Integrated mechanical computing for autonomous soft machines

Material characterization

Dogbone-shaped samples were 3D printed with the soft material that comprises bistable elements (Elastic 50 A, Formlabs), and uniaxial tensile tests were then performed according to ASTM D412-C specifications to measure the linear elastic and viscoelastic property, yield stress, and long-term toughness. We used a universal testing machine (Model 5942, Instron®, Illinois Tool Works Inc.) coupled with a non-contacting video extensometer (Model AVE2, Instron®, Illinois Tool Works Inc.) to perform uniaxial tensile tests at a crosshead loading rate of 25 mm min^-1 while using a 50 N load cell to measure the force. The Elastic 50 A material was found to be linearly elastic with an elastic modulus of 2.04 MPa and a Poisson’s ratio of 0.37 (Supplementary Fig. 2b). A set of step strain experiments with two different strain levels (20%, 40%) was conducted to obtain the viscoelastic behavior of the Elastic 50 A material by applying a fast, instantaneous displacement (4 mm at 120 mm min^-1) for 300 s (Supplementary Fig. 2c, d). We then fitted the Maxwell form of the standard linear solid (SLS) model of linear viscoelasticity to our experimental data; where the instantaneous force (F(t)) can be represented as F(t) = F₀ + F₁e^−t/τ, where F₀ and F₁ are the time independent and time dependent components of the total force, respectively, and τ is the time constant (Supplementary Fig. 2e). The results suggest that the stress relaxation effect was almost negligible after 3 min, and that our SLS model had a high correlation with the experimental data (R = 0.981) and showed a very low F₁/F₀ ratio (= 0.08), indicating that the material demonstrates primarily pure elastic behavior. Apart from the viscoelastic behavior, we also noticed that the yield stress of the Elastic 50 A material is ~1 MPa (at the yield strain of 50–60%; see Supplementary Fig. 2b) and the maximum local stress applied to the beam in its preset buckling state is ~0.7 MPa (Supplementary Fig. 2g), which implies the maximum local strain of ~40%. Based on this investigation, we additionally performed a step strain experiment of 40% strain until the sample was fully damaged (Supplementary Fig. 2f). The result shows that the developed soft bistable beam could sustain its preset buckling state up to ~75 min (~4500 s).

Dynamics and design optimization of a soft bistable unit element

Our design strategy for unit bistable elements aims to achieve the compact and functional assembly of bistable beams with different bistable profiles within a unit lattice, which possibly forms into various mechanical computing components. Previous studies using tilted beams^10,11,12 or 2D triangular tessellation¹³ achieved bistable properties in a relatively simple manner, but the attainable range of their input/output potential energies has been limited solely to the geometries of the unit building blocks. We addressed this challenge by introducing an additional independent parameter, the lattice constant (a) (Supplementary Figs. 2j and 3). According to this approach, even a bistable beam of the same length (L) can have various potential energy profiles depending on the value of a (<L). Interestingly, when being placed into the lattice frame, soft bistable beams with predefined biased beam angles (θ) were buckled into almost the same geometry regardless of θ, thereby enabling easy and symmetric design of mechanical computing units composed of several bistable elements (Supplementary Figs. 2j and 3). For the dynamics and general input/output characteristics, classical beam theories suggest approximate solutions for important parameters in a typical bistable force-displacement profile as shown in Supplementary Fig. 2a, which read^41,42

where c₁ to c₅ are non-dimensional constants, E is the Young’s modulus, I is the moment of inertia, and h is the apex height of the buckled geometry. Importantly, here, h is given as a function of a and L, such that h(a, L). Introducing the detailed expression of I, we obtain:

where b is the depth, and t is the thickness of the beam. The approximation of neglecting high modes of buckling simplifies the force-displacement profile into the shape as shown in Supplementary Fig. 2a, which provides the estimation of potential energies:

