# On the bifurcation analysis of thin multilayer structures by asymptotic numerical method

Hamza Azzayani, Hamid Zahrouni, Norman Mathieu, Pascal Ventura, Michael Brun, Michel Potier-Ferry

► To cite this version:

Hamza Azzayani, Hamid Zahrouni, Norman Mathieu, Pascal Ventura, Michael Brun, et al. On the bifurcation analysis of thin multilayer structures by asymptotic numerical method. 2023. hal-03817307v3

HAL Id: hal-03817307

https://hal.univ-lorraine.fr/hal-03817307v3

Preprint submitted on 28 Apr 2023

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
# On the bifurcation analysis of thin multilayer structures by asymptotic numerical method

Hamza Azzayani, Hamid Zahrouni, Norman Mathieu, Pascal Ventura, Michael Brun, Michel Potier-Ferry

Université de Lorraine, CNRS, Arts et Métiers ParisTech, LEM3, F-57000 Metz, France
Laboratory of Excellence on
Design of Alloy Metals for low-mAss Structures (DAMAS), Université de Lorraine, France
{hamza.azzayani,hamid.zahrouni,norman.mathieu,pascal.ventura,michael.brun,michel.potier-ferry}@univ-lorraine.fr

## Abstract

In this paper, we are interested in the stability analysis of multilayer thin shells using asymptotic numerical method (ANM) associated to Padé approximants. This technique is very efficient in solving nonlinear problems in particular for instability modeling of thin structures thanks to the high order algorithm leading to a high accuracy in computing singular points along the nonlinear solution branches. We present different techniques to detect bifurcation points. The first technique is based on a bifurcation indicator introduced in the nonlinear problem in the form of a scalar function representing the intensity of a fictitious perturbation force which is evaluated along the equilibrium branch and which vanishes exactly at singular points. The second technique is based on Padé approximants that can be used as a bifurcation indicator by analyzing the denominator of rational fractions. The bifurcation corresponds to the first real root of the denominator. A third technique consists in a combination of buckling and linear vibrations which allows the extraction of singular points by analyzing the evolution of natural frequencies along the equilibrium path. Several numerical examples show the efficiency and robustness of the proposed methods.

**Keywords:** Asymptotic numerical method, instability, bifurcation, Padé approximants, thin multilayer structures
# 1. Introduction

Numerical modeling of structural problems has become an indispensable tool since the development of the finite element method, allowing a better understanding of the phenomena that occur during deformation. For reasons of cost and competitiveness, manufacturers tend to use strong and light materials and complex structures. Multilayer composite shell structures have been widely employed in the fields of civil engineering, aerospace, automotive, energy and many others (for example: thermal protection of satellite launch tanks [1], printed circuit boards [2, 3], morphing airfoils [4],...). However, these structures are very sensitive to instabilities (e.g. buckling) [5, 6], and are subjected to a drastic reduction of the critical load in the presence of geometrical or material defects. This instability often presents complex mechanical behaviors with large deformation, strong nonlinearity, multiple bifurcation points and multiple equilibrium paths. In the context of multilayered materials, the instability is very sensitive to the geometry, the number of layers in the thickness, their orientations and the anisotropic material properties [7, 8]. In order to have a better understanding of the complex instability phenomena in various materials, numerous experimental studies as well as analytical and numerical models have been developed [9, 10, 11]. Another application of these numerical models is the design of stacking sequence to obtain multi-stable structures, i.e. structures having several stable states for the same loads. The corresponding load-displacement curves are highly nonlinear with many bifurcations and stability changes, see [12]. In this paper, various numerical procedures are developed to predict these changes of stability for multi-layer shells.

Nonlinear problems in structural mechanics are usually solved by iterative methods based on the predictor-corrector process. The principle is to follow the nonlinear solution branch of the problem in a stepwise manner point by point in two steps. The first step consists in linearizing the initial nonlinear problem and predicting a solution, and the second step consists in correcting the residual by successive iterations to reach the equilibrium. This method is widely used in finite element codes. However, the computation time remains a major obstacle of this method. For instability problems, these algorithms require an adequate path-following strategy to follow complex response branches. In parallel to these methods, a second approach to solve nonlinear problems is the perturbation technique. In this case, the solution branch is determined as a power series with respect to a parameter $a$
40

which can be defined in the same way as the control parameters for classical iterative algorithms. This gives an approximate analytical representation of the solution. This method is based on an analytical development and the applications were limited to particular problems with simple geometry.

45

The asymptotic numerical method (ANM) based on the perturbation technique presents an alternative [13]. It allows to solve nonlinear problems thanks to the combination of asymptotic developments and numerical methods. The unknowns of the problem are developed in power series. The ANM transforms the initial problem (nonlinear) into a succession of linear problems admitting the same tangent stiffness operator. Then, these problems are solved by a numerical method such as the finite element method (FEM). The ANM allows to compute a large part of the nonlinear branch with only one decomposition of the stiffness matrix per step [14, 15]. Many studies have been presented in the literature to compare the number of matrix decompositions required for a given solution path by the asymptotic numerical method or by more classical approaches [15]. However, the ANM has been the most efficient compared to the Newton-Raphson methods. As the power series have a limited radius of convergence, a continuation technique has been proposed to obtain the whole solution branch. Each end of step will serve as the starting point of the next step and the step length is computed a posteriori, requiring that the relative difference between solutions at two consecutive orders is less than a user-defined accuracy parameter. The automatic determination of the step lengths from the terms of the series and according to the local nonlinearity of the response curve is the key point of the efficiency and robustness of these methods compared to classical iterative methods. ANM has been successfully applied in many fields: Plastic beams [16, 17], Navier-Stockes equation [18], contact mechanics [19], instability of fiber reinforced composites [20],...

50

55

60

65

The use of a rational representation (Padé approximants [21]) instead of a polynomial representation improves the validity range of the solution [19]. Najah et al. [22] have shown the efficiency of the Padé approximants over the power series representation and the Rayleigh-Ritz technique [22, 23]. Using rational fractions that share a common denominator can help to minimize the number of poles in a given representation. This can be particularly useful in the context of numerical methods, where a large number of poles can lead to computational inefficiencies.

70

Within the framework of ANM, bifurcation indicators have been proposed to identify critical points. A first indicator is a scalar function computed by

75
introducing a fictitious perturbation force in the equilibrium problem [24]. The indicator can be determined explicitly along the equilibrium branch by the perturbation technique. The roots of this function characterize the singular points [25]. Several applications of this technique have been implemented in solid and fluid mechanics [18, 26]. The second indicator consists in analyzing the poles of the rational representation. Overall, the bifurcation point corresponds to the smallest real pole of the Padé approximant [24]. A third technique for detecting singular points is a combination of buckling and vibration analysis. The method starts with the calculation of the equilibrium branch. Then, we proceed to the determination of the eigenfrequencies with the corresponding vibration modes around this branch. In this method, the natural frequencies are used as indicator of bifurcation. Several studies show that the natural frequencies decrease when the applied load increases and their roots correspond to critical points (bifurcation points or limit points) [27].

In what follows, we present in section 2 the shell formulation for equilibrium conditions and the asymptotic numerical algorithm to solve the resulting nonlinear problem. In section 3, we give details about the three techniques used to detect bifurcation points. In section 4, we present some numerical results and in section 5 we conclude our study.

## 2. Static shell formulation and resolution algorithm

### 2.1. Geometric and kinematic description

In this paper, we consider perfectly bonded multilayer thin structures. The behavior of each layer is described as orthotropic elastic. While this limits the range of practical applications to small deformations, it does allow for large rotations to occur. A shell formulation, adapted for large displacements and large rotations, using only displacement variables is used [28]. Figure 1 shows the geometry and kinematics of the proposed shell in the reference and deformed configuration. The position vector of an arbitrary material point of the initial configuration can be defined as follows:

$$x(\theta_1, \theta_2, \theta_3) = r(\theta_1, \theta_2) + \theta_3 a_3(\theta_1, \theta_2) \quad (1)$$

where $r$ is the projection of this point on the mid-surface, $a_3$ the shell director vector and $(\theta_1, \theta_2, \theta_3)$ the curvilinear coordinates. Then, the covariant
<figure><img src="image_2.png" /><figcaption>Figure 1: Geometric and kinematic description of the shell</figcaption></figure>

base vectors in the undeformed configuration are expressed as:

$$ \left\{ \begin{array}{l} g_{\alpha} = \frac{\partial x}{\partial \theta_{\alpha}} = r_{,\alpha} + \theta_{3} a_{3},_{\alpha} \\ g_{3} = \frac{\partial x}{\partial \theta_{3}} = a_{3} \end{array} \right. \qquad \alpha = 1, 2 \qquad (2) $$

where $,_{\alpha} = \frac{\partial}{\partial \theta_{\alpha}}$

Assuming a linear variation of the displacement in the thickness direction, the displacement field corresponding to the deformed configuration can be written as:

$$ u(\theta_1, \theta_2, \theta_3) = v(\theta_1, \theta_2) + \theta_3 w(\theta_1, \theta_2) \qquad (3) $$

where $v$ and $w$ represent, respectively, the mid-surface displacement and the difference vector between the undeformed and deformed shell director vectors. In this formulation, six degrees of freedom can be distinguished, ($v_1, v_2, v_3$) are relative to the translation of the mid-surface and ($w_1, w_2, w_3$) updating the director vector. Thus, we can define the position vector of the deformed shell by the following expression:

$$ \bar{x} = (r + v) + \theta_3(a_3 + w) = \bar{r} + \theta_3 \bar{a}_3 \qquad (4) $$

In the same way as Eq.(2), the covariant base vectors in the deformed configuration are expressed as:
$$ \left\{ \begin{aligned} \bar{g}_{\alpha} &= \frac{\partial \bar{x}}{\partial \theta_{\alpha}} = \bar{r}_{,\alpha} + \theta_{3} \bar{a}_{3,\alpha} & \alpha = 1, 2 \\ \bar{g}_{3} &= \frac{\partial \bar{x}}{\partial \theta_{3}} = \bar{a}_{3} \end{aligned} \right. \qquad (5) $$

The reader wishing more details on this formulation can refer to the work of Büchter et al. [28].

## 2.2. Strain and stress equations

The Green-Lagrange strain tensor is written in the contravariant basis:

$$ \gamma = \frac{1}{2}(\bar{g}_{ij} - g_{ij})g^i \otimes g^j \qquad g_{ij} = g_i \cdot g_j \qquad (6) $$

where $g_{ij}$ and $\bar{g}_{ij}$ are respectively the covariant components of the metric tensors in the reference and deformed configuration and $g^i$ denote the contravariant basis vectors. it can be defined by the orthogonality condition $g^i \cdot g_j = \delta^i_j$.

From Eqs (2) and (5), the components of the Green-Lagrange strain can be written as:

$$ \gamma_{ij} = \alpha_{ij} + \beta_{ij}\theta_{3} + Q_{ij}(\theta_{3})^{2} \qquad (7) $$

with

$$ \left\{ \begin{aligned} \alpha_{ij} &= \frac{1}{2}(\bar{a}_i \bar{a}_j - a_i a_j) \\ \beta_{\alpha\beta} &= \frac{1}{2}(\bar{a}_{3,\alpha} \bar{a}_{\beta} + \bar{a}_{3,\beta} \bar{a}_{\alpha} - a_{3,\alpha} a_{\beta} - a_{3,\beta} a_{\alpha}) \\ \beta_{\alpha 3} &= \frac{1}{2}(\bar{a}_{3,\alpha} \bar{a}_3 - a_{3,\alpha} a_3) \\ \beta_{33} &= 0 \\ Q_{\alpha\beta} &= \frac{1}{2}(\bar{a}_{3,\alpha} \bar{a}_{3,\beta} - a_{3,\alpha} a_{3,\beta}) \\ Q_{\alpha 3} &= 0 \\ Q_{33} &= 0 \end{aligned} \right. \qquad (8) $$

To avoid numerical locking drawback in shell modeling, the EAS (Enhanced Assumed Strain) concept proposed by Simo and Rifai [29] will be used. This method consists in introducing an additional strain field $\tilde{\gamma}$ incompatible with the displacement and which is chosen, by construction, orthogonal to the stress field $S^m$ relative to each layer $m$. This additional variable allows to use a complete 3D constitutive law without condensation. Since no inter-element continuity is required, the additional deformation is
eliminated at the elementary level, thus preserving the formal structure of a 6-parameter shell theory. The incompatible strain field $\tilde{\gamma}$ is written:

$$ \tilde{\gamma} = \theta_3 \tilde{\beta}_{33} g^3 \otimes g^3 \qquad \int_{\Omega} ( {}^t S^m : \tilde{\gamma} ) d\Omega = 0 \qquad (9) $$

The term $\tilde{\beta}_{33}$ provides a linear variation of the strain field through the thickness. Thus, the total strain field $\gamma$ and the stress field (second Piola-Kirchhoff stress tensor) in each layer $S^m$ are written:

$$ \gamma = \gamma(u) + \tilde{\gamma} = (\gamma_l(u) + \gamma_{nl}(u, u)) + \tilde{\gamma} \qquad (10) $$

$$ S^m = D^m : \gamma \qquad (11) $$

in which $\gamma_l(u)$ and $\gamma_{nl}(u, u)$ denote the linear and nonlinear parts of the compatible strain field respectively and $D^m = D^{ijkl(m)}g_i \otimes g_j \otimes g_k \otimes g_l$ refers to the fourth-order elasticity tensor of the $m$th layer. Since the elasticity tensor is often obtained in cartesian coordinates, a transformation to curvilinear coordinates is used in the numerical implementation

$$ D^{ijkl(m)} = C^{pqrs(m)}(e_p.g^i)(e_q.g^j)(e_r.g^k)(e_s.g^l) \qquad (12) $$

where $C^{pqrs(m)}$ presents the components of elasticity tensor in cartesian coordinates.

105

## 2.3. A three-field variational formulation

Based on the Hu-Washizu functional, we consider a shell formulation with three independent fields $(u, S^m, \tilde{\gamma})$ [15]. It is written in the following form:

$$ \pi_{EAS}(u, \tilde{\gamma}, S^m) = \sum_{m=1}^{N} \int_{\Omega} \{ {}^t S^m : \gamma - \frac{1}{2} {}^t S^m : (D^m)^{-1} : S^m \} d\Omega - \lambda P_e(u) \qquad (13) $$

where $N$ refers to the number of layers of laminated shells and $\lambda P_e(u)$ is the work of the external forces ($\lambda$ is a scalar load parameter).

According to the principle of minimum potential energy, the stationarity condition leads to the virtual work equation (Eq.(14)):

$$ \sum_{m=1}^{N} \int_{\Omega} \{ {}^t \delta S : [(\gamma(u) + \tilde{\gamma}) - (D^m)^{-1} : S^m] + {}^t S^m : \delta \gamma(u) \} d\Omega - \lambda P_e(\delta u) = 0 \qquad (14) $$
By defining the mixed unknown vector $U = (u, \tilde{\gamma}, S^m)$, we can write Eq.(14) as the following simple form:

$$R(U, \lambda) = L(U) + Q(U, U) - \lambda F = 0 \quad (15)$$

where $L(.)$ is a linear operator, $Q(., .)$ a quadratic one, $F$ the external load vector and $R$ the residual vector. The expressions of these operators are given by:

$$\langle L(U), \delta U \rangle = \sum_{m=1}^{N} \int_{\Omega} ({}^{t}\delta S : [(\gamma_l(u) + \tilde{\gamma}) - (D^m)^{-1} : S^m] + {}^{t}S^m : \gamma_l(\delta u)) d\Omega \quad (16)$$

$$\langle Q(U, U), \delta U \rangle = \sum_{m=1}^{N} \int_{\Omega} ({}^{t}\delta S : \gamma_{nl}(u, u) + {}^{t}S^m : 2\gamma_{nl}(u, \delta u)) d\Omega \quad (17)$$

$$\langle F, \delta U \rangle = P_e(\delta u) \quad (18)$$

Thus, the problem is to find the mixed variable $U$ and the loading parameter $\lambda$ verifying Eq.(15).

## 2.4. Asymptotic numerical algorithm

We propose to solve the nonlinear problem (15) by the asymptotic numerical method [28]. First, we use the perturbation technique to transform the nonlinear problem into a succession of linear ones all admitting the same tangent operator. The resulting problems up to order $p$ are discretized and solved by the finite element method. Thus, we analytically obtain a part of the solution path with only one stiffness matrix decomposition. Compared to the predictor-corrector method such as Riks method, this method is automatic and more efficient [22].

Starting from an initially known solution $(U_0, \lambda_0)$, the perturbation technique consists in searching the solution path of the nonlinear problem under an asymptotic expansion form with respect to a path parameter $(a)$:

$$\left\{ \begin{array}{l} U(a) = U_0 + aU_1 + a^2U_2 + \dots + a^nU_n \\ \lambda(a) = \lambda_0 + a\lambda_1 + a^2\lambda_2 + \dots + a^n\lambda_n \end{array} \right. \quad (19)$$
in which $U_p$ are mixed unknown vectors, $\lambda_p$ the unknown coefficients at order $p$ and $n$ is the truncation order of the series. The range of validity of the representation (19) depends on the path parameter $a$. By analogy with arc-length iterative methods, the control parameter $a$ is chosen as the projection of the displacement increment $u - u_0$ and the load increment $\lambda - \lambda_0$ on the tangent direction $(u_1, \lambda_1)$ where $<.,>$ is the Euclidian scalar product:

$$a = <u - u_0, u_1> + (\lambda - \lambda_0)\lambda_1 \quad (20)$$

By substituting the series (19) into Eqs (15) and (20) and identifying the terms according to the power of $a$, we obtain a recurrent sequence of linear mixed problems admitting the same tangent operator.

Order 1:

$$\left\{ \begin{array}{l} L_t^0(U_1) = \lambda_1 F \\ <u_1, u_1> + \lambda_1^2 = 1 \end{array} \right. \qquad (21)$$

130

Order 2:

$$\left\{ \begin{array}{l} L_t^0(U_2) = \lambda_2 F - Q(U_1, U_1) \\ <u_2, u_1> + \lambda_2 \lambda_1 = 0 \end{array} \right. \qquad (22)$$

Order $p$:

$$\left\{ \begin{array}{l} L_t^0(U_p) = \lambda_p F - \sum_{r=1}^{p-1} Q(U_r, U_{p-r}) \\ <u_p, u_1> + \lambda_p \lambda_1 = 0 \end{array} \right. \qquad (23)$$

The tangent operator $L_t^0$ depends only on the initial solution and is defined by:

$$L_t^0(\cdot) = L(\cdot) + Q(U_0, \cdot) + Q(\cdot, U_0)$$

135

After stress condensation, the linear problems above are solved using the finite element method. For discretization, we use the eight node shell element with reduced integration detailed in [15]. The finite element discretisation is described in Appendix. The discretized form of the problem at order $p$ ($p \ge 2$) is given by:
$$ \left\{ \begin{array}{l} [K_t^0]\{q_p\} = \lambda_p\{F\} + \{F_p^{nl}\} \\ {}^t\{q_p\} \cdot \{q_1\} + \lambda_p\lambda_1 = 0 \end{array} \right. \qquad (24) $$

140

where $[K_t^0]$ is the classical tangent stiffness matrix at starting point $(U_0, \lambda_0)$, $\{F\}$ is the external force vector, $\{q_p\}$ and $\lambda_p$ are respectively the discretized form of the displacement $u_p$ and the loading parameter at order $p$. The second member $\{F_p^{nl}\}$ depends on the solutions computed at the previous $(p-1)$ orders.

In the ANM, the variables are expressed as power series expansions. However, these series have a radius of convergence that depends on the problem to be solved, the truncated orders of the series and the required accuracy. A simple continuation procedure consists in considering that the relative difference between the displacement series at two successive orders must remain small with respect to a critical value $\delta_1$ [14]. Thus, the maximal value of the path parameter $a_{max}$ is given by:

$$ a_{max} = \left( \delta_1 \frac{\| u_1 \|}{\| u_n \|} \right)^{\frac{1}{n-1}} \qquad (25) $$

145

150

155

with $\delta_1$ a user-defined accuracy parameter and $\| \cdot \|$ the Euclidean norm of the vector. The solution branch is obtained step by step by considering each end of step as the starting point of the next step. Note that this criterion gives a good order of magnitude of the validity of the solution, while it requires almost no computational time. One of the originalities of this procedure is that the step length is determined a posteriori according to the intrinsic characteristics of the computed series. This criterion avoids the user intervention to define the step length (which is often the case in commercial finite element software). Thus, this continuation algorithm is automatic, more robust and easier to use for predicting complex shell behaviors. We would also like to mention that the step length is automatically adaptive to the non-linearity of the problem.

## 2.5. Convergence improvement by Padé approximants

In the previous section, we presented a continuation method based on the polynomial representation. The solution path is represented by truncated power series. Several tests have shown the efficiency and reliability of this algorithm (thin shells in Lagrangian formulation [14, 30], viscous fluid

160
mechanics and plastic structures [15]). In addition, the use of the rational representation, called Padé approximants, instead of the polynomial representation allows to increase the range of validity of the solution. Thus, we present a new algorithm for the numerical computation of the solution path. For a more detailed study on Padé approximants, we refer the reader to the reference [21].

We denote $P_n$ the representation by rational fractions of the displacement $\mathbf{u}$ and the loading parameter $\lambda$ for a truncation order $n$:

$$ \left\{ \begin{array}{l} P_n(\mathbf{u}(\mathbf{a})) = \mathbf{u}_0 + a \frac{D_{n-2}}{D_{n-1}} \mathbf{u}_1 + a^2 \frac{D_{n-3}}{D_{n-1}} \mathbf{u}_2 + \dots + a^{n-1} \frac{1}{D_{n-1}} \mathbf{u}_{n-1} \\ P_n(\lambda(\mathbf{a})) = \lambda_0 + a \frac{D_{n-2}}{D_{n-1}} \lambda_1 + a^2 \frac{D_{n-3}}{D_{n-1}} \lambda_2 + \dots + a^{n-1} \frac{1}{D_{n-1}} \lambda_{n-1} \end{array} \right. \qquad (26) $$

where $D_i(\mathbf{a})$ are polynomials of degree $i$ with real coefficients $d_i$:

$$ D_i(\mathbf{a}) = 1 + a d_1 + a^2 d_2 + \dots + a^i d_i \qquad (27) $$

These rational approximants have a common denominator to limit the number of poles of the representation. Using the same criterion of the polynomial representation based on displacement, we can define a range of validity ($a_{max.p}$) of the rational representation (25). In this case, we must require that the difference between two rational solutions at consecutive orders remains small at the end of the step. This can be expressed by the following relation with $\delta_2$ a new accuracy parameter:

$$ \frac{\| P_n(\mathbf{u}(\mathbf{a}_{max.p})) - P_{n-1}(\mathbf{u}(\mathbf{a}_{max.p})) \|}{\| P_n(\mathbf{u}(\mathbf{a}_{max.p})) - P_0 \|} = \delta_2 \qquad (28) $$

The range of validity ($a_{max.p}$) is searched in the interval $[a_{max}; \beta a_{max}]$ by the bisection method using the criterion (28) where $\beta$ is a parameter defined by the user ($\beta > 1$). The Padé approximants are computed by a simple Gram-Schmidt orthogonalization which does not require too much computational time compared to the representation by series. Several examples show the efficiency and robustness of the Padé algorithm [19].

# 3. Bifurcation detection methods

## 3.1. Bifurcation indicator

The goal of this section is to develop a reliable algorithm, based on ANM, for the determination of bifurcation points. These points are detected by
evaluating, along the equilibrium branch, a scalar function called bifurcation indicator. This scalar is obtained by introducing a fictitious perturbation force in the equilibrium problem, and which vanishes at the singular points. Several examples are presented to evaluate the efficiency of the proposed method [24].

Let $\Delta\mu\mathbf{f}$ be a fictitious perturbation force applied to the structure in a deformed state $(\mathbf{U}, \lambda)$. $\Delta\mu$ represents the intensity of the force and $\Delta\mathbf{U} = (\Delta\mathbf{u}, \Delta\tilde{\mathbf{\gamma}}, \Delta\mathbf{S})$ its associated response. Under these conditions, the perturbed equilibrium equation is written:

$$L(\mathbf{U} + \Delta\mathbf{U}) + Q(\mathbf{U} + \Delta\mathbf{U}, \mathbf{U} + \Delta\mathbf{U}) = \lambda F + \Delta\mu\mathbf{f} \quad (29)$$

Considering the equilibrium equation (Eq.(15)) and neglecting the quadratic terms, the perturbed problem is written:

$$L_t(\Delta\mathbf{U}) = \Delta\mu\mathbf{f} \quad (30)$$

where $L_t(\cdot) = L(\cdot) + 2Q(\mathbf{U}, \cdot)$ is the tangent operator taken at the equilibrium point $(\mathbf{U}, \lambda)$. An additional condition based on the displacement is imposed:

$$\langle L_t^0(\Delta\mathbf{U} - \Delta\mathbf{U}_0), \Delta\mathbf{U}_0 \rangle = 0 \quad (31)$$

where $L_t^0$ is the tangent operator at the starting point $(\mathbf{U}_0, \lambda_0)$.

Eqs (30) and (31) will be solved by the asymptotic numerical method and the unknowns $\Delta\mathbf{U}$ and $\Delta\mu$ are written as power series:

$$\left\{
\begin{aligned}
\Delta\mathbf{U}(\mathbf{a}) &= \Delta\mathbf{U}_0 + a\Delta\mathbf{U}_1 + a^2\Delta\mathbf{U}_2 + \dots + a^n\Delta\mathbf{U}_n \\
\Delta\mu(\mathbf{a}) &= \Delta\mu_0 + a\Delta\mu_1 + a^2\Delta\mu_2 + \dots + a^n\Delta\mu_n
\end{aligned}
\right.
\quad (32)$$

In the same way as the fundamental branch and by substituting Eq.(32) in Eqs (30) and (31), we obtain a sequence of linear problems

Order 0:

$$L_t^0(\Delta\mathbf{U}_0) = \Delta\mu_0\mathbf{f} \quad (33)$$

$\Delta\mathbf{U}_0$ is computed by imposing $\Delta\mu_0 = 1$.

Order $p \ge 1$:

$$\left\{
\begin{aligned}
L_t^0(\Delta\mathbf{U}_p) &= \Delta\mu_p\mathbf{f} - \sum_{r=1}^p (Q(\mathbf{U}_r, \Delta\mathbf{U}_{p-r}) + Q(\Delta\mathbf{U}_{p-r}, \mathbf{U}_r)) \\
\langle \Delta\mathbf{U}_p, \mathbf{f} \rangle &= 0
\end{aligned}
\right.
\quad (34)$$
The $U_r$ vectors correspond exactly to those determined during the computation of the equilibrium branch. The application of the ANM requires the computation of a second series at each step. Nevertheless, the corresponding computation time remains less important, because the stiffness matrix used for the calculation of the equilibrium branch is the same for the bifurcation indicator.

Such as the equilibrium branch, the expansion terms of the series are determined by the finite element method after a condensation step. The discretization of the problem at order $p$ gives:

$$ \left\{ \begin{array}{l} [K_t^0]\{\Delta q_p\} = \Delta\mu_p f + \{\Delta F_p\} \\ t\{\Delta q_p\}[K_t]\{\Delta q_0\} = 0 \end{array} \right. \qquad (35) $$

where $[K_t^0]$ denotes the tangent stiffness matrix at the starting point, $\{\Delta q_0\}$ and $\{\Delta q_p\}$ are respectively the nodal displacement vectors at orders 0 and $p$ associated to the perturbation force $\{f\}$. The vector $\{\Delta F_p\}$ depends on the solutions $U$ up to order $p$ and $\Delta U$ up to order $(p-1)$. Thus, the system (35) gives:

$$ \Delta\mu_p = - \frac{\langle \Delta F_p, \Delta q_0 \rangle}{\langle f, \Delta q_0 \rangle} \qquad (36) $$

Now that we have computed the terms $\Delta\mu_p$ ($1 \le p \le n$), we can construct the series development of $\Delta\mu$ and study its sign along the equilibrium branch in order to detect the bifurcation.

## 3.2. Poles of the Padé approximants

Another simple method to detect bifurcation points is the analysis of the poles of the Padé approximants. In fact, the bifurcation points correspond to the real roots of the denominator of the rational fraction [24][25]. We note here another interest of the rational representation with respect to the polynomial representation.

## 3.3. Natural frequency of vibration

Buckling and vibration are instability phenomena that can coexist and induce large displacements generally causing damage in structural mechanics. The knowledge of the natural frequencies and the critical load with their corresponding mode shapes allows to predict the resonance and instability regions. A classical technique to combine vibration and buckling analysis is to start first with the computation of the equilibrium branch. Afterwards,
one proceeds to the determination of the frequencies with the corresponding vibration modes around this branch. Several studies show that the natural frequencies decrease as the applied load increases and exactly cancel at the critical points (bifurcation points or limit points) [27]. We use the asymptotic numerical method to solve the equilibrium and the linear vibration problems. The unknowns of the problem (solution branch, frequency and eigenmode) are determined by a perturbation technique whose terms are calculated by the finite element method.

For a given applied load, the structure is assumed to oscillate around a static state $U$. These oscillations are described by the time-dependent mixed vector $V(\theta_1, \theta_2, \theta_3, t)$, so the global response is written $U(\theta_1, \theta_2, \theta_3) + V(\theta_1, \theta_2, \theta_3, t)$. The oscillations of the structure around a static equilibrium state are described by the following equation:

$$M(\ddot{V}) + L(U + V) + Q(U + V, U + V) = \lambda F \quad (37)$$

with

$$<M\ddot{V}, \delta U> = \int_{\Omega} \rho(\ddot{v}_1\delta u_1 + \ddot{v}_2\delta u_2 + \ddot{v}_3\delta u_3)d\Omega$$

where $\ddot{V}$ is the second derivative with respect to time of $V$. $v_i$ are the components of the displacement $V$, $M$ is the mass matrix and $\rho$ is the density.

Neglecting the quadratic term in $V(Q(V, V))$, Eq.(37) becomes:

$$L_t(V) + M(\ddot{V}) = 0 \quad (38)$$

where $L_t = L(.) + 2Q(U,.)$ is the tangent operator at a static deformed state $U$. The general solution of (38) is expressed as:

$$V = Xe^{i\omega t} \quad (39)$$

By inserting Eq.(39) into Eq.(38), we obtain the linear eigenvalue problem:

$$L_t(X) - w^2 M(X) = 0 \quad (40)$$

To study the eigenfrequencies $w$ and their associated modes $X$ in pre-buckling and post-buckling ranges, we need to solve the following coupled problem:
$$
\begin{equation}
\begin{cases}
L(U) + Q(U, U) = \lambda F \\
L_t(X) = w^2 M(X)
\end{cases}
\tag{41a}
\end{equation}
$$

Eq.(41a) corresponds to the static displacement load solution (Eq.(15)).
The obtained displacement **U** is used to determine the tangent operator **L**<sub>t</sub>.
Note that, if the static equilibrium **U**(**λ**) is stable, all the eigenvalues of
(41b) are positive and **w** represents the natural frequencies. On the other
hand, if **U**(**λ**) is unstable, (41b) has negative eigenvalues. The critical points
correspond to the load value for which **w** = 0. This criterion will be used to
detect the bifurcation points.

The static problem (41a) has been solved using the ANM in section 2.4.
The same technique will be used to solve the dynamic problem (41b). The
idea is to determine analytically in power series form the unknowns (X, w)

$$
\begin{equation}
\left\{
\begin{aligned}
\mathbf{X} &= \mathbf{X}_0 + a\mathbf{X}_1 + a^2\mathbf{X}_2 + \dots + a^n\mathbf{X}_n \\
w^2 &= w_0^2 + aw_1 + a^2w_2 + \dots + a^n w_n
\end{aligned}
\right.
\tag{42}
\end{equation}
$$

where (X₀, w₀²) is the solution of the eigenvalue problem at the starting
equilibrium point (U₀, λ₀), Xᵣ = (uₓr, γₓr, Sₓr) is a mixed vector and wᵣ
is a scalar parameter. The parameter a is the same as in problem (41a).
Substituting Eq.(42) into Eq.(41b) and identifying the terms according to
the power of a, we obtain a sequence of linear problems:

Order 1:

$$
\mathbf{L}_t^0(\mathbf{X}_1) - w_0^2 M \mathbf{X}_1 = w_1 M \mathbf{X}_0 - \mathbf{Q}(U_1, \mathbf{X}_0) - \mathbf{Q}(\mathbf{X}_0, U_1) \quad (43)
$$

Order r:

$$
\mathbf{L}_t^0(\mathbf{X}_r) - w_0^2 M \mathbf{X}_r = w_r M \mathbf{X}_0 + \sum_{j=1}^{r-1} w_j M \mathbf{X}_{r-j} - \sum_{j=1}^{r} (Q(U_j, \mathbf{X}_{r-j}) + Q(\mathbf{X}_{r-j}, U_j)) \quad (44)
$$

Again, we see that all these problems have the same linear operator ($L_t^0$ - $w_0^2 M$) which must be decomposed only once for all $\mathbf{X}_r$ ($r = 0, 1, ..., n$). The projection of (43) and (44) into the vibration mode $\mathbf{X}_0$ gives:

Order 1:

$$
w_1 = \frac{\langle \mathbf{X}_0, Q(U_1, \mathbf{X}_0) \rangle + \langle \mathbf{X}_0, Q(\mathbf{X}_0, U_1) \rangle}{\langle M\mathbf{X}_0, \mathbf{X}_0 \rangle} \quad (45)
$$
Order $r$:

$$w_r = \frac{\sum_{j=1}^{r} (\langle \mathbf{X}_0, \mathbf{Q}(U_j, \mathbf{X}_{r-j}) \rangle + \langle \mathbf{X}_0, \mathbf{Q}(\mathbf{X}_{r-j}, U_j) \rangle) - \sum_{j=1}^{r-1} w_j \langle \mathbf{X}_0, M\mathbf{X}_{r-j} \rangle}{\langle M\mathbf{X}_0, \mathbf{X}_0 \rangle} \quad (46)$$

The finite element discretization of the vibration problem gives:

Order 1:

$$\left\{ \begin{array}{l} w_1 = \frac{t\{x_0\}\{F_1\}}{t\{x_0\}[M]\{x_0\}} \\ ([K_t^0] - w_0^2[M])\{x_1\} = w_1[M]\{X_0\} + \{F_1\} \end{array} \right. \quad (47)$$

Order $r$:

$$\left\{ \begin{array}{l} w_r = \frac{t\{x_0\}\{F_r\} - \sum_{j=1}^{r-1} w_j^t\{x_0\}[M]\{x_{r-j}\}}{t\{x_0\}[M]\{x_0\}} \\ ([K_t^0] - w_0^2[M])\{x_r\} = w_r[M]\{X_0\} + \sum_{j=1}^{r-1} w_j[M]\{x_{r-j}\} + \{F_r\} \end{array} \right. \quad (48)$$

where $\{x_r\}$ is the discretized form of $u_{xr}$ and the vectors $\{F_r\}$ depend only on the solutions up to order $(r-1)$.

## 4. Numerical applications

In this section, a bifurcation analysis is conducted for a multilayer structure with nonlinear pre-buckling to evaluate the efficiency and robustness of bifurcation detection algorithms. We consider a cylindrical laminated roof composed of three orthotropic layers (Figure 2). The thickness of each layer is the same. The top and bottom layers are oriented at 45° and the middle layer is at 0° (45°/0°/45°). The material characteristics are given in Table 1. The roof is pinned on the two straight edges and free on the two others edges and subjected to point loading at its center. The structure is discretized by 36 shell elements uniformly distributed.

The response curve of the structure using the rational representation is given in Figure 3. Figure 4 shows the evolution of the deformation of the structure along the equilibrium path of Figure 3. To reach a vertical displacement of the load application point of 30 mm, the ANM requires 35
<figure><img src="image_3.png" /><figcaption>Figure 2: Geometrical parameters of the cylindrical roof.</figcaption></figure>

<table><caption>Table 1: Mechanical properties of each roof layer.</caption><thead><tr><th></th><th>E(MPa)</th><th>μ</th><th>G(MPa)</th></tr></thead><tbody><tr><td>E<sub>1</sub></td><td>3300</td><td>μ<sub>12</sub></td><td>0.25</td><td>G<sub>12</sub></td><td>660</td></tr><tr><td>E<sub>2</sub></td><td>1100</td><td>μ<sub>23</sub></td><td>0.25</td><td>G<sub>23</sub></td><td>440</td></tr><tr><td>E<sub>3</sub></td><td>1100</td><td>μ<sub>13</sub></td><td>0.25</td><td>G<sub>13</sub></td><td>660</td></tr></tbody></table>

steps (35 matrix decompositions) for truncation order $n = 15$ and accuracy parameter $\delta = 10^{-6}$ (the maximum residual of the solution path is equal to $log_{10}(Res) = -2.32$). A summary of the number of steps to have ($u \approx 30$ mm) according to $n$ and $\delta$ is given in Table 2. The same problem is simulated using Abaqus with Newton Raphson algorithm associated to Riks technique. The mesh is the same as the one used for the ANM with an eight node element S8R available in Abaqus implicit. A displacement of 30 mm from the load application point requires 124 steps with 374 matrix decompositions by using a fixed increment equal to 0.05 (or 62 steps with 203 matrix decompositions by using 0.1). This comparison shows the efficiency of the ANM.

Depending on the parameters $n$ and $\delta$ and on the technique used to represent the solution (polynomial or rational representation), the asymptotic numerical method allows to have other equilibrium paths as shown in Figure 5. For $n = 15$ and $\delta = 10^{-6}$ and with a polynomial representation, we notice a deviation of the curve at point C and then another deviation at point B, both of which are bifurcation points. The curve continues to follow
<figure><img src="image_4.png" /><figcaption>Figure 3: Load - vertical displacement of the load application point ($n = 15$ and $\delta = 10^{-6}$). Black line represents the response given by ANM with Padé approximation, circles are relative to the end of each ANM-Padé step and triangles represent 10 points of Abaqus-Riks response (the total number of steps is 124).</figcaption></figure>

the same path between C and B. Keeping the same representation, another equilibrium path is observed with $n = 20$ and $\delta = 10^{-8}$. Points A and F are also bifurcation points of the structure. Two other bifurcation points (points D and E) appear with a rational representation and parameters $n = 15$ and $\delta = 10^{-10}$. Thus, we have six bifurcation points when we try to determine them using the three techniques described in Section 3. It should be noted that when the structure reaches these bifurcation points, it becomes unstable and an important change in its shape may occur.

The bifurcation indicator detects the six bifurcation points (Figure 6). It becomes zero exactly at the critical points. Table 3 shows the critical loading, detection step and residual at the end of the step for each bifurcation point, respectively. The residual remains quite small along the equilibrium branch.
<figure><img src="image_5.png" /><figcaption>Figure 4: The deformations of the structure at steps (a) 20, (b) 50, (c) 90 and (d) 120 in Figure 3.</figcaption></figure>

The second technique for the detection of bifurcation points is the pole analysis of the Padé approximants. Table 4 presents the real roots of the denominator of the Padé fractions in the vicinity of the bifurcation points for different truncation order $n$ and accuracy parameter $\delta$. The poles shown are those belonging to the loading interval and the three configurations analyzed give the same response of Figure 3. We notice that the bifurcation points detected by the bifurcation indicator are among the poles. Taking
<table><caption>Table 2: Summary of the number of ANM-Padé steps to obtain a vertical displacement of the load application point of 30 mm as a function of the truncation order *n* and the accuracy *δ*. The symbol (-) indicates that the corresponding parameters *n* and *δ* do not give the same equilibrium path described in Figure 3.</caption><thead><tr><th>δ<br>n</th><th>10<sup>-6</sup></th><th>10<sup>-8</sup></th><th>10<sup>-10</sup></th><th>10<sup>-12</sup></th></tr></thead><tbody><tr><th>15</th><td>35</td><td>47</td><td>-</td><td>-</td></tr><tr><th>20</th><td>-</td><td>36</td><td>52</td><td>-</td></tr><tr><th>25</th><td>-</td><td>-</td><td>44</td><td>-</td></tr><tr><th>30</th><td>-</td><td>-</td><td>32</td><td>48</td></tr></tbody></table>

<table><caption>Table 3: Critical loading, detection step and relative residual *Res* = $\frac{||F_{int} - F_{ext}||}{||F_{ext}||}$ at each bifurcation point.</caption><thead><tr><th>Bifurcation point</th><th>A</th><th>B</th><th>C</th><th>D</th><th>E</th><th>F</th></tr></thead><tbody><tr><th>critical loading λ<sub>c</sub></th><td>0.252</td><td>0.273</td><td>-0.054</td><td>-0.184</td><td>-0.199</td><td>-0.111</td></tr><tr><th>Step</th><td>2</td><td>4</td><td>9</td><td>23</td><td>26</td><td>34</td></tr><tr><th>Log<sub>10</sub>(Res)</th><td>-4.683</td><td>-4.448</td><td>-3.923</td><td>-4.383</td><td>-4.421</td><td>-2.319</td></tr></tbody></table>

the example of *n* = 15 and *δ* = 10⁻⁶, the second step gives us a single pole (*p*₂ = 0.252) that belongs to the loading interval and coincides exactly with the critical loading relative to the bifurcation point A. The same analysis for the fourth step leads to point B. Nevertheless, in the ninth step we find two poles (*p*₉₁ = -0.032) and (*p*₉₂ = -0.054 = λₑ(*C*)). This can be attributed to defects in the rational representation because the pole *p*₉₁ no longer appears by changing *n* and *δ*. The same observation appears in the other situations where we find more than one pole. Thus, Padé approximants offer a useful tool for identifying bifurcation points without incurring any additional computational costs. To obtain reliable results, it is recommended to explore multiple configurations.

The third technique for bifurcation analysis is to track the zero values of the natural frequencies along the equilibrium path. Figure 7 shows the evolution of the relative frequency *β* for the first vibration mode as a function of the loading parameter *λ*. As with the bifurcation indicator, *β* is zero at the critical points and accurately detects the six bifurcation points.

In multilayer structures, the bifurcation points depend on the number of layers, the orientation, thickness and mechanical properties of each layer. Thus, we tested two other roof configurations by changing the thickness (*h* =
<figure><img src="image_6.png" /><figcaption>Figure 5: Other equilibrium paths of the structure as a function of representation type, truncation order *n* and accuracy parameter *δ*. (a) *n* = 15, *δ* = 10<sup>-6</sup>, series representation (b) *n* = 20, *δ* = 10<sup>-8</sup>, series representation (c) *n* = 15, *δ* = 10<sup>-10</sup>, Padé representation (d) a zoom of Figure 4.c (*u* ∈ [10; 20] and *λ* ∈ [-0.3; 0]).</figcaption></figure>

9 mm) and then the orientation of the layers while keeping constant the other parameters. Figure 8 shows the response curves for the two configurations as well as the bifurcation points, obtained using the techniques discussed previously, with the corresponding critical loads. These tests confirm the effectiveness of our algorithms for following equilibrium paths and detecting bifurcations.

## 5. Conclusion

In this work, we have presented a bifurcation analysis for multilayer structures using the asymptotic numerical method. Three bifurcation detection
<figure><img src="image_7.png" /><figcaption>Figure 6: Bifurcation indicator, along the equilibrium branch, as a function of the loading parameter for $n = 15$ and $\delta = 10^{-6}$ (Padé). Points A, B, C, D, E and F are the six bifurcation points of the structure.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Evolution of the relative natural frequency $\beta = \frac{w^2}{w_2^2}$, along the equilibrium branch of Figure 3, as a function of the loading $\lambda$ for the first vibration mode ($n = 15$ and $\delta = 10^{-6}$ with $\lambda = 0$). Points A, B, C, D, E and F are the six bifurcation points of the structure.</figcaption></figure>
<table><caption>Table 4: Poles of the Padé approximants for three couples of *n* and *δ* in the vicinity of the bifurcation points.</caption><thead><tr><th></th><th colspan="6">n=15</th></tr><tr><th></th><th colspan="6">δ = 10<sup>-6</sup></th></tr><tr><th>Step</th><th>2</th><th>4</th><th>9</th><th>23</th><th>26</th><th>34</th></tr></thead><tbody><tr><th>Loading interval</th><td>[0.221;0.282]</td><td>[0.212;0.306]</td><td>[-0.061;-0.026]</td><td>[-0.191;-0.182]</td><td>[-0.199;-0.197]</td><td>[-0.165;-0.001]</td></tr><tr><th>λ<sub>c</sub></th><td>0.252</td><td>0.273</td><td>-0.032<br/>-0.054</td><td>-0.184</td><td>-0.199</td><td>-0.111</td></tr></tbody><thead><tr><th></th><th colspan="6">n=20</th></tr><tr><th></th><th colspan="6">δ = 10<sup>-10</sup></th></tr><tr><th>Step</th><th>2</th><th>5</th><th>13</th><th>30</th><th>39</th><th>49</th></tr></thead><tbody><tr><th>Loading interval</th><td>[0.185;0.274]</td><td>[0.263;0.299]</td><td>[-0.055;-0.049]</td><td>[-0.186;-0.184]</td><td>[-0.199;-0.182]</td><td>[-0.197;-0.077]</td></tr><tr><th>λ<sub>c</sub></th><td>0.252</td><td>0.273</td><td>-0.054</td><td>-0.184</td><td>-0.184<br/>-0.199</td><td>-0.111<br/>-0.130</td></tr></tbody><thead><tr><th></th><th colspan="6">n=30</th></tr><tr><th></th><th colspan="6">δ = 10<sup>-12</sup></th></tr><tr><th>Step</th><th>2</th><th>6</th><th>18</th><th>33</th><th>37</th><th>46</th></tr></thead><tbody><tr><th>Loading interval</th><td>[0.206;0.265]</td><td>[0.273;0.278]</td><td>[-0.056;-0.047]</td><td>[-0.189;-0.183]</td><td>[-0.199;-0.193]</td><td>[-0.203;-0.049]</td></tr><tr><th>λ<sub>c</sub></th><td>0.252</td><td>0.273<br/>0.274</td><td>-0.054</td><td>-0.184</td><td>-0.199<br/>-0.198</td><td>-0.111</td></tr></tbody></table>

techniques have been recalled and they allowed us to conduct a deep study on the instability of thin shell structures. The first technique is based on the bifurcation indicator which is computed along the equilibrium branch and which vanishes exactly at the critical points. The use of a rational representation (Padé approximants) instead of a polynomial representation (Taylor series) allows both to improve the validity range of the solution and to detect the bifurcation by analyzing the roots of the denominator of the rational fractions. A combination of buckling and vibration describes the third technique and the natural frequencies of vibration allow to detect bifurcations. According to the numerical tests realized, the three techniques are able to detect the critical points with a high accuracy whatever the truncation order $15 \le n \le 30$ and the parameter $10^{-12} \le \delta \le 10^{-6}$. The ANM is therefore a very efficient algorithm to study instabilities and detect bifurcation points in nonlinear pre-buckling branches, which is not the case for classical iterative algorithms.
<figure><img src="image_9.png" /></figure>

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

<table><thead><tr><th>Bifurcation point</th><th>G</th><th>H</th><th>I</th><th>Bifurcation point</th><th>J</th><th>K</th><th>L</th><th>M</th></tr></thead><tbody><tr><td>λ<sub>c</sub></td><td>0.520</td><td>-0.012</td><td>-0.260</td><td>λ<sub>c</sub></td><td>0.354</td><td>0.366</td><td>-0.102</td><td>-0.151</td></tr></tbody></table>

Figure 8: Response curve for two different roof configurations by changing (a) the thickness (h=9mm) (b) the orientation of the layers (90°/0°/90°). The green points represent the bifurcation points and the corresponding table shows the critical loads.

## Appendix: Finite element discretization

The fields $v$ and $w$, which represent respectively the displacement and director difference, are interpolated via the shape functions of the classical eight node serendipity quadrilateral. By collecting the nodal values of $v$ and $w$ in the vector $\{q\}$, the displacement $\{u\}$, its virtual part $\{\delta u\}$ and its gradient $\{\theta(u)\}$ are related to nodal displacements $q$ and $\delta q$ as follows:

$$ \{u\} = [N]\{q\} \qquad \{\delta u\} = [N]\{\delta q\} \qquad \{\theta(u)\} = [G]\{q\} \qquad (49) $$

where $[N]$ is the matrix of shape functions and $[G]$ is their gradient matrix.

The Green-Lagrange strain (compatible part) is expressed in the covariant basis as:

$$ \gamma_{ij}(u) = \frac{1}{2} \left( \frac{\partial u}{\partial \theta^i} g_j + \frac{\partial u}{\partial \theta^j} g_i + \frac{\partial u}{\partial \theta^i} \frac{\partial u}{\partial \theta^j} \right) \qquad (50) $$

where $g_i$ represents the covariant base vectors. $\gamma(u)$ can be decomposed into
a linear part $\gamma_l(u)$ and a non-linear part $\gamma_{nl}(u, u)$:

$$ \gamma(u) = \gamma_l(u) + \gamma_{nl}(u, u) \quad (51) $$

with

$$ \left\{ \begin{array}{l} \gamma_l = [\mathbf{R}]\{\theta(q)\} = [\mathbf{R}][\mathbf{G}]\{q\} \\ \gamma_{nl} = \frac{1}{2}[\mathbf{A}(q)]\{\theta(q)\} = \frac{1}{2}[\mathbf{A}(q)][\mathbf{G}]\{q\} \end{array} \right. \quad (52) $$

The virtual strain $\delta\gamma$ is written as follows:

$$ \delta\gamma = ([\mathbf{R}] + [\mathbf{A}(q)])[\mathbf{G}]\{\delta q\} = [\overline{\mathbf{B}}]\{\delta q\} \quad (53) $$

with

$$ [\overline{\mathbf{B}}] = ([\mathbf{R}] + [\mathbf{A}(q)])[\mathbf{G}] $$

$[\mathbf{R}]$ and $[\mathbf{A}(q)]$ represent the component matrix of the covariant basis and the displacement gradient matrix, respectively.

For the enhanced assumed strain $\tilde{\gamma} = \theta_3 \tilde{\beta}_{33} g^3 \otimes g^3$, the component $\tilde{\beta}_{33}$ does not require any inter-element continuity. Thus, it is discretized as a bilinear polynomial:

$$ \tilde{\beta}_{33} = \alpha_1 + \alpha_2\xi + \alpha_3\eta + \alpha_4\xi\eta \quad (54) $$

where $\xi$ and $\eta$ represent the isoparametric coordinates. The unknown parameters $\alpha_1, \alpha_2, \alpha_3$ and $\alpha_4$ can be eliminated at the elementary level [28]. Thus, the additional deformation is written in the following matrix form:

$$ \{\tilde{\gamma}\} = [\mathbf{B}_{\alpha}]\{\alpha\} \quad (55) $$

345

By substituting Eqs (52), (53) and (55) into Eq. (23), we obtain the discretization form which represents the equilibrium of global structure

$$ \left\{ \begin{array}{l} [\mathbf{K}_t^0]\{\mathbf{q}_p\} = \lambda_p \mathbf{F} + \mathbf{F}_p^{nl} - [\mathbf{K}_{\alpha u}][\mathbf{K}_{\alpha \alpha}]^{-1}\{\mathbf{R}_p^{\alpha}\} \\ [\mathbf{K}_{\alpha \alpha}]\{\alpha_p\} + [\mathbf{K}_{\alpha u}]\{\mathbf{q}_p\} = \{\mathbf{R}_p^{\alpha}\} \\ {}^t\{\mathbf{q}_p\} \cdot \{\mathbf{q}_1\} + \lambda_p \lambda_1 = 0 \end{array} \right. \quad (56) $$

with:

$$ [\mathbf{K}_{\alpha\alpha}] = \sum_{m=1}^{N} \int_{\Omega} ({}^t[\mathbf{B}_{\alpha}][\mathbf{D}]^m[\mathbf{B}_{\alpha}]) dv $$
$$[K_{\alpha u}] = \sum_{m=1}^{N} \int_{\Omega} ({}^{t}[B_{\alpha}][D]^{m}[\overline{B}]) dv$$

$$[K_{uu}] = \sum_{m=1}^{N} \int_{\Omega} ({}^{t}[\overline{B}][D]^{m}[\overline{B}] + {}^{t}[G][M]^{m}[G]) dv$$

$$[K_t^0] = [K_{uu}] - [K_{\alpha u}][K_{\alpha \alpha}]^{-1}[K_{\alpha u}]$$

where $[K_t^0]$ is the tangent stiffness matrix, $\{F\}$ is the external load vector and $[M]$ the initial stress matrix. Compared to Riks method, the new terms are the right-hand nonlinear terms $\{F_p^{nl}\}$ and $\{R_p^{\alpha}\}$ related to solutions at previous orders $p-1$. These two terms are expressed as follows

$$\{F_p^{nl}\} = - \sum_{m=1}^{N} \left\{ \int_{\Omega} ({}^{t}[G] \sum_{r=1}^{p-1} {}^{t}[A(q_{p-r})] \{S_r\})^{m} + {}^{t}[\overline{B}(q_0)][D]^{m} \sum_{r=1}^{p-1} \frac{1}{2} [A(q_{p-r})] \{\theta(q_r)\} \right\} dv$$

$$\{R_p^{\alpha}\} = - \sum_{m=1}^{N} \left\{ \int_{\Omega} ({}^{t}[B_{\alpha}][D]^{m} \sum_{r=1}^{p-1} \frac{1}{2} [A(q_{p-r})] \{\theta(q_r)\}) dv \right\}$$
References

[1] Z. S. Toor, Space applications of composite materials, Journal of Space Technology 8 (1) (2018) 65–70.

350

[2] Y. Wang, K. H. Low, H. L. J. Pang, K. H. Hoon, F. X. Che, Y. S. Yong, Modeling and simulation for a drop-impact analysis of multi-layered printed circuit boards, Microelectronics Reliability 46 (2-4) (2006) 558–573.

[3] D. Kim, S. Joo, D. Kwak, H. Kim, Anisotropic viscoelastic shell modeling technique of copper patterns/photoimageable solder resist composite for warpage simulation of multi-layer printed circuit boards, Journal of Micromechanics and Microengineering 25 (10) (2015) 105016.

355

[4] C. G. Diaconu, P. M. Weaver, F. Mattioni, Concepts for morphing airfoil sections using bi-stable laminated composite structures, Thin-Walled Structures 46 (6) (2008) 689–701.

360

[5] E. Ramm, Buckling of Shells: Proceedings of a State-of-the-Art Colloquium, Universität Stuttgart, Germany, May 6–7, 1982, Springer Science & Business Media, 2012.

[6] K. K. Yadav, S. Gerasimidis, Instability of thin steel cylindrical shells under bending, Thin-Walled Structures 137 (2019) 151–166.

365

[7] M. Brunetti, S. Vidoli, A. Vincenti, Bistability of orthotropic shells with clamped boundary conditions: an analysis by the polar method, Composite Structures 194 (2018) 388–397.

[8] M. Finot, S. Suresh, Small and large deformation of thick and thin-film multi-layers: effects of layer geometry, plasticity and compositional gradients, Journal of the Mechanics and Physics of Solids 44 (5) (1996) 683–721.

370

[9] A. El Chebair, M. P. Païdoussis, A. K. Misra, Experimental study of annular-flow-induced instabilities of cylindrical shells, Journal of Fluids and Structures 3 (4) (1989) 349–364.

375
[10] E. Frulloni, J. M. Kenny, P. Conti, L. Torre, Experimental study and finite element analysis of the elastic instability of composite lattice structures for aeronautic applications, *Composite Structures* 78 (4) (2007) 519–528.

380

[11] T. Fu, X. Wu, Z. Xiao, Z. Chen, Dynamic instability analysis of FG-CNTRC laminated conical shells surrounded by elastic foundations within FSDT, *European Journal of Mechanics/A Solids* 85 (2021) 104139.

[12] Z. Kuang, Q. Huang, W. Huang, J. Yang, H. Zahrouni, M. Potier-Ferry, H. Hu, A computational framework for multi-stability analysis of laminated shells, *Journal of the Mechanics and Physics of Solids* 149 (2021) 104317.

385

[13] N. Damil, M. Potier-Ferry, A new method to compute perturbed bifurcations: application to the buckling of imperfect elastic structures, *International Journal of Engineering Science* 28 (9) (1990) 943–957.

390

[14] B. Cochelin, A path-following technique via an asymptotic-numerical method, *Computers & Structures* 53 (5) (1994) 1181–1192.

[15] H. Zahrouni, B. Cochelin, M. Potier-Ferry, Computing finite rotations of shells by an asymptotic-numerical method, *Computer methods in applied mechanics and engineering* 175 (1-2) (1999) 71–85.

395

[16] H. Abichou, H. Zahrouni, M. Potier-Ferry, Asymptotic numerical method for problems coupling several nonlinearities, *Computer methods in applied mechanics and engineering* 191 (51-52) (2002) 5795–5810.

[17] M. Assidi, H. Zahrouni, N. Damil, M. Potier-Ferry, Regularization and perturbation technique to solve plasticity problems, *International Journal of Material Forming* 2 (1) (2009) 1–14.

400

[18] J. M. Cadou, M. Potier-Ferry, B. Cochelin, N. Damil, ANM for stationary Navier-Stokes equations and with Petrov-Galerkin formulation, *International Journal for Numerical Methods in Engineering* 50 (4) (2001) 825–845.

405
[19] A. Elhage-Hussein, M. Potier-Ferry, N. Damil, A numerical continuation method based on Padé approximants, International Journal of Solids and Structures 37 (46-47) (2000) 6981-7001.

410

[20] R. Xu, Y. Hui, H. Hu, Q. Huang, H. Zahrouni, T. Ben Zineb, M. Potier-Ferry, A Fourier-related FE2 multiscale model for instability phenomena of long fiber reinforced materials, Composite Structures 211 (2019) 530-539.

[21] G. A. Baker, P. Graves-Morris, Padé approximants. part 1: Basic theory, Encyclopedia of Mathematics and Its Applications (1981).

415

[22] A. Najah, B. Cochelin, N. Damil, M. Potier-Ferry, A critical review of asymptotic numerical methods, Archives of Computational Methods in Engineering 5 (1) (1998) 31-50.

[23] A. K. Noor, J. M. Peters, Reduced basis technique for nonlinear analysis of structures, AIAA Journal 18 (4) (1980) 455-462.

420

[24] E. H. Boutyour, H. Zahrouni, M. Potier-Ferry, M. Boudi, Bifurcation points and bifurcated branches by an asymptotic numerical method and Padé approximants, International Journal for Numerical Methods in Engineering 60 (12) (2004) 1987-2012.

[25] A. Tri, H. Zahrouni, M. Potier-Ferry, Bifurcation indicator based on meshless and asymptotic numerical methods for nonlinear poisson problems, Numerical Methods for Partial Differential Equations 30 (3) (2014) 978-993.

425

[26] M. Jamal, H. Elasmar, B. Braikat, E. H. Boutyour, B. Cochelin, N. Damil, M. Potier-Ferry, Bifurcation indicators, Acta Mechanica 139 (1) (2000) 129-142.

430

[27] E. H. Boutyour, L. Azrar, M. Potier-Ferry, Vibration of buckled elastic structures with large rotations by an asymptotic numerical method, Computers & Structures 84 (3-4) (2006) 93-101.

[28] N. Büchter, E. Ramm, D. Roehl, Three-dimensional extension of nonlinear shell formulation based on the enhanced assumed strain concept, International Journal for Numerical Methods in Engineering 37 (15) (1994) 2551-2568.

435
[29] J. C. Simo, M. S. Rifai, A class of mixed assumed strain methods and the method of incompatible modes, International Journal for Numerical Methods in Engineering 29 (8) (1990) 1595–1638.

[30] B. Cochelin, N. Damil, M. Potier-Ferry, The asymptotic-numerical method: an efficient perturbation technique for nonlinear structural mechanics, Revue Européenne des Éléments Finis 3 (2) (1994) 281–297.