# A nine nodes solid-shell finite element with enhanced pinching stress

Mouhamadou Dia, Nahiene Hamila, Mickaël Abbas, Anthony Gravouil

► To cite this version:

Mouhamadou Dia, Nahiene Hamila, Mickaël Abbas, Anthony Gravouil. A nine nodes solid-shell finite element with enhanced pinching stress. Computational Mechanics, 2020, 10.1007/s00466-020-01825-1. [hal-02483993](https://hal.science/hal-02483993)

HAL Id: hal-02483993

https://hal.science/hal-02483993

Submitted on 23 Aug 2021

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
# A nine nodes solid-shell finite element with enhanced pinching stress

Mouhamadou Dia<fnref n="1" />,<fnref n="2" /> · Nahiene Hamila<fnref n="2" /> ID · Mickaël Abbas<fnref n="1" /> · Anthony Gravouil<fnref n="2" />

## Abstract

In this paper we present a low-order solid-shell element formulation—having only displacement degrees of freedom (DOFs), i.e., without rotational DOFs. The element has an additional middle node, that allows efficient and accurate analyses of shell structures using elements at extremely high aspect ratio. The formulation is based on the Hu–Washizu variational principle leading to a novel enhancing strain and stress tensor that renders the computation particularly efficient, with improved in-plane and out-of-plane bending behavior (Poisson thickness locking). The middle-node is endowed with only one degree of freedom, in the thickness direction, allowing the assumption of a quadratic interpolation of the transverse displacement. Unlike solid-shell finite elements reported previously in the literature and formulated under the hypothesis of plane stress or with enhanced assumed strain parameter, the new solid-shell element here mentioned uses a complete three-dimensional constitutive law and gives an enhanced pinching stress, thanks to the middle-node. Moreover, to handle the various locking problems that usually arise on solid-shell formulation, the reduced integration technique is used as well as the assumed shear strain method. Finally to assess the effectiveness and performance of this new formulation, a set of popular benchmark problems, involving geometric non-linear analysis as well as elastic-plastic behavior has been investigated.

**Keywords** Solid-Shell · 3D Constitutive law · Improved normal stress · Robust stabilization

## 1 Introduction

Shell-like structures are largely present in most engineering design and process control. To model such structures, engineers generally use classical shell element based on the degenerated shell concept or classical shell theories. These elements perform well for the simulation of bending problems; in linear as well as for some nonlinear problems. However in certain engineering problems, like sheet metal forming in case of hemming, bottoming, ironing, hydro-

forming or simply V-type bending, those elements can show some inadequacy related to the hypothesis they are embedded with. The main limiting hypothesis being the plane stress state, which literally means that the normal transverse stress is negligible. Under such assumption, one must modify the 3D nonlinear constitutive law to be numerically integrated. This is very often a laborious work particularly when the material law is intrinsically three dimensional.

There are several approaches used to overcome the problem of plane stress state in shell element and they can be divided into two main groups or methodologies. The first group is based on hexahedral elements, the so-called ‘solid-shell’, where modifications are introduced to enforce flexural behavior that permits stretching through a given thickness direction (see, for instance, References [1–14] to name just a few contributions). The second group, which could be called ‘shell-solid’, is based mainly upon Mindlin–Reissner shell elements with three displacements and two local rotational degrees of freedom at each node while the three-dimensional (3D) constitutive behavior is enforced via additional degrees of freedom giving the so-called 5, 6 or 7-parameters shell models, see for example [15–29], among many others. The 5-parameter shell models have been enriched by a desired number of parameters at nodes to permit a representation

<footnote>✉ Nahiene Hamila
Nahiene.hamila@insa-lyon.fr</footnote>

<footnote>Mouhamadou Dia
mouhamadou-m.dia@edf.fr; mouhamadou.dia@insa-lyon.fr</footnote>

<footnote>Mickaël Abbas
mickael.abbas@edf.fr</footnote>

<footnote>Anthony Gravouil
anthony.gravouil@insa-lyon.fr</footnote>

<footnote n="1">¹ EDF R&D ERMES and IMSIA, UMR EDF/CNRS/CEA/ENSTA 9219, 7 boulevard Gaspard Monge, 91120 Palaiseau, France</footnote>

<footnote n="2">² Univ Lyon, INSA-Lyon, CNRS UMR 5259, LaMCoS, 69621 Lyon, France</footnote>
of through-thickness stretching. Several shell elements have been developed which explicitly account for the thickness change as an additional degree of freedom leading to 6-parameter models. But due to the coupling with the Poisson ratio in bending dominated cases, the linear displacement field in thickness direction gives a constant strain which in turn causes an artificial normal stress. As a remedy, either the coupling terms have to be removed from the constitutive law by a plane stress state or the shell formulation has to be extended by a linear normal strain, leading to a 7-parameter model. This is achieved by a quadratic variation of the normal displacement field and the details can be found in References [15–18,20–22,24].

In a recently published paper [30], Sansalone et al. have proposed a third approach where an additional node is introduced in the center of three-node and four-node shell elements with only two through-thickness translational degrees of freedom of the upper and lower surfaces of the shell. Then a full 3D constitutive strain-stress behavior can be used. For triangles in bending state, either based on Kirchhoff's or on Mindlin's assumptions, it has been shown that the results are exactly the same as those given by the initial formulation of these elements using a plane stress hypothesis. For quadrilaterals, the results are slightly different but many numerical examples including nonlinear computations prove that those differences are not significant.

Therefore, the aim of this paper is to extend this methodology by using an additional node at the center of an eight-node brick element in a specific solid-shell formulation. Bassa et al. [31] have proposed a similar approach. Nevertheless the studies have shown a low convergence performance of that formulation in an implicit code, specially in nonlinear problems, due to the three-parameters stabilization procedure. The shear strain adopted in the $SB_{9\gamma}25$ was interpolated following the work of Dvorking and Bathe [32–34]. The Reissner multiplicative function was used as interpolation through the thickness direction to ensure a static admissibility condition. However such interpolation is only valid in linear situations. It is then essential to realize that in this paper, the third approach proposed by Sansalone [30] is here applied to a brick element. An additional node is then introduced with only one through-thickness translational dof with a robust stabilisation procedure that performs well in implicit code. This new method goes beyond the EAS approach. In fact the EAS parameter is often added to solve Poisson thickness locking and volumetric locking by transforming a linear normal strain into a quadratic normal strain. However the value of this parameter is so low that it does not improve the pinching stress that much. In a problem with a shell structure under pressure in one or both side, the solid-shell element with EAS parameter alone gives a normal stress that is almost constant, and do not verify the Neumann boundary conditions. The aim of the solid-shell here proposed is to overcome that limitation

by enhancing the normal stress thank to a additional node. Moreover, only reduced integration is practically interesting for this element, for further non-linear computations. Therefore, all the results given in this paper correspond to this specificity.

Since this approach differs significantly from the previous two ones, only the main features of the solid-shell concept are recalled. The so called solid-shell elements are very attractive and are still the subject of many researches. As standard elements, solid-shell elements can incorporate the normal stress along thickness direction if the pinching is correctly handled. General 3D-constitutive material can be used without any assumption on the normal stress. These elements have only translational degrees of freedom (DOF) which make their formulation very simple. The difficulties associated with complex shell formulation with nodal rotation are then avoided. Moreover, since they have only translational DOF, solid-shell elements can easily be connected with standard solids elements when there is coexistence of three-dimensional and structural zones. All these advantages have lead to develop a set of solid shell element in the last decades: [3,7,10,24,35–41].

Now, for a solid-shell to perform accurately, it must avoid the amount of locking phenomena that low-order elements encounter when modeling shell structures. If not handled conveniently, solid-shell elements suffer from locking pathologies that lead the elements to give poor results, specially in case of out-of-plane bending load analysis or isochoric plasticity.

The first locking pathology that low order solid-shell element encounter is volumetric locking which arises when the material is incompressible or nearly-incompressible, as for rubber like material (see the work of Reese and Wriggers [42]) or for metal plasticity (see Simo and Taylor [43], Miehe [44] to name just a few). In such material the non-vanishing volumetric strain makes the element stiffer, resulting in excessively small deformations. Several solutions have been proposed to solve volumetric locking. There are, among others, The reduced or selective integration technique proposed by Zienkiewicz [45], the B-bar approach of Hughes [46], applied for solid elements in the work of Belytschko [47] and the enhanced assumed strain initially proposed by Simo et al. [15,48,49]. These methods are largely applied in many solid-shell formulations. See for example the work of Mas sud et al. [14], Alves de Souza [37,50,51], Cardoso et al [12,24], Sze and Yao [5], Klinkel et al. [38]. It's worth noting that the reduced integration or selective integration technique are not cost-free. Indeed appropriate hourglass stabilization techniques are necessary to prevent the spurious deformations modes that may arise and ensure a correct rank of the element stiffness matrix. It's well known that the spurious patterns or hourglass modes correspond to the kernel vectors of the stiffness matrix, aside from the rigid body modes. Dif-
ferent techniques have been proposed to deal with the rank deficiency due to selective or reduced integration technique. One can see the remarkable work of Belytschko et al. [52–54], Reese et al. [10,36] as well as Combescure et al. [7,41].

Poisson thickness locking is due to resulting incorrect-constant distribution of the normal stress in the thickness direction. To avoid this locking phenomenon it is necessary to assure that the normal stress varies linearly along the thickness direction in bending situation. It can be done using a quadratic interpolation or the compatible normal displacement e.g. Parisch [1]. Another and most often used method is to enrich the transverse normal strain by mean of the enhanced-assumed stress method [10,37,48].

A further locking effect observed for solid-shell element is the phenomenon of a so called curvature locking or sometime “trapezoidal locking”. This phenomenon is found in structures where the out of plane element edges are not perpendicular to the mid-layer, which is the case for originally curved or heavily deformed structures. This locking effect can be overcomed by using the naturally assumed interpolation of the normal strain as proposed by Bischoff and Ramm [18] or Betsch and Stein [17].

Another commonly encountered locking pathology in solid-shell is the transverse shear locking. This happens because normal strains of linear element are coupled by shear strain. Low order elements do not have pure bending modes to behave correctly for pure bending load cases. Hence, there are parasitic shear strains that appear and become more important than normal strain. This makes the element stiffer than necessary resulting in a poor bending behavior. An interesting remedy is the Assumed Natural Strain (ANS) method first proposed by Hughes and Tezduyar[55], followed by Wempner [56] then by Dvorkin et al. [34,57,58] for shell elements. The approach has since then been applied in solid-shell formulation by many authors. See e.g. [3,12,13,24,38,59].

In this paper the third approach proposed by Sansalone et al. [30] to eliminate plane stress-state in classical shell element is extended to a brick element. A nine nodes solid-shell element is then proposed. The first eight nodes are those of a classical hexahedral element and have three translation degrees of freedom each. The ninth node is an additional node that have only one translation degree of freedom, in the element-thickness direction. The latter node is added to create a quadratic interpolation of the transverse displacement and consequently enrich the element pinching strain. This way the element can accurately work with a three dimensional constitutive law without the common plane stress hypothesis. Poisson thickness locking is then naturally avoided. But not only that, the additional DOF allows to enrich the element pinching stress and further match the upper and lower boundaries condition of the shell. We adopt the ANS technique with 4 tying points for the shear strain as well as the normal transverse strain. Since the reduced

integration technique is adopted we use the stabilization technique following the work of Belytschko [54] and Schwartz [13]. This new element has a wide range of applications, showing very good convergence, robustness and accuracy in nonlinear problems. The element is implemented into the quasi-static implicit software code_aster [60] developed by the French energy and electricity company (eDF).

## 2 Variational formulation

The ninth node of this solid-shell aims to enhance the element displacement gradient and further improves volumetric and Poisson thickness locking. The Hu–Washizu principle is the starting point for this element formulation. As being proposed by Simo and Rifai [61], the displacement field $\mathbf{u} \in \mathcal{V}$, the conjugated Green–Lagrange strain tensor $\mathbf{E} \in \mathcal{E}$ and the second Piola–Kirchhoff stress $\mathbf{S} \in \mathcal{S}$ are treated as independent variables. $\mathcal{V}$, $\mathcal{E}$ and $\mathcal{S}$ being respectively the space of admissible displacement, strain and stress. The Hu–Washizu principle is written as follows:

$$ \pi^{HW}(\mathbf{u}, \mathbf{E}, \mathbf{S}) = \int_{\Omega_0} W(\mathbf{E}) d\Omega_0 + \int_{\Omega_0} \mathbf{S} : \left[ \frac{1}{2}(\mathbf{F}^T\mathbf{F} - \mathbf{I}) - \mathbf{E} \right] d\Omega_0 - \pi_{\text{ext}}(\mathbf{u}) \quad (1) $$

where $W$ is the strain energy, $\mathbf{F}$ the deformation gradient depending on the displacement field $\mathbf{u}$, $\mathbf{I}$ the metric tensor and $\pi^{ext}(\mathbf{u})$ the external force power. As being advocated first by Andelfinger and Ramm [62] and further by Bischoff and Ramm [18] for large deformation with (EAS) method, the approach proposed herein is also based on an enrichment of the Green–Lagrange strain tensor.

$$ \mathbf{E} = \mathbf{E}^H + \mathbf{E}^9 \quad \text{with} \quad \mathbf{E}^H = \frac{1}{2} (\mathbf{F}^T\mathbf{F} - \mathbf{I}) \quad (2) $$

This means that the eight vertices nodes displacement-dependent strain tensor $\mathbf{E}^H$ of the element is enriched by an additional enhanced assumed strain $\mathbf{E}^9 \in \mathcal{E}^9$ thanks to the ninth node. The variation of the above functional is obtained from the directional derivative and leads to

$$ \delta\pi^{HW}(\mathbf{u}, \mathbf{E}^9, \mathbf{S}) = \int_{\Omega_0} \left[ \frac{\partial W}{\partial \mathbf{E}} : \delta\mathbf{E}^H + \frac{\partial W}{\partial \mathbf{E}} : \delta\mathbf{E}^9 \right] d\Omega_0 - \delta\pi_{\text{ext}} - \int_{\Omega_0} \mathbf{S} : \delta\mathbf{E}^9 d\Omega_0 - \int_{\Omega_0} \delta\mathbf{S} : \mathbf{E}^9 d\Omega_0 \quad (3) $$

As in many researches, the three-field functional is reduced to a two field functional as suggested by Simo and Rifai [61], by choosing the interpolation such that $\mathbf{S}$ and $\mathbf{E}^9$ become orthogonal. A slightly different approach is applied in this new formulation. If the element is not pinched, in a
sense that there is no pressure applied above or below the element, the formulation goes with the orthogonality condition. Otherwise, if the element is pinched, let's say by an upper pressure $P_u$ and a lower pressure $P_l$, the second Piola Kirchoff stress $S$ is chosen such that the corresponding Cauchy stress be linearly dependent to the applied pressure:

$$ \sigma_{33} = -\frac{1}{2}(1-\zeta)P_l - \frac{1}{2}(1+\zeta)P_u \quad (4) $$

$\zeta$ being the element thickness parameter, in the covariant frame. This way, the element normal stress is statically enhanced so that a correct pinching stress is derived.

Assuming a nil body force, the Euler–Lagrange equations [18] associated with Eq. (3) are the standard equilibrium equation of the domain:

$$ \begin{aligned} \operatorname{div} \left( \mathbf{F} \frac{\partial W}{\partial \mathbf{E}} \right) &= 0 \\ \mathbf{S} - \frac{\partial W}{\partial \mathbf{E}} &= 0 \\ \mathbf{E}^9 &= 0 \end{aligned} \quad (5) $$

Although $\mathbf{E}^9 = 0$ for the continuum problem, in general $\mathbf{E}_h^9 \neq 0$, when we introduce finite element approximations. By denoting $\mathcal{D}$ the domain occupied by the body, we can observe that the space of enhanced strain field is in $[L_2(\mathcal{D})]^6$ (see Simo [61]). Hence, no inter-element continuity on the $\mathbf{E}^9$ need be enforced when construction finite element approximation. Note also that the enhanced strain interpolation $\mathcal{E}_h^9$ and the standard strain interpolation defined by $\nabla^S \psi_h$ are independent in the sense that: $\mathcal{E}_h^9 \cap \nabla^S \psi_h = \emptyset$

## 3 Kinematics

The SB9 element has nine nodes, one in-plane integration point and the ability to accommodate several integration points along the thickness direction of the element. Figure 1 represents the element topology with the order of node numbering (related to the isoparametric coordinate system of the element, defined by the natural coordinates $\xi, \eta, \zeta$). Also, in the same figure, the distribution of Gauss-Lobatto integration points is given. Compared with the other eight-node ‘solid-shell’ bricks, the presence of a supplementary node has two main aims. First getting a linear normal strain component which, along with a full three-dimensional constitutive strain-stress behavior, allows to achieve similar results in bending cases as those obtained with the usual plane stress state hypothesis. For that, the ninth node DOF plays the role of an extra parameter essential for a quadratic interpolation of the displacement in the thickness direction. The second advantage is that this DOF has a physical meaning and, for

<figure><img src="image_2.png" /><figcaption>Fig. 1 Location of the Gauss-Lobatto integration points for the SB9 element</figcaption></figure>

instance, a strength equivalent to a normal pressure can be prescribed to improve the normal stress when the shell structure is moderately thick. In this section we focus on the interpolation of the 8 vertices nodes, we will give more details about the ninth node in the Sect. 6.

The location of nodes in the isoparametric coordinate system is given by the following isoparametric vectors:

$$ \begin{aligned} a_1^I &= \{-1, 1, 1, -1, -1, 1, 1, -1\} \\ a_2^I &= \{-1, -1, 1, 1, -1, -1, 1, 1\} \\ a_3^I &= \{-1, -1, -1, -1, 1, 1, 1, 1\} \end{aligned} \quad (6) $$

where $I$ means node number. Auxiliary vectors are defined from previous combinations of the isoparametric vectors:

$$ \begin{aligned} h_1^I &= \{-1, -1, 1, -1, 1, -1, 1, -1\} \\ h_2^I &= \{1, 1, -1, -1, -1, -1, 1, 1\} \\ h_3^I &= \{1, -1, -1, 1, -1, 1, 1, -1\} \\ h_4^I &= \{-1, 1, -1, 1, 1, -1, 1, -1\} \end{aligned} \quad (7) $$

These vectors form the basis for the interpolation field of the nodal displacements. The shape functions $N_I(\xi, \eta, \zeta)$ are obtained from a linear combination of the above-mentioned isoparametric vectors:

$$ N_I = \frac{1}{8} (1 + \xi a_1^I + \eta a_2^I + \zeta a_3^I + \xi \eta h_1^I + \eta \zeta h_2^I + \xi \zeta h_3^I + \xi \eta \zeta h_4^I) \quad (8) $$

The derivatives of the shape functions with respect to the isoparametric coordinates are also spanned by the isoparametric vectors in the form:
$$
\begin{align}
N_{I,\xi} &= \frac{1}{8}(a_1^I + h_2^I\zeta + h_3^I\eta + h_4^I\eta\zeta) \nonumber \\
N_{I,\eta} &= \frac{1}{8}(a_2^I + h_1^I\zeta + h_3^I\xi + h_4^I\xi\zeta) \tag{9} \\
N_{I,\zeta} &= \frac{1}{8}(a_3^I + h_1^I\eta + h_2^I\xi + h_4^I\xi\eta) \nonumber
\end{align}
$$

One of the advantages of the SB9 element is the use of only
translational degrees of freedom. Therefore the position of a
point inside the element in the initial (X) or the current (x)
state is obtained from the interpolation of the nodal transla-
tional degrees of freedom:

$$
\mathbf{X}(\xi, \eta, \zeta) = N_I(\xi, \eta, \zeta)\mathbf{X}_I \quad (10)
$$

$$
\mathbf{x}(\xi, \eta, \zeta) = N_I(\xi, \eta, \zeta)\mathbf{x}_I \quad (11)
$$

where **x**<sub>I</sub> and **X**<sub>I</sub> are the nodal coordinates in the current and initial positions. The displacement field is also obtained from nodal values after a proper interpolation with shape functions

$$
\mathbf{u}(\xi, \eta, \zeta) = N_I(\xi, \eta, \zeta)\mathbf{u}_I \quad (12)
$$

where **u**<sub>I</sub> are the nodal displacement. Keep in mind that the
effect of the ninth node will be discussed in a later sec-
tion. The derivatives of displacements in the isoparametric
coordinate system are interpolated from the shape functions
derivatives as

$$
\begin{align}
\mathbf{u}_{,\xi}(\xi, \eta, \zeta) &= N_{I,\xi}(\xi, \eta, \zeta)\mathbf{u}_I \notag \\
\mathbf{u}_{,\eta}(\xi, \eta, \zeta) &= N_{I,\eta}(\xi, \eta, \zeta)\mathbf{u}_I \tag{13} \\
\mathbf{u}_{,\zeta}(\xi, \eta, \zeta) &= N_{I,\zeta}(\xi, \eta, \zeta)\mathbf{u}_I \notag
\end{align}
$$

$$
\mathbf{u}_{,\zeta}(\xi, \eta, \zeta) = N_{I,\zeta}(\xi, \eta, \zeta)\mathbf{u}_I
$$

4 Strain field

The deformation gradient is written as follows:

$$
\mathbf{F} = \frac{\partial \mathbf{x}}{\partial \mathbf{X}} = \frac{\partial x^i}{\partial \xi^k} \mathbf{e}_i \otimes \frac{\partial \xi^k}{\partial X^j} \mathbf{e}_j = \mathbf{g}_i \otimes \mathbf{G}^i \quad (15)
$$

F is a tensor which maps the reference basis **G**<sub>i</sub> to the current one **g**<sub>i</sub> (Fig. 2). **e**<sub>i</sub> is the Cartesian base vector. For simplicity, the superscript **H** will be omitted and **E**<sup>**H**</sup> will be written as **E**. Further the Green Lagrange strain tensor is represented by it's Cartesian and covariant components as follows

$$
\begin{align}
\mathbf{E} &= \frac{1}{2} (\mathbf{g}_i \cdot \mathbf{g}_j - \mathbf{G}_i \cdot \mathbf{G}_j) \mathbf{G}^i \otimes \mathbf{G}^j \nonumber \\
&= E_{ij} \mathbf{e}_i \otimes \mathbf{e}_j = \hat{E}_{ij} \mathbf{G}^i \otimes \mathbf{G}^j \tag{16} \\
&= \frac{1}{2} \left( \mathbf{G}_i \cdot \frac{\partial \mathbf{u}}{\partial \xi^j} + \frac{\partial \mathbf{u}}{\partial \xi^i} \cdot \mathbf{G}_j + \frac{\partial \mathbf{u}}{\partial \xi^j} \cdot \frac{\partial \mathbf{u}}{\partial \xi^i} \right) \mathbf{G}^i \otimes \mathbf{G}^j \nonumber
\end{align}
$$

E<sub>ij</sub> and ˆE<sub>ij</sub> being respectively the components of the Green
Lagrange tensor in the Cartesian and covariant frame. In
Voigt notation the Cartesian and covariant Green–Lagrange
strain are related as follows

$$
\mathbf{E} = \mathbf{T}\hat{\mathbf{E}} \tag{17}
$$

with

$$
\begin{align}
\mathbf{E} &= \{E_{11}, E_{22}, E_{33}, 2E_{12}, 2E_{23}, 2E_{13}\}^t \text{ and} \nonumber \\
\hat{\mathbf{E}} &= \{\hat{E}_{\xi\xi}, \hat{E}_{\eta\eta}, \hat{E}_{\xi\zeta}, 2\hat{E}_{\xi\eta}, 2\hat{E}_{\eta\zeta}, 2\hat{E}_{\xi\zeta}\}^t \tag{18}
\end{align}
$$

and **T** a second-order matrix which contains the terms of the
inverse of the jacobian **J** matrix as follows

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

From Eq. (9) the Jacobian matrix can be decomposed into
constant, linear and bi-linear components depending on
ξ, η, ζ, ξη, ζη, ζξ

$$
\mathbf{J} = \mathbf{J}^0 + \xi \mathbf{J}^\xi + \eta \mathbf{J}^\eta + \zeta \mathbf{J}^\zeta + \xi \eta \mathbf{J}^{\xi\eta} + \eta \zeta \mathbf{J}^{\eta\zeta} + \xi \zeta \mathbf{J}^{\xi\zeta} \quad (14)
$$

This procedure is required for the construction of the
hourglass strain field that stabilizes the SB9 element. The
components of Eq. (14) are given in “Appendix”.

From the variation of the first line of Eq. (16) we can write
the covariant strain-displacement **B̂**<sub>I</sub> at the node *I* (*I* =
1, . . . , 8), such that δ**Ê**<sub>*I*</sub> = **B̂**<sub>*I*</sub>δ**U**<sub>*I*</sub> as follows

$$
\hat{\mathbf{B}}_I = \begin{bmatrix}
N_{I,\xi} \mathbf{g}_1^t \\
N_{I,\eta} \mathbf{g}_2^t \\
N_{I,\zeta} \mathbf{g}_3^t \\
N_{I,\xi} \mathbf{g}_2^t + N_{I,\eta} \mathbf{g}_1^t \\
N_{I,\xi} \mathbf{g}_3^t + N_{I,\zeta} \mathbf{g}_1^t \\
N_{I,\eta} \mathbf{g}_3^t + N_{I,\zeta} \mathbf{g}_2^t
\end{bmatrix} \quad (20)
$$
<figure><img src="image_3.png" /><figcaption>Fig. 2 Configuration of the SB9 element</figcaption></figure>

By considering the Eqs. (8) et (14), $\hat{\mathbf{B}}_I$ can be decomposed into the reduced integrated part $\hat{\mathbf{B}}_I^{\text{ri}}$ and a stabilization part $\hat{\mathbf{B}}_I^{\text{stab}}$ as follows:

$$ \hat{\mathbf{B}}_I = \hat{\mathbf{B}}_I^{\text{ri}} + \hat{\mathbf{B}}_I^{\text{stab}} \qquad (21) $$

With

$$ \hat{\mathbf{B}}_I^{\text{ri}} = \hat{\mathbf{B}}_I^0 + \zeta \hat{\mathbf{B}}_I^{\zeta} \qquad (22) $$

And

$$ \hat{\mathbf{B}}_I^{\text{stab}} = \xi \hat{\mathbf{B}}_I^{\xi} + \eta \hat{\mathbf{B}}_I^{\eta} + \xi \eta \hat{\mathbf{B}}_I^{\xi\eta} + \eta \zeta \hat{\mathbf{B}}_I^{\eta\zeta} + \xi \zeta \hat{\mathbf{B}}_I^{\xi\zeta} \qquad (23) $$

The components of Eqs. (22) and (23) are the basic building blocks of all solid-shell elements and are detailed in “Appendix”.

# 5 Assumed natural strain

## 5.1 Assumed shear strain

One way to reduce shear locking is to use the ANS concept following the work of Bathe and Dvorkin [32,34]. The transverse shear strain is evaluated in the four-middle edges points of the mid-surface of the solid-shell. To interpolate through the volume certain authors multiplied the mid-plane shear strain with the well known Reissner function [30,31]. This kind of interpolation ensures a systematic respect of the statistical condition of shear stress in the element upper and lower faces. However in non-linear problem involving material non linearity and large deformation, that interpolation becomes obsolete since it does not guarantee the satisfaction of statistical condition in the shell faces. Beside, such interpolation does not allow the element to naturally handle frictional contact problems since the shear deformation is assumed to be nil in the element upper and lower faces. Also it has been

<figure><img src="image_4.png" /><figcaption>Fig. 3 Positions of ANS tying points</figcaption></figure>

shown that such interpolation lead to a ill conditioned rigidity matrix. To avoid those inconveniences we applied the ANS technique following the work of Cardoso [12], Alves de Souza [37,51], Schwarze [10]. In the covariant frame, each of the shear deformation $\hat{E}_{\eta\zeta}$ and $\hat{E}_{\xi\zeta}$ are interpolated using the equivalent shear deformation from four different tying points (Fig. 3). E(-1, 0, -1), F(1, 0, -1), G(1, 0, 1) and H(-1, 0, 1), for $\hat{E}_{\eta\zeta}$ and J(0, -1, -1), K(0, 1, -1), L(0, 1, 1) and M(0, -1, 1) for $\hat{E}_{\xi\zeta}$. Hence the assumed shear deformations of the element is computed as follows

$$ \begin{aligned} \hat{E}_{\xi\zeta}^{\text{ANS}} &= \frac{1}{4}(1 \mp \zeta)((1 - \eta)\hat{E}_{\xi\zeta}^{J} + (1 + \eta)\hat{E}_{\xi\zeta}^{K} \\ &\quad + (1 + \eta)\hat{E}_{\xi\zeta}^{L} + (1 - \eta)\hat{E}_{\xi\zeta}^{M}) \\ \hat{E}_{\eta\zeta}^{\text{ANS}} &= \frac{1}{4}(1 \mp \zeta)((1 - \xi)\hat{E}_{\eta\zeta}^{E} + (1 + \xi)\hat{E}_{\eta\zeta}^{F} \\ &\quad + (1 + \xi)\hat{E}_{\eta\zeta}^{G} + (1 - \xi)\hat{E}_{\eta\zeta}^{H}) \end{aligned} \qquad (24) $$

## 5.2 Assumed normal strain

In order to eliminate the curvature thickness locking, we follow the assumed pinching strain as suggested by Bischoff and Ramm [18], Betsch and Stein [17], Schwatz and Reese [13].

$$ \hat{E}_{\xi\zeta}^{\text{ANS}} = N_A \hat{E}_{\xi\zeta}^{A} + N_B \hat{E}_{\xi\zeta}^{B} + N_C \hat{E}_{\xi\zeta}^{C} + N_D \hat{E}_{\xi\zeta}^{D} \qquad (25) $$

$$ N_k = \frac{1}{4}(1 + \xi_k \xi)(1 + \eta_k \eta) \text{ for } k = A, \dots, D \qquad (26) $$

# 6 Introduction of the enhancement

As assessed in Sect. 3 the presence of the additional node has two main aims. First getting a linear normal strain component which, along with a full three-dimensional constitutive strain-stress behavior, allows to achieve similar results in bending cases as those obtained with the usual plane stress state hypothesis. For that, the ninth node DOF plays the role of an extra parameter essential for a quadratic interpolation of the displacement in the thickness direction. In this case the extra DOF is similar to the EAS parameter seen in many solid-shell finite element. The main goal of the EAS parameter is
to enrich the normal deformation so that to avoid Poisson thickness locking. However what we have discovered is that even though this parameter allows to enhanced the element deformations, its contribution to the pinching stress is little to none. When a shell structure under one or double sided pressure load, solid shell elements proposed in the literature can give accurate results for the displacement thanks to the EAS parameter however when we look at the stresses, mainly the normal stress, we can see that they are less accurate. The normal stress do not respect the Neumann boundary condition when the shell is under pressure load. One should put more elements through the shell thickness direction to improve the normal stress accuracy, which is limiting since the main goal of solid-shell element is to put only one element in the shell thickness direction. Hence, The second advantage of the additional DOF is that it has a physical meaning and, for instance, a strength equivalent to a normal pressure can be prescribed to improve the normal stress.

## 6.1 The ninth node

The additional central node is endowed with only one translation degree of freedom in the element thickness direction $\zeta$. The enhanced normal displacement in the thickness direction is written as follows

$$u_{\zeta}(\xi, \eta, \zeta) = u_{\zeta}^{H}(\xi, \eta, \zeta) + (1 - \zeta^2)u_{\zeta}^{9} \quad (27)$$

$u_{\zeta}^9$ being the one and only relative displacement of the ninth node, in through the thickness direction $g_3$ as proposed by ahmad [63], $u_{\zeta}^H$ the normal displacement of the eight node element considered alone. A new column is then added to the covariant matrix $\hat{\mathbf{B}}^H$ making it a 6 by 25 matrix $\hat{\mathbf{B}}$

$$\hat{\mathbf{B}} = [\hat{\mathbf{B}}^H \quad \hat{\mathbf{B}}^9] \quad (28)$$

with

$$\hat{\mathbf{B}}^9 = [0 \quad 0 \quad -2\zeta \quad 0 \quad 0 \quad 0]^T$$

One can then write the variation of the pinching strain in the covariant basis as follows:

$$\delta E_{\zeta\zeta} = \hat{\mathbf{B}}_3^H \delta u_{\zeta}^H - 2\zeta \delta u_{\zeta}^9 \quad (29)$$

$\hat{\mathbf{B}}_3^H$ being the third line of matrix $\hat{\mathbf{B}}^H$.

## 6.2 Equivalent generalized nodal pressure forces

The advantage of having real extra degrees of freedom (DOF) instead of simple parameters is the possibility to physically act on them. For example, the nodal forces equivalent to a

normal pressure are prescribed at the apexes but also on the extra node in order to get the proper normal stress distribution. This has been previously done for the element Q5TTS [30] and a similar method is presented here for the 9-node, 25-DOF, solid-shell element SB9. To easily find the distribution of the required forces at each node, let us consider a hexahedral element (Fig. 4) on which are defined two normal pressures $P_s$ in the upper face and $P_i$ in the lower face. From Eq. (29), assuming a small change of direction of covariant vector $g_3$ we can deduce the variation of the actual normal strain, in the Ahmad base frame as follows

$$\delta E_{33} = \frac{\delta u_9^+ - \delta u_9^-}{h} - 4\zeta \frac{\delta u_9}{h} \quad (30)$$

where $\delta u_9^+$ and $\delta u_9^-$ are respectively the variation of the average normal displacement of the upper and lower face of the element, $h$ the element average thickness and $\delta u_9$ the variation of the ninth node displacement [31]. Assuming the pinching stress $S_{33}$ is linear through the element thickness:

$$S_{33} = \frac{1}{2}(\zeta - 1)P_i - \frac{1}{2}(\zeta + 1)P_s \quad (31)$$

the virtual equilibrium condition in the thickness direction can then be written and the equivalent nodal force distribution deduced as follows:

$$\sum_{i=1}^{8} F_i \delta u_i + F_9 \delta u_9 = \int_{S_0} (P_i - P_s) \delta u_i \, dS_0 + \int_{\Omega_0} \delta E_{33} \cdot S_{33} \, d\Omega_0 \quad (32)$$

$$F_i = \frac{1}{4} P_i A_i, \quad i = 1, 2, 3, 4$$

$$F_j = \frac{1}{4} P_s A_s, \quad j = 5, 6, 7, 8 \quad (33)$$

$$F_9 = \frac{2}{3} (P_s A_s - P_i A_i)$$

with $A_i$ and $A_s$ being the element lower and upper surface area. It is well known that solid-shells with only one layer of elements over the thickness are not able to reproduce a transverse normal stress state which is equal to the applied facials pressures. The new formulation herein (33) improves this lack, so the normal stress gives the accurate values corresponding to the applied pressures in the boundaries with only one element in the thickness direction.
<figure><img src="image_5.png" /><figcaption>Fig. 4 Projected view of the SB9 under pressure</figcaption></figure>

# 7 Stabilization procedure

## 7.1 Green Lagrange strain stabilization

The goal of the stabilization procedure is to correct the rank deficiency of the stiffness matrix coming from the adopted reduced integration scheme. The reader is referred to the important work in this subject by Liu, Belytschko and al [47,52], Cardoso et al. [12,24], Alves de Souza [37,51], schwarze [13]. In some solid-shell formulation, the Jacobian matrix and its inverse are only evaluated at the element center, see for example [31,40,41]. Such approximation assumes that the real element can be represented by its equivalent parallelepiped. And this is quite accurate for thin and not very distorted meshes. However for highly initially distorted meshes that assumption shows lack of accuracy [13]. In order to take into consideration the realistic shape of the element and stabilize properly the stiffness matrix, a polynomial decomposition of the inverse Jacobian matrix is derived following the work of schwarze and Reese [13]. The matrix T of Eq. (19) is decomposed into constant and linear terms as follows

$$T \approx T^0 + \xi T^\xi + \eta T^\eta + \zeta T^\zeta \quad (34)$$

It has been shown in [13] that a very good accuracy is reached just using the constant and linear terms of T. Hence, the bilinear terms are ignored. In the same way, the covariant Green Lagrange strain tensor can be split into its constant, linear and bilinear terms and the Cartesian Green Lagrange strain can then be derived as follows:

$$ 
\begin{aligned}
E & \approx (T^0 + \xi T^\xi + \eta T^\eta + \zeta T^\zeta) \\
& \quad (\hat{E}^0 + \zeta \hat{E}^\xi + \xi \hat{E}^\xi + \eta \hat{E}^\eta + \xi \eta \hat{E}^\xi \eta + \xi \zeta \hat{E}^\xi \zeta + \eta \zeta \hat{E}^\eta \zeta) \\
& \approx \underbrace{E^0 + \zeta E^\xi + \zeta^2 E^\xi \zeta}_{\text{E}^{\text{ri}}} + \underbrace{\xi E^\xi + \eta E^\eta + \xi \eta E^\xi \eta + \xi \zeta E^\xi \zeta + \eta \zeta E^\eta \zeta}_{\text{E}^{\text{stab}}}
\end{aligned} 
\quad (35) $$

and

$$E = E^{\text{ri}} + E^{\text{stab}} \quad (36)$$

The strain-displacement gradient can be split the same way

$$ 
\begin{aligned}
B & \approx (T^0 + \xi T^\xi + \eta T^\eta + \zeta T^\zeta) \\
& \quad (\hat{B}^0 + \zeta \hat{B}^\xi + \xi \hat{B}^\xi + \eta \hat{B}^\eta + \xi \eta \hat{B}^\xi \eta + \xi \zeta \hat{B}^\xi \zeta + \eta \zeta \hat{B}^\eta \zeta) \\
& \approx \underbrace{B^0 + \zeta B^\xi + \zeta^2 B^\xi \zeta}_{\text{B}^{\text{ri}}} + \underbrace{\xi B^\xi + \eta B^\eta + \xi \eta B^\xi \eta + \xi \zeta B^\xi \zeta + \eta \zeta B^\eta \zeta}_{\text{B}^{\text{stab}}}
\end{aligned} 
\quad (37) $$

so that B can be written as

$$B = B^{\text{ri}} + B^{\text{stab}} \quad (38)$$

$E^{\text{stab}}$ and $B^{\text{stab}}$ represent the terms cancelled in the reduced integration. The process of stabilization consists in analytically restoring those terms into the stiffness matrix. Furthermore to eliminate volumetric locking, the B-bar approach [64] is adopted. Hence the hourglass counterpart of the stabilized strain end strain-displacement operators are split into there volumetric and deviatoric components, and only the deviatoric part are kept

$$B^{\text{stab}} = B_{\text{dev}}^{\text{stab}}(\xi, \eta, \zeta) + B_{\text{vol}}^{\text{stab}}(0, 0, 0) \quad (39)$$

Since no constant term is present in the expansion of $B^{\text{stab}}$ (see Eq. 37)

$$B_{\text{vol}}^{\text{stab}}(0, 0, 0) = 0 \quad (40)$$

and

$$B^{\text{stab}} = B_{\text{dev}}^{\text{stab}}(\xi, \eta, \zeta) \quad (41)$$

The stabilization components of the strain and strain-displacement matrices being well identified, a similar approach is used to identify the stabilization stress.

## 7.2 Second Piola Kirchhoff stress stabilization

From Eq. (36), since the two parts of the Green Lagrange strain tensor are orthogonal, the internal energy is split like

$$W(E) = W^{\text{ri}}(E^{\text{ri}}) + W^{\text{stab}}(E^{\text{stab}}) \quad (42)$$

Since second Piola Kirchhoff tensor is the derivative of the internal energy with respect to Green–Lagrange strain, the stabilization counterpart can be identified by simply writing the derivative:

$$S = \frac{\partial W(E)}{\partial E} \quad (43)$$

$$= \frac{\partial W^{\text{ri}}(E^{\text{ri}})}{\partial E^{\text{ri}}} \frac{\partial E^{\text{ri}}}{\partial E} + \frac{\partial W^{\text{stab}}(E^{\text{stab}})}{\partial E^{\text{stab}}} \frac{\partial E^{\text{stab}}}{\partial E} \quad (44)$$

$$= S^{\text{ri}}(E^{\text{ri}}) + S^{\text{stab}}(E^{\text{stab}}) \quad (45)$$
$\mathbf{S}^{\text{ri}}(\mathbf{E}^{\text{ri}})$ is Piola Kirschhoff obtained from the integration of constitutive law and $\mathbf{S}^{\text{stab}}(\mathbf{E}^{\text{stab}})$ Piola Kirchhoff stress for stabilization. The key point of the stabilization is to find an optimal way to evaluate that stabilization stress without integrating it so that the computation time is minimized. Now just as being explained before, to avoid volumetric locking, only the deviatoric part of the Piola Kirchhoff stabilization tensor is kept. Doing so, the stabilization of the stress is given as

$$ \mathbf{S}^{\text{stab}} = \mathbf{C}^{\text{stab}} : \mathbf{E}^{\text{stab}} \qquad (46) $$

$\mathbf{C}^{\text{stab}}$ being the deviatoric part of St Venant Kirchhoff material given as follows

$$ \mathbf{C}^{\text{stab}} = \mu^{\text{stab}} \begin{bmatrix} \frac{4}{3} & -\frac{2}{3} & -\frac{2}{3} & 0 & 0 & 0 \\ -\frac{2}{3} & \frac{4}{3} & -\frac{2}{3} & 0 & 0 & 0 \\ -\frac{2}{3} & -\frac{2}{3} & \frac{4}{3} & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix} \qquad (47) $$

