Robust four-node elements based on Hu–Washizu principle for nonlinear analysis of Cosserat shells

Mixed 4-node shell elements with the drilling rotation and Cosserat-type strain measures based on the three-field Hu–Washizu principle are proposed. In the formulation, apart from displacement and rotation fields, both strain and stress resultant fields are treated as independent. The elements are derived in the framework of a general nonlinear 6-parameter shell theory dedicated to the analysis of multifold irregular shells. The novelty of the developed elements stems from the fact that the measures of assumed strains and stress resultants are asymmetric. The original interpolation of drilling and bending components of strains and stress resultants is proposed. In the formulation of new mixed elements, two variants of the interpolation of membrane components are used and interpolation of the independent fields is defined in the natural or skew coordinates. Accuracy and efficiency of the developed elements are tested in several linear and nonlinear numerical examples. It is shown that smaller number of independent parameters in the interpolation of membrane components gives more accurate results. The proposed mixed 4-node elements enable the use of large load steps in nonlinear computations. Moreover, they require significantly less equilibrium iterations than other shell elements formulated in the 6-parameter shell theory.

the shape and topology optimization of shell structures lead to the shells with more complex geometry. The shell finite element is a crucial tool in accurate analysis of such structures that is still a demanding task.
Four-node shell finite elements are commonly used among engineers nowadays. The majority of commercial FEM systems offer such elements as a default tool in analysis of different mechanical problems. As it is well known, the classical pure displacement formulation with fully integrated matrices gives usually too stiff deformation for general shell problems due to shear and membrane locking. In commercial FEM codes, finite elements rely on reduced (or selective reduced) integration scheme to avoid locking phenomena. However, the problem of non-physical deformations of analysed structures appears in solutions obtained by the elements with reduced integration. Different techniques to mitigate spurious zero energy modes were proposed in the literature [1][2][3], but they only limit this effect. An alternative approach to deal with the locking phenomena is application of mixed hybrid finite elements. These elements are based on multifield variational principles and besides displacement field, stress field and/or strain field are treated as independent. The mixed hybrid elements are characterized by high accuracy and robustness and reduce substantially the locking effect; see, for instance, the review work [4]. For the first time, hybrid finite element with assumed stress field was introduced in the paper [5]. The element formulation was based on the complementary energy principle. In further works, see for instance [6], the Hellinger-Reissner variational principle was used in the formulation of hybrid stress elements. Pian and Sumihara [7] proposed 4-node plane stress element less sensitive to geometric distortion and providing accurate stress results.
The mixed shell elements were extensively developed in 1980s and 1990s. The approach of assumed stress resultants for membrane and bending components and the assumed natural strain (ANS) concept [8] for transverse shear components was used in formulation of mixed shell element in papers [9,10]. The mixed hybrid plate and shell elements with assumed strains based on the modified Hellinger-Reissner functional were described in [11,12]. The concept of enhanced assumed strains (EAS) defined on the element level was proposed in [13] and then developed in [14,15]. The works on the mixed hybrid shell elements were continued in twenty-first century. Wagner and Gruttmann presented robust 4-node shell elements based on three-field Hu-Washizu functional in papers [16,17]. The proposed mixed hybrid elements allow for very large increments (e.g. load steps) in elasto-plastic analysis and require significantly less equilibrium iterations in comparison with other shell elements. The effective shell elements with interpolation of assumed stress resultants and strains defined in skew coordinates were proposed in [18]. In these mixed Hu-Washizu elements, the drilling rotation was included in the formulation through the drilling rotation constraint equation. Four-node hybrid shell element based on the modified complementary energy functional almost insensitive to mesh distortions was described in [19]. Different methods of incorporating the drilling degree of freedom (DOF) in formulation of the EAS shell elements were compared in [20]. Mixed shell elements with the drilling DOF using the co-rotational kinematics concept were presented in [21][22][23]. Recently, some new shell and plate elements have been described in [24][25][26][27][28][29].
The mixed hybrid elements discussed here are developed in the framework of nonlinear 6-parameter shell theory, see for instance [30,31]. A characteristic feature of the nonlinear 6-parameter shell theory is the lack of symmetry of strain and stress resultant measures. The sixth parameter is the drilling rotation that is included in the formulation in the natural way. The theory is especially dedicated to analysis of irregular shell structures, containing, for instance, orthogonal self-intersections. From the theoretical viewpoint, the approach is statically and geometrically exact, since the preliminary postulates and simplifying assumptions are not needed. A kinematic model of the shell reference surface is formally equivalent to the Cosserat surface [32,33] with three rigidly rotating directors. The only assumptions are made during formulation of constitutive equations. Some of the recent results in this respect can be found in [34][35][36][37][38].
For the first time, hybrid plane stress elements with the assumed asymmetric stress field were proposed in papers [39,40]. However, these elements formulated based on the complementary energy principle yielded too stiff response in numerical tests, because of too rich assumed interpolation. Three new variants of interpolation of membrane stresses were described in [41]. However, only the variant with constant interpolation of the membrane shear components was used in the formulation of further two-dimensional assumed stress hybrid elements [42,43]. Sansour and Bednarczyk [44] proposed 4-and 9-node mixed shell elements with asymmetric assumed stress resultants. These elements were defined on the Cosserat surface and based on the two-field Hellinger-Reissner principle. The mixed hybrid shell elements with non-symmetrical interpolation of assumed strains and enhanced assumed strains were formulated in [45].
In the framework of the present shell theory, mixed and semi-mixed elements based on two-field variational principles were proposed in [31,46] and reformulated using the EAS concept in [47,48]. In [48], the ANS approach was applied for transverse shear components. In the shell elements mentioned above [44,45,47,48], only the membrane and transverse shear components were treated as independent on the element level. Moreover, interpolation of assumed strain and stress resultant fields was defined only in the natural coordinates. The present paper introduces new approach to the idea of mixed shell elements within the framework of considered shell theory. The essential features and original aspects of the proposed shell elements are: • For the first time, the mixed hybrid elements with asymmetric independent strain and stress resultant fields are formulated. The elements are based on the three-field Hu-Washizu functional, and here for the simplicity of considerations a linear elastic material model is used. The first part of strain interpolation is the same as for the stress resultants, while the second part is assumed orthogonal to the stress resultant shape functions, like in the enhanced strain formulation. The stress resultant and strain parameters are effectively eliminated on the element level. • In this paper, the interpolation of asymmetric strains and stress resultants is defined in the natural and skew coordinates. The original interpolation matrices for drilling and bending components are proposed. The second part of strain field approximation from paper [17] is modified and generalized to the 6-parameter shell theory. The proper selection of interpolation functions for independent stress and strain fields plays crucial role in the element formulation. Here, two variants of interpolation matrices for asymmetric membrane components are tested. • The accuracy of the proposed mixed shell elements is verified in several shell benchmark problems. The special attention is paid to examples with irregular mesh. The performed analyses show that the developed elements are robustness in nonlinear computations. They allow for large load steps and require significantly less equilibrium iterations than other elements formulated in the 6-parameter shell theory.
The family of mixed hybrid shell elements in the 6-parameter shell theory was introduced in doctoral thesis [49]. The preliminary results for two initial variants of mixed elements were presented in the conference proceedings [50]. Here, enhanced assumed strain approach is used only for membrane components, and new interpolation variant of drilling components is proposed. Moreover, detailed description of formulation of the mixed hybrid elements and the results for several numerical examples are presented.

