ON THE CHOICE OF THE DIRECTOR FOR ISOGEOMETRIC REISSNER-MINDLIN SHELL ANALYSIS

The contribution deals with the treatment of the director vector and the rotational degrees of freedom in isogeometric Reissner-Mindlin shell analysis. In NURBS-based formulations the director vectors and the nodal coordinate systems are defined at the control points. The interpolation with the basis functions maps them to the shell reference surface. For curved geometries the interpolated director at an arbitrary point is in general neither perpendicular with respect to the shell surface nor has the correct length. Thus, basic kinematic assumptions are violated. It is shown that this is one reason for the observed degradation of numerical results for shells when the polynomial order gets higher than five. In the present work a method to overcome these problems is provided. A geometrically exact Reissner-Mindlin shell formulation with five degrees of freedom, three displacements and two rotations, is presented. The two rotations are defined with respect to the orthogonal tangent vectors at the shell reference surface. In standard formulations the nodal coordinate system which also defines the director is constructed from the derivatives in the closest point projection of the control point. It is observed that error norms of angle and length of the interpolated director increase with the order of the basis functions. To improve the accuracy the nodal coordinate systems and directors are chosen in a way, that the interpolated basis systems in all integration points are orthonormal and in correct orientation. The existence of an unique solution is shown. This leads to nodal coordinate systems which are no longer orthonormal. An interpolation, which ensures orthonormality at the integration points, is introduced. The finite element implementation is tested with standard benchmark problems. In contrast to computations with nodal directors being the normal to the shell surface, deformations converge monotonically. In all computed examples the accuracy of the solution for a given mesh increases with the order of the basis functions. Simplifications in the underlying shell theory inhibit exponential convergence rates for curved domains. However, the proposed measures significantly improve the accuracy of computations with high order basis functions.


INTRODUCTION
Isogeometric Analysis is a high order alternative to standard finite element methods.Introduced in 2005 by Hughes et al. [8] it attracts much interest recently.The idea of isogeometric analysis is to use the same mathematical geometry description for modeling and computations.This averts costly conversions.Most approaches base on non-uniform rational B-Splines (NURBS) which are widespread in the graphical industry.These functions posses the required mathematical properties as shown by Bazilevs et al. in [1] for the use as basis functions in a finite element context.Geometries modeled by NURBS are usually described by surfaces.This requires a shell implementation to avoid a conversion to a solid model.The continuity properties of NURBS surfaces permit the use of Kirchhoff-Love shells, entailing a few drawbacks regarding boundary conditions and connection of multiple NURBS patches.An implementation is described in [9].To avoid these drawbacks and to allow for thick shells a Reissner-Mindlin shell formulation is chosen.A formulation based on a degenerated solid with director vectors obtained by a closest point projection is described in [4].A rotation-free isogeometric shell based on a discrete Kirchhoff -type formulation is proposed by Benson et al. in [3].There additional factors computed on the element level to attain an exact interpolation of the normal vectors are introduced.
In this work we present a NURBS-based isogeometric Reissner-Mindlin shell formulation with exact director vectors.The nodal director vectors valid throughout the entire domain are computed in a preprocess with a global system of equations.First, all required NURBS terminology and the employed Reissner-Mindlin shell theory are outlined in section 2 and 3.In section 4 the computation of the director vector and the nodal basis systems is explained in detail.This is followed by the description of the isogeometric finite element formulation.Numerical examples demonstrate the capability of the presented formulation.

