Transverse shear parametrization in hierarchic large rotation shell formulations

Consistent treatment of large rotations in common Reissner–Mindlin formulations is a complicated task. Reissner–Mindlin formulations that use a hierarchic parametrization provide an elegant way to facilitate large rotation shell analyses. This can be achieved by the assumption of linearized transverse shear strains, resulting in an additive split of strain components, which technically simplifies implementation of corresponding shell finite elements. The present study aims at validating this assumption by systematically comparing numerical solutions with those of a newly implemented hierarchic and fully nonlinear Reissner–Mindlin shell element.

Nearly two decades ago, the introduction of isogeometric analysis (IGA) by Hughes et al. 10 triggered a renaissance of the Kirchhoff-Love model, because the usage of splines, particularly non-uniform rational b-splines (NURBS) promises an easier construction of C 1 -continuous function spaces for discretization.Kiendl et al. 11 where the first to exploit this feature for a Kirchhoff-Love type IGA shell formulation.
Soon, however, it has been found that C 1 -continuity can also be beneficial for the Reissner-Mindlin model.Long et al. 12 presented a shell formulation based on subdivision surfaces, in which the total rotation of the shell director is split into a contribution of the rotation of the surface normal and the transverse shear part.Before the development of IGA, the same group had already presented a Kirchhoff-Love shell formulation based on subdivision surfaces in Cirak et al. 13 as early as 2000.A concept similar to Long et al. 12 was developed in Echter et al. 14 for a hierarchic Reissner-Mindlin IGA shell formulation.Because derivatives of the displacement field are involved in the construction of the rotated normal, C 1 -continuity is required for discretization.
As the transverse shear part is described by distinct degrees of freedom, these formulations are intrinsically free from transverse shear locking, without any additional measures.Indeed, a Kirchhoff-Love formulation can directly be extracted by simply constraining the transverse shear degrees of freedom to zero.This particular feature is shared by the mixed element formulation of Auricchio and Taylor. 15,16Therein, the specific construction of the underlying multi-field functional allows to "switch off" the transverse shear part, thus combining the Reissner-Mindlin and the Kirchhoff-Love model in one unique hierarchic formulation.It has been extended to shells by Bischoff and Taylor. 17he formulation of Echter et al. 14 is limited to geometrically linear problems and uses a difference vector to represent the change of direction of the director related to transverse shear deformation.The components of this difference vector can be easily transformed into rotation angles and in fact this formulation is aptly characterized as a formulation with hierarchic rotations.As already mentioned, it is intrinsically free from transverse shear locking.However, if equal order interpolation is used for the mid-surface displacements and the components of the difference vector (a.k.a.hierarchic rotations) oscillations of the transverse shear forces can be observed.
Oesterle et al. 18 succeeded in overcoming this problem by proposing a completely rotation-free Reissner-Mindlin shell formulation, in which the hierarchic part, describing the shear deformation, is obtained from derivatives of additional displacement degrees of freedom that can be interpreted as "shear displacements."This still geometrically linear formulation is free from transverse shear locking and does not exhibit any oscillations.Henceforth, this type of model is denoted as formulation with hierarchic displacements.
Eventually, Oesterle et al. 19 presented geometrically nonlinear versions of the hierarchic rotation and hierarchic displacement formulations.They are capable of dealing with large rotations and large membrane and bending strains.The transverse shear part, however, is linear.In other words, total rotations can be large, but the shear angles are assumed to be small.The reasons for this restriction are mostly technical.In particular, it dramatically simplifies linearization, because there is an additive split of strains emanating from the bending and transverse shear parts.Unlike conventional nonlinear Reissner-Mindlin shell elements, this concept does not require to deal with rotation tensors and algorithms for large rotations.
Recently, Neunteufel and Schöberl 20 adopted the hierarchic rotation approach to enhance their previously presented nonlinear Kirchhoff-Love shell element by transverse shear strains to obtain a Reissner-Mindlin shell element.The formulation is based on the mixed Hellan-Herrmann-Johnson method with a focus on membrane locking.
Quite naturally, the question arises, whether the restriction to small transverse shear angles would compromise the predictive capabilities of the model in certain situations.It may well be possible that for relatively thick shells or very large gradients of bending moments or composite shells with a soft core, large shear angles and corresponding geometrically nonlinear effects are not negligible.Although this problem has already been mentioned to a certain extent, 19,21 it is the major motivation for the present study.
Interestingly, nonlinear Reissner-Mindlin shell elements with linearized shear rotations can also be encountered in commercial finite element software, although this is barely noted in the literature.The Belytschko-Lin-Tsay element 22 is one of the workhorses in the commercial code LS-DYNA 23 for explicit dynamic analyses, for instance in crash and metal forming simulations.As described in Belytschko et al., 22 it uses the assumption that the angle between the director and the normal remains small.For the shear deformable shell element SHELL181 in ANSYS, the underlying mathematical formulation is not disclosed in the theory manual 24 and the authors are not aware of any detailed theoretical description of the element.The documentation indeed mentions that "SHELL181 includes the linear effects of transverse shear deformation." 24A simple numerical experiment, applying prescribed nodal rotations, while fixing all other degrees of freedom, along with the option of a linear elastic material, reveals for both elements a linear relation between shear stress and rotation, regardless of the magnitude of the rotations.
The reasons or the motivation for these simplifications in commercial codes are, to the authors' knowledge, not named.It is, however, well known that consistent treatment of large rotations is an awkward task, because, unlike displacements, large rotations do not live in a linear vector space, see Zienkiewicz, Taylor and Fox 9 for a concise theoretical summary and Müller and Bischoff 25 for a recent contribution and an overview of the state of the art in this area.At least for the concept followed in LS-DYNA, it can be said that it successfully circumvents dealing with rotation tensors and corresponding update algorithms for large rotations.
The present study provides a systematic investigation of the effect of linearized transverse shear parametrization in the context of hierarchic shell elements by means of numerical experiments, devised to trigger potential nonlinear phenomena associated with (large) transverse shear deformations.

Differential geometry and shell kinematics
The presented shell formulations use a convective curvilinear coordinate system with in-plane coordinates   , which span the shell mid-surface, and a thickness coordinate with shell thickness t.Greek indices take on values of 1 or 2, while Latin indices run from 1 to 3. Partial derivatives w.r.t. the curvilinear coordinates are written as ( ) ,i = ( )  i .Quantities of the reference configuration are represented by capital letters, whereas small letters refer to the current configuration, compare Figure 1.
Using dimensional reduction and a total Lagrangian setting, the position of an arbitrary point of the shell body is described by its mid-surface position and the director: The reference director A 3 is defined as the normalized cross product of the in-plane covariant base vectors A  : The definition of the current director a 3 depends on the shell formulation and is given in Sections 2.2 and 2.3.
F I G U R E 1 Shell body in reference configuration and current configuration.
Thus, the displacement field for a point of the shell body u is defined by the difference of position vectors in reference configuration and current configuration.With dimensional reduction, the deformation is split into mid-surface displacements v and the difference between directors: The three-dimensional Green-Lagrange strain tensor is defined as with the covariant base vectors of the shell body Using Equation ( 3), utilizing orthogonality (A  ⋅ A 3 = 0 and a 3, ⋅ a 3 = 0), and neglecting terms quadratic in  3 , the strain components read ) , ) , The current covariant base vectors of the mid-surface can be constructed by While constant parts of E  represent membrane strains, the parts linear in  3 represent curvatures.E 3 are the shear strains, which vanish in the case of shear rigid formulations, in which the current director a 3 is orthogonal to the in-plane base vectors a  .

Hierarchic shell formulation with linearized transverse shear parametrization
The shear rigid Kirchhoff-Love (KL) shell formulation with nonlinear kinematics, which serves as a basis for our hierarchic shear deformable shell in Oesterle et al., 19 describes the director in the current configuration as the current normal, using the in-plane base vectors From the orthogonality of the normal and the in-plane base vectors a  ⋅ a ⊥ 3 = 0, it follows by differentiation w.r.t.  that a derivative of the normal can be moved to the in-plane basis: Using the identity which follows from differentiation of the definition in Equation ( 7) w.r.t.  , the strain components for the Kirchhoff-Love formulation simplify to ) , ) , For our geometrically nonlinear shear deformable Reissner-Mindlin shell formulation with linearized transverse shear parametrization 19 (RM-LS), the director is constructed in a hierarchic manner, by superimposing shear in form of a shear difference vector w to the deformed normal, compare Figure 2: This construction technically violates the inextensibility condition 3 || holds.However, the hypothesis in Oesterle et al. 19 is that in many problems shear rotations are relatively small compared to total rotations.Thus, the idea is to consider shear rotations in a linearized way, while bending and membrane parts of the deformation can be large, that is, they fully include geometrically nonlinear effects.A benefit of the proposed construction is that simplifications due to the orthogonality of a  and a ⊥ 3 , that were made for the components of the strains in Equation (11), can be reused for this formulation.These membrane and bending related strain components are complemented by additional parts, responsible for transverse shear contributions.Consequently, the strain components show an additive structure: ) , The additive structure of the strain components nicely visualizes the ability of this formulation to avoid transverse shear locking a priori, that is, before discretization, compare also Section 2.4.Furthermore, it is beneficial for the derivation of finite elements, as described in Section 2.5.

Hierarchic shell formulation with nonlinear transverse shear parametrization
A fully nonlinear description of the director of the shear deformable Reissner-Mindlin formulation (RM-NL) using a difference vector is proposed by Long et al., 12 compare Figure 2: By normalization of the total director, the inextensibility condition is fulfilled, that is, ||.However, the additive structure of the director is lost, and simplifications are impeded.The strain components remain similar to the ones described in Equation ( 6): ) , ) , ) , Although the director and the strain components do not show a purely additive structure, as opposed to Equations ( 12) and ( 13), the formulation is free from transverse shear locking, as explained in the following Section 2.4.

Parametrization of the hierarchic shear difference vector
The shear difference vector for both linearized and nonlinear shear parametrization is constructed using the in-plane base vectors and the primary variables w  , which are characterized as hierarchic rotations: For the RM-LS formulation, they are directly related to linearized shear angles, if the in-plane base vectors in Equation ( 16) are additionally normalized.For the RM-NL formulation, such geometrical relations for w  are more cumbersome to establish, due to the non-additive construction of the director in Equation (14).With hierarchic rotations as distinct primary variables for transverse shear contributions at hand, it is obvious that bending deformations without transverse shear can easily be achieved (cf.Appendix A).Thus, in contrast to standard Reissner-Mindlin formulations, the hierarchic RM-LS and RM-NL formulations are intrinsically free from transverse shear locking, that is, transverse shear locking is avoided on formulation level instead of removing its effects after discretization on element level.However, as described in depth for the model problem of a Timoshenko beam in Oesterle et al., 18 structural element formulations with hierarchic rotations still show an imbalance of function spaces in the kinematic equations for every equal order interpolation.While shear strains are fully balanced and, thus, shear locking is avoided, curvatures show different derivatives of the primary variables.The effect of this, herein called, bending locking in hierarchic rotation formulations for static analyses is that shear stress resultants still show some oscillations at the domain boundaries for low slenderness.Consequently, an alternative construction of the shear difference vector (16) was introduced in Oesterle et al. 18 and applied to large rotations in Oesterle et al. 19 It uses so called hierarchic displacements v s  instead of hierarchic rotations as primary variables.With these, the hierarchic rotation angles are described by derivatives of the hierarchic displacements w  = v s  , , which removes all imbalances from the kinematic equations and further improves the quality of stress resultants.However, since the RM-NL formulation of Long et al. 12 is originally parametrized by hierarchic rotations, we use the same for the RM-LS formulation to consistently compare the two formulations regarding the effects of linearized transverse shear.

Derivation of element vectors and matrices
Thinking further in derivation of displacement based finite shell elements, the vector of internal element forces is built from the gradient of the strain energy, and the tangent stiffness matrix is the Hessian of the strain energy.Thus, first and second derivatives of the strain components w.r.t. the five degrees of freedom u r are needed: ( ) ,r = ( ) u r and ( ) ,r,s =  2 ( ) u r u s with r = 1 … 5 and s = 1 … 5. Due to the construction of the director or, respectively, the normal in Equations ( 8), (12), and ( 14), this involves the use of the quotient rule.Comparing the strain components in Equations ( 11), (13), and (15), it is seen that the RM-NL formulation is the only one that needs derivatives of a quotient at this stage.So, while for KL and RM-LS formulations the highest needed derivative of the director or, respectively, of the normal, is a ⊥ 3,r,s , which then already consists of five quotients (cf.eq.(5.32) in Kiendl 26 ), the RM-NL formulation needs one further derivative a RM-NL 3,,r,s .Therefore, manual derivation of the element vectors and matrices is more complex and error-prone for the RM-NL formulation.As an alternative, automatic differentiation can be used, as recently done by Oberbichler et al. 27 or Vigliotti and Auricchio 28 for other elements.

