Hamza Azzayani
email: hamza.azzayani@univ-lorraine.fr

Hamid Zahrouni
email: hamid.zahrouni@univ-lorraine.fr

Norman Mathieu
email: norman.mathieu@univ-lorraine.fr

Pascal Ventura
email: pascal.ventura@univ-lorraine.fr

Michael Brun
email: michael.brun@univ-lorraine.fr

Michel Potier-Ferry
email: michel.potier-ferry@univ-lorraine.fr

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

Keywords: Asymptotic numerical method, instability, bifurcation, Padé approximants, thin multilayer structures

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.

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 [START_REF] Toor | Space applications of composite materials[END_REF], printed circuit boards [START_REF] Wang | Modeling and simulation for a drop-impact analysis of multi-layered printed circuit boards[END_REF][START_REF] Kim | Anisotropic viscoelastic shell modeling technique of copper patterns/photoimageable solder resist composite for warpage simulation of multi-layer printed circuit boards[END_REF], morphing airfoils [START_REF] Diaconu | Concepts for morphing airfoil sections using bi-stable laminated composite structures[END_REF],...). However, these structures are very sensitive to instabilities (e.g. buckling) [START_REF] Ramm | Buckling of Shells: Proceedings of a State-of-the-Art Colloquium[END_REF][START_REF] Yadav | Instability of thin steel cylindrical shells under bending[END_REF],

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 [START_REF] Brunetti | Bistability of orthotropic shells with clamped boundary conditions: an analysis by the polar method[END_REF][START_REF] Finot | Small and large deformation of thick and thinfilm multi-layers: effects of layer geometry, plasticity and compositional gradients[END_REF]. 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 [START_REF] Chebair | Experimental study of annular-flow-induced instabilities of cylindrical shells[END_REF][START_REF] Frulloni | Experimental study and finite element analysis of the elastic instability of composite lattice structures for aeronautic applications[END_REF][START_REF] Fu | Dynamic instability analysis of FG-CNTRC laminated conical shells surrounded by elastic foundations within FSDT[END_REF]. 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 [START_REF] Kuang | A computational framework for multi-stability analysis of laminated shells[END_REF]. 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 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.

The asymptotic numerical method (ANM) based on the perturbation technique presents an alternative [START_REF] Damil | A new method to compute perturbed bifurcations: application to the buckling of imperfect elastic structures[END_REF]. 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 [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF][START_REF] Zahrouni | Computing finite rotations of shells by an asymptotic-numerical method[END_REF]. 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 [START_REF] Zahrouni | Computing finite rotations of shells by an asymptotic-numerical method[END_REF]. 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 [START_REF] Abichou | Asymptotic numerical method for problems coupling several nonlinearities[END_REF][START_REF] Assidi | Regularization and perturbation technique to solve plasticity problems[END_REF], Navier-Stockes equation [START_REF] Cadou | ANM for stationary Navier-Stokes equations and with Petrov-Galerkin formulation[END_REF], contact mechanics [START_REF] Elhage-Hussein | A numerical continuation method based on Padé approximants[END_REF], instability of fiber reinforced composites [START_REF] Xu | A Fourier-related FE2 multiscale model for instability phenomena of long fiber reinforced materials[END_REF],... The use of a rational representation (Padé approximants [START_REF] Baker | Padé approximants. part 1: Basic theory[END_REF]) instead of a polynomial representation improves the validity range of the solution [START_REF] Elhage-Hussein | A numerical continuation method based on Padé approximants[END_REF]. Najah et al. [START_REF] Najah | A critical review of asymptotic numerical methods[END_REF] have shown the efficiency of the Padé approximants over the power series representation and the Rayleigh-Ritz technique [START_REF] Najah | A critical review of asymptotic numerical methods[END_REF][START_REF] Noor | Reduced basis technique for nonlinear analysis of structures[END_REF]. 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.

