Toward new scaling laws for wrinkling in biologically relevant fiber-reinforced bilayers Open Access

Corrugation is a very common geometrical feature in nature. This is indeed the case for skin, blood vessel walls, the brain, etc. (see, e.g., Genzer and Groenewold, 2006;Hohlfeld and Mahadevan, 2011;Budday , 2017; and Holland , 2020 among many others).

For instance, as pointed out in Nguyen (2020), a wide number of papers regarding wrinkling, period doubling and quadrupling, creasing, and folding in biological systems, including tissues such as ant’s eyes (see, Fig. 1), have been produced in the last two decades (see, e.g., Genzer and Groenewold, 2006;Ciarletta and Ben Amar, 2012;Ben Amar and Jia, 2013;Ciarletta , 2014;Balbi , 2015;Mottahedi and Han, 2016;Goriely, 2017; Alawiye , 2019; 2020; Nath , 2020;Chen , 2021;Kai , 2022; and Mostafavi Yazdi and Baqersad, 2022 and references cited therein). Furthermore, important results regarding various thin man-made mechanical systems exhibiting corrugation have been largely investigated in parallel (see, e.g., Biot, 1963;Cerda and Mahadevan, 2003;Pocivavsek , 2008; and Cutolo , 2020 and references cited therein, among many others).

Most of the investigations mentioned above have dealt with homogeneous (hyper-) elastic bilayers, perfectly bonded to one another. Those studies have been performed primarily under either an applied prestretch or compressive in-plane external tractions. Occasionally, thermal actions or growth (with reference to biological systems) have also been analyzed as a source of possible instability through corrugation, although not so extensively. For the given action, the aforementioned literature shows that the enabling features for wrinkling are (i) the extreme thinness of one of such layers relative to the thickness of the whole system, and (ii) the mismatch of the elastic moduli of such layers.

Unlike other phenomena, though, very few scaling laws connecting the geometrical features of the exhibited corrugations and the mechanical properties of the bilayers described above are available. In particular, in Allen (1969) (Sec. 8.2), a slightly modified version of the scaling laws (24) and (25), displayed in the sequel, governing the critical strain and the wavenumber at the onset of wrinkling, were obtained in a fairly simple and clever way for an elastic strut bonded to an isotropic elastic core. Such laws have been revisited in more recent times by Sun (2011), where those relationships have been obtained (without showing the actual derivation) as asymptotic expansions of the analytic solutions of the wrinkling problem for isotropic hyperelastic bilayers. With regard to a completely different situation, such as free-standing thin polymeric sheets under tension, a new set of scaling laws has been provided in Cerda and Mahadevan (2003) and analytically validated (with a slight change) in Puntel (2011). More recently, in Goriely and Mihai (2021), generalizations of (24) and (25) were obtained for bilayers made of liquid crystal elastomers (with certain given initial orientations of the domains) bonded with a homogeneous and hyperelastic neo-Hookean material, both in the case in which the thin layer is neo-Hookean and the substrate is made of the liquid crystal elastomer and vice versa.

Primarily due to the presence of fibers, scaling laws for wrinkling occurring in biological tissues are not yet available in the literature. Nonetheless, an ad hoc equation has been recently provided in Nguyen (2020) for the critical strain at the onset of the instability, although it did not come from any mathematical justification.

