偏心料斗颗粒速度和体积分数分布自相似特性

Self⁃similar characteristics of particle velocity and volume fraction distribution in eccentric hopper

摘要

基于离散元法对颗粒流在重力驱动下通过偏心料斗出口的流动特性进行了数值模拟分析，分别使用偏心量e（料斗中心与出口中心的距离）和s（料斗左侧墙壁与料斗出口左端之间的距离）表征偏心料斗卸料过程，系统研究了出口偏心位置对出口处颗粒速度和体积分数分布的影响。结果表明，针对不同出口尺寸，采用偏心量e表征卸料过程会导致不同的颗粒流动模式，而采用偏心量s时，卸料过程中颗粒流动模式具有相似性。进一步基于偏心量s的分析表明，颗粒速度和体积分数分布具有自相似特征，其归一化轮廓曲线可以用指数函数描述，指数参数决定曲线轮廓的变化趋势，具有明确的物理意义。结合出口偏心位置与质量流量，提出了矩形偏心料斗的颗粒流量预测公式，通过该公式发现，出口在料斗壁面（s=\mathrm{0}）相对其他位置具有较高的卸料速度。在微观尺度下，研究料斗卸料过程中的颗粒动力学特征和物理规律，有助于优化料斗结构和颗粒流动，并精确预测卸料流量。

Abstract

Based on the discrete element method， the flow characteristics of granular matter through eccentric hopper outlet driven by gravity were numerically investigated. Two different eccentricities e （distance between the center of the hopper and the center of the outlet） and s （distance between the left wall of the hopper and the left end of the hopper outlet） were used to characterize the discharge process of the eccentric hopper， the influence of the eccentricities on the distribution of particle velocity and volume fraction at the orifice was systematically studied. The results show that for different outlet sizes， when the eccentricity e is used to characterize the discharge process， the granular flow exhibits different flow patterns， while the discharge flow patterns are similar when the eccentricity variable s is used. Further， according to the analysis of eccentricity s， the distribution of particle velocity and volume fraction have self-similar characteristics， their normalized contour curves can be well described by the exponential function， and the exponential parameter determines the trend of curve contour change， which has clear physical significance. Combining the profiles of particle velocity and volume fraction， a normalized formula for predicting flow rate of rectangular eccentric hopper is proposed. It is found that it has the maximum discharge rate when the outlet is located at the wall of the hopper （s=\mathrm{0}）. At the microscopic scale， understanding the dynamic characteristics and physical laws of the particle discharging is helpful to optimize the structure of the hopper and particle flow， as well as accurately predict the discharge rate.

Graphical abstract

关键词

农业工程 / 颗粒速度分布 / 体积分数分布 / 颗粒流量预测 / 离散元

Key words

agricultural engineering / particle velocity distribution / volume fraction distribution / discharge rate prediction / discrete element method

0 引 言

颗粒物质是由宏观尺寸的离散粒子组成的体系，广泛存在于工农业生产和自然环境中。一般来说，颗粒物质的存放通常采用料斗，其从料斗中流出的过程称为卸料^［1-5］。颗粒流的许多问题都和颗粒物质在料斗中的流动特性相关^［6，7］，尤其是在农业物料贮藏与排放领域，精准预测颗粒流量已成为当前研究的热点。目前大多数研究都集中在宏观尺度上，如料斗出口尺寸大小^［8］，料斗的几何形状^［9］、外部驱动力等^［10］因素对流量的影响。颗粒的动力学行为也决定着料斗卸料的宏观表现，因此，在微观尺度下研究颗粒的流动状态和分布情况至关重要。另外，由于卸料过程不可见，目前对料斗内部颗粒流动过程和物理规律的认识仍然不足。

基于以上问题，研究人员使用离散元法对料斗卸料进行仿真模拟，以研究某些物理因素对料斗颗粒流动的影响。离散元法是研究散体物料运动规律的一种数值模拟方法，被广泛应用于核能发电、工业加工、农业贮藏与运输、岩土工程等领域^［11-19］。采用离散元法进行数值模拟时，主要通过时步迭代计算区间内颗粒的受力和位移。根据每一时步颗粒之间的相互作用以及牛顿运动定律进行迭代，进而计算颗粒的运动特征。Gao等^［20］使用离散元法对不同长径比椭圆颗粒通过二维矩形料斗排放进行了模拟，结果显示颗粒流在出口处的速度分布与体积分数分布都遵循相同的自相似定律，并计算了椭圆颗粒的流量。Ji等^［10］采用离散元法模拟了外部压力驱动的颗粒流，研究了不同外部压力下的筒仓流量，并且分析了颗粒的速度分布及法向接触力，证明了外部压力对排放速率的影响。Wang等^［21］采用离散元法对水稻种子从三维料斗不同出口位置的排放进行了模拟，分析了颗粒速度分布，提出了料斗临界偏心率以及流量预测模型。Li等^［22］采用离散元法模拟研究了高温气冷堆（HTGR）吸收球系统中带弯曲出口的垂直管的颗粒排放率，分析讨论了管道几何参数、立管内流量波动以及摩擦因数对颗粒流量的影响，并建立了流量预测方程。以上研究中料斗出口均位于中间位置，然而在实际生产过程时，在料斗结构选型和优化设计时，为调整出口位置以适应工农业生产需求，往往采用偏心料斗。此时颗粒流的相关物理量（如速度、体积分数分布）如何变化，其流量计算方法是否与中心对称料斗保持一致尚需进一步研究。针对上述问题，本文基于自由落体理论^{［23，24］}，采用离散元法对偏心料斗卸料过程进行研究。通过改变偏心率分析颗粒在料斗出口的流动特性，并提出偏心料斗流量预测的归一化公式，为料斗的结构选型和设计提供理论依据。

