# Efficient algorithms for parametric non-linear instability analysis

Antoine Legay, Alain Combescure

► To cite this version:

Antoine Legay, Alain Combescure. Efficient algorithms for parametric non-linear instability analysis. International Journal of Non-Linear Mechanics, 2002, 37 (4-5), pp.709 - 722. 10.1016/S0020-7462(01)00094-4. hal-01422436

HAL Id: hal-01422436

https://hal.science/hal-01422436

Submitted on 25 Dec 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Efficient algorithms for parametric non-linear instability analysis

Antoine Legay, Alain Combescure

LMT-Cachan, E.N.S. de Cachan, C.N.R.S., 61, avenue du Président Wilson, 94235 Cachan Cedex, France

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 fits within the framework of non-linear reliability which requires numerous non-linear mechanical calculations with different sets of parameters. Therefore, we aim for a method which has optimum numerical efficiency.

Keywords: Buckling; Thin structure; Non-linear; Reliability

## 1. 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 affect 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 $\lambda_c$ (Fig. 1). Imperfections on structural parameters are known to lower the critical buckling load from $\lambda_c$ to $\lambda_1$ or $\lambda_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.

Several works have addressed the sensitivity of the non-linear response to geometric defects [1-4].

The study of the sensitivity of the response to the different imperfections can be undertaken in several different ways.

The first 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
<figure><img src="image_2.png" /><figcaption>Fig. 1. Perfect and imperfect equilibrium paths and limit point branches.</figcaption></figure>

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 first 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 modified parameters (see details in Section 2.2.1).

In the second step, in order to calculate a minimum number of points to find 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 modified 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 finite 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 stiffened ring (Section 3.3) under external pressure.

## 2. Description of the method

### 2.1. The Riks' method for the basic curve

The implicit updated Lagrangian algorithm based on Riks' method [9] involves two loops: the first 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 first 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 matrix.

#### 2.1.1. Evolution loop

This part of the algorithm is used to initialize the values for the equilibrium calculation. The external forces are calculated on the configuration of Step n; they are used to obtain a “linear” estimate of the displacement increment for the next step n + 1.

#### 2.1.2. 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 verified are

$$ \Delta U_{n+1}^{i+1} = \underbrace{\Delta \lambda_{n+1}^{i+1} U_{n+1}^{\text{linear}}}_{\text{linear part}} + \underbrace{\Delta U_{\text{NL}}^{i+1}}_{\text{non-linear part}}, \\ \Delta U_{n+1}^{i+1} \Delta U_n = \Delta \beta, \qquad (1) $$

with

$$ \Delta U_{\text{NL}}^{i+1} = \Delta U_{\text{NL}}^{i} + \mathbf{K}_{\mathbf{T}}^{-1} R^{i} \qquad (2) $$

and

$$ R^i = F_{\text{internal}}^i - F_{\text{external}}^i. \qquad (3) $$

In Eq. (2), which increments the non-linear part with the residual, the notation $\mathbf{K}_T^{-1}$ is symbolic: it
means that one actually solves the system by Crout decomposition. System (1) has the exact solution:

$$ \Delta\lambda_{n+1}^{i+1} = \frac{\Delta\lambda_{n+1}^{i}}{2} + \frac{\Delta\beta}{2U_{n+1}^{\text{linear}}\Delta U_{n}} - \frac{\Delta U_{n}^{T}\Delta U_{\text{NL}}^{i+1}}{2U_{n+1}^{\text{linear}}\Delta U_{n}}, $$

$$ \Delta U_{n+1}^{i+1} = \Delta\lambda_{n+1}^{i+1} U_{n+1}^{\text{linear}} + \Delta U_{\text{NL}}^{i+1}. $$

The internal and external forces are calculated in updated Lagrangian mode on the configuration of the previous iteration i.

The arc-length increment is given by the variable $\Delta\beta$. It is initialized using the linear part of the first calculation step.

$$ \Delta\beta = U_{1}^{\text{linear}^{T}} U_{1}^{\text{linear}}. $$

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.

## 2.2. The calculation of one modified curve