Among other issues, the main problem of soft biological tissues is certainly heterogeneity. This can influence the mechanical response of the tissue in terms of inhomogeneity of its pointwise elastic properties and, depending on the shape and functionality of the tissue, its residual stresses (see, e.g., Taber and Humphrey, 2001 and Hayn , 2020, and references cited therein). Nevertheless, averaging and homogenization methods yielding effective mechanical properties for biological tissues have been developed over the last two decades (for a more detailed discussion, see, e.g., Robertson and Watton, 2013;Cyron , 2016; and Braeu , 2017). This leads to overall characterizations of the constitutive behavior of such complex systems (see, e.g., Humphrey and Rajagopal, 2002;Gasser , 2005; and Bellini , 2013 and references cited therein) for which the degree of approximation can, of course, vary significantly (see, e.g., Robertson and Watton, 2013 for a detailed discussion of this aspect). The most utilized approach for such tissues, with particular regard to arterial walls, is certainly the one introduced in Holzapfel (2000) (called OGH in the sequel). Such a constitutive equation has been utilized in Nguyen (2020) and it will also be employed in the sequel together with the Standard Reinforcing Model (called SRM in the sequel). The latter has been studied since the eighties (see, e.g., Kurashige, 1981 and Triantafyllidis and Abeyaratne, 1983), although it was later in Qiu and Pence (1997) that the impact of such a constitutive law on the deformation modes exhibited by this material was investigated. More recently in Melnik (2015) and Sen (2022), the SRM law has also been exploited in relation to the dispersion of the fibers.

The present paper is the first step toward finding a rigorous procedure enabling one to systematically find scaling laws for corrugation starting from the equations governing such a phenomenon. In particular, the work here is organized as follows. In Sec. II the approach undertaken in Nguyen (2020) for the study of wrinkling in bilayers formed by a three-dimensional stiff thin film adhering on top of an OGH (and then SRM) fiber-reinforced, and much softer and thicker, substrate is revisited through a simplified approach. Here, instead of treating the top layer as a three-dimensional solid, a dimensionally reduced formulation [like the plate model in Shield (1994)] is assumed, and the simplified constitutive SRM law for the fiber-reinforced substrate is considered.

In Sec. III, a comparison of the outcomes of choosing SRM instead of the more complex OGH law is performed. Indeed, such a comparison is produced for the analytic results for both the critical strain and the wavenumber at the onset of wrinkling coming from the OGH constitutive law both by considering the top layer as (i) a three-dimensional solid and (ii) as a plate, and (iii) the SRM law for such a dimensionally reduced approach.

Furthermore, asymptotic expansions for both the wrinkling strain and the corresponding wavenumber have been provided for high-contrast elastic mismatches between the thin layer and the substrate in the presence of the reinforcing fibers. This starts from the outcome of the analytic procedure performed to seek the (a) minimum critical strain with respect to the wavenumber among the ones solving the eigenvalue problem characterizing the balance of forces at the interface between film and substrate, and (b) the corresponding wavenumber. The latter is then processed through a suitable sequence of Taylor’s expansions yielding (19), a novel scaling law for the wavenumber itself formed by a product of two terms, a basal one and an amplifying factor. The former term turns out to coincide with (25), namely, the asymptotic law for the wavenumber of a purely neo-Hookean bilayer reported in Sun (2011) and Cao and Hutchinson (2012). In the cases of either the absence of the fibers or their perfect randomness, the amplifying factor goes to one, thereby letting the novel scaling law for the wavenumber degenerate to (25). There, the wavenumber scales like the cubic root of the elastic mismatch of the two layers forming the system in that case. Full novelty is instead in the amplifying factor (26) due to the presence of load-bearing distributed fibers within the matrix of the substrate. That factor turns out to scale with the sixth root of a sum of terms. The latter turns out to be even in the spatial dispersion of the fibers (up to the fourth power of that parameter), and modulated by suitable powers of the modified stiffness ratio between fibers and matrix (accounting for the volume concentration of the former), and on the square of the sin of four times the relative orientation of the fibers themselves. With an analog procedure, the novel scaling law (23) for the critical strain at the onset of wrinkling is also obtained. Not surprisingly, this retrieves (25) (see Sun , 2011 and Cao and Hutchinson, 2012) for isotropic neo-Hookean bilayers, either when fibers are absent or whenever they are randomly distributed. In all of the other cases, the modulating function arising in (23) depends on the presence of the fibers and it is nothing but the square of the amplifying factor previously obtained for the wavenumber. In this same section, diagrams showing the comparisons between the obtained scalings, the analytic results obtained in the previous section, and numerical results performed by using ABAQUS for finite element method (FEM) simulations have been displayed. Such figures relate to results for high elastic contrast between the top thin layer and the substrate and given sets of parameters, carefully discussed in Sec. III.

