Tailoring high precision polynomial architected material constitutive responses via inverse design

As a first step to our general design methodology, we examine how to achieve broad classes of nonlinearity in a structure. Qualitatively, a softening nonlinearity can be achieved by transitioning a beam-like structure via large deformation and rotation from an axial-dominated loading regime to a bending dominated loading regime, while a stiffening nonlinearity can be achieved by reversing this transition (see, e.g., Figure 2 A). By tailoring the shape of the structure, we can adjust both the amount of softening and stiffening and where these regimes exist in the deformation range. This can be done in either a single “spring” (in which, herein, we refer to a beam-like structure undergoing large deformation as our spring), or through a combination of multiple springs in series or in parallel forming an effective spring. We note that within the context of linearized stiffnesses, the stiffness of the overall structure will be dominated by the stiffness of the softer spring for springs in series, while the stiffer spring will dominate for springs in parallel. The resultant behavior can thus be tailored to produce the desired type of nonlinear behavior in the desired displacement or force region.

Parameters of importance can be identified by examining the case of a straight or curved beam being pulled through a fixed displacement range as shown in Figure 2 A. A spring of this general design could easily be incorporated into a spring-mass chain, with a mechanical frame present to enforce rigid boundary conditions, as displayed in Figure 1 E. Considering a system of this design, we take the beam highlighted in the dotted red box in Figure 1 E as our design domain. There exist inherent limits in how stiff or soft this structure can be made, and similarly, the magnitude of forces that can be reached in a given displacement range. Consider first the blue dashed upper line in Figure 2 B, which corresponds to the limit $EAL^{-1}$ and describes linear stiffness corresponding to the axial deformation of a long, thin elastic beam, where $E$, $A$, and $L$ are the beam’s elastic modulus, cross sectional area, and length, respectively. This limit can be interpreted to define the approximate maximum force that could be achieved at a given deformation $\delta$, by pulling a beam that extends normally from a fixed boundary, in an axial fashion. Such a regime could be imagined to be achieved when the beam in Figure 2 A is pulled far to the right ($\delta>>L$). This further serves as a small ($\delta<<t$, where $t$ is the beam width) and very large deformation ($\delta>>L$) limit on the structural stiffness, with implications on for the effect of volume ratio in an inverse design method on achievable stiffnesses. An exception to these limits, as concerns local stiffness at intermediate displacements is the case of snap through and snap back (or snap down), where infinitely positive, infinitely negative, and zero local stiffnesses can be achieved [36]. For instance, considering a structure that exhibits a single snap through or snap back feature, we could expect the following. For both a snap-through and snap-back structure, there is guaranteed to be two points of zero stiffness. For a snap-through structure, there is also guaranteed to be a point of finite negative stiffness which approaches a vertical slope in the limit of snap-back behavior. This vertical slope could be considered as an infinitely positive or negative stiffness. For a snap back structure, there is similarly guaranteed to be two points of vertical slope. Similarly to the blue dashed line forming the upper limit, we can define an approximate lower bound, denoted by the green dashed line in Figure 2 A for stiffness in small deformation ($\delta<<t$) regimes, which corresponds to the $EIL^{-3}$ bending stiffness of a long thin beam, with $I$ being the second moment of area. The exception to this bound would be to include contact effects, or to include a pre-load on a structure exhibiting negative stiffness onset. Without considering contact and pre-loads, this has important implications for the dependence of the small deformation stiffness that can be achieved on the structure’s minimum feature size.

Within the aforementioned bounds, we suggest that a wide range of qualitatively different nonlinear force-displacement paths, which can be fit by a polynomial (as a function of displacement/strain), should be achievable through this exploitation of geometry, even those which elsewhere in literature have been achieved through other mechanisms, mechanical or otherwise (e.g., magnetism, etc.). Further, while these qualitative design guidelines are helpful, we look herein to topology optimization techniques to refine the design process for more exact targeted nonlinear behaviors.

A 2D plane stress, geometrically nonlinear (large deformation), displacement controlled, level set optimization was implemented, following the method described in [37, 38]. For simplicity, we restricted this work to a tension regime. A Kirchhoff material model was implemented, allowing for geometric nonlinearity and large strains while maintaining a linear elastic material model. We choose a level set optimization method because of its compatibility with our pursuit of simply manufacturable, binary (material and no material) geometries.

