Post-buckling of an elastic half-space coated by double layers

Abstract

We investigate the buckling and post-buckling properties of a hyperelastic half-space coated by two hyperelastic layers when the composite structure is subjected to a uniaxial compression. In the case of a half-space coated with a single layer, it is known that when the shear modulus ${\mu }_f$ of the layer is larger than the shear modulus ${\mu }_s$ of the half-space, a linear analysis predicts the existence of a critical stretch and wave number, whereas a weakly nonlinear analysis predicts the existence of a threshold value of the modulus ratio ${\mu }_s/{\mu }_f\approx 0.57$ below which the buckling is super-critical and above which the buckling is sub-critical. It is shown that when another layer is added, a larger variety of behaviour can be observed. For instance, buckling can occur at a preferred wavenumber super-critically even if both layers are softer than the half-space although the top layer would need to be harder than the bottom layer. When the shear modulus of the bottom layer lies in a certain interval, the super-critical to sub-critical transition can happen a number of times as the shear modulus of the top layer is increased gradually. Thus, an extra layer imparts more flexibility in producing wrinkling patterns with desired properties, and our weakly nonlinear analysis provides a road map on the parameter regimes where this can be achieved.

1. Introduction

Buckling of a coated elastic half-space induced by uniaxial compression is a topic that has been much studied in recent decades. Motivated mostly by the desire to suppress buckling as a precursor to structural failure, early studies under the framework of nonlinear elasticity include the linear analyses by Dorris and Nemat-Nasser [1], Shield et al. [2], Ogden and Sotiropoulos [3], Bigoni et al. [4], Steigmann and Ogden [5] and the nonlinear post-buckling analysis by Cai and Fu [6]. The subsequent experimental work by Bowden et al. [7,8] on pattern formation at the micrometre and sub-micrometre scales opened the possibility to use such pattern formation to achieve a variety of useful purposes. Known applications now range from cell patterning [9], optical gratings [10–12], creation of surfaces with desired wetting and adhesion [13, 14] and buckling-based metrology [15]. Pattern formation is also known to play an important role in many biological processes where the main driving mechanism is growth [16]. Motivated by these newly found applications, further linear and nonlinear analyses have been conducted by Cai and Fu [17, 18], Steigmann and Ogden [19], Sun et al. [20], Hutchinson [21], Ciarletta [22], Fu and Ciarletta [23], Ciarletta and Fu [24], Zhang [25], Holland et al. [26], Alawiye et al. [27, 28] and Rambausek and Danas [29]. There also exists a large body of literature that employs approximate plate theories and simplified interfacial conditions in order to provide qualitative understanding of the phenomena observed; see, e.g., Chen and Hutchinson [30], Huang et al. [31], Audoly and Boudaoud [32] and Song et al. [33]. We refer to the reviews by Genzer and Groenewold [34], Yang et al. [35], Li et al. [36] and Dimmock [37] for a comprehensive review of the relevant literature.

A major result arising from a weakly nonlinear analysis [6, 21, 28] is that there exists a critical stiffness ratio that marks the transition from sub-critical bifurcation to super-critical bifurcation. It can then be expected, and indeed confirmed by many recent numerical and experimental studies, that in the sub-critical regime localization is the norm [38] whereas in the super-critical regime stable periodic patterns and further secondary bifurcations are the norm [18, 39–47].

In this article, we extend the analysis of Cai and Fu [6] and consider the effects of adding an extra layer to the film–substrate bilayer structure. The structure of a substrate coated by two layers has been recognized to offer enhanced capabilities in buckling-based metrology [48–50], stretchable electronics [51] and surface pattern switching [49, 52]. Multi-layering is of course the norm in biological systems (e.g. the human skin), whereas in other situations adding extra layers arises out of necessity. For instance, in cell patterning, there may exist the conflicting demands that, on the one hand, well-ordered wrinkling patterns are desired and, on the other hand, the top film is required to be sufficiently soft to dictate a certain cell behaviour [37].

For the structure under consideration, there exist two modulus ratios, namely $r_1={\mu }_s/{\mu }_1$ and $r_2=u_s/{\mu }_2$ , where ${\mu }_s$ , ${\mu }_1$ and ${\mu }_2$ are the shear moduli of the substrate, the first layer (i.e. the lower layer) and the top layer, respectively. In this case, the sign of the nonlinear coefficient $c_1$ in the amplitude equation, that determines the sub-critical to super-critical transition, is a function of $r_1$ and $r_2$ . Our main objective is to determine the dividing curve in the $(r_1,r_2)$ -plane where this coefficient vanishes. It is hoped that such information can help guide the design process when robust wrinkles are desired. We present results for three representative cases when the thickness ratio of the two layers is equal to 1, $0.1$ and 10, respectively, but the methodology is valid for any other thickness ratio (and for any material model).

The rest of this paper is divided into five sections as follows. After formulating the buckling problem in the next section, we present in Section 3 the necessary linear analysis that produces the bifurcation condition. The bifurcation condition exhibits two interesting features. First, even if both layers are softer than the half-space, the stretch can still exhibit a maximum at a finite wavenumber. Second, there is a range of thickness ratios for which multiple stretch maxima exist and mode switching becomes possible as a material parameter is varied. This is in contrast with the situation of a single layer where a stretch maximum is only possible when the layer is stiffer than the half-space and no mode switching is possible. In Section 4 we conduct a weakly nonlinear analysis in order to determine the nature of buckling, namely whether it is super-critical or sub-critical. A richer variety of behaviour than what is possible for a single layer is uncovered. The paper is concluded in Section 5 with a summary, and a discussion of limitations of the current study, and possible future work.

2. Governing equations

We first formulate the governing equations for a general homogeneous elastic body B composed of a non-heat-conducting incompressible elastic material. Three configurations of B are involved in our analysis: the initial unstressed configuration $B_0$ , the finitely stressed equilibrium configuration denoted by $B_{\mathrm{e}}$ , and the current configuration $B_{\mathrm{t}}$ that is obtained from $B_{\mathrm{e}}$ by a small perturbation. The position vectors of a representative material particle relative to a common Cartesian coordinate system are denoted by $\mathit{X}$ , $\mathit{x}$ and $\tilde{\mathit{x}}$ in $B_0$ , $B_{\mathrm{e}}$ and $B_{\mathrm{t}}$ , respectively. The associated coordinates are written as $X_A$ , $x_i$ and ${\tilde{x}}_i$ . We write

$$\tilde{\mathit{x}}=\mathit{x}(\mathit{X})+\mathit{u}(\mathit{x}),$$

(1)

where $\mathit{u}(\mathit{x})$ is the incremental displacement associated with the deformation $B_{\mathrm{e}}\to B_{\mathrm{t}}$ .

The deformation gradients arising from the deformations $B_0\to B_{\mathrm{t}}$ and $B_0\to B_{\mathrm{e}}$ are denoted by $\tilde{\mathit{F}}$ and $\overline{\mathit{F}}$ , respectively, and are defined by their Cartesian components

$${\tilde{F}}_{\mathrm{iA}}=\frac{\partial {\tilde{x}}_i}{\partial X_A},{\overline{F}}_{\mathrm{iA}}=\frac{\partial x_i}{\partial X_A}.$$

(2)

It then follows that

$${\tilde{F}}_{\mathrm{iA}}=({\delta }_{\mathrm{ij}}+u_{i,j}){\overline{F}}_{\mathrm{jA}},$$

(3)

where here and henceforth a comma signifies differentiation with respect to the indicated coordinate, with ‘ $,A$ ’ and ‘ $,j$ ’ meaning differentiation with respect to $X_A$ and $x_j$ , respectively.

In the absence of body forces, the equations of equilibrium and the incompressibility constraint are given by

$${\pi }_{\mathrm{iA},A}=0,\det \mathit{F}=1,$$

(4)

where $\mathit{\pi }$ is the first Piola–Kirchhoff stress which, in component form, is given by

$${\pi }_{\mathrm{iA}}=\frac{\partial W}{\partial F_{\mathrm{iA}}}-{\mathrm{pF}}_{\mathrm{Ai}}^{-1},$$

