# On the fully 3D simulations of thermoelastic models defined in plate and shell geometries

Brice Bognet, Adrien Leygue, Francisco Chinesta

► To cite this version:

Brice Bognet, Adrien Leygue, Francisco Chinesta. On the fully 3D simulations of thermoelastic models defined in plate and shell geometries. Revue Européenne de Mécanique Numérique/European Journal of Computational Mechanics, 2012, 21 (1-2), pp.40-51. 10.1080/17797179.2012.702429. hal-04325809

HAL Id: hal-04325809

https://hal.science/hal-04325809v1

Submitted on 11 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# On the fully 3D simulations of thermoelastic models defined in plate and shell geometries

Brice Bognet, Adrien Leygue and Francisco Chinesta*

EADS Foundation Chair, GeM Ecole Centrale of Nantes, 1 rue de la Noe, BP 92101, F 44321
Nantes cedex 3, France

Many models in polymer processing and composites manufacturing are defined in degenerated three dimensional domains (3D), involving plate or shell geometries. The reduction of models from 3D to two dimensional (2D) is not obvious when complex physics or particular geometries are involved. The hypotheses to be introduced for reaching this dimensionality reduction are unclear, and most of the possible proposals will have a narrow interval of validity. The only gateway is to explore new discretisation strategies able to circumvent or at least alleviate the drawbacks related to mesh based discretisations of fully 3D models defined in plate or shell domains. Appropriate separated representation of the involved fields within the context of the proper generalised decomposition allows solving the fully 3D model by keeping a 2D characteristic computational complexity.

De nombreux modèles rencontrés dans la modélisation des procédés de fabrication de structures composites sont définis dans des domaines avec au moins une dimension très petite par rapport aux autres, comme le sont les géométries plaque ou coque. Dans le cas des physiques ou géométries complexes, les hypothèses capables de réduire la dimensionnalité et qui ont fait leurs preuves dans la théorie de plaques et coques restent incertaines. La seule possibilité est de procéder à une résolution 3D. Une représentation séparée des champs concernés dans le cadre de la Proper Generalized Decomposition permet d'effectuer une résolution 3D avec un coût caractéristique des modélisations 2D.

**Keywords:** proper generalised decomposition;; model reduction;; plate geometries

**Mots-clés:** proper generalized decomposition; réduction de modèles; géométries plaques

## 1. Introduction

Many models in polymer processing and composites manufacturing are defined in degenerated three-dimensional (3D) domains. By degenerated we understand that at least one of the characteristic dimensions of the domain is much lower than the other ones. This situation is particularly common in models defined in plate or shell type geometries.

When computing elastic response of plates, two-dimensional (2D) plate theories are usually preferred to the numerically expensive solution of the full 3D elastic problem. Going from a 3D elastic problem to a 2D plate theory model usually involves some kinematical and mechanical hypotheses on the evolution of the solution through the thickness of the plate.

<footnote>*Corresponding author. Email: Francisco.Chinesta@ec.nantes.fr</footnote>
Despite the quality of existing plate theories, their solution close to the plate edges is usu-
ally wrong as the displacement fields are truly 3D in those regions and do not satisfy the
kinematic hypothesis. Moreover, kinematic hypothesis fails where Saint-Venant’s principle
does not apply. It is well known that some heterogeneous complex plates do not verify the
Saint-Venant’s principle nowhere. In that case, the solution of the 3D model is mandatory
even if its computational complexity could be out of the nowadays calculation capabilities.

Most commercial codes for structural mechanics calculations propose a different type of
plate and shell finite elements, even in the case of multilayered composites plates or shells.
However, in composite manufacturing processes the physics encountered in such a multilay-
ered plate or shell domains is much richer, because it usually involves chemical reactions,
crystallisation and strongly coupled and non-linear thermomechanical behaviours. The com-
plexity of the involved physics makes the introduction of pertinent hypotheses for reducing a
priori the dimensionality of the model from 3D to 2D impossible. In that case a fully 3D
modelling is compulsory, and because the richness of the thickness description (many coupled
physics and many plies with different physical states and directions of anisotropy) the approx-
imation of the fields involved in the models needs thousands of nodes distributed along the
thickness direction. Thus, fully 3D descriptions may involve millions of degrees of freedom
that should be solved many times because of the history dependent thermomechanical behav-
iour. Moreover, when we are considering optimisation or inverse identification, many direct
problems have to be solved in order to reach the minimum of a certain cost function.

