Reissner-Mindlin shell theory based on tangential differential calculus

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

The linear Reissner-Mindlin shell theory is reformulated in the frame of the tangential differential calculus (TDC) using a global Cartesian coordinate system. The rotation of the normal vector is modelled with a difference vector approach. The resulting equations are applicable to both explicitly and implicitly defined shells, because the employed surface operators do not necessarily rely on a parametrization. Hence, shell analysis on surfaces implied by level-set functions is enabled, but also the classical case of parametrized surfaces is captured. As a consequence, the proposed TDC-based formulation is more general and may also be used in recent finite element approaches such as the TraceFEM and CutFEM where a parametrization of the middle surface is not required. Herein, the numerical results are obtained by isogeometric analysis using NURBS as trial and test functions for classical and new benchmark tests. In the residual errors, optimal higher-order convergence rates are confirmed when the involved physical fields are sufficiently smooth.



page 1

page 2

page 3

page 4


Kirchhoff-Love shell theory based on tangential differential calculus

The Kirchhoff-Love shell theory is recasted in the frame of the tangenti...

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

A Consistent Higher-Order Isogeometric Shell Formulation

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

Numerical solution of the div-curl problem by finite element exterior calculus

We are interested in the numerical reconstruction of a vector field with...

Subdivision surfaces with isogeometric analysis adapted refinement weights

Subdivision surfaces provide an elegant isogeometric analysis framework ...

A novel four-field mixed variational approach to Kirchhoff rods implemented with finite element exterior calculus

A four-field mixed variational principle is proposed for large deformati...
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

Shells are curved, thin-walled structures that occur in a wide range of applications in both nature and technology. Because they often feature a high bearing capacity, they are frequently used in engineering applications such as automotive, aerospace, biomedical- and civil-engineering, see e.g., [12, 30, 65]. For the mechanical modelling, one may distinguish two classes of models based on the fact whether traverse shear deformations are considered or not. The Kirchhoff-Love shell theory [45, 47]

neglects such shear deformation resulting in a fourth-order partial differential equation (PDE) for the displacement of the middle surface of thin shells. For the numerical treatment, it is crucial to consider continuity requirements due to the variational index 2

[28]. When taking shear deformations into account, the Reissner-Mindlin shell theory [54] is typically employed and models the deformation of thin and moderately thick shells. The resulting shell equations are a system of second-order PDEs with the unknowns being the displacement of the middle surface and the rotation of the normal vector. An advantage of this model is that the corresponding variational index is and only -continuity in a finite element analysis is required. On the other hand, the model often suffers from locking phenomena for increasingly thin shells which may require further measures in the numerical treatment. An overview of classical shell theory is given, e.g., in [15, 58, 59, 7] or in the text books [12, 8, 2, 62, 65].

The middle surface of a shell is a manifold of codimension 1 embedded in the physical space . One may distinguish two alternatives for the definition of the shells: (1) explicit or (2) implicit. The first approach is typically used in the classical shell theory, where the middle surface of the shell is given through a parametrization, i.e., a map from a two-dimensional parameter space to the real surface embedded in the physical space . Based on this parametrization, curvilinear coordinates may be introduced and geometric quantities (normal vectors, curvatures, etc.) and differential surface operators (gradients, divergence, etc.) are defined. These quantities are the key ingredients of the PDEs obtained in the classical shell theory. There are various implementations of the classical shell models using this approach, a general overview is given, e.g., in [63], Kirchhoff-Love shells may be found, e.g., in [17, 16, 43, 48, 60], Reissner-Mindlin shells, e.g., in [5, 25, 24, 3, 46, 44] and hierarchical shells, e.g., in [28, 49, 6].

The alternative to parametrizations is when the middle surface is defined implicitly, e.g. as the zero-isosurface of a scalar function in three dimensions following the level-set method [33, 34, 52, 57]. The geometric quantities and differential surface operators are defined in the global Cartesian coordinates system of the surrounding physical space as done in the tangential differential calculus (TDC) [22, 36, 39, 56]. When the physical modeling process is based on the TDC, it is applicable to surfaces which are parametrized (explicit definition) or not (implicit definition). In this sense, the TDC-based approach is more general than approaches based on local coordinates, which are restricted to explicit surface descriptions. Models based on the TDC are found in various applications, see [23, 26, 27, 34] for scalar problems such as heat flow and [32, 42] for flow problems on manifolds. In the context of structure mechanics, this approach is used in [19, 20, 21, 61, 38, 56] for Kirchhoff-Love shells, in [40] for curved beams, and in [37, 39, 36] for membranes. In addition to a more general formulation of the PDEs in the frame of the TDC, there is unified and elegant implementation in the finite element context that may be recycled in other situations where PDEs on surfaces are considered. More precisely, typical surface operators in the TDC, in particular, the surface gradient may be applied to the finite element shape functions independently of the concrete application. This enables a shift of significant parts of the implementation needed for shells to the underlying finite element technology.

