DKMQ24 shell element with improved membrane behaviour

TL;DR

This paper introduces DKMQ24$_2$+,

Contribution

It presents an improved shell element with enhanced membrane behavior using new approximation and integration techniques.

Findings

Free of shear and membrane locking.

Lower numerical errors in benchmark tests.

Effective for membrane-dominated problems.

Abstract

An approach to improve the membrane behaviour of the four-node shell element with 24 degrees of freedom DKMQ24 proposed by Katili et al. (2015) is presented. This improvement is based on a different approximation of drilling rotations, based on Allman's shape functions. Further, the element formulation is enhanced by the use of selective reduced integration of shear terms and by the proportional scaling of the penalty constant to the shell thickness. Additionaly, nodal moment corrections for the approximation of distributed normal loads are proposed. This element, called DKMQ24$_2$+, was tested on nine standard benchmark problems for the case of linear elasticity. The element is free of shear and membrane locking and shows lower numerical errors for membrane and membrane-dominated problems.

Figures

Tables

Table 1. Table 1: Basis functions

node $i$	$[r_i,s_i]$	$a_i(r,s)$
$1$	$[-1,-1]$	$\frac{1}{4}\left(1-r\right)\left(1-s\right)$
$2$	$[1,-1]$	$\frac{1}{4}\left(1+r\right)\left(1-s\right)$
$3$	$[1,1]$	$\frac{1}{4}\left(1+r\right)\left(1+s\right)$
$4$	$[-1,1]$	$\frac{1}{4}\left(1-r\right)\left(1+s\right)$

Table 2. Table 2: The supplementary basis functions for midside rotations △ ​ β t ​ k △ subscript 𝛽 𝑡 𝑘 \triangle\beta_{tk} on side k = 5 , … , 8 𝑘 5 … 8 k=5,\dots,8

node $k$	$[r_k,s_k]$	$a_k(r,s)$
$5$	$[0,-1]$	$\frac{1}{2}\left(1-r^2\right)\left(1-s\right)$
$6$	$[1,0]$	$\frac{1}{2}\left(1+r\right)\left(1-s^2\right)$
$7$	$[0,1]$	$\frac{1}{2}\left(1-r^2\right)\left(1+s\right)$
$8$	$[-1,0]$	$\frac{1}{2}\left(1-r\right)\left(1-s^2\right)$

Table 3. Table 3: Quadrature rules used for the stiffness matrix assembling

Term	Gauss integration rule
${\int }_{\text{e}}{𝑩}^{\text{T}}𝑫𝑩\text{d}A$	$2\times 2$
$\frac{5}{6}G{\int }_{\text{e}}h{𝒃}_{\text{s}x}{𝒃}_{\text{s}x}^{\text{T}}\text{d}A$	$1\times 2$
$\frac{5}{6}G{\int }_{\text{e}}h{𝒃}_{\text{s}y}{𝒃}_{\text{s}y}^{\text{T}}\text{d}A$	$2\times 1$
${𝑲}_{\text{stab}}$	$1\times 1$

Table 4. Table 4: Tested elements

Notation	Source	Description
MITC4+	from ref.[32], (2017)	-
​US-ATFHS8	from ref.[15], (2018)	8-nodal solid-shell element, without rotational dofs, 24 dofs
DKMQ24 ​[18]	from ref.[18], (2015)	-
DKMQ24	our implementation	DKMQ24 element ($c=0.001$ as in [18])
DKMQ24($c={10}^{-7}$)	our implementation	DKMQ24 element ($c={10}^{-7}$)
DKMQ241+	our implementation	improved DKMQ24 element, which uses Eq. (409) ($c_1=0.15$)
DKMQ242+	our implementation	improved DKMQ24 element, which uses Eq. (427) ($c_2=0.1$)

Table 5. Table 5: Deflection u Z ​ ( A ) subscript 𝑢 𝑍 𝐴 u_{Z}(A) for Cook’s membrane problem

Mesh	MITC4+		DKMQ24		DKMQ241+		DKMQ242+
	$u_Z\left(A\right)$	err.	$u_Z\left(A\right)$	err.	$u_Z\left(A\right)$	err.	$u_Z\left(A\right)$	err.
	[mm]	[%]	[mm]	[%]	[mm]	[%]	[mm]	[%]
$2\times 2$	$11.729$	$51.06$	$11.845$	$50.57$	$22.353$	$6.72$	$22.396$	$6.55$
$4\times 4$	$18.266$	$23.78$	$18.299$	$23.64$	$23.397$	$2.37$	$23.415$	$2.30$
$8\times 8$	$22.075$	$7.89$	$22.079$	$7.87$	$23.764$	$0.84$	$23.776$	$0.79$
$16\times 16$	$23.430$	$2.23$	$23.430$	$2.23$	$23.892$	$0.30$	$23.897$	$0.29$
$32\times 32$	$23.818$	$0.62$	$23.818$	$0.61$	$23.939$	$0.11$	$23.940$	$0.11$
Reference solution: $u_Z\left(A\right)=23.965\text{mm}$ (DKMQ242+ element on mesh $192\times 128$)

Table 6. Table 6: Deflection u X ​ ( A ) subscript 𝑢 𝑋 𝐴 u_{X}(A) for the hemispherical shell problem

Mesh	MITC4+		DKMQ24 [18]		DKMQ24		DKMQ24		DKMQ241+		DKMQ242+
							($c={10}^{-7}$)
	$u_X\left(A\right)$	err.	$u_X\left(A\right)$	err.	$u_X\left(A\right)$	err.	$u_X\left(A\right)$	err.	$u_X\left(A\right)$	err.	$u_X\left(A\right)$	err.
	[mm]	[%]	[mm]	[%]	[mm]	[%]	[mm]	[%]	[mm]	[%]	[mm]	[%]
$4\times 4$	$94.802$	$2.60$	$78.721$	$14.80$	$49.911$	$45.98$	$85.305$	$7.68$	$40.625$	$56.03$	$48.902$	$47.08$
$8\times 8$	$92.585$	$0.20$	$93.624$	$1.32$	$74.940$	$18.90$	$93.963$	$1.69$	$88.929$	$3.76$	$89.224$	$3.44$
$16\times 16$	$92.030$	$0.40$	$92.246$	$0.17$	$83.775$	$9.33$	$92.289$	$0.12$	$92.356$	$0.05$	$92.372$	$0.03$
$32\times 32$	$92.234$	$0.18$	-	-	$86.654$	$6.22$	$92.146$	$0.28$	$92.405$	$0.01$	$92.406$	$0.01$
Reference solution: $u_X\left(A\right)=-u_Y\left(B\right)=92.400\text{mm}$ ([29])

Table 7. Table 7: Deflection u Z ​ ( B ) subscript 𝑢 𝑍 𝐵 u_{Z}(B) for the Scordelis-Lo roof problem, B ≈ [ 1.928 , 3 , 2.298 ] 𝐵 1.928 3 2.298 B\approx[1.928,3,2.298]

Mesh	DKMQ24 [18]		DKMQ24		DKMQ241+		DKMQ242+
	$u_Z\left(B\right)$	err.	$u_Z\left(B\right)$	err.	$u_Z\left(B\right)$	err.	$u_Z\left(B\right)$	err.
	[m]	[%]	[m]	[%]	[m]	[%]	[m]	[%]
$4\times 4$	$-0.03425$	$5.12$	$-0.034258$	$5.10$	$-0.037036$	$2.59$	$-0.037045$	$2.62$
$8\times 8$	$-0.03528$	$2.26$	$-0.035284$	$2.26$	$-0.036228$	$0.35$	$-0.036233$	$0.37$
$16\times 16$	$-0.03585$	$0.70$	$-0.035858$	$0.67$	$-0.036102$	$0.01$	$-0.036105$	$0.01$
$32\times 32$	-	-	$-0.036070$	$0.08$	$-0.036086$	$0.04$	$-0.036078$	$0.04$
Reference solution: $u_Z\left(B\right)=-0.0361\text{m}$ (theory of deep shell, [34])

Table 8. Table 8: Deflection u Z ​ ( C ) subscript 𝑢 𝑍 𝐶 u_{Z}(C) for the Scordelis-Lo roof problem, C = [ 0 , 3 , 3 ] 𝐶 0 3 3 C=[0,3,3]

Mesh	DKMQ24 [18]		DKMQ24		DKMQ241+		DKMQ242+
	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.
	[m]	[%]	[m]	[%]	[m]	[%]	[m]	[%]
$4\times 4$	$0.00513$	$5.18$	$0.005130$	$5.17$	$0.005599$	$3.49$	$0.005600$	$3.51$
$8\times 8$	$0.00529$	$2.21$	$0.005294$	$2.15$	$0.005447$	$0.67$	$0.005447$	$0.69$
$16\times 16$	$0.00538$	$0.55$	$0.005378$	$0.58$	$0.005419$	$0.17$	$0.005419$	$0.17$
$32\times 32$	-	-	$0.005407$	$0.05$	$0.005414$	$0.07$	$0.005415$	$0.08$
Reference solution: $u_Z\left(C\right)=0.00541\text{m}$ (theory of deep shell, [34])

Table 9. Table 9: Internal moment m x ​ ( C ) subscript 𝑚 𝑥 𝐶 m_{x}(C) for the Scordelis-Lo roof problem

Mesh	DKMQ24 [18]		DKMQ24		DKMQ241+		DKMQ242+
	$m_x\left(C\right)$	err.	$m_x\left(C\right)$	err.	$m_x\left(C\right)$	err.	$m_x\left(C\right)$	err.
	[N]	[%]	[N]	[%]	[N]	[%]	[N]	[%]
$4\times 4$	$1874$	$8.85$	$1740.1$	$15.36$	$1930.7$	$6.09$	$1931.3$	$6.07$
$8\times 8$	$1998$	$2.82$	$1947.8$	$5.26$	$2023.4$	$1.59$	$2023.7$	$1.57$
$16\times 16$	$2041$	$0.73$	$2027.3$	$1.40$	$2048.0$	$0.39$	$2048.1$	$0.38$
$32\times 32$	$2054$	$0.10$	$2050.4$	$0.27$	$2054.5$	$0.07$	$2054.6$	$0.07$
Reference solution: $m_x\left(C\right)=2056\text{N}$ (theory of shallow shell, [34])

Table 10. Table 10: Internal force n y ​ ( B ) subscript 𝑛 𝑦 𝐵 n_{y}(B) for the Scordelis-Lo roof problem

Mesh	DKMQ24 [18]		DKMQ24		DKMQ241+		DKMQ242+
	$n_y\left(B\right)$	err.	$n_y\left(B\right)$	err.	$n_y\left(B\right)$	err.	$n_y\left(B\right)$	err.
	[Nm-1]	[%]	[Nm-1]	[%]	[Nm-1]	[%]	[Nm-1]	[%]