Finally, in Sec. IV, a re-interpretation of the obtained asymptotic expansions for both the critical strain and the associated wavenumber is proposed in terms of the resulting properties of the linearized system about the underformed state obtained in Nguyen (2020). It is worth recalling that the result of such a linearization yields an actual orthotropic material response for the substrate. This innovative way of looking at the newly derived scaling laws illustrates how the modulating function mentioned above is essentially related to the orthotropy of the linearized solid. Indeed, the modulating factor introduced above goes with the sixth root of a term governed by the ratio of the Young moduli evaluated in the principal system of the resulting linearized orthotropic medium, while still depending on the square of the sin of four times the relative angle between the family of the reinforcing fibers.

In Nguyen (2020), an approach to computing the critical strain for which a thin membrane adhering to a soft substrate experiences wrinkling is presented. In that paper, the computation of such a strain is performed by considering the system as composed of two three-dimensional solids and then writing appropriate plane strain balance equations. However, this approach has the computational disadvantage of solving a highly non-linear system. In order to circumvent this drawback, the geometry and the physics of the problem suggest key simplifying assumptions leading, in a much simpler way, to almost the same results obtained from the fully three-dimensional model cited above.

A more efficient approach can be undertaken by focusing the present analysis on:

• bilayers for which the mismatch between the elastic moduli of the layer and of the substrate is very high (i.e., between $10 4 and 10 6$⁠);

• the layer being considered as very thin compared to the substrate (which, in mathematical terms, is in fact assumed infinitely deep).

The outcomes of a dimensionally reduced theory for the top layer, such as the plate one adopted here, and its interactions with an OGH infinite layer have not yet been explored in the literature. Indeed, in the aforementioned recent paper by Nguyen (2020), both the layer and the OGH substrate were treated as fully three-dimensional bodies under plane-strain conditions. No matter the constitutive response of both layers nor how the film is modeled, the balance of tractions at the interface between film and substrate governs the configurations of the bilayer.

The incompressibility condition $det F=1$ is satisfied at the first order, i.e., $λ 1 − 2∂ u 1 ~/∂ X 1+∂ u 2 ~/∂ X 2=0$ (see e.g., Pence and Song, 1991 and Yue , 1994). This last equation is identically satisfied by (3) if $u ~ i= x i− λ i X i(i=1,2)$ is assumed.

In order to facilitate the reader, the following notation is utilized in the sequel: M stands for “matrix,” F for “fibers,” and L for “layer.” In addition, the order reflects the position of the shear modulus of a given system within the ratio: for example, $ρ M L$ means stiffness of the matrix (M) forming the substrate over the one of the layer (L).

Due to the complexity of the OGH model, it is worthy of mention that the amount of calculations required to solve (14) significantly increases relative to the case of neo-Hookean bilayers. Hence, a simpler model than OGH would be especially useful if it would deliver comparable results to its associated optimality conditions.

Solutions of (14) obtained by utilizing the SRM strain density energy are shown in Fig. 3. From there, it is manifest that the critical strains are practically the same by using either constitutive equation, thereby suggesting that the assumption of a simpler constitutive law, such as the SRM, leads indeed to comparable results. This outcome is related to the independence of SRM from the parameter $k 2$⁠. As a confirmation of this circumstance, Nguyen (2020) illustrated that the OGH law does not depend on $k 2$ in the small strain regime: in the sequel, (see Figs. 3–6) the magnitude of the arising strains are shown to be small enough.