With $\mu$ the shear modulus evaluated, in elastic material behavior, as $\mu^{\text{stab}} = E/2(1 + \nu)$, in which $E$ and $\nu$ are respectively the Young's modulus and the Poisson ratio. For inelastic materials, using the elastic shear modulus lead to an overestimation of the stabilization stress. To overcome such problem the secant modulus as defined in Belytschko and Bindeman [54] is used:

$$ 2\mu^{\text{stab}} = \sqrt{\frac{\pi\mathbf{S}}{\pi_E}}, \qquad (48) $$

where

$$ \pi_S = \frac{1}{2} \mathbf{S}_{\text{dev}} : \mathbf{S}_{\text{dev}} \quad , \quad \pi_E = \frac{1}{2} \mathbf{E}_{\text{dev}} : \mathbf{E}_{\text{dev}} \qquad (49) $$

These ways of evaluating the stabilization parameter $\mu^{\text{stab}}$ is very convenient because there is no need for consistent linearization of stabilization constitutive matrix, which is very interesting in terms of computation time. Moreover, the proposed definition of $\mathbf{C}^{\text{stab}}$ ensures an efficient and robust stabilization.

## 8 Linearization of the weak form

In order to solve the nonlinear variational equation, we use a Newton-Raphson iterative method through a sequence of linearization. Assuming that the external energy is displacement independent, we only give details of the linearization