Within the context of the stiffness bounds described as part of Figure 2, we highlight that there exists an inherent coupling between the degree of nonlinearity and the stiffness of a structure (excepting phenomena related to snap through and snap back). We define degree of nonlinearity as the relative magnitude of nonlinear terms in the force-displacement relationship with respect to the linear term at a given displacement. To elucidate the coupling between nonlinearity and stiffness, we consider the transition between axial and bending deformation of a long thin beam. Taking a ratio of those two characteristic stiffnesses (axial by bending) for a beam of rectangular cross section, we obtain a ratio of $(Lt^{-1})^{2}$, demonstrating that thinner structures are capable of greater change in stiffness (i.e., nonlinearity) over a fixed displacement range. However, thinning a structure makes it softer, resulting in a uniformly lower magnitude force displacement curve. As mentioned above, many previous attempts at this problem have employed an objective function which attempts to directly target a force-displacement curve. This severely limits the design space, excluding multiple, and sometimes all, viable solutions. In that case, a better goal, in pursuit of a specific nonlinear behavior, is to search for geometries which will enable the desired strength and shape of nonlinearity, regardless of the stiffness of the structure. We address this by defining the objective function in terms of a polynomial that is fitted to the force-displacement curve, wherein the objective is set to match a desired ratio between polynomial coefficients. This objective is written as

Sensitivities for this objective function can be written as a combination of sensitivities for compliance and sensitivities related to the fitting of the polynomial. We calculate self-adjoint, nonlinear, displacement control compliance sensitivities, $\partial C_{j}/\partial\Omega$, as calculated in [38], where $C$ is the compliance and $\Omega$ is the domain. Calculating our polynomial fit through use of the equation $P=V^{-1}C$, where $V^{-1}$ is the inverse of the Vandermonde matrix, the sensitivity for each polynomial term is written as:

Despite the ability to identify gradients in the design space, we note that this particular inverse design problem is highly non-unique and non-convex. Hence, optimization can get stuck in local minima which may not be sufficiently close to the objective nonlinearity. In order to improve the quality of a solution, we use mechanical intuition to choose an initial condition that is qualitatively “nearby” the desired type of nonlinearity, through the heuristic spring design process described above. This narrows down the number of possible solutions the optimizer can reach, and reduces the chances of becoming trapped in an undesirable local minimum. A brief exploration of the effect initial condition choice can have on the ability to find a satisfactory local optimum is included in the SI. This approach, presented herein, requires some mechanical insight and designer input in choosing an initial condition. Similar approaches have been employed in other works, which have relied upon the input of reduced order training datasets [28] and the use of multiple initial condition guesses [27] to achieve the target constitutive behavior.

To demonstrate the capability of the design method and optimization process, we optimize several spring structures to target specific nonlinear constitutive laws from several classes of nonlinear behavior, i.e. softening, stiffening, softening to stiffening, and stiffening to softening. By optimizing a structure for each class of nonlinearity, we demonstrate the capability of the approach to span the nonlinear design space and target any type of nonlinearity that falls within the range of these classes. These results are depicted in Figure 3, in which the force-displacement equations with the following polynomial coefficients are targeted: Softening, $F=x-40.0x^{3}$; stiffening, $F=x+163.3x^{3}$; softening to stiffening, $F=1.40x-17.332x^{2}+166.875x^{3}$; and stiffening to softening, $F=0.599978x+17.332x^{2}-166.875x^{3}$. The displacement range was chosen arbitrarily, as were the coefficient weights, though in the latter two cases, the target polynomial coefficients were chosen to achieve symmetric stiffness switching springs, as outlined in [35]. Nearly all targeted nonlinearities were matched to within less than 2% of the desired coefficients (with the sole exception being the softening-to-stiffening case, discussed below). In the optimizations shown in Figure 3, the initial conditions were selected to be in the same class of nonlinearity as the desired target (for example, an initial condition which displayed a general stiffening nonlinearity was selected for the stiffening-type nonlinearities), with coarse adjustments being required to tune the initial polynomial coefficient signs. The initial conditions and polynomial fit details for each optimization case can be found in the SI. We suggest that all of the optimized designs can be considered “non-intuitive”, in that they would be difficult or impossible to design by hand to the desired degree of precision using human intuition alone. We draw particular attention to the two stiffness switching optimized designs shown in Figure 3 C and D, as these targeted constitutive laws were more complex in nature, requiring the optimizer to capture a qualitative shift in behavior and to perform optimization for two polynomial ratios simultaneously. The stiffening-to-softening spring design in Figure 3 D, in particular, is of note in that it can be considered to be the combination of two springs (beams) in series (one stiffening spring and one softening spring) in order to achieve the desired behavior. These two springs are highly coupled, with the nonlinear behavior of each affecting not only the overall behavior of the structure, but also the boundary conditions and applied displacement of the other spring, creating a highly complex and nonlinear design scenario. Although the fitness of the match is slightly lower for this case, at an error of 5.242%, this is less surprising considering the greater structural complexity (and presumed design space complexity) seemingly required to achieve this response.