$4\times 4$	$498430$	$22.24$	$498741$	$22.19$	$549772$	$14.23$	$550517$	$14.12$
$8\times 8$	$592680$	$7.54$	$592745$	$7.53$	$607481$	$5.23$	$607709$	$5.19$
$16\times 16$	$620750$	$3.16$	$620909$	$3.13$	$624429$	$2.59$	$624490$	$2.58$
$32\times 32$	$628280$	$1.98$	$628908$	$1.89$	$629006$	$1.87$	$629044$	$1.87$
Reference solution: $n_y\left(B\right)=641000{\text{Nm}}^{-1}$ (theory of shallow shell, [34])

Table 11. Table 11: Deflection u Z ​ ( B ) subscript 𝑢 𝑍 𝐵 u_{Z}(B) for the Scordelis-Lo roof problem with modified parameters according to [ 29 ] ( L = 50 ​ m , R = 25 ​ m , h = 0.25 ​ m , E = 432 ​ MPa , ν = 0 , ρ = 36 ​ kg m − 3 , g = 10 ​ m s − 2 formulae-sequence 𝐿 50 m formulae-sequence 𝑅 25 m formulae-sequence ℎ 0.25 m formulae-sequence 𝐸 432 MPa formulae-sequence 𝜈 0 formulae-sequence 𝜌 36 superscript kg m 3 𝑔 10 superscript m s 2 L=50\,\text{m},\,R=25\,\text{m},\,h=0.25\,\text{m},\,E=432\,\text{MPa},\nu=0,\,\rho=36\,\text{kg\,m}^{-3},\,g=10\,\text{m\,s}^{-2} )

Mesh	MITC4+		US-ATFHS8		DKMQ24		DKMQ241+		DKMQ242+
	$u_Z\left(B\right)$	err.	$u_Z\left(B\right)$	err.	$u_Z\left(B\right)$	err.	$u_Z\left(B\right)$	err.	$u_Z\left(B\right)$	err.
	[m]	[%]	[m]	[%]	[m]	[%]	[m]	[%]	[m]	[%]
$4\times 4$	$-0.31692$	$4.80$	$-0.31767$	$5.05$	$-0.30825$	$1.93$	$-0.30832$	$1.96$	$-0.30367$	$0.42$
$8\times 8$	$-0.30391$	$0.50$	$-0.30443$	$0.67$	$-0.30180$	$0.20$	$-0.30185$	$0.18$	$-0.30069$	$0.56$
$16\times 16$	$-0.30158$	$0.27$	$-0.30173$	$0.22$	$-0.30083$	$0.52$	$-0.30085$	$0.51$	$-0.30056$	$0.61$
$32\times 32$	$-0.30113$	$0.42$	-	-	$-0.30071$	$0.56$	$-0.30072$	$0.56$	$-0.30065$	$0.58$
Reference solution: $u_Z\left(B\right)=-0.3024\text{m}$ ([29])

Table 12. Table 12: Deflection u Z ​ ( C ) subscript 𝑢 𝑍 𝐶 u_{Z}(C) for the pinched cylinder problem, thickness h = 0.03 ​ m ℎ 0.03 m h=0.03\,\text{m} , uniform mesh

Mesh	MITC4+		DKMQ24 [18]		DKMQ24		DKMQ241+		DKMQ242+
	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.
	[mm]	[%]	[mm]	[%]	[mm]	[%]	[mm]	[%]	[mm]	[%]
$4\times 4$	$-0.71240$	$60.03$	$-1.1255$	$39.26$	$-1.12416$	$39.33$	$-1.18399$	$36.10$	$-1.18474$	$36.06$
$8\times 8$	$-1.37736$	$24.15$	$-1.7246$	$6.93$	$-1.72382$	$6.97$	$-1.74588$	$5.78$	$-1.74647$	$5.75$
$16\times 16$	$-1.69944$	$6.77$	$-1.8599$	$0.37$	$-1.85952$	$0.35$	$-1.86512$	$0.65$	$-1.86557$	$0.68$
$32\times 32$	$-1.80254$	$1.20$	-	-	$-1.85688$	$0.21$	$-1.85939$	$0.34$	$-1.85951$	$0.35$
$64\times 64$	$-1.83392$	$0.49$	-	-	$-1.85151$	$0.08$	$-1.85444$	$0.08$	$-1.85445$	$0.08$
Reference solution: $u_Z\left(C\right)=-1.853\text{mm}$ (DKMQ242+ element on mesh $128\times 128$)

Table 13. Table 13: Deflection u Z ​ ( C ) subscript 𝑢 𝑍 𝐶 u_{Z}(C) for the pinched cylinder problem, thickness h = 0.03 ​ m ℎ 0.03 m h=0.03\,\text{m} , distorted mesh

Mesh	MITC4+		DKMQ24 [18]		DKMQ24		DKMQ241+		DKMQ242+
	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.
	[mm]	[%]	[mm]	[%]	[mm]	[%]	[mm]	[%]	[mm]	[%]
$4\times 4$	$-0.69215$	$61.13$	$-0.2397$	$87.06$	$-0.24254$	$86.91$	$-0.24973$	$86.52$	$-0.25396$	$86.29$
$8\times 8$	$-1.37499$	$24.27$	$-1.2326$	$33.46$	$-1.23106$	$33.56$	$-1.19506$	$35.51$	$-1.19549$	$35.48$
$16\times 16$	$-1.70090$	$6.69$	$-1.7144$	$7.48$	$-1.71369$	$7.52$	$-1.70823$	$7.81$	$-1.70857$	$7.79$
$32\times 32$	$-1.80400$	$1.12$	-	-	$-1.82068$	$1.74$	$-1.82087$	$1.73$	$-1.82104$	$1.72$
$64\times 64$	$-1.83392$	$0.49$	-	-	$-1.84221$	$0.58$	$-1.84481$	$0.44$	$-1.84484$	$0.44$
Reference solution: $u_Z\left(C\right)=-1.853\text{mm}$ (DKMQ242+ element on mesh $128\times 128$)

Table 14. Table 14: Deflection u Z ​ ( C ) subscript 𝑢 𝑍 𝐶 u_{Z}(C) for the pinched cylinder problem, thickness h = 0.3 ​ m ℎ 0.3 m h=0.3\,\text{m} , uniform mesh

Mesh	DKMQ24 [18]		DKMQ24		DKMQ241+		DKMQ242+
	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.
	[mm]	[%]	[mm]	[%]	[mm]	[%]	[mm]	[%]
$4\times 4$	$-0.011445$	$19.97$	$-0.011435$	$20.04$	$-0.011622$	$18.73$	$-0.011624$	$18.71$
$8\times 8$	$-0.012180$	$14.83$	$-0.012125$	$15.21$	$-0.012272$	$14.18$	$-0.012273$	$14.18$
$16\times 16$	$-0.012822$	$10.34$	$-0.012638$	$11.62$	$-0.012822$	$10.34$	$-0.012823$	$10.33$
$32\times 32$	-	-	$-0.013099$	$8.40$	$-0.013474$	$5.78$	$-0.013462$	$5.86$
$64\times 64$	-	-	$-0.013505$	$5.56$	$-0.014191$	$0.76$	$-0.014163$	$0.96$
Reference solution: $u_Z\left(C\right)=-0.01430\text{mm}$ (DKMQ242+ element on mesh $128\times 128$)

Table 15. Table 15: Deflection u Z ​ ( C ) subscript 𝑢 𝑍 𝐶 u_{Z}(C) for the pinched cylinder problem, thickness h = 0.3 ​ m ℎ 0.3 m h=0.3\,\text{m} , distorted mesh

Mesh	DKMQ24 [18]		DKMQ24		DKMQ241+		DKMQ242+
	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.
	[mm]	[%]	[mm]	[%]	[mm]	[%]	[mm]	[%]
$4\times 4$	$-0.008617$	$39.74$	$-0.008675$	$39.33$	$-0.008857$	$38.07$	$-0.008858$	$38.05$
$8\times 8$	$-0.011443$	$19.98$	$-0.011444$	$19.97$	$-0.011767$	$17.71$	$-0.011768$	$17.70$
$16\times 16$	$-0.012508$	$12.53$	$-0.012428$	$13.09$	$-0.012748$	$10.85$	$-0.012750$	$10.84$
$32\times 32$	-	-	$-0.013039$	$8.82$	$-0.013472$	$5.79$	$-0.013464$	$5.85$
$64\times 64$	-	-	$-0.013499$	$5.60$	$-0.014190$	$0.77$	$-0.014161$	$0.97$
Reference solution: $u_Z\left(C\right)=-0.01430\text{mm}$ (DKMQ242+ element on mesh $128\times 128$)

Table 16. Table 16: Deflection u Z ​ ( O ) subscript 𝑢 𝑍 𝑂 u_{Z}(O) for the hyperbolic paraboloid shell problem

Mesh	DKMQ24 [18]		DKMQ24		DKMQ241+		DKMQ242+
	$u_Z\left(O\right)$	err.	$u_Z\left(O\right)$	err.	$u_Z\left(O\right)$	err.	$u_Z\left(O\right)$	err.
	[mm]	[%]	[mm]	[%]	[mm]	[%]	[mm]	[%]
$2\times 2$	$0.35900$	$44.52$	$0.35889$	$44.47$	$0.35874$	$44.41$	$0.35891$	$44.48$
$4\times 4$	$0.27074$	$8.99$	$0.26938$	$8.44$	$0.27064$	$8.94$	$0.27055$	$8.91$
$8\times 8$	$0.25023$	$0.73$	$0.25030$	$0.76$	$0.25013$	$0.69$	$0.25019$	$0.71$
$16\times 16$	$0.24864$	$0.09$	$0.24866$	$0.10$	$0.24859$	$0.07$	$0.24868$	$0.10$
$32\times 32$	-	-	$0.24847$	$0.02$	$0.24845$	$0.01$	$0.24846$	$0.02$
Reference solution: $u_Z\left(O\right)=0.248417\text{mm}$ (DKMQ242+ element on mesh $128\times 128$)

Table 17. Table 17: Deflection u Y ​ ( A ) subscript 𝑢 𝑌 𝐴 u_{Y}(A) for the twisted beam problem with h = 0.0032 ​ m , F Y = 0.001 ​ N , F Z = 0 formulae-sequence ℎ 0.0032 m formulae-sequence subscript 𝐹 𝑌 0.001 N subscript 𝐹 𝑍 0 h=0.0032\,\text{m},~{}F_{Y}=0.001\,\text{N},~{}F_{Z}=0

Mesh	MITC4+		US-ATFHS8		DKMQ24 [18]		DKMQ24		DKMQ241+		DKMQ242+
	$u_Y\left(A\right)$	err.	$u_Y\left(A\right)$	err.	$u_Y\left(A\right)$	err.	$u_Y\left(A\right)$	err.	$u_Y\left(A\right)$	err.	$u_Y\left(A\right)$	err.
	[m]	[%]	[m]	[%]	[m]	[%]	[m]	[%]	[m]	[%]	[m]	[%]
