Introduction

The increasing use of composite material, in particular for curved structures, has led to the intensive development of appropriate twodimensional shell theories that can accurately describe the response of multilayered anisotropic shells. In fact, the use of three-dimensional modelisation allows us to compute accurate local solutions, but, unfortunately, it is not well suitable due to prohibitive computational costs.

Two main ways are considered to obtain shell FE: (i) the "pure shell model" in which the displacement is associated with the local curvilinear vectors, and strain and stress are deduced using differential geometry [START_REF] Bernadou | Finite Element Methods for Thin Shell Problems[END_REF][START_REF] Dau | A efficient c 1 finite element with continuity requirements for multilayered/sandwich shell structures[END_REF]; (ii) the shell-like solid approach [START_REF] Zienkiewicz | The Finite Element Method[END_REF] where the displacement vector is defined in the global cartesian frame and jacobian matrix transformation is used to express strain and stress with respect to reference frame defined on the middle surface in order to introduce the constitutive law. In this approach, differentiation is simplified and the curvatures are not directly calculated [START_REF] Vidal | An efficient finite shell element for the static response of piezoelectric laminates[END_REF]. It is widely used in commercial software for simplicity reasons. Considering anisotropic structures where 3D phenomena occur, it is also needed to develop enriched theories. According to published research, various approaches based on the Finite Element (FE) method for composite shells have been carried out. In the following, most of the mentioned works refer to the pure shell model. Thus, two families of models [START_REF] Reddy | Mechanics of Laminated Composite Plates and Shells -Theory and Analysis[END_REF] can be identified:

• the Equivalent Single Layer (ESL) Models: It includes the most common theories, namely the classical Shell Theory (CST/Koïter) and First Order Shear Deformation Theory (FSDT/Nagdhi). The reader can refer to Ref. [START_REF] Leissa | Vibration of Shells[END_REF] for a description of the assumptions on the strain to derive different shell models. CST leads to inaccurate results for composites because both transverse shear and normal strains are neglected. In this way, shallow laminated shells are modelized in Refs. [START_REF] Jeyachandrabose | Explicit formulation of two anisotropic, triangular, thin, shallow shell elements[END_REF][START_REF] Rao | A rectangular laminated anisotropic shallow thin shell finite element[END_REF][START_REF] Qatu | Bending analysis of laminated plates and shells by different methods[END_REF]. FSDT is the most popular model due to the possibility to use a C 0 FE, but it needs shear correction factors and transverse normal strain is always neglected, e.g. Ref. [START_REF] Hossain | A finite element formulation for the analysis of laminated composite shells[END_REF]. So, Higher-order Shear Deformation Theories (HSDT) have been developed to overcome these drawbacks. Various theories are proposed based on 7 ( [START_REF] Sgambitterra | A robust three-node shell element for laminated composites with matrix damage[END_REF][START_REF] Balah | Finite element formulation of a third order laminated finite rotation shell element[END_REF]) or 9 parameters [START_REF] Kant | Estimation of interlaminar stresses in fibre reinforced composite cylindrical shells[END_REF]. Wide range of theories based on the Carrera's Unified Formulation are also addressed in Ref. [START_REF] Cinefra | Shell finite elements with different through-the-thickness kinematics for the linear analysis of cylindrical multilayered structures[END_REF]. In the ESL context, a simple way to improve the estimation of the mechanical quantities consists in adding zig-zag functions (Murakami [START_REF] Murakami | Laminated composite plate theory with improved in-plane responses[END_REF][START_REF] Carrera | Developments, ideas and evaluations based upon the reissner's mixed theorem in the modeling of multilayered plates and shells[END_REF]) in the displacement to introduce the slope discontinuity at the interface between two adjacent layers. It allows to describe the so-called zig-zag effect. It has been carried out by Brank [START_REF] Brank | On composite shell models with a piecewise linear warping function[END_REF] with a 7 parameters model. See also the work of Bhaskar [START_REF] Bhaskar | A higher-order theory for bending analysis of laminated shells of revolution[END_REF] based on a HSDT approach.

stresses can be written over each layer. Thus, the number of unknowns depends on the number of layers. A quadratic triangle element based on a constant shear angle is considered in Ref. [START_REF] Seide | Triangular finite element for analysis of thick laminated shells[END_REF], but a shear correction factor is needed. A three-dimensional shell element is proposed in Ref. [START_REF] Klinkel | A continuum based three dimensional shell element for laminated structures[END_REF]. A LW triangle FE is developed in Ref. [START_REF] Botello | A layer-wise triangle for analysis of laminated composite plates and shells[END_REF] with a condensation technique at the pre-processing level. Ref. [START_REF] Liu | Free vibration analysis of laminated composite shell structures using hybrid strain based layerwise finite elements[END_REF]. deals with a hybrid strain flat triangular FE based on the Hellinger-Reissner variational principle. Note that the trans-verse normal and shear stresses are only taken into account in Ref. [START_REF] Basar | Interlaminar stress analysis of composites: layer-wise shell finite elements including transverse strains[END_REF] where a four-node isoparametric assumed strain is considered. We can mention the eight-node 3D hybrid-EAS solid shell element based on the Hu-Washizu variational principle in Ref. [START_REF] Vu-Quoc | Efficient hybrid-eas solid element for accurate stress prediction in thick laminated beams, plates, and shells[END_REF]. See also the previously mentioned work [START_REF] Cinefra | Shell finite elements with different through-the-thickness kinematics for the linear analysis of cylindrical multilayered structures[END_REF].

As an alternative, refined models have been developed in order to improve the accuracy of ESL models avoiding the additional computational cost of LW approach. Based on physical considerations and after some algebraic transformations, the number of unknowns becomes independent of the number of layers. We can mention the work of Yasin [START_REF] Yasin | An efficient layerwise finite element for shallow composite and sandwich shells[END_REF] including 5 parameters dedicated to shallow shells. Shariyat [START_REF] Shariyat | Non-linear dynamic thermo-mechanical buckling analysis of the imperfect laminated and sandwich cylindrical shells based on a global-local theory inherently suitable for non-linear analyses[END_REF] has also developed a so-called zig-zag model including 15 parameters. It should be also mentioned the work of Dau [START_REF] Dau | A efficient c 1 finite element with continuity requirements for multilayered/sandwich shell structures[END_REF] where a C 1 triangular six-node FE (Argyris-Ganev) based on the Sinus model is considered involving 5 parameters. The approach ensures the continuity conditions of the transverse shear stresses at the interfaces between two adjacent layers.

For the present topics, it should be noted that the mentioned works are based on the Finite Element method for linear elasticity problem in mechanics and applied to laminated composites, knowing that many other approaches (meshless, analytical, semi-analytical …) are involved in open literature. Furthermore, the fundamental subject about the shear and membrane locking of shell is not addressed here. So, this above literature deals with only some aspects of the broad research activity about composite shells. An extensive assessment of different approaches for both various theories and/or finite element applications can be found in Refs. [START_REF] Kapania | A review on the analysis of laminated shells[END_REF][START_REF] Noor | Assessment of computational models for multilayered composite shells[END_REF][START_REF] Gilewski | A survey of finite element models for the analysis of moderately thick shells[END_REF][START_REF] Yang | A survey of recent shell finite element[END_REF][START_REF] Carrera | Theories and finite elements for multilayered, anisotropic, composite plates and shells[END_REF][START_REF] Reddy | Shear deformation plate and shell theories: from stavsky to present[END_REF][START_REF] Hohe | Advances in the structural modeling of elastic sandwich panels[END_REF][START_REF] Qatu | Review of recent literature on static analyses of composite shells: 2000-2010[END_REF].

The main goal of this work consists in assessing two new degenerated shell approaches to model any laminated composite structures. For the transverse coordinate, a fourth-order z-expansion in each layer is used as it is well suited for this type of structures [START_REF] Vidal | Proper generalized decomposition and layer-wise approach for the modeling of composite plate structures[END_REF]. The first approach is based on the classical Layer-Wise model [START_REF] Reddy | Mechanics of Laminated Composite Plates and Shells -Theory and Analysis[END_REF][START_REF] Carrera | Developments, ideas and evaluations based upon the reissner's mixed theorem in the modeling of multilayered plates and shells[END_REF] where the number of unknowns per node depends on the number of layers. The linear problem is directly solved. The originality of this approach is the conjunction of this classical LW model with the degenerated shell model. The second one implies the so-called Proper Generalized Decomposition (PGD) [START_REF] Chinesta | An overview of the proper generalized decomposition with applications in computational rheology[END_REF][START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF][START_REF] Allix | A new multi-solution approach suitable for structural identification problems[END_REF] as it has shown interesting features in the reduction model framework. In particular, it has been applied for composite shells in Refs. [START_REF] Pruliere | 3D simulation of laminated shell structures using the proper generalized decomposition[END_REF][START_REF] Bognet | Separated representations of 3D elastic solutions in shell geometries[END_REF][START_REF] Vidal | Shell finite element based on the proper generalized decomposition for the modeling of cylindrical composite structures[END_REF][START_REF] Vidal | Multi-resolution strategies for the modeling of composite shell structures based on the variable separation method[END_REF]. However, these studies are either dedicated to cylindrical shell geometry or not well suited to model composite structures. The present approach is based on the separation representation where the displacements are written under the form of a sum of products of bidimensional polynomials of (x 1 , x 2 ) and unidimensional polynomials of z. The deduced non-linear problem is solved using a fixed point strategy in which the resolution of two linear problems is carried out alternatively. This process yields to a 2D and a 1D problems in which the number of unknowns is smaller than the previous Layerwise approach.

This article is organized as follows. First, the geometry of the shelllike solid is described including the reference frames. The general problem to be solved is precised. Then, the two involved approches, namely the LW model and PGD method, are described in the particular framework of the degenerated shell model. For the latter, the particular assumption on the displacements yields a non-linear problem solved by a fixed point strategy. The FE discretization is also described and finally, numerical tests are performed. The present approaches are assessed on different shell configurations. Deep and shallow shells, different slenderness ratios and stacking sequences are considered for various boundary conditions. The limitations and advantages of the models are precised and comparison with reference solutions are given.

Description of the shell

This preliminary section is dedicated to the geometric description of the shell-like solid and to the different reference frames that will be used to construct the finite shell element.

Description of the geometry

Let us consider a shell

 = Ω × [-e 2 , e 2 ]
where e is the constant thickness of the shell and Ω is the middle surface. The description of the geometry of the shell is based on the Cartesian coordinates of the nodes and on the finite element approximation over the elementary domain. For this, an eight-node quadrilateral finite element is used. Fig. 1 illustrates the employed finite element approximation and shows the different reference frames to describe the geometry and the mechanics of the shell-like body. These different bases are precised in Table 1. The global coordinates of any arbitrary point in the elementary domain can be expressed in terms of the reduced (curvilinear) coordinates 𝜉, 𝜂 and the rectilinear normal coordinate z:

𝚽(𝜉, 𝜂, z) = ⎡ ⎢ ⎢ ⎢ ⎣ X 1 (𝜉, 𝜂, z) X 2 (𝜉, 𝜂, z) X 3 (𝜉, 𝜂, z) ⎤ ⎥ ⎥ ⎥ ⎦ = 8 ∑ i=1 N q i (𝜉, 𝜂) ⎡ ⎢ ⎢ ⎢ ⎣ X 1 X 2 X 3 ⎤ ⎥ ⎥ ⎥ ⎦ i + z 8 ∑ i=1 N q i (𝜉, 𝜂) ⎡ ⎢ ⎢ ⎢ ⎣ t 31 t 32 t 33 ⎤ ⎥ ⎥ ⎥ ⎦ i (1)
where N q i (𝜉, 𝜂) are the classical Serendipity interpolation functions.

The thickness of the shell is described in terms of the unit vector

t 3 = [ t 31 t 32 t 33
] T normal to the middle surface Ω. We prescribe t 3 = a 3 , i.e., the normal direction of the local curvilinear reference frame on Ω coincides with the normal direction of the local tangent plane to Ω (see Fig. 1). The covariant in-plane base vectors a 𝛼 are usually obtained from the map 𝚽 introduced in Eq. ( 1) to define the shell mid-surface Ω:

