An energy balancing strategy for solution of combined geometrical and material nonlinearity problems

AN ENERGY BALANClNG STRATEGY FOR SOLUTION COMBINED GEOMETRICAL AND MATERIAL NONLINEARITY PROBLEMS RUDOLPH

Institute

of Structural

Mechanics.

OF

SZIL.ARl>+

Universiry of Dortmund. Germany

Dortmund.

Federal Republic of

Abstract-This paper describes an energy based line-search technique for large displacement and nonlinear stability analysis of planar truss and frame structures which. in addition IO geometric oonlinearities. s~mult~lneousiy exhibit nonlinear material behavior. To improve the convergence characteristics of incremental-iterative procedures. a line search seeks a scalar factor which. at fixed load level. scales the displacement vector in such a way that the fundamental energy identity is satisfied. Furthermore. this factor can also be used for guiding the load incrementation. Aiming at computational simplicity and economy. efficient “section-type” beam elements are introduced. 6sing the in/ Loprcrrrgitrn formulation, consideration of geometric nonlinearities yields. however, stiffness matrices which are linear and quadratic functions of the displacements. The material nonlinearity is treated b) bilinear laws. The inelastic behavior is defined for several cross sections: thus. by forming the element stiffness matrices for bending the resulting variable stiffness must be accounted for. The effectiveness of the simple solution procedure has been tested using planar truss. beam. and cable-beam-type structures with known analytical and/or alternative solutions. Good agreement with these existing solutions demonstrates the applicability of the proposed method for routine engineering analysis.

traduced which considers the interaction between axial force and bending moment. Due to the simultaneously occurring geometrical nontinearity, the exact expression in the moment-curvature relationship has been employed in determining the reduced bending stiffness properties for several cross sections of the element. By formulating the element stiffness matrices for bending, the virtual work technique coupled with numerical integration is utilized. The physical law adapted is of bilinear type with strain-hardening portions. While such hyperelastic material can exhibit different uniaxial behavior for tension and compression. consideration of permanent plastic deformations in course of unloading is excluded from the present treatment. in the line-search part of the approach. the pertinent constitutive equations are numerically integrated. Since the present method is limited to conservative loads and to systems which. while undergoing large deflections exhibit moderately large rotations, an extension for elimination of these restrictions is also outlined.

I. lNTRODUCTION

The past decade has witnessed considerable progress in the development of computational methods for geometrically nonlinear problems with linear elastic materials. A real need still exists. however, to analyse structural configurations in which geometrical and material nonlinearities jointly occur. An additional important requirement of the engineering praxis is computational simplicity and economy. The present study aims at freedom from complexity in formulation and at efficiency in solution for a broad spectrum of planar bar and beamtype structural problems exhibiting combined nonlinearities. The powerful-yet straightforward-energy balancing techniques introduced by the author for large deflection analysis1 I ] and. to trace the nonlinear pre- and postbuckling behavior of elastic structuresj?]. has been extended herein to include nonlinear material behavior. The solution procedure involves piecewise linearization at each incremental step (given careful consideration to step-increment strategy). followed by a simple energy based line search which accelerates the convergence of the secant-Newton equilibrium iteration. This proposed incremental-iterative method employs the total Lagrangian formulation. The associated computational strategies require both the tangential and secant material descriptions. For treating material nonlinearities a computationally efficient section-type model has been in-

2.

PRELIMINARY

COSSIDERATIOSS

The nonlinearities in structural analysis arise from two distinct sources which are: geometric and constitutive. Since the primary aim of this paper is the development of computationally efficient algorithms considering both types of nonlinearities. simple formulation in the solution strategy as well as in the element modelling is considered to be mandatory. Consequently. in the present paper we employ the incremental-iterative approach using the

* Professor and Chairman. l-17

I L,?

R. SLlt

finite

element

generality. grangian wideI), tively.

since

frame

order erally

energ)’

equation.

additional

to account directly.

The

necessity

are that

of history

using

terms.

limitations

3.

care since

modelling

remain

(2) Although

plane

and normal

sulting (3) The

strains

are small

bending

curvature

moment

straight

compared

in considerable

terial

nonlinearities.

the achievable

stiffness

for

accuracy[3].

to the

slopes

(Fig.

the

stress

that the history

mensional

restrict

truss

conservative

nonlinear [I!].

As

and frame

static

pendent

pressure

present

treatment.

our

loads. loads

discussion

to

structures Thus. are

incre-

= ,K;::‘_,

(1)

displacement

tangent

modulus

e. The

vector

of elasticity

subscripts

1~ and IJI

(0 = start.

for

the geometrically

expressions Kc:(d)

bar

and

beam

and Kz:

indicated.

these are linear

of the displacements.

(d’)

=

stiffin Ref.

nonlinear

and quadratic

functions

respectively.

the first

displacement

f

element

are given

geometrically

3. I .Z Firs1 r.rfimtrrr,. Using strategy.

in

counter

vector

Euler’s

estimate

simple

for-

of the incremental

is obtained

from

is re-

, = f ,K:II:‘_, )- ‘.Ap ,,I_ I.

A dj):‘_

(2)

has no drift

two-di-

subjected

deformation

excluded

current

matrices

Furthermore.

of the loading

the

pin-jointed

which.

we

load

of the system

+ K$‘L]j,‘,’

is the iteration

Explicit

ward

the ac-

material

law.

parentheses

ma-

affecting

reversal.

Finally.

is ;I

some form

I refer to the load level. while the superscript

stiffness

I) re-

in treating

noticeably of

ma-

diagrams

Consequently.

diagram

matrix

+ K,“.

using

ness matrices

structural

simplification

by such a bilinear

it is assumed

that

used here.

law includes

the end of lrlth

rE,,‘,.‘,. for each element

final).

to unity.

of the stress-strain

without

stress-strain

the re-

law).

lines of different

sults

stress

are large.

is proportional

has shown

representation

placed

At

the tangent

dj,/,,‘,.and the pertinent +

