## Abstract

A thin flat rectangular plate supported on its edges and subjected to in-plane loading exhibits stable post-buckling behaviour. However, the introduction of a nonlinear (softening) elastic foundation may cause the response to become unstable. Here, the post-buckling of such a structure is investigated and several important phenomena are identified, including the transition of patterns from stripes to spots and back again. The interaction between these forms is of importance for understanding the possible post-buckling behaviours of this structural system. In addition, both periodic and some localized responses are found to exist as the dimensions of the plate are increased and this becomes relevant when the characteristic wavelengths of the buckle pattern are small compared with the size of the plate. Potential applications of the model range from macroscopic industrial manufacturing of structural elements to the understanding of micro- and nanoscale deformations in materials.

## 1. Introduction

Instability in elastic bodies is of practical importance to engineers, who are responsible for the design of safe structures. If our sole concern was to understand the initial stages of buckling of such bodies then a suitable model consisting of linearized equations would serve the purpose. On the other hand, a deeper appreciation of the true nature of the buckling together with knowledge of the properties of later stages of the failure mechanisms requires that nonlinear effects be included as well. In the overwhelming majority of post-buckling scenarios exact analytical solutions do not exist and so, necessarily, we have to resort to approximation, be it fully numerical or a simplification of the full system followed by analysis in various limited regimes [1]. Both of these strategies can reveal behaviour which is of interest not only to engineers but also to applied scientists and mathematicians. Instability phenomena are becoming of increasing importance in many fields [2] and some previous studies are being given renewed prominence owing to the discovery of physical phenomena down to the atomic level which these theories seem to model well [3]. It is in this spirit that here we study a problem that at first sight seems to be a simplistic two-dimensional extension of an exhaustively studied one, namely the buckling of a strut resting on an elastic foundation.

In structural engineering, plates are used as load-bearing elements in various guises, including floor slabs, retaining walls, pad foundations, ship hulls and aircraft wings. Of these, perhaps the pad would most naturally be modelled by a plate on a foundation. There are examples where certain parts of structural components can be modelled as plates which are supported; an obvious candidate is a sandwich panel for which a core supports a pair of thin skins. Several of these problems have been examined assuming one-dimensional behaviour and our intention is to extend these ideas to incorporate a second dimension. Moreover, the importance of supported plates is not limited to large-scale structures for similar systems have been proposed as models to look at thin elastic films supported by a substrate. There has been much work on one-dimensional *wrinkling* (another form of buckling) either when a specimen is compressed or if the supporting substrate (foundation) is pre-stretched when the elastic layer is placed on it and is then released. The plate equations for such a configuration have been solved using a finite-element technique and compared with experimental data [4]. Moreover, rather than just being an unfortunate by-product of loading, wrinkling in some physical systems may actually be desirable and its control can be effected by judicious application of load in various directions [5].

At an even smaller scale, studies have shown that wrinkling of nanoscale materials can be sufficient to modify electrical, optical and other physical properties. In particular, compressed graphene behaves in this way and a useful model for this material is that of an isolated elastic layer [6]. We believe that improved understanding of wrinkling in nano-materials might be gained by incorporating a supporting elastic medium as an integral part of the system.

The well-studied one-dimensional strut model has numerous applications from the buckling of submarine pipelines and railway tracks through to sandwich panels and as a phenomenological model for various shell structures [7–9]. Various mathematical techniques have been used to study this system including Rayleigh–Ritz and Galerkin analyses [10], multiple-scale perturbation expansions [11,12] and advanced hyperasymptotic methods [13,14]. Part of the motivation for this flurry of activity has been the fact that the strut model can be regarded as a member of a family of problems governed by a Swift–Hohenberg (SH) equation. The SH equation is a partial differential equation in both time and space that simplifies to an ordinary differential equation (ODE) form when interest is restricted to time-invariant behaviour in one spatial dimension. Static equilibria and their stability can readily be deduced from this reduced problem. A key distinction between the strut and plate models is that the former is of unbounded length—suitable for the study of very long structures—and solutions sought can be periodic or homoclinic (localized in extent) where deflection approaches zero in both positive and negative directions [15]. In reality of course, at least at the macroscopic scale, plates are finite and the choice of boundary conditions potentially has a greater effect on the response of the system to load than if they are modelled as being completely remote from the location of the buckle pattern. Nevertheless, if situations arise for which the deflection is restricted to a small portion of the structure then we could say such buckle patterns are effectively localized in an engineering sense although not homoclinic according to a strict mathematical definition.