A comparison between the outcomes of (i) the dimensionally reduced model coupled with SRM for the substrate and (ii) the solid one developed by Nguyen (2020), is shown in Fig. 3. There, the critical strains and non-dimensional wavenumbers are displayed as functions of the angle formed by the two families of fibers (displayed in Fig. 2) and for different stiffness ratio $ρ F M$⁠. This has been done by setting, as in Nguyen (2020), $k 2=0.8393$ and $κ=0$⁠. Moreover, $ν L=0.5$ has been imposed for the incompressibility of the layer and a stiffness ratio $ρ M L= 10 − 4$ between the matrix and the layer has been assumed. The diagrams show that the results are essentially the same as the original model for exceptionally small stiffness ratios. Furthermore, a symmetric behavior for both the critical strain and the corresponding wavenumber about $θ= 45 °$ is detected in the assumed range $10 − 6− 10 − 4$ for $ρ M L$⁠, for the considered values of $ρ F M$ and $κ$⁠, namely, the modified stiffness ratio between fibers and matrix (accounting for the volume concentration of the former) and the spatial dispersion of the fibers themselves.

The solution given by the system (14) is not generally available in a simple form and, for a given set of the model parameters, it is therefore necessary to solve it numerically. However, it is still possible to simplify its expression under the following assumptions:

• the upper layer is considerably stiffer than the substrate, namely $ρ M L=6 μ M/ E L≤ 10 − 4$⁠;

• the fibers ratio $ρ F M= k 1/ μ M$ assumes values between 0 and 10. This hypothesis, although it may appear limiting as it would bring back to what has been already obtained through the neo-Hookean model, produces reliable results when (a) is fulfilled;

• the critical strain is very small (⁠ $ε c r≤ 10 − 3$⁠, see e.g., Sun , 2011 in the absence of fibers) and therefore, as $ε c r=1− λ 1 c r$⁠, $λ 1 c r≈1$⁠. Because no prestretch is considered, the eigensolutions given by (7) are linearized around $λ 1 c r=1$⁠. In this way, the resulting quantities, namely the critical strain and non-dimensional wavenumber, will have no dependence on the stretch;

• the critical strain and non-dimensional wavenumber are approximately constant when $k 2$ changes (at least for $ρ M L= 10 − 4$⁠), as remarked in the previous section. It is thereby possible to replace the contribution of fibers given by the OGH constitutive law (1) with the SRM (15), removing the variable $k 2$ and the whole exponential part associated to that;

• the Poisson ratio of the layer is $ν L= 1 2$⁠.

By analogy with Allen (1969) and Sun (2011), as pointed out in the previous section, explicit solutions to the optimality problem are sought. In other words, reliable asymptotic expansions for (i) the minimum of the critical strain yielding the onset of wrinkling and (ii) its corresponding wavenumber are the targets of this section.

The same plots display the results coming from the utilized FEM. The numerical simulations have been performed by means of the commercial software ABAQUS/Standard (Lic. n. LKO2211177).

The substrate is modeled as a three-dimensional hyperelastic body under plane strain conditions. In particular, within the plane of interest, a rectangle of length equal to ten times the wavelength $λ=2πh/ k h$ [provided by the asymptotic expansion (27)] is considered. For the sake of computation, the depth of the substrate is taken as the maximum between one hundred times the thickness of the top layer and two times the wavelength. This choice essentially provides a semi-infinite substrate compared to the thin layer, and it relies on a theoretical justification thanks to Pence and Song (1991) and Yue (1994) (see the barreling curves reported in such papers in Fig. 4, same numbering in each paper).

Based on the geometry described above, it is clear that the smaller the ratio $ρ M L$ is, the bigger the size of the substrate will be. Upon utilizing two-dimensional plane strain solid elements, the thickness of the layer governs the characteristic size of the mesh (that is to avoid distorted mesh that could lead to numerical ill-posedness of the governing operator). As a result, the simulations with small $ρ M L$ will be penalized from a computational point of view, due to an extremely high number of nodes.

