Reissner-Mindlin shell theory based on tangential differential calculus

TL;DR

This paper reformulates the Reissner-Mindlin shell theory using tangential differential calculus in a Cartesian framework, enabling analysis on both parametrized and level-set defined surfaces with improved generality.

Contribution

It introduces a TDC-based shell formulation that does not depend on surface parametrization, facilitating applications in modern finite element methods like TraceFEM and CutFEM.

Findings

Achieves optimal higher-order convergence rates in numerical tests.

Applicable to both classical parametrized and level-set surface representations.

Demonstrates the method's effectiveness with isogeometric analysis using NURBS.

Abstract

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.

Tables

Table 1. Tab. 1: Set of common boundary conditions

Clamped edge	$u_{{𝒕}_{\partial \Gamma }}=0$	$u_{{𝒏}_{\partial \Gamma }}=0$	$u_{{𝒏}_{\Gamma }}=0$	${\omega }_{{𝒕}_{\partial \Gamma }}=0$	${\omega }_{{𝒏}_{\partial \Gamma }}=0$
Simply supported edge	$u_{{𝒕}_{\partial \Gamma }}=0$	$u_{{𝒏}_{\partial \Gamma }}=0$	$u_{{𝒏}_{\Gamma }}=0$	$m_{{𝒕}_{\partial \Gamma }}=0$	$m_{{𝒏}_{\partial \Gamma }}=0$
Symmetry support	$p_{{𝒕}_{\partial \Gamma }}=0$	$u_{{𝒏}_{\partial \Gamma }}=0$	$p_{{𝒏}_{\Gamma }}=0$	${\omega }_{{𝒕}_{\partial \Gamma }}=0$	$m_{{𝒏}_{\partial \Gamma }}=0$
Free edge	$p_{{𝒕}_{\partial \Gamma }}=0$	$p_{{𝒏}_{\partial \Gamma }}=0$	$p_{{𝒏}_{\Gamma }}=0$	$m_{{𝒕}_{\partial \Gamma }}=0$	$m_{{𝒏}_{\partial \Gamma }}=0$

Equations

$$\mathbf{P}\left(\mathbf{x}\right)=\mathbb{I}-{\mathbf{n}}_{\Gamma }\left(\mathbf{x}\right)\otimes {\mathbf{n}}_{\Gamma }\left(\mathbf{x}\right)$$

$$\mathbf{Q}\left(\mathbf{x}\right)=\mathbb{I}-\mathbf{P}(\mathbf{x})$$

$${\nabla }_{\Gamma }u(\mathbf{x})$$

$${\nabla }_{\Gamma }u(\mathbf{x}(\mathbf{r}))=\mathbf{J}\left(\mathbf{r}\right)\cdot {\mathbf{G}}^{-1}\left(\mathbf{r}\right)\cdot {\nabla }_{\mathbf{r}}u(\mathbf{r}),$$

$${\nabla }_{\Gamma }^{\text{dir}}\mathbf{u}\left(\mathbf{x}\right)$$

$${\nabla }_{\Gamma }^{\mathrm{cov}}\mathbf{u}\left(\mathbf{x}\right)$$

$${\text{div}}_{\Gamma }\mathbf{u}\left(\mathbf{x}\right)$$

$${\text{div}}_{\Gamma }\mathbf{A}\left(\mathbf{x}\right)$$

$${\int }_{\Gamma }\mathbf{u}\cdot {\text{div}}_{\Gamma }\mathbf{A}\mathrm{d}A=-{\int }_{\Gamma }{\nabla }_{\Gamma }^{\text{dir}}\mathbf{u}:\mathbf{A}\mathrm{d}A+{\int }_{\Gamma }\varkappa \cdot \mathbf{u}\cdot \mathbf{A}\cdot {\mathbf{n}}_{\Gamma }\mathrm{d}A+{\int }_{\partial \Gamma }\mathbf{u}\cdot \mathbf{A}\cdot {\mathbf{n}}_{\partial \Gamma }\mathrm{d}s,$$

$$\Omega =\left\{\mathbf{x}\in {\mathbb{R}}^3:|{\varphi }_{\text{SDF}}(\mathbf{x})|\le \frac{t}{2}\right\}.$$

$$\mathbf{x}={\mathbf{x}}_{\Gamma }+\zeta {\mathbf{n}}_{\Gamma }({\mathbf{x}}_{\Gamma }),$$

$${\mathbf{u}}_{\Omega }(\mathbf{x})$$

$${\mathbf{u}}_{\Omega }({\mathbf{x}}_{\Gamma },\zeta )=\mathbf{u}({\mathbf{x}}_{\Gamma })+\zeta \mathbf{w}({\mathbf{x}}_{\Gamma })$$

$$\mathbf{w}({\mathbf{x}}_{\Gamma })=-[{\nabla }_{\Gamma }^{\text{dir}}\mathbf{u}({\mathbf{x}}_{\Gamma }){]}^{\mathrm{T}}\cdot {\mathbf{n}}_{\Gamma }+\mathbf{\gamma }({\mathbf{x}}_{\Gamma }).$$

$$\begin{array}{rl}{\nabla }_{\Gamma }^{\text{dir}}{\mathbf{u}}_{\Omega }(\mathbf{x})& =\mathbf{P}\cdot \frac{\partial {\mathbf{u}}_{\Omega }(\mathbf{x})}{\partial \mathbf{x}}\cdot \frac{\partial \mathbf{x}}{\partial {\mathbf{x}}_{\Gamma }},\\ & =\mathbf{P}\cdot \left(\nabla \tilde{\mathbf{u}}+\nabla \zeta \otimes \mathbf{w}+\zeta \nabla \tilde{\mathbf{w}}\right)\cdot \left(\mathbb{I}+\zeta \mathbf{H}\right),\\ & =\left({\nabla }_{\Gamma }^{\text{dir}}\mathbf{u}+{\mathbf{n}}_{\Gamma }\otimes \mathbf{w}+\zeta {\nabla }_{\Gamma }^{\text{dir}}\mathbf{w}\right)\cdot \left(\mathbb{I}+\zeta \mathbf{H}\right),\end{array}$$

$${\mathbf{\epsilon }}_{\Gamma }(\mathbf{x})$$

$$={\mathbf{\epsilon }}_{\Gamma }^{\text{P}}(\mathbf{x})+{\mathbf{\epsilon }}_{\Gamma }^{\text{S}}(\mathbf{x})$$

$${\mathbf{\epsilon }}_{\Gamma }^{\text{P}}$$

$$={\mathbf{\epsilon }}_{\Gamma ,\text{Mem}}^{\text{P}}+\zeta {\mathbf{\epsilon }}_{\Gamma ,\text{Bend}}^{\text{P}},$$

$${\mathbf{\epsilon }}_{\Gamma ,\text{Mem}}^{\text{P}}(\mathbf{u})=\frac{1}{2}\left[{\nabla }_{\Gamma }^{\text{cov}}\mathbf{u}+({\nabla }_{\Gamma }^{\text{cov}}\mathbf{u}{)}^{\mathrm{T}}\right],$$

$${\mathbf{\epsilon }}_{\Gamma ,\text{Bend}}^{\text{P}}(\mathbf{u},\mathbf{w})=\frac{1}{2}\left[\mathbf{H}\cdot {\nabla }_{\Gamma }^{\text{dir}}\mathbf{u}+({\nabla }_{\Gamma }^{\text{dir}}\mathbf{u}{)}^{\mathrm{T}}\cdot \mathbf{H}+{\nabla }_{\Gamma }^{\text{cov}}\mathbf{w}+({\nabla }_{\Gamma }^{\text{cov}}\mathbf{w}{)}^{\mathrm{T}}\right].$$

$${\mathbf{\epsilon }}_{\Gamma }^{\text{S}}(\mathbf{u},\mathbf{w})$$