Within the framework of ANM, bifurcation indicators have been proposed to identify critical points. A first indicator is a scalar function computed by introducing a fictitious perturbation force in the equilibrium problem [START_REF] Boutyour | Bifurcation points and bifurcated branches by an asymptotic numerical method and Padé approximants[END_REF]. The indicator can be determined explicitly along the equilibrium branch by the perturbation technique. The roots of this function characterize the singular points [START_REF] Tri | Bifurcation indicator based on meshless and asymptotic numerical methods for nonlinear poisson problems[END_REF]. Several applications of this technique have been implemented in solid and fluid mechanics [START_REF] Cadou | ANM for stationary Navier-Stokes equations and with Petrov-Galerkin formulation[END_REF][START_REF] Jamal | Bifurcation indicators[END_REF]. 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 [START_REF] Boutyour | Bifurcation points and bifurcated branches by an asymptotic numerical method and Padé approximants[END_REF]. 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) [START_REF] Boutyour | Vibration of buckled elastic structures with large rotations by an asymptotic numerical method[END_REF].

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.

Static shell formulation and resolution algorithm

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 [START_REF] Büchter | Three-dimensional extension of nonlinear shell formulation based on the enhanced assumed strain concept[END_REF]. 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(θ 1 , θ 2 , θ 3 ) = r(θ 1 , θ 2 ) + θ 3 a 3 (θ 1 , θ 2 ) ( 1
)
where r is the projection of this point on the mid-surface, a 3 the shell director vector and (θ 1 , θ 2 , θ 3 ) the curvilinear coordinates. Then, the covariant

g α = ∂x ∂θα = r, α +θ 3 a 3 , α α = 1, 2 g 3 = ∂x ∂θ 3 = a 3 (2)
where ,α = ∂ ∂θα Assuming a linear variation of the displacement in the thickness direction, the displacement field corresponding to the deformed configuration can be written as:

u(θ 1 , θ 2 , θ 3 ) = v(θ 1 , θ 2 ) + θ 3 w(θ 1 , θ 2 ) (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:

x = (r + v) + θ 3 (a 3 + w) = r + θ 3 a 3 (4)
In the same way as Eq.( 2), the covariant base vectors in the deformed configuration are expressed as:

g α = ∂x ∂θα = r, α +θ 3 a 3 , α α = 1, 2 g 3 = ∂x ∂θ 3 = a 3 (5)
The reader wishing more details on this formulation can refer to the work 100 of Büchter et al. [START_REF] Büchter | Three-dimensional extension of nonlinear shell formulation based on the enhanced assumed strain concept[END_REF].

Strain and stress equations

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

γ = 1 2 (g i j -g ij )g i ⊗ g j g ij = g i .g j (6)
where g ij and g i j 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 .g j = δ i j . From Eqs (2) and ( 5), the components of the Green-Lagrange strain can be written as:

γ ij = α ij + β ij θ 3 + Q ij (θ 3 ) 2 (7)
with

                         α ij = 1 2 (a i a j -a i a j ) β αβ = 1 2 (a 3,α a β + a 3,β a α -a 3,α a β -a 3,β a α ) β α3 = 1 2 (a 3,α a 3 -a 3,α a 3 ) β 33 = 0 Q αβ = 1 2 (a 3,α a 3,β -a 3,α a 3,β ) Q α3 = 0 Q 33 = 0 (8)
To avoid numerical locking drawback in shell modeling, the EAS (Enhanced Assumed Strain) concept proposed by Simo and Rifai [START_REF] Simo | A class of mixed assumed strain methods and the method of incompatible modes[END_REF] will be used. This method consists in introducing an additional strain field γ 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 γ is written:

γ = θ 3 β 33 g 3 ⊗ g 3 Ω ( t S m : γ)dΩ = 0 (9)
The term β 33 provides a linear variation of the strain field through the thickness. Thus, the total strain field γ and the stress field (second Piola-Kirchhoff stress tensor) in each layer S m are written:

γ = γ(u) + γ = (γ l (u) + γ nl (u, u)) + γ (10)
S m = D m : γ (11)
in which γ l (u) and γ nl (u, u) denote the linear and nonlinear parts of the compatible strain field respectively and D m = D ij kl(m) g i ⊗ g j ⊗ g k ⊗ g l refers to the fourth-order elasticity tensor of the mth layer. Since the elasticity tensor is often obtained in cartesian coordinates, a transformation to curvilinear coordinates is used in the numerical implementation

D ij kl(m) = C pqrs(m) (e p .g i )(e q .g j )(e r .g k )(e s .g l ) (12)
where C pq rs(m) presents the components of elasticity tensor in cartesian coordinates.

A three-field variational formulation

Based on the Hu-Washizu functional, we consider a shell formulation with three independent fields (u, S m , γ) [START_REF] Zahrouni | Computing finite rotations of shells by an asymptotic-numerical method[END_REF]. It is written in the following form:

π E AS (u, γ, S m ) = N m=1 Ω { t S m : γ - 1 2 t S m : (D m ) -1 : S m }dΩ -λP e (u) (13
) where N refers to the number of layers of laminated shells and λP e (u) is the work of the external forces (λ is a scalar load parameter).

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

N m=1 Ω { t δS : [(γ(u)+ γ)-(D m ) -1 : S m ]+ t S m : δγ(u)}dΩ-λP e (δu) = 0 (14)
By defining the mixed unknown vector U = (u, γ, S m ), we can write Eq.( 14) as the following simple form:

R(U , λ) = L(U ) + Q(U, U ) -λF = 0 ( 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:

< L(U ), δU >= N m=1 Ω ( t δS : [(γ l (u)+ γ)-(D m ) -1 : S m ]+ t S m : γ l (δu))dΩ (16) < Q(U, U ), δU >= N m=1 Ω ( t δS : γ nl (u, u)+ t S m : 2γ nl (u, δu))dΩ (17) < F, δU >= P e (δu) (18)
Thus, the problem is to find the mixed variable U and the loading parameter λ verifying Eq.(15).

Asymptotic numerical algorithm

We propose to solve the nonlinear problem ( 15) by the asymptotic numerical method [START_REF] Büchter | Three-dimensional extension of nonlinear shell formulation based on the enhanced assumed strain concept[END_REF]. 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 [START_REF] Najah | A critical review of asymptotic numerical methods[END_REF].

Starting from an initially known solution (U 0 , λ 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):

U (a) = U 0 + aU 1 + a 2 U 2 + ... + a n U n λ(a) = λ 0 + aλ 1 + a 2 λ 2 + ... + a n λ n ( 19
)
in which U p are mixed unknown vectors, λ p the unknown coefficients at order p and n is the truncation order of the series. The range of validity of the representation [START_REF] Elhage-Hussein | A numerical continuation method based on Padé approximants[END_REF] depends on the path parameter a. By analogy with arclength iterative methods, the control parameter a is chosen as the projection of the displacement increment uu 0 and the load increment λ -λ 0 on the tangent direction (u 1 , λ 1 ) where < ., . > is the Euclidian scalar product:

a =< u -u 0 , u 1 > +(λ -λ 0 )λ 1 (20)
By substituting the series [START_REF] Elhage-Hussein | A numerical continuation method based on Padé approximants[END_REF] 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:

L 0 t (U 1 ) = λ 1 F < u 1 , u 1 > +λ 2 1 = 1 (21) 130 Order 2: L 0 t (U 2 ) = λ 2 F -Q(U 1 , U 1 ) < u 2 , u 1 > +λ 2 λ 1 = 0 (22)
Order p:

   L 0 t (U p ) = λ p F - p-1 r=1 Q(U r , U p-r ) < u p , u 1 > +λ p λ 1 = 0 (23)
The tangent operator L 0 t depends only on the initial solution and is defined by: L 0 t (.) = L(.) + Q(U 0 , .) + Q(., U 0 ) After stress condensation, the linear problems above are solved using the 135 finite element method. For discretization, we use the eight node shell element with reduced integration detailed in [START_REF] Zahrouni | Computing finite rotations of shells by an asymptotic-numerical method[END_REF]. The finite element discretisation is described in Appendix. The discretized form of the problem at order p (p ≥ 2) is given by:

[K 0 t ]{q p } = λ p {F } + {F n p l } t {q p }.{q 1 } + λ p λ 1 = 0 (24)
where [K 0 t ] is the classical tangent stiffness matrix at starting point (U 0 , λ 0 ) , {F } is the external force vector, {q p } and λ p are respectively the discretized form of the displacement u p and the loading parameter at order p. The second membre {F n p l } 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 δ 1 [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF]. Thus, the maximal value of the path parameter a max is given by:

a max = δ 1 ∥ u 1 ∥ ∥ u n ∥ 1 n-1 (25)
with δ 1 a user-defined accuracy parameter and ∥ . ∥ 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.

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 [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF]30], viscous fluid mechanics and plastic structures [START_REF] Zahrouni | Computing finite rotations of shells by an asymptotic-numerical method[END_REF]). 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 [START_REF] Baker | Padé approximants. part 1: Basic theory[END_REF].

