Avoiding Membrane Locking with Regge Interpolation

TL;DR

This paper introduces a new interpolation-based method using Regge finite element space to effectively reduce membrane locking in thin shell simulations, improving computational efficiency and ease of implementation.

Contribution

It presents a novel interpolation operator into Regge finite element space that halves the membrane locking constraints compared to traditional methods.

Findings

Reduces membrane locking constraints asymptotically by half.

Easily integrable into existing shell element frameworks.

Demonstrates effectiveness through benchmark examples.

Abstract

In this paper we present a novel method to overcome membrane locking of thin shells. An interpolation operator into the so-called Regge finite element space is inserted in the membrane energy term to weaken the implicitly given kernel constraints. The number of constraints is asymptotically halved on triangular meshes compared to reduced integration techniques. Provided the interpolant, this approach can be incorporated easily into any shell element. The performance of the proposed method is demonstrated by means of several benchmark examples.

Tables

Table 1. Table 5.1: Results for cylinder with free ends × 10 5 absent superscript 10 5 \times 10^{5} , method p2 (unstructured mesh) without Regge interpolation.

	0.1	0.01	0.001	0.0001
10	-2.86758	-1.46367	-1.38530	-1.33837
38	-4.20503	-1.98840	-0.74317	-0.68349
160	-4.57442	-3.43556	-1.51590	-0.37284
644	-4.59672	-4.36907	-2.64805	-0.88707
2592	-4.59775	-4.46503	-3.95876	-2.01521
10370	-4.59781	-4.47108	-4.40993	-3.03108
41564	-4.59781	-4.47147	-4.44835	-4.20861

Table 2. Table 5.2: Results for cylinder with free ends × 10 5 absent superscript 10 5 \times 10^{5} , method p2 (unstructured mesh) with Regge interpolation.

	0.1	0.01	0.001	0.0001
10	-4.58489	-4.37366	-4.21873	-4.21326
38	-4.59430	-4.46137	-4.43641	-4.42058
160	-4.59764	-4.46845	-4.45059	-4.44846
644	-4.59780	-4.47137	-4.44950	-4.44727
2592	-4.59781	-4.47150	-4.45080	-4.44609
10370	-4.59781	-4.47150	-4.45102	-4.44607
41564	-4.59781	-4.47150	-4.45103	-4.44631

Table 3. Table 5.3: Results for cylinder with free ends × 10 5 absent superscript 10 5 \times 10^{5} , method p2 (structured mesh) without Regge interpolation.

	0.1	0.01	0.001	0.0001
8	-2.29946	-1.95947	-1.95357	-1.95351
32	-3.54279	-1.85173	-1.80002	-1.79858
128	-4.48211	-2.32588	-1.78037	-1.77192
512	-4.58991	-3.89444	-1.87506	-1.76749
2048	-4.59731	-4.42511	-2.77129	-1.78288
8192	-4.59778	-4.46852	-4.18701	-2.00808
32768	-4.59781	-4.47131	-4.43271	-3.36587
131072	-4.59781	-4.47149	-4.44987	-4.33599

Table 4. Table 5.4: Results for cylinder with free ends × 10 5 absent superscript 10 5 \times 10^{5} , method p2 (structured mesh) with Regge interpolation.

	0.1	0.01	0.001	0.0001
8	-4.58733	-4.48691	-4.48665	-4.48665
32	-4.59726	-4.46548	-4.45891	-4.45837
128	-4.59782	-4.46894	-4.45131	-4.45074
512	-4.59781	-4.47138	-4.44957	-4.44759
2048	-4.59781	-4.47150	-4.45074	-4.44619
8192	-4.59781	-4.47150	-4.45101	-4.44602
32768	-4.59781	-4.47150	-4.45103	-4.44630
131072	-4.59781	-4.47150	-4.45103	-4.44633

Table 5. Table 5.5: Results for axisymmetric hyperboloid with free ends × 10 5 absent superscript 10 5 \times 10^{5} , method p2 without Regge interpolation.

	0.1	0.01	0.001	0.0001
8	-1.40942	-0.08393	-0.00216	-0.00002
32	-1.99964	-0.41024	-0.01227	-0.00013
128	-2.22356	-1.23888	-0.11740	-0.00182
512	-2.24732	-1.77470	-0.61248	-0.02537
968	-2.24855	-1.85674	-1.03445	-0.07409
2048	-2.24893	-1.88687	-1.45795	-0.20875
8192	-2.24903	-1.89560	-1.83745	-0.91489
32768	-2.24904	-1.89616	-1.88895	-1.64433

Table 6. Table 5.6: Results for axisymmetric hyperboloid with free ends × 10 5 absent superscript 10 5 \times 10^{5} , method p2 with Regge interpolation.

	0.1	0.01	0.001	0.0001
8	-2.24057	-1.92791	-1.92433	-1.92429
32	-2.24716	-1.89856	-1.89522	-1.89517
128	-2.24903	-1.89642	-1.89304	-1.89301
512	-2.24904	-1.89622	-1.89276	-1.89271
968	-2.24904	-1.89621	-1.89271	-1.89266
2048	-2.24904	-1.89620	-1.89270	-1.89262
8192	-2.24904	-1.89620	-1.89271	-1.89259
32768	-2.24904	-1.89620	-1.89271	-1.89259

Table 7. Table 5.7: Results for uniform bending of cylindrical shell × 10 4 absent superscript 10 4 \times 10^{4} , method p2 without Regge interpolation.

	0.1	0.01	0.001	0.0001
