Reduced representations of assumed fields for Hu–Washizu solid-shell element

Mixed eight-node (hexahedron) solid-shell elements based on the standard or partial version of the three-field Hu–Washizu (HW) functionals are developed for Green strain. Three reduced representations of the assumed stress/strain fields are selected. They improve effectiveness, yet retaining good accuracy and convergence properties. At the outset, the standard HW functional and the assumed stress/strain representations of the 3D solid element B8-15P (Weissman in Int J Numer Methods Eng 39:2337–2361, 1996) are used to derive a solid-shell element with 51 parameters. To eliminate locking, the ANS method is applied to the thickness strain (Betsch and Stein in Commun Numer Methods Eng 11:899–909, 1995) and to the transverse shear strain (Dvorkin and Bathe in Eng Comput 1:77–88, 1984). It is a correct element which, however, yields too large displacements for coarse meshes and trapezoidal through-thickness shapes. To improve the above formulation, the ζ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\zeta $$\end{document}-independent reduced representations of the assumed stress/ strain fields are selected and the transformations to Cartesian components are modified. The thickness strain is enhanced by the EAS method. The element with 35 parameters is derived from the standard/enhanced HW functional, but, to further reduce the assumed fields, partial/enhanced HW functionals are constructed from the 3D potential energy by applying the Lagrange multiplier method only to selected strain components. In the element with 27 parameters, this is applied to the constant in-plane strain and to the transverse shear strain while in the element with 19 parameters, to the constant in-plane strain only.Two other modifications are implemented to enhance the behavior of these elements: (A) the skew coordinates are used in the reduced representations of the in-plane stress/strain (Wisniewski and Turska in Int J Numer Methods Eng 90:506–536, 2012), and (B) the Residual Bending Flexibility correction of the transverse shear stiffness (MacNeal in Comput Struct 8(2):175–183, 1978) is adapted. Finally, the performance of the proposed solid-shell HW elements is demonstrated on several linear and non-linear examples for the linear elastic material and the hyper-elastic material. The proposed elements are compared to each other and to the best existing elements of this class.


Introduction
The 8-node hexahedron solid-shell elements have already achieved a considerable level of maturity and have been successfully applied to model shell structures using various 3D constitutive laws, such as the finite strain J 2 -plasticity [44], the Functionally Graded material [45], and others. Nodal degrees of freedom at the bounding surfaces of the solid-shell elements allow for a convenient aggregation of layers and stress analysis of multi-layer composites [43], [36]. Besides the 8-node also 18-node solid-shell elements, based on biquadratic in-plane interpolations, have been developed, e.g. [14].

3D elements
The currently used 8-node solid-shell elements stem from the earlier developed 3D solid elements, such as e.g. the stress hybrid element of Pian and Tong [33], which is based on the Hellinger-Reissner functional and 18 stress parameters. Equivalent to it is the EAS element with 30 strain parameters of Andelfinger et al. [2] and the Hu-Washizu element of Weissman [46]. The latter element uses the above 18parameter assumed stress, and the assumed strain with 33 parameters, where the first 18 have the same form as that used for the assumed stress, while the remaining 15 are associated with the incompressibility condition.
These three elements are of very good accuracy and pass the 3D constant strain patch test but involve a large number of parameters. The subsequent research aimed to reduce this number and to make these elements suitable for finite strain problems; a comprehensive overview can be found in Wriggers [57]. The most commonly used is the Enhanced Assumed Strain (EAS) method of Simo and Rifai [39], which is flexible and can employ the strain-driven constitutive algorithms. Wriggers and Reese [59] identified the problems related to the hourglass instabilities appearing for compression of the enhanced 2D plane strain 4-node elements of a compressible neo-Hookean material. Several methods of controlling the spurious modes were tested afterwards, see the overview in Korelc et al. [23]. A number of very good elements of this class were developed, see e.g. Simo et al. [37], Korelc and Wriggers [24], Glasier and Armero [13], Puso [35] and Korelc et al. [23].
We note that the 8-node 3D solid elements typically pass the 3D and membrane patch tests but fail the bending patch test, which puts in question their use in shell applications with one layer of elements through thickness.

Solid shell elements
The 8-node solid-shell elements, which are the subject of the current paper, use the same interpolation functions and constitutive modules as the 3D solid elements but also additional specialized methods to improve their behavior in thin shell applications involving bending/twisting-dominated problems.
1. The 0th order thickness strain is improperly approximated for curved or trapozoidal through-thickness shape of elements (deformed or undeformed), which causes the so-called curvature thickness locking. The Assumed Natural Strain (ANS) method in the form proposed in Betsch and Stein [6] is used to circumvent it. 2. The out-of-plane bending is impaired by the zero value of the 1st order thickness strain, which is called the volumetric (or dilatational or Poisson's ratio or Poisson's thickness) locking. To remedy this problem, the thickness strain is enhanced either using the EAS method or a specific representation of the assumed thickness strain; this question is discussed in more detail in the sequel. 3. To reduce the transverse shear locking, the Assumed Natural Strain (ANS) method in the form proposed by Dvorkin and Bathe [12] is applied to the 0th order transverse shear strains.
These modifications significantly improve the 8-node solidshell elements' behavior and for this reason are indispensable in this class of elements. Considering the standard shell models (Kirchhoff-Love or Reissner-Mindlin), the normal strain E 33 is equal to zero due to inextensibility of the shell's director. The classical remedy is to recover this strain from the zero-normal-stress condition, S 33 = 0, which for the strains E 11 and E 22 being linear in the thickness coordinate ζ ∈ [−1, 1], yields the E 33 , which is also linear in ζ . Then, however, the 3D constitutive equations must be modified (reduced), which is cumbersome for non-linear materials.
This problem can also be alleviated differently, by enhancing the shell kinematics with parameters describing the extension of a director. In Makowski and Stumpf [30,40], the position vector in the deformed configuration is a quadratic polynomial of the thickness coordinate where d is the current director, λ 0 and λ 1 are the thickness stretches and a 3 is the forward-rotated normal vector. In Simo et al. [38], the extensible-director kinematics with one thickness stretch λ is used, where Ill-conditioning in the thin limit exhibited by this formulation is eliminated by applying the logarithmic thickness stretch μ = ln λ. Equation (2) is used also in the formulation based on the Green strain in Betsch et al. [5], so the thickness strain E 33 = 1 2 (λ 2 − 1). The initial director is defined as D = D/ D , where D is obtained by a linear interpolation of the nodal directors D I . Note that Eq. (2) does not contain the term ζ 2 λ 1 of Eq. (1), the purpose of which is to improve the bending behavior.
For solid-shells, we consider that the thickness strain E 33 is not necessarily normal to the reference surface as is the normal strain. A constant part of the thickness strain E 0 33 is computed from displacements of the top and bottom bounding surfaces, so only the ζ -dependent strain enhancement needs to be added. For 8-node 3D solid elements, a 4-parameter EAS enhancement was used in Andelfinger and Ramm [1] Eq. (34), where ξ, η ∈ [−1, 1] are the in-plane natural coordinates. Moreover, the first 3 terms of Eq. (3) were found to be sufficient to avoid the volumetric locking. The above enhancement was subsequently applied to solid-shells in Büchter et al. [7] Eq. (35) to relieve the Poisson's thickness locking. Later on, Vu-Quoc and Tan [43] established that the first 3 parameters of Eq. (3) suffice to pass the out-of-plane bending patch test for solid-shells. The reference values for this test were worked out in this paper and since then this test became a standard for solid-shell elements. The enhancement of Eq. (3) can also be accounted for in the assumed thickness strain representation of the HW elements, as e.g. in the 3D 8-node solid of Weissman [46] and the solid-shell of Klinkel et al. [20], where the first 3 terms of it are used.
Note that the thickness strain enhancement can also be formulated using the incompatible displacement (1 − ζ 2 ) w 3 d, where w 3 is an additional parameter, see Parisch [31]. This form of the enhancement, however, renders a problem with conditioning of the stiffness matrix in the thin limit.
Finally, we mention the papers on the solid-shell elements discussing the modifications necessary for hyperelastic and elastoplastic material laws. In Hauptmann et al. [16], they include the right Cauchy-Green deformation tensor, and the deformation gradient for elastoplasticity to make them consistent with the modified Green strain. In Hauptmann et al. [15], various methods to avoiding the volumetric locking for a nearly-incompressible material are compared.

