Efficient algorithms for parametric non-linear instability analysis

In this paper, we describe an original method which enables the rapid calculation of the critical buckling pressure of structures with variable parameters. This study ,ts within the framework of non-linear reliability which requires numerous non-linear mechanical calculations with di.erent sets of parameters. Therefore, we aim for a method which has optimum numerical e(cid:29)ciency.


Introduction
This work was performed in the framework of reliability analysis of the buckling strength of a thin structure subjected to external pressure. The critical buckling load is the maximum load admissible before the structure becomes unstable. We attempt to predict how uncertainties on the geometry, the materials or even the loading a ect the response of the structure.
Today, a majority of the reliability calculations with respect to the buckling load of a structure use a linear analysis to determine Euler's critical load c (Fig. 1). Imperfections on structural parameters are known to lower the critical buckling load from c to 1 or 2 .
In order to take into account geometrical defects as well as elastic-plastic constitutive relations, it is necessary to perform a quasi-static non-linear analysis.
The study of the sensitivity of the response to the di erent imperfections can be undertaken in several di erent ways.
The ÿrst class of methods is based on the calculation of the gradients of every variable [5,6]. One performs a mechanical calculation for each new set of parameters.
The second class of methods tries to evaluate the evolution of the limit point as a function of the parameters [7,8]. They consist of following the evolution of this limit point with regard to the parameters using an augmented system (see details in Section 2.2.4).
The approach described here is half-way between the two preceding methods. It is based on Riks' incremental non-linear calculation algorithm [9] (see details in Section 2.1) in which the tangent matrix is used to solve the equilibrium equations. The assembly and Crout decomposition of this matrix for each load increment represents a large percentage of the computation time, particularly for three-dimensional structures. In its ÿrst step, the strategy adopted consists of reutilizing the Crout decompositions of the tangent matrices for problems with mean values of the parameters in order to calculate the solution of the problem with modiÿed parameters (see details in Section 2.2.1).
In the second step, in order to calculate a minimum number of points to ÿnd the limit point, we introduce the total Lagrangian formulation into the code. This enables us to calculate the portion of the curve preceding the critical point in a single load step (Section 2.2.2).
Thus, the coupling with the reliability analysis is not arbitrary, i.e. we use for a modiÿed calculation the data previously calculated and saved for the reference values. If necessary, this reference can be changed in the course of the reliability iterations.
In Section 3.1, we present the original ÿnite elements COMU and COMI which are particularly well-adapted to non-linear buckling analysis of thin axisymmetric structures.
Finally, we present two applications of the methods developed here: a simple ring (Section 3.2) and a sti ened ring (Section 3.3) under external pressure.

The Riks' method for the basic curve
The implicit updated Lagrangian algorithm based on Riks' method [9] involves two loops: the ÿrst to increment the loading, the second to solve the equilibrium equations. The calculation is performed under the assumption of quasi-static evolution.
• n designates the index of the ÿrst loop, • i designates the index of the second loop, • U designates a displacement vector, • designates the Lagrange multiplier of the loading, • is an increment, • ÿ designates the arc length, • K T designates the tangent matrice.

Evolution loop
This part of the algorithm is used to initialize the values for the equilibrium calculation. The external forces are calculated on the conÿguration of Step n; they are used to obtain a "linear" estimate of the displacement increment for the next step n + 1.

Equilibrium solution
The resolution of the equilibrium equations is based on a linear prediction of the displacement increment corrected by a non-linear increment [10].
The load and displacement increments are steered according to an arc-length type strategy [11]. The equations which must be veriÿed are with and In Eq. (2), which increments the non-linear part with the residual, the notation K −1 T is symbolic: it means that one actually solves the system by Crout decomposition. System (1) has the exact solution: The internal and external forces are calculated in updated Lagrangian mode on the conÿguration of the previous iteration i. The arc-length increment is given by the variable ÿ. It is initialized using the linear part of the ÿrst calculation step.
This value is then increased or reduced automatically taking into account the number of iterations performed until convergence. Convergence is reached as soon as the residual ||R i || becomes small enough.
We apply to the structure a geometric imperfection of the same shape as the buckling mode. This imperfection allows us to take into account the geometric defects of the structure and to follow the post-buckling bifurcated branch easily. The response of a modiÿed structure is usually little di erent from that of the reference structure.
Furthermore, the tangent matrix used in the resolution of the equilibrium equations can be approximate because it is not used in the convergence criterion for the equilibrium. These two remarks inspire us to reutilize for all modiÿed calculations the Crout decompositions of the tangent matrices previously stored during the reference calculation (Fig. 2). Thus, we obtain the modiÿed curve without assembling or inverting any system of equations. On the other hand, the number of iterations to reach convergence in the equilibrium loop is larger. The choice of the matrix to use for each loading step is based on the cumulated arc-length of the calculation. One seeks the step of the reference calculation whose cumulated arc length value is closest to that of the current modiÿed calculation step.