Therefore, the following choices have been adopted for the discretization and computational analysis of the system. As pointed out in Sun (2011), instead of modeling the upper layer on the basis of a two-dimensional geometry, the thinness of the top layer relative to its length suggests the use of modified B22 beam elements, by means of the built-in stringer option in ABAQUS. The elements just mentioned above are indeed properly modified to represent the plate behavior of the layer undergoing plane strain conditions. The latter evidently constrains the lateral contraction/expansion of each stripe of the top film during wrinkling. Henceforth, the representative cross section of the thin layer can be taken with unit height, whereas its base must be set equal to $( 1 − ν L 2 ) − 1$⁠. Regarding the substrate, eight-node, hybrid, plane strain elements CPE8H are adopted. In order to properly display the wrinkling mode, the mesh size is calibrated at about $λ c r/40$⁠.

The outcomes of the numerical simulations performed on the model are in excellent agreement with the analytic ones for $ρ M L= 10 − 6, 10 − 5$⁠. Indeed, the error between the analytic and the computational results is of the order of $10 − 3$⁠. Upon exploring cases for which $ρ M L= 10 − 4$⁠, such an error rises to $10 − 2$⁠, thereby suggesting that lowering this discrepancy can be done by modeling the film itself through the CPE8H elements mentioned above.

Concerning the material properties, while the layer is described by a classical linear elastic law with Poisson ratio equal to $1/2$ (because of incompressibility), for the substrate a custom UMAT routine to properly simulate the OGH constitutive relation defined by (1) has been written. This is because the built-in ABAQUS routine neglects the contribution of the compressed fibers, deactivating them once they buckle.

Finally, an extended buckling analysis is performed by making use of an extended number of simulations (actually over 170), due to the need to cover the whole range of parameters. From Figs. 5–7 it is clear how the outcomes of FEM display a full agreement with theoretical ones. In particular, Figs. 5 and 6 show the critical strain and wavenumber with respect to the angle for $ρ M L= 10 − 4$ and different values of mismatch fibers/matrix $ρ F M$⁠. It is worth noting that, the closer to $1/3$ the dispersion is, the flatter the curves are, approaching to the neo-Hookean case when $κ=1/3$ (as well as when $ρ F M=0$⁠). Finally, in Fig. 7 the critical quantities are shown by setting a constant angle $θ= 70 °$ and varying the ratio $ρ M L$⁠, assuming a perfect alignment of the fibers. Furthermore, it is noteworthy that the representation of the critical strain and wavenumber with respect to the stiffness mismatch $ρ M L$ is susceptible to an interesting property. In fact, by using a bi-logarithmic scale as in Fig. 7, it becomes apparent that these quantities arrange along straight lines with slope $2/3$ and $1/3$⁠, respectively. For the sake of completeness, in Fig. 8 the post-buckling phase of a bilayer is shown. For this case, $ρ M L= 10 − 3$⁠, $ρ F M=2$⁠, $k 2=0.8393$⁠, $κ=0$⁠, and $θ= 90 °$ have been assumed. The amplitude has been evaluated as the absolute value of the deviation from the mean height of the surface of two representative points, namely, the global maximum and minimum.

Figure 8 shows that, depending on the contractile strain, the dimensionless amplitude defined by $A/ λ c r=A k c r/(2π)$ changes, and a re-organization of the surface emerges. Indeed, while the upper surface is initially flat, after the onset of the wrinkling the amplitude increases, with the current wavelength being $k c r/ λ 1 c r$⁠. The field (3) reproduces the kinematics that the bilayer has until the onset of the period-doubling. Note that in such a region, which extends up to ten times the wrinkling strain, the amplitude is governed by the well-known relation $A=h ε / ε c r − 1$⁠. This has been obtained in the absence of fibers, as reported by Chen and Hutchinson (2004), Huang (2005a), and Mane and Huang (2022) among others. Moreover, beyond a certain strain, two consecutive crests begin to join, and the valley between them flattens out, causing period-doubling. Finally, by further increasing the compression the waves move closer, until they make contact with one another (folding).

