基于外积求解颗粒平衡问题

文章编号  1000-5269（2024）01-0037-06

DOI：10.15958/j.cnki.gdxbzrb.2024.01.05

收稿日期：2022-07-07

基金項目：国家自然科学基金资助项目（11965007）；贵州省科技计划资助项目（ZK2022148）

作者简介：张兴刚（1980—），男，副教授，研究方向：统计物理和复杂系统的研究，E-mail：xgzhang@gzu.edu.cn.

*通讯作者：张兴刚，E-mail： xgzhang@gzu.edu.cn.

摘  要：提出一种简洁有效的基于外积的方法，对无摩擦的二维圆盘、二维凸颗粒以及三维球体的力学平衡方程进行求解和分析。理论结果表明，圆盘与球体上输出力的大小与输入力的大小成正比，因此引入了力传递的概念，并且导出力传递系数与接触结构的关系。对于凸颗粒，导出了输出力大小与输入力以及力矩的关系式。在已知圆盘上的平均正压力时，导出了各个接触力大小的公式并给出其几何意义。文中提出的方法以及计算的结果，对于求解其他物理问题中的向量方程有意义，对于进一步研究力学网络上的平衡方程有重要价值。

关键词：颗粒物质；力学平衡；张量代数；外积；力传递系数

中图分类号：O183.2；O312

文献标志码：A

静态颗粒体系是颗粒物理学中基础而重要的部分，其研究与许多实际的问题密切相关，如粮仓的设计、地基的承载能力、颗粒材料的力学性质……由于颗粒体系的复杂性，目前还没有一套完整的理论对其各方面的问题进行统一描述［1-2］。早期，研究者主要运用各种简化的理论模型解释有关的实验现象，例如粮仓效应的Janssen模型［3］和OSL模型［4］、力的几率分布的q模型［5］、应力凹陷的FPA模型 ［6］、……。为了从微观上通过颗粒间的接触结构和相互作用理解颗粒体系的宏观性质，一些研究者提出基于力网的方法和理论研究颗粒物质的静力学问题。Edwards等［7-8］提出了力系综理论，由接触网上颗粒的力学平衡方程推导应力张量的平衡方程以及应力几何本构方程。Snoeijer等［9］提出力网系综理论用于研究力的几率分布。自旋玻璃中的复本对称空腔（replica-symmetric cavity）理论也被用于研究无序颗粒体系中的力传递以及力的几率分布问题［8，10］。近年来，颗粒体系中矩阵形式的力学平衡方程［11］，无序颗粒体系中的力链 ［12］，晶格颗粒体系中应力的传递［13］等这些问题的研究加深了我们对于颗粒物质的力学性质和规律的理解。

在基于力网的理论研究中，一个重要的问题是求解和分析网络上的力学平衡方程组。方程组的复杂性来源于多个方面。一个重要的方面来源于方程组的具体形式与颗粒堆积的接触网络密切相关，而接触网络往往非常复杂，这就使得方程组的求解分析在一般情况下非常困难。另一个原因是力学平衡方程组往往是由多个向量方程而不是标量方程构成的方程组。前人在具体的计算中，经常把向量方程在具体坐标系中转化为标量方程进行求解和分析。但在很多时候，我们希望直接对向量方程进行处理，因为这样不仅能保证方程中物理量的整体性，还使得解的形式有容易理解的几何或物理意义。为了简单起见，本文主要考虑单个颗粒的力学平衡问题。提出基于外积的方法对多种典型情况下向量形式的力学平衡方程进行求解，并且对解的几何和物理意义进行讨论。

1  外积的基本性质

在物理学中，向量空间这一数学概念具有基础性的地位和非常重要的意义。比如我们熟悉的二维、三维空间可以用二维向量空间R2与三维向量空间R3描述，分析力学中往往涉及n维向量空间Rn，量子力学与无穷维向量空间密切相关。而张量代数是描述和研究向量空间的数学语言，因此张量对于物理学十分重要［14］。可以通过在多重向量空间上定义线性函数从而定义张量的概念，并引入加法、张量积（又称直积）、缩并、置换这些基本的运算。于是，标量（零阶张量）、向量（一阶张量）、二阶张量等都可以在张量空间的框架下进行理解；点乘、外积、叉乘（适用于三维）等运算也都可以基于上述基本运算进行定义和研究。用张量的语言进行描述有保持几何或物理对象的整体性、便于理解其几何或物理意义、使得关系式与坐标系的选择无关等优点。

