A displacement-based approach to geometric instabilities of a film on a substrate

Abstract

When a thin film adhered to a compliant substrate is growing, it will eventually buckle in order to release the compressive stresses accumulated within the film due to growth. Such geometric instabilities caused by compressive stresses prevail among all living systems in nature and their outcomes range from highly beneficial to destructive. Therefore, understanding compression induced instabilities is of crucial importance. Note that the origin of the “compression” need not necessarily be differential growth, as it may be due to pre-stretch or thermal expansion. A commonly accepted solution strategy for instabilities in bilayer structures dates back to the seminal work of Allen and employs the Airy stress functions. Owing to its reliance on a stress-based approach, the Allen solution is limited to linear two-dimensional problems and its success depends entirely on choosing an appropriate Airy function. The main objective of this contribution is to circumvent these limitations via a displacement-based approach formally suitable for three-dimensional problems, anisotropic materials, and even applicable to finite deformations. Furthermore, the Allen solution in its original form is valid for the plane-stress condition but often it is mistakenly compared with the numerical simulations corresponding to the plane-strain condition. We analyze the subtle difference between the solutions associated with the plane-strain and plane-stress conditions. Next, the analytical solution is compared against the computational results using the finite element method via eigenvalue analysis. Finally, it is briefly explained how the current approach can be utilized beyond the classical bilayer systems.

Keywords

Wrinkling growth-induced instability bilayer structure

1. Introduction

The instabilities that occur in bilayers due to compressive stresses have been the subject of special interest in the past few decades due to their pervasiveness in nature and the wide array of applications of the topic in material design. Of particular interest to this contribution is wrinkling in bilayer systems consisting of a growing thin film adhered to a deep compliant substrate. Following the pioneering works of Allen [1] and Biot [2], such structures have been investigated extensively to understand the behavior of living tissues before and beyond the onset of wrinkling. The literature on the instabilities due to compressive stresses can be categorized into two classes: bilayer wrinkling and extended bilayer wrinkling. Furthermore, it is also necessary to briefly go over some recent works on growth modeling and mechanics of growth. The remainder of this section provides a brief account of the literature on these topics.

Bilayer wrinkling occurs when a thin stiff film attached to the surface of a compliant substrate undergoes a critical amount of compressive stress and, thus, forms sinusoidal patterns in order to release energy. What is deemed as the critical stress and strain is the minimum amount of stress and strain at the onset of wrinkling. Following the thorough expositions of Allen [1] and Biot [2] on the mechanics of bilayer wrinkling, many aspects of this concept have been investigated in depth. In particular, the wrinkling of deposited thin films on compliant substrates has been studied extensively and utilized in various applications [3 –7] using both platinum and gold films on polydimethylsiloxane (PDMS) or other elastomers. In addition, several contributions [8 –12] have detailed on the geometry of bilayer structures or loading conditions as well as stiffness ratios. Recently, Holland et al. [13] have considered a large range of stiffness ratios in bilayers with various loading conditions. They have revealed that the different loading conditions are only distinguishable in the low stiffness regime, and that these differences disappear when measures of effective strain, stiffness, and wavelength are used, see also [14 –21] for more complicated scenarios as well as large deformations. Post-buckling analysis show further types of instabilities apart from the initial wrinkling of the bilayer structure. Examining the range beyond the critical stress for the onset of wrinkling, the morphology and amplitude of geometric instabilities have been studied in [22 –33] among others. A detailed review on wrinkling, creasing, and folding that occur in soft materials with various geometries and under loading can be found in the work of Li et al. [34].

Extended bilayer wrinkling is an extension of bilayer wrinkling where the stiff film is embedded in a medium and has been actively investigated very recently. For instance, Xie et al. [35] have derived an analytical solution for the critical compressive strain and critical wavelength for the wrinkling of a film embedded between two different soft layers, also conducting post-buckling analysis of the structure. The variation of these critical values with changing stiffness ratio has been investigated in [36]. Brau et al. [37, 38] have considered the development of wrinkles and further instabilities with an elastomer or liquid substrate, looking into the symmetry breaking of the system and its effect on post-buckling behavior. In another example, Li et al. [39] have examined the instability patterns formed when periodic fibers are embedded in a soft matrix, showing the softening of the structure after buckling by constructing stress–strain curves, also relating fiber diameter to critical wavelength and amplitude of formed waves. Colin and Holland [40] have studied critical strains in the case where the stiffness of the matrix above and below the embedded layer have different shear moduli, also identifying different modes of wrinkling in a homogeneous matrix. For further extensions of bilayer wrinkling and related studies, see [41 –46] among others. The instabilities of an immersed film has also been studied in [47 –49] with different geometries and from both theoretical and experimental viewpoints.

Growth modeling and mechanics of growth explains the driving force behind the formation of various biological tissues as a result of constrained or differential growth, see the review by Kuhl [50]. A commonly accepted strategy to model the growth in the context of continuum mechanics is to multiplicatively decompose the deformation into its growth and elastic parts as proposed by Rodriguez et al. [51]. The mechanism of growth and its modeling has been a compelling area of study with a focus on the kinematics, continuum mechanics treatment [51 –57], and thermodynamics [58], or explaining growth using mixture theory [59 –61]. In addition, the growth of biological structures such as horns, tusks, rods, or cells have been modeled and studied in [62 –64] while growth models at sub-cellular scales, plant growth, and bone remodeling have been addressed in [65] among others. Hosseini et al. [66, 67] have elaborated on how mechanical forces shape the developing eye, through experiments and computational 3D models, focusing on differential growth in layers of the eye, see also [68, 69]. More specifically, growth-induced instabilities in the bilayers have been used in explaining the formation of skin wrinkles [70], investigating the instabilities formed in artificial elastomeric skins [71], understanding the development and mechanics of brain folding [72 –75], and in exploring the development of patterns in growing tubular tissues [76 –79]. The buckling behavior of cylindrical geometries due to differential growth has been studied extensively by Moulton and Goriely [80] among others. Eskandari et al. [81] have studied mucosal folding in the post-buckling range, with the airway wall having temporally evolving material properties leading to elastosis. Skin growth [82] has also been investigated in several contributions, for example finding applications in tissue expansion [83] as a novel technique to allow for the growth of extra skin for reconstruction. Apart from these cases found in nature, applications of bilayer wrinkling in other fields include optical sensors [84], pressure sensors [85], microfluidic devices [85, 86], stretchable electronics [87 -90], dielectric plates [91, 92], and material behavior [93]. From a mathematical perspective, Yavari [55] have eliminated the use of an intermediary configuration in the approach of Rodriguez et al. [51] by constructing a geometry theory of growing solids, where the growing body remains stress-free in the material manifold, see also [94 -96]. Taber [97] provides an extensive review of some of the works done on biomechanics of growth, along with remodeling and morphogenesis following the multiplicative decomposition approach.

2. Equivalent stiffness of an infinite half-space

The objective of this section is to compute the resistance (equivalent stiffness) of an infinite half-space (substrate) against prescribing a sinusoidal deformation on its surface. A commonly accepted strategy to compute the substrate equivalent stiffness dates back to the seminal work of Allen [1] and employs the Airy stress functions to solve the problem briefly explained in Section 2.2. While the Allen solution is correct, owing to its reliance on Airy functions, it is limited to two-dimensional linear isotropic problems and its success depends entirely on choosing an appropriate Airy function. A key feature of this contribution is to propose a displacement-based approach, elaborated on in Section 2.3, to compute the equivalent stiffness of the substrate without recourse to Airy stress functions and, hence, suitable for three-dimensional as well as anisotropic problems. The current approach can even be formally applied to finite deformations, unlike Allen’s solution.

2.1. Governing equations

Let $u$ denote the displacement field in the substrate and $ε$ be the associated strain field. The Cauchy stress $σ$ is related to the strain $ε$ via a linear mapping $c$ where $c$ is referred to as the constitutive tensor. General field equations of elasticity are

\begin{matrix} ε = \nabla^{sym} u, \nabla \times ε \times \nabla = 0, div σ = 0, σ = c : ε . \end{matrix}

(1)

in which (1)₂ ensures the compatibility of the strain field and (1)₃ is the linear momentum balance in the absence of body forces. The angular momentum balance requires the symmetry of the stress field $σ$ , which is satisfied a priori via the symmetry properties of the constitutive tensor $c$ and, thus, is omitted from the equations. The fourth-order constitutive tensor $c$ for the case of isotropic linear elasticity¹ of interest here reads

c = 2 μ_{s} i^{sym} + 2 λ i^{vol} with i^{sym} = \frac{1}{2} [i \bar{\otimes} i + i \underline{\otimes} i] and i^{vol} = \frac{1}{2} i \otimes i,

(2)

where $μ_{s}$ and $λ_{s}$ are the Lamé parameters and the non-standard tensor products are defined as

[i \bar{\otimes} i]_{abcd} = [i]_{ac} [i]_{bd}, [i \underline{\otimes} i]_{abcd} = [i]_{ad} [i]_{bc}, i : second - order identity tensor .

Note, the strain $ε$ can be computed from the stress $σ$ via the inverse of the constitutive tensor as

ε = c^{- 1} : σ with c^{- 1} = \frac{1}{2 μ_{s}} i^{sym} - \frac{λ_{s}}{2 μ_{s} [μ_{s} + λ_{s}]} i^{vol} .

(3)

While it may not be possible to obtain a closed-form solution for an arbitrary boundary value problem, there exist analytical solutions for various simplified scenarios. The general field equations of elasticity (1) can be reformulated depending on the nature of the problem of interest. Such reformulations fall into two categories of Stress formulation and Displacement formulation briefly addressed below.

Stress formulation also referred to as the Beltrami–Michell equation is obtained through imposing the strain compatibility (1)₂ in terms of stresses using the constitutive relation (1)₄ and then imposing the linear momentum balance (1)₃ as

\nabla \times ε \times \nabla = 0 with ε = c^{- 1} : σ \Rightarrow \nabla \times c^{- 1} : σ \times \nabla = 0,

\nabla \times c^{- 1} : σ \times \nabla = 0 with div σ = 0 \Rightarrow Beltrami - Michell equation,

which for the two-dimensional plane-strain problem of interest here reads

\begin{matrix} Δ (σ : i) = 0 with Δ {•} = div (\nabla {•}) = \frac{\partial^{2} {•}}{\partial ξ^{2}} + \frac{\partial^{2} {•}}{\partial η^{2}}, \end{matrix}

(4)

where $Δ {•}$ denotes the Laplacian operator. Note, the stress formulation is more appropriate to use with traction-type boundary condition problems. The governing equation for the stress formulation even in two-dimensional problems is still rather complex and analytical solutions are commonly determined through an Airy stress function $ϕ$ satisfying the biharmonic equation

Δ Δ ϕ = 0 with Δ Δ {•} = Δ (Δ {•}) = \frac{\partial^{4} {•}}{\partial ξ^{4}} + 2 \frac{\partial^{4} {•}}{\partial ξ^{2} \partial η^{2}} + \frac{\partial^{4} {•}}{\partial η^{4}},

(5)

