赫兹接触

一个二维接触力学问题：一个圆形弹性压头被压入一块矩形弹性体；界面上的接触压力应当与赫兹解析解相吻合。脚本 examples/solid/hertzian_contact/hertzian_contact.py 用 罚函数法 求解——通过向总势能中添加一个二次罚能量来施加接触约束，然后用 L-BFGS 最小化。

问题

两个线弹性物体（\(E = 1000\) Pa，\(\nu = 0.3\)）：

  • 压头 —— 半径为 \(R = 1.0\) 的完整圆盘，用 Tri6 二次三角形单元划分网格，圆心位于 \((-1.0, 0)\)，使其最右端点恰好位于原点。其最左侧的弧段（\(x < -1.8\)）被固定。

  • 弹性体 —— 矩形 \([0, 1] \times [-1, 1]\)，用 Quad9 二次四边形单元划分网格。右边缘（\(x = 1\)）施加规定位移 \(u_x = -0.15\)（向左推，压向压头）。

接触界面是竖直线 \(x = 0\)。随着弹性体向左移动，压头被压缩，并形成接触斑。

罚函数接触

无摩擦法向接触是约束 \(g(\mathbf{x}) \geq 0\)，其中 \(g\) 是 间隙函数（两个表面之间的有符号距离，分离时为正）。罚函数法用接触能量来代替这一不等式约束

\[E_\text{contact} \;=\; \frac{1}{2}\, k_p \int_{\Gamma_c} \langle -g \rangle_+^2 \,\mathrm{d}\Gamma,\]

其中 \(\langle x \rangle_+ = \max(x, 0)\)，\(k_p = 2 \times 10^6\) 是罚刚度。当 \(k_p \to \infty\) 时约束被精确施加，但问题会变得病态；实际中 \(k_p\) 通常取为比弹性刚度大几个数量级。

在每个 L-BFGS 步中被最小化的总能量为

\[\Pi(\mathbf{u}) \;=\; \int_{\Omega_1} \tfrac12 \boldsymbol{\varepsilon}_1 : \mathbb{C} : \boldsymbol{\varepsilon}_1 \;+\; \int_{\Omega_2} \tfrac12 \boldsymbol{\varepsilon}_2 : \mathbb{C} : \boldsymbol{\varepsilon}_2 \;+\; E_\text{contact}.\]

前两项来自 LinearElasticityElementAssembler 的能量形式（小应变二次能量）；接触项则在脚本中单独实现。

点对段检测

界面 \(\Gamma_c\) 被划分为一个物体上的“从节点”和另一个物体上的“主段”。对每个从节点，脚本将其投影到最近的主段上，计算有符号的法向间隙，并累加罚函数贡献。

脚本采用 双向 处理（每个物体的界面节点相对于另一个物体都作为从节点），这给出了对称的表述，并避免了单主单从选择所带来的偏差。

用 L-BFGS 进行能量最小化

这一模式与 超弹性梁（新胡克） 中介绍的相同：两个网格共用一个 u 参数，边界条件在闭包内部施加，由 autograd 提供梯度。由于间隙函数依赖于当前构型，接触项在每次迭代中都会重新计算：

列表 11 examples/solid/hertzian_contact/hertzian_contact.py（示意）
def closure():
    optimizer.zero_grad()
    u_active = apply_bcs(u)
    E_elastic_1 = model_circle.energy(...)
    E_elastic_2 = model_block.energy(...)
    E_contact   = penalty_contact_energy(u_active, k_p)
    total = E_elastic_1 + E_elastic_2 + E_contact
    total.backward()
    return total

for it in range(25):
    optimizer.step(closure)

与赫兹解析解的对比

对于弹性圆柱与半空间之间的二维平面应变接触，赫兹理论给出抛物线形的压力分布：

\[p(x) \;=\; p_\text{max}\, \sqrt{1 - (x / a)^2}, \qquad p_\text{max} = \sqrt{\frac{P\, E^*}{\pi\, R^*}}, \qquad a = \sqrt{\frac{4 P\, R^*}{\pi\, E^*}},\]

其中等效模量 \(E^*\) 和等效半径 \(R^*\) 由两体几何导出。脚本从收敛解中提取接触压力，并将其叠加到 hertzian_contact.png 中的解析曲线上——只要接触斑附近的网格足够细，有限元结果就吻合得很好（chara_length=0.1 是一个合理的起点）。

赫兹接触的冯·米塞斯应力场

图 44 hertzian_contact.py 的输出：可变形压头（圆形）压向被位移的弹性体（右侧）时的冯·米塞斯应力场。应力集中在压头内部一个大致呈泪滴形的区域中，其中心略微位于接触斑内侧——这是经典的赫兹次表面应力分布。左侧固定壁和右侧弹性体上的规定位移箭头展示了边界设置。

运行示例

cd examples/solid/hertzian_contact
python hertzian_contact.py    # writes hertzian_contact.png

控制台输出会报告接触压力分布和最大位移。

下一步