# Smoothed finite element method implemented in a resultant eight-node solid-shell element for geometrical linear analysis

Xavier J-G Elie-Dit-Cosaque, Augustin Gakwaya, Hakim Naceur

► To cite this version:

Xavier J-G Elie-Dit-Cosaque, Augustin Gakwaya, Hakim Naceur. Smoothed finite element method implemented in a resultant eight-node solid-shell element for geometrical linear analysis. Computational Mechanics, 2015, 55 (1), pp.105-126. 10.1007/s00466-014-1085-2. hal-01348306

HAL Id: hal-01348306

https://hal.science/hal-01348306

Submitted on 22 Jul 2016

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
# Smoothed finite element method implemented in a resultant eight-node solid-shell element for geometrical linear analysis

Xavier J-G Élie-Dit-Cosaque · Augustin Gakwaya · Hakim Naceur

**Abstract** A smoothed finite element method formulation for the resultant eight-node solid-shell element is presented in this paper for geometrical linear analysis. The smoothing process is successfully performed on the element mid-surface to deal with the membrane and bending effects of the stiffness matrix. The strain smoothing process allows replacing the Cartesian derivatives of shape functions by the product of shape functions with normal vectors to the element mid-surface boundaries. The present formulation remains competitive when compared to the classical finite element formulations since no inverse of the Jacobian matrix is calculated. The three dimensional resultant shell theory allows the element kinematics to be defined only with the displacement degrees of freedom. The assumed natural strain method is used not only to eliminate the transverse shear locking problem encountered in thin-walled structures, but also to reduce trapezoidal effects. The efficiency of the present element is presented and compared with that of standard solid-shell elements through various benchmark problems including some with highly distorted meshes.

**Keywords** Resultant solid-shell element · Smoothed finite element method (SFEM) · Polygonal element · Strain smoothing · Mesh sensitivity · Accuracy

## 1 Introduction

In the finite element method (FEM), a lot of work dealing with thin shell structures has been accomplished, led by industry needs. This category of three dimensional (3D) shell structures is useful in many industrial activities, like those involved in the aerospace industry. It is characterized by the existence of a small dimension in one direction, identified as the shell thickness, when compared to the two remaining directions.

The many existing shell element formulations can be grouped into three main classes: (i) classical shell elements depicted into plate elements and shell elements. Plate elements combine membrane and bending characteristics, thus including in-plane and out-of-plane theories. For each node, they use both translational and rotational degrees of freedom. They are widely used in practice because of their simple formulation and attractive time efficiency. Shell elements have the same characteristics as plate elements, except for the fact that they also remain accurate in arbitrary shape contexts. The reader can consult Betsch et al. [1], Bischoff and Ramm [2] and Cardoso and Yoon [3]; (ii) degenerated shell elements based on a 3D variational formulation, as for solid elements, have kinematics expressed on the shell mid-surface exhibiting translational and rotational degrees of freedom. For more details on this topic, one can consult, among others, the following authors: [4–13] and [14]; (iii) solid-shell elements, which are among the most recent shell element developments, are 3D shell elements.

Some advantages of solid-shell compared with shell elements formulations are: (i) the capacity of modeling 3D geometries comprising both thin and thick portions without any need for special transition elements; (ii) boundaries modeling without any extra kinematics assumptions. Since physical nodes are located on the bottom and top surfaces of
the shell structure, the contact definition with associated friction phenomenon can be easily dealt with; (iii) the kinematics description remains simple since their formulation uses only displacement degree of freedoms, avoiding the complex update of the rotation vector in non-linear problems. See [15,16] and [17] for details.

Notwithstanding the above classification, shell elements can be separated into two categories: (i) thin shell elements described by Kirchhoff-Love plate theory, which neglects transverse shear deformation, are described in Timoshenko and Woinowsky-Krieger [18]. The elements using this theory require a $C^1$ displacement continuity; (ii) moderate thick shell elements described by Reissner-Mindlin plate theory and used in Bathe and Dvorkin [10], a theory which considers transverse shear, (see Timoshenko and Woinowsky-Krieger [18]). However, solid-shell elements using only displacement degrees of freedom and based in particular on the Reissner-Mindlin theory, use shape functions that are not able to describe the bending induced kinematics in the thickness direction, which may be responsible for locking occurrences. The locking problems are mainly transverse shear locking (cited above), volumetric locking, Poisson thickness locking, curvature or trapezoidal locking and membrane locking. The transverse shear locking may occur when a shell element with a thickness tending towards to zero is subject to bending.

Among the most popular methods to solve transverse shear locking is found the assumed natural strain (ANS) method, used for instance in Dvorkin and Bathe [8] to remedy to transverse shear locking in shell elements. Then, Klinkel et al. [19], Kim et al. [20] and more recently Hannachi [21] introduced it in solid-shell elements. According to ANS method, the natural transverse shear strains are calculated in four sampling points located at the center of the element mid-surface edges. The transverse shear strain is then calculated at element integration points by interpolating the strains evaluated at the element mid-surface. However, to provide a good accuracy, the four sampling points used above are sufficient only in fully integrated elements. An extended ANS method using eight sampling points initially introduced in Cardoso et al. [22] and also used in Schwarze and Reese [23] and Nguyen [24] is however preferred for reduced integration elements.

Another popular technique, used to eliminate transverse shear locking is the enhanced assumed strain (EAS) method, which consists in adding in the deformation field, another field of internal variables that creates additional modes of deformation. The EAS approach is employed in the work of Andelfinger and Ramm [25], Alves de Sousa et al. [26] and in the work of Wriggers et al. [27], Alves de Sousa et al. [28], Li et al. [29] and Fontes Valente et al. [30] for small and large deformations, respectively. The volumetric locking can occur when a given material remains nearly incompressible or incompressible. Simo and Rifai [12] have demonstrated

that applying the EAS method is a way to reduce the volumetric locking issue. The Poisson thickness locking may take place with the hypothesis of a linear displacement variation and, thus, when considering a constant strain along the thickness. This constant strain in the thickness direction may induce inconsistencies since the thickness strain varies linearly as a combination between Poisson's effect and in-plane strain (Büchter et al. [31]). More recently Schwarze and Reese [23] have successfully used this EAS technique in an eight-node solid shell element in a full 3D integration context in order to avoid the Poisson thickness locking. The curvature or trapezoidal locking may appear when an element is used to mesh curved structures and thus when an edge in the thickness direction of a given element is not perpendicular to its mid-surface. Following the same idea as for transverse shear locking treatment, Betsch et al. [1] applied the ANS method to tackle curvature locking. Accordingly, the natural thickness strain is calculated in four sampling points located at the corners of the element mid-surface. It is then interpolated at the integration points. The membrane locking could emanate from a configuration whereby membrane strains remain very small relative to the bending strains in thin structure bending problems. In order to overcome membrane locking issues, Miehe [32] developed a solid type element using one enhanced parameter. Mesh distortion may also be responsible for inaccurate or oscillatory results as proven in Lyly et al. [33] and Kouhia [34].

In spite of the prolific and insightful prior research aiming at improving the FEM, some unresolved issues remain with the technique, such as those related to element distortion and numerical integration. In an effort to enhance the accuracy of numerical solutions for irregular meshes, an approximation method named smoothed finite element method (SFEM) was recently proposed by Liu et al. [35], [36]. This method integrates the conventional FEM technology with the strain smoothing technique that was introduced by Chen et al. [37] while devising a stabilized nodal integration scheme in mesh-free methods (element free Galerkin method). Essentially, the SFEM consists in dividing each element into smoothing cells over which a strain smoothing operation is performed. The strain energy in each smoothing cell is expressed as an explicit form of the smoothed strain. A subsequent use of the divergence theorem then allows the integration to be transformed into one over the cell boundary eliminating the requirements to use shape function derivatives. In Liu et al. [35], the SFEM was applied to a two dimensional static continuum element and implemented in four-node quadrilateral and polygonal elements. Compared with standard FEM, the work of Liu et al. [35], [36] indicate that SFEM has interesting features such as: (i) simplicity of the method, since no Cartesian derivatives of shape functions are used; (ii) elimination of the inverse of the Jacobian matrix; (iii) apparent insensitivity of the SFEM to mesh distortion. This thus makes
finite elements based on this technique more attractive in sit-
uations requiring convergence for very complex mesh and /
or adaptive meshing when standard finite elements are used;
(iv) compared with standard finite elements, the proposed
SFEM elements appear to be more accurate and significantly
more attractive in terms of computer time efficiency. Results
accuracy could also be increased by using more smoothing
cells.

The variational weak form used in SFEM can be derived
from the Hu-Washizu three-fields variational principle or by
the assumed strain method proposed by Simo and Hughes
[38]. Liu et al. [36] then detailed the SFEM theory and found
that it remains equivalent to standard FEM theory for one
dimension as well as for two dimensional (2D) triangular
elements, contrasting with quadrilateral elements. Moreover,
they investigated on the variational consistency of this last
type of element regarding the number of smoothing cells
used. They observed that if the formulation remains varia-
tionally consistent with one cell, this observation does not
hold in cases using multiple cells. However, when using the
same set of shape functions for a given element, the for-
mulation remains energy consistent since the element nodal
displacements remain continuous. Nguyen-Xuan et al. [39]
and Bordas and Natarajan [40], in their theoretical study of
convergence, accuracy and properties of the SFEM tech-
nique, also contributed to a better understanding of SFEM
and gave a good overview of the SFEM driving concept
as well as some interesting results obtained with a SFEM
quadrilateral shell element. The SFEM applied previously
to static problems was extended to 2D dynamic problems in
Dai and Liu [41], [42]. It is demonstrated in these papers that
the approach yields better results than standard quadrilateral
finite elements. The SFEM was initially used with conven-
tional isoparametric elements in Liu et al. [35]. An extension
of the method to polygonal elements is then performed in
Dai et al. [43] and further extended to node-centered (Liu et
al. [44]) and edge-centered (Nguyen-Thoi et al. [45]). It is
observed that SFEM for general n-sided polygonal elements
works well and yields very accurate results, especially in
cases with heavily distorted meshes. However, due to the cor-
responding high computational cost induced by the inherent
higher number of nodes of the technique, they suggest using a
combination of nSFEM and quadrilateral SFEM elements to
respectively mesh the exterior and interior of a part. Nguyen-
Thanh et al. [46] investigated on some methods to solve shear
and incompressible locking problems for quadrilateral SFEM
elements. They successfully presented three selective inte-
gration schemes. Scheme 1 is based on the decomposition of
the material matrix into two parts; one containing the shear
modulus or μ terms and another containing the remaining
terms. Scheme 2 is based on the B-bar approach of Hughes
[47] to overcome volumetric locking. Scheme 3, used to over-
come shear locking, consists in using different smoothing

cell numbers for each stiffness matrix contribution. How-
ever, this last scheme is not that attractive compared with
standard finite elements, since the shear stiffness matrix still
needs Gauss integration points to be determined. Then, var-
ious developments have been made on SFEM quadrilateral
plate and shell elements. For example, Nguyen-Xuan et al.
[48] developed a four nodes plate element using the SFEM
to deal with membrane and bending stiffness matrix compo-
nents. Their proposed element exhibited good accuracy and
time efficiency compared with equivalent standard elements.
Nguyen-Thanh et al. [49] then extended this work to quadri-
lateral flat shell elements and emphasized the high poten-
tial capabilities of the SFEM. Unlike the MITC4 element,
which uses an isoparametric mapping, their proposed ele-
ment calculates the membrane and bending parts of the stiff-
ness matrix directly in the local Cartesian coordinate system.
As a matter of fact, it remains accurate even with a distorted
mesh. Moreover, this new element is found to be insensitive to
mesh distortion when compared to the standard MITC4 ele-
ment. More recently, Wu and Wang [50] improved this work
by revisiting its variational formulation. They demonstrated
that moving the Gauss points to the SFEM cell-based centroid
position has a softening effect on the structure response. A
key point of their work is the smoothing of the assumed shear
strain of the MITC4 element using the edge-based strain
smoothing method presented in Liu et al. [51]. Their pro-
posed element shows improved results accuracy, especially
in cases of a highly distorted mesh. The extension of SFEM to
3D solid continuum elements is proposed by Nguyen-Xuan
[52] and Bordas et al. [53]. In their work, they used a stabiliza-
tion technique to avoid membrane and volumetric locking.
This is consistent with the stabilization technique used in
the meshfree method proposed by Puso and Solberg [54]. In
order to reduce computation time, their 3D stiffness matrix
is calculated using reduced integration on each facet of the
element cells, instead of a standard 2 × 2 Gauss quadrature
formula. However, a stabilization strategy was developed in
order to keep some element efficiency and to eliminate the
zero energy modes. This method consists in considering the
stiffness matrix as a linear combination of the one subcell ele-
ment and the four or eight subcell element. Due to the relative
novelty of the SFEM, its non-linear development and asso-
ciated applications still remain limited. Cui et al. [55] devel-
oped quadrilateral non-linear SFEM plate and shell elements.
They then illustrated their good accuracy as well as their
rapid convergence towards the analytic solutions. Moreover,
Liu et al. [56] successfully compared a quadrilateral shell
element using SFEM with classical finite element through
non-linear examples using plastic materials. More recently,
Élie-Dit-Cosaque et al. [57] developed a non-linear resultant
eight-node solid-shell element. Their proposed element gives
promising results, but some developments are still required
in order to improve its accuracy in non-linear applications.
One should also mention a method called the partition element method (PEM), detailed in Rashid and Sadri [58]. In this method, similarly to SFEM, an element is divided into quadrature cells over which the shape functions are piecewise linear. The averaged strain operators are then calculated in each cell. A major difference between the PEM and the SFEM is that the PEM method allows arbitrary cells to be defined which gives it greater flexibility.

In summary, previous works on the SFEM have observed the following main features of the method: (i) simple implementation; (ii) insensitivity to mesh distortion; (iii) better accuracy as well as higher convergence rates compared with standard equivalent finite element in many cases; (iv) time efficiency.

However, to the best of our knowledge, the following issues are either scarce or did not receive much attention: (i) a full assessment of the theoretical basis of the stabilization techniques proposed in Nguyen-Xuan [52] and Bordas et al. [53]; (ii) the performance of the eight-node solid element using the SFEM in case of a distorted mesh; (iii) application of SFEM shell elements formulation to geometric and material non-linear problems; (iv) implementation of the SFEM approach in solid-shell elements.

In the present paper, the SFEM is applied to an eight node resultant solid-shell element presented in Kim et al. [20] as an extension of the classical resultant shell theory and in Simo et al. [59], in order to evaluate the membrane and bending element stiffness. The strain smoothing method is then applied in a way similar to its application to a quadrilateral shell element (Nguyen-Thanh et al. [49]) without any use of stabilization, since no reduced integration is required. A preliminary simplified version of this work has already been proposed in Élie-Dit-Cosaque et al. [57]. A complete development of the smoothing technique is proposed in the following sections. In order to reduce element locking problems inherent to the proposed solid-shell element extended to SFEM formulation, the well-known ANS method is applied, in a way that is similar to the procedure proposed in Bordas et al. [53]. The transverse normal strain is also approximated using the ANS technique. Special care has been taken to adequately explain how the SFEM is applied to a resultant solid-shell formulation. This paper organization is described in what follows. In Sects. 2, 3 and 4, the geometric description, the kinematics and the finite element approximation of the resultant solid-shell element are presented. The related virtual work principle is then presented in Sect. 5. The SFEM approximation applied to solid-shell element and its parameters are then discussed in Sect. 6. The details of SFEM discretization of the weak form and the computer implementation including detailed element stiffness matrix structures are then given. Various benchmark examples for geometrically linear problems are presented in Sect. 7 to validate the approach. Finally, a conclusion is given in the last section.

