Magnetoelasticity of highly deformable thin films: Theory and simulation

International Journal of Non-Linear Mechanics 47 (2012) 185–196

Contents lists available at ScienceDirect

International Journal of Non-Linear Mechanics journal homepage: www.elsevier.com/locate/nlm

Magnetoelasticity of highly deformable thin films: Theory and simulation$ M. Barham a,b, D.J. Steigmann a,, D. White b a b

Mechanical Engineering, University of California, Berkeley, CA 94720, USA Lawrence Livermore National Laboratory, Livermore, CA 94550, USA

a r t i c l e i n f o

abstract

Available online 19 May 2011

A non-linear two-dimensional theory is developed for thin magnetoelastic films capable of large deformations. This is derived directly from the three-dimensional theory. Significant simplifications emerge in the descent from three dimensions to two, permitting the self-field generated by the body to be computed a posteriori. The model is specialized to isotropic elastomers and numerical solutions are obtained to equilibrium boundary-value problems in which the membrane is subjected to lateral pressure and an applied magnetic field. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Magnetoelasticity Thin films Dynamic relaxation

1. Introduction There is considerable current interest among mechanicians in non-linear magnetoelasticity [1–5]. This is due to the development of highly deformable magnetizable materials synthesized from elastomers infused with micro- or nano-scopic ferrous particles [6]. Such materials are capable of large deformations induced by magnetic fields. This property may be used to facilitate controlled pumping of fluid, for example, via remote actuation. In the present work we continue our development [7] of a membrane theory for thin films composed of such materials. This is used to simulate membrane response to an applied magnetic field and to a pressure transmitted to the material by a confined gas. Section 2 contains a summary of three-dimensional magnetoelasticity and its specialization to isotropic elastomers. A corresponding membrane model is derived in Section 3 directly from the equations of the three-dimensional theory. It incorporates a constraint requiring the magnetization to remain tangential to the film as it deforms. This is motivated by the fact that such states are energetically optimal in thin films [8,9]. Likewise, we impose the constraint of bulk incompressibility, and thus exclude dilational modes of deformation that are energetically unfavorable in typical elastomers. However, unlike incompressibility, the constraint on magnetization is not of the kind that requires a reactive Lagrange multiplier in the relevant constitutive equation. Rather, it is a restriction involving the deformation, allowing local membrane geometry to adjust in response to an applied field. Constraints on the deformation of the Kirchhoff–

Love type are typically imposed at the outset in theories of thin magnetoelastic plates [10]. However, in general such constraints impede the attainment of minima of the overall energy because, by confining attention to states of magnetization that are optimal at any deformation, we effectively eliminate magnetization as an independent variable. The bias induced by an applied field then yields deformations that violate constraints of the Kirchhoff–Love type. Here, this is addressed via a director field which emerges naturally from the underlying three-dimensional theory in the manner described in [11] for the purely mechanical problem, without restricting the nature of the deformation in thin bodies. In Section 4 we use a finite-difference method to discretize the model spatially and discuss the solution of the resulting equations by the method of dynamic relaxation, in which equilibria are obtained as long-time limits of solutions to an artificial dynamical system with viscosity. The method is applied, in Section 5, to determine the deformation, magnetization and magnetic field generated by a thin film in response to an applied dipole field and pressure load. Notation follows standard usage in non-linear continuum mechanics [12]. Thus, boldface is used to denote vectors and tensors, bold subscripts are used to denote derivatives with respect to the indicated tensor or vector variables, the superscript t is used to denote transposition, and the superscripts 1 and t to denote the inverse and inverted transpose. The symbol  refers to the tensor product of vectors. A dot between variables in bold face is used to denote the standard Euclidean inner product, and j  j refers to the induced norm.

2. Three-dimensional magnetoelasticity $

For Ray Ogden, Prager Medalist.  Corresponding author. E-mail address: [email protected] (D.J. Steigmann). 0020-7462/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijnonlinmec.2011.05.004

The background material on continuum electromagnetism underlying this work may be found in [9,13–17]. We apply this

186

M. Barham et al. / International Journal of Non-Linear Mechanics 47 (2012) 185–196

to the description of incompressible magnetic elastomers undergoing large deformations. A summary of the relevant equations is given followed by discussions of restrictions associated with stable equilibria and the specialization of the theory to isotropic materials.

R,

ð1Þ

where T ¼ rðxF ÞF þ m0 ðh  h12jhj2 IÞ þ m0 h  mqI

ð2Þ

is the magnetoelastic Cauchy stress; xðF,mÞ is the free energy per unit mass; r is the mass density (mass per unit current volume); h is the magnetic field; m is the magnetization per unit current volume; F ¼ Dv is the gradient of the deformation function y ¼ vðx,tÞ in which x is the position of a material point in a fixed reference configuration k and D is the gradient with respect to x; superposed dots are used to denote material derivatives; R is the configuration occupied by the body at time t; and m0 ð 4 0Þ is the free-space permeability. Here I is the unit tensor, div is the spatial divergence based on y and q is a Lagrange-multiplier field associated with the incompressibility constraint. Maxwell’s equations may be used [9] to show that divfh  ðh þ mÞ12jhj2 Ig ¼ ðgrad hÞm,

ð3Þ

where grad is the gradient with respect to y, and thus furnish an equivalent equation of motion: € div½rðxF ÞFt qI þ m0 ðgrad hÞm ¼ ry,

ð4Þ

which proves, for reasons discussed below, to be more convenient for our purposes. Here we have suppressed time derivatives in Maxwell’s equations. This is justified in the absence of electric fields if, as assumed here, there are no free charges or currents and the body is not electrically polarized (see [17]). If ta is the applied (i.e, non-electromagnetic) traction acting on a part @Rt of the boundary @R, then [9]

rðxF ÞF

ð12Þ

where @Rt ¼ wð@kt Þ, having used Nanson’s formula: ð13Þ

pa ¼ ata

ð14Þ

is the applied traction measured per unit area of @kt : The magnetic field is the sum [9] h ¼ ha þ hs

t

t

@kt ,

on

where a ¼ jFn Nj is the local areal dilation of @kt : Here,

The local equation of motion in the absence of electric fields or applied (as distinct from electromagnetic) body forces is [9,17] in

PN ¼ pa þ 12m0 ðm  nÞ2 Fn N

an ¼ Fn N,

2.1. Basic equations

div T ¼ ry€

in which J¼ 1 has been imposed. Further, we find the referential form of the boundary condition (5) to be

nqn ¼ ta þ 12 0 ðm

m

2

 nÞ n

on

@Rt :

ð5Þ

Typical boundary-value problems, including those considered here, entail the assignment of y on the complement @R\@Rt . This system is augmented by the incompressibility constraint

rðvðx,tÞ,tÞ ¼ rk ðxÞ; equivalently, J ¼ 1, where J ¼ detF:

ð6Þ

Our further considerations require equations involving a referential divergence operator. For (4), this is easily achieved via the Piola transformation: t

n

n

P ¼ ½rðxF ÞF qIF ¼ WF qF ,

ð7Þ

ð8Þ

ð9Þ

ð10Þ

where Div is the referential divergence based on x; therefore, (4) is equivalent to € Div P þ m0 ðgrad hÞm ¼ rk y,

ð16Þ

in all of three-space, denoted by E, where curl is the spatial curl operation based on y, whereas curl hs ¼ 0

ð17Þ

in E\@R: The self-field and the magnetization are subject to the jump condition [9]: ½hs  ¼ ðn  mÞn

on

@R,

ð18Þ

where ½ is the difference between the exterior and interior limits of the enclosed variable on @R, and to Maxwell’s equation [9]: div hs ¼ div m

in

R ¼ 0,

and

in

E\R:

ð19Þ

The field ha is assumed to be assigned as a function that is continuously differentiable everywhere in E except at a finite number of singularities in E\R: In the examples discussed in Section 5 we study the response of the material to an applied field generated by a dipole source with the poles aligned along a fixed unit vector k. Accordingly [7], ha ðyÞ ¼ ‘D3 ½3ða  kÞak,

ð20Þ

where the (signed) constant D is the dipole strength, ‘ is the distance from the source to the point with position y A E, and ‘a ¼ yyd ,

ð21Þ

in which jaj ¼ 1, is position measured from the source, located at yd . This has an isolated singularity at the source. The associated gradient, required in (11), is [7] grad ha ¼ 3D‘4 f½ða  kÞI þ a  kP½3ða  kÞak  ag, where

P ¼ Ia  a,

ð22Þ

and is symmetric in accordance with (16). From (17)–(19) we have

ð11Þ

ð23Þ

where the scalar-field js satisfies on

@R

ð24Þ

and the magnetostatic equation: divðgrad js Þ ¼ div m