1 离散元模拟

1.1　数学模型

料斗卸料时，会产生料斗壁面与颗粒、颗粒与颗粒之间的接触与碰撞，因此选择合适的接触模型至关重要。本文视散体物料为球形颗粒，假设物料之间没有黏性，颗粒间也不发生相对滑动，因此采用Hertz Mindlin （no slip）接触模型对料斗卸料进行模拟。每个颗粒受到的合力为颗粒相互作用的接触力与重力之和：

m_i\mathit{g}+{\mathit{F}}_i^{\mathrm{c}}=m_i\frac{\mathrm{d}{\mathit{u}}_i}{\mathrm{d}t}（1）

式中：m_i为颗粒i的质量；{\mathit{F}}_i^{\mathrm{c}}为接触力；{\mathit{u}}_i为颗粒i的速度。

颗粒i与颗粒j之间的接触力{\mathit{F}}_{ij}^{\mathrm{c}}可以分为法向分量{\mathit{F}}_{ij}^{\mathrm{n}}和切向分量{\mathit{F}}_{ij}^{\mathrm{t}}。法向力{\mathit{F}}_{ij}^{\mathrm{n}}由赫兹弹性力{\mathit{F}}_{ij}^{\mathrm{n}\mathrm{e}}和黏性力{\mathit{F}}_{ij}^{\mathrm{n}\mathrm{v}}组成^{［10，12，25］}：

{\mathit{F}}_{ij}^{\mathrm{n}\mathrm{e}}=\frac{\mathrm{4}}{\mathrm{3}}{E}^{\mathrm{*}}(R^{\mathrm{*}}{)}^{\frac{\mathrm{1}}{\mathrm{2}}}{{\delta }_{ij}^{\mathrm{n}}}^{\frac{\mathrm{3}}{\mathrm{2}}}（2）