$$=\frac{1}{2}\left[\mathbf{Q}\cdot {\nabla }_{\Gamma }^{\text{dir}}\mathbf{u}+({\nabla }_{\Gamma }^{\text{dir}}\mathbf{u}{)}^{\mathrm{T}}\cdot \mathbf{Q}+{\mathbf{n}}_{\Gamma }\otimes \mathbf{w}+\mathbf{w}\otimes {\mathbf{n}}_{\Gamma }\right].$$

$${\mathbf{\sigma }}_{\Gamma }(\mathbf{x})=2\mu {\mathbf{\epsilon }}_{\Gamma }(\mathbf{x})+\lambda \text{tr}[{\mathbf{\epsilon }}_{\Gamma }(\mathbf{x})]\mathbb{I}$$

$${\mathbf{\sigma }}_{\Gamma }^{\text{P}}(\mathbf{x})=\mathbf{P}\cdot \left[2\mu {\mathbf{\epsilon }}_{\Gamma }^{\text{P}}(\mathbf{x})+\lambda \text{tr}[{\mathbf{\epsilon }}_{\Gamma }^{\text{P}}(\mathbf{x})]\mathbb{I}\right]\cdot \mathbf{P},$$

$$\begin{array}{rl}{\mathbf{\sigma }}_{\Gamma }^{\text{S}}({\mathbf{x}}_{\Gamma })& =2\mu \left[\mathbf{Q}\cdot {\mathbf{\epsilon }}_{\Gamma }^{\text{S}}({\mathbf{x}}_{\Gamma })+{\mathbf{\epsilon }}_{\Gamma }^{\text{S}}({\mathbf{x}}_{\Gamma })\cdot \mathbf{Q}\right]+\lambda \text{tr}\left[{\mathbf{\epsilon }}_{\Gamma }^{\text{S}}({\mathbf{x}}_{\Gamma })\right]\mathbf{Q}\\ & =2\mu {\alpha }_{\text{s}}{\mathbf{\epsilon }}_{\Gamma }^{\text{S}}({\mathbf{x}}_{\Gamma }).\end{array}$$

$${\mathbf{m}}_{\Gamma }$$

$$=\mathbf{P}\cdot {\mathbf{m}}_{\Gamma }^{\text{dir}}\cdot \mathbf{P},$$

$${\left[{\mathbf{m}}_{\Gamma }^{\text{dir}}\right]}_{11}$$

$${\left[{\mathbf{m}}_{\Gamma }^{\text{dir}}\right]}_{22}$$

$${\left[{\mathbf{m}}_{\Gamma }^{\text{dir}}\right]}_{33}$$

$${\left[{\mathbf{m}}_{\Gamma }^{\text{dir}}\right]}_{12}$$

$${\left[{\mathbf{m}}_{\Gamma }^{\text{dir}}\right]}_{13}$$

$${\left[{\mathbf{m}}_{\Gamma }^{\text{dir}}\right]}_{23}$$

$${\tilde{\mathbf{n}}}_{\Gamma }$$

$$=\mathbf{P}\cdot {\mathbf{n}}_{\Gamma }^{\text{dir}}\cdot \mathbf{P},$$

$${\left[{\mathbf{n}}_{\Gamma }^{\text{dir}}\right]}_{11}$$

Full text

Keywords: Shells, Tangential Differential Calculus, Isogeometric analysis, Manifolds

Institute of Structural Analysis

Graz University of Technology

Lessingstr. 25/II, 8010 Graz, Austria

www.ifb.tugraz.at

schoellhammer@tugraz.at

1. 1 Introduction

2. 2 Preliminaries

3. 2.1 Surface gradients and divergence

4. 3 The shell equations

5. 3.1 Kinematics

6. 3.2 Constitutive Equation

7. 3.3 Stress resultants

8. 3.4 Equilibrium in strong form

9. 3.5 Equilibrium in weak form

10. 4 Implementational aspects

11. 4.1 Discrete weak form

12. 5 Numerical results

13. 5.1 Scordelis-Lo roof

14. 5.2 Partly clamped hyperbolic paraboloid

15. 5.3 Clamped flower shaped shell

16. 6 Conclusions and Outlook

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 $1$ and only $C^{0}$-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 $\mathbb{R}^{3}$. 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 $\mathbb{R}^{3}$. 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. The same results are obtained here, however, based on the new TDC-based formulation resulting into a significant different implementation. In particular, the implementation is more intuitive and compact due to a sharper split of FE-technology and application, as already shown in [56]. 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 $C^{0}$-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 $\Gamma$ embedded in the physical space $\mathbb{R}^{3}$. Let the middle surface of the shell be possibly curved, sufficiently smooth, orientable, connected and bounded by $\partial\Gamma$. The surface may also be called a manifold. Surfaces may be defined (explicitly) through a bijective mapping $\mathchoice{\displaystyle\bm{x}}{\textstyle\bm{x}}{\scriptstyle\bm{x}}{\scriptscriptstyle\bm{x}}(\mathchoice{\displaystyle\bm{r}}{\textstyle\bm{r}}{\scriptstyle\bm{r}}{\scriptscriptstyle\bm{r}}):\hat{\Omega}\to\Gamma$ from the parameter space $\hat{\Omega}\subset\mathbb{R}^{2}$ to the real domain $\Gamma\subset\mathbb{R}^{3}$, see Fig. 1(a). In the case of a surface mesh, this is rather an atlas of local, element-wise mappings. In these cases, the surface is given by a parametrization or parametrized. A surface may also be defined implicitly by a level-set function $\phi(\mathchoice{\displaystyle\bm{x}}{\textstyle\bm{x}}{\scriptstyle\bm{x}}{\scriptscriptstyle\bm{x}}):\mathbb{R}^{3}\to\mathbb{R}$ following the level-set method. Then, $\phi$ is a scalar-valued function and the middle surface is implied by its zero-isosurface $\phi(\mathchoice{\displaystyle\bm{x}}{\textstyle\bm{x}}{\scriptstyle\bm{x}}{\scriptscriptstyle\bm{x}})=0$, which might be bounded by additional level-set functions as described in [34], see Fig. 1(b). For the formal proof of equivalence of both cases we refer to, e.g., [27].

In both cases, there exists a unit normal vector $\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}=[n_{x},\,n_{y},\,n_{z}]^{\mathrm{T}}\in\mathbb{R}^{3}$. The representation 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 $\mathchoice{\displaystyle\mathbf{J}}{\textstyle\mathbf{J}}{\scriptstyle\mathbf{J}}{\scriptscriptstyle\mathbf{J}}(\mathchoice{\displaystyle\bm{r}}{\textstyle\bm{r}}{\scriptstyle\bm{r}}{\scriptscriptstyle\bm{r}})=\nicefrac{{\partial\mathchoice{\displaystyle\bm{x}}{\textstyle\bm{x}}{\scriptstyle\bm{x}}{\scriptscriptstyle\bm{x}}}}{{\partial\mathchoice{\displaystyle\bm{r}}{\textstyle\bm{r}}{\scriptstyle\bm{r}}{\scriptscriptstyle\bm{r}}}}$. In the implicit case, the normal vector is given by the normalized gradient of the level-set function $\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\Gamma}=\nicefrac{{\nabla\phi}}{{\|\nabla\phi\|}}$. Along the boundary $\partial\Gamma$, there is an associated tangential vector $\mathchoice{\displaystyle\bm{t}}{\textstyle\bm{t}}{\scriptstyle\bm{t}}{\scriptscriptstyle\bm{t}}_{\partial\Gamma}\in\mathbb{R}^{3}$ pointing along $\partial\Gamma$ and a co-normal vector $\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\partial\Gamma}=\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\Gamma}\times\mathchoice{\displaystyle\bm{t}}{\textstyle\bm{t}}{\scriptstyle\bm{t}}{\scriptscriptstyle\bm{t}}_{\partial\Gamma}\in\mathbb{R}^{3}$ pointing “outwards” and being perpendicular to the boundary $\partial\Gamma$ yet in the tangent plane of the surface $\Gamma$.

