颗粒堵塞：一堆散粒怎样变成固体

半径为

R_i

和

R_j

的两个软圆盘，圆心相距

r

。重叠量

\delta=R_i+R_j-r

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

从两个圆盘，到堵住的一整盒

第一章求两个重叠圆盘之间的力，再看挤满圆盘的盒子怎样让载荷沿少数几条强路径走。第二章数一数，堆积至少要多少个接触才能开始抵抗。第三章把接触力换算成压强，并和壁上的压强核对。第四章从随机起点找出堵塞点，推出过了堵塞点压强怎样增长。第五章加上摩擦，问为什么只靠剪切也能把堆积堵住。第六章把堆积压紧再放开，看它记住了什么。第七章剪切堵塞堆积，检验最软的振动能否预告它在哪里滑移。最后一章逐图读 O’Hern 等人 2003 年的论文。

需要牛顿定律、胡克定律、平面矢量，会读双对数坐标图。应力张量、标度律和计数论证在第一次用到时再介绍。 FR27

01 · L1
两个圆盘怎样互相推，一百个圆盘怎样分担载荷
从重叠量求接触力，解出夹在两壁之间的两个圆盘，再看 128 个圆盘挤紧后，力为什么沿着少数路径走。
02 · L2
堵住的堆积需要多少个接触
把受力平衡写成方程组，数未知量和方程，推出周期性盒子里的最少接触数 2N′ − 1，再用它解释松动盘和平均接触数接近 4 的原因。
03 · L2
怎样从接触力算出压强
推导二维接触应力公式，证明它等于壁上单位长度的力，也等于能量对面积的导数，并用两个圆盘的例子把三种算法对上。
04 · L3
堵塞点在哪里，越过它之后压强怎样增长
用随机初态和二分搜索求出每个堆积的堵塞阈值，看阈值分布怎样随盘数变窄，再推导并检验压强和接触数在阈值之上的标度律。
05 · L3
加上摩擦后，剪切怎样让稀一些的颗粒堵住
写出带摩擦的接触模型，用约束计数说明摩擦为何降低堵塞所需的接触数，再读 Zhao 等人 2019 年怎样用反向剪切区分脆弱态和剪切堵塞态。
06 · L3
压进去再放出来，堆积记住了什么
比较压缩和卸载两条路径，算清一次循环里壁做的功去了哪里，再看反复压缩的堆积怎样落进一个每轮都相同的循环，以及反向剪切怎样检验剪切堵塞。
07 · L4
剪切下的堆积，下一次会在哪里滑移
缓慢剪切堵塞堆积，它会一跳一跳地滑移。推导最软的振动模为什么先变软，再在没见过的堆积上检验它能否提前点出要动的盘。
08 · L4
逐图读 O’Hern 等人 2003 年的论文，找出堵塞点
逐图读这篇论文：构型怎样生成、每个构型的阈值怎样测、阈值分布怎样外推到无限大系统，哪些常数是拟合出来的，哪些我们复现了。

沿途的实验

01挤压一盒圆盘把 128 个软圆盘压过堵塞点再放回，看接触、压强和力网络怎样出现又消失。02让脉冲穿过一排颗粒给一排珠子中的第一颗一份动能，看脉冲怎样穿过只推不拉的 Hertz 接触。03预测下一次滑移缓慢剪切 64 个软圆盘，看最软的模在每次滑移前失去刚度，再统计它事先点中移动盘的比例。04只靠剪切堵住松散颗粒在堵塞密度以下剪切 128 个带摩擦的圆盘，看力链沿被挤压的对角线长出来，再关掉摩擦或反向剪切。

压缩实验模拟两壁之间的无摩擦软圆盘：每一步把盒子压一点，再让圆盘弛豫到所有力平衡。压强和接触数来自模拟；用来比较的发表拟合值单独列在来源里。剪切实验把同样的圆盘放进周期单元；命中率来自模拟，在另外一组堆积上记录的平均值画成参考刻线。剪切堵塞实验给圆盘加上摩擦和滚动；它的压强完全来自模拟，旁边不画任何发表曲线。

来源与复算范围

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。
L. E. Silbert, D. Ertaş, G. S. Grest, T. C. Halsey, D. Levine and S. J. Plimpton, Granular flow down an inclined plane: Bagnold scaling and rheology, arXiv:cond-mat/0105071v1; Physical Review E 64, 051302 (2001)
2001 年 5 月 3 日 v1。第 II 节：带摩擦的接触模型，法向是带阻尼的弹簧，切向是受 Coulomb 条件截断的弹簧，切向与法向刚度之比为 2/7。式 (12)：应力张量等于分支矢量乘接触力之和除以体积，再加动能项。静态堆积只用式 (12) 的接触部分。
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：本课程使用的默认参数。这篇只提供弛豫方法，与堵塞本身无关。
K. Shundyak, M. van Hecke and W. van Saarloos, Force mobilization and generalized isostaticity in jammed packings of frictional grains, arXiv:cond-mat/0610205v1; Physical Review E 75, 010301(R) (2007)
2006 年 10 月 7 日 v1。开头几段：有摩擦时，等静定接触数从 2d 降到 d + 1。式 (1) 与图 1：摩擦力恰好达到 Coulomb 上限的“完全动员”接触，每个多出一个约束，所以每个颗粒上这类接触数至多为 z 减 d + 1。图 2：动员程度的分布分成上限处的一个峰和其下的一段平台。我们没有复算他们的模拟。
D. Bi, J. Zhang, B. Chakraborty and R. P. Behringer, Jamming by shear, Nature 480, 355–358 (2011)
只读了摘要。低于各向同性堵塞密度的带摩擦光弹圆盘，仅靠剪切形变就能经过脆弱态进入剪切堵塞态，力网络变得各向异性。课程只按摘要和 Zhao 等人（2019）的转述引用这一结论，没有读全文，也没有复算其数据。
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。这些是实验数据，我们只读图，没有复算。
C. E. Maloney and A. Lemaître, Amorphous systems in athermal, quasistatic shear, arXiv:cond-mat/0510677v2; Physical Review E 74, 016118 (2006)
2005 年 10 月 28 日 v2。第 IV 节：能量极小值与鞍点合并时，弹性段就到了尽头。式 (9)、(10)：沿临界模的位移按剩余应变的平方根增长。图 5(b)：动力学矩阵的最低本征值随事件前剩余应变按平方根减小。我们在自己的剪切堆积上检验了这个平方根趋近。
M. L. Manning and A. J. Liu, Vibrational modes identify soft spots in a sheared disordered packing, arXiv:1012.4822v1; Physical Review Letters 107, 108302 (2011)
2010 年 12 月 21 日 v1。2500 个 Hertz 圆盘，堆积分数 0.95，准静态剪切。图 2：取最低 30 个振动模中各自位移最大的 20 个颗粒，拼成“软点”。式 (1)：软点与下一次重排的重合分数；图 2 的例子得分 0.64。图 4(a)：最佳软点的得分随事件前剩余应变的变化。正文：软点处接触数少 6%，自由体积多 103%，但单靠这些局部量找不到它们。课程的盲测用更小的堆积和自己选定的阈值。