In order to explain the trend of the critical strain, an approximation has been constructed in Nguyen (2020) by treating the substrate as an orthotropic material. This leads to the derivation of stiffness moduli and Poisson’s ratios through an appropriate linearization. In this section, it is proved that the scaling laws (27) and (28) can actually be rewritten by taking into account the approximation above, thereby expressing them in terms of stiffness parameters of the linearized orthotropic substrate.

Therefore, in the range of parameters examined in this present work, the novel scaling laws (27) and (28) show minor deviations from the analytical model, and they are definitely much closer to the exact results compared to the other formulations available in the literature.

Compressed bilayers made of stiff thin films perfectly bonded to the top of fiber-reinforced deep soft substrates have been studied in this paper. To begin with, it has been shown how assuming a plate behavior for the thin top layer actually allows the reproduction of similar results to the ones obtained in Nguyen (2020), where a three-dimensional elastic behavior for the film itself was adopted. It has then been illustrated that for high stiffness mismatch ratios between the substrate and the adhering layer, both models essentially show the same results, although the dimensional-reduced model naturally entails significantly less computational cost. Upon varying the relative angle between the two families of fibers, the simulations performed for cases in which the substrate is much softer than the layer yield a symmetric and sinusoidal trend both for the critical strain and for the corresponding wavenumber.

Although complex theoretical and numerical analyses can be performed case-by-case, a prompt evaluation of the main wrinkling features for the bilayers at hand is often needed. This is certainly the case when it comes to comparing experimental findings with handy estimates of the topographic features of the corrugation in such systems. Unfortunately, no tools yet exist owing to the rapid evaluations mentioned above. To this end, appropriate scaling laws would certainly fill such a gap. In particular, no rigorously derived simple relations governing the critical strain and the wavenumber at the onset of wrinkling had been provided for highly mismatched hyperelastic fiber-reinforced bilayers before this present work. This drawback did include biological systems, even in the case of estimating the most basic features of the exhibited corrugated topography of organs such as arteries. Although in Nguyen (2020) several interesting analyses were performed (actually by a team of researchers involving three of the coauthors of this present paper) for flat fiber reinforced bilayers, a systematic rationale for mathematically deriving scaling laws for the features highlighted above was not pursued. Indeed, in such a paper only an ad hoc scaling for the strain at the onset of wrinkling was provided for specific cases.

The main result of this work relies upon the novel scaling laws for the critical strain and its corresponding wavenumber characterizing the initial wrinkling of compressed fiber-reinforced bilayers. This has been achieved only thanks to the outcomes of the simplified model developed in this paper. Indeed, close-form solutions to the bifurcation problem governing the balance of forces at the interface between the film (treated as a thin plate) and the fiber-reinforced substrate (modeled as an SRM) enabled one to analytically find explicit asymptotic expansions owing to the new scalings mentioned above. Remarkably, either the absence of fibers or their complete randomness reduce the new scaling laws to the well-known ones valid for neo-Hookean bilayers (see e.g., Sun , 2011 and Cao and Hutchinson, 2012).

In all the other cases, a modulating factor depending on the properties of the fibers turns out to govern the newly obtained laws (see Fig. 4). As expected, for certain fiber orientations and for certain fiber-matrix stiffness ratios, this can amplify the (dimensionless) wavenumber up to a factor $1.80$ and the corresponding wrinkling strain of over $3.2$⁠. This significant amplification would then be missing if the novel scaling laws were erroneously replaced either by existing estimates based on the film and the matrix properties alone or on the ad hoc relations mentioned above. The analyses in Secs. II–III, among other results, provide detailed comparisons between the outcomes of the novel scaling laws against both the results obtained from the dimensionally reduced analytical approach, and from the FE analyses. These results show how truly satisfactory the outcomes of the newly obtained scaling laws are relative to the corresponding analytic and numerical results. Furthermore, Sec. IV, Fig. 9 displays the discrepancies between the results coming from fully linear approaches (i.e., Vonach and Rammerstorfer, 2000) and the ad hoc relation found in Nguyen (2020) vs both the novel scalings and the analytic formulation for certain values of the parameters of the fiber-reinforced bilayers. Unlike the first two sets of results, Fig. 9 displays how the last two methods, which have been shown to essentially agree with the outcomes of the FE analyses in the previous sections, reliably hold throughout the whole range of variability of the fiber’s angle.