Based on this theoretical background, we carried out a parametric study for optimizing the design of a soft bistable unit element. A wide set of numerical studies was performed to provide a possible design space for the parameters, t/L and θ (Supplementary Fig. 2g–i). To reduce the design space, we determined the value of a equal to 9 mm in order to secure a sufficient space in the unit lattice, and the depth (length in the z direction) of all systems was fixed with 4 mm for ease of manufacturing and placement. Aside from the feasible design space obtained from the numerical studies, practical limitations also occurred due to our manufacturing technology. Given the predefined geometric conditions for a and depth, for example, the beam length L needed to be designed within the range between 10 and 13 mm, otherwise it was unable to place the beam into the lattice frame with given a = 9 mm (in case L < 10 mm), or the beam underwent severe deflection together with excessive compressive stresses when initially buckled (in case L > 13 mm). In a similar context, the thickness t should be designed thinner than 0.7 mm. Moreover, the initial beam angle (θ) over 35° induced undesired instabilities in the state of preset buckling. Given all these considerations, therefore, the experimentally feasible range out of the proposed design space was reduced to 10 mm ≤ L ≤ 13 mm, 0.5 mm≤t ≤ 0.7 mm, and θ ≤ 30°, which, in turn, leads to the design space with t/(L/2) = 0.09–0.12 and θ ≤ 30° (Supplementary Fig. 2h, i). Lastly, we set the optimal parameter set of L = 12 mm and t = 0.6 mm to possess a sufficient amount of potential energies in buckled states for the design of mechanical transmission lines and computing units and while avoiding excessive compressive stresses. Supplementary Fig. 2k, l shows representative bistable profiles of beams (L = 12 mm, t = 0.6 mm) with θ equal to 0°, 15°, and 25°. We found that autonomous propagation starting from the preset buckled state was inhibited in the beam with θ smaller than 15°. Finally, the bistable beam with a parameter set of (a, L, t, θ) = (9 mm, 12 mm, 0.6 mm, 25°) was set as an optimized design, which delivers a net transmissive energy ΔU (= U_out ‒ U_in) of ~0.15 mJ.

Characteristic energy analysis for computational propagation

The input (that is, energy barrier, \({U}_{b}^{T}\)) and output (\({U}_{{out}}^{T}\)) characteristics of a mechanical transmission line were analyzed in terms of potential energies applied to each unit element, a soft bistable beam with the given optimized geometry of (a, L, t, θ) = (9 mm, 12 mm, 0.6 mm, 25°). Specifically, considering that the transmission line of finite length consisting of N ( = N₁ + N₂»N_soliton) bistable elements is a series of the transmission line 1 (N₁) and line 2 (N₂) as shown in Fig. 2b, the input energy transferred to the first element of the transmission line 2 (that is, the N₁+1^th element highlighted in green bold in Fig. 2b), \({U}_{{in}}^{{N}_{1}+1}\), is given by the output energy of the transmission line 1, \({U}_{{out}}^{T({N}_{1})}\). Similarly, the output load of the N₁^th beam, \({U}_{l}^{{N}_{1}}\), is determined by the effective energy barrier (or accumulated impedance) caused by the existence of the transmission line 2, \({U}_{b}^{T({N}_{2})}\). An additional point to consider here is that in homogeneous transmission with N»N_soliton, the input/output characteristics of any arbitrary element must be the same. It indicates that the characteristic potential energies of the N₁+1^th beam should be equal to those of N₁^th (highlighted in black, Fig. 2b), N₁-1^th (highlighted in blue, Fig. 2b), and eventually one down to the critical N_c^th beam, which delineates:

Combining the second terms in each parenthesis, we can reduce the equations as (Fig. 2b):

Since N_c can be approximated by N_soliton based on the hypothesis, the input/output characteristic potential energies of the transmission line is finally obtained as (\({U}_{b}^{T}\), \({U}_{{out}}^{T}\)) ≈ (\({U}_{b}^{T({N}_{{soliton}})}\), \({U}_{{out}}^{T({N}_{{soliton}})}\)). Given this expression, a sufficient condition for autonomous propagation within the homogeneous transmission line is obtained as \({U}_{{out}}^{T({N}_{{soliton}})}\) > \({U}_{b}^{T({N}_{{soliton}})}\).