Overview
In this section, we investigate three geometrically nonlinear shell problems in order to study the effect of transverse shear deformation on the structural behavior.In particular, we scrutinize the assumption of small shear rotations.In the following, the shear deformable Reissner-Mindlin shell element with linearized transverse shear parametrization is denoted as RM-LS (in Oesterle et al. 19 this element is found with the abbreviation RM-hr).The newly implemented, geometrically fully nonlinear Reissner-Mindlin shell element on the basis of Long et al. 12 is abbreviated as RM-NL.To verify its correct implementation, we compare solutions to another fully nonlinear Reissner-Mindlin shell element RM-PBFE, 25 which is parametrized by total rotations, and which is available in our in-house research code Ikarus 29 ("PBFE" is for "projection based finite elements").To contextualize any difference between the two shear deformable elements RM-LS and RM-NL, we add solutions for a shear rigid Kirchhoff-Love shell element (KL) on the basis of Kiendl et al. 11 Furthermore, in the third example we draw comparisons with 2D elements.All elements use NURBS-based shape functions and an equal order interpolation of all primary variables with polynomial degrees p and q for the two dimensions.All problems assume linear-elastic, isotropic material behavior using a St. Venant-Kirchhoff material law.The zero transverse normal stress condition S 33 = 0, applied for Kirchhoff-Love and Reissner-Mindlin shell formulations, is incorporated into the material law as usual by elimination of the transverse normal strain E 33 via static condensation.We compute internal force vectors and tangent stiffness matrices of the shear deformable elements RM-LS and RM-NL with automatic differentiation from the implemented strain energy using the C++ library autodiff. 30

Uniaxial bending
The first problem examines uniaxial bending of a beam-like structure.The problem setup is shown in Figure 3.We choose different thicknesses t, which result in different slenderness ratios L t .In addition to L t = 10 and L t = 100 from Oesterle et al., 19 we include a very thin beam with L t = 1000 and a very thick beam with L t = 4.For the latter, the ratio of the analytical solutions for the midpoint displacement by a Timoshenko beam theory w.r.t.Bernoulli beam theory is ) −2 = 1.1, so the difference between shear rigid and shear deformable theory in the geometrically linear case is 10%.For the given geometrically nonlinear case, containing the membrane deformation, which is not triggered in a linear analysis of this problem, this difference is expected to be smaller.However, the chosen statically determinate support keeps membrane effects as small as possible.For all shell elements, sufficiently fine meshes are ensured with the help of mesh convergence studies.For this purpose, a variable number of cubic NURBS elements in x-direction is used for discretization, while in y-direction only one quadratic NURBS element is used.The mesh density is considered sufficient, if further refinement does not change the displacement solution with nine digits precision, as shown in Table 1.Effects from overstiff structural behavior due to transverse shear locking for the RM-PBFE element, bending locking for the RM-LS and RM-NL elements and membrane locking for all elements are thereby excluded.In x-direction the number of required elements is 80, 120, 160, or 200, depending on the type of shell element and the slenderness ratio.With 160 elements, the RM-PBFE element needs a finer mesh than RM-LS and RM-NL (80 elements) for the low slenderness ratio of 10.For the more slender geometry with L t = 100, a coarser mesh with 80 elements is sufficient for RM-PBFE, while RM-LS and RM-NL need 120 and 160 elements.For the more extreme slenderness ratio values of 4 and 1000, the number of required elements is in the same range for all shell elements.
Table 1 shows the maximum vertical displacement for the various slenderness ratios.For the very thin beam, solutions of all shear deformable elements are in perfect agreement, since transverse shear deformations are practically insignificant for the displacement response.The identical results of RM-NL and RM-PBFE for all slenderness ratios verifies correct implementation of RM-NL.Comparing solutions of the KL element to all RM elements, it is seen that with lower slenderness ratios, that is, with larger thicknesses, transverse shear shows a greater influence on the solutions.The difference between shear rigid and shear deformable solutions for L t = 4 with theoretically 10% for the linear case is reduced to 5% in the nonlinear setting here between KL and RM-NL elements.Compared to this, the difference between RM-LS and RM-NL with 0.2% for the very thick beam can be considered as negligible.Thus, in this range, linearization of transverse shear in the RM-LS formulation is legitimized by these results.

