A refined model for the buckling of film/substrate bilayers

Abstract

The classical reduced model for film/substrate bilayers is the one in which the film is governed by the Euler–Bernoulli beam equation, and the substrate is replaced by an array of springs (the so-called Winkler foundation assumption). We derive a refined model in which the normal and shear tractions at the bottom of the film are expressed in terms of the corresponding horizontal and vertical displacements, and the response of the half-space is described by the exact theory. The self-consistency of the refined model is confirmed by showing that it yields a four-term (in the incompressible case) or six-term (in the compressible case) expansion for the critical strain that agrees with the expansion given by the exact theory.

1. Introduction

The buckling and post-buckling of a film/substrate bilayer, which is often idealised as a coated hyperelastic half-space, has received much attention in recent decades due to applications ranging from cell patterning [1], optical gratings [2–4], and creation of surfaces with desired wetting and adhesion properties [5–7], to the deduction of material properties of ultrathin films [8,9]. Under the framework of non-linear elasticity, early studies include the linear analyses by Dorris and Nemat-Nasser [10], Shield et al. [11], Ogden and Sotiropoulos [12], Bigoni et al. [13], Steigmann and Ogden [14], and the non-linear post-buckling analysis by Cai and Fu [15]. More recent studies have addressed the phenomenon of period doubling [16], effects of pre-stretching in the half-space [17], compressibility [18,19], multi-layering [20–22], growth [23–25], and magnetorheology [26]. There also exists a large body of literature based on various approximate models for the layer and substrate (see, for example, [27–30] and the references therein).

With the use of the exact theory of non-linear elasticity, Cai and Fu [15] derived the following asymptotic expansion for the critical strain when a neo-Hookean film bonded to a much softer neo-Hookean substrate, both in a state of plane strain, is subjected to a uni-axial compression:

$$1-\lambda =\frac{r}{2\mathrm{kh}}+\frac{1}{12}(\mathrm{kh}{)}^2-\frac{13}{480}(\mathrm{kh}{)}^4-\frac{3}{8}·\frac{r^2}{{(\mathrm{kh})}^2}+O({(\mathrm{kh})}^5),$$

(1)

where $\lambda$ is the stretch ratio in the direction parallel to the surface (so that $1-\lambda$ is the strain), $k$ is the wave number of the buckling mode, $h$ is the layer thickness, and $r={\mu }_{\mathrm{s}}/{\mu }_{\mathrm{f}}$ with ${\mu }_{\mathrm{f}}$ and ${\mu }_{\mathrm{s}}$ denoting the shear moduli for the film and substrate, respectively. Note that the above expansion is only valid when $r=O((\mathrm{kh}{)}^3)$ , in which case the first two terms on the right-hand side of equation (1) are of the same order, and the last two terms are of the same order. Solving this expansion together with $d\lambda /d(\mathrm{kh})=0$ , we may determine the critical stretch and associated wave number where the strain $1-\lambda$ attains its minimum. This two-term result extends the leading order result of Allen [31] for the plane-stress case.

It can be shown that the leading order version of equation (1), with only the first two terms on the right-hand side retained, can be derived from the classical model

$$\frac{\mathrm{EI}}{1-{\nu }^2}y″″(x)+\frac{\mathrm{Eh}}{1-{\nu }^2}(1-\lambda )y″(x)=-\gamma \mathrm{ky},\gamma =\frac{2E_s(1-{\nu }_s)}{(1+{\nu }_s)(3-4{\nu }_s)},$$

(2)

by taking the incompressibility limit $\nu ,{\nu }_s\to 1/2$ , where $E$ and $\nu$ are Young’s modulus and Poisson’s ratio for the layer, respectively, $E_s$ and ${\nu }_s$ their counterparts for the half-space, and $I$ ( $=h^3/12$ ) is the second moment of area of the layer cross section. It was shown in Cai and Fu [32] that equation (2) can be derived from the exact theory of non-linear elasticity using a self-consistent asymptotic procedure (see also discussions in Audoly and Boudaoud [29]). When the half-space is absent, the right-hand side of equation (2) is zero, and the resulting equation is the well-known Euler–Bernoulli beam equation, which has been justified in many studies (see, for example, [33–35] and the references therein).

For a stand-alone beam that is only subject to a force $P$ along its axis, the Timoshenko beam theory is known to provide a self-consistent refinement on the Euler–Bernoulli beam theory [36,37]. It is then tempting to model the response of the layer using the Timoshenko beam theory but still model the response of the half-space using the Winkler assumption. This possibility has indeed been analysed by Erbas et al. [38] for the dynamic case without prestress (see also [39]). Adapting equation (17) in Goldenveizer et al. [40] to our static case, we obtain

$$\frac{\mathrm{EI}}{1-{\nu }^2}\{y″″+\frac{(8-3\nu )}{40(1-\nu )}{h}^2{y}^{(6)}\}+\frac{\mathrm{Eh}}{1-{\nu }^2}(1-\lambda )y″+\gamma \mathrm{ky}=0,$$

(3)

where $y^{(6)}$ denotes the sixth-order derivative of $y$ . However, on substituting a periodic solution of the form $y(x)={\mathrm{e}}^{\mathrm{i}\mathrm{kx}}$ into this equation, it is found that the resulting expression for $1-\lambda$ agrees with the first three terms in equation (1), but not the fourth term. This means that this reduced model has not taken care of the response of the foundation in a self-consistent manner. An earlier attempt at refining the leading order model (2) was made by Shield et al. [11]. However, as was shown by Cai and Fu [32], their model is not self-consistent either. The main purpose of this study is to modify the Winkler foundation assumption and allow the interface to transmit both vertical and shear tractions so that the resulting refined model can correctly predict the entire right-hand side of equation (5). To this end, we derive a Timoshenko-type plate theory for a pre-stressed plate using two different approaches, a standard expansion method as used by Cai and Fu [32], and an expansion scheme developed by Dai and Song [41] where a consistent plate theory was proposed for compressible hyperelastic materials. The latter methodology has been applied to derive consistent plate models for incompressible materials [42], growing solids [43], and nematic liquid crystal elastomers [44]. Recently, Wang et al. [45] have shown that the consistent plate model in Dai and Song [41] can recover most existing plate equations by assuming appropriate order relations between the applied load and the aspect ratio. In deriving the above-mentioned Timoshenko-type plate theory for a pre-stressed plate, we naturally obtain an additional connection between the shear traction and the two displacement components at the interface, which then enables us to apply traction and displacement continuity exactly. Thus, our reduced model consists of two differential equations for the displacement components at the interface.

Our reduced model will be derived for a general hyperelastic material, but to illustrate our results, we shall consider the case in which both the film and substrate are described by the strain-energy function

$$W=\frac{1}{2}\mu (I_1-2-2\log J)+\frac{\mu \nu }{1-2\nu }{(J-1)}^2,$$

(4)

where $\mu$ and $\nu$ are the ground state shear modulus and Poisson’s ratio, respectively, $I_1$ is the first principal invariant of the Cauchy–Green strain tensors and $J$ is the determinant of the deformation gradient. We denote the shear moduli for the film and substrate by $\mu$ and ${\mu }_{\mathrm{s}}$ , respectively, but assume that $\nu$ takes the same value for both the film and substrate. A calculation with the aid of Mathematica [46], extending the analysis in Cai and Fu [15], reveals that the compressive strain now has the asymptotic expansion

$$1-\lambda =\frac{2{(\nu -1)}^2}{3-4\nu }\frac{r}{\mathrm{kh}}+\frac{1}{12}(\mathrm{kh}{)}^2+d_{0}r+d_1(\mathrm{kh}{)}^4+d_{2}r(\mathrm{kh})+d_3\frac{r^2}{{(\mathrm{kh})}^2}+O({(\mathrm{kh})}^5),$$

(5)

where

$$d_0=\frac{(1-\nu )(1-2\nu )}{3-4\nu },d_1=\frac{31{\nu }^2+33\nu -40{\nu }^3-29}{1440{(\nu -1)}^2},$$

