The Hellan-Herrmann-Johnson Method for Nonlinear Shells

TL;DR

This paper extends the Hellan-Herrmann-Johnson finite element method from linear Kirchhoff plates to nonlinear shells, enabling accurate modeling of large strains and rotations with improved discretization for complex geometries.

Contribution

It introduces a generalized HHJ method for nonlinear shells, incorporating finite strains, large rotations, and geometric degrees of freedom for kinks.

Findings

Demonstrates effectiveness through benchmark examples

Enables straightforward discretization of structures with kinks

Shows improved modeling of nonlinear shell behavior

Abstract

In this paper we derive a new finite element method for nonlinear shells. The Hellan-Herrmann-Johnson (HHJ) method is a mixed finite element method for fourth order Kirchhoff plates. It uses convenient Lagrangian finite elements for the vertical deflection, and introduces sophisticated finite elements for the moment tensor. In this work we present a generalization of this method to nonlinear shells, where we allow finite strains and large rotations. The geometric interpretation of degrees of freedom allows a straight forward discretization of structures with kinks. The performance of the proposed elements is demonstrated by means of several established benchmark examples.

Tables

Table 1. Table 3.1: Number of degrees of freedom per (hybridized and condensed) element for triangles (T) and quadrilaterals (Q).

	p1 T	p1 Q	p3 T	p3 Q
dof/el	12	16	36	48

Table 2. Table 3.2: Horizontal and vertical deflection of cantilever subjected to end shear force with p1 and 16 × 1 16 1 16\times 1 grid.

$P/P_{max}$	-U	V	$P/P_{max}$	-U	V
0.05	0.026	0.664	0.55	1.811	5.210
0.10	0.104	1.311	0.60	2.007	5.452
0.15	0.225	1.926	0.65	2.195	5.669
0.20	0.382	2.498	0.70	2.375	5.864
0.25	0.565	3.021	0.75	2.546	6.040
0.30	0.765	3.494	0.80	2.710	6.199
0.35	0.974	3.919	0.85	2.867	6.344
0.40	1.187	4.299	0.90	3.015	6.476
0.45	1.399	4.638	0.95	3.157	6.597
0.50	1.608	4.940	1.00	3.292	6.708

Table 3. Table 3.3: Horizontal and vertical deflection of cantilever subjected to end moment for p1 and 16 × 1 16 1 16\times 1 grid.

$M/M_{max}$	$U$	$U_{ex}$	W	$W_{ex}$
0.05	-0.196	-0.196	1.870	1.870
0.10	-0.773	-0.774	3.648	3.648
0.15	-1.698	-1.699	5.249	5.248
0.20	-2.916	-2.918	6.600	6.598
0.25	-4.357	-4.361	7.643	7.639
0.30	-5.942	-5.945	8.338	8.333
0.35	-7.582	-7.585	8.671	8.664
0.40	-9.191	-9.194	8.646	8.637
0.45	-10.687	-10.688	8.291	8.281
0.50	-12.000	-12.000	7.652	7.639
0.55	-13.075	-13.073	6.788	6.775
0.60	-13.875	-13.871	5.772	5.758
0.65	-14.384	-14.377	4.678	4.665
0.70	-14.603	-14.595	3.583	3.571
0.75	-14.556	-14.546	2.556	2.546
0.80	-14.280	-14.270	1.656	1.650
0.85	-13.826	-13.818	0.931	0.926
0.90	-13.254	-13.247	0.407	0.405
0.95	-12.625	-12.621	0.099	0.098
1.00	-12.000	-12.000	0.000	0.000

Table 4. Table 3.4: Vertical deflection at point B 𝐵 B at maximal load for slit annular plate with 10 × 80 10 80 10\times 80 grid.

ref	p1	p3
13.7432	13.8224	13.7772

Table 5. Table 3.5: Radial load-deflection at point B 𝐵 B for the hemispherical shell subjected to alternating radial forces at maximal load for P m ​ a ​ x = 1 subscript 𝑃 𝑚 𝑎 𝑥 1 P_{max}=1 .

ref	p1 1	p1 0.5	p1 0.25
0.093	0.110	0.092	0.0927

Table 6. Table 3.6: Radial load-deflection at point B 𝐵 B for the hemispherical shell subjected to alternating radial forces at maximal load for P m ​ a ​ x = 400 subscript 𝑃 𝑚 𝑎 𝑥 400 P_{max}=400 .

$h$	$2$	$1$	$0.5$	$0.25$
p1	4.1218	3.8811	3.8560	3.8735
p3	3.8319	3.8781	3.8796	3.8796

Table 7. Table 3.7: Deflection U A × 10 3 subscript 𝑈 𝐴 superscript 10 3 U_{A}\times 10^{3} for P x = 10 − 6 subscript 𝑃 𝑥 superscript 10 6 P_{x}=10^{-6} , P z = 0 subscript 𝑃 𝑧 0 P_{z}=0 , and t = 0.0032 𝑡 0.0032 t=0.0032 and W A × 10 3 subscript 𝑊 𝐴 superscript 10 3 W_{A}\times 10^{3} for P x = 0 subscript 𝑃 𝑥 0 P_{x}=0 , P z = 10 − 6 subscript 𝑃 𝑧 superscript 10 6 P_{z}=10^{-6} , and t = 0.0032 𝑡 0.0032 t=0.0032 of twisted beam. Reference values are 5.256 5.256 5.256 and 1.294 1.294 1.294 .

	p1	p3	p1	p3
2x12	5.504	5.460	1.332	1.339
4x24	5.470	5.452	1.337	1.338
6x36	5.460	5.451	1.337	1.338
8x48	5.456	5.451	1.337	1.338

