CHAPTER Elasticity and Viscoelasticity 2 CHAPTER 2.1 Introduction to Elasticity and Viscoelasticity JEAN LEMAITRE Universit!e Paris 6, LMT-Cachan, 61 avenue du Pr!esident Wilson, 94235 Cachan Cedex, France For all solid materials there is a domain in stress space in which strains are reversible due to small relative movements of atoms. For many materials like metals, ceramics, concrete, wood and polymers, in a small range of strains, the hypotheses of isotropy and linearity are good enough for many engineering purposes. Then the classical Hooke’s law of elasticity applies. It can be derived from a quadratic form of the state potential, depending on two parameters characteristics of each material: the Young’s modulus E and the Poisson’s ratio n. 1 AijklðE;nÞ sij skl c* ¼ ð1Þ 2r eij ¼ r @c * 1 þ n n sij skk dij ¼ E E @sij ð2Þ E and n are identified from tensile tests either in statics or dynamics. A great deal of accuracy is needed in the measurement of the longitudinal and transverse strains (de 106 in absolute value). When structural calculations are performed under the approximation of plane stress (thin sheets) or plane strain (thick sheets), it is convenient to write these conditions in the constitutive equation. Plane stress ðs33 ¼ s13 ¼ s23 ¼ 0Þ: 3 1 n 07 E 2 3 6 3 72 6E e11 7 s11 6 7 6 1 6 6 7 7 4 e22 5 ¼ 6 07 74 s22 5 6 E 7 6 e12 7 s12 6 5 4 1þn Sym E 2 Handbook of Materials Behavior Models Copyright # 2001 by Academic Press. All rights of reproduction in any form reserved. ð3Þ 71 72 Lemaitre Plane strain ðe33 ¼ e13 ¼ e23 ¼ 0Þ: 2 3 2 l þ 2m s11 6 7 6 4 s22 5 ¼ 4 s12 l l þ 2m Sym with 8 > > >l ¼ < > > > : 32 3 0 e11 7 76 0 54 e22 5 2m e12 nE ð1 þ nÞð1 2nÞ ð4Þ E m¼ 2ð1 þ nÞ For orthotropic materials having three planes of symmetry, nine independent parameters are needed: three tension moduli E1 ; E2 ; E3 in the orthotropic directions, three shear moduli G12 ; G23 ; G31 , and three contraction ratios n12 ; n23 ; n31 . In the frame of orthotropy: 3 3 2 2 e11 3 s11 2 1 n n 12 13 7 7 6 6 0 0 0 7 6 E 7 6 76 E1 E1 7 6 1 7 6 6 76 7 6 e22 7 6 s 7 22 7 1 n23 7 6 6 6 0 0 0 7 7 6 7 6 6 7 7 6 7 6 E2 E2 76 7 6 7 6 6 76 7 6 e 7 6 1 7 s 33 7 6 33 7 6 6 0 0 0 76 7 6 7 6 E 76 3 7¼6 7 6 ð5Þ 76 7 6 7 6 1 76 7 6 7 6 7 0 0 6 e23 7 6 s23 7 76 2G23 7 6 7 6 76 7 6 7 6 6 7 1 7 6 7 6 6 0 7 7 6 Sym 7 6 6 7 6 e31 7 6 6 2G31 76 s31 7 7 4 7 6 5 7 7 6 6 1 5 5 4 4 2G12 e12 s12 Nonlinear elasticity in large deformations is described in Section 2.2, with applications for porous materials in Section 2.3 and for elastomers in Section 2.4. Thermoelasticity takes into account the stresses and strains induced by thermal expansion with dilatation coefficient a. For small variations of temperature y for which the elasticity parameters may be considered as constant: eij ¼ 1þn n sij skk dij þ aydij E E ð6Þ For large variations of temperature, E; n; and a will vary. In rate formulations, such as are needed in elastoviscoplasticity, for example, the 2.1 Introduction to Elasticity and Viscoelasticity 73 derivative of E; n; and a must be considered. 1þn n @ 1þn @ n @a ’ s’ ij s’ kk dij þ aydij þ sij skk dij þ ydij y’ e’ij ¼ E E @y E @y E @y ð7Þ Viscoelasticity considers in addition a dissipative phenomenon due to ‘‘internal friction,’’ such as between molecules in polymers or between cells in wood. Here again, isotropy, linearity, and small strains allow for simple models. Quadratric functions for the state potential and the dissipative potential lead to either Kelvin-Voigt or Maxwell’s models, depending upon the partition of stress or strains in a reversible part and in an irreversible part. They are described in detail for the one-dimensional case in Section 2.5 and recalled here in three dimensions. Kelvin-Voigt model: sij ¼ lðekk þ yl e’ kk Þdij þ 2mðeij þ ym e’ij Þ ð8Þ Here l and m are Lame’s coefficients at steady state, and yl and ym are two time parameters responsible for viscosity. These four coefficients may be identified from creep tests in tension and shear. Maxwell model: 1þn s n skk ð9Þ s’ ij þ s’ kk þ dij e’ij ¼ E t1 E t2 Here E and n are Young’s modulus and Poisson’s ratio at steady state, and t1 and t2 are two other time parameters. It is a fluidlike model: equilibrium at constant stress does not exist. In fact, a more general way to write linear viscoelastic constitutive models is through the functional formulation by the convolution product as any linear system. The hereditary integral is described in detail for the one-dimensional case, together with its use by the Laplace transform, in Section 2.5. Z t n X dskl p eijðtÞ ¼ dt þ Jijkl ðt tÞ Jijkl ðt tÞDskl ð10Þ dt o p¼1 p JðtÞ is the creep functions matrix, and Dskl are the eventual stress steps. The dual formulation introduces the relaxation functions matrix RðtÞ Z t n X dekl p dt þ Rijkl ðt tÞ Rijkl Dekl sijðtÞ ¼ ð11Þ dt o p¼1 When isotropy is considered the matrix, ½ J and ½R each reduce to two functions: either JðtÞ, the creep function in tension, is identified from a creep 74 Lemaitre test at constant stress; JðtÞ ¼ eðtÞ =s and K, the second function, from the creep function in shear. This leads to Dsij Dskk ð12Þ eij ¼ ð J þ KÞ K dij Dt Dt where stands for the convolution product and D for the distribution derivative, taking into account the stress steps. Or MðtÞ, the relaxation function in shear, and LðtÞ , a function deduced from M and from a relaxation test in tension RðtÞ ¼ sðtÞ =e; LðtÞ ¼ MðR 2MÞ=ð3M RÞ Deij Dðekk Þ ð13Þ dij þ 2M sij ¼ L Dt dt All of this is for linear behavior. A nonlinear model is described in Section 2.6, and interaction with damage is described in Section 2.7. CHAPTER 2.2 Background on Nonlinear Elasticity R. W. OGDEN Department of Mathematics, University of Glasgow, Glasgow G12 8QW, UK Contents 2.2.1 Validity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 2.2.2 Deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 2.2.3 Stress and Equilibrium . . . . . . . . . . . . . . . . . . . . . . 77 2.2.4 Elasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 2.2.5 Material Symmetry . . . . . . . . . . . . . . . . . . . . . . . . . . 78 2.2.6 Constrained Materials . . . . . . . . . . . . . . . . . . . . . . . 80 2.2.7 Boundary-Value Problems . . . . . . . . . . . . . . . . . . . . 82 References. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 2.2.1 VALIDITY The theory is applicable to materials, such as rubberlike solids and certain soft biological tissues, which are capable of undergoing large elastic deformations. More details of the theory and its applications can be found in Beatty  and Ogden . 2.2.2 DEFORMATION For a continuous body, a reference configuration, denoted by Br, is identified and @Br denotes the boundary of Br . Points in Br are labeled by their position vectors X relative to some origin. The body is deformed quasistatically from Br so that it occupies a new configuration, denoted B, with Handbook of Materials Behavior Models Copyright # 2001 by Academic Press. All rights of reproduction in any form reserved. 75 76 Ogden boundary @B. This is the current or deformed configuration of the body. The deformation is represented by the mapping w : Br ! B, so that x ¼ vðXÞ X 2 Br ; ð1Þ where x is the position vector of the point X in B. The mapping v is called the deformation from Br to B, and v is required to be one-to-one and to satisfy appropriate regularity conditions. For simplicity, we consider only Cartesian coordinate systems and let X and x, respectively, have coordinates Xa and xi , where a; i 2 f1; 2; 3g, so that xi ¼ wi ðXa Þ. Greek and Roman indices refer, respectively, to Br and B, and the usual summation convention for repeated indices is used. The deformation gradient tensor, denoted F, is given by F ¼ Grad x Fia ¼ @xi [email protected] ð2Þ Grad being the gradient operator in Br . Local invertibility of v and its inverse requires that 05 J det F51 ð3Þ wherein the notation J is defined. The deformation gradient has the (unique) polar decompositions F ¼ RU ¼ VR ð4Þ where R is a proper orthogonal tensor and U, V are positive definite and symmetric tensors. Respectively, U and V are called the right and left stretch tensors. They may be put in the spectral forms 3 3 X X ð5Þ U¼ li uðiÞ uðiÞ V¼ li vðiÞ vðiÞ i¼1 i¼1 where vðiÞ ¼ RuðiÞ ; i 2 f1; 2; 3g, li are the principal stretches, uðiÞ the unit eigenvectors of U (the Lagrangian principal axes), vðiÞ those of V (the Eulerian principal axes), and denotes the tensor product. It follows from Eq. 3 that J ¼ l1 l2 l3 : The right and left Cauchy-Green deformation tensors, denoted C and B, respectively, are defined by C ¼ FT F ¼ U2 B ¼ FFT ¼ V2 ð6Þ 2.2.3 STRESS AND EQUILIBRIUM Let rr and r be the mass densities in Br and B, respectively. The mass conservation equation has the form rr ¼ rJ ð7Þ 77 2.2 Background on Nonlinear Elasticity The Cauchy stress tensor, denoted r, and the nominal stress tensor, denoted S, are related by S ¼ JF1 r ð8Þ The equation of equilibrium may be written in the equivalent forms div r þ rb ¼ 0 Div S þ rr b ¼ 0 ð9Þ where div and Div denote the divergence operators in B and Br , respectively, and b denotes the body force per unit mass. In components, the second equation in Eq. 9 is @Sai þ rr bi ¼ 0 @Xa ð10Þ Balance of the moments of the forces acting on the body yields simply rT ¼ r, equivalently ST FT ¼ FS: The Lagrangian formulation based on the use of S and Eq. 10, with X as the independent variable, is used henceforth. 2.2.4 ELASTICITY The constitutive equation of an elastic material is given in the equivalent forms S ¼ HðFÞ ¼ @W ðFÞ @F r ¼ GðFÞ J1 FHðFÞ ð11Þ where H is a tensor-valued function, defined on the space of deformation gradients F, W is a scalar function of F and the symmetric tensor-valued function G is defined by the latter equation in Eq. 11. In general, the form of H depends on the choice of reference configuration and it is referred to as the response function of the material relative to Br associated with S. For a given Br , therefore, the stress in B at a (material) point X depends only on the deformation gradient at X. A material whose constitutive law has the form of Eq. 11 is generally referred to as a hyperelastic material and W is called a strain-energy function (or stored-energy function). In components, (11)1 has the form Sai ¼ @[email protected] , which provides the convention for ordering of the indices in the partial derivative with respect to F. If W and the stress vanish in Br , so that WðIÞ ¼ 0 @W ðIÞ ¼ O @F ð12Þ where I is the identity and O the zero tensor, then Br is called a natural configuration. 78 Ogden Suppose that a rigid-body deformation x * ¼ Qx þ c is superimposed on the deformation x ¼ vðXÞ, where Q and c are constants, Q being a rotation tensor and c a translation vector. The resulting deformation gradient, F * say, is given by F * ¼ QF: The elastic stored energy is required to be independent of superimposed rigid deformations, and it follows that WðQFÞ ¼ WðFÞ ð13Þ for all rotations Q. A strain-energy function satisfying this requirement is said to be objective. Use of the polar decomposition (Eq. 4) and the choice Q ¼ RT in Eq. 13 shows that WðFÞ ¼ WðUÞ: Thus, W depends on F only through the stretch tensor U and may therefore be defined on the class of positive definite symmetric tensors. We write @W ð14Þ T¼ @U for the (symmetric) Biot stress tensor, which is related to S by T ¼ ðSR þ RT ST Þ=2. 2.2.5 MATERIAL SYMMETRY Let F and F0 be the deformation gradients in B relative to two different reference configurations, Br and B0r respectively. In general, the response of the material relative to B0r differs from that relative to Br , and we denote by W and W 0 the strain-energy functions relative to Br and B0r . Now let P ¼ Grad X0 be the deformation gradient of B0r relative to Br , where X0 is the position vector of a point in B0r . Then F ¼ F0 P: For specific P we may have W 0 ¼ W, and then WðF0 PÞ ¼ WðF0 Þ ð15Þ for all deformation gradients F0 . The set of tensors P for which Eq. 15 holds forms a multiplicative group, called the symmetry group of the material relative to Br . This group characterizes the physical symmetry properties of the material. For isotropic elastic materials, for which the symmetry group is the proper orthogonal group, we have WðFQÞ ¼ WðFÞ ð16Þ for all rotations Q. Since the Q’s appearing in Eqs. 13 and 16 are independent, the combination of these two equations yields WðQUQT Þ ¼ WðUÞ ð17Þ 79 2.2 Background on Nonlinear Elasticity for all rotations Q. Equation 17 states that W is an isotropic function of U. It follows from the spectral decomposition (Eq. 5) that W depends on U only through the principal stretches l1 ; l2 , and l3 and is symmetric in these stretches. For an isotropic elastic material, r is coaxial with V, and we may write r ¼ a0 I þ a1 B þ a2 B2 ð18Þ where a0 ; a1 , and a2 are scalar invariants of B (and hence of V) given by @W 1=2 @W 1=2 @W 1=2 @W a1 ¼ 2I3 þ I1 ð19Þ a0 ¼ 2I3 a2 ¼ 2I3 @I3 @I1 @I2 @I2 and W is now regarded as a function of I1 ; I2 , and I3 , the principal invariants of B defined by I1 ¼ trðBÞ ¼ l21 þ l22 þ l23 ; ð20Þ I2 ¼ 12 ½I21 tr ðB2 Þ ¼ l22 l23 þ l23 l21 þ l21 l22 ð21Þ I3 ¼ det B ¼ l21 l22 l23 ð22Þ Another consequence of isotropy is that S and r have the decompositions S¼ 3 X ti uðiÞ vðiÞ i¼1 r¼ 3 X si vðiÞ vðiÞ ð23Þ i¼1 where si ; i 2 f1; 2; 3g are the principal Cauchy stresses and ti the principal Biot stresses, connected by ti ¼ @W ¼ Jl1 i si @li ð24Þ Let the unit vector M be a preferred direction in the reference configuration of the material, i.e., a direction for which the material response is indifferent to arbitrary rotations about the direction and to replacement of M by M. Such a material can be characterized by a strain energy which depends on F and the tensor M M [2, 4, 5] Thus, we write WðF; M MÞ. The required symmetry (transverse isotropy) reduces W to dependence on the five invariants I1 ; I2 ; I3 ; I4 ¼ M ðCMÞ I5 ¼ M ðC2 MÞ ð25Þ where I1 ; I2 ; and I3 are defined in Eqs. (20)–(22). The resulting nominal stress tensor is given by S ¼ 2W1 FT þ 2W2 ðI1 I CÞFT þ 2I3 W3 F1 þ 2W4 M FM þ 2W5 ðM FCM þ CM FMÞ where Wi ¼ @[email protected] ; i ¼ 1; . . . ; 5. ð26Þ 80 Ogden When there are two families of fibers corresponding to two preferred directions in the reference configuration, M and M0 say, then, in addition to Eq. 25, the strain energy depends on the invariants I6 ¼ M0 ðCM0 Þ I7 ¼ M0 ðC2 M0 Þ I8 ¼ M ðCM0 Þ ð27Þ 0 and also on M M (which does not depend on the deformation); see Spencer [4, 5] for details. The nominal stress tensor can be calculated in a similar way to Eq. 26. 2.2.6 CONSTRAINED MATERIALS An internal constraint, given in the form CðFÞ ¼ 0, must be satisfied for all possible deformation gradients F, where C is a scalar function. Two commonly used constraints are incompressibility and inextensibility, for which, respectively, CðFÞ ¼ detF 1 CðFÞ ¼ M ðFT FMÞ 1 ð28Þ where the unit vector M is the direction of inextensibility in Br . Since any constraint is unaffected by a superimposed rigid deformation, C must be an objective scalar function, so that CðQFÞ ¼ CðFÞ for all rotations Q. Any stress normal to the hypersurface CðFÞ ¼ 0 in the (nine-dimensional) space of deformation gradients does no work in any (virtual) incremental deformation compatible with the constraint. The stress is therefore determined by the constitutive law (11)1 only to within an additive contribution parallel to the normal. Thus, for a constrained material, the stress-deformation relation (11)1 is replaced by S ¼ HðFÞ þ q @C @W @C ¼ þq @F @F @F ð29Þ where q is an arbitrary (Lagrange) multiplier. The term in q is referred to as the constraint stress since it arises from the constraint and is not otherwise derivable from the material properties. For incompressibility and inextensibility we have @W @W ð30Þ S¼ þ qF1 þ 2qM FM S¼ @F @F respectively. For an incompressible material the Biot and Cauchy stresses are given by @W ð31Þ pU1 T¼ det U ¼ 1 @U 2.2 Background on Nonlinear Elasticity 81 and @W ð32Þ pI det F ¼ 1 @F where q has been replaced by p, which is called an arbitrary hydrostatic pressure. The term in a0 in Eq. 18 is absorbed into p, and I3 ¼ 1 in the remaining terms in Eq. 18. For an incompressible isotropic material the principal components of Eqs. 31 and 32 yield @W @W ð33Þ ti ¼ pl1 s i ¼ li p i @li @li r¼F respectively, subject to l1 l2 l3 ¼ 1. For an incompressible transversely isotropic material with preferred direction M, the dependence on I3 is omitted and the Cauchy stress tensor is given by r ¼ pI þ 2W1 B þ 2W2 ðI1 B B2 Þ þ 2W4FM FM þ 2W5 ðFM BFM þ BFM FMÞ ð34Þ For a material with two preferred directions, M and M0 , the Cauchy stress tensor for an incompressible material is r ¼ pI þ 2W1 B þ 2W2 ðI1 B B2 Þ þ 2W4 FM FM þ 2W5 ðFM BFM þ BFM FMÞ þ 2W6 FM0 FM0 þ 2W7 ðFM0 BFM0 þ BFM0 FM0 Þ þ W8 ðFM FM0 þ FM0 FMÞ ð35Þ where the notation Wi ¼ @[email protected] now applies for i ¼ 1; 2; 4; . . . ; 8. 2.2.7 BOUNDARY-VALUE PROBLEMS The equilibrium equation (second part of Eq. 9), the stress-deformation relation (Eq. 11), and the deformation gradient (Eq. 2) coupled with Eq. 1 are combined to give @W Div F ¼ Grad x x ¼ vðXÞ X 2 Br ð36Þ þ rr b ¼ 0 @F Typical boundary conditions in nonlinear elasticity are x ¼ nðXÞ on @Bxr ð37Þ ST N ¼ sðF; XÞ on @Btr ð38Þ where n and s are specified functions, N is the unit outward normal to @Br , 82 Ogden and @Bxr and @Btr are complementary parts of @B. In general, s may depend on the deformation through F. For a dead-load traction s is independent of F. For a hydrostatic pressure boundary condition, Eq. 38 has the form s ¼ JPFT N on @Btr ð39Þ Equations 36–38 constitute the basic boundary-value problem in nonlinear elasticity. In components, the equilibrium equation in Eq. 36 is written Aaibj @ 2 xj þ rr bi ¼ 0 @Xa @Xb ð40Þ for i 2 f1; 2; 3g, where the coefficients Aaibj are defined by Aaibj ¼ Abjai ¼ @2W @Fia @Fjb ð41Þ When coupled with suitable boundary conditions, Eq. 41 forms a system of quasi-linear partial differential equations for xi ¼ wi ðXa Þ. The coefficients Aaibj are, in general, nonlinear functions of the components of the deformation gradient. For incompressible materials the corresponding equations are obtained by substituting the first part of Eq. 30 into the second part of Eq. 9 to give Aaibj @ 2 xj @p þ rr b i ¼ 0 @Xa @Xb @xi detð@xi @Xa Þ ¼ 1 ð42Þ where the coefficients are again given by Eq. 41. In order to solve a boundary-value problem, a specific form of W needs to be given. The form of W chosen will depend on the particular material considered and on mathematical requirements relating to the properties of the equations, an example of which is the strong ellipticity condition. Equations 40 are said to be strongly elliptic if the inequality Aaibj mi mj Na Nb > 0 ð43Þ holds for all nonzero vectors m and N. Note that Eq. 43 is independent of any boundary conditions. For an incompressible material, the strong ellipticity condition associated with Eq. 42 again has the form of Eq. 43, but the incompressibility constraint now imposes the restriction m ðFT NÞ ¼ 0 on m and N. REFERENCES 1. Beatty, M. F. (1987). Topics in finite elasticity: Hyperelasticity of rubber, elastomers and biological tissues } with examples. Appl. Mech. Rev. 40; 1699–1734. 2.2 Background on Nonlinear Elasticity 83 2. Holzapfel, G. A. (2000). Nonlinear Solid Mechanics. Chichester: Wiley. 3. Ogden, R. W. (1997). Non-linear Elastic Deformations. New York: Dover Publications. 4. Spencer, A. J. M. (1972). Deformations of Fibre-Reinforced Materials. Oxford: Oxford University Press. 5. Spencer, A. J. M. (1984). Constitutive theory for strongly anisotropic solids. In Continuum Theory of the Mechanics of Fibre-Reinforced Composites, CISM Courses and Lectures No. 282, pp. 1–32, Spencer, A. J. M., ed., Wien: Springer-Verlag. CHAPTER 2.3 Elasticity of Porous Materials N. D. CRISTESCU 231 Aerospace Building, University of Florida, Gainesville, Florida Contents 2.3.1 Validity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 2.3.2 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 2.3.3 Identification of the Parameters . . . . . . . . . . . . . . 85 2.3.4 Examples. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 References. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 2.3.1 VALIDITY The methods used to determine the elasticity of porous materials and/or particulate materials as geomaterials or powderlike materials are distinct from those used with, say, metals. The reason is that such materials possess pores and=or microcracks. For various stress states these may either open or closed, thus influencing the values of the elastic parameters. Also, the stress-strain curves for such materials are strongly loading-rate-dependent, starting from the smallest applied stresses, and creep (generally any time-dependent phenomena) is exhibited from the smallest applied stresses (see Fig. 2.3.1 for schist, showing three uniaxial stress-strain curves for three loading rates and a creep curve ). Thus information concerning the magnitude of the elastic parameters cannot be obtained: from the initial slope of the stress-strain curves, since these are loadingrate-dependent; by the often used ‘‘chord’’ procedure, obviously; from the unloading slopes, since significant hysteresis loops are generally present. 84 Handbook of Materials Behavior Models Copyright # 2001 by Academic Press. All rights of reproduction in any form reserved. 2.3 Elasticity of Porous Materials 85 FIGURE 2.3.1 Uniaxial stress-strain curves for schist for various loading rates, showing time influence on the entire stress-strain curves and failure (stars mark the failure points). 2.3.2 FORMULATION The elasticity of such materials can be expressed as ‘‘instantaneous response’’ by 1 1 1 ’ T’ ð1Þ ðtrTÞ1 þ D¼ 3K 2G 3 2G % where G and K are the elastic parameters that are not constant, D is the strain rate tensor, T is the stress tensor, tr( ) is the trace operator, and 1 is the unit tensor. Besides the elastic properties described by Eq. 1, some other mechanical properties can be described by additional terms to be added to Eq. 1. For isotropic geomaterials the elastic parameters are expected to depend on stress invariants and, perhaps, on some damage parameters, since during loading some pores and microcracks may close or open, thus influencing the elastic parameters. 2.3.3 IDENTIFICATION OF THE PARAMETERS The elastic parameters can be determined experimentally by two procedures. With the dynamic procedure, one is determining the travel time of the two 86 Cristescu elastic (seismic) extended longitudinal and transverse waves, which are traveling in the body. If both these waves are recorded, then the instantaneous response is of the form of Eq. 1. The elastic parameters are obtained from 4 2 2 K ¼ r vp vS ð2Þ G ¼ rv2S 3 where vS is the velocity of propagation of the shearing waves, vp the velocity of the longitudinal waves, and r the density. The static procedure takes into account that the constitutive equations for geomaterials are strongly time-dependent. Thus, in triaxial tests performed under constant confining pressure s, after loading up to a desired stress state t (octahedral shearing stress), one is keeping the stress constant for a certain time period tc [2, 3]. During this time period the rock is creeping. When the strain rates recorded during creep become small enough, one is performing an unloading–reloading cycle (see Fig. 2.3.2). From the slopes 1 1 1 1 1 1 ð3Þ þ þ 3G 9K 6G 9K of these unloading–reloading curves one can determine the elastic parameters. For each geomaterial, if the time tc is chosen so that the subsequent unloading is performed in a comparatively much shorter time interval, no significant interference between creep and unloading phenomena will take place. An example for schist is shown in Figure 2.3.3, obtained in a triaxial test with five unloading–reloading cycles. FIGURE 2.3.2 Static procedures to determine the elastic parameters from partial unloading processes preceded by short-term creep. 87 2.3 Elasticity of Porous Materials FIGURE 2.3.3 Stress-strain curves obtained in triaxial tests on shale; the unloadings follow a period of creep of several minutes. If only a partial unloading is performed (one third or even one quarter of the total stress, and sometimes even less), the unloading and reloading follow quite closely straight lines that practically coincide. If a hysteresis loop is still recorded, it means that the time tc was chosen too short. The reason for performing only a partial unloading is that the specimen is quite ‘‘thick’’ and as such the stress state in the specimen is not really uniaxial. During complete unloading, additional phenomena due to the ‘‘thickness’’ of the specimen will be involved, including, e.g., kinematic hardening in the opposite direction, etc. Similar results can be obtained if, instead of keeping the stress constant, one is keeping the axial strain constant for some time period during which the axial stress is relaxing. Afterwards, when the stress rate becomes relatively small, an unloading–reloading is applied to determine of the elastic parameters. This procedure is easy to apply mainly for particulate materials (sand, soils, etc.) when standard (Karman) three-axial testing devices are used and the elastic parameters follow from K¼ 1 Dt 3 De1 þ 2De2 G¼ 1 Dt 2 De1 De2 ð4Þ where D is the variation of stress and elastic strains during the unloading– reloading cycle. The same method is used to determine the bulk modulus K in hydrostatic tests when the formula to be used is K¼ Ds Dev ð5Þ with s the mean stress and ev the volumetric strain. Generally, K is increasing with s and reaching an asymptotic constant value when s is increasing very much and all pores and microcracks are closed 88 Cristescu under this high pressure. The variation of the elastic parameters with t is more involved: when t increases but is still under the compressibility– dilatancy boundary, the elastic parameters are increasing. For higher values, above this boundary, the elastic parameters are decreasing. Thus their variation is related to the variation of irreversible volumetric strain, which, in turn, is describing the evolution of the pores and microcracks existing in the geomaterial. That is why the compressibility–dilatancy boundary plays the role of reference configuration for the values of the elastic parameters so long as the loading path (increasing s and=or t) is in the compressibility domain, the elastic parameters are increasing, whereas if the loading path is in the dilatancy domain (increasing under constant s), the elastic parameters are decreasing. If stress is kept constant and strain is varying by creep, in the compressibility domain volumetric creep produces a closing of pores and microcracks and thus the elastic parameters increase, and vice versa in the dilatancy domain. Thus, for each value of s the maximum values of the elastic parameters are reached on the compressibility–dilatancy boundary. 2.3.4 EXAMPLES As an example, for rock salt in uniaxial stress tests, the variation of the elastic moduli G and K with the axial stress s1 is shown in Figure 2.3.4 . The variation of G and K is very similar to that of the irreversible volumetric FIGURE 2.3.4 Variation of the elastic parameters K and G and of irreversible volumetric strain in monotonic uniaxial tests. 2.3 Elasticity of Porous Materials 89 FIGURE 2.3.5 Variation in time of the elastic parameters and of irreversible volumetric strain in uniaxial creep tests. strain eIV . If stress is increased in steps, and if after each increase the stress in kept constant for two days, the elastic parameters are varying during volumetric creep, as shown in Figure 2.3.5. Here D is the ratio of the applied stress and the strength in uniaxial compression sc ¼ 17:88 MPa. Again, a similarity with the variation of eIV is quite evident. Figure 2.3.6 shows for a different kind of rock salt the variation of the elastic velocities vP and vS in true triaxial tests under confining pressure pc ¼ 5 MPa (data by Popp, Schultze, and Kern ). Again, these velocities increase in the compressibility domain, reach their maxima on the compressibility–dilatancy boundary, and then decrease in the dilatancy domain. For shale, and the conventional (Karman) triaxial tests shown in Figure 2.3.3, the values of E and G for the five unloading–reloading cycles shown are: E ¼ 9:9, 24.7, 29.0, 26.3, and 22.3 GPa, respectively, while G ¼ 4:4, 10.7, 12.1, 10.4, and 8.5 GPa. For granite, the variation of K with s is given as  8 > < K0 K1 1 s ; if s s0 s0 KðsÞ :¼ ð6Þ > : K ; if s s 0 0 with K0 ¼ 59 GPa, K1 ¼ 48 GPa, and s0 ¼ 0:344 GPa, the limit pressure when all pores are expected to be closed. 90 Cristescu FIGURE 2.3.6 The maximum of vs takes place at the compressibility–dilatancy boundary (figures and hachured strip); changes of vp and vs as a function of strain (’e ¼ 105 s1, pc ¼ 5 Mpa, T ¼ 308 C), showing that the maxima are at the onset of dilatancy (after Reference ). The same formula for alumina powder is s KðsÞ :¼ K pa exp b pa 1 ð7Þ with K1 ¼ 1 107 kPa the constant value toward which the bulk modulus tends at high pressures, a ¼ 107 , b ¼ 1:2 104 , and pa ¼ 1 kPa. Also for alumina powder we have EðsÞ :¼ E1 pa b expðdsÞ ð8Þ with E1 ¼ 7 105 kPa, b ¼ 6:95 105 , and d ¼ 0:002. For the shale shown in Figure 2.3.3, the variation of K with s for 0 s 45 MPa is KðsÞ :¼ 0:78s2 þ 65:32s 369 ð9Þ REFERENCES 1. Cristescu, N. (1986). Damage and failure of viscoplastic rock-like materials. Int. J. Plasticity 2 (2): 189–204. 2. Cristescu, N. (1989). Rock Rheology, Kluver Academic Publishing. 3. Cristescu, N. D., and Hunsche, U. (1998). Time Effects in Rock Mechanics, Wiley. 4. Ani, M., and Cristescu N. D. (2000). The effect of volumetric strain on elastic parameters for rock salt. Mechanics of Cohesive-Frictional Materials 5 (2): 113–124. 5. Popp, T., Schultze, O., and Kern, H. ( ). Permeation and development of dilatancy and permeability in rock salt, in The Mechanical Behavior of Salt (5th Conference on Mechanical Behavior of Salt), Cristescu, N. D., and Hardy, Jr., H. Reginald, eds., Trans Tech Publ., Clausthal-Zellerfeld. CHAPTER 2.4 Elastomer Models R. W. OGDEN Department of Mathematics, University of Glasgow, Glasgow G12 8QW, UK Contents 2.4.1 Validity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 2.4.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 2.4.3 Description of the Model. . . . . . . . . . . . . . . . . . . . . 93 2.4.4 Identification of Parameters . . . . . . . . . . . . . . . . . . 93 2.4.5 How to Use It . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 2.4.6 Table of Parameters. . . . . . . . . . . . . . . . . . . . . . . . . . 94 References. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 2.4.1 VALIDITY Many rubberlike solids can be treated as isotropic and incompressible elastic materials to a high degree of approximation. Here describe the mechanical properties of such solids through the use of an isotropic elastic strain-energy function in the context of finite deformations. For general background on finite elasticity, we refer to Ogden . 2.4.2 BACKGROUND Locally, the finite deformation of a material can be described in terms of the three principal stretches, denoted by l1 ; l2 ; and l3 . For an incompressible material these satisfy the constraint l1 l2 l3 ¼ 1 ð1Þ The material is isotropic relative to an unstressed undeformed (natural) configuration, and its elastic properties are characterized in terms of a Handbook of Materials Behavior Models Copyright # 2001 by Academic Press. All rights of reproduction in any form reserved. 91 92 Ogden strain-energy function Wðl1 ; l2 ; l3 Þ per unit volume, where W depends symmetrically on the stretches subject to Eq. 1. The principal Cauchy stresses associated with this deformation are given by si ¼ li @W p; @li ð2Þ i 2 f1; 2; 3g where p is an arbitrary hydrostatic pressure (Lagrange multiplier). By regarding two of the stretches as independent and treating the strain energy as 1 # 1 ; l2 Þ ¼ Wðl1 ; l2 ; l1 a function of these through the definition Wðl 1 l2 Þ, we obtain s1 s3 ¼ l1 # @W @l1 s2 s3 ¼ l2 # @W @l2 ð3Þ For consistency with the classical theory, we must have # Wð1; 1Þ ¼ 0; # # # @2W @W @2W ð1; 1Þ ¼ 2m; ð1; 1Þ ¼ 0; ð1; 1Þ ¼ 4m; 2 @la @l1 @l2 @la ð4Þ a 2 f1; 2g where m is the shear modulus in the natural configuration. The equations in Eq. 3 are unaffected by superposition of an arbitrary hydrostatic stress. Thus, # and hence those of W, it suffices to set in determining the characteristics of W, s3 ¼ 0 in Eq. 3, so that s1 ¼ l1 # @W @l1 s2 ¼ l2 # @W @l2 ð5Þ Biaxial experiments in which l1 ; l2 and s1 ; s2 are measured then provide # Biaxial deformation of a thin sheet where the data for the determination of W. deformation corresponds effectively to a state of plane stress, or the combined extension and inflation of a thin-walled (membranelike) tube with closed ends provide suitable tests. In the latter case the governing equations are written 1 P * ¼ l1 1 l2 # @W @l2 F* ¼ # 1 # @W @W l2 l1 1 @l1 2 @l2 ð6Þ where P * ¼ PR=H, P is the inflating pressure, H the undeformed membrane thickness, and R the corresponding radius of the tube, while F * ¼ F=2pRH, with F the axial force on the membrane (note that the pressure contributes to the total load on the ends of the tube). Here l1 is the axial stretch and l2 the azimuthal stretch in the membrane. 93 2.4 Elastomer Models 2.4.3 DESCRIPTION OF THE MODEL A specific model which fits very well the available data on various rubbers is that defined by N X ð7Þ mn ðla1n þ la2n þ la3n 3Þ=an W¼ n¼1 where mn and an are material constants and N is a positive integer, which for many practical purposes may be taken as 2 or 3 . For consistency with Eq. 4 we must have N X ð8Þ mn an ¼ 2m n¼1 and in practice it is usual to take mn an > 0 for each n ¼ 1; . . . ; N. In respect of Eq. 7, the equations in Eq. 3 become N N X X s1 s3 ¼ mn ðla1n la3n Þ s2 s3 ¼ mn ðla2n la3n Þ n¼1 ð9Þ n¼1 2.4.4 IDENTIFICATION OF PARAMETERS Biaxial experiments with s3 ¼ 0 indicate that the shapes of the curves of s1 s2 plotted against l1 are essentially independent of l2 for many rubbers. Thus the shape may be determined by the pure shear test with l2 ¼ 1, so that N N X X s1 s2 ¼ mn ðla1n 1Þ s2 ¼ mn ðla3n 1Þ ð10Þ n¼1 n¼1 for l1 1; l3 1. The shift factor to be added to the first equation in Eq. 10 when l2 differs from 1 is N X ð11Þ mn ð1 la2n Þ n¼1 Information on both the shape and shift obtained from experiments at fixed l2 then suffice to determine the material parameters, as described in detail in References  or . Data from the extension and inflation of a tube can be studied on this basis by considering the combination of equations in Eq. 6 in the form s1 s 2 ¼ l 1 # # @W @W 1 l2 ¼ l1 F * l22 l1 P * @l1 @l2 2 ð12Þ 94 Ogden 2.4.5 HOW TO USE IT The strain-energy function is incorporated in many commercial Finite Element (FE) software packages, such as ABAQUS and MARC, and can be used in terms of principal stretches and principal stresses in the FE solution of boundary-value problems. 2.4.6 TABLE OF PARAMETERS Values of the parameters corresponding to a three-term form of Eq. 7 are now given in respect of two different but representative vulcanized natural rubbers. The first is the material used by Jones and Treloar : a1 ¼ 1:3; a2 ¼ 4:0; a3 ¼ 2:0; m1 ¼ 0:69; m2 ¼ 0:01; m3 ¼ 0:0122 Nmm2 The second is the material used by James et al. , the material constants having been obtained by Treloar and Riding : a1 ¼ 0:707; a2 ¼ 2:9; a3 ¼ 2:62; m1 ¼ 0:941; m2 ¼ 0:093; m3 ¼ 0:0029 Nmm2 For detailed descriptions of the rubbers concerned, reference should be made to these papers. REFERENCES 1. James, A. G., Green, A., and Simpson, G. M. (1975). Strain energy functions of rubber. I. Characterization of gum vulcanizates. J. Appl. Polym. Sci. 19: 2033–2058. 2. Jones, D. F., and Treloar, L. R. G. (1975). The properties of rubber in pure homogeneous strain. J. Phys. D: Appl. Phys. 8: 1285–1304. 3. Ogden, R. W. (1982). Elastic deformations of rubberlike solids, in Mechanics of Solids (Rodney Hill 60th Anniversary Volume) pp. 499–537, Hopkins, H. G., and Sevell, M. J., eds., Pergamon Press. 4. Ogden, R. W. (1986). Recent advances in the phenomenological theory of rubber elasticity. Rubber Chem. Technol. 59: 361–383. 5. Ogden, R. W. (1997). Non-Linear Elastic Deformations, Dover Publications. 6. Treloar, L. R. G., and Riding, G. (1979). A non-Gaussian theory for rubber in biaxial strain. I. Mechanical properties. Proc. R. Soc. Lond. A369: 261–280. CHAPTER 2.5 Background on Viscoelasticity KOZO IKEGAMI Tokyo Denki University, Kanda-Nishikicho 2-2, Chiyodaku, Tokyo 101-8457, Japan Contents 2.5.1 2.5.2 2.5.3 2.5.4 2.5.5 2.5.6 Validity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 Mechanical Models . . . . . . . . . . . . . . . . . . . . . . . . 95 Static Viscoelastic Deformation. . . . . . . . . . . . . . . 98 Dynamic Viscoelastic Deformation . . . . . . . . . 100 Hereditary Integral . . . . . . . . . . . . . . . . . . . . . . . . 102 Viscoelastic Constitutive Equation by the Laplace Transformation . . . . . . . . . . . . . . . . . . . . 103 2.5.7 Correspondence Principle . . . . . . . . . . . . . . . . . . 104 References. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 2.5.1 VALIDITY Fundamental deformation of materials is classified into three types: elastic, plastic, and viscous deformations. Polymetric material shows time-dependent properties even at room temperature. Deformation of metallic materials is also time-dependent at high temperature. The theory of viscoelasticity can be applied to represent elastic and viscous deformations exhibiting timedependent properties. This paper offers an outline of the linear theory of viscoelasticity. 2.5.2 MECHANICAL MODELS Spring and dashpot elements as shown in Figure 2.5.1 are used to represent elastic and viscous deformation, respectively, within the framework of the Handbook of Materials Behavior Models Copyright # 2001 by Academic Press. All rights of reproduction in any form reserved. 95 96 Ikegami FIGURE 2.5.1 Mechanical model of viscoelasticity. linear theory of viscoelasticity. The constitutive equations between stress s and stress e of the spring and dashpot are, respectively, as follows: de ð1Þ s ¼ ke s¼Z dt where the notations k and Z are elastic and viscous constants, respectively. Stress of spring elements is linearly related with strain. Stress of dashpot elements is related with strain differentiated by time t, and the constitutive relation is time-dependent. Linear viscoelastic deformation is represented by the constitutive equations combining spring and dashpot elements. For example, the constitutive equations of series model of spring and dashpot shown in Figure 2.5.2 is as follows: Z ds de ð2Þ sþ ¼Z k dt dt This is called the Maxwell model. The constitutive equation of the parallel model of spring and dashpot elements shown in Figure 2.5.3 is as follows: s ¼ ke þ Z This is called the Voigt or Kelvin model. de dt ð3Þ 97 2.5 Background on Viscoelasticity FIGURE 2.5.2 Maxwell model. There are many variations of constitutive equations giving linear viscoelastic deformation by using different numbers of spring and dashpot elements. Their constitutive equations are generally represented by the following ordinary differential equation: p0 s þ p1 ds d2 s dn s þ p2 2 þ . . . þ pn n dt dt dt ¼ q0 e þ q1 de d2 e dn e þ q2 2 þ . . . þ qn n dt dt dt ð4Þ The coefficients p and q of Eq. 4 give the characteristic properties of linear viscoelastic deformation and take different values according to the number of spring and dashpot elements of the viscoelastic mechanical model. 98 Ikegami FIGURE 2.5.3 Voigt (Kelvin) model. 2.5.3 STATIC VISCOELASTIC DEFORMATION There are two functions representing static viscoelastic deformation; one is creep compliance, and another is the relaxation modulus. Creep compliance is defined by strain variations under constant unit stress. This is obtained by solving Eqs. 2 or 3 for step input of unit stress. For the Maxwell model and the Voigt model, their creep compliances are represented, respectively, by the following expressions. For the Maxwell model, the creep compliance is t 1 1 t þ1 e þ ¼ ð5Þ Z k k t where tM ¼ Z=k, and this is denoted as relaxation time. For the Voigt model, the creep compliance is 1 kt 1 t e ¼ 1 exp ¼ 1 exp ð6Þ k Z k tk where tK ¼ Z=k, and this is denoted as retardation time. Creep deformations of the Maxwell and Voigt models are illustrated in Figures 2.5.4 and 2.5.5, respectively. Creep strain of the Maxwell model 2.5 Background on Viscoelasticity FIGURE 2.5.4 Creep compliance of the Maxwell model. FIGURE 2.5.5 Creep compliance of the Voigt model. 99 100 Ikegami increases linearly with respect to time duration. The Voigt model exhibits saturated creep strain for a long time. The relaxation modulus is defined by stress variations under constant unit strain. This is obtained by solving Eqs. 2 or 3 for step input of unit strain. For the Maxwell and Voigt models, their relaxation moduli are represented by the following expressions, respectively. For the Maxwell model, kt t s ¼ k exp ¼ k exp ð7Þ Z tM For the Voigt model, s¼k ð8Þ Relaxation behaviors of the Maxwell and Voigt models are illustrated in Figures 2.5.6 and 2.5.7, respectively. Applied stress is relaxed by Maxwell model, but stress relaxation dose not appear in Voigt model. 2.5.4 DYNAMIC VISCOELASTIC DEFORMATION The characteristic properties of dynamic viscoelastic deformation are represented by the dynamic response for cyclically changing stress or strain. FIGURE 2.5.6 Relaxation modulus of the Maxwell model. 2.5 Background on Viscoelasticity 101 FIGURE 2.5.7 Relaxation modulus of the Voigt model. The viscoelastic effect causes delayed phase phenomena between input and output responses. Viscoelastic responses for changing stress or strain are defined by complex compliance or modulus, respectively. The dynamic viscoelastic responses are represented by a complex function due to the phase difference between input and output. Complex compliance J of the Maxwell model is obtained by calculating changing strain for cyclically changing stress with unit amplitude. Substituting changing complex stress s ¼ expðiotÞ, where i is an imaginary unit and o is the frequency of changing stress, into Eq. 2, complex compliance J is obtained as follows: J ¼ 1 1 1 1 i ¼ i ¼ J0 iJ00 k oZ k kotM ð9Þ where the real part J0 ¼ 1=k is denoted as storage compliance, and the imaginary part J00 ¼ 1=kotM is denoted as loss compliance. The complex modulus Y of the Maxwell model is similarly obtained by calculating the complex changing strain for the complex changing strain 102 Ikegami e ¼ expðiotÞ as follows: Y ¼ k ðotM Þ2 1 þ ðotM Þ2 þ ik otM 1 þ ðotM Þ2 ¼ Y 0 þ iY 00 ð10Þ where Y 0 ¼ kððotM Þ2 =ð1 þ ðotM Þ2 ÞÞ and Y 00 ¼ kðotM =ð1 þ ðotM Þ2 ÞÞ. The notations Y 0 and Y 00 are denoted as dynamic modulus and dynamic loss, respectively. The phase difference d between input strain and output stress is given by tan d ¼ Y 00 1 ¼ Y 0 otM ð11Þ This is called mechanical loss. Similarly, the complex compliance and the modulus of the Voigt model are able to be obtained. The complex compliance is " # " # 1 1 1 otK J ¼ i ¼ J0 iJ00 ð12Þ k 1 þ ðotK Þ2 k 1 þ ðotK Þ2 " # " # 1 1 1 otK 00 where J ¼ and J ¼ k 1 þ ðotK Þ2 k 1 þ ðotK Þ2 0 The complex modulus is Y ¼ k þ iotK ¼ Y 0 þ iY 00 ð13Þ where Y 0 ¼ k and Y 00 ¼ kotK . 2.5.5 HEREDITARY INTEGRAL The hereditary integral offers a method of calculating strain or stress variation for arbitrary input of stress or strain. The method of calculating strain for stress history is explained by using creep compliance as illustrated in Figure 2.5.8. An arbitrary stress history is divided into incremental constant stress history ds0 Strain variation induced by each incremental stress history is obtained by creep compliance with the constant stress values. In Figure 2.5.8 the strain induced by stress history for t0 5t is represented by the following integral: Z t ds0 ð14Þ Jðt t0 Þ 0 dt0 eðtÞ ¼ s0 JðtÞ þ dt 0 103 2.5 Background on Viscoelasticity FIGURE 2.5.8 Hereditary integral. This equation is transformed by partially integrating as follows: Z t dJðt t0 Þ 0 dt sðt0 Þ eðtÞ ¼ sðtÞJð0Þ þ dðt t0 Þ 0 Similarly, stress variation for arbitrary strain history becomes Z t ds0 Yðt t0 Þ 0 dt0 sðtÞ ¼ e0 YðtÞ þ dt 0 Partial integration of Eq. & gives the following equation: Z t dYðt t0 Þ 0 dt sðtÞ ¼ eðtÞYð0Þ þ sðt0 Þ dðt t0 Þ 0 ð15Þ ð16Þ ð17Þ Integrals in Eqs. 14 to 17 are called hereditary integrals. 2.5.6 VISCOELASTIC CONSTITUTIVE EQUATION BY THE LAPLACE TRANSFORMATION The constitutive equation of viscoelastic deformation is the ordinary differential equation as given by Eq. 4. That is, n X k¼0 pk m dk s X dk e ¼ q k dtk dtk k¼0 ð18Þ 104 Ikegami This equation is written by using differential operators P and Q, Ps ¼ Qe where P ¼ n P k pk k¼0 m P ð19Þ k d d and Q ¼ qk k . dtk dt k¼0 Equation (1?) is represented by the Laplace transformation as follows. n n X X pk sk s% ¼ qk sk e% ð20Þ k¼0 k¼0 where s% and e% are transformed stress and strain, and s is the variable of the Laplace transformation. Equation 20 is written by using the Laplace % as follows: transformed operators of time derivatives P% and Q % Q ð21Þ s% ¼ e% P% n m P P % ¼ pk sk and Q q k sk . where P% ¼ k¼0 k¼0 % P% Comparing Eq. 21 with Hooke’s law in one dimension, the coefficient Q= corresponds to Young’s modulus of linear elastic deformation. This fact implies that linear viscoelastic deformation is transformed into elastic deformation in the Laplace transformed state. 2.5.7 CORRESPONDENCE PRINCIPLE In the previous section, viscoelastic deformation in the one-dimensional state was able to be represented by elastic deformation through the Laplace transformation. This can apply to three-dimensional viscoelastic deformation. The constitutive relations of linear viscoelastic deformation are divided into the relations between hydrostatic pressure and dilatation, and between deviatoric stress and strain. The relation between hydrostatic pressure and dilatation is represented by m n X dk s0ij X dk eii p0k k ¼ q00k k ð22Þ dt dt k¼0 k¼0 P00 sii ¼ Q00 eii ð23Þ n P dk dk and Q00 ¼ q00k k . In Eq. 22 hydrostatic pressure is (1/3) k dt dt k¼0 k¼0 sii and dilatation is eii . where P00 m P p00k 105 2.5 Background on Viscoelasticity The relation between deviatoric stress and strain is represented by m X p0k k¼0 dk s0ij dtk ¼ n X q0k k¼0 P0 s0ij ¼ Q0 e0ij dk e0ij dtk ð24Þ ð25Þ n P dk dk and Q0 ¼ q0k k . In Eq. 24 deviatoric stress and strain k dt dt k0 k¼0 are s0ij and e0ij , respectively. The Laplace transformations of Eqs. 22 and 24 are written, respectively, as follows: where P0 ¼ m P p0k % 00 e%ii P% 00 s% ii ¼ Q ð26Þ % 00 ¼ Q % sðsÞ, and where P% 00 ¼ P% 00 ðsÞ and Q 00 % 0 e% 0ij P% 0 s% 0ij ¼ Q ð27Þ %0 ¼ Q % 0 ðsÞ. where P% 0 ¼ P% 0 ðsÞ and Q The linear elastic constitutive relations between hydrostatic pressure and dilatation and between deviatoric stress and strain are represented as follows: sii ¼ 3Keii ð28Þ s0ij ¼ 2Ge0ii ð29Þ Comparing Eq. 17 with Eq. 19, and Eq. 18 with Eq. 20, the transformed viscoelastic operators correspond to elastic constants as follows: 3K ¼ % 00 Q P% 00 ð30Þ 2G ¼ %0 Q P% 0 ð31Þ where K and G are volumetric coefficient and shear modulus, respectively. For isotropic elastic deformation, volumetric coefficient K and shear modulus G are connected with Young’s modulus E and Poisson’s ratio n as follows: E 2ð1 þ nÞ ð32Þ E 3ð1 2nÞ ð33Þ G¼ K¼ 106 Ikegami Using Eqs. 30–33, Young’s modulus E and Poisson’s ratio are connected with the Laplace transformed coefficient of linear viscoelastic deformation as follows: % 00 % 0Q 3Q E ¼ 0 00 ð34Þ % %0 2P% Q þ P% 00 Q n¼ % 00 P% 00 Q %0 P% 0 Q 0 00 00 % þ P% Q %0 2P% Q ð35Þ Linear viscoelastic deformation corresponds to linear elastic deformation through Eqs. 30–31 and Eqs. 34–35. This is called the correspondence principle between linear viscoelastic deformation and linear elastic deformation. The linear viscoelastic problem is the transformed linear elastic problem in the Laplace transformed state. Therefore, the linear viscoelastic problem is able to be solved as a linear elastic problem in the Laplace transformed state, and then the elastic constants of solved solutions are replaced with the Laplace transformed operator of Eqs. 30–31 and Eqs. 34–35 by using the correspondence principle. The solutions replaced the elastic constants become the solution of the linear viscoelastic problem by inversing the Laplace transformation. REFERENCES 1. 2. 3. 4. 5. 6. Bland, D. R. (1960). Theory of Linear Viscoelasticity, Pergamon Press. Ferry, J. D. (1960). Viscoelastic Properties of Polymers, John Wiley & Sons. Reiner, M. (1960). Deformation, Strain and Flow, H. K. Lewis & Co. Flluege, W. (1967). Viscoelasticity, Blaisdell Publishing Company. Christensen, R. M. (1971). Theory of Viscoelasticity: An Introduction, Academic Press. Drozdov, A. D. (1998). Mechanics of Viscoelastic Solids, John Wiley & Sons. CHAPTER 2.6 A Nonlinear Viscoelastic Model Based on Fluctuating Modes RACHID RAHOUADJ AND CHRISTIAN CUNAT LEMTA, UMR CNRS 7563, ENSEM INPL 2, avenue de la For#et-de-Haye, 54500 Vandoeuvre-l"esNancy, France Contents 2.6.1 Validity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 2.6.2 Background of the DNLR . . . . . . . . . . . . . . . 108 184.108.40.206 Thermodynamics of Irreversible Processes and Constitutive Laws . . . 108 220.127.116.11 Kinetics and Complementary Laws . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 18.104.22.168 Constitutive Equations of the DNLR . . . . . . . . . . . . . . . . . . . . . . . . 112 2.6.3 Description of the Model in the Case of Mechanical Solicitations . . . . . . . . . . . . . . 113 2.6.4 Identification of the Parameters . . . . . . . . . 113 2.6.5 How to Use It . . . . . . . . . . . . . . . . . . . . . . . . . . 115 2.6.6 Table of Parameters. . . . . . . . . . . . . . . . . . . . . 115 References. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 2.6.1 VALIDITY We will formulate a viscoelastic modeling for polymers in the temperature range of glass transition. This physical modeling may be applied using integral or differential forms. Its fundamental basis comes from a generalization of the Gibbs relation, and leads to a formulation of constitutive laws involving control and internal thermodynamic variables. The latter must traduce Handbook of Materials Behavior Models Copyright # 2001 by Academic Press. All rights of reproduction in any form reserved. 107 108 Rahouadj and Cunat different microstructural rearrangements. In practice, both modal analysis and fluctuation theory are well adapted to the study of the irreversible transformations. Such a general formulation also permits us to consider various nonlinearities as functions of material specificities and applied perturbations. To clarify the present modeling, called ‘‘the distribution of nonlinear relaxations’’ (DNLR), we will consider the viscoelastic behavior in the simple case of small applied perturbations near the thermodynamic equilibrium. In addition, we will focus our attention upon the nonlinearities induced by temperature and frequency perturbations. 2.6.2 BACKGROUND OF THE DNLR 22.214.171.124 THERMODYNAMICS OF IRREVERSIBLE PROCESSES AND CONSTITUTIVE LAWS As mentioned, the present irreversible thermodynamics are based on a generalization of the fundamental Gibbs equation to systems evolving outside equilibrium. Note that Coleman and Gurtin , have also applied this postulate in the framework of rational thermodynamics. At first, a set of internal variables (generalized vector denoted z) is introduced to describe the microstructural state. The generalized Gibbs relation combines the two laws of thermodynamics into a single one, i.e., the internal energy potential: e ¼ eðs; e; n; . . . ; zÞ ð1Þ which depends on overall state variables, including the specific entropy, s. Furthermore, with the positivity of the entropy production being always respected, one obtains for open systems: n X dDi s T ¼ Tss ¼ Js : rT Jk : rmk þ A z’ 0 ð2Þ dt k¼1 where the nonequilibrium thermodynamic forces may be separated into two groups: (i) the gradient ones, such as the gradient of temperature gradient rT, and the gradient of generalized chemical potential rmk ; and (ii) The generalized forces A, or affinities as defined by De Donder  for chemical reactions, which characterize the nonequilibrium state of a uniform medium. The vectors Js , Jk , and z’ correspond to the dual, fluxes, or ratetype variables. To simplify the formulation of the constitutive laws, we will now consider the behavior of a uniform representative volume element (RVE without any 2.6 A Nonlinear Viscoelastic Model Based on Fluctuating Modes 109 gradient), thus: Tss ¼ A z’ 0 ð3Þ The equilibrium or relaxed state (denoted by the index r) is currently described by a suitable thermodynamic potential (cr ) obtained via the Legendre transformation of Eq. 1 with respect to the control or state variable (g). In this particular state, the set of internal variables is completely governed by (g): cr ¼ cr ðg; zr ðgÞÞ ¼ cr ðgÞ ð4Þ Our first hypothesis  states that it is always possible to define a thermodynamic potential c only as a function of g and z, even for systems outside equilibrium: ð5Þ c ¼ cðg; zÞ Then, we assume that the constitutive equations may be obtained as functions of the first partial derivatives of this potential with respect to the dual variables, and depend consequently on both control and internal variables; i.e., b ¼ bðg; zÞ and A ¼ Aðg; zÞ. In fact, this description is consistent with the principle of equipresence, as postulated in rational thermodynamics. Therefore, the thermodynamic potential becomes in a differential form: q r X X dck ¼ bm dgm Aj dzj ð6Þ m¼1 j¼1 Thus the time evolution of the global response, b, obeys a nonlinear differential equation involving both the applied perturbation g and the internal variable z (generalized vector): ð7aÞ b’ ¼ au : g’ þ b : z’ A’ ¼ t b : g’ g : z’ ð7bÞ This differential system resumes in a general and condensed form the announced constitutive relationships. The symmetrical matrix au ¼ @2 [email protected]@g is the matrix of Tisza, and the symmetrical matrix g ¼ @2 [email protected]@z traduces the interaction between the dissipation processes . The rectangular matrix b ¼ @2 [email protected]@g expresses the coupling effect between the state variables and the dissipation variables. In other respects, the equilibrium state classically imposes the thermo’ ¼ 0. From Eq. 7b dynamic forces and their rate to be zero; i.e., A ¼ 0 and A we find, for any equilibrium state, that the internal variables’ evolution results directly from the variation of the control variables: z’ r ¼ g1 : t b : g’ ð8Þ 110 Rahouadj and Cunat According to Eqs. 7b and 8, the evolution of the generalized force becomes ’ ¼ g : ð’z z’ r Þ ð9Þ A and its time integration for transformation near equilibrium leads to the simple linear relationship A ¼ gðz zr Þ ð10Þ where g is assumed to be constant. 126.96.36.199 KINETICS AND COMPLEMENTARY LAWS To solve the preceding three equations (7a–b, 10), with the unknown variables being b, z, zr , and A, one has to get further information about the kinetic relations between the nonequilibrium driving forces A and their fluxes z’ . 188.8.131.52.1 First-Order Nonlinear Kinetics and Relaxation Times We know that the kinetic relations are not submitted to the same thermodynamic constraints as the constitutive ones. Thus we shall consider for simplicity an affine relation between fluxes and forces. Note that this wellknown modeling, early established by Onsager, Casimir, Meixner, de Donder, De Groot, and Mazur, is only valid in the vicinity of equilibrium: z’ ¼ L A ð11Þ z’ ¼ L g ðz zr Þ ¼ t1 :ðz zr Þ ð12Þ and hence, with Eq. 10: According to this nonlinear kinetics, Meixner  has judiciously suggested a base change in which the relaxation time operator t is diagonal. Here, we consider this base, which also represents a normal base for the dissipation modes. In what follows, the relaxation spectrum will be explicitly defined on this normal base. To extend this kinetic modeling to nonequilibrium transformations, which is the object of the nonlinear TIP, we also suggest referring to Eq. 12 but with variable relaxation times. Indeed, each relaxation time is inversely proportional to the jump frequency, u, and to the probability þ;r pj ¼ expðDFþ;r j =RTÞ of overcoming a free energy barrier, DFj . It follows that the relaxation time of the process j may be written: trj ¼ 1=u expðDFþ;r j =RTÞ ð13Þ where the symbol (þ) denotes the activated state, and the index (r) refers to the activation barrier of the REV near the equilibrium. 2.6 A Nonlinear Viscoelastic Model Based on Fluctuating Modes 111 The reference jump frequency, u0 ¼ kB T=h, has been estimated from Guggenheim’s theory, which considers elementary movements of translation at the atomic level. The parameters h, kB , and r represent the constants of Plank, Boltzmann, and of the perfect gas, respectively, and T is the absolute temperature. It seems natural to assume that the frequency of the microscopic rearrangements is mainly governed by the applied perturbation rate, g’ , through a shift function að’gÞ: u ¼ u0 =að’gÞ ð14Þ Assuming now that the variation of the activation energy for each process is governed by the evolution of the overall set of internal variables leads us to the following approximation of first order: þ;r þ Kz :ðz zr Þ DFþ j ¼ DFj ð15Þ In the particular case of a viscoelastic behavior, this variation of the free energy becomes negligible. The temperature dependence obviously intervenes into the basic definition of the activation energy as ¼ DEþ;r T DSþ;r DFþ;r j j ð16Þ where the internal energy DEþ;r is supposed to be the same for all processes. It follows that we may define another important shift function, noted aðTÞ, which accounts for the effect of temperature. According to the Arrhenius approximation, DEþ;r being quasi-constant, this shift function verifies the following relation: ln aðT; Tref Þ ¼ DEþ;r ð1=T 1=Tref Þ ð17Þ where Tref is a reference temperature. For many polymers near the glass transition, this last shift function obeys the WLF empiric law developed by William, Landel, and Ferry : lnðaT Þ ¼ c1 ðT Tref Þ=½c2 þ ðT Tref Þ ð18Þ In summary, the relaxation times can be generally expressed as tj ðTÞ ¼ trj ðTref ÞaðT; Tref Þ að’gÞ aðz; zr Þ ð19Þ r and the shift function aðz; z Þ becomes negligible in viscoelasticity. 184.108.40.206.2 Form of the Relaxation Spectrum near the Equilibrium We now examine the distribution of the relaxation modes evolving during the solicitation. In fact, this applied solicitation, g, induces a state of fluctuations which may be approximately compared to the corresponding equilibrium one. According to prigogine , these fluctuations obey the equipartition of the entropy production. Therefore, we can deduce the expected distribution in 112 Rahouadj and Cunat the vicinity of equilibrium as n qﬃﬃﬃﬃ X p0j ¼ B trj with p0j ¼ 1 and B¼1 X n qﬃﬃﬃﬃ trj j¼1 ð20Þ j¼1 where trj is the relaxation time of the process j, p0j its relative weight in the overall spectrum, and n the number of dissipation processes . As a first approximation, the continuous spectrum defined by Eq. 20 may be described with only two parameters: the longest relaxation time corresponding to the fundamental mode, and the spectrum width. Note that a regular numerical discretization of the relaxation time scale using a sufficiently high number n of dissipation modes, e.g., 30, gives a sufficient accuracy. 220.127.116.11 CONSTITUTIVE EQUATIONS OF THE DNLR Combining Eqs. 7a and 12 gives, whatever the chosen kinetics, b’ ¼ au : g’ b ðz zr Þ: t1 ¼ au : g’ ð#a a# r Þ:t1 b z ð21aÞ To simplify the notation, tb will be denoted t. In a similar form and after introducing each process contribution in the base defined above, one has n n b p0 br X X jm j m b’ m ¼ aump g’ p ð21bÞ t j p¼1 j¼1 where the indices u and r denote the instantaneous and the relaxed values, respectively. Now we shall examine the dynamic response due to sinusoidally varying perturbations gn ¼ g0 expðiotÞ, where o is the applied frequency, and i2 ¼ 1, i.e., g’ n ¼ iogn . The response is obtained by integrating the above differential relationship. Evidently, the main problem encountered in the numerical integration consists in using a time step that must be consistent with the applied frequency and the shortest time of relaxation. Furthermore, a convenient possibility for very small perturbations is to assume that the corresponding response is periodic and out of phase: ð22Þ b ¼ b0 expðiot þ jÞ and b’ ¼ iob n n n where j is the phase angle. In fact, such relations are representative of various physical properties as shown by Kramers  and Kronig . The coefficients of the matrices of Tisza, au and ar , and the relaxation times, tj , may be dependent on temperature and=or frequency. In uniaxial 2.6 A Nonlinear Viscoelastic Model Based on Fluctuating Modes 113 tests of mechanical damping, these Tisza’s coefficients correspond to the storage and loss modulus E0 (or G0 ) and E00 (or G00 ), respectively. 2.6.3 DESCRIPTION OF THE MODEL IN THE CASE OF MECHANICAL SOLICITATIONS We consider a mechanical solicitation under an imposed strain e. Here, the perturbation g and the response b are respectively denoted e and s. According ’ may be finally written to Eqs. 19 and 21b, the stress rate response, s, n n X X sj p0j ar : e s’ ¼ p0j au : e’ ð23Þ að’eÞ aðe; er Þ aðT; Tref Þtj ðTref Þ j¼1 j¼1 As an example, for a pure shear stress this becomes n n X X sj 12 p0j Gr e12 s’ 12 ¼ p0j Gu e’ 12 að’eÞ aðe; er Þ aðT; Tref ÞtGj ðTref Þ j¼1 j¼1 ð24Þ In the case of sinusoidally varying deformation, the complex modulus is given by G ðoÞ ¼ Gu þ ðGr Gu Þ n X j¼1 p0j 1 1 þ iotGj It follows that its real and imaginary components are, respectively, n X 1 p0j G0 ðoÞ ¼ Gu þ ðGr Gu Þ 1 þ o2 ðtGj Þ2 j¼1 G00 ðoÞ ¼ ðGr Gu Þ n X j¼1 p0j otj 1 þ o2 ðtGj Þ2 ð25Þ ð26Þ ð27Þ 2.6.4 IDENTIFICATION OF THE PARAMETERS The crucial problem in vibration experiments concerns the accurate determination of the viscoelastic parameters over a broad range of frequency. Generally, to avoid this difficulty one has recourse to the appropriate principle of equivalence between temperature and frequency, assuming implicitly identical microstructural states. A detailed analysis of the literature has brought us to a narrow comparison of the empirical model of Havriliak and 114 Rahouadj and Cunat Negami (HN)  with the DNLR. The HN approach appears to be successful for a wide variety of polymers; it combines the advantages of the previous modeling of Cole and Cole  and of Davidson and Cole . For pure shear stress the response given by this HN approach is G ¼ GuHN þ ðGrHN GuHN Þ 1 ½1 þ ðiotHN Þa b ð28Þ where GuHN ; GrHN ; a; and b are empirical parameters. Thus the real and imaginary components are, respectively, G0 ¼ GuHN þ ðGrHN GuHN Þ G00 ¼ ðGrHN GuHN Þ cosðbyÞ ½1 þ 2oa taHN cosðap=2Þ þ o2a t2a b=2 sinðbyÞ ½1 þ 2oa taHN cosðap=2Þ þ o2a t2a b=2 ð29Þ ð30Þ The function y is defined by y ¼ tan1 oa taHN sinðap=2Þ 1 þ oa taHN cosðap=2Þ ð31Þ Eqs. 28 to 30 are respectively compared to Eqs. 25 to 27 in order to establish a correspondence between the relaxation times of the two models: logðtGr j Þ ¼ logðtHN Þ þ jL=n þ Y ð32Þ where Y, L, and n are a scale parameter, the number of decades of the spectrum, and the number of processes, respectively. A precise empirical connection is obtained by identifying the shift function for the time scale with the relation tanðbyÞ Gr G Gr Gr tj ¼ að’gÞtj ¼ aðoÞtj ¼ ð33Þ tj otHN This involves a progressive evolution of the difference of modulus as a function of the applied frequency: ðGr Gu Þ ¼ ðGrHN GuHN ÞfG ð34Þ The function fG is given by fG ¼ cosðbyÞ ð1 þ tan2 ðbyÞÞ b=2 ½1 þ 2oa taHN cosðap=2Þ þ o2a t2a HN ð35Þ 2.6 A Nonlinear Viscoelastic Model Based on Fluctuating Modes 115 2.6.5 HOW TO USE IT In practice, knowledge of the only empirical parameters of HN’s modeling (and=or Cole and Cole’s and Davidson and Cole’s) permits us, in the framework of the DNLR, to account for a large variety of loading histories. 2.6.6 TABLE OF PARAMETERS As a typical example given by Hartmann et al. , we consider the case of a polymer whose chemical composition is 1PTMG2000=3MIDI=2DMPD*. The master curve is plotted at 298 K in Figure 2.6.1. The spectrum is discretized FIGURE 2.6.1 Theoretical simulation of the moduli for PTMG ( J).* FIGURE 2.6.2 Theoretical simulations of the shift function aðoÞ and of fG for PTMG.* * PTMG: poly (tetramethylene ether) glycol; MIDI: 4,40 -diphenylmethane diisocyanate; DMPD: 2,2-dimethyl-1, 3-propanediol with a density of 1.074 g=cm3, and glass transition Tg ¼ 408C. 116 Rahouadj and Cunat with L ¼ 6, a scale parameter Y equal to 5.6, and 50 relaxation times. The parameters GrHN ¼ 2:14 MPa, GuHN ¼ Gu ¼ 1859 MPa, tHN ¼ 7 1.649 10 s, a ¼ 0:5709 and b ¼ 0:0363 allow us to calculate the shift function aðoÞ and the function fG which is necessary to estimate the difference between the relaxed and nonrelaxed modulus, taking into account the experimental conditions. Figure 2.6.1 illustrates the calculated viscoelastic response, which is superposed to HN’s one. The function fG and the shift function aðoÞ illustrate the nonlinearities introduced in the DNLR modeling (Fig. 2.6.2). REFERENCES 1. Coleman, B. D., and Gurtin, M. (1967). J. Chem. Phys. 47 (2): 597. 2. De Donder, T. (1920). Lecon de thermodynamique et de chimie physique, Paris: Gauthiers, Villars. 3. Cunat, C. (1996). Rev. Gçn. Therm. 35: 680–685. 4. Meixner, J. Z. (1949). Naturforsch., Vol. 4a, p. 504. 5. William, M. L., Landel, R. F., and Ferry, J. D. (1955). The temperature dependence of relaxation mechanisms in amorphous polymers and other glass-forming liquids. J. Amer. Chem. Soc. 77: 3701. 6. Prigogine, I. (1968). Introduction a" la thermodynamique des processus irr!eversibles, Paris: Dunod. 7. Kramers, H. A. (1927). Atti. Congr. dei Fisici, Como, 545. 8. Kronig, R. (1926). J. Opt. Soc. Amer. 12: 547. 9. Havriliak, S., and Negami, S. (1966). J. Polym. Sci., Part C, No. 14, ed. R. F. Boyer, 99. 10. Cole, K. S., and Cole, R. H. (1941). J. Chem. Phys. 9: 341. 11. Davidson, D. W., and Cole, R. H. (1950). J. Chem. Phys. 18: 1417. 12. Hartmann, B., Lee, G. F., and Lee, J. D. (1994). J. Acoust. Soc. Amer. 95 (1). CHAPTER 2.7 Linear Viscoelasticity with Damage R. A. SCHAPERY Department of Aerospace Engineering and Engineering Mechanics, University of Texas, Austin, Texas Contents 2.7.1 2.7.2 2.7.3 2.7.4 Validity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Description of the Model. . . . . . . . . . . . . . . . . . . Identification of the Material Functions and Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7.5 How to Use It . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 118 119 121 123 123 2.7.1 VALIDITY This paper describes a homogenized constitutive model for viscoelastic materials with constant or growing distributed damage. Included are threedimensional constitutive equations and equations of evolution for damage parameters (internal state variables, ISVs) which are measures of damage. Anisotropy may exist without damage or may develop as a result of damage. For time-independent damage, the specific model covered here is that for a linearly viscoelastic, thermorheologically simple material in which all hereditary effects are expressed through a convolution integral with one creep or relaxation function of reduced time; nonlinear effects of transient crack face contact and friction are excluded. More general cases that account for intrinsic nonlinear viscoelastic and viscoplastic effects as well as thermorheologically complex behavior and multiple relaxation functions are published elsewhere . Handbook of Materials Behavior Models Copyright # 2001 by Academic Press. All rights of reproduction in any form reserved. 117 118 Schapery 2.7.2 BACKGROUND As background to the model with time-dependent damage, consider first the constitutive equation with constant damage, in which e and s represent the strain and stress tensors, respectively, e ¼ fSdsg þ eT ð1Þ where S is a fully symmetric, fourth order creep compliance tensor and eT is the strain tensor due to temperature and moisture (and other absorbed substances which affect the strains). The braces are abbreviated notation for a linear hereditary integral. Although the most general form could be used, allowing for general aging effects, for notational simplicity we shall use the familiar form for thermorheologically simple materials, Z t Z x @g 0 @g 0 ð2Þ f fdgg ¼ f ðx x Þ 0 dt ¼ f ðx x0 Þ 0 dx0 @t @x o o where it is assumed f ¼ g ¼ o for t5o and Z t x dt00 =aT ½Tðt00 Þ x0 ¼ xðt0 Þ ð3Þ o Also, aT ðTÞ is the temperature-dependent shift factor. If the temperature is constant in time, then x x0 ¼ ðt t0 Þ=aT : Physical aging  may be taken into account by introducing explicit time dependence in aT ; i.e., use aT ¼ aT ðT; t00 Þ in Eq. 3. The effect of plasticizers, such as moisture, may also be included in aT : When Eq. 2 is used with Eq. 1, f and g are components of the creep compliance and stress tensors, respectively. In certain important cases, the creep compliance components are proportional to one function of time, S ¼ kD ð4Þ where k is a constant, dimensionless tensor and D ¼ DðxÞ is a creep compliance (taken here to be that obtained under a uniaxial stress state). Isotropic materials with a constant Poisson’s ratio satisfy Eq. 4. If such a material has mechanically rigid reinforcements and=or holes (of any shape), it is easily shown by dimensional analysis that its homogenized constitutive equation satisfies Eq. 4; in this case the stress and strain tensors in Eq. 1 should be interpreted as volume-averaged quantities . The Poisson’s ratio for polymers at temperatures which are not close to their glass-transition temperature, Tg , is nearly constant; except at time or rate extremes, somewhat above Tg Poisson’s ratio is essentially one half, while below Tg it is commonly in the range 0.35–0.40 . 119 2.7 Linear Viscoelasticity with Damage Equations 1 and 4 give e ¼ fDdðksÞg þ eT ð5Þ s ¼ kI fEdeg kI fEdeT g ð6Þ The inverse is where kI ¼ k for t > o, 1 and E ¼ EðxÞ is the uniaxial relaxation modulus in which, fDdEg ¼ fEdDg ¼ 1 ð7Þ In relating solutions of elastic and viscoelastic boundary value problems, and for later use with growing damage, it is helpful to introduce the dimensionless quantities 1 1 1 eR fEdeg eRT fEdeT g uR fEdug ð8Þ ER ER ER where ER is an arbitrary constant with dimensions of modulus, called the reference modulus; also, eR and eRT are so-called pseudo-strains and uR is the pseudo-displacement. Equation 6 becomes s ¼ CeR CeRT ð9Þ where C ER kI is like an elastic modulus tensor; its elements are called pseudo-moduli. Equation 9 reduces to that for an elastic material by taking E ¼ ER ; it reduces to the constitutive equation for a viscous material if E is proportional to a Dirac delta function of x. The inverse of Eq. 9 gives the pseudo-strain eR in terms of stress, # þ eR ð10Þ eR ¼ Ss T where S# ¼ C1 ¼ k=ER : The physical strain is given in Eq. 5. 2.7.3 DESCRIPTION OF THE MODEL The correspondence principle (CPII in Schapery [4, 8]) that relates elastic and viscoelastic solutions shows that Eqs. 1–10 remain valid, under assumption Eq. 4, with damage growth when the damage consists of cracks whose faces are either unloaded or have loading that is proportional to the external loads. With growing damage k; C, and S# are time-dependent because they are functions of one or more damage-related ISVs; the strain eT may also depend on damage. The fourth-order tensor k must remain inside the convolution integral in Eq. 5, just as shown. This position is required by the correspondence principle. The elastic-like Eqs. 9 and 10 come from Eq. 5, and thus have the appropriate form with growing damage. However, with 120 Schapery healing of cracks, pseudo-stresses replace pseudo-strains because k must appear outside the convolution integral in Eq. 5 . The damage evolution equations are based on viscoelastic crack growth equations or, in a more general context, on nonequilibrium thermodynamic equations. Specifically, let W R and WCR denote pseudo-strain energy density and pseudo-complementary strain energy density, respectively, 1 ð11Þ W R ¼ CðeR eRT ÞðeR eRT Þ F 2 1# þ eRT s þ F WCR ¼ Sss 2 ð12Þ WCR ¼ W R þ seR ð13Þ so that and @W R @WCR R ð14Þ e ¼ @eR @s The function F is a function of damage and physical variables that cause residual stresses such as temperature and moisture. For later use in Section 2.7.4, assume the damage is fully defined by a set of scalar ISVs, Sp (p ¼1, 2, . . . P) instead of tensor ISVs. Thermodynamic forces, which are like energy release rates, are introduced, s¼ fp @W R @Sp ð15Þ or fp @WCR @Sp ð16Þ where the equality of these derivatives follows directly from the total differential of Eq. 13. Although more general forms could be used, the evolution equations for S’p dSp =dx are assumed in the form ð17Þ S’p ¼ S’p ðSq ; fp Þ in which S’p may depend on one or more Sq (q ¼ 1, . . . P), but on only one force fp . The entropy production rate due to damage is non-negative if X fp S’p O ð18Þ p thus satisfying the Second Law of Thermodynamics. It is assumed that when j fp j is less than some threshold value, then S’p ¼ O. 2.7 Linear Viscoelasticity with Damage 121 Observe that even when the stress vanishes, there may be damage growth due to F. According to Eqs. 12 and 16, fp ¼ @WCR 1 @ S# @eR @F ¼ ss þ T s þ 2 @Sp @Sp @Sp @Sp ð19Þ which does not vanish when r ¼ o, unless @[email protected] ¼ 0. The use of tensor ISVs is discussed and compared with scalar ISVs by Schapery . The equations in this section are equally valid for tensor and scalar ISVs. 2.7.4 IDENTIFICATION OF THE MATERIAL FUNCTIONS AND PARAMETERS The model outlined above is based on thermorheologically simple behavior in that reduced time is used throughout, including damage evolution (Eq. 17). In studies of particle-reinforced rubber , this simplicity was found, implying that even the microcrack growth rate behavior was affected by temperature only through viscoelastic behavior of the rubber. If the damage growth is affected differently by temperature (and plasticizers), then explicit dependence may be introduced in the rate (Eq. 17). In the discussion that follows, complete thermorheological simplicity is assumed. The behavior of particle-reinforced rubber and asphalt concrete has been characterized using a power law when fp > o, ð20Þ S’p ¼ ð fp Þap where ap is a positive constant. (For the rubber composite two ISVs, with a1 ¼ 4:5 and a2 ¼ 6, were used for uniaxial and multiaxial behavior, whereas for asphalt one ISV, with a ¼ 2:5, was used for uniaxial behavior.) A coefficient depending on Sp may be included in Eq. 20; but it does not really generalize the equation because a simple change of the variable Sp may be used to eliminate the coefficient. Only an outline of the identification process is given here, but details are provided by Park et al.  for uniaxial behavior and by Park and Schapery  and Ha and Schapery  for multiaxial behavior. Schapery and Sicking  and Schapery  discuss the model’s use for fiber composites. The effects of eT and F are neglected here. (a) The first step is to obtain the linear viscoelastic relaxation modulus EðxÞ and shift factor aT for the undamaged state. This may be done 122 Schapery using any standard method, such as uniaxial constant strain rate tests at a series of rates and temperatures. Alternatively, for example, uniaxial creep tests may be used to find DðxÞ, after which EðxÞ is derived from Eq. 7. (b) Constant strain rate (or stress rate) tests to failure at a series of rates or temperatures may be conveniently used to obtain the additional data needed for identification of the model. (However, depending on the complexity of the material and intended use of the model, unloading and reloading tests may be needed .) Constant strain rate tests often are preferred over constant stress rate tests because meaningful post-stress peak behavior (prior to significant strain localization) may be found from the former tests. # where For isothermal, constant strain rate, R, tests, the input is Rt ¼ Rx; R# ¼ RaT and x ¼ t=aT . Inasmuch as the model does not depend on temperature when reduced time is used, all stress vs. reduced time response # regardless of temperature. curves depend on only one input parameter R, Thus, one may obtain a complete identification of the model from a series of tests over a range of R# using one temperature and different rates or one rate and different temperatures; both types of tests may be needed in practice for R# to cover a sufficiently broad range. One should, however, conduct at least a small number of both types of tests to check the thermorheologically simple assumption. (c) Convert all experimental values of displacements and strains from step (b) tests to pseudo-quantities using Eq. 8. This removes intrinsic viscoelastic effects, thus enabling all subsequent identification steps to be those for a linear elastic material with rate-dependent damage. If controlled strain (stress) tests are used, then one would employ W R ðWCR Þ in the identification. However, mixed variables may be input test parameters, such as constant strain rate tests of specimens in a test chamber at a series of specified pressures . In this case it is convenient to use mixed pseudoenergy functions in terms of strain and stress variables. Appropriate energy functions may be easily constructed using methods based on linear elasticity theory. (d) The procedure for finding the exponent a and pseudo Young’s modulus in terms of one damage parameter is given by Park et al. . After this, the remaining pseudo-moduli or compliances may be found in terms of one or more ISVs, as described by Park and Schapery  using constant strain rate tests of bar specimens under several confining pressures. The material employed by them was initially isotropic, but it became transversely isotropic as a result of damage. Identification of the full set of five pseudo-moduli and the pseudo-strain energy function, as functions of two ISVs, is detailed by Ha and Schapery . 2.7 Linear Viscoelasticity with Damage 123 2.7.5 HOW TO USE IT Implementation of user-defined constitutive relations based on this model in a finite element analysis is described by Ha and Schapery . Included are comparisons between theory and experiment for overall load-displacement behavior and for local strain distributions. The model employed assumes the material is locally transversely isotropic with the axis of isotropy assumed parallel to the local maximum principal stress direction, accounting for prior stress history at each point. A procedure is proposed by Schapery  that enables use of the same model when transverse isotropy is lost due to rotation of the local maximum principal stress direction. REFERENCES 1. Ha, K., and Schapery, R. A. (1998). A three-dimensional viscoelastic constitutive model for particulate composites with growing damage and its experimental validation. International Journal of Solids and Structures 35: 3497–3517. 2. Hashin, Z. (1983). Analysis of composite materials } a survey. Journal of Applied Mechanics 105: 481–505. 3. Park, S. W., Kim, Y. R., and Schapery, R. A. (1996). A viscoelastic continuum damage model and its application to uniaxial behavior of asphalt concrete. Mechanics of Materials 24: 241–255. 4. Park, S. W., and Schapery, R. A. (1997). A viscoelastic constitutive model for particulate composites with growing damage. International Journal of Solids and Structures 34: 931–947. 5. Schapery, R. A. (1974). Viscoelastic behavior and analysis of composite materials, in Mechanics of Composite Materials, pp. 85–168, vol. 2, Sendeckyi, G. P., ed., New York: Academic. 6. Schapery, R. A. (1981). On viscoelastic deformation and failure behavior of composite materials with distributed flaws, in 1981 Advances in Aerospace Structures and Materials, pp. 5–20, Wang, S. S., and Renton, W. J., eds., ASME, AD-01. 7. Schapery, R. A. (1982). Models for damage growth and fracture in nonlinear viscoelastic particulate composites, in: Proc. Ninth U.S. National Congress of Applied Mechanics, Book No. H00228, pp. 237–245, Pao, Y. H., ed., New York: ASME. 8. Schapery, R. A. (1984). Correspondence principles and a generalized J integral for large deformation and fracture analysis of viscoelastic media, in: International Journal of Fracture 25: 195–223. 9. Schapery, R. A. (1997). Constitutive equations for special linear viscoelastic composites with growing damage, in Advances in Fracture Research, pp. 3019–3027, Karihaloo, B. L., Mai, Y.W., Ripley, M. I., and Ritchie, R. O., eds., Pergamon. 10. Schapery, R. A. (1999). Nonlinear viscoelastic and viscoplastic constitutive equations with growing damage. International Journal of Fracture 97: 33–66. 11. Schapery, R. A., and Sicking, D. L. (1995). On nonlinear constitutive equations for elastic and viscoelastic composites with growing damage, in Mechanical Behavior of Materials, pp. 45–76, Bakker, A., ed., Delft: Delft University Press. 12. Struik, L. C. E. (1978). Physical Aging in Amorphous Polymers and Other Materials, Amsterdam: Elsevier.