(6)

$$d_2=\frac{23-104\nu +156{\nu }^2-80{\nu }^3}{24{(3-4\nu )}^2},d_3=\frac{2{(\nu -1)}^2(32{\nu }^4-60{\nu }^3+47{\nu }^2-21\nu +5)}{{(4\nu -3)}^3}.$$

(7)

As expected, the expansion equation (5) recovers equation (1) in the incompressible limit $\nu \to 1/2$ . Our aim is to derive a refined model that will yield the first six terms in equation (5) exactly. It will also be demonstrated that including the next term in the asymptotic expansion (1) or (5) can significantly improve its accuracy when the materials are incompressible or nearly incompressible.

We note that the term $d_{0}r$ in the expansion equation (5) is of order $(\mathrm{kh}{)}^3$ , which has no counterpart in the incompressible limit. As a quick comparison, we note that the model adopted by Shield et al. [11] would give the following expansion:

$$1-\lambda =\frac{2{(\nu -1)}^2}{3-4\nu }\frac{r}{\mathrm{kh}}+\frac{1}{12}(\mathrm{kh}{)}^2+d_{0}r+\frac{{(\nu -1)}^2}{2(3-4\nu )}r(\mathrm{kh})+O({(\mathrm{kh})}^6).$$

(8)

This expression correctly predicts the term $d_{0}r$ but differs from equation (5) in the $O((\mathrm{kh}{)}^4)$ terms.

The rest of this paper is organised into four sections as follows. After formulating our problem and deriving the exact bifurcation condition in the next section, we devote section 3 to the derivation of the above-mentioned refined model. We first derive a Timoshenko-like model for the film and then use the Stroh formulation to describe the response of the substrate. In section 4, we specialise the model to the buckling analysis of a film/substrate bilayer and show that the reduced model can indeed predict equation (5). The paper is concluded in section 5 with a summary and a discussion of other possible applications of the current reduced model.

2. Formulation and the exact bifurcation condition

We consider a hyperelastic layer bonded to a hyperelastic half-space. To simplify derivations, we assume for the moment that both the layer and half-space are compressible and are in a state of plane strain. The corresponding results for the incompressible case will be obtained by taking the appropriate limit. We first summarise the governing equations valid for both the layer and half-space.

Consider a general, homogeneous elastic body $B$ composed of a non-heat-conducting elastic material. Three configurations are involved in our analysis: an initial unstressed configuration $B_0$ , a finitely deformed configuration $B_e$ , and the current configuration $B_t$ that is obtained by superimposing a small incremental displacement on $B_e$ . The position vectors of a representative particle are denoted by $\mathit{X}$ , $\mathit{x}$ , and $\tilde{\mathit{x}}$ in $B_0,B_e$ , and $B_t$ , respectively, and their Cartesian coordinates are denoted by $X_A,x_i$ , and ${\tilde{x}}_i$ . We write

$$\tilde{\mathit{x}}=\mathit{x}+\mathit{u},$$

(9)

where $\mathit{u}$ , with components $u_i(x_1,x_2)$ , is the incremental displacement associated with the deformation $B_e\to B_t$ .

The deformation gradients arising from the deformations $B_0\to B_t$ and $B_0\to B_e$ are denoted by $\tilde{\mathit{F}}$ and $\overline{\mathit{F}}$ , respectively, and defined by $d\tilde{\mathit{x}}=\tilde{\mathit{F}}d\mathit{X}$ and $d\mathit{x}=\overline{\mathit{F}}d\mathit{X}$ . It then follows that

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

(10)

where here and henceforth a comma indicates differentiation with respect to the implied spatial coordinate.

It is well-known (see, for example, [47]) that in the absence of body forces, the incremental equilibrium equations may be written in the form

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

(11)

where ${\chi }_{\mathrm{ij}}$ are the Cartesian components of the incremental stress tensor, and its linearised formula is given by

$${\chi }_{\mathrm{ij}}={\mathcal{A}}_{\mathrm{jilk}}{u}_{k,l}.$$

(12)

The ${\mathcal{A}}_{\mathrm{jilk}}$ are the first-order instantaneous elastic moduli defined by (Chadwick and Ogden [48])

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

(13)

where $W$ is the strain-energy function per unit volume. For the isotropic solids under consideration, $W$ can be assumed to be a function of the three principal stretches ${\lambda }_1,{\lambda }_2,{\lambda }_3$ . Then, the moduli can be computed according to the following formulae [49]:

$$\begin{array}{c}{\mathcal{A}}_{\mathrm{iijj}}={\lambda }_i{\lambda }_j{W}_{\mathrm{ij}},\hfill \\ {\mathcal{A}}_{\mathrm{ijkl}}=\frac{{\lambda }_i^2}{{\lambda }_i^2-{\lambda }_j^2}({\lambda }_i{W}_i-{\lambda }_j{W}_j){\delta }_{\mathrm{ik}}{\delta }_{\mathrm{jl}}+\frac{{\lambda }_i{\lambda }_j}{{\lambda }_i^2-{\lambda }_j^2}({\lambda }_j{W}_i-{\lambda }_i{W}_j){\delta }_{\mathrm{il}}{\delta }_{\mathrm{jk}},i\ne j,\hfill \end{array}$$

where $W_i=\partial W/\partial {\lambda }_i,W_{\mathrm{ij}}={\partial }^{2}W/\partial {\lambda }_i\partial {\lambda }_j$ , and so on, and the principal stretches are understood to correspond to the primary deformation $B_0\to B_e$ . It is straightforward to verify that

$${\mathcal{A}}_{1212}-{\mathcal{A}}_{1221}={\overline{\sigma }}_1,{\mathcal{A}}_{2121}-{\mathcal{A}}_{2112}={\overline{\sigma }}_2,$$

(14)

where ${\overline{\sigma }}_1$ ( $={\lambda }_1{W}_1$ ) and ${\overline{\sigma }}_2$ ( $={\lambda }_2{W}_2$ ) are the principal Cauchy stresses.

In terms of the displacement components and in the absence of body forces, the equilibrium equation then takes the form

$${\mathcal{A}}_{\mathrm{jilk}}{u}_{k,\mathrm{lj}}=0.$$

(15)

As indicated earlier, the finite deformation associated with $B_0\to B_e$ only appears through the moduli defined by equation (13). Thus, as it stands, the constitutive equation (12) is of the same form as that for an anisotropic elastic solid [47], and so our results may be adapted to describe anisotropic materials.

For plane-strain problems, $u_3$ is zero, and $u_1$ and $u_2$ are only functions of $x_1$ and $x_2$ . The finite deformations that we consider (e.g., uni-axial compression) are always such that the elastic moduli ${\mathcal{A}}_{\mathrm{jilk}}$ is zero whenever there is an odd number of 1, 2, or 3 in the suffices. For instance, ${\mathcal{A}}_{1112}=0,{\mathcal{A}}_{1123}=0$ , and so on. As a result, the equilibrium equations in equation (15) written out for $i=1,2$ take the following form

$${\mathcal{A}}_{1111}{u}_{1,11}+({\mathcal{A}}_{1122}+{\mathcal{A}}_{2112})u_{2,12}+{\mathcal{A}}_{2121}{u}_{1,22}=0,$$

(16)

$${\mathcal{A}}_{1212}{u}_{2,11}+({\mathcal{A}}_{1221}+{\mathcal{A}}_{2211})u_{1,12}+{\mathcal{A}}_{2222}{u}_{2,22}=0.$$

(17)

The traction vector $\mathit{t}$ on any surface with normal pointing in the positive x₂-direction has components given by $t_i={\chi }_{i2}={\mathcal{A}}_{2\mathrm{ilk}}{u}_{k,l}$ . Thus, we have

$$t_1={\mathcal{A}}_{2112}{u}_{2,1}+{\mathcal{A}}_{2121}{u}_{1,2},t_2={\mathcal{A}}_{2211}{u}_{1,1}+{\mathcal{A}}_{2222}{u}_{2,2}.$$