whereby $Δ Δ {•}$ denotes the biharmonic operator. For the Airy stress function $ϕ$ , the stresses

σ_{ξ ξ} = \frac{\partial^{2} ϕ}{\partial η^{2}}, σ_{ξ η} = - 2 \frac{\partial^{2} ϕ}{\partial x \partial y}, σ_{η η} = \frac{\partial^{2} ϕ}{\partial ξ^{2}},

(6)

satisfy the Beltrami–Michell equation (4) a priori. Once the stresses are calculated, the strains can be obtained via the constitutive relation (1)₄and, eventually, integrating the strains renders the displacement field in the domain.

Displacement formulation also referred to as Navier–Lamé equation is obtained via replacing the stress from the constitutive relation (1)₄ in the linear momentum balance (1)₃ and then imposing the strain definition (1)₁ as

div σ = 0 with σ = c : ε \Rightarrow div (c : ε) = 0,

div (c : ε) = 0 with ε = \nabla^{sym} u \Rightarrow Navier - Lam é equation .

Inserting the constitutive tensor (2) into the general format of the Navier–Lamé equation and using the identities

ε : i = div u, div (\nabla^{t} u) = \nabla (div u), div (ε : i i) = \nabla (div u),

renders a simplified format of the Navier–Lamé equation as

[λ_{s} + μ_{s}] \nabla (div u) + μ_{s} Δ u = 0 .

(7)

The derivations in Sections 2.2 and 2.3 correspond to the stress formulation and displacement formulation, respectively. More precisely, the well-established Allen [1] solution elaborated on in Section 2.2 is based on the stress formulation while the current proposition in Section 2.3 is based on the displacement formulation. We emphasize that the surface instability of an infinite half-space is intrinsically suited for the displacement formulation.

2.2. Allen stress-based approach

Let the infinite half-space be identified by the two coordinates $ξ$ and $η$ along its width and height, respectively, such that the medium is extended to infinity in the negative direction of $η$ . Stresses in the elastic medium can be described via a biharmonic function $ϕ (ξ, η)$ of the form

ϕ (ξ, η) = A [1 + γ η] \sin (k ξ) \exp (k η),

(8)

satisfying the biharmonic equation (5). The parameters $A$ and $γ$ can be identified from the boundary conditions and $k > 0$ denotes the wavenumber $k = 2 π / ℓ$ with $ℓ$ being the wavelength.² The stresses associated with the biharmonic function (8) read

\begin{matrix} σ_{ξ ξ} = \frac{\partial^{2} ϕ}{\partial η^{2}} = A [k^{2} [1 + γ η] + 2 k γ] \sin (k ξ) \exp (k η), \\ σ_{η η} = \frac{\partial^{2} ϕ}{\partial ξ^{2}} = - A k^{2} [1 + γ η] \sin (k ξ) \exp (k η) \\ σ_{ξ η} = - 2 \frac{\partial^{2} ϕ}{\partial ξ \partial η} = - 2 A [k^{2} [1 + γ η] + k γ] \cos (k ξ) \exp (k η), \end{matrix}

(9)

from which the strains can be obtained via the constitutive relation (3). The parameter $γ$ is next identified via imposing the vanishing $ε_{ξ ξ}$ on the surface of the infinite half-space as

ε_{ξ ξ} |_{η = 0} = 0 \Rightarrow γ = - \frac{μ_{s} + λ_{s}}{2 μ_{s} + λ_{s}} k .

(10)

Finally, the displacement field $u$ is obtained via integrating the strain field as

\begin{matrix} u_{ξ} = \frac{μ_{s} + λ_{s}}{2 μ_{s} [2 μ_{s} + λ_{s}]} k^{2} η A \cos (k ξ) \exp (k η) & \Rightarrow u_{ξ} |_{η = 0} = 0, \\ u_{η} = [\frac{μ_{s} + λ_{s}}{2 μ_{s} [2 μ_{s} + λ_{s}]} k η - \frac{3 μ_{s} + λ_{s}}{2 μ_{s} [2 μ_{s} + λ_{s}]}] kA \sin (k ξ) \exp (k η) & \Rightarrow u_{η} |_{η = 0} = - \frac{3 μ_{s} + λ_{s}}{2 μ_{s} [2 μ_{s} + λ_{s}]} kA \sin (k ξ) . \end{matrix}

(11)

The detailed derivations are omitted for the sake of space, however, the essential intermediate steps are given in Appendix A. The effective stiffness of the medium is obtained via dividing the traction in the $η$ direction by the vertical displacement on the surface as

\begin{matrix} K_{s}^{eff} = \frac{f_{η} |_{η = 0}}{u_{η} |_{η = 0}} = \frac{σ_{η η} |_{η = 0}}{u_{η} |_{η = 0}} = \frac{- A k^{2} \sin (k ξ)}{- \frac{3 μ_{s} + λ_{s}}{2 μ_{s} [2 μ_{s} + λ_{s}]} k A \sin (k ξ)} \Rightarrow K_{s}^{eff} = \frac{2 μ_{s} [2 μ_{s} + λ_{s}]}{3 μ_{s} + λ_{s}} k . \end{matrix}

(12)

2.3. Current displacement-based approach

Let $u$ denote a hypothetical stationary wave on the surface of the substrate (infinite half-space) as illustrated in Figure 1. The stationary wave lies on the plane spanned by the orthogonal basis ${t, n, m}$ where $m = t \times n$ is the outward direction along which the wave function remains constant corresponding to the associated plane-strain condition. The wave function $u$ is constructed such that it decays through the substrate meaning that its amplitude vanishes asymptotically when the depth approaches infinity. Let $ξ = x \cdot t$ and $η = x \cdot n$ denote the components of an arbitrary position vector $x$ . The stationary wave $u$ can be prescribed as the product of a decay part $w (η)$ and a waviness term $\exp (ik ξ)$ in which $k > 0$ is the wavenumber and is related to the wavelength $ℓ$ via $k = 2 π / ℓ$ . The decay part itself can be decomposed into the amplitude $a$ and decay exponent $\exp (τ η)$ thereby only $τ > 0$ is valid since $η \leq 0$ and, therefore, the wave function $u$ reads

u = w (η) \exp (ik ξ) with w (η) = a \exp (τ η) .

(13)

To proceed, the wave function $u$ is imposed as the deformation of the substrate and its admissible format satisfying the linear momentum balance is derived. The balance of linear momentum in the substrate reads $div σ = 0$ with the stress–strain relation $σ = c : ε$ . Inserting the strain $ε$ in terms of the wave function $u$ and using the symmetry properties of the constitutive tensor $c$ yields

\begin{matrix} div σ = div (c : ε) = div (c : \nabla u) = [\nabla (c : \nabla u)] : i = [c : \nabla \nabla u] : i = 0 . \end{matrix}

(14)

Using the chain rule

\nabla {•} = \frac{\partial {•}}{\partial ξ} \otimes t + \frac{\partial {•}}{\partial η} \otimes n,

the second gradient of the wave function $u$ is obtained as

\begin{matrix} u = w (η) \exp (ik ξ), \\ \nabla u = [w' \otimes n + ik w \otimes t] \exp (ik ξ), \\ \nabla \nabla u = [w ″ \otimes n \otimes n + ik [w' \otimes t \otimes n + w' \otimes n \otimes t] - k^{2} w \otimes t \otimes t] \exp (ik ξ) . \end{matrix}

(15)

Inserting $\nabla \nabla u$ into the linear momentum balance (14) yields the ordinary differential equation

\begin{matrix} A \cdot w ″ + ik B \cdot w' - k^{2} C \cdot w = 0, \end{matrix}

(16)

in which the second-order (localization) tensors $A$ , $B$ , and $C$ are defined as

\begin{matrix} A \cdot {•} : = [c : [{•} \otimes n]] \cdot n, \\ B \cdot {•} : = [c : [{•} \otimes n]] \cdot t + [c : [{•} \otimes t]] \cdot n, \\ C \cdot {•} : = [c : [{•} \otimes t]] \cdot t . \end{matrix}

(17)

Inserting $c$ from (2) and $\nabla \nabla u$ from (15)₃ into the linear momentum balance (14) yields

\begin{matrix} [μ_{s} i + [λ_{s} + μ_{s}] n \otimes n] \cdot w ″ + ik [[μ_{s} + λ_{s}] [n \otimes t + t \otimes n]] \cdot w' - k^{2} [μ_{s} i + [λ_{s} + μ_{s}] t \otimes t] \cdot w = 0 . \end{matrix}

(18)

Comparing the linear momentum balance (18) and (16), renders the localization tensors as

\begin{matrix} A : = μ_{s} i + [λ_{s} + μ_{s}] n \otimes n, \\ B : = [μ_{s} + λ_{s}] [n \otimes t + t \otimes n], \\ C : = μ_{s} i + [λ_{s} + μ_{s}] t \otimes t . \end{matrix}

(19)

Next, the differential equation (16) is solved in order to obtain the admissible forms of the wave function $u$ . In doing so, it proves convenient to rewrite the differential equation (16) as

\begin{matrix} w ″ + ik A^{- 1} \cdot B \cdot w' - k^{2} A^{- 1} \cdot C \cdot w = 0, \end{matrix}

(20)

in which $A^{- 1}$ is obtained via applying the Sherman–Morrison formula as

A^{- 1} = \frac{1}{μ_{s}} i - \frac{λ_{s} + μ_{s}}{μ_{s} [λ_{s} + 2 μ_{s}]} n \otimes n .

To proceed, the second-order differential equation (20) is decomposed into a system of two first-order differential equations as $W' = M \cdot W$ with

W = [\begin{matrix} w \\ w' \end{matrix}] \Rightarrow W' = [\begin{matrix} w' \\ w ″ \end{matrix}] = [\begin{matrix} w' \\ k^{2} A^{- 1} \cdot C \cdot w - ik A^{- 1} \cdot B \cdot w' \end{matrix}] = \underset{M}{\underset{︸}{[\begin{matrix} 0 & i \\ k^{2} A^{- 1} \cdot C & - ik A^{- 1} \cdot B \end{matrix}]}} \cdot \underset{W}{\underset{︸}{[\begin{matrix} w \\ w' \end{matrix}]}} .

(21)

Replacing the decay part $w$ in $W$ in terms of the amplitude $a$ and the decay part as $w = a \exp (τ η)$ yields