of the internal energy.

$$ L[\delta\pi_{\text{int}}] = \int_{\Omega_0} (\Delta\mathbf{S} : \delta\mathbf{E} + \mathbf{S} : \Delta(\delta\mathbf{E}))d\Omega_0 \qquad (50) $$

with

$$ \begin{align} \Delta \mathbf{S} &= \frac{\partial \mathbf{S}}{\partial \mathbf{E}} : \Delta \mathbf{E} = \mathbf{D} : \Delta \mathbf{E}^{\text{ri}} + \mathbf{D} : \Delta \mathbf{E}^{\text{stab}} \nonumber \\ \Delta \mathbf{E} &= \Delta \mathbf{E}^{\text{ri}} + \Delta \mathbf{E}^{\text{stab}} = \mathbf{B}^{\text{ri}} \Delta \mathbf{U} + \mathbf{B}^{\text{stab}} \Delta \mathbf{U} \tag{51} \end{align} $$

$$ \delta \mathbf{E} = \delta \mathbf{E}^{\text{ri}} + \delta \mathbf{E}^{\text{stab}} = \mathbf{B}^{\text{ri}} \delta \mathbf{U} + \mathbf{B}^{\text{stab}} \delta \mathbf{U} $$

$$ \Delta(\delta\mathbf{E}) = \delta\mathbf{U}^t (\mathcal{B}^{\text{ri}} + \mathcal{B}^{\text{stab}}) \Delta\mathbf{U} $$