$2\times 12$	$5.2718$	$0.30$	$5.6197$	$6.92$	$5.200$	$1.07$	$5.1798$	$1.45$	$4.3494$	$17.25$	$5.1128$	$2.72$
$4\times 24$	$5.2381$	$0.34$	$5.2665$	$0.20$	$5.217$	$0.74$	$5.2162$	$0.76$	$4.9987$	$4.90$	$5.2289$	$0.52$
$8\times 48$	$5.2423$	$0.26$	-	-	$5.238$	$0.34$	$5.2409$	$0.29$	$5.2334$	$0.43$	$5.2500$	$0.11$
$16\times 96$	$5.2444$	$0.22$	-	-	-	-	$5.2479$	$0.15$	$5.2499$	$0.12$	$5.2505$	$0.10$
Reference solution: $u_Y\left(A\right)=5.256\text{m}$ ([40])

Table 18. Table 18: Deflection u Z ​ ( A ) subscript 𝑢 𝑍 𝐴 u_{Z}(A) for the twisted beam problem with h = 0.0032 ​ m , F Y = 0 , F Z = 0.001 ​ N formulae-sequence ℎ 0.0032 m formulae-sequence subscript 𝐹 𝑌 0 subscript 𝐹 𝑍 0.001 N h=0.0032\,\text{m},~{}F_{Y}=0,~{}F_{Z}=0.001\,\text{N}

Mesh	MITC4+		US-ATFHS8		DKMQ24 [18]		DKMQ24		DKMQ241+		DKMQ242+
	$u_Z\left(A\right)$	err.	$u_Z\left(A\right)$	err.	$u_Z\left(A\right)$	err.	$u_Z\left(A\right)$	err.	$u_Z\left(A\right)$	err.	$u_Z\left(A\right)$	err.
	[m]	[%]	[m]	[%]	[m]	[%]	[m]	[%]	[m]	[%]	[m]	[%]
$2\times 12$	$1.2777$	$1.26$	$1.3762$	$6.35$	$1.274$	$1.55$	$1.2690$	$1.93$	$1.2644$	$2.28$	$1.2729$	$1.63$
$4\times 24$	$1.2874$	$0.51$	$1.2953$	$0.10$	$1.285$	$0.70$	$1.2859$	$0.62$	$1.2877$	$0.49$	$1.2898$	$0.32$
$8\times 48$	$1.2908$	$0.25$	-	-	$1.290$	$0.31$	$1.2912$	$0.22$	$1.2925$	$0.11$	$1.2926$	$0.11$
$16\times 96$	$1.2917$	$0.18$	-	-	-	-	$1.2926$	$0.11$	$1.2930$	$0.08$	$1.2930$	$0.08$
Reference solution: $u_Z\left(A\right)=1.294\text{m}$ ([40])

Table 19. Table 19: Deflection u Y ​ ( A ) subscript 𝑢 𝑌 𝐴 u_{Y}(A) for the twisted beam problem with h = 0.32 ​ m , F Y = 1000 ​ N , F Z = 0 formulae-sequence ℎ 0.32 m formulae-sequence subscript 𝐹 𝑌 1000 N subscript 𝐹 𝑍 0 h=0.32\,\text{m},~{}F_{Y}=1000\,\text{N},~{}F_{Z}=0

Mesh	MITC4+		DKMQ24 [18]		DKMQ24		DKMQ241+		DKMQ242+
	$u_Y\left(A\right)$	err.	$u_Y\left(A\right)$	err.	$u_Y\left(A\right)$	err.	$u_Y\left(A\right)$	err.	$u_Y\left(A\right)$	err.
	[m]	[%]	[m]	[%]	[m]	[%]	[m]	[%]	[m]	[%]
$2\times 12$	$5.3925$	$0.58$	$5.393$	$0.57$	$5.3762$	$0.88$	$5.4403$	$0.30$	$5.4411$	$0.32$
$4\times 24$	$5.4023$	$0.40$	$5.403$	$0.39$	$5.4013$	$0.42$	$5.4197$	$0.08$	$5.4198$	$0.08$
$8\times 48$	$5.4066$	$0.32$	$5.410$	$0.24$	$5.4116$	$0.23$	$5.4167$	$0.13$	$5.4167$	$0.13$
$16\times 96$	$5.4083$	$0.29$	-	-	$5.4151$	$0.16$	$5.4165$	$0.14$	$5.4165$	$0.14$
Reference solution: $u_Y\left(A\right)=5.424\text{m}$ ([40])

Table 20. Table 20: Deflection u Z ​ ( A ) subscript 𝑢 𝑍 𝐴 u_{Z}(A) for the twisted beam problem with h = 0.32 ​ m , F Y = 0 , F Z = 1000 ​ N formulae-sequence ℎ 0.32 m formulae-sequence subscript 𝐹 𝑌 0 subscript 𝐹 𝑍 1000 N h=0.32\,\text{m},~{}F_{Y}=0,~{}F_{Z}=1000\,\text{N}

Mesh	MITC4+		DKMQ24 [18]		DKMQ24		DKMQ241+		DKMQ242+
	$u_Z\left(A\right)$	err.	$u_Z\left(A\right)$	err.	$u_Z\left(A\right)$	err.	$u_Z\left(A\right)$	err.	$u_Z\left(A\right)$	err.
	[m]	[%]	[m]	[%]	[m]	[%]	[m]	[%]	[m]	[%]
$2\times 12$	$1.7245$	$1.68$	$1.624$	$7.41$	$1.6198$	$7.65$	$1.7486$	$0.31$	$1.7487$	$0.30$
$4\times 24$	$1.7428$	$0.64$	$1.711$	$2.45$	$1.7109$	$2.46$	$1.7503$	$0.21$	$1.7503$	$0.21$
$8\times 48$	$1.7479$	$0.35$	$1.740$	$0.80$	$1.7409$	$0.75$	$1.7514$	$0.15$	$1.7514$	$0.15$
$16\times 96$	$1.7493$	$0.27$	-	-	$1.7491$	$0.28$	$1.7518$	$0.13$	$1.7518$	$0.13$
Reference solution: $u_Z\left(A\right)=1.754\text{m}$ ([40])

Table 21. Table 21: Deflection u Z ​ ( A ) subscript 𝑢 𝑍 𝐴 u_{Z}(A) for Raasch’s hook

Mesh	DKMQ24		DKMQ241+		DKMQ242+
	$u_Z\left(A\right)$	err.	$u_Z\left(A\right)$	err.	$u_Z\left(A\right)$	err.
	[m]	[%]	[m]	[%]	[m]	[%]
$17\times 3$	$2.5852$	$48.57$	$4.6973$	$6.56$	$4.6997$	$6.51$
$34\times 5$	$2.3730$	$52.79$	$4.7895$	$4.72$	$4.8070$	$4.38$
$68\times 10$	$3.5430$	$29.52$	$4.9537$	$1.46$	$4.9629$	$1.27$
$136\times 20$	$5.9169$	$17.70$	$5.0712$	$0.88$	$5.0435$	$0.33$
$272\times 40$	$8.1284$	$61.70$	$5.1054$	$1.56$	$5.0496$	$0.45$
$544\times 80$	$9.1194$	$81.41$	$5.1146$	$1.74$	$5.0436$	$0.33$
Reference solution: $u_Z\left(A\right)=5.027\text{m}$ ([41])

Table 22. Table 22: Deflection u Z ​ ( C ) subscript 𝑢 𝑍 𝐶 u_{Z}(C) for the partly clamped hyperbolic paraboloid shell with h = 0.001 ​ m , ρ = 360 ​ kg ​ m − 3 formulae-sequence ℎ 0.001 m 𝜌 360 kg superscript m 3 h=0.001\,\text{m},~{}\rho=360\,\text{kg}\,\text{m}^{-3} , uniform mesh

Mesh	MITC4+		US-ATFHS8		DKMQ24		DKMQ241+		DKMQ242+
	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.
	[mm]	[%]	[mm]	[%]	[mm]	[%]	[mm]	[%]	[mm]	[%]
$4\times 8$	$-0.27433$	$4.68$	$-0.28875$	$0.33$	$-0.14864$	$48.35$	$-0.17433$	$39.43$	$-0.18961$	$34.12$
$8\times 16$	$-0.27499$	$4.45$	$-0.27896$	$3.07$	$-0.22297$	$22.53$	$-0.28434$	$1.20$	$-0.28502$	$0.97$
$16\times 32$	$-0.28095$	$2.38$	$-0.28230$	$1.91$	$-0.25861$	$10.14$	$-0.28877$	$0.34$	$-0.28883$	$0.36$
$32\times 64$	$-0.28509$	$0.94$	-	-	$-0.27085$	$5.89$	$-0.28799$	$0.07$	$-0.28800$	$0.07$
$64\times 128$	$-0.28691$	$0.31$	-	-	$-0.27472$	$4.55$	$-0.28782$	$0.01$	$-0.28783$	$0.01$
Reference solution: $u_Z\left(C\right)=-0.28780\text{mm}$ ([32])

Table 23. Table 23: Deflection u Z ​ ( C ) subscript 𝑢 𝑍 𝐶 u_{Z}(C) for the partly clamped hyperbolic paraboloid shell with h = 0.0001 ​ m , ρ = 3.6 ​ kg ​ m − 3 formulae-sequence ℎ 0.0001 m 𝜌 3.6 kg superscript m 3 h=0.0001\,\text{m},~{}\rho=3.6\,\text{kg}\,\text{m}^{-3} , uniform mesh

Mesh	MITC4+		DKMQ24		DKMQ241+		DKMQ242+
	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.
	[mm]	[%]	[mm]	[%]	[mm]	[%]	[mm]	[%]
$4\times 8$	$-0.24333$	$2.00$	$-0.01083$	$95.46$	$-0.00758$	$96.82$	$-0.02796$	$88.28$
$8\times 16$	$-0.23443$	$1.73$	$-0.01316$	$94.48$	$-0.09337$	$60.86$	$-0.11258$	$52.81$
$16\times 32$	$-0.23324$	$2.23$	$-0.03860$	$83.82$	$-0.22284$	$6.59$	$-0.22400$	$6.10$
$32\times 64$	$-0.23541$	$1.32$	$-0.09919$	$58.42$	$-0.23892$	$0.15$	$-0.23900$	$0.18$
$64\times 128$	$-0.23730$	$0.53$	$-0.17178$	$27.99$	$-0.23891$	$0.15$	$-0.23892$	$0.15$
Reference solution: $u_Z\left(C\right)=-0.23856\text{mm}$ ([32])