In order to make headway with our stated goal of introducing a second spatial dimension, it may seem natural to simply investigate the appropriate planar version of the SH equation; see, for example, [16,17] for several investigations into two-dimensional localized patterns possible with the planar SH. However, extending in a simple-minded manner hides the complication that a flat plate may respond in other ways when loaded. For instance, there is now the possibility of bending about two axes and twisting which gives the plate extra stiffness and hence offers a good solution for efficient design [18]. Furthermore, compatibility and the constitutive equations arising from the theory of elasticity need to be satisfied and this adds extra terms together with the need to introduce a new stress function which relates deformations and strains to stresses. For moderate to large deflections, the bending of a flat plate is often modelled by the Föppl–von Kármán (FvK) equations for, although it is known that they have some limitations and inconsistencies, they have been shown to usefully model post-buckling plate behaviour at a qualitative level that is normally accurate enough for most engineering purposes [19].

We aim to study here the case where the behaviour of an elastic plate is modified owing to the presence of a supporting elastic foundation or substrate on which it rests. Depending on the detailed properties of the foundation different phenomena can be realized. We will discuss the various different buckling modes that are possible and how they relate to each other as the loading on the plate is varied. In order to demonstrate the potential usefulness of the plate model, we will mention a few examples from various fields that might benefit from our findings.

An important feature of the plate model is that it necessarily involves separate loading parameters from two directions. For the one-dimensional strut, any loading acts at a point and so cannot vary spatially, but in two directions this constraint has to be relaxed; the loading in either direction need not be constant and, in general, would vary with position. Moreover, proper consideration has to be taken of the aspect ratio of the plate. One might naively expect that a long thin plate will respond to loading in a way similar to that of a strut, but if the two dimensions of the plate are comparable, then it is far from clear what one might expect to happen. It would not be surprising if in this case distinct phenomena occur that are peculiar to the strictly two-dimensional case.

The remainder of the paper is organized as follows. In §2, the model is formally introduced and in §3 a Galerkin method is used to obtain approximate buckling solutions which emerge from the flat fundamental state. These Galerkin predictions suggest how to conduct efficient full numerical solutions. A dynamical stability analysis described in §4 establishes criteria that predict when a one-dimensional (striped) solution may be unstable to transverse instability leading to a fully two-dimensional pattern. Fully numerical solutions are obtained and described in §5 and these are compared directly with the earlier results. The paper closes in §6 with some conclusions of the work and a look forward to possible further investigations of the problem.

## 2. Model formulation

The FvK equations for plates [18] have been shown to provide an accurate model for post-buckling behaviour and an analysis demonstrates that the naturally emerging post-buckling solutions are periodic in two orthogonal directions with the formation of a rectangular cellular pattern of undulations [20].

We consider a thin flat plate made of a homogeneous and isotropic material of Young’s modulus *E*, Poisson’s ratio *ν* and with thickness *h* as shown in figure 1. The coordinate axes are aligned so that the mid-surface of the plate sits in the *a*and
*b*in which *D*=*Eh*^{3}/(12(1−*ν*^{2})) denotes the bending stiffness of the plate,

The crucial difference between our model and that for a plate supported only on its edges is the presence of the foundation resistance *w* and no shear interaction is present; then