Snap-through of a ring
The second problem deals with large total rotations.A closed ring is subject to a prescribed rotation about the z-axis at its top point while the bottom is fixed, so that after snap-through at a rotation of 2, a deformed configuration with three overlapping rings results. 31Figure 4 shows on the left a perspective view of the problem setup with the ring in y-z-plane and its width in x-direction.The concrete definition of the applied Dirichlet boundary conditions at the top and bottom control points is shown in Figure 4 on the right with a view onto the shell mid-surface, where the z-axis is placed at the shell boundary.The ring is discretized by two patches with 50 × 1 elements each.In circumferential direction, a polynomial degree of p = 3 is chosen, whereas in width direction q = 2. Coupling of the two patches is realized using the bending-strip method. 32Due to the pronounced imperfection sensitivity of this problem, the choice of stiffness for the bending strip crucially influences numerical stability of the iterative solution process.For the shear rigid KL element, the directional bending stiffness for the bending strip is set to E s = 5C 11 with C 11 being the entry of the linear-elastic material matrix C that relates curvature E lin 11 to stress S 11 for the given material parameters.For the shear deformable elements RM-LS and RM-NL, directional bending stiffness is set to E s = 4C 11 .Different values may lead to solutions on secondary branches, which do not show symmetric snap-through but rather a lateral evasion of the ring.with the ratio of radius to bending stiffness is taken from the literature, where this example is studied using beam elements.The value EI 3 = 1 12 EtL 3 x describes the bending stiffness w.r.t. the axis in thickness direction  3 .The difference between the solutions of KL and the shear deformable elements makes clear that considering transverse shear in this problem is relevant.We suspect that this difference is related to different ability of the models to represent the torsion of the beam-like structure.However, whether transverse shear is considered in a nonlinear manner or in a linearized manner, appears to be irrelevant here, as displayed by the perfectly coinciding solutions of RM-LS and RM-NL elements.Figure 5 (right) shows snapshots of the deformation process for five different angles of rotation.

Simple shear
The motivation for the third problem is to trigger significant shear deformations.The idea is based on a flat shell with its upper bounding surface subject to horizontal static traction while the lower bounding surface is fixed to a rigid base, compare Figure 6.
Since the solution of the problem is homogeneous in x-direction, a discretization with one single shell element is sufficient.A representative element with in-plane dimensions L x = 2 and L y = 2, compare Figure 6, is chosen.The effect of static traction is applied by inhomogeneous Dirichlet boundary conditions on the shear degrees of freedom w 1 = ŵ1 while all mid-surface displacement degrees of freedom are fixed, compare Figure 7 (left).Consequently, any cross-sectional rotation of the element is a pure shear rotation.
Figure 7 (right) shows the geometrical relation between the relative displacement u rel between the upper and lower bounding surface of the shell and the rotation of the cross section  y .For the RM-LS element this geometrical relation is described by the tangent function, whereas for RM-NL it is the sine function.For a relative displacement of u rel ≈ 0.6, the solutions of RM-LS and RM-NL elements begin to significantly differ from each other.The assumption of linearization of transverse shear can therefore be assumed to be invalid for shear rotations larger than approximately 15 • .However, in this deformation range, the underlying assumption of the Reissner-Mindlin shell model, namely cross-sectional fibers remaining straight during deformation, has to be questioned.
In order to study a possible warping of the cross section, we follow the concept of using a cubic ansatz in thickness direction, as it is used for a higher order 3D shell element in Willmann et al. 33 The problem is discretized by one single 2D NURBS element in the x-z-plane, using linear shape functions in x-direction and a cubic ansatz in z-direction.To isolate the warping effects from additional transverse strain effects (thickness change), inhomogeneous Dirichlet boundary conditions are applied in such a way that the cross-sectional length remains unchanged, compare Figure 8 (left).The deformation plot in Figure 8 (right) clearly shows the cross-sectional warping already for shear deformations smaller than u rel = 0.6.Qualitatively similar results for the deformations can also be obtained by a discretization with bi-linear 2D elements when using multiple elements in z-direction.Hence, it can be concluded that during this deformation process  the basic assumption of straight cross sections becomes invalid before the assumption of linearized transverse shear leads to deviations from the three-dimensional solution.In other words: It does not seem to be sensible (or necessary) to take into account geometrical nonlinearity for transverse shear if at the same time the assumption of cross-sectional fibers remaining straight is maintained.