Equation (50) can be split into a material stiffness and a geometrical stiffness.

$$ \begin{align} \int_{\Omega_0} \Delta \mathbf{S} : \delta \mathbf{E} d\Omega_0 & \nonumber \\ &= \delta \mathbf{U} \left( \int_{\Omega_0} (\mathbf{B}^{\text{ri}})^t \mathbf{D} \mathbf{B}^{\text{ri}} d\Omega_0 + \int_{\Omega_0} (\mathbf{B}^{\text{stab}})^t \mathbf{D} \mathbf{B}^{\text{stab}} d\Omega_0 \right) \Delta \mathbf{U} \nonumber \\ &= \delta \mathbf{U}^t \left( \mathbf{K}_m^{\text{ri}} + \mathbf{K}_m^{\text{stab}} \right) \Delta \mathbf{U} \tag{52} \end{align} $$

and

$$ \begin{align} \int_{\Omega_0} \mathbf{S} : \Delta(\delta\mathbf{E}) d\Omega_0 & \nonumber \\ &= \delta\mathbf{U} \left( \int_{\Omega_0} \mathcal{B}^{\text{ri}} \mathbf{S}^{\text{ri}} d\Omega_0 + \int_{\Omega_0} \mathcal{B}^{\text{stab}} \mathbf{S}^{\text{stab}} d\Omega_0 \right) \Delta\mathbf{U} \tag{53} \\ &= \delta\mathbf{U}^t \left( \mathbf{K}_g^{\text{ri}} + \mathbf{K}_g^{\text{stab}} \right) \Delta\mathbf{U} \nonumber \end{align} $$