Let us take the plate to be rectangular of dimensions *a*×*b* with thickness *h*≪*a*,*b*. Appropriate physical boundary conditions can then be of various types. The edges may be unsupported (free), simply supported (which prevents deflection but allows rotation) or fully fixed. Each one of these conditions supplies two kinematic conditions on the deflection or its derivatives. We will study the case of clamped edges in the *a* and width *b*
*symmetric-section* conditions used when seeking localized solutions to the one-dimensional strut model [12]. We will consider plates of various length in the

Let us pause here to consider the choice of the nonlinear foundation,

### (a) Non-dimensionalization

In order to analyse the system, it is helpful to scale variables to remove as many coefficients as possible. Therefore, we set
^{4} is the biharmonic operator with respect to the scaled variables *x* and *y* and *F*(*w*)=*w*−*k*_{2}*w*^{2}+*k*_{3}*w*^{3}. In this scaled state, we suppose that the plate has length *L* and width *γ*.

#### (i) Loading parameters and some simplifications

In the general nonlinear plate problem, the loading parameters in the *x*- and *y*-directions, *P*_{x} and *P*_{y}, should be treated as independent. There is no necessity for either to be constant in space and such variation is likely to be of interest in provoking other behaviours such as bending about a lateral axis. We defer this extra degree of freedom as a refinement to be dealt with in a future study and instead concentrate solely on the case where both *P*_{x} and *P*_{y} are constant. That said, we will allow the possibility of either load being compressive (positive) or tensile (negative). It is conceptually easy to envisage different loading in the two directions: the simplest of course being uniaxial compression. Furthermore, if the edges are prevented from moving in the plane by a pin or clamped support, then tensile stresses can be built up as the plate deforms in the *z*-direction as moderate to large buckling deformation increases. As outlined by Yang *et al.* [24], one-dimensional wrinkles can be induced in nano-materials by pre-stretching the substrate and then depositing the material of interest on it. When the tensile force is released this results in a compression acting in one direction. Forces can then be applied in an orthogonal direction to induce further changes to the material.

In this study, we will limit investigation to the case when, say, *P*_{y} is held constant and *P*_{x} is parametrically varied. Physically, this is equivalent to preloading the plate in one direction (not necessarily up to the buckling capacity in that direction) and then applying another load in the orthogonal direction.

### (b) Total potential energy and Euler–Lagrange equations

We shall use the energy of the system as a measure of the various buckled states. The total potential energy associated with the non-dimensional form of the von Kármán equations in a slightly compacted form is [20]
*P*_{x} and *P*_{y} and we identify the end-shortening of the plate in the *x*- and *y*-directions to be

Of course, the governing equations are obtainable from *V* by obtaining the associated Euler–Lagrange equations of the system but the energy itself plays an important role in the preferential formation of certain buckle patterns.

### (c) Linearization

To understand the critical behaviour of the system, it is necessary to linearize (2.6) and (2.7) about the trivial state. We will first look at the critical instability of an infinite plate so that boundary conditions can be safely ignored. In order to simplify the subsequent notation let us put *P*_{x}=*P* and *P*_{y}=*αP* where the parameter *α* can theoretically be any real value. However, without loss of generality, if we consider a plate that is longer in the *x*-direction than in the *y* then buckling will be dictated by *P*_{x} and so we can legitimately restrict our interest to |*α*|≤1. If we then seek Fourier mode solutions
*S* becomes negative and this then corresponds to *x*-spatially periodic modes. When *P*<*P*^{C} the eigenvalues are complex-valued with the critical values *P*^{C}=2 when *α*=1 and *α*=0; in contrast when *P*>*P*^{C} two of the eigenvalues are complex but the other two are real. Of course we should also linearize equation (2.7) but as the result is independent of *P* it has no bearing on the bifurcation behaviour of the system with respect to load.

## 3. A Galerkin analysis

To obtain numerical solutions and then investigate their evolution under varying parameters, it is necessary to derive an initial approximate two-dimensional solution to the system. Once a successful solution is obtained, we can employ numerical continuation to change *P*_{y} to the desired value.

