Tunable buckling configurations via in-plane periodicity in soft 3D-fiber composites: Simulations and experiments

We study the buckling of soft 3D-fiber composites (FCs) with varying in-plane microstructure periodicity. Through our experiments and simulations, we find that the out-of-plane buckling orientation of fibers is determined by the constituent material properties, volume fractions, and the in-plane periodicity. The examples are given for the FCs with rectangular in-plane arrangement of fibers, i.e., the fibers are periodically situated at distance 𝑎 from each other in one principal direction and distance 𝑏 along the other principal direction ( 𝑏 > 𝑎 ). We provide a buckling configuration map in the design-space of geometric parameters (fiber volume fraction and in-plane periodicity aspect ratio, 𝑏 ∕ 𝑎 ). Pre-determined by the microstructure parameters, the fibers buckle towards (i) the first principal direction along which the fibers are closer to each other, or (ii) the second principal direction, or (iii) towards a non-principal direction. Furthermore, we find that the characteristics of the buckling plane map is governed by the shear modulus contrast between the phases of FCs. We anticipate that the distinct instability patterns reported here will enrich the design-space for building deformation-controlled tunable materials.


Introduction
Soft fiber composites (FCs) simultaneously possess many desirable mechanical properties including flexibility, strength, toughness, rigidity, and fatigue resistance (Zhang et al., 2021;Wang et al., 2019;Xiang et al., 2020b), and are ubiquitously present in natural and biomaterials (Saheb and Jog, 1999;Humphrey, 2002). Their mechanical behavior can be significantly tailored by modifying the microstructure comprising of soft and stiff phases. The performance of soft FCs can be further enhanced by leveraging the pattern transformations associated with the buckling instability. Historically, the buckling of fibers was considered as a failure mode and, therefore, was to be predicted and avoided (Rosen, 1965;Christensen, 1979;Bazant, 2003). Recently, however, elastic instability phenomenon has been embraced to design tunable materials that can undergo sudden and dramatic microstructural transformations (Mullin et al., 2007;Galich et al., 2021;Florijn et al., 2016;Kochmann and Bertoldi, 2017;Babaee et al., 2016;Li et al., 2021Li et al., , 2022Greco et al., 2021). Furthermore, the recent advancements in the multi-material 3D printing (Zhou et al., 2020;Schwartz and Boydston, 2019;Xiang et al., 2020a;Kuang et al., 2019) N. Arora et al. the order of the initial microstructure. To predict the onset of instabilities, the linearized incremental analysis for small displacement superimposed on finite deformations is widely used (Ogden, 1997). The macroscopic instability can be estimated by performing the loss of ellipticity analysis on the effective tensor of elastic moduli, which can be calculated by employing micromechanics (Agoras et al., 2009;Rudykh and deBotton, 2012;Arora et al., 2020;Slesarenko et al., 2017b), phenomenological (Merodio and Ogden, 2002;Merodio and Pence, 2001;Ogden, 2003, 2005), or numerical-based approaches (Bruno et al., 2010;Greco and Luciano, 2011;Greco et al., 2020). The microscopic instability can be predicted by employing the Bloch-Floquet analysis superimposed on finite deformations (Triantafyllidis and Maker, 1985;Nestorovic and Triantafyllidis, 2004;Triantafyllidis et al., 2006). Moreover, the loss of ellipticity analysis is equivalent to the microscopic instability analysis in the long-wave limit (Geymonat et al., 1993).
In this work, we focus on the influence of in-plane periodic microstructure on the instability development in 3D soft FCs. In this regard, Slesarenko and Rudykh (2017) employed the Bloch-Floquet analysis to examine the instabilities and associated buckling patterns in periodic 3D FCs with the square arrangement of fibers. Recently, Galich et al. (2018) numerically studied the instabilities and propagation of waves in FCs with the rectangular arrangement of fibers (as shown in Fig. 1). By analyzing the polarization of eigenmodes in the vicinity of microscopic instability, they predicted that the orientation of the buckling plane in these composites can be switched by tuning the periodic unit cell's aspect ratio. Followed by this, Li et al. (2018a) experimentally realized the buckling behavior, however, only identifying a single buckling direction. In particular, they observed that the fibers always buckle in the direction they are closer to each other. Here, we find the configurations of FCs exhibiting change in the orientation of buckling plane, and also identify the factors influencing the behavior. As we shall show, in addition to the aspect ratio, the fiber volume fraction and the shear modulus contrast between the phases also significantly affect the direction of fiber buckling. Guided by the simulation results, we realize the instability-induced pattern formations and the rotation of buckling plane in the experiments.