(18)

We now specialise the above equations to a coated half-space. In the undeformed configuration $B_0$ , the coating and half-space are defined by $0\le X_2<H,H\le X_2<\infty$ , respectively. We assume that the primary deformation $B_0\to B_e$ corresponds to a uni-axial compression with stretch ${\lambda }_1=\lambda$ in the x₁-direction. We take $\lambda$ to be our loading/bifurcation parameter. The stretch ${\lambda }_2$ in the x₂-direction can be determined by solving ${\overline{\sigma }}_2=0$ . It follows from equation (14)₂ that the identity ${\mathcal{A}}_{2121}={\mathcal{A}}_{2112}$ holds for all values of $\lambda$ and will be used to eliminate the component ${\mathcal{A}}_{2121}$ in the subsequent analysis. In the finitely deformed configuration $B_e$ , the coating and half-space are defined by

$$0\le x_2\le h\mathrm{a}\mathrm{nd}h\le x_2<\infty ,$$

respectively, where $h={\lambda }_{2}H$ is the thickness of the coating layer in $B_e$ .

The linearised buckling problem consists of solving the equilibrium equations (16) and (17) subject to (1) the traction-free boundary condition $\mathit{t}=0$ at $x_2=0$ , (2) continuity of $\mathit{t}$ and $\mathit{u}$ at the interface $x_2=h$ , and (3) the decaying condition $\mathit{u}\to 0$ as $x_2\to \infty$ . If the solution is assumed to take the form

$$\mathit{u}=\mathit{w}\left(k{x}_2\right){\mathrm{e}}^{\mathrm{i}k{x}_1},$$

(19)

where $k$ is the wave number and $\mathit{w}$ is to be determined, then solving the above buckling problem yields the bifurcation condition in the implicit form

$$\omega (\lambda ,\mathrm{kh},r)=0,$$

(20)

where the expression for $\omega$ is obtained with the aid of Mathematica. Under the assumption that $r<<1$ , we expect $\mathrm{kh}$ and $1-\lambda$ to be both small. Writing $r=(\mathrm{kh}{)}^3{r}_0$ with $r_0$ an $O(1)$ constant, and looking for an asymptotic solution of the form

$$1-\lambda ={\xi }_1(\mathrm{kh}{)}^2+{\xi }_2(\mathrm{kh}{)}^3+{\xi }_3(\mathrm{kh}{)}^4+\cdots ,$$

(21)

we may find the coefficients ${\xi }_1,{\xi }_2,...$ by substituting equation (21) into equation (20) and equating the coefficients of $\mathrm{kh}$ at successive orders. This yields the asymptotic expansion (5) when the strain-energy function (4) is used.

3. Derivation of the reduced model

We now proceed to derive a refined model that will yield the same asymptotic expansion (5) without having to solve the fully three-dimensional (3D) problem. Both the procedure outlined in Cai and Fu [32] and the expansion method pioneered by Dai and Song [41] have been used to derive the same equations, but we shall present our derivations corresponding to the latter method since they are more compact than those obtained using the former method. We also wish to promote the use of the latter method since this method was often explained and used in more involved situations than the current one, and as a result, its full potential does not seem to have been fully appreciated in the wider community.

3.1. Refined model for the coating layer

Following the strategy proposed by Dai and Song [41], we first look for a series solution of the form

$$u_1=\sum _{n=0}^{\infty }\frac{1}{n!}{U}_n(x_1)x_2^n,u_2=\sum _{n=0}^{\infty }\frac{1}{n!}{V}_n(x_1)x_2^n,$$

(22)

where the coefficient functions $U_0(x_1),V_0(x_1),U_1(x_1),V_1(x_1),$ and so on are to be determined. We note that the leading terms $U_0(x_1)$ and $V_0(x_1)$ are the displacement components at the traction-free surface $x_2=0$ .

On substituting equation (22) into the traction-free boundary condition $\mathit{t}(x_1,0)=0$ with the two components $t_1$ and $t_2$ given by equation (18), we obtain

$$U_1(x_1)=-V{\prime }_0(x_1),V_1(x_1)=-({\mathcal{A}}_{1122}/{\mathcal{A}}_{2222})U{\prime }_0(x_1).$$

(23)

Note that the traction-free boundary condition $\mathit{t}(x_1,0)=0$ is now satisfied exactly due to the fact that $u_1$ and $u_2$ are expanded around $x_2=0$ .

On the other hand, on substituting equation (22) into equations (16) and (17) and equating the coefficients of like powers of $x_2^n$ , we obtain the recurrence relations

$${\mathcal{A}}_{2121}{U}_n(x_1)=-({\mathcal{A}}_{1122}+{\mathcal{A}}_{1221})V_{n-1}^{\prime }(x_1)-{\mathcal{A}}_{1111}{U}_{n-2}^{″}(x_1)=0,$$

(24)

$${\mathcal{A}}_{2222}{V}_n(x_1)=-({\mathcal{A}}_{1122}+{\mathcal{A}}_{1221})U_{n-1}^{\prime }(x_1)-{\mathcal{A}}_{1212}{V}_{n-2}^{″}(x_1)=0,$$

(25)

which are valid for $n=2,3,...$ . It can then be seen that these two recurrence relations together with equation (23) can be used to express the right-hand sides of equation (22) entirely in terms of the leading terms $U_0(x_1)$ , $V_0(x_1)$ , and their derivatives. This process may be carried out to any desired order. As a result, the traction vector $\mathit{t}(x_1,h)$ at the interface may also be expressed in terms of $U_0(x_1)$ , $V_0(x_1)$ , and their derivatives to any desired order in $h$ . These expressions have the following special structure:

$$t_1(x_1,h)=h{c}_1{U}_0^{(2)}+h^2{c}_2{V}_0^{(3)}+h^3{c}_3{U}_0^{(4)}+h^4{c}_4{V}_0^{(5)}+\cdots ,$$

(26)

$$t_2(x_1,h)=h{d}_1{V}_0^{(2)}+h^2{d}_2{U}_0^{(3)}+h^3{d}_3{V}_0^{(4)}+h^4{d}_4{U}_0^{(5)}+\cdots ,$$

(27)

where the superscript “ $(n)$ ” $(n=2,3,...)$ denotes the nth-order derivative with respect to $x_1$ , and the constant coefficients $c_1,d_1,c_2,d_2$ , and so on only depend on the elastic moduli. For instance, the leading order coefficients are given by

$$c_1={\mathcal{A}}_{1122}^2/{\mathcal{A}}_{2222}-{\mathcal{A}}_{1111},d_1={\mathcal{A}}_{1221}^2/{\mathcal{A}}_{2121}-{\mathcal{A}}_{1212}.$$

Observe the special feature that the nth-order derivatives $U_0^{(n)}$ and $V_0^{(n)}$ in equations (26) and (27) are multiplied by $h^{n-1}$ ( $n=2,3,...$ ).

As the next step in deriving the reduced model, we need to express the traction vector $\mathit{t}(x_1,h)$ in terms of the displacement components at $x_2=h$ which, with the use of equation (22), are given by

$$U(x_1)\equiv u_1(x_1,h)=\sum _{n=0}^{\infty }\frac{1}{n!}{U}_n(x_1)h^n,$$

(28)

$$V(x_1)\equiv u_2(x_1,h)=\sum _{n=0}^{\infty }\frac{1}{n!}{V}_n(x_1)h^n,$$

(29)

where the sign “≡” defines the short notations $U(x_1)$ and $V(x_1)$ that will be the only variables in our reduced model. Note that the sums in these two expressions only contain $U_0(x_1),V_0(x_1)$ and their derivatives. To express $t(x_1,h)$ in terms of $U(x_1)$ and $V(x_1)$ (and their derivatives), we now invert equations (28) and (29) to express $U_0(x_1)$ and $V_0(x_1)$ in terms of $U(x_1)$ and $V(x_1)$ .