In [56]

, the authors proposed a reformulation of the Kirchhoff-Love shell based on the TDC. Herein, the aim is to recast the Reissner-Mindlin shell in the TDC-framework using a difference vector formulation. All relevant mechanical fields such as membrane forces, bending moments and shear forces are expressed in the global Cartesian coordinate system using the TDC and the computation of invariant quantities such as principal moments is shown. Furthermore, a parametrization-free strong form of the force and moment equilibriums is derived and used as a starting point for the derivation of the weak form. Similar to

[56], the concept of residual errors is also used to 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 scarce and standard benchmarks hardly serve for higher-order convergence studies.

For the numerical solution of the shell boundary value problem (BVP), one may employ two fundamentally different approaches of solving PDEs on manifolds. The first approach is a classical finite element method, labelled Surface FEM herein [23, 27, 32, 34]. In the classical Surface FEM the analytical geometry is decomposed into a set of elements. Each element of the discrete surface implies an element-wise parametrization, no matter if the geometry was originally defined explicitly or implicitly. In this case, the classical shell theory based on local coordinates and the TDC-based formulation are suitable choices provided that the geometry is discretized by finite elements. The alternative numerical approach is to embed the implicitly defined shell surface in a three-dimensional background mesh. Boundaries of the shell can be defined by means of the boundary of the background mesh or with additional level-set functions as in [33, 34]. The trial and test functions are the three-dimensional shape functions of the background mesh restricted to the trace of the shell surface. These methods are labelled as CutFEM [10, 11, 13, 29] or TraceFEM [35, 50, 51, 55]. When applying these methods to shell mechanics, it is crucial to use a TDC-based formulation such as proposed in this work, because no parametrization is available nor needed.

The presented numerical results herein are obtained with the Surface FEM [23, 27, 32, 34] using NURBS as trial and test functions as proposed by Hughes et al. [18, 41] in a general sense. In [5, 25, 24, 44], isogeometric analysis (IGA) is also applied to the Reissner-Mindlin shell, however, based on the classical, parameter-based formulation of the governing equations. It is important to note that these results coincide with the proposed formulation, however, the implementation differs significantly. Although IGA is used in the numerical results shown here, we emphasize that the presented TDC-based formulation, being the main aspect of this work, is more general as it may also be used when no parametrization is available such as in CutFEM or TraceFEM.

It is pointed out that continuity requirements would also allow a standard FEM implementation using -continuous shape functions. Nevertheless, we prefer the use of NURBS here, for example, due to the continuous normal vector and improved convergence properties. As mentioned above, a standard difference vector formulation is employed here and different strategies of discretizing tangential vector fields are outlined. We use Lagrange multipliers in order to weakly enforce (i) the tangentiality constraint on the globally defined difference vector and (ii) the boundary conditions. The proposed approach is in the case of parametrized surfaces equivalent to the classical Reissner-Mindlin shell. Therefore, locking phenomena can be expected in the case of thin shells. As shown in the numerical results, order elevation can quite efficiently reduce the locking effects and is easily achieved with the isogeometric approach. Further treatment of locking is not considered herein and is beyond of the scope of this work.

The outline of the paper is as follows: In Section 2, a brief introduction to the tangential differential calculus (TDC) is given. In Section 3, the classical linear Reissner-Mindlin shell equations are recast in terms of the TDC. Stress resultants such as membrane forces, bending moments and transverse shear forces are defined. Based on the stress resultants, the force and moment equilibriums are presented. In Section 4, the discrete weak form and the resulting system of linear equations is shown. Furthermore, different discetization strategies of the difference vector are outlined. In Section 5, numerical results are presented. The first example is the well-known Scordelis-Lo roof proposed in [4], the second example is the partly clamped hyperbolic paraboloid taken from [3], and the last example is based on the clamped flower shaped shell from [56]. 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