Table 8. Table 3.8: Deflection U A subscript 𝑈 𝐴 U_{A} for P x = 10 − 3 subscript 𝑃 𝑥 superscript 10 3 P_{x}=10^{-3} , P z = 0 subscript 𝑃 𝑧 0 P_{z}=0 , and t = 0.0032 𝑡 0.0032 t=0.0032 and W A subscript 𝑊 𝐴 W_{A} for P x = 0 subscript 𝑃 𝑥 0 P_{x}=0 , P z = 10 − 3 subscript 𝑃 𝑧 superscript 10 3 P_{z}=10^{-3} , and t = 0.0032 𝑡 0.0032 t=0.0032 of twisted beam. Reference values are 4.496 4.496 4.496 and 1.227 1.227 1.227 .

	p1	p3	p1	p3
2x12	4.555	4.536	1.337	1.282
4x24	4.538	4.532	1.279	1.280
6x36	4.535	4.532	1.280	1.280
8x48	4.534	4.531	1.280	1.280

Table 9. Table 3.9: Deflection U A × 10 3 subscript 𝑈 𝐴 superscript 10 3 U_{A}\times 10^{3} for P x = 1 subscript 𝑃 𝑥 1 P_{x}=1 , P z = 0 subscript 𝑃 𝑧 0 P_{z}=0 , and t = 0.32 𝑡 0.32 t=0.32 and W A × 10 3 subscript 𝑊 𝐴 superscript 10 3 W_{A}\times 10^{3} for P x = 0 subscript 𝑃 𝑥 0 P_{x}=0 , P z = 1 subscript 𝑃 𝑧 1 P_{z}=1 , and t = 0.32 𝑡 0.32 t=0.32 of twisted beam. Reference values are 5.424 5.424 5.424 and 1.754 1.754 1.754 .

	p1	p3	p1	p3
2x12	5.654	5.598	1.933	1.798
4x24	5.605	5.591	1.822	1.795
6x36	5.597	5.590	1.806	1.795
8x48	5.593	5.589	1.801	1.795

Table 10. Table 3.10: Deflection U A subscript 𝑈 𝐴 U_{A} for P x = 10 3 subscript 𝑃 𝑥 superscript 10 3 P_{x}=10^{3} , P z = 0 subscript 𝑃 𝑧 0 P_{z}=0 , and t = 0.32 𝑡 0.32 t=0.32 and W A subscript 𝑊 𝐴 W_{A} for P x = 10 3 subscript 𝑃 𝑥 superscript 10 3 P_{x}=10^{3} , P z = 1 subscript 𝑃 𝑧 1 P_{z}=1 , and t = 0.32 𝑡 0.32 t=0.32 of twisted beam. Reference values are 4.610 4.610 4.610 and 1.778 1.778 1.778 .

	p1	p3	p1	p3
2x12	4.661	4.621	1.908	1.789
4x24	4.628	4.618	1.810	1.786
6x36	4.622	4.617	1.796	1.785
8x48	4.619	4.617	1.791	1.785

Table 11. Table 3.11: Membrane stress Σ x ​ x × 10 8 subscript Σ 𝑥 𝑥 superscript 10 8 \Sigma_{xx}\times 10^{8} of Z-section cantilever at maximal load.

	[18]	p1	p3
8x6	$-0.953$	$-0.7620$	$-1.0929$
32x15	$-1.063$	$-1.0777$	$-1.0933$
64x30	-	$-1.0989$	$-1.0933$

Equations

$$H^1({\mathcal{T}}_h):=\left\{u\in L^2({\mathcal{S}}_h)|{\nabla }_{\tau }u\in [L^2({\mathcal{S}}_h){]}^3\right\},$$

$$\Sigma ({\mathcal{T}}_h):=\left\{\mathbf{\sigma }\in [C^{\infty }({\mathcal{T}}_h){]}_{sym}^{3\times 3}|[[{\mathbf{\sigma }}_{\mu \mu }]]=0\right\},$$

$$V_h^k({\mathcal{T}}_h):={\Pi }^k({\mathcal{T}}_h)\cap C^0({\mathcal{S}}_h),$$

$${\Sigma }_h^k({\mathcal{T}}_h):=\left\{\mathbf{\sigma }\in [{\Pi }^k({\mathcal{T}}_h){]}_{sym}^{3\times 3}|[[{\mathbf{\sigma }}_{\mu \mu }]]=0\right\},$$

$${\overline{\Gamma }}_h^k({\mathcal{T}}_h):=\left\{u\in [{\Pi }^k({\mathcal{T}}_h){]}^3|[[u_{\mu }]]=0\right\},$$

$$\hat{\Omega }:=\{\hat{x}+z{\hat{\nu }}_{\hat{\mathcal{S}}}(\hat{x}):\hat{x}\in \hat{\mathcal{S}},z\in [-t/2,t/2]\}.$$

$$\mathcal{W}=\frac{t}{8}\Vert I-\hat{I}{\Vert }_{\mathbf{M}}^2+\frac{t^3}{24}\Vert II-\widehat{II}{\Vert }_{\mathbf{M}}^2.$$

$$\Vert \cdot {\Vert }_{\mathbf{M}}^2:=\frac{E}{1-{\nu }^2}{\int }_{{\hat{\mathcal{S}}}_h}\left(\nu \text{tr}(\cdot {)}^2+(1-\nu )\text{tr}({\cdot }^2)\right)dx,$$

$$E_{\text{mem}}:=\frac{t}{8}\Vert I-\hat{I}{\Vert }_{\mathbf{M}}^2=\frac{t}{2}\Vert \mathbf{E}{\Vert }_{\mathbf{M}}^2.$$

$${\tilde{E}}_{\text{bend}}$$

$$=\frac{t^3}{24}\Vert {\mathbf{F}}^T{\nabla }_{\hat{\tau }}(\nu \circ \varphi )-{\nabla }_{\hat{\tau }}\hat{\nu }{\Vert }_{\mathbf{M}}^2.$$

$$E_{\text{bend}}:=$$

