--- title: "逐图读 O’Hern 等人 2003 年的论文，找出堵塞点" description: "逐图读这篇论文：构型怎样生成、每个构型的阈值怎样测、阈值分布怎样外推到无限大系统，哪些常数是拟合出来的，哪些我们复现了。" canonical_url: "https://sandboxphysics.com/zh/stories/granular-jamming/reading-ohern-2003/" language: "zh-CN" page_type: "story" alternate_languages: en: "https://sandboxphysics.com/stories/granular-jamming/reading-ohern-2003/" last_modified: "2026-10-29" --- [颗粒与堵塞](https://sandboxphysics.com/zh/stories/granular-jamming/) L4 · 第 8 章 # 逐图读 O’Hern 等人 2003 年的论文，找出堵塞点 Sandbox Physics · 2026-10-29 这篇论文用一种几乎没有可调参数的模型，把“随机密堆积”变成了一个可以精确定位的点，叫作 J 点。读它时要分清三类数：直接测出来的、拟合得到的、在某个假设下才成立的。本章按论文第 II 节的顺序，一张图一张图地读。 ## 开始前需要什么 前面各章：接触计数（第二章）、接触应力（第三章）、阈值分布与压强标度（第四章）。会读双对数坐标，知道拟合参数带不确定度。 - 说出论文生成零温构型的协议，以及它用什么判断“堵住了”。 - 读懂图 5 到图 8 之间的推导链：堵住比例 → 阈值分布 → 宽度 → 无限大系统的极限。 - 分清表 I 中哪些数依赖假设，并用我们的复算说明哪些复现了、哪些没有。 ## 先看协议：构型从哪里来 第 II A 节和式 (1) 定义了模型：两个颗粒只有重叠时才相互排斥，势能为 $\varepsilon(1-r_{ij}/\sigma_{ij})^{\alpha}/\alpha$。二维用大小比 1.4、各占一半的两种圆盘，防止结晶；粒子数从 4 到 4096。 构型从无穷高温出发：圆心在周期性盒子里完全随机地放，再用共轭梯度一路往能量低处走。停止条件有两个：每个粒子的能量低于 $10^{-16}$，或者相邻两次迭代的能量差小于 $10^{-15}$。第一个条件就是“未堵塞”的判据，对应压强低于 $10^{-10}$。 读到这里要记住：论文里的“堵塞”是这个协议下的结论。换一种生成方法，比如先慢慢降温，得到的阈值会不同。论文后面所有的分布和外推，都是对这种随机起点的统计。 来源：C. S. O’Hern ## 一个构型的阈值怎样测（第 II B 节） 对每个零温构型：压强为零就把所有颗粒按同一比例放大，压强不为零就缩小，每步之后重新求极小，直到压强刚好变为非零。步长在 $10^{-6}$ 到 $10^{-4}$ 之间，越接近阈值越小，并且逐次减小步长，确认最后落在同一个构型上。 这一节的关键结论是：体积模量和剪切模量在同一个堆积分数上变为零。也就是说，“开始抗压”和“开始抗剪”是同一件事，阈值对每个构型都有确切定义。作者还检查了动力学矩阵：阈值之上至少 1e-6 处，除了松动颗粒和整体平移，没有零频模式。 图 3 有一处要小心。我们读的 arXiv v1 中，图注说参考线斜率是 2.0 和 2.5；而表 I 的压强指数是 1.01 和 1.50，第四章从 $p\propto\Delta\phi^{\alpha-1}$ 推出的也是这两个数。2.0 和 2.5 恰好是能量的指数 $\alpha$。我们按表 I 检验压强；正式发表版的图注我们没有看过。 来源：C. S. O’Hern ## 从堵住比例到阈值分布（图 5 到图 8） 图 5 画的是在给定堆积分数下，随机起点中已经堵住的比例 $f_j(\phi)$。它是一条从 0 升到 1 的台阶曲线。图 6 是它对堆积分数的导数，也就是阈值的分布 $P_j(\phi_c)$。作者把三种接触指数的分布叠在一起，在误差内完全重合：阈值只由几何决定，与接触力的软硬无关。 图 7 是分布的半高全宽随粒子数的变化。读图时注意两段：粒子很少时宽度先随粒子数变大，越过最大值后才按式 (3) 下降。式 (3) 的 $w_0=0.16\pm0.04$、$\Omega=0.55\pm0.03$ 只适用于后一段。 图 8 处理的是分布峰的位置 $\phi_0$。它随系统长度 $L=N^{1/d}$ 按式 (4) 趋近极限值 $\phi^{*}$，拟合给出 $\nu=0.71\pm0.08$、$\delta_0=0.12\pm0.03$。可见 $\phi^{*}$ 是外推出来的，不是哪一次模拟直接测到的。 $$ w=w_0\,N^{-\Omega},\qquad \phi_0-\phi^{*}=\delta_0\,L^{-1/\nu},\qquad L=N^{1/d} $$ 论文的式 (3) 与式 (4)。宽度趋于零、峰位置趋于 $\phi^{*}$，合起来说明在无限大系统中几乎所有随机起点都在同一个堆积分数堵住。图 8 画的是两者之差的绝对值。 来源：C. S. O’Hern ## 接触数与表 I：哪些数依赖假设 第 II E 节用式 (6) $Z-Z_c=Z_0(\phi-\phi_c)^{\zeta}$ 拟合图 9，正文给出 $\zeta=0.50\pm0.03$。表 I 里的 $Z_c=3.98$ 是另一回事：作者先去掉松动颗粒（约 5%），再假定 $\zeta=0.5$，才得到这个“精确”的值。$Z_c$ 与 4 相符，并不能反过来独立证明指数是 1/2。 Goodrich、Liu 与 Nagel 2012 年把计数做得更细：周期性盒子里 N 个圆盘，要有正的体积模量至少需要 $2N-1$ 个接触，并且当 $p^{1/2}N$ 接近 1 时，接触数的平方根律会失效。读 2003 年的表 I 时，要把这条有限尺寸限制带上。 我们的复算：在阈值处的承载网络上，超出最少接触数的接触在 32 和 64 个盘时都是 0，在 128 和 256 个盘时最多 1 个；松动盘平均占 4.0% 到 5.0%，与论文的约 5% 一致。接触数指数在 128 个盘时得到 0.577 ± 0.009，没有复现 0.49。 来源：C. S. O’Hern · C. P. Goodrich ## 完整例题：用式 (3) 和式 (4) 检查我们的阈值统计 问题：我们在周期性盒子里对 32、64、128、256 个圆盘各取 24 个随机起点，测得阈值分布的半高全宽为 0.0259、0.0148、0.0129、0.00654，平均阈值为 0.827、0.835、0.838、0.839。它们与式 (3)、式 (4) 相容吗？无限大系统的阈值大约是多少？ 第一步，宽度。式 (3) 的中心值给出 $0.16\times N^{-0.55}$：0.0238、0.0162、0.0111、0.00758。取 $w_0$ 和 $\Omega$ 不确定度的两端，256 个盘允许的范围是 0.0048 到 0.0112，我们的 0.00654 在里面；其余三个尺寸也都在各自的范围里。注意我们是用标准差乘 $2\sqrt{2\ln2}$ 估计全宽，假设了分布是高斯形。 第二步，外推。二维时 $L=\sqrt{N}$，256 个盘 $L=16$，式 (4) 的中心值给出 $0.12\times16^{-1/0.71}=0.0024$。平均阈值随尺寸上升，说明峰从下方趋近极限，所以 $\phi^{*}\approx0.839+0.0024=0.841$。按 $\delta_0$ 和 $\nu$ 的不确定度，这个修正在 0.0011 到 0.0045 之间。 结论：宽度与式 (3) 相容；外推值约 0.841，落在图 5(a) 中大系统 0.84 附近。这个估计有两处近似：式 (4) 讲分布的峰，我们用 24 个样本的平均值代替；256 个盘本身也只是中等大小。 来源：C. S. O’Hern ## 读这类论文的清单 一、先找协议：构型怎样生成、用什么判断堵住、收敛阈值多大。二、分清直接测量和拟合：J 点的位置、接触数的极限都来自外推或假设。三、找适用范围：式 (3) 只对大系统成立；按图 3 和图 9 的图注，二维的标度数据来自 1024 个圆盘。四、对照图注和正文：像图 3 这样的出入，以表和推导为准，并记下自己没有核实的部分。 五、复算时写明什么算相符。我们要求宽度落在式 (3) 不确定度给出的范围里，接触数指数落在表 I 的 $\pm0.04$ 之内。前者通过；后者 0.577 对 0.49，没通过。两个结果都记录在案。 来源：C. S. O’Hern · C. P. Goodrich ## 到实验里试一试 ### [在刚刚堵住的堆积里找表 I 的量](https://sandboxphysics.com/zh/labs/disk-jamming-compression/) 1. 先预测：刚过阈值时，超出等静定的接触多还是少？松动盘大约占几成？ 2. 选“刚刚堵住”预设，播放到最大堆积分数，读出“超出等静定的接触数”；再在“更多数值与计算方法”里读出松动盘和承载盘的数目。 3. 换“深度挤压”预设重复一遍。哪个量随压缩明显变化？把它和式 (6) 的 Z − Zc 对应起来，并说明实验室用的是有壁的盒子，起点和论文的周期性盒子不同。 ## 练习与解答 1. 按式 (3) 的中心值，多少个圆盘时阈值分布的全宽会小于 0.001？ 展开解答 $N=(0.16/0.001)^{1/0.55}=160^{1.82}\approx1.0\times10^{4}$。这是外推，论文的模拟最多到 4096 个粒子。 2. 表 I 给出 Zc = 3.98，与等静定的 4 相差不到不确定度。为什么这不能独立证明 ζ = 1/2？ 展开解答 作者是在假定 ζ = 0.5 并去掉松动颗粒之后才得到这个 Zc 的。指数换一个值，外推出的 Zc 也会变；要独立检验，需要同时拟合 Zc 和 ζ，或者像 Goodrich 等人那样直接数接触。 3. 我们的 128 盘数据给出接触数指数 0.577，论文是 0.49。列出两个可能的原因，并说哪一个我们已经排除。 展开解答 一是有限尺寸：当 $p^{1/2}N$ 接近 1 时平方根律失效；我们只用了 $p^{1/2}N\ge3$ 的数据，指数仍然偏大，所以这个原因没有被完全排除，但不能单独解释。二是计数方式：论文去掉松动颗粒后假定 ζ = 0.5 求 Zc，我们用 Goodrich 等人的最少接触数作参考，两种参考点不同会改变拟合。要分开这两个原因，需要更大的系统。 [上一章：剪切下的堆积，下一次会在哪里滑移](https://sandboxphysics.com/zh/stories/granular-jamming/predicting-rearrangements/) ## 来源与复算范围 1. [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)](https://arxiv.org/abs/cond-mat/0304421v1) 2003 年 4 月 17 日 v1。式 (1)：软球之间有限程的排斥势。第 II 节：从无穷高温的随机构型出发，在周期性盒子里用共轭梯度法求能量极小，并用能量判断是否未堵塞。图 5 至图 7 与式 (3)：堵住构型所占比例、堵塞阈值的分布，以及分布宽度随粒子数按 −0.55 次方减小。表 I：压强与接触数的指数和系数。图 20：二维堆积的力网络。我们用自己的程序复算了阈值宽度、压强指数和接触数指数；接触数指数在 128 个盘时没有复现。 2. [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)](https://arxiv.org/abs/1204.4915v3) 2012 年 6 月 6 日 v3。式 (2)：除整体平移外不留任何零模所需的接触数。式 (3)：要有正的体积模量，还得再多一个接触，所以二维周期性盒子里 N 个盘至少要 2N − 1 个接触。式 (4)、(5)：相应的最少接触数。正文讨论低压时常见幂律失效的有限尺寸交叉。我们在自己生成的阈值堆积上核对了 2N − 1。