Calculation of a portion of the modiÿed non-linear curve
Within the framework of reliability, it is not necessary to know the complete response of the structure: we are interested only in the limit point. In order to save computation time, our approach is to begin the modiÿed calculation at an arbitrary point on the reference curve. For this ÿrst point, we calculate the internal and external forces in total Lagrangian mode with respect to the undeformed modiÿed structure (see Fig. 3).
Once the starting point on the reference curve has been chosen, we steer the algorithm di erently in order to converge toward a point on the modiÿed curve. The choice of the starting point on the reference curve is detailed in Section 2.2.3. We use the following notations: • N : index of the starting point on the reference curve • U N r : total displacement of the reference structure at the starting point • U N r : displacement increment of the reference structure at the starting point • : loading increment to go from the reference structure to the modiÿed structure We initialize the calculation as follows We seek a displacement correction V as well as a loading parameter correction such that the cumulated arc-length of the new displacement is the same as that of Step N of the reference calculation: Thus, the steering during the ÿrst total Lagrangian calculation step can be summarized as The system of equations (4) is a ÿrst-degree equation in V i+1 and second-degree equation There are two solutions to this second-degree equation, of which we always choose the smaller. The non-linear part of the displacement is deÿned by We need to calculate the internal and external forces on the undeformed conÿguration of the modiÿed structure 0 m by the total Lagrangian approach. Then, we can calculate the residual It ensures convergence of the equilibrium point once the residual ||R i || gets small enough. We can summarize this in the following algorithm:

A general strategy to ÿnd the buckling point on the modiÿed curve
The search for the critical buckling point on the modiÿed curve can be performed in several ways.

2.2.3.1.
Complete calculation of the modiÿed curve. We use the method in Section 2.2.1 and test at the end of each step whether the load parameter is decreasing. The calculation is interrupted as soon as unloading is detected (Fig. 4). This simple and numerically reliable method is still calculation-intensive because we recalculate the whole modiÿed curve.

2.2.3.2.
Calculation of a portion of the modiÿed curve. Since the method proposed in Section 2.2.2 allows us to begin the calculation at an arbitrary point on the modiÿed curve, we are looking only for the peak (Fig. 5). How do we manage to begin the calculation before the limit point? 2.2.3.3. Approximate search for the peak. We calculate by dichotomy the three highest points of the curve with 3 di erent starting points on the reference curve (Fig. 6). This gives only an approximation of the maximum, but this can be su cient  for ÿrst reliability calculations where an error on the mechanical calculation is acceptable.

Exact search for the peak.
We can improve the previous method by using the point calculated in total Lagrangian mode which is located just before the maximum, then continuing the calculation in updated Lagrangian steps (Fig. 7).

Application to elastic-plastic buckling cases and comparison with Erikson's and Cochelin's methods
The methods developed by Erikson et al. in [7] and by Cochelin in [8] consist of identifying the limit point of the mean value problem, then following the path of this limit point in function of the parameters using an augmented system: The method in [7] is based on a classical ÿnite element calculation, whereas [8] uses an asymptotic numerical calculation (see [12]). Both methods are e cient for buckling in the elastic domain. Their drawback is that they do not take the loading history into account and, therefore, they cannot be used for plastic buckling.
Moreover, these two methods developed necessitates to develop the analytical derivation of the matrices with respect to every variables. This is very heavy to develop. Our method does not need such a development.