(Bernoulli-Euler

Experience

tual

(Ber-

hypothesis). the displacements

is

balancing.

nonlinear

to the axis of the beam after deformation noulli-Navier

technique

I, I L’ptlirri~g.

ment

is updated

sections

energy

on the

are:

( I) The normal

of the procedure

dependence.

TK(),{’ = [KL

beam elements

equi-

1). In ad-

the

special

solution. in

stability

(Fig.

and gen-

In

“small”

made

of the

to restore

reduced

solution

an

circumstances.


bvhen the material

Iekel.

than in the case

formulation.

in neglecting

hypotheses

load

the conservation

bp the global

incremental

higher

for

of this approach

may put certain

stiff-

of iteruti\e

a given

these

enhanced

should

by two

Cnder

the numerical

holyever.

The

At

to satisfy

further

in the Lagran-

formulation.

terials

seeks

dition.

respec-

analyses.

load steps can be taken

elimination

the loads.

can be considerably

total Lagrangian their

search

librium

can be solved

be taken

energy number

however.

Lagrangian

incrementing

La-

are by far the mo3t

Advantages

larger

of updated

and

hypotheses.

required

the displacements

total

the assumption

geometry

of reference.

effects.

\kith

structural

are

the

and

the undeformed

ness matrices

adopt

the>

formulation

in nonlinear

Using

we coupled

strains.

used

of its simplicit),

because

formulation.

of Green’s

gian

method. Furthermore.

\KI)

from

to

due to the linearization

from

response

process

the true equilibrium (Fig.

de-

3. I .3 A,I

the

the approximate of the total

applied.

position

will

for nonlinear

2).

rr~o;qx

stwrch. Because

hc/.sc~l lick

nature

of

of eqn (2). the first estimate

displacements

, = dj,‘,’ + 2, d!::‘_ ,

d);;‘_

(3)

3. SOI.L’-l’ION PROCEDL’RES

The

solution

nearization

procedures

of the nonlinear

involve

piecewise

li-

equation

by

response

pertinent

to the current

the fundamental load

increment

tween

first

then

in such

estimate

displacements

load level

energy

a way

d!l:‘_,

We

require

,

kV!,\‘_ ,

E
CJ=E,E~+E&E-EJ

for

that the affinity

dj,‘,‘_, is approximately

the be-

set of

maintained.

we can kvrite

dj,\‘_, satisfy o = E,E for

I will violate

If we select

and an improved

d:,:‘_ , = <.,rr_ , Id::’

-E

p,,,_

identity.

E>E~

that

these

= rr,,:;,

where

Cv,,\‘_ , represents up to load

energy.

displacements

of energy

equation:

= 11/L + UT’ + u,:‘-,:,:‘_ ,.

forces

the pertinent

improved

the conservation

(4)

+ 1\ d!;‘- , ).

linear

level

the real p,,,_,

\vork

. while

and nonlinear

(.‘;I

of external

Ul,‘!_

I contain\

parts of the strain

An energy

balancing

srrategt

for solution of combined

promstricai

and materisl

nonlinearity

problem<

I49

load

Line

- search ,

Secant- EULER Iteration

,

P m+l

Structural

level

‘m

d

d’o’

If)

Fig. 7. Accelerated

m+l

d i splacemenk

m+t

increment~i/it~~t~ve procedure.

The nonlinear work of the external 3) can be expressed as

forces (Fig.

Expansion nomial

order

of eqn (5) yields a fourth

d(f)

d(l)

m +I

m

d

solution

(lo,

4.

1.2

=

-

2 ( I!$)!!:‘+,

(.3 ,

(8~)

polyThe linear and nonlinear portions of the strain energy of an element e can be expressed as tU$,!?t

with the following

I = him+ ,.sd!i?,.r

coefficients:

x IK%!!‘+ (uii:)=ty

r.Jld!?+1.e

(9a)

m- t.,&,O?T, .e x ~Ki%sE’,“‘t ,.e,djtj? t.Jld’,o’-1.r (9b)

(vi:‘:)!?+ 1st. approximation

I = i+ym +.,.,dill:,., x DG?‘f;(J%- ,.,,d’,O’: ,.,)I&%. ,.p

(9~)

The multipliers 116 and l/l? in the strain energy expressions follow from the fact that the nonlinear portions of stiffness matrices are linear and quadratic functions of the displacements[4]. The weighing factor

Fig. 3. Work of the external

forces.

expresses in a very simple but effective way the difference between the energies of linear elastic and hyperelastic structural materials. respectively (Fig. 4). To obtain the current axial strain for each pinjointed bar element, we utilize its nonlinear form.

load

incrementarion

i5 possible

indrv

c’,,, - ,. experience

actiLr

computer

especially’

program

true

in cxvz

his. to be discussed Ia\\

(Fig.

deviate Element

ii more

at’tinitb an inter-

efficient.

of nonlinear

This

stabilit),

representation

I) the second

radically

the that

from

El.

of the constitu-

E modul.

EI. does not

the use of

“e” min c,,, _ , = 0.90 +

0.95

and

(13)

, = I .05 ---t I IO

max c,,, _ can be tentatively crement since .-\ssuming

the

validity

of Green‘s

strain.

is

anally-

later.

If in the bilinear tive

via

has shown

we can

Ivrite

Ip,

recommended.

makes

the assumed

placements

affinity

d’,“’ and

d’,” is usually

The first

an exemption

only

between

their by

load in-

in this

the linear

nonlinear

relatively

respect. dis-

counterparts small

load

step

maintained.

E:,:‘_,.<,= I/l., t : lI~.i,!,:‘,.<.

3. I .4 Eyr~i/ihri/,rt~

Cj = 0. 1.2..

.).

described

(II)

true It should .s~tr~l

be noted

modulus

that

in eqns

of elasticity.

(IO)

the

,t’,::‘+ ,,,,. pertinent

(Y) and

to

the element displacements d)!:‘_ ,,(.. hence ,,,I c,,,, ,.<.. are used. Furthermore. the proposed balancing

technique

equilibrium load

as well

path as the stress-strain

level.

Consequently.

equilibrium creases.

integrates

path The

deviation

decreases

scalar

as the

affinity

found

by solving

u’e select

energy

viation

from

librium

iteration

solution

identity unity.

The affinity

index

That

the

load

since they

such

as the

depend

stage

of

analysis-to

ones.

If large

structures

that suitable

determined

deequi-

to guide

functions

the multipliers

of the’di;placekents

secant-Newton

the

nodal

forces

type.

trhich

require

the equilibrium

iteration

a complete

vari-

stiff-

and quadratic

the

the

unbalanced

from

I -

pm - I .

which

yields

additional

Thus.

the final

. . II)

I. 2. 3.

displacement

increments