<figure><img src="image_2.png" /><figcaption>Fig. 1 Solid-shell element parametrization, undeformed configuration</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2 Solid-shell element configurations</figcaption></figure>

## 2 Geometric description

A solid-shell element based on Reissner-Mindlin kinematics is considered and represented parametrically as in Fig. 1. Only three displacement degrees of freedom are introduced at each node as in a standard solid element (Hauptmann and Schweizerhof [16]).

In order to describe the shell-like behavior of the eight-node solid element, it is useful to introduce a set of three coordinates systems as shown in Fig. 2, inspired from Vu-Quoc and Tan [60] and further detailed in Table 1: a global Cartesian coordinate system $\mathbf{X}(X, Y, Z)$, a convected coordinate system $\xi(\xi^1, \xi^2, \xi^3)$ and a local Cartesian coordinate system $\mathbf{x}(x, y, z)$ whose orientation changes with the configuration. The latter is set up in the element mid-surface in order to be able to handle shell features such as membrane, bending and shear behavior of thin shell structures. The relationships between the different coordinate basis are given in several monographs and papers for instance Domissy [61] and Hannachi [21].
<table><caption>**Table 1** Coordinate systems definition</caption><thead><tr><th>Level</th><th>Base vectors</th><th>Coordinates</th></tr></thead><tbody><tr><td>Global Cartesian</td><td>(e<sub>X</sub>, e<sub>Y</sub>, e<sub>Z</sub>)</td><td>X<sup>0</sup>(X<sup>0</sup>, Y<sup>0</sup>, Z<sup>0</sup>) or X<sub>i</sub><sup>0</sup>, i = 1, 2, 3 undeformed configuration<br>X(X, Y, Z) or X<sub>i</sub>, i = 1, 2, 3 deformed configuration</td></tr><tr><td>Natural</td><td>(G<sub>1</sub>, G<sub>2</sub>, G<sub>3</sub>)<br>(g<sub>1</sub>, g<sub>2</sub>, g<sub>3</sub>)</td><td>ξ(ξ<sup>1</sup>, ξ<sup>2</sup>, ξ<sup>3</sup>) or ξ<sup>i</sup>, i = 1, 2, 3</td></tr><tr><td>Local Cartesian</td><td>(R<sub>x</sub>, R<sub>y</sub>, R<sub>z</sub>)<br>(r<sub>x</sub>, r<sub>y</sub>, r<sub>z</sub>)</td><td>x<sup>0</sup>(x<sup>0</sup>, y<sup>0</sup>, z<sup>0</sup>) or x<sub>i</sub><sup>0</sup>, i = 1, 2, 3 undeformed configuration<br>x(x, y, z) or x<sub>i</sub>, i = 1, 2, 3 deformed configuration</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Resultant eight-node solid-shell element geometry and mid-surface concept</figcaption></figure>

The coordinates of an arbitrary point $q_0$ on the solid-shell structure in the undeformed configuration can then be defined in term of the position vector represented parametrically by:

$$ \mathbf{X}_{q_0}^{0} (\xi^1, \xi^2, \xi^3) = \frac{(1+\xi^3)}{2} \cdot \mathbf{X}^{0,U} (\xi^1, \xi^2) \\ + \frac{(1-\xi^3)}{2} \cdot \mathbf{X}^{0,L} (\xi^1, \xi^2) \\ (-1 \le \xi^1, \xi^2, \xi^3 \le 1) \qquad (1) $$

which, after introducing the mid-surface concepts and considering the parameter $\xi^3$ as the thickness direction (see Fig. 3), can also be written as:

$$ \mathbf{X}_{q_0}^{0} (\xi^1, \xi^2, \xi^3) = \bar{\mathbf{X}}^{0} (\xi^1, \xi^2) + \xi^3 \cdot \Delta\mathbf{X}^{0} (\xi^1, \xi^2), \\ (-1 \le \xi^1, \xi^2, \xi^3 \le 1) \qquad (2) $$

where $\mathbf{X}^{0,L}$ and $\mathbf{X}^{0,U}$ are respectively the position vector of a point on the lower and upper surface in the global Cartesian coordinate system $\mathbf{X}^{0}(X^{0}, Y^{0}, Z^{0})$ of the undeformed configuration, $\bar{\mathbf{X}}^{0} = \frac{\mathbf{X}^{0,U} + \mathbf{X}^{0,L}}{2}$ is the position vector of the mid-surface and $\Delta\mathbf{X}^{0} = \frac{\mathbf{X}^{0,U} - \mathbf{X}^{0,L}}{2}$ is a director vector pointing from the lower to the upper surface.

In the deformed configuration, the coordinates of an arbitrary point $q$ in the global Cartesian coordinate system $\mathbf{X}(X, Y, Z)$ can similarly be represented as:

$$ \mathbf{X}_q (\xi^1, \xi^2, \xi^3) = \bar{\mathbf{X}} (\xi^1, \xi^2) + \xi^3 \cdot \Delta\mathbf{X} (\xi^1, \xi^2), \\ (-1 \le \xi^1, \xi^2, \xi^3 \le 1) \qquad (3) $$

By introducing the displacements field vector $\mathbf{U}_q$, the undeformed and the deformed configurations are related by the standard mapping at point $q$, as:

$$ \mathbf{X}_q (\xi^1, \xi^2, \xi^3) = \mathbf{X}_{q_0}^0 (\xi^1, \xi^2, \xi^3) + \mathbf{U}_q (\xi^1, \xi^2, \xi^3) \qquad (4) $$

in the global Cartesian coordinate system $\mathbf{X}(X, Y, Z)$. Similarly, as in Eqs. (2) and (3), by using mid-surface related parameters, the displacement vector in Eq. (4) can be written as:

$$ \mathbf{U}_q (\xi^1, \xi^2, \xi^3) = \bar{\mathbf{U}} (\xi^1, \xi^2) + \xi^3 \cdot \Delta\mathbf{U} (\xi^1, \xi^2), \\ (-1 \le \xi^1, \xi^2, \xi^3 \le 1) \qquad (5) $$

where $\bar{\mathbf{U}}$ and $\Delta\mathbf{U}$ are defined as $\bar{\mathbf{U}} = \frac{\mathbf{U}^U + \mathbf{U}^L}{2}$ and $\Delta\mathbf{U} = \frac{\mathbf{U}^U - \mathbf{U}^L}{2}$, with $\mathbf{U}^L$ and $\mathbf{U}^U$ being respectively the displacement vector of a point on the lower and upper surface of the shell structure.

The covariant base vectors $(G_1, G_2, G_3)$ respectively the contravariant base vectors $(G^1, G^2, G^3)$ take the following forms in the undeformed configuration:

$$ \mathbf{G}_i = \frac{\partial \mathbf{X}^0}{\partial \xi^i}, \quad i = 1, 2, 3 \qquad (6) $$

$$ \mathbf{G}^i = \frac{\partial \xi^i}{\partial \mathbf{X}^0}, \quad i = 1, 2, 3 \qquad (7) $$

In the same idea, the covariant base vectors $(g_1, g_2, g_3)$, respectively the contravariant base vectors $(g^1, g^2, g^3)$ take the following forms in the deformed configuration:

$$ \mathbf{g}_i = \frac{\partial \mathbf{X}}{\partial \xi^i}, \quad i = 1, 2, 3 \qquad (8) $$

$$ \mathbf{g}^i = \frac{\partial \xi^i}{\partial \bar{\mathbf{X}}}, \quad i = 1, 2, 3 \qquad (9) $$
The local Cartesian base vectors are defined in the unde-
formed configuration from the covariant base vectors as:

$$
\mathbf{R}_z = \frac{\mathbf{G}_1 \wedge \mathbf{G}_2}{\|\mathbf{G}_1 \wedge \mathbf{G}_2\|} \qquad (10)
$$

$$
\mathbf{R}_x = \frac{\mathbf{G}_1}{\|\mathbf{G}_1\|} \qquad (11)
$$

$$
\mathbf{R}_y = \mathbf{R}_z \wedge \mathbf{R}_x \qquad (12)
$$

Similarly, the local Cartesian base vectors are defined in the
deformed configuration from the covariant base vectors as:

$$
\mathbf{r}_z = \frac{\mathbf{g}_1 \wedge \mathbf{g}_2}{\|\mathbf{g}_1 \wedge \mathbf{g}_2\|} \qquad (13)
$$

$$
\mathbf{r}_x = \frac{\mathbf{g}_1}{\|\mathbf{g}_1\|} \qquad (14)
$$

$$
\mathbf{r}_y = \mathbf{r}_z \wedge \mathbf{r}_x \qquad (15)
$$

So, the local Cartesian coordinates **x** (x, y, z) are related
to the global Cartesian coordinates **X** (X, Y, Z) through the
transformation matrix **Q**₁ = [ **r**ₓ **r**ᵧ **r**ᵧ] such that the local
components of a position or displacement vector **u**q =
⟨ **u** **v** **w**⟩ of an arbitrary point **q** are given by:

$$
\mathbf{x}_q = \mathbf{Q}_1^T \cdot \mathbf{X}_q \qquad (16)
$$

$$
\mathbf{u}_q = \mathbf{Q}_1^T \cdot \mathbf{U}_q \qquad (17)
$$

For sake of simplify, the subscript ( )<sub>q</sub> is omitted in the
subsequent developments.

3 Kinematics of shell deformation

The Green–Lagrange strain tensor components $\epsilon_{kl}^{\xi}$, oriented
in the natural coordinate system $\xi$ ($\xi^1, \xi^2, \xi^3$) are (see Vu-
Quoc and Tan [60] for more detail):

$$
\epsilon_{kl}^{\xi} = \frac{1}{2} \left( \mathbf{g}_k \frac{\partial \mathbf{u}}{\partial \xi^l} + \frac{\partial \mathbf{u}}{\xi^k} \mathbf{g}_l + \frac{\partial \mathbf{u}}{\xi^k} \frac{\partial \mathbf{u}}{\xi^l} \right), \quad k, l = 1, 2, 3 \quad (18)
$$

Considering only the linear terms, the above Eq. (18)
becomes:

$$
\epsilon_{kl}^{\xi} = \frac{1}{2} \left( \mathbf{g}_k \frac{\partial \mathbf{u}}{\partial \xi^l} + \frac{\partial \mathbf{u}}{\xi^k} \mathbf{g}_l \right), \quad k, l = 1, 2, 3 \qquad (19)
$$

Considering Eq. (19) component by component, the individ-
ual components of $\epsilon_{kl}^{\xi}$ are then:

$$
\begin{align}
\epsilon_{11}^{\xi} &= \mathbf{g}_1 \cdot \frac{\partial \mathbf{u}}{\partial \xi^1} & \epsilon_{22}^{\xi} &= \mathbf{g}_2 \cdot \frac{\partial \mathbf{u}}{\partial \xi^2} \nonumber \\
\gamma_{12}^{\xi} &= \mathbf{g}_1 \cdot \frac{\partial \mathbf{u}}{\partial \xi^2} + \mathbf{g}_2 \cdot \frac{\partial \mathbf{u}}{\partial \xi^1} & \gamma_{13}^{\xi} &= \mathbf{g}_1 \cdot \frac{\partial \mathbf{u}}{\partial \xi^3} + \mathbf{g}_3 \cdot \frac{\partial \mathbf{u}}{\partial \xi^1} \nonumber \\
\gamma_{23}^{\xi} &= \mathbf{g}_2 \cdot \frac{\partial \mathbf{u}}{\partial \xi^3} + \mathbf{g}_3 \cdot \frac{\partial \mathbf{u}}{\partial \xi^2} & \epsilon_{33}^{\xi} &= \mathbf{g}_3 \cdot \frac{\partial \mathbf{u}}{\partial \xi^3} \tag{20}
\end{align}
$$

The Green–Lagrange strain tensor **ε**, can be formulated in
the directions of the local Cartesian base vectors as:

$$
\boldsymbol{\epsilon} = \epsilon_{ij} \mathbf{r}_i \otimes \mathbf{r}_j, \quad i, j = 1, 2, 3 \qquad (21)
$$

It can also be formulated in the direction of the contravariant
base vectors (g¹, g², g³) defined in Eq. (9) as:

$$
\begin{align*}
\boldsymbol{\epsilon} &= \epsilon_{kl}^{\xi} \mathbf{g}_i^k \otimes \mathbf{g}_j^l \\
&= \epsilon_{kl}^{\xi} \frac{\partial \xi^k}{\partial X_i} \mathbf{r}_i \otimes \frac{\partial \xi^l}{\partial X_j} \mathbf{r}_j \\
&= \epsilon_{kl}^{\xi} \frac{\partial \xi^k}{\partial X_i} \frac{\partial \xi^l}{\partial X_j} \mathbf{r}_i \otimes \mathbf{r}_j, \quad i, j, k, l = 1, 2, 3 \tag{22}
\end{align*}
$$

where $\mathbf{X}_i$, $i = 1, 2, 3$ are the global Cartesian coordinates
defined in Table 1.

Identifying the terms of Eqs. (21) and (22), the transforma-
tion between the local Cartesian and the natural components
of the Green–Lagrange strain tensor represented respectively
by $\epsilon_{ij}$ and $\epsilon_{kl}^{\xi}$ is:

$$
\epsilon_{ij} = \epsilon_{kl}^{\xi} \frac{\partial \xi^k}{\partial X_i} \frac{\partial \xi^l}{\partial X_j}, \quad i, j, k, l = 1, 2, 3 \quad (23)
$$

Starting from the deformation gradient **F** defined in term of
the displacements field described in Eq. (5) as **F** = ∂**x** / ∂**X** =
**I** + ∇**X** **u**, the Green–Lagrange strain tensor **ε** oriented in the
local Cartesian coordinate system can be written as:

$$
\begin{align}
\boldsymbol{\epsilon} &= \frac{1}{2} (\mathbf{F}^T \mathbf{F} - \mathbf{I}) \nonumber \\
&= \frac{1}{2} \left\{ (\nabla_X \mathbf{u}) + (\nabla_X \mathbf{u})^T + (\nabla_X \mathbf{u})^T (\nabla_X \mathbf{u}) \right\} \tag{24}
\end{align}
$$

Limiting the following development to the linear strain com-
ponents, as for the natural expression of the Green–Lagrange
strain tensor, gives the expression:

$$
\epsilon = \frac{1}{2} \left\{ (\nabla_X \mathbf{u}) + (\nabla_X \mathbf{u})^T \right\} \quad (25)
$$

Following the concepts introduced in Eq. (5), the displace-
ments vector can also be written as:

$$
\mathbf{u} = \langle \mathbf{u} \mathbf{v} \mathbf{w} \rangle^T = \left\langle \bar{\mathbf{u}} + \xi^3 \cdot \Delta \mathbf{u} \bar{\mathbf{v}} + \xi^3 \cdot \Delta \mathbf{v} \bar{\mathbf{w}} + \xi^3 \cdot \Delta \mathbf{w} \right\rangle^T \quad (26)
$$

with $\bar{\mathbf{u}} = \begin{bmatrix} \bar{u} \\ \bar{v} \\ \bar{w} \end{bmatrix}$ as the displacement mean value on the mid-
surface and $\Delta \mathbf{u} = \begin{Bmatrix} \Delta u \\ \Delta v \\ \Delta w \end{Bmatrix}$ as the through-thickness relative
displacement vector.

Separating the in-plane and the out-of-plane degrees of
freedoms from Eq. (26), the displacements vector for the
resultant theory can then be defined as:

$$
\tilde{\mathbf{u}} = \langle \bar{u} \quad \bar{v} \quad \bar{w} \quad \Delta u \quad \Delta v \quad \Delta w \rangle^T \quad (27)
$$
Then, using the displacement vector defined as in Eq. (5), the linear terms of the Green–Lagrange strain tensor can be expressed in the local Cartesian coordinate system **x** (x, y, z) defined in Table 1 as:

$$ \varepsilon_{xx} = \frac{\partial u}{\partial x} = \frac{\partial \bar{u}}{\partial x} + \xi^3 \cdot \frac{\partial \Delta u}{\partial x}; $$

$$ \varepsilon_{yy} = \frac{\partial v}{\partial y} = \frac{\partial \bar{v}}{\partial y} + \xi^3 \cdot \frac{\partial \Delta v}{\partial y}; $$

$$ \gamma_{xy} = \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} = \frac{\partial \bar{u}}{\partial y} + \frac{\partial \bar{v}}{\partial x} + \xi^3 \cdot \left( \frac{\partial \Delta u}{\partial y} + \frac{\partial \Delta v}{\partial x} \right) $$

$$ \varepsilon_{zz} = \frac{\partial w}{\partial z} = \frac{\partial \bar{w}}{\partial z} + \frac{\partial \xi^3}{\partial z} \cdot \Delta w + \xi^3 \cdot \frac{\partial \Delta w}{\partial z} $$

$$ \gamma_{yz} = \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y} = \frac{\partial \bar{v}}{\partial z} + \frac{\partial \bar{w}}{\partial y} + \xi^3 \cdot \left( \frac{\partial \Delta v}{\partial z} + \frac{\partial \Delta w}{\partial y} \right) + \frac{\partial \xi^3}{\partial z} \cdot \Delta v $$

$$ \gamma_{xz} = \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} = \frac{\partial \bar{u}}{\partial z} + \frac{\partial \bar{w}}{\partial x} + \xi^3 \cdot \left( \frac{\partial \Delta u}{\partial z} + \frac{\partial \Delta w}{\partial x} \right) + \frac{\partial \xi^3}{\partial z} \cdot \Delta u \quad (28) $$

Developing the terms of Eq. (23) and reorganizing the resulting equations in tensor and matrix forms, the transformation of the Green–Lagrange strain tensor components from the natural to the local coordinate system can be expressed as:

$$ \begin{bmatrix} \varepsilon_{xx} \\ \varepsilon_{yy} \\ \varepsilon_{zz} \\ \gamma_{xy} \\ \gamma_{yz} \\ \gamma_{xz} \end{bmatrix} = \mathbf{T} \cdot \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ \gamma_{12} \\ \gamma_{23} \\ \gamma_{13} \end{bmatrix} \quad (29) $$

where

$$ \mathbf{T} = \begin{bmatrix} j_{11}^2 & j_{21}^2 & j_{31}^2 & j_{11} j_{21} & j_{21} j_{31} & j_{11} j_{31} \\ j_{12}^2 & j_{22}^2 & j_{32}^2 & j_{12} j_{22} & j_{22} j_{32} & j_{12} j_{32} \\ j_{13}^2 & j_{23}^2 & j_{33}^2 & j_{13} j_{23} & j_{23} j_{33} & j_{13} j_{33} \\ 2 j_{11} j_{12} & 2 j_{21} j_{22} & 2 j_{31} j_{32} & j_{12} j_{21} + j_{11} j_{22} & 2 j_{22} j_{31} + j_{21} j_{32} & j_{12} j_{31} + j_{11} j_{32} \\ 2 j_{12} j_{13} & 2 j_{22} j_{23} & 2 j_{32} j_{33} & j_{23} j_{12} + j_{22} j_{13} & j_{23} j_{32} + j_{22} j_{33} & j_{13} j_{32} + j_{12} j_{33} \\ 2 j_{11} j_{13} & 2 j_{21} j_{23} & 2 j_{31} j_{33} & j_{23} j_{11} + j_{21} j_{13} & j_{23} j_{31} + j_{21} j_{33} & j_{13} j_{31} + j_{11} j_{33} \end{bmatrix} \quad (30) $$

with, in case of small displacements and small strains formulation, $j_{ij} = J_{ij}^{-1} = (\mathbf{G}_i \cdot \mathbf{r}_j)^{-1}$ is the inverse of the matrix **J** components such that **J** = **G**<sub>i</sub> ⊗ **r**<sub>j</sub>. Equation (30) expressed on the mid-surface $\xi^3 = 0$ simplifies to $\mathbf{T}_{\xi^3=0}$:

$$ \mathbf{T}_{\xi^3=0} = \begin{bmatrix} j_{11}^2 & j_{21}^2 & 0 & j_{11} j_{21} & 0 & 0 \\ j_{12}^2 & j_{22}^2 & 0 & j_{12} j_{22} & 0 & 0 \\ 0 & 0 & j_{33}^2 & 0 & 0 & 0 \\ 2 j_{11} j_{12} & 2 j_{21} j_{22} & 0 & j_{12} j_{21} + j_{11} j_{22} & 0 & 0 \\ 0 & 0 & 0 & 0 & j_{22} j_{33} & j_{12} j_{33} \\ 0 & 0 & 0 & 0 & j_{21} j_{33} & j_{11} j_{33} \end{bmatrix} \quad (31) $$

However, the calculations required for evaluating the terms of the inverse **J**<sup>-1</sup> of the Jacobian matrix may be laborious and are likely to induce matrix transformations that are potentially computational time consuming, as illustrated in Kim et al. [20]. Nevertheless, introducing the strain smoothing technique in the context of finite element or SFEM technique to calculate membrane and bending components helps avoid such calculations. Indeed, their integration is made in the local Cartesian coordinate system **x** (x, y, z) after smoothing operation. This eliminates the calculation of shape function Cartesian derivatives. This technique is presented in the following sections.

Separating the strain tensor components of Eq. (28) into membrane-bending **ε**<sup>mb</sup>, transverse shear **γ** and transverse normal **ε**<sup>zz</sup>, the linear strain vector can be written, in line with the work proposed in Kim et al. [20], in this compact form:

$$ \boldsymbol{\epsilon} = \begin{Bmatrix} \boldsymbol{\epsilon}^{mb} \\ \boldsymbol{\gamma} \\ \boldsymbol{\epsilon}_{zz} \end{Bmatrix} = \begin{Bmatrix} \boldsymbol{\epsilon}^{m} + \xi^3 \cdot \boldsymbol{\epsilon}^{b} \\ \boldsymbol{\gamma} \\ \boldsymbol{\epsilon}_{zz} \end{Bmatrix} \quad (32) $$

In Eq. (32), **ε**<sup>mb</sup> is the combination of the membrane and the bending strains components, respectively **ε**<sup>m</sup> and **ε**<sup>b</sup> such that

$$ \boldsymbol{\epsilon}^{mb} = \begin{Bmatrix} \boldsymbol{\epsilon}_{xx} \\ \boldsymbol{\epsilon}_{yy} \\ \boldsymbol{\gamma}_{xy} \end{Bmatrix} = \boldsymbol{\epsilon}^{m} + \xi^3 \cdot \boldsymbol{\epsilon}^{b}. \quad \text{Writing Eq. (28) in matrix form:} $$

$$ \boldsymbol{\epsilon}^{m} = \begin{Bmatrix} \frac{\partial \bar{u}}{\partial x} \\ \frac{\partial \bar{v}}{\partial y} \\ \frac{\partial \bar{u}}{\partial y} + \frac{\partial \bar{v}}{\partial x} \end{Bmatrix} = \mathbf{B}^{m} \cdot \tilde{\mathbf{u}} \quad (33) $$

with the membrane strain operator defined as:

$$ \mathbf{B}^m = \begin{bmatrix} \frac{\partial}{\partial x} & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{\partial}{\partial y} & 0 & 0 & 0 & 0 \\ \frac{\partial}{\partial y} & \frac{\partial}{\partial x} & 0 & 0 & 0 & 0 \end{bmatrix} \quad (34) $$

and where:

$$ \boldsymbol{\epsilon}^b = \begin{Bmatrix} \frac{\partial \Delta u}{\partial x} \\ \frac{\partial \Delta v}{\partial y} \\ \frac{\partial \Delta u}{\partial y} + \frac{\partial \Delta v}{\partial x} \end{Bmatrix} = \mathbf{B}^b \cdot \tilde{\mathbf{u}} \quad (35) $$
with the bending strain operator defined as:

$$ \mathbf{B}^b = \begin{bmatrix} 0 & 0 & 0 & \frac{\partial}{\partial x} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{\partial}{\partial y} & 0 \\ 0 & 0 & 0 & \frac{\partial}{\partial y} & \frac{\partial}{\partial x} & 0 \end{bmatrix} \qquad (36) $$

From Eq. (28), the matrix expressions of the transverse shear strain $\gamma$ and the transverse normal strain $\epsilon_{zz}$ are also obtained.

The transverse shear strains $\gamma = \begin{Bmatrix} \gamma_{yz} \\ \gamma_{xz} \end{Bmatrix}$ can be written in matrix form as:

$$ \begin{aligned} \{\gamma_{yz} \gamma_{xz}\} &= \begin{bmatrix} 0 & \frac{\partial \bar{u}}{\partial z} & \frac{\partial \bar{u}}{\partial y} & 0 & \xi^3 \cdot \frac{\partial \Delta v}{\partial z} + \frac{\partial \xi^3}{\partial z} \cdot \Delta v & \frac{\partial \Delta w}{\partial y} \\ \frac{\partial \bar{u}}{\partial z} & 0 & \frac{\partial \bar{u}}{\partial x} & \xi^3 \cdot \frac{\partial \Delta u}{\partial z} + \frac{\partial \xi^3}{\partial z} \cdot \Delta u & 0 & \xi^3 \cdot \frac{\partial \Delta w}{\partial x} \end{bmatrix} \\ &= \mathbf{B}^{\gamma} \cdot \tilde{\mathbf{u}} \end{aligned} \qquad (37) $$

where $\mathbf{B}^{\gamma}$ is:

$$ \mathbf{B}^{\gamma} = \begin{bmatrix} 0 & \frac{\partial}{\partial z} & \frac{\partial}{\partial y} & 0 & \xi^3 \cdot \frac{\partial}{\partial z} + \frac{\partial \xi^3}{\partial z} \xi^3 \cdot \frac{\partial}{\partial y} \\ \frac{\partial}{\partial z} & 0 & \frac{\partial}{\partial x} & \xi^3 \cdot \frac{\partial}{\partial z} + \frac{\partial \xi^3}{\partial z} \xi^3 \cdot \frac{\partial}{\partial x} & 0 \end{bmatrix} \quad (38) $$

The transverse normal strain $\epsilon_{zz}$ can be written in matrix form as:

$$ \epsilon_{zz} = \begin{bmatrix} 0 & 0 & \frac{\partial \bar{w}}{\partial z} & 0 & \xi^3 \cdot \frac{\partial \Delta w}{\partial z} + \frac{\partial \xi^3}{\partial z} \cdot \Delta w \end{bmatrix} = \mathbf{B}^{zz} \cdot \tilde{\mathbf{u}} \quad (39) $$

where $\mathbf{B}^{zz}$ is:

$$ \mathbf{B}^{zz} = \begin{bmatrix} 0 & 0 & \frac{\partial}{\partial z} & 0 & \xi^3 \cdot \frac{\partial}{\partial z} + \frac{\partial \xi^3}{\partial z} \end{bmatrix} \quad (40) $$

Considering Eqs. (33), (35), (37) and (39), the strain tensor can be written:

$$ \boldsymbol{\epsilon} = \mathbf{B} \cdot \tilde{\mathbf{u}} \qquad (41) $$

where $\mathbf{B}$ is a strain operator including membrane-bending, transverse shear and transverse normal components (See Kim et al. [20]).

# 4 Finite element approximation

## 4.1 Strain finite element approximation

The standard trilinear shape functions of the hexahedral element can be defined in terms of standard bilinear shape functions, as in Sze and Yao [62] combined with a parameterization in the thickness direction. The resultant solid-shell geometry is thus defined in terms of shape functions defined on the shell mid-surface combined with a thickness parameterization as:

$$ \mathbf{x} (\xi^1, \xi^2, \xi^3) = \sum_{I=1}^{4} N_I (\xi^1, \xi^2) \cdot (\bar{\mathbf{x}}_I + \xi^3 \cdot \Delta \mathbf{x}_I), \quad (-1 \le \xi^3 \le 1) \qquad (42) $$

where $N_I, I = 1, \dots, 4$ are the standard bilinear shape functions of a quadrilateral element, $\bar{\mathbf{x}}_i = (\bar{x}_{i1}, \bar{x}_{i2}, \bar{x}_{i3}, \bar{x}_{i4})$, ($i = 1, \dots, 3$) are the position vectors of the element mid-surface corners and $\Delta \mathbf{x}_i = (\Delta x_{i1}, \Delta x_{i2}, \Delta x_{i3}, \Delta x_{i4})$, ($i = 1, \dots, 3$) are the orientation vectors at the element mid-surface corners.

The displacements field expressed in Eq. (5) can be approximated in finite element form as $\mathbf{u} \approx \mathbf{u}^e$. For sake of simplicity the superscript ( )$^e$ is omitted in this last expression, such that the nodal displacements field $\mathbf{u}$ is given by:

$$ \begin{aligned} \mathbf{u} (\xi^1, \xi^2, \xi^3) &= \sum_{I=1}^{4} N_I (\xi^1, \xi^2) \cdot \left( \frac{1 + \xi^3}{2} \cdot \mathbf{u}_I^U + \frac{1 - \xi^3}{2} \cdot \mathbf{u}_I^L \right) \\ &= \sum_{I=1}^{4} N_I (\xi^1, \xi^2) \cdot (\bar{\mathbf{u}}_I + \xi^3 \cdot \Delta \mathbf{u}_I), \quad (-1 \le \xi^3 \le 1) \end{aligned} \qquad (43) $$

where $\bar{\mathbf{u}}_I = (\bar{u}_{i1}, \bar{u}_{i2}, \bar{u}_{i3}, \bar{u}_{i4})$, ($i = 1, \dots, 3$) are the displacement vectors of the element mid-surface corners and $\Delta \mathbf{u}_I = (\Delta u_{i1}, \Delta u_{i2}, \Delta u_{i3}, \Delta u_{i4})$, ($i = 1, \dots, 3$) are the rotation vectors at the element mid-surface corners. They are defined at each corner $I$ like in Eq. (5).