We use a straightforward Galerkin procedure [25] for this purpose and start with a purely periodic response of the form
*A* is a constant amplitude and *ψ* is an appropriate shape function. The initial stress function, *ϕ*, can be approximated either as a multiple of *w* or as zero if we are seeking to start from close to the bifurcation from the trivial state. We set the length of the plate in the *x*-direction to be *L*=*π* and impose Neumann boundary conditions in that direction. If the width *γ*=2, then for clamped boundary conditions in the *y*-direction, at *y*=±1, we need to identify an appropriate set of functions that satisfy *w*(*x*,±1)=*w*_{y}(*x*,±1)=0. In the *x*-direction we assume Neumann boundary conditions, i.e. *w*_{x}(0,*y*)=*w*_{x}(*π*,*y*)=0 and, furthermore, if *w*(0,*y*)=*w*(*π*,*y*) then we have periodic boundaries in that direction. A suitable family of even functions obeying these boundary conditions is given by
*m* and *n* are integers which characterize the number of half waves occurring in the *x*- and *y*-directions, respectively, and a periodic boundary condition in the *x*-direction requires *m* to be even. A counterpart odd set of functions in the *y*-direction is given by
*x*-direction, *m* should be even. If we denote the left-hand side of first equation (2.6) by *Q*(*x*,*y*) and feed in the assumed form for *w* in (3.1), then the Galerkin formulation tells us that
*A* with a non-trivial even solution
*k*_{3} but the absence of the quadratic coefficient *k*_{2}. This feature is a consequence of taking a single term for the shape of *w* and leads to a symmetric form for *A*. For greater accuracy, we could use more terms, but, since the motivation for this Galerkin method is simply to supply some initial functions that can be used to generate full numerical solutions close to the trivial one, there is little to be gained by deriving a more accurate initial state. For real values of *A*, the term inside the square root must be greater than zero and this gives us a critical value for *P*_{y} is held constant.

## 4. Transverse instability of quasi-one-dimensional states

We next examine the possibility of a transverse instability of quasi-one-dimensional states (by which we mean one-dimensional states that are trivially extended to two spatial dimensions). We will subsequently find that the transverse instability of quasi-one-dimensional states plays a pivotal role for the emergence of non-trivial two-dimensional states.

We start by seeking solutions that vary only in the transverse direction *y*. If *W*(*y*,*t*) is a one-dimensional solution of (2.6) and if *Φ*(*y*) satisfies (2.7) then it is immediately apparent that
*x*, the end-shortening component *Δ*_{x}=0 and the stress function *Φ*≡0 everywhere. It is well known that stationary solutions of (4.1), *W*(*y*,*t*)=*W*(*y*), exist; see, for instance, [8]. The linear stability properties of the stationary *W*(*y*) with respect to perturbations in *y* can be inferred directly from results of the ubiquitous SH equation (given by (4.1) save that the *W*_{tt} is replaced by *W*_{t}). In particular, if a periodic orbit is linearly stable in the standard SH equation, then it is neutrally stable in (4.1). Hence, there do exist periodic orbits that are linearly stable with respect to perturbations in the *y*-direction.

We next consider a small transverse perturbation in the *x*-direction by letting *x*-direction, so that

We solve the spectral problem (4.6) for *η*=*η*(*β*). For small *β*, we write

At

We now have a procedure to predict the dynamical stability of patterns, at least when they start as functions of *y* alone. We therefore need to obtain solutions, *W*(*y*), of the ODE (4.1), subject to appropriate boundary conditions (clamped in this case). Solutions can readily be obtained using Auto [27] and the quantities *N*_{1} and *N*_{2} used to evaluate the transverse instability criterion (4.12). We compute the first four modes of (4.1) on a finite domain in *y* satisfying the boundary conditions (2.3) and the results are shown in figure 2.

