StructuralMind
Start free trial

Shell element, fundamental tests

Updated:

On this page

Results.

ExampleReferenceRecordsLargest deviation
MacNeal-Harder patch, quadrilateral elementsPublished source7< 0.001%
MacNeal-Harder patch, triangular elementsPublished source7< 0.001%
MacNeal-Harder patch, mixed quadrilateral and triangular meshPublished source7< 0.001%
Quadrilateral shell element, eigenvalue testClosed form50%
Triangular shell element, eigenvalue testClosed form50%

Problem. For a finite element to converge to the right answer as the mesh is refined it must meet two preconditions. On a patch of a few irregular elements it must reproduce a constant stress and a constant curvature field exactly, and in a rigid body motion it must produce no internal force at all. This page tests the quadrilateral and triangular shell elements against these conditions.

Reference document. The patch test is the standard patch of MacNeal and Harder (1985), and the exact solution is the constant field imposed on it. For the eigenvalue test the reference values are zero internal force under rigid-body motion and six zero-energy modes.

In the membrane patch test this displacement field is imposed on the outer nodes and these constant section forces are expected in every element. As the stress is σ=N/t\sigma = N/t, this gives σx=σy=1333.3\sigma_x = \sigma_y = 1333.3 and τxy=400\tau_{xy} = 400.

u=10−3(x+y2),v=10−3(y+x2)⇒Nxx=Nyy=1.333,Nxy=0.400u = 10^{-3}\left(x + \tfrac{y}{2}\right),\quad v = 10^{-3}\left(y + \tfrac{x}{2}\right)\quad\Rightarrow\quad N_{xx} = N_{yy} = 1.333,\quad N_{xy} = 0.400

In the bending patch test this deflection field and its slopes are imposed on the outer nodes and these constant moments are expected in every element. As the surface stress is σ=6M/t2\sigma = 6M/t^2, their magnitudes are 0.667 and 0.200.

w=10−3 x2+xy+y22⇒Mxx=Myy=−1.111×10−7,Mxy=−3.333×10−8w = 10^{-3}\,\dfrac{x^2 + xy + y^2}{2}\quad\Rightarrow\quad M_{xx} = M_{yy} = -1.111 \times 10^{-7},\quad M_{xy} = -3.333 \times 10^{-8}

Numerical model. The patch is 0.24 × 0.12 and consists of five elements around four irregular inner nodes. In the triangular mesh every quadrilateral is split into two triangles, and in the mixed mesh two of them are. The eigenvalue test is run on five randomly rotated elements distorted in their plane and on an unsupported 3 × 3 mesh, and the model shown is the modal solution of that mesh.

MacNeal-Harder patch, quadrilateral elements

Parameters. E = 1,000,000 · ν = 0.25 · t = 0.001

Published source. MacNeal and Harder (1985), A proposed standard set of problems to test finite element accuracy, patch test

QuantityReferenceComputedDeviation
Membrane force Nxx (F11)1.333331.333330%
Membrane force Nyy (F22)1.333331.333330%
Membrane shear Nxy (F12)0.40.40%
Bending moment Mxx (M11)−1.111×10−7-1.111 \times 10^{-7}−1.111×10−7-1.111 \times 10^{-7}< 0.001%
Bending moment Myy (M22)−1.111×10−7-1.111 \times 10^{-7}−1.111×10−7-1.111 \times 10^{-7}< 0.001%
Twisting moment Mxy (M12)−3.333×10−8-3.333 \times 10^{-8}−3.333×10−8-3.333 \times 10^{-8}< 0.001%
Largest displacement error at the inner nodes07.603×10−167.603 \times 10^{-16}0%

Note. The bending moments are small because of the plate thickness (t = 0.001), since the moment scales with the cube of the thickness. The published value is the surface stress, which the relation σ=6M/t2\sigma = 6M/t^2 turns into 0.667 and 0.200.

MacNeal-Harder patch, triangular elements

Parameters. E = 1,000,000 · ν = 0.25 · t = 0.001