The finite element approximation of the membrane and bending strain field, defined respectively in Eqs. (33) and (35), then becomes:

$$ \boldsymbol{\epsilon}^m = \sum_{I=1}^{4} \mathbf{B}_I^m \cdot \tilde{\mathbf{u}}_I \qquad (44) $$

where $\tilde{\mathbf{u}}_I = \langle \bar{u}_I, \bar{v}_I, \bar{w}_I, \Delta u_I, \Delta v_I, \Delta w_I \rangle^T$ for each corner $I$,

with

$$ \mathbf{B}_I^m = \begin{bmatrix} N_{I,x} & 0 & 0 & 0 & 0 & 0 \\ 0 & N_{I,y} & 0 & 0 & 0 & 0 \\ N_{I,y} & N_{I,x} & 0 & 0 & 0 & 0 \end{bmatrix}, \quad I = 1, \dots, 4 \qquad (45) $$

and

$$ \boldsymbol{\epsilon}^b = \sum_{I=1}^{4} \mathbf{B}_I^b \cdot \tilde{\mathbf{u}}_I \qquad (46) $$

with

$$ \mathbf{B}_I^b = \begin{bmatrix} 0 & 0 & 0 & N_{I,x} & 0 & 0 \\ 0 & 0 & 0 & 0 & N_{I,y} & 0 \\ 0 & 0 & 0 & N_{I,y} & N_{I,x} & 0 \end{bmatrix}, \quad I=1, \dots, 4 \qquad (47) $$

The detailed smoothing procedure to calculate this membrane and bending element strain fields, respectively represented by $\boldsymbol{\epsilon}^m$ and $\boldsymbol{\epsilon}^b$, as well as the corresponding strain
operators $\mathbf{B}_I^m$ and $\mathbf{B}_I^b$ in their smoothed form are given in Sect. 6.1.

Developing Eq. (37) in the element mid-surface, at $\xi^3 = 0$, the element transverse shear components can be expressed as:

$$ \gamma = \mathbf{T}_{\xi^3=0}^{\gamma} \cdot \gamma^{\xi} \qquad (48) $$

where the natural transverse shear strain $\gamma^\xi$ in matrix form is obtained by developing Eq. (20) as:

$$ \gamma^{\xi} = \begin{cases} \gamma_{23}^{\xi} \\ \gamma_{13}^{\xi} \end{cases} = \sum_{I=1}^{4} \mathbf{B}_{I}^{\gamma, 0} \cdot \tilde{\mathbf{u}}_{I} \qquad (49) $$

with the transverse shear strain operator $\mathbf{B}_I^{\gamma, 0}$, expressed in the element mid-surface $\xi^3 = 0$ by:

$$ \mathbf{B}_I^{\gamma, 0} = \begin{bmatrix} 0 & 0 & \frac{\partial N_I}{\partial \xi^2} \mathbf{G}_3^T & 0 & N_I \mathbf{G}_2^T & 0 \\ 0 & 0 & \frac{\partial N_I}{\partial \xi^1} \mathbf{G}_3^T & N_I \mathbf{G}_1^T & 0 & 0 \end{bmatrix}, \quad I = 1, \dots, 4 \qquad (50) $$

and $\mathbf{T}_{\xi^3=0}^\gamma$ is a sub-transformation matrix extracted from Eq. (31) as:

$$ \mathbf{T}_{\xi^3=0}^{\gamma} = \begin{bmatrix} j_{22}j_{33} & j_{12}j_{33} \\ j_{21}j_{33} & j_{11}j_{33} \end{bmatrix} \qquad (51) $$

Using Eq. (39), the element normal strain finite element approximation evaluated on the element mid-surface $\xi^3 = 0$ then becomes:

$$ \varepsilon_{zz} = j_{33}^2 \cdot \varepsilon_{33}^{\xi} \qquad (52) $$

where the natural normal strain in matrix form is obtained by developing Eq. (20) as in Eq. (53):

$$ \varepsilon_{33}^{\xi} = \sum_{I=1}^{4} \mathbf{B}_{I}^{\xi\xi} \cdot \tilde{\mathbf{u}}_{I} \qquad (53) $$

where the normal strain operator $\mathbf{B}_I^{\xi\xi}$, evaluated on the element mid-surface $\xi^3 = 0$ is expressed as:

$$ \mathbf{B}_I^{\xi\xi} = \langle 00000 N_I \mathbf{G}_3^T \rangle, \quad I = 1, \dots, 4 \qquad (54) $$

## 4.2 Transverse shear locking

In order to reduce the transverse shear locking when evaluating the transverse shear stiffness using first order shell theory (Reissner-Mindlin assumption), the ANS technique initially introduced by Hughes and Tezduyar [63] and then developed by other researchers like Batoz and Dhatt [64], as Q4-gama element, remains an efficient strategy. Since first order shear deformable shell formulation yields constant transverse shear strain through the shell thickness, this may be a good approximation for very thin shells but not necessarily appropriate for thick shells. Then in order to handle the

<figure><img src="image_5.png" /><figcaption>Fig. 4 Eight-node shell element, ANS method for transverse shear locking</figcaption></figure>

thick shell case, the natural shear strains is evaluated at four points ($A_1, A_2, B_1, B_2$) located in the middle of the mid-surface edge as illustrated in Fig. 4. These four points are used to evaluate the actual transverse shear strain anywhere on the element.

This is done by calculating the element shear strain used in actual shear stiffness computation from the interpolation of the element natural strains defined in Eq. (49) at points $A_1, A_2, B_1$ and $B_2$ as:

$$ \gamma_{13}^{\xi, A_1} = \gamma_{13}^{\xi} \Big|_{\xi^1=0, \xi^2=-1, \xi^3=0} \qquad (55) $$

$$ \gamma_{13}^{\xi, A_2} = \gamma_{13} \Big|_{\xi^1=0, \xi^2=1, \xi^3=0} \qquad (56) $$

$$ \gamma_{23}^{\xi, B_1} = \gamma_{23}^{\xi} \Big|_{\xi^1=-1, \xi^2=0, \xi^3=0} \qquad (57) $$

$$ \gamma_{23}^{\xi, B_2} = \gamma_{23}^{\xi} \Big|_{\xi^1=1, \xi^2=0, \xi^3=0} \qquad (58) $$

Equations (55) to (58) are then used to obtain an improved interpolated shear strain using linear interpolation functions:

$$ \gamma_{13}^{\xi, \text{ANS}} = \frac{1 - \xi^2}{2} \gamma_{13}^{\xi, A_1} + \frac{1 + \xi^2}{2} \gamma_{13}^{\xi, A_2} \qquad (59) $$

$$ \gamma_{23}^{\xi, \text{ANS}} = \frac{1 - \xi^1}{2} \gamma_{23}^{\xi, B_1} + \frac{1 + \xi^1}{2} \gamma_{23}^{\xi, B_2} \qquad (60) $$

where $-1 \le \xi^1, \xi^2 \le 1$.

## 4.3 Trapezoidal effect or curvature thickness effect

In FEM approximation, simulation accuracy is dependent on the mesh quality. However, it remains common to use some trapezoidal elements to mesh curved structures such as, for instance, tubes. In case of bending dominated problems, such a trapezoidal shape may be responsible for the trapezoidal effect occurrence, since it may provide additional through-the-thickness strain. Betsch et al. [1] were among the first to apply the ANS method in order to get rid of this type of locking. This technique tends to reduce the additional thickness strain caused by curved shapes in case of bending situation. Consequently, it also tends to reduce potential locking occur-
<figure><img src="image_6.png" /><figcaption>**Fig. 5** Eight-node shell element, ANS method for trapezoidal effect treatment</figcaption></figure>

rence. The present element has the same improved properties as those detailed in the work of Betsch et al. [1], since the natural thickness strains are evaluated at the four corners of the element mid-surface (C, D, E, F) as illustrated in Fig. 5.

The thickness strain used in actual thickness stiffness calculation is then obtained by interpolation of these natural thickness strains.

The normal strains at the four corners of the mid-surface are then obtained from Eq. (53) as:

$$ \epsilon_{33}^{\xi,C} = \epsilon_{33}^{\xi} \Big|_{\xi^1=-1, \xi^2=-1} \qquad (61) $$

$$ \epsilon_{33}^{\xi,D} = \epsilon_{33}^{\xi} \Big|_{\xi^1=1, \xi^2=-1} \qquad (62) $$

$$ \epsilon_{33}^{\xi,E} = \epsilon_{33}^{\xi} \Big|_{\xi^1=1, \xi^2=1} \qquad (63) $$

$$ \epsilon_{33}^{\xi,F} = \epsilon_{33}^{\xi} \Big|_{\xi^1=-1, \xi^2=1} \qquad (64) $$

And Eqs. (61) to (64) are used to obtain an improved interpolated transverse normal strain using bilinear interpolation functions of the standard mid-surface Q4 element:

$$ \epsilon_{33}^{\xi,ANS} = \frac{1}{4} \cdot (1 - \xi^1) \cdot (1 - \xi^2) \cdot \epsilon_{33}^{\xi,C} \\ + \frac{1}{4} \cdot (1 + \xi^1) \cdot (1 - \xi^2) \cdot \epsilon_{33}^{\xi,D} \\ + \frac{1}{4} \cdot (1 + \xi^1) \cdot (1 + \xi^2) \cdot \epsilon_{33}^{\xi,E} \\ + \frac{1}{4} \cdot (1 - \xi^1) \cdot (1 + \xi^2) \cdot \epsilon_{33}^{\xi,F} \qquad (65) $$

$$ \text{where } -1 \le \xi^1, \xi^2 \le 1. $$

## 5 Principle of virtual work

Following the work of Oden and Reddy [65] and more recently of Simo and Hughes [38], the virtual work principle can be written:

$$ W = W_{\text{int}} - W_{\text{ext}} \qquad (66) $$

where the internal and the external work expressions are given by:

$$ W_{\text{int}} = \int_{V} \boldsymbol{\varepsilon}^{T} : \boldsymbol{\sigma} \, dV \qquad (67) $$

$$ W_{\text{ext}} = \int_{V} \tilde{\mathbf{u}}^{T} \cdot \mathbf{f}_{v} dV + \int_{s} \delta \tilde{\mathbf{u}}^{T} \cdot \mathbf{f}_{s} dS \qquad (68) $$

The weak form expression of the virtual work principle can be written as detailed in Hodge [66]:

$$ \delta W = \delta W_{\text{int}} - \delta W_{\text{ext}} = 0 \qquad (69) $$

with the following definition:

$$ \delta W_{\text{int}} = \int_{V} \delta \boldsymbol{\varepsilon}^{T} : \boldsymbol{\sigma} \, dV \qquad (70) $$

$$ \delta W_{\text{ext}} = \int_{V} \delta \tilde{\mathbf{u}}^{T} \cdot \mathbf{f}_{v} dV + \int_{s} \delta \tilde{\mathbf{u}}^{T} \cdot \mathbf{f}_{s} dS \qquad (71) $$

where the stress field $\boldsymbol{\sigma}$ and the constitutive parameter $\mathbf{H}'$ are assumed to be related by the following constitutive equation:

$$ \boldsymbol{\sigma} = \mathbf{H}' \cdot \boldsymbol{\varepsilon} \qquad (72) $$

where, for an isotropic elastic material

$$ \mathbf{H}' = \begin{bmatrix} \bar{\lambda} + 2\mu & \bar{\lambda} & 0 & 0 & 0 & 0 \\ \bar{\lambda} & \bar{\lambda} + 2\mu & 0 & 0 & 0 & 0 \\ 0 & 0 & \mu & 0 & 0 & 0 \\ 0 & 0 & 0 & \mu & 0 & 0 \\ 0 & 0 & 0 & 0 & \mu & 0 \\ 0 & 0 & 0 & 0 & 0 & E \end{bmatrix} \\ = \begin{bmatrix} \mathbf{H}^{\prime mb} & (0) \\ (0) & \mathbf{H}^{\prime \gamma} \\ \mathbf{H}^{\prime \gamma} & (0) \end{bmatrix} \qquad (73) $$

in which

$$ \mathbf{H}^{\prime mb} = \begin{bmatrix} \bar{\lambda} + 2\mu & \bar{\lambda} & 0 \\ \bar{\lambda} & \bar{\lambda} + 2\mu & 0 \\ 0 & 0 & \mu \end{bmatrix} \qquad (74) $$

$$ \mathbf{H}^{\prime \gamma} = \begin{bmatrix} \mu & 0 \\ 0 & \mu \end{bmatrix} \qquad (75) $$

$$ \mathbf{H}^{\prime \mathrm{zz}} = E \qquad (76) $$

and the modified Lamé's coefficients are $\bar{\lambda} = \frac{E \nu}{(1 - \nu^2)}$ and $\mu = \frac{E}{2(1 + \nu)}$ with E the Young's modulus and $\nu$ the Poisson's ratio.

The constitutive parameter $\mathbf{H}'$, based on the generalized plane stress theory, prevents some coupling between membrane, bending and shear strain modes and the thickness mode. This reduces the risk of extra locking in the thickness direction.

Separating the membrane and bending stresses from the transverse shear and the through-the-thickness stresses, the stress-strain relations can be decomposed as follows:
$$
\begin{align}
\boldsymbol{\sigma}^{mb} &= \mathbf{H}^{\prime mb} \cdot \boldsymbol{\varepsilon}^{mb} = \mathbf{H}^{\prime mb} (\boldsymbol{\varepsilon}^m + \xi^3 \cdot \boldsymbol{\varepsilon}^b) = \boldsymbol{\sigma}^m + \xi^3 \cdot \boldsymbol{\sigma}^b \nonumber \\
\boldsymbol{\sigma}^{\gamma} &= \mathbf{H}^{\prime \gamma} \cdot \boldsymbol{\varepsilon}^{\gamma} \nonumber \\
\boldsymbol{\sigma}^{zz} &= \mathbf{H}^{\prime \prime z} \cdot \boldsymbol{\varepsilon}^{zz} \tag{77}
\end{align}
$$

where $\boldsymbol{\sigma}^m = \mathbf{H}^{\prime mb} \boldsymbol{\varepsilon}^m$ and $\boldsymbol{\sigma}^b = \mathbf{H}^{\prime mb} \boldsymbol{\varepsilon}^b$.

According to the adopted kinematics decomposition for
shell formulation and following Domissy [61] and Hannachi
[21], it is possible to decompose the weak form of the inter-
nal virtual work $\delta W_{\text{int}}$, defined in Eq. (70), into membrane-
bending, transverse shear and normal components as:

$$
\begin{align}
\delta W_{\text{int}} &= \delta W_{\text{int}}^{\text{mb}} + \delta W_{\text{int}}^{\gamma} + \delta W_{\text{int}}^{\text{zz}} \nonumber \\
&= \int_V \left(\delta \boldsymbol{\varepsilon}^{\text{mb}}\right)^T : \boldsymbol{\sigma}^{\text{mb}} dV + \int_V \delta \boldsymbol{\gamma}^T : \boldsymbol{\sigma}^{\gamma} dV \nonumber \\
&\quad + \int_V \left(\delta \boldsymbol{\varepsilon}^{\text{zz}}\right)^T : \boldsymbol{\sigma}^{\text{zz}} dV \tag{78}
\end{align}
$$