### 2.2.1. Calculation of one point on the modified curve

The response of a modified structure is usually little different 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 modified calculations the Crout decompositions of the tangent matrices previously stored during the reference calculation (Fig. 2). Thus, we obtain the modified 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

<figure><img src="image_3.png" /><figcaption>Fig. 2. Choice of the reference matrix to be used.</figcaption></figure>

of the calculation. One seeks the step of the reference calculation whose cumulated arc length value is closest to that of the current modified calculation step.

### 2.2.2. Calculation of a portion of the modified 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 modified calculation at an arbitrary point on the reference curve. For this first point, we calculate the internal and external forces in total Lagrangian mode with respect to the undeformed modified structure (see Fig. 3).

Once the starting point on the reference curve has been chosen, we steer the algorithm differently in order to converge toward a point on the modified 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_r^N$: total displacement of the reference structure at the starting point

* $\Delta U_r^N$: displacement increment of the reference structure at the starting point

* $\lambda_r^N$: loading parameter of the reference structure at the starting point

* $V$: total displacement of the modified structure

* $\Delta V$: displacement increment to get from the reference structure to the modified structure
<figure><img src="image_4.png" /><figcaption>Fig. 3. Relation between the reference and modified calculations.</figcaption></figure>

* $\Delta\lambda$: loading increment to go from the reference structure to the modified structure

We initialize the calculation as follows

* displacement: $V = U_r^N$ and $\Delta V^0 = 0$

* loading: $\lambda = \lambda_r^N$ and $\Delta \lambda^0 = 0$

* tangent matrix: $K_T^{-1} = K_{T,ref}^{-1}$

* geometry: $\Omega_0 = \Omega_{modified}^0$

We seek a displacement correction $\Delta V$ as well as a loading parameter correction $\Delta\lambda$ such that the cumulated arc-length of the new displacement is the

same as that of Step N of the reference calculation:

$$||V|| = ||U_r^N||.$$

Thus, the steering during the first total Lagrangian calculation step can be summarized as

$$\Delta V^{i+1} = \underbrace{\Delta \lambda^{i+1} V_0^{\text{linear}}}_{\text{linear part}} + \underbrace{\Delta V_{\text{NL}}^{i+1}}_{\text{nonlinear part}}, \\ ||U_r^N + \Delta V^{i+1}|| = ||U_r^N||. \qquad (4)$$

The system of equations (4) is a first-degree equation in $\Delta V^{i+1}$ and second-degree equation
in $\Delta\lambda^{i+1}$:

$$ \Delta\lambda^2 V_0^{\text{linear}^2} + 2\Delta\lambda V_0^{\text{linear}} (U_r^N + \Delta V_{\text{NL}}^{i+1}) + 2\Delta V_{\text{NL}}^{i+1} U_r^N + \Delta V_{\text{NL}}^{i+1} = 0. $$

There are two solutions to this second-degree equation, of which we always choose the smaller. The non-linear part of the displacement is defined by

$$ \Delta V_{\text{NL}}^{i+1} = \Delta V_{\text{NL}}^{i} + \mathbf{K}_{\mathbf{T}}^{-1} R^{i}. \quad (5) $$

We need to calculate the internal and external forces on the undeformed configuration of the modified structure $\Omega_m^0$ by the total Lagrangian approach. Then, we can calculate the residual $R^i$

$$ V = U_r^N + \Delta V^i, \quad (6) $$

$$ \lambda = \lambda_r^N + \Delta\lambda^i, $$

$$ R^i = F_{\text{internal}}^i(V, \lambda) - F_{\text{external}}^i(V, \lambda). \quad (7) $$

It ensures convergence of the equilibrium point once the residual $||R^i||$ gets small enough.

We can summarize this in the following algorithm:

initializations:

$$ V = U_r^N $$

$$ \Delta V^0 = 0 $$

$$ \lambda = \lambda_r^N $$

$$ \Delta\lambda^0 = 0 $$

$$ \mathbf{K}_{\mathbf{T}}^{-1} = \mathbf{K}_{\mathbf{T}}^{-1 \text{ ref}} $$

