Robin Bouclier
email: robin.bouclier@insa-lyon.fr

| Thomas Elguedj

Alain Combescure

An isogeometric locking-free NURBS-based solid-shell element for geometrically nonlinear analysis

In this work, we develop an isogeometric NURBS-based solid-shell element for the geometrically nonlinear static analysis of elastic shell structures. A single layer of continuous 3D elements through the thickness of the shell is considered and the order of approximation in that direction is chosen to be equal to two. A complete 3D constitutive relation is assumed. The objective is to develop a highly accurate low-order element for coarse meshes. We propose an extension of the mixed method of Bouclier et al. [START_REF] Bouclier | Efficient isogeometric NURBS-based solid-shell elements: Mixed formulation and B-method[END_REF] to deal with locking in the context of large rotations and large displacements. The main idea is to modify the interpolation of the average through the thickness of the stress components. It is also necessary to stabilize the element in order to avoid the occurrence of spurious zero-energy modes. This was achieved, for the quadratic version, through the adjunction of artificial elementary stabilization stiffnesses.

The result is an element of order 2 which is at least as accurate as standard NURBS shell elements of order 4. Linear and nonlinear test calculations have been carried out along with comparisons with other published NURBS and classical techniques in order to assess the performance of the element.

K E Y W O R D S

Isogeometric analysis, Solid-shell element, Locking, Geometrically nonlinear analysis, Mixed methods, Hourglass control

| INTRODUCTION