[ ’ = ( ,K),:‘_, ) - ’ r!/‘c ,

displacement

vector

(15)

(16)

is

combination d:,‘,‘_, = dj,)‘_, +

nonlinearities. most

index

are

test problem

ot

load increments c),cles

to satisfy

control

Ad:,:::‘.

(.j = I. 1. 3. .

it is recom-

can be determined.

2

, =,

important

for the affinity

automatic

113 in this

d[4]. Employing scheme.

are obtained

1 dj,: I

limits

a few iterative

conditions.

(I?)

and lower

In this ~vay. optimum only

.)

is acceptable:

are involved.

limits

(14)

between

on numerous

the

l/2 and

K.l’ls ,md K,vL are linear

ness matrices

2).

how-

on the basis of a small

similar

even

name

. II)

are due to the fact that the nonlinear

In general.

and the degree ofgeometricimaterial

though

cj = I. 2. 3, . . Again,

falls

structure/loading

,I.

JL’,,:‘e,, d:,:‘:,

c.j =

ever. no hard rules for these upper

mended

+ I[Ki“c

, ,]

dj,:‘,

is used.

reducedtl,

size

matrix

+ &[KFLC $,,‘,‘-I.

1,Ij I = ,S!!‘- Id!:‘r,,,

increment

secant stiffness

, = [KL( ,Ey_ , ,I

is

limits

p,,,-, should be reduced.

can be given. ables:

\K!l’&

> c,,,_ , c miix c,,,_ , .

selected

loads.

equation

(/,I = I. 7. 3.

otherwise

to

external

real roots.

the following

is. if c,,,_,

the

the funda-

the smallest

can also be utilized

prescribed

min c,,, _ ,

(5). displays

the updated

total

well

required

2):

in-

, which

it is further

against

the above

quite

in eqn (4) is

by satisfying

In this way.

level

c,,,_,

can be substantially

algorithm.

certain

at each the true

eqn (7). Of all positive

the one which.

mental

curve

forces

For this purpose (Fig.

Although

approximates

position,

internal

energy

from

vector

check

ifcr.trtio~r.

search

equilibrium

as the forcr/

load

index

used to scale the displacement

strains

line

AIof the

The

iteration

is terminated

when

the

. . II).

unbalanced

nodal forces reach negligible magnitude. Each iteration in the above described New.ton

method

sion of the secant

requires stiffness

the assembly matrix

(17)

,K!,j’_

secant-

and inver-

1.

Iii

.\n e~~trgy balancing strategy for solution of combined geometrical and mareriai noniinsarit~ problcmj

This procedure. however. can be very expensive in the nonlinear analysis of .‘real-life” structures. Therefore. the size of the load increments should he selected so that the number of iterations can be restricted to a few. It is. of course. quite possible to alternate global energy corrections with local equilibrium iterations. That is. after the first energy balancing a one step equilibrium correction is introduced. followed by a second energy correction and a subsequent equilibrium iteration.

Basically. the nonlinear load-displacement response of beam and frame structures is obtained as described in the forgoing section. In this case. howevser, the additional effect of axial force on bending should be accounted for. since large deflections of beams invariably introduce axial forces. frequently to the extent that they predominate. We restrict our consideration to prismatic beam elements with longitudinal plane of symmetry and, assume that plane cross sections remain plane after the loads are applied. Under these circumstances. again, a uniaxial state of stress exists allowing the direct use of onedimensional stress-strain relations: (18) A simpli~cation of eqn (18) by means of bilinear constitutive relations will be discussed later. Since the cross sections of the beam remain plane, the elongation and/or contraction of any fiber (Fig. 5) can be obtained from

must be considered ture relationship

in the Bernoulli-Euler

(j

=O.I

curva-

’ s. 3 * . . ). . -.

Cl)

where I denotes the moment of inertia and E,.represents the so-called reduced modulus of elasticity. respectively. Timoshen~o has demonstrated[~j that, in case of pure bending. material nonlinearity can also be treated by eqn (21). if ive replace the linear modulus of elasticity by an equivalent substitute value. This concept is based on the prevailing analogy between stress-strain curve and bending stresses along the depth of the beam It. if it is replaced by /IK in the U---Ecurvef5]. Figure 6 illustrates an extension of this simple, but very effective idea to the case of bending with axial force by introducing a coordinate transformation. That is. if we enter into the U--E diagram using elongation E,, caused by the axial force f I I). then we obtain the origin of the new n’--f’ coordinate system. It can easily be proven that the shaded area in Fig. 6 represents the distribution of bending stresses in the cross section. Denoting by 11, and !I: distances to the upper and lower fibers from the center axis. respectively, elongations of the extreme fibers for each cross section can be computed from eqn ( 19). Thus, in case of bending with axial force, the reduced modulus of elasticity can be defined by

(E,.2f!r:: I.‘. = (E,, 2 %)I,:‘+ ,.P = (l, *.r + .!,(, :)I111 i’c I ’ ..a

(j = 0. I. 2, 3, I

2 (7K)ljl ., 1>! * I .2. cj = 0. I. 2.3..

. .).

(19)

where (c represents the curvature of the central axis of the element, and subscripts (I and h indicate axial and bending strains, respectively. In the elementary beam theory the curvature is linearized by assuming that the square of the slope is small. thus negligible. This assumption. however. no longer holds when the beam is bent with large displacements. Therefore. exact expression for the curvature

(j = 0, I. 7. 3..

Fig. 5. Stress and

. .) (20)

.I.

(x!)

where C = IZ for rectangular cross section. If. instead of a rectangle we have any other symmetrical shape of cross section. values for the constant C can be readily obtained from the pertinent literature[5]. Next, aiming for computational efficiency. we introduce a simple section-type (vs layer-type)

+..-_--

-...-_--_

Fig. 6. Integration of c’-e’

curve.

model for treating combined geometrical and material nonlinearities (Fig. 7). By considering several cross sections. the effect of material nonlinearity is distributed over the length of the beam element. Further simplification is obtained by the use of bilinear constitutive relations. Thus. bve replace the actual stress-strain curbe by two straight-line segments, as shown in Fig. I. This stress-strain idealization can be expressed by

+5/ 1 I z

u = E,c

for

E > E,

(23)

for

(23)

6.

/

.”

Fig. 7. Conventional

beam element

1 with

five sections.

and cedure: u = EIE~ + EJ(E - E,)

E > E,.,

where u, and Edmark the intersection point of these t\vo segments of slopes El and E2, respectively. Such a simple bilinear stress-strain approximation leads not only to very accurate detlections[3]. but can further simplify nonlinear analysis. That is. if both E’/.‘~and #j. for a beam element are within the region of proportionality [Fig. 8(a)], and the constitutive law is identical for tension and compression. the same algorithm is applicable as in the case of linear elastic material properties[Z]. Supposing that ~1.~and ~2.‘. fall within the second straight line of slope E2. representing gradual yielding, then the reduced modulus of elasticity E, can be replaced by 15: [Fig. 8(c)]. Only in the general case [Fig. 8(b)] is it required to determine E,.,,. and the pertinent tlexural stiffnesses. as described above. This simplification can easily be coded. For the simple six degrees-of-freedom beam element (Fig. 7). the commonly used linear and cubic displacement functions

( f$)h.<,=

I

E-l

~_,, K,(E,_I)K; d< h’:‘(E,I,.) IV;‘. dc,

=Z

(i. j = 2. 3. 5. 6).

(27)

where

(O = location of beam section) (see Fig. 7) (28)

(25)

with

bl

Genercl

case

(26)

have been utilized to develop linear and nonlinear element stiffness matrices[2]. In generating stiffness coefficients pertinent to the linear bending stiffness matrix I&.. we apply the virtual work pro-

Fig.

8. Reduced

mduli of elasticity in case of bilinear strebs-strain rrlation.

.-\n energ> balancing

strateg)

for solution

of combIn&

gsomctric;ll

UICI material

nonlinearit!

modulus of elasticity ~.I?!),!,,_,should be utilized by their evaluation. Similar steps should be taken by the equilibrium iterations. Thus. in eqn Ill). the geometrically linear part of the stiffness matrix should be replaced

represents the reduced moment of inertia. In eqn the symbol ( 1” denotes differentiation with respect to 5. Table I contains explicit expressions for duYdE and d’\t,idc’ to facilitate the evaluation of reduced bending stiffnesses E,- I,.. considering five evenly distributed sections along the element axis .V (Fig. 7). Consequently. any standard subroutine for determination of element stiffness matrices for bending Kit,.. which treats variable flexural rigidity can be utilized. To compute the affinity index c,,,_, , we separate the membrane and flexural parts in the geometrically linear element stiffness matrix: 117)

by CKL)!!‘+ I = cK;f + K:):!‘_ ,. cj = I. 2. 3. .

