Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 49, 4, pp. 1029-1048, Warsaw 2011 INTEGRAL APPROACH FOR TIME DEPENDENT MATERIALS USING FINITE ELEMENT METHOD Claude Chazal Limoges University, Civil Engineering and Durability Team, Egletons, France e-mail : chazal@unilim.fr Rostand Moutou Pitti Clermont Universitè, Universitè Blaise Pascal, Clermont-Ferrand, France e-mail: rostand.moutou.pitti@polytech.univ-bpclermont.fr In this work, we present the development of a mathematical approach for the solution of linear, non-ageing viscoelastic materials undergoing mechanical deformation.We use an integral approach based on a discre- te spectrum representation for the creep tensor in order to derive the incremental viscoelastic formulation. Integral operators are discretized using finite difference techniques. The incremental viscoelastic consti- tutive model contains an internal state variable which represents the influence of the whole past history of stress and strain, thus the difficul- ty of retaining the stress-strainhistory in numerical solutions is avoided. A complete general formulation of linear viscoelastic stress-strain analy- sis is developed in terms of increments of stresses and strains. Numerical simulations are included in order to validate the incremental constitutive equations. Key words: discrete creep spectrum, integral approach, incremental viscoelastic formulation 1. Introduction The increasinguseof viscoelasticmaterials for structural applications expected to operate for long periods of time requires understanding of their mechanical behavior. Stress and strain analysis of viscoelastic phenomena can be observed in the behavior ofmost civil engineeringmaterials, such as concrete, wood and bituminous concrete for road constructions. The assumption of viscoelasticity is one of main characterisation of their behavior and time dependence is a 1030 C. Chazal, R. Moutou Pitti complex problemwhich is of great importance in the determination of stresses and strains in viscoelastic structures. The main problem in computational mechanics is to know the response of such viscoelastic materials taking into account its complete past history of stress and strain or thehereditary loading. Research focuses onmaterial evaluation, characterisation and development of analytical as well as numerical techniques that are capable of accounting for historical effects of stress and environment. The objective of the current work is to present a three-dimensional finite element formulation that is suitable for the analysis of non-ageing, linear visco- elasticmedia.This formulationhasbeen incorporated into a three-dimensional FE program. This code is a general purpose tool capable of predicting the response of a structure to complex loading. The phenomena such as creep, relaxation, and creep-and-recovery can all be predicted using this program. The code also includes automated mesh generators which enable convenient grid generation for problems involving complex geometry. The criterion used in the selection of papers for discussion herein was that the focus be on incremental differential or integral formulations. Among the several formulations proposed in the literature, Chazal and Moutou Pitti (2010b) and Ghazlan et al. (1995) have developed incremental formulation using differential operators. Central to the method was the assumption that the creep compliance could be separated into an “elastic part” and a “creep part” and that the stress could be considered to remain constant across a time step. Chazal and Moutou Pitti (2009a,b, 2010a) have proposed incremental constitutive equations for ageing viscoelastic materials in the finite element context. They discussed conditions under which the requirement of storing all previous solutions could be avoided. Taylor et al. (1970) and Christensen (1980a,b) developed integral constitutive models based on the Volterra he- reditary integral equations. The current solution in such methods is history dependent. Krishnaswamy et al. (1990, 1992) represent modifications of the basic initial strainmethodwhich employ variable stiffness. Thesemethods are muchmore stable but the stiffnessmatrix changeswith time. In fracture visco- elastic mechanics, Chazal and Dubois (2001), Dubois et al. (1998, 2002) and recently Moutou Pitti et al. (2008) have applied the incremental viscoelastic formulation, initially proposed by Ghazlan et al. (1995), in order to evalu- ate crack growth process in wood. However, the formulation used is based on the spectral decomposition technique using a generalized Kelvin Voigt model. Applications of the finite element method in the field of viscoelastic fracture mechanics have been presented by Krishnaswamy et al. (1990), Moran and Knauss (1992) and Chazal and Moutou Pitti (2009a). Brinson and Knauss Integral approach for time dependent materials... 1031 (1992) have used the finite element method to conduct an investigation of viscoelastic micromechanics. Analytical solutions, proposed in the literature, are important to the cur- rent discussion. Duenwald et al. (2009) developed constitutive equations for ligament using experimental tests. Chazal andMoutou Pitti (2009b), Stouffer andWineman(1972) developednewconstitutive equations for linear ageingvi- scoelastic materials taking into account environmental dependent viscoelastic behavior. The problem of nonlinear viscoelastic problems has been addressed in the work of Merodio (2006) and Filograna et al. (2009). These analyti- cal techniques are obviously limited the method to the solution of relatively simple problems. Most numerical methods cannot deal with complex visco- elastic problems, because thesemethods require the retaining of the complete past history of stress and strain. In a recent work, Chazal and Moutou Pitti (2010b) have already proposed incremental viscoelastic formulation based on the creep and relaxation differential approach.Differential operatorswereused to discretize the Volterra integral equations. These methods involve a stepwi- se integration through time. The key to accomplishing this was the use of a Dirichlet-Prony series (in this case a Kelvin model) to represent the kernel of the Volterra integral equation. These methods made the solution of ‘large’ viscoelastic problems possible. For a much more in-depth review, the reader is referred to Chazal andMoutou Pitti (2010b). An efficient incremental formulation, based on discrete creep spectrum, using integral equations in the time domain, is presented. The proposed in- cremental stress strain constitutive equations are not restricted to isotropic materials and can be used to resolve complex viscoelastic problems. The for- mulation is developed todealwith threedimensional viscoelastic problemsand the solutions of a particular time are found from those at the previous time, this leads to great savings in the amount of computer storage requirements needed to solve real problems involving three dimensional loading. 2. Incrementalization of the viscoelastic equations The material is assumed as non-ageing which allows a time independence of instantaneous and long term mechanical properties. This is in agreement with constant environmental conditions. In this case, the fourth-order creep tensor J(t) is written in terms of discrete creep spectrum according to the results of Mandel (1978) and Christensen (1980a,b) 1032 C. Chazal, R. Moutou Pitti J(t)= [ J (0)+ M∑ m=1 ( 1− e−λmt ) J (m) ] H(t) (2.1) In the above, the creep tensor is fitted with aWiechert model. J(0) and J(m), m = 1, . . . ,M, denote for tensors of the fourth rank, H(t) is the Heaviside unit step function and λm are positive scalars. The term J (0) is an elastic compliance while J(m) is a creep compliance function. J(0) and J(m) should be determined in order to represent any particular creep function of interest using experimental data from creep tests. With the linearity assumption, the constitutive equation can be expressed in the time domain by the hereditaryVolterra integral equation. It defines the relationship between strain and stress components, respectively, εij and σkl. According to Boltzmann’s principle superposition (Boltzmann, 1878) applied to linear non-ageing viscoelastic material, the constitutive lawmay be expres- sed in a tensor notation as εij(t)= 3∑ k=1 3∑ l=1 t∫ −∞ Jijkl(t− τ) ∂σkl(τ) ∂τ dτ (2.2) in which τ is the time variable; and t is the time since loading. In equation (2.2), it is assumed that the body is in a stress-free state for t < 0 and the initial response (elastic) is made implicit. The term Jijkl re- presents the fourth-order tensor of creep moduli relating stress to mechanical strain and t is referred to as the current time. Equation (2.2) shows that the strain at anygiven timedependsupon the entire “stress history” σkl(τ), τ < t. The integral in equation (2.2) is called the hereditary integral. The reader will recognize from the form of the constitutive relationship that we have assumed the material to be non-ageing and non-homogeneous. The constitutive rela- tionship given by (2.2) is not suitable for numerical calculus in the context of the finite element method because this leads to the requirement of solving a set of Volterra integrals in order to extract the finite element solution. A different approach will be presented in this work in order to simplify the so- lution of numerical equations for the simulation of viscoelastic behavior. For this reason, constitutive equation (2.2) is discretized in order to establish in- cremental constitutive equations and will lead to the requirement of solving a simple set of algebraic equations. Thus the difficulty of retaining the whole past history in computer solution is avoided. A similar approach has been ta- ken by Ghazlan et al. (1995) and Zocher et al. (1997). In preparation for the developments to follow, we assume that the period of loading is subdivided into discrete intervals ∆tn (time step) such that tn+1 = tn +∆tn. Here we Integral approach for time dependent materials... 1033 will describe the solution process of a step-by-step nature in which the loads are applied stepwise at various time intervals. According to equation (2.2), the strain at time tn and tn+1 may be written as εij(tn)= 3∑ k=1 3∑ l=1 tn∫ −∞ Jijkl(tn− τ) ∂σkl(τ) ∂τ dτ εij(tn+1)= 3∑ k=1 3∑ l=1 tn+1∫ −∞ Jijkl(tn+1− τ) ∂σkl(τ) ∂τ dτ (2.3) In order to establish the incremental viscoelastic formulation, equation (2.3)2 is separated into a sum as follows εij(tn+1)= 3∑ k=1 3∑ l=1 tn∫ −∞ Jijkl(tn+1− τ) ∂σkl(τ) ∂τ dτ + 3∑ k=1 3∑ l=1 tn+1∫ tn Jijkl(tn+1− τ) ∂σkl(τ) ∂τ dτ (2.4) The first integral represents the hereditary response, while the second integral deals with the implicit pseudo instantaneous response of the material. Using equations (2.3) and (2.4), the following incremental viscoelastic equation is obtained ∆εij(tn)= 3∑ k=1 3∑ l=1 tn∫ −∞ [Jijkl(tn+1− τ)−Jijkl(tn− τ)] ∂σkl(τ) ∂τ dτ + 3∑ k=1 3∑ l=1 tn+1∫ tn Jijkl(tn+1− τ) ∂σkl(τ) ∂τ dτ (2.5) in which ∆εij(tn) = εij(tn+1)− εij(tn) is the total strain increment during the time step ∆tn. In order to simplify the incremental viscoelastic law, let us define the influence of the complete past history of stress by the hereditary integral as follows ε̃ij(tn)= 3∑ k=1 3∑ l=1 tn∫ −∞ [Jijkl(tn+∆tn− τ)−Jijkl(tn− τ)] ∂σkl(τ) ∂τ dτ (2.6) 1034 C. Chazal, R. Moutou Pitti ε̃ij(tn) is the memory term which involves the whole past solutions. Insertion of equation (2.6) into equation (2.5) results in ∆εij(tn)= 3∑ k=1 3∑ l=1 tn+1∫ tn Jijkl(tn+1−τ) ∂σkl(τ) ∂τ dτ+ ε̃ij(tn) (2.7) In the next Section, the integral in this equation is discretized in order tomake it completely incremental. 3. Incremental constitutive equations The incremental viscoelastic equations will be derived using a linear approxi- mation of the stress. In fact, we assume that the time derivative during each time increment is constant, a staircase function. ∀τ ∈ [tn, tn+∆tn] ⇒ σkl(τ)=σkl(tn)+ τ − tn ∆tn ∆σkl(tn)H(τ − tn) (3.1) in which ∆σkl(tn) = σkl(tn+1)−σkl(tn) is the total stress increment during the time step ∆tn. Substituting into equation (2.7) the linear approximation of the stress in- crement given by equation (3.1) and using the exponential approximation for the creep compliance components given by equation (2.1), we obtain ∆εij(tn)= 3∑ k=1 3∑ l=1 ∆σkl(tn) ∆tn · tn+∆tn∫ tn { J (0) ijkl+ M∑ m=1 J (m) ijkl [ 1− e −λ (m) ijkl (tn+∆tn−τ) ]} dτ + ε̃ij(tn) (3.2) This equation may be integrated in closed form to produce ∆εij(tn)= 3∑ k=1 3∑ l=1 Πijkl(∆tn)∆σkl(tn)+ ε̃ij(tn) (3.3) Πijkl(∆tn) can be interpreted as the viscoelastic compliance tensor and is given by Πijkl(∆tn)= J (0) ijkl+ M∑ m=1 J (m) ijkl [ 1− 1 ∆tnλ (m) ijkl ( 1− e −λ (m) ijkl ∆tn )] (3.4) Integral approach for time dependent materials... 1035 The memory term ε̃ij(tn) in equation (3.3), which involves the whole past solutions, is given by equation (2.6) and defines the influence of the complete past history of stress by a Volterra hereditary integral. Our purpose is now the conversion of this hereditary integral to a more convenient form in order to integrate it in a finite element solution.We now introduceWiechert model (2.1) onto the formulation by way of substitution into (2.6); doing so yields ε̃ij(tn)= 3∑ k=1 3∑ l=1 M∑ m=1 ( 1−e −λ (m) ijkl ∆tn ) tn∫ −∞ J (m) ijkle −λ (m) ijkl (tn−τ)∂σkl(τ) ∂τ dτ (3.5) This may also be written as ε̃ij(tn)= 3∑ k=1 3∑ l=1 Ψijkl(tn) (3.6) where Ψijkl(tn)= M∑ m=1 ( 1− e −λ (m) ijkl ∆tn ) Φ (m) ijkl(tn) (3.7) In the above, Ψijkl(tn) is a fourth-order viscoelastic tensor which represents the influence of thewhole past history of stress and Φ (m) ijkl(tn) is a pseudovisco- elastic tensor corresponding to one of the spectrumrange of the decomposition given by equation (2.1). It is given by Φ (m) ijkl(tn)= tn∫ −∞ J (m) ijkle −λ (m) ijkl (tn−τ)∂σkl(τ) ∂τ dτ (3.8) Let us divide the domain of integration −∞¬ τ ¬ tn in equation (3.8) into two parts: −∞¬ τ ¬ tn−∆tn and tn−∆tn ¬ τ ¬ tn. This leads to Φ (m) ijkl(tn)= tn−∆tn∫ −∞ J (m) ijkle −λ (m) ijkl (tn−∆tn+∆tn−τ)∂σkl(τ) ∂τ dτ + tn∫ tn−∆tn J (m) ijkle −λ (m) ijkl (tn−τ)∂σkl(τ) ∂τ dτ (3.9) Weassumed that thepartial derivative ∂σkl(τ)/∂τ appearing in equation (3.8) can be approximated as ∆σkl/∆tn. The last equation may be integrated in closed form to produce 1036 C. Chazal, R. Moutou Pitti Φ (m) ijkl(tn)=Φijkl(tn−1)e −λ (m) ijkl ∆tn + ∆σkl(tn−1) ∆tn−1 tn∫ tn−∆tn J (m) ijkle −λ (m) ijkl (tn−τ)∂σkl(τ) ∂τ dτ (3.10) The pseudo viscoelastic tensor at time tn may now be expressed as Φ (m) ijkl(tn)=Φ (m) ijkl(tn−1)e −λ (m) ijkl ∆tn + ∆σkl(tn−1)J (m) ijkl ∆tn−1λ (m) ijkl [ 1− e −λ (m) ijkl ∆tn ] (3.11) Finally, the incremental constitutive law given by equation (3.3) can now be inverted to obtain ∆σij(tn)= 3∑ k=1 3∑ l=1 Ξijkl(∆tn)∆εkl(tn)− σ̃ij(tn) (3.12) where Ξijkl =(Πijkl) −1 is the inverse of the compliance tensor and σ̃ij(tn) is a pseudo stress tensorwhich represent the influence of the complete past history of strain. It is given by σ̃ij(tn)= 3∑ k=1 3∑ l=1 [Πijkl(∆tn)] −1ε̃ij(tn) (3.13) The incremental constitutive law represented by equation (3.13) is introdu- ced in a finite element discretisation in order to obtain solutions to complex viscoelastic problems. 4. Finite element discretization The finite element method is used to implement the proposed approach. This allows us to resolve complex viscoelastic problems with real boundary condi- tions. 4.1. Virtual displacement principle in linear viscoelasticity A standard finite element procedure is used as explained in Zienkiewicz (1958) andGhazlan et al. (1995). Incremental equilibriumequations for the vi- scoelastic problemare obtained through theprinciple of virtualwork (Ghazlan Integral approach for time dependent materials... 1037 et al., 1995). For a three dimensional continuumproblem, it yields for any ele- ment ∫ V e σij(tn+1)δεij(tn+1) dV = ∫ V e fvi (tn+1)δui(tn+1) dV + ∫ ∂V e Ti(tn+1)δui(tn+1) dS (4.1) Substituting incremental relations σij(tn+1)=σij(tn)+∆σij εij(tn+1)= εij(tn)+∆εij ui(tn+1)=ui(tn)+∆ui (4.2) into Eq. (4.1) and observing that δεij(tn+1) = δ(∆εij), δui(tn+1) = δ(∆ui) yield ∫ V e ∆σijδ(∆εij) dV = ∫ V e fvi (tn+1)δ(∆ui) dV + ∫ ∂V e Ti(tn+1)δ(∆ui) dS− ∫ V e σij(tn)δ(∆εij) dV (4.3) Note that the virtual workmade by load increments during the time step ∆tn is given by the difference between the twofirst terms and the third termon the right-handsideof equation (4.3). Introducing incremental viscoelastic equation (3.12) into the first hand side of equation (4.3), the following expression can be obtained ∫ V e Ξijkl(∆tn)∆εklδ(∆εij) dV = ∫ V e fvi (tn+1)δ(∆ui) dV (4.4) + ∫ ∂V e Ti(tn+1)δ(∆ui) dS− ∫ V e σij(tn)δ(∆εij) dV + ∫ V e σ̃ij(tn)δ(∆εij) dV where ∆ui is the incremental displacement field between tn and tn+1, fvi (tn+1) are the body forces per unit volume, Ti(tn+1) is the surface tensions perunit surface, δ is thevariation symboland V e is thevolumeof the element. Assuming small displacements, the strains are derived from shape functions using a standard manner in the context of the finite element method. Using matrix notation, the strain increment can be written as {∆εij}= [D]{∆U e} (4.5) where ∆Ue is the local element displacement increment and [D] is the strain- displacement transformation matrix. 1038 C. Chazal, R. Moutou Pitti 4.2. Viscoelastic stiffness matrix computation The equilibrium equations, in matrix notation, can be simplified by intro- ducing the approximation of the strain components given by equation (4.5) into equation (4.4). This leads to ∫ V e [D][δ(∆Ue)]⊤[Θ(∆tn)][D][∆U e] dV = ∫ V e [δ(∆[Ue)]⊤{fvi (tn+1)} dV + ∫ ∂V e [δ(∆Ue)]⊤{Ti(tn+1)} dS− ∫ V e [D][δ(∆Ue)]⊤σij dV + ∫ V e [D][δ(∆Ue)]⊤σ̃ij(tn) dV (4.6) Interpolation functionsmust be used in order to obtain equilibrium equations in discrete form. The displacement field ∆Ue is approximated by {∆Ue}= [Ne]{∆qe} (4.7) where [Ne] is thematrix of shape functions and {∆qe} are the nodal element displacement increment. Substituting equation (4.7) into equation (4.6), we get ∫ V e [BL] ⊤[Θ(∆tn)][BL]{∆q e} dV = ∫ V e [Ne]⊤{fv(tn+1)} dV (4.8) + ∫ ∂V e [Ne]⊤{T(tn+1)} dS− ∫ V e [BL] ⊤{σ(tn)} dV + ∫ V e [BL] ⊤{σ̃(tn)} dV where [BL] = [D][N e] is the total strain-displacement transformation matrix. These equilibrium equations for linear viscoelastic behavior can be rewritten as [KT(∆tn)]{∆q e}= {Fv(tn+1)}+{F s(tn+1)}−{F σ(tn)}+{F vis(tn)} (4.9) Difference between {Fv(tn+1)}+ {F s(tn+1)} and {F σ(tn)} is the external load increment during time step ∆tn while {F vis(tn)} is the viscous load vector corresponding to the complete past history. Equation (4.9) can then be simplified to obtain [KT(∆tn)]{∆q e}= {∆Fext(tn)}+{F vis(tn)} (4.10) Integral approach for time dependent materials... 1039 where the stiffness matrix [KT(∆tn)] is given by [KT∆(tn)] = ∫ V e [BL] ⊤[Θ(∆tn)][BL] dV (4.11) and the nodal force vector is denoted as the sum of two vectors: {∆Fext(tn)} is the external load increment vector, {Fvis(tn)} is the memory load vector and [Θ(∆(tn)] is the viscoelastic constitutive matrix. The viscous load vector increment is given by {Fvis(tn)}= ∫ V e [BL] ⊤{σ̃(tn)} dV (4.12) 4.3. Incremental viscoelastic algorithm The formulation is introduced in the software Cast3m used by the French Energy Atomic Agency. The software can be employed for plane linear visco- elasticity structures. The global incremental procedure for the creep integral approach is described as below: 1. The instantaneous response (elastic solution) to the applied load is first obtained. Thus at the time t= t0, the stress σ0, the displacements q0 and the strains ε0 are all known. 2. Given a solution at the time tn. That is, given the stresses σ(tn) and displacements q(tn) that satisfy the equilibriumequations, our objective is to determine the solution at the time tn+1 = tn+∆tn. Loop over the time interval of study (a) compute the compliance moduli Πijkl(∆tn) from equation (3.4) and the tangent moduli Ξijkl(∆tn) from the relation Ξijkl(∆tn)= = [Πijkl(∆tn)] −1 and then get the viscoelastic constitutive matrix [Θ(∆(tn)]= [Ξijkl(∆tn)] (b) determine the pseudo stress tensor σ̃ij(tn) from equation (3.13) and then compute the viscous load vector {Fvis(tn)} fromequation (4.12) (c) update the linear viscoelastic stiffnessmatrix [KT(∆tn)] from equ- ation (4.11) (d) assemble and solve the system of equilibrium equations (4.10) to obtain the displacement increment ∆q(∆tn) (e) compute the strain increment ∆εe(∆tn) from equations (4.10) (f) use the result of step (b) to compute the stress increment ∆σe(∆tn) from equation (3.12) 1040 C. Chazal, R. Moutou Pitti (g) evaluate the pseudo viscoelastic tensor Φ (m) ijkl(tn) at time tn+1 from equation (3.11) (h) compute thememory term ε̃ij(tn+1) from equations (3.6) and (3.7) (i) update the state (displacement, stress and strain) at the end of the time increment ∆tn: {q(tn+1)}= {q(tn)}+{∆q(∆tn)} {σij(tn+1)}= {σij(tn)}+{∆σij(∆tn)} {εij(tn+1)}= {εij(tn)}+{∆εij(∆tn)} (j) go to step (a) 5. Numerical applications 5.1. Viscoelastic plane stress panel This example is used to check the validity of the creep incremental con- stitutive formulation proposed in this work. The structure analyzed is a plane stress panel subjected to compressive load in the x direction and transverse load in the y direction distributed along the length L. The panel ismade of a Fig. 1. Viscoelastic panel submitted to axial and transverse load homogeneous and isotropic viscoelastic material with dimensions specified in Fig.1 and is pinned at the ends. Horizontal and vertical loads are applied: • Transverse load in the y direction: py(t)= 32kN/m ∀t­ 0 • Compressive load in the x direction (t in hours): σx(t)= { σ0H(t) t¬ 400 and t­ 600 0 400