Viewing the two expressions (28) and (29) as asymptotic expansions for $U$ and $V$ in terms of the small parameter $h$ , we may invert them by looking for an asymptotic solution of the form

$$U_0(x_1)=U(x_1)+h{f}_1(x_1)+h^2{f}_2(x_1)+\cdots ,V_0(x_1)=V(x_1)+h{g}_1(x_1)+h^2{g}_2(x_1)+\cdots ,$$

(30)

where the unknown functions $f_1(x_1),g_1(x_1)$ , and so on can be obtained by substituting equation (30) into equations (28) and (29) and equating the coefficients of like powers of $h$ . It turns out that although the right-hand sides of equations (28) and (29) contain derivatives of $f_1(x_1),g_1(x_1),...$ , the determination of these unknown functions at successive orders only involves the solution of linear algebraic equations, and these unknown functions can all be expressed in terms of $U,V$ and their derivatives. For instance, equating the coefficients of $h$ and $h^2$ , we obtain

$$f_1(x_1)=V^{(1)}(x_1),g_1(x_1)=\frac{{\mathcal{A}}_{1122}}{{\mathcal{A}}_{2222}}{U}^{(1)}(x_1),$$

(31)

$$f_2(x_1)=\left(\frac{{\mathcal{A}}_{1111}{\mathcal{A}}_{2222}+{\mathcal{A}}_{1122}{\mathcal{A}}_{1221}-{\mathcal{A}}_{1122}^2}{2{\mathcal{A}}_{1221}{\mathcal{A}}_{2222}}\right)U^{(2)}(x_1),$$

(32)

$$g_2(x_1)=\left(\frac{{\mathcal{A}}_{1212}-{\mathcal{A}}_{1221}+{\mathcal{A}}_{1122}}{2{\mathcal{A}}_{2222}}\right)V^{(2)}(x_1).$$

(33)

As a result, on substituting equation (30) back into equations (26) and (27), the traction vector $\mathit{t}(x_1,h)$ can be expressed in terms of $U(x_1)$ , $V(x_1)$ , and their derivatives to any desired order in $h$ . These expressions take the form

$$t_1(x_1,h)=h{c}_1{U}^{(2)}+h^2{\hat{c}}_2{V}^{(3)}+h^3{\hat{c}}_3{U}^{(4)}+h^4{\hat{c}}_4{V}^{(5)}+\cdots ,$$

(34)

$$t_2(x_1,h)=h{d}_1{V}^{(2)}+h^2{\hat{d}}_2{U}^{(3)}+h^3{\hat{d}}_3{V}^{(4)}+h^4{\hat{d}}_4{U}^{(5)}+\cdots ,$$

(35)

where the new constant coefficients ${\hat{c}}_2,{\hat{d}}_2,{\hat{c}}_3,{\hat{d}}_3,...$ are also only dependent on the elastic moduli. Again observe the special feature that the nth-order derivatives $U^{(n)}$ and $V^{(n)}$ in these two expressions are multiplied by $h^{n-1}$ ( $n=2,3,...$ ). The above derivation is straightforward. Although the expressions for the higher order coefficients ${\hat{c}}_2,{\hat{d}}_2,{\hat{c}}_3,{\hat{d}}_3,...$ are quite involved, the derivation can easily be implemented on a symbolic manipulation platform such as Mathematica.

As pointed out by Dai and Song [41] (see also [45]), the expansions such as equations (34) and (35) contain all the necessary terms that can be used to recover most of the existing plate theories, although this has not been carried out for higher order plate theories such as the Timoshenko theory. We now use these expansions to derive our refined model.

We use ${\hat{\mathcal{A}}}_{\mathrm{jilk}}$ to denote the moduli for the half-space. We assume that the thickness $h$ , displacement $u_i$ , and coordinate $x_i$ have all been scaled by a typical length scale $L$ . For instance, for analysis of buckling, we shall take $L=1/k$ with $k$ denoting the wavenumber of the buckling mode. We assume the scaled $h$ to be small and consider, as indicated in Section 1, the case when ${\hat{\mathcal{A}}}_{\mathrm{jilk}}=O(h^3{\mathcal{A}}_{\mathrm{jilk}})$ , i.e., the layer is much stiffer than the half-space. It is this parameter regime in which the layer behaviour can be modelled by the Euler–Bernoulli beam theory to leading order. In this case, the layer will deform like a beam in the sense that its transverse displacement will be much larger than its longitudinal displacement; more precisely, $U=O(\mathrm{hV})$ . We expect that such a layer will buckle at small compressive strains, that is with $\lambda \approx 1$ . Thus, we write

$$\lambda =1+h^2\psi ,$$

(36)

where $\psi$ is a constant of $O(1)$ . As a result, the moduli can be expanded as

$${\mathcal{A}}_{\mathrm{jilk}}={{\mathcal{A}}_{\mathrm{jilk}}|}_{\lambda =1}+h^2\psi \mathcal{A}{\prime }_{\mathrm{jilk}}+\frac{1}{2}{h}^4{\psi }^2\mathcal{A}\prime {\prime }_{\mathrm{jilk}}+\cdots ,$$

(37)

where, e.g., $\mathcal{A}{\prime }_{\mathrm{jilk}}$ denotes the derivative of ${\mathcal{A}}_{\mathrm{jilk}}$ with respect to $\lambda$ evaluated at $\lambda =1$ . To simplify notation, we shall write ${{\mathcal{A}}_{\mathrm{jilk}}|}_{\lambda =1}$ simply as ${\mathcal{A}}_{\mathrm{jilk}}$ . Since $\lambda =1$ corresponds to the undeformed state which is isotropic, if we denote the associated Lame constants by ${\lambda }^{*}$ and $\mu$ , we have

$${\mathcal{A}}_{\mathrm{jilk}}={\lambda }^{*}{\delta }_{\mathrm{ji}}{\delta }_{\mathrm{lk}}+\mu ({\delta }_{\mathrm{jl}}{\delta }_{\mathrm{ik}}+{\delta }_{\mathrm{jk}}{\delta }_{\mathrm{il}}).$$

(38)

The derivatives $\mathcal{A}{\prime }_{\mathrm{jilk}}$ are dependent on the strain-energy function used and are not written out here for the sake of brevity.

We look for an asymptotic solution of the form

$$U(x_1)=h{w}_{11}(x_1)+h^2{w}_{12}(x_1)+h^3{w}_{13}(x_1)+h^4{w}_{14}(x_1)+\cdots ,$$

(39)

$$V(x_1)=w_{20}(x_1)+h{w}_{21}(x_1)+h^2{w}_{22}(x_1)+h^3{w}_{23}(x_1)+\cdots ,$$

(40)

where the functions $w_{20},w_{11},w_{21}$ , and so on are to be determined. Since the elastic moduli of the half-space is of order $h^3$ , the tractions exerted by the half-space can at most be of order $h^3$ . It turns out that the vertical traction is indeed of order $h^3$ , but the horizontal traction can only be of order $h^4$ for asymptotic consistency [32]. This, of course, is why the Winkler assumption is valid to leading order. Thus, we write traction continuity in the form

$$t_1(x_1,h)=h^4{\tilde{t}}_1(x_1),t_2(x_1,h)=h^3{\tilde{t}}_2(x_1),$$

with $h^4{\tilde{t}}_1(x_1)$ and $h^3{\tilde{t}}_2(x_1)$ denoting the tractions exerted by the half-space. As a result, equations (34) and (35) may be replaced by

$$h^4{\tilde{t}}_1(x_1)=h{c}_1{U}^{(2)}+h^2{\hat{c}}_2{V}^{(3)}+h^3{\hat{c}}_3{U}^{(4)}+h^4{\hat{c}}_4{V}^{(5)}+\cdots ,$$

(41)

$$h^2{\tilde{t}}_2(x_1)=d_1{V}^{(2)}+h{\hat{d}}_2{U}^{(3)}+h^2{\hat{d}}_3{V}^{(4)}+h^3{\hat{d}}_4{U}^{(5)}+\cdots ,$$

