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

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.


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][5][6][7][8][9][10][11][12][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 inplane 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 meshfree 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 situations 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 variationally 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 formulation 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 technique, 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 conventional 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 corresponding 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 integration 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 overcome shear locking, consists in using different smoothing cell numbers for each stiffness matrix contribution. However, 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, various 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 components. Their proposed element exhibited good accuracy and time efficiency compared with equivalent standard elements. Nguyen-Thanh et al. [49] then extended this work to quadrilateral flat shell elements and emphasized the high potential capabilities of the SFEM. Unlike the MITC4 element, which uses an isoparametric mapping, their proposed element calculates the membrane and bending parts of the stiffness 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 element. 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 proposed 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 stabilization 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 element and the four or eight subcell element. Due to the relative novelty of the SFEM, its non-linear development and associated applications still remain limited. Cui et al. [55] developed 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.
, , ξ ξ ξ X 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 eightnode 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 X (X, Y, Z ), a convected coordinate system ξ ξ 1 , ξ 2 , ξ 3 and a local Cartesian coordinate system 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].

Fig. 3 Resultant eight-node solid-shell element geometry and midsurface concept
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: which, after introducing the mid-surface concepts and considering the parameter ξ 3 as the thickness direction (see Fig. 3), can also be written as: where X 0,L and X 0,U are respectively the position vector of a point on the lower and upper surface in the global Carte- is the position vector of the mid-surface and X 0 = X 0,U −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 X (X, Y, Z ) can similarly be represented as: By introducing the displacements field vector U q , the undeformed and the deformed configurations are related by the standard mapping at point q, as: in the global Cartesian coordinate system 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: where U and U are defined as U = U U +U L 2 and U = , with U L and 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: 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: The local Cartesian base vectors are defined in the undeformed configuration from the covariant base vectors as: Similarly, the local Cartesian base vectors are defined in the deformed configuration from the covariant base vectors as: 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 1 = r x r y r z such that the local components of a position or displacement vector u q = u v w of an arbitrary point q are given by: For sake of simplify, the subscript ( ) q is omitted in the subsequent developments.

