Nonlinear viscoelastic stress analysis—a finite element approach

Km-7949/&l 13 oat .oo Pergamon Press Lfd

Computers & Structures Vol. 18, No. I. pp 133-139, 1984 Printed in Great Britain.

NONLINEAR VISCOELASTIC STRESS ANALYSIS-A ELEMENT APPROACH

FINITE

MOGENSHENRIKSENt Department of Mechanical Engineering, Texas A&M University, College Station, TX 77843,U.S.A. (Receioed 2 August 1982;received for publication 19 January 1983)

Abstract-This paper describes a finite element algorithm developed for analysis of nonlinear viscoelastic materials. A single integral constitutive law proposed by Schapery is used to describe viscoelastic material behavior. Work leading to this paper focused on adhesives, but the FE formulation is general and readily extended to structural systems other than plane strain, plane stress and axisymmetric analysis as described. Cartesian strain components are written in terms of current and past stress states. Thus strains are conveniently defined by a stress operator that includes instantaneous compliance and hereditary strain which is updated by recursive computation. Equilibrium at each time step is insured with a modified Newton Raphson technique, incorporating convergence acceleration. Verification analyses show excellent agreement with experimental data for FM-73 adhesive systems. plane strain analysis of a butt joint is included

NOTATION

shift factor at time f for reduced time strain-displacement transformation matrix Kronecker delta elastic compliance creep compliance instantaneous compliance creep compliance coefficients in exponential series kinematic minus dilatational strain component at t kinematic strain component hereditary strain component stress operator vector of body forces vector of external forces material kernel functions relaxation coefficient for rth term in hereditary strain series material constant in creep compliance series matrix of stiffness coefficients for multiaxial stress matrix of compliance coefficients of multiaxial stress

Poisson’sratio reducedtime at time t reducedtime incrementfrom t -At to t hereditary integral component at time t for rth term exponential series index dummy variable-time Cauchy component of stress time time increment hygrothermal strain volume and differential volume Matrix symbols [ ] matrix I 1 column vector { }T transpose of column vector INTRODUCTION

Nonlinear viscoelastic materials are seeing increasing usage as structural components. Low manufacturing costs and design flexibility enhance selection of such materials. Automotive components, composities and adhesives represent three primary applications of high

tAssociate Professor.

