Separated representations of 3D elastic solutions in shell geometries

The solution of 3D models in degenerated geometries in which some characteristic dimensions are much lower than the other ones -e.g. beams, plates, shells,...- is a tricky issue when using standard mesh-based discretization techniques. Separated representations allow decoupling the meshes used for approximating the solution along each coordinate. Thus, in plate or shell geometries 3D solutions can be obtained from a sequence of 2D and 1D problems allowing fine and accurate representation of the solution evolution along the thickness coordinate while keeping the computational complexity characteristic of 2D simulations. In a former work this technique was considered for addressing the 3D solution of thermoelastic problems defined in plate geometries. In this work, the technique is extended for addressing the solution of 3D elastic problems defined in shell geometries. The capabilities of the proposed approach are illustrated by considering some numerical examples involving different degrees of complexity, from simple shells to composite laminates involving stiffeners. The analyzed examples prove the potentiality and efficiency of the proposed strategy, where the computational complexity was found evolving as reported in our former works, proving that 3D solutions can be computed at a 2D cost.


Background
Plates and shells are very common in nature and thus they inspired engineers that used both from the very beginning of structural mechanics. Shells offer a diversity of possible shapes and geometries, some of them with simple curvature and most of them with double curvature. Many times they are assembled in complex structural systems, in many applications they contain many stiffeners as in the case of aircraft fuselages.
In general the design of such structural elements requires the calculation of stresses, strains and displacements for the design loads. Strains and stresses are related by the so-called constitutive law. The simplest one consists of the linear elasticity. Despite its simplicity many structures are designed for working precisely within the elastic domain. Other designs require considering more complex behaviors (e.g. non-linear elasticity due to material or geometrical non linearities, elastoplastic behaviors usually encountered in material forming -forging, bending, ... -, or complex multiphysics behaviors as the ones encountered in composites manufacturing processes implying http://www.amses-journal.com/content/1/1/4 change of phases, crystallization, polymerization, ... coupled with rich thermomechanical mechanisms).
Design problems always involve the solution of a set of partial differential equations in the degenerate domain of the plate or the shell with appropriate initial and boundary conditions. These domains are degenerated because one of its characteristic dimensions (the thickness in the present case) is much lower that the other characteristic dimensions. We will understand the consequences of such degeneracy later. When analytical solutions are neither available nor possible because the geometrical or behavior complexities, the solution must be calculated by invoking any of the available numerical techniques (finite elements, finite differences, finite volumes, methods of particles, ...).
In the numerical framework the solution will be only obtained in a discrete number of points, usually called nodes, distributed in the domain. From the solution at those points, it can be interpolated at any other point in the domain. In general regular nodal distributions are preferred because they offer the best accuracy. In the case of degenerated plate or shell domains one could expect that if the solution evolves significantly in the thickness direction, a large enough number of nodes must be distributed along the thickness direction to ensure the accurate representation of the field evolution in that direction. In that case, a regular nodal distribution in the whole domain will imply the use of an extremely large number of nodes with the consequent impact on the numerical solution efficiently.
When simple behaviors and domains were considered semi-analytical models can be considered [1]. For addressing more complex scenarios plate and shell theories were developed allowing, through the introduction of some hypotheses, reducing the 3D complexity to the 2D related to the problem now formulated by considering the in-plane coordinates. The use of these theories have been extended gradually for addressing larger and more complex geometries (anisotropic laminates, ...) and behaviors.
There are thousand of papers concerning the proposal and application of plate and shell models (the interested reader can refer to the recent reviews [2,3] and the references therein). Some models are based on the introduction of kinematic hypotheses in the thickness (e.g. [4] among many others). Transverse shear can be also taken into account [5]. Recent zig-zag representations [6,7], layer-wise models [8][9][10] and solid-shell approaches [11,12], allow addressing accurately more complex scenarios, by increasing the computational complexity slightly. Stiffeners require an appropriate coupling of beam and shell models in order to perform calculations at a moderate computational cost [13].
However, as soon as richer physics are involved in the models and the considered geometries differ of those ensuring the validity of the different reduction hypotheses, efficient simulations are compromised. For example in composites manufacturing processes of large parts many reactions and thermal processes inducing significant evolutions on the thermomechanical fields in the thickness occur. These inhomogeneities are at the origin of residual stresses and the associated distortion of the formed parts [14].
In these circumstances as just indicated the reduction from the 3D model to a 2D simplified one is not obvious, and 3D simulations appear many times as the only valid route for addressing such models, that despite the fact of being defined in degenerated geometries (plate or shell) they seem requiring a fully 3D solution. However in order to integrate such calculations (fully 3D and implying an impressive number of degrees of freedom) in usual design procedures, a new efficient (fast and accurate) solution procedure is needed. http://www.amses-journal.com/content/1/1/4 The Saint Venant's principle was extensively used in the Ladeveze's works for defining elegant and efficient 3D simplified models [15]. This technique was then generalized to dynamics [16]. This technique allowed significant reduction of computational complexity.
Later, a new discretization technique based on the use of separated representations was proposed for addressing space-time nonlinear models [17] and then it was generalized for defining general separated representations of solutions involving conformational coordinates [18], space and time and even parameters considered as extra-coordinates. The interested reader can refer to the recent reviews [19][20][21][22] and the references therein.
A direct consequence was the separated representations involving the space coordinates. Thus in plate domains an in-plane-out-of-plane decomposition was proposed for solving flow problems in laminates [20], then for solving thermal problems in extruded geometries [23], elasticity problems [24] and coupled multiphycisc problems [25]. In those cases the 3D solution was obtained from the solution of a sequence of 2D problems (the ones involving the in-plane coordinates) and 1D problems (the ones involving the coordinate related to the plate thickness).
It is important emphasizing the fact that these approaches are radically different to standard plate and shell approaches. We proposed a 3D solver able to compute the different unknown fields without the necessity of introducing any hypothesis. The most outstanding advantage is that 3D solutions can be obtained with a computational cost characteristic of standard 2D solutions. Moreover, as noticed in [24] no locking effects were found, possibly because the fully 3D solution accomplished.
In this work we will generalize the just referred approach considered in the case of plate domains for calculating the fully 3D solution of the elastic problem in shell domains. The 3D solution will be calculated again from the solution of a sequence of 2D and 1D problems thanks to the in-plane-out-of-plane separated representation.
It is important to note that in this paper we are not addressing a new shell modeling, and by this reason we do not need neither establishing a precise state of the art on shell theories nor comparing our approach with the solutions obtained by using shell models. As we are proposing a new procedure for calculating 3D solutions (keeping a computational complexity characteristic of 2D solution procedures) we will compare our solutions with the ones obtained by considering the fully 3D elastic solution in the shell geometries computed with standard 3D solvers (e.g. finite elements).
Before generalizing the technique proposed in [24] for treating elastic problems defined in shell domains we are summarizing it.
In-plane-out-of-plane separated representation of elastic problems defined in plate domains We proposed in [24] and original in-plane-out-of-plane decomposition of the 3D elastic solution in a plate geometry. The elastic problem was defined in a plate domain = ×I with (x 1 , x 2 ) ∈ , ⊂ R 2 and x 3 ∈ I, I =[ 0, H] ⊂ R, being H the plate thickness. The separated representation of the displacement field u = (u 1 , u 2 , u 3 ) reads: where P i k , k = 1, 2, 3, are functions of the in-plane coordinates (x 1 , x 2 ) whereas T i k , k = 1, 2, 3, are functions involving the thickness coordinate x 3 . In [24] we compared the first modes of such separated representations with the kinematic hypotheses usually considered in plate theories. Similar behavior was noticed in the case of elastic solutions in shell domains with respect to classical shell theories.
Expression (1) can be written in a more compact form by using the Hadamard product: where vectors P i and T i contains functions P i k and T i k respectively. Because neither the number of terms in the separated representation of the displacement field nor the dependence on x 3 of functions T i k are assumed a priori, the approximation is flexible enough for representing the fully 3D solution, being obviously more general than theories assuming particular a priori evolutions in the thickness direction x 3 .
Let's consider a linear elasticity problem on a plate domain = × I. The weak formulation reads: where K is the generalized 6 × 6 Hooke tensor, f d represents the volumetric body forces while F d represents the traction applied on the boundary N . The separation of variables introduced in Eq. (1) yields the following expression for the derivatives of the displacement components u i , i = 1, 2, 3: for j = 1, 2; and from which we can obtain the separated vector form of the strain tensor : Depending on the number of non-zero elements in the K matrix, the development of (u * ) T · K · (u) involves different number of terms, 21 in the case of an isotropic material and 41 in the case of general anisotropic behaviors. http://www.amses-journal.com/content/1/1/4 The separated representation constructor proceeds by computing a term of the sum at each iteration. Assuming that the first n−1 modes (terms of the finite sum) of the solution were already computed, u n−1 (x 1 , x 2 , x 3 ) with n ≥ 1, the solution enrichment reads: where both vectors P n and T n containing functions P n i and T n i (i = 1, 2, 3) depending on (x 1 , x 2 ) and x 3 respectively, are unknown at the present iteration. The test function u * reads u * = P * • T n + P n • T * .
The introduction of Eq. (7) into (3) results in a non-linear problem. We proceed by considering the simplest linearization strategy, an alternated directions fixed point algorithm, that proceeds by calculating P n,k from T n,k−1 and then by updating T n,k from the just calculated P n,k where k refers to the step of the non-linear solver. The iteration procedure continues until convergence, that is, until reaching the fixed point P n,k • T n,k − P n,k−1 • T n,k−1 < , that results in the searched functions P n,k → P n and T n,k → T n . Then, the enrichment step continues by looking for the next mode P n+1 • T n+1 . The enrichment stops when the model residual becomes small enough.
When T n is assumed known, we consider the test function u given by P •T n . By introducing the trial and test functions into the weak form and then integrating in I because all the functions depending on the thickness coordinate are known, we obtain a 2D weak formulation defined in whose discretization (by using a standard discretization strategy, e.g. finite elements) allows computing P n .
Analogously, when P n is assumed known, the test function u is given by P n • T . By introducing the trial and test functions into the weak form and then integrating in because all the functions depending on the in-plane coordinates (x 1 , x 2 ) are at present known, we obtain a 1D weak formulation defined in I whose discretization (using any technique for solving standard ODE equations) allows computing T n .
The problems related to the solution of functions P n and T n are defined in Appendix A. As discussed in [24] this separated representation allows computing 3D solutions while keeping a computational complexity characteristic of 2D solution procedures. If we consider a hexahedral domain discretized using a regular structured grid with N 1 , N 2 and N 3 nodes in the x 1 , x 2 and x 3 directions respectively, usual mesh-based discretization strategies imply a challenging issue because the number of nodes involved in the model scales with N 1 · N 2 · N 3 , however, by using the separated representation and assuming that the solution involves N modes, one must solve about N 2D problems related to the functions involving the in-plane coordinates (x 1 , x 2 ) and the same number of 1D problems related to the functions involving the thickness coordinate x 3 . The computing time related to the solution of the one-dimensional problems can be neglected with respect to the one required for solving the two-dimensional ones. Thus, the resulting complexity scales as N · N 1 · N 2 . By comparing both complexities we can notice that as soon as N 3 N the use of separated representations leads to impressive computing time savings, making possible the solution of models never until now solved, and even using light computing platforms.
In [24] we considered the simplest approximations of functions involving the in-plane coordinates P i k (x 1 , x 2 ) by considering bilinear quadrilateral finite elements and piecewise http://www.amses-journal.com/content/1/1/4 linear 1D elements for approximating functions involving the thickness coordinate x 3 , . Richer approximations were analyzed in [26].
In the present work we are generalizing the just described separated representation for solving 3D models defined in shell domains.

3D elastic problem in a shell domain: Shell representation
The shell domain S , assumed with constant thickness, can be described from a reference surface X, that in what follows will be identified to the shell middle surface but that in the general case could be any other one, parametrized by the coordinates ξ , η, that is X(ξ , η), where: Being n the unit vector normal to the middle surface, the shell domain S can be parametrized from: whose expression is given in the Appendix B (Eqs. (73) and (74) that involve Eqs. (42), (43) and (49)). The inverse transformation (x 1 , x 2 , x 3 ) → (ξ , η, ζ ), described byF −1 is also given in the appendix (Eq. (75)).

3D elastic problem in a shell domain: Weak form
The weak form of the elastic problem defined in the shell domain S writes: Now we are considering the coordinates transformation introduced in the previous section and deeply developed in the appendix, mapping The geometric transformation requires to transform the differential operator as well as the different volume and surface elements. Knowing that under the small displacements and strains assumption, the strain tensor consists of the symmetric part of the gradient of displacement tensor, i.e.
it can be transformed taking into account the transformation of the gradient differential operator where ∇ ξ (·) denotes the gradient in the parametric space. http://www.amses-journal.com/content/1/1/4 The volume element involved in the integral in S writes according to Eq. (70) with a is the determinant of the metric tensor related to the middle surface mapping (46) and H and K the curvatures given by Eqs. (56) and (57) respectively. Finally the integrals applying on the domain boundary (where tractions apply) are transformed according to Eq. (68).

3D elastic problem in a shell domain: In-plane-out-of-plane separated representation
With the weak form defined in = ×I the situation is quite similar to the one encountered in the analysis of elastic problems in plate geometries, that was addressed in [24] and we just summarized in section 'Background'. We could perform an in-plane-out-of-plane separated representation of the displacement field, similar to (1) but now involving the coordinates (ξ , η, ζ ) or in a more compact form As explained in section 'Background' the construction of such a separated representation is performed sequentially, thus assuming known the solution at iteration n − 1, the solution at iteration n is sought as By introducing (17) in the weak form and using the alternated directions fixed point algorithm we can calculate P n (ξ , η) by assuming T n (ζ ) known and then updated T n (ζ ) from the just calculated P n (ξ , η). The iteration continues until reaching the convergence (the fixed point) that determines both functions P n (ξ , η) and T n (ζ ).
However, the decomposition in a problem defined in for calculating function P n (ξ , η), obtained by integrating the weak form in I, and in another problem defined in I for calculating function T n (ζ ), obtained by integrating the weak form in , requires the separated representation of all the operators, variables, coefficients and functions involved in the weak form.
For the displacement (the trial u and the test u * displacements) we just indicated the separated in-plane-out-of-plane representation. This representation allows defining a separated representation of the associated strain tensors (u) and (u * ) as illustrated in Eq. (6) but for this purpose we must define a separated representation of the transformation gradient involved in Eq. (13)F −1 .
In the case of small thickness and curvaturesF −1 can be approximated, according to (73), from: where b n and F only depend on the middle surface parametrization. http://www.amses-journal.com/content/1/1/4 Thus, we can define a direct separated representation ofF −1 : where U 3 is the 3 × 3 matrix with unit components, that is: When this approximation (small thickness and curvatures) does not apply we should consider the separated representation ofF −1 by applying a singular value decomposition -SVD-.
The elasticity tensor and the applied forces must be also expressed in a separated form. Most of time this decomposition is direct and in the more complex cases the application of a singular value decomposition allows such decomposition.
Finally the volume and surface elements must be also written in a separated form. In the case of the volume element dx expressed by Eq.
the following separated representation can be assumed with From Eq. (68) we can expect a similar decomposition of the surface element.

3D elastic problem in a shell domain: Composite laminates
In this section we consider composite laminates composed of P anisotropic (generally orthotropic) plies of thickness h (for the sake of simplicity and without loss of generality we assume all the plies with the same thickness). We consider the mechanical behavior of each ply i expressed by its elasticity tensor K i whose form is quite simple when it is expressed in the basis related to its mechanical principal directions. Now, for the sake of simplicity we are also considering a local orthonormal basis defined at each point on the middle surface (t 1 , t 2 , n), where the normal vector n is the one previously considered, and the tangent vectors t 1 and t 2 to the middle surface at the considered location can be chosen arbitrarily under the constraints: t 1 · t 2 = 0, t 1 · n = 0, t 2 · n = 0, t 1 = 1 and t 2 = 1. http://www.amses-journal.com/content/1/1/4 If we define R i as the rotation tensor allowing the expression of the elastic tensor in the orthonormal local system (t 1 , t 2 , n) and Q i the one allowing to express finally its behavior in the cartesian basis (both expressed using the Voigt notation) it results If we assume that the elastic properties of each ply are constant along the ply thickness, then K i only depends on the in-plane coordinates. Thus, the laminate elastic tensor can be written in the separated form or using the Hadamard notation with k i (ζ ) = χ i (ζ ) · U 6 (U 6 being the 6 × 6 matrix with unit components) and χ i (ζ ) the characteristic function related to the i-ply:

Strategy verification
First we verify the solution computed using the proposed strategy by comparing it with the exact solution of a linear and isotropic elastic problem defined in an infinite tube subjected to an internal pressure. Figure 1   We consider the middle surface representation ⎧ ⎪ ⎨ ⎪ ⎩ with R = 2.5, ξ ∈[ 0, 2 · π · R] and η ∈[ 0, 1] (we are considering an arbitrary tube length due to the plane state of strain). The shell geometry is defined by introducing the thickness coordinate ζ ∈[ −0.5, 0.5].
The exact elastic solution only involves radial displacements, u r , depending on the radial coordinate r: In the numerical solution the domain (ξ , η) ∈ =[ 0, 2 · π · R] ×[ 0, 1] is discretized by using a regular mesh of rectangular bilinear elements whereas the thickness interval ζ ∈ I =[ −0.5, 0.5] is discretized by using standard one-dimensional linear elements (the simplest choices, but higher order approximations can be used in and/or I). The displacement is constrained at some appropriate locations in order to avoid rigid body movements. A convergence analysis is performed by considering the 3 different meshes reported in Table 1 where N ξ , N η and N ζ are the number of elements along the angular, axial and thickness coordinates respectively.
The problem is also solved numerically by using a 3D finite element discretization operating on the domain by considering regular meshes of 3D trilinear elements, compatible with the description given by the meshes in Table 1 consisting of 32, 128 and 384 3D-elements respectively.
Radial displacements calculated with finite elements and separated representations for equivalent meshes are compared with the exact solution (30) in Figure 2. A superiority of the separated representation solution with respect to the standard finite element solution can be noticed. This superiority can be associated to a better representation of the metrics and curvature.
Because the solution of the considered problem only involves radial (thickness) displacements depending on the radial coordinate, a 3D elastic solution seems to be the most appropriate and simplest choice, more natural than using computational shell theories. This problem can be solved easily with quadratic finite elements or isogeometric ones by only one element across the thickness. Also with axisymmetric representation only one 8 nodes element gives a nearly perfect solution. However here we prefer using linear finite elements in order to compare with the linear representation of functions involved in the PGD representation here considered, even if we could use higher order representations within the PGD framework. The objective here is proving the convergence and comparing the simplest finite elements and PGD formulations.
As the solution only depends on the radial direction one could expect to capture the solution with a single mode of the separated representation. This expectation is confirmed by the numerical solution that converges with a single term in the decomposition. Figures 3, 4 and 5 depict the computed mode for each component of the displacement field (P 1 1 (ξ , η), T 1 1 (ζ )), (P 1 2 (ξ , η), T 1 2 (ζ )) and (P 1 3 (ξ , η), T 1 3 (ζ )) associated with the decomposition (15) respectively, all them associated with the finest mesh (MESH 3) reported in Table 1. This numerical example served to validate the proposed approach based on the use of separated representations for solving 3D elastic problems in shell domains. It is important to notice that we compute the fully 3D solution with a prescribed accuracy. Thus, the proposed procedure should be viewed as a 3D solver that due to the separated representation of the 3D fields involved in the model allows a reduction of the computational complexity to the one characteristic of 2D solutions.

Elliptic tube
In order to address a problem whose solution consists of more than a single mode, we consider in this section a slightly different problem. It concerns a tube of infinite length but now with an elliptical cross section and again subjected to an internal pressure. The surface representation in given by: where R = 2.5, ξ ∈[ 0, 2 · π · R] and η ∈[ 0, 1] (we are considering an arbitrary tube length due to the plane strain assumption). The tube thickness is defined from ζ ∈[ −0.5, 0.5]. Material parameters are the same than in the previous test. We consider the same discretization than in the previous example, bilinear rectangles in and standard one-dimensional linear elements in the thickness (again the simplest choice).

Tube with finite length
Now in order to remove the plane strain hypothesis we consider a circular tube of finite length (L = 10) with R = 5 and unit thickness, clamped at both ends and again subjected to an internal pressure. The most significant modes in the present case are depicted inFigures 19, 20, 21 and 22 and the 3D solution in Figure 23.

Twisted tourus
In this section we consider a quarter of a twisted torus of thickness h = 0.2 depicted in Figure 24. Again we consider the linear and isotropic elastic problem related to an internal unit pressure and we constraint the displacement of both bases with respect to their normal directions. Again the displacement is constrained at some points for avoiding rigid body displacements. We consider the middle surface representation ⎧ ⎪ ⎨ ⎪ ⎩ X 1 = (r + R · cos(ξ + 2 · η) + e · cos(6 · ξ )) cos(η) X 2 = (r + R · cos(ξ + 2 · η) + e · cos(6 · ξ )) sin(η) X 3 = R · sin(ξ + 2 · η) + e · sin(6 · ξ ) where R = 2, r = 4 and e = 0.15. The parametric space consists of = × I, with The 2D domain is discretized by using a regular mesh consisting of 64 × 30 2D bilinear finite elements whereas the thickness interval I is discretized from 50 one-dimensional linear elements (again higher order approximation are possible). Figure 25 depicts the first component of the strain 11 . Figure 26 depicts the components of the displacement along the normal direction at the different locations shown in Figure 25. It can be noticed that the solution evolves significantly from one point to other. The separated representation captures this rich evolution justifying the interest of performing a fully 3D solution as soon as the computational complexity can be reduced to the one characteristic of 2D solutions.

Scordelis-Lo roof
In this section we address a critical scenario to check convergence and locking issues: the Scordelis-Lo roof depicted in Figure 27. Other than the symmetry conditions indicated in that figure we constraint the first and third components of the displacement, i.e. u 1 = u 3 = 0 on the end sections (in particular BC in the domain of calculation). The problem material and geometrical parameters are the following E = 3 10 6 , ν = 0, a = 600, r = 300, e = 3 and α = π/3. The shell is subjected to a surface load of 0.625 that induces a vertical displacement of 3.6288 at position A in Figure 27 [27].
The computed results obtained by using the in-plane-out-of-plane separated representations are again compared to the ones obtained by using standard 3D linear finite elements (better choices related to higher order approximations exist, but here again we focus in the simplest choice). The comparison is depicted in Figure 28. From this comparison we can conclude that our approach based on a separated representation converges and that it exhibits similar convergence behavior than standard linear 3D finite elements.

Multilayered tube
We consider now a composite tubular part composed of P = 8 plies of thickness 0.125 mm with a stacking sequence [ 90, 45, 0, −45] s in the local orthonormal basis (t 1 , t 2 , n), with t 1 , t 2 and n defining the axial, the circumferential and the thickness direction respectively as depicted in Figure 29. We solved the same elastic problem than in section 'Strategy verification', being the only difference the one related to the material mechanical behavior. The in-plane discretization of was performed by considering a regular mesh composed of 40 × 2 2D bilinear finite elements, whereas 80 one-dimensional linear elements were considered for discretizing the thickness interval I. Thus each ply is represented by 10 elements ensuring a rich description (perhaps too rich in the case of elastic behaviors but that could be compulsory for addressing more complex mechanical behaviors). In this work we would like emphasizing that the proposed strategy allows very rich descriptions of the thickness coordinate without a significant impact on the numerical efficiency.
The displacement field depicted in Figure 30 remains close to the one obtained in section 'Strategy verification' with some slight differences. On the other hand, the radial component of the strain represented in Figure 31 results, as expected, discontinuous across the interplies.

Analysis of large composite structural parts
In this section we address the elastic analysis of a large composite structural part similar to the ones involved in aircraft fuselages. The whole part, 5 m of diameter, and a detail of it are depicted in Figures 32 and 33 respectively. In both figures can be noticed the presence of stiffeners distributed on the laminate surface, all them reinforced along their longitudinal direction.The structure was subjected to a unit internal pressure and an elastic analysis was carried out by assuming the separated representations described in this paper.
The 3D solution was obtained by using an in-plane-out-of-plane separated representation of the elastic fields. The approximation of functions involving the thickness direction was performed by considering 10 elements per ply. An equivalent finite element mesh of the whole part would have implied 7 · 10 7 elements with 2 · 10 8 nodes, each one involving 3 degrees of freedom.  Figure 34 depicts the magnitude of the normal displacement (mm) on the deformed configuration. Even if we considered a laminate consisting of 12 plies, we could consider hundreds of plies without degrading the computing time thanks to the in-plane-out-ofplane separated representation. In fact the mesh considered in the thickness direction only affects the one-dimensional problems defined in that direction that are decoupled of the ones involving the in-plane coordinates. The ones involving the in-plane coordinates being 2D determine the computational complexity of the whole solution procedure.
The main advantage of the procedure that we propose does not concern the possibility of addressing laminates involving several plies, because when elastic behaviors are considered they exist advanced shell elements able to address efficiently such configurations. However, our technique allows considering in an unified description the laminates and the stiffeners placed on it, by using the same in-plane-out-of-plane representation, as was widely discussed in [24]. The consideration of all these stiffeners within a standard shell theory requires special treatments, however in our fully 3D approach, both the shell laminate and the stiffeners are integrated in an unified hypotheses-free description.

Conclusions
We proposed in our former works a procedure based on the separated representation of the displacement field involved in elastic problems defined in plate domains that allowed calculating fully 3D elastic solutions by solving a sequence of 2D and 1D problems, the former ones define in the plane and the last ones in the thickness. Thus, we calculated 3D solutions with a computational cost characteristic of 2D problems, as the ones related to standard plate theories, however here, as we compute directly the 3D solution, we do not need introducing any hypothesis or correction.
In the present paper we extended that methodology for solving elastic linear problems defined in shell geometries. For that purpose we proposed to use a standard mapping of the real geometry to a plate-type parametric domain in which the procedure proposed in our former works was applied. However, the main difficulty when using this procedure concerns the necessity of performing a separated representation of displacement and strain fields, loads, material behavior (elasticity tensor), the differential operators as well as the surface and volume elements. This decomposition is for some problems quite direct and in all cases it can be performed by invoking the singular value decomposition.
The analyzed examples prove the potentiality and efficiency of the proposed strategy, where the computational complexity was found evolving as reported in [24], proving that 3D solutions can be computed at a 2D cost. On the other hand the fact of solving fully 3D models avoids the locking issue that was never found in our analysis.
The extension of this strategy for addressing more complex scenarios involving material or geometrical non-linearities or for considering more complex geometries in which the domain thickness could vary from one point to other represent some of the works in progress.

In-plane and out-of-plane problems
The problem related to the solution of P n is given by  with By integrating Eq. (37) along the coordinate x 3 we obtain a discrete system from which results P n .
By considering now and by integrating in it results the linear problem from with it results T n .

Elements of differential geometry
Consider a shell domain whose middle surface is parametrized by the coordinates (ξ , η), that is X(ξ , η) as depicted in Figure 35. From this mapping we can define the covariant basis a 1 and a 2 , both vectors being tangent to the middle surface: The unit normal vector to the middle surface can be defined from the cross product of a 1 and a 2 according to: Now, we define tensor F containing the Cartesian components of the covariant basis: We can also define the contravariant basis by considering: Figures 36 and 37 show the covariant and the contravariant basis respectively on a torus.

Surface fundamental forms
The length ds defines the first fundamental form I from where a is the metric tensor defined from a = a 1 · a 1 a 1 · a 2 a 2 · a 1 a 2 · a 2 (45) We denote by a the metric tensor determinant, i.e.: The surface element dA writes: dA = a 1 dξ × a 2 dη = a 1 × a 2 · dξ · dη · n (47) and using the cross product property a 1 × a 2 2 = a 1 2 · a 2 2 − (a 1 · a 2 ) 2 it results: In the case of a developable surface, the metric is constant. Figures 38 and 39 depicts a on two surfaces, a developable one in which the metric tensor is constant and another having double curvature. The second fundamental form II is defined by the projection on the normal direction of the second derivative of the mapping: Being n · dX ≡ 0, it results: 0 = d(n · dX) = dn · dX + n · d 2 X = dn · dX + II we can define the mean curvature H: and the Gaussian curvature K :