Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 50, 2, pp. 357-375, Warsaw 2012 50th Anniversary of JTAM THEORETICAL AND NUMERICAL STUDIES OF RELAXATION DIFFERENTIAL APPROACH IN VISCOELASTIC MATERIALS USING GENERALIZED VARIABLES Claude Chazal Heterogeneous Materials Research Group, University of Limoges, Civil Engineering Department, Egletons, France; e-mail: chazal@unilim.fr Rostand Moutou Pitti Clermont Université, Université Blaise Pascal, Laboratoire de Mécanique et Ingénieries, Clermont Ferrand, France; e-mail: rostand.moutou.pitti@lermont.fr The phenomenon of incrementalization in the time domain, for linear non-ageing viscoelasticmaterials undergoingmechanical deformation, is investigated. Analytical methods of solution are developed for linear vi- scoelastic behavior in two dimensions utilizing generalized variables and realistic material properties. This is accomplished by the use of time- dependent material property characterization through a Dirchilet series representation, thus the transformation of the viscoelastic continuum problem from the integral to a differential form is achieved. The beha- vior equations are derived from linear differential equations based on the discrete relaxation spectrum. This leads to incremental constitutive formulations using the finite difference integration, thus the difficulty of retaining the strain history in computer solutions is avoided.A complete general formulation of linear viscoelastic strain analysis is developed in terms of increments of generalized stresses and strains. Key words: linear viscoelasticity, differential approach, incremental con- stitutive law, discrete relaxation spectrum, generalized variables 1. Introduction Viscoelastic materials are characterized by possessing infinite memory. Their actualmechanical response is a function of thewhole past history of stress and strain. In most cases, the behavior of any linear viscoelastic material may be represented by a hereditary approach based on the Boltzmann superposition 358 C. Chazal, R. Moutou Pitti principle (Boltzmann, 1878). This implies that stress and strain analysis of viscoelastic phenomena which can be observed in the behavior of many real materials, presents many difficulties for real problems of complex geometry. The analysis of linear viscoelastic materials is usually obtained by an applica- tion of the correspondence principle to the equations of elasticity (Birnecker, 1992;Chazal andDubois, 2001;Christensen, 1971). This approach is restricted to problems for which it is possible to find an explicit solution to the associa- ted equations of elasticity. In order to obtain solutions to more complicated problems, it is necessary to develop numerical rather than analytical techni- ques. These numerical methods avoid the retaining of the whole past history of stress and strain in the memory of a digital computer and permit to deal with complex viscoelastic structures involving complicated boundary condi- tions. The key to such methods is to incrementalize the hereditary integral equations bymeans of analytical techniques. Thus, the difficulty of computer storage requirements is avoided and the complete past history of stress and strain is represented bymeans of some auxiliary tensors. The problem of finding incremental formulations of linear non-ageing vi- scoelastic problems has been investigated thoroughly by a number of authors (Jurkiewiez et al., 1999); Chazal andMoutou Pitti, 2010a,b, 2011a,b; Moutou Pitti et al., 2011). The interest for such a formulation lies in its help both in understanding theoretical aspects, especially in the mathematical treatment of the integral equations (for instance, they might turn out to be useful as a tool for deriving approximate governing equations for the behavior of visco- elastic continua) and in developing efficient numerical integration methods on a clearly stated theoretical basis. Several recent studies have addressed the subject of incremental consti- tutive laws for linear non-ageing viscoelastic materials undergoingmechanical deformation (Chazal andMoutouPitti, 2009a,b, 2010a, 2011a,b;MoutouPitti et al., 2011). Early works byChazal andMoutou Pitti (2010b) considered the theory of linear viscoelasticity to establish incremental constitutive equations using creep or relaxation functions. By assuming a strong form of the creep or relaxation function, and by defining the behavior of the material in differen- tial or integral approach that previouslywas proposed by Zocher et al. (1997), incremental equations in the time domain for a linear non-ageing viscoelastic material can be constructed (Chazal, 2000; Chazal and Moutou Pitti, 2009c, 2010a, 2011a,b; Moutou Pitti et al., 2011). However, these formulations deal with local integration for the Volterra equations. Chazal and Moutou Pitti (2009b) established a general method for fin- ding the incremental formulation for any linear ageing viscoelastic problem on Theoretical and numerical studies of relaxation... 359 the basis of the choice of an integrating operator, while Chazal and Moutou Pitti (2009d), Moutou Pitti et al. (2010a,b,c,d) established a similar method for the specific case of the linear non-ageing problem. These methods were successfully applied to general linear material continua in deriving an exten- ded formulation of the classical principle of Boltzmann in the static case (see Chazal 2000; Chazal and Dubois, 2001; Chazal andMoutou Pitti, 2009c; Mo- utou Pitti et al., 2009). Bozza and Gentili (1995) used the theory of linear viscoelasticity to establish constitutive equations using relaxation functions. They sought solutions to the inversion problem of the constitutive equations. Drozdov andDorfmann (2004) derived constitutive equations for the nonline- ar viscoelastic behavior after performing tensile relaxation tests. Kim andLee (2007) and Theocaris (1964) have proposed an incremental formulation and constitutive equations in the finite element context (see also Jurkiewiez et al., 1999; Zocher et al., 1997; Chazal and Poutou Pitti, 2010b, 2011a). In fracture viscoelastic mechanics, Kim and Lee (2007), Moutou Pitti et al. (2009, 2011), Chazal and Dubois (2001) applied the incremental formulation in order to evaluate creep crack growth process in wood. Krempl (1979) andKujawski et al. (1980) performed an experimental study of creep and relaxation in steel at room temperature. However, the formulation used is based on the spectral decomposition using the generalized Maxwell model. All the above procedures adopted in order to transform the original visco- elastic formulation into a new onewith a differential or integral operator were based on the idea of local integration in the global operator of the problem with various techniques. In this paper, a different approach is adopted taking into account that in the presence of a general viscoelastic constitutive law the behavior equations are first integrated over the thickness of the structure in terms of generalized variables. This new formulation is based on a discrete relaxation spectrum and the finite difference method using generalized integral equations in the time domain. The incremental stress and strain constitutive equations are not restricted to isotropic materials and can be used to resolve complex boundary viscoelastic problemswithout retaining the past history of thematerial in the memory of a digital computer. In the following, a formal statement of the viscoelastic initial/boundary value problem is provided. The one dimensional linear viscoelastic behavior is used to account for three dimensional responses. After that, we present the development of generalized differential equations in terms of one dimensional stress and strain components. This is followed by a discussion of the conver- sion through incrementalization (essentially, a finite difference procedure) of 360 C. Chazal, R. Moutou Pitti the linear viscoelastic constitutive equations into a form suitable for imple- mentation in a finite element formulation. Finally, the incremental viscoelastic constitutive equations of the model are established. 2. Problem statement This section concentrates on the viscoelastic response of time dependent ma- terials at isothermal deformationwith small strains. According toChristensen (1971),Mandel (1978), Salençon (1983) andChazal andMoutou Pitti (2010b, 2011a), the components of the relaxation tensor J(t) can be represented in terms of exponential series Jαβγδ(t)= { J∞αβγδ + M∑ m=1 Jmαβγδe −tλm αβγδ } H(t) (2.1) where λmαβγδ,m=1, . . . ,M, are strictly positive scalars and repeated indices do not imply summation convention. J∞αβγδ and J m αβγδ represent the equili- brium and the differed part of the relaxation tensor respectively, and H(t) is the Heaviside unit step function. According to Boltzmann’s superposition principle in linear non-ageing viscoelasticity (Boltzmann, 1878), the constitu- tive equations between the components σαβ(t) of the stress tensor and the components of the strain tensor eαβ(t) for non-ageing linear viscoelastic ma- terials can be expressed in the time domain by hereditary Volterra’s integral equation σαβ(t)= ∑ γ ∑ δ t∫ −∞ Jαβγδ(t− τ) ∂eγδ(τ) ∂τ dτ (2.2) We introduce stresses and strains in generalized variables according to Love’s first-order shell theory. The strain at any point of the thin structure may be given as eαβ(t)= εαβ(t)+ ζχαβ(t) α,β =1,2 (2.3) where εαβ(t) and χαβ(t) are the middle surface extensional strain and cu- rvature, respectively. If we consider a plane stress state, the non-vanishing resultant of stresses is then defined by Nαβ(t)= +h/2∫ −h/2 σαβ(ζ,t) dζ Mαβ(t)= +h/2∫ −h/2 ζσαβ(ζ,t) dζ (2.4) Theoretical and numerical studies of relaxation... 361 Nαβ(t) and Mαβ(t) are the generalized stresses and h is the thickness of the structure assumed to be constant. Note that the radii of the curvature for the middle surface do not enter into equation (2.4) because of the thin shell assumption. In order to determine the constitutive equation in terms of generalized stresses and strains, we introduce generalized strains, given by equation (2.3), into Volterra integral equation (2.2). One finds σNαβ(t)+σ M αβ(t)= ∑ γ ∑ δ t∫ −∞ Jαβγδ(t− τ) ∂ ∂τ [εγδ(τ)+ ζχγδ(τ)] dτ (2.5) Note that the total stress σαβ(t) is separated into two parts: normal stress σNαβ(t) due to extensional strain and bending stress σ M αβ(t) due to curvature. The constitutive equations in generalized variables can now be obtained from behavior equation (2.5). Using equation (2.4) and integrating equation (2.5) over the thickness, we find Nαβ(t)= +h/2∫ −h/2 σNαβ(ζ,t) dζ = ∑ γ ∑ δ h t∫ −∞ Jαβγδ(t− τ) ∂ ∂τ εγδ(τ) dτ Mαβ(t)= +h/2∫ −h/2 ζσMαβ(ζ,t) dζ = ∑ γ ∑ δ h3 12 t∫ −∞ Jεβγδ(t− τ) ∂ ∂τ χγδ(τ) dτ (2.6) Let us consider the two pseudo fourth order tensors N and M of compo- nents Nαβγδ(t) and Mαβγδ(t) respectively. These tensors are defined by { a1Nαβγδ(t) a2Mαβγδ(t) } = t∫ −∞ Jαβγδ(t− τ) ∂ ∂τ { εγδ(τ) χγδ(τ) } (2.7) a1 and a2 are geometric constants and are defined by: a1 = 1/h, and a2 = 12/h 3. Nαβγδ(t) and Mαβγδ(t) are pseudo mono-dimensional stresses obtained from Volterra’s integral equation as given by equation (2.7). The- se components can be interpreted as the contribution of the strain history {eγδ(τ),τ ¬ t} of the components eγδ(t) of the strain tensor to the stress components σαβ(t). Introducing equation (2.7) into equations (2.6), it can be shown that the pseudo mono-dimensional stress components satisfy the following equations 362 C. Chazal, R. Moutou Pitti { Nαβ(t) Mαβ(t) } = 3∑ γ=1 3∑ δ=1 { Nαβγδ(t) Mαβγδ(t) } = 3∑ γ=1 3∑ δ=1 t∫ −∞ Jαβγδ(t− τ) ∂ ∂τ    1 a1 εγδ(τ) 1 a2 χγδ(τ)    (2.8) Each equation of relation (2.8) represents a one-dimensional non-ageing linear viscoelastic material defined by its relaxation function J(t) given by equation (2.1). 3. Analysis of the proposed model When we apply the mechanical strain defined by the strain history {eαβ(τ),τ ∈ ℜ}, the response of the material is then given by the history of stresses {Nαβγδ(t),Mαβγδ(t), t∈ℜ} defined by the behavior equation (2.7) in which the relaxation function is given by equation (2.1). If the generalized strain {εγδ(t),χγδ(t)} is applied to the material at ti- me t, then the response in stresses can be obtained using the finite relaxation spectrum representation given by equation (2.1). This leads to { a1Nαβγδ(t) a2Mαβγδ(t) } = t∫ −∞ [ J∞αβγδ + M∑ m=1 Jmαβγδe −λm αβγδ (t−τ) ] ∂ ∂τ { εγδ(τ) χγδ(τ) } (3.1) Thus the pseudo mono-dimensional stresses given by the last equation, and written as a function of equilibrium and a differed part of the relaxation spectrum, can be rewritten in the following form { Nαβγδ(t) Mαβγδ(t) } =    N ∞ αβγδ(t)+ ∑M m=1N m αβγδ(t) M ∞ αβγδ(t)+ ∑M m=1M m αβγδ(t)    (3.2) with N ∞ αβγδ(t)= 1 a1 t∫ −∞ J∞αβγδ ∂εγδ(τ) ∂τ dτ = 1 a1 J∞αβγδεγδ(t) N m αβγδ(t)= 1 a1 t∫ −∞ Jmαβγδe −λm αβγδ (t−τ)∂εγδ(τ) ∂τ dτ (3.3) M ∞ αβγδ(t)= 1 a2 t∫ −∞ J∞αβγδ ∂χγδ(τ) ∂τ dτ = 1 a2 J∞αβγδχγδ(t) Theoretical and numerical studies of relaxation... 363 M m αβγδ(t)= 1 a2 t∫ −∞ Jmαβγδe −λm αβγδ (t−τ)∂χγδ(τ) ∂τ dτ It should be noted that N ∞ αβγδ(t) and M ∞ αβγδ(t) represent the equilibrium part of thepseudomono-dimensional stress of thematerialwhile N m αβγδ(t) and M m αβγδ(t) represent the differed part of the same pseudo mono-dimensional stress. As we mentioned in the above section, a differential approach is used in order to establish the differential equations of themechanicalmodel. Thus,we need to express the viscoelastic response of thematerial as a function of stress and strain derivatives. For this reason, let us use equation (2.8), the rate of the total stress is determined by ∂ ∂t { Nαβ(t) Mαβ(t) } = 3∑ γ=1 3∑ δ=1 ∂ ∂t { Nαβγδ(t) Mαβγδ(t) } = ∂ ∂t 3∑ γ=1 3∑ δ=1    N ∞ αβγδ(t)+ M∑ m=1 N m αβγδ(t) M ∞ αβγδ(t)+ M∑ m=1 M m αβγδ(t)    (3.4) The rate of the equilibrium part of the pseudo one-dimensional stress N ∞ αβγδ(t) and M ∞ αβγδ(t) is easy to be evaluated. According to equations (3.3)1 and (3.3)3, and after applying a time derivative operator, one find a1 ∂N ∞ αβγδ(t) ∂t = J∞αβγδ ∂εγδ(t) ∂t a2 ∂M ∞ αβγδ(t) ∂t = J∞αβγδ ∂χγδ(t) ∂t (3.5) In other words, the equilibrium part of the pseudo one-dimensional stress is directly proportional to the total strain at time t. It is completely defined by the history of the applied strain. However, the rate of the differed part of the pseudoone-dimensional stresses N m αβγδ(t) and M m αβγδ(t) ismoredifficult to be determined. Using equations (3.3)2 and (3.3)4, and applying a time derivative operator, we can write a1 ∂ ∂t N m αβγδ(t)= J m αβγδe −λm αβγδ (t−t) ∂ ∂t εγδ(t) − t∫ −∞ Jmαβγδλ m αβγδe −λm αβγδ (t−τ) ∂ ∂τ εγδ(τ) (3.6) 364 C. Chazal, R. Moutou Pitti a2 ∂ ∂t M m αβγδ(t)=J m αβγδe −λm αβγδ (t−t) ∂ ∂t χγδ(t) − t∫ −∞ Jmαβγδλ m αβγδe −λm αβγδ (t−τ) ∂ ∂τ χγδ(τ) These integral equations give the total rate of the differed part of the pseudomono-dimensional stresses. The main purpose of our development in this section is to establish dif- ferential equations between the total rate of the pseudo mono-dimensional stresses and the rate of the total strain. For this reason, wewill transform the last equations in a differential type. Let us introduce behavior equations (3.3)2 and (3.3)4 in integral equations (3.6). This leads to linear differential equations with constant coefficients and can be integrated analytically ∂ ∂t N m αβγδ(t)+λ m αβγδN m αβγδ(t)= 1 a1 Jmαβγδ ∂ ∂t εγδ(t) ∂ ∂t M m αβγδ(t)+λ m αβγδM m αβγδ(t)= 1 a2 Jmαβγδ ∂ ∂t εγδ(t) (3.7) The solution to these linear differential equations gives the rate of the pseudo one-dimensional stresses N m αβγδ(t) and M m αβγδ(t). Finally, the general differential equations governing the non-ageing linear viscoelastic behavior can be obtained from equation (3.4) after summation on γ and δ indices. One finds ∂ ∂t { Nαβ(t) Mαβ(t) } = M∑ m=1 ∂ ∂t { Γmαβ(t) Ψmαβ(t) } + 3∑ γ=1 3∑ δ=1 J∞αβγδ ∂ ∂t    1 a1 εγδ(t) 1 a2 χγδ(t)    (3.8) where Γmαβ(t) and Ψ m αβ(t), α,β ∈ {1,2,3}, m ∈ {1, . . . ,M} are the solutions to the following equations { Γmαβ(t) Ψmαβ(t) } = 3∑ γ=1 3∑ δ=1 { N m αβγδ(t) M m αβγδ(t) } (3.9) Note that Γmαβ(t) and Ψ m αβ(t) can be interpreted as pseudo stresses and repre- sent the influence of the past history of strain on thematerial behavior. They are given by the solution to linear differential equations (3.7). It also should be mentioned that the non-ageing linear viscoelastic behavior is completely defined by differential equations (3.8). We note that this formulation, written in terms of generalized stresses and strains rates, is easily adapted to temporal discretization methods such as the finite difference method. Theoretical and numerical studies of relaxation... 365 4. Conversion to incremental equations Here we will describe the solution process of a step-by-step nature in which loads are applied stepwise at various time intervals. Let us consider the time step ∆tn = tn+1 − tn. The subscript n and n+1 refer to the values at the beginning and end of the time step, respectively. This technique is successfully used by Chazal and Dubois (2001) in the case of viscoelastic structures. We assume that the time derivative during each time increment is constant and is expressed by ∂ζij ∂t = ζij(tn+1)− ζij(tn) ∆tn = ∆(ζij)n ∆tn (4.1) where ζij represents generalized strains or stresses. The following expressions can then bewritten for the rate of pseudo stresses at the beginning of the time step ∂ ∂t { Γmαβ(tn) Ψmαβ(tn) } = 1 ∆tn { Γmαβ(tn+1)−Γ m αβ(tn) Ψmαβ(tn+1)−Ψ m αβ(tn) } = 1 ∆tn { ∆Γmαβ(tn) ∆Ψmαβ(tn) } (4.2) A linear approximation is used for strains, and is expressed by { εγδ(τ) χγδ(τ) } = { εγδ(tn) χγδ(tn) } + τ− tn ∆tn H(τ− tn) { εγδ(tn+1)−εγδ(tn) χγδ(tn+1)−χγδ(tn) } (4.3) This linear approximation leads to very accurate results in finite element discretization as it is shown by Chazal and Dubois (2001). Thus we do not need higher approximations for the strain during a finite increment of the time load. This leads to a constant rate during each time increment: ∂ ∂t { εαβ(tn) χαβ(tn) } = 1 ∆tn { εαβ(tn+1)−εαβ(tn) χαβ(tn+1)−χαβ(tn) } = 1 ∆tn { ∆εαβ(tn) ∆χαβ(tn) } (4.4) By integrating equation (3.8) between tn and tn+1, it can be written in the following form { ∆Nαβ(tn) ∆Mαβ(tn) } = M∑ m=1 { ∆Γmαβ(tn) ∆Ψmαβ(tn) } + 3∑ γ=1 3∑ δ=1 J∞αβγδ    1 a1 ∆εγδ(tn) 1 a2 ∆χγδ(tn)    (4.5) In order to determine the generalized stress increments from this equation, we have to determine the generalized pseudo stress increments ∆Γmαβ(tn) and 366 C. Chazal, R. Moutou Pitti ∆Ψmαβ(tn) during the time step ∆tn. First, let us consider differential equation (3.7). The analytical solution to this differential equation can be expressed as N m αβγδ(tn+1)−N m αβγδ(tn)= ( e −λm αβγδ ∆tn −1 ) N m αβγδ(tn) + 1 a1 Jmαβγδ ∆εγδ(tn) λmαβγδ∆tn ( 1− e −λm αβγδ ∆tn ) M m αβγδ(tn+1)−M m αβγδ(tn)= ( e −λm αβγδ ∆tn −1 ) M m αβγδ(tn) + 1 a2 Jmαβγδ ∆χγδ(tn) λmαβγδ∆tn ( 1− e −λm αβγδ ∆tn ) (4.6) Consequently, when we substitute equations (4.6) into equation (3.9), we obtain the generalized pseudo stress increments ∆Γmαβ(tn) and ∆Ψ m αβ(tn)    M∑ m=1 ∆Γmαβ(tn) M∑ m=1 ∆Ψmαβ(tn)    = M∑ m=1 3∑ γ=1 3∑ δ=1 { N m αβγδ(tn+1)−N m αβγδ(tn) M m αβγδ(tn+1)−M m αβγδ(tn) } (4.7) or M∑ m=1 ∆Γmαβ(tn)= M∑ m=1 3∑ α=1 3∑ β=1 ( e −λm αβγδ ∆tn −1 ) N m αβγδ(tn) + 1 a1 Jmαβγδ∆εγδ(tn) λmαβγδ∆tn ( 1− e −λm αβγδ ∆tn ) M∑ m=1 ∆Ψmαβ(tn)= M∑ m=1 3∑ α=1 3∑ β=1 ( e −λm αβγδ ∆tn −1 ) M m αβγδ(tn) + 1 a2 Jmαβγδ∆εγδ(tn) λmαβγδ∆tn ( 1− e −λm αβγδ ∆tn ) (4.8) The incremental constitutive equations can now be obtained from consti- tutive equation (4.5). Substituting equations (4.8) into (4.5), we find { ∆Nαβ(tn) ∆Mαβ(tn) } = 3∑ γ=1 3∑ δ=1 [ Παβγδ(tn) 0 0 Ξαβγδ(tn) ]{ ∆εγδ(tn) ∆χγδ(tn) } − { Ñαβ(tn) M̃αβ(tn) } (4.9) Theoretical and numerical studies of relaxation... 367 where Παβγδ(tn) and Ξαβγδ(tn) are fourth-order tensors which can be inter- preted as rigidity tensors in the extensional and bending state respectively, they are given by Παβγδ(tn)= J ∞ αβγδ + 1 a1 M∑ m=1 Jmαβγδ∆εγδ(tn) λmαβγδ∆tn ( 1− e −λm αβγδ ∆tn ) Ξαβγδ(tn)= J ∞ αβγδ + 1 a2 M∑ m=1 Jmαβγδ∆χγδ(tn) λmαβγδ∆tn ( 1− e −λm αβγδ ∆tn ) (4.10) Ñαβ(tn) and M̃αβ(tn) are pseudo generalized stresses which represent the in- fluence of the complete past history of extensional and bending generalized stresses. They are given by { Ñαβ(tn) M̃αβ(tn) } = 3∑ γ=1 3∑ δ=1 M∑ m=1 ( 1− e −λm αβγδ ∆tn ){Nmαβγδ(tn) M m αβγδ(tn) } (4.11) Finally, the incremental constitutive law given by equation (4.9) can now be inverted to obtain { ∆εγδ(tn) ∆χγδ(tn) } = 3∑ γ=1 3∑ δ=1 [ Θαβγδ(tn) 0 0 Λαβγδ(tn) ]{ ∆Nγδ(tn) ∆Mγδ(tn) } + { ε̃αβ(tn) χ̃αβ(tn) } (4.12) where Θαβγδ(tn) and Λαβγδ(tn) are compliance fourth-order tensors corre- sponding to extensional and bending state of deformation respectively, they are given by the inverse of the rigidity matrix [ Θαβγδ(tn) 0 0 Λαβγδ(tn) ] = [ Παβγδ(tn) 0 0 Ξαβγδ(tn) ] −1 (4.13) ε̃αβ(tn) and χ̃αβ(tn) are pseudo strains tensors which represent the influence of the complete past history of strain. They are given by { ε̃αβ(tn) χ̃αβ(tn) } = 3∑ γ=1 3∑ δ=1 [ Θαβγδ(tn) 0 0 Λαβγδ(tn) ]{ Ñγδ(tn) M̃γδ(tn) } (4.14) The incremental constitutive law represented by equation (4.9) can be introduced in a finite element discretization in order to obtain solutions to complex viscoelastic problems. Finally, in order to use the incremental viscoelastic formulation presented in this paper, we need to identify the relaxation components of the relaxation 368 C. Chazal, R. Moutou Pitti tensor. The experimental identification of viscoelastic properties is treated in details by Jäger and Lackner (2007) and Müllner and Jäger (2008). The viscoelastic solution is obtained by the application of themethod of functional equations to the elastic solution to the indentation problem and by means of torsional rheometry. 5. Finite element discretization The governing equations of the discretized system, using the finite element model, are derived from the principle of virtual displacements. Let us consider a linear quasi-static non-ageing viscoelastic structure. The principle of virtual displacements implies that the increment in external virtual work is equal to the increment in internal virtual work ∫ eA 〈δt∆εαβ)n,δt∆(χαβ)n〉 { ∆(Nαβ)n ∆(Mαβ)n } deA= ∫ eV ∆(fvi )nδt∆(ui)nd eV (5.1) where ∆(ui)n is the incremental displacement field between tn and tn+1, ∆(fvi )n is the incremental body forces per unit volume, eA and eV are the area and the volume of the element, and δt is the variation symbol. For the sake of simplicity, the surface traction term is omitted in the last equation. Assuming small displacements, strains are derived from shape functions using a standardmanner in the context of the finite elementmethod.Using amatrix notation, the strain increment can be written as { ∆εαβ)n ∆(χαβ)n } = [BL]{∆(U e)n} (5.2) where ∆(Ue)n is the local element displacement increment and [BL] is the strain-displacement transformation matrix. Introducing incremental viscoela- stic constitutive equations (4.9) into equilibrium equations (5.1) and using finite element approximation (5.2), the equilibrium equations for linear visco- elastic behavior can be rewritten as [KT ]n{∆(U e)}n = {∆F ext}n+{∆F vis}n (5.3) where [KT ]n = ∫ eA [BL] ⊤[Ωn][BL] d eA (5.4) and Theoretical and numerical studies of relaxation... 369 {∆Fvis}n = ∫ eV [BL] ⊤    1 a1 (Ñαβ)n 1 a2 (M̃αβ)n    deA (5.5) {∆Fext}n = {F ext}n+1 − {F ext}n is the external load vector increment, {∆Fvis}n is the viscous load vector increment corresponding to the complete past history and [Ωn] is the viscoelastic constitutive matrix. The formulation is introduced in the software Cast3m used by the French EnergyAtomicAgency.The software canbe employed for linear viscoelasticity structures using triangular elements. The global incremental procedure for the relaxation differential approach is described as: 1. At time tn, compute the tangentmoduli Παβγδ(tn) and Ξαβγδ(tn) from equations (4.10) 2. Compute the viscoelastic constitutive matrix [Ωn] [Ωn] = [ Παβγδ(tn) 0 0 Ξαβγδ(tn) ] 3. Compute the pseudo generalized stresses {Ñαβ(tn)} and {M̃αβ(tn)} from equation (4.11) 4. Determine the increment of the viscous load vector {∆Fvis}n from equ- ation (5.5) 5. Update the viscoelastic stiffness matrix [KT ]n from equation (5.4) 6. Assemble and solve viscoelastic equilibrium equations (5.3) in order to obtain the displacement increment vector {∆(U)}n 7. Compute the generalized strain increment {∆εαβ(tn)} and {∆χαβ(tn)} from equation (5.2) 8. Compute thegeneralized stress increment {∆Nαβ(tn)}and {∆Mαβ(tn)} from equation (4.9) 9. Using the results of step 3, compute the generalized pseudo-strains {ε̃αβ(tn+1)} and {χ̃αβ(tn+1)} from equation (4.14) 10. Update the state {U}n+1 = {U}n+{∆(U)}n {Nαβ}n+1 = {Nαβ}n+{∆(Nαβ)}n {Mαβ}n+1 = {Mαβ}n+{∆(Mαβ)}n {εαβ}n+1 = {εαβ}n+{∆(εαβ)}n {χαβ}n+1 = {χαβ}n+{∆(χαβ)}n 11. Go to step 1 370 C. Chazal, R. Moutou Pitti 6. Numerical example This example will be illustrated by a viscoelastic circular cylindrical shell fixed at one end and loaded at the free end. The applied load is a unit radial loading at the free end while the other end is built on. The geometrical and loading details are given in Fig.1. It should be noted that the mesh of the cylinder is graded so that there aremore elements near the loaded point, since in this region the stresses and deflections change most rapidly. The material properties used for the cylinder are given in Table 1. Fig. 1. Circular cylindrical shell using axisymetric shell elements Table 1.Constants used for material properties J0 J1 J2 λ1 λ2 1.45 ·10−5 20 ·10−5 33.33 ·10−5 0.001 0.01 The results of the numerical process are shown in Figs. 2-5. In Figure 2, the numerical results for the radial displacement are plotted versus the axial position measured from the free end, while in Fig.3, the meridional moment is plotted for 20-element idealization. Both the radial displacement and the meridionalmoment are comparedwith the theoretical solution given inTimo- shenko et al. (1959). A very good agreement with the theoretical results can be observed. Theoretical and numerical studies of relaxation... 371 The results of the viscoelastic analysis are given in Figs. 4 and 5. Figure 4 shows how the radial deflection varies versus time, while in Fig.5, we plotted the variation of themeridional strain versus time. It can be shown that strains keep on building up leading to the strain failure. Fig. 2. Variation in the radial deflection of the shell undergoing radial load Fig. 3. Variation in the meridional moment in the shell undergoing radial end load Fig. 4. Free end radial displacements versus time in the shell with radial load applied 372 C. Chazal, R. Moutou Pitti Fig. 5. Free endmeridional strains versus time in the shell loaded radially 7. Conclusions The transformation in differential terms of the integral formulation of the viscoelastic continuumproblemhas been successfully achieved through the in- troduction of a discrete spectrum representation of the relaxation tensor. This leads to a new linear incremental formulation in the time domain for non- ageing viscoelastic materials undergoingmechanical deformation. The formu- lation is based on a differential approach using the discrete spectrum repre- sentation for the relaxation components. The governing equations are then obtained using a discretized form of the Boltzmann superposition principle (Boltzmann, 1878). The analytical solution to the differential equations is ob- tained using finite difference discretization in the time domain. In this way, the incremental constitutive equations of the linear viscoelastic material use a pseudo fourth order rigidity tensor and the influence of the whole past hi- story on the behavior of the material at time t is given by a pseudo second order tensor. The generality allowed by this approach has been established by findingan incremental formulation through simple choices of the tensor relaxa- tion components. This approach appears to open awide horizon(to explore) of new incremental formulations according to particular relaxation components. Remarkable incremental constitutive laws, for which the above technique is applied, are given. Among the numerous applications of the incremental formulations presen- ted in this paper, there is numerical implementation in finite element software, thus the behavior of complex boundaryviscoelastic problems can be obtained. References 1. Birnecker R., 1992, Crack tip parameters for growing cracks in linear visco- elastic materials, Proceedings of the 1st International Conference on Localised Damage, Berlin, 85-98 Theoretical and numerical studies of relaxation... 373 2. Boltzmann L., 1878, Zur Theorie der elastischen Nachwirkung Sitzungsber, Mat Naturwiss. Kl. Kaiser. Akad. Wiss., 70, 275 3. Bozza A., Gentili G., 1995, Inversion and quasi-static problems in line- ar viscoelasticity, International Journal of Theoretical and Applied Mechanics: Meccanica, 30, 321-335 4. Chazal C., 2000, A three dimensional incremental formulation for the linear analysis of viscoelastic media,Proceedings of the 3rd International Conference on Mechanics of Time Dependent Materials S.E.M., Erlangen, Germany 5. Chazal C., Dubois F., 2001, A new incremental formulation in the time do- main for crack initiation in an orthotropic linearly viscoelastic solid,Mechanics of Time Dependent Materials, 5, 3, 229-253DOI: 10.1023/A:1017922223094 6. Chazal C., Moutou Pitti R., 2009a, A new incremental formulation for linear viscoelastic analysis: creep differential approach, Journal of Theoretical and Applied Mechanics, 47, 2, 397-409 7. Chazal C., Moutou Pitti R., 2009b, An incremental constitutive law for ageing viscoelastic materials: a three-dimensional approach, Comptes Rendus de l’Académie des Sciences: C.R. Mécanique, 337, 30-33 DOI: 10.1016/j.crme.2008.12.002 8. Chazal C., Moutou Pitti R., 2009c, An incremental formulation for linear viscoelastic analysis: creep differential approach, SEM Annual Conference and Exposition on Experimental and Applied Mechanics, Albuquerque, ISBN 978- 1-935116-03-5 9. Chazal C., Moutou Pitti R., 2009d, Modèle mathématique incrémental par décomposition spectrale pour lesmatériauxviscoélastiques linéaires,27ème Rencontres Universitaires de Génie Civil, AUGC, Saint-Malo, 03-05 10. Chazal C., Moutou Pitti R., 2010a, Modelling of ageing viscoelastic materials in three dimensional finite element approach, International Jour- nal of Theoretical and Applied Mechanics: Meccanica, 45, 3, 439-441 DOI: 10.1007/s11012-009-9244-9 11. Chazal C., Moutou Pitti R., 2010b, Viscoelastic incremental formulation using creep and relaxation differential approaches,Mechanics of Time Depen- dent Materials, 14, 2, 173-190DOI: 10.1007/s11043-009-9101-1 12. Chazal C., Moutou Pitti R., 2011a, Incremental constitutive formulation for time dependentmaterials: creep integral approach, accepted for publication inMechanics of Time DependentMaterials DOI: 10.1007/s11043-011-9135-z 13. Chazal C., Moutou Pitti R., 2011b, Integral approach for time depen- dentmaterials using finite elementmethod, Journal of Theoretical and Applied Mechanics, 49, 4, 1029-1048 14. Christensen R.M., 1971,Theory of Viscoelasticity: an Introduction, Acade- mic Press, NewYork, ISBN 0-12-174250-4 374 C. Chazal, R. Moutou Pitti 15. DrozdovA.D., DorfmannA., 2004,A constitutivemodel in finite viscoela- sticity of particle-reinforced rubbers, International Journal of Theoretical and Applied Mechanics: Meccanica, 39, 245-270 16. Jäger A., Lackner R., et al., 2007, Identification of viscoelastic proper- ties by means of nanoindentation taking the real tip geometry into account, Meccanica, 42, 293-306 17. Jurkiewiez B., Destrebecq J.F, Vergne A., 1999, Incremental analysis of time-dependent effects in composite Structures, Computers and Structures, 73, 425-435 18. Kim K.S, Lee H.S., 2007, An incremental formulation for the prediction of two-dimensional fatigue crack growthwith curved paths, International Journal for Numerical Methods in Engineering, 72, 697-721 19. Krempl E., 1979, An experimental study of uniaxial room-temperature rate- sensitivity, creep and relaxation of AISI type 304 stainless steel, Journal of Mechanics and Physics of Solids, 27, 363-375 20. Kujawski D., Kallianpur V., Krempl E., 1980,An experimental study of uniaxial creep, cyclic creep and relaxation of AISI type 304 stainless steel at room temperature, Journal of Mechanics and Physics of Solids, 28, 129-148 21. Mandel J., 1978, Dissipativité normale et variables caches, Mechanics Rese- arch Communications, 5, 225-229 22. Moutou Pitti R., Chateauneuf A., Chazal C., 2010a,Couplage fiabilité et comportementdifférédes structures enbétonprécontraint,28èmeRencontres Universitaires de Génie Civil, AUGC, La Bourboule, France 23. Moutou Pitti R., Chateauneuf A., Chazal C., 2010b, Fiabilité des structures en béton précontraint avec prise en compte du comportement vi- scoélastique, Fiabilité desMatériaux et des Structures, 6èmes Journées Natio- nales de Fiabilité Toulouse, France 24. Moutou Pitti R., Chazal C., Bouchäır H., 2010c, Modélisation numérique du comportement différé des matériaux du génie civil, ConserBA- TI2010, France 25. Moutou Pitti R., Chazal C., Chateauneuf A., 2010d, Sur une appli- cation numérique du modèle mathématique de la formulation viscoélastique incrémentale en fluage adaptée auxmatériauxdu génie civil, 28ème Rencontres Universitaires de Génie Civil, AUGC, La Bourboule, France 26. Moutou Pitti R., Chazal C., Labesse F., Lapusta Y., 2011, A gene- ralization of Mv integral to axisymmetric problems for viscoelastic materials, accepted for publication inActa Mechanica, DOI: 10.1007/s00707-011-0460-8 Theoretical and numerical studies of relaxation... 375 27. MoutouPittiR., DuboisF., ChazalC., 2009, Initiation and crackgrowth process in viscoelastic orthotropic materials, SEM Annual Conference and Exposition on Experimental and Applied Mechanics, Albuquerque, ISBN 978- 1-935116-03-5 28. Salençon J., 1983, Viscoélasticité, Presse de l’école nationale des ponts et chaussées, Paris 29. Müllner H.W., Jäger A., et al., 2008, Experimental identification of vi- scoelastic properties of rubber compoundsbymeans of torsional rhemetry,Mec- canica, 43, 327-337 30. Theocaris P.S., 1964, Creep and relaxation contraction ratio of linear visco- elastic materials, Journal of Mechanics and Physics of Solids, 12, 125-138 31. Timoshenko S.P., Woinowsky-Krieger S., 1959, Theory of Plates and Shells, McGraw-Hill, NewYork 32. Zocher M.A., Groves S.E., Hellen D.H., 1997, A three-dimensional fini- te element formulation for thermoviscoelastic orthotropic media, International Journal for Numerical Methods in Engineering, 40, 2267-2288 Teoretyczna i numeryczna analiza lepko-sprężystych właściwości materiałów za pomocą różniczkowej metody relaksacji opartej na zmiennych uogólnionych Streszczenie W pracy zbadano zagadnienie inkrementalizacji czasowej zjawisk zachodzących w liniowych, niestarzejących sięmateriałach lepko-sprężystychpoddanychmechanicz- nej deformacji. Zastosowanometody analityczne do określenia lepko-sprężystego za- chowania sięmateriałuw dwuwymiarowej przestrzeni, używając zmiennych uogólnio- nych i realistycznychparametrówokreślającychwłaściwości próbki. Badania przepro- wadzono poprzez zdefiniowanie tych właściwości w postaci szeregu Dirichleta umoż- liwiającego transformację całkowej reprezentacji badanego kontinuum w formę róż- niczkową. Równania stanu zaczerpnięto z liniowego modelu opisanego równaniami różniczkowymibazującymina zdyskretyzowanymwidmie relaksacji.Pozwoliło to uzy- skać konstytutywne wyrażenia przyrostowe poprzez całkowanie różnic skończonych, co z kolei wyeliminowało konieczność zachowywania historii odkształceń w pamięci komputera. Pełna analiza przebiegu deformacji liniowo lepko-sprężystego materiału zostałaprzeprowadzonawdziedzinie przyrostówuogólnionychnaprężeń i odkształceń. Manuscript received February 25, 2011; accepted for print July 6, 2011