Hu-Washizu functional in the nonlinear 6-parameter shell theory
The general nonlinear 6-parameter shell theory is developed in the mixed approach that is neither purely derived nor direct and was described, for example, in [51,52]. The theory accommodates naturally finite displacements and finite rotations and has the same kinematical structure as Cosserat rods, see [53]. The strong ellipticity condition and its weak form known as Hadamard inequality were discussed for the nonlinear 6-parameter shell theory in [54]. The formulation, implementation and computational aspects of the present shell theory have been extensively studied in [30,31,46,48,[55][56][57] and references given there. Therefore, in this section only some aspects of the theory necessary to formulate the mixed hybrid elements are discussed. The special emphasis is placed on the Hu-Washizu principle serving as the basis for the weak statement of the developed elements.
The shell reference surface of Cosserat type is denoted by M in the undeformed configuration (Fig. 1). The shell boundary ∂ M consists of the part ∂ M f with prescribed boundary tractions n * and couples m * , and the part ∂ M d with imposed boundary displacements. Typical point on the shell base surface M is described by the position vector x and the structure tensor T 0 ∈ SO (3). The orthonormal basis t 0 i is defined by the tensor T 0 , as a result of the orthogonal transformation of some global fixed basis where ζ β are arc-length coordinates associated with the unit vectors t 0 i . The deformation of the shell reference surface u = (u, Q) is described by the translation vector u(x) and the proper orthogonal tensor Q(x) ∈ SO(3) specifying independent rotations of the vectors t 0 i . Hence, the position vector y(x) and the structure tensor T(x) ∈ SO(3) at the current base surface of shell are computed in the following way (see [30,31,46]) (3) Fig. 1 The shell kinematics in the nonlinear 6-parameter theory During deformation, the directors t 0 i rotate to its current position We assume that the surface M (excluding intersection edges) is smooth enough to perform necessary mathematical operations, such as computation of the metric and curvature tensors. Correspondingly to the field of generalized displacements u = (u, Q), there exists the field of virtual displacements w = (v, w). The vector fields v and w may be interpreted as the kinematically admissible virtual translations and rotations, respectively, such that v = w = 0 on ∂ M d . They are defined as (see, e.g. [58]) where δ is the symbol of virtual change (variation) and ax(W) is the axial vector of the skew tensor W. In order to parametrize the rotations, the canonical parametrization [59] is used for the tensor Q(x) ∈ SO(3) and the vector w is used in the spatial representation. The natural strain measures on the Cosserat surface characteristic for the nonlinear 6-parameter shell theory [30,31,46] are defined as The stretching vector ε β and the bending vector κ β in symbolic notation form together a column vector ε(u) of the compatible shell strains computed from the generalized displacement field u The virtual measures of the compatible strains (6), (7) are calculated as (see [30,31,46]) They are equal to the co-rotational variations of the strain measures [31,46,60] For convenience, equations (9) are put into the symbolic matrix notation where the matrix operator B describes the relation between the virtual strains and the virtual displacements. For the proposed mixed elements beside the vector of compatible strains ε(u) also a column vector of independent strains ε is defined on the element level. The following formal decompositions of the strain vector and the strain-displacement operator are introduced where subscripts m, s, b, d denote the membrane (in-plane), transverse shear, bending and drilling components, respectively. Characteristic features of Cosserat-type theories [61] are: asymmetric membrane shear strains ε 12 , ε 21 , asymmetric bending strains κ 12 , κ 21 and the presence of drilling strains κ 1 , κ 2 . The vector of independent stress resultants and couples, that is energy conjugated with the strain vector ε, is assumed in the form The vectors of external loads in the matrix notation read where f (x) and c(x) are the external resultant force vector and resultant couple vector, respectively (see Fig. 1).
To simplify presentation and interpretation of further derivations of formulae, the existence of a 2D strain energy density Φ(ε) = 1 2 ε T Cε on the reference surface is postulated. Under this assumption, the constitutive relation reads The formulation of the proposed mixed hybrid shell elements is based on the three-field Hu-Washizu principle [62] where V (u) is the assumed potential of the external loads (15) given by formula The weak form of the boundary value problem may be posed as follows: find kinematically admissible field u and independent fields ε, s, such that W (u, ε, s) → stationar y. In order to simplify notation, the following vectors of element variables are introduced θ = [u, ε, s] T , δθ = [w, δε, δs] T and θ = [ u, ε, s] T . Then, the first variation of the Hu-Washizu functional (17) may be written in the form Linearization of (19) gives the second variation of the Hu-Washizu functional (17) Fig. 2 The image of the standard element π (e) into shell element Π (e) on the reference surface In the above equation, the external loads are assumed as dead loads, i.e. δ 2 V [u, w, u] = 0, so the stiffness matrix associated with loads does not appear.
In the reminder of paper, linear elastic constitutive equations, described for instance in [30,31,[46][47][48], are assumed, because a linear elastic material model was used in the majority of benchmark problems proposed as the tests of shell element formulations.
The interpolation scheme described above assures the required C 0 continuity of the displacement field between elements; therefore, the shell reference surface may be spatially approximated as a sum of the finite elements where N e is the total number of shell finite elements. Relation between the parent (natural) coordinates ξ = (ξ 1 , ξ 2 ) ∈ π (e) = [− 1, +1] × [− 1, +1] and the arc-length coordinates ζ = ζ 1 , ζ 2 is described by the Jacobian matrix J Matrix J can be interpreted as the transformation matrix between the basis vectors t 0 α given by (2) and the basis vectors g α of the parent coordinate system ξ , see Derivatives of the shape functions L a (ξ ) with respect to the arc-length coordinates are computed considering (22) as follows The finite element area A (e) is defined in the coordinates ζ as where j (ξ ) is a determinant of the Jacobian matrix J and can be expanded as follows Finally, the element area may be computed based on the value of the Jacobian determinant j 0 in the element centre A (e) = 4 j 0 . For SO(3)-valued variables, the interpolation scheme described above cannot be used. Hence, in the paper these values are interpolated using the multistep approach, described for instance in [46,47]. The virtual rotation vector w ∈ so(3) is interpolated directly in the global fixed basis {e i }. Using the transformation T : e i → t i , the vector w is written in the form Then, the interpolation of w is performed as follows where T(ξ ) is a transformation matrix corresponding to the tensor T (3).
To minimize shear locking effect, the ANS concept [8] is employed for the transverse shear part ε s (u) of the compatible strains. At selected points A(0,1), B(− 1,0), C(0, − 1), D(1,0) defined in the standard element ( Fig. 2), the following sets of values are computed, cf. (12), (13) The label DI indicates the direct interpolation, i.e. the compatible strains. Then, these values are transformed to the parent element space In the next step, interpolation is performed using the ANS approach Finally, the resulting values are transformed back to the current basis t i in the spatial representation, and appropriate terms in Eqs. (8) and (11) are substituted by the new values.

Skew coordinates
The use of an appropriate coordinate system during approximation at the element level is one of a key aspect in the formulation of a robust 4-node mixed element. Specific modification of the natural coordinates may improve the accuracy of the mixed elements of irregular shape, see for instance results of numerical tests in [63,64]. The skew coordinates were proposed in [65] and then applied in the formulation of different mixed elements, see, e.g. [18,63,66]. These coordinates are defined in the fixed basis at the element's centre In the above formulae, the label a denotes the node and x is the position vector expressed in the Cartesian basis. The position vectorx computed relative to the element's centre x c may be written using the arc-length coordinates ζ 1 , ζ 2 and the skew coordinates x s , y s as The relation between the arc-length coordinates and the skew coordinates can be expressed as where J 0 (ξ 1 = 0, ξ 2 = 0) is the Jacobian matrix (22) calculated at the element's centre. After some further transformations, see, for example, [18,63], the skew coordinates are expressed in terms of the natural coordinates as follows The coefficients j 0 , j 1 , j 2 appearing in Eqs. (26) and (35) are computed as where a i and b i are calculated based on the arc-length coordinates of the element nodes. The coefficients a i , b i appear also in the Jacobian matrix (22) written in the following way The natural and skew coordinates will be applied in interpolation of the assumed strain and stress resultant fields defined in the subsequent paragraphs. In order to simplify the notation, we introduce common symbols ξ * 1 , ξ * 2 for these two types of coordinates. For finite elements in the shape of parallelogram j 1 = 0 and j 2 = 0, thus the natural and skew coordinates are equal for such elements. Consequently, the influence of choice of coordinates on the results will be analysed only in the numerical examples with irregular, distorted mesh.

Interpolation of the stress resultants
The independent stress resultants (14) are approximated by the relations where vector α (e) contains 12 real parameters for the constant part and 7 or 9 real parameters for the varying parts of the stress resultant field, depending on the variant of mixed element. Here 2 options of the membrane components interpolation matrix S mi are analysed Similar interpolation scheme defined in the natural coordinates was proposed in [41]. The varying parts of the transverse shear forces and the bending moments are defined by the following interpolation matrices Matrix S b is an extension of the interpolation method of the bending stress resultants from [16,17] on the nonlinear 6-parameter shell theory. In case of the drilling stress resultants, 4-parameter interpolation variant was introduced in [50]. However, authors' own analyses show that the application of 3-parameter interpolation variant allows us to obtain the results with the same precision. Hence, here to reduce the number of the independent parameters the following interpolation matrix S d is proposed In Eqs. (40)- (42), the transformation matrices are used where J 0 αβ = J αβ (ξ = 0, η = 0) are the components of the Jacobian matrix (22) computed at the element's centre. The matrices T 0 σ ,T 0 σ describe the transformation of the contravariant stress resultant components to the physical coordinate system t 0 i at the element's centre.

Interpolation of the independent shell strains
Similarly as in [17,67], we assume that the interpolation of the independent strains (12) consists of two parts where R is the set of real numbers. The approach discussed above for the stress resultants is used also in the first part of interpolation of the assumed strain field where Taking into consideration the results from paper [68] obtained for the Hu-Washizu mixed elements, the contravariant transformation rule is used also in the interpolation of the first part of the independent strain field As a consequence of assumption (48), the matrices given by (46) and (47) have the same form as the matrices specified in (40) and (41). The covariant transformation rule is used only in interpolation of the assumed drilling strains that allows to satisfy the following linearized compatibility condition, see [49,69] The second part of the independent strain field approximation (44) is assumed L 2 orthogonal to the stress resultant field. In this context, the interpolation scheme corresponds to the EAS formulation [13,48] and mixed-enhanced methods [70,71]. The following form of the interpolation matrix P 2 is assumed In comparison with the previous papers [17,49,50,67], the nonzero interpolation scheme is used only for the membrane components. This approach reduces the number of β 2 (e) parameters. The second part of the strains interpolation is defined only in the natural coordinates to assure fulfilment of the following L 2 orthogonality condition also for the constant stress resultants. We propose the matrix N 2 (52) as the generalization of the matrix M 2 from paper [17] to the theory with asymmetric measures of strains. The analysis performed in [49] indicates that variants of matrix N i with greater number of parameters β 2 (e) give almost the same results as the matrix N 2 .

Linearized variational formulation
In the present shell theory, the vector of nodal degrees of freedom defined in the global coordinate system {e i } at single node a takes the following form where v i are virtual translations andw i are virtual parameters of rotations (28). The vectors δq a are collected into the displacement vector of a 4-node element In the finite element, the subsequent interpolation scheme for the vector of virtual displacements w and the vector of displacement increments u is used whereL is an element matrix of the shape functions L a (ξ ) The vectors of virtual compatible strains δε(u) and strain increments ε(u) after taking into account relations (11), (56), (57) are written as follows The linearized stationary condition for the Hu-Washizu principle (17) may be written as where the first and the second variation of the Hu-Washizu functional are given by (19) and (20), respectively. The substitution of Eqs. (38), (44), (58) describing interpolation of the stress resultant, strain and displacement fields into the stationary condition (59), and given the arbitrariness of the virtual fields, yields the finite element approximation In the set of Eq. (60), the following element matrices and vectors are defined Here K

(e)
G is a geometric matrix described in more detail in "Appendix A". The vector of element external loads p (e) (15) is computed like in the standard displacement formulation.
After static condensation of the parameters β 2 (e) , the set of Eq. (60) takes the form where the following new element matrices and vectors are introduced The integrals in (61), (62), (64) and (65) are computed numerically using the full (2 × 2) integration scheme.
In the next step, the independent parameters α (e) and β 1 (e) are eliminated on the element level, since the assumed strains and stress resultants are approximated discontinuously between finite elements. The effective method of elimination of the strain and stress resultant parameters (described in Appendix D of paper [17]) is used, since equal number of parameters α (e) and β 1 (e) is assumed here. In this order, the vectors α (e) and β 1 (e) are computed from the second and third equation of system (63), respectively and then substituted to the first equation of (63). Finally, the equilibrium equation for the mixed finite element may be written in the standard form where tangent element stiffness matrix K

(e)
T and element residual vector r (e) are introduced The global stiffness matrix K T and the global vectors p, r are obtained from the element matrices and vectors in the standard aggregation procedure. Then, the global set of equations K T q = p − r is solved using an incremental-iterative scheme. The displacements u and the independent parameters α (e) , β 1 (e) are updated in each increment and iteration. On the global level, the translations u are updated additively, while the rotation tensor Q is multiplicatively accumulated in the spatial representation, see for instance [46,47] The independent strain and stress resultant parameters are updated on the element level making use of Eqs. (66) and (67). Since in this study a linear elastic material is used, thus in each iteration it holdsr 1 β = 0. This facilitates the update process of the independent parameters β i+1 where the following matrices and vector are introduced Taking into account (73) and (74) for linear elasticity problems, only the matrices P βq , P αβ and the vectors β 1 ,p β have to be saved in the memory of the authors' FEM code.

Proposed mixed hybrid shell elements
Taking into account two variants of the membrane components interpolation, two types of hybrid mixed shell elements are proposed in Table 1. The element code is assigned based on the total number of the independent strain and stress resultant parameters. The 8-parameter (8p) interpolation of the membrane components means that 4 parameters for the constant part and 4 parameters for the varying part are assumed. While in the 6-parameter interpolation only 2 parameters for the varying part are used. The 8p interpolation is described by the matrices P m4 , S m4 and the 6p interpolation by the matrices P m2 , S m2 , respectively. In the MIX40, MIX44 elements, the natural coordinates are used in the interpolation of the assumed strain and stress resultants. The elements with interpolation defined in the skew coordinates are denoted by letter 's': MIX44s, MIX40s.

Numerical examples
The linear elastic constitutive relations, described in detail in other papers [30,31,47,48], are used in numerical examples. A characteristic for the nonlinear 6-parameter shell theory, the drilling couples M β are calculated from the following constitutive equations where D is the shell bending stiffness and α t is the drilling correction factor, that may be interpreted also as a material parameter, see, e.g. [72]. Here, the following values (see [73]) of the shear correction factors α s = 5/6 and α t = 0.7 are assumed. In this paper, results for the developed mixed hybrid element are compared with the solutions from literature, especially with the mixed 4-node HW47, HW29 elements [18] and the semi-EAS-ANS 4-node EANS4 element [48]. The authors' results obtained with 16-node displacements-rotations-based elements CAMe16 [30] are assumed as basic reference solutions, since they are almost free of the locking phenomena.

Eigenvalues of the tangent matrix
Prior to nonlinear calculations, the spectrum of the stiffness matrix is studied. The eigenvalues and eigenvectors were computed for a single finite element without imposed boundary conditions to avoid masking of the presence of zero eigenvalues. The finite elements of a square and of a deltoid shape are analysed, see Fig. 3. Similarly as in [66], the following material properties E = 10 6 , v = 0.3 and the shell thickness h 0 = 0.1 are assumed.
All developed mixed hybrid elements have a correct rank, since six zero eigenvalues, corresponding to the rigid body motions, were obtained for each element. The nonzero eigenvalues of different mixed shell elements of deltoid shape are compared in Fig. 4. A variant of the membrane components approximation has significant influence on the eigenvalues (see the range 6-8 in Fig. 4), while the influence of type of coordinates is relatively small. The minimally greater eigenvalues were determined for the elements with interpolation of strains and stress resultants defined in the natural coordinates (MIX44, MIX40) for the majority of eigenmodes. Moreover, all nonzero eigenvalues are positive.
which generate the constant strain field and the constant stress field in the shell plane. The proposed mixed elements pass the membrane patch test, since the values of displacements computed at the internal nodes agree with the values from Eq. (76). In passing it is worthy to note that special variants of the membrane patch test for elements with the drilling rotation were proposed in [66]. They verify the presence of additional nonphysical drilling rotation field. In the most demanding variant of this test, the drilling rotations are free at the nodes 1-4. The proposed mixed elements fulfil this test, as the values of the drilling rotation are zero at all nodes.
Consequently, the state of the constant strains and bending moments was obtained in all developed elements. Moreover, the values of deflections at the internal nodes (A, B, C, D)  Table 2 Comparison of values of displacement v A and drilling rotation ϕ A computed in node A of the Cook's membrane (see Fig. 6) for different meshes and finite elements

Cook's membrane
The example was proposed in [75] as a demanding test of 2D finite elements under in-plane shear loading.
Here, the Cook's test is used to evaluate a quality of the developed shell elements in the in-plane bending. A characteristic feature of this example is the skew and tapered shape of finite elements, see Fig. 6. The material parameters are: E = 1, v = 1/3. The membrane of the thickness h 0 = 1 is fixed on the left edge and uniformly loaded on the right edge by the total force P = 1. A discretization 4 × 4 and geometrical dimensions are presented in Fig. 6. The computations were performed also for a coarse mesh 2 × 2 and a fine mesh 32 × 32. The values of the vertical displacement v A and the drilling rotation ϕ A at the middle node A on the right edge (see Fig. 6) computed by the authors' mixed elements are presented in Table 2. The greatest value of v A was determined for the MIX40s element. This value is in a very good agreement with the reference values from papers [7,18,76]. Application of the skew coordinates in the interpolation of the assumed stress resultants and strains gives more accurate results compared to the natural coordinates. A slight locking effect is observed in the solutions obtained by the elements with the 8-parameter interpolation of the membrane components (MIX44, MIX44s). The variant of interpolation of the membrane shear components has the greatest influence on the value of ϕ A . Table 2 shows that also the type of coordinates affects the drilling rotation value.

L-shaped frame
The L-shaped frame is the first example analysed here in a geometrically nonlinear range. It was used originally in [77] and then applied in many papers (e.g. [16,18,48,78]) as a test of the bending part of shell element formulations. The following geometrical dimensions: L = 240, b = 30, h 0 = 0.6 (see Fig. 7) and the material parameters: E = 7.124 × 10 4 , v = 0.31 are assumed. The frame is clamped at the left end and loaded by in-plane force P at point C.
The introduction of the imperfection load P impf = P/1000 at point D (Fig. 7) allows to compute the nonlinear buckling loads and the postbuckling paths. Three FE meshes: (16+16)×2 (16 FE in the leg axis, 2 FE on the width of leg), (32+32)×4 (Fig. 7) and (64+64)×8 are used in computations. Similarly as in [78], two variants of load are investigated: variant A presented in Fig. 7 P A = P and variant B with the forces P B = − P and P B impf = − P impf acting in the opposite direction. The critical buckling forces P cr computed in nonlinear stability analyses using the developed mixed elements are presented in Table 3 for different FE meshes. The values of P cr are almost the same for all authors' mixed elements and are in a very good agreement with the reference values from papers [16,78]. Compared to the basic reference solution (CAMe16), the critical loads for the mixed elements are slightly more accurate than values for the EANS4 element. The nonlinear postbuckling paths computed for a (16 + 16) × 2 mesh are presented in Fig. 8. The increment P = 1 is used in the first step and the arc-length method in the following steps.
A lower number of large steps (see Fig. 8) than in papers [16,18] show the efficiency of the proposed mixed elements. The computations performed with the MIX44 element needed 5 steps (17 iterations), whereas the MIX40 element required only 4 steps (16 iterations). The developed mixed elements are more efficient in comparison with the HW29, EANS4 elements that needed 67 and 73 equilibrium iterations, respectively.

Twisted beam
The example of the clamped beam twisted about 90 o and loaded by concentrated force at the tip was originally proposed in [74]. This test is used to evaluate the influence of warping on the performance of the proposed mixed shell elements. We analyse the most demanding version of this test [18] with a very small shell thickness h 0 = 0.0032; hence, shell elements are more prone to the membrane locking. The other geometric dimensions are as follows: L = 12, b = 1.1, see Fig. 9. The material properties: E = 2.9 × 10 7 , v = 0.22 and a 4 × 24 finite element mesh (Fig. 9) are assumed in geometrically nonlinear analyses. Despite the warped shape of finite elements, the authors' analyses show that type of coordinates has almost no influence on results; thus, only solutions for the natural coordinates are presented below.
The nonlinear equilibrium paths for the developed mixed elements are compared with the reference curves in Fig. 10. Additionally, the computed values of translations u A , w A at point A are reported in Table 4 for the three levels of load. Some discrepancies between the authors' and the reference solutions (see Fig. 10, Table 4) are caused by more than 16 times greater number of nodes in a 6 × 36 mesh of the CAMe16 elements. These differences almost disappear after double mesh refinement, see curves for MIX40 (8 × 48 mesh) in Fig. 10. The load-deflection paths P x (u A ) (Fig. 10) show that the MIX40 element yields more accurate results than the MIX44 element, which gives almost the same results as the EANS4 element.
Then, the convergence rate of the developed mixed elements is evaluated. The maximum load increment P (with accuracy 0.001) and the total number of iterations were determined in the nonlinear analysis up to the force P x = 0.01. The values for different shell element formulations are compared in Table 5. The developed mixed elements allowed to use significantly greater P in comparison with the EANS4 and CAM16 elements.   They require also meaningfully less equilibrium iterations than all reference elements. Table 5 shows that the elements with 6-parameter interpolation of the membrane components are slightly more effective, because of the lower number of iterations.

Pinched clamped cylinder
In this example proposed in [79], the influence of mesh distortion on the results is investigated for shell with nonplanar geometry. In this order, regular and irregular 8 × 16, 16 × 32 (Fig. 11) meshes with two times greater number of elements along cylinder axis than in radial direction were generated. The cylindrical shell is loaded by two opposite forces P at free end and fixed on the other end. The double symmetry of the geometry, boundary conditions and load is exploited in computations, see Fig. 11. The following material properties: E = 2.0685 × 10 7 , v = 0.3 and geometric dimensions R = 1.016, L = 3.048, h 0 = 0.03 are assumed. The values of deflection w A computed in the linear analysis (P = 1000) for different shell elements and FE meshes are presented in Table 6. The type of coordinates has no influence on the value of w A for regular mesh. In the case of irregular discretization (Fig. 11b), application of the MIX44 elements gives slightly more accurate deformation than the MIX44s elements, while the closest value of w A to the reference solution for 8×16 discretization was obtained for the MIX40s element.    Fig. 12. Some discrepancies between the curves obtained for the developed mixed elements and the reference curves almost disappear after double mesh refinement, see Fig. 13. The equilibrium paths λ(u B ), λ(−w A ) (Fig. 12) show that the MIX40 element gives slightly more accurate deformation than the MIX44 and EANS4 elements.
The influence of mesh distortion on the values of load factor λ for chosen values of deflection w A is presented in Table 7. The values computed by the proposed mixed elements for irregular discretizations 8×16  and 16×32 are compared with the reference values obtained for regular 8×16 mesh of the CAMe16 elements. The locking effect and greater differences between values of λ for specific shell elements are observed for irregular mesh compared to the regular mesh. The smallest values of λ, the closest to the reference solution, are obtained for the MIX40 element, except value of λ computed for mesh 8 × 16 and w A = − 0.5. Table 7 shows that mixed elements MIX44 and MIX40 give slightly more precise solution for mesh 16×32 than the corresponding elements with interpolation defined in the skew coordinates.

Pinched hemisphere with a 18 • hole
The hemispherical shell with a 18 • hole is analysed as well-known example (see, e.g. [16,18,44,47,48,78]) of doubly curved shell. The characteristic features of this benchmark test are large almost in-extensional deformation, small strains and trapezoidal-shaped finite elements. Compared to the original version of this example proposed by MacNeal and Harder [74], here a four times smaller shell thickness h 0 = 0.01 is assumed to make shell even more prone to the membrane locking. We analyse only the quarter of the hemisphere with a radius R = 10 (see Fig. 14), due to the double symmetry of the task. The boundary conditions and the external forces are presented in Fig. 14. The material data are as follows: E = 6.825 × 10 7 , v = 0.3. The values of displacement w A at point A computed in the linear analysis for force P = 1 are presented in Table 8. The translations obtained for the developed mixed elements are compared with the reference values for different meshes. The elements with the interpolation of the assumed stress resultants and strains defined in the natural coordinates yield more accurate deformation for coarse meshes than the elements formulated in the skew coordinates. Table 8 shows also that deformation computed by the MIX40 elements is almost not corrupted by the locking effect, despite the very thin shell thickness.  Then, the mesh convergence analysis is performed in a geometrically nonlinear range. The values of translation w A computed for the total force P = 4 and different shell element formulations are compared in Fig. 15. The developed mixed elements enable the faster convergence to the reference solution than the CAMe16 element. The 16×16 FE mesh (see Fig. 14) is assumed in the further nonlinear analyses to compare results with the reference solution from the literature and to expose differences between the proposed elements.
The equilibrium paths obtained for the authors' mixed elements are compared with the reference solutions in Fig. 16. The load-deflection curves for the MIX40, MIX40s elements are in a very good agreement with   Table 9. The results of the linear and nonlinear analyses (Tables 8,9) indicate that the elements with poorer interpolation of the membrane components (MIX40, MIX40s) give more accurate results than other mixed elements. The differences between the solutions for the mixed elements and the CAMe16 element are caused by the very low shell thickness, and they almost disappear for 32×32 FE mesh, see curves for the MIX44 element in Fig. 16.
To evaluate the overall efficiency (robustness) of the proposed mixed elements, the radius and rate of convergence of the Newton method is measured. Namely, the maximum load increment P 1 in the first step of nonlinear analysis is identified for each element with the accuracy of 0.05. The value of load increment P 1 and number of iterations in the first step are compared with the values for the reference shell elements in Table 10. The maximum load increment is observed for the MIX40, MIX40s elements. Then, the efficiency of the mixed elements is tested in the nonlinear analysis for the total force P = 8.8. For each element, the analysis is performed with minimal number of increments and the corresponding maximum step P. The value of P and number of iterations in all increments are presented in Table 11. The values in Tables 10  and 11 show that all developed elements allow for slightly greater step increments and require less Newton iterations than the HW47 element. Moreover, the convergence rate of the proposed elements is significantly better compared to the CAMe16 and EANS4 elements.

Channel section beam
The short-channel section cantilever beam is analysed as an example of the shell with orthogonal intersections that was originally introduced in paper [30]. A characteristic feature of this benchmark is the issue of local and global loss of stability and coupling of the membrane stress state with the bending and torsional stress state. The beam is subjected to the concentrated force P = λP ref (P ref = 100) at point A (see Fig. 17a) and fully clamped at the other end. The following geometric dimensions: L = 36, b = 2, d = 6 (Fig. 17a) and the shell thickness h 0 = 0.05 are assumed. The material data for the channel section beam are as follows: E = 10 7 , v = 0.333. Figure 17a shows the assumed (2 + 6 + 2) × 36 FE mesh consisting of square elements of size 1×1. The computations were performed also for irregular mesh generated by introduction of in-plane distortion d = 0.25 in the case of internal nodes, see Fig. 17b. The nonlinear load-deflection curves for the displacement u A at point A computed for regular mesh are presented in Fig. 18. The equilibrium paths are shown only for the MIX40 and MIX44 elements, since the results for the natural and skew coordinates are equivalent for square finite elements. The displacement control analysis ceased to converge for translation u A = 4.4 in the case of MIX44 element and u A = 5.6 for the MIX40 element. The solutions obtained for the mixed elements are slightly locked in comparison with the reference curves for the HW47 and CAMe16 elements. On the other hand, they are minimally more accurate than solution for the EANS4 element. A discrepancy between the results is the effect of the coarse mesh that does not allow to model precisely local buckling form of deformation at the channel top flange, see Fig. 19. For the same number of degrees of freedom, a very good agreement with the basic reference solution is obtained, see curve obtained for the mesh (6 + 18 + 6) × 108 of the MIX40 elements in Fig. 18.
Extreme values λ extr of the load factor and values of λ for given values of displacement u A are compared in Table 12 for two types of discretization and different shell elements. In the initial phase of deformation   (u A ≤ 1.0), the mesh distortion has relatively small impact on value of λ. In the postbuckling range, the elements MIX44s and MIX40s yield slightly more accurate results compared to the MIX44 and MIX40 elements, respectively. This indicates that in this example the skew coordinates should be used when distorted mesh is considered.

Conclusions
In the paper, original mixed hybrid shell elements with the asymmetric assumed strains and stress resultants are formulated based on the Hu-Washizu principle. The new elements extend the library of elements implemented in the nonlinear 6-parameter shell theory. The accuracy and correctness of the proposed shell elements are positively verified and tested in various numerical examples. The obtained results show that the MIX40s and MIX40 elements are less prone to the membrane locking than the MIX44s and MIX44 elements. Moreover, the MIX40s and MIX40 elements are more accurate than the EANS4 element and give comparable results to the mixed HW47 element. The computations performed for irregular meshes show that the elements with interpolation defined in the skew coordinates give more accurate results in the examples with planar geometry (Cook's test, channel section beam). The application of the natural coordinates is beneficial in the examples of curved shells: pinched hemisphere and pinched clamped cylinder, but only for 16×32 mesh. However, the type of coordinates has relatively small influence on the results, because the difference between solutions obtained for the natural and skew coordinates is greater than 1% only in the pinched hemisphere test for coarse meshes 6×6, 8×8. The specific correction of the parent coordinates applied, for example, in [10,16,67] was not used here, because for the proposed interpolation scheme of the assumed stress resultants and strains it gives the same results as the natural coordinates. The high efficiency of the proposed mixed elements is presented in the examples of twisted beam and pinched hemisphere. They allow for significantly larger load steps than the previous EANS4 and CAMe16 elements formulated in the nonlinear 6-parameter shell theory. They require also less equilibrium iterations than the mixed HW47 element. The convergence rate of the MIX40s and MIX40 elements is slightly better than the MIX44s and MIX44 elements. Taking into account theirs high accuracy and robustness, the shell elements MIX40, MIX40s are selected as the best performers from the developed mixed hybrid elements. Theirs advantage is also the lowest number of the independent parameters which is associated with the shortest time of a single iteration.  N β1 (δ β1 + ε β1 ) + Q β ε β N β1 (δ β2 + ε β2 ) −Q β (δ β2 + ε β2 ) N β2 (δ β1 + ε β1 ) N β2 (δ β2 + ε β2 ) + Q β ε β Q β (δ β1 + ε β1 ) −N β2 ε β N β1 ε β N αβ (δ αβ + ε αβ ) The matrix K (e) G is not entirely symmetric; therefore, the symmetrization procedure is applied to asymmetric terms. In contrast to [16,17], the ANS approach is not applied during calculation of the geometric matrix, because such method yields slightly worse results in the considered shell theory.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.