湿颗粒系统是密相固体颗粒和填隙流体组成的多相结构，与干颗粒系统相比较，接触颗粒间的液体形成液桥，产生黏性力，从而影响颗粒运动。在很多工程实际中，颗粒间含水不可避免且会对颗粒材料宏观特性产生显著影响，如漏斗流，颗粒结拱，料口阻塞，仓壁磨损等^[1]。

湿颗粒材料按照含液量可将间隙水分为摆动状态、链索状态、毛细状态以及浸渍状态^[2]。当湿颗粒材料的固相采用离散颗粒模型时，间隙水的影响可基于液桥模型^[3-4]或连续介质模型^[5-6]考虑。目前已有不少学者通过液桥模型研究湿颗粒材料的力学特性，如杜友耀等^[7]研究二维条件下两相邻颗粒间液桥临界断裂距离及液桥力的计算过程，王辉等^[8]利用液桥理论研究了湿颗粒分离过程中在液桥力作用下颗粒速度的变化规律，Zhang等^[9]建立了考虑滑动摩擦和恢复系数的颗粒接触模型，Louati等^[10]通过定性和定量的方法研究了颗粒尺寸及含液量对液桥力的影响。同时也有一些学者对筒仓卸料过程中储料的流态和侧压力开展研究，如刘克瑾等^[11]基于离散元方法模拟筒仓储料卸料成拱过程并研究筒仓壁压力分布情况，陈磐等^[12]基于离散元方法研究料仓卸料流型与湿颗粒摩擦特性、料仓半顶角、料仓开口尺寸等之间的相互关系，Madrid等^[13]通过离散元方法对比研究了筒仓强制卸料与自由卸料时的速度分布差异。

近些年的研究学者对筒仓卸料压力的离散元模拟研究比较充分，但是针对湿颗粒卸料的力学研究目前还不够充分。Soulie等^[14]通过gorge method计算了大小不等球形颗粒之间的液桥力，但是计算过程比较复杂。笔者尝试将半径不等的2个颗粒通过几何平均的方式简化为2个等效半径相等颗粒计算液桥力，提供一种既能考虑半径多样性，也能保证求解效率的近似途径。同时，在楚锡华等^[15]模拟干颗粒材料筒仓卸料过程的基础上，介绍湿颗粒液相的液桥模型和基于离散颗粒模型计算液桥力的主要途径，用自定义接触本构的方式对颗粒间含液量较少情况下的液桥作用进行数值实现，并且在此基础上模拟筒仓卸料过程，比较干颗粒和湿颗粒在筒仓卸料过程中宏观特性的不同。

1 湿颗粒液桥理论模型

液桥力包括静态的毛细力和动态的黏性力，根据情况筒仓卸料过程的模拟只考虑毛细力。基于Fisher提出的喇叭形及Laplace-Young方程^[16]，用一种简单而便于实现的液桥计算模型。颗粒间隙液桥的几何表征如图 1所示，R₁，R₂为颗粒半径，D为颗粒间距，δ₁，δ₂为半填充角或圆心角，θ为气液固三相交接处的颗粒表面切向与液体表面切向的夹角。该液桥模型计算液桥力可基于如下假定：1)颗粒为圆形；2)颗粒表面光滑不考虑粗糙度；3)液桥在摆动状态下含水量很少，液桥体积较小不考虑重力；4)颗粒材料的准静态加载过程中液桥力可忽略液体黏性的影响。基于以上假设，当颗粒半径相等且颗粒接触(D=0)，θ=0时，有R₁=R₂=R，δ₁=δ₂=δ，可以按公式(1)(2)计算间隙水吸引力s(suction)及液桥力F_co。

$ \begin{array}{c} s=P_{\mathrm{a}}-P_{\mathrm{w}}=T\left(\frac{1}{r}-\frac{1}{b}\right) , \end{array} $

(1)

$ F_{\mathrm{co}}=\left(P_{\mathrm{a}}-P_{\mathrm{w}}\right) \pi b^{2}+2 \pi b T。$

(2)

式中：P_a，P_w分别代表孔隙气体压力和孔隙水压力；T代表毛细水表面张力，取室温20 ℃条件下水的表面张力值为0.072 8 N/m; r，b的值可按公式(3)(4)计算：

$ \begin{array}{c} r=R(\sec \delta-1) , \end{array} $

(3)

$ b=R(\tan \delta+1-\sec \delta)。$

(4)

对于半径不同的颗粒，则以等效半径(R₁+R₂)/2替代R代入公式(3)(4)计算。

2 含液桥力的颗粒接触模型的数值实现 2.1 接触绑定模型

颗粒接触点处接触力与相对位移的关系的本构特性见图 2所示。

通过命令指定颗粒接触点处的法向黏结强度和切向黏结强度，其大小可以根据液桥力的计算公式求得。

2.2 用户自定义接触本构模型

如图 3所示，将液桥力简化为线性变化，颗粒生成之后将颗粒接触点处的弹簧黏壶模型转换为自定义的含液桥力的颗粒接触模型，然后进行卸料计算。

考虑间隙水影响的法向接触力为

$ F_{i}^{\mathrm{n}}=K^{\mathrm{n}} U^{\mathrm{n}} n_{i}-F_{\mathrm{cmax}}。$

(5)

切向接触力的计算公式为

$ \Delta F_{i}^{\mathrm{s}}=K^{\mathrm{s}} \Delta U_{i}^{\mathrm{s}}。$