Objectives of the current paper
The objective of the current paper is to develop several new eight-node (hexahedron) solid-shell elements based on either the standard or partial version of the three-field Hu-Washizu (HW) functionals. The partial HW functionals are constructed from the 3D potential energy by applying the Lagrange multiplier method to the difference between compatible and assumed strain for selected components only. The purpose is to obtain an effective element of high accuracy and good convergence properties using a minimal number of the assumed stress/strain parameters.
1. We start with developing the solid-shell HW element based on the standard HW functional with 51 parameters (HW51). The assumed stress/strain representations are as in the 3D solid element B8-15P of Weissman [46], but are treated here as contravariant, hence, a different transformation to the Cartesian components is used. This transformation is also different compared to this of Klinkel et al. [20]. Additionally, to eliminate locking for shell-type structures, we apply the ANS method in two forms, one proposed by Betsch and Stein [6] for the thickness strain E 0 33 and the other by Dvorkin and Bathe [12] for the transverse shear strain E 0 α3 . This element is stable, passes the membrane/bending patch tests, and yields quite accurate results in the benchmark tests. Nonetheless, it uses a large number of parameters and is slightly too flexible, i.e. the mesh convergence is from above, so the displacements for coarse meshes are excessive in some tests.
Next we aim at obtaining elements of better effectiveness and accuracy than HW51, the standard or the partial version of the HW functional is used. Both the earlier mentioned ANS methods are also applied to the below-described three elements, which we refer to in this paper as the reduced representation elements.
For the standard HW functional, the formulation used in HW51 is modified in two ways: (a) the assumed stress/strain representations are reduced by omitting the ζ -dependent terms, and (b) the thickness strain E 33 is enhanced by the EAS method using the ζ -dependent representation while the assumed thickness stress/strain representations are omitted. The above changes yield the element with 35 parameters (HW35).
Subsequently, the assumed representations are further reduced and two elements based on the partial HW functionals are developed. The Lagrange multiplier method is applied to the in-plane strain E 0 αβ and the transverse shear strain E 0 α3 in the element with 27 parameters (HW27B), and only to the in-plane strain E 0 αβ in the element with 19 parameters (HW19).
Two other modifications are implemented in the tested solid-shell elements to enhance their behavior: (A) the skew coordinates are used in the representations of the assumed in-plane stress/strain in the reduced representation elements, and (B) the Residual Bending Flexibility (RBF) correction of the transverse shear stiffness of MacNeal [26] is adapted to solid-shell elements. This correction is implemented also in the reference solid-shell elements.
3. The above four solid-shell elements are tested together with the two reference solid-shell elements, which are currently considered as the best in this class. This is the HSEE element of Klinkel et al. [20] and the EAS10 element, which is characterized in Sect. 5. Several other 2D and shell elements are also used in comparisons.
Very good performance of the developed solid-shell HW elements is demonstrated on several linear and non-linear examples for the linear elastic (SVK) material, and such aspects are considered as: accuracy, effects of skew coordinates and the RBF correction, and convergence properties of the Newton method and the arc-length method. Two examples are concerned with nearly-incompressible materials; the dilatation for the SVK material and ν ≈ 0.5 is tested and the hourglassing in the large strain compression for the modified neo-Hookean hyper-elastic material is checked.

Outline of the paper
The outline of the paper is as follows: the general characteristics of the solid-shell elements are provided in Sect. 2. The three-field Hu-Washizu functionals and the assumed stress/strain representations of the developed mixed/enhanced HW solid-shell finite elements are described in detail in Sect. 3. Specialized treatment of thickness and transverse shear strains is described in Sect. 4. Numerical tests in Sect. 5 demonstrate the performance of the developed reduced representation solid-shell elements. The paper is closed by final remarks in Sect. 6.
Notation: "parameter" is abbreviated to "p". In the elements' designations, the total number of elemental parameters is included, e.g. HW35 has 35p. Besides, regarding the components, "COV" stands for "covariant", "CTV" for Numbering of nodes and the reference elemental basis of 8node solid-shell element. The element can be flat or warped. ζ is the through-thickness coordinate "contravariant" and "CART" for "Cartesian". The elemental parameters are denoted as q i , i = 1, . . . , N q .

General characteristics of solid-shell elements
In this section, first the general characteristics of the solidshell 8-node elements is given and, next, their kinematics is described.

Basic definitions for solid-shell element
Consider a 8-node isoparametric solid-shell element with the nodes numbered as shown in Fig. 1. The element can be either flat or warped, similarly as the standard 4-node shell element, see e.g. [48]. The nodal "directors" are defined as the vectors linking the corresponding nodes at the bottom and top surfaces, i.e. 1-5, 2-6, 3-7 and 4-8. If these vectors are mutually skew (non-parallel), the element has a trapezoidal through-thickness shape. Consider the following vectors associated with the solidshell element: the initial position X, the current position x, and the displacement u. The above vectors are interpolated as follows: where the standard tri-linear shape functions are, ξ, η, ζ ∈ [−1, +1] are the natural coordinates and {ξ I , η I , ζ I } are the natural coordinates of nodes I = 1, . . . , 8. We assume that the ζ -coordinate is associated with the thickness direction; the reference surface is at ζ = 0 while the bounding top/bottom surfaces at ζ = ±1. The thickness direction (ξ = const., η = const., ζ ) is in general not normal to the reference surface. For a solid-shell element, the reference Cartesian basis {i k } (k = 1, 2, 3) is constructed as follows. First, the vectors of the natural basis on the reference surface ζ = 0 are defined as where in the denominator k is a selector of the coordinate w.r.t. which we differentiate. In general, g k are neither unit nor orthogonal to each other. The co-vectors g k are defined in the standard way by g k · g l = δ k l (k, l = 1, 2, 3). Let us assume that g 1 , g 2 are the tangent vectors and define in the tangent plane, mutually orthogonal vectors i 1 , i 2 equally distant from g 1 , g 2 , see Fig. 2 in [52], where the auxiliary vectors arẽ The normal unit vector i 3 is perpendicular to the tangent vectors g 1 and g 2 , i.e.
As it is also perpendicular to i 1 and i 2 , the vectors i k form a Cartesian basis {i k }. We construct the elemental reference basis at the element's , and transform X and u to this basis. Then, in the element, the local position vector is X = X i 0 1 + Y i 0 2 + Z i 0 3 , and the Jacobian matrix where i 0 1 = [1, 0, 0] T , i 0 2 = [0, 1, 0] T and i 0 3 = [0, 0, 1] T ; for more details on Jacobians see [48] Sect. 10.3. The Jacobian matrix at the element's center is J 0 Remark For solid-shells, the "thickness" is defined as h = 2 g 3 and it is different from the standard thickness used for shells. The latter can be defined for solid-shells as a projection of the natural vector g 3 on the normal vector i 3 , i.e. h = 2 (g 3 · i 3 ), where g 3 is obtained from Eq. (6) for k = 3.

Kinematics of solid-shells
The configuration space of the Cauchy continuum is defined as: where χ is the deformation function defined on the reference configuration of the body B. The deformation function χ : x = χ (X) maps the reference (non-deformed) configuration onto the current (deformed) one. The deformation gradient is defined as where X is the position vector in the initial (non-deformed) configuration and x in the current (deformed) one. The right Cauchy-Green deformation tensor and the Green strain are defined as where C 0 . = C| x=X . We can linearize the strain E in the thickness coordinate ζ at the reference surface ζ = 0 as follows: where the 0th and 1st order strains are We use the convective coordinates ξ , and parameterize the position vectors as X(ξ ) and x(ξ ). The convective coordinates are associated with the natural coordinates ξ, η, ζ ∈ [−1, +1], which are used in Eq. (4) to parameterize the domain of a finite element, and then ξ . = {ξ, η, ζ }. For the components in the Cartesian reference basis {i k }, we obtain where x ,ξ . = ∂x/∂ξ and J . = ∂X/∂ξ is the Jacobian matrix. Besides, x = X + u where u is the displacement vector. From Eqs. (12) and (15), we see that the CART → COV and COV → CART transformations of components of the Green strain E are where the covariant (COV) components of strain can be obtained as either = ∂x/∂ζ are vectors of the natural basis {ḡ i } in the deformed (current) configuration.
Finally, we note that if the components E COV are known then the 0th and 1st order strains of Eq. (14) are expressed in terms of the CART components as, The last simplified form of E 1 is obtained using two assumptions:

Mixed/enhanced HW solid-shell finite elements
In this section, first the standard and partial Hu-Washizu functionals are described, and next the assumed stress/strain representations of the developed solid-shell HW elements are characterized. Let us distinguish the tangent (in-plane) components α, β = 1, 2, the through-thickness component 33 and the transverse shear components α3. The 0th order and the 1st order parts of stress/strain are designated by the superscripts "0" and "1", respectively, i.e where ζ ∈ [−1, 1] is the through-thickness coordinate and i, j = 1, 2, 3.