(42)

where we have divided equation (35) by $h$ to ease our presentation.

On substituting equations (39) and (40) together with equations (36) and (37) into equations (41) and (42) and equating the coefficients of like powers of $h$ , we can find $w_{11},w_{21}$ , and so on at successive orders in terms of the traction functions ${\tilde{t}}_1(x_1)$ and ${\tilde{t}}_2(x_1)$ . The two equations at order $h$ are automatically satisfied, whereas the two equations at order $h^2$ are linear algebraic equations and can be solved to give

$$w_{11}^{(2)}=-\frac{1}{2}{w}_{20}^{(3)},$$

(43)

$$\frac{\mu ({\lambda }^{*}+\mu )}{3({\lambda }^{*}+2\mu )}{w}_{20}^{(4)}-\psi {\overline{\sigma }}_1^{\prime }{w}_{20}^{(2)}={\tilde{t}}_2(x_1),$$

(44)

where ${\overline{\sigma }}_1$ has been defined in equation (14) and ${\overline{\sigma }}_1^{\prime }$ stands for the derivative of ${\overline{\sigma }}_1$ with respect to $\lambda$ evaluated at $\lambda =1$ . The leading order result (44) corresponds to an Euler–Bernoulli beam that is only supported by normal forces.

At order $h^3$ , the two linear algebraic equations can be solved to give

$$w_{12}^{(2)}=-\frac{1}{2}{w}_{21}^{(3)},$$

(45)

$$\frac{\mu ({\lambda }^{*}+\mu )}{3({\lambda }^{*}+2\mu )}{w}_{21}^{(4)}-\psi {\overline{\sigma }}_1^{\prime }{w}_{21}^{(2)}=0.$$

(46)

At order $h^4$ , we solve the two equations to obtain

$$w_{13}^{(2)}=-\frac{1}{2}{w}_{22}^{(3)}-\frac{({\lambda }^{*}+\mu )}{6({\lambda }^{*}+2\mu )}{w}_{20}^{(5)}+\frac{{\lambda }^{*}{\overline{\sigma }}_1^{\prime }\psi }{8\mu ({\lambda }^{*}+\mu )}{w}_{20}^{(3)}-\frac{({\lambda }^{*}+2\mu )}{4\mu ({\lambda }^{*}+\mu )}{\tilde{t}}_1(x_1),$$

(47)

$$\frac{\mu ({\lambda }^{*}+\mu )}{3({\lambda }^{*}+2\mu )}{w}_{22}^{(4)}-\psi {\overline{\sigma }}_1^{\prime }{w}_{22}^{(2)}=\frac{1}{2}\tilde{t}{\prime }_1(x_1)-\frac{\mu ({\lambda }^{*}+\mu )}{15({\lambda }^{*}+2\mu )}{w}_{20}^{(6)}+k_4\psi w_{20}^{(4)}+\frac{1}{2}{\overline{\sigma }}_1^{″}{\psi }^2{w}_{20}^{(2)},$$

(48)

where

$$k_4=-\frac{1}{12{({\lambda }^{*}+2\mu )}^2}\{2{\lambda }^{*}({\lambda }^{*}+2\mu ){\overline{\sigma }}_1^{\prime }+{({\lambda }^{*}+2\mu )}^2\mathcal{A}{\prime }_{1111}-2{\lambda }^{*}({\lambda }^{*}+2\mu )\mathcal{A}{\prime }_{1122}+{\lambda }^{*2}\mathcal{A}{\prime }_{2222}\}.$$

Multiplying equations (46) and (48) by $h$ and $h^2$ , respectively, and adding the two resulting equations to equation (44), we obtain, after making use of equation (40),

$$\begin{array}{c}\frac{\mu ({\lambda }^{*}+\mu )}{3({\lambda }^{*}+2\mu )}{V}^{(4)}-\psi {\overline{\sigma }}_1^{\prime }{V}^{(2)}={\tilde{t}}_2(x_1)+\frac{h^2}{2}\tilde{t}{\prime }_1(x_1)+\frac{1}{2}{\overline{\sigma }}_1^{″}{\psi }^2{h}^2{V}^{(2)}\hfill \\ +k_4{h}^2\psi V^{(4)}-\frac{\mu ({\lambda }^{*}+\mu )}{15({\lambda }^{*}+2\mu )}{h}^2{V}^{(6)}+O(h^3).\hfill \end{array}$$

(49)

On the other hand, multiplying equations (43) and (45) by $h^2$ and $h^3$ , respectively, and adding the two resulting equations together, we obtain, after making use of equation (39) and rearranging,

$$\begin{array}{c}{h}^4{\tilde{t}}_1(x_1)=-\frac{4\mu ({\lambda }^{*}+\mu )}{{\lambda }^{*}+2\mu }(h{U}^{(2)}(x_1)+\frac{1}{2}{h}^2{V}^{(3)}(x_1))-\frac{2\mu {({\lambda }^{*}+\mu )}^2}{3{({\lambda }^{*}+2\mu )}^2}{h}^4{V}^{(5)}(x_1)\hfill \\ +\frac{{\lambda }^{*}{\overline{\sigma }}_1^{\prime }}{2({\lambda }^{*}+2\mu )}{h}^2(\lambda -1)V^{(3)}(x_1)+O(h^5),\hfill \end{array}$$

(50)

where we have used equation (36) to eliminate $\psi$ . It can be shown that the same equations for $U$ and $V$ are obtained even after solutions at the next order are incorporated. Thus, the error terms in equations (49) and (50) may be replaced by $O(h^4)$ and $O(h^6)$ , respectively.

Finally, on eliminating ${\tilde{t}}_1(x_1)$ from equation (49) with the use of equation (50) and solving the resulting equation for $t_2(x_1)$ , we obtain

$$\begin{array}{c}{h}^3{\tilde{t}}_2(x_1)=\frac{\mu ({\lambda }^{*}+\mu )}{3({\lambda }^{*}+2\mu )}{h}^3{V}^{(4)}-(\lambda -1)h{\overline{\sigma }}_1^{\prime }{V}^{(2)}+\frac{2\mu ({\lambda }^{*}+\mu )}{{\lambda }^{*}+2\mu }(h^2{U}^{(3)}+\frac{1}{2}{h}^3{V}^{(4)})\hfill \\ +\frac{h^5\mu (6{\lambda }^{*2}+13{\lambda }^{*}\mu +7{\mu }^2)}{15{({\lambda }^{*}+2\mu )}^2}{V}^{(6)}+{\hat{k}}_4{h}^3(\lambda -1)V^{(4)}-\frac{1}{2}h{(\lambda -1)}^2{\overline{\sigma }}_1^{″}{V}^{(2)}+O(h^7),\hfill \end{array}$$

(51)

where

$${\hat{k}}_4=\frac{1}{12{({\lambda }^{*}+2\mu )}^2}\{-{\lambda }^{*}({\lambda }^{*}+2\mu ){\overline{\sigma }}_1^{\prime }+{({\lambda }^{*}+2\mu )}^2\mathcal{A}{\prime }_{1111}-2{\lambda }^{*}({\lambda }^{*}+2\mu )\mathcal{A}{\prime }_{1122}+{\lambda }^{*2}\mathcal{A}{\prime }_{2222}\}.$$

Note that the first two terms on the right-hand side of equation (51) corresponds to the leading order theory, and we have not combined them with the other higher order terms. In particular, although the second term involving $h^3{V}^{(4)}$ is of the same order as the leading terms, the sum $h^2{U}^{(3)}+\frac{1}{2}{h}^3{V}^{(4)}$ together is of $O(h^5)$ , as can be seen from equation (50) where all terms must necessarily be of the same order.

In the next section, alternative expressions for $h^4{\tilde{t}}_1(x_1)$ and $h^3{\tilde{t}}_2(x_1)$ will be derived by considering the equilibrium of the half-space. Equating the expressions given by equations (50) and (51) with these alternative expressions, respectively, will then yield the reduced model.

3.2. Response of the half-space