$$\displaystyle+\sum_{\hat{E}\in\hat{\mathcal{E}}_{h}}\|\sphericalangle(\nu_{L},\nu_{R})\circ\phi-\sphericalangle(\hat{\nu}_{L},\hat{\nu}_{R})\|_{\bm{M},\hat{E}}^{2}\big{)}.$$

$$\tilde{\mathcal{W}}(u):=E_{\text{mem}}+E_{\text{bend}}.$$

$$\tilde{\mathcal{L}}(u,\mathbf{\sigma }):=\frac{t}{2}\Vert \mathbf{E}{\Vert }_{\mathbf{M}}^2-\frac{6}{t^3}\Vert \mathbf{\sigma }{\Vert }_{{\mathbf{M}}^{-1}}^2+\tilde{B}(\mathbf{\sigma },u),$$

$$\tilde{B}(\mathbf{\sigma },u):=$$

$$-\sum _{\hat{E}\in {\hat{\mathcal{E}}}_h}⟨\sphericalangle ({\nu }_L,{\nu }_R)\circ \varphi -\sphericalangle ({\hat{\nu }}_L,{\hat{\nu }}_R),{\mathbf{\sigma }}_{\hat{\mu }\hat{\mu }}{⟩}_{\hat{E}},$$

$$\mathcal{L}(u,\mathbf{\sigma })=\frac{t}{2}\Vert \mathbf{E}{\Vert }_{\mathbf{M}}^2-\frac{6}{t^3}\Vert \mathbf{\sigma }{\Vert }_{{\mathbf{M}}^{-1}}^2-B(\mathbf{\sigma },u),$$

$$B(\mathbf{\sigma },u)=$$

$$\sum _{\hat{T}\in {\hat{\mathcal{T}}}_h}{\int }_{\hat{T}}\mathbf{\sigma }:({\mathcal{H}}_{\nu \circ \varphi }+(1-\hat{\nu }\cdot \nu \circ \varphi ){\nabla }_{\hat{\tau }}\hat{\nu })dx$$

$$-\sum _{\hat{E}\in {\hat{\mathcal{E}}}_h}{\int }_{\hat{E}}(\sphericalangle ({\hat{\nu }}_L,{\hat{\nu }}_R)-\sphericalangle ({\nu }_L,{\nu }_R)\circ \varphi ){\mathbf{\sigma }}_{\hat{\mu }\hat{\mu }}ds.$$

$$\nu \circ \varphi =\frac{1}{\Vert \text{cof}(\mathbf{F})\hat{\nu }{\Vert }_2}\text{cof}(\mathbf{F})\hat{\nu }=\frac{1}{J}\text{cof}(\mathbf{F})\hat{\nu },$$

$${\tau }_e\circ \varphi =\frac{1}{\Vert \mathbf{F}{\hat{\tau }}_e{\Vert }_2}\mathbf{F}{\hat{\tau }}_e=\frac{1}{J_b}\mathbf{F}{\hat{\tau }}_e,$$

$$\mu \circ \varphi =\pm (\nu \circ \varphi )\times ({\tau }_e\circ \varphi )=\frac{1}{\Vert ({\mathbf{F}}^{†}{)}^T\hat{\mu }{\Vert }_2}({\mathbf{F}}^{†}{)}^T\hat{\mu },$$

$$\Vert \cdot {\Vert }_{{\mathbf{M}}^{-1}}^2:=\frac{1+\nu }{E}{\int }_{{\hat{\mathcal{S}}}_h}(\text{tr}({\cdot }^2)-\frac{\nu }{2\nu +1}\text{tr}(\cdot {)}^2)dx,$$

$$B(\mathbf{\sigma },u)=$$

$$-\sum _{\hat{E}\in {\hat{\mathcal{E}}}_h}{\int }_{\hat{E}}\sphericalangle ({\nu }_L,{\nu }_R)\circ \varphi {\mathbf{\sigma }}_{\hat{\mu }\hat{\mu }}ds.$$

$$\frac{1}{2}\sphericalangle ({\nu }_L,{\nu }_R)=\{\nu \}\cdot {\mu }_L+\mathcal{O}(|\{\nu \}\cdot {\mu }_L{|}^3),$$

$$\mathcal{L}(u,\mathbf{\sigma },\hat{\alpha })=\frac{t}{2}\Vert \mathbf{E}{\Vert }_{\mathbf{M}}^2-\frac{6}{t^3}\Vert \mathbf{\sigma }{\Vert }_{{\mathbf{M}}^{-1}}^2-B(\mathbf{\sigma },u,\hat{\alpha }),$$

$$B(\mathbf{\sigma },u,\hat{\alpha })$$

$$=\sum _{\hat{T}\in {\hat{\mathcal{T}}}_h}{\int }_{\hat{T}}\mathbf{\sigma }:({\mathcal{H}}_{\nu \circ \varphi }+(1-\hat{\nu }\cdot \nu \circ \varphi ){\nabla }_{\hat{\tau }}\hat{\nu })dx$$

$$-\sum _{\hat{E}\in {\hat{\mathcal{E}}}_h}{\int }_{\hat{E}}(\sphericalangle ({\hat{\nu }}_L,{\hat{\nu }}_R)-\sphericalangle ({\nu }_L,{\nu }_R)\circ \varphi )⟨⟨{\mathbf{\sigma }}_{\hat{\mu }\hat{\mu }}⟩⟩ds$$

$$+{\int }_{\hat{E}}{\hat{\alpha }}_{\hat{\mu }}[[{\mathbf{\sigma }}_{\hat{\mu }\hat{\mu }}]]ds,$$

$$\sum _{\hat{E}\in {\hat{\mathcal{E}}}_h}{\int }_{\hat{E}}\sphericalangle ({\hat{\nu }}_L,{\hat{\nu }}_R)-\sphericalangle ({\nu }_L,{\nu }_R)\circ \varphi ds=$$