2.1 Surface gradients and divergence

On the manifold $\Gamma$, the orthogonal projection operator $\mathchoice{\displaystyle\mathbf{P}}{\textstyle\mathbf{P}}{\scriptstyle\mathbf{P}}{\scriptscriptstyle\mathbf{P}}\left(\mathchoice{\displaystyle\bm{x}}{\textstyle\bm{x}}{\scriptstyle\bm{x}}{\scriptscriptstyle\bm{x}}\right)\in\mathbb{R}^{3\times 3}$ is defined based on the normal vector as

[TABLE]

and projects arbitrary vectors $\mathchoice{\displaystyle\bm{v}}{\textstyle\bm{v}}{\scriptstyle\bm{v}}{\scriptscriptstyle\bm{v}}(\Gamma)\in\mathbb{R}^{3}$ into the tangent space $T_{P}\Gamma$. There holds $\mathchoice{\displaystyle\mathbf{P}}{\textstyle\mathbf{P}}{\scriptstyle\mathbf{P}}{\scriptscriptstyle\mathbf{P}}\cdot\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\Gamma}=\mathchoice{\displaystyle\bm{0}}{\textstyle\bm{0}}{\scriptstyle\bm{0}}{\scriptscriptstyle\bm{0}}$, $\mathchoice{\displaystyle\mathbf{P}}{\textstyle\mathbf{P}}{\scriptstyle\mathbf{P}}{\scriptscriptstyle\mathbf{P}}=\mathchoice{\displaystyle\mathbf{P}}{\textstyle\mathbf{P}}{\scriptstyle\mathbf{P}}{\scriptscriptstyle\mathbf{P}}^{{\mathrm{T}}}$, and $\mathchoice{\displaystyle\mathbf{P}}{\textstyle\mathbf{P}}{\scriptstyle\mathbf{P}}{\scriptscriptstyle\mathbf{P}}\cdot\mathchoice{\displaystyle\mathbf{P}}{\textstyle\mathbf{P}}{\scriptstyle\mathbf{P}}{\scriptscriptstyle\mathbf{P}}=\mathchoice{\displaystyle\mathbf{P}}{\textstyle\mathbf{P}}{\scriptstyle\mathbf{P}}{\scriptscriptstyle\mathbf{P}}$. Instead, the projection operator

[TABLE]

projects an arbitrary vector in normal direction of $\Gamma$ and there holds $\mathchoice{\displaystyle\mathbf{Q}}{\textstyle\mathbf{Q}}{\scriptstyle\mathbf{Q}}{\scriptscriptstyle\mathbf{Q}}\cdot\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\Gamma}=\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\Gamma}$, $\mathchoice{\displaystyle\mathbf{Q}}{\textstyle\mathbf{Q}}{\scriptstyle\mathbf{Q}}{\scriptscriptstyle\mathbf{Q}}=\mathchoice{\displaystyle\mathbf{Q}}{\textstyle\mathbf{Q}}{\scriptstyle\mathbf{Q}}{\scriptscriptstyle\mathbf{Q}}^{\mathrm{T}}$, and $\mathchoice{\displaystyle\mathbf{Q}}{\textstyle\mathbf{Q}}{\scriptstyle\mathbf{Q}}{\scriptscriptstyle\mathbf{Q}}\cdot\mathchoice{\displaystyle\mathbf{Q}}{\textstyle\mathbf{Q}}{\scriptstyle\mathbf{Q}}{\scriptscriptstyle\mathbf{Q}}=\mathchoice{\displaystyle\mathbf{Q}}{\textstyle\mathbf{Q}}{\scriptstyle\mathbf{Q}}{\scriptscriptstyle\mathbf{Q}}$.

The tangential gradient operator $\nabla_{\Gamma}$ of a differentiable scalar function $u:\Gamma\to\mathbb{R}$ on the manifold is given by

[TABLE]

where $\nabla$ is the standard gradient operator, and $\tilde{u}$ is a smooth extension of $u$ in a neighborhood $\mathcal{U}$ of the manifold $\Gamma$. For parametrized surfaces defined by the map $\mathchoice{\displaystyle\bm{x}}{\textstyle\bm{x}}{\scriptstyle\bm{x}}{\scriptscriptstyle\bm{x}}\left(\mathchoice{\displaystyle\bm{r}}{\textstyle\bm{r}}{\scriptstyle\bm{r}}{\scriptscriptstyle\bm{r}}\right)$, and a given scalar function $u\left(\mathchoice{\displaystyle\bm{r}}{\textstyle\bm{r}}{\scriptstyle\bm{r}}{\scriptscriptstyle\bm{r}}\right):\mathbb{R}^{2}\to\mathbb{R}$, the tangential gradient may be determined without explicitly computing an extension $\tilde{u}$ using

[TABLE]

with $\mathchoice{\displaystyle\mathbf{G}}{\textstyle\mathbf{G}}{\scriptstyle\mathbf{G}}{\scriptscriptstyle\mathbf{G}}=\mathchoice{\displaystyle\mathbf{J}}{\textstyle\mathbf{J}}{\scriptstyle\mathbf{J}}{\scriptscriptstyle\mathbf{J}}^{\mathrm{T}}\cdot\mathchoice{\displaystyle\mathbf{J}}{\textstyle\mathbf{J}}{\scriptstyle\mathbf{J}}{\scriptscriptstyle\mathbf{J}}$ being the metric tensor (first fundamental form). This is relevant in the context of the Surface FEM, where $u\left(\mathchoice{\displaystyle\bm{r}}{\textstyle\bm{r}}{\scriptstyle\bm{r}}{\scriptscriptstyle\bm{r}}\right)$ 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 $\nabla_{\Gamma}u$ is in the tangent space of $\Gamma$, i.e., $\nabla_{\Gamma}u\in T_{P}\Gamma$, and, thus, $\mathchoice{\displaystyle\mathbf{P}}{\textstyle\mathbf{P}}{\scriptstyle\mathbf{P}}{\scriptscriptstyle\mathbf{P}}\cdot\nabla_{\Gamma}u=\nabla_{\Gamma}u$, $\mathchoice{\displaystyle\mathbf{Q}}{\textstyle\mathbf{Q}}{\scriptstyle\mathbf{Q}}{\scriptscriptstyle\mathbf{Q}}\cdot\nabla_{\Gamma}u=\mathchoice{\displaystyle\bm{0}}{\textstyle\bm{0}}{\scriptstyle\bm{0}}{\scriptscriptstyle\bm{0}}$ and $\nabla_{\Gamma}u\cdot\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\Gamma}=0$. It is straightforward to determine second order derivatives of scalar functions, see, e.g., [22, 56].

When surface gradients of vector functions $\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}\left(\mathchoice{\displaystyle\bm{x}}{\textstyle\bm{x}}{\scriptstyle\bm{x}}{\scriptscriptstyle\bm{x}}\right):\Gamma\to\mathbb{R}^{3}$ are considered, it is important to distinguish directional and covariant gradients:

[TABLE]

Concerning the surface divergence of vector functions $\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}\left(\mathchoice{\displaystyle\bm{x}}{\textstyle\bm{x}}{\scriptstyle\bm{x}}{\scriptscriptstyle\bm{x}}\right):\Gamma\to\mathbb{R}^{3}$ and tensor functions $\mathchoice{\displaystyle\mathbf{A}}{\textstyle\mathbf{A}}{\scriptstyle\mathbf{A}}{\scriptscriptstyle\mathbf{A}}\left(\mathchoice{\displaystyle\bm{x}}{\textstyle\bm{x}}{\scriptstyle\bm{x}}{\scriptscriptstyle\bm{x}}\right):\Gamma\to\mathbb{R}^{3\times 3}$, there holds

[TABLE]

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

[TABLE]