Having derived the relation between traction and displacement at the interface from the side of the coating layer, we now employ the Stroh formulation [50] to derive its counterpart from the side of half-space. The following derivations are adapted from [51].

The equilibrium equations for the half-space are

$${\hat{\mathcal{A}}}_{1i1k}{\hat{u}}_{k,11}+({\hat{\mathcal{A}}}_{1i2k}+{\hat{\mathcal{A}}}_{2i1k}){\hat{u}}_{k,12}+{\hat{\mathcal{A}}}_{2i2k}{\hat{u}}_{k,22}=0,$$

or equivalently,

$$Q{\hat{\mathit{u}}}_{,11}+(R+R^T){\hat{\mathit{u}}}_{,12}+T{\hat{\mathit{u}}}_{,22}=0,$$

(52)

where the three matrices $T,R$ , and $Q$ are defined by

$$T_{\mathrm{ik}}={\hat{\mathcal{A}}}_{2i2k},R_{\mathrm{ik}}={\hat{\mathcal{A}}}_{1i2k},Q_{\mathrm{ik}}={\hat{\mathcal{A}}}_{1i1k}.$$

(53)

The traction vector $\hat{\mathit{t}}$ on any surface with unit normal $\mathit{n}=({\delta }_{i2})$ is given by

$${\hat{t}}_i={\hat{\chi }}_{i2}={\hat{\mathcal{A}}}_{2\mathrm{ilk}}{\hat{u}}_{k,l}={\hat{\mathcal{A}}}_{2i1k}{\hat{u}}_{k,1}+{\hat{\mathcal{A}}}_{2i2k}{\hat{u}}_{k,2},$$

(54)

or equivalently,

$$\hat{\mathit{t}}=R^T{\hat{\mathit{u}}}_{,1}+T{\hat{\mathit{u}}}_{,2}.$$

(55)

We define the Fourier transform (FT) of $\hat{\mathit{u}}$ and its inverse by:

$$\mathit{z}(k,x_2)\equiv \mathcal{F}[\hat{\mathit{u}}]={\int }_{-\infty }^{\infty }\hat{\mathit{u}}(x_1,x_2){\mathrm{e}}^{-\mathrm{i}k{x}_1}d{x}_1,$$

(56)

$$\hat{\mathit{u}}(x_1,x_2)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\mathit{z}(k,x_2){\mathrm{e}}^{\mathrm{i}k{x}_1}\mathrm{dk},$$

(57)

where the first expression defines the vector function $\mathit{z}(k,x_2)$ . Applying the FT to equation (52), we obtain:

$$T\mathit{z}″+\mathrm{i}k(R+R^T)\mathit{z}\prime -k^{2}Q\mathit{z}=0,$$

(58)

where a prime denotes differentiation with respect to $x_2$ . On substituting a trial solution of the form

$$\mathit{z}(k,x_2)=\mathit{a}{\mathrm{e}}^{\mathrm{i}\mathrm{kp}{x}_2},$$

(59)

into equation (58), we find that the constant scalar $p$ and vector $\mathit{a}$ satisfies the eigenvalue problem

$$(p^{2}T+p(R+R^T)+Q)\mathit{a}=0\mathbf{.}$$

(60)

Under the assumption that ${\hat{\mathcal{A}}}_{\mathrm{jilk}}$ satisfies the strong ellipticity condition, the eigenvalues of $p$ in equation (60) cannot be pure real. For the current plane-strain problem, we denote by $p^{(1)},p^{(2)}$ the two eigenvalues of $p$ with positive imaginary parts and ${\mathit{a}}^{(1)},{\mathit{a}}^{(2)}$ the associated eigenvectors. Then, for $k>0$ , a general solution that satisfies the decaying condition is

$$\mathit{z}(k,x_2)=\sum _{j=1}^2{c}_j{\mathit{a}}^{(j)}{\mathrm{e}}^{\mathrm{i}k{p}^{(j)}{x}_2}=A〈{\mathrm{e}}^{\mathrm{i}\mathrm{kp}{x}_2}〉\mathit{c}=A〈{\mathrm{e}}^{\mathrm{i}\mathrm{kp}{x}_2}〉A^{-1}\mathit{d},$$

(61)

where $c_1,c_2$ are constants,

$$A=[{\mathit{a}}^{(1)},{\mathit{a}}^{(2)}],\mathit{c}={[c_1,c_2]}^T,$$

$\mathit{d}=A\mathit{c}$ , and $〈{\mathrm{e}}^{\mathrm{i}\mathrm{kp}{x}_2}〉$ denotes the diagonal matrix

$$\mathrm{d}\mathrm{iag}\{{\mathrm{e}}^{\mathrm{i}k{p}^{(1)}{x}_2},{\mathrm{e}}^{\mathrm{i}k{p}^{(2)}{x}_2}\}.$$

For $k<0$ , we must choose the conjugates of $p^{(1)},p^{(2)}$ in constructing the general solution in order to satisfy the decaying condition. Thus, we have

$$\mathit{z}(k,x_2)=\overline{A}〈{\mathrm{e}}^{\mathrm{i}k\overline{p}{x}_2}〉{\overline{A}}^{-1}\overline{\mathit{d}},k<0,$$

(62)

where an overbar signifies complex conjugation. It then follows that

$$\mathit{z}(k,x_2)=\overline{\mathit{z}(-k,x_2)},k<0.$$

(63)

With the aid of equations (61) and (62), we deduce that

