Modeling of cylindrical composite shell structures based on the Reissner’s Mixed Variational Theorem with a variable separation method

This work deals with the modeling of laminated composite and sandwich shells through a variable separation approach based on a Reissner’s Variational Mixed Theorem (RMVT). Both the displacement and transverse stress fields are approximated as a sum of products of separated functions of the in-plane coordinates and the transverse coordinate. This approach yields to a non-linear problem that is solved by an iterative process, in which 2D and 1D problems are separately considered at each iteration. In the thickness direction, a fourth-order expansion in each layer is used. For the in-plane description, classical Finite Element method is used. Numerical examples involving several representative shell configurations (deep/shallow, thick/thin) are addressed to show the accuracy of the present method. It is shown that it can provide quasi-3D results less costly than classical LW computations. In particular, the estimation of the transverse stresses, which are of major importance for damage analysis, is very good.


Introduction
Composite shells are widely used in the industrial field (aerospace, automotive, marine, medical industries...) due to their excellent mechanical properties, especially their high specific stiffness and strength. For composite design, accurate knowledge of displacements and stresses is required. One way consists in considering the three-dimensional modelisation. However, due to the complexity of such numerical simulations, it is desirable to take advantage of the geometric ratios and to represent the problem as a two-dimensional model by referring to shell theories. There are two ways to define the approximation of the displacement field. A "pure shell" model can be considered in which the displacement is associated with the local curvilinear vectors, and strain and stress are deduced using the differential geometry [1]. Alternatively, the shell-like solid approach [2] is widely used to formulate shell Finite Element (FE), in particular in commercial software. In this case, the displacement vector is defined in the global cartesian frame. The jacobian matrix transformation is used to express strain and stress with respect to the reference frame defined • The Equivalent Single Layer Models (ESLM), to which belong the Classical Shell Theory (CST/Koiter) and First Order Shear Deformation Theory (FSDT/Nagdhi). The reader can refer to [5] for a detailed description of the assumptions on the strain field to derive different shell models. CST leads to inaccurate results for composites because both transverse shear and normal strains are neglected, see [6][7][8] for shallow laminated shells. FSDT is the most popular model due to the possibility to use a C 0 FE, but it needs shear correction factors and the transverse normal strain is still neglected (cf. [9][10][11][12]). So, Higher-order Shear Deformation Theories (HSDT) have been developed to overcome these drawbacks. Different kinematics including 7 [13,14] or 9 parameters [15] have been proposed. Reddy has derived a 5 parameter model starting from a third-order theory by considering the homogeneous top/bottom conditions for the transverse shear stresses [16]. Note also the variable kinematics approach based on the Carrera's Unified Formulation developed in [17,18]. In the ESLM context, a simple way to improve the estimation of the mechanical quantities consists in adding one zig-zag function, called Murakami's function, in the expression of the displacement. In this way, the slope discontinuity at the interface between two adjacent layers is introduced. It allows to describe the so-called zig-zag effect. It has been carried out in conjunction with the FSDT in [19] and [20] based on dedicated mixed formulation. It has been also used with the HSDT in [21] and [22] including 9 and 13 parameters, respectively. • The Layer-Wise Models (LWM), in which the expansion of the mechanical quantities is defined over each layer independently. Some of these works are based on a linear distribution of the in-plane displacements through each layer, without taking into account the transverse normal stress. The transverse displacement can be constant across the whole thickness, such as in [23,24], or in each layer separately, as in [25]. But, this type of approach fails to predict accurate transverse stresses, unless using dedicated post-processing steps [24]. Thus, higher-order approaches taking into account the transverse normal effect have been developed. The threedimensional constitutive law is used. Second-, third-and fourth-order expansions are discussed in [4,26,27]. In this framework, Kulikov [28] has developed the sampling surfaces method, see also the previously mentioned work [18]. In all the aforementioned LWM, the number of unknowns depends on the number of layers, which may thus affect the performance in terms of computational cost.
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. Starting from a LW expansion, enforcing the interlaminar continuity of the displacement and transverse shear stresses, as well as the homogeneous top/bottom conditions, the number of unknowns becomes independent of the number of layers. The deduced model can be derived from the CST [29], as well as from an HSDT with a third-order theory [30][31][32][33] or the sinus model [1]. The number of parameters of the resulting displacement model varies from 5 to 15. We can also cite two other approaches where only the displacement continuity [34] or the top/bottom conditions [16] are satisfied. While the above references propose models still relying on the 2D constitutive law, the papers [35,36] extend the zig-zag approach to include the transverse normal stress. 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 various theories and/or finite element applications can be found in [37][38][39][40][41][42][43][44][45] and more recently in [46].
Nevertheless, in the framework of the failure analysis of composite structures, involving the free edge effect for instance, the prediction of the interlaminar stresses is of major interest. In particular, the difficulty is to well-describe the interlaminar continuous transverse stresses. Most of the ESLM fail to represent these in the most critical cases, unless using a post-processing treatment [11,24,[47][48][49]. To overcome these drawbacks, alternative formulations to the displacement-based approach have been developed. On the one hand, several techniques have been devised to correct the transverse shear locking pathology affecting FSDT-based plate/shell elements, most of which can be stated from hybrid-mixed approaches [50]. For composites, an assumed strain approach has been adapted in a FSDT [51] or in a LWM [52]. On the other hand, assumed partial/complete stress field over the laminate thickness independently aims at increasing the accuracy of this one. Without considering the transverse normal stress, some authors [53,54] have adopted a partial hybrid stress approach based on the HellingerReissner Variational Principle (HRVP). Using a HSDT approach for displacements, Yong [53] has developed a generalized stress assumption for the transverse shear stresses only, whereas in-plane stresses are also involved in [54]. Alternative hybrid approaches take into account the transverse normal effect. The Fraeijs de VeubekeHuWashizu multifield variational principle [55,56] and the Reissner Mixed Variational Theorem [57] are used assuming the transverse shear stresses only. Note that Sgambitterra et al. [14] have introduced a mixedfield assumptions (γ α3 and σ 33 ) to derive a hybrid formulation enforcing the compatibility straindisplacement relations to be least-squares compatible through the shell thickness. As a complement to these hybrid methods, mixed formulations are addressed in conjunction with FSDT [58] (HRVP) or including the Murakami's zig-zag function [20] (Jing and Liao's functional [59]). An advanced method is proposed in [60] by considering all interlaminar stresses between two adjacent layers as primary variables and also a higher-order LW displacement.
In the same way, the partially Reissner's Mixed Variational Theorem (RMVT) assuming two independent fields for all displacement and transverse stress variables allows to ensure a priori interlaminar continuous transverse stress fields. The consideration of the transverse normal stress as an independent assumed variable seems to be important to obtain accurate distributions through the thickness of the shell (see e.g. [55]). The RMVT approach comes from the works of Reissner, see [61,62]. It was first applied for multilayered structures in [63] and then, in [64] with higher order displacement field and [65] with a Layerwise approach for both displacements and transverse stress fields. Afterwards, the approach was widely developed with a systematic approach based on the Carrrera's Unified Formulation to provide a large panel of 2D models for composite structures based on ESL and/or LW descriptions of the unknowns [66][67][68]. It has been also applied for shell structures in [69][70][71]. Note that the reader can refer to [27] for a survey on RMVT in this framework, and to [37,72] for a further discussion.
In the present work, the RMVT is considered because it is a natural variational tool for composite structures as it allows to formulate independent approximations for mechanical variables that are required to be continuous across the stacking direction. It is used in conjunction with a variable separation method, namely the Proper Generalized Decomposition (PGD). The main purpose consists in taking advantage of this promising approach to decrease the high computational cost of a LW approach. In fact, interesting features have been shown in the model reduction framework [73][74][75]. So, the aim of the present paper is to assess this particular representation of the unknowns in the framework of a mixed formulation to model cylindrical laminated and sandwich shells. Both displacements and transverse stresses are written under the form of a sum of products of bidimensional polynomials of (ξ 1 , ξ 2 ) and unidimensional polynomials of z. A piecewise fourth-order Lagrange polynomial of z is chosen across the thickness of each layer. As far as the variation with respect to the in-plane coordinates (ξ 1 , ξ 2 ) is concerned, a 2D eight-node quadrilateral FE is employed. With the PGD approach, each unknown function of (ξ 1 , ξ 2 ) is approximated using one degree of freedom (dof) per node of the mesh and the LW unknown functions of z are global for the whole shell. This process yields to two separate linear problems, i.e., a 2D problem in (ξ 1 , ξ 2 ) and a 1D problem in z, in which the number of unknowns is much smaller compared to a classical Layerwise approach. These two problems are solved successively within an iterative scheme. The interesting feature of this approach lies on the possibility to have a higher-order z-expansion and to refine the description of the mechanical quantities through the thickness without substantially increasing the computational cost. This is particularly suitable for the modeling of composite structures. Note that this method has been successfully applied to displacement-based plate/shell models in [76][77][78][79] and also to RMVT plate models in [80].
We now outline the remainder of this article. First, the shell definition and the differential geometry are recalled. The RMVT formulation is described and the separation of the in-plane and out-of-plane displacements/ transverse stresses is introduced. The principles of the PGD are precised in the framework of our study. The particular assumption on the assumed variables yields to a non-linear problem, that is solved within an iterative process. The FE discretization is also described and finally, numerical tests are addressed. A preliminary convergence study is performed. Then, the influence of classical shell assumptions on the strains and the number of numerical layers are studied. The approach is also assessed for a wide range of applications: deep/shallow shells, cross-ply/angle-ply configurations, different slenderness ratios and different degrees of anisotropy for a sandwich are considered. The accuracy of the results is evaluated by comparison with a 2D elasticity solution from [81], or results available in the open literature [79,82,83].