With

$$ \begin{align} \mathbf{K}_m^{\text{ri}} &= 4J_0 \int_{-1}^{1} (\mathbf{B}^0 + \zeta \mathbf{B}^{\zeta} + \zeta^2 \mathbf{B}^{\zeta^2})^t \mathbf{D}(\mathbf{B}^0 + \zeta \mathbf{B}^{\zeta} + \zeta^2 \mathbf{B}^{\zeta^2}) d\zeta \nonumber \\ \mathbf{K}_m^{\text{stab}} &= J_0 [\mathbf{B}^{\zeta\eta t} \mathbf{D}\mathbf{B}^{\zeta\eta} + \mathbf{B}^{\zeta t} \mathbf{D}\mathbf{B}^{\zeta\zeta}] \frac{8}{3} \nonumber \\ &\quad + J_0 [\mathbf{B}^{\zeta\eta t} \mathbf{D}\mathbf{B}^{\zeta\eta} + \mathbf{B}^{\zeta t} \mathbf{D}\mathbf{B}^{\zeta\zeta} + \mathbf{B}^{\eta\zeta t} \mathbf{D}\mathbf{B}^{\eta\zeta}] \frac{8}{9} \tag{54} \\ \mathbf{K}_g^{\text{ri}} &= J_0 \int_{-1}^{1} (\mathcal{B}^0 + \zeta \mathcal{B}^{\zeta} + \zeta^2 \mathcal{B}^{\zeta^2}) \mathbf{S}^{\text{ri}} d\zeta \nonumber \\ \mathbf{K}_g^{\text{stab}} &= J_0 [\mathcal{B}^{\zeta\eta} \mathbf{S}^{\zeta\eta} + \mathcal{B}^{\eta\zeta} \mathbf{S}^{\eta\zeta} + \mathcal{B}^{\zeta\zeta} \mathbf{S}^{\zeta\zeta}] \frac{8}{9} \nonumber \end{align} $$

$$ \begin{align} \mathbf{K}_g^{\text{stab}} &= J_0 [\mathcal{B}^{\xi\eta}\mathbf{S}^{\xi\eta} + \mathcal{B}^{\eta\xi}\mathbf{S}^{\eta\xi} + \mathcal{B}^{\xi\zeta}\mathbf{S}^{\xi\zeta}] \frac{8}{9} \nonumber \\ &\quad + J_0 [\mathcal{B}^{\xi\eta}\mathbf{S}^{\xi\eta} + \mathcal{B}^{\eta\xi}\mathbf{S}^{\eta\xi} + \mathcal{B}^{\xi\zeta}\mathbf{S}^{\xi\zeta}] \frac{8}{9} \end{align} $$

$\mathbf{K}_m^{\text{ri}}$ and $\mathbf{K}_m^{\text{stab}}$ are respectively the integrated stiffness matrix and its stabilization counterpart, $\mathbf{K}_g^{\text{ri}}$ and $\mathbf{K}_g^{\text{stab}}$ respectively the geometrical stiffness matrix and its counterpart. $J_0$ is the determinant of the Jacobian matrix. Note that the stabilization matrices are computed analytically without numerical integration. In the same way, one can write the internal force and its stabilization counterpart as
$$ 
\begin{aligned}
\mathbf{F}_{int}^{ri} &= 4J_0 \int_{-1}^{1} (\mathbf{B}^0 + \zeta \mathbf{B}^{\zeta} + \zeta^2 \mathbf{B}^{\zeta^2})^T \mathbf{S} \mathrm{d}\zeta \\
\mathbf{F}_{int}^{stab} &= J_0 [\mathbf{B}^{\xi t} \mathbf{S}^{\xi} + \mathbf{B}^{\eta t} \mathbf{S}^{\eta}] \frac{8}{3} \\
&\quad + J_0 [\mathbf{B}^{\xi \eta t} \mathbf{S}^{\xi \eta} + \mathbf{B}^{\xi \zeta t} \mathbf{S}^{\xi \zeta} + \mathbf{B}^{\eta \zeta t} \mathbf{S}^{\eta \zeta}] \frac{8}{9}
\end{aligned}
\quad (55) $$

## 9 Numerical implementation

In this section, the main features of the element implementation is briefly described. The element has been implemented in the implicit nonlinear finite element code Code_aster. In this process, the total Lagrangian formulation is adopted making it obsolete the use of a co-rotational formulation. The equilibrium equations are solved step by step using an iterative procedure based on the Newton-Raphson scheme. These iterations are performed until the residual load vector is sufficiently small, using a constant tangent stiffness matrix built at the beginning of the current time step. It's worth noting that the stabilization process borrowed from the work of Schwarze and Reese [13] shows a very good robustness allowing some structural instabilities problems involving either a load-limit point ('snap-through') or a deflection-limit point ('snap-back') to be overcome simply with a Newton-Raphson algorithm without adopting arc-length control parameter. The first two tests intend to show the importance of the enhancement of the normal stress which is possible thank to the additional middle node. The following tests are classical tests available in the literature to show the performance of the element in bending problems, in large deformation with or without material non-linearity. The different finite elements used as reference for comparison are:

- Q1STs: Reduced integrated solid-shell element with EAS and ANS by Schwarze [13]

- Q5TTS: 5-nodes quadrilateral with through thickness stress by Sansalone [30]

<figure><img src="image_6.png" /><figcaption>Fig. 5 Circular clamped plate</figcaption></figure>

- SHB8PS: Reduced integrated solid-shell with ANS by Abed Meraim and Combescure [40]

- S4R: ABAQUS's four-node shell element models [65]

## 9.1 Circular clamped plate

The goal of this test is to show the main advantage of the ninth node compared to a EAS parameter alone. We consider a disk of radius $R = 100$ mm and of thickness $t = 1$ mm. The material is isotropic. The Young modulus is $E = 2 \times 10^5$ MPa and the Poisson ratio $\nu = 0.3$. Two studies have been done. In the first study a uniform normal pressure $P_u = 0.01172$ N/mm² is applied in the upper side of the plate. In the second study two normal pressures $P_u = 2 \times 0.01172$ N/mm² and $P_l = 0.01172$ N/mm² have been respectively applied in the upper and lower faces of the circular plate. The theoretical displacement of the plate center (C in Fig. 5) is $(U)_C = 1$mm according to Kirchhoff's theory. Tables 1 and 2 show respectively the results for study 1 and study 2 from different finite elements. It is interesting to notice that all different finite elements give good result if we look at the plate center displacement. However the interesting thing to notice is the pinching stress given by the different elements. The SHB8PS and the S4R of abaqus give a nil value for the pinching stress, which is normal since they work with plane stress state. The EAS-only (here used with 5 Gauss points) and the SB9 are both using a full 3D constitutive law but we can see a non negligible difference in the pinching stress. In Tables 3 and 4 the pinching stress for the two elements are detailed in order to have a broad view of what is happening in the shell thickness. The pinching stress is almost constant in the thickness direction for the EAS solid-shell, while for the SB9 element we have a better normal stress which satisfy the boundary condition. This is the main benefice of utilizing a additional node compared to a simple EAS parameter. It allows to split the applied pressure so that the pinching stress is more accurate.

