- Email: [email protected]
Isoparametric finite-element thermoelasto-plastic creep analysis of shells of revolution
Printed ELSEVIER
in Northern
0308-0161(95)00061-5
Ireland.
All rights reserved 030X-0l6l/96/$15.00
Isoparametric finite-element thermoelastoplastic creep analysis of shells of revolution Mechanical
Engineering
Dept.,
M. Shariyat & M. R. Eslami Amirkabir University of Technology, Tehran
Polytechnic.
No. 424, Tehran.
lrtrn
(Received 1 May 1995:accepted 10 June 1995). First the strain-displacementrelations of a third order isoparametric shell element are derived. Using the variational approach the stiffnessand force matricesof a shell under distributed and concentrated loads and temperature gradient on distributed or concentrated elastic foundations are presented. Considering the factors that influence the material properties and the yield function, the thermoelasto-plasticand thermoelasto-plasticcreep constitutive equationsare introduced. Finally, a solution algorithm is proposedand applied to investigate the nature and variations of stressesproduced in cylindrical shellswith sphericaland torospherical heads. Copyright 0 1996Elsevier ScienceLtd.
1 INTRODUCTION
isoparametric shell element is widely considered to be the most efficient and economic.’ In some applications of multishell structures or composite shells, the inner and outer facings can be best approximated by a shell on elastic foundations. The elastic foundations can generally impose distributed or concentrated forces or moments on the shell. Using the variational approach,4*’ the stiffness and force matrices of a shell under distributed and concentrated loads and temperature gradient supported by elastic foundations, are derived. Considering isotropic strain-hardening behaviour for the material and taking into account all of the factors that affect the material properties and the yield function (particularly in situations which involve drastic temperature and strain rate changes), the thermoelasto-plastic and creep constitutive equations corresponding to the present analysis are proposed. Finally, an efficient solution procedure based on the employed consititutive equation is proposed. The method is applied to study the nature and behaviour of stresses produced in two commonly used types of cylindrical shells, one with a spherical and the other with a torospherical head. The main advantage of the solution procedure presented in this paper is that no iteration procedure is necessary during each time increment. This is in contrast to the initial stress, initial strain, and variable stiffness methods which are appropriate for the materially-nonlinear
Many engineering structures are vulnerable to inelastic deformations especially at elevated temperatures. A fundamental observation comparing elastic and inelastic analysis is that in elastic solutions, the total stress is evaluated from the tofal strain alone, whereas in an inelastic response calculation the total stress at time c depends on the stress and strain history. Typical inelastic phenomena are elasto-plasticity, creep and viscoelasticity. Considering the incremental analysis of inelastic response, there are basically three kinematical formulations that are effectively employed: materially-nonlinear analysis, total Lagrangian formulations (TL) and updated Lagrangian Jumann stress rate (ULJ) or updated Lagrangian Hencky (ULH) formulations.’ For small displacements and rotations, the results obtained using each of these approaches are approximately the same.’ Almost all of the solution procedures proposed based on the mentioned methods, suffer from the fact that they require an iterative procedure in each time increment. In the present paper, considering the continuity of the nodal variables and their first derivatives between the adjacent elements, a shell element is third-order isoparametric adopted. It is found that of all elements developed for plate and shell structures the 249
250
M. Sharivat, M. R. Eslami
analysis. It is also in contrast to the methods employed for the linearized TL, UJ and ULH formulations 2 ELASTIC
FORMULATION
To establish the stiffness and force matrices, a third-order isoparametric element as shown in Fig. 1 is used. The use of a dimensionless natural coordinate 5, which is defined as: 5 = 2 (s/L -O-S), 1‘n p1ace of the meridional coordinate s, is preferred. Now, if the vertical and horizontal displacements and their first derivatives are required to satisfy the continuity conditions between any adjacent elements, they should have the following form:
where: N(=$(E”-3&+2); N;=$(-t3+3t+2) NY= a([’ - 5’ - 5 + 1);
is the rotation matrix. For an isoparametric shell element, the variations of the element geometry and element displacements are assumed to be similar. So that, one can write: r = g [N,‘r, + NY($) 3 z = 2 [N:z; + NY ;-,
] I
(5)
From the element geometry, we have:
Thus, from the above equation and definition of 5, it can be easily verified that: (2)
N;=i(e3+t2-[-l)
d+ z= d2+ -= d(’
and:
{S}’ = (IT;.wj, $),, zj,N;, I I
where:
w,, 2) (5)) I I
(7) ~
-__
~
The principal curvature of the element in the meridional cross-section can be calculated as:
The meridional and lateral displacements in local coordinates are related to the vertical and horizontal displacements in global coordinates through the following equation: -dr
d’z __-_
dz, dc$‘d[’ d[’ The matrix relationship between displacements is:’
1
-d2r -’
(8) dtJ2 strains and
{~}*=[;]=[-!jg~;L.J (9 Therefore, after substituting the foregoing equations into the above equation, the matrix relationship between strains and nodal displacements reduces to the following: Fig. 1.
1~1”= [WI
(10)
Creep analysis of shells of revolution
The matrix [B] is a 4 X 8 matrix defined as:
The stress-strain relationship for thermoelastic deformation of the shell can be expressed as: {o-l=
where:” (B,)=[-sin~$/~+$cos@(~+~~~
31
{ ;j
= [D]k’)
= [D](k)
- (~~17‘)
(13)
in which {c} and {E’} are the total strain and the mechanical strain vectors, respectively. The total strain vector can be redefined as:
rer={ E‘+z.xj={E,,,)+~{x} EH + z . XH
X[N]+[cos$sin$1$[N].E
(14)
where {E,,,} is the strain vector of the middle surface of the shell and (x} is the curvature vector. Assuming plane-stress condition t’or thin shells, the elasticity matrix is:
w = [0
(15) Therefore according to eqns (9) and (IO), eqn (13) can be rewritten in the following form:
b> = PlW@~ - {a)T)
(16)
where:
(B,) + L. (B,) PI = [ @h)+ ,7.(I?,)1
(17)
The total potential energy of the shell element under mechanical and thermal loads and the elastic foundation may be obtained from the following general equation:J,i
+ [sin C#J-cos ~$1; [N] . $ + f [cos 4 sin 0][N]] , here:
[iv] =
N; 0 N$ 0
N(
0
0 N$
where the first term in eqn (18) is the strain energy, the second and third terms are the strain energies of the distributed and concentrated elastic foundations, and the fourth and fifth terms are the potential energies of the applied distributed and ring concentrated lateral forces acting on the shell. It is important to notice that while the unknown vector (6) has four nodal values U, W, dii/dt, and dG/dt, the reaction of the elastic foundation can be assumed to be
M. Shariyat.
-.‘FJ .
proportional to 14, W, and p. Here, p, the cross-section rotation, is expressed as:
M. R. Eslami
where:
tP1
WI
[D*] =
$I
PI + [-sin Cpcos +] $. [N]) . (6)
= WW
1
I
(251
1
c rc,1 (19)
It is further assumed that the distributed and concentrated ring forces include only forces in normal and tangential directions and moment about meridian, that is:
According to eqns (3) and (19), it is noticed that:
and:
r/7 T(s, z) dz: I - r/2 II? 7, = T(s, z )z dz i---l/2 to =
(26)
where t is the shell thickness.
h> = vx%where:
3 THERMOELASTO-PLASTIC
[Z] = [‘qy, calling:
P
1000 1 r 0100
0 0
000 000
.]=l- - - - _ _ _ -1 (A)
(22)
and upon substitution from eqns (16), (19) and (22) into eqn (18), the finite-element equilibrium equation is obtained using the variational principal, i.e.,
WI-WI= E,I + UC.1
{de} = {de,} + {de,,.} + {de,;} + {de,-}
(24)
where {de,} is the strain increment due to elastic deformation; {de,,) is the strain increment due to temperature dependent material properties; {dr,,;} is the strain incre,ment due to strain rate dependent material properties; {der) is the thermal strain increment. Therefore according to the generalized Hooke’s law for elastic stresses and strains the total stress increments can be computed from eqn (27) as:
[K](e) = ZE{ 11, ($1 . ( [B]~‘]D*][B]
+ PI7 IK,,lF+. d5 + c,[Rl’ [K,I[Rl} j-’
-I
WW~~ dt
+ c cr[Rl’(FWJ) {&}
= XL j-’
I
[B17-{M,}r . d,$
In the derivation of the following equations it is assumed that the material is isotropic, exhibits the same properties in tension and compression and the material properties are not history dependent. Also, the material is considered to have an isotropic strain-hardening behaviour. The plastic analysis is usually accomplished by incremental methods. For situations involving drastic temperature and strain rate changes, the total incremental strain components may be assumed to consist of the following components:”
(23)
where for the base element (e):
{&,I = 24
ANALYSIS
(21)
kW = [WW = [DIWI - k&I - kkd - kk,;) - id+)
(27)
(28)
The plastic potential function F is generally a function of variation of temperature, strain rate and amount of the plastic work: F = F({o},
K, T, Z)
(29)
253
Creep analysis of shells of revolution
where K is the work-hardening parameter of the material. For material obeying the Von Mises yield criterion under a multiaxial stress state the function F has the form: F = Jz - a: = :(Gz - a;)
(30)
Equation (34) takes the following form: dA=;({g}[D]({dc}-(a)d7’-F X{u}dT-=$JdC)
where a? is the yielding stress: (T is the effective stress and in tensor form:
(37) The plasticity and elasto-plasticity matrices are defined as:
For the present case the effective stress and strain are: a=-\/ 0-f + a’e - U,,U% (31) E=+&f+E;Es6% For Von Mises-type materials the following relation holds for the deviatoric stress vector: {a’} = {aF/au}
(32)
The plastic strain in eqn (27) can be expressed in terms of the Prandtl-Reuss relation as: {de,,} = dh . {a’} = dh{aF/au}
PI,, = $ [W4b’YPl FL = PI - PIP
(39) After substituting eqn (37) into eqn (33) and substituting the resulted equation and eqns (36), (38) and (39) into eqn (28), the thermoelastoplastic constitutive equation is obtained: 1 {da} = [D],{de} - [D]({a}T + s
(33)
When plastic yield is occurring the stresses are on the yield surface: F = 0 and hence dF is zero. Establishing dF = 0 from eqn (29) and substituting eqns (28) and (33) into the resulting equation and noting that K is a function of the plastic strains, the following expression is obtained for the proportionality factor dh:
-___
(40)
One may find that [D] = 2G during the plastic deformation of a solid, due to the fact that the first deviatoric strain invariant vanishes.h On the other hand according to the definition of the work hardening factor K: {dK
‘[D,({dr} dh =
- {de,.} - {d+}
{de,,;}) + %dT
+ gdt
(34) Recognizing the following relationships:6 {d+} = {a}T;
{d+}
=F
{a}
(38)
/de,,}’
=
(41) is easily
{U}’
and eqn (30), the following relation obtained: aF 2 dF 2 -=---a5 A=-H’ (42) 3 dK 3”’ t$,,’ dK where H’ is the plastic modulus of the material. Then, substituting eqns (32), (41) and (42) into eqns (36) and (38) leads to the following: S = $a2(3G + H’) (43) IT 9G2 u:u; (44) ‘D1p = (3G + H’)S ;;:A
1
dT;
Determination
{de<,;} = “[D_]m’ {a} d: ai (35) and letting: (36)
of
the
plastic
rnodrlius
According to the effective stresstw effective strain curve of the materials the effective strain can be divided into two components: dF = dC<,+ dC,, (45) Based on the tangent modulus definition,
M. Sharivat, M. R. Eslami
254
E, = dG/dZ. relation:
eqn
(45) leads to the following
To determine E,, we introduce following constitutive equation:
and use the
2 = $ [a + a,,(a - q,,)“‘j
E a,, . m,,(G - cTJ”q
4. FINITE-ELEMENT FORMULATION THERMOELASTIC-PLASTIC CREEP ANALYSIS
OF
(47)
This equation has a relatively simple and accurate form, specially for steels. Therefore, I?, can be expressed as:
ET=[l+
where k, and fz are material constants, Q is the activation energy, R is the Boltzman’s constant and T is the absolute temperature.
(48)
3.1 Coustitutive equations for themoelastoplastic creep analysis The total incremental strain components can be expressed as: {W = WcJ + k-k) (49) where {de,} represent the thermoelasto-plastic strain components given in eqn (27) and {de<} are the incremental creep strains. According to the Prandtl-Reuss relation for creep deformation of materials obeying the Von Mises yield criterion, the incremental creep strains are expressed as:
(50) The effective creep strain rate is defined as: gc = (;{~“}‘{~“})I/’ (51) since {de} in eqn (40) is equivalent to {de+,} in the present analysis [eqn (49)], eqn (40) can be written in the following finite incremental form:
During small load increments, it can be assumed that a piecewise linear relation between the stress and strain exists. Thus, if all physical quantities including the displacements, strains, etc. are dealt with in the incremental sense, the elastic formulation presented earlier can still be used for the thermoelasto-plastic deformation induced by small load increments. For this purpose, [D] and (6) should be substituted respectively by [D],,, and {ha}. Therefore, eqn (23) changes to the following form: (54) [Kl,,Wl = WitI + {A&l The incremental strain energy can be computed from: AU = f /,{AL}‘{Ag} dV (55) I Thus, the additional incremental load matrix induced by the creep effect can be shown to take the form: {AF,.} = 1 [B’]“[D],,{A~, >dV C’ in other words, the general form equilibrium equation is:
of
the
[K],,,Wl
= Wd + WI > + @r;: > (57) form of the which is the incremental thermoelasto-plastic creep state of equilibrium condition. 5 SOLUTION
For most engineering applications the general Norton’s law for multidimensional analysis is used. This equation can be expressed in the form:’ E = k,G” exp(-QlRT) (53)
(56)
PROCEDURE
The following solution algorithm is applicable for both thermoelasto-plastic and thermoelastoplastic creep analysis. The main advantage of the method presented here is that in contrast to the intitial stress, initial strain and variable stiffness approaches, no iterative procedure (in each time increment) is necessary. However, to ensure the convergence of the solution, very small load increments must be used.
Creep analysis of shells of revolution
In the following algorithm, it is assumed that the overall equilibrium equation is represented by [Kl,Wl = @RI1. Initialize {E} and {k} for all elements. 2. Select thermal and mechanical load increments. 3_. Perform thermal analysis and calculate nodal temperatures {T} at time t. 4. Select material properties based on the average temperature of the element and the strain rate determined from the previous step of loading. 5. Calculate the effective stress corresponding to the mid-point of the element or predict it using an interpolation method such as Lagrange interpolation method based on the nodal point stresses of the present and adjacent elements. Then determine H’ from eqns (46) and (48). 6. Establish matrix [B] from eqn (11). 7. Form the elastoplastic matrix [D], from eqns (39) and (44). 8. For stability of the results in the thermoelasto-plastic creep anlaysis, choose the time increment to be: At,,, = 4( 1 + Y)/3Ek,n&“-‘)
exp( - Q/RT)
9. After the full loads are applied, calculate {AE,.} from eqns (50) and (51), otherwise set it to zero. 10. Form the stiffness and force matrices of the element using eqns (57), (56) and (24). 11. Assemble the overall stiffness and force matrices, [K], and {AR}. 12. Modify {AhR} and [K], corresponding to the applied boundary conditions. 13. After solving the resulting system of equations, {Aa} is found. Then, the corresponding displacement and strain increment is obtained by: {aI= iSI+ -WI
255
Spherical \ Toroidal
172.0’ \
-++-h
Fig. 2.
18. In the case of creep analysis, proceed to the next time increment and repeat Steps 4-16 until the final time is reached. 6 APPLICATIONS
AND RESULTS
The foregoing method was applied to investigate and compare the stress variations in two cylindrical shells with different heads, subjected to mechanical and thermal loadings. The dimensions of the compared shells and their boundary conditions are shown in Figs 2 and 3. These forms are commonly used in large tanks and pressure vessels. As may be noticed from these figures, the length of the cylindrical portion and the boundary conditions of the shells are the same. The shells were considered to be made of a material whose properties are: E = 22.4 X 10hpsi; Y = 0.3: CT,.= 28000 psi: m,, = 4: a = 10-h; % = 1.31479 x lo-“; n,, = 3; k, = 2.6568 X lo-’
{Ad = PIW 14. Compute {Aa} from eqn (52). 15. Compute:
16. Update nodal coordinates by letting: r=r+G;
z=z+ii
17. Repeat the foregoing steps from Step 2, until the full loads are applied.
t 24’
I
I
h
c
Fig. 3.
M. Shariyat,
256
M. R. Eslami
Effective Torospherical,
Circumferential Stress Variations Torospherical, Elastoplast . & Creep Circumferential
Stress
(Psi)
5. Effective /__._.__..
(Thousands)
Strerr
Stress Variations Elastoplast. 8 Creep
(Psi)
(Thousands)
------
307-i
f-
-40
’ 0
I
I
I
I
I
I
I
I
I
20
40
60
60
100
120
140
160
Global
Meridional
Coordinate
-
Ehatlc
-
Elaatopla81
- --
cremp1
- -
creep2
40
i
I
0’
(in)
20
0
40
Global -
Elutlc
I
I
60
60
Meridional -
I
/
120
140
I
100
Coordinate -
Ehrtophrt
Creep1
160
(in) Croop2
Fig. 4.
Fig. 6.
Figures 4 to 15 illustrate the stress variations at the middle surface. The final internal pressure and temperature in Figs 4-9 are 40 psi and 100°F and in Figs lo-l.5 are 60 psi and 100°F
respectively. The elastic curves in Figs 4-9 and 13-15 correspond to application of one half of the loads. Investigation of Figs 4 and 7 shows that due to presence of local moments, abrupt
Meridional Torospherical, 26 Meridional
Stress
Stress Variations Elastoplast . 8, Creep (Pri)
Circumferential Stress Variations Spherical Cap, Elastuplastic h Creep Circumferential
(Thousands)
Stre88
(P8i)
(Thousandr)
1 25.
-15’ 0
I
1
I
I
,
I
I
20
40
60
60
100
120
140
Global
Meridional
Coordinate
(in)
I 160
0’ 0
I
I
50
100
Global -.
ElmtIc _ Cmop1
-
200
150
Coordinate
(in)
Elartoplast croop2
Fig. 5.
Meridional
-
Elwtlo
-
Ehtopkat
- --crewp1
Fig. 7.
-crap2
Creep analysis of shells of revolution
Circumferential (Torospherical
Meridional Stress Variations Spherical Cap, Elastoplastic 8, Creep Meridional
Stress
(Psi)
251
Circumferential
(Thousands)
40
121
Stress Variations Cap , lncr. Press. )
Streae
(Psi)
(Thousands)
r
20
6 -60
1
-60 -100 ,
/
,
50
100
150
0’ 0
Global -
Eta~tlc
Meridional -
-120
Coordinate - -
Elmoplmt
t
200
(in)
cmep1
-
-140
,i
I
0
20
I
/
40
60
Global
crmp2
I
80
Meridional
100
I
120
Coordinate
140
180
(in)
Fig. 8.
Fig. 10.
changes in the circumferential stress have occurred at the junction points (particularly at the spherical-toroidal junction). The elastic results of the torospherical head are
in good agreement with the experimental results presented in Ref. 8. It is seen that the stress concentration in the spherical head is markedly less than that of the torospherical head. Study of
Meridional Stress Variations (Torospherical Cap , Incr. Press. )
Effective Stress Variations Spherical Cap, Elastoplastic h Creep
Meredional Effective
Streso
(P8i)
Strew
(Pail
(Thousands)
150
(Thousands)
20; 100
-
50 Pi = 60 [psi] ==%I z=24 -12
0
-50
0’ 0
I
I
I
I
50
100
150
200
Global -
ElBBtlc
Meridional -
Coordinate --
Ehrtople*t
Fig. 9.
- Craep1
(in) - --
crrrp2
-
-loo0
20
Global
40
80
80
Meridional Fig. 11.
100
120
Coordinate
140
(in)
160
M. Shariyat, M. R. Eslami
258
Effective Stress Variations (Torospherical Cap , Incr. Press. ) Effective 200
r
Stress
(Psi)
Circumferential Stress Variations Spherical Cap, Elastoplastic 8 Creep
(Thousands)
35 Circumferential (--.
Stress
(Psi)
(Thousands) _--__--
__._ .-
30
150
f7-1
2.5
I
/ /
20’ 1
/ ,.
15
-4
t
50
Global
A
20
Global
40
60
80
Meridional
100
120
Coordinate
140
160
-
(in)
‘J :
.
0
0
-~~
Ela8tic
100
Meridional -
/--
150
200
Coordinate -
Elaatopla*t
Creep1
(in) -
creep2
Fig. 12.
Fig. 13.
the effective stress variations of these shells, Figs 6 and 9, shows that while some regions of the torospherical shell undergo plastic deformations, the other shell is elastic everywhere. Figures lo-15 illustrate the results for increased
internal pressure. As can be noticed, failure takes place in the spherical-toroidal junction of the torospherical shell due to excessive plastic deformations, whereas the shell with a spherical head is still elastic.
Meridional
-
Elastic
-
Global Elastoplaat
Stress
Meridional
Variations
Cowdinote
(in) ----
Fig. 14.
creep1
----
creep2
Creep
of shells of revolution
analysis
Effective Stress Variations Spherical Cap, Elastoplastic 8 Creep Effective
Stress
(Psi)
(Thousands)
259
and secondary stress concepts for prediction of the inelastic behaviour of the shells are described in detail in Refs 9, 10 and are in agreement with the present results.
30,
25’ t
REFERENCES
20 1
5t 0’
1 0
/
50
Global Meridional -
Elartic
-
I
1
150
100
200
Coordinate
(in)
Creep1
- -- Creep2
Elartoplaat
Fig. 15.
The thermoelasto-plastic creep curves, are designated by Creep1 and Creep2. The Creep1 curves correspond to t = 500 hr and t = 1000 hr and Creep2 curves correspond to t = 1000 hr and t = 2000 hr, for the shell with torospherical and spherical head, respectively. Investigation of the creep curves demonstrates that as time elapses, the stress variations in the neighbourhood of the geometrical changes become smoother. The few exceptions observed are due to the fact that in the finite-element method it is assumed that the variables do not change abruptly and are due to presence of some primary stresses. The primary
1. Bathe, K. J., Finite elementelastic-plasticanalysisusing updated Lagrangian-Hencky method, in: Finite Element Methods for Nonlinear Problems. edited by Bergman,P. G., Bathe. K. J., Wundelich, W., Springer-Verlag (1986). 2. Bathe, K. J.. Finite Element Procedures in Engineering Analysis, Prentice-Hall (1982) 386-396. 3. Zienkiewicz, 0. C., The Finite Element Method in Engineering Science, McGraw-Hill (1983). 4. Eslami, M. R. & Shariyat, M., Isoparametric finite element of shell of revolution on elastic foundation, Proc. ASME-PVP Conf. (1991) Melbourne-Astre. 5. Eslami, M. E. & Shariyat. M., A third-order isoparametricfinite element of shell of revolution, to be publishedin Int. J. Etzg. Anal. Des. 6. Hsu, T. R., Finite Element Methods in Thermomomechanics, Prentice-Hall, Inc. (1986). 7. Penny, R. K. & Marriot, D. L., Design for Creep, McGraw-Hill (1971). 8. Gouid, P. L., Fin&e Element Analysis of Shell of Revolution, Pitman Adv. Pub. (1988)73. 9. Eslami, M. R. & Shariyat, M., The primary and secondary stressesin pressure vessels,Proc. ASMEPVP Co@, San Diego (June 23-27, 1991). 10. Eslami, M. R. & Shariyat, M., A technique to distinguish the primary and secondary stresses,to be publishedin Trans. ASME J. Press. Vess. 11. Eslami, M. R., Ahmadian, H. & Ayatollahi, M. R., Static and dynamic finite element solution of plate and shell of revolution on elastic foundation, Proc. ASME-PVP Conf. (June 17-21, 1990)Nashville. 12. Eslami, M. R. & Alizadeh, S. H., Galerkin finite element formulation of spherical shells under nonaxisymmetric loading, Proc. ASME-PVP Corzf. (June 17-21, 1990)Nashville. 13. Huang, H. C., Static and Dynamic Analysis of Plates and Shells-Theory, Software and Applications. SpringerVerlag, Berlin (1989).
in Northern
0308-0161(95)00061-5
Ireland.
All rights reserved 030X-0l6l/96/$15.00
Isoparametric finite-element thermoelastoplastic creep analysis of shells of revolution Mechanical
Engineering
Dept.,
M. Shariyat & M. R. Eslami Amirkabir University of Technology, Tehran
Polytechnic.
No. 424, Tehran.
lrtrn
(Received 1 May 1995:accepted 10 June 1995). First the strain-displacementrelations of a third order isoparametric shell element are derived. Using the variational approach the stiffnessand force matricesof a shell under distributed and concentrated loads and temperature gradient on distributed or concentrated elastic foundations are presented. Considering the factors that influence the material properties and the yield function, the thermoelasto-plasticand thermoelasto-plasticcreep constitutive equationsare introduced. Finally, a solution algorithm is proposedand applied to investigate the nature and variations of stressesproduced in cylindrical shellswith sphericaland torospherical heads. Copyright 0 1996Elsevier ScienceLtd.
1 INTRODUCTION
isoparametric shell element is widely considered to be the most efficient and economic.’ In some applications of multishell structures or composite shells, the inner and outer facings can be best approximated by a shell on elastic foundations. The elastic foundations can generally impose distributed or concentrated forces or moments on the shell. Using the variational approach,4*’ the stiffness and force matrices of a shell under distributed and concentrated loads and temperature gradient supported by elastic foundations, are derived. Considering isotropic strain-hardening behaviour for the material and taking into account all of the factors that affect the material properties and the yield function (particularly in situations which involve drastic temperature and strain rate changes), the thermoelasto-plastic and creep constitutive equations corresponding to the present analysis are proposed. Finally, an efficient solution procedure based on the employed consititutive equation is proposed. The method is applied to study the nature and behaviour of stresses produced in two commonly used types of cylindrical shells, one with a spherical and the other with a torospherical head. The main advantage of the solution procedure presented in this paper is that no iteration procedure is necessary during each time increment. This is in contrast to the initial stress, initial strain, and variable stiffness methods which are appropriate for the materially-nonlinear
Many engineering structures are vulnerable to inelastic deformations especially at elevated temperatures. A fundamental observation comparing elastic and inelastic analysis is that in elastic solutions, the total stress is evaluated from the tofal strain alone, whereas in an inelastic response calculation the total stress at time c depends on the stress and strain history. Typical inelastic phenomena are elasto-plasticity, creep and viscoelasticity. Considering the incremental analysis of inelastic response, there are basically three kinematical formulations that are effectively employed: materially-nonlinear analysis, total Lagrangian formulations (TL) and updated Lagrangian Jumann stress rate (ULJ) or updated Lagrangian Hencky (ULH) formulations.’ For small displacements and rotations, the results obtained using each of these approaches are approximately the same.’ Almost all of the solution procedures proposed based on the mentioned methods, suffer from the fact that they require an iterative procedure in each time increment. In the present paper, considering the continuity of the nodal variables and their first derivatives between the adjacent elements, a shell element is third-order isoparametric adopted. It is found that of all elements developed for plate and shell structures the 249
250
M. Sharivat, M. R. Eslami
analysis. It is also in contrast to the methods employed for the linearized TL, UJ and ULH formulations 2 ELASTIC
FORMULATION
To establish the stiffness and force matrices, a third-order isoparametric element as shown in Fig. 1 is used. The use of a dimensionless natural coordinate 5, which is defined as: 5 = 2 (s/L -O-S), 1‘n p1ace of the meridional coordinate s, is preferred. Now, if the vertical and horizontal displacements and their first derivatives are required to satisfy the continuity conditions between any adjacent elements, they should have the following form:
where: N(=$(E”-3&+2); N;=$(-t3+3t+2) NY= a([’ - 5’ - 5 + 1);
is the rotation matrix. For an isoparametric shell element, the variations of the element geometry and element displacements are assumed to be similar. So that, one can write: r = g [N,‘r, + NY($) 3 z = 2 [N:z; + NY ;-,
] I
(5)
From the element geometry, we have:
Thus, from the above equation and definition of 5, it can be easily verified that: (2)
N;=i(e3+t2-[-l)
d+ z= d2+ -= d(’
and:
{S}’ = (IT;.wj, $),, zj,N;, I I
where:
w,, 2) (5)) I I
(7) ~
-__
~
The principal curvature of the element in the meridional cross-section can be calculated as:
The meridional and lateral displacements in local coordinates are related to the vertical and horizontal displacements in global coordinates through the following equation: -dr
d’z __-_
dz, dc$‘d[’ d[’ The matrix relationship between displacements is:’
1
-d2r -’
(8) dtJ2 strains and
{~}*=[;]=[-!jg~;L.J (9 Therefore, after substituting the foregoing equations into the above equation, the matrix relationship between strains and nodal displacements reduces to the following: Fig. 1.
1~1”= [WI
(10)
Creep analysis of shells of revolution
The matrix [B] is a 4 X 8 matrix defined as:
The stress-strain relationship for thermoelastic deformation of the shell can be expressed as: {o-l=
where:” (B,)=[-sin~$/~+$cos@(~+~~~
31
{ ;j
= [D]k’)
= [D](k)
- (~~17‘)
(13)
in which {c} and {E’} are the total strain and the mechanical strain vectors, respectively. The total strain vector can be redefined as:
rer={ E‘+z.xj={E,,,)+~{x} EH + z . XH
X[N]+[cos$sin$1$[N].E
(14)
where {E,,,} is the strain vector of the middle surface of the shell and (x} is the curvature vector. Assuming plane-stress condition t’or thin shells, the elasticity matrix is:
w = [0
(15) Therefore according to eqns (9) and (IO), eqn (13) can be rewritten in the following form:
b> = PlW@~ - {a)T)
(16)
where:
(B,) + L. (B,) PI = [ @h)+ ,7.(I?,)1
(17)
The total potential energy of the shell element under mechanical and thermal loads and the elastic foundation may be obtained from the following general equation:J,i
+ [sin C#J-cos ~$1; [N] . $ + f [cos 4 sin 0][N]] , here:
[iv] =
N; 0 N$ 0
N(
0
0 N$
where the first term in eqn (18) is the strain energy, the second and third terms are the strain energies of the distributed and concentrated elastic foundations, and the fourth and fifth terms are the potential energies of the applied distributed and ring concentrated lateral forces acting on the shell. It is important to notice that while the unknown vector (6) has four nodal values U, W, dii/dt, and dG/dt, the reaction of the elastic foundation can be assumed to be
M. Shariyat.
-.‘FJ .
proportional to 14, W, and p. Here, p, the cross-section rotation, is expressed as:
M. R. Eslami
where:
tP1
WI
[D*] =
$I
PI + [-sin Cpcos +] $. [N]) . (6)
= WW
1
I
(251
1
c rc,1 (19)
It is further assumed that the distributed and concentrated ring forces include only forces in normal and tangential directions and moment about meridian, that is:
According to eqns (3) and (19), it is noticed that:
and:
r/7 T(s, z) dz: I - r/2 II? 7, = T(s, z )z dz i---l/2 to =
(26)
where t is the shell thickness.
h> = vx%where:
3 THERMOELASTO-PLASTIC
[Z] = [‘qy, calling:
P
1000 1 r 0100
0 0
000 000
.]=l- - - - _ _ _ -1 (A)
(22)
and upon substitution from eqns (16), (19) and (22) into eqn (18), the finite-element equilibrium equation is obtained using the variational principal, i.e.,
WI-WI= E,I + UC.1
{de} = {de,} + {de,,.} + {de,;} + {de,-}
(24)
where {de,} is the strain increment due to elastic deformation; {de,,) is the strain increment due to temperature dependent material properties; {dr,,;} is the strain incre,ment due to strain rate dependent material properties; {der) is the thermal strain increment. Therefore according to the generalized Hooke’s law for elastic stresses and strains the total stress increments can be computed from eqn (27) as:
[K](e) = ZE{ 11, ($1 . ( [B]~‘]D*][B]
+ PI7 IK,,lF+. d5 + c,[Rl’ [K,I[Rl} j-’
-I
WW~~ dt
+ c cr[Rl’(FWJ) {&}
= XL j-’
I
[B17-{M,}r . d,$
In the derivation of the following equations it is assumed that the material is isotropic, exhibits the same properties in tension and compression and the material properties are not history dependent. Also, the material is considered to have an isotropic strain-hardening behaviour. The plastic analysis is usually accomplished by incremental methods. For situations involving drastic temperature and strain rate changes, the total incremental strain components may be assumed to consist of the following components:”
(23)
where for the base element (e):
{&,I = 24
ANALYSIS
(21)
kW = [WW = [DIWI - k&I - kkd - kk,;) - id+)
(27)
(28)
The plastic potential function F is generally a function of variation of temperature, strain rate and amount of the plastic work: F = F({o},
K, T, Z)
(29)
253
Creep analysis of shells of revolution
where K is the work-hardening parameter of the material. For material obeying the Von Mises yield criterion under a multiaxial stress state the function F has the form: F = Jz - a: = :(Gz - a;)
(30)
Equation (34) takes the following form: dA=;({g}[D]({dc}-(a)d7’-F X{u}dT-=$JdC)
where a? is the yielding stress: (T is the effective stress and in tensor form:
(37) The plasticity and elasto-plasticity matrices are defined as:
For the present case the effective stress and strain are: a=-\/ 0-f + a’e - U,,U% (31) E=+&f+E;Es6% For Von Mises-type materials the following relation holds for the deviatoric stress vector: {a’} = {aF/au}
(32)
The plastic strain in eqn (27) can be expressed in terms of the Prandtl-Reuss relation as: {de,,} = dh . {a’} = dh{aF/au}
PI,, = $ [W4b’YPl FL = PI - PIP
(39) After substituting eqn (37) into eqn (33) and substituting the resulted equation and eqns (36), (38) and (39) into eqn (28), the thermoelastoplastic constitutive equation is obtained: 1 {da} = [D],{de} - [D]({a}T + s
(33)
When plastic yield is occurring the stresses are on the yield surface: F = 0 and hence dF is zero. Establishing dF = 0 from eqn (29) and substituting eqns (28) and (33) into the resulting equation and noting that K is a function of the plastic strains, the following expression is obtained for the proportionality factor dh:
-___
(40)
One may find that [D] = 2G during the plastic deformation of a solid, due to the fact that the first deviatoric strain invariant vanishes.h On the other hand according to the definition of the work hardening factor K: {dK
‘[D,({dr} dh =
- {de,.} - {d+}
{de,,;}) + %dT
+ gdt
(34) Recognizing the following relationships:6 {d+} = {a}T;
{d+}
=F
{a}
(38)
/de,,}’
=
(41) is easily
{U}’
and eqn (30), the following relation obtained: aF 2 dF 2 -=---a5 A=-H’ (42) 3 dK 3”’ t$,,’ dK where H’ is the plastic modulus of the material. Then, substituting eqns (32), (41) and (42) into eqns (36) and (38) leads to the following: S = $a2(3G + H’) (43) IT 9G2 u:u; (44) ‘D1p = (3G + H’)S ;;:A
1
dT;
Determination
{de<,;} = “[D_]m’ {a} d: ai (35) and letting: (36)
of
the
plastic
rnodrlius
According to the effective stresstw effective strain curve of the materials the effective strain can be divided into two components: dF = dC<,+ dC,, (45) Based on the tangent modulus definition,
M. Sharivat, M. R. Eslami
254
E, = dG/dZ. relation:
eqn
(45) leads to the following
To determine E,, we introduce following constitutive equation:
and use the
2 = $ [a + a,,(a - q,,)“‘j
E a,, . m,,(G - cTJ”q
4. FINITE-ELEMENT FORMULATION THERMOELASTIC-PLASTIC CREEP ANALYSIS
OF
(47)
This equation has a relatively simple and accurate form, specially for steels. Therefore, I?, can be expressed as:
ET=[l+
where k, and fz are material constants, Q is the activation energy, R is the Boltzman’s constant and T is the absolute temperature.
(48)
3.1 Coustitutive equations for themoelastoplastic creep analysis The total incremental strain components can be expressed as: {W = WcJ + k-k) (49) where {de,} represent the thermoelasto-plastic strain components given in eqn (27) and {de<} are the incremental creep strains. According to the Prandtl-Reuss relation for creep deformation of materials obeying the Von Mises yield criterion, the incremental creep strains are expressed as:
(50) The effective creep strain rate is defined as: gc = (;{~“}‘{~“})I/’ (51) since {de} in eqn (40) is equivalent to {de+,} in the present analysis [eqn (49)], eqn (40) can be written in the following finite incremental form:
During small load increments, it can be assumed that a piecewise linear relation between the stress and strain exists. Thus, if all physical quantities including the displacements, strains, etc. are dealt with in the incremental sense, the elastic formulation presented earlier can still be used for the thermoelasto-plastic deformation induced by small load increments. For this purpose, [D] and (6) should be substituted respectively by [D],,, and {ha}. Therefore, eqn (23) changes to the following form: (54) [Kl,,Wl = WitI + {A&l The incremental strain energy can be computed from: AU = f /,{AL}‘{Ag} dV (55) I Thus, the additional incremental load matrix induced by the creep effect can be shown to take the form: {AF,.} = 1 [B’]“[D],,{A~, >dV C’ in other words, the general form equilibrium equation is:
of
the
[K],,,Wl
= Wd + WI > + @r;: > (57) form of the which is the incremental thermoelasto-plastic creep state of equilibrium condition. 5 SOLUTION
For most engineering applications the general Norton’s law for multidimensional analysis is used. This equation can be expressed in the form:’ E = k,G” exp(-QlRT) (53)
(56)
PROCEDURE
The following solution algorithm is applicable for both thermoelasto-plastic and thermoelastoplastic creep analysis. The main advantage of the method presented here is that in contrast to the intitial stress, initial strain and variable stiffness approaches, no iterative procedure (in each time increment) is necessary. However, to ensure the convergence of the solution, very small load increments must be used.
Creep analysis of shells of revolution
In the following algorithm, it is assumed that the overall equilibrium equation is represented by [Kl,Wl = @RI1. Initialize {E} and {k} for all elements. 2. Select thermal and mechanical load increments. 3_. Perform thermal analysis and calculate nodal temperatures {T} at time t. 4. Select material properties based on the average temperature of the element and the strain rate determined from the previous step of loading. 5. Calculate the effective stress corresponding to the mid-point of the element or predict it using an interpolation method such as Lagrange interpolation method based on the nodal point stresses of the present and adjacent elements. Then determine H’ from eqns (46) and (48). 6. Establish matrix [B] from eqn (11). 7. Form the elastoplastic matrix [D], from eqns (39) and (44). 8. For stability of the results in the thermoelasto-plastic creep anlaysis, choose the time increment to be: At,,, = 4( 1 + Y)/3Ek,n&“-‘)
exp( - Q/RT)
9. After the full loads are applied, calculate {AE,.} from eqns (50) and (51), otherwise set it to zero. 10. Form the stiffness and force matrices of the element using eqns (57), (56) and (24). 11. Assemble the overall stiffness and force matrices, [K], and {AR}. 12. Modify {AhR} and [K], corresponding to the applied boundary conditions. 13. After solving the resulting system of equations, {Aa} is found. Then, the corresponding displacement and strain increment is obtained by: {aI= iSI+ -WI
255
Spherical \ Toroidal
172.0’ \
-++-h
Fig. 2.
18. In the case of creep analysis, proceed to the next time increment and repeat Steps 4-16 until the final time is reached. 6 APPLICATIONS
AND RESULTS
The foregoing method was applied to investigate and compare the stress variations in two cylindrical shells with different heads, subjected to mechanical and thermal loadings. The dimensions of the compared shells and their boundary conditions are shown in Figs 2 and 3. These forms are commonly used in large tanks and pressure vessels. As may be noticed from these figures, the length of the cylindrical portion and the boundary conditions of the shells are the same. The shells were considered to be made of a material whose properties are: E = 22.4 X 10hpsi; Y = 0.3: CT,.= 28000 psi: m,, = 4: a = 10-h; % = 1.31479 x lo-“; n,, = 3; k, = 2.6568 X lo-’
{Ad = PIW 14. Compute {Aa} from eqn (52). 15. Compute:
16. Update nodal coordinates by letting: r=r+G;
z=z+ii
17. Repeat the foregoing steps from Step 2, until the full loads are applied.
t 24’
I
I
h
c
Fig. 3.
M. Shariyat,
256
M. R. Eslami
Effective Torospherical,
Circumferential Stress Variations Torospherical, Elastoplast . & Creep Circumferential
Stress
(Psi)
5. Effective /__._.__..
(Thousands)
Strerr
Stress Variations Elastoplast. 8 Creep
(Psi)
(Thousands)
------
307-i
f-
-40
’ 0
I
I
I
I
I
I
I
I
I
20
40
60
60
100
120
140
160
Global
Meridional
Coordinate
-
Ehatlc
-
Elaatopla81
- --
cremp1
- -
creep2
40
i
I
0’
(in)
20
0
40
Global -
Elutlc
I
I
60
60
Meridional -
I
/
120
140
I
100
Coordinate -
Ehrtophrt
Creep1
160
(in) Croop2
Fig. 4.
Fig. 6.
Figures 4 to 15 illustrate the stress variations at the middle surface. The final internal pressure and temperature in Figs 4-9 are 40 psi and 100°F and in Figs lo-l.5 are 60 psi and 100°F
respectively. The elastic curves in Figs 4-9 and 13-15 correspond to application of one half of the loads. Investigation of Figs 4 and 7 shows that due to presence of local moments, abrupt
Meridional Torospherical, 26 Meridional
Stress
Stress Variations Elastoplast . 8, Creep (Pri)
Circumferential Stress Variations Spherical Cap, Elastuplastic h Creep Circumferential
(Thousands)
Stre88
(P8i)
(Thousandr)
1 25.
-15’ 0
I
1
I
I
,
I
I
20
40
60
60
100
120
140
Global
Meridional
Coordinate
(in)
I 160
0’ 0
I
I
50
100
Global -.
ElmtIc _ Cmop1
-
200
150
Coordinate
(in)
Elartoplast croop2
Fig. 5.
Meridional
-
Elwtlo
-
Ehtopkat
- --crewp1
Fig. 7.
-crap2
Creep analysis of shells of revolution
Circumferential (Torospherical
Meridional Stress Variations Spherical Cap, Elastoplastic 8, Creep Meridional
Stress
(Psi)
251
Circumferential
(Thousands)
40
121
Stress Variations Cap , lncr. Press. )
Streae
(Psi)
(Thousands)
r
20
6 -60
1
-60 -100 ,
/
,
50
100
150
0’ 0
Global -
Eta~tlc
Meridional -
-120
Coordinate - -
Elmoplmt
t
200
(in)
cmep1
-
-140
,i
I
0
20
I
/
40
60
Global
crmp2
I
80
Meridional
100
I
120
Coordinate
140
180
(in)
Fig. 8.
Fig. 10.
changes in the circumferential stress have occurred at the junction points (particularly at the spherical-toroidal junction). The elastic results of the torospherical head are
in good agreement with the experimental results presented in Ref. 8. It is seen that the stress concentration in the spherical head is markedly less than that of the torospherical head. Study of
Meridional Stress Variations (Torospherical Cap , Incr. Press. )
Effective Stress Variations Spherical Cap, Elastoplastic h Creep
Meredional Effective
Streso
(P8i)
Strew
(Pail
(Thousands)
150
(Thousands)
20; 100
-
50 Pi = 60 [psi] ==%I z=24 -12
0
-50
0’ 0
I
I
I
I
50
100
150
200
Global -
ElBBtlc
Meridional -
Coordinate --
Ehrtople*t
Fig. 9.
- Craep1
(in) - --
crrrp2
-
-loo0
20
Global
40
80
80
Meridional Fig. 11.
100
120
Coordinate
140
(in)
160
M. Shariyat, M. R. Eslami
258
Effective Stress Variations (Torospherical Cap , Incr. Press. ) Effective 200
r
Stress
(Psi)
Circumferential Stress Variations Spherical Cap, Elastoplastic 8 Creep
(Thousands)
35 Circumferential (--.
Stress
(Psi)
(Thousands) _--__--
__._ .-
30
150
f7-1
2.5
I
/ /
20’ 1
/ ,.
15
-4
t
50
Global
A
20
Global
40
60
80
Meridional
100
120
Coordinate
140
160
-
(in)
‘J :
.
0
0
-~~
Ela8tic
100
Meridional -
/--
150
200
Coordinate -
Elaatopla*t
Creep1
(in) -
creep2
Fig. 12.
Fig. 13.
the effective stress variations of these shells, Figs 6 and 9, shows that while some regions of the torospherical shell undergo plastic deformations, the other shell is elastic everywhere. Figures lo-15 illustrate the results for increased
internal pressure. As can be noticed, failure takes place in the spherical-toroidal junction of the torospherical shell due to excessive plastic deformations, whereas the shell with a spherical head is still elastic.
Meridional
-
Elastic
-
Global Elastoplaat
Stress
Meridional
Variations
Cowdinote
(in) ----
Fig. 14.
creep1
----
creep2
Creep
of shells of revolution
analysis
Effective Stress Variations Spherical Cap, Elastoplastic 8 Creep Effective
Stress
(Psi)
(Thousands)
259
and secondary stress concepts for prediction of the inelastic behaviour of the shells are described in detail in Refs 9, 10 and are in agreement with the present results.
30,
25’ t
REFERENCES
20 1
5t 0’
1 0
/
50
Global Meridional -
Elartic
-
I
1
150
100
200
Coordinate
(in)
Creep1
- -- Creep2
Elartoplaat
Fig. 15.
The thermoelasto-plastic creep curves, are designated by Creep1 and Creep2. The Creep1 curves correspond to t = 500 hr and t = 1000 hr and Creep2 curves correspond to t = 1000 hr and t = 2000 hr, for the shell with torospherical and spherical head, respectively. Investigation of the creep curves demonstrates that as time elapses, the stress variations in the neighbourhood of the geometrical changes become smoother. The few exceptions observed are due to the fact that in the finite-element method it is assumed that the variables do not change abruptly and are due to presence of some primary stresses. The primary
1. Bathe, K. J., Finite elementelastic-plasticanalysisusing updated Lagrangian-Hencky method, in: Finite Element Methods for Nonlinear Problems. edited by Bergman,P. G., Bathe. K. J., Wundelich, W., Springer-Verlag (1986). 2. Bathe, K. J.. Finite Element Procedures in Engineering Analysis, Prentice-Hall (1982) 386-396. 3. Zienkiewicz, 0. C., The Finite Element Method in Engineering Science, McGraw-Hill (1983). 4. Eslami, M. R. & Shariyat, M., Isoparametric finite element of shell of revolution on elastic foundation, Proc. ASME-PVP Conf. (1991) Melbourne-Astre. 5. Eslami, M. E. & Shariyat. M., A third-order isoparametricfinite element of shell of revolution, to be publishedin Int. J. Etzg. Anal. Des. 6. Hsu, T. R., Finite Element Methods in Thermomomechanics, Prentice-Hall, Inc. (1986). 7. Penny, R. K. & Marriot, D. L., Design for Creep, McGraw-Hill (1971). 8. Gouid, P. L., Fin&e Element Analysis of Shell of Revolution, Pitman Adv. Pub. (1988)73. 9. Eslami, M. R. & Shariyat, M., The primary and secondary stressesin pressure vessels,Proc. ASMEPVP Co@, San Diego (June 23-27, 1991). 10. Eslami, M. R. & Shariyat, M., A technique to distinguish the primary and secondary stresses,to be publishedin Trans. ASME J. Press. Vess. 11. Eslami, M. R., Ahmadian, H. & Ayatollahi, M. R., Static and dynamic finite element solution of plate and shell of revolution on elastic foundation, Proc. ASME-PVP Conf. (June 17-21, 1990)Nashville. 12. Eslami, M. R. & Alizadeh, S. H., Galerkin finite element formulation of spherical shells under nonaxisymmetric loading, Proc. ASME-PVP Corzf. (June 17-21, 1990)Nashville. 13. Huang, H. C., Static and Dynamic Analysis of Plates and Shells-Theory, Software and Applications. SpringerVerlag, Berlin (1989).