- Email: [email protected]
A formulation of elasticity and viscoelasticity for fibre reinforced material at small and finite strains
Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
www.elsevier.com/locate/cma
A formulation of elasticity and viscoelasticity for ®bre reinforced material at small and ®nite strains M. Kaliske Institut f ur Statik, Universit at Hannover, Appelstraûe 9a, 30167 Hannover, Germany Accepted 26 March 1999
Abstract The combination of unidirectional ®bres embedded in matrix yields a novel material which is heterogeneous on the micro level and exhibits new mean characteristics on the macro level. Due to the fact that the length scale of a composite structure is several orders larger compared to the micro scale, the composite material is assumed to be homogeneous on the structural level. Therefore, a homogeneous anisotropic constitutive formulation has to be introduced in order to represent the material and the composite structure eciently, for example by a ®nite element discretization. In case of unidirectional reinforced composites a transversely isotropic idealization is meaningful. One typical class of composites is manufactured, for example, from glass or carbon ®bres combined with an epoxy matrix. Usually, structures produced from these materials are subjected to small strains. A geometrically nonlinear description is essential when modelling steel or nylon cord reinforced elastomeric material. The investigation of biomaterials like muscles or bones, which also have a ®brous microstructure, is accomplished on the basis of ®nite strains, too. The static and the dynamic response of a structure can be simulated with the ®nite element method. The stiness is characterized by anisotropic elastic properties. For the dynamic response, the dissipative behaviour, originated in the viscoelastic properties of the polymeric material, is of great importance. Depending on the class of composites investigated, the elastic and the viscoelastic behaviour of ®bre reinforced material has to be formulated by a transversely isotropic approach at small and at ®nite strains respectively. With these models at hand, realistic ®nite element simulations may be carried out. Ó 2000 Elsevier Science S.A. All rights reserved. Keywords: Elasticity; Viscoelasticity; Fibre reinforced material
1. Introduction The status of composite material in engineering application is increasing steadily. Fibre matrix combinations are employed in the aircaft industry and for the construction of space structures. Moreover, composites are also used for civil engineering structures like footbridges. This class of `sti composites', made for example from glass ®bres in epoxy resin is enlarged by so-called `soft composites'. This second type of material is formed by cords with high tensile strength reinforcing an elastomer. Typical areas of application are tires, conveyer belts or pneumatic shock absorbers. The physical description of the industrial products mentioned has to be extended to biomaterials which often exhibit an anisotropic microstructure. Therefore, constitutive models for bones or soft living tissues belong to the same class as the arti®cially produced material. The manufacturing of composites allows the creation of speci®c characteristics of the combined material through variation of geometrical parameters like ®bre volume fraction, ®bre angle and even ®bre shape or through the combination of dierent components. This way, a material design is accomplished to optimize structural characteristics. 0045-7825/00/$ - see front matter Ó 2000 Elsevier Science S.A. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 9 9 ) 0 0 2 6 1 - 3
226
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
A numerical tool to investigate structures eciently is the ®nite element method. By this approach, it is possible to solve complex initial boundary value problems. In order to compute realistic results a constitutive formulation has to be employed that incorporates all main features of the material. Principle properties to be taken into account are anisotropic elasticity and anisotropic viscoelasticity. The description of linear elastic behaviour of ®bre reinforced material is shown in standard text books on composites. Assuming linear elasticity and small strains Jones [7], Christensen [1] or Gibson [4] for example present the constitutive relations for dierent types of anisotropy based on engineering constants. For the geometrically nonlinear regime the linear elastic approach is often used as to be seen in the literature, i.e., modelling is accomplished by the so-called de St. Venant Kirchho material. Klee and Stein [13] show that the formulation may be employed just for a small range of deformation, i.e., up to 5% of strain. General considerations on anisotropic elasticity at ®nite strains are given in Truesdell and Noll [28], Suhubi [25] and Lurie [15]. Spencer [24] developed a constitutive strain energy function in terms of invariants. For this purpose, he extended the basis of three characteristic quantities of isotropy to ®ve invariants of transverse isotropy. Schr oder [18] and Weiss et al. [30] investigate ®nite and anisotropic nonlinear elasticity. Moreover, Weiss et al. [30] show a speci®c strain energy function applied to ®nite element simulations of biological material. Time- and frequency-dependent behaviour is observed not only for an isotropic polymer like epoxy resin or rubber, but the composite itself is also viscoelastic. Biomaterials, which are mostly anisotropic, exhibit viscoelastic behaviour as well (see for example [23,26]). Dierent theoretical models and numerical realizations for the formulation of anisotropic linear viscoelastic properties of thermo-rheologically simple materials (see [19]) at small strains have been published. The approach shown below is motivated by its applicability to a linear as well as to a geometrically nonlinear description. An overview of the present literature on modelling of anisotropic viscoelastic material reponse reveals that most of the authors follow similar ideas. The article of Lin and Hwang [14] is one of the ®rst publications dealing with the mentioned topic in the context of the ®nite element method. They introduce a thermo-viscoelastic model for plane strain and its algorithmic formulation based on the recurrence relation of Taylor et al. [27]. Yi [33] computes the solution of the linear quasi-static anisotropic problem in the Laplace- and in the Fourier-domain. A time-dependent simulation is accomplished by Yi according to the algorithm of Taylor et al. [27] employing a symbolic generalized Maxwell-element. In a number of papers ([33±35]) the authors start with a viscoelastic variation principle, discretize it with the ®nite element method and ®nally end up with integral equations, which can be solved recursively. Wang [29] and Kennedy and Wang [12] developed a solution procedure for nonlinear integral equations utilizing the viscoelastic approach of Lou and Schapery [10] in order to describe the long-term nonlinear viscoelastic response of composite material. The formulation is applied to ®nite element computations. Zocher et al. [36] present an algorithmic approach to orthotropic linear thermo-viscoelasticity. The integral equations of the material history are developed in an incremental form and solved recursively similar to the method of Taylor et al. [27]. In a recent publication Poon and Ahmad [17] also introduce an eective and robust integration algorithm for anisotropic linear viscoelaticity which can considered to be similar to the model of Zocher et al. [36]. An extension of the small strain approach to large strain viscoelasticity is provided in the articles of Simo [21], Govindjee and Simo [5] and Kaliske and Rothert [8] for isotropic material. Basics of ®nite strain viscoelastic formulations are developed by Sidoro [20] and Lubliner [11]. Le Tallec et al. [26] derived a ®nite and nonlinear viscoelastic model and used the approach for ®nite element simulations of the human eye. The outline of this paper is the following: In Section 2 transversely isotropic formulations for elastic and viscoelastic characteristics at small strains are given. Subsequently, Section 3 extends the general structure of the approaches to the ®nite strain regime. Numerical examples of Section 4 illustrate the application of the models to ®nite element computations of a dynamic structural response, reinforced rubber material and simulations of biological tissues. Main features are summarized in Section 5. Finally, Appendix A provides some useful results of tensor analysis, which are employed in the following considerations.
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
227
2. Linear theory 2.1. Elastic constitutive model Hyperelastic material is de®ned by the strain energy function W. In the case of a geometrically linear and isotropic theory, it is a function of the small strain tensor e. A unidirectional ®bre reinforced composite is characterized by a unit vector in ®bre direction a, where kak 1. Transverse isotropy is obtained for ®bres in a hexagonal order, i.e., material properties are invariant with respect to a rotation about the ®bre axis a. Due to this second constitutive characteristic a, the strain energy function W of a composite is formulated in terms of the strain tensor e and the dyadic product a a, which is a second order tensor of the microstructure, i.e., the ®bre orientation.The invariance requirement of the strain energy W e; a a W QeQT ; Qa aQT
1
has to be ful®lled, if e and a are subjected to a rotation by a proper orthogonal tensor Q e0 QeQT
a0 Qa;
and
2
where QQT QT Q 1 with the unit tensor 1 of second order and det Q 1. Spencer [24] derived the basis tr e; tr e2 ; tr e3 ; a : e : a; a : e2 : a
3
of ®ve irreducible invariants of transverse isotropy. Subsequently, the most general quadratic strain energy function k b W tr2 e lT tr e2 a aeatr e 2 lL ÿ lT ae2 a aea2 2 2
4
is formulated using the ®ve independent material parameters k, lT , a, lL , b. Standard evaluations of continuous mechanics and tensor analysis lead to the accompanying stress tensor r
oW ktr e a aea1 2lT e atr e b aea a a 2 lL ÿ lT a ae ea a oe
5
as derivative of the potential function W with respect to the strain tensor e. The linear elastic material tensor a
o2 W k 1 1 a a a 1 1 a a 2lT I b a a a a 2 lL ÿ lT Ia oe2
6
is computed as second derivative of the potential function with respect to e, where 1 is a second order and I is a fourth order unit tensor. The quantity Ia stands for a fourth order tensor formed by the components of the unit vector a (see Eq. (A. 13) in Appendix A for an explicit representation). Supposing a ®bre orientation in ~x1 direction, i.e., aT 1; 0; 0, the constant material tensor (6) reduces to 2 6 6 6 6 6 6 6 ~a 6 6 6 6 6 6 6 6 4
k 2a 4lL ÿ 2lT b
ka
ka
0
0
ka
k 2lT
k
0
0
ka
k
k 2lT
0
0
0
0
0
2lL
0
0
0
0
0
2lL
0
0
0
0
0
0
3
7 7 0 7 7 7 7 0 7 7 7 7 0 7 7 7 0 7 7 5 2lT
7
228
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
in matrix form. The inverse ~ aÿ1 of the 2 1 m12 m12 ÿ ÿ 6 E11 E E 11 11 6 6 m12 1 m23 6ÿ ÿ 6 E11 E22 E22 6 6 m12 m23 1 6ÿ 6 E11 ÿ E22 E22 ~aÿ1 6 6 6 0 0 0 6 6 6 6 0 0 0 6 6 4 0 0 0
classical formulation employing engineering constants yields 3 0 0 0 7 7 7 0 0 0 7 7 7 7 0 0 0 7 7 7 7 1 0 0 7 7 G12 7 7 1 0 7 0 7 G12 7 1 5 0 0 G23
8
for a transversely isotropic composite in material coordinates ~xi (see for example [7]). The shear modulus G23 in Eq. (8) is not an independent parameter and may be computed according to G23
E22 2 1 ÿ m23
9
from Young's modulus E22 and Poisson's ratio m23 . The strain energy function of the classical approach 1 W e:a:e 2
10
is given directly by the strain e and the material tensor a. Since the shown linear formulations in Eqs. (7) and (8) both describe transversely isotropic material, the engineering constants (Eq. (8)) and the coecients of the invariant based approach (Eq. (7)) are equivalent (see [18] for the conversion). 2.2. Viscoelastic constitutive model The so-called generalized Maxwell-model has proved to be a useful rheological structure to formulate time-dependent and frequency-dependent characteristics. The symbolic material model is composed of an elastic spring in parallel with n separate Maxwell-elements j (Fig. 1). This approach results in a stress tensor n X hj 11 r r0 j1
formed additively from purely elastic stresses r0 a0 e and internal viscoelastic stress variables hj , which are determined by the linear rate-equation 1 h_ j hj cj r_ 0 sj
12
for a symbolic material element j. Eq. (12) de®nes the viscoelastic stress variables hj as function of the stress rate r_ 0 . Material parameters are the relaxation time sj , which is chosen to be independent of the ®bre
Fig. 1. Generalized Maxwell-element.
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
229
orientation and, consequently, results in a scalar property. The model is charaterized by the fourth order relaxation tensor cj . This quantity represents the anisotropic viscosity of the material. A similar isotropic formulation employing a scalar factor cj , is discussed by Simo [21] and Kaliske and Rothert [8]. The relaxation time (Eqs. (12) and (46)) used in the mathematical model do not necessarily have to be identical for the different components of the material tensor. In this respect Yi and Hilton [35], Zocher et al. [36] and Poon and Ahmad [17] consider the relaxation time of the anisotropic material as a direction-controlled quantity. The transversely isotropic relaxation tensor for one Maxwell-element is given in material coordinates ~xi by the matrix 2
c11
6 6 c21 6 6 6 c21 ~cj 6 6 6 0 6 6 4 0 0
c12
c12
0
0
c22
c23
0
0
c23
c22
0
0
0
0
c44
0
0
0
0
c44
0
0
0
0
0
3
7 0 7 7 7 0 7 7: 7 0 7 7 7 0 5
13
c66
In contrast to the symmetric elasticity tensor a, the relaxation tensor is unsymmetric for the components c12 and c21 . Due to the speci®c symmetries of transverse isotropy c23 c32 is found. In case of isotropy, the approach reduces ~cj cj I
14
to a scalar factor ci and the fourth order identity tensor I. Eq. (12) carries over to the well known isotropic rate-equation. The fourth order relaxation tensor (13) is computed in a formal step according to ~cj ~aj a~ÿ1 0 ;
15
analogously to the isotropic model, where cj lj =l0 de®nes the contribution of one Maxwell-element j to the total stiffness. The scalar parameters l0 and lj are the elastic constants of the symbolic spring and the Maxwell-element j. Equivalently, ~ a0 in Eq. (15) stands for the purely elastic behaviour of the ®bre reinforced material and ~ aj is the transversely isotropic elasticity tensor of Maxwell-element j. The total timedependent stiffness of the generalized Maxwell-element is given by ( ) n X t cj exp ÿ 16 a0 c ta0 ; a t I sj j1 where c t is the tensorial relaxation function. The relaxation tensor ~cj is transformed either by using the ®bre direction a or by T ÿT ÿT aj T Tÿ1 ~ aÿ1 aÿ1 TT ~ aj ~ TT~cj TÿT cj aj aÿ1 0 T ~ 0 T 0 T
from local material coordinates to a global reference system (see Eq. (A.14) in Appendix A). Starting from the linear rate-equation (12) the hereditary integral Z t t ÿ t or0 cj exp ÿ dt hj sj ot 0
17
18
is derived. It describes the history of the internal stress variables hj . The integration with respect to time is accomplished eciently by a recurrence relation hjn1
exp
ÿ 1 ÿ exp ÿ Dt=sj n1 Dt n r0 ÿ rn0 ÿ hj c j sj Dt=sj
19
230
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
for time step Dt tn1 ÿ tn . This approach assumes a linear stress rate r_ 0
dr0 Dr0 rn1 ÿ rn0 lim lim 0 Dt!0 Dt Dt!0 dt Dt
20
during time step Dt. The recurrence formula is obtained exploiting the features of the exponential function as shown by Herrmann and Peterson [6] and Taylor et al. [27] in a similar way. For the evaluation of the recurrence relation (19) the known solutions hnj , rn0 at time tn are needed. This information has to be stored in a data base for one time step, i.e., for the evaluation at time tn1 .The algorithmic material tensor is crucial in the context of a ®nite element implementation. The fourth order tensor ( ÿ ) n X 1 ÿ exp ÿ Dt=sj orn1 n1 cj 21 an1 a : n1 I 0 oe Dt=s j j1 is determined as derivative of the algorithmic stresses rn1 with respect to strain en1 . Tensor (21) is symmetric. It shows an interesting structure because the viscous part fg and the elastic characteristics a0n1 are fully separated. This formal split is based on the introduction of the relaxation tensor cj . The advantage of the proposed formulation will be obvious for the extension to ®nite strains. The nonlinearities of a ®nite anisotropic linear viscoelasticity are associated with the nonlinear elastic behaviour and with the introduction of a ®nite strain measure. Viscoelasticity is still modelled by a linear rate-equation. Therefore, the separation of linear viscous and nonlinear elastic contributions is useful. 3. Finite strain formulation 3.1. Elastic constitutive model The invariant based elastic approach in Section 2.1 can be extended to a nonlinear transversely isotropic formulation at ®nite strains. On the other hand, an application of the classical linear description (Eq. (8)) to geometrically nonlinear considerations using a constant material tensor A yields a strain energy function 1 22 W E : A : E; 2 where E is the Green Lagrangian strain tensor. This physically linear so-called de St. Venant Kirchho material is valid only in a small range of strains. The strain energy function W of a nonlinear elastic ®bre reinforced material with respect to the reference con®guration is a function of the Right Cauchy Green tensor C and the ®bre direction A, i.e., the tensor A A of the microstructure of the material, where kAk 1. A push forward operation of the unit vector A to the current con®guration a U A FA
23
is accomplished by the deformation gradient F. The strain energy function for ®nite strains analogously to Eq. (1) is given by W W C; A A W QCQT ; QA AQT
24
postulating invariance with respect to a rotation Q. Spencer [24] derives the basis of irreducible invariants 1 25 I tr C; II tr2 C ÿ tr C2 ; III det C; IV A : C : A; V A : C2 : A 2 equivalently to the small strain formulation. The three invariants I, II, III of isotropy are extended by IV and V for transverse isotropy. Subsequently, the elastic potential W W I; II; III; IV; V
26
is formulated as function of the ®ve invariants.The second Piola Kirchho stress tensor S relative to the reference con®guration
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
S2
oW oC
231
27
and the accompanying fourth order material tensor A4
o2 W oC2
28
are computed according to standard methods (see [2]). The general form of the second Piola Kircho stress tensor for the elastic potential as given by Eq. (26) is derived S 2 WI IWII 1 ÿ 2 WII C 2III WIII Cÿ1 2 WIV A A 2 WV M
29
M : A CA CA A:
30
using The tensorial quantities are multiplied by scalar factors consisting of derivatives of the potential function with respect to the invariants WI : oW=oI; WII : oW=oII; WIII : oW=oIII etc. and by invariants I, III. The entire structure of the general material tensor A 4 WI;I WII 2IWI;II I2 WII;II 1 1 ÿ 4 WI;II ÿ 2IWII;II 1 C C 1 4 WII;II C C
ÿ 4 IIIWIII III2 WIII;III Cÿ1 Cÿ1 ÿ 4 IIIWI;III I IIIWII;III 1 Cÿ1 Cÿ1 1 ÿ ÿ 4 IIIWII;III C Cÿ1 Cÿ1 C 4 WI;IV IWII;IV 1 A A A A 1 4 WI;V IWII;V 1 M M 1 ÿ 4 WII;IV C A A A A C ÿ 4 WII;V C M M C ÿ 4 IIIWIII;IV Cÿ1 A A A A Cÿ1 ÿ 4 IIIWIII;V Cÿ1 M M Cÿ1 4 WIV;V M A A A A M 4 WIV;IV A A A A 4 WV;V M M ÿ 4 WII I 4 WV IA ÿ 4 IIIWIII ICÿ1
31
is large. The structure is built up by scalar factors and fourth order tensor quantities. The factors consist of ®rst and second order derivatives WI;I : oW2 =oIoI, WI;II : oW2 =oIoII, WI;III : oW2 =oIoIII etc. of the elastic potential W with respect to the ®ve invariants. The tensors are formed by combinations of the second order unit tensor 1, the Right Cauchy Green tensor C, its inverse Cÿ1 and the unit vector in ®bre direction A. Moreover, Eq. (31) contains the fourth order quantities I, IA , ICÿ1 (see Appendix A). The results (29) and (31) are transformed to the current con®guration by a push forward a U A using the deformation gradient F. The relations may also be derived directly in the current con®guration utilizing a potential function W W b; a a given with respect to the deformed body, where b : FFT is the Left Cauchy Green tensor and a FA is the current ®bre direction. The spatial approach is not considered further in this work because the ®nite viscoelastic formulation (Section 3.2) has to be developed with respect to the undeformed body.
232
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
p A lot of materials show dierent behaviour for volumetric J det F III and for isochoric states of T deformation C : F F. In order to take this observation into account, an additive division of the strain energy function W U J W C; A A
32
is introduced, where det C J ÿ2=3 det C 1. In contrast to function (24), which describes a compressible material, a separate determination of the volumetric term U J and the incompressible part W C; A A is carried out. This representation is based on the multiplicative split of the deformation gradient F J 1=3 1F into a volumetric J 1=3 1 and an incompressible component F. Subsequently, the invariants 1 2 2 I tr C; II tr2 C ÿ tr C ; IV A : C : A; V A : C : A 2
33
are computed analogously to Eq. (25), however, the incompressible Right Cauchy Green tensor C is employed. The third invariant, standing for the volume change, results in III det C 1. Therefore, it is meaningless for the constitutive relation W C; A A. With the formulation shown (32), a material can be modelled being compressible J 6 1 or nearly incompressible J 1 by choosing the material parameters. Special numerical techniques allow the stimulation of exact incompressibility (see [22] for details). An uncoupled structure of the second Piola Kirchho stress tensor is derived S Svol Siso
34
and, equivalently, A Avol Aiso
35
for the tangential material tensor founded on the assumption of a splitted strain energy (32). A distinction between the volume deformation J and the isochoric contribution C is strictly considered. In the reference con®guration the quantities for the hydrostatic pressure Svol JU 0 Cÿ1 ;
36
where U 0 : oU =oJ , and the deviatoric stress component ~ Siso J ÿ 2=3 DEV S
37
using the result ~ : 2 oW 2 W IW 1 ÿ 2 W C 2 W A A 2 W M S I II II IV V oC
38
form the second Piola Kirchho stress tensor. Eq. (38) employs M : A CA CA A
39
and Eq. (37) the de®nition 1 DEV : ÿ C : Cÿ1 ; 3 which gives the equivalent quantity of a second order deviator in the reference con®guration. For the volumetric tangent the standard expression ÿ ÿ Avol J 2 U 00 Cÿ1 Cÿ1 JU 0 Cÿ1 Cÿ1 ÿ 2ICÿ1
40
41
is calculated from the strain energy function and lengthy considerations of tensor analysis yield the isochoric tangent operator
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
ÿ ÿ 2 ~ : C I ÿ1 ÿ 1 Cÿ1 Cÿ1 ÿ 2 Siso Cÿ1 Cÿ1 Siso J ÿ4=3 DEV A ~ iso : Aiso J ÿ2=3 S C 3 3 3
233
42
~ iso is given by the material tensor The fourth order deviatoric quantity in the reference con®guration DEVA ~ iso : 4 W WII 2IW I2 W 1 1 A I;I I;II II;II ÿ 4 WI;II ÿ 2IWII;II 1 C C 1 4 WII;II C C 4 WI;IV IWII;IV 1 A A A A 1 4 WI;V IWII;V 1 M M 1 ÿ 4 WII;IV C A A A A C ÿ 4 WII;V C M M C 4 WIV;V M A A A A M 4 WIV;IV A A A A 4 WV;V M M ÿ 4 WII I 4 WV IA
43
similar to the compressible model (31). The general structure of the material tensor (43) is reduced by the terms, which contain the third invariant III. Properties associated with the volumetric deformation are described merely by the function U J . Parallel to the second order quantity (40) a kind of fourth order deviator is de®ned 1 1 1 DEV : ÿ Cÿ1 : C ÿ : C Cÿ1 C : : CCÿ1 Cÿ1 3 3 9
44
in order to compute (42). The general results (Eqs. (41) and (42)), which formulate the decoupled material behaviour, can be found, for example, in Simo and Hughes [37]. 3.2. Viscoelastic constitutive model The extension of the constitutive approach to ®nite strains introduced in Section 2.2 is accomplished using the second Piola Kirchho stress tensor S C; A A S0
n X
Hj :
45
j1
The generalized Maxwell-element is the symbolic rheological structure used and, therefore, it yields an additive composition of the stresses. The linear rate-equation _ j 1 Hj Cj S_ 0 H sj
46
introduces the constitutive relation between the internal viscoelastic stress variables Hj and the elastic stress rate S_ 0 . Time-dependency is considered by the hereditary integral Z t t ÿ t oS0 C; A A Cj A A exp ÿ dt 47 Hj C; A A ot sj 0
234
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
developed from Eq. (46). Equivalent to the small strain case, the stress rate dS0 DS0 Sn1 ÿ Sn0 lim lim 0 S_ 0 Dt!0 Dt Dt!0 dt Dt is discretized for a time increment Dt. The resulting formula ÿ 1 ÿ exp ÿ Dt=sj n1 Dt n1 n S0 ÿ Sn0 Hj Cj Hj exp ÿ sj Dt=sj
48
49
is an ecient approach for an approximate update of the anisotropic viscoelastic stress variables. The quantities Hnj and Sn0 of time tn have to be supplied for the computation of the present stresses. Relaxation time sj and the fourth order relaxation tensor Cj are material parameters to be determined. The constant tensor Cj A A Aj C 1; A AAÿ1 0 C 1; A A const:
50
is evaluated for the underformed material, i.e., C 1. The constitutive dierential equation (46) de®nes a linear relation between Hj and S_ 0 by the relaxation tensor Cj independent of the deformation. Nonlinearity is introduced by the ®nite strain measure C and the application of nonlinear elastic characteristics using an appropriate potential function W W C; A A. This way, the application of the relaxation tensor allows a separation of linear and nonlinear parts of the model. Moreover, this elastic and viscous split is also important for parameter identi®cation. The dierentiation of the algorithmic stress tensor with respect to strain yields the consistent tangent operator ( ÿ ) n X 1 ÿ exp ÿ Dt=sj oSn1 n1 Cj A A 51 An1 A : 2 n1 I 0 C; A A; Dt=s oC j j1 which is determined by the constant tensorial relaxation function n X t Cj exp ÿ C t I sj j1
52
and the nonlinear elastic material tensor A0n1 . For a time step Dt the algorithmic tensor (51) is divided up into a constant viscous term fg and the nonlinear elastic contribution An1 0 . This fourth order tensor A C 1; A A is symmetric for the undeformed material. In case of loading, i.e., C 6 1, an unsymmetric quantity A C; A A is obtained. Therefore, the entire ®nite element matrices of the discrete problem have to be evaluated and assembled. The implementation of (51) in a ®nite element code ensures a quadratic rate of convergence of the Newton procedure at least near the ®nal solution. The tangential material tensor in the current con®guration a is calculated by a push forward operation a U A
53
utilizing the deformation gradient. A direct analytical determination of the spatial tensor a is not possible in contrast to the isotropic formulation, where the tensorial relaxation function C t reduces to a scalar function. 4. Numerical examples 4.1. Dynamic response of U-pro®le The ®rst application of the introduced models simulates the geometrically linear dynamic response of a composite pro®le. The geometry (U 200 60 10 mm3 ) is taken from the design manual of Fiberline Composites A/S [3]. Glass ®bres, which make up 50% of the material's volume, are arranged in x-direction and combined with a polymer matrix. Consequently, material characteristics in y- and z-orientation are
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
235
Fig. 2. Discretization and deformed con®guration of U-pro®le.
Fig. 3. Dynamic response of U-pro®le.
isotropic. The pro®le l 600 mm is represented by 380 linear brick elements (see Fig. 2). All nodes of the rear cross-section are ®xed and the tip of this cantilever beam is subjected to a distributed load of q 125 N=mm constant with regard to time. Fig. 3 shows the horizontal and the vertical response of marked point P as function of time. Firstly, the pro®le is simulated using a viscoelastic composite material and, secondly, compared to an isotropic steel construction. In the latter case, dissipative behaviour of the structure is not given by the constitutive
236
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
Fig. 4. Rayleigh-damping of steel.
relation but de®ned utilizing Rayleigh-damping. The two damping parameters are taken to yield a nearly frequency-independent dissipation for high modes and a structural damping ratio of n 0:001273 (Fig. 4) according to Petersen [16] for the dominant torsional mode at xsteel 1314 rad=s. The comparison of the two types of pro®le point out signi®cant dierences. The steel pro®le is less deformed staticly and dynamically because of a larger Young's modulus. The vibrations of the composite cantilever are strongly damped owing to the viscoelastic polymer resin. The lower stiness of the composite material results in an oscillation at a reduced main frequency xcomposite 752 rad=s compared to the steel pro®le xsteel 1314 rad=s. The dominant eigenmode is a torsional deformation (Fig. 2) because the resultant of the exciting load is not located at the shear center. 4.2. Elastic analysis of a human cornea In the introduction it was pointed out that a lot of biological tissues are characterized by a ®brous microstructure and, therefore, anisotropic constitutive formulations apply to it. Parts of the human eye, investigated in this numerical example, fall into this class of materials. Fig. 5 depicts a partial discretization with eight noded 2304 mixed Q1/P0-brick elements (according to [22]). The inner part (cornea) is represented by a ®ne mesh and the outer part (sclera) is given by a coarse discretization. The state of deformation of the cornea is of primary interest in this simulation and the sclera is merely considered as elastic boundary condition for the cornea.
Fig. 5. Finite element discretization of the human eye.
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
237
The cornea, which is very important part of the optical system, is built of some hundreds of transparent collagen ®bre layers and it has a total thickness of 520 lm. One collagen ®bre is relatively sti in ®bre direction and very soft perpendicular to it. The ®bre layers are arranged with arbitrary direction in the cornea plane. The superposition of collagen layers with dierent orientation yields a stier isotropic plane and a reduced stiness rectangular to it. In order to obtain an ecient description, a transverse isotropic potential function W U vol J W iso C W aniso C; A A
54
is introduced. The potential is split up into a volumetric and into an isochoric contribution to model quasiincompressibility. The isochoric part consists of a small isotropic basic stiness W iso which is superimposed by an isochoric anisotropic part W aniso
m X n 1 X i1 C a IVj ÿ 1 : m j1 i1 i
55
The ®bre stiness (55) is de®ned as an expansion of the fourth invariant IVj IV C; A aj
56
using n terms, where aj j ÿ 1 p=m for j 1; . . . ; m. The ®bre is rotated in the cornea plane in m steps and superimposed, i.e., the mean stiffness for m ®bre angles aj is computed on the material scale. A ®rst simulation investigates the state of deformation of the cornea at inner pressure pi . Fig. 6 depicts a cross-section of the undeformed ideal cornea and the state of deformation at pressure pi 0:007 MPa (deformation plot is scaled 2.5 times). The accompanying radius of curvature, which is a very important parameter in opthalmology, is shown in Fig. 8. The deformed con®guration (Fig. 6) as well as the plot of the radius of curvature (Fig. 8) clearly show a ¯attening of the cornea. The increase of the radius of curvature, which is constant at r 7:86 mm for the undeformed structure, is larger at the outer part of the cornea compared to the center. The postprocessing algorithm evaluating the curvature from the state of deformation yields a result that is erroneously not axisymmetric, as can be seen from Fig. 8. The algorithm is very sensitive and we compute a result in¯uenced by the ®nite element discretization, which also is not axisymmetric. A second ®nite element simulation (Fig. 7) analyses the in¯uence of a single stitch on the deformation behaviour, i.e., on the radius of curvature (Fig. 8). The nearly rigid stitch is discretized by one-dimensional
Fig. 6. Cross-section of the undeformed and the deformed con®guration of the cornea.
Fig. 7. Position of a single stitch.
238
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
Fig. 8. Radius of curvature (ideal cornea subjected to internal pressure).
Fig. 9. Radius of curvature (cornea with a single stitch subjected to internal pressure).
truss elements implemented at the place where a stitch is located if a new cornea was implanted by eye surgery. Fig. 9 shows the resulting curvature at pressure pi 0:007 MPa. The diagram depicts the local character of the pertubation. About three quarters of the cornea do not show any in¯uence of the stitch. The ®nite element results computed in this study compare qualitatively well with experimental investigations to be published in the future (see [31]).
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
239
4.3. Viscoelastic simulation of a two layer laminate A numerical study investigating a two layer tensile test is carried out in order to illustrate on the one hand the eect of an unsymmetric stacking order on the state of deformation and on the other hand the ®nite viscoelastic model. Geometry and boundary conditions for the two samples (®bre angle: a 90 = 0 , a 45 = ÿ 45 ) are identical. The specimen 400 100 16 mm3 is ®xed at one crosssection and loaded displacement controlled at the other end. The simulation is based on polyester ®bres and elastomeric material. Load-strain characteristics of the material components are given by Fig. 10. The diagram is evaluated for a single ®bre which can be loaded up to e 16:5%. The second graph (Fig. 10) depicts a typical nonlinear material curve of rubber: softening at the beginning followed by stiening. The composite material is de®ned by a constitutive transversely isotropic potential according to Eq. (54). The volumetric term U vol
j 2
J2 ÿ 1 ÿ ln J 2
57
is associated with the quasi-incompressible volume deformation J 1. The characteristics of the isochoric part are determined by the isotropic function 2
W iso C1 I ÿ 3 C2 I ÿ 3 C3 I ÿ 3
3
58
introduced by Yeoh [32], which is a realistic phenomenological model for rubber elasticity (see [9]). The expansion W aniso
n X i1
Cia IV ÿ 1
i1
59
in n terms of the fourth invariant IV describes the nonlinear stiffness in ®bre direction. The numerical model of the specimen (Fig. 11) consists of a discretization of 459 Q1/P0-brick elements. Fig. 11 shows the deformed con®gurations at an elongation of 66 mme16:5%. For layers a 90 = 0 the state of deformation is scaled three times while for a 45 = ÿ 45 the deformation plot is not enlarged. The deformed specimen shows either a bending-extension coupling a 90 = 0 or a torsion±extension coupling a 45 = ÿ 45 . The latter state of deformation results in a twisting of about 90 . In a displacement controlled numerical cyclic test the reaction force F is computed as function of the displacement u at the cross-section subjected to deformation (Fig. 12). The ®nite transversely isotropic viscoelastic model is employed. The load±displacement curves form a hysteresis because of the application of the inelastic material model. The area described by the curves stands for the dissipated energy per loading cycle. The difference between the viscoelastic stiffness of the two ®bre con®gurations can be seen clearly.
Fig. 10. Experimental investigation of polyester ®bre and rubber.
240
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
Fig. 11. Deformed con®gurations of two layer laminate (e 16:5%).
Fig. 12. Cyclic extension test.
5. Conclusions This paper has presented constitutive formulations for material that is idealized as transversely isotropic. The approaches cover material with ®brous microstructure like biological tissues and all kinds of undirectional ®bre matrix combinations. Elasticity and viscoelasticity are considered as two main aspects of the material characteristics. Depending on application and composite components, it has to be decided whether in®nitesimal or ®nite strains apply. The transversely isotropic concepts shown can be extended to fully orthotropic description if necessary. A reduction of the stiness depending on loading, summarized as `damage phenomena', is not considered in the present paper. Future work will add the micromechanical change of the material structure due to ®bre or matrix rupture and molecular debonding like Mullins-eect for rubber material. An anisotropic damage model is capable of predicting this behaviour. The present article focuses on the homogeneous macroscopic material level in order to provide a formulation suitable for structural ®nite element simulations. The ®eld of parameter identi®cation and determination of material characteristics from the material's heterogenous micro structure is excluded and will be discussed in a second study. In principle, the mean properties can be evaluated analytically, numerically and experimentally. For small strains, very ecient homogenization procedures either in the timedomain, or in the Laplace-domain, depending on elastic or viscoelastic characterization, have been derived and veri®ed by ®nite element simulations on the micro scale. At ®nite strains homogenization is quite dicult and still an area of research.
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
241
Appendix A A.1. Tensor invariants Based on an arbitrary second order tensor B and a unit vector A ®ve irreducible invariants IB tr B;
A:1
1 IIB tr2 B ÿ tr B2 ; 2
A:2
IIIB det B;
A:3
IVB A : B : A;
A:4
VB A : B2 : A
A:5
are formulated, which can be used for example to de®ne a transversely isotropic strain energy function (Eq. (26)). In order to develop the complete material formulation, partial derivatives oIB 1; oB
A:6
oIIB IB 1 ÿ B; oB
A:7
oIIIB IIIB Bÿ1 ; oB
A:8
oIVB A A; oB
A:9
oVB A BA BA A oB
A:10
of the invariants with respect to tensor B are needed (see Section 3.1). The speci®c tensors shown below are given with respect to a Cartesian basis in symbolic form and they are also depicted with its components. These quantities are employed to compute the tangential material tensor. The derivative of tensor B with respect to itself yields a fourth order unit tensor oB 1 I Iijkl ei ej ek el dik djl dil djk ei ej ek el ; oB 2
A:11
where dij is the Kronecker symbol. The derivative of the inverse tensor Bÿ1 with respect to B oBÿ1 1 ÿ1 ÿIBÿ1 ÿIBÿ1 ;ijkl ei ej ek el ÿ Bÿ1 Bÿ1 Bÿ1 il Bjk ei ej ek el ; 2 ik jl oB
A:12
results in a fourth order tensor equivalent to the unit tensor but utilizing the components of Bÿ1 . Computing the derivative of M A B A B A A with respect to tensor B, i.e., evaluating the second derivative of the ®fth invariant with respect to B, leads to a constant fourth order tensor
242
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
oM 1 IA IA;ijkl ei ej ek el Ak Ai djl Ai Al djk Aj Ak dil Aj Al dik ei ej ek el ; oB 2
A:13
assembled from the components of vector A. A.2. Transformation operator A transformation from Cartesian material coordinates ~xi the transformation operator 2 2 m21 n21 l1 m 1 l1 n1 l1 6 l2 2 2 m n l m l2 n2 6 2 2 2 2 2 6 2 2 2 6 l3 m3 n3 l3 m 3 l3 n3 T6 6 2l l 2m m 2n n l m l m l n l n 1 2 1 2 1 2 2 1 1 2 2 1 6 1 2 6 4 2l1 l3 2m1 m3 2n1 n3 l1 m3 l3 m1 l1 n3 l3 n1 2l2 l3 2m2 m3 2n2 n3 l2 m3 l3 m2 l2 n3 l3 n2
to global coordinates xi is accomplished with 3 m 1 n1 7 m 2 n2 7 7 7 m 3 n3 7; m1 n2 m2 n1 7 7 7 m1 n3 m3 n1 5 m2 n3 m3 n2
A:14
where li ; mi ; ni are the direction cosines.
References [1] R.M. Christensen, Mechanics of Composite Materials, Wiley, New York, 1979. [2] T.C. Doyle, J.L. Ericksen, Nonlinear elasticity, in: Advances in Applied Mechanics IV, Academic Press, New York, 1956. [3] Fiberline, Design- und Konstruktionshandbuch f ur Konstruktionspro®le aus Verbundwerkstoen, Fiberline Composites A/S, Kolding, 1995. [4] R.F. Gibson, Principles of Composite Material Mechanics, McGraw-Hill, New York, 1994. [5] S. Govindjee, J.C. Simo, Mullins' eect and the strain amplitude dependence of the storage modulus, International Journal of Solids and Structures 29 (1992) 1737±1751. [6] L.R. Herrmann, F.E. Peterson, A numerical procedure for viscoelastic stress analysis, in: Proceedings of the Seventh Meeting of ICRPG Mechanical Behaviour Working Group, Orlando, 1968. [7] R.M. Jones, Mechanics of Composite Materials, McGraw-Hill, Tokyo, 1975. [8] M. Kaliske, H. Rothert, Formulation and implementation of three-dimensional viscoelasticity at small and ®nite strains, Computational Mechanics 19 (1997) 228±239. [9] M. Kaliske, H. Rothert, On the ®nite element implementation of rubberlike materials at ®nite strains, Engineering Computations 14 (1997) 217±233. [10] Y.C. Lou, R.A. Schapery, Viscoelastic characterization of a nonlinear ®ber reinforced plastic, Journal of Composite Materials 5 (1971) 208±234. [11] J. Lubliner, A model of rubber viscoelasticity, Mechanics Research Communications 12 (1985) 93±99. [12] T.C. Kennedy, M. Wang, Three-dimensional, nonlinear viscoelastic analysis of laminated composites, Journal of Composite Materials 28 (1994) 902±925. [13] K.-D. Klee, E. Stein, Discretization in nonlinear continuum mechanics and comparison of some ®nite element algorithms, in: Proceedings of the Pressure Vessels and Piping Conference, Orlando, 1982. [14] K.Y. Lin, I.H. Hwang, Thermo-viscoelastic analysis of composite materials, Journal of Composite Materials 23 (1989) 554±569. [15] A.I. Lurie, Nonlinear theory of elasticity, in: J.D. Achenbach et al. (Eds.), North-Holland Series in Applied Mathematics and Mechanics, North-Holland, Amsterdam, 1990. [16] C. Petersen, Dynamik der Baukonstruktionen, Vieweg, Braunschweig, 1996. [17] H. Poon, F.M. Ahmad, A material point time integration procedure for anisotropic thermo rheologically simple, viscoelastic solids, Computational Mechanics 21 (1998) 236±242. [18] J. Schr oder, Theoretische und algorithmische Konzepte zur ph anomenologischen Beschreibung anisotropen Materialverhaltens, Dissertation, Universit at Hannover, 1996. [19] F.R. Schwarzl, A.J. Staverman, Time±temperature dependence of linear viscoelastic behaviour, Journal of Applied Physics 23 (1952) 838±843. [20] F. Sidoro, Un modele viskoelastique non lineaire avec con®guration intermediaire, Journal de Mecanique 13 (1974) 679±713. [21] J.C. Simo, On a fully three-dimensional ®nite-strain viscoelastic damage model: formulation and computational aspects, Computer Methods in Applied Mechanics and Engineering 60 (1987) 153±173. [22] J.C. Simo, R.L. Taylor, Quasi-incompressible ®nite elasticity in principal stretches, continuum basis and numerical algorithms, Computer Methods in Applied Mechanics and Engineering 85 (1991) 273±310.
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
243
[23] F. Soergel, Biomechanische Charakterisierung der menschlichen Augenhornhaut mit dynamisch-mechanischer Spektroskopie, Dissertation, Universit at Ulm, 1994. [24] A.J.M. Spencer, Continuum theory of the mechanics of ®bre-reinforced composites, CISM Courses and Lectures, No. 282, Springer, Wien, 1984. [25] E.S. Suhubi, Thermoelastic solids, in: A.C. Eringen (Ed.), Continuum Physics II, Academic Press, New York. [26] P. Le Tallec, C. Rahier, A. Kaiss, Three-dimensional incompressible viscoelasticity in large strains: formulation and numerical approximation, Computer Methods in Applied Mechanics and Engineering 109 (1993) 233±258. [27] R.L. Taylor, K.S. Pister, G.L. Goudreau, Thermomechanical analysis of viscoelastic solids, International Journal for Numerical Methods in Engineering 2 (1970) 45±59. [28] C. Truesdell, W. Noll, The nonlinear ®eld theories, in: S. Fl ugge (Ed.), Handbuch der Physik, Band III/3, Springer, Berlin, 1965. [29] M. Wang, Three-dimensional, nonlinear viscoelastic analysis of laminated composites: a ®nite element approach, Dissertation, Oregon State University at Corvallis, 1993. [30] J.A. Weiss, N.M. Bradley, S. Govindjee, Finite element implementation of incompressible, transversely isotropic hyperelasticity, Computer Methods in Applied Mechanics and Engineering 135 (1996) 107±128. [31] T. Wolf, Druckinduzierte Anderungen der Hornhautkr ummungen bei Schweinebulbie mit Einzelknopfnaht und fortlaufender Keratoplastiknaht, Dissertation, Medizinische Hochschule Hannover (to be published). [32] O.H. Yeoh, Characterization of elastic properties of carbon-black-®lled rubber vulcanizates, Rubber Chemistry and Technology 63 (1990) 792±805. [33] S. Yi, Finite element analysis of anisotropic viscoelastic composite structures and analytical determination of optimum viscoelastic material properties, Dissertation, University of Illinois at Urbana-Champaign, 1992. [34] S. Yi, H. Hilton, Dynamic ®nite element analysis of viscoelastic composite plates in the time domain, International Journal for Numerical Methods in Engineering 37 (1994) 4081±4096. [35] S. Yi, H. Hilton, Hygrothermal eects on viscoelastic responses of laminated composites, Composites Engineering 5 (1995) 183± 193. [36] M.A. Zocher, S.E. Groves, A three-dimensional ®nite element formulation for thermoviscoelastic orthotropic media, International Journal for Numerical Methods in Engineering 40 (1997) 2267±2288. [37] J.C. Simo, T.J.R. Hughes, Computational Inelasticity, Interdisciplinary Applied Mathematics 7, Springer, New York, 1998.
www.elsevier.com/locate/cma
A formulation of elasticity and viscoelasticity for ®bre reinforced material at small and ®nite strains M. Kaliske Institut f ur Statik, Universit at Hannover, Appelstraûe 9a, 30167 Hannover, Germany Accepted 26 March 1999
Abstract The combination of unidirectional ®bres embedded in matrix yields a novel material which is heterogeneous on the micro level and exhibits new mean characteristics on the macro level. Due to the fact that the length scale of a composite structure is several orders larger compared to the micro scale, the composite material is assumed to be homogeneous on the structural level. Therefore, a homogeneous anisotropic constitutive formulation has to be introduced in order to represent the material and the composite structure eciently, for example by a ®nite element discretization. In case of unidirectional reinforced composites a transversely isotropic idealization is meaningful. One typical class of composites is manufactured, for example, from glass or carbon ®bres combined with an epoxy matrix. Usually, structures produced from these materials are subjected to small strains. A geometrically nonlinear description is essential when modelling steel or nylon cord reinforced elastomeric material. The investigation of biomaterials like muscles or bones, which also have a ®brous microstructure, is accomplished on the basis of ®nite strains, too. The static and the dynamic response of a structure can be simulated with the ®nite element method. The stiness is characterized by anisotropic elastic properties. For the dynamic response, the dissipative behaviour, originated in the viscoelastic properties of the polymeric material, is of great importance. Depending on the class of composites investigated, the elastic and the viscoelastic behaviour of ®bre reinforced material has to be formulated by a transversely isotropic approach at small and at ®nite strains respectively. With these models at hand, realistic ®nite element simulations may be carried out. Ó 2000 Elsevier Science S.A. All rights reserved. Keywords: Elasticity; Viscoelasticity; Fibre reinforced material
1. Introduction The status of composite material in engineering application is increasing steadily. Fibre matrix combinations are employed in the aircaft industry and for the construction of space structures. Moreover, composites are also used for civil engineering structures like footbridges. This class of `sti composites', made for example from glass ®bres in epoxy resin is enlarged by so-called `soft composites'. This second type of material is formed by cords with high tensile strength reinforcing an elastomer. Typical areas of application are tires, conveyer belts or pneumatic shock absorbers. The physical description of the industrial products mentioned has to be extended to biomaterials which often exhibit an anisotropic microstructure. Therefore, constitutive models for bones or soft living tissues belong to the same class as the arti®cially produced material. The manufacturing of composites allows the creation of speci®c characteristics of the combined material through variation of geometrical parameters like ®bre volume fraction, ®bre angle and even ®bre shape or through the combination of dierent components. This way, a material design is accomplished to optimize structural characteristics. 0045-7825/00/$ - see front matter Ó 2000 Elsevier Science S.A. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 9 9 ) 0 0 2 6 1 - 3
226
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
A numerical tool to investigate structures eciently is the ®nite element method. By this approach, it is possible to solve complex initial boundary value problems. In order to compute realistic results a constitutive formulation has to be employed that incorporates all main features of the material. Principle properties to be taken into account are anisotropic elasticity and anisotropic viscoelasticity. The description of linear elastic behaviour of ®bre reinforced material is shown in standard text books on composites. Assuming linear elasticity and small strains Jones [7], Christensen [1] or Gibson [4] for example present the constitutive relations for dierent types of anisotropy based on engineering constants. For the geometrically nonlinear regime the linear elastic approach is often used as to be seen in the literature, i.e., modelling is accomplished by the so-called de St. Venant Kirchho material. Klee and Stein [13] show that the formulation may be employed just for a small range of deformation, i.e., up to 5% of strain. General considerations on anisotropic elasticity at ®nite strains are given in Truesdell and Noll [28], Suhubi [25] and Lurie [15]. Spencer [24] developed a constitutive strain energy function in terms of invariants. For this purpose, he extended the basis of three characteristic quantities of isotropy to ®ve invariants of transverse isotropy. Schr oder [18] and Weiss et al. [30] investigate ®nite and anisotropic nonlinear elasticity. Moreover, Weiss et al. [30] show a speci®c strain energy function applied to ®nite element simulations of biological material. Time- and frequency-dependent behaviour is observed not only for an isotropic polymer like epoxy resin or rubber, but the composite itself is also viscoelastic. Biomaterials, which are mostly anisotropic, exhibit viscoelastic behaviour as well (see for example [23,26]). Dierent theoretical models and numerical realizations for the formulation of anisotropic linear viscoelastic properties of thermo-rheologically simple materials (see [19]) at small strains have been published. The approach shown below is motivated by its applicability to a linear as well as to a geometrically nonlinear description. An overview of the present literature on modelling of anisotropic viscoelastic material reponse reveals that most of the authors follow similar ideas. The article of Lin and Hwang [14] is one of the ®rst publications dealing with the mentioned topic in the context of the ®nite element method. They introduce a thermo-viscoelastic model for plane strain and its algorithmic formulation based on the recurrence relation of Taylor et al. [27]. Yi [33] computes the solution of the linear quasi-static anisotropic problem in the Laplace- and in the Fourier-domain. A time-dependent simulation is accomplished by Yi according to the algorithm of Taylor et al. [27] employing a symbolic generalized Maxwell-element. In a number of papers ([33±35]) the authors start with a viscoelastic variation principle, discretize it with the ®nite element method and ®nally end up with integral equations, which can be solved recursively. Wang [29] and Kennedy and Wang [12] developed a solution procedure for nonlinear integral equations utilizing the viscoelastic approach of Lou and Schapery [10] in order to describe the long-term nonlinear viscoelastic response of composite material. The formulation is applied to ®nite element computations. Zocher et al. [36] present an algorithmic approach to orthotropic linear thermo-viscoelasticity. The integral equations of the material history are developed in an incremental form and solved recursively similar to the method of Taylor et al. [27]. In a recent publication Poon and Ahmad [17] also introduce an eective and robust integration algorithm for anisotropic linear viscoelaticity which can considered to be similar to the model of Zocher et al. [36]. An extension of the small strain approach to large strain viscoelasticity is provided in the articles of Simo [21], Govindjee and Simo [5] and Kaliske and Rothert [8] for isotropic material. Basics of ®nite strain viscoelastic formulations are developed by Sidoro [20] and Lubliner [11]. Le Tallec et al. [26] derived a ®nite and nonlinear viscoelastic model and used the approach for ®nite element simulations of the human eye. The outline of this paper is the following: In Section 2 transversely isotropic formulations for elastic and viscoelastic characteristics at small strains are given. Subsequently, Section 3 extends the general structure of the approaches to the ®nite strain regime. Numerical examples of Section 4 illustrate the application of the models to ®nite element computations of a dynamic structural response, reinforced rubber material and simulations of biological tissues. Main features are summarized in Section 5. Finally, Appendix A provides some useful results of tensor analysis, which are employed in the following considerations.
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
227
2. Linear theory 2.1. Elastic constitutive model Hyperelastic material is de®ned by the strain energy function W. In the case of a geometrically linear and isotropic theory, it is a function of the small strain tensor e. A unidirectional ®bre reinforced composite is characterized by a unit vector in ®bre direction a, where kak 1. Transverse isotropy is obtained for ®bres in a hexagonal order, i.e., material properties are invariant with respect to a rotation about the ®bre axis a. Due to this second constitutive characteristic a, the strain energy function W of a composite is formulated in terms of the strain tensor e and the dyadic product a a, which is a second order tensor of the microstructure, i.e., the ®bre orientation.The invariance requirement of the strain energy W e; a a W QeQT ; Qa aQT
1
has to be ful®lled, if e and a are subjected to a rotation by a proper orthogonal tensor Q e0 QeQT
a0 Qa;
and
2
where QQT QT Q 1 with the unit tensor 1 of second order and det Q 1. Spencer [24] derived the basis tr e; tr e2 ; tr e3 ; a : e : a; a : e2 : a
3
of ®ve irreducible invariants of transverse isotropy. Subsequently, the most general quadratic strain energy function k b W tr2 e lT tr e2 a aeatr e 2 lL ÿ lT ae2 a aea2 2 2
4
is formulated using the ®ve independent material parameters k, lT , a, lL , b. Standard evaluations of continuous mechanics and tensor analysis lead to the accompanying stress tensor r
oW ktr e a aea1 2lT e atr e b aea a a 2 lL ÿ lT a ae ea a oe
5
as derivative of the potential function W with respect to the strain tensor e. The linear elastic material tensor a
o2 W k 1 1 a a a 1 1 a a 2lT I b a a a a 2 lL ÿ lT Ia oe2
6
is computed as second derivative of the potential function with respect to e, where 1 is a second order and I is a fourth order unit tensor. The quantity Ia stands for a fourth order tensor formed by the components of the unit vector a (see Eq. (A. 13) in Appendix A for an explicit representation). Supposing a ®bre orientation in ~x1 direction, i.e., aT 1; 0; 0, the constant material tensor (6) reduces to 2 6 6 6 6 6 6 6 ~a 6 6 6 6 6 6 6 6 4
k 2a 4lL ÿ 2lT b
ka
ka
0
0
ka
k 2lT
k
0
0
ka
k
k 2lT
0
0
0
0
0
2lL
0
0
0
0
0
2lL
0
0
0
0
0
0
3
7 7 0 7 7 7 7 0 7 7 7 7 0 7 7 7 0 7 7 5 2lT
7
228
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
in matrix form. The inverse ~ aÿ1 of the 2 1 m12 m12 ÿ ÿ 6 E11 E E 11 11 6 6 m12 1 m23 6ÿ ÿ 6 E11 E22 E22 6 6 m12 m23 1 6ÿ 6 E11 ÿ E22 E22 ~aÿ1 6 6 6 0 0 0 6 6 6 6 0 0 0 6 6 4 0 0 0
classical formulation employing engineering constants yields 3 0 0 0 7 7 7 0 0 0 7 7 7 7 0 0 0 7 7 7 7 1 0 0 7 7 G12 7 7 1 0 7 0 7 G12 7 1 5 0 0 G23
8
for a transversely isotropic composite in material coordinates ~xi (see for example [7]). The shear modulus G23 in Eq. (8) is not an independent parameter and may be computed according to G23
E22 2 1 ÿ m23
9
from Young's modulus E22 and Poisson's ratio m23 . The strain energy function of the classical approach 1 W e:a:e 2
10
is given directly by the strain e and the material tensor a. Since the shown linear formulations in Eqs. (7) and (8) both describe transversely isotropic material, the engineering constants (Eq. (8)) and the coecients of the invariant based approach (Eq. (7)) are equivalent (see [18] for the conversion). 2.2. Viscoelastic constitutive model The so-called generalized Maxwell-model has proved to be a useful rheological structure to formulate time-dependent and frequency-dependent characteristics. The symbolic material model is composed of an elastic spring in parallel with n separate Maxwell-elements j (Fig. 1). This approach results in a stress tensor n X hj 11 r r0 j1
formed additively from purely elastic stresses r0 a0 e and internal viscoelastic stress variables hj , which are determined by the linear rate-equation 1 h_ j hj cj r_ 0 sj
12
for a symbolic material element j. Eq. (12) de®nes the viscoelastic stress variables hj as function of the stress rate r_ 0 . Material parameters are the relaxation time sj , which is chosen to be independent of the ®bre
Fig. 1. Generalized Maxwell-element.
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
229
orientation and, consequently, results in a scalar property. The model is charaterized by the fourth order relaxation tensor cj . This quantity represents the anisotropic viscosity of the material. A similar isotropic formulation employing a scalar factor cj , is discussed by Simo [21] and Kaliske and Rothert [8]. The relaxation time (Eqs. (12) and (46)) used in the mathematical model do not necessarily have to be identical for the different components of the material tensor. In this respect Yi and Hilton [35], Zocher et al. [36] and Poon and Ahmad [17] consider the relaxation time of the anisotropic material as a direction-controlled quantity. The transversely isotropic relaxation tensor for one Maxwell-element is given in material coordinates ~xi by the matrix 2
c11
6 6 c21 6 6 6 c21 ~cj 6 6 6 0 6 6 4 0 0
c12
c12
0
0
c22
c23
0
0
c23
c22
0
0
0
0
c44
0
0
0
0
c44
0
0
0
0
0
3
7 0 7 7 7 0 7 7: 7 0 7 7 7 0 5
13
c66
In contrast to the symmetric elasticity tensor a, the relaxation tensor is unsymmetric for the components c12 and c21 . Due to the speci®c symmetries of transverse isotropy c23 c32 is found. In case of isotropy, the approach reduces ~cj cj I
14
to a scalar factor ci and the fourth order identity tensor I. Eq. (12) carries over to the well known isotropic rate-equation. The fourth order relaxation tensor (13) is computed in a formal step according to ~cj ~aj a~ÿ1 0 ;
15
analogously to the isotropic model, where cj lj =l0 de®nes the contribution of one Maxwell-element j to the total stiffness. The scalar parameters l0 and lj are the elastic constants of the symbolic spring and the Maxwell-element j. Equivalently, ~ a0 in Eq. (15) stands for the purely elastic behaviour of the ®bre reinforced material and ~ aj is the transversely isotropic elasticity tensor of Maxwell-element j. The total timedependent stiffness of the generalized Maxwell-element is given by ( ) n X t cj exp ÿ 16 a0 c ta0 ; a t I sj j1 where c t is the tensorial relaxation function. The relaxation tensor ~cj is transformed either by using the ®bre direction a or by T ÿT ÿT aj T Tÿ1 ~ aÿ1 aÿ1 TT ~ aj ~ TT~cj TÿT cj aj aÿ1 0 T ~ 0 T 0 T
from local material coordinates to a global reference system (see Eq. (A.14) in Appendix A). Starting from the linear rate-equation (12) the hereditary integral Z t t ÿ t or0 cj exp ÿ dt hj sj ot 0
17
18
is derived. It describes the history of the internal stress variables hj . The integration with respect to time is accomplished eciently by a recurrence relation hjn1
exp
ÿ 1 ÿ exp ÿ Dt=sj n1 Dt n r0 ÿ rn0 ÿ hj c j sj Dt=sj
19
230
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
for time step Dt tn1 ÿ tn . This approach assumes a linear stress rate r_ 0
dr0 Dr0 rn1 ÿ rn0 lim lim 0 Dt!0 Dt Dt!0 dt Dt
20
during time step Dt. The recurrence formula is obtained exploiting the features of the exponential function as shown by Herrmann and Peterson [6] and Taylor et al. [27] in a similar way. For the evaluation of the recurrence relation (19) the known solutions hnj , rn0 at time tn are needed. This information has to be stored in a data base for one time step, i.e., for the evaluation at time tn1 .The algorithmic material tensor is crucial in the context of a ®nite element implementation. The fourth order tensor ( ÿ ) n X 1 ÿ exp ÿ Dt=sj orn1 n1 cj 21 an1 a : n1 I 0 oe Dt=s j j1 is determined as derivative of the algorithmic stresses rn1 with respect to strain en1 . Tensor (21) is symmetric. It shows an interesting structure because the viscous part fg and the elastic characteristics a0n1 are fully separated. This formal split is based on the introduction of the relaxation tensor cj . The advantage of the proposed formulation will be obvious for the extension to ®nite strains. The nonlinearities of a ®nite anisotropic linear viscoelasticity are associated with the nonlinear elastic behaviour and with the introduction of a ®nite strain measure. Viscoelasticity is still modelled by a linear rate-equation. Therefore, the separation of linear viscous and nonlinear elastic contributions is useful. 3. Finite strain formulation 3.1. Elastic constitutive model The invariant based elastic approach in Section 2.1 can be extended to a nonlinear transversely isotropic formulation at ®nite strains. On the other hand, an application of the classical linear description (Eq. (8)) to geometrically nonlinear considerations using a constant material tensor A yields a strain energy function 1 22 W E : A : E; 2 where E is the Green Lagrangian strain tensor. This physically linear so-called de St. Venant Kirchho material is valid only in a small range of strains. The strain energy function W of a nonlinear elastic ®bre reinforced material with respect to the reference con®guration is a function of the Right Cauchy Green tensor C and the ®bre direction A, i.e., the tensor A A of the microstructure of the material, where kAk 1. A push forward operation of the unit vector A to the current con®guration a U A FA
23
is accomplished by the deformation gradient F. The strain energy function for ®nite strains analogously to Eq. (1) is given by W W C; A A W QCQT ; QA AQT
24
postulating invariance with respect to a rotation Q. Spencer [24] derives the basis of irreducible invariants 1 25 I tr C; II tr2 C ÿ tr C2 ; III det C; IV A : C : A; V A : C2 : A 2 equivalently to the small strain formulation. The three invariants I, II, III of isotropy are extended by IV and V for transverse isotropy. Subsequently, the elastic potential W W I; II; III; IV; V
26
is formulated as function of the ®ve invariants.The second Piola Kirchho stress tensor S relative to the reference con®guration
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
S2
oW oC
231
27
and the accompanying fourth order material tensor A4
o2 W oC2
28
are computed according to standard methods (see [2]). The general form of the second Piola Kircho stress tensor for the elastic potential as given by Eq. (26) is derived S 2 WI IWII 1 ÿ 2 WII C 2III WIII Cÿ1 2 WIV A A 2 WV M
29
M : A CA CA A:
30
using The tensorial quantities are multiplied by scalar factors consisting of derivatives of the potential function with respect to the invariants WI : oW=oI; WII : oW=oII; WIII : oW=oIII etc. and by invariants I, III. The entire structure of the general material tensor A 4 WI;I WII 2IWI;II I2 WII;II 1 1 ÿ 4 WI;II ÿ 2IWII;II 1 C C 1 4 WII;II C C
ÿ 4 IIIWIII III2 WIII;III Cÿ1 Cÿ1 ÿ 4 IIIWI;III I IIIWII;III 1 Cÿ1 Cÿ1 1 ÿ ÿ 4 IIIWII;III C Cÿ1 Cÿ1 C 4 WI;IV IWII;IV 1 A A A A 1 4 WI;V IWII;V 1 M M 1 ÿ 4 WII;IV C A A A A C ÿ 4 WII;V C M M C ÿ 4 IIIWIII;IV Cÿ1 A A A A Cÿ1 ÿ 4 IIIWIII;V Cÿ1 M M Cÿ1 4 WIV;V M A A A A M 4 WIV;IV A A A A 4 WV;V M M ÿ 4 WII I 4 WV IA ÿ 4 IIIWIII ICÿ1
31
is large. The structure is built up by scalar factors and fourth order tensor quantities. The factors consist of ®rst and second order derivatives WI;I : oW2 =oIoI, WI;II : oW2 =oIoII, WI;III : oW2 =oIoIII etc. of the elastic potential W with respect to the ®ve invariants. The tensors are formed by combinations of the second order unit tensor 1, the Right Cauchy Green tensor C, its inverse Cÿ1 and the unit vector in ®bre direction A. Moreover, Eq. (31) contains the fourth order quantities I, IA , ICÿ1 (see Appendix A). The results (29) and (31) are transformed to the current con®guration by a push forward a U A using the deformation gradient F. The relations may also be derived directly in the current con®guration utilizing a potential function W W b; a a given with respect to the deformed body, where b : FFT is the Left Cauchy Green tensor and a FA is the current ®bre direction. The spatial approach is not considered further in this work because the ®nite viscoelastic formulation (Section 3.2) has to be developed with respect to the undeformed body.
232
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
p A lot of materials show dierent behaviour for volumetric J det F III and for isochoric states of T deformation C : F F. In order to take this observation into account, an additive division of the strain energy function W U J W C; A A
32
is introduced, where det C J ÿ2=3 det C 1. In contrast to function (24), which describes a compressible material, a separate determination of the volumetric term U J and the incompressible part W C; A A is carried out. This representation is based on the multiplicative split of the deformation gradient F J 1=3 1F into a volumetric J 1=3 1 and an incompressible component F. Subsequently, the invariants 1 2 2 I tr C; II tr2 C ÿ tr C ; IV A : C : A; V A : C : A 2
33
are computed analogously to Eq. (25), however, the incompressible Right Cauchy Green tensor C is employed. The third invariant, standing for the volume change, results in III det C 1. Therefore, it is meaningless for the constitutive relation W C; A A. With the formulation shown (32), a material can be modelled being compressible J 6 1 or nearly incompressible J 1 by choosing the material parameters. Special numerical techniques allow the stimulation of exact incompressibility (see [22] for details). An uncoupled structure of the second Piola Kirchho stress tensor is derived S Svol Siso
34
and, equivalently, A Avol Aiso
35
for the tangential material tensor founded on the assumption of a splitted strain energy (32). A distinction between the volume deformation J and the isochoric contribution C is strictly considered. In the reference con®guration the quantities for the hydrostatic pressure Svol JU 0 Cÿ1 ;
36
where U 0 : oU =oJ , and the deviatoric stress component ~ Siso J ÿ 2=3 DEV S
37
using the result ~ : 2 oW 2 W IW 1 ÿ 2 W C 2 W A A 2 W M S I II II IV V oC
38
form the second Piola Kirchho stress tensor. Eq. (38) employs M : A CA CA A
39
and Eq. (37) the de®nition 1 DEV : ÿ C : Cÿ1 ; 3 which gives the equivalent quantity of a second order deviator in the reference con®guration. For the volumetric tangent the standard expression ÿ ÿ Avol J 2 U 00 Cÿ1 Cÿ1 JU 0 Cÿ1 Cÿ1 ÿ 2ICÿ1
40
41
is calculated from the strain energy function and lengthy considerations of tensor analysis yield the isochoric tangent operator
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
ÿ ÿ 2 ~ : C I ÿ1 ÿ 1 Cÿ1 Cÿ1 ÿ 2 Siso Cÿ1 Cÿ1 Siso J ÿ4=3 DEV A ~ iso : Aiso J ÿ2=3 S C 3 3 3
233
42
~ iso is given by the material tensor The fourth order deviatoric quantity in the reference con®guration DEVA ~ iso : 4 W WII 2IW I2 W 1 1 A I;I I;II II;II ÿ 4 WI;II ÿ 2IWII;II 1 C C 1 4 WII;II C C 4 WI;IV IWII;IV 1 A A A A 1 4 WI;V IWII;V 1 M M 1 ÿ 4 WII;IV C A A A A C ÿ 4 WII;V C M M C 4 WIV;V M A A A A M 4 WIV;IV A A A A 4 WV;V M M ÿ 4 WII I 4 WV IA
43
similar to the compressible model (31). The general structure of the material tensor (43) is reduced by the terms, which contain the third invariant III. Properties associated with the volumetric deformation are described merely by the function U J . Parallel to the second order quantity (40) a kind of fourth order deviator is de®ned 1 1 1 DEV : ÿ Cÿ1 : C ÿ : C Cÿ1 C : : CCÿ1 Cÿ1 3 3 9
44
in order to compute (42). The general results (Eqs. (41) and (42)), which formulate the decoupled material behaviour, can be found, for example, in Simo and Hughes [37]. 3.2. Viscoelastic constitutive model The extension of the constitutive approach to ®nite strains introduced in Section 2.2 is accomplished using the second Piola Kirchho stress tensor S C; A A S0
n X
Hj :
45
j1
The generalized Maxwell-element is the symbolic rheological structure used and, therefore, it yields an additive composition of the stresses. The linear rate-equation _ j 1 Hj Cj S_ 0 H sj
46
introduces the constitutive relation between the internal viscoelastic stress variables Hj and the elastic stress rate S_ 0 . Time-dependency is considered by the hereditary integral Z t t ÿ t oS0 C; A A Cj A A exp ÿ dt 47 Hj C; A A ot sj 0
234
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
developed from Eq. (46). Equivalent to the small strain case, the stress rate dS0 DS0 Sn1 ÿ Sn0 lim lim 0 S_ 0 Dt!0 Dt Dt!0 dt Dt is discretized for a time increment Dt. The resulting formula ÿ 1 ÿ exp ÿ Dt=sj n1 Dt n1 n S0 ÿ Sn0 Hj Cj Hj exp ÿ sj Dt=sj
48
49
is an ecient approach for an approximate update of the anisotropic viscoelastic stress variables. The quantities Hnj and Sn0 of time tn have to be supplied for the computation of the present stresses. Relaxation time sj and the fourth order relaxation tensor Cj are material parameters to be determined. The constant tensor Cj A A Aj C 1; A AAÿ1 0 C 1; A A const:
50
is evaluated for the underformed material, i.e., C 1. The constitutive dierential equation (46) de®nes a linear relation between Hj and S_ 0 by the relaxation tensor Cj independent of the deformation. Nonlinearity is introduced by the ®nite strain measure C and the application of nonlinear elastic characteristics using an appropriate potential function W W C; A A. This way, the application of the relaxation tensor allows a separation of linear and nonlinear parts of the model. Moreover, this elastic and viscous split is also important for parameter identi®cation. The dierentiation of the algorithmic stress tensor with respect to strain yields the consistent tangent operator ( ÿ ) n X 1 ÿ exp ÿ Dt=sj oSn1 n1 Cj A A 51 An1 A : 2 n1 I 0 C; A A; Dt=s oC j j1 which is determined by the constant tensorial relaxation function n X t Cj exp ÿ C t I sj j1
52
and the nonlinear elastic material tensor A0n1 . For a time step Dt the algorithmic tensor (51) is divided up into a constant viscous term fg and the nonlinear elastic contribution An1 0 . This fourth order tensor A C 1; A A is symmetric for the undeformed material. In case of loading, i.e., C 6 1, an unsymmetric quantity A C; A A is obtained. Therefore, the entire ®nite element matrices of the discrete problem have to be evaluated and assembled. The implementation of (51) in a ®nite element code ensures a quadratic rate of convergence of the Newton procedure at least near the ®nal solution. The tangential material tensor in the current con®guration a is calculated by a push forward operation a U A
53
utilizing the deformation gradient. A direct analytical determination of the spatial tensor a is not possible in contrast to the isotropic formulation, where the tensorial relaxation function C t reduces to a scalar function. 4. Numerical examples 4.1. Dynamic response of U-pro®le The ®rst application of the introduced models simulates the geometrically linear dynamic response of a composite pro®le. The geometry (U 200 60 10 mm3 ) is taken from the design manual of Fiberline Composites A/S [3]. Glass ®bres, which make up 50% of the material's volume, are arranged in x-direction and combined with a polymer matrix. Consequently, material characteristics in y- and z-orientation are
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
235
Fig. 2. Discretization and deformed con®guration of U-pro®le.
Fig. 3. Dynamic response of U-pro®le.
isotropic. The pro®le l 600 mm is represented by 380 linear brick elements (see Fig. 2). All nodes of the rear cross-section are ®xed and the tip of this cantilever beam is subjected to a distributed load of q 125 N=mm constant with regard to time. Fig. 3 shows the horizontal and the vertical response of marked point P as function of time. Firstly, the pro®le is simulated using a viscoelastic composite material and, secondly, compared to an isotropic steel construction. In the latter case, dissipative behaviour of the structure is not given by the constitutive
236
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
Fig. 4. Rayleigh-damping of steel.
relation but de®ned utilizing Rayleigh-damping. The two damping parameters are taken to yield a nearly frequency-independent dissipation for high modes and a structural damping ratio of n 0:001273 (Fig. 4) according to Petersen [16] for the dominant torsional mode at xsteel 1314 rad=s. The comparison of the two types of pro®le point out signi®cant dierences. The steel pro®le is less deformed staticly and dynamically because of a larger Young's modulus. The vibrations of the composite cantilever are strongly damped owing to the viscoelastic polymer resin. The lower stiness of the composite material results in an oscillation at a reduced main frequency xcomposite 752 rad=s compared to the steel pro®le xsteel 1314 rad=s. The dominant eigenmode is a torsional deformation (Fig. 2) because the resultant of the exciting load is not located at the shear center. 4.2. Elastic analysis of a human cornea In the introduction it was pointed out that a lot of biological tissues are characterized by a ®brous microstructure and, therefore, anisotropic constitutive formulations apply to it. Parts of the human eye, investigated in this numerical example, fall into this class of materials. Fig. 5 depicts a partial discretization with eight noded 2304 mixed Q1/P0-brick elements (according to [22]). The inner part (cornea) is represented by a ®ne mesh and the outer part (sclera) is given by a coarse discretization. The state of deformation of the cornea is of primary interest in this simulation and the sclera is merely considered as elastic boundary condition for the cornea.
Fig. 5. Finite element discretization of the human eye.
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
237
The cornea, which is very important part of the optical system, is built of some hundreds of transparent collagen ®bre layers and it has a total thickness of 520 lm. One collagen ®bre is relatively sti in ®bre direction and very soft perpendicular to it. The ®bre layers are arranged with arbitrary direction in the cornea plane. The superposition of collagen layers with dierent orientation yields a stier isotropic plane and a reduced stiness rectangular to it. In order to obtain an ecient description, a transverse isotropic potential function W U vol J W iso C W aniso C; A A
54
is introduced. The potential is split up into a volumetric and into an isochoric contribution to model quasiincompressibility. The isochoric part consists of a small isotropic basic stiness W iso which is superimposed by an isochoric anisotropic part W aniso
m X n 1 X i1 C a IVj ÿ 1 : m j1 i1 i
55
The ®bre stiness (55) is de®ned as an expansion of the fourth invariant IVj IV C; A aj
56
using n terms, where aj j ÿ 1 p=m for j 1; . . . ; m. The ®bre is rotated in the cornea plane in m steps and superimposed, i.e., the mean stiffness for m ®bre angles aj is computed on the material scale. A ®rst simulation investigates the state of deformation of the cornea at inner pressure pi . Fig. 6 depicts a cross-section of the undeformed ideal cornea and the state of deformation at pressure pi 0:007 MPa (deformation plot is scaled 2.5 times). The accompanying radius of curvature, which is a very important parameter in opthalmology, is shown in Fig. 8. The deformed con®guration (Fig. 6) as well as the plot of the radius of curvature (Fig. 8) clearly show a ¯attening of the cornea. The increase of the radius of curvature, which is constant at r 7:86 mm for the undeformed structure, is larger at the outer part of the cornea compared to the center. The postprocessing algorithm evaluating the curvature from the state of deformation yields a result that is erroneously not axisymmetric, as can be seen from Fig. 8. The algorithm is very sensitive and we compute a result in¯uenced by the ®nite element discretization, which also is not axisymmetric. A second ®nite element simulation (Fig. 7) analyses the in¯uence of a single stitch on the deformation behaviour, i.e., on the radius of curvature (Fig. 8). The nearly rigid stitch is discretized by one-dimensional
Fig. 6. Cross-section of the undeformed and the deformed con®guration of the cornea.
Fig. 7. Position of a single stitch.
238
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
Fig. 8. Radius of curvature (ideal cornea subjected to internal pressure).
Fig. 9. Radius of curvature (cornea with a single stitch subjected to internal pressure).
truss elements implemented at the place where a stitch is located if a new cornea was implanted by eye surgery. Fig. 9 shows the resulting curvature at pressure pi 0:007 MPa. The diagram depicts the local character of the pertubation. About three quarters of the cornea do not show any in¯uence of the stitch. The ®nite element results computed in this study compare qualitatively well with experimental investigations to be published in the future (see [31]).
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
239
4.3. Viscoelastic simulation of a two layer laminate A numerical study investigating a two layer tensile test is carried out in order to illustrate on the one hand the eect of an unsymmetric stacking order on the state of deformation and on the other hand the ®nite viscoelastic model. Geometry and boundary conditions for the two samples (®bre angle: a 90 = 0 , a 45 = ÿ 45 ) are identical. The specimen 400 100 16 mm3 is ®xed at one crosssection and loaded displacement controlled at the other end. The simulation is based on polyester ®bres and elastomeric material. Load-strain characteristics of the material components are given by Fig. 10. The diagram is evaluated for a single ®bre which can be loaded up to e 16:5%. The second graph (Fig. 10) depicts a typical nonlinear material curve of rubber: softening at the beginning followed by stiening. The composite material is de®ned by a constitutive transversely isotropic potential according to Eq. (54). The volumetric term U vol
j 2
J2 ÿ 1 ÿ ln J 2
57
is associated with the quasi-incompressible volume deformation J 1. The characteristics of the isochoric part are determined by the isotropic function 2
W iso C1 I ÿ 3 C2 I ÿ 3 C3 I ÿ 3
3
58
introduced by Yeoh [32], which is a realistic phenomenological model for rubber elasticity (see [9]). The expansion W aniso
n X i1
Cia IV ÿ 1
i1
59
in n terms of the fourth invariant IV describes the nonlinear stiffness in ®bre direction. The numerical model of the specimen (Fig. 11) consists of a discretization of 459 Q1/P0-brick elements. Fig. 11 shows the deformed con®gurations at an elongation of 66 mme16:5%. For layers a 90 = 0 the state of deformation is scaled three times while for a 45 = ÿ 45 the deformation plot is not enlarged. The deformed specimen shows either a bending-extension coupling a 90 = 0 or a torsion±extension coupling a 45 = ÿ 45 . The latter state of deformation results in a twisting of about 90 . In a displacement controlled numerical cyclic test the reaction force F is computed as function of the displacement u at the cross-section subjected to deformation (Fig. 12). The ®nite transversely isotropic viscoelastic model is employed. The load±displacement curves form a hysteresis because of the application of the inelastic material model. The area described by the curves stands for the dissipated energy per loading cycle. The difference between the viscoelastic stiffness of the two ®bre con®gurations can be seen clearly.
Fig. 10. Experimental investigation of polyester ®bre and rubber.
240
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
Fig. 11. Deformed con®gurations of two layer laminate (e 16:5%).
Fig. 12. Cyclic extension test.
5. Conclusions This paper has presented constitutive formulations for material that is idealized as transversely isotropic. The approaches cover material with ®brous microstructure like biological tissues and all kinds of undirectional ®bre matrix combinations. Elasticity and viscoelasticity are considered as two main aspects of the material characteristics. Depending on application and composite components, it has to be decided whether in®nitesimal or ®nite strains apply. The transversely isotropic concepts shown can be extended to fully orthotropic description if necessary. A reduction of the stiness depending on loading, summarized as `damage phenomena', is not considered in the present paper. Future work will add the micromechanical change of the material structure due to ®bre or matrix rupture and molecular debonding like Mullins-eect for rubber material. An anisotropic damage model is capable of predicting this behaviour. The present article focuses on the homogeneous macroscopic material level in order to provide a formulation suitable for structural ®nite element simulations. The ®eld of parameter identi®cation and determination of material characteristics from the material's heterogenous micro structure is excluded and will be discussed in a second study. In principle, the mean properties can be evaluated analytically, numerically and experimentally. For small strains, very ecient homogenization procedures either in the timedomain, or in the Laplace-domain, depending on elastic or viscoelastic characterization, have been derived and veri®ed by ®nite element simulations on the micro scale. At ®nite strains homogenization is quite dicult and still an area of research.
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
241
Appendix A A.1. Tensor invariants Based on an arbitrary second order tensor B and a unit vector A ®ve irreducible invariants IB tr B;
A:1
1 IIB tr2 B ÿ tr B2 ; 2
A:2
IIIB det B;
A:3
IVB A : B : A;
A:4
VB A : B2 : A
A:5
are formulated, which can be used for example to de®ne a transversely isotropic strain energy function (Eq. (26)). In order to develop the complete material formulation, partial derivatives oIB 1; oB
A:6
oIIB IB 1 ÿ B; oB
A:7
oIIIB IIIB Bÿ1 ; oB
A:8
oIVB A A; oB
A:9
oVB A BA BA A oB
A:10
of the invariants with respect to tensor B are needed (see Section 3.1). The speci®c tensors shown below are given with respect to a Cartesian basis in symbolic form and they are also depicted with its components. These quantities are employed to compute the tangential material tensor. The derivative of tensor B with respect to itself yields a fourth order unit tensor oB 1 I Iijkl ei ej ek el dik djl dil djk ei ej ek el ; oB 2
A:11
where dij is the Kronecker symbol. The derivative of the inverse tensor Bÿ1 with respect to B oBÿ1 1 ÿ1 ÿIBÿ1 ÿIBÿ1 ;ijkl ei ej ek el ÿ Bÿ1 Bÿ1 Bÿ1 il Bjk ei ej ek el ; 2 ik jl oB
A:12
results in a fourth order tensor equivalent to the unit tensor but utilizing the components of Bÿ1 . Computing the derivative of M A B A B A A with respect to tensor B, i.e., evaluating the second derivative of the ®fth invariant with respect to B, leads to a constant fourth order tensor
242
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
oM 1 IA IA;ijkl ei ej ek el Ak Ai djl Ai Al djk Aj Ak dil Aj Al dik ei ej ek el ; oB 2
A:13
assembled from the components of vector A. A.2. Transformation operator A transformation from Cartesian material coordinates ~xi the transformation operator 2 2 m21 n21 l1 m 1 l1 n1 l1 6 l2 2 2 m n l m l2 n2 6 2 2 2 2 2 6 2 2 2 6 l3 m3 n3 l3 m 3 l3 n3 T6 6 2l l 2m m 2n n l m l m l n l n 1 2 1 2 1 2 2 1 1 2 2 1 6 1 2 6 4 2l1 l3 2m1 m3 2n1 n3 l1 m3 l3 m1 l1 n3 l3 n1 2l2 l3 2m2 m3 2n2 n3 l2 m3 l3 m2 l2 n3 l3 n2
to global coordinates xi is accomplished with 3 m 1 n1 7 m 2 n2 7 7 7 m 3 n3 7; m1 n2 m2 n1 7 7 7 m1 n3 m3 n1 5 m2 n3 m3 n2
A:14
where li ; mi ; ni are the direction cosines.
References [1] R.M. Christensen, Mechanics of Composite Materials, Wiley, New York, 1979. [2] T.C. Doyle, J.L. Ericksen, Nonlinear elasticity, in: Advances in Applied Mechanics IV, Academic Press, New York, 1956. [3] Fiberline, Design- und Konstruktionshandbuch f ur Konstruktionspro®le aus Verbundwerkstoen, Fiberline Composites A/S, Kolding, 1995. [4] R.F. Gibson, Principles of Composite Material Mechanics, McGraw-Hill, New York, 1994. [5] S. Govindjee, J.C. Simo, Mullins' eect and the strain amplitude dependence of the storage modulus, International Journal of Solids and Structures 29 (1992) 1737±1751. [6] L.R. Herrmann, F.E. Peterson, A numerical procedure for viscoelastic stress analysis, in: Proceedings of the Seventh Meeting of ICRPG Mechanical Behaviour Working Group, Orlando, 1968. [7] R.M. Jones, Mechanics of Composite Materials, McGraw-Hill, Tokyo, 1975. [8] M. Kaliske, H. Rothert, Formulation and implementation of three-dimensional viscoelasticity at small and ®nite strains, Computational Mechanics 19 (1997) 228±239. [9] M. Kaliske, H. Rothert, On the ®nite element implementation of rubberlike materials at ®nite strains, Engineering Computations 14 (1997) 217±233. [10] Y.C. Lou, R.A. Schapery, Viscoelastic characterization of a nonlinear ®ber reinforced plastic, Journal of Composite Materials 5 (1971) 208±234. [11] J. Lubliner, A model of rubber viscoelasticity, Mechanics Research Communications 12 (1985) 93±99. [12] T.C. Kennedy, M. Wang, Three-dimensional, nonlinear viscoelastic analysis of laminated composites, Journal of Composite Materials 28 (1994) 902±925. [13] K.-D. Klee, E. Stein, Discretization in nonlinear continuum mechanics and comparison of some ®nite element algorithms, in: Proceedings of the Pressure Vessels and Piping Conference, Orlando, 1982. [14] K.Y. Lin, I.H. Hwang, Thermo-viscoelastic analysis of composite materials, Journal of Composite Materials 23 (1989) 554±569. [15] A.I. Lurie, Nonlinear theory of elasticity, in: J.D. Achenbach et al. (Eds.), North-Holland Series in Applied Mathematics and Mechanics, North-Holland, Amsterdam, 1990. [16] C. Petersen, Dynamik der Baukonstruktionen, Vieweg, Braunschweig, 1996. [17] H. Poon, F.M. Ahmad, A material point time integration procedure for anisotropic thermo rheologically simple, viscoelastic solids, Computational Mechanics 21 (1998) 236±242. [18] J. Schr oder, Theoretische und algorithmische Konzepte zur ph anomenologischen Beschreibung anisotropen Materialverhaltens, Dissertation, Universit at Hannover, 1996. [19] F.R. Schwarzl, A.J. Staverman, Time±temperature dependence of linear viscoelastic behaviour, Journal of Applied Physics 23 (1952) 838±843. [20] F. Sidoro, Un modele viskoelastique non lineaire avec con®guration intermediaire, Journal de Mecanique 13 (1974) 679±713. [21] J.C. Simo, On a fully three-dimensional ®nite-strain viscoelastic damage model: formulation and computational aspects, Computer Methods in Applied Mechanics and Engineering 60 (1987) 153±173. [22] J.C. Simo, R.L. Taylor, Quasi-incompressible ®nite elasticity in principal stretches, continuum basis and numerical algorithms, Computer Methods in Applied Mechanics and Engineering 85 (1991) 273±310.
M. Kaliske / Comput. Methods Appl. Mech. Engrg. 185 (2000) 225±243
243
[23] F. Soergel, Biomechanische Charakterisierung der menschlichen Augenhornhaut mit dynamisch-mechanischer Spektroskopie, Dissertation, Universit at Ulm, 1994. [24] A.J.M. Spencer, Continuum theory of the mechanics of ®bre-reinforced composites, CISM Courses and Lectures, No. 282, Springer, Wien, 1984. [25] E.S. Suhubi, Thermoelastic solids, in: A.C. Eringen (Ed.), Continuum Physics II, Academic Press, New York. [26] P. Le Tallec, C. Rahier, A. Kaiss, Three-dimensional incompressible viscoelasticity in large strains: formulation and numerical approximation, Computer Methods in Applied Mechanics and Engineering 109 (1993) 233±258. [27] R.L. Taylor, K.S. Pister, G.L. Goudreau, Thermomechanical analysis of viscoelastic solids, International Journal for Numerical Methods in Engineering 2 (1970) 45±59. [28] C. Truesdell, W. Noll, The nonlinear ®eld theories, in: S. Fl ugge (Ed.), Handbuch der Physik, Band III/3, Springer, Berlin, 1965. [29] M. Wang, Three-dimensional, nonlinear viscoelastic analysis of laminated composites: a ®nite element approach, Dissertation, Oregon State University at Corvallis, 1993. [30] J.A. Weiss, N.M. Bradley, S. Govindjee, Finite element implementation of incompressible, transversely isotropic hyperelasticity, Computer Methods in Applied Mechanics and Engineering 135 (1996) 107±128. [31] T. Wolf, Druckinduzierte Anderungen der Hornhautkr ummungen bei Schweinebulbie mit Einzelknopfnaht und fortlaufender Keratoplastiknaht, Dissertation, Medizinische Hochschule Hannover (to be published). [32] O.H. Yeoh, Characterization of elastic properties of carbon-black-®lled rubber vulcanizates, Rubber Chemistry and Technology 63 (1990) 792±805. [33] S. Yi, Finite element analysis of anisotropic viscoelastic composite structures and analytical determination of optimum viscoelastic material properties, Dissertation, University of Illinois at Urbana-Champaign, 1992. [34] S. Yi, H. Hilton, Dynamic ®nite element analysis of viscoelastic composite plates in the time domain, International Journal for Numerical Methods in Engineering 37 (1994) 4081±4096. [35] S. Yi, H. Hilton, Hygrothermal eects on viscoelastic responses of laminated composites, Composites Engineering 5 (1995) 183± 193. [36] M.A. Zocher, S.E. Groves, A three-dimensional ®nite element formulation for thermoviscoelastic orthotropic media, International Journal for Numerical Methods in Engineering 40 (1997) 2267±2288. [37] J.C. Simo, T.J.R. Hughes, Computational Inelasticity, Interdisciplinary Applied Mathematics 7, Springer, New York, 1998.