Repository navigation
Shell analysis
extensions/structure/shell.mojo gives a three-node flat shell for linear-elastic analysis. It adds a constant-strain membrane, the discrete Kirchhoff plate triangle of Batoz, Bathe and Ho, and a small drilling stiffness. Shells and frame members share nodes in one StructuralModel. Frame analysis covers the model, the solvers and the units.
Add three nodes, then a shell over them. The normal follows the right hand rule over the three corners.
var a = model.add_node(Vec3d(0, 0, 3))
var b = model.add_node(Vec3d(4, 0, 3))
var c = model.add_node(Vec3d(0, 4, 3))
var slab = model.add_shell(a, b, c, Length64(0.2, METER), concrete())
model.add_pressure(live, slab, Pressure64(2.4, KILOPASCAL))A positive pressure pushes against the normal. Each corner takes a third of the resultant, with no moments. add_self_weight puts a third of ρ t A g at each corner, along -z.
After a static solve, result.shell_resultants[s] gives the shell's local axes, membrane stresses and moments. model.shell_element(s) gives the element itself, with its stiffness, transformation and corner coordinates.
The local x axis runs from the first corner to the second. The local z axis is the normal. The local y axis is z × x. Each corner has six degrees of freedom, the same as a frame node.
| Part | Formulation |
|---|---|
| Membrane | The constant-strain triangle, A t Bᵀ D B, with D for plane stress |
| Plate | The discrete Kirchhoff triangle (DKT). The rotations are quadratic, and the Kirchhoff constraint holds at the corners and the side midpoints. |
| Plate integration | Three interior Gauss points, at (1/6, 1/6), (2/3, 1/6) and (1/6, 2/3), with weight A / 3. They integrate the DKT exactly. |
| Drilling | A penalty ties each corner's rotation about the normal to the membrane rotation ω = (∂v/∂x - ∂u/∂y) / 2. Its modulus is 0.001 G t. |
| Transformation | Six copies of the 3 by 3 rotation. A global matrix is Tᵀ K T. |
| Mass | Lumped: ρ t A / 3 on each translation of each corner, and that mass times t² / 12 on each rotation |
The rotation about the local x axis is ∂w/∂y. The rotation about the local y axis is -∂w/∂x. These are right-handed rotations, the same as the frame element's.
The drilling penalty follows Hughes and Brezzi. A rigid rotation gives no drilling strain, so the patch tests stay exact. A flat mesh has no free drilling mode, so a coplanar slab is not a mechanism.
ShellResultants gives the values in the shell's local axes.
| Field | Meaning |
|---|---|
sigma_x, sigma_y, tau_xy
|
The membrane stresses, constant over the element |
m_x |
-D (∂²w/∂x² + ν ∂²w/∂y²), at the centroid, in newton meters per meter |
m_y |
-D (∂²w/∂y² + ν ∂²w/∂x²), at the centroid |
m_xy |
-D (1 - ν) ∂²w/∂x∂y, at the centroid |
D is E t³ / (12 (1 - ν²)). The moment fields carry the Force64 type, because a moment per length has the dimension of a force.
A flat plate model with only out-of-plane loads still has membrane degrees of freedom. Fix UX, UY and RZ at enough nodes to stop the in-plane rigid motions, or at every node for a pure bending problem. For a simply supported edge, fix UZ at the edge nodes. For a clamped edge, also fix RX and RY.
Each row is a test in tests/test_structure_shell.mojo. The error is the measured relative difference from the reference.
| Case | Reference | Error |
|---|---|---|
| Membrane patch test, irregular patch, uniform tension | σx = σ, σy = τxy = 0 in every element; u = σ x / E, v = -ν σ y / E | Below 1e-6 |
| Plate patch test, irregular patch, constant curvature | Interior nodes from the field; mx, my and mxy constant | Below 1e-9 |
| Rigid translation and rotation of a tilted element | No force, no stress | Below 1e-9 |
| Tilted shell under pressure and its own weight | Reactions sum to p A n plus ρ t A g | Below 1e-9 |
| Simply supported square plate, uniform load, 4, 8 and 16 divisions per side | w = 0.00406 q a⁴ / D, ν = 0.3 | -4.2%, -1.3%, -0.29% |
| Clamped square plate, uniform load, 4, 8 and 16 divisions per side | w = 0.00126 q a⁴ / D, ν = 0.3 | +10.2%, +3.0%, +1.1% |
| Simply supported square plate, first frequency, 12 divisions | ω = 2 π² / a² sqrt(D / (ρ t)) | -0.33% |
Each plate mesh splits each square into two triangles. The plate tests fix UX, UY and RZ at every node.
- The element is flat. A curved shell needs many small elements.
- The membrane is a constant-strain triangle. It is stiff in in-plane bending, so a shear wall or a deep beam needs a fine mesh.
- The material is isotropic and in plane stress. Layers of a construction are not modeled one by one.
- The plate has no transverse shear deformation. It is for thin plates.
- A pressure is uniform over each element and lumped at the corners.
- The mass is lumped, not consistent.
- Batoz, Bathe and Ho, "A study of three-node triangular plate bending elements", International Journal for Numerical Methods in Engineering 15, 1980.
- Hughes and Brezzi, "On drilling degrees of freedom", Computer Methods in Applied Mechanics and Engineering 72, 1989.
- Irons and Razzaque, "Experience with the patch test for convergence of finite elements", 1972.
- Timoshenko and Woinowsky-Krieger, "Theory of Plates and Shells", 2nd edition, McGraw-Hill, 1959: tables 8 and 35.
- Leissa, "Vibration of Plates", NASA SP-160, 1969.
- Cook, Malkus, Plesha and Witt, "Concepts and Applications of Finite Element Analysis", 4th edition, 2002: chapter 15.