8	5.99797	5.20307	2.43784	2.23646
16	6.00048	5.95953	4.44977	2.89527
32	6.00005	5.99788	5.88840	4.62870
128	6.00005	5.99970	5.96252	5.18276
384	6.00009	6.00007	5.99700	5.80409
1024	6.00010	6.00010	5.99987	5.98020
2560	6.00010	6.00010	6.00008	5.99852
6144	6.00010	6.00010	6.00010	5.99998
14336	6.00010	6.00010	6.00010	6.00009

Table 8. Table 5.8: Results for uniform bending of cylindrical shell × 10 4 absent superscript 10 4 \times 10^{4} , method p2 with Regge interpolation.

	0.1	0.01	0.001	0.0001
8	5.99838	6.00183	6.00186	6.00186
16	5.99925	6.00010	6.00011	6.00011
32	5.99986	6.00007	6.00007	6.00007
128	6.00004	6.00009	6.00009	6.00009
384	6.00009	6.00010	6.00010	6.00010
1024	6.00010	6.00010	6.00010	6.00010
2560	6.00010	6.00010	6.00010	6.00010
6144	6.00010	6.00010	6.00010	6.00010
14336	6.00010	6.00010	6.00010	6.00010

Table 9. Table 5.9: Results for hyperbolic paraboloid × 10 2 absent superscript 10 2 \times 10^{2} , method p2 without Regge interpolation.

	0.1	0.01	0.001	0.0001
8	1.72969	0.23110	0.03778	0.00235
32	1.80316	0.38329	0.04234	0.00063
128	1.83255	0.72759	0.12137	0.00630
512	1.83461	0.87677	0.34994	0.04454
2048	1.83476	0.89757	0.59515	0.13695
8192	1.83479	0.89897	0.67297	0.38658
32768	1.83481	0.89906	0.68057	0.55735

Table 10. Table 5.10: Results for hyperbolic paraboloid × 10 2 absent superscript 10 2 \times 10^{2} , method p2 with Regge interpolation.

	0.1	0.01	0.001	0.0001
8	1.78963	0.97481	0.90386	0.90289
32	1.83907	0.89488	0.74326	0.71977
128	1.83561	0.89766	0.68867	0.64244
512	1.83481	0.89928	0.68089	0.60693
2048	1.83477	0.89911	0.68098	0.59784
8192	1.83479	0.89907	0.68111	0.59649
32768	1.83481	0.89907	0.68110	0.59660

Table 11. Table 5.11: Results for open hemisphere with clamped ends × 10 5 absent superscript 10 5 \times 10^{5} , method p2 without Regge interpolation.

	0.1	0.01	0.001	0.0001
11	1.40550	1.41500	1.41510	1.41510
20	1.17835	1.14513	1.14461	1.14460
84	1.13244	1.12607	1.12551	1.12551
342	1.12707	1.12036	1.12034	1.12034
1368	1.12762	1.11990	1.11917	1.11919
5420	1.12765	1.12028	1.11895	1.11899
21482	1.12765	1.12031	1.11917	1.11897

Table 12. Table 5.12: Results for open hemisphere with clamped ends × 10 5 absent superscript 10 5 \times 10^{5} , method p2 with Regge interpolation.

	0.1	0.01	0.001	0.0001
11	1.31232	1.31583	1.31587	1.31587
20	1.20520	1.20922	1.20927	1.20927
84	1.12825	1.13285	1.13302	1.13303
342	1.12726	1.12282	1.12369	1.12370
1368	1.12762	1.12028	1.12045	1.12054
5420	1.12765	1.12030	1.11942	1.11964
21482	1.12765	1.12031	1.11923	1.11926

Equations

$${\text{Reg}}^k:=\{\mathbf{\sigma }\in [{\Pi }^k({\mathcal{T}}_h){]}_{\text{sym}}^{2\times 2}:[[{\mathbf{\sigma }}_{tt}]]=0\},$$

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

$${\Psi }_{E_{ij},l}:$$

$${\Psi }_{T,l}:$$

$${\mathbf{\phi }}_{E_{ij},0}:=\nabla {\lambda }_i\odot \nabla {\lambda }_j,$$

$$p_n^{(\alpha ,\beta )}(x)$$

$$\times \frac{d^n}{d{x}^n}((1-x{)}^{\alpha }(1+x{)}^{\beta }(x^2-1{)}^n),n\in {\mathbb{N}}_0,$$

$${\hat{p}}_n^{(\alpha ,\beta )}(x)$$

$${\int }_{-1}^1(1-x{)}^{\alpha }{p}_j^{\alpha }(x)p_l^{\alpha }(x)dx=\frac{2^{\alpha +1}}{2j+\alpha +1}{\delta }_{jl},$$

$${\int }_{-1}^1(1-x{)}^{\alpha }{\hat{p}}_j^{\alpha }(x){\hat{p}}_l^{\alpha }(x)dx=0\text{for }|j-l|>2.$$

$${\mathbf{\phi }}_{E_{ij},l}:={\hat{p}}_l^0\left(\frac{{\lambda }_i-{\lambda }_j}{{\lambda }_i+{\lambda }_j}\right)({\lambda }_i+{\lambda }_j{)}^l\nabla {\lambda }_i\odot \nabla {\lambda }_j.$$

$${\mathbf{\phi }}_{T^1,l_1,l_2}:=w^{(l_1,l_2)}{\lambda }_1(\nabla {\lambda }_2\odot \nabla {\lambda }_3),l_1+l_2\le k-1,$$

