Snap-through of a shallow cylindrical roof¶
The hinged shallow cylindrical roof under a central point load is a classical verification problem for geometrically nonlinear analysis. Its equilibrium path has a limit point: beyond a critical load the structure snaps through to a far, inverted configuration, and the branch in between has negative stiffness.
A load-controlled solve cannot follow that branch, because above the limit load there is no equilibrium near the current state. Viscous stabilization supplies the missing stiffness through the alpha argument of solve(), which adds a damping force $\alpha \mathbf{M} \Delta \mathbf{u} / \Delta \lambda$ to the residual and its tangent. The solve then crosses the valley instead of tracing it, which is what the Abaqus *STATIC, STABILIZE option does.
Geometry and material follow the benchmark of Sabir and Lock, taken here as a plane strain section so the model is a shallow arch rather than the full shell panel. Loads are therefore per unit depth.
import matplotlib.pyplot as plt
import torch
from torchfem import Planar
from torchfem.materials import HyperelasticPlaneStrain
torch.set_default_dtype(torch.float64)
Material¶
We use a Neo-Hookean strain energy density in a hyperelastic plane-strain material.
def psi(F, params):
"""Neo-Hookean strain energy density."""
G, D = params[0], params[1]
C = F.transpose(-1, -2) @ F
J = torch.exp(0.5 * torch.logdet(C))
C_bar = C * J ** (-2.0 / 3.0)
return G * (torch.trace(C_bar) - 3.0) + 1 / D * (J - 1) ** 2
# Neo-Hookean parameters matched to E and nu
E = 3102.75
NU = 0.3
mu = E / (2 * (1 + NU))
kappa = E / (3 * (1 - 2 * NU))
material = HyperelasticPlaneStrain(psi, params=[mu / 2, 2.0 / kappa])
Geometry¶
# Geometry and mesh parameters
R = 2540.0
THETA = 0.1
T = 6.35
N_ARC = 40
N_T = 4
# Nodes on concentric arcs through the thickness, crown at phi=0
phi = torch.linspace(0.0, THETA, N_ARC + 1)
y0 = R * torch.cos(torch.tensor(THETA))
nodes = torch.cat(
[
torch.stack([r * torch.sin(phi), r * torch.cos(phi) - y0], dim=-1)
for r in torch.linspace(R - T / 2, R + T / 2, N_T + 1)
]
)
npl = N_ARC + 1
elements = torch.tensor(
[
[j * npl + i, j * npl + i + 1, (j + 1) * npl + i + 1, (j + 1) * npl + i]
for j in range(N_T)
for i in range(N_ARC)
]
)
# Node columns: crown at phi=0, hinge on the mid-surface at phi=THETA
col = torch.arange(len(nodes)) % npl
row = torch.arange(len(nodes)) // npl
crown = col == 0
hinge = (col == N_ARC) & (row == N_T // 2)
roof = Planar(nodes, elements, material)
roof.thickness = torch.ones(len(elements))
roof.constraints[hinge] = True # hinged support
roof.constraints[crown, 0] = True # symmetry plane
roof.plot(title="Half section of the cylindrical roof")
Reference equilibrium path¶
Prescribing the crown displacement makes the problem stable, so the full path including the descending branch can be traced. This is the reference the load-controlled solve is compared against.
# Apply boundary conditions: crown node is displaced downward by 30 mm
roof.constraints[crown, 1] = True
roof.displacements[crown, 1] = -30.0
# Solve
increments = torch.linspace(0, 1, 31)
u_ref, f_ref, _, _, _ = roof.solve(
increments=increments, nlgeom=True, return_intermediate=True, verbose=True
)
# Postprocess results
w_ref = -u_ref[:, crown, 1].mean(dim=1)
P_ref = -f_ref[:, crown, 1].sum(dim=1) * 2.0 # half model -> full arch
# Limit point: first point where the path stops rising
descending = torch.nonzero((P_ref[1:] - P_ref[:-1]) < 0).ravel()
i_limit = int(descending[0])
P_limit = float(P_ref[i_limit])
─── torch-fem · solve ──────────────────────────────────────────────────────────────────
model Planar · 160 elem · 410 dof · float64
machine AMD EPYC 7763 64-Core Processor · 2 threads · 16 GB RAM
solver spsolve · direct · scipy · cpu
newton rtol 1e-08 · atol 1e-06 · ≤10 it · nlgeom
────────────────────────────────────────────────────────────────────────────────────────
Increment Load factor Steps Iterations Residual Wall time
1 0.03333 1 4 3.00e-11 1.23 s
2 0.06667 1 2 2.55e-07 0.11 s
3 0.1 1 2 3.69e-07 0.11 s
4 0.1333 1 2 9.61e-07 0.11 s
5 0.1667 1 3 3.01e-11 0.14 s
6 0.2 1 3 2.73e-11 0.14 s
7 0.2333 1 3 2.89e-11 0.14 s
8 0.2667 1 3 2.96e-11 0.15 s
9 0.3 1 3 2.76e-11 0.15 s
10 0.3333 1 3 2.47e-11 0.14 s
11 0.3667 1 3 2.99e-11 0.14 s
12 0.4 1 3 3.71e-11 0.14 s
13 0.4333 1 3 1.10e-10 0.14 s
14 0.4667 1 3 5.32e-10 0.14 s
15 0.5 1 3 4.14e-09 0.14 s
16 0.5333 1 3 5.71e-08 0.14 s
17 0.5667 1 4 3.32e-11 0.17 s
18 0.6 1 4 1.78e-09 0.17 s
19 0.6333 1 5 3.31e-11 0.20 s
20 0.6667 1 4 8.78e-11 0.17 s
21 0.7 1 5 2.68e-11 0.20 s
22 0.7333 1 4 2.87e-11 0.18 s
23 0.7667 1 3 8.71e-09 0.14 s
24 0.8 1 3 5.61e-11 0.14 s
25 0.8333 1 3 2.45e-11 0.15 s
26 0.8667 1 3 3.22e-11 0.14 s
27 0.9 1 3 2.49e-11 0.15 s
28 0.9333 1 2 2.93e-07 0.12 s
29 0.9667 1 2 9.13e-08 0.11 s
30 1 1 2 3.48e-08 0.11 s
────────────────────────────────────────────────────────────────────────────────────────
converged · 30 increments · 93 iterations · 5.44 s
Load control with stabilization¶
The same model is reused with the crown displacement released and a point load applied instead, ramped to twice the limit load. Viscous damping keeps the tangent stiffness invertible while the arch traverses the unstable region, and the solve ends on the far branch.
Setting alpha=0 here makes the solve fail: without the damping term there is no equilibrium near the limit point for Newton-Raphson to find, and subdividing the increment only moves closer to the singular tangent.
# Ramp up to twice the limit load
load = 2.0 * P_limit
# Swap displacement control for load control on the same model
roof.constraints[crown, 1] = False
roof.displacements[crown, 1] = 0.0
roof.forces[crown, 1] = -(load / 2.0) / int(crown.sum()) # half model
# Solve with stabilization
u, _, _, _, _ = roof.solve(
increments=increments,
nlgeom=True,
alpha=1e-6,
return_intermediate=True,
verbose=True,
)
# Postprocess results
w = -u[:, crown, 1].mean(dim=1)
P = increments * load
─── torch-fem · solve ──────────────────────────────────────────────────────────────────
model Planar · 160 elem · 410 dof · float64
machine AMD EPYC 7763 64-Core Processor · 2 threads · 16 GB RAM
solver spsolve · direct · scipy · cpu
newton rtol 1e-08 · atol 1e-06 · ≤10 it · nlgeom · stabilized α=1e-06
────────────────────────────────────────────────────────────────────────────────────────
Increment Load factor Steps Iterations Residual Wall time
1 0.03333 1 3 1.27e-10 0.14 s
2 0.06667 1 2 2.68e-08 0.11 s
3 0.1 1 2 8.35e-09 0.11 s
4 0.1333 1 2 1.25e-08 0.11 s
5 0.1667 1 2 2.19e-08 0.11 s
6 0.2 1 2 4.06e-08 0.11 s
7 0.2333 1 2 7.96e-08 0.11 s
8 0.2667 1 2 1.67e-07 0.11 s
9 0.3 1 2 3.84e-07 0.11 s
10 0.3333 1 2 9.87e-07 0.11 s
11 0.3667 1 3 2.84e-11 0.14 s
12 0.4 1 3 2.77e-11 0.14 s
13 0.4333 1 3 5.39e-11 0.15 s
14 0.4667 1 3 4.64e-09 0.14 s
15 0.5 1 4 2.91e-11 0.17 s
16 0.5333 3 18 1.22e-08 0.67 s ↓1 ↑2
17 0.5667 2 16 2.38e-09 0.57 s ↑2
18 0.6 2 8 3.50e-08 0.34 s ↑2
19 0.6333 2 5 1.93e-11 0.25 s ↑2
20 0.6667 1 3 2.19e-11 0.14 s
21 0.7 1 2 2.73e-11 0.11 s
22 0.7333 1 2 2.33e-11 0.11 s
23 0.7667 1 2 2.40e-11 0.11 s
24 0.8 1 2 2.14e-11 0.11 s
25 0.8333 1 2 2.27e-11 0.11 s
26 0.8667 1 2 2.31e-11 0.11 s
27 0.9 1 2 2.15e-11 0.11 s
28 0.9333 1 2 2.19e-11 0.11 s
29 0.9667 1 2 2.23e-11 0.11 s
30 1 1 2 2.59e-11 0.11 s
────────────────────────────────────────────────────────────────────────────────────────
converged · 30 increments · 107 iterations · 5.01 s
Comparison¶
The stabilized solve follows the reference path up to the limit point, cuts across the valley at nearly constant load instead of descending into it, and rejoins the reference on the far branch. The descending branch itself is never visited, which is the price of load control. How far the damping perturbs the result is measured by roof.stabilization_energy, which should stay small compared to the work done by the external load.
plt.plot(w_ref, P_ref, "o-", color="gray", label="Displacement control")
plt.plot(w, P, "o-", color="black", label="Stabilized load control")
plt.xlabel("Crown displacement $w$ in mm")
plt.ylabel("Load $P$ in N/mm")
plt.title("Shallow roof snap-through")
plt.legend()
plt.grid()