$$\sum _{\hat{T}\in {\hat{\mathcal{T}}}_h}{\int }_{\partial \hat{T}}\sphericalangle (\{\hat{\nu }\},\hat{\mu })-\sphericalangle (\{\nu \},\mu )\circ \varphi ds.$$

$$\{\nu \}\approx \frac{1}{\Vert {\mathbf{P}}_{{\tau }_e}^{\perp }(\{\nu {\}}^n){\Vert }_2}{\mathbf{P}}_{{\tau }_e}^{\perp }(\{\nu {\}}^n)=:{\overline{\{\nu \}}}^n.$$

$$\text{div}(\text{div}({\nabla }^{2}w))=f,$$

$$\mathcal{L}(w,\mathbf{\sigma }):=$$

Full text

Keywords: nonlinear shells; structural mechanics; discrete differential geometry; mixed finite elements; Kirchhoff hypothesis

1 Introduction

The difficulty of constructing simple $C^{1}$-conforming Kirchhoff–Love shell elements led to the development of the well-known discrete Kirchhoff (DKT) elements [30, 46, 3], where the Kirchhoff constraint was inforced in a discrete way along the edges. The class of rotation-free (RF) elements eliminate the rotational degrees of freedom by using out-of-plane translation degrees of freedom (dofs) [33, 9, 19]. Alternative approaches are discontinuous Galerkin (DG) methods [17, 21, 47] and Isogeometric Analysis (IGA) [26, 40, 28, 16].

The HHJ method for fourth order Kirchhoff plates has been developed and analyzed in [22, 23, 27]. Later work has been done in the 80s [13, 1], 90s [44] and recently after 20 years [11, 25, 7]. It overcomes the issue of $C^{1}$-conformity by introducing the moment tensor as an additional tensor field leading to a mixed method. The tangential displacement and normal-normal stress method (TDNNS) developed for linear elasticity and Reissner–Mindlin plates in [42, 34, 35, 36] follows the idea of mixed methods, where the stress tensor gets interpolated in the reinvented $H(\text{divdiv})$ space from the HHJ method.

In this paper modern coordinate-free differential geometry, see e.g. [14, 43], is used to define the shell energy. The aim of this work is to find a (high-order) finite element shell element, consisting of $H^{1}$-conforming finite elements for the displacement and $H(\text{divdiv})$ elements for the moments. It turns out, that this model can be seen as a generalization of the HHJ method to nonlinear shells. Furthermore, the method can handle surfaces with kinks in a natural way without additional treatment. Numerical results are shown to confirm the model.

2 Methodology

2.1 Notation and finite element spaces

Let $\mathcal{S}$ be a 2-dimensional surface in $\mathbb{R}^{3}$, and let $\mathcal{S}_{h}=\bigcup_{T\in\mathcal{T}_{h}}T$ be its approximation by a triangulation $\mathcal{T}_{h}$ consisting of possibly curved triangles or quadrilaterals. The set of all edges in $\mathcal{T}_{h}$ is denoted by $\mathcal{E}_{h}$. Further, let $L^{2}(\mathcal{S}_{h})$ and $C^{0}(\mathcal{S}_{h})$ be the set of all square-integrable and continuous functions on $\mathcal{S}_{h}$, respectively.

For each element in $\mathcal{T}_{h}$ we denote the surface normal vector by $\nu$ and the normalized edge tangent vector between elements by $\tau_{e}$. The outgoing element-normal vector $\mu$ is defined as $\mu=\pm\nu\times\tau_{e}$ depending on the orientation of $\tau_{e}$, cf. Figure 2.1.

The set of all piece-wise polynomials of degree $k$ on $\mathcal{T}_{h}$ is denoted by $\Pi^{k}(\mathcal{T}_{h})$. With this, we define the following function and finite element spaces

[TABLE]

where we used the notations $\bm{\sigma}_{\mu\mu}:=\mu^{T}\bm{\sigma}\mu$ and $u_{\mu}:=u\cdot\mu$ with $\llbracket\cdot\rrbracket$ denoting the jump over elements. Note, that $\nabla_{\tau}u$ denotes the surface gradient of $u$, which can be introduced in weak sense [15], or directly as Fréchet-derivative.

For the construction of finite element spaces and explicit basis functions of (2.3), (2.4) and (2.5) we refer to [6, 42, 48, 32, 8, 49]. With a hierarchical basis of (2.5) we can define the finite element space $\Gamma_{h}^{k}(\mathcal{T}_{h})$ as the space $\overline{\Gamma}_{h}^{k}(\mathcal{T}_{h})$, where the inner degrees of freedom are neglected. $\Sigma_{h}^{k}(\mathcal{T}_{h})$ will be called the Hellan–Herrmann–Johnson finite element space.

In the following we denote the Frobenius scalar product of two matrices $\bm{A}$, $\bm{B}$ by $\bm{A}:\bm{B}:=\sum_{i,j}\bm{A}_{ij}\bm{B}_{ij}$, $\|\bm{A}\|_{F}:=\sqrt{\bm{A}:\bm{A}}$ and $\sphericalangle(a,b):=\arccos(\frac{a\cdot b}{\|a\|_{2}\|b\|_{2}})$ measures the angle between two vectors $a$, $b$, $\|\cdot\|_{2}$ denoting the Euclidean norm.

2.2 Shell model

Let $\hat{\Omega}\subset\mathbb{R}^{3}$ be an undeformed configuration of a shell with thickness $t$, described by the mid-surface $\hat{\mathcal{S}}$ and the according orientated normal vector $\hat{\nu}_{\hat{\mathcal{S}}}$

[TABLE]

