Straightening wrinkles

Journal of the Mechanics and Physics of Solids 65 (2014) 1–11 Contents lists available at ScienceDirect Journal of the Mechanics and Physics of Solids journal homepage: www.elsevier.com/locate/jmps Straightening wrinkles M. Destrade a,b,n, R.W. Ogden c, I. Sgura d, L. Vergori a,c a School of Mathematics, Statistics & Applied Mathematics, National University of Ireland, Galway, Ireland School of Mechanical & Materials Engineering, University College Dublin, Belfield, Dublin 4, Ireland c School of Mathematics & Statistics, University of Glasgow, UK d Dipartimento di Matematica e Fisica, Università del Salento, Lecce, Italy b a r t i c l e i n f o abstract Article history: Received 4 September 2013 Received in revised form 21 December 2013 Accepted 3 January 2014 Available online 15 January 2014 We consider the elastic deformation of a circular cylindrical sector composed of an incompressible isotropic soft solid when it is straightened into a rectangular block. In this process, the circumferential line elements on the original inner face of the sector are stretched while those on the original outer face are contracted. We investigate the geometrical and physical conditions under which the latter line elements can be contracted to the point where a localized incremental instability develops. We provide a robust algorithm to solve the corresponding two-point boundary value problem, which is stiff numerically. We illustrate the results with full incremental displacement fields in the case of Mooney–Rivlin materials and also perform an asymptotic analysis for thin sectors. & 2014 Elsevier Ltd. All rights reserved. Keywords: Soft solids Straightening Instability Asymptotic analysis Impedance Matrix method 1. Introduction The back of the elbow on an extended arm and the region of contact between the road and a tyre are two examples of a straightening deformation. In this paper we revisit and take further the analysis of this somewhat neglected exact solution of nonlinear elasticity. We extend it to include a study of incremental instability: in the case of the elbow, instability could model the appearance of wrinkles; in that of the tyre, it would correspond to the onset of the so-called Schallamach (1971) stripes; see Fig. 1 for photographs. There are only six known families of large deformations which are universal to all nonlinear incompressible isotropic materials (see for instance the textbook Tadmor et al., 2012 for a recent exposé). They are       Family Family Family Family Family Family 0: 1: 2: 3: 4: 5: Homogeneous deformations. Bending, stretching and shearing of a rectangular block. Straightening, stretching and shearing of a sector of a hollow cylinder. Inflation, bending, torsion, extension and shearing of an annular wedge. Inflation or eversion of a sector of a spherical shell. Inflation, bending, extension and azimuthal shearing of an annular wedge. There are countless stability studies for Families 0, 1, 3, and 4; however, there exists no work on the stability of a straightened sector. In fact, the literature on the large straightening deformation itself is very sparse and we have been able n Corresponding author at: School of Mathematics, Statistics & Applied Mathematics, National University of Ireland, Galway, Ireland. E-mail address: michel.destrade@nuigalway.ie (M. Destrade). 0022-5096/$ - see front matter & 2014 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jmps.2014.01.001 2 M. Destrade et al. / J. Mech. Phys. Solids 65 (2014) 1–11 Fig. 1. Straightening instabilities. (a) As the arm extends, wrinkles may appear on the surface of the elbow, the number of which seems to depend mostly on the age of the subject. (b) Single frame of a film of the contact between perspex and a butyl sphere sliding over it at 0.043 cm/s, showing the appearance of instability stripes (taken from Schallamach, 1971. to identify only five contributions on the subject (Ericksen, 1954; Truesdell and Noll, 2004; Aron et al., 1998; Aron, 2000, 2005). Although in practice, bending (Family 1) and straightening (Family 2) are the opposite of each other, their respective theoretical modelling differs completely. For instance, we shall see in this paper that the large straightening deformation occurs in plane strain and in plane stress (λ3 ¼ 1 and s1 ¼ 0 throughout the block), while large bending is accompanied by plane strain only (the normal stress is zero only point-wise on the bent faces). That no result at all can be deduced from one deformation to apply to the other is particularly true for the stability analysis, which proves to be extremely difficult to conduct. Both problems can be formulated in terms of a linear ordinary differential system with variable coefficients which is stiff numerically. For the bending problem, this numerical stiffness can be smoothed out by the Compound Matrix method, which has proved to be most successful in the past in stability analyses for Family 0 (e.g. Haughton, 2011), Family 1 (e.g. Haughton, 1999; Coman and Destrade, 2008; Roccabianca et al., 2011), Family 3 (e.g. Destrade et al., 2010) and Family 4 (e.g. Fu, 1998; Fu and Lin, 2002). Here we tried several numerical methods in turn for the straightening problem: Determinantal method, Compound Matrix method, Surface Impedance method, and only the latter one turned out to be robust enough to handle this very stiff problem (see also Destrade et al., 2009, 2014). In his seminal paper on “Deformations Possible in Every Isotropic, Incompressible, Perfectly Elastic Body”, Ericksen (1954) showed that a circular sector can always be “straightened” into a rectangular block, by mapping the reference cylindrical polar coordinates ðR; Θ; ZÞ of the material points in the region 0 o R1 r Rr R2 ; Θ0 rΘ rΘ0 ; 0 r Z rH; ð1Þ to the rectangular Cartesian coordinates ðx1 ; x2 ; x3 Þ of points in the region a r x1 r b; l r x2 r l; 0 r x3 r H; ð2Þ in the deformed configuration through the deformation x1 ¼ 1 2 AR ; 2 x2 ¼ Θ ; A x3 ¼ Z; ð3Þ see Fig. 2 and also Truesdell and Noll (2004). Here 0 o2Θ0 o2π is the angle spanning the solid sector, and 2π 2Θ0 is what we will call the “open angle” of the original, undeformed, sector (in contrast to the term “opening angle” of a deformed sector used, for example, in Destrade et al., 2010). Also, R ¼ R1 and R ¼ R2 are the inner and outer faces of the open sector, respectively, while x1 ¼ a and x1 ¼ b are their counterparts in the deformed rectangular block. The quantity A, which has the dimensions of the inverse of a length, is to be determined from the boundary conditions. In general, its value depends on the loads imposed on the end faces x2 ¼ 7 l of the straightened block to sustain the deformation. Examples of such loads (in the absence of body forces) include a resultant normal force N, a resultant moment M, or a mixture of the two. In this paper we take the point of view that the final length 2l of the block is found from the stability analysis and we calculate the corresponding A, N and M required to attain it. In particular, it is clear from (3) that A ¼ Θ0 =l; ð4Þ and that the faces normal to the x1 axis are at x1 ¼ a ¼ 1 2 Θ0 2 AR ¼ R ; 2 1 2l 1 x1 ¼ b ¼ 1 2 Θ0 2 AR ¼ R : 2 2 2l 2 ð5Þ Section 2 covers the characteristics of the large straightening deformation. In particular, we show there that the principal Cauchy stress component s1 is zero throughout the block. 3 M. Destrade et al. / J. Mech. Phys. Solids 65 (2014) 1–11 2l H H b-a R2 R1 x1 x3 x2 Fig. 2. Straightening of a circular cylindrical sector with initial radii ratio R1 =R2 ¼ 0:2 and open angle 4π=3. The line elements along x2 (originally circumferential) on the inner straightened face x1 ¼ a (originally at R ¼ R1 ) are extended; those on the outer face x1 ¼ b (originally at R ¼ R2 ) are contracted. Eventually, at a critical stretch of contraction, the outer face buckles and wrinkles appear. Fig. 3. When a sector made of Mooney–Rivlin material with angle Θ0 ¼ 2π=3 and radii ratio R1 =R2 ¼ 0:1 is straightened, it buckles with the formation of k ¼ 3 wrinkles, with amplitude decaying rapidly away from the face x1 ¼ b. Prescription of the final length 2l will prove to be most useful in the subsequent instability analysis (Section 3), especially for making connections with known results of surface instability. Indeed, here the straightening instability threshold corresponds to the possible existence of static wrinkles on the face (5)2 of the straight block, see Fig. 3 for a 3D representation of such wrinkles. Consider the deformations on these surfaces: the face R ¼ R2 in the circumferential direction is deformed into the face x1 ¼ b, and this face is compressed in the x2 direction if R2 Θ0 4 l. This will be the case provided any axial tension applied to the faces x2 ¼ 7l is not too large. Similarly, the face R ¼ R1 will be stretched by the deformation if R1 Θ0 o l. This suggests, in particular, that on the face x1 ¼ b there should exist a critical compressive stretch, λcr say, at which wrinkles with sinusoidal variations in the x2-direction should appear, and we shall use λcr as our bifurcation parameter. The corresponding bifurcation criterion should be dispersive because of the presence of inherent characteristic lengths and, in particular, in the thick block/ small wavelength limit it should yield the surface stability criterion of plane strain compression. For instance, for a block made of neo-Hookean, Mooney–Rivlin, or generic third-order incompressible elastic material, the asymptotic critical stretch should be λcr ¼ 0:544, as established by Biot (1963). We conduct the full numerical analysis of straightening wrinkles in Section 5, where we treat the example of the Mooney–Rivlin material (which here coincides with the neo-Hookean material since the deformation is restricted to plane strain). In Section 4, we provide an asymptotic stability analysis for thin sectors of a general strain-energy density and establish the expression for the critical compressive stretch up to second order in the thickness to radius ratio of the sector. 4 M. Destrade et al. / J. Mech. Phys. Solids 65 (2014) 1–11 2. Large straightening deformation We compute the deformation gradient F from (3) as F ¼ ARe1 ER þ 1 e2 AR EΘ þe3 EZ ; ð6Þ where ER ; EΘ ; EZ and e1 ; e2 ; e3 are the cylindrical polar and Cartesian unit basis vectors in the reference and deformed configurations, respectively. It follows that the Eulerian principal directions of the deformation (defined as the directions of the eigenvectors of FFT ) are the Cartesian basis vectors and that the principal stretches (the square roots of the eigenvalues of FFT ) are λ1 ¼ AR; λ2 ¼ 1 ; AR λ3 ¼ 1; ð7Þ showing that the deformation (3) is plane strain and isochoric. In this paper, we consider homogeneous, incompressible isotropic hyperelastic materials, for which the straightening deformation is universal (Ericksen, 1954). We denote by W ¼ Wðλ1 ; λ2 ; λ3 Þ the strain energy per unit volume, so that the Cauchy stress tensor is r ¼ s1 e 1 e1 þ s2 e2 e2 þ s3 e3 e3 ; ð8Þ where s1 ; s2 ; s3 are the principal Cauchy stresses, given by s1 ¼ λ1 ∂W ∂λ1 p; s2 ¼ λ2 ∂W ∂λ2 p; s3 ¼ ∂W ∂λ3 p; ð9Þ and p is a Lagrange multiplier associated with the constraint of incompressibility (det F ¼ λ1 λ2 λ3  1). From now on it is convenient to introduce the notation λ1 ¼ λ 1 ; λ2 ¼ λ; ð10Þ and to define the strain energy in terms of a single deformation variable, the circumferential stretch λ. Thus, we define the ^ ¼W ^ ðλÞ by function W ^ ðλÞ ¼ Wðλ W 1 ; λ; 1Þ; ð11Þ and it then follows that s2 ^ 0 ðλÞ; s1 ¼ λW ð12Þ where the prime denotes the derivative with respect to λ. We take the undeformed configuration to be stress free, so that ^ 0 ð1Þ ¼ 0: W ð13Þ Clearly the deformation depends only on the single variable R (or equivalently x1) and hence the second and third components of the equilibrium equation div r ¼ 0 in the absence of body forces show that p is independent of x2 and x3, while the first component yields simply ds1 =dx1 ¼ 0. Hence the normal principal stress component s1 is constant throughout the block. Next we assume that the boundaries x1 ¼ a; b are free of traction, and we deduce that 0 ^ ðλÞ; s2 ¼ λW s1  0; ð14Þ while s3 has to be calculated from (9)3 with p ¼ λ1 ∂W=∂λ1 . We require the stress s2 to be positive when corresponding to stretching and negative when corresponding to contraction, and we therefore impose the conditions ^ 0 ðλÞ⪌0 W according as λ⪌1: ð15Þ ^ is a convex function of λ, and we assume that this is the case henceforth. In view of (13), these conditions hold if W Now we calculate the resultant normal force N and moment M (about the origin of the Cartesian coordinate system) on a face x2 ¼ constant of the block. They are independent of x2 and are given by Z b Z b N¼H s2 dx1 ; M ¼ H s2 x1 dx1 : ð16Þ a a On use of the connections λ ¼ 1=ðARÞ and x1 ¼ AR2 =2, it is straightforward to show that Z Z λa H λa H ^ 0 ðλÞ dλ; M ¼ ^ 0 ðλÞ dλ; N¼ λ 2W λ 4W A λb 2A2 λb ð17Þ where λa ¼ 1 1 4 λb ¼ ; AR1 AR2 ð18Þ M. Destrade et al. / J. Mech. Phys. Solids 65 (2014) 1–11 5 are the stretches in the x2-direction on the faces x1 ¼ a and x1 ¼ b of the straightened block, respectively. In general, the resultant force N is not zero. Note that the possibility of straightening a circular cylindrical sector by applying a moment alone and no normal force is treated in Destrade et al. (2014). 3. Small wrinkles In order to study the (linearized) stability of the deformed rectangular configuration, we investigate the possible existence of infinitesimal static solutions in the neighbourhood of the large straightening deformation. Hence we consider superimposed displacements u ¼ uðx1 ; x2 ; x3 Þ. We let L ¼ grad u denote the displacement gradient, so that the incremental incompressibility condition reads tr L  Lii ¼ ui;i ¼ 0 ð19Þ in the usual summation convention for repeated indices, where a subscript i following a comma signifies differentiation with respect to xi. The incremental nominal stress s_ 0 is given by (Ogden, 1997) s_ 0 ¼ A 0 L þ pL _ pI; ð20Þ where a superposed dot signifies an increment, the zero subscript indicates evaluation in the deformed configuration, I is the identity tensor, and A 0 is the fourth-order tensor of instantaneous elastic moduli. In components s_ 0ij ¼ A0ijkl ul;k þ pui;j _ ij ; pδ ð21Þ where δij is the Kronecker delta. Referred to the (Eulerian) principal axes of the underlying deformation, the only non-zero components of A 0 are A0iijj ¼ λi λj W ij ; i; j ¼ 1; 2; 3;  A0ijij ¼ A0ijji þλi W i ¼ λi W i λj W j  λ2i λ2i λ2j ; i aj; λi a λj ; ð22Þ where W i  ∂W=∂λi , W ij  ∂2 W=∂λi ∂λj (note the major symmetry A0piqj ¼ A0qjpi ). The corresponding formulas for λj ¼ λi ; i a j, in (22)2 may be obtained by a limiting procedure, but are not needed here. In the absence of body forces, the incremental equation of equilibrium has the form div s_ 0 ¼ 0: ð23Þ Here, for simplicity, we focus on plane two-dimensional incremental deformations such that u3 ¼ 0 and u1, u2 are independent of x3. Further, we seek solutions which have sinusoidal variations in the x2-direction, i.e. of the form _ ¼ fU 1 ðx1 Þ; U 2 ðx1 Þ; Pðx1 Þgeinx2 ; fu1 ; u2 ; pg ð24Þ where n ¼ kπA=Θ0 ¼ kπ=l; ð25Þ is the wavenumber, and the integer k is the mode number, identifying the number of wrinkles on the faces x1 ¼ constant, with amplitude decaying from x1 ¼ b to x1 ¼ a. Clearly, the wavelength L of the wrinkles is L ¼ 2π=n ¼ 2l=k; ð26Þ i.e. it is equal to the length of the deformed block divided by the number of wrinkles. These two quantities are to be determined from the numerical treatment of the stability analysis. It follows from (21) that the components of the incremental nominal stress tensor have a similar form s_ 0ij ¼ Sij ðx1 Þeinx2 ; i; j ¼ 1; 2; ð27Þ say. Then, by using a standard procedure (see, e.g., Destrade and Scott, 2004; Destrade and Ogden, 2005; Destrade et al., 2009, 2010), we cast the above equations as a first-order differential system for the four-component displacement–traction vector η  ½U 1 ; U 2 ; iS11 ; iS12 ŠT . This is the so-called Stroh formulation: ! G1 ðx1 Þ G2 ðx1 Þ d ηðx1 Þ ¼ i Gðx1 Þηðx1 Þ ¼ i ð28Þ ηðx1 Þ: G3 ðx1 Þ G1 ðx1 Þ dx1 Here the real matrix G has the 0 0 n B 0 B n G¼B B n2 s2 0 @ ^ ″ðλÞ 0 n2 λ2 W form 1 0 0 0 0 C 1=α C C; n C A n 0 ð29Þ 6 M. Destrade et al. / J. Mech. Phys. Solids 65 (2014) 1–11 where ^ 0 ðλÞ=ðλ4 α ¼ A01212 ¼ λW 1Þ; ð30Þ and its 2  2 sub-matrices G1 ; G2 ; G3 are symmetric. Now, if the incremental equations of equilibrium can be solved, subject to the traction-free conditions S11 ¼ S12 ¼ 0 on x1 ¼ a; b; ð31Þ then possible equilibrium states exist in a neighbourhood of the straightened configuration, signalling the onset of instability. When such a configuration is determined, we refer to the value of the circumferential stretch λb ¼ 1=ðAR2 Þ at this point as the critical value for compression, which we denote by λcr . Simple calculations show that Aðb aÞ ¼ R22 R21 2 2 2λcr R2 ; ð32Þ which allows for the complete determination of the straightened geometry just prior to instability. In particular, the dimensions of the rectangular block in the x1 and x2 directions are then, respectively, b a¼ R22 R21 ; 2λcr R2 l ¼ R2 Θ0 λcr : ð33Þ This latter equation, coupled to (26), shows that the critical wavelength of bifurcation is the length of the external curved boundary of the undeformed sector, multiplied by the critical stretch ratio and divided by the number of wrinkles. Take, for instance, the sector in Fig. 2(a). To draw it we scaled lengths with respect to R2 and picked R1 ¼0.2, R2 ¼ 1, Θ0 ¼ π=3. The stability analysis of Section 5 will reveal that if the material is modelled by the Mooney–Rivlin potential, then the corresponding buckled state occurs at λcr ¼ 0:5710, with k ¼1 wrinkle forming over the surface; see Fig. 4. It follows from (33) that just prior to instability, at λb ¼ 0:5711 say, the dimensions of the rectangular block are b a ¼ 0:840 and l¼0.598; see Fig. 2(b). Then A can be found from λcr , and N and M from (17). By (26), the critical wavelength is L ¼ 1:20. 4. Asymptotic analysis for thin sectors In this section we derive the approximation to the critical threshold λcr when the undeformed cylindrical sector is thin compared to the radius of the outer face, that is when ɛ1 R1 51: R2 ð34Þ To conduct a perturbation analysis of the incremental equations of equilibrium (28)–(31) in terms of the small parameter ɛ, we need to expand the strain energy in terms of the strain. All expansions are equivalent, and here we choose to use the Green–Lagrange strain tensor E ¼ ðFT F IÞ=2. For incompressible isotropic elastic solids, the expansion begins as   A   tr E3 þ DðtrðE2 ÞÞ2 þ ⋯; ð35Þ W ¼ μ tr E2 þ 3 where μ, A, D, the respective second-, third-, and fourth-order elastic constants, are all of the same order of magnitude and of the same dimensions (see Destrade and Ogden, 2010 and references therein). For the straightening deformation (6), indeed for any plane strain deformation, we obtain for future reference the results ^ 0 ð1Þ ¼ 0; W ^ ″ð1Þ ¼ 4μ; W ^ ‴ð1Þ ¼ W ^ ⁗ð1Þ ¼ 156μþ 48A þ96D; W 12μ; αð1Þ ¼ μ; α0 ð1Þ ¼ 2μ; ð36Þ where α ¼ αðλÞ is defined by (30). We now introduce the following dimensionless quantities: sffiffiffiffiffiffiffiffiffiffiffiffiffi! 1 2x1 λcr kπ α 1 ; α¼ ; ξ¼ ; n¼ ɛ Θ0 μ R2 s2 ¼ s2 ; μ Ui ¼ Ui ; R2 S 1i ¼ S1i ; μ W¼ ^ W : μ ð37Þ The space variable is thus ξ, which has the advantage of spanning the thickness of the block with the fixed range ½0; 1Š. In terms of the new variable ξ, the stretch λ reads as λ¼ 1 λcr : ɛξ ð38Þ M. Destrade et al. / J. Mech. Phys. Solids 65 (2014) 1–11 7 Fig. 4. Some examples of the critical value of the stretch λcr as a function of the initial radii ratio R1 =R2 for a Mooney–Rivlin material. The number of wrinkles k depends on the angle Θ0 and on R1 =R2 ; see Table 1 for details. For instance, for Θ0 ¼ 2π=3 (curve 2) the number of wrinkles is k ¼3, 2, or 1, depending on the value of R1 =R2 . For Θ0 ¼ π=3; π=4; π=5; π=6 (curves 4; 5; 6; 7) there is only one wrinkle, for any radii ratio. Substituting (37) into Eqs. (28)–(29) yields dη ¼ dξ iɛλ 1 Gη; ð39Þ where η ¼ ½U 1 ; U 2 ; iS 11 ; iS 12 ŠT and the dimensionless Stroh matrix is 0 B n=λ cr B G ¼B B n 2 s 2 =λ2 cr @ 0 0 0 n=λcr 0 0 0 0 n 2 W ″=ð1 ɛξÞ2 n=λcr 0 1 1=α C C C: n=λcr C A 0 ð40Þ Finally, we introduce the power series expansions of λcr and of η in ɛ: 1 λcr ¼ ∑ ɛ i λðiÞ cr ; i¼0 1 η ¼ ∑ ɛi ηðiÞ : ð41Þ i¼0 The dimensionless displacement–traction vector η must satisfy the boundary conditions (31) at all orders, so that ðiÞ ηðiÞ 3 ¼ η4 ¼ 0; for ξ ¼ 0; 1: ð42Þ 8 M. Destrade et al. / J. Mech. Phys. Solids 65 (2014) 1–11 We start the asymptotic analysis at order Oð1Þ, where Eq. (39) gives dηð0Þ ¼ 0; dξ ð43Þ for which the solutions satisfying (42) are of the form ηð0Þ ¼ ½a0 ; b0 ; 0; 0ŠT ; ð44Þ where a0 and b0 are integration constants, to be determined at the next order of the asymptotic analysis. At order OðɛÞ, Eq. (39) gives dηð1Þ ¼ dξ 1 iðλð0Þ cr Þ Gð0Þ ηð0Þ ; ð45Þ with 0 ð0Þ n=λcr 0 B ð0Þ B n=λcr B Gð0Þ ¼ B 2 0 ð0Þ B n W =λcr @ 0 0 0 0 0 0 n2W ″ 0 ð0Þ n=λcr 0 1 C 1=α C C C; n=λð0Þ cr C A 0 ð46Þ 0 ð0Þ where α, W and W ″ are evaluated at λcr (note that we have used the connection (14)2 to express s 2 in terms of W ). By direct integration of (45) we obtain   1 a1 ðn=λð0Þ η1ð1Þ ¼ iðλð0Þ cr Þ cr Þb0 ξ ;   1 ð0Þ b1 ðn=λcr Þa0 ξ ; η2ð1Þ ¼ iðλð0Þ cr Þ 1 0 ðn 2 W =λð0Þ cr Þa0 ξ; η3ð1Þ ¼ iðλð0Þ cr Þ η4ð1Þ ¼ iðn 2 W ″Þb0 ξ; ð47Þ η1ð1Þ where a1 and b1 are constants of integration for and respectively, while the constants of integration for and ηð1Þ 4 were taken as zero in order to satisfy (42) at ξ ¼ 0. To satisfy the boundary condition at ξ ¼ 1 we must have, simultaneously, 0 ð0Þ feither W ðλð0Þ cr Þ ¼ 0 or a0 ¼ 0g and feither W ″ðλcr Þ ¼ 0 or b0 ¼ 0g. As a consequence of (15) and of the convexity of W , the ð0Þ ð0Þ unique value of λcr satisfying this is λcr ¼ 1, which is in line with our expectation that sectors of vanishing thickness should buckle as soon as deformed (see Coman and Destrade, 2008 or Goriely et al., 2008 for similar conclusions for the bending or compression of thin elements). Using (36), we thus obtain the following reduced expressions for the displacement–traction vector: ηð0Þ ¼ ½a0 ; 0; 0; 0ŠT ; and for the matrix: 0 0 B n B ð0Þ G ¼B @ 0 0 n 0 0 4n 2 ηð1Þ ¼ 0 0 0 n i½a1 ; b1 0 η2ð1Þ , na0 ξ; 0; 0ŠT ; 1 1C C C: nA ηð1Þ 3 ð48Þ ð49Þ 0 Moving on to the collection of the terms of order Oðɛ2 Þ in Eq. (39), we obtain h i  ð1Þ  dηð2Þ ¼ i Gð1Þ ηð0Þ þ Gð0Þ ηð1Þ λcr þ ξ Gð0Þ ηð0Þ ; dξ where the matrix Gð0Þ is given by (49) from now on, and 0 1 ð1Þ 0 nλcr 0 0 B C ð1Þ nλcr 0 0 2ðλð1Þ B cr þ ξÞ C B C ð1Þ Gð1Þ ¼ B 4n 2 ðλð1Þ þ ξÞ C; 0 0 nλcr B C cr   @ A ð1Þ ð1Þ 2 0 4n ξ þ 3λcr nλcr 0 ð50Þ ð51Þ where again, we have used (36). We then integrate these equations directly and apply the boundary conditions to deduce ð1Þ λcr , and repeat the process at the next order. The results are obtained simply as  1 1  2 ð2Þ ; λcr 2n λð1Þ ¼ 3 : ð52Þ cr ¼ 2 24 Here we see that the correction at order one to the critical stretch is the same for any cylindrical sector irrespective of its open angle, and is valid for any value of the mode number k. The effects of the open angle and of the mode number on the critical threshold λcr reveal themselves at order two in the thickness, although that correction accounts only for geometrical M. Destrade et al. / J. Mech. Phys. Solids 65 (2014) 1–11 9 ð2Þ effects since no physical elastic constant is present. Clearly, the largest value of λcr corresponds to n ¼ π=Θ0 , that is when k ¼1 and only one vertical wrinkle occurs on the outer face of the straightened block. We have thus established the following universal asymptotic expansion of the critical compression value of the circumferential stretch: !   1 1 π2 ɛ 2 2 3 ɛ2 þ O ɛ3 : λcr ¼ 1 ð53Þ 2 24 Θ0 5. General stability analysis: numerics The inhomogeneous differential system (28) is numerically stiff and calls for a robust algorithm. Here, we favour the Impedance Matrix method (see Destrade et al., 2009 and Norris and Shuvalov, 2012 for successful implementation of this algorithm in solids with cylindrical symmetry, and see Shuvalov (2003) for a rigorous and exhaustive treatment of the underlying theory). To summarise, the conditional impedance matrix Za is zero on the “inner” face x1 ¼ a and singular on the “outer” face x1 ¼ b: det Za ðbÞ ¼ 0: Za ðaÞ ¼ 0; ð54Þ By definition it relates the displacement U to the traction S through Sðx1 Þ ¼ Za ðx1 ÞUðx1 Þ; ð55Þ and by construction it satisfies the matrix differential Riccati equation dZa ¼ iðG1 Za dx1 Za G1 Þ þ Za G2 Za þG3 ; ð56Þ where the G0 s are 2  2 sub-blocks of the Stroh matrix (29). To solve numerically the boundary value problem given by (56) with boundary conditions (54) and to identify the critical value λcr , we apply a shooting-like technique that combines the bisection (dichotomy) method with an initial value solver for (56) with (54)1 . Hence, for a given value of λcr in some starting guess interval, we integrate (56) subject to (54)1 by using the ode45 function of Matlab (based on Runge–Kutta methods) and, if the target value (54)2 is not met at x1 ¼ b, then by dichotomy programming we adjust λcr until the target is satisfied with a given precision (tol ¼ 10 7 ). Then, once the critical value of the circumferential stretch is found, we have SðbÞ ¼ Za ðbÞUðbÞ ¼ 0, which means that U 2 ðbÞ ¼ U 1 ðbÞ Z a11 ðbÞ ¼ Z a12 ðbÞ Z a21 ðbÞ : Z a22 ðbÞ ð57Þ This ratio determines the shape of the wrinkles on the outer face of the straightened block. Moreover, the displacement U satisfies the differential equation dU ¼ iG1 U dx1 G2 Za U; ð58Þ and we may thus use (57) as the starting point for its backward (from face x1 ¼ b to face x1 ¼ a) numerical integration in order to compute the entire displacement field (a similar equation also exists for the traction field). For our numerical experiments, we first non-dimensionalise the equations (note that this is a different nondimensionalisation process from that conducted in the asymptotic analysis of the previous section). Specifically, we use the dimensionless quantities x1 h 2 2 i kπ U b xn1 ¼ A R1 =R2 ; 1 ; nn ¼ ; U ni ¼ i ; Zna ¼ Za ; b b 2Θ0 μ ^ S α s W 1j 2 n n ; S1j ¼ ; ð59Þ ; sn2 ¼ αn ¼ ; W ¼ μ μ μ μ ^ ″ð1Þ=4 is the infinitesimal shear modulus. Note that by using (3), (4) and λcr , it becomes a simple matter to where μ ¼ W establish the useful links xn1 ¼ λ 2 2 λcr ; nn ¼ bλ2cr n: ð60Þ Then, substitution of (59) into Eq. (28) yields a starred version of the governing equations (28), (56), (57) and (58), where the non-dimensionalised Stroh matrix Gn is made up of the blocks ! ! 0 0 0 nn λcr 2 n n ; G2 ¼ μG2 ¼ ; G1 ¼ bG1 ¼ 0 1=αn nn λcr 2 0 10 M. Destrade et al. / J. Mech. Phys. Solids 65 (2014) 1–11 Table 1 Description of the results displayed in Fig. 4: For each numbered curve, the corresponding angle and number of wrinkles is given. Curve number Angle Number of wrinkles 1 Θ0 ¼ π 2 Θ0 ¼ 2π=3 3 Θ0 ¼ π=2 4 5 6 7 Θ0 ¼ π=3 Θ0 ¼ π=4 Θ0 ¼ π=5 Θ0 ¼ π=6 k ¼4 k ¼1 k ¼3 k ¼2 k ¼1 k ¼2 k ¼1 k ¼1 k ¼1 k ¼1 k ¼1 2 b G3 ¼ G3 ¼ μ n nn2 λcr 4 sn2 0 0 nn2 λcr 4 λ2 W n″ ! ; for for for for for for for for for for for 0 r R1 =R2 r 0:15 0:15 r R1 =R2 r 1 0 r R1 =R2 r 0:11 0:11 r R1 =R2 r 0:17 0:17 r R1 =R2 r 1 0 r R1 =R2 r 0:18 0:18 r R1 =R2 r 1 0 r R1 =R2 r 1 0 r R1 =R2 r 1 0 r R1 =R2 r 1 0 r R1 =R2 r 1 ð61Þ and the boundary conditions (54) are replaced by Zna ðR21 =R22 Þ ¼ 0; det Zna ð1Þ ¼ 0: ð62Þ Next, we focus on a specific material for the constitutive modelling, namely the Mooney–Rivlin material, with strain energy density: W ¼ C 1 ½trðFFT Þ 3Š þ C 2 ½trðFFT Þ 1 3Š; ð63Þ where C1, C2 are positive constants, with C 1 þ C 2 ¼ 2μ. Then we find that αn ¼ λcr 2 xn1 ; sn2 ¼ λ2cr xn1 xn1 λ2cr ; λ 2 W n″ ¼ xn λ2cr þ 3 21 ; xn1 λcr ð64Þ and the non-dimensional equations can therefore be solved independent of the physical constants C1 and C2 (see Destrade et al., 2009 for a study of the influence of constitutive parameters on the stability of bent blocks). The quantities left at our disposal are geometric: the angle Θ0 and k, the number of wrinkles on the face x1 ¼ b (i.e. xn1 ¼ 1). Each choice of these quantities determines nn and then, for each initial thickness ratio R1 =R2 , we seek the corresponding critical stretch λcr . For practical purposes here, we choose a given angle Θ0 and vary k to ensure that the value of λcr corresponds to the earliest onset of buckling (i.e. λcr is as close to 1 as possible). Take for instance a sector with angle Θ0 ¼ π=2. We find that for 0 rR1 =R2 r0:1833, the straightened block buckles with k¼2 wrinkles, and for 0:1833 rR1 =R2 r1, it buckles with k ¼1 wrinkle. Several different examples are listed in Table 1 and displayed in Fig. 4. In the thick sector (i.e. R1 =R2 close to 0), small wavelength (i.e. Θ0 close to 0) limit, all curves tend to λcr ¼ 0:544, the value found by Biot (1963) for surface instability in plane strain of a Mooney–Rivlin half-space (note that for plane strain the Mooney–Rivlin energy function reduces to the neo-Hookean one). In the thin sector limit, all curves are approximated by the quadratic (53). The Impedance Matrix method has the advantage that it provides not only the critical compressive stretch but also the full incremental displacement field as obtained by solving (58). Hence we have obtained the exact incremental fields shown in Figs. 3 (in 3D) and 4 (in 2D). Note that the incremental analysis gives access to the amplitude of the incremental displacements up to an arbitrary scalar factor, and that we have chosen an exaggerated scale to display them. In the examples shown in Fig. 4, we see that the displacement fields for the modes k ¼ 2; 3; 4 are highly localized near the surface of the compressed face, their amplitudes decaying rapidly with depth. For the mode k ¼1 (only one wrinkle), the displacements are not as localized, and propagate somewhat through the thickness of the block. This phenomenon is particularly visible in the Θ0 ¼ 2π=3, R1 =R2 ¼ 0:2 case. Acknowledgments Partial funding from the Royal Society of London (International Joint Project grant for MD, RWO, LV), from the Istituto Nazionale di Alta Matematica (Marie Curie COFUND Fellowship for LV; Visiting Professor Scheme for MD, IS), and from the Italian Ministry of Education, Universities and Research (PRIN-2009 project “Matematica e meccanica dei sistemi biologici e dei tessuti molli” for IS) is gratefully acknowledged. We are also indebted to Jerry Murphy (Dublin City University) for helpful discussions and to Artur Gower (National University of Ireland Galway) for technical support. M. Destrade et al. / J. Mech. Phys. Solids 65 (2014) 1–11 11 References Aron, M., Christopher, C., Wang, Y., 1998. On the straightening of compressible, nonlinearly elastic, annular cylindrical sectors. Math. Mech. Solids 3, 131–145. Aron, M., 2000. Some remarks concerning a boundary-value problem in non-linear elastostatics. J. Elast. 60, 165–172. Aron, M., 2005. Combined axial shearing, extension, and straightening of elastic annular cylindrical sectors. IMA J. Appl. Math. 70, 53–63. Biot, A.M., 1963. Surface instability of rubber in compression. Appl. Sci. Res. A 12, 168–182. Coman, C., Destrade, M., 2008. Asymptotic results for bifurcations in pure bending of rubber blocks. Q. J. Mech. Appl. Math. 61, 395–414. Destrade, M., Ni Annaidh, A., Coman, C.D., 2009. Bending instabilities of soft biological tissues. Int. J. Solids Struct. 46, 4322–4330. Destrade, M., Ogden, R.W., 2005. Surface waves in a stretched and sheared incompressible elastic material. Int. J. Non-Linear Mech. 40, 241–253. Destrade, M., Ogden, R.W., 2010. On the third- and fourth-order constants of incompressible isotropic elasticity. J. Acoust. Soc. Am. 128, 3334–3343. Destrade, M., Ogden, R.W., Sgura, I., Vergori, L., 2014. Straightening: existence, uniqueness and stability. Proc. R. Soc. Lond. A, in press. Destrade, M., Murphy, J.G., Ogden, R.W., 2010. On deforming a sector of a circular cylindrical tube into an intact tube: existence, uniqueness, and stability. Int. J. Eng. Sci. 48, 1212–1224. Destrade, M., Scott, N.H., 2004. Surface waves in a deformed isotropic hyperelastic material subject to an isotropic internal constraint. Wave Motion 40, 347–357. Ericksen, J.L., 1954. Deformations possible in every isotropic, incompressible, perfectly elastic body. Z. Angew. Math. Phys. 5, 466–489. Fu, Y.B., 1998. Some asymptotic results concerning the buckling of a spherical shell of arbitrary thickness. Int. J. Non-Linear Mech. 33, 1111–1122. Fu, Y.B., Lin, Y.P., 2002. A WKB analysis of the buckling of an everted neo-Hookean cylindrical tube. Math. Mech. Solids 7, 483–501. Goriely, A., Vandiver, R., Destrade, M., 2008. Nonlinear Euler buckling. Proc. R. Soc. Lond. Ser. A 464, 3003–3019. Haughton, D.M., 1999. Flexure and compression of incompressible elastic plates. Int. J. Eng. Sci. 37, 1693–1708. Haughton, D.M., 2011. A practical method for the evaluation of eigenfunctions from compound matrix variables in finite elastic bifurcation problems. Int. J. Non-Linear Mech. 46, 795–799. Norris, A.N., Shuvalov, A.L., 2012. Elastodynamics of radially inhomogeneous spherically anisotropic elastic materials in the Stroh formalism. Proc. R. Soc. Lond. Ser. A 468, 467–484. Ogden, R.W., 1997. Nonlinear Elastic Deformations. Dover, New York. Roccabianca, S., Bigoni, D., Gei, M., 2011. Long-wavelength bifurcations and multiple neutral axes in elastic multilayers subject to finite bending. J. Mech. Mater. Struct. 6, 511–512. Schallamach, A., 1971. How does rubber slide? Wear 17 (191). Shuvalov, A.L., 2003. A sextic formalism for three-dimensional elastodynamics of cylindrically anisotropic radially inhomogeneous materials. Proc. R. Soc. Lond. Ser. A 459, 1611–1639. Tadmor, E.B., Miller, R.E., Elliott, R.S., 2012. Continuum Mechanics and Thermodynamics. University Press, Cambridge. Truesdell, C., Noll, W., 2004. The Non-linear Field Theories of Mechanics. Springer, Berlin.