There have been many attempts at introducing 3D behaviour in models defined in plate and shells not only in the case of homogeneous and isotropic behaviours but also in the case of composite laminates. An exhaustive review is out of the scope of this work, but the interested reader can refer to different families of strategies, ranging in the case of elastic simple thermo-mechanical behaviours from the pure 3D solid models, the 3D type models, the quasi-3D models or the 3D-shell elements (Khalili, Soroush, Davar, & Rahmani, 2011; Matsunaga, 2004; Parente, Fontes Valente, Natal Jorge, Cardoso, & Alves de Sousa, 2006; Qatu, Sullivan, & Wang, 2010; Xue, Cheng, & Hu, 2003; Ye, 2003; Zhang & Yang, 2009). In the context of plate and shell finite elements, there are 1000s of works devoted to the topic.

Today, the solution of such fully 3D models in large structures remains intractable despite
the impressive progresses reached in mechanical modelling, numerical analysis, discretisation
techniques and computer science during the last decade. New numerical techniques are
needed for approaching such complex scenarios, to proceed to the solution of fully 3D multi-
physics models in geometrically complex parts (e.g. a whole aircraft). The well-established
mesh-based discretisation techniques fail because the excessive number of degrees of freedom
involved in the fully 3D discretisations, where very fine meshes are required in the thickness
direction (despite its reduced dimension) and also in the in-plane directions to avoid too dis-
torted meshes or to capture rich in-plane descriptions.

In this manuscript, we propose the application of the model reduction method known as
proper generalised decomposition – PGD – to the simulation of 3D thermomechanical models
defined in plate geometries. This technique was proposed in two recent papers (Ammar, Mok-
dad, Chinesta, & Keunings, 2006, 2007) for circumventing, or at least alleviating, the curse
of dimensionality, then applied in panoply of scenarios (Chinesta, Ammar, & Cueto, 2010;
Chinesta, Ammar, Falco, & Laso, 2007; Chinesta, Ammar, & Joyot, 2008; Mokdad, Ammar,
Normandin, Chinesta, & Clermont, 2010; Mokdad, Pruliere, Ammar, & Chinesta, 2007; Pru-
liere, Ammar, El Kissi, & Chinesta, 2009).

The fully 3D solution of models defined in degenerated domains is also an appealing field
of application of the PGD as proved in Bognet et al. (2012). Consider the unknown field
$u(x, y, z)$ defined in a plate domain $\Xi$. We proposed in Bognet et al. (2012) an in-plane-out-of-plane separated representation

$$u(x, y, z) \approx \sum_{i=1}^{N} X_i(\mathbf{x}) \cdot Z_i(z) \quad (1)$$

with $\mathbf{x} = (x, y) \in \Omega \subset \mathbb{R}^2$ and $z \in I \subset \mathbb{R}$. This strategy is particularly suitable when $\Xi = \Omega \times I$.

Because such decomposition involves the calculation of 2D functions $X_i(\mathbf{x})$ and one-dimensional (1D) functions $Z_i(z)$ (these ones with a computational complexity negligible with respect to the computation of the 2D functions), we can conclude that the computational complexity of the fully 3D solution is of the same order of magnitude than the solution of 2D models.

The main contributions of this work are: (i) the analysis of the first PGD mode and its connection with the usual plate hypotheses; (ii) the extension to profiles in which the separated representation operates by separating the transverse section and the profile axis and (iii) address the extension towards a fully 3D solution of models defined in shell domains. Before addressing the three questions just referred we recall the in-plane-out-of-plane separated representation that we are considering in the displacement representation.

## 2. In-plane-out-of-plane separated representation

When we consider the elastic behaviour of a plate, it suffices considering a separated representation of each component of the displacement vector:

$$\mathbf{u}(x, y, z) = \begin{pmatrix} u(x, y, z) \\ v(x, y, z) \\ w(x, y, z) \end{pmatrix} \approx \sum_{i=1}^{N} \begin{pmatrix} u_{xy}^{i}(x, y) \cdot u_{z}^{i}(z) \\ v_{xy}^{i}(x, y) \cdot v_{z}^{i}(z) \\ w_{xy}^{i}(x, y) \cdot w_{z}^{i}(z) \end{pmatrix} \quad (2)$$

Details on the constructor of such separated representation were given in Bognet et al. (2012).

In order to highlight the interest of such decomposition, we are comparing the complexity of PGD-based solvers with respect to the standard finite element method. For the sake of simplicity, we will consider a hexahedral domain discretised using a regular structured grid with respectively, $N_x$, $N_y$ and $N_z$ nodes in the $x$, $y$ and $z$ directions, respectively. Even if the domain thickness is much lower than the other characteristic in-plane dimensions, the physics in the thickness direction could be quite rich, mainly when we consider composite plates composed of hundreds of anisotropic plies in which the complex physics involved requires fully 3D descriptions. In that case, 1000s of nodes in the thickness direction could be required to represent accurately the solution behaviour in that direction. In usual mesh-based discretisation strategies, this fact induces a challenging issue because the number of nodes involved in the model scales with $N_x \times N_y \times N_z$, however, if one applies a PGD-based discretisation, and the separated representation of the solution involves $N$ modes (terms in the finite sum decomposition), one should solve $N$ 2D problems related to the functions involving the in-plane coordinates and $N$ 1D problems related to the functions involving the thickness coordinate. The computing time related to the solution of the one-dimensional (1D) problems can be neglected with respect to the one required for solving the 2D ones. Thus, the PGD complexity scales as $N \times N_x \times N_y$, $N$ being the number of terms in the decomposition and
$N_x \times N_y$ being the number of nodes for describing each function defined in the $\Omega$. The amount of information in the PGD solution is $N \times (N_x \times N_y + N_z)$, taking into account both the representation of 2D functions defined in $\Omega$ and 1D functions defined in $I$, with $\Xi = \Omega \times I$.

By comparing both complexities, $N_x \times N_y \times N_z$ and $N \times N_x \times N_y$, we can notice that as soon as $N_z \gg N$ the use of PGD-based discretisation leads to impressive computing time savings, making possible even the solution of models never until now solved, even using low performance computing platforms. In our numerical experiments we realise that $N$ is in general of the order of few 10s.

### 3. Analysis of the separated representation of an elastic plate

To validate the proposed separated representation, we consider the square homogeneous plate depicted in Figure 1 and we compare the classical 3D linear elastic finite element solution and the one obtained by using the PGD with an equivalent discretisation, that is, the 2D functions involving the in-plane coordinates in the PGD are approximated using the same mesh that the finite element considered on the plate surface and the 1D functions involving the thickness coordinate when using the PGD were approximated by using the same number of nodes that was considered in the thickness finite element approximation.

The applied load consists of a uniform pressure applied on the upper face. The finite element solution was performed by considering a uniform mesh composed of 100 × 100 × 50 eight-nodes hexahedral elements. The PGD solution was performed by using the uniform mesh composed of 100 × 100 four-nodes elements for approximating the functions involving the in-plane coordinates, whereas a uniform 1D mesh composed of 50 two-nodes 1D linear elements were used for approximating the functions involving the thickness coordinate. Nine modes were needed for approximating the solution when using the PGD, most of them describe the 3D effects that appear in the neighbourhood of the boundaries where the displacement was prescribed. In order to compute these nine terms involved in the separated representation 165 2D and other 165 1D problems were solved. The error (considering the $L^2$ norm) related to the solution performed by using the PGD method, considering as reference solution the one computed by using the FEM, is lower than 0.3% everywhere except in the vicinity of the plate corners where it reaches a value of 0.57%. This error can be reduced by considering more terms in the separated representation, i.e. higher $N$ in the finite sum decomposition.

<figure><img src="image_2.png" /><figcaption>Figure 1. Problem geometry.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2. Comparison of the PGD and FEM based 3D discretisations.</figcaption></figure>