{\mathit{F}}_{ij}^{\mathrm{n}\mathrm{v}}=C_n{\left[\mathrm{8}{m}^{\mathrm{*}}{E}^{\mathrm{*}}{\left(R^{\mathrm{*}}{\delta }_{ij}^{\mathrm{n}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\right]}^{\frac{\mathrm{1}}{\mathrm{2}}}{\mathit{v}}_{ij}^{\mathrm{n}}（3）

式中：{\delta }_{ij}^{\mathrm{n}}为颗粒法向重叠度；C_{\mathrm{n}}为法向阻尼系数；{\mathit{v}}_{ij}^{\mathrm{n}}为法向相对速度；E^{\mathrm{*}}为颗粒等效弹性模量；R^{\mathrm{*}}为等效半径；m^{\mathrm{*}}为等效质量。

\frac{\mathrm{1}}{E^{\mathrm{*}}}=\frac{\mathrm{1}-{\nu }_i^{\mathrm{2}}}{E_i}+\frac{\mathrm{1}-{\nu }_j^{\mathrm{2}}}{E_j}（4）

\frac{\mathrm{1}}{R^{\mathrm{*}}}=\frac{\mathrm{1}}{R_i}+\frac{\mathrm{1}}{R_j}（5）

\frac{\mathrm{1}}{m^{\mathrm{*}}}=\frac{\mathrm{1}}{m_i}+\frac{\mathrm{1}}{m_j}（6）

式中：E_i、E_j分别为粒子i、j的弹性模量；{\nu }_i、{\nu }_j分别为粒子i、j的泊松比；R_i、R_j分别为粒子i、j的半径；m_j为粒子j的质量。

切向力{\mathit{F}}_{ij}^{\mathrm{t}}基于Mohr-Coulomb定律Mindlin理论^［26］，由黏性力{\mathit{F}}_{ij}^{\mathrm{t}\mathrm{v}}和弹性力{\mathit{F}}_{ij}^{\mathrm{t}\mathrm{e}}组成：

\left\{\begin{array}{l}{\mathit{F}}_{ij}^{\mathrm{t}\mathrm{e}}={\mu }_{\mathrm{p},\mathrm{p}}\left|{\mathit{F}}_{ij}^{\mathrm{n}\mathrm{e}}\right|\left(\mathrm{1}-{\left(\mathrm{1}-\frac{\mathrm{m}\mathrm{i}\mathrm{n}({\delta }_{\mathrm{t}},{\delta }_{\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}})}{{\delta }_{\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}}\right)}^{\frac{\mathrm{3}}{\mathrm{2}}}\right)\overline{\mathit{t}}\\ {\mathit{F}}_{ij}^{\mathrm{t}\mathrm{v}}=C_{\mathrm{t}}\left(\mathrm{6}{\mu }_{\mathrm{p},\mathrm{p}}{m}^{\mathrm{*}}\left|F_{ij}^{\mathrm{n}\mathrm{e}}\right|\sqrt[]{\frac{\mathrm{1}-\mathrm{m}\mathrm{i}\mathrm{n}({\delta }_{\mathrm{t}},{\delta }_{\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}})}{{\delta }_{\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}}}\right)v_{ij}^{\mathrm{t}}\end{array}\right.（7）

式中：{\mu }_{\mathrm{p},\mathrm{p}}、C_{\mathrm{t}}分别为滑动摩擦因数、切向阻尼系数；{\delta }_{\mathrm{t}}和v_{ij}^{\mathrm{t}}分别为切线方向的位移和相对速度，{\delta }_{\mathrm{t},\mathrm{m}\mathrm{a}\mathrm{x}}={\mu }_{\mathrm{p},\mathrm{p}}\left(\mathrm{2}-\nu \right)/\mathrm{2}\left(\mathrm{1}-\nu \right)\bullet {\delta }_{ij}^{\mathrm{n}}，\nu为泊松比；\overline{\mathit{t}}为切线方向的单位向量。

当颗粒间存在相对转动时，滚动摩擦力矩{\mathit{M}}_R可表示为：

{\mathit{M}}_R={\mu }_{\mathrm{r}\mathrm{⁃}\mathrm{p},\mathrm{p}}{R}_i\left|{\mathit{F}}_{ij}^{\mathrm{n}}\right|{\widehat{\omega }}_{ij}（8）

式中：{\mu }_{\mathrm{r}\mathrm{⁃}\mathrm{p},\mathrm{p}}为滚动摩擦因数，{\widehat{\omega }}_{ij}为两个接触颗粒的相对角速度，由{\widehat{\omega }}_{ij}={\omega }_{ij}/\left|{\omega }_{ij}\right|给出。

在Hertz-Mindlin接触模型中，粒子间临界时间步长计算公式为^［27］：

t_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{\pi }{r}_{\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\sqrt[]{{\rho }_{\mathrm{p}}/G}}{\mathrm{0.163}\nu +\mathrm{0.876}\mathrm{66}}（9）

式中：r_{\mathrm{m}\mathrm{i}\mathrm{n}}与{\rho }_{\mathrm{p}}分别为颗粒的最小半径和密度；G为剪切模量。

在本文的离散元法模拟中，时间步长取小于\mathrm{0.1}{t}_{\mathrm{m}\mathrm{a}\mathrm{x}}，以避免数值方法的不稳定性。

1.2　物理模型

本文采用上述接触模型对偏心料斗卸料进行仿真模拟。偏心料斗模型如图1所示（初始体积分数{\varphi }_{\mathrm{B}}=\mathrm{0.569}），其中料斗宽度为L，出口尺寸为D。在仿真中，采用两个偏心量来定义出口的位置，分别为e（料斗中心与出口中心的距离）和s（料斗左侧墙壁与料斗出口左端之间的距离），满足s+e=(L-D)/\mathrm{2}。料斗的几何尺寸以及颗粒的仿真参数如表1所示。

整个卸料过程分为加载堆积和卸料两个步骤：初始时，颗粒在料斗中自由落体随机堆积，此时卸料口为关闭状态；当堆积达到稳定状态（t=\mathrm{5}\mathrm{s}）后，卸料口打开，颗粒在重力作用下流出料斗。单一粒径的颗粒容易互相聚集，形成晶格化现象，进而影响流动^［28-31］。然而，本文对10%粒径标准差与单一粒径颗粒在稳定状态下的瞬时流量进行了对比，发现平均流量误差在2%以内，表明粒径离散性对结果的影响可忽略。为方便统计，本文采用单一粒径颗粒进行仿真研究。

2 结果及讨论

图2为料斗出口位于底部最左侧时的颗粒排出过程。在重力驱动下，颗粒流动由出口向周围逐步扩散，同时出口右侧颗粒向左侧挤压汇聚，最终通过出口排出。由图2可知，料斗上方颗粒流动较缓慢，速度向出口方向逐渐增大，并在卸料口处达到最大。由于出口处颗粒受到更大的压力作用且运动约束减小，颗粒以较大的速度从料斗排出；而料斗底部出口对侧的颗粒运动较为缓慢，甚至形成颗粒停滞区。随着时间的推移，大部分颗粒最终流出料斗。

对于偏心料斗卸料过程，分别用偏心量e和s 表征流动状态。图3（a）（b）展现了两种不同的颗粒流动模式。当e=\mathrm{27}\mathrm{m}\mathrm{m}时，图3（a）中颗粒从口径比D/d=\mathrm{6}的料斗排出，颗粒从两侧向出口汇聚；而图3（b）中颗粒从口径比D/d=\mathrm{15}的料斗排出，颗粒均从一侧向出口流动。这表明在相同偏心量e下，颗粒排出的流动模式不同。对比图3（b）和（c），在相同s=\mathrm{0}情况下，颗粒从口径比D/d=\mathrm{15}和D/d=\mathrm{6}的偏心料斗流出时，均向最左侧汇聚且从料斗最左侧流出，其速度变化规律基本一致，流动模式相同。因此，当使用偏心量s来描述颗粒流动时，尽管出口尺寸不同，但都能得到相似的流动状态。

2.1　偏心量\mathit{e}表征颗粒速度和体积分数

基于以上分析，首先讨论以偏心量e表征料斗排放颗粒的速度与体积分数分布。如图4所示，在料斗出口处测量颗粒的速度与体积分数。其中，虚线为颗粒速度测量位置，方框为颗粒体积分数测量区域，通过稳定流动时间段内颗粒通过的平均数量进行计算。

当颗粒通过料斗出口时，记录其位置坐标和速度垂直分量v。图5（a）展示了出口位于料斗中心处时颗粒速度与颗粒位置的变化关系，出口尺寸影响颗粒速度分布，其曲线呈轴对称状。随着偏心量e的增加，颗粒速度分布发生明显偏移，最大值不再位于料斗出口中心处，而是向e的改变方向偏移，如图5（b）（c）所示。另外，研究发现不同口径比的颗粒速度均呈抛物线分布，具有一定的相似性。因此，本文对速度分布进行归一化处理。根据自由落体理论^［24］和颗粒速度分布特点，采用式（10）对颗粒速度与位置进行归一化处理：

\overline{v}=\frac{v}{\sqrt[]{\mathrm{2}gR}}=\epsilon {\left[\mathrm{1}-{\left(\frac{x-e}{R}\right)}^{\mathrm{2}}\right]}^m（10）

式中：v为颗粒速度；g为重力加速度；R为出口半径（R=D/2）；\left(x-e\right)/R为颗粒水平位置在出口处的归一化；\epsilon和m为拟合参数。

速度归一化结果如图5（d）~（f）所示，其中 图5（d）（e）中颗粒速度数据点集中在同一条曲线上，而图5（f）中，颗粒速度归一化之后，口径比D/d=\mathrm{10}和D/d=\mathrm{13}的颗粒速度散点落在归一化曲线上，而其他出口的颗粒速度均在出口偏心一侧发散，尤其是口径比D/d=\mathrm{15}的颗粒速度与归一化曲线偏移较大。通过式（10）得到的拟合参数如表2所示。参数\epsilon在e=\mathrm{27}\mathrm{m}\mathrm{m}时的值大于e=\mathrm{0}时的值，表明不同偏心量e下的速度无法用该公式表征。

对出口处颗粒的体积分数与位置进行分析，结果如图6（a）~（c）所示。在e=\mathrm{0}时，颗粒体积分数分布保持轴对称；随着出口向左偏移至e=\mathrm{27}\mathrm{m}\mathrm{m}时，颗粒体积分数分布呈现左高右低的特征，但变化趋势不显著。采用与速度曲线归一化相同的方法对颗粒体积分数进行归一化处理：

\frac{\varphi }{{\varphi }_{\mathrm{c}}}={\left[\mathrm{1}-{\left(\frac{x-e}{R}\right)}^{\mathrm{2}}\right]}^n（11）

式中：{\varphi }_{\mathrm{c}}为出口中心颗粒的体积分数；n为颗粒体积分数分布平滑度参数，其值见表2。

随着偏心量e的增大，拟合参数n呈现微小增长趋势，表明在偏心量e表征下，颗粒体积分数受其影响较小。图6（d）~（f）为采用式（11）拟合的归一化曲线。在e=\mathrm{27}\mathrm{m}\mathrm{m}时，颗粒体积分数并没有出现发散现象，表明不同偏心量e下的颗粒体积分数分布仍符合自相似定律。

当e=\mathrm{0}时，图6（d）中得到的颗粒体积分数指数n=\mathrm{0.386}， 比其他研究者得出的实验值略大（Zhou等^［32］得出n=\mathrm{0.19}，Janda等^［33］得出n=\mathrm{0.22}）。研究发现，料斗壁厚和颗粒-壁面摩擦因数会影响颗粒体积分数。因此，本文研究了3种壁厚（Z=\mathrm{1.01}d、\mathrm{1.05}d、\mathrm{1.10}d）和3种颗粒-壁面滑动摩擦因数（{\mu }_{\mathrm{p},\mathrm{w}}=\mathrm{0.01}、0.2、0.5）对颗粒流动的影响，结果如图7所示。结果表明，料斗质量流量W随Z和{\mu }_{\mathrm{p},\mathrm{w}}的增大而减小（见图7（a）~（b）），进而影响出口处颗粒速度和体积分数分布。然而，当Z=\mathrm{0.01}d时，n值不再改变，如图7（c）所示；当{\mu }_{\mathrm{p},\mathrm{w}}=\mathrm{0.01}时，颗粒体积分数指数n=\mathrm{0.2}，此时的结果与实验值接近，如图7（d）所示。由此可知，调节颗粒-壁面滑动摩擦因数可得到与实验一致的结果，与Blanco-Rodríguez等^［34］得出的结论一致。

2.2　偏心量\mathit{s}表征颗粒速度与体积分数

对出口处颗粒速度与位置的关系（用偏心量s表征）进行分析，发现所有s值下各出口位置的颗粒速度变化趋势均相同，其速度呈抛物线分布，如图8（a）~（c）所示。当出口位于料斗最左侧（s=\mathrm{0}）时，出口中心左侧颗粒速度均大于右侧；随着出口向料斗中心处偏移，其最大速度接近出口中心。在极端偏心情况下（s=\mathrm{0}），颗粒速度剖面表现出不对称性，最大速度不再位于料斗出口中心，式（10）已无法表征其特征，需要引入额外项描述偏心的影响，进而对颗粒速度分布进行归一化：

\overline{v}={\alpha }_{\mathrm{1}}\left(\frac{x-e}{R}\right)+{\beta }_{\mathrm{1}}{\left[\mathrm{1}-{\left(\frac{x-e}{R}\right)}^{\mathrm{2}}\right]}^a（12）

式中：{\alpha }_{\mathrm{1}}、{\beta }_{\mathrm{1}}均为拟合参数；a为颗粒归一化指数，取值见表3。{\alpha }_{\mathrm{1}}的平均值为1.3，{\alpha }_{\mathrm{1}}\left[(x-e)/R\right]项可以描述偏心量s对速度分布的影响，如图8（d）~（f）所示。所有颗粒速度散点均分布在归一化曲线上，当偏心量s从0开始增加时，料斗出口在边界（s=\mathrm{0}）处的速度归一化曲线与其他偏心量s下的曲线不同。图8（d）显示，颗粒最大速度不再位于出口中心，而是向左偏移（出口位于料斗最左侧），且其值大于其他偏心量下的最大速度。

由于偏心量s的增加，速度分布由不对称向对称过渡，中心两侧颗粒最大速度向料斗出口中心处颗粒靠近。参数{\alpha }_{\mathrm{1}}\left[(x-e)/R\right]中的{\alpha }_{\mathrm{1}}调节偏心料斗下颗粒速度归一化曲线与最大速度的偏移，指数a决定速度归一化曲线的平滑度。这两个参数均与偏心量s有关，因此综合考虑偏心量s、料斗宽度L和出口尺寸D的影响，对偏心量s进行无量纲化处理得到\mathrm{2}s/\left(L-D\right)，结果如图9所示。当出口从料斗边界（s=\mathrm{0}）向中心移动时，\mathrm{2}s/\left(L-D\right)增加，{\alpha }_{\mathrm{1}}\left[(x-e)/R\right]的值减小，速度曲线的不对称性逐渐减弱（{\alpha }_{\mathrm{1}}\to \mathrm{0}）。当出口逐渐接近料斗中心位置时，{\alpha }_{\mathrm{1}}呈指数增长，从一个负值趋近于0，其变化规律符合关系式：{\alpha }_{\mathrm{1}}=-\mathrm{0.3}{\mathrm{e}}^{-\mathrm{5.25}\left[\mathrm{2}s/\left(L-D\right)\right]}。图9（b）显示了指数a随\mathrm{2}s/\left(L-D\right)的变化，其变化趋势与{\alpha }_{\mathrm{1}}相似，最终趋近于0.39，可表示为a=\mathrm{0.39}-\mathrm{0.12}{\mathrm{e}}^{-\mathrm{3.74}\left[\mathrm{2}s/\left(L-D\right)\right]}。偏心量s=0时的参数{\alpha }_{\mathrm{1}}和a与其他偏心量下的值差异显著。尤其是附加项{\alpha }_{\mathrm{1}}\left[(x-e)/R\right]中的{\alpha }_{\mathrm{1}}在偏心量s=\mathrm{0}处值得探究：由于料斗出口位于边界处，出口处料斗壁面滑移导致颗粒速度的增加以及最大速度偏移，从而产生不对称分布。而在其他偏心情况下，颗粒速度分布均趋于对称，最大速度均分布在料斗出口中心，{\alpha }_{\mathrm{1}}均趋近于0。指数a反映出口从偏心料斗向中心料斗过渡时，颗粒分布平滑度呈指数增长，并趋近于中心料斗的值。结果表明，当\mathrm{2}s/\left(L-D\right)从0增加至1时，速度分布对称轴由料斗边界向中心移动，颗粒速度分布曲线由平滑转变为“弯曲”，表示颗粒出口两侧的速度开始减小，证实偏心量s不仅会影响颗粒最大速度与速度分布，也会改变出口中心两侧的速度大小。

以偏心量s表征颗粒体积分数，图10展示了颗粒通过偏心出口的体积分数剖面及其归一化曲线。图10（a）中，当s=\mathrm{0}时，颗粒体积分数分布呈现明显的不对称性。随着s逐渐增大，所有出口尺寸下的体积分数剖面对称性逐渐恢复，见图10（b）和（c）。为了验证体积分数分布的自相似特性，采用出口处最大体积分数{\varphi }_{\mathrm{c}}进行归一化处理。与预期一致，这些点在出口中心和边界均能很好地拟合为一条拟合曲线（图10（d）~（f））。采用与速度曲线归一化相似的方法，体积分数剖面可表示为：

\frac{\varphi }{{\varphi }_{\mathrm{c}}}={\alpha }_{\mathrm{2}}\left(\frac{x-e}{R}\right)+{\beta }_{\mathrm{2}}{\left[\mathrm{1}-{\left(\frac{x-e}{R}\right)}^{\mathrm{2}}\right]}^b（13）

式中：{\alpha }_{\mathrm{2}}和{\beta }_{\mathrm{2}}为拟合参数；指数b决定剖面形状，其参数值见表3。

在图10（d）中，所有体积分数散点通过式（13）归一化为一条曲线，相比偏心量e的体积分数归一化结果，并没有出现发散现象。图10（e）（f）中，随着s的增大，体积分数归一化曲线从向左偏移过渡至出口中心处对称，表明{\alpha }_{\mathrm{2}}\left[(x-e)/R\right]参数项的作用逐渐减弱。式（13）可以对偏心量s表征的偏心料斗出口处体积分数分布精准拟合，同时也考虑了偏心量s的影响。

如图11所示，式（13）中的参数{\alpha }_{\mathrm{2}}和b均可以由无量纲参数\mathrm{2}s/\left(L-D\right)描述。图11（a）中，{\alpha }_{\mathrm{2}}与无量纲参数\mathrm{2}s/\left(L-D\right)呈指数递增关系。偏心量s=\mathrm{0}时，参数{\alpha }_{\mathrm{2}}小于其他偏心量下的参数值，其分布向左偏移，呈现不对称性，此时参数项{\alpha }_{\mathrm{2}}\left[(x-e)/R\right]达到最大值，对颗粒体积分数分布影响最大，这是由于料斗壁面的滑移对颗粒流动有较大影响。但是这种影响有限，当出口接近料斗中心时，参数值{\alpha }_{\mathrm{2}}趋近于0，颗粒体积分数分布恢复对称特性，偏心量s对体积分数分布的影响消失。参数b被用来描述颗粒体积分数分布的平滑度。在图11（b）中，该参数值在s=\mathrm{0}时小于其他偏心情况；当偏心量s为6~30 mm时，b在0.42左右波动。在此范围内，颗粒体积分数曲线分布几乎具有相同的曲率。

2.3　偏心料斗的质量流量预测模型

为建立偏心料斗的流量预测模型，将质量流量W(s)与偏心量s无量纲化，得到图12。所有质量流量数据均由仿真后处理中的质量流量监测器获取，该监测器置于出口下方。其中，W_{\mathrm{0}}为料斗无偏心时的流量（即e=\mathrm{0}），W(s)为不同出口尺寸料斗在不同偏心量s下的质量流量。通过图中的拟合曲线来预测不同偏心量s下的流量，其公式为：

\frac{W(s)}{W_{\mathrm{0}}}=A{\mathrm{e}}^{-\frac{\mathrm{2}s/(L-D)}{B}}+C（14）

式中：A、B、C为拟合参数。

当料斗为小出口（D/d=\mathrm{6}）时，A=1.086，B=0.325，C=0.128；料斗出口为其他尺寸时，A=0.994，B=0.264，C=0.143。

在图12中，随着偏心量增大，各出口流量逐渐接近于非偏心料斗流量。对于D/d=\mathrm{6}的出口，其流量比均大于其他出口尺寸。在出口极度偏心（s=\mathrm{0}）情况下，颗粒在壁面处的碰撞加剧且受到壁面滑移作用，其最大速度增大并且向偏心方向偏移，从而增大颗粒排出流量。在其他偏心情况下，各出口流量比均接近1，这是由于非极度偏心时颗粒最大速度波动较小且偏移不明显，对料斗流量的影响较弱。

当\mathrm{2}s/\left(L-D\right)=\mathrm{0}时，料斗流量比是非偏心出口的1.2~1.4倍；在0.2<\mathrm{2}s/\left(L-D\right)<1.0范围内，料斗的流量几乎与非偏心料斗的流量相同。对于相同偏心率\mathrm{2}s/\left(L-D\right)，出口较小的料斗可提供更大的流量增量。这一发现可为优化出口位置以控制颗粒流速提供理论依据。

3 结 论

（1）采用偏心量s表征出口位置时，颗粒在不同出口尺寸下呈现出相似的流动模式，而偏心量e则导致不同的流动模式，尤其是在出口位于料斗壁面时。由此判断，采用s表征料斗出口偏心距离对颗粒流动的影响更为合适。

（2）针对矩形偏心料斗卸料过程，采用s表征偏心距离时，出口处颗粒速度与体积分数分布曲线仍然符合自相似定律，这与非偏心料斗的分布特性一致。

（3）通过将偏心量s与质量流量无量纲化，提出了一项考虑偏心位置的归一化公式来预测偏心料斗的流量。结果表明，料斗出口在壁面处（s=\mathrm{0}）的质量流量始终高于其他位置。

参考文献

[1] Lu F, Zhang C, Wei F. Compressibility of granular fluids [J]. Physics of Fluids, 2022, 34(3): No. 033307.

[2] Huang Y, Zhu C, Xiang X. Granular flow under microgravity: A preliminary review[J]. Microgravity Science and Technology, 2014, 26(2): 131-138.

[3] Meng F, Pang M, Liu K. Constitutive relation and mechanical mechanism analysis in granular lubrication [J]. Industrial Lubrication and Tribology, 2020, 72(5): 621-628.

[4] Campbell C S. Granular material flows—An overview[J]. Powder Technology, 2006, 162(3): 208-229.

[5] Zheng N, Zhu H W, Huang D C, et al. The effects of physical factors on the granular hopper flow[J]. Scientia Sinica: Physica, Mechanica et Astronomica, 2020, 50(9): No.090005.

[6] 孟凡净, 刘焜. 密集剪切颗粒流中速度波动和自扩散特性的离散元模拟[J]. 物理学报, 2014, 36(13): 262-270.

[7] Meng Fan-jing, Liu Kun. Velocity fluctuation and self diffusion character in a dense granular sheared flow studied by discrete element method[J]. Acta Physica Sinica, 2014,63(13): 262-270.

[8] 胡国琦, 张训生, 鲍德松, 等. 二维颗粒流通道宽度效应的分子动力学模拟[J]. 物理学报, 2004, 53(12):4277-4281.

[9] Hu Guo-qi, Zhang Xun-sheng, Bao De-song, et al. The molecular dynamics simulation of the effect of channel width on two-dimensional granular flow[J]. Acta Physica Sinica, 2004, 53(12): 4277-4281.

[10] Anand A, Curtis J S, Wassgren C R, et al. Predicting discharge dynamics from a rectangular hopper using the discrete element method (DEM)[J]. Chemical Engineering Science, 2008, 63(24): 5821-5830.

[11] Beverloo W A, Leniger H A, Van de velde J. The flow of granular solids through orifices[J]. Chemical Engineering Science, 1961, 15: 260-269.

[12] Ji S, Wang S, Peng Z. Influence of external pressure on granular flow in a cylindrical silo based on discrete element method[J]. Powder Technology, 2019, 356: 702-714.

[13] Ma H, Zhao Y. Influence of external pressure on granular flow in a cylindrical silo based on discrete element method[J]. Chemical Engineering Science, 2017, 172: 636-651.

[14] Wang S, Zhuravkov M, Ji S. Granular flow of cylinder-like particles in a cylindrical hopper under external pressure based on DEM simulations[J]. Soft Matter, 2020, 16(33): 7760-7777.

[15] Wu S M, Zhu H P, Yu A B, et al. DEM simulation of granular flows on a heap[J]. AIP Conference Proceedings, 2009, 1145: 621-624.

[16] Zhou G, Ng C. Numerical investigation of reverse segregation in debris flows by DEM[J]. Granular Matter, 2010, 12(5): 507-516.

[17] Fan J, Luu L H, Philippe P, et al. Discharge rate characterization for submerged grains flowing through a hopper using DEM-LBM simulations[J]. Powder Technology, 2022, 404: No.117421.

[18] Fan J, Luu L H, Noury G, et al. DEM-LBM numerical modeling of submerged cohesive granular discharges[J]. Granular Matter, 2020, 22(3): No.66.

[19] 原方,杜乾,徐志军,等. 基于卸料流态模拟与观测的储粮仓壁动态压力增大机理研究[J]. 农业工程学报, 2019, 35(5): 286-293.

[20] Yuan Fang, Du Qian, Xu Zhi-jun, et al. Mechanism of dynamic pressure increase of grain silo wall based on simulation of discharge flow states[J]. Transactions of the Chinese Society of Agricultural Engineering, 2019, 35(5): 286-293.

[21] 徐贵玲, 李梦慧, 卢平, 等. 烘培生物质与煤不同配比混合物的流动及下料特性[J]. 农业工程学报, 2018, 34(1): 186-192.

[22] Xu Gui-ling, Li Meng-hui, Lu Ping, et al. Flowability and discharge characteristics of mixtures with different ratio of torrefied biomass and pulverized coal[J]. Transactions of the Chinese Society of Agricultural Engineering, 2018, 34(1): 186-192.

[23] 刘克瑾, 姚辉江, 袁铮. 平底筒仓中心卸料流态演化过程及仓壁压力波动性[J]. 农业工程学报,2023,39(15):15-24.

[24] Liu Ke-jin, Yao Hui-jiang, Yuan Zheng. Evolution process of center discharge flow pattern and pressure fluctuation of silo wall in a flat-bottomed silo[J]. Transactions of the Chinese Society of Agricultural Engineering, 2023, 39(15): 15-24.

[25] Gao Q, Chen Y, Zhao C. Self-similarity of density and velocity profiles in a 2D hopper flow of elliptical particles: discrete element simulation[J]. Powder Technology, 2022, 402: No.117338.

[26] Wang Y, Jia F, Zhang J, et al. Model construction method of discharge rate of eccentric silo[J]. Powder Technology, 2022, 405: No.117555.

[27] Li T, Zhang H, Liu M, et al. DEM study of granular discharge rate through a vertical pipe with a bend outlet in small absorber sphere system[J]. Nuclear Engineering and Design, 2017, 314: 1-10.

[28] Hilton J E, Cleary P W. Granular flow during hopper discharge[J]. Phys Rev E Stat Nonlin Soft Matter Phys, 2011, 84(1): No.011307.

[29] Tang Y, Chan D H, Zhu D Z. The modeling of free-fall arch formation in granular flow through an aperture[J]. Frontiers in Physics, 2022, 10: No.963495.

[30] Zhang Z, Liu Y, Zheng B, et al. Discharge characteristics of binary particles in a rectangular hopper with inclined bottom[J]. Computational Particle Mechanics, 2021, 8(2): 315-324.

[31] Cundall P, Strack O. A discrete numerical model for granular assemblies[J]. Géotechnique, 1979, 29(1): 47-65.

[32] Di Renzo A, Paolo Di Maio F. An improved integral non-linear model for the contact of particles in distinct element simulations[J]. Chemical Engineering Science, 2005, 60(5): 1303-1312.

[33] Daniels K E, Behringe R P. Hysteresis and competition between disorder and crystallization in sheared and vibrated granular flow[J]. Physical Review Letters, 2005, 94(16): No.168001.

[34] Jing L, Yang G C, Kwok C Y, et al. Dynamics and scaling laws of underwater granular collapse with varying aspect ratios[J]. Physical Review E, 2018, 98(4): No.042901.

[35] Yang G C, Jing L, Kwok C Y, et al. Pore-scale simulation of immersed granular collapse: implications to submarine landslides[J]. Journal of Geophysical Research: Earth Surface, 2020, 125(1): e2019JF005044.

[36] Yang G C, Jing L, Kwok C Y, et al. Size effects in underwater granular collapses: Experiments and coupled lattice Boltzmann and discrete element method simulations[J]. Physical Review Fluids, 2021, 6(11): No.114302.

[37] Zhou Y, Ruyer P, Aussillous P. Discharge flow of a bidisperse granular media from a silo: discrete particle simulations[J]. Physical Review E, 2015, 92(6): No.062204.

[38] Janda A, Zuriguel I, Maza D. Flow rate of particles through apertures obtained from self-similar density and velocity profiles [J]. Physical Review Letters, 2012, 108(24): No.248001.

[39] Blanco-Rodríguez R, Hidalgo R C, Pérez-Ángel G, et al. Critical numerical analysis of quasi-two-dimensional silo-hopper discharging[J]. Granular Matter, 2021, 23(4): 86-96.

基金资助

吉林省科技厅重点研发项目(20240303040NC)

国家自然科学基金重点项目(52130001)