Numerical analysis of wrinkled, anisotropic, nonlinearly elastic membranes

Mechanics Research Communications 57 (2014) 1–5

Contents lists available at ScienceDirect

Mechanics Research Communications journal homepage: www.elsevier.com/locate/mechrescom

Numerical analysis of wrinkled, anisotropic, nonlinearly elastic membranes A. Atai a , D.J. Steigmann b,∗ a b

Department of Mechanical Engineering, University of Tehran, Tehran, Iran Department of Mechanical Engineering, University of California, Berkeley, CA 94720, USA

a r t i c l e

i n f o

Article history: Received 6 November 2013 Received in revised form 1 January 2014 Accepted 4 January 2014 Available online 13 January 2014

a b s t r a c t The method of dynamic relaxation is used to simulate wrinkling in anisotropic sheets modelling biotissues and structural membranes. © 2014 Elsevier Ltd. All rights reserved.

Keywords: Anisotropic membranes Tension-field theory Dynamic relaxation

1. Introduction In this work we extend the dynamic relaxation method, applied to the study of wrinkling in isotropic membranes in Haseganu and Steigmann (1994), to the analysis of partly wrinkled anisotropic membranes. The considered model is based on tension-field theory (Pipkin, 1994), which is derived from conventional membrane theory via quasiconvexification of the associated strain-energy function (Dacarogna, 1982), yielding the minimum strain energy that can be attributed to a given state of strain. The relaxed strain-energy function automatically excludes the destabilizing compressive stresses predicted by the original, unrelaxed, energy. It is derived by constructing an energy-minimizing sequence of deformations characterized by ever-more finely spaced wrinkles and effectively encodes the average energy density of a membrane containing many wrinkles. The same model emerges rigorously from three-dimensional nonlinear elasticity via the method of gamma convergence (Le Dret and Raoult, 1995), a technique for extracting the leading-order variational problem in the small-thickness limit. Whereas the expression for the potential energy in tension-field theory has been rigorously justified, the question of the existence of energy minimizers remains open. This is due to the failure of the relaxed strain-energy function to satisfy coercivity conditions underlying the hypotheses of available existence theorems based on the direct method of the calculus of variations (Dacarogna, 1982). This circumstance stems from the fact that the absence of

∗ Corresponding author. Tel.: +1 510 684 5380; fax: +1 510 642 6144. E-mail address: [email protected] (D.J. Steigmann). 0093-6413/$ – see front matter © 2014 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.mechrescom.2014.01.002

stress in connected regions of strain space – a feature of the relaxed energy – implies an absence of stiffness also. This situation furnishes impetus for methods such as that adopted here, in which the equilibrium problem is embedded in an artificial dynamical system constructed in such a way that the desired equilibria are globally asymptotically stable with respect to arbitrary initial data. Indeed, such an approach may prove fruitful in devising constructive existence theorems; however, this issue is beyond our present scope. We discretize the surrogate dynamical system using a finitedifference method based on Green’s theorem, and forward integrate in (artificial) time using a simple difference scheme. Temporal accuracy is not an issue, as it is only the asymptotic states that are of interest. This allows for the use of simple explicit finite difference operators to achieve an efficient vectorized system for computations. The basic theory and associated numerical analysis are well developed, and discussed in detail in the references cited. For this reason we forego their detailed description, emphasizing instead their application to anisotropic bio-elastic or artificial structural membranes. The strain-energy function adopted is inspired by the observed response of bio-tissues, in which an initially soft in-plane tensile response is followed by strain-stiffening as the underlying collagen fibers straighten and stretch on the micro-scale to accommodate an overall finite deformation. A similar mechanism characterizes woven PVC-clad polyester structural fabric used in tension structures, in which the interlaced fibers of the weave are initially crimped and hence relatively soft in response to inplane tension, stiffening as the fibers straighten in the transition to a stretching mode (Nadler et al., 2006). Alternative analyses of

2

A. Atai, D.J. Steigmann / Mechanics Research Communications 57 (2014) 1–5

wrinkling in anisotropic membranes are discussed in Woo et al. (2004), Barsotti and Vannucci (2013) for the case of small strains. In the present work we model finite strains, using Pipkin’s algorithm (Pipkin, 1994) to construct relaxed membrane theory. This is applicable if the underlying (unrelaxed) energy is a convex function of strain – a condition that is usually satisfied in applications. An interesting alternative treatment is described in Epstein (1999).