Finally, in Sec. V it has also been shown how material parameters obtained in Nguyen (2020), Appendix 2, providing orthotropic linearized moduli for the substrate can be used to express the newly obtained scaling laws.

Providing tools like appropriate and rigorously derived scaling laws is key for the analysis of geometrically complex situations involving compressed highly mismatched bilayers containing fibers, with whatever degree of dispersion relative to a main orientation they have. Indeed, having the availability of scaling laws for the main features of wrinkling, i.e., the critical strain at its onset and the associated wavenumber, becomes definitely useful when it comes to performing a firsthand comparison with experimental measurements even before running computational/analytical analyses.

The strategy undertaken in this paper to obtain such laws is currently under investigation for low-mismatch fiber-reinforced bilayers, that are even more amenable for soft biological tissues. Furthermore, new scaling laws in the presence of possible sources of inelasticity, such as growth, or even yielding and plasticity of the top layer in the presence of metallic films may be accounted for upon generalizing the methodology proposed in this paper to such situations. In other words, the systematic approach introduced in this paper to asymptotically expand complex representation formulas of the main wrinkling features paves the way for many other related problems, such as pre-stretch induced corrugation both in flat and curved fiber-reinforced bilayers, the latter being much closer to the geometry of the cross sections of wrinkled arteries.

A.M. and L.D. gratefully acknowledge the Italian Ministry of Universities and Research (MUR) in the framework of the project DICAM-EXC, Departments of Excellence 2023–2027 (Grant No. DM 230/2022).

L.D., A.C., and M.F. gratefully acknowledge partial support from the grants: Italian Ministry of Universities and Research (MUR) through Grant Nos. PRIN-20177TTP3S and PON “Stream” —ARS01 01182.

L.D. and A.C. gratefully acknowledge partial support from Grant No.PRIN-2022XLBLRX.

M.F. gratefully acknowledges partial support from Grant No. PRIN-2022ATZCJN.

L.D. also gratefully acknowledges partial support from the ERC through (1) FET Open “Boheme” Grant No. 863179, (2) LIFE GREEN VULCAN LIFE19 ENV/IT/000213, (3) ERC-ADG-2021-101052956-BEYOND, (4) the Italian Group of Theoretical Mechanics GNFM-INDAM of the National Institute of High Mathematics, and (5) the Italian Government through the 2023–2025 PNRR_CN_ICSC_Spoke 7_CUP E63C22000970007 grant, awarded to the University of Trento, Italy.

The authors wish to thank Dr. Thomas A. Witten, Homer J. Livingston Professor Emeritus within the Department of Physics and the James Franck Institute of the University of Chicago for the fruitful online conversations about the wrinkling of thin systems.

The authors have no conflicts to disclose.

A. Mirandola: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (lead); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). A. Cutolo: Conceptualization (supporting); Data curation (supporting); Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Software (equal); Validation (supporting); Writing – review & editing (equal). A. R. Carotenuto: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (supporting); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (supporting). N. Nguyen: Formal analysis (equal); Investigation (supporting); Methodology (equal); Supervision (supporting); Writing – review & editing (equal). L. Pocivavsek: Conceptualization (equal); Formal analysis (supporting); Investigation (supporting); Methodology (equal); Software (supporting); Validation (equal); Writing – original draft (supporting); Writing – review & editing (supporting). M. Fraldi: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Software (supporting); Supervision (equal); Validation (supporting); Writing – original draft (equal); Writing – review & editing (equal). L. Deseri: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (lead); Writing – original draft (equal); Writing – review & editing (lead).

The data that support the findings of this study are available within the article.