Published source. MacNeal and Harder (1985), A proposed standard set of problems to test finite element accuracy, patch test

QuantityReferenceComputedDeviation
Membrane force Nxx (F11)1.333331.333330%
Membrane force Nyy (F22)1.333331.333330%
Membrane shear Nxy (F12)0.40.40%
Bending moment Mxx (M11)−1.111×10−7-1.111 \times 10^{-7}−1.111×10−7-1.111 \times 10^{-7}< 0.001%
Bending moment Myy (M22)−1.111×10−7-1.111 \times 10^{-7}−1.111×10−7-1.111 \times 10^{-7}< 0.001%
Twisting moment Mxy (M12)−3.333×10−8-3.333 \times 10^{-8}−3.333×10−8-3.333 \times 10^{-8}< 0.001%
Largest displacement error at the inner nodes01.359×10−161.359 \times 10^{-16}0%

Note. The bending moments are small because of the plate thickness (t = 0.001), since the moment scales with the cube of the thickness. The published value is the surface stress, which the relation σ=6M/t2\sigma = 6M/t^2 turns into 0.667 and 0.200.

MacNeal-Harder patch, mixed quadrilateral and triangular mesh

Parameters. E = 1,000,000 · ν = 0.25 · t = 0.001

Published source. MacNeal and Harder (1985), A proposed standard set of problems to test finite element accuracy, patch test

QuantityReferenceComputedDeviation
Membrane force Nxx (F11)1.333331.333330%
Membrane force Nyy (F22)1.333331.333330%
Membrane shear Nxy (F12)0.40.40%
Bending moment Mxx (M11)−1.111×10−7-1.111 \times 10^{-7}−1.111×10−7-1.111 \times 10^{-7}< 0.001%
Bending moment Myy (M22)−1.111×10−7-1.111 \times 10^{-7}−1.111×10−7-1.111 \times 10^{-7}< 0.001%
Twisting moment Mxy (M12)−3.333×10−8-3.333 \times 10^{-8}−3.333×10−8-3.333 \times 10^{-8}< 0.001%
Largest displacement error at the inner nodes08.146×10−168.146 \times 10^{-16}0%

Note. The bending moments are small because of the plate thickness (t = 0.001), since the moment scales with the cube of the thickness. The published value is the surface stress, which the relation σ=6M/t2\sigma = 6M/t^2 turns into 0.667 and 0.200.

Quadrilateral shell element, eigenvalue test

Parameters. element count = 5

Closed form. Properties every correct finite element must have

QuantityReferenceComputedDeviation
Symmetry error of the stiffness matrix max⁡∣K−KT∣/max⁡∣K∣\max\lvert K - K^{\mathsf T}\rvert / \max\lvert K\rvert01.037×10−161.037 \times 10^{-16}0%
Internal force in rigid body motion max⁡∣f∣/max⁡∣K∣\max\lvert f\rvert / \max\lvert K\rvert08.299×10−168.299 \times 10^{-16}0%
Number of zero-energy modes, single element660%
Number of zero-energy modes, free mesh660%
Number of zero-frequency modes, modal solution of the free mesh660%

Triangular shell element, eigenvalue test

Parameters. element count = 5

Closed form. Properties every correct finite element must have

QuantityReferenceComputedDeviation
Symmetry error of the stiffness matrix max⁡∣K−KT∣/max⁡∣K∣\max\lvert K - K^{\mathsf T}\rvert / \max\lvert K\rvert06.091×10−176.091 \times 10^{-17}0%
Internal force in rigid body motion max⁡∣f∣/max⁡∣K∣\max\lvert f\rvert / \max\lvert K\rvert01.233×10−151.233 \times 10^{-15}0%
Number of zero-energy modes, single element660%
Number of zero-energy modes, free mesh660%
Number of zero-frequency modes, modal solution of the free mesh660%

Note. In a single triangular element the seventh eigenvalue is about 10⁻⁸ of the largest. This mode comes from the softened in-plane twisting stiffness of the MITC3+ element. It cannot pass from one element to the next, so the free mesh has exactly six zero-energy modes.