- Email: [email protected]
Orthogonal trajectory accession to the nonlinear equilibrium curve
COMPUTER METHODS NORTH-HOLLAND
IN APPLIED
MECHANICS
AND ENGINEERING
47 (1984) 283-297
ORTHOGONAL TRAJECTORY ACCESSION TO THE NONLINEAR EQUILIBRIUM CURVE
Isaac FRIED Department of Mathematics, Boston University, Boston, MA 02215, U.S.A. Received 12 December 1983 Revised manuscript received 14 March 1984
An alternative to the Riks-Wempner-Crisfield iterative correction scheme is presented that does not require an explicit displacement-load accession path to the nonlinear equilibrium curve, nor a known equilibrium point. Its symmetry with respect to the displacement and load assures success in rounding limit points as well as turning points.
1. Introduction
Consider the single, implicit equilibrium
curve
r(x, A)= 0
(1)
in which x is the displacement and A and load. To trace A versus x we shall have to compute close root pairs (x, A) to (1). Since iterative methods must invariably be used to solve the nonlinear t(~, A) = 0, two strategy decisions concerning the solution procedure have to be made: (i) how to advance from an established equilibrium point, say (x0, ho), to a next initial guess, say (x1, A& and (ii) what corrective method to use if r(xi, A1) is untolerably large. The simplest and cheapest tracking, or continuation, procedure for r(x, A) = 0 consists of increasing A0 to Ai for a new initial guess (x1 = x0, Al) and a (modified) Newton-Raphson iterative solution of r(x, Ai), with a consman? load Al. That a monotonic A sequence can miss sections of the equilibrium curve has long been pointed out [l, 21. To remedy this, Riks [3,5], Wempner [4], Crisfield [6,7] and others [8,9] suggest an iterative approach to the equilibrium curve on a variable A = A(x) load-displacement constraint. The difference between the different methods being that while Riks and Wempner advocate a linear (planar) constraint, Crisfield promotes a spherical one. Continuation [lo, 111 of the equilibrium curve is based on linearization. Let (x,, A,) be a point in the (x, A) plane not necessarily satisfying the equilibrium equation, r(x,, A,) Z 0. Let further Sx = x”+~- x, and SA = A,+i - A, denote the displacement and load corrections. Linearization of r(x, + Sx, A, + 6,) = 0 yields r, + (,6x + i,SA = 0 0045-7825/84/$3.00 @ 1984, Elsevier Science Publishers B.V. (North-Holland)
(2)
284
I. Fried, On nonlinear equilibrium curves
where I’ = dr/ax and i = arlah. Henceforth we shall omit the subscript II when referring to the current nth point. As for the prediction, that is, for the move from the already computed equilibrium point (x0, A,) to the next initial guess (x,, A,), it is commonly agreed to supplement (2) with the elliptical constraint 8x2+ Sh2= s2
(3)
where s is the step size. Combining 6ho = Al - A0 = k
(2) and (3) we get
‘I’
sx() = x, -
(i’ + rr2Y2 ’
x0
=
si
+
V’ + w2
since r(xO, ho) = 0. The choice of sign in (4) determines the direction of travel. But instead of a tangential prediction we may use the chord. For this we need two equilibrium points, say (x_~, Adl) and (x0, A,,), and we linearly predict that x1 = x0(1 + CX)- x-lcu )
Al = Ao(1 + a) - A_lc-u,
(9
which with (3) produces a =
s[(xo- x-J2 + (ho - A_,)2]-“2 .
(6)
In this paper we present and experiment with a correction procedure, based on the Newton-Raphson method as expressed in (2) whereby the initial guess (xi, A,) is iteratively improved on an orthogonal trujectory to the equilibrium curve, but with no analytical expression for it and without the involvement of the previously computed (x0, Ao).
2. The Riks-Wempner
method
The cleverness of this method lies in the replacement linear relationship crx+pA = y
of the constant A correction
by the
(7)
so that a8x + /?6A = 0. Adding this constraint to (2) yields SA = y+
3
6X=-&,
or if the constraint line in (7) is chosen to be orthogonal Riks and Wempner suggest, then SA =
r.SA: .’ sxo r
6x = -
r.Sx . r’- r$
to the tangent vector (8x,, aho), as
(9)
285
I. Fried, On nonlinear equilibrium curves
In preparation rewritten as
for things to come we further observe that since rhSxO+ iOoSho = 0, (9) may be
If T(X,A) = f(x)-
A, then r’ = f’, i = -1, and (10) becomes
that is shown to be directly obtainable from the Newton-Raphson method applied to the intersection of r = T(X,A) and A = -x/f I,+ j3, Indeed, since Sx = -r/(dr/dx), we have that 6x = --r/v’- A’), but A’ = -l/f;, and hence (11). For the success and failure of the Riks-Wempner method in turning a limit point see Fig. 1.
3. The Crisfield method Crisfield replaces the linear constraint in (7) by the circular (elliptical) (x - &J2 + (A - hJ2 =
s2.
(12)
Using the linearization r + r’Sx + 64 = 0 and the constraint (x + Sx - xJ2 + (A + SA - A0)2= s2, $A is eliminated and we are left with the quadratic equation [12] 8x2@” + i’) + 284~’ + (x - x0)i2 - (A - A&‘] + r2 - 2(A - A&i = 0
(13)
for Sx. The load increment 6A is obtained either from r + r’Sx + %A = 0 or directly from the circular constraint. A direct Newton-Raphson application to r = T(X,A), (x - x0)2+ (A - ho)’ = s2 produces
Fig. 1. Success (a) and failure
(b) of the linear constraint
method.
I. Fried, On nonlinear equilibn’um curves
286 A
I
0
0 I
fis
a
b X
Fig. 2. Success
Sx = -r/(dr/dx),
(a) and failure
(b) of the circular
constraint
method.
drldx = r’ + i dA/dx, (x - x0) dx + (A - h,) dA = 0, and
6x=-
;_,
(14)
r~_r--d2
A - A0
or alternatively 6A=-
;_*.
(15)
r-r!2
x-
x0
It is interesting to compare (14) (15) with (9). Note that (14) and (15) do not correspond to Crisfield’s method. For the success and failure of the Crisfield method in rounding limit and turning points see Fig. 2. 4. Orthogonal trajectory What we consider now is an orthogonal trajectory approach to the equilibrium curve without the need for an explicit constraint, and without the involvement of an equilibrium point (x0, A,). All we have to do for that is apply the Newton-Raphson method to r = r(x, A) = 0 so as to have Sx = -r/(dr/dx), add to this dr = r’ dx + i dh, and impose the orthogonality condition dh r’ --_= dxi
1
(16)
’
and we get
SA=-gp,
~x=-L
rr2 + i’*
(17)
I, Fried, On nonlinear equilibrium curves
287
which is symmetric in x and A. Comparing this with (10) for the Riks-Wempner method reveals how the orthogonal constraint is updated here. In fact, Ramm 1131, and apparently also Haselgrove [21], suggested a corrections method that includes updated normals. When applied to-r(x, A) = 0, Ramm’s method produces the load increment SA = -
(18)
i - r’(A1 - rh,),(X - X0)’
involving not only the previous equilibrium point (x0, Ao) but also the new initial guess (x1, Al). It is interesting to compare (18) with (9) and (15).
5. The vector form
Our main interest is in systems of nonlinear stiffness equation as produced by finite elements applied to elastic problems, and we turn our attention in this section to the application of the orthogonal trajectory method to problems with many degrees of freedom. Let T&X,A) be the total potential energy of the discretization, with x being the displacement vector and h the scalar load. Here r(x, A) = &r/ax is the gradient vector of T. If (x0, A,,) is, again, an equilibrium point so that r(xo, A,) = 0, then a linear expansion around it is written as (19) where ar/ax is the stiffness matrix, say K, and where adah is the (negative) load vector, say p. For how to compute K and p see [14-191. A prediction with step size s subject to the constraint (xl - xo)f(xl - x,,) + (A1 - Ao)~ = s2,
(20)
together with the linearization
(21)
Ko(xl-xo)+(Al-Ao)po=O, leads to the initial guess Al = A,,+ ~(1 + p;K,‘K,‘po)-1’2, where SAo= A, - Ao. A chordial prediction a =
x1 =
x,, - SA,K,'po
(22)
is here essentially the same as in (5) and (6),
s[(xo- xmI)‘(xo- x_~)+ (A0 -
A_1)2]-1'2
(23)
Al = A,,(1+ a) - A_la,
(24)
and x1 = x0(1 + a) - X-I(Y) with no sign ambiguity.
288
I. Fried, On nonlinear equilibrium curves
Newton-Raphson’s 6x =
correction applied to the system of linear equations r(x, A) = 0 becomes
-K-l@ + SAP).
(25)
A constant load iteration with 6A = 0 reduces (25) to SX = -K’r. constraint that relates Sh to 6x by
A linear load-displacement
SA = at8x
(26)
with the given vector a produces the correction a’K_‘r
SA = -
1 + a’K-‘p
For orthogonal
trajectory
’
Sx = -K-l@ + hip),
(27)
accession we set a = Kelp, and have 6x = - K-‘(r + 6Ap).
(28)
Observe that both the Riks-Wempner method and the present one call for the computation of K-‘r and K-‘p, and otherwise entail comparable computational effort.
6. Midpoint correction A Runge-Kutta type correction in the form of a midpoint reevaluation of the stiffness matrix and load vector is also often used [20] in nonlinear computational elasticity. In this method an initial guess is predicted along the tangent to the equilibrium curve with (20) or (24), but then K and p are reformed at x = x0 + $8~ and A = A0+ $A, and the new data are resubstituted into (20). Each step requires in this way two assemblies and two inversions. An alternative of equal computational effort would be a linear prediction with ol~e Newton-Raphson correction. It happens that this procedure is considerably more accurate than a midpoint correction. It is sufficient that we analyze the methods on the parabolic stiffness equation A = crx’, where CY= $A”, with the circular constraint x2 + A2 = s*. Their exact intersection point (x,, A,) is at x, = k[-? + i(l + 4a*syy*
)
A, = &[-
1 + (1
+ 4a*s*)“*] .
(29)
Without correction point p in Fig. 3, (x = s, A = 0) constitutes the approximation to equilibrium. Since both the exact solution (x,, A,) and all other approximations are at an equal distance s from the origin we prefer to use the directional error cp = cos
cxxe :*^^‘J
-l
(30)
289
I. Fried, On nonlinear equilibrium curves
s/2
0 Fig. 3. Midpoint
of (x, A). For point
p
S
and Newton-Raphson
corrections
we have that cp = cos-‘(x,/s),
we get that for tangential prediction
constraint.
or with the expansions
A, = &(l
X,=S(l-_t(yZS2+~(y4S4?...),
with circular
-
(y*s*+ 2cu4S4T . . . )
(31)
only
cp = cys
(32)
if s is small. Midpoint correction reaches point m of Fig. 3 with
cY2s2
s x
=
(1
+
a2s2)l/2
)
*
=
(1
+
a2s2)l/2
(33)
’
Using (31) and with similar expansions in powers of s of x and A in (33) we have for small s that cp=
&3S3
(34)
which is a substantial improvement over (32). One Newton-Raphson correction with a circular constraint leads to 6x = -2cu2s3 - s + s(1 + 3a2s2)“2 1+4a2s2 ’ and a directional
.X=,X+&X,
h = as2 +
2CkS6x
(35)
error
(p = &# which is two orders better than (34).
(36)
290
I. Fried, On nonlinear equilibrium curves
7. Numerical examples The ability of orthogonal trajectory accession to succeed around limit points, bifurcation points, turning points, snap-throughs and detachments is numerically demonstrated on two simple structural elastic problems. First to be considered is the two-member hinged truss shown in Fig. 5 in its original undeformed state, As the load A is increased the two elastic bars shorten symmetrically until a snap-through occurs and the triangle is inverted. A displacement x causes the strain e= 1-(1+~~-2hx)~‘~
(37)
in the bars. With each bar being of the stiffness $k the total potential system is n(x) = ke2 - hx, and hence r(x, A) = kee’ - A ,
energy of the two bar
r’ = k(e’2 + ee”)
i=-1,
(38)
where h-x
e’= l_ey The exact equilibrium
I
l
starf
err=
-l+ef2 l_e
-
curve r(x, A) = 0 is shown in Fig. 4.
Iexact 1.0 b
I *.5-
:.
s=o. 05
k=25
h=0.5
NR=3
‘,. ‘4
-0.5
Fig. 4. Orthogonal trajectory convergence to the equilibrium curve of the truss in Fig. 5 from different starting points.
-
:
S
:
Fig. 5. Pointwise tracking of the equilibrium curve with a tangential corrector and orthogonal trajectory corrections.
I. Fried, On nonlinear equilibrium curves
291
Fig. 4 shows also the convergence paths of orthogonal trajectory accession given by (17) for different starting points that are far from equilibrium. To observe the orthogonal trajectory one must look close to the end of the iterative approach. A starting point on the normal to a limit point converges in one step. Notice also the special attraction of the limit points. Fig. 5 shows the point by point continuation of the equilibrium curve of Fig. 4 with the predictor of (4) and the corrector of (17), with 3 Newton-Raphson corrections per step. The two-member truss shown in Fig. 6 is flat when A = 0. Its behavior under load includes bifurcation, turning points and detachment. Two degrees of freedom, x and 4, determine the configuration of the deformed system, whose total potential energy reads now r(x, cp) = 2k,(p2 + kle2 - Ax
(40)
where k, is the stiffness of each bar, k2 the stiffness of the torsional spring at the vertex, and e the strain in the bars,
,+1-x
(41)
cos cp -
The two stiffness equations
are obtained from n(x, 4) through differentiation, %=2k,e$--A=O.
!$=4kZg+2kle$-,
If k, = $ and k2 = &E, then after eliminating
(42)
x between the pair of (42) we are left with
(43)
: :
>=I . .
.
.‘.e=1.023~ . ,e=l.l027 .'., . . . .. ‘. '.
'. . . . . . . . . ...* ,_ . : : : : : ::::::.:::::: r=0.9 ". . . . . . . . .
"
:
k, =I/2 kr =e/16
1 0
1
~18
Fig. 6. Computed
T/4 equilibrium
points
I
I
3rr/8
n-/2
for the two-member
truss.
+
292
I. Fried, On nonlinear equilibrium curves
When cp = 0 we have either A = x or A = g1 k (1 - #‘2] )
(44)
so that for E < 1 two bifurcation points occur on the A axis, and they coalesce for E = 1. As E surpasses 1 the nonlinear bifurcation curve detaches from the A axis and a turning point (i.e., dA/dx = a) is created. The critical cp* and A* at which turning occurs are given by tan +* +* -e,
1 h*=2cos$*.
(45)
For a given E > 1 no nonlinear solution exists for cp< cp”, whatever A. Tangential prediction plus orthogonal trajectory accession has no difficulty following the equilibrium curves past the turning points, as can be seen in Fig. 6. The two values of E = 1.0235 and E = 1.1027 in Fig. 6 correspond to q* = &n and cp* = HIT,respectively. Fig. 7 uses C$and A from Fig. 6 to trace A resp. x using x= 1+cos~(Acos~-1).
(46)
Finally, we compare the performance of the various continuation methods discussed in this paper in rounding a limit point (dA/dx = 0) and a turning point (dh/dx = ~0). The parabolic I = ax(1 - x) - A = 0 has a limit point at x = 4, A = $(Yand is of the shape of Fig. 1. We choose (Y= 8, x0 = 0.4375 (at which r’ = l), and A0= 1.96875. A tangential predictor of step size s = 0.2 lands us at x1 = 0.57892136 and AI = 2.1101714. This is approximately the situation shown in Fig. l(b) and Fig. 2(a). Both the linear constraint method and Ramm’s method fail to find an equilibrium point starting with these x1 and AI. Only the circular constraint method and the l.O- h
l=l.:.E=O.g l
bifurcation
X I
0
0.5 Fig. 7. As in Fig. 6 but x vs. A.
I
1.0
293
I. Fried, On nonlinear equilibrium curves Table 1 Convergence r=8x(l-x)-h Step
of the circular
Circular
constraint
constraint
s = 0.2
X
0 1 2 3 4
and orthogonal
Orthogonal
A
0.57892136 0.62723680 0.62287506 0.62283738 0.62283737
2.1101714 1.8891614 1.8793660 1.8792878 1.8792878
trajectory
upon
trajectory
X
0.57892136 0.50104992 0.50023522 0.50023520 0.50023520
methods
A
2.1101714 2.0485029 2.0000049 1.9999996 1.9999996
orthogonal trajectory method converge here-to different points, though, as can be seen in Table 1. Reversing the roles of A and X: r = mh (1 - A) - x creates an equilibrium curve as in Fig. 2(b) with a turning point. Starting from the reversed AI = 0.57892136 and x1 = 2.1101714, the circular constraint method fails, but the orthogonal trajectory method, because of its inherent symmetry to x and A, converges to the reversed values of x and A in Table 1.
8. Bending of a rod with a yielding support The paper will remain incomplete without a finite element multi degrees of freedom example. For this purpose we shall apply the predictor-corrector method of Section 5 to follow the large deformation of a unit elastic rod shown in Fig. 9. The rod is not built-in but is rather
2
1
M/a M/a=+((l+
I~/p)2)/(li(#/p)4))p’2
I. 3 I. I 0.9 0.7 0.5 0.3 P=O.l 4
0
I Fig. 8. Constitutive
2 behavior
of a yielding
spring.
I. Fried, On nonlinear equilibrium curves
294
hinged on a yielding torsional spring with the constitutive ship
moment-twist,
M = M(4), relation-
1 t (c#J/@)’ d2 (r#J/py>
M=QI4 (1t
(47)
where (Y,j3 and p are parameters. Fig. 8 shows M/a versus C$ for different values of /3 and 0.9. We use this value of p in all the following computations. The total potential energy of the inextensible, tip loaded, unit rod can be written in terms of the slope 0 only as
p =
?T@) = 01(;(Y2- P sin 8 t Q cos 6) ds + E(4) I
(48)
where s is arch length and E(4) the elastic energy stored in the root spring due to a rotation C#J = e(0). Following [15] we discretize the rod with quadratic elements and integrate the element total potential energy n= with a twcD-1point Gauss method. For a typical element of size h we have
3
I
0 Fig. 9. Tip
deflection
vs. load in a rod with a built-in
I (a) and a yielding spring (b) support.
295
I. Fried, On nonlinear equilibrium curves
that ne =
$h i 4h-28;j=I
where j is the integration ei = ehf+,
P sin 0, + Q cos 0,
(49)
index and bj
=
(50)
fl$j
where 0: = (e,, &, 0,) is the element nodal values vector and 4: = :(l+ 4: =
4’2 = k(l -
V/3,4,1 - V/3))
$(-2A&
3,4V\/3, -2X6+
3))
&
=
G,
4,1+ V/3) )
$(2V/3 - 3, -4d/3,2V\/3 + 3) .
The element gradient r, and element stiffness matrix k, are obtained repeated integrations with respect to &,
from rre in (49) by
r,=~=~2h-‘B,~i-fh(Pcos8,+QSinB,)k, e j=l k,
=
3
=
c
i
2h-‘&d;
-
$
(51)
(52)
h(Q cos t3j- P sin Oi)&#j .
(53)
j=l
Since we plan to keep the tip force P constant and vary Q only we derive the element vector corresponding to the global p in (25) from
pe
pc = g
= $ - 4h sin e&j. I-1
Fig. 10. Computed
equilibrium
configurations
of a rod with a yielding
support.
296
I. Fried, On nonlinear equilibrium curves
Q/Qcr
0
Fig. 11. Lowest
eigenvalue
A? of the global stiffness
matrix
K VS. Q/Q,, for rod (b) in Fig. 9.
The element vectors pe and r, are routinely assembled into the global p and r, and k, is likewise assembled into the global stiffness matrix K. To account for the root spring we add M = dE/&$ of (47) to the first entry of r that corresponds to the degree of freedom C$= O(O), and dM/&j = d*E/d4* to the first diagonal term of K. In the computations we normalize P and Q with respect to Q,, = $rr*, choose the small P = 0.01 Q,, to avoid the symmetry, and take the root spring with CY= 10, p = 0.2 and p = 0.9. To start the predictor-corrector algorithm of (23) (24) and (28) we compute the two equilibrium configurations for Q = 0 and Q/Q,, = 0.1 and proceed cyclically with the step size s = 0.1. Fig. 9 traces the thus computed y(l) as a function of Q/Q,, for the built-in rod (a), and for the rod with a yielding root (b). Fig. 10 shows the computed bent configuration of the rod at the end of each cycle, while Fig. 11 shows the variation with Q/Q,, of the lowest eigenvalue h? of the global stiffness matrix K at equilibrium for the rod (b) in Fig. 9. All discretizations are done with 7 quadratic elements.
References J.H. Argyris, Continua and discontinua in matrix methods of structural analysis, AFFDL-TR-66-80 (1966) 11. G.A. Thurston, Continuation of Newton’s method through bifurcation points, J. Appl. Mech. 36 (1969) 425430. of Newton’s method to the problem of elastic stability, J. Appl. Mech. 39 (1972) [31 E. Riks, The application 1060-1066. related to nonlinear theories of solids, Internat. J. Solids Structures 7 141 G.A. Wempner, Discrete approximations (1971) 1581-1599. approach to the solution of snapping and buckling problems, Internat. J. Solids [51E. Riks, An incremental Structures 15 (1979) 329-551.
[ll PI
I. Fried, On nonlinear equilibrium curves
297
[6] M.A. Crisfield, A fast incremental/iterative solution procedure that handles “snap-through”, Comput. & Structures 13 (1981) X-62. [7] M.A. Crisfield, An arc-length method including line searches and accelerations, Internat. J. Numer. Meths. Engrg. 19 (1983) 1269-1289. [S] J. Padorean and S. Tovichakchaikul, Self adaptive predictor-corrector algorithms for static nonlinear structural analysis, Comput. & Structures 15 (1982) 365-377. [9] J. Padorean and S. Tovichakchaikul, Operating characteristics of hyperbolically and elliptically constrained self-adaptive incremental Newton-Raphson algorithms, J. Franklin Institute 316 (1983) 197-223. [lo] E.L. Allgower and K. George, Simplicial and continuation methods for approximating fixed points and solutions to systems of equations, SIAM Rev. 22 (1980) 28-85. [ll] E.L. Allgower, A survey of homotopy methods for smooth mappings, in: E.L. Allgower et al., eds., Numerical Solution of Nonlinear Equations (Springer, Berlin, 1981). [12] L.T. Watson and S.M. Holzer, Quadratic convergence of Crisfield’s method, Comput. & Structures 17 (1983) 69-72. [13] E. Ramm, Strategies for tracing the nonlinear response near limit points, in: W. Wunderlich et al., eds., Non-linear Finite Element Analysis in Structural Mechanics (Springer, Berlin, 1981). [14] I. Fried, Numerical Solution of Differential Equations (Academic Press, New York 1979). [15] I. Fried, Stability and equilibrium of the straight and curved elastica-finite element computation, Comput. Meths. Appl. Mech. Engrg. 28 (1981) 49-61. [16] I. Fried, Nonlinear finite element computation of the equilibrium and stability of the circular plate, Internat. J. Numer. Meths. Engrg. 17 (1981) 1427-1440. [17] I. Fried, Finite element computation of large rubber membrane deformations, Internat. J. Numer. M&c_. Engrg. 18 (1982) 653-660. [18] I. Fried, Nonlinear finite element computation of the equilibrium, stability and motion of the extensional beam and ring, Comput. Meths. Appl. Mech. Engrg. 38 (1983) 29-44. [19] I. Fried, Large deformation static and dynamic finite element analysis of extensible cables, Comput. & Structures 15 (1982) 315-319. [20] T.Y. Yang, Matrix displacement solution Structures 9 (1973) 829-842. [21] C.B. Haselgrove, The solution of non-linear conditions, Comput. J. 4 (1961) 255-259.
to elastica equations
problems
of beams
and of differential
and
equations
frames,
Internat.
with two-point
J. Solids boundary
IN APPLIED
MECHANICS
AND ENGINEERING
47 (1984) 283-297
ORTHOGONAL TRAJECTORY ACCESSION TO THE NONLINEAR EQUILIBRIUM CURVE
Isaac FRIED Department of Mathematics, Boston University, Boston, MA 02215, U.S.A. Received 12 December 1983 Revised manuscript received 14 March 1984
An alternative to the Riks-Wempner-Crisfield iterative correction scheme is presented that does not require an explicit displacement-load accession path to the nonlinear equilibrium curve, nor a known equilibrium point. Its symmetry with respect to the displacement and load assures success in rounding limit points as well as turning points.
1. Introduction
Consider the single, implicit equilibrium
curve
r(x, A)= 0
(1)
in which x is the displacement and A and load. To trace A versus x we shall have to compute close root pairs (x, A) to (1). Since iterative methods must invariably be used to solve the nonlinear t(~, A) = 0, two strategy decisions concerning the solution procedure have to be made: (i) how to advance from an established equilibrium point, say (x0, ho), to a next initial guess, say (x1, A& and (ii) what corrective method to use if r(xi, A1) is untolerably large. The simplest and cheapest tracking, or continuation, procedure for r(x, A) = 0 consists of increasing A0 to Ai for a new initial guess (x1 = x0, Al) and a (modified) Newton-Raphson iterative solution of r(x, Ai), with a consman? load Al. That a monotonic A sequence can miss sections of the equilibrium curve has long been pointed out [l, 21. To remedy this, Riks [3,5], Wempner [4], Crisfield [6,7] and others [8,9] suggest an iterative approach to the equilibrium curve on a variable A = A(x) load-displacement constraint. The difference between the different methods being that while Riks and Wempner advocate a linear (planar) constraint, Crisfield promotes a spherical one. Continuation [lo, 111 of the equilibrium curve is based on linearization. Let (x,, A,) be a point in the (x, A) plane not necessarily satisfying the equilibrium equation, r(x,, A,) Z 0. Let further Sx = x”+~- x, and SA = A,+i - A, denote the displacement and load corrections. Linearization of r(x, + Sx, A, + 6,) = 0 yields r, + (,6x + i,SA = 0 0045-7825/84/$3.00 @ 1984, Elsevier Science Publishers B.V. (North-Holland)
(2)
284
I. Fried, On nonlinear equilibrium curves
where I’ = dr/ax and i = arlah. Henceforth we shall omit the subscript II when referring to the current nth point. As for the prediction, that is, for the move from the already computed equilibrium point (x0, A,) to the next initial guess (x,, A,), it is commonly agreed to supplement (2) with the elliptical constraint 8x2+ Sh2= s2
(3)
where s is the step size. Combining 6ho = Al - A0 = k
(2) and (3) we get
‘I’
sx() = x, -
(i’ + rr2Y2 ’
x0
=
si
+
V’ + w2
since r(xO, ho) = 0. The choice of sign in (4) determines the direction of travel. But instead of a tangential prediction we may use the chord. For this we need two equilibrium points, say (x_~, Adl) and (x0, A,,), and we linearly predict that x1 = x0(1 + CX)- x-lcu )
Al = Ao(1 + a) - A_lc-u,
(9
which with (3) produces a =
s[(xo- x-J2 + (ho - A_,)2]-“2 .
(6)
In this paper we present and experiment with a correction procedure, based on the Newton-Raphson method as expressed in (2) whereby the initial guess (xi, A,) is iteratively improved on an orthogonal trujectory to the equilibrium curve, but with no analytical expression for it and without the involvement of the previously computed (x0, Ao).
2. The Riks-Wempner
method
The cleverness of this method lies in the replacement linear relationship crx+pA = y
of the constant A correction
by the
(7)
so that a8x + /?6A = 0. Adding this constraint to (2) yields SA = y+
3
6X=-&,
or if the constraint line in (7) is chosen to be orthogonal Riks and Wempner suggest, then SA =
r.SA: .’ sxo r
6x = -
r.Sx . r’- r$
to the tangent vector (8x,, aho), as
(9)
285
I. Fried, On nonlinear equilibrium curves
In preparation rewritten as
for things to come we further observe that since rhSxO+ iOoSho = 0, (9) may be
If T(X,A) = f(x)-
A, then r’ = f’, i = -1, and (10) becomes
that is shown to be directly obtainable from the Newton-Raphson method applied to the intersection of r = T(X,A) and A = -x/f I,+ j3, Indeed, since Sx = -r/(dr/dx), we have that 6x = --r/v’- A’), but A’ = -l/f;, and hence (11). For the success and failure of the Riks-Wempner method in turning a limit point see Fig. 1.
3. The Crisfield method Crisfield replaces the linear constraint in (7) by the circular (elliptical) (x - &J2 + (A - hJ2 =
s2.
(12)
Using the linearization r + r’Sx + 64 = 0 and the constraint (x + Sx - xJ2 + (A + SA - A0)2= s2, $A is eliminated and we are left with the quadratic equation [12] 8x2@” + i’) + 284~’ + (x - x0)i2 - (A - A&‘] + r2 - 2(A - A&i = 0
(13)
for Sx. The load increment 6A is obtained either from r + r’Sx + %A = 0 or directly from the circular constraint. A direct Newton-Raphson application to r = T(X,A), (x - x0)2+ (A - ho)’ = s2 produces
Fig. 1. Success (a) and failure
(b) of the linear constraint
method.
I. Fried, On nonlinear equilibn’um curves
286 A
I
0
0 I
fis
a
b X
Fig. 2. Success
Sx = -r/(dr/dx),
(a) and failure
(b) of the circular
constraint
method.
drldx = r’ + i dA/dx, (x - x0) dx + (A - h,) dA = 0, and
6x=-
;_,
(14)
r~_r--d2
A - A0
or alternatively 6A=-
;_*.
(15)
r-r!2
x-
x0
It is interesting to compare (14) (15) with (9). Note that (14) and (15) do not correspond to Crisfield’s method. For the success and failure of the Crisfield method in rounding limit and turning points see Fig. 2. 4. Orthogonal trajectory What we consider now is an orthogonal trajectory approach to the equilibrium curve without the need for an explicit constraint, and without the involvement of an equilibrium point (x0, A,). All we have to do for that is apply the Newton-Raphson method to r = r(x, A) = 0 so as to have Sx = -r/(dr/dx), add to this dr = r’ dx + i dh, and impose the orthogonality condition dh r’ --_= dxi
1
(16)
’
and we get
SA=-gp,
~x=-L
rr2 + i’*
(17)
I, Fried, On nonlinear equilibrium curves
287
which is symmetric in x and A. Comparing this with (10) for the Riks-Wempner method reveals how the orthogonal constraint is updated here. In fact, Ramm 1131, and apparently also Haselgrove [21], suggested a corrections method that includes updated normals. When applied to-r(x, A) = 0, Ramm’s method produces the load increment SA = -
(18)
i - r’(A1 - rh,),(X - X0)’
involving not only the previous equilibrium point (x0, Ao) but also the new initial guess (x1, Al). It is interesting to compare (18) with (9) and (15).
5. The vector form
Our main interest is in systems of nonlinear stiffness equation as produced by finite elements applied to elastic problems, and we turn our attention in this section to the application of the orthogonal trajectory method to problems with many degrees of freedom. Let T&X,A) be the total potential energy of the discretization, with x being the displacement vector and h the scalar load. Here r(x, A) = &r/ax is the gradient vector of T. If (x0, A,,) is, again, an equilibrium point so that r(xo, A,) = 0, then a linear expansion around it is written as (19) where ar/ax is the stiffness matrix, say K, and where adah is the (negative) load vector, say p. For how to compute K and p see [14-191. A prediction with step size s subject to the constraint (xl - xo)f(xl - x,,) + (A1 - Ao)~ = s2,
(20)
together with the linearization
(21)
Ko(xl-xo)+(Al-Ao)po=O, leads to the initial guess Al = A,,+ ~(1 + p;K,‘K,‘po)-1’2, where SAo= A, - Ao. A chordial prediction a =
x1 =
x,, - SA,K,'po
(22)
is here essentially the same as in (5) and (6),
s[(xo- xmI)‘(xo- x_~)+ (A0 -
A_1)2]-1'2
(23)
Al = A,,(1+ a) - A_la,
(24)
and x1 = x0(1 + a) - X-I(Y) with no sign ambiguity.
288
I. Fried, On nonlinear equilibrium curves
Newton-Raphson’s 6x =
correction applied to the system of linear equations r(x, A) = 0 becomes
-K-l@ + SAP).
(25)
A constant load iteration with 6A = 0 reduces (25) to SX = -K’r. constraint that relates Sh to 6x by
A linear load-displacement
SA = at8x
(26)
with the given vector a produces the correction a’K_‘r
SA = -
1 + a’K-‘p
For orthogonal
trajectory
’
Sx = -K-l@ + hip),
(27)
accession we set a = Kelp, and have 6x = - K-‘(r + 6Ap).
(28)
Observe that both the Riks-Wempner method and the present one call for the computation of K-‘r and K-‘p, and otherwise entail comparable computational effort.
6. Midpoint correction A Runge-Kutta type correction in the form of a midpoint reevaluation of the stiffness matrix and load vector is also often used [20] in nonlinear computational elasticity. In this method an initial guess is predicted along the tangent to the equilibrium curve with (20) or (24), but then K and p are reformed at x = x0 + $8~ and A = A0+ $A, and the new data are resubstituted into (20). Each step requires in this way two assemblies and two inversions. An alternative of equal computational effort would be a linear prediction with ol~e Newton-Raphson correction. It happens that this procedure is considerably more accurate than a midpoint correction. It is sufficient that we analyze the methods on the parabolic stiffness equation A = crx’, where CY= $A”, with the circular constraint x2 + A2 = s*. Their exact intersection point (x,, A,) is at x, = k[-? + i(l + 4a*syy*
)
A, = &[-
1 + (1
+ 4a*s*)“*] .
(29)
Without correction point p in Fig. 3, (x = s, A = 0) constitutes the approximation to equilibrium. Since both the exact solution (x,, A,) and all other approximations are at an equal distance s from the origin we prefer to use the directional error cp = cos
cxxe :*^^‘J
-l
(30)
289
I. Fried, On nonlinear equilibrium curves
s/2
0 Fig. 3. Midpoint
of (x, A). For point
p
S
and Newton-Raphson
corrections
we have that cp = cos-‘(x,/s),
we get that for tangential prediction
constraint.
or with the expansions
A, = &(l
X,=S(l-_t(yZS2+~(y4S4?...),
with circular
-
(y*s*+ 2cu4S4T . . . )
(31)
only
cp = cys
(32)
if s is small. Midpoint correction reaches point m of Fig. 3 with
cY2s2
s x
=
(1
+
a2s2)l/2
)
*
=
(1
+
a2s2)l/2
(33)
’
Using (31) and with similar expansions in powers of s of x and A in (33) we have for small s that cp=
&3S3
(34)
which is a substantial improvement over (32). One Newton-Raphson correction with a circular constraint leads to 6x = -2cu2s3 - s + s(1 + 3a2s2)“2 1+4a2s2 ’ and a directional
.X=,X+&X,
h = as2 +
2CkS6x
(35)
error
(p = &# which is two orders better than (34).
(36)
290
I. Fried, On nonlinear equilibrium curves
7. Numerical examples The ability of orthogonal trajectory accession to succeed around limit points, bifurcation points, turning points, snap-throughs and detachments is numerically demonstrated on two simple structural elastic problems. First to be considered is the two-member hinged truss shown in Fig. 5 in its original undeformed state, As the load A is increased the two elastic bars shorten symmetrically until a snap-through occurs and the triangle is inverted. A displacement x causes the strain e= 1-(1+~~-2hx)~‘~
(37)
in the bars. With each bar being of the stiffness $k the total potential system is n(x) = ke2 - hx, and hence r(x, A) = kee’ - A ,
energy of the two bar
r’ = k(e’2 + ee”)
i=-1,
(38)
where h-x
e’= l_ey The exact equilibrium
I
l
starf
err=
-l+ef2 l_e
-
curve r(x, A) = 0 is shown in Fig. 4.
Iexact 1.0 b
I *.5-
:.
s=o. 05
k=25
h=0.5
NR=3
‘,. ‘4
-0.5
Fig. 4. Orthogonal trajectory convergence to the equilibrium curve of the truss in Fig. 5 from different starting points.
-
:
S
:
Fig. 5. Pointwise tracking of the equilibrium curve with a tangential corrector and orthogonal trajectory corrections.
I. Fried, On nonlinear equilibrium curves
291
Fig. 4 shows also the convergence paths of orthogonal trajectory accession given by (17) for different starting points that are far from equilibrium. To observe the orthogonal trajectory one must look close to the end of the iterative approach. A starting point on the normal to a limit point converges in one step. Notice also the special attraction of the limit points. Fig. 5 shows the point by point continuation of the equilibrium curve of Fig. 4 with the predictor of (4) and the corrector of (17), with 3 Newton-Raphson corrections per step. The two-member truss shown in Fig. 6 is flat when A = 0. Its behavior under load includes bifurcation, turning points and detachment. Two degrees of freedom, x and 4, determine the configuration of the deformed system, whose total potential energy reads now r(x, cp) = 2k,(p2 + kle2 - Ax
(40)
where k, is the stiffness of each bar, k2 the stiffness of the torsional spring at the vertex, and e the strain in the bars,
,+1-x
(41)
cos cp -
The two stiffness equations
are obtained from n(x, 4) through differentiation, %=2k,e$--A=O.
!$=4kZg+2kle$-,
If k, = $ and k2 = &E, then after eliminating
(42)
x between the pair of (42) we are left with
(43)
: :
>=I . .
.
.‘.e=1.023~ . ,e=l.l027 .'., . . . .. ‘. '.
'. . . . . . . . . ...* ,_ . : : : : : ::::::.:::::: r=0.9 ". . . . . . . . .
"
:
k, =I/2 kr =e/16
1 0
1
~18
Fig. 6. Computed
T/4 equilibrium
points
I
I
3rr/8
n-/2
for the two-member
truss.
+
292
I. Fried, On nonlinear equilibrium curves
When cp = 0 we have either A = x or A = g1 k (1 - #‘2] )
(44)
so that for E < 1 two bifurcation points occur on the A axis, and they coalesce for E = 1. As E surpasses 1 the nonlinear bifurcation curve detaches from the A axis and a turning point (i.e., dA/dx = a) is created. The critical cp* and A* at which turning occurs are given by tan +* +* -e,
1 h*=2cos$*.
(45)
For a given E > 1 no nonlinear solution exists for cp< cp”, whatever A. Tangential prediction plus orthogonal trajectory accession has no difficulty following the equilibrium curves past the turning points, as can be seen in Fig. 6. The two values of E = 1.0235 and E = 1.1027 in Fig. 6 correspond to q* = &n and cp* = HIT,respectively. Fig. 7 uses C$and A from Fig. 6 to trace A resp. x using x= 1+cos~(Acos~-1).
(46)
Finally, we compare the performance of the various continuation methods discussed in this paper in rounding a limit point (dA/dx = 0) and a turning point (dh/dx = ~0). The parabolic I = ax(1 - x) - A = 0 has a limit point at x = 4, A = $(Yand is of the shape of Fig. 1. We choose (Y= 8, x0 = 0.4375 (at which r’ = l), and A0= 1.96875. A tangential predictor of step size s = 0.2 lands us at x1 = 0.57892136 and AI = 2.1101714. This is approximately the situation shown in Fig. l(b) and Fig. 2(a). Both the linear constraint method and Ramm’s method fail to find an equilibrium point starting with these x1 and AI. Only the circular constraint method and the l.O- h
l=l.:.E=O.g l
bifurcation
X I
0
0.5 Fig. 7. As in Fig. 6 but x vs. A.
I
1.0
293
I. Fried, On nonlinear equilibrium curves Table 1 Convergence r=8x(l-x)-h Step
of the circular
Circular
constraint
constraint
s = 0.2
X
0 1 2 3 4
and orthogonal
Orthogonal
A
0.57892136 0.62723680 0.62287506 0.62283738 0.62283737
2.1101714 1.8891614 1.8793660 1.8792878 1.8792878
trajectory
upon
trajectory
X
0.57892136 0.50104992 0.50023522 0.50023520 0.50023520
methods
A
2.1101714 2.0485029 2.0000049 1.9999996 1.9999996
orthogonal trajectory method converge here-to different points, though, as can be seen in Table 1. Reversing the roles of A and X: r = mh (1 - A) - x creates an equilibrium curve as in Fig. 2(b) with a turning point. Starting from the reversed AI = 0.57892136 and x1 = 2.1101714, the circular constraint method fails, but the orthogonal trajectory method, because of its inherent symmetry to x and A, converges to the reversed values of x and A in Table 1.
8. Bending of a rod with a yielding support The paper will remain incomplete without a finite element multi degrees of freedom example. For this purpose we shall apply the predictor-corrector method of Section 5 to follow the large deformation of a unit elastic rod shown in Fig. 9. The rod is not built-in but is rather
2
1
M/a M/a=+((l+
I~/p)2)/(li(#/p)4))p’2
I. 3 I. I 0.9 0.7 0.5 0.3 P=O.l 4
0
I Fig. 8. Constitutive
2 behavior
of a yielding
spring.
I. Fried, On nonlinear equilibrium curves
294
hinged on a yielding torsional spring with the constitutive ship
moment-twist,
M = M(4), relation-
1 t (c#J/@)’ d2 (r#J/py>
M=QI4 (1t
(47)
where (Y,j3 and p are parameters. Fig. 8 shows M/a versus C$ for different values of /3 and 0.9. We use this value of p in all the following computations. The total potential energy of the inextensible, tip loaded, unit rod can be written in terms of the slope 0 only as
p =
?T@) = 01(;(Y2- P sin 8 t Q cos 6) ds + E(4) I
(48)
where s is arch length and E(4) the elastic energy stored in the root spring due to a rotation C#J = e(0). Following [15] we discretize the rod with quadratic elements and integrate the element total potential energy n= with a twcD-1point Gauss method. For a typical element of size h we have
3
I
0 Fig. 9. Tip
deflection
vs. load in a rod with a built-in
I (a) and a yielding spring (b) support.
295
I. Fried, On nonlinear equilibrium curves
that ne =
$h i 4h-28;j=I
where j is the integration ei = ehf+,
P sin 0, + Q cos 0,
(49)
index and bj
=
(50)
fl$j
where 0: = (e,, &, 0,) is the element nodal values vector and 4: = :(l+ 4: =
4’2 = k(l -
V/3,4,1 - V/3))
$(-2A&
3,4V\/3, -2X6+
3))
&
=
G,
4,1+ V/3) )
$(2V/3 - 3, -4d/3,2V\/3 + 3) .
The element gradient r, and element stiffness matrix k, are obtained repeated integrations with respect to &,
from rre in (49) by
r,=~=~2h-‘B,~i-fh(Pcos8,+QSinB,)k, e j=l k,
=
3
=
c
i
2h-‘&d;
-
$
(51)
(52)
h(Q cos t3j- P sin Oi)&#j .
(53)
j=l
Since we plan to keep the tip force P constant and vary Q only we derive the element vector corresponding to the global p in (25) from
pe
pc = g
= $ - 4h sin e&j. I-1
Fig. 10. Computed
equilibrium
configurations
of a rod with a yielding
support.
296
I. Fried, On nonlinear equilibrium curves
Q/Qcr
0
Fig. 11. Lowest
eigenvalue
A? of the global stiffness
matrix
K VS. Q/Q,, for rod (b) in Fig. 9.
The element vectors pe and r, are routinely assembled into the global p and r, and k, is likewise assembled into the global stiffness matrix K. To account for the root spring we add M = dE/&$ of (47) to the first entry of r that corresponds to the degree of freedom C$= O(O), and dM/&j = d*E/d4* to the first diagonal term of K. In the computations we normalize P and Q with respect to Q,, = $rr*, choose the small P = 0.01 Q,, to avoid the symmetry, and take the root spring with CY= 10, p = 0.2 and p = 0.9. To start the predictor-corrector algorithm of (23) (24) and (28) we compute the two equilibrium configurations for Q = 0 and Q/Q,, = 0.1 and proceed cyclically with the step size s = 0.1. Fig. 9 traces the thus computed y(l) as a function of Q/Q,, for the built-in rod (a), and for the rod with a yielding root (b). Fig. 10 shows the computed bent configuration of the rod at the end of each cycle, while Fig. 11 shows the variation with Q/Q,, of the lowest eigenvalue h? of the global stiffness matrix K at equilibrium for the rod (b) in Fig. 9. All discretizations are done with 7 quadratic elements.
References J.H. Argyris, Continua and discontinua in matrix methods of structural analysis, AFFDL-TR-66-80 (1966) 11. G.A. Thurston, Continuation of Newton’s method through bifurcation points, J. Appl. Mech. 36 (1969) 425430. of Newton’s method to the problem of elastic stability, J. Appl. Mech. 39 (1972) [31 E. Riks, The application 1060-1066. related to nonlinear theories of solids, Internat. J. Solids Structures 7 141 G.A. Wempner, Discrete approximations (1971) 1581-1599. approach to the solution of snapping and buckling problems, Internat. J. Solids [51E. Riks, An incremental Structures 15 (1979) 329-551.
[ll PI
I. Fried, On nonlinear equilibrium curves
297
[6] M.A. Crisfield, A fast incremental/iterative solution procedure that handles “snap-through”, Comput. & Structures 13 (1981) X-62. [7] M.A. Crisfield, An arc-length method including line searches and accelerations, Internat. J. Numer. Meths. Engrg. 19 (1983) 1269-1289. [S] J. Padorean and S. Tovichakchaikul, Self adaptive predictor-corrector algorithms for static nonlinear structural analysis, Comput. & Structures 15 (1982) 365-377. [9] J. Padorean and S. Tovichakchaikul, Operating characteristics of hyperbolically and elliptically constrained self-adaptive incremental Newton-Raphson algorithms, J. Franklin Institute 316 (1983) 197-223. [lo] E.L. Allgower and K. George, Simplicial and continuation methods for approximating fixed points and solutions to systems of equations, SIAM Rev. 22 (1980) 28-85. [ll] E.L. Allgower, A survey of homotopy methods for smooth mappings, in: E.L. Allgower et al., eds., Numerical Solution of Nonlinear Equations (Springer, Berlin, 1981). [12] L.T. Watson and S.M. Holzer, Quadratic convergence of Crisfield’s method, Comput. & Structures 17 (1983) 69-72. [13] E. Ramm, Strategies for tracing the nonlinear response near limit points, in: W. Wunderlich et al., eds., Non-linear Finite Element Analysis in Structural Mechanics (Springer, Berlin, 1981). [14] I. Fried, Numerical Solution of Differential Equations (Academic Press, New York 1979). [15] I. Fried, Stability and equilibrium of the straight and curved elastica-finite element computation, Comput. Meths. Appl. Mech. Engrg. 28 (1981) 49-61. [16] I. Fried, Nonlinear finite element computation of the equilibrium and stability of the circular plate, Internat. J. Numer. Meths. Engrg. 17 (1981) 1427-1440. [17] I. Fried, Finite element computation of large rubber membrane deformations, Internat. J. Numer. M&c_. Engrg. 18 (1982) 653-660. [18] I. Fried, Nonlinear finite element computation of the equilibrium, stability and motion of the extensional beam and ring, Comput. Meths. Appl. Mech. Engrg. 38 (1983) 29-44. [19] I. Fried, Large deformation static and dynamic finite element analysis of extensible cables, Comput. & Structures 15 (1982) 315-319. [20] T.Y. Yang, Matrix displacement solution Structures 9 (1973) 829-842. [21] C.B. Haselgrove, The solution of non-linear conditions, Comput. J. 4 (1961) 255-259.
to elastica equations
problems
of beams
and of differential
and
equations
frames,
Internat.
with two-point
J. Solids boundary