E ciency considerations
The time saved for the complete calculation of the modiÿed curve is on the order of 50% for structures with many degrees of freedom. It is about 70% if one calculates only the peak of the curve (see Section 3.3.5). The number of iterations to obtain a point on the modiÿed curve starting from the reference curve is approximately equal to that needed to get from one point to the next in the updated Lagrangian scheme.
However, this is very dependent on the structure studied: the more non-linear the portion of the curve preceding the limit point, the higher the gain.

The COMU and COMI ÿnite elements
These two ÿnite elements are general conical axisymmetric shell elements. They were developed to predict the non-linear response of axisymmetric shells subjected to any type of loading (axisymmetric or not) and either geometric initial imperfections or non-axisymmetric thickness initial imperfections. These elements have been proved to be very e cient if the initial imperfection can easily be decomposed into a few Fourier harmonics (say 5 -8). When the imperfection is highly localized, a 3D shell representation is more e cient. We will now brie y describe the formulation of the elements. Both elements are conical Mindlin shell elements with two nodes. The displacement ÿeld is decomposed into a Fourier series. For the  COMI element, the geometry is assumed to be imperfect in terms of coordinates (the mean radius is non-constant around the circumference) as well as thickness [13]. For the COMU element [14], the initial imperfection is decomposed into a Fourier series. For the COMI element, the imperfections (for the geometry as well as for the thickness) are assumed to be given for each element on a set of points -designated by Â i -around the circumference. In principle, the COMI element can handle any type of imperfect axisymmetric structure and, therefore, should enable one to avoid any 3D computation at all for axisymmetric structures. However, when the imperfections or the loads are localized, the e ciency of the element is poor compared to a full 3D shell analysis due to the Fourier decomposition of the displacement ÿeld.
3.1.1. Formulation 3.1.1.1. Notations and reference frames. Let us introduce a cylindrical coordinate system: each point M of the shell has three coordinates (r; Â; z). The local reference frame at point M is deÿned by (ñ;s;t), whereñ is the normal to the shell,s is the tangent along the meridian direction andt the tangent along the circumferential direction. In this local system, the displacement vector will be represented by a 4-degrees-of-freedom vector U with 3 translations (w; u; v) and one rotation ÿ around directiont.  Let us designate by the angle (−z;s). As usual for axisymmetric structures, we decompose the displacement ÿeld U into a Fourier series.
3.1.1.2. Strains. We will now express the strains in the local reference frame as a function of the displacement ÿeld in the same frame. For a conical shell with no initial imperfections, the linear strains are, as usual, separated into membrane ( lm ), bending ( lb ) and shear ( s ) strains given by r 2 (− @ 2 w @Â 2 − r sin @w @Â + cos @v @Â ) lb sÂ = 1 2r 2 [ − 2r @ 2 w @s@Â + 2r cos @v @s − v cos sin + sin @w @Â ] The quadratic membrane strains q (U ; U ) used in large displacement analysis as well as in buckling analysis are given by @w @s + ( @u @Â − v sin ) @u @s + ( @v @Â + w cos + u sin ) @v @s ] 3.1.1.3. Strains for the geometrically imperfect shell element. We now introduce a geometric initial imperfection represented by the displacement vector d 0 from the perfect structure to the imperfect one. The displacement ÿeld can be discretized either by its expression in a Fourier basis for the COMU element or on the set of discrete points Â i for the COMI element. The strain ÿeld is now lm (d 0 ; U ) = lm (U ) + q (d 0 ; U ); In this expression, q (d 0 ; U ) is obtained simply by replacing one of the vectors U by d 0 in the expression of the quadratic strains.
3.1.1.4. Stress-strain law. We simply use Hooke's plane stress law to deduce the membrane and bending stresses from the corresponding strains. The shear stress is obtained by s = 5 6 G s ; where G is the shear modulus. In the elastic-plastic case, an Illiushin "global" model is chosen [15].
3.1.2. The COMU and COMI ÿnite elements 3.1.2.1. Geometry. The element is a two-node conical element. The discretization of the geometrical imperfection di ers between the two ÿnite elements: for the COMU element, the imperfection is given as coe cients of the Fourier series at each node, whereas for the COMI element it is given at each node as a set of local displacements for each point Â i . For the COMI element, the thickness ÿeld is also deÿned on the set of points Â i . In the case of large displacement analysis, the geometry is updated at each load increment.  with the derivative of the normal displacement w with respect to s.

