两个圆盘怎样互相推，一百个圆盘怎样分担载荷

先约定对象、单位和方向

研究对象是平面上的圆盘：一半直径为 $d$ ，另一半直径为 $1.4\,d$ 。两种大小混在一起，是为了让圆盘排不成整齐的晶格。长度以 $d$ 为单位，能量以接触刚度 $\varepsilon$ 为单位，于是力的单位是 $\varepsilon/d$ 。

取圆盘 $i$ 和 $j$ ，半径为 $R_i$ 、 $R_j$ ，圆心距离为 $r$ 。两盘刚好相切时 $r=\sigma_{ij}=R_i+R_j$ 。软圆盘可以稍微压进对方，重叠量记作 $\delta=\sigma_{ij}-r$ ，只有 $\delta>0$ 时才有力。

半径为

R_i

和

R_j

的两个软圆盘，圆心相距

r

。重叠量

\delta=R_i+R_j-r

画得夸大了；箭头是沿圆心连线、大小相等方向相反的推力。原创示意图。

来源：C. S. O’Hern

重叠越多，推得越狠

O’Hern 等人用的接触势能只在重叠时不为零。指数 $\alpha=2$ 时它就是弹簧的势能； $\alpha=5/2$ 时是 Hertz 接触，三维弹性球就是这样，二维圆盘沿用它只是一种建模选择。

对 $r$ 求导再取负号，得到沿圆心连线、把两盘推开的力。谐振情形下 $f=\varepsilon\delta/\sigma_{ij}^2$ ，刚度是常数 $\varepsilon/\sigma_{ij}^2$ ；Hertz 情形下力按 $\delta^{3/2}$ 增长，越压越硬。

没有摩擦，所以接触力只有沿圆心连线的这一个分量。这一点在第五章加入摩擦之前一直成立，它决定了第二章里要数的未知量有多少。

V(r)=\frac{\varepsilon}{\alpha}\left(1-\frac{r}{\sigma_{ij}}\right)^{\alpha},\qquad f=-\frac{dV}{dr}=\frac{\varepsilon}{\sigma_{ij}}\left(\frac{\delta}{\sigma_{ij}}\right)^{\alpha-1}\quad(\delta>0)

O’Hern 等人式 (1)，Goodrich 等人用的是同一势能。

\delta\le0

时势能和力都为零：圆盘不相碰就互不影响。

来源：C. S. O’Hern · C. P. Goodrich

壁像一面镜子

实验室里的盒子有四面平直、无摩擦的壁。圆心到壁的距离 $h$ 小于半径 $R$ 时，壁对圆盘的作用就像圆盘碰到了自己的镜像：势能为 $(\varepsilon/\alpha)(1-h/R)^{\alpha}$ ，谐振情形下刚度是 $\varepsilon/R^2$ 。

这条壁的定律是本课程自己的选择。O’Hern 等人用周期性盒子，没有壁；我们在第四章读他们的数据时，会用周期性盒子重新算一遍。

完整例题：夹在两壁之间的两个圆盘

把一个小盘（ $R_1=0.5$ ）和一个大盘（ $R_2=0.7$ ）放在一条水平线上，左右各一面壁。两盘的自然长度合计 $2R_1+2R_2=2.4$ ，把两壁的间距压到 $2.39$ ，总共要吸收 $\Delta=0.01$ 的重叠。取谐振接触。

先建模。每个圆盘只受左右两个力，要平衡，这两个力必须大小相等。三个接触传递的是同一个力 $f$ ，这正是串联弹簧。三处的刚度分别是 $1/R_1^2=4$ 、 $1/\sigma^2=1/1.44$ 和 $1/R_2^2=1/0.49$ ，单位 $\varepsilon/d^2$ 。

再计算。三处重叠相加等于 $\Delta$ ，解出 $f=0.01/2.18=4.59\times10^{-3}\,\varepsilon/d$ 。三处重叠依次是 $1.15\times10^{-3}$ 、 $6.61\times10^{-3}$ 和 $2.25\times10^{-3}$ ，相加正好是 $0.01$ 。两盘之间的接触最软，吸收了约三分之二的压缩。

再和数值解比。我们的程序用 FIRE 弛豫求出同一构型的力，与这个闭式解的相对误差不超过 $2.2\times10^{-10}$ ，三种重叠、两种接触定律都检查过。

\Delta=\delta_{1w}+\delta_{12}+\delta_{2w}=\frac{f}{\varepsilon}\left(R_1^2+\sigma^2+R_2^2\right)\quad\Longrightarrow\quad f=\frac{\varepsilon\,\Delta}{R_1^2+\sigma^2+R_2^2}

串联弹簧的柔度相加。这一结果只对谐振接触成立；Hertz 接触的力不与重叠成正比，要用数值方法解。

来源：C. S. O’Hern · J. Guénolé

一百多个圆盘：力沿着少数路径走

把 128 个圆盘装进盒子挤紧，每个接触都按上面的定律给出一个力，但这些力大小悬殊。实验室默认循环挤到最紧时（堆积分数 0.86），有 303 个承力接触：约 25% 的接触上力不到平均值一半，约 6% 超过平均值两倍，最强的约为平均值的 3.1 倍。总的接触力有一半由 29% 的接触承担。

把每个接触画成一条线，线宽正比于力，强接触会连成穿过整个盒子的链条，这就是力链。O’Hern 等人图 20 用同样的画法画出二维双分散堆积的力网络。

在实验里，力链可以直接看到。光弹材料受力后会改变偏振光的颜色和亮度，Zhao 等人把光弹圆盘放在正交偏振光下拍照，从亮度的空间变化读出每一帧的压强。