Furthermore, let $\Phi:\hat{\Omega}\rightarrow\Omega$ be the deformation from the initial to the deformed configuration of the shell and $\phi:\hat{\mathcal{S}}_{h}\rightarrow\mathcal{S}_{h}$ the deformation of the approximated mid-surface. I.e., let $\phi\in[V_{h}^{k+1}(\hat{\mathcal{T}}_{h})]^{3}$ with $\hat{\mathcal{T}}_{h}$ and $\mathcal{T}_{h}=\phi(\hat{\mathcal{T}}_{h})$ the according triangulations of $\hat{\mathcal{S}}_{h}$ and $\mathcal{S}_{h}$. Then, we define $\bm{F}:=\nabla_{\hat{\tau}}\phi$ and $J:=\|\,\text{cof}(\bm{F})\|_{F}=\|\,\text{cof}(\bm{F})\hat{\nu}\|_{2}$ as the deformation gradient and the deformation determinant, respectively. Here, $\,\text{cof}(\bm{F})$ denotes the cofactor matrix of $\bm{F}$. We can split the deformation into the identity function and the displacement, $\phi=\text{id}+u$, and thus, $\bm{F}=\bm{P}_{\hat{\tau}}+\nabla_{\hat{\tau}}u$ with the projection onto the tangent plane $\bm{P}_{\hat{\tau}}:=\bm{I}-\hat{\nu}\otimes\hat{\nu}$, $\otimes$ denoting the dyadic outer product.

We consider the Kirchhoff–Love assumption, where the deformed normal vector has to be orthogonal to the deformed mid-surface $\mathcal{S}_{h}$. With Steiner’s formula, asymptotic analysis in the thickness parameter $t$ and using the plane strain assumption for the material norm, we obtain for the according shell energy functional

[TABLE]

$\mathcal{W}$ is given in terms of differential forms, see [41], and (2.7) is comparable to the classical formulations [12, 5, 10].

The material norm is given by

[TABLE]

with $E$ the Young’s modulus and $\nu$ the Poisson’s ratio. $\hat{{I}}$, ${I}$ and $\hat{{II}}$, ${II}$ denote the (pull-backed) first and second fundamental form of the reference and deformed configuration, respectively. With the Green strain tensor $\bm{E}:=1/2(\bm{C}-\bm{P}_{\hat{\tau}})$ restricted on the tangent space, $\bm{C}=\bm{F}^{T}\bm{F}$ denoting the Cauchy-Green tensor, we obtain

[TABLE]

This corresponds to the membrane energy of the shell. The difference between the curvature of the deformed and initial second fundamental form describes the bending energy, for which holds

[TABLE]

Motivated by discrete differential geometry, see [20] and references therein, and DG methods [2] we add also distributional contributions to the bending energy

[TABLE]

Thus, with the notation of (2.9) and (2.11), we have to minimize

[TABLE]

To reduce this fourth order problem to a second order one, we introduce a new variable $\bm{\sigma}$ which leads to a mixed saddle point problem. Hence, we have to find the critical points of the following Lagrange functional, which is equivalent to minimize (2.12), see Appendix A,

[TABLE]

where

[TABLE]

with $\langle\cdot,\cdot\rangle$ denoting the $L^{2}$-scalar product on an element $\hat{T}$ or on an edge $\hat{E}$. With some computations, see Appendix A, we finally obtain the following Lagrange functional

[TABLE]

with

[TABLE]

In (2.16) $\bm{\mathcal{H}}_{\nu\circ\phi}:=\sum_{i}(\nabla_{\hat{\tau}}^{2}u_{i})\nu_{i}\circ\phi$, where $\nabla_{\hat{\tau}}^{2}$ denotes the surface Hessian [15]. For the deformed normal and tangent vectors the following identities hold

[TABLE]

where $\bm{F}^{\dagger}$ denotes the Moore–Penrose pseudo-inverse of $\bm{F}$.

The Lagrange multiplier $\bm{\sigma}$ has the physical meaning of the moment. Note, that the thickness parameter $t$ appears now also in the denominator and the inverse material tensor

[TABLE]

is used.

In case of a flat plane (2.16) becomes

[TABLE]

A possible simplification of (2.16) can be achieved by the approximation

[TABLE]

where $\{\nu\}:=\frac{1}{\|\nu_{L}+\nu_{R}\|_{2}}(\nu_{L}+\nu_{R})$ denotes the averaged normal vector.

The resulting system is a saddle point problem, which would lead to an indefinite matrix after assembling. To overcome this problem, we can use complete discontinuous elements for the moment $\bm{\sigma}$ and introduce a hybridization variable $\hat{\alpha}\in\Gamma_{h}^{k}(\hat{\mathcal{T}}_{h})$ to reinforce the normal-normal continuity of $\bm{\sigma}$:

[TABLE]

where (2.16) is now given by

[TABLE]

with $\langle\langle\bm{\sigma}_{\hat{\mu}\hat{\mu}}\rangle\rangle:=1/2(\bm{\sigma}_{\hat{\mu}_{L}\hat{\mu}_{L}}+\bm{\sigma}_{\hat{\mu}_{R}\hat{\mu}_{R}})$. Due to the hybridization variable $\hat{\alpha}$, we can use static condensation to eliminate the moment $\bm{\sigma}$ locally, which leads to a positive definite problem again. The new unknown $\hat{\alpha}$ has the physical meaning of the changed angle, the rotation, between two elements.

For the computation of the jump term we use that

[TABLE]

To compute the deformed averaged normal vector $\{\nu\}$ on an edge, information of the two neighbored elements is needed at once, which would need e.g. Discontinuous Galerkin techniques. Instead, one can use the information of the last (load-step) solution $\{\nu\}^{n}$, see Figure 2.2. To measure the correct angle, we have to project $\{\nu\}^{n}$ to the plane orthogonal to the tangent vector $\tau_{e}$ by using the projection $\bm{P}_{\tau_{e}}^{\perp}=\bm{I}-\tau_{e}\otimes\tau_{e}$, and then re-normalize it