where $\nabla_{\Gamma}^{\text{dir}}\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}:\mathchoice{\displaystyle\mathbf{A}}{\textstyle\mathbf{A}}{\scriptstyle\mathbf{A}}{\scriptscriptstyle\mathbf{A}}=\mathrm{tr}\left(\nabla_{\Gamma}^{\text{dir}}\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}\cdot\mathchoice{\displaystyle\mathbf{A}}{\textstyle\mathbf{A}}{\scriptstyle\mathbf{A}}{\scriptscriptstyle\mathbf{A}}^{{\mathrm{T}}}\right)$. For tangential (in-plane) tensor functions with $\mathchoice{\displaystyle\mathbf{A}}{\textstyle\mathbf{A}}{\scriptstyle\mathbf{A}}{\scriptscriptstyle\mathbf{A}}=\mathchoice{\displaystyle\mathbf{P}}{\textstyle\mathbf{P}}{\scriptstyle\mathbf{P}}{\scriptscriptstyle\mathbf{P}}\cdot\mathchoice{\displaystyle\mathbf{A}}{\textstyle\mathbf{A}}{\scriptstyle\mathbf{A}}{\scriptscriptstyle\mathbf{A}}\cdot\mathchoice{\displaystyle\mathbf{P}}{\textstyle\mathbf{P}}{\scriptstyle\mathbf{P}}{\scriptscriptstyle\mathbf{P}}$, the term involving the mean curvature $\varkappa=\text{tr}(\nabla_{\Gamma}^{\text{dir}}\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\Gamma})$ vanishes and one finds $\nabla_{\Gamma}^{\text{dir}}\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}:\mathchoice{\displaystyle\mathbf{A}}{\textstyle\mathbf{A}}{\scriptstyle\mathbf{A}}{\scriptscriptstyle\mathbf{A}}=\nabla_{\Gamma}^{\text{cov}}\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}:\mathchoice{\displaystyle\mathbf{A}}{\textstyle\mathbf{A}}{\scriptstyle\mathbf{A}}{\scriptscriptstyle\mathbf{A}}$.

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 $\Omega$ of thickness $t$ can be defined implicitly by a special level-set function, called signed-distance function $\phi_{\text{SDF}}$. This scalar-valued function gives the shortest Euclidean distance from $\mathchoice{\displaystyle\bm{x}}{\textstyle\bm{x}}{\scriptstyle\bm{x}}{\scriptscriptstyle\bm{x}}$ to the middle surface of the shell and is positive in direction of the normal vector $\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\Gamma}$ and negative otherwise. Then,

[TABLE]

Alternatively, when the middle surface $\Gamma$ is parametrized with a map $\mathchoice{\displaystyle\bm{x}}{\textstyle\bm{x}}{\scriptstyle\bm{x}}{\scriptscriptstyle\bm{x}}_{\Gamma}(\mathchoice{\displaystyle\bm{r}}{\textstyle\bm{r}}{\scriptstyle\bm{r}}{\scriptscriptstyle\bm{r}})$, the domain of the shell is defined by

[TABLE]

where $\zeta$ is the coordinate in thickness direction $|\zeta|\leq\nicefrac{{t}}{{2}}$.

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 $\mathchoice{\displaystyle\bm{P}}{\textstyle\bm{P}}{\scriptstyle\bm{P}}{\scriptscriptstyle\bm{P}}\in\Omega$ is the difference between the spatial and reference configuration

[TABLE]

which takes the form

[TABLE]

with $\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}(\mathchoice{\displaystyle\bm{x}}{\textstyle\bm{x}}{\scriptstyle\bm{x}}{\scriptscriptstyle\bm{x}}_{\Gamma}):\Gamma\to\mathbb{R}^{3}$ being the displacement of the middle surface and $\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}(\mathchoice{\displaystyle\bm{x}}{\textstyle\bm{x}}{\scriptstyle\bm{x}}{\scriptscriptstyle\bm{x}}_{\Gamma}):\Gamma\to T_{P}\Gamma$ being the difference vector, describing the rotation of the normal vector. In contrast to the Kirchhoff-Love shell, transverse shear deformations $\mathchoice{\displaystyle\bm{\gamma}}{\textstyle\bm{\gamma}}{\scriptstyle\bm{\gamma}}{\scriptscriptstyle\bm{\gamma}}$ occur, which results in an additional rotation of the normal vector $\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\Gamma}$, as illustrated in Fig. 2.

The difference vector $\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}$ expressed in terms of the TDC is then defined as in [21, 56] with additional transverse shear deformations

[TABLE]

Note that the difference vector is a tangential vector as in the classical theory. The surface gradient of $\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}_{\Omega}(\mathchoice{\displaystyle\bm{x}}{\textstyle\bm{x}}{\scriptstyle\bm{x}}{\scriptscriptstyle\bm{x}})$ is given by

[TABLE]

with $\mathchoice{\displaystyle\mathbf{H}}{\textstyle\mathbf{H}}{\scriptstyle\mathbf{H}}{\scriptscriptstyle\mathbf{H}}=\nabla_{\Gamma}^{\text{dir}}\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\Gamma}$ being the Weingarten map [42, 22]. Note that the gradient of the thickness parameter $\zeta$ can be identified as the normal vector $\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\Gamma}$. The second term $(\mathbb{I}+\zeta\mathchoice{\displaystyle\mathbf{H}}{\textstyle\mathbf{H}}{\scriptstyle\mathbf{H}}{\scriptscriptstyle\mathbf{H}})$ 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 $\sum_{i=1}^{3}\mathchoice{\displaystyle\bm{G}}{\textstyle\bm{G}}{\scriptstyle\bm{G}}{\scriptscriptstyle\bm{G}}^{i}\otimes\mathchoice{\displaystyle\bm{A}}{\textstyle\bm{A}}{\scriptstyle\bm{A}}{\scriptscriptstyle\bm{A}}_{i}$, where $\mathchoice{\displaystyle\bm{G}}{\textstyle\bm{G}}{\scriptstyle\bm{G}}{\scriptscriptstyle\bm{G}}^{i}$ are the contra-variant base vectors in the shell continuum and $\mathchoice{\displaystyle\bm{A}}{\textstyle\bm{A}}{\scriptstyle\bm{A}}{\scriptscriptstyle\bm{A}}_{i}$ the co-variant base vectors on the surface [7].

The linearised strain tensor $\bm{\varepsilon}_{\Gamma}$ is defined by the symmetric part of the surface gradient of $\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}_{\Omega}$

[TABLE]

and is split into an in-plane strain $\bm{\varepsilon}_{\Gamma}^{\text{P}}$ and a transverse shear strain $\bm{\varepsilon}_{\Gamma}^{\text{S}}$. Neglecting higher-order terms in thickness direction, as usual in the classical theory [7], the in-plane strain is defined by [56]

[TABLE]

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

[TABLE]

and the bending strain is

[TABLE]

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

[TABLE]

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 $\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}_{t}=\mathchoice{\displaystyle\mathbf{P}}{\textstyle\mathbf{P}}{\scriptstyle\mathbf{P}}{\scriptscriptstyle\mathbf{P}}\cdot\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}$. Since the curvature is zero in this case, the Weingarten map $\mathchoice{\displaystyle\mathbf{H}}{\textstyle\mathbf{H}}{\scriptstyle\mathbf{H}}{\scriptscriptstyle\mathbf{H}}$ vanishes. Therefore, the bending strain is only a function of the difference vector $\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}$ and the transverse shear strain becomes a function of the normal displacement $u_{n}=\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}\cdot\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\Gamma}$ and $\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}$, 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,

[TABLE]

where $\mu=\frac{E}{2(1+\nu)}$ and $\lambda=\frac{E\nu}{1-\nu^{2}}$. The stress tensor is decomposed in a similar manner as above into in-plane (membrane and bending) stresses

[TABLE]

and transverse shear stresses