Thanks to their additional nodes the Q5TTS and the SB9 use a 3D constitutive law and give a pinching stress that is very accurate, on the contrary of the other ele-
<table><caption>**Table 1** Center displacement and normal stress</caption><thead><tr><th></th><th>Q5TTS</th><th>S4R</th><th>SHB8PS</th><th>SB9</th><th>EAS only</th></tr></thead><tbody><tr><th>Uc</th><td>0.9874</td><td>0.995</td><td>0.995</td><td>0.995</td><td>0.995</td></tr><tr><th>σ<sup>top</sup><sub>bottom</sub></th><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>-0.0058</td></tr><tr><th>σ<sup>bottom</sup></th><td>-0.01172</td><td>0.0</td><td>0.0</td><td>-0.01172</td><td>-0.0058</td></tr></tbody></table>

<table><caption>**Table 2** Center displacement and normal stress</caption><thead><tr><th></th><th>Q5TTS</th><th>S4R</th><th>SHB8PS</th><th>SB9</th><th>EAS only</th></tr></thead><tbody><tr><th>Uc</th><td>0.9874</td><td>0.995</td><td>0.995</td><td>0.995</td><td>0.995</td></tr><tr><th>σ<sup>top</sup><sub>bottom</sub></th><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>-0.0058</td></tr><tr><th>σ<sup>bottom</sup></th><td>-0.01172</td><td>0.0</td><td>0.0</td><td>-0.01172</td><td>-0.0058</td></tr></tbody></table>

<table><caption>**Table 3** Pinching stress at Gauss point: study 1</caption><thead><tr><th>σ at Gauss points</th><th>EAS-only</th><th>SB9</th></tr></thead><tbody><tr><td>1</td><td>-0.00585</td><td>-0.00000</td></tr><tr><td>2</td><td>-0.00585</td><td>-0.002021</td></tr><tr><td>3</td><td>-0.00585</td><td>-0.005857</td></tr><tr><td>4</td><td>-0.00585</td><td>-0.009693</td></tr><tr><td>5</td><td>-0.00585</td><td>-0.011717</td></tr></tbody></table>

<table><caption>**Table 4** Pinching stress at Gauss point: study 2</caption><thead><tr><th>σ at Gauss points</th><th>EAS-only</th><th>SB9</th></tr></thead><tbody><tr><td>1</td><td>-0.01757</td><td>-0.01171</td></tr><tr><td>2</td><td>-0.01757</td><td>-0.01373</td></tr><tr><td>3</td><td>-0.01757</td><td>-0.01757</td></tr><tr><td>4</td><td>-0.01757</td><td>-0.02140</td></tr><tr><td>5</td><td>-0.01757</td><td>-0.02343</td></tr></tbody></table>

ments that use the plan stress assumption or only the EAS parameter.

## 9.2 Under pressure cylinder

This test case, like the previous one, is designed to study the pinch response of the SB9 element but this time with a curved structure. We consider a cylinder under internal and external pressure and evaluate the normal stress (pinching stress). For symmetry reasons only one eighth of the cylinder is meshed, the Fig. 6 being obtained by several reflections relative to the axes of symmetries. The geometric and material characteristics of the test are displayed in Table 5. The results are compared to an approximate analytical solution. The results are compared to an approximate analytical solution:

<table><caption>**Table 5** Geometric and material characteristics</caption><tr><td>Radius</td><td>R</td><td>1</td></tr><tr><td>Thickness</td><td>t</td><td>0.01</td></tr><tr><td>Young modulus</td><td>E</td><td>1 × 10<sup>7</sup></td></tr><tr><td>Poisson ratio</td><td>ν</td><td>0.0</td></tr><tr><td>Internal pressure</td><td>Pi</td><td>2</td></tr><tr><td>External pressure</td><td>Pe</td><td>1</td></tr></table>

$$ \left\{ \begin{array}{l} \sigma_{rr}(r) = \left( P \frac{R}{2t} - \bar{P} \right) - P \frac{R}{2t} \left( \frac{R}{r} \right)^2 \\ P = P_i - P_e \\ \bar{P} = \frac{1}{2} (P_i + P_e) \end{array} \right. \qquad (56) $$

$P_i$ and $P_e$ are the pressures applied on the inner and outer walls of the cylinder. $R$ and $t$ represent the radius and thickness of the cylinder. In Table 6 we can observe the values of the normal stress, compared to those given by the analytical solution of Eq. 56.

In this test also, we find the static admissibility $\sigma_{rr} = [\sigma \cdot n] \cdot n = -p$, where $\sigma$ is the stress field, $n$ the normal on the face and $p$ the applied pressure. We can thus see the interest of the additional term compared with a EAS parameter alone without enhancing the pinching stress.

## 9.3 Bending of a cantilever beam

In this test we investigate the out-of-plane bending of a cantilever beam under a tip load. The problem has been analyzed by many authors [8,10,48,59,66,67]. As shown in Fig. 7, the beam has a side length of $L = 10$ mm and a rectangular cross section of width $B = 1$ mm and a thickness $t = 0.1$ mm. The material is isotropic and has a Young modulus $E = 10^7$ N/mm$^2$ and a Poisson ratio $\nu = 0.3$. The total load is $F_{max} = 40$ N introduced in ten time increments. The Beam is discretized by one element along the width and thickness direction and 16 elements along the length. The numerical integration is performed with Gauss points through the thickness direction. A first study has been made to assess the convergence of the SB9 for this test, using a regular mesh and then a highly distorted mesh. Figure 8a gives the load displacement path for the two different meshes. The displacement at convergence is in agreement with the analytic solution given in [68]. The load displacement path is identical for both the regular and distorted mesh, proving that the SB9 is very insensitive to mesh distortion. In a second time the behaviour of the element in the near quasi-incompressible situations is investigated. The regular mesh is kept but this time the Poisson ratio is varied between $\nu = 0.0$ and $\nu = 0.4999$. The element still gives good result showing no sign of volumetric locking, see Fig. 8b.
<table><caption>**Table 6** Pinch stress on the 5 Gauss points</caption><thead><tr><th>Pts Gauss</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>σ<sub>rr</sub> Analyt</td><td>-2.00378</td><td>-1.82894</td><td>-1.5</td><td>-1.17427</td><td>-1.00373</td></tr><tr><td>σ<sub>rr</sub> SB9</td><td>-1.99915</td><td>-1.82648</td><td>-1.49915</td><td>-1.17183</td><td>-0.99915</td></tr><tr><td>Érreur (%)</td><td>0.04</td><td>0.00</td><td>0.05</td><td>0.02</td><td>0.45</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 6** Cylinder's mesh</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7** Bending of a cantilever beam</figcaption></figure>

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

## 9.4 Spherical shell with 18 degree hole

In this example the hemispherical shell with an 18 degree circular cut out at its pole is studied. The shell is loaded by an alternating radial point forces F at 90 degree intervals, see Fig. 9. This is a well known problem as it has been considered by many authors [10,15,69–71], to name just a few. Thanks to symmetry, only one-quater of the shell is modeled with a 16 × 16 mesh. The geometrical parameters used are the same used by [71] with a radius of $R = 10.0$ mm and a thickness of $t = 0.04$ mm which give a ratio $R/t = 250$. So this is a relatively severe test more likely to exhibit locking effects. The structure is subdued to two concentrated force $F = 100N$ as shown in Fig. 9. The material parameters are the Young modulus $E = 6.825 \times 10^7$ N/mm² and the Poisson ratio $\nu = 0.3$. The total number of increment to reach the

<figure><img src="image_10.png" /><figcaption>(a) Distortion analysis</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) near-incompressibility analysis</figcaption></figure>

**Fig. 8** Evolution of the cantilever beam's free-edge displacement

maximum load is 10. In Fig. 10, the deformation path of the points A and B are plotted and compared to the values given by Sze [71]. The results given by the SB9 element is in good agreement with the reference solutions.
<figure><img src="image_12.png" /><figcaption>**Fig. 9** Pinched hemisphere with hole</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 10** Evolution of displacement of hemisphere</figcaption></figure>

## 9.5 Stretched cylinder with free edge

Loaded by two opposite singles forces, see Fig. 11, the cylinder undergoes a significantly large rotation combining bending and membrane effects. The geometric is defined

by a length $L = 10.35$, a radius $R = 4.953$ and a thickness $T = 0.094$ mm. The material properties are given by the Young modulus $E = 10.5 \times 10^6$ N/mm$^2$ and the Poisson ratio $\nu = 0.3125$ and the applied load in each side is $F = 40$ KN. Due to the symmetry of the problem, only one eighth of the system is discretized. In order to investigate the convergence several meshes refinements are tested including 8 × 12, 16 × 24 and 20 × 30, with only a single element through the thickness. Using the present formulation of the SB9 the problem is modeled using a 16 × 24 mesh and the results are compared to those tabulated in Sze et al [71]. Figure 12 shows the results in term of load versus radial displacements at points A, B and C. Point A corresponds to the point under the loading, while point B and C are on the side of the cylinder and undergo horizontal displacements see Fig. 12. As it can be seen in the load-displacement curves, the overall response exhibits two regimes: a first stage dominated by bending effect and characterized by large displacements and rotations, and a second phase dominated by membrane effects, which may cause locking. It's also important to note the snap-through phenomenon arising when the loading reaches the critical value of around 20 kN. This can be seen through the displacement reversal that occurs on the

<figure><img src="image_14.png" /><figcaption>**Fig. 11** Stretched cylinder with free edge</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. 12** Evolution of displacement of stretched cylinder</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 14** Evolution of displacement of pinched cylinder</figcaption></figure>

load displacement curve of point C. This new formulation of the SB9 element gives results in very good agreement with the reference solution, with no locking effect.

## 9.6 Pinched cylindrical shell mounted over rigid diaphragms

This problem and its variations have been considered by many authors [3,69–71], among others. The cylinder shell is mounted on a rigid diaphragms and is subjected to two pinching forces load $P_{max}$ = 12000 as shown is the Fig. 13a. The dimensions of the cylinder are a length of $L$ = 200, a radius $R$ = 100 and a thickness $h$ = 1. The material properties are Young Modulus $E$ = $30 \times 10^3$ and the Poisson ratio $\nu$ = 0.3. Thanks to symmetry, one-eighth of the shell is modeled. A commonly employed mesh of 40 × 40 solid-shell elements is used in this test. Figure 14 shows the displacement of the

points A and B, compared with those given in [71] with the Abaqus S4R element.

## 9.7 Simply supported plate

In this example the simply supported plate under pressure is addressed. The problem has been treated by many authors [8, 10, 59, 72, 73] among others. The geometric parameters of the plate are represented by the side length $L$ = 508 mm and the thickness $t$ = 2.54. Thanks to symmetry only one quarter of the plate is represented, Fig. 15b being obtained by reflexion. As boundary conditions, the outer lower side edge of the plate is fixed in the vertical direction ($u_3$ = 0), hence the rotation around these axes is possible. A uniform pressure $p = \lambda p_0$ with $p = 0.01 \text{N/mm}^2$ is applied. The maximum of the load factor is $\lambda_{max}$ = 40. The material is taken ideally elastic-plastic with $E$ = $6.9 \times 10^4 \text{N/mm}^2$, $\nu$ = 0.3 and the yield stress $\sigma_y$ = 248N/mm$^2$. Figure 16 gives the evolution of the normal displacement of the plate center compared to the results given in [10].