d,l::!,,-,

+ 4 W%%.,

I,)

dd!!,-, .

(30)

where (d,T?j, - I = Id, &%‘!!n- I

(31)

represents only those element displacement components which are connected with the external axial forces, while

denotes their flexural counterparts (Fig. 7). Equations (9b) and (9~) are still valid, provided the nonlinear element stiffness matrices for beams are used. Again, explicit expressions for K?:(d) and Kzf;.td’) are given in Ref. [?I. Since these matrices are functions of the cross section area A, rather than that of the moment of inertia I, the pertinent secant

I.

Explicit

det ,K$’

[det if

Equation

expressions

buckling.

d2w

5

3 L

1

-&d&d,+

d,-

+

d,+

Furthermore,

as well as for the derivative

I 51

+

5

d,+

-$

d,-

+

d,

f2

+xd3+161

- &

d,-

A

d 6

fed,-

+ds+

_Ld6 lfj

+

-+d,-&d,+$d,-+ f

73

+d6

d,

d,

_-

d,

+d,+ L

;

d3

d 2+sd3-

+

+

+d6

+d5+&

d,-

Ad,+ r’

(35)

the equilibrium is unstable cc).

2

-

l8

3

,

for IV anJ W”

dE

d3

-$d2

0

at the critical points tb)

(35b) is valid for bifurcation

snap-through

5

1

(34)

the equilibrium is stable (a)

TK!? = 0 1 ,

w”=

z-

+ K,$tE,. I,.,:./::,,_,

,

> 0

det TK’r’< 111

!+

0

similar

3.3 iVorf/i/retrr sl~thilir_YtrlIcc!\sis The previously described straightforward incremental iterative solution technique is applied to trace the load-detlection path in the prebuckling region which is often referred to as primary curve. In the “predictor” phase the current tangent stiffness matix of the system ( I) is computed. The behavior of the determinant of this matrix can give valuable insight into the state of equilibrium pertinent to the load level already reached by incrementation. The following classification is well known: if

Table

(33)

which also contains both types of nonlinearities.