(5)

with W denoting the strain-energy function (per unit volume in the reference configuration) and p a Lagrange multiplier enforcing the incompressibility constraint. We denote by $\overline{p}$ and $\overline{p}+\tilde{p}$ the values of p associated with the deformations $B_0\to B_{\mathrm{e}}$ and $B_0\to B_{\mathrm{t}}$ , respectively. Thus, $\tilde{p}$ is the increment of p which is an additional field induced by the constraint of incompressibility.

We define the incremental stress tensor $\mathbf{\chi }$ through

$$\mathbf{\chi }={\overline{J}}^{-1}(\tilde{\mathit{\pi }}-\overline{\mathit{\pi }}){\overline{F}}^{\mathrm{T}},$$

(6)

where $\overline{J}=\det \overline{F}$ and $\overline{\mathit{\pi }}$ and $\tilde{\mathit{\pi }}$ are the first Piola–Kirchhoff stresses in $B_{\mathrm{e}}$ and $B_{\mathrm{t}}$ , respectively. It can be shown, see e.g. [53], that $\mathbf{\chi }$ satisfies the equilibrium equation

$${\chi }_{\mathrm{ij},j}=0,$$

(7)

and has a Taylor expansion given by

$$\begin{array}{c}{\chi }_{\mathrm{ij}}={\mathcal{A}}_{\mathrm{jilk}}^1{u}_{k,l}+\frac{1}{2}{\mathcal{A}}_{\mathrm{jilknm}}^2{u}_{k,l}{u}_{m,n}+\frac{1}{6}{\mathcal{A}}_{\mathrm{jilknmqp}}^3{u}_{k,l}{u}_{m,n}{u}_{p,q}\hfill \\ +\overline{p}(u_{j,i}-u_{j,k}{u}_{k,i}+u_{j,k}{u}_{k,l}{u}_{l,i})-\tilde{p}({\delta }_{\mathrm{ji}}-u_{j,i}+u_{j,k}{u}_{k,i})+O({ϵ}^4),\hfill \end{array}$$

(8)

where $ϵ$ is a small parameter characterizing the magnitude of $u_{i,j}$ and $\tilde{p}$ , and the tensors ${\mathcal{A}}^1$ , ${\mathcal{A}}^2$ and ${\mathcal{A}}^3$ are the first-, second- and third-order tensors of instantaneous elastic moduli in $B_{\mathrm{e}}$ whose expressions are typified by [54]

$${\mathcal{A}}_{\mathrm{jilk}}^1={\overline{J}}^{-1}{\overline{F}}_{\mathrm{jA}}{\overline{F}}_{\mathrm{lB}}\frac{{\partial }^{2}W}{\partial F_{\mathrm{iA}}\partial F_{\mathrm{kB}}}{|}_{\mathit{F}=\overline{\mathit{F}}}.$$

(9)

On substituting (8) into (7) and simplifying with the use of the identity $({\mathrm{JF}}_{\mathrm{Ai}}^{-1}{)}_{,A}\equiv 0$ , we obtain

$$\begin{array}{c}{\mathcal{A}}_{\mathrm{jilk}}^1{u}_{k,\mathrm{lj}}+{\mathcal{A}}_{\mathrm{jilknm}}^2{u}_{m,n}{u}_{k,\mathrm{lj}}+\frac{1}{2}{\mathcal{A}}_{\mathrm{jilknmqp}}^3{u}_{m,n}{u}_{p,q}{u}_{k,\mathrm{lj}}\hfill \\ -{\tilde{p}}_{,j}({\delta }_{\mathrm{ji}}-u_{j,i}+u_{j,m}{u}_{m,i})+O({ϵ}^4)=0.\hfill \end{array}$$

(10)

These incremental equilibrium equations are supplemented by the incompressibility condition in the form

$$u_{i,i}=\frac{1}{2}{u}_{m,n}{u}_{n,m}-\frac{1}{2}(u_{i,i}{)}^2-\det (u_{m,n}).$$

(11)

We now specialize the above equations to the structure of an incompressible elastic half-space coated by two incompressible elastic layers. Each component (layer or half-space) in this structure is a homogeneous elastic body to which the above equations apply. We choose our common coordinate system such that the half-space, the first layer and second layer are defined by $-\infty <x_2\le 0$ , $0\le x_2\le h_1$ and $h_1\le x_2\le h_1+h_2$ , respectively, where $h_1$ and $h_2$ are constant thicknesses of the two layers to be specified.

To simplify analysis, we assume that the bonded structure is in a state of plane strain so that $u_3=0$ and $u_1$ and $u_2$ are independent of $x_3$ . We also assume that the principal axes of stretch corresponding to the finite deformation are aligned with the coordinate axes.

Our problem is then to solve (11) and (10) in $-\infty <x_2\le h_1+h_2$ subject to the following auxiliary conditions:

$(\mathrm{i})$ traction-free boundary conditions,

$${\chi }_{i2}=0,\mathrm{on}{x}_2=h_1+h_2,$$

(12)

$(\mathrm{ii})$ interfacial continuity conditions,

$$[u_i]=0,[{\chi }_{i2}]=0,\mathrm{at}{x}_2=0,\mathrm{or}{h}_1,$$

(13)

$(\mathrm{iii})$ decay conditions,

$$u_i\to 0\mathrm{as}{x}_2\to -\infty ,$$

(14)

where the notation $[f]$ denotes the jump of f when an interface is crossed. The trivial solution $u_i=0$ , $\tilde{p}=0$ is clearly one solution. Our aim is to find the conditions under which the above boundary value problem has a non-trivial solution.

Our numerical calculations will be carried out for the case when the layers and half-space are all composed of either neo-Hookean or Gent materials, and the prestress in $B_{\mathrm{e}}$ takes the form of a uniaxial compression along the $x_1$ -direction with stretch $\lambda$ . For neo-Hookean and Gent materials, the strain energy function is given by

$$W=\frac{1}{2}\mu (I_1-3)\mathrm{and}W=-\frac{1}{2}\mu J_m\ln (1-\frac{I_1-3}{J_m}),$$

(15)

respectively, where $I_1=\mathrm{tr}\mathit{B}$ , $\mathit{B}=\mathit{F}{\mathit{F}}^{\mathrm{T}}$ , $\mu$ is the shear modulus and $J_m$ is a constitutive constant characterizing material extensibility with the limit $J_m\to \infty$ recovering the neo-Hookean model. We denote the shear moduli of the substrate and the two layers by ${\mu }_s$ , ${\mu }_1$ and ${\mu }_2$ , respectively, and define the dimensionless parameters $r_1={\mu }_s/{\mu }_1$ , $r_2={\mu }_s/{\mu }_2$ . When the Gent material model is used, we shall assume that $J_m$ takes the same value in all the three components, but our methodology can accommodate different values of $J_m$ and is, in fact, valid for any strain-energy function. All of our symbolic manipulations and numerical integrations are carried out with the aid of Mathematica [55].

3. Linear theory

In preparation for the later nonlinear analysis, we first consider the linearized version of the problem specified by (10), (11) and (12) – (14). The linearized governing equations are

$$u_{i,i}=0,{\mathcal{A}}_{\mathrm{jilk}}^1{u}_{k,\mathrm{lj}}-{\tilde{p}}_{,i}=0,-\infty <x_2\le h_1+h_2,$$

(16)

and the auxiliary conditions are

$$T_i^{(l)}=0\mathrm{on}{x}_2=h_1+h_2,$$

(17)

$$[u_i]=0,[T_i^{(l)}]=0,\mathrm{on}{x}_2=0\mathrm{or}{h}_1,$$

(18)

$$u_i\to 0\mathrm{as}{x}_2\to -\infty ,$$

(19)

where the linearized traction $T_i^{(l)}$ is given by

$$T_i^{(l)}={\mathcal{A}}_{2\mathrm{ilk}}^1{u}_{k,l}+\overline{p}{u}_{2,i}-\tilde{p}{\delta }_{2i}.$$