[TABLE]

Note that $\tau_{e}$ itself depends on the unknown deformation. By using (2.26) we have to ensure that $\{\nu\}^{n}$ lies between the two element-normal vectors, see Figure 2.2. For smooth manifolds the angle between the element-normal vectors tends to $180$ degree as $h\rightarrow 0$. Hence, this assumption is fulfilled, if the elements do not rotate more than half of their included angle during one load-step, which is an acceptable and realistic assumption.

2.3 Relation to the HHJ-method

If we assume to have a plate which lies in the x-y-plane and a force $f$ is acting orthogonal on it, we can compute the linearized bending energy by solving the following fourth order scalar equation on $\hat{\mathcal{S}}_{h}$

[TABLE]

where the thickness $t$ and all material parameters are hidden in the right-hand side $f$.

Therefore, the HHJ-method [22, 23, 27] introduces the linearized moment tensor $\bm{\sigma}$ and solves the following saddle point problem instead, given by the Lagrange functional

[TABLE]

If we now consider our shell model (2.15), neglect the membrane energy term and the material parameters and linearize the bending energy, see Appendix B, we obtain (2.3). Thus, (2.15) can be seen as a generalization of the HHJ-method (2.3) from linear plates to nonlinear shells.

2.4 Boundary conditions and kink structures

For $H^{1}$ the Dirichlet boundary condition $u=u_{D}$ can be used to prescribe the displacement on the boundary, whereas the do-nothing condition is used for free boundaries. For $\bm{\sigma}\in H(\text{divdiv})$ we can prescribe the normal-normal component, $\bm{\sigma}_{\mu\mu}$, on the boundary. Homogeneous Dirichlet data, $\bm{\sigma}_{\mu\mu}=0$, are used for free boundaries. By setting non-homogeneous data one can prescribe a moment. The do-nothing Neumann boundary condition $\bm{\sigma}_{\mu\tau_{e}}=0$ is used for clamped boundaries.

In the case of a complete discontinuous moment tensor and the hybridization variable $\hat{\alpha}$, the boundary conditions for $\bm{\sigma}$ have to be incorporated in terms of $\hat{\alpha}$. Note that the essential and natural boundary conditions swap, i.e. the clamped boundary condition is now set directly as homogeneous Dirichlet data and the prescribed moment is handled natural as a right-hand side.

If we compute the variations of (2.15) with respect to $\bm{\sigma}$, we obtain in strong form that the angle from the initial configuration gets preserved, see (C.2). The hidden interface condition for the displacement $u$ in strong form are not needed for the method itself. However, if one uses e.g. Residual error estimators, the boundary conditions are crucial, see Appendix C for the calculations.

The method can also handle non-smooth surfaces with kinks and branching shells, where one edge is shared by more than two elements, in a natural way, without any extra treatments. Due to the normal-normal continuity of $\bm{\sigma}$ the moment gets preserved over the kinks and as the angle is the same on the initial and deformed configuration, the kink itself gets also preserved. Note, that in this case simplification (2.22) cannot be used any more, as $|\{\nu\}\cdot\mu|\nrightarrow 0$ as $h\rightarrow 0$ at the kinks.

2.5 Shell element

Combining the displacement $u\in[V_{h}^{k}(\hat{\mathcal{T}}_{h})]^{3}$, the moment tensor $\bm{\sigma}\in\Sigma^{k-1}_{h}(\hat{\mathcal{T}}_{h})$ and, eventual, the hybridization space $\Gamma_{h}^{k-1}(\hat{\mathcal{T}}_{h})$ leads to our shell element. For polynomial order $k$, the method will be denoted by pk, i.e. p1 is the lowest order method consisting of piece-wise linear displacements and piece-wise constant moments. In Figure 2.3 the hybridized p1 and p2 element for quadrilaterals can be seen. Note, that the hybridized lowest order triangle shell element is equivalent to the Morley element [30]. If we use the lowest order elements on triangles for (2.15) then the Hessian term vanishes, as only linear polynomials are used. For quadrilaterals the Hessian is constant on each element in this case.

To solve (2.15) we have to assemble the according matrix. As it is formulated in terms of a Lagrange functional, the first variations must be computed, which is a bit challenging due to the non-linearity but doable, see Appendix C. If, however, the finite element software supports energy based integrators where the variations are calculated automatically, one can use directly the Lagrange functional (2.15).

2.6 Membrane locking

We observed that the lowest order elements do not suffer from locking, but for the higher order methods membrane locking, cf. [37], may occur, e.g. in the benchmark cantilever subjected to end moment, section 3.2. To overcome this problem one can interpolate the membrane stress tensor by a $L^{2}$-projection into a space of reduced dimension, $\|\mathcal{I}^{h}_{L^{2}}\bm{E}\|_{\bm{M}}^{2}$. The projection can be incorporated to (2.15) by introducing an auxiliary variable $\bm{R}$ and adding for the displacement $u\in[V^{k}_{h}(\hat{\mathcal{T}}_{h})]^{3}$ and $\bm{R}\in[\Pi^{k-1}(\hat{\mathcal{T}}_{h})]^{3\times 3}_{sym}$

[TABLE]

to the Lagrange functional. As $\bm{R}$ is discontinuous, we can use static condensation to eliminate it locally. This works well for structured quadrilateral meshes and is similar to reduced integration order methods. For triangles, however, the locking is reduced, but still has an impact to the solution. Here, other interpolation operators and spaces have to be used, which is topic of further research.

3 Numerical results

The method is implemented in the NGS-Py interface, which is based on the finite element library Netgen/NGSolve222www.ngsolve.org [38, 39].

We will use the lowest order elements p1and also the p3 method as an high-order example. Table 3.1 lists the according number of degrees of freedom for each element.

3.1 Cantilever subjected to end shear force

