- Email: [email protected]
Analysis of geometrically nonlinear anisotropic membranes: theory and verification
Finite Elements in Analysis and Design 41 (2005) 963 – 988 www.elsevier.com/locate/finel
Analysis of geometrically nonlinear anisotropic membranes: theory and verification Wenqing Zhang, John L. Leonard, Michael L. Accorsi∗ Department of Civil and Environmental Engineering, University of Connecticut, Storrs, CT 06269-2037, USA Received 30 August 2004; received in revised form 16 November 2004; accepted 29 December 2004
Abstract A new anisotropic constitutive model is introduced to simulate the geometrically nonlinear dynamic behavior of general anisotropic membranes experiencing large deformations. This model relates the second Piola–Kirchhoff stress with the Green–Lagrange strain in a 3-D convected curvilinear coordinate frame. The construction of the anisotropic constitutive law from the predefined principal material coordinate system to the final convected curvilinear coordinate system using coordinate transformations is presented. The proposed theory is implemented into a finite element code and several numerical examples are given for validation. An application problem is described to demonstrate the importance of the new model in real applications. 䉷 2005 Elsevier B.V. All rights reserved. Keywords: Anisotropic/orthotropic constitutive law; Geometrically nonlinear dynamic behavior; Membrane and fabric structures; Finite element method
1. Introduction Membrane and fabric structures are used in a wide variety of applications. Several examples are airbags in the automobile industry [1], parachutes and parafoils in airdrop applications [2], and pneumatic muscle actuators that are used in the robotics industry [3,4] and more recently for parachute control [5]. In many such applications, the membrane structure undergoes large rotations and displacements, but not necessarily large strains. For these cases, a constitutive law that accounts for geometrically nonlinear motion is required. ∗ Corresponding author. Tel.: +1 860 486 5642; fax: +1 860 486 2298.
E-mail address: [email protected] (M.L. Accorsi). 0168-874X/$ - see front matter 䉷 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2004.12.003
964
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Much research has been performed in modeling the nonlinear anisotropic behavior of fabric membranes using a micromechanics approach. For example, Stubbs and Thomas developed a nonlinear elastic constitutive model for coated fabrics in 1984 [6]. In 1997, Sun et al. studied the stress and strain of each individual yarn in a weave [7] and Kato et al. conducted research to model the constitutive properties of weaves from the concept of a fabric lattice [8] in 1999. In these micromechanical fabric models, detailed mechanics of the fabric, such as the effects of coatings, decrimping, friction and slip are modeled. Other researchers have concentrated on macroscopic models for orthotropic materials. Kyriacou et al., for example, presented a finite element code to analyze nonlinear orthotropic and hyper-elastic membranes in large deformation problems in 1996 [9]. The membrane strain energy is expressed in terms of the invariants of the Cauchy–Green strain tensor and the orthotropic properties are pre-defined in the “preferred directions”. In 1999, Manch and Rio established another orthotropic constitutive model in a convected coordinate frame [10] where the strain energy is calculated from the Almansi strain tensor expressed in the convected configuration. In this work, a new macroscopic model is introduced to simulate the geometrically nonlinear dynamic behavior of general anisotropic membranes using a Total Lagrange formulation. A convected curvilinear coordinate frame is used to describe the large deformation and the geometrically nonlinear behavior. General convected anisotropic constitutive laws, which relate the second Piola–Kirchhoff stress to the Green–Lagrange strain, are established in the curvilinear coordinate system on the basis of the material properties predefined in the membranes’ principal material directions. A finite element code has also been developed to implement the anisotropic constitutive matrix for membranes in a general structural dynamics analysis. Several numerical examples are given to validate the capabilities of the model proposed in this research by comparing the dynamic behavior of initially orthotropic/anisotropic membrane structures with mechanically identical woven/braided fiber networks. For woven networks, the fibers are initially orthogonal, whereas for braided networks, the fibers are initially nonorthogonal. Finally, an application of anisotropic materials in parachute design is also presented. In these example problems, a very simple fabric micromechanical model is used to facilitate verification of the results, as discussed later. Specifically, the effects of fiber undulation and decrimping, fabric coating, and friction and slippage are not modeled in these examples. The inclusion of these fabric micromechanical effects into the macroscopic anisotropic model presented here will be considered in future work.
2. Methodology Anisotropic constitutive laws and three-dimensional kinematics to be used in a finite element model [11] are developed in this section from considerations of the generalized Hooke’s law relating stresses and strains [12]. When studying the dynamic behavior of anisotropic membranes, the following two assumptions are made: (1) the membrane material is assumed to be linear elastic and (2) only plane stresses are considered, that is, the stresses in the normal direction of the membranes are assumed to be zeros. In this part, the theoretical kinematic formulations are first developed, followed by the construction of the anisotropic constitutive laws from the principal material coordinate system to the convected curvilinear coordinate system by using coordinate transformation. The incorporation of the anisotropic constitutive theory into an existing finite element code is also presented. In the remaining part of this paper, Latin indices take on the values ranging over 1, 2, 3, while Greek indices take on the values of 1, 2. Repeated indices, unless otherwise noted, imply summation over the
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
965
Fig. 1. Lagrangian description of continuum membrane element.
range of the indices. The arrow above a variable indicates the variable is a vector. The tilde symbol above a variable denotes the variable is a tensor. 2.1. Kinematics A total Lagrangian description [13] is employed in this work to describe the deformation of a continuum solid membrane element. That is, every material particle position is referred to the undeformed configuration. A material particle in the initial undeformed configuration is identified in a fixed global Cartesian coordinate system as = Xi ei X
(1)
and its position in the deformed configuration at time t can be described as x = xi ei ,
(2)
where ei are the unit vectors along the coordinate axis of the global Cartesian coordinate system, Xi the initial coordinate of the particle along the ith axis, and xi the current coordinate of the particle along the ith axis, as shown in Fig. 1. Because of the assumption of plane stress, a 2-D curvilinear coordinate system is defined for every element at time t to describe the membrane element deformation in the convected frame. The covariant base vectors are g =
jx j
=
jxs j
es = gs es ,
(3)
966
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
= jX = jXs es = Gs es . G j
j
(4)
for the initial configuration. The Obviously, g are the base vectors for the current configuration and G covariant components of the metric tensor with respect to the initial and the current configuration are expressed as the following equations, respectively: m · G = Gm G = G G ,
(5)
g = g · g = gm gm ,
(6)
where the dot symbol represents dot product operation. According to Green and Zerna [14], the Green– Lagrange strain tensor is therefore expressed as ⊗ G ˜ = 21 (g − G )G
(7)
with the components = 21 (g − G )
(8)
and ⊗ denotes the dyad of two vectors. The Green–Lagrange strain is defined with respect to the initial configuration, and relates to the second Piola–Kirchhoff stress by the following equation [15]: ˜ , S˜ = C˜
(9)
where S˜ denotes the second Piola–Kirchhoff stress tensor, and C˜ is a fourth-order tensor representing the stress–strain relationships. Once the appropriate expressions of C˜ for the anisotropic material are obtained, the internal forces can be calculated by using Eq. (9), and a stiffness matrix can be also constructed. More detailed explanations are presented in a later section. 2.2. Description of anisotropic material property This section presents the transformation procedure for the anisotropic constitutive laws in the finite element analysis. Three coordinate systems, namely, a principal material coordinate system, a rectangular Cartesian coordinate system and a curvilinear coordinate system, are used to describe the anisotropic material properties of membrane elements. At first, the anisotropic constitutive matrix is defined in the principal material coordinate system, assuming the initial state of the plane sheet of material is flat. Then, a transformation of the constitutive matrix from the principal material coordinate system to the three-dimensional frame of Cartesian coordinate system is performed. After another transformation of constitutive matrix from the Cartesian coordinate system to the curvilinear coordinate system of the element is done, the anisotropic material properties expressed in the element’s curvilinear coordinate system are finally obtained, which are used in the calculations of the finite element analysis later. 2.2.1. Principal material coordinate When studying membranes, it is possible to divide a large membrane into sectors that are arbitrarily oriented in three-dimensional space, with every sector containing one or many isoparametric quadrilateral or triangular membrane elements [16]. These various sectors are assumed to be flat initially, and can be described with anisotropic or orthotropic materials oriented in the principal material directions within
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
967
Fig. 2. Three coordinate systems.
them. Therefore, in this study, each sector is considered as a plane sheet in its initial flat shape and a principal material coordinate system is defined in it. The anisotropic constitutive matrix is then developed in this coordinate system to describe the anisotropic material properties of the elements in the plane sheet sector. 2.2.1.1. Definition of the principal material coordinate system The material properties of a plane sheet are initially described in the principal material coordinate system, assuming that the initial state of the plane sheet is flat. The principal material coordinate system is defined as an orthogonal coordinate system such that two of its axes are mapped to the principal material directions in the plane sheet. As shown in Fig. 2, the unit axial vector m 1 and m 2 are orthogonal vectors in the principal directions on the initial flat sheet. They are determined by the global coordinates of three pre-selected nodes P1 , P2 and P3 on j the sheet, with position vectors R0i . The orthonormal direction cosines i relate m i to the Cartesian unit axis vectors ej by the following equation (no sum on i): j
m i = i ej = di /|di |
(10)
with d1 = R02 − R01 ,
(11.1)
968
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
d3 = d1 × (R03 − R01 ),
(11.2)
d2 = d3 × d1 ,
(11.3)
where the symbol × denotes cross-product operation. j The inverse directional cosines i are defined to relate unit vectors ei to unit vectors m j by j
j, ei = i m j
(12) j
where i are the inverse of i and are determined by ik kj = ij
(13) j
and ij is the Kronecker delta. Because i are orthonormal, Eq. (13) yields j
i = ij .
(14)
2.2.1.2. Description of anisotropic properties in principal direction coordinate system According to Jones [12], the 2-D anisotropic constitutive relations described in the principal material coordinate system can be expressed as
= S t ,
(15)
where and t represent the covariant components of the Green–Lagrange strains and the second Piola–Kirchhoff stresses in the principal material directions, respectively. S are the corresponding stress–strain coefficients with the following symmetry properties:
S = S = S = S = S .
(16)
In terms of physically measurable quantities E1 , E2 , 12 , 21 , G, 1 and 2 , the anisotropic compliances are 11 S11 = 1/E1 ,
(17.1)
22 = 1/E2 , S22
(17.2)
12 12 21 21 = S21 = S12 = S21 = 1/G, S12
(17.3)
11 S22 = − 12 /E1 ,
(17.4)
12 = − 21 /E2 , S21
(17.5)
12 21 11 11 S11 = S11 = S12 = S21 = 1 /E1 ,
(17.6)
12 21 22 22 = S22 = S12 = S21 = 2 /E2 , S22
(17.7)
1 where E1 and E2 stand for the Young’s modulus in the first and second principal material direction (m and m 2 directions), respectively. ij represents the Poisson’s ratio of transverse strain in thej th direction 11 equals to S 12 . G is the shear when stressed in the ith direction. Note that because 12 E2 = 21 E1 , S22 21
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
969
modulus in 1–2 plane. Since the plane studied is assumed in plane stress state, the shear modulus in 1–3 and 2–3 planes are ignored. i is the Lekhnitski coefficient [17] of mutual influence, which is defined as i = 12 /i .
(18)
It characterizes shearing in the plane caused by the normal stress in the ith direction. In some cases where G is difficult to determine, the Firt approximation can be utilized [18]: √
12 / 21 . G = E1 (19) √ 2 1 + 12 21 As a result, the contracted stress–strain relations for the plane sheet in the principal coordinate system can be obtained by inversion of Eq. (6), and be expressed as t = Q .
(20)
Thus, Q11 11 = E1 (1 − G 2 2 /E2 )/ ,
(21.1)
Q22 22 = E2 (1 − G 1 1 /E1 )/ ,
(21.2)
12 21 21 Q12 12 = Q21 = Q12 = Q21 = G(1 − 12 21 )/ ,
(21.3)
22 Q11 22 = Q11 = (E2 12 + G 1 2 )/ ,
(21.4)
21 11 11 Q12 11 = Q11 = Q12 = Q21 = −G( 1 + 12 2 )/ ,
(21.5)
21 22 22 Q12 22 = Q22 = Q12 = Q21 = −G( 2 + 21 1 )/ ,
(21.6)
= (1 − 12 21 ) − G/E1 [2 12 1 2 + 1 1 + (E1 /E2 ) 2 2 ].
(21.7)
where
For the special case of an orthotropic material sheet, 1 and 2 both equal zero. 2.2.2. Cartesian coordinate system Once the stress–strain relations are obtained in the principal material coordinate system, it is necessary to transform them to the Cartesian coordinate system. The Cartesian coordinate system employed here is identical to the fixed global Cartesian coordinate system mentioned before and has the unit axis vector ei . The 3-D relations between stresses and strains in this system can be expressed as ij
Trs = Crs Eij ,
(22)
where Trs represent the second Piola–Kirchhoff stress, Eij the Green–Lagrange strain components, both in the reference Cartesian coordinate frame. Considering Eq. (22) together with Eqs. (9) and (20), the components of the constitutive matrix expressed in the Cartesian coordinate system are given by ij
j
p
t Crs = im n Qmn pt r s .
(23)
970
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Because of the plane-stress state assumed, only the Q terms are non-zeros. Thus, Eq. (23) can be rewritten as ij
j
Crs = i Q r s .
(24)
2.2.3. Curvilinear coordinate system In order to obtain the final representation of constitutive laws in the current curvilinear coordinate system, another transformation is needed. The convected curvilinear coordinate system described previously based on every element integration point is also employed here at time t with the natural base vectors as described in Eq. (3). The contravariant base vectors are g = hs es ,
(25)
whose components hs are determined by g · g = hs gs =
(26)
and the contravariant components of the metric tensor are h = g · g = hm hm .
(27)
Finally, the constitutive relations in the element curvilinear system are expressed as S = E ,
(28)
where S denote the contravariant components of the second Piola–Kirchhoff stress in the element curvilinear coordinate system, the covariant components of the Green–Lagrange strain components. By comparing Eq. (28) with Eqs. (22) and (9), we see that the components of the constitutive tensor in the “integration-point-wise” curvilinear coordinate system can be calculated by
E = E h h
(29)
with
mn p t E = h m h n Cpt g g .
(30)
This completes the transformation of constitutive components from Q given in Eq. (20) to E in the curvilinear coordinate system at time t. The results obtained here are implemented in the finite element code to perform the anisotropic material analysis. 2.3. Finite element implementation The anisotropic constitutive laws derived previously have been incorporated into an existing finite element code to execute the nonlinear dynamic analysis of anisotropic membranes. A FORTRAN program named TENSION [11], which is developed for the 3-D nonlinear dynamic analysis of membrane, cable and payload structures, is used. In this section, the procedures of the implementation are described. As mentioned previously, once the current curvilinear coordinate system is determined, the Green– Lagrange strains are calculated from Eqs. (7) and (8). Note that the components gs of the covariant base
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
971
vectors are obtained from the derivatives of the element shape functions [16] in the program. Then, the components of the second Piola–Kirchhoff stress S are determined from Eq. (28) and are used in the calculation of the internal restoring force RN i and tangent stiffness matrix KMN ij [11]. RN i =
1 2
√ S (HM, HN, + HN, HM, )yMi A dX1 dX2
(31)
KMN ij =
1 j 1 i S (HM, HN, + HM, HN, ) + E (HR, HN, + HR, HN, )yRi 2 4 √ ×(HQ, HM, + HM, HQ, )yQj A dX1 dX2 , (32)
where RN i denotes the ith component of the internal force on node N, yMi the current coordinate of node M in the ith direction. KMN ij are the components of the tangent stiffness matrix, H represents the shape functions and comma in the subscript denotes the partial derivatives of the variable with respect to the subscript after the comma. These equations express the current internal force and stiffness at time t, but referred to the initial undeformed state. A more detailed study of the dynamic analysis method used in TENSION can be found in [11]. Therefore, a brief flow chart of the program using an implicit time integration scheme and full Newton–Raphson iteration with the geometrically nonlinear anisotropic constitutive laws can be expressed as follows: (i) Set up the model, input the anisotropic material properties in principal material coordinate system (ii) Transform the anisotropic material properties from the principal material coordinate system to Cartesian coordinate system (Transformation 1) (iii) Begin the time step loop (iv) Begin Newton–Raphson iteration (v) Element loop (vi) Calculate the curvilinear coordinate system and Green–Lagrange strain (vii) Transform the anisotropic material properties from Cartesian coordinate system to curvilinear coordinate system (Transformation 2) (viii) Calculate the second Piola–Kirchhoff stress, restoring forces and tangent stiffness matrix (ix) End element loop (x) Construct and solve system equations (xi) Convergence check, if not convergent, go to step (iv) (xii) End time step loop ij
The first transformation takes place only once and the calculation of Crs is the same for all the elements in a particular sheet of material. The second transformation takes place for every membrane element in all the iterations of every time step, because the deformation changes and therefore the curvilinear coordinate system of each element changes in time. Since the constitutive relation is evaluated within the Newton–Raphson iteration, additional nonlinear effects (such as material nonlinearities) could also be included in the current formulation.
972
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
3. Numerical examples Two numerical examples are presented to verify the constitutive model developed previously. In each example, the dynamic behavior of a membrane structure is compared to a simple woven/braided cable network that is mechanically identical to the membrane structure. The model is initially orthotropic in the first example and initially anisotropic in the second example. 3.1. Comparison of an initially orthotropic membrane structure with a mechanically identical cable network Consider the woven cable network shown in Fig. 3. The width and height of the network are w and h, respectively. The number of rows of horizontal cables and columns of vertical cables are both n. The vertical cables are woven with the horizontal cables according to the pattern shown in Fig. 4. A simple fabric micromechanical model is assumed for the cable network to facilitate comparison with the macroscopic anisotropic membrane model. It is assumed that the horizontal and vertical cables are flat (i.e. no undulation). It is further assumed that the horizontal and vertical cables are pinned at their junctions and that no rotational resistance or slippage occurs at the junctions. Under these assumptions, the properties of the mechanically equivalent anisotropic membrane can be determined analytically. If the cross-section areas of the vertical cables and the horizontal cables are (w/n)t and (h/n)t, respectively, then the cable network is equivalent to a membrane structure with the same dimensions, w
Fig. 3. A woven cable network.
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
973
Fig. 4. Initial weave pattern of the cable network.
Fig. 5. An initially orthotropic cable network and membrane model.
and h, and thickness t. The Young’s modulus of the horizontal cables and vertical cables are set as E1 and E2 , respectively. Assuming that the density of all cables is and there is no shear between cables, then the membrane structure, which has the same mass and material properties as the cable network, should
974
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Fig. 6. Deformed shape of the initially orthotropic membrane structure.
have a density of 2 and the following stress–strain relations in the Cartesian coordinate system: E1 0 0 C= 0 E2 0 . 0 0 0
(33)
Obviously, this membrane structure is initially orthotropic. The methodology of how to derive the corresponding material property for a membrane structure to be identical to a braided cable network is given in Appendix A. If the boundary conditions of the membrane structure and the cable network are set the same, these two models are then mechanically identical and they should present the same deformation when subjected to the same loads. In this example, the left edge and the bottom edge of both models are fixed, and a point force P with a diagonal magnitude of 0.4 pound is applied to the node on the upper right corner, as shown in Fig. 5. The width w and height h for these two models are both 1 in. In the membrane structure, the thickness t is 0.4 in and the density is 0.8 pound/in3 . In the cable network, n is equal to 50; the density is 0.4 pound/in3 ; the cross-section area for all the cables is 0.005 in2 ; and E1 and E2 are both 106 pound/in2 . The deformation shapes for both models are tracked to verify the theory. The change of the angles between the horizontal
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
975
Fig. 7. Deformed shape of the initially orthotropic cable network.
cable and the vertical cable at two randomly selected points in the cable network, points A and B, as in Fig. 5, is also compared with the corresponding change of the angles obtained from the membrane structure at the same points. Figs. 6 and 7 show the deformation shapes of the membrane structure and the cable network after the point force is applied for 0.4 s, respectively. The deformation shapes of these two models agree very well. Fig. 8 shows the angle changes obtained from both the membrane structure and the braided cable network at points A and B. The results demonstrate that, over a large range of deformation, the initially orthotropic membrane structure behaves the same as the mechanically identical cable network. 3.2. Comparison of an initially anisotropic membrane structure with a mechanically identical cable network In this example, the dynamic behavior of a membrane structure is compared with a mechanically identical braided cable network. Fig. 9 shows the initial weave pattern of the cable network in this example. The same simplifying assumptions used in the previous example for the fabric micromechanical behavior are also used here to facilitate comparison with the membrane model. As a result, a mechanically identical membrane structure not only has the same dimensions and mass as the cable network, but also has the
976
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Fig. 8. Angle change versus time for the initially orthotropic membrane structure and cable network.
Fig. 9. Initial weave pattern of the anisotropic cable network.
stress–strain relations in the Cartesian coordinate system as
0.0045(E1 + E2 ) C = 0.0625(E1 + E2 ) 0.0335(E1 − E2 )
0.0625(E1 + E2 ) 0.8705(E1 + E2 ) 0.4665(E1 − E2 )
0.0335(E1 − E2 ) 0.4665(E1 − E2 ) . 0.1250(E1 + E2 )
(34)
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
977
Fig. 10. An initially anisotropic membrane structure and cable network model.
Fig. 11. Deformed shape of the initially anisotropic membrane at 0.26 s.
Assuming the Young’s modulus for cables in directions 1 and 2, which are marked in Fig. 9, are E1 and E2 , respectively. More detailed derivation for this constitutive matrix can be found in Appendix A. Eq. (34) indicates that the membrane structure initially possesses anisotropic material properties. Fig. 10 shows the initial shapes and the boundary conditions of the membrane structure and the cable network model used in this example. As can be seen, the left edges of both models are fixed, and the right edges are subject to a uniform edge force F with a magnitude of 8.1 lb/in. One point, point C, is also randomly selected to compare the dynamic behavior of the membrane structure and the braided cable network. During the analysis, the length of the horizontal edge w is 1 in and the height h is 0.5 in for both models. In the membrane structure, the thickness t is 0.1 in and the density is 0.8 pound/in3 . In the
978
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Fig. 12. Deformed shape of the initially anisotropic cable network at 0.26 s.
Fig. 13. X Displacement comparison of Point C.
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Fig. 14. Y Displacement comparison of Point C.
Fig. 15. Angle change at Point C.
979
980
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Fig. 16. Initial mesh of a quarter-scale C9 parachute model.
cable network, n is equal to 80; the density is 0.4 pound/in3 ; the cross-section area for all the cables is 0.00125 in2 and E1 and E2 are both 100 pound/in2 . Figs. 11 and 12 show the deformed shapes at a time of 0.26 s for the membrane structure and the cable network, respectively. The displacements of point C obtained from the two models are given in Figs. 13 and 14, and the change of the angles between intersecting cables/fibers at point C is also given in Fig. 15. All these results indicate that the two models demonstrate the same dynamic behavior.
4. Applications The anisotropic theory developed here has been applied to the simulation of dynamic parachute problem. In this example, a quarter-symmetric model of a round C9 parachute is used and the inflation of the parachute from the initially unstressed state to the totally inflated state is simulated. The finite element mesh of the initial configuration is shown in Fig. 16. This model consists 88 three-node cable elements
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
981
Fig. 17. Fully inflated shape of a quarter-scale C9 parachute model.
for the cables in canopy, 32 two-node cable elements the suspension lines, 1 payload element and 154 nine-node membrane elements for the canopy, which is made up of 7 gores. Three kinds of materials, i.e., isotropic, orthotropic and anisotropic materials, are used separately as the canopy fabrics, and the stress distributions along the centerline of one gore are compared. For all cases, the true (Cauchy) stress is presented. In this example, the isotropic material used has the Young’s modulus as 4.32 × 105 lb/ft2 , Poisson’s ratio as 0.3; the orthotropic material used has the Young’s modulus of 4.32 × 105 lb/ft 2 in one direction and 2.16 × 105 lb/ft2 in another direction; the anisotropic material is set to be weaved up by two directional fibers with the Young’s modulus as 4.32 × 105 lb/ft2 , and the intersecting angle is 60◦ . The parachute canopy undergoes 0.5 lb/ft 2 pressure. The canopy’s thickness is 0.1 ft and the payload is 50 lb for all cases. The membrane thickness is 0.0001 ft and the cable cross-section area is 0.0001 ft 2 . Fig. 17 represents the inflated configuration of the model at 6.0 s. Although three different materials are used for the canopy fabrics, the deformation shape of the parachute does not change much. However, the maximum and the minimum principal stress distributions differ significantly. Figs. 18 and 19 give the magnitudes of the maximum principal stress and minimum principal stress along the centerline of a canopy gore. Figs. 20–25 show the top view of distributions for the maximum and minimum principal
982
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Fig. 18. Maximum principal stress for different materials.
Fig. 19. Minimum principal stress for different materials.
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
983
Fig. 20. Maximum principal stress vector of isotropic canopy.
stress vectors, respectively. As can be seen, when anisotropic membrane is used, the maximum principal stresses on the canopy are all in circumferential direction and the minimum principal stresses are all meridian stresses. When isotropic material or orthotropic material is used, the maximum principal stresses are circumferential stress around the canopy center and change to meridian stresses with the increase of radius, while the minimum principal stresses are in radial direction around the canopy center and change to hoop stresses with radius increase. Note that the radius where the principal stresses begin to change direction is different for isotropic canopy and orthotropic canopy. Figs. 18 and 19 show that the magnitudes of principal stresses are different when different materials are used.
5. Conclusions In this research, the theoretical formulation, finite element implementation and numerical validation of a new anisotropic constitutive law are presented to describe the geometrically nonlinear dynamic behavior of general anisotropic membranes with large deformations. The theory is based on a Total Lagrange formulation using the second Piola–Kirchhoff stress and the Green–Lagrange strain. A 3-D convected curvilinear coordinate frame is employed to describe the large deformation and the geometrically nonlinear behavior of membranes. The methodology of how to derive the constitutive law from the initial predefined material principal coordinate frame to the final curvilinear convected coordinate frame using coordinate transformations is explained, and the theory is implemented into a finite element code.
984
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Fig. 21. Minimum principal stress vector of isotropic canopy.
Fig. 22. Maximum principal stress vector of anisotropic canopy.
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Fig. 23. Minimum principal stress vector of anisotropic canopy.
Fig. 24. Maximum principal stress vector of orthotropic canopy.
985
986
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Fig. 25. Minimum principal stress vector of orthotropic canopy.
Several numerical examples are presented to verify the theory by comparing the dynamic behavior of initially orthotropic/anisotropic membrane structures with mechanically identical braided cable networks. These results show that the new anisotropic constitutive model correctly describes the geometrically nonlinear behavior of anisotropic membranes experiencing large deformations. Finally, an application problem of parachute inflation is analyzed. The stress distributions in the membranes obtained by applying different materials (i.e., isotropic, orthotropic and anisotropic materials) are significantly different.
Acknowledgements The research work described in this work was sponsored by the US Air Force of Scientific Research under Contract number F49620-98-1-02/4 to the University of Connecticut.
Appendix A. When a membrane structure is comprised of braided fibers, its material properties can be considered as a combination of all these fibers’ material properties. As an example, a membrane structure composed of two directional fibers is shown in Fig. 26, and the angles between these two fiber directions and the x-axis are and , respectively. The Young’s modulus of the fibers in the first direction is set toE1 and the second direction as E2 . Two ortho-normal coordinate systems, namely, the x y -coordinate system and
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
987
Fig. 26. General weave pattern of a fabric membrane.
the x y -coordinate system, are defined as shown in Fig. 26. Therefore, the material properties expressed in the x y -coordinate system for the fibers in first direction are E1 0 0 C = 0 0 0 . (35) 0 0 0 After the coordinate transformation, the material properties for the fibers in the first direction can be expressed in the xy-coordinate system as [12] cos2 sin2 2 cos3 sin cos4 C1 = cos2 sin2 (36) sin4 2 cos sin3 E1 . 3 3 2 cos sin 2 cos sin 2 cos2 sin2 Similarly, the material properties for the fibers in the second direction are described in the x y -coordinate system as E2 0 0 C = 0 0 0 (37) 0 0 0
988
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
and in the xy-coordinate system as cos2 sin2 cos4 2 2 C2 = cos sin sin4 3 2 cos sin 2 cos sin3
2 cos3 sin 2 cos sin3 E2 . 2 cos2 sin2
(38)
Therefore, assuming the shear effect is negligible among the fibers, the membrane’s material properties expressed in the xy-coordinate system are C = C1 + C2 .
(39)
References [1] J.J. Nieboer, J. Wismans, P.J.A. De Coo, Airbag modeling techniques, Soc. Automotive Eng. Trans. 99 (6) (1990) 1855–1870. [2] R.J. Benney, J.W. Leonard, A 3-D finite element structural parachute model, Proceedings of the 13th AIAA Aerodynamic Decelerator Conference, Clearwater, FL, 1995. [3] S.D. Prior, A.S. White, Measurements and simulation of a pneumatic muscle actuator for a rehabilitation robot, Simul. Pract. Theory 3 (1995) 81–117. [4] B. Tondu, P. Lopez, Modeling and control of Mckibben artificial muscle robot actuators, IEEE Control Syst. Mag. 20 (2000) 15–38. [5] G. Brown, R. Haggard, R.J. Benney, N. Rosato, A new pneumatic actuator and its use in airdrop applications, AIAA Conference Paper, AIAA-99-1719, Toulouse, France, 1999. [6] N. Stubbs, S. Thomas, A nonlinear elastic constitutive model for coated fabrics, Mech. Mater. 3 (1984) 157–168. [7] F. Sun, A.M. Seyam, B.S. Gupta, A generalized model for predicting load–extension properties of woven fabrics, Textile Res. J. 67 (1997) 866–874. [8] S. Kato, T. Yoshino, H. Minami, Formulation of constitutive equations for fabric membranes based on the concept of fabric lattice model, Eng. Struct. 21 (1999) 691–708. [9] S.K. Kyriacou, C. Schwab, J.D. Humphrey, Finite element analysis of nonlinear orthotropic hyperelastic membranes, Comput. Mech. 18 (1996) 269–278. [10] P.Y. Manach, G. Rio, Analysis of orthotropic behavior in convected coordinate frames, Comput. Mech. 23 (1999) 510–518. [11] M.L. Accorsi, J.W. Leonard, R.J. Benney, K. Stein, Structural modeling of parachute dynamics, AIAA J. 38 (2000) 139–146. [12] R.M. Jones, Mechanics of Composite Materials, Hemisphere Publishing Corporation, Washington, DC, 1975. [13] L.E. Malvern, Introduction to the Mechanics of a Continuous Medium, Prentice-Hall, Inc., Englewood Cliffs, NJ, 1969. [14] A.E. Green, W. Zerna, Theoretical Elasticity, Dover Publication Inc., New York, 1992. [15] K.J. Bathe, Finite Element Procedures, Prentice-Hall, Englewood Cliffs, NJ, 1996. [16] J.W. Leonard, Tension Structures, McGraw-Hill Co., New York, 1988. [17] S.G. Lekhnitski, Theory of Elasticity of an Anisotropic Elastic Body, Holden-Day Publishing Corporation, San Francisco, CA, 1963. [18] V. Firt, Statics. Formfinding and Dynamics of Air-Supported Membrane Structures, Martin Nijhoff Publishers, The Hague, 1983.
Analysis of geometrically nonlinear anisotropic membranes: theory and verification Wenqing Zhang, John L. Leonard, Michael L. Accorsi∗ Department of Civil and Environmental Engineering, University of Connecticut, Storrs, CT 06269-2037, USA Received 30 August 2004; received in revised form 16 November 2004; accepted 29 December 2004
Abstract A new anisotropic constitutive model is introduced to simulate the geometrically nonlinear dynamic behavior of general anisotropic membranes experiencing large deformations. This model relates the second Piola–Kirchhoff stress with the Green–Lagrange strain in a 3-D convected curvilinear coordinate frame. The construction of the anisotropic constitutive law from the predefined principal material coordinate system to the final convected curvilinear coordinate system using coordinate transformations is presented. The proposed theory is implemented into a finite element code and several numerical examples are given for validation. An application problem is described to demonstrate the importance of the new model in real applications. 䉷 2005 Elsevier B.V. All rights reserved. Keywords: Anisotropic/orthotropic constitutive law; Geometrically nonlinear dynamic behavior; Membrane and fabric structures; Finite element method
1. Introduction Membrane and fabric structures are used in a wide variety of applications. Several examples are airbags in the automobile industry [1], parachutes and parafoils in airdrop applications [2], and pneumatic muscle actuators that are used in the robotics industry [3,4] and more recently for parachute control [5]. In many such applications, the membrane structure undergoes large rotations and displacements, but not necessarily large strains. For these cases, a constitutive law that accounts for geometrically nonlinear motion is required. ∗ Corresponding author. Tel.: +1 860 486 5642; fax: +1 860 486 2298.
E-mail address: [email protected] (M.L. Accorsi). 0168-874X/$ - see front matter 䉷 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2004.12.003
964
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Much research has been performed in modeling the nonlinear anisotropic behavior of fabric membranes using a micromechanics approach. For example, Stubbs and Thomas developed a nonlinear elastic constitutive model for coated fabrics in 1984 [6]. In 1997, Sun et al. studied the stress and strain of each individual yarn in a weave [7] and Kato et al. conducted research to model the constitutive properties of weaves from the concept of a fabric lattice [8] in 1999. In these micromechanical fabric models, detailed mechanics of the fabric, such as the effects of coatings, decrimping, friction and slip are modeled. Other researchers have concentrated on macroscopic models for orthotropic materials. Kyriacou et al., for example, presented a finite element code to analyze nonlinear orthotropic and hyper-elastic membranes in large deformation problems in 1996 [9]. The membrane strain energy is expressed in terms of the invariants of the Cauchy–Green strain tensor and the orthotropic properties are pre-defined in the “preferred directions”. In 1999, Manch and Rio established another orthotropic constitutive model in a convected coordinate frame [10] where the strain energy is calculated from the Almansi strain tensor expressed in the convected configuration. In this work, a new macroscopic model is introduced to simulate the geometrically nonlinear dynamic behavior of general anisotropic membranes using a Total Lagrange formulation. A convected curvilinear coordinate frame is used to describe the large deformation and the geometrically nonlinear behavior. General convected anisotropic constitutive laws, which relate the second Piola–Kirchhoff stress to the Green–Lagrange strain, are established in the curvilinear coordinate system on the basis of the material properties predefined in the membranes’ principal material directions. A finite element code has also been developed to implement the anisotropic constitutive matrix for membranes in a general structural dynamics analysis. Several numerical examples are given to validate the capabilities of the model proposed in this research by comparing the dynamic behavior of initially orthotropic/anisotropic membrane structures with mechanically identical woven/braided fiber networks. For woven networks, the fibers are initially orthogonal, whereas for braided networks, the fibers are initially nonorthogonal. Finally, an application of anisotropic materials in parachute design is also presented. In these example problems, a very simple fabric micromechanical model is used to facilitate verification of the results, as discussed later. Specifically, the effects of fiber undulation and decrimping, fabric coating, and friction and slippage are not modeled in these examples. The inclusion of these fabric micromechanical effects into the macroscopic anisotropic model presented here will be considered in future work.
2. Methodology Anisotropic constitutive laws and three-dimensional kinematics to be used in a finite element model [11] are developed in this section from considerations of the generalized Hooke’s law relating stresses and strains [12]. When studying the dynamic behavior of anisotropic membranes, the following two assumptions are made: (1) the membrane material is assumed to be linear elastic and (2) only plane stresses are considered, that is, the stresses in the normal direction of the membranes are assumed to be zeros. In this part, the theoretical kinematic formulations are first developed, followed by the construction of the anisotropic constitutive laws from the principal material coordinate system to the convected curvilinear coordinate system by using coordinate transformation. The incorporation of the anisotropic constitutive theory into an existing finite element code is also presented. In the remaining part of this paper, Latin indices take on the values ranging over 1, 2, 3, while Greek indices take on the values of 1, 2. Repeated indices, unless otherwise noted, imply summation over the
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
965
Fig. 1. Lagrangian description of continuum membrane element.
range of the indices. The arrow above a variable indicates the variable is a vector. The tilde symbol above a variable denotes the variable is a tensor. 2.1. Kinematics A total Lagrangian description [13] is employed in this work to describe the deformation of a continuum solid membrane element. That is, every material particle position is referred to the undeformed configuration. A material particle in the initial undeformed configuration is identified in a fixed global Cartesian coordinate system as = Xi ei X
(1)
and its position in the deformed configuration at time t can be described as x = xi ei ,
(2)
where ei are the unit vectors along the coordinate axis of the global Cartesian coordinate system, Xi the initial coordinate of the particle along the ith axis, and xi the current coordinate of the particle along the ith axis, as shown in Fig. 1. Because of the assumption of plane stress, a 2-D curvilinear coordinate system is defined for every element at time t to describe the membrane element deformation in the convected frame. The covariant base vectors are g =
jx j
=
jxs j
es = gs es ,
(3)
966
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
= jX = jXs es = Gs es . G j
j
(4)
for the initial configuration. The Obviously, g are the base vectors for the current configuration and G covariant components of the metric tensor with respect to the initial and the current configuration are expressed as the following equations, respectively: m · G = Gm G = G G ,
(5)
g = g · g = gm gm ,
(6)
where the dot symbol represents dot product operation. According to Green and Zerna [14], the Green– Lagrange strain tensor is therefore expressed as ⊗ G ˜ = 21 (g − G )G
(7)
with the components = 21 (g − G )
(8)
and ⊗ denotes the dyad of two vectors. The Green–Lagrange strain is defined with respect to the initial configuration, and relates to the second Piola–Kirchhoff stress by the following equation [15]: ˜ , S˜ = C˜
(9)
where S˜ denotes the second Piola–Kirchhoff stress tensor, and C˜ is a fourth-order tensor representing the stress–strain relationships. Once the appropriate expressions of C˜ for the anisotropic material are obtained, the internal forces can be calculated by using Eq. (9), and a stiffness matrix can be also constructed. More detailed explanations are presented in a later section. 2.2. Description of anisotropic material property This section presents the transformation procedure for the anisotropic constitutive laws in the finite element analysis. Three coordinate systems, namely, a principal material coordinate system, a rectangular Cartesian coordinate system and a curvilinear coordinate system, are used to describe the anisotropic material properties of membrane elements. At first, the anisotropic constitutive matrix is defined in the principal material coordinate system, assuming the initial state of the plane sheet of material is flat. Then, a transformation of the constitutive matrix from the principal material coordinate system to the three-dimensional frame of Cartesian coordinate system is performed. After another transformation of constitutive matrix from the Cartesian coordinate system to the curvilinear coordinate system of the element is done, the anisotropic material properties expressed in the element’s curvilinear coordinate system are finally obtained, which are used in the calculations of the finite element analysis later. 2.2.1. Principal material coordinate When studying membranes, it is possible to divide a large membrane into sectors that are arbitrarily oriented in three-dimensional space, with every sector containing one or many isoparametric quadrilateral or triangular membrane elements [16]. These various sectors are assumed to be flat initially, and can be described with anisotropic or orthotropic materials oriented in the principal material directions within
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
967
Fig. 2. Three coordinate systems.
them. Therefore, in this study, each sector is considered as a plane sheet in its initial flat shape and a principal material coordinate system is defined in it. The anisotropic constitutive matrix is then developed in this coordinate system to describe the anisotropic material properties of the elements in the plane sheet sector. 2.2.1.1. Definition of the principal material coordinate system The material properties of a plane sheet are initially described in the principal material coordinate system, assuming that the initial state of the plane sheet is flat. The principal material coordinate system is defined as an orthogonal coordinate system such that two of its axes are mapped to the principal material directions in the plane sheet. As shown in Fig. 2, the unit axial vector m 1 and m 2 are orthogonal vectors in the principal directions on the initial flat sheet. They are determined by the global coordinates of three pre-selected nodes P1 , P2 and P3 on j the sheet, with position vectors R0i . The orthonormal direction cosines i relate m i to the Cartesian unit axis vectors ej by the following equation (no sum on i): j
m i = i ej = di /|di |
(10)
with d1 = R02 − R01 ,
(11.1)
968
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
d3 = d1 × (R03 − R01 ),
(11.2)
d2 = d3 × d1 ,
(11.3)
where the symbol × denotes cross-product operation. j The inverse directional cosines i are defined to relate unit vectors ei to unit vectors m j by j
j, ei = i m j
(12) j
where i are the inverse of i and are determined by ik kj = ij
(13) j
and ij is the Kronecker delta. Because i are orthonormal, Eq. (13) yields j
i = ij .
(14)
2.2.1.2. Description of anisotropic properties in principal direction coordinate system According to Jones [12], the 2-D anisotropic constitutive relations described in the principal material coordinate system can be expressed as
= S t ,
(15)
where and t represent the covariant components of the Green–Lagrange strains and the second Piola–Kirchhoff stresses in the principal material directions, respectively. S are the corresponding stress–strain coefficients with the following symmetry properties:
S = S = S = S = S .
(16)
In terms of physically measurable quantities E1 , E2 , 12 , 21 , G, 1 and 2 , the anisotropic compliances are 11 S11 = 1/E1 ,
(17.1)
22 = 1/E2 , S22
(17.2)
12 12 21 21 = S21 = S12 = S21 = 1/G, S12
(17.3)
11 S22 = − 12 /E1 ,
(17.4)
12 = − 21 /E2 , S21
(17.5)
12 21 11 11 S11 = S11 = S12 = S21 = 1 /E1 ,
(17.6)
12 21 22 22 = S22 = S12 = S21 = 2 /E2 , S22
(17.7)
1 where E1 and E2 stand for the Young’s modulus in the first and second principal material direction (m and m 2 directions), respectively. ij represents the Poisson’s ratio of transverse strain in thej th direction 11 equals to S 12 . G is the shear when stressed in the ith direction. Note that because 12 E2 = 21 E1 , S22 21
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
969
modulus in 1–2 plane. Since the plane studied is assumed in plane stress state, the shear modulus in 1–3 and 2–3 planes are ignored. i is the Lekhnitski coefficient [17] of mutual influence, which is defined as i = 12 /i .
(18)
It characterizes shearing in the plane caused by the normal stress in the ith direction. In some cases where G is difficult to determine, the Firt approximation can be utilized [18]: √
12 / 21 . G = E1 (19) √ 2 1 + 12 21 As a result, the contracted stress–strain relations for the plane sheet in the principal coordinate system can be obtained by inversion of Eq. (6), and be expressed as t = Q .
(20)
Thus, Q11 11 = E1 (1 − G 2 2 /E2 )/ ,
(21.1)
Q22 22 = E2 (1 − G 1 1 /E1 )/ ,
(21.2)
12 21 21 Q12 12 = Q21 = Q12 = Q21 = G(1 − 12 21 )/ ,
(21.3)
22 Q11 22 = Q11 = (E2 12 + G 1 2 )/ ,
(21.4)
21 11 11 Q12 11 = Q11 = Q12 = Q21 = −G( 1 + 12 2 )/ ,
(21.5)
21 22 22 Q12 22 = Q22 = Q12 = Q21 = −G( 2 + 21 1 )/ ,
(21.6)
= (1 − 12 21 ) − G/E1 [2 12 1 2 + 1 1 + (E1 /E2 ) 2 2 ].
(21.7)
where
For the special case of an orthotropic material sheet, 1 and 2 both equal zero. 2.2.2. Cartesian coordinate system Once the stress–strain relations are obtained in the principal material coordinate system, it is necessary to transform them to the Cartesian coordinate system. The Cartesian coordinate system employed here is identical to the fixed global Cartesian coordinate system mentioned before and has the unit axis vector ei . The 3-D relations between stresses and strains in this system can be expressed as ij
Trs = Crs Eij ,
(22)
where Trs represent the second Piola–Kirchhoff stress, Eij the Green–Lagrange strain components, both in the reference Cartesian coordinate frame. Considering Eq. (22) together with Eqs. (9) and (20), the components of the constitutive matrix expressed in the Cartesian coordinate system are given by ij
j
p
t Crs = im n Qmn pt r s .
(23)
970
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Because of the plane-stress state assumed, only the Q terms are non-zeros. Thus, Eq. (23) can be rewritten as ij
j
Crs = i Q r s .
(24)
2.2.3. Curvilinear coordinate system In order to obtain the final representation of constitutive laws in the current curvilinear coordinate system, another transformation is needed. The convected curvilinear coordinate system described previously based on every element integration point is also employed here at time t with the natural base vectors as described in Eq. (3). The contravariant base vectors are g = hs es ,
(25)
whose components hs are determined by g · g = hs gs =
(26)
and the contravariant components of the metric tensor are h = g · g = hm hm .
(27)
Finally, the constitutive relations in the element curvilinear system are expressed as S = E ,
(28)
where S denote the contravariant components of the second Piola–Kirchhoff stress in the element curvilinear coordinate system, the covariant components of the Green–Lagrange strain components. By comparing Eq. (28) with Eqs. (22) and (9), we see that the components of the constitutive tensor in the “integration-point-wise” curvilinear coordinate system can be calculated by
E = E h h
(29)
with
mn p t E = h m h n Cpt g g .
(30)
This completes the transformation of constitutive components from Q given in Eq. (20) to E in the curvilinear coordinate system at time t. The results obtained here are implemented in the finite element code to perform the anisotropic material analysis. 2.3. Finite element implementation The anisotropic constitutive laws derived previously have been incorporated into an existing finite element code to execute the nonlinear dynamic analysis of anisotropic membranes. A FORTRAN program named TENSION [11], which is developed for the 3-D nonlinear dynamic analysis of membrane, cable and payload structures, is used. In this section, the procedures of the implementation are described. As mentioned previously, once the current curvilinear coordinate system is determined, the Green– Lagrange strains are calculated from Eqs. (7) and (8). Note that the components gs of the covariant base
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
971
vectors are obtained from the derivatives of the element shape functions [16] in the program. Then, the components of the second Piola–Kirchhoff stress S are determined from Eq. (28) and are used in the calculation of the internal restoring force RN i and tangent stiffness matrix KMN ij [11]. RN i =
1 2
√ S (HM, HN, + HN, HM, )yMi A dX1 dX2
(31)
KMN ij =
1 j 1 i S (HM, HN, + HM, HN, ) + E (HR, HN, + HR, HN, )yRi 2 4 √ ×(HQ, HM, + HM, HQ, )yQj A dX1 dX2 , (32)
where RN i denotes the ith component of the internal force on node N, yMi the current coordinate of node M in the ith direction. KMN ij are the components of the tangent stiffness matrix, H represents the shape functions and comma in the subscript denotes the partial derivatives of the variable with respect to the subscript after the comma. These equations express the current internal force and stiffness at time t, but referred to the initial undeformed state. A more detailed study of the dynamic analysis method used in TENSION can be found in [11]. Therefore, a brief flow chart of the program using an implicit time integration scheme and full Newton–Raphson iteration with the geometrically nonlinear anisotropic constitutive laws can be expressed as follows: (i) Set up the model, input the anisotropic material properties in principal material coordinate system (ii) Transform the anisotropic material properties from the principal material coordinate system to Cartesian coordinate system (Transformation 1) (iii) Begin the time step loop (iv) Begin Newton–Raphson iteration (v) Element loop (vi) Calculate the curvilinear coordinate system and Green–Lagrange strain (vii) Transform the anisotropic material properties from Cartesian coordinate system to curvilinear coordinate system (Transformation 2) (viii) Calculate the second Piola–Kirchhoff stress, restoring forces and tangent stiffness matrix (ix) End element loop (x) Construct and solve system equations (xi) Convergence check, if not convergent, go to step (iv) (xii) End time step loop ij
The first transformation takes place only once and the calculation of Crs is the same for all the elements in a particular sheet of material. The second transformation takes place for every membrane element in all the iterations of every time step, because the deformation changes and therefore the curvilinear coordinate system of each element changes in time. Since the constitutive relation is evaluated within the Newton–Raphson iteration, additional nonlinear effects (such as material nonlinearities) could also be included in the current formulation.
972
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
3. Numerical examples Two numerical examples are presented to verify the constitutive model developed previously. In each example, the dynamic behavior of a membrane structure is compared to a simple woven/braided cable network that is mechanically identical to the membrane structure. The model is initially orthotropic in the first example and initially anisotropic in the second example. 3.1. Comparison of an initially orthotropic membrane structure with a mechanically identical cable network Consider the woven cable network shown in Fig. 3. The width and height of the network are w and h, respectively. The number of rows of horizontal cables and columns of vertical cables are both n. The vertical cables are woven with the horizontal cables according to the pattern shown in Fig. 4. A simple fabric micromechanical model is assumed for the cable network to facilitate comparison with the macroscopic anisotropic membrane model. It is assumed that the horizontal and vertical cables are flat (i.e. no undulation). It is further assumed that the horizontal and vertical cables are pinned at their junctions and that no rotational resistance or slippage occurs at the junctions. Under these assumptions, the properties of the mechanically equivalent anisotropic membrane can be determined analytically. If the cross-section areas of the vertical cables and the horizontal cables are (w/n)t and (h/n)t, respectively, then the cable network is equivalent to a membrane structure with the same dimensions, w
Fig. 3. A woven cable network.
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
973
Fig. 4. Initial weave pattern of the cable network.
Fig. 5. An initially orthotropic cable network and membrane model.
and h, and thickness t. The Young’s modulus of the horizontal cables and vertical cables are set as E1 and E2 , respectively. Assuming that the density of all cables is and there is no shear between cables, then the membrane structure, which has the same mass and material properties as the cable network, should
974
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Fig. 6. Deformed shape of the initially orthotropic membrane structure.
have a density of 2 and the following stress–strain relations in the Cartesian coordinate system: E1 0 0 C= 0 E2 0 . 0 0 0
(33)
Obviously, this membrane structure is initially orthotropic. The methodology of how to derive the corresponding material property for a membrane structure to be identical to a braided cable network is given in Appendix A. If the boundary conditions of the membrane structure and the cable network are set the same, these two models are then mechanically identical and they should present the same deformation when subjected to the same loads. In this example, the left edge and the bottom edge of both models are fixed, and a point force P with a diagonal magnitude of 0.4 pound is applied to the node on the upper right corner, as shown in Fig. 5. The width w and height h for these two models are both 1 in. In the membrane structure, the thickness t is 0.4 in and the density is 0.8 pound/in3 . In the cable network, n is equal to 50; the density is 0.4 pound/in3 ; the cross-section area for all the cables is 0.005 in2 ; and E1 and E2 are both 106 pound/in2 . The deformation shapes for both models are tracked to verify the theory. The change of the angles between the horizontal
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
975
Fig. 7. Deformed shape of the initially orthotropic cable network.
cable and the vertical cable at two randomly selected points in the cable network, points A and B, as in Fig. 5, is also compared with the corresponding change of the angles obtained from the membrane structure at the same points. Figs. 6 and 7 show the deformation shapes of the membrane structure and the cable network after the point force is applied for 0.4 s, respectively. The deformation shapes of these two models agree very well. Fig. 8 shows the angle changes obtained from both the membrane structure and the braided cable network at points A and B. The results demonstrate that, over a large range of deformation, the initially orthotropic membrane structure behaves the same as the mechanically identical cable network. 3.2. Comparison of an initially anisotropic membrane structure with a mechanically identical cable network In this example, the dynamic behavior of a membrane structure is compared with a mechanically identical braided cable network. Fig. 9 shows the initial weave pattern of the cable network in this example. The same simplifying assumptions used in the previous example for the fabric micromechanical behavior are also used here to facilitate comparison with the membrane model. As a result, a mechanically identical membrane structure not only has the same dimensions and mass as the cable network, but also has the
976
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Fig. 8. Angle change versus time for the initially orthotropic membrane structure and cable network.
Fig. 9. Initial weave pattern of the anisotropic cable network.
stress–strain relations in the Cartesian coordinate system as
0.0045(E1 + E2 ) C = 0.0625(E1 + E2 ) 0.0335(E1 − E2 )
0.0625(E1 + E2 ) 0.8705(E1 + E2 ) 0.4665(E1 − E2 )
0.0335(E1 − E2 ) 0.4665(E1 − E2 ) . 0.1250(E1 + E2 )
(34)
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
977
Fig. 10. An initially anisotropic membrane structure and cable network model.
Fig. 11. Deformed shape of the initially anisotropic membrane at 0.26 s.
Assuming the Young’s modulus for cables in directions 1 and 2, which are marked in Fig. 9, are E1 and E2 , respectively. More detailed derivation for this constitutive matrix can be found in Appendix A. Eq. (34) indicates that the membrane structure initially possesses anisotropic material properties. Fig. 10 shows the initial shapes and the boundary conditions of the membrane structure and the cable network model used in this example. As can be seen, the left edges of both models are fixed, and the right edges are subject to a uniform edge force F with a magnitude of 8.1 lb/in. One point, point C, is also randomly selected to compare the dynamic behavior of the membrane structure and the braided cable network. During the analysis, the length of the horizontal edge w is 1 in and the height h is 0.5 in for both models. In the membrane structure, the thickness t is 0.1 in and the density is 0.8 pound/in3 . In the
978
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Fig. 12. Deformed shape of the initially anisotropic cable network at 0.26 s.
Fig. 13. X Displacement comparison of Point C.
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Fig. 14. Y Displacement comparison of Point C.
Fig. 15. Angle change at Point C.
979
980
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Fig. 16. Initial mesh of a quarter-scale C9 parachute model.
cable network, n is equal to 80; the density is 0.4 pound/in3 ; the cross-section area for all the cables is 0.00125 in2 and E1 and E2 are both 100 pound/in2 . Figs. 11 and 12 show the deformed shapes at a time of 0.26 s for the membrane structure and the cable network, respectively. The displacements of point C obtained from the two models are given in Figs. 13 and 14, and the change of the angles between intersecting cables/fibers at point C is also given in Fig. 15. All these results indicate that the two models demonstrate the same dynamic behavior.
4. Applications The anisotropic theory developed here has been applied to the simulation of dynamic parachute problem. In this example, a quarter-symmetric model of a round C9 parachute is used and the inflation of the parachute from the initially unstressed state to the totally inflated state is simulated. The finite element mesh of the initial configuration is shown in Fig. 16. This model consists 88 three-node cable elements
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
981
Fig. 17. Fully inflated shape of a quarter-scale C9 parachute model.
for the cables in canopy, 32 two-node cable elements the suspension lines, 1 payload element and 154 nine-node membrane elements for the canopy, which is made up of 7 gores. Three kinds of materials, i.e., isotropic, orthotropic and anisotropic materials, are used separately as the canopy fabrics, and the stress distributions along the centerline of one gore are compared. For all cases, the true (Cauchy) stress is presented. In this example, the isotropic material used has the Young’s modulus as 4.32 × 105 lb/ft2 , Poisson’s ratio as 0.3; the orthotropic material used has the Young’s modulus of 4.32 × 105 lb/ft 2 in one direction and 2.16 × 105 lb/ft2 in another direction; the anisotropic material is set to be weaved up by two directional fibers with the Young’s modulus as 4.32 × 105 lb/ft2 , and the intersecting angle is 60◦ . The parachute canopy undergoes 0.5 lb/ft 2 pressure. The canopy’s thickness is 0.1 ft and the payload is 50 lb for all cases. The membrane thickness is 0.0001 ft and the cable cross-section area is 0.0001 ft 2 . Fig. 17 represents the inflated configuration of the model at 6.0 s. Although three different materials are used for the canopy fabrics, the deformation shape of the parachute does not change much. However, the maximum and the minimum principal stress distributions differ significantly. Figs. 18 and 19 give the magnitudes of the maximum principal stress and minimum principal stress along the centerline of a canopy gore. Figs. 20–25 show the top view of distributions for the maximum and minimum principal
982
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Fig. 18. Maximum principal stress for different materials.
Fig. 19. Minimum principal stress for different materials.
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
983
Fig. 20. Maximum principal stress vector of isotropic canopy.
stress vectors, respectively. As can be seen, when anisotropic membrane is used, the maximum principal stresses on the canopy are all in circumferential direction and the minimum principal stresses are all meridian stresses. When isotropic material or orthotropic material is used, the maximum principal stresses are circumferential stress around the canopy center and change to meridian stresses with the increase of radius, while the minimum principal stresses are in radial direction around the canopy center and change to hoop stresses with radius increase. Note that the radius where the principal stresses begin to change direction is different for isotropic canopy and orthotropic canopy. Figs. 18 and 19 show that the magnitudes of principal stresses are different when different materials are used.
5. Conclusions In this research, the theoretical formulation, finite element implementation and numerical validation of a new anisotropic constitutive law are presented to describe the geometrically nonlinear dynamic behavior of general anisotropic membranes with large deformations. The theory is based on a Total Lagrange formulation using the second Piola–Kirchhoff stress and the Green–Lagrange strain. A 3-D convected curvilinear coordinate frame is employed to describe the large deformation and the geometrically nonlinear behavior of membranes. The methodology of how to derive the constitutive law from the initial predefined material principal coordinate frame to the final curvilinear convected coordinate frame using coordinate transformations is explained, and the theory is implemented into a finite element code.
984
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Fig. 21. Minimum principal stress vector of isotropic canopy.
Fig. 22. Maximum principal stress vector of anisotropic canopy.
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Fig. 23. Minimum principal stress vector of anisotropic canopy.
Fig. 24. Maximum principal stress vector of orthotropic canopy.
985
986
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
Fig. 25. Minimum principal stress vector of orthotropic canopy.
Several numerical examples are presented to verify the theory by comparing the dynamic behavior of initially orthotropic/anisotropic membrane structures with mechanically identical braided cable networks. These results show that the new anisotropic constitutive model correctly describes the geometrically nonlinear behavior of anisotropic membranes experiencing large deformations. Finally, an application problem of parachute inflation is analyzed. The stress distributions in the membranes obtained by applying different materials (i.e., isotropic, orthotropic and anisotropic materials) are significantly different.
Acknowledgements The research work described in this work was sponsored by the US Air Force of Scientific Research under Contract number F49620-98-1-02/4 to the University of Connecticut.
Appendix A. When a membrane structure is comprised of braided fibers, its material properties can be considered as a combination of all these fibers’ material properties. As an example, a membrane structure composed of two directional fibers is shown in Fig. 26, and the angles between these two fiber directions and the x-axis are and , respectively. The Young’s modulus of the fibers in the first direction is set toE1 and the second direction as E2 . Two ortho-normal coordinate systems, namely, the x y -coordinate system and
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
987
Fig. 26. General weave pattern of a fabric membrane.
the x y -coordinate system, are defined as shown in Fig. 26. Therefore, the material properties expressed in the x y -coordinate system for the fibers in first direction are E1 0 0 C = 0 0 0 . (35) 0 0 0 After the coordinate transformation, the material properties for the fibers in the first direction can be expressed in the xy-coordinate system as [12] cos2 sin2 2 cos3 sin cos4 C1 = cos2 sin2 (36) sin4 2 cos sin3 E1 . 3 3 2 cos sin 2 cos sin 2 cos2 sin2 Similarly, the material properties for the fibers in the second direction are described in the x y -coordinate system as E2 0 0 C = 0 0 0 (37) 0 0 0
988
W. Zhang et al. / Finite Elements in Analysis and Design 41 (2005) 963 – 988
and in the xy-coordinate system as cos2 sin2 cos4 2 2 C2 = cos sin sin4 3 2 cos sin 2 cos sin3
2 cos3 sin 2 cos sin3 E2 . 2 cos2 sin2
(38)
Therefore, assuming the shear effect is negligible among the fibers, the membrane’s material properties expressed in the xy-coordinate system are C = C1 + C2 .
(39)
References [1] J.J. Nieboer, J. Wismans, P.J.A. De Coo, Airbag modeling techniques, Soc. Automotive Eng. Trans. 99 (6) (1990) 1855–1870. [2] R.J. Benney, J.W. Leonard, A 3-D finite element structural parachute model, Proceedings of the 13th AIAA Aerodynamic Decelerator Conference, Clearwater, FL, 1995. [3] S.D. Prior, A.S. White, Measurements and simulation of a pneumatic muscle actuator for a rehabilitation robot, Simul. Pract. Theory 3 (1995) 81–117. [4] B. Tondu, P. Lopez, Modeling and control of Mckibben artificial muscle robot actuators, IEEE Control Syst. Mag. 20 (2000) 15–38. [5] G. Brown, R. Haggard, R.J. Benney, N. Rosato, A new pneumatic actuator and its use in airdrop applications, AIAA Conference Paper, AIAA-99-1719, Toulouse, France, 1999. [6] N. Stubbs, S. Thomas, A nonlinear elastic constitutive model for coated fabrics, Mech. Mater. 3 (1984) 157–168. [7] F. Sun, A.M. Seyam, B.S. Gupta, A generalized model for predicting load–extension properties of woven fabrics, Textile Res. J. 67 (1997) 866–874. [8] S. Kato, T. Yoshino, H. Minami, Formulation of constitutive equations for fabric membranes based on the concept of fabric lattice model, Eng. Struct. 21 (1999) 691–708. [9] S.K. Kyriacou, C. Schwab, J.D. Humphrey, Finite element analysis of nonlinear orthotropic hyperelastic membranes, Comput. Mech. 18 (1996) 269–278. [10] P.Y. Manach, G. Rio, Analysis of orthotropic behavior in convected coordinate frames, Comput. Mech. 23 (1999) 510–518. [11] M.L. Accorsi, J.W. Leonard, R.J. Benney, K. Stein, Structural modeling of parachute dynamics, AIAA J. 38 (2000) 139–146. [12] R.M. Jones, Mechanics of Composite Materials, Hemisphere Publishing Corporation, Washington, DC, 1975. [13] L.E. Malvern, Introduction to the Mechanics of a Continuous Medium, Prentice-Hall, Inc., Englewood Cliffs, NJ, 1969. [14] A.E. Green, W. Zerna, Theoretical Elasticity, Dover Publication Inc., New York, 1992. [15] K.J. Bathe, Finite Element Procedures, Prentice-Hall, Englewood Cliffs, NJ, 1996. [16] J.W. Leonard, Tension Structures, McGraw-Hill Co., New York, 1988. [17] S.G. Lekhnitski, Theory of Elasticity of an Anisotropic Elastic Body, Holden-Day Publishing Corporation, San Francisco, CA, 1963. [18] V. Firt, Statics. Formfinding and Dynamics of Air-Supported Membrane Structures, Martin Nijhoff Publishers, The Hague, 1983.