Benchmark: MacNeal-Harder standard problems for the shell elements#
Date: 22 of October 2026
Author: Saullo G. P. Castro
MacNeal and Harder [1] proposed a set of small problems that excite, with coarse and distorted meshes, the failure modes of shell and solid elements: spurious mechanisms and locking. This notebook runs the problems of that set that a flat shell element can model with the four shell elements of pyfe3d, Quad4, Quad4R, Tria3R and Tria3DSG, every element with its default attributes, and compares the normalised results with the theoretical values of the paper and with the results
of the MSC/NASTRAN QUAD4 element tabulated there.
Problem of [1] |
Figure, tables |
Here |
|---|---|---|
Patch test for plates |
Fig. 2, Tables 2 and 7 |
membrane and bending |
Straight cantilever beam |
Fig. 4, Tables 3 and 8 |
regular, trapezoidal and parallelogram meshes |
Curved beam |
Fig. 5, Tables 3 and 9 |
in-plane and out-of-plane load |
Twisted beam |
Fig. 6, Tables 3 and 10 |
in-plane and out-of-plane load |
Rectangular plate |
Fig. 7, Tables 4, 11 and 12 |
the two cases of Tables 11 and 12 |
Scordelis-Lo roof, spherical shell |
Figs. 8 and 9, Tables 5, 13 and 14 |
in the notebook of the shell obstacle course |
Patch test for solids, thick-walled cylinder |
Figs. 3 and 10 |
not modelled: solid and plane strain problems |
The triangular meshes are obtained by cutting every quadrilateral of the mesh of [1] along its diagonal 1-3, as the paper suggests for triangles (p. 7). Unless said otherwise, a displacement is normalised by the theoretical value of [1], and the letter grades follow the rule of [1], p. 14: A for an error of at most 2%, B up to 10%, C up to 20%, D up to 50% and F above, the absolute errors being averaged when several cases contribute to one grade.
Common code#
stiffness() assembles the linear stiffness matrix KC0 of a mesh of quadrilaterals for one of the four elements, solve() solves the static problem with prescribed degrees-of-freedom, tip_forces() distributes a unit end load as in Fig. 11(a) of [1], half of the load of each element edge going to each of its two nodes.
[1]:
import numpy as np
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import spsolve
from pyfe3d.shellprop_utils import isotropic_plate
from pyfe3d import (Quad4, Quad4Data, Quad4Probe, Quad4R, Quad4RData,
Quad4RProbe, Tria3R, Tria3RData, Tria3RProbe, Tria3DSG,
Tria3DSGData, Tria3DSGProbe, INT, DOUBLE, DOF)
ELEMENTS = {
'Quad4': (Quad4, Quad4Data, Quad4Probe),
'Quad4R': (Quad4R, Quad4RData, Quad4RProbe),
'Tria3R': (Tria3R, Tria3RData, Tria3RProbe),
'Tria3DSG': (Tria3DSG, Tria3DSGData, Tria3DSGProbe),
}
def split(quads):
# two triangles per quadrilateral, cut along the diagonal 1-3
return [t for (n1, n2, n3, n4) in quads for t in ((n1, n2, n3), (n1, n3, n4))]
def stiffness(ncoords, quads, name, prop):
Elem, Data, Probe = ELEMENTS[name]
conn = split(quads) if name.startswith('Tria') else quads
data, probe = Data(), Probe()
ncf = np.ascontiguousarray(ncoords, dtype=DOUBLE).flatten()
N = DOF*ncoords.shape[0]
size = data.KC0_SPARSE_SIZE*len(conn)
KC0r = np.zeros(size, dtype=INT)
KC0c = np.zeros(size, dtype=INT)
KC0v = np.zeros(size, dtype=DOUBLE)
areas = np.zeros(len(conn))
for k, nodes in enumerate(conn):
e = Elem(probe)
for i, n in enumerate(nodes):
setattr(e, 'n%d' % (i + 1), n)
setattr(e, 'c%d' % (i + 1), DOF*n)
e.init_k_KC0 = k*data.KC0_SPARSE_SIZE
e.update_rotation_matrix(ncf)
e.update_probe_xe(ncf)
e.update_area()
e.update_KC0(KC0r, KC0c, KC0v, prop)
areas[k] = e.area
K = coo_matrix((KC0v, (KC0r, KC0c)), shape=(N, N)).tocsc()
return K, conn, areas
def solve(K, bk, f, uk=None):
# static solution with the degrees-of-freedom of the mask bk prescribed to uk
u = np.zeros(K.shape[0])
if uk is not None:
u[bk] = uk[bk]
bu = ~bk
u[bu] = spsolve(K[bu][:, bu], f[bu] - K[bu][:, bk] @ u[bk])
return u
def area_load(ncoords, conn, areas, q):
# nodal forces of a load per unit area q, lumped equally to the nodes
f = np.zeros(DOF*ncoords.shape[0])
for nodes, A in zip(conn, areas):
for n in nodes:
f[DOF*n:DOF*n + 3] += np.asarray(q)*A/len(nodes)
return f
def grid(nx, ny):
# node ids of a structured grid and its quadrilaterals
ids = np.arange((nx + 1)*(ny + 1)).reshape(nx + 1, ny + 1)
quads = [(ids[i, j], ids[i + 1, j], ids[i + 1, j + 1], ids[i, j + 1])
for i in range(nx) for j in range(ny)]
return ids, quads
def tip_forces(N, tip, comp, total=1.):
f = np.zeros(N)
for a, b in zip(tip[:-1], tip[1:]):
f[DOF*a + comp] += total/(len(tip) - 1)/2
f[DOF*b + comp] += total/(len(tip) - 1)/2
return f
def grade(normalised):
# letter grade of [1], p. 14, from the mean absolute error
err = 100*np.mean(np.abs(np.atleast_1d(normalised) - 1.))
for letter, limit in (('A', 2.), ('B', 10.), ('C', 20.), ('D', 50.)):
if err <= limit:
return letter
return 'F'
Patch tests for plates#
The patch of Fig. 2 of [1] is a rectangle of 0.24 by 0.12 divided into five quadrilaterals by four interior nodes at (0.04, 0.02), (0.18, 0.03), (0.16, 0.08) and (0.08, 0.08), with \(t = 0.001\), \(E = 10^6\) and \(\nu = 0.25\). The displacements of the exact field of Table 2 of [1] are prescribed at the four corner nodes and the interior nodes are free:
membrane: \(u = 10^{-3}(x + y/2)\), \(v = 10^{-3}(y + x/2)\), a state of constant strain \(\epsilon_x = \epsilon_y = \gamma_{xy} = 10^{-3}\);
bending: \(w = 10^{-3}(x^2 + xy + y^2)/2\), with the rotations of zero transverse shear, which in the convention of
pyfe3dare \(r_x = \partial w/\partial y\) and \(r_y = -\partial w/\partial x\), a state of constant curvature.
The out-of-plane degrees-of-freedom are fixed in the membrane test, and \(u\), \(v\) and \(r_z\) in the bending test. The table of [1] prescribes \(u\) and \(v\) only, for elements without drilling rotations. Quad4 and Tria3R enrich their in-plane displacements with the drilling rotations [2], so that the displacement of a boundary edge depends on \(r_z\) at its nodes, and the drilling rotation of the boundary nodes is prescribed as well, to the in-plane rotation
of the exact field, \((v_{,x} - u_{,y})/2 = 0\). Left free, it lets the boundary edges deform beyond the prescribed field.
[1] records the maximum error of the stresses (Table 7). Stresses are recovered differently by each element of pyfe3d, so the measures here are the maximum error of the displacements of the interior nodes, relative to the largest value of the exact field, and the error of the strain energy \(U = \frac{1}{2} \boldsymbol{u}^\top \boldsymbol{K} \boldsymbol{u}\) against the exact \(\frac{1}{2} \boldsymbol{\epsilon}^\top \boldsymbol{A} \boldsymbol{\epsilon}\) (membrane) or
\(\frac{1}{2} \boldsymbol{\kappa}^\top \boldsymbol{D} \boldsymbol{\kappa}\) (bending) times the area. Both are zero, to round-off, when the element passes.
[2]:
PATCH_XY = np.array([[0., 0.], [0.24, 0.], [0.24, 0.12], [0., 0.12],
[0.04, 0.02], [0.18, 0.03], [0.16, 0.08], [0.08, 0.08]])
PATCH_QUADS = [(4, 5, 6, 7), (0, 1, 5, 4), (1, 2, 6, 5), (2, 3, 7, 6), (3, 0, 4, 7)]
def patch_test(name, form):
ncoords = np.column_stack((PATCH_XY, np.zeros(8)))
prop = isotropic_plate(thickness=0.001, E=1.e6, nu=0.25)
K, conn, areas = stiffness(ncoords, PATCH_QUADS, name, prop)
N = K.shape[0]
x, y = PATCH_XY.T
boundary = np.arange(8) < 4
uex = np.zeros(N)
bk = np.zeros(N, dtype=bool)
if form == 'membrane':
uex[0::DOF] = 1e-3*(x + y/2)
uex[1::DOF] = 1e-3*(y + x/2)
checked, fixed, prescribed = (0, 1), (2, 3, 4), (0, 1, 5)
C = np.array([[prop.A11, prop.A12, prop.A16],
[prop.A12, prop.A22, prop.A26],
[prop.A16, prop.A26, prop.A66]])
else:
uex[2::DOF] = 1e-3*(x**2 + x*y + y**2)/2
uex[3::DOF] = 1e-3*(y + x/2)
uex[4::DOF] = -1e-3*(x + y/2)
checked, fixed, prescribed = (2, 3, 4), (0, 1, 5), (2, 3, 4)
C = np.array([[prop.D11, prop.D12, prop.D16],
[prop.D12, prop.D22, prop.D26],
[prop.D16, prop.D26, prop.D66]])
for d in fixed:
bk[d::DOF] = True
for d in prescribed:
bk[d::DOF] |= boundary
u = solve(K, bk, np.zeros(N), uex)
err_u = max(np.abs(u[d::DOF][~boundary] - uex[d::DOF][~boundary]).max()
/np.abs(uex[d::DOF]).max() for d in checked)
# the three strains, or the three curvatures, are all equal to 1e-3 in magnitude
e = np.full(3, 1e-3)
U_exact = 0.5*e @ C @ e*0.24*0.12
err_U = 0.5*u @ (K @ u)/U_exact - 1.
return err_u, err_U
patch = {}
print('element test interior displ. error energy error check')
for name in ELEMENTS:
for form in ('membrane', 'bending'):
err_u, err_U = patch[name, form] = patch_test(name, form)
ok = err_u < 1e-8 and abs(err_U) < 1e-8
print('%-10s %-10s %20.2e %14.2e %s' % (name, form, err_u, err_U,
'PASS' if ok else 'FAIL'))
element test interior displ. error energy error check
Quad4 membrane 2.82e-16 -1.11e-16 PASS
Quad4 bending 7.39e-01 5.80e+00 FAIL
Quad4R membrane 4.52e-16 -1.11e-16 PASS
Quad4R bending 7.24e-01 1.34e+01 FAIL
Tria3R membrane 6.78e-16 -4.44e-16 PASS
Tria3R bending 8.41e-02 2.43e-01 FAIL
Tria3DSG membrane 1.08e-15 -1.11e-16 PASS
Tria3DSG bending 3.82e-14 -2.38e-11 PASS
All four elements pass the membrane patch test. Only Tria3DSG passes the bending patch test. [1] reports 0.0% for the QUAD4 of MSC/NASTRAN in both tests (Table 7).
The failure of Quad4 and Quad4R in bending is in the transverse shear. The exact transverse displacement is quadratic and both quadrilaterals interpolate it bilinearly, with the transverse shear sampled at the centroid; on these distorted elements the gradient of the bilinear interpolant at the centroid differs from the prescribed rotation there, on every element of the patch, which leaves a transverse shear strain and, with \(\ell/t \approx 100\), a large parasitic shear energy.
Tria3R takes its transverse shear at the centroid as well, softened by its alpha_shear_locking factor. The discrete shear gap of Tria3DSG makes the shear strain of a state of constant curvature vanish exactly, which is the property this test checks.
Straight cantilever beam#
Fig. 4 of [1]: length 6.0, width 0.2, thickness 0.1, \(E = 10^7\), \(\nu = 0.3\), a 6 by 1 mesh clamped at one end, with regular, trapezoidal and parallelogram elements, the interior node lines of the distorted meshes being slanted at 45 degrees. The unit loads are applied at the free end, the twist as a couple of transverse forces \(\pm 1/0.2\) on the two tip nodes (Fig. 11(b) of [1]). The extension and the shear results are the mean tip displacement in the direction of the load, the twist result is the tip rotation \((w_2 - w_1)/0.2\).
[3]:
STRAIGHT_REF = {'extension': 3.0e-5, 'in-plane shear': 0.1081,
'out-of-plane shear': 0.4321, 'twist': 0.03208} # Table 3 of [1]
MH_QUAD4_STRAIGHT = { # Table 8 of [1]
'regular': {'extension': 0.995, 'in-plane shear': 0.904, 'out-of-plane shear': 0.986, 'twist': 0.941},
'trapezoidal': {'extension': 0.996, 'in-plane shear': 0.071, 'out-of-plane shear': 0.968, 'twist': 0.951},
'parallelogram': {'extension': 0.996, 'in-plane shear': 0.080, 'out-of-plane shear': 0.977, 'twist': 0.945},
}
def straight_beam(name, shape, load):
L, w, nx, ny = 6., 0.2, 6, 1
ids, quads = grid(nx, ny)
ys = np.linspace(0., w, ny + 1)
ncoords = np.zeros(((nx + 1)*(ny + 1), 3))
for i in range(nx + 1):
d = 0.
if 0 < i < nx and shape == 'parallelogram':
d = 0.1
elif 0 < i < nx and shape == 'trapezoidal':
d = 0.1 if i % 2 else -0.1
ncoords[ids[i, :], 0] = i*L/nx + d*(2*ys/w - 1)
ncoords[ids[i, :], 1] = ys
prop = isotropic_plate(thickness=0.1, E=1.e7, nu=0.3)
K, conn, areas = stiffness(ncoords, quads, name, prop)
N = K.shape[0]
bk = np.zeros(N, dtype=bool)
for d in range(DOF):
bk[d::DOF] |= np.isclose(ncoords[:, 0], 0.)
tip = ids[-1, :]
if load == 'twist':
f = np.zeros(N)
f[DOF*tip[0] + 2] = -1/w
f[DOF*tip[-1] + 2] = +1/w
u = solve(K, bk, f)
return (u[DOF*tip[-1] + 2] - u[DOF*tip[0] + 2])/w
comp = {'extension': 0, 'in-plane shear': 1, 'out-of-plane shear': 2}[load]
u = solve(K, bk, tip_forces(N, tip, comp))
return u[DOF*tip + comp].mean()
straight = {}
print('normalised tip displacement')
print('%-14s %-19s' % ('mesh', 'load') + ''.join('%10s' % n for n in ELEMENTS) + ' MH QUAD4')
for shape in ('regular', 'trapezoidal', 'parallelogram'):
for load, ref in STRAIGHT_REF.items():
for name in ELEMENTS:
straight[name, shape, load] = straight_beam(name, shape, load)/ref
print('%-14s %-19s' % (shape, load)
+ ''.join('%10.3f' % straight[name, shape, load] for name in ELEMENTS)
+ '%10.3f' % MH_QUAD4_STRAIGHT[shape][load])
normalised tip displacement
mesh load Quad4 Quad4R Tria3R Tria3DSG MH QUAD4
regular extension 0.999 0.998 0.992 0.992 0.995
regular in-plane shear 0.904 0.119 0.145 0.027 0.904
regular out-of-plane shear 0.980 0.973 0.876 0.971 0.986
regular twist 377.368 4.780 1.669 0.675 0.941
trapezoidal extension 0.998 0.999 0.992 0.991 0.996
trapezoidal in-plane shear 0.761 0.033 0.149 0.013 0.071
trapezoidal out-of-plane shear 0.981 0.233 0.589 0.959 0.968
trapezoidal twist 109.715 4.907 3.244 0.658 0.951
parallelogram extension 0.998 0.999 0.992 0.991 0.996
parallelogram in-plane shear 0.838 0.068 0.102 0.010 0.080
parallelogram out-of-plane shear 0.982 0.772 0.831 0.958 0.977
parallelogram twist 137.968 5.258 2.709 0.638 0.945
Extension is reproduced by every element on every mesh. The other cases show the following failures:
Twist of ``Quad4``. The one-point transverse shear of
Quad4, from Hughes et al. [3], leaves the hourglass mode \(w = (x - x_c)(y - y_c)\) without stiffness, and a strip one element wide loaded in torsion deforms in that mode, so the twist is two orders of magnitude too large.Quad4Rcontrols that mode with its hourglass stiffness, and the result depends on the amount of it, here about five times too flexible. Both are documented in the Notes of the two modules, and the teststest_quad4_spurious_shear_mode.pyandtest_drilling_benchmarks.pyshow that two elements across the width recover the twist withQuad4.In-plane shear, i.e. in-plane bending.
Quad4, with the drilling enrichment of Allman [2], reaches 90% of the reference on the regular mesh, like the QUAD4 of [1], and keeps 76% to 84% on the distorted meshes, where the QUAD4 of [1] locks (7% to 8%).Quad4Rlocks on all three meshes: its one-point membrane gives no strain to the in-plane bending pattern of the bilinear element, which is resisted only by its hourglass stabilisation and by its drilling term, and on this one-element-deep strip these make it far too stiff. The triangles lock too:Tria3DSGhas the constant-strain membrane of the linear triangle, andTria3R, which carries the Allman enrichment, reaches 10% to 15%.Out-of-plane shear on distorted meshes.
Quad4Rdrops to 23% on the trapezoidal mesh andTria3Rto 59%.Twist of the triangles.
Tria3Ris 1.7 to 3.2 times too flexible andTria3DSG32% to 36% too stiff.
Curved beam#
Fig. 5 of [1]: a 90 degree arc of inner radius 4.12 and outer radius 4.32, thickness 0.1, \(E = 10^7\), \(\nu = 0.25\), a 6 by 1 mesh clamped at one end and loaded by unit forces at the tip, in the plane of the arc (radially at the tip) or normal to it.
[4]:
CURVED_REF = {'in-plane': 0.08734, 'out-of-plane': 0.5022} # Table 3 of [1]
MH_QUAD4_CURVED = {'in-plane': 0.833, 'out-of-plane': 0.951} # Table 9 of [1]
def curved_beam(name, load):
nx, ny = 6, 1
th = np.linspace(0., np.pi/2, nx + 1)
r = np.linspace(4.12, 4.32, ny + 1)
T, Rr = np.meshgrid(th, r, indexing='ij')
ncoords = np.column_stack(((Rr*np.cos(T)).ravel(), (Rr*np.sin(T)).ravel(),
np.zeros(T.size)))
ids, quads = grid(nx, ny)
prop = isotropic_plate(thickness=0.1, E=1.e7, nu=0.25)
K, conn, areas = stiffness(ncoords, quads, name, prop)
N = K.shape[0]
bk = np.zeros(N, dtype=bool)
for d in range(DOF):
bk[d::DOF] |= np.isclose(ncoords[:, 1], 0.)
tip = ids[-1, :]
comp = 1 if load == 'in-plane' else 2
u = solve(K, bk, tip_forces(N, tip, comp))
return u[DOF*tip + comp].mean()
curved = {}
print('normalised tip displacement')
print('%-14s' % 'load' + ''.join('%10s' % n for n in ELEMENTS) + ' MH QUAD4')
for load, ref in CURVED_REF.items():
for name in ELEMENTS:
curved[name, load] = curved_beam(name, load)/ref
print('%-14s' % load + ''.join('%10.3f' % curved[name, load] for name in ELEMENTS)
+ '%10.3f' % MH_QUAD4_CURVED[load])
normalised tip displacement
load Quad4 Quad4R Tria3R Tria3DSG MH QUAD4
in-plane 0.908 0.103 0.119 0.022 0.833
out-of-plane 26.828 1.587 1.100 0.760 0.951
The in-plane case repeats the in-plane bending result of the straight beam: Quad4 reaches 91%, better than the 83% of the QUAD4 of [1], and the other three elements lock. The out-of-plane load twists the strip, and Quad4 is again far too flexible through the hourglass mode of its transverse shear.
Twisted beam#
Fig. 6 of [1]: length 12.0, width 1.1, thickness 0.32, twisted by 90 degrees from root to tip, \(E = 29 \times 10^6\), \(\nu = 0.22\), a 12 by 2 mesh clamped at the root, with unit forces at the tip. The reference surface is a helicoid, so every quadrilateral is warped, by 7.5 degrees according to [1] (p. 11). At the tip the width lies along \(z\): the in-plane load is along \(z\) and the out-of-plane load along \(y\).
The quadrilaterals of pyfe3d take their normal from the cross product of the two diagonals and work with the coordinates of the nodes in that frame (update_rotation_matrix() and update_probe_xe()), with no correction for warping; a triangle is flat by construction. The mesh is refined beyond the one of [1] to show the effect.
[5]:
TWISTED_REF = {'in-plane': 0.005424, 'out-of-plane': 0.001754} # Table 3 of [1]
MH_QUAD4_TWISTED = {'in-plane': 0.993, 'out-of-plane': 0.985} # Table 10 of [1]
def twisted_beam(name, load, nx=12, ny=2):
L, w = 12., 1.1
X, S = np.meshgrid(np.linspace(0., L, nx + 1), np.linspace(-w/2, w/2, ny + 1),
indexing='ij')
beta = np.pi/2*X/L
ncoords = np.column_stack((X.ravel(), (S*np.cos(beta)).ravel(),
(S*np.sin(beta)).ravel()))
ids, quads = grid(nx, ny)
prop = isotropic_plate(thickness=0.32, E=29.e6, nu=0.22)
K, conn, areas = stiffness(ncoords, quads, name, prop)
N = K.shape[0]
bk = np.zeros(N, dtype=bool)
for d in range(DOF):
bk[d::DOF] |= np.isclose(ncoords[:, 0], 0.)
tip = ids[-1, :]
comp = 2 if load == 'in-plane' else 1
u = solve(K, bk, tip_forces(N, tip, comp))
return u[DOF*tip + comp].mean()
twisted = {}
print('normalised tip displacement')
print('%-9s %-14s' % ('mesh', 'load') + ''.join('%10s' % n for n in ELEMENTS) + ' MH QUAD4')
for nx, ny in ((12, 2), (48, 4), (192, 8), (768, 16)):
for load, ref in TWISTED_REF.items():
for name in ELEMENTS:
twisted[name, nx, load] = twisted_beam(name, load, nx, ny)/ref
print('%-9s %-14s' % ('%dx%d' % (nx, ny), load)
+ ''.join('%10.3f' % twisted[name, nx, load] for name in ELEMENTS)
+ ('%10.3f' % MH_QUAD4_TWISTED[load] if nx == 12 else ''))
normalised tip displacement
mesh load Quad4 Quad4R Tria3R Tria3DSG MH QUAD4
12x2 in-plane 0.252 0.234 0.910 0.969 0.993
12x2 out-of-plane 0.332 0.325 0.884 0.817 0.985
48x4 in-plane 0.659 0.691 0.967 0.994
48x4 out-of-plane 0.711 0.743 0.965 0.954
192x8 in-plane 0.899 0.946 0.994 0.998
192x8 out-of-plane 0.917 0.957 0.994 0.990
768x16 in-plane 0.973 0.992 0.998 0.999
768x16 out-of-plane 0.978 0.994 0.998 0.997
On the mesh of [1] the two quadrilaterals are three to four times too stiff, while the QUAD4 of [1] is within 2%. The triangles, which are flat by construction, are within 20%. The error of the quadrilaterals decreases as the warping of each element decreases with the refinement, and all four elements approach the reference of [1] on the finest mesh. MacNeal [4], one of the two references given for the QUAD4 in Table 6 of [1], modifies the matrices of the plane element for warped quadrilaterals;
nothing equivalent is implemented in the quadrilaterals of pyfe3d.
Rectangular plate#
Fig. 7 of [1]: \(a = 2.0\), \(b = 2.0\) or \(10.0\), thickness 0.0001, \(E = 1.7472 \times 10^7\), \(\nu = 0.3\), an \(N \times N\) mesh on one quarter of the plate with symmetry conditions on the two inner edges. The two cases tabulated by [1] for the QUAD4 are computed: simply supported edges with a uniform pressure \(q = 10^{-4}\) (Table 11) and clamped edges with a central load \(P = 4 \times 10^{-4}\) (Table 12), the quarter model carrying \(P/4\). Simply supported means \(w = 0\) on the edge, with free rotations; clamped means \(w = r_x = r_y = 0\). The in-plane degrees-of-freedom \(u\), \(v\) and \(r_z\) are fixed, the plate being loaded in bending only. The pressure is lumped to the nodes, one quarter (one third for triangles) of the load of each element going to each node. The result is the deflection at the centre normalised by Table 4 of [1].
[6]:
PLATE_REF = {('simple', 1.): 4.062, ('simple', 5.): 12.97, # Table 4 of [1]
('clamped', 1.): 5.60, ('clamped', 5.): 7.23}
MH_QUAD4_PLATE = {('simple', 1.): [0.981, 1.004, 1.003, 1.002], # Table 11 of [1]
('simple', 5.): [1.052, 0.991, 0.997, 0.998],
('clamped', 1.): [0.934, 1.010, 1.012, 1.010], # Table 12 of [1]
('clamped', 5.): [0.519, 0.863, 0.940, 0.972]}
NS = (2, 4, 6, 8)
def rect_plate(name, n, ratio, support):
a, t, E, nu = 2., 0.0001, 1.7472e7, 0.3
b = ratio*a
X, Y = np.meshgrid(np.linspace(0., a/2, n + 1), np.linspace(0., b/2, n + 1),
indexing='ij')
ncoords = np.column_stack((X.ravel(), Y.ravel(), np.zeros(X.size)))
ids, quads = grid(n, n)
prop = isotropic_plate(thickness=t, E=E, nu=nu)
K, conn, areas = stiffness(ncoords, quads, name, prop)
N = K.shape[0]
x, y = ncoords[:, 0], ncoords[:, 1]
bk = np.zeros(N, dtype=bool)
for d in (0, 1, 5):
bk[d::DOF] = True
edge = np.isclose(x, 0.) | np.isclose(y, 0.)
bk[2::DOF] |= edge
if support == 'clamped':
bk[3::DOF] |= edge
bk[4::DOF] |= edge
bk[4::DOF] |= np.isclose(x, a/2) # symmetry, rotation about y
bk[3::DOF] |= np.isclose(y, b/2) # symmetry, rotation about x
centre = ids[n, n]
if support == 'simple':
f = area_load(ncoords, conn, areas, (0., 0., 1.e-4))
else:
f = np.zeros(N)
f[DOF*centre + 2] = 4.e-4/4
u = solve(K, bk, f)
return u[DOF*centre + 2]
plate = {}
for (support, ratio), ref in PLATE_REF.items():
print('%s, b/a = %g, normalised deflection at the centre'
% ('simple supports, uniform load' if support == 'simple'
else 'clamped supports, central load', ratio))
print('%-4s' % 'N' + ''.join('%10s' % n for n in ELEMENTS) + ' MH QUAD4')
for i, n in enumerate(NS):
for name in ELEMENTS:
plate[name, support, ratio, n] = rect_plate(name, n, ratio, support)/ref
print('%-4d' % n + ''.join('%10.3f' % plate[name, support, ratio, n] for name in ELEMENTS)
+ '%10.3f' % MH_QUAD4_PLATE[support, ratio][i])
print()
simple supports, uniform load, b/a = 1, normalised deflection at the centre
N Quad4 Quad4R Tria3R Tria3DSG MH QUAD4
2 1.065 0.781 1.410 0.883 0.981
4 1.016 0.930 1.069 0.989 1.004
6 1.007 0.967 0.995 1.009 1.003
8 1.004 0.981 0.965 1.017 1.002
simple supports, uniform load, b/a = 5, normalised deflection at the centre
N Quad4 Quad4R Tria3R Tria3DSG MH QUAD4
2 0.826 0.875 4.583 0.922 1.052
4 0.978 0.974 1.868 1.115 0.991
6 0.991 0.988 1.391 1.155 0.997
8 0.995 0.993 1.225 1.170 0.998
clamped supports, central load, b/a = 1, normalised deflection at the centre
N Quad4 Quad4R Tria3R Tria3DSG MH QUAD4
2 0.865 0.533 2.286 0.000 0.934
4 0.965 0.819 1.401 0.000 1.010
6 0.985 0.907 1.170 0.002 1.012
8 0.992 0.944 1.074 0.009 1.010
clamped supports, central load, b/a = 5, normalised deflection at the centre
N Quad4 Quad4R Tria3R Tria3DSG MH QUAD4
2 0.318 0.287 8.366 0.000 0.519
4 0.825 0.666 3.942 0.000 0.863
6 0.911 0.799 2.591 0.001 0.940
8 0.946 0.867 1.994 0.003 0.972
Quad4 converges in all four cases, more slowly than the QUAD4 of [1] for the clamped plate with \(b/a = 5\). Quad4R converges from below, more slowly. The two triangles are the interesting cases at this thickness, \(\ell/t\) between 1250 and 25000:
Tria3Ris too flexible on the coarse meshes, up to eight times on the elongated elements of \(b/a = 5\), and is not converged at \(N = 8\). This is the transverse shear stiffness divided by \(1 + \alpha \ell^2/t^2\) with the defaultalpha_shear_locking = 0.7, whose error behaves as \(O(1/N^2)\) (Notes ofpyfe3d.tria3r).Tria3DSGlocks with clamped supports and the central load, at less than 1% of the reference up to \(N = 8\), and with simple supports and the uniform load it is too flexible on the elongated elements of \(b/a = 5\), by 17% at \(N = 8\). Both converge under further refinement, computed below, but slowly.
[7]:
print('%-36s %4s %10s %10s' % ('case', 'N', 'Quad4', 'Tria3DSG'))
for support, ratio, ns in (('clamped', 1., (16, 32)), ('simple', 5., (16, 32, 64))):
for n in ns:
print('%-36s %4d' % ('%s, b/a = %g' % (support, ratio), n)
+ ''.join('%10.3f' % (rect_plate(name, n, ratio, support)/PLATE_REF[support, ratio])
for name in ('Quad4', 'Tria3DSG')))
case N Quad4 Tria3DSG
clamped, b/a = 1 16 0.999 0.328
clamped, b/a = 1 32 1.001 0.913
simple, b/a = 5 16 0.999 1.177
simple, b/a = 5 32 1.000 1.033
simple, b/a = 5 64 1.000 1.000
Summary#
The grades of Table 16 of [1] for the tests modelled here, the absolute errors of the cases of each test being averaged. For the patch tests the error is the larger of the two measures above, in per cent. As a check of the grading, the same function applied to the QUAD4 results of Tables 8 to 12 of [1] must give the QUAD4 column of their Table 16.
[8]:
def summary(r):
# r(test) returns the normalised results of the cases of each test
tests = [
('(1) Patch test, membrane', r('patch', 'membrane')),
('(2) Patch test, bending', r('patch', 'bending')),
('(3) Straight beam, extension', r('straight', 'extension', ('regular', 'trapezoidal', 'parallelogram'))),
('(4) Straight beam, in-plane, regular', r('straight', 'in-plane shear', ('regular',))),
('(5) Straight beam, in-plane, irregular', r('straight', 'in-plane shear', ('trapezoidal', 'parallelogram'))),
('(6) Straight beam, out-of-plane, regular', r('straight', 'out-of-plane shear', ('regular',))),
('(7) Straight beam, out-of-plane, irregular', r('straight', 'out-of-plane shear', ('trapezoidal', 'parallelogram'))),
('(8) Straight beam, twist', r('straight', 'twist', ('regular', 'trapezoidal', 'parallelogram'))),
('(9) Curved beam, in-plane', r('curved', 'in-plane')),
('(10) Curved beam, out-of-plane', r('curved', 'out-of-plane')),
('(11) Twisted beam', r('twisted')),
('(12) Rectangular plate (N = 4)', r('plate')),
]
return [(label, grade(v) if v is not None else '-') for label, v in tests]
def ours(name):
def r(test, *args):
if test == 'patch':
# error in per cent turned into a normalised value
return 1. + max(abs(e) for e in patch[name, args[0]])
if test == 'straight':
load, shapes = args
return [straight[name, s, load] for s in shapes]
if test == 'curved':
return curved[name, args[0]]
if test == 'twisted':
return [twisted[name, 12, l] for l in TWISTED_REF]
if test == 'plate':
return [plate[name, s, ra, 4] for (s, ra) in PLATE_REF]
return r
def macneal_harder_quad4(test, *args):
if test == 'patch':
return 1. # Table 7 of [1]: 0.0% in both tests
if test == 'straight':
load, shapes = args
return [MH_QUAD4_STRAIGHT[s][load] for s in shapes]
if test == 'curved':
return MH_QUAD4_CURVED[args[0]]
if test == 'twisted':
return list(MH_QUAD4_TWISTED.values())
if test == 'plate':
return [MH_QUAD4_PLATE[k][1] for k in PLATE_REF]
# QUAD4 column of Table 16 of [1], tests (1) to (12)
TABLE16_QUAD4 = list('AAABFABBCBAB')
check = [g for _, g in summary(macneal_harder_quad4)]
print('grading of the QUAD4 results of [1] reproduces their Table 16:',
'PASS' if check == TABLE16_QUAD4 else 'FAIL %s' % check)
print()
grades = {name: summary(ours(name)) for name in ELEMENTS}
labels = [label for label, _ in grades['Quad4']]
print('%-44s' % 'test' + ''.join('%10s' % n for n in ELEMENTS) + ' MH QUAD4')
for i, label in enumerate(labels):
print('%-44s' % label + ''.join('%10s' % grades[n][i][1] for n in ELEMENTS)
+ '%10s' % TABLE16_QUAD4[i])
print('%-44s' % 'failed tests (D and F)'
+ ''.join('%10d' % sum(g in 'DF' for _, g in grades[n]) for n in ELEMENTS)
+ '%10d' % sum(g in 'DF' for g in TABLE16_QUAD4))
grading of the QUAD4 results of [1] reproduces their Table 16: PASS
test Quad4 Quad4R Tria3R Tria3DSG MH QUAD4
(1) Patch test, membrane A A A A A
(2) Patch test, bending F F D A A
(3) Straight beam, extension A A A A A
(4) Straight beam, in-plane, regular B F F F B
(5) Straight beam, in-plane, irregular D F F F F
(6) Straight beam, out-of-plane, regular A B C B A
(7) Straight beam, out-of-plane, irregular A D D B B
(8) Straight beam, twist F F F D B
(9) Curved beam, in-plane B F F F C
(10) Curved beam, out-of-plane F F B D B
(11) Twisted beam F F C C A
(12) Rectangular plate (N = 4) B C F F B
failed tests (D and F) 5 8 7 6 1
With their default attributes, the four elements of pyfe3d fail between five and eight of these twelve tests, against one for the QUAD4 of [1]. The failures of Quad4 are the bending patch test, the in-plane bending of the distorted straight beam (D, where the QUAD4 of [1] gets an F), the two tests that twist a strip one element wide, and the warped elements of the twisted beam. Quad4R and the two triangles lock in in-plane bending on these one-element-deep meshes, and Tria3DSG is
the only element that passes both patch tests. Several failures are confined to the coarse meshes that [1] prescribes on purpose: the refinements above show the twisted beam converging for every element and the locked plates of Tria3DSG recovering.
The Scordelis-Lo roof and the spherical shell of [1] are run, with Tables 13 and 14 of [1], in the notebook of the shell obstacle course.
References#
[1] MacNeal, R. H., and Harder, R. L., A proposed standard set of problems to test finite element accuracy, Finite Elements in Analysis and Design, 1(1), 3-20, 1985. https://doi.org/10.1016/0168-874X(85)90003-4
[2] Allman, D. J., A compatible triangular element including vertex rotations for plane elasticity analysis, Computers & Structures, 19(1-2), 1-8, 1984. https://doi.org/10.1016/0045-7949(84)90197-4
[3] Hughes, T. J. R., Taylor, R. L., and Kanoknukulchai, W., A simple and efficient finite element for plate bending, International Journal for Numerical Methods in Engineering, 11(10), 1529-1543, 1977. https://doi.org/10.1002/nme.1620111005
[4] MacNeal, R. H., A simple quadrilateral shell element, Computers & Structures, 8(2), 175-183, 1978. https://doi.org/10.1016/0045-7949(78)90020-2