Post-buckling behaviors of thin-film soft-substrate bilayers with finite-thickness substrate

Surface buckling behaviors of thin-film soft-substrate bilayers have important research value. Recent research has focused on bilayers with infinite-thickness substrates. However, bilayers with finite-thickness substrates widely exist. To study this problem more comprehensively, we extended the stability theory of a beam on an elastic foundation to bilayers and then established a finite element method of bilayers using the neo-Hookean hyperelastic constitutive model. A self-contact analysis method was coupled to the finite element method so that the further buckling evolution of the film surface after folding could be simulated. Based on our analysis of various modulus ratios and thickness ratios, the evolution of the buckling path was significantly influenced by the thickness ratio. Without considering the situation of a prestressed substrate, four new buckling paths were found. Thus, we extended the single buckling path (under infinite thickness substrate) to five types. Our study also found that for path four, the substrate with a certain thickness exhibited a special final stable surface morphology. That is, regardless of the friction, a new periodic morphology after film folding appeared due to the contact slip of the film surface. Finally, further analysis showed that these five buckling paths were all dependent on different modulus ratios and thickness ratios.


Scientific Reports
| (2022) 12:4074 | https://doi.org/10.1038/s41598-022-08136-w www.nature.com/scientificreports/ hyperelastic constitutive model to simulate the buckling behaviors under various modulus ratios and thickness ratios. We also coupled the self-contact analysis method to our finite element model, so that the further buckling behaviors after folding could be determined. All of the simulation examples were calculated for surface contacted slip after folding, and the final stage of the buckling path was considered to be symmetry breaking of the surface morphology of the film. We first compared our analytical solution and the finite element model results for film wrinkling to verify the accuracy of our finite element model. The change of the buckling path of the film surface with the thickness ratio was then studied by more post-buckling simulations. The results showed that the buckling path of the film surface was significantly influenced by the thickness ratio, and with the increase in the thickness ratio, five type buckling paths were found. The typical buckling path was path five, corresponding to an infinite-thickness substrate. This was consistent with the results of Xu 16 , and we further compared with other previously reported experimental data 15 . Path one was a local ridged buckling evolution, which was similar to the surface ridge for prestressed substrates reported previously 13,30 . However, the local ridge of path one with a very small thickness ratio and very large modulus ratio exhibited a different behavior, i.e., with increasing load, the local ridge of the film surface stopped increasing. In contrast, the adjacent wave crest began to increase and continued until all of the wave crests increased. Finally, the surface morphology became a sinusoidal periodic mode, similar to wrinkling. This mode was similar to the telephone cord buckling reported previously [31][32][33] , in which the variations of the buckling behaviors of annealed telephone cord films under high compressive stresses with temperature were reported. Path two also appeared when the substrate was very thin, but unlike path one, after slightly increasing the substrate thickness, the surface of the film retained a sinusoidal periodic mode after wrinkling, and new buckling modes did not appear. Path three involved a double period mode after wrinkling in the film surface, but unlike path five, with increasing load, the period quadrupling did not appear, but presented a new buckling mode, i.e., the upward wave trough of the formed double period mode disappeared, and a new morphology of the film surface that resembled a double folding appeared. In path four, the complete buckling evolution of Wrinkling-period doubling-period quadrupling-quadruple folding occurred. Upon further loading, the quadruple folding diminished due to contact slip of the film surface, and after a period of asymmetric evolution, a stable and symmetric triple folding appeared. This final buckling stage of path four was similar to that reported previously 30,34 , in which the triple folding of bilayers occurred under substrate prestressing. All of the paths are shown in Fig. 1.
Finally, the buckling path changes with different modulus ratios and thickness ratios was given in our paper.