Standard/partial HW functionals for solid-shells
Standard 3D HW functional Let us consider the 3D formulation based on the 2nd Piola-Kirchhoff stress S and the Green strain E. The standard Hu-Washizu (HW) functional is mixed by definition and in addition to the compatible displacements u, it involves also two independent fields of stresses and strains, designated as S * and E * , respectively. In some elements, the enhanced strain E(u) + E enh is used, where E enh is the EAS enhancement of Sect. 4.1.
The standard form of the three-field HW functional is as follows: where the strain energy density W(E * ) is expressed by the independent strain E * . The independent stress S * plays the role of the Lagrange multiplier of the constraint linking the independent strain E * and the compatible strain E(u). Also, F ext is the potential of the external loads, the body force, and the displacement boundary conditions and V is the volume of the 3D body in the initial configuration. In the current paper we use the strain energy density W for the St.Venant-Kirchhoff (SVK) material and for the modified neo-Hookean hyper-elastic material (Sect. 5.3.6). We develop the following two solid-shell elements based on the standard HW functional: 1. For the solid-shell element HW51, the standard HW functional is where 51 elemental parameters are used. In this element, the assumed stress/strain representations depend also on the thickness coordinate ζ , see Sect. 3.4.1. 2. For the solid-shell element HW35, the standard/enhanced HW functional is where 35 elemental parameters are used, and E 1 enh i j results from the enhancement of the 1st order thickness strain, see Eq. (60). In this element, all the assumed stress/strain representations are ζ -independent, and are defined in Sect. 3.4.2.
In both the above elements, the ANS method is used for the covariant 0th order thickness strain and the transverse shear strain, see Sects. 4.1 and 4.2 Partial HW functional for solid-shell elements To obtain an efficient finite element, a minimal number of elemental parameters, which guarantee the element's stability and an acceptable accuracy, should be used. Typically, the standard HW functional is reduced to either the Hellinger-Reissner (HR) functional or the Enhanced Assumed Strain (EAS) functional. Another possibility exists for shells and solidshells. We can introduce the partial HW functionals, which are obtained by the application of the Lagrange multiplier method to selected components of strain only.
We develop two solid-shell elements based on the partial shell HW functional, with the Lagrange multiplier method applied to the 0th order strain E 0 i j as specified below: 1. For the solid-shell element HW27B, the partial/enhanced HW functional is where the Lagrange multiplier method is applied to E 0 αβ and E 0 α3 while E 0 33 and E 1 αβ appear in the strain energy in the standard way. The underlined term in W is the sum of the assumed 0th order strain E 0 * αβ and the compatible 1st order strain E 1 αβ multiplied by ζ , and, S 0 * αβ and S 0 * α3 are the Lagrange multipliers. In this element, all the assumed stress/strain representations are ζ -independent, and are defined in Sect. 3.4.3. 2. For the solid-shell element HW19, the partial/enhanced HW functional is where the Lagrange multiplier method is applied only to E 0 αβ while E 0 33 , E 0 α3 and E 1 αβ appear in the strain energy in the standard way. S 0 * αβ is the Lagrange multiplier. Note that the same underlined term appears in Eqs. (23) and (24). In this element, all the assumed stress/strain representations are ζ -independent, and are defined in Sect. 3.4.4. In both the above functionals, E 1 enh αβ , E 1 enh 33 and E 1 enh α3 result from the enhancement of the 1st order thickness strain, see Eq. (60). We note that, compared to the element HW51, the elements HW35, HW27B and HW19 use the reduced assumed stress/strain representations, and in the sequel they are referred to as the reduced representation elements.
Finally, taking a variation of each of the governing functionals of Eqs. (21)- (24) and performing a consistent linearization of it, we obtain the tangent stiffness matrix K for a particular element. This procedure is standard and, therefore, not discussed here.

Skew coordinates
Consider two bases associated with the element's center: the reference Cartesian basis {i k } and the natural basis {g 0 k } (k = 1, 2, 3). The position vector relative to the element center X 0 is given by X − X 0 , and it can be expressed in these two bases as follows: where the skew coordinates ξ S , η S , ζ S are associated with the natural basis {g 0 k }, and x, y, z are the Cartesian coordinates. By taking a scalar product of this equation with the vectors i k , we obtain three equations, which can be solved for the skew coordinates, ⎡ where J 0 is the Jacobian of Eq. (10) at the element's center. As x, y, z depend on the natural coordinates ξ, η, ζ , hence, for a 3D 8-node (tri-linear) element, Eq. (26) becomes where A kl (l = 1, . . . , 4) are the scalar coefficients. For a 2D 4-node (bi-linear) element, Eq. (27) is reduced to where A 11 = (J ,η ) 0 /J 0 , A 21 = (J ,ξ ) 0 /J 0 and J = det J. For parallelogram shapes, A i j = 0, so the difference between the natural and skew coordinates vanishes. Note that the skew coordinates are introduced above quite generally, without resort to any particular approximation problem. In Yuan et al. [60], these coordinates appear as a by-product, when the equilibrium equations for 2D 4-node element are enforced for the 12-parameter representation of the assumed stress in terms of ξ, η.
Remark An important property of the skew coordinates is that the homogeneous equilibrium equations and the conditions of compatibility of strains are satisfied point-wise regardless of the element's shape. For the natural coordinates, they are satisfied point-wise only for parallelogram shape elements, as verified for 2D in [50,51]. In effect, the HR and HW elements with the assumed stress/strain representations expressed by the skew coordinates have an improved accuracy for irregular element's shapes. This is also true for 4-node HW shell elements with 6 dofs/node, see [52].
For the 8-node solid-shell elements, we have to choose which of the two above defined forms of skew coordinates to use. In this paper, we selected the 2D version of Eq. (28) because it is simpler and its properties are well established and tested. Besides, the effects of skewness of g 0 3 and g 0 α are well circumvented by the ANS method for the thickness strain of Betsch and Stein [6], see Sect. 4.1.
We stress that the skew coordinates are used only in the reduced representation elements, to define the in-plane components (11, 22 and 12) of the assumed stress/strain representations, while the natural coordinates are applied in the other components.

Transformation operators for strain/stress vectors
Because of the symmetry of strain and stress tensors, instead of matrices we can use the vectors of their components and define the transformation matrices to obtain these components with respect to another basis. Let us define the transformation matrix where J ik = g i ·i 0 k (i, k = 1, 2, 3) are the components of the Jacobian J and a, b are scalars. To perform the transformations between the contravariant (CTV) components and the Cartesian (CART) components of vectors (29), we define two operators The use of these operators is equivalent to the matrix operations, as specified below.
1. The CTV → CART transformation of components of strain can be written either as We can check their equivalence, i.e.
where (·) v designates the operation of taking out the components of a matrix to obtain the strain vector of Eq. (29). Analogous equivalent relations can be written for components of stress, The modified versions of the above transformations are used in Klinkel et al. [20]; we use them also in our HW solid-shell elements but they are modified in a different way than in the cited paper. 2. The COV → CART transformation of components of strain can be written either as where the inverse of T S is used to transform the strain vectors E v . This transformation is used e.g. in Weissman [46] for the assumed strain together with that of Eq. (33) for the assumed stress. We use it only for comparison in some examples. In the sequel, we additionally indicate where the operators T E and T S are computed. Then the superscript "0" stands for the element's center (ξ = η = ζ = 0), while "ζ = 0" for the reference surface.

Assumed representations of strain/stress
The independent strains and stresses E * and S * are treated as the assumed fields, i.e. E * → E a and S * → S a , which is indicated by the superscript "a". These fields are selected separately for each solid-shell element described below.

Element HW51
This element is obtained from the standard three-field F H W functional of Eq. (21). It is based on the assumed representations which were selected earlier for 8-node 3D elements: the 18p representation for stress of Pian and Tong [33] and the 33p for strain used in the 3D element B8-15P of Weissman [46]. Totally, the HW51 element has N q = 18 + 33 = 51 additional parameters q i .
Assumed stress The CTV → CART transformation of components of stress of Eq. (33) is applied in the following form: where T 0 S is computed at the element's center. The assumed CTV representations are where C 11 = q 7 η + (q 8 + q 9 η)ζ, C 22 = q 10 ξ + (q 11 + q 12 ξ)ζ, Note that only 3 components, i.e. 11, 22 and 12, are ζdependent. Totally, 18 parameters q i are used for the assumed stress. This 18p representation was proposed in Pian and Tong [33] for the assumed stress of 3D 8-node hybrid element. Subsequently it was also used in other 3D elements, e.g. in Andelfinger and Ramm [1] Eq. (34) and Weissman [46] Eqs. (41)(42). Klinkel et al. [20] used it in the solid-shell element HSEE, where the transformation operators for particular terms of this representation are computed at different locations within an element.
Assumed strain The CTV → CART transformation of components of strain of Eq. (32) is applied as follows: where T 0 E is the operator at the element's center. The assumed CTV representation consists of 2 parts, where the first one (underlined) has the same form as specified for stress in Eqs. (36)(37) and involves 18p. The second part contains 15p and is as follows: where Totally, 33 parameters q i are used for the assumed strain. Note that C 1 33 involves the 3p which are used in Eq. (59), therefore, the EAS enhancement of E 33 is not required.

Remark
Note the difference between the transformation rule for the assumed strains used in Weissman [46] Eq. (46) and ours of Eq. (38). In that paper, the COV → CART rule is used (it is identical to Eq. (34) in the current paper) while we use the CTV → CART rule of Eq. (32). These rules are tested and compared for 4-node 2D HW elements in Wisniewski et al. [56]. The CTV → CART rule is found to be in general more accurate, see e.g. the two-element distortion test of Sect. 3.3 therein. It allows for the reduction of the number of parameters used for the assumed strain.

