Elasto-plastic analysis of axisymmetric structures subject to arbitrary loads by hybrid-stress finite elements

cQ45-7949/M s3.w + .al Pergamon Press Ltd.

Computers& StructuresVol. 19, No. 3, pp. 447-465, 1984 Printed in the U.S.A.

ELASTO-PLASTIC ANALYSIS OF AXISYMMETRIC STRUCTURES SUBJECT TO ARBITRARY LOADS BY HYBRID-STRESS FINITE ELEMENTS SOM P. SINGHt Applied Technology Division, Association of American Railroads, 3140 S. Federal, Chicago, IL 60616, U.S.A

and R. L. SPILKER~ Department of Civil Engineering, Mechanics and Metallurgy, University of Illinois at Chicago, IL 60680, U.S.A. (Received 26 July 1983; receivedfor

publication 26 September

1983)

Abstract-A semi-analytic finite element method in conjuction with a hybrid-stress functional based on the initial stress approach is presented for elasto-plastic analysis. A three-dimensional solid element based on hybrid-stress model is also implemented for the elasto-plastic analysis. A procedure to compute the non-linear effects in terms of Fourier series in the hybrid-stress model is described. The accuracy and

efficiencyof the semi-analytic method is evaluated via numerical examples by comparing the solution with a full 3-D solution. The semi-analytic method is observed to be a viable alternative to 3-D analysis in elasto-plastic

analysis of axisymmetric structures subject to arbitrary loads.

INTRODUCTION

The finite element method has become an accepted technique to perform linear and nonlinear analysis of structural problems involving complex loading and geometry. A special class of three dimensional structures is one for which geometry, material properties, boundary conditions and loading are independent of one direction; e.g. axisymmetric structures have stresses, strains and displacements which are independent of the circumferential direction, 6, when loading is independent of 6. In such cases a two dimensional (e.g. axisymmetric) analysis can be used. If an axisymmetric structure is subjected to arbitrary loading (i.e. dependent on O), one possible recourse is three dimensional (3-D) analysis. A computationally effective alternative to 3-D analysis in these cases is the semi-analytic finite element approach. In this approach, load, displacement, stresses and strains are expanded in a Fourier (harmonic) Series and the complete analysis reduces to the summation of a series of 2-D analyses. The accuracy and computational efficiency of the semi-analytic method in linear analysis is demonstrated in the works of Wilson[l], Zienkiewicz and Cheung[2], ArgyrisU] and Ahmad and Zienkiewicz[4]. In each of these studies assumed displacement finite elements were used. The hybrid-stress model, originated by Pian[S], can also be used to develop the element matrices. Elements based on the hybrid-stress model can exhibit superior displacement convergence, and stress and strain distribution compared with the analogous assumed displacement elements as shown for example

TAssistant Manager (Analytical). SAssociate Professor of Structural Mechanics.

in earlier studies by Pian and Tong[6], Fian[7], Ahmad and Irons[8]. Spilker and Daugirda[9] have developed an element for semi-analytic analysis (AXH14H), shown in Fig. 1, using the hybrid-stress model. In nonlinear analyses the solution accuracy is governed by the method of solution as well as the type of finite element used. In material nonlinear analyses such as elasto-plastic, visco-plastic and creep the structural response is governed by the stress state of the material. Thus it is important to use elements which predict accurate stress distribution and hybrid-stress based elements are good candidates for such analyses. The material nonlinear analysis of axisymmetric structures subject to arbitrary loading is usually performed by using 3-D elements. Since com-

, 1,-. c-

/’

__--__-

--__

-__“4w4

yy3

Y

/ 7. (r3.23) H. u4 ,L____* , I r , (r4,z.J “3 \ / I \ v2 w2 I’ , Wl 1 “2 ‘-__ “1 --__ ’ b2’Z2) AC----. ,D I (‘,.Z,) Ul

Fig. 1. Geometry and node sequence for the 4 node isoparametric element AXH14H.

448

S. P. SINGHand R. L. SPILKER