为方便起见，以二维向量空间R2为例简要说明外积的定义和基本性质，为颗粒平衡问题的数学推导奠定一些基础。设向量a，b∈R2，则a与b的外积定义为

a∧b=ab－ba（1）

在本文中直积符号‘可以省略。由外积得到的结果是一个二阶完全反对称张量，或称为2-形式。类似地，r-形式就是r阶完全反对称张量；R2上的所有r-形式构成r-形式空间Λr（R2）。外积满足线性性、结合律、分配律等良好的性质，特别值得说明的是它满足反交换律

a∧b=－b∧a（2）

如图1所示，{e1，e2}是R2上的一组单位正交基底，规定逆时针绕向为正向。在R2上取两个单位向量

n（θ1）=cosθ1e1+sinθ1e2

n（θ2）=cosθ2e1+sinθ2e2 （3）

通过计算容易得

n（θ1）∧n（θ2）=sin（θ12）τ0（4）

其中，θ12表示由从n（θ1）转向n（θ2）的角度，τ0=e1∧e2叫做容积元，表示面积大小为1的有向正方形。在式（4）令θ12=π/2，可以看出容积元在转动操作下是不变的，因此可以把τ0理解为单位面积张量。显然a∧b=absin（θ）τ0，其中θ是从a到b夹角；因此a∧b表达了向量a与b所形成的平行四边形的面积张量。上述关于外积的基本性质，对于n维向量空间也一般是适用的。

2  颗粒力学平衡的求解

2.1  二维光滑圆盘

在光弹性颗粒、泡沫等体系的研究中，我们经常运用二维光滑圆盘组成的颗粒体系作为物理模型研究这些体系的力学性质。当体系在一定加载条件下达到力学平衡态后，可从中取出任意一个圆盘进行受力分析和求解。对于重力作用下的粮仓效应、力的几率分布等这类问题的求解，按力的传递这种方式讨论问题会带来很大的方便［5，15］。

如图2中的左图所示，在圆盘上选择两个接触点C1和C2；把作用在C1、C2上的接触力f1和f2的负值定义为输出力。这两个输出力可分别表示为－f1=f1n1，－f2=f2n2；其中n1和n2是单位向量，是已知的；f1和f2是输出力的大小，是未知的。除了f1和f2之外，圆盘一般还受到其他接触力或重力的作用，把这些力统称为输入力。设输入力的合力是f=fn，可等效地将它作用于圆盘中心，一般是已知的。为了由输入力获得输出力的大小f1和f2，由于光滑圆盘上的力矩平衡自动满足，因此只需考虑力平衡方程

f1n1+f2n2=fn（5）

即可。这是一个向量方程，在二维情况下等价于两个标量方程；当n1∧n2=sin（θ12）τ0≠0時它有唯一解。利用a∧a=0，以及R2上的各个2-形式之间只相差一个系数，将n2从左边外积上方程（5）两端，再进行整理得

f1=fn∧n2n1∧n2=fsinθf，2sinθ12（6）

类似地，将n1从左边外积上方程（5）两端，可导出

f2=fn1∧nn1∧n2=fsinθ1，fsinθ12（7）

力平衡方程（5）的含义是总的输入力等于总的输出力，因此可以用力传递的观点理解输入力与输出力的关系，不过这种传递是矢量传递而不是标量传递。由（6）、（7）两式，可以分别定义力传递系数

λf，2=sinθf，2sinθ12

λ1，f=sinθ1，fsinθ12 （8）

显然，力传递系数只由力向量间的夹角（或接触结构）确定，表征了某个输入力对某个输出力的传递能力。

对于堵塞态［8］颗粒体系的力学性质等问题，可按图2中的右图分析颗粒上的受力和平衡，这时一般不考虑重力。在有3个接触点的情况下，设接触力分别为f1=－f1n1，f2=－f2n2，f3=－f3n3，则由力学平衡条件有