[TABLE]

As readily seen, the transverse shear stress is only a function of $\mathchoice{\displaystyle\bm{x}}{\textstyle\bm{x}}{\scriptstyle\bm{x}}{\scriptscriptstyle\bm{x}}_{\Gamma}$, 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 $\alpha_{\text{s}}$ is introduced [7]. A common choice of the shear correction factor is $\alpha_{\text{s}}=\nicefrac{{5}}{{6}}$. Note that due to the double projection with $\mathchoice{\displaystyle\mathbf{P}}{\textstyle\mathbf{P}}{\scriptstyle\mathbf{P}}{\scriptscriptstyle\mathbf{P}}$ 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 $\bm{\sigma}_{\Gamma}(\mathchoice{\displaystyle\bm{x}}{\textstyle\bm{x}}{\scriptstyle\bm{x}}{\scriptscriptstyle\bm{x}})$ is only a function of the deflection of the middle surface $\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}$, the difference vector $\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}$ 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 $\mathchoice{\displaystyle\mathbf{m}}{\textstyle\mathbf{m}}{\scriptstyle\mathbf{m}}{\scriptscriptstyle\mathbf{m}}_{\Gamma}$ is defined as

[TABLE]

with $j=1,\,2,\,3$ and $D_{\text{B}}=\frac{Et^{3}}{12(1-\nu^{2})}$ being the flexural rigidity of the shell. The two non-zero eigenvalues of $\mathchoice{\displaystyle\mathbf{m}}{\textstyle\mathbf{m}}{\scriptstyle\mathbf{m}}{\scriptscriptstyle\mathbf{m}}_{\Gamma}$ are the principal moments $m_{1}$ and $m_{2}$. The effective membrane (normal) force tensor $\tilde{\bm{n}}_{\Gamma}$ is defined as

[TABLE]

where $D_{\text{M}}=\frac{Et}{1-\nu^{2}}$ 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 $\mathchoice{\displaystyle\mathbf{n}}{\textstyle\mathbf{n}}{\scriptstyle\mathbf{n}}{\scriptscriptstyle\mathbf{n}}^{\text{real}}_{\Gamma}$ is defined by

[TABLE]

and is, in general, not symmetric, but features one zero eigenvalue just as $\tilde{\mathchoice{\displaystyle\mathbf{n}}{\textstyle\mathbf{n}}{\scriptstyle\mathbf{n}}{\scriptscriptstyle\mathbf{n}}}_{\Gamma}$. With Eq. 3.16, the resulting transverse shear force tensor is

[TABLE]

where $D_{\text{Shear}}=\alpha_{\text{s}}\mu t=\alpha_{\text{s}}\frac{Et}{2(1+\nu)}$ 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

[TABLE]

with $\mathchoice{\displaystyle\bm{f}}{\textstyle\bm{f}}{\scriptstyle\bm{f}}{\scriptscriptstyle\bm{f}}$ being the load vector per area and $\mathchoice{\displaystyle\bm{c}}{\textstyle\bm{c}}{\scriptstyle\bm{c}}{\scriptscriptstyle\bm{c}}$ being a distributed moment vector on the middle surface $\Gamma$. One may split the force equilibrium into the tangential and normal direction

[TABLE]

Alternatively, Eq. 3.21 can be rewritten in terms of the effective normal force tensor by substituting $\mathchoice{\displaystyle\mathbf{n}}{\textstyle\mathbf{n}}{\scriptstyle\mathbf{n}}{\scriptscriptstyle\mathbf{n}}^{\text{real}}_{\Gamma}$ with Eq. 3.19

[TABLE]

Assuming a bounded shell with boundary $\partial\Gamma$, there exist for each field (deflection of the middle surface $\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}$ and difference vector $\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}$) two non-overlapping parts, the Dirichlet boundary $\partial\Gamma_{\text{D,i}}$ and the Neumann boundary $\partial\Gamma_{\text{N,i}}$, with $i=\{\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}},\,\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}\}$. The corresponding boundary conditions for the displacement $\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}$ are

[TABLE]

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

[TABLE]

In Fig. 3, the possible boundary conditions are illustrated. In Fig. 3(a), the displacement field $\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}$ at the boundary is expressed in terms of the local triad $(\mathchoice{\displaystyle\bm{t}}{\textstyle\bm{t}}{\scriptstyle\bm{t}}{\scriptscriptstyle\bm{t}}_{\partial\Gamma},\,\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\partial\Gamma},\,\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\Gamma})$ and, since the difference vector is tangential, the rotation of the normal vector may be written in terms of $(\mathchoice{\displaystyle\bm{t}}{\textstyle\bm{t}}{\scriptstyle\bm{t}}{\scriptscriptstyle\bm{t}}_{\partial\Gamma},\,\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\partial\Gamma})$

[TABLE]

The conjugated forces $(\mathchoice{\displaystyle\bm{p}}{\textstyle\bm{p}}{\scriptstyle\bm{p}}{\scriptscriptstyle\bm{p}}_{\mathchoice{\displaystyle\bm{t}}{\textstyle\bm{t}}{\scriptstyle\bm{t}}{\scriptscriptstyle\bm{t}}_{\partial\Gamma}},\,\mathchoice{\displaystyle\bm{p}}{\textstyle\bm{p}}{\scriptstyle\bm{p}}{\scriptscriptstyle\bm{p}}_{\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\partial\Gamma}},\,\mathchoice{\displaystyle\bm{p}}{\textstyle\bm{p}}{\scriptstyle\bm{p}}{\scriptscriptstyle\bm{p}}_{\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\Gamma}})$ and bending moments $(\mathchoice{\displaystyle\bm{m}}{\textstyle\bm{m}}{\scriptstyle\bm{m}}{\scriptscriptstyle\bm{m}}_{\mathchoice{\displaystyle\bm{t}}{\textstyle\bm{t}}{\scriptstyle\bm{t}}{\scriptscriptstyle\bm{t}}_{\partial\Gamma}},\,\mathchoice{\displaystyle\bm{m}}{\textstyle\bm{m}}{\scriptstyle\bm{m}}{\scriptscriptstyle\bm{m}}_{\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\partial\Gamma}})$ at the boundary are shown in Fig. 3(b).

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

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 $\mathchoice{\displaystyle\bm{v}}{\textstyle\bm{v}}{\scriptstyle\bm{v}}{\scriptscriptstyle\bm{v}}_{\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}}$ and the divergence theorem Eq. 2.5 is applied to the terms $-\int_{\Gamma}\mathchoice{\displaystyle\bm{v}}{\textstyle\bm{v}}{\scriptstyle\bm{v}}{\scriptscriptstyle\bm{v}}_{u}\cdot\text{div}_{\Gamma}\mathchoice{\displaystyle\mathbf{n}}{\textstyle\mathbf{n}}{\scriptstyle\mathbf{n}}{\scriptscriptstyle\mathbf{n}}^{\text{real}}_{\Gamma}\ \mathrm{d}A$ and $-\int_{\Gamma}\mathchoice{\displaystyle\bm{v}}{\textstyle\bm{v}}{\scriptstyle\bm{v}}{\scriptscriptstyle\bm{v}}_{u}\cdot\mathchoice{\displaystyle\mathbf{Q}}{\textstyle\mathbf{Q}}{\scriptstyle\mathbf{Q}}{\scriptscriptstyle\mathbf{Q}}\cdot\text{div}_{\Gamma}\mathchoice{\displaystyle\mathbf{q}}{\textstyle\mathbf{q}}{\scriptstyle\mathbf{q}}{\scriptscriptstyle\mathbf{q}}_{\Gamma}\ \mathrm{d}A$ resulting into

[TABLE]