Figure 2 compares the CPU time of both the PGD and the FEM based discretisations for solving the linear elasticity problem previously described as a function of the number of in-plane degrees of freedom, $N_x \times N_y$ and of the number of degrees of freedom in the thickness, $N_z$. We can notice the linear evolution of the computational complexity with the number of in-plane or out-of-plane degrees of freedom when using the PGD instead of the exponential growing when using the finite element discretisation.

In this simple problem, the edge effects are already present and confirm the necessity of several modes to correctly describe the solution in the boundary neighbourhood. In Figures 3

<figure><img src="image_4.png" /><figcaption>Figure 3. First mode of the PGD solution. Left: $v_{xy}^1(x, y)$ top and $v_z^1(z)$ down. Right: $w_{xy}^1(x, y)$ top and $w_z^1(z)$ down.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4. Second mode of the PGD solution. Left: $v_{xy}^2(x, y)$ top and $v_z^2(z)$ down. Right: $w_{xy}^2(x, y)$ top and $w_z^2(z)$ down.</figcaption></figure>

and we depict respectively, the first and the second mode of the PGD solution. The first mode seems to represent classical plate theory solutions because the first mode of the displacements $u_z^1(z)$, $v_z^1(z)$ and $w_z^1(z)$ shows a linear evolution in the thickness direction. The second mode of the PGD solution shows a more complex z-dependence but it should be noticed that in the xy-plane it essentially contributes to the solution in the plate edges neighbourhood where one expects to observe a truly 3D displacement field. The subsequent modes of the PGD solution gradually improve the solution quality close to the plate edges and corners.

<figure><img src="image_6.png" /><figcaption>Figure 5. Evolution of the angle between the cross section and the mid surface with the plate thickness.</figcaption></figure>
As just pointed out, the first PGD mode seems to agree with classical plate theory where the z-displacement is constant in the thickness whereas the in-plane displacements evolve linearly in the thickness. It is well known that depending on the thickness Kirchhoff-Love or Reissner-Mindlin theories apply. Thus for thick enough plates, planar cross-sections orthogonal to the mid-surface, remain after deformation planar but no more orthogonal to the mid-surface. However, when the plate thickness reduces, cross-sections after deformation remain orthogonal to the mid-surface. Figure 5 depicts the evolution of the angle between the deformed cross section and the mid-surface with the plate thickness when considering the displacement field related to the first mode of the separated decomposition calculated by applying the PGD.

These results open numerous perspectives and the analysis of the successive modes could offer a way to define enriched plate theories. This analysis constitutes a work in progress.

## 4. Alternative separated representations

The efficient solution process by applying the PGD requires a separated representation of differential operators and material parameters as illustrated in Bognet et al. (2012). Thus, the in-plane-out-of-plane separated representation of the displacement field requires the following representation of the elasticity tensor:

$$ \mathbf{K}_{ijkl}(x, y, z) \approx \sum_{l=1}^{l} \sum_{l=1}^{M} k_{l}^{ijkl}(\mathbf{x}) \cdot \kappa_{l}^{ijkl}(z) \quad (3) $$

In Bognet et al. (2012), we illustrate the separated representation of complex geometries. Thus, if we consider the structure depicted in Figure 6 involving a homogeneous and isotropic material, the Young's modulus could be written from:

$$ E(x, y, z) \approx E \cdot \sum_{l=1}^{l} \sum_{l=1}^{2} e_{l}(\mathbf{x}) \cdot \varepsilon_{l}(z) \quad (4) $$

where the functions involved in the separated representation (4) are depicted in Figure 7.

More complex geometries can be addressed, however, when thickness evolves a change of coordinates should be applied, even if there are other possibilities that are being explored. Figure 8 illustrates the displacement fields computed by using an in-plane-out-of-plane separated representation by considering an elastic constitutive equation, of two structures subjected to a uniform pressure in the bottom surface with two of its lateral sides clamped and the other two sides free.

In what follows we are considering a more complex geometry, illustrated in Figure 9, that cannot be easily written in an in-plane-out-of-plane separated form.

<figure><img src="image_7.png" /><figcaption>Figure 6. Separable geometry.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7. Functions involved in the separated representation (4).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8. Thermo elastic solution using an in plane out of plane decomposition.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9. Angle geometry.</figcaption></figure>

