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

Wrinkling, creasing and folding are frequent phenomena encountered in biological and man-made bilayers made by thin films bonded to thicker and softer substrates often containing fibers. Paradigmatic examples of the latter are the skin, the brain, and arterial walls, for which wiggly cross-sections are detected. Although experimental investigations on corrugation of these and analog bilayers would greatly benefit from scaling laws for prompt comparison of the wrinkling features, neither are they available nor have systematic approaches yielding to such laws ever been provided before. This gap is filled in this paper, where a uniaxially compressed bilayer formed by a thin elastic film bonded on a hyperelastic fiber-reinforced substrate is considered. The force balance at the film-substrate interface is here analytically and numerically investigated for highly mismatched film-substrates. The onset of wrinkling is then characterized in terms of both the critical strain and its corresponding wavenumber. Inspired by the asymptotic laws available for neo-Hookean bilayers, the paper then provides a systematic way to achieve novel scaling laws for the wrinkling features for fiber-reinforced highly mismatched hyperelastic bilayers. Such novel scaling laws shed light on the key contributions defining the response of the bilayer, as it is characterized by a fiber-induced complex anisotropy. Results are compared with Finite Element Analyses and also with outcomes of both existing linear models and available adhoc scalings. Furthermore, the amplitude, the global maximum and minimum of ruga occurring under increasing compression spanning the wrinkling, period doubling and folding regimes are also obtained.