Note that, in order to obtain Eq. 3.29 the identities $(\mathchoice{\displaystyle\mathbf{H}}{\textstyle\mathbf{H}}{\scriptstyle\mathbf{H}}{\scriptscriptstyle\mathbf{H}}\cdot\nabla_{\Gamma}^{\text{dir}}\mathchoice{\displaystyle\bm{v}}{\textstyle\bm{v}}{\scriptstyle\bm{v}}{\scriptscriptstyle\bm{v}}_{u}):\mathchoice{\displaystyle\mathbf{m}}{\textstyle\mathbf{m}}{\scriptstyle\mathbf{m}}{\scriptscriptstyle\mathbf{m}}_{\Gamma}=\nabla_{\Gamma}^{\text{dir}}\mathchoice{\displaystyle\bm{v}}{\textstyle\bm{v}}{\scriptstyle\bm{v}}{\scriptscriptstyle\bm{v}}_{u}:(\mathchoice{\displaystyle\mathbf{H}}{\textstyle\mathbf{H}}{\scriptstyle\mathbf{H}}{\scriptscriptstyle\mathbf{H}}\cdot\mathchoice{\displaystyle\mathbf{m}}{\textstyle\mathbf{m}}{\scriptstyle\mathbf{m}}{\scriptscriptstyle\mathbf{m}}_{\Gamma})$ and $[\mathchoice{\displaystyle\mathbf{Q}}{\textstyle\mathbf{Q}}{\scriptstyle\mathbf{Q}}{\scriptscriptstyle\mathbf{Q}}_{,x}^{\Gamma}\cdot\mathchoice{\displaystyle\bm{v}}{\textstyle\bm{v}}{\scriptstyle\bm{v}}{\scriptscriptstyle\bm{v}}_{u}\quad\mathchoice{\displaystyle\mathbf{Q}}{\textstyle\mathbf{Q}}{\scriptstyle\mathbf{Q}}{\scriptscriptstyle\mathbf{Q}}_{,y}^{\Gamma}\cdot\mathchoice{\displaystyle\bm{v}}{\textstyle\bm{v}}{\scriptstyle\bm{v}}{\scriptscriptstyle\bm{v}}_{u}\quad\mathchoice{\displaystyle\mathbf{Q}}{\textstyle\mathbf{Q}}{\scriptstyle\mathbf{Q}}{\scriptscriptstyle\mathbf{Q}}_{,z}^{\Gamma}\cdot\mathchoice{\displaystyle\bm{v}}{\textstyle\bm{v}}{\scriptstyle\bm{v}}{\scriptscriptstyle\bm{v}}_{u}]:\mathchoice{\displaystyle\mathbf{q}}{\textstyle\mathbf{q}}{\scriptstyle\mathbf{q}}{\scriptscriptstyle\mathbf{q}}_{\Gamma}=\mathchoice{\displaystyle\bm{v}}{\textstyle\bm{v}}{\scriptstyle\bm{v}}{\scriptscriptstyle\bm{v}}_{u}^{\mathrm{T}}\cdot\mathchoice{\displaystyle\mathbf{H}}{\textstyle\mathbf{H}}{\scriptstyle\mathbf{H}}{\scriptscriptstyle\mathbf{H}}\cdot(\mathchoice{\displaystyle\mathbf{q}}{\textstyle\mathbf{q}}{\scriptstyle\mathbf{q}}{\scriptscriptstyle\mathbf{q}}_{\Gamma}\cdot\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\Gamma})$ 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 $\mathchoice{\displaystyle\bm{v}}{\textstyle\bm{v}}{\scriptstyle\bm{v}}{\scriptscriptstyle\bm{v}}_{w}$ and the divergence theorem of Eq. 2.5 is applied to the first term, leading to

[TABLE]

Suitable trial and the test function spaces are subspaces of the $\mathcal{H}^{1}(\Gamma)^{3}$ Sobolev space, where $\mathcal{H}^{1}$ 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 $C^{0}$-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 $T$ is the middle surface of the shell and the elements $\tau_{i}\ (i=1,\,\ldots,\,n_{\text{Elem}})$ 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 $N_{j}^{k}(\mathchoice{\displaystyle\bm{x}}{\textstyle\bm{x}}{\scriptstyle\bm{x}}{\scriptscriptstyle\bm{x}}_{\Gamma}(\mathchoice{\displaystyle\bm{r}}{\textstyle\bm{r}}{\scriptstyle\bm{r}}{\scriptscriptstyle\bm{r}}))$ employed for all (physical) fields involved are NURBS of order $k$ with $j=1,\,\ldots,\,n_{c}$ 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 $\nabla_{\Gamma}N(\mathchoice{\displaystyle\bm{x}}{\textstyle\bm{x}}{\scriptstyle\bm{x}}{\scriptscriptstyle\bm{x}}_{\Gamma}(\mathchoice{\displaystyle\bm{r}}{\textstyle\bm{r}}{\scriptstyle\bm{r}}{\scriptscriptstyle\bm{r}}))$ 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 $k$ is now defined by

[TABLE]

where the degrees of freedom $\hat{u}_{j}$ are stored at the control points.

The discrete displacement of the middle surface results as $\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}^{h}=u^{h,i}\mathchoice{\displaystyle\bm{E}}{\textstyle\bm{E}}{\scriptstyle\bm{E}}{\scriptscriptstyle\bm{E}}_{i}$, with $\mathchoice{\displaystyle\bm{E}}{\textstyle\bm{E}}{\scriptstyle\bm{E}}{\scriptscriptstyle\bm{E}}_{i}$ being Cartesian base vectors, with $i=1,\,2,\,3$ and $u^{h,i}=\mathchoice{\displaystyle\bm{N}}{\textstyle\bm{N}}{\scriptstyle\bm{N}}{\scriptscriptstyle\bm{N}}_{\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}}^{\mathrm{T}}\cdot\hat{\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}}^{i}$. In contrast to $\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}$, the difference vector $\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}$ 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 $\mathchoice{\displaystyle\bm{A}}{\textstyle\bm{A}}{\scriptstyle\bm{A}}{\scriptscriptstyle\bm{A}}_{\alpha},\ \alpha=1,\,2$, which are by construction tangential, may be used to define the difference vector $\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}^{h}=w^{h,\,\alpha}\mathchoice{\displaystyle\bm{A}}{\textstyle\bm{A}}{\scriptstyle\bm{A}}{\scriptscriptstyle\bm{A}}_{\alpha}$, where $w^{h,\,\alpha}=\mathchoice{\displaystyle\bm{N}}{\textstyle\bm{N}}{\scriptstyle\bm{N}}{\scriptscriptstyle\bm{N}}_{\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}}^{\mathrm{T}}\cdot\hat{\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}}^{\alpha}$. 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 $\mathchoice{\displaystyle\mathbf{H}}{\textstyle\mathbf{H}}{\scriptstyle\mathbf{H}}{\scriptscriptstyle\mathbf{H}}$, 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 $\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}^{h}=w^{h,i}\mathchoice{\displaystyle\bm{E}}{\textstyle\bm{E}}{\scriptstyle\bm{E}}{\scriptscriptstyle\bm{E}}_{i}$, with $i=1,\,2,\,3$ and $w^{h,i}=\mathchoice{\displaystyle\bm{N}}{\textstyle\bm{N}}{\scriptstyle\bm{N}}{\scriptscriptstyle\bm{N}}_{\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}}^{\mathrm{T}}\cdot\hat{\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}}^{i}$ and the constraint $\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}^{h}\cdot\mathchoice{\displaystyle\bm{n}}{\textstyle\bm{n}}{\scriptstyle\bm{n}}{\scriptscriptstyle\bm{n}}_{\Gamma}=0$ 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 $\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}^{h}=\mathchoice{\displaystyle\mathbf{P}}{\textstyle\mathbf{P}}{\scriptstyle\mathbf{P}}{\scriptscriptstyle\mathbf{P}}\cdot w^{h,i}\mathchoice{\displaystyle\bm{E}}{\textstyle\bm{E}}{\scriptstyle\bm{E}}{\scriptscriptstyle\bm{E}}_{i}$. 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 $\lambda^{h}_{n}=\mathchoice{\displaystyle\bm{N}}{\textstyle\bm{N}}{\scriptstyle\bm{N}}{\scriptscriptstyle\bm{N}}^{\mathrm{T}}_{\lambda_{n}}\cdot\hat{\mathchoice{\displaystyle\bm{\lambda}}{\textstyle\bm{\lambda}}{\scriptstyle\bm{\lambda}}{\scriptscriptstyle\bm{\lambda}}}_{n}$ 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 ${N}_{\lambda_{\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}}}=\{{N}_{{\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}}}\rvert_{\partial\Gamma_{\text{D},\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}}}\}$ and for the difference vector as ${N}_{\lambda_{\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}}}=\{{N}_{{\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}}}\rvert_{\partial\Gamma_{\text{D},\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}}}\}$.

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