Table 24. Table 24: Deflection u Z ​ ( C ) subscript 𝑢 𝑍 𝐶 u_{Z}(C) for the partly clamped hyperbolic paraboloid shell with h = 0.001 ​ m , ρ = 360 ​ kg ​ m − 3 formulae-sequence ℎ 0.001 m 𝜌 360 kg superscript m 3 h=0.001\,\text{m},~{}\rho=360\,\text{kg}\,\text{m}^{-3} , distorted mesh

Mesh	MITC4+		DKMQ24		DKMQ241+		DKMQ242+
	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.
	[mm]	[%]	[mm]	[%]	[mm]	[%]	[mm]	[%]
$4\times 8$	$-0.21755$	$24.41$	$-0.06171$	$78.56$	$-0.07755$	$73.05$	$-0.07956$	$72.35$
$8\times 16$	$-0.27692$	$3.78$	$-0.18916$	$34.27$	$-0.22173$	$22.96$	$-0.22324$	$22.43$
$16\times 32$	$-0.28504$	$0.96$	$-0.26034$	$9.54$	$-0.28385$	$1.37$	$-0.28403$	$1.31$
$32\times 64$	$-0.28670$	$0.38$	$-0.27294$	$5.16$	$-0.28781$	$0.00$	$-0.28784$	$0.01$
$64\times 128$	$-0.28745$	$0.12$	$-0.27582$	$4.16$	$-0.28784$	$0.02$	$-0.28785$	$0.02$
Reference solution: $u_Z\left(C\right)=-0.28780\text{mm}$ ([32])

Table 25. Table 25: Deflection u Z ​ ( C ) subscript 𝑢 𝑍 𝐶 u_{Z}(C) for the partly clamped hyperbolic paraboloid shell with h = 0.0001 ​ m , ρ = 3.6 ​ kg ​ m − 3 formulae-sequence ℎ 0.0001 m 𝜌 3.6 kg superscript m 3 h=0.0001\,\text{m},~{}\rho=3.6\,\text{kg}\,\text{m}^{-3} , distorted mesh

Mesh	MITC4+		DKMQ24		DKMQ241+		DKMQ242+
	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.
	[mm]	[%]	[mm]	[%]	[mm]	[%]	[mm]	[%]
$4\times 8$	$-0.12648$	$46.98$	$-0.00120$	$99.50$	$-0.00161$	$99.33$	$-0.00173$	$99.27$
$8\times 16$	$-0.21096$	$11.57$	$-0.00751$	$96.85$	$-0.01770$	$92.58$	$-0.01807$	$92.42$
$16\times 32$	$-0.23703$	$0.64$	$-0.03676$	$84.59$	$-0.10256$	$57.01$	$-0.10314$	$56.77$
$32\times 64$	$-0.23856$	$0.00$	$-0.11378$	$52.31$	$-0.21383$	$10.37$	$-0.21411$	$10.25$
$64\times 128$	$-0.23844$	$0.0005$	$-0.18645$	$21.84$	$-0.23724$	$0.55$	$-0.23726$	$0.54$
Reference solution: $u_Z\left(C\right)=-0.23856\text{mm}$ ([32])

Table 26. Table 26: Deflection u Z ​ ( C ) subscript 𝑢 𝑍 𝐶 u_{Z}(C) for the plate with distorted mesh

Mesh	DKMQ24 [18]		DKMQ24		DKMQ241+		DKMQ242+
	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.	$u_Z\left(C\right)$	err.
	[m]	[%]	[m]	[%]	[m]	[%]	[m]	[%]
$2\times 2$	$-1920.3$	$28.37$	$-1910.8$	$27.73$	$-1725.3$	$15.33$	$-1725.3$	$15.33$
$4\times 4$	$-1616.1$	$8.03$	$-1906.6$	$7.39$	$-1560.0$	$4.27$	$-1560.0$	$4.27$
$8\times 8$	$-1538.0$	$2.81$	$-1529.2$	$2.22$	$-1517.6$	$1.45$	$-1517.6$	$1.45$
$16\times 16$	$-1519.0$	$1.54$	$-1510.6$	$0.98$	$-1507.7$	$0.79$	$-1507.7$	$0.79$
$32\times 32$	-	-	$-1506.1$	$0.68$	$-1505.4$	$0.63$	$-1505.4$	$0.63$
Reference solution: $u_Z\left(C\right)=-1495.95\text{m}$ ([42])

Equations

$${\beta }_x$$

$${\beta }_y$$

$$\mathbf{X}(r,s,t)=\sum _{i=1}^4{a}_i(r,s)\left[{\mathbf{X}}_i+t\frac{h}{2}{\mathbf{n}}_i\right],$$

$${\mathbf{X}}_{,r}(r,s,t)$$

$${\mathbf{X}}_{,s}(r,s,t)$$

$${\mathbf{X}}_{,t}(r,s,t)$$

$${\mathbf{X}}_r(r,s)$$

$${\mathbf{X}}_s(r,s)$$

$${\mathbf{v}}_3(r,s)$$

$${\mathbf{v}}_1(r,s)$$

$${\mathbf{v}}_2(r,s)$$

$$\displaystyle\bm{T}(r,s)=\left[\begin{array}[]{c}\bm{v}^{\text{T}}_{1}(r,s)\\ \bm{v}^{\text{T}}_{2}(r,s)\\ \bm{v}^{\text{T}}_{3}(r,s)\end{array}\right],~{}~{}\bm{T}^{-1}(r,s)=\bm{T}^{\text{T}}(r,s)=[\bm{v}_{1}(r,s)~{}\bm{v}_{2}(r,s)~{}\bm{v}_{3}(r,s)].$$

$$\displaystyle\left[\begin{array}[]{c}x\\ y\\ z\end{array}\right]=\bm{T}\left[\begin{array}[]{c}X\\ Y\\ Z\end{array}\right],~{}\left[\begin{array}[]{c}u_{x}\\ u_{y}\\ u_{z}\end{array}\right]=\bm{T}\left[\begin{array}[]{c}u_{X}\\ u_{Y}\\ u_{Z}\end{array}\right],~{}\left[\begin{array}[]{c}\varphi_{x}\\ \varphi_{y}\\ \varphi_{z}\end{array}\right]=\bm{T}\left[\begin{array}[]{c}\varphi_{X}\\ \varphi_{Y}\\ \varphi_{Z}\end{array}\right],$$

$$\displaystyle\left[\begin{array}[]{c}X\\ Y\\ Z\end{array}\right]=\bm{T}^{\text{T}}\left[\begin{array}[]{c}x\\ y\\ z\end{array}\right],~{}\left[\begin{array}[]{c}u_{X}\\ u_{Y}\\ u_{Z}\end{array}\right]=\bm{T}^{\text{T}}\left[\begin{array}[]{c}u_{x}\\ u_{y}\\ u_{z}\end{array}\right],~{}\left[\begin{array}[]{c}\varphi_{X}\\ \varphi_{Y}\\ \varphi_{Z}\end{array}\right]=\bm{T}^{\text{T}}\left[\begin{array}[]{c}\varphi_{x}\\ \varphi_{y}\\ \varphi_{z}\end{array}\right].$$

$$\displaystyle\left[\begin{array}[]{c}u_{X}\\ u_{Y}\\ u_{Z}\end{array}\right](r,s,t)=$$

$$t\frac{h}{2}\sum _{k=5}^8{a}_k(r,s)△{\beta }_{tk}{\mathbf{t}}_k,t\in [-1,1],$$

$$\displaystyle\bm{E}_{i}=-\text{Spin}(\bm{v}^{i}_{3})=\left[\begin{array}[]{ccc}~{}~{}0&~{}~{}v^{i}_{3Z}&-v^{i}_{3Y}\\ -v^{i}_{3Z}&~{}~{}0&~{}~{}v^{i}_{3X}\\ ~{}~{}v^{i}_{3Y}&-v^{i}_{3X}&~{}~{}0\\ \end{array}\right].$$

$$\displaystyle\left[\begin{array}[]{c}\triangle u_{iX}\\ \triangle u_{iY}\\ \triangle u_{iZ}\\ \end{array}\right]=\,\,$$

$$\displaystyle z\!\left[\begin{array}[]{ccc}-v^{i}_{2X}&v^{i}_{1X}&0\\ -v^{i}_{2Y}&v^{i}_{1Y}&0\\ -v^{i}_{2Z}&v^{i}_{1Z}&0\\ \end{array}\right]\underbrace{\left[\begin{array}[]{ccc}v^{i}_{1X}&v^{i}_{1Y}&v^{i}_{1Z}\\ v^{i}_{2X}&v^{i}_{2Y}&v^{i}_{2Z}\\ v^{i}_{3X}&v^{i}_{3Y}&v^{i}_{3Z}\\ \end{array}\right]}_{\bm{T}(r_{i},s_{i})}\left[\begin{array}[]{c}\varphi_{iX}\\ \varphi_{iY}\\ \varphi_{iZ}\end{array}\right]=$$

$$\displaystyle z\!\underbrace{\left[\begin{array}[]{ccc}0&~{}~{}v^{i}_{3Z}&-v^{i}_{3Y}\\ -v^{i}_{3Z}&0&~{}~{}v^{i}_{3X}\\ ~{}~{}v^{i}_{3Y}&-v^{i}_{3X}&0\\ \end{array}\right]}_{\bm{E}_{i}}\left[\begin{array}[]{c}\varphi_{iX}\\ \varphi_{iY}\\ \varphi_{iZ}\end{array}\right],~{}~{}z=t\frac{h}{2},~{}~{}t\in[-1,1],~{}~{}i=1,2,3,4,$$

$${\beta }_t(t)$$

$${\beta }_m(t)$$

$$\mathbf{J}=\frac{\partial (X,Y,Z)}{\partial (r,s,t)}$$

$${\mathbf{X}}_r$$

$${\mathbf{X}}_s$$

$$\mathbf{m}$$

$${\mathbf{m}}_r$$

$${\mathbf{m}}_s$$

$$\overline{\mathbf{J}}:=\underset{t\to 0}{\lim }\frac{\partial (r,s,t)}{\partial (X,Y,Z)}$$

$$\frac{\text{d}{\mathbf{A}}^{-1}}{\text{d}t}=-{\mathbf{A}}^{-1}\frac{\text{d}\mathbf{A}}{\text{d}t}{\mathbf{A}}^{-1},$$

$${\overline{\mathbf{J}}}^{\prime }$$

$$-\overline{\mathbf{J}}\left[\underset{t\to 0}{\lim }\frac{\text{d}\mathbf{J}}{\text{d}t}\right]\overline{\mathbf{J}}=-\overline{\mathbf{J}}\left[{\mathbf{m}}_r,{\mathbf{m}}_s,0\right]\overline{\mathbf{J}},$$

$$\displaystyle\overline{\bm{J}}^{\prime}=\left[\begin{array}[]{c}\bm{o}_{r}^{\text{T}}\\ \bm{o}_{s}^{\text{T}}\\ \bm{o}_{t}^{\text{T}}\\ \end{array}\right].$$