and the normal vector t 3 = a 3 is finally obtained from the perpendicularity condition:

t 3 = a 3 = a 1 ∧ a 2 ‖a 1 ∧ a 2 ‖ (3)
The construction of the local base vectors t 1 and t 2 follows the procedure suggested by Ref. [START_REF] Zienkiewicz | The Finite Element Method[END_REF].

The change of bases

The displacement field vector is expressed with respect to the global reference frame e i as

u(X 1 , X 2 , X 3 ) = 3 ∑ i=1 u G i (X 1 , X 2 , X 3 )e i (4)
where the superscript G indicates that the components are taken in the global reference frame. The displacement vector u = [u 1 u 2 u 3 ] T is expressed in the local orthonormal basis (L) by

u L = [ T LG ] u G with [ T LG ] = ⎡ ⎢ ⎢ ⎢ ⎣
In the same way, the following expression can be obtained:

u G = [ T GL ] u L with [ T GL ] = [ T LG ] -1 = [ T LG ] T = ⎡ ⎢ ⎢ ⎢ ⎣
⎤ ⎥ ⎥ ⎥ ⎦ (6)
The displacement vector u R defined in the reduced local basis is constructed as

u R = [ T RG ] u G with [ T RG ] = ⎡ ⎢ ⎢ ⎢ ⎣
⎤ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎣ a 1 • e 1 a 1 • e 2 a 1 • e 3 a 2 • e 1 a 2 • e 2 a 2 • e 3 a 3 • e 1 a 3 • e 2 a 3 • e 3 ⎤ ⎥ ⎥ ⎥ ⎦ (7)
Furthermore, the inverse relation holds:

u G = [ T GR ] u R with [ T GR ] = [ T RG ] -1 (8)
Similar relations can be established between the local orthonormal L and reduced R reference frames by substituting G by L.

Reference problem description

The weak form of the boundary value problem

The shell is submitted to a surface force density t defined over a subset 𝜕 F of the boundary. We assume that a prescribed displacement u = u d is imposed on Γ D = 𝜕 -𝜕 F . The mechanical problem is based on the variational principle which is given by: find u(M) ∈ U (space of admissible displacements) such that

a(u, 𝛿u) = b(𝛿u), ∀𝛿u ∈ 𝛿U (9) with a(u, 𝛿u) = ∭  𝛿𝜺 L T 𝝈 L d b(𝛿u) = ∬ 𝜕 F 𝛿u • t d𝜕
where U is the space of admissible displacements, i.e. U = {u ∈ (H 1 ())

3 ∕u = u d on Γ D }. We have also 𝛿U = {u ∈ (H 1 ()) 3 ∕u = 0 on Γ D }.
For the present work, 𝜕 F is considered as the top or bottom surface of the shell, that is z = z F with z F = ±e∕2. For sake of clarity, the body forces are neglected in the developments.

The constitutive relation

The stress tensor is obtained from the strain tensor using the constitutive equations. All these tensors will be referred to the local basis vectors associated with the middle surface of the shell. In case of laminated shells composed of NC perfectly bonded orthotropic plies, the threedimensional constitutive law of the kth layer in the local orthonormal basis is expressed by: 𝝈 L(k) = C (k) 𝜺 L(k) [START_REF] Hossain | A finite element formulation for the analysis of laminated composite shells[END_REF] For convenience reasons, stresses 𝝈 and strains 𝜺 are split into three groups: [START_REF] Sgambitterra | A robust three-node shell element for laminated composites with matrix damage[END_REF] where the subscripts p, n and s denote in-plane, transverse normal and shear values, respectively.

𝝈 T p = [𝜎 11 𝜎 22 𝜏 12 ], 𝝈 T s = [𝜏 13 𝜏 23 ], 𝝈 n = 𝜎 33 𝜺 T p = [𝜀 11 𝜀 22 𝛾 12 ], 𝜺 T s = [𝛾 13 𝛾 23 ], 𝜺 n = 𝜀 33
Using the separation between the different components described above, the three dimensional constitutive law of the kth layer is given by:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝝈 (k) p = Q (k) pp 𝜺 (k) p + Q (k) pn 𝜺 (k) n 𝝈 (k) n = Q (k) np 𝜺 (k) p + Q (k) 33 𝜺 (k) n 𝝈 (k) s = Q (k) s 𝜺 (k) s ( 12
)
where

Q (k) pp = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ Q (k) 11 Q (k) 12 Q (k) 16 Q (k) 12 Q (k) 22 Q (k) 26 Q (k) 16 Q (k) 26 Q (k) 66 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ Q (k) pn = Q (k) np T = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ Q (k) 13 Q (k) 23 Q (k) 36 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ Q (k) s = ⎡ ⎢ ⎢ ⎣ Q (k) 55 Q (k) 45 Q (k) 45 Q (k) 44 ⎤ ⎥ ⎥ ⎦ (13)
where Q (k) ij are the three-dimensional stiffness coefficients of the layer (k).

The strain

Thecompatiblestrainfieldisobtained from the linear straindisplacement relations, which in the global reference frame read

𝜀 G ij ( u G ) = 1 2 ( u G ,i + u G ,j ) (14)
where i, j = X 1 , X 2 , X 3 . The components of the strain tensor can be calculated with respect to the three bases G, L and R:

𝜀(u G ) = 𝜀 G ij (e i ⊗ e j ) = 𝜀 L ij (t i ⊗ t j ) = 𝜀 R ij (a i ⊗ a j ) = 𝜀 R ij (a i ⊗ a j ) (15)
It must be noted that covariant and contravariant components are the same in an orthonormal basis, for example e i and t i . This is not the case when the strain components are expressed in the local reduced reference frame a i . Therefore, we can express the strain tensor in the reduced basis in either covariant or contravariant components.

In this work, strain components are first calculated with respect to the local reduced base vectors, using FE approximations. Then, using tensorial transformation, they are expressed with respect to the local cartesian base vectors, in order to take into account the constitutive equations. Therefore, the local cartesian base vectors must be chosen as reference frame and introduced in the variational formulation.

For this purpose, we have

𝜖 L ij = 𝜖 R kl (a k .t i ) (a l .t j ) (16)
This tensorial transformation can be written under matrix product form, for the following components [START_REF] Sgambitterra | A robust three-node shell element for laminated composites with matrix damage[END_REF][START_REF] Liu | Free vibration analysis of laminated composite shell structures using hybrid strain based layerwise finite elements[END_REF][START_REF] Balah | Finite element formulation of a third order laminated finite rotation shell element[END_REF], using the transformation matrix denoted

[ TT LR ]
. The transformation between local and reduced for the third index (3 or z) is 1, because same unit vector is used.

The expression of this transformation is obtained from:

[ T LR ] = [ T RL ] -1 with [ T RL ] = [ a 1 .t 1 a 1 .t 2 a 2 .t 1 a 2 .t 2 ] (17)
and using tensorial matrix of order two which have the following expression:

[ TT LR ] = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ T LR (1, 1) 2 T LR (1, 2) 2 T LR (1, 1) T LR (1, 2) T LR (2, 1) 2 T LR (2, 2) 2 T LR (2, 1) T LR (2, 2) 2 T LR (1, 1) T LR (2, 1) 2 T LR (1, 2) T LR (2, 2) T LR (1, 1) T LR (2, 2) + T LR (1, 2) T LR (2, 1) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (18)
Finally, we have to calculate the strain in the reduced local basis. In this local basis, coordinates denoted [START_REF] Bernadou | Finite Element Methods for Thin Shell Problems[END_REF][START_REF] Dau | A efficient c 1 finite element with continuity requirements for multilayered/sandwich shell structures[END_REF] are in fact (𝜉, 𝜂) due to the use of FE approximation.

We have

𝜖 R ij = 1 2 ( u L ,i . a j + u L ,j . a i ) (19)
and

u L ,i . a j = u L k,i (t k . a j ) + u L k (t k,j . a i ) (20)
In this last expression, the second term needs to use differential geometry of shell in order to express derivatives of the base vectors. In the following, it is neglected in order to simplify the numerical implementation.

First approach: classical LW model with degenerated shell approach

In this section, the classical LW approach is briefly recalled as it is widely developed in open literature [START_REF] Reddy | Mechanics of Laminated Composite Plates and Shells -Theory and Analysis[END_REF][START_REF] Carrera | Developments, ideas and evaluations based upon the reissner's mixed theorem in the modeling of multilayered plates and shells[END_REF]. Hereafter, the developments are focused on the framework of our degenerated shell approach.

The fourth-order expansion of the displacement in the kth layer is expressed as

u L(k) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ u (k) 1 (x 1 , x 2 , z) u (k) 2 (x 1 , x 2 , z) u (k) 3 (x 1 , x 2 , z) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = 5 ∑ j=1 f j LW (z) ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ u (k) j (x 1 , x 2 ) v (k) j (x 1 , x 2 ) w (k) j (x 1 , x 2 ) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (21)
where

f j LW (z) are Lagrangian interpolation functions. u (k) j , v (k) j , w (k) j are defined in Ω.
Using this expression to compute the strains in conjunction with the changes of basis described in Section 3.3 and introducing them in Eq. ( 9) with a FE approximation, it drives to a classical linear system to be solved as [START_REF] Liu | Free vibration analysis of laminated composite shell structures using hybrid strain based layerwise finite elements[END_REF] where K LW and  LW are the rigidity matrix and the force vector of the LW approach, q LW contains all the unknown degrees of freedom. The expression of the classical elementary matrices before the assembly process is not given for brevity reason.

K LW q LW =  LW
At each node, the number of unknowns depends on the number of layers, N LW dofnode = 3.(4.N z + 1), where N z is the number of numerical layers. Thus, the size of the system increases with N z as

N LW = 3.(4.N z + 1).(3.N x .N y + 2(N x + N y ) + 1)
. N x xN y are the number of elements in the 𝜉 and 𝜂 directions, respectively.

Second approach: application of the Proper Generalized Decomposition to degenerated shell model

In this section, the application of the PGD is developed for the present degenerated shell analysis. This work is an extension of the previous studies on shell structures ( [START_REF] Vidal | Shell finite element based on the proper generalized decomposition for the modeling of cylindrical composite structures[END_REF][START_REF] Vidal | Multi-resolution strategies for the modeling of composite shell structures based on the variable separation method[END_REF]) to take into account any geometrical composite shells.

The displacement

The displacement solution is constructed as the sum of N products of functions of in-plane coordinates and transverse coordinate (N ∈ ℕ is the order of the representation) in the global basis:

u G = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ u 1 (x 1 , x 2 , z) u 2 (x 1 , x 2 , z) u 3 (x 1 , x 2 , z) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = N ∑ j=1 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ f j 1 (z) v j 1 (x 1 , x 2 ) f j 2 (z) v j 2 (x 1 , x 2 ) f j 3 (z) v j 3 (x 1 , x 2 ) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (23)
where

(f j 1 , f j 2 , f j 3 ) are defined in Ω z and (v j 1 , v j 2 , v j
3 ) are defined in Ω. Ω and Ω z are bi-dimensional domain associated with the mid-surface of the shell, and the unidimensional domain associated with the normal fiber respectively.

In this paper, a classical eight-node FE approximation is used in Ω and a LW description is chosen in Ω z as it is particularly suitable for the modeling of composite structure.

The problem to be solved

Now, we briefly recall the process to build the different terms (couples) involving in Eq. ( 23). Classically, the resolution of Eq. ( 9) is based on a greedy algorithm. If we assume that the first m functions have been already computed, the trial function for the iteration m + 1 is written as

u m+1 = u m + ⎡ ⎢ ⎢ ⎢ ⎣ f 1 v 1 f 2 v 2 f 3 v 3 ⎤ ⎥ ⎥ ⎥ ⎦ = u m + f ∘ v (24)
where the "•" operator stands for Hadamard's element-wise product, (v 1 , v 2 , v 3 ) and (f 1 , f 2 , f 3 ) are the functions to be computed and u m is the associated known set at iteration m defined by

u m = m ∑ i=1 ⎡ ⎢ ⎢ ⎢ ⎣ f i 1 v i 1 f i 2 v i 2 f i 3 v i 3 ⎤ ⎥ ⎥ ⎥ ⎦ (25)
The test functions are