Integrations and derivations.
Numerical integration with one Gauss point at the center of the element is used along the meridian direction. The integration around the circumference is analytical in the case of the COMU element and numerical (using Simpson's rule) for the COMI element. The derivation along the circumferential direction (coordinate Â) is analytical for the COMU element and numerical for the imperfection ÿeld and for the COMI element.

Non-linear calculations
The geometrically non-linear strategy is simply the updated Lagrangian technique. Riks' arc-length control [16] is used to proceed beyond snap-through or bifurcation points. The ring is meshed simply with a COMU element on a Fourier basis decomposed on modes (0,2,4). The perturbations applied in the calculation of the modiÿed curve are −10% on the Young's modulus, the thickness, the amplitude of the defect and the yield stress. This example has been treated by several authors. Wang [17] calculated a ring with a circumferential hinge. He showed that the critical buckling load drops signiÿcantly compared to the conÿguration without hinge. Fu and Waas [18] studied a ring made of a composite material. More recently, Huddleston and Sivaselvan [19] analyzed the buckling of a compressible ring. Fig. 9 shows the result obtained by calculating the modiÿed curve completely with the method in Section 2.2.1.

Part of the modiÿed curve
We applied the method in Section 2.2.2. Only nine points su ce to ÿnd the limit point of the modiÿed curve (Fig. 10).

Deÿnition of the problem
The results presented here were calculated on a sti ened ring (see Figs. 11 and 12), a structure whose non-linear behavior is more complex. This example was carried out by Bourinet et al. [20]. The mesh comprises 80 COMU elements (see Section 3.1). The Fourier basis chosen is (0,2,4,6,8). The geometric defect is deÿned on the second harmonic, its amplitude is 9:952 mm. The geometry and the

Complete modiÿed curve
The perturbation is introduced on the exterior thickness: −17%. Fig. 13 shows, on the example of the sti ened ring, the shape of the total radial displacement of the shell node. Fig. 14 shows the shape of the total radial displacement after a total Lagrangian calculation of the ÿrst step (see Section 2.2.2).

Quality of the response
As can be seen from the curves, the results of the calculation with the e cient method are identical to those of the direct calculation of the modiÿed curve. One can thus appreciate the quality of the method when applied to this example.

E ciency
On this example, there is saving of time on the order of 40% when the modiÿed curve is calculated by the method given in Section 2.2.1. When only the upper part of the curve is calculated in order to obtain the limit point, the gain is 60%. The order of magnitude of the reference computation time for this example is 5 -10 min.

Conclusion and perspectives
We have introduced an original method of parametric non-linear calculation in order to provide, with reliability analysis in mind, a capability to perform successive calculations for di erent, but similar, sets of parameters. The key to this method is the reutilization of inverted tangent matrices saved from the mean value problem. We have also introduced the possibility of using a total Lagrangian formulation in performing the ÿrst step of the calculation of the modiÿed curve.
These ideas are inapplicable if the responses of the reference structure and of the modiÿed structure are very di erent, particularly in cases where the primary buckling modes have di erent shape. If so, one can simply change reference calculation.
For plastic buckling, the most e cient method, i.e. calculation of the peak only, cannot be used because the non-linear path depends on the loading history. One must choose a starting point on the reference curve such that the target point on the modiÿed structure is only very slightly plastic.
The use of COMU and COMI elements enables one to mesh axisymmetric structures and to seek the displacement in a given, non-axisymmetric Fourier basis. The geometric defect is given in the Fourier basis (COMU) or point-by-point on the circum-ference (COMI). By using these elements one can achieve a further gain in e ciency for modal defects.
In certain cases these coupled methods allow a gain of more than an order of magnitude in the computation time in non-linear buckling compared to a conventional ÿnite element code with three-dimensional shell elements. This evaluation is in progress.
The e cient method developed here can be applied directly to the three-dimensional cases where elements COMU and COMI cannot be used. It should prove particularly advantageous because the weight of the calculation and triangulation of the tangent matrix is, in this case, very important.