NURBS-BASED ISOGEOMETRIC ANALYSIS
Geometry description and interpolation of the presented shell formulation is based on NURBS basis function in an isogeometric framework.NURBS terminology needed in the following is shortly introduced.For more details see [13] and [6].NURBS surface patches are spanned by a tensor-product combination of two knot vectors In ξ 1 -direction the order of the NURBS basis functions is referred to as p and the number of control points in this direction is termed n.Analogously, q and m define the ξ 2 -direction.Thus, a control point net with n np = n × m control points arises.The control points are projected to a 3-dimensional space with the NURBS basis functions N I , where the index I is a specified function of the indices i and j.The NURBS basis functions are computed with the weighting factor w ij of the control point.The B-Splines basis functions N p i (ξ 1 ) and N q j (ξ 2 ) are determined by the recursive Cox-de Boor formula, see [13].The number of nonzero basis function in every element is n en = (p + 1) (q + 1).The number of elements per patch is n el = (n − p) (m − q).For the determination of a physical point on the NURBS surface only n en specified basis functions have to be considered.See [13] for the derivatives with respect to the parametric coordinates ξ α with α = 1, 2.

REISSNER-MINDLIN SHELL THEORY
The proposed shell element incorporates a linear Reissner-Mindlin shell theory with an inextensible director field to account for transverse shear strains.The shell formulation is directly derived from the continuum theory as proposed in [11].The variational formulation is pure displacement based and no measures against locking are taken.As the focus of this contribution lies on the choice of the director vectors a linear theory for small strains and small rotations with linear elastic material is implemented.

Kinematic description
Location X and deformations û of a physical point on the shell are given by where ξ α are the parametric coordinates of the shell mid-surface and ξ 3 is the thickness coordinate.Here Greek indices range from 1 to 2 and Latin indices from 1 to 3. In the linear theory reference and current configuration coincide.The position vector of the mid-surface is denoted by X (ξ α ).The deformed configuration computes with the deformations u (ξ α ) to x (ξ α ) = X (ξ α ) + u (ξ α ).The director vector is perpendicular to the shell mid-surface and is in accordant orientation with the thickness coordinate ξ 3 with the range − h 2 ≤ ξ 3 ≤ h 2 .The vector d is the difference vector between the director vector and the deformed director vector.It can be computed by d = ω × D. The covariant basis vectors G α are defined as the derivatives of the position vector X with respect to the coordinates ξ α .This results in where subscript commas denote derivation with respect to the parametric coordinates ξ i .The contravariant basis vectors G i can be computed by with the Kornecker-delta δ j i .The mid-surface area of an element is

Stresses and strains
The strains are measured with the linearized Green-Lagrange strain tensor For an efficient formulation the strains are split up into shell strains with under disregard of second order curvatures.Thus the resulting membrane strains ε αβ , the curvatures κ αβ and the shear strains γ α are with the derivatives and is work-conjugate to the resulting stress vector σ = n 11 , n 22 , n 12 , m 11 , m 22 , m 12 , q 1 , q 2 T .(13)

Variational formulation
The presented shell is based on a pure displacement formulation with static loading, which is considered as Neumann conditions t on the boundary Γ t and surface loads p.The potential reads with the independent displacement vector v = [u, ω] T containing the displacements u and the rotational parameters ω of the shell mid-surface.The variation of the potential leads to the weak form of equilibrium and is given as where C is the linear elastic material tensor.The virtual shell strains are arranged akin to equation (12).