Element HW35
This element is obtained from the standard/enhanced functional F HW35 of Eq. (22). Totally, it has N q = 13+19+3 = 35 additional parameters q i , where 13p are used for the assumed stress, 19p for the assumed strain and 3p to enhance the thickness strain.
Assumed stress The CTV → CART transformation of components of stress of Eq. (33) is applied as follows: where T 0 S is computed at the element's center and T ζ =0 S at the reference surface (ζ = 0). The assumed CTV representations are defined by A v and C v of Eq. (36), where The above 13p representation is obtained from the 18p representation of Eq. (37) by omitting the ζ -dependent terms and replacing the natural coordinates ξ, η by the skew ones ξ S , η S in C 11 and C 22 . Besides, the operator T S in the second term of Eq. (35) is at the element's center while here in Eq. (41) at the reference surface ζ = 0. Note that C 12 = 0 but, because of the presence of q 4 in A v , the matrix which is inverted during the condensation out of additional parameters is non-singular. Assumed strain The CTV → CART transformation of components of strain of Eq. (32) is applied as follows: where T 0 E is computed at the element's center and T ζ =0 E at the reference surface (ζ = 0). The assumed CTV representations are defined by A v and C v of Eq. (36), where Totally, 19 parameters q i are used for the assumed strain. The 0th order assumed thickness strain E a 33 involves 4 parameters; q 3 in A v and 3p in C 33 . Also, 3p are used by the EAS method of Eq. (59) to enhance the 1st order thickness strain.
Comparison of HW35 and HSEE. The main differences between our HW35 and the reference HSSE (43p) of Kinkel et al. [20] are as follows: 1. HSEE uses 43 parameters; 18 for the assumed stress, 18 for the assumed strain, and 7 for the enhancement of the 11, 22 and 33 strain components. HW35 uses 35 parameters; 13 for the assumed stress, 19 for the assumed strain and 3 to enhance the 33 strain component. HW35 uses skew coordinates in the in-plane assumed representations. 2. The enhancing strain components of HSEE are where E COV 11 and E COV 22 are the enhancements of the 0th order strain, while E COV 33 of the 1st order strain. In HW35, only the last one is used. 3. Different transformations and representation are used for the assumed strain. In HSSE, the assumed strain is obtained via the CTV → CART transformation of Eq. (32) applied as follows: where T 0 E is computed at the element's center and T E (ξ, η, ζ ) at Gauss integration points. Comparing to Eq. (43) for HW35, we see that the second term T 0 (46) is not used, and T E in the third term is replaced by T ζ =0 E computed at the reference surface. In the assumed representations, the components 11, 22 and 12 are different, and this difference does not vanish if the natural coordinates are replaced by the skew ones. 4. In HSEE, the CTV → CART transformation of components of stress of Eq. (33) is applied in the form analogous to that used for strain, where T 0 S is computed at the element's center and T S (ξ, η, ζ ) at Gauss integration points. Hence, the differences between the transformations used in these two elements are analogous to these for the assumed strain.

Element HW27B
This element is obtained from the partial/enhanced functional F HW27B of Eq. (23), in which only the membrane and transverse shear strain components are treated by the Lagrange multiplier method. Totally, it has N q = 9 + 15 + 3 = 27 additional parameters q i , where 9p are used for the assumed stress, 15p for the assumed strain and 3p to enhance the thickness strain.
Assumed stress The CTV → CART transformation of Eq. (41) is applied, and the assumed CTV representation is defined as follows: the components of C v of Eq. (36) are The skew coordinates ξ S , η S are used in C 11 and C 22 . Totally, 9 parameters q i are used for the assumed stress. Five of them correspond to the Pian and Sumihara membrane stress of [32], but expressed in skew coordinates, and 4 are added for the transverse shear stress.
Totally, 15 parameters q i are used for the assumed strain. Additionally the 3p EAS enhancement of Eq. (59) is used for the thickness strain.

Element HW19
This element is obtained from the partial functional F HW19 of Eq. (24), in which only the in-plane components are treated by the Lagrange multiplier method. Totally, it has N q = 5 + 11 + 3 = 19 additional parameters q i , where 5p are used for the assumed stress, 11p for the assumed strain and 3p to enhance the thickness strain.
Assumed stress The CTV → CART transformation of Eq. (41) is applied, and the assumed CTV representations are defined as follows: and the components of C v of Eq. (36) are where the skew coordinates ξ S , η S are used in C 11 and C 22 . Totally, 5 parameters q i are used for the assumed stress, which corresponds to the Pian and Sumihara stress of [32], but expressed in skew coordinates. Assumed strain The CTV → CART transformation of Eq. (43) is applied, and the assumed CTV representations are defined as follows: and the components of C v of Eq. (36) Totally, 11 parameters q i are used for the assumed strain. Additionally the 3p EAS enhancement of Eq. (59) is used for the thickness strain.
Remark If the bi-linear (underlined) terms in the assumed strain is omitted then we obtain the 9p representation, which is identical to that selected for the 2D element HW14-S [51]. The inverse constitutive equation and the Pian-Sumihara's 5p stress representation was instrumental in selecting this representation. In effect, it is accurate and insensitive to transformations, see [56]. In the 8-node solid-shell element, however, it yields 2 large eigenvalues for the nearly-incompressible material, see Sect. 5.1. The underlined terms render that only one large eigenvalue is obtained for the HW19 element and the volumetric locking is avoided.

Treatment of thickness and transverse shear strains
In this section we discuss the specialized methods used to improve the thickness strain and the transverse shear strain. The covariant strain components are modified.

ANS method and EAS enhancement for the thickness strain
Two problems related to the thickness strain E 33 of solidshell elements can be identified.
A. For trapezoidal through-thickness shape of solid-shell elements, the initial natural vector g 3 .
= ∂X/∂ζ is not normal to the reference surface ζ = 0. Also the current natural vectorḡ 3 .
= ∂x/∂ζ becomes not normal to the reference surface ζ = 0 for certain types of deformation, e.g. the out-of-plane bending. In effect, the covariant 0th order thickness strain E COV 33 is improperly approximated, which causes the so-called curvature thickness locking. We denote the covariant thickness strain at the reference surface as ε 33 .
. d and h are used for the initial configuration,. Then the thickness strain at the reference line, obtained from Eq. (17) for i, j = 3, is For the parameterization of x of Eq. (4), we obtaind(ξ ) = where A and B are the end-points (ξ = ±1) on the reference line. Whend A andd B are nonparallel, such an interpolation ofd(ξ ) does not preserve its length, i.e.
The second term causes the thickness straining in bending, and to circumvent it Betsch and Stein [6] propose the ANS method. Within the ANS method, the product is interpolated using appropriate products at the end-points, It is always equal to 1 becaused A andd B are unit vectors. For the initial "director" d and the product d · d of Eq. (54), we proceed in a similar way. Effectively, we can interpolate ε 33 using the values at the end-points A and B, where (ε 33 ) A and (ε 33 ) B are directly computed using not Eq. (54) but Eq. (17).
For the solid-shells we proceed similarly. First, the values of ε 33 are computed at 4 corner points (ξ = ±1, η = ±1) at the reference surface ζ = 0; they are denoted as (ε 33 ) I , I = 1, 2, 3, 4. Next, this strain is interpolated within an element with the bi-linear shape functions N I (ξ, η), Finally, the Cartesian components of this strain are obtained using the COV → CART transformation of Eq. (34).

B.
In the solid-shell kinematics, the 1st order thickness strain E 1 33 is equal to zero, which impairs element's accuracy in the out-of-plane bending. An overview of the methods which have been used to enhance this strain in shells, 3D elements and solid-shell elements, is given in Sect. 1.
In the current paper, we use the 3-parameter EAS enhancement of the covariant 1st order thickness strain and the standard transformation to Cartesian components, where T 0 S is computed at the element's center and J . = det J. This enhancement is used in all our solid-shells elements except HW51, in which the terms of Eq. (59) are included in the assumed strain representation, see Eq. (40).

ANS method for transverse shear strains
Let us denote the covariant components of the transverse shear strains at the reference surface ζ = 0 as To reduce the transverse shear locking, the Assumed Natural Strain (ANS) method in the form proposed by Dvorkin and Bathe [12] is applied to γ 13 and γ 23 as follows: 1. First, the γ α3 (α = 1, 2) are computed (sampled) at the middle points of element sides, at two points for each component, see Fig. 2. These values are denoted as γ 5 13 The interpolated components (designated by a tilde) are constant in the direction in which the derivative is calculated and linear in the other direction, which means that the ANS method is orientation-dependent. 3. Finally, at the Gauss integration points, the Cartesian transverse shear strains are obtained by the COV → CART transformation of Eq. (34).
Finally, the ANS modified strains of Eqs. (58) and (62) are transformed to the Cartesian components simultaneously with the compatible strain components,