来源：C. S. O’Hern · Y. Zhao

到实验里试一试

在堆积里找力链

先预测：挤到最紧时，力直方图会集中在平均值附近，还是拖着一条长尾？
播放默认预设，等盒子挤到最紧时暂停。看舞台上最粗的几条线连成什么形状，再看舞台下方的力直方图。
把接触定律从谐振换成 Hertz，重复一次。比较两种定律下直方图的宽度，并说明你用什么量来比较。

和一维颗粒链对照

打开一维 Hertz 颗粒链，看冲击怎样沿唯一的一条路径传到墙上。
回到二维堆积，数一数一个承力圆盘平均有几个接触。说明为什么二维堆积里的力可以绕路，而一维链不行。

练习与解答

两个同样大小的圆盘（ $R=0.5$ ）夹在两壁之间，总重叠 $\Delta=0.01$ ，谐振接触。求力和三处重叠。
展开解答
柔度之和为 $0.25+1+0.25=1.5$ ，所以 $f=0.01/1.5=6.67\times10^{-3}\,\varepsilon/d$ 。壁处重叠各为 $fR^2=1.67\times10^{-3}$ ，两盘之间为 $f\sigma^2=6.67\times10^{-3}$ ，合计 0.01。
对 Hertz 接触，把重叠加倍，力变成几倍？刚度（力对重叠的导数）变成几倍？谐振接触呢？
展开解答
Hertz 接触 $f\propto\delta^{3/2}$ ，重叠加倍时力变为 $2^{3/2}\approx2.83$ 倍，刚度 $\propto\delta^{1/2}$ ，变为约 1.41 倍。谐振接触的力加倍，刚度不变。
一个无摩擦圆盘只有两个接触，它能承受载荷吗？至少要几个接触才行？
展开解答
无摩擦的接触力都沿圆心连线，指向圆心。两个这样的力要相加为零，必须大小相等、方向相反，也就是两个接触点正好在一条直径的两端。一般的构型做不到，稍一扰动就失去平衡。三个方向不全相同的接触才能平衡一般的力，所以实验室把承力接触少于三个的圆盘算作松动盘。
用例题的数据，求两盘系统储存的弹性能，并用两种方法计算。
展开解答
串联弹簧整体： $U=f^2(R_1^2+\sigma^2+R_2^2)/(2\varepsilon)=2.29\times10^{-5}\,\varepsilon$ 。分开算：三个弹簧各为 $k\delta^2/2$ ，分别是 $2.63\times10^{-6}$ 、 $1.52\times10^{-5}$ 、 $5.16\times10^{-6}$ ，相加也是 $2.29\times10^{-5}\,\varepsilon$ 。

来源与复算范围

C. S. O’Hern, L. E. Silbert, A. J. Liu and S. R. Nagel, Jamming at zero temperature and zero applied stress: the epitome of disorder, arXiv:cond-mat/0304421v1; Physical Review E 68, 011306 (2003)
2003 年 4 月 17 日 v1。式 (1)：软球之间有限程的排斥势。第 II 节：从无穷高温的随机构型出发，在周期性盒子里用共轭梯度法求能量极小，并用能量判断是否未堵塞。图 5 至图 7 与式 (3)：堵住构型所占比例、堵塞阈值的分布，以及分布宽度随粒子数按 −0.55 次方减小。表 I：压强与接触数的指数和系数。图 20：二维堆积的力网络。我们用自己的程序复算了阈值宽度、压强指数和接触数指数；接触数指数在 128 个盘时没有复现。
C. P. Goodrich, A. J. Liu and S. R. Nagel, Finite-size scaling at the jamming transition, arXiv:1204.4915v3; Physical Review Letters 109, 095704 (2012)
2012 年 6 月 6 日 v3。式 (2)：除整体平移外不留任何零模所需的接触数。式 (3)：要有正的体积模量，还得再多一个接触，所以二维周期性盒子里 N 个盘至少要 2N − 1 个接触。式 (4)、(5)：相应的最少接触数。正文讨论低压时常见幂律失效的有限尺寸交叉。我们在自己生成的阈值堆积上核对了 2N − 1。
J. Guénolé, W. G. Nöhring, A. Vaid, F. Houllé, Z. Xie, A. Prakash and E. Bitzek, Assessment and optimization of the fast inertial relaxation engine (FIRE) for energy minimization in atomistic simulations and its implementation in LAMMPS, arXiv:1908.02038v4 (2020)
2020 年 1 月 30 日 v4。算法 2：FIRE 2.0，含半隐式 Euler 积分，以及走上坡后退回半步的修正。表 1：本课程使用的默认参数。这篇只提供弛豫方法，与堵塞本身无关。
Y. Zhao, J. Barés, H. Zheng, J. E. S. Socolar and R. P. Behringer, Shear-jammed, fragile, and steady states in homogeneously strained granular materials, arXiv:1904.10051v2; Physical Review Letters 123, 158001 (2019)
2019 年 8 月 14 日 v2。图 1(b)：21 个同心环组成的 Couette 剪切台，以每步 0.6% 的应变剪切 1447 到 2101 个双分散光弹圆盘（直径 1.59 cm 和 1.27 cm，摩擦系数 0.9）。图 2(d)、2(e)：用反向剪切区分脆弱态和剪切堵塞态，压强噪声阈值为 0.3 N/m。正文：能被剪切堵住的最小堆积分数为 0.745 ± 0.005，外推得到的无摩擦堵塞密度为 0.84 ± 0.02。这些是实验数据，我们只读图，没有复算。