We show in figure 3 a two-parameter plot of the transverse (in)stability region when *k*_{2}=2 and *k*_{3}=0.5. In figure 3, the values of *N*_{1} and *N*_{2} (see (4.11)) have been calculated numerically for a one-dimensional pattern and the boundary where the value of *P*_{x} above which the solutions become unstable in time. We see that the location of the lower boundary for *P*_{y} is almost independent of *P*_{x}; indeed, when *P*_{x} is pushed beyond some critical value, after which it becomes unstable and the pattern evolves to some other state.

A word of caution is appropriate. Our analysis described here is strictly only valid for a plate of infinite extent in the *x*-direction. Thus, it is unlikely to provide a complete picture of the stability of plates of finite length for which the choice of boundary conditions ought to be relevant. That said, the analysis should provide some clue as to the types of behaviours that may be seen when a striped solution of a finite plate (with deflection independent of *x*) encounters a perturbation which varies in the transverse planar direction.

## 5. Numerical solutions

The above approximations allow us to speculate on the form of some of the possibilities of solutions of the plate-on-foundation model. However, given the experience gained with related problems, such as the (one-dimensional) strut model [13], it is highly probable that there is an intricate solution structure which, initially at least, can best be explored numerically. This was done by discretizing the boundary value problem using a standard fourth-order scheme in the *x*-direction and then imposing Neumann boundary conditions using ghost points [28]. A pseudo-spectral Chebyshev discretization was employed in the *y*-direction with the clamped boundary conditions enforced by choosing the polynomial interpolant, *p*(*y*), of a function in *y* to be *p*(*y*)=(1−*y*^{2})*q*(*y*), where *q*(*y*) is a polynomial with *q*(±1)=0 ([29], ch. 14). The two-dimensional differentiation operators were constructed using Kronecker products as described in [29] and the boundary value problem solved using a Newton trust region-based algorithm [30]. We also embedded the boundary value problem into a secant continuation formulation to facilitate path following of the solutions in parameter space [31].

The scheme was implemented in Matlab (v. 2015a) with typical discretizations in *x* of *N*_{x}=41 mesh points on *x*∈[0,*π*] and *N*_{y}=21 Chebyshev collocation points in *y*.

### (a) Periodic solutions

As mentioned earlier, it is frequently the case that for plate problems focus is concentrated on periodic behaviour, particularly in the macroscopic example of engineering structures. In such instances, the wavelength of the buckle pattern is of the same order of magnitude as the overall dimensions of the plates so only a relatively small number of waves are expected to span the whole domain. The aspect ratio is also important; a plate that is very much thinner in one dimension than the other might conceivably behave somewhat like a strut [10], but wider plates may well exhibit very different behaviour. We take some tentative numerical steps towards investigating this variation in §5c.

The earlier Galerkin analysis suggested that solutions which are periodic in both *x* and *y* can bifurcate from the trivial flat state and they do so in pairs. However, we now appreciate a limitation of the simple-minded Galerkin approximation used above, because it is evident that, since the foundation is asymmetric with respect to deflection, the pair of modes that evolve are not mirror images of each other. In order to model this asymmetric behaviour, more than a single mode would be required to break the implicit symmetry [11]. Nevertheless, the single-mode approach gives us a firm connection between the one- and two-dimensional studies. A typical example of the periodic post-buckling regime is depicted in figure 4. Here the two branches emerging from the flat state are labelled *A* and *B* and shown in different colours. In *P*–*Δ*_{x} space, both branches overlap exactly and it is therefore more informative to present the plot with one extra dimension (here, deflection at the origin) to prise the two branches apart and show that they really are distinct. The buckled shape corresponds to *m*=3, *n*=1, as defined in (3.2). With *P*_{y}=2, the prediction from the one-dimensional approach in §4 is that at bifurcation *P*_{x}≈17.03 (see (4.12)) while the numerical solution of the full equation reveals that the bifurcation occurs at around *P*_{x}=22.7. Each point on the curve is a solution of the static equilibrium equations, although the energy *V* does vary considerably. The value of the energy can be loosely thought of as an informal relative measure of how distorted the deflected surface is, starting at zero at the bifurcation point (Bif) and dropping off as the striped solution is approached (St) (although as the striped solution is obviously not flat, the value of the energy is not zero).