This geometry constitutes a profile that can be viewed as the cross-section extruded along the axis $z$. Now, the natural separated representation writes again:

$$ \mathbf{u}(x, y, z) = \begin{pmatrix} u(x, y, z) \\ v(x, y, z) \\ w(x, y, z) \end{pmatrix} \approx \sum_{i=1}^{N} \begin{pmatrix} u_{xy}^{i}(x, y) \cdot u_{z}^{i}(z) \\ v_{xy}^{i}(x, y) \cdot v_{z}^{i}(z) \\ w_{xy}^{i}(x, y) \cdot w_{z}^{i}(z) \end{pmatrix} \quad (5) $$

but now $\mathbf{X} = (x, y)$ defines the cross-section coordinates and $z$ the one defined in the axis direction. Functions defined in the cross-section are described from a usual finite element representation as depicted in Figure 10.

Figure 11 depicts the stress field on the deformed part resulting after a cooling of 180 degrees of an angle-laminate initially right composed of eight unidirectional plies, where the fibres were aligned in the sequence [0,90,90,0,0,90,90,0]. The mesh of the cross-section depicted in Figure 11, contains five elements per ply. The resulting 2D mesh contains 2680 bi-linear elements whereas 100 1D linear elements are considered for approximating functions
<figure><img src="image_11.png" /><figcaption>Figure 10. Mesh defined on the cross section.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11. Stress field on the deformed geometry.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12. Coordinate system considered for defining the separated representation in shell domains.</figcaption></figure>

defined in the axis direction. The number of degrees of freedom involved in an equivalent 3D finite element description results 836,400. The use of the PGD techniques allows the solution of such rich model in 77 s using a standard laptop and Matlab.
In those conditions, when we consider an initial angle $\alpha = 90^\circ$, the thermo-elastic deformation induces a final angle slighted lower $\alpha = 89.6^\circ$. Figure 11 depicts one component of the stress tensor on the deformed geometry. If we desire a deformed part characterised by $\alpha = 90^\circ$, we should consider a slightly larger angle. In order to determine such an optimal initial angle we consider an optimisation strategy where the desired parameter consists of the angle $\alpha$, and the cost function is the deviation of the final angle with respect to the desired one. A simple secant strategy converges in five iterations, being the optimal initial angle $\alpha = 90.4^\circ$. The whole optimisation procedure was performed in a few minutes.

## 5. In-plane-out-of-plane extension to shell geometries

When models are defined in shell geometries the decomposition applies by considering again the middle surface and its normal direction. For this purpose, we consider the natural system of coordinates depicted in Figure 12.

In this case, the separated representation writes:

$$ \mathbf{u}(\xi_1, \xi_2, \xi_3) = \begin{pmatrix} u(\xi_1, \xi_2, \xi_3) \\ v(\xi_1, \xi_2, \xi_3) \\ w(\xi_1, \xi_2, \xi_3) \end{pmatrix} \approx \sum_{i=1}^{N} \begin{pmatrix} u_{xy}^{i}(\xi_1, \xi_2) \cdot u_{z}^{i}(\xi_3) \\ v_{xy}^{i}(\xi_1, \xi_2) \cdot v_{z}^{i}(\xi_3) \\ w_{xy}^{i}(\xi_1, \xi_2) \cdot w_{z}^{i}(\xi_3) \end{pmatrix} \quad (6) $$

Then, after some technical manipulations similar to the ones performed in Bognet et al. (2012), in the case of plate geometries, we can construct such separated representation. Figure 13 depicts the elastic 3D solution (displacement normal to the middle surface) defined in a section of aircraft fuselage subjected to an internal pressure.

<figure><img src="image_14.png" /><figcaption>Figure 13. Section of aircraft fuselage subjected to an internal pressure: normal displacement.</figcaption></figure>
## 6. Conclusions

