In the deformation of layered materials such as geological strata, or stacks of paper, mechanical properties compete with the geometry of layering. Smooth, rounded corners lead to voids between the layers, while close packing of the layers results in geometrically induced curvature singularities. When voids are penalized by external pressure, the system is forced to trade off these competing effects, leading to sometimes striking periodic patterns. In this paper, we construct a simple model of geometrically nonlinear multi-layered structures under axial loading and pressure confinement, with non-interpenetration conditions separating the layers. Energy minimizers are characterized as solutions of a set of fourth-order nonlinear differential equations with contact-force Lagrange multipliers, or equivalently of a fourth-order free-boundary problem. We numerically investigate the solutions of this free-boundary problem and compare them with the periodic solutions observed experimentally.
(a) Folds in multi-layer structures: rocks and paper
Geological folds arise from the in-plane compression of layers of rocks from tectonic plate movement, and provide one of the overarching themes of this Theme Issue (see also [1,2]). Geometric constraints, following the need for layers to fit together, generate a number of distinctive forms of folding, not seen in single-layered systems (see Price & Cosgrove ). A couple of illustrative examples are shown: figure 1a is an example of parallel folding where, as the eye moves from bottom left to top right, the curvature of the layers increases until eventually it becomes too tight and forms a distinctive swallowtail singularity , after which the layers are seen to crumble. In figure 1b, the geometric constraints have been overcome by the rocks adopting straight limbs and sharp corners, in what is known as chevron folding . This phenomenon is strongly reminiscent of kink banding, seen also in other contributions in this issue [6,7].
The sequence of events leading to such structures is of course lost in geological time, but small-scale experiments on layered materials, such as papers, do provide some suggestions. The outcome of one such experiment  is shown in figure 2. About 800 half sheets of A4 size cut longitudinally were squashed together between two steel plates by a transverse load; subsequently, the axial displacement was increased by a corresponding axial load, whilst the transverse displacement was kept constant. Load cells were used to record both loads at 1 s intervals, and an in-plane transducer registered the axial end displacement. To aid visualization, a black-edged sheet was introduced approximately every 25 layers.
Output from a typical experiment is shown in figure 2b. Kink bands formed in sequence from the loaded end, starting at the point B, followed in turn by D, F and G, as seen on the output diagram. Each band forms with a clear-cut instability, observed by a corresponding instantaneous fall in applied axial load and followed by subsequent re-stabilization or lock-up. By comparison, figure 3 shows various stages in the loading history of a similar experiment, but with soft foam layers replacing the upper and lower (stiff) supporting foundation, such that the overburden pressure and accompanying stiffness are much reduced. The folds still develop and lock-up sequentially, but the sharp corners and straight limbs of the kink bands are replaced by curved layers, as in the parallel folding example of figure 1; lock-up is here associated with the tightening of the curvature on the inside of a fold .
The test shown sequentially in figure 4 starts under conditions similar to parallel folding, but evolves to a situation more like kink banding. Foam is again used as a supporting medium, but at the onset of loading, only a small overburden pressure is applied, so effectively an Euler buckle initially appears over one half wave; this is just apparent in the left-hand image of figure 4. The long wavelength means that deflections into the foundation grow relatively rapidly, and quickly bring the outer steel plates into play, stiffening the foundation. This strongly nonlinear effect both flattens the crest of the wave and instigates two outlying kink bands orientated across the sample at about 30° to the transverse direction, seen most clearly in the final picture. Although these appear similar to those of figure 2, they are formed under quite different circumstances, as follows.
When a kink band is orientated obliquely across the sample, as it displaces from its undeflected (flat) state, the layers initially dilate within the band . Those seen in figure 2 have used this property to select their angles of orientation across the sample, such that maximum available dilation exactly cancels the pre-compression induced by the overburden pressure, and friction is eliminated. They propagate across the sample in this optimum configuration, in which work done against friction is minimized; this orientation is then maintained as the band develops and eventually locks up.
The same is not true for the system of figure 4. The development of the band proceeds slowly from the onset of loading; its orientation is not freely chosen but imposed from the start, and dictated by factors such as sample length and overall thickness. If the angle of orientation of the band is greater than its optimum value, voids will result, as seen here. Note that, although the introduction of a black-edged sheet clearly helps to induce voids, they appear as a periodic sequence every five black sheets, or about once every 125 layers. Very similar effects are also found in paperboard creasing and folding, as seen in fig. 8 of Beex & Peerlings  (this issue).
(b) A model of multi-layer folding with voids
In this paper, we focus on how bending stiffness, overburden pressure and geometric constraints combine to cause or restrict the appearance of voids. The intriguing periodicity of the voids in figure 4 serves as an example: is it possible to explain, or even predict, the way voids concentrate at certain layers?
We do so by formulating a simple model of a multi-layer material in which we eliminate all but the essential aspects. The layers are described by (linear elastic) Euler–Bernouilli beams, which are inextensible and laterally incompressible; the contact condition is ‘hard’, and the overburden pressure is modelled by a simple global energetic expression (see §3). Geometric nonlinearity is preserved, and the system will undergo large rigid-body rotations.
Such a model is highly simplified in comparison to other studies in the literature (see Hudleston & Treagus  for a recent review), and differs from a large part of the literature by assuming elasticity rather than viscosity. By studying such a simplified model, however, the appearance and size of voids are precisely defined, and we have detailed access to their properties. As for the elasticity, we believe that for the type of question that we pose here, elastic and viscous materials behave sufficiently similarly on the relevant time scales. In other work, we have seen how the assumption of elasticity allows for precise analysis [4,8,11,12].
Our previous study of the behaviour of voids concentrated on a single infinite elastic inextensible beam forced into a V-shaped corner under uniform pressure q . We generalize this model to multi-layer folding in two steps. First, we recast the model of voiding of Dodwell et al.  to embrace finite-length beams and axial loads. The introduction of axial load negates the favourable properties found for minimizers of the purely pressured case: minimizers are no longer necessarily convex, symmetric or have a single void space. It is shown that, under various loading conditions, minimum energy can lie with buckling on a local length scale [13,14], demonstrating higher order accommodation structures.
The main section of this paper then combines the multi-layered setup presented in Boon et al.  with the elasticity properties of §2 to capture the periodic behaviour of figure 4. A first experiment in this direction is to stack identical copies of elastic layers; an elementary geometric argument suggests that reduction in the total void space could force such a stack into a straight-limbed, sharp-cornered configuration, as seen in figure 5.
In §3, we construct this experiment mathematically by deriving a potential energy function for a system of N identical layers that are forced into a V-shaped corner by a uniform axial load P and overburden pressure q, to produce the so-called 1-period solution. The energy for this specific case may be simply expressed in terms of a single layer, which leads to a similar analysis to that of Dodwell et al. , as presented in §2. The results match our initial elementary argument and show that the total void space reduces as the overburden pressure increases, giving deformation with straight limbs and sharp corners. If we then consider m-period solutions, where one void forms every m layers and then repeats, we are led to an interesting question. Which m-periodic solution has minimal energy? Constructing an energy for m-periodic solutions and the subsequent variational analysis is non-trivial, and we therefore leave this to a future publication. We do, however, offer some preliminary energy calculations, which suggest that the energy typically is minimized at an m that is larger than 1 but still finite. When we choose reasonable values for the parameters in the experiments described in figure 4, the predicted value of m is of the same order of magnitude as the observed value, giving us hope that there is more interesting research to follow.
2. Single-layer model for void formation
The potential energy functional, V (w), comprises strain energy stored in bending, UB, and work done against overburden pressure in voiding, UV, so that 2.1where B is the elastic bending stiffness of the layer and q the overburden pressure per unit length. Note that UB and UV are measured from different zero-energy configurations: the former is zero in the flat state and the latter zero in the V-shape. For finite-length beams, the two can be reconciled by noting that the work done by q in moving from the flat beam configuration to the displaced state is a constant (the area of the triangle between the V and the flat beam) minus the work done against q in voiding: constants of energy are known not to affect equilibrium equations. The same idea extends in the limit to beams of infinite length.
By considering stationary solutions of (2.1) with respect to vertical displacement, equilibrium for the system is described by the following Euler–Lagrange equation: 2.2where μ is a Lagrange multiplier, arising from the pointwise obstacle constraint w≥f. Physically, this quantity can be interpreted as the vertical component of the reactive contact load with the boundary f; therefore, μ is zero wherever there is no contact, matches the overburden pressure q along stretches of contact, and equals the discontinuous shear force at points of delamination. Dodwell et al.  prove that there exists a unique minimizing solution to the Euler–Lagrange equation (2.2), and that such solutions are convex and symmetric, and have a single finite interval in which they are not in contact with the boundary.
Equation (2.2) is written in Cartesian coordinates (x,w), but may be rewritten in intrinsic coordinates parameterized by (s,ψ), where s is the arc length and ψ=dw/ds. As an aside, both representations have their uses: intrinsic coordinates give simpler equations, but for multi-layered problems, they carry the added complication of a different coordinate frame for each layer. Although not shown here, it is often useful to switch freely between coordinate systems, using the benefits of either when it suits (Dodwell et al. , §3 provides some details on how this can be done).
The formulation of Dodwell et al.  is readily extended to embrace the different configurations of a long beam of finite length L, with dead load P acting down the slope (figure 7). In this work, we are primarily interested in the bottom picture, where the combination of P and q has pushed the beam into a similar configuration as figure 6, such that some portion of the length is in contact with the supporting boundary. The shortening Δ, compared with the flat state wx=0, and measured along the limb, is given by 2.3where ℓ is the horizontal non-contact length (figure 6). The work done by the compressive load P is then PΔ and after ignoring a constant, the total potential energy over the non-contact region can be written as 2.4This has the Euler–Lagrange equation 2.5We note that the effect of the dead load P can be replaced by rigid loading by introducing the extra constraint Δ=constant, where the load P becomes a free parameter. It can be argued that rigid loading more closely resembles the experimental set-up; however, both forms of loading share the same equilibrium solutions.
(a) Numerical solutions
We now consider solutions that minimize V subject to the constraint Δ=constant. With the additional loading constraint, the favourable properties of minimizers for the purely pressured case (unconstrained minimization) may cease to hold: constrained minimizers of V are not necessarily unique, symmetric or convex. Indeed, although there remains the possibility of higher order asymmetric solutions [15,16], we shall here restrict ourselves to symmetric solutions that have a single finite interval for which they are not in contact with the obstacle f.
The Lagrange multiplier μ can be eliminated by writing equation (2.5) as a free-boundary problem . Numerically, this problem can then be treated as an initial boundary problem with at most two unknown shooting parameters, ℓ and P, depending on the form of loading. Solutions of (2.5) are found by shooting from the free boundary at x=−ℓ with the initial conditions 2.6towards the target symmetric section at x=0 given by 2.7These initial conditions (2.6) can be found by integrating (2.5) once, and imposing continuity of w, wx and wxx at points of delamination x=±ℓ.
By dividing (2.5) by the elastic stiffness B, it becomes clear that solutions depend on the two variables q/B and P/B. Figures 8 and 9 show a number of numerically obtained solutions found by varying these two parameters. Figure 8 gives the equilibrium solutions for P=1 and increasing values of q/B, whereas figure 9 displays plots of one of the equilibrium solutions for q=1 and increasing P/B. The important feature to note is that as either q/B or P/B increases, then the size of the delamination decreases, giving an equilibrium profile with straight limbs and a sharper corner.
Figure 10 shows numerical solutions of (2.5), on a plot of load P against its corresponding deflection, as measured by Δ, at constant values of q, B, L and k. This plot shows two of the many possible types of solution. From P=0, we first follow the solid line corresponding to the simple convex solution of figure 9. We note that as the layer is forced into the corner singularity, , then . Figure 10 marks the Maxwell displacement, ΔM, where energy levels in the two equilibrium states are the same . If Δ>ΔM, minimum energy rests with the up-buckled shape. In particular, Δ=L/2 gives an upper bound for the existence of simple convex solutions beyond which only higher order solutions, such as that indicated by the dashed-line curve, exist. It is worth noting that such higher modes apparently do exist in geological structures, known as accommodation structures , p. 323. Figure 11 shows one example from a rock outcrop in North Wales. The possibility of energy minimizers involving buckling on a smaller scale is particularly interesting, but we leave the investigation of such convoluted solutions to a future contribution.
3. Multi-layered model for void formation
(a) Multi-layered formulation
We now extend the single-layered model to a multi-layered formulation. Consider a stack of N identical layers of uniform thickness h and length L forced into a V-shaped singularity by an overburden pressure defined by q on the top layer, and axial load P per layer, as seen in figure 12.
The centreline of the ith layer from the bottom is labelled as wi(x), and the top and bottom boundary of that layer as and , respectively. A constant axial load P acts on the stacks, which shortens the ith layer by a distance parameterized by Δi. As before, Δi is measured from the centre of the fold, and here additionally ξi gives the horizontal length of the layer. We also introduce the obstacle constraint for i=1 to N−1 and w1(x)≥f, so that adjacent layers cannot interpenetrate or pass through the boundary f=k|x|.
We now construct a potential energy function for the system, by adding together the contributing parts.
(b) Bending energy
The total bending energy of the layers is the sum of the bending energy of each layer, given by The quadratic dependence on implies that sharp corners have infinite bending energy and, therefore, for finite overburden pressure, layers will void rather than producing singular corners.
(c) Work done against overburden pressure
In experiments with paper, where voids occur, there is clearly a subtle interplay between release in overburden pressure in these areas, and some compensatory increase in internal strain energy elsewhere. Rather than attempting to model this complex situation directly, we simply assume here that the creation of all voids does work against a constant pressure q, regardless of position in the sample. Further development is clearly necessary for application to experiments like that of figure 4. However, from the experiment presented in figure 2, it appears that the transverse load in the post-buckling range may be a good indicator of the work done in forming voids. This provides motivation to carry out a series of more accurate experiments of type presented in figure 4 where transverse load is also measured, but we leave this for future work.
The overburden pressure acting on the layers is q per unit length, and therefore the total work done against overburden pressure is q multiplied by the sum of all voids between the layers, given by The first term gives the energy associated with the void formed between the bottom layer w1 and the boundary f, and the second calculates the size of the voids between each of the remaining layers. If q is large, UV implies a severe energy penalty.
The fundamental assumption of engineering bending theory that plane sections remain plane and normal to the centreline means the top and bottom boundaries of a layer can be calculated parametrically by propagating the centerline forwards and backwards a distance h/2. In this case, it is necessary for the bottom of one layer, , and the top of the next, , to be defined according to the common fixed coordinate system (x,w). At a given x, we wish to calculate the void between the ith and (i+1)th layer, defined as . Figure 13a shows that can be determined by propagating wi+1 backwards from the point (i), and by propagating wi forwards from (ii). These centreline positions are not known a priori, but can be reverse-engineered from the common position x by the first-order approximation 3.1Figure 13b demonstrates this approximation, where (iii) and (iv) are found by constructing right-angled triangles. This is a thin-layer approximation, in which we assume that the layer is much thinner than the radius of curvature of the layer. Whether this is justified or not depends on the situation. We also note that the approximation is exact for straight segments, whereas for the convex segment shown in figure 13b, (iii) is an underestimate and (iv) an overestimate; and vice versa for a concave segment. With this approximation, the total energy in voiding becomes
(d) Work done by load
If we consider a constant axial load P acting on each layer of the multi-layered stack, the total work done by the load is the sum of the work done by the load in the shortening of each layer. Therefore, up to a constant, the total work done by the load is
(e) Total potential energy function
The total potential energy for the complete system can be written as the sum of work done in bending by each layer, total work done against overburden pressure in voiding and work done on the system by the axial load P in shortening, as follows: 3.2Appropriate differentiation leads to a system of N nonlinear fourth-order equations, coupled by contact conditions.
4. Periodic solutions
The formulation above gives a general energy description for N identical layers forced into a V-shaped singularity. Now, motivated by figure 4, we restrict ourselves to minimizing V over the class of symmetric periodic functions with one void per period. We also make the assumption that layers are symmetric, with straight limbs of fixed angle k connected by a convex curved section. From §2a, we recognize that higher order solutions and asymmetric solutions might exist, but the experiments shown in figure 4 suggest that the assumption of this simplified geometry is justified.
(a) 1-periodic solutions
We start with the simplest periodic solution where the structure repeats every layer, hereafter named a 1-periodic solution. Since each layer deforms identically, the total potential energy may be described in terms of a single layer. Each layer does the same work in bending and shortening, and the size of the void can be calculated by shifting the bottom boundary of a layer to fit above the top of the same layer. Therefore, The total potential energy (3.2) simplifies to 4.1where the shortening Δ is again given by (2.3). Minimizers of (4.1) solve the Euler–Lagrange equation over the non-contact region, given by 4.2As in §2, boundary conditions can be found by integrating (4.2) with respect to the spatial coordinate x and imposing continuity of wx and wxx at points of delamination x=±ℓ, so that 4.3coupled with the symmetric boundary condition at x=0. Figure 14 shows an example of solution profiles obtained by this numerical shooting method, as well as the hinge zone of a chevron fold in a recent paper experiment that demonstrates a 1-period solution. We note that the results match physical intuition in that as overburden pressure q or load P increases, then the void size decreases, leaving in the limit deformations with straight limbs and sharp corners. This is demonstrated in figure 15, which shows how the size of the void varies with angle of the limbs k and increasing values of q/B.
(b) m-periodic solutions
Section 4a naturally extends to m-periodic solutions, where a void forms after m layers and then the structure repeats. This leads to various interesting questions. Is there a value of m, not always 1 or , which has minimum energy among the class of such periodic solutions? Can we predict the value of m=125 seen in the periodic solutions presented in figure 4? As we have seen, 1-periodic solutions can be solved using previously developed analysis . This is not the case for more general m-periodic solutions, where deriving the energy and the subsequent variational analysis is a non-trivial exercise. We therefore leave the details of this work to a future publication. To test the basic ideas, however, we perform a very rough calculation for the experiment presented in figure 4, to see if a simple energy analysis, with reasonable values of the physical constants, gives a value of m of the correct order of magnitude.
The central premise is that in equilibrium, the bending energy in a single m-periodic packet is approximately equal to the energy required to form one of the voids. For this analysis, we assume that each layer has straight limbs connected by a circular arc (figure 16), where the circular arc of the outside layer is of radius R swept through an angle 2θ. The angle θ is fixed for a given gradient of the limbs k. The bending energy stored in an m-periodic packet of layers is therefore where κi is the curvature and si is the arc length of the ith layer. Here, we number the layers as 1 to m, from the outside inwards. Following a result presented in Boon et al. , κi=κ1/(1−(i−1)hκ1), and it therefore follows that For large m, we consider f(κ1) to be a Riemann sum, and approximate it by the corresponding integral, i.e. setting xi=h(i−1), by which For a circular arc, κ1=1/R is constant over θ and therefore 4.4Estimating the work done to produce a void as simply the overburden pressure q times the size of the void, we can equate the forms of energy, so that 4.5Table 1 displays estimates of the physical constants for the experiment and gives a brief description of their source.
By substituting the values of the parameters from the table into (4.5), we obtain giving a rough estimate of m≈17.
The fact that an intermediate value of m is chosen suggests that the basic premise of the calculation (balance of bending and voiding energy) might be reasonable. The fact that the predicted number m≈17 is in the same ballpark as the observed value of 125 in figure 4 further reinforces this suggestion.
5. Concluding remarks
This paper extends a single-layered model for folding and void formation  to the simplest periodic multi-layered structure, including the vital effects of axial loading. The model extends our insights into the interplay between material and loading properties on the one hand, and the geometric constraints of layers fitting together on the other. The results match physical intuition in that as overburden pressure increases, void size decreases, in the limit leading to deformations with straight limbs and sharp corners. Although an oversimplification in terms of the paper experiments, and even more so of their true geological counterparts, it does capture many of their important hallmarks, as well as demonstrating the possibility of periodic solutions arising naturally in an elastic context.
Studying these 1-periodic solutions has led to interesting questions about the possibility of m-periodic solutions. In particular, can we find which m-periodic solution has minimum energy? The simple energy argument presented in the previous section gives a rough estimate for m for the experiment presented in figure 4. This gives an interesting path of research to follow, and one we wish to consider in a future contribution.
One contribution of 15 to a Theme Issue ‘Geometry and mechanics of layered structures and materials’.
- This journal is © 2012 The Royal Society