Kirchhoff-Love shell theory based on tangential differential calculus

The Kirchhoff-Love shell theory is recasted in the frame of the tangential differential calculus (TDC) where differential operators on surfaces are formulated based on global, three-dimensional coordinates. As a consequence, there is no need for a parametrization of the shell geometry implying curvilinear surface coordinates as used in the classical shell theory. Therefore, the proposed TDC-based formulation also applies to shell geometries which are zero-isosurfaces as in the level-set method where no parametrization is available in general. For the discretization, the TDC-based formulation may be used based on surface meshes implying element-wise parametrizations. Then, the results are equivalent to those obtained based on the classical theory. However, it may also be used in recent finite element approaches as the TraceFEM and CutFEM where shape functions are generated on a background mesh without any need for a parametrization. Numerical results presented herein are achieved with isogeometric analysis for classical and new benchmark tests. Higher-order convergence rates in the residual errors are achieved when the physical fields are sufficiently smooth.


Introduction
The mechanical modeling of shells leads to partial differential equations (PDEs) on manifolds where the manifolds are curved surfaces in the three-dimensional space. An overview in classical shell theory is given, e.g., in [4,9,32,44,45] or in the textbooks [1,5,41,49]. When modeling physical phenomena on curved surfaces, definitions for geometric quantities (normal vectors, curvatures, etc.) and differential surface operators (gradients, divergence, etc.) are key ingredients. These quantities may be either defined based on two-dimensional, curvilinear local coordinates living on the manifold or on global coordinates of the surrounding, three-dimensional space.
In the first case, the curved surface is parametrized by two parameters, i.e., there is a given map from the two-dimensional parameter space to the three-dimensional physical space, see Fig. 1(a). For the definition of geometrical quantities and surface operators, co-and contra-variant base vectors and Christoffel-symbols naturally occur. It is important to note that a parametrization of a surface is not unique, hence, there are infinitely many maps which result in the same curved surface. Obviously, the physical modeling must be independent of a concrete parametrization, which suggests the existence of a parametrization-free formulation.
In the second case, the geometric quantities and surface operators are based on global coordinates as done in the tangential differential calculus (TDC) [15,25,28]. Then, a model may also be defined even if a parametrization of a curved surface does not exist, for example, when it is a zero-isosurface of a scalar function in three dimensions following the level-set method [21,22,39,43]. When the physical modeling is based on the TDC, i.e., on global coordinates, it is applicable to surfaces which are parametrized or not. In this sense, the TDC-based approach is more general than approaches based on local coordinates. Models based on the TDC are found in various applications, see [16,17,18,22] for scalar problems such as heat flow and [20,31] for flow problems on manifolds. In the context of structure mechanics, this approach is used in [29] for curved beams, in [26,28,25] for membranes, and in [27] for flat shells embedded in R 3 .
Herein, we apply the TDC for the reformulation of the classical Kirchhoff-Love shell theory which is typically formulated based on a given parametrization. Based on the TDC, it is possible to also formulate the boundary value problem (BVP) for shell geometries where no parametrization is given as for the example in Fig. 1 given by the zero-isosurface of φ(x) = x − r with x ∈ [−r, r] 2 × [0, r] and the mechanical response to the force F is sought. As mentioned before, the TDC-based formulation is also valid when a parametrization is available; it is then equivalent to the classical formulation based on local coordinates.
Other attempts to parametrization-free formulations of the Kirchhoff-Love shell theory are found, e.g., in [11,12,13,14] with a mathematical focus and in [33,47,50] from an engineering perspective, however, only with focus on displacements. Herein, the Kirchhoff-Love shell theory is recasted in the frame of the TDC including all relevant mechanical aspects. For the first time, the parametrization-free strong form of the Kirchhoff-Love shell is given and taken as the starting point to derive the weak form. Then, boundary terms for the relevant boundary conditions of Kirchhoff-Love shell theory are naturally achieved. Furthermore, mechanical quantities such as moments, normal and shear forces are defined based on global coordinates and it is shown how (parametrization-)invariant quantities such as principal moments are computed. Finally, the strong form of Kirchhoff-Love shells is also found highly useful to define residual errors in the numerical results. Of course, evaluating this error in the strong form requires up to forth-order derivatives on the surface, which is implementationally quite some effort. The advantage, however, is that one may then confirm higher-order convergence rates in the corresponding error norm for suitable shell test cases. This is, otherwise, very difficult as exact solutions for shells are hardly available and classical benchmark tests typically give only selected scalar quantities, often with moderate accuracy.
For the numerical solution of shells, i.e., the approximation of the shell BVP based on numerical methods, we distinguish two fundamentally different approaches. The first is a classical finite element analysis based on a surface mesh, labelled Surface FEM herein [16,18,20,22]. Once a surface mesh is generated, it implies element-wise parametrizations for the shell geometry, see Fig. 1(c), no matter whether the underlying (analytic) geometry was parametrized or implied by level sets. In this case, classical shell theory based on parametrizations is suitable at least for the discretized geometry. The proposed TDC-based formulation is suitable as well which shall be seen in the numerical results. The other numerical approach is to use a three-dimensional background mesh into which the curved shell surface is embedded, cf. Fig. 1(d). Then, the shape functions of the (three-dimensional) background elements are only evaluated on the shell surface and no parametrization (and surface mesh) is needed to furnish basis functions for the approximation. For these methods, e.g., labelled CutFEM [6,7,8,19] or TraceFEM [23,37,38,42], applied to the case of shell mechanics, it is no longer possible to rely on classical parametrization-based formulations of the shell mechanics, however, the proposed TDC-based formulation is still applicable.
For the numerical results presented herein, the continuous weak form of the BVP is discretized with the Surface FEM [16,18,20,22] using NURBS as trial and test functions as proposed by Hughes et al. [10,30] due to the continuity requirements of Kirchhoff-Love shells. The boundary conditions are weakly enforced via Lagrange multiplies [51]. The situation is similar to [2,32,35,36], however, based on the proposed view point, the implementation is quite different. In particular, when PDEs on manifolds from other application fields than shell mechanics are also of interest (e.g., when transport problems [16,17,18] or flow problems [20,31] on curved surfaces are considered), there is a unified and elegant way to handle this by computing surface gradients applied to finite element shape functions which simplifies the situation considerably. In that sense one may shift significant parts of the implementation needed for shells to the underlying finite element technology and recycle this in other situations where PDEs on surfaces are considered.
We summarize the advantages of the TDC-based formulation of Kirchhoff-Love shells: (1) The definition of the BVP does not need a parametrization of the surface (though it can also handle the classical situation where a parametrization is given), (2) the TDC-based formulation is also suitable for very recent finite element technologies such as CutFEM and TraceFEM (though the typical approach based on the Surface FEM or IGA is also possible and demonstrated herein), (3) the implementation is advantagous in finite element (FE) codes where other PDEs on manifolds are considered as well due to the split of FE technology and application. From a didactic point of view, it may also be advantageous that troubles with curvilinear coordinates (co-and contra-variance, Christoffel-symbols) are avoided in the TDC-based approach where surface operators and geometric quantities are expressed in tensor notation.
The outline of the paper is as follows: In Section 2, important surface quantities are defined, and an introduction to the tangential differential calculus (TDC) is given. In Section 3, the classical linear Kirchhoff-Love shell equations under static loading are recast in terms of the TDC. Stress resultants such as membrane forces, bending moments, transverse shear forces and corner forces are defined. In Section 4, implementational aspects are considered. The element stiffness matrix and the resulting system of linear equations are shown. The implementation of boundary conditions based on Lagrange multipliers is outlined. Finally, in Section 5, numerical results are presented. The first example is a flat shell embedded in R 3 , where an analytical solution is available. The second and third example are parts of popular benchmarks as proposed in [2]. In the last example, a more general geometry without analytical solution or reference displacement is considered. The error is measured in the strong form of the equilibrium in order to verify the proposed approach and higherorder convergence rates are achieved.