Shell definitions and differential geometry
A shell C with a middle surface S and a constant thickness e, see Fig. 1, is defined by [84]: where the middle surface can be described by a map from a parametric bidimensional domain as: In Fig. 1, the map describing the shell middle surface (in grey) and the local basis vectors are presented. The basis vectors a i are defined for a point on S and the basis vectors g i are defined for a generic point of the shell.
For a point on the shell middle surface, the covariant basis vectors defining the tangent plane to the middle surface are usually obtained as follows: where a 3 is the unit normal vector to the surface S, see The contravariant vectors are constructed from the covariant ones using the following equations: where δ α β is the Kronecker symbol. For a generic point of the shell, covariant basis vectors must be defined and we have where b β α is the mixed form of the second fundamental form. This basis g i , illustrated in Fig. 1, must be used to define quantities for any point of the shell. The form μ β α (z) introduced in Eq. (5) defines the transport from the shell middle surface to any point of the shell and is associated with the curvature variation along the thickness direction z of the shell. The inverse tensor of the mixed tensor μ β α is denoted m β α and is defined as: where we have introduced the determinant of the mixed tensor μ = det(μ β α ) = 1 − 2Hξ 3 + (ξ 3 ) 2 K and the invariants of the second fundamental form H = 1 2 tr(b β α ) and K = det(b β α ). Finally, the surface element dS and the volume element dC are given by: where a is the determinant of the first fundamental form a αβ . The geometry of a shell can also be defined using contravariant or mixed forms. Furthermore, covariant and contravariant differentiation involving Christoffel symbols are not detailed here and readers can refer to the book [84].