$$\frac{\partial (u_X,u_Y,u_Z)}{\partial (r,s,t)}$$

$${\mathbf{u}}_r$$

$${\mathbf{u}}_s$$

$$\mathbf{\Phi }$$

$${\mathbf{\Phi }}_r$$

$${\mathbf{\Phi }}_s$$

Full text

††journal: Finite Elements in Analysis and Design

Nomenclature

Abbreviations

[TABLE]

Latin symbols

[TABLE]

[TABLE]

Greek symbols

[TABLE]

1 Introduction

Finite element modeling of shell structures is of great practical importance in engineering practice. Therefore there is a need for an efficient shell element, which would have an optimal convergence rate for different kinds of problems (membrane dominated, bending dominated and mixed shell problems), which would be at the same time free of shear locking, free of membrane locking, without zero energy modes and which would keep the optimal convergence rate also on distorted meshes. Especially, it is desired to find an optimal low-order shell element, which would provide the highest benefit-cost ratio when used in practical engineering applications. To find such an element is a tremendous task, not solved until today.

Let us restrict ourselves in this paper to 4-nodal quadrangle shell elements, based on the linear approximation of the transversal shear (the Reissner-Mindlin concept), which use both translational and rotational degrees of freedom in a node. Moreover, let us concentrate on the case of small deformations, however, an extension to the case of large deformations, keeping the small strain assumption, can be easily made in the corotational approach [1, 2, 3], in which no change of the element formulation is needed.

In order to construct a quadrangle shell element, it is possible to either consider a flat element geometry, by which however it is necessary to use the warp correction term suggested by Taylor [4, 5] in the element formulation, or a warped element geometry. Formulations of flat quadrangles can be found in [6, 3, 7, 8, 9]. Formulations using a non-flat geometry include the QUAD4 element suggested by MacNeal [10] and the famous MITC4 (Mixed Interpolation Tensorial Components) element proposed by Dvorkin and Bathe [11], which was further extended to the so called MITC4+ element [12, 13]. Another extension of the MITC4 element called MITC4S was introduced by Niemi [14]. As alternative, shells can be modelled using solid-shell elements, as for example the US-ATFHS8 [15], which is a 8-nodal element with 24-dofs (displacement dofs only). Other recent shell formulations can be found in Gruttmannn [16], Kim [17] or Katili [18]. The latter one introduces the DKMQ24 shell element, based on the Naghdi/Reissner/Mindlin shell theory with six degrees of freedom per node. This shell element is an extension of the DKMQ (Discrete Kirchhoff-Mindlin Quadrangle) plate element suggested by the same author, which proved to have good numerical properties on a wide set of problems [19]. The DKMQ24 shell element has full coupling of membrane-bending energy, is free of shear locking, free of membrane locking and proved to converge on a wide set of standard benchmark problems [18]. It was successfully applied to orthotropic multilayered shells [20] and to a beam analysis [21].

The DKMQ24 shell element can be, due to its good numerical properties, considered as a promising candidate for an efficient low order quadrangle shell element. However, despite its numerical efficiency, it suffers from one drawback – the drilling rotations are independent of the remaining degrees of freedom. The stiffness matrix of the DKMQ24 shell element has a full rank due to stabilization, however, the drilling rotations do not enrich the displacement interpolation space while increasing computational cost at the same time. One possibility to avoid this inefficiency would be the use of five degrees of freedom per node only, which was proposed by Irpanni [22]. However, from a practical point of view six degrees of freedom at a node is of benefit in order to easily connect more arbitrary oriented shell (or beam) elements to a single node.

The aim of our contribution is to modify the formulation of the DKMQ24 element in order to overcome this imperfection. We propose the incorporation of drilling rotations with the help of Allman’s shape functions [23] in the element formulation. This enriches the in-place displacement field with incomplete quadratic functions and improves the membrane behaviour of the element. The element formulation contains a penalty term, forcing the drilling rotations to be equal to an in-plane rotation of the displacement field. This term is multiplied by the penalty coefficient, which is either constant or scaled to the shell thickness $h$. A further improvement of the element formulation is achieved by the selective reduced integration scheme of shear terms, as described in Section 3.4. Moreover, nodal moment corrections for the case of distributed normal loading are derived in Section 3.5. The modified element formulation is named either DKMQ241+ (for a constant penalty coefficient) or DKMQ242+ (for a penalty coefficient scaled to the shell thickness $h$). The formulation was tested on nine well-established benchmark problems, which show that it is free of shear locking, free of membrane locking, converges in all considered test cases, and provides a reduced numerical error in the membrane-dominated problems.

The practical implementation of the DKMQ24, DKMQ241+ and DKMQ242+ shell elements was done in the C# programming language using the Microsoft .NET Framework technology as a part of the finite element program femCalc [24].

In Section 2 we give a full and detailed derivation of the DMKQ24 shell element introduced in [19]. In Section 3 modifications of the DKMQ24 element, which construct the DKMQ241+ and DKMQ242+ elements, are introduced. The dependence of the relative error at the selected benchmarks on the penalty constant of the DKMQ24, DKMQ241+ and DKMQ242+ elements is discussed in Section 4. In Section 5, convergence behaviour of the DKMQ24, DKMQ241+ and DKMQ242+ elements on pure membrane, pure bending and on mixed membrane-bending benchmark problems is tested and discussed. Closing remarks are given in Section 6.

2 Formulation of the DKMQ24 shell element

In this Section a detailed derivation of the DKMQ24 shell element, introduced in [18], is given. The reason for this decision is that in reference [18] some parts of the derivation are not presented in detail.

2.1 Element characterization

The presented shell element is four-nodal with six unknowns at each node (three displacements $[u_{X},u_{Y},u_{Z}]$ and three rotations $[\varphi_{X},\varphi_{Y},\varphi_{Z}]$ in the global c.s.). The element is considered to be warped and considers full coupling between membrane and bending terms. Linear elasticity theory is assumed throughout the whole article.

2.2 Definition of rotations

Rotations $\varphi_{x},\varphi_{y},\varphi_{z}$ are defined using the standard definition, i.e. as rotations about axes $x,y,z$ using the right-hand-rule. Rotations $\varphi_{x},\varphi_{y}$ can be equivalently defined as rotations of axis $z$ towards the axis $x$ or $y$ respectively, denoted by $\beta_{x},\beta_{y}$. The relation between these definitions is the following (see Figure 1):

[TABLE]

We use rotations $\beta_{x},\beta_{y}$ during the derivation of the element, because it is more convenient. However, the final formulation of the DKMQ24 shell element is given in the standard nodal rotations $\varphi_{x},\varphi_{y},\varphi_{z}$.

2.3 Coordinate system

The geometry of an isoparametric curvilinear shell element in the reference configuration is described by

[TABLE]

where $\bm{X}_{i}=[X_{i},Y_{i},Z_{i}]$ is the location of the $i$-th node in the reference (undeformed) configuration, $h$ is a shell thickness and $\bm{n}_{i}$ is a given normal vector at node $i$, see Figure 2. The shell thickness $h$ in Eq. (3) is assumed to be constant in derivation of the strain matrices (similar to the approach used in [25]). The actual thickness and its variation within the element is introduced into the constitutive relations only (Eq. (217, 221, 224)).

The basis functions are summarized in Table 1. At this point we can calculate a covariant vector basis with respect to all parameters

[TABLE]

where we have used the notation $\bullet_{,r}:=\frac{\partial\bullet}{\partial r},~{}\bullet_{,s}:=\frac{\partial\bullet}{\partial s}$. These vectors are in general neither orthogonal nor orthonormal. Let us denote

[TABLE]

Let us define a normal vector at any point $(r,s)$ by

[TABLE]

with the help of which we define normal vectors at nodes $\bm{n}_{i}=\bm{v}_{3}(r_{i},s_{i}),~{}i=1,\dots,4$. These vectors are used in Eq. (3). The remaining orthonormal vectors $\bm{v}_{1}(r,s)$ and $\bm{v}_{2}(r,s)$ can be obtained simply by

[TABLE]

where the vectors $\bm{v}_{1},~{}\bm{v}_{2}$ and $\bm{v}_{3}$ comprise the local coordinate system. If an orthotropic material of the shell is considered, the vector $\bm{v}_{1}$ is usually chosen to be collinear with the direction of fibers.

2.4 Transformation between global and local coordinate systems

The local orthonormal coordinate system $\bm{v}_{1}(r,s),\bm{v}_{2}(r,s),\bm{v}_{3}(r,s)$ forms the transformation matrix $\bm{T}(r,s)$

[TABLE]

The same transformation rules hold for coordinates, displacements and rotations

[TABLE]

where upper-case letters refer to the reference coordinate system and lower-case letters refer to the local coordinate system. The matrix which comprises the orthonormal basis is orthogonal, therefore inversion of above relations is trivial:

[TABLE]

2.5 Displacements

Let us denote the orthonormal basis vectors $\bm{v}_{1}(r,s),\bm{v}_{2}(r,s),\bm{v}_{3}(r,s)$ evaluated at $i$-th node simply by $\bm{v}_{1}^{i},\bm{v}_{2}^{i},\bm{v}_{3}^{i}$. The displacement field of the DKMQ24 shell element is then defined by

[TABLE]

where

[TABLE]

We have used the definition of the spin matrix $\text{Spin}(\bm{v})$ taken from [1]. The supplementary quadratic basis functions $a_{5},a_{6},a_{7},a_{8}$, which enrich the bending behaviour of the element, are summarized in Table 2.

The second term in Eq. (64) represents the displacement increment caused by nodal rotations. Let us derive this term for the $i$-th node:

[TABLE]

where we have used the definition of the vector multiplication $\bm{v}_{3}=\bm{v}_{1}\times\bm{v}_{2}$. The last term in Eq. (64) contains supplementary degrees of freedom $\triangle\beta_{tk}$ (see Figure 3). The supplementary degrees of freedom enrich the rotation approximation in the tangential $t$-direction, where the quadratic approximation is considered. In the perpendicular $m$-direction to side $k$ only a linear approximation of rotation $\beta_{m}$ is used:

[TABLE]

The supplementary degrees of freedom $\triangle\beta_{tk}$ will be eliminated later on the element level, so that the only element unknowns remain three rotations and three displacements at each node (see Figure 4).

2.6 Jacobi matrix

The derivative of the nodal location in the global coordinate system (3) with respect to parameters $r,s,t$ yields

[TABLE]

where

[TABLE]

The Jacobi matrix required for the calculation of the membrane stiffness term yields

[TABLE]

For the calculation of the bending stiffness term the following term is also needed: $\lim_{t\to 0}\frac{\textrm{d}}{\textrm{d}t}\frac{\partial(r,s,t)}{\partial(X,Y,Z)}$. By recalling the matrix identity

[TABLE]

which follows from the derivation of the identity matrix $\bm{I}=\bm{A}\bm{A}^{-1}$, we get