𝛿(f ∘ v) = 𝛿f ∘ v + 𝛿v ∘ f (26) with v = ⎡ ⎢ ⎢ ⎢ ⎣ v 1 v 2 v 3 ⎤ ⎥ ⎥ ⎥ ⎦ f = ⎡ ⎢ ⎢ ⎢ ⎣ f 1 f 2 f 3 ⎤ ⎥ ⎥ ⎥ ⎦ (27)
Introducing the test function defined by Eq. ( 26) and the trial function defined by Eq. ( 24) into the weak form Eq. ( 9), the two following equations to be solved can be deduced:

• for the test function 𝛿f a(v ∘ f, v ∘ 𝛿f) = b(v ∘ 𝛿f) -a(u m , v ∘ 𝛿f) ∀𝛿f (28)
• for the test function 𝛿v

a(f ∘ v, f ∘ 𝛿v) = b(f ∘ 𝛿v) -a(u m , f ∘ 𝛿v) ∀𝛿v (29)
The bilinear form a can be written under the following form using Eq. ( 12)

a(u, 𝛿u) = ∫  𝛿𝜺 L p T ( Q (k) pp 𝜺 L p + Q (k) pn 𝜺 L n ) + 𝛿𝜺 L n T ( Q (k) np 𝜺 L p + Q (k) 33 𝜺 L n ) + 𝛿𝜺 L s T Q (k) s 𝜺 L s d (30)
This expression has to be written in the suitable form depending on the two problems Eq. ( 28) and Eq. ( 29). This coupled non-linear problem is solved using a classical strategy based on a fixed point method. For each problem, only unknown 1D or 2D function has to be found. So, the approach leads to the process given in Algorithm 1. The fixed point algorithm is stopped when the distance between two consecutive terms are sufficiently small (Cf. [START_REF] Vidal | Proper generalized decomposition and layer-wise approach for the modeling of composite plate structures[END_REF]).

Algorithm 1 Classical algorithm

for m = 0 to N max do Initialize ṽ(0) for k = 1 to k max do Compute f(k) from Eq. ( 28) (linear equation on Ω z ), ṽ(k-1) being known Compute ṽ(k) from Eq. ( 29) (linear equation on Ω), f(k)

being known Check for convergence end for

Set f m+1 = f(k) , v m+1 = ṽ(k) Set u m+1 = u m + f m+1 •v m+1
Check for convergence end for Fig. 2. discretization of the displacements for the PGD method.

Finite element discretization

In this section, the shell finite element approximation is built. For the present approach, the computation of the discretized displacements relies on 2 meshes, as shown in Fig. 2. The first one is 1D with 3 unknowns par node related to the discretization of the z-

functions f 1 (z), f 2 (z), f 3 (z).
It is possible to choose different interpolation orders. That will determine the number of nodes of the 1D mesh. For instance, a fourth-order expansion with NC layers implies 4.NC + 1 nodes (considering one element per physical layer). The second one is 2D with 3 unknowns per node related to the discretization of the in-plane functions

v 1 (x 1 , x 2 ), v 2 (x 1 , x 2 ), v 3 (x 1 , x 2 )
. Thus, considering any points on the middle plane of the shell structure, the computations of the displacements through the thickness depend on the same N z-functions. The unknown functions of z are global for the whole shell.

In this work, a discrete representation of the functions (v, f) is introduced, using (i) a classical eight-node shell finite element approximation in Ω, (ii) Lagrange interpolations in Ω z . The elementary vector of degrees of freedom (dof) associated with one element Ω e of the mesh in Ω is denoted q v e . The vector of dofs associated with the expansion in Ω z is denoted q f . The displacement field is determined from the values of q v e and q f e by

v e = N 𝜉 q v e , f e = N z q f e ( 31
)
The matrices N 𝜉 , N z contain the interpolation functions.

Finite element problem to be solved on Ω

For the sake of simplicity, the function f(k) which is assumed to be known, will be denoted f, and the function ṽ(k) to be computed will be denoted v. Using the particular form of the displacement Eq. ( 23) in Eq. ( 19) and the FE approximation for v and f, the reduced in-plane, normal and shear strains can be deduced under the following form:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝜺 R p ( f v) = 𝚺 p z ( f ) BaseTA p q v e 𝜺 R n ( f v) = 𝚺 n z ( f ) BaseTA n q v e 𝜺 R s ( f v) = [ 𝚺 s1 z ( f ) BaseTA s1 + 𝚺 s3 z ( f ) BaseTA s3 ] q v e ( 32
)
where BaseTA p , BaseTA n and BaseTA s1 , BaseTA s3 contain the shape functions, their derivatives, the scalar products t i • a i and t i • e i . Their expressions and those of 𝚺 p z , 𝚺 n z , 𝚺 s1 z , 𝚺 s3 z are given in A.1.

Using Eq. ( 32) and Eq. ( 16) in the expression Eq. ( 30), the discretization of the bilinear form a in the weak formulation given in Eq. ( 29) can be deduced:

a(u, 𝛿u) = ∑ e ∫ Ω e ∫ Ω z 𝛿q v e T { BaseTA p T 𝚺 p z ( f ) T TT T LR Q (k) pp TT LR 𝚺 p z ( f ) BaseTA p + BaseTA p T 𝚺 p z ( f ) T TT T LR Q (k) pn 𝚺 n z ( f ) BaseTA n + BaseTA n T 𝚺 n z ( f ) T Q (k) np TT LR 𝚺 p z ( f ) BaseTA p + BaseTA n T 𝚺 n z ( f ) T Q (k) 33 𝚺 n z ( f ) BaseTA n + [ 𝚺 s1 z ( f ) BaseTA s1 + 𝚺 s3 z ( f ) BaseTA s3 ] T T T LR Q (k) ss T LR [ 𝚺 s1 z ( f ) BaseTA s1 + 𝚺 s3 z ( f ) BaseTA s3
] } q v e dΩ z dΩ e [START_REF] Hohe | Advances in the structural modeling of elastic sandwich panels[END_REF] To split the integration with respect to the midsurface and the zcoordinate, it is needed to write the three matrix products TT LR 𝚺 p z ( f ),

T LR 𝚺 s1 z ( f ) and T LR 𝚺 s3 z ( f ) in Eq. ( 33
) under a specific form M 1 ( f )M 2 .
Thus, we introduce these new expressions as

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ TT LR 𝚺 p z ( f ) = 𝚺 p new ( f ) MM p LR T LR 𝚺 s1 z ( f ) =𝚺 s1 new ( f ) M s1 LR T LR 𝚺 s3 z ( f ) =𝚺 s3 new ( f ) M s3 LR ( 34
)
Only the expressions of 𝚺

+ BaseTA s1 T M s1 LR T k s11 z ( f ) M s1 LR BaseTA s1 + BaseTA s1 T M s1 LR T k s13 z ( f ) M s3 LR BaseTA s3 + BaseTA s3 T M s3 LR T k s31 z ( f ) M s1 LR BaseTA s1 + BaseTA s3 T M s3 LR T k s33 z ( f ) M s3 LR BaseTA s3 } q v e dΩ e ( 35
)
with

k pp z ( f ) = ∫ Ω z 𝚺 p new ( f ) T Q (k) pp 𝚺 p new ( f )dz k pn z ( f ) = ∫ Ω z 𝚺 p new ( f ) T Q (k) pn 𝚺 n z ( f )dz k np z ( f ) = k pn z ( f ) T k 33 z ( f ) = ∫ Ω z 𝚺 n z ( f ) T Q (k) 33 𝚺 n z ( f )dz k s11 z ( f ) = ∫ Ω z 𝚺 s1 new ( f ) T Q (k) ss 𝚺 s1 new ( f )dz k s13 z ( f ) = ∫ Ω z 𝚺 s1 new ( f ) T Q (k) ss 𝚺 s3 new ( f )dz k s33 z ( f ) = ∫ Ω z 𝚺 s3 new ( f ) T Q (k) ss 𝚺 s3 new ( f )dz k s31 z ( f ) = k s13 z ( f ) T (36)
Note that the calculation of Eq. ( 36) is performed using an analytical integration. Using Eq. ( 35), the variational problem defined on Ω from Eq. ( 29) leads to the linear system

K z ( f )q v =  v ( f , u n ) (37)
where K z ( f ) is deduced from Eq. ( 35) and q v is the vector of the nodal displacements in the global basis. The right hand side of Eq. ( 29) is not given for brevity reason. The same discretization process is used to deduce it.

Finite element problem to be solved on Ω z

For the sake of simplicity, the function ṽ(k-1) which is assumed to be known, will be denoted ṽ, and the function f(k) to be computed will be denoted f.

Using the particular form of the displacement Eq. ( 23) in Eq. ( 19) and the FE approximation for v and f, the three parts of the strains can be conveniently written under the following form:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝜺 R p (ṽ f ) =𝚺 p xy (ṽ) BaseNz q f e 𝜺 R n (ṽ f ) =𝚺 n xy (ṽ) BaseBz q f e 𝜺 R s (ṽ f ) =𝚺 s1 xy (ṽ) BaseNz q f e + 𝚺 s3 xy (ṽ) BaseBz q f e ( 38
)
where BaseNz and BaseBz contain the shape functions and their derivatives. Their expressions and those of 𝚺 p xy (ṽ), 𝚺 n xy (ṽ), 𝚺 s1 xy (ṽ), 𝚺 s3 xy (ṽ)

are given in B.

Starting from the problem to be solved Eq. ( 28), Eq. ( 38) and Eq. ( 16) can be introduced to deduce the expression of a:

a(v ∘ f, v ∘ 𝛿f) = ∑ e ∫ Ω z e ∫ Ω 𝛿q f e T { BaseNz T 𝚺 p xy (ṽ) T TT T LR Q (k)
pp TT LR 𝚺 p xy (ṽ) BaseNz + BaseNz T 𝚺 p xy (ṽ) T TT T LR Q (k) pn 𝚺 n xy (ṽ) BaseBz

+ BaseBz T 𝚺 n xy (ṽ) T Q (k) np TT LR 𝚺 p xy (ṽ) BaseNz + BaseBz T 𝚺 n xy (ṽ) T Q (k) 33 𝚺 n xy (ṽ) BaseBz + [ 𝚺 s1 xy (ṽ) BaseNz + 𝚺 s3 xy (ṽ) BaseBz ] T T T LR Q (k) ss T LR [ 𝚺 s1 xy (ṽ) BaseNz + 𝚺 s3 xy (ṽ) BaseBz ] } q f e dΩdΩ z e ( 39
)
By integrating with respect to the in-plane coordinate, Eq. ( 39) becomes

a(v ∘ f, v ∘ 𝛿f) = ∑ e 𝛿q f e T ∫ Ω z e { BaseNz T k pp xy (k) (ṽ) BaseNz + BaseNz T k pn xy (k) (ṽ) BaseBz + BaseBz T k np xy (k) (ṽ) BaseNz + BaseBz T k 33 xy (k) (ṽ) BaseBz + BaseNz T k s11 xy (k) (ṽ) BaseNz + BaseNz T k s13 xy (k) (ṽ) BaseBz + BaseBz T k s31 xy (k) (ṽ) BaseNz + BaseBz T k s33 xy (k) (ṽ) BaseBz } dΩ z e q f e (40)
with

k pp xy (k) (ṽ) = ∫ Ω 𝚺 p xy (ṽ) T TT T LR Q (k) pp TT LR 𝚺 p xy ( f )dΩ k pn xy (k) (ṽ) = ∫ Ω 𝚺 p xy (ṽ) T TT T LR Q (k) pn 𝚺 n xy (ṽ)dΩ k np xy (k) (ṽ) = k pn z (k) ( f ) T k 33 xy (k) (ṽ) = ∫ Ω 𝚺 n xy (ṽ) T Q (k) 33 𝚺 n xy (ṽ)dΩ k s11 xy (k) (ṽ) = ∫ Ω 𝚺 s1 xy (ṽ) T T T LR Q (k)
ss T LR 𝚺 s1 xy (ṽ)dΩ

k s13 xy (k) (ṽ) = ∫ Ω 𝚺 s1 xy (ṽ) T T T LR Q (k)
ss T LR 𝚺 s3 xy (ṽ)dΩ