RBF correction of transverse shear stiffness
Low-order finite elements, such as the two-node Timoshenko beam, the four-node Reissner-Mindlin shell and the eightnode solid-shell, lock for the sinusoidal bending because such a form of deformation is not properly represented by linear shape functions. To resolve this problem, the Residual Bending Flexibility (RBF) correction was proposed for beams and shells, see [26,28]. This correction is beneficial also for very thin elements. Below the RBF correction is adapted to the 8-node solid-shell elements.
For a 2D beam element with a rectangular cross-section, we define the ratio of the shear energy W γ to the total energy as follows: where W κ is the bending energy. For the Discrete Kirchhoff (DK) beam element based on cubic displacements (see e.g. [48] Sect. 13.1), the deflection for the sinusoidal bending is where θ is the nodal rotation and l is the element's length. Let us define the RBF correction coefficient c RBF . = c and the corrected shear modulus where the underlined term is called the Residual Bending Flexibility. Note that G is multiplied by the shear correction factor k to account for a parabolic distribution of the transverse shear stress through thickness. Its use is not indicated in [26] but is obvious. Note that the c RBF and G * depend on the material and geometrical properties of the element but not on the rotation θ .
Although the c RBF and G * are derived for the DK beam element, they can also be applied to the Timoshenko beam element based on linear interpolations. Similarly, a generalization of G * can be applied to bi-linear 4-node plate and shell elements.
For the 8-node solid-shell element, we need to account for the transverse shear strains at the reference surface ζ = 0, hence we proceed similarly as for a 4-node shell element. The corrected G * of Eq. (66) is applied separately for each direction, i.e.
where l 1 and l 2 are the lengths of vectors connecting opposite mid-side points at the reference surface. (The shear correction factor k = 5/6 is used in numerical examples of Sect. 5 unless defined differently.) Let us assume that h, G * 1 and G * 2 are constant over the reference surface and express the transverse shear strain energy for a single element as follows: To pass the bending patch test, the ANS method of Sect. 4.2 should be applied to the transverse shear strain. Then the element performs well for bending but yields excessive results for twist. For instance, in the linear test of a slender straight cantilever modeled by one row of 4-node shell elements and twisted by a pair of forces, the rotation r x of the tip is too large by about 38%, see [48] Sect. 15.3.1.
To avoid the excessive twist, the full RBF correction is applied to the average values of the transverse shear strains ("av") and a fraction of it to the difference of them (designated by "d"), see [27]. To separate these parts, we re-write Eqs. (62) of the ANS method as follows: Taking for instance only the two last rows of Eq. (63), we have where each component is a linear function of ξ and η, e.g. γ CART 13 =γ av 13 +γ d1 13 ξ +γ d2 13 η. To separate these terms in the strain energy W γ of Eq. (68), we assume that J (ξ, η) ≈ J 0 , where J 0 is the Jacobian determinant at the element's center. Then Finally, the integrand of the strain energy (68) is modified and the full RBF correction is applied only to the "av" strains and a fraction of it to the "d" strains, with the corrected shear modulus of Eq. (67) additionally modified as follows: where is a corrective coefficient; = 0.04 is selected in [26].Ṅote that too small values of cause problems with the conditioning of the tangent stiffness matrix. We proceed similarly with G * Comparison The shear correction factor for a Timoshenko beam element is defined in Tessler and Hughes [42] Eq. (4.14) as follows: where k 2 * = π 2 /12 is the analytic shear correction factor, which corresponds to our k = 5/6. Besides, r . = √ I /A, and for a rectangular cross-section, when A = bh and I = bh 3 /12, it yields r 2 = h 2 /12. Hence, the inverse of Eq. (74) divided by G is We see that k 2 TH in the above formula is identical to c RBF of Eq. (66) when the analytic shear correction factors (k 2 * and k) are defined in the same way.

Numerical tests
In this section, we describe the numerical tests of four 8node Hu-Washizu (HW) solid-shell elements developed in the current paper, see Table 1. These elements are based on the Green strain and developed from either the standard HW functional of Eqs. (21) and (22), or the partial HW functionals of Eqs. (23) and (24).
The assumed representations of stresses/strains are described in Sects. 3 HW27B and HW19). In these elements, the in-plane components of the assumed strain/stress are expressed by the skew coordinates of Sect. 3.2. The number of parameters used for particular components of stress/strain is given in Table 2. The 0th and 1st order parts of stress/strain are designated by the superscripts 0 and 1, respectively. In the sequel, "parameter" is abbreviated to "p".
Several methods of improving the behavior of solid-shells are applied and they are described as follows: the ANS methods for the thickness strain and the transverse shear strains in Sects The multipliers q i of additional modes are eliminated on the element's level and updated by the scheme U2, see [48]. All solid-shell elements are integrated using the 2 × 2 × 2 Gauss integration rule.
Two reference solid-shell elements are used: 1. HSEE with 43p of Klinkel et al. [20], which is an HW element; it is described in detail and compared to our HW35 element in Sect. 3.4.2. 2. EAS10, in which the 7p enhancement of Wilson et al. [47] is applied to the in-plane strain, and the 3p enhancement of Eq. (59) to the thickness strain E 33 . The first version of this element is designated Q1A3E5 in [19]. In [20], two terms are added to the thickness strain enhancement to pass the bending patch test, as proposed in [43], and it is designated Q1A3E7.
These elements are derived using the automatic differentiation program AceGen described in [22,25], and are tested within the finite element program FEAP developed by R.L. Taylor [41,61]. The use of these programs is gratefully acknowledged. Our parallel multithreaded (OMP) version of FEAP is described in [18].
We tacitly assume that any consistent set of units is used for the data defined in numerical examples.

Eigenvalues of single element
The eigenvalues of the tangent matrix are computed for a single unsupported element, and for the material data: the Young's modulus E = 1 and the Poisson's ratio ν = 0.3.
First, a single element of size 1 × 1 × h and the thickness h = 0.1 is checked. Then, several other element shapes are Ref.
examined, such as a truncated pyramid, a rhomboid with parallel corner "directors" and a warped element with elevated nodes 3 and 7 in Fig 1. For all the shapes tested, our solidshell elements have a correct number of zero eigenvalues (6). For the triangular prism, e.g when nodes 3 and 4, and nodes 7 and 8, have the same coordinates, the number of zero eigenvalues remains correct. In the so-called thin limit, an additional zero eigenvalue appears for all tested solid-shell elements.
For the nearly-incompressible material, i.e. ν = 0.499999999, one large eigenvalue 8.5 × 10 8 and 6 zero eigenvalues is obtained for all tested solid-shell elements. For the distorted shapes and the triangular prism, the number of large and zero eigenvalues remains the same. The absence of volumetric locking was additionally confirmed by the "Pure bending of a cantilever in plane strain" test of MacNeal [28].

Patch tests
The solid-shell elements fail the 3D patch test with nonparallel top and bottom bounding sides. Hence, the shell-type patch tests are adapted and performed.
The five-element patch of elements described in MacNeal and Harder [29] is used, with the material data: the Young's modulus E = 10 6 and the Poisson's ratio ν = 0.25. The displacements are prescribed at the external nodes of the patch and computed at the internal nodes. The reference results for the membrane patch test are the same as for the 4-node shell elements but additionally the values of u 3 and the thickness strain/stress (E 33 , S 33 ) must be checked. For the out-of-plane bending/twisting patch test, the reference results are given in VuQuoc and Tan [43] Table 2.
All our solid-shell elements yield correct displacements at the internal nodes and correct compatible strains and stresses at the Gauss Points.

Dilatation for out-of-plane bending
This test is used to detect the dilatational locking for the nearly-incompressible linear elastic material., The dilatation is the volume change for small strains only. Pure out-of-plane bending for plane strain is defined by the following displacements: for which the small strains are 11 = z, 22 = 0, For the Poisson's ratio ν = 0.5, i.e. for the incompressible material, we obtain 33 = −z so the dilatation e = 11 + 22 + 33 = 0. This test is an adaptation of the pure in-plane bending test of [28], p. 216.
To obtain the nearly-incompressible material, we use the Poisson's ratio ν = 0.499999999, for which the dilatation should be e ≈ 0. Besides, the Young's modulus E = 1, the thickness h = 0.2, the half-length L = 5 and the width w = 1. The RBF correction of Sect. 4.3 is not used.
The 10 × 1 × 1-element mesh is used, with 1 layer of elements in the through-thickness 0Z direction, see Fig. 3. The values of u x and u z of Eq. (77) are applied to 8 corner nodes and at x = z = 0 to nodes 6, 17, 28 and 39. To render the plane strain state, the boundary condition u y = 0 is set for all nodes at y = ±w/2. The obtained results are presented in Table 3. Strains and the dilatation are computed in element no. 5 at point "C", i.e. at (ξ, η, ζ ) = (0, 0, −0.5). The analytical solution at point "C" is 11 = 0.05. For all elements, 22 ≈ 0 and 33 ≈ − 11 are obtained as required, hence these components are not reported. The values of 11 only negligibly differ from the reference analytical value.
The transverse shear strain 13 and the dilatation e are close to zero, which indicates that all the tested and reference solid-shell elements are free from the transverse shear locking and the dilatational locking.

Cook's membrane
In Cook's test [8], the in-plane shear deformation dominates and the elements are skew and tapered. The membrane is clamped at one end, while at the other end, the uniformly distributed tangent load P = 1 is applied, see Fig. 4. The nodes shown in this figure are doubled to obtain the mesh for 8-node solid-shell elements. The data is as follows: E = 1, ν = 1/3, and the thickness h = 1. Two meshes are used in computations; a coarse 2×2-element mesh and a fine 32×32element mesh. For the reference 9-node 2D element the same number of nodes is used, and the 1 × 1 and 16 × 16-element meshes.
The vertical displacements u y at point A obtained in the linear analysis are presented in Table 4. The relative errors [in %] are also provided. All the tested solid-shell elements perform identically for the fine 32×32-element mesh. Regarding the coarse 2 × 2-element mesh, our elements based on the reduced representations and the skew coordinates (HW35, HW27 and HW19) are the most accurate ones but the difference between them and the other elements is small.