The optimization results demonstrate that the method is capable of designing a variety of nonlinear springs. The high degree of precision in matching a target behavior also makes it possible to target nonlinearities which vary from one another only slightly, sweeping through the design space for the desired nonlinear behavior. To demonstrate this capability, we chose a target which consisted of a linear plus cubic stiffening type nonlinearity ($F=x+Ax^{3}$), and varied the target coefficient weights (A is varied). This type of sweep is fairly straightforward to perform once the first design has been optimized, as the optimum from one design can be used as the initial condition for the next step, in many cases (but not all) eliminating the need to design the initial conditions. The results are shown in Figure 4. We first draw attention to the distinct shift in design between ratios 40 and 80. This transition represents a local minimum which led to a relatively poor optimization of the ratio 40 target, and prevented an acceptable optimization of the ratio 20 design from an initial condition drawn from this poor ratio 40 optimization. Intuitively, this makes sense, as the hinge-like structure present in the more nonlinear designs is qualitatively different from the type of structure which might be expected to exhibit more linear behavior. Selecting a new starting design, which in this case consisted of straight edges and a thicker beam to promote a weaker nonlinearity, provided a successful workaround which avoided the local minimum and found a satisfactory design for the ratio 20 target, as well as an improved optimal design for the ratio 40 target. This difficulty highlights the challenges of the highly nonlinear design space, and the ease with which such local minima can be circumvented with educated initial design. We also note that the coupling of stiffness and nonlinearity is clearly visible in the optimized curves of the sweep designs (considering ratios 40-163 to be grouped together, and ratios 10-20 to be grouped together) in Figure 4B, in which the more linear designs display stiffer overall behavior, while the more nonlinear designs are softer. This change in nonlinearity (and stiffness) is achieved by increasing or decreasing the thickness. We highlight this in the zoomed-in images of the ratio 40-163 optimized structures in Figure 4, noting that the ratio 163 design is the optimized stiffening design from Figure 3, and that each subsequent ratio uses the optimum from the previous ratio as the initial condition.

With respect to the coupling between stiffness and nonlinearity, while, as in Figure 1, it is simple to normalize or scale the linear stiffness (or resulting force magnitude) computationally, we also note that there are several mechanisms to achieve this in practice. Perhaps the most straightforward manner to scale the resulting stiffness or effective modulus would be to adjust the elastic modulus of the constituent material forming the design. However, in practice, materials occupy discrete moduli ranges and come with other practical limitations such as small elastic ranges. As an alternative approach, consider first the structural stiffness of the design. To keep the same degree of nonlinearity, while increasing the magnitude of the force (thus compensating for the normalized stiffness), one could change the depth of the 2D geometry, as shown by the varied colors in Figure 1C and the resulting force-displacement curves in Figure 1D. One could also scale the area of the 2D geometry, thus increasing the force, however, the applied displacement must also be proportionally scaled to achieve the same degree of nonlinearity. We next consider scaling the response from the perspective of an effective stress-strain response. Consider the 2D design with a nonlinear stiffness $k(\delta)$ is assigned some depth $d$, length $l$, and width $w$, as shown in Figure 1C. This structure is then repeated within a larger unit-cell domain of depth $D$, length $L$, and width $W$, wherein the designed region is attached to a rigid frame as shown in Figure 1E. The resulting stress strain relationship would thus be $\sigma=\epsilon kL(DW)^{-1}$, such that the size of the unit cell becomes essentially a single scaling factor equivalent to scaling the modulus of the constituent material. There is a key difference from scaling the constituent modulus, however, where the lengths defining the design domain cannot exceed that of the larger, unit cell domain.