The above investigation can be equally extended to the case with the inclusion of a computing unit, mechanologic X (Fig. 2f, g). Any mechanical computing system driven by transition wave-based computational propagations can be readily divided into two parts: (i) the transmission line 1 (N₁) before the computing unit (mechanologic X), and (ii) a series of the logic X and the transmission line 2 (N₂) connected after it. As explored in the transmission line case, the input energy transferred to the logic X, \({U}_{{in}}^{X}\), is given by the output energy of the transmission line 1, \({U}_{{out}}^{T({N}_{1})}\), such that \({U}_{{in}}^{X}\) = \({U}_{{out}}^{T\left({N}_{1}\right)}\) ≈ \({U}_{{out}}^{T\left({N}_{{soliton}}\right)}\). Similarly, the output load of the N₁^th beam, which is shared by the transmission line 1 and the logic X, is determined by the effective energy barrier induced by the logic X and the transmission line 2, \({U}_{b}^{X+T({N}_{2})}\) ≈ \({U}_{b}^{X+T\left({N}_{{soliton}}\right)}\). Consequently, a sufficient condition for autonomous computational propagation via the mechanologic X is obtained as \({U}_{{out}}^{T\left({N}_{{soliton}}\right)}\) > \({U}_{b}^{X+T({N}_{{soliton}})}\). Note that this condition is valid under our assumption that the mechanologic has the characteristic dimension (lattice constant, a) shorter than the soliton width; otherwise, dynamics at the output interface of the logic should be additionally considered (Supplementary Fig. 9b, d).

Estimation of N
_soliton

The N_soliton values in our approach were estimated in two different ways: (i) First, we experimentally measured or numerically obtained the time-evolving normalized displacement data of all bistable elements involved in computational propagation, and then counted the number of moving beams at any arbitrary time frame. The state of moving beams was defined as when the normalized displacement had a value between 0.05 and 0.95. The subdivided steps involve counting the number of intersections between the displacement profile of each beam and arbitrary vertical lines on the time axis (Supplementary Fig. 8j). The obtained N_soliton data were plotted along with the normalized time (Supplementary Fig. 8k), which produced the average N_soliton of each computing event with calculable error bounds represented by SD (Supplementary Fig. 8l). (ii) Second, we obtained N_soliton by experimentally measuring the input/output characteristic potential energies of the transmission line with N varying from 0 to 9 (Fig. 2c, d and Supplementary Fig. 5). As well described in the main text, we found that, with the increase in N, the force-displacement and energy-displacement curves asymptotically approach the saturated profiles, supporting the existence of N_soliton (Supplementary Fig. 5). These two methods yielded exactly the same result showing the N_soliton value of ~2.5.

Design of mechanologics

All mechanologics take the form of bistable square structures, with designs fit compactly into a unit lattice (a = 9 mm) (Fig. 3a,b and Supplementary Fig. 10). We designed their computing functions using at most four bistable building blocks, which were arranged perpendicular to each other via linear elastic linkages. NOT gates, which need to have two different phases (ϕ=0, π) depending on the phase of the input signals, consist of four identical bistable elements with a parameter set of (L, t) = (12 mm, 0.6 mm) together with θ = 15° (for ϕ=0) or 25° (for ϕ=π). Note that, within the same unit lattice, we could easily tune the potential energy profiles developed in each of the four bistable elements by changing L, t and θ (Fig. 3a). Given the practical design space (Supplementary Fig. 2k, l), the values we proposed, especially for θ, indicate the minimum values allowing for autonomous computing, that is, computational propagation without additional environmental or human intervention. All the bistable elements in NOT gates are connected with their neighboring components via a linear elastic linkage. The geometry of the linear linkage is divided into two parts having different width, 0.4 and 1 mm, respectively. The difference in bending stiffness between these two parts, which are ~15.6 times [=(1/0.4)³], makes the linkages behave like a rigid four-bar linkage so as to enable the inversion of transmitted mechanical bits, while the thinner parts permit the rotation in the proximity of bistable beams (Supplementary Movie 3). AND gates basically take a similar form as in NOT gates. Differences include the presence of two input terminals, which are arranged mirror-symmetrically, and additional geometric engineering of bistable beams to carry different potential energies for the input and output terminals (the output one should have a larger energy barrier; see Supplementary Fig. 10). To achieve this, we designed the beam angles of the input (θ₁) and output (θ₂) terminals to have the values of 20° and 5°, respectively, which delivered successful computational propagation (Supplementary Movie 4). Distinct from other mechanologic designs, OR gates consist of three bistable elements connected via a softer and longer elastic linkage. This open-ended design is helpful to decouple the actuations of two independent input beams. The geometry of the spring-like linkage can be replaced by alternative structures, only if its total length matches the distance between the input and output terminals so that the reversion of the output terminal does not affect the non-triggered input. The beam angle (θ) was designed to have the value of 30° for successful computation (Supplementary Movie 5).