A shell is a thin-walled structure, which can be modelled as a surface embedded in the physical space . Let the middle surface of the shell be possibly curved, sufficiently smooth, orientable, connected and bounded by . The surface may also be called a manifold. Surfaces may be defined (explicitly) through a bijective mapping from the parameter space to the real domain , see Fig. (a). Hence, the surface is given by a parametrization or parametrized. A surface may also be defined implicitly by a level-set function following the level-set method. Then, is a scalar-valued function and the middle surface is implied by its zero-isosurface , which might be bounded by additional level-set functions as described in [34], see Fig. (b). For the formal proof of equivalence of both cases we refer to, e.g., [27].

((a)) explicit
((b)) implicit
Fig. 1: Examples of bounded surfaces embedded in the physical space : (a) explicitly defined surface with a parametrization , (b) implicitly defined surface with a master level-set function (yellow) and slave level-set functions for the boundary definition (grey).

In both cases, there exists a unit normal vector . The definition of the normal vector depends on whether the surface is based on a parametrization or not. In the parametrized case, the normal vector is determined by the normalized cross-product of the tangent vectors living in the shell surface and given by the columns of the Jacobi-matrix . In the implicit case, the normal vector is defined by the normalized gradient of the level-set function . Along the boundary , there is an associated tangential vector pointing along and a co-normal vector pointing “outwards” and being perpendicular to the boundary yet in the tangent plane of the surface .

2.1 Surface gradients and divergence

On the manifold , the orthogonal projection operator is defined based on the normal vector as


and projects arbitrary vectors into the tangent space . There holds , , and . Instead, the projection operator


projects an arbitrary vector in normal direction of and there holds , , and .

The tangential gradient operator of a differentiable scalar function on the manifold is given by


where is the standard gradient operator, and is a smooth extension of in a neighborhood of the manifold . For parametrized surfaces defined by the map , and a given scalar function , the tangential gradient may be determined without explicitly computing an extension using



being the metric tensor (first fundamental form). This is relevant in the context of the Surface FEM, where

may be the shape functions in the reference element and tangential gradients are to be determined in the physical surface elements. It is noteworthy that is in the tangent space of , i.e., , and, thus, , and . It is straightforward to determine second order derivatives of scalar functions, see, e.g., [22, 56].

When surface gradients of vector functions are considered, it is important to distinguish directional and covariant gradients:

Concerning the surface divergence of vector functions and tensor functions , there holds

To derive the weak form of the governing equations, the following divergence theorem on manifolds is needed [21, 22, 56],


where . For tangential (in-plane) tensor functions with , the term involving the mean curvature vanishes and one finds .

3 The shell equations

In this section, we derive the linear Reissner-Mindlin shell theory or first-order shear deformation theory in the frame of tangential operators based on a global Cartesian coordinate system. As mentioned before, this has the advantage over classical shell theory that the resulting model is valid no matter whether a parametrization is available or not.

We restrict ourselves to infinitesimal deformations and rotations, which means that the reference and spatial configuration are indistinguishable. For simplicity, a linear elastic material governed by Hooke’s law is assumed. In contrast to the Kirchhoff-Love shells, the additional constraint on the shell director is omitted, thus allowing transverse shear strains, leading to the well-known 5-parameter shell models, e.g., [7]. Furthermore, we assume a constant shifter in the material law, which enables an analytical pre-integration in thickness direction.

The shell continuum of thickness can be defined implicitly with a signed distance function ,


Alternatively, when the middle surface is parametrized with a map , the domain of the shell is defined by


where is the coordinate in thickness direction .

3.1 Kinematics

For Reissner-Mindlin shells, the cross section remains straight after the deformation, but not necessarily normal to the middle surface due to transverse shear deformations. Herein, the rotation of the normal vector is modelled with a standard difference vector formulation [28, 44]. Other approaches such as exponential maps, rotation tensors, etc. have been proposed, e.g., in [2, 58, 59, 7, 25], but are not considered here. The overall displacement of a point is the difference between the spatial and reference configuration

which takes the form


with being the displacement of the middle surface and being the difference vector, describing the rotation of the normal vector. In contrast to the Kirchhoff-Love shell, transverse shear deformations occur, which results in an additional rotation of the normal vector , as illustrated in Fig. 2.

Fig. 2: Displacement field of the Reissner-Mindlin shell.

The difference vector expressed in terms of the TDC is then defined as in [21, 56] with additional transverse shear deformations


Note that the difference vector is a tangential vector as in the classical theory. The surface gradient of is given by


with being the Weingarten map [42, 22]. Note that the gradient of the thickness parameter can be identified as the normal vector . The second term in Eq. 3.5 is the inverse of the shell shifter [7], which is a second order tensor. The equivalent expression in local coordinates is , which maps the metric form the shell body to the middle surface, where are the contra-variant base vectors in the shell continuum and the co-variant base vectors on the surface.

