Evaluation of solution schemes for nonlinear structures

EVALUATION OF SOLUTION SCHEMES FOR NONLINEAR STRUCTURES D. P. Division

of Structural

Engineering

MONDKAR~

and Structural

and G. H. POWELLS

Mechanics,

Department

of

Civil Engineering, University of

California, Berkeley,CA 94720,U.S.A. (Received

I February 1976; received for publication

16 September

1977)

Abstract-This paper investigates solution schemes of the type that can be implemented in general purpose computer codes for nonlinear finite element analysis. A very flexible solution strategy is discussed, in which a variety of solution schemes can be implemented by specifying the values of certain solution control parameters. A number of nonlinear structures with diverse nonlinear characteristics are selected and each is analyzed using different solution schemes. The results and solution efficiency are then compared to provide guidance in the selection of the most appropriate solution procedures for various types of nonlinear behavior. 1. INTRODUCTION

number of theoretical and computational advances in the analysis of nonlinear structures have been made in recent years. Many computer analyses have been carried out, and a number of general purpose computer programs for nonlinear elastic and inelastic analysis have been developed. Nevertheless, a great deal of investigation remains to be carried out, particularly to determine the most economical solution schemes. The purpose of this paper is to investigate solution strategies of the type which can be implemented in general purpose computer codes for nonlinear finite element analysis. A preliminary version of a general purpose code for static and dynamic nonlinear analysis has been developed by the authors. This code incorporates a very flexible solution strategy, in which any one of a variety of solution schemes can be implemented by specifying the values of certain solution control parameters. The available schemes include step-by-step solution. Newton-Raphson iteration, constant stiffness iteration, and an adaptable scheme in which the structure stiffness is reformed only when needed. The flexible strategy has been designed to give the experienced analyst the opportunity to select the most effective scheme. because no single solution scheme appears to be optimal for all types of nonlinear structure. In this paper a number of structures with a variety of nonlinear characteristics have been chosen, each has been analyzed using a number of different solution schemes, and the results and solution efficiency have been compared. These comparisons demonstrate that different schemes may differ substantially on the basis of accuracy and efficiency. They provide guidance for the selection of the most appropriate solution scheme for problems of different types. It is believed that the results reported in this paper will supplement results obtained by other investigators to determine economical solution schemes for nonlinear structures (for example, see Ref. [13. 161). A

large

equations of motion for a structure undergoing large displacement, finite strain deformation has been presented in [l], using a Lagrangian frame of reference in which all static and kinematic variables are referred to the initial undeformed state. These discrete incremental equations for an undamped system are M.ii+[‘~~+‘~~1.q_=2p-(M.‘~+‘R).

(1)

Here &f is the structure mass matrix; ‘& is the linear stiffness matrix and ‘&. is the structure geometric stiffness matrix, both of which are evaluated in the current deformed configuration C, (time f); the state of motion in configuration C!, is defined by displacement ‘a velocity ‘4 and acceleration ‘d; the corresponding incremental quantities between the current configuration C, and the neighboring configuration C2 (time I + Ar) are e cj and 4; ‘11 is the load applied in configuration C2 and ‘6 is the load in equilibrium with the state of stress in the current configuration C,. It should be noted that eqn (1) includes equilibrium corrections to compensate for error due to linearization assumptions in deriving the equations of motion, and it is assumed that the prescribed surface tractions and body forces are not dependent on deformation, i.e. the load ‘p is independent of the deformation of the body. Material nonlinearity may be included while evaluating IKE, by specifying a tangent constitutive relationship between increment of stress and increment of strain; both the stress and strain are measured with respect to the undeformed configuration C, at time t = 0. Appropriate modifications to eqn (1) need be made if viscous damping is included. Details of the theoretical formulation are presented in [I], and may also be found elsewhere[2-4]. The equations of equilibrium for static loadings, or the equations of motion when integrated using an implicit, single step method[5] reduce to the following form: K*. q=f.*

(2)

2. COMPUTATIONAL ASPECTS 2.1

Discrete

incremental

equations

The finite element formulation

where &* is the tangent stiffness matrix in static analyor the effective stiffness matrix (which includes contribution of the mass and damping matrices) in dynamic analysis; similar interpretation can be attached. to the vector p. The forms of K* and f* are presented in [ 11,using Newmark’s method for dynamic analysis. sis,

of the incremental

+Assistant Research Engineer. $Professor of Civil Engineering. CIS\0l

9Llo

5-4

223

D. P. MONDKARand G. H. POWELL

Typicaliy, the structure response will be computed by applying the load in small steps. However. because eqn (2) will give only an approximate solution for the displacement increment. equilibrium iterations may have to be carried out in some cases to obtain results with a sufficient degree of accuracy.

An iterative procedure can be assumed to have converged when the unbalanced load (right hand side of eqn 1) becomes acceptably small, judged by the Euclidean norm or some other property of the vector. Convergence may also be based on the increments of nodal displacement or element strain.

2.2 Step-by-step and iteration procedures In the step-by-step solution procedure the load is applied in several small steps and the structure is assumed to respond linearly within each step, the response being obtained without iteration (Fig. la). This procedure is simple to apply and has been widely used, parti~~~iy in elasto-plastic analysis. It is essentially a Euler method, and errors are likely to accumulate after several steps unless very fine steps are adopted. The solution may, therefore, diverge considerably from the true response. The accuracy of the computed response can be improved by applying equilibrium corrections (Fig. lb). Two types of iterative procedure are commonly used, namely Newton-Raphson iteration (Fig. lc) and Constant Stiffness iteration (Fig. id). In Newton-Raphson iteration the structure matrix is reformulated at every iteration, and a disadvantage of this procedure is that a large amount of computational effort may be required to form and decompose the stiffness matrix. In Constant Stiffness itera~on the stiffness matrix is formed only once. However, Constant Stiffness iteration will typically converge more slowly than Newton-Raphson iteration, and schemes to accelerate convergence may be desirable. It may also be advantageous to devise mixed iteration schemes combining the features of both iteration procedures.

2.3 Computational steps A ^solution scheme for the nonlinear analysis of a structure involves three major computational steps. namely (I) linearization, (2) equation solution and (31 state dete~i~ation. The evaluation of the tangent stiffness matrix in static analysis, or of the effective stiffness matrix in dynamic analysis, constitutes the Iine~i~ation phase. Typically. in a finite element form~iation the structure matrix will be obtained by assembly of the element stiffness matrices. and requires the material constitutive relationship. For eiasto-plastic materials, this relationship varies with the state of stress and strain, and hence will vary over the element volume. If numerical integration is used to compute the element stiffness, it is convenient to evaluate constitutive relationship at each integration point. The computation of an effective load vector, and the solution of a set of linear equations constitute the equation solution phase. Form&as for the load vectors in static and dynamic analysis are presented in [I]. Techniques for solving linear equations are we11known and need not be discussed here 161. After eqn (2) has been solved for the displacement increment e it is necessary to compute the corresponding increments in stress and strain, and hence obtain a new state, This is the state determination phase. This phase is at least as important as the linearization phase. but has received considerably less attention. Computation of strain increments from displacement increments involves only kinematics, specifically the straindisplacement transformation. The problem of computing a stress increment from a given strain increment involves the material constitutive relationship. For hypoelastic materials, such as elast~plastic materials, the stress increment from an elastic to a yielded state, or during continued plastic loading, needs to be computed using numerical integration. Various algorithms can be designFd for this purpose, some of which have been presented in [ 1,7]. The question of path dependence assumes an added significance for these materials. A “true” loading path can be followed for such materials by using the step-by-step procedure with very small steps. If. however, an iterative procedure is used, a strain increment is calculated for each iteration and the increment for the load step is obtained by summation. In this case, it is important to distinguish between “path dependent” and “path inde~ndent” state determinations, as described in the following section. It should be noted that the compu~tional cost of the state determination phase can be significant compared with that of the linearization and equation solution phases, and therefore this cost shouId be included when evaluating the computational efficiency of a solution scheme.

P

COMPUTED RESPONSE

COMPUTED RESPONSE

P t

3P

3P

2P

2P

‘P

‘P RESPONSE

(0)

STEP-BY-STEP WITHOUT EPUlLleRlUY CORRECTION

(b)

STEP-BY-STEP WITH EPUILIBRVJY CORRECTlON

2.4 Path dependent and path independent 9

(cl

NEWTON-RAPHSON ITERATION

9 (d)

CONSTANT ITERATION

STIFFNESS

Fig. 1. Step-by-stepand iterative procedures.

state deter-

mination

With path dependent state determination, the stresses are updated at the end of each iteration, based on the strain increments computed for that iteration. In con-

Evaluation

of solution schemesfor nonlinear structures

frast, with path independent state determination the stress increment is computed for all strain increments accumulated up to any iteration, and the stresses are updated only after the iteration process has converged. For cases in which the strains increase monotonically during iteration, the results with the two types of state determination can be expected to be in close agreement (Fig. 2a). However, if the strains do not increase monotonically. the two methods may give significantly different results (Fig. 2b). Path dependent state determination appears to be more consistent for NewtonRaphson iteration, whereas for Constant Stiffness iteration it is more logical to use path independent state determination. 2.5 Imposed loads us imposed displacements For limit load analysis of a yielding structure, it can be expected that an iteration solution will tend to converge slowly at a load approaching the ultimate capacity of the structure, and for loads exceeding the ultimate capacity convergence can obviously not be obtained. This problem of sensitivity can be overcome in some cases by imposing increments of displacement, rather than force. It is important to note that with imposed disptacements, the changes in strain during iteration will not typically be in the form of progressively increasing strain over the entire structure, as in the case of imposed loads, but in the form of redistributions of strain involving strain increase in soft parts and corresponding decrease in stiff parts of the structure. For this reason, it may be advisable to use constant stiffness iteration with path independent state determination in most imposed displacement situations, because path dependent stress computation may lead to erroneous results. 2.6 Iteration with acceleration An acceleration scheme to improve the convergence rate in Constant Stiffness iteration has been described by Nayak and Zienkiewicz@l. In this scheme, called the alpha-constant acceleration scheme, the disptacement increments during any iteration are scaIed in an attempt to obtain the same result as if Newton-Raphson iteration were employed. For each iteration the scheme involves two steps of equation solution and two steps of state determination. A modi~ed version of the scheme has been implemented in the computer program, in which each iteration involves only one step of equation solution and one step of state determination. Further, the scaling matrix is reinitialized to a unit matrix if the iteration fails 10 converge within a specified number of iterations or if the Euclidean norm of the unbalanced load vector increases from one iteration to the next. Thess ~estr~ctions are intended to prevent possible divergence. Details of the modified version may be found in [I]. The scheme can be expected to improve the rate of A

iai

W,,

PATH

lb)

Fig. 1. Path dependence.

INOEPENDENI

22s

convergence for moderately nonlinear structures in which the load-displacement curve is of a softening type. However, for structures with stiffening load-displacement response, or for loadings which produce stress reversals. the acceleration scheme has been found to be unreliable and should be used with caution. The scheme can be used in both static and dynamic analysis. 3. SOLUTION STRATEGY the general case, the nonlinear response will be computed by a combination of the step-by-step and iterative procedures. Depending on the degree of nonlinearity, equilibrium iterations may or may not be needed, and the tangent or effective dynamic stiffness matrix may or may not be formulated in every step. Because of the large computational effort typically required for decomposition of the tangent stiffness matrix, it will usually be desirable to seek a strategy in which the number of stiffness matrix reformulation is minimized. Nevertheless, it is also necessary to consider the computational effort for the state determination calculations. Basically, there are only two solution procedures included in the computer program, namely (a) NewtonRaphson iteration and (b) Constant Stiffness iteration, the step-by-step procedure being treated as a special case of Newton-Raphson iteration. It is possible to obtain a variety of solution schemes by specifying values for a number of control parameters, as described in the following section. In

3.1 Control parameters of the solution strategy The control parameters are as follows: (a) Number of steps (NSTEP). For static .analysis, this parameter is the number of equal load steps. For dynamic analysis, it equals the number of time steps. (b) Type of ~~e~a~~onp~cedure ffTYP). This parameter takes a value of zero for Newton-Raphson iteration and a value of one for constant stiffness iteration. If a value greater than one is assigned, constant stiffness iteration with the alpha-constant acceleration scheme is used. The value of the parameter is then used to control reinitialization of the scaling matrix. Reinitialization to a unit matrix is carried out if the number of iterations before convergence exceeds the value of the parameter. (c) Type of slate dererminarion (KPATH). This parameter controls the state determination process for elasto-plastic materials. For Newton-Raphson iteration, only path dependent state determination can be used, and the parameter is ignored. For constant stiffness iteration the parameter takes a value of zero for path independent and one for path dependent state determination. (d) Stiffness reformulation code (KRUSE). If the load step or time step is sufficiently small, then the tangent stiffness matrix may be kept constant over several steps. The value of this parameter can be used to control the frequency of stiffness reformulation. If the parameter is assigned a value of zero, then the tangent stiffness matrix or the effective matrix from the preceding step is retained. If a value equal to n (n 2 1) is specified, the stiffness matrix is updated every n steps. (e) M~x~m~rn number of iteration cycles within a step ~~A~CYC~. This parameter specifies the number of stiffness updates (“cycles” of iteration) permitted during

D. P. MONDKAR and G. H. POWELL

226

(MAXIT).

(k) Limit on displacement increment (DISLIM). It may be desirable to limit increments of displacement during any step of iteration. In most structures. it is possible to identify a principal component of displacement which may be controlled using this parameter. (I) E.recution termination code (IQUIT). If convergence is not obtained within the specified number of iteration “cycles” and iterations per cycle. then it may be desirable to terminate execution to prevent waste of computer time. This parameter is assigned a value of one to terminate execution. or a value of zero to proceed to the next load step or time step.

TOLF.

3.2 Solution schemes By assigning appropriate values to the parameters described in the preceding section. it is-possible to construct a variety of solution schemes. The basic procedures and some of the mixed solution procedures are identified in this section, and values of the parameters for these schemes are suggested in Table 1. The tolerances are rather loosely identified as “small”. “moderate” or “large”. Specific values of these tolerances will be influenced by the particular structure being analyzed. Some solution schemes are as follows: (a) Step-by-step procedure without iteration. with stiffness reformulation every step. (b) Newton-Raphson iteration. (c) Constant stiffness jteration with stiffness reformulation every step. (d) Constant stiffness iteration using initial elastic stiffness throughout. (e) Step-by-step procedure with stiffness reformulation every step, but with constant stiffness iteration in the last load step for static analysis, or at specified time step intervals for dynamic analysis. (f) Automatic stiffness reformulation procedure, in which the stiffness is reformulated only if convergence is not obtained in a specified number of constant stiffness iterations. (g) Mixed iteration procedure in which Newton-

constant stiffness iteration within ‘a single load step or time step. If convergence is not attained within a specified number of iterations, the stiffness is updated, and a new “cycle” of iteration, with a new constant stiffness begins. If the solution does not converge after a specified number of iteration cycles, the solution will typically be terminated, although continuation with the next load step may be specified through use of the parameter IQUIT. (f) Maximum number of iterations per cycle This parameter controls the maximum number of iterations permitted in each iteration cycle. tolerance (TOLF). This convergence (g) fine parameter defines the convergence tolerance in the last load step of any series for static analysis. For dynamic analysis it specifies the convergence tolerance to be used in certain time steps, the frequency of which is controlled by the parameter NZTE The convergence criterion is based on the residual load in any iteration. (h) Coarse convergence tolerance (TOL). This convergence tolerance is used in all steps except the last in static analysis, and in all time steps except those at the specified frequency in dynamic analysis. This tolerance can be specified to be equal to the fine tolerance, TOLF, if desired, but will commonly be larger. (i) StifJess reformulation tolerance (TOLK). In certain cases it may be desirable to perform, in any load step, Newton-Raphson iterations initially, but as the changes in stiffness become progressively smaller, to retain the same stiffness and convert the solution procedure to constant stiffness iteration. With NewtonRaphson iteration, if the residual load vector at the end of any iteration satisfies this tolerance criterion, the previous stiffness will be retained until convergence. (j) Frequency of time steps for fine convergence tolerance (NITF). This parameter controls the time step interval, in dynamic analysis, at which iteration will be performed to satisfy the fine convergence tolerance,

Table I. Control parametersfor some solution schemes ’ Control I Parameter

I

NSTEP

/

I (a) 1 (b) (c) (d) As As 1 As I As ! / appropriate j appropriate( appropriate! appropriate

I I

/

ITYP

I

KPATH

~

KRUSE

!

0 needed Not

!

i

0 needed Not

111

MAXCYC

i

MAXIT

I ,

TOLF

'

TOL

I 1

large

TOLK

1

large

I 1 large

’

11

i'l

0 or

/ 0 or 1

j I /

j

0

l

I /

1

As As appropriate, appropriate

j small 1 small or 1 moderate

I I small I small or ~ moderate

large

I

DISLIH

( 1 1 As i appropriate

IQUIT

(

/1~1;1~0~1j1I

0

1

1

1 = TOLF large As ) As I i appropriate appropriate / As As 1 appropriate1 appropriate II

NITF

!

SOLUTION SCHEME

I

small small or moderate large As appropriate

As

appropriate

I 1 I

/ / 1

(e) j !f) As I AS appropriate: appropriate 1

!

/ Oorl

1

1

i i I : ( 1 :

1

large large

( E NSTEP

0 or 1 0

j

large small

1

say,3

say, 5 I

small

small or 1 moderate large

I

1 appropriate / j

II / 1 / 1

/

needed Nt O

I

1

j I

1 AS appropriate

1 small / small or

I ~ I

moderate

I ,

AS appropriate

AS As appropriate/ appropriate

0

j

moderate AS appropriate AS appropriate

j

Evaluation

of solution

schemes

Raphson iterations are followed by constant stiffness iterations. All of these schemes can be used for either static or dynamic analysis. A limit on the displacement increment in any iteration may be specified, and for schemes with constant stiffness iteration the alpha-constant acceleration scheme may be used. Different convergence tolerances may be specified to obtain results more accurately in some steps than in others. 4. CASE

for nonlinear structures

227

IO

6

STUJMES

4.1 Two-bay

plane truss: large displacemenf elastoplastic static and dynamic response (a) General description. A two-bay plane truss, shown

in Fig. 3, has been analyzed considering both large and small displacement effects. In each case, the truss members were assigned an elasto-plastic constitutive relationship of Ramberg-Osgood type, with the parameters indicated in Fig. 3. For the dynamic analysis, a lumped mass idealization was adopted, with the mass of each truss bar lumped at the bar ends. The transient integration was carried out usmg Newmark’s average constant acceleration operator (fi = l/4. y = l/2, S = 0). (b) Solution schemes. The following solution schemes, from Section 3.2, were selected: (I) Scheme (a) (2) Scheme (c) (3) Scheme (c) with stiffness formulation whenever the structure yielded or unloaded. (c) Static response. Figure 4 shows the static response of the truss with and without geometric nonlinearity. For large displacements the response was computed with a load step of 1 kip. using solution schemes 1 and 2. For small displacements. the response was obtained with a load step of 1 kip and using solution scheme 3. For iteration the convergence tolerance was specified to be 0.05 kip. Convergence was rapid requiring an average of

b

30

I

in

30

I”

I 40

1.

Stress-Strain Law (Ramberg-Osgood Curve)

= 40.52 = i

K = 3:;

KS1

C COj2) lC_? ~n:lr

r=,

= 1SO"C is, = ;?.j' ::“

It, seci,,,4

i :undamentdi Ferlcc

620C .seis

Fig. 3. Plane

truss

;il

04 VERTICAL

Fig. 4. Static

1.2

06

DISPLACEMENT

response

16

20

w(in)

of plane truss

4 iterations per step beyond a load of 6 kips, with almost linear behavior below 6 kips. The number of stiffness formulations in each of the three solution procedures were 11, 11 and 5, respectively, for 11 load steps; whereas the numbers of state determinations were Il.26 and 26, respectively. The computational efforts in each of the three solution schemes can be considered to be moderate and comparable, so that any one of the three schemes may be used for this example. With large displacement effects, the results show very close agreement with the results reported by Noor[91. Also the effect of geometric nonlinearity on the elastoplastic response of the truss is small for the range of loading considered. Table 2 compares the present results (at P = 10 kips) with those of Noor[9] and Goldberg et aL[lO]. Again, the agreement is very close. (d) Dynamic response. The response of the truss to a step load of 10 kips is shown in Fig. 5. The combined effect of geometric and material nonlinearities was included in computing this response. Three different time steps weere selected, namely At, = 62 psec (T,/lOO), At, = 124 Fsec (TJSO), and Ar, = 248 psec (TJ25). where T, equals the fundamental period of the elastic truss. The results given by the solution scheme 1 using each of the three time steps are displayed in Fig. 5 for a duration of 6200 psec. The linear response is also shown, for a time step of 124gsec only. The results of the three nonlinear response analyses are nearly identical. The effect of using equilibrium iterations (solution scheme 2) with a larger time step of 248~sec is also shown. An average of 1.6 constant stiffness iterations per step were required, and the response closely matches those obtained without iteration. Computationally, the three analyses using the step-bystep scheme without iteration required 100, 50 and 25 stiffness formulations. and equal numbers of state determinations; whereas the analysis with iteration required 25 stiffness formulations and 40 state determinations. It is interesting to note that the analysis using iteration and a larger time step of 248/~sec is more efficient compared with the step-by-step analysis without iteration and a smaller time step of 124 psec.

D. P. MONDKAR

(b) I

*

-

i

and G. i-f. P~~LL.

DISPLACEMENT U, INCHES

1.0511

DISPLACEMENT METHOD

1.0574

O"923

**

I

MIXED (FORCEAND DISPLACEMENT) METHOD

40

i

24

E Y fj d

16

NONLINEAR RESPONSE EOUIL. ITERATIONS 41 e 24a (ISECS

NLINEAR RESPONSE EOUIL. ITERATIONS

NONLINEAR RESPONSE NO EQUIL ITERATIONS LINEAR

ii E 08

NONLINEAR RESPONSE NO EQUC. ITERATIONS At 1 248 pSECS

4.2 Cantilever beam: lage displacement elastic static and dynamic response (a) General description. The response of a cantilever

beam subjected to a uniformly distributed load has been computed, considering large displacement effects but assuming the. material to be isa~o~ic, IineatIy elastic. The cantilever dimensions and material properties are given in Fig 6(a). The finite element model of the cantilever is shown in Fig. 6fb). F&e &node plane stress elements were used along the length, with 2x2 Gauss quadrature integration. A lumped mass idealization was used for the dynamic analysis, with one-fourth of the mass of each element lumped at each of the four corner nodes. Newmark’s average constant acceleration operator (B = l/4, y = l/2, 8 = 0) was used for numerical time integration. (b) Solution schemes. The following solution schemes, from Section 3.2, were selected. (1) Scheme (a). (2) Scheme (b). (31 Scheme (cf.

RESPONSE

‘. ‘I

\

(4) Scheme (f), with stiffness reformulation only if convergence was not obtained in five iterations. (c) Static response. The computed tip displacement under static loading is shown in Fig. 7. These results were obtained using solution schemes I and 2. The total pressure of IO~b/~~~was applied in 10. 20. 30 and I00 equal steps for the step-by-step scheme, and in 1, 5 and 20 equal steps for the Newto~~aphson iteration. A convergence tolerance of 0.01 Ib was specified for iteration. Nearly identical results were obtained with 100 steps using the step-by-step scheme and with 20 steps using Newton-Raphson iteration. and both resuits compare well with the analytical results given by Holden[ll]. With the load applied in 1 and 5 steps. solution scheme 2 gives results which are almost the same as those obtained with 20 steps. With a large load step (10 equal steps) and solution scheme I’(step-by-step procedure) the results show cons~de~ble deviation from the analyticat resufts of Hotden. The results using smaller load steps (i.e. 20 and 30 equal steps) show improvements the

Evaluation of solution schemes for nonlinear structures

fl P/2

in

lb/in* L

I-

I

(01 BEAM

PROPERTIES

P/2

(b)

FINITE

ELEMENT

Fig. 6. Cantilever

P/4

MODEL

beam.

229

These analyses indicate a need to use small load steps with the step-by-step solution scheme. For this stiffening structure the Newton-Raphson iteration scheme performed well for all load steps. (d) Dynamic response. The dynamic response of the beam subjected to a step load of 2.85 lb/in’ is shown in Fig. 8 for a duration of 0.0135 sec. Two different time steps were considered, namely AI, = 45 psec (TJl20) and At, = 135psec (T,,/40), where T0 is the fundamental period of the elastic cantilever. The results shown in Fig. 8 were obtained by solution schemes 1 and 3 using both the time steps in each case, and by scheme 4 using a time step At,. A convergence tolerance of 0.1 lb was specified for iteration. The response given by scheme 3 with the short time step AI, will be assumed to be “exact” for the purpose of comparison. The linear response (considering only small displacement effects) obtained using a time step At, is also shown. As is to be expected, the nonlinear response is stiffer than the linear response. With a time step of 45 psec, the results obtained with and without iteration are very close, indicating that the equilibrium iterations can be omitted with such a small time step. With a time step of 135gsec, the results with iteration are very close to the “exact” results, whereas those without iteration are grossly in error. This emphasizes the need for equilibrium iterations when larger time steps are used. The automatic stiffness reformulation with iteration (solution scheme 4) using a time step Af, gives results very close to the “exact” response. These analyses indicate that a time step equal to approximately T,,/40 is needed to obtain accurate response for this example. The computational efficiency of the various analyses is compared in Table 3. It can be seen that among the procedures giving accurate results, scheme 4 (automatic stiffness reformulation) is computationally most efficient for this example. This efficiency can be attributed to the fact that ihe stiffness is reformulated only when necessary, rather than arbitrarily at each time step. The exact response compares reasonably well with the results reported by Bathe et af.[12], with maximum difference not exceeding 5%. 4.3 Clamped beam: large displacement elastic static and dynamic response (a) General description. A clamped beam subjected to

TIP

DISPLACEMENT

Fig. 7. Large displacement

RATIO

w/L

elastic static response of cantilever beam.

results with 30 load steps being accurate for practical purposes. Computationally. the step-by-step scheme required IO. 20. 30 and 100 state determinations, and equal numbers of stiffness formulations for 10,20,30 and 100 loading steps. respectively; whereas the NewtonRaphson iteration scheme involved 8, 25 and 71 state determinations and equal numbers of stiffness formulations for I. 5 and 20 loading steps, respectively. The Newton-Raphson iteration required. therefore, an average of 8. 5 and 3.5 iterations per step. respectively.

a central concentrated load has been analyzed. Only large displacement effects were considered, and the material was assumed to be isotropic, linearly elastic. The geometry of the beam and the material properties are shown in Fig. 9. The finite element model of the beam consisted of five I-node plane stress elements in the half span, with 2 x 2 Gauss quadrature integration. For dynamic analysis, a lumped mass idealization was used and the time integration was performed with Newmark’s operator (/3 = l/4, y = l/2, S = 0). The dynamic response of this beam has also been studied by McNamara[l3], and a similar problem was solved by Weeks[l4]. In the studies by McNamara the beam half span was modelled using five beam bending elements, whereas Weeks used a simpler model with approximate nonlinear effects. (b) Solution schemes. The following solution schemes. from Section 3.2. were selected. (I) Scheme (a). (2) Scheme (b). (3) Scheme (f). with stiffness reformulation only if

D. P. MOSDKAR and G. H. POWELL

230

0.8

NONLINEAR

RESPONSE MQ EQUlL ITER~PTIOHS AI:OM IO+SEC

NONLINEARRESPONSE WlTn EQ",L ITERATIORS At 2 135 SEC ALSO FORUUL4TlOM

At

Fig. 8. Large displacement Table 3. Cantilever

(1) with at = At,

!

~am~ompar~son

100

(3) with bt = At2

100

with

I w?th At = At2 ~‘Aut~tic” Stiffness'

fotnulation

*

71.

etlkiency

No. of State Determinations 300 I

100

Average No. of Iterations/Step

/ 1

776 529

No iteration Vo iteration

7.8

I

/

!

5.3

Procedure gives inaccurate response

convergence was not obtained in five iterations. (4) Scheme (c). (c) Sfofic ~~p~#~~. The variation of central displacement with load is shown in Fig, 9 for both Linear and nonlinear behavior. It can he observed that the beam exhibits a highly stiffening behavior, the linear displacement being several times larger than the nonlinear displacement. The total load of 700 lb was incremented using two different stepping schemes, namely (1) unequaS stepping scheme in which the load was applied in steps of 10ib up to 100lb, and then in steps of 100lb up to a total of 700 lb and (2) equal stepping scheme in which the load was applied in seven equal steps of IOOIb each. Using the unequal stepping scheme, the response was computed by both the step-by-step procedure (solution scheme 1) and Newton-Raphson iteration (solution scheme 2). The results of the two analyses were nearly identical for the entire range of Ioading, and are shown in Fig. 9 (the results with the step-by-step procedure are plotted only up to a load of IOOIb). Newton-Raphson iteration required 47 state determinations and an equal number

of computational

300

bt = At2

*(l)

: 13sIO-‘SEC

elastic dynamic responseof cantilever beam.

No. of Stiffness Formulations

SOLUTION PROCEDURE

IO-’

“&UTObtATIC” STIFF WIT” ITER6ilONS

of stiffness

formulations.

It involved

an average

of 3 iterations per step to obtain convergence to a tolerance of I fb. On the other hand, the step-by-step procedure (solution scheme if required only 16 stiffness

formulations and an equal number of state determinations. Using the equal load stepping scheme. the displacement was calcufated by Newton-Raphson iteration with a Iimit of 0.5 in, imposed on the displacement increment in any iteration. The results. shown in Fig. 9, are almost the same as those obtained with the unequal load stepping. An.average of 4 iterations per step, and a total of 27 stiffness formulations (with an equal number of state determination~~ were required. Newton-Raphson iteration without a displacement limit failed to converge within 5 iterations in the first load step. The step-by-step procedure (solution scheme 1) with equal load steps gave totally incorrect results, which are not shown, because of the highly nonlinear load-displacement response of the beam in the first load step. These analyses indicate that a small load step is needed to obtain accurate results in the initial stiffening part, and that a larger step may be used subsequently. The displacement limit in the Newton-Raphson scheme has the effect of reducing the load increment in any iteration. Solution scheme 1 (step-by-step procedure) with unequal load steps, and solution scheme 2 (NewtonRaphson iteration) with a displacement limit and equal load steps performed best for this example. (d) l)ynamic response. The response of the beam to a

Evaluation of solution schemes for nonlinear structures

0

0.1

0.2

0.3

0.4

0.5

DISPLACEMENT, Y (INCHES)

Fig. 9. Large displacement elastic static response of clamped beam.

dynamic step load of 640 lb was studied for a duration of XKM~sec, Because of the highly stiffening behavior at a static load of 44OIb, it can be expected that the beam subjected to a dynamic load will vibrate with a period considerably shorter than the linear period of vibration, and that this will influence the choice of time step. Figure 10 shows a comparison of the linear and nonlinear responses of the beam. The linear responses were obtained with two different time steps, namely 50 psec and lOOpsee. The nonlinear response was computed with a time step of 50 gsec, using solution scheme 1. The linear responses obtained with the two time steps are almost identical, and the period of vibration compares well with the value of To = 9056 Fsec obtained using a beam formula[I5]. On the other hand, the period of the first cycle of nonlinear vibration is seen to be approximately 2300 bsec. The considerable difference in the maximum displacements of the linear and nonlinear solutions can also be noted.

231

Figure 11 gives the non~in~ responses of the beam, computed using various time steps and solution schemes. Three different time steps, namely Ai, = 100 &sec, At, = 50 *see and Ar, = 25psec were used, corresponding to ratios of about T&JO, T&80 and Td360, respectively. With each of these time steps, the nonlinear analyses were carried out using the step-by-step procedure (solution scheme 1). The analysis with the smallest time step, A&, will be assumed to be “exact” for comparison purposes. The response was also computed using the automatic stiffness reformulation procedure (solution scheme 3) and time steps of At, and At,. The convergence tolerance for iteration was specified to be 10% of the applied load. The results of the step-by-step analysis with time step At, are reasonably close to the “exact” results; whereas those with time step At, show considerable differences, indicating a need for equilibrium iterations to obtain accurate response. The results obtained with the automatic stiffness reformulation (solution scheme 3) using time steps At, and AI, both compare well with the “exact” response. Computationally, the step-by-step analyses with time steps, Al,, Ahr,and At, required 50, 100 and 200 stiffness formulations, respectively, and the same numbers of state determinations. The analyses with the automatic stiffness reformuIation scheme required 45 stiffness formulations and 489 state determi~tions for a time step of 50 psec and 40 stiffness fo~ulations and 297 state determinations for a time step of 1OO~ec; that is, these analyses involved an average of 4.9 and 6 iterations per step for the two time steps, respectively. The analysis of the beam was repeated using solution scheme 4 and a time step of 1OOpsec. However, the iteration failed to converge in the second time step within a specified limit of 10 iterations and a convergence tolerance of 1046 of the applied load. From computational standpoint it can be seen that the automatic stiffness reformulation procedure with iteration (solution scheme 3) worked better among solution schemes giving accurate results. McNamara[ 131 analyzed the beam using a central difference operator with a time step of 5 met, and obtained a maximum displacement of 0.9Oin. compared to 0.77 in. in the present analysis. The period of the first cycle of nonlinear vibration as reported by McNamara was 2883 psec, compared to a period of about 2300 psec

I I

l LINEAR

At

1000

2000 3000 TIME ( pSEC 1

RESPDN6E * 100 p6EC6

4000

5000

Fig. IO. Comparison of linear and nonlinear dynamic responses for clamped beam.

232

D. P. MONDKAR and G. H. POWELL I

000 .

g

0.675

NOULINER RESPONSE “AUTOMATIC“ STIFFNESS FORMULATION WlTH ITERATIONS ~50 ~SECS.

A+

5 5

0.750

l c

0.625

z z !: 0.500 ii y 0.375 0

i

x

0.250

NONLLHEAR RESPONSE “AUTOYATlC* STlfF,,ESS FORMULATION WIT!+

Fig. II. Large displacement elastic dynamic response of clamped beam. P

obtained by the present analysis. These discrepancies in the two investigations can be attributed to the fact that McNamara used a beam bending element which is more flexible than the two dimensional plane stress element used in this analysis.

(a) General dewiprion. The behavior of a spheticai cap under concentrated apex load has been studied for both static and dynamic loads. The geometry of the shell and the material properties are shown in Fig. 12. Only geometric nonlinearity was considered, and the material was assumed to be isotropic, linearly elastic. The finite element model of the cap consisted of ten 8-node axisymmetric elements, with 2 x 2 Gauss quadrature integration. For dynamic analysis a lumped mass ideahzation was used, and the transient integration was carried out using Newmark’s average constant acceleration scheme (j3 = l/4, y = 112,6 = 0). (b) So~~f~onschemes. The following solution schemes, from Section 3-2, were selected. (1) Scheme (a). (2) Scheme (b). (3) Scheme (c). (4) Scheme (f), with stiffness formulation only if convergence was not obtained in five iterations. (c) Static response. The apex displacement response of the cap is shown in Fig. 12. The results were computed using solution schemes 1 and 2, up to a total load of 1001b. For solution scheme I, the load was incremented using steps of i lb and 5 lb, and for solution scheme 2 the load was incremented using steps of 5 Ib and 20 lb. The cap behavior (shell parameter A = 6 where A’= t2( 1- ~z)~4~~2~2)as shown in Fig. 12 is highly nonlinear, with initial softening and subsequent stiffening. With a load step of 1 lb, the step-by-step procedure (solution scheme 1) yields results which closely agree with those obtained by Newton-Raphson iteration with a 5 lb load step. Newton-Raphson iteration with steps of 201b also gives nearly identical results (these results are not plotted in Fig. 12). However, with a load step of 5 lb the step-by-step scheme gives response showing considerabte drift from that obtained with a 1 lb step, except in the stiffening region beyond a load of about 70 lb. This

R=476in -0soin

0

HI

006569tn

t = 001576 E = ioooo

VI t*i

Y -03 p = 2.4Sx10*41b-,~cZ/ln*

0

r

/

/

0.4

0.8

8.2

DISPLACEMENT

Fig. 12. Large displacement

t .6

RATIO

2.0

r/H

elastic static response of shallow

cap.

indicates a need to use a small load increment in the initial softening region if the step-by-step procedure is used. Because the cap behavior is almost linear in the stiffening region, larger load steps may be taken subsequently. ~ompu~tion~y, the step-by-step procedure required 100 and 20 stiffness formuIations for load steps of 1lb and 5 lb, respectivety, whereas Newton-Raphson iteration involved 83 stiifness formulations with an average of 4.2 iterations per step for 5 lb toad increment, and 36 stiffness formulations with an average of 7.2 iteration per step for 20lb load increment. The Newton-Raphson iteration scheme therefore performed better for this example in terms of accuracy and efficiency. Figure 12 shows a comparison of these analyses with the numerical results of Haisler et ol.[ 161, and the experimental studies of Evan-Iwanovski et a1.1171.The comparison with the results of Haisler is very close, whereas the experimental results show considerable

Evaluation of solution schemesfor nonlinear structures

233 LINEAR

RESPONSE

At : z pacts

_

At : 4 @ECS

l

tKXNClWCAR

___

24

At

‘i 0:

‘@

‘8.

AI

O 0

Q

r

:\

0

NO

_ A ‘\

16

$ $

”

0

I 50

I 100

I

I

I

I

150

200

250

350

TIME

U

ITERATIONS. ~SECS

EOUIL.

ITERATIOWS.

“AUTOUATIC” STIFCNESS FORYULATiON WITH 1TERATION:A~ : 4 pSECS

0

9

1.2

: z

ITERATiONS.

~SECS

At = 4 ~~SFCS

D

: ‘+ 5

: 2

_._ NOEOUIL

$.@,?O

t

RESPONSE

W1WEOUtL

0

. I

I

350

400

I 450

I 500

QiSECSf

Fig. 13. Large displacementelastic dynamic response of shallow cap. differences. However, it should he noted that the experimental cap model had a shell parameter A = 6.23. The present analysis also agrees fairly closely with the finite difference solution of Mescallll8]. (d) ~nu~i~ response. The response of the cap sub jetted to a dynamic step load of 100 lb is shown in Fig, 13 for a duration of 500 psec. The linear and nonlinear solutions are both shown. The response was computed using two different time steps, namely At, = 2ksec (approximately r&O) and At, = 4psec (approximately Ta/30), where r, is the fun~en~l period of the elastic cap. The linear response was calculated using both time steps. For nonlinear analysis, the results shown in Fig. 13 were obtained by solution scheme 1 using both time steps, by solution scheme 3 using a time step AtI and by solution scheme 4 using a time step A&. A convergence tolerance of 10% of the applied load (i.e. IOlb) was specified for the iteration. The results of the analysis using constant stiffness iteration (solution scheme 3) with a time step of 2psec are assumed to be “exact” for comparison purposes. The linear responses obtained using the two,time steps are practically the same, except that some differences can be observed after the third cycle of vibration. With a time step of 2 gsec, the results of the nonlinear analysis without iteration differ little from those with iteration (“exact” response). A need to perform equilibrium iterations when a larger time step is used is demonstrated in Fig. 13. in which it is seen that with a time step of 4 gsec the computed response without iteration (solution scheme 1) differs considerably from the “exact” response: whereas the results using the automatic stiffness reformulation with iteration (solution scheme 4) are close to the “exact” resuits. The analysis was repeated using solution scheme 3 (constant stiffness iteration) and a time step of 4 gsec. but the iteration failed to converge in the second time step within a specified limit of 10 iterations per step. The step-by-step scheme involved 250 and 125

stiffness formulations (and equal numbers of state determinations) for time steps of 2~sec and 4psec, respectively. The constant stiffness iteration (solution scheme 3) required 250 stiffness formulations and 1025 state determinations (i.e. an average of 4.1 iterations per step) for a time step of Zibsec, The analysis using automatic stiffness reformulation scheme required only 117 stiffness formulations, but 735 state determinations (i.e. an average of 5.9 iterations per step) for a time step of 4~sec. For this example, the automatic stiffness reform~atio~ scheme proved to be ~ompu~~o~lIy attractive. It is also interesting to note that in this scheme the stiffness was reformulated nearly each step, but that plain constant stiffness iteration failed to converge. This indicates the advantage of reformulating the stiffness only when necessary rather than automatically at each time step. The “exact” response shown in Fig. 13 compares closely (within 5% error) with Lheresults reported in [ 121, in which the analysis was carried out using the Wilson-B methodi and a time step of approximately T&O. 4.5 Shallow arch: large displacement elastic static response

(a) Ge~e~af desc~pt~o~. The large displacement response of a shallow elastic arch subjected to a concentrated apex load has been studied. The geometry of the arch and the material properties are indicated in Fig. 14. The material is assumed to be elastic and isotropic. The finite element model of the arch consists of eight &node plane stress elements between the fixed end and the apex (i.e. symmetry is taken into account), with 2 x 2 Gauss quadrature integration. (b) Solution schemes. The following solution schemes, from Section 3.2, were selected. (1) Scheme (a). (2) Scheme (c). (c) Static response. Figure 14 shows the static. response of the arch. The analyses were carried out with the step-by-step procedure using load steps of 1lb and 5 lb,

D. P. MONDKAR and G. H. POWELL

234

scheme) diverged considerably. Comp~tat~onally. the constant stiffness solution scheme required an average of 3.8 and 6.0 iterations per step for 1lb and 2.5 lb, load steps, respectively. The results of the present study are also compared in Fig. 14 with the numerical results of Mallet et al.[20] and with the experimental results of Cijelsvik ef al.[21].

R: 133.1138Inches 8: 73373 dsgrats t = 0.188 Inch Q = 170 Inches H = 1.09 Inchas

.

4.6 Plane strain punch problem: small displacement elasto-plastic static response (a) Generai description. The behavior of a solid plane strain specimen subjected to indentation by a rigid punch has been studied. The materiai of the specimen was assumed to be elastic-perfectly plastic and only small displacement effects were included. The specimen dimensions and mate~a~ properties are given in Fig. 15. For the finite element model, higher order @-node) plane strain elements with 2 x 2 Gauss point integration were used in the region of the specimen directly under the punch, where yielding of the material is expected, whereas lower order (7-node and 4-node) elements were used in regions of the specimen away from the punch. Two different types of analysis were carried out, namely an imposed load analysis and an imposed displacement analysis. For the imposed toad analysis, the punch was modelled as a number of finite elements with a very large value of Young’s modulus. For the imposed displacement analysis, displacements were imposed through very stiff ve~ic~.sp~ngs. In both cases friction at the interface was assumed to be zero. (b) Solution schemes. The following solution schemes, from Section 3.2, were selected. (1) Scheme (c). (2) Scheme (b). (c) Static response. Figure 15 shows a plot of mean pressure vs punch indentation (both quantities are nondimensionalized). For the imposed displacements case, the total punch indentation of 0.008 was applied in three ways, namely (1) 16 equai steps of 0.0005 each, (2) 4

CORMEO EACH STEP PRESENT YIUDY. LOAD STEP I ID COI)nECTIoN, NO ITERATION LOAD STEP I lb‘.

EQUU..

A

IP

Of

i

Q

PRESENT &tUOY. EOUIL. CORRIICTION. NO ITERATION PlESENT STUDY, LOiOSTEP 2.3 lb. CONSTANT ST,FF.

ITE”ATIOY

I

I

0.1

0.2

DlSPLACEMENf

0.3 v

0.4

0.5

(tNCHES)

Fig. 14. Large displace~enF elastic static response of shallow arch.

and with constant stiffness iteration using load steps of I lb, 2.5 lb and 5 lb. With a load step of 1 lb, the result without iteration and with iteration (convergence tolerance 0.1 lb) were almost identical. With a load step of 2.5 lb and iteration (tolerance 0.25 lb) the results were again very close. However, with 5 lb load step the iteration procedure failed to converge within 15 iterations and a tolerance of 0.5 lb beyond a load of 2Olb. whereas for the same load step the results with solution scheme I (step-by-step

lYPOSED DISPLACEHEHTS CDNST. STIFF. ITERATION, STIFF. REFORM. EACH STEP 16 EPUAL STEPS, 0~0.0003~0.008 0 SINGLE STEP OF 0.008 N-R ITERATION 4 EOUAL STEPS, 0~0.002)0.008 A PATH-DEPENDENT STATE DETERMINATION IMPOSE0 LOADS 0 N-R ITERATION, 14 STEPS, 0(0.2)0.8(0.04)1.2

2.0 -

1.6 -

.aP 0

/ 1.2 NONLINEAR

t

2b \

\--_------/

I

* F

h

2L

RESPONSE

L/b 92.7 h/b = 1.7 ELASTO-PLASTIC MATERIAL E = 10,000 kri Y = 0.33 *Y =13.0 ksl PO =2 UY /Jj

~~, _....A_

.ooe

PUNCH t~OEt+tAtt0~

RATIO w/b

Fig. 15. Static response of plane strain punch problem.

Evaluation

of solution schemes for nonlinear structures

posed displacements case. In both analyses, collapse was obtained at a pressure of 1.2. Essentially identical results were obtained by Nayak and ZienkiewiczI221, and the collapse load closely agrees with that predicted by a slip-line solution of the probleml23). The yielded integration points for the imposed load case at a pressure of 1.2 are shown in the left half of Fig. 16. Yielded zones are shown for the imposed displacement case in the right half of this figure, at displacements of 0.002, 0.005 and 0.008, respectively.

equal steps of 0.002 each and (3) a single step of 0.008. in each case the response was computed using ~oos~nt stiffness iteration (solution scheme I) with path independent state determination. The results of the analysis with 16 equal steps are shown in Fig. IS. Nearly identical results were obtained with 4 equal steps (these results are not plotted). With the displacement applied in a single step very similar results were also obtained. The analysis with I6 steps required I3 stiffness formulations, 24 state determinations and an average of 1.6 iterations per step with a convergence tolerance of 1% of the yield stress. The corresponding numbers were 4,26 and 6.5 for 4 equal steps, and I,25 and 25 for one step of imposed inden~tion. A larger convergence tolerance, equal to 5% of the yield stress, was specified for the single step case. The analysis was also carried out using NewtonRaphson iteration (solution scheme 2) and 4 equal steps. In this analysis path dependent state determination was used. However, the results obtained were totally incorrect, as shown in Fig. 15. Also, at the center of the most severely strained element a horizontal strain of 0.00915 and an effective stress of 8.47 were computed at the imposed displacement of 0.008, compared with a strain of 0.00979 and an effective stress of 13.0 (indicating yield) obtained using constant stiffness iteration with path independent state determination. This illustrates the need to use path inde~ndent state determination for imposed displacement cases, as was pointed out in Section 2.5. For the imposed load computations, the total pressure of 1.2 was applied in 4 equal steps of 0.20 each followed by 10 equal steps of 0.04 each. The response was computed using Newton-Raphson iteration with path dependent state determination. The results show very close agreement with 16 steps of imposed displacement This analysis required 18 stiffness formulations, 25 state determinations and an average of 2.6 iterations per step beyond a pressure of about 0.92, up to which the response is seen to be linear. Because the strains increase progressively for imposed loads, Newton-Raphson iteration with a path dependent state determination yielded correct results, in contrast to the totally incorrect results obtained with this solution scheme for the imLINE RIGID

PLANE

OF

235

5. COMXUDiNG REMARKS The principal intent of this paper has been to discuss a solution strategy implements in a general purpose finite element computer code for nonlinear structural analysis. The use of a flexible solution strategy such as this is believed to be desirable because it is unlikely that a single solution scheme provides optimal computational efficiency and accuracy for all types of nonlinear behavior. The numerical experiments described in this paper cover a variety of nonlinear structures, and confirm the need for such a solution strategy. The results presented in this paper show that it is necessary to perform equilibrium iterations to obtain accurate results when large or moderately large load steps or time steps are used. However, it is generally useful to compute the response using the step-by-step scheme without iteration when the nature and degree of the nonlinear behavior is not known d priori. More sophisticated solution schemes can then be used to refine the initial analysis. For the dynamic analysis of geometrically nonlinear structures, in which a continuous change in stiffness occurs, the automatic stiffness reformulation scheme with iteration has been shown to be both accurate and computationally efficient. The advantage of this scheme is that it updates the stiffness matrix only when necessary, rather than at every step or at specified time step intervals. For static analyses of elasto-plastic or geome~cally nonlinear s~ctures, Newton-~phson iteration always performed well, even for large load steps. However, with

SYMMETRY

PUNCH

-

w/b

= 0.002

____

w/b

= 0.005

--

w/b = 0 008

STRAIN

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\~\\\\\\\\\\\\\\\\\\\\\\\\ (a)

IMPOSED

LOADS,

p/p,

= 1.2

!

(b)

IMPOSED

DISPLACEMENTS

Fig. 16. Plastic zone profiles for plane strain punch problem

D. P. MONDKAR and G. H. POWELL

236

a stiffening load displacement response it may be necessary to use small steps in the low load region, or to specify limits on the displacement increment in any iteration. The step-by-step scheme has been shown to require small load steps, particularly if there is an initial stiffening region The study of the punch indentation problem has shown the need to use consent stiffness iteration with path independent state detestation in the case of imposed displacements, and that path dependent state determination may lead to erroneous results.

REFERENCES

I. D. P. Mondkar and G. H. Powell, Static and dynamic analp sis of nonlinear structures. Report No. EERC 75-10, Earthquake Engineering Research Center, University of CaIifornia, Berkeley (1975). 2. P. K. Larsen. Large displacement analysis of shells of revolution including creep, phsticity and viscoelasticity.Report No. SESM 71-22. ~u~rnent of Civil Enaineerina, University of C~~o~~, &rkefey (1971). 3. J. T. Gden, Finite ~ft~~~~ of ~oal~ear ~aat~au. McGrawHi& New York (1972). 4. H. D. Hibbitt, P. V. Marcal and J. R. Rice. A finite element formulation for problems of large strain and large displacement. Inf. I. Solids Strucr. 6, l06!%1@6 (1970). 5. N. M. Newmark, A method of computation for structural dynamics. 1. Engng Mech. Div. AXE 85(EM3),67-94 (1959). 6. D. P. Mondkar and G. H. Powell, Towards optimal in-core equation solving. hr. I: Compuf. !Wcr. 4.531-548 (1974). 7. G. H. Powell, Finite element analysis of elasto-plastic tee

joints. Report No. UC SEW 74-14, Department of Civil Engineering, University of California, Berkeley (1974). 8. G. C. Nayak and 0. C. Zienkiewicz. Note on the alphaconstant stiftness method for the analysis of nonlinear problems. int. 1. for NWW. M&h. Eagttg 4,5?9-582 (1972).

9. A. K. Noor, Nonline~ aaalysis of space trusses. I. Srrucf. Div,, ASCE Hl&ST3),533-546 (1974).

IO. J. E. Goldberg and R. M. Richard. Analysis of nonhnear structures. J. Struct. Dir.. ASCE 89lST41.333-35I (19631. 11. J. T. Holden, On the finite deflections of thin beams. Irrr. j. Solids Sfrucr. 8, IOSI-1055(1972). 12. K. J. Bathe, H, Gzdemir and E. L. Wilson. Static and dynamic geometric and material nonlinear analysis. SESXI Report No. 74-4. Department of Civil Engineering. L’niversity of California, Berkeley (1974). 13. J. F. McNamara, Solution schemes for problems of nonlinear structural dynamics, J. Pressure Vessel Tech.. .-1S&3E96 96102 (1974). 14. G. Weeks, Temporal operators for nonlinear structural dynamics. J. Engng Meek. Dir., ASCE 98(EM41. 1087-I104 (1972). IS. W. C. Hurty and hf. F. Rubinstein. Dynamics of Srrwfures. Prentice Hall, Englewood Cliffs, New Jersey (19671. 16. W. E. Haisler. J. A. Stricklin and F. J. Stebbins. Drvelooment and evaluation of solution procedures for geometricali!, nonlinear structural analysis. AIAA J. lo(J). X4-272 (1971). 17. R. M. Evan-Iwanovski. H. S. Cheng and T. C. Loo. Experimental investigation of deformations and stability of spherical shells subjected to concentrated load at the apex. Proc. Foutih U.S. ~atioaaf Congress of ‘~~~~ied j~erkanics. ASME pp. 563-575 (1952).

18. J. F. Mescall. Large deflections of spherical shells under concentrated loads, 1. A&. Meek. 32.936-938 (19651. 19. K, J. Bathe and E. L. Wilson, Stability and accuracy analysis of direct integration methods. Inr. J. Earthquake Engng nnd SMIK~.Dynamics. 1, 283-291 (1973). 20. R. H. Mallet and L. Berke, Automated method for the large deflection and instability analysis of three dimensional truss and frame assemblies, AFFDL-TR-66-102(1966). 21. A. Gjelsvik and S. R. Bodner. The energy criterion and snap buckling of arches. J. Engng Meek. Die.. ASCE 88fEMS). 87-134 (1%2). 22. G. C. Nayak and 0. C. Zienkiewicz, Elasto-plastic stress analysis. A generalization for various constitutive relations including strain softening. Int. J. /or Namer. Meth. Engng 5. 113-13511972). 23. R. Hi& Tfre ~a~~~~utic~i Theory of P~ast~cit~~ Clarendon Press. Oxford (19.50).