The definition of the strain field
For geometrically linear elastic analysis, the components of the strain tensor ε ij expressed in the contravariant basis a i are obtained through covariant differentiation, denoted |, as follows: The mixed tensor m β λ carries out the transport from any point of the shell to the shell middle surface, that is from g i to a i .

Constitutive relation
The stress tensor is obtained from the strain tensor using the constitutive equations. For this purpose, all these tensors must be referred to the covariant and contravariant basis vectors, a i and a i respectively, associated with the middle surface of the shell. In case of laminated shells composed of orthotropic plies, this reference frame ensures to consistently take into account the different material orientations of the layers. The tensorial relation is: It should be noted that the stress tensor is defined in the covariant reference frame, whereas the strain components are defined in the contravariant frame. If the frame is assumed to be orthonormal then covariant and contravariant components are equal, that is super-script and sub-script are interchangeable.
In the framework of the RMVT formulation, the Hooke's law has to be rewritten under a convenient mixed form.
Firstly, stresses σ and strains ε are split into two groups: where the subscripts n and p denote transverse and in-plane values, respectively. The shell can be made of NC perfectly bonded orthotropic layers. Using the previous assumptions and the separation between transverse and in-plane components, the three dimensional constitutive law of the kth layer is given by: In this equation, the subscript G indicates that the strain is issued from the geometrical relations, while H means that the quantities are calculated from the Hooke's law. Q ij are the three-dimensional stiffness coefficients of the layer (k). Finally, RMVT is associated to a mixed form of Hooke's law, which can be written as The assumed transverse stresses are denoted as σ nM . The superscript (k) and the coordinates (ξ , z) are omitted for clarity reason.
The relations between the coefficients of the classical Hooke's law Eq. (12) and the mixed one Eq. (13) are: The weak form of the boundary value problem The formulation of the problem is based on the Reissner's partially Mixed Variational Theorem [61], denoted RMVT. In this formulation, the Principle of Virtual Displacement is modified by introducing the constraint equation to enforce the compatibility of the transverse strain components. This term also depends on the assumed transverse stresses, see also [17,72]. Thus, the problem can be formulated as follows: Find u(M) ∈ U (space of admissible displacements) and σ nM such that where t is the prescribed surface forces applied on ∂C F . The body force is not considered in this expression.