We have provided a brief resumé of the theoretical background in a companion paper (Atai and Steigmann, 2012). This is summarized here for the benefit of the reader. The equilibrium equation for membranes has a simple divergence structure identical to that for conventional bulk continua. We use a plane  as reference configuration, taken here to be unstressed. Let e˛ ; ˛ = 1, 2, be orthonormal vectors spanning (the translation space of) , and let e3 = e1 × e2 . We use the usual summation convention with Greek indices taking values in {1, 2} and Latin indices in {1, 2, 3}, while Greek subscripts preceded by commas are used to denote partial derivatives with respect to initial Cartesian coordinates x˛ . The equilibrium equation in the absence of lateral loads is or Ti˛,˛ = 0,

(1)

where T = Ti˛ ei ⊗ e˛

(2)

is the Piola stress (resultant) and div is the two-dimensional divergence operator on  . In the case of elasticity the stress is determined by the deformation gradient F = Fi˛ ei ⊗ e˛ ,

where Fi˛ = ri,˛

∂ri ≡ ∂x˛

(3)

in which ri are the Cartesian coordinates of a material point after deformation. The relationship is (Haseganu and Steigmann, 1994) T = WF ,

or Ti˛ =

∂W , ∂Fi˛

(4)

where W is the strain energy per unit area of  . These equations are augmented by traction data ti = Ti˛ ˛ or position data ri = Ri on complementary parts of the boundary, where ˛ are the components of the exterior unit normal to an edge in the reference configuration. For the mixed zero-load/position boundary-value problems considered here, equilibria furnish energy minimizers only of their gradients F satisfy the Legendre-Hadamard condition a ⊗ b · M(F)[a ⊗ b]≥0,

(5)

pointwise in , for all non-zero three-vectors a and two-vectors b ∈, where M(F) = WFF

(6)

is the tensor of elastic moduli. We shall also make use of the straindependent elastic moduli C(E), where E=

1 t (F F − 1); 2

E˛ˇ =

M(F)[A] = AS +

1 FC(E)[At F + Ft A] 2

1 (F F − ı˛ˇ ) 2 i˛ iˇ

(7)

is the strain in which 1 is the identity for 2-space, and C(E) = UEE

(8)

are the plane-stress elastic moduli, in which U(E) = W (F) is the associated strain-energy function.

(9)

(10)

for any tensor A of the form A = Ai˛ ei ⊗ e˛ , where S = UE

(11)

is the symmetric (plane) second Piola-Kirchhoff (S = S˛ˇ e˛ ⊗ eˇ ), related to the Piola stress by

2. Membrane theory

divT = 0,

An application of the chain rule furnishes the useful connection

T = FS;

Ti˛ = Fiˇ Sˇ˛ .

stress

(12)

We observe, using the minor symmetries of C, that A · M(F)[A] = At A · S + At F · C(E)[At F].

(13)

In terms of components, Mi˛jˇ Ai˛ Ajˇ = Ai˛ Aiˇ Sˇ˛ + C˛ˇ Ai˛ Fiˇ Aj Fj .

(14)

The strain-energy function W is locally convex (as a function of F) if and only if A · M(F)[A]≥0 for all A. Similarly, the strainenergy function U is locally convex (as a function of E) if and only if B · C(E)[B]≥0 for all B. It follows that if U is convex in E and the stress is positive semi-definite, then W is convex. This can be seen imme2 diately by using the spectral representation S = S u ⊗ u˛ , ˛=1 ˛ ˛ where S˛ are the principal stresses and u˛ are the associated principal directions; thus, if the principal stresses are non-negative, then At A · S =

2   2 S˛ Au˛  ≥0

(15)

˛=1

for all non-zero A. In the present context this observation is due to Pipkin (1993). In fact it is well known (Haseganu and Steigmann, 1994) in membrane theory that positive semi-definiteness of the 2nd PiolaKirchhoff stress is a necessary condition for the Legendre-Hadamard inequality and hence necessary for the energy to be minimized. This observation is the basis of tension-field theory, in which the membrane strain-energy function is relaxed so as to ensure that a negative-definite or indefinite stress never arises in any configuration (Pipkin, 1994). 3. Solution procedure We use a discrete version of the Green-Stokes theorem to discretize the equations directly on the computational plane . This method, incorporating position and traction boundary conditions, is described comprehensively by Haseganu and Steigmann (1994), to which reference may be made for a detailed discussion. The associated nodal points are distributed at the intersections of a curvilinear grid of boundary-fitted coordinates computed using the grid-mapping procedure discussed in (Wang and Steigmann, 1997) and so again we omit the details. This scheme is used to discretize the artificial dynamical system ˙ divT = ¨r + c r,