Material fabrication and experiments
In this section, we provide the details of specimen fabrication, experimental device and setup. The 3D fiber composites are printed using the Stratasys Objet260 Connex 3D printer. Circular fibers are vertically aligned along the direction 2 and are periodically distributed in transparent soft matrix with the separation of distance along the 1 direction, and along the 3 direction, as per the dimensions of unit cell (Fig. 1a); note that > . The schematic of an example specimen is shown in Fig. 1b. The stiff fibers are printed using the digital material (mixture of TangoPlus and VeroBlack) DM-95 (shear modulus, ≈ 4.9 MPa Meng et al., 2020), while transparent soft matrix is printed using TangoPlus ( ≈ 0.21 MPa Slesarenko et al., 2017a;. Here, and thereafter, the fields and parameters corresponding to fibers and matrix are denoted as (•) and (•) , respectively. Geometrically, the FCs are defined by the aspect ratio of unit cell, ∕ , and the volume fraction of fiber phase, = 2 ∕(4 ), where is the diameter of each fiber (Fig. 1a). We print four different FCs with the geometric parameters, number of fibers, and dimensions provided in Table 1. To reduce the influence of boundary effects on the periodic FCs, the specimens are printed with an additional TangoPlus material layer of thickness 5 mm on the boundary. Uniaxial compression tests are performed using MTS Sintech 10/GL testing machine (maximum load 50 kN). The specimens are compressed quasi-statically along the direction of fibers, 2 , at the strain rate of 10 −3 s −1 . Since the soft matrix's materials is transparent, the embedded fibers and their buckling patterns are visible from the outside. The deformation of FCs is recorded using two high resolution cameras: The first camera is facing towards the ⟨ 1 , 2 ⟩ plane referred to as Plane 1 (Fig. 1c) and the second camera captures the deformation in plane ⟨ 3 , 2 ⟩ or Plane 2 (Fig. 1d).

Numerical modeling
In this section, we discuss the numerical analysis performed to predict the onset of instability and the associated buckling modes in FCs with rectangular arrangement of fibers. We employ Bloch-Floquet analysis superimposed on finite deformations to investigate the microscopic instability. We performed the analysis on the representative volume element (RVE) of the periodic FC with dimensions as shown in Fig. 1a; the height of the RVE is ℎ = 0.05 . The FC is subjected to the uniaxial deformation along the fibers defined by the following macroscopic deformation gradient: where is the applied macroscopic stretch along the direction of fibers and are basis vectors as shown in Fig. 1. Numerical simulations are conducted by the means of finite element code COMSOL Multiphysics 5.4. Fibers and matrix are defined by the neo-Hookean strain energy density function, which is integrated in COMSOL as where is the shear modulus and is the first Lame's parameter; 1 = tr( ) is the first invariant of the right Cauchy-Green deformation tensor, = ⊤ . Here, is the deformation gradient and ≡ det . To model the nearly-incompressible behavior of both phases, we set a very high ratio between the first Lame's parameter and shear modulus, ∕ = 1000. The instability analysis is performed in two steps: First, we apply macroscopic deformation (1) by imposing periodic boundary conditions on the faces of the unit cell. Second, Bloch-Floquet conditions are imposed on the unit cell as ( + ) = ( ) exp(− ⋅ ), where is position vector, is displacement vector, is the Bloch wave vector, and defines the spatial periodicity in the reference configuration. The eigenvalue problem is solved until a non-trivial zero eigenvalue is detected. The corresponding stretch ratio and wavenumber are identified as the critical stretch and critical wavenumber , respectively, and the orientation of buckling plane is determined by the polarization of the associated eigenmode. For a more detailed description of the numerical method implementation, see Slesarenko and Rudykh (2017) and Galich et al. (2018).

Results and discussion
We start by investigating the effect of the in-plane periodicity aspect ratio, ∕ , on the instability parameters and the orientation of buckling plane in the FCs. Fig. 2a and b show the critical stretch and normalized critical wavenumber ∕(2 ) as the functions of ∕ , respectively. The results are presented for the composites with fiber volume fraction = 0.05, ≈ 4.9 MPa, and ≈ 0.21 MPa, such that the shear modulus contrast is ∕ ≈ 23. The triangles with dotted curves represent the simulation results, and the circular symbols denote the experimental findings. The experimental value of critical wavenumber in the reference configuration is evaluated using the buckling wavelength in the current configuration. The relation between them is = 2 ∕ , where stretch ratio corresponds to the deformation level at which the wavelength is measured. Blue filled symbols mark the composites in which the fibers buckle in the principal Plane 1, whereas red hollow symbols show buckling in principal Plane 2. The green half-filled symbols and the estimated region marked in the green color denote the buckling of fibers in a non-principal plane. Fig. 2c shows the images of instability patterns observed experimentally: ∕ = 4 (top left), ∕ = 10 (top right), and ∕ = 7 (bottom). The horizontal bar is shaded in color to show the orientation of buckling at a given ∕ value: blue for Plane 1, red represents Plane 2, and the green region denotes buckling in a non-principal plane.
We observe that the FCs develop instability at smaller deformation levels ( increases) with an increase in the periodicity aspect ratio, ∕ . Moreover, the wavelength of instability patterns (inverse of ) increases as the periodicity aspect ratio increases. The values of the numerically predicted critical wavenumber values are in a good agreement with the experimental results. In experiments, however, the instability develops at higher deformation levels than predicted by simulations. The discrepancy observed here can be attributed to multiple factors, such as, for example, boundary effects or the stabilization effect due to the interphasial zone between the fibers and matrix (Arora et al., 2019). Nevertheless, the variation of critical parameters with aspect ratio follow the same trend in simulations and experiments.
Interestingly, the buckling plane of the FC rotates with a change in its aspect ratio. In particular, we observe that at smaller values of ∕ , the fibers buckle in the direction in which they are closer to each other, i.e., in Plane 1. For instance, at ∕ = 4, we find in our simulations and experiments, that the fibers undergo cooperative buckling in ⟨ 1 , 2 ⟩ principal plane (Fig. 2c-top left). On the other hand, at higher values of ∕ , Plane 2 is preferred over the first one for buckling. For example, in the FCs with ∕ = 10, the buckling is observed in the principal plane ⟨ 3 , 2 ⟩ ( Fig. 2c-top right).
Note that the orientation of the buckling plane does not discretely switch from Plane 1 to Plane 2. Instead, there exists a transition region (shaded in green color) where the fibers buckle in a non-principal plane. In simulations, this region can be identified by observing the polarization of the associated eigenmode(s) in the following two ways: (i) if the eigenmode for which the zero non-trivial eigenfrequency is obtained has a polarization along the non-principal direction, or (ii) for the eigenmodes with polarization in the principal directions ( 1 and 3 ), if the eigenvalues become (non-trivially) zero at almost same deformation levels. For the FCs considered here, the transition region approximately lies in the range 6 ≲ ∕ ≲ 7. In experiments, for the composite with aspect ratio ∕ = 7, the instability patterns are observed in both Plane 1 and Plane 2 (Fig. 2c-bottom), which indicates that the buckling occurs in a non-principal plane. This is in contrast to experimental results for the FCs with ∕ = 4 and ∕ = 10: for these composites the buckling is only observed in a single principal plane.
Next, we investigate the role of fiber composite microstructure geometry on the buckling plane orientation. To this end, Fig. 3 shows the buckling plane map as the function of geometric parameters: fiber volume fraction and aspect ratio ∕ , for the composites with ∕ = 23. The regions are color shaded according to the buckling plane orientations, similar to Fig. 2. The arrows in the schematic of respective regions show the direction of fiber buckling. The dotted curve separates the gray region that is geometrically inadmissible, given by the equation, ≥ ∕(4 ). The circular symbols represent the experimental results.
It is evident from Fig. 3 that the FCs can be classified into four geometric parameter groups, based on their buckling plane orientations. First, in the composites with small aspect ratios, i.e. with nearly square arrangement of fibers, the buckling plane orientates towards the N. Arora et al. non-principal direction (denoted by the green shaded region). These composites are not influenced by the unequal distribution of fibers along the two principal directions and their instability behavior is close to that of the square periodic FCs. With increase in the aspect ratio, however, the buckling plane switches to Plane 1. This transition, approximately marked by the dash-dotted curve, occurs at comparatively higher values of ∕ for the composites with smaller fiber volume fractions. For instance, in the FC with ∕ = 2, the fibers buckle along the non-principal direction when = 0.03, and along 1 -direction for = 0.05. In the composites that show the effects of unequal in-plane distribution of fibers, we observe that the volume fraction of fibers has a similar influence as that of periodicity aspect ratio on the buckling plane orientation. Specifically, in the FCs with smaller fiber volume fractions and smaller aspect ratios, the fibers are more likely to buckle along the direction in which they are closer to each other (or Plane 1). In contrast, the buckling along the second principal plane is preferred at higher values of and ∕ . The switch in buckling plane orientation may be attributed to the decrease in the distance between the periphery of fibers with an increase in and/or ∕ : This leads to a significantly higher reinforcement along the 1 -direction in comparison to the 3 -direction and the FCs tend to behave as laminates with layers periodically distributed along 3 -direction. Hence, ⟨ 3 , 2 ⟩ plane (or Plane 2) becomes energetically favorable for the instability development. These two regions (shaded in red and blue) are separated by a transition region in the middle (shaded in green color). In this region, similar to the nearly square periodic FCs, the fibers buckle in neither Plane 1 nor Plane 2, instead buckle in a non-principal plane. The black solid curve in Fig. 3 is approximately drawn at the center of the transition region; it is termed as the transition curve.
The experimental results for FCs with = 0.05 having aspect ratios ∕ = 4, ∕ = 7, and ∕ = 10 are in good agreement with the simulations (as was also previously shown in Fig. 2). Moreover, the fibers in FC with = 0.06 and ∕ = 10 are experimentally observed to buckle in Plane 2 (hollow red symbol), further validating the numerical predictions.
Finally, we examine the dependence of fiber buckling orientation on the shear modulus contrast between the phases. In Fig. 4, we plot the transition curves, which correspond to the transition region between Plane 1 and Plane 2 on the buckling plane map (see the black solid curve in Fig. 3). We consider the FCs with shear modulus contrast: ∕ = 10 (blue), ∕ = 15 (red), ∕ = 23 (green), and ∕ = 50 (orange). We observe that with an increase in the shear modulus contrast ( ∕ ), the transition curve shifts towards the N. Arora et al.  higher periodicity aspect ratios and higher fiber volume fractions. For instance, consider the FCs with geometric parameters = 0.03 and ∕ = 8, they develop buckling in Plane 2 for ∕ = 10, whereas fibers buckle in Plane 1 when the shear modulus contrast is increased to ∕ = 50. Recall that the area between the transition curve and the dotted curve (separating the geometrical inadmissible region) correspond to the FCs that develop buckling in Plane 2. Therefore, as the transition curve shifts with an increase in the shear modulus contrast, FCs are more likely to develop buckling in the first principal plane for a wide range of geometric parameters. Hence, higher values of ∕ favors the buckling of fibers in Plane 1.

Conclusion
In summary, we have experimentally and numerically investigated the instability development and associated buckling patterns in 3D FCs with rectangular in-plane arrangement of fibers. We have observed that the composites with higher values of the periodicity aspect ratio develop instability at smaller deformation levels and their buckling patterns are characterized by longer wavelengths. We have provided a buckling orientation map dividing the microstructural geometric parameters into four groups. First, the FCs with small aspect ratios almost behave like square periodic ones, and their fibers buckle in the nonprincipal plane. The increase in the aspect ratio and/or fiber volume fraction, however, leads to the switch in the buckling plane towards the principal direction in which the fibers are closer to each other (or Plane 1). At very high values of these geometric parameters, the FCs are more likely to develop instability in the second principal plane (or Plane 2). Lastly, we have found a narrow parameter-space, transition region, between the Plane 1 and Plane 2 regions. In the transition region, similar to nearly square periodic FCs, the fibers buckle towards a non-principal direction. Moreover, the buckling plane orientation of composites with identical geometric parameters can vary from Plane 2 to Plane 1 with an increase in the shear modulus contrast. This dependence on the stiffness of phases may be further utilized to enrich the design-space through, for example, introducing viscoelastic phases; such inelastic behavior may allow rotation of the buckling plane based on the loading history.
The results presented here can be utilized in designing systems that exploits the targeted buckling plane; for instance, an acoustic metamaterial that can filter waves based on their polarization. Moreover, these controlled microstructure transformations in FCs can be merged with other material systems at different length scales to build hierarchical elastic metamaterials (Hussein, 2019;Bacquet et al., 2018;Arora et al., 2022;Ostoja-Starzewski et al., 2016;. The dynamic properties of such systems can be studied by applying the recently developed adaptive isogeometric methods (Yu et al., 2022;Jansari et al., 2022;Chen et al., 2022Chen et al., , 2020. Furthermore, we note that the onset of instabilities can be affected by the uncertainties or imperfections in the material or geometry of the composite. The effect of such uncertainties can be quantified and integrated into the numerical framework through stochastic analysis (Ding et al., , 2021Hauseux et al., 2017Hauseux et al., , 2018.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.