Some solutions at various values of *P*_{x} are shown in figures 5 and 6. In figure 5*a*, we see the buckle pattern, which emerges from the flat state *w*=0. Of course, this solution relates to just one side of the bifurcation point as a second pattern, which is not quite the mirror image of the first, evolves on the other branch. Figure 5*b*,*c* depicts solutions further into the post-buckling regime, and it is much easier to see that these solutions are asymmetric about the flat state owing to the bias in the foundation introduced by the quadratic term.

Solutions on each branch continue to maintain the same structure as *P*_{x} is varied except that the amplitudes of both the peaks and troughs increase steadily until the branches near the striped solution St which is shown in figure 6. This last state corresponds to the situation when the two branches eventually coalesce, thus forming a closed bifurcation curve (isola) on which the solution appears to be trapped.

#### (i) Striped solution

It is evident that a special example of a periodic solution is that discussed in §4 where the deflection is purely a function of *y*. The solution has been shown to be stable to a two-dimensional perturbation only if the axial load *P*_{x} lies within a certain range. Our computations above indicate that the solution represents a transition from solutions of one parity to another. The numerical results also confirm that *Δ*_{x} is zero both at the bifurcation from the trivial state *and* at the striped state (figure 4*a*). A very important distinction between these two limiting configurations is that the potential energy is quite different (figure 4*b*). The outcome of this is that the two states are not interchangeable and a straightforward transition from one to the other is impossible.

For the same value of *P*_{y} (=2), two further families of solutions can be found that both bifurcate from the fundamental state but that do so at different values of *P*_{x} (figure 7). The emerging solutions are distinct, having a different number of waves in the *y*-direction but, as these solutions evolve, they converge to the same striped solution presented earlier. Though kinematically similar, it should be remembered that *P*_{x} is different in the two instances although the final striped solution has potential energy *V* ≈14.649 (equation (2.8)). The potential energy of the striped solution is lower than other equilibrium states in the vicinity, which satisfies the axiom for static stability [1] and numerically agrees with our analysis of dynamical stability in §4 that such solutions should be statically stable.

Two further bifurcation diagrams are presented for a different value of transverse load: *P*_{y}=7 (figure 8). It can be seen once more that, although the solutions emergent from the trivial state have separate forms, they again coalesce onto the same striped solution. In this case, *V* ≈73.821.

#### (ii) Transverse stability of the striped solutions

It is evident that the same striped solution exists at the limit point of several solutions bifurcating from the flat state (e.g. figures 6*e* and 7*c*, where three distinct two-dimensional bifurcating patterns share the same striped pattern, or figure 8*c*, for two others for a different value of *P*_{y}). The forms of these sets of one-dimensional solutions are dictated solely by the value of the load parameter *P*_{y} only because the end-shortening in the *x*-direction, *Δ*_{x}, is zero and so there is no contribution to the potential energy from *P*_{x}. On the other hand, the form of the solutions emerging from the trivial state is dependent on both *P*_{x} and *P*_{y}. So, the question that inevitably arises is: how can the same striped solution be linked to multiple two-dimensionally periodic forms? The answer is guided both by the approximate analysis carried out in §4 and by the numerical simulations.

Theoretical analysis reveals that the striped solution remains dynamically stable as long as
*P*_{x} is larger than this limit then *η*_{1} (4.12) is real and the solution will grow exponentially in time. The striped solution under consideration closely resembles mode 1 in figure 2, and it is determined that the critical value of *P*_{x}≈19.7. This compares reasonably well with the numerically computed limit values of *N*_{1} and *N*_{2} for the full equations which lead to the prediction *P*_{x}≈22.7.