k s33 xy (k) (ṽ) = ∫ Ω 𝚺 s3 xy (ṽ) T T T LR Q (k) ss T LR 𝚺 s3 xy (ṽ)dΩ k s31 xy (k) (ṽ) = k s13 xy (k) (ṽ) T (41)
After an assembling procedure, Eq. ( 40) can be reported in Eq. ( 28) to deduce a linear system to be solved:

K xy (ṽ)q f =  f (ṽ, u n ) (42)
where K xy (ṽ) is deduced from Eq. ( 40) and q f is the vector of the nodal unknowns associated to f. The right hand side of Eq. ( 28) is not given for brevity reason. The same discretization process is used to obtain it.

Numerical results

In this section, an eight-node quadrilateral FE based on the Serendipity interpolation functions is used for the unknowns depending on the in-plane coordinates. The geometry of the shell is approximated by this classical FE in the parametric space. A Gaussian numerical inte-gration with 3 × 3 points is used to evaluate the elementary matrices.

We have developed two different degenerated shell approaches assuming a fourth-order z-expansion through each layer: (i) the first one, denoted LD4-S, is based on the LD4 model, as developed in the systematic work of Carrera and his "Carrera's Unified Formulation" (CUF), see Refs. [START_REF] Carrera | Developments, ideas and evaluations based upon the reissner's mixed theorem in the modeling of multilayered plates and shells[END_REF][START_REF] Carrera | Historical review of zig-zag theories for multilayered plates and shells[END_REF][START_REF] Ottavio | Considerations on higher-order finite elements for multilayered plates based on a unified formulation[END_REF]. In the latter, 12NC + 3 unknown functions per node are used in this kinematic; (ii) the second one, denoted VS-LD4-S, is based on the variable separation approach (PGD), the number of unknowns is reduced to three per node for the 1D and 2D models. The results are compared with an approach based on an exact geometry description using LD4 model. It is denoted LD4-M. In the subsequent examples, this model can be considered as a reference solution. As the same assumptions on the shell geometry are used for both LD4-S and VS-LD4-S models, the comparison between these two models allows us to evaluate only the error associated to the variable separation part.

Several static tests are presented in order to validate our approaches and evaluate their efficiencies. For each of them, a preliminary convergence study has been carried out to fix the suitable mesh. Only the first one is presented for conciseness reason. Different boundary conditions, geometries and number of layers are addressed to assess the present approaches and highlight the possibilities and limitations.

For the numerical simulations related to the PGD model, two simultaneous couples are first computed and the following ones are built one by one. This strategy is needed to achieve the convergence of the method.

Simply-supported cylindrical shell

The first test case is based on a configuration proposed by Ren [START_REF] Ren | Exact solutions for laminated cylindrical shells in cylindrical bending[END_REF]. It concerns a semi-infinite simply-supported cylindrical shell submitted to a constant pressure. It is described below: Fig. 3). boundary conditions: simply-supported shell along its straight edges, subjected to a constant pressure along the curvature:

convergence study -S

= 4 -𝜑 = 𝜋∕2-3 layers [0 • ∕90 • ∕0 • ].
q(𝜉 1 ) = q 0 material properties: E L = 25 GPa , E T = 1 GPa , G LT = 0.2 GPa, G TT = 0.5 GPa , 𝜈 LT = 𝜈 TT = 0.25
where L refers to the fiber direction, T refers to the transverse direction. mesh: only a quarter of the structure is meshed. A space ratio is considered in these two directions (ratio between the size of the larger and the smaller element). Different meshes are considered.

See Fig. 3 for instance (N x = 44, N y = 10 with a space ratio of 8).

numerical layers: N z is the total number of numerical layers.

number of dofs for VS-LD4-S:

Ndof xy = 3(N x + 1)(N y + 1)
and Ndof z = 12 × 𝛼NC + 3 are the number of dofs of the two problems associated with v i j and f i j respectively. 𝛼 is the number of numerical layers per physical layer. So the total number of dofs is

Ndof xy + Ndof z .
number of dofs for LD4-S:

N LW = 3.(4.N z + 1).(3.N x .N y + 2(N x + N y ) + 1)

results:

The results are made nondimensional using:

u = u 1 (a  ∕2, b  ∕2, z) E T eq 0 S 3 , w = u 3 (0, b  ∕2, z) 100E T eq 0 S 4 .

𝜎 𝛼𝛼 =

𝜎 𝛼𝛼 (0,b  ∕2,z)

q 0 S 2 ,. 𝜎 13 = 𝜎 13 (a  ∕2,b  ∕2,z) q 0 S , 𝜎 33 = 𝜎 33 (0,b  ∕2,z) q 0
An error rate is defined as

Δ AV = 100 |x -x AV | |x AV |
where x = w, 𝜎 13 [START_REF] Vidal | Multi-resolution strategies for the modeling of composite shell structures based on the variable separation method[END_REF] the subscript AV stands for the asymptotic value (corresponding to the refined mesh). First of all, a convergence study is performed for the three approaches (Figs. 4 and5). A three-layer configuration with two slenderness ratios (thick and thin shells) is considered. The variations of the transverse displacement, the transverse shear stress and the error rate Δ AV (Eq. ( 43)) with respect to the number of elements along the 𝜉 axis are reported in these figures. A mesh with a space ratio equal to 8 in the Fig. 7. distribution of 𝜎 11 along the thickness -convergence study on the couples -S = 10 -

𝜑 = 𝜋∕2-3 layers [0 • ∕90 • ∕0 • ] -VS-LD4-S model.
= 𝜋∕2-3 layers [0 • ∕90 • ∕0 • ] -VS-LD4-S model.
two directions is used. It can be inferred from these results that the convergence rate of the shell theory (LD4-M model) is better than the other one (the degenerated shell). As it will be discussed in the subsequent examples, the asymptotic value is different for the thick case between the LD4-M and LD4-S models due to the approximation of the geometry of the degenerated shell approach. We can also notice that the variable separation drives to the same results for the displacement. As far as the transverse shear stress is concerned, the convergence rate is different, but the asymptotic value is very close.

Based on this study, a 44 × 10 mesh will be used in the following.

Then, the study is focused on the VS-LD4-S model in order to illustrate the behavior of the PGD method. A three-layer shell with S = 10 is considered. The behavior of the method is given in Figs. 678. These figures the convergence of the PGD process on both displacements and stresses. The two first couples allow us to build a solution very close to the converged one. The following couples bring only small corrections on the transverse stresses.

further illustration, the functions f i 1 (z) are given in Fig. 9. The two first functions are constant and linear, corresponding to classical models (FSDT). Then, the distribution of the two other couples is more complex, they represent small contributions to the whole solution.

The present method is also assessed for different slenderness ratios, from thick to thin shell. In Table 2, the results are compared with the LD4-S approach for displacements, in-plane and transverse stresses.

Table 2