Wrinkling mode of bilayers in small deformation
The thin-film soft-substrate bilayers can be regarded as an incompressible film attached to an elastic half space, and we can use the traditional elastic space method to solve this problem. We assume that the film is incompressible with a length of 2L and a thickness of h, and it is attached to an elastic substrate of thickness H. The film and substrate do not undergo any relative sliding and separation (Fig. 2).  The deflection of thin film is denoted as w(x), and we assumed that the wrinkling of the bulk structure was periodic with period l 0 under a symmetric displacement load ε n at both ends. The load ε n satisfies the following relations: The force of the elastic substrate acting on the film was assumed to be P y (any point P y along the vertical direction). We considered the differential equation of the deflection of a beam on an elastic foundation. In the case of plane strain, this differential equation is as follows: where D and Q are the bending stiffness and in-plane load, respectively, for a rectangular section, , the ε is the in-membrane compressive strain, E f is the film elastic modulus, and µ f is the film Poisson's ratio. P y can be obtained by coordinating the deformation of the film and substrate. We consider the situation of a small deformation. According to linear elasticity theory, the displacement in two directions, u(x, y), v(x, y), of the substrate is given as follows 15 : where µ s is the substrate Poisson ratio.
The substrate stress σ i (x, y) (i = x, y) satisfies the geometric and physical equations, and we obtain the following: where the E s is substrate elastic modulus.
At y = 0 , the vertical displacement f(x) of the substrate is coordinated with the film deflection, i.e., v(x, 0) = f (x) = w(x) , and other boundary conditions can be set to u( To solve Eq. (3), we use the Fourier transform. The substrate displacements u, v were transformed along the x-axis, yielding where ũ and ṽ are the frequency transforms of u and v, respectively, and ω is the Fourier parameter. The frequency transforms of Eq. (3) are as follows: (2) Dw (4) + Qw (2) + P y = 0, where The solution is the displacement of the finite-thickness substrate. In the general case, we can obtain an accurate analytical solution of any initial buckling mode using Eq. (9). We next solve Eq. (2) in the Fourier space. The expression is rewritten as where F P y =σ y (x, 0) is the Fourier transform of the substrate stress function. By Eq. (4), we obtain Using Eq. (9), the effect of the substrate F P y is All the effects of the substrate constitutive relations and boundary conditions are included in the function �(ω) , as follows: Using Eqs. (9), (11), and (13), we obtain the governing equation of the frequency space: In the case of a small deformation, a presumptive displacement function f (ω) is included, and the vicinity of the initial buckling solution can be obtained by Eq. (15). To study this problem more comprehensively, we used a function Z to represent Eq. (15): We define an intermediate function g(ω) to determine the relationship of the buckling parameters. g(ω) is assumed to be Using the Eq.(16) and Eq. (17), combined with integral operation, easy to get .  (18) can be simplified, and the following compressive strain is obtained: . The right side of Eq. (19) will no longer be related to parameter ω after the integral operation. We now analyze the wrinkling mode of the film. We assumed the Poisson's ratios µ f and µ s were both 0.433 15 to simulate a soft material. Meanwhile, we were more interested in the relationship between the buckling mode and the thickness ratio. Thus, the modulus ratio was assumed to 600. We selected the displacement functions f (x) = a cos(kx) and f (ω) = π(δ(ω + k) + δ(ω − k)) , we set the critical compression strain ε 0 of the initial buckling stage to an extreme value of ε , and the initial wave number set to k 0 .
This coincided with the experimental observation that the film surface morphology retained a sinusoidal period during surface wrinkling. Further results of Eq. (19) for different thickness ratios are given in Fig. 3, which shows that the compression strain of wrinkling changed with the thickness ratio. The case with H h = 1 often occurs when the substrate is attached to a harder object, such as bones or wires 26,27 . Meanwhile, with the increase in H h , the extreme point of the curve moved downward such that the critical wavenumber and strain were reduced. When H h = 50 , further increases in H h produced negligible changes of the curves, so the substrate could be regarded as an infinite-thickness substrate. The extreme point of the curve provides the critical compression strain. In the general case, the thinner the substrate is, the smaller the critical wave number becomes, and the amplitude from Eq. (1) is larger.

Finite element model of thin-film soft-substrate bilayers
The initial buckling stage of the thin-film soft-substrate bilayers appeared as a small deformation. Although we have given an analytical solution, the post-buckling of large deformation is still difficult to solve. Brau et al. 15 analyzed the infinite substrate buckling behavior used the Green strain. A governing equation of large displacement was provided, the PDEs were solved by a perturbation method and verified with experimental data, and an approximate method to predict the film buckling of wrinkling and period doubling was presented. However, due to the high nonlinearity of the problem, this method cannot incorporate the hyperelastic constitutive relation of a soft material. We consider the post-buckling of a finite substrate. It is difficult to obtain a mathematical solution, but the finite element method (FEM) provides good approximations 16,23,35 .
Based on the powerful FEM, to obtain a more profound understanding of the physical process of film deformation, we considered the contact-sliding of a film after the film was folding. This process used the augmented Lagrangian method based on a penalty function 36,37 that is widely used in finite element analysis of adhesive friction contact, sliding friction contact, and frictionless contact problems. To ensure the correctness of the numerical simulation, a hexahedral mesh was selected, the element size of the film was set to 0.4h along x-axis and 0.2h along the y-axis. For the substrate, the element size was same as film thickness along the x-axis, and the element size of the y-direction was gradually increased from 0.2 h to 8 h. The film and substrate were compatibly deformed.
The two-sided displacement load of the finite element model was applied using the rigid columns shown in Fig. 4. The rigid column and model were set with a contact element with zero friction (to speed up convergence, the frictional force of the film and the rigid column can be set to smaller values initially and then set to zero after Figure 3. Relationship between the compression strain and the wavenumber given by Eq. (19). The modulus ratio was assumed to be 600, and µ f and µ s were assumed to be 0.433. The different curves represent the different thickness ratios. www.nature.com/scientificreports/ wrinkling, but it is not mandatory). We considered the contact slip of the film surface, a self-contact element was set on top of film. The y-direction displacement of the element model was limited at the bottom. Meanwhile a displacement load was applied to the rigid columns. Otherwise, the large deformation of the buckling was considered. We used the Solid 185 element of the ANSYS software and the Neo-Hookeam hyperelastic constitutive model. All the simulations were completed using the ANSYS software, and part of the theoretical calculations were aided by the MAPLE software. We adopted the default convergence criterion of ANSYS to ensure the reliability of the calculation results.

Initial stage buckling of analytical solution and FEM.
We investigated the wrinkling of the finite element model under small loads. In this case, the Neo-Hookean hyperelastic constitutive model was similar to a linear elastic constitutive model, allowing us to compare our finite element model with the analytical solution proposed in Section 2. We considered different modulus ratios and thickness ratios, critical wavelengths 0 , and critical loads ε n0 of wrinkling obtained by the finite element model and Eq. (19), as shown in Fig. 5. When the thickness ratio was large, the FEM and theoretical solutions of the wrinkling wavelength agreed closely. However, when the substrate thickness was reduced, our FEM solution was larger than that of Eq. (19), which may possibly because of that there is a competition between sinusoidal mode and local ridge mode. In the case of thin substrate, the influence of local ridge mode is amplified synchronously. So, with the comprehensive influence, although the film is still wrinkled in a sinusoidal mode, there are still systematic differences of the wavelength between the idealized theoretical solution and the FEM simulation. For the case of a thick substrate, the FEM simulation is in good agreement with the theoretical results. The calculations also show that the FEM  Buckling path of local ridge. Due to physical constraint of the substrate in a bilayer system, the inital buckling mode shows a sinusoidal period mode. The film surface morphology also generally exhibited a sinusoidal period after wrinkling. However, when the modulus ratio was very large and the thickness ratio was very small, our simulation showed a different buckling evolution in which the sinusoidal period mode became a local ridged mode. With increasing load, the adjacent wave crests began to increase until all of the wave crests increased, and the surface of the film exhibited a sinusoidal period mode again. Based on the comparison of the FEM solution and previously reported results 32 , for the conditions in our FEM simulations, the elastic contribution of the substrate can be neglected due to the very large modulus ratio and very small thickness ratio. However, the bottom edge restriction cannot be ignored. Thus, the buckling path of the bilayers in this case is similar to the buckling of a film in contact with an air layer and fixed on a stiff substrate. Our FEM simulation results were also similar to the experimental results reported by 32 . In their experiment, the confined air layer enclosed in the thin layer can be assumed as an elastic substrate with very small stiffness (the elastic stiffness here is very small, and the influence of shear stiffness is ignored). Accordingly, the antisymmetric axial plane of the warped thin layer is similar to the fixed lower surface of our finite element model. So, the thin layer and the air layer of the experiment can form a bilayers with small thickness ratio and large modulus ratio as shown in Fig. 7. The FEM solution of thin-film soft-substrate bilayers with a modulus ratio of 600 and a thickness ratio of 5 was analyzed. The buckling path was wrinkling, followed by local ridge formation, and finally, sinusoidal period structure formation. As shown in Fig. 7b and c, our FEM simulation results were similar to previously reported experimental results 32 . We selected a symmetric model and observed the amplitude change of seven wave crests (for dimensionless purposes, A i / 0 was chosen), and all of the solutions were entered and shown in Fig. 7d (the several lines prior to bifurcation in the figure are overlapping, and amplitudes are following the same paths).
Buckling path of bilayers with large thickness ratio. We next analyzed the post-buckling of thin-film soft-substrate bilayers. To compare with experimental data 15 , we set the modulus ratio to 120, and the thickness ratio was 60. The simulation results are shown in Fig. 8.
The observed amplitude changes of A 1 , A 2 , and A 3 are shown in Fig. 9 as the load was increased (the several lines prior to bifurcation in the figure are overlapping, and amplitudes are following the same paths). The comparison shows that our FEM method provides a suitable method for predicting the wrinkling-period doubling mode. As for the quadruple period, there were some errors between the FEM solution and the experimental data. However, the simulation results using a lattice model method showed the same error trend as that of previously reported data 16 .
The buckling path shown in Fig. 8 is the typical evolution of the thin film with a very large substrate thickness. To distinguish this from other paths, it is denoted as path five. In this path, the final buckling stage was a quadruple folding of the film surface. A further increased load could not form any new symmetric periods but broke the quadrupling mode, causing the surface morphology to be disordered.
We next reduced the thickness ratio, and our FEM simulation gave the path four mode shown in Fig. 10. In this example, the thickness ratio was set to 45, and the modulus ratio was set to 400. The result of our simulation revealed a new morphological transformation, i.e., the quadruple folding to the triple folding (Fig. 11). www.nature.com/scientificreports/ Paths two and three are formed by further reducing the substrate thickness. To simulate these situations, we selected example thickness ratios of 20 and 10 with a modulus ratio of 240. Figure 12 shows path three for a thickness ratio of 20. The upward wave trough after the period doubling disappeared but formed a double folding mode.
We are most interested in the amplitude changes of A 1 and A 2 in Fig. 12. Similarly, Fig. 13 shows the FEM simulation of path two, in which the film retained a sinusoidal period structure after wrinkling . The amplitude changes of paths two and three are shown in Fig. 14a and b (the several lines prior to bifurcation in the figure are overlapping, and amplitudes are following the same paths), respectively.

Results and analysis.
We analyzed the FEM simulations with different modulus ratios and thickness ratios, and all of the results are presented as the two-dimensional image shown in Fig. 15. Paths one to five show the buckling evolution of the thin-film soft-substrate bilayers, and different modulus ratios and thickness ratios produced different buckling paths.

Concluding remarks
This paper considers the actual physical limitations of thin-film soft-substrate bilayers with the finite-thickness substrates. The buckling evolution path simulations of the bilayers were simulated by the FEM method with different modulus ratios and thickness ratios with compression loads imposed on both sides. We also extended the buckling theoretical analysis to bilayers with finite thicknesses. The general solution of an arbitrary surface deflection with a linear elasticity model was given to verify the correctness of the FEM method. The post-buckling behavior of large deformation was simulated by our FEM model.  www.nature.com/scientificreports/ Our calculations showed that the film buckling went through five paths with the change of the thickness ratio of the film and substrate. We hope our results will play a role in predicting or adjusting the film surface morphology after buckling.
Finally, there is some error between our FEM results and the experimental data in the final buckling stage, which may be related to the selection of our constitutive relationship. Perhaps a new constitutive equation such as consider of the strain stiffening model of Gent or Ogden may be contribute to reduce the error of simulations.