Closest point projection
Reissner-Mindlin shell formulations require the definition of nodal basis systems for the interpolation of the director vector and for the rotational directions.These systems have to be given as an input.In standard finite element shell formulations these nodal basis systems are e.g.determined by analytical functions for regular geometries or they are interpolated from the neighboring elements.For more details about the choice of directors in standard finite element shell formulations see [5].In isogeometric analysis the control points, which correspond to the nodes in finite element formulations, are in general not located on the physical domain.The control point net is disjunct from the element mesh.Thus in [4] the use of the normal from the control point with respect to the shell reference surface is proposed.This normal vector has to be computed in a preprocess by a closest point projection.The characteristics of NURBS surfaces preclude a closed solution, an iterative solution process, e.g. with the Newton-Raphson method is necessary.See [13] for details.With the normal vector and the local derivatives in the projected point orthonormal nodal basis systems can be determined.Computations with an isogeometric Reissner-Mindlin shell formulation along the lines of [7] do not show correct convergence behavior for order elevation (k-refinement).This is also visible in the results of Benson et al. in [4].As an example the Scordelis-Lo roof, part of the shell obstacle course proposed by Belytschko et al. in [2], is examined.System and applied load are exactly chosen as in [2].In Figure 1 the h-refinement convergence behavior for three orders p of NURBS basis functions is given.A system sketch is also provided for convenience.Figure 1 clearly shows that the accuracy deteriorates for rising orders p.The distance of control points from the relevant element and the number of control points influ-encing one element n en = (n − p) (m − q) grows with the order p respective q.In the case of curved domains this leads to an addition of vectors with distinct differing directions which has a negative impact on the interpolation of the director vector.This can be clearly seen in Figure 2, where the error norms for deviation of length and angle of the director vector are plotted for three different orders p of NURBS basis functions.This deterioration in interpolation quality averts proper convergence behavior.

Computation of exact basis systems
In the following section a method is introduced, which ensures the orthonormality and the correct orientation of the director vector at the integration points.The orthonormal basis vectors A i entail nine equations in each of the n gp integration points in each of the n el elements of the relevant patch.As all nine directions can be treated independently this yields (n gp × n el ) equations of the the type which can be arranged in a system of equations Here A GP ij is the j-th component of the i-th basis vector of the integration point basis system.Analogously A CP ijI are the unknown values for the nodal point basis systems at the control point I.The matrix N with (n gp × n el ) rows and n np columns contains the basis functions N I and is identical for all nine directions.This requires only one solution of equation ( 18) with multiple right-hand-sides.The number of integration points per element n gp is chosen as for the structural analysis and must lead to a determined or overdetermined system of equations.The condition n np ≥ n gp × n el holds in all cases where the structural analysis is stable.As the columns of N contain the linear independent basis functions N I the columns of the matrix are linear independent.This implies that holds, which guarantees an unique solution, see e.g.[10].This allows to solve the determined normal equation instead of the overdetermined equation (18).

ISOGEOMETRIC FINITE ELEMENT IMPLEMENTATION
The shell mid-surface, the deformations and the director vector are interpolated by with the NURBS basis functions given in equation ( 2).The superscript h identifies approximated quantities.The director vectors D I = A 3I are computed in the preprocess described in section 4. The derivatives of the NURBS basis functions N I,ξ α with respect to the parametric coordinates ξ α have to be transformed to a local Cartesian basis system t.With the resulting derivatives N I,α the quantities are interpolated.For the variation of the potential the interpolation of the variated shell strains is necessary.It results in with The interpolation of δd h ,α and δd h is given in the following.The proposed approach does not interpolate the nodal variation of the director, but the local rotations δβ.It is somehow similar to the full SO(3) update proposed in [12], where the interpolation of the rotational parameter δω is discussed.However the proposed approach is more accurate for curved domains as account is taken for the differing orientation of the nodal rotations.
The variation of the director vector can be rewritten to with W = skew D. The equivalence δw = δω holds for small rotations.With the interpolations of the nodal basis systems the matrix establishes the connection δω = T 3 δβ.Thus follows as relation between variated director and variated nodal rotations.The interpolation of δβ requires a transformation from local rotations in the control points δβ I to local rotations in the integration point.From can be derived.The drilling rotation δβ 3 is fixed for a 5-parameter shell.Therefore a statical condensation is possible and δβ 3I can be eliminated from equation (30).As result the interpolation with is established.The combination of equation ( 28) and (31) yields as interpolated variated director vector.The derivatives δd ,α are more lengthy to attain.They read with and the derivatives of the transformation matrix Inserting equations ( 24), ( 33) and (34) into equation ( 23) yields the B I matrix, which establishes the relation between virtual strains and the variations of rotations and displacements with the submatrices b Iα = M T I T T X ,α and bαβI = TT I,α X ,β .The matrix TI,α is defined by The strains ε h are interpolated akin to equation (34).The approximation of the weak form (15) is given as (40)