f1n1+f2n2+f3n3=0（9）

我们一般还知道所加载的平均应力，这在典型情况下表现为已知平均正压力f，也就是说

f1+f2+f3=3f（10）

联立方程（9）与（10），在通常情况下应该可以确定唯一的f1、f2以及f3。利用（10）将f3代入（9）可得

f1（n1－n3）+f2（n2－n3）=－3fn3（11）

再利用（6）和（7）以及外积的性质可解出

f1=3fsinθ23sinθ12+sinθ23+sinθ31（12）

f2=3fsinθ31sinθ12+sinθ23+sinθ31（13）

将它们代入（10）可解得

f3=3fsinθ12sinθ12+sinθ23+sinθ31（14）

设圆盘的半径为r，三角形C1C2C3的面积张量为

S△C1C2C3τ0=12r2（n1∧n2+n2∧n3+n3∧n1）=

12r2（sinθ12+sinθ23+sinθ31）τ0 （15）

类似地，再考虑3个三角形的面积S△OC1C2、S△OC2C3、S△OC3C1，由上述结果容易得到

f1=3fS△OC2C3S△C1C2C3

f2=3fS△OC3C1S△C1C2C3

f3=3fS△OC1C2S△C1C2C3 （16）

这个公式以具有明显几何意义的方式给出了圆盘上3个接触力大小的解；显然，接触力的大小按照三角形面积的百分比进行分配。对于图2中右图所示的圆盘，如果接触点超过3个，各个接触点处力的大小将不能唯一确定，就会出现力的不确定性。

2.2  二维光滑凸颗粒

在图3中，给出了一个处于力学平衡状态的二维光滑凸颗粒，其中O点是该颗粒的质心。

类似于光滑圆盘的情形，为了考虑力的传递问题，在凸颗粒上选择3个接触点C1、C2和C3。以接触力f1、f2和f3的负值作为输出力，它们分别表示为－f1=f1n1、－f2=f2n2、－f3=f3n3，其中n1、n2和n3分别是凸颗粒在这3个接触点处的法向单位向量。设输入力的合力是f=fn，它等效地作用于质心O。与光滑圆盘不同的是，接触力还可能对质心产生力矩；为方便起见，这里采用反对称的力矩张量表达力矩，并且设输入力的总力矩张量为M=Mτ0。由颗粒的力平衡与力矩平衡得方程组

f1n1+f2n2+f3n3=f

f1r1∧n1+f2r2∧n2+f3r3∧n3=M  （17）

它等价于3个标量方程，因而由它一般可确定唯一的f1、f2以及f3。由（17）中的第2个方程可得

f3=M－f1r1∧n1+f2r2∧n2r3∧n3（18）

将其代入（17）中的第1个方程可得到类似于（5）的向量方程，利用（6）和（7）对其进行求解并且对结果进行化简，可得

f1=（f∧n2）（r3∧n3）+M（n2∧n3）+（n3∧f）（r2∧n2）（n1∧n2）（r3∧n3）+（n2∧n3）（r1∧n1）+（n3∧n1）（r2∧n2）

f2=（n1∧f）（r3∧n3）+（f∧n3）（r1∧n1）+M（n3∧n1）（n1∧n2）（r3∧n3）+（n2∧n3）（r1∧n1）+（n3∧n1）（r2∧n2）

f3=M（n1∧n2）+（n2∧f）（r1∧n1）+（f∧n1）（r2∧n2）（n1∧n2）（r3∧n3）+（n2∧n3）（r1∧n1）+（n3∧n1）（r2∧n2） （19）

可以看到，力的平衡与力矩平衡方程耦合在一起，使得解的结果要复杂很多。这时输出力的大小不只由输入力的大小及接触结构决定，还与输入力对质心产生的力矩有关，因此不能像光滑圆盘那样获得简单的力传递系数的定义。

2.3  三维光滑球体

外积的概念可用于任意n维向量，因此可以将上述二维情况下颗粒平衡的讨论推广到三维。一种简单且重要的情形是考虑若干光滑球体组成的静态颗粒体系。如图4所示，从中选取一个处于力学平衡状态的球体。