and revise eqn (9a) accordingly. Thus, the pertinent part of the element-strain energy can be written as = i Ym-,.e[d:K;(,E)

.I.

where, on the element level. a relationship to that given by eqn (29) exists: (Ii,4 )!!I, I = [Kf;(,E)

L’f’“’

153

problem>

d6

$

d,

I’4

K. S/Ii.

of the determinant bvith respect to the load factor h can further help to distinguish between snapthrough and bifurcation buckling. That is. if

we are dealing with snap-through other hand. if

;

buckling.

det TKj;,’ = finite negative

On the

value.

(37)

the buckling is of bifurcation t)pe. Since eqn (2) is normally solved by some decomposition method. the determinant. dct ,K!‘,‘. is obtained as a by-product of this procedure. Using the Gaussian elimination. for instance. product of all diagonal elements in the triangularited matrix gives the determinant of the current tangent stiffness matrix. Consequently. by monitoring the determinant (including the signs of the diagonal terms) we have a very good measure for detecting critical points. That is. observed low determinant values indicate an approaching singularity. In the near vicinity of the critical point. however, the solution should be brought to a temporary halt to consider a possible change in the subsequent solution strategy. to be discussed later. Another by-product of the energy balancing is a stability indicator. described in Ref. [Z]. This scalar quantity can also be used for advanced detection of limit points, in the case of snap-through buckling. In spite of its pronounced similarity with the “current-stiffness parameter.+’ introduced by Bergan[6J, they are not identical. Although both are very useful in prediction of limit points, their application is restricted to snap-through buckling. 3.3. I Bifilrc~rriorr hrwklirq. Experience and experiments show that pure-bifurcation bucklings of real-life structures are relatively rare phenomena[7]. In spite of this fact. we briefly discuss here how the present method is applied in the case such an stability problem occurs.

\HI)

The smooth primary load-detlscrion curve is traced using the present incremental-iterative procedure, as described above. Since the determinant of the tangent stiffness matrix is positive for the stable part of the equilibrium path and zero at the bifurcation point (351. continuous observation of the determinants and that of other stability indicators is considered to be mandatory. This should be coupled with the simultaneous control of the affinity index. In this way. the analyst is warned in advance of the approaching bifurcation point. A simple extrapolation of ths determinant curve is all that is required to establish the critical load (Fig. 9). Determination of the secondary unstable branch, however. requires an indirect technique. whereby the original structure is slightly modified by introducing initial imperfections. This well known approach transforms the bifurcation problem into an “equivalent” snap-through buckling. as shown in Fig. 9. Consequently. we compute a smooth load-deflection curve of the slightly imperfect structure. The postbuckljng region of this deformation path can be used for approximation of the two intersecting branches (Fig. 9). Of course, this approach is only applicable for simple (vs multiple) bifurcation points. In any case. the equivalent limit point usually gives a safe estimate of the critical load. The technique to be used fo determine the limit point and the postbuckling path is discussed next. 3.3.2 Stop-rhro/r~h hrrc,klin~. Again. the nonlinear prebuckling path is determined as described above. Since at the limit points the tangent stiffness matrix becomes singular the “predictor”-eqn (2) cannot be applied to obtain the first estimate of the incremental displacements. The energy balancing strategy can, however, handle this. otherwise difficult. situation with relative ease and permits the tracing of the postbuckling path without changing the algorithm. As the determinants and the other stability indicators signal the approach of the limit point. the load increments are subsequently reduced. In the

bifurction point

displacements ai Load- displacement

curve

det bi values of determinants

Fig. 9. Bifurcation buckling.

vicimty brought already estimate

of the limit point. the computation is to a temporary halt to slightly reduce the achieved load level p,,
deflection curve in the postbuckling region. ibe use a linear combination of these tuo displacement curves (10)

d:;:‘_, = d,, + @.

by d!:;‘_, = kd:,‘;’ .

p. = (I.ol-l.usI.

(38)

Equation (38) now replaces eqn (3) in the algorithm. The line-search phase. followed by the equilibrium iteration. remains unchanged. In tracing further portions of the postbuckling path, the solution strategy described in Sets. !.I and 3.2 is used until the next limit point is encountered. Here. again. we apply the above given strategy. The method of solution for combined nonlinear instability probIems, discussed above. can be applied successfully as long as the buckled mode is similar to the one obtained prior to reaching the limit point. In some cases. however. after passing this point. the deflected configuration of the structure suddenly changes. This phenomenon is usually signaled by one or more llnproportionally large displacement components. Such a situation often occurs when in the prebuckling range the geometric nonlinearities dominate but at the limit point a “plastic” buckling takes place: since. at this stage, the unstable state of equilibrium is mainly due to their geometric counterparts. In this case, the continuation process beyond the limit point must be temporarly changed. This. however. may require some insight into the failure mode of the structure. For this purpose, we compute the lowest eigenvalue w, and the pertinent eigenvector 6r, using the tangent stiffness matrix rKcrat the critical point. Again. simple extrapolation of the determinants or one of the stability indicators helps to locate d,,. Since the inverse vector interation simultaneously yields wI and 8,. this simple technique should be applied. The new buckling mode is

5=

cy =r:II dcr - &-I

CA,.

/l~l II

To establish the first approximate

Lmt

.

(39)

point of the load-

point

Fig. IO. Passing the limit

tl

point.

where p is a ueiphing factor indicating the degree of participation of the new buckling mode in the superimposed deflection vector. The weighing factor p is computed from the incremental energy loss .I Wh I?,- I = Au* 01- I which expanded

gives the following

b,, + b,P

+ 6:-S= f

h$”

(41) polynomial: = 0

(41)

with the coefficients

(43)

The lowest real root of the above polynomial is the weighing factor required in eqn (40). Using this first estimate of displacements, the line search and the equilibrium iteration parts of the present method are carried out as previously described. Needless to say that this special procedure to overcome the limit point is more expensive that its regular counterpart. In addition, this special phase of computation requires a high degree of man/cornputer interaction. Fortunately. the need for such a procedure is an exception rather than the rule. Furthermore, in the engineering practice the stability analysis is usually terminated after the critical load has been determined: i.e. tracing the postbuckling path is only in isolated cases required.

.8. COSSIDERATIOY

OF L.ARCE ROT,ATlOSS

The total Lagrangian formulation coupled with the use of Green’s strain, is by far the most widely used approach in application of finite element methods of nonlinear structural problems, because of its simplicity. However, this combination also entails some disadvantages. In Green’s simplification of the strain-displacement relation ( I I ). for instance, all higher order terms are arbitrarily neglected except for the one quadratic rotation term. As the structure undergoes excessive deformations. these higher order terms gain in importance: thus. they cannot be further neglected. In addition. in case of large rotations. relation of external loads to the deformed structure might undergo considerable changes making the structure nonconservative. Hence, application of the present method to such structural problems requires the following exten-

1%

K. S/IL.

sions: ii) Additional con5idrration of burns of the neglected cmali terms in the strain-displacement expression. Consequently,. the stift’ne3, matrix of the system must be augmented by K,)~Id’) and Kkctd’ 1. tii)

Since the load vector is no longer constant but depends on the gradients of the deformations, it should be continuously modified during the computation. (iii) In the case of displacement dependent loads. introduction of an additional “load correction matrix ” is also requiredf8-1 I]. Work is currently being made to include all these extensions needed for more exact consideration of large rotations. It is hoped that the results of this effort can be published soon.