Two-element distortion test
This test is used to verify either the in-plane part or the through-thickness part of solid-shell elements. Bending always takes place in the X0Y plane, but the nodal "directors" are either aligned with the 0Z-axis (in-plane bending) or parallel to the X0Y-plane (out-of-plane bending).
The cantilever is modeled by two solid-shell elements, and a tilt of their common side is defined by the parameter d, see Fig. 5. The data is as follows: E = 1500, ν = 0, h = 1, and P = 10. The nodes shown in Fig. 5 are doubled to obtain the mesh for 8-node solid-shell elements, and the pair of forces ±P is replaced by four forces ±P/2. Boundary conditions are applied as described in [20].
In-plane bending In Fig. 6, the u y displacement at the tip is shown for a varying distortion d. Up to d = 1, for which the tilt of the common nodal "directors" is 45 • , all the curves almost coincide. For d > 1, the curves for our HW51 and the reference EAS10 are identical, and the curves for all our reduced representations elements (HW35, HW27B and HW19) coincide.
Note that for larger distortions less accurate solutions are expected, see, e.g., [1] Fig. 8 and [56] Fig. 5. Besides, the curve for a very accurate 2D 4-node element HW14-S [51] has a plateau for d > 1. QE2 of [34] and HR5-S of [50] perform as HW14-S. As the curve for our reduced representations elements is the closest to this 2D solution, it seems to be better than the other curves.
In Fig. 6, we also show the results for two sets of transformation for stress-strain: (a) CTV-CTV of Eqs. (33) and (32) and (b) CTV-COV of Eqs. (33) and (34). For HW51, the curves for both transformations coincide, but for the reduced representations elements there is a significant difference between them. The CTV-CTV transformations yields more accurate solutions, hence, they are used in our solidshell elements.
Out-of-plane bending In Fig. 7, the u y displacement at the tip is shown for a varying distortion d. The results obtained without the RBF correction are presented. Our reduced representations elements (HW35, HW27B and HW19) and the reference EAS10, yield the analytical solution u y = 1 in the whole range of d.
Comparing our HW51 and the reference HSEE, the curves for them coincide and the error of them grows fast for increasing d. For d = 1, i.e. for the 45 • tilt of the common nodal "directors", the relative error is 8.2%. Since the thickness strain is treated in these elements by the ANS method of Sect. 4.1, the cause of the error is not this strain. Note that this error is eliminated in HW35, where the ζ -dependent terms are omitted in the assumed representations and the thickness strain is enhanced by the EAS method.
Finally, all our elements are tested using the two sets of transformations for the assumed stress-strain, CTV-CTV and CTV-COV, and the curves obtained are identical.

Straight cantilever of trapezoidal elements
This classical test of MacNeal and Harder [29] is applied to the transverse deformation of solid-shell elements, as e.g. in Harnau et al. [14]. The accuracy of displacements for a trapezoidal through-thickness shape of solid-shell elements and the effects of the RBF correction of Sect. 4.3 are assessed. The trapezoidal through-thickness shape of solid-shell elements, see Fig. 8, is obtained by tilting the common nodal "directors" by ± 45 • . The vertical force P = 1 is applied at the right end, equally to the top and bottom nodes, causing out-of-plane bending and transverse shearing. The mesh consists of 6 elements, and the aspect ratio for the rectangular elements is 5. The nodes shown in Fig. 8 are doubled in the 0Y direction to obtain the mesh for 8-node solid-shell elements. The boundary conditions imposed at the left end of the cantilever eliminate the rigid body motion but enable the deformation in the 0Y and 0Z directions. The data is as follows: E = 10 7 , ν = 0.3, the length L = 6, the thickness h = 0.2 and the width in the 0Y direction w = 0.1.
The vertical displacement u z at node A is presented in The effects of the RBF correction are generally positive, except for our HW51 and the reference HSEE, for which the accuracy worsens. Nonetheless, for these elements the displacements are larger than the reference solution, excluding locking.
Note that all the reference 2D plane stress elements show locking for the trapezoidal shape but to different degrees. For the 9-node elements it is over 3 times smaller than for the 4node ones, for which the error exceeds − 77%. For both, the errors are much bigger than for the solid-shell elements.
In Table 5, we also evaluate the effects of the ANS method for the thickness strain of Betsch and Stein [6]. Comparing the results obtained with and without this method for the trapezoidal mesh (the latter are designated "EAS10 w/o ANS [6]"), we conclude that it reduces the relative error from − 71.5 to − 0.84% for computations without the RBF, and from − 59.68 to − 0.15% for computations with the RBF correction; a significant reduction indeed.

Remark
The mesh of parallelogram (through-thickness) shape elements is obtained by tilting the common nodal "directors" by 45 • , see Fig. 8. For this mesh and computations with the RBF correction, the maximum relative error of u z is 0.78% for HSEE and 0.5% for HW51. For the other elements, it is much smaller. Generally, the error of parallelogram shape solid-shell elements is much smaller than for the trapezoidal shapes.

Curved 3D cantilever
The element's distorted shape (Case 1, Fig. 9) and the skew (non-parallel) nodal "directors" (Case 2, Fig. 11) affect the accuracy of a solution, especially when the thickness h diminishes. By the nodal "directors" we designate the through-thickness vectors linking the corresponding bottom and top nodes; for a single element they are shown in Fig. 1.
The curved 3D cantilever is fixed at one end and loaded by a moment M z at the other, see Fig. 9. The nodes shown in this figure are doubled in the radial direction to obtain the mesh for 8-node solid-shell elements. The external moment is assumed as M z = (R/h) −3 , so a solution of a linear problem should remain constant. This moment is applied to solid-shell elements as two pairs of opposite tangent forces The data is as follows: E = 2 × 10 5 , ν = 0, width b = 0.025 and radius of curvature R = 0.1. The reference analytical solution for a curved beam is where I is the moment of inertia. The FE mesh consists of 6 solid-shell elements.    Fig. 9; the distortion is defined as in Koschnick et al. [21] p. 2454. For the distorted mesh, the 8-node solid-shell elements are warped (non-planar), which makes this case very demanding and for which the RBF correction can be beneficial.
The scaled displacement u y /u ana y at point A obtained in the linear analysis are shown in Fig. 10, where two sets of results are presented, without and with the RBF correction. (a) Without the RBF correction, the best accuracy provides the reference EAS10, next are all the remaining elements, which perform nearly identically. (b) With the RBF correction, the best accuracy provides the reference EAS10 and ours HW51. Next are the remaining tested and reference elements, which nearly coincide. The curves for our HW35, HW27B and HW19 almost coincide.
For the slenderness log 10 (R/h) = 3, the RBF correction significantly improves the accuracy of displacements; reducing the error from 10 to 2% for EAS10, from 30 to 5% for our HW51 and from 30 to 10% for the remaining elements. For log 10 (R/h) = 4, the accuracy is unacceptable for all elements, without/with the RBF correction.
Case 2: Radial and skew nodal "directors" For the curved cantilever, the nodal "directors" can be either radial or skew (non-parallel) as shown in Fig. 11 for the thickness h = 10 −2 (R/h = 10). For the skew ones, the angle ϕ ≈ 35 • ; it would decrease if h were decreased further. In the circumferential direction, the regular mesh of Fig. 9a is used.
The displacements u y at point A (node 1) obtained in the linear analysis are shown in Table 6. The relative errors [in %] also are presented. The most accurate is the reference EAS10, then our reduced representation HW elements (HW19, HW27B and HW35). The latter elements perform similarly and the RBF correction has a small negative effect for them. For EAS10, this effect is more pronounced. The elements with a large number of parameters, such as our HW51 and the reference HSEE, perform worse than the other elements and effect of the RBF correction is negligible.

Pinched cylinder with diaphragms
The purpose of this example is to show mesh convergence for the tested and reference HW solid-shell elements. This is one of the most severe tests for both inextensional bending and complex membrane states, see Belytschko et al. [4] p. 239.
A cylindrical shell is closed at both ends by rigid diaphragms and is pinched by two opposite forces P applied at the middle section, see Fig. 12a. The material and geometry data of [4] are as follows: E = 3 × 10 6 , ν = 0.3, R = 300, L = 300, h = 3 and P = 1. Because of symmetries, only one-eighth of the cylinder is analyzed, see The vertical displacement w F E M under load at node A on the external surface is monitored, and its value obtained for the element HW35 and N = 128 is used as w R E F = 1.8757 × 10 −5 to obtain the scaled displacements. The results for various N are presented in Fig. 13; the reference curves (Ref. shell) are obtained using the 4-node shell element HW47 with 6 dofs/node of [52] (thick broken lines).
We see that the solutions for the tested and reference solid-shell elements almost coincide and are very close to the reference shell solutions. The positive effects of the RBF correction for coarse meshes are clear; the curves obtained without the RBF correction can be compared to Fig. 13 of [17]. For the densest mesh, the difference between the solutions without and with the RBF correction becomes very small, as expected.
The solutions obtained without the RBF correction converge monotonically from below and can be used to find the convergence rate for varying the number of elements per side N . Because the differences between the curves for our solid-shell HW elements are very small in Fig. 13; we further use only HW35. The relative error is defined as e w = |(w HW35 − w REF )|/w REF and the log 10 e w is plotted against log 10 (L/N ) in Fig. 14, where L/N is the element's size. Next the line is fitted to these 5 points, and y(x) = −3.04157 + 1.51914 x is obtained. Therefore the a)

