1 导读

两只钢球相撞，“叮”的一声，前后不过几十微秒。在这几十微秒里，接触处被压出一个直径约半毫米的小圆，压强高达数 GPa，然后球又几乎毫无损失地弹开。一排紧挨着的钢球（牛顿摆）被一端的球撞击，另一端恰好只飞出一个球。把成百上千个球排成一列，敲击一端，会出现一种宽约 5 个球、形状不变、速度随振幅变化的“孤立波”。这三件事由同一个力学规律串起来：两个弹性曲面接触时，压力与压缩量之间是 \(F\propto\delta^{3/2}\) 的非线性关系。本讲从这个规律出发，一步一步走到颗粒晶体中的强非线性波。

学习目标

学完本讲，你应当能够：

1. 用量纲分析和标度论证得到 \(a\sim\sqrt{R\delta}\)、\(F\sim E^*\sqrt{R}\,\delta^{3/2}\)，并说明精确系数 \(4/3\) 的来源；

2. 由能量守恒推出碰撞的最大压缩量与接触时间 \(t_c\approx2.94\,\delta_m/v\propto v^{-1/5}\)，判断准静态近似与弹性假设何时成立；

3. 解释牛顿摆为什么近似“进一出一”，理解色散与非线性在其中的作用；

4. 从离散运动方程出发，经长波近似和行波代入，得到 Nesterenko 孤立波的 \(\cos^4\) 剖面、宽度与 \(V_s\propto F_m^{1/6}\propto v_m^{1/5}\) 的标度，并能用数值模拟检验；

5. 处理“电路通断法测接触时间”一类实验数据：线性化、最小二乘、不确定度和系统误差分析。

建议学时

共 2 学时（约 90 分钟）：前置模块与赫兹接触 25 分钟；碰撞与恢复系数 15 分钟；牛顿摆与色散 15 分钟；Nesterenko 孤立波与数值验证 25 分钟；应用 10 分钟。例题与实验题供课后练习，附录供深入阅读。

推导路线

胡克定律 + Boussinesq 点力解 \(\Rightarrow\) 赫兹接触律 \(F=\tfrac43E^*\sqrt{R^*}\delta^{3/2}\) \(\Rightarrow\)（能量守恒）碰撞时间 \(\Rightarrow\)（多球）牛顿摆 \(\Rightarrow\)（长链 + 长波近似）孤立波 \(\Rightarrow\) 界面与器件。

2 名词与符号

两个容易混淆的名称：孤立波（solitary wave）指形状不变、局域在有限区域内传播的行波；孤子（soliton）原指某些特殊方程（如描述浅水波的 KdV 方程）中、彼此碰撞后能完全复原的孤立波。颗粒链中的孤立波彼此碰撞后会留下少量“残波”，严格说不是孤子，本讲统一称“孤立波”。

3 前置模块：弹性体的基本描述

3.1 应力、应变与胡克定律

3.2 量纲分析与标度论证

标度论证比量纲分析多走一步：先根据物理图像判断哪些量“同一数量级”（例如“接触区内的应变约为压缩量除以接触半径”），再用这些粗略关系推出幂律。它给出指数和数量级，精确系数要靠完整求解。本讲反复使用这两种方法，并在每处用精确结果加以核对。

3.3 弹性半空间受点力：Boussinesq 解

由于方程是线性的，分布载荷产生的位移等于各个面元 \(\dd S'\) 上的力 \(p\,\dd S'\) 产生的位移之和（叠加原理）： \[\tag{3.5} u_z(\vect{r})=\frac{1-\nu^2}{\pi E}\iint\frac{p(\vect{r}')\,\dd S'}{\abs{\vect{r}-\vect{r}'}} .\] 这个积分就是附录 C 中求赫兹压强所产生位移的出发点。注意积分形式与“电荷面分布产生的静电势”完全相同（\(1/\abs{\vect{r}-\vect{r}'}\) 核），这个类比在附录中很有用。

3.4 曲面的局部几何：间隙函数

4 赫兹接触

4.1 问题的提法

两个弹性体在一点相切，再用力 \(F\) 把它们压紧。远离接触区的参考点（例如两球心）相互靠近 \(\delta\)，叫压缩量。我们要求 \(F(\delta)\)、接触圆半径 \(a\) 和接触面上的压强分布 \(p(r)\)。

图 1(a) 画出了几何关系。如果两物体可以互相“穿透”，未形变的轮廓会重叠，重叠深度正好是 \(\delta\)。真实情况是两物体在接触区都被压下去，表面重合。设物体 \(i\) 表面在 \(r\) 处被压下的法向位移为 \(u_i(r)\)。未形变时，两表面在 \(r\) 处的间距是 \(h(r)\)；两物体整体靠近 \(\delta\) 后，若不变形，间距变为 \(h(r)-\delta\)（负值表示穿透）；两表面再各被压下 \(u_1\)、\(u_2\)，间距变为 \(h-\delta+u_1+u_2\)。接触区内这个间距为零，即 \[\tag{4.1} u_1(r)+u_2(r)=\delta-h(r)=\delta-\frac{r^2}{2R^*},\qquad r\le a ,\] 接触区外则要求两表面不重叠：\(u_1+u_2>\delta-h(r)\)（\(r>a\)）。

4.2 标度论证

先不求解，估计各量的数量级。

4.3 精确解：赫兹压强分布

现在求系数。思路是“猜测—检验”：找一个压强分布 \(p(r)\)，使它在接触区内产生的位移恰好满足式 (4.1)，即为常数减去 \(r^2\) 的项，并满足 \(p\ge0\)、在边界 \(r=a\) 处 \(p\to0\)（表面光滑地分离）。

4.4 弹性势能与刚度

接触力只依赖于压缩量 \(\delta\)，且是弹性的（加载、卸载走同一条曲线），所以它是保守力，存在势能： \[\tag{4.11} U(\delta)=\int_0^{\delta}K\delta'^{3/2}\dd\delta'=\frac25K\delta^{5/2}=\frac{8}{15}E^*\sqrt{R^*}\,\delta^{5/2}=\frac25F\delta .\] 同样的压缩量与终态力下，它是线性弹簧储能 \(\frac12F\delta\) 的 \(4/5\)：赫兹“弹簧”一开始很软。刚度（切线斜率）为 \[\tag{4.12} k(\delta)=\dv{F}{\delta}=\frac32K\delta^{1/2}=2E^*a .\] 最后一个等号有简单的物理意义：刚度只由接触圆半径决定，与平底圆柱压头（半径 \(a\)）的 \(2E^*a\) 相同。

4.5 何时开始塑性屈服

赫兹理论假定纯弹性。由式 (4.10)，\(p_0\) 可达 GPa 量级，必须检查材料是否屈服。

接触区表面上三个方向都受压，切应力反而不大；最大切应力出现在表面下方。弹性力学给出对称轴上（\(r=0\)，深度 \(z\)）的主应力（引用结果[4]）： \[\frac{\sigma_z}{p_0}=-\frac{1}{1+z^2/a^2},\qquad \frac{\sigma_r}{p_0}=\frac{\sigma_\theta}{p_0}=-(1+\nu)\Big(1-\frac za\arctan\frac az\Big)+\frac{1}{2(1+z^2/a^2)} ,\] 负号表示压应力。于是 \(\tau_{\max}=\frac12\abs{\sigma_z-\sigma_r}\)。对 \(z\) 数值求极大（\(\nu=0.3\)）：\(\tau_{\max}\approx0.31p_0\)，位于深度 \(z\approx0.48a\) 处。因此屈服起始条件为 \[\tag{4.13} p_0\approx1.6\,Y .\] 塑性首先在表面下方约半个接触半径深处萌生，被周围的弹性材料包围，表面上看不出来。这就是轴承“次表面疲劳”剥落的原因[4]。低碳结构钢 \(Y\) 只有约 0.25 GPa，调质钢约 1 GPa，淬火轴承钢可达约 2 GPa。由上面的数值例，两只 10 mm 钢球静压 10 N 时 \(p_0\approx{1\,}\mathrm{GPa}\)，仍在弹性范围；碰撞时的压强则要高得多（见第 5.6 节）。

5 两球碰撞

5.1 相对运动方程

两球质量 \(m_1,m_2\)，沿连心线以相对速度 \(v\) 正碰。设球心位置为 \(x_1<x_2\)，压缩量 \(\delta=(x_1-x_2)+(R_1+R_2)\)（未接触时 \(\delta<0\)，接触时 \(\delta>0\)）。接触力 \(F(\delta)\) 沿连心线，大小相等方向相反： \[m_1\ddot x_1=-F(\delta),\qquad m_2\ddot x_2=+F(\delta).\] 第一式除以 \(m_1\)，第二式除以 \(m_2\)，相减： \[\tag{5.1} \ddot\delta=\ddot x_1-\ddot x_2=-\Big(\frac1{m_1}+\frac1{m_2}\Big)F \quad\Longrightarrow\quad m^*\ddot\delta=-K\delta^{3/2},\qquad m^*=\frac{m_1m_2}{m_1+m_2}.\] 这是一个质量为 \(m^*\)（约化质量）、受非线性弹力的单自由度振子，初始条件 \(\delta(0)=0\)、\(\dot\delta(0)=v\)。这里用到了准静态假设：每一时刻接触力由该时刻的 \(\delta\) 按静态赫兹公式给出，球体内部的弹性波已经“来得及”把力传开（第 5.5 节检验）。球落在固定的大质量平板上，取 \(m_2\to\infty\)，\(m^*=m_1\)，\(R^*=R_1\)。

5.2 最大压缩量

式 (5.1) 乘以 \(\dot\delta\) 并对时间积分，得能量守恒： \[\frac12m^*\dot\delta^2+\frac25K\delta^{5/2}=\frac12m^*v^2 .\] 最大压缩时 \(\dot\delta=0\)： \[\tag{5.2} \delta_m=\Big(\frac{5m^*v^2}{4K}\Big)^{2/5}=\Big(\frac{15\,m^*v^2}{16\,E^*\sqrt{R^*}}\Big)^{2/5}.\] 相应的最大力与最大接触半径为 \(F_m=K\delta_m^{3/2}\)、\(a_m=\sqrt{R^*\delta_m}\)。

5.3 接触时间

由能量守恒解出 \(\dot\delta=v\sqrt{1-(\delta/\delta_m)^{5/2}}\)。压缩阶段从 \(\delta=0\) 到 \(\delta_m\)，恢复阶段是它的时间反演（式 (5.1) 不含 \(\dot\delta\)，在 \(t\to-t\) 下不变），两段用时相等。于是 \[\tag{5.3} t_c=2\int_0^{\delta_m}\frac{\dd\delta}{\dot\delta}=\frac{2\delta_m}{v}\int_0^1\frac{\dd x}{\sqrt{1-x^{5/2}}},\qquad x=\delta/\delta_m .\] 附录 D 把这个积分化成 B 函数（\(B(x,y)\equiv\int_0^1t^{x-1}(1-t)^{y-1}\dd t\)，一种可查表的标准积分）：\(\int_0^1(1-x^{5/2})^{-1/2}\dd x=\frac25B(\frac25,\frac12)=1.4716\)。所以 \[\tag{5.4} \boxed{t_c=2.943\,\frac{\delta_m}{v}=2.943\Big(\frac{15\,m^*}{16\,E^*\sqrt{R^*}}\Big)^{2/5}v^{-1/5}}.\]

5.4 数值例：钢球

两只半径 \(R={10\,}\mathrm{mm}\) 的钢球（\(m={32.9\,}\mathrm{g}\)，\(m^*=m/2\)，\(R^*=R/2\)，\(E^*={115.4\,}\mathrm{GPa}\)）相撞，结果列于表 2。

几个值得注意的数：以 1 m/s 相撞，峰值力约 1 kN，是球重的约 3000 倍，而压缩量只有 20 μm。接触时间几十微秒，正是实验中用计数器或示波器能测的量级（第 10 节）。赫兹的例子（\(R={25\,}\mathrm{mm}\) 钢球，\(v={10\,}\mathrm{mm/s}\)）用式 (5.4) 算得 \(t_c={378\,}\mathrm{\mu s}\)、\(a_m={0.127\,}\mathrm{mm}\)，与原文的 380 μs、0.13 mm 一致；地球大小的钢球按 \(t_c\propto R\) 放大得 26.7 小时[1]。

5.5 准静态近似何时成立

式 (5.1) 假设接触力随 \(\delta\) 瞬时调整。实际上，接触区的应力要靠弹性波传遍整个球，球才作为整体减速。只有碰撞持续时间远大于弹性波穿过球所需的时间，才能把每一时刻看成静态平衡： \[\tag{5.5} t_c\gg\frac{2R}{c_0}.\]

准静态成立时，碰撞中转化为球内弹性振动和波的能量很少。Hunter（1957）与 Reed（1985）专门估算了这部分能量[22], [23]，对于 \(v\ll c_0\) 的碰撞它只占入射动能的很小一部分（Reed 指出 Hunter 的估计偏低约 4.5 倍，但结论的量级不变）。

5.6 恢复系数

理想赫兹碰撞完全弹性

式 (5.1) 中的力是保守力（只依赖于 \(\delta\)、单值），且方程具有时间反演对称性，所以分离时 \(\dot\delta=-v\)，\(e=1\)。这是在线弹性、准静态和无内耗三条假设下才成立的结论。真实碰撞中 \(e<1\)，原因有三类：

1. 塑性屈服。由式 (4.10)，\(p_0=\frac{2E^*}{\pi}(\delta_m/R^*)^{1/2}\)，而 \((\delta_m/R^*)^{5/2}=15m^*v^2/(16E^*R^{*3})\)。令 \(p_0=1.6Y\)（式 (4.13)），对球（半径 \(R\)、密度 \(\rho\)）撞刚性平板（\(m^*=\frac43\pi R^3\rho\)，\(R^*=R\)）解出屈服速度： \[\tag{5.6} \frac{\rho v_Y^2}{Y}=\frac{4}{5\pi}(0.8\pi)^5\Big(\frac{Y}{E^*}\Big)^4\approx26\Big(\frac{Y}{E^*}\Big)^4 .\] 与 Johnson 书中给出的系数 26 一致[4]。式中没有 \(R\)：屈服速度与球的大小无关。钢球撞钢板，\(Y={1\,}\mathrm{GPa}\) 时 \(v_Y\approx{0.14\,}\mathrm{m/s}\)，\(Y={2\,}\mathrm{GPa}\) 时约 0.77 m/s；两等球相撞时 \(m^*/R^{*3}\) 是球撞平板时的 4 倍，相对屈服速度减半。所以日常的钢球碰撞往往已经轻微超过初始屈服。不过塑性区很小，被弹性材料包围，能量损失有限，恢复系数仍接近 1，并随速度增大而缓慢减小[4]。

2. 粘弹性内耗。材料内部的分子摩擦使接触力多出一项与形变速率有关的耗散力，Brilliantov 等由准静态粘弹性理论推出了它的形式（\(\propto\sqrt\delta\,\dot\delta\)）[24]。Kuwabara 与 Kono 由此得到 \(1-e\propto v^{1/5}\)[25]：速度越大，\(e\) 越小，但依赖很弱。

3. 表面粗糙与振动激发。表面微凸体先接触、先屈服；球体振动模式会带走少量能量（第 5.5 节）。

6 牛顿摆：为什么只弹出一个球

6.1 守恒律不够

\(N\) 个质量都为 \(m\) 的球排成一列，第 1 个以 \(v_0\) 撞入，其余静止。设末态速度为 \(v_1,\dots,v_N\)。动量与能量守恒给出 \[\tag{6.1} \sum_{n=1}^Nv_n=v_0,\qquad \sum_{n=1}^Nv_n^2=v_0^2 .\] \(N\) 个未知数、两个方程。\(N=2\) 时，除平凡解（未碰撞）外只有“交换速度”一解。\(N\ge3\) 时有无穷多解。例如 \(N=3\)，\((v_1,v_2,v_3)=(0,0,v_0)\) 满足式 (6.1)，\(\big(-\tfrac13,\tfrac23,\tfrac23\big)v_0\) 也满足。所以“为什么只飞出一个”不是守恒律能回答的，答案取决于球之间的相互作用力律[6]。

6.2 有间隙的球：二体碰撞序列

如果相邻球之间有微小间隙 \(g\)，且 \(g\) 远大于碰撞的最大压缩量 \(\delta_m\)，那么第 1、2 球的碰撞在第 2 球碰到第 3 球之前就已结束。每次碰撞都是等质量的完全弹性二体碰撞，速度完全交换：\(1\to2\to\dots\to N\)，最后只有第 \(N\) 个球以 \(v_0\) 飞出，其余全部静止。这是教科书上“进一出一”的理想情形。由第 5 节，\(\delta_m\) 只有几到几十微米，所以几十微米的间隙就已足够。

6.3 紧贴的球：多体问题

真实的牛顿摆通常是紧贴的（零间隙、零预压）。这时第 1 个球还没停下，压缩已传到第 3、4 个球，多个接触同时受压，必须解 \(N\) 体运动方程 \[\tag{6.2} m\ddot u_n=F(u_{n-1}-u_n)-F(u_n-u_{n+1}),\qquad F(\psi)=\begin{cases}A\psi^{\beta}, & \psi>0\\ 0, & \psi\le0\end{cases}\] 其中 \(u_n\) 是第 \(n\) 个球离开初始位置的位移，\(\psi=u_{n-1}-u_n\) 是相邻两球的压缩量。\(F=0\)（\(\psi\le0\)）表示球只能互相推、不能互相拉。赫兹接触 \(\beta=3/2\)；作为对照，\(\beta=1\) 是“只压不拉”的线性弹簧。在压缩区内，\(\beta=1\) 的单边链与下面式 (6.3) 的双边线性链方程完全相同，区别只在于接触分离后不再有拉力。

图 3(a) 给出 5 球牛顿摆的数值结果。赫兹接触下末球获得 \(0.989v_0\)，第 4 球 \(0.127v_0\)，前三个球以 \(-0.071v_0\)、\(-0.030v_0\)、\(-0.014v_0\) 后退；动量和、动能和分别守恒到 \(10^{-10}\) 与 \(10^{-7}\)。末球带走 \(97.8\%\) 的能量，其余球速度很小，肉眼看来就是“只飞出一个”。线性弹簧链的末球只有 \(0.942v_0\)，第 4 球却有 \(0.296v_0\)，“多球飞出”明显得多。Hinch 与 Saint-Jean 对长链的计算给出末端两球 \(0.986V\)、\(0.149V\)，撞入球反弹 \(-0.07V\)[7]，与我们 5 球的结果定性一致（末两球分别为 \(0.989v_0\) 与 \(0.127v_0\)，撞入球 \(-0.071v_0\)）。

6.4 色散：线性链的对照

为什么线性链会“散开”？考虑线性弹簧（刚度 \(k\)，允许拉伸）连接的无限长等质量链，运动方程为 \[\tag{6.3} m\ddot u_n=k(u_{n+1}-2u_n+u_{n-1}).\] （第 41 届决赛考过双原子链，这里只用单原子链作为对照。）

由式 (6.4)，\(v_g=d\sqrt{k/m}\cos(qd/2)\)，短波（\(qd\to\pi\)）的群速度趋于零（图 3(c)）。任何脉冲都可以写成许多不同波数的正弦波之和（傅里叶分解），脉冲越窄，所需的短波成分越多。一次碰撞产生的脉冲很尖，短波成分很多，这些成分跑得慢，于是脉冲展宽、后面拖着振荡尾巴（图 3(b) 紫线）。对于只压不拉的单边线性链（橙线），压缩部分的传播与双边链相同，但尾巴中本应出现拉伸的地方接触分离了，振荡被截掉，只剩下展宽的主脉冲和后方的缓慢后退区；色散导致展宽这一结论在定性上不变。到达链末端时，能量分摊到好几个球上，就会有几个球一起飞出。Hinch 与 Saint-Jean 的计算表明，线性力律下飞出的球数约为 \(1.5N^{1/3}\)[7]。

6.5 赫兹链：声真空

赫兹链的情况根本不同。先看有静态预压力 \(F_0\) 的链：每个接触有静压缩 \(\psi_0=(F_0/A)^{2/3}\)。在它附近做小振动，等效线性刚度为切线斜率 \(k_0=\frac32A\psi_0^{1/2}\)（式 (4.12)），长波声速为 \[\tag{6.5} c_{\mathrm{lin}}=d\sqrt{\frac{k_0}{m}}=d\sqrt{\frac{3A}{2m}}\Big(\frac{F_0}{A}\Big)^{1/6}\propto F_0^{1/6}.\] 当 \(F_0\to0\) 时线性声速趋于零。无预压的赫兹链里不存在线性声波，任何扰动都是“强非线性”的。Nesterenko 把这种介质称为声真空[11]。

没有线性声速，就没有“所有长波都一样快”的基准。脉冲的传播速度取决于它自身的振幅：振幅大处更“硬”、跑得快，非线性有使波前变陡的趋势；离散性（色散）又使脉冲展宽。两者在某个固定宽度上达到平衡，形成一个形状不变的孤立波。这正是图 3(b) 中赫兹链脉冲始终紧凑的原因，也是 Herrmann 与 Schmälzle 所说的“无色散传递”的物理实质[6]。下一节定量地求出这个孤立波。

7 Nesterenko 孤立波

7.1 离散方程与标度对称性

考虑一列相同的弹性球（直径 \(d=2R\)，质量 \(m=\frac{\pi}{6}\rho d^3\)），初始时恰好相切、无预压。第 \(n\) 个球的位移为 \(u_n\)，第 \(n-1\) 与第 \(n\) 个球之间的接触压缩量为 \(\psi_n=u_{n-1}-u_n\)。由式 (4.9)，等径同材料球 \(R^*=R/2=d/4\)、\(E^*=E/[2(1-\nu^2)]\)，接触系数 \[\tag{7.1} A=\frac43E^*\sqrt{R^*}=\frac{E\sqrt{d}}{3(1-\nu^2)}=\frac{E\sqrt{2R}}{3(1-\nu^2)} .\] 运动方程为（括号下标 \(+\) 表示只取正值） \[\tag{7.2} m\ddot u_n=A\big[(\psi_n)_+^{3/2}-(\psi_{n+1})_+^{3/2}\big].\] 用 \(\psi\) 表示更简洁：\(\ddot\psi_n=\ddot u_{n-1}-\ddot u_n\)，代入两次式 (7.2)。在孤立波所在的区域，各接触都处于压缩状态（\(\psi>0\)），可以去掉“取正值”的记号： \[\tag{7.3} m\ddot\psi_n=A\big[\psi_{n-1}^{3/2}-2\psi_n^{3/2}+\psi_{n+1}^{3/2}\big].\] 右边是 \(f=\psi^{3/2}\) 的二阶差分，与线性链式 (6.3) 结构相同，只是对象从位移换成了“力”。

7.2 长波近似

定义压缩应变 \(w(x,t)=-\pdv{u}{x}>0\)（压缩为正）。附录 E 逐步完成如下三个展开：

1. 第 \(n\) 个球受的合力是两侧接触力之差，接触位于 \(x\mp d/2\)，所以 \(m\ddot u=f(x-\frac d2)-f(x+\frac d2)=-d\,f'-\frac{d^3}{24}f'''+\cdots\)；

2. 接触压缩量 \(\psi(x)=u(x-\frac d2)-u(x+\frac d2)=d\,w+\frac{d^3}{24}w''+\cdots\)；

3. 把 2 代入 \(f=\psi^{3/2}\) 并展开（要求 \(d^2w''/w\ll1\)），再代入 1，整理后对 \(x\) 求一次导数，得到应变场满足的方程。

结果是 Nesterenko 的长波方程[9], [11]： \[\tag{7.6} w_{tt}=c^2\Big[w^{3/2}+\frac{d^2}{10}\,w^{1/4}\big(w^{5/4}\big)_{xx}\Big]_{xx},\qquad c^2=\frac{Ad^{5/2}}{m}=\frac{2E}{\pi\rho(1-\nu^2)} .\] 方括号中第一项是连续介质的非线性“应力”（赫兹律），第二项来自离散性（色散）。如果把赫兹律换成线性律，第一项变成 \(w\)，就回到普通的波动方程加色散修正。

7.3 行波解：余弦四次方剖面

7.4 数值验证

我们用速度 Verlet 法（附录 F）积分 50 球链的式 (7.2)。第 1 球以 \(v_0\) 撞入（相当于一个与链中球相同的撞击器），能量相对误差 \(<10^{-7}\)。结果见图 4、5 与表 3。

8 应用：界面、冲击防护与声学器件

8.1 异质球链的界面

线性波在两种介质的界面上部分反射、部分透射，反射系数由声阻抗（密度与声速之积 \(Z=\rho c\)）之比决定，与振幅无关。强非线性的颗粒链则不同。我们模拟了等径（\(d={5\,}\mathrm{mm}\)）PTFE 球链与钢球链的界面（图 8）。PTFE 参数取 Daraio 等使用的动态值 \(E={1.46\,}\mathrm{GPa}\)、\(\nu=0.46\)、\(\rho={2200\,}\mathrm{kg/m^3}\)[13]。钢球质量是同径 PTFE 球的 3.6 倍，接触刚度约大两个数量级。

物理解释

两种情形都可以用“撞击器”的图像理解。(b) 中界面处最后一个钢球像一个重撞击器压向轻的 PTFE 链。重撞击器减速很慢，持续推挤，就会依次发出一串孤立波；它自身最终几乎静止，所以钢中没有反射压缩波。(a) 中轻的 PTFE 球撞向重的钢球，像轻球撞重球一样被弹回，形成反射的压缩波。这些结论与 Nesterenko、Daraio 等的实验和计算一致[15]。他们还发现，若对界面施加预压，从钢到 PTFE 的方向也会出现“反常”的反射压缩波，这一现象可用来设计非对称的声传输[15]。

8.2 冲击力的衰减与能量捕获

由图 8(b)，一个强脉冲进入软而轻的链后分裂成一串较弱的孤立波，能量在时间上被“摊开”，单个峰值力大大降低。Daraio 等把钢球段与 PTFE 段交替排列，发现软段之间会“捕获”能量，再以弱脉冲的形式缓慢放出，从而显著降低传到末端的冲击力[16]。这为冲击防护材料提供了一种不同于“吸能”（塑性、断裂）的思路：通过非线性动力学重新分配能量。

8.3 颗粒晶体声透镜与声二极管

声弹（非线性声透镜）

由式 (6.5) 与图 7，预压越大，波传得越快。Spadoni 与 Daraio 把 21 条平行的钢球链（每条 21 个直径 9.5 mm 的球）排成阵列，对各链施加不同的静态预压，使中间链的波慢、外侧链的波快。各链出口的孤立波在下游介质中同时到达一点，形成高度集中的“声弹”；改变预压即可调节焦点位置与强度[17]。这与光学透镜“中间厚、走得慢”的原理相同，只是“折射率”由预压控制。

声开关与声整流

Boechler 等在预压的颗粒链靠近一端处放一个轻的“缺陷”球，用特定频率驱动，当驱动振幅超过某个阈值时，体系的稳态响应突然改变性质（称为分岔），从单一频率的振动跳到含许多频率的宽带振动（准周期或混沌态）。这种跳变很陡，可用作声开关；正反两个方向的传输差别可达 \(10^4\) 倍以上，形成“声二极管”[18]。这些器件都利用了赫兹接触的非线性：小振幅下近似线性，大振幅下强非线性，以及只压不拉带来的不对称。

9 决赛风格例题

10 实验与数据处理：电路通断法测钢球接触时间

数据说明：表 4 是模拟数据。生成方法：由式 (5.4) 取 \(E^*={115.4\,}\mathrm{GPa}\) 算出 \(t_c\)，扣除阈值 \(\delta_{\mathrm{th}}={0.10\,}\mathrm{\mu m}\) 造成的漏计时间（按 \(2\int_0^{\delta_{\mathrm{th}}}\dd\delta/\dot\delta\) 精确计算），再乘以 \(1+\epsilon\)，\(\epsilon\) 为标准差 0.8% 的高斯随机数，每个角度 5 次。随机种子 20261003（生成程序随讲义源文件提供）。

11 小结：推导链

1. 弹性基础：胡克定律 + 泊松效应 \(\Rightarrow\) 半空间点力解 \(u_z=(1-\nu^2)F/(\pi Es)\)（式 (3.4)），由量纲分析确定形式；叠加原理给出分布载荷的位移积分。

2. 赫兹接触：几何间隙 \(r^2/2R^*\) + 接触条件 + 半椭球压强 \(\Rightarrow\) \(a=\sqrt{R^*\delta}\)、\(F=\frac43E^*\sqrt{R^*}\delta^{3/2}\)、\(p_0=3F/(2\pi a^2)\)、\(U=\frac25F\delta\)（式 (4.9)–(4.11)）。非线性来自接触面积随压缩增大。

3. 碰撞：约化质量 + 能量守恒 \(\Rightarrow\) \(\delta_m\propto v^{4/5}\)、\(t_c=2.943\delta_m/v\propto v^{-1/5}\)（式 (5.4)）。准静态判据 \(t_c/(2R/c_0)\propto(c_0/v)^{1/5}\gg1\)。理想赫兹碰撞 \(e=1\)；屈服速度与尺寸无关（式 (5.6)）。

4. 牛顿摆：守恒律不足以决定末态 \(\Rightarrow\) 需要力律。线性链有色散（式 (6.4)），脉冲散开；赫兹链无线性声速（式 (6.5)，声真空），脉冲保持紧凑。

5. 孤立波：离散方程的标度对称性（式 (7.4)）\(\Rightarrow\) 精确指数 \(V_s\propto F_m^{1/6}\propto v_m^{1/5}\)、宽度与振幅无关；长波近似（式 (7.6)）+ 行波 + 幂函数代换 \(\Rightarrow\) \(\cos^4\) 剖面、宽度约 \(5d\)、前因子 0.6802（式 (7.10)，对应 \(C=2/\sqrt5=0.894\)），离散数值解给出 0.689（\(C=0.906\)）。

6. 应用：撞击器图像解释界面的反射与分裂；预压调速 \(\Rightarrow\) 声透镜；非线性与分岔 \(\Rightarrow\) 声开关与整流。

参考文献

附录 A 弹性半空间受点力：Boussinesq 解的推导

本附录推导第 3.3 节的式 (3.4)，即表面位移 \(u_z=\frac{(1-\nu^2)F}{\pi Es}\)。路线是：

1. 写出三维弹性体的平衡方程（Navier 方程）；

2. 用一个“把弹性问题化为找调和函数”的表示把它化简；

3. 仿照静电学中的点电荷与像电荷，选出两个势函数；

4. 用表面无外力与合力为 \(F\) 两个条件定出待定常数。

所有中间式都用 sympy 逐一核验过。注意：本附录中 \(\rho=\sqrt{r^2+z^2}\) 表示到受力点的距离，\(\vect\psi\)、\(\beta\) 是势函数，与正文中的密度 \(\rho\)、接触压缩量 \(\psi\) 无关。

A.1 三维弹性力学的基本方程

A.2 问题的提法

柱坐标下的轴对称形变

问题绕 \(z\) 轴对称，位移只有 \(u_r(r,z)\)、\(u_z(r,z)\) 两个分量。各应变分量为 \[\varepsilon_{rr}=\pdv{u_r}{r},\quad \varepsilon_{\theta\theta}=\frac{u_r}{r},\quad \varepsilon_{zz}=\pdv{u_z}{z},\quad \varepsilon_{rz}=\frac12\Big(\pdv{u_r}{z}+\pdv{u_z}{r}\Big).\] 其中 \(\varepsilon_{\theta\theta}=u_r/r\) 来自周向：半径 \(r\) 的一圈质点整体向外移动 \(u_r\)，周长由 \(2\pi r\) 变为 \(2\pi(r+u_r)\)，相对伸长为 \(u_r/r\)。在任一点取局部的 \((r,\theta,z)\) 正交方向，式 (A.2) 照样成立（各向同性材料不分方向），于是 \(\sigma_{zz}=2G\varepsilon_{zz}+\lambda\operatorname{tr}\varepsilon\)，\(\sigma_{rz}=2G\varepsilon_{rz}\)，\(\operatorname{tr}\varepsilon=\varepsilon_{rr}+\varepsilon_{\theta\theta}+\varepsilon_{zz}\)。

A.3 Papkovich–Neuber 表示：化为找调和函数

满足 \(\nabla^2f=0\) 的函数叫调和函数，静电学中无电荷区域的电势就是调和函数。下面的表示把求解 Navier 方程化为寻找调和函数。

线弹性问题的解是唯一的：满足方程和全部边界条件的位移场只有一个。所以只要找到满足所有条件的 \(\vect\psi\)、\(\beta\)，得到的就是答案。这与静电学中“猜”出满足边界条件的电势（如像电荷法）是同一个思路。

A.4 选择势函数

第 1 步：点源。力集中作用在原点，位移在原点附近应当发散。仿照点电荷的电势，取 \[\vect\psi=\frac{A}{\rho}\,\uvec z,\qquad \rho=\sqrt{r^2+z^2}=\sqrt{x^2+y^2+z^2},\] \(A\) 为待定常数。\(1/\rho\) 是点电荷的电势，除原点外处处调和。

第 2 步：“像”源。只用第 1 步的 \(\vect\psi\) 时，表面上会留下切应力（下面第 A.5 节会算出），需要再加一个奇点在物体外面的调和函数来抵消。取 \[\beta=B\ln(\rho+z),\] \(B\) 为待定常数。它的电学图像是：沿负 \(z\) 轴（\(z<0\)）从原点到无穷远均匀分布的线电荷，其电势为 \(\int_0^L\frac{\dd s}{\sqrt{r^2+(z+s)^2}}=\ln(2L)-\ln(\rho+z)+\order{1/L}\)，常数 \(\ln(2L)\) 不影响导数。线电荷由点电荷叠加而成，每个点电荷的电势都调和，所以 \(\ln(\rho+z)\) 在线电荷以外处处调和（也可以直接求导验证）。它的奇点 \(\rho+z=0\) 只在负 \(z\) 轴上，即物体外部（加上原点），不会在物体内部引入新的力源。

第 3 步：位移。此时 \(\vect x\cdot\vect\psi=zA/\rho\)，记 \(\phi=zA/\rho+B\ln(\rho+z)\)。轴对称时梯度为 \(\nabla\phi=\pdv{\phi}{r}\uvec r+\pdv{\phi}{z}\uvec z\)。用 \(\pdv{\rho}{r}=\frac r\rho\)、\(\pdv{\rho}{z}=\frac z\rho\) 以及 \(\pdv{}{z}\ln(\rho+z)=\frac{z/\rho+1}{\rho+z}=\frac1\rho\)，由式 (A.6) 得 \begin{align*} u_r&=-\pdv{\phi}{r}=A\frac{rz}{\rho^3}-B\frac{r}{\rho(\rho+z)},\tag{A.7}\\ u_z&=\frac{4(1-\nu)A}{\rho}-\pdv{\phi}{z}=A\frac{z^2}{\rho^3}+\frac{(3-4\nu)A-B}{\rho}.\tag{A.8} \end{align*}

A.5 用边界条件定常数

第 1 步：表面无切应力定 \(B\)。在 \(z=0\)（\(r>0\)）处对式 (A.7)、(A.8) 求导。\(\pdv{}{z}\big[\rho(\rho+z)\big]_{z=0}=r\)，所以 \(\pdv{}{z}\frac{1}{\rho(\rho+z)}\Big|_{z=0}=-\frac{1}{r^3}\)，于是 \[\pdv{u_r}{z}\Big|_{z=0}=\frac{A+B}{r^2},\qquad \pdv{u_z}{r}\Big|_{z=0}=-\frac{(3-4\nu)A-B}{r^2},\] \[\sigma_{rz}(r,0)=G\Big(\pdv{u_r}{z}+\pdv{u_z}{r}\Big)_{z=0}=\frac{G}{r^2}\big[2B-2(1-2\nu)A\big]=0\quad\Longrightarrow\quad B=(1-2\nu)A .\] 如果没有 \(\beta\)（\(B=0\)），表面切应力为 \(-2(1-2\nu)GA/r^2\ne0\)，这就是第 A.4 节第 2 步要加“像源”的原因。

第 2 步：求 \(\sigma_{zz}\)。代入 \(B=(1-2\nu)A\)，\(u_z=A\big[\frac{z^2}{\rho^3}+\frac{2(1-\nu)}{\rho}\big]\)。求导得 \[\varepsilon_{zz}=\pdv{u_z}{z}=A\Big[\frac{2\nu z}{\rho^3}-\frac{3z^3}{\rho^5}\Big].\] 由第 A.3 节第 3 步，\(\operatorname{tr}\varepsilon=\nabla\cdot\vect u=2(1-2\nu)\pdv{\psi_z}{z}=-2(1-2\nu)\frac{Az}{\rho^3}\)。于是 \[\sigma_{zz}=2G\varepsilon_{zz}+\lambda\operatorname{tr}\varepsilon=2GA\Big[\frac{2\nu z}{\rho^3}-\frac{3z^3}{\rho^5}\Big]-\frac{2G\nu}{1-2\nu}\cdot\frac{2(1-2\nu)Az}{\rho^3}=-\frac{6GAz^3}{\rho^5}.\] \(\nu\) 的项恰好抵消。在 \(z=0\)、\(r\ne0\) 处 \(\sigma_{zz}=0\)，表面法向无外力的条件自动满足。同理可得 \(\sigma_{rz}=-6GA\,rz^2/\rho^5\)，在表面也为零。

第 3 步：合力定 \(A\)。代入式 (A.5)，令 \(t=r^2+c^2\)，\(\dd t=2r\,\dd r\)： \[\int_0^\infty-\frac{6GAc^3}{(r^2+c^2)^{5/2}}\,2\pi r\,\dd r=-6\pi GAc^3\int_{c^2}^\infty t^{-5/2}\dd t=-6\pi GAc^3\cdot\frac{2}{3c^3}=-4\pi GA .\] 结果与 \(c\) 无关，正是整体平衡所要求的（任一平面传递的合力都等于 \(F\)）。令其等于 \(-F\)： \[A=\frac{F}{4\pi G}.\]

A.6 结果与检验

附录 B 间隙函数的推导

本附录补全第 3.4 节的几何推导：先约定坐标与符号，再给出球面的精确式与误差估计；接着处理一般光滑曲面和主轴不重合的两个物体；最后说明间隙怎样进入接触条件。

B.1 坐标与符号约定

1. 公共切平面。两个未形变的光滑物体在点 O 相切。O 处两个表面各有一个切平面；两表面相切而不相交，这两个切平面必须重合，称为公共切平面。取 O 为原点，\(x\)、\(y\) 轴在公共切平面内，\(z\) 轴沿公共法线、指向物体 1（上方的物体）。

2. 两个高度函数。物体 1 的表面在 \((x,y)\) 处离开切平面的高度记为 \(z_1(x,y)\)，向上为正；物体 2 的表面离开切平面的深度记为 \(z_2(x,y)\)，向下为正。两种约定都使“表面离开切平面、远离对方”时取正值。于是两未形变表面在 \((x,y)\) 处的竖直间距为 \[\tag{B.1} h(x,y)=z_1(x,y)+z_2(x,y).\]

3. 凹面的符号。若物体 2 是一个凹窝（表面向上弯、包住物体 1），离开 O 后它的表面高于切平面，按约定 \(z_2<0\)。下文的曲率半径由 \(z=\frac{r^2}{2R}\) 定义，所以凹面的曲率半径取负值。

4. 为什么没有一次项。把 \(z_1\) 在 O 点作二元泰勒展开： \[z_1(x,y)=z_1(0,0)+\pdv{z_1}{x}\Big|_O x+\pdv{z_1}{y}\Big|_O y+\frac12\Big(\frac{\partial^2z_1}{\partial x^2}x^2+2\frac{\partial^2z_1}{\partial x\partial y}xy+\frac{\partial^2z_1}{\partial y^2}y^2\Big)_O+\cdots\] O 在表面上，所以 \(z_1(0,0)=0\)；切平面就是在 O 点与表面坡度相同的平面，而它就是 \(z=0\) 平面，所以两个一阶导数（表面在 O 点沿 \(x\)、\(y\) 方向的坡度）都为零。因此展开从二次项开始。\(z_2\) 同理。

B.2 球面：精确式、展开与误差

半径 \(R\) 的球面（以 O 为最低点）上，离对称轴距离 \(r=\sqrt{x^2+y^2}\) 处的高度为（勾股定理） \[\tag{B.2} z=R-\sqrt{R^2-r^2}.\] 第 1 步：有理化。分子分母同乘 \(R+\sqrt{R^2-r^2}\)： \[z=\frac{R^2-(R^2-r^2)}{R+\sqrt{R^2-r^2}}=\frac{r^2}{R+\sqrt{R^2-r^2}} .\] 第 2 步：提出抛物线近似。记 \(x\equiv r^2/R^2\)，则 \[\tag{B.3} z=\frac{r^2}{2R}\cdot\frac{2}{1+\sqrt{1-x}} .\] \(x\to0\) 时第二个因子趋于 1，所以领头项是 \(z\approx r^2/(2R)\)。

第 3 步：级数。由二项式展开 \(\sqrt{1-x}=\sum_{k\ge0}\binom{1/2}{k}(-x)^k=1-\frac x2-\frac{x^2}{8}-\frac{x^3}{16}-\cdots\)，其中 \(k\ge1\) 的通项为 \(-\frac{(2k)!}{(2k-1)\,4^k(k!)^2}x^k\)。代入式 (B.2)： \[z=R\big(1-\sqrt{1-x}\big)=R\Big(\frac x2+\frac{x^2}{8}+\frac{x^3}{16}+\cdots\Big)=\frac{r^2}{2R}+\frac{r^4}{8R^3}+\frac{r^6}{16R^5}+\cdots\] 第 4 步：丢掉高阶项的误差。定义相对误差 \(\varepsilon\equiv z\big/\frac{r^2}{2R}-1\)。由式 (B.3)，记 \(s=\sqrt{1-x}\)， \[\varepsilon=\frac{2}{1+s}-1=\frac{1-s}{1+s}=\frac{(1-s)(1+s)}{(1+s)^2}=\frac{x}{(1+\sqrt{1-x})^2}.\] 这是精确式。因为 \(1-x\le\sqrt{1-x}\le1\)（\(0\le x\le1\)），所以 \[\tag{B.4} \frac x4\le\varepsilon\le\frac{x}{(2-x)^2},\qquad \varepsilon\approx\frac{x}{4}=\frac{r^2}{4R^2}\quad(x\ll1).\] 接触区内 \(r\le a\)，所以丢掉 \(r^4\) 及更高次项造成的相对误差不超过约 \((a/R)^2/4\)。用正文的数值：两只 \(R={10\,}\mathrm{mm}\) 钢球静压 10 N，\(a/R^*=0.014\)，\(\varepsilon\approx5\times10^{-5}\)；以 1 m/s 相撞，\(a_m/R^*=0.064\)，\(\varepsilon\approx1.0\times10^{-3}\)；即使以 3 m/s 相撞，\(a_m/R^*=0.099\)，\(\varepsilon\approx2.5\times10^{-3}\)。（这里保守地用 \(R^*\) 估计；对每个球应当用它自己的 \(R_i\ge R^*\)，误差更小。）而接触区的应变本身就是 \(a/R\) 量级，线弹性理论忽略的非线性修正也是这个量级，比 \((a/R)^2/4\) 大得多。所以在赫兹理论的精度内，只需保留二次项。

B.3 一般光滑曲面：二次型与主曲率

第 1 步：二次型。由第 B.1 节第 4 条，任何光滑表面在切点附近都可以写成 \[\tag{B.5} z=\frac12\big(\kappa_{xx}x^2+2\kappa_{xy}xy+\kappa_{yy}y^2\big)+\order{r^3},\] 其中 \(\kappa_{xx}=\partial^2z/\partial x^2\) 等是 O 点的二阶导数（量纲为 1/长度）。球面 \(\kappa_{xx}=\kappa_{yy}=1/R\)，\(\kappa_{xy}=0\)。

第 2 步：转动坐标轴。把坐标轴绕 \(z\) 轴转过 \(\theta\)，新坐标 \((x',y')\) 与旧坐标的关系为 \(x=x'\cos\theta-y'\sin\theta\)，\(y=x'\sin\theta+y'\cos\theta\)。代入式 (B.5) 并按 \(x'^2\)、\(x'y'\)、\(y'^2\) 合并： \begin{align*} \kappa'_{xx}&=\kappa_{xx}\cos^2\theta+2\kappa_{xy}\sin\theta\cos\theta+\kappa_{yy}\sin^2\theta,\\ \kappa'_{yy}&=\kappa_{xx}\sin^2\theta-2\kappa_{xy}\sin\theta\cos\theta+\kappa_{yy}\cos^2\theta,\\ \kappa'_{xy}&=-\frac{\kappa_{xx}-\kappa_{yy}}{2}\sin2\theta+\kappa_{xy}\cos2\theta . \end{align*} 第 3 步：两个组合量。前两式相加得 \(\kappa'_{xx}+\kappa'_{yy}=\kappa_{xx}+\kappa_{yy}\)：平均曲率 \(S\equiv\frac12(\kappa_{xx}+\kappa_{yy})\) 与转角无关。前两式相减再用倍角公式得 \(\frac12(\kappa'_{xx}-\kappa'_{yy})=\frac{\kappa_{xx}-\kappa_{yy}}2\cos2\theta+\kappa_{xy}\sin2\theta\)。把它与 \(\kappa'_{xy}\) 合成一个复数： \[\tag{B.6} D\equiv\frac{\kappa_{xx}-\kappa_{yy}}2+\ii\,\kappa_{xy}\qquad\Longrightarrow\qquad D'=\ee^{-2\ii\theta}D .\] 即坐标轴转过 \(\theta\) 时，\(D\) 只是在复平面上转过 \(-2\theta\)，模 \(\abs D\) 不变。

第 4 步：消去交叉项。交叉项系数 \(\kappa'_{xy}=\Im D'\)。选 \(2\theta\) 等于 \(D\) 的辐角，\(D'\) 就成为正实数，交叉项消失： \[\tag{B.7} \tan2\theta=\frac{2\kappa_{xy}}{\kappa_{xx}-\kappa_{yy}} .\] 此时 \(\kappa'_{xx}=S+\abs D\)，\(\kappa'_{yy}=S-\abs D\)，即两个主曲率 \[\tag{B.8} \kappa_\pm=\frac{\kappa_{xx}+\kappa_{yy}}2\pm\sqrt{\Big(\frac{\kappa_{xx}-\kappa_{yy}}2\Big)^2+\kappa_{xy}^2},\qquad R'=\frac1{\kappa_+},\quad R''=\frac1{\kappa_-},\] \(R'\)、\(R''\) 叫主曲率半径，对应的方向叫主轴。在主轴坐标下，\(z=\frac12\big(x'^2/R'+y'^2/R''\big)\)。这就是“对称矩阵 \(\begin{pmatrix}\kappa_{xx}&\kappa_{xy}\\ \kappa_{xy}&\kappa_{yy}\end{pmatrix}\) 的本征值与本征方向”，这里只用到了三角恒等式和复数。

B.4 两个物体：主轴斜交

设物体 1 的主轴沿 \(x\)、\(y\)，主曲率半径 \(R_1'\)、\(R_1''\)；物体 2 的主曲率半径 \(R_2'\)、\(R_2''\)，其主轴相对 \(x\) 轴转过 \(\omega\)。

第 1 步：同一坐标系中相加。两个高度函数都是 \((x,y)\) 的二次型，由式 (B.1)，间隙的系数是两者系数之和：\(K_{xx}=\kappa^{(1)}_{xx}+\kappa^{(2)}_{xx}\) 等。相应地，平均曲率与复数 \(D\) 也直接相加。

第 2 步：写出各自的 \(S\)、\(D\)。物体 1 在自己的主轴坐标中 \(\kappa_{xy}=0\)： \[S_1=\frac12\Big(\frac1{R_1'}+\frac1{R_1''}\Big),\qquad D_1=\frac12\Big(\frac1{R_1'}-\frac1{R_1''}\Big).\] 物体 2 在自己的主轴坐标中 \(D_2^{(\text{自})}=\frac12\big(\frac1{R_2'}-\frac1{R_2''}\big)\) 为实数。它的主轴坐标相对 \((x,y)\) 转过 \(\omega\)，由式 (B.6) 反推，在 \((x,y)\) 坐标中 \(D_2=\ee^{2\ii\omega}D_2^{(\text{自})}\)；而 \(S_2=\frac12\big(\frac1{R_2'}+\frac1{R_2''}\big)\) 与坐标无关。

第 3 步：合成。\(S=S_1+S_2\)，\(D=D_1+\ee^{2\ii\omega}D_2^{(\text{自})}\)。由余弦定理（两复数之和的模）， \[\abs D^2=D_1^2+\big(D_2^{(\text{自})}\big)^2+2D_1D_2^{(\text{自})}\cos2\omega .\] 第 4 步：转到合成主轴。按第 B.3 节第 4 步转到 \(D\) 的主轴 \((x',y')\)，间隙成为 \(h=\frac12(S-\abs D)x'^2+\frac12(S+\abs D)y'^2\equiv Ax'^2+By'^2\)（取 \(B\ge A\)）。于是 \begin{gather*} A+B=S=\frac12\Big(\frac1{R_1'}+\frac1{R_1''}+\frac1{R_2'}+\frac1{R_2''}\Big),\tag{B.9}\\ B-A=\abs D=\frac12\sqrt{\Big(\frac1{R_1'}-\frac1{R_1''}\Big)^2+\Big(\frac1{R_2'}-\frac1{R_2''}\Big)^2+2\Big(\frac1{R_1'}-\frac1{R_1''}\Big)\Big(\frac1{R_2'}-\frac1{R_2''}\Big)\cos2\omega} .\tag{B.10} \end{gather*} 这正是接触力学文献中常用的形式[4]。我们用 200 组随机曲率（含凹面）与矩阵本征值逐一比较，二者一致到 \(10^{-12}\)；sympy 也验证了 \((B-A)^2=(A+B)^2-4\det\) 这一恒等式。

第 5 步：等效半径。把 \(h=Ax'^2+By'^2\) 写成 \(h=\frac{x'^2}{2R_a}+\frac{y'^2}{2R_b}\)，\(R_a=\frac1{2A}\)，\(R_b=\frac1{2B}\)，取几何平均，定义等效半径 \[\tag{B.11} R_e=\sqrt{R_aR_b}=\frac{1}{2\sqrt{AB}} .\] \(A=B\) 时 \(R_e=1/(2A)\)，就是正文的 \(R^*\)。\(h=\) 常数的等间隙线是以 \(x'\)、\(y'\) 为主轴的椭圆（图 11(a)）。\(A\ne B\) 时接触区也是椭圆，压强分布是半椭球面，求解要用到椭圆积分[4]，超出本讲范围。本讲只处理 \(A=B\) 的轴对称情形。

B.5 特例

1. 两球：\(R_i'=R_i''=R_i\)，\(D_1=D_2=0\)，所以 \(A=B=\frac12\big(\frac1{R_1}+\frac1{R_2}\big)=\frac1{2R^*}\)，\(h=r^2/(2R^*)\)，回到正文式 (3.6)。

2. 两等半径圆柱正交：\(R_1'=R_2'=R\)，\(R_1''=R_2''=\infty\)，\(\omega=90^\circ\)。\(D_1=\frac1{2R}\)，\(D_2=\ee^{\ii\pi}\frac{1}{2R}=-\frac1{2R}\)，相加为零，所以 \(A=B=\frac{1}{2R}\)，\(h=r^2/(2R)\)：与半径 \(R\) 的球压在平面上完全相同。测量表面力的实验正是利用这一点，把两片云母贴在两个轴线互相垂直的圆柱面上，Tabor 与 Winterton 在论文中明确指出“这种接触与球–平面接触相似”[33]。一般夹角时，由式 (B.10)，\(B-A=\abs{\cos\omega}/R\)，只有 \(\omega=90^\circ\) 时才轴对称。

3. 平行圆柱：\(\omega=0\)，\(\abs D=S\)，所以 \(A=0\)，\(B=\frac12\big(\frac1{R_1}+\frac1{R_2}\big)\)。间隙与沿轴线的坐标无关，接触区是一条窄带（线接触），二维问题的力–压缩关系与 \(F\propto\delta^{3/2}\) 不同，赫兹点接触的结果不适用[4]。

4. 球在凹球窝中：\(R_2<0\)。由 \(\frac1{R^*}=\frac1{R_1}-\frac1{\abs{R_2}}\)，当 \(\abs{R_2}>R_1\) 时 \(R^*>R_1\)：球窝越贴合，\(R^*\) 越大，同样载荷下接触圆越大、压强越低。\(\abs{R_2}\to R_1\) 时 \(R^*\to\infty\)，间隙处处为零，“小接触区”假设失效。

B.6 间隙怎样进入接触条件

第 1 步：刚性趋近量。在两个物体内部、远离接触区（距离 \(\gg a\)）的地方各取一个参考点，例如两球心。压紧过程中两参考点沿 \(z\) 轴相互靠近的距离记为 \(\delta\)。若两物体是刚体，靠近 \(\delta\) 后，\((x,y)\) 处两表面的竖直间距由 \(h\) 变为 \(h(x,y)-\delta\)；在 \(h<\delta\) 的区域，这个量为负，表示两刚体会“互相穿透”。

第 2 步：加上弹性位移。实际上两物体都会形变。记物体 \(i\) 表面在 \((x,y)\) 处的法向位移为 \(u_{zi}\)，并规定指向物体 \(i\) 内部为正（即离开对方）。参考点离接触区很远，那里的弹性位移已衰减到可以忽略（由 Boussinesq 解，位移随距离按 \(1/s\) 衰减，第 3.3 节）。所以形变后的实际间距为 \[\tag{B.12} g(x,y)=h(x,y)-\delta+u_{z1}(x,y)+u_{z2}(x,y).\] 第 3 步：接触的两种情形。实际间距不能为负（不能穿透）：\(g\ge0\)。在接触区内两表面贴合，\(g=0\)；在接触区外，表面分开，\(g>0\)。代入式 (B.12)： \[u_{z1}+u_{z2}=\delta-h\quad(\text{接触区内}),\qquad u_{z1}+u_{z2}>\delta-h\quad(\text{接触区外}).\] 另外，接触区内压强 \(p\ge0\)（无粘附），接触区外 \(p=0\)，两种情形可以合写为 \(g\ge0\)、\(p\ge0\)、\(p\,g=0\)。

第 4 步：与正文衔接。对轴对称情形 \(h=r^2/(2R^*)\)，这就是正文式 (4.1) 及其后的不等式。剩下的问题，即什么样的 \(p(r)\) 使 \(u_{z1}+u_{z2}\) 恰好等于 \(\delta-r^2/(2R^*)\)，由第 4.3 节与附录 C 解决。

附录 C 赫兹压强产生的表面位移

本附录证明式 (4.6)。由式 (4.4)，需计算接触区内一点 B（到中心 O 的距离为 \(r<a\)）处的积分 \[\tag{C.1} I(r)=\iint_{\rho\le a}\frac{p(\rho)\,\dd S}{s},\qquad p(\rho)=\frac{p_0}{a}\sqrt{a^2-\rho^2},\] 其中 \(\rho\) 是面元到 O 的距离，\(s\) 是面元到 B 的距离。

第 1 步：以 B 为极点的极坐标

以 B 为原点、以 \(\overrightarrow{OB}\) 方向为极轴建立极坐标 \((s,\varphi)\)。面元 \(\dd S=s\,\dd s\,\dd\varphi\)，于是被积函数中的 \(1/s\) 被消去： \[I=\int_0^{2\pi}\dd\varphi\int_0^{s_1(\varphi)}p\,\dd s .\] 这就是选 B 为极点的好处：点力解的 \(1/s\) 奇点被面元的 \(s\) 抵消。

第 2 步：沿射线的压强

射线上距 B 为 \(s\) 的点 C，满足 \(\overrightarrow{OC}=\overrightarrow{OB}+\overrightarrow{BC}\)，所以 \(\rho^2=r^2+s^2+2rs\cos\varphi\)。记 \(\alpha^2=a^2-r^2>0\)，\(\beta=r\cos\varphi\)，则 \[a^2-\rho^2=\alpha^2-2\beta s-s^2=(\alpha^2+\beta^2)-(s+\beta)^2 .\] 射线在 \(a^2-\rho^2=0\) 处离开接触圆：\(s_1=-\beta+\sqrt{\alpha^2+\beta^2}\)。

第 3 步：沿射线积分

令 \(s+\beta=\sqrt{\alpha^2+\beta^2}\,\sin\theta\)。\(s=0\) 对应 \(\theta_0=\arcsin\big(\beta/\sqrt{\alpha^2+\beta^2}\big)\)，\(s=s_1\) 对应 \(\theta=\pi/2\)： \[J(\varphi)\equiv\int_0^{s_1}\sqrt{\alpha^2-2\beta s-s^2}\,\dd s=(\alpha^2+\beta^2)\int_{\theta_0}^{\pi/2}\cos^2\theta\,\dd\theta =\frac{\alpha^2+\beta^2}{2}\Big(\frac\pi2-\theta_0\Big)-\frac{\alpha\beta}{2},\] 最后一步用了 \(\sin\theta_0\cos\theta_0=\alpha\beta/(\alpha^2+\beta^2)\)。

第 4 步：方向配对

射线 \(\varphi\) 与反方向射线 \(\varphi+\pi\) 的 \(\beta\) 互为相反数，因而 \(\theta_0\) 也互为相反数。两者相加，含 \(\theta_0\) 的项和 \(\alpha\beta\) 项都抵消： \[J(\varphi)+J(\varphi+\pi)=\frac{\pi}{2}(\alpha^2+\beta^2)=\frac\pi2\big(a^2-r^2+r^2\cos^2\varphi\big).\]

第 5 步：对方向积分

\[\int_0^{2\pi}J\,\dd\varphi=\int_0^{\pi}\big[J(\varphi)+J(\varphi+\pi)\big]\dd\varphi=\frac\pi2\Big[\pi(a^2-r^2)+\frac{\pi r^2}{2}\Big]=\frac{\pi^2}{4}(2a^2-r^2).\] 于是 \(I=\frac{p_0}{a}\cdot\frac{\pi^2}{4}(2a^2-r^2)\)，代入式 (4.4)： \[\tag{C.2} u_1+u_2=\frac{1}{\pi E^*}I=\frac{\pi p_0}{4aE^*}(2a^2-r^2),\qquad r\le a ,\] 即式 (4.6)。

接触区外

对 \(r>a\) 的点，用同样方法（此时只有一部分方向的射线穿过接触圆）可得[4] \[\tag{C.3} u_1+u_2=\frac{p_0}{2aE^*}\Big[(2a^2-r^2)\arcsin\frac ar+a\sqrt{r^2-a^2}\Big],\qquad r>a .\] 在 \(r=a\) 处与式 (C.2) 连续，均为 \(\pi p_0a/(4E^*)\)。图 1(b) 的圆圈是对式 (C.1) 的直接数值积分（对方位角的积分可化为第一类完全椭圆积分这一标准特殊函数，再对径向用自适应数值求积），与式 (C.2)、(C.3) 在 \(0\le r\le3a\) 内的最大相对偏差为 \(4\times10^{-9}\)。

附录 D 接触时间积分与 B 函数

式 (5.3) 中的积分，令 \(t=x^{5/2}\)，即 \(x=t^{2/5}\)，\(\dd x=\frac25t^{-3/5}\dd t\)： \[\tag{D.1} \int_0^1\frac{\dd x}{\sqrt{1-x^{5/2}}}=\frac25\int_0^1t^{-3/5}(1-t)^{-1/2}\dd t=\frac25B\Big(\frac25,\frac12\Big)=\frac25\cdot\frac{\Gamma(\frac25)\Gamma(\frac12)}{\Gamma(\frac{9}{10})}.\] 数值：\(\Gamma(0.4)=2.21816\)，\(\Gamma(0.5)=1.77245\)，\(\Gamma(0.9)=1.06863\)，积分 \(=1.47164\)，所以 \(t_c=2\times1.47164\,\delta_m/v=2.9433\,\delta_m/v\)。用 30 位精度直接数值积分，结果一致。

附录 E 长波展开与行波解的细节

(1) 记号

球 \(n\) 的平衡位置 \(x_n=nd\)，位移 \(u_n(t)=u(x_n,t)\)。接触 \((n-1,n)\) 位于两球之间 \(x_n-d/2\) 处，其压缩量 \[\psi(x_n-\tfrac d2)=u(x_n-d)-u(x_n).\] 定义连续函数 \(\psi(x)=u(x-\frac d2)-u(x+\frac d2)\) 与接触力 \(f(x)=A\psi(x)^{3/2}\)。

(2) 力的差分

式 (7.2) 为 \(m\ddot u(x)=f(x-\frac d2)-f(x+\frac d2)\)。把两项在 \(x\) 处泰勒展开，偶数阶项相消： \[\tag{E.1} f(x-\tfrac d2)-f(x+\tfrac d2)=-2\Big[\frac d2f'+\frac{1}{3!}\Big(\frac d2\Big)^3f'''+\cdots\Big]=-d\,f'-\frac{d^3}{24}f'''+\order{d^5}.\]

(3) 压缩量的展开

同理，\(\psi=u(x-\frac d2)-u(x+\frac d2)=-d\,u'-\frac{d^3}{24}u'''+\cdots\)。用 \(w=-u'\)（\(u'''=-w''\)）： \[\tag{E.2} \psi=d\,w+\frac{d^3}{24}w''+\order{d^5}.\]

(4) 3/2 次幂的展开

写 \(\psi=dw\,(1+\epsilon)\)，\(\epsilon=\frac{d^2w''}{24w}\)。长波近似要求 \(\abs\epsilon\ll1\)。用 \((1+\epsilon)^{3/2}=1+\frac32\epsilon+\order{\epsilon^2}\)： \[\tag{E.3} f=A(dw)^{3/2}\Big(1+\frac{d^2w''}{16w}\Big)=Ad^{3/2}\Big[w^{3/2}+\frac{d^2}{16}w^{1/2}w''\Big],\]

略去的是 \(\order{\epsilon^2}\)，即相对主项 \(d^4\) 阶的项。

(5) 合并

把式 (E.3) 代入式 (E.1)，保留到 \(d^2\) 的相对阶（第二项中只需用主阶 \(f\approx Ad^{3/2}w^{3/2}\)）： \[m\ddot u=-Ad^{5/2}\,\partial_x\Big[w^{3/2}+\frac{d^2}{16}w^{1/2}w''+\frac{d^2}{24}\big(w^{3/2}\big)''\Big].\] 展开 \((w^{3/2})''=\frac32w^{1/2}w''+\frac34w^{-1/2}w'^2\)，方括号中的 \(d^2\) 项为 \(\frac{d^2}{8}w^{1/2}w''+\frac{d^2}{32}w^{-1/2}w'^2\)。另一方面 \(w^{1/4}(w^{5/4})''=\frac54w^{1/2}w''+\frac{5}{16}w^{-1/2}w'^2\)，乘以 \(d^2/10\) 后正好是 \(\frac{d^2}{8}w^{1/2}w''+\frac{d^2}{32}w^{-1/2}w'^2\)。所以方括号 \(=w^{3/2}+\frac{d^2}{10}w^{1/4}(w^{5/4})''\)（sympy 逐项核对）。最后对 \(x\) 求导并用 \(\ddot w=-\partial_x\ddot u\)，即得式 (7.6)，其中 \(c^2=Ad^{5/2}/m\)。代入式 (7.1) 与 \(m=\frac\pi6\rho d^3\) 得 \(c^2=2E/[\pi\rho(1-\nu^2)]\)。

(6) 以 \(\psi\) 为变量的另一种近似（近似 II）

也可以不经过 \(w\)，直接对式 (7.3) 的二阶差分展开：\(f_{n-1}-2f_n+f_{n+1}=d^2f''+\frac{d^4}{12}f''''+\cdots\)，得 \[\tag{E.4} m\psi_{tt}=Ad^2\Big[\psi^{3/2}+\frac{d^2}{12}\big(\psi^{3/2}\big)_{xx}\Big]_{xx}.\] 它与式 (7.6) 在 \((d/L)^2\) 阶上同样正确，区别只在于 \(\psi\) 与 \(dw\) 本身相差 \(\order{d^3}\)（式 (E.2)）。行波解（例 2 中 \(\beta=3/2\)）为 \(\psi=\psi_m\cos^4\big(z/\sqrt3d\big)\)，宽度 \(\sqrt3\pi d=5.44d\)，\(V_s^2=\frac56(Ad^2/m)\psi_m^{1/2}\)，即 \(C=\sqrt{5/6}=0.913\)。两种近似的差别（约 2%）反映了被略去的高阶项的大小。

(7) 积分常数与“紧支撑”

积分两次时令常数为零，用的是孤立波前方 \(w\to0\)。\(\cos^4\) 解在 \(\abs{\kappa z}=\pi/2\) 处与零解光滑衔接（\(w\)、\(w'\)、\(w''\)、\(w'''\) 均为零），所以分段函数是长波方程的合法解。这种“紧支撑”（\(w\) 只在有限区间内不为零）是长波近似的产物。原因可以从展开的有效性条件看出：第 (4) 段要求 \(\abs\epsilon=\abs{d^2w''/(24w)}\ll1\)。对 \(w=\xi_m\cos^4\theta\)（\(\theta=\kappa z\)），\(w''/w=\kappa^2(12\tan^2\theta-4)\)，于是 \[\epsilon(\theta)=\frac{\kappa^2d^2}{24}(12\tan^2\theta-4)=\frac{\tan^2\theta}{5}-\frac1{15}.\] 在中心 \(\epsilon=-0.067\)，\(\theta=\pi/4\) 处 \(\epsilon=0.13\)，到 \(\theta\approx1.16\)（约为半宽 \(\pi/2\) 的 74%）处 \(\epsilon=1\)，再往外发散。所以在剖面外侧约四分之一的区域里，展开已经失效，“恰好降到零”的边界并不可信；离散链的真实孤立波尾部以双指数方式衰减[26]，但衰减极快，在图 4(b) 上与 \(\cos^4\) 的差别已难以分辨。

附录 F 数值方法

积分器

所有球链模拟使用速度 Verlet 格式： \[v_{n}^{k+1/2}=v_n^k+\tfrac{h}{2}\ddot u_n(u^k),\quad u_n^{k+1}=u_n^k+h\,v_n^{k+1/2},\quad v_n^{k+1}=v_n^{k+1/2}+\tfrac h2\ddot u_n(u^{k+1}).\] 它是二阶精度的“辛”格式（symplectic，即离散后的演化仍保持相空间体积，像真实的哈密顿系统一样），长时间积分中能量误差有界、不会系统漂移。步长取特征时间 \(\tau=(m/A)^{2/5}v_0^{-1/5}\) 的 \(1/300\)–\(1/1000\)，各模拟的总能量相对误差 \(<10^{-5}\)（多数 \(<10^{-7}\)），总动量守恒到舍入误差。两球碰撞（图 2）用 SciPy 的 RK45 积分，相对容差 \(10^{-11}\)，并用事件检测确定分离时刻。

测量方法

孤立波速度：记录第 \(j_1=20\) 与 \(j_2=40\) 个接触力达到最大的时刻 \(t_1\)、\(t_2\)，\(V_s=(j_2-j_1)d/(t_2-t_1)\)。宽度：把第 30 个接触的 \(\psi(t)\) 用 \(z=V_s(t-t_{\mathrm{peak}})\) 换成空间坐标后取半高全宽。由于标度对称性，无量纲结果对所有振幅相同，7 个振幅的计算只用来演示和检验程序。

交互演示

网页演示用同样的格式（JavaScript 实现）。hertz-collision 用 RK4 积分（步长 \(t_c/4000\)），读数中数值与解析 \(t_c\) 一致到 \(10^{-4}\)；hertz-soliton 在全钢链中测得的 \(V_s\) 与式 (7.10) 之比约为 1.012，与表 3 一致。

附录 G 最小二乘公式的推导

设 \(y_i=\alpha+bx_i+\epsilon_i\)，\(\epsilon_i\) 相互独立、方差都为 \(\sigma^2\)。令 \(Q(\alpha,b)=\sum(y_i-\alpha-bx_i)^2\) 取极小： \[\pdv{Q}{\alpha}=-2\sum(y_i-\alpha-bx_i)=0\Rightarrow \alpha=\bar y-b\bar x;\qquad \pdv{Q}{b}=-2\sum x_i(y_i-\alpha-bx_i)=0 .\] 把 \(\alpha\) 代入第二式，利用 \(\sum x_i(y_i-\bar y)=S_{xy}\)、\(\sum x_i(x_i-\bar x)=S_{xx}\)，得 \(b=S_{xy}/S_{xx}\)。由于 \(b=\sum_i\frac{x_i-\bar x}{S_{xx}}y_i\) 是 \(y_i\) 的线性组合，独立变量的方差相加： \[\sigma_b^2=\sum_i\frac{(x_i-\bar x)^2}{S_{xx}^2}\sigma^2=\frac{\sigma^2}{S_{xx}} .\] 同理 \(\alpha=\bar y-b\bar x\)，\(\bar y\) 与 \(b\) 不相关（因为 \(\sum(x_i-\bar x)=0\)），\(\sigma_\alpha^2=\sigma^2(\frac1n+\frac{\bar x^2}{S_{xx}})\)。\(\sigma\) 未知时用残差估计 \(s^2=Q_{\min}/(n-2)\)：\(n\) 个数据中两个自由度被 \(\alpha\)、\(b\) 用掉，分母取 \(n-2\) 才使 \(s^2\) 的期望等于 \(\sigma^2\)。这就是式 (10.1)。由函数 \(g(C_t)\) 求不确定度时用一阶泰勒展开：\(u_g=\abs{g'(C_t)}u_{C_t}\)，如 \(E^*\propto C_t^{-5/2}\) 给出 \(u_{E^*}/E^*=\frac52u_{C_t}/C_t\)。

附录 H 名词解释