\begin{matrix} W = [\begin{matrix} w \\ w' \end{matrix}] = [\begin{matrix} a \\ τ a \end{matrix}] \exp (τ η) \Rightarrow W' = [\begin{matrix} w' \\ w ″ \end{matrix}] = [\begin{matrix} τ a \\ τ^{2} a \end{matrix}] \exp (τ η) = τ W . \end{matrix}

(22)

Comparing (21) and (22), renders the eigenproblem $M \cdot W = τ W$ or alternatively $[M - τ 1] \cdot W = 0$ whose non-trivial solutions correspond to the condition $det (M - τ 1) = 0$ , which can be represented in the matrix format as

[M - τ 1] = [\begin{matrix} - τ & 0 & 1 & 0 \\ 0 & - τ & 0 & 1 \\ k^{2} [1 + α_{s}] & 0 & - τ & - ik α_{s} \\ 0 & k^{2} [1 - β_{s}] & - ik α_{s} [1 - β_{s}] & - τ \end{matrix}] and α_{s} = \frac{λ_{s} + μ_{s}}{μ_{s}}, β_{s} = \frac{λ_{s} + μ_{s}}{λ_{s} + 2 μ_{s}},

which eventually yields

\begin{matrix} det (M - τ 1) = 0 \Rightarrow {[k^{2} - τ^{2}]}^{2} = 0 \Rightarrow τ = \pm k \overset{τ > 0}{= = = \Rightarrow} τ = k . \end{matrix}

(23)

The generic bounded solution of the differential equation (20) obtained by solving for the eigenvectors corresponding to $τ > 0$ , see Appendix B, reads

\begin{matrix} w (η) = \frac{1}{k} [[Ci - D \frac{2 μ_{s} i}{k [μ_{s} + λ_{s}]} - Di η] t + [C + D \frac{1}{k} - D η] n] \exp (k η), \end{matrix}

(24)

in which the constants $C$ and $D$ can be specified via imposing the boundary conditions. Inserting the decay function $w$ into the wave function $u$ reads

\begin{matrix} u = w (η) \exp (ik ξ) = \frac{1}{k} [[Ci - D \frac{2 μ_{s} i}{k [μ_{s} + λ_{s}]} - Di η] t + [C + D \frac{1}{k} - D η] n] \exp (k [i ξ + η]) . \end{matrix}

(25)

The relation $τ = k$ is significant since it indicates that the decay coefficient $τ$ must be identical to the wavenumber $k$ in order to have an admissible wave function (13). Therefore, the surface waves with larger wavenumbers decay faster along the depth of the substrate. In other words, the waves with larger wavelength decay slower and, thus, their effect vanishes at large depth $η$ , as shown in Figure 2. This observation is particularly remarkable from a computational perspective. In finite element simulations, the infinite half-space must be approximated as a domain with a prescribed depth. The finding (23) indicates that the finite depth of the substrate in numerical simulations can have a considerable influence on the results depending on the wavelength itself. Based on the relation $τ = k$ , from a computational perspective, one can adjust the depth of the substrate according to the wavelength so as to obtain comparable solutions.

At the surface of the substrate corresponding to $η = 0$ , the displacement is identical to that of the wrinkling film. Note, the wrinkling of the film is associated with its internal growth and the compressive residual stress thereof but no actual geometrical compression is applied. Exactly on the onset of instabilities of the film, the substrate displacement in the $ξ$ -direction is zero and in the $η$ -direction is a sinusoidal wave. The displacement in the $ξ$ -direction and $η$ -direction correspond to the $t$ -component and $n$ -component of $u$ , respectively. The hypothetical wave $u$ on the surface reads

\begin{matrix} u |_{η = 0} = \frac{1}{k} [[Ci - D \frac{2 μ_{s} i}{k [μ_{s} + λ_{s}]}] t + [C + D \frac{1}{k}] n] \exp (ik ξ), \end{matrix}

(26)

or, more specifically,

\begin{matrix} v_{ξ} : = {[u |_{η = 0}]}_{ξ} \equiv \frac{1}{k} [- C + D \frac{2 μ_{s}}{k [μ_{s} + λ_{s}]}] \sin (k ξ) and v_{η} : = {[u |_{η = 0}]}_{η} \equiv \frac{1}{k} [C + D \frac{1}{k}] \cos (k ξ) . \end{matrix}

(27)

Imposing the boundary condition $v_{ξ} = 0$ , yields the key relation

\frac{D}{C} = k \frac{μ_{s} + λ_{s}}{2 μ_{s}},

between the constants $C$ and $D$ rendering the wave function on the surface as

\begin{matrix} u |_{η = 0} = C \frac{1}{k} [\frac{3 μ_{s} + λ_{s}}{2 μ_{s}}] n \exp (ik ξ), \end{matrix}

(28)

or more specifically

\begin{matrix} v_{ξ} = {[u |_{η = 0}]}_{ξ} \equiv 0 and v_{η} = {[u |_{η = 0}]}_{η} \equiv C \frac{3 μ_{s} + λ_{s}}{2 μ_{s} k} \cos (k ξ) . \end{matrix}

(29)

Hence, the wave function $u$ reads

\begin{matrix} u = w (η) \exp (ik ξ) = \frac{1}{k} [[Ci - D \frac{2 μ_{s} i}{k [μ_{s} + λ_{s}]} - Di η] t + [C + D \frac{1}{k} - D η] n] \exp (k [i ξ + η]) \\ = C \frac{1}{k} [[i - \frac{D}{C} \frac{2 μ_{s} i}{k [μ_{s} + λ_{s}]} - \frac{D}{C} i η] t + [1 + \frac{D}{C} \frac{1}{k} - \frac{D}{C} η] n] \exp (k [i ξ + η]) \\ = C \frac{1}{k} [[- \frac{μ_{s} + λ_{s}}{2 μ_{s}} ik η] t + [\frac{3 μ_{s} + λ_{s}}{2 μ_{s}} - \frac{μ_{s} + λ_{s}}{2 μ_{s}} k η] n] \exp (k [i ξ + η]), \end{matrix}

(30)

or, alternatively,

\begin{matrix} u = C [[- \frac{μ_{s} + λ_{s}}{2 μ_{s}} i η] t + [\frac{3 μ_{s} + λ_{s}}{2 μ_{s} k} - \frac{μ_{s} + λ_{s}}{2 μ_{s}} η] n] \exp (k [i ξ + η]) . \end{matrix}

(31)

Equipped with the admissible wave function (31), the next step is to compute the traction on the surface of the substrate. The traction $f$ on the surface $η = 0$ is obtained via the Cauchy relation $f = σ \cdot n$ whereby $σ$ is the Cauchy stress on the surface as

\begin{matrix} f = σ \cdot n = σ |_{η = 0} \cdot n = [c : ε |_{η = 0}] \cdot n = [c : [\nabla u {] |}_{η = 0}] \cdot n . \end{matrix}

(32)

Therefore, to proceed, we compute the displacement gradient on the surface. The gradient of the admissible wave function (31) reads

\begin{matrix} \nabla u & = \frac{\partial u}{\partial ξ} \otimes t + \frac{\partial u}{\partial η} \otimes n = ik u \otimes t + k u \otimes n + C [- \frac{μ_{s} + λ_{s}}{2 μ_{s}} i t \otimes n - \frac{μ_{s} + λ_{s}}{2 μ_{s}} n \otimes n] \exp (k [i ξ + η]), \end{matrix}

(33)

and reduces on the surface to

\begin{matrix} [\nabla u] |_{η = 0} & = C [\frac{3 μ_{s} + λ_{s}}{2 μ_{s}} i n \otimes t - \frac{μ_{s} + λ_{s}}{2 μ_{s}} i t \otimes n + n \otimes n] \exp (ik ξ) . \end{matrix}

(34)

The surface stress is then calculated as

\begin{matrix} σ |_{η = 0} = c : [\nabla u {] |}_{η = 0} = C [μ_{s} i [n \otimes t + t \otimes n] + [2 μ_{s} + λ_{s}] n \otimes n + λ_{s} t \otimes t] \exp (ik ξ), \end{matrix}

(35)

which is notably symmetric, as expected. Finally, the surface traction $f$ reads

\begin{matrix} f = σ |_{η = 0} \cdot n = C [μ_{s} i t + [2 μ_{s} + λ_{s}] n] \exp (ik ξ), \end{matrix}

(36)

or, more specifically,

\begin{matrix} f_{ξ} : = [f]_{ξ} \equiv - C μ_{s} \sin (k ξ) and f_{η} : = [f]_{η} \equiv C [2 μ_{s} + λ_{s}] \cos (k ξ) . \end{matrix}

(37)

The effective stiffness of the substrate can be understood as the ratio $f_{η} / v_{η}$ indicating the resistance of the substrate against prescribing the wave function $v_{η}$ on its surface. Inserting the vertical surface wave $v_{η}$ from (29)₂ and the vertical surface traction $f_{η}$ from (37)₂ yields

\begin{matrix} K_{s}^{eff} = \frac{f_{η}}{v_{η}} = \frac{C [2 μ_{s} + λ_{s}] \cos (k ξ)}{C \frac{3 μ_{s} + λ_{s}}{2 μ_{s} k} \cos (k ξ)} \Rightarrow K_{s}^{eff} = \frac{2 μ_{s} [2 μ_{s} + λ_{s}]}{3 μ_{s} + λ_{s}} k . \end{matrix}

(38)

Note, the equivalent stiffness of an infinite half-space (38) using out displacement-based approach here is identical to (12) obtained from the stress-based approach of Allen employing the Airy stress functions. This conclusion is especially remarkable since the displacement fields (11) and (31) from the two approaches are different, but the equivalent stiffnesses thereof are not. Table 1 gathers the key quantities obtained from the two approaches and highlights the similarities and dissimilarities of the two formulations, see also Figure 3 for a graphical representation.

Table 1.

Analytical solution of the main fields obtained from both Allen and current approaches. Note, Allen’s solution relies on two-dimensional isotropic behavior of the domain. However, the current approach is displacement-based and can be extended to three-dimensional and anisotropic cases. Obviously, the current approach recovers the solution of Allen by setting $C = i k^{2} A / [2 μ_{s} + λ_{s}]$ .

	Allen	Current
$u_{ξ}$	$k^{2} A \frac{[μ_{s} + λ_{s}]}{2 μ_{s} [2 μ_{s} + λ_{s}]} η \cos (k ξ) \exp (k η)$	$- i C \frac{[μ_{s} + λ_{s}]}{2 μ_{s}} η \cos (k ξ) \exp (k η)$
$u_{η}$	$k A [\frac{[μ_{s} + λ_{s}]}{2 μ_{s} [2 μ_{s} + λ_{s}]} k η - \frac{3 μ_{s} + λ_{s}}{2 μ_{s} [2 μ_{s} + λ_{s}]}] \sin (k ξ) \exp (k η)$	$- i C [\frac{[μ_{s} + λ_{s}]}{2 μ_{s}} η - \frac{3 μ_{s} + λ_{s}}{2 μ_{s} k}] \sin (k ξ) \exp (k η)$
$ε_{ξ ξ}$	$- k^{2} A [\frac{[μ_{s} + λ_{s}]}{2 μ_{s} [2 μ_{s} + λ_{s}]} k η] \sin (k ξ) \exp (k η)$	$i C [\frac{[μ_{s} + λ_{s}]}{2 μ_{s}} k η] \sin (k ξ) \exp (k η)$
$ε_{η η}$	$k^{2} A [\frac{[μ_{s} + λ_{s}]}{2 μ_{s} [2 μ_{s} + λ_{s}]} k η - \frac{1}{2 μ_{s} + λ_{s}}] \sin (k ξ) \exp (k η)$	$- i C [\frac{[μ_{s} + λ_{s}]}{2 μ_{s}} k η - 1] \sin (k ξ) \exp (k η)$
$σ_{ξ ξ}$	$- k^{2} A [\frac{[μ_{s} + λ_{s}]}{2 μ_{s} + λ_{s}} k η + \frac{λ_{s}}{2 μ_{s} + λ_{s}}] \sin (k ξ) \exp (k η)$	$i C [[μ_{s} + λ_{s}] k η + λ_{s}] \sin (k ξ) \exp (k η)$
$σ_{η η}$	$k^{2} A [\frac{[μ_{s} + λ_{s}]}{2 μ_{s} + λ_{s}} k η - 1] \sin (k ξ) \exp (k η)$	$- i C [[μ_{s} + λ_{s}] k η - [2 μ_{s} + λ_{s}]] \sin (k ξ) \exp (k η)$

2.4. Effective stiffness under plane-strain and plane-stress condition

All the derivations so far correspond to the plane-strain condition. Nevertheless, we can adopt the current results for plane-stress instead of plane-strain, however, the detailed derivations are omitted for brevity. In particular, the effective stiffnesses of the substrate $K_{s}^{eff}$ against a sinusoidal wave with the wavenumber $k = 2 π / ℓ$ on its surface for the two cases read

\begin{matrix} Plane Strain : K_{s}^{eff} = \frac{2 μ_{s} [2 μ_{s} + λ_{s}]}{3 μ_{s} + λ_{s}} k, Plane Stress : K_{s}^{eff} = \frac{8 μ_{s} [μ_{s} + λ_{s}]}{6 μ_{s} + 5 λ_{s}} k, \end{matrix}

(39)

in terms of the Lamé parameters $μ_{s}$ and $λ_{s}$ . It is particularly notable that in the fully compressible limit associated with $λ_{s} = 0$ , the effective stiffnesses $K_{s}^{eff}$ for both cases coincide as

\begin{matrix} λ_{s} = 0 \Rightarrow Plane Strain : K_{s}^{eff} = \frac{4}{3} μ_{s} k, Plane Stress : K_{s}^{eff} = \frac{4}{3} μ_{s} k . \end{matrix}

(40)

On the other hand, in the incompressible limit associated with $λ_{s} \to \infty$ , the effective stiffnesses $K_{s}^{eff}$ for the two cases reach their maximum difference as

\begin{matrix} λ_{s} \to \infty \Rightarrow Plane Strain : K_{s}^{eff} = 2 μ_{s} k, Plane Stress : K_{s}^{eff} = \frac{8}{5} μ_{s} k . \end{matrix}

(41)

In the incompressible limit, the effective stiffness associated with the plane-stress case underestimates its counterpart in the plane-strain case by 20%. This observation is particularly important and relevant to the applications of the current study since such instabilities usually occur in soft materials and biological tissues that are often nearly incompressible. Nevertheless, this subtle and yet important nuance is repeatedly overlooked in the literature to date. Figure 4 illustrates the profile of the effective stiffness versus Lamé parameters. It can be seen that the effective stiffness under plane-stress condition is consistently less than its counterpart under the plane-strain condition. This can be explained due to the decreased resistance of the substrate in plane-stress condition due to its additional flexibility in the lateral direction. Exactly for the same reason, in the fully compressible limit both plane-strain and plane-stress formulations render identical results.

The problem formulation so far has been in terms of the Lamé parameters $μ_{s}$ and $λ_{s}$ . Nonetheless, it is possible to carry out the derivations in terms of other pairs of parameters such as the elastic modulus $E_{s}$ and Poisson’s ratio $ν_{s}$ . A remarkable outcome is that while $K_{s}^{eff}$ in terms of the Lamé parameters differs for plane-strain and plane-stress, translating $K_{s}^{eff}$ in terms of $E_{s}$ and $ν_{s}$ renders identical forms for both cases as

\begin{matrix} K_{s}^{eff} (E_{s}, ν_{s}) = \frac{2 E_{s}}{[3 - ν_{s}] [1 + ν_{s}]} k . \end{matrix}

(42)

This finding is significant since it provides a unified formulation for both the plane-strain and plane-stress conditions. Consequently, we formulate the remainder of this manuscript in terms of the elastic modulus and Poisson’s ratio instead of the Lamé parameters. Note, the use of the Lamé parameters until this point was advantageous since otherwise the format of the localization tensors (19) would have been more complicated.

We mention in passing that care should be taken when calculating the material parameters in two dimensions since the relations between the material parameters are not necessarily identical to those in three dimensions and even the physical bounds for the same parameter could be different. For instance, the Poisson’s ratio could reach 1 in the incompressibility limit for the plane-strain case unlike 0.5 for the three-dimensional case. Table 2 gathers the relations between different material parameters for both the plane-strain as well as plane-stress scenarios and provides further insight. The conventional three-dimensional elastic modulus and Poisson’s ratio are denoted ${{}^{3 D}E}_{s}$ and ${{}^{3 D}ν}_{s}$ , respectively, so as to clearly distinguish them from their two-dimensional counterparts.

Table 2.

Effective stiffness of the substrate $K_{s}^{eff}$ against a sinusoidal wave with the wavenumber $k = 2 π / ℓ$ on its surface. The two different formats of $K_{s}^{eff}$ in terms of the Lamé material parameters $μ_{s}$ and $λ_{s}$ coincide in terms of $E_{s}$ and $ν_{s}$ . The elastic modulus and Poisson’s ratio associated with conventional three-dimensional tests are denoted ${{}^{3 D}E}_{s}$ and ${{}^{3 D}ν}_{s}$ , respectively.

Plane-strain		Plane-stress
$K_{s}^{eff} (μ_{s}, λ_{s}) = \frac{2 μ_{s} [2 μ_{s} + λ_{s}]}{3 μ_{s} + λ_{s}} k$		$K_{s}^{eff} (μ_{s}, λ_{s}) = \frac{8 μ_{s} [μ_{s} + λ_{s}]}{6 μ_{s} + 5 λ_{s}} k$
$E_{s} = \frac{4 μ_{s} [μ_{s} + λ_{s}]}{2 μ_{s} + λ_{s}}$	$ν_{s} = \frac{λ_{s}}{2 μ_{s} + λ_{s}}$	$E_{s} = \frac{μ_{s} [2 μ_{s} + 3 λ_{s}]}{μ_{s} + λ_{s}}$	$ν_{s} = \frac{λ_{s}}{2 [μ_{s} + λ_{s}]}$
$E_{s} = \frac{{{}^{3 D}E}_{s}}{1 -^{3 D} ν_{s}^{2}}$	$ν_{s} = \frac{{{}^{3 D}ν}_{s}}{1 -^{3 D} ν_{s}}$	$E_{s} =^{3 D} E_{s}$	$ν_{s} =^{3 D} ν_{s}$
$μ_{s} = \frac{E_{s}}{2 [1 + ν_{s}]}$	$λ_{s} = \frac{E_{s} ν_{s}}{[1 - ν_{s}] [1 + ν_{s}]}$	$μ_{s} = \frac{E_{s}}{2 [1 + ν_{s}]}$	$λ_{s} = \frac{E_{s} ν_{s}}{[1 - 2 ν_{s}] [1 + ν_{s}]}$
$K_{s}^{eff} (E_{s}, ν_{s}) = \frac{2 E_{s}}{[3 - ν_{s}] [1 + ν_{s}]} k$		$K_{s}^{eff} (E_{s}, ν_{s}) = \frac{2 E_{s}}{[3 - ν_{s}] [1 + ν_{s}]} k$

3. Wrinkling of a growing layer

This section details on growth-induced instabilities of a thin growing film. To begin with, the thin film is assumed to be on top of a deep (infinite half-space) compliant substrate, as illustrated in Figure 5. The analytical solution of this problem is derived in Section 3.1 based on the equivalent stiffness of the substrate (42). The analytical solution is illustrated next using numerical examples in Section 3.2 and its key features are discussed. Section 3.3 compares the analytical solution with computational simulations using the finite element method. Finally, Section 3.4 investigates the extended bilayer wrinkling and elaborates on how to generalize the proposed analytical solution to more complicated scenarios.

Note that the origin of the compressive stresses leading to instabilities of the film need not necessarily be differential growth. For instance, pre-stretch or thermal expansion could also result in such geometric instabilities. In fact, the nature of the compressive stresses could influence geometric instabilities in bilayers and this issue has been carefully analyzed in [13, 14, 98] very recently. Based on these studies, the origin of the compressive stresses has more significant effects on post-buckling behavior and secondary instabilities but it has very little influence on the primary buckling modes that are of particular interest here. More precisely, one underlying assumption of this manuscript is that at the onset of wrinkling the substrate is stress-free, which shall be revisited if other scenarios are to be considered. Various extensions of the proposed approach to more complicated cases and to introduce pre-stretch as well as nonlinearities into the picture are possible, but shall be pursued in separate contributions.

In the subsequent derivations, the thickness of the domain in the out-of-plane direction is assumed to be uniform and denoted by $b$ . The problem can be modeled as being plane-strain or plane-stress in two dimensions. Based on the discussion in Section 2.4 both formulations would formally coincide if the substrate parameters $E_{s}$ and $ν_{s}$ are defined according to Table 2 and used instead of the Lamé parameters $μ_{s}$ and $λ_{s}$ . Therefore, in what follows, the derivations are carried out in terms of the elastic modulus and Poisson’s ratio instead of the Lamé parameters. The boundary conditions are illustrated in Figure 5 (right) indicating that the horizontal displacements along the length of the domain are prohibited but the vertical displacements remain free.

3.1. Analytical solution

In order to derive the analytical solution for wrinkling of a growing thin film on an infinite half-space, we decompose the domain into two subdomains, namely, the film and the substrate, as shown in Figure 6. Then each subdomain is treated separately and their corresponding deformations are superimposed according to the geometrical constraint of perfect bonding between the film and the substrate. That is, the film is constantly adhered to the substrate and cannot detach from it and, thus, the deformation of the film must be identical to that of the substrate. The elastic modulus and Poisson’s ratio of the film are denoted by $E_{f}$ and $ν_{f}$ , respectively. The growing elastic film fixed at both ends resembles a thin beam under compression wherein the compression results from the accumulation of the residual stresses due to growth. Furthermore, the interaction of the film with the substrate manifests itself as a distributed force on the film. If $v$ denotes the deflection of a thin beam, the classical governing equation of the beam reads

\begin{matrix} E_{f} I v ″ ″ + σ b h v ″ = f, \end{matrix}

(43)

in which $E_{f}$ and $I$ are the elastic modulus and the moment of inertia of the beam, respectively. The compressive stress along the film due to growth is denoted by $σ$ and $f$ on the right-hand side denotes the distributed force on the beam in the orthogonal direction. Note, positive values for $σ$ correspond to compression based on our definition. The fundamental solution of the film equation (43) for the deflection $v$ is of sinusoidal form. On the other hand, for a sinusoidal deflection function on the surface of the substrate, the force distribution reads

\begin{matrix} f = - b K_{s}^{eff} v with K_{s}^{eff} = \frac{2 E_{s}}{[3 - ν_{s}] [1 + ν_{s}]} k, \end{matrix}

(44)

where $k = 2 π / ℓ$ is the wavenumber of the deflection $v$ . Clearly, the relation (44) establishes link between the findings of Section 2 and the derivations in this section. Furthermore, the minus sign in front of $f$ is due to the fact that the force $f$ in (44) is the (reaction) force exerted from the substrate on the film hence, negative of the force on the surface of the substrate.

Next, we insert the distributed force (44)₁ into the governing equation of the film (43). Also, the moment of inertia for a film with a rectangular cross-section of the width $b$ and height $h$ is $I = b h^{3} / 12$ . Finally, considering the sinusoidal nature of the deflection, the relation $v ″ = - k^{2} v$ holds and yields

\begin{matrix} \frac{1}{12} E_{f} h^{3} k^{4} v - σ h k^{2} v = - \frac{2 E_{s}}{[3 - ν_{s}] [1 + ν_{s}]} k v \Rightarrow σ = \frac{1}{12} E_{f} h^{2} k^{2} + \frac{2 E_{s}}{[3 - ν_{s}] [1 + ν_{s}]} \frac{1}{hk}, \end{matrix}

(45)

for the compressive stress $σ$ in the film due to growth. In order to calculate the critical wavenumber $k_{cr}$ at the onset of the growth-induced instability, we set the derivative of the stress $σ$ with respect to $k$ to zero and after some manipulation $k_{cr}$ reads

\begin{matrix} \frac{\partial σ}{\partial k} = 0 \Rightarrow k_{cr} = \frac{1}{h} {[\frac{12 E_{s}}{E_{f} [3 - ν_{s}] [1 + ν_{s}]}]}^{\frac{1}{3}}, \end{matrix}

(46)

from which the critical induced stress in the film $σ_{cr}$ and the associated critical strain $ε_{cr} = σ_{cr} / E_{f}$ can be obtained by inserting $k_{cr}$ into (45) as

\begin{matrix} σ_{cr} = E_{f} {[\frac{3 E_{s}}{2 E_{f} [3 - ν_{s}] [1 + ν_{s}]}]}^{\frac{2}{3}}, ε_{cr} = {[\frac{3 E_{s}}{2 E_{f} [3 - ν_{s}] [1 + ν_{s}]}]}^{\frac{2}{3}} . \end{matrix}

(47)

Eventually, the critical wavelength $ℓ_{cr}$ and the critical growth $g_{cr}$ can be computed as $ℓ_{cr} = 2 π / k_{cr}$ and $g_{cr} = ε_{cr} / [1 - ε_{cr}]$ thereby a positive strain is assumed to be compressive. Note, the final results here are formally identical to those given by Allen [1]. Nonetheless, the material parameters in (46) and (47) must be related to the conventional three-dimensional material parameters according to Table 3. Thus, the pioneering Allen’s relations in their original form are valid only for the plane-stress case, but not for the plane-strain case unless the transformations according to Table 3 are taken into account. This is particularly remarkable since the associated computational simulations using the finite element method to date are often based on the plane-strain condition but are mistakenly compared with the Allen’s solution.

Table 3.

Critical wavelength $ℓ_{cr}$ and critical stress $σ_{cr}$ in the film to initiate growth-induced instabilities at the critical growth $g_{cr}$ . The elastic modulus and Poisson’s ratio associated with conventional three-dimensional tests are denoted by ${{}^{3 D}E}_{s}$ and ${{}^{3 D}ν}_{s}$ , respectively.

$ℓ_{cr} = 2 π h {[\frac{E_{f} [3 - ν_{s}] [1 + ν_{s}]}{12 E_{s}}]}^{\frac{1}{3}}$ , $ε_{cr} = {[\frac{3 E_{s}}{2 E_{f} [3 - ν_{s}] [1 + ν_{s}]}]}^{\frac{2}{3}}$ , $g_{cr} = \frac{ε_{cr}}{1 - ε_{cr}}$
plane-strain	plane-stress
$E_{s} = \frac{{{}^{3 D}E}_{s}}{1 -^{3 D} ν_{s}^{2}}$ , $ν_{s} = \frac{{{}^{3 D}ν}_{s}}{1 -^{3 D} ν_{s}}$	$E_{s} =^{3 D} E_{s}$ , $ν_{s} =^{3 D} ν_{s}$
$E_{f} = \frac{{{}^{3 D}E}_{f}}{1 -^{3 D} ν_{f}^{2}}$	$E_{f} =^{3 D} E_{f}$

We emphasize that the relation between the compressive strain in the film and the film growth reads $g = ε / [1 - ε]$ , or alternatively $ε = g / [1 + g]$ , in which $g$ is the (anisotropic) film growth coefficient along the film itself with $g > 0$ indicating growth and $g < 0$ corresponding to shrinkage. In addition, throughout the previous derivations it is assumed that a positive strain is compressive, which complies with our sign convention for growth $g$ . In other words, shrinkage of the film results in negative (tensile) strain and growth of the film leads to positive (compressive) strain in the film. In view of a finite deformation setting, the growth tensor can be constructed via $F_{g} = [1 + g] \bar{I}$ in which $\bar{I} = I - N \otimes N$ is the identity tensor along the film with $N$ being the unit normal vector to the film.

3.2. Numerical illustrations

The analytical expressions obtained in Section 3.1 are illustrated here numerically to better interpret them and use them to their fullest extent in numerous applications. Figure 7 gathers the graphs for the critical growth $g_{cr}$ and the critical wavelength $ℓ_{cr}$ at the incompressibility limit and provides a thorough comparison between the plane-stress and plane-strain conditions. The critical growth is plotted against the film to substrate stiffness ratio $E_{f} / E_{s}$ as well as the film thickness $h$ . Note that throughout the derivations in the previous section, it was assumed that the growing thin film is considerably stiffer than the substrate and, therefore, our results are more accurate for $E_{f} / E_{s} > 10$ . In other words, the validity of our approach for the stiffness ratios less than 10 is debatable since the compressive strain $ε$ can no longer be regarded as infinitesimal, see [13, 99] for further details and discussions. However, for the sake of completeness, we provide the graphs for stiffness ratios as low as $E_{f} / E_{s} = 1$ . As is illustrated in Figure 9, even for such small stiffness ratios, we observe a very good agreement between the analytical solution and the numerical results using the finite element method. This surprising agreement between the numerical and analytical results can be justified since the compressive stresses within the film result from growth only and not compression on the whole domain. This in turn corresponds to a stress-free substrate at the onset of wrinkling behaving according to small-strain linear elasticity. See [99] for a thorough and systematic study on instabilities of (soft) films on (compliant) substrates.

The surface plots clearly illustrate that the critical growth $g_{cr}$ is independent of the film thickness and is only a nonlinear function of the stiffness ratio. While the difference between plane-stress and plane-strain for the effective stiffness of the substrate could reach up to 20%, see Figure 4, this difference reduces to less than 2% for the critical growth $g_{cr}$ , especially for higher stiffness ratios. Nonetheless, it is crucial to recognize the difference between the two approaches even though for this specific problem they happen to render nearly identical results. This difference is subtle yet noteworthy since the discrepancies encountered when comparing analytical and computational results are generally due to inconsistent plane-stress or plane-strain assumptions. However, such mismatches are to date frequently dismissed as being numerical errors. The same discussion holds for the critical wavelength $ℓ_{cr}$ , as well.

The critical wavelength $ℓ_{cr}$ is plotted against film thickness $h$ and the third root of the stiffness ratio $E_{f} / E_{s}$ . The linear relationship between the critical wavelength $ℓ_{cr}$ and the third root of the stiffness ratio is portrayed clearly by the plots. The subtle difference in the critical wavelength arising from plane-stress or plane-strain assumption is highlighted in the center panel by plotting the critical wavelength $ℓ_{cr}$ normalized by thickness $h$ against the third root of the stiffness ratio for both conditions. As the film thickness $h$ increases for a fixed stiffness ratio, the number of waves at the onset of wrinkling decreases or, more precisely, the critical wavelength varies proportionally to the film thickness, as expected. In addition, the correlation between the critical wavelength $ℓ_{cr}$ and the stiffness ratio is observed by moving down the stiffness ratio axis at a certain thickness. Such thorough representations are extremely helpful in determining the design parameters to fulfill the requirements in various applications and to better understand the interplay between the growth rate and the associated instability patterns.

3.3. Comparison with finite element method

The objective of this section is to compare briefly the analytical solution gathered in Table 3 with the numerical results obtained from computational simulations using the finite element method via an eigenvalue analysis. The comparisons given here are associated with the plane-strain condition. Nonetheless, we have performed a similar set of examples for the plane-stress condition in accordance with Section 3.2 but omitted them here since they mainly lead to the same observations and conclusions without providing any further insight. In the analytical approach, the geometrical instabilities are implicitly accounted for via buckling analysis of the film. However, the computational simulations are based on the nonlinear finite strain theory of continuum mechanics where the geometrical nonlinearity is absolutely crucial to capture geometric instabilities.

Figure 8 illustrates the kinematics of growth and the associated instabilities within the framework of nonlinear continuum mechanics. In this framework, the material (reference) configuration $B_{0}$ is mapped onto the spatial (current) configuration $B_{t}$ via the nonlinear deformation map $φ$ . That is, the placement of the spatial points $x$ are related to their material counterparts $X$ via $x = φ (X)$ . The corresponding linear deformation map $F : = Grad φ$ relates the line element $d x$ in the spatial configuration to its material counterpart $d X$ via the relation $d x = F \cdot d X$ .

To account for growth, we employ the commonly accepted concept of the multiplicative decomposition of the deformation gradient into its elastic and growth part as $F = F_{e} \cdot F_{g}$ , see [51, 97, 100, 101] among others. Clearly, the grown (intermediate) configuration $B_{g}$ need not be necessarily compatible. The stored energy in the material is not altered due to morphogenetic nature of the growth and, thus, the free energy of a growing continuum would only be a function of the elastic part of the deformation as $ψ = ψ (F) = ψ_{e} (F_{e})$ where $ψ$ and $ψ_{e}$ denote the total free energy density and the elastic free energy density, respectively. Thus, the total Piola stress $P = \partial ψ / \partial F$ can be calculated from the elastic Piola stress $P_{e} = \partial ψ_{e} / \partial F_{e}$ via the relation $P = P_{e} \cdot {F_{g}}^{- t}$ . Next, similar to conventional finite deformation elasticity and upon prescribing an elastic free energy density $ψ_{e}$ , the Piola stress field satisfying the momentum balance and subsequently the associated deformation field is calculated. Obviously, additional care should be taken due to the geometrical instabilities inherent in this problem. For the details of the finite element implementation for growth-induced instabilities, see [102] and references therein.

Figure 9 gathers an exhaustive comparison between the numerical results using the finite element method and the analytical solution for a broad range of film to substrate stiffness ratios denoted $E_{f} / E_{s}$ . The comparative studies are carried out for different film thicknesses and also for three different Poisson’s ratios. For all the examples, it is assumed that the Poisson’s ratio of the film and substrate are identical and, therefore, the stiffness ratio $E_{f} / E_{s}$ remains independent of the Poisson’s ratio. To ensure a meaningful comparison, identical three-dimensional parameters are chosen for both the analytical solution as well as the computational simulations, see Table 3. For each Poisson’s ratio and film thickness, both the critical growth $g_{cr}$ and the critical wavelength $ℓ_{cr}$ are plotted. Motivated from the analytical solution, the horizontal axis for the critical wavelength $ℓ_{cr}$ is chosen as the third root of the stiffness ratio so as to recover the linear dependence. The solid lines indicate the analytical solution while the points are the results of the computational simulations using the finite element method in accordance with [102], see also [103]. Unlike for the critical growth $g_{cr}$ , the numerical results for the critical wavelength $ℓ_{cr}$ exhibit a step-wise behavior due to the fact that the underlying finite element discretization can only comply with multiples of half-wavelength but cannot be infinitely refined to continuously conform with the larger wavelengths corresponding to larger stiffness ratios. Overall, the analytical solution and computational results are in excellent agreement and notably so for larger Poisson’s ratio. This is particularly important since most of the applications of the current study deal with biological soft tissues that are usually nearly incompressible.

3.4. Extended bilayer wrinkling

A key feature of this contribution is to decompose a bilayer structure into a substrate and a film and subsequently derive the equivalent stiffness of the substrate against the displacement of the film and, therefore, regard the substrate as a nonlinear spring with the effective stiffness $K_{s}^{eff}$ . Introducing the notion of effective stiffness for the substrate does not fundamentally change the derivations for a simple bilayer structure, however, it unravels the nature of the problem and paves the way to extend the analytical solution to more complicated scenarios such as extended bilayer wrinkling [40, 41], trilayer wrinkling [21], or multilayer wrinkling [46], as schematically illustrated in Figure 10. See also the extended reference list in the introduction for further details.

For instance, the extended bilayer wrinkling shown in Figure 10 (left) is composed of a substrate and superstrate with the material properties and geometries shown in Figure 11 for which a similar analysis to that in Section 2.3 can be carried out. However, employing the notion of equivalent stiffness allows us to readily interpret the extended bilayer system as a combination of two springs with the stiffnesses $K_{s}^{eff}$ and $K_{u}^{eff}$ replacing the substrate and superstrate, respectively. In this case, the displacements for both springs are identical but their effective forces on the film are added leading to a spring combination in parallel. Thus, the total stiffness of the system against a sinusoidal displacement of the film reads $K_{tot}^{eff} = K_{s}^{eff} + K_{u}^{eff}$ . On the other hand, the trilayer wrinkling depicted in Figure 10 (center) is different in that the intermediate layer with the effective stiffnesses $K_{i}^{eff}$ lies between the film and the substrate corresponding to a spring combination in series that yields $1 / K_{tot}^{eff} = 1 / K_{s}^{eff} + 1 / K_{i}^{eff}$ . Obviously, this approach can be applied to multilayer wrinkling of Figure 10 (right) or other combinations of multiple layers.

4. Concluding remarks

Growth often plays a crucial role in the behavior of living systems. Here we have presented our first attempt to provide a displacement-based approach to analyze geometric instabilities in bilayer structures. The most commonly accepted solution strategy to identify the critical conditions to initiate such instabilities dates back to the seminal work of Allen [1] and is based on Airy stress functions. The Allen solution is limited to two-dimensional linear isotropic problems and its success depends entirely on choosing an appropriate Airy function. This contribution circumvents the limitations associated with the Allen solution via a displacement-based approach that is intrinsically suitable for this problem. Moreover, the subtle and yet frequently overlooked difference between the solutions corresponding to the plane-strain and plane-stress conditions has been carefully analyzed. Through a series of numerical examples, the analytical solution has been compared against computational simulations using the finite element method via eigenvalue analysis and, overall, an excellent agreement has been observed. Finally, it has been briefly explained how the current approach can be utilized in various scenarios such as extended bilayer wrinkling and trilayer wrinkling. In particular, introducing the notion of effective stiffness here provides a great insight into instabilities of bilayer systems and equips us with a powerful methodology to deal with more complicated cases beyond the classical bilayer structures. Our next immediate plan is to extend the proposed strategy to study geometric instabilities accounting for nonlinearities and secondary patterns. In summary, this manuscript presents our first attempt to shed light on geometric instabilities in bilayers via a displacement-based approach inherently suited for this problem, instead of the stress-based approach of Allen. We believe that this framework can significantly enhance our understanding of geometric instabilities greatly pertinent to soft bio-materials and living systems.

Footnotes

Appendix A. Derivations of stress formulation

This section provides detailed derivations and steps regarding the approach presented in Section 2.2. While the main goal here is to include the intermediate steps in the derivations, only the crucial steps are listed and the trivial ones are omitted for the sake of brevity. The derivatives of the stress function

\begin{matrix} ϕ (ξ, η) = A [1 + γ η] \sin (k ξ) \exp (k η), \end{matrix}

read

(48)

\begin{matrix} \frac{\partial ϕ}{\partial η} = A [k [1 + γ η] + γ] \sin (k ξ) \exp (k η), \\ \frac{\partial^{2} ϕ}{\partial η^{2}} = A [k^{2} [1 + γ η] + 2 k γ] \sin (k ξ) \exp (k η), \\ \frac{\partial ϕ}{\partial ξ} = A k [1 + γ η] \cos (k ξ) \exp (k η), \\ \frac{\partial^{2} ϕ}{\partial ξ^{2}} = - A k^{2} [1 + γ η] \sin (k ξ) \exp (k η), \\ \frac{\partial^{2} ϕ}{\partial ξ \partial η} = A [k^{2} [1 + γ η] + k γ] \cos (k ξ) \exp (k η), \end{matrix}

from which the stresses (6) can be obtained. Next, the strains are computed from the stress field and, eventually, the displacement field $u$ is derived via integrating the strain field as

(49)

\begin{matrix} u_{ξ} = \int ε_{ξ ξ} d ξ and u_{η} = \int ε_{η η} d η, \end{matrix}

subjected to the boundary conditions. The strains are related to the stresses according to the relations

(50)

\begin{matrix} ε_{ξ ξ} = \frac{2 μ_{s} + λ_{s}}{4 μ_{s} [μ_{s} + λ_{s}]} σ_{ξ ξ} - \frac{λ_{s}}{4 μ_{s} [μ_{s} + λ_{s}]} σ_{η η}, \\ ε_{η η} = \frac{2 μ_{s} + λ_{s}}{4 μ_{s} [μ_{s} + λ_{s}]} σ_{η η} - \frac{λ_{s}}{4 μ_{s} [μ_{s} + λ_{s}]} σ_{ξ ξ}, \end{matrix}

in which the stresses are

(51)

\begin{matrix} σ_{ξ ξ} = A [k^{2} [1 + γ η] + 2 k γ] \sin (k ξ) \exp (k η), \\ σ_{η η} = - A k^{2} [1 + γ η] \sin (k ξ) \exp (k η) . \end{matrix}

After several mathematical steps, the displacement field $u$ reads

(52)

\begin{matrix} u_{ξ} = - [\frac{1}{2 μ_{s}} k [1 + γ η] + \frac{2 μ_{s} + λ_{s}}{2 μ_{s} [μ_{s} + λ_{s}]} γ] A \cos (k ξ) \exp (k η), \\ u_{η} = - [\frac{1}{2 μ_{s}} k [1 + γ η] - \frac{1}{2 [μ_{s} + λ_{s}]} γ] A \sin (k ξ) \exp (k η) . \end{matrix}

The restriction of the displacement field to the surface of the infinite half-space $η = 0$ yields

(53)

\begin{matrix} u_{ξ} = - [\frac{1}{2 μ_{s}} k + \frac{2 μ_{s} + λ_{s}}{2 μ_{s} [μ_{s} + λ_{s}]} γ] A \cos (k ξ), \\ u_{η} = - [\frac{1}{2 μ_{s}} k - \frac{1}{2 [μ_{s} + λ_{s}]} γ] A \sin (k ξ) . \end{matrix}

On the other hand, the stresses on the surface $η = 0$ read

(54)

\begin{matrix} σ_{ξ ξ} |_{η = 0} = A [k^{2} + 2 k γ] \sin (k ξ), \\ σ_{η η} |_{η = 0} = - A k^{2} \sin (k ξ), \\ σ_{ξ η} |_{η = 0} = - 2 A [k^{2} + k γ] \cos (k ξ), \end{matrix}

resulting in the surface strain along the surface

(55)

\begin{matrix} ε_{ξ ξ} |_{η = 0} = \frac{2 μ_{s} + λ_{s}}{4 μ_{s} [μ_{s} + λ_{s}]} σ_{ξ ξ} |_{η = 0} - \frac{λ_{s}}{4 μ_{s} [μ_{s} + λ_{s}]} σ_{η η} |_{η = 0} = \frac{k}{2 μ_{s} [μ_{s} + λ_{s}]} [[μ_{s} + λ_{s}] k + γ [2 μ_{s} + λ_{s}]] A \sin (k ξ), \end{matrix}

which must vanish identically and, thus, furnishing the parameter $γ$ as

(56)

\begin{matrix} ε_{ξ ξ} |_{η = 0} = 0 \Rightarrow γ = - \frac{μ_{s} + λ_{s}}{2 μ_{s} + λ_{s}} k . \end{matrix}

Using $γ$ from (56), the stresses read

(57)

\begin{matrix} σ_{ξ ξ} = - [\frac{μ_{s} + λ_{s}}{2 μ_{s} + λ_{s}} k η + \frac{λ_{s}}{2 μ_{s} + λ_{s}}] A k^{2} \sin (k ξ) \exp (k η), \\ σ_{η η} = [\frac{μ_{s} + λ_{s}}{2 μ_{s} + λ_{s}} k η - 1] A k^{2} \sin (k ξ) \exp (k η), \\ σ_{ξ η} = 2 [\frac{μ_{s} + λ_{s}}{2 μ_{s} + λ_{s}} k η - \frac{μ_{s}}{2 μ_{s} + λ_{s}}] A k^{2} \cos (k ξ) \exp (k η) . \end{matrix}

Furthermore, it is possible to compute the surface strain normal to the surface as well as the surface shear strain as

(58)

\begin{matrix} ε_{η η} |_{η = 0} = \frac{2 μ_{s} + λ_{s}}{4 μ_{s} [μ_{s} + λ_{s}]} σ_{η η} |_{η = 0} - \frac{λ_{s}}{4 μ_{s} [μ_{s} + λ_{s}]} σ_{ξ ξ} |_{η = 0} = - \frac{1}{2 μ_{s} + λ_{s}} k^{2} A \sin (k ξ), \\ ε_{ξ η} |_{η = 0} = \frac{1}{2 μ_{s}} σ_{ξ η} |_{η = 0} = \frac{1}{[2 μ_{s} + λ_{s}]} k^{2} A \cos (k ξ) . \end{matrix}

Finally, the parameter $γ$ from (56) is replaced in the displacement field (52) resulting in

(59)

\begin{matrix} u_{ξ} = \frac{μ_{s} + λ_{s}}{2 μ_{s} [2 μ_{s} + λ_{s}]} k^{2} η A \cos (k ξ) \exp (k η), \\ u_{η} = [\frac{μ_{s} + λ_{s}}{2 μ_{s} [2 μ_{s} + λ_{s}]} k η - \frac{3 μ_{s} + λ_{s}}{2 μ_{s} [2 μ_{s} + λ_{s}]}] k A \sin (k ξ) \exp (k η), \end{matrix}

whose restriction on the surface reads

(60)

\begin{matrix} u_{ξ} |_{η = 0} = 0 and u_{η} |_{η = 0} = - \frac{3 μ_{s} + λ_{s}}{2 μ_{s} [2 μ_{s} + λ_{s}]} k A \sin (k ξ) . \end{matrix}

Appendix B. Derivation of the general solution

The general solution to (16) may be obtained by solving for the eigenvectors corresponding to the eigenvalue obtained in (23). Solving for $w$ in the following system renders the eigenvector corresponding to the eigenvalue $k$

[M - k 1] = 0 \Rightarrow [\begin{matrix} - k & 0 & 1 & 0 \\ 0 & - k & 0 & 1 \\ k^{2} [1 + α_{s}] & 0 & - k & - ik α_{s} \\ 0 & k^{2} [1 - β_{s}] & - ik α_{s} [1 - β_{s}] & - k \end{matrix}] [\begin{matrix} v_{1} \\ v_{2} \\ v_{3} \\ v_{4} \end{matrix}] = 0 .

This system can be solved using three of the four equations. All the possible solutions to this system are linearly dependent, resulting in the relation $v_{1} = i v_{2} = v_{3} / k = i v_{4} / k$ for which setting $v_{1} = 1$ renders the eigenvector

v^{*} = {[\begin{matrix} 1 & - i & k & - ik \end{matrix}]}^{t} .

This gives the first fundamental solution $W_{1} = v \exp (k η)$ of the differential equation $W' = M \cdot W$ . Since the multiplicity of the eigenvalue $k$ is 2, but only one corresponding eigenvector is obtained (i.e. the geometric and algebraic multiplicities of the eigenvalue are different), another vector is needed for the second fundamental solution. This second fundamental solution is of the form $W_{2} = v η \exp (k η) + v^{*} \exp (k η)$ . Combined with the first fundamental solution, the general solution for the system $W' = M \cdot W$ can, hence, be written as

(61)

\begin{matrix} W = W_{1} + W_{2} = v \exp (k η) + v η \exp (k η) + v^{*} \exp (k η) . \end{matrix}

Placing the proposed $W$ from (61) into $W' = M \cdot W$ , after some manipulations, yields the relations $[M - k 1] v = 0$ and $[M - k 1] v^{*} = v$ . The first equation has already been solved for when computing the eigenvector corresponding to the eigenvalue of $k$ . The second equation presents the system

[\begin{matrix} - k & 0 & 1 & 0 \\ 0 & - k & 0 & 1 \\ k^{2} [1 + α_{s}] & 0 & - k & - ik α_{s} \\ 0 & k^{2} [1 - β_{s}] & - ik α_{s} [1 - β_{s}] & - k \end{matrix}] [\begin{matrix} v_{1}^{*} \\ v_{2}^{*} \\ v_{3}^{*} \\ v_{4}^{*} \end{matrix}] = [\begin{matrix} 1 \\ - i \\ k \\ - ik \end{matrix}] .

Solving for $v^{*}$ in this system results in a new eigenvector, with $v_{1}^{*}$ set to 1, as

v^{*} = {[\begin{matrix} 1 & \frac{k α_{s} - 2 - α_{s}}{ik α_{s}} & 1 + k & \frac{k α_{s} - 2 - α_{s}}{i α_{s}} - i \end{matrix}]}^{t} .

The solution to the system (61) can, thus, be written as

(62)

\begin{matrix} W = D_{1} W_{1} + D_{2} W_{2} = D_{1} v \exp (k η) + D_{2} [v η \exp (k η) + v^{*} \exp (k η)] . \end{matrix}

Referring to (22), the solution for the differential equation (16) can be written as the first two components of the general solution tensor $W$ , which reads

(63)

\begin{matrix} w = D_{1} [\begin{matrix} 1 \\ - i \end{matrix}] \exp (k η) + D_{2} η [\begin{matrix} 1 \\ - i \end{matrix}] \exp (k η) + D_{2} [\begin{matrix} 1 \\ \frac{k α_{s} - 2 - α_{s}}{ik α_{s}} \end{matrix}] \exp (k η), \end{matrix}

expressed in terms of the basis $t$ and $n$ as

(64)

\begin{matrix} w = [[D_{1} + D_{2} + η D_{2}] t + [- i D_{1} + \frac{[k α_{s} - 2 - α_{s}]}{i k α_{s}} D_{2} - i η D_{2}] n] \exp (k η) . \end{matrix}

Finally, the constants $D_{1}$ and $D_{2}$ can be expressed in terms of two other constants $C$ and $D$ defined as

(65)

\begin{matrix} D_{1} = \frac{i}{k} C + [\frac{k α_{s} - 2}{k α_{s}}] \frac{i}{k} D and D_{2} = - \frac{i}{k} D, \end{matrix}

(66)

\begin{matrix} C = \frac{k}{i} D_{1} + [\frac{k α_{s} - 2}{α_{s} i}] D_{2} and D = - \frac{k}{i} D_{2}, \end{matrix}

leading to an alternative expression

(67)

\begin{matrix} w = \frac{1}{k} [[i C - \frac{2}{k α_{s}} i D - i η D] t + [C + \frac{1}{k} D - η D] n] \exp (k η), \end{matrix}

which using the definition $α_{s} = [λ_{s} + μ_{s}] / μ_{s}$ eventually reads

(68)

\begin{matrix} w (η) = \frac{1}{k} [[i C - \frac{2 μ_{s}}{k [μ_{s} + λ_{s}]} i D - i η D] t + [C + \frac{1}{k} D - η D] n] \exp (k η) . \end{matrix}

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

Notes

References

Allen

. Analysis and Design of Structural Sandwich Panels. Oxford: Pergamon Press, 1969.

Biot

. Bending of an infinite beam on an elastic foundation. J Appl Math Mech 1937; 203: A-1–A-7.

Iacopi

Brongersma

Maex

. Compressive stress relaxation through buckling of a low-k polymer-thin cap layer system. Appl Phys Lett 2003; 82: 1380.

Yoo

Suh

Park

, et al. Physical self-assembly of microstructures by anisotropic buckling. Adv Mater 2002; 14: 1383–1387.

Ohzono

Shimomura

. Ordering of microwrinkle patterns by compressive strain. Phys Rev B 2004; 69: 132202.

Volynskii

Bazhenov

Lebedeva

, et al. Mechanical buckling instability of thin coatings deposited on soft polymer substrates. J Mater Sci 2000; 35: 547–554.

Huang

Suo

. Instability of a compressed elastic film on a viscous layer. Int J Solids Struct 2002; 39: 1791–1802.

Potier-Ferry

Belouettar

, et al. Multiple bifurcations in wrinkling analysis of thin films on compliant substrates. Int J Non-Linear Mech 2015; 76: 203–222.

Groenewold

. Wrinkling of plates coupled with soft elastic media. Phys A 2001; 298: 32–45.

10.

Stoop

Lagrange

Terwagne

, et al. Curvature-induced symmetry breaking determines elastic surface patterns. Nat Mater 2015; 14: 337–342.

11.

Potier-Ferry

. A multi-scale modeling framework for instabilities of film/substrate systems. J Mech Phys Solids 2016; 86: 150–172.

12.

Zhao

. Multimodal surface instabilities in curved film-substrate structures. ASME J Appl Mech 2017; 84: 081001.

13.

Holland

Feng

, et al. Instabilities of soft films on compliant substrates. J Mech Phys Solids 2017; 98: 350–365.

14.

Andres

Steinmann

Budday

. The origin of compression influences geometric instabilities in bilayers. Proc R Soc A Math Phys Eng Sci 2018; 474: 20180267.

15.

Lee

Triantafyllidis

Barber

, et al. Surface instability of an elastic half space with material properties varying with depth. J Mech Phys Solids 2008; 56: 858–868.

16.

Huang

Suo

. Very thin solid-on-liquid structures: the interplay of flexural rigidity, membrane force, and interfacial force. Thin Solid Films 2003; 429: 273–281.

17.

Huang

Stafford

Vogt

. Effect of surface properties on wrinkling of ultrathin films. J Aerospace Eng 2007; 20: 38–44.

18.

Song

Jiang

Liu

, et al. Buckling of a stiff thin film on a compliant substrate in large deformation. Int J Solids Struct 2008; 45: 3107–3121.

19.

Yin

Cao

, et al. Surface wrinkling of anisotropic films bonded on a compliant substrate. Int J Solids Struct 2018; 141–142: 219–231.

20.

Jia

Cao

Liu

, et al. Wrinkling of a bilayer resting on a soft substrate under in-plane compression. Phil Mag 2012; 92: 1554–1568.

21.

Lejeune

Javili

Linder

. Understanding geometric instabilities in thin films via a multi-layer model. Soft Matter 2016; 12: 806.

22.

Huang

Hong

Suo

. Nonlinear analyses of wrinkles in a film bonded to a compliant substrate. J Mech Phys Solids 2005; 53: 2101–2118.

23.

Audoly

Boudaoud

. Buckling of a stiff film bound to a compliant substrate - Part I: Formulation, linear stability of cylindrical patterns, secondary bifurcations. J Mech Phys Solids 2008; 56: 2401–2421.

24.

Chen

Hutchinson

. Herringbone buckling patterns of compressed thin films on compliant substrates. ASME J Appl Mech 2004; 71: 597–603.

25.

Huang

Hong

Suo

. Evolution of wrinkles in hard films on soft substrates. Phys Rev E 2004; 70: 030601.

26.

Wang

Zhao

. A three-dimensional phase diagram of growth-induced surface instabilities. Sci Rep 2015; 5: 8887.

27.

Jin

Chen

Hayward

, et al. Creases on the interface between two soft materials. Soft Matter 2014; 10: 303.

28.

Cao

Hutchinson

. From wrinkles to creases in elastomers: the instability and imperfection-sensitivity of wrinkling. Proc Royal Soc A 2012; 468: 94–115.

29.

Hong

Zhao

Suo

. Formation of creases on the surfaces of elastomers and gels. Appl Phys Lett 2009; 95: 111901.

30.

Jin

Cai

Suo

. Creases in soft tissues generated by growth. Epl 2011; 95: 64002.

31.

Jin

Auguste

Hayward

, et al. Bifurcation diagrams for the formation of wrinkles or creases in soft bilayers. J Appl Mech 2015; 82: 061008.

32.

Zang

Zhao

Cao

, et al. Localized ridge wrinkling of stiff films on compliant substrates. J Mech Phys Solids 2012; 60: 1265–1279.

33.

Potier-Ferry

Belouettar

, et al. 3D finite element modeling for instabilities in thin films on soft substrates. Int J Solids Struct 2014; 51: 3619–3632.

34.

Cao

Feng

, et al. Mechanics of morphological instabilities and surface wrinkling in soft materials: a review. Soft Matter 2012; 8: 5728.

35.

Xie

Huang

Cao

, et al. Buckling and postbuckling of stiff lamellae in a compliant matrix. Compos Sci Technol 2014; 99: 89–95.

36.

Kaynia

Rudykh

, et al. Wrinkling of interfacial layers in stratified composites. Adv Eng Mater 2013; 15: 921–926.

37.

Brau

Vandeparre

Sabbah

, et al. Multiple-length-scale elastic instability mimics parametric resonance of nonlinear oscillators. Nat Phys 2011; 7: 56–60.

38.

Brau

Damman

Diamant

, et al. Wrinkle to fold transition: influence of the substrate response. Soft Matter 2013; 9: 8177.

39.

Slesarenko

Galich

, et al. Instabilities and pattern formations in 3D-printed deformable fiber composites. Compos B Eng 2018; 148: 114–122.

40.

Colin

Holland

. Layer wrinkling in an inhomogeneous matrix. Int J Solids Struct 2019; 156-157: 119–125.

41.

Slesarenko

Rudykh

. Microscopic instabilities and elastic wave propagation in finitely deformed laminates with compressible hyperelastic phases. Eur J Mech A Solid 2019; 73: 126–136.

42.

Zhao

Cao

, et al. Buckling of an elastic fiber with finite length in a soft matrix. Soft Matter 2016; 12: 2086.

43.

. A mechanics model of microtubule buckling in living cells. J Biomech 2008; 41: 1722–1729.

44.

Gao

. Tuning the wrinkling patterns of an interfacial/coating layer via a regulation interphase. Int J Solids Struct 2017; 104-105: 92–102.

45.

Zhang

. Symplectic analysis for wrinkles: a case study of layered neo-hookean structures. ASME J Appl Mech 2017; 84: 071002.

46.

Lejeune

Javili

Linder

. An algorithmic approach to multi-layer wrinkling. Extreme Mech Lett 2016; 7: 10–17.

47.

O’Keeffe

Moulton

Waters

, et al. Growth-induced axial buckling of a slender elastic filament embedded in an isotropic elastic matrix. Int J Non-Linear Mech 2013; 56: 94–104.

48.

Rudykh

Boyce

. Transforming wave propagation in layered media via instability-induced interfacial wrinkling. Phys Rev Lett 2014; 112: 034301.

49.

Liu

Terwagne

, et al. Buckling of an elastic rod embedded on an elastomeric matrix: planar vs. non-planar configurations. Soft Matter 2014; 10: 6294.

50.

Kuhl

. Growing matter: A review of growth in living systems. J Mech Behav Biomed Mater 2014; 29: 529–543.

51.

Rodriguez

Hoger

McCulloch

. Stress-dependent finite-growth in soft elastic tissues. J Biomech 1994; 27: 455–467.

52.

Garikipati

Arruda

Grosh

, et al. A continuum treatment of growth in biological tissue: The coupling of mass transport and mechanics. J Mech Phys Solids 2004; 52: 1595–1625.

53.

Lubarda

Hoger

. On the mechanics of solids with a growing mass. Int J Solids Struct 2002; 39: 4627–4664.

54.

Skalak

Zargaryan

Jain

, et al. Compatibility and the genesis of residual stress by volumetric growth. J Math Biol 1996; 34: 889–914.

55.

Yavari

. A geometric theory of growth mechanics. J Nonlinear Sci 2010; 20: 781–830.

56.

Humphrey

. Continuum biomechanics of soft biological tissues. Proc R Soc Lond A 2003; 459: 3–46.

57.

Menzel

Kuhl

. Frontiers in growth and remodeling. Mech Res Commun 2012; 42: 1–14.

58.

Epstein

Maugin

. Thermomechanics of volumetric growth in uniform bodies. Int J Plast 2000; 16: 951–978.

59.

Ateshian

. On the theory of reactive mixtures for modeling biological growth. Biomech Model Mechanobiol 2007; 6: 423–445.

60.

Grillo

Wittum

Giaquinta

, et al. A multiscale analysis of growth and diffusion dynamics in biological materials. Int J Eng Sci 2009; 47: 261–283.

61.

Humphrey

Rajagopal

. A constrained mixture model for growth and remodeling of soft tissues. Math Models Methods App Sci 2002; 12: 407–430.

62.

Garikipati

. The kinematics of biological growth. Appl Mech Rev 2009; 62: 030801.

63.

Ciarletta

Preziosi

Maugin

. Mechanobiology of interfacial growth. J Mech Phys Solids 2013; 61: 852–872.

64.

Moulton

Lessinnes

Goriely

. Morphoelastic rods. Part I: A single growing elastic rod. J Mech Phys Solids 2013; 61(2): 398–427.

65.

Ambrosi

Ateshian

Arruda

, et al. Perspectives on biological growth and remodeling. J Mech Phys Solids 2011; 59: 863–883.

66.

Hosseini

Taber

. How mechanical forces shape the developing eye. Prog Biophys Mol Biol 2018; 137: 25–36.

67.

Hosseini

Beebe

Taber

. Mechanical effects of the surface ectoderm on optic vesicle morphogenesis in the chick embryo. J Biomech 2014; 47: 3837–3846.

68.

Araujo

McElwain

DLS

. A history of the study of solid tumour growth: The contribution of mathematical modelling. Bull Math Biol 2004; 66: 1039–1091.

69.

Hosseini

Garcia

Taber

. A new hypothesis for foregut and heart tube formation based on differential growth and actomyosin contraction. Development 2017; 144: 2381–2391.

70.

Genzer

Groenewold

. Soft matter with hard skin: From skin wrinkles to templating and material characterization. Soft Matter 2006; 2: 310.

71.

Efimenko

Rackaitis

Manias

, et al. Nested self-similar wrinkling patterns in skins. Nat Mater 2005; 4: 293–297.

72.

Hohlfeld

Mahadevan

. Unfolding the sulcus. Phys Rev Lett 2011; 106: 105702.

73.

Budday

Steinmann

Kuhl

. The role of mechanics during brain development. J Mech Phys Solids 2014; 72: 75–92.

74.

Verner

Garikipati

. A computational study of the mechanisms of growth-driven folding patterns on shells, with application to the developing brain. Extreme Mech Lett 2018; 18: 58–69.

75.

Budday

Steinmann

. On the influence of inhomogeneous stiffness and growth on mechanical instabilities in the developing brain. Int J Solids Struct 2018; 132-133: 31–41.

76.

Ciarletta

Balbi

Kuhl

. Pattern selection in growing tubular tissues. Phys Rev Lett 2014; 113: 248101.

77.

Cao

Feng

, et al. Surface wrinkling of mucosa induced by volumetric growth: Theory, simulation and experiment. J Mech Phys Solids 2011; 59: 758–774.

78.

Wiggs

Hrousis

Drazen

, et al. On the mechanism of mucosal folding in normal and asthmatic airways. J Appl Physiol 1997; 83: 1814–1821.

79.

Ciarletta

Destrade

Gower

, et al. Morphology of residually stressed tubular tissues: Beyond the elastic multiplicative decomposition. J Mech Phys Solids 2016; 90: 242–253.

80.

Moulton

Goriely

. Circumferential buckling instability of a growing cylindrical tube. J Mech Phys Solids 2011; 59: 525–537.

81.

Eskandari

Javili

Kuhl

. Elastosis during airway wall remodeling explains multiple co-existing instability patterns. J Theor Biol 2016; 403: 209–218.

82.

Limbert

Kuhl

. On skin microrelief and the emergence of expression micro-wrinkles. Soft Matter 2018; 14: 1292.

83.

Buganza Tepole

Joseph Ploch

Wong

, et al. Growing skin: A computational model for skin expansion in reconstructive surgery. J Mech Phys Solids 2011; 59: 2177–2190.

84.

Bowden

Brittain

Evans

, et al. Spontaneous formation of ordered structures in thin films of metals supported on an elastomeric polymer. Nature 1998; 393: 146–149.

85.

Huck

WTS

Bowden

Onck

, et al. Ordering of spontaneously formed buckles on planar surfaces. Langmuir 2000; 16: 3497–3501.

86.

Sun

Xia

Moon

, et al. Folding wrinkles of a thin stiff layer on a soft substrate. Proc R Soc A 2012; 468: 932–953.

87.

Hutchinson

. The role of nonlinear substrate elasticity in the wrinkling of thin films. Phil Trans R Soc A 2013; 371: 20120422.

88.

Watanabe

Shirai

Hirai

. Wrinkled polypyrrole electrode for electroactive polymer actuators. J Appl Phys 2002; 92: 4631.

89.

Lacour

Wagner

Huang

, et al. Stretchable gold conductors on elastomeric substrates. Appl Phys Lett 2003; 82: 2404.

90.

Lacour

Jones

Suo

, et al. Design and performance of thin metal film interconnects for skin-like electronic circuits. IEEE Electron Device Lett 2004; 25: 179–181.

91.

Conroy Broderick

Chen

, et al. Wrinkles in soft dielectric plates. J Mech Phys Solids 2018; 119: 298–318.

92.

Seifi

Wang

Park

. Surface tension effects on surface instabilities of dielectric elastomers Preprint 2016; arXiv:1611.06419.

93.

Chung

Nolte

Stafford

. Surface wrinkling: A versatile platform for measuring thin-film properties. Adv Mater 2011; 23: 349–368.

94.

DiCarlo

Quiligotti

. Growth and balance. Mech Res Commun 2002; 29: 449–456.

95.

Javili

Steinmann

Kuhl

. A novel strategy to identify the critical conditions for growth-induced instabilities. J Mech Behav Biomed Mater 2014; 29: 20–32.

96.

Lü

Chen

, et al. Modified multiplicative decomposition model for tissue growth: Beyond the initial stress-free state. J Mech Phys Solids 2018; 118: 133–151.

97.

Taber

. Biomechanics of growth, remodeling, and morphogenesis. Appl Mech Rev 1995; 48: 487.

98.

Dortdivanlioglu

Javili

Linder

. Computational aspects of morphological instabilities using isogeometric analysis. Comput Meth Appl Mech Eng 2017; 316: 261–279.

99.

Cao

Hutchinson

. Wrinkling phenomena in neo-Hookean film/substrate bilayers. ASME J Appl Mech 2012; 79: 031019.

100.

Goriely

Ben Amar

. Differential growth and instability in elastic shells. Phys Rev Lett 2005; 94: 198103.

101.

Ben Amar

Goriely

. Growth and instability in elastic tissues. J Mech Phys Solids 2005; 53: 2284–2319.

102.

Javili

Dortdivanlioglu

Kuhl

, et al. Computational aspects of growth-induced instabilities through eigenvalue analysis. Comput Mech 2015; 56: 405–420.

103.

Kuhl

Menzel

Steinmann

. Computational modeling of growth. Comput Mech 2003; 32: 71–88.