在球体上选择3个接触点C1、C2和C3，以它们的接触力的负值－f1=f1n1、－f2=f2n2、－f3=f3n3为输出力。设输入力的合力为f=fn，力矩平衡自动满足，所以只需考虑力平衡方程

f1n1+f2n2+f3n3=fn（20）

它相当于3个标量方程。为了求f1，将n2∧n3从右端乘以方程（20）的兩端，利用外积的运算性质以及R3上各个3-形式之间只相差一个系数，容易导出

f1=fn∧n2∧n3n1∧n2∧n3（21）

用类似的方法，可以导出

f2=fn1∧n∧n3n1∧n2∧n3

f3=fn1∧n2∧nn1∧n2∧n3 （22）

类似于光滑圆盘的情形，利用（21）与（22）的结果，我们可以定义与f1、f2、f3相对应的力传递系数λf，2，3、λ1，f，3、λ1，2，f。于是有

f1=λf，2，3f，f2=λ1，f，3f，f3=λ1，2，ff（23）

为了求出力传递系数与单位向量n、n1、n2、n3之间的夹角的关系，我们以n1∧n2∧n3为例计算这种三阶完全反对称张量。设Ω0是R3上的容积元，那么3-形式n1∧n2∧n3=VΩ0，其中V是一个实数，其绝对值代表由n1、n2、n3构成的平行六面体的体积。为讨论方便起见，假设n1∧n2∧n3与Ω0的定向是一致的，且n1、n2、n3线性无关，也就是说V>0。由（4）式易知

n1∧n2=sin（θ12）τ12（24）

其中，τ12代表由n1、n2确定的有向平面P（n1，n2）的容积元。设n3在平面P（n1，n2）上的投影向量可表示为l//n//，其中n//是单位投影向量；那么n3可表示为

n3=l//n//+l⊥n⊥（25）

其中，n⊥是垂直于P（n1，n2）的单位向量。设t//是P（n1，n2）上的与n//相垂直的单位向量，那么容积元Ω0可表示为Ω0=n//∧n//∧n⊥。另外，我们有τ12=n//∧t//，于是

n1∧n2∧n3=sin（θ12）l⊥Ω0（26）

这样，只需求出l⊥，就可以得出n1∧n2∧n3的结果。通过计算容易得到

τ12·n3=－l//t//（27）

因此，τ12作用于n3的效果是把n3投影到平面P（n1，n2）上，然后逆向旋转90度。利用这一几何意义，可以构造二阶对称张量

E12=－τ12·τ12（28）

显然，E12·n3=l//n//；因此，E12就是与平面P（n1，n2）相对应的投影张量，而且我们有E12·E12=E12。有了E12之后，可以得到

l2//=（n3-E12·n3）·（n3-E12·n3）=

1+n3·τ12·τ12·n3（29）

由（24）式可得

τ12·τ12=cosθ12（n1n2+n2n1）－（n1n1+n2n2）sin2θ12（30）

综合（26）、（29）和（30），可以算出

n1∧n2∧n3=

1+2cosθ12cosθ23cosθ31－cos2θ12－cos2θ23－cos2θ31Ω0 （31）

将这个公式中角度的指标进行适当替换，便可写出n∧n2∧n3的结果。于是力传递系数λf，2，3与n、n1、n2、n3之间的夹角的关系式为

λf，2，3=

1+2cosθf，2cosθ23cosθ3，f－cos2θf，2－cos2θ23－cos2θ3，f1+2cosθ12cosθ23cosθ31－cos2θ12－cos2θ23－cos2θ31（32）

类似地，可以写出λ1，f，3、λ1，2，f的公式。由此可见，三维情况下的力传递系数比二维情况要复杂很多。

3  总结

本文提出基于张量代数中的外积的方法对几种典型情况下颗粒的力学平衡方程进行求解和分析。可以看到，利用外积的方法可以以简洁明了的方式对一些向量方程进行求解，而且解的形式只由已知量之间的内在关系确定，与坐标系无关；这使得解的结果往往有良好的几何或物理含义，便于理解和分析。具体的求解表明，只有对于光滑圆盘和光滑球体这样简单的情形，才能严格地引入力传递的概念以及力的传递系数。对于其他形状的颗粒以及有摩擦的情形，这时力平衡方程与力矩平衡方程往往耦合在一起，使得输出力不只与输入力有关，还与输入力的力矩有关。不过在讨论颗粒的力学平衡时，如果可以忽略力矩平衡条件带来的约束，也可以近似地使用力传递系数简化问题的研究。论文中提出的方法以及计算的结果，对于求解其他物理问题中的向量方程有意义，对于进一步研究网络上的力学平衡方程组有价值。

感谢贵州大学实验室开放项目资助。

参考文献：

[1]陆坤权， 刘寄星. 颗粒物质（上）［J］. 物理， 2004， 33（9）： 629-635.

［2］ 孫其诚， 厚美瑛， 金峰， 等. 颗粒物质物理与力学［M］. 北京： 科学技术出版社， 2011.

［3］ DE GENNES P G. Granular matter： a tentative view［J］. Reviews of Modern Physics， 1999， 71（2）： 374-382.

［4］ VANEL L， CLAUDIN P， BOUCHAUD J P， et al. Stresses in silos： comparison between theoretical models and new experiments［J］. Physical Review B： Condensed Matter & Materials Physics， 2000， 84（7）： 1439-1442.

［5］ COPPERSMITH S N， LIU C， MAJUMDAR S， et al. Model for force fluctuations in bead packs［J］. Physical Review， 1996， 53（5）： 4673-4685.

［6］ WITTMER J P， CLAUDIN P， CATES M E， et al. An explanation for the contral stress minimum in sand piles［J］. Nature， 1996， 382（25）： 336-338.

［7］ EDWARDS S F， OAKESHOTT R B S. Theory of powders［J］. Physica A， 1989， 157（3）： 1080-1090.

［8］ BAULE A， MORONE F， HERRMANN H J， et al. Edwards statistical mechanics for jammed granular matter［J］. Reviews of Modern Physics， 2018， 90（1）： 15001-15006.

［9］ SNOEIJER J H， VLUGT T J H， VAN HECKE M， et al. Force network ensemble： a new approach to static granular matter［J］. Physical Review Letters， 2004， 92（5）： 054302.

［10］BO L， MARI R， SONG C M， et al. Cavity method for force transmission in jammed disordered packings of hard particles［J］. Soft Matter， 2014， 10（37）： 7379-7392.

［11］GENDELMAN O， POLLACK Y G， PROCACCIA I， et al. What determines force chains in granular media？［J］. Physical Review Letters， 2016， 116（2）： 078001.

［12］KRISHNARAJ K P， NOTT P R. Coherent force chains in disordered granular materials emerge from a percolation of quasilinear clusters［J］. Physical Review Letters， 2020， 124（19）： 198002.

［13］ZHANG X G， DAI D. Exact solutions for stress distribution in rhombic disk packings［J］. Granular Matter， 2021， 23（2）： 51.

［14］郭仲衡. 张量（理论和应用）［M］. 北京： 科学技术出版社，1988.

［15］张兴刚， 胡林， 隆正文. 粮仓效应的圆盘堆积模型［J］. 计算物理， 2006， 23（6）： 642-646.

（责任编辑：曾  晶）

Solving Particle Equilibrium Problem Based on Wedge Product

ZHANG Xinggang*

（College of Physics， Guizhou University， Guiyang 550025， China）

Abstract：

A simple and effective method based on exterior product is proposed to solve and analyze the mechanical equilibrium equations for frictionless particles with different shapes such as two-dimensional disk， two-dimensional convex particle and three-dimensional sphere. The theoretical results show that the output force on a disk or a sphere is proportional to the input force. Therefore， the concept of force transmission is introduced， and the relationship between the coefficient of force transmission and the contact structure is derived. For case of convex particle， the relationship between the output force and the input force and torque is derived. When the average normal pressure on a disk is known， the formulas of the contact forces are derived and their geometric significances are given. The method and results presented in this paper are also significant for solving vector equations in other physical problems and for further studying the equilibrium equations on mechanical networks.

Key words：

granular matter; mechanical equilibrium; tensor algebra; wedge product; coefficient of force transmission