Fabrication of integrated mechanical computing systems

All bistable beams, transmission lines, and logic units were fabricated monolithically by 3D printing using an elastic material (Elastic 50 A) in a stereolithography printer (Form 3B, Formlabs). All 3D printed parts were cleaned by ultrasonicating (USC 1200 T, VWR) in isopropyl alcohol for 30 min, and subsequently post-processed in an ultraviolet (UV) oven (Form Cure, Formlabs) for 90 min at 60 °C. The lattice frames with the lattice constant (a) equal to 9 mm were 3D printed using a rigid material (Clear, Formlabs) and spray painted black to provide contrast for better visualization (Supplementary Fig. 3a).

Mechanical characterization

The uniaxial tensile testing setup same as that for the material characterization above was used to obtain the force-displacement characteristics of the unit bistable beams, transmission lines, and mechanologics (Fig. 3a, right, and Supplementary Fig. 5a). All the prepared samples were placed in the printed lattice frame, which was affixed to the bottom clip. Specifically, two holes spaced 2.5 mm apart were patterned onto every sample and anchored by two separate cylindrical pins (diameter = 0.4 mm) to allow only linear displacement by preventing undesired rotation of the load point. A displacement-controlled loading was applied by specifying a crosshead loading rate of 4 mm min^-1 and while using a 5 N load cell to measure the force. For most measurements, at least three different samples were prepared for each design and measured at least three times for each sample. All plots represent the average values with error bounds indicating standard deviation (SD).

Experimental analysis of computational propagation of mechanical signals

Mechanical computing systems with bistable elements being marked with color dots were placed in the lattice frame. The behavior of computational propagation was recorded and then analyzed by a free video analysis and modeling tool (Tracker 6). The obtained displacement data were analyzed in their normalized forms which were defined by the net displacement of each bistable unit in the propagation direction divided by the average maximum displacement (~4.6 mm) of all the reverted beams.

Finite element analysis

Finite element (FE) analysis was performed to model the dynamic behavior of individual bistable elements as well as the propagation of transition waves along a lattice of multiple bistable elements connected by elastic links. All FE simulations were performed with a commercial finite element solver (Abaqus/Standard, version 2020, Simulia, Dassault Systèmes). Plane strain conditions were assumed to increase computational efficiency and elements of type CPE8 were used to construct the mesh. All simulations involved three steps: (i) An initial static step to compress the beam so that it matches the lattice length. (ii) A following static step with force-controlled loading to invert (preset) the beam to its second stable state. (iii) A final dynamic implicit step with displacement-controlled loading to revert the beam to its first stable state while recording the reaction force. Using a dynamic analysis enabled us to accurately capture the instabilities and snap-through behavior of the bistable beams. A displacement-controlled loading was necessary to obtain the complete equilibrium path in the force-displacement plane. For the final dynamic step, the two ends of the beams were fixed with an encastre boundary condition and the center of the beam displaced vertically. The energy-displacement data were obtained by integrating the recorded force-displacement data. We ensured quasi-static conditions in the simulations by introducing a small damping factor and by monitoring the kinetic energy of the simulations. Transition wave simulations required the same steps, with all the beams being subjected to steps (i) and (ii), and only the first (end) beam being provided a prescribed displacement in the final step. All other beams are allowed to deform freely as the transition wave propagates through the system and the displacement of the midpoint of each beam is recorded to monitor the wave propagation.

Fabrication of soft hydrogel actuators

A solution of photoinitiator was prepared by dissolving 150 mg of 2,2-dimethoxy-2-phenylacetophenone (99%, Sigma Aldrich) in 1 mL of dimethyl sulfoxide (Sigma Aldrich). Ethylene glycol dimethacrylate (100 µL) and the photoinitiator solution (150 µL) were added in 2 ml of poly(ethylene glycol) diacrylate (PEGDA, Mn 575, Sigma Aldrich). An 8 µL of the prepared PEGDA solution was dropped onto the surface of 3D printed actuator body made up of the Elastic 50 A material (25 mm × 4 mm) and then covered with a slide glass to adjust its thickness to 25 ± 5 µm. The confined solution was polymerized by irradiation of UV light (Bio-Link 365, wavelength of 365 nm, Vilber Lourmat GmbH) for 30 min.