Introducing Eq. (77) in Eq. (78) gives:

$$
\begin{align}
\delta W_{\text{int}} &= \delta W_{\text{int}}^{\text{mb}} + \delta W_{\text{int}}^{\gamma} + \delta W_{\text{int}}^{\text{zz}} \nonumber \\
&= \int_V \left(\delta \boldsymbol{\varepsilon}^{\text{mb}}\right)^T \cdot \mathbf{H}^{\prime \text{mb}} \cdot \boldsymbol{\varepsilon}^{\text{mb}} dV + \int_V \delta \boldsymbol{\gamma}^T \cdot \mathbf{H}^{\prime \gamma} \cdot \boldsymbol{\varepsilon}^{\gamma} dV \nonumber \\
&\quad + \int_V \left(\delta \varepsilon_{zz}\right)^T \cdot \mathbf{H}^{\prime \prime z} \cdot \varepsilon_{zz}^d V \tag{79}
\end{align}
$$

Then, using Eq. (79) and separating integration into thickness
and mid-surface integrations gives:

$$
\begin{align}
\delta W_{\text{int}} &= \delta W_{\text{int}}^{\text{mb}} + \delta W_{\text{int}}^{\gamma} + \delta W_{\text{int}}^{\text{zz}} \nonumber \\
&= \int_S (\delta \boldsymbol{\varepsilon}^m)^T \cdot \mathbf{A} \mathbf{H} \cdot \boldsymbol{\varepsilon}^m dS + \int_S (\delta \boldsymbol{\varepsilon}^b)^T \cdot \mathbf{D} \mathbf{H} \cdot \boldsymbol{\varepsilon}^b dS + \dots \nonumber \\
&\quad + \int_S \delta \boldsymbol{\gamma}^T \cdot \bar{\mathbf{A}} \mathbf{H} \cdot \boldsymbol{\gamma} dS + \int_S (\delta \varepsilon_{zz})^T \cdot \mathbf{E} \mathbf{H} \cdot \varepsilon_{zz}^d S \nonumber \\
&= \delta \tilde{\mathbf{u}}^T \cdot \left[ \int_S (\mathbf{B}^m)^T \cdot \mathbf{A} \mathbf{H} \cdot \mathbf{B}^m dS + \int_S (\mathbf{B}^b)^T \cdot \mathbf{D} \mathbf{H} \cdot \mathbf{B}^b dS + \dots \right. \nonumber \\
&\qquad \left. + \int_S (\mathbf{B}^\gamma)^T \cdot \bar{\mathbf{A}} \mathbf{H} \cdot \mathbf{B}^\gamma dS + \int_S (\mathbf{B}^{\prime \prime z})^T \cdot \mathbf{E} \mathbf{H} \cdot \mathbf{B}^{\prime \prime z} dS \right] \cdot \tilde{\mathbf{u}} \tag{80}
\end{align}
$$

where *h* is the shell thickness and **A**H, **D**H, $\bar{\mathbf{A}}\mathbf{H}$ and **E**H are
through-the-thickness integrated material parameters thus
defining the new resultant constitutive law expressed as:

$$
\mathbf{A} \mathbf{H} = \int_{-h/2}^{h/2} \mathbf{H}_{ij}^{\prime mb} d\xi^3 = h \cdot \mathbf{H}^{\prime mb} \quad (81)
$$

$$
\mathbf{D} \mathbf{H} = \int_{-h/2}^{h/2} (\xi^3)^2 \cdot \mathbf{H}^{\prime mb} d\xi^3 = \frac{h^3}{12} \cdot \mathbf{H}^{\prime mb} \quad (82)
$$

$$
\bar{\mathbf{A}}\mathbf{H} = \int_{-h/2}^{h/2} \mathbf{H}^{\prime \gamma} d\xi^3 = h \cdot \mathbf{H}^{\prime \gamma} \quad (83)
$$

$$
EH = \int_{-h/2}^{h/2} H'^{zz} d\xi^3 = h \cdot H'^{zz} \quad (84)
$$

Introducing the finite element approximation, the element
internal virtual work $\delta W_{\text{int}}^e$ in Eq. (78) can be expressed as:

$$
\begin{align}
\delta W_{\text{int}}^{e} &= \delta W_{\text{int}}^{e,mb} + \delta W_{\text{int}}^{e,\gamma} + \delta W_{\text{int}}^{e,zz} \nonumber \\
&= \int_{-1}^{1} \int_{-1}^{1} \sum_{I=1}^{4} (\delta \boldsymbol{\varepsilon}_{I}^{m} + \xi^{3} \cdot \delta \boldsymbol{\varepsilon}_{I}^{b})^{T} : \hat{\boldsymbol{\sigma}}^{mb} |\mathbf{J}| d\xi^{1} d\xi^{2} + \dots \nonumber \\
&\quad + \int_{-1}^{1} \int_{-1}^{1} \sum_{I=1}^{4} \delta \boldsymbol{\gamma}_{I}^{T} : \hat{\boldsymbol{\sigma}}^{\gamma} |\mathbf{J}| d\xi^{1} d\xi^{2} \nonumber \\
&\quad + \int_{-1}^{1} \int_{-1}^{1} \sum_{I=1}^{4} (\delta \boldsymbol{\varepsilon}_{zz})^{T} : \hat{\boldsymbol{\sigma}}^{zz} |\mathbf{J}| d\xi^{1} d\xi^{2} \nonumber \\
&= \delta \tilde{\mathbf{u}}_{I}^{T} \cdot \int_{-1}^{1} \int_{-1}^{1} \left[ \sum_{I=1}^{4} (\mathbf{B}_{I}^{m} + \xi^{3} \cdot \mathbf{B}_{I}^{b})^{T} \cdot \hat{\boldsymbol{\sigma}}^{mb} |\mathbf{J}| + \dots \right. \nonumber \\
&\qquad \left. + \sum_{I=1}^{4} (\mathbf{B}_{I}^{\gamma})^{T} \cdot \hat{\boldsymbol{\sigma}}^{\gamma} |\mathbf{J}| + \sum_{I=1}^{4} (\mathbf{B}_{I}^{zz})^{T} \cdot \hat{\boldsymbol{\sigma}}^{zz} |\mathbf{J}| \right] d\xi^{1} d\xi^{2} \nonumber \\
&= \delta \tilde{\mathbf{u}}_{I}^{T} \cdot \mathbf{f}_{\text{int}} \tag{85}
\end{align}
$$

where $|\mathbf{J}|$ is the determinant of $\mathbf{J}$, the resultant stresses
are $\hat{\boldsymbol{\sigma}}^{mb} = \begin{bmatrix} \mathbf{A} \mathbf{H} \\ \mathbf{D} \mathbf{H} \end{bmatrix} \cdot \begin{bmatrix} \boldsymbol{\varepsilon}_{I}^{m} \\ \boldsymbol{\varepsilon}_{I}^{b} \end{bmatrix}, \hat{\boldsymbol{\sigma}}^{\gamma} = \bar{\mathbf{A}} \mathbf{H} \cdot \boldsymbol{\gamma}_{I}, \hat{\boldsymbol{\sigma}}^{zz} = \mathbf{E} \mathbf{H} \cdot \boldsymbol{\varepsilon}_{zz, I}$ and

$$
\mathbf{f}_{\text{int}} = \int_{-1}^{1} \int_{-1}^{1} \left[ (\mathbf{B}_{I}^{m} + \xi^3 \cdot \mathbf{B}_{I}^{b})^{T} \cdot \hat{\boldsymbol{\sigma}}^{mb} + (\mathbf{B}_{I}^{\gamma})^{T} \cdot \hat{\boldsymbol{\sigma}}^{\gamma} + (\mathbf{B}_{I}^{zz})^{T} \cdot \hat{\boldsymbol{\sigma}}^{zz} \right] |\mathbf{J}| d\xi^1 d\xi^2 \quad (86)
$$

The finite element approximation of the weak form of the
internal virtual work $\delta W_{\text{int}}^e$ expressed in Eq. (80) then yields
(See Kim et al. [20]):

$$
\begin{align}
\delta W_{\text{int}}^{e} &= \delta \tilde{\mathbf{u}}_{I}^{T} \cdot \int_{-1}^{1} \int_{-1}^{1} \left[ \sum_{I=1}^{4} (\mathbf{B}_{I}^{m})^{T} \cdot \mathbf{A} \mathbf{H} \cdot \mathbf{B}_{I}^{m} \cdot |\mathbf{J}| \right. \nonumber \\
&\qquad \left. + \sum_{I=1}^{4} (\mathbf{B}_{I}^{b})^{T} \cdot \mathbf{D} \mathbf{H} \cdot \mathbf{B}_{I}^{b} \cdot |\mathbf{J}| + \dots \right. \nonumber \\
&\qquad \left. + \sum_{I=1}^{4} (\mathbf{B}_{I}^{\gamma})^{T} \cdot \bar{\mathbf{A}} \mathbf{H} \cdot \mathbf{B}_{I}^{\gamma} \cdot |\mathbf{J}| \right. \nonumber \\
&\qquad \left. + \sum_{I=1}^{4} (\mathbf{B}_{I}^{zz})^{T} \cdot \mathbf{E} \mathbf{H} \cdot \mathbf{B}_{I}^{zz} \cdot |\mathbf{J}| \right] d\xi^{1} d\xi^{2} \cdot \Delta \tilde{\mathbf{u}}_{I} \nonumber \\
&= \delta \tilde{\mathbf{u}}_{I}^{T} \cdot \mathbf{k} \cdot \Delta \tilde{\mathbf{u}}_{I} \tag{87}
\end{align}
$$
where the element stiffness matrix **k** is:

$$ \mathbf{k} = \mathbf{k}^{mb} + \mathbf{k}^{\gamma} + \mathbf{k}^{zz} \qquad (88) $$

with

$$ \mathbf{k}^{mb} = \int_{-1}^{1} \int_{-1}^{1} \sum_{I=1}^{4} \left[ \left( \mathbf{B}_{I}^{m} \right)^{T} \cdot \mathbf{A} \mathbf{H} \cdot \mathbf{B}_{I}^{m} + \left( \mathbf{B}_{I}^{b} \right)^{T} \cdot \mathbf{D} \mathbf{H} \cdot \mathbf{B}_{I}^{b} \right] \cdot |\mathbf{J}| d\xi^{1} d\xi^{2} \qquad (89) $$

$$ \mathbf{k}^{\gamma} = \int_{-1}^{1} \int_{-1}^{1} \sum_{I=1}^{4} (\mathbf{B}_{I}^{\gamma})^{T} \cdot \bar{\mathbf{A}} \mathbf{H} \cdot \mathbf{B}_{I}^{\gamma} \cdot |\mathbf{J}| d\xi^{1} d\xi^{2} \qquad (90) $$

$$ \mathbf{k}^{zz} = \int_{-1}^{1} \int_{-1}^{1} \sum_{I=1}^{4} (\mathbf{B}_{I}^{zz})^{T} \cdot \mathbf{E} \mathbf{H} \cdot \mathbf{B}_{I}^{zz} \cdot |\mathbf{J}| d\xi^{1} d\xi^{2} \qquad (91) $$

# 6 Smoothed finite element method

## 6.1 Smoothed strain field FEM formulation

At this stage, the standard FEM and the SFEM formulations share the same shell kinematics and variational principles. A fundamental difference will however result from the introduction, in the SFEM, of the smoothed strain obtained from a strain smoothing operation introduced in Liu et al. [35], Dai and Liu [42] or Nguyen-Thanh et al. [49]. This smoothed strain is defined as:

$$ \tilde{\varepsilon}_{ij}(\mathbf{x}_c) = \int_{\Omega_c^h} \varepsilon_{ij}(\mathbf{x}) \cdot \varphi(\mathbf{x} - \mathbf{x}_c) d\Omega \qquad (92) $$

where $\varphi$ is a smoothing function with $\varphi(\mathbf{x} - \mathbf{x}_c) = \begin{cases} 1/S_c & \mathbf{x} \in \Omega_c \\ 0 & \mathbf{x} \notin \Omega_c \end{cases}$ and $S_c = \int_{\Omega_c} d\Omega$ is the area of the smoothing cell $\Omega_c$ with $\mathbf{x}_c$ the position vector of a point on the shell structure. This operation is very similar to the mean dilatation procedure used to deal with the incompressibility in non-linear mechanics. It has also been used in the weak-form meshless method based on nodal integration or in the moving particle finite element method (MPFEM) developed in Hao et al. [67].

Applying the strain smoothing concepts to the deformation gradient, the smoothed deformation gradient can be written as:

$$ \tilde{F}_{ij}(\mathbf{x}_c) = \frac{1}{S_c} \int_{\Omega_c} \left[ \frac{\partial u_i}{\partial x_j} + \delta_{ij} \right] d\Omega \qquad (93) $$

By applying the divergence theorem, $\int_{\Omega_c} \left[ \frac{\partial u_i}{\partial x_j} \right] d\Omega = \int_{\Gamma_c} (u_i n_j) d\Gamma$, the smoothed strain formulation is then expressed as an integration of the strain matrix on the cell boundary and Eq. (93) becomes:

$$ \tilde{F}_{ij}(\mathbf{x}_c) = \tilde{\varepsilon}_{ij}(\mathbf{x}_c) + \delta_{ij} \qquad (94) $$

<figure><img src="image_7.png" /><figcaption>Fig. 6 Smoothing concept: mid-surface identification</figcaption></figure>

where

$$ \tilde{\varepsilon}_{ij}(\mathbf{x}_c) = \frac{1}{S_c} \int_{\Gamma_c} (u_i \cdot n_j) d\Gamma \qquad (95) $$

Hence the associated smoothed Green–Lagrange strain $\tilde{\varepsilon}$ becomes:

$$ \begin{aligned} \tilde{\varepsilon}_{ij} &= \frac{1}{2} (\tilde{F}_{ik} \tilde{F}_{kj} - I_{ij}) \\ &= \frac{1}{2} [(\tilde{\varepsilon}_{ik} + \delta_{ik}) (\tilde{\varepsilon}_{kj} + \delta_{kj}) - \delta_{ij}] \end{aligned} \qquad (96) $$

In the previous equation, $\Gamma_c$ is the cell boundary and $n_j$ the outward normal of the considered cell boundary.

In the concept presented here, the mid-surface of the shell is first identified. Then, the length of its boundary $\Gamma_c$, the normal vectors of its boundary $\mathbf{n}$ and its area $S_c$ are determined as illustrated in Fig. 6.

Considering only the linear part and after simplification, the smoothed strain becomes:

$$ \tilde{\varepsilon}_{ij}^{l} = \frac{1}{2} [\tilde{\varepsilon}_{ij} + \tilde{\varepsilon}_{ji}] \qquad (97) $$

Introducing the local displacements field $\mathbf{u}(\xi^1, \xi^2, \xi^3) = \bar{\mathbf{u}}(\xi^1, \xi^2) + \xi^3 \cdot \Delta\mathbf{u}(\xi^1, \xi^2)$, the deformation gradient can then be written as:

$$ \begin{aligned} \tilde{F}_{ij}(\mathbf{x}_c) &= \frac{1}{S_c} \int_{\Omega_c} \left[ \frac{\partial [\bar{u}_i(\xi^1, \xi^2) + \xi^3 \cdot \Delta u_i(\xi^1, \xi^2)]}{\partial x_j} \right] + \delta_{ij} d\Omega \\ &= \tilde{\varepsilon}_{ij}(\mathbf{x}_c) + \delta_{ij} \end{aligned} \qquad (98) $$

with

$$ \tilde{\varepsilon}_{ij}(\mathbf{x}_c) = \frac{1}{S_c} \int_{\Gamma_c} \left( \bar{u}_i \cdot \mathbf{n}_j + \xi^3 \Delta u_i \mathbf{n}_j \right) d\Gamma \qquad (99) $$

Here, the smoothed element formulation is applied for the numerical integration of the membrane and bending parts of the strain field on the mid-surface boundary for the resultant eight-node solid-shell element.

Separating the strain field into membrane-bending, transverse shear and transverse normal components then gives:

$$ \tilde{\varepsilon}_{ij} = \tilde{\varepsilon}_{ij}^{m} + \tilde{\varepsilon}_{ij}^{b} + \varepsilon_{ij}^{\gamma} + \varepsilon_{ij}^{zz} \qquad (100) $$

From Eqs. (89) and (99) it is possible to deduce the continuous smoothed form of the membrane and bending stiffness matrix $\tilde{\mathbf{k}}^{mb}$:
<figure><img src="image_8.png" /><figcaption>**Fig. 7** Mid-surface partition in cells</figcaption></figure>

$$ \tilde{\mathbf{k}}^{mb} = \sum_{c=1}^{n_c} \left[ \begin{aligned} & (\tilde{\mathbf{B}}_c^m(x_c))^T \cdot \mathbf{A} \mathbf{H} \cdot \tilde{\mathbf{B}}_c^m(x_c) \\ & + (\tilde{\mathbf{B}}_c^b(x_c))^T \cdot \mathbf{D} \mathbf{H} \cdot \tilde{\mathbf{B}}_c^b(x_c) \end{aligned} \right] \cdot S_c \quad (101) $$

where the smoothed strain-displacement matrix operators are determined similarly as in Nguyen-Thanh et al. [49]:

$$ \tilde{\mathbf{B}}_c^m(x_c) = \frac{1}{S_c} \cdot \sum_{b=1}^{4} \sum_{I=1}^{4} \begin{bmatrix} N_I \cdot n_x & 0 & 0 & 0 & 0 \\ 0 & N_I \cdot n_y & 0 & 0 & 0 \\ N_I \cdot n_y & N_I \cdot n_x & 0 & 0 & 0 \\ 0 & 0 & N_I \cdot n_x & 0 & 0 \\ 0 & 0 & 0 & N_I \cdot n_y & 0 \end{bmatrix} \cdot l_b^c \quad (102) $$

$$ \tilde{\mathbf{B}}_c^b(x_c) = \frac{1}{S_c} \cdot \sum_{b=1}^{4} \sum_{I=1}^{4} \begin{bmatrix} 0 & 0 & 0 & N_I \cdot n_x & 0 \\ 0 & 0 & 0 & 0 & N_I \cdot n_y \\ 0 & 0 & 0 & N_I \cdot n_y & N_I \cdot n_x \\ 0 & 0 & 0 & N_I \cdot n_y & N_I \cdot n_x \\ 0 & 0 & 0 & N_I \cdot n_x & 0 \end{bmatrix} \cdot l_b^c \quad (103) $$

In Eqs. (101) to (103), *I* is the mid-surface node indices from 1 to 4; *x*<sub>c</sub> is the unique gauss point along the element edge and *l*<sub>b</sub><sup>c</sup> is the length of Γ<sub>b</sub><sup>c</sup>; *n*<sub>x</sub> and *n*<sub>y</sub> are the components of the cell border outward normal vector **n**; *S*<sub>c</sub> is the area of cell *c*; *n*<sub>c</sub> is the total number of cells; *N*<sub>I</sub> is the shape functions at the integration point and *b* is the edge number of cell *c* going from 1 to 4. One of the main advantages of the SFEM is that *N*<sub>I</sub>, *l*<sub>b</sub><sup>c</sup>, **n** and *S*<sub>c</sub> are directly calculated from the real structure geometry as presented in Fig. 7.

As a result, the expression of the total element stiffness matrix expressed in Eq. (88), can take the smoothed form in the local Cartesian coordinate system **x** (*x*, *y*, *z*):

$$ \tilde{\mathbf{k}} = \tilde{\mathbf{k}}^{mb} + \mathbf{k}^{\gamma} + \mathbf{k}^{zz} \quad (104) $$

The approximation of the shear strain field presented in section 4.2 is used to determine the shear strain matrix **B**<sup>ANS,γ</sup>. This matrix is then used in the numerical integration of the transverse shear stiffness at four classical Gauss points. A 2 × 2 integration scheme is adopted in the element mid-surface. The transverse shear stiffness, evaluated at the geometry mid-surface and free of transverse shear locking, is then calculated introducing Eqs. (59) and (60) into Eq. (49) to give:

$$ \mathbf{k}^{\text{ANS}, \gamma} = \sum_{I=1}^{4} (\mathbf{B}_{I}^{\text{ANS}, \gamma})^{T} \cdot \bar{\mathbf{A}} \mathbf{H} \cdot \mathbf{B}_{I}^{\text{ANS}, \gamma} \cdot |\mathbf{J}| \quad (105) $$

Then, the approximation of the transverse normal strain field presented in Sect. 4.3 is used to determine the transverse normal strain matrix **B**<sup>ANS,zz</sup>. It is this latter which is used in the numerical integration of the through-the-thickness normal stiffness Eq. (65) defined as:

$$ \mathbf{k}^{\text{ANS,zz}} = \sum_{I=1}^{4} (\mathbf{B}_{I}^{\text{ANS,zz}})^{T} \cdot \mathbf{E} \mathbf{H} \cdot \mathbf{B}_{I}^{\text{ANS,zz}} \cdot |\mathbf{J}| \quad (106) $$
<figure><img src="image_9.png" /><figcaption>Fig. 8 Membrane/bending integration scheme per element; a classical and b SFEMs</figcaption></figure>

The adopted integration scheme is $2 \times 2$.

The total local element stiffness matrix then becomes:

$$ \tilde{\mathbf{k}}^{ANS} = \tilde{\mathbf{k}}^{m} + \tilde{\mathbf{k}}^{b} + \mathbf{k}^{ANS,\gamma} + \mathbf{k}^{ANS,zz} \quad (107) $$

The local to global stiffness matrix transformation is similar to the method used in Kim et al. [20].

## 6.2 Operation count

A quick overview of the integration scheme used in the SFEM is given in this section. It is compared with the integration scheme used in a classical eight-node resultant shell element to deal with the membrane and bending effects as illustrated in Fig. 8.

In the classical scheme, the intrinsic shape function and shape function derivatives are calculated at the four integration points located in the element mid-surface. Then the Jacobian matrix **J** and its determinant |**J**| are calculated. As the element is integrated in local Cartesian coordinate system **x** (x, y, z), the Cartesian shape function derivatives, required to define the strain operators **B**<sup>m</sup> and **B**<sup>b</sup>, are computed by multiplying the intrinsic derivative with the inverse of the jacobian matrix **J**<sup>-1</sup>. Then the stiffness matrix **K**<sup>mb</sup> is calculated at each integration point and assembled with the corresponding element vector and matrices.

In the SFEM, there is no need to calculate the intrinsic shape function derivatives. At the beginning of the integration, each element is subdivided into $n_c$ cells and over each cell a smoothed strain is used, which is a constant integrand,

and the integral over each cell is evaluated over the cell boundary. Geometric parameters such as the length of the four edges of each cells $l_b^c$, the cell areas $S_c$ and the cell edge normals **n** are calculated depending on the number of partitioning cells used within the element. In the present study, even if the number of cells remains limited to one, two or four, it is possible to increase this number depending on the required level of accuracy. In order to get familiar with the calculated geometric parameters, let us consider as an exercise, one cell for the element integration.

As presented in Sect. 6.1, only lengths $l_1^1, l_2^1, l_3^1, l_4^1$ and normals $\mathbf{n}_1^1, \mathbf{n}_2^1, \mathbf{n}_3^1, \mathbf{n}_4^1$, of four edges are calculated as well as one area $A_1$, all of them calculated in the local Cartesian coordinate system **x** (x, y, z). Mathematically, the length of the four edges of each cell $l_b^c$ corresponds to the Jacobian used in the loop on cell edges integration. The intrinsic shape functions are calculated in the center of each cell boundary. Then, the Jacobian matrix **J** and its determinant |**J**| are evaluated in the center of each cell located in the element mid-surface. As no Cartesian derivatives of the shape function are calculated, there is no need to calculate the inverse of the Jacobian matrix **J**<sup>-1</sup> to determine the smoothed operator $\tilde{\mathbf{B}}^{mb}$ as presented in Eqs. (102) and (103). Then the smoothed stiffness matrix $\tilde{\mathbf{K}}^{mb}$ is calculated and summed on each cell if more than one cell is used for the integration. The smoothing method simplifies the expression of the integration terms and is very attractive, especially when only one cell is used, since the number of operations remains small.

## 6.3 Element consistency

Taking advantage of the divergence theorem presented in Sect. 6.1, the relations between the classic and the smoothed strain-displacement matrix operators are easily determined (See Liu et al. [36] for details):

$$ \tilde{\mathbf{B}}_{cl}^{m} = \frac{1}{S_C} \int_{\Omega_C} \mathbf{B}_{l}^{m}(\mathbf{x}) \, d\Omega \quad (108) $$

$$ \tilde{\mathbf{B}}_{cl}^{b} = \frac{1}{S_C} \int_{\Omega_C} \mathbf{B}_{l}^{b}(\mathbf{x}) \, d\Omega \quad (109) $$

Indeed, the smoothed strain operators of Eqs. (102) and (103) are equivalent to the strain operators of Eqs. (45) and (47), smoothed over a given smoothing cell. The solution of the present resultant solid-shell element is then different than the solution of a resultant element using a classical FEM. This can be observed in Sect. 7.

Considering the smoothed strain field $\tilde{\mathbf{e}}$, the orthogonal condition takes an analogous form to that in Liu et al. [36]:

$$ \begin{aligned} \int_V \delta \bar{\sigma}^T : (\nabla_X \tilde{\mathbf{u}} - \tilde{\mathbf{e}}) dV + \int_V \bar{\sigma}^T : (\nabla_X \delta \tilde{\mathbf{u}} - \delta \tilde{\mathbf{e}}) dV &= 0 \\ \Rightarrow \int_V \delta \tilde{\mathbf{e}}^T \cdot \mathbf{H}' \cdot \nabla_X \tilde{\mathbf{u}} dV - \int_V \delta \tilde{\mathbf{e}}^T \cdot \mathbf{H}' \cdot \tilde{\mathbf{e}} dV \end{aligned} $$
$$
\begin{align*}
&+ \int_V \tilde{\boldsymbol{\epsilon}}^T \cdot \boldsymbol{H}' \cdot \nabla_X \delta \tilde{\boldsymbol{u}} \, dV - \int_V \tilde{\boldsymbol{\epsilon}}^T \cdot \boldsymbol{H}' \cdot \delta \tilde{\boldsymbol{\epsilon}} \, dV = 0 \\
&\Rightarrow \delta \tilde{\boldsymbol{u}}^T \left( \int_V \tilde{\boldsymbol{B}}^T \cdot \boldsymbol{H}' \cdot \boldsymbol{B} \, dV \right) \tilde{\boldsymbol{u}} - \delta \tilde{\boldsymbol{u}}^T \left( \int_V \tilde{\boldsymbol{B}}^T \cdot \boldsymbol{H}' \cdot \tilde{\boldsymbol{B}} \, dV \right) \tilde{\boldsymbol{u}} \\
&\quad + \tilde{\boldsymbol{u}}^T \left( \int_V \tilde{\boldsymbol{B}}^T \cdot \boldsymbol{H}' \cdot \boldsymbol{B} \, dV \right) \delta \tilde{\boldsymbol{u}} - \dots \\
&\quad - \tilde{\boldsymbol{u}}^T \left( \int_V \tilde{\boldsymbol{B}}^T \cdot \boldsymbol{H}' \cdot \tilde{\boldsymbol{B}} \, dV \right) \delta \tilde{\boldsymbol{u}} = 0 \\
&\Rightarrow \int_V \tilde{\boldsymbol{B}}^T \cdot \boldsymbol{H}' \cdot \boldsymbol{B} \, dV = \int_V \tilde{\boldsymbol{B}}^T \cdot \boldsymbol{H}' \cdot \tilde{\boldsymbol{B}} \, dV \tag{110}
\end{align*}
$$

with $\tilde{\mathbf{B}}$, a smoothed strain operator similar to that of Eq. (41), but replacing the membrane and bending components with their smoothed form and where $\tilde{\mathbf{e}} = \tilde{\mathbf{B}} \cdot \tilde{\mathbf{u}}$.

From Eq. (110) it can be easily demonstrated that the SFEM formulation remains consistent if one smoothing cell is used within the element. (See proof of Theorem 2 in Liu et al. [36]). However, in such a situation, they have shown that the solution of the SFEM has the same properties as for a reduced integration element, initially developed by Zienkiewicz et al. [68]. In this type of element, only one integration point is used. It is located at the centroid of the isoparametric element, where $\xi^1 = \xi^2 = \xi^3 = 0$. The strain operator $\mathbf{B}$ is then equivalent to an average value over the isoparametric element mid-surface. Similarly, the smoothed strain operator $\tilde{\mathbf{B}}$ is an average value over the element mid-surface. This could result in the apparition of zero energy modes called “hourglass modes” in the stiffness matrix.

If the number of smoothing cells is at least two, the SFEM element does not stay variationally consistent as demonstrated in Theorem 4 of Liu et al. [36], but energy consistent. In this case neither the stress continuity at the cell interfaces nor the equilibrium within the element are guaranteed. However, these configurations have the advantage of not having hourglass modes.

## 7 Benchmark problems

The proposed eight-node solid-shell element is named RH8s-X (Resultant hexahedra eight nodes smoothed with X smoothing cells) where X represents the number of smoothing cells dividing its mid-surface, 1, 2 or 4. Its capabilities have been compared in typical benchmark problems with other existing elements presented in the following list:

*   **SH8**—eight-node element with ANS and EAS equivalent to the SCH8 element from Hannachi [21].

*   **SC8R**—eight-node element with reduced integration and plane stress assumption Abaqus user's manual, ver 6.8 [69].

*   **RH8**—resultant eight-node element with ANS.

<figure><img src="image_10.png" /><figcaption>Fig. 9 Cook's membrane problem</figcaption></figure>

*   **Xsolid85**—resultant eight-node element with ANS and plane stress assumption Kim et al. [20].

*   **RH8s-X** (X = 1, 2 or 4)—present eight-node resultant element.

*   **Sch2009**—Schwarze and Reese [23] : reduced integration eight—node solid-shell element with ANS and EAS.

*   **MISTk** (k = 1, 2 or 4) – Nguyen-Xuan [52] : four—node shell element with smoothed membrane and bending strains.

All results have been normalized to the exact reference solution of each problem.

### 7.1 Cook's membrane problem

This benchmark problem is generally used to verify the ability of an element to handle membrane problems. The irregular mesh also helps verifying the accuracy of the isoparametric Jacobian transformation matrix calculation at the same time. A plate, with the geometry presented in Fig. 9, is clamped on the left side and subjected to a shear load $F = 1$ on the right side. The Young's modulus is $E = 1$ and the Poisson's ratio is $\nu = 0.33$. The vertical displacement of point A, located at the right side center, is compared to the exact reference solution $U_y = 23.91$ in Fig. 10 and Table 2.

From Table 2 it is seen that when the number of smoothing cells increases from 1 to 4, the displacement of point A decreases progressively from an overestimated to an underestimated solution approaching the result obtained with the classical RH8 element. For a very coarse mesh, at least two cells seem to be necessary to arrive at a sufficiently accurate solution to the problem. The current element shows good performance. Even if the SC8R element results have the best convergence rate towards the reference result, the RH8s-1 and RH8s-2 remain the most accurate elements.
<figure><img src="image_11.png" /><figcaption>**Fig. 10** Cook's membrane problem normalized results</figcaption></figure>

## 7.2 Pinched cylinder with end diaphragm problem

In the pinched cylinder problem, a cylinder with an end diaphragm is subjected to a concentrated load at its center

top surface A. Its geometry parameters, presented in Fig. 11a are a length L = 600, a radius R = 25 and a thickness t = 0.25. Due to the symmetry of the model, only a quarter of the cylinder is considered in simulations. The Young's modulus is $E = 3e-6$ and the Poisson's ratio is $\nu = 0.3$. This typical benchmark problem is very useful to determine the capacity of an element to deal with inextentional bending modes and complex membrane states. The vertical displacement of point A is compared to the reference solution $U_y = 1.8248 e^{-5}$. Irregular meshes, as in Fig. 11b, created using a technique presented in Nguyen-Thanh et al. [49] are also studied. Results for different mesh sizes are presented in Table 3 and in Fig. 12.

From Table 3, it is seen that when the number of smoothing cells is increased from 1 to 4 in the RH8s-X element, the displacement of point A decreases progressively, approaching the result obtained with the classical RH8 element. The same observation can be made for both regular and irregular meshes. In terms of result accuracy, the proposed element remains the best for regular and irregular meshes among the resultant shell elements. Thus, RH8s-1 is here again the best

<table><caption>**Table 2** The normalized results of Cook membrane problem</caption><thead><tr><th>Element number</th><th>SC8R</th><th>Xsolid85</th><th>SH8</th><th>RH8</th><th>RH8s-1</th><th>RH8s-2</th><th>RH8s-4</th></tr></thead><tbody><tr><td>2 × 2</td><td>1.112</td><td>0.495</td><td>0.492</td><td>0.492</td><td>1.259</td><td>0.784</td><td>0.547</td></tr><tr><td>4 × 4</td><td>0.963</td><td>0.765</td><td>0.763</td><td>0.763</td><td>1.032</td><td>0.928</td><td>0.810</td></tr><tr><td>8 × 8</td><td>0.974</td><td>0.953</td><td>0.922</td><td>0.922</td><td>1.008</td><td>0.980</td><td>0.941</td></tr><tr><td>16 × 16</td><td>0.987</td><td>0.98</td><td>0.979</td><td>0.979</td><td>1.003</td><td>0.995</td><td>0.985</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Fig. 11** Pinched cylinder with end diaphragm problem</figcaption></figure>
<table><caption>**Table 3** The normalized results of pinched cylinder with end diaphragm problem</caption><thead><tr><th colspan="11">Regular mesh</th></tr><tr><th>Element number</th><th>SC8R</th><th>Mist1</th><th>Mist2</th><th>Mist4</th><th>Xsolid85</th><th>SH8</th><th>RH8</th><th>RH8s-1</th><th>RH8s-2</th><th>RH8s-4</th></tr></thead><tbody><tr><td>4 × 4</td><td>0.453</td><td>0.475</td><td>0.442</td><td>0.388</td><td>0.382</td><td>0.383</td><td>0.375</td><td>0.478</td><td>0.417</td><td>0.392</td></tr><tr><td>8 × 8</td><td>0.785</td><td>0.809</td><td>0.788</td><td>0.755</td><td>0.751</td><td>0.75</td><td>0.746</td><td>0.81</td><td>0.782</td><td>0.758</td></tr><tr><td>16 × 16</td><td>0.945</td><td>0.957</td><td>0.948</td><td>0.935</td><td>0.932</td><td>0.935</td><td>0.929</td><td>0.955</td><td>0.945</td><td>0.934</td></tr><tr><td>32 × 32</td><td>0.995</td><td>—</td><td>—</td><td>—</td><td>0.991</td><td>0.989</td><td>0.989</td><td>0.997</td><td>0.993</td><td>0.99</td></tr></tbody><thead><tr><th colspan="11">Irregular mesh</th></tr><tr><th>Element number</th><th>SC8R</th><th>SH8</th><th>RH8</th><th>RH8s-1</th><th>RH8s-2</th><th>RH8s-4</th></tr></thead><tbody><tr><td>4 × 4</td><td>0.458</td><td>0.369</td><td>0.361</td><td>0.461</td><td>0.4</td><td>0.377</td></tr><tr><td>8 × 8</td><td>0.784</td><td>0.742</td><td>0.739</td><td>0.811</td><td>0.775</td><td>0.751</td></tr><tr><td>16 × 16</td><td>0.953</td><td>0.928</td><td>0.927</td><td>0.963</td><td>0.948</td><td>0.933</td></tr><tr><td>32 × 32</td><td>0.999</td><td>0.988</td><td>0.988</td><td>0.997</td><td>0.994</td><td>0.989</td></tr></tbody></table>

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

<figure><img src="image_14.png" /><figcaption>**Fig. 12** Pinched cylinder with end diaphragm problem normalized results</figcaption></figure>

configuration of the proposed eight-node elements. It is worth noting that the RH8 element remains slightly more sensitive to mesh distortion than the RH8s-X elements providing a lowest accuracy in the case of an irregular mesh. The convergence rate of the Mistk elements is similar to that of the RH8s-X elements.

### 7.3 Scordelis–Lo roof problem

The Scordelis–Lo roof problem was first introduced by Scordelis and Lo [70]. Since then, many researchers have used it as a standard benchmark problem to check membrane as well as bending shell elements performance. The geometry parameters of the simple curved structure illustrated in Fig. 13a are a length L = 50, a radius R = 25 and a thickness t = 0.25. The material used has a Young’s modulus E=4.32e8 and a Poisson’s ratio v = 0.0. The end diaphragms, located at each extremity of the structure, are taken into account using the boundary conditions $u_x = u_y = 0$. As was the case for the pinched cylinder with end diaphragm problem, the symmetry of the model allows one to consider only a quarter of the structure in the calculation. The vertical displacement of point A is compared to the reference solution $U_y = -0.3024$. The proposed element capability is tested using a regular as well as an irregular mesh configuration, as illustrated in Fig. 13b. The convergence of the numerical results towards the analytical solution is then presented in Table 4 and in Fig. 14.

From Table 4 it is seen that when the number of smoothing cells increases from 1 to 4 in the RH8s-X, the displacement of point A decreases progressively, approaching the reference analytical solution and the result obtained with the classical RH8 element. The same observation could be made for both regular and irregular meshes. Among the eight-node elements, the classical RH8 and the SH8 are the most accurate when the mesh is regular and their results offer a better convergence rate towards the reference solution than the results of the proposed RH8s-X elements. However, when an irregular mesh is used, the proposed RH8s-X elements converge faster towards the reference solution than the other resultant elements. Moreover, RH8s-4 remains the most accurate
<figure><img src="image_15.png" /><figcaption>**Table 4** The normalized results of Scordelis-Lo roof problem</figcaption></figure>

<table><thead><tr><th colspan="10">Regular mesh</th></tr><tr><th>Element number</th><th>SC8R</th><th>Mist1</th><th>Mist2</th><th>Mist4</th><th>SH8</th><th>RH8</th><th>RH8s-1</th><th>RH8s-2</th><th>RH8s-4</th></tr></thead><tbody><tr><td>4 × 4</td><td>1.206</td><td>1.168</td><td>1.06</td><td>0.977</td><td>1.035</td><td>1.031</td><td>1.297</td><td>1.179</td><td>1.089</td></tr><tr><td>8 × 8</td><td>1.041</td><td>1.028</td><td>1.001</td><td>0.976</td><td>1.019</td><td>1.018</td><td>1.096</td><td>1.68</td><td>1.041</td></tr><tr><td>16 × 16</td><td>1.010</td><td>1.008</td><td>1.001</td><td>0.995</td><td>1.009</td><td>1.009</td><td>1.035</td><td>1.028</td><td>1.021</td></tr><tr><td>32 × 32</td><td>1.000</td><td>—</td><td>—</td><td>—</td><td>1.000</td><td>1.000</td><td>1.012</td><td>1.010</td><td>1.008</td></tr></tbody><tbody><tr><th colspan="10">Irregular mesh</th></tr><tr><th>Element number</th><th>SC8R</th><th>SH8</th><th>RH8</th><th>RH8s-1</th><th>RH8s-2</th><th>RH8s-4</th><th></th><th></th><th></th></tr><tr><td>4 × 4</td><td>1.281</td><td>1.075</td><td>1.071</td><td>1.382</td><td>1.245</td><td>1.135</td><td></td><td></td><td></td></tr><tr><td>8 × 8</td><td>1.027</td><td>0.95</td><td>0.949</td><td>1.086</td><td>1.005</td><td>0.973</td><td></td><td></td><td></td></tr><tr><td>16 × 16</td><td>1.010</td><td>1.000</td><td>1.000</td><td>1.036</td><td>1.022</td><td>1.013</td><td></td><td></td><td></td></tr><tr><td>32 × 32</td><td>0.998</td><td>0.995</td><td>0.995</td><td>1.010</td><td>1.007</td><td>1.004</td><td></td><td></td><td></td></tr></tbody></table>

resultant element emphasizing the ability of SFEM technique to deal with membrane and bending problems without any locking. It is also noteworthy that the convergence rate of the Mistk elements is better than that of the RH8s-X elements.

## 7.4 Pinched hemispherical with 18° hole problem

This problem was initially proposed to examine the element capacity to deal with inextentional bending modes and rigid body movements. The double curvature geometry used in this problem is a severe test and a great way to estimate the Jacobian transformation accuracy. A hemispherical with a hole of 18° is subjected to two opposed concentrated loads $Fx = 1$ and $Fy = -1$, as presented in Fig. 15. Its radius is $R = 10$ and its thickness is $t = 0.04$. The Young's modulus is $E = 6.825\text{e}7$ and the Poisson's ratio is $\nu = 0.3$. Once again, given the problem symmetry, only a quarter of the geometry

was studied. The reference analytical point A displacement proposed in Macneal and Harder [71] was 9.4e-2. Results for different mesh sizes are presented in Table 5 and in Fig. 16.

From Table 5 it is seen that when the number of smoothing cells increases from 1 to 4, the displacement of point A decreases progressively from an overestimated to an underestimated solution approaching the result obtained with the classical RH8 element. The RH8s-1 element remains more accurate than the other elements as well as the latest Schwarze and Reese [23] element for refined mesh. However, the other RH8s-2 and RH8s-4 elements have a better convergence rate to the reference solution.

## 8 Concluding remarks

A review of existing developments using the SFEM is proposed in the introduction of this paper. Even though some
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14 Scordelis–Lo roof problem normalized results</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 15 Pinched hemispherical with 18° hole problem</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 16 Pinched hemispherical with 18° hole problem normalized results</figcaption></figure>

issues are still opened, it should be noted that the SFEM is already used in both shell and solid element formulations. Open issue examples are, for instance, found in Nguyen-Xuan [52] who have emphasized the necessity to introduce in their solid element formulation based on the SFEM a stabilization procedure to compensate for the induced rank deficiency and to increase its accuracy. However, this last point still needs some development. The present work has demonstrated that solid-shell elements can also successfully take advantage of this method through the resultant stress theory. Following some recent contributions, such as in Nguyen-Thanh et al. [49], a 3D variational formulation has been developed yielding in a set of resultant solid-shell elements with smoothed membrane and bending strains. In the proposed element the membrane and bending stiffness matrix integration is transferred to the boundary of cells defined on the mid-surface. No intrinsic shape function derivative and no transformation Jacobian matrix inverse calculation are required which makes the element easier to implement than a classical resultant shell element. Moreover, only the normal vectors of four cell boundaries are required per cell to integrate the membrane and bending part of the proposed element. The corresponding theoretical aspects of this method are presented. An accurate calculation of the Jacobian matrix allows the proposed element to represent very well plate and double curved structures. This element remains accurate even

<table><caption>Table 5 The normalized results of pinched hemispherical problem</caption><thead><tr><th>Element number</th><th>SC8R</th><th>Xsolid85</th><th>Sch2009</th><th>SH8</th><th>RH8</th><th>RH8s-1</th><th>RH8s-2</th><th>RH8s-4</th></tr></thead><tbody><tr><td>4</td><td>1.138</td><td>1.058</td><td>1.043</td><td>1.055</td><td>1.024</td><td>1.073</td><td>1.037</td><td>1.030</td></tr><tr><td>8</td><td>1.035</td><td>1.005</td><td>1.002</td><td>1.001</td><td>0.994</td><td>1.016</td><td>1.003</td><td>0.997</td></tr><tr><td>16</td><td>1.005</td><td>-</td><td>0.993</td><td>0.992</td><td>0.990</td><td>1.001</td><td>0.996</td><td>0.992</td></tr><tr><td>20</td><td>0.995</td><td>-</td><td>-</td><td>0.992</td><td>0.991</td><td>0.999</td><td>0.996</td><td>0.992</td></tr></tbody></table>
with highly distorted meshes and even provides the best accu-
racy among the tested resultant shell elements. Several 3D
numerical benchmark problems have been used as evaluation
tools, revealing the efficiency of the SFEM, when compared
to the classical FEM in membrane and bending problems. In
accordance with *Theorem 3* of Liu et al. [36], it is observed
that the solution obtained with SFEM elements will approach
the standard compatible displacement FEM model when the
number a smoothing cells increases. The RH8s-1 configu-
ration generally gave the closest results regarding the ana-
lytical reference solution in the tested benchmark problems.
However, this configuration could be exposed to hourglass
modes. Much like in Nguyen-Thanh et al. [49], the authors
of the present paper believe that the RH8s-2 configuration is
the most acceptable configuration. With two surface smooth-
ing cells, the present eight-node resultant solid-shell element
does not require any stabilization to remain stable and accu-
rate either for regular or distorted meshes. Only the linear
elastic aspect is considered in this contribution. However,
the authors believe that the SFEM could remain efficient
in non-linear geometric context given the work initiated in
Élie-Dit-Cosaque et al. [57]. In addition, the classical ANS
technique used in the present element to approximate the
transverse shear and transverse normal strains has been suf-
ficient enough to avoid trapezoidal as well as shear locking
in the test problems. Nevertheless, the accuracy of this tech-
nique could be significantly reduced in case of highly dis-
torted meshes and very small thicknesses. Inspired by Wu &
Wang [50], investigations on the extension of the SFEM to
the two other effects, i.e the transverse shear and through-
the-thickness effects, could be a good axis of development
for future work.

**Acknowledgments** The authors are willing to thank the following institutions for their supports: the Department of Mechanical Engineering of Laval University, the Consortium for Research and Innovation in Aerospace in Quebec (CRIAQ), the Martinique Regional Council, the European Social fund and the Aluminum Research Center (REGAL).

**References**

1. Betsch P, Gruttmann F, Stein E (1996) A 4-node finite shell element for the implementation of general hyperelastic 3D-elasticity at finite strains. Comput Methods Appl Mech Eng 130:57-79. doi:10.1016/0045-7825(95)00920-5

2. Bischoff M, Ramm E (1997) Shear deformable shell elements for large strains and rotations. Int J Numer Methods Eng 40(23):4427-4449

3. Cardoso RPR, Yoon JW (2005) One point quadrature shell element with through-thickness stretch. Comput Methods Appl Mech Eng 194(9-11):1161-1199. doi:10.1016/j.cma.2004.06.017

4. Ahmad S, Irons BM, Zienkiewicz O (1970) Analysis of thick and thin shell structures by curved finite elements. Int J Numer Methods Eng 2:419-451

5. Hughes TJR, Liu WK (1981) Nonlinear finite element analysis of shells: Part I. three-dimensional shells. Comput Methods Appl Mech Eng 26(3):331-362. doi:10.1016/0045-7825(81)90121-3

6. Hughes TJR, Liu WK (1981) Nonlinear finite element analysis of shells: Part II. two-dimensional shells. Comput Methods Appl Mech Eng 27(2):167-181. doi:10.1016/0045-7825(81)90148-1

7. Parisch H (1978) Geometrical nonlinear analysis of shells. Comput Methods Appl Mech Eng 14(2):159-178. doi:10.1016/0045-7825(78)90091-9

8. Dvorkin EN, Bathe K-J (1984) A continuum mechanics based four-node shell element for general non-linear analysis. Eng Comput 1:77-88. doi:10.1108/eb023562

9. Belytschko T, Lin JI, Chen-Shyh T (1984) Explicit algorithms for the nonlinear dynamics of shells. Comput Methods Appl Mech Eng 42(2):225-251. doi:10.1016/0045-7825(84)90026-4

10. Bathe K-J, Dvorkin EN (1985) A four-node plate bending element based on Mindlin/Reissner plate theory and a mixed interpolation. Int J Numer Methods Eng 21(2):367-383. doi:10.1002/nme.1620210213

11. Simo JC, Fox DD (1989) On a stress resultant geometrically exact shell model. Part I: formulation and optimal parametrization. Comput Methods Appl Mech Eng 72(3):267-304. doi:10.1016/0045-7825(89)90002-9

12. Simo JC, Rifai MS (1990) A class of mixed assumed strain methods and the method of incompatible modes. Int J Numer 29:1595-1638. doi:10.1002/nme

13. Belytschko T, Liu WK, Moran B (2000) Nonlinear finite elements for continua and structures, 1st edn. Wiley, New York

14. César de Sá JMA, Natal Jorge RM, Fontes Valente RA, Almeida Areias PM (2002) Development of shear locking-free shell elements using an enhanced assumed strain formulation. Int J Numer Methods Eng 53(7):1721-1750. doi:10.1002/nme.360

15. Parisch H (1995) A continuum-based shell theory for non-linear applications. Int J Numer Methods Eng 38(11):1855-1883. doi:10.1002/nme.1620381105

16. Hauptmann R, Schweizerhof K (1998) A systematic development of “solid shell” element formulations for linear and non-linear analyses employing only displacement degrees of freedom. Int J Numer Methods Eng 42:49-69

17. Hauptmann R, Schweizerhof K, Doll S (2000) Extension of the “solid shell” concept for application to large elastic and large elastoplastic deformations. Int J Numer Methods Eng 49(9):1121-1141

18. Timoshenko S, Woinowsky-Krieger S (1959) Theory of plates and shells, 2nd edn. Mcgraw-Hill College, New York

19. Klinkel S, Gruttmann F, Wagner W (2006) A robust non-linear solid shell element based on a mixed variational formulation. Comput Methods Appl Mech Eng 195(1-3):179-201. doi:10.1016/j.cma.2005.01.013

20. Kim KD, Liu GZ, Han SC (2005) A resultant 8-node solid-shell element for geometrically nonlinear analysis. Comput Mech 35(5):315-331. doi:10.1007/s00466-004-0606-9

21. Hannachi M (2007) Formulation d'éléments finis volumiques adaptés à l'analyse, linéaire et non linéaire, et à l'optimisation de coques isotropes et composites. Université de technologie de Compiègne, UTC

22. Cardoso RPR, Yoon J-W, Mahardika M, Choudhry S (2008) Enhanced assumed strain (EAS) and assumed natural strain (ANS) methods for one point quadrature solid shell elements. Int J Numer Methods Eng 75:156-187. doi:10.1002/nme.2250

23. Schwarze M, Reese S (2009) A reduced integration solid-shell finite element based on the EAS and the ANS concept-Geometrically linear problems. Int J Numer Methods Eng 80(10):1322-1355. doi:10.1002/nme.2653

24. Nguyen NH (2009) Development of solid-shell elements for large deformation simulation and springback prediction. Université de Liège, Liège
25. Andelfinger U, Ramm E (1993) EAS-elements for two-dimensional, three-dimensional, plate and shell structures and their equivalence to HR-elements. Int J Numer Methods Eng 36(8):1311-1337. doi:10.1002/nme.1620360805

26. Alves de Sousa RJ, Natal Jorge RM, Fontes Valente RA, César De Sá JMA (2003) A new volumetric and shear locking-free 3D enhanced strain element. Eng Comput 20(7):896-925. doi:10.1108/02644400310502036

27. Wriggers P, Eberlein R, Reese S (1996) A comparison of three-dimensional continuum and shell elements for finite plasticity. Int J Solids Struct 33(20-22):3309-3326

28. Alves de Sousa RJ, Cardoso RPR, Yoon J-W, Grácio JJ (2005) A new one-point quadrature enhanced assumed strain (EAS) solid-shell element with multiple integration points along thickness: Part I - geometrically linear applications. Int J Numer Methods Eng 62(7):952-977. doi:10.1002/nme.1226

29. Li LM, Peng YH, Li DY (2011) A stabilized underintegrated enhanced assumed strain solid-shell element for geometrically non-linear plate/shell analysis. Finite Elem Anal Des 47(5):511-518

30. Fontes Valente RA, Alves de Sousa RJ, Natal Jorge RM (2004) An enhanced strain 3D element for large deformation elastoplastic thin-shell applications. Comput Mech 34:38-52. doi:10.1007/s00466-004-0551-7

31. Büchter N, Ramm E, Roehl D (1994) Three-dimensional extension of non-linear shell formulation based on the enhanced assumed strain concept. Int J Numer Methods Eng 37(15):2551-2568. doi:10.1002/nme.1620371504

32. Miehe C (1998) A theoretical and computational model for isotropic elastoplastic stress analysis in shells at large strains. Comput Methods Appl Mech Eng 155(3-4):193-233. doi:10.1016/S0045-7825(97)00149-7

33. Lyly M, Stenberg R, Vihinen T (1993) A stable bilinear element for the Reissner-Mindlin plate model. Comput Methods Appl Mech Eng 110(3-4):343-357. doi:10.1016/0045-7825(93)90214-I

34. Kouhia R (2008) On stabilized finite element methods for the Reissner-Mindlin plate model. Int J Numer Methods Eng 74(8):1314-1328. doi:10.1002/nme.2211

35. Liu GR, Dai KY, Nguyen TT (2007) A smoothed finite element method for mechanics problems. Comput Mech 39(6):859-877. doi:10.1007/s00466-006-0075-4

36. Liu GR, Nguyen TT, Dai KY, Lam KY (2007) Theoretical aspects of the smoothed finite element method (SFEM). Int J Numer Methods Eng 71(8):902-930. doi:10.1002/nme.1968

37. Chen J-S, Wu C-T, Yoon S, You Y (2001) A stabilized conforming nodal integration for Galerkin mesh-free methods. Int J Numer Methods Eng 50(2):435-466

38. Simo JC, Hughes TJR (1986) On the variational foundations of assumed strain methods. J Appl Mech 53:51-54. doi:10.1115/1.3171737

39. Nguyen-Xuan H, Bordas S, Nguyen-Dang H (2008) Smooth finite element methods: convergence, accuracy and properties. Int J Numer Methods Eng 74(2):175-208. doi:10.1002/nme.2146

40. Bordas SPA, Natarajan S (2009) On the approximation in the smoothed finite element method (SFEM). Int J Numer Methods Eng 81(5):660-670. doi:10.1002/nme.2713

41. Dai KY, Liu GR (2007) Smoothed finite element method, CE006. http://hdl.handle.net/1721.1/35825

42. Dai KY, Liu GR (2007) Free and forced vibration analysis using the smoothed finite element method (SFEM). J Sound Vib 301(3-5):803-820. doi:10.1016/j.jsv.2006.10.035

43. Dai KY, Liu GR, Nguyen TT (2007) An n-sided polygonal smoothed finite element method (nSFEM) for solid mechanics. Finite Elem Anal Des 43(11-12):847-860. doi:10.1016/j.finel.2007.05.009

44. Liu GR, Nguyen-Thoi T, Nguyen-Xuan H, Lam KY (2009) A node-based smoothed finite element method (NS-FEM) for upper bound

solutions to solid mechanics problems. Comput Struct 87(1-2):14-26. doi:10.1016/j.compstruc.2008.09.003

45. Nguyen-Thoi T, Liu GR, Nguyen-Xuan H (2011) An n-sided polygonal edge-based smoothed finite element method (nES-FEM) for solid mechanics. Int J Numer Methods Biomed Eng 27(9):1446-1472. doi:10.1002/cnm.1375

46. Nguyen-Thanh T, Liu GR, Dai KY, Lam KY (2007) Selective smoothed finite element method. Tsinghua Sci Technol 12(5):497-508. doi:10.1016/S1007-0214(07)70125-6

47. Hughes TJR (1980) Generalization of selective integration procedures to anisotropic and nonlinear media. Int J Numer Methods Eng 15(9):1413-1418. doi:10.1002/nme.1620150914

48. Nguyen-Xuan H, Rabczuk T, Bordas S, Debongnie J (2008) A smoothed finite element method for plate analysis. Comput Methods Appl Mech Eng 197(13-16):1184-1203. doi:10.1016/j.cma.2007.10.008

49. Nguyen-Thanh N, Rabczuk T, Nguyen-Xuan H, Bordas S (2008) A smoothed finite element method for shell analysis. Comput Methods Appl Mech Eng 198(2):165-177. doi:10.1016/j.cma.2008.05.029

50. Wu CT, Wang HP (2013) An enhanced cell-based smoothed finite element method for the analysis of Reissner-Mindlin plate bending problems involving distorted mesh. Int J Numer Methods Eng 95(4):288-312. doi:10.1002/nme.4506

51. Liu GR, Nguyen-Thoi T, Lam KY (2009) An edge-based smoothed finite element method (ES-FEM) for static, free and forced vibration analyses of solids. J Sound Vib 320(4-5):1100-1130. doi:10.1016/j.jsv.2008.08.027

52. Nguyen-Xuan H (2008) A strain smoothing method in finite elements for structural analysis. Université de Liège, Liège

53. Bordas SPA, Rabczuk T, Hung N-X, Nguyen VP, Natarajan S, Bog T, Quan DM, Hiep NV (2010) Strain smoothing in FEM and XFEM. Comput Struct 88(23-24):1419-1443. doi:10.1016/j.compstruc.2008.07.006

54. Puso MA, Solberg J (2006) A stabilized nodally integrated tetrahedral. Int J Numer Methods Eng 67(6):841-867. doi:10.1002/nme.1651

55. Cui X, Liu G, Li GY, Zhao X, Nguyen TT, Sun GY (2008) A smoothed finite element method (SFEM) for linear and geometrically nonlinear analysis of plates and shells. Comput Model Eng Sci 28(2):109-125. doi:10.3970/cmes.2008.028.109

56. Liu SJ, Wang H, Zhang H (2010) Smoothed finite elements large deformation analysis. Int J Comput Methods 7(3):513-524. doi:10.1142/S0219876210002246

57. Élie-Dit-Cosaque X, Gakwaya A, Lévesque J, Guillot M (2011) Smoothed finite element method for the resultant eight-node solid shell element analysis. In: Simulia customer conference, Barcelona

58. Rashid MM, Sadri A (2012) The partitioned element method in computational solid mechanics. Comput Methods Appl Mech Eng 237-240:152-165. doi:10.1016/j.cma.2012.05.014

59. Simo JC, Armero F, Taylor RL (1993) Improved versions of assumed enhanced strain tri-linear elements for 3D finite deformation problems. Comput Methods Appl Mech Eng 110(3-4):359-386. doi:10.1016/0045-7825(93)90215-J

60. Vu-Quoc L, Tan XG (2003) Optimal solid shells for nonlinear analyses of multilayer composites. I. Statics. Comput Methods Appl Mech Eng 192(9-10):975-1016. doi:10.1016/S0045-7825(02)00435-8

61. Domissy E (1997) Formulation et évaluation d'éléments finis volumiques modifiés pour l'analyse linéaire et non linéaire des coques. Université de technologie de Compiègne, UTC

62. Sze KY, Yao L (2000) A hybrid-stress ANS solid-shell element and its generalization for smart structure modelling—part I: solid-shell element formulation. Int J Numer Methods Eng 48:545-564

63. Hughes TJR, Tezduyar TE (1981) Finite elements based upon mindlin plate theory with particular reference to the four-node bilin-
ear isoparametric element. J Appl Mech 48(3):587-596. doi:10.1115/1.3157679

64. Batoz J-L, Dhatt G (1992) Modélisation des structures par éléments finis: coques, Les Presses de l'Université Laval

65. Oden JT, Reddy JN (1974) On dual-complementary variational principles in mathematical physics. Int J Eng Sci 12:1-29. doi:10.1016/0020-7225(74)90073-1

66. Hodge PG (1970) Continuum mechanics: an introductory text for engineers. McGraw-Hill, New York

67. Hao S, Liu WK, Belytschko T (2004) Moving particle finite element method with global smoothness. Int J Numer Methods Eng 59(7):1007-1020. doi:10.1002/nme.999

68. Zienkiewicz OC, Taylor RL, Too JM (1971) Reduced integration technique in general analysis of plates and shells. Int J Numer Methods Eng 3(2):275-290. doi:10.1002/nme.1620030211

69. Abaqus user's manual, ver 6.8, (2007) Abaqus user's manual. Dassault Systèmes, Simulia

70. Scordelis AC, Lo KS (1964) Computer analysis of cylindrical shells. J Am Concr Inst 69:539-559

71. Macneal RH, Harder RL (1985) A proposed standard set of problems to test finite element accuracy. Finite Elem Anal Des 1:3-20. doi:10.1016/0168-874X(85)90003-4