putational cost and effort associated with a full 3-D analysis are prohibitive, some researchers have successfully extended the semi-analytic method to nonlinear analysis: works due to Stricklin[lO], Kotanchik et al. [l 11, Meissner[12] and Winnicki and Zienkiewicz [13] are based on the assumed displacement elements. In this paper the hybrid-stress element developed by Spilker and Daugirda[9] is selected for extension of the semi-analytic method into elasto-plastic analysis with the hybrid-stress formulation. To evaluate the accuracy and relative computational effort of the semi-analytic method compared with a full 3-D analysis, a hybrid-stress 20-node solid element (called SH69Q) developed by Spilker and Singh[ 141 is also implemented for elasto-plastic analysis. Node numbering and geometry for this isoparametric quadratic element are shown in Fig. 2. Generally, incremental theory of plasticity is used in the elasto-plastic analysis. In conjunction with this theory, three solution techniques are most popular: Tangent Stiffness, Initial Stress and Initial Strain methods. The latter two require assembly and factoring of the global stiffness matrix only once, thus leading to computational savings. In the hybridstress approach the initial stress method is prefered since it does not require determination of strains from incremental stresses during plastic flow in elastic, perfectly-plastic materials. In this paper the general hybrid-stress based formulation for elasto-plastic analysis is briefly summarized, followed by the procedure for extension of the method to semi-analytic finite elements. Procedures used in the calculation of equivalent load vectors due to initial stresses are explained and the convergence criterion used is discussed. Finally two example problems solved by a full 3-D solution method and the semi-analytic method are presented. The elasto-plastic analysis is based on a hybrid-stress formulation in conjunction with the initial stress approach. Alternate functionals of this type were developed by Spilker and Plan [151and used in the analysis of axisymmetric solids. Spilker and MunirI16] also expanded the alternate functionals for plate bending analysis. The hybrid functional n, labekd xfik, was apparently prefered in both studies and is adopted for the present study. Briefly, hybrid functional n: is distinguished in that the elastic

stress increment during plastic flow is assumed to satisfy equilibrium in the interior of an element. The primary objective of this study is the application of this functional to the semi-analytic analysis, A brief summary of the functional, matrix definitions, and procedures are presented; further details are found in Ref. [15] VARIATIONALPRINCIPLE

The functional nz

for a 3-D solid is given by [15]: &,,,,,A~;~ds~ dv

V” + Iv:{,Atk,du

+

-

jvn

As;,&,dv

i=kAukds + R,*

s &m

(1) 1

tensor; where: S,, = linear elastic strain-stress s; = elastic stress tensor; & = fictitious initial stress tensor; ek[= Strain tenSOr determined from displacement field; u, = kinematically admissible displacement field in the interior of the element; r, = prescribed boundry traction; A = increment in quantity ( ); (s& = elastic-plastic stress tensor = s;,- S&. The stress and strain quantities are illustrated for a hypothetical one-dimensional case in Fig. 3. At a load increment (step i), the total incremental stress Asj is known and the strain de, is determined through the linear stress-strain relations. On the basis of strain be, and current stress (s” + As;), As,, can be determined via plasticity relations. Since As’ and AS+, are known, .rf is computed as the difference between As’ and As,. The term R,+ in eqn (1) represents an equilibrium imbalance in element “n” at load step ‘7” and is given by

Uniaxial strain(e) Fig. 2. Geometry and node sequence for the 20 node isoparametric element SH69Q.

Fig. 3. Schematics of stress quantities in initial stress approach for a loading increment.

Elasto-plastic analysis of axisymmetric structures where ( )’ = total value of the quantity ( ); (F& = the total body force. Basically, eq (2) is a virtual work statement of the equilibrium imbalance expression. MATRIX FORMULATION (SH69Q ELEMENT)

In hybrid formulation rr:, stresses As’ are interpolated in V, in terms of stress parameter Aa in the form As’ = PAa such that equilibrium

(3)

is satisfied, i.e.

As;,,, = 0

(9b)

Equation (9a) is written on element level since the da’s are independent from element to element. When eqn (9a) is solved for Aa and substituted into eqn (9b), the summation over all the elements leads to the following global system of equations for the ith loading step: K*Aq:=AQ,++Q;*+Ff*-e* =AQ:+R;*+Ff*

in I’,.

(10)

(4)

The intra-element displacements, u, are interpolated in V,, in terms of nodal displacements, Aq, of the element as: u=NAq

(5)

such that the appropriate interelement displacement continuity is satisfied. The strain-displacement relations in conjunction with eqn (5) yield: i = BAq.

(6)

Equations (2)-(6) are substituted into eqn (1) and the following matrices are defined: H= r PrSPdv

(7a)

0) Fd=

~~d-GTAa-F+AQ+Qo)i=O.

449

PrSfdv sY

(7c)

where K* = global stiffness matrix of the structures; Aq: = incremental global nodal displacement vector at load step i; A QF = incremental load vector due to applied loads at step i; RF* = global nodal load vector due to equilibrium imbalance correction; and P’* = global nodal load vector due to fictitious initial stresses. For the ith increment (from applied load AQ,*_, to AQf) the applied load increment is known. The equilibrium imbalance correction term applies to the state of stress/strain at the beginning of the step and is therefore also known. However, the equivalent nodal loads corresponding to the fictitious initial stresses are dependent on the stress state at the end of the ith increment and therefore are not known. For a small increment, it may be possible to approximate c* by Ff? ,, so that eqn (10) can be written as: K*Aqt=AQ;+R::,+c:,.

(11)

In the solution of eqn (1 l), the global stiffness matrix K* is factored prior to the first loading step.

Ud) MATRIXFORMULATlON FOR THE SEMIANALYTIC METHOD (AXH14H ELEMENT)

(74 Q"=

s

m

NT&.

sn

The functional the form:

rrf (eqn 1) then can be expressed in

nk(Aa,Aq) =c

n

{l/2 AarHAa + AqW-

aTGAq

-Aq%“+AqTAQT+AqTQo}.

(8)

Since nk is stationary with respect to arbitrary nonzero variations of Aa and Aq, the variation S&, (for 6a # 0 and 6q # 0) leads to the following equations: wAali = [GAqJ

(94

tNo distinction is made here between symmetric and anti-symmetric harmonics, and summation first extends over the elements and then harmonics.

The efficiency of the semi-analytic method is a consequence of the representation of stresses, strains, displacements and loads by a Fourier series. The harmonics present in the series are dependent on the arbitrary load distribution. In linear analysis, the total structure response is obtained by superposition of contributions of individual (uncoupleo, harmonics. In elastoplastic analysis, nonlinear effects (plastic flow) are due to the cummulative effect of the harmonics. Thus some procedure to extract the contribution of each harmonic to the nonlinear effects is required. In extending the semi-analytic approach to nonlinear analysis, an assumption is made that nonlinear effects, which are cummulative effects of the harmonics used in the analysis, are expandable in a Fourier series. The series in general contains the same harmonics as used to approximate the load distribution. With this assumption, the variational principle used may be expressed as a sum of variational principles for individual harmonics present in the analysis. The functional, s&, of eqn (8) can then be written as:?

S. P. SINGH and R. L. SPILKER

450 ~~[~

follows:

a,, A Gl

KZ(A qZ)i = [(A Q*)m f (QO*)m- F’*>m + (Fd*)mli

=~{~[1/2Aa~~~Aa,+Aq”L,l,d

= [(A Q*>m+ (R”*)m + (E”*)mli

Aa,TG,,,LL,A~T - A$f,%,w

-

+ A s*T,,J+=A Q, + A q*TJ,TQ,,,“l

(12)

where definitions of a,,,, H,,,, G,, q,, F,,,‘, F,,,’ and Q, are the same as defined in eqns (lb(7), except that these quantities are now defined in conjunct;on with harmonic “m” (for detailed definition of H,,,, G, and Q,, the reader is advised to see [9]). The quantities Fmd, F,,,’ and Q,,,’ are defined as: F,d = C,,, B,,,=(r,z)s,‘dr dz ISr I

(13a)

F,,,O= C,

(13b)

BmT(r,z)s,,,Odr dz

s

NT,,,‘d.s

Q; = C,,,

where Kz = global stiffness matrix; Aq*m= incremental nodal displacements vector; AQZ = incremental nodal load vector due to applied loads; R?=nodal load vector due to equilibrium imbalance correction; and Fi* = nodal load vector due to fictitious initial stress. Following the discussion of eqn (lo), eqn (15) can be written in the form; K.Z(Ati), = (A Q3i + i1+ (F”,*)t- 1. (16) At each increment, eqn (16) is solved for each harmonic. The total displacements, stresses and strains are then computed by superposition. At the first load step, Kz for each harmonic is factored and stored (on disc in our case) for further use.

(13c)

S”

CALCULATION

OF FICTITIOUS

STRESS AND EQUIVALENT

where

INITIAL

LOADS VECTOR

FOR THE SH69Q AND AXH14H ELEMENTS

211 m=O II m>l’

c,=

and B, = strain-displacement matrix for harmonic “m”; s, ’ = fictitious initial stress vector for harmonic “m” determined from the Fourier series expansion of total 9/ for a given (r, z) point along the circumference; s,” = total stress vector for harmonic “m” determined from so (as for smf); T,,,O= total traction for harmonic “m”. The vectors Fd and F” are computed from fictitious initial stress sf and total stress so. In the semi-analytic method, these stresses must correspond to the cummulative effect of all the harmonics used. In order to write the functional, nk, in the form of eqn (12), Fmd and F,’ must be computed from s,f and s,,,‘, which in turn, should be obtained from the total quantities; sf and so. A procedure for computing these stresses will be discussed later. Since nonlinear effects can be decoupled into harmonic components (assumption) and eqn (12) represents a superposition of functionals for individual harmonics, the stationary property of rr& can be used. By taking the first variation of xz for each harmonic with respect to a, and q,,, (for da, # 0 and 6q, # 0), the following set of equations can be obtained; (for mth harmonic) w,,,Aa, - G,L,Aq,,J = 0 c [L,qmd - L,=G,,,=Aa, - L,qmo n + L,=AQ, + L,=Q,q = 0.

(144

(14b)

Solving for Aa,,, from eqn (14a) and substituting into (14b) and carrying out the summation over the number of elements leads to a global system of eouations for harmonic “m”. (ith loadine sten) as -&--~~---

(15)

~-~

~~~

-~~~~

I

\

In the initial stress approach for elasto-plastic analysis, the accuracy of the analysis is dependent on accurate estimation of the fictitious initial stress, and the computaion of equivalent nodal loads. A method to compute these quantities in conjunction with the 2-D (axisymmetric) analysis and in plate bending application has been described by Spilker and Pian[ 151 and Spilker and Munir [ 161, respectively. The procedure can be easily adapted for the 3-D (SH69Q) analysis. Here, a description of the extension of this method in conjunction with the semianalytic (AXH14H) model is presented; additional details are given by Singh [17]. At a load step ‘i” eqn (16) is solved for the incremental displacements Aq,,, for an element for m th harmonic, and the incremental stress parameters Aa, are computed via eqn (14a). The incremental (elastic) stress AsA and strains Ae; at a point (r, z, 0) are computed through the relations; Ask(r, z, 0) = A,(@)P,(r,z)Aa,

(17)

Aeh(r, z, 6) = SAsk.

(18)

The total incremental stress (As’)~and strain (Ae’)i are determined by superposition of the contributions of all the harmonics as: (As’)~= CA s:,

(1%

(Ae’)i = 1 A eh

(20)

and

The total stress so and total strains e” at a point (r, z, 13)are updated by adding the incremental quan-

Elastoplastic analysis of axisymmetric structures tities to the current states as: so@,z, 6) = sy- i(r, z, 6) + As’(r, z, e)

(21)

e”(r, Z, 6) = ey_,(r, z,e) + de’@, z, e).

(22)

451

integrating the corresponding stress coefficients s,‘(r, z) and s,‘(r, z) over the plane (r, z) for every element, as: (Fnt?i

=

BmT(r, z)s{(r, z)r dr dz

cm

(28a)

ss ,z

Thus using eqn (21) and (22), stresses and strains in a given element are computed at the discrete circumferential locations for all the (r, z) points selected as stress/strain points. For the AXH14H element, the centroid of the element was found to be the optimal stress-strain point [9], and is used in this study. The yield criterion is checked at every (r, z, e) sampling station using the total stress so and the elastic-plastic stress s,Jr, z, 0) is determined via the plasticity relations [ 181; Ase,, = D,Ae’.

(23)

The total stress so and fictitious stress $ are determined as follows. sf(r,z,e)=s~_,(r,z,e)+As,(r,z,e)

(F,O)i=cm

B,,,r(r, z)sf(r, z)r dr dz I z

t28b)

IS

(F,,,d)i and (F,“)i are assembled to obtain global quantities (Fi*)i and @m*)P This procedure as described from eqns (17)-(28) has been presented in Fig. 4 as a flow-chart. From this discussion it is apparent that for an accurate approximation of Fourier coefficients, a suflicient number of sampling stations in the 0 direction must be used. For the problem in this study where more than a single harmonic is used, the sampling stations are located at 5” intervals. Only two points are used for the case where a single harmonic is present. The distribution of the load along the 8 direction governs the number of these points. CONVERGENCE CRITRRION

(24)

In incremental solution schemes, accuracy of the solution is governed by the size of the load increment. s[(r, Z, 0) = As’(r, Z, 0) - As&r, Z, 0). (25) Since the use of smaller load increment for convergence would lead to prohibitive solution cost, Equations (24) and (25) represent the quantities some convergence criterion should be used. In this which are the cummulative effect of all the harstudy, the incremental-iterative scheme used in Ref. monics. To extract the harmonic components (s@)[ [15] is used for the 3-D (SH69Q) analysis and a and (s/)~, the distributions of total quantities along modified form is used in the semi-analytic (AXH14H) the 6 direction are used to compute Fourier analysis. coefficients s,,$, z) and s,‘(r, z). Since values of the For illustrative purposes, consider eqn (11) for the quantities are given at discrete 8 locations, some kind ith load step and kth iteration, which can be written of numerical method must be used to determine these as: coefficients. In order to do so, (s”)~and (4/)i at a point (r, z) can be written as a Fourier series, i.e. K*A@=AQ:+R$+F$*. (29) sO(r, z, 6 ) = [SO@,zllo+

-%A~ W’(r, z X, + MU%% 41;

si(r, z, 0) =
In this scheme eqn (29) is solved repetitively for the “i” until the equivalent load Fp corresponding to the fictitious initial stress is sufficiently small, i.e. until:

(26) load step

(27)

where [ lo = zero coefficient of the series (4); [ If,,= mth symmetric coefficient of the series (a,,,); and [ 1; = Ith antisymmetric coefficient of the series (0 4, a, and 6, correlate the coefficients to a convenient Fourier series g(6) = a, + a,,,cos (me) + b,sin(Z@) and eqns (26) and (27) are written for a general case where symmetric and antisymmetric harmonics are present. A,(B) and AX@) are square matrices of order 6 with nonzero terms only on the diagonals. These are taken from Ref. [9] and are defined in the Appendix of this paper. As mentioned [ lo, [ ]s, and [ 1; are to be determined from a given distribution of the quantities along the 8 direction. Here, the scheme used by Zienkiewicz et al. [ 131, Bodie’s Rule [ 191, was adopted for computing these coefficients. When the coefficients for all harmonics are known, the load vectors F,,,d and F,O can be computed by

(30) where 1 1 denotes the magnitude (squared) of the vector and TOL is a small non-negative (e.g. TOL < 0.01) parameter. The extension of this scheme to the semi-analytic method, eqn (16) can be rewritten in a form similar to eqn (29) as: %(A qz) = (A Q3i

+ (E*)k + (Fd,*)t.

(31)

To apply the criterion of eqn (30) to the semi-analytic method, convergence is determined by using the cummulative magnitude of the fictitious initial stress load vector due to all harmonics, i.e. (32) Here, it4 is the total number of harmonics. Now eqn

452

S. P. SING&I and

&In

LOOP ON STRESS

POINTS n.n*1 t

R. L.

SPILKER

= 2 [A~‘il,

i roil,

= f~oi_ll,

+ [A%],

f 2'"

= r$,_,l,

+ iA&

ALL POINTS CHECKED FOR YIELDING

YES

-L-J.-NO

F [ btll,

NO

I

kl

1 LOOP ON I

30

STRESS POINTS

YES

IA.jpll

= [ D&

[A$,

I

I

1

I

I

I

LOOP OVER HARHONICS h.h+l

NO

ALL HARXONIC DONE

GAUSS POINTS psp+1

ALL GAUSS WINTS NO

Fig. 4. Flow chart for computing equivalent nodal loads for AXH14H (semi-analytic) model.

Elasto-plastic

analysis of axisymmetric structures

453

(31) is repetitively solved for every harmonic until the criterion of eqn (32) is satisfied. In this study a value of TOL = 0.005 was used for the problem requiring more than one harmonic. For the 3-D analysis (SH69Q model) a value of TOL = 0.01 was used. A smaller value of TOL was used in the semi-analytic analysis to avoid gross violation equilibrium (eqn 3 1) by any harmonic. EXAMPLR PROBLEMS

Two examples problems are selected to investigate the accuracy and convergence of elements SH69Q and AXHl4H in elasto-plastic analysis based on the formulation presented in the previous section. Both examples involve elastic perfectly-plastic materials. The first example is a thick circular cylinder subjected to uniform internal pressure, which has been selected because an exact solution to the problem is available. The second example is a thick cylinder subject to a nonuniform pressure load so that multiple harmonics are required in the semi-analytic analysis. In the first example a quarter of the cylinder is modelled for the 3-D analysis. The finite element mesh used for the SH69Q analysis is shown in Fig. 5 (i.e. 4 elements circumferencially, 2 elements radially and one element along the axis of the cylinder). In the AXHl4H model, five elements (5 elements along the radius and one element along the cylinder axis) are used as shown in Fig. 6. Due to symmetry of the load only the single zeroth harmonic is used in the AXHl4H model. The Young’s M~ulus, Poissons Ratio and Yield stress used for the elasticperfectly plastic analysis are 106psi, 0.3 and 20,000 psi, respectively. The exact solution of this problem is given by Venkatraman and Patelf20J. For the geometry used in this example (inner radius = 5 in, outer radius = loin) initial yielding occurs at the inner surface at an applied pressure of 8600 psi and a fully plastic state occurs at 16,000 psi. In both analyses, the pressure was increased from 8000 to 16,000 psi in 40 equal load increments. The convergence criterion of TOL = 0.01 was used for both models and both methods required iterations once plastic flow began. The last load increment which converged was number 39 in both models and corresponds to =98.75x of the collapse load. To compare the predictions of the two models, the

Y

IL t,

X

)\

0”

X

Fig. 5. Geometry and &rite element grid for the thick cylinder (SH69Q) model.

Fig. 6. Finite element grid for the thick cylinder (AXHl4H model.

radial displacement y, at the inner surface and radial stress s, and hoop stress s, in the innermost element are selected. For the SH69Q analysis, the stresses are evaluated at R = 5.527in and for AXHl4H at R = 5.5 in. These points represent the optimal stress-strain Gauss points for the two elements. The results obtained for u,, s, and s, versus applied pressure along with the exact solution, are plotted in Figs. 7(a-c), respectively. It can be observed from Figs. 7(ac) that both models produce results for the displacement and stresses which closely follow the exact solutions. Since only a single harmonic is used in the AXHl4H model, the computation time required by the AXHl4H model and SH69Q differ by an order of magnitude. The CPU time used by the SH69Q model on a DEC-2050 Computer system was 33.8 min, whereas the AXHl4H model required only 3.4min. All computations were carried out in the double precision. The second example is a thick cylinder subjected to a harmonically distributed, inner surface pressure loading as shown in Fig. 8. This example is designed to demonstrate the application of the semi-analytic method in elasto-plastic analysis where the applied load distribution is approximated by more than one harmonic, and to provide a more realistic comparison of the relative efficiency of the semi-analytic and 3-D analyses. The expression for pressure distribution in the lirst quadrant (0 < 6 < 90”) is: p(e) = I

q4ej7t P(4eln

+ sin(2e))

0 < e < 7~14

-sin(20))

7r/4Geerrn/2

where P is 8OOOpsi and the peak value of p(0) is 16,000 psi at 8 = 45”. Similar pressure dist~bution is applied on the inner surface of the three other quadrants of the cylinder. This pressure distribution is designed to minimize the 3-D modelling effort (due to memory and computational constraints); i.e. due to symmetry about two axes, only a quarter of the cylinder is modelled in the SH69Q analysis. In the semi-analytic analysis (AXH14H model), the cylinder is modelled by ten elements radially and one element along the cylinder axis. The load distribution is approximated by using six symmetric harmonics (0,4,8, 12, 16,20). The 3-D model (SH69Q model) consisting of 32 elements (four elements radially, eight elements circumferencially and one etement along the cylinder axis not shown) is similar to one shown in Fig. 5, except the new grid

S.

454

0 AXHl4H

P.

SINGH and R. L. SPILKER

Model

SH690

l

3.2-

-

Model

Exact

2.62 & ZJ 24g El E 58 2.01.6-

. P

d .L” D 1.2 G -0 s

0.6-

1

0

16.0

Internal

Fig. 7(a). Radial displacement

0 AXHl4H

1

pressure (KS1

(u,) vs internal pressure.

Model

l

16.0-

SH69Q

-

Model

Exact

lO.Ozi c -

co-

fr I 2

m iii

6.0 -

Internal

pressure

(KS11

Fig. 7(b). Radial stress sr vs internal pressure.

(used here) has four additional elements along the circumferential direction and two additional elements along the radial direction. The material is assumed to be elastic perfectly plastic; material properties used are identical to the first example. Since there is no solution available for this problem, trial runs were made to determine the load level at which the solutions started diverging. The pressure was increased from 16,000 (peak value at ,$I= 45”) to 27,520 psi in 12 load increments. In both models it was found that the solutions diverged (corresponding to a fully-plastic load) at a pressure smaller than 27,520psi. Finally, a loading sequence was used in

which

the

pressure

was

increased

from

16,000

to

25,856 psi (higher of the two collapse loads predicted by the two models in trial runs) through 12 load increments. During iterations within a load step, a convergence criterion of TOL = 0.005 was chosen. During the collapse load step of each model, the solution did not converge within maximum of 30 iterations permitted. Since computational times for both AXH14H and SH69Q models were quite large, it was decided not to increase the number of load increments to establish the collapse load more exactly. The SH69Q model failed to converge in the 8th

Elasto-plastic 0 AXHl4H

455

analysis of axisymmetric structures Model

l

16.07

SH69Q

-

Model

Exact

14.0-

12.0-

2

lO.O-

Y I,

e.o-

E .a 07

6.0 -

4.0-

2.0-

0

/

1.6

3.2

4.6

6.4

Internal

6.0

9.6

11.2

14.4

16.0

pressure (KS11

Fig. 7(c). Hoop stress S, vs internal pressure.

1.3

2.5

3.6

5:O I 6.3

7.5

5.6

10.0

Radius (in.)

Fig. 8. Geometry and pressure distribution for the cylinder subject to harmonic internal pressure.

I

S. P. SINGH and R. L. SPILKER

456

load increment, and the AXH14H model at the 1lth load step. The two collapse loads correspond to pressure levels of 24,960 and 26,045 psi, respectively. The load deflection (radial displacement at 0 = 45”) curves obtained from two models are shown in Fig. 9. It can be seen that in early stages of yielding, both models are in close agreement. As load is increased, the AXH14H model predicts stiffer behavior than the SH69Q model. However, it is worth noting here that the collapse load predicted by the AXHl4H model is only 4.357” more than the SH69Q model. To further investigate the results from the two models, various stress distributions and the yielded region of the cylinder are plotted for two load increments. Note that the stresses predicted by the two models at the first load increment are essentially the same since the load did not cause yielding. Load increment five and the load increment corresponding to collapse (increment 10 for AXHl4H and 8 for SH69Q) are chosen for comparison. The radial stress, s,, tangential stress, s, and axial stress, s, have been plotted as a function of circumferential coordinate (0) in Figs. lO(a-c). The reported stresses are for the mnermost elements in both models. In the SH69Q model, 2 point Gauss stations are used to evaluate the stresses, and stresses plotted in the figures corresponds to a radial location of R = 5.264 in. In the AXH14H model, stresses are evaluated at the centroids of the elements and stresses plotted in the figures correspond to a radial location of R = 5.25 in. The yielded regions predicted by the two models for the chosen load steps have been plotted in Fig. 10(d). The solid line enclosing the interior of the cylinder quadrant denotes the plastic zone predicted by the SH69Q model, whereas, dashed lines represent the prediction of the AXH14H model. The distribution of radial stress, s,, (Fig. lOa), predicted by the two models at increment number 5 and 8 is in good agreement. The distribution of tangential stress, s,, (Fig. lob), from the two models are also in good agreement except near the bound-

o

AXHl4H

aries (0 = 0” and 8 = 90”). The same behavior is exhibited by the axial stress, s,. (Fig. 101~).In Fig. 10(d), yielded regions predicted by the two models are shown. It is apparent that the AXH14H model has a tendency to spread the plasticity along the tangential direction at a slightly faster rate than the SH69Q model. Though the boundaries of the yielded zones in the two models differ, the areas enclosed by them do not differ by a significant amount, which is indicative of only 4.35% difference in the predicted collapse loads. The distributions of stresses, s,, s, and s, at radii, R = 7.75 in (for AXH14H model) and R = 7.764in (for SH69Q model), for the two load increments have been plotted in Fig. 1l(a-c). Again it can be observed that at increment number five, the model predictions are close. At increment number eight, the SH69Q model results produce an approximate fit through the results of AXH14H model. It should be noted that AXH14H model results will have a smoother distribution than SH69Q model due to Fourier series expansion of the stresses and strains during the analysis. It was not possible to increase the number of elements in the SH69Q model due to constraints on memory space and increased computational cost, but a more refined 3-D model is apparently needed here to establish more accurately the “exact” solution. A final point should be made here regarding the Fourier series expansion of the plasticity effects (fictitious stress) in the AXH14H model. As is evident from the results of the SH69Q model and the applied load distribution, yielding would begin at a plane f3= 45”. Since a numerical scheme is used to compute the Fourier coefficients from a distribution of stresses at finite 0 locations, the nonlinear effects may be underestimated. This could cause the AXH14H model to average sharp peaks in stress over an increased interval in 0. This may be the explanation for the observed spreading of plasticity along the circumference.

l

Model

SH690

Model

6.00 r

Normalized

internal

pressure (P/PI)

Fig. 9. Radial displacement ur vs normalized internal pressure for (harmonic loading).

Elasto-plastic

0

451

analysis of axisymmetric structures

AXH14H

Model

0 Incramant

Radius

SH690

Modal

No. 5 Radius

= 5.25in.

= 3.234

in.

IncromontNo. S Radius

= 5.2Sin

Radius

Theta

= 5.254

in

(00)

Fig. 10(a). Radial stress s, vs B location for load increment No. 5 and 8 (near inner surface).

458

S. P. SINGHand R. L. SPILKJZR

0

AXH14H

.

Model Incremant

31.0

Radius

SH690

Model

NO 5 Radius

- 5.25 tn.

= 5.254

in.

26.5

22.0

17.5 13.0

6.5

4.0

Increment Radius

No. 6

= 5 25in.

Radius - 5.254in

26 .S

t

17.5

13.0

6.5

4.0

-0.5

-5.01 910

’ 16.0

I 27.0

I 36.0

I 45.0 Theta

I 54.0

I 63.0

I 72.0

, 81.0

J 90.0

IDG)

Fig. 10(b). Hoop stress s, vs tI location for load increment No. 5 and 8 (near inner surface).

Elasto-plastic

0

analysis of axisymmetric structures

AXHI4H

Model

. Incramont

Radius = 5.25 in

Incremsnt Radius

2 5.25in.

St+690

Model

No, 6 Radius

: 5.254

Radius:

5.254in

in.

Na6

Fig. 10(c). Axial stress s, vs 0 location for load increment No. 5 and g (near inner

surface).

S. P. SINGHand R. L. SPIJXEIZ

---

AXHI4H

Modok

-

SH69Q

Modal

Plastic-regime Incremsnt

No. 5

8.8 -

7.5-

2.5-

lncremartt

Radius

No.8

(in.)

Fig. 10(d). Yielded region for load increment No. 5 and collapse load increment.

Elasto-plastic

0

AXHI4H

Modrl

0 Incremant

7.0

Radius

461

analysis of axisymmetric structures

6~690

Modat

No.5

- 7.75in.

Radius

- 7.7542

Radius

= 7.754

in.

3.0

-17.0

-21.0 i

L E” .? m

Increment 7.0

r

Radius

No.6

I 7.75 in.

in

-17.0

-21.0 i -25.0

I 9.0

I 16.0

1 27.0

I 36.0

I 45.0

Thata

1 54.0

I 63.0

I 72.0

I GID

I SO.0

LDG)

Fig. 11(a). Radial stress s, vs 0 location for load increment No. 5 and 8 (near mid-height).

462

S.P.

0

SINGH

AXHl4H

and R. L. SPILKER

Modal

. Increment

Radius

= 7.75

I

1

SH690

Model

No.5 Radius:

in.

7.7542

26.5

6.5 t 4.0

-0.5 1

i

-6.01

I

I

1

I

I

I

I

I

Y

c 8 D ;i

Incremmt 31.0

Radius

= 7.75

No.6 Radius - 7.754

in.

In.

26.5

Thoto

(001

Fig. 11(b). Hoop stress S, vs 0 location for load increment No. 5 and 8 (near mid-height).

Elasto-plastic

0

analysis of axisymmetric structures

AXttl4H

Mod.1

0 Incrrm.nt

Radius

= 7.75

E”

Radius:

7.7542

in.

No.8

IncrNnant Radius

14.0

Hod.1

No.5

in.

0

;?

SH69Q

463

: 7.75in.

Radius

: 7.754in

-7.01

-10.0

I

9.0

I

18.0

1

27.0

I

36.0

1

4s.G

Theto

I

540

I

93.0

I

720

,

El.0

I 9G.G

IDG)

Fig. 1l(c). Axial stress s, vs 0 location for load increment No. 5 and 8 (near mid-height).

R. L.

The computational times required by the two models in the complete analysis were not as different as in the previous example. The total time used by the SH69Q model for eight load steps was 77 min. The AXH14H model required 80min for eleven steps. Here, it is worth noting that this problem had been designed with symmetry about two axes for the 3-D model, which contributed to substantial reduction in computational time for the SH69Q model. Still, the 3-D model required 32% more time per load increment compared with the semi-analytic analytic analysis.

4.

5.

6.

CONCLUDING REMARK!3

The elastic-plastic analysis of axisymmetric structures subject to arbitrary loading has been examined by using an initial stress approach in conjunction with the hybrid-stress model. The semi-analytic approach, in which quantities of interest are expanded in a Fourier series in the circumferential direction, has been used and the appropriate hybrid-stress functional for elastic-plastic semi-analytic analysis has been defined. This allows the three-dimensional nonlinear problem to be solved as a series of twodimensional analyses. To evaluate the accuracy and efficiency of the semi-analytic method, a 3-D 20-node hybrid-stress element developed earlier by the authors has also been implemented for elastic-plastic analysis. The accuracy of the two, semi-analytic element AXH14H, and 3-D element SH69Q, were first estab-

lished by comparing the solutions with the exact solution. Another numerical example has been presented which exploits the use of the semi-analytic method. Six harmonics were used in the semi-analytic model and a full 3-D solution was obtained by the SH69Q element model. The collapse load for the structure was obtained along with the detailed stress distribution and yielded region at different load levels. A comparison of the semi-analytic and the full 3-D solution demonstrated close collapse load predictions. Further comparison of the stress distribution and yielded region during loading and just prior to collapse, showed a good match between the semianalytic and 3-D solution. In addition substantial savings in the effort of model preparation and computational time are observed through the use of the semi-analytic method. Acknowledgement-The authors wish to thank the Association of American Railroads, Chicago, Illinois, which permitted the (first author) use of the computational facilities for the work reported here. REFERENCES 1. E. L. Wilson, Structural analysis of axisymmetric solids. AIAA J. 3, 2265-2214 (1965). 2. 0. C. Zienkiewicz and Y. K. Cheung, Stress in shafts. The Engineer 24 (1967). 3. J. H. Argyris, K. E. Buck, I. Grieger and G. Marcezek,

7.

SPILKER

Application of the matrix displacement method to analysis of pressure vessels. Trans. ASME, 92B, 317-329 (1970). S. Ahmad, B. M. Irons and 0. C. Zienkiewicz, Curved thick shell and membrane elements with particular reference to axisymmetric problems. Proc. 2nd Conf. Matrix Methods ’ Structural Mechanics, Wright-Patterson A. F.rgase Ohio (1968). T. H. H. Pian, Hybrid Models in Numerical and Computer Methods in Applied Mechanics (Edited by S. J. Fenves, N. Perrone, A. R. Robinson and W. C. Schnobrich). Academic Press, New York (1973). T. H. H. Pian and P. Tong, Basis of finite element methods for solid continua. Inl. J. Num. Merh. Enann _Y I, 3-28 (1969). T. H. H. Pian, Variational principles for incremental finite element methods. J. Franklin Inst.. 302., 473-488 (1976).

8. S. Ahmad and B. M. Irons, An assumed stress approach to refined isoparametric finite elements in three dimensions. Finite Element Methods in Engineering, DD. -- __ 85-100. University of New South Wales, U.K. (1974). 9. R. L. Snilker and Diane M. Dauairda. Analvsis of azisymmetric structures under arbitrary loading using the hybrid-stress model. Int. J. Num. Meth. Engng 17, 801-828 (1981). 10. J. A. Stricklin et al. Linear and non-linear analysis of shells of revolution with asymmetric stiffness properties. Proc. 2nd Co& Matrix Methods in Structural Mechanics. Wright-Patterson A. F. Base, Ohio (1968). 11. E. A. Witmer and J. J. Kotanchik, Progress report on discrete element elastic and elastioplastic analyses of shells of revolution subjected to axisymmetric and asymmetric loading Proc. 2nd Conf Matrix Methods in Structural Mechanics, Wright-Patterson A. F. Base, Ohio (1968). 12. H. E. Meissner, Laterally loaded pipe piles in cohesionless soil. Proc. 2nd Cont Numerical Methods in Geo-Mechanics, (Edited by C. S. Desai), Vol. III, Virginia (1976). 13. L. A. Winnicki and 0. C. Zienkiewicz, Plastic (or visco-plastic) behavior of axisymmetric bodies subjected to non-symmetric loading-semi-analytical finite element solution. Int. J. Num. Meth. Engng 14, 1399 1412 (1979). 14. R. L. Spilker and S. P. Singh, Three dimensional hybrid-stress isoparametric quadratic displacement elements. Int. J. Num. Meth. Engng 18, 445466 (1982). 15. R. L. Spilker and T. H. H. Pian, Hybrid-stress models for elasto-plastic analysis by the initial stress approach. Int. J. Num. Meth. Engng 14, 359-378 (1979). 16. R. L. Spilker and N. I. Munir, Elasto-plastic analysis of plates by hybrid-stress model and initial stress approach. Znt. J. Num. Merh. Engng 17, 1791-1810 (1981). 17. S. P. Singh, Material nonlinear analysis of 3-D and axisymmetric structures (under arbitrary loads) using hybrid-stress finite elements, pp. 8489. Ph. D. Dissertation, University of Illinois, Chicago (May, 1982). 18. 0. C. Zienkiewicz, The Finite Element Methods, 3rd Edn, p. 493. McGraw-Hill, New York (1977). 19. M. Abramodvitz and I. A. Stegun, Mathematical Function, 9th Edn. Dover, New York (1970). 20. 6. Venkatraman and S. A. Pate& Structural Mechanics with Introduction to Elasticity and Plasticity. McGrawHill, New York (1970).

APPENDIX

Symmetric harmonic “m” -oos(mQ 0

0 cos(me)

0 0

0 0

0 0

0

0 80

0 00

0 cos(m0) 00

0 cos(me) 00

sin(mf7) 00

0 sin (me) 0 1

0

1

I

A,,,@)=

Elasto-plastic

analysis of axisymmetric structures

465

Antisymmetric harmonic “I”

A,@) =

1 sin(N) 0

sin0(re)

sin0(10)

sin(B) 0

-cOs(le) 0

-c0qe) 0 !