Kinematics of shell deformation
The Green-Lagrange strain tensor components ε ξ kl , oriented in the natural coordinate system ξ ξ 1 , ξ 2 , ξ 3 are (see Vu-Quoc and Tan [60] for more detail): Considering only the linear terms, the above Eq. (18) becomes: Considering Eq.
The Green-Lagrange strain tensor ε, can be formulated in the directions of the local Cartesian base vectors as: It can also be formulated in the direction of the contravariant base vectors g 1 , g 2 , g 3 defined in Eq. (9) as: where X i , i = 1, 2, 3 are the global Cartesian coordinates defined in Table 1.
Identifying the terms of Eqs. (21) and (22), the transformation between the local Cartesian and the natural components of the Green-Lagrange strain tensor represented respectively by ε i j and ε ξ kl is: 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: Limiting the following development to the linear strain components, as for the natural expression of the Green-Lagrange strain tensor, gives the expression: Following the concepts introduced in Eq. (5), the displacements vector can also be written as: as the displacement mean value on the mid- 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: 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: 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: 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 with, in case of small displacements and small strains formulation, (30) expressed on the mid-surface ξ 3 = 0 simplifies to T ξ 3 =0 : However, the calculations required for evaluating the terms of the inverse J −1 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 ε mb , transverse shear γ and transverse normal ε zz , the linear strain vector can be written, in line with the work proposed in Kim et al. [20], in this compact form: In Eq. (32), ε mb is the combination of the membrane and the bending strains components, respectively ε m and ε b such that with the membrane strain operator defined as: and where: with the bending strain operator defined as: From Eq. (28), the matrix expressions of the transverse shear strain γ and the transverse normal strain ε zz are also obtained.
The transverse shear strains γ = γ yz γ xz can be written in matrix form as: where B γ is: The transverse normal strain ε zz can be written in matrix form as: where B zz is: Considering Eqs. (33), (35), (37) and (39), the strain tensor can be written: where B is a strain operator including membrane-bending, transverse shear and transverse normal components (See Kim et al. [20]).

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: where N I , I = 1, . . . , 4 are the standard bilinear shape functions of a quadrilateral element, 3) are the position vectors of the element midsurface corners and 3) are the orientation vectors at the element midsurface corners. The displacements field expressed in Eq. (5) can be approximated in finite element form as u ≈ u e . For sake of simplicity the superscript ( ) e is omitted in this last expression, such that the nodal displacements field u is given by: where 3) are the displacement vectors of the element mid-surface corners and 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: where u I = u I v I w I u I v I w I T for each corner and The detailed smoothing procedure to calculate this membrane and bending element strain fields, respectively represented by ε m and ε b , as well as the corresponding strain operators B m I and B b I in their smoothed form are given in Sect. 6.1.
Developing Eq. (37) in the element mid-surface, at ξ 3 = 0, the element transverse shear components can be expressed as: where the natural transverse shear strain γ ξ in matrix form is obtain by developing Eq. (20) as: with the transverse shear strain operator B γ, 0 I , expressed in the element mid-surface ξ 3 = 0 by: and T γ ξ 3 =0 is a sub-transformation matrix extracted from Eq. (31) as: Using Eq. (39), the element normal strain finite element approximation evaluated on the element mid-surface ξ 3 = 0 then becomes: where the natural normal strain in matrix form is obtained by developing Eq. (20) as in Eq. (53): where the normal strain operator B ζ ζ I , evaluated on the element mid-surface ξ 3 = 0 is expressed as:

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  A 2 , B 1 , B 2 ) located in the middle of the midsurface 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: Equations (55) to (58) are then used to obtain an improved interpolated shear strain using linear interpolation functions: where −1 ≤ ξ 1 , ξ 2 ≤ 1.

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 throughthe-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- 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: 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: where −1 ≤ ξ 1 , ξ 2 ≤ 1.

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: where the internal and the external work expressions are given by: The weak form expression of the virtual work principle can be written as detailed in Hodge [66]: with the following definition: where the stress field σ and the constitutive parameter H are assumed to be related by the following constitutive equation: where, for an isotropic elastic material in which and the modified Lamé's coefficients are λ = E.ν (1−ν 2 ) and μ = E 2.(1+ν) with E the Young's modulus and ν the Poisson's ratio.
The constitutive parameter 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: where σ m = H mb ε m and σ b = H mb · ε 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 internal virtual work δW int , defined in Eq. (70), into membranebending, transverse shear and normal components as: Introducing Eq. (77) in Eq. (78) gives: Then, using Eq. (79) and separating integration into thickness and mid-surface integrations gives: where h is the shell thickness and AH , DH , AH and E H are through-the-thickness integrated material parameters thus defining the new resultant constitutive law expressed as: Introducing the finite element approximation, the element internal virtual work δW e int in Eq. (78) can be expressed as: where where the element stiffness matrix k is: Smoothed finite element method

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: where ϕ is a smoothing function with ϕ (x − x c ) = 1/S c x ∈ c 0 x / ∈ c and S c = c d is the area of the smoothing cell c with 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 weakform 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: By applying the divergence theorem, C u · i n j d , the smoothed strain formulation is then expressed as an integration of the strain matrix on the cell boundary and Eq. (93) becomes: Hence the associated smoothed Green-Lagrange strainε becomes: In the previous equation, 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 c , the normal vectors of its boundary 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: Introducing the local displacements field u ξ 1 , ξ 2 , ξ 3 = u ξ 1 , ξ 2 + ξ 3 · u ξ 1 , ξ 2 , the deformation gradient can then be written as: 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: In Eqs. (101) to (103), I is the mid-surface node indices from 1 to 4; x c is the unique gauss point along the element edge and l c b is the length of c b ; n x and n y are the components of the cell border outward normal vector n; S c is the area of cell c; n c is the total number of cells; N I 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 I , l c b ,n and S c 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): The approximation of the shear strain field presented in section 4.2 is used to determine the shear strain matrix B AN S,γ . 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 midsurface. 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: Then, the approximation of the transverse normal strain field presented in Sect. 4.3 is used to determine the transverse normal strain matrix B AN S,zz . It is this latter which is used in the numerical integration of the through-the-thickness normal stiffness Eq. (65) defined as: The adopted integration scheme is 2 × 2.
The total local element stiffness matrix then becomes: The local to global stiffness matrix transformation is similar to the method used in Kim et al. [20].

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 m and B b , are computed by multiplying the intrinsic derivative with the inverse of the jacobian matrix J −1 . Then the stiffness matrix K mb 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 c b , 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 1 2 , l 1 3 ,l 1 4 and normals n 1 1 , n 1 2 , n 1 3 , n 1 4 , 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 c b 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 −1 to determine the smoothed operatorB mb as presented in Eqs. (102) and (103). Then the smoothed stiffness matrix 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.

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): 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ε, the orthogonal condition takes an analogous form to that in Liu et al. [36]:  [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 ξ 1 = ξ 2 = ξ 3 = 0. The strain operator B is then equivalent to an average value over the isoparametric element mid-surface. Similarly, the smoothed strain operatorB is an average value over the element midsurface. 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.

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. 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] : fournode shell element with smoothed membrane and bending strains.
All results have been normalized to the exact reference solution of each problem.

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 ν = 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.

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 ν = 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   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.

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 ν = 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  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.

Pinched hemispherical with 18 o 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 o 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.825e7 and the Poisson's ratio is ν = 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.

Concluding remarks
A review of existing developments using the SFEM is proposed in the introduction of this paper. Even though some     issues are still opened, it should be noted that the SFEM is already used in both shell and solid element formulations.

Number of elements
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 with highly distorted meshes and even provides the best accuracy 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 configuration generally gave the closest results regarding the analytical 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 smoothing cells, the present eight-node resultant solid-shell element does not require any stabilization to remain stable and accurate 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 sufficient enough to avoid trapezoidal as well as shear locking in the test problems. Nevertheless, the accuracy of this technique could be significantly reduced in case of highly distorted 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 throughthe-thickness effects, could be a good axis of development for future work.