Improved nine-node shell element MITC9i with reduced distortion sensitivity

The 9-node quadrilateral shell element MITC9i is developed for the Reissner-Mindlin shell kinematics, the extended potential energy and Green strain. The following features of its formulation ensure an improved behavior: and transverse shear strains, which regular straight element sides middle positions of midside a To reduce shape distortion effects, so-called corrected shape functions of Celia and Gray 1984) are extended to shells and used instead of the standard ones. In all patch tests are passed additionally for shifts of the midside nodes along straight element sides and for arbitrary shifts of the central node. Several extensions the are proposed of non-ﬂat shells. In particular, a criterion is put forward to determine the shift parameters associated with the central node for non-ﬂat elements. Additionally, the method is presented to construct a parabolic side for a shifted midside node, which improves accuracy for symmetric curved edges. Drilling rotations are included by using the drilling Rotation Constraint equation, in a way consistent with the additive/multiplicative rotation update scheme for large rotations. We show that the corrected shape functions reduce the sensitivity of the solution to the regularization parameter γ of the penalty method for this constraint. The MITC9i shell element is subjected to a range of linear and non-linear tests to show passing the patch tests, the absence of locking, very good accuracy and insensitivity to node shifts. It favorably compares to several other tested 9-node elements.


Introduction
It has been clear since the earliest implementations of a basic 9-node element, that it has several limitations, i.e. it is excessively stiff (locks) and its accuracy is diminished by shape distortions.
To alleviate the problem of an excessive stiffness, five types of modifications to the standard formulation of the 9node element have been proposed: 1. Uniform Reduced Integration (URI) in [53], with the 2 × 2 instead of the 3 × 3 Gauss integration scheme. This yields a rank-deficient stiffness matrix and requires a stabilization, which was developed e.g., in [4,5,31] and for shells in [50]. 2. Selective Reduced Integration (SRI) in [32], with different integration schemes for various parts of the strain energy. This method yields the correct rank of the stiffness matrix but strain components are computed at different integration points, which limits the range of application of the element. 3. Two-level approximations of strains, which are developed as either the Assumed Strain (AS) method or the Mixed Interpolation of Tensorial Components (MITC) method, see, e.g., two books [11,17]. These were invented to overcome the above mentioned limitation of the SRI. (Note that in the 4-node element, this technique is applied to transverse shear strains and designated the Assumed Natural Strain (ANS) method.) We will show in the current paper that the MITC method still has the potential for improvement. 4. The enhanced element based on the Enhanced Assumed Strain (EAS) technique, which uses a set of 11 modes for membrane behavior [7]. They are incorporated using an identical formula to that used for 4-node elements, see [39,40,42,51]. We use these modes in our version of the 9-EAS11 shell element. 5. Partly hybrid formulations based on the Hellinger-Reissner functional, e.g., [34,35]. Hybrid stress, hybrid strain and enhanced strain 9-node shell elements are developed and assessed in [36]. We are not aware of any 9-node shell element based on the Hu-Washizu functional, which was a basis of several 4-node shell elements of very good accuracy and exceptional robustness in nonlinear applications, see e.g. [16,41,45,47,48].
The technique of two-level approximations of strains is applied either to Cartesian strain components (the AS method) or to covariant strain components (the MITC method); in the current paper we focus on the latter one. The strain components are sampled at the points where they are of good accuracy, and, next, extrapolated over the element. In effect, all strain components are available at each of the 3 × 3 Gauss integration points.
For sampling of the shell strain components (11,13) and (22,23), the 2 × 3 and 3 × 2-point schemes are used in the literature respectively. Several sets of sampling points and interpolation functions were proposed for the (12) strain components, i.e. the in-plane shear strain ε 12 and the twisting strain κ 12 , 1. The 2 × 3 and the 3 × 2-point schemes in [19,21,31]. Note that the last paper uses different approximation functions from the first two. 2. The 2 × 2-point scheme in [8,28]. This scheme uses exactly the same points for the shear component as the SRI of [32] (Table 1, p. 580).
In our tests, the 2 × 2-point scheme for the (12) strain components provided the correct rank of the tangent matrix and was the most accurate.
Another problem with the 9-node elements is the sensitivity of solutions to shifts of the midside nodes and the central node from the middle positions, which causes loss of accuracy for all the element formulations listed above. Some earlier studies, e.g., [26] suggest that only regular meshes, i.e. with straight element boundaries and evenly spaced nodes, give satisfactory results, while the accuracy is poor for curved boundaries or unevenly spaced nodes. Since then, however, research has been conducted to reduce this sensitivity.
The best remedy found so far, are the Corrected Shape Functions (CSF) of [9] used instead of the standard shape functions based on isoparametric transformations of [15]. In [9], the CSF are tested for an 8-node (serendipity) element for the Laplace equation (heat conduction) and the 4 × 4 integration rule. In [30], we examine several 9-node 2D elements for plane stress elasticity and the 3 × 3 Gauss integration: QUAD9** [19], MITC9 [2] and ours: 9-AS [28] and MITC9i [43]. Summarizing the performed tests, because of the CSF, the midside nodes can be shifted along straight element sides without rendering errors while for perpendicular shifts of these nodes errors are smaller and negative values of the Jacobian determinant are often avoided. Also the central node can be shifted in an arbitrary direction without causing errors. All of the tested 2D elements benefited from using the CSF, hence, in the current paper, we further develop this technique to make it applicable to shells.
Objectives and scope of the paper The objective of the current paper is to develop an improved 9-node MITC9i shell element with drilling rotations.
1. The original 9-node MITC9 element has good accuracy but does not pass the patch test. In [43], we have scrutinized its formulation to find the source of this problem. By an alternative but equivalent formalism, the formulation of the MITC9 element was seen from a different perspective, and the modified transformations were proposed. As a result of this change, the membrane MITC9i element passed the patch test for a regular mesh, i.e. with straight element sides and middle positions of the midside nodes and the central node.
In the current paper this technique is extended to the bending and transverse shear strains of a shell element MITC9i, for which proper transformations between Cartesian and covariant components are devised. In particular, we verify whether this change suffices to pass the patch tests for middle positions of the midside nodes and the central node. 2. The problem of sensitivity to node shifts is particularly important for 9-node elements as it affects a parametrization of the element domain and element accuracy. In the current paper, we improve the shell element robustness to node shifts extending the corrected shape functions of [9]. In [30], we tested several membrane 9-node elements integrated by the 3 × 3 Gauss rule, and these functions have proven beneficial to all of them. In the current paper, we further develop this technique to make it applicable to shells. For a 9-node shell element located in 3D space, we have to extend the method of computing the node shift parameters due to a presence of the third coordinate. Additionally, a new criterion is needed to determine the shift parameters associated with the central node for nonflat elements. Next, we have to verify whether the shell element MITC9i using the corrected shape functions passes patch tests for the midside nodes shifted along straight element sides and for the central node in an arbitrary position. Additionally, we noticed that the original method to determinate the shift parameters of [9] yields a non-symmetric side curve when the midside node is shifted. Therefore, we propose an alternative method, which does not rely on the proportion of arc-lengths, but uses a parametric equation of a parabola to construct a symmetric curve also when the midside node is shifted. We describe this method and evaluate its accuracy for an example with curved and symmetric shell edges. 3. The drilling rotation is incorporated into the formulation of the 9-node shell element via the drilling Rotation Constraint (RC) and the penalty method. We implemented it consistently with the additive/multiplicative rotation update scheme for large rotations, which combines a 3parameter canonical rotation vector and a quaternion.
Regarding the drilling RC, we analyze its form for equalorder bi-quadratic interpolations of displacements and the drilling rotation to check whether it contains a faulty term, analogous to the ξη-term in 4-node elements, see [48], and to evaluate how it affects the solution. Besides, we verify how the corrected shape functions change the sensitivity of a solution to the regularization parameter γ of the penalty method, in particular in the case of curved or distorted elements.
Finally, the 9-node MITC9i shell element with drilling rotations is tested on a range of linear and non-linear numerical examples, which are performed to check passing the patch tests, an absence of locking, accuracy, an insensitivity to node shifts, and correctness of implementation of the drilling rotation; a selection of these tests is presented in Sect. 6. Performance of the MITC9i element is mostly compared to other 9-node elements, some of them with similar improvements implemented as in the tested element. Besides, reference results obtained by our 4-node shell elements are provided; the corresponding ones for triangular elements can be found e.g., in [10].