Introduction
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.Budday et al. (2017), Genzer and Groenewold (2006), Hohlfeld andMahadevan (2011), andHolland et al. (2020) among many others).For instance, as pointed out in Nguyen et al. (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 Figure 1), have been produced in the last two decades (see e.g.Alawiye et al. (2019Alawiye et al. ( , 2020)), Balbi et al. (2015), Ben Amar and Jia (2013), Y. Chen et al. (2021), Ciarletta and Ben Amar (2012), Ciarletta Figure 1: An original example of wrinkling in biological tissues is displayed in the images above.LEFT: A full-scale black-and-white SEM (HITACHI TM 4000 PLUS) image of the Ant's eye is reported.CENTER: A red-framed zoomed-in area from the left image is blown-up at the center of the figure: at the resolution reported in that frame, radial wrinkles are visible around each unit forming the eye's compound.Images are all original and taken by the coauthors of this paper affiliated at the LIMITS Laboratory, within the University of Napoli "Federico II".RIGHT: Nonlinear FE simulation.The top right of the figure displays the projection onto the X-Z plane of the resulting displacement field representing the zoom of the blue-framed inset taken from the center of this figure: details simulated at this level of observation reveal a few azimuthal wrinkled crowns separating the central undisturbed zone from the radial wrinkles observed before at a coarser resolution.The bottom right of the figure displays a 3D image of the vertical displacement resulting from the FE analysis, thereby reliably reproducing the experimental observation reported above.et al. (2014), Genzer and Groenewold (2006), Goriely (2017), Kai (2022), Mostafavi Yazdi and Baqersad (2022), Mottahedi and Han (2016), and Nath et al. (2020) 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), Cutolo et al. (2020), and Pocivavsek et al. (2008) 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) (Sect.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 et al. (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 et al. (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 et al. (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.Hayn et al. (2020) and Taber and Humphrey (2001), 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.Braeu et al. (2017), Cyron et al. (2016), and Robertson and Watton (2013)).This leads to overall characterizations of the constitutive behavior of such complex systems (see e.g.Bellini et al. (2013), Gasser et al. (2005), and Humphrey and Rajagopal (2002) and ref.s 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 et al. (2000) (called OGH in the sequel).Such a constitutive equation has been utilized in Nguyen et al. (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 et al. (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 towards 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 Sect. 2 the approach undertaken in Nguyen et al. (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 et al. (1994)) is assumed, and the simplified constitutive SRM law for the fiber-reinforced substrate is considered.In Sect.3, 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 wave number 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 Cao and Hutchinson (2012) and Sun et al. (2011).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 Cao and Hutchinson (2012) and Sun et al. (2011)) 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 Methods (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 Sect.3. Finally, in Sect. 4 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 et al. (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.

A simplified model
In Nguyen et al. (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: (a) bilayers for which the mismatch between the elastic moduli of the layer and of the substrate is very high (i.e. between 10 4 ÷ 10 6 ); (b) the layer being considered as very thin compared to the substrate (which, in mathematical terms, is in fact assumed infinitely deep).
Item (b) allows for considering thin plate behavior for the top layer, and this inspired many studies on fully isotropic, homogeneous and elastic bilayers already present in the literature (see e.g.Cao and Hutchinson (2012) and Sun et al. (2011) and ref.s cited therein).In the present analysis, a thin plate theory to model the thin film bonded to the fiber-reinforced substrate is adopted.As previously mentioned, the latter here is modeled through the OGH constitutive equation.Such a material has a strain energy that is additively composed of two terms.The first one is due to the classical neo-Hookean matrix.The second term is due to the presence of fibers, organized in families, dispersed in the matrix, and reciprocally oriented with one another at a certain angle 2θ (Figure 2).Finally, the total strain energy density of the substrate W s is given by the sum of those two contributions, i.e.
W s,OGH = W s,matrix + W s,fibers , with W s,matrix = µ M (I 1 − 3), I 1 := tr (C) , where N is the number of families of fibers in the matrix and F is the deformation gradient.The term 2µ M stands for the shear stiffness of the matrix, k 1 is a parameter related to the stiffness of the fibers and Figure 2: Schematics of the plane strain bilayer system.The substrate is composed by two families of fibers with relative angle 2θ embedded in a neo-Hookean matrix.In (a) the bilayer is undeformed, with thickness h and length L 0 .This configuration is assumed to be the reference one, with the material coordinates system X 1 − X 2 , with a substrate much deeper than the layer (h/H → 0).In (b) the deformed configuration is shown.There u L is an imposed contractile displacement, ϵ L = u L /L 0 is the corresponding strain, λ 1 = 1 + ϵ L is the resulting stretch and, hence, the bilayer's deformed length is λ 1 L 0 .This geometry remains valid for higher values of the stretch λ 1 ≥ λ cr 1 , where λ cr 1 is the longitudinal stretch at the onset of wrinkling (see Eqn. ( 13)), spanning the whole wrinkling regime until period-doubling starts (see Fig.

in Sect. 3). During deformation the angle between the two families of fibers takes the form
, where l m , m = 1, 2 are defined below (2).The upper left corner displays a magnification of the reactive tractions arising in the bilayer due to the imposed displacement (the tractions between the layer and the substrate are not shown to scale relative to one another).
k 2 is a non-dimensional parameter determined experimentally.The term I 1 is the first invariant, i.e. the trace of any of the two Cauchy-Green tensors.The argument of the exponential defined above, besides k 2 is given by where l m = (cos θ, sin θ, 0) T is the unit vector representing the m-th fibers family with respect to the horizontal axis.It is worth noting that I 4m is the magnitude (squared) of the extension/contraction of the fibers.
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 et al. (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.
Here the configurational changes of such an interface are analyzed through a small-on-large approach, consistent with the existing literature on compressed bilayers formed by stiff films on non fiber-reinforced soft substrates (see e.g.Cao and Hutchinson (2012), Hutchinson (2013), Jiang et al. (2007), and Wang et al. (2023) and references cited therein, among many others).To this aim, following Nguyen et al. (2020), Eqn.(4) (see e.g. also Shield et al. (1994), Eqn.
(3) for the sole displacement field, and Sun et al. ( 2011), Eqn.s (2.1) ÷(2.3)) a sinusoidal perturbation (of amplitude δ ≪ 1) is given to a homogeneous planestrain, volume-preserving finite deformation of the substrate induced through a longitudinal shortening (imposed in the direction 1, as specified below), i.e. x (3) It is worth noting that λ 2 = 1/λ 1 , p is the hydrostatic pressure needed to maintain incompressibility (namely the reactive action needed to keep isochoricity of the substrate), whereas p 0 , p 1 and α are constants (to be determined through boundary conditions).Of course, in (3) the pairs (x 1 , x 2 ) and (X 1 , X 2 ) give the coordinates of a generic point in the deformed and in the reference configuration, respectively.As the wrinkling of the interface occurs, the corrugation will be characterized by a space wavelength λ, or by its corresponding wavenumber k, linked by the relation k = 2π/λ.The deformation and pressure fields (3) satisfy some basic considerations about the nature of the problem.Along the direction X 2 , the perturbation must fade at a great distance from the interface, and horizontally the motion must be periodic.Moreover, as shown further, the fields are solutions of the equilibrium of the substrate, for a suitable choice of α and p 1 .The incompressibility condition det F = 1 is satisfied at the first order, i.e.Pence and Song (1991) and Yue et al. (1994)).This last equation is identically satisfied by (3) if Finally, as shown in Figure 6(e) the bifurcation is characterized by a sinusoidal profile and, by observing Figure 8, this extends up to ten times the strain at the onset of the wrinkling.Tractions acting on the film coming from the substrate must be evaluated in order to characterize which superimposed deformations are admissible for the bilayer.To do this, the first Piola-Kirchhoff stress tensor for the substrate P s can be computed as follows: When the wrinkling has not yet occurred, the constant p 0 can be obtained noting that the normal traction at the interface P s 22 (δ = 0) vanishes.By letting δ = 0 the following expression for p 0 is determined, leading to which, introduced into (3) and then into (4), allows us to write the equilibrium equations for the substrate where ( • ) ,j indicates ∂( • )/∂X j and the repeated index means summation.It is worth noting that a full analytic proof of the fact that (3) is a representation formula for the solution of the boundary value problem at hand for purely neo-Hookean bilayers (with no fibers) formed by stiff films on softer substrates could be provided by generalizing the approaches utilized in Pence and Song (1991), Qiu and Pence (1997), and Wang et al. (2023) to account for the presence of the fibers reinforcing the substrate.Of course, unlike the case of neo-Hookean materials, for fiber-reinforced substrates the constants α and k characterizing the eigenmodes are expected to be influenced both by the fiber and by the matrix parameters.
In particular, solving equation ( 6) for the systems under consideration yields four different pairs of solutions (α, p 1 (α)).Note that the solution of α is formed by two complex conjugate pairs, which differ from one another with the sign of their positive part.Nevertheless, only two of those pairs (α, p 1 ) can be used, more precisely the ones that have an α with a strictly non-negative part, as the perturbation effects must vanish at long distances from the interface (it is worth noting that going inward deep into the substrate entails negative values for X 2 , see fig. (2)).After labeling α 1 and α 2 the values satisfying equations ( 6), the resulting quantities in the reference configuration are given by a linear combination of the respective eigenfunctions, i.e.
In the case that the stiffness of the fibers approaches zero, it is worth noting that ( 5) and ( 6) give the same results found in Sun et al. (2011) (note that Q used in Sun et al. (2011) is equal to 2µ M ), namely: It is worthy of mention that the approach followed by Eqn.(3) to Eqn. ( 8) is analogous to the one introduced in Nguyen et al. (2020).The assumption (b) introduced above, i.e. the layer is assumed to be very thin compared to the substrate, is now useful.In this case it appears reasonable to assume that the layer starts to wrinkle with a wavelength that is large compared to the thickness of the upper layer.Henceforth, a plate behavior with a single bending axis which lies on a semi-infinite space can be assumed (see e.g.(Shield et al., 1994)).Upon utilizing the balance equation at the interface between the top layer and the substrate (in the reference configuration) the following expressions are obtained (see Shield et al. (1994) and Sun et al. ( 2011)) where is the longitudinal strain in the absence of prestretch, P L is the corresponding longitudinal stress arising across the layer (Shield et al., 1994), and E L and ν L are the Young modulus and Poisson ratio of the layer, respectively.Furthermore, P s ij and u i = x i − λ i X i , (i, j = 1, 2), are the stresses exchanged between the substrate and the layer and the displacements at the interface (hence evaluated at X 2 = 0 and obtained from Eqn. ( 7)), respectively.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, ρ ML means stiffness of the matrix (M) forming the substrate over the one of the layer (L).Recalling that λ = 2π/k denotes the spatial wavelength of periodic wrinkles, and by introducing k h = 2πh/λ, namely its corresponding non-dimensional wavenumber, by setting ρ FM = k 1 /µ M the ratio between stiffness information about both the fibers and the matrix, and by noting that ρ ML = 6µ M /E L is the stiffness ratio between the substrate and the layer, the substitution of expressions ( 7) into (9) leads to the following homogeneous linear system in the amplitudes C 1 and C 2 appearing in (7): where M is the resulting coefficients matrix and the pair C 1 , C 2 characterize the wrinkling eigenmodes.
For the sake of brevity, the explicit form of M is omitted, although available upon request.Of course, the amplitude modes C 1 and C 2 are associated with the values of ε for which a bifurcation of equilibrium occurs, i.e. such that where It is worth noting that ε L appears only in the first row of M, hence its determinant is linear, as well.
Following the findings of Pence and Song (1991) and Yue et al. (1994) (see e.g.Fig. ( 4) in both papers), and employed later by Sun et al. (2011), among the possible values satisfying Eqn. ( 12), only the ones corresponding to the smallest wavenumber are of interest.This leads to writing the following optimality conditions, namely the ones governing both the minimum strain at which the onset of wrinkling occurs and the corresponding wavenumber: Indeed, in full analogy with Sun et al. (2011), the conditions above can be shown to deliver the critical stretch and the corresponding non-dimensional wavenumber at which wrinkling occurs.
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.
To this end, the OGH constitutive equation here is replaced by the Standard Reinforcing Model, SRM, mentioned above (for details see e.g.(Kurashige, 1981;Qiu & Pence, 1997;Triantafyllidis & Abeyaratne, 1983)).The SRM energy density reads as follows: where γ has the dimension of an elastic modulus.Of course, (15) must be added to the neo-Hookean term, accounting for the hyperelasticity of the matrix.It should be noted that, as the parameter k 2 approaches zero, the derivative of W fiber,OGH with respect to the strain invariant E m coincides with the one of W fiber,SRM when γ = k 1 .Solutions of ( 14) obtained by utilizing the SRM strain density energy are shown in Figure 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 et al. (2020) illustrated that the OGH law does not depend on k 2 in the small strain regime: in the sequel, (see Figures 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 et al. (2020), is shown in Figure 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 Figure 2) and for different stiffness ratio ρ FM .This has been done by setting, as in Nguyen et al. (2020), k 2 = 0.8393 and κ = 0.Moreover, ν L = 0.5 has been imposed for the incompressibility of the layer and a stiffness ratio ρ ML = 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 ρ ML , for the considered values of ρ FM 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.

A new asymptotic law for fiber-reinforced bilayers
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: (a) the upper layer is considerably stiffer than the substrate, namely ρ ML = 6µ M /E L ≤ 10 −4 ; (b) the fibers ratio ρ FM = 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; (c) the critical strain is very small (ε cr ≤ 10 −3 , see e.g.Sun et al. (2011) in the absence of fibers) and therefore, as Because no prestretch is considered, the eigensolutions given by ( 7) are linearized around λ cr 1 = 1.In this way, the resulting quantities, namely the critical strain and non-dimensional wavenumber, will have no dependence on the stretch; (d) the critical strain and non-dimensional wavenumber are approximately constant when k 2 changes (at least for ρ ML = 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; (e) the Poisson ratio of the layer is ν L = 1 /2.
By analogy with Allen (1969) and Sun et al. (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.
In order to do so, one can start by taking advantage of the assumption (a), i.e. ρ ML ≈ 0. Henceforth, a Taylor expansion of (∂ ε cr )/(∂ k h,cr ) (the second equation appearing in ( 14)) with respect to ρ ML around zero can be considered.Upon equating the obtained expression to zero, it is not difficult to check that a closed-form solution of an algebraic third-order equation in k h,cr can be obtained.The only possible physically admissible root of such an equation reads as follows: where It is worth noting that the zero-order term of the expansion ( 16) vanishes.From the physical viewpoint, this can be interpreted as the substrate becoming extremely soft relative to the thin top layer, when the matrix-to-film stiffness ratio ρ ML → 0. Hence, the buckling strain of such layer (for a finite depth of one unit length) of thickness h turns out to be ϵ cr ∼ (h/L 0 ) 2 , as the wavelength tends to the physical length of the film.The thinness of the latter implies h/L 0 → 0, hence both the critical strain and the corresponding dimensionless wave number, k h,cr ρ ML →0 ∼ h/L 0 , tend to zero.It is worth noting that g 1 depends upon ρ FM = k 1 /µ M , relating the stiffness of the fibers, weighted against their volume concentration, and the shear modulus of the substrate.For low densities of fiber reinforcements ρ FM tend to zero.Therefore, g 1 can be replaced by a suitable expansion obtained as follows: , where where h 1j (θ, κ), j = 0, 1, 2 are the terms of the expansions.Note that such expressions are valid only if g i are positive functions for every value of ρ FM and for 0 ≤ θ ≤ π/2.This is reasonable since the wavenumber is a positive quantity.Finally, carrying out the computations of the previous expressions, one has Henceforth, the resulting critical (non-dimensional) wavenumber takes the following form: and, as was previously pointed out, this value is unique.
A corresponding asymptotic expansion for the wrinkling strain can also be obtained.Indeed, by substituting ( 19) in the first equation of ( 14), and by computing the Taylor expansion of the resulting expressions up to second order, the following form for ε cr is achieved: . (20) Similarly to the previous case, the zero-order term in the Taylor expansion for the argument in ( 20) is zero and, furthermore, in this specific case, even the first-order one identically vanishes.In particular, the zeroorder term corresponds to ε cr ρ ML =0 , which is again consistent with having a compressed free-standing (because ρ ML = 0 would essentially mean to have a substrate with zero stiffness relative to the top layer) infinitely thin film with a finite length.Furthermore, in order to achieve an irreducible representation for the critical strain, a Taylor expansion of t 2 can be provided.By expanding that with respect to ρ FM around 0, the following expression follows: where Furthermore, by carrying out the computations for q 2i , i = 0, 1, 2, their values read as follows: Upon substituting ( 22) into ( 21), and the obtained result into (20), the following asymptotic expression for the critical strain is delivered: For the sake of consistency, ( 19) and ( 23) are explored for the simpler case of neo-Hookean bilayers.As expected, those expressions reduce to what already found in Sun et al. (2011) and Cao and Hutchinson (2012) by assuming ν L = 1/2, i.e.
k h,nh = 3 3ρ ML (24) It is evident by inspections of ( 19) and ( 23) that the presence of the fibers turns out to significantly influence both the wrinkling strain and the corresponding wavenumber.Indeed, the asymptotic expressions above involve the following modulating factor: Henceforth, the quantities mentioned above can be written as where ζ(ρ FM , κ, θ) (see Figure 4) is defined by the expressions above and its square has the meaning of amplitude factors for the neo-Hookean values of k h and ε cr , respectively.It is just as simple to verify that the modulating factor is really an amplification function.This reduces to 1 in both cases in the absence of fibers, as shown in Figure 4, and when κ approaches to 1/3.This demonstrates that in the case of a total dispersion of fibers within the matrix these give no contribution to the critical dimensionless wavenumber and strain.
Using expressions ( 27) and ( 28), and the results achieved in ( 14), the plots displayed in Figures 5 ÷ 7 for different angles, stiffness ratio ρ ML , ρ FM and dispersion factors are obtained.It is worth noting that Eqn.s ( 27) and ( 28) work particularly well for small ρ ML ratios: this can be achieved even for stiff fibers, provided that concentration per unit volume is adequately low.Another aspect is related to the behavior of the curves: for small ρ ML , the trends of both critical strain and wavenumber are symmetrical with respect to π/4.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)  Based on the geometry described above, it is clear that the smaller the ratio ρ ML 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 ρ ML 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 et al., 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 − ν 2 L ) −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 λ cr /40.The outcomes of the numerical simulations performed on the model are in excellent agreement with the analytic ones for ρ ML = 10 −6 , 10 −5 .Indeed, the error between the analytic and the computational results are of the order of 10 −3 .Upon exploring cases for which ρ ML = 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 Figures 5 ÷ 7 it is clear how the outcomes of FEM display a full agreement with theoretical ones.In particular, Figures 5 and 6 show the critical strain and wavenumber with respect to the angle for ρ ML = 10 −4 and different values of mismatch fibers/matrix ρ FM .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 ρ FM = 0).Finally, in Figure 7 the critical quantities are shown by setting a constant angle θ = 70 • and varying the ratio ρ ML , 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 ρ ML is susceptible to an interesting property.In fact, by using a bi-logarithmic scale as in Figure 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 Figure 8 the post-buckling phase of a bilayer is shown.For this case ρ ML = 10 −3 , ρ FM = 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/λ cr = Ak cr /(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 cr /λ cr 1 .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 √ ε/ε cr − 1.This has been obtained in the absence of fibers, as reported by X. Chen and Hutchinson (2004), Huang et al. (2005), 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).

Approximation based on linearized orthotropic properties
In order to explain the trend of the critical strain an approximation has been constructed in Nguyen et al. (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.In a fiber-reinforced material, the principal directions P 1 − P 2 are identified by the orientation of fibers and by their normal, which may not coincide with the "natural" reference system X 1 − X 2 used to solve the equilibrium problem.Using the expressions found by Nguyen et al. (2020), the stiffness moduli with respect to the natural directions result as follows Henceforth, their value in the principal system is obtained by placing θ = 0, by assuming the axis P 1 and    As the contractile strain increases different patterns emerge on the surface.Initially, the film is flat but, after reaching the critical strain ε cr wrinkling occurs.In such a region the amplitude turns out to scale like h √ ε/ε cr − 1.In the post-wrinkling regime period-doubling and creases can be observed.
X 1 aligned, so that Noting that their ratio is it is clear it can be substituted into the amplitude function ( 27), obtaining In this way the scaling laws ( 27) and ( 28) can be rewritten as follows from which one can see how the response depends on the ratio of the orthotropic stiffness moduli in the principal system and on the angle that fibers have with respect to the natural one.By assuming a bilayer with a geometry similar to what is considered in this paper, though with a linear elastic orthotropic substrate instead of a fiber-reinforced one, Vonach and Rammerstorfer (2000) obtained scaling laws for such systems.In Eqns.( 22) and ( 23), such authors provided explicit formulas for the semi-wavelength and the critical longitudinal load at the onset of wrinkling of an isotropic thin plate resting on an elastic foundation.By adapting these results to write down the actual wavelength and the critical strain it follows that where is the bending stiffness of the top layer, while k s is the substrate stiffness.It is noteworthy that these relations are particularly similar to the ones valid for isotropic bilayers.Depending on the specific problem under consideration, k s can assume different forms.For instance, the outcomes of Eqn.s (34) and (35) displayed in Figure 9, arise by choosing the substrate's stiffness k s as in Eqn.(21) of Vonach and Rammerstorfer (2000).It must be mentioned that for the evaluation of such a quantity the linearized orthotropic moduli obtained in Appendix 2 of Nguyen et al. (2020) have been used.Furthermore, in Eqn.(7) of that same paper, the authors constructed a scaling law for the critical strain that fits the results obtained for low substrate/layer mismatches with good agreement.Such an expression is neither based on an asymptotic expansion nor on a rigorous derivation, and it makes use of an ad-hoc elastic module, E s eff , such that the following relations hold: Finally, unlike Eqn. ( 22) and ( 23) by Vonach and Rammerstorfer (2000) and the Eqn.( 7) by Nguyen et al. (2020), it is worthy of mention that Figure 9 shows how scaling laws ( 27) and (28) derived in this paper well capture the quasi-symmetric trend of the critical wavenumber and critical strain at the onset of wrinkling.Indeed, the case ρ ML = 10 −3 portrayed in Figure 9 presents a slight asymmetry when ρ FM = 10 that is not present whenever ρ ML ≤ 10 −4 and 0 ≤ ρ FM ≤ 10.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.

Conclusions
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 et al. (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 wave number 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 et al. (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, closeform 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.(Cao & Hutchinson, 2012;Sun et al., 2011)).In all the other cases, a modulating factor depending on the properties of the fibers turns out to govern the newly obtained laws (see Figure 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 Sect.s 2 ÷ 3, 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, Sect.4, Figure 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 et al. (2020) versus 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, Figure 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 Sect. 5 it has also been shown how material parameters obtained in Nguyen et al. (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 lowmismatch 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.

Figure 3 :
Figure3: Comparison between the critical strain and the non-dimensional wavenumber between the plate model (dots) and the 3D-solid one (solid line) with respect to the angle θ and for three different values of ρ FM .The considered ratio ρ ML = 10 −5 is in the middle of the range of interest.The diamonds are the critical strains obtained from the Standard Reinforcing Model (15).

Figure 4 :
Figure 4: The amplitude function ζ(ρ FM , κ, θ) varying ρ FM and θ, assuming κ = 0 (perfect alignment).In particular, two projections are shown: the first one, within the plane ρ FM − ζ shows the amplification for a fixed angle, namely θ = 70 • .The second one, in the plane θ − ζ represents the sinusoidal amplification by setting ρ FM = 8.Note that when the dispersion factor κ approaches zero the function is an horizontal plane with ζ = 1.
and Yue et al. (1994) (see the barreling curves reported in such papers in Fig.s 4, same numbering in each paper).

Figure 6 :
Figure 6: Comparison of the critical dimensionless wavenumber between FEM, the plate model (dots) and the asymptotic expansion (19) (solid line) with respect to the angle θ, ρ ML = 10 −4 and ρ FM = 2, 5, 10 ((a), (c), (d)).In (b) the wrinkling mode, with normalized amplitude, of a representative set of values is shown.Finally, in (e) the magnitude of the normalized displacements resulting from the FE analysis is plotted.

Figure 7 :
Figure 7: Comparison of the critical strain and wavenumber for κ = 0, θ = 70 • and k 2 = 0.8393, as a function of the stiffness mismatch between the matrix and the film ρ ML .Note that the plots are on a bi-logarithmic scale.

Figure 8 :
Figure8: Dimensionless amplitude A/λ cr for the case ρ ML = 10 −3 , ρ FM = 2, k 2 = 0.8393, κ = 0 and θ = 90 • .As the contractile strain increases different patterns emerge on the surface.Initially, the film is flat but, after reaching the critical strain ε cr wrinkling occurs.In such a region the amplitude turns out to scale like h √ ε/ε cr − 1.In the post-wrinkling regime period-doubling and creases can be observed.