A FOKTRAS program covering the present method has been first implemented on a PRIME 650 minicomputer. and subsequently on a VICTOR 9000 microcomputetj I?]. The fact that such reasonably priced upper-class micros are nowadays widely available with virtually no runnjngcosts. makes this alternative very attractive for smaft engineering consulting firms. Furthermore. considering the actual needs of practicing engineers, the program has been so developed that it can also treat. in addition to the here discussed combined nonfinearities. separately occctring geometrical or material nonlinearities, using the very same approach. The basic limitation in using microcomputers for nonlinear structural analysis, is the small rapid access memory (RAIM) available. This compels the programmer to the extensive use of slow access memory in the form of floppy disks. with the consequent penalty in length of computational time. Although, the VICTOR microcomputer can be expanded up to 896 KB core memory. this relatively large R4ivl is, however, divided into I4 addressable segments of 64 KB each. Consequently, an extensive moduiarization of even a moderately sized finite element program is mandatory. Experience shows that the most critical and time consuming part of nonlinear structural analyses using small computers. is the solution of a large system of simultaneous equations. The computational time percentage devoted to this task lies between 70% and 80% of the total time. The algorithm developed for solution of a large systems of linear equations is a hypermatrix variety of Cholesky’s method. The reason for this selection is manifold. First of all. the common characteristic of the stiffness matrices used in the present method is that all are symmetric, positive definite. and banded. Thus. making use of these characteristics. the number of matrix operations could be substantially reduced. Secondly. Cholesky’s method yields the determi-

\KI>

nant of the coefficient matrix as a by-product ofthe triangularization process. Since determinant of the tangent stiffness matrix ptays (as mentioned above) a dominant role in the nonlinear stability analysis. this item alone can be a decisive factor in the case of very large systems. Last but not least. the method is transparent and compares favorably, with other similar type of solution procedures. such ai Iron’s frontal method, for instance. All stiffness matrices are divided into equal blocks. TO perform the factorisation. however. two of these blocks and the pertinent program must occupy one of the 64 KB segments at the same time. This determines the maximum size of a hypermatris. By large structural systems some of the hypermatrices are null. These need not be stored OI operated on. Even inside of nonnnull submatrices. the operations with zeros are avoided. Using such segmentation and overlay techniques, nonlinear analyses of large structttr~~i systems can be carried out with microcomputers. This appears to be one of the most appealing features ofthe present method.

I‘o test the accuracy and cornptlt~Iti~)n~~i efficiency of the proposed method. cable. truss. beam, cable-beam. and arch structures have been analysed and the obtained results compared with their analytical or alternative solutions. Due to the extreme mathematical complexity in achieving quantitative results for structures exhibiting combined geometrical and material nonlinearities. few exact solutions to such problems appear in the literature. Those solutions which are available deal with relatively simple structural contigurations, loads, and constitutive laws. respectively. The test problems discussed here are of similar nature. Erc,mplc I. A simple truss problem, shown in Fig. I l(a) has been chosen to test the effect of the proposed energy based line-search technique in connection with the accuracy of weighing factors given in eqn t 10). in this case. the effect of material nonlinea~ty predominates. For purposes of comparison. an analytic solution of this relatively simple problem has been developed. Complete agreement of the computed results vvith those of the analytical solution has been achieved. The obtained load-deflection curve for nodal point 0 is illustrated in Fig. I l(b). One interesting point deserves special attention, namely in the linear elastic region (up to load level 140 kN), for all practical purposes no equilibrium iterations are required. Ev,en up to load level 160 kN. where the truss already exhibits pronounced inelastic behavior limited cycles of secant-Newton iteration suffice to restore compieie nodal equilibrium. Since in the region of gradual yielding relatively large load increments have been utilized. in this region the number of rcquirrtl it-

tp, al

Apphed load

ido

Structure

[kNI

deflection

b) Load-deflection Fig.

I

I.

Large

de&ctions

erative cycles became large. By proper reduction of the load-increment sizes. such computationally undesirable effects can be eliminated or. at least. mitigated. Esamplr 2. The suspension bridge, shown in Fig. I?, exhibits both geometrical and material nonlinearities. An alternative solution of this cablebeam problem is given in Ref. 1131. The original Ramsberg-Osgood type constitutive law has been replaced here by an equivalent bilinear stressstrain approximation. A comparison of the exact expression for strain energy[l3] with eqn (9) has shown practically no deviation. Similarly. excellent agreement has been obtained for all cable forces (Table 2). The apparent small errors in the girder bending moments are due to the fact that the

300cm

I 300cm

I 300cm Fig.

diagram of il truss.

present method considers large deflections in girders, while Ref. [13] does not. E.rcu~rplc 3. Next, inelastic deflections of a simple supported beam, shown in Fig. 13. have been investigated. In this case only material nonlinearity has been considered. The purpose of this investigation was to test the applicability and accuracy of the proposed method and, hence. the validity of the pertinent computer code developed for the general case. Again. the hyperelastic material behavior has been approximated by a bilinear stress-strain relationship. Comparison of the center deflections with the exact load-displacement curve[3] indicates that the proposed method not only leads to accurate deflections, but to considerably reduced numerical work. That is. convergence to the exact solution

I 300cm

12. Suspension

I

300 cm 1 300 cm bridge.

I’S

CABLE

t

Ilemen

FORCES

Present

I

method

I

216

q

I

L86

iL

, i

02

Fig.

anal

00

q

216

216 71

001

L86

BENDING

OL

Error %

0

26 632

67

0

II.

036

39

298

35

0

1L 036

27

0

26.632

68

-Ll

13. Simple

--__ 08

f

text

E. :

106

10’k,ps/,n2

ET=

2 53

lO”k~ps/m~

oy :

52 k@,n’

b I

21n

L =

101n

h:

51n

I31

pose of this study was twofold. First. the general validity of the computer code has been tested by assuming a linear elastic stress-strain relationship. In this case. even with relatively coarse discretization and large load steps excellent agreement with the exact solution could be obtained. It is significant that without any equilibrium iterations the error in the maximum deflection was merely - 0.5%. Using only one cycle in the iterative steps after energy balancing, errors related to the internal stress resultants became also insignificant. The second application of the computer code involves determination of the additional nonlinear effects caused by an assumed bilinear stress-strain relationship. The results for the center deflection indicate a considerable softening effect due to the gradual yielding of the beam. To check the validity of these results. an approximate analytical solution

D : 3LL 25 N/cm 15

10 :

method

5:

06