Shell element characteristic
Reissner-Mindlin shell kinematics The position vector of an arbitrary point of a shell in the initial configuration is expressed as where X 0 is a position of the reference surface and t 3 is the shell director, a unit vector normal to the reference surface. Besides, ζ ∈ [−h/2, +h/2] is the coordinate in the direction normal to the reference surface, where h denotes the initial shell thickness. In a deformed configuration, the position vector is expressed by the Reissner-Mindlin kinematical assumption, where x 0 is a position of the reference surface and Q 0 ∈ SO(3) is a rotation tensor.
= ∂X/∂ξ . The right Cauchy-Green deformation tensor is and the Green strain is where C 0 . = C| x=X = I. The Green strain can be expressed as a series of the ζ coordinate, where the higher order terms are neglected. The 0th order strain E 0 includes the membrane components ε and the transverse shear components γ /2 while the 1st order strain E 1 includes the bending/twisting components κ. We assume that the transverse shear part of E 1 is negligible, i.e. κ α3 ≈ 0 (α = 1, 2). Note that the normal strains ε 33 and κ 33 are equal to zero because of Eq. (2) and must be either recovered from an auxiliary condition, such as the plane stress condition or the incompressibility condition, or be introduced by the EAS method. In the current paper, we use the plane stress condition for this purpose.
Incremental form of rotation Q 0 To enable computation of large rotations of a shell, we use an incremental approach and an update scheme involving a rotation vector ψ and a quaternion q n , where n is the increment (step) index. Consistently with this scheme, the rotation Q 0 ∈ SO(3) of Eq. (2) is assumed in the following form where Q 0 ( ψ) which is the Rodrigues's formula for the canonical parametrization of a rotation tensor. Besides, = ψ × I is the skew-symmetric tensor and I is the second-rank identity tensor. Within an increment, we use the Newton method and additively update the rotation vector after each iteration, ψ i+1 = ψ i +δψ, where δψ is the rotation vector increment for an iteration and i is the index of iterations. When the Newton method have converged, ψ i+1 is converted to a quaternion for an increment q, and the total quaternion is multiplicatively updated q n+1 = q • q n . This scheme works very well for statical computations in the range of large rotations.

Rotations in mechanics
In mechanics, the rotations can be treated in two ways: either as independent variables, e.g. in Cosserat continuum, or as dependent variables, e.g. in Cauchy continuum; the latter case is of interest in the current paper. In the Cosserat shells, see e.g. [12,34] and [35], the drilling rotation is naturally present but the constitutive equations are complicated.
For Cauchy continuum and the known deformation gradient F, the rotations can be obtained by the polar decomposition F = RU = VR, where the rotation R ∈ SO(3) and U, V is the right and left stretching tensor, respectively. But this is merely a post-processing which is not a subject of the current paper; our goal is to include the rotations into the formulation.
To include rotations as additional variables into 3D equations, first, we define the extended configuration space where C . = {χ : B → R 3 } is the classical configuration space, and, next, we impose the Rotation Constraint (RC) to obtain the Cauchy continuum, see [37]. Here Q ∈ SO(3) is the rotation, which at a solution of the extended system of equations (equilibrium equations plus the RC) is equal to R yielded by the polar decomposition of F; this issue is considered in detail in [44].
For shells, we can use an analogous approach to that for the 3D Cauchy continuum, but with the following alterations: (a) the rotation tensor is constant in ζ , i.e. Q ≈ Q 0 , where Q 0 is defined in Eq. (7), and (b) only the drilling component of ψ is considered because the tangent ones are introduced by the term Q 0 t 3 in Eq. (2).
Drilling rotation constraint The drilling rotation for an increment is defined as an elementary rotation about the forward-rotated director a n 3 , i.e. ω . = ψ · a n 3 , where a n 3 = Q n 0 t 3 , see Fig. 1. Note that Q 0 t 3 = Q 0 a n 3 in Eq. (2), so the tangent components of ψ in the {a n k } basis are present in shell equations while the normal one, i.e. the drilling rotation, must be included in a different way.
Let us define the drilling Rotation Constraint (drilling RC) as the (1,2) component of the RC of Eq. (10) in the local basis {t k } (k = 1, 2, 3), It can be written in the incremental form consistent with the form of Q 0 of Eq. (7) and the multiplicative decompositions where A . = Q T 0 F, the convected vectors are t n α . = F n t α and the forward-rotated vectors are a n α . = Q n 0 t α (α = 1, 2); a derivation of the above formula is given in Sect. 5. To obtain the shell equations with drilling rotations, we construct an extended potential energy functional including the above scalar equation.
Physical interpretation of drilling rotation For simplicity, consider the planar (2D) deformation, so only the drilling rotation ω matters in the definition of Q 0 while the tangent rotation components are equal to zero. To obtain an interpretation of ω, we consider a pair of orthogonal unit vectors t 1 and t 2 , associated with the initial configuration, which are rotated and stretched by F, so we can write where λ 1 , λ 2 > 0 are stretches and β 1 , β 2 are rotation angles. Using these formulas in the drilling RC in the standard (non-incremental) form (see the derivation in Eq. (59)), we obtain for cos ω = 0, λ 1 c 1 + λ 2 c 2 = 0 and λ α ≈ 1, where c 1 .
= cos β 2 . Hence, the drilling angle ω is an average of rotations of vectors t 1 and t 2 ; for details see [46] "Appendix".
Extended potential energy functional for shells with drilling rotation To obtain the shell equations with drilling rotations, we define the extended shell potential energy, This functional includes the shell strain energy, the work of external forces F sh ext , and an additional term for the drilling RC, F sh drill . Besides μ .
= det Z, where Z is the shifter tensor, and A is the shell reference surface. We discuss below the first and third term. (A) The strain energy W sh (E) = W sh (E 0 , E 1 ) by using the shell form of Green strain E(ζ ) ≈ E 0 +ζ E 1 , see Eq. (6). For instance, for a linear material and symmetry of its properties w.r.t. the reference surface, upon integration over thickness, we obtain the well-known: The drilling rotation term F sh drill is assumed in the penalty form, where c is defined by Eq. (11) and γ ∈ (0, ∞) is the regularization parameter. The upper bound on γ equal to the shear modulus G was found for an isotropic elastic material in [20], but a reduced value can be beneficial for non-planar shell elements; in the current paper we test also γ = G/1000. Note the thickness h in this formula. Another question is how to select γ for non-isotropic and/or inelastic materials for which theoretical results are not available. For an orthotropic elastic material the shear modulus for the in-plane shear G 12 should be used instead of G. For multilayer shells composed of orthotropic elastic layers of different orientation, the effective material is anisotropic. We treat this problem approximately, expressing G 12 of a substitute orthotropic material in terms of the effective inplane stiffness matrix D 0 . The elasto-plastic tangent matrix C ep can be used similarly.