Furthermore, we include several additional examples of optimizations in order to demonstrate the range and versatility of the method. While the examples above highlight distinct nonlinear behaviors, several of the targeted curves are weakly nonlinear. To demonstrate that the method is not constrained to weak nonlinearities, Figure 5A highlights a targeted stiffening curve displaying a greater degree of nonlinearity (and a purely quadratic constitutive law). The target in this case is $F=x+200.0x^{2}+0x^{3}$. We note that, when the maximum displacement of 0.07 is substituted into $x$ in this target, the ratio of the quadratic term to the linear term becomes 14, demonstrating an order of magnitude difference between linear and nonlinear terms. In Figure 5C, a constitutive law containing a fourth order term is targeted in order to show that the above cubic examples are an arbitrary choice of nonlinearity, and there is nothing constraining the method to third order terms. In addition to targeting the distinct nonlinear behaviors outlined above, we desire the capability of targeting nonlinearities which lie in the negative stiffness, as well as fully bistable regime (a capability which is highly useful in pursuing superior impact-resistant structures.) [2] One benefit of our displacement-controlled approach is the capability to solve problems which exhibit negative stiffness and snap-through type behavior with a Newton-Raphson solver. In Figure 5B below we show optimizations which target a fully bistable forced-displacement curve. Results of an optimized negative stiffness (but not fully bistable) structure are shown in Figure 5D, for a targeted nonlinear law of $F=x-10.0x^{2}+30.004x^{3}$. For both of these cases, we optimize over a greater displacement range in order to more easily realize the more nonlinear negative stiffness behavior without resorting to a significantly increased mesh density (as would be required for a thinner structure). As above, the displacement range and targeted nonlinearity are otherwise arbitrary. This optimization requires the simultaneous optimization of two ratios. For the initial optimization run, both ratios were matched to within 1% (0.638% for the cubic and 0.636% for the quadratic). However, in this case, the negative stiffness portion of the curve is heavily dependent upon only one of these two ratios (the quadratic, which is responsible for the transition to negative stiffness due to the negative sign associated with the term), and the behavior in the negative stiffness region is accordingly much more sensitive to this ratio. We therefore add a weight (a 10x multiplier) to the quadratic ratio for this optimization, prioritizing this term over the cubic term. This result (shown in Figure 5D) achieved a lower match to the cubic term (1.712%), but a much closer match to the prioritized quadratic term (0.180%), resulting in a closer normalized match of the negative stiffness portion of the curve.

In order to confirm the accuracy of the method and the validity of the optimized springs, several of the springs shown herein were chosen for manufacturing and experimental verification of the nonlinear force-displacement curves. Polycarbonate was initially chosen to fabricate the springs from, due to its relatively linear behavior at lower strains and high strain-to-break value. However, preliminary tests indicated that for thicker structures, plasticity onset (which is not accounted for in the optimization’s Kirchhoff material model) would introduce unaccounted nonlinear behaviors. The thinnest structure, the bistable spring from Figure 5B, was therefore manufactured and tested in polycarbonate, with the results shown in Figure 6A indicating close agreement between optimization outputs and experiments. It is of note that such experimental validation proved to be incredibly sensitive to misalignment issues (such as any out of plane offset or the direction of displacement application not being precisely horizontal) and quality of the boundary condition approximation. Regarding the latter, the imposition of the y-direction being fixed on the top and bottom edges via the frame containing the optimized structure, in particular, was found to be very important, due to a great deal of the nonlinearity stemming from the structure pressing against these surfaces during deformation. Indeed, insufficient stiffness of the frame, or lateral forces applied to it, were observed to modify the resulting nonlinear force-displacement response of the designed structure. Careful construction of the system was thus required in order to achieve good experimental match, which may require additional attention when considering real-world applications. We also note that, due to the FE simulations being run on a 1x1 m domain, often with a modulus of E = 1 Pa, simulation results required scaling to match up to experimental results. This is to be expected, and highlights the ability to manipulate the scale of the structure in order to achieve desired force ranges. Simulation results were scaled along the x-axis by the size of the unit cell (i.e., experimental unit cell length/simulation unit cell length), and along the y-axis by the modulus of the constituent material used in the experiment (found via scaling max and min values of nonlinear curves, as described in the SI), the edge length of the unit cell, the depth, and a factor of two to account for two springs being present in the experiment(i.e. experimental unit cell length/simulation unit cell length), and along the y-axis by the modulus of the constituent material used in the experiment (found via scaling max and min values of nonlinear curves, as described in the SI), the edge length of the unit cell, the depth, and a factor of two to account for two springs being present in the experiment (i.e. ($E_{Exp}/E_{Sim})*(L_{Exp}/L_{Sim})*(D_{Exp}/D_{Sim})*2$, in which E refers to Young’s Modulus, L refers to the edge length of the unit cell, D refers to the depth of the unit cell, subscript Exp refers to the value in the experiment, and subscript Sim refers to the value used in the optimization).

As the ability to test any optimized structure was desired, an updated material model was introduced to the modeling for optimization, which would allow for the use of fabrication materials that exhibit a larger elastic range. A Neo-Hookean (hyperelastic) material model was introduced, with a plane strain formulation. Several of the above optimization cases were re-run with the updated material model, beginning from the optimal configuration identified with the Kirchhoff material model. Small adjustments were made to these designs by the optimizer to account for the new material model. These updated structures were then fabricated out of silicone (see SI for material and fabrication details) and tested. The results, shown in Figure 6B,C, demonstrate that the updated material model has similar levels of agreement to the original Kirchhoff material model, and allows for the fabrication and testing of thicker structures while avoiding the onset of plasticity. We note that analogous scaling was performed to map between simulation and experimental results, as described above in the elastic case. Further information on experimental validation is included in the SI.