WHILE $||R^i||$ is too big DO

compute $\Delta\lambda^i$

compute $\Delta V^i$

compute the new displacement: $V = U_r^N + \Delta V^i$

compute the new loading: $\lambda = \lambda_r^N + \Delta\lambda^i$

compute $F_{\text{internal}}^i(V, \lambda)$

compute $F_{\text{external}}^i(V, \lambda)$

compute the residual: $R^i = F_{\text{internal}}^i(V, \lambda) - F_{\text{external}}^i(V, \lambda)$

compute the new non-linear part
$\Delta V_{\text{NL}}^{i+1} = \Delta V_{\text{NL}}^{i} + \mathbf{K}_{\mathbf{T}}^{-1} R^i$

END WHILE

<figure><img src="image_5.png" /><figcaption>Fig. 4. Complete calculation of the modified curve.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Calculation of a portion of the modified curve.</figcaption></figure>

### 2.2.3. A general strategy to find the buckling point on the modified curve

The search for the critical buckling point on the modified curve can be performed in several ways.

#### 2.2.3.1. Complete calculation of the modified 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 modified curve.

#### 2.2.3.2. Calculation of a portion of the modified curve.
Since the method proposed in Section 2.2.2 allows us to begin the calculation at an arbitrary point on the modified 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 different starting points on the reference curve (Fig. 6). This gives only an approximation of the maximum, but this can be sufficient
<figure><img src="image_7.png" /><figcaption>Fig. 6. Approximate search for the peak.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Exact search for the peak.</figcaption></figure>

for first reliability calculations where an error on the mechanical calculation is acceptable.

2.2.3.4. 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).

### 2.2.4. 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:

* $u$: displacement

* $\lambda$: loading

* $X$: parameters

$G(u, \lambda, X)$

$$= \begin{cases} F(u, \lambda, X) & \text{equilibrium equation} \\ g(u, \lambda, X) & \text{limit point criterion} \\ N(u, \lambda, X) & \text{steering} \end{cases} = 0.$$

The method in [7] is based on a classical finite element calculation, whereas [8] uses an asymptotic numerical calculation (see [12]).

Both methods are efficient 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.

### 2.2.5. Efficiency considerations

The time saved for the complete calculation of the modified 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 modified 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.

## 3. Application to axisymmetric structures with non-axisymmetric initial imperfections

### 3.1. The COMU and COMI finite elements

These two finite 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 efficient 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 efficient. We will now briefly describe the formulation of the elements. Both elements are conical Mindlin shell elements with two nodes. The displacement field is decomposed into a Fourier series. For the
<figure><img src="image_9.png" /><figcaption>Fig. 8. Ring under external pressure.</figcaption></figure>

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 $\theta_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 efficiency of the element is poor compared to a full 3D shell analysis due to the Fourier decomposition of the displacement field.

### 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, \theta, z)$. The local reference frame at point $M$ is defined by $(\vec{n}, \vec{s}, \vec{t})$, where $\vec{n}$ is the normal to the shell, $\vec{s}$ is the tangent along the meridian direction and $\vec{t}$ 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 $\beta$ around direction $\vec{t}$.

$$ U = \begin{pmatrix} w \\ u \\ v \\ \beta \end{pmatrix}. \qquad (8) $$
<figure><img src="image_10.png" /><figcaption>Fig. 9. Complete modified curve.</figcaption></figure>

Let us designate by $\Phi$ the angle $(-\vec{z}, \vec{s})$. As usual for axisymmetric structures, we decompose the displacement field $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 field in the same frame. For a conical shell with no initial imperfections, the linear strains are, as usual, separated into membrane ($\varepsilon^{lm}$), bending ($\varepsilon^{lb}$) and shear ($\varepsilon^s$) strains given by