We denote P n the representation by rational fractions of the displacement u and the loading parameter λ for a truncation order n:

     P n (u(a)) = u 0 + a D n-2 D n-1 u 1 + a 2 D n-3 D n-1 u 2 + ... + a n-1 1 D n-1 u n-1 P n (λ(a)) = λ 0 + a D n-2 D n-1 λ 1 + a 2 D n-3 D n-1 λ 2 + ... + a n-1 1 D n-1 λ n-1 (26)
where D i (a) are polynomials of degree i with real coefficients d i :

D i (a) = 1 + ad 1 + a 2 d 2 + ... + a i d i (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 [START_REF] Tri | Bifurcation indicator based on meshless and asymptotic numerical methods for nonlinear poisson problems[END_REF]. 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 δ 2 a new accuracy parameter:

∥ P n (u(a max.p )) -P n-1 (u(a max.p )) ∥ ∥ P n (u(a max.p )) -P 0 ∥ = δ 2 (28)
The range of validity (a max.p ) is searched in the interval [a max ; βa max ] by the bisection method using the criterion [START_REF] Büchter | Three-dimensional extension of nonlinear shell formulation based on the enhanced assumed strain concept[END_REF] where β is a parameter defined by the user (β > 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 [START_REF] Elhage-Hussein | A numerical continuation method based on Padé approximants[END_REF].

Bifurcation detection methods

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 [START_REF] Boutyour | Bifurcation points and bifurcated branches by an asymptotic numerical method and Padé approximants[END_REF].

Let ∆µf be a fictitious perturbation force applied to the structure in a deformed state (U , λ). ∆µ represents the intensity of the force and ∆U = (∆u, ∆ γ, ∆S) its associated response. Under these conditions, the perturbed equilibrium equation is written:

L(U + ∆U ) + Q(U + ∆U, U + ∆U ) = λF + ∆µf ( 29
)
Considering the equilibrium equation (Eq.( 15)) and neglecting the quadratic terms, the perturbed problem is written:

L t (∆U ) = ∆µf (30)
where L t (.) = L(.) + 2Q(U, .) is the tangent operator taken at the equilibrium point (U , λ). An additional condition based on the displacement is imposed:

< L 0 t (∆U -∆U 0 ), ∆U 0 >= 0 (31)
where L 0 t is the tangent operator at the starting point (U 0 , λ 0 ). Eqs (30) and (31) will be solved by the asymptotic numerical method and the unknowns ∆U and ∆µ are written as power series:

∆U (a) = ∆U 0 + a∆U 1 + a 2 ∆U 2 + ... + a n ∆U n ∆µ(a) = ∆µ 0 + a∆µ 1 + a 2 ∆µ 2 + ... + a n ∆µ n (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 185 Order 0:

L 0 t (∆U 0 ) = ∆µ 0 f ( 33
)
∆U 0 is computed by imposing ∆µ 0 = 1. Order p ≥ 1:

   L 0 t (∆U p ) = ∆µ p f - p r=1 (Q(U r , ∆U p-r ) + Q(∆U p-r , U r )) < ∆U p , f >= 0 (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:

[K 0 t ]{∆q p } = ∆µ p f + {∆F p } t {∆q p }[K t ]{∆q 0 } = 0 (35)
where [K 0 t ] denotes the tangent stiffness matrix at the starting point, {∆q 0 } and {∆q p } are respectively the nodal displacement vectors at orders 0 and p associated to the perturbation force {f }. The vector {∆F p } depends on the solutions U up to order p and ∆U up to order (p -1). Thus, the system (35) gives:

∆µ p = - < ∆F p , ∆q 0 > < f, ∆q 0 > (36)
Now that we have computed the terms ∆µ p (1 ≤ p ≤ n), we can construct the series development of ∆µ and study its sign along the equilibrium branch in order to detect the bifurcation.

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 [START_REF] Boutyour | Bifurcation points and bifurcated branches by an asymptotic numerical method and Padé approximants[END_REF] [START_REF] Tri | Bifurcation indicator based on meshless and asymptotic numerical methods for nonlinear poisson problems[END_REF]. We note here another interest of the rational representation with respect to the polynomial representation.

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) [START_REF] Boutyour | Vibration of buckled elastic structures with large rotations by an asymptotic numerical method[END_REF]. 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 (θ 1 , θ 2 , θ 3 , t), so the global response is written U (θ 1 , θ 2 , θ 3 ) + V (θ 1 , θ 2 , θ 3 , t). The oscillations of the structure around a static equilibrium state are described by the following equation:

M ( V ) + L(U + V ) + Q(U + V, U + V ) = λF (37) with < M V , δU >= Ω ρ(v 1 δu 1 + v2 δu 2 + v3 δu 3 )dΩ
where 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 ρ is the density.

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

L t (V ) + M ( V ) = 0 ( 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 iwt (39)
By inserting Eq.(39) into Eq.( 38), we obtain the linear eigenvalue problem:

L t (X) -w 2 M (X) = 0 (40)
To study the eigenfrequencies w and their associated modes X in prebuckling and post-buckling ranges, we need to solve the following coupled problem:

L(U ) + Q(U, U ) = λF (41a) L t (X) = w 2 M (X) (41b)
Eq.(41a) corresponds to the static displacement load solution (Eq.( 15)). The obtained displacement U is used to determine the tangent operator L t . 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)

X = X 0 + aX 1 + a 2 X 2 + ... + a n X n w 2 = w 2 0 + aw 1 + a 2 w 2 + ... + a n w n (42)
where (X 0 , w 2 0 ) is the solution of the eigenvalue problem at the starting equilibrium point (U 0 , λ 0 ), X r = (u xr , γ xr , S xr ) is a mixed vector and w r 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:

L 0 t (X 1 ) -w 2 0 M X 1 = w 1 M X 0 -Q(U 1 , X 0 ) -Q(X 0 , U 1 ) (43)
Order r:

L 0 t (X r )-w 2 0 M X r = w r M X 0 + r-1 j=1 w j M X r-j - r j=1 (Q(U j , X r-j )+Q(X r-j , U j ))
(44) Again, we see that all these problems have the same linear operator (L 0 tw 2 0 M ) which must be decomposed only once for all X r (r = 0, 1, ..., n). The projection of (43) and (44) into the vibration mode X 0 gives: Order 1:

w 1 = < X 0 , Q(U 1 , X 0 ) > + < X 0 , Q(X 0 , U 1 ) > < M X 0 , X 0 > (45)
Order r:

w r = r j=1 (< X 0 , Q(U j , X r-j ) > + < X 0 , Q(X r-j , U j ) >) - r-1 j=1 w j < X 0 , M X r-j > < M X 0 , X 0 > (46)
The finite element discretization of the vibration problem gives: Order 1:

     w 1 = t {x 0 }{F 1 } t {x 0 }[M ]{x 0 } ([K 0 t ] -w 2 0 [M ]){x 1 } = w 1 [M ]{X 0 } + {F 1 } (47)
Order r:

             w r = t {x 0 }{Fr}- r-1 j=1 w t j {x 0 }[M ]{xr -j } t {x 0 }[M ]{x 0 } ([K 0 t ] -w 2 0 [M ]){x r } = w r [M ]{X 0 } + r-1 j=1 w j [M ]{x r-j } + {F r } (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).

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 a(mm) R(mm) h(mm) α(rad) 508 2540 6 0.2 steps (35 matrix decompositions) for truncation order n = 15 and accuracy parameter δ = 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 ≈ 30 mm) according to n and δ 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 δ 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 δ = 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 the same path between C and B. Keeping the same representation, another equilibrium path is observed with n = 20 and δ = 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 δ = 10 -10 . Thus, we have six bifurcation points when we try to determine 275 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, the example of n = 15 and δ = 10 -6 , the second step gives us a single pole (p 2 = 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 91 = -0.032) and (p 92 = -0.054 = λ c (C)). This can be attributed to defects in the rational representation because the pole p 91 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 = 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.

Conclusion

In this work, we have presented a bifurcation analysis for multilayer structures using the asymptotic numerical method. Three bifurcation detection Step 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 ≤ n ≤ 30 and the parameter 10 -12 ≤ δ ≤ 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.

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 {δu} and its gradient {θ(u)} are related to nodal displacements q and δq as follows: The Green-Lagrange strain (compatible part) is expressed in the covariant basis as:

γ ij (u) = 1 2 ( ∂u ∂θ i g j + ∂u ∂θ j g i + ∂u ∂θ i ∂u ∂θ j ) ( 50
)
where g i represents the covariant base vectors. γ(u) can be decomposed into

Figure 1 :

 1 Figure 1: Geometric and kinematic description of the shell

  105

  180

Figure 2 :

 2 Figure 2: Geometrical parameters of the cylindrical roof.

Figure 3 :

 3 Figure 3: Load -vertical displacement of the load application point (n = 15 and δ = 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).

Figure 4 : 3 .

 43 Figure 4: The deformations of the structure at steps (a) 20, (b) 50, (c) 90 and (d) 120 in Figure 3.

Figure 5 :

 5 Figure 5: Other equilibrium paths of the structure as a function of representation type, truncation order n and accuracy parameter δ. (a) n = 15, δ = 10 -6 , series representation (b) n = 20, δ = 10 -8 , series representation (c) n = 15, δ = 10 -10 , Padé representation (d) a zoom of Figure 4.c (u ∈ [10; 20] and λ ∈ [-0.3; 0]).

FFigure 6 :FFigure 7 : 2 w 2 0, along the equilibrium branch of Figure 3 ,

 67223 Figure 6: Bifurcation indicator, along the equilibrium branch, as a function of the loading parameter for n = 15 and δ = 10 -6 (Padé). Points A, B, C, D, E and F are the six bifurcation points of the structure.

Figure 8 :

 8 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.

  {u} = [N ]{q} {δu} = [N ]{δq} {θ(u)} = [G]{q}(49) where [N ] is the matrix of shape functions and [G] is their gradient matrix.

Table 1 :

 1 Mechanical properties of each roof layer.

		E(MPa)	µ		G(MPa)
	E 1	3300	µ 12	0.25	G 12	660
	E 2	1100	µ 23	0.25	G 23	440
	E 3	1100	µ 13	0.25	G 13	660

Table 2 :

 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 Figure3.

	Bifurcation point	A	B	C	D	E	F
	critical loading λ c	0.252	0.273	-0.054	-0.184	-0.199	-0.111
	Step	2	4	9	23	26	34
	Log 10 (Res)	-4.683	-4.448	-3.923	-4.383	-4.421	-2.319

Table 3 :

 3 Critical loading, detection step and relative residual Res = ∥Fint-Fext∥

	∥Fext∥	at each

Table 4 :

 4 Poles of the Padé approximants for three couples of n and δ in the vicinity of the bifurcation points.

		2	5	13	30	39	49
	Loading	[0.185;0.274]	[0.263;0.299]	[-0.055;-0.049]	[-0.186;-0.184]	[-0.199;-0.182]	[-0.197;-0.077]
	interval
	λ c	0.252	0.273	-0.054	-0.184	-0.184	-0.111
						-0.199	-0.130
				n=30	δ = 10 -12
	Step	2	6	18	33	37	46
	Loading	[0.206;0.265]	[0.273;0.278]	[-0.056;-0.047]	[-0.189;-0.183]	[-0.199;-0.193]	[-0.203;-0.049]
	interval
	λ c	0.252	0.273	-0.054	-0.184	-0.199	-0.111
			0.274			-0.198

detection step and residual at the end of the step for each bifurcation point, respectively. The residual remains quite small along the equilibrium branch.

a linear part γ l (u) and a non-linear part γ nl (u, u):

with

The virtual strain δγ is written as follows:

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

For the enhanced assumed strain γ = θ 3 β 33 g 3 ⊗ g 3 ,the component β 33 does not require any inter-element continuity. Thus, it is discretized as a bilinear polynomial:

where ξ and η represent the isoparametric coordinates. The unknown parameters α 1 , α 2 , α 3 and α 4 can be eliminated at the elementary level [START_REF] Büchter | Three-dimensional extension of nonlinear shell formulation based on the enhanced assumed strain concept[END_REF]. Thus, the additional deformation is written in the following matrix form:

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

with:

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 n p l } and {R α p } related to solutions at previous orders p -1. These two terms are expressed as follows

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.