An end shear force $P$ on the right boundary is applied to a cantilever, which is fixed on the left. The material and geometrical properties are $E=1.2\times 10^{6}$, $\nu=0$, $L=10$, $W=1$, $t=0.1$ and $P_{max}=4$, see Figure 3.1. A structured $16\times 1$ rectangular grid is used. The reference values are taken from [45]. In Figure 3.2 one can see the initial and deformed mesh and in Figure 3.3 and Table 3.2 the results.

3.2 Cantilever subjected to end moment

A cantilever is clamped on the left side and a moment $M$ is applied on the right. On the other boundaries we use the symmetry-condition. The material and geometrical properties are $E=1.2\times 10^{6}$, $\nu=0$, $L=12$, $W=1$, $t=0.1$ and $M_{max}=50\pi/3$, see Figure 3.4. The results can be found in Figure 3.6 and Table 3.3, and the initial and final mesh in Figure 3.5.

3.3 Slit annular plate

The material and geometrical properties are $E=2.1\times 10^{8}$, $\nu=0$, $R_{i}=6$, $R_{o}=10$, $t=0.03$ and $P_{max}=4.034$, see Figure 3.7. We used structured quadrilateral meshes. The quantity of interest is the transverse displacement at point B. The reference value of $13.7432$ is taken from [24]. The initial and deformed mesh can be seen in Figure 3.8 and the results in Figure 3.9 and Table 3.4.

3.4 Hemispherical shell subjected to alternating radial forces

The material and geometrical properties are $E=6.825\times 10^{7}$, $\nu=0.3$, $R=10$, $t=0.04$, see Figure 3.10. A non-structured triangulation is used with different mesh-sizes. For $P_{max}=1$ [41] gives the reference value of the vertical deflection at point $B$ with $0.093$ at maximal load. In Table 3.5 the results for p1 and three different meshes can be found. For the large displacement case we used $P_{max}=400$, see Figure 3.11 and 3.12. The results shown in Table 3.6 are convenient with [24].

3.5 Twisted beam

A beam is twisted by $90$ degrees and clamped on the left side, whereas a point load is applied on the middle of the right boundary. The material and geometrical properties are $E=2.9\times 10^{7}$, $\nu=0.22$, $L=12$, $b=1.1$, $t=0.0032$, $0.32$, see Figure 3.13.

Different forces, $P_{max}\in\{10^{-6}$, $10^{-3}$, $1$, $10^{3}\}$, are applied in x- and z-direction. Some combinations of thickness and force parameters led to a solution in a linear regime, see Table 3.7 and 3.9, where the reference solutions are taken from [4] and [29], respectively. Others are already in the nonlinear regime, see Table 3.8 and 3.10. Therefore, the full three-dimensional model is used with a $150\times 14\times 2$ structured cubic grid and standard Lagrangian elements of polynomial order 3, i.e. 162 dofs/cube, to generate a reference solution.

3.6 Z-section cantilever

A moment $M=1.2\times 10^{6}$ is applied at the right end of a Z-section, which is fixed on the left side. Therefore, two shear forces $P=6\times 10^{5}$ are involved, see Figure 3.14. The material and geometrical properties are $E=2.1\times 10^{11}$, $\nu=0.3$, $t=0.1$, $L=10$, $W=2$ and $H=1$. The quantity of interest is the membrane stress $\Sigma_{xx}$ at point $A$. The reference value $-1.08\times 10^{8}$ is taken from NAFEMS [31]. The results are compared with rotation-free elements [18] and can be found in Table 3.11.

3.7 T-section cantilever

We propose an example where more than two elements share an edge. The material and geometrical properties are $E=6\times 10^{6}$, $\nu=0$, $t=0.1$, $L=1$, $W=1$ and $H=1$. The structure is clamped on the bottom and a shear force $P_{max}=1000$ is applied on the left boundary, see Figure 3.15.

The moment induced by the shear force $P$ on the left top branch goes over the kink to the bottom branch where the structure is fixed without inducing moments on the right top one. Thus, it only rotates and the curvature is zero also after the deformation. The deflections of the point $A$ are given in Figure 3.16.

Acknowledgments

Michael Neunteufel has been funded by the Austrian Science Fund (FWF) project W1245.

Appendix A Lagrange functional

We compute the variations of the Lagrange functional in (2.13), neglecting the sums over $\hat{T}$ and $\hat{E}$,

[TABLE]

Expressing $\bm{\sigma}$ from (A.1) and inserting it into (A.2) yields to the same expression as the variation of (2.12) with respect to the displacement $u$. We conclude that (2.13) and (2.12) are equivalent.

Equivalence of (2.15) and (2.13) follows by differentiating the identity $\bm{F}^{T}\nu\circ\phi=0$ and some computations

[TABLE]

where $\bm{H}_{i}:=\nabla_{\hat{\tau}}^{2}u_{i}+\nabla_{\hat{\tau}}((\bm{P}_{\hat{\tau}})_{i})$, $(\bm{P}_{\hat{\tau}})_{i}$ denoting the i-th column of $\bm{P}_{\hat{\tau}}$ and $\partial_{\hat{x}_{i}}\bm{F}$ the i-th partial derivative of $\bm{F}$. With $\bm{P}_{\hat{\tau}}=\bm{I}-\hat{\nu}\otimes\hat{\nu}$, neglecting $\phi$, and sum convention for $i$ we obtain

[TABLE]

where we used that $\nabla_{\hat{\tau}}\hat{\nu}_{i}\otimes\hat{\nu}:\bm{\sigma}\equiv 0$.

Appendix B Linearization

