Computing critical points and secondary paths in nonlinear structural stability analysis by the finite element method

Compurers & Strucrures Vol.58,No. l.pp. 203-220, 1996

Pergamon

Copyright 0 1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0045.7949/96 $9.50 + 0.W

0045-7949(95)00114-x

COMPUTING CRITICAL POINTS AND SECONDARY PATHS IN NONLINEAR STRUCTURAL STABILITY ANALYSIS BY THE FINITE ELEMENT METHOD J. Shi Department of Aeronautical Engineering, Queen’s University of Belfast. David Keir Building, Stranmillis Road, Belfast BT9 5AG, U.K. (Received 11 May 1994) Abstract-Stability analysis is a crucial aspect of structural design, particularly for beam, plate and shell type structures. The stability of a structure is often governed by singular points in its equilibrium path. These points usually pose serious problems to ordinary finite element solution methods and require special treatment. In this paper, the state-of-the-art of finite element path-following techniques is critically reviewed with special attention to bracketing procedures for the singular points and branch switching algorithms. Some practical aspects involved will be explained in detail and demonstrated through numerical examples at the end

the numerical difficulties associated with critical points and branching techniques, in other words the post-critical response. In this paper, only static stability is studied. Moreover we will confine ourselves to the finite element analysis of elastic structural stability. For a broader treatment, the reader is referred to Refs [1, 1621 and 23-291 and also the references therein. In the following, we first discuss the bracketing of singular points, then we explain different branching procedures.

1. INTRODUCTION

Nowadays structures often work in a nonlinear range in order to meet demands like low structural weight but high strength. Since stability is an integral part of structural design, a nonlinear stability analysis is hence inevitable, which often poses a serious challenge to finite element analysts. The two important approaches of stability analysis are perturbation method and continuation method [l]. Although there exist a few finite element versions of the former method [2-51, for various reasons it is the latter method that is predominant in present finite element analysis and thus is the only concern of this paper. The continuation method is usually built on a Newton type iteration scheme [6-91, a two step predictor-corrector approach. As a result the efficiency of the scheme depends on both predictor and corrector. Stability analysis is characterized by two aspects: (1) path following techniques when away from the critical points, see Rheinboldt [lo] for a mathematical treatment and Clarke and Hancock [111,for a recent comparative study of different finite element solution techniques. More up-to-date reviews can be found in Refs [12-l 51. (2) Near the critical points, the determination of critical points and branching of critical points is a bifurcation point. Critical points are divided into limit and bifurcation points. While a solution method for limit points is quite well developed [16], there is still a lot to be done about the bifurcation point, especially for the case of cluster and multiple bifurcation points [ 171. Except, directly at or very close to the critical points, where singularity and large curvature might cause problems, the arc-length type method [l, 1 l-141 can perform quite satisfactorily. So in this paper we will focus on

2.

DETERMINATION

OF

SINGULAR

POINTS

The defining conditions for the critical points are well established, i.e. the tangent stiffness matrix K, becomes singular [ 1,7, 171. In this paper we will concentrate on the case where the nullity of K, is one. In other words the critical point is either a limit or simple bifurcation point. For the case of nullity greater than one, we have multiple bifurcation, which is quite complicated and will only be discussed briefly here. It is also assumed that the lowest eigenvalues are “reasonably” separated so that we are not dealing with clustered bifurcation. The singularity of K, can be interpreted in many different ways, e.g. zero pivot of the factorised K,, zero determinant, zero eigenvalue and so on, which can then be employed as a criterion or test function (Seydel[30], p. 132) to locate the singular point. The two basic types of bracketing schemes are addressed below. 2.1. The indirect method The bi-section related methods. When moving along an equilibrium path, we check the positive 203

204

J. Shi

definiteness of stiffness matrix K, by looking at its smallest pivot, determinant or the smallest eigenvalue: test functions r = Det(K,) = smallest pivot = lowest eigenvalue of K, = lowest frequency of natural vibrations = current stiffness parameter (limit point only) = . Normally we start with a stable equilibrium with a positive r. If we have passed a singular point the opposite is true, provided that the post-critical primary path is unstable. As a result, to locate a singular point is simply a matter of reducing the uncertain interval between two consecutive stable and unstable equilibrium states. To put this in another way, we want to pinpoint a neutrally stable equilibrium. This is actually closely related to the concept of line searches in the theory of optimization [3 1, 321, which find the minimum/maximum of a function while in our case we want to find the zeros of a function. Both of the minimum/maximum and zero can be simply deemed as special solution points of a function. As a matter of fact if we work on the square of the test function T, then it is the minimum that we are looking for. Built on the basic idea of bracketing in the previous various bracketing schemes can be paragraph, devised. The simplest is the bi-section method that just tries to find a solution at half the value of the previous solution control parameter increment A?, whether it is loading factor, controlled displacement or arc-length. From this solution and the previous ones, we can tell the uncertain interval closed by the two solutions having opposite signs of the test function and consequently where the singular point lies. The above process is repeated until the uncertain interval is smaller than a pre-set value (see Fig. la), which can be decided in a number of ways [33]. In general it is not wise to control test function only, as load factor or displacement may change very little even though test function itself is still quite large. Shi and Crisfield [33] used the geometrical average of test function and load factor (if load control) or arc length (if arc length control) as bracketing tolerance. For the case of multiple bifurcation, more than one pivot or eigenvalue are zero, practically it is

acceptable to make the smallest eigenvalue less than a tolerance. To avoid unreasonable excessive values, they are often normalized by their initial values. However if critical point is far away from initial solution, such scaling may be meaningless as the magnitude of test functions can be drastically different. As a result, Shi and Crisfield [33] advocate the geometrical average of test functions just before and after critical point as the scaling factor. Although a zero pivot, determinant and eigenvalue all mean a singular K,, obviously it is computationally most efficient to use the smallest pivot as the test function, because it is readily available from the factorization of K,, while extra computational work is needed to evaluate the determinant and even more work for the lowest eigenvalue. One should take care for large systems whose Det(K,) may be bigger or smaller (depending on structure, units and element used) than the largest or smallest real number a computer can represent. Another drawback with Det(K,) is that it does not change sign if its nullity is even. To improve convergence rate, the relative magnitude of test functions can be employed as a guide to the position of the singular point [7, 381. This introduces the interpolation or weighted average to better approximate the critical point.

‘1,=‘1,

T, -

T?

See Fig. 1b for the definition of the z, and vs. As extrapolation sometimes leads to ironic results, so interpolation is preferred. In fact if we store the closest pair of positive and negative pivots and their associated control parameter, interpolation between these two solutions is guaranteed. Equation (1) tells us that to expect a decent answer, (q2 - PI,) and (r2 - T,) should be neither too small nor too large, which cannot always be guaranteed. So to safeguard against such cases, some sort of cut-off mechanism is necessary when (n2 - n,) is either greater than 80% or less than 20% of its previous value [33]. Instead of

(b)

(a) Fig. I. The bi-section

+=r2.

related

method

205

Computing critical points and secondary paths one can resort to higher linear interpolation, order interpolations, which however may not be optimum, because more computation is necessary to evaluate the extreme solution of a high order polynomial. When two or more pivots become negative due to the presence of a multiple singular point, one of them has to be selected for interpolation. But it is impossible to keep using the same pivot for the subsequent interpolations in order to be consistent. So bi-section may be preferred. There may be more than one singular point along an equilibrium path and it may be necessary to compute all of them to get a global picture of structural response. The nullity of a singular point is equal to the increase in the number of negative pivots or eigenvalues just after passing it. For example when continuing from point 1 to 5 (see Fig. 2a), one will notice changes in the number of negative pivots three times, i.e. between points 1 and 2, 3 and 4, 3 and 4. So there are three singular points. Also the differences of number of negative pivots are 1, 2 and 1, respectively, since Bl and B3 are simple bifurcation points (nullity = 1), while B2 is a double bifurcation point. In practice it is quite possible to step over more than one singular point in one step, even though they are reasonably separated. In the previous example, it is quite likely that increment two passes both Bl and B2 (Fig. 2b). This has two immediate implications: (1) there is a singular point of nullity three between point 1 and 2. If a bracketing is initiated here it could subsequently be found that there are actually two bifurcation points instead of one. (2) If it is the second singular point that is specified to be calculated, then one gets B3 rather than B2. In general, to expect a satisfactory performance, test functions should ideally be monotonic and have a large non-zero slope when crossing zero. Although this is to a large extent decided by the parameterization [38], Tabadd [39] has shown that it is also linked with element nonlinearity effects. This conforms with the author’s own experience that 7 depends on the particular type of element being used. More theoretical investigation on the properties of

Fig. 2. Bracketing

test functions has been done by Mang [40,41]. Also Tao and Ramm [42] have set up a relation between the characteristic of a full system and that of its counterpart in the reduced bases, thus making bracketing in reduced bases possible. Previous bracketing procedures have been built on a scalar quantity-test function. It is also possible to calculate singular point using K, directly [35, 361, where interpolation is based on load factor and displacement, respectively. To handle both bifurcation and limit point it may be easier to use arc length: Det[K: + Al”+‘/Al”(K~+’ -q)]

= 0

(2)

which can be transformed into an eigenvalue problem K; = Al” +‘/Al”(c - K: + ‘)

(3)

to find the next arclength Al”+ ‘. We can see that the indirect method is not efficient in the sense that we have to iterate at a number of intermediate, otherwise unwanted equilibrium points (q, , q2, q3, . . .), to finally reach an acceptable interval. This issue of efficiency is particularly acute when one just wants to decide the critical boundary. In such a case it is the critical points only that are of interest and thus it is highly preferable to have an algorithm that minimizes effort in calculating these bracketing solutions. Shi and Crisfield [33] relax equilibrium tolerance at the beginning of bracketing, thus fewer iterations are required, but gradually tighten it up as the solution approaches a singular point. Chen [37] also uses nonequilibrium states for bracketing purpose. 2.2. Direct method An alternative way to increase computational efficiency is to calculate these points directly [43]. Accuracy is another reason for developing the direct method [30,44,45, 1011.Critical points are also equilibrium solutions, the only thing special about them is that the test functions are zero at such points as well. Consequently if a solution satisfies equilibrium and its test functions are zero at the same time, then

with more than one singular

point present.

206

J. Shi

it must be a critical solution. This is exadty what Abbott [46,47] did more than 10 years ago. He actually used zero determinant as his constraint, which in turn requires one to know analytically the derivative of G [see eqn (4)]. In order to be able to apply the Newton’s method, Moore and Spence [48] worked on the second constraint with non-zero 4. As far as computational efficiency is concerned, Moore’s approach is more appealing. Because of its partitioned formulation, storage requirement and the amount of computation are greatly reduced despite working on a 2n + 1 system.

Equilibrium:

G(u, 2) = R(u) - Iq = 0 Fig. 3. Unfolding (structural imperfection).

Constraints:

Det(K(u)) = 0 K(u)4 = 0

.

II$II = 1 or

&=

1,

(4)

where G(u, 2) is out of balance force, R(u) internal force, q load pattern, 1 load factor. The first contraint on the eigenmode leads to a quadratic equation for the load factor, which means that we have to select a root. This can be avoided by working with its linearization [44,45]. Another possibility is constraining a specific component of the eigenmode. However, we will see later that this will predetermine the singular point that is to be computed. A finite element version of Moor’s scheme has recently been proposed by Wriggers et al. [44,45]. It claimed some success for both elastic and plastic problems. For cases where the analytical form of the second derivatives are not available, an approximation scheme was given. Generally speaking the direct method needs a good starting vector to achieve convergence. It rarely converges to the critical point from anywhere in the load displacement space. Usually one just follows the equilibrium path sufficiently close to the critical point, then he can switch on the direct method to locate it (denoted as direct 1). A by-product of the above scheme is an eigenvector associated with zero eigenvalue, which can be used for branching purpose. Direct bracketing does not differentiate simple critical point (nullity = 1) and multiple critical point (nullity > 1). This makes it difficult to control which singular point to bracket if there are a number of them. Despite all the good things mentioned about the direct method, it should be made clear that for practical problems its applicatoin is limited [30, p. 1611, because the high accuracy of the direct method may not be needed or achievable. This is especially true when the discretization error exceeds a tight tolerance, say, IL, - /1I or when branching point is sensitive to perturbations.

be low and the solution may lose its accuracy, which surprisingly is not often true. According to Seydel[30, p. 1491, this is due to the fact that iteration vectors may lie in the symmetric and anti-symmetric sub-domains. where the Jacobian is no longer singular [49]. Mathematical studies on Newton’s with singularity can be found in Refs method [50-52, 1001. For a general treatment of ill conditioned problems one may consult Rothwell and Drachman [53]. However if convergence is affected by the singularity of K,, two possible remedies can be made, based on the techniques that have been suggested by Riks and Rankin [54] and Thurston et al. [55]. Riks’s method is related to the bordered scheme which relies on the selection of a nonsingular minor of K,. He developed two choices of the minor from information on the lowest eigenvector and path tangent vector, respectively. Although this method is restricted to the limit point only, it has potential to be extended to the bifurcation point. In a more general way Thurston employed an equivalent transformation T which transforms the original iteration into AX=TY, where T is made of the m lowest eigenvectors of K,

and n -m elementary permutation vectors. Unlike the linear Newton iteration method the higher order terms in the Taylor series expansion are retained in Thurston’s adapted approach. Moreover only those elements of Y associated with the eigenvectors are kept in the higher order terms. If m is chosen to be larger than the nullity of K,, a consistent system of equations

can always be guaranteed.

Both of these

2.3. Near singularity In the vicinity

the K, is nearly the convergence rate could

of a critical point,

singular. Consequently

Fig. 4. A simple one bar spring system

Computing

critical

Fig. 5. Unfolding

points

and secondary

by structurai

two methods involve extra computation at each load step, but this may be justified on the efficiency side as fewer load steps are required to reach the desired load level.

207

paths

or load perturbation

A more straight forward approach has been suggested by Wriggers and Simo [45], which comes from the penalty formulation by Fellipa [99]. With this technique, the smallest

(a) 26.Or

m

Analytical solution FEM solution I

-26.0L

-

.

Analytical solution FEM solution

Fig. 6. (a) Branching by structural perturbation-load u horizontal displacement u. (b) Branching by structural perturbation-load u vertical displacement r.

208

J. Shi

(a) 26.C

I .

Analytical solution FEM solution

I

5500

-26.0

(b)

.

Analytical solution FEM solution

Fig. 7. (a) Branching by load perturbation-load u horizontal diplacement U. (b) Branching by load perturbation-load u vertical displacement v.

leading diagonal of K, is replaced by a large machine precision dependent constant. Since equilibrium is checked directly with the internal and

external forces, a modified K, may only affect the convergence rate without the structure itself being changed.

Fig. 8. A simple asymmetrical bifurcation problem.

Computing critical points and secondary paths

209

Table I. Bracketing results for the two bar truss with asymmetric bifurcation Bi-section

Determinant

Smallest pivot

Direct 1

Direct 2

11

4

4

3

5

Bracketing inc.

3. BRANCH SWITCHING

For various reasons it is sometimes desirable to trace secondary branches in order to get a general view of the entire branching diagram. Roughly speaking, there are again two ways of applying the direct and indirect approaches, the relative merits of which will be evaluated in this section. 3.1. Indirect method-unfolding

It is well known now that only a limit point is generic and a bifurcation point is not [30, p. 781. In other words the latter only exists when structure and loading systems are perfect without any imperfections from manufacturing, assembling, etc. This makes it possible to change a bifurcation problem into that of a limit point one by perturbing the original perfect system, the so-called unfolding [ 11. Provided that the right perturbation is applied, the perturbed solution will be close to the perfect solution. As the limit point is no longer a problem. we can employ any existing standard path following method. To be more precise, what one needs to do first is an analysis on the perfect system along the fundamental path OAB (Fig. 3). On passing from A to B, he or she should notice the presence of a bifurcation point between A and B, e.g. by checking the diagonal pivots and the current stiffness parameter at the same time. Then some initial imperfection is introduced and the perturbed system is reanalysed delivering a different path OC. Once a “considerable” departure from the perfect fundamental path is observed, say at point C, he can stop the program and save current equilibrium state to the disk. He may then restart from C with the initial imperfection switched off. With a little bit of luck, one can hopefully expect to converge to a point D on the secondary path, from where he can then carry on in either direction. Although this sort of approach does not require any new solution procedures, it is apparently costly because of the reanalysis. Moreover, applying the right perturbation is by no means a trivial task and the success of the method is highly dependent on it. Generally speaking the perturbation can be either in the loading or in the structures. Its relevant position, orientation and relative intensity are governed by the buckling modes [56]. One way to avoid resolution is calculating the bifurcation point first without any perturbation. A fictitious force is then applied to the structure to force it into a secondary path. The force is adaptively adjusted during branching and it becomes zero once convergence is achieved. To make it work, one may

apply some well [102].

constraint

on

the

deformation

as

3.2. Direct method Branching by direct method is normally a two step process. A guess is made first to a solution on the secondary branch, it is then improved by some correction procedures. Convergence of such a process is apparently determined by both initial guesspredictor and iteration scheme-corrector. A good corrector may compensate for a poor predictor, or vice versa, but a bad predictor combined with a bad corrector is bound to fail. In this section the existing ways of constructing predictors and correctors are reviewed. Predictors. Based on the first eigenvectors: 4,, &, &, . . One simple approach that does not require any higher order differentiation than the K, is the eigenmode injection method [17,57], in which once the bifurcation point uc is located with its smallest eigenvector(s) 4, the displacement field is perturbed by scaled eigenvector(s).

The constants 5: represent the contribution by each eigenmode to the solution path. Generally they have to be tuned to get the correct value and converge on the secondary paths. An ad hoc selection procedure of c, has been mentioned in Wagner and Wriggers [57]. A similar approach is also adopted by Thurston et al. [55]. As a rule of thumb, {, should not be too small. so iteration does not converge back to the fundamental path. At the same time they cannot be too large so the predictor is still within the ball of convergence. Based on line searches. An interesting predictor based on line search has been suggested by Fujii et al. [58-601. Essentially he attempts to find an equilibrium state of a stationary energy along a certain direction at a fixed load level. To get the direction of the line search, he started with the incremental form of the equilibrium equation:

$dx=O

Table 2. Branching residuals for asymmetric bifurcation Eigenmode injection Path tangents

Iteration 0

Iteration 1

0.3874E-1 0.26OOE-1

0.2417E-4 0.6057E-5

J. Shi

Fig. 9. (a) Structural

response of the two bar truss-load u horizontal displacement. response of the two bar truss-load v vertical displacement.

where x=u, forj=l,..., n: x=1, forj=n+l. He then replaced the dth row of the above equation by dl = 0, which results in:

(7)

Table Bi-section Bracketing

inc.

13

3. Bracketing

where K,=(k,,,..;..k+

qi=(-$

Once this search direction tn = (dx, , dx,, . . . dx,) is found, a trial and error procedure is initiated to make

results of the van Mises truss

Determinant 3

(b) Structural

Smallest 3

pivot

Direct 4

1

Direct 4

2

211

Computing critical points and secondary paths G(u)ti = 0 along the line uc + {tl. uc is a numerical solution close to the bifurcation point and { is a constant that has to be adjusted to find a stationary energy. This is effectively making dII = 0, where II is the potential energy. As numerically it is very difficult to force C(u)tl to be zero, so in practice a small C(u)t, will be accepted. If d is such a point (see Fig. 3), then bd will form a good predictor, b is the bifurcation point. Notice well that along curve ac, structure is not in equilibrium except at point a and c. Obviously for the above method, choosing the right equation to be replaced by I = c needs some intuition as no solid guide is available. One drawback is that the initial value of 5 and its increment can only be guessed. If they are of the wrong magnitude, apparently line search will not get anywhere near the stationary point. Fujii’s approach is essentially the same as the one advocated by Seydel[30, p. 1691, except for the fact that Seydel does not perform any further line search and lets the corrector take care of a possibly bad predictor. As a consequence, Seydel employs an empirical formula to decide the magnitude of his predictor. The formula does not seem to be well grounded because it solely relies on the information of the primary path, which does not neccessarily give a clue of the secondary path. A similar energy search concept has been adopted by Kroplin et al. [61-64]. In his energy perturbation method, the line search direction (or perturbation pattern as he calls it) is the lowest eigenvector. Another major difference is that the amount of perturbation required to achieve a convergence on another equilibrium path can be decided beforehand. Based on the tangents to the branches. In concept, branch switching can be done as long as the bifurcation point and the tangents to different branches are known [l, 171. In theory this information can be obtained by continuously differentiating the equilibrium equation, but in practice the differentation process is expensive, if possible, otherwise it may be replaced by some approximation schemes, as given by Riks [l], Kouhai and Mikkola [65] and Eriksson [66]. To derive the tangents we first differentiate the original equilibrium equation (4) with respect to some path parameter, which will give:

c,li+cJ

=o.

Equation (8) can be differentiated again to yield the relation G,ii + G,l = -[G&i

+ 2GUlliA + G,,iA].

Since at the bifurcation point G,4 = 0, which implies G,+ = 0 as in eqn (8), hence eqn (10) left multiplied by @ will produce: @‘[G,,ii + 2GUlij + G,,M] = 0. If we substitute eqn (9) into eqn (ll),

where

a = 4TG,,~4 b = b’(G,,o

+ G&

c = dT(GU,~o + 2G,a

+ GA,).

The quadratic eqn (12) will give two sets of solutions of y, and yz. which provide us with the two tangents at the bifurcation point. One of them, say t, is the tangent to the primary path, the other t, to the secondary path:

f? = (ilz, 1,).

(13)

The tangent to the secondary path can consequently be employed to construct a predictor: u=Ar/t,.

(14)

The magnitude of the constant Aq is not a known a priori. Basically it has to be decided by trial and

error. A tentative guidance can be found in Riks [17]. If the bifurcation is symmetric, i.e. a = 0, then t2

(8)

fY20

(9)

we have:

where

y2 = 2 G,a +Gn=O.

(11)

(12)

Suppose (u, 1) is the critical solution and 4 is the zero eigenvector of K, then from eqn (8) we have: ;=r,+

(10)

Fig. 10. The von Mises Truss.

212

J. Shi

2x 6 t c OU : -I c) I

-3

Vertical

displacement

-2

(b)

‘.’ T

Load factor A

An&ytical 0.8

-0.8

Fig. 11. (a) Structural response of the von Mises truss-load c vertical displacement; (b) structural response of the von Mises truss--load u horizontal displacement.

reduces

to 4. This is in fact the justification

of the

eigenmode injection method mentioned earlier. The above derivation is general whether GU is not symmetric or loading is not proportional. However if the opposite is true, the previous equations can be greatly simplified. For example (0, 1) is in fact the tangent to the primary branch, hence it must satisfy eqn (11). So the constant c is actually zero. Moreover G, is only a function of u or G, is the constant loading pattern, consequently G, and G,, are zero. These effectively mean that

required to approximate the bifurcation point with a high accuracy and the evaluation of the second order derivatives [30, p. 1661. An approximation of these higher order terms has been devised [ 171.To simplify the process of calculating the constants, Riks [7] made t, orthogonal to t,: t;ts=t&J,fF+yp7)=0.

(15)

From which we can get: YJYI = -U(&)

4YG,,.,) and is zero. The major shortcoming of constructing a predictor with tangents is the excessive amount of computation

t, = yI(t,-- o/W)),

where yI is the magnitude of the predictor. t, can be approximated by the difference of the two solutions

Table 4. Bracketing results of the dome-first bifurcation point Bracketinn inc.

(16)

B§ion

Determinant

Smallest pivot

Direct 1

Direct 2

9

6

5

4

12 (limit point)

Computing critical points and secondary paths

213

Table 5. Bracketina results of the dome-second bifurcation noint Bracketing inc.

B§ion

Determinant

Smallest pivot

Direct 1

Direct 2

5

failed

6

5

12 (limit point)

on the fundamental path that enclose the bifurcation point. Duffeat [67] generalized the above approach by making tB at an angle to t, tT,t, = k.

(17)

Apparently this constant k is not known without any further complicated calculations. In practice it can be adjusted by trial and error to make t, closer to the real tangent to the secondary path. In his example problems, Duffeat claimed that k = 0 was always effective. Correctors with selective properties. In the last

section we have outlined some of the predictors. To guarantee a convergence onto the secondary path, we also need a good corrector as well as a close initial guess. In some senses a corrector can be regarded as a constraint that forces the iteration process in a certain direction or within a desired subspace. Guided mostly by geometrical argument, a number of correctors have been formulated. For example one can iterate “parallel” to the fundamental path [30, p. 1701:

f(u,1)= Au,-yk = const.

(18)

Here, one of the displacement components Au, is different from its counterpart y, on the primary path of a fixed constant. Since the corrector is “parallel” to the primary path, the chances of going back to it must be slim. Another heuristic constraint was proposed by Riks [17]:

f(u91) =

{u- YKU

- Y)T41t,f’

{u - y - [(u - y)Tt,]t,} - Aq*= 0, (19) which is basically a “cylinder” whose axis is the tangent t, of the primary path at the bifurcation point and whose radius is Aq. With a good predictor, we can expect two intersections on the secondary path with the cylinder. Another possible corrector: L = constant.

(20)

However it will fail in a situation when the secondary path does not intersect with Iz = constant. It is well known that for symmetric bifurcation, the lowest eigenvector Cp,is orthogonal to the fundamental path and parallel to the secondary path. As a result any increment Au from the bifurcation point with a nonzero projection on 4, must lead away from the primary path and possibly to the secondary path. Such a constraint was devised by AH-Chan [68]:

(AuTq5,)(c$;Au)-(Al)*= 0.

(21)

Other techniques. There are some other techniques which can not be easily fitted into the two categories [30,69]. Their importance should not be overlooked. Especially the paper by Kearfott [70] whose method is so general that it can handle both primary and secondary bifurcation points, whether they are simple or multiple. In addition, it does not require the computation of the second derivatives. However because of the limit on space, they will not be described here.

EA=104

Fig. 12. A three-dimensional dome.

214

J. Shi 4. OTHER RELATED RESEARCH AREAS

Multi-parameter system So far in this paper we have confined ourselves to single parameter system. In practice it is more likely that a nonlinear system has several control parameters. A multi-parameter system has been well addressed by applied mathematicians [23,91], a finite element formulation should be possible, in theory, for real engineering problems [92,93]. For such problems, the direct approach to calculate the singular point is of great value because the pragmatic way of varying only one parameter with the rest fixed, is simply too expensive to be of any practical importance.

Dynamic bifurcation, chaos? Dynamic stability is of tremendous importance to structural engineers. A large amount of work has been done on the theoretical side (see Herrmann [71] and Smites [72, 731 for instance). Some attempt has been made by the finite element analysts (for some of latest work the reader is referred to Refs [74-811). Because of the nature of commonly used solution schemes, test functions are not readily available. This is particularly true when it comes to the detection of Hopf bifurcation [82]. Recently dynamic stability has gained more and more attention because of its link with chaos, which is now deemed to be nonlinear the fundamental nature of many systems [30, 26-291.

Plastic bifurcation Bifurcation does not have to be elastic and conservative. It can also happen in plastic dissipative system like the necking and shear banding problems [45]. Since solution is now history dependent, bracketing must be done in one direction only along the solution path to avoid artificial unloading [94]. To find recent more general treatment on plastic bifurcation one is referred to Kubicek and Marck [95], Petryk and Thumann [96,97] and Nguyen and Triantafyllidis [98].

Multiple or clustered bifurcation Nearly all the theories presented in this report are only applicable to simple bifurcation problems. In practice multiple bifurcation may also occur. In such cases one can no longer use Det(K,), etc. as they may not change sign after passing the bifurcation point. Further to get the tangents to the bifurcated paths, one needs higher and higher derivatives than K, depending on the multiplicity of the bifurcation. Rabier and Rheinboldt [83] has recently done some work in this area (see also Kearfott [70] for some earlier work). The bifurcation points may also not necessarily be well separated. On the contrary they can be close to each other resulting in a scenario of clustered bifurcation. This can pose serious numerical difficulty as the bracketing and branching procedure have to be able to distinguish nearly identical bifurcation points.

5.

In the following examples trusses made of bars are considered for simplicity. Other types of structures with bifurcation can be found in Kouhia and Mikkola [65], Shi and Crisfield [33], Wagner and Wriggers [57] and Wriggers et al. [44,45]. One bar spring system We use a one bar spring system (see Fig. 4) to illustrate the idea of branching by unfolding. We perturb the original bifurcation problem shown in Fig. 4 by either an initial structure imperfection 6 or loading perturbation 6,. The first test was done on structural perturbation with 6 = 100, S, = 0. As one could expect because of the unfolding, we encountered a limit point which can be handled quite nicely by the arc-length method. After passing the limit point, the program was stopped with displacement, load factor, etc. written to a restart file. The program was then restarted without any imperfection, we can see clearly the “jump” from the imperfect solution to the perfect secondary path (Fig. 6). Although we only traced half of the secondary path, obviously the other

Symmetry and asymmetr) Even with the rapid growth of computer power, a full nonlinear finite element analysis is still sometimes too expensive. On the other hand there is often hidden (a)symmetry in structural responses [84]. which can be exploited to reduce computer storage requirement and increase solution efficiency. In this regard Carstensen and Wagner [90], Noor and Peters [85,86]. Werner and Spence [49] and Chang and Healy [87] have done some work. More recently Healy [88] and Ikeda and Murata [89] have proposed a more schematic approach in exploiting symmetry by group-theoretic and block diagonalization. Table 6. Branching

First bifurcation Second Third

bifurcation bifurcation

point point point

NUMERICAL EXAMPLES

parameters

of the dome

Step length

5”

i’,

52

2.0

0.0

I.0

0.0

1.0 1.0

0.0 0.0

1.0 0.2

5.0 1.0

1.0 1.0

0.1 0.0

0.2 0.1

1.0 5.0

215

Computing critical points and secondary paths

0

at -5.000

\

L

Fig. 13. Structural response of the dome--load

half can be followed equally well. A similar thing has been done with load perturbation 6 = 0, 6, = 100 (Fig. 7) although the “jump” is less obvious. Despite that all these worked nicely for simple test problems, they may very well break down for large complicated cases.

Two bar truss with asymmetric bifurcation This is a two-dimensional truss system exhibiting asymmetric bifurcation invented by the author to test branching procedures. Asymmetric bifurcation is essentially a reflection of asymmetry in structure or loading. In our case we can see that the inclined bar gives different stiffness when the joint is moving up and down. Firstly we show bracketing results of different schemes in Table 1. It is obvious that the direct methods took far fewer steps to compute the asymmetric bifurcation point. The direct 1, whereby the direct approach is invoked

L:vertical displacement at the apex

once the bifurcation is passed, is less efficient than the direct 2, which computed the bifurcation point from unloaded state without any normal path following. There were about three increments (depending on the initial load increment size, desired number of iterations, maximum increment size, etc.) before the bifurcation point was passed. All indirect methods and direct 1 had these steps and each step has about two iterations, but not direct 2. Strictly speaking, the “3” and “5” corresponding to direct 1 and 2, respectively, in Table 1 are iteration numbers rather than increment numbers. They should be treated differently from “11 ” , “4” and “4” of the indirect method, when comparing efficiency, as these are actual increment numbers with iterations in each increment. We will see later that it is advantageous to adopt direct 1, instead of direct 2, if there is more than one singular point present. Clearly, interpolation on either determinant or smallest pivot is more efficient than the bisection method. In general interpolations

Fig. 14. A three-dimensional

arch truss.

216

J. Shi

0.100 0.080 0.060 0.0&O 0.020 0.0

0

-0.020 -0.040 -0.060 -0.080 -0.100

t X E3

(b)

0.100

r

-6

0

-0.100

L

X E3

(c)

O.l.QOn

<

Fig. IS. (a) Structural response of the three-dimensional arch truss-load o vertical displacement at the top. (b) Structural response of the three-dimensional arch truss-load u x displacement at the top. (c) Structural response of the three-dimensional arch truss-load u y displacement at the top.

Computing critical points and secondary paths are superior to the bisection method if the singular point is simple. It should be emphasized here that the number of bracketing steps should be taken as relative, since they are influenced by the bracketing tolerance. After the bracketing, we can use either the eigenmode injection or path tangents to follow the secondary paths. Because this is a simple structure, there is very small difference between the two. The residual forces at the first two iterations are as in Table 2. Again the above numbers should be regarded as relative as they vary with step length. Finite element results are shown in Fig. 9, from which we can clearly see the asymmetry at bifurcation. The eon Mises truss

The geometry of this truss is made, such that its response shows both a limit and a bifurcation point. Bracketing results for the first bifurcation point is given in Table 3. Again the superiority of bracketing with interpolation is apparent. One special thing about this example is that the indirect bracketing results depend on the equilibrium convergence tolerance. A loose one (10m6) led to bracketing failure because the number of negative pivots was wrongly based. However a tighter tolerance (lo-“) overcame this problem and gave the results shown in Table 3. As for direct bracketing, we adopt the constraint in eqn (4) on the eigenmode as & = 1, if k = 3, we got the bifurcation point while k = 4 gave us the limit point. This is not surprising, as k = 3, we are forcing the eigenmode to have a unit component at the top node in x direction-a bifurcation mode; k = 4, we are making the eigenmode in the direction of primary path-y direction at the top node. So if we have some idea of the eigenmodes associated with different singular points, by constraining proper eigenmode components, we can compute, directly, one particular singular point. Branching was conducted with eigenmode injection, no difficulty was encountered. Results of load against vertical and horizontal displacement are shown in Fig. 11.

v

217

This example [58-601 has four singular pointsthree bifurcation points (including one double bifurcation) and a limit point on the primary path. They are fairly closely located, 1, = 8.68, 2.r= 10.26, A3= 15.67, 1, = 18.40. To avoid overstepping, the initial load step was taken as 3 and subsequent increment to be 1. Bracketing was performed on the first bifurcation point and bracketing performance is given in Table 4. The advantage of using interpolation or direct 1 is clearly demonstrated in the above table. It should be noted that because of the existence of four singular points, the direct 2 actually computed the limit point instead of the first bifurcation point. Due to the same reason, one has to ensure that the increment size should be small enough that it will not step over two singular points in one increment. The same excise has been repeated on the second singular point, a double symmetric bifurcation point, which has made the interpolation methods inferior to the bi-section method, as we could not consistently interpolate on the same negative pivot out of two. In this particular case, since the nullity is even, the determinant does not change sign, which made it impossible to interpolate with the determinant. The performance of interpolation deteriorated more when attempting to compute the third bifurcation point. There are now a maximum of five negative pivots, so confusing that bracketing failed for pivot based interpolation. The third bifurcation point is a double one, so determinant based interpolation could not cope with it. The direct 2 always gave the limit point, if we use the constraint in eqn (4) on the eigenmode as 4, = 1, where k is fixed. However if we constrain another k we then get a different singular point, say the third bifurcation point. Branch switching at the first bifurcation point is easy, while for the second and third double bifurcation points, a lot of “tuning” has to be done to find the &s in eqn (5) and step length. These are given in Table 6.

I Fig. 16. A two-dimensional circular arch.

218

J. Shi

Vertical

Fig. 17. Structural

response

of the two-dimensional

arch-load

The dome response, load vs vertical displacement at the apex is given in Fig. 13. Three-dimensional

rectangular

arch

This is another example that has several singular points (four bifurcation points R, = 12.28, d, = 69.13, i., = 99.26, I, = 22.94) on the equilibrium curve. but they are all simple bifurcation points. As a result no particular problem was encountered in either bracketing or branching (except at the last bifurcation point, it is difficult to follow the whole secondary path). Similarly regarding the dome, one should make sure that the increment size is not too big. The branching step lengths for the first three bifurcation points are: 1.0,2.0 and 0.5, respectively. The three displacements at one of the top nodes under vertical loading are given in Fig. 15. At the first and fourth bifurcation points, the structure buckle in YZ plane (Fig. 15b), i.e. no displacement in x direction; at the second and third bifurcation, it buckles in xz plane (Fig. 15~) i.e. _r displacements are symmetric at the top nodes. This is the only example that has bifurcation points after limit point. TWO-dimensional

circular

3.

4.

5.

6.

7.

8.

9.

arch

This is a more realistic problem with a limit and bifurcation point. The geometric and material data are given in Fig. 16. Bracketing on the bifurcation point (iv = 6.29) and branching had been conducted without any difficulty. Its response (Fig. 17) is different from the von Mises’s truss in that secondary paths do not rejoin the primary path.

10. II.

12.

REFERENCES

13. 1. E.

Riks, Some computational aspects of stability analysis of nonlinear structures. Comput. Meth. appl. Mech. Engng 47, 219-260 (1984) 2. L. Azrar et al. Asymptotic numerical methods to

L; vertical

displ. at apex

displacement

at the apex

compute bifurcating branches. In: Proc. Eur. Conf. on New Adtlances in Computer Slructure and Mechanics, Giens, France (Edited by P. Ladeveze and 0. C. Zienkiewicz) (1991). R. C. Batista. R. C. Antonini and R. V. Alves, An asymptotic model approach to nonlinear strucutral elastic instability. Comput. Struct. 38, 475484 (1991). R. T. Haftka, R. H. Mallet and W. Nachbar, Adaptation of Koiter’s methods to finite element analysis of snap-buckling behaviour. Int. J. Solids. Strucl. 7, 1427-1447 (1971). R. H. Mallet and R. T. Haftka, Progress inponlinear finite element analysis using asymptotic solution technique. In: Advances in Computer Methods in Structures anh Mechanics (Edited by J. T. Oden er al.), DD. 357-374. UAH Press. The Universitv of Alabama in Huntsville, AL (1972). E. Riks. Application of Newton’s method to the problem of elastic stability. J. appl. mech. 39, 1060-1065 (1972). E. Riks, An incremental approach to the solution of snapping and buckling problems. Inl. J. Solids Strucf. 15, 529-551 (1979). M. A. Crisfield, New solution procedures for linear and nonlinear finite element analysis. In: The Mathematics of Finite Elements and Applications V (Edited by J. R. Whiteman), pp. 49981. Academic Press, London (1985). R. Kouhia, Generalised Newton-Raphson Techniques. In: Nonlinear Dynamics, Chap. 7 (Edited by K. S. Virdi). Blackie Academic and Professional, London (1994). W. C. Rheinboldt, Numerical Analysis of Parametericed nonlinear Equalions. Wiley, New York (1986). J. M. Clarke and G. J. Hancock. A study of incremental-iterative strategies for nonlinear analyses. Int. J. numer. Meth. Engng 29, 1365-1391 (1990). M. A. Crisfield and J. Shi. 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 P. Wriggers and W. Wagner), pp. 47-68. Springer, Berlin-71991). E. Riks. On formulations of oath-following techniques for structural stability analysis. In: Nena Advances in Computational Structural Mechanics (Edited by P. Ladeveze and 0. C. Zienkeiwicz). Elsevier, Oxford (1992).

Computing

critical

points

14. G. Skeie and C. A. Fehppa, Detecting and traversing bifurcation points in nonlinear structural analysis. Int. J. Space Struct. 6 (1991). 15a. E. L. Allgower and K. Georg, Continuation and path following. Acta Numer. 164 (1993). 15b. M. Papadrakakis, Solving large-scale nonlinear problems in solid and structural mechanics. In: Solving Large-scale Problems in Mechanics (Edited by M. _ Papadrakakis). Wiley, New York (1992). 16. E. Riks. Proeress in collaose analvsis. _ NLR MP 84031 U (1984). and stability, a numerical 17. E. Riks, 1984, Bifurcation approach. In: Innovative Methods for Nonlinear Problems (Edited by W. K. Liu et al.), pp. 313-344. Pineridge. Swansea. U.K. (1984). 18. D. J. Allman, On the general theory of the stability of equilibrium of discrete conservative system. Aero. J. 93, 29935 (1989). Theory of buckling and post-buckling 19. B. Budiansky. behaviour of elastic structures. In: Advances in Applied Mechanics, Vol. 14. Academic Press (1974). (Ed.). Buckling of structures. In: IU20. B. Budiansky TAM Symposium Cambridge. U.S.A. June 1974. Springer, Berlin (1976). and W. T. Koiter, Postbuckling 21. J. W. Hutchinson theory. Appl. Mech. Rev. (1970). Plastic buckling. In: Advances in 22. J. W. Hutchinson, Applied Mechanics, Vol. 14. pp. 67--141. Academic Press, London (1974). 23. K. Huseyin, Multiple Parameter Stability Theory and its Applications. Clarendon Press, Oxford (1986). 24. A. Komarakul-Na-Nakorn and J. S. Arora, Stability criteria: a review. Comput. Sfruct. 37, 3549 (1990). 25. E. Ramm (Ed.). Buckling of Shells. In: Proc. Stale-ofthe-Ar/ Colloquium, University of Stuttgart. Springer, Berlin ( 1982). 26. J. M. T. Thompson, Instabilities and Catastrophies in Science and Engineering. Wiley, New York (1982). 27. J. M. T. Thompson and G. W. Hunt (Eds), COLLAPSE: the buckling of structures in theory and practice. In: Proc. IUTAM Symp. Cambridge University Press (1982). 28. J. M. T. Thompson and G. W. Hunt, Elasric Instability Phenomena. Wiley. New York (1984). 29. J. M. T. Thompson and H. B. Stewart, Nonlinear Dynamics and Chaos. Wiley. New York (1986). 30. R. Seydel, From Equilibrium to Chaos. Elsevier, Oxford (1988). 31. D. M. Greig, Optimisation. Longman, Harlow (1980). 32. G. R. Walsh. Method of Optimization. Wiley, New York (I 975). 33. J. Shi and M. Crisfield, A semi-direct approach in the computation of singular points. Inr. J. Comput. Struct. (in press). 34. R. Kouhia. private communication. 35. M. Fujikake, A simple approach to bifurcation and limit point calculations, 21, 183-191 (1985). 36. S. H. Lee and D. N. Herting, Comments on “A simple approach to bifurcation and limit point calculations. Inr. J. numer. Meth. Engng 21, 193551937 (1985). 37. S. L. Chen, A non-linear numerical method for accurate determination of limit and bifurcation points. Int. J. numer. Meth. Engng. 36, 277992790 (1993). 38. Z. Waszczyszyn. Numerical problems of nonlinear stability analysis of elastic structures. Comput. Strucr. 17, 13-24 (1983). 39. F. Tabaddor and J. Padovan, Element nonlinearity effects on local and global stability. Compur. Struct. 40, 401407 (1991). 40. H. Mang. On special points on loadchsplacement paths in the prebuckling domain of thin shells. Int. J. numer. Meth. Engng. 31, 207-228 (1991). 41. H. Mang, On bounding properties of eigenvalues from

and secondary

42.

43.

44.

45.

46.

47.

48.

49.

50. 51.

52.

53.

54.

55.

56.

57.

58.

59.

60.

61.

62.

paths

219

linear initial FE stability anlaysis of thin, elastic shells with respect to stability limits from geometrically non-linear prebuckling analysis. Int. J. numer. Meth. Engng 31, 1087-I 111 (1991). D. H. Tao and E. Ramm, Characteristics of subspace in reduced basis technique. Int. J. numer. Merh. Engng 31, 1567-1583 (1991). W. C. Rheinboldt and E. Riks, Solution techniques for nonlinear finite element equations. In: State-of-lheArt surveys on Finite Element Technology (Edited by A. K. Noor and W. D. Pilkey). ASME, New York (1983). P. Wriggers, W. Wagner and C. Miehe, A quadratically convergent procedure for the calculation of stability points in finite element anlaysis. Comput. Mefh. appl. Mech. Engng 70, 3299347 (1988). P. Wriggers and J. C. Simo. A general procedure for the direct computation of turning and bifurcation points. Im. mimer. Meth. Engng 30, 155-176 (1990). J. P. Abbott, Numerical continuation methods for nonlinear equations and bifurcation problems. Ph.D. dissertation. Australian National University (1977). J. P. Abbott, An efficient algorithm for the determination of certain bifurcation points. J. Comput. appl. Math. 4, 19-27 (1978). G. Moore and A. Spence, The calculation of turning points of nonlinear equations. SIAM J. numer. Anal. 17 (1980). B. Werner and A. Spence. The computation of symmetry-breaking bifurcation points. SIAM J. numer. Anal. 21 (1984). G. W. Reddien, On Newton’s method for singular problems. SIAM J. numer. Anal. 15, 993-996 (1978). D. W. Decker. H. B. Keller and C. T. Kellv. Convergence rates for Newton’s method at singmar points. SIAM J. numer. Anal. 20, 296314 (1983). A. Griewank, On solving nonlinear equations with simple singularities or near singular solutions. S1.4 M Rev. 27, 5377563 (1985). E. Rothwell and B. Drachman, A unified approach to solving ill-conditioned matrix problems. Inr. J. numer. Melh. Engng 28, 6099620 (1989). E. Riks and C. C. Rankin, Bordered equations in continuation methods: an improved solution technique. NLR MP 87057 U (1987). G. A. Thurston. F. A. Brogan and P. Stehlin, Postbuckling analysis using a general-purpose code. J. AIAA 24, p. 1013 (1986). M. Talbot and G. Dhatt. Three discrete Kirchoff elements for shell analysis with large geometrical nonlinearities and bifurcations. Engng Compur. 4 (1987). W. Wagner and P. Wriggers. A simple method for the calculation of postcritical branches. Engng. Comput. 5, 103 (1988). K. K. Choong and F. Fujii, Line search for asymmetric bifurcation path of elastic structures. In: Proc. of Symp. on Computational Methods in Structural Engineering and Related Fields, Vol. 14. pp. 73378. Tokyo (1990). F. Fujii and K. K. Choong. Branch switching in bifurcation of structures. J. Engng Mech. ASCE 118, (1992). F. Fujii and K. Asada, Branch switching in simple spatial bifurcation models. Proc. Seiken-IASS s.ymp.. 19-22 October, Tokyo, Japan. B. Kroplin, D. Dinkler, and J. Hillman, An energy perturbation applied to nonlinear structural analysis, Comput. Meth. appt. Mech. Engng 52, 885-897 (1985). B. Krophn and D. Dinkler, Dynamic versus static buckling analysis of thin walled steel structures. In: Finite Element Method for Plate and Shell Structure, Vol. 2, Formulation and Algorithms (Edited by T. J, R. Hughes and E. Hinton). Pineridge Press, Swansea (1986).

220

3. Shi

63. B. Kroplin and D. Dinkler, On static and dynamic instability analysis of thin-walled structures. In: Discretizatih Methods in Structural Mechanics (Edited by G. Kuhn and H. Mana) IUTAM/IACM Svmoosium. Vienna, Austria, Springer, Berlin’ (1989). - * 64. B. Kroplin and D. Dinkler. An Energy-based concept of structural stability. In: Proc. Eur. Conf. on New Advances in Computer Structure and Mechanics, Giens, France (Edited by P. Ladeveze and 0. C. Zienkiewicz) (1991). 65. R. Kouhia and M. Mikkola, Tracing the equilibrium path beyond simple critical points. Int. J. numer. Meth. Engng 28, 2923-2941 (1989). 66. A. Eriksson, Derivatives of tangential stiffness matrices for equilibrium path descriptions, Int. J. numer. Merh. Engng 32, 1093-l 113 (1991). 61. G. A. Duffeat. Some aspects of the numerical solution of equilibrium problems in finite elasticity. Ph.D thesis, Department of Civil Engineering, University of Cape Town (1985). 68. H. E. AH-Chan, The large displacement behaviour of thin plates, Ph.D. thesis, Department of Structural Engineering, University of New South Wales (1986). 69. W. C. Rheinboldt. Numerical methods for a class of finite dimensional bifurcation problems. SIAM J. numer. Anal. 15 (1978). 70. R. B. Kearfott. Some general Bifurcation Techniques. SIAM J. Sri. Stat. Comput. 4 (1983). 71. G. Herrmann (Ed.), Dynamic Sta6iIit.y of Structures. Pergamon Press, Oxford (1965). 72. G. J. Simitses, Instability of dynamically loaded structures. Appl. Mech. Rev. 40 (1987). 73. G. J. Simitses, Dynamic Stability qf Suddenly Loaded Structure. Springer, Berlin (1989). 74. M. Kleiber, W. Kotula and M. Saran, Numerical analysis of dynamic quasi-bifurcation. Engng Comput. 4, p. 48 (1987). 15. M. Kleiber et al., Finite element anlaysis of dynamic quasi-bifurcation in elastic-plastic structures. Comput. Struct. 34, 6499652 (1990). On the dynamic 76. W. Auli and F. G. Rammerstorfer, instability of shell structures-riteria and algorithms. In: Finite Element Method for PIate and Shell Structure, Vol. 2, Formulation and Algorithms (Edited by T. J. R. Hughes and E. Hinton). Pineridge Press, Swansea (1986). 77. M. A. Crisfield, J. Shi and K. L. Lim, Finite element and nonlinear dynamics. In: Modern Practice in Stress und Vibraiton Analysis. Stress Analysis Group of the Institute of Phvsics. Sheffield (1993). finite element 78. K. L. Lim and ‘M. A. Crisfield,‘Dynamic analysis applied to a simple model exhibiting dynamic instability. Engng Compu/. (submitted). 79. P. Wriggers and C. Carstensen, An efficient algorithm for the computation of stability points of dynamical systems under step loading. Engng Comput. 19, 6699679 ( 1992). and P. Wriggers, On perturbation be80. C. Carstensen haviour in non-linear dynamics. Commun. numer. Meth. Ennng 9, 165-175 (1993). 81. K. H. Wmters; K. A. Cliffe and C. P. Jackson, The prediction of instabilities uisng bifurcation theory. In: Numerical Methods in Transient and Coupled Problems (Edited by R. W. Lewis). Wiley. New York (1987).

82. T. J. Garratt, G. Moore and A. Spence, Two methods for the numerical detection of Hopf bifurcations (1991). 83. P. J. Rabier and W. C. Rheinboldt, On a computational method for the second fundamental tensor and its application to bifurcation problems. Numer. Math. 57, 681694 (1991). 84. G. W. Hunt, Hidden (a)symmetries of elastic and plastic bifurcation. Appl. Mech. Rm. 39 (1986). 85. A. K. Noor and J. M. Peters, Tracing post-limit-point paths with reduced basis technique. Comput. Meth. appl. Mech. Engng 28, 217-240 (1981). 86. A. K. Noor and J. M. Peters, Recent advances in reduction methods for instability analysis of structures. Comput. Struct. 16, 6780 (1983). 87. P. Chang and T. Healey, Computation of symmetry modes and exact reduction in nonlinear structural anlaysis. Compuf. Struct. 28, 135-142 (1988). 88. T. J. Healey, A group-theoretic approach to computational bifurcation problems with symmetry. Comput. Meth. appl. Mech. Engng 67, 257-295 (1988). 89. K. Ikeda and K. Murota, Bifurcation analysis of symmetric structures using block-diagonalization. Cornnut. Meth. appt. Mech. Engng 86, 215-243 (199 1). 90. C. darstensen and W. Wagner, Detecting symmetrybreaking bifurcation points of symmetric structures. Int. J. Numer. Meth. Engng 36, 3019-3039 (1993). 91. A. Spence and A. Jepson, Numerical Techniques for Multi-Prameter Problems. Numerical Analysis, Lecture Notes in Mathematics 1066. Springer, Berlin (1984). 92. W. Wagner and P. Wriggers, Calculation of bifurcation points via fold curves. In: Nonlinear Computational Mechanics: State of the Art (Edited by P. Wriggers and W. Wagner), pp. 69-84. Springer, Berlin (1991). Fold lines for sensitivity analyses in 93. A. Eriksson, structural instability. Comput. Meth. appl. Mech. Engng (submitted). and H. Zhang, On the 94. W. F. Chen, E. Yamaguchi loading criteria in the theory of plasticity. Comput. Struct. 39, 679683 (1991). 95. M. Kubicek and M. Marck, Computational Methods in Bifurcation Theory and Dissipative Structure. Springer, Berlin (1984). 96. H. Petryk, General Theory of Bifurcation and Stability in Time-Indeuendent Plasficity. Lecutre notes, CISMcourse, Udine. Springer, Berlin (in press). 97. H. Petryk and K. Thermann, On discretized plasticity problems with bifurcations. 29, 745-765 (1992). 98. S. Q. Nguyen and N. Triantafyllidis, Plastic bifurcation and- postbifurcation analysis for generalised standard continua. J. Mech. Phvs. Solids 37, 545-566 (1989). 99. C. A. Fellippa, Transversing critical points with penalty springs. In: Transient/Dynamic Analysis and Constitutive Laws for Engineering Material (Edited by Pande and Middelton), pp. C2/1lC2/8. Nijhoff, Dordrecht, (1987). 100. C. T. Helley and R. Suresh, A new acceleration method for Newton’s method at singular points. SIAM J. Numer. Anal. 20, 1001-1009 (1983). of second order 101. J. Shi and M. A. Crisfield, Application derivatives in path following (in preparation). by specific displacement 102. J. Shi, Branch switching control and adaptive fictious stiffness. Note (1992).