Compressive properties of silicone Bouligand structures

This article presents an examination of silicone, Bouligand lattices in compression. Appearing frequently in biological organisms and manufacturing design, Bouligand structures comprise layers of parallel strands or fibers oriented in a helical fashion. They can exhibit exceptional fracture resistance when composed of rigid or composite materials. The behavior of elastomeric Bouligand structures, however, is less well understood. Additively manufactured (AM) elastomeric lattices have applications in stress mitigation, medical devices, and soft robotics. This article demonstrates that Bouligand structures are a useful addition to the design space of AM elastomers. By adjusting the layer-rotation parameters, lattice stress can increase by more than 300% without altering the porosity. Additionally, we introduce path length metrics that help explain the observed relationship between layer rotation and compression response. Additive manufacturing (AM) continues to push the boundary of manufacturable structures and enhance the ability to robustly design for specific properties and behaviors. The more we understand the design space of a novel AM microstructure, the greater its application range. In this article, we describe the mechanical behavior of helicoidal, elastomeric lattices and introduce path length metrics to help explain their stress response. We show that these structures can exhibit a large range of mechanical behaviors in compression, making them well suited for applications such as stress mitigation and impact absorption. Additionally, the path length metrics could become useful design tools and may be applicable to a larger set of cellular structures. These findings expand our ability to rapidly design materials with highly specific and customizable properties to meet the needs of modern engineering challenges.


Introduction
Direct-ink-writing is an additive manufacturing (AM) technique that enables the fabrication of complex three-dimensional shapes that traditional manufacturing methods cannot achieve. 1 Researchers have direct written many materials, including elastomers, 2,3 metals, 4 ceramics, 5 and their composites. 6 The resultant novel structures have advanced fields such as soft robotics, 7 stretchable electronics, 8 biomedical engineering, 9 optics, 10 and stress mitigation, where precise deposition of elastomeric, cellular microstructures enables tunable mechanical properties and superior aging performance than stochastic foams. 11 Two common direct-ink-write (DIW) structures, "simple cubic," and "facecentered tetragonal," exhibit a range of mechanical behaviors with changes in porosity; 2,12,13 however, alternative microstructures exist that can further expand the application space of DIW. [14][15][16] In this work, we add to this expansion by modeling and characterizing elastomeric, Bouligand structures, which comprise layers of straight, evenly spaced filaments that are rotated by a fixed angle relative to each other. Many organisms employ these architectures-typically with rigid materials-to achieve superior fracture toughness. [17][18][19][20][21] Inspired by these organisms' capabilities, engineering researchers have aimed to reproduce their structures and properties. [22][23][24][25][26] The mechanical properties of elastomeric, Bouligand structures, however, are not well researched. We aim to fill this gap in knowledge.
Cushioning and stress mitigation applications often employ elastomeric cellular materials for their ability to endure compressive loads; therefore, we aimed to understand how elastomeric Bouligand Impact Letter Impact statement Additive manufacturing (AM) continues to push the boundary of manufacturable structures and enhance the ability to robustly design for specific properties and behaviors. The more we understand the design space of a novel AM microstructure, the greater its application range. In this article, we describe the mechanical behavior of helicoidal, elastomeric lattices and introduce path length metrics to help explain their stress response. We show that these structures can exhibit a large range of mechanical behaviors in compression, making them well suited for applications such as stress mitigation and impact absorption. Additionally, the path length metrics could become useful design tools and may be applicable to a larger set of cellular structures. These findings expand our ability to rapidly design materials with highly specific and customizable properties to meet the needs of modern engineering challenges.
structures respond in compression. We performed parametric analyses in simulation, focusing on the parameters that define rotation between layers: the angle of rotation and the position of the point about which all filaments are rotated. To compare the nonlinear stress-strain behavior of the different parameters, we report the lattice stress at multiple strains. (We define lattice stress to be the reaction force divided by the crosssectional area.) We validated the simulation results by printing and mechanically testing a subset of the simulated structures.
To help explain the observed relationship between compression response and the rotation parameters, we defined path length metrics representative of the structure and found a strong correlation with lattice stress. Future work will involve performing dimensional analysis to seek an analytical or reduced order model relating lattice parameters to stress-strain response, which could enable more rapid design of DIW microstructures.