By numerically finding the striped solution for various values of *P*_{x}, bifurcations can be detected by numerically estimating when an eigenvalue of the approximate solution becomes zero. The structure of the eigensolution closely resembles the deflection pattern in the vicinity of critical buckling and a bifurcation is detected at around *P*_{x}=10.17. The unstable eigensolution at this load is plotted in figure 9*a*, and it can be seen that, although the actual deflected shape is one dimension, the eigensolution reveals the direction in which any disturbance would grow. Thus, we conclude that, in the vicinity of this load, the solution would naturally evolve towards a pattern indicated by that unstable solution. Further transverse instabilities can be detected at various values of *P*_{x} (while holding *P*_{y} constant) and some are shown in figure 9*b*,*c*. These unstable eigensolutions resemble the buckling patterns in figure 7. The numerically determined values of *P*_{x} correlate very well with the values at the limit points of the bifurcation curves. Thus, we conclude that the same striped pattern can indeed be linked to various two-dimensional patterns in a systematic way.

We have found that for each given *P*_{y}, there is an associated striped solution. However, we can subsequently induce buckle patterns of different types by then loading the plate in the orthogonal direction, *P*_{x}. Theoretically, this could be used as a way of producing the desired deflection pattern in a repeatable manner. A consequence of the phenomenon of multiple patterns being able to develop from a particular one-dimensional striped mode could be used to manufacture a layer of material with a desired pattern in a two-step process. Firstly, the specimen could be loaded in one in-plane direction to give a striped buckle shape and then judicious loading in the orthogonal direction would give rise to the required two-dimensional pattern. In particular, carefully varying the secondary load would, in principle, be a possible method of manufacturing materials with various tuneable patterns and hence physical properties [24].

### (b) Localized behaviour

As we indicated in the previous section, periodicity is the predominant form of behaviour to be expected in most plated structures, because the finite domain precludes genuinely localized solutions which are only evident over an infinite domain. However, if the domain is significantly larger than the natural wavelength of the buckle pattern then a (periodic) solution can arise which comes to resemble its truly localized counterpart in the sense that, over the whole extent of the plate, only a small fraction undergoes notable deflection. Localization of patterns can be favourable under certain conditions, particularly where a softening nonlinearity makes such a pattern energetically more favourable to jump to than its periodic counterpart [11] even though both are kinematically admissible. Even in the case that such localized solutions are unstable, they can form critical (mountain pass) modes [32] to other forms. Thus, a localized solution may develop in a structure of finite extent if at one end the solution and all its derivatives approach zero. Such is the case for the buckling of a thin cylindrical shell [33]. In other applications, such as contemporary material science, periodic buckling has been the primary focus [34]. However, if the extent of the domain is much larger than the wavelength of buckle patterns, then isolated (localized) patterns can coexist with periodic solutions.

To demonstrate the relationship between periodic and localized forms, figure 10 shows the bifurcation structure of two plates of non-dimensional lengths 10 and 20. These plots can be interpreted as localized patterns with symmetry about the position *x*=0. As can be seen, a pair of solutions bifurcate from the trivial (zero-energy) state; the other is a subsidiary solution forming a multi-humped pattern [13]. We infer that as length increases the limiting behaviour will be localized states that bifurcate from the unbuckled state and coexist with their periodic counterparts just as in the one-dimensional strut model [11]. For a finite-length plate, the way in which a localized pattern may be realized is by an initial periodic pattern which is then disturbed, perhaps by a secondary bifurcation, into a localized state. Such multi-staged buckling phenomena are common in plates [22,23].

### (c) Aspect ratio

