## Abstract

This study investigates the buckling of a uni-axially compressed neo-Hookean thin film bonded to a neo-Hookean substrate. Previous studies have shown that the elastic bifurcation is supercritical if *r*≡*μ*_{f}/*μ*_{s}>1.74 and subcritical if *r*<1.74, where *μ*_{f} and *μ*_{s} are the shear moduli of the film and substrate, respectively. Moreover, existing numerical simulations of the fully nonlinear post-buckling behaviour have all been focused on the regime *r*>1.74. In this paper, we consider instead a subset of the regime *r*<1.74, namely when *r* is close to unity. Four near-critical regimes are considered. In particular, it is shown that, when *r*>1 and the stretch is greater than the critical stretch (the subcritical regime), there exists a localized solution that arises as the limit of modulated periodic solutions with increasingly longer and longer decaying tails. The evolution of each modulated periodic solution is followed as *r* is decreased, and it is found that there exists a critical value of *r* at which the deformation gradient develops a discontinuity and the solution becomes a static shock. The semi-analytical results presented could help future numerical simulations of the fully nonlinear post-buckling behaviour.

## 1. Introduction

There exists a large number of bifurcation phenomena in Nature, engineering and everyday life that are fully nonlinear and cannot be described by the traditional near-critical weakly nonlinear analysis. Examples include the cusp-shaped patterns that appear on the inner surface of an axially compressed hollow cylinder [1] or an everted hollow cylinder (A. Juel 2006, personal communication), kinks in a highly deformed solid cylinder or tube [2,3], wrinkling patterns in film/substrate bilayers (e.g. [4,5]), and creases in swollen cellular foams and gels (e.g. [6]). Because of their fully nonlinear nature, such phenomena are usually studied either numerically using finite-element packages or with the aid of simplified model equations (e.g. [7]). In this paper, we study one problem belonging to this class of elastic bifurcations using a semi-analytical approach. In particular, we aim at investigating the nonlinear behaviour of the Biot instability by taking the limit from a coated to a homogeneous elastic half-space under uni-axial compression. Taking this limit enables us to bring in all Fourier modes in a weakly nonlinear analysis, and thus to demonstrate the existence of patterns that can usually only be observed in the fully nonlinear regime.

Taking the limit from a coated to a homogeneous half-space can be viewed as one scheme to unfold the nonlinear characteristics of the Biot instability. Linear buckling analysis for a homogeneous half-space was first carried out by Biot in 1963, but its post-buckling behaviour has eluded full understanding for many decades. The difficulty lies in the fact that a homogeneous half-space does not have a natural length scale, and consequently there does not exist a distinguished mode number—all modes having equal status. A naive weakly nonlinear post-buckling analysis would be to write the solution as a Fourier integral or Fourier series, and then impose a solvability condition at the second order of successive approximations. This procedure would result in an infinite system of quadratic equations which do not seem to have convergent non-trivial solutions [8]. This indicates that, if post-buckling solutions existed, they should contain some form of discontinuities, e.g. in the form of static shocks. The next natural step is then to approach such non-smooth solutions by a limiting process. One such scheme is to assume that the surface of the half-space is corrugated and to follow the evolution of the surface profile as lateral compression is increased gradually. It was shown by Fu [9] that the evolution is indeed terminated by the formation of a static shock wave, and that the bifurcation is likely to be subcritical. A similar scheme was employed by Cao & Hutchinson [10] augmented also by a fully numerical simulation. Another scheme is to first assume that the half-space is coated by a thin layer with bending stiffness and then let the bending stiffness tend to zero. This scheme was carried out by Hohlfeld & Mahadevan [11], where analysis of a simplified model was combined with a fully numerical simulation.

The above search for a non-trivial post-buckling solution is closely associated with another challenging problem in continuum mechanics, namely the problem of the existence of surface acoustic waves of permanent form [12]. Since Biot's surface wrinkling mode is simply a standing surface wave with zero speed induced by pre-stress, it is not surprising that similar problems were encountered in the determination of nonlinear surface wave solutions. We refer to Fu & Hill [13] for a critical review of the relevant literature.

The problem of pattern formation on a coated half-space, or a thin film bonded to a substrate, has received much attention in recent years due to its potential application in a variety of situations (see, for example, Li *et al.* [14] for a comprehensive review). Given *μ*_{f} and *μ*_{s} as the shear moduli of the film and substrate, respectively, the post-bucking analysis by Cai & Fu [15] demonstrated that the nature of the elastic bifurcation is determined by the ratio *r*=*μ*_{f}/*μ*_{s}, being supercritical if *r*>1.74 and subcritical if *r*<1.74 (note, however, that the *r* in that paper corresponds to 1/*r* here). This fact was confirmed by a more recent study by Hutchinson [16], using a slightly different procedure, under a more general framework where a prestretch in the substrate before the film is attached is also considered. Seemingly unaware of this result, previous authors have always focused their attention on the regime *r*>1.74 (e.g. [6]). Thus, much remains unknown about the regime *r*<1.74, which is expected to be harder to simulate numerically due to sensitivity to imperfections. In this paper, we provide some semi-analytical solutions for the case when *r* is close to unity. The latter assumption enables us to incorporate the effects of all modes in a self-consistent manner, and to offer some insight on the post-buckling solution as the coated half-space reduces to a single homogeneous half-space. It is hoped that our (weakly nonlinear) results will provide a basis for future numerical simulations of the fully nonlinear post-buckling behaviour when *r*<1.74.

The rest of the paper is organized as follows. After formulating the buckling problem in the next section, we describe in §3 our asymptotic procedure and summarize the solution to the leading-order buckling problem. The nonlinear amplitude equations are derived in §4 by two alternative procedures; one of them is to express the total energy in terms of the Fourier amplitudes and then apply energy extremization. Numerical solutions are presented in §5, and we conclude in §6 with a summary of the main results.

## 2. Governing equations

Let us consider a general homogeneous elastic body *B* that is composed of a non-heat-conducting incompressible elastic material. We denote by *B*_{0} its initial unstressed state and by *B*_{e} a finitely deformed configuration. We choose a rectangular coordinate system relative to which the position vector of a material particle in *B*_{e} is denoted by ** x** with components (

*x*

_{i}). Suppose that

*B*

_{e}is further subjected to a small-amplitude static perturbation. The corresponding incremental displacement and pressure fields are denoted by

**and**

*u**p*, both of which are functions of

**. The governing equations for**

*x***and**

*u**p*consist of the incremental equilibrium equation

*χ*

_{ij}are components of the incremental stress tensor defined by

*B*

_{0}to

*B*

_{e}, and

*ϵ*is a small parameter characterizing the amplitude of (

*u*

_{i,j}) and

*p*. To simplify the analysis, we assume that both the primary and incremental deformations are plane-strain so that

*u*

_{3}≡0 and

*u*

_{1}and

*u*

_{2}are independent of

*x*

_{3}. Throughout this paper, the usual summation convention on repeated suffixes is observed and the range of summation is from 1 to 2.

In terms of the stress tensor ** χ**, the incremental surface traction

**on a material surface with outward unit normal**

*t***in**

*n**B*

_{e}is given by

*μ*is a shear modulus and

*I*

_{1}is the first principal invariant of the left Cauchy–Green strain tensor. We then have

**is the left Cauchy–Green strain tensor associated with the finite deformation from**

*B**B*

_{0}to

*B*

_{e}. For a neo-Hookean material, the second and higher order elastic moduli are all zero. Thus, a neo-Hookean material behaves as a linear material under incremental deformations, so that only geometrical nonlinearity is taken into account.

We now specialize the above general equations to the structure of a neo-Hookean elastic layer with shear modulus *μ*_{f} bonded to another neo-Hookean elastic half-space with shear modulus *μ*_{s}. We choose our coordinate system such that the layer and half-space occupy in *B*_{e} the regions −*h**≤*x*_{2}≤0 and *h** is a constant. The first-order elastic moduli are given by
*B*_{0} to *B*_{e} is a uni-axial compression. We then have ** B**=diag{λ

^{2}, λ

^{−2}}, where λ is the principal stretch along the

*x*

_{1}-direction. Applying the traction-free boundary condition on

*x*

_{2}=0 and the traction continuity condition on the interface

*x*

_{2}=−

*h**, we obtain

*B*

_{0}to

*B*

_{e}and the associated stress field are determined by a single parameter: the principal stretch λ. It is expected that, if λ reaches a certain critical value, an inhomogeneous solution may bifurcate from the above homogeneous state. The bifurcation value and the post-buckling states are determined by solving (2.1), (2.2) in

We non-dimensionalize all the governing equations and auxiliary conditions by using *μ*_{s} as the stress scale (for *L* as a typical length scale (for *x*_{i} and *u*_{i}), where *L* is a positive constant that will be specified later. For the sake of clarity we shall use the same letters to denote their non-dimensional counterparts except that the non-dimensional film thickness is denoted by *h*. Then the non-dimensional governing equations and auxiliary conditions are obtained from their dimensional counterparts by setting *μ*_{s}=1, *μ*_{f}=*r*, *h**=*h*, where *r*=*μ*_{f}/*μ*_{s} is the only non-dimensional material constant that characterizes the stiffness of the coating layer relative to that of the half-space.

## 3. Asymptotic expansions and the leading-order problem

Let us now consider the limit *ϵ* is the same small parameter as in (2.3) and *r*_{0} is an arbitrary *O*(1) constant. The deviation (*r*−1) is chosen to be *O*(*ϵ*) so that weakly dispersive and nonlinear effects appear at the same order. The case of a homogeneous half-space is recovered by taking the limit *r*_{0}→0.

We shall focus on the near-critical regime by writing
_{cr 0} and λ_{0} are constants. The critical stretch λ_{cr 0} will be determined shortly, whereas λ_{0} is to be specified. Note that we have reserved λ_{cr} to denote the mode-number-dependent critical stretch for an arbitrary coated half-space later. Correspondingly, the strain tensor ** B** has the expansion

*p*

^{(1)},

*p*

^{(2)},… are functions of

*x*

_{1}and

*x*

_{2}. On substituting (3.1)−(3.7) into the nonlinear eigenvalue problem and equating the coefficients of

*ϵ*, we obtain at leading order

*h*<

*x*

_{2}≤0. We note that, since the coating and half-space behave as the same material at leading order, the interfacial continuity conditions (2.9) are automatically satisfied. The leading-order problem (3.8)−(3.11) is then the same as the one for a single homogeneous half-space, which was first solved by Biot [18]. Next we summarize the main results which will be used later.

Equation (3.8*b*) implies the existence of a function *ψ*(*x*_{1},*x*_{2}), such that
*a*) and the boundary conditions can then be expressed entirely in terms of *ψ*. The resulting equations admit a single-mode surface-wave-type solution of the form
*H*(*x*_{2},*k*) is the shape function which can be determined by substituting (3.13) into the equilibrium equations (3.8)_{1}. It can be shown with the use of the decay condition (3.10) that, to within a multiplicative constant,
_{cr 0} must satisfy

Substituting (3.13) and (3.14) into (3.12), we obtain
*W*_{1} and *W*_{2} defined by
*s*_{1} and *s*_{2} are defined by *s*_{1}=1, *s*_{2}=λ^{2}_{cr 0} and are introduced in order to apply the summation rule.

Corresponding to (3.17), the pressure *p*^{(1)} must necessarily take the form
*a*) with *i*=1, we find
*F*(*s*_{2})≡0 in the present case, we choose not to use the simplified form in order to facilitate comparison with the more general cases.

We finally note that the above leading-order solution is valid for any mode number *k*.

## 4. Amplitude equations

For the leading-order problem which is ‘non-dispersive’, we look for a solution of the form
*A*_{m} are complex constants and the summations exclude *m*=0 (or equivalently we assume *A*_{0}=0). Without loss of generality, we have taken the length-scale *L* to be the inverse of the *dimensional* fundamental mode number, *k** say, so that the *non-dimensionalized* mode number of the fundamental mode in the above expression is unity. Correspondingly, we have *h*=*k***h**. Thus, for example, if the original (dimensional) coating thickness *h** were fixed, *h*→0 would correspond to the small mode number limit.

It can be checked that, if only a finite number of modes are included in (4.1), the second-order problem will become unsolvable. This is a unique feature of non-dispersive bifurcation or wave propagation problems.

To ensure that the expressions (4.1) for *p*^{(1)} are real, we impose the condition that

In principle, the evolution equations for the Fourier amplitudes *A*_{m}(*m*=1,2,…) can be determined from a solvability condition imposed on the governing equations for *p*^{(2)}; these governing equations are obtained by substituting (3.1)−(3.7) into the nonlinear eigenvalue problem (specified in §2) and equating the coefficients of *ϵ*^{2}. However, this approach is algebraically cumbersome. A more efficient method is the virtual work method proposed in Fu [19]; this method was used in Fu & Hill [13] to derive the evolution equations for nonlinear travelling waves in a coated elastic half-space having similar material properties as the one considered in this paper. Since we would need to assess the energy of different solutions anyway, as a by-product we may also derive the amplitude equations from the stationarity of the total energy. In the following, both approaches are employed and are used to provide a useful check on each other.

We start with the general (incremental) energy expression
*S*_{t} is the part of ∂*B*_{0} where traction *B*_{e} and then using the divergence theorem to eliminate the surface integral, we obtain
*ϵ*^{3} for other material models because of material nonlinearities.

We now specialize the above general energy expression to the coated elastic half-space with *B*_{e} replaced by
_{2} and (3.7) into the above expression, we obtain
*ϵ*^{2} term can be shown to vanish identically, whereas terms of order *ϵ*^{4} or higher have been neglected. We observe that in terms of unscaled coordinates the actual energy per unit wavelength is equal to the above *G* divided by *h*, a fact that will be used when we compare energies of different wavelengths when the coating thickness is fixed.

To evaluate (4.6) explicitly, we follow Fu & Hill [13] and write
*x*_{1} will survive the integration. We thus obtain
*ω*^{(1)}(*r*) and *ω*^{(2)}(*p*,*r*) are defined by
*G* with respect to *A*_{k} for any *k*, and then setting the resulting derivative to zero. Their reduced form has been checked against the amplitude equations derived by adapting the procedure used in Fu & Hill [13] and given by
*a*,*b*,*c*,*m*,*n* in (4.12) are carried out beforehand. It turns out that the resulting expression is not too long, as written in appendix A.

Neglecting the nonlinear terms in the amplitude equations (4.10), we would obtain *c*(*kh*)=0, which is the leading-order bifurcation condition. This condition gives a relation between *r*_{0} and λ_{0}, or equivalently a relation between *r*−1 and the bifurcation value λ_{cr}−λ_{cr 0} of λ−λ_{cr 0}. On the other hand, by taking the limit *r*→1 in the (exact) bifurcation condition of Cai & Fu [15], we may obtain

The asymptotic relation (4.13) together with (3.1) and (3.2) implies that the leading-order bifurcation condition must be given by λ_{0}=*a*_{1}(*kh*)*r*_{0}. We have checked numerically to verify that λ_{0}=*a*_{1}(*kh*)*r*_{0} with *a*_{1}(*kh*) given above does indeed satisfy *c*(*kh*)=0. Furthermore, since *c*(*kh*) must necessarily be a multiple of (λ_{0}−*a*_{1}(*kh*)*r*_{0}) and the first two terms in the expression for *c*(*kh*) are independent of λ_{0} and *r*_{0}, we may also write *c*(*kh*) in the alternative form
*k*>0 this expression can be simplified to
*a*_{1}(*kh*) defined by (4.14); the latter is monotone increasing for 0<*kh*<*k*_{cr}, and monotone decreasing for *kh*>*k*_{cr}, with the maximum *a*_{1 max}(≈0.2418) attained at *kh*=*k*_{cr}, where *k*_{cr}≈1.73. Figure 1 compares the exact bifurcation condition and the two-term asymptotic expression given by (4.13) for the cases when *r*=1.05 and *r*=0.95, respectively. It shows that the asymptotic expression provides a uniformly valid approximation for the exact bifurcation condition over the entire range of *kh*.

## 5. Post-buckling solutions

Let us now investigate the behaviour of post-buckling solutions. Since the kernel *D*_{−k}=−*D*_{k}, and the amplitude equations (4.10) then reduce to the following system of real equations:
*ω*^{(1)}(*k*)=−λ_{0}*kc*(*kh*). The Hessian matrix of G, formed from the second-order derivatives above, will be used to assess the stability of the solutions obtained. A solution is said to be unstable if this matrix has at least one negative eigenvalue, and be stable if all the eigenvalues are positive (e.g. [21]).

We also observe that the kernel *a*. Thus, if we obtain a solution such that the only non-zero Fourier coefficients are *D*_{K},*D*_{2K},*D*_{3K},… with *K* a fixed integer, then the amplitude equations (5.1) may be replaced by
*m*=1,2,…. It then follows that, when *h* is replaced by *Kh*, there is a solution in which the *m*th Fourier coefficient is given by *KD*_{Km} (and, in particular, the first Fourier coefficient is non-zero). In fact, it can be shown that these two solutions correspond to the same solution in terms of the unscaled variables, and so in our interpretation of numerical results solutions with *D*_{1}=*D*_{2}=⋯=*D*_{K−1}=0 for some integer *K*>1 are not counted as independent solutions.

### (a) Numerical procedure for solving the amplitude equations

We first observe that the infinite system (5.1) depends on *r*_{0} and λ_{0} through the combination *r*_{0}/λ_{0}. We shall therefore take *r*_{0} to be ±1 and only vary λ_{0} in our numerical calculations, bearing in mind that increasing λ_{0}, for instance, is equivalent to decreasing *r*_{0}. Once a solution for *D*_{k} is found, various quantities can be evaluated. For instance, the vertical displacement *h** in *B*_{e}, instead of 1/*k**, as the length scale. This is adopted in the rest of this paper, and we then have
*h*=*k***h** denotes the non-dimensional fundamental mode number.

We first replace the infinite system (5.1) by the truncated system
*M* is a suitably chosen positive integer. This truncated system is augmented with the assumption
*M* quadratic equations with *M* unknowns.

We start with *M*=2, in which case the two quadratic equations can be solved explicitly to give
*M* in unit steps. At each step, *M*=*N*+1 say, the solution from the previous step together with *D*_{N+1}=0 is used as the initial guess in the solution of *N*+1 quadratic equations. If the convergence criterion
^{−4} and 10^{−6}. If the above criterion is not satisfied for *N* as large as 300, we consider that the procedure does not produce a valid solution. The procedure usually has one of four outcomes: (i) the trivial solution, (ii) a nicely converged non-trivial solution, (iii) a situation where convergence seems possible but an infinite number of modes would be needed (the case of a static shock), and (iv) a situation where convergence is impossible because the Fourier coefficients do not decay at all. This procedure is validated using a simple problem in appendix B.

We may also start with *M*=3, in which case the three simultaneous equations can be reduced to

The coefficient of _{1} is found to be a positive constant. It then follows that the associated solutions are real only if *c*(*h*)*c*(2*h*)≥0, that is,
_{3} are real if *c*(3*h*)[*c*(3*h*)+5.02969*c*(*h*)]≥0, that is,

### (b) Numerical results

In view of the fact that the neutral case corresponds to λ_{0}=*r*_{0}*a*_{1}(*h*) and the maximum of *r*_{0}*a*_{1}(*h*) is *r*_{0}*a*_{1 max} when *r*_{0}>0 and zero when *r*_{0}<0 (figure 2), we shall discuss the following four different cases separately:

*Case* (*i*): *r*_{0}=1, λ_{0}>*r*_{0}*a*_{1 max}.

This is the subcritical regime when the coating is slightly stiffer than the substrate, and there exists a critical mode number *k*_{cr}=1.73. In this case, it is found that the two starting solutions given by (5.8) lead to two solutions which are equivalent in the sense that one is simply a translation of the other (by *π*) in the *x*_{1}-direction. We have also tried starting with *M*=3. It is found that two of the four solutions given by (5.10) do not lead to a convergent solution, whereas the other two lead to the same solution as when *M*=2.

Recalling that *h* denotes the non-dimensional mode number of the fundamental mode, we first consider the case when *h*=*k*_{cr}, which is usually the concern of a weakly nonlinear analysis. An unmodulated sinusoidal solution is found for each value of λ_{0} less than a certain cut-off value (to be defined shortly). The dashed line in figure 3*a* shows how the amplitude of surface elevation varies with respect to 1/λ_{0} (we have used 1/λ_{0} because it increases with the load), where here and hereafter we use *u*_{2} to denote _{0} is gradually increased, more and more modes are needed to obtain a convergent solution, resulting in a steeper and steeper surface profile. The evolution is terminated with the formation of a static shock in which *u*_{2,1}(*x*_{1},0) becomes discontinuous (figure 4). To show this behaviour more clearly, we have shown in figure 3*b* the variation of |*u*_{2,11}(0,0)| with respect to 1/λ_{0}. It is found that this gradient tends to infinity as 1/λ_{0} approaches a cut-off value approximately equal to 3.82. Secondly, although as 1/λ_{0} is decreased higher modes contribute to the steepening of the surface elevation profile, they have very little effect on the amplitude |*u*_{2}(0,0)|: this quantity when calculated using the three-mode approximation (5.10) is almost indistinguishable from the dashed curve presented in figure 3*a*.

For values of *h* in a small neighbourhood of *k*_{cr}, other periodic solutions are also found. Two such solutions corresponding to *h*/*k*_{cr}=1.15,0.8 are shown in figure 5 (the top two plots). To demonstrate its unmodulated periodic nature, we have shown each profile over the interval [−*π*/*h*,3*π*/*h*], which is twice the period. The rest of the plots in this figure show how the periodic solution changes character as *h* is reduced further: it gradually evolves into a modulated solution in which the fast oscillation seems to have a mode number equal to *k*_{cr}, whereas the slow modulation has a mode number equal to *h*. With regard to the two profiles corresponding to *h*/*k*_{cr}=0.2,0.1, which typify those obtained with other small values of *h*, it is found that the central part of each profile over one period, e.g. the part consisting of the two main peaks and the single trough at *x*_{1}=0, is almost identical to the central part of the other (i.e. it is almost independent of the value of *h*). The main difference between the two profiles is that the oscillatory tail is much longer when *h*/*k*_{cr} is smaller. In the limit *h*/*k*_{cr}→0, the tail becomes infinitely long and the modulated periodic solution becomes a fully localized solution. Therefore, as far as the calculation of *u*_{2}(0,0) and *u*_{2,11}(0,0) for the localized solution is concerned, we may use the profile obtained with any *h* that is small enough. It is found, however, that the maximum *h* that can be used to serve this purpose depends on the value of λ_{0}: it decreases as λ_{0}−*a*_{1 max}*r*_{0} decreases. For instance, taking _{0}−*a*_{1 max}*r*_{0} is greater than 0.0025, but it only returns an unmodulated periodic solution when λ_{0}−*a*_{1 max}*r*_{0} becomes less than 0.0025. This is why the solid curves in figure 3 are obtained by taking two different values of *a*, demonstrating the fact that the two solutions are not exactly identical, but both provide a good approximation for the localized solution.

In trying to understand what solutions are returned when the numerical procedure does not return a modulated periodic solution, we consider small values of λ_{0}−*a*_{1 max}*r*_{0} and find the following ‘mode-locking’ behaviour when *h* is not small enough: the *D*_{j}'s are non-zero only for *j*≥*K* and *j* an integer multiple of *K*, where the integer *K* depends on the value of *h*. For instance, when λ_{0}−*a*_{1 max}*r*_{0}=0.001, we have the following correspondence:
*Kh* always lies between 1 and 2, the two integers bracketing the critical value 1.73 of *h*. In view of the comments made below (5.5), these solutions are the same solutions as when *h* is replaced by *hK* and the fundamental mode is non-zero. Thus, when λ_{0} is sufficiently close to its critical value *a*_{1 max}*r*_{0}, the fundamental mode selected in the solutions is always the mode nearest to the critical mode. The associated surface profile in each case is unmodulated and convergence is achieved with only a few modes. The mode-locking exhibited in (5.13) stops when *h* is approximately smaller than _{0} approaches *a*_{1 max}*r*_{0}, the modulation takes place over an increasingly greater and greater length scale and this can only be achieved by decreasing the fundamental mode number.

As we may expect, the longer the tail of the modulated solution is, the smaller the averaged energy over one period is. This is also verified numerically using expression (5.2). Thus, the limiting localized solution has the smallest energy among all modulated periodic solutions. For each modulated periodic solution obtained, the associated Hessian matrix of the energy function *G* has exactly one negative eigenvalue. If we order the eigenvalues *ω*_{i}(*i*=1,2,…) such that |*ω*_{1}|<|*ω*_{2}|<|*ω*_{3}|<⋯ and denote the negative eigenvalue by *ω*_{N}, then the index *N* increases as *h* decreases. For instance, for *N* is equal to 1,2,3,8,18 and 27, respectively, and with the tolerance in (5.9) set to be 10^{−4} convergence of solution is achieved with 12,26,34,72,128 and 180 modes, respectively. It has also been verified that the index *N* is independent of the tolerance set and hence of the number of modes included.

As in the case of periodic solutions with *h*=*k*_{cr}, the surface elevation profile also steepens up as λ_{0} increases (figure 6). In this case, the value of *u*_{2,11}(0,0) tends to infinity when 1/λ_{0} tends to 3.89, slightly higher than the corresponding value for periodic solutions (figure 3*b*).

We finally observe that the fully localized solution is similar to the localized solution predicted by the simple model of a linearly elastic beam on a nonlinearly elastic softening foundation (e.g. [22,23]). In fact we believe that the localized solution determined here can be matched with the localized solution determined by Cai & Fu [15] under appropriate limits.

*Case* (*ii*): *r*_{0}=1, λ_{0}<*a*_{1 max}*r*_{0}.

In this case, both conditions (5.11) and (5.12) are satisfied only for values of *h* in some intervals (figure 2). For instance, when λ_{0}−*a*_{1 max}*r*_{0}=−0.001, these intervals are
*K* has the same meaning as in (5.13) and *K*=0 signifies the fact that only the zero solution is obtained. All of the non-zero solutions are unmodulated periodic solutions, and, in view of the comments made below (5.5), these are the same solutions as the ones when *h* is replaced by the product *Kh*. Thus, only solutions with mode numbers close to *k*_{cr} are selected. In particular, no non-trivial solutions are obtained when *h* is larger than 2.5, and in contrast with the previous case no modulated periodic solutions are obtained no matter how small the value of *h* is.

*Case* (*iii*): *r*_{0}=−1, λ_{0}>0.

We first observe that when *r*_{0}=−1 the bifurcation value of λ_{0}, namely *r*_{0}*a*_{1}(*h*), approaches zero from below as _{0} and *h* are small enough but modulated periodic solutions are never found. For each unmodulated periodic solution found, the associated Hessian matrix of *G* has one negative eigenvalue, and so the solution is unstable. In figure 7, we have shown a typical solution corresponding to *r*_{0}=−1, _{0}=0.001,0.005,0.01,0.012. It is observed that, as λ_{0} is increased, the solution has larger and larger amplitude but we believe that it does not tend to a static shock based on two reasons. Firstly, the gradient of surface elevation does increase but it does not change significantly as λ_{0} is increased. Secondly, as λ_{0} becomes large enough, it becomes impossible to obtain a convergent solution since the solution becomes increasingly more and more sensitive to the tolerance imposed and the truncation number. In figure 8, we have shown the Fourier amplitude *D*_{k} against *k* when three typical truncation numbers are used. It is found that, although *D*_{k} for large values of *k* are sufficiently small to have a negligible effect on the profile of *u*_{2}(*x*_{1},0), they do not decay as *k* increases and make *u*_{2,1}(*x*_{1},0) and *u*_{2,11}(*x*_{1},0) diverge. This divergent behaviour becomes more and more pronounced (through increase in the amplitudes of *D*_{k} for large *k*) as λ_{0} or *h* is increased. The amplitudes of *D*_{k} for large *k* may even grow with respect to *k* when λ_{0} or *h* is large enough. Reliable convergent solutions could be obtained only for *h* approximately smaller than _{0} is. Thus, in this case only (unmodulated) periodic solutions with wavelengths much larger than the coating thickness are obtained.

*Case* (*iv*): *r*_{0}=−1, λ_{0}<0.

Finally, we consider the case when *r*_{0}<0 and λ_{0} is reduced from zero. As in case (ii) both conditions (5.11) and (5.12) are satisfied only for values of *h* in some intervals (figure 2). For instance, when λ_{0}=−0.001, these intervals are (0,0.0198) and (0.0396,4.0660). No solutions are found for values of *h* in the first interval. For values of *h* in the second interval, unmodulated periodic solutions are only found for *h* between 0.0396 and 0.14, but modulated periodic solutions are never found. For *h*>0.14 and no matter how small |λ_{0}| is, the procedure produces either a trivial solution or a divergent solution that has the same behaviour as in figure 8. For each choice of *h* for which an unmodulated periodic solution could be found when λ_{0} is reduced from zero, the solution can be viewed as a continuation of the solution found in the previous case at the same value of *h*. The solution can be continued all the way to when the equality in (5.11) holds, in which case the solution has zero amplitude.

## 6. Conclusion

In this paper, we have investigated the post-buckling behaviour of a coated elastic half-space in which the coating and substrate have almost identical properties. This parameter regime has previously not been examined in the literature, either numerically or analytically. A theoretical analysis is much desired since this problem would be very hard to study by only using numerical tools, such as finite elements. In fact, since the bifurcation in this parameter regime is subcritical, the post-buckling behaviour would be sensitive to imperfections. Such imperfection sensitivity would in turn introduce major issues in the characterization and continuation of post-buckling solutions. Thus, our semi-analytical results and the phase map of the solution behaviour over the entire parameter space would provide a useful guide in any future analytical or numerical investigations.

When the coating and half-space have almost identical material properties, the coated half-space becomes ‘almost non-dispersive’ in the sense that all the modes are near-critical so that solvability conditions need to be imposed at second order of a successive asymptotic analysis, while for ‘dispersive’ problems involving a single critical mode the solvability is imposed at third order to obtain the amplitude/evolution equation. We have adopted the approach of using Fourier series expansion of the incremental displacement and pressure fields. The mode number of the fundamental mode is not fixed to be the critical mode number and is allowed to be arbitrary. If the fundamental mode number were small and all subharmonics (relative to the critical mode) were found to be non-zero, we would obtain modulated or localized solutions. Unmodulated periodic solutions would be obtained if all the subharmonics were found to be zero. Both situations did arise in our numerical calculations.

Our numerical calculations divide the near-critical parameter space into four regimes. The most interesting regime seems to be where the coating is slightly stiffer and the uni-axial stretch is slightly above its critical value, which corresponds to case (i) in §5b. In this case, there is a distinguished critical mode given by *kh*=1.73, although all the other modes can interact resonantly. This is also the regime where our numerical calculation met fewest difficulties because each solution can be followed in a well-understood manner, and divergence or convergence is always clear-cut. As a consequence, we have been able to understand this regime completely. Two main results may be highlighted. Firstly, among all solutions possible, the localized solution always has the least energy. Secondly, under the assumptions made, our solution for the displacement field is of the form
_{cr 0} and the Fourier coefficients depend on the loading parameter λ and material constant *r* through the combination *r*_{0}/λ_{0}. It is found that as *r*_{0}/λ_{0} approaches 3.82, or equivalently,

We further note that the localized solution reported is the limit of modulated periodic solutions with increasingly longer tails. This limit is achieved when the fundamental mode number tends to zero. For each such nearly localized solution, which is typically formed from 200 modes, the total energy *G* is a cubic function of the Fourier amplitudes, and its Hessian matrix always has a single negative eigenvalue, which indicates that the solution is unstable. However, it is observed that, as *h* is decreased, the negative eigenvalue is associated with an increasingly higher and higher mode. It seems that, in the limit when the fundamental mode number tends to zero, the single negative eigenvalue would be associated with a mode with an infinite mode number. In the present context, we have not been able to deduce stability or instability for the localized solutions. However, it might be possible to draw an analogy with perhaps a similar situation, namely the bifurcation of a pressurized long membrane tube into a localized bulge (e.g. [24]). In that case although the localized bulging configuration is unstable in its early stage of development, it does eventually lead to a stable bulging state, and analysis of the intermediate (unstable) bulging configurations is an essential part in the understanding of the whole inflation process.

## Data accessibility

This manuscript does not contain primary data and as a result has no supporting material associated with the results presented.

## Funding statement

The work of Y.B.F. was supported by the National Natural Science Foundation of China (grant no. 11372212), and by a Visiting Professorship of the Universit Pierre et Marie Curie, Paris 6. The work of P.C. was partially supported by the INSERM grant AAP PhysiCancer, and by the ‘Start-up Packages and PhD Program’ project, co-funded by Regione Lombardia through the ‘Fondo per lo sviluppo e la coesione 2007–2013’ (formerly FAS program).

## Author contributions

The two authors contributed equally in developing the mathematical model and method of analysis. Y.B.F. performed the calculations and wrote the paper. Both authors gave final approval for publication.

## Conflict of interests

The authors have no competing interests.

## Appendix A. Simplified expression for the kernel K ( p , q ) given by (4.12)

## Appendix B. Validation of the numerical procedure

The procedure proposed in §5a is validated by considering the simple equation

If we expand *u* as a Fourier series
*a* is a constant characterizing the mode number of the fundamental mode, then the Fourier coefficients satisfy the following system of quadratic equations:
*a*→0, the procedure reproduces the exact solution (B 2) and no other spurious solutions.

The numerical scheme proposed above has also previously been validated by considering the more complicated equation *u*_{xxxx}+*Pu*_{xx}+*u*−*u*^{2}=0, which has a richer family of solutions; this equation describes the deflection of an elastic beam on a softening foundation that is subjected to a compressive force *P* (e.g. [25]).

- Received December 18, 2014.
- Accepted April 1, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.