Remark on PL method
We have also tested the Perturbed Lagrange (PL) form of the drilling rotation term, where T * is the Lagrange multiplier. The second term under the integral provides a regularization in T * , and a small perturbation of the tangent matrix which is needed when it is singular, i.e. for the drilling rotation ω = 0. Note that for 9node elements, we have to use 9 parameters for T * to avoid singularity of the tangent matrix. They can be either values at nodes or coefficients of interpolation functions, and can be eliminated (condensed out) on the element level. For 4node mixed/enhanced shell elements, the PL form implies a reduced sensitivity to distortions and a larger radius of convergence in non-linear problems, see [48,49]. For 9-node elements, the PL form also works well but is slightly less beneficial, hence the penalty form of Eq. (18), for which the element is faster, is tested in Sect. 6.
Taking a variation of the governing functional of Eq. (16) and performing a consistent linearization of it, we obtain the tangent stiffness matrix K for the Newton method; this procedure is standard and, therefore, not described here. The standard procedure is also used to compute the stress and couple resultants.
Basic definitions for 9-node shell element We consider a 9node isoparametric element with the nodes numbered and named as in Fig. 2, where 1, 2, 3, 4 are the corner nodes, 5, 6, 7, 8 are the midside nodes, and 9 is the central node. For a shell element, we use a Cartesian basis associated with it as a reference basis.
Consider the following vectors associated with the reference surface of a shell: the initial position X, the current The above vectors are interpolated as follows: where the natural coordinates ξ, η ∈ [−1, +1] and I is the node number. Note that the rotation vectors ψ I are used in a step, and converted to quaternions at the end of the step. The shape functions N I (ξ, η) are presented in Sect. 4. The tangent vectors of the natural basis on the reference surface are defined as The vectors g k of the co-basis {g 1 , g 2 } are defined in the standard way by g α · g β = δ α β , α, β = 1, 2. In general, g 1 and g 2 as well as g 1 and g 2 are neither unit nor mutually orthogonal. The director t 3 is defined as a unit vector perpendicular to vectors g 1 and g 2 , Tangent vectors of the Cartesian basis {i k } can be defined as where the auxiliary vectors arẽ The normalized natural vectors are denoted by a tilde, i.e. g k = g k / g k . Note that i 1 and i 2 are equally distant from g 1 and g 2 , see Fig. 2 in [48]. The above vectors for the element center are designated by the letter "c".

MITCi method for bending/twisting and transverse shear strains
Standard and improved 2D MITC9 element A lot of research has been devoted to various aspects of the Mixed Interpolation of Tensorial Components (MITC) method; for a comprehensive review see [11], [8]. Here we focus on the transformations used in this method. The MITC9 is a 9-node element formulated using the MITC method and fully integrated (3 × 3 Gauss rule). This element uses the covariant components of the Green strain, which e.g. for 2D are as follows For X and x interpolated as in Eq. (20) and g α of Eq. (21), the above components can be directly computed. Then, in the MITC method, these strains are sampled at selected points, interpolated, and transformed to the local Cartesian basis at each Gauss point. Generally, the MITC9 element performs quite well but it does not pass the patch test even for a regular mesh shown in Fig. 7, i.e. with straight elements' sides and middle positions of the midside nodes and the central node. This motivated our work on the formulation of this element in [43] and [30], and below we shortly describe the improvement which we proposed.
First, we note that the matrix ε ξ of covariant components of Eq. (25) can be obtained from the Cartesian components ε re f by the transformation where j .
. We noticed that j varies over the element, and, in consequence, the covariant ε ξ is associated with a different co-basis at each sampling point. We put forward the hypothesis that this was the cause of failing the patch test by the classical MITC9 element, and proposed an improvement described below.
Because the Jacobian j varies over the element, it can be decomposed j . = j c j, where j designates the relative Jacobian between the element center and a Gauss point. Then, Hence, we proposed to neglect j and use the Jacobian at the element center j c instead of j. Note that the components ε c ξ yielded by such a transformation are not exactly the covariant components of Eq. (25); for clarity, we designate them by "COVc". Remark The Jacobian at center was used in [40] to improve the method of incompatible modes for a 4-node element.
However, no such an improvement was proposed for the MITC method until our paper [43], where the MITC9i element was developed. The main difficulty was to notice that the covariant components of Eq. (25) can be obtained from the Cartesian ones by Eq. (26), and apply this formula in the context of sampling/interpolation.