strength polymers. Increased use of high strength composites by aeronautical and automotive industries has created renewed interest in the mechanics of such materials. Presently, research focuses on material evaluation, characterization and development of analytical techniques that are capable of accounting for historical effects of stress and environment. This paper is concerned with development of a nonlinear viscoelastic finite element code that is based on Schapery’s single integral constitutive equation. A stress operator that defines uniaxial strain as a function of current and past stress is developed. Extension to multiaxial stress states is accomplished by incorporating Poisson effects, resulting in a constitutive matrix that consists of instantaneous compliance, Poisson’s ratio and a vector of hereditary strains. The constitutive equations are of a form suitable for nonlinear finite element analysis. Plane stress, plane strain and axisymmetric formulations are included. The finite element method is the most powerful technique presently available for analysis of materially and geometrically nonlinear mechanical systems. Bathe et al. [I] provided a thorough review of non-linear time dependent analysis; their work emphasized the rigorous foundation on which nonlinear finite element analysis should be based. Oden[2] was among the first to suggest FE formulations for elastic, hypo-elastic and hyper-elastic materials systems. Others including Zienkiewicz[3], Gallagher[4], and Bathe and Wilson[5], Becker et al.[6] gave a clear presentation of the finite element method in general. Although not specifically addressing the FEM, Fung[71 presented the physical and mathematical foundation on which nonlinear analysis should build. The single integral law for viscoelastic materials systems was first proposed by Schapery[8]. The law is based on rational thermodynamics and, in its simplest form, defines linearly elastic stress strain behavior. Sanders and Haisler[9] evaluated the feasibility of predicting cyclic creep via the single integral law. Their results show some promise although refinements, for cyclic creep, are required. Flaggs and Crossman[lO] investigated effects of thermal and moisture transients in viscoelastic lap joints, using a linear elastic code in 133

134

M. HENRIKSEN

which elastic constants were determined from the single integral law. Adams and Schaffer[ll, 121 reported on residual process induced stress in fibrous composites. Their analysis[ll, 121 was based on the single integral equation and showed good agreement between experiments and analytical predictions. Zienkiewicz[3] suggested that linear viscoelastic materials could be modeled by “equivalent elastic constants”. Varying loads were permitted through time incrementation; iterative solutions of equilibrium equations led to acceptable results for linear viscoelasticity. A lack of supportive materials data, however, precluded choosing this alternative for the present work. Weitsman has published extensively[l3-151 on the mechanics of viscoelastic adhesive joints. Using variational methods and finite difference solution techniques, Weitsman analyzed a variety of idealized adhesive systems. Elasticity theory was useful in accounting for the presence of stress singularities but could not easily include relaxation phenomena. Peretz and Weitsman[ 161 recently released data that resulted from a materials characterization program of the FM-73 adhesive. Material properties and experimental data reported by Weitsman and Peretz are used in verification of this code for uniaxial loads and provide the data required in analysis of multi-axial stress states. CONSTITUTIVE

Stress

- MPa (Psi/ 145)

Fig. I. Kernel functions RO’.~I’and fi’ vs stress at 303°KFM-73 adhesive.

EQUATIONS

The single integral constitutive equation may, for uniaxial stress-strain, be stated as: Ef = ‘QID”U’ + g,’

I’

D, (+~‘,&‘n‘)ds. 0

,/,, I ,,,, I ,,,, (,,,, /

16

(I)

0

10

20

Stress

40

30

50

MPa (psi/ 145)

In eqn (I), the left hand term represents uniaxial kine-

matic strain at time t. The term Do is elastic compliance while D,(9) is a creep compliance function. go’ defines stress and temperature effects on elastic compliance and is a measure of state dependent reduction (or increase) of stiffness. Transient (or creep) compliance factor g,’ has an identical purpose, operating on the creep compliance component. Superscript, t, denotes current time, u* is Cauchy stress at t while s is a dummy variable of integration. gz’ accounts for influence of load rate on creep; g,” is also stress and temperature dependent. Functions %” and Ys are reduced time at t and s and are defined by

Fig. 2. Shift factor, a,,~‘. vs stress-FM-73 adhesive.

work, the exponential form is chosen since it permits hereditary effects to be computed recursively. Data reduced in either form may, of course, be converted to the other. Equation (1) can be expressed in the following operational form in which the subscripts denote operation on a Cartesian component of stress, u
ei;=F(uii’)

(4)

where F(oi;) is a stress operator, defined by eqn (5): where a:T is a time “shift factor”. This function physically modifies viscosity as a function of temperature and stress. Mathematically, a,,Ts shifts the creep data parallel to the time axis relative to a master curve for creep strain u. time. Figure I depicts, for the FM-73 system, the functional relationship between gO’,g,‘,gz’, and stress at a reference temperature of 303°K. These data were reported by Weitsman and Peretz[l6]. Similarly, Fig. 2 depicts shift factor avTSas a function of stress at three different temperatures. Transient creep compliance, D=(Y) may be expressed in either of two forms D,(T)= D,Q"= c D, [I - exp D,(q)

from creep recovery data.

this

F(a,j’) = D:aiif t E,i’ (i and j not summed). It may be shown that instantaneous given by

compliance

Dr’ = go’Do+ g,‘gz’c D,(l - I’,‘) r

(5) Dr’ is

(6)

where r,’ is a relaxation coefficient in the creep compliance series, defined as r,’ = (1 - exp [ - A,AV’])/A,A\Y

(7)

135

Nonlinearviscoelasticstress analysis-a finite element approach Hereditary strain component Eij’ may be written as

{a’} = [M']{e'}+{E'}

Ei/ = g,’ 2 (gJ-A’r;a;,:A’ - g,oai; I I

where [M’] is

- q :;r exp [ - LA’?‘])

(9)

I

1 [“‘l = Dr’(1t v’)(l -2~‘)

and the rth component, q:2A’, of the hereditary integral series at the end of the previous load step is 1-M 4 ,,ij =

r

{l - exp [ - h,(qt-A’ -V)]}

(10)

q :,ii = q :.iA’exp [ - h,Aq’] + (gz’giiif- g:-“’ aii-*‘)r,‘. (11) Equation (11) is derived in Appendix 1, using the assumption that product g2’uijS varies linearly over the current time step. This assumption limits time step size when stress states or temperatures change rapidly.

VL

1

$(g,‘ccj‘) ds.

Note that the integral on the r.h.s. of the equality sign of eqn (10) is identical to the rth term of eqn (1) when creep compliance is expressed in exponential series form. Only the integration limits differ. Integration is accomplished by updating each q :yiA’ as follows;

'

l-v’

1-A’ I -0

(17)

0

At

0

0

0

V’

V’

1-2v’ 2

0 (18)

For plane strain, kinematic strain component eX3’is identically zero; stresses may be computed by setting kinematic strain component E,,’ = 0. Hereditary strains are evaluated from eqn (9). For plane stress, only three independent strain components exist; kinematic strain esst can be evaluated by modification of eqn (18) as suggested by Bathe et af.[5]. SOLUTIONALGORITHM

The FE equilibrium equations may be established using the principle of virtual work as indicated in eqn (19):

hlULTIAXIALSTRESSFORMULATION

In order to formulate a stress strain relationship for multiaxial stress states, each strain component is assumed to be a linear function of the stress operators. Therefore, as in linear elastic analysis, the Poisson effect is incorporated. Coupling effects between stress components are assumed to be primary only and fully defined by the matrix relationship {e’}= D;[N’]{a’} t [N]{E’}.

(12)

In eqn (12), we strive to incorporate possible states of plane stress, plane strain and rotational symmetry, hence {e’} is a vector containing the algebraic difference of kinematic strains {E’}and thermal strains {a,#} {e’}T= {(El,’- e’),(E*a’- L0,Y,&33’ - 0’))

(13)

{su}T(~y~BIT~uT~d~-~Fcx’~)=IO~

in which {Su} is a vector of virtual displacements, and [B] is the strain displacement transformation matrix, defined by the relationship {de}= [B]{du}.

{8u}=(I

”

@‘IT = IE,,‘Ezz’E’~‘E~~‘}

(19

I

1

-ok’

;

-v’

-v*

0

-v*

l/2 0 0 1

.

(16)

Inversion of eqn (16) gives a constitutive relationship which, written in matrix form, is

v~[Blt{E’}dV-(f’..‘}}={O} I

(21)

or simply {SU}~{R'}={O}

0-v' [N'] =I-v* 1

[Nl is a matrix (4 X4), defined by eqn (16), in which V’is Poisson’s ratio:

[B]T[M’]{e’}dV-jy[B]T[M’]{e’}dV -

(14)

and {E’} is a vector of hereditary strains, components of which are defined by eqn (9).

(20)

The vector {u’} consists of Cauchy stresses, evaluated at internal integration points while CfteX’}is a vector of externally applied boundary forces. Application of eqns (13) and (17) lead to the following expression for equilibrium:

while {a’} contains four components of Cauchy stress {a’}T = {~,,‘u**‘u,z’u~s’ 1

(19)

(22)

where W’J=j

”

DITb’W-{f’}.

(23)

as expanded in eqn (21). In order for eqn (22) to be satisfied for arbitrary virtual displacements, the vector of residual forces, {R'} must be equal to zero. Since 4 is a function of reduced time, and therefore dependent on stress and temperature at the current time step, a direct solution cannot be obtained. Hence, an iterative solution, approximately satisfying (19) is sought.

M. HENRIKSEN

136

In determining a solution to (19), the following approach is followed: if {Rf,,} is the vector of unbalanced forces at iteration it 1 for time t, then {R:+,) = {Ri’}+

[KT’lIAui’}

where tangent stiffness [KT’] is given by

The solution is started by assuming a stress free state to exist. Iteration continues until the convergence tolerance is met. The problem is then restarted for the next load step, assuming a zero displacement state at the beginning of that load step. Thus a series of nonlinear problems is solved. Other integration schemes were investigated and found to offer no advantage over the chosen method which is unconditionally stable.

[KT’l=~ in which {u’} is a vector of nodal displacement at t. It can be easily shown, from eqn (21), that [KG] is given by [K,‘]=[

Y

[B]T[M*][B]dV.

(26)

Since we seek for a solution in which {Ri,,} = {0}, incremental displacements {Aui’}are evaluated from {Au;} = - [K;]{Ri’}

(27)

{u;+,}={ui’}+{Au,‘}.

(28)

and

The solution process for time step t proceeds as follows: (1) At the beginning of the time step, evaluate [KG]; (2) Decompose [K:] via an [L] [D] [LIT skyline structured equation solver; (3) Iterate through eqns (27) and (28) until the ratio of the norms of incremental displacements and total displacement satisfies < TOL

Material properties used in sample analysis are those reported by Weitsman and Peretz[l6]. Creep compliance data are given in Table 1 along with other pertinent material data. Thermal and stress dependent material functions (g0’,glt,g2’) as well as shift factor (a,<) are reported in the literature [ 161and repeated in Appendix 2. In order to account for multiaxial stress states, the data published by Weitsman and Peretz[ 161is assumed valid in octahedral shear stress form. Three uniaxial test cases (plane stress), corresponding with laboratory tests conducted by Weitsman and Peretz, were analyzed. In the first of these, a uniaxial stress of 10 MPa (1450psi) was applied to a specimen for 1200set, followed by a step increase to 26.6 MPa (3860psi) lasting an addition 1200sec. As may be seen in Fig. 3, excellent agreement between experiment and finite element prediction results. Figure 4 depicts stress vs strain behavior in uniaxial tension as both stress and temperature increase (from 0 to 20MPa and 303-333°K) linearly over a period of 3600 sec. Excellent agreement between finite element prediction and experimental data may be seen.

(29)

where TOL is a small value, typically 0.000 01. Solution accuracy does not depend on how well the tangent stiffness matrix represents the actual system but convergence rates are strongly affected. Therefore, occasional updates of the stiffness matrix are desirable. The modified Newton-Raphson technique requires a large number of iterations if rapid compliance changes in the system are encountered. In order to reduce computational effort a convergence acceleration scheme proposed by Crisfield[17] is used. A loadstep rarely requires more than ten iterations for convergence to a tolerance of 0.000 01.

Time (Seconds)

Fig. 3. Strain vs time-uniaxialtension at 323°K.

Table 1. Materialdata used in sample analyses FM-73 Adheeive - Unsctinmed (All data at reference temperature of 303OK. (66%) Property

SI units

Elastic compliance Do Creep con~llanee D

c

English

unita

360(&6)/KPa

2.48(E-b)/pal

29(&6)/MPa

2.00(E-7)/psi

Power law exponent "

.li

.12

Poisson's Ratio

.38

.30

Coefficient of thermal

expansion

6.6(E-5) mlml°K 11.9(E-5) in/idF

137

Nonlinear viscoelastic stress analysis-a finite element approach

Strain-Percent

Fig. 4. Stress vs per cent strain-uniaxial tension.

(a)

(b)

Fig. 7. Stress state at

0

2000

1000

3000

4000

Time-Seconds

Fig. 5. Per cent strain vs time-uniaxialtension stress: 10MPa

(1450).

Figure 5 shows strain vs time in uniaxial tension as temperature increases from 303 to 333°K over a period of 3600 sec. Satisfactory agreement between experiment and analysis is noted. The algorithm developed in this paper is designed to evaluate residual stresses in nonlinear viscoelastic adhesive joints. Such joints are generally cured near the glass transition temperature and subsequently cooled to room temperature. In assessing joint efficiencies, residual stresses must be evaluated. Example 3 treats a butt joint in which post cure stress relaxation occurs. Figure 6

Cc)

Cd)

= 120sec. for shaded sector of Fig. 6 (T = 303°K).

depicts the joint. Initial conditions assume the joint to be stress free at a curing temperature of 333°K. After the cure cycle, the butt joint is subjected to a temperature drop from 333 to 303°K that is assumed to occur linearly over a period of 120sec. In the analysis, the mechanical properties for FM-73 adhesives are used as listed in Table 1. The boundary conditions are imposed such that adherends are immobile and rigid. The latter assumption is reasonable in light of the low compliance of aluminum relative to that of FM-73 adhesive. The focus of the analysis is on assessment of stress relaxation that occurs over a period of approx. 3 yr. Figure 7 depicts the stress state existing immediately after the temperature drop to 303°K. Each of the cartesian stress components are depicted in Fig. 8 also, but at the end of the 3-yr period. It can be seen that stresses relax to approx. 50% of their previous magnitude in that time. Table 2 summarizes the data depicted in the post processing plots. In order for joint efficiencies to improve, post cure treatment is required. In this example, a relatively fine mesh of 6 x 12 &node isoparametric elements was used. The element adjacent to the singularity point is modified as suggested by Barsoum[l7] in order that the model might more adequately account for the singularity that is encoun-

x R~lon

Modeled

-I

I-

O.1Snlrn

0.006

(a)

(b)

Cd)

I”

Fig. 6. Butt joint subjected to cooling from 333 to 303°K.

Fig. 8. Relaxed stress state at t = IO’sec. for shaded sector of Fig. 6.

M. HENRIKSEN

138

Table 2. Summary of stress-states-Figs. 7 and 8

FigUKe

stress

COIltOUr MPa

Component*

Increment psi

la

20.9

3,030

200

lb

44.1

6,400

400

20.3

2,950

150

Jd

27.9

4,050

200

Je

51.7

7,500

500

12.4

1,800

100

25.5

3,100

200

11.4

8b 8c

0e

* Please

refer

co

1,650

100

20.0

2,900

100

29.6

4,300

250

Figure

6 for

tered in elastic analysis. The irregular stress contour near the singularity is produced by the post processor which does not include a singularity element. Effects of post cure heat treatment are not included in this sample analysis. The code, however, allows quantitative evaluation of stress reduction effects produced by various thermal cycles. CONCLUSIONS

The nonlinear material model incorporated in the finite element code can be extended to account for the finite strains that frequently exist in adhesive joints. The single integral law is held to be valid for finite as well as infinitesimal strains. Work is presently under way to extend the formulation to include high temperature cyclic creep. The algorithm described yields results that are in excellent agreement with experimental data. Usage is relatively economical. The code is implemented on a 32 bit minicomputer with a virtual memory operating system. Thus, large arrays required for storage of hereditary data are handled with a minimum of programming effart. Evaluation of alternate methods of integration of eqn (1) suggests the method described herein to be most economical among those studied.

Acknowledgement-The author wishes to express his gratitude to Profs. Y. Weitsman and R. A. Schapery. Our many discussions were extremely helpful in the execution of this project. REFERENCES 1. K. J. Bathe, H. Odzemir and E. L. Wilson, Static and dynamic geometric and material nonlinear analysis. Rep. UCSESM 74-4, Structural Engineering Laboratory, University of California, Berkeley (Feb. 1974). 2. J. T. Oden, Finite Elements of Nonlinear Continua. McGrawHill, New York (1972). 3. 0. C. Zienkiewicz, The Finite Element Method, 3rd Edn. McGraw-Hill, New York (1977).

coordinate

system

R. H. Gallagher, Finite Element Analysis, Fundamentals. Prentice-Hall, Englewood Cliff (1975). K. J. Bathe and E. L. Wilson, Numerical Methods in Finite Element Analysis. Prentice-Hall, Englewood Cliffs (1976). E. B. Becker, G. F. Carey and J. T. Oden, Finite Element Analysis-An Introduction. Prentice-Hall, Englewood Cliffs (1981). 7. Y. C. Fung, Foundations of Solid Mechanics. Prentice-Hall, Englewood Cliffs (1965). 8. R. A. Schapery, Further development of a thermodynamic constitutive theory: stress formulation. Purdue Res. Foundation, Project No. 4958,(Feb. 1969). 9. D. R. Danders and W. E. Haisler, An incremental form of the single integral nonlinear viscoelastic theory for elastic-plastic creep finite element analysis. ASME 79-PVP-I 14 (Aug. 1979). 10. D. L. Flaggs and F. W. Crossman, Viscoelastic response of a bonded joint due to transient hygrothermal exposure. Modem tures,

Deoelopments

in Composite

Materials

and

Struc-

pp. 299-314.ASME, New York (1978). 11. B. G. Schaffer and D. F. Adams, Nonlinear viscoelastic behavior of a composite material using a finite element micromechanical analysis. University of Wyoming, UWMEDR-COI-101-l(June 1980). 12. B. G. Schaffer and D. F. Adams, Nonlinear viscoelastic analysis of a unidirectional composite material. Trans ASME J. Appl. Mech. 48 (1981). 13. Y. Weitsman, Stresses in adhesive joints due to moisture and temperature. J. Composite Malls 11, 378-394(1977). 14. Y. Weitsman, Effects of fluctuating moisture and temperature on the mechanical response of resin plates. J. Appl. Mech. 44(4),571-576(1977). IS. Y. Weitsman, Interfacial stresses in viscoelastic adhesive layers due to moisture sorption. Int. J. Solids Structures 15 (1979). 16. D. Peretz and Y. Weitsman, The non-linear thermo-viscoelastic characterizations of FM-73 adhesives. J. Rheology M,245-261(1983). 17. M A. Crisfield, A faster modified Newton-Raphson iteration. In Computer Methods in Applied Mechanics and Engineering. pp. 267-278.North-Holland, Amsterdam (1979). 18. R. S. Barsoum, On the use of isoparametric finite elements in linear fracture mechanics. Int. I. Num. Meth. Engrs 10,25-37 (1976).

139

Nonlinear viscoelastic stress analysis-a finite element approach series, given by

APPENDIX1

Derivation of recursive equations Equation (1) is restated in exponential form for component (no summation is implied).

(1.6) ij

,

d Ei/ = go’Doailtt g,’ HD,[1 - exp {- A,(Y’-Y’“)}]- (g2’gijs)ds. ds I -0 (1) Let the product gzsui; be expressed by variable Gi/. Equation (I) may be rewritten as:

I

I -0

exp [ - A@ - W)] 9

ds.

The remaining components represent hereditary strains, Eii which are given by

,-A1 1 _. exp [ - A,(@-OS)]‘2 I

After manipulations of eqns (l), (1.4) and (1.5) it can be shown that terms in the hereditary integral at the end of the current time step is given by q~,~~=exp[-A,AY’]q~~ALt(G~~-G~,~A’)~~,

exp [ - h,(V’ -‘V)] 2

ds

ds.

(1.2)

Shift factors and material Kernel functions For the sake of convenience, shift factor equations and kernel functions, as reported by Weitsman and Peretz[l6] are included here: Material kernel functions are go’= (I t @)(I t jf,

Note now that the first integral on the r.h.s. eqn (1.2) may be rewritten by the use of basic properties of logarithms as I+Pf exp [ - A# -0 I where q:,iA’is given by I-& I -0

exp [ - A#PmA’

- Ur”)]

qds ds

6’ = ;{(c~;)*t (a;)* t (oi’)*- 6(7,,1)2}“2/a,,

(2.2) (2.3)

exp[-A,(Y’-Ys)]~ds=(Gi~-G~,~A’)~,’ 1-81

T’ = T’ - To TO

where T, is temperature in “K and stress is given in MPa. Constants in the material kernel functions are as follows: .

(1.4)

j=O.lS,

The second integral of eqn (1.2) is integrated by assuming Gii varies linearly over the current time step. Therefore,

I’

exp [ - @‘I

ds = exp ( - A,AY’)q:LF’ (1.3)

,-A, 4 r.ij =

gl’ = {I t 4(#4}

(2.1)

g2’= {l t fi($)‘} exp [fit’]

- 9”)] 2

(1.9)

APPENDIX2

t I ,-A,

4 i;y,

(1.8)

exp [ - h,(q’ - ‘us)]2 ds

t

(1.7)

(1.1)

Now break a series term in integral (1.1) into two parts, the first part having limits from - 0 to t -At, the second integral spanning only the current load step.

I’-0

D:=golDotg,‘gl’CD,(l-rr:).

Eii’= g; x D,{ g imA’ r!Ui,‘Af - g2’uij@j - exp[ - &A@] r

E/ = g,,‘D& t g,’ 2 D, _$ Gii’ds , It

-g,‘cD,

Collecting those terms in eqn (1) that are multiplied by current stress yields instantaneous compliance:

p^=O.91, q=1.435,

4=8.5,

p=O.75, fi=12.12

~0 = 23.6MPa. The shift factor a,T’ is reported as

(1.5)

aoTt= exp [ - (7~7’ t it’)] where

where Tr’ is a relaxation coefficient in the transient comdiance

(2.4