NUMERICAL EXAMPLES
The improved accuracy of the proposed approach is shown with the help of two numerical examples.Both benchmarks are part of the shell obstacle course proposed by Belytschko et al. in [2].The Scordelis-Lo roof illustrates the difference between the proposed approach and a standard approach using closest point projection and orthonormal nodal basis systems.In the second example the present shell formulation is compared to the isogeometric Reissner-Mindlin shell described in [4].All computations are linear.

Scordelis-Lo roof
In section 4 the Scordelis-Lo roof is used to illustrate the shortcomings of Reissner-Mindlin shell analysis with directors attained by closest point projection and interpolated variated nodal directors.The latter one is called standard approach and is compared to the present approach.Figure 3 shows the convergence behavior of a characteristic deformation with respect to mesh refinement (h-refinement) for different polynomial orders.The symmetry of the system is used and only one quarter of the system is computed.The length of the full system is L = 50 with a radius R = 25 .The wall thickness is t = 0.25 with a Young's modulus E = 4.32 • 10 8 and a Poisson's ratio ν = 0 .The system is loaded with a uniform gravity load of g = 90 per unit area.For more details about the geometry of the system see [2].The deformations computed with the proposed approach clearly show monotonously growing convergence behaviour for all order of NURBS basis functions.The standard approach overestimates the defomations for a order p = 4 for very coarse meshes, then shows oscillatory behaviour and finally converges from above.For an order p = 6 the standard approach converges slower than for p = 4, which indicates wrong p-convergence behaviour for the used k-refinement.For low order NURBS basis functions only small differences between both approaches occur.

Pinched cylinder
The pinched cylinder sketched in Figure 4 is used to compare the accuracy of the proposed formulation with existing shell formulations.As a reference the isogeometric Reissner-Mindlin shell proposed by Benson et al. in [4] is used.The system consists of a cylinder with a radius R = 300 and is L = 600 long.The material is specified by a Young's modulus E = 3.0 • 10 6 , Poisson's ratio ν = 0.3 and t = 3 is the wall thickness.The system is loaded with two opposing radial forces F = 1 .Both ends of the cylinder are constrained by a rigid diaphragm.Due to symmetry only one eighth of the system is modeled.More details can be found in [2].The deformations given in Figure 4   Computations with the present formulation show improved accuracy with rising order p of the NURBS basis functions, whereas in the reference computations order p = 4 outperforms p = 5.Furthermore the accuracy of the present approach is superior to the reference results when the same order of basis functions and mesh discretization is used.

CONCLUSION
In this contribution a Reissner-Mindlin shell formulation employing NURBS basis functions for the description of geometry and director vectors and the interpolation of deformations is introduced.The formulation allows for an exact evaluation of the initial director at an arbitrary point of the shell geometry.The proposed interpolation of the variated local rotations is more accurate than the common interpolation of the variated director vector.Whereas the difference in standard finite element formulations is negligible, the provided numerical examples have shown that the accuracy of isogeometric formulations with higher order is massively impaired by such an approach.The present approach yields superior results and thus allows the use of coarse meshes with high orders of basis functions.
The extension of the presented formulation to nonlinear computations is straightforward and will be treated in a subsequent paper.Further research will be necessary to handle shell with kinks.Measures against the various effects of locking and the consideration of second order curvatures could further improve the results for coarse meshes with high order basis functions.

Figure 1 .
Figure 1.Scordelis-Lo roof: System and deformations for director vectors determined by closest point projection.

I
CB K dΩ .

Figure 3 .
Figure 3. Scordelis-Lo roof: Comparison of deformations in point A.
are normalized to a value of 1.83 • 10 −5 .

Figure 4 .
Figure 4. Pinched cylinder: System sketch and comparison of deformations.