MITCi method for bending/twisting and transverse shear strains
For the membrane strains ε of the shell element MITC9i, we apply the transformations which we developed in [43]. In the present paper, we propose and test analogous transformations for the bending/twisting strains κ and the transverse shear strains γ . Note that the transformation formulas between Cartesian and covariant components are different for the in-plane strains, i.e. ε and κ, and for the transverse shear strain γ , see "Appendix".
Steps defining the improved shell element MITC9i are as follows: 1. The representations in the reference Cartesian basis are transformed to the co-basis at the element center, to obtain the COVc components, 2. The two-level approximations of the COVc components are performed, This involves sampling and interpolation, which are described in detail in the next paragraph. Results are designated by the tilde " .
3. The approximated COVc components in the co-basis at the element center are transformed back from the co-basis at the element center to the reference Cartesian basis, The transformations of the first and third steps are reciprocal and without the second step we would simply obtain κ re f = κ re f and γ re f = γ re f .

Sampling points and interpolation functions
The two-level approximations are applied to the COVc components of shell strains, as symbolically given by Eq. (29). The COVc strain components are sampled at points shown in Fig. 3, where a = √ 1/3 and b = √ 3/5. Let us group the sampled shell strain components, Then the particular strain groups are interpolated, using the sampled values, as follows: where R l are the interpolation functions defined below, and l is the index of sampling points: l = A, B, C, D, E, F for the 6-point schemes and l = A, B, C, D for the 4-point scheme. The sets of interpolation functions are defined as follows: 1. for ε 11 , κ 11 , γ 31 , the points of Fig. 3a are used (2 points in the ξ direction), 2. for ε 22 , κ 22 , γ 32 , the points of Fig. 3b are used (2 points in the η direction), 3. For ε 12 , κ 12 , the points of Fig. 3c are used (2 points in ξ and η directions), Remark The sampling points of Fig. 3 correspond to the integration points proposed for the Selective Reduced Integration (SRI) of strain energy terms in [32]. Several attempts were made to apply the 6-point schemes of Fig. 3a and b to ε 12 and κ 12 but without a significant improvement, see [19,21,31].
Improved sampling strategy Note that the schemes of Fig. 3a and b involve 6 sampling points each, which implies a con-siderable number of evaluations. To reduce it, we use the so-called sampling lines instead of points, which yields a more efficient implementation. This method is not fully equivalent to the 6-point scheme due to the transformations of Eqs. (28) and (30) but works very well.
To explain the method, we consider only one strain component E c ξξ , for which sampling points and integration points are shown in Fig. 4a. Let us, for a moment, neglect -sampling points -Gauss points the effect of transformations of Eqs. (28) and (30). Then we see that both these types of points are located at the same η ∈ {−b, 0, +b}, where b = √ 3/5. Hence, the 3×3 integration does evaluate E c ξξ at correct η-coordinates, and, therefore, no separate sampling and interpolation in the η-direction is needed. In consequence, we sample and interpolate E c ξξ only in the ξ -direction, where This formula involves two sampling lines shown in Fig. 4b, where L1 is located at ξ = −a and L2 at ξ = a. Besides, note that the interpolation functions R L1 and R L2 of Eq. (37) replace six much more complicated functions R l of Eq. (33). The sampling lines are transformed back by Eq. (30).

Corrected shape functions for shell element
In this section, we present several modifications to the method of calculation the shift parameters for the corrected shape functions, which: (a) are necessary for 9-node shell elements located in 3D space, and (b) improve accuracy of a solution for non-flat shell elements and (c) enable to construct a symmetric side curve for a shifted midside node. The last modification is also applicable to 2D problems. The standard isoparametric shape functions are obtained by the assumption that the midside nodes (5,6,7,8) are located at the middle positions between the respective corner nodes and the central node 9 is located at the element center. When these nodes are shifted from the middle positions then the physical space parameterized by the standard shape functions is distorted, see e.g. Figs. 13a and 20 in [30].
To alleviate this problem, the corrected shape functions (CSF) were proposed in [9], with six additional parameters introduced in the local coordinates space, α, β, γ , , θ, κ ∈ [−1, +1], see Fig. 5. The distance of these nodes from the middle positions in the local coordinates space is computed as proportional to their distance in the physical space. To determine these 6 parameters, nonlinear equations must be solved: 4 equations with 1 unknown each and 2 equations with 2 unknowns. This is done only once, so the time overhead is insignificant.
The corrected shape functions for the 9-node element are defined in two steps. First, the shape functions of the 8-node (serendipity) element are defined, , , , and they account for shifts of the midside nodes from the middle positions. Next, the basis function for the central node 9 is added hierarchically to the shape functions for the 8-node element. The obtained shape functions for the 9-node element are [9], Eq. (20). These functions account for shifts of the midside nodes from the middle positions and the central node from the central position. When the shift parameters are equal to zero, then the CSF of Eq. (39) are reduced to the standard shape functions.

Computation of shift parameters in 3D for midside nodes
For the midside nodes (5,6,7,8), the procedure to determine the shift parameters α, β, γ , can be extended from 2D to 3D space as follows. For instance, let us determine the shift parameter α of node 5 in local coordinates. Consider the parametric form of the side 1-5-2, which, in general, can be curved. For 2D problems, it is where the vector of corrected shape functions for a 3-node element is For shells in 3D space, we have to write additionally a similar expression for the third coordinate, In the method proposed in [9], the shift parameter α is obtained from the proportion of arc-lengths, i.e. the fractional distance of node 5 along the curved side relative to nodes 1 and 2 is required to be identical in a physical space and in a local space, where the arc-length of the side in a physical space from In the above form, we included the term (d Z/dξ ) 2 , i.e. the contribution of the third coordinate, which is needed for shells. The definition of L +1 −1 is analogous; we have to replace α by +1. Besides, for the right-hand-side, we obtain where a, b, c depend on differences in node positions and on α. We solve this equation using the Newton method, and the term added for shells, (d Z/dξ ) 2 , contributes to the tangent matrix and the residual vector. An initial value of α, denoted as α 0 , is obtained by a proportion and a change of variables as follows: where L i j is the distance between nodes i and j in 3D space, computed as a length of a vector connecting these points. When the side 1-5-2 is straight then, L 15 + L 52 = L 12 and the so-defined α 0 is a solution of the nonlinear Eq. (43). The parameters for the other midside nodes, i.e. β, γ and , are obtained in an analogous way.

Computation of shift parameters in 3D for central node
For the central node 9 and the (θ, κ) parameters, the generalization to 3D space is more complicated. For 2D space, we use an assumption that the Jacobian of transformation from physical to local coordinates should not be affected by the central node 9. For shells, we have additionally the third coordinate Z , and we can consider two methods of treating it: M1. Neglect a contribution of the Z -coordinate, and solve a system of 2 equations identical to that used for 2D space, This method is an option when the elemental reference basis is used, e.g. the tangent basis attached at node 9, as then Z describes elevation or warping of an element. M2. Account for a contribution of the Z -coordinate. We propose to incorporate it using a definition of the square error where r(θ, κ) . = ⎡ ⎢ ⎣ The error e is a scalar, and minimization of it w.r.t. θ and κ yields 2 equations This method can be applied when either the elemental or global reference basis is used for an element. For a flat shell element, M2 yields exactly the same results as M1, so the proposed generalization is correct indeed.
In both methods, the set of nonlinear equations is solved for θ and κ using the Newton method. Initial values can be calculated, e.g., as follows: where L i j is the distance between nodes i and j in 3D space. When nodes 8, 9 and 6 lie on one straight line then θ = θ 0 . Similarly, when nodes 5, 9 and 7 lie on one straight line then κ = κ 0 . The proposed method M2 is tested in Sect. 6.4.

Alternative computation of shift parameters for curved sides
When the element side is straight, then the shift α = α 0 , where α 0 is given by Eq. (46). Computations are more complicated when the side is curved; e.g. in the method of [9], the proportion of arc-lengths yields a nonlinear Eq. (45) from which α is calculated. We have noticed two features of that procedure: 1. when the midside node is shifted from the middle position then the so-computed α yields a non-symmetrical side curve, 2. any particular shape of the side curve is not assumed in this method so we cannot control it. It is a drawback, e.g., for domain boundaries, the shape of which we try to represent as exactly as possible.
We propose an alternative method by which a parabola is constructed in such a way that a symmetric shape of the side is obtained even for a shifted midside node. The steps of the proposed method for the given nodes P 1 , P 2 and P 5 are as follows: 1. Check the position of P 5 relative to P 1 and P 2 where p 1 , p 2 and p 5 are position vectors of P 1 , P 2 and P 5 , respectively, and tol is a prescribed tolerance. If d < tol then P 5 is on the axis of symmetry so α = 0 and this ends the calculations; otherwise we proceed further. 2. Construct vectors of a local Cartesian basis in 3D using P 1 , P 2 and P 5 Note that s t and s n belong to the plane defined by P 1 , P 2 and P 5 . 3. Calculate coordinates of P 1 , P 2 and P 5 in the local basis {s t , s n , s r }, where p 0 = (p 1 +p 2 )/2 and O . = [s t | s n | s r ] ∈ SO(3). The applied transformation involves a shift by p 0 and a rotation by O. Associate the local coordinates x, y with the vectors s t and s n , respectively, so let them have zero at p 0 . The local coordinates of nodes are designated: r 1 = [x 1 , y 1 ], r 2 = [x 2 , y 2 ] and r 5 = [x 5 , y 5 ], where y 1 = y 2 = 0. 4. Calculate coefficients of a parabola in the local basis, assuming that it passes through P 1 and P 2 , and is symmetrical w.r.t. x = 0. Writing y = ax 2 + c for P 1 and P 5 , we obtain a system of two equations, from which we can calculate 5. The parametric equation of parabola y = ax 2 + c is where t is a parameter. The value of t for P 1 can be obtained using the condition x = x 1 , and the value of t for P 5 using the condition x = x 5 , Then, the shift parameter α ∈ [−1, 1] can be defined as a ratio The new method is simpler than that of [9] and its accuracy is better when the shell edge is curved and symmetric indeed, see Sect. 6.5. Note that the midside shift parameters, such as α, affect also the (θ, κ) parameters for the central node, see Eq. (50). Finally, we note that e.g. a semi-circle can be implemented in a similar manner, but a parabola is simpler.
Example Consider three points P 1 (0,0), P 2 (4,0) and P 5 (1,2) shown in Fig. 6, where the side curves are obtained for the standard shape functions ("standard") and the corrected shape functions. For the latter, α is obtained in two ways: the method of [9] yields α = − 0.264405 and the nonsymmetric curve "CSF,CG", while the proposed method yields α = − 0.5 and the symmetric curve "CSF,new".

Treatment of drilling RC
The purpose of using the drilling Rotation Constraint (drilling RC) is to provide an additional equation for the = ψ · a n 3 , which is not present in the shell kinematical assumption (2) and in strains.
We define the drilling RC as the (1,2) where a 1 . = Q 0 t 1 and a 2 . = Q 0 t 2 are the forward-rotated tangent vectors t 1 and t 2 .
The drilling RC can also be obtained in the forward-rotated basis {a k }, where a k . = Q 0 t k , which is associated with the current (deformed) configuration. Let us denote A . = Q T 0 F and its forward-rotated form A * .
= Q 0 AQ T 0 = FQ T 0 . We see that where the last form uses the forward-rotated basis vectors. When we approach the solution, i.e. Q 0 → R, then A .
where the rotation R ∈ SO(3) and the right and left stretching tensors, U and V, respectively, are obtained by the polar decomposition of F.
Incremental form of the drilling RC An incremental form of the drilling RC can be obtained from the last form of Eq. (59) using the multiplicative decompositions F = FF n and Q 0 = Q 0 Q n 0 , where the superscript "n" designates the last converged configuration. Introducing the convected vectors t n α . = F n t α and the forward-rotated vectors a n α . = Q n 0 t α (α = 1, 2), after straightforward transformations, we obtain c = 1 2 a n 1 · At n 2 − t n 1 · A T a n 2 , where A . = Q T 0 F. This form of c complies with the additive/multiplicative rotation update scheme, in which, we use Q 0 parametrized by ψ for the increment, see Eq. (8). This scheme works correctly for | ω| < π/2, and, when combined with a quaternion, allows for arbitrarily large drilling rotations. For the initial configuration, F n = I and Q n 0 = I, so we have t n α = t α , a n α = t α , and the above formula reduces to where a α . = Q 0 t α , of an apparent similarity to the last form of Eq. (59).
Note that the last ξ 2 η 2 -term contains only the rotational coefficient 9 but no displacement coefficients U i , V i , which is incorrect, because both these types of coefficients should be linked in each term of the drilling RC. This problem is caused by equal-order approximations of displacements and the drilling rotation; a similar one appears for 4-node (bilinear) elements, but then the ξη-term is faulty, see [48].
Let us designate the faulty term as c 9 . = 9 ξ 2 η 2 . For a bi-unit square element and the linearized drilling RC, we can express c 9 in terms of nodal drilling rotations ω k (k = 1, . . . , 9) as follows: Note that for the rigid element rotation ω rr , all nodal ω k = ω rr and the term in brackets vanishes, so the faulty term does not cause any problem then. However, in general, = 0, which means that the penalty form involving c 9 is non-zero and affects the solution.
To establish a magnitude of the effect caused by c 9 , we: (1) removed c 9 from c, which implied 1 additional zero eigenvalue of the tangent matrix, and (2) stabilized this matrix by the scaled-down c 9 . Then, Eq. (18) valid for an arbitrary element (not only for the bi-unit one), has the following modified form where 1/10 3 is the scaling factor. This form was tested on the Cook's membrane example of Sect. 6.3, and the effect of c 9 term was not strong. Hence, in further tests, we use the form of F sh drill given by Eq. (18), which is simpler.

Numerical tests
In this section, we present numerical tests of the developed 9-node shell element MITC9i. The element uses the modified transformations of Sect. 3 and the corrected shape functions (CSF) in the extended version for shells of Sect. 4, which enables calculation of shift parameters for non-flat element geometry and generation of a symmetric side curve for a nonsymmetric position of a midside node. An implementation of the drilling rotation is described in Sects. 2 and 5. The element is integrated using the 3 × 3 Gauss rule. The tested and reference shell elements are listed in Table 1. All "ours" shell elements are of the Reissner-Mindlin type and with drilling rotations (6 dofs/node).
The elements were derived using the automatic differentiation program AceGen described in [23,24], and were tested within the finite element program FEAP developed by R.L. Taylor [52]. The use of these programs is gratefully acknowledged. Our parallel multithreaded (OMP) version of FEAP is described in [22].
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 a tangent matrix are computed for a single unsupported element. First, a single bi-unit square element with regularly placed mid-side and central nodes is checked. Then, several other flat and non-flat (h-p, one-and twocurvature) element shapes with shifted mid-side and central nodes are examined. This test is performed for the standard as well as corrected shape functions, which implies different element shapes. For all these cases, the tested shell element MITC9i has a correct number of zero eigenvalues (6).

Patch tests
The patch tests are crucial to demonstrate the benefits of using the new MITC9i shell element. The five-element patch of elements proposed in [33] of Fig. 7 and its distorted form of Fig. 8 are used. The data E = 10 6 , ν = 0.25 and a thickness h = 0.001. The stretching/shearing (membrane) and bending/twisting patch tests are performed as described in [27]; displacements and rotations are prescribed at external nodes and the computed values are checked at internal nodes. The transverse shear test is performed for the load case defined for a 9-node plate in [18], see "Shearing case" in Fig. 2b therein. The particular types of constant strains are verified as follows: 1. For the stretching/shearing (membrane) patch test, the displacement and rotation fields are This test verifies that ε x x = ε yy = 0.001 and ε xy = 0.0005 while the bending/twisting strains and the transverse shear strains are equal to zero.
The drilling rotation, ψ z is unconstrained, and we should obtain the solution ψ z = 0, where zero results from the linearized drilling RC, ψ z = ω . = 1 2 (u x,y − u y,x ).
2. For the bending/twisting patch test, the displacement and rotation fields are This test verifies that the bending strains κ x x = κ yy = 0.001 and the twisting strain κ xy = −0.0005 while the membrane strains and the transverse shear strains are equal to zero. Note that passing this test is not required for consistency; it is merely a test of ability to model constant bending/twisting strains. 3. For transverse shear patch, the uniform transverse force P z = 1 is applied at the right boundary, the left boundary is clamped. Rotations ψ x and ψ y are set to zero at all nodes of the mesh to prevent the bending moments developing. The analytical solution is and at the right boundary u z = 0.006. This test verifies that the transverse shear strains are: ε xz = 0.0125 and ε yz = 0. The standard patch tests are performed using the regular mesh of Fig. 7, which has straight element sides, the side nodes are placed at the middle positions and the central node of each element is located at the intersection of lines connecting opposite midside nodes.
We additionally perform the patch tests using the distorted meshes obtained from the regular mesh by shifts of the selected five nodes indicated in Fig. 8. The shifts are kept within the radius r = 0.00673145, which is a half of the minimal distance of nodes in the entire mesh. Pseudorandom real numbers are generated 100 times in the range [−r, r ], and the relative error is computed over all nodes, where u is a selected component and the analytical solution u ana is computed using either Eq. (67) or (68) or (69). The rounded logarithms of these errors are shown in Table 2, and we see that for the regular mesh (Case A), the errors are of order 10 −16 , which means that the MITC9i element passes all the standard patch tests. For the distorted meshes, additionally the corrected shape functions are activated, and the results can be summarized as follows: 1. For the shifted central node 25 (Case B) and for parallel shifts of nodes 21, 22, 23, 24 (Case C), the errors are of order 10 −5 , so the element passes all the patch tests for these cases. Note that to obtain such a performance the standard shape functions are not sufficient and the corrected shape functions are required. 2. For perpendicular shifts of nodes 21, 22, 23 and 24 (Case D), we observe a drop of accuracy in all three patch tests.
The biggest error occurs for bending/twisting, and this patch test is not passed. We checked bending and twisting separately, and they both contribute to the error. (To shed additional light on this case, we replaced the transverse shear part of the 9-node element by an assembly of the transverse shear parts of four 4-node elements, in which the ANS method of [3] was applied to eliminate locking. Then the error dropped to log 10 e u = −3.1.) Because the bending/twisting patch test is not necessary for consistency, the convergence of the FE solution is possible despite failing this test.
Finally, note that the classical MITC9 element does not pass these patch tests even for the regular mesh (Case A), so the MITC9i element provides an improvement.

Cook's membrane
In this test, we demonstrate that the use of corrected shape functions reduces sensitivity of a solution to the value of regularization parameter γ in the drilling RC term. The implementation of F sh drill of Eq. (18) is tested. The membrane is clamped at one end (all degrees of freedom are restrained, including drilling rotations), while at the other end, the uniformly distributed shear load P = 1 is applied, see Fig. 9. The data is as follows: E = 1, ν = 1/3 and a thickness h = 1.
Two meshes are used, the skew one of Fig. 9a, and the irregular one, where the central node 9 and three midside nodes (5,7,8) are shifted by d, see Fig. 9b. We use two mesh densities: 1 × 1 and 16 × 16 elements.
We test the MITC9i element using: (a) two values of the regularization parameter: γ = G and γ = G/1000, and (b) two types of shape functions, standard and corrected. The results are given in Table 3, and we summarize them as follows: 1. Solutions for the coarse 1 × 1-element mesh are more sensitive to the value of γ than those for the dense 16 × 16-element mesh. For both tested types of shape functions, the reduced value γ = G/1000 improves the accuracy. 2. For the irregular mesh of Fig. 9b and both mesh densities, the use of the corrected shape functions is beneficial and the solution is almost insensitive to the node shifts. The MITC9i element is more accurate than the other elements.
We see that for the corrected shape functions, the solution is less sensitive to the value of regularization parameter γ , which is another positive effect of using these functions.

