- Email: [email protected]
A time-integration algorithm for thermo-rheologically complex polymers
Available online at www.sciencedirect.com
Computational Materials Science 41 (2008) 576–588 www.elsevier.com/locate/commatsci
A time-integration algorithm for thermo-rheologically complex polymers Anastasia Muliana *, Kamran A. Khan Department of Mechanical Engineering, Texas A&M University, College Station, TX 77843-3123, United States Received 30 November 2006; received in revised form 28 May 2007; accepted 30 May 2007 Available online 20 July 2007
Abstract This study presents nonlinear thermo-viscoelastic analyses of isotropic polymers that belong to a class of thermo-rheologically complex materials (TCM). The nonlinear viscoelastic constitutive model is expressed with an integral form of a creep function, whose initial, long-term, and transient properties change with stresses and temperatures. A combined recursive–iterative method is formulated to solve the viscoelastic integral equation. The recurrence formula allows bypassing the need to store entire strain histories. Two types of iterative procedures, fixed point (FP) and Newton–Raphson (NR), are examined within the recursive algorithm. Furthermore, a consistent tangent stiffness matrix is formulated to accelerate convergence and avoid divergence. The algorithm is compatible with displacement based finite element (FE) structural analyses for small deformation gradients and uncoupled thermo-mechanical problems. Verification of the numerical algorithm is performed using the thermo-viscoelastic experimental data available in the literature. The influences of temperature dependent material properties on the material’s viscoelastic responses are discussed for general thermo-mechanical loadings. The capability of the algorithm in predicting multi-axial viscoelastic responses is then verified using creep data of adhesive bonded joints at various temperatures. Finally, numerical simulations of time-dependent crack propagations in adhesive bonded joints are presented for the opening and shearing modes. 2007 Elsevier B.V. All rights reserved. Keywords: Viscoelastic; Thermo-mechanical; Polymers; Finite element; Recursive; Iterative
1. Introduction Polymers exhibit linear and nonlinear viscoelastic behaviors under combined mechanical and environmental loadings. The nonlinear viscoelastic responses, which often occur at high load levels and elevated temperatures, are indicated by non-constant material properties. Based on their rheological responses at different temperatures, polymers can be classified as thermo-rheologically simple materials (TSM) or thermo-rheologically complex materials (TCM). If temperature affects mainly the time-dependent material properties and the effect of temperature can be incorporated through a time-scale shift factor only, materials are categorized as TSM. On the other hand, for *
Corresponding author. Tel.: +1 979 458 3579; fax: +1 979 845 3081. E-mail address: [email protected] (A. Muliana).
0927-0256/$ - see front matter 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2007.05.021
the TCM, temperature influences the material’s initial (elastic), long-term (equilibrium), and time-dependent (transient) properties. Further discussion regarding the TSM and the TCM can be found in Harper and Weitsman [13], Caruthers and Cohen [4], Caruthers et al. [5], Wineman and Rajagopal [33], Tschoegl et al. [31]. Numerical algorithms that are compatible with finite element (FE) analyses have been developed mainly for analyzing viscoelastic constitutive models of the TSM. Zienkiewicz et al. [36] and Bazant [2] presented incremental algorithm for simulating viscoelastic behaviors derived directly from a governing differential equation of a spring-dashpot mechanical model. The strain rate was modeled as functions of the current stresses and the accumulated strains. Other types of constitutive viscoelastic models, i.e. hereditary integral forms, have been formulated for predicting time-dependent material responses
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
under more general loading histories. A recurrence numerical algorithm was first used by Taylor et al. [30] for solving a linear viscoelastic integral equation. The recursive approach allows formulating current stress tensor as a function of current strain increments and history variables stored in the previous time step. This approach bypasses the need to store entire history variables, which minimizes the storage required in performing the integration. Similar recursive approach was used by Feng [8] for linear viscoelastic materials described by the Voltera integral equation. The recursive method has also been used for solving nonlinear viscoelastic constitutive equations of the TSM typed. Henriksen [10], Roy and Reddy [25], and Lai and Bakker [17] developed recursive integration algorithms for solving the Schapery’s [27] nonlinear constitutive model of isotropic materials under combined mechanical and environmental loadings. A Prony exponential series was used in the transient compliance. The computer storage requirement for using this method was proportional to the number of terms in the Prony series. The effects of temperature, moisture, and aging were carried through the time-shift factors only. Haj-Ali and Muliana [9] formulated a recursive– iterative algorithm of the Schapery stress-dependent viscoelastic constitutive models. An iterative scheme based on the Newton–Raphson method was added for correcting the trial linearized stresses and minimizing the residual strains. Simo and Hughes [29] presented several numerical algorithms of inelastic and isotropic materials with internal state variables including viscoelasticity. The hereditary integral models have been extended for anisotropic viscoelastic materials. Poisson’s effect is often assumed to be constant. Zocher et al. [37], Poon and Ahmad [21], Poon and Ahmad [22], Hilton and Yi [11], and Yi et al. [35] developed a FE scheme for linear and nonlinear viscoelastic materials of anisotropic media using integration point update algorithm. To enhance convergence rate, consistent tangent stiffness matrix was formulated. Idealizing thermo-viscoelastic responses of polymers as TSM is sufficient when the materials are subjected to moderate temperature changes and temperature does not vary with time or when the materials are under constant temperature rates and constant mechanical loads. There are many applications when the deformed/stressed bodies are subjected to transient temperatures or when thermo-mechanical coupling due to internal energy dissipation is prominent. Holzapfel and Reiter [12], Reese and Wriggers [24], Johnson and Chen [14] presented time-integration procedures for coupled thermo-viscoelastic behaviors for small and large deformation problems, in which the materials deviate from the TSM. Constitutive equations were derived from the Maxwell mechanical analog model with constant material properties. Reese and Wriggers [24] and Johnson and Chen [14] have pointed out that the constitutive equations should incorporate material degradation due to the temperature changes to accurately predict the thermoviscoelastic behaviors especially when the thermo-viscoelastic coupling is not negligible, however, constitutive
577
models for thermo-viscoelastic responses of the TCM are still lacking. Major studies on the thermo-viscoelastic behaviors of polymers have been focused on a class of TSM. There are many engineering applications that involve transient temperatures and a thermo-mechanical coupling, under which the viscoelastic responses of a materials deviate from the TSM. Thus, there is a need to develop a constitutive equation for the TCM that can be integrated with general structural analysis tools. This study presents a time-integration algorithm for predicting thermo-viscoelastic behaviors of isotropic polymers that belong to a class of TCM. This algorithm is designed to be compatible with the FE framework to perform large scale structural analyses under general thermo-mechanical loadings. For this purpose, the effects of stress and temperature should be incorporated in the initial, long-term, and time-dependent material properties. The theoretical basis of the studied viscoelastic constitutive equation is given in Harper and Weitsman [13]. In Section 2, a recursive–iterative algorithm is formulated. Two types of iterative procedures, fixed point (FP) and Newton–Raphson (NR), are examined within the recursive scheme. Furthermore, a consistent tangent stiffness matrix is formulated at the material level to accelerate convergence and avoid divergence at the structural level. Verifications of the numerical algorithm with thermo-mechanical viscoelastic experimental data available in the literature are presented in Section 3. Convergence behaviors using the FP and NR methods are examined for several thermo-mechanical loading histories. In addition, performances of the TCM under general thermo-viscoelastic behaviors are discussed based on the numerical results. In Section 4, examples of the integrated material model with the FE analyses are given on time-dependent and crack propagation analyses of polymeric adhesive bonded joints. Creep data on the adhesive bonded lap joints at different temperatures and humidity are used for comparisons. Numerical predictions of a crack initiation and a complete debonding time for mode I (opening) and II (shearing) failures are also presented. 2. Formulation of the time-integration algorithm Schapery [28] and Sadkin and Aboudi [26] have derived constitutive material models for the TCM that involves two time-scales for incorporating temperature and stress variations. The use of two time-scales in the constitutive material models makes development of numerical algorithm very challenging and is often difficult to integrate these constitutive material models in general FE structural analyses. In this study, the time–stress–temperature effects on the material behaviors are linked to one time-scale without any restriction for modeling temperature or stress variations with time. The Schapery [27] nonlinear single integral constitutive model is modified to include the effects of stress and temperature on the initial, long-term, and time-dependent material properties for non-aging materials, which is given as:
578
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
Z t t s et eðtÞ ¼ g0 ðrt ; T t ÞD0 rt þ g1 ðrt ; T t Þ DDðw w Þ 0 Z t s d½g0 ðrs ; T s Þrs s dDT ds þ ds aðT Þ ds ds 0
ð1Þ
Here D0 and DD are the uniaxial instantaneous elastic and transient compliances. The uniaxial transient compliance is expressed as: t
DDw ¼
N X
Dn ð1 exp½kn wt Þ
ð2Þ
n¼1
w is the reduced-time (effective time) given by: Z t Z s dn dn t s w wðsÞ ¼ w wðtÞ ¼ n n n n 0 a r ;T 0 aðr ; T Þ
where D0 ¼
1 E0
Dð1Þ ¼ g0 ðrt ; T t ÞD0 þ g1 ðrt ; T t Þg2 ðrt ; T t Þ
ð4Þ N X
Dn
n¼1
1 where Dn ¼ En
ð5Þ
The uniaxial viscoelastic behavior in Eq. (1) is generalized for multi-axial (3D) constitutive relations of isotropic materials by separating the deviatoric and volumetric strain–stress relations as: 1 etij ¼ etij þ etkk dij þ aðT t T 0 Þdij 3 1 t eij ¼ g0 ð rt ; T t ÞJ 0 S tij 2 Z t o½g2 ð rs ; T s ÞS sij 1 wt ws Þ t t ð þ g1 ð r ;T Þ DJ ds 2 os 0
ð6Þ
The parameters J0 and B0 are the instantaneous elastic shear and bulk compliances, respectively. The terms DJ and DB are the time-dependent shear and bulk compliances, respectively. The nonlinear parameters of the multi-axial behaviors are modeled as a function of current t and temperature Tt. The matrix Poisson’s effective stress r ratio, m, is assumed to be time-independent, which allows expressing the shear and bulk compliances as: J 0 ¼ 2ð1 þ mÞD0
ð3Þ
The constitutive model in Eq. (1) can be visualized as series arrangement of N Kelvin–Voigt elements and an elastic spring. Each spring has stress and temperature dependent material constants. The parameters g0 and g1 measure changes in the elastic and the transient compliances, respectively, with stresses and temperatures. The parameter g2 accounts for the loading rate effect on the time-dependent response. The parameter aðrt ; T t Þ is a time-shift (interchange) factor measured with respect to the reference stress and temperature. If the thermal effects are carried through the shift factor a(Tt) only, the thermo-rheologically simple material (TSM) is exhibited. The parameters T and T0 are the current and the reference temperatures, respectively. The superscript denotes the time-dependent variable. In general, the coefficient of thermal expansion a depends on the temperature; in this study a is assumed to be constant. The initial and the long-term nonlinear creep compliances can be determined from a step stress history by setting time t = 0 and 1, respectively, which are expressed by: Dð0Þ ¼ g0 ðrt ; T t ÞD0
1 rt ; T t ÞB0 rtkk etkk ¼ g0 ð 3 Z t 1 t s o½g 2 ð rs ; T s Þrskk þ g1 ð rt ; T t Þ DBðw w Þ ds 3 os 0
DJ
wt
B0 ¼ 3ð1 2mÞD0 t
¼ 2ð1 þ mÞDDw
t
t
DBw ¼ 3ð1 2mÞDDw
ð7Þ
For simplicity, the nonlinear material parameters will be written as: gb ð rt ; T t Þ ¼ gtb að rt ; T t Þ ¼ at
b ¼ 0; 1; 2 rffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 t t t ¼ S S r 2 ij ij
ð8Þ
2.1. Recursive formulation A recursive method is used to solve the integrals in Eq. (6). The algorithm is designed to be compatible with general displacement based FE structural analyses for uncoupled thermo-mechanical loadings. The numerical algorithm of the TCM will be called at every material (Gaussian) integration point within elements at each structural iteration to achieve structural and material convergences simultaneously. Thus, an efficient and accurate numerical algorithm for solving the constitutive material model becomes necessary. Using the transient compliance in Eq. (2) and the relations in Eq. (7), the strains in Eq. (6) can be written as: N N X 1 1 1 X etij ¼ gt0 J 0 S tij þ gt1 S tij J n gt1 J n qtij;n 2 2 2 n¼1 n¼1 Z t o gs2 S sij where qtij;n ¼ ds exp½kn ðwt ws Þ os 0 N N X 1 1 1 X etkk ¼ gt0 B0 rtkk þ gt1 rtkk Bn gt1 Bn qtkk;n 3 3 3 n¼1 n¼1 s s Z t t s o g 2 rkk t where qkk;n ¼ ds exp½kn ðw w Þ os 0
ð9Þ
ð10Þ
A recursive integration method, which was derived by Taylor et al. [30] for a linear viscoelastic integral model and extended by Haj-Ali and Muliana [9] for a nonlinear stress-dependent viscoelastic constitutive model, is used to solve the integral parts in Eqs. (9) and (10). These are expressed as:
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
qtij;n
¼
Z
tDt s
t
exp½kn ðw w Þ
o gs2 S sij
Detij ¼ etij eijtDt ¼ Jt S tij JtDt S tDt ij
ds os Z t o gs2 S sij ds þ exp½kn ðwt ws Þ os tDt t t tDt tDt ¼ exp½kn Dwt qtDt S ij ij;n þ g2 S ij g2 0
qtkk;n
1 exp½kn Dwt kn Dwt
ð11Þ
tDt
t
1 exp½kn Dwt kn Dwt
s
ð12Þ
The first integral includes the limits (0, t Dt), i.e. up to the previous converged time step, which are stored as history variables of the hereditary strain integrals. The limits of the second integral are taken as (t Dt, t), which is the current incremental step. Detailed derivation of the recursive– iterative algorithm for the nonlinear viscoelastic model can tDt be found in Haj-Ali and Muliana [9]. The parameters qij;n tDt and qkk;n , n = 1, . . ., N are the hereditary integral for every term in the Prony series in the form of deviatoric and volumetric strains. These parameters are history state variables stored from the last converged step at time (t Dt). The total number of history variables is 7 · N · NGauss; where 7 is the total number of deviatoric and volumetric components, N is the number of terms in the Prony series, and NGauss is the number of Gaussian integration points. The incremental reduced-time is expressed as Dwt wt wtDt. Substituting the integral Eqs. (11) and (12) into Eqs. (9) and (10) gives the deviatoric and volumetric total strains: etij
" # N N X X 1 t 1 exp ½kn Dwt t t t t t g J 0 þ g1 g2 ¼ J n g1 g 2 Jn S ij 2 0 kn Dwt n¼1 n¼1
N t 1 X tDt 1 exp½kn Dw tDt gt1 J n exp½kn Dwt qtDt S ij t ij;n g 2 2 n¼1 kn Dw
Jt S tij d tij ð13Þ " # N N X X 1 t 1 exp½kn Dwt t rkk etkk ¼ Bn gt1 gt2 Bn g0 B0 þ gt1 gt2 3 kn Dwt n¼1 n¼1
N t 1 X tDt 1 exp½kn Dw tDt gt1 Bn exp ½kn Dwt qtDt g r kk;n 2 kk 3 n¼1 kn Dwt t rt V t ð14Þ B kk
N tDt 1X J n gt1 exp ½kn Dwt g1tDt qij;n 2 n¼1 " ! tDt N 1 tDt X tDt 1 exp kn Dw g2 J n g1 2 kn DwtDt n¼1 # t t 1 exp ½kn Dw g1 S ijtDt kn Dwt
o gs2 rskk ds ¼ exp½kn ðw w Þ os 0 Z t oðgs2 rskk Þ ds þ exp½kn ðwt ws Þ os tDt t t tDt tDt ¼ exp½kn Dwt qtDt rkk kk;n þ g2 rkk g 2 Z
579
kk
The parameters J t and bulk Bt are the current time-dependent deviatoric and volumetric compliances, respectively, while d tij and vtkk are the deviatoric and volumetric strains that carry the history variables. The incremental forms of the deviatoric and volumetric strains are:
t rt B tDt rtDt ¼B Detkk ¼ etkk etDt kk kk kk N t 1X Bn g1 exp ½kn Dwt g1tDt qtDt kk;n 3 n¼1 " ! tDt N 1 tDt X tDt 1 exp kn Dw g2 Bn g 1 3 kn DwtDt n¼1 # t t 1 exp ½kn Dw g1 rtDt kk kn Dwt
ð15Þ
ð16Þ
Eqs. (13)–(16) define complete solutions for the total and incremental strain tensors, respectively. The nonlinear parameters in Eqs. (13)–(16) are expressed as functions of the current temperature and total stresses; however, the total stresses at the current time (t) are not known. Linearized trial stresses are used as a starting point to obtain the current stresses using either Eqs. (13) and (14) or (15) and (16). An iterative scheme is added in order to find the correct stress tensor. The chosen trial values (initial approximation solutions) can significantly determine the convergence solutions. Poor initial approximation values often lead to divergence. Two iterative methods, which are fixed point (FP) and Newton–Raphson’s (NR) are used in this study. The FP method gives a slower but stable convergence rate and requires less computational effort during each iteration step. The NR method can increase convergence process rapidly, however, linearized tangent functions often lead to divergence solutions if during the iteration process the convergence is not monotonic or if the magnitude of the tangent functions is very small (approaching zero). 2.2. Fixed point (FP) method The iteration procedure for solving a variable x that satisfies x = f(x) in the FP method is summarized as: xðkþ1Þ ¼ f ðxðkÞ Þ R ¼ kxðkþ1Þ xðkÞ k Converged solutions:
ð17Þ R!0
where k indicates the iteration counter with an initial approximation solution x(0) and R is the measured error. The current stress tensor at the iteration k from a given mechanical strain is: i 1 1 1h rtij ¼ S tij þ rtkk dij ¼ t etij þ d tij þ t etkk þ vtkk dij ð18Þ 3 J B
580
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
The effective shear J t and bulk Bt compliances, and the nonlinear parameters in the Eq. (18) are functions of the t and temperatures Tt. The curcurrent effective stresses r rent temperature is defined by Tt = TtDt + DTt. The initial approximation (trial) stress tensor is determined using the following approximation of nonlinear parameters: gt;tr b a
t;tr
¼
t;ð0Þ gb
¼a
t;ð0Þ
¼ gb ð r ¼ að r
tDt
tDt
;T
;T
tDt
tDt
Þ b ¼ 0; 1; 2
Þ
ð19Þ
The trial stress tensor is formed based on the given (known) variables and history variables: t;ð0Þ
rt;tr ij ¼ rij
t;ð0Þ
¼ rijtDt þ Drij t;ð0Þ
¼ rijtDt þ DS ij
1 t;ð0Þ þ Drkk dij 3
ð20Þ
where the trial incremental deviatoric and volumetric stresses are given by: " # N X 1 1 t;ð0Þ t;ð0Þ t tDt DS ij ¼ t;ð0Þ Detij þ g1 J n ðexp½kn Dw 1Þqij;n 2 J n¼1 " # N 1 1 t;ð0Þ X t;ð0Þ t t tDt Drkk ¼ t;ð0Þ Dekk þ g1 Bn ðexp½kn Dw 1Þqkk;n 3 B
defined by using either the incremental strains or the total strains. The residual strain tensor is expressed by: t;ðkÞ
Rij
t;ðkÞ
¼ Deij
1 t;ðkÞ þ Dekk dij D^etij 3
ð25Þ
The incremental deviatoric and volumetric strain tensors at the iteration k in the material level are defined using Eqs. (15) and (16). The current incremental mechanical strains t;ðmÞ D^etij D^eij is obtained from the structural level at the mth global iteration. The initial approximation (trial) stress tensor is determined using Eqs. (19)–(21). A Jacobian tensor is formed by taking the derivative of the residual tensor in t Eq. (25) with respect to the incremental stress tensor oRij . Finally, the consistent tangent stiffness matrix C tijkl oDrtkl at the converged state is: t 1 oDrtij oRij t C ijkl ¼ ; kRtij k ! 0 ð26Þ oD^etkl oDrtkl Once convergence has been achieved, the stress, the hereditary integrals for each Prony series, and the consistent tangent stiffness are updated. The recursive–iterative algorithm is summarized in Fig. 1.
n¼1
ð21Þ During the iteration, the component of the residual tensor is expressed in terms of the total stress: t;ðkþ1Þ Rij
¼
t;ðkþ1Þ rij
t;ðkÞ rij
ð22Þ
The nonlinear parameters in Eq. (18) is defined by t;ðkÞ ð t;ðkÞ ; T t Þ and at;ðkÞ ¼ að gb ¼ g b r rt;ðkÞ ; T t Þ. Once convergence is achieved, the consistent tangent stiffness matrix in Eq. (23), the current stress tensor, and the history variables in Eqs. (11) and (12) are updated and sent to the structural level. In the coupled thermo-mechanical analyort ses, it is also necessary to determine oTijt .
ortij 1 1 1 1 t ð23Þ C ijkl ¼ t ¼ t dik djk þ dij 3 Bt J t o^ekl J The variable ^etij is the component of the mechanical strain tensor. 2.3. Newton–Raphson’s (NR) method The iterative procedure using the NR’s method for solving a variable x that satisfies F(x) = 0 is summarized as: xðkþ1Þ ¼ xðkÞ
F ðxðkÞ Þ F 0 ðxðkÞ Þ
dF ðxÞ dx Converged solutions:
ð24Þ
F 0 ðxÞ ¼
kRk ¼ kF k ! 0
This procedure requires that F(x) and F 0 (x) are continuous near the solution x and the magnitude of F 0 (x) should not approach zero or infinite. The residual tensor can be
2.4. Implementation of the recursive–iterative algorithm in the FE framework The constitutive thermo-viscoelastic material model with linearized and corrector schemes (Fig. 1) is implemented in the user material (UMAT) ABAQUS FE (2005). At the FE structural level, an iterative equation solution is also performed. During this structural iteration the recursive–iterative algorithm will be called at each material (Gaussian) point. In every iteration within an incremental time step Dt(m), trial incremental strain tensor t;ðmÞ Deij and temperature DTt,(m) are prescribed, which are the input variables to the subroutine UMAT. The superscript (m) denotes the global iteration counter within the current incremental time step. The output from the UMAT t;ðmÞ is the total stresses rij and the consistent tangent stiffness t;ðmÞ t;ðMÞ C ijkl . The converged C ijkl after M global iteration at the current time t will be used to provide incremental trial strains for the next time step (t + Dt). In the uncoupled thermo-mechanical analyses, the trial incremental temperature is directly linked to the incremental time step. By adding the corrector scheme at the material level, it is not necessary to have a tight incremental time Dt to achieve convergence at the structural level. In fact, only by tightening Dt without adding a corrector scheme at the material level can lead to divergence, which was discussed in the previous studies by Haj-Ali and Muliana [9] and Muliana and Kim [18]. At the material level, the convergence criterion can be defined using the residual vectors in Eqs. (22) or (25), which are expressed in terms of stress and strain vectors, respectively. In this study, the given tolerances allow for the maximum residual stress to be 0.01% of the total stress and the maximum residual strain to be 1 micro-
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
581
Fig. 1. Recursive–iterative algorithm for the nonlinear isotropic viscoelastic material.
strain. The above tolerances are set as the limits since reducing the tolerances below these values does not change the overall responses within our numerical values of interest, while increasing the tolerance beyond these values can lead to divergence in the overall responses. The users should also realize that the stress vector at the local (material) iteration (k + 1) in the NR iterative method is defined as an intersection between the stress vector and the tangent oRt to the residual vector oDrijt at the iteration k. Thus, if at the kl
iteration k the value of
oRtij oDrtkl
is approaching zero or infinite
convergence at the material level will not be achieved, which indicates by either an extremely large stress correction or a zero correction. In this situation, adding convergence criteria such as a maximum number of iteration and a minimum correction rate can be used to terminate the
iteration at the material level and a smaller time increment Dt can be chosen to possibly achieve a convergence.
3. Verification of the thermo-mechanical viscoelastic constitutive model The proposed numerical algorithm of the TCM is verified using the thermo-viscoelastic experimental data available in the literature. Test data on FM73 adhesive polymers reported by Peretz and Weitsman [19], Peretz and Weitsman [20] and epoxy Hercules 3502 of Harper and Weitsman [13] are used. The algorithm is implemented in the 3D continuum elements. The 3D nonlinear viscoelastic material response can also be integrated with other typed of elements e.g: 1D (truss, beam) or 2D (plane stress,
582
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
plane strain) by imposing uniaxial or in-plane stress–strain relations to the 3D constitutive model.
where T0 is the reference temperature, which is 303 K; the stress r and temperature T are in MPa and Kelvin, respectively. Convergence behaviors using the FP and NR methods are presented first. Fig. 2 illustrates the nonlinear elastic responses under a linear stress and temperature ramp loading. Iteration processes are performed both at the material and structural levels. Iteration procedures at the material level from the FP and NR methods are monitored at two stress levels: 30 and 40 MPa, as shown in Fig. 3. It is seen that the NR method results in less iteration numbers than the ones using the FP method, while at the structural level both methods require two iterations. The total CPU time required for this case using the NR and FP methods are 3.68 and 3.75 s, respectively. With significantly smaller number of iterations, the CPU time of the NR method is Table 1 Elastic properties and Prony series coefficients for the FM73 adhesive polymer n
kn (s1)
Dn · 106 (MPa)1
1 2 3 4 5 6
1 101 102 103 104 105
21.0 21.6 11.8 15.9 21.6 20.0
E = 2710 MPa m = 0.35
Stress (MPa)
The efficiency and accuracy of the recursive–iterative algorithm with the fixed point (FP) and Newton–Raphson (NR) iterative schemes are examined using experimental data of FM73 polymers. Peretz and Weitsman [19], Peretz and Weitsman [20] performed 15 min creep followed by 15 min recovery tests on FM73 adhesive under various uniaxial stresses: 10–30 MPa and temperatures: 30–60 C. The calibrated time-dependent parameters (Prony series) and elastic properties for FM73 adhesive is given in Table 1. From the experimental results, Peretz and Weitsman [19], Peretz and Weitsman [20] established that the nonlinear stress and temperature dependent material parameters can be fitted using the following functions: r T T0 g0 ¼ 1 þ 0:15 1 þ 0:91 50 T0
r 2:4 T T0 g1 ¼ 1 þ 1:435 exp 8:5 50 T0
ð27Þ r 2 T T0 g2 ¼ 1 þ 0:75 exp 12:12 50 T0
r T T0 a ¼ exp 1:75 40 50 T0
Loading: C r = 1 o sec k = 1 MPa sec
T = rt + T0; σ = kt;
Temperature (ºC)
FM 73 adhesive
3.1. Responses of FM73 adhesive polymers
Response with NR method Response with FP method
Strain (%) Fig. 2. Stress–strain responses on FM73 adhesive under stress and temperature ramp loading using the recursive algorithm with the Newton– Raphson (NR) and fixed point (FP) iterative methods.
FP method, σ=30 MPa FP method, σ=40 MPa NR method, σ=30 MPa NR method, σ=40 MPa
log
[TRol] Convergence line
Iteration number Fig. 3. Residual at the material level from ramp loading using the NR and FP methods.
reduced by only 2%. This is due to the extra computational effort, which requires calculating the residual strain tensor, the derivative of the residual tensor, and its inverse, which is needed when using the NR method. Convergence behaviors of the NR and FP methods are also presented for creep loading under constant stress (30 MPa) and temperature (60 C), shown in Fig. 4. The responses under these loads and temperatures show high nonlinear behaviors. Fig. 5 illustrates the iteration processes at the material levels during the initial and final creep tests. The CPU time using the NR and FP methods are 4.64 and 5.3 s, respectively. The NR method gives more efficient solutions provided that the residual function is continuous; its derivative exists; and is nonsingular within the loading range. The time-dependent responses from the recursive–iterative algorithm are then compared with available experimental data. Figs. 6 and 7 show creep and stress–strain responses under stress and temperature loading histories illustrated in the figures. Both the NR and the FP iterative methods are used and compared with the experimental
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
583
FM 73 adhesive Response with NR method Response with FP method
Response with NR method Response with FP method
Stress (MPa)
Strain (%)
FM 73 adhesive Creep σ=30 MPa and T = 60oC
Experimental data (Peretz and Weitsman, 1983) Input:
σ( MPa) 20
60 3600
Time (sec)
time
30 3600
time
Strain (%)
Fig. 4. Creep strain responses on FM73 adhesive using the recursive algorithm with the NR and FP iteration methods.
FP method, t=10 −6 sec FP method, t = 900 sec NR method, t =10 −6 sec NR method, t = 900 sec log
Τ(C)
[TRol]
Fig. 7. Creep strain responses under linear stress and temperature ramp.
method. The temperature recovery tests were performed by reducing the testing temperature to the reference temperature at the constant load. The stress recovery tests were performed by removing the applied stresses at the fix temperature. (See Figs. 8 and 9.) 3.2. Responses of epoxy Hercules 3502 polymers
Convergence line
Iteration number Fig. 5. Residual at the material level from creep loading using the NR and FP methods.
The capability of the recursive–iterative algorithm in predicting thermo-mechanical creep behaviors is also verified for responses under nonlinear ramp loadings (non-constant loading rates). For this purpose, thermo-mechanical experimental tests on the epoxy Hercules 3502 of Harper and Weitsman [13] are used. The time and temperature dependent material properties were characterized using the creep recovery tests under a constant stress 10 MPa and several temperatures ranging from 50 to 130 C. Table 2 gives the elastic and time-dependent material properties
Experimental data (Peretz and Weitsman, 1983)
FM 73 adhesive σo σult = 0.2 = 10 MPa
Input: σ (MPa)
Strain (%)
Total strain (%)
FM 73 adhesive Response with NR method Response with FP method
T (C) 60
10
30 3600
time
3600
time
Experimental data (Peretz and Weitsman, 1983) Numerical model Recovery from 60oC to 30oC Recovery from 50oC to 30oC Recovery from 40oC to 30oC
Time (sec) Fig. 6. Creep strain responses under linear temperature ramp.
Time (sec)
data. Good predictions are shown. The rest of this study will use the recursive algorithm with the NR iterative
Fig. 8. Temperature recovery responses of FM 73 adhesive at constant stress. (Numerical prediction using the NR method; experimental data of [20].)
584
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
σ (MPa)
FM 73 adhesive σ= 3.7; 10; 15; 18; 22; 26; 31 MPa Τ= 30 º C
127 600
Strain (%)
Experimental data Numerical prediction
Strain (%)
Τ(C)
10 time
1200
27
Hercules 3502
600
1200
time
Experimental data (Harper and Weitsman, 1985) Numerical model (NR method) Time (sec)
Time (sec)
Fig. 10. Strain responses for the Hercules epoxy 3502.
Fig. 9. Stress recovery responses of FM 73 adhesive at constant temperature. (Numerical prediction using the NR method; experimental data of [19].)
Experimental data (Harper and Weitsman,1985) Numerical model (NR method)
n
kn (s1)
Dn · 106 (MPa)1
1 2 3 4 5 6 7 8 9
1 101 102 103 104 105 106 107 108
21 4.7 2.5 9.9 26.6 20.1 104.2 496.5 562.2
Strain (%)
Table 2 Elastic properties and Prony series coefficients for the Hercules 3502 polymer
Hercules 3502 σ (MPa)
T(C) 127
10
2000
E = 4292 MPa m = 0.35
time
27 2000
time
Time (sec) Fig. 11. Strain responses for the Hercules epoxy 3502.
of the Hercules 3502. The temperature dependent parameters are: "
" 0:787 # 2:48 # T T0 T T0 g2 ¼ exp 4:5 g0 ¼ exp 0:675 T0 T0 " " 1:11 # 0:985 # T T0 T T0 g1 ¼ exp 5:79 aT ¼ exp 13:7 T0 T0
ð28Þ
where T0 is the reference temperature, which is 303 K; the temperature T is in Kelvin. Figs. 10 and 11 illustrate creep strain responses for the epoxy Hercules 3502 under two different combinations of stress and temperature histories. The recursive numerical algorithm with the NR iterative method is used. Good predictions are shown for all loading histories. 3.3. Numerical assessment on the performance of TCM This section examines the responses of the TSM and the TCM under various thermo-mechanical loading histories. For this purpose, the FM73 material reported by Peretz
and Weitsman [20] is used. The first study examines the behaviors of the TSM and the TCM under constant temperature rates (slow and fast loading). The responses from the viscoelastic constitutive models having constant (temperature independent) material properties are also presented for comparison. Fig. 12 shows the lower loading rate leads to higher accumulation of the strains, while higher loading rates show smaller strain buildup with increase of the temperature. For the TSM, the difference in strains between the two temperature rates is constant, while for the TCM, the strain differences increase with increasing temperatures. This is due to the additional temperature dependent parameters (g0, g1, and g2). The responses of the TSM and temperature independent materials show nearly equal strain differences under the two temperature rates. These differences are mainly due to the time-dependent material properties. The deviation of the TSM response from the time-independent material is caused by the shift factor aT. The second study presents the material responses under multiple temperature histories, illustrated in Fig. 13.
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588 350 Thermorheologically Complex Material
Τ(ºC)
Thermorheologically Simple Material Temperature independent materials
340
60
Temperature (ºK)
30 100 6000
330
σ(MPa)
(**) (**) (*)
(*)
(**)
t
(*)
15
320
t
6000
310
300 0.005
0.006
0.007
0.008
(*)
0.005 K/sec
(**)
0.3 K/sec
0.009
0.01
Strains Fig. 12. Strain responses under two constant temperature rates and a constant stress.
0.9 Thermorheologically complex material
Strains (%)
Thermorheologically simple material Temperature independent materials
0.6 Τ(C)
σ(MPa)
585
rithms of time-dependent constitutive material models that can be integrated with the FE framework. FE crack propagation analyses of bonded joints having linear viscoelastic adhesives have also been performed by several researchers [16,32,6,7,23]. Most studies are done with adhesive polymers of the TSM typed. In this study, time-dependent and crack propagation analyses are performed on adhesive bonded joints, in which the adhesives follow the TCM. Nonlinear time-dependent responses of the FM73M adhesive single lap joints, tested by Jurf and Vinson [15] at different temperatures and moistures, are predicted first. The geometry of the lap joint is given in Fig. 14b. The FE model for the lap joint is generated using 3D continuum elements (C3D8). Five layers of these elements are generated across the bond thickness. The recursive–iterative algorithm is implemented at the material points within the elements. Short-term creep tests at different temperatures ranging: 22–106 C and at dry, 63% and 95% moisture contents are measured for the FM73M adhesive. The time temperature superposition technique was employed to create a long-term response at the reference temperature for each humidity condition. The calibrated temperature dependent parameters from the shifting processes are then used to predict viscoelastic responses at other non-reference temperatures. Fig. 15a and b show the shear creep compliances from the experimental data and FE analyses at several temperatures and moisture contents. The shear
60 10 30
t
1000 1800 2400
0.3 0
500
1000
2400
1500
t
P 2000
2500
Aluminum
Time (sec) Fig. 13. Strain responses under multiple temperature histories and a constant stress.
Fig. 13 shows strain responses from the TSM and the TCM. Higher strain responses have been observed for the TCM due to the presence of all nonlinear temperature dependent parameters. For the TSM, the temperature dependent property is only carried through a time-shift factor; therefore the transient temperature recovery is not prominent. It is also seen that the initial strains for the TSM and the temperature independent materials are equal. The presence of the time-shift factor in the TSM results in a constant strain difference between the TSM and the temperature independent materials.
Adhesive
P
A
A
Aluminum
4. Nonlinear time-dependent and failure analyses of adhesive bonded joints Adhesive
FE structural analyses have been employed to predict time-dependent behaviors of adhesive bonded lap joints, e.g. Roy and Reddy [25], Yadagirz et al. [34], and Carpenter [3]. The analyses require formulating numerical algo-
SECTION A-A Fig. 14. Geometry models for Mode I and Mode II time-dependent crack propagation analyses on adhesive bonded joints. (a) Double cantilever beam. (b) Lap joint (not to scale, unit in mm).
586
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
crack initiation prediction should also incorporate plastic deformation and thermal stresses due to the different coefficient of thermal expansions in the adherend and adhesive.
FE results Experimental results (Jurf & Vinson, 1985) 79 ºC
4.1. Mode I (opening) time-dependent crack propagation
0.008
71 ºC
0.004
65 ºC
52 ºC 20 ºC
0
0
200
400
600
800
1000
log time (sec)
Shear creep compliance (1/MPA)
0.04 FE results Experimental results (Jurf & Vinson, 1985)
0.03
f ¼
72 ºC
0.02
0.01
68 ºC 65 º C
0
51 º C
0
200
400
600
800
Convergence studies are first performed to determine suitable number of elements used in the crack propagation analyses. Finer meshes are required for the adhesive layers. Total of 4 and 64 elements are generated along the adhesive’s thickness and bonded length, respectively. Cohesive elements are applied along the predetermined failure path, which is located at the mid-height of the adhesive layers. The last four nodes near the fix end are tied, in which the crack will no longer allow propagating. Thus, complete debonding is determined once the current crack tip reaches the tied node. The crack opening displacement (COD) failure criterion is used to evaluate the crack initiation and propagation, which is given as:
1000
log time (sec) Fig. 15. Shear creep compliances calculated from the FE analyses and the experimental data of the bonded lap joint at different temperatures and at humidity; (a) 63% and (b) 95%.
compliance in the FE analysis is calculated by dividing the shear strains by the constant shear stress in the mid-layer of the adhesive. Good agreement is observed in all cases. The second study deals with crack propagation simulations of bonded joints under the opening and the shearing modes. The goal is to predict a crack initiation and a complete failure time. The geometries of the mode I and II are taken from the studies by Allen and Searcy [1] and Jurf and Vinson [15], respectively, which are shown in Fig. 14. FE models for the mode I (opening) and the mode II (shearing) bonded joint failures are generated using plane strain elements (CPE4). For this purpose, a plane strain condition is imposed to the 3D recursive–iterative algorithm. The nonlinear viscoelastic constitutive model is incorporated for the adhesive layers made of FM73M polymer. The aluminum adherends in both lap joints are modeled as a linear elastic material and the temperature effects in the adherends’ elastic properties are ignored. It is noted that in order to simulate more accurate time-dependent crack propagation analyses under thermo-mechanical loadings, detailed constitutive model of the adherends under elevated temperature should be incorporated. In addition, a more accurate
dtn dn;cr
ð29Þ
where dtn is the current normal opening displacement across the cohesive surface, dn,cr is the critical opening displacement. The critical displacement dn,cr is first determined using failure analyses of a quasi-static problem with a critical stress failure criterion: ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s 2 2 r s f ¼ ð30Þ þ rcr scr where rcr and scr are the normal and shear critical stresses, respectively. Thus, the critical displacement dn,cr is the opening displacement when the failure criterion in Eq. 250
P ult.
200
150
Load (N)
Shear creep compliance (1/MPA)
0.012
100
50
0 0
0.04
0.08
0.12
Displacement (mm) Fig. 16. Load–displacement response for the double cantilever beam (mode I) with a displacement rate of 0.2 mm/min.
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
(30) reaches one. The load-displacement response of the quasi-static double cantilever beam (DCB) specimen under a displacement control with a rate of 0.2 mm/min at 295 K is shown in Fig. 16. The failure criterion in Eq. (29) is monitored at a distance of 3 mm ahead of the current crack tip. When the failure criterion reaches one, crack is initiated and the fracture energy is released to create two traction free surfaces. The energy released rate needs to be determined from the experimental crack propagation. The load, at which crack initiates, is around 200 N. This load corresponds to a 5 MPa normal critical stress. The critical opening displacement dn,cr is found to be equal to 0.002 mm. These parameters correspond to the opening failure within the adhesive materials. Next, the time-dependent crack propagation analysis of the double cantilever is performed with the COD failure criterion. Two constant loads: 140 N and 160 N, which correspond to 70% and 80% failure load,
587
Table 3 Crack initiation and propagation time for mode I failure analyses Temperature (K)
Crack initiation time (s)
Debonding time (s)
140 N
160 N
140 N
160 N
295 305 310 315 320 323
– – – 3500 1050 262
2097 52.4 1.64 4e–4 3e–4 2e–4
– – – – – –
– – 323 10.9 11.2 0.5
are applied under several isothermal temperatures for one hour. Fig. 17a and b present the cumulative crack length for the two loads at several temperatures. As expected,
10
12
Cumulative Crack Length (mm)
8 T=295 ºK
Cumulative Crack Length (mm)
T=310 ºK
8
T=315 ºK T=320 ºK T=323 ºK
4
6
4 T=295 ºK
2
T=311 ºK T=323 ºK
0 0
0
6
12
18
24
Time (sec) -4
0
600
1200
1800
2400
3000
10
3600
Time (sec) 24
Cumulative Crack Length (mm)
Cumulative Crack Length (mm)
8 20
T=295 ºK T=305 ºK
16
T=310 ºK T=315 ºK T=320 ºK
12 8 4
6
4 T=295 ºK T=311 ºK
2 0 -4 -600
T=323 ºK
0
600
1200
1800
2400
3000
3600
Time (sec) Fig. 17. Time-dependent cumulative crack length under load levels: (a) 140 N and (b) 160 N.
0 0
3
6
9
12
15
Time (sec) Fig. 18. Time-dependent cumulative crack length under displacement rate: (a) 0.1 mm/min and (b) 0.2 mm/min.
588
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
Table 4 Crack initiation and propagation time for mode II failure analyses
ble of predicting thermo-viscoelastic responses under general stress–temperature histories.
Temperature (K)
Crack initiation time (s)
Debonding time (s)
0.1 mm/ min
0.2 mm/ min
0.1 mm/ min
0.2 mm/ min
Acknowledgement
295 311 323
6.16 6.70 7.41
3.07 3.34 3.68
15.2 19.1 23.0
7.37 9.31 12.1
This research is sponsored by the National Science Foundation under Grant No. 0546528. References
higher load and temperatures accelerate the crack initiation and the complete debonding time. Table 3 presents crack initiation and debonding time at 70% and 80% failure load. 4.2. Mode II (shearing) time-dependent crack propagation A similar procedure as in the mode I is adopted for the mode II failure analysis. In this case, 40 meshes are generated along the bonded length for the adhesive layers and total 5 elements are used along the adhesive’s thickness. Cohesive elements are applied along the predetermined failure paths, which are located at the top and bottom interfaces between the adhesive layers and the aluminum adherend. A stress based failure criterion, Eq. (30), is used. This failure criterion corresponds to the separation between adhesive and aluminum. Fig. 18a and b show the crack propagation analysis for the mode II based on the critical stress. Analyses are run at displacement rates of 0.1 mm/min and 0.2 mm/min and at different temperatures. Softer material reaches the critical stress level slower than the stiffer material does, which cause slower crack propagation. The crack propagates equally in the anti-symmetric failure paths at the top and bottom interfaces. Table 4 shows the crack initiation and debonding time at different displacement rates and temperatures. 5. Conclusions A numerical algorithm is developed for the nonlinear viscoelastic material models of the TCM. The effects of stress and temperature are incorporated in the initial (elastic), long-term (equilibrium), and time-dependent (transient) compliances. The algorithm is compatible with the displacement based FE model for small deformations and uncoupled thermo-mechanical analyses. The recursive algorithm along with iterative schemes is performed at the material level, which enhances convergence and accuracy both at the material and structural levels. The numerical algorithm is verified with the available experimental data under general stress–temperature loading histories. It is shown that the proposed time-integration algorithm of a nonlinear viscoelastic model of the TCM typed is capa-
[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37]
David H. Allen, Chad R. Searcy, Int. J. Fract. 107 (2001) 159–176. Z.P. Bazant, Int. J. Numer. Meth. Eng. 4 (1) (1972) 11–15. W.C. Carpenter, Comput. Struct. 36 (6) (1990) 1141–1152. J.M. Caruthers, R.E. Cohen, Rheol. Acta 19 (1980) 606–613. J.M. Caruthers, D.B. Adolf, R.S. Chambers, P. Shrikhande, Polymer 45 (2004) 4577–4597. B. Chen, D.A. Dillard, Int. J. Solids Struct. 38 (2001) 6907–6924. F. Dubois, C. Chazal, C. Petit, Mech. Time-Depend. Mater. 2 (1999) 269–286. W.W. Feng, Int. J. Nonlinear Mech. 27 (4) (1992) 675–678. R. Haj-Ali, A. Muliana, Int. J. Numer. Meth. Eng. 59 (1) (2004) 25– 45. M. Henriksen, Comput. Struct. 18 (1) (1984) 133–139. H. Hilton, S. Yi, Compos. Eng. 3 (1993) 123–135. G.A. Holzapfel, G. Reiter, Int. J. Eng. Sci. 33 (1995) 1037–1058. B.D. Harper, Y. Weitsman, J. Rheol. 29 (1985) 49–66. A.R. Johnson, T.K. Chen, Comput. Meth. Appl. Mech. Eng. 194 (2005) 313–325. R.A. Jurf, J.R. Vinson, J. Mater. Sci. 20 (1985) 2979–2989. P. Krishnaswamy, M.E. Tuttle, A.F. Emery, J. Ahmad, Int. J. Numer. Meth. Eng. 30 (1990) 371–387. J. Lai, A. Bakker, Comput. Mech. 18 (1996) 182–191. A.H. Muliana, J.S. Kim, Int. J. Solids Struct., in press. D. Peretz, Y. Weitsman, J. Rheol. 26 (3) (1982) 245–261. D. Peretz, Y. Weitsman, J. Rheol. 27 (2) (1983) 97–114. H. Poon, F. Ahmad, Comput. Mech. 21 (1998) 236–242. H. Poon, F. Ahmad, Int. J. Numer. Meth. Eng. 46 (1999) 2027–2041. P. Rahulkumar, A. Jagota, S.J. Bennison, S. Saigal, Int. J. Solids Struct. 37 (2000) 1873–1897. S. Reese, P. Wriggers, ZAMM Appl. Math. Mech. 78 (1998) S157. S. Roy, J.N. Reddy, Comput. Struct. 29 (6) (1988) 1011–1031. Y. Sadkin, J. Aboudi, Compos. Sci. Technol. 36 (1989) 351–365. R.A. Schapery, Polym. Eng. Sci. 9 (1969) 295–310. R.A. Schapery, Compos. Mater. 5 (1974) 85–168. J.C. Simo, T.J.R. Hughes, Comput. Inelast., Springer-Verlag, New York, 1998. R.L. Taylor, K.S. Pister, G.L. Goudreau, Int. J. Numer. Meth. Eng. 2 (1970) 45–59. N.W. Tschoegl, W.G. Knauss, I. Emri, Mech. Time-Dept. Mater. 6 (2002) 53–99. M.K. Warby, J.R. Walton, J.R. Whiteman, Comput. Meth. Appl. Mech. Eng. 97 (1992) 375–397. A.S. Wineman, K.R. Rajagopal, Mechanical Response of Polymers, Cambridge University Press, 2000. S. Yadagirz, C. Reddy, T. Reddy, Comput. Struct. 27 (4) (1987) 445W. S. Yi, H.H. Hilton, M.F. Ahmad, J. Appl. Mech. 63 (1996) 218–224. O.C. Zienkiewicz, M. Watson, I.P. King, Int. J. Mech. Sci. 10 (1968) 807–827. M.A. Zocher, S.E. Groves, D.H. Allen, Int. J. Numer. Meth. Eng. 40 (1997) 2267–2280.
Computational Materials Science 41 (2008) 576–588 www.elsevier.com/locate/commatsci
A time-integration algorithm for thermo-rheologically complex polymers Anastasia Muliana *, Kamran A. Khan Department of Mechanical Engineering, Texas A&M University, College Station, TX 77843-3123, United States Received 30 November 2006; received in revised form 28 May 2007; accepted 30 May 2007 Available online 20 July 2007
Abstract This study presents nonlinear thermo-viscoelastic analyses of isotropic polymers that belong to a class of thermo-rheologically complex materials (TCM). The nonlinear viscoelastic constitutive model is expressed with an integral form of a creep function, whose initial, long-term, and transient properties change with stresses and temperatures. A combined recursive–iterative method is formulated to solve the viscoelastic integral equation. The recurrence formula allows bypassing the need to store entire strain histories. Two types of iterative procedures, fixed point (FP) and Newton–Raphson (NR), are examined within the recursive algorithm. Furthermore, a consistent tangent stiffness matrix is formulated to accelerate convergence and avoid divergence. The algorithm is compatible with displacement based finite element (FE) structural analyses for small deformation gradients and uncoupled thermo-mechanical problems. Verification of the numerical algorithm is performed using the thermo-viscoelastic experimental data available in the literature. The influences of temperature dependent material properties on the material’s viscoelastic responses are discussed for general thermo-mechanical loadings. The capability of the algorithm in predicting multi-axial viscoelastic responses is then verified using creep data of adhesive bonded joints at various temperatures. Finally, numerical simulations of time-dependent crack propagations in adhesive bonded joints are presented for the opening and shearing modes. 2007 Elsevier B.V. All rights reserved. Keywords: Viscoelastic; Thermo-mechanical; Polymers; Finite element; Recursive; Iterative
1. Introduction Polymers exhibit linear and nonlinear viscoelastic behaviors under combined mechanical and environmental loadings. The nonlinear viscoelastic responses, which often occur at high load levels and elevated temperatures, are indicated by non-constant material properties. Based on their rheological responses at different temperatures, polymers can be classified as thermo-rheologically simple materials (TSM) or thermo-rheologically complex materials (TCM). If temperature affects mainly the time-dependent material properties and the effect of temperature can be incorporated through a time-scale shift factor only, materials are categorized as TSM. On the other hand, for *
Corresponding author. Tel.: +1 979 458 3579; fax: +1 979 845 3081. E-mail address: [email protected] (A. Muliana).
0927-0256/$ - see front matter 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2007.05.021
the TCM, temperature influences the material’s initial (elastic), long-term (equilibrium), and time-dependent (transient) properties. Further discussion regarding the TSM and the TCM can be found in Harper and Weitsman [13], Caruthers and Cohen [4], Caruthers et al. [5], Wineman and Rajagopal [33], Tschoegl et al. [31]. Numerical algorithms that are compatible with finite element (FE) analyses have been developed mainly for analyzing viscoelastic constitutive models of the TSM. Zienkiewicz et al. [36] and Bazant [2] presented incremental algorithm for simulating viscoelastic behaviors derived directly from a governing differential equation of a spring-dashpot mechanical model. The strain rate was modeled as functions of the current stresses and the accumulated strains. Other types of constitutive viscoelastic models, i.e. hereditary integral forms, have been formulated for predicting time-dependent material responses
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
under more general loading histories. A recurrence numerical algorithm was first used by Taylor et al. [30] for solving a linear viscoelastic integral equation. The recursive approach allows formulating current stress tensor as a function of current strain increments and history variables stored in the previous time step. This approach bypasses the need to store entire history variables, which minimizes the storage required in performing the integration. Similar recursive approach was used by Feng [8] for linear viscoelastic materials described by the Voltera integral equation. The recursive method has also been used for solving nonlinear viscoelastic constitutive equations of the TSM typed. Henriksen [10], Roy and Reddy [25], and Lai and Bakker [17] developed recursive integration algorithms for solving the Schapery’s [27] nonlinear constitutive model of isotropic materials under combined mechanical and environmental loadings. A Prony exponential series was used in the transient compliance. The computer storage requirement for using this method was proportional to the number of terms in the Prony series. The effects of temperature, moisture, and aging were carried through the time-shift factors only. Haj-Ali and Muliana [9] formulated a recursive– iterative algorithm of the Schapery stress-dependent viscoelastic constitutive models. An iterative scheme based on the Newton–Raphson method was added for correcting the trial linearized stresses and minimizing the residual strains. Simo and Hughes [29] presented several numerical algorithms of inelastic and isotropic materials with internal state variables including viscoelasticity. The hereditary integral models have been extended for anisotropic viscoelastic materials. Poisson’s effect is often assumed to be constant. Zocher et al. [37], Poon and Ahmad [21], Poon and Ahmad [22], Hilton and Yi [11], and Yi et al. [35] developed a FE scheme for linear and nonlinear viscoelastic materials of anisotropic media using integration point update algorithm. To enhance convergence rate, consistent tangent stiffness matrix was formulated. Idealizing thermo-viscoelastic responses of polymers as TSM is sufficient when the materials are subjected to moderate temperature changes and temperature does not vary with time or when the materials are under constant temperature rates and constant mechanical loads. There are many applications when the deformed/stressed bodies are subjected to transient temperatures or when thermo-mechanical coupling due to internal energy dissipation is prominent. Holzapfel and Reiter [12], Reese and Wriggers [24], Johnson and Chen [14] presented time-integration procedures for coupled thermo-viscoelastic behaviors for small and large deformation problems, in which the materials deviate from the TSM. Constitutive equations were derived from the Maxwell mechanical analog model with constant material properties. Reese and Wriggers [24] and Johnson and Chen [14] have pointed out that the constitutive equations should incorporate material degradation due to the temperature changes to accurately predict the thermoviscoelastic behaviors especially when the thermo-viscoelastic coupling is not negligible, however, constitutive
577
models for thermo-viscoelastic responses of the TCM are still lacking. Major studies on the thermo-viscoelastic behaviors of polymers have been focused on a class of TSM. There are many engineering applications that involve transient temperatures and a thermo-mechanical coupling, under which the viscoelastic responses of a materials deviate from the TSM. Thus, there is a need to develop a constitutive equation for the TCM that can be integrated with general structural analysis tools. This study presents a time-integration algorithm for predicting thermo-viscoelastic behaviors of isotropic polymers that belong to a class of TCM. This algorithm is designed to be compatible with the FE framework to perform large scale structural analyses under general thermo-mechanical loadings. For this purpose, the effects of stress and temperature should be incorporated in the initial, long-term, and time-dependent material properties. The theoretical basis of the studied viscoelastic constitutive equation is given in Harper and Weitsman [13]. In Section 2, a recursive–iterative algorithm is formulated. Two types of iterative procedures, fixed point (FP) and Newton–Raphson (NR), are examined within the recursive scheme. Furthermore, a consistent tangent stiffness matrix is formulated at the material level to accelerate convergence and avoid divergence at the structural level. Verifications of the numerical algorithm with thermo-mechanical viscoelastic experimental data available in the literature are presented in Section 3. Convergence behaviors using the FP and NR methods are examined for several thermo-mechanical loading histories. In addition, performances of the TCM under general thermo-viscoelastic behaviors are discussed based on the numerical results. In Section 4, examples of the integrated material model with the FE analyses are given on time-dependent and crack propagation analyses of polymeric adhesive bonded joints. Creep data on the adhesive bonded lap joints at different temperatures and humidity are used for comparisons. Numerical predictions of a crack initiation and a complete debonding time for mode I (opening) and II (shearing) failures are also presented. 2. Formulation of the time-integration algorithm Schapery [28] and Sadkin and Aboudi [26] have derived constitutive material models for the TCM that involves two time-scales for incorporating temperature and stress variations. The use of two time-scales in the constitutive material models makes development of numerical algorithm very challenging and is often difficult to integrate these constitutive material models in general FE structural analyses. In this study, the time–stress–temperature effects on the material behaviors are linked to one time-scale without any restriction for modeling temperature or stress variations with time. The Schapery [27] nonlinear single integral constitutive model is modified to include the effects of stress and temperature on the initial, long-term, and time-dependent material properties for non-aging materials, which is given as:
578
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
Z t t s et eðtÞ ¼ g0 ðrt ; T t ÞD0 rt þ g1 ðrt ; T t Þ DDðw w Þ 0 Z t s d½g0 ðrs ; T s Þrs s dDT ds þ ds aðT Þ ds ds 0
ð1Þ
Here D0 and DD are the uniaxial instantaneous elastic and transient compliances. The uniaxial transient compliance is expressed as: t
DDw ¼
N X
Dn ð1 exp½kn wt Þ
ð2Þ
n¼1
w is the reduced-time (effective time) given by: Z t Z s dn dn t s w wðsÞ ¼ w wðtÞ ¼ n n n n 0 a r ;T 0 aðr ; T Þ
where D0 ¼
1 E0
Dð1Þ ¼ g0 ðrt ; T t ÞD0 þ g1 ðrt ; T t Þg2 ðrt ; T t Þ
ð4Þ N X
Dn
n¼1
1 where Dn ¼ En
ð5Þ
The uniaxial viscoelastic behavior in Eq. (1) is generalized for multi-axial (3D) constitutive relations of isotropic materials by separating the deviatoric and volumetric strain–stress relations as: 1 etij ¼ etij þ etkk dij þ aðT t T 0 Þdij 3 1 t eij ¼ g0 ð rt ; T t ÞJ 0 S tij 2 Z t o½g2 ð rs ; T s ÞS sij 1 wt ws Þ t t ð þ g1 ð r ;T Þ DJ ds 2 os 0
ð6Þ
The parameters J0 and B0 are the instantaneous elastic shear and bulk compliances, respectively. The terms DJ and DB are the time-dependent shear and bulk compliances, respectively. The nonlinear parameters of the multi-axial behaviors are modeled as a function of current t and temperature Tt. The matrix Poisson’s effective stress r ratio, m, is assumed to be time-independent, which allows expressing the shear and bulk compliances as: J 0 ¼ 2ð1 þ mÞD0
ð3Þ
The constitutive model in Eq. (1) can be visualized as series arrangement of N Kelvin–Voigt elements and an elastic spring. Each spring has stress and temperature dependent material constants. The parameters g0 and g1 measure changes in the elastic and the transient compliances, respectively, with stresses and temperatures. The parameter g2 accounts for the loading rate effect on the time-dependent response. The parameter aðrt ; T t Þ is a time-shift (interchange) factor measured with respect to the reference stress and temperature. If the thermal effects are carried through the shift factor a(Tt) only, the thermo-rheologically simple material (TSM) is exhibited. The parameters T and T0 are the current and the reference temperatures, respectively. The superscript denotes the time-dependent variable. In general, the coefficient of thermal expansion a depends on the temperature; in this study a is assumed to be constant. The initial and the long-term nonlinear creep compliances can be determined from a step stress history by setting time t = 0 and 1, respectively, which are expressed by: Dð0Þ ¼ g0 ðrt ; T t ÞD0
1 rt ; T t ÞB0 rtkk etkk ¼ g0 ð 3 Z t 1 t s o½g 2 ð rs ; T s Þrskk þ g1 ð rt ; T t Þ DBðw w Þ ds 3 os 0
DJ
wt
B0 ¼ 3ð1 2mÞD0 t
¼ 2ð1 þ mÞDDw
t
t
DBw ¼ 3ð1 2mÞDDw
ð7Þ
For simplicity, the nonlinear material parameters will be written as: gb ð rt ; T t Þ ¼ gtb að rt ; T t Þ ¼ at
b ¼ 0; 1; 2 rffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 t t t ¼ S S r 2 ij ij
ð8Þ
2.1. Recursive formulation A recursive method is used to solve the integrals in Eq. (6). The algorithm is designed to be compatible with general displacement based FE structural analyses for uncoupled thermo-mechanical loadings. The numerical algorithm of the TCM will be called at every material (Gaussian) integration point within elements at each structural iteration to achieve structural and material convergences simultaneously. Thus, an efficient and accurate numerical algorithm for solving the constitutive material model becomes necessary. Using the transient compliance in Eq. (2) and the relations in Eq. (7), the strains in Eq. (6) can be written as: N N X 1 1 1 X etij ¼ gt0 J 0 S tij þ gt1 S tij J n gt1 J n qtij;n 2 2 2 n¼1 n¼1 Z t o gs2 S sij where qtij;n ¼ ds exp½kn ðwt ws Þ os 0 N N X 1 1 1 X etkk ¼ gt0 B0 rtkk þ gt1 rtkk Bn gt1 Bn qtkk;n 3 3 3 n¼1 n¼1 s s Z t t s o g 2 rkk t where qkk;n ¼ ds exp½kn ðw w Þ os 0
ð9Þ
ð10Þ
A recursive integration method, which was derived by Taylor et al. [30] for a linear viscoelastic integral model and extended by Haj-Ali and Muliana [9] for a nonlinear stress-dependent viscoelastic constitutive model, is used to solve the integral parts in Eqs. (9) and (10). These are expressed as:
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
qtij;n
¼
Z
tDt s
t
exp½kn ðw w Þ
o gs2 S sij
Detij ¼ etij eijtDt ¼ Jt S tij JtDt S tDt ij
ds os Z t o gs2 S sij ds þ exp½kn ðwt ws Þ os tDt t t tDt tDt ¼ exp½kn Dwt qtDt S ij ij;n þ g2 S ij g2 0
qtkk;n
1 exp½kn Dwt kn Dwt
ð11Þ
tDt
t
1 exp½kn Dwt kn Dwt
s
ð12Þ
The first integral includes the limits (0, t Dt), i.e. up to the previous converged time step, which are stored as history variables of the hereditary strain integrals. The limits of the second integral are taken as (t Dt, t), which is the current incremental step. Detailed derivation of the recursive– iterative algorithm for the nonlinear viscoelastic model can tDt be found in Haj-Ali and Muliana [9]. The parameters qij;n tDt and qkk;n , n = 1, . . ., N are the hereditary integral for every term in the Prony series in the form of deviatoric and volumetric strains. These parameters are history state variables stored from the last converged step at time (t Dt). The total number of history variables is 7 · N · NGauss; where 7 is the total number of deviatoric and volumetric components, N is the number of terms in the Prony series, and NGauss is the number of Gaussian integration points. The incremental reduced-time is expressed as Dwt wt wtDt. Substituting the integral Eqs. (11) and (12) into Eqs. (9) and (10) gives the deviatoric and volumetric total strains: etij
" # N N X X 1 t 1 exp ½kn Dwt t t t t t g J 0 þ g1 g2 ¼ J n g1 g 2 Jn S ij 2 0 kn Dwt n¼1 n¼1
N t 1 X tDt 1 exp½kn Dw tDt gt1 J n exp½kn Dwt qtDt S ij t ij;n g 2 2 n¼1 kn Dw
Jt S tij d tij ð13Þ " # N N X X 1 t 1 exp½kn Dwt t rkk etkk ¼ Bn gt1 gt2 Bn g0 B0 þ gt1 gt2 3 kn Dwt n¼1 n¼1
N t 1 X tDt 1 exp½kn Dw tDt gt1 Bn exp ½kn Dwt qtDt g r kk;n 2 kk 3 n¼1 kn Dwt t rt V t ð14Þ B kk
N tDt 1X J n gt1 exp ½kn Dwt g1tDt qij;n 2 n¼1 " ! tDt N 1 tDt X tDt 1 exp kn Dw g2 J n g1 2 kn DwtDt n¼1 # t t 1 exp ½kn Dw g1 S ijtDt kn Dwt
o gs2 rskk ds ¼ exp½kn ðw w Þ os 0 Z t oðgs2 rskk Þ ds þ exp½kn ðwt ws Þ os tDt t t tDt tDt ¼ exp½kn Dwt qtDt rkk kk;n þ g2 rkk g 2 Z
579
kk
The parameters J t and bulk Bt are the current time-dependent deviatoric and volumetric compliances, respectively, while d tij and vtkk are the deviatoric and volumetric strains that carry the history variables. The incremental forms of the deviatoric and volumetric strains are:
t rt B tDt rtDt ¼B Detkk ¼ etkk etDt kk kk kk N t 1X Bn g1 exp ½kn Dwt g1tDt qtDt kk;n 3 n¼1 " ! tDt N 1 tDt X tDt 1 exp kn Dw g2 Bn g 1 3 kn DwtDt n¼1 # t t 1 exp ½kn Dw g1 rtDt kk kn Dwt
ð15Þ
ð16Þ
Eqs. (13)–(16) define complete solutions for the total and incremental strain tensors, respectively. The nonlinear parameters in Eqs. (13)–(16) are expressed as functions of the current temperature and total stresses; however, the total stresses at the current time (t) are not known. Linearized trial stresses are used as a starting point to obtain the current stresses using either Eqs. (13) and (14) or (15) and (16). An iterative scheme is added in order to find the correct stress tensor. The chosen trial values (initial approximation solutions) can significantly determine the convergence solutions. Poor initial approximation values often lead to divergence. Two iterative methods, which are fixed point (FP) and Newton–Raphson’s (NR) are used in this study. The FP method gives a slower but stable convergence rate and requires less computational effort during each iteration step. The NR method can increase convergence process rapidly, however, linearized tangent functions often lead to divergence solutions if during the iteration process the convergence is not monotonic or if the magnitude of the tangent functions is very small (approaching zero). 2.2. Fixed point (FP) method The iteration procedure for solving a variable x that satisfies x = f(x) in the FP method is summarized as: xðkþ1Þ ¼ f ðxðkÞ Þ R ¼ kxðkþ1Þ xðkÞ k Converged solutions:
ð17Þ R!0
where k indicates the iteration counter with an initial approximation solution x(0) and R is the measured error. The current stress tensor at the iteration k from a given mechanical strain is: i 1 1 1h rtij ¼ S tij þ rtkk dij ¼ t etij þ d tij þ t etkk þ vtkk dij ð18Þ 3 J B
580
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
The effective shear J t and bulk Bt compliances, and the nonlinear parameters in the Eq. (18) are functions of the t and temperatures Tt. The curcurrent effective stresses r rent temperature is defined by Tt = TtDt + DTt. The initial approximation (trial) stress tensor is determined using the following approximation of nonlinear parameters: gt;tr b a
t;tr
¼
t;ð0Þ gb
¼a
t;ð0Þ
¼ gb ð r ¼ að r
tDt
tDt
;T
;T
tDt
tDt
Þ b ¼ 0; 1; 2
Þ
ð19Þ
The trial stress tensor is formed based on the given (known) variables and history variables: t;ð0Þ
rt;tr ij ¼ rij
t;ð0Þ
¼ rijtDt þ Drij t;ð0Þ
¼ rijtDt þ DS ij
1 t;ð0Þ þ Drkk dij 3
ð20Þ
where the trial incremental deviatoric and volumetric stresses are given by: " # N X 1 1 t;ð0Þ t;ð0Þ t tDt DS ij ¼ t;ð0Þ Detij þ g1 J n ðexp½kn Dw 1Þqij;n 2 J n¼1 " # N 1 1 t;ð0Þ X t;ð0Þ t t tDt Drkk ¼ t;ð0Þ Dekk þ g1 Bn ðexp½kn Dw 1Þqkk;n 3 B
defined by using either the incremental strains or the total strains. The residual strain tensor is expressed by: t;ðkÞ
Rij
t;ðkÞ
¼ Deij
1 t;ðkÞ þ Dekk dij D^etij 3
ð25Þ
The incremental deviatoric and volumetric strain tensors at the iteration k in the material level are defined using Eqs. (15) and (16). The current incremental mechanical strains t;ðmÞ D^etij D^eij is obtained from the structural level at the mth global iteration. The initial approximation (trial) stress tensor is determined using Eqs. (19)–(21). A Jacobian tensor is formed by taking the derivative of the residual tensor in t Eq. (25) with respect to the incremental stress tensor oRij . Finally, the consistent tangent stiffness matrix C tijkl oDrtkl at the converged state is: t 1 oDrtij oRij t C ijkl ¼ ; kRtij k ! 0 ð26Þ oD^etkl oDrtkl Once convergence has been achieved, the stress, the hereditary integrals for each Prony series, and the consistent tangent stiffness are updated. The recursive–iterative algorithm is summarized in Fig. 1.
n¼1
ð21Þ During the iteration, the component of the residual tensor is expressed in terms of the total stress: t;ðkþ1Þ Rij
¼
t;ðkþ1Þ rij
t;ðkÞ rij
ð22Þ
The nonlinear parameters in Eq. (18) is defined by t;ðkÞ ð t;ðkÞ ; T t Þ and at;ðkÞ ¼ að gb ¼ g b r rt;ðkÞ ; T t Þ. Once convergence is achieved, the consistent tangent stiffness matrix in Eq. (23), the current stress tensor, and the history variables in Eqs. (11) and (12) are updated and sent to the structural level. In the coupled thermo-mechanical analyort ses, it is also necessary to determine oTijt .
ortij 1 1 1 1 t ð23Þ C ijkl ¼ t ¼ t dik djk þ dij 3 Bt J t o^ekl J The variable ^etij is the component of the mechanical strain tensor. 2.3. Newton–Raphson’s (NR) method The iterative procedure using the NR’s method for solving a variable x that satisfies F(x) = 0 is summarized as: xðkþ1Þ ¼ xðkÞ
F ðxðkÞ Þ F 0 ðxðkÞ Þ
dF ðxÞ dx Converged solutions:
ð24Þ
F 0 ðxÞ ¼
kRk ¼ kF k ! 0
This procedure requires that F(x) and F 0 (x) are continuous near the solution x and the magnitude of F 0 (x) should not approach zero or infinite. The residual tensor can be
2.4. Implementation of the recursive–iterative algorithm in the FE framework The constitutive thermo-viscoelastic material model with linearized and corrector schemes (Fig. 1) is implemented in the user material (UMAT) ABAQUS FE (2005). At the FE structural level, an iterative equation solution is also performed. During this structural iteration the recursive–iterative algorithm will be called at each material (Gaussian) point. In every iteration within an incremental time step Dt(m), trial incremental strain tensor t;ðmÞ Deij and temperature DTt,(m) are prescribed, which are the input variables to the subroutine UMAT. The superscript (m) denotes the global iteration counter within the current incremental time step. The output from the UMAT t;ðmÞ is the total stresses rij and the consistent tangent stiffness t;ðmÞ t;ðMÞ C ijkl . The converged C ijkl after M global iteration at the current time t will be used to provide incremental trial strains for the next time step (t + Dt). In the uncoupled thermo-mechanical analyses, the trial incremental temperature is directly linked to the incremental time step. By adding the corrector scheme at the material level, it is not necessary to have a tight incremental time Dt to achieve convergence at the structural level. In fact, only by tightening Dt without adding a corrector scheme at the material level can lead to divergence, which was discussed in the previous studies by Haj-Ali and Muliana [9] and Muliana and Kim [18]. At the material level, the convergence criterion can be defined using the residual vectors in Eqs. (22) or (25), which are expressed in terms of stress and strain vectors, respectively. In this study, the given tolerances allow for the maximum residual stress to be 0.01% of the total stress and the maximum residual strain to be 1 micro-
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
581
Fig. 1. Recursive–iterative algorithm for the nonlinear isotropic viscoelastic material.
strain. The above tolerances are set as the limits since reducing the tolerances below these values does not change the overall responses within our numerical values of interest, while increasing the tolerance beyond these values can lead to divergence in the overall responses. The users should also realize that the stress vector at the local (material) iteration (k + 1) in the NR iterative method is defined as an intersection between the stress vector and the tangent oRt to the residual vector oDrijt at the iteration k. Thus, if at the kl
iteration k the value of
oRtij oDrtkl
is approaching zero or infinite
convergence at the material level will not be achieved, which indicates by either an extremely large stress correction or a zero correction. In this situation, adding convergence criteria such as a maximum number of iteration and a minimum correction rate can be used to terminate the
iteration at the material level and a smaller time increment Dt can be chosen to possibly achieve a convergence.
3. Verification of the thermo-mechanical viscoelastic constitutive model The proposed numerical algorithm of the TCM is verified using the thermo-viscoelastic experimental data available in the literature. Test data on FM73 adhesive polymers reported by Peretz and Weitsman [19], Peretz and Weitsman [20] and epoxy Hercules 3502 of Harper and Weitsman [13] are used. The algorithm is implemented in the 3D continuum elements. The 3D nonlinear viscoelastic material response can also be integrated with other typed of elements e.g: 1D (truss, beam) or 2D (plane stress,
582
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
plane strain) by imposing uniaxial or in-plane stress–strain relations to the 3D constitutive model.
where T0 is the reference temperature, which is 303 K; the stress r and temperature T are in MPa and Kelvin, respectively. Convergence behaviors using the FP and NR methods are presented first. Fig. 2 illustrates the nonlinear elastic responses under a linear stress and temperature ramp loading. Iteration processes are performed both at the material and structural levels. Iteration procedures at the material level from the FP and NR methods are monitored at two stress levels: 30 and 40 MPa, as shown in Fig. 3. It is seen that the NR method results in less iteration numbers than the ones using the FP method, while at the structural level both methods require two iterations. The total CPU time required for this case using the NR and FP methods are 3.68 and 3.75 s, respectively. With significantly smaller number of iterations, the CPU time of the NR method is Table 1 Elastic properties and Prony series coefficients for the FM73 adhesive polymer n
kn (s1)
Dn · 106 (MPa)1
1 2 3 4 5 6
1 101 102 103 104 105
21.0 21.6 11.8 15.9 21.6 20.0
E = 2710 MPa m = 0.35
Stress (MPa)
The efficiency and accuracy of the recursive–iterative algorithm with the fixed point (FP) and Newton–Raphson (NR) iterative schemes are examined using experimental data of FM73 polymers. Peretz and Weitsman [19], Peretz and Weitsman [20] performed 15 min creep followed by 15 min recovery tests on FM73 adhesive under various uniaxial stresses: 10–30 MPa and temperatures: 30–60 C. The calibrated time-dependent parameters (Prony series) and elastic properties for FM73 adhesive is given in Table 1. From the experimental results, Peretz and Weitsman [19], Peretz and Weitsman [20] established that the nonlinear stress and temperature dependent material parameters can be fitted using the following functions: r T T0 g0 ¼ 1 þ 0:15 1 þ 0:91 50 T0
r 2:4 T T0 g1 ¼ 1 þ 1:435 exp 8:5 50 T0
ð27Þ r 2 T T0 g2 ¼ 1 þ 0:75 exp 12:12 50 T0
r T T0 a ¼ exp 1:75 40 50 T0
Loading: C r = 1 o sec k = 1 MPa sec
T = rt + T0; σ = kt;
Temperature (ºC)
FM 73 adhesive
3.1. Responses of FM73 adhesive polymers
Response with NR method Response with FP method
Strain (%) Fig. 2. Stress–strain responses on FM73 adhesive under stress and temperature ramp loading using the recursive algorithm with the Newton– Raphson (NR) and fixed point (FP) iterative methods.
FP method, σ=30 MPa FP method, σ=40 MPa NR method, σ=30 MPa NR method, σ=40 MPa
log
[TRol] Convergence line
Iteration number Fig. 3. Residual at the material level from ramp loading using the NR and FP methods.
reduced by only 2%. This is due to the extra computational effort, which requires calculating the residual strain tensor, the derivative of the residual tensor, and its inverse, which is needed when using the NR method. Convergence behaviors of the NR and FP methods are also presented for creep loading under constant stress (30 MPa) and temperature (60 C), shown in Fig. 4. The responses under these loads and temperatures show high nonlinear behaviors. Fig. 5 illustrates the iteration processes at the material levels during the initial and final creep tests. The CPU time using the NR and FP methods are 4.64 and 5.3 s, respectively. The NR method gives more efficient solutions provided that the residual function is continuous; its derivative exists; and is nonsingular within the loading range. The time-dependent responses from the recursive–iterative algorithm are then compared with available experimental data. Figs. 6 and 7 show creep and stress–strain responses under stress and temperature loading histories illustrated in the figures. Both the NR and the FP iterative methods are used and compared with the experimental
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
583
FM 73 adhesive Response with NR method Response with FP method
Response with NR method Response with FP method
Stress (MPa)
Strain (%)
FM 73 adhesive Creep σ=30 MPa and T = 60oC
Experimental data (Peretz and Weitsman, 1983) Input:
σ( MPa) 20
60 3600
Time (sec)
time
30 3600
time
Strain (%)
Fig. 4. Creep strain responses on FM73 adhesive using the recursive algorithm with the NR and FP iteration methods.
FP method, t=10 −6 sec FP method, t = 900 sec NR method, t =10 −6 sec NR method, t = 900 sec log
Τ(C)
[TRol]
Fig. 7. Creep strain responses under linear stress and temperature ramp.
method. The temperature recovery tests were performed by reducing the testing temperature to the reference temperature at the constant load. The stress recovery tests were performed by removing the applied stresses at the fix temperature. (See Figs. 8 and 9.) 3.2. Responses of epoxy Hercules 3502 polymers
Convergence line
Iteration number Fig. 5. Residual at the material level from creep loading using the NR and FP methods.
The capability of the recursive–iterative algorithm in predicting thermo-mechanical creep behaviors is also verified for responses under nonlinear ramp loadings (non-constant loading rates). For this purpose, thermo-mechanical experimental tests on the epoxy Hercules 3502 of Harper and Weitsman [13] are used. The time and temperature dependent material properties were characterized using the creep recovery tests under a constant stress 10 MPa and several temperatures ranging from 50 to 130 C. Table 2 gives the elastic and time-dependent material properties
Experimental data (Peretz and Weitsman, 1983)
FM 73 adhesive σo σult = 0.2 = 10 MPa
Input: σ (MPa)
Strain (%)
Total strain (%)
FM 73 adhesive Response with NR method Response with FP method
T (C) 60
10
30 3600
time
3600
time
Experimental data (Peretz and Weitsman, 1983) Numerical model Recovery from 60oC to 30oC Recovery from 50oC to 30oC Recovery from 40oC to 30oC
Time (sec) Fig. 6. Creep strain responses under linear temperature ramp.
Time (sec)
data. Good predictions are shown. The rest of this study will use the recursive algorithm with the NR iterative
Fig. 8. Temperature recovery responses of FM 73 adhesive at constant stress. (Numerical prediction using the NR method; experimental data of [20].)
584
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
σ (MPa)
FM 73 adhesive σ= 3.7; 10; 15; 18; 22; 26; 31 MPa Τ= 30 º C
127 600
Strain (%)
Experimental data Numerical prediction
Strain (%)
Τ(C)
10 time
1200
27
Hercules 3502
600
1200
time
Experimental data (Harper and Weitsman, 1985) Numerical model (NR method) Time (sec)
Time (sec)
Fig. 10. Strain responses for the Hercules epoxy 3502.
Fig. 9. Stress recovery responses of FM 73 adhesive at constant temperature. (Numerical prediction using the NR method; experimental data of [19].)
Experimental data (Harper and Weitsman,1985) Numerical model (NR method)
n
kn (s1)
Dn · 106 (MPa)1
1 2 3 4 5 6 7 8 9
1 101 102 103 104 105 106 107 108
21 4.7 2.5 9.9 26.6 20.1 104.2 496.5 562.2
Strain (%)
Table 2 Elastic properties and Prony series coefficients for the Hercules 3502 polymer
Hercules 3502 σ (MPa)
T(C) 127
10
2000
E = 4292 MPa m = 0.35
time
27 2000
time
Time (sec) Fig. 11. Strain responses for the Hercules epoxy 3502.
of the Hercules 3502. The temperature dependent parameters are: "
" 0:787 # 2:48 # T T0 T T0 g2 ¼ exp 4:5 g0 ¼ exp 0:675 T0 T0 " " 1:11 # 0:985 # T T0 T T0 g1 ¼ exp 5:79 aT ¼ exp 13:7 T0 T0
ð28Þ
where T0 is the reference temperature, which is 303 K; the temperature T is in Kelvin. Figs. 10 and 11 illustrate creep strain responses for the epoxy Hercules 3502 under two different combinations of stress and temperature histories. The recursive numerical algorithm with the NR iterative method is used. Good predictions are shown for all loading histories. 3.3. Numerical assessment on the performance of TCM This section examines the responses of the TSM and the TCM under various thermo-mechanical loading histories. For this purpose, the FM73 material reported by Peretz
and Weitsman [20] is used. The first study examines the behaviors of the TSM and the TCM under constant temperature rates (slow and fast loading). The responses from the viscoelastic constitutive models having constant (temperature independent) material properties are also presented for comparison. Fig. 12 shows the lower loading rate leads to higher accumulation of the strains, while higher loading rates show smaller strain buildup with increase of the temperature. For the TSM, the difference in strains between the two temperature rates is constant, while for the TCM, the strain differences increase with increasing temperatures. This is due to the additional temperature dependent parameters (g0, g1, and g2). The responses of the TSM and temperature independent materials show nearly equal strain differences under the two temperature rates. These differences are mainly due to the time-dependent material properties. The deviation of the TSM response from the time-independent material is caused by the shift factor aT. The second study presents the material responses under multiple temperature histories, illustrated in Fig. 13.
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588 350 Thermorheologically Complex Material
Τ(ºC)
Thermorheologically Simple Material Temperature independent materials
340
60
Temperature (ºK)
30 100 6000
330
σ(MPa)
(**) (**) (*)
(*)
(**)
t
(*)
15
320
t
6000
310
300 0.005
0.006
0.007
0.008
(*)
0.005 K/sec
(**)
0.3 K/sec
0.009
0.01
Strains Fig. 12. Strain responses under two constant temperature rates and a constant stress.
0.9 Thermorheologically complex material
Strains (%)
Thermorheologically simple material Temperature independent materials
0.6 Τ(C)
σ(MPa)
585
rithms of time-dependent constitutive material models that can be integrated with the FE framework. FE crack propagation analyses of bonded joints having linear viscoelastic adhesives have also been performed by several researchers [16,32,6,7,23]. Most studies are done with adhesive polymers of the TSM typed. In this study, time-dependent and crack propagation analyses are performed on adhesive bonded joints, in which the adhesives follow the TCM. Nonlinear time-dependent responses of the FM73M adhesive single lap joints, tested by Jurf and Vinson [15] at different temperatures and moistures, are predicted first. The geometry of the lap joint is given in Fig. 14b. The FE model for the lap joint is generated using 3D continuum elements (C3D8). Five layers of these elements are generated across the bond thickness. The recursive–iterative algorithm is implemented at the material points within the elements. Short-term creep tests at different temperatures ranging: 22–106 C and at dry, 63% and 95% moisture contents are measured for the FM73M adhesive. The time temperature superposition technique was employed to create a long-term response at the reference temperature for each humidity condition. The calibrated temperature dependent parameters from the shifting processes are then used to predict viscoelastic responses at other non-reference temperatures. Fig. 15a and b show the shear creep compliances from the experimental data and FE analyses at several temperatures and moisture contents. The shear
60 10 30
t
1000 1800 2400
0.3 0
500
1000
2400
1500
t
P 2000
2500
Aluminum
Time (sec) Fig. 13. Strain responses under multiple temperature histories and a constant stress.
Fig. 13 shows strain responses from the TSM and the TCM. Higher strain responses have been observed for the TCM due to the presence of all nonlinear temperature dependent parameters. For the TSM, the temperature dependent property is only carried through a time-shift factor; therefore the transient temperature recovery is not prominent. It is also seen that the initial strains for the TSM and the temperature independent materials are equal. The presence of the time-shift factor in the TSM results in a constant strain difference between the TSM and the temperature independent materials.
Adhesive
P
A
A
Aluminum
4. Nonlinear time-dependent and failure analyses of adhesive bonded joints Adhesive
FE structural analyses have been employed to predict time-dependent behaviors of adhesive bonded lap joints, e.g. Roy and Reddy [25], Yadagirz et al. [34], and Carpenter [3]. The analyses require formulating numerical algo-
SECTION A-A Fig. 14. Geometry models for Mode I and Mode II time-dependent crack propagation analyses on adhesive bonded joints. (a) Double cantilever beam. (b) Lap joint (not to scale, unit in mm).
586
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
crack initiation prediction should also incorporate plastic deformation and thermal stresses due to the different coefficient of thermal expansions in the adherend and adhesive.
FE results Experimental results (Jurf & Vinson, 1985) 79 ºC
4.1. Mode I (opening) time-dependent crack propagation
0.008
71 ºC
0.004
65 ºC
52 ºC 20 ºC
0
0
200
400
600
800
1000
log time (sec)
Shear creep compliance (1/MPA)
0.04 FE results Experimental results (Jurf & Vinson, 1985)
0.03
f ¼
72 ºC
0.02
0.01
68 ºC 65 º C
0
51 º C
0
200
400
600
800
Convergence studies are first performed to determine suitable number of elements used in the crack propagation analyses. Finer meshes are required for the adhesive layers. Total of 4 and 64 elements are generated along the adhesive’s thickness and bonded length, respectively. Cohesive elements are applied along the predetermined failure path, which is located at the mid-height of the adhesive layers. The last four nodes near the fix end are tied, in which the crack will no longer allow propagating. Thus, complete debonding is determined once the current crack tip reaches the tied node. The crack opening displacement (COD) failure criterion is used to evaluate the crack initiation and propagation, which is given as:
1000
log time (sec) Fig. 15. Shear creep compliances calculated from the FE analyses and the experimental data of the bonded lap joint at different temperatures and at humidity; (a) 63% and (b) 95%.
compliance in the FE analysis is calculated by dividing the shear strains by the constant shear stress in the mid-layer of the adhesive. Good agreement is observed in all cases. The second study deals with crack propagation simulations of bonded joints under the opening and the shearing modes. The goal is to predict a crack initiation and a complete failure time. The geometries of the mode I and II are taken from the studies by Allen and Searcy [1] and Jurf and Vinson [15], respectively, which are shown in Fig. 14. FE models for the mode I (opening) and the mode II (shearing) bonded joint failures are generated using plane strain elements (CPE4). For this purpose, a plane strain condition is imposed to the 3D recursive–iterative algorithm. The nonlinear viscoelastic constitutive model is incorporated for the adhesive layers made of FM73M polymer. The aluminum adherends in both lap joints are modeled as a linear elastic material and the temperature effects in the adherends’ elastic properties are ignored. It is noted that in order to simulate more accurate time-dependent crack propagation analyses under thermo-mechanical loadings, detailed constitutive model of the adherends under elevated temperature should be incorporated. In addition, a more accurate
dtn dn;cr
ð29Þ
where dtn is the current normal opening displacement across the cohesive surface, dn,cr is the critical opening displacement. The critical displacement dn,cr is first determined using failure analyses of a quasi-static problem with a critical stress failure criterion: ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s 2 2 r s f ¼ ð30Þ þ rcr scr where rcr and scr are the normal and shear critical stresses, respectively. Thus, the critical displacement dn,cr is the opening displacement when the failure criterion in Eq. 250
P ult.
200
150
Load (N)
Shear creep compliance (1/MPA)
0.012
100
50
0 0
0.04
0.08
0.12
Displacement (mm) Fig. 16. Load–displacement response for the double cantilever beam (mode I) with a displacement rate of 0.2 mm/min.
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
(30) reaches one. The load-displacement response of the quasi-static double cantilever beam (DCB) specimen under a displacement control with a rate of 0.2 mm/min at 295 K is shown in Fig. 16. The failure criterion in Eq. (29) is monitored at a distance of 3 mm ahead of the current crack tip. When the failure criterion reaches one, crack is initiated and the fracture energy is released to create two traction free surfaces. The energy released rate needs to be determined from the experimental crack propagation. The load, at which crack initiates, is around 200 N. This load corresponds to a 5 MPa normal critical stress. The critical opening displacement dn,cr is found to be equal to 0.002 mm. These parameters correspond to the opening failure within the adhesive materials. Next, the time-dependent crack propagation analysis of the double cantilever is performed with the COD failure criterion. Two constant loads: 140 N and 160 N, which correspond to 70% and 80% failure load,
587
Table 3 Crack initiation and propagation time for mode I failure analyses Temperature (K)
Crack initiation time (s)
Debonding time (s)
140 N
160 N
140 N
160 N
295 305 310 315 320 323
– – – 3500 1050 262
2097 52.4 1.64 4e–4 3e–4 2e–4
– – – – – –
– – 323 10.9 11.2 0.5
are applied under several isothermal temperatures for one hour. Fig. 17a and b present the cumulative crack length for the two loads at several temperatures. As expected,
10
12
Cumulative Crack Length (mm)
8 T=295 ºK
Cumulative Crack Length (mm)
T=310 ºK
8
T=315 ºK T=320 ºK T=323 ºK
4
6
4 T=295 ºK
2
T=311 ºK T=323 ºK
0 0
0
6
12
18
24
Time (sec) -4
0
600
1200
1800
2400
3000
10
3600
Time (sec) 24
Cumulative Crack Length (mm)
Cumulative Crack Length (mm)
8 20
T=295 ºK T=305 ºK
16
T=310 ºK T=315 ºK T=320 ºK
12 8 4
6
4 T=295 ºK T=311 ºK
2 0 -4 -600
T=323 ºK
0
600
1200
1800
2400
3000
3600
Time (sec) Fig. 17. Time-dependent cumulative crack length under load levels: (a) 140 N and (b) 160 N.
0 0
3
6
9
12
15
Time (sec) Fig. 18. Time-dependent cumulative crack length under displacement rate: (a) 0.1 mm/min and (b) 0.2 mm/min.
588
A. Muliana, K.A. Khan / Computational Materials Science 41 (2008) 576–588
Table 4 Crack initiation and propagation time for mode II failure analyses
ble of predicting thermo-viscoelastic responses under general stress–temperature histories.
Temperature (K)
Crack initiation time (s)
Debonding time (s)
0.1 mm/ min
0.2 mm/ min
0.1 mm/ min
0.2 mm/ min
Acknowledgement
295 311 323
6.16 6.70 7.41
3.07 3.34 3.68
15.2 19.1 23.0
7.37 9.31 12.1
This research is sponsored by the National Science Foundation under Grant No. 0546528. References
higher load and temperatures accelerate the crack initiation and the complete debonding time. Table 3 presents crack initiation and debonding time at 70% and 80% failure load. 4.2. Mode II (shearing) time-dependent crack propagation A similar procedure as in the mode I is adopted for the mode II failure analysis. In this case, 40 meshes are generated along the bonded length for the adhesive layers and total 5 elements are used along the adhesive’s thickness. Cohesive elements are applied along the predetermined failure paths, which are located at the top and bottom interfaces between the adhesive layers and the aluminum adherend. A stress based failure criterion, Eq. (30), is used. This failure criterion corresponds to the separation between adhesive and aluminum. Fig. 18a and b show the crack propagation analysis for the mode II based on the critical stress. Analyses are run at displacement rates of 0.1 mm/min and 0.2 mm/min and at different temperatures. Softer material reaches the critical stress level slower than the stiffer material does, which cause slower crack propagation. The crack propagates equally in the anti-symmetric failure paths at the top and bottom interfaces. Table 4 shows the crack initiation and debonding time at different displacement rates and temperatures. 5. Conclusions A numerical algorithm is developed for the nonlinear viscoelastic material models of the TCM. The effects of stress and temperature are incorporated in the initial (elastic), long-term (equilibrium), and time-dependent (transient) compliances. The algorithm is compatible with the displacement based FE model for small deformations and uncoupled thermo-mechanical analyses. The recursive algorithm along with iterative schemes is performed at the material level, which enhances convergence and accuracy both at the material and structural levels. The numerical algorithm is verified with the available experimental data under general stress–temperature loading histories. It is shown that the proposed time-integration algorithm of a nonlinear viscoelastic model of the TCM typed is capa-
[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37]
David H. Allen, Chad R. Searcy, Int. J. Fract. 107 (2001) 159–176. Z.P. Bazant, Int. J. Numer. Meth. Eng. 4 (1) (1972) 11–15. W.C. Carpenter, Comput. Struct. 36 (6) (1990) 1141–1152. J.M. Caruthers, R.E. Cohen, Rheol. Acta 19 (1980) 606–613. J.M. Caruthers, D.B. Adolf, R.S. Chambers, P. Shrikhande, Polymer 45 (2004) 4577–4597. B. Chen, D.A. Dillard, Int. J. Solids Struct. 38 (2001) 6907–6924. F. Dubois, C. Chazal, C. Petit, Mech. Time-Depend. Mater. 2 (1999) 269–286. W.W. Feng, Int. J. Nonlinear Mech. 27 (4) (1992) 675–678. R. Haj-Ali, A. Muliana, Int. J. Numer. Meth. Eng. 59 (1) (2004) 25– 45. M. Henriksen, Comput. Struct. 18 (1) (1984) 133–139. H. Hilton, S. Yi, Compos. Eng. 3 (1993) 123–135. G.A. Holzapfel, G. Reiter, Int. J. Eng. Sci. 33 (1995) 1037–1058. B.D. Harper, Y. Weitsman, J. Rheol. 29 (1985) 49–66. A.R. Johnson, T.K. Chen, Comput. Meth. Appl. Mech. Eng. 194 (2005) 313–325. R.A. Jurf, J.R. Vinson, J. Mater. Sci. 20 (1985) 2979–2989. P. Krishnaswamy, M.E. Tuttle, A.F. Emery, J. Ahmad, Int. J. Numer. Meth. Eng. 30 (1990) 371–387. J. Lai, A. Bakker, Comput. Mech. 18 (1996) 182–191. A.H. Muliana, J.S. Kim, Int. J. Solids Struct., in press. D. Peretz, Y. Weitsman, J. Rheol. 26 (3) (1982) 245–261. D. Peretz, Y. Weitsman, J. Rheol. 27 (2) (1983) 97–114. H. Poon, F. Ahmad, Comput. Mech. 21 (1998) 236–242. H. Poon, F. Ahmad, Int. J. Numer. Meth. Eng. 46 (1999) 2027–2041. P. Rahulkumar, A. Jagota, S.J. Bennison, S. Saigal, Int. J. Solids Struct. 37 (2000) 1873–1897. S. Reese, P. Wriggers, ZAMM Appl. Math. Mech. 78 (1998) S157. S. Roy, J.N. Reddy, Comput. Struct. 29 (6) (1988) 1011–1031. Y. Sadkin, J. Aboudi, Compos. Sci. Technol. 36 (1989) 351–365. R.A. Schapery, Polym. Eng. Sci. 9 (1969) 295–310. R.A. Schapery, Compos. Mater. 5 (1974) 85–168. J.C. Simo, T.J.R. Hughes, Comput. Inelast., Springer-Verlag, New York, 1998. R.L. Taylor, K.S. Pister, G.L. Goudreau, Int. J. Numer. Meth. Eng. 2 (1970) 45–59. N.W. Tschoegl, W.G. Knauss, I. Emri, Mech. Time-Dept. Mater. 6 (2002) 53–99. M.K. Warby, J.R. Walton, J.R. Whiteman, Comput. Meth. Appl. Mech. Eng. 97 (1992) 375–397. A.S. Wineman, K.R. Rajagopal, Mechanical Response of Polymers, Cambridge University Press, 2000. S. Yadagirz, C. Reddy, T. Reddy, Comput. Struct. 27 (4) (1987) 445W. S. Yi, H.H. Hilton, M.F. Ahmad, J. Appl. Mech. 63 (1996) 218–224. O.C. Zienkiewicz, M. Watson, I.P. King, Int. J. Mech. Sci. 10 (1968) 807–827. M.A. Zocher, S.E. Groves, D.H. Allen, Int. J. Numer. Meth. Eng. 40 (1997) 2267–2280.