An improved quadrilateral shell element based on the Hu–Washizu functional

In this paper a robust and effective 4-node shell element for the structural analysis of thin structures is described. A Hu–Washizu functional with independent displacements, stress resultants and shell strains is the variational basis of the theory. Based on a previous paper an additional interpolation part using quadratic shape functions is introduced for the independent shell strains. Especially for unstructured meshes this leads to an improved convergence behavior. The expanded element formulation proves to be insensitive to mesh distortion. Another well-known feature of the mixed hybrid element is the robustness in nonlinear applications with large deformations.


Introduction
Nonlinear structural analysis of thin structures requires effective and robust element formulations. Especially the possibility of large solution steps and high accuracy when using reasonable unstructured meshes are desired properties.
To bypass the difficulties caused by C 1 -requirements of the Kirchhoff-Love theory many of the shell models consider transverse shear deformations within a Reissner-Mindlin theory. Low order elements like quadrilaterals using a standard displacement interpolation are characterized by locking phenomena and lead to unacceptable stiff results when reasonable finite element meshes are employed. In shells two types of locking occur: transverse shear locking in which bending modes are excluded and nearly all energy is stored in transverse shear terms, and membrane locking in which bending energy is restrained and energy is stored in membrane terms.
An effective method to avoid transverse shear locking is based on assumed shear strain fields first proposed in Ref. [1], and subsequently extended among others in [2,3]. The assumed strain method has also been applied to approximate the membrane strains, e.g. [4][5][6][7][8][9][10]. The papers show that locking is relieved.
The basis for assumed strain methods are multi-field variational principles. Especially for linear elasticity the Hellinger-Reissner functional is adequate as variational founda-tion for mixed interpolated elements, e.g. [11][12][13]. In case of a nonlinear material law a local iteration for the determination of the physical strains is necessary. Hence, a Hu-Washizu functional with independent displacements, stresses and strains seems to be more appropriate, e.g. [9,10,[14][15][16][17][18]. Within the so-called enhanced strain formulation the independent stresses are eliminated from the set of equations using orthogonality conditions and a two field formulation remains [19]. This approach has been successfully applied for shell problems in a multiplicity of publications.
An important issue within the context of developing a finite shell element is the number and type of rotation parameters on the element. Mostly general shell theories exclude explicit dependence of a rotational field about the normal to the shell surface which leads to a five parameter model (three displacements and two local rotations). Use of five degrees-of-freedom frame requires construction of special coordinate systems for the rotational parameters. Considering the so-called drilling degree-of-freedom leads to a finite element discretization with six nodal parameters. This has some advantages since both displacement and rotation parameters are associated with a global coordinate frame, e.g. [20].
The essential features and new aspects of present formulation are as follows: i. Reissner-Mindlin kinematic assumptions considering finite deformations are applied. The variational formulation is based on a Hu-Washizu functional with independent displacements, stress resultants and shell strains. We choose a strain energy density as a quadratic form in terms of the independent shell strains. The finite element formulation for 4-node elements is specified, where the approximation of the displacements and stress resultants is unaltered in comparison to our previous publication [16]. The following amendments are included. ii. In comparison with Ref. [16] the interpolation functions for the shell strains are expanded. Additionally, quadratic shape functions are introduced which are orthogonal to the constant part of the strains. The orthogonality is important for convergence against the correct solution. Furthermore, a shape factor which considers the deviation of the element geometry from a square is incorporated. It leads to an essentially improved convergence behavior especially for unstructured meshes. iii. The derived mixed hybrid quadrilateral element has 5 or 6 degrees of freedom at the nodes, possesses the correct rank and fulfills the membrane and bending patch test. The element formulation is tested by means of several nonlinear shell problems. It is shown that the expanded interpolation of the shell strains with quadratic shape functions relieves membrane locking in an effective way.

Hu-Washizu variational formulation
Let B be the three-dimensional Euclidean space occupied by the shell of thickness h in the reference configuration. With ξ i we denote a convected coordinate system of the body. The coordinate in thickness direction ξ 3 is bounded by h − ≤ ξ 3 ≤ h + , where h − and h + are the coordinates of the outer surfaces. In the following the summation convention is used for repeated indices, where Latin indices range from 1 to 3 and Greek indices range from 1 to 2. Commas denote partial differentiation with respect to the coordinates ξ α . The coordinate on the boundary = u ∪ σ of the initial reference surface is denoted by s.
The position vectors of the initial and current shell reference surface are denoted as X(ξ 1 , ξ 2 ) and x(ξ 1 , ξ 2 ), respectively. Hence, the displacement vector of the reference surface is defined with u = x−X. A vector field D(ξ 1 , ξ 2 ) with |D(ξ 1 , ξ 2 )| = 1, associated with the initial configuration, is introduced. The unit director d of the current configuration is obtained by an orthogonal transformation of the initial vector D. With x, α ·d = 0 shear deformations are accounted for within the Reissner-Mindlin theory.
The shell is loaded statically by surface loadsp on as well as by boundary loadst and couple loadsm on the boundary σ . The loads are assumed to be independent of the displacements. Hence, the variational foundation using the Hu-Washizu functional with dA = j dξ 1 dξ 2 and j = |X, 1 ×X, 2 | is given. Here, v = [u, ϕ] T contains the displacements u and rotational parameters ϕ, as well as ε and σ denote the independent shell strains and stress resultants, respectively. We assume a strain energy density as a quadratic form W (ε) = 1 2 ε T C ε using the constant elasticity matrix C. The geometric shell strains are organized in the vector where the membrane strains ε αβ , curvatures κ αβ and transverse shear strains γ α follow from the Green-Lagrangian strain tensor The vector of independent stress resultants σ = [n 11 , n 22 , n 12 , m 11 , m 22 , m 12 , q 1 , is defined with membrane forces n αβ = n βα , bending moments m αβ = m βα and shear forces q α .
Introducing θ := [v, σ, ε] T and admissible variations δθ := [δv, δσ, δε] T the stationary condition associated with functional (1) reads With integration by parts and application of standard arguments of variational calculus one obtains the associated Euler-Lagrange equations. These are the static field equations, the geometric field equations and the constitutive equations in , as well as the static boundary conditions on σ , see [15]. The associated finite element equations are iteratively solved applying Newton's method. For this purpose the linearization of the stationary condition (5) is derived with C = ∂ 2 ε W as L [g(θ, δθ), θ] := g(θ, δθ) + Dg · θ Finally, the geometric boundary conditions v =v on u have to be fulfilled as constraints.

Finite element equations
We refer to Refs. [15,16] where details of the finite element formulation are specified. The isoparametric concept for 4-node elements using bilinear functions N I (ξ , η) is applied. For the coordinates of the unit square holds −1 ≤ {ξ , η} ≤ 1. The constant orthonormal element coordinate system is denoted by [t 1 , t 2 , t 3 ] and is computed with the nodal vectors X I , I = 1, 2, 3, 4 as follows The superscript h refers to the finite element approximation of the particular quantity, and commas denote the partial derivative with respect to ξ or η. The matrices 22 (9) cause a transformation of contravariant tensor components to the constant element base system t i . The entries J 0 αβ are the components of J evaluated at the element center. The factors a and b are specified below. Detailed investigations on the use of ansatz functions for contravariant stress and strain components in the framework of a Hu-Washizu functional are contained in Ref. [18].
The finite element approximation of the vector δθ h := [δε h g , δσ h , δε h ] T can be written as ⎡ To avoid transverse shear locking, ansatz functions of the assumed strain method [3] are incorporated in B, see Ref. [16].
The matrix N σ for the interpolation of σ h = N σσ as well as δσ h is chosen as follows where 1 n denotes a unit matrix of order n. The coefficient matrices read T 0 σ = T 0 with a = 2 and b = 1 as well asT 0 σ =T 0 . The constantsξ andη are the coordinates of the center of gravity of the particular element. For rectangular elements holdsξ =η = 0. The parameter vectorsσ and δσ contain 8 parameters for the constant part and 6 parameters for the varying part of the stress field. The interpolation of the membrane forces and bending moments corresponds to the membrane part in Ref. [21]. The original approach for plane stress problems was published in Ref. [22]. Regarding requirements on the interpolation functions to fulfill the patch test and to ensure stability of the discrete system of equations we refer to the discussion in Ref. [15]. The matrix N ε for the interpolation of the independent strains ε h = N εε as well as δε h = N ε δε is subdivided in two parts The number of parameters l = n + k of the second part is specified below. The submatrices N 1 ε and N 2 ε read Here, where the index n ∈ {0, 2, 4, 6, 7, 9, 11} has the meaning that optionally the first n columns are taken. With n = 0 the matrix M m n is omitted. The shape factor c considers the deviation of the element geometry from a square. For this purpose the metric coefficients G αβ of the initial reference surface are evaluated at the element center Hence, c is obtained with the ratio of the eigenvalues of G 0 The factor c has a geometrical meaning. One can show that √ c is the ratio of the semiaxes a and b (a ≥ b) of an ellipse which can be inscribed in a distorted element. For a square element holds c = 1. The matrix M b k , associated with the curvatures, reads Again, the index k refers to the number of columns that optionally are taken. The meaningful parameters are k ∈ {0, 2, 4, 6}. In this work we investigate k = 0 (M b k and associated parameters are omitted) and the complete matrix with k = 6.
Remark Due to the factor j 0 /j, the constant coefficient matrix T 0 ε and the functions (ξ , η, ξ η, ξ 2 η, η 2 ξ ) the integral of N 2 ε over the element domain e vanishes. Thus, the functions are orthogonal to the constant part of the shell strains. The orthogonality is important for convergence against the correct solution. This is shown by means of a numerical example in the next section.
The use of transformation matrix T 0 ε in N 2 ε is in contrast to Ref. [16], where with (T 0 σ ) −T transformations of covariant tensor components are described. The numerical tests show that both versions lead with mesh refinement to the same converged solution. The application of T 0 ε yields for coarse meshes a slightly softer behavior.
The finite element approximation of the external virtual work ofp,t andm leads to Here, numel denotes the total number of finite shell elements to discretize the problem and f a corresponds to the element load vector of a standard displacement method. Furthermore, it holds where k g is specified in detail in Ref. [15]. We insert δθ h = N θ δθ according to Eq. (10) and the corresponding equation θ h = N θ θ into the linearized variational equation (6), which now reads with The integrals over an element domain e of a particular element e are computed numerically using a 2×2 Gauss integration scheme. With incorporation of the quadratic functions in Eq. (15) a 3 × 3 Gauss integration is necessary.
Matrix F is expressed with (12) The last four columns with quadratic shape functions in (15) are not orthogonal to column 9 and 10 of N σ according to (11) and thus lead to entries in F 2 . They are consistently omitted when setting F 2 = 0 in F, f e and f s .
where r denotes the vector of element nodal forces. Since the stress resultants and shell strains are interpolated discontinuously across the element boundaries the parameters σ and ε can be eliminated from the set of equations. This is done applying a standard Gaussian elimination procedure to the system of equations (24), see Ref. [23]. One obtains the tangential element stiffness matrix k e T , the element residual vectorf and (21) reduces to The shell elements possess 5 or 6 degrees of freedom (dofs) at the nodes. At nodes on intersections 6 dofs (3 global displacements and 3 global rotations) and at the remaining nodes 5 dofs (3 global displacements and 2 local rotations) are present. The linear element stiffness matrix possesses with six zero eigenvalues the correct rank. The derived element formulation has been implemented in an extended version of the general purpose finite element program FEAP [24].

Eigenvalue analysis of the element stiffness matrix
At first, we compute eigenvalues and eigenvectors of the linear element stiffness matrix considering different parameters. Following Ref. [10] we examine a square element and a distorted (warped) element with a = 2, h = 0.02, E = 10 8 , ν = 0.3, see Fig. 1. For the parameter n we choose n = 0, 7, 11, whereas k = 0 is set in Fig. 2. The 4-node element has 5 degrees of freedom at each node, thus the element stiffness matrix is of order 20. All versions lead to six zero eigenvalues corresponding to the six rigid body modes. The remaining 14 nonzero eigenvalues are depicted in Fig. 2. For comparison we add results of the +HW element (taken from [10]). Lower eigenvalues relate to the bending modes, and higher eigenvalues relate to the stiffer membrane-and shear modes. Both are divided by a pronounced jump, see   The influence of the parameter k, associated with the curvatures, is depicted in Fig. 4. We present results for n = 11 combined with k = 0 and k = 6. As expected, k = 6 leads to a reduction of the eigenvalues of the bending modes 7 and 8 and additionally of mode 9 for the distorted case.
Finally, the pure membrane case considering a flat element (z ≡ 0) is investigated. The distortion in the x − y plane corresponds to Fig. 1. The out-of plane displacements and the rotations are fixed. The remaining degrees of freedom are 8 in-plane displacements. As Fig. 5 shows, there are 3 zero eigenvalues associated with the 3 rigid body movements of a flat sheet. Using the parameter n = 11 and c according to Eq. (17) yields two further eigenvalues which are almost zero. It leads to an unstable element behavior. This is not the case for c = 0.
Summarizing, the parameter n = 11 leads in comparison with n = 0 and n = 7 to lower eigenvalues 12-15. The parameter k = 6 reduces the eigenvalues of some bending modes. This is the reason for the improved convergence behavior in the subsequent depicted test examples. It holds especially for element geometries which deviate notably from a square. Pure membrane problems can be computed with n = 11 and c = 0.

Hemispherical shell
The next problem is the hemispherical shell with an 18 • cutout subjected to alternating radial point loads P at its equator, shown in Fig. 6a. This geometrically non-linear example is often cited as a benchmark problem for shell elements. It is a test for the ability to   [25]. Geometrical and material data are R = 10, ϕ = 18 • , thickness h = 0.04 and E = 6.825 · 10 7 , ν = 0.3. Considering symmetry one quarter of the structure corresponding to the region ABCD in Fig. 6a is discretized using 8 × 8 and 12 × 12 regular meshes. We employ the boundary conditions u y = β = 0 on AD, u x = β = 0 on BC and u z = 0 at a point on AB, e.g. at A. Figure 7 shows the load displacement curves for the regular meshes. The defined converged solution is computed with a 128 × 128 regular mesh. Results are only presented for P − u xA ; similar output can be obtained for P − u yB . In addition, Fig. 8 depicts results for distorted meshes. The principal mesh distortion is described in Fig. 6b for a 4 × 4 mesh. Each edge is discretized using the aspect ratios L 1 : L 2 : L 3 : ... : L N = 1 : 2 : 3 : ... : N , where N denotes the number of elements per direction. The 12 × 12 distorted mesh is illustrated in Fig. 6c. As can be seen in Figs. 7 and 8, significant improvements can be achieved along with the quadratic terms in Eq. (15) (n = 11), especially for distorted meshes. For comparison we add results from Ref. [8] using the MITC4+ element. A plot of the factor c according to Eq. (17) is shown for the distorted 8 × 8 mesh in Fig. 9. The factor may take values much larger than 1. The importance of the orthogonality of the used functions in (15) and (18) is visualized in Fig. 10. The function ξ 2 instead of ξ 2 η is not orthogonal over the unit square with respect to a constant. In Fig. 10 non-orthogonal 1 means use of ξ 2 and non-orthogonal 2 use of ξ 2 − c as well as in an analogues way for η 2 . One can see that both versions lead to convergence against a wrong solution. Also with c = 4/3, whereby the integrals of ξ 2 − c and η 2 − c over the unit square vanish, one obtains likewise wrong displacements with mesh refinement. The results are computed with regular meshes, where the defined converged solution u xA = −8.1546 is obtained with a mesh of 128 × 128 elements. The final deformed mesh is depicted in Fig. 11.

Cylindrical shell segment
In this subsection we examine a cylindrical shell segment, e.g. [26], subjected to a uniform bending moment M = M 0 · h 3 along BC. The shell segment is fully clamped at DE, see     h = R/10, 000 and E = 2.1 · 10 6 , ν = 0. Figure 13 depicts load displacement curves for point A and regular meshes. A 128 × 128 regular mesh is utilized for the defined converged solution. In addition, Fig. 14 shows results for distorted meshes. The principal mesh distortion is described in Fig. 12b for a 4 × 4 mesh. The curved edges are discretized using the aspect ratios L 1 : L 2 : L 3 : ... : L N = 1 : 2 : 3 : ... : N , where again N denotes the number of elements per direction. A 12 × 12 distorted mesh is presented in Fig. 12c in a perspective view. Again, improvements can be achieved for distorted meshes when using the quadratic terms in Eq. (15) (n = 11). For comparison we add results of the element formulation [10] denoted as +HW. The performance of the MITC4+ element [8] is similar. The final deformed mesh is depicted in Fig. 15.

Twisted beam
We consider the twisted beam problem shown in Fig. 16, originally introduced in [25]. Geometrical and material data are L = 12, b = 1.1, thickness h = 0.0032 and E = 29 · 10 6 , ν = 0.22, respectively. The cantilever beam is clamped at one end and is loaded by an out-of-plane acting load P at point A. A regular 4 × 24 mesh is chosen for the solution. Figure 17 depicts the convergence behavior of the displacements of point A  for different parameters n and results using the MITC4+ element [8]. The converged solution is obtained employing a 32 × 192 regular mesh. Furthermore mesh distortion is investigated. The first distorted mesh is shown in Fig. 18a together with a flat projection in Fig. 18b, both in a perspective view. A ratio L max /L min = 2 is chosen, where L max and L min denote the longest and shortest element length in the flat projection, respectively. Figure 19 depicts the resulting load displacement curves of point A. Very good results can be seen, even for n = 0. In addition, we investigate a second distorted mesh, where the distortion is introduced in the opposite direction, see Fig. 20. The associated load displacement curves of point A for a different choice of n are contained in Fig. 21. Again, the quadratic terms in Eq. (15) (n = 11) are necessary to produce accurate results. The convergence behavior of displacement u yA for the second distorted mesh versus the number of elements N in width direction is presented in Fig. 22. Again, n = 11 leads to a significant improvement of the element behavior. This is only achieved with the shape factor c according to (17), as the comparison with the curve n = 11 (c = 0) shows. The deformed beam using the distorted mesh 2 for P = 4 · 10 −2 is depicted in Fig. 23.

Hook problem
Next, we consider the hook problem shown in Fig. 24, referred to in linear analysis as the Raasch challenge [27]. For the FE-discretization we use N × 2 N × 3 N elements with N   The structure is fully clamped at one end and is loaded by a shear load P applied as a uniformly distributed traction at the free end. For the solution, we use a regular 4 × 8 × 12 mesh. Figure 25 shows the resulting load-displacement curves of point A, where curves for the MITC4+ element [8] are included. Similar results can be found for the +HW element (see Figs. 12b, 13b of Ref. [10]). The defined converged solutions are obtained with a 32 × 64 × 96 regular mesh. The principal distorted mesh is shown in Fig. 26 with respect to a flat projection together with a perspective view of the structure. Here, L max /L min = 1.5 is chosen for the first arch and L max /L min = 2.0 for the second arch according to [10]. Figure 27 shows the convergence behavior of displacement u zA of point A versus the number of elements N in width direction. Results for MITC4, MITC4+ and +HW are taken from Figs. 12a and 13a in Ref. [10]. The superior behavior of the MITC4+ and +HW elements as well as of present element formulation with n = 11 is shown. The deformed regular mesh is depicted in Fig. 28.

Cook's problem
Here, we discuss the influence on the nonlinear analysis of the well-known Cook's membrane, first introduced in [28] for the nonlinear case. It is a tapered panel clamped on one end and uniformly loaded with a resultant P = 1 on the other end, see Fig. 29. Geometrical   and material data are h = 1 and E = 2, ν = 1/3. The problem provides a pure membrane test including element distortions and is a test for handling the in-plane bending dominated by shear. The discretization is performed with a N ×N mesh. Using present element the total load P = 1 can be applied in one step with six iterations. For the displacement u yA we depict the performance in dependence of N in Table 2 and Fig. 30. Results for MITC4, MITC4+ and +HW are taken from Table 3 in Ref. [10]. The output for HW in [10] is identical with our results using n = 0. The convergence behavior of the MITC4 and MITC4+ elements is relative slow. Solutions for +HW [10] and n = 7, 11 exhibit a fast convergence. As shown by means of above computed element eigenvalues the parameter n = 11 can only be used with c = 0. Otherwise the element formulation leads in a pure membrane case to hourglassing. The deformed mesh for P = 1 is depicted in Fig. 31.

Annular plate
This example, which is shown in Fig. 32, has been introduced in [29]. We have presented results for this problem in [16] using also the element formulation [30]. The annular plate is loaded at its free edge with a constant load p = 1; the other edge is clamped.   [16]. Here, 0 • and 90 • refers to the circumferential direction and the radial direction, respectively.
For the material behavior transversal isotropy is assumed. The material and geometrical data are: E 1 = 40 · 10 6 E 2 = 1 · 10 6 G 12 = G 23 = 0.6 · 10 6 ν 12 = 0.25 The analysis is based on a N × 5 N mesh, where N denotes the number of elements in radial direction. Furthermore mesh distortion is investigated. Meshes are generated in polar coordinates using four corner nodes. For the distorted meshes in circumferential direction two intermediate nodes are added. With the nodes C-r and D-r, see Fig. 33, a regular mesh and with nodes C-d and D-d a distorted mesh is generated. The position of C-d and D-d is described in the following. We introduce a rotation angle ϕ in the x − y plane, which is defined via ϕ(ξ ) =  Fig. 34. In both cases we define a converged solution using a 24 × 120 regular mesh.    The convergence behavior of displacement u zB versus the number of elements N in width direction is presented in Fig. 35 for p = 1. One can see that the matrix M b k with k = 6, see Eq. (18), leads to minor improvements of the element behavior. The deformed plate using the regular mesh is depicted in Fig. 36 for p = 1.

Stiffened cylindrical shell
With standard degrees of freedom at the nodes also shell intersection problems can be computed. As an example the stiffened cylindrical shell according to Fig. 37 is considered. The figure shows a cross-section of the shell and a finite element mesh of half the structure accounting for symmetry conditions. Radius and length of the cylinder are R = 1000 mm, L = 2000 mm and the shell thickness is h = 10 mm. The shell is free at y = z = 0 and clamped at y = L. A concentrated load F acts at the coordinates (x, y, z) = (0, 0, R). The skin of the structure consists of a [0 • /90 • /0 • ] lay-up, where 0 • refers to the circumferential direction and 90 • to the y-direction. The stiffeners with measurements d = 50 mm and h = 10 mm are arranged in radial direction. In the symmetry axis a thickness 2 h is present. The stiffeners are homogeneous and the fibre direction coincides with the length direction. The material parameters assuming transversal isotropy are chosen as follows E 1 = 125000 N/mm 2 G 12 = 4800 N/mm 2 E 2 = 7400 N/mm 2 G 23 = 2700 N/mm 2 ν 12 = 0.34 .
The mesh density is denoted by N × M × K , where in Fig. 37 N = 12 is the number of elements in circumferential direction, M = 8 the number in length direction and K = 2 the number of elements for a stiffener in radial direction. Thus, the relations M = 2 3 N and K = 1 12 N ≥ 1 are used. A distorted mesh is depicted in Fig. 38. Here, we apply the same distortion technique for the skin as described before for the annular plate. Now, boundary values y 1 = y(−1) = 0 , y 2 = y(0) = L/2(1 ± d L ) and y 3 = y(1) = L with d L = 0.3 are chosen at each stiffener intersection.
The geometrical nonlinear computations are performed with displacement control and a step size w = 20 mm. The load F is computed as reaction on the prescribed displacement. The load step size w can be enlarged up to w = 100 mm in the first increments. In Fig. 39 load F is plotted versus the prescribed displacement w for the chosen parameters. The convergence behavior for the external load for the prescribed displacement w = 400 mm with respect to the number of elements N is depicted in Fig. 40. Again, the matrix M b k with k = 6 leads to minor improvements of the element behavior. Figure 41 shows that the final configuration is characterized by finite deformation.

Conclusions
Based on a previous paper on a mixed hybrid quadrilateral shell element the interpolation matrix for the shell strains is expanded by quadratic shape functions. Thereby membrane locking can be significantly relieved. The new developments lead to a considerable improvement of the approximation behavior especially when the element form deviates from a square. Based on the performed numerical tests we recommend n = 11 for the membrane part with shape function matrix (15) and k = 0 for matrix (18) as the additional terms for the curvatures lead to minor improvements. For pure membrane problems the shape factor has to be used with c = 0. A well-known feature of present element formulation is the remarkable robustness in nonlinear applications. It allows very large load steps in comparison to element formulations based on the displacement method or enhanced strain elements.