(16)

obtained by appending an explicit viscous damping term to the inertial term in the actual equations of motion. In turn, the spatially discretized system is used together with an explicit central-difference time integration scheme designed for efficient vectorization of the equations. This is not the actual equation of motion. It is an artificial system introduced solely to expedite the computation of equilibria.

A. Atai, D.J. Steigmann / Mechanics Research Communications 57 (2014) 1–5

3

Eq. (4) may be used to show that the total mechanical energy E, defined by E=

   1  2  r˙ + W da 

(17)

2

for the null traction/assigned position boundary-value problem, satisfies (Atai and Steigmann, 1998) d E=− dt



 2

c r˙  da

(18)



and hence decays on solution trajectories of (16), reaching its minimum, assuming one exists, in a static configuration. As shown in Atai and Steigmann (2012) the convexity of W(F), which is identically satisfied by the present relaxed formulation of membrane theory (Section 4), ensures that any static configuration is in fact a minimizer of the total strain energy. Accordingly, the total mechanical energy furnishes a Lyapunov function for the dynamical system. The method is implemented using quiescent initial conditions (assigned initial position field and vanishing initial velocity), yielding equilibria in the long-time limit of the solution to the artificial dynamical system. Thus temporal accuracy is not an issue, whereas temporal stability is ensured by repeating the computations, as necessary, using successively smaller time steps. Some of the examples discussed in Section 4 entail suturing of parts of the membrane boundary to form a curve. This is achieved using an algorithm developed in detail in (Atai and Steigmann, 1998). In all examples the equations are non-dimensionalized by a (positive) modulus C and the solution is advanced in   (artificial) time until the maximum value of local norm C −1 divT falls below a specified tolerance, on the order of 10−6 , corresponding to the approximate attainment of mechanical equilibrium. Here the divergence operator is non-dimensionalized by a problem-dependent length scale. 4. Pipkin’s theorem and its implementation ˆ Let S(E) be the 2nd Piola-Kirchhoff stress computed from a conventional plane-stress strain-energy function U(E) (cf. 11). In typical applications this function is convex, so that B · C(E)[B]≥0 for all B; and stationary at zero strain, so that S(E) = 0. Then, according to Pipkin’s theorem (Pipkin, 1994), the constitutive equation of relaxed membrane theory is delivered by the following algorithm, based on a partition of strain space into sets S, T and U defined by S = {E : E is negative semi − definite}, ˆ is positive semi − definite}, T = {E : S(E)

(19)

U = {E : E not in S or T }. If the value of E delivered by the explicit time-wise numerical integration procedure satisfies E ∈ S; i.e., if E11 ≤ 0, E22 ≤ 0 and 2 , then set S = 0; the membrane is slack. If E ∈ T; i.e., if E11 E22 ≥E12 ˆS11 ≥0, Sˆ 22 ≥0 and Sˆ 11 Sˆ 22 ≥Sˆ 2 , then set S = S(E); ˆ the membrane is 12

biaxially stressed and the original constitutive function is operative. Otherwise; i.e. for E ∈ U, find (> 0), a and a unit vector t = t˛ e˛ such that ˆ ∗ ), where E∗ = E + a2 n ⊗ n t ⊗ t = S(E

and

n = t × e3 .

(20)

Fig. 1. Initial configuration of annular membrane.

The implementation of this algorithm is illustrated for an orthotropic membrane having material axes that coincide with e˛ . It is straightforward to apply standard coordinate transformation rules to accommodate axes of orthotropy that are rotated relative to these axes (see Barsotti and Vannucci, 2013), but for the sake of brevity we refrain from doing so here. 4.1. A quadratic form in the strains Our model relies on a dimensionless, homogeneous, quadratic and convex function Q(E) of the strain that exhibits orthotropic symmetry with respect to the e˛ -axes. Thus, 2 2 2 Q = b1 E11 + b2 E22 + 2b3 E11 E22 + b4 E12 ,

where b1 –b4 are dimensionless constants, which is convex if and only if b1 > 0,

b2 > 0,

b1 b2 > b23

and

b4 > 0.

(23)

4.2. Exponential strain-energy function Biological tissues and structural fabric share the feature that their tensile responses are quite soft at small strains and stiffen dramatically as the underlying fibrous microstructure unbends and subsequently stretches in response to continued macroscopic extension. To model this behavior we adopt a strain-energy function of Fung type (Holzapfel et al., 2000; Alhayani et al., 2013; Haughton and Merodio, 2009), given by U=

1 C[exp(Q ) − 1], 2

(24)

where C is a positive modulus and Q is given by (22); the former is used to non-dimensionalize the equations (Section 3). We observe that U (Q) > 0 and U  (Q ) > 0, so that U is a convex and increasing function of Q . Because Q(E) is convex, it follows that U is a convex function of E and hence that (24) satisfies the hypotheses of Pipkin’s theorem. To see this we compute

For strains in U the membrane is wrinkled in the direction orthogonal to the (positive) uniaxial stress along the line panned by t (in the reference plane). Pipkin (1994) has proved that (20) is uniquely solvable for , a2 and t (apart from the irrelevant sense of t) under the stated hypotheses on U(E). Finally, for E ∈ U we set

UE = U  (Q )QE

S = t ⊗ t.

in which both terms on the right are non-negative.

(21)

(22)

and UEE = U  (Q )QEE + U  (Q )QE ⊗ QE ,

(25)

yielding B · C(E)[B] = U  (Q )B · QEE [B] + U  (Q )(B · QE )2 ,

(26)

4

A. Atai, D.J. Steigmann / Mechanics Research Communications 57 (2014) 1–5

Fig. 4. Initial configuration of a membrane with circular external boundary containing an elliptical hole. Fig. 2. Deformed configuration with a 20-degree counterclockwise rotation of the inner boundary. The material parameters are b1 = 0.2655, b2 = 0.7966, b3 = 0.0664 and b4 = 0.13.

Eqs. (11 and 25) furnish (Sˆ 11 /C)e−Q = b1 E11 + b3 E22 , and

(Sˆ 22 /C)e−Q = b3 E11 + b2 E22

(Sˆ 12 /C)e−Q = b4 E12 .

(27)

Using n1 = cos ,

n2 = sin ;

t1 = − sin

and t2 = cos ,

(28)

we cast (20) as the system 2

2

(/C) exp(−Q )sin

= b1 (E11 + a2 cos2 ) + b3 (E22 + a2 sin ),

(/C) exp(−Q )cos2

= b3 (E11 + a2 cos2 ) + b2 (E22 + a2 sin ),

(29)

2

5. Examples

−(/C) exp(−Q ) sin cos = b4 (E12 + a sin cos ). 2

a2

To solve this system we eliminate (/C) exp(− Q) and from the first two equations and substitute into the third. This yields the single nonlinear equation f( ) = 0, where f ( ) = b1 b4 E12 cos4 − [b3 (b3 + b4 )E22 + b1 (b4 E11 − b2 E22 )] × cos3 sin − [(b1 b2 − b23 − b3 b4 )E11 − b2 b4 E22 ] × cos sin3 − b2 b4 E12 sin4 ,

in which the E˛ˇ are fixed by the numerical integration procedure. Similar transcendental equations emerge in the small-strain models described in Barsotti and Vannucci (2013), Woo et al. (2004). In fact, because of the properties of the exponential function, precisely the same equation arises if the present strain-energy function is replaced by Q itself, yielding a model that is appropriate if the strains are sufficiently small. Because of the irrelevance of the sign of t in (20), a solution is physically the same as + , and so we confine attention to the domain 0 ≤ < . Using a simple interval-splitting method, we find that there is typically more than one root of f( ); but, in accordance with Pipkin’s theorem, that there is only one that delivers positive values of /C and a2 .

(30)

Fig. 3. Deformed configuration with a 20-degree counterclockwise rotation of the inner boundary. The material parameters are b1 = 1.001, b2 = 0.1001, b3 = 0.01002 and b4 = 0.1.

We consider two basic problems. In the first, an annular membrane (Fig. 1) is deformed by fixing it at the outer boundary and assigning a displacement at the inner boundary corresponding to a counterclockwise rotation of 20◦ . The radius of the outer boundary is 2.5 that of the inner boundary and the latter is used as the length scale for non-dimensionalization. The predicted deformation is seen to inherit the four-fold symmetry of the underlying axes of orthotropy. Wrinkled regions, corresponding to the subdomain U of strain space, are indicated by dashed segments aligned

Fig. 5. Deformed configuration of the sutured membrane. The material parameters are b1 = 0.2655, b2 = 0.7966, b3 = 0.0664 and b4 = 0.13.

A. Atai, D.J. Steigmann / Mechanics Research Communications 57 (2014) 1–5

5

symmetry of the problem yields a configuration of the suture line, which is not specified a priori, that coincides with the horizontal axis of material symmetry. The solution exhibits a combination of tense, wrinkled and slack (unstressed) regions; the latter corresponding to the subdomain S of strain space. Wrinkled and tense regions are indicated as before, whereas slack regions are identified by small open circles. Figs. 5 and 6 depict the solution in the vicinity of the right-hand end of the suture line, and indicate that the distributions of wrinkled and slack states are strongly sensitive to material parameters. The largest principal stretch in the example of Fig. 5 is 1.12976 and occurs at the material point labelled A in Fig. 4. In Fig. 6 the largest principal stretch is 1.15912 and occurs at point B. These results indicate that the material is deformed in the moderate finite-strain range. Acknowledgments DJS gratefully acknowledges the support of the Powley Fund for ballistics research. References

Fig. 6. Deformed configuration of the sutured membrane. The material parameters are b1 = 1.001, b2 = 0.1001, b3 = 0.01002 and b4 = 0.1.

locally with the principal axes of stress in the current configuration (Figs. 2 and 3). These directions are parallel to the vector field Ft, where t is the axis of (positive) uniaxial 2nd Piola-Kirchhoff stress (cf. (21)). Elsewhere the stress is positive definite, corresponding to the subdomain T of strain space. Two sub-examples are shown, based on different sets of the material parameters b1 –b4 in (22). These indicate that the wrinkling pattern is somewhat sensitive to material parameters. Our results are in qualitative agreement with those obtained in Woo et al. (2004), Barsotti and Vannucci (2013) for the small-strain theory and in (Haseganu and Steigmann, 1996) for the nonlinear theory of orthotropic networks. In the second problem we analyze the suturing of a membrane containing an elliptical hole (Fig. 4). This problem was studied in Atai and Steigmann (2012) for the case of isotropic membranes. The exterior boundary is fixed while the upper and lower edges of the interior elliptical boundary are sutured together using the algorithm developed by Atai and Steigmann (1998). The

Alhayani, A.A., Giraldo, J.A., Rodríguez, J., Merodio, J., 2013. Computational modeling of bulging of inflated cylindrical shells applicable to aneurysm formation and propagation in arterial wall tissue. Finite Elem. Anal. Des. 73, 20–29. Atai, A., Steigmann, D.J., 1998. Coupled deformations of elastic curves and surfaces. Int. J. Solids Struct. 35, 1915–1952. Atai, A., Steigmann, D.J., 2012. Modeling and simulation of sutured biomembranes. Mech. Res. Commun. 46, 34–40. Barsotti, R., Vannucci, P., 2013. Wrinkling of orthotropic membranes: an analysis by the polar method. J. Elast. 113, 5–26. Dacarogna, B., 1982. Quasiconvexity and relaxation of nonconvex problems in the calculus of variations. J. Funct. Anal. 46, 102–118. Epstein, M., 1999. On the wrinkling of anisotropic elastic membranes. J. Elast. 55, 99–109. Haseganu, E., Steigmann, D.J., 1994. Analysis of partly wrinkled membranes by the method of dynamic relaxation. Comput. Mech. 14, 596–614. Haseganu, E., Steigmann, D.J., 1996. Equilibrium analysis of finitely deformed elastic networks. Comput. Mech. 17, 359–373. Haughton, D., Merodio, J., 2009. The elasticity of arterial tissue affected by Marfan’s syndrome. Mech. Res. Commun. 36, 659–668. Holzapfel, G.A., Gasser, T.C., Ogden, R.W., 2000. A new constitutive framework for arterial wall mechanics and a comparative study of material models. J. Elast. 61, 1–48. Le Dret, H., Raoult, A., 1995. The nonlinear membrane model as a variational limit of nonlinear three-dimensional elasticity. J. Math. Pures Appl. 75, 551–580. Nadler, B., Papadopoulos, P., Steigmann, D.J., 2006. Convexity of the strain-energy function in a two-scale model of ideal fabrics. J. Elast. 84, 223–244. Pipkin, A.C., 1993. Convexity conditions for strain-dependent energy functions for membranes. Arch. Ration. Mech. Anal. 121, 361–376. Pipkin, A.C., 1994. Relaxed energy densities for large deformations of membranes. IMA J. Appl. Math. 52, 297–308. Wang, M., Steigmann, D.J., 1997. Small oscillations of finitely deformed elastic networks. J. Sound Vib. 202, 619–631. Woo, K., Igawa, H., Jenkins, C.H., 2004. Analysis of wrinkling behavior of anisotropic membrane. Comput. Model. Eng. Sci. 6, 397–408.