Orientational order in buckling elastic membranes

Physica D 205 (2005) 267–274

Orientational order in buckling elastic membranes Nariya Uchida∗ Department of Physics, Tohoku University, Sendai 980-8578, Japan Available online 10 February 2005

Abstract We study pattern formation in buckling membranes under lateral compression, using the F¨oppl–von K´arm´an (FvK) approximation for nonlinear elasticity. To identify the role of shear strain, we numerically compare the patterns on elastic and fluid membranes as well as their coarsening dynamics. By mapping the FvK equations onto a vector spin model with nonlocal interactions, we find anisotropic correlation of ridges that characterizes the topography of crumpled elastic membranes. © 2005 Elsevier B.V. All rights reserved. PACS: 46.25+x; 05.45−a; 68.55−a Keywords: Nonlinear elasticity; Pattern formation; Phase ordering

1. Introduction Since Euler, instability of thin elastic bodies has been a topic of long-lasting interest, both from academic and practical points of view. Straight rods and plates compressed along their axes undergo a buckling instability at critical load. While the equilibrium shape of buckled rods are well known as elastica, the patterns on buckled plates or membranes, as seen in crumpled papers and crushed cans, are far more complicated. Recently, some progress has been made in the analysis of singular structures in a strongly deformed membrane. Witten and coworkers [1] have developed a scaling theory for the geometry of a ridge, where two slopes meet and the bending energy concentrates. Moldovan and Golubovi´c [2] performed an overdamped molecular dynamics simulation and found power-law type growth of buckling patterns, which supports the scaling theory. A ridge terminates at the peak of a cone, which is another kind of singularity in a crumpled membrane. The geometry and dynamics of a developable cone, created by a point force acting on the membrane, have been studied in detail [3–5]. However, still little is known about the interactions between ridges and/or cones, which control the topography of crumpled membranes. The difficulty lies ∗

Tel.: +81 22 217 7756; fax; +81 22 217 7756. E-mail address: [email protected] (N. Uchida).

0167-2789/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.physd.2005.01.009

268

N. Uchida / Physica D 205 (2005) 267–274

in the analysis of nonlinear elasticity involving coupled fourth-order differential equations, known as the F¨oppl–von K´arm´an (FvK) equations [6]. For instance, it is only recently that secondary (undulatory) instabilities of a straight ridge under nearly uni-directional compression are classified [7]. In this article, we attempt a new approach seeking an analogy with orientational phase ordering. The translational invariance of our system naturally lead us to regard the gradient of the membrane height as the order parameter. By projecting out the in-plane displacement from the FvK elastic energy, we obtain effective interactions between the local gradient vectors. Of our particular interest is the role of shear strain, which we discuss by comparing the patterns on elastic and fluid membranes. The organization of this paper is as follows. In Section 2, we briefly review the FvK approximation for large deformation of membranes. In Section 3, we describe a numerical simulation of buckling dynamics and interpret the results within the scaling theory. In Section 4, we introduce a mapping onto a vector spin model, with which we characterize the topography of elastic membranes. We conclude in Section 5.

2. The F¨oppl–von K´arm´an approximation The FvK approximation for deformation of an initially flat membrane assumes the following two conditions: (i) the thickness of the membrane is negligibly small in comparison to the characteristic lengthscale of the deformation; this allows us to treat the membrane as a curved surface of vanishing thickness; (ii) the gradient m = ∇h of the membrane, h being its height from a reference plane, is much smaller than unity in magnitude. Let us take the Cartesian coordinate system in which the undeformed membrane lies in the xy plane. The membrane’s configuration is specified by the three-dimensional position R = R(x, y) of each material element whose initial position is R0 = (r, 0) = (x, y, 0). We decompose the displacement into horizontal and vertical components as R − R0 = (u, h) = (ux , uy , h). The FvK elastic energy consists of strain and bending contributions as F = Fs + Fb . The strain part is given by  
λ 2 Fs = dx dy (1) Eii + µEij2 , 2 where λ and µ are the Lam´e constants and Eij = 21 [∂i uj + ∂j ui + (∂i h)(∂j h)]

(i, j = x, y)

(2)

is the strain tensor linearized with respect to the in-plane displacement u [6]. The linearization is justified by the assumption m = |∇h|  1, because the in-plane strain Uij = 21 (∂i uj + ∂j ui ) is expected to be an O(m2 ) quantity in mechanical equilibrium. The bending part has the simple form
κ Fb = dx dy (∇ 2 h)2 , 2