[TABLE]

4.1 Discrete weak form

The discrete weak form of the Reissner-Mindlin shell with Lagrange multipliers is the following. Given Young’s modulus $E\in\mathbb{R}^{+}$, Poisson’s ratio $\nu\in[0,\,0.5]$, surface load and moment $\mathchoice{\displaystyle\bm{f}}{\textstyle\bm{f}}{\scriptstyle\bm{f}}{\scriptscriptstyle\bm{f}}^{h},\,\mathchoice{\displaystyle\bm{c}}{\textstyle\bm{c}}{\scriptstyle\bm{c}}{\scriptscriptstyle\bm{c}}^{h}$ on $\Gamma$, traction $\hat{\mathchoice{\displaystyle\bm{p}}{\textstyle\bm{p}}{\scriptstyle\bm{p}}{\scriptscriptstyle\bm{p}}}^{h}$ on $\partial\Gamma_{\text{N},\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}}$, bending moments $\hat{\mathchoice{\displaystyle\bm{m}}{\textstyle\bm{m}}{\scriptstyle\bm{m}}{\scriptscriptstyle\bm{m}}}_{\partial\Gamma}^{h}$ on $\partial\Gamma_{\text{N},\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}}$ and boundary conditions $\hat{\mathchoice{\displaystyle\bm{g}}{\textstyle\bm{g}}{\scriptstyle\bm{g}}{\scriptscriptstyle\bm{g}}}_{\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}}^{h}$ in $\partial\Gamma_{\text{D},\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}}$, $\hat{\mathchoice{\displaystyle\bm{g}}{\textstyle\bm{g}}{\scriptstyle\bm{g}}{\scriptscriptstyle\bm{g}}}_{\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}}^{h}$ on $\partial\Gamma_{\text{D},\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}}$, find the displacement field $\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}^{h}\in\mathcal{S}_{\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}}^{h}$, the difference vector $\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}^{h}\in\mathcal{S}_{\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}}^{h}$, and the Lagrange multiplier fields $(\lambda^{h}_{n},\,\mathchoice{\displaystyle\bm{\lambda}}{\textstyle\bm{\lambda}}{\scriptstyle\bm{\lambda}}{\scriptscriptstyle\bm{\lambda}}_{\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}}^{h},\,\mathchoice{\displaystyle\bm{\lambda}}{\textstyle\bm{\lambda}}{\scriptstyle\bm{\lambda}}{\scriptscriptstyle\bm{\lambda}}_{\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}}^{h})\in\mathcal{L}_{\lambda_{n}}^{h}\times\mathcal{L}_{\mathchoice{\displaystyle\bm{\lambda}}{\textstyle\bm{\lambda}}{\scriptstyle\bm{\lambda}}{\scriptscriptstyle\bm{\lambda}}_{\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}}}^{h}\times\mathcal{L}_{\mathchoice{\displaystyle\bm{\lambda}}{\textstyle\bm{\lambda}}{\scriptstyle\bm{\lambda}}{\scriptscriptstyle\bm{\lambda}}_{\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}}}^{h}$ such that for all test functions $(\mathchoice{\displaystyle\bm{v}}{\textstyle\bm{v}}{\scriptstyle\bm{v}}{\scriptscriptstyle\bm{v}}_{u}^{h},\,\mathchoice{\displaystyle\bm{v}}{\textstyle\bm{v}}{\scriptstyle\bm{v}}{\scriptscriptstyle\bm{v}}_{w}^{h},\,v^{h}_{\lambda_{n}},\,\mathchoice{\displaystyle\bm{v}}{\textstyle\bm{v}}{\scriptstyle\bm{v}}{\scriptscriptstyle\bm{v}}^{h}_{\mathchoice{\displaystyle\bm{\lambda}}{\textstyle\bm{\lambda}}{\scriptstyle\bm{\lambda}}{\scriptscriptstyle\bm{\lambda}}_{\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}}},\,\mathchoice{\displaystyle\bm{v}}{\textstyle\bm{v}}{\scriptstyle\bm{v}}{\scriptscriptstyle\bm{v}}^{h}_{\mathchoice{\displaystyle\bm{\lambda}}{\textstyle\bm{\lambda}}{\scriptstyle\bm{\lambda}}{\scriptscriptstyle\bm{\lambda}}_{\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}}})\in\mathcal{V}_{\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}}^{h}\times\mathcal{V}_{\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}}^{h}\times\mathcal{V}_{\lambda_{n}}^{h}\times\mathcal{V}_{\mathchoice{\displaystyle\bm{\lambda}}{\textstyle\bm{\lambda}}{\scriptstyle\bm{\lambda}}{\scriptscriptstyle\bm{\lambda}}_{\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}}}^{h}\times\mathcal{V}_{\mathchoice{\displaystyle\bm{\lambda}}{\textstyle\bm{\lambda}}{\scriptstyle\bm{\lambda}}{\scriptscriptstyle\bm{\lambda}}_{\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}}}^{h}$, there holds in $\Gamma$

[TABLE]

The usual element assembly yields a linear system of equations (if displacements and rotations are prescribed) in the form

[TABLE]

with $[\underline{\hat{\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}}},\ \underline{\hat{\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}}},\ \hat{\mathchoice{\displaystyle\bm{\lambda}}{\textstyle\bm{\lambda}}{\scriptstyle\bm{\lambda}}{\scriptscriptstyle\bm{\lambda}}}_{n},\ \underline{\hat{\mathchoice{\displaystyle\bm{\lambda}}{\textstyle\bm{\lambda}}{\scriptstyle\bm{\lambda}}{\scriptscriptstyle\bm{\lambda}}}}_{\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}},\ \underline{\hat{\mathchoice{\displaystyle\bm{\lambda}}{\textstyle\bm{\lambda}}{\scriptstyle\bm{\lambda}}{\scriptscriptstyle\bm{\lambda}}}}_{\mathchoice{\displaystyle\bm{w}}{\textstyle\bm{w}}{\scriptstyle\bm{w}}{\scriptscriptstyle\bm{w}}}]^{\mathrm{T}}$ being the sought displacements and rotations of the control points and Lagrange multiplier fields. As usual in the context of Lagrange multiplier methods, Eq. 4.5 has a saddle point structure and the well-known Babuška-Brezzi condition [1, 9, 31] must be satisfied in order to obtain useful solutions in all involved fields.

We have conducted a number of numerical experiments with respect to the orders of the approximations of the displacements, rotations, and the Lagrange multiplier fields enforcing the boundary conditions and constraint on the difference vector. Results presented in the next section are obtained using the same order of these fields, which is implementationally the most convenient setting. Choosing the order for the constraint on the difference vector one order higher or lower did not noticeably change these results. We have also studied the penalty method for this constraint and obtained optimal convergence rates for the penalty parameter $\alpha$ in the range between $10^{6}$ and $10^{10}$.

5 Numerical results