Shift parameters for parabolic cylindrical element
In this test, we evaluate the method M2 of calculation the shift parameters θ and κ of the corrected shape functions for shells, which was proposed in Sect. 4. The parabolic cylindrical element shown in Fig. 10a is obtained from a bi-unit square X, Y ∈ [−1, +1] by shifting nodes 5, 7 and 9 by the vector [d x , 0, d z ]. (In the notation of Sect. 4, Z 5 = Z 7 = Z 9 = d z .) The remaining nodes are at Z = 0 and the standard X, Y positions.
For this element, we first calculate the shift parameters α, β, γ , and next θ and κ using methods M1 and M2 of Sect. 4. For the initial values given by Eq. (51), the Newton method always converged in less than 5 iterations. The values of θ obtained for the horizontal shift d x ∈ [0, 0.5] and the elevation d z ∈ [0, 0.7] are plotted in Fig. 10b, and we see that: 1. M1 yields always the value θ = d x , which is correct for flat elements and acceptable for shallow elements, when 0 ≤ d z ≤ 0.1, but is incorrect beyond this range. 2. For M2, the value of θ varies with d z and the bigger d z the θ differs more from d x , which is correct as the arc-length changes then. Note that, e.g., for d z = 0.7 and d x = 0.5, the error of M1 relative to M2 is about 17.7%. Concluding, the method M2 is appropriate for flat, shallow and non-flat elements, and can be used as a general method of calculation θ and κ for 9-node shell elements.

Single curved element
In this test, we evaluate the method of computation of the shift parameters for curved sides proposed in Sect. 4, which constructs a symmetrical side curve when the midside node is shifted from the middle position. The curved cantilever is clamped at one end and at the other end is loaded by either "Load 1" causing in-plane bending by a pair of in-plane forces, or "Load 2" causing twisting by a pair of transverse forces, see Fig. 11. The material data and geometry are: E = 1.0 × 10 5 , ν = 0.0, h = 0.025, inner radius = 0.095, outer radius = 0.105, force P = 0.1, as in [25]. One 9-node element is used and it's nodes are shifted as follows: (1) the midside nodes 6 and 8 by d = n * 0.001 along straight boundaries, (2) the midside nodes 5 and 7 by φ = n * 1 • in the circumferential direction along curved boundaries, and (2) the central node 9 in the radial direction by d = n * 0.001, see Fig. 11. The element shape for n = 3 in shown in Fig. 12, and we see that the standard shape functions very poorly approximate the geometry while the corrected ones much better.
The displacements at node 1 obtained for the increasing distortion parameter n are shown in Figs. 13 and 14. The curves are designated as follows: "standard" for the standard shape functions, while for the corrected shape functions: "CSF,CG" for the method of calculation of the shift parameter of [9], and "CSF,new" for the method of calculation of the shift parameter proposed in Sect. 4.
Next, we can calculate errors of the displacement at node 1 for the distortion parameter n = 4 relative to the displacement for n = 0, see Table 4. We see that the corrected shape functions are beneficial for both tested elements and significantly reduce errors. The proposed method of calculation of the shift parameter "CSF,new" is more accurate than the original method of [9] "CSF,CG". It renders that the element MITC9i is practically insensitive to nodes' shifts for both load cases.

Curved cantilever
In this test, we assess the effects of elements' curvature, skew shape and a varying shell thickness h on the accuracy of a solution.
The curved cantilever is fixed at one end and loaded by a moment M z at the other, see Fig. 15. The data is as follows: E = 2 × 10 5 , ν = 0, width b = 0.025 and radius of curvature R = 0.1. The FE mesh consists of 6 nine-node elements, which have either rectangular or skew shape, see [25]. The analytical solution for a curved beam is where I is the moment of inertia. Note that for the rectangular mesh, coefficients of the corrected shape functions are equal to zero while for the skew mesh non-zero values are obtained; the values for the 'middle' elements (θ and κ are obtained by M2 of Sect. 4) are given in Table 5.
The displacement u y at point A obtained by the linear analysis are shown in Fig. 16. The plots of the rotation ψ z are similar and for this reason not shown here. We can conclude this test as follows: 1. For the rectangular mesh, the solutions for all elements are are close to the analytical value, and we provide for reference the curve designated "MITC9i,rectang". The only exception is the element '9', which severely locks; then the curve is similar to the one shown for this element and the skew mesh. 2. For the skew mesh, a. all the tested elements, except '9', yield accurate solution for R/ h ≤ 100. Beyond this range, the accuracy of all elements worsens. The most accurate is S9R5, second is MITC9i but the difference between it and the next 9-EAS11, 9-AS and 9-SRI is not big. The solution for MITC9 was obtained only for log(R/ h) = 1, 2, see [29] p. 93. b. MITC9i behaves slightly better than the elements 9-AS and 9-SRI of [29]; this can be attributed to the corrected shape functions implemented in MITC9i. MITC9i performs also slightly better than 9-EAS11 when both use the corrected shape functions.

Pinched hemispherical shell with hole
In this test, the shell undergoes an almost in-extensional deformation and membrane locking can strongly manifest itself, see [4]. The elements are curved and we test two values of the regularization parameter γ for the drilling rotation term. Linear analyses Two values of the regularization parameter γ for the drilling RC term are tested, γ = G and γ = G/1000, and the results for the 8 × 8 element mesh are given in Tables 6 and 7. A displacement at the point where a force is applied and in a direction of this force is reported. The conclusions are as follows: We see that the solution for the basic element 9 is very locked, and such techniques as EAS, MITC, AS or SRI improve accuracy about threefold for h = 0.04 and over 30 times for h = 0.01.

Non-linear analyses
The non-linear analyses are performed using the Newton method and P = 0.2. We use the mesh of 8 × 8 elements, unless differently indicated. The solution curves for the inward displacement under the force are shown in Fig. 18.
1. The curves MITC9i, 9-SRI0 and 9-EAS11 are obtained for γ = G/1000. We see that they almost coincide with the solution for the 4-node element HW47. 2. For MITC9i, we additionally provide the solution for γ = G ("MITC9i,G") and it is locked compared to the one for γ = G/1000. The solution for a 16 × 16element mesh ("MITC9i, 16x16") is softer than that for the 8 × 8-element mesh. Load Inward displacement at force P Concluding the linear and non-linear analyses, we see that the accuracy of the MITC9i shell element is very good compared to the 9-SRI0 element and that the reduced value γ = G/1000 for the drilling rotation term is beneficial in the case of curved elements.