Application of the proper generalized decomposition to the cylindrical shell
In this section, we develop the application of the PGD for shell analysis with a mixed formulation. This work is an extension of a previous study on composite cylindrical shell structures [79].

The cylindrical geometry
A cylindrical shell (see Fig. 2) is described using the following map: Following Eq. (3), the non-zero terms for the covariant and mixed forms are: and μ = 1 + z R .

The displacement and transverse stress field
Let us denote u i (ξ 1 , ξ 2 , ξ 3 = z) and σ i3 (ξ 1 , ξ 2 , ξ 3 = z) the curvilinear components of the displacement field and the transverse stresses respectively, associated with the contravariant basis vectors a i . Let and z be the bidimensional domain associated with the mid-surface of the shell [see Eq. (1)], and the unidimensional domain associated with the normal fiber, respectively. The displacement and transverse stress fields are constructed as the sum of N products of functions of in-plane coordinates and transverse coordinate (N ∈ N is the order of the representation) according to are defined in . The "•" operator is Hadamard's element-wise product.
In this paper, a classical eight-node FE approximation is used in and a LW description is chosen in z as it is particulary suitable for the modeling of composite structures.

The strain field for the cylindrical composite structure
The strain field in Eq. (8) is free of any approximated shell kinematics. These strain components are simplified using Eq. (17) and we recover the following relations: where covariant derivative becomes classical derivative for the case of a cylindrical shell. Equation (18) must be introduced at this level in order to obtain the compatible strain expansion along the normal coordinate z of the shell. So, we obtain: where the prime stands for the derivative with respect to z (f i = df i dz ), () ,α indicates the partial derivative with respect to ξ α (α = 1, 2).

The problem to be solved
The resolution of Eq. (15) 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 are the functions to be computed and u m , σ m nM are the associated known sets at iteration m defined by The test functions are with The test functions defined by Eqs. (26,27), the trial functions defined by Eqs. (23,24), and the mixed constitutive relation Eq. (13) are introduced into the weak form Eq. (15) to obtain the two following equations: For the present work, ∂C F is considered as the partial top or bottom surface of the shell, that is z = z F with z F = ±e/2.
From Eqs. (29) and (30), a coupled non-linear problem is derived whose solution is seeked by means of a fixed point method for simplicity reason. An initial function f (0) , f (0) σ is set, and at each step, the algorithm computes two new pairs (v (k+1) , f (k+1) ), These two equations are linear. The first one is solved on S, while the second one is solved on z . As in [80], the fixed point algorithm is stopped when the distance between two consecutive solutions is smaller than a fixed value ε 0 .

Finite element discretization
To build the shell finite element approximation, a discrete representation of the functions (v, τ, f , f σ ) must be introduced. In this work, a classical finite element approximation in and z is used. The elementary vector of degrees of freedom (dof) associated with one element e of the mesh in is denoted q vτ e , while the elementary dof vector associated with one element ze of the mesh in z is denoted q where and The matrices N ξ , B ξ , N z , B z , N σ ξ , N σ z contain the interpolation functions, their derivatives and the jacobian components, respectively.