Simulation parameter sweeps
A fully defined Bouligand lattice requires six parameters: strand diameter ( d ), horizontal strand spacing ( s ), vertical offset between layers ( vo ), rotation angle between layers (pitch angle, θ), rotation center position ( rc ), and number of layers ( n ) (Figure 1). Using theoretical porosity calculations based on simulated geometry data, we found that porosity increased significantly with increases in vo/d and s/d , while changes in rc and θ affected the porosity by less than 4% (Supplementary Figure 1). Because the effect of structure porosity on lattice stress is well understood, 12 we focused our parametric analyses on the porosity-independent parameters: rotation center position and pitch angle.
Note: To calculate the theoretical porosity ( ), we used the following equation: where V is the volume of the lattice geometry, w equals the geometry width, and th equals the geometry thickness: To understand how pitch angle affects lattice stress, we simulated compression processes for structures with rc = 1 and pitch angles θ {20°, 25°, 30°, …, 85°, 90°} and calculated lattice stress, F/w 2 , where F equals the reaction force output. In general, as the pitch angle approached 90°, the lattice stress increased; however, there were exceptions, most notably at 60° and 45°. We also analyzed the relationship between elastic modulus (for strains 0.07 to 0.12) and pitch angle and found that-for 15% strain-it was similar to the lattice stress versus pitch angle relationship (Figure 2a-b). For larger strains, however, the modulus versus angle relationship deviated from the stress versus angle relationship due to the nonlinearity of the stress-strain curve Figure 1. (a) Lattice parameters and renderings of simple cubic, face-centered tetragonal, and helicoidal lattices, (b) adjacent-layer geometry for 30° and 90° rotation angles, (c) bird's-eye view of printed lattices. Figure 2). To fully capture this nonlinear behavior, we focused our analyses on stresses at multiple strains.