is the cofactor of the deformation gradient. Thus, J div½rðxF ÞFt qI ¼ Div P,

curl ha ¼ 0

½grad js  ¼ ðn  mÞn

is the referential strain-energy density, and Fn ¼ JFt

of an applied field ha , generated by remote sources, and the selffield hs generated by the magnetized body. In the present circumstances both satisfy the relevant Maxwell equation without time derivatives; thus

hs ¼ grad js ,

where WðF,mÞ ¼ rk x

ð15Þ

in

R ¼ 0,

and

in

E\R:

ð25Þ

At any given time the unique solution satisfying js  jyj1 as jyj-1 is [14,16] Z Z mðy0 Þ  nðy0 Þ div mðy0 Þ 0 daðy dvðy0 Þ for y= 4pjs ðyÞ ¼ Þ 2@R: 0j 0 jyy @R R jyy j ð26Þ

M. Barham et al. / International Journal of Non-Linear Mechanics 47 (2012) 185–196

187

The magnetization and magnetic field are related constitutively by [9]

where L is a suitable load potential. We then have the conservation law

Wm ¼ m0 h ¼ m0 ðha grad js Þ:

d 0 E dt

ð27Þ

Thus, if the constitutive function WðF,mÞ is known, Eqs. (11)–(14) and (23), (26), (27) yield a coupled integro-differential system to be solved for the deformation and magnetization. This presents considerable analytical and numerical challenges [18]. In [7] these were avoided by considering the limit of a weakly magnetized body in the presence of a strong applied field. In this limit the self-field may be generated from (26) a posteriori, and plays only a passive role in the analysis. Alternatively, a direct simulation of the field may be based on a discretization of Maxwell’s equations in the space surrounding the body [4,19]. In Section 3 we use a result derived in [8] for thin films to show that the tractability of the formulation adopted in [7] is retained when the magnetization and applied fields are comparable in magnitude. This yields a conventional differential-algebraic system to be solved on a reference surface associated with the thin film. To facilitate subsequent analysis, we use a pull-back M of m defined by Z Z m  n da ¼ M  N dA, ð28Þ s

S

in which S  k is an arbitrary orientable surface and s ¼ wðS,tÞ  R is its image in the current configuration. Nanson’s formula then furnishes M ¼ JF1 m:

ð29Þ

In particular, this yields the convenient connections

am  n ¼ M  N and J div m ¼ Div M,

ð30Þ

which enable us to use, in place (12) and (26), respectively, the equivalent expressions PN ¼ pa þ 12m0 a2 ðM  NÞ2 Fn N

on

@kt

ð31Þ

and 4pjs ðy,tÞ ¼

Z

Z Mðx,tÞ  NðxÞ Div Mðx,tÞ dAðxÞ dVðxÞ, @k jyvðx,tÞj k jyvðx,tÞj

for

x= 2@k, ð32Þ

in which the role of time has been made explicit and incompressibility has been imposed.

2.2. Stability and strong ellipticity A magneto-mechanical energy balance may be derived from (11), (12), (18), (19). Thus [9,20],   Z Z Z d K þ rx dvþ Mm0 ha  m dv ¼ ta  y_ da, ð33Þ dt R R @Rt where Z 1 _ 2 dv K¼ rjyj 2 R is the conventional kinetic energy and [3,9,15,16] Z 1 M ¼  m0 hs  m dv 2 R

ð34Þ

ð35Þ

is the magnetostatic energy of the self-field. In this work we consider conservative applied tractions for which Z d L, ð36Þ ta  y_ da ¼ dt @Rt

¼ 0,

where

E0 ¼ K þE

is the total magneto-mechanical energy in which Z Z E ¼ rx dvþ Mm0 ha  m dvL R

ð37Þ

ð38Þ

R