[ kNm 1

MOMENTS

method

solutfon

present

00

Present

.

0

!

Pofnts

I

-

197L2 197 LS

I

/

00

Nodal

5 *

2 .

LB

197.39L

could be achieved by using merely IO simple beam elements of variable flexural rigidity, as described in Sec. 3.2. It is of interest to note that 4 times as many beam elements are required to obtain the same results if the moment-curvature relationship, given in eqn (2 I ). is determined only for the middle of elements. Eirrmpk 4. A uniformly loaded beam of rectangular cross section with freely rotating but not sliding ends, as shown in Fig. 14. has been considered for the next investigation. An exact solution of this geometrically nonlinear problem assuming the validity of Hooke’s law is given in Ref. [ 141. The pur-

3

%

a

+ see

‘

1 Error

1

I

0

I

13

216

197 39L

GIRDER

p, 1 kips I m2

[ kN 1 Reference

10

>upportcd beam Hirh nonlinrarit~.

E. E-: E, b L h

:

136MNlcm’ 02E. = L 10.’ : 2Ocm ___ 1 BOOcm : 30cm

/

,,si’

_,;/ / _Y

--iE.E,l’ -0.

appr

method Onal

linear

WC Im I material Fi_e. I-!. Effects

of combined

nonlinsarities

Solutlan

PlkNI

Exact iv.

I

L 100 E,= 206%

I 10gkNlm2

Ez: OOlE* Es” 1%. b -1000m

Fig,

Ifi.

Deflections

of a beam with fixed ends.

has been developed using current secant moduli for stretching and bending. Considering the approximate nature of the analytical solution. the deviations are small and in the right direction. To restore nodal equilibrium after line search, at most, six iterative cycles are required. E.~rnpie 5. The next problem involves a ctamped beam under a concentrated load at midspan. Dimensions of the beam are given in Fig. IS. The nonlinear material properties are, again, represented by a bilinear stress-strain relationship. Because of symmetry, one-half of the beam was

subdivided into 10 finite beam elements. Again. in the first part of the analysis only geometrical nonlinearities have been considered. As the computed load-deflection curve indicates an excellent agreement with the analytical solution of this problemjIll could be achieved. in spite of the relatively large load increments. At all load levels. the proposed line search has approximated the true equilibrium condition quite closely. Consequently. one step equilibrium correction was sufficient. The combined effects of geometrical and material nonlinearities are similar to those discussed above. with the exception that gradual yielding of the beam takes place at a somewhat higher load level. In this case, the geometrical nonlinearity predominates. in spite of the fact that the selected bilinear stressstrain relationship with E2 = 0.01.15, quite closely approximates the elastic-plastic material behavior. In the region of gradual yielding, the load increments are reduced to obtain the correct ioad-deflection curve without excessive equilibrium iterations. Escmple 6. Figure I6 shows a very shallow clamped arch, loaded with a concentrated load at the apex. In modelling the structure. 20 straight beam elements have been used for half the span length. Considering only geometrical nonlinearities, the equilibrium path is always stable for this structural configuration. Consequently, there are no limit points on the load-deflection curve. Comparing the obtained results with the exact solution

?

F

E,= IO7 F : 297OL I_ L: L838 R : 100 h I 2.0 brl0

PS’ I” in I” I” m

10000

20 Fig.

16. Very

30 shallow

LO arch.

50

w,linl

Fig.

17. Snap through

of the problem[lj], the agreement appears to be very good. Next, we consider the effects of combined geometric and material nonlinearities on the equilibrium path. As can be seen from Fig. 16. the response of the structure is now completel) changed, since the pertinent load-displacement curves have limit points. Furthermore. decreased E2 reduces the critical load. Extunplr 7. We next consider a shallow arch whose nonlinear response is characterized by limit points, as shown in Fig. 17. In the present analysis. the arch was subdivided into 20 straight beam elements over half the span. The obtained results for the central deflection involving only geometrical nonlinearities are in excellent agreement with the exact solution given by Timoshenko[l6]. Introduction of material nonlinearities lowers the critical load, as shown in Fig. 17. In tracing the complete F= 109in Ls 31 I" R = 133111 I"

50

15

E,Z

:

-

2%.

E = E,= constI171

.,

)

LO 1 35: 30 i

10.'

Fig.

02

OL

06

08

18. Load-deflection

10

12

curve

1L

16

16

of a shallo\v

20

w,ltnl

arch

of ;I shallow

arch.

postbuckling paths no numerical difficulties have been encountered. Er~/t~p/e 8. Two shallow circular arches subjected to concentrated loads at their apex have been selected as test problems for elastic and inelastic stability analyses. The parameters for both arches are taken from the paper of Belytschko and Glaum[ 171. First. the arch shown in Fig. I8 has been investigated involving only geometric nonlinearities. This structure was modelled by 40 straight elements. The prebuckling part of the obtained loaddeflection curve (including the limit point) shows good agreement with results published in Ref. [ 17). The agreement in the postbuckling region can be considered as reasonable. The discrepancy accounts for the different formulation and elements. That is, in Ref. 1171 higher order curved beam elements have been employed. whereas the present method applies simple straight beams. In addition. the deformed configuration in the postbuckling range involves significant rotations which can be better accounted for by a “higher order corrotational formulation”[ 171 than by the total Lagrangian description coupled with Green’s strain. Since in this case the geometric nonlinearities govern. additional consideration of material nonlinearity with E, = 0,04.E, in the bilinear stress-strain relation does. not significantly change the full load-deflection curve (Fig. IS). By increasing the depth of the arch from tl = 0.1875 in. to tl = 0.59 in., the nonlinear response of the structure was completely changed as well with elastic as bvith elastic-plastic materials. as shown in Fig. 19. With constant modulus of elasticity the equilibrium path is stable in its full length. The pertinent deflection curve compares reasonably well with that obtained by Belytschko(l71. The

balanccng >trateg) for dution

An cnere)

of combined

$+&

and m‘lIerl;il

nonlinzarir~

Ihl

problems

; i ;;Y‘;:

~

pl

grometrkll

b : 1001n E,= IO' PSI

I

Ilbl

Ez= OOL E, E, = 2%0

1600; !

-

E = E;: const[171

--A--E : E,: const.

05

discrepancy

is mostly

differences

Comparison the results Again, the

good

prebuckling

Similarly.

agreement range

including

the

method

corresponding

with

the critical

shows

results

suited

struc-

close agreement published

for

tracing

without

The various in

load.

of this

solutions.

proposed relying

in Ref.

nonlinear

method upon on

however.

been

problems

Some

a meaningful

This paper describes simple of

modelling

planar

achieve

technique

framed

with

deflection sociated

dof straight

element

beam

An energy much than

methods. cycles

in nonlinear

these tions ularly Final

for

selection

the

of load

down

load

final

level

displace-

it

balancing

provides

a rational

increments.

In spite

excessive

required.

This flow”

is that the proposed point

techglobally

The

such

should

to perform the present its re-

modelling

tech-

restrictions

of the

are:

to large deflections large

(ii)

The

loading

(iii)

No

unloading

relation

formulation with

used.

limits

in the

the method

small strains

and mod-

rotations. must

be conservative.

from

the plastic

state

can take

place. Although,

these

handicaps

for a large percentage

nate them.

with

work

IS]

such or

the present

do not present

is presently

Furthermore.

to combine

techniques. proach[

limitations

as

some energy

under

it appears

other

arc-length

extrapolation balancing

way

used

to elimi-

to be quite

convergence

the

serious

of structures

lu-

accelerating control technique[

ap191

strategy.

of

itera-

is particis exited.

method

when

understand

the employed law.

imple-

difficulties,

using

must clearly

Lagrangian

erately

solutions.

is. in order

strain-displacement

total

crative

is significantly methods

approach

of the

been

analysis

analysis

and constitutive

in practice,

of the iterative

energy

large “plastic

at the limit

six

incremental-iterative

improvements.

the case when observation

to the

Newton-type

may be occasionally

not break

at fixed

the number

Furthermore.

significant

section-type

local convergence

the

The as-

That

nonlinear

and justify

(i) The

dof bar.

estimate

In addition.

geometr-

of a

herein

to large scale projects

the analyst

present

for large

involving

line search

way,

in

nique

problems:

nonlinearities.

traditional

In this

makes

convergent. way

the

to achieve

reduced. nique

closer

is to

and economy

is a simple

or a four

based

aim

nonlinear

analysis

and

analyses

main

is applicable

and combined

finite

ments

of static

introduced

and stability

ical. material,

Its

simplicity

all types

i.e. the algorithm

strategy

for nonlinear

structures.

computational

dealing

gives

a new solution

known

a microcomputer.

not be underestimated. method.

with

has even

is applied

strictions

capabilities

demonstrated

code

F0KTK.G

involved

a computer

7. CONCLUSlONS

path

the basic computational analysis

have

test

The pertinent mented

[l71.

the postbuckling

changing

strategy.

interest.

be achieved

analysis

well

structure

solution

could

the postbuckling

ture by the present

discussed

[ 171 is of particular

in Ref.

w,linl

of elastic and elastic-plastic

of the problem.

of the elastic-plastic

given

very

with

due to the above

in formulation

20

1.5

1.0

Fig. 19. Comparison

does

and is. in general.

Ac~no~l~/ecl,~~tlc~,lf.s-The author gratefully acknowledges support for this research from the Deutsche Forschungsgemeinschaft (German Research Foundation) and wishes 10 thank Mr R. Jahn and hlr D. Ziesiq for their work on the numerical examples.

REFEREWES

I. R. Szilard. olaccment

.An energ! bal;mc~ng method for large diyanal~ii~ of \truiture\. Coltrp//t. .\/~~Iif.

d.‘Szilard. Critlcul load and postbuckling analLA b) FE)1 using energy balancing technique. Ccj/rrprtr. Srr~cr. 20. 277-236 ( 1985). A. Cajes. Inelastic detlections of beams. J. Srr~cf. L)il.. .-\SCE 91, li4Y-1565 (1968). R. Mallet and P. V. Marcal. Finite element analysis of nonlinear structures. J. .Sfrccci. D/t,. A.‘YCE 9-l. ZOHI-2105 (1968). S. Timoshenko. Srwrr~r/~ l!f’ .\/trr‘,/.ici/.\--~tr~r II. 3rd Edn. Van Nostrand. Princeton I lYi6l. P. G. Bergan et (I/.. Solution techniques for nonlinear finite element problems. //cr. J. ,A’lf/,rc,r. .C/e!lr. En,v~g I?, 1677-1696 ! 1978,. D. Bushnell. Bucklinnofshells-Pittfall fordesieners. AM.4 J. 19. 1183-1~~6 (IYYII. J. T. Oden and J. E. Kev. Numerical analysis for finite axisymmetric deformat;ons for incompressible elastic solids of revolution. 1111.J. So/it/s Srrlrc.!. 5, 407-j IY (1970). H. D. Hlbbit, P. V. Marcal and J. R. Rice. A finite element formulation for problems of large strains and large displacments. Irrr. J. So/it/.\ 5trclc.l. 6. 106% 1086 t 1970). M. S. Gadala and A. E. Oravas. Numerical solutions

I I. K. II.

Schueitzerhofand E. Ramm. Dizplacmcnt Jepsndent prc$surs loads in nonlinear finite element anal~jis. Co,rrprtr. S~UCC(. 18. IOYY-I I I-l I IYX-I). R. Srilard (‘I II/.. Energy Balancing II: ;\nalys~s Geometrically Physically Nonlinear Beamr. Frame>. Cables and Arches Iin German). Report No. Sr .35;1. Institute of Structural Mechanics. Cniversit) of Dortmund. Dortmund. Federal Republic of German!

of

i I%GI. Static analysis ofgeometrically and phys13. J. Pistrrak. ically nonlint‘ar cable-beam 5tructure5. Ucl,. Rorr/r~. Sci. Tctiw. .I/cc~/r. .App/. 18, 45-53 ( lYX3). II. K. J. Roark. ~~wr~~r/tr.\ 1;~ .S/rv.vs tljr(/ .Srrcii//. -4th Edn. RlcGra\\-Hill. New York 11965). and E. F. hlasur. Buckling shallow IS. H. L. Schreyer arches. J. .Mec~lr. Di\,. ,-tSCE 92, I - 19 I IY661. and J. XI. Gel-e, Ilrc,or~ r!f’t‘/o.v/ic~ Srcl16. S. Timoshenko hi/it\. 2nd Edn. 1lcGraw-Hill. tiew York (1961). and L. Gaum. Application of higher 17. T. Belytschko order corrotational stretch theories to nonlinear finite element analysis. C~~,rrpc/l. Strrrc,l. IO, 175-182 ( IY7Y). method including line 18. hl. A. CrisPield. An arc-length seaces and accelerations. I/r!. J. :\‘rorrc,r. ,I/crt/r. Eyctrc 19. I’-W-l%9 11983,. wlution to elastica I’). T. Y. Yang. IL1atri.x displacement problems of beams and frames. Irr!. J. So/it/s .Swwf. 9, 829-841 (1973).

of