- Email: [email protected]
Computational bifurcation theory: path-tracing, pinpointing and path-switching
Engineering Structures. Vol.
19. No. 5, pp. 385-392, 1997 (,~) 1997 Elsevier Science Lid Printed in Great Britain. All rights reserved Ill41 0296/97 $17.00 + ().IX)
PII: S0141-0296(96)00094-6
ELSEVIER
Computational bifurcation theory: path-tracing, pinpointing and pathswitching F. Vujii Department of Civil Engineering, GiJh Universio,, Gifu 501-11, Japan
E. Ramm lnstitut fiir Baustatik, Universi~, of Stuttgart, Pfaffenwaldring 7, D-70550 Stuttgart, Germany
The stability points in nonlinear elasticity may be classified into limit points and bifurcation points. For limit points, the path-tracing scheme to successively compute the regular equilibrium points on the equilibrium path and the pinpointing scheme to precisely locate the singular equilibrium points are sufficient for the computational stability analysis. For bifurcation points, however, a specific procedure for path-switching is also necessary to detect the branching paths to be traced in the postbuckling region. The present paper briefly describes these fundamental strategies, path-tracing, pinpointing and path-switching, in computational bifurcation theory. © 1997 Elsevier Science Ltd. All rights reserved. Keywords: stability, bifurcation, nonlinear solution procedures
1.
finite element context, equation (1) may usually be written as
Introduction
The general stability theory developed by Koiter j and Thompson and Hunt 2 is useful in theoretical buckling analysis. The highly mathematical stability theory seems, however, not always feasible in computational practice. More specifically, use of the higher order derivatives of the equilibrium equations, for example, is not realistic in existing finite element codes. Computational bifurcation theory based on recent computational technology such as finite element analysis, finite displacement theory 3, the arc-length method4,5 and the extended system6,7, will approach bifurcation problems from a more practical and computational viewpoint. An excellent review on the bifurcation analysis of geometrically nonlinear structures has been given by Choong and Hangai 8. The present paper describes the fundamental strategies, path-tracing, pinpointing and pathswitching9 ~2, in computational bifurcation theory.
2.
R(u ) - pe = 0
with the internal resistance R(u), loading parameter p and its reference load vector e. The linearized equations of (2) are Kdu -dpe=O
(3)
where K is the tangent stiffness matrix and ~ , 0r Or '
K-, . . .z...,.
j= 1
(4)
)Li
using the eigenvalue Aj and the normalized eigenvector 0r. At the regular equilibrium point A close to a simple bifurcation point B with A~--0, the term including the critical eigenpair (A], 01) is dominant on the right-hand side in equation (4), so that it holds approximately that
Mathematical background
For a nonlinear elastic structure with N degrees of freedom and subject to a single proportional loading, let the equilibrium equations be E(u,p) = 0
(2)
0t0T AI
K~'~- . . . .
( 1)
with nodal displacements u and load parameter p. In the
(5)
and any column d (1 -- d -< N) of K~~ approximates the
385
Computational bifurcation theory: F. Fujii and E. Ramm
386
direction of the critical eigenvector 0~, unless the component d of 0t is equal to zero. (column d OfKA~) ((component d of 0,)] \ hi 1 01
(6)
and (column d of K~ ~) = KA~f
(7)
with f being a unit base vector I
d
When equation (15) is used in Newton-Raphson steps, the observed eigenpair (A,0) should be such that it will be continuously changing and become critical at the singular equilibrium point to be attained. To compute (A,0) at the current point, the Rayleigh quotient iteration with cubic convergent property is recommended. An accidental switch to another eigenpair during Rayleigh quotient iteration may cause unsuccessful pinpointing, or attain a different target which is not aimed at. This, however, can be avoided by using the updated eigenpair as the initial approximation for the Rayleigh quotient iteration during the Newton-Raphson procedure. All these computational steps will be used later in the present study to pinpoint the stability point.
N
f r = (0,...,1,...,0)
(8)
f may also be any arbitrary vector in equation (7) representing a 'one-step' inverse iteration method. This is, however, true only in the vicinity of a simple bifurcation point. For multiple bifurcation, a more reasonable choice is necessary for f, as will be shown later. The displacement increment Au may generally be split into
Au = Auo + 61) AUp + 6q Au.
t')
AKo = el, {K(u + e,,Au,) - K(u)}
11)
AKp = ( 2 p ) { K ( u + et,A u i , ) - K ( u ) }
12) 13)
To predict the eigenvalue increment 6A due to AK, the formula (Chatelin ~3) 6A = 0TAK 0
(14)
will be used to pinpoint the stability point. 0 is the observed eigenvector computed at the current point. From equations (10) and (14), it follows that 6A = 6Ao + 6p6Ap + 6qSAu
(19)
duTdu + dp 2 = A 2
(20)
K6u-6pe=-E
(21)
dUT6U + dp6p = 0
( 22 )
for (6u, 6p). An interesting alternative to the constraint during iteration might be to minimize the magnitude of the displacement residual 6u with respect to 6,o~4. Fried ~5 also proposed the orthogonal trajectory accession to the equilibrium path from an arbitrary nonequilibrium point.
4. AKq = ( 2 u ) { K ( u + e q A u q ) - K ( u ) }
Kdu -dpe=0
with a specified magnitude A of the tangent predictor (du, dp). The corrector (Su, 6/)), for example, being normal to the computed predictor (du, dp), will be calculated by solving
10)
with the finite difference approximations
Tracing the equilibrium path
The arc-length control is the most generally applicable method to trace the equilibrium path e (Figure 2) defined by equation(I). A number of schemes have been established in the latter half of the 1970s and the first half of the 1980s 3 5. Only a fundamental scheme is outlined below due to the limited space. The tangent predictor (du, dp), may be computed by solving
(9)
where (6p, &l) represent the correctors of two loading parameters (p, q), as they will be used later in the globally convergent nonlinear schemes, and (Auo, Au F, Auq) are the incremental displacement vectors for the imbalance load vector ( - E ) and for (6p= 1, ~q= 1), respectively. The change z~K in the tangent stiffness matrix due to Au also may be split into A K = AK,, + 61) AKp + 6q AK, t
3.
Pinpointing the stability point
To precisely compute the stability point, the present paper proposes locally as well as globally convergent nonlinear schemes. In the local scheme, the starting point must be close to the stability point. In the global schemes, however, it may be far from the target point. 4.1. Local scheme (Newton-Raphson procedure) Chan ~6proposed a numerical procedure to accurately determine the singular equilibrium point. For direct computation of the stability points, the 'extended system approach' is available 6'7, in which
(15)
with 6Ao = Or AK,,O
(16)
6Ap = Or AKp 0
(17)
6Aq = 0 T Z]kKqO
(18)
E(u~) = 0
(23)
KO = 0
(24)
101 : l
(25)
will be solved for (2N+ 1) unknowns, (u,O,p), by the Newton-Raphson procedure.
Computational bifurcation theory: F. Fujii and E. Ramm
387
E(u~)
In the present study, an alternative to the extended system is used as follows. E(u,p) = 0
(26)
h = 0
(27)
A - qAA = 0
(28)
6A = - h
(29)
The eigenvalue increment 6A is predicted by (30)
6A = 6A,, + 6p6hp
8u = 8u,, + 8pSul,
(32)
KAOA= hAOA
(38)
- dq E A
K IE
(33)
6 u , = + K Je
(34)
The iterative steps will be repeated, until the solution becomes convergent. The Newton-Raphson iteration is locally convergent and works only when the assumed predictor provides an initial approximation close to the target solution. To obtain a reliable predictor, the following globally convergent schemes are proposed.
(39) (40)
duTdu + d p 2 + dq 2 = A 2
(41)
E A = -(E
-
qEa )
(42)
6)t - 6q~ A = - ( )t - qAa )
(43)
du ~ 8u + dp ~p + dq 6q = 0
(44)
, respectively. A is a step-size to be specified. When q changes its sign on the homotopy path in the (u,p,q)-space, one should switch to the local pinpointing scheme, as described in Section 4.1. The stability point may also be attained by tracing the detour (Figure 1) defined by
Global s c h e m e s 9
The local procedure described above can be combined with an arc-length method to come close to a stability point by tracing the equilibrium path. When it is determined that the current equilibrium point is close to or beyond the stability point, the procedure can simply switch from path-tracing to the local pinpointing scheme. Two more global schemes, homotopy and the 'detour' method (Figure 1 ), are possible. The homotopy path (Figure 1) is defined by
0
dA - d q A A = 0
K 6 u - 61) e - 8q
6u,,=-
=
and
with
4.2.
(37)
Kdu -dpe
A + 8A,, (31)
EA = E ( U A ~ A )
The starting point A with q = 1 and the target point with q = 0 are obviously on the homotopy path, so that by tracing the homotopy path the stability point with A = 0 is attainable. The predictor (du,dp,dq) and corrector (Su,6p,6q) to trace the homotopy path may be computed by
with 6q = 0 in equation (15), and equations (28) and (29) can be solved to obtain the corrector (6u,Sp). 8p = -
(36)
where E a and A.A a r e the imbalance load vector and the observed eigenvalue at the starting point A (nonequilibrium point), respectively.
These equations may be linearized for a Newton-Raphson iteration as g 6 u - 6p e = - E
(35)
-QU A = 0
E(u,p) - (q - q2)f= 0
(45)
A-(1-q2)A
(46)
a=0
to connect the current equilibrium point A with q = 0 and the stability point with q = 1. The v e c t o r f will be identified later in equations(59)-(61). The switch to N e w t o n Raphson iterations must be initiated once the path point with q = 1 is passed. In the predictor and corrector steps, K d u - dp e - dq( l - 2q ) f : 0
(47)
dA - d q ( - 2 q ) A a = 0
(48)
du T du + dp 2 + dq 2= A 2
(49)
and q=l
K 6 u - ~p e - ~q( l - 2q L'~
~ q=0
q=l
= - {E - (q - q2b¢} ~A - 8q(-2q)AA = - {A - (1 duT 6u + dp 6p + dq 8 q = 0
A
q=
Figure 1 Global schemes
(50) - q2)AA}
(51) (52)
must be solved for (du,dp,dq) and (6u,6p,&/), respectively. The path-tracing scheme for the homotopy path and detour is almost the same as the arc-length control to compute the equilibrium path, except that in the global schemes, there are two load parameters (p, q) and the eigenvalue constraint is included in the path equations.
Computational bifurcation theory: F. Fujii and E. Ramm
388
5.
f = 0~
Path.switchingg-]2
for hi = 0 at B
f = cosa0~ + sina02 for hi = A2 = 0 at B
Path-switching is computationally nothing else than a nonlinear solution of the equilibrium equations,
(59) (60)
and
E(U,pA ) = 0
(53)
to find equilibria other than UA, when the load p is fixed at PA (Figure 2). Newton-Raphson iterations are applicable only after a reliable 'branching predictor' is obtained. The present study introduces a globally convergent nonlinear solution procedure to reliably switch from the primary path to the branching path. Modify the equilibrium equations in such a way that
F(u,q) = 0
(54)
F(u,q) = E(u,pA) - q f
(55)
with
in which q and f are a disturbing load parameter and its associated load vector, respectively. Equation (54) defines a solution path ~ in the (u,q)-space and q will be changing along c. When q = 0, ~ will intersect the equilibrium path (Figure 2). Linearizing equation (54) at point A as (56)
g A du A = f ( f o r dq = 1)
provides the initial tangent direction of deformation away from the primary path. Before identifying the vector f in equation (55), recall here that the branching predictor is usually constructed by
f = cosacos/30~ + sinacos/302 + sin/303 lbr h~ = h~
=
h 3 =
0 at B
(61)
with (a,/3) being the longitude and latitude to identify a point on a unit sphere surface in (o-], 0-2, 0-3)-space, as shown in Figure 3. Care must be taken here, that ( 0~, 02, 03) used in equations (59)-(61) should not be computed at B, but at the point A with h~ # 0, A2 # 0 and h 3 ~ 0. Equations (59)-(61) are the best choices of vectorfto push the structure in the most prospective direction to the target C on the branching path, because for example the structure will deform in the direction of
(',)
du,~= h
(')
cosa0p + A2 sina02
(62)
away from A to the target, when the vector f is chosen as in equation (60) for two zero eigenvalues, At = A~ = 0, at B. Let ff and rl be the projection of ~c on the u-space and the direction of the unit tangent vector du A
tA = ia, r
t 63)
du A
duH=dut+ o'~01 + 0"202 + 0"303
(57)
dp = ~ (specified)
(58)
using 'eigenvector injection'. (du~, dp) represents the tangent vector of the primary path at the bifurcation point B and may approximately be replaced by the tangent vector computed at the regular equilibrium point A close to B. (0~, 02, 03) are the critical eigenvectors identified at a threefold bifurcation point B with A l = /~2 = /~3 ~--0. (O-l, Or2, 0"3) are indefinite variables in the random access to the bifurcation paths. From equation (57), it may be assumed that
, respectively (Figures 2 and 4). The paths ~: and ~/ serve as a guideline to the target. The following two practical strategies are recommended for path-switching.
5.1. Line search (Figure 4) Conduct a line search in the "0-direction to detect a stationary point of the potential by assuming that u ( ~ ) = UA + ~Jta
The termination criterion of the line search is
ErtA = 0
(65)
The detected stationary point C' in the line search direction is usually nearby the target C and provides a reliable branching predictor. No equations need to be solved during
0" 3
._-o
q /
Figure2 Path-switching
(64)
Figure3 a and /3
Computational bifurcation theory: F. Fujii and E. Ramm
389
-q
TJ'(~'~I'1)2,P) = (1) 1 + ~2) 2 + 0.5k{sin(~o + ~2) - sin(~o + ~ , )}2 sin2~, c'
c
-p{2cos1)o - cos(l)o + l~t )sin(l~o + ~11) - cos(D.o + ~L-)}
(66)
k = (KL2)/C
(67)
with
Figure 4 Line search line search and only the imbalance vector E must be computed.
5.2. ~-Tracing (Figure 2) The arc-length control is applicable to trace ~: from A to C. q = 0 is the termination criterion of ~-tracing. The control of q does not work due to a turning point of q on ~: between A and C. ~-tracing is more time-consuming compared to line search. However, not only C, but other equilibria far away from point A may be detected. This is a significant advantage of ~:-tracing as a globally convergent nonlinear scheme. 6.
Numerical examples
The computational strategies described above are all tested by numerical examples. Although seemingly simple and small, the computed bifurcation models and a shell example are in reality difficult from the computational point of view. The finite element program 'CARAT' which was developed at the Institut ftir Baustatik, University of Stuttgart, Germany, was used in the shell analysis.
(1)1, ~ 2 , ~0) are two degrees of freedom and an imperfection measure, respectively. (K, C) in equation (67) represent the spring stiffness for elongation and rotation, respectively, and k stands for the relative stiffness. The equilibrium equations and tangent stiffness matrix may easily be derived analytically from equation (66). For the model with (~-)~o,1)k,k) (5°,90°,2), the careful inspection of the number of the negative pivots in the LDLr-decomposed stiffness matrix revealed two bifurcation points, BPI and BP2 (Figure 6), and the local scheme precisely computed their locations as BP1 : (~,,122,p) = (+0.4296,+0.4296,+3.477) and BP2: (~l,~2,p) = (+ 1.483,+ 1.483,+5.934 ). The eigenvector 0T=(+0.7071,--0.7071) is critical at both bifurcation points. The Newton-Raphson procedure, when initiated at (['~,,['~2,P) = (0,0,0), pinpointed BP2, at which (1),, + 1~,1)o + 1)2) = (w/2,¢r/2), as Figure 6 shows. The detour method also worked well starting from (flt,[~2,P) = (0,0,0) and successfully pinpointed BPI, when the corresponding eigenpair was observed. The homotopy path, starting from a nonequilibrium point, (1)1,~2,P) -(+0.1,-0.1,+2.0), could also guide the nonlinear solution to BP1. =
P
|~z2 . . . . "'-.
°/t
6. I. Bench mark for pinpointing The simple bench model in Figure 5 was designed to test the pinpointing strategies. The potential for the twodegrees-of-freedom model is written in dimensionless form as 1)o= f , l~k = 90~, k-- 2"0
/0
~
i l
u,'
!
!
C~
I
Figure 5 Bench mark model
!
0~=(-0.7071 ,-0.7071 )
!
!
/
lie+ I'll
,~o + IQ2
01=(+0.7071 ,-0.7071 ) I
I
I
I
I
I
"","~
2 i
a,ioo
~" _ . ~ j
!
Deformed configuration
D.L
Figure6 BP1 and BP2
I
I
I
I
Computational bifurcation theory: F. Fujfi and E. Ramm
390
6.2. Multiple bifurcation model The inverted-L frame with two elastic hinge connections in Figure 7 was designed to validate the path-switching scheme in multiple bifurcation. For the initial postbuckling region, in which two degrees of freedom (~,f~2) are sufficiently small, the potential may be remarkably simplified in the following dimensionless form.
a P=4.2 -0.4 -0.6 - 0 . ~ ~
-0.2
0.4
(
0
0.8 0 8 6 ~2 0 -0.2 -0.4 -0.6 -0.8
7r(~Ql,~Q2,p ) = +2(fl2 + f12) _ p
(112 + 1)2) + 4 2 - a 2 + 4 2 - a ~ J
(68)
This simple potential expression enables an analytical treatment of the model. The potential map is plotted in Figure 8 at the load level higher (P = 4.2) and lower (P = 3.7) than the buckling load (Per = 4.0). Beyond the multibifurcation point, there is no equilibrium point on the potential surface at the load level with P--4.2, except on the primary path, while at the lower load level with P = 3.7, eight distinct equilibrium points exist on four bifurcation paths. All the branching paths are symmetric and unstable. The equilibrium points in one-quarter of the symmetric potential map could be detected by line search and ~-tracing (Figure 9). The two critical eigenvectors were used in the linear combination with one indefinite parameter cK0 --< c~ ~ ~). For an increment z~a = ~, vector f was assumed in seven different directions in one-quarter surface and a stationary point was detected in each line search direction. By Newton-Raphson iterations, the target point nearby could be easily attained from the stationary point detected. In the ~-tracing also, v e c t o r f need not be initially in the direction to the target. The tangent direction of the ~:curve has been corrected stepwise during E-tracing, finally hitting the target, Line search and ~-tracing have been compared for o~= ~r/6 and ~/3. For the same vector f, different bifurcation paths were attained. This is due to the fact that the line search always follows a straight line, while during ~-tracing, the path ~ to the target will be such that along the steepest descent E of the potential is parallel to f due to E = q f, as shown by equations (54) and (55).
6.3. Toggle frame The third example is a two-dimensional clamped-clamped toggle frame j7.~8, as shown in Figure 10. The system paraL 2
I.
:
L 2
P
k 2///"'; Figure 7 Inverted L-frame
-
I
1
1) P=3.7
Per=4.0
-0.4 -0.2
0.2 0.4
e
fle 0 0.2
Figure 8 Potential map
meters to identify the geometry and stiffness of the frame are as follows. L (half span) = 328.56 cm, H (height) = 38.6 cm, E (Young's modulus)= 100 kN/cm 2, A (cross-sectional area) = 84.16cm 2 and 1 (cross-sectional inertia) = 265.74 cm 4 A total of 10 straight beam elements are used in discretization (27 degrees of freedom). A vertical load P is assumed at the midpoint. The characteristic feature in this example is a number of stability points on the primary path (Figure lO), which exhibits snap-through and bifurcation. Two limit points and six bifurcation points were inspected during tracing of the primary path and all these stability points could be precisely located by the local pinpointing scheme. The eigenvectors, which are critical at the first four singular points, are illustrated in Figure 10. The orthogonality condition (0re = 0) of the critical eigenvector at the bifurcation points and the reference load vector can easily be examined. An accidental switch to another eigenpair during iterations might cause unsuccessful pinpointing or we might attain a different stability point which is not expected. This is, however, unlikely to occur during iteration in the vicinity of the target point. During path-switching, there was no particular problem in attaining the bifurcation paths. Path-tracing, pinpointing and path-switching are three indispensable strategies in computational bifurcation theory.
Computational bifurcation theory: F. Fujii and E. Ramm
391 P
,,,/,,/,/,,/
/
o.5_..,,~% o s i,)l 0
.
~
0.4
P=3.7
/ o,
/
o
f
/
f
a
/
/
~p
R=83.33m a= 10m h=0.1m E = 34 000 000 kN/m 2 v =0.2 p = h x 1512 kN/m
Figure 11 Shell panel
P
f=cos ~01 + sin a 02, Aa=15 ~
2 .00E+00' 1.90 1.80 1.70 1.60
1.50
0
.
~
~
displacement
1.40 1.30 1.20 1.10 1.00 0.90
0.5
0.80t o.7ot 0.60t 0.50t
(/.401(I.3010.20 t
W
O. I 0~"
-4.81`)
-4.00
-3.20
-2.40
-1.61,)
0.81,)
0.80E-01
Figure 12 Path-switching in asymmetric bifurcation
~-tracing Figure 9 Line search and ~C-tracing '
I
'
I
'
LPf'~BP3
/ ~.
BP1
\
~
I
~
__
.
~
LPI
LP21BP5
I 6O
Eigenvectors BP4 ~
I 0
I 20
I
A
I 40
Figure 10 Toggle frame
6.4. Shell bifurcation The cylindrical shell panel under axial load (Figure 11) was computed to test the path-switching strategy in shell bifurcation ~°. The membrane displacements are not constrained along the four hinged boundary lines. One quarter of the shell panel was discretized by 2 x 2 elements. In Figure 12, the load parameter, scaled by the buckling load, is plotted against the normal deflection W(=A) at the panel centre. Except for the negligible Poisson effect, the primary path coincides with the load axis. In the vicinity of the asymmetric bifurcation point, both line search and ~-tracing
worked well for path-switching. When the panel buckled inwards, the secondary path was unstable with one negative diagonal element in the LDLr-decomposed stiffness matrix. Beyond the lower limit point, no negative diagonal element was monitored and the path became stable again. ~-tracing was then successfully initiated to attain another equilibrium point on the secondary branch, which was associated with the outward buckling of the panel. After computing several equilibrium points on the attained equilibrium branch, an attempt to jump back to the previous secondary branch was also successful. During ~-tracing in these two attempts, the unstable primary branch beyond the asymmetric bifurcation point could be crossed. In this example, it is clear that line search works well when the target equilibrium point is close to the starting point. For a long-distance path-switching, ~tracing with a globally convergent property is more reliable.
7.
Conclusions
The computational bifurcation theory was briefly outlined along with its fundamental strategies. So far as the simple branching is concerned, there seem to be no problems in both theory and computation. Therefore, research efforts in the future may be directed towards multiple-bifurcation and hill-top branching. The major task is to find all multi-equilibria in the initial postbuckling region. The group theoretic bifurcation analysis might be a promising research direction in the investigation of multiple-bifurcation due to symmetry. Hill-top branching at a limit point may frequently occur but remains unexplored in the computational bifurcation theory.
392
Computational bifurcation theory: F. Fujii and E. Ramm
It may be worth noting here the applicability of the proposed computational procedures to large structures with thousands of degrees of freedom. All the global schemes for path-tracing, pinpointing and path-switching may be considered as a procedure to track the equilibrium path, homotopy path, detour and the trajectory ~ by simply repeating a local procedure. There seems to be no particular problem in successively solving the linearized equilibrium equations even for large structures. The computational bifurcation theory, however, requires the eigenvalues and eigenvectors of the stiffness matrix for pinpointing and path-switching. Especially, the eigenpair must be computed precisely in the local pinpointing procedure. An eigenanalysis of large structures might cause numerical difficulties in the procedures. Numerical tests for this issue are necessary in the future.
5
6 7
8
9
10 II
Acknowledgments The present work is based upon the research activity of the first author from September 1991 to August 1992 at the Institut ftir Baustatik, University of Stuttgart, Germany. The research scholarship provided by the Alexander-yon-Humboldt Stiftung is gratefully acknowledged. References 1 Koiter, W. T. 'On the stability of elastic equilibrium', Thesis, Polytechnic Institute, Delft H. T., Paris, Amsterdam, 1945 2 Thompson, J. M. T. and Hunt. G. W. A general theory of elastic stability John Wiley & Sons, Chichester, UK, 1973 3 Yang, Y.-B. and Kuo, S.-R. Theory and analysis of nonlinear J?amed structures Prentice Hall, Englewood Cliffs, NJ, 1994 4 Crisfield, M. A. and Shi, J. 'A review of solution procedures and path-following techniques in relation to the non-linear finite element analysis of structures', in Nonlinear computational mechanics: state of the art, edited by Wriggers, P. and Wagner, W. (Eds), SpringerVerlag, Berlin, pp. 47-68
12
13 14
15 16
17
18 19
Ramm, E. 'Strategies for tracing the non-linear response near limit points', in Nonlinear.finite element analysis in structural mechanic,~. Wunderlich, W. et al. (Eds), Springer-Verlag. Berlin, 1981, pp. 63-89 Wagner, W. and Wriggers, P. "A simple method for the calculation of postcritical branches' Engng Comput. 1988, 5, 103-109 Wriggers. P. and Simo, J. C. "A general procedure l~r the direct computation of turning and bifurcation points', Int. J. Numer. Meth. Engng, 1990, 30, 155-176 Choong, K. K. and Hangai, Y. 'Review on methods of bilurcation analysis for geometrically nonlinear structures', Bull. Int. Assoc. Shells Spatial Struct. (issue dedicated to SE1KEN-IASS Symposiunl on Nonlinear Analysis and Design for Shell and Spatial Structures) 1993, 34 (2), 133-149 Fujii, F. 'Bifurcation and path-switching in nonlinear elasticity', 3rd World Congress on Computational Mechanics, I 5 August 1994, Makuhari, Japan, Extended abstracts, Vol ll, pp 1246-1247 Fujii, F. and Choong, K. K. 'Branch-switching in bifurcation of structures', J. Engng Meeh., ASCE 1992, 118 (8), 1578-1595 Fujii, F., Choong, K. K. and Kitagawa, T. 'Branch-switching in muhibifurcation of structures', in Computing methods in applied .~cienee,~ and engineering, R. Glowinski (Ed), Nova Science Publishers, Inc. NY, pp. 457-466 Fujii, F. and Asada, K. 'Branch-switching in simple spatial bifur cation models', SEIKEN-1ASS Svmp. on Nonlinear Analysis and Design j b r Shells and Spatial Structures, Tokyo, 1993, pp. 515-522 Chatelin, F. Valeurs Propres de Matrices, Masson, Paris, 1988 Chan, S. L. "Geometric and material nonlinear analysis of beam-columns and frames using the minimum residual displacement method', Int. J. Numer. Meth. Engng 1988, 26, 2657 2669 Fried, 1. "Orthogonal trajectory accession to the nonlinear equilibrium curve', Comput. Meth. Appl. Mech. Engng 1984, 47, 283 297 Chan, S. L. 'A non-linear numerical method for accurate determination of limit and bifurcation points', Int. J. Numer. Meth. Enqn~ 1993, 36, 2779-2790 Eriksson, A. 'On accurate descriptions for primary and secondary paths in equilibrium problems', Comput. Struet. 1992. 44 (1/2). 229-242 Fujii, F., Okazawa, S. and Gong, S.-X. 'Computational stability theory - its strategy', ICES-95, 30 July 3 August 1995, Hawaii Eckstein, U. 'Nichtlineare Stabilitfitsberechnung elastischer Schalentragwerke', Mitteilungen Nr. 83-3, Institut ffir Konstruktiven lngc nieurbau, Ruhr-Universitfit Bochum, November 1983
19. No. 5, pp. 385-392, 1997 (,~) 1997 Elsevier Science Lid Printed in Great Britain. All rights reserved Ill41 0296/97 $17.00 + ().IX)
PII: S0141-0296(96)00094-6
ELSEVIER
Computational bifurcation theory: path-tracing, pinpointing and pathswitching F. Vujii Department of Civil Engineering, GiJh Universio,, Gifu 501-11, Japan
E. Ramm lnstitut fiir Baustatik, Universi~, of Stuttgart, Pfaffenwaldring 7, D-70550 Stuttgart, Germany
The stability points in nonlinear elasticity may be classified into limit points and bifurcation points. For limit points, the path-tracing scheme to successively compute the regular equilibrium points on the equilibrium path and the pinpointing scheme to precisely locate the singular equilibrium points are sufficient for the computational stability analysis. For bifurcation points, however, a specific procedure for path-switching is also necessary to detect the branching paths to be traced in the postbuckling region. The present paper briefly describes these fundamental strategies, path-tracing, pinpointing and path-switching, in computational bifurcation theory. © 1997 Elsevier Science Ltd. All rights reserved. Keywords: stability, bifurcation, nonlinear solution procedures
1.
finite element context, equation (1) may usually be written as
Introduction
The general stability theory developed by Koiter j and Thompson and Hunt 2 is useful in theoretical buckling analysis. The highly mathematical stability theory seems, however, not always feasible in computational practice. More specifically, use of the higher order derivatives of the equilibrium equations, for example, is not realistic in existing finite element codes. Computational bifurcation theory based on recent computational technology such as finite element analysis, finite displacement theory 3, the arc-length method4,5 and the extended system6,7, will approach bifurcation problems from a more practical and computational viewpoint. An excellent review on the bifurcation analysis of geometrically nonlinear structures has been given by Choong and Hangai 8. The present paper describes the fundamental strategies, path-tracing, pinpointing and pathswitching9 ~2, in computational bifurcation theory.
2.
R(u ) - pe = 0
with the internal resistance R(u), loading parameter p and its reference load vector e. The linearized equations of (2) are Kdu -dpe=O
(3)
where K is the tangent stiffness matrix and ~ , 0r Or '
K-, . . .z...,.
j= 1
(4)
)Li
using the eigenvalue Aj and the normalized eigenvector 0r. At the regular equilibrium point A close to a simple bifurcation point B with A~--0, the term including the critical eigenpair (A], 01) is dominant on the right-hand side in equation (4), so that it holds approximately that
Mathematical background
For a nonlinear elastic structure with N degrees of freedom and subject to a single proportional loading, let the equilibrium equations be E(u,p) = 0
(2)
0t0T AI
K~'~- . . . .
( 1)
with nodal displacements u and load parameter p. In the
(5)
and any column d (1 -- d -< N) of K~~ approximates the
385
Computational bifurcation theory: F. Fujii and E. Ramm
386
direction of the critical eigenvector 0~, unless the component d of 0t is equal to zero. (column d OfKA~) ((component d of 0,)] \ hi 1 01
(6)
and (column d of K~ ~) = KA~f
(7)
with f being a unit base vector I
d
When equation (15) is used in Newton-Raphson steps, the observed eigenpair (A,0) should be such that it will be continuously changing and become critical at the singular equilibrium point to be attained. To compute (A,0) at the current point, the Rayleigh quotient iteration with cubic convergent property is recommended. An accidental switch to another eigenpair during Rayleigh quotient iteration may cause unsuccessful pinpointing, or attain a different target which is not aimed at. This, however, can be avoided by using the updated eigenpair as the initial approximation for the Rayleigh quotient iteration during the Newton-Raphson procedure. All these computational steps will be used later in the present study to pinpoint the stability point.
N
f r = (0,...,1,...,0)
(8)
f may also be any arbitrary vector in equation (7) representing a 'one-step' inverse iteration method. This is, however, true only in the vicinity of a simple bifurcation point. For multiple bifurcation, a more reasonable choice is necessary for f, as will be shown later. The displacement increment Au may generally be split into
Au = Auo + 61) AUp + 6q Au.
t')
AKo = el, {K(u + e,,Au,) - K(u)}
11)
AKp = ( 2 p ) { K ( u + et,A u i , ) - K ( u ) }
12) 13)
To predict the eigenvalue increment 6A due to AK, the formula (Chatelin ~3) 6A = 0TAK 0
(14)
will be used to pinpoint the stability point. 0 is the observed eigenvector computed at the current point. From equations (10) and (14), it follows that 6A = 6Ao + 6p6Ap + 6qSAu
(19)
duTdu + dp 2 = A 2
(20)
K6u-6pe=-E
(21)
dUT6U + dp6p = 0
( 22 )
for (6u, 6p). An interesting alternative to the constraint during iteration might be to minimize the magnitude of the displacement residual 6u with respect to 6,o~4. Fried ~5 also proposed the orthogonal trajectory accession to the equilibrium path from an arbitrary nonequilibrium point.
4. AKq = ( 2 u ) { K ( u + e q A u q ) - K ( u ) }
Kdu -dpe=0
with a specified magnitude A of the tangent predictor (du, dp). The corrector (Su, 6/)), for example, being normal to the computed predictor (du, dp), will be calculated by solving
10)
with the finite difference approximations
Tracing the equilibrium path
The arc-length control is the most generally applicable method to trace the equilibrium path e (Figure 2) defined by equation(I). A number of schemes have been established in the latter half of the 1970s and the first half of the 1980s 3 5. Only a fundamental scheme is outlined below due to the limited space. The tangent predictor (du, dp), may be computed by solving
(9)
where (6p, &l) represent the correctors of two loading parameters (p, q), as they will be used later in the globally convergent nonlinear schemes, and (Auo, Au F, Auq) are the incremental displacement vectors for the imbalance load vector ( - E ) and for (6p= 1, ~q= 1), respectively. The change z~K in the tangent stiffness matrix due to Au also may be split into A K = AK,, + 61) AKp + 6q AK, t
3.
Pinpointing the stability point
To precisely compute the stability point, the present paper proposes locally as well as globally convergent nonlinear schemes. In the local scheme, the starting point must be close to the stability point. In the global schemes, however, it may be far from the target point. 4.1. Local scheme (Newton-Raphson procedure) Chan ~6proposed a numerical procedure to accurately determine the singular equilibrium point. For direct computation of the stability points, the 'extended system approach' is available 6'7, in which
(15)
with 6Ao = Or AK,,O
(16)
6Ap = Or AKp 0
(17)
6Aq = 0 T Z]kKqO
(18)
E(u~) = 0
(23)
KO = 0
(24)
101 : l
(25)
will be solved for (2N+ 1) unknowns, (u,O,p), by the Newton-Raphson procedure.
Computational bifurcation theory: F. Fujii and E. Ramm
387
E(u~)
In the present study, an alternative to the extended system is used as follows. E(u,p) = 0
(26)
h = 0
(27)
A - qAA = 0
(28)
6A = - h
(29)
The eigenvalue increment 6A is predicted by (30)
6A = 6A,, + 6p6hp
8u = 8u,, + 8pSul,
(32)
KAOA= hAOA
(38)
- dq E A
K IE
(33)
6 u , = + K Je
(34)
The iterative steps will be repeated, until the solution becomes convergent. The Newton-Raphson iteration is locally convergent and works only when the assumed predictor provides an initial approximation close to the target solution. To obtain a reliable predictor, the following globally convergent schemes are proposed.
(39) (40)
duTdu + d p 2 + dq 2 = A 2
(41)
E A = -(E
-
qEa )
(42)
6)t - 6q~ A = - ( )t - qAa )
(43)
du ~ 8u + dp ~p + dq 6q = 0
(44)
, respectively. A is a step-size to be specified. When q changes its sign on the homotopy path in the (u,p,q)-space, one should switch to the local pinpointing scheme, as described in Section 4.1. The stability point may also be attained by tracing the detour (Figure 1) defined by
Global s c h e m e s 9
The local procedure described above can be combined with an arc-length method to come close to a stability point by tracing the equilibrium path. When it is determined that the current equilibrium point is close to or beyond the stability point, the procedure can simply switch from path-tracing to the local pinpointing scheme. Two more global schemes, homotopy and the 'detour' method (Figure 1 ), are possible. The homotopy path (Figure 1) is defined by
0
dA - d q A A = 0
K 6 u - 61) e - 8q
6u,,=-
=
and
with
4.2.
(37)
Kdu -dpe
A + 8A,, (31)
EA = E ( U A ~ A )
The starting point A with q = 1 and the target point with q = 0 are obviously on the homotopy path, so that by tracing the homotopy path the stability point with A = 0 is attainable. The predictor (du,dp,dq) and corrector (Su,6p,6q) to trace the homotopy path may be computed by
with 6q = 0 in equation (15), and equations (28) and (29) can be solved to obtain the corrector (6u,Sp). 8p = -
(36)
where E a and A.A a r e the imbalance load vector and the observed eigenvalue at the starting point A (nonequilibrium point), respectively.
These equations may be linearized for a Newton-Raphson iteration as g 6 u - 6p e = - E
(35)
-QU A = 0
E(u,p) - (q - q2)f= 0
(45)
A-(1-q2)A
(46)
a=0
to connect the current equilibrium point A with q = 0 and the stability point with q = 1. The v e c t o r f will be identified later in equations(59)-(61). The switch to N e w t o n Raphson iterations must be initiated once the path point with q = 1 is passed. In the predictor and corrector steps, K d u - dp e - dq( l - 2q ) f : 0
(47)
dA - d q ( - 2 q ) A a = 0
(48)
du T du + dp 2 + dq 2= A 2
(49)
and q=l
K 6 u - ~p e - ~q( l - 2q L'~
~ q=0
q=l
= - {E - (q - q2b¢} ~A - 8q(-2q)AA = - {A - (1 duT 6u + dp 6p + dq 8 q = 0
A
q=
Figure 1 Global schemes
(50) - q2)AA}
(51) (52)
must be solved for (du,dp,dq) and (6u,6p,&/), respectively. The path-tracing scheme for the homotopy path and detour is almost the same as the arc-length control to compute the equilibrium path, except that in the global schemes, there are two load parameters (p, q) and the eigenvalue constraint is included in the path equations.
Computational bifurcation theory: F. Fujii and E. Ramm
388
5.
f = 0~
Path.switchingg-]2
for hi = 0 at B
f = cosa0~ + sina02 for hi = A2 = 0 at B
Path-switching is computationally nothing else than a nonlinear solution of the equilibrium equations,
(59) (60)
and
E(U,pA ) = 0
(53)
to find equilibria other than UA, when the load p is fixed at PA (Figure 2). Newton-Raphson iterations are applicable only after a reliable 'branching predictor' is obtained. The present study introduces a globally convergent nonlinear solution procedure to reliably switch from the primary path to the branching path. Modify the equilibrium equations in such a way that
F(u,q) = 0
(54)
F(u,q) = E(u,pA) - q f
(55)
with
in which q and f are a disturbing load parameter and its associated load vector, respectively. Equation (54) defines a solution path ~ in the (u,q)-space and q will be changing along c. When q = 0, ~ will intersect the equilibrium path (Figure 2). Linearizing equation (54) at point A as (56)
g A du A = f ( f o r dq = 1)
provides the initial tangent direction of deformation away from the primary path. Before identifying the vector f in equation (55), recall here that the branching predictor is usually constructed by
f = cosacos/30~ + sinacos/302 + sin/303 lbr h~ = h~
=
h 3 =
0 at B
(61)
with (a,/3) being the longitude and latitude to identify a point on a unit sphere surface in (o-], 0-2, 0-3)-space, as shown in Figure 3. Care must be taken here, that ( 0~, 02, 03) used in equations (59)-(61) should not be computed at B, but at the point A with h~ # 0, A2 # 0 and h 3 ~ 0. Equations (59)-(61) are the best choices of vectorfto push the structure in the most prospective direction to the target C on the branching path, because for example the structure will deform in the direction of
(',)
du,~= h
(')
cosa0p + A2 sina02
(62)
away from A to the target, when the vector f is chosen as in equation (60) for two zero eigenvalues, At = A~ = 0, at B. Let ff and rl be the projection of ~c on the u-space and the direction of the unit tangent vector du A
tA = ia, r
t 63)
du A
duH=dut+ o'~01 + 0"202 + 0"303
(57)
dp = ~ (specified)
(58)
using 'eigenvector injection'. (du~, dp) represents the tangent vector of the primary path at the bifurcation point B and may approximately be replaced by the tangent vector computed at the regular equilibrium point A close to B. (0~, 02, 03) are the critical eigenvectors identified at a threefold bifurcation point B with A l = /~2 = /~3 ~--0. (O-l, Or2, 0"3) are indefinite variables in the random access to the bifurcation paths. From equation (57), it may be assumed that
, respectively (Figures 2 and 4). The paths ~: and ~/ serve as a guideline to the target. The following two practical strategies are recommended for path-switching.
5.1. Line search (Figure 4) Conduct a line search in the "0-direction to detect a stationary point of the potential by assuming that u ( ~ ) = UA + ~Jta
The termination criterion of the line search is
ErtA = 0
(65)
The detected stationary point C' in the line search direction is usually nearby the target C and provides a reliable branching predictor. No equations need to be solved during
0" 3
._-o
q /
Figure2 Path-switching
(64)
Figure3 a and /3
Computational bifurcation theory: F. Fujii and E. Ramm
389
-q
TJ'(~'~I'1)2,P) = (1) 1 + ~2) 2 + 0.5k{sin(~o + ~2) - sin(~o + ~ , )}2 sin2~, c'
c
-p{2cos1)o - cos(l)o + l~t )sin(l~o + ~11) - cos(D.o + ~L-)}
(66)
k = (KL2)/C
(67)
with
Figure 4 Line search line search and only the imbalance vector E must be computed.
5.2. ~-Tracing (Figure 2) The arc-length control is applicable to trace ~: from A to C. q = 0 is the termination criterion of ~-tracing. The control of q does not work due to a turning point of q on ~: between A and C. ~-tracing is more time-consuming compared to line search. However, not only C, but other equilibria far away from point A may be detected. This is a significant advantage of ~:-tracing as a globally convergent nonlinear scheme. 6.
Numerical examples
The computational strategies described above are all tested by numerical examples. Although seemingly simple and small, the computed bifurcation models and a shell example are in reality difficult from the computational point of view. The finite element program 'CARAT' which was developed at the Institut ftir Baustatik, University of Stuttgart, Germany, was used in the shell analysis.
(1)1, ~ 2 , ~0) are two degrees of freedom and an imperfection measure, respectively. (K, C) in equation (67) represent the spring stiffness for elongation and rotation, respectively, and k stands for the relative stiffness. The equilibrium equations and tangent stiffness matrix may easily be derived analytically from equation (66). For the model with (~-)~o,1)k,k) (5°,90°,2), the careful inspection of the number of the negative pivots in the LDLr-decomposed stiffness matrix revealed two bifurcation points, BPI and BP2 (Figure 6), and the local scheme precisely computed their locations as BP1 : (~,,122,p) = (+0.4296,+0.4296,+3.477) and BP2: (~l,~2,p) = (+ 1.483,+ 1.483,+5.934 ). The eigenvector 0T=(+0.7071,--0.7071) is critical at both bifurcation points. The Newton-Raphson procedure, when initiated at (['~,,['~2,P) = (0,0,0), pinpointed BP2, at which (1),, + 1~,1)o + 1)2) = (w/2,¢r/2), as Figure 6 shows. The detour method also worked well starting from (flt,[~2,P) = (0,0,0) and successfully pinpointed BPI, when the corresponding eigenpair was observed. The homotopy path, starting from a nonequilibrium point, (1)1,~2,P) -(+0.1,-0.1,+2.0), could also guide the nonlinear solution to BP1. =
P
|~z2 . . . . "'-.
°/t
6. I. Bench mark for pinpointing The simple bench model in Figure 5 was designed to test the pinpointing strategies. The potential for the twodegrees-of-freedom model is written in dimensionless form as 1)o= f , l~k = 90~, k-- 2"0
/0
~
i l
u,'
!
!
C~
I
Figure 5 Bench mark model
!
0~=(-0.7071 ,-0.7071 )
!
!
/
lie+ I'll
,~o + IQ2
01=(+0.7071 ,-0.7071 ) I
I
I
I
I
I
"","~
2 i
a,ioo
~" _ . ~ j
!
Deformed configuration
D.L
Figure6 BP1 and BP2
I
I
I
I
Computational bifurcation theory: F. Fujfi and E. Ramm
390
6.2. Multiple bifurcation model The inverted-L frame with two elastic hinge connections in Figure 7 was designed to validate the path-switching scheme in multiple bifurcation. For the initial postbuckling region, in which two degrees of freedom (~,f~2) are sufficiently small, the potential may be remarkably simplified in the following dimensionless form.
a P=4.2 -0.4 -0.6 - 0 . ~ ~
-0.2
0.4
(
0
0.8 0 8 6 ~2 0 -0.2 -0.4 -0.6 -0.8
7r(~Ql,~Q2,p ) = +2(fl2 + f12) _ p
(112 + 1)2) + 4 2 - a 2 + 4 2 - a ~ J
(68)
This simple potential expression enables an analytical treatment of the model. The potential map is plotted in Figure 8 at the load level higher (P = 4.2) and lower (P = 3.7) than the buckling load (Per = 4.0). Beyond the multibifurcation point, there is no equilibrium point on the potential surface at the load level with P--4.2, except on the primary path, while at the lower load level with P = 3.7, eight distinct equilibrium points exist on four bifurcation paths. All the branching paths are symmetric and unstable. The equilibrium points in one-quarter of the symmetric potential map could be detected by line search and ~-tracing (Figure 9). The two critical eigenvectors were used in the linear combination with one indefinite parameter cK0 --< c~ ~ ~). For an increment z~a = ~, vector f was assumed in seven different directions in one-quarter surface and a stationary point was detected in each line search direction. By Newton-Raphson iterations, the target point nearby could be easily attained from the stationary point detected. In the ~-tracing also, v e c t o r f need not be initially in the direction to the target. The tangent direction of the ~:curve has been corrected stepwise during E-tracing, finally hitting the target, Line search and ~-tracing have been compared for o~= ~r/6 and ~/3. For the same vector f, different bifurcation paths were attained. This is due to the fact that the line search always follows a straight line, while during ~-tracing, the path ~ to the target will be such that along the steepest descent E of the potential is parallel to f due to E = q f, as shown by equations (54) and (55).
6.3. Toggle frame The third example is a two-dimensional clamped-clamped toggle frame j7.~8, as shown in Figure 10. The system paraL 2
I.
:
L 2
P
k 2///"'; Figure 7 Inverted L-frame
-
I
1
1) P=3.7
Per=4.0
-0.4 -0.2
0.2 0.4
e
fle 0 0.2
Figure 8 Potential map
meters to identify the geometry and stiffness of the frame are as follows. L (half span) = 328.56 cm, H (height) = 38.6 cm, E (Young's modulus)= 100 kN/cm 2, A (cross-sectional area) = 84.16cm 2 and 1 (cross-sectional inertia) = 265.74 cm 4 A total of 10 straight beam elements are used in discretization (27 degrees of freedom). A vertical load P is assumed at the midpoint. The characteristic feature in this example is a number of stability points on the primary path (Figure lO), which exhibits snap-through and bifurcation. Two limit points and six bifurcation points were inspected during tracing of the primary path and all these stability points could be precisely located by the local pinpointing scheme. The eigenvectors, which are critical at the first four singular points, are illustrated in Figure 10. The orthogonality condition (0re = 0) of the critical eigenvector at the bifurcation points and the reference load vector can easily be examined. An accidental switch to another eigenpair during iterations might cause unsuccessful pinpointing or we might attain a different stability point which is not expected. This is, however, unlikely to occur during iteration in the vicinity of the target point. During path-switching, there was no particular problem in attaining the bifurcation paths. Path-tracing, pinpointing and path-switching are three indispensable strategies in computational bifurcation theory.
Computational bifurcation theory: F. Fujii and E. Ramm
391 P
,,,/,,/,/,,/
/
o.5_..,,~% o s i,)l 0
.
~
0.4
P=3.7
/ o,
/
o
f
/
f
a
/
/
~p
R=83.33m a= 10m h=0.1m E = 34 000 000 kN/m 2 v =0.2 p = h x 1512 kN/m
Figure 11 Shell panel
P
f=cos ~01 + sin a 02, Aa=15 ~
2 .00E+00' 1.90 1.80 1.70 1.60
1.50
0
.
~
~
displacement
1.40 1.30 1.20 1.10 1.00 0.90
0.5
0.80t o.7ot 0.60t 0.50t
(/.401(I.3010.20 t
W
O. I 0~"
-4.81`)
-4.00
-3.20
-2.40
-1.61,)
0.81,)
0.80E-01
Figure 12 Path-switching in asymmetric bifurcation
~-tracing Figure 9 Line search and ~C-tracing '
I
'
I
'
LPf'~BP3
/ ~.
BP1
\
~
I
~
__
.
~
LPI
LP21BP5
I 6O
Eigenvectors BP4 ~
I 0
I 20
I
A
I 40
Figure 10 Toggle frame
6.4. Shell bifurcation The cylindrical shell panel under axial load (Figure 11) was computed to test the path-switching strategy in shell bifurcation ~°. The membrane displacements are not constrained along the four hinged boundary lines. One quarter of the shell panel was discretized by 2 x 2 elements. In Figure 12, the load parameter, scaled by the buckling load, is plotted against the normal deflection W(=A) at the panel centre. Except for the negligible Poisson effect, the primary path coincides with the load axis. In the vicinity of the asymmetric bifurcation point, both line search and ~-tracing
worked well for path-switching. When the panel buckled inwards, the secondary path was unstable with one negative diagonal element in the LDLr-decomposed stiffness matrix. Beyond the lower limit point, no negative diagonal element was monitored and the path became stable again. ~-tracing was then successfully initiated to attain another equilibrium point on the secondary branch, which was associated with the outward buckling of the panel. After computing several equilibrium points on the attained equilibrium branch, an attempt to jump back to the previous secondary branch was also successful. During ~-tracing in these two attempts, the unstable primary branch beyond the asymmetric bifurcation point could be crossed. In this example, it is clear that line search works well when the target equilibrium point is close to the starting point. For a long-distance path-switching, ~tracing with a globally convergent property is more reliable.
7.
Conclusions
The computational bifurcation theory was briefly outlined along with its fundamental strategies. So far as the simple branching is concerned, there seem to be no problems in both theory and computation. Therefore, research efforts in the future may be directed towards multiple-bifurcation and hill-top branching. The major task is to find all multi-equilibria in the initial postbuckling region. The group theoretic bifurcation analysis might be a promising research direction in the investigation of multiple-bifurcation due to symmetry. Hill-top branching at a limit point may frequently occur but remains unexplored in the computational bifurcation theory.
392
Computational bifurcation theory: F. Fujii and E. Ramm
It may be worth noting here the applicability of the proposed computational procedures to large structures with thousands of degrees of freedom. All the global schemes for path-tracing, pinpointing and path-switching may be considered as a procedure to track the equilibrium path, homotopy path, detour and the trajectory ~ by simply repeating a local procedure. There seems to be no particular problem in successively solving the linearized equilibrium equations even for large structures. The computational bifurcation theory, however, requires the eigenvalues and eigenvectors of the stiffness matrix for pinpointing and path-switching. Especially, the eigenpair must be computed precisely in the local pinpointing procedure. An eigenanalysis of large structures might cause numerical difficulties in the procedures. Numerical tests for this issue are necessary in the future.
5
6 7
8
9
10 II
Acknowledgments The present work is based upon the research activity of the first author from September 1991 to August 1992 at the Institut ftir Baustatik, University of Stuttgart, Germany. The research scholarship provided by the Alexander-yon-Humboldt Stiftung is gratefully acknowledged. References 1 Koiter, W. T. 'On the stability of elastic equilibrium', Thesis, Polytechnic Institute, Delft H. T., Paris, Amsterdam, 1945 2 Thompson, J. M. T. and Hunt. G. W. A general theory of elastic stability John Wiley & Sons, Chichester, UK, 1973 3 Yang, Y.-B. and Kuo, S.-R. Theory and analysis of nonlinear J?amed structures Prentice Hall, Englewood Cliffs, NJ, 1994 4 Crisfield, M. A. and Shi, J. 'A review of solution procedures and path-following techniques in relation to the non-linear finite element analysis of structures', in Nonlinear computational mechanics: state of the art, edited by Wriggers, P. and Wagner, W. (Eds), SpringerVerlag, Berlin, pp. 47-68
12
13 14
15 16
17
18 19
Ramm, E. 'Strategies for tracing the non-linear response near limit points', in Nonlinear.finite element analysis in structural mechanic,~. Wunderlich, W. et al. (Eds), Springer-Verlag. Berlin, 1981, pp. 63-89 Wagner, W. and Wriggers, P. "A simple method for the calculation of postcritical branches' Engng Comput. 1988, 5, 103-109 Wriggers. P. and Simo, J. C. "A general procedure l~r the direct computation of turning and bifurcation points', Int. J. Numer. Meth. Engng, 1990, 30, 155-176 Choong, K. K. and Hangai, Y. 'Review on methods of bilurcation analysis for geometrically nonlinear structures', Bull. Int. Assoc. Shells Spatial Struct. (issue dedicated to SE1KEN-IASS Symposiunl on Nonlinear Analysis and Design for Shell and Spatial Structures) 1993, 34 (2), 133-149 Fujii, F. 'Bifurcation and path-switching in nonlinear elasticity', 3rd World Congress on Computational Mechanics, I 5 August 1994, Makuhari, Japan, Extended abstracts, Vol ll, pp 1246-1247 Fujii, F. and Choong, K. K. 'Branch-switching in bifurcation of structures', J. Engng Meeh., ASCE 1992, 118 (8), 1578-1595 Fujii, F., Choong, K. K. and Kitagawa, T. 'Branch-switching in muhibifurcation of structures', in Computing methods in applied .~cienee,~ and engineering, R. Glowinski (Ed), Nova Science Publishers, Inc. NY, pp. 457-466 Fujii, F. and Asada, K. 'Branch-switching in simple spatial bifur cation models', SEIKEN-1ASS Svmp. on Nonlinear Analysis and Design j b r Shells and Spatial Structures, Tokyo, 1993, pp. 515-522 Chatelin, F. Valeurs Propres de Matrices, Masson, Paris, 1988 Chan, S. L. "Geometric and material nonlinear analysis of beam-columns and frames using the minimum residual displacement method', Int. J. Numer. Meth. Engng 1988, 26, 2657 2669 Fried, 1. "Orthogonal trajectory accession to the nonlinear equilibrium curve', Comput. Meth. Appl. Mech. Engng 1984, 47, 283 297 Chan, S. L. 'A non-linear numerical method for accurate determination of limit and bifurcation points', Int. J. Numer. Meth. Enqn~ 1993, 36, 2779-2790 Eriksson, A. 'On accurate descriptions for primary and secondary paths in equilibrium problems', Comput. Struet. 1992. 44 (1/2). 229-242 Fujii, F., Okazawa, S. and Gong, S.-X. 'Computational stability theory - its strategy', ICES-95, 30 July 3 August 1995, Hawaii Eckstein, U. 'Nichtlineare Stabilitfitsberechnung elastischer Schalentragwerke', Mitteilungen Nr. 83-3, Institut ffir Konstruktiven lngc nieurbau, Ruhr-Universitfit Bochum, November 1983