Large displacement analysis of viscoelastic shells of revolution

COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING 3 (1974) 237-253 Q NORTI~-HOLL~D PUBLISHING COMPANY

LARGE DISPLACEMENT ANALYSIS OF VISCOELASTIC SHELLS OF REVOLUTION P.K. LARSEN and E.P. POPOV Department

of CivilEngineering, Universityof California Berkeley, USa

Received 16 April 197 3 Revised manuscript received 10 August 1973

The incremental equilibrium equations for a body subjected to iarge deformations are given in a Lagrangian description. The loading terms due to both conservative and nonconservative types of loading are derived using virtual work considerations. The simple integral constitutive laws of the infinitesimal theory of linear viscoelasticity are generalized to small strain, large rotation applications through the use of invariant strain and stress measures. The viscoelastic material is characterized using a Prony series expansion for the relaxation modulus, leading to a simple algorithm for the evaluation of the convolution integrals. The shell geometry and displacement fields are discretized by the “degenerate” isoparametric shell element. The numerical examples include the post-buckling behavior of an elastic shallow spherical shell and the viscoelastic creep buckling of a spherical cap subjected to sustained loading.

1. Introduction Plastics and other viscoelastic materials have been used for nonstructural purposes for a number of years. During the last decade, however, their use in load carrying structures has increased significantly. High polymeric materials are extensively used today as light weight cores in sandwich structures and in many highly stressed products serving in exposed environments. An example of such an application is the use of plastic windows in deep sea exploration pressure vessels. Here the viscoelastic character of the material must be considered both from safety and functional aspects, since the optical characteristics of the window are deformation dependent. The theory of linear viscoelasticity for infinitesimal deformations is well known, and solutions to boundary value problems for such materials have been presented in the literature. The finite element method has also been used extensively to obtain solutions to this class of problems [ 1,2l. The nonlinear theory of viscoelasticity is significantly more difficult, and few problems have been solved for this case. Here the nonlinear stress-strain transformation is commonly given in the form of multiple convolution integrals, and the identification of these kernels is difficult. Furthermore, such laws are expensive computationally since few effective recursive algorithms are available. These complications were the motivation for the development of the theory of finite linear viscoelasticity, which is a generalization of linear theory suitable for small strain, large rotation problems 13, 41, In this theory the Stress is assumed linear in the strain, but the kinematics of the deformation is nonlinear. Such generalized theories have in a few cases been used for the large displacement analysis of shells f 5, 61. The response of a viscoelastic material depends on the stress history, and the governing equations must hence be integrated along the deformation path. Incremental solution techniques must therefore be used. Such techniques for a body undergoing finite deformations have been derived

238

I? K. Larsen, E.!Ip.Popov, Analysis of viscoelastic shells of revolution

previously and have been implemented for the finite element analysis of both geometric and material nonlinea~ties [ 7, 8, 91. Herein the incremental formulation for the nonlinear finite element analysis using the Lagrangian description is presented. The traction bounda~ conditions for both ~onse~ative and nonconservative loading are given. The unbalanced load due to the lack of equilibrium between the externally applied load and the internal stress Geld is included. Various generalizations of the linear theory of viscoelasticity are discussed and evaluated in view of their computational feasibility. The form chosen for this application is given in terms of Lagrangian strain rate and the 2nd Piola-Kirchhoff (P-K) stress tensor. Isotropic materials are characterized by their bulk and shear relaxation moduli, here represented by Prony series expansions. The shell is discretized using an axisymmet~~ finite element of the “degenerate” so-called isoparametric type [ lo]. Parametric interpolations are used both for geometry and displacement fields, and the kinematics of the deformation is determined by the theory of a three-dimensional solid. The derivation of the linearized tangent stiffness of the element, consisting of the incremental and geometric stiffness matrices, is presented in matrix form. The load vectors of applied loading and creep pseudo-loading are likewise given. The numerical examples include the post-buckling analysis of a linearly elastic shallow spherical cap and the viscoelastic creep buckling analysis of a viscoelastic spherical shell. Some problems inherent in the use of the correction for “unbalanced” load in viscoelastic boundary value problems are discussed.

2. Incremental equilibrium equations During its deformation a body will describe a continuous motion from an initial configuration B, to the final ~onfi~ration B. Let an intermediate ~~~~ent) configuration be B, and its neighboring one B, (Fig. 1). The equations governing the motion of this body can be written in two alternative ways. The direct formulation is given in terms of the initial configuration B, and B without considering the deformation path. The resulting equilibrium equations are nonlinear and must be solved by itera-

Fig. 1

P.K. Larsen, E.P. Popov, Analysis of viscoelastic shells of revolu?ion

239

tion. The incremental formulation, on the other hand, describes the motion from B, to B, , and the final solution is then obtained by forward integration along the deformation path. For boundary value problems involving materials with memory, the deformation path must be considered, and only the incremental formulation is computationally feasible. Depending upon the choice of kinematic description, the resulting equilibrium equation will take on different forms. The most common descriptions used in solid mechanics may be referred to as Lagrangian and Eulerian. The former refers to the initial confi~ration B, , and the latter most commonly to the current configuration B, . The choice of description depends upon the solution method to be used and the mathematical structure of the const~tutive laws for the given material. In the context of the finite element method using refined elements, the Lagrangian mode is the more advantageous [ 11 I. This is normally the case even when the constitutive laws are more conveniently formulated in a different description. The incremental equilibrium equations in a Lagrangian description are well known [ 8, 91, and only a few steps in their derivation are shown here. For quasistatic problems the inertia terms may be neglected, and this leads to the following virtual work expression .for the body in configuration B 1 . 6W, =

$

(624)’(9) da +

34

Jp(W(“f) du.

Similarly one has in confi~ration SW, = f a4

(1)

4

tsuY c2t)fzz +

B,

~~(~~}j~2~~ du,

(2)

4

where ‘t, %, ‘f, and 2f are traction vectors and body forces in B, and B, f and 3Bi refers to the surface of B, where surface tractions are prescribed. Moreover, &a, du, p are the infinitesimal surface efement, volume element, and mass density in B, , and the barred quantities are the equivalent elements in 3,. Using the principle of conservation of mass, the volume integrals in Eqs. ( 1) and (2) can be transformed to integrals over the volume in configuration 43,. Cauchy’s theorem in terms of the P-K stress tensor ‘S and defo~ation gradient ‘F in 3, is given [ 121 by

Ot=‘SN and

‘tda, = fF"tdA, where da, takes the value du or d5, and N and 64 are the normal vector and infinitesimal area element in B, respectively. Using Nanson’s formula and the Green-Gauss theorem, these relatiosnhips facilitate the transformation of the surface integrals in Eqs. (1) and (2) to configuration B, I1 1,181. For the sake of simplicity, rectangular Cartesian coordinates are used when component forms are exhibited; capital indices indicate that components are given relative to B,.

240

I? K. Larsen, E.P. Popov, Analysis of viscoelasticshells of revolution

and 6W, = j$$u) 30

(2F2S)dV=

j” =SJL!&dV, BO

(4)

where (a,) is the right divergence operator and d V is the volume element in B, . Also E is the strain increment between B, and B, , and e is the part of E that is linear in u

EfJ = “EIJ - ”ErJ= e, + qIJ, where 7 is the nonlinear part of E. Using rectangular Cartesian coordinates,

and

Combination stress field

of Eqs. (3) and (4) yields the incremental virtual work equation for the internal

The increment in P-K stress is given by S = 2S - “S and can be expressed in terms of the strain increment E by the linear tensor transformation C S=CE.

(7)

S is nonlinear in the displacement increment ~1,and substitution of Eq. (7) in Eq. (6) will render the latter nonlinear in u. A linearized version of Eq. (6) is obtained when E is replaced by e in Eq. (7). This leads to the well-known expression for the system stiffness consisting of the incremental and geometric stiffness. The incremental virtual work expression for the prescribed tractions and body forces may be obtained directly by subtracting Eq. (1) from Eq. (2). However, the use of Eq. (7) in stress computations combined with the linearized version of Eq. (6) will give an internal stress field that does not necessarily equilibrate the applied loading. This nonequilibrated force (residual load correction) can be accounted for in the load increment by combining Eqs. (2) and (3)

J? K. Larsen, E.P. Popov, Analysis of viscoelasticshelIs of revolution

241

This general expression must be specialized for conservative and nonconservative loading. For nonconse~ative loading of the pressure type the traction vector 2t is given by 2t&j=

po ( 2 F-‘)‘NdA. --2p~d2 = -2p T fJ

where 5 is the outward normal to the surface element dB in B, , and 2p is the pressure on this element. The transformation from B, to B, is accomplished by using Nanson’s formula. 2F-1 is given in terms of gradients with respect to confi~ration B, and hence cannot be determined directly. Furthermore, it is linear in the displacement increment U, and substitution of Eq. (9) into Eq. (8) will hence give rise to an additional stiffness term [9]. It should be noted that this term is nonsymmetric. An approximation to Eq. (9) is obtained when the displacement increment u is assumed to be small, and one can assume P+=P

and

a . a -=----* ax, . ax,

This leads to the approximate virtual work expression

The second term in the parenthesis in this equation represents the additional stiffness term. For most en~nee~n~ applications the norm of this term is small compared to the norm of the stiffness terms of Eq. (6). Since the residual load correction term is included, the benefit of computing the tangent stiffness matrix to a very high degree of accuracy by including this term is therefore negligible compared to the added computational effort needed to solve the equilibrium equations with an unsymmetric stiffness matrix. For conservative loading the applied load is usually given in terms of vectors with fixed directions in space. Since there are no rules for the transformation of such vectors from one configuration to another, the load must be redefined in terms of vectors that do transform. A convenient set of such vectors is the convected base vectors 2G, to the coordinate syrface X1= const. in B, which transfroms as “Gz = (2F-1)JI”GJ,

where “GI is the set of base vectors to X1 = const. in 13,. For gravity loading this gives

The norm and direction of the gravity load vector 2q are known, but its components relative to 2G are deformation dependent and must be computed in configuration B, . As before, the deformaiion gradient 2F-1 cannot be evaluated directly, and approximations must be introduced, as

242

for the nonconse~ative arly to Eq. ( t 0).

R K. Larsen, E.P. Popov, Analysis of viscoelastic shells

of revolution

loading, The final external virtual work expression can be obtained simif-

3. Constitutive theory for viscoelastic materials The theory of linear viscoelasticity for infinitesimal deformations is well known and has in the past been applied to both isothermal and nonisothermal boundary value porblems. The response oii to a given strain history e&t) is here given by the convolution integral

(12) where E.. and 0.. are the infinitesimal strain and stress. The integrating function G,, reduces to two independent components for isotropic materials and exhibits the fading memory of the material. The nonlinear theory of viscoelasticity has many features in common with the linear theory, but more precise measures for stress and strains must be used. This also means that a suitable mode of description must be chosen for the constitutive law. The most general nonlinear theory gives the free energy as a functional of the strain history, usually in the form of multiple integrals, and derives the stress response from this functional. This leads to a formulation where the response cannot easily be given by a recursive algortihm and is hence expensive computationally. The theory of finite linear viscoelasticity was derived for the case where the stress is linearly dependent on the strain history but where the kinematics of the deformation is nonlinear. For this case a generalization of Eq. (I 2) would simplify the formulation and reduce the computational effort. Truesdell and No11 [ 131 give such a generalization using the current configuration as reference; however, their formulation requires the polar decomposition of the deformation gradient in order to compute an invariant stress measure and is computationally unattractive. An alternative formulation using convected coordinates has been used for such problems [ 14, 151. Here the mathematical structure of Eq. (12) is retained, but oii and :ii are referred to convected base vectors. This formulation is energy conjugate and satisfies the invanance requirements. However, it should be noted that the integrating function G,, gives a relationship between physical components of stress and strain as determined by experiments. Since the physical components of the convected base vectors are deformation dependent, additional transformations must be introduced when the tensor components in Eq. (I 2) are replaced by their physical components. In order to avoid this problem a pure Lagrangian generalization of Eq. (12) is chosen here. The symmetric P-K stress tensor and the Lagrangian strain tensor are used, and Eq. ( 12) is interpreted as a functional relationship between the physical components of these. (13) where ‘SIJ is the stress in configuration B, at time t,. Note that E,,(T) is the total Lagrangian strain at time r and not the strain increment as defined in Eq. (5).

P. K. Larsen, E.P. Popov, Analysis

of viscoelastic skells of revohtion

243

G~~~~ should be identified from experiments at finite strains. However, for small strain, large rotation applications, the components of GIJKL may be approximated with sufficient accuracy by the co~esponding components G,, as determined from infinitesimal strain tests.

4. Material characterization For nonaging viscoelastic materials the simplest possible response algorithm is obtained for generalized Maxwell-type materials. Here the relaxation modulii may be represented by a Prony series am

= G,” +

z

Gkemtfxi,

(Y= 1,2,

f=1

where in the isotropic case G! and G, are the bulk and shear moduli, Gz are the equilibrium values, Xi the discrete relaxation times, and f the number of Maxwell elements in the representation. The advantage of this formulation is that only the solution at the previous time step is needed for the evaluation of the convolution integral in Eq. ( 13). Using the recursive algorithm of Taylor et al. [ 23, the mean stress @{t,> at time t, is given by i$t,)

=

G’:+

6 G;hi(At,)] + da(t,)

,f=1

G;eltn_,>

+6 g,(t,), f=l

where

and

The mean strain is assumed to vary linearly within the time interval (fj_, , tj), and Atj = ti - tj_l , dZ($) = P(fj) -

z(tj_l).

The increment in mean stress going from configuration B,_r to Bn is hence dii(t,) = @(t,) - rQ&

=

(1.9

244

P. K. Larsen, BP. Popov, Ann&sis of viscoelasticshells of revolution

The stress deviator is given by a similar equation obtained by replacing G, by G, and mean strain and stress by deviator&z stress and strain. The stress increment can now be obtained easily from the mean and deviatoric stress. For materials where Poisson’s ratio remains constant during deformation, v(a) = v(O), the above procedure is simplified, since E)may be kept outside the time integration. In this case Eq. (14) will represent Young’s modulus. Eq. (16) gives a relationship between increments of stress and strain in the form (17) where sIJ is the delayed response. From Eq. (16) it can be seen that CrJKL is a function of the time steps which will be constant only when equal time steps are used in the analysis.

Finite Element Formulation 5. Representation

of geometry and displacement field

Various plate and shell theories were originally developed in order to reduce the complexity of the general three-dimensional theory of elasticity. By considering the structure to be thin in the normal direction and introducing Kirchhoff s assumption, the general theory was reduced to a two-dimensional problem. However, for finite element formulations this procedure is of less benefit, as it severely complicates the construction of admissible elements that satisfy the continuity requirements. In the present application the “degenerate” isoparametric shell element [ 101 was therefore chosen. Here the element is considered as part of a solid, allowing the use of the simple straindisplacement relationships from the three-dimensional theory of elasticity. Complicated kinematic relationships of the shell theories are thereby bypassed. These elements give good results for moderately thick shells but are known to be in error for thin shells unless special inte~ation schemes are used [ 16, 171. This is caused by the appearance of shear deformations in pure bending modes for the lower order elements. This effect may be alleviated by using a higher order element with a cubic displacement field. For a cubic axisymmetric element the shell geometry is represented by a coordinate transformation (Fig. 2).

(18) where cg,({) are interpolation polynomials, h, the thickness, and 8, the angle between the r-axis and the shell “normal” at node “i”. In addition, g and q are the local natural coordinates defined by the shape of the element such that $ = T 1 at nodes 1 and 4, and q = + 1 on the outer and inner face of the element respectively. The kinematics of the element are based on the following assumptions: Plane sections initially normal to the middle surface remain plane but hot necessarily normal. The displacement normal to the middle surface is constant through the thickness.

P.K. Lmsen, IS? Popov, Analysis of viscoelastic she& ofrevolution

245

Fig. 2.

The same displacement field representation ment zcand the total displacement *u

may now be used for both the displacement incre-

(19) where h, is the thickness at node i. The three degrees of freedom at each node are displacements ui and wi and the rotation of the shell normal (Ye(Fig. 2). For the cubic elements the nodes are located at ($ = -11 -l/3, +1/3, +l), and the interpolation polynominals are given by

6. Kinematics of the deformation In a local orthogonal coordinate system (s, t) (Fig. 2) the linear part of the Lagrangian strain increment between configurations B, and B, can be written as

5s’

t =

eee e ~_ St

1+-

aQ4,

as

)

0

a%$

0

_“~

as 0

--l

a;

1_

(I+% )

r au,

at

au2 as

t

(20)

246

P. K. Larsen, E.P. Popov, Analysis of viscoelastic shells of revolution

where u1 and uz are the displacements relative to this local system, and u is the global displacement in the r-direction. From the assumed displacement field it follows that &.l, at

% _=:

and

= 0,

at

0.

For a state of plane stress, the normal strains ‘e,* and ett wili be determined from the constitutive law and not from kinematics. The term (1 + 3’uz/at> is therefore replaced by (I + let,), where reit is the physical component of the normal strain as computed from the plane stress condition. In matrix form this gives e=Au.&,

(211

where the vector of displacement gradients u, is transformed from the displacement increment u by the matrix operator B u, =Bg,

B=W,

??I.

(22)

Combining Eqs. (2 1) and (22) e = A&i.

(23

The matrix B consists of gradient operators with respect to ~onfi~urat~on B, and remains unchanged during deformation, while the effect of the previous deformation is contained in A. The displacement increment u consists of the global displacements u and w. The nonlinear part of the strain increment relative to the same local coordinate system is given as _ _ rl SS

u 2

0

(241

r

rl%%

au, au, rl St

_

-

_

“cat

_

Symbolicahy this can be written in matrix form as

cm where

P. K. Larsen, E.P Popov. Analysis of viscoelastic sheIfs of revoltrtion

I$

= H&

q) =

B*R,B

,

247 (26)

with

The Hli are given only in terms of B and will hence remain unchanged throughout the deformation path. Due to the parametric representation of both the geometry and the displacement field for the element used here, the elements of B cannot be given explicitly and must be computed numerically.

7. Stiffness matrices and load vectors The incremental equilibrium equations governing the deformation from configuration B, to B, are obtained from the linearized version of Eq. (6) and Eq. (8) using the stress-strain law given by Eq. (17). (KO

+ K,)u=Rp-RE-RD

.

(27)

When the load “stiffness” in Eq, (10) is neglected, the incremental and geometric stiffness matrices K, and Ko are given by

%

=

2n j’ s’ -1 -1

‘S’H’

r(r$, q)

dqd&

(29)

where

H*

=

He6

I %t

+ HL

.

I

The matrix operation in Eq. (29) implies that the submatrix U;, should be multiplied by the scalar element ‘Sss of ‘S and so forth. The nodal load needed to equilibrate the internal stress field is % = 2n jl -1

j’ -1

“PAI? r(& rj)dqd[

,

(30)

P. K. Larsen, BP. Popov, Analysis of viscodastic shells of revolution

248

and the pseudo-loading due to the delayed stress response is jr %iBr(&

R, = 277 j1 -1

rj)dgdt

.

(31)

-i

Finally, the applied pressure-type loading is given by +1

RP = 2x

s

-2$$-

(32).

[‘F-“]?Vr(t)dE,

-I

where N is the vector containing the components of the outward normal. Both stress and kinematic quantities in Eqs. (27) through (32) are functions of both t and q. However, by proper decomposition of these quantities, the integration over the shell thickness can be performed separately saving substantial effort in computation.

8. Solution method The repeated solution of Eq. (27) by Gaussian elimination constitutes the forward integration of the quasistatic version of the equations of motion. The load term R is obtained from Eq. (32) includes by giving the load at any time as a percentage of the full load, The difgrence R -R both the increment in external load and the nonequilibrated load at the end of f;he p&ious step. The formulation hence implies that a one-step Newton iteration is performed on the latter components The load factors and length of time steps are originally given as input. For elasticSly applied loads the load factors are modified as the determinant of the stiffness matrix changes near the critical loads. Likewise, the step length in the forward integration in time is changed as the rate-ofchange of displacement varies with time where a characteristic displacement component or displacement norm is used to monitor the behavior. Numerical examples. 1. Post-~u~~~~~~behavior

ofshallow ~p~er~ca~shell

The accuracy of the large displacement aspect of the program developed using the preceding formulation is illustrated by the post-buckling analysis of a shallow spherical cap. The unstable part of the load-deflection path is established using load reversals. Since the local extremum values are not know a priori, a simple test must be devised for the load reversal decision. In the present example the loading procedure is governed by the positive definiteness of the tangent stiffness matrix K, + K,. The determinant of this is the trace of the triangularized matrix obtained during the Gaussian elimination process, and it will indicate whether the confi~ration is stable or not. For more details regarding the loading/unloadin~ strategies, see Ref. [ 181. The geometry of the shell is as follows: R = 50.0 inch,

h = 0.1 inch,

CY = 4.48”,

P. K. Larsen, E. P. Popov. Analysis of viscoelastic shells of revolu tion

I

05

I

I

1.0

1.5

Dlsplocement

249

b I

2.0

I

2.5

I

30

w/h

Fig. 3.

and the linearly elastic material has E = lo6 psi and v = l/3. The applied load was 40 psi, and the normalized load parameter is

where ??I4= 12(1-V2). The shell was analyzed using 5 elements and a total of 36 load increments. Of these 8 were used to determine the upper critical load and the remaining to determine the post-buckling behavior. Even with these few load increments the results agree quite well with Weinitschke’s solution [ 191, as shown in Fig. 3. The computer time needed for the complete analysis was 0.3 set (CP time, CDC 6400) per element and load step. 2. Viscoelastic creep buckling of spherical shell The behavior of a viscoelastic shell subjected to sustained pressure loading is illustrated for a shallow spherical shell. The objective of this study is to determine the life expectancy of the shell when subjected to various pressure levels.

P K. Larsen, E. P. Popov, Analysis of viscoelastic shells of revolu tion

250

The geometry of the shell is identical to that of a plexiglas window for a deep submergence pressure vessel tested by the U.S. Navy: R = 31.57 inch,

h =

2.5 inch,

CY = 36.0”.

The material is a PMMA, Plexiglas grade G with instantaneous Young’s modulus E = 4.5 X lo5 psi and Y= 0.35. Poisson’s ratio is assumed time-independent for this study. Relaxation data obtained from creep tests are given in Ref. [20] . The data are given only for a,stress level of 1000 psi, and this is assumed to be representative for the present study. A four-term Prony series is used to represent the relaxation modulus G(t) = G, +

5 Gieetixi; i=l

choosing hi one decade apart, the Gi are determined from the relaxation data using a least square fit. A simple search was made to pick Xi that minimized the least square error in G(t). This procedure will in most cases give a sufficiently smooth approximation for G(t). Based on this analysis the following relaxation times and moduli are used. i

Gi (psi)

hi (hr)

0

1.25

x lo5

0

1 2 3 4

0.9765 0.8102 0.5776 0.8939

x 10’ x lo5 X 10’ X lo5

20 200 2000 20000

The behavior of the shell under instantaneous

APEX

pressure loading is given in Fig. 4. The shell ex-

DISPLACEMENT

Fig. 4.

w iin

i

PK. Larsen, E. P. Popov, Analysis

of viscoelastic

shells of revolution

25 I

hibited snap-through behavior with upper and lower critical pressures of 1885 psi and 1780 psi respectively. Six elements and 23 increments were used to obtain the load-deflection relationship. The creep analysis was made with pressure levels of 1 130 psi, 1320 psi, and 15 10 psi. This represents respectively 60%, 70%, and 80% of the upper critical loading. Pressure was applied instantaneously using 12 load increments and then sustained during the creep process. Fig. 5 shows the apex displacement versus time for the three cases. As seen, only the case p = p/p,, = 0.8 exhibits snap-through with sudden loss in stiffness. For p = 0.7 and p = 0.6 th e “snapped” configuration is reached after a gradual increase in displacement over a substantial period of time. For the case p = 0.8 oscillations in the solution were observed starting at time r = 60 hr. These oscillations are due to the use of the out-of-balance force RP - R, in the load vector. As the critical configuration is approached this force may cause the lack of equilibrium to be overcorrected, and at the end of the time step the stresses obtained are too small. In the next step these stresses must again be corrected for equilibrium; and when combined with the fading memory of the material, oscillations are set up. The use of smaller time steps will not solve the problem, since this may give even more erroneous variations in stress with time. Instead the out-of-balance force should be applied instantaneously (i.e., elastically) using the tangent stiffness from the previous time step and the corresponding strain increment added to that caused by the creep pseudo-load. This gives a smoother variation in the stress and reduces oscillation. 20 ,

I

I

Time ( hours ) Fig. 5.

I

I

252

P. K. Larsen, E. P. Popov, Analysis of vi~~oeiastic shells of revolution

9. Summary and ~onc~usious

The basic virtual work expressions for a body undergoing finite deformations have been used to obtain the incremental equilibrium equations. Using well-known transformation rules it has been shown that both conservative and nonconservative types of loading give rise to nonsymmetric stiffness terms in a Lagrangian formulation. The tangent stiffness matrix for the step forward inte~ation of the equations of motion is obtained by linearization of the incremental equations. A convenient generalization of the integral ~onstitutive laws of infinitesjmal theory of linear viscoelasticity is given in terms of the symmetric Piola-Kirchhoff stress and Lagrangian strain increments. Furthermore, by using a Prony series approximation for the relaxation modulus a very efficient recursive algorithm is obtained for the evaluation of the convolution integrals. The formulation as applied to axissymmetric shells of revolution is shown to give excellent results for the post-buckling behavior of a spherical cap. When applied to viscoelastic materials, it has been observed that the use of the out-of-balance load may lead to oscillations in the solution unless special care is used in the computation of the strain increment for each step.

Addendum

In connection with the last sentance in the first paragraph of the section on Finite Element Formulation, the author’s attention has been called to some pertinent references. Successful admissible elements satisfying the continuity requirements have been ~onst~cted; these are the t~an~lar shell elements SHEBA with arbitrary geometry in space, and the axisymmetri~ shell elements SABA [ 2 l-241 . A generalization including transverse shear strain effects is given in Ref. 1251. References [I] OX. Zienkie~cz, M. Watson and I.P. Ring, A numerical method for visco-elastic stress analysis, Int. J. Mech. Sci. 10 (1968). 807-827. 121 R.L. Taylor, KS. Pister and G.L. Goudreau, Thermomechanical analysis of viscoealstic solids, fnt. J. Num. Meth. Eng. 2, (1970) 45-59. [3] B.D. Coleman and W. NOR, Foun~tion of linear ~scoelasti~ty, Review of Modern Physics 33 (1961) 239-249. 141 A.C. Pipkin, Small finite deformations of viscoelastic solids, Review of Modern Physics (1964) 1034-1041. [5] N.C. Huang, Ax~y~etric creep buckling of clamped shailow spherical shells, Journ. of Applied Mechanics 32 (1965) 323-330. [6] Z. Bychawski, Some problems of creep bending and creep buckling of viscoelastic sheet panels in the range of large deflections, Non-Classical Shell-Problems, Proc. of IASS Symposium, Warsaw, 1964 (North-Ho~and, Amsterdam, 1964) 368-383. [ 71 S. Yaghmai, fncremental analysis of large deformations in mechanics of solids with applications to axisymmetric shells of revolution, Ph.D. Dissertation, University of California, Berkeley, 1968. NASA CR-1350, June 1969. [ 81 J.T. Oden, Finite element formulation of problems of finite deformation and irreversible thermodynamics of nonlinear continuua - a survey and extension of recent developments, Japan - U.S. Seminar and Mathematical Methods of Structural Analysis and Design, Tokyo, August 1969. R.H. Gallagher, Y. Yamada and J.T. Oden (eds) (Univ. of Alabama Press, 1970). (91 H.D. Hibbit, P.V. Marcal and J.R. Rice, A finite element formulation for problems of large strain and large disp~a~ments, Int. J. Solids Struct. 6 (1970) 1069-1086. [ 1O] S. Ahmad, B.M. Irons and O.C. Zienkieticz, Curved thick shell and membrane elements with particular reference to axisymmetric problems, Proc. 2nd Conf. on Matrix Meth. in Struct. Mech., AFFDL-TR-68-150 (1968) 5399572. [ 111 P.K. Larson and E.P. Popov, A note on incremental equilibrium equations and approximate constitutive relations in large inelastic deformations. To be published in Acta Mechanica.

PK. Larsen, E.P. Popou, Analysis of viseoeiasticshells of revolution

253

[ 121 L.E. Malvern, Introduction to the mechanics of a continuum medium (Prentice Hall, 1969). [ 131 C. Truesdell and W. Noll, The nonlinear field theories of mechanics, in: S. Fliigge (ed), Handbuch der Physik III/3 (1965). [ 141 A.G. Fredrickson, Principles and applications of rheology (Prentice Hall, 1964). ]15] J.T. Oden, A generalization of the finite element concept and its application to a class of problems in nonlinear viscoelasticity, in: Developments in Theoretical and Applied Mechanics IV, D. Frederick (ed.) (Pergamon, 1968) 581-592. [ 161 SF. Pawsey and R.W. Clough, Improved numerical integration of thick shell finite elements, Int. J. Numer. Methods Eng. 3 (1971) 576-586. [ 17) 0.Z. Zienkiewicz, R.L. Taylor and J.M. Too, Reduced inte~ation technique in general analysis of plates and shell, Int. J. Num. Meth. Eng. 3 (1971) 275-290. [ 181 P.K. Larsen, Large displacement analysis of shells of revolution including creep, plasticity and viscoelasticity, Ph.D. Dissertation, University of California, Berkeley, Report No. UCSESM 71-22, December, 1971. [ 191 H. Weinitschke, On the nonlinear theory of shallow spherical shells, SIAM 6 (1958) 209-232. [ 2O] Modern Plastic Encyclopedia (McGraw-Hill, 1970). [ 21f J.H. Argyris and D.W. Scharf, “The SHEBA of shell elements for the matrix displacement method, Parts 1 and 1 I ,” The Aeronautical Journal of the Royal Aeronautical Society 72 (1968) 873-883. [22] J.H. Argyris and D.W. Scharpf, “The SHEBA family of shell elements for the matrix displacement methods, Part III: Large displacements,” The Aeronautical Jounrla of the Royal Aeronautical Society 73 (1969) 423-426. [ 231 A.S.L. Chan and A. Firmin, “The analysis of cooling towers by the matrix displacement finite element method, Part I: Small displacement,” The Aeronautical Journal of the Royal Aronautical Society 74 (1970) 826-835. [ 241 A.S.L. Chan and A. Firmin, “Large displacement analysis of axisymmetric shells, Part II: Formulation of the geometric stiffness matrix.” The Aeronautical Journal of the Royal Aeronautical Society 74 (1970) 971-982. [ 251 J.H. Argyris and D.W. Scharpf, “Matrix displacement analysis of shells and plates including transverse shear strain effects,” Computer Methods in Applied Mechanics and Engineering 1 (1972) 81-139.