004s7949/8753.

00
+ 0.00
Q 1987RrgamonJournalsLtd.

IMPROVED ARC LENGTH

ORTHOGONALITY METHODS FOR
NONLINEAR FINITE ELEMENT ANALYSIS
BRUCE W. R. FORDE and SIEGFRIED F. STIEMER
Department of Civil Engineering, University of British Columbia, 2324 Main Mail, Vancouver, B.C.,
Canada V6T 1WS

(Received 20 January 1987)

Ahstraet-The path following technique known as The Arc Length Method’ has evolved over the past
decade into a commonly used tool in nonlinear finite element analysis. Variations of this method have
orginated from a number of other workers. This paper presents a new general formulation for all arc
length procedures. A derivation is given based on orthogonality principles which provides a new
perspective illustrating the relationship between the existing and proposed path following techniques. A
simplified procedure, obtained directly from the new general formulation, provides the same results as
Crisgeld’s explicit iteration procedure (CompuJ. Strucr. 13, 5562 (1981)) with a reduction in com-
putational effort. Practical appii~tion of this theory is demonstrated using a numerical example of a
reticulated shell structure.

1. INTRODUCTION displacement update for a conventional ‘load con-

trolled’ Newton procedure due to the unbalanced
The general goal of all arc length procedures is the loads G(j).
control of iteration in the numerical solution of
complex nonlinear problems. Application in finite Au = AR’61 + Au” (6)
element analysis involves the tracing of a complex
path in the load/displacement space into the post- Au'=K-'P (7)
critical range. Procedures of this nature require a
A,,"
= _ K-I@'). (8)
reformulation of the traditional ~themati~l prob-
lem into a different form.
The standard equilibrium equation for propor- This procedure can be applied with either a modified
tional loading is Newton, quasi Newton, or pure Newton iteration
scheme. The relationship between these components
1P-F(x)=O. (1) in the load/displacement space is displayed in Fig. 1.

This is traditionally solved using a Newton in- 2. ARC LENGTH ORDINALS

METHODS
cremental iterative solution.
A general arc length procedure can be derived from
K’“Au = 1P - F(x”‘) (2) orthogonality principles. An arbitrary update direc-

x(1+ 1)= x”’ + Au. (3)

These equations must be modified for use with arc

length procedures. Writing the proportional loading
factor I in incremental form yields a general formu-
lation:

K”‘Au - PAl = (“A + i”‘)P - P(x”)) = - Gti) (4)

x(’+ 1)= mX+ u(‘)+ Au. (5)

The incremental displacement can be written in two

components. The first is the displacement Au’ due to
a unit load factor multipii~ by the incremen~l
variation in the load level AR. The second is the Fig. 1. Generalized incremental formulation
625
626 BRUCEW. R. FORDEand SIU~FRIED
F. Srtausr~

tion II”) can be chosen with reference to the tangent

t”’ of the current incremental load-displacement A
configuration. The scalar product of these two vec-
tors yields a residual R “1. Normal and updated
normal plane [ 11, consistently linearized spherical h

hyperplane [2], and explicit spherical [3] iteration pro-

cedures all can be obtained from the results of this
general method.
The tangential and normal vectors consists of m
dimensions from the displacement vector and one
dimension from the load parameter. These com-
ponents are combined using a scaling factor fl (with
units of displacement) to form vectors with m + 1
dimensions as shown in Fig. 2.
Fig. 2. Orthogonality relationship.
$0 _ u’0 + fly”)
(9)

a”’ = Au + gti. (10)

A simple correction can be applied to the residual
The scalar product of these vectors yields a residual once divergence from the arc has taken place. The
R’“. This can be visualized as the projection of vector difference between the length of the current tangent
II”’ onto vector t”) (the length of vector r”’ multiplied vector and the desired length can be projected onto
by the length of vector t”)). the current tangent vector to provide the residual for
the orthogonality expression.
t@).n(‘)= 1e) 11
n(‘)1corj Q
t”‘.r”‘= _ ]t”‘]]r”‘]
= Itqlq
= - t’qt”) - s)
= RQ. (11)
= R’@_ (14)
Forming the scalar product of to and II@)using eqns
(6), (9) and (lo), a general expression can be derived
Substitution of this residual into the general expres-
for A1:
sion yields an update for linearized iteration on
t”‘. ,‘I) = @Au + /?zJ”‘Al spherical hyperplanes as obtained by Schweizorhof
and Wriggers [2].
= II”)‘{ Au’Al, + Au”)
Consistently linearized method
+ /PPM _ tW(to _ s) _ @‘&+I
A.4 = @*A”’+ &?+A,+ ’ (15)
= u’0rAu” + M(u”‘rAu’

The update direction provided by this method also

+ /Pa”‘).
tends initially away from the spherical path; however,
General expression