$$ \varepsilon^{\text{lm}} = \begin{cases} \varepsilon_{ss}^{\text{lm}} = \frac{\partial u}{\partial s} \\ \varepsilon_{\theta\theta}^{\text{lm}} = \frac{1}{r}(\frac{\partial v}{\partial \theta} + w \cos \Phi + u \sin \Phi) \\ \varepsilon_{s\theta}^{\text{lm}} = \frac{1}{2}[\frac{1}{r}(\frac{\partial u}{\partial \theta} - v \sin \Phi) + \frac{\partial v}{\partial s}] \end{cases}, \quad (9) $$

$$ \varepsilon^{\text{lb}} = \begin{cases} \varepsilon_{ss}^{\text{lb}} = -\frac{\partial^2 w}{\partial s^2} \\ \varepsilon_{\theta\theta}^{\text{lb}} = \frac{1}{r^2} \left( -\frac{\partial^2 w}{\partial \theta^2} - r \sin \Phi \frac{\partial w}{\partial \theta} + \cos \Phi \frac{\partial v}{\partial \theta} \right) \\ \varepsilon_{s\theta}^{\text{lb}} = \frac{1}{2r^2} \left[ -2r \frac{\partial^2 w}{\partial s \partial \theta} + 2r \cos \Phi \frac{\partial v}{\partial s} - v \cos \Phi \sin \Phi + \sin \Phi \frac{\partial w}{\partial \theta} \right] \end{cases}, \quad (10) $$

$$ \varepsilon^s = \frac{1}{2} \left( \beta + \frac{\partial w}{\partial s} \right). \quad (11) $$

The quadratic membrane strains $\varepsilon^q(U, U)$ used in large displacement analysis as well as in buckling analysis are given by

$$ \varepsilon^q = \begin{cases} \varepsilon_{ss}^q = \left(\frac{\partial w}{\partial s}\right)^2 + \left(\frac{\partial u}{\partial s}\right)^2 + \left(\frac{\partial v}{\partial s}\right)^2 \\ \varepsilon_{\theta\theta}^q = \frac{1}{r^2} \left[ \left(\frac{\partial w}{\partial \theta}\right)^2 - v \cos \Phi \right]^2 \\ + \left(\frac{\partial u}{\partial \theta}\right)^2 - v \sin \Phi \right]^2 \\ + \left(\frac{\partial v}{\partial \theta}\right)^2 + w \cos \Phi + u \sin \Phi \right]^2 \\ \varepsilon_{s\theta}^q = \frac{1}{r} \left[ \left(\frac{\partial w}{\partial \theta}\right)^2 - v \cos \Phi \right] \\ \frac{\partial w}{\partial s} + \left(\frac{\partial u}{\partial \theta}\right)^2 - v \sin \Phi \frac{\partial u}{\partial s} + \left(\frac{\partial v}{\partial \theta}\right)^2 + w \cos \Phi \\ + u \sin \Phi \frac{\partial v}{\partial s} \end{cases}. \quad (12) $$

### 3.1.1.3. Strains for the geometrically imperfect shell element.
We now introduce a geometric initial imperfection represented by the displacement
<figure><img src="image_11.png" /><figcaption>Fig. 10. Approximation of the peak of the modified curve.</figcaption></figure>

vector $d_0$ from the perfect structure to the imperfect one. The displacement field can be discretized either by its expression in a Fourier basis for the COMU element or on the set of discrete points $\theta_i$ for the COMI element. The strain field is now

$$ \varepsilon^{\text{lm}}(\underline{d}_0, \underline{U}) = \varepsilon^{\text{lm}}(\underline{U}) + \varepsilon^q(\underline{d}_0, \underline{U}), \quad (13) $$

$$ \varepsilon^{\text{lb}}(\underline{d}_0, \underline{U}) = \varepsilon^{\text{lb}}(\underline{U}), \quad (14) $$

$$ \varepsilon^q(\underline{d}_0, \underline{U}) = \varepsilon^q(\underline{U}, \underline{U}). \quad (15) $$