<figure><img src="image_17.png" /><figcaption>(a) Geometric description</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(b) Visual results</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>(a) Geometric description</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(b) Visual results</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 16** Evolution of the normal displacement of the plate center</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 18** Hardening law</figcaption></figure>

## 9.8 Stamping of a sheet by a cylindrical punch

This test represents a calculation of the drawing of a square sheet metal by a rigid cylindrical punch in the presence of large plastic deformations. This test is useful in the simulation of sheet metal forming. The sheet metal is modelled in SB9 elements. The punch and die are rigid elements. The contact between the different elements is modelled by a continuous formulation. The side length of the sheet is $L = 160$ mm and thickness $e = 2$ mm. The die is $h = 35$ mm deep with

a diameter of $D = 100$ mm. The geometric details of the tools are visible in the Fig. 17. The elastic parameters are the Young modulus $E = 200,000$ MPa and the Poisson ratio $v = 0.3$ and the hardening law is plotted in Fig. 18. To stamp the sheet, the punch is imposed with a vertical displacement from 0 to 35 mm. Figure 19 shows the deformation state of the sheet. This concludes that for these kind of problems the SB9 element is very robust, gives excellent results and is quite efficient.

<figure><img src="image_23.png" /><figcaption>**Fig. 17** Stamping of a sheet by a cylindrical punch, tools geometry</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Fig. 19 Stamping of a sheet by a cylindrical punch, deformation state</figcaption></figure>

## 10 Conclusion

In this paper we have presented a nine nodes solid-shell finite element with a improved pinching stress. A ninth node, with one translational degree of freedom through the element thickness, is added in the middle of the classical hexahedral solid-shell to allow the unrestricted use of 3D constitutive law without encountering Poisson thickness locking. Furthermore the adjunction of this middle node allows a redistribution of a pinching pressure force between a surface contribution and a volume contribution. This improves considerably the normal stress and gives better results than elements with the EAS parameter alone, without pinching stress enrichment. In addition to eliminating pinching and volumetric locking thanks to the middle node, transverse and trapezoidal shear locking are reduced by the assumed strain method applied to both transverse shear and pinch strain. To avoid zero energy modes, the stabilization terms are isolated and analytically integrated considering only their deviatoric part following the B-bar method. The finite element formulation has been subjected to a number of representative numerical assessments and the results are very accurate. The SB9 finite element is implemented into the quasi-static implicit software code_aster developed by the French energy and electricity company (eDF).

**Acknowledgements** The authors gratefully acknowledge the National Association of Research and Technology (ANRT) and électricité de France (eDF) for the funding of this work.

## 11 Appendix

$$ \mathbf{J}^0 = \frac{1}{8} \begin{bmatrix} a_1^I x_1^I & a_1^I x_2^I & a_1^I x_3^I \\ a_2^I x_1^I & a_2^I x_2^I & a_2^I x_3^I \\ a_3^I x_1^I & a_3^I x_2^I & a_3^I x_3^I \end{bmatrix} \quad (57) $$

$$ \mathbf{J}^\xi = \frac{1}{8} \begin{bmatrix} 0 & 0 & 0 \\ h_1^I x_1^I & h_1^I x_2^I & h_1^I x_3^I \\ h_3^I x_1^I & h_3^I x_2^I & h_3^I x_3^I \end{bmatrix} \quad (58) $$

$$ \mathbf{J}^\eta = \frac{1}{8} \begin{bmatrix} h_1^I x_1^I & h_1^I x_2^I & h_1^I x_3^I \\ 0 & 0 & 0 \\ h_2^I x_1^I & h_2^I x_2^I & h_2^I x_3^I \end{bmatrix} \quad (59) $$

$$ \mathbf{J}^\zeta = \frac{1}{8} \begin{bmatrix} h_3^I x_1^I & h_3^I x_2^I & h_3^I x_3^I \\ h_2^I x_1^I & h_2^I x_2^I & h_2^I x_3^I \\ 0 & 0 & 0 \end{bmatrix} \quad (60) $$

$$ \mathbf{J}^{\xi\eta} = \frac{1}{8} \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ h_4^I x_1^I & h_4^I x_2^I & h_4^I x_3^I \end{bmatrix} \quad (61) $$

$$ \mathbf{J}^{\xi\zeta} = \frac{1}{8} \begin{bmatrix} 0 & 0 & 0 \\ h_4^I x_1^I & h_4^I x_2^I & h_4^I x_3^I \\ 0 & 0 & 0 \end{bmatrix} \quad (62) $$

$$ \mathbf{J}^{\eta\zeta} = \frac{1}{8} \begin{bmatrix} h_4^I x_1^I & h_4^I x_2^I & h_4^I x_3^I \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \quad (63) $$

Here are given the in-plane part of the parameter of equations (22) and (23). The normal and shear part being sorted from Sect. 5.
$$ \hat{\mathbf{B}}_{mI}^{0} = \begin{bmatrix} a_{1I} \mathbf{J}_{1}^{0} \\ a_{2I} \mathbf{J}_{2}^{0} \\ a_{1I} \mathbf{J}_{2}^{0} + g_{2I} \mathbf{J}_{1}^{0} \end{bmatrix} \qquad (64) $$

$$ \hat{\mathbf{B}}_{mI}^{\xi} = \begin{bmatrix} h_{3I} \mathbf{J}_{1}^{0} + a_{1I} \mathbf{J}_{1}^{\xi} \\ h_{2I} \mathbf{J}_{2}^{0} + a_{2I} \mathbf{J}_{2}^{\xi} \\ h_{2I} \mathbf{J}_{1}^{0} + h_{3I} \mathbf{J}_{2}^{0} + a_{2I} \mathbf{J}_{1}^{\xi} + a_{1I} \mathbf{J}_{2}^{\xi} \end{bmatrix} \qquad (65) $$

$$ \hat{\mathbf{B}}_{mI}^{\xi\xi} = \begin{bmatrix} h_{3I} \mathbf{J}_{1}^{\xi} \\ h_{2I} \mathbf{J}_{2}^{\xi} \\ h_{2I} \mathbf{J}_{1}^{\xi} + h_{3I} \mathbf{J}_{2}^{\xi} \end{bmatrix} \qquad (66) $$

$$ \hat{\mathbf{B}}_{mI}^{\xi} = \begin{bmatrix} 0 \\ h_{1I} \mathbf{J}_{2}^{0} + a_{2I} \mathbf{J}_{2}^{\xi} \\ h_{1I} \mathbf{J}_{2}^{0} + a_{1I} \mathbf{J}_{2}^{\xi} \end{bmatrix} \qquad (67) $$

$$ \hat{\mathbf{B}}_{mI}^{\eta} = \begin{bmatrix} h_{1I} \mathbf{J}_{1}^{0} + a_{1I} \mathbf{J}_{1}^{\eta} \\ 0 \\ h_{1I} \mathbf{J}_{2}^{0} + a_{2I} \mathbf{J}_{1}^{\eta} \end{bmatrix} \qquad (68) $$

$$ \hat{\mathbf{B}}_{mI}^{\eta\xi} = \begin{bmatrix} h_{4I} \mathbf{J}_{1}^{0} + h_{3I} \mathbf{J}_{1}^{\eta} + h_{1I} \mathbf{J}_{1}^{\xi} + a_{1I} \mathbf{J}_{1}^{\eta\xi} \\ 0 \\ h_{4I} \mathbf{J}_{2}^{0} + h_{2I} \mathbf{J}_{1}^{\eta} + h_{1I} \mathbf{J}_{2}^{\xi} + g_{2I} \mathbf{J}_{1}^{\eta\xi} \end{bmatrix} \qquad (69) $$

$$ \hat{\mathbf{B}}_{mI}^{\xi\xi} = \begin{bmatrix} 0 \\ h_{4I} \mathbf{J}_{2}^{0} + h_{2I} \mathbf{J}_{2}^{\xi} + h_{1I} \mathbf{J}_{2}^{\xi} + a_{2I} \mathbf{J}_{2}^{\xi\xi} \\ h_{4I} \mathbf{J}_{1}^{0} + h_{3I} \mathbf{J}_{2}^{\xi} + h_{1I} \mathbf{J}_{1}^{\xi} + a_{1I} \mathbf{J}_{2}^{\xi\xi} \end{bmatrix} \qquad (70) $$

$\mathbf{J}_1, \mathbf{J}_2, \mathbf{J}_3$ being the rows of the Jacobian matrix in the current configuration. Note that the rows of the jacobian matrix constrain the coordinates of the covariant $g_i$ vectors.

## 11.1 Decomposition of T

The Inverse of the Jacobian matrix is decomposed keeping only the constant and linear terms as follows

$$ \mathbf{J}^{-1} \approx \mathbf{J}^{-1}|_{\xi=0} + \sum_{i=1}^{3} \mathbf{J}_{,\xi_i}^{-1}|_{\xi=0} \xi_i \qquad (71) $$

The constant term being easily determined, the work will be to determine the linear terms with very limited resources. To do so, the Eq. (72) is simply derived with respect to the corresponding convective parameter and gives the Eq. (73).

$$ \mathbf{J}\mathbf{J}^{-1} \approx (\mathbf{J}\mathbf{J}^{-1})|_{\xi=0} + \sum_{i=1}^{3} (\mathbf{J}\mathbf{J}^{-1})_{,\xi_i}|_{\xi=0} \xi_i \qquad (72) $$

$$ 0 = (\mathbf{J}\mathbf{J}^{-1})_{,\xi_i} |_{\xi=0} = \mathbf{J}_{,\xi_i} |_{\xi=0} \mathbf{J}^{-1} |_{\xi=0} + \mathbf{J} |_{\xi=0} \mathbf{J}_{,\xi_i}^{-1} |_{\xi=0} \qquad (73) $$

From Eq. (73) one can easily determine the linear Jacobian terms, as follows, all terms being known or determined easily.

$$ \mathbf{J}_{,\xi_i}^{-1}|_{\xi=0} = -(\mathbf{J}^0)^{-1} \mathbf{J}^{\xi_i} (\mathbf{J}^0)^{-1} \qquad (74) $$

Hence, a good representation of the inverse Jacobian matrix is known and one can simply insert these terms into Eq. (19) to sort the taylor decomposition of the matrix T

$$ \mathbf{J}^{-1} \approx \mathbf{J}^{-1}|_{\xi=0} - \sum_{i=1}^{3} (\mathbf{J}^0)^{-1} \mathbf{J}^{\xi_i} (\mathbf{J}^0)^{-1} \xi_i \qquad (75) $$

## References