Preliminaries
Shells are geometrical objects, where one dimension is significantly smaller compared to the other two dimensions. In this case, the shell can be reduced to a surface Γ embedded in the physical space R 3 . In particular, the surface is a manifold of codimension 1. Let the surface be possibly curved, sufficiently smooth, orientable, connected and bounded by ∂Γ. There are two alternatives for defining the shell geometry. One is through a parametrization, i.e., a (bijective) mapping from the parameter spaceΩ ⊂ R 2 to the real domain Γ ⊂ R 3 . The other approach is based on the level-set method. Then, a level-set function φ(x) : R 3 → R with x ∈ Ω ⊂ R 3 exists and the shell is implicitly given by Additional level-set functions may restrict the zero-isosurface to the desired, bounded shell as described in [22]. In Fig. 2(a) and (b) the two different approaches are schematically shown. The definition of the normal vector depends on whether the shell geometry is based on a parametrization or not. In the first case (cf. Fig. 2(a)), the shell geometry results from a map x(r). Then, the normal vector n Γ of the shell surface is determined by a cross-product of the columns of the Jacobi-matrix J(r) = ∂x /∂r. The resulting geometric quantities, surface operators, and models in this case are parametrization-based.
In the case where the shell geometry is implied by the zero-isosurface of a level-set function φ(x) (cf. Fig. 2(b)) and no parametrization is available, the normal vector may be determined by n Γ = ∇φ / ∇φ . All resulting quantities including the BVP of the Kirchhoff-Love shell are parametrization-free in this case. Of course, when in the wake of discretizing the BVP, the Surface FEM is used for the approximation, then a surface mesh of the shell geometry is needed and the surface elements do imply a parametrization again. It was already mentioned above, that other numerical methods such as the TraceFEM and Cut-FEM do not rely on a surface mesh. In this case, the countinuous and discrete BVP for the shell are truly parametrization-free.
In addition to the normal vector on the surface, along the boundary ∂Γ there is an associated tangential vector t ∂Γ ∈ R 3 pointing in the direction of ∂Γ and a co-normal vector n ∂Γ = n Γ × t ∂Γ ∈ R 3 pointing "outwards" and being perpendicular to the boundary yet in the tangent plane of the surface Γ. For the proof of equivalence of both cases we refer to, e.g., [18].

Tangential Differential Calculus
The TDC provides a framework to define differential operators avoiding the use of classical differential geometric methods based on local coordinate systems and Christoffel symbols. In the following, an overview of the operators and relations in the frame of the TDC are presented. For simplicity, we restrict ourselves to the case of surfaces embedded in the three dimensional space. However, the shown relations and definitions may be adopted to other situations accordingly (e.g., curved lines embedded in 2D or 3D). An introduction from a more mathematical point of view is given in [15,25,31].
Orthogonal projection operator P The orthogonal projection operator or normal projector P ∈ R 3×3 is defined as The operator ⊗ is the dyadic product of two vectors. The normal projector P projects a vector v onto the tangent space T P Γ of the surface. Note that P is idempotent (P·P = P), symmetric (P = P ) and obviously in the tangent space T P Γ of the surface, i.e., P · n Γ = n Γ · P = 0).
The projection of a vector field v : Γ → R 3 onto the tangent plane is defined by where v t is tangential, i.e. v t · n Γ = 0. The double projection of a second-order tensor function A(x) : Γ → R 3×3 leads to an in-plane tensor and is defined as with the properties A t = P · A t · P and A t · n Γ = n Γ · A t = 0.

Tangential gradient of scalar functions
The tangential gradient ∇ Γ of a scalar function u : Γ → R on the manifold is defined as where ∇ is the standard gradient operator in the physical space andũ is a smooth extension of u in a neighbourhood U of the manifold Γ. Alternatively,ũ is given as a function in global coordinatesũ(x) : R 3 → R and only evaluated at the manifoldũ | Γ = u.
For parametrized surfaces defined by the map x(r), and a given scalar function u(r) :Ω → R, the tangential gradient can be determined without explicitly computing an extensionũ using with J(r) = ∂x /∂r ∈ R 3×2 being the Jacobi-matrix, G = J · J is the metric tensor or the first fundamental form and the operator ∇ r is the gradient with respect to the reference coordinates. The components of the tangential gradient are denoted by representing first-order partial tangential derivatives. An important property of ∇ Γ u is that the tangential gradient of a scalar-valued function is in the tangent space of the surface ∇ Γ u ∈ T P Γ, i.e., ∇ Γ u · n Γ = 0. When using the Surface FEM to solve BVPs on surfaces, one may use Eq. 2.7 to compute tangential gradients of the shape functions. If, on the other hand, TraceFEM or CutFEM is used, one may use Eq. 2.6.

Tangential gradient of vector-valued functions
Consider a vector-valued function v(x) : Γ → R 3 and apply to each component of v the tangential gradient for scalars. This leads to the directional gradient of v defined as Note that the directional gradient is not in the tangent space of the surface, in general. A projection of the directional gradient to the tangent space leads to the covariant gradient of v and is defined as which is an in-plane tensor, i.e., ∇ cov Γ v ∈ T P Γ. The covariant gradient often appears in the modelling of physical phenomena on manifolds, i.e., in the governing equations. In contrast the directional gradient appears naturally in product rules or divergence theorems on manifolds.
In the following, partial surface derivatives of scalar functions are denoted as ∂ Γ x i u or u Γ ,i with i = 1, 2, 3. Partial surface derivatives of vector or tensor components are denoted as v dir i,j for directional and v cov i,j for covariant derivatives with i, j = 1, 2, 3.

Tangential gradient of tensor functions
For a second-order tensor function A(x) : Γ → R 3×3 , the partial directional gradient with respect to x i is defined as The directional gradient of the tensor function is then defined as The covariant partial derivative is determined by projecting the partial directional derivative onto the tangent space

Second-order tangential derivatives
Next, second-order derivatives of scalar functions are considered. The directional second order gradient of a scalar function u is defined by where He dir is the tangential Hessian matrix which is not symmetric in the case of curved manifolds [15], i.e., u dir ,ij = u dir ,ji . For the case of parametrized surfaces and a given scalar function in the reference space, the tangential Hessian-matrix can be determined by where, Q = J · G −1 , and Q ,r i denotes the partial tangential derivative of Q with respect to r i . The covariant counterpart is In contrast to He dir , He cov is symmetric and an in-plane tensor [48]. In the special case of flat surfaces embedded in R 3 the directional and covariant Hessian matrix are equal.

Tangential divergence operators
The divergence operator of a vector-valued function v(x) : Γ → R 3 is given as 17) and the divergence of a matrix or tensor function Note that div Γ A is, in general, not a tangential vector. It would only be tangential if the surface is flat and A is an in-plane tensor.

Weingarten map and curvature
The Weingarten map as introduced in [31,15] is defined as and is related to the second fundamental form in differential geometry. The Weingarten map is a symmetric, in-plane tensor and its two non-zero eigenvalues are associated with the principal curvatures The minus in Eq. 2.20 is due to fact that the Weingarten map is defined with the "outward" unit normal vector instead of the "inward" unit normal vector, which leads to positive curvatures of a sphere. The third eigenvalue is zero, because H is an in-plane tensor. The corresponding eigenvectors t 1 , t 2 and n Γ are perpendicular as H is symmetric. In Fig. 3, the osculating circles with the radii r i = 1 /κ i and the eigenvectors at a point P are shown.
The Gauß curvature is defined as the product of the principal curvatures K = 2 i=1 κ i and the mean curvature is introduced as κ = κ 1 + κ 2 = tr(H).

Divergence theorems in terms of tangential operators
The divergence theorem or Green's formula for a scalar function f ∈ C 1 (Γ) and a vector valued function v ∈ C 1 (Γ) 3 are defined as in [13,15] (2.21) The term with the mean curvature κ is vanishing if the vector v is tangential, then v · n Γ = 0. In extension to Eq. 2.21, Green's formula for second order tensor functions In the case of in-plane tensors, e.g., A t = P · A t · P, the term with the mean curvature κ vanishes due to A t · n Γ = 0 and we also have

The shell equations
In this section, we derive the linear Kirchhoff-Love shell theory in the frame of tangential operators based on a global Cartesian coordinate system. We restrict ourselves to infinitesimal deformations, which means that the reference and spatial configuration are indistinguishable. Furthermore, a linear elastic material governed by Hooke's law is assumed. As usual in the Kirchhoff-Love shell theory, the transverse shear strains and the change of curvature in the material law are neglected, which restricts the model to thin shells (tκ max 1).
With these assumptions, an analytical pre-integration with respect to the thickness leads to stress resultants such as normal forces and bending moments. The equilibrium in strong form is then expressed in terms of the stress resultants. Finally, the transverse shear forces may be identified via equilibrium considerations.

Kinematics
The middle surface Γ of the shell is a sufficiently smooth manifold embedded in the physical space R 3 . A point on the middle surface is denoted as x Γ ∈ Γ ⊂ R 3 and may be obtained explicitly or implicitly, see Section 2. With the unit-normal vector n Γ a point in the domain of the shell Ω of thickness t is defined by with ζ being the thickness parameter and |ζ| ≤ t /2. Alternatively, if the middle surface is defined implicitly with a signed distance function φ(x) the domain of the shell Ω is defined by In this case the middle surface Γ is the zero-isosurface of φ(x), see Eq. 2.2. The displacement field u Ω of a point P (x Γ , ζ) in the shell continuum Ω takes the form being the displacement field of the middle surface and w(x Γ ) being the difference vector, as illustrated in Fig. 4. Without transverse shear strains, the difference vector w expressed in terms of TDC is defined as in [13] w As readily seen in the equation above, the difference vector w is tangential. Alternatively, the difference vector w may also be re-written in terms of partial tangential derivatives of u and the normal vector n Γ Consequently, the displacement field of the shell continuum is only a function of the middle surface displacement u, the unit normal vector n Γ and the thickness parameter ζ.
The linearised, in-plane strain tensor ε Γ is defined by the symmetric part of the directional gradient of the displacement field u Ω , projected with P [26] (3.6) Finally, the whole strain tensor may be split into a membrane and bending part, as usual in the classical theory with Note that in the linearised bending strain tensor ε Γ,B , the term (∇ dir Γ u) · H is neglected as in classical theory [45,Remark 2.2] or [49]. The resulting membrane and bending strain in Eq. 3.7 are equivalent compared to the classical theory, e.g., [1]. In the case of flat shell structures as considered in [27] the membrane strain is only a function of the tangential displacement u t = P · u and the bending strain only depends on the normal displacement u n = u · n Γ , which simplifies the whole kinematic significantly. Moreover, the normal vector n Γ is then constant and the difference vector simplifies to w(x Γ ) = −∇ Γ u n .

Constitutive Equation
As already mentioned above, the shell is assumed to be linear elastic and, as usual for thin structures, plane stress is presumed. The in-plane stress tensor σ Γ is defined as where µ = E 2(1+ν) and λ = Eν (1−ν 2 ) are the Lamé constants and ε dir Γ is the directional strain tensor from Eq. 3.6. With this identity the in-plane stress tensor can be computed only with the directional strain tensor which is from an implementational point of view an advantage, because covariant derivatives are not needed explicitly. In comparison to the classical theory, the in-plane stress tensor expressed in terms of TDC does not require the computation of the metric coefficients in the material law. Therefore, the resulting stress tensor does not hinge on a parametrization of the middle surface and shell analysis on implicitly defined surfaces is enabled.

Stress resultants
The stress tensor is only a function of the middle surface displacement vector u, the difference vector w(u) and the thickness parameter ζ. This enables an analytical preintegration with respect to the thickness and stress resultants can be identified. The following quantities are equivalent to the stress resultants in the classical theory [45,1], but they are expressed in terms of the TDC using a global Cartesian coordinate system.
The moment tensor m Γ is defined as with is the flexural rigidity of the shell. The moment tensor m Γ is symmetric and an in-plane tensor. Therefore, one of the three eigenvalues is zero and the two non-zero eigenvalues of m Γ are the principal bending moments m 1 and m 2 . The principal moments are in agreement with the eigenvalues of the moment tensor in the classical setting, see [1]. For the effective normal force tensorñ Γ we havẽ with Similar to the moment tensor, the two non-zero eigenvalues ofñ Γ are in agreement with the effective normal force tensor expressed in local coordinates. Note that for curved shells this tensor is not the physical normal force tensor. This tensor only appears in the variational formulation, see Section 4. The physical normal force tensor n real Γ is defined by and is, in general, not symmetric and also has one zero eigenvalue. The occurrence of the zero eigenvalues in m Γ ,ñ Γ and n real Γ is due to fact that these tensors are in-plane tensors, i.e. m Γ · n Γ = n Γ · m Γ = 0 . The normal vector n Γ is the corresponding eigenvector to the zero eigenvalue and the other two eigenvectors are tangential.

Equilibrium
Based on the stress resultants from above, one obtains the equilibrium for a curved shell in strong form as with f being the load vector per area on the middle surface Γ. A summation over the indices i, k = 1, 2, 3 has to be performed. The obtained equilibrium does not rely on a parametrization of the middle surface but is, otherwise, equivalent to the equilibrium in local coordinates [1,49]. From this point of view, the reformulation of the linear Kirchhoff-Love shell equations in terms of the TDC may be seen as a generalization, because the requirement of a parametrized middle surface is circumvented. With boundary conditions, as shown in detail in Section 3.3.2, the complete fourth-order boundary value problem (BVP) is defined.
Based on the equilibrium in Eq. 3.13, the transverse shear force vector q is defined as (3.14) Note that in the special case of flat Kirchhoff-Love structures embedded in R 3 the divergence of an in-plane tensor is a tangential vector, as already mentioned in Section 2.1. Therefore, the definition of the transverse shear force vector in [27] is in agreement with the obtained transverse shear force vector herein.

Equilibrium in weak form
The equilibrium in strong form is converted to a weak form by multiplying Eq. 3.13 with a suitable test function v and integrating over the domain, leading to With Green's formula from Section 2.1, we introduce the continuous weak form of the equilibrium: The corresponding function spaces are where ∂Γ D is the Dirichlet boundary and ∂Γ N is the Neumann boundary. The advantage of this procedure is that the boundary terms naturally occur and directly allow to consider for mechanically meaningful boundary conditions.

Boundary conditions
As well known in the classical Kirchhoff-Love shell theory, special attention needs to be paid to the boundary conditions. In the following, the boundary terms of the weak form in Eq. 3.16 are rearranged in order to derive the effective boundary forces.
Using Eqs. (3.12) and (3.4), we have As already mentioned above, the difference vector w is a tangential vector. Consequently, the difference vector at the boundary may be expressed in terms of the tangential vectors t ∂Γ and n ∂Γ where ω t ∂Γ = ω t ∂Γ t ∂Γ may be interpreted as rotation along the boundary and ω n ∂Γ = ω n ∂Γ n ∂Γ is the rotation in co-normal direction, when the test function v is interpreted as a displacement, see Fig. 5(a). Analogously to the difference vector, the expressions n real Γ · n ∂Γ and m Γ · n ∂Γ in Eq. 3.19 are decomposed in a similar manner Next, the term p n Γ = P · div Γ m Γ · n ∂Γ represents the resultant force in normal direction. In Fig. 5(b) the forces and bending moments along a curved boundary are illustrated.
(b) forces and moments Fig. 5: Decomposition of the difference vector w, in-plane normal forces n real Γ · n ∂Γ and bending moments m Γ · n ∂Γ along the boundary ∂Γ in terms of t ∂Γ and n ∂Γ : (a) Rotations at the boundary, (b) normal force tensor and bending moments at the boundary.
Inserting these expressions in Eq. 3.19, the integral along the Neumann boundary simplifies to As discussed in detail, e.g., in [1], the rotation in co-normal direction ω n ∂Γ is already prescribed with v| ∂Γ . Therefore, the term ω n ∂Γ m n ∂Γ is expanded and with integration by parts we obtain where +C and −C are points close at a corner C. The new boundary term are the Kirchhoff forces or corner forces. Note that if the boundary of the shell is smooth, the corner forces vanish. Finally, the integral over the Neumann boundary in Eq. 3.16 is expressed in terms of the well-known effective boundary forces and the bending moment along the boundary The obtained effective boundary forces and moments are in agreement with the given quantities in local coordinates [1,49]. The prescribeable boundary conditions are the conjugated displacements and rotations to the effective forces and moments at the boundarỹ In Tab. 1, common support types are given. Other boundary conditions (e.g., membrane support, . . . ) may be derived, with the quantities above, accordingly.

Implementational aspects
The continuous weak form is discretized using isogeometric analysis as proposed by Hughes et al. [30,10]. The NURBS patch T is the middle surface of the shell and the elements τ i (i = 1, . . . , n Elem ) are defined by the knot spans of the patch. The mesh is then defined by the union of the elements Γ = τ ∈T τ .
There is a fixed set of local basis functions {N k i (r)} of order k with i = 1, . . . , n k being the number of control points and the displacements {û i ,v i ,ŵ i } stored at the control points i are the degrees of freedom. Using the isoparametric concept, the shape functions N k i (r) are NURBS of order k. The surface derivatives of the shape functions are computed as defined in Section 2 , similar as in the Surface FEM [18,16,20,22] using NURBS instead of Lagrange polynomials as ansatz and test functions. The shape functions of order k ≥ 2 are in the function space V, see Eq. 3.17. In fact, the used shape functions are in the The resulting element stiffness matrix K Elem is a 3 × 3 block matrix and is divided into a membrane and bending part (4.1) The membrane part is defined by summation over a and b. The matrixK is determined by directional first-order derivatives of the shape functions N . One may recognize that the structure of the matrixK is similar to the stiffness matrix of 3D linear elasticity problems. For the bending part we have A summation over a, b on the one hand and c, d on the other has to be performed. The first term ofK is the contraction of the covariant Hessian matrix He cov Γ and the second term may be identified as the Bi-Laplace operator. Note that for the Bi-Laplace operator also Implementational aspects directional derivatives may be used, due to the fact that the trace of second order derivatives is invariant, although the components differ. This suggests a further rearrangement of the contraction of the covariant Hessian matrix in order to use only directional derivatives, which is preferred from an implementational point of view. The equivalent expression of K using only second-order directional derivatives is with summation over a, b, e and c, d as above. When the shell is given through a parametrization, the resulting element stiffness matrix in the classical theory, e.g., [9] is equivalent to the element stiffness matrix from above, but in the classical setting the computation may be found more cumbersome due to fact that the local basis vectors and the metric tensor in co-and contra-variant form has to be computed. In contrast, herein, the surface gradients and second-order derivatives are first applied to the shape functions (NURBS or classical finite element functions) to obtain N Γ ,i , N dir ,ij and N cov ,ij , which is independent of the application.
In this sense, a significant part of the complexity of implementing shells is shifted to finite element technology and may be recycled for any kind of boundary value problems on curved surfaces in R 3 . Examples are transport problems [17,16,18] or flow problems [20,31] on curved surfaces. We expect that future implementations in finite element software will provide frameworks for solving PDEs on manifolds and, based, e.g., on this work will also apply to shells. In order to emphasize the differences in the implementation, example Matlab®codes for the proposed TDC-based formulation and the classical parametrization-based formulation are given in Appendix A, clearly highlighting the differences.
The boundary conditions are weakly enforced by Lagrange multipliers [51]. The shape functions of the Lagrange multipliers are NURBS of the same order than the shape functions of the displacements. Therefore, the shape functions of the Lagrange multiplier i is defined as For the test cases shown in Section 5, bounded condition numbers and unique solutions are observed. The usual assembly yields a linear system of equations in the form with [û,λ] = [û,v,ŵ,λ] being the sought displacements of the control points and Lagrange multipliers. With the shape functions of the Lagrange multipliers N L , the constraint matrix C for simply supported edges is defined by for clamped edges and for symmetry supports Note that all constraint matrices have three block-columns refering to the unknownŝ u,v,ŵ.

Numerical results
The numerical results are achieved using NURBS functions for the geometry and shape function definition, following the methodology of isogeometric analysis (IGA) [30,32,3,34,24]. The definition of NURBS is omitted here for brevity but is found at numerous references in the frame of IGA, e.g., [40,10].
The obtained shell equations are carefully verified and compared to the classical approach with different test cases. As already mentioned above the proposed approach leads to an equivalent stiffness matrix for arbitrary curved and non-curved shells. Consequently, the same convergence properties as shown, e.g., in [32,9] are expected. In the following, the results of the convergence analyses of a flat shell embedded in R 3 , the Scordelis-Lo roof, and the pinched cylinder test (part of the shell obstacle course proposed by Belytschko et al. [2]) are shown. Furthermore, a new test case with a challenging geometry is proposed which features smooth solutions enabling higher-order convergence rates. These rates are confirmed in the residual error as no analytic solution exists, see Section 5.4. Other examples (e.g., pinched hemispherical shell, shells of revolution, etc.) have been considered but are omitted here for brevity.
In the convergence studies, NURBS patches with different orders and numbers of knot spans in each direction are employed. This is equivalent to meshes with higher-order elements and n =

Flat shell embedded in R 3
Following a similar rationale as in [27], as a first test case, we consider a simple quadrilateral, flat shell with the normal vector n Fig. 6. The shell is simply supported at all edges. For verification, the load vector f is split into tangential f t and normal f n loads. The tangential loads are obtained with the method of manufactured solution for a given displacement field u t (x) = [ [1, 1] · 1 /4 · sin(πr) sin(πs)] • χ −1 . In normal direction, a sinusoidal load f n (x) = [−D sin(πr) sin(πs)] • χ −1 is applied to the shell. Herein, χ is an affine mapping function (rigid-body rotation) from the horizontal parameter space to the real domain. An analytic solution for the normal displacements is easily obtained with u n (x) = [−(sin(πr) sin(πs))/(4π 4 )] • χ −1 , [46]. The shell is defined with L = 1 and the thickness is set to t = 0.01. The material parameters are: Young's modulus E = 10 000 and the Poisson's ratio ν = 0.3.
In Fig. 7, the solution of the shell is illustrated. The displacements are scaled by two orders of magnitude. The colours on the deformed surface indicate the Euclidean norm of the displacement field u .
The results of the convergence analysis are shown in Fig. 8. The curves are plotted as a function of the element size 1 /n (which is rather a characteristic length of the knot spans). The dotted lines indicate the theoretical optimal order of convergence. In Fig. 8(a), the relative L 2 -error of the primal variable (displacements) is shown. Optimal higher-order convergence rates O(p + 1) are achieved. In the figures Fig. 8(b) to Fig. 8(d), the relative L 2 -errors of the normal forces (membrane forces), bending moments and transverse shear forces are plotted. For all stress resultants the theoretical optimal orders of convergence are achieved. It is clear that the same results were obtained if the results are computed for the purely two-dimensional case as, e.g., in [9].

Scordelis-Lo roof
The Scordelis-Lo roof is a cylindrical shell and is supported with two rigid diaphragms at the ends. The shell is loaded by gravity forces, see In Fig. 10(a), the numerical solution of the Scordelis-Lo roof is illustrated. The displacements are magnified by one order of magnitude.
In Fig. 10(b), the convergence of the maximum displacement u z,max is plotted up to polynomial order of p = 6 as a function of the element (knot span) size. It is clearly seen that the expected results are achieved, with increasing accuracy for higher-order NURBS. Due to the lack of a more accurate reference solutions, it is not useful to show these results in a double-logarithmic diagram as usual for error plots. The style of presentation follows those of many other references such as, e.g., in [2,9,32].

Pinched cylinder
The next test case is a cylindrical shell pinched with two diametrically opposite unit loads located within the middle of the shell, see Fig. 11. The cylinder is defined with L = 600, R = 300. The thickness is set to t = 3. The material properties are: Young's modulus E = 3 × 10 6 and the Poisson's ratio ν = 0.3. The reference displacement at the loading points are u Ref = 1.824 88 × 10 −5 as given in reference [2]. Due to symmetry only one eighth of the geometry is modelled.
In Fig. 12(a), the numerical solution of the pinched cylinder is illustrated with scaled displacements by a factor of 5 × 10 6 .
As in the example before, in Fig. 12(b), the convergence to a normalized reference displacement as a function of the element size is plotted. The results confirm with the expected convergence behaviour as shown in [9,32]. It is noted that due to the singularity in some mechanical quantities due to the single force, higher-order convergence rates are not possible here. However, the improvement for increasing the order of the NURBS is still seen in the figure. An additional grading of the elements in order to better resolve the singularity would have further improved the situation but is omitted here.

Flower shaped shell
As a last example, a more complex geometry is considered, which enables smooth mechanical fields and thereby enables higher-order convergence rates. The geometry of the middle surface is given with and illustrated in Fig. 13. On the right side of the figure, the boundary conditions and material parameters are defined. The middle surface of the shell features varying principal curvatures and curved boundaries. The curved boundaries are clamped and the corresponding conditions (from Tab. 1) have to be properly enforced. An analytical solution or reference displacement is not available. Therefore, the error is measured in the strong form of the equilibrium from Eq. 3.13 and may be called residual error. In particular, the residual error is the summed element-wise relative L 2 -error The computation of the residual error requires the evaluation of fourth-order surface derivatives. It is noteworthy that the implementation of these higher-order derivatives is not without efforts. For example, recall that mixed directional surface derivatives are not symmetric. That is, there are 3 4 = 81 partial fourth-order derivatives. Nevertheless, if the displacement field is smooth enough this error measure is a suitable quantity for the convergence analysis.
In Fig. 14(a), the deformed shell is illustrated. The displacement field is scaled by one order of magnitude. In Fig. 14(b), the results of the convergence analysis are plotted. Due to the fact that fourth-order derivatives need to be computed, at least fourth-order shape functions are required. The theoretical optimal order of convergence is O(p − 3) if the solution is smooth enough. One may observe that higher-order convergence rates are achieved, however, rounding-off errors and the conditioning may slightly influence the convergence. Nevertheless, the results are excellent also given the fact that higher-order accurate results for shells (given in double-logarithmic error plots) are the exception.
The stored elastic energy at the finest level with a polynomial order p = 8, which may be seen as an overkill solution, is e = 1.7635958 ± 1 × 10 −7 kN m. This stored elastic energy may be used for future benchmark tests, without the need to implement fourth-order

Conclusions and Outlook
The linear Kirchhoff-Love shell theory is reformulated in terms of the TDC using a global Cartesian coordinate system and tensor notation. The resulting model equations apply to shell geometries which are parametrized or not. For example, a parametrization may not be available when shell geometries are implied by the level-set method. Because the TDC-based formulation holds in both cases, it may be seen as a generalization to the classical shell theory which is based on parametrizations and curvilinear coordinates.
The TDC-based strong form is used as the starting point to consistently obtain the weak form including all boundary terms well-known in the Kirchhoff-Love theory. Mechanical stress-resultants such as moments, normal and shear forces are defined in global coordinates. Furthermore, the strong form may be used in the numerical results to compute residual errors and thus enable convergence analyses even without the knowledge of exact solutions which, for shells, are scarce.
For the discretization, the Surface FEM is used with NURBS as trial and test functions. That is, an isogeometric approach is demonstrated due to continuity requirements. In this case, the presence of a surface mesh (i.e., a NURBS patch), implies a parametrization and although the involved equations and the resulting implementations vary significantly, it is seen that the classical, parametrization-based and the proposed TDC-based formulation are equivalent. For a generic finite element framework enabling various implementations for PDEs on manifolds (in addition to only shells), the TDC-based approach is benefitial, because surface gradients of shape functions may be computed beforehand and are independent of the application.
The numerical results confirm higher-order convergence rates. As mentioned, based on the residual errors, a framework for the verification of complex test cases is presented. There is a large potential in the parametrization-free reformulation of shell models, because the obtained PDEs may be discretized with new finite element techniques such as TraceFEM or CutFEM based on implicitly defined surfaces. In this case, neither the problem statement nor the discretization is based on a parametrization. % Add contribution at integration point to element matrix and rhs. ElemMat = ElemMat + ipRef.ww(i)*A3bar * (KKmemb + KKbend); ElemRhs = ElemRhs + ipRef.ww(i)*A3bar * Rhs; end function [fx, fy, fz] = EvaluateLoad(xx, yy, zz) nn = length(xx); fx = zeros(nn, 1); fy = zeros(nn, 1); fz = zeros(nn, 1); end 38 Element stiffness matrix Code 2: Element contribution in terms of local coordinates