$${\mathbf{\phi }}_{T^2,l_1,l_2}:=w^{(l_1,l_2)}{\lambda }_2(\nabla {\lambda }_3\odot \nabla {\lambda }_1),l_1+l_2\le k-1,$$

$${\mathbf{\phi }}_{T^3,l_1,l_2}:=w^{(l_1,l_2)}{\lambda }_3(\nabla {\lambda }_1\odot \nabla {\lambda }_2),l_1+l_2\le k-1,$$

$$w^{(l_1,l_2)}:=p_{l_1}^0\left(\frac{x}{1-y}\right)(1-y{)}^{l_1}{p}_{l_2}^{2l_1+1}(2y-1).$$

$${\mathcal{I}}_{h,k}^R$$

$$\mathbf{\sigma }\mapsto \sum _{i=0}^{N_k}{\alpha }_i{\mathbf{\phi }}_i,$$

$$\sum _{E\in {\mathcal{E}}_h}{\int }_E{\mathbf{\sigma }}_h:q_{\mathcal{E}}{t}_E\otimes t_{E}ds=\sum _{E\in {\mathcal{E}}_h}{\int }_E\mathbf{\sigma }:q_{\mathcal{E}}{t}_E\otimes t_{E}ds,$$

$$\sum _{T\in {\mathcal{T}}_h}{\int }_T{\mathbf{\sigma }}_h:{\mathbf{q}}_{\mathcal{T}}dx=\sum _{T\in {\mathcal{T}}_h}{\int }_T\mathbf{\sigma }:{\mathbf{q}}_{\mathcal{T}}dx,$$

$${\mathbf{M}}^D⟨{\alpha }_i⟩=f,{\mathbf{M}}^D=\left(\begin{array}{cc}{\mathbf{M}}_{EE}^D& {\mathbf{M}}_{ET}^D\\ {\mathbf{M}}_{TE}^D& {\mathbf{M}}_{TT}^D\end{array}\right),f=\left(\begin{array}{c}{f}_E\\ f_T\end{array}\right),$$

$${\mathbf{\sigma }}_h\circ \Phi :={\mathbf{F}}^{-\top }{\hat{\mathbf{\sigma }}}_h{\mathbf{F}}^{-1}$$

$${\mathbf{\sigma }}_h\circ \Phi :=({\mathbf{F}}^{†}{)}^{\top }{\hat{\mathbf{\sigma }}}_h{\mathbf{F}}^{†}.$$

$$q_{\mathcal{E}}\circ \Phi :=J_b{\hat{q}}_{\mathcal{E}},$$

$${\mathbf{q}}_{\mathcal{T}}\circ \Phi :=\frac{1}{J}\mathbf{F}{\hat{\mathbf{q}}}_{\mathcal{T}}{\mathbf{F}}^{\top },$$

$${\int }_E{\mathbf{\sigma }}_h:q_{\mathcal{E}}t\otimes tds$$

$$={\int }_{\hat{E}}{\hat{\mathbf{\sigma }}}_h:{\hat{q}}_{\mathcal{E}}\hat{t}\otimes \hat{t}d\hat{s},$$

$${\int }_T{\mathbf{\sigma }}_h:{\mathbf{q}}_{\mathcal{T}}dx$$

$$={\int }_{\hat{T}}{\hat{\mathbf{\sigma }}}_h:{\hat{\mathbf{q}}}_{\mathcal{T}}d\hat{x}.$$

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

$$\mathcal{W}(u)=\frac{t}{2}{E}_{\text{mem}}(u)+\frac{t^3}{2}{E}_{\text{bend}}(u)+\frac{t}{2}{E}_{\text{shear}}(u)-f(u),$$

$$E_{\text{mem}}(u):={\int }_S\Vert {\mathbf{E}}_{\tau }{\Vert }_{\mathbf{M}}^{2}dx,$$

$$\Vert \cdot {\Vert }_{\mathbf{M}}^2:=\frac{\bar{E}}{1-{\bar{\nu }}^2}\left(\bar{\nu }\text{tr}(\cdot {)}^2+(1-\bar{\nu })\text{tr}({\cdot }^2)\right),$$

$$E_{\text{mem}}^{\text{lin}}(u):={\int }_S\Vert \text{sym}({\mathbf{P}}_{\tau }{\nabla }_{\tau }u){\Vert }_{\mathbf{M}}^{2}dx,$$

$$(E_{\text{mem}}^{\text{lin}}(u){)}_{\alpha \beta }={\int }_S\Vert \frac{1}{2}(u_{\alpha |\beta }+u_{\beta |\alpha })-b_{\alpha \beta }{u}_3{\Vert }_{\mathbf{M}}^{2}dx,$$

$${\int }_S\Vert {\mathcal{I}}_{h,k-1}^R{\mathbf{E}}_{\tau }{\Vert }_{\mathbf{M}}^{2}dx.$$

$$\mathcal{L}(u,\mathbf{R},\mathbf{Q}):={\int }_S\Vert \mathbf{R}{\Vert }_{\mathbf{M}}^{2}dx+⟨\mathbf{R}-{\mathbf{E}}_{\tau },\mathbf{Q}{⟩}_{{\mathcal{T}}_h},$$

