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

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.


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 'solidshell', where modifications are introduced to enforce flexural behavior that permits stretching through a given thickness direction (see, for instance, References [1][2][3][4][5][6][7][8][9][10][11][12][13][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][16][17][18][19][20][21][22][23][24][25][26][27][28][29], among many others. The 5-parameter shell models have been enriched by a desired number of parameters at nodes to permit a representation 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 6parameter 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][16][17][18][20][21][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 a implicit code, specially in nonlinear problems, due to the three-parameters stabilization procedure. The shear strain adopted in the SB9γ 25 was interpolated following the work of Dvorking and Bathe [32][33][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 threedimensional and structural zones. All theses advantages have lead to develop a set of solid shell element in the last decades: [3,7,10,24,[35][36][37][38][39][40][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 Massud 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][53][54], Reese et al. [10,36] as well as Combescure et al. [7,41].
Poisson thickness locking is due to resulting incorrectconstant 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 a 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).

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 u ∈ V , the conjugated Green-Lagrange strain tensor E ∈ E and the second Piola-Kirchhoff stress S ∈ S are treated as independent variables. V , E and S being respectively the space of admissible displacement, strain and stress. The Hu-Washizu principle is written as follows: where W is the strain energy, F the deformation gradient depending on the displacement field u, I the metric tensor and π ext (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.
This means that the eight vertices nodes displacementdependent strain tensor E H of the element is enriched by an additional enhanced assumed strain E 9 ∈ E 9 thanks to the ninth node. The variation of the above functional is obtained from the directional derivative and leads to 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 S and 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 Kirchhoff stress S is chosen such that the corresponding Cauchy stress be linearly dependent to the applied pressure: ζ 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: Although E 9 = 0 for the continuum problem, in general E 9 h = 0, when we introduce finite element approximations. By denoting D the domain occupied by the body, we can observe that the space of enhanced strain field is in [L 2 (D)] 6 (see Simo [61]). Hence, no inter-element continuity on the E 9 need be enforced when construction finite element approximation. Note also that the enhanced strain interpolation E 9 h and the standard strain interpolation defined by S V h are independent in the sense that:

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 ξ , η, ζ ). 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 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: where I means node number. Auxiliary vectors are defined from previous combinations of the isoparametric vectors: These vectors form the basis for the interpolation field of the nodal displacements. The shape functions N I (ξ, η, ζ ) are obtained from a linear combination of the above-mentioned isoparametric vectors: The derivatives of the shape functions with respect to the isoparametric coordinates are also spanned by the isoparametric vectors in the form: 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 translational degrees of freedom: where x I and X I 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 where u I are the nodal displacement. Keep in mind that the effect of the ninth node will be discussed in a later section. The derivatives of displacements in the isoparametric coordinate system are interpolated from the shape functions derivatives as From Eq. (9) the Jacobian matrix can be decomposed into constant, linear and bi-linear components depending on ξ, η, ζ, ξ η, ζ η, ζ ξ 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".

Strain field
The deformation gradient is written as follows: F is a tensor which maps the reference basis G i to the current one g i (Fig. 2). e i is the Cartesian base vector. For simplicity, the superscript H will be omitted and E H will be written as E. Further the Green Lagrange strain tensor is represented by it's Cartesian and covariant components as follows E i j andÊ i j 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 and T a second-order matrix which contains the terms of the inverse of the jacobian J matrix as follows From the variation of the first line of Eq. (16) we can write the covariant strain-displacementB I at the node I (I = 1, . . . , 8), such that δÊ I =B I δU I as followŝ  (14),B I can be decomposed into the reduced integrated partB ri I and a stabilization part B stab I as follows: And The components of Eqs. (22) and (23) are the basic building blocks of all solid-shell elements and are detailed in "Appendix".

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 statical 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 statical 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 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Ê ηζ andÊ ξζ are interpolated using the equivalent shear deformation from four different tying points (Fig. 3)

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

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

The ninth node
The additional central node is endowed with only one translation degree of freedom in the element thickness direction ζ . The enhanced normal displacement in the thickness direction is written as follows u 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 H ζ the normal displacement of the eight node element considered alone. A new column is then added to the covariant matrixB H making it a 6 by 25 matrixB witĥ One can then write the variation of the pinching strain in the covariant basis as follows: B H 3 being the third line of matrixB H .

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 where δu + 9 and δ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 δu 9 the variation of the ninth node displacement [31]. Assuming the pinching stress S 33 is linear through the element thickness: the virtual equilibrium condition in the thickness direction can then be written and the equivalent nodal force distribution deduced as follows: 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.

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 refered 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 it's 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 it's 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 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 it's constant, linear and bilinear terms and the Cartesian Green Lagrange strain can then be derived as follows: and The strain-displacement gradient can be split the same way so that B can be written as E stab and B 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 Since no constant term is present in the expansion of B stab (see Eq. 37) The stabilization components of the strain and straindisplacement matrices being well identified, a similar approach is used to identify the stabilization stress.

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 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 ri (E ri ) is Piola Kirschhoff obtained from the integration of constitutive law and S stab (E 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 C stab being the deviatoric part of St Venant Kirchhoff material given as follows With μ the shear modulus evaluated, in elastic material behavior, as μ stab = E/2(1 + ν), in which E and ν 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: where These ways of evaluating the stabilization parameter μ 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 C stab ensures an efficient and robust stabilization.

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.
Equation (50) can be split into a material stiffness and a geometrical stiffness.
K ri m and K stab m are respectively the integrated stiffness matrix and its stabilization counterpart, K ri g and K stab g respectively the geometrical stiffness matrix and it's 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 it's stabilization counterpart as 8 9 (55)

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 overcomed simply with a Newton-Raphson algorithm without adopting arclength 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] -SHB8PS: Reduced integrated solid-shell with ANS by Abed Meraim and Combescure [40] -S4R: ABAQUS's four-node shell element models [65]

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×10 5 MPa and the Poisson ratio ν = 0.3. Two studies have been done. In the first study a uniform normal pressure P u = 0.01172 N/mm 2 is applied in the upper side of the plate. In the second study two normal pressures P u = 2 × 0.01172N/mm 2 and P l = 0.01172 N/mm 2 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 = 1mm 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-

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: 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 σ rr = [σ · n] · n = −p, where σ 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.

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 ν = 0.3. The total load is F max = 40N 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 ν = 0.0 and ν = 0.4999. The element still gives good result showing no sign of volumetric locking, see Fig. 8b.

Spherical shell with 18 degree hole
In this example the hemispherical shell with an 18 degree circular cut out at it's 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][70][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×10 7 N/mm 2 and the Poisson ratio ν = 0.3. The total number of increment to reach the   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.

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 × 10 6 N/mm 2 and the Poisson ratio ν = 0.3125 and the applied load in each side is F = 40K N. 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 20K N. This can be seen through the displacement reversal that occurs on the

Pinched cylindrical shell mounted over rigid diaphragms
This problem and its variations have been considered by many authors [3,[69][70][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×10 3 and the Poisson ratio ν = 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

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 = λ p 0 with p = 0.01N/mm 2 is applied. The maximum of the load factor is λ max = 40. The material is taken ideally elastic-plastic with E = 6.9 × 10 4 N/mm 2 , ν = 0.3 and the yield stress σ y = 248N/mm 2 . Figure 16 gives the evolution of the normal displacement of the plate center compared to the results given in [10].    Fig. 17 The elastic parameters are the Young modulus E = 200,000 MPa and the Poisson ratio ν = 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.

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).
J 1 , J 2 , J 3 being the rows of the Jacobian matrix in the current configuration. Note that the rows of the jacobian matrix contrain the coordinates of the covariant g i vectors.

Decomposition of T
The Inverse of the Jacobian matrix is decomposed keeping only the constant and linear terms as follows 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).
From Eq. (73) one can easily determine the linear Jacobian terms, as follows, all terms being known or determined easily.
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