1. Parisch H (1995) A continuum-based shell theory for non-linear applications. Int J Numer Meth Eng 38(11):1855–1883

2. 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

3. 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 Meth Eng 42(1):49–69

4. Hauptmann R, Schweizerhof K, Doll S (2000) Extension of the ‘solid-shell’ concept for application to large elastic and large elasto-plastic deformations. Int J Numer Meth Eng 49(9):1121–1141

5. Sze KY, Yao LQ (2000) A hybrid stress ans solid-shell element and its generalization for smart structure modelling. part i-solid-shell element formulation. Int J Numer Meth Eng 48(4):545–564

6. Vu-Quoc L, Tan XG (2003) Optimal solid shells for non-linear analyses of multilayer composites. i. statics. Comput Methods Appl Mech Eng 192(9–10):975–1016

7. Legay A, Combescure A (2003) Elastoplastic stability analysis of shells using the physically stabilized finite element SHB8PS. Int J Numer Meth Eng 57(9):1299–1322

8. Fontes Valente RA, Parente MPL, Natal Jorge RM, César de Sá JMA, Grácio JJ (2005) Enhanced transverse shear strain shell formulation applied to large elasto-plastic deformation problems. Int J Numer Meth Eng 62(10):1360–1398

9. Kim KD, Liu GZ, Han SC (2005) A resultant 8-node solid-shell element for geometrically nonlinear analysis. Comput Mech 35(5):315–331

10. Reese S (2007) A large deformation solid-shell concept based on reduced integration with hourglass stabilization. Int J Numer Meth Eng 69(8):1671–1716

11. Alves De Sousa RJ, Yoon JW, Cardoso RPR, Fontes Valente RA, Gracio JJ (2007) On the use of a reduced enhanced solid-shell (RESS) element for sheet forming simulations. Int J Plast 23(3):490–515

12. Cardoso RPR, Yoon JW, Mahardika M, Choudhry S, Alves de Sousa RJ, Fontes Valente RA (2008) Enhanced assumed strain (EAS) and assumed natural strain (ANS) methods for one-point quadrature solid-shell elements. Int J Numer Methods Eng 75(2):156–187

13. 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 Meth Eng 80(10):1322–1355

14. Masud A, Tham CL (2000) Three-dimensional corotational framework for elasto-plastic analysis of multilayered composite shells. AIAA J 38(12):2320–2327
15. Simo JC, Fox DD, Rifai MS (1990) On a stress resultant geometrically exact shell model. Part iii: computational aspects of the nonlinear theory. Comput Methods Appl Mech Eng 79(1):21-70

16. Braun M, Bischoff M, Ramm E (1994) Nonlinear shell formulations for complete three-dimensional constitutive laws including composites and laminates. Comput Mech 15(1):1-18

17. 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(1-2):57-79

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

19. Kemp BL, Cho C, Lee SW (1998) A four-node solid shell element formulation with assumed strain. Int J Numer Meth Eng 43(5):909-924

20. Bischoff M, Ramm E (2000) On the physical significance of higher order kinematic and static variables in a three-dimensional shell formulation. Int J Solids Struct 37(46-47):6933-6960

21. El-Abbasi N, Meguid SA (2000) A new shell element accounting for through-thickness deformation. Comput Methods Appl Mech Eng 189(3):841-862

22. Cardoso RPR, Yoon J-W, Grácio JJ, Barlat F, César de Sá JMA (2002) Development of a one point quadrature shell element for nonlinear applications with contact and anisotropy. Comput Methods Appl Mech Eng 191(45):5177-5206

23. Brank B, Korelc J, Ibrahimbegović A (2002) Nonlinear shell problem formulation accounting for through-the-thickness stretching and its finite element implementation. Comput Struct 80(9-10):699-717

24. Cardoso RPR, Yoon JW (2005) One point quadrature shell element with through-thickness stretch. Comput Methods Appl Mech Eng 194(9):1161-1199

25. Brank B (2005) Nonlinear shell models with seven kinematic parameters. Comput Methods Appl Mech Eng 194(21-24):2336-2362

26. Sussman T, Bathe K-JR (2013) 3D-shell elements for structures in large strains. Comput Struct 122:2-12

27. Kulikov GM, Plotnikova SV (2008) Finite rotation geometrically exact four-node solid-shell element with seven displacement degrees of freedom. Computer Modeling in Engineering & Sciences 28(1):15-38

28. Kim D-N, Bathe K-J (2008) A 4-node 3D-shell element to model shell surface tractions and incompressible behavior. Comput Struct 86(21-22):2027-2041

29. Klinkel S, Gruttmann F, Wagner W (1999) A continuum based three-dimensional shell element for laminated structures. Comput Struct 71(1):43-62

30. Sansalone M, Sabourin F, Brunet M (2011) A new shell formulation using complete 3D constitutive laws. Int J Numer Methods Eng 86(6):688-716

31. Bassa B, Sabourin F, Brunet M (2012) A new nine-node solid-shell finite element using complete 3D constitutive laws. Int J Numer Methods Eng 92(7):589-636

32. 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

33. Bathe K-J, Dvorkin EN (1986) A formulation of general shell elements-the use of mixed interpolation of tensorial components. Int J Numer Methods Eng 22(3):697-722

34. Dvorkin EN, Pantuso D, Repetto EA (1995) A formulation of the mitc4 shell element for finite strain elasto-plastic analysis. Comput Methods Appl Mech Eng 125(1-4):17-40

35. 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

36. Reese S, Wriggers P, Reddy BD (2000) A new locking-free brick element technique for large deformation problems in elasticity. Comput Struct 75(3):291-304

37. Alves de Sousa RJ, Cardoso RPR, Fontes Valente RA, Yoon J-W, Grácio JJ, Natal Jorge RM (2006) A new one-point quadrature enhanced assumed strain (EAS) solid-shell element with multiple integration points along thickness-part ii: nonlinear applications. Int J Numer Methods Eng 67(2):160-188

38. Klinkel S, Wagner W (1997) A geometrical non-linear brick element based on the EAS-method. Int J Numer Methods Eng 40(24):4529-4545

39. Wriggers P, Korelc J (1996) On enhanced strain methods for small and finite deformations of solids. Comput Mech 18(6):413-428

40. Abed-Meraim F, Combescure A (2002) Shb8ps-a new adaptive, assumed-strain continuum mechanics shell element for impact analysis. Comput Struct 80(9):791-803

41. Abed-Meraim F, Combescure A (2009) An improved assumed strain solid-shell element formulation with physical stabilization for geometric non-linear applications and elastic-plastic stability analysis. Int J Numer Meth Eng 80(13):1640-1686

42. Reese S, Wriggers P (1997) A material model for rubber-like polymers exhibiting plastic deformation: computational aspects and a comparison with experimental results. Comput Methods Appl Mech Eng 148(3-4):279-298

43. Simo JC, Taylor RL (1985) Consistent tangent operators for rate-independent elastoplasticity. Comput Methods Appl Mech Eng 48(1):101-118

44. Miehe C, Apel N (2004) Anisotropic elastic-plastic analysis of shells at large strains. A comparison of multiplicative and additive approaches to enhanced finite element design and constitutive modelling. Int J Numer Methods Eng 61(12):2067-2113

45. Zienkiewicz OC, Taylor RL, Zienkiewicz OC, Taylor RL (1977) The finite element method, vol 3. McGraw-hill, London

46. 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

47. Liu WK, Guo Y, Tang S, Belytschko T (1998) A multiple-quadrature eight-node hexahedral finite element for large deformation elastoplastic analysis. Comput Methods Appl Mech Eng 154(1-2):69-132

48. Simo JC, Laursen TA (1992) An augmented Lagrangian treatment of contact problems involving friction. Comput Struct 42(1):97-116

49. 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

50. Alvesde Sousa RJ, Natal Jorge RM, FontesV alente RA, César de Sá JMA (2003) A new volumetric and shear locking-free 3D enhanced strain element. Eng Comput 20(7):896-925

51. Alves de Sousa RJ, Cardoso RPR, Fontes Valente RA, Yoon J-W, Grácio JJ, Natal Jorge RM (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 Meth Eng 62(7):952-977

52. Belytschko T, Tsay C-S (1983) A stabilization procedure for the quadrilateral plate element with one-point quadrature. Int J Numer Meth Eng 19(3):405-419

53. Belytschko T, Ong JS-J, Liu WK, Kennedy JM (1984) Hourglass control in linear and nonlinear problems. Comput Methods Appl Mech Eng 43(3):251-276

54. Belytschko T, Bindeman LP (1993) Assumed strain stabilization of the eight node hexahedral element. Comput Methods Appl Mech Eng 105(2):225-260
55. 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

56. Wempner G, Talaslidis D, Hwang C-M (1982) A simple and effi-
cient approximation of shells via finite quadrilateral elements. J
Appl Mech 49(1):115-120

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

58. Dvorkin EN (1995) Nonlinear analysis of shells using the mitc formulation. Arch Comput Methods Eng 2(2):1

59. Hauptmann R, Doll S, Harnau M, Schweizerhof K (2001) Solid-shell'elements with linear and quadratic shape functions at large deformations with nearly incompressible materials. Comput Struct 79(18):1671-1685

60. Electricité de France. Finite element *code_aster*, structures and thermomechanics analysis for studies and research. https://www.code-aster.org, 1989-2019

61. Simo JC, Rifai MS (1990) A class of mixed assumed strain methods and the method of incompatible modes. Int J Numer Methods Eng 29(8):1595-1638

62. 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

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

64. Hughes TJR (1980) Generalization of selective integration proce-
dures to anisotropic and nonlinear media. Int J Numer Methods
Eng 15(9):1413–1418

65. Hibbett K, Sorensen (1998) ABAQUS/standard: user's manual, vol
1. Hibbitt, Karlsson & Sorensen, Providence

66. 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

67. 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

68. Frisch-Fay R (1962) Flexible bars. Butterworths, New York

69. Saleeb AF, Chang TY, Graf W, Yingyeunyong S (1990) A hybrid/mixed model for non-linear shell analysis and its applications to large-rotation problems. Int J Numer Methods Eng 29(2):407-446

70. Sansour C, Bufler H (1992) An exact finite rotation shell theory, its mixed variational formulation and its finite element implementation. Int J Numer Methods Eng 34(1):73-115

71. Sze KY, Liu XH, Lo SH (2004) Popular benchmark problems for geometric nonlinear analysis of shells. Finite Elem Anal Des 40(11):1551-1569

72. Betsch P, Stein E (1999) Numerical implementation of multiplica-
tive elasto-plasticity into assumed strain elements with application
to shells at large strains. Comput Methods Appl Mech Eng 179(3-
4):215-245

73. Eberlein R, Wriggers P (1999) Finite element concepts for finite elastoplastic strains and isotropic stress response in shells: theoretical and computational analysis. Comput Methods Appl Mech Eng 171(3-4):243-279