NGSolve Technical Note 2026-04 \(\bullet\) October 2026
A flat circular sheet \(S \subset \mathbb{R}^3\) of radius \(1.2\), with a slightly ragged hem \(r(\theta) = 1.2\,(1 + 0.04\sin 11\theta + 0.025\sin(23\theta+1))\) and three elliptic holes (eyes and mouth), starts at rest at height \(z_0 = 0.6\) above a sphere (“head”) of radius \(R_h = 0.35\). The head is at rest until \(t = 2.5\), then follows a smoothly ramped path \[\boldsymbol c(t) = \bigl(1-e^{-s/0.7}\bigr)\bigl(0.55\sin 1.3s,\ 0.15\sin 0.9s,\ 0.45\sin^2 0.7s\bigr), \qquad s = t - 2.5 .\] The sheet is a surface mesh in 3D with 3804 triangles (\(h = 0.05\)). The unknown is the displacement \(\boldsymbol u: S \to \mathbb{R}^3\); the deformed position is \(\boldsymbol X+ \boldsymbol u\).
With the unit normal \(\boldsymbol n\) of the reference surface and the tangential projection \(\boldsymbol P= I - \boldsymbol n\otimes\boldsymbol n\), the surface deformation gradient and Green strain are \[F = \nabla_\tau \boldsymbol u+ \boldsymbol P, \qquad E = \tfrac12 (F^T F - \boldsymbol P),\] and the deformed normal is \(\boldsymbol n_{\mathrm{def}} = \operatorname{cof}(F)\boldsymbol n/ |\operatorname{cof}(F)\boldsymbol n|\).
To avoid membrane locking, \(E\) is interpolated into lowest-order Regge elements (\(H(\mathrm{curl\,curl})\), order 1), \(E_r = \mathcal I^{\mathrm{Regge}} E\). The membrane energy is compressible neo-Hooke, \[W_{\mathrm{mem}} = \frac{\mu}{2}(\operatorname{tr} C - 2) - \mu \ln J + \frac{\lambda}{2} (\ln J)^2, \qquad C = 2E_r + \boldsymbol P,\quad J^2 = \det(C + \boldsymbol n\otimes\boldsymbol n),\] with \(\mu = 500\), \(\lambda = 200\) (per unit area). St. Venant–Kirchhoff was used first, but it loses convexity for stretches below \(1/\sqrt3\), which happens at the tips of the holes when the sheet is compressed; Newton then stalled.
We use the nonlinear HHJ shell element of Neunteufel and Schöberl. The bending moment \(\sigma\) is an independent unknown in \(H(\mathrm{div\,div})\) on the surface. The bending part of the Lagrangian is \[\sum_{T} \int_T \Bigl( -\tfrac12\, \mathcal M^{-1}\sigma : \sigma + \sigma : H_\nu(\boldsymbol u) \Bigr) + \sum_T \int_{\partial T} \Bigl( \sphericalangle(\boldsymbol\mu_{\mathrm{def}}, \bar{\boldsymbol\mu}) - \tfrac{\pi}{2} + \alpha\cdot\boldsymbol\mu \Bigr)\, \sigma_{\mu\mu},\] where
\(H_\nu(\boldsymbol u) = \sum_i (\nabla^2 \boldsymbol u_i)\,(\boldsymbol n_{\mathrm{def}})_i\) is the Hessian of the deformation contracted with the deformed normal;
\(\boldsymbol\mu\) is the outer conormal of the element edge; \(\boldsymbol\mu_{\mathrm{def}} = \boldsymbol n_{\mathrm{def}} \times \boldsymbol t_{\mathrm{def}}\) is the deformed one;
\(\bar{\boldsymbol\mu}\) is an averaged normal on each edge (a vector facet field, order 2), updated once per time step from the current deformed normal and projected orthogonal to the deformed edge tangent. The edge angle measures the kink of the deformed surface across edges; \(\pi/2\) is its value for the flat reference sheet;
the material law is \(\mathcal M^{-1}\sigma:\sigma = \frac{1+\nu}{D}\bigl(|\sigma|^2 - \frac{\nu}{1+\nu}(\operatorname{tr}\sigma)^2\bigr)\), i.e. the Kirchhoff plate law with \(D = Et^3/12\); here \(D = 0.05\), \(\nu = 0.3\).
\(\sigma\) is taken
discontinuous (HDivDivSurface with
discontinuous=True), and the normal-normal continuity of
\(\sigma\) is restored by the Lagrange
multiplier \(\alpha\) in
NormalFacetSurface (order 1), which plays the role of the
rotation across edges. On the free boundary (hem and hole rims) \(\alpha\) has no boundary condition, which
gives \(\sigma_{\mu\mu} = 0\) weakly.
Since \(\sigma\) is now element-local,
it is eliminated by static condensation; the condensed system in \((\boldsymbol u, \alpha)\) is symmetric
positive definite and is solved with a sparse Cholesky factorization.
(The non-hybridized mixed form is a saddle point problem and needed a
general LU.)
Sanity check: a cantilever strip under small load gives a tip deflection of \(-0.024988\) vs. \(-qL^4/(8D) = -0.025\).
Contact with the analytic sphere by a penalty on the signed distance: \(\frac{k_c}{2} [R_h - |\boldsymbol X+ \boldsymbol u- \boldsymbol c(t)|]_+^2\), \(k_c = 2\cdot 10^4\).
Contact friction as viscous drag relative to the head velocity \(\boldsymbol v_h\), active where the sheet was within \(0.02\) of the head in the previous step (\(\chi^n \in \{0,1\}\)), \(c_f = 200\).
Air drag everywhere, \(c_a = 3\). Without it the heavy hem swings under the head, folds through itself and an element inverts (\(\ln J\) undefined).
Gravity \(g = 9.81\); areal density \(\rho = 1\), and \(\rho = 2\) in the hem ring \(r > 1\).
Implicit Euler written as an incremental energy: in each step we find the stationary point of \[\begin{aligned} \Pi^{n+1}(\boldsymbol u,\sigma,\alpha) = \int_S \Bigl[\, & \frac{\rho}{2\Delta t^2}\,|\boldsymbol u- \boldsymbol u^n - \Delta t\, \boldsymbol v^n|^2 + W_{\mathrm{mem}}(\boldsymbol u) + \rho g\, u_z + \frac{k_c}{2} [R_h - |\boldsymbol X+ \boldsymbol u- \boldsymbol c^{n+1}|]_+^2 \\ & + \frac{c_f \chi^n + c_a}{2\Delta t}\, |\boldsymbol u- \boldsymbol u^n - \Delta t\, \boldsymbol v_h|^2 \Bigr] + \text{(HHJ bending terms)}, \end{aligned}\] then set \(\boldsymbol v^{n+1} = (\boldsymbol u^{n+1} - \boldsymbol u^n)/\Delta t\). The step is \(\Delta t = 0.01\), 800 steps for \(t \in [0, 8]\). Each step uses Newton on the condensed system, starting from the predictor \(\boldsymbol u^n + \Delta t\,\boldsymbol v^n\). A Newton update is halved while the residual is not finite (inverted element); a failed step would be retried with \(\Delta t/2\), recursively down to \(\Delta t/16\). In the final run no step needed this: 2–6 Newton iterations per step.
| displacement \(\boldsymbol u\) | VectorH1, order 2 |
| moment \(\sigma\) | HDivDivSurface, order 1,
discontinuous (condensed) |
| multiplier \(\alpha\) | NormalFacetSurface, order
1 |
| membrane strain | Regge interpolation,
HCurlCurl order 1 |
| averaged normal \(\bar{\boldsymbol\mu}\) | VectorFacetSurface, order
2 |
| unknowns | 69 426 total, 35 190 after condensation |
| linear solver | sparse Cholesky |
Each energy term is its own Variation; a single
summed Variation linearized about \(20\times\) slower. The membrane, bending
and edge terms are Compile()d.
Assembling the linearization takes about 0.7 s per Newton iteration and dominates; the Cholesky factorization takes about 0.05 s.
The total run is roughly 30–80 minutes on an Apple M4 (depending on threads).
Visualization: the order-2 displacement is drawn with the NGSolve webgui as curved Bernstein triangles; frames are captured with headless Chrome and encoded to H.264.
M. Neunteufel, J. Schöberl, The Hellan–Herrmann–Johnson method for nonlinear shells, Computers & Structures 225 (2019); M. Neunteufel, J. Schöberl, Avoiding membrane locking with Regge interpolation, CMAME 373 (2021); see also the ESI minicourse https://jschoeberl.github.io/ESI-Elasticity-Geometry-NGSolve/ESI-minicourse.html.