堵塞点在哪里，越过它之后压强怎样增长

怎样求一个堆积的堵塞阈值

O’Hern 等人的协议是这样的：在周期性盒子里随机撒下圆心，相当于无穷高温下的构型；然后求能量的局部极小值。若每个盘平均能量低于 $10^{-16}\,\varepsilon$ ，圆盘之间没有重叠，堆积没有堵住；否则就堵住了。

他们用共轭梯度法求极小值，我们用 FIRE 2.0，两者找的是同一种东西：力处处平衡的构型。没堵住就把所有圆盘按同一比例放大，堵住了就缩小，每次切换方向时把步长减半，直到步长小于十万分之一。最后停下的堆积分数就是这个初态的阈值。

阈值是一个分布，盘数越多越窄

每个随机初态给出一个阈值。O’Hern 等人图 5 画出已堵住初态的比例怎样随堆积分数变化，对它求导就得到图 6 的阈值分布。分布的半高全宽在大系统里按式 (3) 减小，图 7 给出 $w_0=0.16\pm0.04$ 、 $\Omega=0.55\pm0.03$ 。

我们对 32、64、128、256 个盘各取 24 个随机种子，用样本标准差乘 $2\sqrt{2\ln2}$ 估计半高全宽，这一步假设分布是高斯形。结果依次为 0.0259、0.0148、0.0129、0.00654，都落在式 (3) 按误差给出的范围内。平均阈值依次为 0.827、0.835、0.838、0.839，随盘数增加，向图 5(a) 中大系统约 0.84 的位置靠近。

盘数趋于无穷时，宽度趋于零：几乎所有随机初态都在同一个堆积分数堵住。这就是“堵塞点”有明确位置的依据。

w=w_0\,N^{-\Omega},\qquad w_0=0.16\pm0.04,\quad \Omega=0.55\pm0.03

O’Hern 等人式 (3)，只适用于大系统。他们指出，盘数很少时宽度先随盘数增大，过了最大值才按这个幂律减小。

来源：C. S. O’Hern

越过阈值，压强怎样增长

记 $\Delta\phi=\phi-\phi_c$ 。挤过阈值一点，典型的重叠量与 $\Delta\phi$ 成正比，因为多出来的面积要由重叠吸收。第一章说过 $f\propto\delta^{\alpha-1}$ ；接触数在阈值附近几乎不变，第三章的压强是接触力乘距离之和，所以 $p\propto\Delta\phi^{\alpha-1}$ 。谐振接触时压强线性增长，Hertz 接触时按 3/2 次方增长。

O’Hern 等人表 I 给出二维双分散系统：谐振 $\psi=1.01$ 、 $p_0=0.34$ ，Hertz $\psi=1.50$ 、 $p_0=0.27$ 。我们用 128 个盘、8 个种子，在 $\Delta\phi$ 从 0.001 到 0.1 的范围拟合，得到谐振 $1.039\pm0.021$ 、 $p_0=0.35$ ，Hertz $1.52\pm0.02$ 、 $p_0=0.233$ 。指数复现了；Hertz 的系数低于表 I 的误差范围。

p=p_0\,(\phi-\phi_c)^{\psi},\qquad \psi=\alpha-1

表 I 中

p_0

的误差为 ±0.05，

\psi

的误差为 ±0.03。上面的论证只说明指数，系数要靠计算。

压强对“超过各自阈值的压缩量”作图，两轴均为对数坐标。圆点：本课程复算的 128 个盘、8 个种子的平均值。虚线：O’Hern 等人表 I 的拟合，谐振接触

0.34\,\Delta\phi^{1.01}

（青色），Hertz 接触

0.27\,\Delta\phi^{1.50}

（琥珀色），作为单独的参照。

来源：C. S. O’Hern

接触数的标度没有复现

表 I 还给出 $Z-Z_c=Z_0\,\Delta\phi^{\zeta}$ ，二维谐振 $\zeta=0.49$ 、 $Z_0=3.6$ 。我们在 128 个盘时得到 $\zeta=0.577\pm0.009$ 、 $Z_0=4.41$ ，明显偏陡。

Goodrich、Liu 和 Nagel 解释了小系统的问题：当 $p^{1/2}N$ 只有 1 的量级时，接触数的平方根律失效。我们只拟合 $p^{1/2}N\ge3$ 的数据，指数仍然偏陡，所以只能说它在 128 个盘时没有复现，要验证需要更大的系统。

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

完整例题：挤过阈值百分之一

问题：二维谐振双分散圆盘，挤到 $\Delta\phi=0.01$ ，压强和平均接触数约是多少？

用表 I： $p=0.34\times0.01^{1.01}=3.2\times10^{-3}\,\varepsilon/d^2$ ； $Z=3.98+3.6\times0.01^{0.49}=3.98+0.38=4.36$ 。

与我们的数值比较。在 64、128、256 个盘的周期性盒子里各取 12 个种子，挤到 $\phi_c+0.01$ ，平均压强分别为 $(2.62\pm0.22)\times10^{-3}$ 、 $(2.73\pm0.13)\times10^{-3}$ 、 $(2.45\pm0.23)\times10^{-3}$ ，比表 I 的拟合值低约 15% 到 25%，且不随盘数系统地变化。表 I 的二维数值来自 1024 个盘在一段压缩范围上的整体拟合，单个点差这么多不算矛盾。表 I 的系数应当看作准确到两成左右的估计。

来源：C. S. O’Hern

排除假的堵塞：结晶

同样大小的圆盘容易排成六角晶格，堵住的位置会高得多。我们给每个盘算局部六重序，超过 0.85 的算作“像晶体”。挤过阈值 0.01 时，双分散堆积里这类盘平均占 14.5%，最多 21.1%；等大圆盘平均占 64.6%，平均阈值升到 0.859。

双分散堆积只要像晶体的盘超过 35%，就被标记为结晶，不计入统计。所有基准种子都没有超过这个线。实验室的“等大圆盘”预设专门展示这个反例。

来源：C. S. O’Hern

到实验里试一试

比较谐振和 Hertz 的起点

先预测：刚过起点时，哪种接触定律的压强曲线更平？
播放默认预设，看压强图离开零点后是直线还是向上弯。再选“Hertz 接触”预设，比较两条曲线在起点附近的形状。
选“等大圆盘”预设，读出起点堆积分数和平均局部六重序，与默认预设比较，说明为什么基准统计要排除结晶。

练习与解答

把挤压量从阈值之上 0.01 加到 0.02，谐振和 Hertz 接触的压强各变成几倍？
展开解答
按表 I 的指数，谐振约 $2^{1.01}=2.01$ 倍，Hertz 约 $2^{1.50}=2.83$ 倍。系数 $p_0$ 在比值里约掉了。
按式 (3)，要让阈值分布的半高全宽小于 0.001，大约需要多少个圆盘？
展开解答
$N=(w_0/w)^{1/\Omega}=160^{1/0.55}\approx1.0\times10^{4}$ 。考虑到 $w_0$ 和 $\Omega$ 的误差，答案可能相差几倍，只能说是一万个左右。
32 个盘的周期性盒子里，单个随机初态的阈值大致会落在什么范围？只做一次能说明大系统的 $\phi_c$ 吗？
展开解答
我们测得平均 0.827、标准差 0.011，单次结果大致落在 0.81 到 0.84 之间。它比大系统约 0.84 的值系统地偏低，单次结果的涨落又和这个偏差差不多大。要得到大系统的 $\phi_c$ ，必须用更大的系统。

来源与复算范围

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：本课程使用的默认参数。这篇只提供弛豫方法，与堵塞本身无关。