Comparison VS-LD4-S/LD4-S -[0

• ∕90 • ∕0 • ] -𝜑 = 𝜋∕2 -N z = NC. S u(e∕2) w(0) 𝜎 11 (-

Table 3

Comparison Two, three or four couples are built depending on the slenderness ratio. The excellent agreement with the LD4-S model shows that the variable separation does not introduce additional errors for the solution regardless of the value of S.

VS-LD4-S/LD4-M -[0 • ∕90 • ∕0 • ] -𝜑 = 𝜋∕2 -N z = NC. S u(h∕2) w(0) 𝜎 11 (-
The results are also compared with the LD4-M model (cf. Table 3). It can be inferred from this table that the slenderness ratio has an influence on the results. For very thick case, the simplified assumption on the geometry of the shell drives to a significant error rate (from 10 to 21%). For further illustration, Fig. 10 shows the distribution of the inplane and transverse normal stresses through the thickness at the center of the structure. The main difference between the two approaches occurs at the bottom of the shell. Nevertheless, we can notice that the present method can be applied for semi-thick to thin structures as the error rate becomes acceptable. Finally, the approach is assessed on a semi-thick shell (S = 10) for different numbers of layers (NC = 6, 9, 24). The choice of this configuration comes from the previous analysis where the error rate remains good. From Table 4, we can notice that the error rate is less than 4.5% for the displacements and 8.7% for the stresses when compared with the reference solution. Through-thickness distribution of stresses for a 9-layer configuration are also plotted in Fig. 11. The present distributions are in excellent agreement with the LD4-S model and the accuracy 12. shell panel with rigid diaphragm. of this approach is very good with respect to the reference solution. Equivalent Single models fail to model shell structures with these stacking sequences.

Table 5

Comparison VS-LD4-S/LD4-S -rigid diaphragm -NC = 3.

S

= u 3 (a  ∕2, b  ∕2, z) 100E T eq 0 S 4 . 𝜎 𝛼𝛼 = 𝜎 𝛼𝛼 (a  ∕2,b  ∕2,z) q 0 S 2 ,. 𝜎 33 = 𝜎 33 (a  ∕2,b  ∕2,z) q 0 reference solution: LD4-M model
For this test case, different values of the parameter R∕a  are considered referring to deep or shallow shells. The convergence in terms of number of couples is given in Fig. 13 for S = 20, R∕a  = 2. A high convergence rate is found for the transverse displacement, while more couples are needed to obtain the transverse normal stress with accuracy. For this configuration, eight couples are built. Further illustrations are shown in Fig. 14 where different values of R∕a  are considered. The previous remark remains true. The convergence rate of the different mechanical quantities is not the same. A small error rate is obtained for the transverse displacement with few couples, that is not the case for 𝜎 33 . Moreover, it can be noticed that the configuration with R∕a  = 2 seems to be the most severe case for this example.

Numerical results are also summarized in Table 5 where the two present approaches are compared for two slenderness ratios. The results are in very good agreement for the transverse displacements and the inplane stress. But, a substantial error can be noted for the transverse stress. As in Ref. [START_REF] Vidal | Multi-resolution strategies for the modeling of composite shell structures based on the variable separation method[END_REF], it can be reduced by updating the z-functions f i . This is achieved efficiently (see Table 6), and the maximum error rate for 𝜎 33 decreases from 19% to 1.8%. The results can be also compared with the reference solution (LD4-M). The displacement drives to a maximum error rate of 4.9%, while the maximal error rate on the stresses is 9%. It is related to the case R∕a  = 2. This error comes from the assumptions on the geometry of the shell. To show the wide range of configurations involved in this test, Fig. 15 R∕a  = 1, 2, 10, 20. Despite the complexity of these distributions, the present approach drives to satisfactory results when compared with the reference solution.

Another aspect is addressed in Fig. 16. For S = 10 and R∕a  = 20, oscillations in the first and third layers occur in the distribution of the transverse normal stress (Fig. 16 left). It can be overcome by introducing numerical layers in each physical layer (N z = 6). The new results are given in Fig. 16 right. As the problem to be solved is split into an inplane/out-of-plane problem, it does not affect the computational cost of the method as in Refs. [START_REF] Vidal | Shell finite element based on the proper generalized decomposition for the modeling of cylindrical composite structures[END_REF][START_REF] Vidal | Multi-resolution strategies for the modeling of composite shell structures based on the variable separation method[END_REF]. Note that the same phenomenon occurs for the in-plane stress, but it is not shown for brevity reason.

Spherical shell

A doubly-curved laminated composite structure is considered in this section. The test is described as follows geometry: spherical square composite cross-ply shell, R = 3, R∕a  = 5, with 3/4/7/10 layers. S = 25. All the layers have the same thickness. boundary conditions: clamped shell on all edges, subjected to a constant pressure q 0 material properties: same materials as in Section 6.1 mesh: only a quarter of the structure is meshed. The distribution of 𝜎 11 along the thickness at the center of the shell is given in Fig. 17 for different numbers of layers. Seven couples are computed for VS-LD4-S model. It can be inferred from this figure that the results of the present approaches are in very good agreement with the reference solution. For such structures, the use of variable separation does not affect the accuracy of the results (Fig. 17 left). Moreover, the number of unknowns involved in the VS-LD4-S model remains low despite the increase of the number of layers.

Conclusion

In this article, two degenerated shell approaches are applied in conjunction with a LayerWise approach and a variable separation method for the modeling of composite structures. A fourth-order LW description of the displacements is used. The two models are assessed through a wide range of configurations involving different geometries, stacking sequences and boundary conditions. The results are compared with a reference solution which comes from a model available in open literature. Accurate results have been obtained despite the assumptions on the geometry of the shell. The approaches have the capability to model semi-thick to thin structures. Moreover, if needed, numerical layers in each physical layer can be introduced without increasing significantly the computational cost for the model based on the variable separation.

Based on these promising results, the membrane and shear locking phenomenon will be addressed in future investigation. , i = 1, 2, 3

𝚺 p z ( f ) = [ f 1 Id 3×3 ⋮ f 2 Id 3×3 ⋮ f 3 Id 3×3 ] (A.5) 𝚺 n z ( f ) = [ f ′ 1 ⋮ f ′ 2 ⋮ f ′ 3 ] (A.6) 𝚺 s1 z ( f ) = [ f 1 Id 2×2 ⋮ f 2 Id 2×2 ⋮ f 3 Id 2×2 ] (A.7) 𝚺 s3 z ( f ) = [ f ′ 1 Id 2×2 ⋮ f ′ 2 Id 2×2 ⋮ f ′ 3 Id 2×2 ] (A.8)
where Id i×i is the i × i identity matrix.

Fig. 1 .

 1 Fig. 1. Description of the geometry of the 8-node shell element.

Fig. 3 .

 3 Fig. 3. Simply supported cylindrical shell (left) -mesh of a quarter of the structure (right).

Fig. 5 .

 5 Fig. 5. convergence study -S = 50 -𝜑 = 𝜋∕2-3 layers [0 • ∕90 • ∕0 • ].

Fig. 6 .

 6 Fig. 6. distribution of u w (right) the thickness -convergence study on the couples -S = 10 -𝜑 = 𝜋∕2-3 layers [0 • ∕90 • ∕0 • ] -VS-LD4-S model.

Fig. 8 .

 8 Fig. 8. distribution of 𝜎 13 (left) and 𝜎 33 (right) along the thickness -convergence study on the couples -S = 10 -𝜑 = 𝜋∕2-3 layers [0 • ∕90 • ∕0 • ] -VS-LD4-S model.

Fig. 9 .

 9 Fig. 9. distribution of f i 1 (i = 1, 2, 3, 4) along the thickness -S = 10 -𝜑 = 𝜋∕2-3 layers [0 • ∕90 • ∕0 • ] -VS-LD4-S model.

Fig. 10 .

 10 Fig. 10. distribution of 𝜎 11 (left) and 𝜎 33 (right) along the thickness -S = 4 -3 layers.

Fig. 11 .

 11 Fig. 11. distribution of 𝜎 11 (left), 𝜎 13 (middle) and 𝜎 33 (right) along the thickness -S = 10 -9 layers.

Fig. 13 . 2 6. 2 .Fig. 12 . 3 Fig. 14 .

 132212314 Fig. 13. distribution of the error with respect to the number of couples -NC = 3 -rigid diaphragm -S = 20 -R∕a  = 2

Fig. 15 .

 15 Fig. 15. distribution of 𝜎 11 through the thickness -NC = 3 -rigid diaphragm -S = 20.

Fig. 16 . 20 Fig. 17 .

 162017 Fig. 16. distribution of 𝜎 33 along the thickness -N z = 3 (left)/N z = 6 (right) numerical layers -NC = 3 -rigid diaphragm -S = 10 -R∕a  = 20

S 2 ,

 2 N x = 24, N y = 24 number of dofs for VS-LD4-S: Ndof xy = 1875 and Ndof z = 12 × NC + 3. Ndof z = 39 to Ndof z = 123 (3/10 layers) number of dofs for LD4-S: N LW = 71175 for N z = 3, N LW = 224475 for N z = 10 results: The results are made nondimensional using: 𝜎 11 = 𝜎 11 (a  ∕2,b  ∕2,z) q 0 reference solution: LD4-M model

Appendix A. 2 .

 2 Expression of 𝚺 □ new ( f ) and MM p LR 𝚺 p new ( f ) = 𝚺 p z ( f ) element problem to be solved on Ω z Appendix B.1. Expression of BaseNz and BaseBz

Table 1

 1 Reference frames used for the description of the shell.

	(e 1 , e 2 , e 3 )	Direct orthonormal Cartesian basis, (G)
	(X 1 , X 2 , X 3 )	Global Cartesian coordinates, (G)
	(t 1 , t 2 , t 3 )	Direct orthonormal local basis, plane tangent to Ω, local (L)
	(x 1 , x 2 , x 3 = z)	Local curvilinear coordinates on Ω, local (L)
	(a 1 , a 2 , a 3 )	Reduced local basis, plane tangent to Ω, Reduced (R)
	(𝜉, 𝜂)	Reduced curvilinear coordinates on Ω, Reduced (R)
	a 𝛼 = 𝚽(𝜉, 𝜂, z) ,𝛼 for 𝛼 = 𝜉, 𝜂

e 1 t 1 • e 2 t 1 • e 3 t 2 • e 1 t 2 • e 2 t 2 • e 3 t 3 • e 1 t 3 • e 2 t 3 • e 3

	t 11 t 12 t 13 t 21 t 22 t 23 t 31 t 32 t 33	⎤ ⎥ ⎥ ⎥ ⎦	=	⎡ ⎢ ⎢ ⎢ ⎣	t 1 • ⎤ ⎥ ⎥ ⎥ ⎦

t 11 t 21 t 31 t 12 t 22 t 32 t 13 t 23 t 33

Table 4

 4 Comparison VS-LD4-S/LD4-M -𝜑 = 𝜋∕2 -S = 10 -N z = NC.

	NC		u(e∕2)	w(0)	𝜎 11 (-e∕2)	𝜎 13 max	𝜎 33max
	6	VS-LD4-S	70.8408 (2.7%)	1.6747 (2.6%)	-4.4873 (7.3%)	1.5219 (4.6%)	-5.9020 (5.6%)
		LD4-M	72.7697	1.7194	-4.8381	1.5954	-6.2529
	9	VS-LD4-S	54.1182 (4.4%)	1.2318 (4.5%)	-3.7848 (8.7%)	1.3624 (3.3%)	-5.2850 (6.4%)
		LD4-M	56.6165	1.2895	-4.1451	1.4096	-5.6479
	24	VS-LD4-S	67.3724 (4.0%)	1.5552 (4.0%)	-4.7467 (8.6%)	1.4873 (3.1%)	-5.8211 (4.0%)
		LD4-M	70.1935	1.6203	-5.1951	1.5355	-6.0636

Table 6

 6

	Comparison VS-LD4-S/LD4-S/LD4-M -rigid diaphragm -NC = 3 -f i (z) updated.
	S	R∕a 			w(0)			𝜎 11 (-h∕2)			𝜎 33 max
			1	2	10	20	1	2	10	20	1	2	10	20
	10	VS-LD4-S	2.6971	9.6868	13.4184	13.4266	0.5306	-0.1826	0.1293	-0.2365	1.2152	1.1658	1.5307	1.6107
		Err-LD4-S	0.04%	0.02%	0.01%	0.00%	0.44%	2.85%	0.13%	0.40%	1.82%	0.51%	0.79%	0.92%
		Err-LD4-M	4.80%	4.84%	4.78%	4.78%	3.30%	6.71%	3.74%	2.42%	1.50%	2.66%	0.22%	3.59%
		LD4-S	2.6959	9.6852	13.4173	13.4259	0.5330	-0.1775	0.1292	-0.2374	1.1934	1.1718	1.5187	1.6256
		LD4-M	2.8331	10.1792	14.0916	14.1006	0.5487	0.1711	0.1344	-0.2424	1.1972	1.1356	1.5274	1.6706
	20	VS-LD4-S	0.7778	5.1588	13.2326	13.2695	0.4247	0.1159	0.1336	0.1257	1.9206	1.2228	2.3452	2.4346
		Err-LD4-S	0.03%	0.02%	0.00%	0.00%	0.35%	1.86%	0.02%	0.02%	1.44%	1.48%	1.88%	0.01%
		Err-LD4-M	2.44%	2.49%	2.47%	2.46%	2.46%	9.41%	1.90%	1.97%	0.33%	6.66%	3.98%	1.18%
		LD4-S	0.7776	5.1581	13.2324	13.2693	0.4262	0.1128	0.1336	0.1257	1.8934	1.1757	2.3019	2.4348
		LD4-M	0.7972	5.2905	13.5671	13.6046	0.4354	0.1050	0.1362	0.1282	1.9142	1.1185	2.2555	2.4061

  TA pi = t 1 . e i . a 1 N 𝜉,𝜉 t 1 . a 2 N 𝜉,𝜂 t 1 . a 1 N 𝜉,𝜂 + t 1 . a 2 N 𝜉,𝜉 . a 1 N 𝜉,𝜉 t 2 . a 2 N 𝜉,𝜂 t 2 . a 1 N 𝜉,𝜂 + t 2 . a 2 N 𝜉,𝜉

	with
						[	]
	TA s1i = t 3 . e i		N 𝜉,𝜂 N 𝜉,𝜉	, i = 1, 2, 3
	BaseTA s3 =	⎡ ⎢ ⎢ ⎢ ⎣	TA s31 0 2×1 0 2×1	0 2×1 TA s32 0 2×1	0 2×1 0 2×1 TA s33	⎤ ⎥ ⎦ ⎥ ⎥	(A.4)
	with
	TA s3i = t 1 . e i	[	t 1 . a 2 N 𝜉 t 1 . a 1 N 𝜉	]	+ t 2 . e i	[	t 2 . a 1 N 𝜉 t 2 . a 2 N 𝜉	]
					⎡ ⎢ ⎢ ⎢ ⎣		t 1 ⎤ ⎥ ⎥ ⎥ ⎦	+ t 2 . e i	⎡ ⎢ ⎢ ⎢ ⎣	t 2 ⎤ ⎥ ⎦ ⎥ ⎥	, i = 1, 2, 3
	BaseTA n =	⎡ ⎢ ⎢ ⎢ ⎣	t 3 . e 1 N 𝜉 0 0	0 t 3 . e 2 N 𝜉 0	0 0 t 3 . e 3 N 𝜉	⎤ ⎥ ⎦ ⎥ ⎥	(A.2)
	BaseTA s1 =	⎡ ⎢ ⎢ ⎢ ⎣	TA s11 0 2×1 0 2×1	0 2×1 TA s12 0 2×1	0 2×1 0 2×1 TA s13	⎤ ⎥ ⎦ ⎥ ⎥	(A.3)

  tta 12 ṽ2,𝜉 tta 13 ṽ3,𝜉 tta 21 ṽ1,𝜂 tta 22 ṽ2,𝜂 tta 23 ṽ3,𝜂 tta 11 ṽ1,𝜂 + tta 21 ṽ1,𝜉 tta 12 ṽ2,𝜂 + tta 22 ṽ2,𝜉 tta 13 ṽ3,𝜂 + tta 23 ṽ3,𝜉 with tta ij = t 1 .e j t 1 .a i + t 2 .e j t 2 .a i . e 1 ṽ1 t 3 . e 2 ṽ2 t 3 . e 3 ṽ3 . e 1 ṽ1,𝜂 t 3 . e 2 ṽ1,𝜂 t 3 . e 3 ṽ1,𝜂 t 3 . e 1 ṽ1,𝜉 t 3 . e 2 ṽ1,𝜉 t 3 . e 3 ṽ1,𝜉 [ tta 21 ṽ1 tta 22 ṽ2 tta 23 ṽ3 tta 11 ṽ1 tta 12 ṽ2 tta 13 ṽ3

	𝚺 n xy (ṽ) =	[ t 3 ]	(B.4)
	𝚺 s1 xy (ṽ) =	[	t 3 ]	(B.5)
	𝚺 s3 xy (ṽ) =
			⎤ ⎥ ⎦ ⎥ ⎥	(B.3)

] (B.6)

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Finite element problem to be solved on Ω

Appendix A.1. Expression of BaseTA □ and 𝚺 □ z ( f )

with


Keywords: Functionally graded materials, Layerwise approach, Separation of variables, Finite Element

This work deals with the modeling of plate structures constituted of Functionally Graded Material through a variable separation approach. The displacement field is approximated as a sum of separated functions of the in-plane coordinates x, y and the transverse coordinate z. This choice yields to a non-linear problem that can be solved by an iterative process.

Two 2D and 1D linear problems are solved successively. In the thickness direction, a fourthorder expansion in each layer is considered. For the in-plane description, classical Finite Element method is used. The number of unknowns is reduced compared to classical layerwise approach.

Classical numerical tests encountered in the literature are provided to show the accuracy of the present LayerWise (LW) method. The plate can be constituted of a FGM layer or combined FGM/orthotropic layers. Different slenderness ratios and anisotropy indexes are also considered. Accurate results with a low computational cost are obtained when compared with reference solutions and models available in literature.

Introduction

In this study, Functionally Graded Materials (FGM) in which the mechanical properties vary smoothly and continuously from one surface to the other are considered. To obtain that, different ways of manufacturing process can be performed (see [START_REF] Naebe | Functionally graded materials: A review of fabrication and properties[END_REF][START_REF] Parihar | Recent advances in the manufacturing processes of functionally graded materials: a review[END_REF]). One of the advantages is the elimination of the transverse stresses discontinuity which can be encountered in laminated composite. Thus, it allow us to avoid one of the most critical failure mode, i.e. the delamination process. Other attractive characteristics can be mentionned such as: enhanced thermal properties, higher fracture toughness, ... [START_REF] Birman | Modeling and analysis of functionally graded materials and structures[END_REF]. Note that FGM can be used for a wide field of applications (see [START_REF] Gupta | Recent development in modeling and analysis of functionally graded materials and structures[END_REF][START_REF] Parihar | Recent advances in the manufacturing processes of functionally graded materials: a review[END_REF]). That justifies the increasing interest for this type of materials. Thus, it will be important to have efficient numerical tools so as to compute accurate displacements and stresses for design purpose.

These topics have been the focus of attention of many researchers during the past few years.

Different three-dimensional models have been carried out in [START_REF] Kashtalyan | Three-dimensional elasticity solution for bending of functionally graded rectangular plates[END_REF][START_REF] Kashtalyan | Three-dimensional elasticity solution for sandwich panels with a functionally graded core[END_REF][START_REF] Burlayenko | Free vibrations and static analysis of functionally graded sandwich plates with three-dimensional finite elements[END_REF] for FGM and sandwich plates, but the computational cost becomes high. Thus, a wide range of 2D models with a low number of generalized unknowns has been developed. For that, the displacements are expressed through the whole thickness of the plate (Equivalent Single Layer model -ESL).

Classical Laminate Theory and First-order Shear Deformation Theory (FSDT) have been considered in [START_REF] Chi | Mechanical behavior of functionally graded material plates under transverse load -part i: Analysis[END_REF] and [START_REF] Cheng | Deflection relationships between the homogeneous kirchhoff plate theory and different functionally graded plate theories[END_REF], respectively. But, the first one fails to predict accurate results for thick cases and the second one needs a shear correction factor. So, Higher-order Shear Deformation Theory (HSDT) has been proposed considering different types of shear strain shape functions: polynomial [START_REF] Reddy | Analysis of functionally graded plates[END_REF][START_REF] Cheng | Deflection relationships between the homogeneous kirchhoff plate theory and different functionally graded plate theories[END_REF][START_REF] Ferreira | Static analysis of functionally graded plates using third-order shear deformation theory and a meshless method[END_REF][START_REF] Kumar | New transverse shear deformation theory for bending analysis of FGM plate under patch load[END_REF], trigonometric [START_REF] Kumar | New transverse shear deformation theory for bending analysis of FGM plate under patch load[END_REF], inverse trigonometric [START_REF] Nguyen | A new inverse trigonometric shear deformation theory for isotropic and functionally graded sandwich plates[END_REF], specified [START_REF] Nguyen-Xuan | Isogeometric analysis of functionally graded plates using a refined plate theory[END_REF][START_REF] Li | A novel shear deformation theory for static analysis of functionally graded plates[END_REF] forms. But, as it has been highlighted in [START_REF] Carrera | Effects of thickness stretching in functionally graded plates and shells[END_REF], the stretching effect is of major importance for the modeling of this type de materials. Thus, these HSDT models have been extended to take into account this phenomenon. They can be again classified depending of the strain shape functions as previously: polynomial [START_REF] Qian | Static and dynamic deformations of thick functionally graded elastic plates by using higher-order shear and normal deformable plate theory and meshless local petrov-galerkin method[END_REF][START_REF] Natarajan | Bending and vibration of functionally graded material sandwich plates using an accurate theory[END_REF][START_REF] Thai | Improved refined plate theory accounting for effect of thickness stretching in functionally graded plates[END_REF], trigonometric [START_REF] Zenkour | Generalized shear deformation theory for bending analysis of functionally graded plates[END_REF][START_REF] Neves | A quasi-3d sinusoidal shear deformation theory for the static and free vibration analysis of functionally graded plates[END_REF][START_REF] Mantari | Generalized hybrid quasi-3d shear deformation theory for the static analysis of advanced composite plates[END_REF][START_REF] Mantari | Optimized sinusoidal higher order shear deformation theory for the analysis of functionally graded plates and shells[END_REF][START_REF] Thai | A simple quasi-3d sinusoidal shear deformation theory for functionally graded plates[END_REF], exponential [START_REF] Mantari | Five-unknowns generalized hybrid-type quasi-3d hsdt for advanced composite plates[END_REF], hyperbolic [START_REF] Belabed | An efficient and simple higher order shear and normal deformation theory for functionally graded material (FGM) plates[END_REF][START_REF] Mantari | Four-unknown quasi-3d shear deformation theory for advanced composite plates[END_REF], hybrid [START_REF] Mantari | Generalized hybrid quasi-3d shear deformation theory for the static analysis of advanced composite plates[END_REF][START_REF] Mantari | Five-unknowns generalized hybrid-type quasi-3d hsdt for advanced composite plates[END_REF] forms. Note that most of these approaches are built such that the free boundary conditions on the top and bottom surfaces are fullfiled. It should be also mentioned that some of these works needs the determination of one or two parameters involved in the expression of the kinematics, see [START_REF] Mantari | Optimized sinusoidal higher order shear deformation theory for the analysis of functionally graded plates and shells[END_REF][START_REF] Mantari | Four-unknown quasi-3d shear deformation theory for advanced composite plates[END_REF][START_REF] Nguyen | A new inverse trigonometric shear deformation theory for isotropic and functionally graded sandwich plates[END_REF][START_REF] Kumar | New transverse shear deformation theory for bending analysis of FGM plate under patch load[END_REF][START_REF] Li | A novel shear deformation theory for static analysis of functionally graded plates[END_REF]. This type of models can be well-suited for monolayered structure, but limitations appear considering multi-layered or sandwich panels. An alternative and simple way consists in including the so-called Murakami's zig-zag function as in [START_REF] Natarajan | Bending and vibration of functionally graded material sandwich plates using an accurate theory[END_REF], but more precised theories have been developed based on a LayerWise approach (LW). Note the extensive works based on the systematic approach developed by Carrera with his Carrera's Unified Formulation (CUF) and its extension to FGM material in [START_REF] Carrera | Variable kinematic model for the analysis of functionally graded material plates[END_REF]. Displacement-based and Reissner's Mixed Variational Theorem (RMVT) approaches are addressed in [START_REF] Brischetto | Advanced mixed theories for bending analysis of functionally graded plates[END_REF][START_REF] Carrera | Effects of thickness stretching in functionally graded plates and shells[END_REF]. Discrete layers have been also considered in [START_REF] Ramirez | Static analysis of functionally graded elastic anisotropic plates using a discrete layer approach[END_REF]. Only partial works on the modeling of FGM plate structures have been given herein and interested readers can refer to the reviews provided in [START_REF] Jha | A critical review of recent research on functionally graded plates[END_REF][START_REF] Swaminathan | Stress, vibration and buckling analyses of FGM plates -a state-of-the-art review[END_REF][START_REF] Gupta | Recent development in modeling and analysis of functionally graded materials and structures[END_REF][START_REF] Thai | A review of theories for the modeling and analysis of functionally graded plates and shells[END_REF][START_REF] Ghatage | On the numerical modelling and analysis of multidirectional functionally graded composite structures: A review[END_REF].

Over the past years, some methods based on the separation of variables have shown interesting features to model composite structures. For FGM material, the extended Kantorovich method has been applied in [START_REF] Aghdam | A semi-analytical solution for bending of moderately thick doubly curved functionally graded panels[END_REF]. Another interesting way is based on the so-called Proper Generalized Decomposition (PGD). It has been successfully developed for the modeling of laminated and sandwich structures [START_REF] Savoia | A variational approach to three-dimensional elasticity solutions of laminated composite plates[END_REF][START_REF] Bognet | Advanced simulation of models defined in plate geometries: 3D solutions with 2D computational complexity[END_REF][START_REF] Vidal | Composite beam finite element based on the proper generalized decomposition[END_REF][START_REF] Vidal | Proper generalized decomposition and layer-wise approach for the modeling of composite plate structures[END_REF] and will be extended hereafter for the analysis of FGM plates. In this purpose, the displacements are written under the form of a sum of products of bidimensional functions of (x,y) and unidimensional functions of z. A piecewise fourth-order Lagrange polynomial of z is chosen and a 2D eight-node quadrilateral FE is employed for the in-plane coordinates. Each unknown function of (x,y) is classically approximated using one degree of freedom (dof) per node of the mesh and the LW unknown functions of z are global for the whole plate. Finally, the deduced non-linear problem implies the resolution of 2D and 1D problems alternatively, in which the number of unknowns is smaller than a classical Layerwise approach. Considering the models developed in open literature, it calls for a well-suited kinematics, i.e. a relevant shear strain shape function.

The present work allows us to determine automatically the z-functions. Moreover, any configurations could be considered, namely structures constituted of only FGM or combined FGM/classical (laminated) layers.

The article is organized as follows. Firstly, the classical mechanical problem involving a FGM is given. Then, the formulation of the PGD is described in this particular framework.

The assumption on the displacements yields a non-linear problem that is solved using an iterative process. The FE discretization is also described. The separation variables can be efficiently applied as the in-plane and out-of-plane integrations can be split. Finally, numerical results are provided to assess the accuracy of the method. Classical examples are addressed and the results are compared with both reference solutions and results available in open literature. The influence of the slenderness ratio and a characteristic parameter of the FGM (volume fraction exponent) is studied. The present approach also allows us to modelize the FGM with piecewise constant mechanical properties through the thickness at a low computational cost.

Reference problem description: the governing equations

Let us consider a composite plate structure occupying the domain

V = Ω × Ω z where Ω = [0, a] × [0, b] (a, b being the dimensions of the plate) and Ω z = [-h 2 , h 2 ]
in a Cartesian coordinate (x, y, z). h is the thickness of the plate, see Fig. 1.

[Figure 1 about here.]

Constitutive relation

As we are interested in the analysis of FGM plates, the continuous variations of the material characteristics along the thickness have to be well described. The FGM layers can be single or embedded in other isotropic/orthotropic/FGM layers. Classically, exponential and/or polynomial functions are applied to either engineering constants (Young Modulus, shear Modulus, Poisson ratio) or directly to material stiffnesses C ij .

The plate can be made of NC perfectly bonded classical/FGM layers. The constitutive equations for a layer k can be written as

σ (k) = C (k) (z)ε z Ω k z with Ω z = ∪ NC k=1 Ω k z ( 1
)
where we denote the stress vector by σ, the strain vector via ε.

We have

C (k) (z) =               C (k) 11 C (k) 12 C (k) 13 0 0 C (k) 16 C (k) 22 C (k) 23 0 0 C (k) 26 C (k) 33 0 0 C (k) 36 C (k) 44 C (k) 45 0 sym C (k) 55 0 C (k) 66               (2)
where C

(k) ij are the three-dimensional stiffness coefficients of the layer (k).

In the present work, without loss of generality, we assume that we can write the constitutive law as

C (k) (z) = g (k) (z) C (k) 0 (3)
where C

(k) 0 is constant in each layer. If needed, the decomposition of the matrix C (k) (z) can be splited into several terms to keep the efficiency of the method.

Note that g (k) (z) is constant for a classical isotropic or orthotropic layer.

The weak form of the boundary value problem

In the present work, a static problem is studied. For conciseness reason, the body force density is not considered. The surface force t is applied on ∂V F = ∂V Fxy ×{z F }. We assume that a prescribed displacement u = 0 is imposed on Γ D . Thus, the classical problem to be solved can be formulated as follows:

Using the above matrix notations and for admissible displacement δu ∈ δU, the variational principle is given by:

find u ∈ U such that - V ε(δu) T σ dV + ∂V F δu T t d∂V = 0, ∀δu ∈ δU ( 4
)
where U is the space of admissible displacements, i.e. U = {u ∈ (H 1 (V)) 3 /u = 0 on Γ D }.

We have also δU

= {u ∈ (H 1 (V)) 3 /u =0 on Γ D }.
This problem can be expressed as

a(u, δu) = b(δu) ∀δu ∈ δU (5) with a(u, δu) = V ε(δu) T C(z)ε(u) dΩ dΩ z b(δu) = ∂V F δu T t d∂V (6)

Application of the separated representation to the FGM plate

The present section focuses on the procedure that allows us to expand the PGD to the analysis of FGM layers. This specific separation has shown interesting features in the framework of static analysis [START_REF] Bognet | Advanced simulation of models defined in plate geometries: 3D solutions with 2D computational complexity[END_REF][START_REF] Vidal | Proper generalized decomposition and layer-wise approach for the modeling of composite plate structures[END_REF]. The separated formulation is briefly introduced in the subsequent sections.

The displacement and the strain field

In the framework of the PGD, the displacement solution u is classically written as the sum of N products of separated functions (N ∈ N + is the order of the representation)

u(x, y, z) = N i=1 f i (z) • v i (x, y) ( 7 )
where f i (z) and v i (x, y) are unknown functions which must be computed during the resolution process. f i (z) and v i (x, y) are defined on Ω z and Ω respectively. The "•" operator in Eq. ( 7) is Hadamard's element-wise product. We have:

f i • v i = v i • f i =      f i 1 (z)v i 1 (x, y) f i 2 (z)v i 2 (x, y) f i 3 (z)v i 3 (x, y)      with v i =      v i 1 (x, y) v i 2 (x, y) v i 3 (x, y)      f i =      f i 1 (z) f i 2 (z) f i 3 (z)      (8)
The strain can be expressed with respect to the reference frame in which the dependance with respect to the space coordinates is omitted as follows:

ε(u) = N i=1               f i 1 v i 1,1 f i 2 v i 2,2 (f i 3 ) v i 3 (f i 2 ) v i 2 + f i 3 v i 3,2 (f i 1 ) v i 1 + f i 3 v i 3,1 f i 1 v i 1,2 + f i 2 v i 2,1               (9)
where the prime stands for the classical derivative (f i = df i dx ), and () ,α for the partial derivative.

Formulation of the problem to be solved

An iterative procedure is performed to solve Eq. ( 5) (greedy algorithm). If we assume that the first n functions have been already computed, the trial function for the iteration n + 1 is written as

u n+1 = u n + f • v = u n +      f 1 v 1 f 2 v 2 f 3 v 3      (10)
where (f 1 , f 2 , f 3 ) and (v 1 , v 2 , v 3 ) are the functions to be computed and u n is the associated known set at iteration n defined by

u n = n i=1      f i 1 v i 1 f i 2 v i 2 f i 3 v i 3      (11)
The test functions are

δ(f • v) = δf • v + δv • f (12) with v =      v 1 v 2 v 3      f =      f 1 f 2 f 3      (13)
Introducing the test function defined by Eq. ( 12) and the trial function defined by Eq.

(10) into the weak form Eq. ( 5), the two following equations to be solved can be deduced:

• for the test function δf a(v • f, v • δf) = b(v • δf) -a(u n , v • δf) ∀δf ( 14
)
• for the test function δv

a(f • v, f • δv) = b(f • δv) -a(u n , f • δv) ∀δv ( 15
)
From Eq. ( 14) and Eq. ( 15), a coupled non-linear problem is derived. A fixed point method is a simple way to solve it. Starting from an initial function ( f(0) , ṽ(0) ), we construct a sequence ( f(l) , ṽ(l) ) which satisfy Eq. ( 14) and Eq. ( 15) respectively. Only 1D or 2D functions have to be found, as the other one is assumed to be known (from the previous step of the fixed point strategy). The fixed point algorithm is stopped when the distance between two consecutive terms are sufficiently small.

Finite element discretization

In this section, a classical finite element approximation in Ω and Ω z for (v, f) is introduced. The elementary vector of degrees of freedom (dof) associated with one element Ω e of the mesh in Ω is denoted q v e . The elementary vector of dofs associated with one element Ω ze of the mesh in Ω z is denoted q f e . The displacement and strain fields are determined from the values of q v e and q f e by

v e = N xy q v e , E e v = B xy q v e , f e = N z q f e , E e f = B z q f e ( 16
)
where

E e v T = v 1 v 1,1 v 1,2 v 2 v 2,1 v 2,2 v 3 v 3,1 v 3,2 E e f T = f 1 f 1 f 2 f 2 f 3 f 3
The matrices N xy , B xy , N z , B z contain the interpolation functions, their derivatives and the jacobian components.

Finite element problem to be solved on Ω

For the sake of simplicity, the functions f(l) which are assumed to be known, will be denoted f. And the function ṽ(l) to be computed will be denoted v. The strains included in Eq. ( 15) are defined as

ε( f • v) = Σ z ( f )E v (17)
with

Σ z ( f) =               0 f1 0 0 0 0 0 0 0 0 0 0 0 0 f2 0 0 0 0 0 0 0 0 0 f 3 0 0 0 0 0 f 2 0 0 0 0 f3 f 1 0 0 0 0 0 0 f3 0 0 0 f1 0 f2 0 0 0 0               (18)
For convenience reasons, the displacement u is written under the form

u = D z ( f )v =      f1 0 0 0 f2 0 0 0 f3      v
The variational problem defined on Ω from Eq. ( 15) is

Ω δE v T k z ( f )E v dΩ = ∂V Fxy δv T F z ( f ) d∂V - Ω δE v T σ z ( f , u n )dΩ (19)
where

k z ( f) = NC k=1 Ω k z g (k) (z) Σ z ( f ) T C (k) 0 Σ z ( f )dz σ z ( f , u n ) = NC k=1 Ω k z g (k) (z) Σ z ( f ) T C (k) 0 ε(u n )dz F z ( f ) = D z ( f ) T t z=z F ( 20
)
The introduction of the finite element approximation Eq. ( 16) in the variational Eq. ( 19) leads to the linear problem:

(21) K z (f ˜) q v = R vz ( ˜f, u n )
where • q v is the vector of the nodal displacements, associated with the finite element mesh in Ω,

• K z ( f ) is the mechanical stiffness matrix obtained by summing the elements' stiffness

matrices K e z ( f) = Ωe B T xy k z ( f )B xy dΩ e • R vz ( f , u n ) is the equilibrium residual obtained by summing the elements' residual load vectors R e vz ( f , u n ) = ∂V e Fxy N T xy F z ( f )d∂V e - Ωe B T xy σ z ( f, u n )dΩ e

Finite element problem to be solved on Ω z

As in the previous section, the known functions ṽ(l-1) will be denoted ṽ and the functions f(l) to be computed will be denoted f. The strain included in Eq. ( 14) is defined as

ε(ṽ • f ) = Σ xy (ṽ)E f (22)
where

Σ xy (ṽ) =              
ṽ1,1 0 0 0 0 0 0 0 ṽ2,2 0 0 0 0 0 0 0 0 ṽ3 0 0 0 ṽ2 ṽ3,2 0 0 ṽ1 0 0 ṽ3,1 0 ṽ1,2 0 ṽ2,1 0 0 0

              (23)
and the displacement is

u = D xy (ṽ)E f =              
ṽ1 0 0 0 0 0 0 0 0 0 0 0 0 0 ṽ2 0 0 0 0 0 0 0 0 0 0 0 0 0 ṽ3 0 0 0 0 0 0 0

              E f
The variational problem defined on Ω z from Eq. ( 14) is

NC k=1 Ω k z g (k) (z) δE f T k (k) xy (ṽ)E f dz = δE f T z=z F F xy (ṽ) - NC k=1 Ω k z g (k) (z) δE f T σ xy (k) (ṽ, u n )dz (24)
where k

(k)
xy (ṽ), σ xy (k) (ṽ, u n ) and F xy (ṽ) can be expressed under the following separated form:

k (k) xy (ṽ) = Ω Σ xy (ṽ) T C (k) 0 Σ xy (ṽ)dΩ σ xy (k) (ṽ, u n ) = Ω Σ xy (ṽ) T C (k) 0 ε(u n )dΩ F xy (ṽ) = ∂V Fxy D xy (ṽ) T t d∂V xy ( 25
)
Note that the computation of k

(k)
xy depends on the type of involved layers, i.e. FGM or classical layer.

The introduction of the finite element discretization Eq. ( 16) in the variational Eq. ( 24) leads to the linear problem:

K xy (ṽ) q f = R f (ṽ, u n ) ( 26
)
where

• q f is the vector of degree of freedom associated with the F.E. approximations in Ω z .

• K xy (ṽ) is obtained by summing the elements' stiffness matrices:

K xy e (ṽ) = Ω k ze g (k) (z) B T z k (k) xy (ṽ)B z dz e ( 27
)
• R f (ṽ, u n ) is obtained by summing the residual' vectors:

R e f (ṽ, u n ) = B T z | z=z F F xy (ṽ) - Ω k ze g (k) (z) B T z σ xy (k) (ṽ, u n ) dz e (28)

Numerical results

In the numerical examples, an eight-node quadrilateral FE based on the classical Serendipity interpolation functions is used for the unknowns depending on the in-plane coordinates.

For the unknowns depending on the z-coordinate, the displacement is described by a fourthorder interpolation as it is justified in [START_REF] Vidal | Proper generalized decomposition and layer-wise approach for the modeling of composite plate structures[END_REF]. A Gaussian numerical integration with 3 × 3 points is used to evaluate the elementary matrices. As far as the integration with respect to the transverse coordinate is concerned, 5 integrations points are chosen in each numerical layer. The present approach, denoted VS-LD4, is compared with both reference solutions and other models available in open literature.

Hereafter, two classical numerical examples given by [START_REF] Zenkour | Generalized shear deformation theory for bending analysis of functionally graded plates[END_REF] are addressed to assess the accuracy of the present method. A one-layered plate and a sandwich structure are considered with different slenderness ratios and values of volume fraction exponent. Moreover, a third example from [START_REF] Wu | An rmvt-based third-order shear deformation theory of multilayered functionally graded material plates[END_REF] is studied involving laminates with both FGM and orthotropic layers. Unless otherwise mentioned, the fourth-order layerwise model LD4, referring to the systematic work of Carrera ("Carrera's Unified Formulation" (CUF)) will be considered as the reference solution. This approach is based on a displacement formulation, and each displacement component is expanded until the fourth-order. 12NC + 3 unknown functions are used in this kinematic.

Note that only one couple is built for each test due to the type of load and boundary conditions.

one-layered FGM plate

In this section, one-layered FGM plate is considered with different slenderness ratios.

The data is given as follows:

geometry: square FGM plate with length-to-thickness ratio S = a/h = 4, 10 boundary conditions: simply-supported plate subjected to a bi-sinusoidal pressure: p(x, y) = p 0 sin( πx a ) sin( πy b ) material properties: elastic properties varying along the thickness direction z by a polynomial law, as proposed by Zenkour [START_REF] Kashtalyan | Three-dimensional elasticity solution for sandwich panels with a functionally graded core[END_REF]; The plate is made of aluminum on the bottom and alumina on the top, and the following functional relationship is considered for convergence. We also notice that the in-plane mesh has to be rather fine so that the boundary conditions on the transverse normal stress be fulfilled and it seems to be more difficult for k = 10.

E(z): E(z) = E m + (E c -E m ) 2z + h 2h
[Table 1 about here.]

Considering a mesh with N x = N y = 32 elements, the influence of the numerical layers is addressed in Tab. 2. For the displacements, the convergence rate is high as above. For the stresses, the influence of the value of k on the convergence rate is more pronounced.

For k = 10, a higher Young modulus gradient through the thickness occurs and 8 numerical layers are needed. For k = 1, only 2 are sufficient. The transverse normal stress is the most sensitive value to the mesh.

[Table 2 about here.]

For further illustration, the distributions of σ ¯13 and σ ¯33 through the thickness are given in Fig. 2 for k = 10. It can be inferred from this figure that one numerical layer is insufficient to achieve accurate results. The boundary conditions are satisfied with at least 4 numerical layers.

[Figure 2 about here.]

Piecewise constant Young Modulus: Convergence study

As in [START_REF] Carrera | Variable kinematic model for the analysis of functionally graded material plates[END_REF] and [START_REF] Ramirez | Static analysis of functionally graded elastic anisotropic plates using a discrete layer approach[END_REF], the FGM plate is modelized using a piecewise constant Young modulus. The influence of the number of numerical layers is studied hereafter. E is constant in each numerical layer, the value is taken at the middle of the layer. This type of computations can be performed in a commercial software considering numerous materials. Here, the PGD framework is advantageous as the computational cost is reduced even using high number of layers, as only 1D problems are involved (see [START_REF] Vidal | Proper generalized decomposition and layer-wise approach for the modeling of composite plate structures[END_REF]). To illustrate that, Tab. 3 gives the comparison between the present method (VS-LD4) vs the fourth-order LW expansion approach (LD4) in terms of the number of unkowns. The gain becomes very high when the number of layers increases. Thus, it is interesting to investigate this issue with the present numerical tool. For that, the error with respect to the mesh refinement is given in a log-log curve in Fig. 3 for both displacements and stresses. The convergence rate is the same for all the quantities excepted for the in-plane stress (divided by two). For this latter, numerical layers imply discontinuities between each element in the thickness direction, so high number of numerical layers is required to circumvent this disadvantage. An illustration is shown in Fig. 4. It can be also noted that oscillations occurs even with N z = 64 numerical layers.

[Table 3 about here.]

[Figure 3 about here.]

[Figure 4 about here.]

Different configurations for the one-layer test case

Based on the previous convergence studies, the subsequent computations are performed usinga32×32 mesh with 8 numerical layers. Note that the reference solution (LD4) involves [Table 4 about here.]

[Table 5 The results are summarized in Tab. 7. The reference solution, denoted Q3D, is issued from [START_REF] Carrera | Variable kinematic model for the analysis of functionally graded material plates[END_REF]. For all models, the accuracy of the displacements is excellent. Nevertheless, when k increases, the HSDT model drives to an error rate of more than 8% for the transverse shear stress. Layerwise type models are required to compute accurate transverse stresses for all cases. Moreover, in [START_REF] Carrera | Effects of thickness stretching in functionally graded plates and shells[END_REF], the authors have highlighted the importance to take into account the transverse stretching in the model. This phenomenon is included in the present approach.

[Table 6 about here.]

[Table 7 about here.]

The distributions of the in-plane and transverse shear stress through the thickness are also provided in Fig. 5 for k = 1 and k = 10. Results issued from the LD4 model are not given for clarity reason because the agreement with the present model is excellent. The satisfaction of the free boundary conditions is clear. We also notice the influence of the volume fraction exponent on the in-plane stress. A high stress gradient occurs on the top of the plate when k increases.

For this test case, the number of dofs of the LD4 model (316899 dofs) is 33 times greater than the VS-LD4 model.

[Figure 5 about here.]

Sandwich plate

In this section, a sandwich plate is considered. The test case is described as follows:

geometry: square sandwich plate (a = b = 1 m) with length-to-thickness ratio S = 4, 20, constituted of three layers, the thickness of the faces and the core is h f = 0.1h and h c = 0.8h, respectively. First of all, the importance of the numerical layers is highlighted to compute the transverse stresses. As illustrated in Fig. 6, a high discontinuity occurs at the interface between the core and the upper face when using one numerical layer per physical layer. It can be circumvented by introducing three numerical layers.

In the subsequent example, the results are given following this recommandation.

[Figure 6 about here.]

Different slenderness ratios and values of k are considered to assess the accuracy of the method. The results are summarized in Tab. 8. For both thin and thick cases, whatever the values of the volume fraction exponent, the results are in excellent agreement with the reference solution. The distribution of the in-plane and transverse shear stresses are also given in Fig. 7 for k = 1 and k = 10. The complex distribution of σ11 through the thickness is well represented. The traction free condition is fulfilled on the top and bottom surfaces.

[Table 8 about here.]

[Figure 7 about here.]

Multilayered FGM plate

In this section, the static behavior of a multilayered FGM plate is studied. The test case

given in [START_REF] Wu | An rmvt-based third-order shear deformation theory of multilayered functionally graded material plates[END_REF] is such as where L refers to the fiber direction, T refers to the transverse direction.

The FGM material is such that :

g (1) (z) = exp(k exp z+h/6 -h/3 ) for -h/2 ≤ z ≤ -h/6 g (3) (z) = exp(k exp z-h/6 h/3 ) for h/6 ≤ z ≤ h/2
with the material-property gradient index: k exp = 5 mesh: N x = N y = 32 results: The results are made nondimensional using:

ū = u 1 (0, b/2, z) E T hq 0 S 3 , σαα = σ αα (a/2, b/2, z) q 0 S 2 σ13 = σ 13 (0, b/2, z) q 0 S , σ33 = σ 33 (a/2, b/2, z) q 0
reference values exact elasticity results are given in [START_REF] Wu | An rmvt-based third-order shear deformation theory of multilayered functionally graded material plates[END_REF] and obtained from [START_REF] Wu | Three-dimensional static behavior of functionally graded magnetoelectro-elastic plates using the modified pagano method[END_REF].

This example seems to be a severe test case. The results will be compared with two Thirdorder Shear Deformation Theories with a displacement-based (denoted TSDT-Wu) and a RMVT-based (denoted RMVT-Wu) approach available in [START_REF] Wu | An rmvt-based third-order shear deformation theory of multilayered functionally graded material plates[END_REF]. Note that the transverse normal strain is assumed to be nil in these works. For the present approach, three numerical layers per physical layer is required to improve the accuracy of the transverse shear and normal stresses. The in-plane displacement and stress, and the transverse shear / normal stresses are shown in Fig. 8 and Fig. 9. It can be inferred from these figures that the VS-LD4 method is in excellent agreement with the exact solution. The RMVT approach is not sufficient to recover accurate results. Indeed, as previously observed in open literature, the influence of the stretching effect is of major importance.

[Figure 8 about here.]

[Figure 9 about here.]

Conclusion

In this article, an approach based on a variable separation method has been extended

k

  where E m = 70 GPa, E c = 380 GPa are the elastic modulus of the aluminium and alumina, respectively. k is the volume fraction exponent (k > 0). ν = 0.3 is considered as constant. mesh: N x = N y = 32 where N x and N y are the number of elements along the x and y directions, respectively. number of dofs: Ndof xy = 3.(3.N x .N y + 2(N x + N y ) + 1) and Ndof z = 12.N z + 3 are the number of dofs of the two problems associated with v i j and f i j respectively. N z is the number of numerical layers. So, the total number of dofs is Ndof xy + Ndof z . results displacements and stresses are made non-dimensional according to ū(z) = u 1 (0, b/2, z) 100 E c h 3 p 0 a 4 ; w(z) = u 3 (a/2, b/2, z) 10 E c h 3 p 0 a 4 ; σ11 (z) = σ 11 (a/2, b/2, z) h a p 0 ; σ13 (z) = σ 13 (0, b/2, z) h a p 0 ; σ33 (z) = σ 33 (a/2, b/2, z) 4.1.1. Convergence study A convergence study is carried out for two values of k corresponding to two different distributions of properties through the thickness. In-plane and out-of-plane mesh refinements are considered independently. First, Tab. 1 shows the convergence of the displacements and stresses with a fixed value of N z = 4. Classically, the convergence rate of the displacement is fast. Only N x = N y = 4 is sufficient. For stresses, 32 elements are needed to achieve

  also 8 numerical layers. Two slenderness ratios (thick and semi-thick configurations) with two values of k are considered. Results are summarized in Tab. 4 and Tab. 5. Results are in excellent agreement with the LD4 model for both displacements and stresses, regardless the values of S and k.The error rate is less than 0.1%. As the LD4 model gives quasi-3D results, the present results are very satisfactory.

  boundary conditions: simply-supported plate on all sides subjected to a bi-sinusoidal pressure p(x, y) = p 0 sin πx a sin πy b material properties: The two external faces are in aluminium at the bottom and in alumina at the top. bottom face sheets: isotropic material with E = 70 GPa, ν = 0.3 (aluminium) top face sheets: isotropic material with E = 380 GPa, ν = 0.3 (alumina) Core : FGM as in Section 4.1 mesh: N x = N y = 32 number of dofs: Ndof xy = 9603 and Ndof z = 12.N z + 3 = 111 number of dofs for LD4: N LW = 3.(4.N z + 1).(3.N x .N y + 2(N x + N y ) + 1) = 355311 results: non-dimensional results as in Section 4.1

geometry:

  square sandwich plate (a = b = 1 m) with length-to-thickness ratio S = 5, constituted of three layers [F GM/0 • /F GM], the thickness of the 3 layers are the same. boundary conditions: simply-supported plate on all sides subjected to a bi-sinusoidal pressure p(x, y) = p 0 sin πx a sin πy b material properties: E L = 25 GPa , E T = 1 GPa , G LT = 0.2 GPa , G T T = 0.5 GPa , ν LT = ν T T = 0.25

Figure 1 :Figure 2 :

 12 Figure 1: Geometry of the plate

Figure 3 :Figure 4 :Figure 5 :Figure 6 :

 3456 Figure 3: Convergence study -S=10 -one-layered -polynomial law -k = 1

Figure 7 :Figure 8 :Figure 9 :

 789 Figure 7: distribution of σ11 (left) and σ13 (right) along the thickness -S=4 -sandwich -polynomial lawk = 1/10

Table 1 :

 1 one-layered FGM plate -polynomial law -b = a -S = 10 -N z = 4

	k	Nx = Ny	ū(-h/4)	w(0)	σ11 (h/3)	σ33 (h/2)	σ13 (h/6)
		1	0.5790	0.4478	1.2655	2.3411	0.8612
		2	0.6461	0.5798	1.5361	1.6133	0.4561
		4	0.6438	0.5870	1.5256	1.2594	0.3123
	1	8	0.6436	0.5875	1.5121	1.0753	0.2672
		16	0.6436	0.5875	1.5077	1.0196	0.2550
		32	0.6436	0.5875	1.5066	1.0049	0.2519
		64	0.6436	0.5875	1.5063	1.0012	0.2511
		1	1.0066	0.8069	0.8090	4.1083	0.5376
		2	1.0874	0.9967	0.9205	2.2401	0.3270
		4	1.0844	1.0067	0.9069	1.4554	0.2531
	10	8	1.0842	1.0074	0.8983	1.1313	0.2313
		16	1.0842	1.0074	0.8957	1.0381	0.2256
		32	1.0842	1.0074	0.8950	1.0140	0.2241
		64	1.0842	1.0074	0.8948	1.0079	0.2238

Table 2 :

 2 one-layered FGM plate -polynomial law -b = a -S = 10 -N x = N y = 32

Zenkour06 [START_REF] Zenkour | Generalized shear deformation theory for bending analysis of functionally graded plates[END_REF] it refers to Generalized shear deformation theory. An ESL model including a sinus function with 5 unknowns is used. The transverse displacement is constant through the thickness. It satisfies the zero tangential traction boundary conditions on the surfaces of the plate. A Navier type solution is developed. Mantari12 [START_REF] Mantari | A new trigonometric layerwise shear deformation theory for the finite element analysis of laminated composite and sandwich plates[END_REF] a High-order Shear Deformation Theory (HSDT) model is used with 5 generalized unknowns (ESL model). The tangential stress-free boundary conditions on the plate surfaces is satisfied. A Navier-type solution is carried out.