Lagrangian Crumpling Equations 9 0 0 Mark A. Peterson 2 Mount Holyoke College n a January 27, 2009 J 7 2 Abstract ] t A concise method for following the evolving geometry of a moving f o surface using Lagrangian coordinates is described. All computations can s be done in the fixed geometry of the initial surface despite the evolving . complexity of the moving surface. The method is applied to three prob- t a lems in nonlinear elasticity: the bulging of a thin plate under pressure m (theoriginal motivationforF¨oppl-vonKarmantheory),thebucklingofa - spherical shell underpressure, and the phenomenon of capillary wrinkles d induced by surface tension in a thin film. In this last problem the inclu- n sion ofagravitational potentialenergyterm inthetotalenergyimproves o the agreement with experiment. c [ 1 1 Introduction v 6 The elasticity theory of thin shells is largely differential geometry by another 6 name. Inthis paper I describe a method for followingthe differentialgeometric 1 4 dataofasurfaceasitmoves,andillustrateitsapplicationtonon-linearelasticity . theory. Theequationsofthemethodarecompletelygeneralforsmoothsurfaces, 1 0 and so could in principle describe the complex motions of crumpling up to the 9 formation of singularities. 0 Problemsinvolvingelasticmembraneshavebeenapproachedinseveralways, v: including numericalsimulationby triangulatedsurfaces,using a polyhedralap- i proximation to differential geometry [1]. Another approach has been to use X differential geometry and scaling laws to understand the line and point singu- r larities of crumpled surfaces analytically [2, 3, 4, 5, 6, 7], and numerically[8]. a The method of this paper generalizes familiar methods of mechanical engineer- ing for the non-linear elasticity theory of thin shells [9, 10, 11] in going beyond second order, and in treating initially curved surfaces in a unified way. Section 2 establishes notation for the differential geometry of a moving sur- face and shows how to use Lagrangian coordinates to simplify its description, themainideaofthispaper. Section3summarizestheobservationsoftheprevi- ous section in a system of differential equations for the evolving surface and its strains. Section 4 compares this approachto F¨oppl-von Karman(FvK) theory, 1 andsolvesthe motivating problemfor thattheory,the bulging ofathin rectan- gular plate subject to pressure, by integrating the evolution equations forward in time. Section 5 uses second order expansions of the crumpling equations to describe the buckling of a sphere under pressure. Section 6 uses the insights of CerdaandMahadevan[12]togiveamoredetaileddescriptionofaphenomenon recentlydiscussedin[13],capillarywrinklesinducedbysurfacetensioninathin film. A previously unnoticed discrepancy with experiment is partially resolved with the inclusion of the gravitationalpotential energy of the system. 2 Geometrical Methods In terms of smooth coordinates (x1,x2,x3) in space one can describe the de- formation of a material object by the trajectories of its constituent particles, solutions of equations of motion dxi =Vi(x,t) (1) dt where V =Vi∂ (2) i isthe vectorfieldgeneratingthe flow,andtisaparameteralongthe flow. Inte- gratingthe systemforwardtot=1,onecanalsothink ofVi asadisplacement, a slight abuse of notation that should be clear from context. Metric relations among the particles are given by ds2 =g dxidxj (3) ij where g is a Riemannian metric tensor, perhaps, but not necessarily, the Eu- ij clidean metric. I coordinatize the material object by Lagrangiancoordinates, convected by the flow, i.e., every material point keeps the same coordinates that it had orig- inally. In this case the changing metric relationship of material points, namely the changein the expressionEq.(3), is due entirely to the changein the metric components g , because dxi, which for this purpose simply assigns to a line ij segment the coordinate difference of its endpoints, is invariant. The rate of changeas a consequenceofthis deformationinthe components ofthe metric g, or of any second rank tensor G, expressed in convectedcoordinates, is givenby the Lie derivative[14][15] £ G(∂ ,∂ )=VG +G([∂ ,V],∂ )+G(∂ ,[∂ ,V]) (4) V j k jk j k j k Here [ , ] is the Lie bracket of vector fields. It is more common to express objectslikethis,derivativesoftensorswhicharethemselvestensors,intermsof the covariant derivative with respect to the metric connection, and to employ the conventions of raising and lowering indices with g and its matrix inverse ij gij, such that, for example the covector with components V =g Vj (5) i ij 2 has covariantderivative with respect to xk (denoted V ), in terms of the ordi- i;k nary partial derivative (denoted V ) given by i,k V =V +Γj V (6) i;k i,k ik j where the coefficients of connection Γ are 1 Γj = gjm(g g g ) (7) ik 2 ik,m− mi,k− km,i It is straightforwardto verify for any second rank tensor G that µν £ G(∂ ,∂ )=gij(V G +V G +V G ) (8) V k ℓ j kl;i j;k iℓ j;ℓ ki In particular, if G is the metric tensor g, which is a covariant constant, we recover the well known result £ g(∂ ,∂ )=V +V =2U (9) V k ℓ ℓ;k k;ℓ kℓ where U is the rate of strain tensor of the flow V (or the first order strain of the displacement V). Nothing said above was specific to three dimensions, and therefore every statement can be interpreted as referring to a surface with a Riemannian structure if the indices take only two values and not three. From now on I shall use Latin indices for tensors in three-space, and Greek indices for tensors on a surface. Now consider a smooth material surface M, so thin that one may regard it as 2-dimensional,andlet (x2,x3) be coordinatesin this surface, while x1 =z is displacementalongthenormaltothesurface,withthepositivedirectionchosen conventionally, such that the surface M is z = 0. Such a coordinate system exists for a neighborhood of M such that z <1/C, where C is the supremum | | over M of both principal curvatures in absolute value. The metric tensor in these coordinates takes the form 1 0 g = . (10) 0 g +2zh +z2k µν µν µν (cid:18) (cid:19) The tensor g , with Greek indices taking values (2,3), is the first fundamental µν form of M, h is the second fundamental form, and k = hλh is the third µν µν µ λν fundamental form. All these tensors are associated with the surface M, and not with the ambient space. They do not depend on z, i.e. all z dependence in Eq. (10) is explicit. The plus sign on the middle term is a conventional choice. On a sphere, for example, one could take the positive direction for z to be the outernormaldirection,andtheprincipalcurvaturesofthespheretobepositive. Now let a vector field (a,Vµ) be prescribed on M with normal component a(x2,x3) and tangential components Vµ(x2,x3), and extend it to a neighbor- hood of M as W =a∂ +Vµ∂ zG′µνa ∂ , (11) z µ ,µ ν − whereinitiallythetensorG′µν =gµν. Inashorttime ∆t,the flowgeneratedby the velocity field W changes the metric tensor components by approximately ∆g =∆t£ g (12) W 3 The tensor g+∆g regarded as a tensor on 3-space expresses the ambient Eu- clideangeometryinLagrangiancoordinates. Ifg+∆gisrestrictedtothesurface z = 0 and indices (2,3), one has the slightly altered first fundamental form of M G =(g +∆g ) , (13) µν µν µν |z=0 expressingthe non-Euclideangeometry ofthe slightly alteredM induced by its embedding in the ambient Euclidean space. The term linear in z in Eq. (11) was chosen to maintain the block diagonal form of Eq. (10) to first order in z. Therefore, taking the z-derivative, one has the slightly altered second funda- mental form of M, ∂ H = (g +∆g ) , (14) µν µν µν ∂z (cid:20) (cid:21)|z=0 The third fundamentalform couldnotbe computed in this way,but it is deter- mined by H , µν K =H Hλ . (15) µν µλ ν I now imagine taking a sequence of such small steps, and I will continue to denote by G andH the evolvingfirstandsecondfundamentalforms giving µν µν the Riemannian structure on M induced by the embedding in Euclidean space. I will not make use of this Riemannian structure for computations, however. There is anothernaturalRiemannianstructureonM, namely thatgivenby theoriginal,undeformedfirstfundamentalformg ,togetherwithitsassociated µν connection, etc., which I shall continue to use, being careful not to give it erro- neous interpretations. This Riemannian structure, unlike G , has no obvious µν geometrical meaning on the deformed surface, but it is still useful in a formal way. Another possible interpretation, deliberately suppressing the geometrical meaningofG ,is toimagineasurfacethatis notdeformedbythe flowW but µν carries tensor fields G and H , initially coinciding with g and h , that µν µν µν µν aredeformed by W. That these happen to be the first andsecondfundamental formsofanevolvingsurfaceisforgotten. Inthispicturetheundeformedg has µν an obvious geometrical meaning as the metric on the underlying, unchanging surface which is the arena for the evolving G and H . µν µν In Eq. (11) I introduced the tensor G′µν, initially gµν. More generally G′µν istheinverseofG asamatrix. ItisatensorfieldonM,butitisnotobtained µν from G by raising indices. Raising indices is an operation accomplished by µν gµν, my chosen Riemannian structure, not by G′µν. The prime on G′ is a reminder that it is not some version of the tensor G. I have shown how G and H change, to first order, under a deformation µν µν (a,Vµ) of M, assumed now always to be extended off M as in Eq (11). In turn, (a,Vµ) might evolve so as to reduce at each step a free energy functional depending on G and H . In this way I will arrive at crumpling equations, µν µν a system of differential equations for (a,Vµ), G , and H , describing the µν µν evolution of M. Before considering the equation for (a,Vµ), though, there is another issue to consider. 4 This formulation leaves implicit what the evolving surface actually looks like, since mere knowledge of G and H is not a convenient description of µν µν M. To keeptrackofthe positions ofpoints onthe surface,one shouldintegrate Eq.(1)usingcomponentsofW(0,xµ)=(a,Vµ)withrespecttofixedCartesian coordinate axes. Let XA(x1,x2,x3,t) be a Cartesian coordinate function in space. It is time independent in the physical sense, but its functional form depends on time because the xi evolve in time. The 1-form dXA = XA dxj ,j assigns the XA component WA to the vector W. This 1-form evolves in time at the rate given by the Lie derivative £ dXA(∂ ) = W dXA(∂ )+dXA([∂ ,W]) (16) W i i i = WjXA +Wj XA (17) ,ij ,i ,j = (XA Wj) (18) ,j ,i Thus UA = UiXA, the XA-component of any vector field U = Ui∂ at time t, i i can be found using XA(x1,x2,x3,t) solving i ∂XA ∂(XAWj) i = j (19) ∂t ∂xi with appropriate initial conditions. By the definition of the coordinate x1 =z, the Cartesian coordinate XA is an affine linear function of z. It is essential thereforetoexpandXAWj onlytofirstorderinz inEq.(19). Tobecompletely j explicit, XA is independent of z and we can represent 1 XA =YA(x2,x3)+zZA(x2,x3). (20) µ µ µ Then Eq. (19) says ∂X1A = ZAVµ YAa G′µν (21) ∂t µ − µ ,ν ∂YA µ = (XAa+XAVν) (22) ∂t 1 ν ,µ ∂ZA µ = (ZAVν YAa G′νλ) (23) ∂t ν − ν ,λ ,µ The linear approximation I have made in the neighborhood of M obscures the fact that if W were made to carry affine normal lines to affine normal lines exactly, as one could always require by a suitable nonlinear extension W of (a,Vµ) off M, then XAWj would be exactly an affine linear function of z j without approximation. The evolution of M is the same for any extension, however, so what looks like a linear approximation in the method is actually exact. As a special case, I describe motion at constant velocity, i.e., ∂WA/∂t = 0 for each component A. Then ∂ XAWj ∂Wj 0= j =(XAWk) Wj +XA (24) ∂t k ,j j ∂t (cid:0) (cid:1) 5 Thus the components of W must evolve according to ∂Wk = Xk(XAWj) Wi (25) ∂t − A ,j ,i Here Xk is the inverseofXA, consideredas a matrix. Eq.(25)for straightline A j motion is recognizable as ∂Wk +Wj Wk =0 (26) j ∂t ∇ where isthecovariantderivativewithrespecttothemetricconnectionofthe k ∇ Euclidean metric in 3-space expressed in the evolving Lagrangian coordinates. I emphasize that I have chosen, however, not to use the evolving geometry but rather the fixed initial geometry of M for all computations, a great simplifica- tion. 3 Evolution Equations By the arguments of the previous section the surface M evolves according to ∂G κλ = VµG +Vµ G +Vµ G +2aH (27) ∂t κλ;µ ;κ µλ ;λ µκ κλ ∂Hκλ = aK a + 1a G′µν( G +G +G ) κλ ,λ;κ ,µ κλ;ν νλ;κ νκ;λ ∂t − 2 − +VµH +Vµ H +Vµ H (28) κλ;µ ;κ µλ ;λ µκ (29) Using these relations one can find how other geometric quantities change, for example the area element √Gdx2dx3, involving the determinant of the first fundamental form G=G G G G (30) 22 33 23 32 − The result is ∂√G =(Vµ√G) +aG′µνH √G (31) ,µ µν ∂t Integrating one finds √G and hence dilation strain. The strain tensor 1 (G g ) (32) µν µν 2 − can be found by integrating Eq. (27). A natural definition for nonlinear shear strain S is µν ∂S 1 ∂G 1 ∂√G µν µν = G . (33) µν ∂t 2 ∂t − √G ∂t ! The subtracted term removes the contribution of dilation strain. S is not µν traceless, in general, beyond first order. 6 4 Comparison with F¨oppl-von Karman Approach A simple example illustrates the use of this formalism and points out its rela- tionshiptoF¨oppl-vonKarman(FvK)theory[9]. FvKconsiderstheequilibrium state of a thin membrane subject to external forces and boundary conditions. Since the metric strain within a membrane is typically small, even for large normal displacements, it makes sense to continue to use linear stress-strain re- lationships. The strain may, however, be a nonlinear function of displacement, and therefore displacement may be nonlinearly related to stress. FvK thus produces nonlinear equations for the equilibrium shape of an elastic membrane subject to external stress. Historically this idea was implemented by expanding the strain tensor to first orderin tangentialdisplacement but secondorder in normaldisplacement. I derive the FvK strain by solving the evolution equations to first order in Vµ and secondorder in a, continuing to use the notation of previous sections, with the initial velocity vector W(0) =a∂ +V(0)µ∂ zG′µνa ∂ (34) z µ ,µ ν − of Eq. (11). I am using the superscript(0) to indicate the initial value, which is also the zeroth approximationfor an iterative solution. Other initial values are g = G(0) = δ and h = H(0) = 0. I use Picard’s method to generate the µν µν µν µν µν solution to the differential system Eqs. (19), (25, (27), and (28) iteratively as a power series in t, taking M to be the Euclidean plane with the usual Cartesian coordinates. In this case there is no distinction between indices up and indices down,andcovariantderivativesareordinarypartialderivatives. Iteratingonce, and ignoring quadratic terms except in a gives G(1) = δ +t(V(0)+V(0)) (35) κλ κλ κ,λ λ,κ H(1) = ta (36) κλ − ,λκ V(1) = taa (37) µ ,µ Iterating a second time, still ignoring quadratic terms except in a, gives G(2) =δ +t(V(0)+V(0)) 2ta +t2a a (38) κλ κλ κ,λ λ,κ − ,κλ ,κ ,λ Finally, evaluating at t=1, gives the FvK metric strain 1 1 (G(2) δ )= (V(0)+V(0)+a a ) a (39) 2 κλ − κλ 2 κ,λ λ,κ ,κ ,λ − ,κλ ThisisthecomputationalstartingpointforFvKtheory. Therestofthattheory follows from minimizing the elastic energy, expressed as a quadratic functional of this strain and the first order bending strain H(1), to find the equilibrium µν shape. The approachofthis paper isto developthe nonlinearstrainasthe solution toadifferentialsystem. FromthatpointofviewthederivationofEq.(39)isnot 7 very natural, since to obtain it one must artificially impose the condition that the trajectories of the particles are straight lines, a condition that introduces, viaEq.(37),asecondordercorrectionintothestrainthatisnecessarytoobtain Eq. (39). Although one can certainly parameterize the possible final shapes of M by displacementofparticles alongstraightlines, it is a differentthing to say that particles actually move along straight lines. FvK theory does not claim this, and in that sense it is not a dynamical theory. A dynamical theory would determine the evolution of the velocity vector (a,Vµ) by some local physical law, replacing Eq. (25) in the differential system. It would be a simpler theory, both conceptually and computationally, in that solving it would only require integrating a differential system forward in time. I will do the obvious thing and choose W to reduce the elastic energy at each step, seeking the minimum. A typical phenomenological elastic energy functional is E =E +E +E (40) d s c where 2 Λ √G E = 1 √gdx2dx3 (41) d 2 ZM √g − ! E = µ SκλS √gdx2dx3 (42) s κλ ZM E = κ G′µνH gµνh 2√gdx2dx3 (43) c µν µν 2 − ZM(cid:16) (cid:17) andwhereΛ,µ,κarethe2Dcompressionmodulus,shearmodulus,andbending modulus respectively. The area element involves √g, not √G, because the energy due to metric strain is better understood to be per unit mass, not per unitarea,andthe massisconvectedwiththematerialcoordinates. Thesystem will move,if possible, to lower its energy, so one must compute the variationof E with respect to a small normal displacement δa and tangential displacement δVµ δE δE δE = δVµ+ δa √gdx2dx3 (44) δVµ δa ZM(cid:20) (cid:21) The work done on M in deforming it represents energy given up by some other partof the system, so this workshould be added with a minus signto the total change in energy. Work done by pressure P in a small normal deformation δa, for instance, is W =P δa√Gdx2dx3 (45) ZM where now one must use the physical area element √Gdx2dx3 on M. A small displacement in the direction opposite to this “gradient,” i.e. δE √G a = L +P (46) a −δa √g! δE Vµ = L (47) V −δVµ (cid:18) (cid:19) 8 willlowerthe energyandmovethe systemtowarda localminimum. The linear operators L and L include a projection onto the space of admissible veloc- a V ity vector fields. They must define positive semi-definite quadratic forms with respect to the inner product givenby integrationoverM. Apart from these re- quirements,they will varywith the application. This is just the familiarnotion of conjugate gradient. One could also think of L and L together as defin- a V ing a generalized mobility tensor, because it transforms generalized force into velocity. If one only wants to know the final state, one could try to choose L a and L so as to reach equilibrium in the most efficient way. In any case, the V dynamics of the system is not completely determined by the elastic energies, andadditionalphysicalconsiderationsmustbeaddedto completethetheoryin a specific application. Eqs. (46) and (47), together with the evolution equations of Section 3, are what I mean by Lagrangian crumpling equations. The original problem ad- dressed by FvK theory, the bulging of a square plate fixed on the boundary and subject to pressure,can be solvedstraightforwardlyin this way. Represent all geometric data by discretization on a square grid of points of the original square. Spectral methods (fast Fourier transform with anti-aliasing) make the computation efficient, and the gradient flow converges quickly to a solution. 5 Buckling of a sphere under pressure I consider an elastic spherical shell subject to pressure P, described by the phenomenological energies of Eqs. (41), (42), (43), and (45). For small enough pressurethesphereisuniformlycompressed,butaspressureincreasesitbuckles. I will describe the buckling by using expansions of strain to second order in displacement where necessary, not the FvK expansion, but the “dynamic” one ofthis paper,foundby solvingthe crumplingequationsiteratively. Itturns out that the expansion must include more terms than FvK. For a sphere of radius R, in terms of spherical polar coordinates (θ,φ), g = diag(R2,R2sin2θ) (48) µν h = g /R (49) µν µν k = g /R2 (50) µν µν TakingR=1,andregardingallquantities nowasdimensionless,the perturbed geometric quantities in a general displacement (a,Vµ) are G = g +V +V +2ag (51) µν µν µ;ν ν;µ µν H = g +ag a +V +V (52) µν µν µν ,µ;ν µ;ν ν;µ − 1 √G = √g[1+(Vµ +2a)+ (Vµ Vν +VµVν ;µ 2 ;µ ;ν ;ν;µ + 4aVµ +2Vµa aaµ +2a2)] (53) ;µ ,µ− ;µ Theareaelement√Ghadtobe foundtosecondorderindisplacement. Tofirst 9 order in displacement the shear strain in the sphere is 1 S = (V +V Vλ g ) (54) µν 2 µ;ν ν;µ− ;λ µν Parameterize the displacement by coefficients (a ,b ,c ), such that ℓm ℓm ℓm a = a Y (55) ℓm ℓm ℓm X Vµ = gµν b Y +ǫµν c Y (56) ℓm ℓm,ν ℓm ℓm,ν ℓm ℓm X X where the Y are spherical harmonics and ǫ = ǫ =sinθ, ǫ =ǫ =0 is ℓm 32 23 22 33 − the antisymmetric tensor. Then for example the change in the mean curvature of the perturbed sphere is δH =GµνH gµνh = [ℓ(ℓ+1) 2]a Y (57) µν µν ℓm ℓm − − ℓm X so that the curvature energy is κ E = [ℓ(ℓ+1) 2]2 a 2 (58) c ℓm 2 − | | ℓm X Itvanishes for ℓ=1, asit must by Galileaninvariance,andit is independent of the tangential displacement Vµ. The other energy expressions are W = 4πPa Y P ℓ(ℓ+1)a b +2P a 2 00 00 ℓm ℓm ℓm − | | ℓm ℓm X X 1 + Pa Y [ 2ℓ(ℓ+1)a b +ℓ(ℓ+1)a 2+2a 2] (59) 00 00 ℓm ℓm ℓm ℓm 2 − | | | | ℓm X Λ E = [ ℓ(ℓ+1)b +2a ]2 d ℓm ℓm 2 − ℓm X + Λa Y [ 2ℓ(ℓ+1)a b +ℓ(ℓ+1)a 2+2a 2] (60) 00 00 ℓm ℓm ℓm ℓm − | | | | ℓm X µ E = ℓ(ℓ+1)[ℓ(ℓ+1) 2](b 2+ c 2) (61) s ℓm ℓm 2 − | | | | ℓm X These expansions have been carried out to second order in all coefficients, but they anticipate that a is the same order as a 2 for ℓ > 1, so that some 00 ℓm | | terms quadratic in a appear to be third order. I also anticipate that the first 00 response to pressure is a uniform compression P a (62) 00 | |∼ Λ sothatconsistencyrequiresP <<Λ. Nowseektheminimumofthetotalenergy E =W +E +E +E (63) tot d s c 10