Convoluted accommodation structures in folded rocks

February 26, 2012 19:10 Philosophical Magazine PhilMagFinal Philosophical Magazine Vol. 00, No. 00, 00 Month 2011, 1–18 RESEARCH ARTICLE Convoluted accommodation structures in folded rocks T. J. Dodwell and G. W. Hunt∗ Department of Mechanical Engineering, University of Bath, Bath, BA2 7AY, UK. (Received 00 Month 200x; final version received 00 Month 200x) A simplified variational model for the formation of convoluted accommodation structures, as seen in the hinge zones of larger-scale geological folds, is presented. The model encapsulates some important and intriguing nonlinear features, notably: infinite critical loads, formation of plastic hinges, and buckling on different length scales. An inextensible elastic beam is forced by uniform overburden pressure and axial load into a V-shaped geometry dictated by formation of a plastic hinge. Using variational methods developed in [1], energy minimization leads to representation as a fourth-order nonlinear differential equation with free boundary conditions. Equilibrium solutions are found using numerical shooting techniques. Under the Maxwell stability criterion, it is recognized that global energy minimizers can exist with convoluted physical shapes. For such solutions, parallels can be drawn with some of the accommodation structures seen in exposed escarpments of real geological folds. Keywords: accommodation structures, localization, obstacle problem, multilayered geometry, geological folding, Maxwell stability criterion. 1. 1.1. Introduction Accommodation structures in multilayered folds Multilayered structures found in the Earth’s crust may arise from a sequence of sedimentary deposits, subsequently held together by additional overlying layers. Tectonic plate movement causes in-plane compression of these layers of rocks that can fold them into convoluted patterns. The resulting deformation is not only dependent on the complex mechanical and material properties of each layer, but on the highly nonlinear geometric constraints of them fitting together. It is widely accepted that a layered system buckles to a wavelength determined primarily by the bulk material properties of the stack [2, 3], and as a result less competent layers can be forced to accommodate the overall geometry imposed. In some cases layers can be obliged to fit into very tight or even singular geometries [1, 3–5], and as a result various patterns, named accommodation structures [6], may form on a localized length scale. These structures include accommodation thrusts, keel-like hinges [3, Ch. 12], void or saddle reef formation [1, 5] and in particular higher-order localised buckles. Two examples of such localised buckles are seen in Fig. 1. Developing a simplified structural model to understand this particular type of accommodation structure is the primary purpose of this paper. Dodwell et al. [1] model the behaviour of an elastic layer forced by a constant overburden pressure into a simplified V-shaped geometry, where the layer is freely allowed to slip on its limbs. Constrained energy minimization shows that any such solution is convex, symmetric and unique [1]. Since the model gives only convex ∗ Corresponding author. Email: [email protected] ISSN: print/ISSN online c 2011 Taylor & Francis DOI: http://www.informaworld.com February 26, 2012 19:10 Philosophical Magazine 2 PhilMagFinal T. J. Dodwell and G. W. Hunt Figure 1. Convoluted accommodation structures in the hinge region of larger multilayered folds; (left) outcrop from N.Wales (photo courtesy of J. Cosgrove); (right) Porthleven, South Cornwall. solutions we see that pressure, by itself, cannot initiate the formation of higherorder folds like those in Fig. 1. Any valid model must therefore embrace not only the complex topological effects of layers fitting together, but also those of slip relative to each other [7, pp. 167-168]. In the inevitable presence of friction, axial loads and/or corresponding constraints hence become part of the story, in the formation of such convoluted accommodation structures. 1.2. The model To encapsulate the formulation of these structures at singular geometries, we construct the simple mechanical model of Fig. 2. A knee joint is made up of two rigid bars of length L/2, hinged by a plastic rotational spring with the static response characteristic shown at the bottom left; no rotation θ can take place at the hinge until the moment in the spring reaches the plastic moment Mp , whereupon it rotates against this constant resistive moment. This knee joint represents the mechanics and geometry of a larger fold or competent layer. A weaker accommodating beam of bending stiffness EI is pinned to it at each end, where in-line springs of stiffness k2 , which are unstressed in the flat state, undergo deflections Q representing relative slip between the components. The accommodating structure is taken as inextensional but is allowed to buckle, as seen at the bottom right, with displacements described by the angle ψ(s) where s represents measures along the buckled length as shown. The system is loaded horizontally by a load P , displacing by an amount ∆, and there is a further spring k1 at the point of loading to give an initial, finite (pre-buckled) stiffness. The entire system is placed in a bath of pressure q, against which work needs to be done to create voids. Dodwell et al [1] simultaneously define layer geometry in cartesian coordinates, characterised by the vertical displacement u(x) with respect to fixed horizontal coordinate x, both being measured from the center of the hinge. We note that the ability to switch freely between coordinate frames is an important part of these analytical developments, each description having its benefits. Intrinsic coordinates generally give concise governing equations, but a common fixed coordinate system using cartesian coordinates proves useful in describing voiding between layers [5]. The linear elastic springs of stiffness k2 (see Fig. 2), control the amount of slip between the layers, and are taken to mimic axial compression in extended lengths of the accommodating layer which are able to unload into the buckle. The modelling February 26, 2012 19:10 Philosophical Magazine PhilMagFinal 3 Convoluted accommodation structures in folded rocks L cos θ ℓ k2 k2 pressure q Q k1 ∆ P Q B θ Mp M Mp θ Figure 2. This figure shows the setup of the model and parameters discussed in this paper. Solution profile demonstrates a typical up buckled profile. of such lengths to put a value to k2 would depend heavily on frictional properties at the interface and would certainly be non-trivial in practice, mirroring that of the selection of so-called active lengths in pipeline buckling problems [8]. We define u(x) over the set X = [−L/2 cos θ, L/2 cos θ], as well as the contact set Γ = {x ∈ X : u(x) = f (x)} , together with its complement the non-contact set Γc . Throughout this work we restrict ourselves to consider both convex and nonconvex solutions of u, which are symmetric and have a single interval non-contact set. Finally the parameter ℓ = inf{x > 0 : u(x) = f (x)} characterises the length of the separation. 1.3. Variational toolkit and comments on modelling Using variational techniques developed in [1, 5], we construct a potential energy function for the system described above, and derive necessary equilibrium conditions for each of the free parameters u, ∆, θ and Q. This leads to representation of u as a four-order nonlinear free boundary problem, which once recast in intrinsic coordinates can be solved as an initial value problem with 3 unknown shooting parameters. These solutions not only have clear similarities with the geological examples presented in Fig. 1, but demonstrate some characteristic features of such nonlinear systems, infinite critical loads and localization together with interactive buckling on different scales for instance. Similar interactive behaviour has also been identified with various forms of failure in manufactured materials, for example mode interaction in the buckling of sandwich structures [9], local delamination in composites [10] and local and global mode interaction in I-beams under uniform bending [11]. A recent interdisciplinary special issue of Philosophical Transactions of the Royal Society A has been devoted to highlighting similarities and differences February 26, 2012 19:10 Philosophical Magazine 4 PhilMagFinal T. J. Dodwell and G. W. Hunt between layered structural geometry and behaviour in geology, carbon-fibre composites, and paperboard used for packaging [12]. In the work presented, we see a departure from the trend of modelling advanced geological features with large complex finite element simulations (see [13] for a recent review). As well as having clear visual appeal, such studies have also been geared towards answering a host of important geological questions. However, with the vast range of parameters often involved in such formulations, it is sometimes unclear how their geological counterparts might accurately be represented. Given a sufficiently complex setup and the ability to tweek parameters, it would seem churlish to suggest that any given result might be achieved, but there could nevertheless be some degree of subjectivity in the choices of parameter values. While we fully appreciate the value of such approaches, our objective here is to complement them with a simplified toy model in which the parameters influencing the localised bucking process can be closely defined and carefully monitored. The model is a part of a sequence of such systems [1, 5], moving towards a full description of multilayered buckling under the highly nonlinear constraints imposed by layers being required to fit together [14]. We believe that such simplifications are allowing new interpretations of the development of localized buckles in singular geometries, where voiding and related depressurization effects may be important. 2. Total potential energy and equilibrium solutions We now construct each contributing energy in turn for the system described in section 1.2, before deriving equilibrium conditions with respect to perturbations of each of the free parameters. 2.1. Total potential energy function, V The potential energy stored in the layer-parallel spring UIS , the axial springs UA , and the work done in the plastic rotational spring UKJ , follow from standard elasticity and plasticity theories and simple geometric arguments, so that 1 UIS = k1 (∆ − L(1 − cos θ))2 , 2 UA = k 2 Q2 and UKJ = 2Mp θ. (1) The elastic layer does work in bending, UB , from a flat configuration and against overburden pressure in forming a void UV . Classical bending theory [15, Ch.1] gives the bending energy over a small segment of the beam ds as dUB = 21 Bκ(s)2 ds, where κ is the curvature and B is the elastic stiffness of the beam. Integrating over all s we find that 1 UB = B 2 Z 1 κ(s) ds = B 2 X 2 Z X 1 ds u2xx dx = B 2 3 (1 + ux ) dx 2 Z X u2xx dx. (1 + u2x )5/2 If an overburden pressure q per unit length acts on the system, the work done in creating a void is given by q(u − f )dx, integrating over all x gives, UV = q Z (u − f ) dx. X If q is large then UV imposes a severe energy penalty. February 26, 2012 19:10 Philosophical Magazine PhilMagFinal 5 Convoluted accommodation structures in folded rocks We assume the elastic beam u is inextensible. This imposes a length condition on the solution u, which in turn implies compression or extension, Q, of the axial springs. Therefore the inextensibility condition can be written as J (u) = Z p X 1 + u2x dx − L = Q. (2) The layer u is under hard unilateral constraint from the arms of the hinge; therefore the contact condition u ≥ f , is imposed. The total potential energy for the system is the sum of all contributing potential energies derived above, minus the work done by the load in end shortening, P ∆, therefore 1 V (u, ∆, θ, Q) = k1 (∆ − L(1 − cos θ))2 + 2Mp θ + k2 Q2 + W (u, θ, Q) − P ∆. (3) 2 Here W (u, θ, Q) defines the potential energy functional of the layer u, given by 1 W (u, θ, Q) = B 2 Z X u2xx dx + q (1 + u2x )5/2 Z (u − f ) dx, (4) X which is constrained by the contact condition u ≥ f and the inextensibility constraint (2). 2.2. Equilibrium solutions We now derive a series of equilibrium equations, by finding stationary conditions for the potential energy with respect to perturbations of ∆, θ, Q and u. 2.2.1. Stationarity with respect to ∆ and θ Two equilibrium equations are derived by considering stationarity of V with respect to ∆ and θ, i.e. ∂V = k1 (∆ − L(1 − cos θ)) − P = 0, ∂∆ (5) ∂W ∂V = −Lk1 (∆ − L(1 − cos θ)) sin θ + 2Mp + . ∂θ ∂θ (6) and Solving these two equations simultaneously provides two solution paths. Firstly the trivial solution for which θ = 0, P = k1 ∆ (7) and secondly various non-trivial equilibrium solutions (θ 6= 0) 1 P = L sin θ  ∂W 2Mp + ∂θ  . (8) We note that ∂W/∂θ depends on the solutions profiles u, for which we now derive equilibrium conditions. February 26, 2012 19:10 Philosophical Magazine 6 PhilMagFinal T. J. Dodwell and G. W. Hunt 2.2.2. The Euler-Lagrange and Hamiltonian equations We now apply Kuhn-Tucker theory [16, pp.249] and standard Lagrange multiplier theory [16], to derive necessary conditions for a stationary solutions of V with respect to perturbations of u, subject to the contact constraint u ≥ f and the inextensibility constraint J (u) = Q. A standard variational approach considers stationarity of V with respect to vertical perturbations, i.e. u(x) → u(x)+εϕ(x), for some small constant ε and suitable test function ϕ. As a result the corresponding Euler-Lagrange equations gives the condition for a vertical force balance of the system [1]. Equally, one may insist on stationarity with respect to horizontal perturbations, i.e u(x) → u(x + εϕ). In this case the resulting differential equation is the static Hamiltonian for the system [1], and can be physically interpreted as a horizontal force balance. Firstly, applying Kuhn-Tucker theory, stationarity with respect vertical perturbations leads to the Euler-Lagrange equation  u3xx ux uxx uxxx 5 35 u3xx u2x uxxxx + − 10 − + B 2 (1 + u2x )9/2 (1 + u2x )5/2 (1 + u2x )7/2 2 (1 + u2x )7/2   u2x uxx uxx − + q = µ, (9) +λ (1 + u2x )1/2 (1 + u2x )3/2  where µ, a non negative measure, and λ ∈ R are Lagrange multipliers, arising from the contact constraint u ≥ f and the inextensibility constraint (2) respectively. Physically µ can be interpreted as the vertical component of the reactive contact load with the boundary f . Kuhn-Tucker theory represents this by the complementarity condition, Z (u − f ) dµ = 0. (10) X Therefore µ is zero wherever there is no contact, and equals a discontinuous shear force at points of delamination. The constant Lagrange multiplier λ can be interpreted as the in-plane load required to enforce the inextensibility constraint (2), and as with µ is to be found as part of the solution. There is also a first integral of (9), B   ux 5 ux uxx uxxx +λ − + qx = αδℓ . (1 + u2x )5/2 2 (1 + u2x )7/2 (1 + u2x )1/2 (11) where α=B tan θ uxxx (ℓ) +λ + qℓ, 2 5/2 (1 + tan θ) (1 + tan2 θ)1/2 which follows by imposing continuity of first and second derivatives at x = ℓ. Secondly we derive the Hamiltonian for the system by finding the stationary solutions with respect to horizontal perturbations, i.e. u(x) → u(x + εϕ). This February 26, 2012 19:10 Philosophical Magazine PhilMagFinal Convoluted accommodation structures in folded rocks 7 yields a second fourth order nonlinear differential equation,   uxxxx ux u2x uxx uxxx 5 u3xx ux 35 u3xx u3x B − 10 − + + 2 (1 + u2x )9/2 (1 + u2x )5/2 (1 + u2x )7/2 2 (1 + u2x )7/2   u3x uxx uxx ux +λ + qux = µux , (12) − (1 + u2x )1/2 (1 + u2x )3/2 which describes the horizontal force balance on the beam. This also has a first integral, given by  u2xx u2x 5 u2x u2xx 1 1 uxxx ux −λ − − +qu = qu0 − BΛ2 −λ. B 5/2 7/2 5/2 1/2 2 2 2 2 2 (1 + ux ) 2 (1 + ux ) 2 (1 + ux ) (1 + ux ) (13) Here the right hand side follows from applying the fixed symmetric boundary conditions,  u(0) = u0 , ux (0) = 0, uxx (0) = Λ and uxxx (0) = 0, (14) where u0 and Λ ∈ R are free parameters to be found as part of the solution. 2.2.3. Intrinsic representation Equations (11) and (13) can be written in terms of intrinsic coordinates, characterized by the arc length s, measured from the point of symmetry, and the tangent angle ψ with the horizontal. Recalling the relationship between both sets of coordinates [14], cos ψ = ψs = κ = uxx /(1 + u2x )3/2 , ds/dx = (1 + u2x )1/2 , dx = 1/(1 + u2x )1/2 ds sin ψ = and du = ux /(1 + u2x )1/2 ; ds (15) (16) equation (9) can be rewritten as #    " d u2xx uxx ux 1 B B +λ p + + qx (1 + u2x ) = 0, 2 dx (1 + u2x )3/2 2 (1 + u2 )3 1 + ux which is equivalent to 1 Bψss cos ψ + ( Bψs2 + λ) sin ψ + qx = 0. 2 (17) Similarly we rewrite (12), so that 1 1 Bψss sin ψ − ( Bψs2 + λ) cos ψ + qu = qu0 − BΛ2 − λ. 2 2 (18) Multiplying (17) by cos ψ and (18) by sin ψ and adding together, yields the equilibrium equation for u normal to the layer, given by 1 Bψss + qx cos ψ + qu sin ψ = (qu0 − BΛ2 − λ) sin ψ. 2 (19) February 26, 2012 19:10 Philosophical Magazine 8 PhilMagFinal T. J. Dodwell and G. W. Hunt Similarly multiplying (17) by sin ψ and (18) by cos ψ and subtracting, one obtains an equilibrium equation for u parallel to the layer, 1 1 Bψs2 + λ + qx sin ψ − qu cos ψ = (λ + BΛ2 − qu0 ) cos ψ. 2 2 (20) Substituting the free boundary condition at x = ℓ, u(ℓ) = ℓ tan θ, ux (ℓ) = tan θ and uxx (ℓ) = 0, into the above we derive the following condition 1 1 Bψs2 + λ + qℓ sin θ − qℓ tan θ cos θ = (λ + BΛ2 − qu0 ) cos θ, 2 2 (21) rearranging for u0 1 u0 = q   1 2 (1 − sec θ)λ + BΛ . 2 (22) This condition now implies an additional fixed boundary condition at x = 0. 2.3. Stationarity with respect to Q The condition for stationarity with respect to perturbations of Q, is obtained in a similar manner, so that ∂W ∂V = 2k2 Q + = 0. ∂Q ∂Q 2.4. (23) Physical interpretation of equilibrium equations The equilibrium equations give conditions under which the energy V is stationary with respect to perturbations or displacements of each of the free parameters ∆, θ, u and Q. Each of these equations describes the corresponding force or moment balance for equilibrium with respect to displacements of each of these parameters. Mathematically they represent the dual elements associated with each of these free parameters, loads with displacements and moments with rotations. The equilibrium equation P = k1 (∆ − L(1 − cos θ)), describes the force balance associated with perturbations of ∆, and is interpreted as a balance of the reactive force exerted by the horizontal spring against the in-plane force P . The Euler-Lagrange equations and Hamiltonian equation describe force balances of u vertically and horizontally respectively. We can identify each term in both of these equations as follows. Firstly for the Euler-Lagrange equation, shear force F axial load A z }| { z }| { 1 2 = (Bψs )s cos ψ + ( Bψs + λ) sin ψ + qx |{z} {z } | |2 {z } total vertical vertical component of F vertical component of A pressure αδℓ |{z} vertical contact force at delamination , (24) February 26, 2012 19:10 Philosophical Magazine PhilMagFinal 9 Convoluted accommodation structures in folded rocks and then for the Hamiltonian (Bψ ) sin ψ } | s{zs horizontal component of F 1 1 − ( Bψs2 + λ) cos ψ + q(u − u(0)) + BΛ2 + λ = {z } | 2 | |2 {z } {z } horizontal component of A total horizontal pressure horizontal load at s = 0 α tan θδℓ | {z } horizontal contact force at delamination (25) The equilibrium equation for Q is, ∂W = −2k2 Q, ∂Q Noting that Q is also written as a constraint equation on u, (2), it follows that the above equilibrium equation is equivalent to λ(Q) = −k2 Q. (26) This makes physical sense as λ(Q) is the axial load in the beam so k2 Q is the force in each spring undergoing a displacement Q. Finally the equilibrium equation (8), 1 P = L sin θ  ∂W 2Mp + ∂θ  , gives the stationary condition with respect to rotational perturbations about the hinge θ; it can therefore be interpreted as the moment balance about the hinge. Additionally ∂W/∂θ is the contribution to the total moment made by the elastic beam. This observation is important, as it allows us to derive an analytical expression for ∂W/∂θ, which substantially simplifies the following analysis. θ Figure 3. Free body diagram of whole system plus block of pressurized fluid. Figure 3 shows the moments acting about the hinge. The total moments in the anticlockwise direction are given by, 1 P L sin θ + qL2 , 2 (27) . February 26, 2012 19:10 Philosophical Magazine 10 PhilMagFinal T. J. Dodwell and G. W. Hunt whereas moments in the clockwise direction are 1 1 2Mp + M2 + P2 u0 + qL2 cos2 θ + q(L sin θ + u0 )(L sin θ − u0 ), 2 2 (28) Comparing each of these forces and moments (24) and (25), we note that P2 = 1 2 2 BΛ + λ and M2 = Bκ(x = 0) = BΛ; therefore it follows that 1 1 P L sin θ = 2Mp + BΛ + ( BΛ2 + λ)u0 − qu20 . 2 2 (29) Comparing this moment balance with equation (8), we observe that ∂W 1 1 = BΛ + ( BΛ2 + λ)u0 − qu20 . ∂θ 2 2 3. 3.1. (30) Numerical solutions Shooting for stationary solutions of u We now obtain numerical solutions of u, for a range of values of θ and physical parameters k2 , q and B. We note that for particularly convoluted solutions, ux becomes infinite, which means such solutions are unattainable in a cartesian framework. This problem is eliminated by solving the equilibrium equation normal to the layer (19). This emphasises the importance of calculating both the Euler-Lagrange and Hamiltonian equations. Numerical results can be obtained by reducing this free boundary problem to an initial value shooting problem, and by rewriting (19) as a four-dimensional first order system,    ψs ψ 1 1 2   q  ψs    =  − B (x cos ψ + w sin ψ) + B ( 2 Λ + λ − qw) sin ψ  .   x cos ψ u s sin ψ  (31) A suitable numerical integrator (e.g. MATLAB’s ode45) can then be used to shoot from the symmetric section ψ(x = 0) = 0 ψs (x = 0) = Λ and 1 u0 = q   1 2 (1 − sec θ)λ + Λ , 2 (32) towards the free boundary at x = −ℓ with boundary conditions ψ(x = −ℓ) = −θ, ψs (x = −ℓ) = 0, and u = −ℓ tan(θ). (33) so that the length constraint (2) holds. Here u0 is derived from the condition (22) and Λ, ℓ and λ are free shooting parameters to be found numerically. A three dimensional search routine was constructed [14] by nesting MATLAB’s built-in function fminsearch. This proved efficient, given a reasonable initial guess. Several different outcomes of this shooting method, from different initial guesses, are shown in Fig. 4. Unlike the earlier model [1] solutions are no longer unique, and up-buckled shapes of increasing complexity start to appear. February 26, 2012 19:10 Philosophical Magazine PhilMagFinal 11 Convoluted accommodation structures in folded rocks (a) (b) (c) Figure 4. Possible modeshapes found from (31), (32) and (33), for fixed values of q, B, θ, k1 and k2 = 2. (a) simple convex solution. (b) up-buckled solution. (c) higher-order up-buckle. 3.2. Plots of λ(Q) vs. Q The shapes of Fig. 4 are drawn at a particular value of θ and most importantly at a fixed value of Q determined by the length constraint (2). The variation of λ with Q is as shown in Fig. 5, for the convex and the first-order up-buckled profiles of Fig. 4 (a) and (b) respectively. Also shown is the straight line λ = −k2 Q. Positions a, b and c on Fig. 5 mark points where this straight line intercepts λ(Q) and so equilibrium equation (26) is satisfied. Here the load exerted by the spring exactly matches the required axial load for that solution. Profiles at each of these positions are as shown. We also show, in Fig. 5(d) a possible solution close to the limit of validity of the up-buckled shape, where the accommodating beam is about to self-intersect; we shall refer to this as the teardrop profile. The solid blue curve represents the convex solution of Fig. 4(a), and intersects the red straight line just once, at point a. There is thus a unique convex deflection profile for every value of θ. However the same is not true for the up-buckled solution, seen as the dot-dash curve in Fig. 5. For the given value of θ the straight line −k2 Q intercepts λ(Q) twice, at b and c, meaning that two different up-buckled shapes exist for the same value of θ. The shapes of the curves also suggest that, whereas the convex (solid blue) curve is always going to be crossed once by a straight line passing through the origin at Q = λ = 0, the up-buckled curve typically is crossed either twice or not at all. This is seen in the plots of Fig. 6, showing the straight and curved loci for the up-buckled shape when θ is (a) small; (b) moderate, and (c) relatively large. It is clear that up-buckled shapes only exist for θ values greater than the threshold value θmin of Fig. 6(b). This significant point is highlighted by the load-deflection plots of the next section. Once equilibrium solutions u to equation (26) are attained for each value of θ, an approximation for ∂W/∂θ can be constructed. Numerically this is done using the trapezoidal rule to calculate W for each θ, and then a finite difference approximation for ∂W/∂θ. Equation (8) can then be used to find the corresponding value of P and the following response curves constructed. 4. Overall system response The differences between the convex and up-buckled equilibrium solutions can be seen most clearly in plots of load P against end-shortening ∆, as in Fig. 7, which embrace all possible values of θ as the displacement evolves. Fig. 7(a) shows the equilibrium paths for the convex shape of Fig. 4(a). Two possible equilibrium states exist, the fundamental path seen as a red solid line where the arms are completely aligned (θ = 0), and buckled path in the bent configuration (θ 6= 0) shown as February 26, 2012 19:10 Philosophical Magazine 12 PhilMagFinal T. J. Dodwell and G. W. Hunt b. a. c. d. (a) (b) (c) (d) Figure 5. Top: numerically obtained solutions λ(Q), for simple convex solution profiles (solid blue line) and first-order up-buckled solution (dot-dash blue line). Linear relation λ = −k2 Q is also shown (straight red line). Equilibrium condition (26) holds at points a, b and c. Bottom: solution profiles at these three locations, along with d indicating the teardrop profile close to self contact. Drawn for q = B = k1 = tan θ = 1 and k2 = 50. blue/dashed. With increasing P , θ → 0 in the buckled state and these paths approach each other asymptotically, but never meet. As the load falls from this position, θ grows, until at some positive value of load the path reverses direction and follows a route with P increasing and θ continuing to grow. While falling the buckled path is unstable, but under conditions of controlled end-shortening it would restabilize at the point where ∆ is a minimum. Thus the possibility arises of a dynamic jump at constant ∆, from the fundamental state to the now-stable rising post-buckled path. To initiate this an energy hump needs to be overcome, but at high loads this would be expected to be relatively small. Because an elastic layer can only bend to finite curvature, to adopt this deflected state there must be significant slip between the layers taken up by the springs k2 . The equilibrium path for the up-buckled shape, shown in Fig. 7(b), is similar but February 26, 2012 19:10 Philosophical Magazine PhilMagFinal 13 Convoluted accommodation structures in folded rocks λ λ Q (a) Q (b) λ Q (c) Figure 6. Plots like Fig. 5 for three different values of θ. (a) No equilibrium solution at θ = 0.2 rad. (b) Limit point (one equilibrium solution) at θ = 0.67 rad. (c) Two equilibrium solutions at θ = 0.8 rad. Drawn for q = B = k1 = 1 and k2 = 50. has important differences. In particular, although it is gets close to the fundamental path over some of its length, this is only over a finite load range. The solution reaches a limit load (seen here as a cusp) at the θmin of Fig. 6(b), where it separates into the simple up-buckled shape of Fig. 5(b) and its counterpart 5(c). The simple up-buckled shape eventually restabilizes with increasing θ much like the convex mode, whereas the more open configuration is on the way to contacting with itself and becoming the teardrop of Fig. 5(d). Here the path is stopped before this state is reached, because the delamination has extended to the ends of the rigid arms. It should be noted that, although the paths of Fig. 7(a) and (b) occupy the same region of P − ∆ space, they would in fact be far from each other in the full displacement space, separated in particular by different values of Q. In the absence of any further parametric variation such as a changeable imperfection there is, we think, no likelihood of any bifurcation linking the two sets of paths. 4.1. Maxwell stability criterion For load versus corresponding deflection plots like those of Fig. 7, the total energy stored in the system can be inferred from the area under the appropriate graph [15]. The schematic of Fig. 8 shows relevant areas for three possible equilibrium states that can arise at a typical, constant, value of applied end-shortening ∆. Energy on the fundamental path, at A, simply involves the triangular area shown in 8(b). Energy at B adds to this the area E1 between the fundamental path and the high load/small θ region of the buckled path (see 8(c)), and for the energy at C area E2 is similarly subtracted, as 8(d). Note that, even when the paths approach asymptotically but never meet as in Fig. 7(a), these areas can often be found by direct integration (see for example [17]). E1 represents the input required to get February 26, 2012 19:10 Philosophical Magazine 14 PhilMagFinal T. J. Dodwell and G. W. Hunt P 150 100 50 1 2 3 4 5 ∆ 4 5 ∆ (a) P 150 100 50 1 2 3 (b) Figure 7. Equilibrium paths for the model of Fig. 2, plotted for B = 0.1, k1 = 20, k2 = 200, Mp = 20, q = 1. (a) Convex profile: the buckled (blue/dashed) solution and the fundamental (red/solid) paths approach each other asymptotically as P → ∞ but never meet. (b) First-order up-buckled profile: the buckled solution reaches a limit point (seen here as a cusp) at a finite load value where the response separates into shapes with growing and diminishing void sizes. over the energy hump at B, while E2 is the difference between this maximum and C, which under constrained end-shortening would be an energy minimum. There is of course no bifurcation point on the fundamental path, and so A is an energy minimum at all load values and therefore technically always stable. High loads however, would push the system into a dangerous state of metastability; as P increases, E1 becomes increasingly thin and shrinks in area. In such situations it is common to adopt the Maxwell stability criterion [18, p. 53] where stability rests with the global energy minimum rather than local minima. This marks out the Maxwell displacement ∆M – the point at which the energy at A equals the energy at C and E1 = E2 – as a lower bound for the point of buckling. For applied ∆ < ∆M , the global energy minimum rests with the fundamental equilibrium path: for ∆ > ∆M it rests with the buckled state. Depending perhaps on the amount of available energy from outside disturbances relative to the size of the energy hump E1 , this may or may not be a good estimate of the actual buckling load; February 26, 2012 19:10 Philosophical Magazine PhilMagFinal Convoluted accommodation structures in folded rocks 15 E1 E2 C (a) (b) (c) (d) Figure 8. (a) Typical form of load/displacement plot, P vs ∆. (b) Total energy at A. (c) Total energy at B. (d) Total energy at C. however in the absence of any other ballpark estimate it provides some kind of framework within which to assess the buckling process. The Maxwell displacement itself can rest with the convex or the up-buckled solution, depending on the choice of parameters. In practical terms the actual buckling mode is also likely to be heavily influenced by knock-on effects of initial imperfections [12]. 5. Parametric variations The plots of Figs. 7 and 8 are of the load P against its corresponding deflection (end-shortening ∆), and give a good overview of the system response in just two dimensions. Whether P or ∆ is chosen as the parameter to be controlled the same equilibrium curves are produced, although their stability is likely to depend on this choice [15]. Other parameters of the system could also be varied, and we now review the effect on the system response of two such variations. Consider for example the plots of Fig. 9, showing the response of the up-buckled shape of Fig. 4(b) at three different values at Mp . While the convex solution always takes the same form whatever the chosen Mp , this is clearly not the case here. Focussing on the limit point (seen here as a cusp) where the load is locally a maximum and θ = θmin , increasing Mp has the effect of significantly increasing its corresponding load in comparison with the remainder of the response. We note that ∂W/∂θ of equation (8) is independent of Mp , and therefore comes up with the same contribution whatever its chosen value; θmin is thus the same for all three cases. But the response for the higher values of Mp begins to look much more like that of the convex shape than for low values; the relative energy hump to be overcome to initiate buckling clearly drops as Mp is increased. We might therefore expect the Maxwell stability criterion to carry more relevance for high than for low values of February 26, 2012 19:10 Philosophical Magazine 16 PhilMagFinal T. J. Dodwell and G. W. Hunt Mp = 125 Mp = 50 Mp = 5 Figure 9. Change to up-buckled response for different values of Mp . For each, q = 0.1, B = 0.5, k1 = 30, k2 = 200 and L = 5. Also θmin = 0.17 rad at the limiting (cusp) point in each case. Mp . Much the same effect is found if k2 is changed, as shown in Fig 10. Again an k2 = 200 k2 = 400 θmin = 0.17 θmin = 0.112 k2 = 800 θmin = 0.064 Figure 10. Change to up-buckled response for different values of k2 . For each, q = 0.1, B = 0.5, k1 = 30, Mp = 5 and L = 5. Note that, unlike in Fig. 9, θmin reduces as k2 increases. increase in k2 raises the load, but this time simultaneously reduces the value of θmin at the limit point. It also causes the post-buckling curve to stiffen more quickly, quickly becoming almost parallel to the fundamental path as shown. Movement in the k2 springs is a necessary component of the up-buckled shape, equivalent of slip between the layers in the geological context, and as these springs stiffen up more work is being done in them to accommodate the same shape. For a given February 26, 2012 19:10 Philosophical Magazine PhilMagFinal Convoluted accommodation structures in folded rocks 17 θ, the higher the value of k2 the more quickly a buckling load will develop in the accommodating beam and hence the corresponding value of θmin is lowered. 6. Concluding remarks This paper presents a simplified structural model to describe the formation of convoluted accommodation structures in tight geometries that arise naturally in layered systems. It has been formulated specifically to highlight the important nonlinear geometric constraints imposed by the need for layers to fit together, and demonstrates mechanisms by which the less competent layers accommodate such tight geometries imposed by their more competent counterparts. The work extends previous work by the authors [1], to include the more general formation of the V-shaped singularity, as well as important effects of axial constraint. The pin joints at the ends of the structure, along with the introduction of springs k2 and the inextensibility, are useful simplifications but remain some way from reality. Relative slip is concentrated into deflection Q at the ends of the beam rather than being spread out over a finite, extensible, contact length as would be the case in more practical circumstances. The situation is not unlike that of a straight long pipeline on a flat bed carrying hot oil or gas, which can release heat-induced axial compression by up-buckling against gravity over a finite length [8]. Lengths of pipe at each end outside the buckled region stretch and slip along the bed to feed the buckle. In the absence of friction, these so called active lengths could be infinitely long and feed ad infinitum into the buckle. With realistic friction between the bed and the pipe this mechanism is restricted to finite lengths, over which the drop in compressive load from the unbuckled to buckled region is offset by friction from the bed. For pipelines these active lengths can be long in comparison with the buckled length, but in the geological setting friction would be expected to dominate and corresponding active lengths be relatively short. We note that the inclusion of inter layer friction is important for future models for multilayered structures; however it does offer a difficult challenge, because the normal contact forces change as the response develops. One parametric change not pursued here is that of pressure q. For the present parameter values the teardrop shape of Fig. 5(d) never fully develops, as the separated region grows to the full length of the model beforehand. However, with relative increases in q this limitation may well disappear. Pressure forces acting normal to the surface of layer should tend to reduce the separation length involved, and could perhaps even generate tension in the accommodating beam. We leave this development for future work. The paper draws on a philosophy found also in a companion paper in this volume [19], in selecting wherever possible linear elastic constitutive laws to express material behaviour. Thus a small element of the accommodating beam has linear bending stiffness B, and frictional effects are modelled by linear elastic springs k2 . The exception is the rigid-plastic rotational spring Mp . This could have been taken as a stiff linear elastic spring; bifurcation at a high rather than infinite value of P would then occur, but the post-buckling response remains much the same [14]. Linear constitutive laws may need replacing if good experimental comparisons are to be achieved, but they do have certain modelling advantages over more complex rheologies. First, they are usually more straightforward to analyse. But most importantly, they highlight the roles that nonlinear geometry and the associated potential for many alternative solutions can play in such systems. Modelling has sometimes failed in the past not because of the formulation, nor the resulting gov- February 26, 2012 19:10 Philosophical Magazine 18 PhilMagFinal REFERENCES erning equations, but because of poor choice of solution. Extra criteria for selection are sometimes necessary or useful. The buckling of the long axially compressed cylindrical shell gives a good example of just such a circumstance [20]. Because nonlinear geometry is central to the formulation, the system is capable of undergoing large deflections and rotations. Governing equations for the layer are derived in both cartesian and intrinsic coordinates, the latter meaning that for convoluted profiles as seen in Fig. 4, numerical solutions do not blow up. We again recognise the purely elastic assumptions are not realistic, but more complex rheological models might likewise be subject to error. Real geological systems would involve an intricate balance of nonlinear constitutive laws, strain history and temperature, pressure, chemical coupling etc. [21], such that precise analysis would appear an impossible dream. Real progress is perhaps best made by a combination of modelling on many levels, with proper interaction and cross-fertilization of ideas across all elements of the spectrum. References [1] T.J. Dodwell, M.A. Peletier, C.J. Budd and G.W. Hunt, SIAM J. Appl. Math. (SIAP) 72 (2012) p. 444–463. [2] M.A. Biot, Mechanics of Incremental Deformations, Wiley, New York, 1965. [3] N.J. Price and J.W. Cosgrove, Analysis of geological structures, 1st Edition Cambridge University Press, Cambridge, UK, 1990. [4] J.A. Boon, C. Budd and G.W. Hunt, Proc R Soc A 463 (2007) p. 1447–1466. [5] T.J. Dodwell, G. Hunt, M.A. Peletier and C.J. Budd, Multilayered folding with voids, (2012), for Geometry and Mechanics of Layered Structures and Materials, Special Issue of Phil. Trans. R. Soc. Lond. A (to appear). [6] J.G. Ramsay, Geol. Soc. Am. Bul. 85 (1974) p. 1741–54. [7] B.E. Hobbs, W.D. Means and P.F. Williams, An outline of structural geology, Wiley, New York, 1976. [8] R.E. Hobbs, ASCE J. Transp. Eng. 110 (1984) p. 175–189. [9] G.W. Hunt and M.A. Wadee, Proc. R. Soc. Lond. A 454 (1998) p. 1187–1216. [10] R. Butler, A.T. Rhead, W. Liu and N. Kontis, Compressive strength of delaminated aerospace composites, (2012), for Geometry and Mechanics of Layered Structures and Materials, Special Issue of Phil. Trans. R. Soc. Lond. A (to appear). [11] M.A. Wadee and L. Gardner, Proc. R. Soc. Lond. A 468 (2012) p. 245–268. [12] C.J. Budd, R. Butler and G.W. Hunt (eds.), Geometry and Mechanics of Layered Structures and Materials, Preface to Special Issue of Phil. Trans. R. Soc. Lond. A (to appear), 2012. [13] P.J. Hudleston and S.H. Treagus, J. Struct. Geol. 32 (2010) p. 2042–2071. [14] T.J. Dodwell, Multilayered folding with constraints, PhD Thesis, University of Bath, UK, 2011. [15] J.M.T. Thompson and G.W. Hunt, A General Theory of Elastic Stability, Wiley, London, 1973. [16] D.G. Luenberger, Optimization by Vector Space Methods, Wiley-Interscience, 1968. [17] G.W. Hunt, M.A. Peletier and M.A. Wadee, J. Struct. Geol. 22 (2000) p. 669–681. [18] E.C. Zeeman, Catastrophe Theory: Selected papers, 1972-1977, Addison-Wesley, 1977. [19] Giles Hunt and Jill Hammond, Philos. Mag. (2012) (this special issue). [20] G.W. Hunt, IMA J. Appl. Math. 76 (2011) p. 2–26. [21] B.E. Hobbs, K. Regenauer-Lieb and A. Ord, J. Struct. Geol. 30 (2008) p. 1572–1592.