Finite element problem to be solved on
For the sake of simplicity, the functions f (k) , f (k) σ which are assumed to be known, will be denotedf ,f σ , respectively, and the functions v (k+1) , τ (k+1) to be computed will be denoted v and τ, respectively. The strains and the assumed transverse stresses in Eq. (29) are defined as The variational problem defined on from Eq. (29) is Note that the units of C nn , C np , C pn and C pp are different. The introduction of the finite element approximation Eq. (31) in the variational Eq. (36) leads to the linear system where • q vσ is the vector of the nodal displacements / transverse stresses associated with the finite element mesh in , • K z (f ,f σ ) is the stiffness matrix obtained by summing the elements' stiffness

Finite element problem to be solved on z
For the sake of simplicity, the functions v (k+1) , τ (k+1) which are assumed to be known, will be denotedṽ,τ and the functions f (k+1) , f (k+1) σ to be computed will be denoted f , f σ . The strain in Eq. (30) is defined as The variational problem defined on z from Eq. (30) is and The introduction of the finite element discretization Eq. (31) in the variational Eq. (44) leads to the linear system where • q ff σ is the vector of degree of freedom associated with the F.E. approximations in z , • K ξ (ṽ,τ ) is obtained by summing the elements' stiffness matrices: is obtained by the summation of the elements' residual vectors given by

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. The geometrical transformation is based on an explicit map . A Gaussian numerical integration with 3 × 3 points is used to calculate the elementary matrices. Several static tests are presented with the aim of validating our approach and evaluating its efficiency. First, a convergence study is carried out to determine the suitable mesh for the further analysis. Then, the influence of the approximation on the factor 1/μ is addressed [1,85]. Three orders of expansion are considered. The influence of the numerical layers is also studied to illustrate the possibility to refine the transverse description of the displacements and the stresses. The present approach is assessed for deep/shallow and thick/thin cross-ply/angle-ply/sandwich shells to show the wide range of validity.
Unless otherwise mentioned, the numerical assessments are based on the test case proposed by Ren [81]. It concerns a semi-infinite simply-supported cylindrical shell submitted to a sinusoidal pressure, as detailed below:  Fig. 2). Boundary conditions: Simply-supported shell along its straight edges (transverse and tangential displacements are fixed on the (ξ 1 , ξ 2 ) mesh), sinusoidal pressure along the curvature: q(ξ 1 ) = q 0 sin πξ 1 Rφ Material properties: where L refers to the fiber direction, T refers to the transverse direction. Mesh: Only a quarter of the structure is meshed. The mesh is constituted of N x × N y elements in the ξ 1 and ξ 2 directions respectively. A space ratio is considered in these two directions (ratio between the size of the larger and the smaller element). Numerical layers: N z is the total number of numerical layers. Number of dofs: Ndof xy = 6(3N x N y + 2(N x + N y ) + 1) and Ndof z = 24 × αNC + 6 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 . Results: The results are made nondimensional using: q 0 Reference values: The 2D exact elasticity results are obtained as in [81].
Few iterations are required to reach the convergence of the fixed point algorithm (cf. Section ). Moreover, for this test case, only one couple is needed to obtain the solution.
The present approach, denoted VS-LM4, is compared with results available in open literature: 33   LM4: It refers to the systematic work of Carrera and his "Carrera's Unified Formulation" (CUF), see [37,72,86]. A LayerWise model based on a RMVT approach where each component is expanded until the fourth-order, is given. 24NC +6 unknown functions per node are used in this kinematic. VS-LD4: It refers to the work on the proper generalized decomposition with a spatial separation between (x, y) and z. A fourth-order expansion is used for the 1D problem associated to the z-direction. The formulation is based on a displacement approach (see [79]).

Convergence study
Firstly, a convergence study is performed to determine the suitable mesh to be used for the following test cases. A one-layered shell with S = 40 and φ = π/3 is considered. For the mesh, a space ratio is taken equal to 50 in the ξ 2 -direction, and a regular one is used in ξ 1 -direction. The results issued from the Navier approach which is the asymptotic value of the present model, are also given. The latter are quite similar to the reference 2D elasticity solution given by Ren [81]. Table 1 shows that a 26 × 10 mesh (see Fig. 3) is adequate to model the structure. These results can be considered as converged values. The maximum error rate remains less than 0.9% for both displacements and stresses with respect to the reference solution. This mesh is used for the following investigations.