In this paper, we have demonstrated that PGD-based 3D modelling using an in-plane-out-of-plane separated representation is an appealing choice to perform complex 3D simulations of thermo-mechanical models defined in degenerated domains. This approach is particularly suitable for a very wide class of classical problems such as thick plates, laminate geometries or any other complex geometry because of the high impact of edge effects, which cannot be described using classical simplified theories. Even if enriched plates theories exists, enrichments have to be introduced ‘a priori’ to capture edge effects. Using fully 3D PGD separated representation, edge effects or any other 3D local or non-local effect are automatically captured, because of its fully 3D character. Finally, this new in-plane-out-of-plane approach allows solving different 3D modes while keeping a very rich description of the unknown evolution throughout the thickness which cannot be solved by using standard mesh-based discretisations because of their prohibitive computational cost.

## References

* Ammar, A., Mokdad, B., Chinesta, F., & Keunings, R. (2006). A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. *Journal of Non Newtonian Fluid Mechanics, 139*, 153–176.

* Ammar, A., Mokdad, B., Chinesta, F., & Keunings, R. (2007). A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. Part II: transient simulation using space time separated representations. *Journal of Non Newtonian Fluid Mechanics, 144*, 98–121.

* Bognet, B., Leygue, A., Chinesta, F., Poitou, A., & Bordeu, F. (2012). Advanced simulation of models defined in plate geometries: 3D solutions with 2D computational complexity. *Computer Methods in Applied Mechanics and Engineering, 201*, 1–12.

* Chinesta, F., Ammar, A., Falco, A., & Laso, M. (2007). On the reduction of stochastic kinetic theory models of complex fluids. *Modeling and Simulation in Materials Science and Engineering, 15*, 639–652.

* Chinesta, F., Ammar, A., & Joyot, P. (2008). The nanometric and micrometric scales of the structure and mechanics of materials revisited: An introduction to the challenges of fully deterministic numerical descriptions. *International Journal for Multiscale Computational Engineering, 6*(3), 191–213.

* Chinesta, F., Ammar, A., & Cueto, E. (2010). On the use of proper generalized decompositions for solving the multidimensional chemical master Equation. *European Journal of Computational Mechanics, 19*, 53–64.

* Khalili, S.M.R., Soroush, M., Davar, A., & Rahmani, O. (2011). Finite element modeling of low velocity impact on laminated composite plates and cylindrical shells. *Composite Structures, 93*, 1363–1375.

* Matsunaga, H. (2004). A comparison between 2D single layer and 3D layerwise theories for computing interlaminar stresses of laminated composite and sandwich plates subjected to thermal loadings. *Composite Structures, 271*, 651–670.

* Mokdad B., Pruliere, E. Ammar, A., & Chinesta, F. (2007). On the simulation of kinetic theory models of complex fluids using the Fokker Planck approach. *Applied Rheology, 17*(2), 26494, 1–14.

* Mokdad, B., Ammar, A., Normandin, M., Chinesta, F., & Clermont, J.R. (2010). A fully deterministic micro macro simulation of complex flows involving reversible network fluid models. *Mathematics and Computer in Simulation, 80*, 1936–1961.

* Parente, M.P.L., Fontes Valente, R.A., Natal Jorge, R.M., Cardoso, R.P.R., & Alves de Sousa, R.J. (2006). Sheet metal forming simulation using EAS solid shell finite elements. *Finite Elements in Analysis and Design, 42*, 1137–1149.

* Pruliere, E., Ammar, A., El Kissi, N., & Chinesta, F. (2009). Recirculating flows involving short fiber suspensions: Numerical difficulties and efficient advanced micro macro solvers. *Archives of Computational Methods in Engineering, State of the Art Reviews, 16*, 1–30.

* Qatu, M.S., Sullivan, R.W., & Wang, W. (2010). Recent research advances on the dynamic analysis of composite shells: 2000–2009. *Composite Structures, 93*, 14–31.
Xue, M., Cheng, L., & Hu, N. (2003). The stress analysis of sandwich shells faced with composite sheets based on 3D FEM. *Composite Structures*, **60**, 33–41.

Ye, J. (2003). *Laminated composite plates and shells: 3D modeling*. London: Springer Verlag.

Zhang, Y.X., & Yang, C.H. (2009). Recent developments in finite element analysis for laminated composite plates. *Composite Structures*, **88**, 147–157.