Kirchhoff-Love shell theory based on tangential differential calculus

by   D. Schöllhammer, et al.
TU Graz

The Kirchhoff-Love shell theory is recasted in the frame of the tangential differential calculus where differential operators on surfaces are formulated without the need for a parametrization, i.e., local coordinates. The governing equations are presented in strong and weak form including a detailed discussion of the boundary conditions and mechanical quantities such as moments, normal and shear forces. Although the parameter-free formulation leads to identical results than the classical formulation, it is more general as it applies also to shells whose geometry is implied by level-sets and, hence, no parametrization is available. Furthermore, it enables a different and unified view point on shell mechanics in general and leads to an elegant implementation. Numerical results are achieved based on 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.



page 32

page 35


Reissner-Mindlin shell theory based on tangential differential calculus

The linear Reissner-Mindlin shell theory is reformulated in the frame of...

A Consistent Higher-Order Isogeometric Shell Formulation

Shell analysis is a well-established field, but achieving optimal higher...

A general higher-order shell theory for compressible isotropic hyperelastic materials using orthonormal moving frame

The aim of this study is three-fold: (i) to present a general higher-ord...

A Unified Finite Strain Theory for Membranes and Ropes

The finite strain theory is reformulated in the frame of the Tangential ...

Geometric Formulation for Discrete Points and its Applications

We introduce a novel formulation for geometry on discrete points. It is ...

A C^1-continuous Trace-Finite-Cell-Method for linear thin shell analysis on implicitly defined surfaces

A Trace-Finite-Cell-Method for the numerical analysis of thin shells is ...

Optimization and Learning With Nonlocal Calculus

Nonlocal models have recently had a major impact in nonlinear continuum ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 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. (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 .

((a)) parametric via map
((b)) Shell-BVP on zero-isosurface
((c)) maps implied by FEM
((d)) TraceFEM, CutFEM
Fig. 1: (a) In classical shell mechanics, the middle surface is defined by a parametrization, i.e., a map . (b) The cupola is given by the zero-isosurface of and the mechanical response to the force is sought. (c) The surface mesh implies element-wise, approximate parametrizations even if the initial geometry is defined by level-sets. (d) For implicitly defined shells in the context of TraceFEM and CutFEM, no parametrization is needed at all.

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. (b): The cupola with radius r is given by the zero-isosurface of with and the mechanical response to the force 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. (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. (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 , 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 higher-order convergence rates are achieved.

2 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 . 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 to the real domain . The other approach is based on the level-set method. Then, a level-set function with 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. (a) and (b) the two different approaches are schematically shown.

((a)) explicit
((b)) implicit
Fig. 2: Examples of bounded surfaces embedded in the physical space : (a) Explicitly defined surface with a map , (b) Implicitly defined surface with a master level-function (yellow) and slave level-set functions for the boundary definition (gray).

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. (a)), the shell geometry results from a map . Then, the normal vector of the shell surface is determined by a cross-product of the columns of the Jacobi-matrix . 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 (cf. Fig. (b)) and no parametrization is available, the normal vector may be determined by . 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 CutFEM 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 pointing in the direction of and a co-normal vector 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].

2.1 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 D or D). An introduction from a more mathematical point of view is given in [15, 25, 31].

Orthogonal projection operator

The orthogonal projection operator or normal projector is defined as


The operator is the dyadic product of two vectors. The normal projector projects a vector onto the tangent space of the surface. Note that is idempotent , symmetric and obviously in the tangent space of the surface, i.e., .

The projection of a vector field onto the tangent plane is defined by


where is tangential, i.e. . The double projection of a second-order tensor function leads to an in-plane tensor and is defined as


with the properties and .

Tangential gradient of scalar functions

The tangential gradient of a scalar function on the manifold is defined as


where is the standard gradient operator in the physical space and is a smooth extension of in a neighbourhood of the manifold . Alternatively, is given as a function in global coordinates and only evaluated at the manifold .

For parametrized surfaces defined by the map , and a given scalar function , the tangential gradient can be determined without explicitly computing an extension using


with being the Jacobi-matrix, is the metric tensor or the first fundamental form and the operator 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 is that the tangential gradient of a scalar-valued function is in the tangent space of the surface , i.e., . 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 and apply to each component of the tangential gradient for scalars. This leads to the directional gradient of 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 and is defined as


which is an in-plane tensor, i.e., . 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 or with . Partial surface derivatives of vector or tensor components are denoted as for directional and for covariant derivatives with .

Tangential gradient of tensor functions

For a second-order tensor function , the partial directional gradient with respect to 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 is defined by


where is the tangential Hessian matrix which is not symmetric in the case of curved manifolds [15], i.e., . For the case of parametrized surfaces and a given scalar function in the reference space, the tangential Hessian-matrix can be determined by


where, , and denotes the partial tangential derivative of with respect to . The covariant counterpart is


In contrast to , is symmetric and an in-plane tensor [48]. In the special case of flat surfaces embedded in the directional and covariant Hessian matrix are equal.

Tangential divergence operators

The divergence operator of a vector-valued function is given as


and the divergence of a matrix or tensor function , is


Note that is, in general, not a tangential vector. It would only be tangential if the surface is flat and 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

is an in-plane tensor. The corresponding eigenvectors

and are perpendicular as is symmetric. In Fig. 3, the osculating circles with the radii and the eigenvectors at a point are shown.

Fig. 3: Osculating circles (blue, red) and eigenvectors of at point on a surface embedded in .

The Gauß curvature is defined as the product of the principal curvatures and the mean curvature is introduced as .

Divergence theorems in terms of tangential operators

The divergence theorem or Green’s formula for a scalar function and a vector valued function are defined as in [13, 15]


The term with the mean curvature is vanishing if the vector is tangential, then . In extension to Eq. 2.21, Green’s formula for second order tensor functions , is


where . In the case of in-plane tensors, e.g., , the term with the mean curvature vanishes due to and we also have .

3 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 .

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.

3.1 Kinematics

The middle surface of the shell is a sufficiently smooth manifold embedded in the physical space . A point on the middle surface is denoted as and may be obtained explicitly or implicitly, see Section 2. With the unit-normal vector a point in the domain of the shell of thickness is defined by


with being the thickness parameter and . Alternatively, if the middle surface is defined implicitly with a signed distance function the domain of the shell is defined by


In this case the middle surface is the zero-isosurface of , see Eq. 2.2. The displacement field of a point in the shell continuum takes the form


with being the displacement field of the middle surface and being the difference vector, as illustrated in Fig. 4.

Fig. 4: Displacements and of the shell.

Without transverse shear strains, the difference vector expressed in terms of TDC is defined as in [13]


As readily seen in the equation above, the difference vector is tangential. Alternatively, the difference vector may also be re-written in terms of partial tangential derivatives of and the normal vector


Consequently, the displacement field of the shell continuum is only a function of the middle surface displacement , the unit normal vector and the thickness parameter .

The linearised, in-plane strain tensor is defined by the symmetric part of the directional gradient of the displacement field , projected with [26]


Finally, the whole strain tensor may be split into a membrane and bending part, as usual in the classical theory


Note that in the linearised bending strain tensor , the term 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 and the bending strain only depends on the normal displacement , which simplifies the whole kinematic significantly. Moreover, the normal vector is then constant and the difference vector simplifies to .

3.2 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 and are the Lamé constants and 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.

3.2.1 Stress resultants

The stress tensor is only a function of the middle surface displacement vector , the difference vector and the thickness parameter . This enables an analytical pre-integration 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 is defined as


where is the flexural rigidity of the shell. The moment tensor is symmetric and an in-plane tensor. Therefore, one of the three eigenvalues is zero and the two non-zero eigenvalues of are the principal bending moments and . 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 have


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 is defined by


and is, in general, not symmetric and also has one zero eigenvalue. The occurrence of the zero eigenvalues in , and is due to fact that these tensors are in-plane tensors, i.e.  . The normal vector is the corresponding eigenvector to the zero eigenvalue and the other two eigenvectors are tangential.

3.3 Equilibrium

Based on the stress resultants from above, one obtains the equilibrium for a curved shell in strong form as


with being the load vector per area on the middle surface . A summation over the indices 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 is defined as


Note that in the special case of flat Kirchhoff-Love structures embedded in 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.

3.3.1 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 and integrating over the domain, leading to


With Green’s formula from Section 2.1, we introduce the continuous weak form of the equilibrium:

Find such that


The corresponding function spaces are


where is the Dirichlet boundary and 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.

3.3.2 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 is a tangential vector. Consequently, the difference vector at the boundary may be expressed in terms of the tangential vectors and


where may be interpreted as rotation along the boundary and is the rotation in co-normal direction, when the test function is interpreted as a displacement, see Fig. (a). Analogously to the difference vector, the expressions and in Eq. 3.19 are decomposed in a similar manner


Next, the term represents the resultant force in normal direction. In Fig. (b) the forces and bending moments along a curved boundary are illustrated.

((a)) rotations
((b)) forces and moments
Fig. 5: Decomposition of the difference vector , in-plane normal forces and bending moments along the boundary in terms of and : (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 is already prescribed with . Therefore, the term is expanded and with integration by parts we obtain


where and are points close at a corner . 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 boundary

In Tab. 1, common support types are given. Other boundary conditions (e.g., membrane support, …) may be derived, with the quantities above, accordingly.

Clamped edge
Simply supported edge
Symmetry support
Free edge
Tab. 1: Set of common boundary conditions

4 Implementational aspects

The continuous weak form is discretized using isogeometric analysis as proposed by Hughes et al. [30, 10]. The NURBS patch is the middle surface of the shell and the elements are defined by the knot spans of the patch. The mesh is then defined by the union of the elements

There is a fixed set of local basis functions of order with being the number of control points and the displacements stored at the control points

are the degrees of freedom. Using the isoparametric concept, the shape functions

are NURBS of order . 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 are in the function space , see Eq. 3.17. In fact, the used shape functions are in the Sobolev space iff .

The resulting element stiffness matrix is a block matrix and is divided into a membrane and bending part


The membrane part is defined by


summation over and . The matrix is determined by directional first-order derivatives of the shape functions . One may recognize that the structure of the matrix is similar to the stiffness matrix of D linear elasticity problems. For the bending part we have