$$\mathit{z}\prime =\mathrm{i}k·\{\begin{array}{cc}A〈p〉A^{-1}z(k,x_2),\hfill & \mathrm{w}\mathrm{hen}k>0,\hfill \\ \multicolumn{1}{c}{\overline{A}〈\overline{p}〉{\overline{A}}^{-1}z(k,x_2),}& \mathrm{w}\mathrm{hen}k<0.\hfill \end{array}$$

(64)

Thus, we have

$$\mathit{z}\prime =(\mathrm{R}\mathrm{e}G+\mathrm{i}\mathrm{s}\mathrm{gn}(k)\mathrm{I}\mathrm{m}G)\mathcal{F}\left[{\hat{\mathit{u}}}_{,1}\right],$$

(65)

where

$$G=A〈p〉A^{-1},$$

(66)

Re and Im denote the real and imaginary parts, respectively, and sgn is the sign function.

To invert equation (65), we make use of the result that

$$\mathrm{s}\mathrm{gn}(k)=-\frac{1}{\mathrm{i}\pi }\mathcal{F}\left[\frac{1}{x_1}\right],$$

(67)

where the FT is understood to take its principal value. It then follows by inverting equation (65), followed by the use of the convolution theorem, that

$${\hat{\mathit{u}}}_{,2}=\left(\mathrm{R}\mathrm{e}G\right){\hat{\mathit{u}}}_{,1}-\left(\mathrm{I}\mathrm{m}G\right)\frac{1}{\pi x_1}*{\hat{\mathit{u}}}_{,1},$$

(68)

where the star signifies integral convolution. Alternatively, we may write this result as

$${\hat{\mathit{u}}}_{,2}=(\mathrm{R}\mathrm{e}G){\hat{\mathit{u}}}_{,1}+(\mathrm{I}\mathrm{m}G)\mathcal{H}[{\hat{\mathit{u}}}_{,1}],$$

(69)

where $\mathcal{H}$ denotes the Hilbert transform defined by

$$\mathcal{H}[g(x_1)]=\frac{1}{\pi }\mathrm{p}.\mathrm{v}.{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\frac{\mathrm{g}(\mathrm{y})}{\mathrm{y}-{\mathrm{x}}_1}\mathrm{dy}=-\frac{1}{\mathrm{\pi }{\mathrm{x}}_1}*\mathrm{g}({\mathrm{x}}_1).$$

(70)

On eliminating ${\hat{\mathit{u}}}_{,2}$ from equation (55) with the use of equation (69), we obtain

$$\hat{\mathit{t}}=\{R^T+T(\mathrm{R}\mathrm{e}G)\}{\hat{\mathit{u}}}_{,1}+T(\mathrm{I}\mathrm{m}G)\mathcal{H}[{\hat{\mathit{u}}}_{,1}].$$

(71)

Note that the continuity condition $\hat{\mathit{u}}(x_1,h)=\mathit{u}(x_1,h)$ at the interface may be differentiated to yield ${\hat{\mathit{u}}}_{,1}(x_1,h)={\mathit{u}}_{,1}(x_1,h)$ . As a result, when evaluated at the interface $x_2=h$ , the $\hat{\mathit{u}}$ in equation (71) may be replaced by $\mathit{u}(x_1,h)=\{U(x_1),V(x_1){\}}^T$ .

It was shown in Fu [51] that the matrix $G$ is related to the surface impedance matrix $M$ [52] by

$$M=-\mathrm{i}(R^T+\mathrm{TG}),$$

(72)

so that equation (71) may also be rewritten as

$$\hat{\mathit{t}}=-(\mathrm{I}\mathrm{m}M){\hat{\mathit{u}}}_{,1}+(\mathrm{R}\mathrm{e}M)\mathcal{H}[{\hat{\mathit{u}}}_{,1}].$$

(73)

The surface impedance matrix $M$ plays an important role in the surface wave theory [53], and many useful results about $M$ are known. In particular, it is positive definite under the convexity assumption for the elastic moduli, and satisfies the matrix Riccati equation (see Biryukov [54] and Fu and Mielke [55])

$$(M-\mathrm{i}R)T^{-1}(M+\mathrm{i}{R}^{\mathrm{T}})-Q=0.$$

(74)

This matrix equation can be solved analytically when there are enough symmetries in the problem. For the current plane-strain problem where the coordinate axes are in the principal directions of deformation, the three matrices $Q,R,T$ have the simple forms

$$T=\left(\begin{array}{cc}{T}_1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& T_2\hfill \end{array}\right),R=\left(\begin{array}{cc}0\hfill & R_1\hfill \\ \multicolumn{1}{c}{R_2}& 0\hfill \end{array}\right),Q=\left(\begin{array}{cc}{Q}_1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& Q_2\hfill \end{array}\right),$$

(75)

where

$$\begin{array}{c}{T}_1={\hat{\mathcal{A}}}_{2121},T_2={\hat{\mathcal{A}}}_{2222},R_1={\hat{\mathcal{A}}}_{1122},\hfill \\ R_2={\hat{\mathcal{A}}}_{2112},Q_1={\hat{\mathcal{A}}}_{1111},Q_2={\hat{\mathcal{A}}}_{1212}.\hfill \end{array}$$

As a result, equation (74) can be solved explicitly to yield

$$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),M_i\mathrm{r}\mathrm{eal},$$

(76)

with

$$\begin{array}{c}{M}_1=\sqrt{T_1{Q}_1-\frac{T_1}{T_2}{\left(\frac{R_1+R_2}{1+\zeta }\right)}^2},\zeta =\sqrt{\frac{T_1{Q}_2}{T_2{Q}_1}},\hfill \\ M_2=\zeta \frac{T_2}{T_1}{M}_1,M_4=\frac{\zeta R_1-R_2}{1+\zeta }.\hfill \end{array}$$

(77)

Returning to equation (73) and evaluating it at the interface $x_2=h$ , we may replace $\hat{\mathit{t}}$ by $\{h^4{\tilde{t}}_1,h^3{\tilde{t}}_2{\}}^T$ and $\hat{\mathit{u}}$ by $\{U,V{\}}^T$ to obtain

$$-M_4{V}^{(1)}(x_1)+M_1\mathcal{H}[U^{(1)}(x_1)]=h^4{\tilde{t}}_1,$$

(78)

$$M_4{U}^{(1)}(x_1)+M_2\mathcal{H}[V^{(1)}(x_1)]=h^3{\tilde{t}}_2.$$

(79)

With $h^4{\tilde{t}}_1$ and $h^3{\tilde{t}}_2$ given by equations (50) and (51), equations (78) and (79) are two ordinary differential equations for $U(x_1)$ and $V(x_1)$ , and are the refined film/substrate model that we set out to derive.

The classical model corresponds to balancing the term $M_2\mathcal{H}[V^{(1)}(x_1)]$ with the first two terms in $h^3{\tilde{t}}_2$ and is given by

$$M_2\mathcal{H}[V^{(1)}(x_1)]=\frac{\mu ({\lambda }^{*}+\mu )}{3({\lambda }^{*}+2\mu )}{h}^3{V}^{(4)}-(\lambda -1)h{\overline{\sigma }}_1^{\prime }{V}^{(2)}.$$

(80)

The $M_2$ should be evaluated at $\lambda =1$ and so we have

$$M_2=\frac{2{\mu }_{\mathrm{s}}({\lambda }_{\mathrm{s}}+2{\mu }_{\mathrm{s}})}{{\lambda }_{\mathrm{s}}+3{\mu }_{\mathrm{s}}},{\overline{\sigma }}_1^{\prime }=\frac{4\mu ({\lambda }^{*}+\mu )}{{\lambda }^{*}+2\mu },$$

(81)

where ${\lambda }_{\mathrm{s}}$ and ${\mu }_{\mathrm{s}}$ represent the Lame constants for the substrate. For sinusoidal solutions of the form $V(x_1)=A\sin (k{x}_1)$ with $A$ a constant, we have

$$\mathcal{H}[V^{(1)}(x_1)]=\mathrm{Ak}\mathcal{H}[\cos (k{x}_1)]=-\mathrm{Ak}\sin (k{x}_1)=-\mathrm{kV}(x_1).$$

Equation (80) then recovers equation (2).

4. Self-consistency of the refined model

On substituting a periodic buckling solution of the form

$$U(x_1)=A\sin (k{x}_1),V(x_1)=B\cos (k{x}_1),$$

(82)

into the reduced model equations (78) and (79) and cancelling $\sin (k{x}_1)$ and $\cos (k{x}_1)$ , we obtain two homogeneous linear equations for $A$ and $B$ . For a non-trivial solution, we must set the determinant of the coefficient matrix to zero, which then yields the bifurcation condition. By assuming that $\mathrm{kh}$ is small, and both ${\lambda }_{\mathrm{s}}/{\lambda }^{*}$ and ${\mu }_{\mathrm{s}}/\mu$ are of order $(\mathrm{kh}{)}^3$ , we find from the bifurcation condition that $1-\lambda$ has an asymptotic expansion given by

$$\begin{array}{c}1-\lambda =\frac{2{(\nu -1)}^2}{3-4\nu }\frac{r}{\mathrm{kh}}+\frac{1}{12}(\mathrm{kh}{)}^2+d_{0}r+d_1(\mathrm{kh}{)}^4+d_{2}r(\mathrm{kh})\hfill \\ +d_3\frac{r^2}{{(\mathrm{kh})}^2}+d_4(\mathrm{kh}{)}^{2}r+d_5\frac{r^2}{\mathrm{kh}}+O({(\mathrm{kh})}^6),\hfill \end{array}$$

(83)

where the constants $d_0,d_1,d_2$ , and $d_3$ are the same as in equation (5), and the new constants $d_4$ and $d_5$ are given by

$$d_4=\frac{192{\nu }^4-512{\nu }^3+508{\nu }^2-218\nu +33}{24{(4\nu -3)}^3},$$

(84)

$$d_5=\frac{256{\nu }^6-992{\nu }^5+1608{\nu }^4-1380{\nu }^3+644{\nu }^2-147\nu +11}{2{(4\nu -3)}^3}.$$

(85)

Compared with equation (5), the expansion equation (83) contains two extra terms of order $(\mathrm{kh}{)}^5$ . We have checked to verify that the expansion equation (83) is the same as that given by the exact bifurcation condition, thus verifying the self-consistency of the refined model. We now show that the two extra order $(\mathrm{kh}{)}^5$ terms in equation (83) can significantly improve its accuracy.

The maximum stretch is attained when $d\lambda /d(\mathrm{kh})=0$ . On substituting equation (83) and

$$\mathrm{kh}=g_1{r}^{1/3}+g_2{r}^{3/3}+g_3{r}^{4/3}+O(r^{5/3}),$$

(86)

into this equation and then equating the coefficients of like powers of $r$ to zero, we obtain

$$g_1=2^{2/3}\sqrt[3]{3}\sqrt[3]{\frac{{(\nu -1)}^2}{3-4\nu }},g_2=\frac{-1920{\nu }^4+3776{\nu }^3-2504{\nu }^2+500\nu +33}{60{(3-4\nu )}^2},$$

(87)

$$g_3=\frac{256{\nu }^5-640{\nu }^4+592{\nu }^3-208{\nu }^2-2\nu +11}{23^{2/3}{(3-4\nu )}^2\sqrt[3]{8{\nu }^2-14\nu +6}}.$$

(88)

On substituting the expansion equation (86) back into equation (83), we obtain the corresponding expression for the critical stretch

$${\lambda }_{\mathrm{cr}}=1-\frac{3^{2/3}{(1-\nu )}^{4/3}}{{(6-8\nu )}^{2/3}}{r}^{2/3}+\frac{(\nu -1)(2\nu -1)}{4\nu -3}r+g_4{r}^{4/3}+g_5{r}^{5/3}+O(r^{6/3}),$$

(89)

where

$$g_4=\frac{{(1-\nu )}^{2/3}(480{\nu }^4-556{\nu }^3+199{\nu }^2-115\nu +72)}{203^{2/3}\sqrt[3]{6-8\nu }{(3-4\nu )}^2},$$

(90)

$$g_5=-\frac{\sqrt[3]{3-4\nu }(256{\nu }^6-1040{\nu }^5+1796{\nu }^4-1670{\nu }^3+863{\nu }^2-227\nu +22)}{22^{2/3}\sqrt[3]{3}{(1-\nu )}^{2/3}{(4\nu -3)}^3}.$$

(91)

In the incompressibility limit $\nu \to 0.5$ , the expressions for the critical stretch and wavenumber reduce to

$${\lambda }_{\mathrm{cr}}=1-\frac{1}{4}(3r{)}^{2/3}+\frac{11}{160}(3r{)}^{4/3}-\frac{1}{24}(3r{)}^{5/3}+O(r^2),$$

(92)

$$(\mathrm{kh}{)}_{\mathrm{cr}}=(3r{)}^{1/3}+\frac{3}{20}r+O(r^{5/3}).$$

(93)

Since care needs to be taken in using the truncated asymptotic expansion (83) to compute $d\lambda /d(\mathrm{kh})$ , we have checked to verify that all the terms displayed in equations (87)–(93) are the same as those given by the exact theory, and those remainder terms represented by the $O$ symbol cannot be derived using the current reduced model.

We observe that in the incompressible case, we have $g_3=0$ so that adding the $O((\mathrm{kh}{)}^5)$ terms in equation (83) does not result in any change in $(\mathrm{kh}{)}_{\mathrm{cr}}$ but does give rise to an extra $O(r^{5/3})$ term in equation (92). In Figure 1(a) and (b), we have shown the effect of adding the $O(r^{5/3})$ term in equation (92) by comparing the asymptotic results with their exact counterparts. It is seen that adding the $O(r^{5/3})$ term in equation (92) significantly improves the accuracy of the resulting asymptotic expansion for ${\lambda }_{\mathrm{c}r}$ .

Figure 1.

Comparison of asymptotic results (dashed lines) with exact results (black solid lines) in the incompressible case. The red dashed lines in (a) and (b) correspond to equations (92) and (93), respectively, and the blue line in (a) corresponds to equation (92) with the $O(r^{5/3})$ term neglected.

Having demonstrated the importance of the extra $O((\mathrm{kh}{)}^5)$ terms in equation (83), we now investigate their effect for the compressible case. We note that the extra $O((\mathrm{kh}{)}^5)$ terms in equation (83) give rise to the $g_3{r}^{4/3}$ term in equation (86) and the $g_5{r}^{5/3}$ term in equation (89). In Figures 2–4, we have shown the counterparts of Figure 1 for $\nu =0.45,0.3$ and $0.1$ , respectively, where the red and blue dashed lines are asymptotic results with the these extra terms added or neglected, respectively. It is seen that adding the extra terms continues to improve the accuracy for $\nu$ close to $0.5$ (e.g., $\nu =0.45$ ), but it has the opposite effect when $\nu =0.3$ and a mixed effect when $\nu =0.1$ (in the sense that it improves the accuracy in ${\lambda }_{\mathrm{c}r}$ but worsens the accuracy in $(\mathrm{kh}{)}_{\mathrm{cr}}$ ). Thus, we conclude that including the $O((\mathrm{kh}{)}^5)$ terms in equation (83) improves the accuracy of the asymptotic results significantly only when the materials are incompressible or nearly incompressible. For other values of $\nu$ , the optimal truncation of the asymptotic expansion (83) will be dependent on $\nu$ .

Figure 2.

Comparison of asymptotic results (dashed lines) with exact results (black solid lines) when $\nu =0.45$ . The red dashed lines in (a) and (b) correspond to equations (89) and (86), respectively, and whereas the blue lines represent the lower order approximations with $g_3$ and $g_5$ set to zero.

Figure 3.

Comparison of asymptotic results (dashed lines) with exact results (black solid lines) when $\nu =0.3$ . The red dashed lines in (a) and (b) correspond to equations (89) and (86), respectively, and whereas the blue lines represent the lower order approximations with $g_3$ and $g_5$ set to zero.

Figure 4.

Comparison of asymptotic results (dashed lines) with exact results (black solid lines) when $\nu =0.1$ . The red dashed lines in (a) and (b) correspond to equations (89) and (86), respectively, and whereas the blue lines represent the lower order approximations with $g_3$ and $g_5$ set to zero.

5. Conclusion

Reduced plate and beam models have played an important role in engineering applications since they allow the most important information to be extracted without having to solve the fully 3D elasticity problem. In particular, the Euler–Bernoulli beam theory has frequently been used in studying pattern formation in film–substrate bilayers. Together with the Winkler assumption for the response of the substrate, this classical model is self-consistent and yields a leading order asymptotic expansion for the critical strain that agrees with that given by the exact 3D theory. In this paper, we have derived a refined model for the film–substrate interaction. The main motivation for deriving such a model is our planned study of pattern formation in a coated half-space where the coating has periodic but piecewise homogeneous material properties. For this problem, it is much harder to use the exact 3D theory to obtain analytical results than for the much studied case when the coating is homogeneous. Our preliminary study has shown that using the classical model is an attractive option, but its effectiveness is seriously hampered by the limited range of validity in which the theory is valid. It is hoped that the current refined model will increase the range of validity, and the relevant results will be published in a sequel to the current paper.

We used the methodology first proposed by Dai and Song [41] to derive the relationship between the traction and displacement vectors at the interface. The derivation consists of two steps. The first step is to express the traction and displacement vectors at the interface, $\mathit{t}(x_1,h)$ and $\mathit{u}(x_1,h)$ , in terms of the displacement vector at the traction-free surface, $\mathit{u}(x_1,0)$ . The second step is to treat the expression for $\mathit{u}(x_1,h)$ as an asymptotic expansion in $h$ and invert it to express $\mathit{u}(x_1,0)$ in terms of $\mathit{u}(x_1,h)$ . The most attractive feature of this methodology is that all manipulations only involve the solutions of linear algebraic equations, and the expansions can be carried out to any order in $h$ on a symbolic manipulation platform. Also, our treatment of both the coating and half-space can easily be generalised to deal with anisotropic and/or inhomogeneous materials [56,57], non-linear effects [58], liquid crystal elastomers [59,60], an infinite beam embedded in an infinite elastic solid [61], and the case when $U$ and $V$ depend on both of the in-plane coordinates.