The linearised strain tensor is defined by the symmetric part of the surface gradient of


and is split into an in-plane strain and a transverse shear strain . Neglecting higher-order terms in thickness direction, as usual in the classical theory [7], the in-plane strain is defined by [56]


which is divided into a membrane and bending strain. The in-plane membrane strain becomes


and the bending strain is


The transverse shear strain is defined in a similar manner as in [40]


When the shell surface is parametrized, the resulting strain components are equivalent compared to the classical theory in local coordinates, see, e.g., [2, 7, 28, 44]. In the case of flat shell structures, the membrane strain is only a function of the tangential part of the middle surface displacement . Since the curvature is zero in this case, the Weingarten map vanishes. Therefore, the bending strain is only a function of the difference vector and the transverse shear strain becomes a function of the normal displacement and , resulting into the well-known Reissner-Mindlin plate, see, e.g., [53].

3.2 Constitutive Equation

The shell is assumed to be linear elastic and, as usual for thin structures, the Lamé constants are chosen such that the normal stress in thickness direction is eliminated, hence,


where and . The stress tensor is decomposed in a similar manner as above into in-plane (membrane and bending) stresses


and transverse shear stresses


As readily seen, the transverse shear stress is only a function of , which results in a constant transverse shear stress in thickness direction within the Reissner-Mindlin shell theory. In order account for this, a shear correction factor is introduced [7]. A common choice of the shear correction factor is . Note that due to the double projection with in Eq. 3.15 of the in-plane stress, also directional gradients can be used, which is beneficial from an implementational point of view.

3.3 Stress resultants

Assuming a constant shifter in the material law, the stress tensor is only a function of the deflection of the middle surface , the difference vector and linear in thickness direction. This enables an analytical pre-integration with respect to the thickness and stress resultants, such as effective membrane forces, bending moments and transverse shear forces can be identified. The following quantities are expressed in terms of the TDC using a global Cartesian coordinate system and are equivalent to the stress resultants in the classical theory using curvilinear coordinates, e.g., [2, 7].

The symmetric, in-plane moment tensor is defined as

resulting in the components

with and

being the flexural rigidity of the shell. The two non-zero eigenvalues of

are the principal moments and . The effective membrane (normal) force tensor is defined as

with components

where is the membrane stiffness. Analogously to the moment tensor, the effective normal force tensor is also a symmetric, in-plane tensor. Note that for curved shells this tensor is not the physical normal force tensor, but occurs in the weak form, see Section 3.5. In case of a curved shell, the physical normal force tensor is defined by


and is, in general, not symmetric, but features one zero eigenvalue just as . With Eq. 3.16, the resulting transverse shear force tensor is

with components

where is the transverse shear stiffness.

3.4 Equilibrium in strong form

Based on the stress resultants from above, one obtains the force and moment equilibrium for a curved Reissner-Mindlin shell in terms of the TDC using a global Cartesian coordinate system in strong form as


with being the load vector per area and being a distributed moment vector on the middle surface . One may split the force equilibrium into the tangential and normal direction


Alternatively, Eq. 3.21 can be rewritten in terms of the effective normal force tensor by substituting with Eq. 3.19


Assuming a bounded shell with boundary , there exist for each field (deflection of the middle surface and difference vector ) two non-overlapping parts, the Dirichlet boundary and the Neumann boundary , with . The corresponding boundary conditions for the displacement are


For the rotation of the normal vector, the boundary conditions are


In Fig. 3, the possible boundary conditions are illustrated. In Fig. (a), the displacement field at the boundary is expressed in terms of the local triad and, since the difference vector is tangential, the rotation of the normal vector may be written in terms of


The conjugated forces and bending moments at the boundary are shown in Fig. (b).

((a)) displacements and rotations
((b)) forces and moments
Fig. 3: Decomposition of the middle surface displacement , difference vector , forces and bending moments along the boundary in terms of and : (a) displacements and rotations at the boundary, (b) forces and bending moments at the boundary.

A set of common support types is given in Tab. 1. Other boundary conditions (e.g., membrane support, etc.) can be found, e.g., in [2].

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

With the boundary conditions for the displacements and rotations, the complete second-order boundary value problem (BVP) is defined. The obtained BVP in terms of the TDC is valid in the case of implicitly and explicitly defined surfaces. In the case of parametrized shells, the equilibrium in strong form is equivalent to the strong form formulated in local coordinates [2, 62, 44]. However, because the obtained BVP does not rely on a parametrized middle surface of the shell, the formulation in the frame of TDC is more general.