Today, due to both the emergence of numerical simulation for analysis and the increasing use of Computer-Aided Design (CAD) software in design, it seems judicious to address the calculation of shell structures in the framework of IsoGeometric Analysis (IGA). Indeed, in this framework which was introduced in Hughes et al. [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry, and mesh refinement[END_REF], the calculations are carried out using NURBS functions, which are the most common CAD technology, instead of Lagrange polynomials.

More precisely, since CAD uses surface representations and boundaries, this approach is especially suited for surfacebased geometries, which is the case of shells. Thus, the IGA approach enables one to address both the design and the analysis of shells using exactly the same geometric models. The meshing stage, which was sometimes delicate and expensive, becomes immediate, and the geometry of the analysis model can be represented exactly regardless of the refinement of the mesh. In addition to that geometric aspect, NURBS functions appear to have higher-order continuity at knots. The use of the well-known k -refinement (smooth order elevation, see Cottrell et al. [START_REF] Cottrell | Studies of refinement and continuity in isogeometric structural analysis[END_REF]), in which C p-1 continuity is achieved using discretizations of order p, leads to improved accuracy and robustness compared to traditional finite elements. Based on this principle, numerous works have been initiated to assess the effectiveness of the method in the case of shells. First, due to the ease with which C 1 continuity can be implemented with these functions, Kiendl et al. [START_REF] Kiendl | Isogeometric shell analysis with Kirchhoff-Love elements[END_REF][START_REF] Kiendl | The bending strip method for isogeometric analysis of Kirchhoff-Love shell structures comprised of multiple patches[END_REF] developed NURBS elements of the Kirchhoff-Love type. At the same time, Benson et al. [START_REF] Benson | Isogeometric shell analysis: the Reissner-Mindlin Shell[END_REF][START_REF] Benson | A large deformation, rotation-free, isogeometric shell[END_REF] focused more on elements of the Reissner-Mindlin type, while Hughes et al. [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry, and mesh refinement[END_REF] and Echter and Bischoff [START_REF] Echter | Numerical efficiency, locking and unlocking of NURBS finite elements[END_REF] were able to demonstrate the effectiveness of the approach for structures calculated using 3D NURBS solid elements. Subsequently, Dornisch et al. [START_REF] Dornisch | Isogeometric Reissner-Mindlin shell analysis with exactly calculated director vectors[END_REF] proposed a means of describing the normal to the shell exactly in the Reissner-Mindlin model. Finally, Echter et al. [START_REF] Echter | A hierarchic family of isogeometric shell finite elements[END_REF], Hosseini et al. [START_REF] Hosseini | An isogeometric solid-like shell element for nonlinear analysis[END_REF], Bouclier et al. [START_REF] Bouclier | Efficient isogeometric NURBS-based solid-shell elements: Mixed formulation and B-method[END_REF] and Caseiro et al. [START_REF] Caseiro | On the Assumed Natural Strain Method to alleviate locking in solid-shell NURBS-based finite elements[END_REF] developed finite elements which are capable of taking into account the behavior through the thickness of the shell.

In practice, in spite of the advantage of a reduction in the number of degrees of freedom, the implementation of a shell structural model based on a discretization of the mid-surface alone is difficult. For example, coupling these shell elements with 3D solid elements, introducing general 3D constitutive relations, extracting information about the behavior through the thickness and describing contact conditions are delicate endeavors. In classical finite element analysis, elements known as solid-shell elements have been widely considered (see, for example, [START_REF] Legay | Elastoplastic stability of shells using the physically stabilized finite element SHB8PS[END_REF][START_REF] Abed-Meraim | An improved assumed strain solid-shell element formulation with physical stabilization for geometric nonlinear applications and elastic-plastic stability 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 II -Non linear applications[END_REF][START_REF] Klinkel | A geometrical non-linear brick element based on the EAS-method[END_REF][START_REF] Wriggers | A note on enhanced strain methods for large deformations[END_REF][START_REF] Reese | A large deformation solid-shell concept based on reduced integration with hourglass stabilization[END_REF][START_REF] Bassa | A new nine-node solid-shell finite element using complete 3D constitutive laws[END_REF]).

Here, the idea is to take a continuous 3D solid element and, in order to minimize the computation cost, consider a single element layer through the thickness. The advantage of these elements is that they have only nodal displacements as degrees of freedom, which makes their implementation easier for complex simulations (coupling with solid elements, varying shell thickness, geometric, material and contact nonlinearities). However, there is a difficulty in dealing with locking, which these elements are very sensitive to. Since the element relies on a full 3D model, all the locking mechanisms are possible: transverse and in-plane shear locking, membrane locking and thickness locking. Resorting to the NURBS context does not eliminate this problem. Indeed, it has been shown that NURBS elements are sensitive to the same types of locking as classical finite elements based on Lagrange-type polynomials (see, in particular, Echter and Bischoff [START_REF] Echter | Numerical efficiency, locking and unlocking of NURBS finite elements[END_REF] for shells and Elguedj et al. [START_REF] Elguedj | B-bar an F-bar projection methods for nearly incompressible linear and non linear elasticity and plasticity using higher order NURBS element[END_REF] for nearly incompressible materials). The initial strategy in early works on NURBS shells consisted in increasing the order of interpolation (see, again, [START_REF] Kiendl | Isogeometric shell analysis with Kirchhoff-Love elements[END_REF][START_REF] Benson | Isogeometric shell analysis: the Reissner-Mindlin Shell[END_REF][START_REF] Dornisch | Isogeometric Reissner-Mindlin shell analysis with exactly calculated director vectors[END_REF]). However, doing so merely reduces locking without eliminating the issue completely: even high-order elements are not completely locking-free and the accuracy depends on such critical parameters as the slenderness or the curvature of the structure. In addition, increasing the order of approximation leads to an increase in the number of Gauss points for calculating the integrals, thus making these elements costly, especially since the objective is to perform nonlinear calculations.

Going back to solid-shell elements in the context of the traditional finite element method, many authors have focused on the development of locking-free 8-node elements (i.e. elements with linear interpolation functions). In most of these cases, the locking issue could be resolved by using assumed strain methods [START_REF] Bathe | A four node plate bending element based on Mindlin/Reissner theory and a mixed interpolation[END_REF][START_REF] Bucalem | Higher-order MITC general shell elements[END_REF][START_REF] Belytschko | Assumed strain stabilization of the eight node heaxahedral element[END_REF] coupled with reduced-integration schemes [START_REF] Zienkiewicz | Reduced integration technique in general analysis of plates and shells[END_REF][START_REF] Hughes | A simple and efficient finite element for plate bending[END_REF] with hourglass control [START_REF] Flanangan | A Uniform strain hexahedron and quadrilateral with orthogonal hourglass control[END_REF][START_REF] Belytschko | Hourglass control in linear and nonlinear problems[END_REF]. One can mention the solid-shell elements of Combescure et al. [START_REF] Legay | Elastoplastic stability of shells using the physically stabilized finite element SHB8PS[END_REF][START_REF] Abed-Meraim | An improved assumed strain solid-shell element formulation with physical stabilization for geometric nonlinear applications and elastic-plastic stability analysis[END_REF], Alves de Sousa et al. [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 -Non linear applications[END_REF], Wriggers and Reese [START_REF] Wriggers | A note on enhanced strain methods for large deformations[END_REF], Reese [START_REF] Reese | A large deformation solid-shell concept based on reduced integration with hourglass stabilization[END_REF] and Bassa et al. [START_REF] Bassa | A new nine-node solid-shell finite element using complete 3D constitutive laws[END_REF]. In the NURBS context, however, far fewer studies have been undertaken in order to eliminate locking in shells, and even fewer using the solid element approach. Some works were carried out for the simplified case of beams (see, for instance, Echter and Bischoff [START_REF] Echter | Numerical efficiency, locking and unlocking of NURBS finite elements[END_REF] using a DSG method, Beirão da Veiga et al. [START_REF] Beirão Da Veiga | Avoiding shear locking for the Timoshenko beam problem via isogeometric collocation methods[END_REF][START_REF] Auricchio | Locking-free isogeometric collocation methods for spatial Timoshenko rods[END_REF] using a collocation method and Bouclier et al. [START_REF] Bouclier | Locking free isogeometric formulations of curved thick beams[END_REF] using reduced-integration techniques and B-projection strategies). In the context of shell structural models, one of the rare contributions was that of Echter et al. [START_REF] Echter | A hierarchic family of isogeometric shell finite elements[END_REF], who developed a hierarchical family of NURBS shell elements. Closer to the solid-shell context, one can mention the work of Hosseini et al. [START_REF] Hosseini | An isogeometric solid-like shell element for nonlinear analysis[END_REF], who developed a solidlike Bézier shell element. However, only Bouclier et al. [START_REF] Bouclier | On the development of NURBS-based isogeometric solid shell elements: 2D problems and preliminary extension to 3D[END_REF][START_REF] Bouclier | Efficient isogeometric NURBS-based solid-shell elements: Mixed formulation and B-method[END_REF], using a mixed approach, and more recently Caseiro et al. [START_REF] Caseiro | On the Assumed Natural Strain Method to alleviate locking in solid-shell NURBS-based finite elements[END_REF], using an Assumed Natural Strain method, seem to have applied the pure solid-shell approach to the NURBS context. Up until now, these works were limited to linear elastic structural analysis.

In this context, our aim is to develop a solid-shell NURBS element for the analysis of geometrically nonlinear problems. The shell is assumed to be elastic, linear and isotropic, but is subjected to large rotations and large displacements. A full 3D constitutive relation is considered, and functions of degree 2 are assumed through the thickness.

The objective is to build a highly accurate low-order element for coarse meshes. In order to do that, we propose an extension of the mixed method presented in [START_REF] Bouclier | Efficient isogeometric NURBS-based solid-shell elements: Mixed formulation and B-method[END_REF] to handle the locking issue in the context of large rotations and large displacements. We also stabilize the element by adding artificial stiffnesses to prevent the occurrence of any spurious zero-energy modes.

The paper is organized as follows: after this introduction, Section 2 reviews the construction of the element of [START_REF] Bouclier | Efficient isogeometric NURBS-based solid-shell elements: Mixed formulation and B-method[END_REF] in the context of small deformations. Then its stabilization in the case of small deformations is presented in Section 3. Section 4 extends these considerations to the geometrically nonlinear context. Linear and nonlinear numerical tests are presented in Section 5 in order to assess the accuracy of the element. Finally, Section 6 makes concluding remarks on this work.

| THE MIXED SOLID-SHELL NURBS ELEMENT FOR ELASTICITY UNDER SMALL DEFORMATIONS

In this section, we briefly summarize the construction of the mixed solid-shell NURBS element developed in Bouclier et al. [START_REF] Bouclier | Efficient isogeometric NURBS-based solid-shell elements: Mixed formulation and B-method[END_REF]. The interested reader is referred to that reference for further details. This element is suitable for the linear elastic calculation of shell structures under small deformations. Prior to focusing on the mixed element itself, we begin by introducing the shell problem and its discretization using the IGA concept based on NURBS functions.

| The reference problem

We undertake to study the elastic behavior of a thick or thin shell of thickness h described in coordinate system (x, y, z)

(see Figure 1). The domain occupied by the shell is denoted Ω. At any point M , we consider the local coordinate system (e 1 , e 2 , e 3 ), where e 1 and e 2 are two vectors which are tangent to the mid-surface and e 3 is the normal vector. We will use the terms "thickness" in reference to direction e 3 and "lengths" for directions e 1 and e 2 . The local coordinate through the thickness is denoted z ( ∈ -h 2 , h 2 ). We will also use ξ, η and χ, the normalized coordinates (∈ [0, 1]) along the lengths and the thickness respectively. The associated domain (ξ, η, χ) is called the isoparametric domain.

The shell is assumed to be subjected to a body force f in Ω and we prescribe a pressure F over boundary Γ F and a displacement U d over boundary Γ u . These boundaries are such that Γ u ∩ Γ F = ∅ and Γ u ∪ Γ F = ∂ Ω. Finally, we consider a homogeneous, isotropic material with Young's modulus E and Poisson's coefficient ν. We want to solve this problem using the recent concept of IGA along with the solid-shell element approach [START_REF] Bouclier | On the development of NURBS-based isogeometric solid shell elements: 2D problems and preliminary extension to 3D[END_REF][START_REF] Bouclier | Efficient isogeometric NURBS-based solid-shell elements: Mixed formulation and B-method[END_REF][START_REF] Caseiro | On the Assumed Natural Strain Method to alleviate locking in solid-shell NURBS-based finite elements[END_REF]. In order to do that, we will discretize the shell using a single layer of solid 3D NURBS elements through the thickness.

| Description of the geometry using NURBS

As is usually the case in the IGA approach, the geometry of the shell is described using the NURBS technology. In this section, we present a brief review of the bases of that technique. For further details, the interested reader should consult the references listed below. The concept was first introduced in Hughes et al. [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry, and mesh refinement[END_REF] and, more recently, was formalized in the book by Cottrell et al. [START_REF] Cottrell | Isogeometric analysis: Toward Integration of CAD and FEA[END_REF], where a detailed description can be found. In recent years, this has led to a series of many publications (see, for instance, [START_REF] Cottrell | Studies of refinement and continuity in isogeometric structural analysis[END_REF][START_REF] Kiendl | Isogeometric shell analysis with Kirchhoff-Love elements[END_REF][START_REF] Kiendl | The bending strip method for isogeometric analysis of Kirchhoff-Love shell structures comprised of multiple patches[END_REF][START_REF] Benson | Isogeometric shell analysis: the Reissner-Mindlin Shell[END_REF][START_REF] Echter | Numerical efficiency, locking and unlocking of NURBS finite elements[END_REF][START_REF] Benson | A large deformation, rotation-free, isogeometric shell[END_REF][START_REF] Dornisch | Isogeometric Reissner-Mindlin shell analysis with exactly calculated director vectors[END_REF][START_REF] Echter | A hierarchic family of isogeometric shell finite elements[END_REF][START_REF] Hosseini | An isogeometric solid-like shell element for nonlinear analysis[END_REF][START_REF] Bouclier | Efficient isogeometric NURBS-based solid-shell elements: Mixed formulation and B-method[END_REF][START_REF] Caseiro | On the Assumed Natural Strain Method to alleviate locking in solid-shell NURBS-based finite elements[END_REF][START_REF] Elguedj | B-bar an F-bar projection methods for nearly incompressible linear and non linear elasticity and plasticity using higher order NURBS element[END_REF][START_REF] Beirão Da Veiga | Avoiding shear locking for the Timoshenko beam problem via isogeometric collocation methods[END_REF][START_REF] Auricchio | Locking-free isogeometric collocation methods for spatial Timoshenko rods[END_REF][START_REF] Bouclier | Locking free isogeometric formulations of curved thick beams[END_REF][START_REF] Bouclier | On the development of NURBS-based isogeometric solid shell elements: 2D problems and preliminary extension to 3D[END_REF][START_REF] Hughes | Efficient quadrature for NURBS-based isogeometric analysis[END_REF][START_REF] Auricchio | A Simple algorithm for obtaining nearly optimal quadrature rules for NURBS-based isogeometric analysis[END_REF]).

NURBS functions, which are a generalized version of B-spline functions, have become a standard for geometric modeling in CAD and computer graphics (see, for instance, Cohen et al. [START_REF] Cohen | Discrete B-spline and subdivision techniques in computer aided geometric design and computer graphics[END_REF], Piegl and Tiller [START_REF] Piegl | The NURBS Book (Monographs in Visual Communication)[END_REF] and Farin [START_REF] Farin | Curves and Surfaces for CAGD, A Practical Guide[END_REF]). These functions provide exact representations of many shapes used in engineering, such as conical sections. They can be viewed as rational projections of high-order B-splines and, thus, possess many of the properties associated with Bsplines, including the most interesting one which is their high degree of continuity. Further details can be found in Cohen et al. [START_REF] Cohen | Discrete B-spline and subdivision techniques in computer aided geometric design and computer graphics[END_REF], Rogers [START_REF] Rogers | An introduction to NURBS With Historical Perspective[END_REF] and Piegl and Tiller [START_REF] Piegl | The NURBS Book (Monographs in Visual Communication)[END_REF], as well as in Cottrell et al. [START_REF] Cottrell | Isogeometric analysis: Toward Integration of CAD and FEA[END_REF].

If N A , A ∈ {1, 2, .., n } denotes the n 3D NURBS functions, ω A , A ∈ {1, 2, .., n } the associated weights and P A , A ∈ {1, 2, .., n } the associated control points with coordinates x A in coordinate system (x, y, z), the geometry of the shell is described through the position vector M, which is defined as:

M = n A=1 N A x A , (1)
where the NURBS functions N A , A ∈ {1, 2, .., n } are obtained from the B-spline functions N A through:

N A = N A w A n A=1 N A w A . (2)
Then in order to define the 3D B-spline functions N A at control point P A , all one needs to do is use the tensor product, which consists simply in multiplying the 1D B-spline functions associated with that point in the three spatial directions. k , k ∈ 1, 2, .., n e 3 denote the n e 1 , n e 2 and n e 3 1D B-spline functions associated with directions e 1 , e 2 and e 3 respectively, at control point P A which corresponds to the i th , j th and k th control points in directions e 1 , e 2 and e 3 respectively, one has:

N A = M e 1 i M e 2 j M e 3 k . (3)
The 1D B-spline functions are built using a knot vector. Each knot vector associated with a direction is defined in the isoparametric domain (ξ, η, χ). For example, for the first direction, one uses knot vector Ξ = ξ 1 , ξ 2 , .., ξ ne 1 +p+1 , where ξ l ∈ is the l th knot, l is the knot index (l = 1, 2, .., n e 1 + p + 1) and p is the polynomial degree of the functions

M e 1 i , i ∈ 1, 2, .., n e 1 .
There can be more than one knot in the same location of the parametric space. If m is the multiplicity of a given knot, the functions are of continuity class C p-m at that location. If the knots are regularly spaced, the knot vector is said to be uniform. A knot vector whose first and last knots are of multiplicity p + 1 is said to be open. In that case, the basis is interpolating at the end knots of the interval, which simplifies the application of the boundary conditions. For the sake of simplicity, our work considers only knot vectors which are both uniform and open. The 1D B-spline basis functions for a given order p are defined recursively from the knot vector using the Cox-de Boor recursion formula (see, for example, Cohen et al. [START_REF] Cohen | Discrete B-spline and subdivision techniques in computer aided geometric design and computer graphics[END_REF]). One starts from functions which are constant by parts (p = 0):

M e 1 i ,0 (ξ) = 1 if ξ i ≤ ξ < ξ i +1 , 0 otherwise (4)
Then, for p = 1, 2, 3, .., one builds:

M e 1 i ,p (ξ) = ξ -ξ i ξ i +p -ξ i M e 1 i ,p-1 (ξ) + ξ i +p+1 -ξ ξ i +p+1 -ξ i +1 M e 1 i +1,p-1 . (5)
Contrary to classical finite elements whose shape functions are usually chosen to be linear, here, in order to benefit from the superior approximation properties of the NURBS functions, we choose them to be of a polynomial degree at least equal to two in the three spatial directions. Regarding continuity, the mesh refinement is of the k -refinement type, which means that elements are added while maintaining the highest order of continuity of the NURBS functions, i.e. C p-1 at the knots. The positions and associated weights of the control points can be set out to give an exact representation of conical shell geometries (cylinder, cone, sphere...). Subsequently, these geometries are preserved during mesh refinement. For a good overview of mesh generation and refinement, see Cottrell et al. [START_REF] Cottrell | Studies of refinement and continuity in isogeometric structural analysis[END_REF]. One should note that because of the structured nature of the NURBS functions local refinement cannot be applied directly. Today, it seems that a remedy for this can be to use T-splines (see, for example, Scott et al. [START_REF] Scott | Local refinement of analysis-suitable T-splines[END_REF]) or hierarchical B-splines (see, for example, Vuong et al. [START_REF] Vuong | A hierarchical approach to adaptive local refinement in isogeometric analysis[END_REF]). These techniques will not be considered in this work. Finally, one can note that these CAD functions are incapable of representing helixes or helicoidal functions exactly. Geometry specialists are working on this problem (see, for example, Letcher and Shook [START_REF] Letcher | NURBS considered harmful for gridding (Alternative Offered[END_REF]).

| The mixed finite element NURBS formulation

Now that we have defined the shell problem and its discretized geometry using NURBS functions, let us proceed to its resolution. In order to carry out this step effectively, i.e. while keeping the adverse locking phenomenon at a minimum, the solid-shell element of [START_REF] Bouclier | Efficient isogeometric NURBS-based solid-shell elements: Mixed formulation and B-method[END_REF] is based on a mixed method.

| Preliminaries and notations

The displacement solution U of the problem at point M is expressed in coordinate system (x, y, z) as:

U = {u v w } T (x,y,z) . ( 6
)
This displacement is sought in functional space U, which, along with functional space V of the test displacement fields, is defined as follows:

U = V ∈ H 1 (Ω) 3 , V | Γu = U d ; V = V ∈ H 1 (Ω) 3 , V | Γu = 0 . (7)
Then, one can write l (U), which corresponds to the work of the external loads, as:

l (U) = ∫ Ω U T fdΩ + ∫ Γ F U T FdΓ F . (8)
In continuous variational formulations, strains and stresses are expressed using a tensor-type notation. Throughout the paper, this will be the preferred notation in the sections labeled "Continuous version". Thus, in (x, y, z), we use the linearized strain written as ε x and the Cauchy stress written as σ x . For the implementation part, we prefer to use a notation of the Voigt type, in which the same quantities are written in vector form:

ε x = ε x x ε y y ε z z √ 2ε x y √ 2ε y z √ 2ε x z T ; σ x = σ x x σ y y σ z z √ 2σ x y √ 2σ y z √ 2σ x z T . (9)
This, rather than the tensor notation, will be used in the sections labeled "Discrete version" in the rest of the paper.

The strain is obtained classically from the displacement as follows:

ε x x = u, x ; ε y y = v , y ; ε z z = w , z ; √ 2ε x y = 1 √ 2 (u, y +v , x ) ; √ 2ε y z = 1 √ 2 (v , z +w , y ) ; √ 2ε x z = 1 √ 2 (u, z +w , x ) ; (10)
Let us mention again that a constitutive relation of the elastic, linear, homogeneous and isotropic type will be assumed throughout this work. In order to do that, a complete 3D behavior law is used. The relation on the tensor level is written as: σ x = C ε x , where C is Hooke's operator. The relation is obtained, using Voigt's notations [START_REF] Echter | A hierarchic family of isogeometric shell finite elements[END_REF], through the introduction of operator D such that σ x = D ε x .

In order to make the shell's strains and stresses meaningful (i.e. allocate the membrane, shear and thickness contributions properly) and, thus, describe the shell in the same manner as with a structural model, it is preferable to work in the local coordinate system (e 1 , e 2 , e 3 ) (see [START_REF] Bouclier | Efficient isogeometric NURBS-based solid-shell elements: Mixed formulation and B-method[END_REF] for details regarding the construction of this local system). In that system, the strain and stress tensors are denoted ε t and σ t and the associated vectors are:

ε t = ε 11 ε 22 ε 33 √ 2ε 12 √ 2ε 23 √ 2ε 13 T ; σ t = σ 11 σ 22 σ 33 √ 2σ 12 √ 2σ 23 √ 2σ 13 T . (11)
These are obtained through a change of basis starting from global quantities [START_REF] Echter | A hierarchic family of isogeometric shell finite elements[END_REF]. This transformation, denoted R, is defined on the tensor level and is such that:

ε t = P T ε x P = R
The same relation can also be expressed using Voigt's notations, which boils down to calculating operator R as a function of the components P i j of P so that ε t = R ε x and σ t = R σ x .

| The continuous version

In mixed methods, the stress field and the usual displacement field are sought simultaneously. This is interesting when it comes to relieving locking because one is given the possibility of seeking the stress field in a space which is different from that of the strain field (which derives from the displacement field). Consequently, the displacement-strain relation [START_REF] Hosseini | An isogeometric solid-like shell element for nonlinear analysis[END_REF] is enforced in a weak sense, contrary to basic displacement-type formulations in which it is enforced in a strong sense. This is how, in Herrmann et al. [START_REF] Herrmann | Elasticity equations for incompressible and nearly incompressible materials by a variational theorem[END_REF], the mixed formulation for nearly incompressibility problems was introduced with the pressure as an additional unknown to be sought in a space of a lesser degree than the displacements, leading to elements which are free from volumetric locking.

In our case of a model intended for the mechanics of continuous media, the basic mixed problem consists in finding both U ∈ U and σ t ∈ S, where S is the stress space defined such that:

S = σ t ∈ L 2 (Ω) 6 , (13)
which minimizes the following mixed functional over V ⊗ S:

Π mi x ed U, σ t = ∫ Ω σ t : ε t dΩ - 1 2 ∫ Ω σ t : D -1 σ t dΩ -l (U) . (14)
Here and from now on, we use the notation σ : ε to designate the twice-contracted product of stresses and strains.

The minimization of ( 14) leads to the expression of the problem in the weak form: find U ∈ U and σ t ∈ S such that:

∫ Ω σ t * : ε t dΩ + ∫ Ω ε t * : σ t dΩ - ∫ Ω σ t * : C -1 σ t dΩ = l (U * ), U * ∈ V and σ t * ∈ S. (15)
In order to extract the contributions which induce locking in the solid-shell element, we propose to use the average stresses through the shell's thickness. Indeed, one would expect the locking mechanisms involved to be (i) membrane locking, (ii) shear locking and (iii) thickness locking. Globally, the terms which are responsible for these mechanisms seem to be the averages through the thickness of (i) σ 11 and σ 22 , (ii) σ 12 , σ 23 and σ 13 , and (iii) σ 33 . This idea of making use of an averaging operator to extract the locking contributions was initially given in [START_REF] Bouclier | On the development of NURBS-based isogeometric solid shell elements: 2D problems and preliminary extension to 3D[END_REF]. Thus, one writes:

σ t = σ t mid ⊥ + σ t mid , (16)
where σ t mid is the average stress through the thickness and σ t mid ⊥ is the complement to obtain the true stress. In order to treat locking, one modifies the interpolation of σ t mid . Thus, that contribution, hereafter denoted simply σ, becomes the only unknown stress part of the mixed formulation. The other contribution σ t mid ⊥ comes from the displacement and is equal to:

σ t mid ⊥ = C ε t (U) -M I D C ε t (U) = C ε t -M I D ε t with M I D ε t = 1 h ∫ h/2 -h/2 ε t (U( z )) d z . ( 17
)
Since the material is homogeneous, the constitutive relation is constant through the thickness and, therefore, C can be taken out of the averaging operator M I D ( •) in Equation [START_REF] Wriggers | A note on enhanced strain methods for large deformations[END_REF]. With these assumptions, one can modify 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 -Non linear applications[END_REF] in order to obtain the weak formulation of the mixed element: find U ∈ U and σ ∈ S such that:

∫ Ω ε t * : C ε t -M I D ε t * : C M I D ε t dΩ + ∫ Ω M I D ε t * : σdΩ + ∫ Ω σ * : M I D ε t dΩ - ∫ Ω σ * : C -1 σdΩ = l (U * ), U * ∈ V and σ * ∈ S. (18)

| The discrete version

Now let us define consistent approximation spaces for the two unknowns U and Û σ of the problem. Following the principle of isoparametric elements, the functions N A , A ∈ {1, 2, .., n } of Section 2.2 are used to interpolate the displacement. One uses degree 2 along χ, which is sufficient to obtain quasi-optimal accuracy in the case of slender structures. For the loading part, one uses one polynomial degree less for the lengthwise functions. Since Û σ is constant through the thickness, constant functions of χ suffice. In summary, if, for the approximation in displacement U h , one takes the space Q p,q,2 which corresponds to the functions of polynomial degree p in ξ, q in η and 2 in χ, one ends up with the functions ÑC , C ∈ {1, 2, .., ñ } which generate the space Q p-1,q -1,0 which is used to interpolate the approximation of the stress part σh (see [START_REF] Elguedj | B-bar an F-bar projection methods for nearly incompressible linear and non linear elasticity and plasticity using higher order NURBS element[END_REF][START_REF] Bouclier | Locking free isogeometric formulations of curved thick beams[END_REF] for examples of the construction of such spaces).

Thus, the approximation U h is sought in the form:

U h = n A=1 N A (ξ, η, χ)U A i.e.            u h v h w h            = n A=1 N A            u A v A w A            , ( 19
)
where

U A = u A v A w A T is the vector of the control variables of U h = u h v h w h T at control point A. Ex- pression (19
) can be written in matrix form as:

U h = n A=1 N A U A with N A =         N A 0 0 0 N A 0 0 0 N A         . ( 20
)
Based on that, one can express the strains according to [START_REF] Hosseini | An isogeometric solid-like shell element for nonlinear analysis[END_REF] and derive the approximation ε h x of ε x :

ε h x = n A=1 B L A U A where B L A =                 N A,x 0 0 0 N A,y 0 0 0 N A,z 1 √ 2 N A,y 1 √ 2 N A,x 0 0 1 √ 2 N A,z 1 √ 2 N A,y 1 √ 2 N A,z 0 1 √ 2 N A,x                 . ( 21
)
B L A is the classical strain-displacement matrix at control point A. Regarding approximation σh , its expression is:

σh = ñ C =1 Ñ C σC t ex t r mw i t h Ñ C =                 ÑC 0 0 0 0 0 0 ÑC 0 0 0 0 0 0 ÑC 0 0 0 0 0 0 √ 2 ÑC 0 0 0 0 0 0 √ 2 ÑC 0 0 0 0 0 0 √ 2 ÑC                 . ( 22
)
σC = σC 19) and σC , C ∈ {1, 2, .., ñ } of ( 22) respectively, one gets the discrete version of Problem [START_REF] Reese | A large deformation solid-shell concept based on reduced integration with hourglass stabilization[END_REF]. Remembering that the local coordinate system is constant through the shell's thickness:

MID ε h t = MID R ε h x = R MID ε h x , (23)
and that the material is isotropic:

D = R T D R = R D R T , (24)
one arrives at the following linear system to be solved:

      K L uu K L σu T K L σu -K L σ σ              {U } { σ }        =        {F } {0}        , (25)
with, depending on the various operators introduced:

K L uu = ∫ Ω       B L T D B L - 1 h ∫ h/2 -h/2 B L ( z )d z T D 1 h ∫ h/2 -h/2 B L ( z )d z       dΩ ; K L σu = ∫ Ω Ñ T R 1 h ∫ h/2 -h/2 B L ( z )d z dΩ ; K L σ σ = ∫ Ω Ñ T D -1 Ñ dΩ ; {F } = ∫ Ω N T fdΩ + ∫ Γ F N T FdΓ F . (26)
Global operators B L , N and Ñ contain respectively each of the nodal contributions

B L A , N A , A ∈ {1, 2, .., n } and Ñ C , C ∈ {1, 2, .
., ñ }. The numerical integrations are carried out according to the classical strategy for NURBS. One uses standard Gauss quadrature with (p + 1) * (q + 1) * (2 + 1) integration points. This is the same standard Gauss integration scheme as that which would be used in the case of the full integration of classical finite elements. Recently, different and more optimal strategies have been proposed in this case of NURBS functions. One can mention the contributions of Hughes et al. [START_REF] Hughes | Efficient quadrature for NURBS-based isogeometric analysis[END_REF] and Auricchio et al. [START_REF] Auricchio | A Simple algorithm for obtaining nearly optimal quadrature rules for NURBS-based isogeometric analysis[END_REF]. These more sophisticated integration procedures were not considered in our work.

In practice, on this level, one carries out a static condensation of the stress unknowns over the whole structure in order to reduce System [START_REF] Hughes | A simple and efficient finite element for plate bending[END_REF] to the kinematic unknowns alone, leading to a more convenient form of the mixed element. The expression of this static condensation is:

K L σu {U } -K L σ σ { σ } = {0} . (27)
which leads to the equilibrium as a function of the displacement alone:

K L mi x ed {U } = {F } with K L mi x ed = K L uu + K L σu T K L σ σ -1 K L σu . (28)
Operator K L σ σ is symmetric and positive definite. Therefore, it is invertible. K L mi x ed is the global stiffness matrix of the mixed solid-shell element. This is a fully-populated matrix because of the inversion of K L σ σ .

| STABILIZATION OF THE ELEMENT UNDER SMALL DEFORMATIONS

In fact, one disadvantage of the mixed element of [START_REF] Bouclier | Efficient isogeometric NURBS-based solid-shell elements: Mixed formulation and B-method[END_REF] presented above is that it possesses hourglass modes. Among the first works attempting to address these hourglass modes, one can mention that of Belytschko et al. [START_REF] Flanangan | A Uniform strain hexahedron and quadrilateral with orthogonal hourglass control[END_REF][START_REF] Belytschko | Hourglass control in linear and nonlinear problems[END_REF].

One should bear in mind that, still from a classical finite element standpoint, dealing with such spurious modes is a recurring issue in the development of high-performance solid-shell elements. Therefore, effective techniques for controlling these modes have been proposed (see, for example [START_REF] Legay | Elastoplastic stability of shells using the physically stabilized finite element SHB8PS[END_REF][START_REF] Abed-Meraim | An improved assumed strain solid-shell element formulation with physical stabilization for geometric nonlinear applications and elastic-plastic stability 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 II -Non linear applications[END_REF][START_REF] Reese | A large deformation solid-shell concept based on reduced integration with hourglass stabilization[END_REF][START_REF] Bassa | A new nine-node solid-shell finite element using complete 3D constitutive laws[END_REF]). The same situation arises in the NURBS context. Far fewer works exist dealing with hourglass control in the NURBS context, but one could mention Bouclier et al. [START_REF] Bouclier | Locking free isogeometric formulations of curved thick beams[END_REF], which controlled hourglass modes resulting from a reduced integration of NURBS functions.

Hourglass modes are modes whose energy is nearly zero and which do not correspond to rigid body movements.

In a classical (e.g. displacement-based) formulation using a sufficient number of integration points, these modes would not occur because they would generate a nonzero energy. In our case, it is because of the mixed formulation and the choice of the approximation spaces that these modes appear. Their presence is also linked to the displacement boundary conditions which are applied to the shell problem. In Bouclier et al. [START_REF] Bouclier | Efficient isogeometric NURBS-based solid-shell elements: Mixed formulation and B-method[END_REF], the test cases considered had boundary conditions which prevented the development of such modes. The extension of the mixed element to geometrically nonlinear cases has led to the study of problems with different displacement boundary conditions, and it is in this context that the occurrence of these hourglass modes was observed. In order to control these modes, we propose to use the same approach as with classical finite elements, which consists in introducing an artificial stabilizing stiffness. This is done on order 2 lengthwise (p = q = 2). In fact, as noted in [START_REF] Bouclier | Efficient isogeometric NURBS-based solid-shell elements: Mixed formulation and B-method[END_REF], the quadratic restriction of the mixed element seems to be sufficient to calculate most of the usual shell problems accurately. The resulting stabilized mixed quadratic element enables one to deal with any problem under small deformations effectively, regardless of the displacement boundary conditions considered.

| Identification of the hourglass modes

In order to control hourglass modes, the first step is to identify them. More precisely, one must seek the strain field or the stress field of the hourglass mode in order to, subsequently, derive an associated energy which can be used to eliminate it. In order to do that, we studied the simplified problem of a plate of slenderness 10 (see Figure 2). The lengths L were equal to 100 and the thickness h was equal to 10. For the material parameters, we chose a usual configuration (E = 200000 and ν = 0.3). We discretized the plate using 8 elements per side. The objective was only to calculate and study the stiffness matrix of the problem. No boundary condition was introduced. That plate geometry was chosen because it is easy to visualize. It was obtained through linear mapping in the three spatial directions of the elements in the isoparametric domain (ξ, η, χ) (see again Figure 2). Thus, the hourglass modes had the same shapes (apart from the scaling coefficients) in the physical domain (x, y , z ) and in the isoparametric domain.

F I G U R E 2 The plate problem used to visualize the hourglass modes

We calculated the condensed global stiffness matrix of the plate problem using Equation ( 28) and studied its eigenvalues, which are shown in Table 1 in order of increasing absolute values. As expected, the first six eigenvalues (numbered 1 to 6) are very small compared to the others. The corresponding modes represent the rigid body motions which, in practice, would be eliminated by the displacement boundary conditions. Then, following these rigid body modes, one can observe four eigenvalues (numbered 7, .., 10) whose absolute values are small compared to the rest.

These are the eigenvalues associated with the hourglass modes. Their presence implies the existence of a rank deficiency in the stiffness matrix which generated an instability of the system (28) to be solved. In other words, this means that in the mixed formulation proposed too many constraints were released. The constraints which are responsible for the locking effects were reduced, but, at the same time, additional constraints were also eliminated, resulting in the rank problem affecting the stiffness matrix. Figure 3 shows the deformed configurations of the four hourglass modes in the physical domain.

F I G U R E 3
The deformed configurations of the hourglass modes for the plate problem One can observe oscillations in the mode shapes. This is a characteristic of hourglass modes. If these oscillations are not controlled, the deformed shape of the shell problem also presents "bumps" (see Section 5 "Numerical results", and especially Subsection 5.1). At this point, it is important to note that these hourglass modes in the IGA context have the same regularity as the NURBS functions. In this case, the displacement field is quadratic lengthwise and has C 1 continuity at the interior knots. Such a regularity can be observed in the deformed shapes (Figure 3). Consequently, the resulting strain and stress are linear within each element and C 0 at the element's edges. Therefore, controlling one of these hourglass modes in a single mesh element makes it disappear from the whole structure. The hourglass control is carried over, so to speak, from one element to the next by continuity. For further details concerning this point, the interested reader is urged to refer to Bouclier et al. [START_REF] Bouclier | Locking free isogeometric formulations of curved thick beams[END_REF].

Remark One can also note that due to the shape of the hourglass modes these can also disappear as a result of the displacement boundary conditions. For example, one can see that clamping one edge of the plate would be sufficient to eliminate them. By continuity, the hourglass control would be carried over from one edge to the next.

Thus, in order to control the spurious modes, it suffices to identify their shapes within an element. Here, we prefer to take a central element (see Figure 4), i.e. an element with four C 1 edges in terms of displacement continuity.

The modes in such an element are simpler than in a border element. In addition, in order to improve the accuracy of the calculations, one tends to perform the hourglass control in an element which is as close as possible to the center of the structure. In other words, if the mesh has n el 1 elements lengthwise along e 1 and n el 2 elements lengthwise along e 2 , one performs the hourglass control in the element whose position is n o hg :

n o hg = n el 1 × E n el 2 2 + E n el 1 2 + 1, (29)
where E( •) denotes the integer part. Thus, in our case, we are looking at the representation of the modes in Element 37. Describing the modes based on their deformed shape, i.e. describing the associated displacement field, is complicated. In fact, recalling that what is actually being sought is the mode's stress field or strain field, it is easier to take advantage of the mixed formulation and obtain the stress directly. In order to do that, one would consider the uncondensed stiffness matrix of Element 37 and calculate the stress unknowns σe hg i associated with the i th hourglass mode. Since the corresponding stress field, denoted σe h hg i

, is linear lengthwise and constant through the thickness [START_REF] Bucalem | Higher-order MITC general shell elements[END_REF], it can be defined very easily everywhere within the element. Reverting back to the element's parent domain ( ξ, η, χ) ∈ [-1; 1] 3 in which the Gaussian integration is to be carried out (see again Figure 4), one can describe the four spurious modes in terms of stresses analytically as follows:

σe h hg 1 =                            0 0 0 0 η ξ                            ; σe h hg 2 =                            - ξ ξ 0 - η 0 0                            ; σe h hg 3 =                            - η η 0 ξ 0 0                            ; σe h hg 4 =                            1 -1 0 0 0 0                            . ( 30
)
Hourglass modes 1, 2, 3 and 4 correspond respectively to eigenmodes 7, 8, 9 and 10 of Table 1. These stress fields, which seem to be consistent with the deformed shapes in Figure 3, will be used for hourglass control.

| Control of the hourglass modes

Now that we have the stress field associated with the hourglass mode in the central element of the structure, we must, in order to control that mode, calculate an elementary stabilization stiffness K L e st ab i which produces energy when this i th hourglass mode is activated. One can see that the modes given by ( 30) produce energy for a classical

∫ Ω ε h * t T D ε h t dΩ = ∫ Ω σ h * t T D -1 σ h t dΩ. ( 31
)
Using that bilinear form and the approach of Bouclier et al. [START_REF] Bouclier | Locking free isogeometric formulations of curved thick beams[END_REF], one expresses the elementary stabilization stiffness in the form:

K L e st ab i = K L e u a i K L e a i a i -1 K L e a i u , (32)
where, using the notations of Section 2, the various operators are given by:

K L e a i a i = ∫ Ωe ( σe h hg i ) T D -1 σe h hg i dΩ e ; K L e a i u = ∫ Ωe ( σe h hg i ) T R B L dΩ e ; K L e u a i = K L e a i u T . (33)
Ω e represents the elementary domain of the central element of the structure which was chosen for the hourglass control. Then, these elementary stabilization stiffnesses are added to the global stiffness of System [START_REF] Beirão Da Veiga | Avoiding shear locking for the Timoshenko beam problem via isogeometric collocation methods[END_REF]. Thus, the expression of the stabilized equilibrium becomes:

K L st ab mi x ed {U } = {F } with K L st ab mi x ed = K L mi x ed + 4 i =1 γ i K L e st ab i . ( 34
)
(γ i ) i ∈{1,2,3,4}
are real coefficients which are the penalty parameters of the stabilization. They can be chosen so that they balance the amplitudes of the hourglass modes to be controlled. The amplitude a i of the i th mode can be inferred from [START_REF] Cottrell | Isogeometric analysis: Toward Integration of CAD and FEA[END_REF] and is equal to:

{a i } = K L e a i a i -1 K L e a i u {U } , (35)
where {U } is the solution of the stabilized equilibrium [START_REF] Auricchio | A Simple algorithm for obtaining nearly optimal quadrature rules for NURBS-based isogeometric analysis[END_REF]. In practice, it has been observed that a wide range of values of the penalty parameters can result in an effective control of the hourglass modes: for example, for the plate problem considered, proper stabilization is obtained with γ i ∈ 10 -3 , 10 i ∈ {1, 2, 3, 4}). For the sake of simplicity, we usually take these coefficients equal to 1. The global stiffness matrix K L st ab mi x ed of the stabilized mixed solid-shell element expressed according to [START_REF] Auricchio | A Simple algorithm for obtaining nearly optimal quadrature rules for NURBS-based isogeometric analysis[END_REF] has correct rank. Table 2 shows the eigenvalues of K L st ab mi x ed for the initial plate problem. One can see that only the eigenvalues associated with the hourglass modes are different. Now they are positive and sufficiently large to generate some energy. The element resulting from this operation is no longer subject to hourglass modes.

Remark Since the hourglass modes have been identified on a central element of the structure (element with 4 C 1 edges), the control is consistent for meshes presenting such a central element, i.e. composed of at least 3 elements per side. However, we have been able to notice through the numerical experiments of section 5 that the control happens to be effective most of the time for coarser meshes. Thus, we propose to keep the strategy for all mesh sizes : the hourglass control is performed in the central element (Equation ( 29)).

Remark Since, in the case of NURBS, hourglass control is necessary only in a single element, the stabilization process is inexpensive. This is another advantage compared to classical finite elements based on Lagrange-type polynomials. Indeed, since these Lagrange functions are C 0 along element edges, the strains and stresses are C -1 and, therefore, the hourglass control must be performed in every element of the mesh, which increases considerably the number of operations required to stabilize the stiffness matrix.

Remark It is important to note that in our case of spurious modes reduced integration is never envisaged. The integration scheme chosen is the classical NURBS scheme as presented in Section 2.3. Thus, the modes expressed according to [START_REF] Bouclier | Locking free isogeometric formulations of curved thick beams[END_REF] produce energy with bilinear form [START_REF] Bouclier | On the development of NURBS-based isogeometric solid shell elements: 2D problems and preliminary extension to 3D[END_REF].

| EXTENSION OF THE ELEMENT TO GEOMETRICALLY NONLINEAR PROB-LEMS

Let us now consider geometrically nonlinear problems. The shell remains elastic, but is subjected to large rotations and large displacements which induce major changes in configuration Ω during the calculations. Therefore, the configuration can no longer be assumed to be identical to the reference configuration, now denoted Ω 0 . In this section, we propose an extension of the mixed solid-shell NURBS element of Section 2 which addresses this geometric nonlinearity. In order to do that, we use a total Lagrangian strategy, which means that all the calculations are performed based on the reference configuration Ω 0 . Initially, in Sections 4.1, 4.2 and 4.3, we will not attempt to control the hourglass modes, but only express the equilibrium of the structure in the geometrically nonlinear case and then study its stability, leading to the element's tangent stiffness matrix and finally to the algorithm for the resolution of the nonlinear problem. It is only afterward, in Section 4.4, that the principles of Section 3 will be adapted to the new situation in order to control the hourglass modes in the geometrically nonlinear context.

| Equilibrium equations in the geometrically nonlinear case

| The continuous version

The same technique that was presented in Section 2.3 can be easily reused in the geometrically nonlinear case to extend the weak formulation of the mixed problem [START_REF] Reese | A large deformation solid-shell concept based on reduced integration with hourglass stabilization[END_REF] to the geometrically nonlinear case. Throughout this section, we use the same notations as in Section 2. Thus, the mixed variational formulation which expresses the geometrically nonlinear equilibrium becomes: find U ∈ U and Π ∈ S such that:

∫ Ω 0 E t * : C E t -M I D E t * : C M I D E t dΩ 0 + ∫ Ω 0 M I D E t * : ΠdΩ 0 + ∫ Ω 0 Π * : M I D E t dΩ 0 - ∫ Ω 0 Π * : C -1 ΠdΩ 0 = l (U * ), U * ∈ V and Π * ∈ S. (36)
As in the small-deformation case, the stresses and strains are expressed in the shell's local coordinate system. Since a total Lagrangian strategy is used, the local system (e 1 , e 2 , e 3 ) is the one of the reference initial undeformed configuration, i.e. that it is constant along the nonlinear simulation and it is the same as the one in the linear case. ε t and σ are replaced by E t and Π respectively. E t is the Green-Lagrange strain tensor and Π is the average through the thickness of the second Piola-Kirchhoff stress tensor. The work of the external loads l (U * ) is calculated as in [START_REF] Dornisch | Isogeometric Reissner-Mindlin shell analysis with exactly calculated director vectors[END_REF], but on the reference geometry Ω 0 . In fact, everything turns out to be very similar to the previous small-deformation case. ∫ Ω has simply become ∫ Ω 0

, except that ε t and σ were replaced respectively by E t and Π. This is not surprising because the behavior between σ t and ε t in Section 2 was the same as the behavior between Π t and E t in Ω 0 now:

Π t = C E t .
The Green-Lagrange strain tensor in the local basis (e 1 , e 2 , e 3 ) is expressed as a function of the displacement field as follows:

E t = R E x with E x = ε x (U) + 1 2 ε x Q (U, U) . ( 37
)
ε x is the linearized strain corresponding to small deformations as in [START_REF] Echter | A hierarchic family of isogeometric shell finite elements[END_REF], and ε x Q (U, U) is the quadratic strain, which is a function of the displacement gradient in (x, y, z):

ε x Q (U, U) = x U T x U or, in indexed form, ε Q i j (U, U) = 3 k =1 ∂U k ∂x i ∂U k ∂x j . ( 38
)
As for the virtual strain E t * , its expression is:

E t * = R E x * with E x * = ε x U * + ε x Q U, U * . (39)

| The discrete version

As previously in Section 2, one seeks an approximation of the displacement solution U h in Q p,q,2 and an approximation of the loading part Û ˝h in Q p-1,q -1,0 . The discrete expression is:

E h t (U h ) = R B L + 1 2 B N L (U h ) {U } , E h t * = R B L + B N L (U h ) U * and Πh = Ñ Π . (40)
Here, operator B N L (U h ) was introduced to express Relation [START_REF] Rogers | An introduction to NURBS With Historical Perspective[END_REF] in discrete form. The mixed system to be solved is obtained by substitution of the approximations (40) into the variational expression [START_REF] Piegl | The NURBS Book (Monographs in Visual Communication)[END_REF] and takes the form:

       K N L uu U h K N L u π U h K N L πu U h -K N L π π               {U } Π        =        {F } {0}        , (41)
with the following expressions of the various operators:

K N L uu U h = ∫ Ω 0 B L + B N L (U h ) T D B L + 1 2 B N L (U h ) dΩ 0 - ∫ Ω 0 1 h ∫ h/2 -h/2 B L + B N L (U h ) d z T D 1 h ∫ h/2 -h/2 B L + 1 2 B N L (U h ) d z dΩ 0 ; K N L u π U h = ∫ Ω 0 1 h ∫ h/2 -h/2 B L + B N L (U h ) d z T R T Ñ dΩ 0 ; K N L πu U h = ∫ Ω 0 Ñ T R 1 h ∫ h/2 -h/2 B L + 1 2 B N L (U h ) d z dΩ 0 ; K N L π π = ∫ Ω 0 Ñ T D -1 Ñ dΩ 0 . (42)
Then, as in [START_REF] Belytschko | Hourglass control in linear and nonlinear problems[END_REF], a static condensation is carried out:

K N L πu U h {U } -K N L π π Π = {0} . ( 43
)
which, eventually, leads to the following geometrically nonlinear global stiffness matrix:

K N L mi x ed U h = K N L uu U h + K N L u π U h K N L π π -1 K N L πu U h . (44)
This matrix, as was the case under small deformations [START_REF] Beirão Da Veiga | Avoiding shear locking for the Timoshenko beam problem via isogeometric collocation methods[END_REF], is a fully-populated matrix because of the inversion of operator K N L π π . Finally, the discrete equilibrium can be expressed as a function of the displacement alone:

K N L mi x ed U h {U } = {F } . ( 45
)
4.2 | Derivation of the mixed tangent matrix : stability study in the geometrically nonlinear case

| The continuous version

We now study the stability of the structure in order to construct the tangent matrix of the mixed solid-shell element.

To study the stability, one must test the existence of a second displacement field belonging to U and a second stress field belonging to S which satisfy the mixed equilibrium [START_REF] Piegl | The NURBS Book (Monographs in Visual Communication)[END_REF]. This would mean that the structure has reached an unstable equilibrium state and is likely to buckle, as illustrated in Figure 5. More precisely, let us assume that U and Π satisfy the equilibrium equation ( 36) and then seek other equilibrium solutions of the form U+αU 1 (with the associated test field U * ) and Π + α Π1 (with the associated test field Π * ), α being a real coefficient which tends toward zero. In this case, one can observe that U 1 belongs to V. For this second displacement field, the strain can be rewritten as: F I G U R E 5 Buckling : the three configurations Then, substituting these expressions [START_REF] Gruttman | Theory and numerics of thin elastic shells with finite rotations[END_REF] into variational form [START_REF] Piegl | The NURBS Book (Monographs in Visual Communication)[END_REF], one gets the expression of the geometrically nonlinear equilibrium for these second fields. The first terms (which are constant with respect to α) cancel each other out because U and Π are already solutions of Equation [START_REF] Piegl | The NURBS Book (Monographs in Visual Communication)[END_REF]. There remain the terms in α and α 2 . For a sufficiently small α, it suffices that the α term be zero for the mixed equilibrium to be satisfied for U + αU 1 and Π + α Π1 . Thus, the stability of equilibrium U, Π boils down to: given U ∈ U and Π ∈ S, find U 1 ∈ V and Π1 ∈ S such that:

         E t → R E x (U + αU 1 ) with E x (U + αU 1 ) = E x (U) + α ε x (U 1 ) + ε x Q (U, U 1 ) + α 2 2 ε x Q (U 1 , U 1 ) ; E t * → R ε x (U * ) + ε x Q (U, U * ) + α ε x Q (U 1 , U * ) . (46)
∫ Ω 0 ε x Q U 1 , U * : C ε x (U) + 1 2 ε x Q (U, U) dΩ 0 + ∫ Ω 0 ε x U * + ε x Q U, U * : C ε x (U 1 ) + ε x Q (U, U 1 ) dΩ 0 - ∫ Ω 0 M I D ε x Q U 1 , U * : C M I D ε x (U) + 1 2 ε x Q (U, U) dΩ 0 + ∫ Ω 0 M I D ε x U * + ε x Q U, U * : C M I D ε x (U 1 ) + ε x Q (U, U 1 ) dΩ 0 + ∫ Ω 0 R M I D ε x U * + ε x Q U, U * : Π1 dΩ 0 ∫ Ω 0 R M I D ε x Q U 1 , U * : ΠdΩ 0 + ∫ Ω 0 Π * : R M I D ε x (U 1 ) + ε x Q (U, U 1 ) dΩ 0 - ∫ Ω 0 Π * : C -1 Π1 dΩ 0 = 0, U * ∈ V and Π * ∈ S. (47)
This formulation will be used later to derive the tangent matrix of the mixed solid-shell element. Right now, one can already identify the terms which participate in the geometric stiffness of the structure. Because of the fact that in our element the stress is divided into two parts, the geometric stiffness appears more precisely to be composed of two terms: the double-bracket term { {•} } which characterizes the geometric stiffness of the average stress through the thickness, and the single-bracket term {•} which is related to the geometric stiffness of the orthogonal part of the average stress through the thickness.

| The discrete version

In practice, U h 1 and Πh 1 are sought in the same spaces as U h and Πh :

U h 1 = N {U 1 } ; Πh 1 = Ñ Π1 . ( 48
)
The discrete version of Equation [START_REF] Sze | Popular Benchmark problems for geometric nonlinear analysis of shells[END_REF] through which these can be determined leads to the resolution of a system of the form:

       K 1 N L Tuu U h + K 2 N L Tuu U h + K 3 N L Tuu Πh K N L T u π U h K N L T πu U h -K N L T π π               {U 1 } Π1        =        {0} {0}        . ( 49
)
where the classical operators take the following expressions:

K 1 N L Tuu U h = ∫ Ω 0 B L + B N L U h T D B L + B N L U h dΩ 0 - ∫ Ω 0 1 h ∫ h/2 -h/2 B L + B N L U h T D 1 h ∫ h/2 -h/2 B L + B N L U h dΩ 0 ; K N L T u π U h = K N L T π u U h T = ∫ Ω 0 1 h ∫ h/2 -h/2 B L + B N L U h T R T Ñ dΩ 0 ; K N L T π π = ∫ Ω 0 Ñ T D -1 Ñ dΩ 0 , (50)
and the geometric stiffnesses are obtained as:

U * T K 2 N L Tuu U h {U 1 } = ∫ Ω 0 ε x Q U h 1 , U h * : C ε x U h + 1 2 ε x Q U h , U h dΩ 0 - ∫ Ω 0 M I D ε x Q U h 1 , U h * : C M I D ε x U h + 1 2 ε x Q U h , U h dΩ 0 ; U * T K 3 N L Tuu Πh {U 1 } = ∫ Ω 0 R M I D ε x Q U h 1 , U h * : Πh dΩ 0 . (51)
The construction of these geometric stiffness matrices follows the same principle as for classical finite elements (see, for example, Legay et al. [START_REF] Legay | Elastoplastic stability of shells using the physically stabilized finite element SHB8PS[END_REF]). Then a static condensation is carried out, leading to the global tangent matrix:

K N L T mi x ed U h , Πh = K 1 N L Tuu U h + K 2 N L Tuu U h + K 3 N L Tuu Πh + K N L T u π U h K N L T π π -1 K N L T πu U h . (52)
Again, this matrix is fully populated due to the calculation of

K N L T π π -1
.

Finally, one can see that if there is a state such that the tangent matrix defined by Equation ( 52) is singular this means that we found a nonzero displacement field U 1 which satisfies Equation [START_REF] Sze | Popular Benchmark problems for geometric nonlinear analysis of shells[END_REF]. This field corresponds to the buckling mode. The structure has reached an unstable equilibrium condition.

| Calculation of the geometrically nonlinear equilibrium of the structure

In order to solve the discrete geometrically nonlinear equilibrium equation [START_REF] Buechter | Shell theory versus degeneration -a comparison in large rotation finite element analysis[END_REF], one sets up a classical Newton algorithm in which the tangent matrix is updated at each iteration. Because of the static condensations of the various operators, one can follow the usual finite element approach: starting from an undeformed state F ext 1 = {F 1 } = {0} i.e. U h 1 = 0 and Π1 one seeks the state which corresponds to the maximum loading F ext max . One defines a loading step ˊF, a target accuracy and a maximum number of iterations nb max , then one applies the procedure described in Algorithm 1:

| Stabilization: control of the hourglass modes in the geometrically nonlinear case

There remains the question of the stabilization of the mixed solid-shell element in the geometrically nonlinear case. Let us proceed as in Section 3. Here again, we limit ourselves to functions which are quadratic lengthwise (p = q = 2).This stabilization will result in the adjunction of artificial elementary stabilization stiffnesses to the equilibrium stiffness matrix as in Equation [START_REF] Simo | On stress resultant geometrically exact shell model, Part III: computational aspects of nonlinear theory[END_REF] and to the tangent matrix as in Equation (52). In order to carry out the geometrically

| Stabilization of the equilibrium

Under large rotations and large displacements, the classical bilinear form (31) used to produce energy from the spurious modes becomes:

∫ Ω 0 (E h * t ) T D E h t dΩ 0 = ∫ Ω 0 (Π h * t ) T D -1 Π h t dΩ 0 . (54)
With that expression, repeating the approach of Section 3.2, the elementary stabilization stiffness for the i th hourglass mode is calculated as:

K N L e st ab i U h = K N L e u a i U h K N L e a i a i -1 K N L e a i u U h , (55)
with

K N L e u a i U h = ∫ Ω 0e B L + B N L U h T R T Πe h hg i dΩ 0e ; K N L e a i u U h = ∫ Ω 0e ( Πe h hg i ) T R B L + 1 2 B N L U h dΩ 0e ; K N L e a i a i = ∫ Ω 0e ( Πe h hg i ) T D -1 Πe h hg i dΩ 0e . (56)
Ω 0e represents the elementary domain in the reference configuration of the central element chosen for the hourglass control. The stabilized geometrically nonlinear equilibrium is expressed by:

K N L st ab mi x ed U h {U } = {F } with K N L st ab mi x ed U h = K N L mi x ed U h + 4 i =1 γ i K N L e st ab i U h , (57)
where K N L mi x ed U h is the unstabilized stiffness matrix of Equation [START_REF] Simo | On stress resultant geometrically exact shell model, Part III: computational aspects of nonlinear theory[END_REF]. Thus, in the resolution algorithm of Section 4.3, K N L st ab mi x ed U h replaces K N L mi x ed U h for the calculation of the internal loads. The amplitudes of the hourglass modes are also calculated in that step in the same way as in Equation [START_REF] Cohen | Discrete B-spline and subdivision techniques in computer aided geometric design and computer graphics[END_REF]:

{a i } = K N L e u a i U h K N L e a i a i -1 {U } , i ∈ {1, 2, 3, 4} , (58)
where {U } is the solution of the stabilized equilibrium (57).

| Stabilization of the tangent matrix

In order to obtain the elementary stabilization stiffness of the tangent matrix, one goes back to the expression of the stability of the geometrically nonlinear equilibrium (see the details in Section 4.2). In addition to the new displacement and stress fields of the forms U + αU 1 (with the associated test field U * ) and Π + α Π1 (with the associated test field Π * ), one verifies the existence of new hourglass modes whose elementary stresses are of the form Πe , i ∈ {1, 2, 3, 4} are already solutions of the stabilized equilibrium), one ends up with an elementary stabilization stiffness matrix for the i th mode of the form:

hg i + α Πe 1 hg i , i ∈ {1, 2,
K N L e T st ab i U h , Πe h hg i = K N L e Tuu Πe h hg i + K N L e Tua i U h K N L e
Ta i a i -1

K N L e Ta i u U h . ( 59
)
The expressions of the classical operators are:

K N L e Tua i U h = K N L e Ta i u U h T = ∫ Ω 0e B L + B N L U h T R T Πe h hg i dΩ 0e ; K N L e Ta i a i = ∫ Ω 0e ( Πe h hg i ) T D -1 Πe h hg i dΩ 0e , (60)
and the single-bracket term {•}, which corresponds to a geometric hourglass stiffness, is calculated as:

U * T K N L e Tuu Πe h hg i {U 1 } = ∫ Ω 0e R ε x Q U h 1 , U h * : Πe h hg i dΩ 0e . (61)
Finally, the stabilized tangent matrix is given as a function of the unstabilized tangent matrix (52) by:

K N L st ab T mi x ed U h , Πh , Πe h hg i = K N L T mi x ed U h , Πh + 4 i =1 γ i K N L e T st ab i U h , Πe h hg i . ( 62
)
This is the matrix which was built in the resolution algorithm of Section 5.2. With this modification and the modification of the equilibrium, the solid-shell mixed element becomes stable and can be used to solve any geometrically nonlinear problem.

| NUMERICAL RESULTS

In this section, in order to assess the performance of the stabilized mixed solid-shell element developed in the previous sections, we present the results of several common test cases. The objective of these tests is to evaluate the robustness and the accuracy of the element both in the linear case and in the geometrically nonlinear case. This presentation is divided into two parts: first, Section 5.1 presents the results of two linear problems; then, Section 5.2 reanalyzes, this time in the geometrically nonlinear context, the two previous test cases and considers an additional case.

Concerning the discretization, we started from a patch composed of a single element, to which we applied the k -refinement strategy. Thus, the continuity across the interior knots was C p-1 , p being the polynomial degree of the NURBS function. From here on, the mesh of the structure composed of N elements along the first length ξ and M elements along the second length η will be denoted N × M . Concerning the mixed element, we used the quadratic version p = q = 2 whose hourglass control was studied in Sections 3 and 4.4. This element will be referred to as "Mixed 2". Once the element has been subjected to hourglass control, "hg" will be added to its designation which will become "Mixed 2 hg".

Our goal is to show that this low-order mixed element is unaffected by locking and, therefore, is very accurate for the resolution of difficult problems. In order to do that, we compared our element to standard NURBS solid-shell elements of a higher order. According to Bouclier et al. [START_REF] Bouclier | Efficient isogeometric NURBS-based solid-shell elements: Mixed formulation and B-method[END_REF], the behavior of these standard solid-shell elements seems to be quite similar to that of standard Reissner-Mindlin or Kirchhoff-Love NURBS shell elements [START_REF] Kiendl | Isogeometric shell analysis with Kirchhoff-Love elements[END_REF][START_REF] Benson | Isogeometric shell analysis: the Reissner-Mindlin Shell[END_REF]. Therefore, such a comparison is a good test of the superiority of our approach over any standard NURBS shell element, be it a structural model or a solid model. Here, what is meant by "standard" is that the usual NURBS modeling approach is applied to the shell elements directly: a classical displacement-based formulation is assumed and no specific treatment in order to handle locking is considered. The standard solid-shell elements are denoted "basic", with the adjunction of "2" for quadratic elements (p = q = 2), "3" for cubic elements (p = q = 3), and "4" for quartic elements (p = q = 4). In addition, to guarantee the validity of our simulations, the solutions obtained are compared to reference solutions, denoted "Reference", which were taken from the finite element literature. Finally, since the main objective of the paper is to extend the mixed method of [START_REF] Bouclier | Efficient isogeometric NURBS-based solid-shell elements: Mixed formulation and B-method[END_REF] to the geometrically nonlinear case, the nonlinear results are compared to the results of other high-performance finite element techniques, both NURBS and classical, taken from the literature.

| The linear case: elasticity under small deformations

The two problems considered for the linear case were the free-edged pinched cylinder and the twisted beam under out-of-plane loading. These problems were initially proposed in MacNeal and Harder [START_REF] Macneal | A proposed standard set of problems to test finite element accuracy[END_REF] and, subsequently, were widely used in the literature to the point that it is now customary to use these test cases to validate the performance of elements. In particular, these problems enable one to test the locking behavior of elements. Reference solutions can be found in [START_REF] Macneal | A proposed standard set of problems to test finite element accuracy[END_REF]. These are the reference solutions we used in the studies reported below. The pinched cylinder problem possesses displacement boundary conditions which let hourglass modes of the mixed element develop, thus making its study useful to validate the stabilization method proposed in Section 3. In the case of the twisted beam, these spurious modes are eliminated by the displacement boundary conditions applied, so stabilization would be unnecessary. For a more in-depth study, the interested reader is encouraged to refer to Bouclier et al. [START_REF] Bouclier | Efficient isogeometric NURBS-based solid-shell elements: Mixed formulation and B-method[END_REF]. In that reference, other linear test cases in which the hourglass modes cannot develop were also studied.

| The free-edged pinched cylindrical shell problem

In this problem, a free-edged cylinder is subjected to two diametrally opposed concentrated loads (see Figure 6).

There is also a similar test case with diaphragms at the level of the cylinder's ends, but the free-edged version is more interesting in our case because it lets hourglass modes develop. Besides, the version with diaphragms has already been studied in [START_REF] Bouclier | Efficient isogeometric NURBS-based solid-shell elements: Mixed formulation and B-method[END_REF]. Due to the symmetry of the problem, the model is restricted to one-eighth of the cylinder. This problem involves significant bending energy, which makes its resolution difficult due to its propensity to locking. In the linear case, the cylinder is assumed to be pinched. A load P with a magnitude of 25 is applied downward at point A. The reference solution from [START_REF] Macneal | A proposed standard set of problems to test finite element accuracy[END_REF], defined as the displacement along z at the point of application of the load, is equal to -W r ef (A) = 0.1139. appears to be perfectly smooth. These oscillations are characteristic of the presence of hourglass modes (whose shapes were shown in Figure 3). Thus, for this problem, the stabilization method proposed in Section 3 seems to be effective. It enables one to get rid of these oscillations, which means that there are no more hourglass modes. In addition, in Figure 8 which shows the convergence of the displacement solution as a function of the mesh refinement, one can observe that these oscillations have an adverse effect on the displacement of interest. Indeed, due to these oscillations the Mixed 2 vertical displacement at point A (Figure 8(a)) converges much more slowly and is oddly stiffer than the corresponding Mixed 2 hg solution (Figure 8(b)). In this case, the latter converges as rapidly as the Basic 4 solution. This is due to the presence of locking in the basic element. Increasing the order of interpolation of the basic element reduces locking, but does not eliminate it completely as the mixed strategy does.

Remark

It is important to note that the reference solution of [START_REF] Macneal | A proposed standard set of problems to test finite element accuracy[END_REF] which is used here is provided from a pure shell model (of Kirchhoff-Love type). In our case of a full 3D model, the displacement under the load is expected to converge to infinity because of the singularity of the pointwise force. Nevertheless, for the coarse meshes of interest, the effect of the singularity is negligible that is why the study to assess locking remains valid. Such a comment also hold for the pinched hemispherical shell with a hole studied in section 5.2.

| The twisted beam problem under out-of-plane loading

The second problem considered is that of a helicoidal shell which is clamped at one end and subjected to an out-ofplane load at the other end (see Figure 9). The load is applied in our model using an out-of-plane uniform surface force on the end cross-section of the beam. This problem is well-known for the evaluation of element performance in the context of twisted structures. The distortion of the elements makes the test especially discriminating regarding locking phenomena. For the linear case, we considered two beam thicknesses: a thick beam with h = 0.32 and a thin beam with h = 0.0032. The associated reference solutions from [START_REF] Macneal | A proposed standard set of problems to test finite element accuracy[END_REF] in terms of the displacement V (A) for an applied load P = 1 are V r ef (A) = 0.00175 and V r ef (A) = 1296 respectively.

x y z P

U(A) V(A) W(A) A clamped (x,y,z) b L F I G U R E 9
The twisted beam: description and data of the problem The structure's mid-surface is defined by parameters s and θ as follows:

x = 2θL/π ; y = s cos θ ; z = -s sin θ ;

0 ≤ θ ≤ π/2 ; -b 2 ≤ s ≤ b 2 . (63)
Since this surface is helicoidal, it cannot be represented exactly using NURBS functions. Here, in order to calculate the problem, we approximated the geometry using a single NURBS element. Then, in strict compliance with the IGA approach, we kept the same geometry while refining the NURBS mesh. The construction of the approximate geometry is described in Annex A.

The analysis of the results follows the same pattern as in the first test case. First, Figure 10 shows the deformed configurations obtained with the Mixed 2 hg element for an 8 × 2 mesh. The thick beam is shown in Figure 10(a In this case, it has to be noted that the Mixed 2 element would have provided very similar results. The deformed configurations would have not been subjected to oscillations and the displacement plots would have been very close. This can be explained by the clamped end on the right-hand side of the beam which eliminates hourglass modes. Stabilization as proposed in Section 3 is actually not necessary here. Nevertheless, in order to show the ability of the proposed element to be really locking-free, whatever the displacement boundary conditions are, we decided to keep the stabilized element to perform this study. The important thing to observe on the solutions shown in Figure 11 is that the Mixed 2-hg solution does not seem to vary as a function of the slenderness of the beam (the shapes in Figures 11(a) and 11(b) are similar). This is not the case for the basic solutions which, for identical meshes, become much less accurate when the slenderness of the beam increases. For the thick beam, one can note again that Basic 4 is quite comparable to Mixed 2 hg, whereas for the thin beam Basic 4 is too stiff. This is a consequence of locking. An increase in the order of approximation enables one to reduce locking, but does not eliminate it altogether. The more slender the beam, the more pronounced this phenomenon.

(a) Thickness h=0.32 (scale factor 500) (b) Thickness h=0.0032 (scale factor 0.0015)

F I G U R E 1 0 Initial configuration (mesh) and deformed shape (displacements along y) using the stabilized mixed element for the twisted beam in elasticity under small deformations

| Geometrically nonlinear calculations

In this second section, we test the extension of the element proposed in Section 4 to large rotations and large displacements. In order to do that, the test cases considered were the nonlinear counterparts of the test cases of the previous section and the additional case of a pinched hemispherical shell with a hole.This time, the magnitudes of the loads were greater (see the values of P nonl i n max in Figures 6,9 and 12). For these calculations, we implemented the usual nonlinear incremental algorithm of Section 4.3 along with the hourglass control of Section 4.4 when necessary. The values used for the parameters of the algorithm were = 10 -5 (for the accuracy) and nb max = 25 (for the maximum number of iterations). The loading step ∆F was chosen, depending on the test case, such that the convergence of the algorithm would occur after less than 10 iterations. In addition to comparing the mixed element and higher-order basic elements, we were also interested in the reliability of the proposed strategy compared to other high-performance shell techniques available today. In the NURBS context, our comparison included the solid-like shell Bézier element (SLSBEZ) of Hosseini et al. [START_REF] Hosseini | An isogeometric solid-like shell element for nonlinear analysis[END_REF]. Indeed, even though that technique does not have as its main objective the elimination of locking, it uses a solid-like model and is claimed to be efficient for low orders. Today, this is the geometrically nonlinear NURBS shell technique which seems to be closest to ours. The corresponding cubic element will be denoted "SLSBEZ 3". For comparison, in the classical finite element context, we chose the SHB8PS element of Abed-Meraim and Combescure [START_REF] Abed-Meraim | An improved assumed strain solid-shell element formulation with physical stabilization for geometric nonlinear applications and elastic-plastic stability analysis[END_REF]. This element is an 8-node brick which, even though it involves a modified 3D constitutive relation, appears to be competitive against other traditional solid-shell finite element techniques. Since its interpolation functions are linear, this element will be denoted "SHB8PS 1".

| Case of a pinched hemispherical shell with a hole

The first test case is the hemispherical shell problem of Figure 12. In this problem, a thin hemispherical shell with a hole at its apex is subjected to concentrated loads at four points of its free base. Because of the symmetry of the problem, we modeled only one-fourth of the structure and, in order to eliminate the remaining rigid body movement, we fixed the topmost left corner of the one-eighth sphere (point C ) in the z direction. This is a very popular test case when it comes to verifying the absence of membrane locking and the good representation of rigid body modes (see, for example, among the earliest works, [START_REF] Simo | On stress resultant geometrically exact shell model, Part III: computational aspects of nonlinear theory[END_REF][START_REF] Buechter | Shell theory versus degeneration -a comparison in large rotation finite element analysis[END_REF]). The maximum load P max was set to 100. We used a load step ∆F = 5 for the calculations. Since hourglass modes are likely to occur in this test case, it was necessary to apply the hourglass control of Section 4.4 for the mixed element. Usually, throughout the calculation, one studies the loading curve P as a function of the displacement u at the point of application A and in the direction of the load. The corresponding solution given in Hosseini et al. [START_REF] Hosseini | An isogeometric solid-like shell element for nonlinear analysis[END_REF], based on the traditional S4R shell finite element of ABAQUS and a 16 × 16 mesh, was used as the reference for our comparisons. observe the same trends as in the small-deformation studies, which proves the good behavior of the element in the geometrically nonlinear context for this test case. Once again one can note that:

1. An increase in the order of the basic elements leads to an improvement in the solution, but for orders 3 and 4 a good approximation requires no less than 16 × 16 elements. The coarser meshes are too stiff because of the locking phenomenon;

2. The choice of a mixed strategy enabled us to get a good approximation for order 2 and a coarser mesh (4 × 4).

With that strategy, there appears to be no locking, which explains the mixed element's very good accuracy even with very coarse meshes.

Next, Figure 15 shows the solutions of element Mixed 2 hg compared to the solutions of element SLSBEZ 3. With SLSBEZ 3, a 16 × 16 mesh appears to be necessary to obtain a correct solution. Therefore, that is the mesh which should be compared to the 8 × 8 mesh (or even to the 4 × 4 mesh) using the Mixed 2 hg element. This comparison is given in Table 3, which shows the numbers of degrees of freedom (DOFs) and the numbers of Gauss points associated with these meshes. These results confirm the good performance of the Mixed 2 hg element developed here, since it appears to be more accurate for the same number of DOFs and the same number of Gauss points than element SLSBEZ 3. F I G U R E 1 5 Comparison of the Mixed 2 hg solutions with the SLSBEZ 3 solutions of Hosseini et al. [START_REF] Hosseini | An isogeometric solid-like shell element for nonlinear analysis[END_REF] for the pinched hemispherical shell

| Pull-out of an open-ended cylindrical shell

In the nonlinear version, the cylinder of Figure 6 is pulled out until a maximum load P max = 10, 000 is reached for one-eighth of the structure. With this loading and the boundary conditions, the cylinder is subjected to considerable rotations combining bending and membrane effects, which make the numerical resolution very difficult. Because of that, this test case has become very popular (let us mention [START_REF] Gruttman | Theory and numerics of thin elastic shells with finite rotations[END_REF][START_REF] Sze | Popular Benchmark problems for geometric nonlinear analysis of shells[END_REF] to cite only two references). For our calculation, we considered a loading step ∆F = 100. Since hourglass modes could appear, we had to use the stabilization method of Section 4.4 for the mixed element. The relations of the load P to the radial displacements W (A), -U (B) and -U (C ) during the calculation are interesting to follow. (See again Figure 6 for the positions of points A, B and C .) For the convergence study, we considered 8 × 8, 16 × 8 and 16 × 16 meshes (the first direction ξ being the circumferential direction and the second direction η being the longitudinal direction). Since element Mixed 2 hg seems to have reached convergence with 8 × 8 elements, the corresponding results are presented only for the coarsest mesh. These solutions are also compared to the reference solutions given in Sze et al. [START_REF] Sze | Popular Benchmark problems for geometric nonlinear analysis of shells[END_REF] using the traditional S4R element, particularly with 24 × 16 and 36 × 24 meshes. The deformed configurations obtained during the calculations with element Mixed 2 hg are shown in Figure 17 and the (normalized load)-(displacement) curves are shown in Figure 17. Again, no oscillations of the deformed shapes can be observed which accounts for the stabilization strategy proposed. With regard to the behavior of this test case, the results indicate the presence of two regimes: the first regime is dominated by bending effects and is characterized by large displacements and large rotations, whereas in the second regime mostly membrane effects are present. In addition, one can observe the occurrence of a "snap-through" effect once the loading reaches a critical value close to 5, 000. Then, the direction of the displacement at point C changes. Regarding the quality of the solution, the accuracy with Basic 4 and Mixed 2 hg appear to be similar. For this coarse mesh, these two solutions seem to be reasonably close to the reference solution of the problem.

Finally, in Figure 18, the results of element Mixed 2 hg are compared to those of elements SLSBEZ 3 and SHB8PS 1.

The convergence mesh for SLSBEZ 3 seems to be composed of 16 elements per side and, according to [START_REF] Abed-Meraim | An improved assumed strain solid-shell element formulation with physical stabilization for geometric nonlinear applications and elastic-plastic stability analysis[END_REF], the convergence with element SHB8PS 1 requires a 30 × 20 mesh. Table 4 gives the characteristics of these meshes along with those of the 8 × 8 mesh for element Mixed 2 hg. Again, one can observe that element Mixed 2 hg is attractive compared to element SLSBEZ 3. One can also note that the NURBS strategies perform better than the classical finite element strategy in terms of the number of DOFs, but are less interesting in terms of the number of Gauss points. This is a classical observation. Today, in IGA, the adverse effect of too many integration points in NURBS functions is well-known. The rule which consists in taking (p + 1) Gauss points per element for functions of degree p is far from optimal. A few techniques aiming to alleviate this problem have already been proposed. One could certainly use the more optimal integration rules of [START_REF] Hughes | Efficient quadrature for NURBS-based isogeometric analysis[END_REF][START_REF] Auricchio | A Simple algorithm for obtaining nearly optimal quadrature rules for NURBS-based isogeometric analysis[END_REF], or consider reduced-integration rules as proposed is the recent work of Schillinger et al. [START_REF] Schillinger | Reduced bezier element quadrature rules for quadratic and cubic splines in isogeometric analysis[END_REF], including hourglass control if necessary, following the formalism proposed in Sections 3 and 4.4.

For example, the computations for the mixed element have also been carried out here using a reduced Gaussian rule which consists of taking p Gauss points instead of (p + 1) for functions of degree p. This strategy has initially been investigated in [START_REF] Benson | Isogeometric shell analysis: the Reissner-Mindlin Shell[END_REF][START_REF] Benson | A large deformation, rotation-free, isogeometric shell[END_REF] before being considered again in [START_REF] Schillinger | Reduced bezier element quadrature rules for quadratic and cubic splines in isogeometric analysis[END_REF] where a different way to treat the boundary elements is pointed out. For the Mixed 2 hg meshes considered on the cylinder problem, no additional instabilities have been reported using 2 × 2 × 2 = 8 Gauss points per element everywhere in the mesh. The obtained solutions look identical to the Mixed 2 hg solutions with 27 Gauss points per element. This mixed element of degree 2 with reduced Gaussian integration is denoted "Mixed 2 hg reduced" in what follows. Its properties for the 8 × 8 mesh are added in Table 4.

We observe that such a strategy already enables to greatly reduce the numerical integration cost. In terms of Gauss points, the NURBS mesh becomes competitive with respect to the classical finite elements mesh. F I G U R E 1 8 Comparison of the Mixed 2 hg solutions with the SLSBEZ 3 solutions of Hosseini et al. [START_REF] Hosseini | An isogeometric solid-like shell element for nonlinear analysis[END_REF] and the SHB8PS 1 solutions of Abed-Meraim and Combescure [START_REF] Abed-Meraim | An improved assumed strain solid-shell element formulation with physical stabilization for geometric nonlinear applications and elastic-plastic stability analysis[END_REF]

| The twisted beam under out-of-plane loading

Finally, we present the nonlinear calculation of the twisted beam under out-of-plane loading. In this version, the beam is assumed to have an intermediate thickness h = 0.032 and the load goes up to P max = 40 (see Figure 9). This test case is very popular in the context of small strains, but is considered less often in geometrically nonlinear calculations.

However, it does show the capability of an element to undergo large displacements and large rotations. For the calculations using NURBS, we followed the same strategy as for the linear case: we approximated the geometry with a single NURBS element, then kept the same geometry during the refinement. In this case, since there are no spurious modes, the Mixed 2 element could have been used. Nevertheless, to be consistent with the study carried out in the linear case, we preferred considering the stabilized mixed element. The results with the Mixed 2 element would have been very similar. The load step was chosen to be ∆F = 1. We studied the normalized load curves as functions of the displacements of the center of the section at end point (V (A), -U (A),W (A)) (see Figure 9 again for the displacement directions). For the convergence study, 4×1, 8×2 and 16×4 meshes were considered. In the case of element Mixed 2 hg, the convergence was obtained for the 8 × 2 mesh. Thus, the results are shown for that intermediate mesh alone. The solutions are compared to the reference solutions given in Smolenski [START_REF] Smolenski | Statically and kinematically exact nonlinear theory of rods and its numerical verification[END_REF].

F I G U R E 1 9
The deformed configuration (scale factor 1, displacement along y) during the geometrically nonlinear calculation for the twisted beam with Mixed 2 hg and 8 × 2 elements

Using the same approach as before, we first show the deformed configurations for the mixed element (see Fig-

ure 19), then we move on to the (normalized load)-(displacements) responses during the calculation for the higher- F I G U R E 2 0 The load-displacement curves (P /P max = f (V (A)), P /P max = f (-U (A)) and P /P max = f (W (A))) for the twisted beam under out-of-plane loading in the geometrically nonlinear case (8 × 2 mesh) order basic elements and the mixed element (see Figure 20). The Mixed 2 hg element seems to converge more rapidly than the quartic basic element (which is the only basic element which leads to a good approximation of the solution at the end of the calculation with an 8 × 2 mesh). Finally, we compare the Mixed 2 hg element to the SHB8PS 1 element. It was necessary to use a 24 × 4 mesh for the 8-node brick to converge. This mesh and the 8 × 2 mesh for element Mixed 2 hg are compared in Table 5. Again, the behavior of our proposed strategy is better from the point of view of the number of DOFs. However, from the point of view of the number of integration points, the NURBS integration techniques need improvements in order to become attractive compared to the solid-shell strategy using traditional finite elements. This can for instance be done using the reduced Gaussian rule given above in the cylinder's example. Again, no additional instabilities have been reported on this test case with such a strategy. The associated solution Mixed 2 hg reduced is still very close to the Mixed 2 hg solution but, is reached for a number of Gauss points equivalent to the one of the SHB8PS 1 mesh (see again Table 5).

| CONCLUSION

In this paper, we developed a solid-shell NURBS element which is suitable for the geometrically nonlinear calculation of slender structures. The principle of the solid-shell element approach is that the structure is described using a single layer through the thickness of 3D continuous elements. This approach is interesting for the engineer because only displacement DOFs at the control points need to be considered, which makes the implementation of the element relatively easy for complex simulations (large rotations and large displacements, large strains, plasticity, coupling with a solid model, contact, etc). In addition, such an element gives access to the strain through the thickness of the shell directly, which can be interesting for some applications (e.g. forming problems). We proposed an extension of the mixed method presented in a previous contribution (see Bouclier et al. [11]) in order to deal with locking in the case of large rotations and large displacements. The main idea of the mixed method is to modify the interpolation of the average stress components through the thickness. We derived the element's mixed tangent matrix and ended up with a high-performance low-order isogeometric element which is extremely accurate for coarse meshes. Thus, this element provides an interesting alternative to the basic strategy of increasing the order of interpolation, whose drawback is that is leads to very large numbers of Gauss points. Since, in most cases, the NURBS technology enables the exact geometry of the shell to be represented with only a few low-order elements, we were able to show that a NURBS mesh composed of these few low-order elements is sufficient to calculate the shell problem very accurately.

More precisely, this mixed element is of order 2 through the thickness and a complete 3D constitutive relation is assumed. The element also needs to be stabilized in order to prevent the occurrence of spurious zero-energy modes. This was achieved with the introduction of artificial stabilization stiffnesses. Because of the higher continuity of the NURBS functions, we were able to limit the calculation of these artificial stiffnesses to a single element of the mesh while still controlling the hourglass modes throughout the whole structure. Thus, this stabilization is very inexpensive. We limited the approach to the quadratic version of the mixed element. The result is a solid-shell NURBS element of order 2 in the three spatial directions which is at least as accurate as standard NURBS shell elements of order 4 for a given number of DOFs. Test cases for small deformations and, large rotations and large displacements were studied in order to assess the performance of the element. As far as we know, this is one of the first times that a locking-prevention method for shell elements is proposed in the NURBS context for large rotations and large displacements. Here, this technique was developed in the general framework of solid-shell elements in order to take into account as many types of contributions as possible (bending, membrane and shear, and also variation of thickness).

A comparison of our element to an element derived from another NURBS technique with a solid element approach enabled us to confirm the attractiveness of the element developed in this work. Finally, we compared the NURBS solid-shell element with one of the most accurate solid-shell elements based on classical finite element functions of the Lagrange polynomial type. This comparison has led to classical conclusions regarding the comparison of NURBS and Lagrange polynomials for the development of finite elements: the NURBS element is much more accurate than the classical finite element for a given number of DOFs, but it requires a larger number of Gauss points to calculate the integrals. Since the search for more optimal and practical integration rules for NURBS is a very active field of development today, one can hope to see a drastic reduction in the number of integration points in the near future thanks to new rules which will be more suitable for NURBS, including, if necessary, hourglass control following the formalism proposed in this work. For example, it seems at first glance that a simple reduced Gaussian rule would already enable to sufficiently reduce the number of Gauss points to be competitive with respect to techniques proper to classical finite elements.

A | CONSTRUCTION OF THE APPROXIMATE NURBS GEOMETRY OF THE TWISTED BEAM

The approximate geometry of the mid-surface we undertake to develop is defined with a NURBS element of order 2 lengthwise and order 1 across. The control points in (x, y, z) are such that: To illustrate and justify that construction, we show in Figure 22 the method which would be chosen to approximate a helix. In fact, one proceeds in the same way as for building a quarter-circle, except that the control points also move lengthwise.

2.

3.

2 w w w

F I G U R E 2 2
The approximate construction chosen for a helix using NURBS functions Figure 23 shows the superposition of the exact mid-surface (see Formula (63)) and the approximate mid-surface (see the control points defined by (64)). One can verify that the geometries are rather similar. There seems to be an exact match at the edges and in the center, and the error remains small everywhere else. By calculating the distance, in the least squares sense, between the two mid-surfaces, one finds a 3.8% error. Since the geometry is not exact, the problem solved is slightly different from the original problem. Nevertheless, from the calculation point of view, the test case thus modified is just as interesting as the original: the same phenomena (locking, twisted structures) are present. The converged solution obtained is 0.00181 for thickness h = 0.32 (compared to 0.00175 for the finite element solution, i.e. a 3.3% error) and 1357 for thickness h = 0.0032 (compared to 1296 i.e. a 4.7% error). These values were obtained with the Basic 4 element of our NURBS calculation code. Rather than using the solutions of [START_REF] Macneal | A proposed standard set of problems to test finite element accuracy[END_REF], we chose these new values as the reference for this work.

1

 1 The reference problem

  ε x and σ t = P T σ x P = R σ x with the transition matrix P = [e 1 e 2 e 3 ] .

13 T

 13 is the vector of the control variables of σh at control point C . If {U } and { σ } denote the vectors of the components U A , A ∈ {1, 2, .., n } of (

  (a) The mode associated with eigenvalue n o 7 (displacement along z) (b) The mode associated with eigenvalue n o 8 (displacement along y) (c) The mode associated with eigenvalue n o 9 (displacement along x) (d) The mode associated with eigenvalue n o 10 (displacement along x)

  mode in the parent domain Parent domain F I G U R E 4 Hourglass control in the central element: one seeks to describe the modes in the associated parent domain displacement-based formulation whose bilinear form is of the type:

6

 6 Figure 7 shows the deformed shapes obtained with element Mixed 2 (Figure 7(a)) and with element Mixed 2 hg (Figure 7(b)) using an 8 × 8 mesh. There are clearly oscillations in Figure 7(a), whereas in Figure 7(b) the structure

8

 8 (a) Without hourglass control (b) With hourglass controlF I G U R E 7 Initial mesh configuration and deformed shape for the mixed element (scale factor 500, displacements along z) for the pinched cylinder in elasticity under small deformations Convergence of the displacements in elasticity under small deformations for the free-edged pinched cylinder

  ) and the thin beam is shown in Figure10(b). As for the pinched cylinder, no oscillations can be observed in the deformed configurations. Then, we examine the convergence of the displacement of interest as a function of the number of control points lengthwise (See Figure11(a) for the thick beam and Figure 11(b) for the thin beam.) The meshes considered for the calculations were: 1 × 1, 4 × 1, 8 × 2 and 16 × 4.

1

 1 Convergence of the displacements for the twisted beam in elasticity under small deformations

F

  Figure13shows the deformed configuration of the hemispherical shell at the end of the calculation with element

Mesh SLSBEZ 3 ,TA B L E 3 4

 334 Comparison of the numbers of DOFs and Gauss points for the different convergence meshes Load-displacement curves (P /P max = f (u)) for the pinched hemispherical shell with a hole in the geometrically

F I G U R E 1 6 F I G U R E 1 7

 67 (a) The deformed shape for P = 4, 000 (b) The deformed shape for P = 6, 000(c) The deformed shape for P = 10, 000 The deformed configuration (scale factor 1, displacement along z) for the pulled-out cylinder in a geometrically nonlinear calculation with Mixed 2 hg and 8 × 8 elements The P /P max = f (W (A)), P /P max = f (-U (B)) and P /P max = f (-U (C )) load-displacement curves for the pulled-out cylinder (8 × 8 mesh) in the geometrically nonlinear calculation

2 F 4 5

 245 I G U R E 2 1 Comparison of the Mixed 2 hg solutions with the SHB8PS 1 solutions of Abed-Meraim and Combescure [14] for the twisted beam Mesh Mixed 2 hg, 8 × 2 Mixed 2 hg reduced, 8 × 2 SHB8PS 1, 24 × Number Comparisons for the twisted beam: number of DOFs and Gauss points for the various convergence meshes

P 1 =P 4 = 2 .

 142 (0, -b/2, 0) ; P 2 = (L/2, -b/2, b/2) ; P 3 = (L, 0, b/2) ; (0, b/2, 0) ; P 5 = (L/2, b/2, -b/2) ; P 6 = (L, 0, -b/2) ; (64) with the associated weights: ω 1 = ω 3 = ω 4 = ω 6 = 1 and ω 2 = ω 5 = 1/ √ The control points are numbered first lengthwise, then across the width (length n o 1: points 1 to 3 ; length n o 2: points 5 to 6).

  The eigenvalues of the mixed element's condensed global stiffness matrix for the plate problem

	Eigenvalue n o	1 to 6	7	8	9	10	11	12	< • • • <	900
	Eigenvalue	10 -10	-0.0086	-0.0193	-0.0193	-0.0212	538	1.42 × 10 3	< • • • <	7.88 × 10 6
	Rigid body modes	×
	Hourglass modes		×	×	×	×
	TA B L E 1

  The eigenvalues of the mixed element's stabilized condensed global stiffness matrix for the plate problem

	Eigenvalue n o	1 to 6	7	8	9	10	11	12	< • • • <	900
	Eigenvalue	10 -10	37	504	556	572	538	1.42 × 10 3	< • • • <	7.88 × 10 6
	Rigid body modes	×
	TA B L E 2

  for the pinched cylinder Comparisons for the cylinder: number of DOFs and Gauss points for the various convergence meshes

	Mesh	SLSBEZ 3, 16 × 16	Mixed 2 hg, 8 × 8	Mixed 2 hg reduced, 8 × 8	SHB8PS 1, 30 × 20
	Number of DOFs	2527	900	900	3906
	Number of Gauss points	4096	1728	512	1200
	TA B L E 4

Acknowledgements

Alain Combescure's contribution was partially supported by the AREVA-SAFRAN chair on life extension and manufacturing processes. This support is gratefully acknowledged.

Calculation of the displacement increment:

Calculation of the displacement of State 2: U h n 2 = U h n 1 + ˊUh n ; Calculation of the inner loading of State 2 using Equation ( 45):

Calculation of the stress of State 2 using Equation ( 43):

Calculation of the residual: R n = F ext 2 -F int 2 ; Updating of the displacement and of the unknown stress part:

Zeroing of the loading increment: ˊF0 = 0;

Storage of the solution;

Updating of the external loading: F ext 1 = F ext 2 ; end while nonlinear hourglass control, we consider that the shapes of the hourglass modes are the same as those which were observed under small deformations, see Equation [START_REF] Bouclier | Locking free isogeometric formulations of curved thick beams[END_REF]. Now, these are expressed as functions of the second Piola-Kirchhoff stresses, which can be rewritten as follows: