Farid Abed-Meraim
email: farid.abed-meraim@metz.ensam.fr

Alain Combescure

An improved assumed strain solid-shell element formulation with physical stabilization for geometric non-linear applications and elastic-plastic stability analysis

Keywords: SHB8PS solid-shell, hourglass, shear and membrane locking, assumed strain method, orthogonal projection, elastic-plastic stability

In this paper, the earlier formulation of the SHB8PS finite element is revised in order to eliminate some persistent membrane and shear locking phenomena. This new formulation consists of a solid-shell element based on a purely three-dimensional approach. More specifically, the element has eight nodes, with displacements as the only degrees of freedom, as well as an arbitrary number of integration points, with a minimum number of two, distributed along the 'thickness' direction. The resulting derivation, which is computationally efficient, can then be used for the modeling of thin structures, while providing an accurate description of the various through-thickness phenomena. A reduced integration scheme is used to prevent some locking phenomena and to achieve an attractive, low-cost formulation. The spurious zero-energy modes due to this in-plane one-point quadrature are efficiently controlled using a physical stabilization procedure, whereas the strain components corresponding to locking modes are eliminated with a projection technique following the assumed strain method. In addition to the extended and detailed formulation presented in this paper, particular attention has been focused on providing full justification regarding the identification of hourglass modes in relation to rank deficiencies. Moreover, an attempt has been made to provide a sound foundation to the derivation of the co-rotational coordinate frame, on which the calculations of the stabilization stiffness matrix and internal load vector are based. Finally to assess the effectiveness and performance of this new formulation, a set of popular benchmark problems is investigated, involving geometric non-linear analyses as well as elastic-plastic stability issues.

INTRODUCTION

Large-scale finite element simulations are extensively used in engineering design and process control. In various fields of engineering (civil, aerospace, automotive, nuclear, and aerospace), three-dimensional non-linear problems still tax the resources of most computers. Indeed, meshes with sufficient resolution to achieve reasonable accuracy often require many hours of computer time, even with explicit methods. Therefore, the efficiency of finite elements is of crucial importance to speed up the design process and reduce the computational cost of these simulations. Over the past 20 years, considerable progress has been achieved in developing fast and reliable elements. In this regard, efficient (e.g. reduced integration) elements were developed for continuum mechanics problems [START_REF] Hughes | Generalization of selective integration procedures to anisotropic and nonlinear media[END_REF][START_REF] Flanagan | A uniform strain hexahedron and quadrilateral with orthogonal hourglass control[END_REF][START_REF] Belytschko | Hourglass control in linear and nonlinear problems[END_REF][START_REF] Liu | Finite element stabilization matrices-a unification approach[END_REF][START_REF] Liu | Multiple-quadrature underintegrated finite elements[END_REF][START_REF] Koh | New improved hourglass control for bilinear and trilinear elements in anisotropic linear elasticity[END_REF][START_REF] Belytschko | Assumed strain stabilization of the eight node hexahedral element[END_REF][START_REF] Hu | A one-point quadrature eight-node brick element with hourglass control[END_REF][START_REF] Masud | A stabilized 3D co-rotational formulation for geometrically nonlinear analysis of multi-layered composite shells[END_REF]. Concurrently, for structural mechanics problems, numerous efficient plate and shell elements have been developed based on mixed formulations or enhanced assumed strain (EAS) methods in order to avoid locking problems. Among these are Bathe and Dvorkin [START_REF] Bathe | Four-node plate bending element based on Mindlin/Reissner plate theory and a mixed interpolation[END_REF][START_REF] Bathe | Formulation of general shell elements-the use of mixed interpolation of tensorial components[END_REF], Onate and Castro [START_REF] Onate | Derivation of plate based on assumed shear strain fields[END_REF], Cheung and Chen [START_REF] Cheung | Refined hybrid method for plane isoparametric element using an orthogonal approach[END_REF], Ayad et al. [START_REF] Ayad | Un e ´le ´ment quadrangulaire de plaque base ´ sur une formulation mixte-hybride avec projection en cisaillement[END_REF][START_REF] Ayad | A new hybrid-mixed variational approach for Reissner-Mindlin plates: the MiSP model[END_REF], Chapelle and Bathe [START_REF] Chapelle | Fundamental considerations for the finite element analysis of shell structures[END_REF][START_REF] Chapelle | The Finite Element Analysis of Shells-Fundamentals[END_REF], Cardoso et al. [START_REF] Cardoso | Development of a one-point quadrature shell element for nonlinear applications with contact and anisotropy[END_REF], Fontes Valente et al. [START_REF] Valente | On the use of an enhanced transverse shear strain shell element for problems involving large rotations[END_REF][START_REF] Valente | Enhanced transverse shear strain shell formulation applied to large elasto-plastic deformation problems[END_REF], Gruttmann and Wagner [START_REF] Gruttmann | A stabilized one-point integrated quadrilateral Reissner-Mindlin plate element[END_REF], Cardoso and Yoon [START_REF] Cardoso | One-point quadrature shell elements for sheet metal forming analysis[END_REF][START_REF] Cardoso | One-point quadrature shell elements: a study on convergence and patch tests[END_REF], and Cardoso et al. [START_REF] Cardoso | A new approach to reduce membrane and transverse shear locking for one-point quadrature shell elements: linear formulation[END_REF].

However, in real-life structures, coexistence of three-dimensional and structural zones is quite common, and both types of elements must be used simultaneously. Elements that behave well in both continuum and structural applications considerably simplify the modeling of such structures, and avoid both arbitrary definitions of separation zones (e.g. continuum/structural) and the intricacies of connecting different types of elements (e.g. shell/continuum). Furthermore, continuum-based elements have many other advantages: the avoidance of complex shell-type kinematics, the use of general three-dimensional constitutive models, direct calculation of thickness (strain) variations, easy treatment of large rotations along with simple updating of configurations, straightforward connection with three-dimensional elements since displacements are the only degrees of freedom, and natural contact conditions on both sides of the structure.

As a consequence, much effort has been devoted to the development of solid-shell elements for use in finite element models of thin structures [START_REF] Domissy | Formulation et e ´valuation d'e ´le ´ments finis volumiques modifie ´s pour l'analyse line ´aire et non line ´aire des coques[END_REF][START_REF] Cho | Stability analysis using a geometrically nonlinear assumed strain solid shell element model[END_REF][START_REF] Hauptmann | A systematic development of solid-shell element formulations for linear and non-linear analyses employing only displacement degrees of freedom[END_REF][START_REF] Lemosse | ´ments finis isoparame ´triques tridimensionnels pour l'e ´tude des structures minces[END_REF][START_REF] Sze | A hybrid stress ANS solid-shell element and its generalization for smart structure modelling, part I-solid-shell element formulation[END_REF][START_REF] Hauptmann | Solid-shell elements with linear and quadratic shape functions at large deformations with nearly incompressible materials[END_REF][START_REF] Abed-Meraim | SHB8PS a new intelligent assumed strain continuum mechanics shell element for impact analysis on a rotating body[END_REF][START_REF] Abed-Meraim | SHB8PS-a new adaptive, assumed-strain continuum mechanics shell element for impact analysis[END_REF][START_REF] Legay | Elastoplastic stability analysis of shells using the physically stabilized finite element SHB8PS[END_REF][START_REF] Vu-Quoc | Optimal solid shells for non-linear analyses of multilayer composites[END_REF][START_REF] Chen | A mixed 8-node hexahedral element based on the Hu-Washizu principle and the field extrapolation technique[END_REF][START_REF] Kim | A resultant 8-node solid-shell element for geometrically nonlinear analysis[END_REF][START_REF] Alves De Sousa | A new one-point quadrature enhanced assumed strain (EAS) solid-shell element with multiple integration points along thickness: part I-geometrically linear applications[END_REF][START_REF] Alves De Sousa | A new one-point quadrature enhanced assumed strain (EAS) solid-shell element with multiple integration points along thickness: part II-nonlinear applications[END_REF][START_REF] Alves De Sousa | On the use of a reduced enhanced solid-shell (RESS) element for sheet forming simulations[END_REF][START_REF] Klinkel | A robust non-linear solid shell element based on a mixed variational formulation[END_REF][START_REF] Reese | A large deformation solid-shell concept based on reduced integration with hourglass stabilization[END_REF][START_REF] Cardoso | Enhanced assumed strain (EAS) and assumed natural strain (ANS) methods for one-point quadrature solid-shell elements[END_REF]. This growing interest has been motivated by the above-mentioned requirements and constraints that are common in many industrial applications. Most of the methods developed earlier were based on EAS fields, and consisted of either the use of a conventional integration scheme with appropriate control of all locking phenomena or the application of reduced integration with hourglass control. Both approaches have been extensively investigated and evaluated in various structural applications, as reported in the work of Dvorkin and Bathe [START_REF] Dvorkin | Continuum mechanics based four-node shell element for general non-linear analysis[END_REF], Belytschko and Bindeman [START_REF] Belytschko | Assumed strain stabilization of the eight node hexahedral element[END_REF], Zhu and Cescotto [START_REF] Zhu | Unified and mixed formulation of the 8-node hexahedral elements by assumed strain method[END_REF], Wriggers and Reese [START_REF] Wriggers | A note on enhanced strain methods for large deformations[END_REF], Klinkel and Wagner [START_REF] Klinkel | A geometrical non-linear brick element based on the EAS-method[END_REF], Klinkel et al. [START_REF] Klinkel | A continuum based three-dimensional shell element for laminated structures[END_REF], Wall et al. [START_REF] Wall | A deformation dependent stabilization technique, exemplified by EAS elements at large strains[END_REF], Reese et al. [START_REF] Reese | A new locking-free brick element technique for large deformation problems in elasticity[END_REF], Puso [START_REF] Puso | A highly efficient enhanced assumed strain physically stabilized hexahedral element[END_REF], Alves de Sousa et al. [START_REF] Alves De Sousa | A new volumetric and shear locking-free 3D enhanced strain element[END_REF], and Fontes Valente et al. [START_REF] Valente | An enhanced strain 3D element for large deformation elastoplastic thin-shell applications[END_REF].

Among the pioneering research dealing with thin structure modeling by means of threedimensional elements without rotational degrees of freedom, the work of Graf et al. [START_REF] Graf | On the numerical performance of three-dimensional thick shell elements using a hybrid/mixed formulation[END_REF] is notable for developing 8, 16, and 18-node three-dimensional elements based on hybrid/mixed formulations. Xu and Cai [START_REF] Xu | A new plate shell element of 16 nodes and 40 degrees of freedom by relative displacement method[END_REF] proposed a 16-node displacement-based isoparametric element with 40 degrees of freedom and plane-stress assumptions. Sze and Ghali [START_REF] Sze | Hybrid hexahedral element for solids, plates, shells and beams by selective scaling[END_REF] modified the 8-node hexahedral hybrid element first proposed by Pian and Tong [START_REF] Pian | Relations between incompatible model and hybrid stress model[END_REF] by introducing adjustable parameters in order to avoid excessively stiff behavior and to recover shell, plate, and beam solutions. Kim and Lee [START_REF] Kim | A solid element formulation for large deflection analysis of composite shell structures[END_REF] developed an 18-node hexahedral element for the analysis of large deflections of composite shell structures, in which the constitutive law was modified in order to uncouple the normal transverse

• The ability to model thin, three-dimensional structures using only a single layer of elements along the thickness, while accurately describing the various through-thickness phenomena (e.g. bending and elasto-plasticity). • Simplified meshing of complex structural forms, where shell and solid elements must coexist without any compatibility problems between different families of elements (continuum and structural elements for instance). • Easy treatment of large rotations and a straightforward procedure of updating configurations (with no rotational degrees of freedom involved) when compared with conventional shell element formulations. • Computational efficiency due to large admissible aspect ratios (allowing for optimal meshes), the use of reduced integration, and the elimination of shear and membrane locking by appropriate techniques. • A simple and attractive formulation (hexahedral geometry, eight nodes, only three translational degrees of freedom per node) thus avoiding complex and tedious pure-shell element formulations.

In this work, the formulation of the SHB8PS element is enhanced with new projections in order to eliminate some membrane and shear locking phenomena that were still present in the original formulation [START_REF] Abed-Meraim | SHB8PS-a new adaptive, assumed-strain continuum mechanics shell element for impact analysis[END_REF][START_REF] Legay | Elastoplastic stability analysis of shells using the physically stabilized finite element SHB8PS[END_REF]. Despite the geometry of the element (eight-node hexahedron with only displacement degrees of freedom), several modifications are introduced in order to incorporate shell features. Among them, a shell-like behavior for the element is achieved by modifying the three-dimensional constitutive law so that plane-stress conditions are approached and by aligning all of the integration points along a preferential direction, called the thickness.

The reduced integration scheme, initiated by the early contributions of Zienkiewicz et al. [START_REF] Zienkiewicz | Reduced integration technique in general analysis of plates and shells[END_REF], Hughes et al. [START_REF] Hughes | Reduced and selective integration techniques in finite element analysis of plates[END_REF], and Hughes [START_REF] Hughes | Generalization of selective integration procedures to anisotropic and nonlinear media[END_REF], is used in order to improve the computational efficiency and to alleviate some membrane and shear locking phenomena. The spurious zero-energy deformation modes due to this in-plane reduced integration are efficiently controlled by a stabilization technique following the approach given in Reference [START_REF] Belytschko | Assumed strain stabilization of the eight node hexahedral element[END_REF]. First, the corresponding hourglass modes are shown to be the vectors of the kernel of the stiffness matrix aside from the rigid body modes. To circumvent this stiffness matrix rank deficiency, the hourglass modes are explicitly derived using a basis of the vector space of the discretized displacements, and then are efficiently stabilized. It is worth noting that the proposed formulation is valid for any set of integration points located along the thickness direction and comprising at least two integration points.

In order to eliminate the various locking effects (transverse shear, membrane), the discrete gradient operator is projected onto an appropriate sub-space. This projection technique can be derived from the formalism of the assumed strain method based on the pioneering work of Simo and Rifai [START_REF] Simo | A class of mixed assumed strain methods and the method of incompatible modes[END_REF], which was successfully applied in subsequent contributions [START_REF] Simo | Geometrically non-linear enhanced strain mixed methods and the method of incompatible modes[END_REF][START_REF] Simo | Improved versions of assumed enhanced strain tri-linear elements for 3D finite deformation problems[END_REF][START_REF] Andelfinger | EAS-elements for 2D, 3D, plate and shell structures and their equivalence to HRelements[END_REF][START_REF] Kasper | A mixed-enhanced strain method: part I-linear problems[END_REF]. This approach is also shown to be justified within the framework of the Veubeke-Hu-Washizu mixed variational principle [START_REF] De Veubeke | Diffusion des inconnues hyperstatiques dans les voilures a ` longeron couple[END_REF]. It is well known that the procedure for choosing an assumed strain field is substantially complex since each term of the discrete gradient operator has to be handled separately in order to eliminate the components responsible for membrane and transverse shear locking.

The SHB8PS element was first developed within an explicit formulation and implemented into an explicit dynamic code (EUROPLEXUS) in order to simulate impact problems [START_REF] Abed-Meraim | SHB8PS a new intelligent assumed strain continuum mechanics shell element for impact analysis on a rotating body[END_REF][START_REF] Abed-Meraim | SHB8PS-a new adaptive, assumed-strain continuum mechanics shell element for impact analysis[END_REF]. This explicit version was also used to simulate bird ingestion by aircraft gas turbine engines, as well as other accidental situations suggested by the aeronautical company SNECMA. Next, an implicit version of the element was formulated and implemented into the quasi-static implicit code Stanlax-INCA for elastic-plastic stability applications [START_REF] Legay | Elastoplastic stability analysis of shells using the physically stabilized finite element SHB8PS[END_REF]. More recently, this version was implemented into the quasi-static implicit code ASTER, which was developed by the energy and electricity company EDF, due to its good performance in various applications.

In spite of the built-in projection aimed at eliminating locking phenomena, the former (explicit and implicit) formulations of SHB8PS showed a relatively slow convergence rate in the case of the pinched hemispherical test problem. The driving force behind the present development was the persistence of some locking modes in certain applications, as revealed in Reference [START_REF] Legay | Elastoplastic stability analysis of shells using the physically stabilized finite element SHB8PS[END_REF] and pointed out by Reese [START_REF] Reese | A large deformation solid-shell concept based on reduced integration with hourglass stabilization[END_REF][START_REF] Reese | On a physically stabilized one-point finite element formulation for three-dimensional finite elastoplasticity[END_REF]. This work focuses on the projection techniques in order to better eliminate various locking phenomena. The newly developed version of SHB8PS is presented in this paper and its very good performance is demonstrated. Through numerous well-known benchmark problems, this new formulation proves to be free of locking and exhibits good convergence toward analytical or numerical reference solutions.

The remainder of this paper is outlined as follows. In Section 2, the extended and detailed formulation is presented, including the new assumed strain fields. Through this formulation, particular attention has been focused on providing full justification regarding the identification and isolation of hourglass modes in relation to element rank deficiencies. Moreover, an attempt has been made to provide a sound foundation for the derivation of the co-rotational coordinate frame on which the calculation of the stabilization stiffness matrix and internal load vector is based. To assess the effectiveness and performance of this new formulation, Section 3 is devoted to numerical experiments involving geometric non-linear analyses as well as elastic-plastic stability problems. Finally, Section 4 presents some concluding remarks.

NEW IMPROVED FORMULATION OF THE SHB8PS ELEMENT

This section details the newly developed formulation of the SHB8PS element. Several features regarding the formulation (kinematics and interpolation) are common to the previously published

node 1 -1 -1 -1 2 1 -1 -1 3 1 1 -1 4 -1 1 -1 5 -1 -1 1 6 1 -1 1 7 1 1 1 n …
version [START_REF] Abed-Meraim | SHB8PS a new intelligent assumed strain continuum mechanics shell element for impact analysis on a rotating body[END_REF][START_REF] Abed-Meraim | SHB8PS-a new adaptive, assumed-strain continuum mechanics shell element for impact analysis[END_REF][START_REF] Legay | Elastoplastic stability analysis of shells using the physically stabilized finite element SHB8PS[END_REF]; the main modification here lies in the choice of the assumed strain field in relation to the associated orthogonal projection technique. The new assumed strain field has been specifically aimed at eliminating some residual shear and membrane locking phenomena. As a result, the stabilization fields (stiffness matrix and internal load vector), which differ from those of the previous formulation, will be accordingly derived.

Finite element interpolation

SHB8PS is a hexahedral, eight-node, and isoparametric element with linear interpolation. It is provided with a set of nint integration points spread along the 0 direction in the local coordinate frame. Figure 1 shows the reference geometry of the element, the nodal coordinates, as well as the location of its integration points. The coordinates xi , i = 1, 2, 3, of a point in the element are related to the nodal coordinates xiI using the classical linear isoparametric shape functions NI (I = 1, . . . , 8) and the relations:

8 xi = xiI NI (Ç, 1, 0) = }. xiI NI (Ç, 1, 0) ( 1
)
I =1
The convention of implied summation for repeated subscripts will be used hereafter, unless specified otherwise. The lowercase subscripts i vary from one to three and represent the directions of the spatial coordinates. The uppercase subscripts I vary from one to eight and correspond to the nodes of the element. With this convention, the interpolation of the displacement field ui inside the element in terms of the nodal displacements uiI is similar:

ui = uiI NI (Ç, 1, 0) (2)

Strain-displacement relation and discrete gradient operator

The displacement field interpolation, Equation (2), allows the strain field to be related to the nodal displacements. The linear part of the strain tensor is written as

εij = 1 (ui, j + u j,i ) = 1 (uiI NI, j + ujI NI,i ) (3) 2 2 5 8 6 7 int 2 1 1 4 2 3 8 -1 1 1
Figure 1. SHB8PS reference geometry, integration point location, and nodal coordinates. Then, the classical tri-linear shape functions for eight-node hexahedral elements are considered:

d T x T h T h T h T h T
NI (Ç, 1, 0) = 1 (1+Ç I Ç)(1+1 I 1)(1+0 I 0) 8 Ç, 1, 0 ∈ [-1, 1], I = 1,..., 8 (4)
These shape functions transform a unit cube in the reference space (Ç, 1, 0) into a general hexahedron in the (x1, x2, x3) space, as illustrated in Figure 2.

Combining Equations ( 1), [START_REF] Flanagan | A uniform strain hexahedron and quadrilateral with orthogonal hourglass control[END_REF], and (4) leads to the expansion of the displacement field as a constant term, linear terms in xi , and some terms depending on the ha functions:

ui = a0i +a1i x +a2i y +a3i z +c1i h1 +c2i h2 +c3i h3 +c4i h4, i = 1, 2, 3 h1 = 10, h2 = 0Ç, h3 = Ç1, h4 = Ç10 (5)
When this equation, for which the constants aji and cai will be subsequently defined, is evaluated at the element nodes, the following three eight-equation systems are obtained:

di = a0i s+a1i x1 +a2i x2 +a3i x3 +c1i h1 +c2i h2 +c3i h3 +c4i h4, i = 1, 2, 3 (6)
In the above equation, the di and xi vectors indicate the nodal displacements and coordinates, respectively, and are defined as

i = (ui 1, ui 2, ui 3, . . . , ui 8) i = (xi 1, xi 2, xi 3, . . . , xi 8)
The vectors s and ha (a = 1,..., 4) are given by

s T = (1, 1, 1, 1, 1, 1, 1, 1) 1 = (1, 1, -1, -1, -1, -1, 1, 1) 2 = (1, -1, -1, 1, -1, 1, 1, -1) 3 = (1, -1, 1, -1, 1, -1, 1, -1) 4 = (-1, 1, -1, 1, 1, -1, 1, -1) (7) (8) ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ dz
The unknown constants aji and cai given in Equations ( 5) and ( 6) are determined by introducing the bi (i = 1,..., 3) vectors from Hallquist [START_REF] Hallquist | Theoretical manual for DYNA3D[END_REF], defined as

aN bi = N,i (0) = ax i |Ç=1=0=0 , i = 1, 2, 3 Hallquist Form (9)
Explicit expressions for the derivatives of the shape functions evaluated at the origin of the (Ç, 1, 0) frame are derived in Appendix A together with some useful orthogonality properties leading to

aji = b T • di , cai = c T • di , i, j = 1, . . . , 3 with j a 1 3 T l ( 10
)
c a = 8 ha - }. (h a •x j )b j j =1 , a = 1, . . . , 4
This allows us to express the discrete gradient operator relating the strain field to the nodal displacements as

∇s (u) = B• d ⎡ b T T ⎤ x + ha,x c a 0 0 ⎢ T T ⎢ ⎢ 0 b y + ha,y c a 0 ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 0 0 b T + ha,z c T ⎢ (11) B = ⎢ ⎢ b T z a ⎢ T T T ⎢ ⎢ y + ha,y c a b x + ha,x c a 0 ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 0 b T + ha,z c T b T + ha,y c T ⎢ ⎢ z ⎣ b T T a y a ⎢ ⎦ T T where z + ha,z c a 0 b x + ha,x c a ⎡ u x ,x ⎤ ⎢ u y,y ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ u z,z ⎢ ⎡ d x ⎤ ∇s (u) = ⎢ ⎢ , d = ⎢ dy ⎢ (12) ⎢ ⎢ ⎢ ⎢ ⎢ u x ,y +u y,x ⎢ ⎣ ⎦ ⎢ ⎢ ⎢ ⎢ ⎢ u y,z +u z,y ⎢ ⎢ ⎢ ⎣ ⎦ u x ,z + u z,x
This form of the discrete gradient operator is useful since it allows each of the non-constant strain modes to be handled separately, so that an assumed strain field can be easily and conveniently built. Moreover, some orthogonality conditions involving the c a vectors which enter the expression of matrix B, are demonstrated in Appendix A. These properties will be useful in the subsequent hourglass stability analysis of the SHB8PS element. They will also help in choosing an appropriate assumed strain field and in evaluating the stabilization stiffness.

Hourglass analysis for the SHB8PS

The hourglass modes of the SHB8PS element are analyzed following the approach first introduced in Reference [START_REF] Belytschko | Assumed strain stabilization of the eight node hexahedral element[END_REF]. For the SHB8PS element, these spurious modes are shown to originate in the particular location of the integration points (along a line). They are characterized by a vanishing energy, while they should induce a non-zero strain. This singular behavior is explained by the difference between the kernel of the discrete and continuous stiffness operators. Recall that the shelllike behavior of the SHB8PS element is obtained by modifying its three-dimensional constitutive law to approach plane-stress conditions and by aligning the integration points of the element along a particular direction, called the thickness. This in-plane reduced integration also aims to increase the computational efficiency and to avoid some shear locking phenomena in bending-dominated problems. In a standard displacement-based formulation, the elastic stiffness is obtained using the integration points as follows:

r Ke = Oe n int B T •C•B dO = }. w(0 I ) J (0 I )B T (0 I )•C•B(0 I ) ( 13
)
I =1
where J (0 I ) is the Jacobian of the transformation between the unit reference configuration and the current configuration of an arbitrary hexahedron. It is important to underline that, although undertaken within a displacement-based approach, the investigation of hourglass modes conducted in this section and the associated stabilization procedure given in Section 2.4 are quite general as long as at least two integration points are used. This is further detailed in Appendix B where it is shown that the proposed stabilization still applies after the projection technique is implemented within the assumed strain framework, for which the underlying variational principle is presented in Section 2.5. Note also that two integration points are sufficient for both providing a rank sufficient element and dealing with elastic problems, as will be shown through the numerical examples given in Section 3. It has also been revealed, from illustrative test problems, that a minimum of five integration points should be used when dealing with elastic-plastic applications. Table I gives the coordinates and the associated weights of the Gauss points, which represent the roots of the Gauss-Legendre polynomial, in the case of five integration points along the thickness direction. For a set of nint integration points (I = 1,..., nint), with coordinates Ç I = 1 I = 0, 0 I 1= 0, the derivatives ha,i (a = 3, 4; i = 1, 2, 3) vanish. Consequently, for nint�2 operator B defined by Table I. Coordinates and weights of the Gauss points for five integration points along thickness. Equation [START_REF] Bathe | Formulation of general shell elements-the use of mixed interpolation of tensorial components[END_REF] reduces to B12, where the sum on the index a only goes from 1 to 2:

2 T x + ⎢ a 1 ⎢ ⎢ ⎢ ⎤ ha,x c a 0 0 ⎢ ⎢ 2 ⎢ ⎢ 0 b T + }. ha,y c T 0 ⎢ ⎢ a=1 ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ T }. ⎢ ⎢ B12 = ⎢ 0 0 b z + a=1 h a,z c a ⎢ ⎢ ⎢ (14) ⎢ 2 2 ⎢ ⎢ b T }. }. ⎢ ⎢ y + ⎢ ⎢ ⎢ ⎢ a=1 h a,y c T b T + a=1 ha,x c T 0 ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 0 b T + }. ⎢ a=1,2 ⎢ ⎢ 2 h a,z T b T + }. h a=1 a,y c T ⎢ a ⎢ ⎢ ⎢ ⎢ ⎣ T z + a=1 h a,z c a 0 b x + }. a=1,2 h a,x c T ⎦
In order to identify the kernel of the stiffness matrix, a basis for the vector space of the discretized displacements is built. Then, the reduced integration is shown to diminish the rank of the discrete stiffness. Indeed, according to Equation ( 13), the rank of the stiffness matrix Ke is closely related to that of the B matrix. In other words, the zero-strain modes d that verify at each integration point the equation following should be found:

∇s (u) = B(0 I )• d = 0 (15)
A detailed analysis of hourglass modes is given in Appendix B; hereafter, only the main results are reported. Using expression [START_REF] Ayad | Un e ´le ´ment quadrangulaire de plaque base ´ sur une formulation mixte-hybride avec projection en cisaillement[END_REF] for the discrete gradient operator computed at the integration points and making use of the orthogonality relations (A13) and (A16), see Appendix A, the kernel of the stiffness matrix can be explicitly derived. This naturally reveals six rigid body modes of which consists the kernel of a fully integrated stiffness:

⎛ s ⎞ ⎛ 0 ⎞ ⎛ 0 ⎞ ⎛ y ⎞ ⎠ ⎛ z ⎞ , ⎛ 0 ⎞ z (16) ⎠ ⎝ ⎠ , ⎝ ⎠ , ⎝ ⎠ , ⎝-x , ⎝ 0 ⎠ ⎝ ⎢ 0 ⎢ 0 ⎢ s ⎢ 0 ⎢ 0 ⎢ s ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 0 -x -y
The first three column vectors correspond to the translations along the Ox, Oy and Oz axes, respectively. The three remaining vectors refer to the rotations about the Oz, Oy and Ox axes, respectively. For nint�2, in addition to these six rigid body modes, the following six vectors are also found in the kernel of the stiffness matrix Ke: The hourglass modes corresponding to the Ox axis are shown in Figure 3 for a hexahedron with a single integration point, nint = 1, located at the origin of the reference frame. Similar modes are obtained for the Oy and Oz axes by axis permutation. Unlike the one-point quadrature hexahedron (see Reference [START_REF] Belytschko | Assumed strain stabilization of the eight node hexahedral element[END_REF]) comprising 12 hourglass modes as shown in Figure 3, only 6 hourglass modes are found for the SHB8PS provided that at least two integration points are considered. They are composed of the h3 and h4 vectors as expressed in Equation [START_REF] Chapelle | The Finite Element Analysis of Shells-Fundamentals[END_REF] and illustrated in Figure 4.

⎛ h 3 ⎞ ⎢ 0 ⎢ ⎛ 0 ⎞ ⎢ h ⎢ ⎛ 0 ⎞ ⎢ 0 ⎢ ⎛ h 4 ⎞ ⎢ 0 ⎢ ⎛ 0 ⎞ ⎢ h ⎢ ⎛ 0 ⎞ ⎢ 0 ⎢ ⎝ ⎠ , ⎝ 0 3⎠ , ⎝ 0 ⎠ , ⎝ h3 ⎠ , ⎝ 0 4⎠ , ⎝ 0 (17) ⎠ h4

Stabilization of spurious zero-energy modes

The control of the six hourglass modes of the SHB8PS element, as revealed by Equation [START_REF] Chapelle | The Finite Element Analysis of Shells-Fundamentals[END_REF], is achieved by adding a stabilization stiffness to the stiffness matrix Ke. This is drawn from the approach of Reference [START_REF] Belytschko | Assumed strain stabilization of the eight node hexahedral element[END_REF], in which an efficient stabilization technique was applied along with an assumed strain method for the eight-node hexahedral element with uniform reduced integration. The stabilization forces are deduced in the same way. It is important to note that this stabilization

⎢ ⎢ ⎢ 4 ⎢ a T ⎢ 4 ⎢ a
part is treated completely independently of the assumed strain projection part, since the latter is intended to eliminate the locking phenomena. This projection technique will be applied in Section 2.5. The starting point consists in decomposing the discrete gradient operator B into two parts as follows:

B = B12 +B34 (18)
The first term in this additive decomposition is given by Equation ( 14). The second term B34 is precisely the one that vanishes at the integration points, and is given by the following matrix form:

⎡ 4 ⎤ }. h a,x c T ⎢ a 3 ⎢ ⎢ ⎢ 0 ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 4 ⎢ }. h a,y c T ⎢ a=3 ⎢ ⎢ }. ⎢ ⎢ 0 0 ⎢ ⎢ ⎢ 4 4 a=3 h a,z c a ⎢ ⎢ ⎢ ⎢ (19) ⎢ }. h a,y c T ⎢a=3 ⎢ ⎢ }. a=3 ha,x c a 0 ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 0 }. ⎢ a=3,4 ⎢ ⎢ 4 h a,z T }. h a=3 a,y c T ⎢ a ⎢ ⎢ ⎢ ⎢ ⎣ }. h a,z c T h a,x c T ⎦ a=3 a=3,4
In the standard displacement approach, the stiffness matrix and the internal forces are defined as r Ke =

f int = B T • C• B dO Oe r B T • r dO Oe (20)
By introducing the additive decomposition [START_REF] Cardoso | Development of a one-point quadrature shell element for nonlinear applications with contact and anisotropy[END_REF] of the B operator, the stiffness matrix becomes

Ke = r r B T •C•B12 dO+ Oe Oe r B T •C•B34 dO+ Oe r B T •C•B12 dO+ Oe B T •C•B34 dO ( 21
)
which can be simply written as

Ke = K12 +KSTAB (22)
The first term, K12, is the only one taken into account when the stiffness is evaluated at the integration points as defined previously: The second term, KSTAB, represents the stabilization stiffness since it vanishes if evaluated at the integration points:

r K 12 = n int B T •C•B12 dO = }. w(0 I ) J (0 I )B T (0 I )•C•B12(0 I ) ( 23
)
r 12 • C•B Oe r dO+ Oe r B T •C•B12 dO+ Oe B T •C•B34 dO (24)
In a similar way, the internal forces of the element can be written as

f int = f int STAB ( 25
)
The first term, f int , is the only one taken into account when the forces are evaluated at the integration points:

int

r 12 = Oe B T • r dO = n int }. I =1 w(0 I ) J (0 I )B T (0 I )• r(0 I ) ( 26
)
The second term f STAB of Equation ( 25) represents the stabilization forces and should be consistently calculated according to the stabilization stiffness given by Equation ( 24), see Reference [START_REF] Belytschko | Assumed strain stabilization of the eight node hexahedral element[END_REF].

Since the stabilization stiffness matrix and internal load vector cannot be calculated properly at the integration points, we will calculate them in the co-rotational coordinate system proposed in Reference [START_REF] Belytschko | Assumed strain stabilization of the eight node hexahedral element[END_REF] in order to prevent the hourglass mode phenomena in case of geometric nonlinearities. Some justification for the derivation of this co-rotational frame is also provided in Appendix C. An intermediate stage of this approach consists in projecting B onto a B matrix in order to eliminate the remaining locking problems.

Hu-Washizu variational principle and assumed strain field

The discrete gradient operator is projected onto an appropriate sub-space in order to eliminate shear and membrane locking. This projection technique can be derived from the formalism of the assumed strain method. This approach can also be justified within the framework of the Hu-Washizu non-linear mixed variational principle (see for instance Korelc and Wriggers [START_REF] Korelc | Consistent gradient formulation for a stable enhanced strain method for large deformations[END_REF]). Indeed, this three-field variational principle reads

r n(v, ε ˙, r ¯ ) = Oe T bε ˙ r •r dO+b Oe r ¯ T •(∇s (v)-ε ˙) dO-bd ˙ T •f ext = 0 ( 27
)
where b denotes a variation, v the velocity field, ε ˙ the assumed strain rate, r the interpolated stress, r the stress evaluated by the constitutive law, d ˙ the nodal velocities, f ext the external nodal forces and ∇s (v) the symmetric part of the velocity gradient. The assumed strain formulation used to construct the SHB8PS element is a simplified form of the Hu-Washizu variational principle as described by Simo and Hughes [START_REF] Simo | On the variational foundations of assumed strain methods[END_REF]. In this simplified form, the interpolated stress is chosen to be orthogonal to the difference between the symmetric part of the velocity gradient and the assumed strain rate. Consequently, the second term of Equation ( 27) vanishes, yielding r

T n(ε ˙) = bε ˙ Oe •r dO-bd ˙ T •f ext = 0 (28)
In this form, the variational principle is independent of the stress interpolation, since the interpolated stress is eliminated and no longer needs to be defined. The discrete equations then only require

O y a B ˆ = y ⎢ ⎢ ˆ ⎢ T x T ⎢
the interpolation of the velocity and of the assumed strain field. The assumed strain rate ε ˙ is expressed in terms of a B matrix, projected starting from the classical discrete gradient B defined by Equation ( 11):

ε ˙(x , t ) = B ¯ (x )• d ˙ (t ) (29)
Once this expression is substituted into the variational principle [START_REF] Lemosse | ´ments finis isoparame ´triques tridimensionnels pour l'e ´tude des structures minces[END_REF], new expressions for the elastic stiffness and internal forces are obtained:

r T r T Ke = B •C•B dO, f int = Oe B •r(ε ˙) dO ( 30
)
Oe Before defining the projected B operator, let us replace in the previous equations the Hallquist form of the bi vectors, Equation ( 9), with the mean form b ˆ i from Flanagan and Belytschko [START_REF] Flanagan | A uniform strain hexahedron and quadrilateral with orthogonal hourglass control[END_REF]:

1 b ˆ i = e r N,i (Ç, 1, 0) dO, i = 1, 2, 3 (31)
Oe Accordingly, the vectors c a are replaced by the vectors c ˆa defined as

1 3 T l c ˆa = 8 ha - }. (h a • x j )b ˆ j j =1 (32)
Finally, matrix B, defined by Equation [START_REF] Bathe | Formulation of general shell elements-the use of mixed interpolation of tensorial components[END_REF], is replaced by the B ˆ operator, defined as

x + ha,x c ˆT 0 0 ⎡ b ˆ T ⎢ 0 ⎢ ⎢ a b ˆ T + h a,y ⎤ c T 0 ⎢ ⎢ ⎢ ⎢ 0 0 b ˆ T + ha,z c ˆT ⎢ ⎢ ⎢ ⎢ b ˆ T ⎢ ⎢ + h a,y c ˆa b ˆ T + h a,x c ˆa z a ⎢ ⎢ 0 ⎢ ⎢ ⎢ (33)
⎢ z T y T ⎢ ⎢ 0 b ˆ T + ha,z c ˆa ⎣ b ˆ T + ha,y c ˆa ⎢ ⎦ b ˆ T T T T z + ha,z c ˆa 0 b ˆ x + ha,x c ˆa
The approach developed earlier still applies, as well as the expressions of the stabilization stiffness and internal forces, as long as the same additive decomposition is adopted:

B ˆ = B ˆ 12 +B ˆ 34 ( 34
)
It is noteworthy that in the former version of the SHB8PS element, the Hallquist forms bi were only replaced by the mean expressions b ˆ i of Flanagan-Belytschko in the stabilization terms B ˆ 34 and thus in KSTAB.

It is also important to note that both forms bi and b ˆ i have been tested on a large number of test problems and that Flanagan-Belytschko's mean form performed better in all cases. The better convergence of this latter form is most clear when few, highly distorted elements are used. Similar results have been reported in Reference [START_REF] Belytschko | Assumed strain stabilization of the eight node hexahedral element[END_REF] with an assumed strain, eight-node solid element with one-point quadrature.

a 0 0 ⎢ a ⎢ 4 ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎢
At this stage, operator B ˆ in Equation ( 34) can be projected onto a B ˆ operator such that

B ˆ = B ˆ 12 +B ˆ 34 (35)
Only the second term B ˆ 34 from Equation ( 34) is projected; the first term B ˆ 12 remains unchanged and is given by Equation ( 14) where vectors bi are replaced by b ˆ i . The operator B ˆ 34 is projected onto B ˆ 34, given by

⎡ 4 ⎤ }. h a,x c ˆT ⎢ a=3 ⎢ ⎢ ⎢ 0 ⎢ B ˆ 34 = ⎢ 4 }. a=3 ⎢ ⎢ ⎢ ha,y c ˆT 0 ⎢ ⎢ T ⎢ (36) ⎢ 0 0 h3,z c ˆ3 ⎢ ⎢ ⎢ ⎢ 0 0 0 ⎢ ⎢ 0 0 0 ⎢ ⎦ 0 0 h4,x c ˆT
The elastic stiffness is then given by Equation ( 22) as the sum of the following two contributions:

r K 12 = n int B ˆ T • C• B ˆ 12 dO = }. w(0 I ) J (0 I )B ˆ T (0 I )• C• B ˆ 12(0I ) (37) 12 12
Oe

I =1 r r T r T K STAB = B ˆ T ˆ 34 ˆ 34 12 ˆ 34 ˆ 34 12 • C•B Oe dO+ B Oe •C•B ˆ dO+ B Oe • C•B dO (38)
The stabilization stiffness, Equation [START_REF] Alves De Sousa | A new one-point quadrature enhanced assumed strain (EAS) solid-shell element with multiple integration points along thickness: part II-nonlinear applications[END_REF], is calculated in a co-rotational coordinate system given in Reference [START_REF] Belytschko | Assumed strain stabilization of the eight node hexahedral element[END_REF]. This orthogonal co-rotational system, which is embedded in the element and rotates with it, is chosen to be aligned with the referential coordinate system (see Figure 5). This choice

F z y x z S R y x Figure 5
. Schematic representation of the co-rotational coordinate system.

K T K T 1 a 1 ax i i a T a2 =
is justified here by the rotation extracted from the polar decomposition of the deformation gradient as reported in Reference [START_REF] Abed-Meraim | Stabilite ´ des e ´le ´ments finis sous-inte ´gre[END_REF] and discussed in Appendix C. As noted in Reference [START_REF] Belytschko | Assumed strain stabilization of the eight node hexahedral element[END_REF], such a co-rotational approach has numerous advantages, including simplified expressions for the above stabilization stiffness matrix, whose first two terms vanish, and a more effective treatment of shear locking in this frame. In addition, the co-rotational system assures a frame-invariant element.

The main equations defining the adopted co-rotational coordinate system are given in the following. First, the components of the column vectors forming the rotation matrix are computed:

a1i = K T • xi , a2i = K T • xi , i = 1, 2, 3 (39) 1 2
with

K T 1 = (-1, 1, 1, -1, -1, 1, 1, -1) 2 = (-1, -1, 1, 1, -1, -1, 1, 1) 3 = (-1, -1, -1, -1, 1, 1, 1, 1) (40)
Then, the correction term ac is calculated so that the orthogonality relation a T •(a2 +ac) = 0 is verified:

T ac = -1 • a1 (41)
1 •a1 The third base vector a3 is then obtained by the cross-product:

a3 = a1 ×(a2 +ac) (42)
The rotation matrix R that maps a vector in the global coordinate system to the co-rotational system is finally given, after normalization, by

R1i = a 1i , R2i = ±a1± a2i + aci ±a2 +ac± , R3i = a 3i , i = 1, 2, 3 (43

) ±a3±

The stabilization terms (stabilization stiffness and internal forces; Equation [START_REF] Alves De Sousa | A new one-point quadrature enhanced assumed strain (EAS) solid-shell element with multiple integration points along thickness: part II-nonlinear applications[END_REF]) are computed in this co-rotational coordinate system, where several terms can be simplified. Because this corotational coordinate system is chosen to be aligned with the reference frame, the relationship between the two coordinate systems can be approximated as

ax ˜i 1 K T x , ax ˜i aÇ j 0 if i j (44) aÇ i = aÇ i /ax ˜i = 8 i • ˜i aÇ j = ax ˜i = 1=
in which vector x ˜i denotes the nodal coordinates expressed in the co-rotational system, and repeated subscripts do not indicate a summation. Equation ( 44) allows the following simplifications to be made:

ahi 8Ç k 8Ç j Ç k K T •x ˜1 K T • x ˜2 K T • x ˜3 h i,i = ˜i = 0, h j,i = K T • x ˜i , h4,i = , K T • x ˜i J ˜ 1 2 8 8 3 (45) 8
where J ˜ denotes the determinant of the Jacobian matrix. Note also that in these last formulas as well as in the subsequent equation, there is no sum on repeated subscripts; moreover, subscripts

¯ ¯ +2 0 0 0 0 0 E 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 i 0 C = ⎢ 0 3 ⎢ i , j
, and k are two by two distinct and take values 1, 2, and 3 with all of the possible permutations. Then, Equation ( 45) simply leads to r hi, j dO = 0

Oe

r H ii = Oe (h j,i ) 2 dO = r (hk,i ) 2 dO = 3 Oe r (h4,i ) 2 dO Oe (46) 1 (K T •x ˜ j )(K T •x ˜k ) j k = 3 (K T •x ˜i ) r 1 T Hij = hi, j h j,i dO = K k • x ˜k
Oe Using these explicit expressions, the stabilization stiffness given in Equation ( 38) is obtained completely analytically in this co-rotational system as

⎡ k 11 k 12 k 13 ⎤ K STAB = ⎢ k 21 k 22 k 23 ⎢ ( 47
)
where the 8×8 matrices kij are given by

⎣ ⎦ k 31 k 32 k 33 k11 = ( ¯ +2 )H11[c ˆ3c ˆT + 1 c ˆ4c ˆT] 3 3 4 k22 = ( ¯ +2 )H22[c ˆ3c ˆT + 1 c ˆ4c ˆT] k 33 = H 11 1 T 3 3 4 ( 48
) 3 c ˆ4c ˆ4 kij = 0, i 1= j
In a linear analysis with perfectly rectangular elements, the above third block matrix k33 should be rather taken as k33 = H11[cc ˆ3c ˆT + 1 c ˆ4c ˆT], with c = 0.01. Note also that an improved, plane-stress-type constitutive law is adopted for the SHB8PS element, in order to enhance its immunity with regard to thickness locking. This specific law, which uncouples the response in terms of in-plane and transverse normal stress versus normal strain is given by

⎡ +2 ¯ ⎢ ⎢ ⎢ ⎢ 0 0 0 0 ⎤ ⎢ ⎢ ⎢ 0 ⎢ E E v ⎢ ⎢ ⎢ ⎢ = , ¯ = (49) ⎢ 0 ⎢ ⎢ ⎢ 2(1+v) 1-v 2 ⎢ ⎢ ⎢ ⎢ ⎣ ⎦ f f f a
where E is Young's modulus and v is Poisson's ratio. Note that as usually adopted, these material properties are specified with respect to a local physical coordinate system, in which the x -y plane corresponds to the element mid-plane defined by the 0-coordinate of the considered integration point. To illustrate this, a possible choice of such a local physical coordinate system, such as adopted here, is described in Appendix D. The choice of this constitutive matrix avoids the locking encountered with a full three-dimensional law. Moreover, in contrast to the commonly adopted plane-stress assumption, this modified stiffness matrix allows the deformation energy associated with the strains normal to the mean surface of the element to be taken into consideration.

For the computation of the internal forces of the element, the same approach is adopted (see also References [START_REF] Abed-Meraim | SHB8PS a new intelligent assumed strain continuum mechanics shell element for impact analysis on a rotating body[END_REF][START_REF] Abed-Meraim | SHB8PS-a new adaptive, assumed-strain continuum mechanics shell element for impact analysis[END_REF]). The additive decomposition [START_REF] Chen | A mixed 8-node hexahedral element based on the Hu-Washizu principle and the field extrapolation technique[END_REF] and the projection [START_REF] Kim | A resultant 8-node solid-shell element for geometrically nonlinear analysis[END_REF] allow us to calculate the stabilization forces:

n int f int = }. w(0 I ) J (0 I )B ˆ T (0 I )• r(0 I )+ f STAB ( 50
)
where

I =1 ⎛ STAB ⎞ 1 f STAB = ⎢ f STAB ⎢ (51)
and

⎢ 2 ⎢ ⎝ ⎠ STAB 3 STAB i = 4 }. a=3 Qi ac ˆa, i = 1, 2, 3 (52)
Qi a, called the generalized stresses and entering the expressions of the stabilization forces, are related to the so-called generalized strains qi a by the following incremental equations:

Q ˙ 13 = ( ¯ +2 )H11q ˙13 Q ˙ 14 = 1 11 14 3 ( ¯ +2 )H q˙ Q ˙ 23 = ( ¯ +2 )H22q ˙23 (53) Q ˙ 24 = 1 22 24 3 ( ¯ +2 )H q˙ Q ˙ 33 = 0 Q ˙ 34 = 1 11 34
3 H q ˙ The generalized strain rates q ˙i a are given by

q ˙i a = c ˆT •d ˙ i , i = 1, 2, 3, a = 3, 4 (54)
Once the stabilization terms, i.e. stabilization stiffness and internal load vector, are computed in the co-rotational coordinate system through Equations ( 47)-( 48) and ( 51)- [START_REF] Xu | A new plate shell element of 16 nodes and 40 degrees of freedom by relative displacement method[END_REF], respectively, they have to be transformed back to the global coordinate system. Note also that the previous expressions for the stabilization stiffness and forces hold for elastic behavior. In the case of elastic-plastic behavior, Young's modulus E is replaced by the mean tangent modulus (i.e. the average of the tangent moduli at the integration points across the thickness). This choice avoids an overstiff response that would correspond to a purely elastic hourglass stabilization scheme. Moreover, this strategy results in an adaptive element provided with a stabilization technique that automatically adjusts to the physical situation of the element, whether elastic or elastic-plastic.

Numerical implementation for non-linear analyses

In this section, the main features of the implementation of the SHB8PS element are briefly described. For this purpose, the incremental, non-linear, and implicit finite element code Stanlax-INCA has been used. In this process, the updated Lagrangian strategy is adopted. For the stress and internal variable updates, the well-known co-rotational formulation given in Reference [START_REF] Liu | A multiple-quadrature eight-node hexahedral finite element for large deformation elastoplastic analysis[END_REF] is used. 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. For structural instability problems involving either a load-limit point ('snap-through') or a deflection-limit point ('snap-back'), as well as for material instability (softening behavior), the path-following Riks algorithm [START_REF] Riks | An incremental approach to the solution of snapping and buckling problems[END_REF], which is based on an arc-length control parameter, is adopted.

It is worth noting that the proposed formulation is based on simple ideas, which makes it easy to incorporate into implicit and explicit non-linear programs, but in turn it requires the use of smaller load steps than more sophisticated solid-shell elements based on mixed methods (see, e.g. Klinkel et al. [START_REF] Klinkel | A robust non-linear solid shell element based on a mixed variational formulation[END_REF]). On the other hand, this simplicity is such that the CPU time is only 20% higher for the SHB8PS element than for the standard DKT shell finite element for instance. This is also attributable, in part, to the fast computation of the stabilization stiffness matrix, for which the running time is about the same as that of the assumed strain stabilization of Belytschko and Bindeman [START_REF] Belytschko | Assumed strain stabilization of the eight node hexahedral element[END_REF], which was shown to be nearly as fast as the perturbation-type hourglass control of Flanagan and Belytschko [START_REF] Flanagan | A uniform strain hexahedron and quadrilateral with orthogonal hourglass control[END_REF]. Indeed, the use of the co-rotational coordinate system allows great simplifications to be made, and makes it possible to integrate the stabilization stiffness matrix in closed form so that numerical integration is not required.

For coupling with non-linear behavior models, an elastic-plastic constitutive law with isotropic hardening and an associated plastic flow rule has been used. As previously mentioned, the standard three-dimensional elastic constitutive law has been specifically modified for this element formulation, and this must accordingly be taken into account for the time integration of the set of constitutive equations. This is the main modification with respect to the classical radial return mapping algorithm based on Newton-Raphson's iterative procedure. Note also that this first choice of a relatively simple non-linear behavior model has been adopted for the sake of simplicity; more advanced, physically-based elastic-plastic models can be used in the same way.

The associated yield criterion is defined by

F = aeq -ay (ε p )�0 ( 55
)
where aeq is the von Mises equivalent stress and ay is the yield stress, which can be described by a non-linear function of the equivalent plastic strain ε p . For isotropic hardening, Equation ( 55) can be regarded as a geometric transformation for the yield surface, in which this surface, whose current size is ay , expands homogenously without distortion in stress space.

NUMERICAL EXAMPLES

In order to validate the new version of the SHB8PS element, its performance has been assessed based on the analysis of a variety of benchmark problems frequently used in the literature. For each test problem, the results were compared with the reference solutions and with those given by the earlier version of the SHB8PS element [START_REF] Legay | Elastoplastic stability analysis of shells using the physically stabilized finite element SHB8PS[END_REF]. Note that several projections have been formulated in this study and extensively tested over a wide range of benchmark problems. The projection presented here is the one that showed the best accuracy and convergence rate and exhibited no transverse shear or membrane locking phenomena. This projection had better results than the former version of the SHB8PS element in all situations, especially in the test of the pinched hemispherical shell, where the improvement is particularly significant. The results given hereafter include this above-mentioned test case as well as a set of representative popular benchmark problems commonly used to test finite element performance. While the first two numerical examples are linear elastic, all the remaining benchmark problems are non-linear, involving geometric or material non-linear computations as well as stability analyses. For elastic problems, two Gauss points have been considered along the thickness direction. Recall that this choice of a minimum number of two integration points is primarily dictated by rank deficiency considerations. Hence, the default number of two integration points is recommended for the proposed solid-shell element formulation in elastic problems as no effect of increasing this number has been experienced. The third example problem is specifically designed to address the choice of the adequate number of through-thickness integration points in the context of plasticity. It is revealed that a minimum of five integration points should be used when dealing with elastic-plastic applications. The last numerical test, in which a non-linear elastic-plastic buckling analysis is undertaken, is also specifically intended to evaluate the SHB8PS element in the context of plastic instability.

Pinched hemispherical shell

This test problem, which is often used to assess the three-dimensional inextensional bending behavior of shells, has become very popular and has been adopted by many authors since it was proposed by MacNeal and Harder [START_REF] Macneal | A proposed standard set of problems to test finite element accuracy[END_REF]. This test is severe because the transverse shear and membrane locking phenomena are predominant and are further accentuated by the particular geometry of the problem (distorted, skewed elements). This problem was studied in detail by Belytschko et al. [START_REF] Belytschko | Assumed strain stabilization procedure for the 9-node Lagrange shell element[END_REF], who showed that since all the elements are incurved, the intensity of membrane and shear locking is increased. They also showed that in this doubly curved shell problem, membrane locking is much more severe than shear locking. Figure 6 shows the geometry, loading, and boundary conditions for this elastic thin shell problem (R/t = 250). The radius is R = 10, the thickness is t = 0.04, Young's modulus is E = 6.825×10 7 , and the Poisson ratio is v = 0.3. Owing to the symmetry of the problem (i.e. planes (XZ) and (YZ)), only one quarter of the hemisphere is meshed using a single layer of elements through the thickness and with two unit loads along the directions Ox and Oy. Except for the symmetry, the boundary conditions are free; nevertheless, the displacement of one point in the z-direction is fixed in order to prevent rigid body motions. According to the reference solution [START_REF] Macneal | A proposed standard set of problems to test finite element accuracy[END_REF], the displacement of point A along the x--direction is equal to 0.0924 (see Figure 6). The convergence results are reported in Table II in terms of normalized displacement at point A in the x -direction versus the number of elements. The new version of the SHB8PS element is compared with the former version and with the three elements HEX8, HEXDS, and H8-ct-cp. The HEX8 element is the standard, eight-node, full integration solid element (eight integration points).

The HEXDS element is an eight-node, four-point quadrature solid element (see Liu et al. [START_REF] Liu | A multiple-quadrature eight-node hexahedral finite element for large deformation elastoplastic analysis[END_REF]).

The H8-ct-cp element is formulated and described in Reference [START_REF] Lemosse | ´ments finis isoparame ´triques tridimensionnels pour l'e ´tude des structures minces[END_REF]. Table II demonstrates that the new version of the SHB8PS element provides excellent convergence and shows no locking. This represents a significant improvement on the poor convergence showed by the previous version, as revealed in Reference [START_REF] Legay | Elastoplastic stability analysis of shells using the physically stabilized finite element SHB8PS[END_REF] and pointed out in [START_REF] Reese | A large deformation solid-shell concept based on reduced integration with hourglass stabilization[END_REF][START_REF] Reese | On a physically stabilized one-point finite element formulation for three-dimensional finite elastoplasticity[END_REF].

Patch tests

The patch tests were introduced by MacNeal and Harder [START_REF] Macneal | A proposed standard set of problems to test finite element accuracy[END_REF] in order to assess the ability of newly developed finite elements to reproduce constant stress states with meshes using few distorted elements. For plate and shell elements, two patch tests have been defined for investigating the membrane and the out-of-plane bending behavior, respectively. These consist of a patch of five elements subjected to two different deformation states by prescribing particular displacements at the exterior nodes of the mesh. These patch tests have been investigated for solid-shell formulations in several recent papers (see e.g. Vu-Quoc and Tan [START_REF] Vu-Quoc | Optimal solid shells for non-linear analyses of multilayer composites[END_REF], Klinkel et al. [START_REF] Klinkel | A robust non-linear solid shell element based on a mixed variational formulation[END_REF], Reese [START_REF] Reese | A large deformation solid-shell concept based on reduced integration with hourglass stabilization[END_REF], and Cardoso et al. [START_REF] Cardoso | Enhanced assumed strain (EAS) and assumed natural strain (ANS) methods for one-point quadrature solid-shell elements[END_REF]). In both loading situations, we consider the five-element mesh described in Figure 7.

In this figure, the geometric and material parameters are specified along with the coordinates of Table III. Membrane patch test: displacements of the interior nodes obtained with the current SHB8PS formulation.

E = 1×10 6 = 0.25 L = 0.24 l = 0.
Node

u 1 u 2 1 5.00×10 -5 4.00×10 -5 2 1.95×10 -4 1.20×10 -4 3 2.00×10 -4 1.60×10 -4 4 1.20×10 -4 1.20×10 -4
the four top interior nodes; the coordinates of the eight exterior nodes being simply deduced from the geometric data.

Membrane patch test.

In this first loading case, a pure membrane strain state is investigated that must result in a constant membrane stress field. To this end, the following displacements:

u1 = 10 -3 (x1 + 1 x2), u2 = 10 -3 ( 1 x1 + x2) ( 56
) 2 2
are prescribed at the eight exterior nodes (four on the bottom layer and four on the top layer). The u3 displacement of the four bottom exterior nodes is set to zero. The numerical results obtained with the current SHB8PS formulation in terms of displacements at the interior nodes are reported in Table III. These confirm with the analytical solution given by Equation [START_REF] Pian | Relations between incompatible model and hybrid stress model[END_REF]. In addition, a plane-stress state is obtained with a constant in-plane membrane stress field:

a11 = a22 = 1333, a12 = 400 (57)
which is in agreement with the theoretical solution given in Reference [START_REF] Macneal | A proposed standard set of problems to test finite element accuracy[END_REF]. This shows that the proposed element formulation passes the membrane patch test exactly.

Out-of-plane bending patch test.

In this second loading case, a bending strain state is investigated for which the stress analytical solution also reveals a constant stress field. For this purpose, the following displacements:

10 -3 t ( 1 \ 10 -3 t ( 1 \ 10 -3 2 2 u1 = ± 2 x1 + 2 x2 , u2 = ± 2 x1 + x2 , u3 = 2 (x 1 + x1 x2 + x 2 ) ( 58
)
are applied to the bottom and top exterior nodes, with t the thickness of the plate (see Figure 7). The numerical results obtained with the proposed formulation are compared with the analytical solution (see Reference [START_REF] Macneal | A proposed standard set of problems to test finite element accuracy[END_REF]). The latter analytical solution consists of the displacements of the interior nodes, which satisfy Equation ( 58), and the stress values at the top and bottom surfaces given by

a11 = a22 = ±0.667, a12 = ±0.200 (59)
It is well known that passing the out-of-plane bending patch test is more difficult than the membrane counterpart (see the analysis conducted in Reese [START_REF] Reese | A large deformation solid-shell concept based on reduced integration with hourglass stabilization[END_REF]). A detailed discussion on this issue can also be found in Vu-Quoc and Tan [START_REF] Vu-Quoc | Optimal solid shells for non-linear analyses of multilayer composites[END_REF], who reported that without the use of the Assumed Natural Strain (ANS) method for the transverse shear strains originated by Dvorkin and Bathe [START_REF] Dvorkin | Continuum mechanics based four-node shell element for general non-linear analysis[END_REF], it turns out to be difficult to fulfill the analytical solution for the bending patch test. They combined both the EAS and the ANS methods and proposed a minimal number of EAS parameters that is required to satisfy the out-of-plane bending patch test. What is remarkable in their investigation is that even with such a combination, an element formulation may fail to pass the bending patch test if a smaller number of EAS parameters is employed.

Following the same analysis as in Reference [START_REF] Reese | A large deformation solid-shell concept based on reduced integration with hourglass stabilization[END_REF], each of the five original elements is successively subdivided into 1 ×1, 2×2, 4×4 ... sub-elements until convergence is reached. In these successive discretizations, we use the mesh nomenclature 1 ×1 (original patch), 2×2, 4×4, and 8×8 elements. The numerical results are normalized with respect to the analytical solution and reported in Table IV for the displacements and in Table V for the stresses. For comparison purposes, the results yielded by the HEX8 and the Q1SPs elements are also reported when available. The latter represents the solid-shell formulation proposed in Reference [START_REF] Reese | A large deformation solid-shell concept based on reduced integration with hourglass stabilization[END_REF]. Note that in order to analyze the stress state on the top and bottom surfaces of the plate, the trapezoidal rule with three integration points along the thickness has been used. As can be observed in Tables IV andV, although the original five-element patch does not fulfill the analytical bending solution, it converges toward the analytical solution with mesh refinement in the same way as the Q1SPs solid-shell element proposed in [START_REF] Reese | A large deformation solid-shell concept based on reduced integration with hourglass stabilization[END_REF]. The performance of the proposed formulation on the original bending patch test is quite predictable since the element does not incorporate any EAS or ANS fields to enrich its transverse shear strain variation. Indeed, recall that the transverse shear strain is constant through the thickness within the current formulation, while it has a linear through-thickness variation in the bending patch test.

Elastic-plastic bending and number of through-thickness integration points

As previously discussed, two integration points are sufficient for linear elasticity applications. However, when elastic-plastic constitutive behavior is considered, the non-linear stress variation through the thickness requires more integration points for an accurate response. In addition, in order to correctly track the progress of the elastic-plastic front, the constitutive equations need to be evaluated at several points within the element. This issue has been discussed in several ), which generally revealed that at least five to seven integration points are required to capture the non-linear effects characteristic of elasto-plasticity.

u 1 /

Pure bending of an elastic-perfectly plastic beam.

Pure bending is investigated here to illustrate the need for more than two integration points when plasticity occurs. This simple illustrative problem, for which an analytical solution is available, is very suitable to emphasize the non-linear effects induced by plasticity in contrast to elasticity. The main objective of this preliminary example, for which no finite element solution is required, is to demonstrate the impact of the number of through-thickness integration points on the accuracy of the results. Figure 8 provides the material data for this problem as well as the computed responses in terms of moment versus curvature for different numbers of Gauss points through the thickness. This clearly confirms that while the elastic part of the loading curve is exactly described, two or three Gauss points fail to predict the entire curve properly, when plasticity occurs. These results suggest that a minimum number of five integration points are required to approach the analytical solution.

Elastic-plastic bending of a cantilever beam.

The trends revealed by the previous example can be confirmed here using the proposed solid-shell element. In this test problem, an elastic-plastic cantilever beam is investigated. Figure 9(a) shows the geometry, loading, and boundary conditions for this thin beam problem (L/ h = 100). The material properties, elastic-plastic behavior with linear isotropic hardening of constant parameter H , are those of References [START_REF] Eberlein | Finite element concepts for finite elastoplastic strains and isotropic stress response in shells: theoretical and computational analysis[END_REF] and [START_REF] Tan | Efficient and accurate multilayer solid-shell element: nonlinear materials at finite strain[END_REF], which serve as reference solutions for the current computations. These reference solutions were obtained with a mesh of 20×1×1 shell, respectively, solid-shell elements using five Gauss points through the thickness, and the former reference solution has also been validated with respect to numerical results by Dvorkin et al. [START_REF] Dvorkin | A formulation of the MITC4 shell element for finite strain elasto-plastic analysis[END_REF]. The load-deflection results given in Figure 9(b) correspond to a mesh of 10×1×1 SHB8PS elements. As can be observed in this figure, the onset of plasticity appears at a load of about P = 5 and the calculations using five Gauss points through the thickness are in good agreement with the reference results. When only two or three integration points are considered, it is clearly shown, once again, that the corresponding load-deflection curves are far away from the reference solution.

In summary, throughout the paper, two integration points are adopted for all elasticity applications, while five integration points are used for the cases involving plasticity. It should be noted, however, that the proposed formulation is quite general and the number of through-thickness integration points can be increased whenever required for particular applications.

Slit annular plate subjected to lifting line force

The following series of benchmark problems involve geometric non-linear effects (large rotations and displacements). The first test problem in this category was originally considered by Basar and Ding [START_REF] Basar | Finite rotation elements for the nonlinear analysis of thin shell structures[END_REF][START_REF] Basar | Finite-rotation shell elements for the analysis of finite-rotation shell problems[END_REF], and has since been adopted by many other authors (see References [START_REF] Domissy | Formulation et e ´valuation d'e ´le ´ments finis volumiques modifie ´s pour l'analyse line ´aire et non line ´aire des coques[END_REF][START_REF] Lemosse | ´ments finis isoparame ´triques tridimensionnels pour l'e ´tude des structures minces[END_REF][START_REF] Sze | Popular benchmark problems for geometric nonlinear analysis of shells[END_REF][START_REF] Buechter | Shell theory versus degeneration-a comparison in large rotation finite element analysis[END_REF][START_REF] Wriggers | Thin shells with finite rotations formulated in Biot stress: theory and finite element formulation[END_REF][START_REF] Brank | On implementation of a nonlinear four node shell finite element for thin multilayered elastic shells[END_REF][START_REF] Barut | Nonlinear analysis of laminates through a Mindlin-type shear deformable shallow shell element[END_REF][START_REF] Sansour | On hybrid stress, hybrid strain and enhanced strain finite element formulations for a geometrically exact shell theory with drilling degrees of freedom[END_REF][START_REF] Sansour | Families of 4-node and 9-node finite elements for a finite deformation shell theory. An assessment of hybrid stress, hybrid strain and enhanced strain elements[END_REF][START_REF] Hong | An assumed strain triangular curved solid shell element formulation for analysis for plates and shells undergoing finite rotations[END_REF][START_REF] Sze | An eight-node hybrid-stress solid-shell element for geometric nonlinear analysis of elastic shells[END_REF][START_REF] Kim | Curved quadratic triangular degenerated-and solid-shell elements for geometric nonlinear analysis[END_REF][100]). The main interest of this problem is that it is well suited for testing shell elements under significant finite rotations. This elastic thin slit annular plate (Ri /t = 200) is shown in Figure 10(a), which provides its geometric and material parameters as well. The initially circular ring has a slit cut along the radial direction A-B, at which a vertical line force P is applied along its free edge, while the Different meshes, including 6×48, 8×64, and 10×80, all of which had a single layer of elements along the thickness, were tested to investigate the convergence. Very accurate reference results were tabulated by Sze et al. [START_REF] Sze | Popular benchmark problems for geometric nonlinear analysis of shells[END_REF] using the ABAQUS shell element S4R with two different meshes, 6×30 and 10×80. Therefore, this reference solution was taken for comparison, and a similar intermediate mesh with 8×64 SHB8PS elements was chosen. The results are given in Figure 10(c) in terms of normalized load versus vertical displacement at the tips A and B, which reveal very good agreement with the reference solution.

Pull-out of an open-ended cylindrical shell

This test problem consists of an elastic thin cylindrical shell with free edges subjected to a pair of diametrically opposite radial forces (Figure 11). These particular boundary conditions allow the cylinder to undergo significantly large rotations, combining bending and membrane effects and thus providing a severe test for finite element formulations. For this reason, this benchmark has become very popular, as indicated by the number of studies that have used it [START_REF] Masud | A stabilized 3D co-rotational formulation for geometrically nonlinear analysis of multi-layered composite shells[END_REF][START_REF] Sze | Popular benchmark problems for geometric nonlinear analysis of shells[END_REF][START_REF] Brank | On implementation of a nonlinear four node shell finite element for thin multilayered elastic shells[END_REF][START_REF] Sansour | On hybrid stress, hybrid strain and enhanced strain finite element formulations for a geometrically exact shell theory with drilling degrees of freedom[END_REF][START_REF] Sansour | Families of 4-node and 9-node finite elements for a finite deformation shell theory. An assessment of hybrid stress, hybrid strain and enhanced strain elements[END_REF][START_REF] Hong | An assumed strain triangular curved solid shell element formulation for analysis for plates and shells undergoing finite rotations[END_REF][START_REF] Sze | An eight-node hybrid-stress solid-shell element for geometric nonlinear analysis of elastic shells[END_REF][START_REF] Kim | Curved quadratic triangular degenerated-and solid-shell elements for geometric nonlinear analysis[END_REF][100][101][102][103][104][105]. The geometric and material properties, as well as the boundary conditions and loading, are described in Figure 11(a). Only one octant of the cylinder is modeled due to the symmetry of the geometry, loading, and boundary conditions. In order to investigate the convergence, several meshes were tested, including 16×24, 20×30, and 24×36, with only a single element along the thickness. Again, Sze et al. [START_REF] Sze | Popular benchmark problems for geometric nonlinear analysis of shells[END_REF] provided, in tabulated form, highly accurate reference results for this test, using the ABAQUS shell element S4R with converged meshes, namely 16 ×24 and 24×36 elements. This reference solution was used for comparison, and an intermediate mesh of 20×30 SHB8PS elements (20 elements along the longitudinal direction and 30 along the circumferential direction) was considered. Figure 11(c) shows the results in terms of normalized load versus radial displacements at points A, B, and C. Point A corresponds to the point under loading, while points B and C are on the side of the cylinder and undergo horizontal displacements, with point C also on the free edge. As revealed by the load-displacement curves, the overall response exhibits two regimes: a preliminary stage dominated by bending effects and characterized by large displacements and rotations, and a later phase dominated by membrane effects, which may cause locking. Another interesting feature of this benchmark problem is associated with a snap-through phenomenon exhibited when the loading reaches a critical value of around 20×10 3 . This can clearly be seen through the displacement reversal that occurs on the load-displacement curve of point C, which is also visible on the deformed mesh under loading (see Figure 11(b)). The results of the current SHB8PS formulation are in good agreement with the reference solution, and no locking effects are experienced.

P/Pmax E = 21×10 6 = 0 R i = 6,

Twisted beam under out-of-plane loading

The geometric and material parameters for this test problem, as well as the boundary conditions and loading, are shown in Figure 12(a). This thin elastic beam (L/t = 375), twisted by an angle of 90 • between its two ends, is clamped at its right edge and subjected to an out-of-plane loading at its free left edge. The linear version of this benchmark problem has been used extensively to test the performance of finite elements in the context of warped configurations. This test is now considered as a reference shell test, since the initial element distortion increases the severity of the test with regard to locking. The linear convergence of SHB8PS has been previously demonstrated in References [START_REF] Abed-Meraim | SHB8PS-a new adaptive, assumed-strain continuum mechanics shell element for impact analysis[END_REF][START_REF] Legay | Elastoplastic stability analysis of shells using the physically stabilized finite element SHB8PS[END_REF], for both in-plane and out-of-plane loadings. The non-linear case under consideration has been investigated much less frequently in the literature. Different meshes were tested in order to achieve convergence, including 12×2, 24×4, and 48×4, with a single element in the thickness direction. The current SHB8PS results, obtained with a mesh of 24×4 elements, are shown in Figure 12(c), while Figure 12(b) displays the deformed configurations under successive loading states. The reference solution given by Smolenski [106] employing shell elements was used for comparison. Figure 12(c) illustrates this comparison by providing the plots of normalized load versus displacement at the tip in the three directions, which prove to be in good agreement. [106] and current SHB8PS results.

Lateral buckling of a cantilever beam

This test problem is illustrated below, with the geometric parameters, material properties, boundary conditions, and loading all specified in Figure 13(a). The lateral buckling of an elastic cantilever beam is considered. As shown in Figure 13(a), one side is clamped, while an in-plane loading is applied to the free edge. The beam undergoes large rotations, and after a linear pre-buckling regime along the fundamental path, a bifurcation point is detected, resulting in a bifurcated branch with symmetric stable bifurcation. A regular 10×2 mesh, with only a single layer of elements along the thickness, was used as in Reference [107], in which an EAS shell formulation was employed in conjunction with an efficient and accurate technique to pinpoint bifurcation points and to track the associated bifurcation branches. For this reason, those results were taken as reference solutions for the sake of comparison. The deformed configurations, characteristic to the lateral post-buckling deformation, are shown in Figure 13(b) for several successive loading states. Figure 13(c) gives the plots of the load-deflection curves in terms of load versus vertical displacement at the tip. Both the linear part of the fundamental path (pre-buckling) and the non-linear bifurcated branch are represented. The Euler critical load is predicted quite well by this current SHB8PS formulation. Likewise, the non-linear post-buckling branch (bifurcated path) is also correctly predicted. Through this example, the present finite element formulation proves its capability of predicting such types of instability by detecting singular points and the associated post-buckling behavior. The geometric parameters, material properties, boundary conditions, and loading are all described in this figure. Owing to the snap-through behavior exhibited by this test problem, it has been frequently used in the literature with different types of boundary conditions. The clamped-hinged version studied here is deemed more interesting and challenging, since its asymmetric boundary conditions allow non-symmetric buckling to occur. In order to follow the load-deflection curves beyond the limit points, the path-following Riks method based on an arc-length control parameter is adopted. The entire post-buckling response was provided in Reference [START_REF] Klinkel | A geometrical non-linear brick element based on the EAS-method[END_REF] using a robust strategy. Because of the lack of symmetry, the arch is modeled here entirely by means of a 40×1 mesh with only a single element along the thickness. In the literature, similar mesh densities have been adopted; in Reference [108] a 16×1 mesh with eight-node shell elements was employed, while Reference [START_REF] Domissy | Formulation et e ´valuation d'e ´le ´ments finis volumiques modifie ´s pour l'analyse line ´aire et non line ´aire des coques[END_REF], taken as the reference solution, used a 32×1 mesh with eight-node solid elements specifically modified for the efficient analysis of shell structures. Figure 14(c) shows the results obtained in terms of load versus displacements of the load point in both the x and the y directions. Figure 14(b) shows the evolution of the shape of the arch under different loading levels.

Clamped-hinged deep circular arch subjected to a concentrated load

The current SHB8PS calculations match the reference solution very well, confirming once again the ability of the present formulation to reproduce such unstable snap-through behavior.

Elastic-plastic buckling of a stiffened cylindrical ring

In this last benchmark problem, the SHB8PS element is again tested in the framework of elasticplastic behavior. This problem also allows the verification of the formulation of the geometric stiffness matrix Kr (see Appendix E), as well as that associated with the follower pressure K p . Furthermore, this problem emphasizes the ability of the present element formulation to take geometric imperfections into account. The test consists of a portion of a submarine hull subjected to external pressure, as illustrated in Figure 15(a). This test was previously used in Reference [109], while in Reference [110] the submarine was modeled in its entirety. The following analysis only considers a single ring, and is described in Figure 15(a) along with its geometric and material data. Note that due to its symmetry, only one quarter of the ring is modeled and subjected to the corresponding symmetry boundary conditions (see Figure 15). Table VI provides the stress-strain curve defining the uniaxial tensile curve for the adopted elastic-plastic model with non-linear isotropic hardening. The calculations take into account the follower pressure and are carried out in two stages. In the first step, a linear analysis of buckling is performed to determine the Euler critical pressure, based on the eigenvalue analysis of the global stiffness matrix. The second stage consists in performing a fully non-linear analysis after slightly altering the structure by adding a small geometric imperfection to the initial shape along the first Euler buckling mode (first linear eigenmode revealed by the preliminary Euler buckling analysis).

Linear buckling analysis (first Euler's critical load).

In this linear analysis, the Euler critical pressure is determined as well as the corresponding buckling mode. This critical state is associated with the lowest pressure that makes the global stiffness matrix singular, and is classically obtained by solving the following eigenvalue problem:

(Ke + c (Kr +K p ))•Xc = 0 (60)
in which c is the critical buckling load and Xc is the associated buckling mode. The first critical load obtained with the current SHB8PS calculations was 5.49 MPa; the associated buckling mode is shown in Figure 15(b). This value is in very good agreement with the value of 5.52 MPa provided by the COMU element. The formulation of this latter element is based on an axisymmetric geometry together with the Fourier series decomposition for the displacement [111].

3.9.2. Non-linear elastic-plastic buckling analysis. The fully non-linear buckling analysis accounts for both the elastic-plastic constitutive law, described in Table VI, and the geometric imperfection. This imperfection was embedded in the initial mesh, which used 440 elements, in the form of a small amplitude of the Euler eigenmode (in this example, the amplitude was 0.1 times the thickness). The path-following Riks strategy was used to carry out the calculations, and the obtained results are shown in Figure 15(c) in terms of normalized pressure versus radial displacement of point A. As can be expected, the non-linear constitutive model in conjunction with the geometric imperfection decreases the elastic critical load, which is now about 60% of its initial value. These results are in good agreement with those of Reference [109]; the latter results were also validated by comparison with results performed with the COMU finite element.

CONCLUSIONS

The earlier formulation of the SHB8PS solid-shell element has been revised with regard to locking, and a new, improved version has been developed and implemented into the implicit, non-linear finite element code Stanlax-INCA. This resulting derivation has been assessed, based on the analysis of a variety of popular benchmark problems frequently used in the literature. The key idea of this development is the adequate combination of a reduced integration rule with the wellknown assumed strain method. Another interesting feature of this approach is the convenient fully three-dimensional framework on which this solid-shell element is based (eight-node hexahedron with only three translational degrees of freedom per node).

In order to achieve an attractive, low-cost formulation, the computational efficiency of the element has been enhanced by adopting an in-plane one-point quadrature scheme. The resulting spurious zero-energy modes are controlled using an effective and physical stabilization procedure. For the out-of-plane integration, a set consisting of an arbitrary number of integration points, with a minimum number of two, has been chosen along a particular direction, designated as the 'thickness'. These choices, along with the use of a modified elastic constitutive law specifically aimed at alleviating thickness locking, all contribute to the accurate analysis of bending-dominated structural problems using only a single layer of elements along the thickness. Moreover, the projection technique adopted in the current formulation eliminates the various locking phenomena much better. Indeed, the excellent accuracy and convergence properties of the element have been clearly demonstrated through numerous linear and non-linear benchmark problems. All of these tests revealed that no residual locking (membrane, shear) has been experienced. In particular, the

1 • ⎢ aNI |0
improvement is shown to be significant on the pinched hemispherical shell problem, in which the amount of locking observed in the earlier formulation was completely eliminated. Several non-linear structural analysis test problems have been treated as well, involving geometric non-linear effects and/or snap-through or bifurcation-type instabilities. The SHB8PS element performed well in all of these examples and proved its capability for predicting such types of instability by detecting singular points and the associated post-buckling behavior. The coupling with elasto-plasticity, which thus far has used only a relatively simple constitutive law based on isotropic hardening, confirms the benefits of such a solid-shell formulation in dealing with nonlinear elastic-plastic stability problems for shell structures. This development should be pursued further, since it shows promise in prospective applications to sheet metal forming processes, which would be of great interest.

APPENDIX A: ORTHOGONALITY CONDITIONS AND VECTORS bi

In this appendix, some orthogonality relations necessary for the derivation of the discrete gradient operator, expressed in Section 2.2, are demonstrated. These properties are also needed for the analysis of the element rank deficiency (see Section 2.3), as well as for the identification of hourglass modes, which is carried out in Appendix B. The first step consists of providing explicit expressions for the derivatives of the shape functions, evaluated at the origin of the reference coordinate system, as given by Equation ( 9). This particular form of vectors bi , referred to as Each of the eight components of the above vector can be expressed as

⎛ aN I ⎞ ⎢ aÇ ⎢ ⎢ ⎢ a NI 1 ( a NI aÇ a NI a1 a NI a0 \ ( aÇ a1 a0 \ ⎢ aN ⎢ 1 ⎢ ⎢ axi 1 0 = aÇ axi + a1 axi + a0 axi = axi axi axi ⎢ ⎢ |0 ⎢ ⎢ I ⎢ a1 ⎢ ⎢ (A2) ⎢ ⎢ ⎝ ⎠ a0 |0
From Equation (4), the derivatives aNI /an, where n T = (Ç 1 0), are straightforward:

⎛ Ç I (1+1 I 1)(1+0 I 0) ⎞ a NI 1 ⎢ ⎢ an = 8 ⎢ 1 I (1+Ç I Ç)(1+0 I 0) ⎢ (A3) ⎝ ⎠ 0 I (1+Ç I Ç)(1+1 I 1) K 8 K ⎢ ⎢ ⎢ ⎢ ⎠ b ⎢ ⎢ ⎢ ⎢ ⎢ ⎠ = = j which gives ( aN \ T ⎛ -1 1 1 -1 -1 1 1 -1 ⎞ 1 1 = ⎢ -1 -1 1 1 -1 -1 1 1 ⎢ = ⎛ T ⎞ 1 ⎢ K T ⎢ (A4) an |0 8 ⎝ ⎠ -1 -1 -1 -1 1 1 1 1 ⎢ 2 ⎢ ⎝ ⎠ T 3
For the remaining terms, aÇ j /axi , further derivations are needed. This requires the calculation of matrix J such that Once this matrix J is known, the Hallquist form of vectors bi is simply obtained by

⎛ b T ⎞ 1 ⎢ b T ⎢ ( aN \T ⎢ 2 ⎢ = J|0 • ⎝ ⎠ T 3 an |0 (A6)
Using the Jacobian matrix and its inverse, matrix J can be rewritten as Making use of Equation ( 1), the element coordinate interpolation can be expanded as

xi = 1 (s T +ÇK T +1K T +0K T + h1h T + h2h T + h3h T + h4h T )•xi , i = 1, 2, 3 (A8) 8 1 2 3 1 2 3 4
The components of the Jacobian matrix are then obtained by differentiation: axi 1 T T which gives

F ij = aÇ = 8 (K j + h a,Ç j h a )• xi , i, j = 1, . . . , 3, a = 1, . . . , 4 (A9) ⎛ K T T T ⎞ 1 •x1 K 2 •x1 K 3 •x1 1 T T T F|0 = 8 ⎢ K 1 •x2 K 2 • x2 K 3 •x2 ⎢ (A10) ⎢ ⎢ ⎝ ⎠ K T T T 1 •x3 K 2 •x3 K 3 •x3 1 1 -1 -1 1 1 -1 1 1 -1 -1 1 |0 2 2 2 j a 8 and hence ⎛ K T T T ⎞ -1 1 •x1 K 1 •x2 K 1 •x3 T -1 ⎢ T T T ⎢ J|0 = (F ) = 8 ⎢ K •x1 K • x2 K • x3 ⎢ (A11) ⎝ ⎠ K T T T 3 •x1 K 3 •x2 K 3 •x3 The Hallquist form of vectors bi is finally given by ⎛ b T ⎞ ⎛ K T T T ⎞ -1 1 1 •x1 K 1 •x2 K 1 •x3 ⎛ -1 -1 ⎞ ⎢ b T ⎢ ⎢ K T T T ⎢ ⎢ ⎢ ⎢ 2 ⎢ = ⎢ 2 •x1 K 2 •x2 K 2 •x3 ⎢ •⎝-1 1 (A12) ⎠ ⎝ ⎠ ⎝ ⎠ b T K T T T 3 3 •x1 K 3 •x2 K 3 •x3 -1 -1 -1 -1 1 1 1 1
Using this expression (A12) of vectors bi , most of the following orthogonality conditions become straightforward by simple algebra:

b T T T i •ha = 0, b i •s = 0, b i •x j = bij h T T a •s = 0, h a • h fJ = 8b afJ i, j = 1,..., 3
•x j = N T (A13) •x j = x j,i 0 = b . Finally, i ,i |0 | ij
the constants aji and cai given in Equation ( 10) are obtained by multiplying Equation ( 6) by b T and h T , respectively, and by using the above orthogonality relations:

1 3 l aji = b T •di , cai = c T •di with c = ha - }. (h T • x j )b j (A14) j a a 8 a j =1
Although not required, constant a0i can be deduced by multiplying Equation ( 6) by s T :

1 3 l a0i = b T •di with b = s- }. (s T • x j )b j j =1 (A15)
Moreover, it is easy to show, using Equation (A13), that the c a vectors that are defined in Equation (A14) verify the following useful orthogonality conditions:

c T T a •x j = 0, c a •h fJ = b afJ (A16)
Replacing the constants a0i , aji, and cai in Equation ( 5), the displacement field can be expressed in the following convenient form: ui

= (b T + x1b T + x2b T + x3b T + h1c T + h2c T + h3c T + h4c T )•di (A17) 1 2 3 1 2 3 4
Differentiating this last equation with respect to x j yields the displacement gradient as

1 4 j + h \ c T •di = (b T + ha, j c T )•di (A18) a j a a=1 ⎢ 0 ⎢ , e = ⎢ h ⎢ , ⎝ ⎠ 14 ⎝ 1⎠ e = ⎢ 0 ⎢ , e = ⎢ 0 ⎢ , e = ⎢ h ⎢ , e = ⎢ 0 ⎢ 0 0 h1 0 0 h2 ⎛ h 3 ⎞ ⎛ 0 ⎞ ⎛ 0 ⎞ ⎛ h 4 ⎞ ⎛ 0 ⎞ ⎛ 0 ⎞ i i i i i i

APPENDIX B: IDENTIFICATION OF HOURGLASS MODES

As mentioned in Section 2.3, the hourglass patterns of the element correspond to zero-energy modes, i.e. eigenvectors associated with zero-eigenvalues, aside from rigid body modes. In order to properly identify the kernel of the stiffness matrix, the orthogonality conditions stated in Equation (A13) will be used together with some algebraic derivations. Since the discrete stiffness matrix is given by Equation ( 13), a zero-energy mode is a vector h g that satisfies

B(0 I )•h g = 0, I = 1,..., nint (B1)
It can be easily shown that a basis for the vector space of the discretized displacements is given by the following vectors:

⎛ s ⎞ e1 = ⎢ 0 ⎢ , e ⎛ 0 ⎞ = ⎢ s ⎢ , e ⎛ 0 ⎞ = ⎢ 0 ⎢ , e ⎛ x ⎞ = ⎢ 0 ⎢ , e ⎛ 0 ⎞ = ⎢ x ⎢ , e ⎛ 0 ⎞ = ⎢ 0 ⎢ ⎝ ⎠ 2 ⎝ ⎠ 3 ⎝ ⎠ 4 ⎝ ⎠ 5 ⎝ ⎠ 6 ⎝ ⎠ 0 ⎛ y ⎞ e7 = ⎢ 0 ⎢ , e 0 ⎛ 0 ⎞ = ⎢ y ⎢ , e s ⎛ 0 ⎞ = ⎢ 0 ⎢ , e 0 ⎛ z ⎞ = ⎢ 0 ⎢ , e 0 ⎛ 0 ⎞ = ⎢ z ⎢ , e x ⎛ 0 ⎞ = ⎢ 0 ⎢ ⎝ ⎠ 8 ⎝ ⎠ 9 ⎝ ⎠ 10 ⎝ ⎠ 11 ⎝ ⎠ 12 ⎝ ⎠ e13 = 0 ⎛ h 1 ⎞ 0 ⎛ 0 ⎞ y ⎛ 0 ⎞ 0 ⎛ h 2 ⎞ 0 z ⎛ 0 ⎞ ⎛ 0 ⎞ (B2) 15 ⎝ ⎠ 16 ⎝ ⎠ 17 ⎝ 2⎠ 18 ⎝ ⎠ e19 = ⎢ 0 ⎢ , e = ⎢ h ⎢ , e = ⎢ 0 ⎢ , e = ⎢ 0 ⎢ , e = ⎢ h ⎢ , e = ⎢ 0 ⎢ ⎝ ⎠ 20 0 ⎝ 3⎠ 0 21 ⎝ ⎠ h3 22 ⎝ ⎠ 0 23 ⎝ 4⎠ 0 24 ⎝ ⎠ h4
To demonstrate this, let us assume that there exists a combination of coefficients ai such that 24

}. ai ei = 0 (B3) i =1
Multiplying Equation (B3) by (b T 0 0), (0 b T 0), and (0 0 b T ), i = 1, 2, 3, successively, and making use of Equation (A13) yields

1 24 \ (b T 0 0)• }. ai ei i =1 1 24 \ = a3i +1 = 0 (0 b T 0)• }. ai ei i =1 1 24 \ = a3i +2 = 0, i = 1, 2, 3 (B4) (0 0 b T )• }. ai ei i =1 = a3i +3 = 0 a a a a a a ⎢ ⎢ ⎢
Repeating this operation with (h T 0 0), (0 h T 0), and (0 0 h T ), a = 1,..., 4, respectively:

(h T 0 0)• 1 24 \ }. ai ei i =1 1 24 \ = 8a3a+ 10 = 0 (0 h T 0)• }. ai ei i =1 1 24 \ = 8a3a+ 11 = 0, a = 1, . . . , 4 (B5) (0 0 h T )• }. ai ei i =1 = 8a3a+ 12 = 0
Combining Equations (B2)-(B5), it becomes obvious that

a1 = a2 = a3 = 0 (B6)
This shows that vectors (ei , i = 1,..., 24) are linearly independent, and hence they form a basis for the vector space of discretized displacements.

Let us now assume that vector h g belongs to the kernel of the stiffness matrix. This vector can be expanded in terms of the above base vectors as

24 h g = }. ci ei (B7) i =1
Using Equation (B1) together with Equation ( 14) and orthogonality condition (A13) yields

⎛ c 4 + h 1,x (Ç I )c 13 + h 2,x (Ç I )c 16 ⎢ c8 + h1,y (Ç I )c14 + h2,y (Ç I )c17 ⎢ ⎢ c 12 + h 1,z (Ç I )c 15 + h 2,z (Ç I )c 18 ⎢ ⎞ ⎢ ⎢ ⎢ ⎢ ⎢ = 0, I = 1,..., nint (B8) ⎢ c 5 ⎢ ⎢ +c7 + h 1,y (Ç I )c13 + h 1,x (Ç I )c14 + h 2,y (Ç I )c16 + h 2,x (Ç I )c 17 ⎢ ⎢ ⎢ ⎢ c9 +c11 + h1,z (Ç )c + h (Ç )c + h (Ç )c + h (Ç )c ⎢ I 14 ⎝ 1,y I 15 2,z I 17 2,y I 18 ⎠ c6 +c10 + h1,z (Ç I )c13 + h1,x (Ç I )c15 + h2,z (Ç I )c16 + h2,x (Ç I )c18
Evaluating the above equation at the nint different integration points implies that

c4 = c13 = c16 = 0 c8 = c14 = c17 = 0 c12 = c15 = c18 = 0 c5 +c7 = 0 c9 +c11 = 0 c6 +c10 = 0 (B9) ⎡ b T z ⎢ ⎢ ⎢ ⎢ ⎢ 8 ⎢ ⎢ and hence ⎛ s ⎞ ⎛ 0 ⎞ ⎛ 0 ⎞ ⎛ -y ⎞ ⎛ -z ⎞ ⎛ 0 ⎞ h g = c1 ⎢ 0 ⎢ +c ⎢ s ⎢ +c ⎢ 0 ⎢ +c ⎢ x ⎢ +c ⎢ 0 ⎢ +c ⎢ -z ⎢ ⎝ ⎠ 2 ⎝ ⎠ 0 0 3 ⎝ ⎠ s 5 ⎝ ⎠ 0 6 ⎝ ⎠ x 9 ⎝ ⎠ y ⎛ h 3 ⎞ +c19 ⎢ 0 ⎢ +c ⎛ 0 ⎞ ⎢ h ⎢ +c ⎛ 0 ⎞ ⎢ 0 ⎢ +c ⎛ h 4 ⎞ ⎢ 0 ⎢ +c ⎛ 0 ⎞ ⎢ h ⎢ +c ⎛ 0 ⎞ ⎢ 0 ⎢ (B10) ⎝ ⎠ 20 ⎝ 3⎠ 0 0 21 ⎝ ⎠ h3 22 ⎝ ⎠ 0 23 ⎝ 4⎠ 0 24 ⎝ ⎠ h4
This reveals that the kernel of the stiffness matrix for the SHB8PS element consists of the usual six rigid body modes (i.e. the first six vectors in Equation (B10)), and six hourglass modes as given by the last six vectors in Equation (B10) and illustrated in Figure 4. It should be noted that this formulation of the SHB8PS is valid for any set of nint integration points located along the same line Ç I = 1 I = 0, I = 1,..., nint, and comprising at least two integration points (nint�2). Note also that if a one-point quadrature rule had been used (i.e. nint = 1 and Ç I = 1 I = 0 I = 0), it would have resulted in 12 hourglass modes, as illustrated in Figure 3. Indeed, with a reduced integration scheme employing only a single integration point located at the center of the reference frame, the discrete gradient operator given by Equation ( 11) reduces to its constant part:

x 0 0 ⎤ ⎢ ⎢ ⎢ 0 b T 0 ⎢ ⎢ y ⎢ ⎢ ⎢ ⎢ 0 0 b T ⎢ Bc = ⎢ T T ⎢ (B11) ⎢b y b x 0 ⎢ ⎢ ⎢ 0 b T ⎢ b T ⎢ ⎢ z y ⎢ ⎣ b T T ⎦ z 0 b x
In this case, the same analysis as above would have resulted in

⎛ c 4 ⎢ c ⎢ ⎢ c 12 ⎢ ⎢ c +c ⎞ ⎢ ⎢ ⎢ ⎢ ⎢ = 0 (B12) ⎢ ⎢ 5 7 ⎢ ⎢ ⎢ ⎢ c9 +c11 ⎢ Hence ⎝ ⎠ c6 +c10 ⎛ s ⎞ h g = c1 ⎢ 0 ⎢ +c ⎛ 0 ⎞ ⎢ s ⎢ +c ⎛ 0 ⎞ ⎢ 0 ⎢ +c ⎛ -y ⎞ ⎢ x ⎢ +c ⎛ -z ⎞ ⎢ 0 ⎢ +c ⎛ 0 ⎞ ⎢ -z ⎢ ⎝ ⎠ 2 ⎝ ⎠ 0 0 3 ⎝ ⎠ s 5 ⎝ ⎠ 0 6 ⎝ ⎠ x 9 ⎝ ⎠ y ⎛ h 1 ⎞ +c13 ⎢ 0 ⎢ +c ⎛ 0 ⎞ ⎢ h ⎢ +c ⎛ 0 ⎞ ⎢ 0 ⎢ +c ⎛ h 2 ⎞ ⎢ 0 ⎢ +c ⎛ 0 ⎞ ⎢ h ⎢ +c ⎛ 0 ⎞ ⎢ 0 ⎢ ⎝ ⎠ 14 ⎝ 1⎠ 0 0 15 ⎝ ⎠ h1 16 ⎝ ⎠ 0 17 ⎝ 2⎠ 0 18 ⎝ ⎠ h2 ⎛ h 3 ⎞ +c19 ⎢ 0 ⎢ +c ⎛ 0 ⎞ ⎢ h ⎢ +c ⎛ 0 ⎞ ⎢ 0 ⎢ +c ⎛ h 4 ⎞ ⎢ 0 ⎢ +c ⎛ 0 ⎞ ⎢ h ⎢ +c ⎛ 0 ⎞ ⎢ 0 ⎢ (B13) ⎝ ⎠ 20 ⎝ 3⎠ 0 0 21 ⎝ ⎠ h3 22 ⎝ ⎠ 0 23 ⎝ 4⎠ 0 24 ⎝ ⎠ h4
As can be seen from Equation (B13) above, with a uniform reduced integration scheme, there are 12 hourglass modes (i.e. the last 12 vectors in Equation (B13)).

On the contrary, with a full integration scheme, functions ha,i (a = 1,..., 4; i = 1, 2, 3) would not have vanished at the integration points, and repeating the same analysis as before would give

c4 = c13 = c16 = c19 = c22 = 0 c8 = c14 = c17 = c20 = c23 = 0 c12 = c15 = c18 = c21 = c24 = 0 c5 +c7 = 0 c6 +c10 = 0 c9 +c11 = 0 (B14)
In this case, Equation (B14) clearly shows that only rigid body modes belong to the kernel of the stiffness matrix, and thus no rank deficiency is observed.

APPENDIX C: CO-ROTATIONAL COORDINATE FRAME

For the calculation of the stabilization terms (i.e. the stabilization stiffness matrix and stabilization internal load vector), several choices of co-rotational coordinate system are possible. The co-rotational frame proposed in Reference [START_REF] Belytschko | Assumed strain stabilization of the eight node hexahedral element[END_REF] was adopted in the current formulation, and some justification for the derivation of the rotation matrix associated with this local coordinate system is provided in this appendix. Note that this choice was motivated by several advantages, as stated in Section 2.5, and in particular, it allows simplified expressions for the above-mentioned stabilization stiffness matrix. The key property of this orthogonal co-rotational system, which is embedded in the element and rotates with it, is that it is chosen to be aligned with the reference coordinate system as illustrated in Figure 5. In the following, attention is focused on the extraction of the rotation matrix R that transforms a vector x in the global coordinate system (x , y, z) into a vector x ˜ in the co-rotational system (x ˜, y ˜, z ˜), using the vector transformation x ˜ = R• x.

The starting point of this derivation is the multiplicative decomposition of the deformation gradient (i.e. polar decomposition theorem), illustrated in Figure 5, which can be expressed as dx = F•dn = R T •S•dn (C1) This multiplicative decomposition of the deformation gradient ensures a local coordinate system that is aligned with the reference frame. Combining Equation (C1) with the relation dx ˜ = S•dn provides the desired transformation matrix:

dx = R T •dx ˜ ⇔ dx ˜ = R• dx (C3)
In order to extract this rotation matrix from the deformation gradient tensor, let us rewrite Equation ( 1) in an equivalent form: xi = 1 (s T +ÇK T +1K T +0K T + h1h T + h2h T + h3h T + h4h T )•xi (C4)

8 1 2 3 1 2 3 4
The components of the deformation gradient can then be expressed as axi

1 T T Fij = aÇ = 8 (K j + h a,Ç j h a )•xi (C5)
For the sake of simplicity, this deformation gradient is evaluated at the origin of the reference frame (Ç = 1 = 0 = 0):

⎡ K T T T ⎤ 1 •x1 K 2 •x1 K 3 •x1 1 T T T F(0) = 8 ⎢ K 1 •x2 K 2 •x2 K 3 •x2 ⎢ (C6) ⎢ ⎢ ⎣ ⎦ K T T T 1 •x3 K 2 •x3 K 3 •x3
This simplified form allows the calculation of the first two column vectors of the rotation matrix R T after adding a correction term to make them orthogonal (see Equations ( 39)-( 41)). The third column vector is simply obtained by the cross-product of the above normalized vectors (Equations ( 42)-( 43)).

APPENDIX D: LOCAL COORDINATE SYSTEM FOR ELASTICITY MATRIX

The local physical coordinate system, in which the material properties are specified, i.e. the modified elastic constitutive matrix given in Equation [START_REF] Reese | A new locking-free brick element technique for large deformation problems in elasticity[END_REF], is described here. For each integration point, such a local coordinate system is associated, in which the x -y plane corresponds to the element mid-plane defined by the 0-coordinate of the considered integration point. Several choices are possible; the adopted local coordinate system is depicted in Figure D1. The starting point is to define the element mid-plane corresponding to a given integration point K , 1�K �nint. This is determined by four points P K , 1�L�4, of coordinates xi ( P K ) such that In the above equation, K is defined as

K = 1 (1-0 K ) (D2)
where 0 K is the third coordinate of the current integration point, and xiI (i = 1,..., 3, I = 1,..., 8) are the nodal coordinates expressed in the global coordinate system. Then, the coordinates of the center O K of the local coordinate system are defined as In the quadrangle shown in Figure D1, M K , M K , M K , and M K are the barycentres (i.e. mid- points) of ( P K P K ), ( P K P K ), ( P K P K ), and ( P K P K ), respectively. This allows us to define the

xi (O K ) = 1 (xi ( P K )+ xi ( P K )+ xi ( P K )+ xi ( P K )), i = 1, 2, 3 (D3)
⎡ a xx ⎤ ⎡ e xx ⎤ ⎢ ⎢ ⎢ e Q ⎢ ⎢ayy⎢ ⎢ yy ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ a zz ⎢ ⎢ Q ⎢ ⎢ zz ⎢ r = ⎢ ⎢ , e Q = ⎢ ⎢ (E3) ⎢ ⎢ ⎢ Q Q ⎢ ⎢ a xy⎢ ⎢e xy +e yx ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ Q Q ⎢ ⎢ a yz ⎢ ⎢ ⎢ ⎣ ⎦ ⎢e yz +e zy ⎢ a xz ⎣ e Q Q ⎦
xz +e zx and the components of the quadratic part of the strain tensor are given by The components of the quadratic part of the strain tensor can then be discretized as (bd

T •Bi )(B T •�dk ) = bd T • B Q • �d where k j ij k=1 ⎡ Bi B T 0 0 ⎤ ⎡ bd 1 ⎤ ⎡ �d 1 ⎤ B Q ⎢ 0 Bi B T 0 ⎢ , bd = ⎢ bd2 ⎢ , �d = ⎢ �d2 ⎢ (E6) ij = ⎢ ⎣ j ⎢ ⎦ 0 0 Bi B T ⎣ ⎦ bd3 ⎣ ⎦ �d3
With these quadratic discrete gradient operators B Q , the contribution kr(0 I ) at integration point 0 I to the overall geometric stiffness matrix is given by kr(0 I ) = axx(0 I )B Q (0 )+ayy(0 )B Q (0 )+azz(0 )B Q (0 )+axy(0 )(B Q (0 )+B Q (0 )) The geometric stiffness matrix is finally obtained using the integration points as

n int Kr = }.
w(0 I ) J (0 I )kr(0 I ) (E7)

I =1

6 Figure 2 .

 62 Figure 2. Reference space (Ç, 1, 0) and physical space (x1, x2, x3) of the element.

Figure 3 .

 3 Figure 3. Hourglass modes in the x -direction for a one-point quadrature hexahedron.

Figure 4 .

 4 Figure 4. Hourglass modes in the x -direction for the SHB8PS element.

Figure 6 .

 6 Figure 6. Pinched hemispherical shell test: (a) geometric and material data as well as boundary conditions and loading specifications and (b) initial and deformed configurations.

12

 12

Figure 7 .

 7 Figure 7. The five-element mesh of the patch test; geometric and material data as well as the coordinates of the interior nodes 1, 2, 3, and 4 located on the top surface.

Figure 8 .

 8 Figure 8. Moment-curvature responses for an elastic-perfectly plastic beam bending problem: effect of through-thickness integration points and comparison with the analytical solution.

Figure 9 .

 9 Figure 9. Elastic-plastic bending of a cantilever beam: (a) geometric, material, loading, and boundary condition data and (b) load-displacement curves for different numbers of through-thickness integration points and comparison with reference results.

Figure 10 .

 10 Figure 10. Slit annular plate lifted by a line force P: (a) geometric, material, and loading data; (b) initial and deformed configurations; and (c) load-deflection curves for the slit annular plate: reference solution [90] and current SHB8PS results.

Figure 11 .

 11 Figure 11. Open-ended cylindrical shell subjected to radial pulling forces: (a) problem data; (b) deformed mesh; and (c) load-displacement curves of the open-ended cylinder subjected to pulling forces: reference data [90] and current SHB8PS results.

  point A in the three directions

Figure 12 .

 12 Figure 12. Twisted beam under out-of-plane loading: (a) problem data; (b) successive deformed configurations; and (c) load-displacement curves at the load point A of the twisted beam: reference solution [106] and current SHB8PS results.

Figure 13 .

 13 Figure 13. Lateral buckling of a cantilever beam: (a) problem data; (b) successive deformed configurations; and (c) load-displacement curves (fundamental and bifurcated paths): reference solution [107] and current SHB8PS results.

Figure 14 (

 14 Figure14(a) depicts an elastic deep circular arch subjected to a point load at the apex of the arch. The geometric parameters, material properties, boundary conditions, and loading are all described in this figure. Owing to the snap-through behavior exhibited by this test problem, it has been frequently used in the literature with different types of boundary conditions. The clamped-hinged version studied here is deemed more interesting and challenging, since its asymmetric boundary conditions allow non-symmetric buckling to occur. In order to follow the load-deflection curves beyond the limit points, the path-following Riks method based on an arc-length control parameter is adopted. The entire post-buckling response was provided in Reference[START_REF] Klinkel | A geometrical non-linear brick element based on the EAS-method[END_REF] using a robust strategy. Because of the lack of symmetry, the arch is modeled here entirely by means of a 40×1 mesh with only a single element along the thickness. In the literature, similar mesh densities have

Figure 14 .

 14 Figure 14. Clamped-hinged deep circular arch under central point load: (a) problem data; (b) successive deformed configurations; and (c) load-displacement curves of the point under vertical load: reference solution [25] and current SHB8PS results.

Figure 15 .

 15 Figure 15. Elastic-plastic buckling of a stiffened cylinder: (a) problem data specification; (b) first Euler's buckling mode; and (c) pressure to Euler's critical load ratio versus radial displacement at point A: reference solution [109] and current SHB8PS results.

  Hallquist's form, reads b T T ( a N1 a N2 a N3 a N4 a N5 a N6 a N7 a N8 \ i = N ,i |Ç=1=0=0 = axi axi axi axi axi axi axi axi |Ç=1=0=0 (A1)

  , a, fJ = 1,...,[START_REF] Liu | Finite element stabilization matrices-a unification approach[END_REF] The above third consistency condition follows simply from b T

j

  in which S is a symmetric matrix and R is an orthogonal matrix representing a pure rotation. The latter does not contribute to the deformation since e= 1 (F T •F-I) = 1 (S T •R•R T •S-I) = 1 (S 2 -I)

L 8 L, 3 Figure D1 .

 83D1 Figure D1. Schematic representation of the local material coordinate system associated with the K th integration point of the element.

  , e1, of the local coordinate system, parallel to (O K M K ), while e2 is defined parallel to (O K M K ), see FigureD1. A correction term ec is then added to vector e2, so that e1 and (e2 + ec) are orthogonal, which gives The third base vector is simply obtained by the cross-producte3 = e1 ×(e2 +ec) (D5)Finally, the rotation matrix R mat that maps a vector in the global coordinate system to the local coordinate system is derived, after normalization, in terms of the components of the base vectors byR mat e 1i mat e 2i + e ci e 3i 1i = ±e1± , R 2i = ±e2 +ec± , R mat = ±e3± , i = 1, 2, 3(D6)APPENDIX E: GEOMETRIC STIFFNESS MATRIX FOR BUCKLING ANALYSISIn this appendix, the geometric stiffness matrix is derived for the SHB8PS element. This geometric stiffness matrix Kr is to be added to the regular tangent stiffness matrix K in a usual structural see Equation (60), for instance. Recall that the geometric stiffness matrix originates from the linearization of the virtual work principle and is due to the non-linear (quadratic) part of the strain tensor. In its continuum form, it reads: r Kr(bu, �u) = Oe r r : ∇bu T • ∇�u dO = Oe r : e Q (bu, �u) dO (E1) Using the vector form of the stress tensor and the quadratic part of the strain tensor, respectively, Equation (E1) can be rewritten as where r Kr(bu, �u) = Oe r T • e Q (bu, �u) dO (E2)

3 ij

 3 (bu, �u) = }. k=1 buk,i �uk, j = buk,i �uk, j (E4)Using the discrete form of the displacement gradient, as given in Equation (A18), we obtain buk,i= (b T + ha,i c T )•bdk = B T •bdk i a i (E5) �uk, j = (b T + ha, j c T )• �dk = B T • �dk j a j

3 ij

 3 (bu, �u) = }.

  (0 I )(B Q (0 )+B Q (0 ))+axz(0 )(B Q (0 )+B Q (0 ))

Table II .

 II Normalized displacement at point A of the pinched hemispherical shell.

		SHB8PS				SHB8PS
	Number of	previous formulation	HEX8	HEXDS	H8-ct-cp	current formulation
	elements	ux /u ref	ux /u ref	ux /u ref	ux /u ref	ux /u ref
	12	0.0629	0.0005		0.05	0.8645
	27	0.0474	0.0011			1.0155
	48	0.1660	0.0023	0.408	0.35	1.0098
	75	0.2252	0.0030	0.512	0.58	1.0096
	192	0.6332	0.0076	0.701	0.95	1.0008
	363	0.8592	0.0140	0.800		1.0006
	768	0.9651	0.0287			1.0006
	1462	0.9910	0.0520			1.0009

Table IV .

 IV u ref Out-of-plane bending patch test: normalized displacements of interior nodes.

				1		u 2	ref /u 2	u 3	ref /u 3
	Node Mesh layout	HEX8 SHB8PS Q1SPs HEX8 SHB8PS HEX8 SHB8PS
	1	1×1	256.25	5.90	4.77	10.17	3.49	3.50	4.67
		2×2	2.24	1.47	1.08	1.45	0.98	2.08	1.31
		4×4	1.27	1.08	1.05	1.24	0.91	1.47	1.17
		8×8	1.23	1.01	1.01	1.23	0.98	1.36	1.02
	2	1×1	23.90	0.87	0.53	54.59	0.64	1.26	1.82
		2×2	0.68	0.87	0.99	0.84	0.92	1.08	1.03
		4×4	0.91	0.99	1.00	0.81	1.00	1.03	1.01
		8×8	0.92	1.00	1.00	0.84	1.00	1.02	1.00
	3	1×1 2×2	22.21 0.85	0.18 0.74	-0.06 0.95	103.29 0.82	0.15 0.79	1.22 1.07	0.81 1.05
		4×4	1.00	0.98	1.00	1.09	0.99	1.03	1.02
		8×8	1.01	1.00	1.01	1.07	1.00	1.02	1.00
	4	1×1	347.48	1.26	2.23	371.28	2.66	1.28	2.07
		2×2	1.10	1.24	1.09	0.44	0.70	1.22	1.12
		4×4	0.97	1.03	1.01	0.87	1.02	1.10	1.04
		8× 8	0.95	1.01	1.01	0.87	1.00	1.07	1.01

Table V .

 V Out-of-plane bending patch test: normalized stress values on the top surface.

		a 11 /a ref 11	ref a 22 /a 22	ref a 12 /a 12
	Mesh layout	HEX8	SHB8PS	HEX8	SHB8PS	HEX8	SHB8PS
	1×1	133.04	0.63	263.39	2.53	435.35	1.22
	2×2	0.88	1.20	1.28	1.28	0.42	1.09
	4×4	0.92	0.90	0.94	0.99	1.21	1.00
	8× 8	0.78	0.99	0.69	0.99	1.18	1.01
	contributions (see References [7, 37-39, 41, 42, 79-86]

Table VI .

 VI Stress-strain curve for the elastic-plastic behavior with isotropic hardening.

	ε (%)	a (MPa)	ε(%)	a (MPa)	ε (%)	a (MPa)
	0.1300	260	0.2560	360	2.7043	480
	0.1470	285	0.3362	380	4.2872	500
	0.1506	290	0.3951	390	10.5084	540
	0.1783	320	0.7003	420	16.1733	560
	0.2081	340	1.0805	440