The problem as formulated in this study consists of an investigation into the theoretical buckling aspects of a structure that extends into two dimensions. The equivalent one-dimensional problem (the strut model) has been extensively studied but the plate on a foundation has attracted far less interest. There are many potential avenues to explore for this model and we have hinted at some possible ideas here. In many ways, it is perhaps the role played by the aspect ratio that is most fundamental, so we introduce the aspect ratio, *l*_{y}≡*γ*/*L*, as a continuous parameter. In figure 11, we illustrate the evolution of a basic periodic solution as *l*_{y} is varied from its original value *l*_{0}=2/*π*. As the ratio is increased so the pattern spreads to fill the extra width and yet the buckled cells approximately maintain their own aspect ratio. As the width gets larger, there are a few transitions when the number of cells increases by one each time; this is clearly a form of mode jumping [20]. Further evidence that lends credence to this conclusion is provided in the plot of potential energy versus *l*_{y} as the curve exhibits undulations at the point where the number of buckle cells changes.

## 6. Conclusion and further study

The current work has examined the fundamental behaviour of a structural system in order to investigate the effect of adding an extra spatial dimension to a well-studied one-dimensional problem. As expected, the introduction of a finite width opens up a new richness of possibilities of solutions available after buckling ensues. The loading regime is seen to be enhanced with the ability to apply different forces in different directions and indeed to allow the loading to vary spatially, although that has been left for future work.

We have seen that some conventional techniques, such as the Galerkin method, are able approximately to capture buckling behaviour emergent from the flat state but fully numerical methods are needed to follow what happens to these solutions further into post-buckling. We found that there are special one-dimensional (striped) patterns which eventually connect both the branches which emerge from the critical state and furthermore that the same striped pattern (governed by the transverse load, *P*_{y}) can lie at the end of multiple branches. The unstable eigensolutions of the striped solution reveal which two-dimensional (spot) solution it links to. An interesting possibility arises from our findings in that we have outlined a plausible procedure by which material could be manufactured with pre-determined characteristics by careful application of load in one in-plane direction followed by loading at right angles to it. This technique could be a simpler alternative to current methods used to manufacture highly tuned nanomaterials.

A straightforward analysis of the striped solution using an adaptation of the strut model [10] reveals information about the stability of solutions and this, coupled to numerics, has enabled us to locate the boundary between stable and unstable behaviour. The boundary gives a reasonable approximation for the longitudinal load at which the striped solution gives way to spots or cellular buckle patterns.

Although we have considered the loading as two independent parameters in this study, it may be instructive in future to allow the load to vary along either or both edges in order to simulate various conditions that may arise in reality. An example of such varying stress may arise through thermal expansion or contraction of a plate embedded in an elastic medium. Some consideration has been given to localized behaviour and the effects of aspect ratio. These features would increase in importance were we to look at the case where the buckling cells are small compared with the length and width of the plate.

The plate–foundation model has been used to mimic many physical situations ranging from material wrinkling in industrial manufacturing [35], theoretical neo-Hookean materials [36] and fold hierarchies in thin films [37] through to the morphological development of tissue and organs such as the brain [38]. However, most modelling usually starts with some simplifying assumptions including the reduction of the problem into one-dimensional phenomena. We offer here a way to begin systematically to study two-dimensional behaviour, which holds the possibility of uncovering a much richer set of competing and interacting solutions. Although resolution is presently far off, progress in this direction may conceivably go part way to explaining some of the complex behaviour observed in practice.

## Data accessibility

Calculations were performed using equations which are all published in this paper and numerical computations were carried out using published algorithms with Matlab.

## Authors' contributions

M.K.W. conceived of the problem; M.K.W. and D.J.B.L. designed the programme of study and M.K.W. undertook the Galerkin calculations and computations based on numerical algorithms developed by D.J.B.L. M.K.W. and A.P.B. worked on the reduced one-dimensional stability problem and D.J.B.L. devised a numerical verification method for these results. All authors contributed to the writing of the paper and approved its final submission.

## Competing interests

The authors had no competing interests in performing this work.

## Funding

The research was funded as part of the authors’ standard employment contracts.

## Acknowledgements

We express our thanks to Prof. Geoffrey Nash at the University of Exeter for helpful discussions regarding the manufacture and properties of modern two-dimensional materials and undergraduate David Frampton, also at the University of Exeter, for helping to develop the Matlab code to find the Galerkin approximations in §3. Finally, we are grateful to the referees for their helpful comments, which considerably enhanced the presentation of this work.

- Received December 21, 2015.
- Accepted March 11, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.