[TABLE]

where upper indices denote dual vectors111For the set of linearly independent vectors $\bm{u}_{1},\dots,\bm{u}_{N}$ one can uniquely define a set of dual vectors $\bm{u}^{1},\dots,\bm{u}^{N}$, which satisfy $\bm{u}_{i}^{\text{T}}\bm{u}^{j}=\delta_{ij},~{}i,j=1,\dots,N$, where $\delta_{ij}$ is the Kronecker delta symbol.. Rows of the $\overline{\bm{J}}^{\prime}$ matrix are explicitly denoted by vectors $\bm{o}_{r}^{\text{T}},\bm{o}_{s}^{\text{T}},\bm{o}_{t}^{\text{T}}$:

[TABLE]

2.7 Displacement derivatives

The derivative of the displacement field (64) with respect to parameters $r,s,t$ yields

[TABLE]

where

[TABLE]

2.8 Membrane and bending strain matrices

The vector of element unknowns $\bm{q}$ of size $24\times 1$ is composed in the following way

[TABLE]

The membrane and bending strain matrices can be calculated in the local coordinate system using the strain definition

[TABLE]

Using the chain rule we get

[TABLE]

and analogously for the remaining derivatives $\frac{\partial u_{y}}{\partial x}$ and $\frac{\partial u_{y}}{\partial y}$. The membrane strain $\bm{\varepsilon}$ yields

[TABLE]

where

[TABLE]

where the columns are indexed by index $j$. The bending strain $\bm{\kappa}$ yields

[TABLE]

where

[TABLE]

[TABLE]

where the columns are indexed by index $j$ and $k$ respectively (no summation over these indices is considered).

2.9 Strain energy and constitutive relations

The expression for the strain energy of the element reduces to the standard form:

[TABLE]

where the vector of unknowns $\bm{q}$ is given in Eq. (143). The stiffness matrix of an element $e$ is given by the formula

[TABLE]

where

[TABLE]

and $A$ is an area of the element $e$, $\bm{B}_{\text{b}}$ is the bending strain matrix defined later in Eq. (270), $\bm{B}_{\text{m}}$ is the membrane strain matrix defined in Eq. (171) and $\bm{K}_{\text{stab}}$ is a matrix stabilizing drilling rotations specified later in Eq. (351). The elasticity matrices in Eq. (213) equal in case of single-layered shells made of isotropic material to

[TABLE]

Actual thickness and its variation within the element is accounted for in the calculation of $\bm{D}_{\text{b}},~{}\bm{D}_{\text{m}},~{}\bm{D}_{\text{s}}$ matrices. The strain resultants evaluated in the Gauss quadrature points are

[TABLE]

and analogously the stress resultants evaluated in the Gauss quadrature points are

[TABLE]

The stress resultants $\bm{m}$, $\bm{n}$ are defined in the local coordinate system defined by vectors $\bm{v}_{1},~{}\bm{v}_{2}$. If it is desired to get the stress resultants in other coordinate system $\bm{v}_{1}^{\prime},~{}\bm{v}_{2}^{\prime}$, we use the transformation

[TABLE]

where $T_{ij}=(\bm{v}_{i}^{\prime})^{\text{T}}\bm{v}_{j},~{}i=1,2$. The bending strain matrix $\bm{B}_{\text{b}}$ and the shear strain matrix $\bm{B}_{\text{s}}$ are derived to depend also on supplementary degrees of freedom $\triangle\beta_{tk}$. These supplementary degrees of freedom are then eliminated with the help of the matrix $\bm{A}_{\text{n},4\times 24}$, which relates them to nodal degrees of freedom as follows:

[TABLE]

The assembly of the bending stiffness matrix is then given by

[TABLE]

where

[TABLE]

The shear stiffness matrix takes the form

[TABLE]

2.10 Evaluation of shear matrix

Matrix $\bm{B}_{\text{s}\triangle\beta,2\times 4}$, which relates the transversal shear strains to the supplementary rotations $\triangle\beta_{tk}$, is identical to the same matrix used in the DKMQ plate element [19]. We recap the derivation for completeness.

Let us evaluate an average transversal shear strain on each side. Taking an arbitrary side $k$, the directional parameter is denoted by $t$. The transversal shear strain in direction $t$ has the form

[TABLE]

where $\mathfrak{q}_{t}$ is a shear internal force and $D_{\text{s}}$ is the shear stiffness, which in case of single-layered isotropic shells equals to $\frac{5}{6}Gh$, where $G$ is the shear modulus and $h$ is an average shell thickness on side $k$. The stress equilibrium yields

[TABLE]

and Hooke’s law has the form

[TABLE]

where $D_{\text{b}}$ is an average bending stiffness on the side $k$, which in case of single-layered isotropic shell equals to $\frac{Eh^{3}}{12(1-\nu^{2})}$, where $E$ is Young’s modulus and $h$ the shell thickness. We now consider a linear approximation of angle $\beta$ in normal $m$-direction and quadratic in tangential $t$-direction according to Eq. (96). Both approximations do not depend on the $m$-direction, which immediately yields

[TABLE]

where

[TABLE]

We see that the transversal shear approximation $\overline{\gamma}_{tz}$ is constant on the side. If we liked to transform these shear strains back to the reference plane, we would have to incorporate the determinant of the transformation and also the parametrization sign change, see to Figure 5:

[TABLE]

As in the MITC4 element [11] the transverse shear strains are approximated with

[TABLE]

where

[TABLE]

so we finally get

[TABLE]

2.11 Evaluation of matrix A${}_{\text{n}}$

We use the Hu-Washizu functional according to which

[TABLE]

where $\gamma_{tz}$ is a shear approximated from displacement and $\overline{\gamma}_{tz}$ is an average shear approximated by Eq. (285). We have

[TABLE]

By rewriting the scalar unknowns to vector ones we get

[TABLE]

In the matrix form the Eq. (311) takes the form

[TABLE]

where $\bm{n}_{k}$ is the normal vector corresponding to the side $k$ starting with node $i$ and ending with node $j$. Finally, the supplementary rotations $\triangle\beta_{tk}$ can be evaluated by the following equation

[TABLE]

2.12 Drilling rotations stabilization

The shell element must be stabilized with respect to drilling rotations $\varphi_{z}$ in order to avoid a singular stiffness matrix, which appears if a node connects coplanar elements. We follow the stabilization procedure given in [18]. The added energy, which stabilizes drilling rotation has the form

[TABLE]

where $c$ is the penalty constant. The value $c=0.001$ is considered in our implementation of the DKMQ24 element as well as in [18]. The corresponding stiffness matrix in the global coordinate system has the form:

[TABLE]

The first term is called the MacNeal stabilization and has rank 1, the second has rank 3. It can be proven that such a matrix removes all zero energy modes connected with drilling rotations [18].

2.13 Numerical integration

Elements are integrated by means of the standard $2\times 2$ Gauss quadrature except of the stabilization term (351), for which a reduced one-point integration rule is used.

2.14 Load vector

The work of external forces for an arbitrary element $e$ in the discretized form yields

[TABLE]

In order to evaluate the consistent load vector $\bm{f}$ let us evaluate the deformation of the middle surface using Eq. (64) by setting $t=0$, which yields

[TABLE]

Let us further consider a general loading case given by tangential pressures $p_{x}(r,s)$ and $p_{y}(r,s)$, by the normal pressure $p(r,s)$ and by the element loading force $\bm{F}_{\text{e}}$. The consistent nodal force load vector in the global coordinate system is given by the integration of the basis functions given by Eq. (371) and yields

[TABLE]

where

[TABLE]

where $\bm{v}_{1}(r,s),\bm{v}_{2}(r,s),\bm{v}_{3}(r,s)\equiv\bm{n}(r,s)$ are the basis vectors, $A$ is the element area and $w_{q},~{}q=1,2,3,4$ are the weights and $r_{q},\,s_{q},~{}q=1,2,3,4$ are the quadrature points of the Gauss quadrature $2\times 2$.

3 Improvements to the DKMQ24 shell element

In this Section we introduce modifications of the DKMQ24 shell element, derived in Section 2. These modifications apply both to the DKMQ241+ and DKMQ242+ shell elements. The difference between these two variants is explained in Section 3.3.

3.1 The treatment of drilling rotations

In the original derivation drilling rotations are stabilized according to formula (350), which causes the local element stiffness matrix to have a full rank even for the coplanar element configuration. However, this treatment keeps drilling rotations independent from the remaining degrees of freedom and as a consequence the element behaviour is not improved although we pay for these additional degrees of freedom by an increase of computational costs. Therefore, we use the approach proposed in [27], where the drilling rotations are considered as derivatives of the in-plane displacement field. The space of basis functions is enriched by incomplete quadratic polynomials in this case.

Let us start with calculation of a rotation around a normal vector $\mathbf{n}$ at an arbitrary point $(r,s)$

[TABLE]

The displacement field \bm{u}=\left[\begin{array}[]{c}u_{X}(r,s)\\ u_{Y}(r,s)\\ u_{Z}(r,s)\end{array}\right] given by formula (64) is then enriched by the following term

[TABLE]

where the starting node $I(k)$ and the ending node $J(k)$ for the side $k=5,6,7,8$ are given by

[TABLE]

and $\bm{l}_{k}$ is the in-plane normal vector corresponding to the side $k$, see Figure 4, given by $\bm{l}_{k}=\frac{\bm{n}_{k}\times\bm{t}_{k}}{||\bm{n}_{k}\times\bm{t}_{k}||}$, where $\bm{n}_{k}=\frac{\bm{n}_{I(k)}+\bm{n}_{J(k)}}{2},~{}\bm{t}_{k}=\frac{P_{J(k)}-P_{I(k)}}{|P_{J(k)}-P_{I(k)}|}$, see Figure 6. The membrane strain field \bm{\varepsilon}=\left[\begin{array}[]{c}\varepsilon_{x}(r,s)\\ \varepsilon_{y}(r,s)\\ \gamma_{xy}(r,s)\end{array}\right] given by formula (167) is enriched by the term \left[\begin{array}[]{c}F_{11}\\ F_{22}\\ F_{12}+F_{21}\end{array}\right], where

[TABLE]

where

[TABLE]

where the vector $\mathbf{e}_{i}$ is a zero vector with 1 at index $i$. The bending strain field $\bm{\kappa}$ given by Eq. (180) stays unchanged. To conclude this Section, let us remark that the DKMQ24 shell element has improved bending behaviour, due to quadratic enrichment of the displacement field with the help of out-of-plane rotations $\varphi_{x},\varphi_{y}$ and in the same manner, the DKMQ241+ and DKMQ242+ shell elements have improved membrane behaviour due to quadratic enrichment of the displacement field with the help of drilling rotations $\varphi_{z}$.

3.2 The static condensation

We enrich the displacement field by a bubble mode according to the approach given in [27]. Two new degrees of freedom $u_{1},u_{2}$ corresponding to the mid-element node $9$ are considered:

[TABLE]

The corresponding strain vector has the form

[TABLE]

These additional degrees of freedom are statically condensed out on the element level, therefore the shell element has again 24 degrees of freedom. The static condensation is based on a standard procedure applied on the membrane part only (Appendix A). The statically condensed strain matrix $\bm{B}$ takes the form

[TABLE]

where

[TABLE]

3.3 Penalty term

The drilling rotations stabilization term (350), used in the case of the DKMQ24 shell element, is not used in case of the DKMQ241+ and DKMQ242+ shell elements. However, a penalty term which forces the drilling rotations $\varphi_{z}(r,s)$ to be equal to a rotation of the displacement field $\psi$ at the same point, derived by Hughes and Brezzi [27], is considered. This energy term takes the form

[TABLE]

where

[TABLE]

and $c_{1}$ is the penalty constant and is taken as $c_{1}=0.15$ in our implementation. The corresponding stiffness matrix in the global coordinate system has the form:

[TABLE]

where the vector $\bm{g}$ has length 24 and takes the form

[TABLE]

The one point integration rule is used to integrate the stiffness matrix (425) in order to relax the prescribed condition (409).

We use either a constant penalty parameter $c_{1}$ (let us call this element variant DKMQ241+) or we propose the following scaling

[TABLE]

where $h$ is the element thickness and $A$ is the element area. This makes the scaling dependent only on the local information and independent of length units. Moreover, this scaling increase the penalty constant for thicker and smaller elements (let us call this element variant DKMQ242+).

3.4 Numerical integration of shear terms

By using the following explicit notation of the shear strain matrix rows

[TABLE]

we can decompose the element stiffness matrix $\bm{K}_{\text{e}}$, given by Eq. (208), to the form

[TABLE]

where the stabilisation stiffness matrix $\bm{K}_{\text{stab}}$ is given by Eq. (425). We have numerically tested different quadrature rules for each term in Eq. (431) and we suggest the reduced integration of shear terms, summarized in Table 3. The proposed reduced integration provides smaller numerical errors when compared to the $2\times 2$ quadrature rule, mainly at the pinched cylinder benchmark in variant $h=0.3\,\text{m}$ (Section 5.4). Also, the stiffness matrix assembling is slightly faster, because the number of quadrature points for the shear terms is reduced to half.

3.5 Nodal moment corrections

Note that due to absence of rotations in Eq. (371) the load vector does not contain any nodal moment corrections. The reason is the Mindlin–Reissner hypothesis of the formulation, which separates high order terms of the displacement to the rotations. The real (physical) displacement of the element middle shell surface can be, analogously to Eq. (384), considered in the following way

[TABLE]

Using the normal projection of the pressure and force loading introduced in Eq. (373), denoted by $\lambda(r,s)$, given by formula

[TABLE]

and considering the work of the normal forces $\delta W=\lambda(r,s)\delta u_{z}(r,s)$, the second term in Eq. (444) yields the nodal moment corrections corresponding to the distributed load described by Eq. (373) of the form

[TABLE]

where $w_{q},~{}q=1,2,3,4$ are the weights and $r_{q},\,s_{q},~{}q=1,2,3,4$ are the quadrature points of the Gauss quadrature $2\times 2$. In the component form we finally get

[TABLE]

These nodal moment corrections are used in the shell elements DKMQ241+ and DKMQ242+.

4 Dependence on penalty constants

In this Section, the dependence of the relative error at the selected benchmark problems on the penalty constant $c$ (in case of the DKMQ24 shell element), on the penalty constant $c_{1}$ (in case of the DKMQ241+ shell element) and the penalty constant $c_{2}$ (in case of the DKMQ242+ shell element) is investigated. In general, the benchmark problems can be divided in three categories: The first category contains problems whose error is higher for higher values of the penalty constant. The second category, which is in our case presented only by Raasch’s hook, is opposite – a sufficient value of the penalty constant is needed in order to get a precise solution. The last and the largest category of benchmark problems exhibits no or negligible dependence on the penalty constant. We have selected some representative benchmark problems, so that all three categories are included in our comparison.

Let us start with the DKMQ24 shell element in our implementation, see Fig. 7. Raasch’s hook allows only a small range of the penalty constant, approx. $c\in[10^{-2},10^{-1}]$, in which the relative error is acceptible. On the other hand, the hemispherical shell benchmark with $16\times 16$ mesh belongs to the first category – the penalty constant needs to satisfy $c<10^{-4}$ in order to get acceptible accuracy at this benchmark. In conclusion, there is no satisfactory choice of the penalty constant $c$ for this element. Let us now focus on the DKMQ241+ shell element, see Fig. 8. Raasch’s hook benchmark problem needs the penalty constant $c_{1}$ to be higher then approx. $10^{-1}$ to get a reasonable error. The twisted beam benchmark problem (with setting $h=0.0032\,\text{m},F_{Y}>0$) provides the lower error for the lower values of the penalty constant $c_{1}$. Therefore, the value $c_{1}=0.15$, used in our calculation is near the optimal value. Let us now focus on the shell element DKMQ242+, see Fig. 9. The situation here is qualitatively similar to the situation of the shell element DKMQ241+, only the relative errors are lower, see Fig. 9. Raasch’s hook benchmark problem needs the penalty constant $c_{2}$ to be higher then approx. $0.05$. The twisted beam benchmark problem (with setting $h=0.0032\,\text{m},F_{Y}>0$) needs the penalty constant $c_{2}$ not to be higher then approx. $0.1$. Therefore, the value $c_{2}=0.1$ used in our calculation is near the optimal value.

By comparison of the DKMQ241+ and DKMQ242+ shell elements at all benchmark problems given in Section 5, we see that they provide almost an identical error except for two cases: the twisted beam benchmark problem ($h=0.0032$) and the Raasch’s hook benchmark problem, at both of which the DKMQ242+ element is superior. Therefore these two benchmarks were selected into our comparison in this Section. Moreover, we therefore recommend the DKMQ242+ shell element over the DKMQ241+ shell element.

5 Convergence tests

The behaviour of the DKMQ24, DKMQ241+ and DKMQ242+ shell elements has been tested on nine well-established benchmark problems, including one pure membrane benchmark problem suggested by Cook [28], six shell and one plate benchmark problems suggested by MacNeal and Harder [29] and Belytschko [30], which were reused by Katili [18], and one challenging benchmark proposed by Raasch [31]. We compare the results of the DKMQ241+ and DKMQ242+ shell elements with the recent state-of-the-art quadrangle shell/solid-shell elements, which are all summarized in Table 4. In all cases, geometric and material linearity is considered and there is no initial deformation.

5.1 Cook’s membrane

The standard Cook’s membrane problem tests an element behaviour for the case of in-plane shear loading. The problem also tests the effect of mesh distortion. The membrane is fully fixed at $X=0$ and loaded by a distributed force at $X=L_{1}$ (see Figure 10). The $Z$-component of the displacement is evaluated at the test point A. The analytical solution for this example is not available, therefore the reference solution is obtained numerically using the DKMQ242+ element on the refined mesh $192\!\times\!128$. The numerical results are summarized in Table 3. We observe a massive reduction of the relative error in case of the DKMQ241+ and DKMQ242+ shell elements in comparison to the MITC4+ and DKMQ24 shell elements.

5.2 Hemispherical shell

The hemispherical shell problem, consisting of a hemisphere loaded by four opposite single forces, is depicted in Figure 12. Due to symmetries, only one quarter of this hemisphere needs to be considered. The following symmetry conditions are applied: $u_{Y}=\varphi_{X}=\varphi_{Z}=0$ on the side AC and $u_{X}=\varphi_{Y}=\varphi_{Z}=0$ on the side BC. Additionally, the following boundary condition is applied at point C: $u_{Z}=0$. The analytical solution to this problem is given by [29] as $u_{X}(A)=-u_{Y}(B)=92.400\,\text{mm}$.

Note that this example is sensitive to the mesh quality and it is also the only benchmark in our set for which the mesh is not exactly specified by its depiction in Figure 12. In our case we have used the Laplacian smoothing of internal nodes [33] until we have converged to a unique mesh for the calculation.

We report that our implementation of the DKMQ24 shell element gives significantly higher error then the one reported in [18]. A similar error was obtained for the reduced penalty constant $c=10^{-7}$. In comparison to the DKMQ24 element from [18] and MITC4+ element, our elements DKMQ241+ and DKMQ242+ provide higher errors on coarser meshes and lower errors on finer meshes.

5.3 Scordelis-Lo roof

The Scordelis-Lo roof is a cylindrical surface with radius $R$, length $L$ and thickness $h$, reinforced at positions $Y=0$ and $Y=L$ by two rigid diaphragms (see Figure 14). Because of symmetry only one quarter is calculated. The shell is loaded by self-weight, where gravitational acceleration is considered to be $g=10\,\text{m}\,\text{s}^{-2}$, which yields the total force acting on the roof $f_{Z}=-h\rho gA=-12500\pi\,\text{N}$. The following boundary and symmetric conditions are applied: $u_{X}=u_{Z}=\varphi_{Y}=0$ on the side AD, $u_{X}=\varphi_{Y}=\varphi_{Z}=0$ on the side CD and $u_{Y}=\varphi_{X}=\varphi_{Z}=0$ on the side CB. The analytical solution is given in [34]. We note that two versions of this benchmark, having the same geometry setting differing only in the parameter setting, are given in the literature. The one version is taken from [34] (described in Fig. 14 and used in Tables 7-10 for comparison with the DKMQ24 [18] element), the second version taken from [29] (used in Table 11 for comparison with the MITC4+ and US-ATFHS8 elements).

With this benchmark problem we test not only nodal displacements, but also values of the nodal stress-resultants. We point out that nodal stress-resultants are bilinearly extrapolated from the Gauss quadrature points to nodes.

If compared to the DKMQ24 element from [18], our modified elements DKMQ241+ and DKMQ242+ exhibit lower errors at all test cases. If compared to the MITC4+ and US-ATFHS8 elements, our modified elements DKMQ241+ and DKMQ242+ exhibit lower errors on both coarser meshes, but higher errors on the finest mesh $16\times 16$.

5.4 Pinched cylinder shell