(Supplementary
To analyze the effect of rotation center position on lattice stress, we simulated all combinations of rc ∈ {0.0, 0.25, 0.5, 0.75, 1.0} and θ {20°, 25°, 30°, …, 90°} and found that for some pitch angles, the lattice stress varied by up to 220%, whereas for others, the stress changed by less than 8% (Figure 2c-d). This result reveals the importance of precisely defining the rotation center for certain pitch angles, which becomes especially relevant when printing lattices on nonplanar surfaces, where the local rc values can deviate from the intended.
The increasing trends in lattice stress versus pitch angle may be explained by inspecting the geometry of two adjacent layers ( Figure 1b). As the pitch angle approaches 90°, intersection points are pushed closer together, increasing their frequency. Therefore, lattices with pitch angles that are closer to 90° have more intersection points in a given volume. Strands at intersection points experience more stretching than unsupported strand segments (which can bend more), and research has shown that stretch-dominated lattices exhibit higher lattice stresses. 27, 28 Although the frequency of intersections helps us understand why lattice stress increases as pitch angle increases, it does not explain the spikes in stress at 45°, 60°, and 90°. To elucidate the origin of these spikes, we introduce new lattice metrics ( L p and C p ) in the "Lattice vertical path length" section.

Lattice vertical path length
We discussed in the "Simulation parameter sweeps" section that intersection point frequency plays an important role in lattice stress; however, relative intersection point position also has a significant impact. Consider structures with a pitch angle of 90°: when rc = 0 , all intersection points align vertically, vertical stress columns form in compression, and the structure has a high stress response (Figure 3a). When rc = 0.5 , however, the intersection points are maximally misaligned, and the structure is among the softest of the structures we simulated.
To quantify the relative position of intersection points in a structure, we defined a variable called "path length," L p , to be the length of the minimum path needed to traverse along filaments from an intersection point in the bottom layer, to the top. To understand this metric, imagine an ant on the lattice. To take the shortest path from bottom to top, it would traverse a filament in the first layer until it reached an intersecting filament in the second layer. It would climb up to that layer and repeat this process until it reached the top. Movie 1 in the supplements is a visual representation of an L p calculation.
We calculated L p for all intersection points between the bottom two layers in a six-strand-wide, square boundary and defined a variable called "path length coefficient," C p , which equals the mean of all L p values divided by the geometry thickness (Supplementary Figure 3). The equation for C p is as follows:

COMprEssiVE prOpErtiEs Of siLiCOnE BOULigAnd strUCtUrEs
Here, n p equals the total number of paths. As C p increased, lattice stress decreased (Figure 3b-c, Supplementary Figure 4). This correlation can be understood by considering that C p is a surrogate for the extent to which compressive stresses are distributed within the structure. For example, at a 90° pitch angle with rc ∈ {0,1}-where intersection points are laterally aligned and C p is low-stresses are concentrated to the columns of aligned intersection points, resulting in a stiffer response. Conversely in structures where intersection point alignment is poor and C p is high, the stress is more distributed, and the overall lattice stress is low.
The correlation between C p and lattice stress helps to explain why some structures with pitch angles of 45°, 60°, or 90° exhibited dramatically different stress responses than others. These pitch angles are large factors of 180°, and therefore their structures form short helicoidal periods (four, three, and two layers, respectively). When helicoidal periods are short, the lateral positions of intersection points repeat often, increasing the frequency of their alignment or misalignment, which can minimize or maximize C p , pushing the mechanical response to the extremes of the design space. This observation is consistent with research that has shown that structures with short helicoidal periods can exhibit unique responses when compared to structures with longer helicoidal periods. 18 Having observed that lattice stress decreases with an increase in both C p and intersection point spacing (~1sinθ, Figure 1b), we combined the variables into one expression, sinθ/C p , and plotted those values as a function of pitch angle (Figure 3c, dashed line, Supplementary Figure 5). This curve, plotted on a separate y-axis, overlaps well with the simulation stress curve. Future work will seek an analytical expression relating structure parameters to lattice stress by deriving a scaling factor for this expression.

Experimental validation of models
We printed and mechanically tested a subset of the simulated architectures (Supplementary Table I) and found that, in general, the experimental stresses were higher (Figure 4a). To understand why the experimental stresses exceeded the simulations, we observed that the printed lattices experienced sagging, and therefore, their measured thicknesses were below the commanded value (Supplementary Figure 6). As a result, the printed samples had lower porosities than their simulated counterparts, and we know that lower porosity structures have higher stresses in compression. 12 To determine if the disagreement could be fully explained by the thickness mismatch, we ran updated simulations using new vertical offset values, which we calculated by applying the measured thickness values to Equation 2. (For instance, to compute the new vo for the lattice with rc = 0 , θ = 50°, and a printed thickness of 1.5 mm, we used the expression 1.5 mm = 0.2 mm + vo· 9 and found that vo = 0.14 mm.) The updated simulation results matched the experiments much more closely ( Figure 4b); however, agreement remained less good at higher strains. Lattices with pitch angles 45°, 60°, and 90° had the largest error, where the simulation stresses exceeded the experiments for strains above ~20 percent. These are structures with short helicoidal periods, and we hypothesize that frequent, repeated layers might experience more buckling, which-being a high-acceleration, short time-scale process-cannot be captured by the implicit computational approaches we used. Comparing stress-strain curves supports this hypothesis (Supplementary Figure 7). For the 60° and 90° lattices, deviation between simulation and experiment increased dramatically at the strains where the experimental curves began to plateau, and plateaus of this kind often indicate buckling. 28,29 To improve the simulation accuracy, we could employ a computational model that captures buckling, such as explicit modeling. Another possible source of error is that the simulations were modeled as confined compression, whereas the experiments were unconfined. We assumed that larger structures (such as the ones we mechanically tested) would experience negligible lateral movement in a six-strand-wide region; however, this assumption might not hold for larger strains. Additionally, some aspects of the simulation geometries were idealized: the strands were modeled as perfect cylinders, while the printed parts sagged, the bottom layers' strands flattened against the print surface, and the models omitted the excess material that exists at the intersections of the printed parts (Supplementary Figure 8). Each of these features could be included in future models to reduce error.

Discussion and conclusion
In this work, we analyzed two porosity-independent variables-pitch angle and rotation center position-and found that they affect the lattice compressive stress by up to 326 percent. We also introduced path length metrics, L p and C p , that help to explain these results. The correlation between the path length coefficient, C p , and lattice stress suggests that it is possible to derive an analytical or reduced order model that predicts lattice stress at a given strain from structure parameters and material properties. Because Bouligand structures are parametrically defined, it may be possible to derive an analytical form for L p that could serve as a basis for the C p expression. Such a model would greatly improve the ability to rapidly design structures for a given application, because it would reduce or eliminate the need to run FEA simulations.

COMprEssiVE prOpErtiEs Of siLiCOnE BOULigAnd strUCtUrEs
To the best of our knowledge, this is the first examination of how rotational structure parameters affect the compressive properties of elastomeric Bouligand structures. These results facilitate the design of additively manufactured, elastomeric microarchitectures, and the methods could be transferable to rigid materials. By elucidating the compressive responses of these structures, this work enables additive manufacturing design to access more mechanical behaviors, enabling a larger application space.

Compression simulations
The main goal of this work was to expand our knowledge of the Bouligand structure design space for elastomers. Average structure simulation time was 3.2 h, enabling us to gather data for many different structures relatively quickly. By contrast, mechanically testing one sample required 1 h to print, 17 h to cure, and 10 min to test.

Simulation software
We chose Ansys as our finite element analysis (FEA) software because it allowed us to automate the entire simulation process, from generating the geometries, to gathering the desired output variables. Moreover, the static structural algorithm achieved higher compressive strains for DIW architectures than other commercial or custom codes.

Geometry definition
We defined the outer boundary of the structure as a square with width equal to 6 s. We chose this value by drawing from research that shows that to minimize edge effects when modeling lattice structures in compression, the width of the boundaries must be at least as wide as six strands/cell units. 30 We meshed the lattice geometries using linear-element tetrahedrons with four degrees of freedom at each node (displacement and hydrostatic pressure). The edge length could be as large as 0.25 × d . To simulate the compression, we defined two square plates directly above and below the geometry and meshed them using hexahedrons. Their length and width equaled those of the geometry, and their height equaled d/2.

Material modeling
To model the silicone mechanical properties, we used a second-order Ogden model of a custom DIW silicone whose parameters we obtained through mechanical testing of the bulk material. The effective Young's modulus was 3 MPa and the bulk modulus was 1 GPa. In tensile elongation-to-break experiments, the material failed at 5 MPa. We used Ansys' built-in material model for structural steel for the plates' material.

Compression setup
We used implicit mechanical models. For applications where compression speeds are low (mm/min), the systems can be described as quasi-static, making implicit modeling a feasible choice. Because, in general, implicit models take less time to converge on a solution, and we wanted to simulate many structures, implicit models fit our purposes well. However, some of the structures we simulated experienced local buckling, where explicit modeling could have worked better.
Because we used an implicit model, the results were time-independent, meaning that we could use an arbitrary time interval over which to simulate the compression. We set the compression to ramp at a constant rate from 0 to 50% strain over one second. As output, we gathered the reaction force values on the moving plate and von Mises stresses at the mesh nodes.

Experimental validation of models
Using a custom ink whose cured properties we applied to the simulations, we printed circular samples 50 mm in diameter with the lattice parameters outlined in Supplementary Table I. The material was formulated as a two-part (A:B) thermal cure Pt-catalyzed siloxane (10:1) having a vinyl terminated poly(dimethylsiloxane)-co-(diphenylsiloxane), reinforcing fumed silica, a thixotropic additive, a reaction inhibitor, a vinyl MQ siloxane resin, and poly(dimethylsiloxane)co(methylhydrosiloxane) as cross-linker. The samples were printed at room temperature using a custom DIW machine built on a three-axis Aerotech platform and were cured at 150°C.
We performed compression tests using an Instron 5969 with a 5 kN load cell. The upper platen was 50 mm in diameter with a spherical seat, and the lower, fixed platen was 28.68 mm in diameter. We performed compliance correction up to 1.5 kN prior to running any tests. For each sample, we performed three load-unload cycles up to 1.5 MPa with a compression rate of 1.26 mm/min and an unload speed of 2.54 mm/min. We collected displacement and reaction force data at 10-20 Hz.
To measure the printed lattice parameters and compare them to the nominal values, we cut square-shaped sections out of the printed samples (in a region that had not been compressed during testing) and imaged them using a Keyence VHX-6000 microscope.

Funding
This work was performed under the auspices of the US Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.

Data availability
The data sets generated during the current study are available from the corresponding author on reasonable request.

Code availability
The custom code used to compute L p is available from the corresponding author on reasonable request.

Conflict of interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.

Open access
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.