Influence of the expansion order for the factor 1/μ
In this section, the influence of the approximation on the factor 1/μ = 1/(1 + z/R) is addressed for very thick to thin homogeneous shells and for deep or shallow shells, respectively. For this purpose, the results involving the zero, one and two-order expansions are compared in Tables 2 and 3. The zero-order approximation (thin shell hypothesis or Love's hypothesis) seems to be suitable only for the thin structure with S ≥ 40. For S ≤ 10, the order of the development of 1/μ should be increased to improve significantly the accuracy of the solution, including both the displacements and the stresses. Nevertheless, the use of the second-order expansion has a significant influence only for the thick shell.
Considering different opening angles of the structure in Table 3, it can be inferred from these results that the use of the two-order expansion allows us again to decrease significantly the error rate for both displacements and stresses. Nevertheless, the influence on the transverse normal stress does not appear for the shallow shells (φ ≤ π/3).
It can be also noticed that the error rate on the transverse shear and normal stresses remains high for both the thick structures and the shallow shells (φ ≤ π/3), regardless of the degree of the approximation of 1/μ. Even for φ = π/8, the error rates associated to the displacement u and the normal stress σ 11 are still high. Thus, a residual error remains. This problem is investigated in the following section.

Influence of the number of numerical layers
The influence of the number of numerical layers is studied for the homogeneous case.
In Tables 4 and 5, four numerical layers per physical layer are considered for moderately thick to very thick cases and shallow shells, respectively. These configurations are the most representative ones. First, it can be observed that the accuracy of results with different opening angles (Table 5) are improved when compared with Table 3, regardless of the degree of approximation. In particular, the effect on the stresses is more pronounced. Then, the use of the second order expansion of 1/μ drives to a very good agreement with the reference solution (Tables 4, 5). The maximum error rate is 1%, even for the very thick case. Finally, a refinement in the z-direction in conjunction with a high-order expansion of 1/μ is required to obtain accurate results for the shallow shells (φ ≤ π/3) or for the moderately to very thick cases (S ≤ 10). In Table 6, the suitable number of numerical layers is given for different slenderness ratios so as to obtain an error rate of about 1% for φ = π/8 (the most critical test). As  expected, it can be seen that the number of numerical layers increases with the thickness of the structure. Even considering a structure with only one physical layer, the refinement of the description of mechanical quantities is required as their distributions through the thickness could be very complex for a shell structure. For illustration, the distributions of the displacementū and the stressesσ 11 ,σ 13 through the thickness are provided in Fig. 4. An oscillating behavior occurs for the displacement and the in-plane stress, which is quite different than the plate case. Moreover, an asymmetrical distribution can be observed. The maximum value of the transverse shear stress is obtained near the top of the homogeneous shell.
For the present approach, it should be also noted that the additional computational cost inducing by the z-refinement is negligible as only the number of dofs of the 1D problem increases (Ndof z = 24 × αNC + 6, see Table 6). The size of the 2D problem remains unchanged. This is a main difference with respect to a LW approach where the total number of unknowns would be Ndof xy × Ndof z = Ndof xy × (24 × αNC + 6).
Thus, despite the simplicity of the stacking sequence, it is observed that the use of a higher-order theory involving numerical layers is required for thick or shallow shells.

Cross-ply test case
Firstly, a thick three-layer [0 • /90 • /0 • ] shell is considered, referring to the test case proposed by Ren [81] described in the preliminaries of the present section. Only distributions of the transverse normal and shear stresses through the thickness are provided for deep and shallow shells (see Fig. 5), as those are the most difficult to obtain. Only 1 couple is built and only one element per layer is used for the problem in z . It is inferred from this figure that the present model gives excellent results when compared to the exact solution. It can be noticed that the top/bottom surface conditions are fulfilled. Moreover, the maximum value of the transverse normal stress inside the structure increases for the deepest shell. It becomes greater than the applied load value on the top surface. For further illustrations, a 24-layer case with S = 4 and φ = π/3 is addressed. The distributions of displacements and stresses are given in Figs. 6 and 7. This example shows the capability of the approach to provide accurate results for a significant number of layers. It is to be noted that the agreement of the present solutions with the exact one is very good. To achieve a solution with the same accuracy, it is needed to use a LW approach. To illustrate the interest of the present method in terms of computational cost, we can compare the number of dofs involved for each model. For the present one, the 2D and 1D problems imply Ndof xy = 5118 and Ndof z = 582, respectively, while, a LW approach with a fourth-order expansion induces Ndof LW = 496, 446.

Angle-ply test case
In this section, a test case proposed by Bhaskar [82] involving a three-layer [45 • /−45 • /45 • ] with S = 4 and φ = 1 is given. The material properties and the loads are the same as the Ren's test. Numerical results for displacements and stresses are summarized in Table  7 and are compared with the exact solution. Again, the accuracy of the results is very satisfactory. The continuity of the transverse stresses is naturally fulfilled owing to the mixed formulation. We also note that the free boundary conditions are satisfied.

Bending analysis of a sandwich shell under a sinusoidal pressure
The approach is now assessed on a sandwich shell. The test is based on the Ren test and is proposed by Carrera in [83]. It is detailed below: Geometry: Cylindrical shell with R = 10, S = 4. The thickness of each face sheet is e 10 . The panel is supposed infinite along the x 2 = ξ 2 direction. Boundary conditions: Simply-supported shell along its straight edges, subjected to a sinusoidal pressure along the curvature: q(ξ 1 ) = q 0 sin πξ 1 Rφ .
Reference values: the LM4 results are given in [83].
The parameter α allows us to define different face-to-core stiffness ratios, i.e.
E face E Tcore = 73.10 5 and E face E Tcore = 73.10 7 for α = 1 and α = 10 −2 , respectively. This test case is discriminating for the assessment of composite models. In [83], the authors have shown that LayerWise models are needed to obtain accurate distributions of the transverse stresses through the thickness for this soft core sandwich shell. Moreover, the use of a mixed formulation with an order expansion greater than two, also increases the accuracy of the results. The present approach can be assessed for this discriminating test case. The distributions through the thickness of the transverse displacement and the transverse shear/normal stresses are given in Figs. 8 and 9. It can be noticed that the results are in very good agreement with the LM4 model. The localisation of the transverse stress in the faces of the sandwich is well-represented. Very high values are reached, contrary to the plate case where the maximum value occurs on the loading surface. The results issued from the VS-LD4 model are also provided for further comparison. For both cases, it shows that the transverse normal stress is very difficult to obtain by a displacement based formulation even with a higher-order theory.

Bending analysis of laminated shells under a constant pressure
In this section, the Ren configuration involving three layers is considered with a constant global pressure. To assess the present mixed model, the results of the displacementbased model with variable separation, denoted VS-LD4, are given. In this simulation, four numerical layers per physical layer with a fourth-order expansion through the thickness is  Fig. 10 as these are the most difficult quantities to compute. It can be inferred from this figure that five couples are needed to obtain accurate results. Nevertheless, the convergence rate is higher for the transverse normal stress. For this later, only two couples allow us to recover the load value on the top surface of the shell.

Conclusion
In this paper, a variable separation method in the framework of Reissner's Mixed Variational Theorem is proposed for the modeling of laminated composite and sandwich shells. A 8-node FE for the in-plane unknown approximation and a fourth-order LW description for the thickness unknown approximation are used. In this formulation, all interface conditions are exactly satisfied. The approach has been assessed through different benchmarks proposed in the open literature. The influence of classical strain assumptions (approximation of 1/μ) is discussed. We have shown the importance of the expansion of this term for thick shells. At the same time, the accuracy of the results could be also improved by using numerical layers in each physical layer without increasing significantly the computational cost. In fact, the number of layers has no influence on this cost as only the cost of the 1D problem is affected. This is particularly interesting in the framework of a mixed approach, where the number of unknowns involving both displacements and transverse stresses becomes very important in a classical LW method. Comparisons with exact reference solutions, results available in open literature have shown the very good accuracy of the method for a wide range of applications. Deep/shallow, very thick/thin shell structures with various stacking sequences and high anisotropy configurations can be modeled with an excellent accuracy. So, the present work can provide quasi-3D results avoiding expensive 3D FEM or LW computations. Therefore, this method seems to have very attractive features.