In this section the proposed reformulation of the classical Reissner-Mindlin shell equations in the frame of the TDC is applied to a set of benchmark examples, consisting of the well-known Scordelis-Lo roof from [4], the partly clamped hyperbolic paraboloid from [3, 14] and the clamped flower shaped shell from [56]. In all examples, the shells are rather thin and locking phenomena can be expected, especially in the case of low ansatz orders. However, when increasing the order $p$, locking phenomena decrease significantly and, therefore, no further measures against locking phenomena are considered herein.

In the convergence studies, uniform 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=\{2,\,4,\,8,\,16,\,32,\,64,\,128\}$ elements per side are used. The orders are varied as $p=\{2,\,3,\,4,\,5,\,6\}$.

5.1 Scordelis-Lo roof

The Scordelis-Lo roof is part of the so-called shell obstacle course [4]. The shell is a cylindrical section and is supported with two rigid diaphragms at the ends and loaded by gravity forces, see Fig. 4. The cylinder is defined with a length $L=$50\text{,}, a radius of $R=$25\text{\,} and the angle subtended by the roof is $\phi=$80\text{,}\mathrm{\SIUnitSymbolDegree}. The thickness of the shell is set to $t=$0.25\text{\,}. The material parameters are Young’s modulus $E=$4.32\text{\times}{10}^{8}\text{,} and the Poisson’s ratio $\nu=$0.0\text{\,}. In the convergence study, the maximum vertical displacement $u_{z,\max}$ is compared with the reference solution $u_{z,\max,\text{Ref}}=$0.3024\text{,}$$ as given in [4]. The largest vertical displacement occurs in the midpoint of the free edges $[\pm R\cos($50\text{,}\mathrm{\SIUnitSymbolDegree}$),\,25,\,R\sin($50\text{,}\mathrm{\SIUnitSymbolDegree}$)]^{\mathrm{T}}$.

In Fig. 5(a), the numerical solution of the Scordelis-Lo roof is illustrated. The displacements are magnified by one order of magnitude. The colors on the deformed surface indicate the Euclidean norm of the displacement field $\mathchoice{\displaystyle\bm{u}}{\textstyle\bm{u}}{\scriptstyle\bm{u}}{\scriptscriptstyle\bm{u}}$. In Fig. 5(b), the normalized 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 results improve upon increasing the order of the NURBS. It is noted, that due to the lack of smoothness of the solution for this classical benchmark test (at the corners of the shell), optimal higher-order convergence rates (as typically visualized in double-logarithmic error plots) can not be expected. Hence, the visualization as in Fig. 5(b) follows the usual style of many other references such as, e.g., in [4, 17, 43]. Except for the order $p=2$, the locking phenomena are not very pronounced. Hence, graded meshes as used, e.g., in [44] in order to resolve the boundary layers, are not used here for the sake of simplicity.

5.2 Partly clamped hyperbolic paraboloid

The next test case is a partly clamped hyperbolic paraboloid and is taken from [3, 14]. The shell is defined by $z=x^{2}-y^{2}$ with $(x,\,y)\in[-\nicefrac{{1}}{{2}},\,\nicefrac{{1}}{{2}}]^{2}$, the thickness is set to $t=0.01$ and is loaded by gravity forces, see Fig. 6. The edge at $x=-\nicefrac{{1}}{{2}}$ is clamped and the other three edges are free. The material parameters are Young’s modulus $E=$2.0\text{\times}{10}^{11}\text{,}, Poisson’s ratio $\nu=$0.3\text{\,}. Similar to the example before, the displacements are compared with a reference solution. In particular, the vertical displacement at point $i=(0.5,\,0,\,0.25)^{\mathrm{T}}$ is compared with the reference solution $u_{z,\text{Ref}}=$-9.3355\text{\times}{10}^{-5}\text{,}$$ given in [3].

In Fig. 7(a), the undeformed domain (grey) and the deformed shell is presented, with displacements scaled by a factor of $2000\text{\,}$. In Fig. 7(b), analogously to the example before, the normalized convergence of the vertical displacement $u_{z,i}$ at point $i$ is plotted up to polynomial order of $p=6$ as a function of the element (knot span) size. For the lower orders $p=2,\,3$, the expected locking phenomena is more pronounced compared to the example before. Nevertheless, it is clearly seen that the accuracy for higher-order NURBS increases significantly and the behaviour of convergence is in agreement with the results shown e.g., in [3, 44].

5.3 Clamped flower shaped shell

The geometry for the last example is taken from [56]. The surface is rather complex, but suitable to feature smooth solutions of all physical fields and, consequently, higher-order convergence rates can be achieved. The geometry of the middle surface is given by

[TABLE]

and illustrated in Fig. 8. The shell with the thickness $t=0.1$ is loaded in all three directions and the material parameters are Young’s modulus $E=$10\text{,}, Poisson’s ratio $\nu=$0.3\text{\,}. The middle surface features varying principal curvatures and the curved boundaries are clamped.

Following the same rationale as in [56], the force equilibrium of Eq. 3.21 and moment equilibrium of Eq. 3.22 are computed in strong form and may be called residual errors. In particular, the $L_{2}$-norms of the residual errors are calculated as

[TABLE]

The computation of the residual errors requires second-order surface derivatives which is another argument to employ IGA in this work. The theoretical optimal order of convergence in the residual errors is $\mathcal{O}(p-1)$ due to the presence of second order derivatives.

The numerical solution is presented in Fig. 9(a) and scaled by one order of magnitude. The results of the convergence study are shown in Fig. 9(b) and Fig. 9(c). The polynomial orders are varied up to $p=6$. Due to the complex geometry and boundary layer effects, the pre-asymptotic range is rather pronounced. Therefore, the results of the coarser levels $n=\{2,\,4\}$ are omitted in the results. With a sufficiently small knot span size (element size), the expected higher-order convergence rates are achieved in both residual errors Eq. 5.3 and Eq. 5.4.

The stored elastic energy at the finest level with a polynomial order $p=6$ is $\mathfrak{e}=$5.052,979,16\text{\times}{10}^{-04}\text{,}$$, which may be seen as an overkill solution. This stored elastic energy may be used for future benchmarking when the computation of the residual errors is not desired.

6 Conclusions and Outlook

The linear Reissner-Mindlin shell theory is reformulated in terms of the TDC using a global Cartesian coordinate system. The total deformation is split into the deformation of the middle surface and the rotation of the normal vector. The latter is expressed with a difference vector. The resulting shell equations are applicable to explicitly (parametrized) or implicitly defined shell geometries. Therefore, the proposed reformulation may seen as a generalization of the classical shell theory which is based on local coordinates resulting from a parametrization of the middle surface.

The TDC-based strong form of the force and moment equilibriums is expressed in terms of stress resultants such as membrane forces, bending moments and transverse shear forces using tensor notation. The weak form is consistently derived from the strong form including boundary terms. In the numerical results, the error in the strong form of the shell BVP is integrated over the domain. The resulting residual errors converge with the expected higher orders provided that the involved physical fields are sufficiently smooth. An advantage of investigating the residual errors is that they may be computed for general shell geometries, loadings and boundary conditions without the need for analytical solutions.

The weak form is discretized with an isogeometric approach which may be seen as a realization of the Surface FEM using NURBS as trial and test functions. Despite the fact that standard $C^{0}$-continuous FEM shape functions are sufficient, the high continuity of the isogeometric approach is preferred herein due to the continuous normal vector of the discrete surface, improved convergence properties, and the ability to compute the residual errors. It is clear that in this particular case, the surface mesh, i.e., a NURBS patch, implies a parametrization. This enables a comparison with the classical shell theory and it is seen that the two approaches are equivalent in case of parametrized surfaces, nevertheless, the implementations vary significantly. For a general purpose finite element framework, the TDC-based formulation is beneficial, because significant parts of the implementation needed for shells can be shifted to the underlying finite element technology and may be recycled when other PDEs on surfaces are considered.

In the numerical results, classical and new benchmark tests are presented and optimal higher-order convergence rates are achieved. 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.