Roll-up of a clamped beam
This is a standard test to verify finite rotation procedures for beams and shells, but here it is used to asses the quality of solid-shell elements in pure bending. The planar straight beam is clamped at one end and loaded by a bending moment M at the other end, see Fig. 15a. The nodes shown in this figure are doubled to obtain the mesh of 25 8-node solid-shell elements. The moment M is replaced by two pairs of opposite forces ±P, which are tangent to the free edge in the current configuration. The data is as follows: E = 12 × 10 6 , ν = 0, the thickness h = 0.1, the width w = 1, the length L = 10. The RBF correction is irrelevant to this test.
A full circle is obtained for the bending moment M * = 2π E I /L = 628.319. The beam's tip displacements up to  More precise are the values of the tip's displacements at M * given in Table 7, where also the relative error of u x [in %] is presented. All the solid-shell elements slightly lock, but the range of relative error is narrow [1.14-1.28]%. Least locks the reference EAS10 without the ANS for the thickness strain of [6], as then the error is 0.41%.

Pinched hemispherical shell with hole
The hemispherical shell with an 18 • hole is loaded by two pairs of equal but opposite external forces, see Fig. 16. The shell undergoes an almost in-extensional deformation, so a membrane locking of a solid-shell element can be detected by this test.
Due to the double symmetry, a quarter of the shell is modeled. One element is used through-thickness and three meshes with 8×8, 16×16 and 64×64 elements over the sur-face. The nodes shown in Fig. 16 are doubled in the thickness direction to obtain the mesh for 8-node solid-shell elements, which have the shape of truncated pyramids. For reference, the 4-node and 9-node shell elements are used; for the latter the 4 × 4, 8 × 8 and 32 × 32-element meshes are used.
The data is as follows: E = 6.825 × 10 7 , ν = 0.3, thickness h = 0.04 of [29]. Non-linear solutions are also provided for a smaller thickness h = 0.01, for which the elements are more prone to locking.
The results of linear analyses for the thickness h = 0.04 are given in Table 8, where the inward displacement u y under the force P = 1 at the inner node is reported.
In general, the RBF correction renders that the displacements are larger than without using it. The relative error is reported for the 64 × 64-element mesh and the RBF correction reduces its range from [− 0.47, − 0.48]% to [− 0.27, − 0.28]%, so this correction is beneficial. Besides, it causes that all the solid-shell elements converge from above for the increasing mesh density. The non-linear analyses are performed using the Newton method and the solution curves for the 16×16-element mesh   The solution curves coincide for all solid-shell elements in the whole range of load. They exactly follow the reference curve up to a certain level from which they start departing from it. Then the displacements become slightly either excessive or locked. It depends on whether the RBF correction is used or not.
We also include non-linear solution curves for the smaller thickness h = 0.01, see Fig. 17. Comparing Figs. 17 and 18, we see that the RBF correction is certainly more beneficial for the smaller thickness h = 0.01.
One-step non-linear test An overall effectiveness of a finite element in non-linear analyses depends not only on the time of a single Newton iteration but also on the radius of convergence and the rate of convergence, which can be   characterized by the maximum P for which the Newton method converges and the number of iterations required. This test is performed for the pinched hemisphere of h = 0.01 and the 16 × 16-element mesh. The maximum load for which the Newton method converges is found for each solidshell element separately, increasing P by 0.1. The results obtained without and with the RBF correction are presented in Table 9.
We see that the HW elements with a large number of parameters, such as the reference HSEE (43 p) and ours HW51 and HW35, converge for bigger max P and require less iterations. Those with less parameters, such as our HW27B and HW19, converge for a smaller max P. For them the RBF correction is particularly beneficial. The reference EAS10 requires the smallest max P and the RBF correction does not improve this.
Comparing the performance of our reduced representation HW elements, we draw the following conclusions: 1. Our HW51 and HW35 perform similarly, but HW35 uses less parameters. HW27B and HW19 perform likewise, hence the assumed transverse shear stress and strain (S 0a α3 and E 0a α3 ), which are used in HW27B but not in HW19, are not required. 2. Another formulation with 27 parameters has been tested, in which the assumed thickness stress and strain (S 0a 33 and E 0a 33 ) are retained and the assumed transverse shear stress and strain (S 0a α3 and E 0a α3 ) are omitted. Such an element performed worse than HW19, especially without the RBF correction.

Twisted beam
This test belongs is from MacNeal and Harder [29]. The initial geometry of the beam is twisted so all the elements are warped (non-flat), but the initial strain is equal to zero. The beam is clamped at one end and loaded by a force P y at the other, see Fig. 19. The data is as follows: E = 2.9× 10 7 , ν = 0.22, the length L = 12, the width w = 1.1 and the twist angle is 90 • .
In our computations, the 4 × 24-element mesh of the 8node solid-shell elements and a small thickness h = 0.0032 are used. One element is used through thickness. For the reference shell elements, the 4 × 24-element mesh is used for the 4-node HW47 and the 2 × 12-element mesh for the 9node MITC9i. This example is also computed without/with the RBF correction.
The results of the linear analysis for P y = 1 × 10 −6 are shown in Table 10, where the u y × 10 3 displacement at point A and its relative error are given. Without the RBF correction, our HW35 and HW27B are the most accurate, next are the reference elements EAS10 and HSEE, and finally HW51. The relative errors are in the range [0.07, − 0.46]%. With the RBF correction, the most accurate is the reference EAS10, next ours HW27B and HW19, then our HW51 and HW35 and the reference HSEE. The relative errors are in the range [0.10, − 0.67]%.
The non-linear load-deflection curves obtained for P y = 1 × 10 −4 by the arc-length method are shown in Fig. 20, and the displacement u y at point A is monitored. Without the RBF correction, the curves for all tested and reference elements coincide. With the RBF correction, all the curves except this for EAS10 coincide, and are slightly 'softer' than those obtained without the RBF correction. The softest is the reference EAS10, for which the solution curve coincides with that for the reference 4-node HW47.
Summarizing, the beam is extremely thin but the results are not corrupted by the membrane locking and are close to the reference ones.