This expression can be simplifkd for particular cases

of o~hogonality. For the case of R’OSE0, the tangent
t(l) and update II”) are mutually normal at each
iteration step, yielding the expression for updated
normal planes [l]. This formula provides a series of
updated iteration directions which converge as shown
in Fig. 3.
I
Updated normal plane method I X
I
I +
in’
x
(13) Fig. 3. Updated normal plane method.
 Arc length method for nonlinear finite element analysis 627

The length of r(‘+ ‘) can be calculated in terms of the

old length and the incremental update vector.

t(‘+ ‘1= (t(‘)’ + L\urAu + f12A12)o.s. (18)

Substitution of this expression into the general
formula for Al leads directly to the quadratic equa-
tion obtained by Crisfield [3]. The problems encoun-
tered in the application of Crisfield’s method can be
Sm Path ’ avoided by solving for I(’+ ‘) using the updated nor-
mal plane expression for AI, and then calculating a
residual which returns the iteration direction to a
point directly on the sphere. This provides exactly the
same result as in [3], wifhour the solution of a
Fig. 4. Consistently linearized method. quadratic equation and the selection of one root.

Algorithm for explicit iteration on spheres

the application of the residual draws the path back to 1. Use the update for orthogonal iteration.
the arc once convergence is achieved.
Explicit iteration in a spherical path requires the
formulation of a residual based on the error that
would be obtained using an orthogonal iteration
path. This error can be corrected in advance to
provide the desired arc. The residual vector rc’+ I) is in 2. Calculate the associated displacement vector.
line with the tangent t(‘+‘) and acts in the opposite
direction. The deviation in iteration i + 1 is the Au = Au’M -I- Au”.
difference between the length of the tangent in this
configuration and the desired length.
3. Find the length of the tangent in the potential
configuration.

r”+” = (r(‘p + AuTAu + j2AA2)0s.

= (t(‘+ 1) - s).
(16)
The orthogonality condition of eqn (11) yields the 4. Calculate the required residual for explicit iter-
scalar residual R(j). ation on a sphere.

f(r).r(l+I)= - (t(‘)l)t(‘+‘)lcosa
R”’ = - _&(I”+‘)-s).
s
= _ s(t(i+l)-s)-
10 + I’
5. Return to the general formula for iteration path
direction.
= - & (t(‘+ 1’- s)

= R(“. (17)

6. Calculate the desired displacement vector.

Au = Au’M + Au”.

This procedure involves only a small amount of

extra work in comparison to that of the consistently
linearized procedure. Steps l-5 require the identical
amount of computation. Step 6 introduces one addi-
tional vector multiplication and addition. The larger
domain of attraction offered by this algorithm, in
comparison to the other methods based on orthogo-
nality, may be important in the vicinity of limit points
Fig. 5. Tangent vectors for orthogonal iteration. with sharp gradients.
 BRUCE W. R. FORDEand SIEGFRIEDF. STIEMER

3. ffuMEmC.4L EXAMPLE
Application of path following procedures to the
large problems found in engineering practice requires
a comprehensive comparison of the various methods.
An extensive series of examples were used to rank the
updated plane, linearized, and C&field methods in
[4], so only one example is given here to demonstrate
the use of the new procedure for explicit spherical
iteration.
A reticulated shell structure previously analysed by
Papadrakakis [5] provides the basis for comparison
of the implemented path following procedures in the
finite element analysis program NISA. The complex
behaviour of this structure is caused by the multiple
Fig. 6. Explicit iteration on spheres.
snapping models present in the analysis of reticulated
shells with many levels. The present study deals with
only the first snapping mode.
The complex displacement response curves found
in this problem provide a typical application for arc
length procedures. The results obtained by Papa-
In contrast, the approach of C&field [3] can be
drakakis agree completely with solutions obtained in
derived starting from the definition of the current and
the present study. This problem is restricted to elastic
following load/displacement configurations:
behaviour (even though member buckling may be
ff’)Z= /32A(tlr)2
+ p(‘)r”(‘f prevalent) simply to provide a means of comparison
(19) to an established solution.
A variety of analyses were performed using
I(‘+ ‘P = 82@(i) + Aa)2
different path following methods. Four arc length
procedures were used:
+ {I&‘)+ Au’ Al + Au”} r {II(‘)
1. Updated normal plane iteration [l].
+ Au’AA + Au”j. (20) 2. Consistently linearized iteration [2].
3. Explicit spherical iteration using Crisfield’s
Crisfield [3] utilizes a further simplification of /I = 0
method [3].
at this stage; however, this is included here for 4. Explicit spherical iteration using ortho-
consistency of the derivation. The two configurations
gonality [4].
t(‘) and t(‘+ ‘) are identically equal to the prescribed
value of the arc length. Combination of the two Each of these procedures were combined with three
expressions yields a quadratic equation: Newton update methods:
1. Modified Newton.
aM2+bAl+c=0, (21) 2. BFGS.
3. Pure Newton.
where
The total number of iterations required for equi-
a = @* =L A#Au’ librium convergence using two different arc lengths
and the total CPU time required for computation are
b = 2@*1”’ + u”“Au’ + Au”~Au’) listed in Tables 1 and 2.
The results of the analyses indicate that con-
c = ~u”‘~&J’* + Au”~Au”, vergence is relatively independent of the applied
arc length procedure as found in 141, and that the
The solution yields two roots for AR: new procedure for explicit iteration provides results
identical to that of Crisfield’s method with an 18%
reduction in the time spent computing ArZ.

4. SUMMARY

The selection of the appropriate root depends on Arc length procedures can be categorized using one
the current tendency of the load/displa~ment curve. general formula from orthogonafity principles. Sub-
Precautions must also be taken in the case of complex stitution of a variety of assumed or calculated scalar
roots. The solutions of these problems are explained residual quantities into this formula yields the well
by Crisfield [3]. known ‘updated plane’, ‘consistently linearized’, and
 Arc length method for not&ear finite clement analysis 629

P=l ER-IO’

Fig. 7. Geometry of a reticulated shell

Fig. 9. Vertical displacement of joint 2.

Fig. 8. Vertical displacement of joint I.

Fig. 10. Horizontal displacement of joint 2.

Table 1. Arc length = 5.0

Total number of iterations Total CPU time (see)
for seven arc length steps Vax 1l/730

Modified BFGS quasi Pure Modified BFGS quasi Pure

Method Newton Newton Newton Total Newton Newton Newton Total
Updated
planes 32 32 21 85 340 333 to31 1704
Consistently
linearized 32 30 85 341 330 1106 1777
Crisfield 32 31 22: 84 341 334 1038 1713
Explicit
spherical 32 31 21 84 339 333 1030 1702
630 BRUCEW. R. FORDEand SJEGFRIED
F. STIEMER

Table 2. Arc length = 10.0

Total number of iterations Total CPU time (set)
for five arc length steps Vax I l/730

Modified BFGS quasi Pure Modified BFGS quasi Pure

Method Newton Newton Newton Total Newton Newton Newton Total
Updated
planes 25 IS 0 262 739 I248
Consistently
linearized - 27 17 - - 254 811 -
Crisfield - 25 15 - - 256 748 -
Explicit
spherical - 25 15 - - 248 748 -

‘Crisfield’ methods or the new simplified method for finite clement analysis in structural mechanics, Ruhr-
‘explicit spherical iteration’. Universitiit Bochum, Germany (Edited by W. Wun-
derlich, E. Stein, and K. J. Bathe), pp. 63-89 (1980).
Convergence velocity is relatively independent of
2. K. Schweizerhoof and P. Wriggen, Consistent linear-
arc length procedures, so the choice of method is ization for path following methods in nonlinear FE
dependent on the desired overait convergence charac- analysis. Comput. Meth. appl. Mech. Engng. 59.261-279
teristics. For the solution of highly nonlinear prob- (1986).
3. M. A. Crisfield, A fast incremental/iterative solution
lems, with severe local alteration of stiffness, the most
procedure that handles snap through. Compur. Srr~cr.
robust arc length procedures are Crisfield’s method 13, 55-62 (1981).
and the new simplified method for ‘explicit spherical 4. B. W. R. Forde, Iteration procedures for sudden local
iteration’. Pure Newton or BFGS quasi Newton alteration of structural stiffness. Mitteilung Nr. 6, Insti-
updates can be used in conjunction with these tut fur Baustatik der Unive~i~t Stuttgart, Professor
Dr.-fng., E. Ramm (1986).
methods to improve convergence velocity and com- 5. M. Papadrakakis, Post-buckling analysis of spatial
putational efficiency. The other arc length procedures structures by vector iteration methods. Compur. Srrucr.
can be used for the solution of all problems where 14, 393-402 (1981).
convergency is not severely strained. 6. E. Riks, The application of Newton’s method to the
problem of elastic stability. J. appl. Neck 39,
10%1066 (1972).
REFERENCES 7. G. A. Wempner, Discrete approximations related to
I. E. Ramm, Strategies for tracing the nonlinear response nonlinear theories of solids. Inr. j. Sofidr Srruct. 7,
near limit points. In Europe-U.S. Workshop, nonlinear 1581-1599 (1971).