is the magnetoelastic potential energy. We remark that our energy balance excludes certain terms that are present in the balance discussed in [20]. These vanish collectively when the applied field is assigned as a stationary function of y, as assumed here; that is, as a function which is independent of t in the spatial description [9]. Further, the results of [9] may be used to show that the static specialization of (11), in which inertia is suppressed, furnishes an Euler–Lagrange equation for E. In this work we consider pressure acting on a part @Rt of the boundary formed by the union of two surfaces, @Rtþ and @R t , having no points in common. Uniformly distributed pressures, P þ and P , respectively, are acting on these surfaces. Let S be a fixed orientable surface such that @S ¼ C, the curve bounding @R t : We choose S such that its closure, and that of @R t , intersect only in C,  so that S [ @R t encloses a well-defined volume V  E: In the þ  applications of interest here, @Rt and @Rt , respectively are the ‘upper’ and ‘lower’ lateral surfaces of a thin sheet which, together with S, contains a compressible gas that transmits a pressure P to the lower surface. In Section 5 we identify S with the reference plane for the sheet. The upper surface is subjected to a fixed pressure P þ supplied by a large reservoir. This loading is conservative, and the associated potential, modulo an unimportant constant, is [9] Z V L¼ P ðvÞ dvP þ ðV þV  Þ, ð39Þ where P  ðV  Þ is the pressure–volume relation for the compressible gas and V is the volume of the body in configuration R. In the present context, the incompressibility of the magnetoelastic material allows us to suppress V on the right-hand side. Further, Z 1 V ¼  y  Fn N dA, ð40Þ 3 @kt  where @k t is the pre-image of @Rt in the reference configuration with exterior unit normal N [9]. In a full thermodynamic treatment accounting for dissipative effects, the energy balance (37) is replaced by the imbalance 0 dE =dt r 0 [20], so that if a state with vanishing initial velocity tends asymptotically to an equilibrium state, then the latter minimizes the potential energy E [5,20]; i.e., it furnishes a value of the potential energy not exceeding that supplied by the initial state. Because K is a positive-definite function of the velocity, it follows that E0 delivers a Lyapunov function for the dynamical system provided that the potential energy is strictly minimized at the equilibrium state. The considered equilibrium state is then stable. Without further qualification, this claim applies rigorously only to finite-dimensional systems [21]. Thus, we apply it only to the system that has been discretized for the purpose of numerical analysis. This is the basis of a dynamic relaxation method for computing equilibria (Section 5). In particular, then, an asymptotically stable equilibrium state minimizes the potential energy. In the purely mechanical setting, it is well known that a minimizing deformation necessarily satisfies the (local) strong-ellipticity inequality pointwise (see, for example, [22]). In the present setting this is replaced by the magnetoelastic strong-ellipticity inequalities [5]:

a  AðbÞa 4 0

and

c  ðWmm Þc 4 0,

ð41Þ

188

M. Barham et al. / International Journal of Non-Linear Mechanics 47 (2012) 185–196

where AðbÞ is the acoustic tensor defined by a  AðbÞa ¼ a  b  fWFF WFm ðWmm Þ

1

WmF g½a  b:

ð42Þ

These inequalities apply for all non-zero vectors a,b,c, with a and b subject to the restriction:

2

ð43Þ

associated with incompressibility. The second inequality implies that Wmm is invertible, as required by the first inequality. In terms of Cartesian components, inequalities ð41Þ1,2 are and

ð@2 W=@mi @mj Þci cj 4 0,

ð44Þ

where 2 Aij ðbÞ ¼ f@2 W=@FiA @FjB ð@2 W=@FiA @mk ÞðWmm Þ1 kl ð@ W=@ml @FjB ÞgbA bB :

ð45Þ They furnish pointwise restrictions on energy-minimizing states of deformation and magnetization jointly, which in turn play a central role in reducing the three-dimensional theory to a twodimensional membrane model (Section 3). 2.3. Reduced constitutive equations and isotropic materials Balance of moment of momentum requires that the Cauchy stress tensor be symmetric in the absence of electric fields [14,17]. Using (2), the symmetry requirement may be reduced to the statement: ðWF ÞFt þWm  m

is symmetric,

2

U ¼ m2fðC10 þC11 J1 =M s ÞðI1 3Þ þ ðC20 þ C21 J1 =M s ÞðI2 3Þ 2

a  Fn b ¼ 0

Aij ðbÞai aj 4 0

To use this formalism we adopt a magnetoelastic extension of the classical Mooney–Rivlin strain-energy function proposed in [5] and defined by

2

2

n ½coshðJ1 =M s Þ1g, þ C01 J1 =M s þC02 J2 =M s þ C01

ð56Þ

where J1 ¼ I5 I1 I4 þ I2 I6

and

J2 ¼ I6 ,

ð57Þ

in which detF ¼ 1 has been imposed. Here m is the ground-state shear modulus, M s is the saturation value of magnetization per unit volume, and the Cij are dimensionless constants. Numerical values of m, Cij and m0 M s are given in Ref. [5, Table 2], where m0 is the free-space permeability. The symbols J1,2 are used in [5] to denote invariants based on magnetization per unit mass. These are recovered on dividing our invariants by r2k , and (56) takes this adjustment into account. Further, we have used ð49Þ2 and (53) to express the invariants adopted in [5], here based on magnetization per unit volume, in terms of the Ik. In [5] it is claimed that (56) satisfies ð41Þ1 without qualification. This comports with the fact that the standard Mooney–Rivlin model satisfies the purely mechanical strong-ellipticity condition at all deformations [24]. Inequality ð41Þ2 was also shown in [5] to be satisfied over a substantial range of strain.

3. Membrane approximation

ð46Þ

and so W may be written as a (different) function of C and M, if desired. We make use of this function in Section 3. We assume the material to be isotropic, with a center of symmetry, relative to the reference configuration k: Then [7]

We consider a body whose reference configuration k is a prismatic region generated by the parallel translation of a simply connected plane O with piecewise-smooth boundary curve @O: The closure of k is O  ½h=2,h=2, where O ¼ O [ @O and h is the (uniform) thickness. Let l be another length scale such as the diameter of a hole in O or a typical spanwise dimension. We assume that E6h=l 51, and, in the theoretical development, adopt l as the unit of length (l ¼1). We derive a two-dimensional membrane model by estimating the equations of the threedimensional theory to leading order in E: Further, we suppose the deformation to be C2 and the magnetization to be C1 in the interior of the body, so that the local equations of the foregoing theory apply almost everywhere. With minor loss of generality we assume the dipole in (20) to be orthogonal to the plane O, which is thus oriented by the unit vector k: The projection onto the plane is

W ðC,MÞ ¼ W ðRt CR,Rt MÞ

1 ¼ Ik  k

which is found, following [7], to be equivalent to the requirement: WðF,mÞ ¼ WðQF,QmÞ

for all rotations Q ,

ð47Þ

and this in turn is satisfied if and only if [7] WðF,mÞ ¼ W ðC,MÞ

ð48Þ

for some function W , where C ¼ Ft F

and

M ¼ Ft m:

ð49Þ

The latter is related to the pull-back M by (cf. (29)) JM ¼ CM,

ð50Þ

for all orthogonal R:

ð51Þ

For R ¼ I this yields W ðC,MÞ ¼ W ðC,MÞ, which is satisfied if ^ ðC,M  MÞ for some function W ^ and only if [9] W ðC,MÞ ¼ W subject to the restriction: ^ ðC,M  MÞ ¼ W ^ ðRt CR,Rt ðM  MÞRÞ W

for all orthogonal R:

ð52Þ

For incompressible materials, standard representation theory ^ ¼ UðI1 ,I2 ,I4 I6 Þ for some function U, where [23] implies that W I1 ¼ trC,

I2 ¼ 12½I12 trðC2 Þ,

I5 ¼ C  M  M, I6 ¼ M  M:

ð53Þ

and

Wm ¼ FW M ,

ð54Þ

with

ð59Þ

of the Piola transform (7). Let B be a linear coordinate in the direction of k, and suppose B ¼ 0 on O. Eq. (11) is then equivalent to € DivJ ðP1Þ þ P0 kþ m0 ðgrad hÞm ¼ rk y,

W M ¼ 2ðU4 C þ U5 C2 þ U6 IÞM,

where Uk ¼ @U=@Ik :

ð60Þ

where ðÞ ¼ @ðÞ=@B and DivJ is the (referential) two-dimensional divergence with respect to position u on O, where ð61Þ

This holds at all points in the interior of the body and therefore at

B ¼ 0 in particular. Thus, DivJ ðP0 1Þ þ P00 k þ m0 ðgrad hÞ0 m0 ¼ rk0 y€ 0 ,

Sym W C ¼ ðU1 I1 U2 ÞI þ U2 Cþ U4 M  M þ U5 ½CðM  MÞ þðM  MÞC

and

P ¼ P1 þ Pk  k

x ¼ u þ Bk:

Proceeding as in [9] we then obtain WF ¼ 2FðSym W C Þ þ m  W M

and generates the orthogonal decomposition

0

I4 ¼ C  M  M,

2

ð58Þ

ð55Þ

ð62Þ

where the subscript 0 identifies the values of functions at B ¼ 0; i.e., on the plane O: For example, P0 ¼ WF ðF0 ,m0 Þq0 Fn0

and

m0 h0 ¼ Wm ðF0 ,m0 Þ

ð63Þ

M. Barham et al. / International Journal of Non-Linear Mechanics 47 (2012) 185–196

which yields

in which [11] F0 ¼ f þ d  k,

ð64Þ

where f ¼ rr,

ð65Þ

rðu,tÞð ¼ y0 Þ and dðu,tÞ are the restrictions to O of v and v0 , respectively, and r is the two-dimensional gradient on O; i.e., the gradient with respect to u: We note that rðu,tÞ maps O to the deformed membrane surface o ¼ wðO,tÞ: Accordingly, f maps O0 , the translation space associated with the plane O, to To , the tangent plane to o at the material point u A O. To accommodate the constraint of bulk incompressibility we impose 1 ¼ detF0 ¼ Fn0 k  F0 k ¼ an  d,

ð66Þ

where (64) and Nanson’s formula (13) have been used in the final equality. Here, a is the local areal dilation of O and n is the orientation of the surface onto which O is deformed; i.e., the unit normal to To : The general solution is d ¼ a1 n þfe,

ð67Þ

where eA O0 is arbitrary. Further, Eqs. (13) and (64) yield

an ¼ fi1  fi2 ,

ð68Þ

where ia A O0 are subject only to the requirement that fi1 ,i2 ,kg be a positively oriented orthonormal set. Thus F0 is determined by f and e, regarded as independent variables. The associated Cauchy–Green deformation tensor, C0 ¼ Ft0 F0 , is C0 ¼ c þce  kþ k  ce þ ða2 þe  ceÞk  k, where

t

c¼f f

ð69Þ

and a is obtained by evaluating the norm of (68), yielding pffiffiffiffiffiffiffiffiffiffi a ¼ detc:

ð70Þ

Div M ¼ DivJ ð1 MÞ þM 0 :

3.1. The leading-order model The foregoing equations, holding on O, are exact consequences of the three-dimensional theory. Approximations arise in using them to represent material response in O  ½E=2, E=2: Let P 7 be the interior limits of P as B- 7 E=2, where the exterior unit normals are N 7 ¼ 7k: Their Taylor expansions yield þ





P N þP N ¼ E

P00 k þ oð

þ

þ





EÞ and P N P N ¼ 2P0 k þ oðEÞ: ð71Þ

On the left-hand sides we use (12) together with the estimates ðFn NÞ 7 ¼ 7ðFn Þ 7 k ¼ 7Fn0 k þ ðE=2ÞðFn Þ00 kþ oðEÞ

ð72Þ

ð78Þ

It follows from (71) and (75) that (62) yields a well-defined differential equation in the limit of small E only if P00 k remains bounded. Further, (63) implies that the deformation gradient and magnetization are bounded on O only if P0 is bounded. From (75) and (76) it is therefore necessary that paþ þ p a ¼ Ep þ oðEÞ

and

paþ p a ¼ 2q þoð1Þ,

ð79Þ

where p and q are of order unity in magnitude. It follows that, to leading order in E, n n 0 0 0 1 P00 k ¼ p þ m0 a2 0 M0 f½M 0 ða0 =a0 ÞM0 F0 k þ 2M0 ðF Þ0 kg

and

2 n P0 k ¼ q þ 12m0 a2 0 M0 F0 k:

ð80Þ

3.2. Estimate of the self-field Before proceeding we obtain an estimate of the leading-order self-field potential (32). An elementary calculation based on (77) and (78) gives Z  Z ½DivJ ð1 M0 Þ þ M 00  1 M0  m dS 4pjs ðy,tÞ ¼ E dA jyrj @O jyrj O Z Z Mþ M dA dAþ oðEÞ, ð81Þ þ þ  @k þ jyv j @k jyv j where the superscripts 7 identify the values of functions at the upper and lower lateral surfaces @k 7 ¼ O  f 7 E=2g and m A O0 is the unit normal exterior to O: This is valid provided that y arðu,tÞ for any u A O : To estimate the associated integrals we compute jvj0 ¼ jvj1 v  v0 , where v ¼ yvðx,tÞ and the derivative is with respect to B at fixed y. Accordingly, v0 ¼ Fk, and (64) gives jyvj00 ¼ 

þ

189

ðyrÞ  d: jyrj

For y ar this yields   1 1 E ðyrÞ 17  d þoðEÞ, 7 ¼ 2 jyrj 2 jyrj jyvj

ð82Þ

ð83Þ

which, when combined with M 7 ¼ M0 7 ðE=2ÞM 00 þoðEÞ,

ð84Þ

results in Z  1 M0  m M0 dSþ ðyrÞ  d 2 @O jyrj O jyrj  DivJ ð1 M0 Þ dA þoðEÞ=E,  jyrj

4pjs ðy,tÞ=E ¼

Z

ð85Þ

and

a 7 ¼ a0 7ðE=2Þa00 þ oðEÞ,

ð73Þ

where n n 0 a00 ¼ a1 0 F0 k  ðF Þ0 k,

ð74Þ n

which follows on differentiation of a ¼ jF kj: After some algebra we obtain 2 0 0 P þ N þ þP N ¼ paþ þp a þ Em0 a0 M0 f½M 0 ða0 =a0 ÞM0 

Fn0 k þ 12M0 ðFn Þ00 kg þoðEÞ

ð75Þ

and 2 2 n P þ N þ P N ¼ paþ p a þ m0 a0 M0 F0 k þ OðEÞ,

ð76Þ

pa7

are the applied tractions at the lateral surfaces and where M ¼ M  k: The role of the latter suggests the decomposition M ¼ 1 Mþ Mk,

ð77Þ

provided that (83) is uniformly valid over the domain. This limitation effectively restricts the use of (85) to points y whose distances from the membrane are of order unity compared to E; that is, to points in space whose minimum distances from the deforming membrane surface are large compared to membrane thickness. Accordingly, it may not be used to describe the selffield in the interior of the material. To characterize the magnetic state inside the film, we estimate (32) at an interior point x A k. For points x near x, the presumed differentiability of the deformation implies that jyvðx,tÞj ¼ OðjnjÞ, where n ¼ xx and y ¼ vðx,tÞ: The self-field is obtained by computing the gradient of js with respect to y and evaluating the result at y; thus, for x= 2@k, Z Z ðM  NÞu ðDiv MÞu 4phs ðyÞ ¼ dA dV, 2 jy jy v ðxÞj vðxÞj2 @k k ð86Þ where u ¼ ðyvðxÞÞ=jyvðxÞj,

190

M. Barham et al. / International Journal of Non-Linear Mechanics 47 (2012) 185–196

in which t has been suppressed. The singularity is of order jnj2 , which is integrable in k. Therefore the volume integral makes a contribution of order E. The boundary integral includes a contribution from the surface @O  ðE=2, E=2Þ, on which jyvðxÞj is strictly bounded away from zero for any E: Accordingly, this too contributes at order E, leaving Z ðM  NÞu dAþ OðEÞ, ð87Þ 4phs ðyÞ ¼ þ  jy vðxÞj2 @k [@k in which M  N ¼ 7M 7 on @k 7 , respectively. Thus, Z  Z   ðM  NÞu jM 7 j   dA r dA   @k þ [@k jyvðxÞj2  þ  jy vðxÞj2 @k [@k Z 1 jM 7 j dA: r max 2 @k þ [@k @k þ [@k jyvðxÞj

P00 k ¼ p

ð88Þ

The integrand in the final integral is dominated by its asymptotic behavior near x; i.e., by jnj2 . For small thickness, the integral may then be shown to be OðjlnEjÞ in magnitude. In view of (84), the upper bound remains finite in the limit only if maxO jM0 j ¼ 0, in which case it is of order jElnEj: This guarantees that jhs ðyÞj is finite and vanishes with E. In particular, then, hs

vanishes on

O, at leading order:

ð89Þ

The alternative (M0 a 0Þ yields an upper bound of order jlnEj, which allows the self-field to grow without bound as thickness tends to zero. In this case the magnetostatic energy, and therefore the potential energy, may become unbounded. However, this alternative does not require the self-field to become unbounded, and so our analysis, while suggestive, is not conclusive. In other words, we have only shown that the constraint: M0 ¼ 0

on

O

ð90Þ

is sufficient for (89) and for boundedness of the magnetostatic energy. Using (64) and (77), we then derive m0 A To

on

o, at every u A O:

ð91Þ

To explore this issue further, consider the part of the potential energy involving magnetization. This is (cf. (35), (36) and (38)) Z 1 ð92Þ Emag ¼ ½W ðh þha Þ  m dV 2 k in which y and F are fixed, and reduces to Z Emag ¼ ðWha  mÞ dV

function of y. The claim then follows from (16), which is equivalent to the symmetry of grad ha : Strictly, these results are known to be necessary only for optimal (energy-minimizing) states of magnetization and so may not apply to dynamical states. However, in this work we use dynamics solely to facilitate the computation of equilibria. We do not model actual dynamic interactions. Accordingly, we restrict attention to states of magnetization that are energetically optimal at any deformation, equilibrated or otherwise. Eq. (90) affords the important simplifications: and

P0 k ¼ q:

ð94Þ

For points remote from the deforming film (85) is applicable and simplifies, by virtue of (90), to  Z 1 dA þ oðEÞ=E, ð95Þ 4pjs ðy,tÞ=E ¼ 1 M0  r jyrj O where r is the two-dimensional gradient with respect to u A O and Green’s theorem has been used to combine terms. Proceeding as in the calculation leading to (85), we put vðu,tÞ ¼ yrðu,tÞ and use (65) to derive dðjvj1 Þ ¼ jvj3 v  dv,

where

dv ¼ drðuÞ ¼ fðduÞ,

yielding  1 t ¼ jyrj3 f ðyrÞ r jyrj

ð96Þ

ð97Þ

and 4pjs ðy,tÞ=E ¼

Z

fM0  ðyrÞ dA þoðEÞ=E: 3 jyrj O

ð98Þ

A straightforward computation based on (23) generates the scaled self-field in the surrounding space: Z 4phs ðy,tÞ=E ¼ GðfM0 Þ dAðuÞ þ oðEÞ=E, O

3 1 where G ¼ ðyrÞ  ðyrÞ I jyrj5 jyrj3

ð99Þ

in which rðu,tÞ is the membrane position field at time t. Thus, the leading-order model generates the dominant part of the self-field in the surrounding space (which is of order EÞ a posteriori. 3.3. Loading

ð93Þ

k

if the self-field is negligible; i.e., if (91) holds. Here the magnetization is obtained by solving (27) in which the self-field is suppressed, so that Emag is controlled entirely by the deformation. This effectively eliminates the magnetization as an independent variable. In the work of Gioia and James [8] on non-deforming films it is proved that minimizers of (92) furnish energies that converge to (93) as film thickness tends to zero. It was also proved that optimal states of magnetization necessarily satisfy (91) and that the residual self-field vanishes, in accordance with (89) (see also [9]). Further, in [8] it is shown that there is no residual magnetostatic equation to leading order; indeed, the solution (26) to (25) has already been used in the course of obtaining (91) and therefore plays no further role. These results imply that (89) and (91) characterize optimal states of magnetization in a sufficiently thin film, at any fixed deformation. In particular, the magnetostatic energy is negligible at leading order. The Euler equation for the deformation that emerges from this leading-order approximation is given by (11) [9], but with grad h replaced by grad ha . This follows from the fact that the variational derivative of ha , identified by a superposed dot, is purely _ if the applied field is a stationary convective; i.e., h_ a ¼ ðgrad ha Þy,

Turning now to the loading, suppose the lateral surfaces are subjected to pressures P 7 . The applied tractions are pa7 ¼ 8ðP 7 ÞðFn Þ 7 k,

ð100Þ

and we assume that P 7 ¼ Ep 7 þoðEÞ,

ð101Þ

where p 7 are of order unity. In this case q ¼ 0 and p ¼ aðDpÞn

ð102Þ 

þ

in (94), where Dp ¼ p p is the net lateral pressure across the membrane. Invoking (89) and the foregoing thin-film approximations, we find that (62) reduces to DivJ ðP0 1Þ þ aðDpÞn þ m0 ðgrad ha Þ0 m0 ¼ rk0 r€ ,

ð103Þ

to leading order, where ðgrad ha Þ0 is evaluated using (21) and (22) with y replaced by r, and P0 ¼ WF ðF0 ,m0 Þq0 Ft 0 :

ð104Þ

This is augmented by the algebraic constraints ð63Þ2 and ð94Þ2 . Using (100), (101) and y0 ¼ vðu,tÞ ¼ r, the leading-order

M. Barham et al. / International Journal of Non-Linear Mechanics 47 (2012) 185–196

constraints are found to be Wm ðF0 ,m0 Þ ¼ m0 ha ðrÞ

and

P0 k ¼ 0:

ð105Þ

Together with (64), (67), (77), (90) and (104), these furnish a system for the determination of rðu,tÞ, eðu,tÞ, M0 ðu,tÞð ¼ 1 M0 Þ and q0 ðu,tÞ: In practice we solve the equation obtained on multiplying (103) through by E: This yields the equation of motion for the membrane, which in turn furnishes the leading-order approximation to that for a thin sheet. Our preference for (4) over (1) is due the availability of an explicit formula for the gradient of the applied field (cf. (22)). From ð105Þ1 and (64) it is clear that the constraint (91) imposes a restriction, not only on the magnetization, but also on the deformation and director fields r and e, and thus on the geometry of the film in the presence of an applied field. In particular, this allows the orientation of the tangent plane to the membrane to adjust in response to the applied field. The literature on magnetoelasticity in thin structures is typically based on an a priori constraint of the Kirchoff–Love type (i.e., e ¼ 0Þ on the director field (e.g. [10]). However, in Section 5 we find that solutions deviate substantially from Kirchhoff–Love kinematics. Because we have confined attention to states of magnetization that are optimal at any deformation, and thereby eliminated magnetization as an independent variable, it follows by relaxation of constraints that restrictions of the Kirchhoff–Love type impede the attainment of minima of the overall potential energy. Indeed, the analysis of [7] indicates that the Kirchhoff– Love constraint is generally incompatible with (91). Therefore the present model is optimal relative to formulations in which such constraints are imposed at the outset. Kirchhoff–Love kinematics obtain if the effects of deformation and magnetization are uncoupled in the expression for the strain-energy function, as in weakly magnetized materials subjected to applied fields of sufficient intensity [7]. Standard mixed traction/position problems consist of the specification of r and the traction:

s ¼ P0 1m,

ð106Þ

on complementary parts of the boundary curve @O: Here s is the value on @O of the exact traction field acting on a part of the cylindrical generating surface of the body where tractions are assigned. In this work we assume position to be prescribed on @O  ½E=2, E=2 and thus assign r everywhere on @O:

191

We regard this as an equation for M in which r, f and e (hence FÞ are assigned. To investigate its solvability we compute another derivative, again at fixed F, obtaining ~ MM ÞM _ ¼ f t ðWmm Þm _ _ ¼ ½f t ðWmm ÞfM: ðW

ð110Þ

Therefore, _ ¼ fM _  ðWmm Þf M, _ _  ðW ~ MM ÞM M

ð111Þ

_ by virtue of ð41Þ . Accordwhich is positive for all non-zero M 2 ~ MM is positive-definite and W ~ ðF,Þ is strictly convex. ingly, W ^ which miniEq. (109) therefore possesses a unique solution M ~ at fixed F: This in turn determines the magnetization mizes W ^ which furnishes the unique solution to m ¼ f M, Wm ¼ m0 ha ðrÞ:

ð112Þ

Next, we fix f and define Gðe,mÞ ¼ Wðf þ dðf,eÞ  k,mÞ,

ð113Þ

where dðf,eÞ is the function defined by (67), in which a and n are determined by f via (68) and (70). Consider one-parameter families eðuÞ and mðuÞ: The former induces the one-parameter family FðuÞ of deformation gradients with derivative F_ ¼ f e_  k (cf. (64) and (67)). Accordingly, t _  Wm , G_ ¼ e_  f ðWF Þk þ m t

Ge ¼ f ðWF Þk

yielding

and

Gm ¼ Wm :

ð114Þ

t

Using (104) and the invertibility of F , we find the constraint ð105Þ2 to be equivalent to Ge ¼ 0

and

q ¼ d  ðWF Þk:

ð115Þ

To address the first of these equations, we keep f fixed and compute: t _ ðGe Þ ¼ f fWFF ½f e_  k þ ðWFm Þmgk:

ð116Þ

Here we regard mðuÞ as the magnetization induced by eðuÞ via ^ (112); that is, we use the unique solution M ¼ MðeÞ to (109), ^ associated with fixed r and f, to generate mðuÞ ¼ f MðeðuÞÞ: This satisfies (112) identically for all eðuÞ with u in some open interval. It follows that _ ¼ ðWmm Þ1 ðWmF Þðf e_  kÞ, m

ð117Þ

so that (116) reduces to t

_ ðGe Þ ¼ f fAðkÞgðf eÞ, 3.4. Solvability of the constraints We demonstrate the solvability of the constraints ð105Þ1,2 for M and e at a given deformation rðu,tÞ of O. To ease the notation, here and henceforth the subscript 0 is suppressed on the understanding that all fields discussed are the restrictions to O of threedimensional fields identified by the same symbol. We impose (90) at the outset, and thus find it more convenient to work with a formulation based on the use of F and M, rather than F and m, as independent variables. To this end we invoke bulk incompressibility and use (29) to define the function: ~ ðF,MÞ ¼ WðF,FMÞ W

for

M A O0 :

for all

_ A O0 , M

ð107Þ

ð108Þ

~ M A O0 . It follows that W ~ M ¼ 1Ft ðWm Þ, where 1Ft ¼ f t wherein W by virtue of (64); Eq. (108) then implies that ~ M ¼ m f t ha ðrÞ: W 0

where AðÞ is defined by (42). With these results in hand, we define a function GðeÞ by ^ ^ GðeÞ ¼ Gðe,f MðeÞÞ m0 ha ðrÞ  f MðeÞ,

ð119Þ

at the same r and f: Inserting eðuÞ and evaluating the derivative, we find from (112) and ð114Þ2 that

G_ ¼ Ge  e_

ð120Þ

for all u in some open interval. Then, € G€ ¼ ðGe Þ  e_ þ Ge  e:

ð121Þ 0

Consider a one-parameter family of magnetizations MðuÞ A O0 . Using a superposed dot to denote the derivative with respect to ~ MM _ ¼ Wm  FM _ at fixed F, and the parameter, we derive W therefore _ ¼0 ~ M Ft ðWm Þ  M ½W

ð118Þ

ð109Þ

The domain of GðÞ is the linear space O , a convex set. If e1,2 belong to this domain, then so do all points on the straight-line path eðuÞ ¼ ue2 þ ð1uÞe1 ;

u A ð0,1Þ,

ð122Þ

on which (121) reduces to _ G€ ¼ f e_  fAðkÞgðf eÞ;

e_ ¼ e2 e1 a 0:

ð123Þ

Setting a ¼ f e_ A To and b ¼ k, we find that (43) is satisfied because Fn k ¼ an is orthogonal to To (cf. (68)). Accordingly, the strong€ 40. ellipticity inequality ð41Þ1 is applicable and implies that G Integration of this result over ð0,uÞ and then again over (0,1)

192

M. Barham et al. / International Journal of Non-Linear Mechanics 47 (2012) 185–196

yields the conclusion that GðeÞ is strictly convex; i.e.,

Gðe2 ÞGðe1 Þ 4 Ge ðe1 Þ  ðe2 e1 Þ

ð124Þ

for all unequal pairs e1,2 , wherein Ge ¼ Ge by virtue of (120). Because strictly convex functions possess unique stationary points, we conclude that ð115Þ1 has a unique solution en : In particular, this solution satisfies GðeÞ 4 Gðen Þ for all e aen and therefore furnishes the unique minimizer of GðeÞ: With this solution in hand, the unique magnetization field associated with a given deformation rðu,tÞ, and attendant gradient f, is given by ^ n Þ: m ¼ f Mðe 3.5. Lyapunov functions We have shown that the constraints (109) and ð115Þ1 possess unique solutions e and M at fixed r and f: To obtain them, use may be made of the Newton–Raphson method, for example. The convexity conditions established in the foregoing ensure that the associated iterates converge to a unique solution. Alternatively, we may embed (109) and ð115Þ1 into the artificial dynamical problems: € þ cM _ þW ~ M ¼ m f t ha ðrÞ mM 0

and

me€ þ ce_ þGe ¼ 0,

ð125Þ

respectively, in which m and c are positive constants and the superposed dots in the two equations now identify derivatives with respect to time-like parameters t1,2 , respectively. Equilibria of this system are precisely the unique solutions to (109) and ð115Þ1 . Further, solutions to this system satisfy the energy balances:   d 1 d 1 _ 2 þW ~ ¼ cjMj _ 2 and _ 2 þ G ¼ cjej _ 2: mjMj mjej dt1 2 dt2 2 ð126Þ Standard theory for ordinary differential equations then ensures the existence of trajectories of ð125Þ1,2 for arbitrary initial data on which _ 2 þW ~ L1 ¼ 12 mjMj

and

_ 2 þG L2 ¼ 12mjej

ð127Þ

are strictly decreasing. Our results concerning the minimizing properties of equilibria then imply that L1,2 furnish Lyapunov functions for ð125Þ1,2 , respectively. All trajectories tend asymptotically to solutions of the constraints (109) and ð115Þ1 , and these are stable equilibria of the dynamical system [21]. The implementation of these results is discussed in Section 4. Finally, we use the energy balance (37) to construct a Lyapunov function for the motion rðu,tÞ: To this end we observe, using (34)–(36), (38) and (39) that K ¼ EKM þ oðEÞ,

L ¼ ELM þoðEÞ

and

E ¼ EEM þ oðEÞ,

ð128Þ

where KM ¼

1 2

Z

and

O

rk jr_ j2 dA, LM ¼

EM ¼

Z O

Z

V

p ðvÞ dvp þ V Z WðF,mÞ dAm0 ha ðrÞ  m dALM ,

ð129Þ

O

respectively, are the leading-order (membrane) approximations to the kinetic energy, pressure potential, and potential energy, in which Z 1 V¼ an  r dA ð130Þ 3 O is the volume of the compressible gas contained by the membrane. From the leading-order equation of motion (103), we obtain Z ð131Þ K_ M ¼ r_  ½DivJ ðP1Þ þ aðDpÞn þ m0 ðgrad ha Þm dA: O

Using [9] Z L_ M ¼ aðDpÞn  r_ dA

ð132Þ

O

this is reduced to Z Z Z P1m  r_ dS P1  f_ dA, K_ M ¼ L_ M þ m0 m  h_ a dA þ O

@O

ð133Þ

O

where we have used h_ a ¼ ðgrad ha Þr_ for stationary applied fields, together with the symmetry of grad ha . In this work we assume r to be fixed on @O and accordingly suppress the integral over @O. We now use ð105Þ1 and combine the result with ð129Þ3 to derive Z d _ ðKM þ EM Þ ¼ ðWF  FP1  f_ Þ dA: ð134Þ dt O Using (104) with the constraint of bulk incompressibility in the form Ft  F_ ¼ 0, together with (59) and (64), we find that WF  F_ ¼ P  F_ ¼ P1  f_ þPk  d_ and thereby reduce (134) to Z d ðKM þ EM Þ ¼ Pk  d_ dA, dt O

ð135Þ

ð136Þ

which vanishes by virtue of (105). Thus the leading-order model yields the conservation law ðd=dtÞðKM þ EM Þ ¼ 0, which is replaced, in the presence of dissipation, by the imbalance ðd=dtÞ ðKM þ EM Þ r 0: This observation suggests that a dissipative numerical scheme may be based on a discretization of the artificial dynamical equation: DivJ ðP1Þ þ aðDpÞn þ m0 ðgrad ha Þm ¼ rk r€ þ cr_ ,

ð137Þ

where c is a suitable constant. It is straightforward to show that if this equation is used in place of (103), then the leading-order energy balance is replaced by Z d ðKM þ EM Þ ¼ c jr_ j2 dA: ð138Þ dt O Our earlier observation that stable equilibria minimize E implies that EM is minimized, to leading order in thickness. Indeed, it is easily verified that (112) and the static specialization of (103) furnish the Euler–Lagrange equations for EM. Consequently, KM þ EM decays on trajectories of (137), provided that c 40, and achieves a strict minimum at a stable equilibrium. It therefore yields a Lyapunov function for (137), whose equilibria coincide with those of (103). This conclusion applies strictly only to a finite-dimensional projection of the problem associated with a spatial discretization of the equations on O: It also presumes that equilibria are minimizers of EM : Here, however, we have only imposed necessary conditions for a minimum of the energy. In particular, in the purely mechanical specialization of the theory it is known that these conditions are insufficient to preclude compressive stresses in equilibrium, which violate the Legendre–Hadamard necessary condition for minimizers of EM [25]. In such circumstances the existence of minimizers may be restored by replacing the membrane energy with a suitable relaxation [25–28] which excludes compressive stress a priori via the mechanism of fine-scale wrinkling. This is the subject of tension-field theory [29]. In this work we apply the theory to problems that do not exhibit wrinkling and therefore do not require the explicit relaxation. We emphasize the fact that (137) does not describe the actual dynamics of the membrane. Rather, it is used here solely to expedite the computation of equilibria by embedding the equilibrium problem into an artificial (finite-dimensional) dynamical system whose equilibria coincide with those of the physical problem. As such, it furnishes a convenient regularization of the equations. The strictly dissipative nature of this system is a feature shared by actual equations of motion that account for

M. Barham et al. / International Journal of Non-Linear Mechanics 47 (2012) 185–196

dissipation through constitutive equations rather than through modification of the equation of motion. However, (137) proves more convenient for the purpose of generating equilibria because it allows the discrete equations associated with the temporal evolution to be decoupled, affording a more efficient solution procedure. This is discussed in the next section.

4. Finite-difference scheme Equation (137) is discretized by using a finite-difference scheme derived from Green’s theorem. The application of this scheme to plane-strain problems in non-linear elasticity theory is described by Silling [30]. Its adaptation to membrane theory is developed in [26,27]. Here, we summarize the method and describe its extension to magnetoelasticity. The reference plane O is covered by a grid consisting of cells of the kind depicted in Fig. 1. Nodes are labeled using integer superscripts ði,jÞ: Thus, ui,j a are the referential coordinates of node ði,jÞ, where ua ¼ u  ia ; a ¼ 1,2: The four regions formed by a node, together with its nearest-neighbor nodes, are called zones. Zonecentered points, identified by open circles in Fig. 1, are labeled using half-integer superscripts. Green’s theorem may be stated in the form: Z Z sa, a da ¼ eab sa dub , ð139Þ D

where þ1 þ1 Ai,j ¼ 14½ðui1,j ui2þ 1,j Þðui,j ui,j1 Þðui1,j ui1þ 1,j Þðui,j ui,j1 Þ 2 1 1 1 2 2

ð141Þ is one-half the area of the quadrilateral. We also have need of gradients of various functions at zonecentered points. First, we apply (139) with sa ¼ ca sðu1 ,u2 Þ, where s is a smooth scalar-field and ca are arbitrary constants. This yields Z Z s, a da ¼ eab s dub : ð142Þ D

@D

We now identify D with the shaded region in the figure. The lefthand side is approximated by the product of the shaded area with the integrand, evaluated at the zone-centered point, and the four edge contributions to the right-hand side are approximated by replacing the integrand in each with the average of the nodal values at the endpoints. This gives [26] i þ 1=2,j þ 1=2

2Ai þ 1=2,j þ 1=2 ðs, a

i þ 1=2,j þ 1=2

2Ai,j ðsa, a Þi,j ¼ eab ½sa

i1=2,j þ 1=2

þ sa

i1=2,j1=2

þ sa

þ1 ðui1,j ui,j Þ b b

ðuibþ 1,j ui,j1 Þ, b

ð143Þ

where þ1 i,j þ 1 Ai þ 1=2,j þ 1=2 ¼ 12½ðui,j ui2þ 1,j Þðui1þ 1,j þ 1 ui,j ui1þ 1,j Þ 2 1 Þðu1

ðui2þ 1,j þ 1 ui,j 2 Þ:

ð144Þ

The magnetoelastic equilibrium equation is given by (137) in which the right-hand side is suppressed. To facilitate its discretization, we exploit the fact that the term an associated with the pressure load may be expressed as a divergence on O [31]. Thus, n ¼ nk ik , where i3 ¼ k,

ank ¼ 12eijk eab fia fjb ¼ Gkb, b ,

ð145Þ

and Gkb ¼ 12eijk eab fia rj

ð146Þ

in which eijk and eab , respectively, are the three- and twodimensional unit alternators (e123 ¼ e12 ¼ þ 1Þ: This result is useful in the present circumstances because the net lateral pressure on the membrane is uniformly distributed. Thus, the equilibrium equation is equivalent to the system

þ1 ðui,j uibþ 1,j Þ b

ðui,j1 ui1,j Þ b b

i þ 1=2,j1=2

þ sa

þ1 Þ ¼ eab ½ðsi þ 1,j þ 1 si,j Þðui,j uibþ 1,j Þ b

ðsi,j þ 1 si þ 1,j Þðuibþ 1,j þ 1 ui,j b Þ,

@D

where sa ðu1 ,u2 Þ is a smooth two-dimensional vector field, D is an arbitrary simply connected subregion of O and commas followed by subscripts identify partial derivatives with respect to the indicated coordinates. To approximate the divergence sa, a at node (i,j) we identify D with the quadrilateral contained within the dashed contour of Fig. 1. The left-hand side of (114) is estimated as the nodal value of the integrand multiplied by the area of D; the right-hand side as the zone-centered values of the integrand on each of the four edges of @D multiplied by the appropriate edge length. Thus [26],

193

ð140Þ

Tka, a ¼ Rk ,

where

Tka ¼ Pka þ ðDpÞGka

and

ðaÞ Rk ¼ m0 hk,i fia Ma :

ð147Þ hkðaÞ

¼ ik  ha are the Here Pka ¼ P ik  ia are the components of P1, components of the applied field, fka ¼ f  ik  ia ¼ rk, a are the components of the surface deformation gradient, rk ¼ ik  r are the Cartesian coordinates of a material point on the deformed surface and Ma ¼ M  ia are the magnetization components. Each of Eqs. (147) is of the form:

sa, a ¼ f ,

ð148Þ

where sa ¼ Tka and f ¼ Rk ; k ¼ 1,2,3: The sa , in turn, depend on the magnetization and on the gradients of s ¼ rk ðk ¼ 1,2,3Þ. To solve (148) at node (i,j) we integrate it over the region enclosed by the dashed quadrilateral of Fig. 1 with vertices at the nearestneighbor nodes, obtaining

Si,j ¼ F i,j ,

ð149Þ

where

Si,j ¼ 2Ai,j ðsa, a Þi,j

ð150Þ

and Fig. 1. Finite-difference mesh.

F i,j ¼ 2Ai,j f i,j :

ð151Þ

194

M. Barham et al. / International Journal of Non-Linear Mechanics 47 (2012) 185–196

In (150) the right-hand side is evaluated in terms of the zonecentered values of sa via (140). The latter are determined constitutively by corresponding zone-centered values of magnetization together with the gradients s, a which, in turn, are given via (143) by the values of s at the nodes located at the vertices of the shaded region of Fig. 1. The scheme is seen to require one degree less differentiability than that required by the local differential equations. Discussions of the associated truncation errors are given by Silling [30] and Hermann and Bertholf [32]. To solve Eq. (149) we introduce a regularization based on the artificial dynamical system (cf. (137)):

Si,j,n ¼ mi,j s€ i,j,n þ ci,j s_ i,j,n þ F i,j,n , i,j

i,j

ð152Þ i,j

i,j

where m ¼ 2A r is the nodal mass, c ¼ 2A c is the nodal damping coefficient, n is the time step, and superposed dots refer to derivatives with respect to (artificial) time. This is not the discrete form of the actual dynamical equations. Rather, it is an artificial system introduced solely to expedite the computation of equilibria. The basic method, known as dynamic relaxation [33], is a powerful tool for generating equilibria in a wide variety of nonlinear problems. It was developed for membrane theory in a purely mechanical setting in [26–28] and extended to coupled thermoelasticity in [31]. We observe that the matrix Gkb associated with lateral pressure is evaluated at zone-centered points (cf. (147) and (150)). However, this involves the deformation rk (cf. (146)), a nodal variable. The required evaluation is based on the average of the deformations at the four adjacent nodes. Similarly, (147) requires nodal values of fka Ma , which are obtained by averaging values at the four adjacent zone-centered points. In the case of volume-dependent pressure loading it is necessary to evaluate the volume enclosed by the deformed membrane and the plane O: This is obtained from (130) in which an  r ¼ 1 2eijk eab fia fjb rk : The domain is divided into zones—the shaded regions in Fig. 1—and the integral over each is estimated as the zone-centered value of the integrand multiplied by the shaded area, given by (144). Similarly, the scaled self-field at a given position y in the surrounding space is obtained by using (99), in which the integral is replaced by the sum of the integrals over the zones. Each of these is approximated by multiplying the value of the integrand at the relevant zone-centered point by the shaded area. The integrand is formed from zone-centered values of f and M and the averaged values of the nodal membrane position r. The time derivatives in (152) are approximated by the central differences: 1 s_ n ¼ ðs_ n þ 1=2 þ s_ n1=2 Þ, 2

1 s€ ¼ ðs_ n þ 1=2 s_ n1=2 Þ, h n

1 h

s_ n1=2 ¼ ðsn sn1 Þ,

ð153Þ

where h is the time increment and the node label (i,j) has been suppressed. Substitution into (152) furnishes the explicit, decoupled system: ðh1 þc=2Þmi,j s_ i,j,n þ 1=2 ¼ ðh1 c=2Þmi,j s_ i,j,n1=2 þ Si,j,n F i,j,n ,

si,j,n þ 1 ¼ si,j,n þ hs_ i,j,n þ 1=2 ,

ð154Þ

which is used to advance the solution in time node-by-node. The stress at zone-centered points is updated by using (64), (65), (104), in which the reactive constraint pressure q is computed from (67), (68) and ð115Þ2 . The starting procedure is derived from the quiescent initial conditions: _ i,j,0 ¼ 0, si,j,0 ¼ s0 ðui,j a Þ, s

ð155Þ

where s0 ðua Þ is assigned. Thus, from (154) we obtain ð2=hÞmi,j s_ i,j,1=2 ¼ Si,j,0 F i,j,0 ,

ð156Þ

in which the right-hand side is determined by the function s0 : The system is non-dimensionalized and the solution is advanced i,j,n to the first tn such that max F i,j,n j o d, a suitable tolerance. i,j jS We remark that because only long-time limits of solutions are relevant, temporal accuracy is not an issue. Stability is addressed by using sufficiently small time steps selected on the basis of successive trials based on a sequence of values of h. A similar temporal discretization is used to update the magnetization and director fields M and e at zone-centered points. Consistency with the derivation of the Lyapunov functions L1,2 of Section 3 requires the use of a staggered scheme in which the predicted position field at time step n þ1 is fixed while integrating (126). We then start the integration of ð126Þ2 using the value of e at step n as the initial condition (with the initial value e_ ¼ 0Þ. This calculation proceeds in increments of the time-like variable t2 . We fix the predicted value of e at the subsequent step and use this value to integrate ð126Þ1 with respect to t1 , using the value of M generated by the previous value of e as the initial condition (with _ ¼ 0Þ. This continues until convergence is achieved, initial value M yielding the magnetization associated with the predicted value of e: The integration with respect to t2 then resumes and the cycle is repeated until convergence is achieved, yielding the values of e and M associated with the position field at step n þ 1: The process is repeated until the deformation field converges, yielding the final equilibrium position, magnetization and director fields over all nodes and zone-centered points. However, numerical experiments indicate that this computationally intensive double-staggered scheme is not required in practice. Instead, we find that equilibrium states may be achieved by treating all fields on an equal basis as far as temporal integration is concerned. The magnetization at step n¼0 is set to zero. This is the unique solution to (112) if the applied field vanishes. Accordingly, the applied field intensity is first set to a small value and the equilibrium fields are obtained by the foregoing procedure. Successive equilibria are then computed for a sequence of increasing field intensities, using the equilibria associated with each member of the sequence as initial values for the next member.

5. Examples In this final section we discuss the results of some numerical experiments. All examples pertain to a membrane that is initially square, of side 8 mm and thickness h ¼ 50  103 mm: The latter is used in place of E in the formula (99) for the self-field, which was derived using a scheme in which E is interpreted as (dimensionless) thickness. The mass density is r ¼ 1750 kg=m3 ; the freespace permeability is m0 ¼ 4p  107 N=A2 (Newton per square Ampere) [14]; and the dipole source is centrally located above the plane at yd ¼ ð8 mmÞk. We find that convergence is achieved in all cases using a regular 33  33 mesh. Material parameters are taken to be those suggested in [34]. Thus, the saturation magnetization is M s ¼ m0 =2, the shear modulus is m ¼ 1:0  106 N=m2 , and the remaining parameters in (56) are C10 ¼ 1:0, C20 ¼ 0:625, C11 ¼ 2

0:0791, C21 ¼ 0:0, C01 ¼ b=6 and C02 ¼ b=2, where b ¼ m0 M s =2: Fig. 2 depicts the deformation of the membrane under zero pressure in response to a dipole of strength D ¼ 160  106 A m2 (cf. (20)). The vertical and in-plane dimensions are scaled differently to aid in visualization. We have used the data generated by the simulation, together with (69), to verify that the threedimensional principal stretches on the membrane surface are

M. Barham et al. / International Journal of Non-Linear Mechanics 47 (2012) 185–196

195

Fig. 2. Deformed membrane at D ¼ 160  106 A m2 .

Fig. 4. In-plane part, e, of director field, at D ¼ 160  106 A m2 .

Fig. 3. Referential magnetization at D ¼ 160  106 A m2 on the reference plane.

well within the limits required for the validity of ð41Þ2 . The referential in-plane magnetization field M is shown in Fig. 3. This field is directed everywhere toward the center of the membrane, where differentiability requires that it diminish to zero in intensity. This and the constraint (91) cause the interaction with the applied field to weaken near the center, resulting in a deformed surface that is relatively flat under the dipole source. Fig. 4 shows the variation of the in-plane part, e, of the director field with respect to position on the reference plane. The deformation deviates from Kirchhoff–Love kinematics wherever this is non-zero. This reflects the bias induced by the dipole source at points lying off the dipole axis, causing the director d on the deformed surface to tilt relative to the tangent plane as the membrane adjusts to the applied field. The effect diminishes near the corners of the membrane where the field is relatively weak, and near the center where the field lines intersect the membrane orthogonally and the associated bias vanishes; in either case the kinematics revert to the Kirchhoff–Love mode. Fig. 5 illustrates the self-field generated by the membrane, computed post-facto using (99), in a plane of symmetry obtained by fixing a reference coordinate at the value zero. Finally, the effects of pre-stretch and pressure are displayed in Fig. 6, in which the height of the deformed surface, at a point on the dipole axis, is plotted against dipole strength. The open circles and crosses correspond to zero applied pressure; the former corresponding to no pre-stretch and the latter to a uniform prestretch of 1.2 induced by an outward displacement of nodes on the boundary; these are subsequently fixed in the course of the simulation. Pre-stretch is seen to stiffen the membrane

Fig. 5. Self-field in space at D ¼ 160  106 A m2 , in the plane defined by u2 ¼ 0.

Fig. 6. Membrane displacement under the dipole source, as a function of dipole strength. Effect of pre-stretch indicated by circles (3) and crosses (); effect of fixed or volume-dependent pressure is indicated by dots () and stars (*), respectively (see text).

dramatically, resulting in a much smaller deflection at any given field strength. The effect of pressure (at no pre-stretch) is illustrated by the dotted and starred data, the former corresponding to a fixed inflation pressure P ¼ 2:0  105 Pa acting on the

196

M. Barham et al. / International Journal of Non-Linear Mechanics 47 (2012) 185–196

interior of the membrane; the external pressure is assumed to vanish. This is regarded as being supplied by a large reservoir with an opening on the reference plane. The stars correspond to a volume-dependent pressure in which the product of the pressure and the enclosed volume remains constant, as in an ideal gas at fixed temperature. The constant is derived by using (130) to compute the contained volume generated in response to the fixed pressure at zero field strength. As expected, pressure has a significant effect on deformation at small field intensities, but its relative importance diminishes with increasing intensity. Moreover, at any value of field intensity the volume-dependent pressure yields a smaller displacement than that produced by the fixed pressure. The discrepancy increases with field intensity due to the attendant increase in volume, which causes the volumedependent pressure to be reduced in magnitude. In all cases an upper limit is predicted for the deformation that can be maintained in equilibrium. Such limits are identified by the failure of the dynamic relaxation method to generate equilibria when the field intensity is increased above a critical value. Our results thus establish the existence of a limit-point instability at sufficiently high field intensities. This corroborates the analysis of [35], based on a low-order finite-dimensional projection of the model developed in [7].

Acknowledgments This work was performed under the auspices of the U.S. Department of Energy by the University of California, Lawrence Livermore National Laboratory under contract DE-AC5207NA27344. It was conducted during the tenure of a Lawrence Scholarship held by M. Barham from 2008 to 2011 in support of graduate studies at the University of California, Berkeley.

References [1] A. Dorfmann, R.W. Ogden, Magnetoelastic modeling of elastomers, Euro. J. Mech. A Solids 22 (2003) 497–507. [2] A. Dorfmann, R.W. Ogden, Nonlinear magnetoelastic deformations of elastomers, Acta Mech. 167 (2004) 13–28. [3] R. Bustamante, A. Dorfmann, R.W. Ogden, On variational formulations in nonlinear magnetoelastostatics, Math. Mech. Solids 13 (2008) 725–745. [4] R. Bustamante, A. Dorfmann, R.W. Ogden, Numerical solution of finite geometry boundary-value problems in nonlinear magnetoelasticity, Int. J. Solids Struct. 48 (2011) 847–883. [5] S.V. Kankanala, N. Triantafyllidis, On finitely strained magnetorheological elastomers, J. Mech. Phys. Solids 52 (2004) 2869–2908. ˜ oz, A model of the behaviour of magnetor[6] M.R. Jolly, J.D. Carlson, B.C. Mun heological materials, Smart Mater. Struct. 5 (1996) 607–614.

[7] M. Barham, D.J. Steigmann, M. McElfresh, R.E. Rudd, Finite deformation of a pressurized magnetoelastic membrane in a stationary dipole field, Acta Mech. 191 (2007) 1–19. [8] G. Gioia, R.D. James, Micromagnetics of very thin films, Proc. R. Soc. London A 453 (1997) 213. [9] D.J. Steigmann, Equilibrium theory for magnetic elastomers and magnetoelastic membranes, Int. J. Non-linear Mech. 39 (2004) 1193–1216. [10] G.A. Maugin, Continuum Mechanics of Electromagnetic Solids, NorthHolland, Amsterdam, 1988. [11] D.J. Steigmann, A concise derivation of membrane theory from three-dimensional nonlinear elasticity, J. Elasticity 97 (2009) 97–101. [12] M.E. Gurtin, An Introduction to Continuum Mechanics, Academic Press, Orlando, 1981. ¨ [13] C. Truesdell, R. Toupin, The classical field theories, in: S. Flugge (Ed.), Handbuch der Physik, vol. III/1, Springer, Berlin, 1960. [14] A. Kovetz, Electromagnetic Theory, Oxford University Press, 2000. [15] W.F. Brown, Magnetoelastic Interactions, Springer, Berlin, 1966. [16] A. DeSimone, P. Podio-Guidugli, On the continuum theory of deformable ferromagnetic solids, Arch. Ration. Mech. Anal. 136 (1996) 201. [17] D.J. Steigmann, On the formulation of balance laws for electromagnetic continua, Math. Mech. Solids 14 (2009) 390–402. [18] A. Nobili, A.M. Tarantino, Magnetostriction of a hard ferromagnetic and elastic thin-film structure, Math. Mech. Solids 13 (2009) 95–123. [19] M. Barham, D. White, D.J. Steigmann, Finite-element modeling of the deformation of magnetoelastic film, J. Comput. Phys. 229 (2010) 6193–6207. [20] R.D. James, Configurational forces in magnetism with application to the dynamics of a small-scale ferromagnetic shape memory cantilever, Continuum Mech. Thermodyn. 14 (2002) 55–86. ¨ [21] R.J. Knops, E.W. Wilkes, Theory of elastic stability, in: S. Flugge (Ed.), Handbuch der Physik, vol. VIa/3, Springer, Berlin, 1963. [22] R.W. Ogden, Non-linear Elastic Deformations, Dover, NY, 1997. [23] Q.-S. Zheng, Theory of representations for tensors: a unified invariant approach to constitutive equations, Appl. Mech. Rev. 47 (1994) 545. [24] R.W. Ogden, Waves in isotropic elastic materials of Hadamard Green or harmonic type, J. Mech. Phys. Solids 18 (1970) 149–163. [25] A.C. Pipkin, The relaxed energy density for isotropic elastic membranes, IMA J. Appl. Math. 36 (1986) 85. [26] E. Haseganu, D.J. Steigmann, Analysis of partly wrinkled membranes by the method of dynamic relaxation, Comput. Mech. 14 (1994) 596–614. [27] A. Atai, D.J. Steigmann, Coupled deformations of elastic curves and surfaces, Int. J. Solids Struct. 35 (1998) 1915–1952. [28] D.J. Steigmann, Eliza Haseganu’s analysis of wrinkling in pressurized membranes, in: A. Guran, A.L. Smirnov, D.J. Steigmann, R. Vaillancourt (Eds.), Advances in the Mechanics of Solids, in memory of E.M. Haseganu, vol. 15, World Scientific, Singapore, 2006, pp. 3–16. [29] D.J. Steigmann, Tension-field theory, Proc. R. Soc. Lond. A 429 (1990) 141–173. [30] S.A. Silling, Phase changes induced by deformation in isothermal elastic crystals, J. Mech. Phys. Solids 37 (1989) 293–316. [31] M. Taylor, D.J. Steigmann, Simulation of laminated thermoelastic membranes, J. Therm. Stresses 32 (2009) 448–476. [32] W. Hermann, L.D. Bertholf, Explicit Lagrangian finite-difference methods, in: T. Belytschko, T.J.R. Hughes (Eds.), Computational Methods for Transient Analysis, Elsevier, Amsterdam, 1983, pp. 361–416. [33] P. Underwood, Dynamic relaxation, in: T. Belytschko, T.J.R. Hughes (Eds.), Computational Methods for Transient Analysis, Elsevier, Amsterdam, 1983, pp. 245–265. [34] S.V. Kankanala, N. Triantafyllidis, Magnetoelastic buckling of a rectangular block in plane strain, J. Mech. Phys. Solids 56 (2008) 1147–1169. [35] M. Barham, D.J. Steigmann, M. McElfresh, R.E. Rudd, Limit-point instability of a magnetoelastic membrane in a stationary magnetic field, J. Smart Mater. Struct. 17 (2008) 6–11.