The pinched cylinder, depicted in Figure 20, is loaded by two opposite pinching forces $F$. The cylinder has a rigid diaphragm at both ends. Only one eighth of the model is calculated due to symmetries (to which the force $\frac{F}{4}$ is applied only). The boundary conditions are as follows: $u_{X}=u_{Z}=\varphi_{Y}=0$ on the side AD, on the remaining sides symmetric conditions are prescribed: $u_{Z}=\varphi_{X}=\varphi_{Y}=0$ on side AB, $u_{Y}=\varphi_{X}=\varphi_{Z}=0$ on the side BC and $u_{X}=\varphi_{Y}=\varphi_{Z}=0$ on side CD. Two different shell thicknesses are considered in our calculation, $h=0.03\,\text{m}$ and $h=0.3\,\text{m}$. We use Young’s modulus $E=3\,\text{MPa}$ instead of the value $30000\,\text{MPa}$ used in [18]. The uniform and distorted meshes are analysed. For the distorted mesh a bias factor, i.e. the ratio of the largest and smallest element edge, of 10 is used. The analytical solution for the case $h=0.03\,\text{m}$ is based on the Kirchhoff theory proposed by Lindberg, et al. [35], which follows the solution of Flügge [36], and equals to $u_{Z}(C)=-1.825\,\text{mm}$. The analytical solution for the case $h=0.3\,\text{m}$, taken from [37], is based on the Mindlin theory and equals to $u_{Z}(C)=-0.01261\,\text{mm}$. However, as reported in [18] and according to our numerical results, the numerical results converge to slightly larger values. Therefore, numerical results of the DKMQ242+ shell element on mesh $128\!\times\!128$ are taken as reference values. The numerical results presented in Tables 9-12 reveal slow convergence of all the considered elements, which is in concordance with other authors [30, 38] and is a consequence of the single point loading.

When compared to the DKMQ24 [18] element, our modified elements DKMQ241+ and DKMQ242+ have almost identical results (for $h=0.03\,\text{m}$ at some cases slightly better, at some slightly worse, for $h=0.3\,\text{m}$ slightly better at all test cases). When compared to the MITC4+ element, where we only have data for the case $h=0.03\,\text{m}$, the modified elements DKMQ241+ and DKMQ242+ exhibit lower errors in case of the uniform mesh, and slightly higher errors in case of the distorted mesh.

5.5 Hyperbolic paraboloid shell

The hyperbolic paraboloid shell (hypar shell) problem is depicted in Figure 25. The shell middle surface is described by $Z=\frac{b}{a^{2}}XY$. The shell is loaded by self-weight, where the gravitational acceleration is considered to be $g=10\,\text{m}\,\text{s}^{-2}$. The shell is fully clamped on all sides ABCD, i.e. $u_{X}=u_{Y}=u_{Z}=\varphi_{X}=\varphi_{Y}=\varphi_{Z}=0$. The analytical solution is not available ([39], p.245), therefore we use the numerical result on the fine mesh as reference. Our modified elements DKMQ241+ and DKMQ242+ provide almost identical results at this benchmark when compared to the shell element DKMQ24 [18].

5.6 Twisted beam

A shell strip of length $L$, width $b$ and thickness $h$ is twisted in the undeformed configuration by angle $\frac{\pi}{2}$ between points $O$ and $A$ (see Figure 27). The strip is clamped at $X\!=\!0$ and loaded by a single force $[0,F_{Y},F_{Z}]$ at point $A$. The shell middle surface is described by: $[X(\varphi,s),Y(\varphi,s),Z(\varphi,s)]=\left[\frac{2\varphi L}{\pi},~{}s\sin\varphi,~{}s\cos\varphi\right],$ $0\leq\varphi\leq\frac{\pi}{2},-\frac{b}{2}\leq s\leq\frac{b}{2}$. The example was introduced by MacNeal and Harder [29] and represents an important test example at which many flat quadrilateral shell elements fail to converge. This is in accordance with our numerical experience, which shows that the DKMQ24/DKMQ241+/DKMQ242+ elements fail to converge if the coupling term, including derivatives of the element normal with respect to reference coordinates in Eq. (193, 206), which is identicaly zero for flat elements, is missing. The analytical solution based on the beam theory (which ignores the deformation in the transverse direction) is taken from [40]. Four different settings of loading forces and shell thicknesses were tested, see Tables 17-20.

For the case $h=0.0032\,\text{m}$, when compared to the DKMQ24 [18] element, our modified elements DKMQ241+ and DKMQ242+ provide smaller errors at all mesh sizes except at the coarsest one. When comparing the same case $h=0.0032\,\text{m}$ with the MITC4+ element, our modified elements DKMQ241+ and DKMQ242+ provide smaller errors on the finer meshes. The smallest error at mesh size $4\times 24$ is given by the US-ATFHS8 element. For the case $h=0.0032\,\text{m}$ our modified elements DKMQ241+ and DKMQ242+ are superior to all other elements in all test cases.

5.7 Raasch’s hook

Raasch’s hook problem was first reported by Ingo Raasch of the BMW corporation in Germany (see Figure 32). It is a challenging problem because of the inherent coupling between three modes of deformation: bending, extension, and twist. Many shell elements with drilling dofs exhibit divergence of the tip deflection at this benchmark ([31]). Except of the coarsest mesh $17\times 3$, which has uniformly distributed elements, all other meshes we used have, in accordance with [32], the same number of elements before and after the inflection point, see Figure 32.

The results are shown in Table 21. We note that a sufficient size of the penalty constant is needed in order to obtain convergence at this benchmark.

Our implementation of the DKMQ24 element exhibits large errors at this benchmark and, moreover, diverges to infinity when the mesh gets finer. An increase of the penalty constant $c$ does not yield expected convergence, moreover, an increase of the penalty constant negatively influences the error at other benchmarks. For further discussion on this topic see Section 4. The DKMQ242+ element provides lower errors when compared to the DKMQ241+ element. Moreover, the DKMQ242+ element provides a fast convergence to a fixed value, which is however slightly higher then the reference value (by $0.3$%). The DKMQ242+ element and the MITC4+ element provide similar errors, however, we note that the MITC4+ element is tested with a slightly different setting ($\nu=0.3$ instead of $\nu=0.35$) and, therefore, converging to a different target value ($4.82482\,\text{m}$ instead of $5.027\,\text{m}$).

5.8 Partly clamped hyperbolic paraboloid shell

The partly clamped hyperbolic paraboloid shell benchmark, taken from [32], whose geometry is described by equation

[TABLE]

is loaded by self-weight (see Figure 34). On side CD, symmetry conditions are prescribed ($u_{Y}=\varphi_{X}=\varphi_{Z}=0$), the side DA is fully clamped ($u_{X}=u_{Y}=u_{Z}=\varphi_{X}=\varphi_{Y}=\varphi_{Z}=0$) and other sides are free. Two different shell thicknesses are considered: $h/L=0.001$ and $h/L=0.0001$. The corresponding load magnitudes are scaled by using different densities $\rho=360\,\text{kg}\,\text{m}^{-3}$ and $\rho=3.6\,\text{kg}\,\text{m}^{-3}$, respectively. The reference solution is taken from [32].

Let us first comment the results on uniform meshes. For the case $h=0.001\,\text{m}$ the modified elements DKMQ241+ and DKMQ242+ are superior to all other elements, except of the coarsest mesh $4\times 8$. For the case $h=0.0001\,\text{m}$ the modified elements DKMQ241+ and DKMQ242+ provide smaller errors when compared to the MITC4+ element only at finer meshes ($32\times 64$ and $64\times 128$). Let us now comment the results on distorted meshes. For the case $h=0.001\,\text{m}$ the modified elements DKMQ241+ and DKMQ242+ when compared to the MITC4+ element, are superior at finer meshes ($32\times 64$ and $64\times 128$), comparable at the middle-sized mesh $16\times 32$ and worse on coarser meshes. For the case $h=0.0001\,\text{m}$ the MITC4+ elements shows lowest errors at all mesh sizes, which demonstrated its good ability to resist to a mesh distorsion. In all the tested cases we see that the modified elements DKMQ241+ and DKMQ242+ provide significantly lower errors when compared to the DKMQ24 element, which demonstrated better behaviour of both elements at membrane-dominated problems.

5.9 Square clamped plate under uniform loading with distorted mesh

A fully clamped pressure loaded square plate with distorted mesh is considered (see Figure (40)). Only one quarter is computed due to symmetries. Boundary and symmetry conditions are as follows: $u_{X}=u_{Y}=u_{Z}=\varphi_{X}=\varphi_{Y}=\varphi_{Z}=0$ on sides AB and AD, $u_{X}=\varphi_{Y}=\varphi_{Z}=0$ on side BC and $u_{Y}=\varphi_{X}=\varphi_{Z}=0$ on side CD. This benchmark problem tests sensitivity to mesh distortion. The reference solution is taken from [42]222The value $62.83$ given in [18], taken from [42], is not the deformation $u_{Z}$ as stated, but $\frac{Gu_{Z}}{hp}$.. This benchmark is taken as an example of pure bending, where our modifications of the membrane part of the DKMQ24 element does not play any role. Despite of that, our modified elements DKMQ241+ and DKMQ242+ provide significantly lower errors in all test cases, which is due to the nodal moment corrections.

6 Conclusions

The new quadrangle shell element DKMQ242+ with four nodes and six degrees of freedom per node, based on the DKMQ24 quadrangle shell element introduced by Katili et al. [18], with improved membrane behaviour was introduced. The element DKMQ242+ uses quadratic enrichment of the in-plane approximations with the help of drilling rotations, a selective reduced integration of shear terms, a proportional scaling of the penalty coefficient to the shell thickness $h$ and nodal moment corrections. Another variant of the element formulation called DKMQ241+ uses a constant penalty coefficient without any scaling and otherwise is identical with the DKMQ242+ variant. According to our numerical comparisons, the DKMQ242+ element variant is superior to the DKMQ241+ element variant.

The new quadrangle shell element DKMQ242+, when compared to the DKMQ24 element by Katili [18], exhibits lower errors on majority of the tested benchmark problems, in particular on pure membrane problems and on shell problems where in-plane loading is dominant. Moreover, the presented element provides convergent behaviour at Raasch’s hook benchmark.

The DKMQ242+ shell element, when compared to the MITC4+ and US-ATFHS8 elements, exhibits in many cases a lower level of error. The error is in general higher on coarser meshes and smaller on finer meshes. The resistance to distored meshes of the DKMQ242+ shell element is not as high as of the MITC4+ and US-ATFHS8 shell elements.

The DKMQ242+ shell element passes all our convergence tests, without presence of neither shear locking nor membrane locking and converges in all considered cases to the reference solution. Moreover, all our modifications, except of the static condensation, do not introduce any additional computational costs when compared to the DKMQ24 shell element. To conclude, the DKMQ242+ shell element could be considered as a replacement of the DKMQ24 shell element.

7 Acknowledgements

The authors gratefully acknowledge the support of the FWF Austrian Science Fund through the START project Y 1093 as well as the support through the doctoral college Computational Design at TU Wien.

8 Appendix A – Static condensation

We consider the following system of equations defined on a single element $e$

[TABLE]

where

[TABLE]

If we eliminate the unknown $\bm{q}_{\text{n}}$, we get the eliminated system of equations of the form

[TABLE]

where

[TABLE]

From the programming point of view the following equivalent form can be advantageous

[TABLE]

because the condensation is implemented only on the strain matrix level.