CONCLUSIONS AND OUTLOOK
In this article, we critically questioned the assumption of small shear angles, which was made earlier for geometrically nonlinear hierarchic Reissner-Mindlin shell formulations in Oesterle et al. 19 This assumption leads to an additive split of strains in Equation ( 13), where parts related to shear deformation can be added to the strain components of a Kirchhoff-Love formulation.This simplifies the derivation of nonlinear finite shell elements compared to using a nonlinear representation of transverse shear.We have studied the effect of linearized transverse shear compared to nonlinear shear using NURBS based elements on the basis of hierarchic rotation formulations.We emphasize that the same conclusions hold for the hierarchic displacement Reissner-Mindlin formulation also derived in Oesterle et al., 19 which also utilizes the assumption of linearized transverse shear.
For the numerical problem of uniaxial bending, we concluded that even for a very thick beam the influence of linearization of transverse shear is negligible when compared to the difference between shear deformable and shear rigid solutions.The problem of a ring snapping through showed that in a large total rotation problem, in which transverse shear deformations are in general relevant, there is no visible difference between solutions of the linearized transverse shear elements and the fully nonlinear ones.The simple shear problem, motivated by a static surface traction, showed that the assumption of linearized transverse shear is justified for problems with transverse shear deformations that are in a range where the assumption of cross-sectional fibers remaining straight during deformation is met.
In summary, we conclude that linearization of transverse shear strains in large rotation shell elements based on the Reissner-Mindlin model is a valid method.
When larger cross-sectional deformations evolve, which are in particular due to nonlinear constitutive behavior, for example, in large strain elasto-plasticity, further investigations are necessary extending the kinematics to shells with thickness changes, here named as 3D shells, or to shells including cross-sectional warping, here named as higher order 3D shells, as presented for example by Willmann et al. 33 for sheet metal forming.Linearization of transverse shear inside a 3D shell formulation leads to a consequent continuation towards nonlinear kinematics of the hierarchic family of isogeometric shell finite elements, which was introduced in Echter et al. 14 for linear kinematics.The application of 3D shells or higher order 3D shells to laminated composites or sandwich structures with stiff outer layers and soft cores potentially gives also rise to large transverse shear deformations.Building on the results from the simple shear problem in Section 3.4, however, the authors assume that warping effects will be more relevant than the representation of nonlinear transverse shear also in these cases.

F I G U R E 2
Mid-surface of the shell in reference configuration and current configurations for hierarchic shear deformable shell formulations with different shear parametrizations.

F
I G U R E 4 Snap-through of a ring, problem setup.F I G U R E 5 Snap-through of a ring, load-displacement diagram for scaled reaction moment over rotation angle (left) and five stages of the deformation process (right).The load-displacement diagram in Figure 5 (left) shows the scaled reaction moment plotted versus the rotation angle φz .The reaction moment M is calculated from the pair of reaction forces of the control points at the bottom.Its scaling MR EI 3

F I G U R E 6
Simple shear, motivation.F I G U R E 7 Simple shear, problem setup for shells (left) and relative displacement over cross-sectional rotation (right).

F I G U R E 8
Simple shear, problem setup for 2D element (left) and deformations at various stages (right).
Uniaxial bending, maximum vertical displacement for various slenderness ratios.
Uniaxial bending, problem setup.TA B L E 1Note: Element names highlighted in color for better recognition.