Short C-beam
This example includes the right angle intersections of elements and is typically used to test shell elements with 3 rotational dofs/node, see [9]. It is useful also for solid-shell elements because of the through-thickness element's shapes, for which the nodal "directors" are skew (non-parallel), see Fig. 21b.
The beam is fully clamped at one end and loaded by two vertical forces P/2 at the other end, see Fig. 21a, b. The data is as follows: E = 10 7 , ν = 0.333, the length L = 36, the thickness h = 0.05. The web is modeled by 36 × 6 elements and each flange by 36 × 2 elements, so the total number of elements is 360. The vertical displacement u z is a) Fig. 19 Twisted beam. a Initial mesh and load. b Deformed mesh at P y = 0.1 The linear reference solution u z = 1.1544 × 10 3 is computed using the 16-node shell element in [10]. This value is confirmed by our results of Table 11, where the error is: (a) 0.1% or less for 4-node and 9-node shell elements, and (b) 0.24% for the 3D 8-node element and a dense mesh of 23040 elements, i.e. for the mesh density multiplier equal to 4, also through-thickness.
The linear solutions for P = 1 and the mesh of 360 elements are given in Table 11, where the displacement −u z × 10 3 and its relative error are shown. Without the RBF correction, the range of relative errors is [  The non-linear solutions for the 360-element mesh are computed using the arc-length method with the initial P = 5, see Fig. 22. The reference curve is obtained using the 4-node shell element HW47 of [52] (thick continuous line). The solution curves for the solid-shell elements with the RBF correction almost coincide and are slightly more stiff than the reference curve. The solution curve for the 3D 8-node solid element specified in Table 11 and the 360-element mesh (dotted line) is stiff.  It was also verified that the curves for HW51 and the two sets of transformation for the assumed stress/strain: (a) CTV-CTV of Eqs. (33) and (32) and (b) CTV-COV of Eqs. (33) and (34), coincide.

L-shaped plate
The purpose of this non-linear test is to characterize the convergence of the solid-shell elements by comparing the last solution points reached in 10 steps.
The L-shaped plate is clamped at one end and the in-plane force P is applied at the other end. The data is as follows: E = 71240, ν = 0.31, w = 30, L = 240, and the  thickness h = 0.6, see [3]. A mesh of 64 square elements of size 15 × 15 is used and two forces P/2 are applied to the top and bottom nodes to replace the force P, see Fig. 23. The RBF correction of Sect. 4.3 is used. The deformation associated with the primary solution path takes place in the X0Y plane and has a bifurcation point at which an out-of-plane deformation occurs. We add a small out-of-plane load P z at point B to initiate this deformation. (More details related to a precise estimation of the bifurcation load are given in [48] p. 440.) The equilibrium problem is solved using the arc-length method for the prescribed initial load increment P and the requested number of iterations per step I req . The arc-length is modified at the beginning of each load step using the scaling parameter where I prev is the number of iterations in the previous (convergent) load step. First, the non-linear solutions for all the tested and reference solid-shell elements are computed and they coincide with the reference curve (Ref.), obtained for the 4-node shell element HW47. Next, 10 steps are performed using the arc-length method for P = 0.5, P 3 = P/1000, and I req = 25. In Fig. 24, the final solution points for particular elements are indicated. Two elements, our HW51 and the reference HSEE, cover the longest distance in 9 steps but diverge in the 10th step. In 10 steps, the longest distance covers our HW35, then our HW27B and HW19, and finally the reference EAS10.

Compression of a nearly-incompressible block
For the enhanced low-order solid elements and nearly incompressible materials, the hourglassing (zero-energy) modes can occur for large compressive (uniform) strains, see, e.g., [11,13,23,24,58,59]. It is important to check this aspect for 2D 4-node plane strain elements and 3D 8-node solid elements as it indicates an instability of their formulation.
For the solid-shell elements, hourglassing does not seem to be recognized as a problem and, e.g. in [15], which deals with nearly-incompressible materials, this issue is not addressed. This can be explained by the fact that the solid-shell elements are mainly applied to shell structures, in which large compressive strains are prevented by design to avoid the loss of structural stability.
The purpose of this test is to establish whether hourglassing appears or not in our solid-shell HW elements and for this purpose we find the first zero eigenvalue and the associated eigenvector of their tangent matrix. Two orientations of the solid-shell elements are checked so they are compressed in either their in-plane or normal (through thickness) direction.  The block of size 1 × 1 × 1 is modeled by a mesh of 10 × 10 × 1 elements, see Fig. 25, where one element is used in the thickness direction (normal to the plane of figure). The block is supported in the vertical direction at the bottom and vertically compressed by a sequence of increments of displacement v = 0.001 applied at the top. The Newton method is used to solve the equilibrium equations for each increment and 10 lowest eigenvalues of the tangent stiffness matrix are computed at the converged configurations using ARPACK. For the lowest eigenvalue equal to zero, the corresponding scaled eigenvector is plotted, so it can be seen directly whether hourglassing appears.
The modified neo-Hookean hyper-elastic material model is used, with the strain energy function W defined in terms of the volumetric/deviatoric split of the deformation gradient F, where C = F T F is the right Cauchy-Green deformation tensor and J = det F. The material constants are as follows: the shear modulus μ = 20, the bulk modulus K = 4 × 10 5 and the dimensionless parameter β = −2. The RBF correction of Sect. 4.3 is not used and the shear correction factor k = 1. The obtained critical strains for all the tested and reference solid-shell elements are given in Table 12. The reduced representation elements and the reference elements behave similarly and display hourglassing at the first (lowest) zero eigenvalue. The HW51 behaves differently; the zero eigenvalues appear at lower strains (two of them are reported) but there is no hourglassing for the in-plane compression. Nonetheless hourglassing appears for the normal compression of HW51.
For reference, this example is also computed using two 3D 8-node solid elements; our mGu-EADG12B and TSCG12 of [23], both integrated using the 2 × 2 × 2 Gauss rule. Our element uses the mean gradient of displacements (mGu), the γ -vectors and the Enhanced Assumed Displacement Gradient (EADG12B) with 12 parameters. It is slightly faster than TSCG12 but less accurate. Both these elements, are suitable for the material of Eq. (81).
The two lowest eigenvalues of the solid-shell element HW35 and the reference 3D solid elements are shown in Fig. 26. Interesting is the difference between the curves for the normal and in-plane compression of HW35 (thick lines), which is a result of using the techniques of Sect. 4 and characteristic for solid-shell elements. Note that the curves for HW35 and 3D mGu-EADG12B coincide while the curves for 3D TSCG12 are different.
To detect hourglassing at the first zero eigenvalues, scaled eigenvectors are superimposed on the deformed configuration, see Fig. 27 for the in-plane compression. The HW35 and our 3D mGu-EADG12B display similar hourglassing while no hourglassing appears for the 3D TSCG12.
Finally, the reduced representation solid-shell HW elements show hourglassing at the first zero eigenvalue but it appears at the large strain equal to 0.37 for the in-plane compression. This is the level of strain for which the shell structures are unlikely to be designed.

Summary of results related to the RBF correction
The RBF correction was implemented to enhance the elements' behavior in sinusoidal bending and for significant elements' thinness. It was tested on several linear and nonlinear examples, the results of which are summarized below.
For the rectangular mesh, all elements show a small improvement when the RBF correction is used. For the trapezoidal mesh, the errors of the reduced representation HW elements are diminished. The errors of HW51 and the reference HSSE are the largest; for the RBF correction they slightly increase to exceed 5%. 2. Curved 3D cantilever of Sect. 5.2.6. Case 1 with warped elements. For the slenderness log 10 (R/h) = 3, the RBF correction significantly improves the accuracy, reducing the error from 10 to 2% for EAS10, from 30 to 5% for HW51 and from 30 to 10% for the reduced representation HW elements and the reference HSEE, see Fig. 10. Case 2 with skew nodal "directors". All the reduced representation HW elements perform similarly, and the RBF correction has a small negative effect on them. HW51 and the reference HSSE yield excessive displacements of an over 12% error, which remains practically unchanged by the RBF correction. 3. Pinched hemisphere with hole of Sect. 5.3.2. In the linear test, the RBF correction reduces the relative errors and they are smaller than 0.48% for the dense mesh, see Table 8. In the nonlinear test, all the elements perform almost identically. The solution curves follow the reference curve up to a certain load, and then gradually depart from it. The displacements obtained with the RBF correction are slightly excessive for the thickness h = 0.04 and more exact for h = 0.01, see Figs. 17 and 18. In the one-step non-linear test, the RBF correction is particularly beneficial for the reduced representation HW elements, and they all have a much bigger radius of convergence than the reference EAS10, see Table 9.
4. Twisted beam of Sect. 5.3.3. In the linear test, the errors are smaller than 0.67%, and the effects of the RBF correction are ambiguous (positive or negative), see Table 10. In the non-linear test without/with the RBF correction, the curves for all HW elements coincide, and these obtained with the RBF correction are slightly 'softer'. The softest is the reference EAS10, which exactly matches the reference 4-node HW47, see Fig. 20. 5. Short C-beam of Sect. 5.3.4. In the linear test without the RBF correction, the errors are below 3.35%, which slightly increases to 3.46% when this correction is applied. All HW elements perform well, and HW19 is the most accurate, see Table 11. In the non-linear test, the solution curves obtained with the RBF nearly coincide for all elements, and are slightly more stiff than the reference curve. 6. L-shaped plate of Sect. 5.3.5. In this non-linear test, the solution curves are obtained with the RBF correction and they coincide for all elements. Among the reduced representation HW elements, the biggest radius of convergence has HW35, then HW27B and HW19, see Fig. 24; In conclusion, we see that the performance of the reduced representation HW elements is in the majority of tests improved by the RBF correction. The radius of convergence is definitely improved, and is larger than the one of the reference element EAS10.

Final remarks
The paper concerns the development of eight-node Hu-Washizu (HW) solid-shell elements, with particular emphasis on the reduction of the number of internal parameters. Four HW elements are developed and tested. They have correct rank, pass the membrane and bending patch tests and are free from the curvature thickness, transverse shear and volumetric locking.
1. At the outset, an element based on a 51-parameter representation of the assumed stress/strain (HW51) is developed, next also three others (HW35, HW27B and HW19) using the reduced representations with 35, 27 and 19 parameters, respectively, including 3 parameters of the EAS enhancements of the thickness strain. (43 parameters are used in the reference HSEE element [20].) The latter two elements are based on the so-called partial HW functionals, which differs from the approach used in other papers, where the HW functional is needed for all strain components. 2. All the reduced representations of the assumed stress/strain are ζ -independent, where ζ is the thickness coordinate, with the skew coordinates in the in-plane (11, 22 and 12) components. The assumed stress/strain representations involve respectively 13/19, 9/15 and 5/11 parameters, which is a reduction from 18/33 parameters of HW51. 3. For the elements based on the reduced representations, the transformations used for the assumed stress/strain have a significant effect for distorted shapes, see Fig. 6. The contravariant (CTV-CTV) transformations yield more accurate solutions for such elements and, therefore, are applied to all of them. The transformation operators T S and T E are computed at the element center and at the reference surface (ζ = 0), which is slightly different than in the HSEE element, see Sect. 3.4.2. 4. An interesting outcome of the paper is that the elements with a larger number of parameters, such as HW51 and the reference HSEE (43 p), are less accurate than the elements with the reduced representations for trapezoidal throughthickness shapes of elements, see Fig. 7, and Tables 5 and  6. This feature provides a strong argument for using the reduced representation elements.

Fig. 27
Compression of a nearly-incompressible block. Eigenvectors superimposed on the configuration corresponding to the lowest zero eigenvalue: a Solid-shell HW35 (in-plane compression) and 3D mGu-EADG12B. b 3D TSCG12 5. The convergence properties of the reduced representation elements indicate that HW35 has a larger radius of convergence than HW27B and HW19, see Sect. 5.3.2 (Table 9) and Sect. 5.3.5 (Fig. 24). The elements with a large number of parameters HW51 and the reference HSEE (43p) converge only slightly better. All HW elements converge much better than the enhanced strain element EAS10. 6. By comparing the reduced representation HW elements with each other, we obtain that they are nearly equally accurate in the linear tests of Sects. 5 In summary, the developed HW elements with the reduced representations of the assumed stress/strain perform very well, especially for trapezoidal through-thickness shapes of elements (skew nodal "directors"). When a large radius of convergence is required then HW35 is advised. Otherwise, considering efficiency, we recommend the element with the smallest number of parameters HW19.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.