To show that (2.15) simplifies to (2.3) in the linear regime we use that the gradient of the displacement of the full three-dimensional model $\nabla U=\nabla(u+z\nu\circ\phi)$ is small, $\nabla U=\mathcal{O}(\varepsilon)\ll 1$. Thus, we immediately obtain that $\nabla_{\hat{\tau}}u=\mathcal{O}(\varepsilon)$, $\bm{F}=\bm{I}+\mathcal{O}(\varepsilon)$, $J_{b}=1+\mathcal{O}(\varepsilon)$ and $\bm{\sigma}=\mathcal{O}(\varepsilon)$. Furthermore, there holds $\nu\circ\phi-\hat{\nu}=-\hat{\nu}^{T}\nabla_{\hat{\tau}}u+\mathcal{O}(\varepsilon^{2})$ for $\varepsilon\rightarrow 0$. We neglect all terms of order $\mathcal{O}(\varepsilon^{2})$ or higher. For simplicity we will also neglect the $\phi$ dependency, e.g. we write $\nu$ instead of $\nu\circ\phi$.

Starting from (2.21), we obtain on each $\hat{T}\in\hat{\mathcal{T}}_{h}$

[TABLE]

For the jump term we use (2.22) and (2.26), such that

[TABLE]

For ease of presentation we neglect $\bm{\sigma}_{\hat{\mu}\hat{\mu}}$ in (B.2), employ that $\frac{1}{\|\bm{P}_{\tau_{e}}^{\perp}(\{\hat{\nu}\})\|_{2}}\bm{P}_{\tau_{e}}^{\perp}(\{\hat{\nu}\})=\bm{P}_{\tau_{e}}^{\perp}(\{\hat{\nu}\})+\mathcal{O}(\varepsilon^{2})$ and that $\{\hat{\nu}\}=\hat{\nu}$ on a flat plane to obtain

[TABLE]

For $\int_{\partial\hat{T}}\{\nu\}\cdot\mu\,ds$ the linearization is done analogously and leads to the same result as (B.3).

If we now use (B.1), (B.2), (B.3) and (2.21) and apply it to (2.15), neglect the membrane energy term and the constants, and employ that $\hat{\nu}=\begin{pmatrix}0\\ 0\\ 1\end{pmatrix}$, we finally obtain

[TABLE]

which is indeed (2.3).

Appendix C Variations

We compute the variations of (2.15) to deduce the bilinear form of the according variational equations. Then we will (partly) integrate by parts to find the hidden boundary conditions in strong form.

For simplicity, we will neglect the material tensor $\bm{M}$ and write only $\nu$ instead of $\nu\circ\phi$. The same holds for $\mu$ and $\tau_{e}$. We will consider only the formulation (2.15), the case with the hybridization variable $\hat{\alpha}$ in (2.23) can be done analogously.

Computing the first variation of problem (2.15) with respect to $\bm{\sigma}$ gives

[TABLE]

for all permissible directions $\delta\bm{\sigma}$. Testing (C.1) with functions which have only support on one edge $\hat{E}$ of the triangulation $\hat{\mathcal{T}}_{h}$ yields in strong

[TABLE]

For the first variation of the membrane energy term of (2.15) in direction $v:=\delta u$ we immediately obtain for every $\hat{T}\in\hat{\mathcal{T}}_{h}$

[TABLE]

The other variations are more involved. We define the operator $\bar{(\cdot)}_{ij}:\mathbb{R}^{3\times 3}\rightarrow\mathbb{R}^{2\times 2}$, which maps $3\times 3$ matrices to its $2\times 2$ sub-matrix where the $i$-th row and $j$-th column are canceled out. Further, let $\mathcal{A}^{ij}(\cdot):\mathbb{R}^{2\times 2}\rightarrow\mathbb{R}^{3\times 3}$ denotes the operator which embeds $2\times 2$ matrices into $3\times 3$ matrices, such that $\overline{\mathcal{A}^{ij}(\bm{A})}_{ij}=\bm{A}$ and the $i$-th row and the $j$-th column of $\mathcal{A}^{ij}(\bm{A})$ are zero. Thus, $\mathcal{A}^{ij}(\cdot)$ is the right-inverse of $\bar{(\cdot)}_{ij}$.

With this, we define for $i,j\in\{1,2,3\}$

[TABLE]

Then, the following identity holds for all smooth matrix valued functions for $i,j,k\in\{1,2,3\}$

[TABLE]

With the notation $\bm{A}^{\#}_{\hat{\nu},\nu}:=\nu_{i}\bm{A}^{\#}_{ij}\hat{\nu}_{j}$ and $v:=\delta u$ there further holds

[TABLE]

Note, that (C.7) has the form of a covarient derivative.

By using $(\bm{F}\hat{\tau}_{e})\times(\bm{F}\hat{\mu})=\,\text{cof}(\bm{F})\hat{\nu}$, $\bm{F}^{\#}_{\hat{\nu},i}:\nabla_{\hat{\tau}}v$ can be simplified to

[TABLE]

Now, the volume term of (2.16) is split into two terms depending on $u$

[TABLE]

and

[TABLE]

For the boundary integral of (2.16) we use the averaged normal vector $\{\nu\}=\frac{1}{\|\nu^{o}+\nu\|_{2}}(\nu^{o}+\nu)$, with $\nu^{o}$ denoting the element normal vector on the neighbored element. This yields

[TABLE]

which can be computed exploiting (C.7). Using (2.26) instead of $\{\nu\}$ yields to a similar expression.

To obtain the boundary conditions of $u$ in strong form, which are hidden naturally in the weak form of the equation, we have to integrate by parts until no derivatives of $v$ appear.

E.g., (C.3) yields

[TABLE]

For (C.10) we have to integrate twice by parts obtaining

[TABLE]

where $\partial\partial T$ are the vertices of the element $T$ and $dss$ denotes point evaluation.

For (C.11) we get

[TABLE]

Finally, one has to use integration by parts for (C.12) to obtain the last boundary terms. Adding up all boundary terms, taking care of the constants and material parameters, one obtain the natural boundary conditions in strong form with respect to the displacement $u$.