(6)

考虑间隙水影响的切向接触力的极值为

$ F_{\max }^{\mathrm{s}}=\mu\left|F_{i}^{\mathrm{n}}+F_{\mathrm{cmax}}\right|+F_{\mathrm{cw}}, $

(7)

其中F_cw=μF_cmax。

笔者主要探究湿颗粒对筒仓内颗粒流动的影响，考虑FarhangRadjai和Vincent Richefeu^[17]对湿颗粒材料的模拟建议含水量由断裂距离δ_n取值控制，不同的含水量对应不同的δ_n。根据文献[14]给出的断裂距离δ_n和液桥体积的关系为

$ \delta_{n}=(1+0.5 \theta) V^{(1 / 3)}。$

(8)

首先取δ_n为定值0.02 mm，对应的最大液桥力F_cmax由公式(2)计算得到，液桥断裂时的作用力大小设定为0.1F_cmax。

3 筒仓卸料数值算例

筒仓为平面模型，如图 4所示，筒仓内颗粒视作圆球。仓高H=0.5 m，顶部宽W=0.25 m，底部开口直径D和倾角α取60°。

颗粒样本参数参考了文献[15]，如下：颗粒半径为1 mm，筒仓内颗粒数目为1 000，孔隙度n=0.10，颗粒与墙的刚度均取值kⁿ=1.0×10⁶，k^s=1.0×10⁶，其中kⁿ为法向刚度，k^s为切向刚度。颗粒与颗粒间摩擦系数取μ_B=0.5，颗粒与墙之间摩擦系数取μ_W=0.5。黏性系数β_N=β_T=0.02。颗粒密度ρ=2 000 kg/m³，重力加速度取g=9.8 m/s²。

3.1 仓壁压力

筒仓内仓壁压力的变化主要以右侧仓壁为研究对象，将仓壁从下到上10等分，取10个数据提取点，提取点分布如图 5所示，统计计算过程中每点处所受到的压力。筒仓内样本为平均半径1.0mm的颗粒，图 6(a)(b)是对仓壁底端压力提取点1处在一次完整卸料过程中的压力值(时间步长用p表示)。

从图 6中可以得到，仓壁底端的压力峰值发生在卸料过程中，而不是发生在物料堆积阶段。干颗粒弹簧黏壶模型在物料堆积的平衡过程中压力波动较小，开始卸料后波动变大，且峰值大于湿颗粒自定义模型。

3.2 阻塞概率

为模拟阻塞概率，在一次卸料完成或发生阻塞后，按照完全相同的粒度分布在筒仓内重新生成相同数目颗粒样本，再次模拟卸料过程。对任意一组特定的筒仓形状参数或是物性参数，反复模拟卸料过程100次，记录其中发生阻塞的次数N_J，阻塞概率P_J=N_J%。

取样本，仓内颗粒集合体平均半径取1.00 mm，粒径取均匀分布，R_max=1.25×10^-3 m，R_min=0.80×10^-3 m，R_ave=1.00×10^-3 m，其余参数不变。

笔者用颗粒流方法分别做了粒度单一分布和服从均匀分布的颗粒材料阻塞概率P_J与底部开口直径与颗粒平均直径的比值ζ=D/2R_ave的关系，计算结果如图 7所示。

数据分析可知，相同条件下，均匀分布粒径的样本阻塞概率要大于单一粒径分布的样本。在卸料的过程中，大颗粒扮演了重要的角色，对阻塞有重要影响。对于绑定模型，单一粒径分布时湿颗粒阻塞概率要大于干颗粒，但是颗粒粒径均匀分布时，湿颗粒的作用规律不明显, 用自定义接触本构对液桥作用进行模拟，计算结果如图 8所示。

图 8是筒仓中颗粒接触本构分别取干颗粒弹簧黏壶模型和湿颗粒自定义模型时的阻塞概率，其中最大液桥力分别取0.000 5 N和1.000 0 N。湿颗粒自定义模型的阻塞概率与弹簧黏壶模型相比有所提高，但整体趋势一致。最大液桥力取0.000 5 N时与干颗粒的差别不是很大，在4.5＜ζ＜5.5时甚至出现了干颗粒样本的阻塞概率更高的情况，这也与计算样本的局限性有关，颗粒数目较少，循环次数较少。但最大液桥力取1.000 0 N时湿颗粒样本明显比干颗粒样本容易阻塞，反映了一定含液量对颗粒流动性的影响。

4 结论

笔者基于湿颗粒液桥理论模型对含液量较少的湿颗粒材料进行了筒仓卸料的离散元模拟，主要工作和结论如下：

1) 采用接触点绑定和自定义本构两种途径模拟液桥力，对颗粒间含液量较少情况下的液桥作用进行数值实现。

2) 通过对筒仓卸料的数值模拟，比较了干颗粒和湿颗粒在筒仓卸料过程中宏观特性的主要不同，得出仓壁底端的压力峰值发生在卸料过程中，且干颗粒的压力峰值高于湿颗粒压力峰值。

3) 筒仓卸料过程中，在相同条件下，粒径均匀分布的颗粒其阻塞概率要大于单一粒径的颗粒。湿颗粒阻塞概率大于干颗粒，液桥力整体上对颗粒流动有阻碍作用。

本文的研究仍然有一定的局限性，如对筒仓卸料过程中样本颗粒数量偏少，后续研究可以采用更大的样本量进行定量分析，使此模型模拟更为完善。