Equation (16)₁ implies the existence of a ‘stream’ function $\psi (x_1,x_2)$ such that

$$u_1={\psi }_{,2},u_2=-{\psi }_{,1}.$$

(20)

Substituting (20) into (16)₂ and eliminating $\tilde{p}$ through cross differentiation, we obtain [56]

$$\gamma {\psi }_{,2222}+2\beta {\psi }_{,1122}+\alpha {\psi }_{,1111}=0,$$

(21)

where

$$\alpha ={\mathcal{A}}_{1212}^1,2\beta ={\mathcal{A}}_{1111}^1+{\mathcal{A}}_{2222}^1-2{\mathcal{A}}_{1122}^1-2{\mathcal{A}}_{1221}^1,\gamma ={\mathcal{A}}_{2121}^1.$$

(22)

We look for a periodic wrinkling solution in the form

$$\psi =H(k{x}_2){\mathrm{e}}^{\mathrm{i}k{x}_1},$$

(23)

where k is the wave number. On inserting (23) into (21) and solving the resulting fourth-order ordinary differential equation for H, we obtain

$$H(k{x}_2)=\{\begin{array}{cc}{A}_1\exp (k{s}_1{x}_2)+A_2\exp (k{s}_2{x}_2),\hfill & x_2\in (-\infty ,0),\hfill \\ \multicolumn{1}{c}{\sum _{j=1}^4{\hat{A}}_j\exp (k{s}_j{x}_2),}& x_2\in (0,h_1),\hfill \\ \multicolumn{1}{c}{\sum _{j=1}^4{\tilde{A}}_j\exp (k{s}_j{x}_2),}& x_2\in (h_1,h_1+h_2),\hfill \end{array}$$

(24)

where $A_1,A_2,\dots ,{\tilde{A}}_3,{\tilde{A}}_4$ are disposable constants, and $s_1$ and $s_2$ are the two roots of

$$\gamma s^4-2\beta s^2+\alpha =0,$$

(25)

that have positive real parts (so that $H\to 0$ as $x_2\to -\infty$ ), and $s_3=-s_1$ , $s_4=-s_2$ . In the above expressions, we have assumed that $s_1$ and $s_2$ are independent of the elastic moduli so that they take the same values in the layers and half-space. For the neo-Hookean material model, we have $s_1=1$ , $s_2={\lambda }^2$ .

On substituting the above solutions into the auxiliary conditions (17) – (19), we obtain a matrix equation of the form $K\mathit{c}=\mathit{0}$ , where K is a $10\times 10$ matrix and $\mathit{c}$ is the column vector formed from the 10 unknowns $A_1,A_2,\dots ,{\tilde{A}}_3,{\tilde{A}}_4$ . The existence of a non-trivial solution requires

$$\det K=0,$$

(26)

which is the bifurcation condition relating the pre-stretch $\lambda$ to the wavenumber k. The wavenumber appears in the bifurcation condition through $k{h}_1$ and $k{h}_2$ . Thus, in terms of the thickness ratio $h=h_2/h_1$ , the bifurcation condition (26) takes the form $F(k{h}_1,h,\lambda ,r_1,r_2)=0$ .

For each choice of h, $r_1$ and $r_2$ , the above bifurcation condition define $\lambda$ as a function of $k{h}_1$ . It has two important limits, $k\to 0$ or $k\to \infty$ , for which the corresponding value of $\lambda$ is well-known when the material is neo-Hookean. When $k\to 0$ , which is equivalent to $h_1\to 0$ and $h_2\to 0$ , the effect of the two layers becomes negligible and the critical value of $\lambda$ is that associated with a half-space, that is $0.54$ as first obtained by Biot [57]. When $k\to \infty$ , the wavelength of the modes tends to zero and the critical modes will localize either near the traction-free surface or at one of the two interfaces (between the two layers or between the first layer and the half-space). The mode localized near the traction-free surface is again a surface wave mode with critical value of $\lambda$ equal to $0.54$ , and this mode corresponds to the same branch/curve that tends to $0.54$ when $k\to 0$ . The interfacial modes have previously been studied by Dowaikh and Ogden [56] and Fu [58]. It was shown in [58] that the surface impedance tensorM for a pre-stressed half-space occupying the region $x_2>0$ , with principal axes of stretched aligned with the coordinate axes, has the explicit expression

$$M=\left(\begin{array}{cc}{M}_1\hfill & \mathrm{i}{M}_4\hfill \\ \multicolumn{1}{c}{-\mathrm{i}{M}_4}& M_2\hfill \end{array}\right),$$

(27)

with

$$\begin{array}{c}{M}_1=\sqrt{T_{11}(Q_{11}+T_{22}+2\sqrt{Q_{22}{T}_{11}}-2R_{12}-2R_{21})},\hfill \\ M_2=\frac{M_1^3}{2T_{11}^2}-\frac{M_1}{2T_{11}}(Q_{11}+T_{22}-2R_{12}-2R_{21}),\hfill \\ M_4=\frac{M_1^2}{2T_{11}}-\frac{1}{2}(Q_{11}+T_{22}-2R_{12}+2\overline{p}),\hfill \end{array}$$

and

$$T_{\mathrm{ij}}={\mathcal{A}}_{2i2j}^1,R_{\mathrm{ij}}={\mathcal{A}}_{1i2j}^1,Q_{\mathrm{ij}}={\mathcal{A}}_{1i1j}^1.$$

The bifurcation condition for a single half-space may be written as $\det M=0$ that reproduces the well-known equation ${\lambda }^6+{\lambda }^4+3{\lambda }^2-1=0$ for the neo-Hookean material. This condition yields the Biot value $\lambda =0.54$ . The bifurcation condition for an interfacial mode localized at the interface of two half-spaces is given by $\det (M+{\hat{M}}^{*})=0$ , where M has the same meaning as above, $\hat{M}$ can be computed using (27) but with material constants replaced by those for the other half-space, and the superscript ‘*’ signifies complex conjugation. With the ratio of the two shear moduli denoted by r such that $0\le r\le 1$ (if r is larger than 1 we use its inverse instead because the two half-spaces can be exchanged), the bifurcation condition $\det (M+{\hat{M}}^{*})=0$ for the case of uniaxial compression reduces to

$$(r^2+1)({\lambda }^6+{\lambda }^4+3{\lambda }^2-1)+2r({\lambda }^6+3{\lambda }^4-2{\lambda }^2+1)=0;$$

(28)

see also [56]. It is known that a necessary condition for $\det (M+{\hat{M}}^{*})=0$ to be satisfied is that either $\det M<0$ or $\det {\hat{M}}^{*}<0$ , that is, the interfacial mode occurs at a value of $\lambda$ smaller than the Biot value $0.54$ . The solution of (28) for $\lambda$ as a function of r is a monotonically decreasing function of r, with maximum and minimum given by $0.54$ and 0, respectively. It is found that only one interfacial mode can exist as $k\to \infty$ and this mode is always localized at the interface with greater contrast in stiffness. As this mode lies below the curve that tends to $0.54$ when $k\to 0$ or $\infty$ , it is of little interest in the current context and, hence, will not be displayed. However, the two limits discussed above are used as useful checks on our numerical results for intermediate values of k.

The bifurcation condition discussed above is next used to generate plots in the $(r_1,r_2)$ -plane showing domains where a stretch maximum exists and where mode switching takes place. By mode switching we refer to situations where at a particular set of material parameters two equal maxima of $\lambda$ occur at two different values of k (corresponding to two modes with short and long wavelengths, respectively), and a small perturbation of such parameter values can make one or the other the only preferred mode. In other words, with a slight change of the material parameters, the bifurcation mode can switch from a short mode to a long mode or vice versa [49]. We shall present illustrative results for three representative cases, and for the neo-Hookean material model. Before presenting results for each case, however, we first observe some general features shared by all cases. These general features correspond to the three limits, $r_1\to 1$ , $r_2\to r_1$ or $r_2\to \infty$ , under which the structure under consideration reduces to a half-space coated by a single layer. As for the latter reduced case a stretch maximum exists only when the single layer is stiffer than the substrate, the boundary of domain where a stretch maximum exists for the current problem must contain the point $(r_1,r_2)=(1,1)$ and must approach the asymptote $r_1=1$ as $r_2\to \infty$ . It may also be deduced that the semi-infinite straight line $r_2=r_1$ with $r_1>1$ must lie in the domain of non-existence of a stretch maximum. It turns out that this line is actually part of the boundary of this domain.

In Figures 1(a), 2(a) and 3(a), we have shown the domain in the $(r_1,r_2)$ -plane where the bifurcation condition gives a stretch maximum for the three representative cases $h=1,0.1$ or 10, respectively (unshaded region). These numerical results are all consistent with the general observation made above. In addition, with regard to Figure 1(a), we note that for each fixed $r_1$ in the interval $(a_1,1)$ , as $r_2$ is increased from zero, a stretch maximum first exists, then disappears, and finally emerges again at a large value of $r_2$ . This behaviour is displayed in Figure 4 for $r_1=0.9$ . Similar behaviour is observed for fixed $r_2>1$ and variable $r_1$ . Another important result seen in Figure 1(a) is that a stretch maximum exists for all values of $r_2$ if $r_1<a_1$ , or for all values of $r_1$ if $r_2<1$ . The former scenario corresponds to the fact that provided the lower layer is sufficiently harder than the substrate (more precisely ${\mu }_1>a_1^{-1}{\mu }_s$ ), a stretch maximum exists no matter how soft the top layer is. Also, when both layers are softer than the half-space, a stretch maximum is still possible provided the first layer is softer than the second layer (i.e. $r_1>r_2$ ).

Figure 1.

Linear and nonlinear bifurcation behaviour for $h=1$ . (a) Parameter domain (unshaded) where the bifurcation stretch has a maximum. The left-most point of the bounding curve has coordinates $(a_1,2.69)$ where $a_1=0.89$ . (b) Dependence of the sign of $c_1$ on $r_1$ and $r_2$ with $G_1$ having coordinates $(0.57,0.57)$ and the blue line approaching $r_1=0.57$ as $r_2\to \infty$ .

Figure 2.

Results for $h=0.1$ . (a) Parameter domain (unshaded) where the bifurcation stretch has a maximum. The red line is where mode switching takes place and it terminates at $(a_2,0.008)$ where $a_2=0.16$ . (b) Dependence of the sign of $c_1$ on $r_1$ and $r_2$ with $G_2$ having coordinates $(0.57,0.57)$ and the blue line approaching $r_1=0.57$ as $r_2\to \infty$ . The dashed red curve is where mode switching takes place, and the green line coinciding with the red dashed line is where the sign change of $c_1$ is due to mode jumping, instead of via $c_1=0$ (i.e. $c_1\ne 0$ on the green line except at the two ends). The values of $r_2$ at the ends of the dotted vertical lines starting at $f_2=0.46$ and $g_2=0.59$ are $0.19$ and $0.34$ , respectively.

Figure 3.

Results for $h=10$ . (a) Parameter domain (unshaded) where the bifurcation stretch has a maximum. The red curves are where mode switches from short modes to long modes. As $r_2\to \infty$ , the red and black curves intersect at $r_1=0.75$ approximately, after which the red curve coincides with the black curve (so no mode switching takes place for $r_1$ large enough). The values of $r_2$ at the ends of the dotted vertical lines starting at $a_3=0.004$ , $b_3=0.42$ and $d_3=0.56$ are $0.24$ , $1.15$ and $1.21$ , respectively. (b) Dependence of the sign of $c_1$ on $r_1$ and $r_2$ with $G_3$ having coordinates $(0.57,0.57)$ and the blue line approaching $r_1=0.57$ as $r_2\to \infty$ . The dashed red curves are where mode switching takes place, and the green line coinciding with the red dashed lines is where the sign change of $c_1$ is due to mode jumping, instead of via $c_1=0$ (i.e. $c_1\ne 0$ on the green line except at the two ends). The values of $r_2$ at the ends of the dotted vertical lines starting at $f_3=0.28$ and $g_3=0.53$ are $0.84$ and $1.75$ , respectively.

Figure 4.

Evolution of the bifurcation curve when $h=1,r_1=0.9$ showing the fact that as $r_2$ is increased, the stretch maximum first disappears and then re-merges at a larger value of $r_2$ .

As h is varied around 1, it is found that multiple stretch maxima may occur when $|h-1|$ is sufficiently large. Figure 2(a) shows the counterpart of Figure 1(a) when $h=0.1$ (the lower layer is now 10 times as thick as the upper layer). The left boundary of the shaded domain is not a straight line although it looks that way. It actually has a similar shape to its counterpart in Figure 1(a) but the $r_1$ in the current case varies in a smaller interval, namely between $0.998$ and 1. The red line is where mode switching takes place (short modes are preferred below the line and long modes are preferred above the line). A typical example of mode switching is displayed in Figure 5 for $r_1=0.6$ .

Figure 5.

Evolution of the bifurcation curve when $h=0.1,r_1=0.6$ showing the fact that as $r_2$ is increased, mode switching may take place once.

When h is sufficiently larger than 1 (the lower layer becoming thinner than the top layer), we observe the novel phenomenon that mode switching may take place more than once when $r_2$ is varied whereas $r_1$ is fixed in a finite interval. This is due to the fact that the bifurcation curve now may have three maxima. Figure 3(a) shows the counterpart of Figure 2(a) when $h=10$ . The red curve is again where mode switching takes place, but now it has two branches. Thus, if $r_1$ is fixed in the interval $(b_3,d_3)$ and $r_2$ is increased gradually from 0, a long mode is preferred first. This mode gives way to a mode with intermediate wavelength as the short red line is crossed. Finally, this intermediate mode jumps to a short mode as the upper part of the long red line is crossed. An example of this mode switching is displayed in Figure 6 for $r_1=0.43$ .

Figure 6.

Evolution of the bifurcation curve when $h=10,r_1=0.43$ showing the fact that as $r_2$ is increased, mode switching may take place among three modes.

It is clear from the above representative results that the addition of a second layer may provide more flexibility in generating stable periodic patterns. Of course, the existence of a stretch maximum is only one of the necessary conditions for the existence of stable periodic patterns. Another necessary condition is that the associated bifurcations need to be super-critical. This will be discussed in the next section by referring to Figures 1(b), 2(b) and 3(b).

4. Post-buckling analysis

In this section, we shall first derive the amplitude equation for a single near-critical mode for general prestress and general material models. We then present numerical results for uniaxial compression and for neo-Hookean or Gent materials.

4.1. General prestress and material models

For general prestress, we assume that the finitely stressed state $B_{\mathrm{e}}$ is determined by a single parameter $\lambda$ , which, for instance, in the case of a uniaxial compression is the stretch in the $x_1$ -direction. We shall only consider the case when the bifurcation curve has a single maximum for $\lambda$ , which we denote by ${\lambda }_0$ . For an analysis of multiple mode interaction, we refer to Fu [59] and Fu and Ogden [60]. We denote by $B_{\mathrm{cr}}$ the critical configuration where the stretch in the $x_1$ -direction is equal to ${\lambda }_0$ . Guided by the analysis in [6], we assume that

$$\lambda ={\lambda }_0+{ϵ}^2{\lambda }_1,\overline{p}={\overline{p}}_0+{ϵ}^2{\lambda }_1{\overline{p}}_1,$$

(29)

where ${\lambda }_1$ is a constant, and ${\overline{p}}_0$ and ${\overline{p}}_1$ can be expressed in terms of ${\lambda }_0$ (and also the elastic moduli). We denote the uniformly deformed configuration associated with (29) by $B_{\mathrm{e}}$ . We observe that the expansions (3) – (11) are valid for any finitely deformed state, e.g. the critical state $({\lambda }_0,{\overline{p}}_0)$ or the perturbed state $(\lambda ,\overline{p})$ . As an alternative to the approach adopted in [6] where the independent variables $x_1$ and $x_2$ are defined in $B_{\mathrm{e}}$ , in the current post-buckling analysis these variables are defined in $B_{\mathrm{cr}}$ . It is then convenient to identify $\overline{F}$ and $\overline{p}$ in (3) – (11) with $\mathrm{diag}\{{\lambda }_0,{\lambda }_0^{-1}\}$ and ${\overline{p}}_0$ , respectively. In other words, we assume that the expansions (3) – (11) are around the critical configuration $B_{\mathrm{cr}}$ . We shall show that the current approach and the approach of [6] give the same expression for the nonlinear coefficient.

We look for a solution of the form

$$u_j={\epsilon }^2{u}_i^{(0)}{\lambda }_1+ϵu_j^{(1)}(x_1,x_2)+{ϵ}^2{u}_j^{(2)}(x_1,x_2)+{ϵ}^3{u}_j^{(3)}(x_1,x_2)+\cdots ,$$

(30)

$$\tilde{p}={ϵ}^2{\overline{p}}_1{\lambda }_1+ϵp^{(1)}(x_1,x_2)+{ϵ}^2{p}^{(2)}(x_1,x_2)+{ϵ}^3{p}^{(3)}(x_1,x_2)+\cdots ,$$

(31)

where the first terms represent the uniform perturbation from $B_{\mathrm{cr}}$ to $B_{\mathrm{e}}$ , and are given by $u_1^{(0)}=x_1/{\lambda }_0$ , $u_2^{(0)}=-x_2/{\lambda }_0$ , whereas the leading-order solution takes the form

$$u_1^{(1)}={\psi }_{,2},u_2^{(1)}=-{\psi }_{,1},\psi =\mathrm{AH}(x_2)E+\mathrm{c}.\mathrm{c}.,E={\mathrm{e}}^{\mathrm{i}{x}_1},$$

(32)

and

$$p^{(1)}=A{P}_1(x_2)E+\mathrm{c}.\mathrm{c}.,P_1(x_2)=-\mathrm{i}\gamma H″\prime +\mathrm{i}({\mathcal{A}}_{1111}^1-{\mathcal{A}}_{1122}^1-{\mathcal{A}}_{1221}^1)H\prime .$$

(33)

In the last two expressions, c.c. denotes the complex conjugate of the preceding term, A is the unknown amplitude that is to be determined and $H(x_2)$ is given by (24). The expression for E indicates that we have chosen the mode number of the near-critical mode to be unity. This is without loss of generality because the bifurcation condition takes the form $F(k{h}_1,h,\lambda ,r_1,r_2)=0$ and so we can obtain different values of k by varying $h_1$ alone. Alternatively, for each fixed $h_1$ , we may use the inverse of the critical wave number to scale $x_1$ and $x_2$ , and then the wave number becomes unity relative to the scaled coordinates. In the following analysis, $h_1$ is the thickness of the first layer in the critical configuration $B_{\mathrm{cr}}$ and is equal to the critical value of $k{h}_1$ in the linear analysis. We have separated the $O({\epsilon }^2)$ uniform perturbations in (30) and (31) from the other $O({\epsilon }^2)$ terms in order to facilitate comparisons with the expressions in [6].

As the above leading-order solution is simply the linearized solution with the undetermined amplitude chosen to be A, when we substitute (30) – (31) into the nonlinear governing equations and auxiliary conditions, we find that the system of equations obtained by equating the coefficients of $ϵ$ is automatically satisfied. By equating the coefficients of ${ϵ}^2$ , we obtain the second-order governing equations

$$u_{1,1}^{(2)}+u_{2,2}^{(2)}=\frac{1}{2}{u}_{i,j}^{(1)}{u}_{j,i}^{(1)},$$

(34)

$${\mathcal{A}}_{\mathrm{jilk}}^1{u}_{k,\mathrm{lj}}^{(2)}-p_{,i}^{(2)}=-{\mathcal{A}}_{\mathrm{jilknm}}^2{u}_{m,n}^{(1)}{u}_{k,\mathrm{lj}}^{(1)}-p_{,j}^{(1)}{u}_{j,i}^{(1)},$$

(35)

and the corresponding auxiliary conditions

$$T_i^{(2)}=0\mathrm{on}{x}_2=h_1+h_2,$$

(36)

$$[u_i^{(2)}]=0,[T_i^{(2)}]=0,\mathrm{on}{x}_2=0\mathrm{or}{h}_1,$$

(37)

$$u_i^{(2)}\to 0\mathrm{as}{x}_2\to -\infty ,$$

(38)

where

$$T_i^{(2)}={\mathcal{A}}_{2\mathrm{ilk}}^1{u}_{k,l}^{(2)}+{\overline{p}}_0{u}_{2,i}^{(2)}-p^{(2)}{\delta }_{2i}+\frac{1}{2}{\mathcal{A}}_{2\mathrm{ilknm}}^2{u}_{k,l}^{(1)}{u}_{m,n}^{(1)}-{\overline{p}}_0{u}_{2,k}^{(1)}{u}_{k,i}^{(1)}+p^{(1)}{u}_{2,i}^{(1)}.$$

(39)

In writing down the last expression, we have made use of the result ${\mathcal{A}}_{2\mathrm{ilk}}^1{u}_{k,l}^{(0)}+{\overline{p}}_0{u}_{2,i}^{(0)}-{\overline{p}}_1{\delta }_{2i}\equiv 0$ which is the traction-free boundary condition associated with the infinitesimal homogeneous perturbation from $B_{\mathrm{cr}}$ to $B_{\mathrm{e}}$ . As a result, ${\lambda }_1$ does not appear in (34) – (39) which take the same form as in [6].

The second-order problem specified by (34) – (39) can be solved once the form of prestress and the elastic moduli are given. In the next section we explain how this problem can be solved and present explicit results for the special case of uniaxial compression and neo-Hookean materials. However, our following derivation of the amplitude equation does not depend on the explicit solution of this problem.

To derive the amplitude equation which must be satisfied by A, we follow [59] and [61] and make use of the virtual work principle

$${\int }_{-\infty }^{h_1+h_2}\mathrm{d}{x}_2{\int }_0^{2\pi }{\chi }_{\mathrm{ij}}{u}_{i,j}^{*}\mathrm{d}{x}_1=0,$$

(40)

where $u_i^{*}$ is a linear solution corresponding to $k=-1$ . More precisely, we take

$$u_1^{*}=H^{\prime }(x_2){\mathrm{e}}^{-\mathrm{i}{x}_1},u_2^{*}=\mathrm{i}H(x_2){\mathrm{e}}^{-\mathrm{i}{x}_1},-\infty <x_2\le h_1+h_2.$$

(41)

The identity (40) can be proved by integration by part followed by an application of the divergence theorem [6]. The expansions (30) – (31) can now be substituted into (8) and the resulting expression into (40). On equating the coefficients of $ϵ$ , ${ϵ}^2$ and ${ϵ}^3$ , it is found that the two equations obtained from equating the coefficients of $ϵ$ and ${ϵ}^2$ are automatically satisfied. From equating the coefficient of ${ϵ}^3$ , we obtain

$${\int }_{-\infty }^{h_1+h_2}\mathrm{d}{x}_2{\int }_0^{2\pi }\{({\mathcal{L}}_{\mathrm{ij}}[{\mathit{u}}^{(3)},p^{(3)}]+{\sigma }_{\mathrm{ij}}^{(3)})u_{i,j}^{*}\}\mathrm{d}{x}_1=0,$$

(42)

where

$${\mathcal{L}}_{\mathrm{ij}}[{\mathit{u}}^{(3)},p^{(3)}]={\mathcal{A}}_{\mathrm{jilk}}^1{u}_{k,l}^{(3)}+{\overline{p}}_0{u}_{j,i}^{(3)}-p^{(3)}{\delta }_{\mathrm{ji}},$$

(43)

$$\begin{array}{c}{\sigma }_{\mathrm{ij}}^{(3)}={\mathcal{A}}_{\mathrm{jilknm}}^2{u}_{k,l}^{(1)}(u_{m,n}^{(2)}+{\lambda }_1{u}_{m,n}^{(0)})+\frac{1}{6}{\mathcal{A}}_{\mathrm{jilknmqp}}^3{u}_{k,l}^{(1)}{u}_{m,n}^{(1)}{u}_{p,q}^{(1)}\hfill \\ +p^{(1)}(u_{j,i}^{(2)}+{\lambda }_1{u}_{j,i}^{(0)})+({\lambda }_1{\overline{p}}_1+p^{(2)})u_{j,i}^{(1)}.\hfill \end{array}$$

(44)

In writing down (44) we have made use of the identities that

$$(u_{j,k}^{(1)}{u}_{k,i}^{(2)}+u_{j,k}^{(2)}{u}_{k,i}^{(1)}-u_{j,k}^{(1)}{u}_{k,l}^{(1)}{u}_{l,i}^{(1)})u_{i,j}^{*}=0,u_{j,k}^{(1)}{u}_{k,i}^{(1)}{u}_{i,j}^{*}=0,$$

which can be verified by expanding the summations and then making use of the properties $u_{j,j}^{(1)}=u_{j,j}^{*}=0$ and (34).

Only the first term in the integrand of (42) now contains the unknown third-order solution $({\mathit{u}}^{(3)},p^{(3)})$ . By integrating ${\mathcal{L}}_{\mathrm{ij}}[{\mathit{u}}^{(3)},p^{(3)}]u_{i,j}^{*}$ repeatedly by parts and making use of the fact that $u_i^{*}$ given by (41) is a linear solution, this term can be expressed in terms of the first- and second-order solutions, and (42) then reduces to

$${\int }_{-\infty }^{h_1+h_2}\mathrm{d}{x}_2{\int }_0^{2\pi }\{p^{*}{u}_{j,i}^{(1)}{u}_{i,j}^{(2)}+{\sigma }_{\mathrm{ij}}^{(3)}{u}_{i,j}^{*}\}\mathrm{d}{x}_1=0,$$

(45)

where $p^{*}$ is the pressure field corresponding to $u_i^{*}$ , and is given by

$$p^{*}=\{\mathrm{i}\gamma H″\prime +\mathrm{i}({\mathcal{A}}_{1122}^1+{\mathcal{A}}_{1221}^1-{\mathcal{A}}_{1111}^1)H\prime \}{\mathrm{e}}^{-\mathrm{i}{x}_1}.$$

In obtaining (45), we have also made use of the result $u_{i,i}^{(3)}=u_{i,j}^{(1)}{u}_{j,i}^{(2)},$ which is obtained from equating the coefficients of ${ϵ}^3$ in (11).

To facilitate the remaining presentation, we write

$$u_{i,j}^{(1)}=A{\Gamma }_{\mathrm{ij}}^{(1)}E+\mathrm{c}.\mathrm{c}.,u_{i,j}^{*}={\overline{\Gamma }}_{\mathrm{ij}}^{(1)}\overline{E},$$

(46)

$$p^{(1)}=A{P}_{1}E+\mathrm{c}.\mathrm{c}.,p^{(2)}=A\overline{A}{P}_0+A^2{P}_2{E}^2+\mathrm{c}.\mathrm{c}.,$$

(47)

$$u_{i,j}^{(2)}=A\overline{A}{\Gamma }_{\mathrm{ij}}^{(m)}+A^2{\Gamma }_{\mathrm{ij}}^{(2)}{E}^2+\mathrm{c}.\mathrm{c}.,$$

(48)

where the bars on A and ${\Gamma }_{\mathrm{ij}}$ signify complex conjugation and the expressions for ${\Gamma }_{\mathrm{ij}}^{(1)}$ , ${\Gamma }_{\mathrm{ij}}^{(m)}$ , ${\Gamma }_{\mathrm{ij}}^{(2)}$ , $P_1$ , $P_0$ and $P_2$ can be obtained from the leading-order and second-order solutions.

On substituting (46) – (48) into (45) and evaluating the integral with respect to $x_1$ , we obtain the amplitude equation

$$c_0{\lambda }_{1}A+c_1|A{|}^{2}A=0,$$

(49)

where the linear and nonlinear coefficients $c_0$ and $c_1$ are given, respectively, by

$$c_0={\int }_{-\infty }^{h_1+h_2}\{{\mathcal{A}}_{\mathrm{jilknm}}^2{u}_{m,n}^{(0)}{\Gamma }_{\mathrm{kl}}^{(1)}+P_1{u}_{j,i}^{(0)}+{\overline{p}}_1{\Gamma }_{\mathrm{ji}}^{(1)}\}{\overline{\Gamma }}_{\mathrm{ij}}^{(1)}\mathrm{d}{x}_2,$$

(50)

$$c_1={\int }_{-\infty }^{h_1+h_2}\{{\overline{P}}_1({\Gamma }_{\mathrm{ij}}^{(1)}{\Gamma }_{\mathrm{ji}}^{(m)}+2{\overline{\Gamma }}_{\mathrm{ij}}^{(1)}{\Gamma }_{\mathrm{ji}}^{(2)})+K_{\mathrm{ij}}{\overline{\Gamma }}_{\mathrm{ij}}^{(1)}\}\mathrm{d}{x}_2.$$

(51)

In the above expression for $c_1$ , $K_{\mathrm{ij}}$ is given by

$$\begin{array}{c}{K}_{\mathrm{ij}}={\mathcal{A}}_{\mathrm{jilknm}}^2({\Gamma }_{\mathrm{kl}}^{(1)}{\Gamma }_{\mathrm{mn}}^{(m)}+{\overline{\Gamma }}_{\mathrm{kl}}^{(1)}{\Gamma }_{\mathrm{mn}}^{(2)})+\frac{1}{2}{\mathcal{A}}_{\mathrm{jilknmqp}}^3{\overline{\Gamma }}_{\mathrm{kl}}^{(1)}{\Gamma }_{\mathrm{mn}}^{(1)}{\Gamma }_{\mathrm{pq}}^{(1)}\hfill \\ +P_1{\Gamma }_{\mathrm{ji}}^{(m)}+P_0{\Gamma }_{\mathrm{ji}}^{(1)}+P_2{\overline{\Gamma }}_{\mathrm{ji}}^{(1)}.\hfill \end{array}$$

(52)

The amplitude equation (49) admits the non-trivial post-buckling solution

$$|A{|}^2=-\frac{c_0}{c_1}{\lambda }_1.$$

It can be shown [6] that $c_0$ is always positive. Thus, the above solution can be obtained only if ${\lambda }_1/c_1<0$ . It then follows that the bifurcation is super-critical if $c_1>0$ and sub-critical if $c_1<0$ .

On comparing the above expressions for $c_0$ and $c_1$ with those in Cai and Fu [6], we see that the two expressions for $c_1$ are identical, but those for $c_0$ are different. The discrepancy can be explained by noting that if we were to make a variable transformation from the coordinates in $B_{\mathrm{e}}$ to $B_{\mathrm{cr}}$ in the expression for ${\chi }_{\mathrm{ij}}$ in [6], the following extra term would be produced to order ${\epsilon }^3$ :

$${\epsilon }^3\frac{{\lambda }_1}{{\lambda }_0}\{{(-1)}^l{\mathcal{A}}_{\mathrm{jilk}}^1{u}_{k,l}^{(1)}+{(-1)}^i{\overline{p}}_0{u}_{j,i}^{(1)}\},(\mathrm{no}\mathrm{summation}\mathrm{on}i\mathrm{but}\mathrm{summation}\mathrm{on}l).$$

(53)

We have verified numerically that when this term is added in the evaluation of the virtual work principle in [6], the approach used by Cai and Fu [6] gives the same result for $c_0$ as the approach adopted in the current paper. We note, however, that this discrepancy is immaterial because it is the sign of $c_1$ that determines whether the bifurcation is super-critical or sub-critical.

4.2. Uniaxial compression and neo-Hookean materials

In this subsection, we calculate the coefficients in the amplitude equation (49) for the special case when the layers and half-space are made of different neo-Hookean materials and the prestress takes the form of a uniaxial compression. We shall present numerical results for the three cases considered in Figures 1(a), 2(a) and 3(a).

We assume that the maximum ${\lambda }_0$ in (29) is attained at $k{h}_1=h_{1\mathrm{cr}}$ . As we have taken $k=1$ , this implies that $h_1=h_{1\mathrm{cr}}$ . We note that both ${\lambda }_0$ and $h_{1\mathrm{cr}}$ depend on $r_1$ and $r_2$ , and ${\overline{p}}_1$ is related to ${\lambda }_1$ by ${\overline{p}}_1=-2\mu {\lambda }_1/{\lambda }_0^3$ . Our aim in this subsection is to determine the dependence of $c_1$ on $r_1$ and $r_2$ .

With the aid of (32) – (33) and (46)₁, we obtain

$${\Gamma }_{11}^{(1)}=\mathrm{i}{H}^{\prime },{\Gamma }_{12}^{(1)}=H^{″},{\Gamma }_{21}^{(1)}=H,{\Gamma }_{22}^{(1)}=-\mathrm{i}{H}^{\prime },$$

(54)

and

$$P_1(x_2)=\mathrm{i}\mu ({\lambda }^2{H}^{\prime }-{\lambda }^{-2}{H}^{″\prime }).$$

(55)

The governing equations (34) and (35) for the second-order solution reduce to

$$u_{1,1}^{(2)}+u_{2,2}^{(2)}=(u_{1,1}^{(1)}{)}^2+u_{1,2}^{(1)}{u}_{2,1}^{(1)},\mu {\overline{B}}_{\mathrm{jl}}{u}_{i,\mathrm{jl}}^{(2)}-p_{,i}^{(2)}=-p_{,j}^{(1)}{u}_{j,i}^{(1)},$$

(56)

where $\overline{\mathit{B}}=\mathrm{diag}\{{\lambda }_0^2,{\lambda }_0^{-2}\}$ . These equations are to be solved subject to the auxiliary conditions (36) – (38), where the expression (39) for $T_i^{(2)}$ now reduces to

$$T_i^{(2)}=\mu {\overline{B}}_{2l}{u}_{i,l}^{(2)}+{\overline{p}}_0{u}_{2,i}^{(2)}-p^{(2)}{\delta }_{2i}-{\overline{p}}_0{u}_{2,k}^{(1)}{u}_{k,i}^{(1)}+p^{(1)}{u}_{2,i}^{(1)}.$$

(57)

Owing to quadratic interaction, the right-hand side of (56)₁ is a linear combination of $E^0$ , $E^2$ and $E^{-2}$ . Thus, the solution for $u_i^{(2)}$ takes the form

$$u_1^{(2)}=A\overline{A}{U}_0(x_2)+A^2{U}_2(x_2)E^2+\mathrm{c}.\mathrm{c}.,$$

(58)

$$u_2^{(2)}=A\overline{A}{V}_0(x_2)+A^2{V}_2(x_2)E^2+\mathrm{c}.\mathrm{c}.,$$

(59)

and $p^{(2)}$ takes the form (47)₂. On substituting (58) – (59) and (47)₂ into (56), we find, after some manipulation,

$$U_0=0,V_0=2H{H}^{\prime },2\mathrm{i}{U}_2=-V_2^{\prime }+H{H}^{″}-H^{\prime 2},$$

(60)

$$P_0=2\mu {\lambda }^{-2}(H{H}^{\prime }{)}^{\prime }+2\mathrm{i}{P}_1{H}^{\prime },2\mathrm{i}{P}_2=\mu {\lambda }^{-2}{U}_2^{″}-4\mu {\lambda }^2{U}_2+P_1^{\prime }H-P_1{H}^{\prime },$$

(61)

$$V_2^{″″}-4(1+{\lambda }^4)V_2^{″}+16{\lambda }^4{V}_2=3({\lambda }^4-1)(H^{\prime }{H}^{″}-H{H}^{″\prime }).$$

(62)

It was shown in [6] that (62) has a particular integral given by

$$V=\Gamma (x_2)\equiv \frac{{\lambda }^4-1}{9{\lambda }^8-82{\lambda }^4+9}(9H{H}^{″\prime }-21H^{\prime }{H}^{″}).$$

(63)

Thus, the general solution to (62) is given by

$$V_2=\{\begin{array}{cc}{B}_1\exp (2s_1{x}_2)+B_2\exp (2s_2{x}_2)+\Gamma (x_2),\hfill & x_2\in (-\infty ,0)\hfill \\ \multicolumn{1}{c}{\sum _{j=1}^4{\hat{B}}_j\exp (2s_j{x}_2)+\Gamma (x_2),}& x_2\in (0,h_1)\hfill \\ \multicolumn{1}{c}{\sum _{j=1}^4{\tilde{B}}_j\exp (2s_j{x}_2)+\Gamma (x_2),}& x_2\in (h_1,h_1+h_2),\hfill \end{array}$$

(64)

where $s_1=1,s_2={\lambda }_0^2$ , $s_3=-1,s_4=-{\lambda }_0^2$ , and $B_1$ , $B_2$ , ${\hat{B}}_1$ – ${\hat{B}}_4$ and ${\tilde{B}}_1$ – ${\tilde{B}}_4$ are constants that are determined by the auxiliary conditions (36) – (38).

With the aid of (58) – (59), we may calculate $u_{i,j}^{(2)}$ . Comparing the resulting expressions with (48) then yields

$$\begin{array}{c}{\Gamma }_{11}^{(m)}=0,{\Gamma }_{12}^{(m)}=U_0^{\prime },{\Gamma }_{21}^{(m)}=0,{\Gamma }_{22}^{(m)}=V_0^{\prime },\hfill \\ {\Gamma }_{11}^{(2)}=2\mathrm{i}{U}_2,{\Gamma }_{12}^{(2)}=U_2^{\prime },{\Gamma }_{21}^{(2)}=2\mathrm{i}{V}_2^{\prime },{\Gamma }_{22}^{(2)}=V_2^{\prime }.\hfill \end{array}$$

(65)

The expressions for $c_0$ and $c_1$ can be simplified further by noting that the second- and third-order elastic moduli are all zero.

We now investigate the dependence of $c_1$ on $r_1$ and $r_2$ for the three representative cases shown in Figures 1(a), 2(a) and 3(a). As in the linear analysis, some general results may be deduced by referring to the three limits, $r_1\to 1$ , $r_2\to r_1$ and $r_2\to \infty$ , under which the structure under consideration reduces to a half-space coated by a single layer. As for the latter reduced case $c_1$ vanishes when the modulus ratio is equal to $0.57$ , we may deduce for the current problem that the curve corresponding to $c_1=0$ must contain the point $(r_1,r_2)=(0.57,0.57)$ and must approach the asymptote $r_1=0.57$ as $r_2\to \infty$ . These facts are used to validate the Mathematica code that is used to compute $c_1$ for any choice of $r_1,r_2$ and h for which a stretch maximum exists.

Figure 1(b) shows the sign of $c_1$ in the $(r_1,r_2)$ -plane when $h=1$ . The plane is divided into three regions by two solid curves, and the three regions correspond to $c_1>0$ , $c_1<0$ , and non-existence of a stretch maximum, respectively. In addition to the general observations made above, the blue solid curve also tends to an asymptote as $r_1\to \infty$ . Fitting the numerical results for $2.1<r_1<10$ to a straight line, we obtain $r_2=0.38r_1+0.24$ which is displayed in Figure 1(b). Our results conform with the expectation that the bifurcation will be super-critical if both layers are much stiffer than the half-space (corresponding to the area near the origin in Figure 1(b)), but there are also two novel aspects. First, even if both layers are softer than the half-space (i.e. $r_1>1,r_2>1$ ), the bifurcation can still be super-critical (and so robust wrinkling patterns can be observed) provided $r_2<0.38r_1+0.24$ , that is if the top layer is sufficiently harder than the first layer. Second, if $r_2$ is fixed to lie in the interval $(0.414,1)$ and $r_1$ is increased from zero, then $c_1$ changes sign twice: it is positive for sufficiently small or large values of $r_1$ , but is negative in between. Similarly, when $r_1$ is fixed to be between $0.57$ and $0.648$ and $r_2$ is increased from zero, the $c_1$ also changes sign twice. Thus, adding an extra layer enables robust wrinkling patterns to be achieved over a larger parameter regime.

Figure 2(b) displays the corresponding results when $h=0.1$ . The asymptotes associated with the limits $r_1\to \infty$ and $r_2\to \infty$ are similar to those in the previous case, but now the curve corresponding to $c_1=0$ splits into two branches owing to the presence of mode switching. These two branches are connected by the green line across which the sign of $c_1$ changes abruptly owing to mode switching. On the other two segments of the red dashed line across which mode jumping takes place, the sign of $c_1$ remains unchanged when the line is crossed; see later discussion related to Figure 7. Note that the vertical asymptote $r_2=0.57$ lies between $f_2$ and $g_2$ . Thus, when $r_1$ is fixed to lie in the interval $(f_2,0.57)$ and $r_2$ is increased from zero gradually, the nature of bifurcation changes according to super-critical → sub-critical → super-critical, whereas when $r_1$ is fixed to lie in the interval $(0.57,0.60)$ it evolves as super-critical → sub-critical → super-critical → sub-critical.

Figure 7.

Comparison of the dependence of $c_1/c_0$ on $r_1$ when the Gent material model with $J_m=30$ or the neo-Hookean model is used. Two typical mode jumping behaviours are shown: (a) $r_2=0.5,h=0.1$ , here $c_1$ does not change sign when mode jumping takes place at $r_1=0.72$ , ${\lambda }_0=0.63$ , with associated wave numbers given by $(k{h}_1{)}_{\mathrm{cr}}=1.39,14.57$ ; (b) $r_2=0.25,h=0.1$ , here $c_1$ changes sign when mode jumping takes place at $r_1=0.52$ , ${\lambda }_0=0.70$ , with associated wave numbers given by $(k{h}_1{)}_{\mathrm{cr}}=1.17,12.34$ .

Finally, in Figure 3(b) we display the results for $h=10$ . In the limit $r_2\to \infty$ , the black and blue lines asymptote to $r_1=1$ and $r_1=0.57$ , respectively, whereas the red dotted line intersects the black line at $r_1=0.75$ approximately, after which the red curve stays on the black curve (so no mode switching takes place for $r_1$ large enough). The curve corresponding to $c_1=0$ again splits into two branches owing to the presence of mode switching. If we now fix $r_1$ and increase $r_2$ from zero gradually, the bifurcation behaviour is again dependent on the fixed value of $r_1$ but even more complicated than in the previous case. In both cases when mode switching is possible, the mode switching lines (the red dashed lines) consist of three distinctive parts: a part that is entirely in the domain of $c_1>0$ , a part that is entirely in the domain of $c_1<0$ , and a part across which the sign of $c_1$ changes. As bifurcation with $c_1<0$ is sensitive to imperfections, we expect that mode switching can only take place in a predictable fashion across the first part.

4.3. Uniaxial compression and general material model

When more general strain energy functions are used, solution of the second order problem follows the same procedure as in the previous section. The only major difference is that a particular integral as simple as (63) does not seem to be possible. Instead the required particular integral in (64) can be found using the method of variation of parameters [42]. If the equation for $V_2$ is written in the form

$$V_2^{″″}-(1+s^2)a^2{V}_2^{″}+s^2{a}^4{V}_2=\omega (x_2),$$

(66)

where a and s are known constants, then the particular integral is given by

$$\begin{array}{c}\Gamma (x_2)=\frac{1}{2a^{3}s(1-s^2)}\{s{\mathrm{e}}^{a{x}_2}\int {\mathrm{e}}^{-a{x}_2}\omega (x_2)d{x}_2-{\mathrm{e}}^{\mathrm{as}{x}_2}\int {\mathrm{e}}^{-\mathrm{as}{x}_2}\omega (x_2)d{x}_2\hfill \\ +{\mathrm{e}}^{-\mathrm{as}{x}_2}\int {\mathrm{e}}^{\mathrm{as}{x}_2}\omega (x_2)d{x}_2-s{\mathrm{e}}^{-a{x}_2}\int {\mathrm{e}}^{a{x}_2}\omega (x_2)d{x}_2\},\hfill \end{array}$$

(67)

with the understanding that the arbitrary constants in the indefinite integrals are all set to zero (this is necessary in order to satisfy the decay condition as $x_2\to -\infty$ ). This expression is valid for the two layers as well as the substrate although $\omega (x_2)$ takes different expressions in the three different regions.

We have written a separate Mathematica code based on the above formula to compute $c_1$ for any material model. The program is used to validate the program used in the previous subsection which is written specifically and independently for the neo-Hookean material model. As an illustrative example to show the effect of the material extensibility $J_m$ , in Figure 7(a) and (b) we have shown $c_1/c_0$ against $r_1$ when $h=0.1$ and $r_2$ is fixed to be $0.5$ and $0.25$ , respectively. The leading-order solution is normalized such that $u_2=A\cos x_1$ at $x_2=h_1+h_2$ . The results are displayed against their counterparts when the neo-Hookean model is used for which the results are already given in Figure 2(b). It is seen that decreasing $J_m$ has little effect on the value of $r_2$ at which mode jumping takes place and only slightly widens the interval of $r_1$ where $c_1$ is negative.

5. Conclusion

In this paper, we have investigated the linear and nonlinear buckling properties of a hyperelastic half-space coated with two layers. At macro-scales, buckling usually undermines a structure’s integrity and should be avoided. When a structure is sensitive to imperfections, any imperfection, material or geometrical, will significantly reduce the critical load at which bifurcation takes place. Thus, from a practical point of view, it is important to find the parameter regime in which the structure is imperfection sensitive. At micrometre and sub-micrometre scales, robust wrinkling patterns can be harnessed to serve useful purposes. As only super-critical bifurcations may be observable/realizable in practice, results from our weakly nonlinear analysis provide a road map on how to choose a variety of combinations of material parameters to achieve robust wrinkling patterns.

Our analysis is conducted with the aid of the exact theory of nonlinear elasticity and for general strain energy functions. A Mathematica code is written for computing the coefficient $c_1$ the sign of which determines whether the bifurcation is super-critical or not. For the current two-layers/substrate structure, $c_1$ depends on the modulus ratios $r_1$ and $r_2$ as well as the thickness ratio h. For each fixed h, we may display the sign of $c_1$ in the $(r_1,r_2)$ -plane, covering all the possibilities. Illustrative results are presented for the case when the material is modelled by the neo-Hookean model or the Gent model, and the prestress takes the form of a uniaxial compression.

When the neo-Hookean model is used, three sets of representative results are presented corresponding to $h=0.1,1$ and 10, respectively. They illustrate the three possibilities of no mode switching (when $h=1$ ), mode switching occurring once (when $h=0.1$ ) and mode switching occurring twice (when $h=10$ ), respectively. For each case, we display in the $(r_1,r_2)$ -plane domains where the stretch has a maximum and where $c_1$ is positive. One important finding is that when mode switching is theoretically possible based on the linear analysis, it may not be observable/realisable/controllable if it occurs on a part of the mode switching line where $c_1$ is negative or changes sign. When the Gent model is used, we determine the effects of varying the extensibility parameter $J_m$ and it is found that changing $J_m$ does not seem to change our results in any qualitative way, and the quantitative differences it makes are still insignificant when $J_m$ has become as small as 30.

We remark that although we have presented some representative behaviours, our numerical calculations are by no means intended to be exhaustive. For instance, we cannot conclude whether three or more stretch maxima can occur or not for other parameter combinations. Neither have we considered the effects of allowing a pre-stretch in the substrate [21]. The main aim of this paper has been to demonstrate that the sign of $c_1$ can be computed semi-analytically, with the aid of Mathematica, without making any approximations even for the most general material model. Our Mathematica code is freely available to any interested reader upon request.