$$\displaystyle\langle\bm{R}-\bm{E}_{\tau},\bm{Q}\rangle_{\mathcal{T}_{h}}:=\sum_{T\in\mathcal{T}_{h}}\Big{(}\int_{T}(\bm{R}-\bm{E}_{\tau}):\bm{Q}_{T}\,dx$$

Full text

Keywords: locking; shells; finite element method; Regge calculus

1 Introduction

In mathematical formulations of plates and shells a small parameter, the thickness $t$, is involved. The lack of finite element approximations fulfilling the implicitly given constraints of the physical model leads to so-called locking phenomena [6, 17]. As the thickness becomes small –depending heavily on the geometry and boundary conditions– the shell falls in one of two different categories: the membrane dominated or bending dominated case [12]. For shells, shear and membrane locking can be observed in the case of non-inhibited pure bending. The former, induced by the Kirchhoff constraint in the limit case, has been extensively discussed and a variety of shear locking free plate and shell elements have been proposed. For membrane locking, also called inextensional locking, where the (curved) elements fail to represent pure bending, only little numerical analysis has been done [3, 30, 24, 18, 14]. A framework based on discrete models for constructing low order membrane locking free elements has been proposed in [40]. In practice, mostly reduced integration schemes combined with stabilization techniques [51, 44, 45, 36] and assumed strain methods [32, 36, 28, 29, 13] are used. Therein, the membrane constraints are weakened due to under-integration, without deteriorating the membrane stability in the membrane dominated case, or the strain components are evaluated at certain points and extrapolated, respectively. The discrete strain gap method [30], related to assumed strain procedures, modifies the normal strains eliminating parasitic membrane parts. Further, mixed methods introducing the membrane force tensor as additional unknown have been proposed [3, 14, 23]. It is well known that p and hp-refinement strategies [38, 46, 25] may overcome the problem of membrane locking, but in the case of low order triangular elements only little influence of the reduced integration techniques has been observed [18].

Tullio Regge derived in [41] a geometric discretization of the Einstein field equations by approximations with a piece-wise constant metric. In theoretical, and later also numerical, physics so-called Regge calculus was applied e.g., in fields of relativity and quantum mechanics and has been further developed the last fifty years, see [48] for an overview. An analytical perspective of Regge calculus was given in [16, 15]. It has been observed that Regge’s approach is equivalent to specify lengths at all edges of a mesh, analogical to Whitney-forms [47]. In the context of finite element exterior calculus (FEEC) [5, 4] a finite element point of view was given in [19, 20] and the resulting Regge elements have been generalized to arbitrary polynomial order on triangles and tetrahedrons [31].

In this work the resulting Regge interpolant is used to construct membrane locking free methods for shells, staying stable in the case of inhibited pure bending. It can be incorporated into any existing method and finite element code, provided the interpolation operator. In a variety of numerical examples the performance of the method is tested.

This paper is structured as follows: In the next section we will give an overview of the construction of Regge elements and the corresponding interpolation operator. Section 3 is devoted to the description of the proposed method in the general setting of shells. In Section 4 the method is discussed and in the last section we apply the resulting algorithm to several established membrane locking benchmark examples.

2 Regge elements

For the reader’s convenience we give a brief introduction in the construction of Regge finite elements. Then the implementation of the Regge interpolation operator is discussed.

As the Regge elements approximate symmetric tensor fields, we seek for a matrix valued finite element space. To prescribe the edge lengths globally only the tangential-tangential components need to be continuous. Therefore, let $\mathcal{T}_{h}$ a triangulation of a bounded domain $\Omega\subset\mathbb{R}^{2}$ and $\mathcal{E}_{h}$ the corresponding skeleton, i.e., the set of all edges. The set of all piece-wise polynomials up to degree $k$ on $\mathcal{T}_{h}$ and $\mathcal{E}_{h}$ is denoted by $\Pi^{k}(\mathcal{T}_{h})$ and $\Pi^{k}(\mathcal{E}_{h})$, respectively. For each element $T\in\mathcal{T}_{h}$ the tangential and outer normal vector on the boundary $\partial T$ are given by $t$ and $n$, respectively. The outer dyadic product of two vectors is denoted by $\otimes$.

The Regge finite element space is given by

[TABLE]

where $\bm{\sigma}_{tt}:=(\bm{I}-n\otimes n)\bm{\sigma}(\bm{I}-n\otimes n)$ is the tangential-tangential component of $\bm{\sigma}$, $\bm{I}$ the identity matrix, and $\llbracket\bm{\sigma}_{tt}\rrbracket$ denotes the corresponding jump over elements.

Further, we define the Lagrangian nodal finite element space as

[TABLE]

where $C^{0}(\Omega)$ denotes to set of all continuous functions on $\Omega$. For a (high-order) construction of $H^{1}$-conforming finite elements we refer to [10, 50, 49].

In the context of [27] we define functionals, the degrees of freedom (dofs), on the reference or physical triangle $T$ with the local space $[\Pi^{k}(T)]^{2\times 2}_{\text{sym}}$. Therefore, let $\{q_{E,l}\}$ and $\{\bm{q}_{T,l}\}$ denote a polynomial basis of $\Pi^{k}(E_{ij})$ on the edge $E_{ij}$ between two vertices, $i\neq j$, and $[\Pi^{k-1}(T)]^{2\times 2}_{\text{sym}}$ on the triangle $T$, respectively. Then the functionals read

[TABLE]

where $\bm{A}:\bm{B}:=\sum_{ij}\bm{A}_{ij}\bm{B}_{ij}$ denotes the Frobenius scalar product and $t_{E}$ the tangent vector of the edge $E_{ij}$. Note that $t_{E}\otimes t_{E}$ is single valued, i.e., does not depend on the orientation of $t_{E}$.

In one dimension the Regge elements coincide with $L^{2}$-conforming discontinuous finite elements. On triangular elements the dofs are associated with the edges (2.3), analogical to $H(\text{curl})$-conforming [33, 34] elements, and inner bubbles (2.4) for higher polynomial degrees. The lowest order (polynomial order $k=0$) and first order Regge elements on segments and triangles are illustrated in Figure 2.1.

We now give an explicit basis for the corresponding shape functions $\{\bm{\varphi}_{i}\}$ to (2.3)–(2.4). Therefore, let $\lambda_{i}\in\Pi^{1}(\hat{T})$ denote the barycentric coordinates of the vertices $V_{1}=(-1,0)$, $V_{2}=(1,0)$, and $V_{3}=(0,1)$ of the reference triangle $\hat{T}$, i.e. $\lambda_{i}(V_{j})=\delta_{ij}$, where $\delta_{ij}$ denotes the Kronecker delta. Then, the lowest order basis functions are given by

[TABLE]

where $\odot$ denotes the symmetric dyadic product of two vectors. For a generalization to arbitrary order $k$ let for $\alpha,\beta>-1$, $p_{n}^{(\alpha,\beta)}$ and $\hat{p}_{n}^{(\alpha,\beta)}$ denote the $n$-th Jacobi and $n$-th integrated Jacobi polynomial [1, 2, 8], respectively,

[TABLE]

As we are only interested in the case $\beta=0$ the notation can be simplified by $p_{n}^{(\alpha,0)}(x)=p_{n}^{\alpha}(x)$ and $\hat{p}_{n}^{(\alpha,0)}(x)=\hat{p}_{n}^{\alpha}(x)$. The latter fulfill the following orthogonality properties

[TABLE]

Note that with $\alpha=0$ the (integrated) Jacobi polynomials reduce to the (integrated) Legendre polynomials. Hence, the high order edge basis functions read, $l=1,\dots,k$,

[TABLE]

The (high order) cell basis functions of order $k>0$ are given by, $l_{1},l_{2}\geq 0$,

[TABLE]

with the Dubiner basis

[TABLE]

The shape functions (2.5), (2.10), and (2.11) build a basis of $[\Pi^{k}(\mathcal{T}_{h})]^{2\times 2}_{\text{sym}}$, which has a dimension of $3(k+1)(k+2)/2$: The edge shape functions are linearly independent as $(\bm{\varphi}_{E_{i},l})_{t_{E_{j}}t_{E_{j}}}=\delta_{ij}$, where $E_{i}$ and $t_{E_{i}}$ denotes the i-th edge and corresponding tangent vector. Further, $(\bm{\varphi}_{T^{i},l_{1},l_{2}})_{t_{E_{j}}t_{E_{j}}}=0$ for $i,j=1,2,3$ and thus, the inner shapes are independent of the edge basis. The claim follows together with the independence of $\lambda_{i}(\nabla\lambda_{j}\odot\nabla\lambda_{k})$ for $i\neq j\neq k$ and counting all shape functions. A different basis for triangular Regge elements is constructed in [31].

Given functionals (2.3)–(2.4) $\{\Psi_{i}\}$ and the corresponding shape functions $\{\bm{\varphi}_{i}\}$, one can define the following Regge interpolation operator

[TABLE]

where $N_{k}\in\mathbb{N}$ denotes the number of degrees of freedom (the number of shape functions), $C^{\infty}(\Omega)$ the set of all smooth functions on $\Omega$, and the coefficients $\alpha_{i}$ are obtained by the following consideration:

Let $\bm{\sigma}_{h}=\mathcal{I}^{R}_{h,k}(\bm{\sigma})$ be the interpolation of a given function $\bm{\sigma}$. Then, with (2.3)–(2.4), $\bm{\sigma}_{h}$ is the solution of the following system of equations

[TABLE]

for all $\bm{q}_{\mathcal{T}}\in[\Pi^{k-1}(\mathcal{T}_{h})]^{2\times 2}_{\text{sym}}$ and $q_{\mathcal{E}}\in\Pi^{k}(\mathcal{E}_{h})$. In matrix form we obtain the linear equation

[TABLE]

with $\langle\alpha_{i}\rangle$ denoting the corresponding coefficient vector of $\bm{\sigma}_{h}$ and $\bm{M}^{D}_{EE}=\Psi_{E_{ij}}(\bm{\varphi}_{E_{ij}})$, $\bm{M}^{D}_{ET}=\Psi_{E_{ij}}(\bm{\varphi}_{T})$, $f_{E}=\Psi_{E_{ij}}(\bm{\sigma})$ and analogously for the other components. There holds by construction $\Psi_{E_{ij}}(\bm{\varphi}_{T})=0$ and thus, $\bm{M}^{D}_{ET}=0$, resulting in a lower block triangular dual mass matrix, see Figure 2.2 (a). This structure can be exploited for inverting $\bm{M}^{D}$. Note that the matrix is not symmetric, as in general $\Psi_{T}(\bm{\varphi}_{E_{ij}})\neq 0$. As will be discussed later the tangential-tangential continuity of the elements can be broken leading to a block diagonal matrix depicted in Figure 2.2 (b). Note that in this case the dofs are ordered element-wise instead of splitting them into edge and inner dofs.

Given a Regge element $\hat{\bm{\sigma}}_{h}$ on the reference triangle $\hat{T}$, one needs to map it to a physical element $T=\Phi(\hat{T})$, $\Phi:\hat{T}\rightarrow\mathbb{R}^{2}$, in such a way that the tangential-tangential continuity is preserved. For $H(\text{curl})$-conforming elements the covariant transformation $u\circ\Phi:=\bm{F}^{-\top}\hat{u}$, $\bm{F}:=\nabla\Phi$, is used to preserve the tangential continuity. Thus, by applying the covariant transformation on both sides

[TABLE]

leads to a symmetric tangential-tangential continuous function $\bm{\sigma}_{h}$ on $T$. For shells one needs to map the reference triangle $\hat{T}$ to a (possibly curved) surface triangle $T\subset\mathbb{R}^{3}$. The deformation gradient $\bm{F}\in C^{\infty}(\Omega,\mathbb{R}^{3\times 2})$ is not invertible and thus, the transformation rule (2.16) needs to be adopted by using the Moore–Penrose pseudo-inverse instead

[TABLE]

It is possible to obtain a geometry free dual mass matrix $\bm{M}^{D}$ by transforming the basis $q_{\mathcal{E}}$ and $\bm{q}_{\mathcal{T}}$ in (2.14) with

[TABLE]

where $J_{b}:=\|\bm{F}t\|_{2}$ denotes the boundary determinant and $J:=\det(\bm{F})$, leading to ($T=\Phi(\hat{T})$, $E=\Phi(\hat{E})$)

[TABLE]

Therefore, one may exploit this property to assemble only (permutations of) one element and “fill” the whole matrix $\bm{M}^{D}$. This procedure works also for transformation (2.17), where one has to choose the surface determinant $J=\sqrt{\det(\bm{F}^{\top}\bm{F})}$ for $\bm{F}\in\mathbb{R}^{3\times 2}$ or $J=\|\text{cof}(\bm{F})\nu\|_{2}$ for the extended matrix $\bm{F}\in\mathbb{R}^{3\times 3}$ in (2.19) with $\text{cof}(\bm{F})$ denoting the cofactor matrix and $\nu$ the surface normal vector.

3 Methodology

3.1 Shells and membrane energy

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

[TABLE]

Furthermore, let $\Phi:\Omega\rightarrow\mathbb{R}^{3}$ be the deformation from the initial to the deformed configuration of the shell and $\phi:S\rightarrow\mathbb{R}^{3}$ the deformation of the approximated mid-surface. The corresponding triangulation consisting of possibly curved triangles of $S$ is denoted by $\mathcal{T}_{h}$. Then, we define $\bm{F}:=\nabla_{\tau}\phi$ as the deformation gradient. Here, $\nabla_{\tau}\phi$ denotes the surface gradient of $\phi$, which can be introduced in weak sense [22], or directly as Fréchet-derivative. We can split the deformation into the identity function and the displacement, $\phi=\text{id}+u$, and thus, $\bm{F}=\bm{P}_{\tau}+\nabla_{\tau}u$ with the projection onto the tangent plane $\bm{P}_{\tau}:=\bm{I}-\nu\otimes\nu$.

The shell energy functional can be split into a membrane, bending and shear energy part, cf. [21, 13, 9],

[TABLE]

where $f$ denotes the external forces.

We focus on the membrane energy and consider the full nonlinear term

[TABLE]

with $\bm{E}_{\tau}:=1/2(\bm{C}_{\tau}-\bm{I})$ denoting the Green strain tensor restricted on the tangent plane, $\bm{C}_{\tau}:=\bm{F}^{\top}\bm{F}$ the Cauchy–Green strain tensor. The material norm is given by

[TABLE]

with the material tensor $\bm{M}$, the Young’s modulus $\bar{E}$, and the Poisson’s ratio $\bar{\nu}$, respectively.

The linearization of (3.3) is given by

[TABLE]

where $u_{\alpha|\beta}$ denotes the covariant derivative, $b_{\alpha\beta}$ the second fundamental form, $\alpha,\beta\in\{1,2\}$, and $u_{3}$ the displacement component in normal direction, see e.g. [13] for the notation of curvilinear coordinates.

3.2 Usage of Regge interpolant

In what follows let the discrete displacements $u_{h}\in[V_{h}^{k}]^{3}$. For the proposed method we insert the Regge interpolation operator of order $k-1$ into the membrane energy (3.3)

[TABLE]

Due to the tangential-tangential continuity of $\bm{E}_{\tau}$ –the discrete Jacobian $\nabla_{\tau}u_{h}$ is tangential-continuous– it is mathematically equivalent to apply the projection operator $\mathcal{I}^{R}_{h,k-1}$ only element-wise. This makes the method cheap, as no additional global system has to be solved. Further, properties as symmetry and positivity gets preserved. Therefore, only small problems of the form (2.14) on one element have to be solved on each integration point.

Further, (3.6) is equivalent to a three-field formulation by introducing the (locally and thus discontinuous) Regge interpolant $\bm{R}\in\text{Reg}^{k-1,\mathrm{dc}}$ and the corresponding local shape functionals $\bm{Q}\in[\text{Reg}^{k-1,\mathrm{dc}}]^{*}$ as additional unknowns. Note that the functionals (2.3)–(2.4) span the topological dual space $[\text{Reg}^{k-1,\mathrm{dc}}]^{*}$. The corresponding Lagrangian reads

[TABLE]

where, according to (2.3)–(2.4),

[TABLE]

and thus, with (2.14), $\bm{R}=\mathcal{I}^{R}_{h,k-1}\bm{E}_{\tau}$. In numerical experiments we observed that if the full nonlinear Green strain tensor $\bm{E}_{\tau}$ is used in (3.7) less Newton iterations are needed than for the direct interpolation procedure (3.6).

4 Discussion

Let $u$ be the exact solution of the shell problem (3.2) in the case of non-inhibited pure bending such that $E^{\text{lin}}_{\text{mem}}(u)=0$. Interpolating $u$ into the Lagrangian finite element space $[V_{h}^{k}]^{3}$, $u_{h}:=\mathcal{I}_{h,k}u$, $\mathcal{I}_{h}$ denoting the standard nodal interpolation operator, does not guarantee in general that $E^{\text{lin}}_{\text{mem}}(u_{h})=0$ for the discrete displacements. I.e., the interpolation operator does not preserve the kernel of the membrane operator. Therefore, pure bending modes induce discrete membrane energy modes due to the discrete constraints. This effect dominates for small thickness parameters $t$, the shell element is called to be too stiff and locking occurs.

By using the Regge interpolant $\mathcal{I}^{R}_{h,k-1}E^{\text{lin}}_{\text{mem}}(u_{h})$ we weak the discrete constraints. Reduced integration schemes follow the same idea, using less Gauß-integration points, which corresponds to an $L^{2}$ instead of a Regge interpolation. When we compare the number of dofs, which can be interpreted as the number of constraints, one can observe that on a single triangle $T$ the number of constraints are equal, as the dimension of both spaces are the same, $\text{dim}=3(k+1)(k+2)/2$.

For a triangulation $\mathcal{T}_{h}$, however, the number of constraints differ already in the lowest order case significantly. For Regge elements we have one degree of freedom per edge, whereas in the reduced integration scheme one has three per element. Asymptotically there holds

[TABLE]

where $\#T$, $\#E$, and $\#V$ denote the number of triangles, edges, and vertices of the triangulation $\mathcal{T}_{h}$, respectively. Therefore,

[TABLE]

and thus, the Regge interpolation reduces the number of constraints asymptotically by a factor of two compared to the $L^{2}$-projection. Furthermore, on a triangulation $\mathcal{T}_{h}$ of a flat two-dimensional domain or a surface described by one single embedding (and thus not closed) there holds

[TABLE]

where $\#V_{B}$ and $\#V_{I}$ denote the number of vertices on the boundary and in the inner domain of the surface, respectively. The discrepancy of three corresponds to the number of rigid-body motions in two dimensions, two translations and one rotation. Therefore, for given displacements at the vertices one can find a unique value per edge describing the (tangential-tangential) stretching between two vertices. This fits perfectly to the following (linear) exact sequence

[TABLE]

where $RB:=\{Ax+b\,|\,A\in\mathbb{R}^{2\times 2},\,A^{\top}=-A,\,b\in\mathbb{R}^{2}\}$ denotes the set of linearized rigid body motions.

In [20, 26] they used this sequence in three dimensions as a part of a larger complex and proofed in the lowest order case commuting and exactness properties. For a nonlinear complex one has to replace the symmetric gradient by the Green strain tensor and $RB=\{Ax+b\,|\,A\in SO(2),\,b\in\mathbb{R}^{2}\}$, where $SO(2)$ denotes the set of all orthogonal $2\times 2$ matrices with determinant one.

In case of the full nonlinear membrane energy term (3.3) the Green strain operator $\bm{E}_{\tau}:[\Pi^{k}(\mathcal{T}_{h})]^{d}\rightarrow[\Pi^{2k-2}(\mathcal{T}_{h})]^{d\times d}_{\text{sym}}$ doubles the polynomial degree asymptotically element-wise, with the exception $k=1$. This may lead to even worse discrete kernel conservation. Thanks to the Regge interpolant, the Green strain tensor gets projected back to polynomial degree $k-1$ and again the number of constraints are reduced.

The idea of inserting an interpolation operator has already been successfully applied to avoid shear locking. E.g. for the mixed interpolated tensorial components (MITC) elements [7, 11] an $H(\text{curl})$ interpolant is inserted into the shear energy term. Also methods where the rotations get directly approximated by $H(\text{curl})$-conforming finite elements overcome shear locking [37].

In the lowest order case $k=1$ for the displacements, membrane locking is not observed as long as an isoperimetric mapping for the shell geometry is considered. Curving the geometry by a higher polynomial degree as the displacements leads to enormous membrane locking in the lowest order case. However, using the Regge interpolation $\mathcal{I}^{R}_{h,0}$ reduces this locking phenomena too.

5 Numerical examples

To avoid shear locking effects we use the Kirchhoff–Love shell model introduced in [35], where the Regge interpolation has been successfully used for triangular meshes. The method is implemented in the open source finite element library Netgen/NGSolve222www.ngsolve.org [42, 43].

For the benchmarks we use second order finite elements for the displacements, where the geometry is mapped isoperimetrically, i.e. curved elements are used, called method p2.

The forces are chosen such that the deformations are in the linear regime. Therefore, the differences between the linearized (3.1) and full nonlinear (3.3) membrane energy is marginal. Further, the forces are scaled appropriately with the thickness parameter $t$ ($t^{3}$ in the bending dominated and $t$ in the membrane dominated case) such that the deformations are in the same magnitude. Due to the nonlinear membrane and bending energy, however, the results may vary little with respect to the thickness parameter. The reference values are computed by using fourth order finite elements for the displacement on the finest mesh, called method p4, and the relative error is computed by |result - reference|/|reference|.

5.1 Cylinder with free ends

A cylinder with free ends is loaded with a periodic force [39, 13], see Figure 5.1. By symmetries the computational domain is one eighth of the original and symmetry boundary conditions are used, see Figure 5.2 and 5.3. The material and geometric parameters are $R=1$, $E=3\times 10^{4}$, $\nu=0.3$, $t\in\{0.1,0.01,0.001,0.0001\}$ and the cylinder is loaded by the normal pressure distribution $P=t^{3}\cos(2\zeta)\hat{\nu}$, $\zeta$ and $\hat{\nu}$ denoting the circumferential arc-length and the normal vector on the reference configuration, respectively, cf. Figure 5.1.

The radial deflection at point $A$ is measured and listed in Tables 5.1-5.4. The relative error for unstructured meshes can be found in Figure 5.4 and for structured meshes in Figure 5.5. There the classical locking behavior can be observed if the Regge interpolant is not used, as the pre-asymptotic range increases rapidly for smaller thicknesses. Using Regge interpolation avoids this pre-asymptotic behavior. Further, for a small amount of elements the relative errors start already with $<10$ percent, also for thick parameters $t$. The results on the unstructured meshes are comparable to the structured one.

5.2 Axisymmetric hyperboloid with free ends

An axisymmetric hyperboloid is described by the equation

[TABLE]

with free boundaries is loaded by a force, see [13]. Due to symmetries it is sufficient to use one eighth of the geometry and symmetry boundary conditions, see Figure 5.6 for the geometry and Figure 5.7 for a coarse and fine mesh. The material and geometric parameters are $R=1$, $E=2.85\times 10^{4}$, $\nu=0.3$, $t\in\{0.1,0.01,0.001,0.0001\}$, $P=\frac{t^{3}}{\sqrt{x^{2}+y^{2}}}\cos(2\zeta)\begin{pmatrix}x\\ y\\ 0\end{pmatrix}$ similar to the previous benchmark.

The radial deflection at point $A$ is listed in Table 5.5-5.6 and the relative error in Figure 5.8. Again the results improve using Regge interpolation and we emphasize that for $t=0.0001$ with 8 elements the difference with a factor of $10^{5}$ is immensely ($-2\times 10^{-10}$ vs $-2\times 10^{-5}$ with the reference value $-1.89\times 10^{-5}$).

5.3 Uniform bending of cylindrical shell

A moment $M$ is applied to a cylindrical shell, which is fixed at the top [30]. The material and geometric parameters are $R=0.1$, $b=0.025$, $E=2\times 10^{5}$, $\nu=0$, $t\in\{0.1,0.01,0.001,0.0001\}$, $M_{0}=(t/R)^{3}$, see Figure 5.9 and 5.10.

This time the deflection orthogonal to the radial direction is computed at point $A$. The results can be found in Table 5.7-5.8 and Figure 5.11. In this benchmark the method without interpolation operator does not produce a strong pre-asymptotic regime for small thicknesses. However, the initial relative error gets larger. In contrast, the errors with the Regge interpolation start all at nearly the same value and show a uniform convergence behavior.

5.4 Hyperbolic paraboloid

A hyperbolic paraboloid, which is described by the embedding

[TABLE]

is clamped at the bottom and subjected to a surface force $f$ [18]. On the right side symmetry boundary conditions are used, the other boundaries are free. The material and geometric parameters are $\alpha=0.2$, $E=2.85\times 10^{4}$, $\nu=0.3$, $t\in\{0.1,0.01,0.001,0.0001\}$, $f=8t^{3}\hat{\nu}$. Here, $\hat{\nu}$ denotes the normal vector on the reference configuration, see Figure 5.12 and 5.13.

The deflection in z-direction at point $A$ can be seen in Table 5.9-5.10 and the relative error is depicted in Figure 5.14.

5.5 Open hemisphere with clamped ends

As a membrane dominated example an 18∘ open hemisphere with clamped top and bottom edges is used [13]. Due to symmetry only one fourth of the hemisphere is considered with appropriate symmetry boundary conditions, see Figure 5.15 and 5.16. The material and geometric parameters are $R=10$, $E=6.825\times 10^{7}$, $\nu=0.3$, $t\in\{0.1,0.01,0.001,0.0001\}$ the volume force density is $P=\frac{t}{10}\cos(2\zeta)\hat{\nu}$, where $\zeta$ denotes the angle between the $x$ and $y$ component.

The deflection in $x$-direction at point $A$ is listed in Table 5.11-5.12 and Figure 5.17 shows the relative error. As expected the method without Regge interpolation does not lock in the case of inhibited pure bending. Using the interpolation operator yields to only slight deterioration in the convergence rates for smaller thicknesses. However, this effect is marginal compared to the improvements in the bending dominated regime benchmarks. Further, also in this membrane dominated example the interpolation yields to better results for a small amount of elements.

6 Conclusions

In this work the Regge interpolation operator was inserted into the membrane energy part relaxing the kernel constraints and avoiding membrane locking for thin shells. For triangular meshes the number of constraints is significantly reduced without deteriorating the membrane stability in the membrane dominated regime. The performance was demonstrated by benchmark examples including membrane and bending dominated cases. A rigorous mathematical proof of uniform convergence independently of the thickness parameter is topic of further research. In contrast to shear locking, which can be already observed for simple plate problems, membrane locking occurs only if curved elements are used and is thus more involved.

Acknowledgements

The authors acknowledge support from the Austrian Science Fund (FWF) through grant number W 1245.