超弹性梁（新胡克）¶
示例库中最简单的有限应变示例：一根 \(1.0 \times 0.4 \times 0.4\) m 的橡胶梁，一端固定，另一端由扭转力场扭转。采用可压缩新胡克应变能密度和 L-BFGS 能量最小化求解。脚本位于 examples/solid/hyperelastic_beam/hyperelastic_beam.py。
为什么在学习下面更难的示例之前先研究这个示例：它在一个干净的几何上单独呈现了 能量最小化方法，没有孔洞、没有接触、也没有辅助场——因此本示例与 cantilever_beam.py 之间的差异，恰好就是从“线性求解”到“在非凸能量上做 L-BFGS”的转变。
问题¶
可压缩新胡克应变能：
其中 \(\mathbf{F} = \mathbf{I} + \nabla\mathbf{u}\) 是变形梯度，\(I_1 = \mathrm{tr}(\mathbf{F}^T \mathbf{F})\)，\(J = \det\mathbf{F}\)。拉梅参数 \((\mu, \lambda)\) 来自 tensormesh.material 中的 Rubber.lame_params（\(E \approx 10\) MPa，\(\nu \approx 0.48\)，近乎不可压缩）。
边界条件：
在 \(x = 0\) 处固定：\(\mathbf{u} = \mathbf{0}\)，
在 \(x = 1\) 处的扭转载荷：诺伊曼面力 \(\mathbf{t} = C\, (0, -z, y)\)，其中 \(C = 2.4 \times 10^7\)，通过
FacetAssembler在端面上积分，得到一致节点载荷 \(f_i = \int_\Gamma N_i\, \mathbf{t}\, \mathrm{d}A\)，并在 10 个载荷步内逐步加载，其余所有表面：无面力（自由）。
TensorMesh 设置¶
NeoHookeanModel 是一个自定义的 ElementAssembler，其 element_energy 在每个求积点返回 \(\Psi(\mathbf{F})\)——仅此而已。TensorMesh 在各单元上积分，求和得到一个标量总能量，随后 PyTorch 的 autograd 为 L-BFGS 提供梯度：
class NeoHookeanModel(ElementAssembler):
def __post_init__(self, mu, lam):
self.mu = mu
self.lam = lam
def element_energy(self, gradu):
dim = gradu.shape[-1]
I = torch.eye(dim, device=gradu.device, dtype=gradu.dtype)
F = I + gradu
J = torch.clamp(torch.det(F), min=1e-6)
I1 = (F**2).sum()
logJ = torch.log(J)
return (0.5 * self.mu * (I1 - 3)
- self.mu * logJ
+ 0.5 * self.lam * logJ**2)
mesh = gen_cube(chara_length=0.05, order=2,
left=0, right=1.0, bottom=0, top=0.4,
front=0, back=0.4)
model = NeoHookeanModel.from_mesh(mesh, mu=mu, lam=lam)
u = torch.zeros_like(mesh.points, requires_grad=True)
optimizer = optim.LBFGS([u], lr=1.0, max_iter=50,
history_size=50,
line_search_fn="strong_wolfe")
for step in range(1, n_load_steps + 1):
lam_load = step / n_load_steps # incremental load
def closure():
optimizer.zero_grad()
u_active = apply_bcs(u, lam_load)
E = model.energy(point_data={"displacement": u_active})
E.backward()
return E
optimizer.step(closure)
几个值得注意的细节：
二次四面体单元（
order=2）。线性（P1）四面体在近乎不可压缩极限（\(\nu \to 0.5\)）下会发生 闭锁。对本问题而言，P2 单元是开销最小的解决办法。对 J 进行截断。 当 L-BFGS 过冲进入翻折构型时，
torch.clamp(J, min=1e-6)能让能量保持有限。一种更优雅但更慢的替代方案是添加一个障碍项，把 \(J\) 推回 \((0, \infty)\)。通过掩码施加边界条件。 脚本使用
u_active = u * mask + val这一模式在闭包内部施加狄利克雷边界条件——被掩码的自由度上的梯度自动为零，因此 L-BFGS 不会去动它们。这就无需单独创建一个Condenser实例，并且能与 autograd 顺畅协作。增量加载是必须的。 否则，从静止构型热启动的 L-BFGS 通常无法在满扭转下找到对应的能量盆地。对这根梁来说 10 步就足够了；更刚硬的材料或更剧烈的加载则需要更多步。
图 43 hyperelastic_beam.py 的输出：橡胶梁在其最终扭转构型下，按位移大小着色。右端（蓝色立方体）固定；左端的红色箭头表示在 10 个 LBFGS 载荷步内逐步施加的规定扭转载荷。与上面的线性悬臂梁不同，变形以 1.0 倍显示——新胡克响应无需线性化即可处理有限转动，不会产生线性化伪影。¶
运行示例¶
cd examples/solid/hyperelastic_beam
python hyperelastic_beam.py # writes hyperelastic_rubber.png
控制台输出会在每个载荷步报告最大节点位移。