Short C-beam
Drilling rotations are mainly introduced to allow for analyses of intersecting shells and in this example we have two 90degree intersections. A short C-beam is fully clamped at one end and loaded by a vertical force P at the other, see Fig. 19a. At the clamped end, displacements and rotations are constrained to zero, as proposed in [12]. The data is as follows: E = 10 7 , ν = 0.333, the thickness h = 0.05. The web is modeled by 18 × 3 elements and each flange by 18 × 1 elements, so the total number of elements is 90. The vertical displacement u z at the point where the force P is applied is monitored.
The mesh is regular so coefficients of the corrected shape functions (CSF) are equal to zero in this example.
We see that the solutions for MITC9i, 9-AS, 9-SRI and 9-EAS11 are similar and slightly above the reference value. Besides, MITC9i is insensitive to the value of γ , and yields almost the same results for γ = G/1000 and γ = G.   The solution by MITC9 is slightly excessive. Surprisingly, the solution for a basic element 9 does not show much locking.
The non-linear solution computed using the arc-length method with the initial P = 20 is shown in Fig. 20. All elements are tested for γ = G/1000. The curves for 9-AS and 9-SRI coincide with that for MITC9i and the solution for MITC9 is quite close to them. The 4-node HW47 element is Note that the reference elements 9 and S9R5 behave erroneously; the solution for 9 is too stiff while for S9R5 is too flexible, compared also to the solutions from, e.g., [6,12]. Note that S9R5 lost convergence at about u z ≈ −2.7.

L-shaped plate
In this test we compare the convergence properties of the MITC9i element to the 9-node element 9-EAS11 and two 4-node elements, including the mixed/enhanced HW29 element of [48]. This example was also analyzed using quadrilateral and triangular shell elements in [38] and [10].
The L-shaped plate is clamped at one end and the in-plane force P is applied at the other end, see Fig. 21. The data is as follows: E = 71,240, ν = 0.31, w = 30, L = 240, and thickness h = 0.6. The value of the regularization parameter for the drilling RC term γ = G. A mesh of 17 nine-node elements (102 nodes) is used.
The solution of this problem has a bifurcation point at which an out-of-plane deformation occurs. We add a small out-of-plane load and solve the equilibrium problem using the arc-length method. The in-plane force P and the outof-plane load P 3 = P/1000 are applied at point B, at which the out-of-plane displacement u 3 is monitored.
In this test, 10 steps are performed by the arc-length method for the initial load increment P = 0.5 and the requested number of iterations per step I req = 20. The non-linear solutions are presented in Fig. 22 curve was obtained in [48] using 64 4-node elements and P = 0.1. We see in Fig. 22, that the longest steps were made by the 4-node HW29 element, next by the 9-node MITC9i and 9-EAS11, and the shortest ones by the 4-node EADG5A. The total number of Newton iterations used in 10 steps was as follows: (a) 4-node elements: HW29 -67, EADG5A -90, (b) 9-node elements: MITC9i -91, 9-EAS11 -91. We see that the 4-node HW29 outperforms the other 4and 9-node elements in this test. We recall that a treatment of the drilling RC was pivotal to improve convergence properties of this element, see [48]. In this element, we applied the Perturbed Lagrange method and the 3-parameter Lagrange multiplier, which yielded one spurious mode, for which we developed the γ -stabilization. Finally, we note that this analysis can be continued further without any problem, and, e.g., for MITC9i and 20 steps, we end up at the load level of about 130.

Final remarks
The improved nine-node shell element MITC9i with drilling rotations was developed and tested in this paper. The element is based on the Reissner-Mindlin kinematics, Green strain and the extended potential energy functional for the plane stress condition. Its formulation includes several improvements, which we summarize as follows: 1. We propose to formulate the MITC method for the bending/twisting strains κ and the transverse shear strains γ using the 'COVc' components of Green strain, see Sect. 3, which is different from in the classical MITC9 element. This modification renders that the patch tests are passed by the MITC9i element for the regular mesh of Fig. 7. 2. The corrected shape functions are applied instead of the standard shape functions, see Sect. 4. They eliminate distortions of a local coordinate space caused by shifted nodes, which results in a better interpolation accuracy. These functions enable the MITC9i element to pass the patch tests for meshes distorted by parallel shifts of midside nodes and arbitrary shifts of the central node, see Fig. 8. The method of computation of the shift parameters of [9] was generalized from 2D to shells located in 3D space in Sect. 4. To compute the shift parameters of a central node, we propose the method based on the minimization of a square error, Eq. (50), which is suitable for flat, shallow and non-flat elements, see Sect. 6.4. Additionally, we propose an alternative method to calculate the shift parameters for midside nodes, Eq. (58), different from the one of [9]. It does not rely on the proportion of arc-lengths but uses a parametric equation of a parabola. In effect, a symmetric curved side is generated even for a shifted midside node, which improves accuracy when the shell boundary is curved and symmetric, see Sect. 6.5. 3. The drilling rotation is incorporated using the drilling part of the Rotation Constraint equation, Eqs. (10) and (11). We found that the corrected shape functions reduce sensitivity of solutions to the value of regularization parameter γ of the penalty method. Besides, we confirmed that the reduced value γ = G/1000 improves accuracy for curved or distorted element shapes, see Sects. 6.3 and 6.7.
Finally, the shell element MITC9i with drilling rotations not only passes the patch tests described above, but its accuracy is very good in linear tests for coarse distorted meshes. In non-linear tests, see Sects. 6.7, 6.8 and 6.9 (we have also run several other tests, not reported in this paper), it is very robust and its convergence properties are good. Nonetheless, its radius of convergence is smaller than of our best 4-node element HW29, see Sect. 6.9. Therefore, further studies are needed to fully compare its performance to other existing very good quadrilateral and triangular elements.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. j −1 = g 1 · t 1 g 1 · t 2 g 2 · t 1 g 2 · t 2 , and g α are the co-basis vectors.

Transformation of in-plane and transverse components
Transformation formulas for covariant and Cartesian components of the second-rank tensor A can be obtained in the following way.
In-plane (αβ) components First, we calculate the Cartesian components A C αβ = t α · (At β ) and, next, transform them using Eq. (75). Finally, we identify the covariant components A αβ . = g α ·(A g β ) of the tensor representation A = A αβ g α ⊗ g β . The transformations are performed as follows: We can rewrite the above formulas in the vector-matrix form where Alternatively, in terms of 2 × 2 matrices of components A C and A ξ , we have and the covariant components in terms of Cartesian ones, Transverse (α3 and 3α) components The derivation is analogous to that for the in-plane components. First, we calculate the Cartesian components A C α3 = t α · (At 3 ) and, next, transform them using Eq. (75). Finally, we identify the covariant components A α3 . = g α · (A g 3 ) of the tensor representation A = A α3 g α ⊗ g 3 . (Note that for the normal basis, We can rewrite the above formula in the vector-matrix form The formulas for A C 3α components are analogous.