3.5 Equilibrium in weak form

We first formulate the weak form of the force equilibrium. Eq. 3.21 is multiplied with a suitable test function and the divergence theorem Eq. 2.5 is applied to the terms and resulting into


Note that, in order to obtain Eq. 3.29 the identities and are used. As previously mentioned, in the weak form of the force equilibrium only the effective normal force tensor appears instead of the non-symmetric, physical normal force tensor.

The weak form of the moment equilibrium is obtained in a similar manner by multiplying Eq. 3.22 with a suitable, tangential test function and the divergence theorem of Eq. 2.5 is applied to the first term, leading to


Suitable trial and the test function spaces are subspaces of the Sobolev space, where is the space of functions with square integrable first derivatives.

4 Implementational aspects

The previously derived continuous weak forms can be discretized with different finite element approaches such as the classical Surface FEM or more recent approaches such as the CutFEM [10, 11, 13, 29] and TraceFEM [35, 50, 51, 55]. Herein, the weak form of the BVP is discretized using isogeometric analysis (IGA) as proposed by Hughes et al. [41, 18], being closely related to Surface FEM when applied to solve PDEs on manifolds. It is pointed out that continuity requirements would also allow a standard FEM implementation using -continuous shape functions. Nevertheless, we prefer the use of NURBS here, for example, due to the improved convergence properties and higher smoothness of the results (including smooth forces and moments). Furthermore, the improved smoothness enables us to compute errors in the strong form of the BVP, see the numerical results in Section 5.

The NURBS patch is the middle surface of the shell and the elements are defined by the knot spans of the patch. Linking isogeometric analysis to standard FE terminology, one may naturally refer to (i) the NURBS patch as the “mesh”, (ii) the knot spans as the “elements”, and (iii) the NURBS-functions as the shape, test, and/or trial functions. The shape functions employed for all (physical) fields involved are NURBS of order with being the number of control points. We avoid the mathematical definition of NURBS and the resulting patches here because of the abundance of literature devoted to IGA and consider this as common state of the art.

The surface derivatives of the shape functions are computed as defined in Eq. 2.4 , similar to the Surface FEM [27, 23, 32, 34] using NURBS instead of Lagrange polynomials as trial and test functions. A general finite element space of order is now defined by

where the degrees of freedom

are stored at the control points.

The discrete displacement of the middle surface results as , with being Cartesian base vectors, with and . In contrast to , the difference vector is a tangential vector, describing the rotation of the normal vector. The discretization of a tangential vector is, in general, not straightforward and, in the following, different strategies are examined:

(1) In the case of a parametrized surface, the co-variant base vectors , which are by construction tangential, may be used to define the difference vector , where . This approach is used in the classical 5-parameter models [7].

(2) Alternatively, the directions of the principal curvatures, which are eigenvectors of the Weingarten map

, can be used as basis vectors. These vectors are perpendicular and also tangential by construction. This might be a reasonable choice in the case of curved, implicitly defined surfaces, where a parametrization is not available. Compared to the first approach, the crucial requirement of a parametrization is circumvented and the number of degrees of freedom per control point is equal.

(3) Another possibility is to define the difference vector in the global Cartesian coordinate system , with and and the constraint is weakly enforced with a Lagrange multiplier or with the penalty method.

(4) A variant of (3), is to project the difference vector onto the tangent space of the middle surface . An advantage of this approach is, that the additional Lagrange multiplier field is not needed. On the other hand, due to the projection, conditioning issues occur, which may be addressed with an additional stabilization term.

Herein the third approach, where the difference vector is globally defined and the constraint is enforced with a Lagrange multiplier, is chosen. The shape functions of the discrete Lagrange multiplier for the constraint on the difference vector is defined in the same manner as the components of the middle surface displacement. Furthermore, the boundary conditions shall be enforced weakly with Lagrange multipliers [64]. The shape functions of the discrete Lagrange multiplier field for the displacements is defined as and for the difference vector as .

Based on this, the following discrete trial and test functions spaces are defined:


4.1 Discrete weak form

The discrete weak form of the Reissner-Mindlin shell with Lagrange multipliers is the following. Given Young’s modulus , Poisson’s ratio , surface load and moment on , traction on , bending moments on and boundary conditions in , on , find the displacement field , the difference vector , and the Lagrange multiplier fields such that for all test functions