(3)

(4)

where κ is the bending modulus and can be expressed in terms of λ, µ, and the membrane’s thickness [6]. We also introduce the 2D bulk elastic modulus K = λ + µ for later use. We control the externally applied strain Uij (bars denote spatial averages throughout this paper), and assume periodic boundary conditions for δu = u − U · r and h. 3. Numerical simulation of buckling dynamics Now let us assume that the membrane is suddenly compressed at t = 0 after which the applied strain Uij (<0) is kept constant. The simplest dynamical equation to model the resultant buckling process is the overdamped

N. Uchida / Physica D 205 (2005) 267–274

269

one, ∂R δF = −Γ , ∂t δR

(5)

where Γ is the inverse of the friction constant. This can be a realistic model for membranes in a liquid environment [8]. We choose this form not for its reality, but because it produces dynamical scaling properties from which we can extract useful information on the mechanics of structural elements. The above continuum model was implemented on a 512 × 512 square lattice and we chose (K, µ, κ, Γ ) = (2, 1, 0.01, 0.01) as the standard parameter set, with the time step t = 1 and grid size x = 1. Isotropic compression with Uxx = Uyy = −0.01 was imposed at t = 0 and then evolution of the patterns were followed by integrating Eq. (5) with the Euler scheme. For comparison, we also simulated the buckling of a fluid membrane with µ = 0 and the other parameters unchanged. In Fig. 1, we show snapshots of the height field h(r), mean curvature ∇ 2 h, and square gradient (∇h)2 for both elastic and fluid membranes. For the elastic case, the plot of (∇h)2 clearly shows ridges as black lines, while their junctions correspond to cones. The cones have large absolute values of the mean curvature and appear as bright and dark spots in the plot of ∇ 2 h. The fluid membrane shows a very different pattern composed of nearly isotropic mounds and basins, the top and bottom (resp.) of which appear as dark spots in the plot of (∇h)2 . Coarsening of the patterns is monitored through the bending energy Fb  and square height h2 , plotted versus time in Fig. 2. We also show the orientational correlation length ξm , defined as the first zero of the orientational correlation function: ˆ + r) · m(r), ˆ Gm (|x|) = m(x

ˆ = m

m ∇h = . |m| |∇h|

(6)

Evolution of these quantities are excellently fitted by the power laws Fb  ∝ t −α , h2  ∝ t β , and ξm ∝ t γ , with (α, β, γ) being (0.47, 0.54, 0.27) for the elastic membrane and (0.47, 0.45, 0.24) for the fluid one. The former set of exponents agrees well with the previous scaling argument [2], according to which (α, β, γ) = (5/11, 6/11, 3/11) (0.45, 0.54, 0.27). On the other hand, the exponents for fluid membranes can be understood by a reduced dynamical equation which we will derive in Section 4. After projecting out the horizontal displacement u as a fast variable, we will have the dynamical equation for m in the form:   ∂m = Γ Uii + 21 m2 ∇ 2 m − Γκ∇ 2 ∇ 2 m. ∂t

(7)

We solved it numerically, and the resultant growth curves are also plotted in Fig. 2. In the late stage, they are in good agreement with those from the original model. Noting that the equation is locally linear (i.e., the average m2 −4 in the late is insensitive to changes in m at a given point) and that the two terms in the RHS should scale as ξm stage, we can apply a simple dimensional analysis to obtain (α, β, γ) = (1/2, 1/2, 1/4). This is fairly close to the numerical result, and is also in consistent with the fact that buckling rods in viscous media has a lengthscale that grows as t 1/4 [9]. A more direct check of the dynamical scaling comes from the plots of the orientational correlation function, given in Fig. 3. A marked difference between the elastic and fluid membranes is the oscillatory correlation in the latter. It has at least four minima in Gm (x), only two of which are within the range of plot (Fig. 3(b)).

4. Effective interactions and orientational correlation Having seen the differences between buckling patterns of elastic and fluid membranes, we are in a position to analyze the role of shear rigidity. We proceed by eliminating the in-plane displacement u from the original FvK theory, assuming the mechanical equilibrium condition δF/δu = 0. The resultant effective strain energy depends on the height field h only through its gradient, which is a natural consequence of the translational invariance of our

270

N. Uchida / Physica D 205 (2005) 267–274

Fig. 1. Snapshots of the height h, mean curvature ∇ 2 h, and square gradient (∇h)2 , for (a) an elastic membrane (K = 2, µ = 1) and (b) a fluid membrane (K = 2, µ = 0). In the greyscale, brighter colors are used for larger values; dark colors correspond to small positive values of (∇h)2 and negative values of h and ∇ 2 h.

system. To be explicit, the force balance equation can be written as 0=−

δF = (λ + µ)∇(∇ · u) + µ∇ 2 u + g, δu

g = λ2 ∇(m2 ) + µ∇ · (mm).

(8) (9)

From the divergence of Eq. (8), we obtain ∇ ·u=−

1 g + const, (λ + 2µ) ∇ 2

(10)

N. Uchida / Physica D 205 (2005) 267–274

271

Fig. 2. (a) Bending energy, (b) square height, and (c) orientational correlation length vs. time. Solid lines show the fitting by power laws.

with which the in-plane displacement is determined up to a linear function of r. After some calculation, we arrive at the expression for the strain tensor: Eij = Pij + Cij ,

(11)

Fig. 3. Orientational correlation function Gm (x) for (a) elastic and (b) fluid membranes. The distance is scaled by the correlation length ξm .

272

N. Uchida / Physica D 205 (2005) 267–274

where Cij is spatially constant and    1 λ + µ ∇ i ∇j ∇ k ∇l 2 Pij = δij − m − m k ml , 2 λ + 2µ ∇ 2 ∇2

(12)

where 1/∇ 2 is the inverse operator of Laplacian. Noting that Eij = Uij + mi mj /2 + δPij with δPij = Pij − Pij , we obtain the effective strain energy in the form:       
λ 1 1 (13) Uii + m2 m2 + δPii2 + µ Uij + mi mj mi mj + δPij2 + const. Fs = dr 2 4 4 For fluid membranes (µ = 0), Pij does not contribute to the effective energy because Pii = 0. The reduced dynamical Eq. (7) can be derived by converting the original equation ∂h/∂t = −ΓδF/δh into the form: δH ∂m = Γ ∇2 , ∂t δm where H = H[m] is the effective Hamiltonian given by   
1 1 dr λ Uii + m2 m2 + κ(∇m)2 . H= 2 2

(14)

(15)

This resembles the Hamiltonian for the 2D XY magnet but an important difference is the quartic term now replaced by the global coupling ∝ m2 m2 .

Fig. 4. Snapshots of the orientational order parameter field Qij (r) for (a) elastic and (b) fluid membranes. In (a), ridges (where Qij 0) are seen as grey lines.

N. Uchida / Physica D 205 (2005) 267–274

273

Shear rigidity gives rise to anisotropic long-range interactions between the local gradient vectors. Since Pij contains the gradient only through the combination mi mj , it is useful to introduce the orientational order parameter tensor: 1 Qij = mi mj − m2 δij . 2

(16)

It is represented by the unitary matrix:

 Q cos 2θ Q sin 2θ , {Qij } = Q sin 2θ −Q cos 2θ

(17)

where θ = arctan(my /mx ) and Q = m2 /2. Reexpressing the effective elastic energy (13) in terms of the Fourier components of Qij and Q, we obtain terms that depend on the direction of the wavevector. They are combined into the form:
dq µ(λ + µ) q |ˆqi qˆ j Qij − Qq |2 , (18) Faniso = 2(λ + 2µ) q (2π)2 where qˆ = q/|q| gives the angle φ of the wavevector as qˆ = (cos φ, sin φ). This long-range interaction q q q q q q breaks the symmetry between the variables Q1 = 2ˆqi qˆ j Qij = cos(2φ)Qxx + sin(2φ)Qxy and Q2 = sin(2φ)Qxx − q q q q cos(2φ)Qxy , which appear symmetrically in local terms such as |Qij |2 = |Q1 |2 + |Q2 |2 . In real space, Faniso causes difference in the gradients of Qxx in the x- and y-directions and those in the oblique directions. The simulated patterns in Fig. 4(a) indeed shows that ∇Qxx favors the oblique directions. Note that ridges run in perpendicular to m and gradients of Qij correspond to deflection of ridges. The observed anisotropy means that the correlation of ridge directions are weaker in the parallel and perpendicular directions than in the oblique directions. For comparison, we give the Qij fields of a fluid membrane in Fig. 4(b). They show that the correlation anisotropy is opposite to that of elastic membranes. Its origin is not apparent from the effective free energy (13), and should be sought for in the constraint ∂x my = ∂y mx (= ∂x ∂y h). This constraint can be approximately imposed by adding a free energy term fν = ν(∂x my − ∂y mx )2 , and is made rigid in the limit ν → ∞. The additional free energy suppresses orientational correlation parallel to the local gradient m. This is seen by temporarily assuming |m| = 1 and rewriting fν as ν(m · ∇θ)2 ; spatial variation in |m| would not qualitatively change the anisotropy. In order to quantify the correlation anisotropy, now we introduce the relative orientational correlation function: Grel (X, Y ) = Qij (r)Qij (r + R),

R = (X cos θ, Y sin θ),

(19)

where θ is the angle of m as above. This function on the X- and Y-axes gives correlation in the directions parallel and perpendicular to the local gradient, respectively. In Fig. 5, we plot the function for elastic and fluid membranes.

Fig. 5. Relative orientational correlation function Grel (X, Y ) for (a) elastic and (b) fluid membranes. For each, statistical average is taken over 10 snapshots at t = 160, 000. The vertical axis is rescaled such that Grel (0) = 1.

274

N. Uchida / Physica D 205 (2005) 267–274

For elastic membranes, we see that the correlation is suppressed on the X- and Y-axes while it is enhanced in the oblique directions, except near the origin where two satellite peaks appear on the X-axis. On the other hand, Grel for fluid membranes is oscillatory in every direction but negative peaks are found only in the oblique directions. These clearly show that the shear-mediated long-range interaction Faniso changes the correlation anisotropy. It should be mentioned that the anisotropy of fluid membranes remains for when the ratio µ/λ is small. The crossover from the fluid-like to elastic behaviors as well as the effect of anisotropic external strain will be discussed elsewhere.

5. Summary To summarize, we have characterized the patterns on buckling elastic membranes through mapping the FvK theory onto a model of orientational ordering. The nonlocal interactions mediated by in-plane shear strain cause anisotropic correlation of ridges. It may worth mention that similar nonlocal interactions arise from strain–orientation couplings in some soft matter complexes, such as liquid–crystalline elastomers, membranes, and crosslinked polymer blends [10]. Here the concept has been applied to elastic membranes under the assumption of weak nonlinearity. It remains for future work to extend the analysis to a broader class of crumpling phenomena that involve large and often irreversible deformations. This work is supported by the Grant-in-Aid for Scientific Research from Japan Ministry of Education, Culture, Sports, Science, and Technology.

References [1] T.A. Witten, H. Li, Europhys. Lett. 23 (1993) 51; A.E. Lobkovsky, et al., Science 270 (1995) 1482; E.M. Kramer, T.A. Witten, Phys. Rev. Lett. 78 (1997) 1303; B.A. Didonna, T.A. Witten, Phys. Rev. Lett. 87 (2001) 206105. [2] D. Moldovan, L. Golubovi´c, Phys. Rev. Lett. 82 (1999) 2884; D. Moldovan, L. Golubovi´c, Phys. Rev. E 60 (1999) 4377. [3] M. Ben Amar, Y. Pomeau, Proc. R. Soc. London A 453 (1997) 729. [4] S. Chaieb, F. Melo, Phys. Rev. Lett. 80 (1998) 2354; E. Cerda, L. Mahadevan, Phys. Rev. Lett. 80 (1998) 2358; E. Cerda, et al., Nature 401 (1999) 46. [5] A. Boudaoud, et al., Nature 407 (2000) 718. [6] L.D. Landau, E.M. Lifshitz, Theory of Elasticity, Pergamon, Oxford, 1986. [7] B. Audoly, B. Roman, A. Poceau, Eur. Phys. J. B 27 (2002) 7. [8] L. Bordieu, et al., Phys. Rev. Lett. 72 (1994) 1502; A. Saint-Jalmes, F. Gallet, Eur. Phys. J. B 2 (1998) 498. [9] L. Golubovi´c, D. Moldovan, A. Peredera, Phys. Rev. Lett. 81 (1998) 3387; O. Hallatschek, E. Frey, K. Kroy, cond-mat/0402082 (preprint). [10] P.D. Olmsted, J. Phys. II 4 (1994) 2215; N. Uchida, A. Onuki, Europhys. Lett. 45 (1999) 341; N. Uchida, Phys. Rev. E 66 (2002) 040902; N. Uchida, J. Phys.: Condens. Matter 16 (2004) L21.