In this expression, $\varepsilon^q(\underline{d}_0, \underline{U})$ is obtained simply by replacing one of the vectors $\underline{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

$$ \sigma^s = \frac{5}{6} G \varepsilon^s, \quad (16) $$

where $G$ is the shear modulus. In the elastic–plastic case, an Illushin “global” model is chosen [15].

### 3.1.2. The COMU and COMI finite elements

#### 3.1.2.1. Geometry.
The element is a two-node conical element. The discretization of the geometrical imperfection differs between the two finite elements: for the COMU element, the imperfection is given as coefficients 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 $\theta_i$. For the COMI element, the thickness field is also defined on the set of points $\theta_i$. In the case of large displacement analysis, the geometry is updated at each load increment.

#### 3.1.2.2. Degrees of freedom.
At each node, the displacement field $\underline{U}$ has four components for each Fourier mode chosen. It is expressed in the cylindrical reference frame.

#### 3.1.2.3. Shape functions.
The shape functions are linear for each of the four degrees of freedom. The shear energy is used to connect the rotation
<figure><img src="image_12.png" /><figcaption>Fig. 11. Stiffened ring: geometry.</figcaption></figure>

with the derivative of the normal displacement w
with respect to s.

3.1.2.4. *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 (co-ordinate θ) is analytical for the *COMU* element and numerical for the imperfection field and for the *COMI* element.

3.1.3. *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.

## 3.2. Example of the parametric non-linear buckling of a ring under external pressure

3.2.1. *Definition of the problem*

The characteristics of the ring are (Fig. 8):

• radius $R = 100$ mm,

• thickness $h = 1$ mm,
<figure><img src="image_13.png" /><figcaption>Fig. 12. Stiffened ring: material behavior.</figcaption></figure>

* buckling mode studied $n=2$,

* amplitude of the defect on mode $n$: $\delta = 0.1$ mm,

* Young's modulus $E = 200 000$ MPa,

* Poisson's ratio $\nu = 0.28$,

* yield stress 500 MPa,

* Euler's critical pressure: $h^3 E(n^2 - 1)/12vR^3$.

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 modified 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 significantly compared to the configuration without hinge. Fu and Waas [18] studied a ring made of a composite material. More recently, Hudleston and Sivaselvan [19] analyzed the buckling of a compressible ring.

### 3.2.2. Complete modified curve

Fig. 9 shows the result obtained by calculating the modified curve completely with the method in Section 2.2.1.

### 3.2.3. Part of the modified curve

We applied the method in Section 2.2.2. Only nine points suffice to find the limit point of the modified curve (Fig. 10).

## 3.3. Buckling of a stiffened cylinder

### 3.3.1. Definition of the problem

The results presented here were calculated on a stiffened 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 defined on the second harmonic, its amplitude is 9.952 mm. The geometry and the
<figure><img src="image_14.png" /><figcaption>Fig. 13. Stiffened ring: complete modified curve.</figcaption></figure>

elastic-plastic constitutive law of the material are
defined on Figs. 11 and 12. The young's modulus
is $E = 200\,000$ MPa, the conventional yield stress
is $\sigma_y = 160$ MPa.

3.3.2. *Complete modified curve*

The perturbation is introduced on the exterior
thickness: -17%. Fig. 13 shows, on the example of
the stiffened ring, the shape of the total radial dis-
placement of the shell node.

3.3.3. *Portion of the modified curve*

Fig. 14 shows the shape of the total radial dis-
placement after a total Lagrangian calculation of
the first step (see Section 2.2.2).

3.3.4. *Quality of the response*

As can be seen from the curves, the results of the
calculation with the efficient method are identical
to those of the direct calculation of the modified

curve. One can thus appreciate the quality of the
method when applied to this example.

3.3.5. *Efficiency*

On this example, there is saving of time on the
order of 40% when the modified 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.

**4. Conclusion and perspectives**

We have introduced an original method of para-
metric non-linear calculation in order to provide,
with reliability analysis in mind, a capability to per-
form successive calculations for different, but sim-
ilar, sets of parameters.
<figure><img src="image_15.png" /><figcaption>Fig. 14. Stiffened ring: approximation of the peak of the modified curve.</figcaption></figure>

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 first step of the calculation of the modified curve.

These ideas are inapplicable if the responses of the reference structure and of the modified structure are very different, particularly in cases where the primary buckling modes have different shape. If so, one can simply change reference calculation.

For plastic buckling, the most efficient 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 modified 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 circumfer-

ference (COMI). By using these elements one can achieve a further gain in efficiency 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 finite element code with three-dimensional shell elements. This evaluation is in progress.

The efficient 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.

## References

[1] G. Cederbaum, J. Arbocz, On the reliability of imperfection-sensitive long isotropic cylindrical shells, Struct. Safety 18 (1) (1996) 1–9.
[2] J. Arbocz, J.M.A.M. Hol, Collapse of axially compressed cylindrical shells with random imperfections, AIAA J. 29 (12) (1991) 2247-2256.

[3] G. Cederbaum, J. Arbocz, Reliability of shells via Koiter formulas, Thin-Walled Struct. 24 (2) (1996) 173-187.

[4] L.A. Godoy, E.O. Taraco, Design sensitivity of post-buckling states including material constraints, Comp. Meth. Appl. Mech. Eng. 188 (2000) 665-679.

[5] M. Lemaire, A. Mohamed, Finite element and reliability: a happy marriage? in: Nineth IFIP WG 7.5 Working Conference, Reliability and optimization of structural systems, September 25-27, 2000, Ann Arbor, Michigan, USA, 2000.

[6] T. Abdo, R. Rackwitz, A new beta-point algorithm for large time-invariant and time variant reliability problems, Reliability and optimization of structures, 3rd WG 7.5 IFIP Conference, Vol. 1, Berkeley, 1990, p. 11.

[7] A. Eriksson, C. Pacoste, A. Zdunek, Numerical analysis of complex instability behaviour using incremental-iterative strategies, Comp. Meth. Appl. Mech. Eng. 179 (1999) 265-305.

[8] B. Cochelin, A path following technique via an asymptotic numerical method, Comput. Struct. 29 (1994) 1181-1192.

[9] E. Riks, An incremental approach to the solution of snapping and buckling problems, Int. J. Solids Struct. 15 (1979) 524-551.

[10] A. Combescure, Strategies for non-linear static calculations, Irdi/dedr/demt/smts/lams, CEA.

[11] M.A. Crisfield, A fast incremental/iterative solution procedure that handles "snap-through", Comput. Struct. 13 (1980) 55-62.

[12] B. Cochelin, N. Damil, M. Potier-Ferry, The asymptotic-numerical-method: an efficient perturbation technique for non-linear structural mechanics, Rev. Eur. Éléments Finis 3 (1994) 281-297.

[13] G. Gusic, A. Combescure, J.F. Jullien, Influence of circumferential thickness variations on buckling of cylindrical shells under external pressure, Comput. Struct. 74 (4) (2000) 461-477.

[14] G.D. Galletly, A. Combescure, Plastic buckling of complete toroidal shells of elliptical cross-section subjected to internal pressure, Thin-Walled Struct. 34 (2) (1999) 135-146.

[15] M.A. Crisfield, Nonlinear Finite Element Analysis of Solids and Structures, Vol. I, Wiley, New York, 1988.

[16] E. Riks, Some computational aspects of the stability of non-linear structures, Comp. Meth. Applied Mech. Eng. 47 (1984) 219-259.

[17] C.Y. Wang, Postbuckling of an externally pressurized ring with a hinge, Int. J. Solids Struct. 27 (10) (1991) 1287-1293.

[18] L. Fu, A.M. Waas, Initial post-buckling behavior of thick rings under uniform external hydrostatic pressure, J. Appl. Mech. 62 (1995) 338-345.

[19] J.V. Huddleston, M.V. Sivaselvan, Buckling and postbuckling of compressible circular rings under hydrostatic pressure, AIAA J. 38 (12) (2000) 2361-2363.

[20] J.M. Bourinet, N. Gayton, M. Lemaire, A. Combescure, Reliability analysis of stability of shells based on combined finite element and response surface methods, in: M. Padadrakis, A. Samartin, E. Onate (Eds.), Computational Methods for Shell and Spatial Structures, IASS-IACM 2000, ISASTR-NTUA, Athens, Greece, 2000.