摘要:应用基于局部本构理论的连续数值模拟方法, 研究出口在底部和侧面的颗粒物在类三维矩形容器内的卸载现象. 重点是容器厚度W和出口高度D对颗粒物压强与速度的影响. 受力分析和数值模拟结果均表明, 距离出口较近区域的颗粒物压强与W及D呈现如下相关性: 当D/W足够小时, 压强只与D相关; 当D/W足够大时, 压强只与W相关. 且出口在底部和侧面时均有上述结果. 模拟结果还显示, 当出口在底部时, 对于模拟中所有D/W值, 出口中心处法向速度只和D相关; 当出口在侧面时, 颗粒物出口中心处法向速度则与压强变化规律一致. 由此可见, 出口处的压强并不控制颗粒物的出口法向速度. 另外, 与出口在侧面相比, 出口在底部时, 造成流量相关性规律改变的D/W临界值较大, 一般实际情况无法满足, 因此出口中心处法向速度只与D相关, 始终满足Beverloo定律.
关键词: 局部本构关系/
连续数值模拟/
颗粒物质/
筒仓卸载

English Abstract

全文HTML

--> --> -->

科学界通常认为颗粒介质(granular media)是大于100 μm的颗粒集合体^[1], 其布朗效应和颗粒物质之间的凝聚力相较于重力及碰撞力可以忽略. 颗粒介质在自然界以及工业生产中广泛存在, 因此对颗粒体系的深入研究将会对全球工业与经济的发展有极大助益^[2]. 早期颗粒流的研究手段主要是实验和离散颗粒介质的数值模拟(DEM). 离散元数值模拟方法主要是考虑每个个体颗粒物质之间的碰撞^[3]. 该方法发展已较为成熟, 但计算速度慢. 近年来, 法国颗粒物质研究组提出的局部本构理论^[4]得到了广泛使用, 这一局部流变理论使得将流体力学方法应用于颗粒流模拟成为可能, 这极大程度地降低了运算时间. 另外, 该理论加深了人们对颗粒流内在物理机理的理解. 该方法主要描述的是与固体态相差甚远的液体状态, 其本质是将颗粒物质看成带有摩擦性质的非牛顿流体. 之后Jop等^[5]将此局部流变学从二维向量形式推广至三维张量形式. 该理论被证实能够应用在许多密集颗粒流构型当中, 例如颗粒物坍塌^[6]、筒仓卸载^[7,8]等. 但是存在一定的局限性, 第一个局限性是无法描述颗粒物开始流动和停止流动的现象, 第二个局限性是无法描述堵塞状态的运动^[9].
颗粒物在筒仓或沙漏中的流动是一个经典的颗粒介质基础研究现象. 人们很早就发现沙漏比水漏作为计时器更加简便和准确. 对于出口在底部的筒仓卸载问题已有大量的研究. 当容器尺寸足够大时, 颗粒物的流量始终保持恒定, 并且只与出口尺寸有关, 这是它的重要特征, 但其原因和物理机理仍有待研究. 普遍为人们所接受的解释是Janssen效应^[10,11], 即容器底部压强达到一定值后不再随颗粒高度增加而增加. 这是由于颗粒与筒壁间存在着摩擦力, 筒壁承担了部分颗粒重量. 由于筒仓出口处的压强与颗粒物高度无关, 因此流量也与高度无关. 然而近来一些实验证明颗粒物的流量与出口处的压强不具有相关性^[12,13]. 因此Janssen效应并不能完整地解释此现象. 流量Q随孔径D变化的经验公式, 早在60年代就由Beverloo 等^[14]利用量纲分析和实验方法给出, 并被大量文献证实相当精确可靠. 近几年人们主要开展了筒仓倾斜情况下的研究^[15—17], 即出口方向与重力不平行, 研究发现流量也保持恒定, 并得到了与倾斜角度相关的流量经验公式. 此外Zhou等^[18]研究了出口位置在侧面(即与地面垂直, 倾角为90°)的颗粒物卸载情况(实验所用筒仓如图1所示), 该工作重点研究了容器厚度W和出口高度D对颗粒物流量的影响. 该工作表明, 出口在侧面的情况下, 颗粒物流量仍旧保持恒定, 并且作者发现了与出口在底部情况不同的流量规律. 当D/W 足够小时, 流量规律与出口在底部时的情况一样, 满足Beverloo经验公式, 即流量$ Q\propto S\sqrt{gD} $, 其中$ S $为出口面积; 而当$ D/W $足够大时, 流量规律则为$ Q \propto S\sqrt{gW} $. 可见, 在容器厚度很小时, 卸载流量呈现了与出口在底部情况下不同的规律. 作者使用上述基于局部本构模型的连续模拟方法得到了与实验结果相同的流量规律, 并且证明了导致该结果的主要原因是出口处流线的倾斜角受前后容器壁的摩擦力影响很大. 可见, 上述连续模拟方法能够较好地模拟筒仓卸载实验. 本文将在已有工作的基础上, 利用上述基于本构模型的连续数值模拟方法深入分析筒仓卸载现象. 首先, 分析距离出口近处压强和速度与W及D之间的关系. 进而, 由此研究压强与速度之间的关系, 并利用上述数值模拟解释为何出口在底部时流量规律始终为$ Q\propto S\sqrt{gD} $.
图 1 Zhou等实验所用矩形筒仓示意图(取自文献[18])
Figure1. Schematic apparatus of the rectangular silo used by Zhou et al. (extracted from Ref.[18]).

本文采用基于局部本构模型的连续数值模拟方法, 即将颗粒物看成非牛顿流体, 并且考虑颗粒物之间的摩擦力对颗粒流动的影响. 非牛顿不可压缩的纳维-斯托克斯公式可以写为:

$ {\nabla} \cdot {u} = 0, $

$ \rho \left( \frac{\partial {u}}{\partial t} + {u} \cdot {\nabla} {u} \right) = - {\nabla} p + {\nabla} \cdot (2 \eta {D}) + \rho {g} + {f}_{\rm w}, $

其中$ {D} $是应变速率张量$ (\nabla{u} +\nabla{u}^{\rm T}) / 2 $. 根据文献[5]中提出的理论, $ \eta $为有效黏性系数, 它与摩擦系数 $ \mu(I) $以及局部压强$ p $相关. $ D_2 $为应变速率张量的第二不变量, $ D_2 = \sqrt{D_{ij}D_{ij}} $).

$\begin{split}&\eta = \frac{\mu(I)}{\sqrt{2}D_2}p\; , I = \frac{d\sqrt{2}D_2}{\sqrt{p/\rho}}, \\&\quad \text{另外}, \mu(I) = \mu_{\rm s} + \frac{\Delta \mu}{I_0/I+1}. \end{split}$

取流变学常数$ \mu_{\rm s} = 0.4 $, $ \Delta \mu = 0.28 $和 $ I_0 = 0.4 $, 但不考虑文献[5]提出的体积分数$ \phi $与$ I $的关系式, 此处取$ \phi = 1 $. 方程(2)中$ {f}_{\rm w} $代表前后两个壁面的总体摩擦力. 本文实际采用的是二维模拟, 为考虑容器前后壁对颗粒物的摩擦以及模拟三维效果, 像赫尔-肖(Hele-Shaw)流动一样, 将动量方程沿着$ y $方向进行平均^[19,20]. 根据颗粒流的特性, 假设速度沿着$ y $方向几乎保持不变, 因此动量方程中的非线性项前的系数为$ 1 $. 而黏滞力沿着$ y $方向的积分之和则等于前后壁对颗粒物的摩擦力, 假设该力为库仑力(在每个壁上值为$ -\mu_{\rm w} p $). 该力的方向应当与颗粒物运动方向相反, 因此前后两个壁面的总体摩擦力大小为

$ {{f}}_{\rm w} = -2 \frac{\mu_{\rm w} {\rm p}}{W} \frac{{u}}{ |{u}|}. $

取$ \mu_{\rm w} = 0.1 $, 通过改变容器宽度$ W $的值来改变摩擦力的大小. 可见本文中容器厚度$ W $的影响等效于容器前后壁摩擦力大小的影响. 本文利用上述类三维方法模拟图1所示实验, 具体计算方法参见文献[18]和文献[21]. 在固体边界加上非滑移(no-slip)的边界条件, 在出口处压强设为零. 数值模拟中容器及颗粒物的尺寸如表1所示(容器相关尺寸参数的定义见图1).

$h_p$	4L
$D$	[0.3125, 0.40125, 0.4375, 0.5, 0.5938, 0.625, 0.6562, 0.6875, 0.75]L
$W$	[0.16, 0.2, 0.25, 0.5, 1, 2]L
$d$	$L/90$

表1容器及颗粒物的尺寸
Table1.Size of silo and particle.

2

-->

3.1.出口在底部

选取筒仓卸载过程中流量稳定情况下的某一时刻$ t/\sqrt{L/g} = 4 $(其中$ L $为容器宽度), 该时刻颗粒物压强场、速度场以及无量纲数$ I $的分布在图2中画出. 从容器顶部至底部看压强场云图, 压强变化呈现如下几个特点: 1) 在颗粒物球床较高的位置($ z \approx 3L $), 颗粒物压强随着$ z $降低增大; 2) 很快压强不再升高, 而是在$ 1<z/L<2 $区域几乎保持不变; 3) 对于$ z<L $区域, 压强迅速下降直至出口处压强为$ 0 $; 4) 靠近容器壁处压强比中心处更大, 可见容器壁承担了部分颗粒物的重量, 体现了Janssen 效应. 从速度场云图可以看出, 在$ z>L $的区域内, 竖直方向速度大小几乎相等, 而位于容器底部, 在尺寸约为$ D $的区域内, 速度随着$ z $减小显著增大. 另外, 颗粒物停滞区域(即速度大小为$ 0 $)的边界线与竖直方向所成夹角几乎保持不变. 从无量纲常数$ I $的分布图可以看出, 其数值在整个容器区域内都比较小, 沿着容器壁以及停滞区域的颗粒物剪切力较大, 此时$ I $的值也较大(浅蓝色区域), 但$ I $值最大的区域是颗粒物加速度较大并且压强较小的出口处. $ I $的最大值为$ 0.2 $, 表示临近出口处颗粒物仍旧为密集颗粒流.
图 2 $ D = 0.3125L $以及$ W = 0.25L $情况下$ t/\sqrt{L/g} = 4 $时刻连续数值模拟结果, 从左至右: 容器内压强与其最大值之比; 竖直方向速度与其最大值之比, 其中黑色实线表示颗粒物流线; 无量纲常数$ I = d\sqrt{2}D_2/(\sqrt{p/\rho}) $
Figure2. Continuum simulation results with $ D = 0.3125L $ and $ W = 0.25L $ at time $ t/\sqrt{L/g} = 4 $, from the left to the right: pressure $ p^{\rm p} $ normalized by it's maximum value within the silo; the vertical velocity $ v^{\rm p} $ normalized by it's maximum value within the silo, the black lines represent the streamlines; dimensionless number $ I = d\sqrt{2}D_2/(\sqrt{p/\rho}) $.

可见在靠近出口处, 动能项开始逐渐增大并且可以部分补偿压强降低, 但由于摩擦力的存在, 二者之间的关系并不满足伯努利方程. 接下来重点分析距离出口较近处的颗粒物压强$ p^p $和速度$ v^p $与容器厚度$ W $及出口高度$ D $的关系.
3 -->

3.1.1.压强结果分析

由前面分析可知, 颗粒物在距离出口近处压强和速度变化明显, 在远离出口的区域内压强变化较小, 速度基本保持恒定, 因此本文主要对距离出口近的区域进行分析. 下面首先取该区域一个微元薄层进行受力分析. 定义一个高度为$ \partial \lambda $的微元薄层区域(见图3中的蓝色虚线区域), 且若取容器后壁面出口中点处为坐标原点, 该区域坐标满足$ x\in[-D/2, D/2] $, $ y \in [0, W] $, $ z \in [\lambda, \lambda+\partial \lambda] $. 假设这个区域的边界处所受摩擦力与重力平衡, 可得:

$ 2 \bar{F}_{\rm g}+ 2 \bar{F}_{\rm w} + P = 0, $

其中$ \bar{F}_{\rm g} $代表颗粒物之间的摩擦力(图3中实心箭头), $ \bar{F}_{\rm w} $代表颗粒物与前后容器壁之间的摩擦力(图3中空心箭头), $ P $为重量. 假设$ p^{\rm p} $在这个区域内几乎无变化.
图 3 出口在底部情况下容器内不同区域受力图
Figure3. Force diagram of different zones within the silo for the case with orifice at the bottom.

在$ y $和$ z $方向上积分, 得到$ \bar{F}_{\rm g} $为

$\begin{split} \bar{F}_{\rm g} = &\int_{y = 0}^{W} \int_{z = \lambda}^{\lambda + \partial \lambda}\mu(I) p^{\rm p}(x = D/2)\\&\times \partial y \partial z \propto \mu_{\rm s} p^{\rm p}(\lambda)W \partial \lambda \; ;\end{split}$

在$ x $和$ z $方向积分, 可得$ \bar{F}_{\rm w} $为

$ \begin{split}\bar{F}_{\rm w} =\;& \int_{x = -D/2}^{D/2} \int_{z = \lambda}^{\lambda + \partial \lambda}\mu_{\rm w} p^{\rm p}(y = 0)\\&\times \partial x\partial z\propto \mu_{\rm w} p^{\rm p}(\lambda)D \partial \lambda \; .\end{split} $

由于重量$ P = \rho g D \partial\lambda W $, 取常数$ \alpha_1 $和$ \alpha_2 $, 方程(5)可写成

$ 2 \alpha_1 \mu_{\rm s} p^{\rm p}(\lambda) W \partial \lambda + 2 \alpha_2 \mu_{\rm w} p^{\rm p} (\lambda) D \partial \lambda = \rho g D \partial \lambda W. $

对颗粒物压强进行无量纲化操作, 可得到:

$ \frac{p^{\rm p}(\lambda) }{\rho g D} = \frac{ \gamma_1 }{1 + \gamma_2 D/W}, $

其中$ \gamma_1 = {1}/({2 \alpha_1 \mu_{\rm s}}) $, $ \gamma_2 = ({\alpha_2 \mu_{\rm w}})/({\alpha_1 \mu_{\rm s}}) $. (9)式表明颗粒物在靠近出口处压强仅仅与$ D/W $相关. 我们认为$ z = D $处且位于容器中心线上的点也属于该区域, 图4(a)显示了不同$ W $值时$ z = D $处无量纲化压强$ p^{\rm p}(D)/(\rho g L) $随$ D/L $的变化曲线(模拟过程中容器宽度$ L $始终不变). 可见, 当其他物理量不变时, $ D $增大, 出口处压强随之增大, 这与方程$ (9) $的增减性相符合. 且随着$ W $降低, 压强与$ D $的相关性逐渐降低. 当其他物理量不变, $ W $变化时, 颗粒物压强与$ W $强相关. 显见当$ W $降低, 颗粒物在$ z = D $处的压强降低, 同样与方程$ (9) $的增减性相符合. 图4(b)显示无量纲化压强$ p^{\rm p}(D)/(\rho g D) $随$ D/W $的变化曲线, 模拟数据几乎重叠并落到一条曲线上, 图中黑色虚线是利用最小二乘法将模拟结果与方程(9)拟合得到的, 可见模拟结果与方程(9) 符合良好. 可得到如下关系:

$ \frac{p^{\rm p}(D) }{\rho g D} = \frac{ 1.22 }{1 + 0.49 D/W}. $

由方程(10)中$ \gamma_2 = 0.49 $可推出, 当$ D/W\ll 2 $时, 可得到$ p^{\rm p}(D) \propto \rho g D $, 压强仅与$ D $相关; 而当$ D/W\gg 2 $时, 可得到$ p^p(D) \propto \rho g W $, 压强只与$ W $相关. 因此根据$ D/W $的大小, 压强呈现出两种截然不同的相关性. 可见出口在底部情况下出口附近的压强与Zhou 等^[18]得到的出口在侧面的颗粒物流量有相似的规律.
图 4 在$ W $不同的情况下距离出口$ D $处的颗粒物压强$ p^{\rm p}(D) $结果　(a) 无量纲化压强$ p^{\rm p}(D)/(\rho g L) $随无量纲化出口尺寸$ D/L $的变化; (b) 无量纲化压强$ p^{\rm p}(D)/(\rho g D) $随无量纲化出口尺寸$ D/W $的变化
Figure4. Results of granular pressure at the distance $ D $ from the orifice $ p^{\rm p}(D) $ for various $ W $: (a) Dimensionless pressure $ p^{\rm p}(D)/(\rho g L) $ vs dimensionless orifice size $ D/L $; (b) dimensionless pressure $ p^{\rm p}(D)/(\rho g D) $ vs dimensionless orifice size $ D/W $.

3 -->

3.1.2.出口速度结果分析

大量的实验表明^[22,23], 对于出口在底部的情况, 流量始终满足Beverloo方程, $ Q\propto S \sqrt {gD} $, 可见出口处速度没有像压强一样与$ W $相关. 文献[22,23]指出, 颗粒物流量与出口中心处法向速度呈正相关, 因此图5(a)给出了无量纲出口中心处法向速度(即竖直方向速度)$ {v^p_0}^2/(2 gL) $随着$ D/L $的变化曲线. 由图5(a)可见, 当$ D/L $值较大时, 速度随着$ W $变化会发生轻微改变. 图5(b)画出了无量纲化出口中心处竖直速度$ {v^{\rm p}_0}^2/(2gW) $随无量纲化出口尺寸$ D/W $的变化, 图中黑色虚线为通过最小二乘法得到的方程, 可推出出口中心处法向速度与容器厚度及出口尺寸满足如下关系:
图 5 在$ W $不同的情况下位于出口处中心线的竖直方向速度$ v^{\rm p}_0 $结果　(a) 无量纲化速度$ {v^{\rm p}_0}^2/(2gL) $随无量纲化出口尺寸$ D/L $的变化; (b) 无量纲化速度$ {v^{\rm p}_0}^2/(2gW) $随无量纲化出口尺寸$ D/W $的变化
Figure5. Results of vertical velocity on the central streamline $ v^{\rm p}_0 $ for various $ W $: (a) Dimensionless vertical velocity $ {v^{\rm p}_0}^2/(2gL) $ vs dimensionless orifice size $ D/L $; (b) dimensionless vertical velocity $ {v^{\rm p}_0}^2/(2gW) $ vs dimensionless orifice size $ D/W $.

$ \frac{{v^{\rm p}_0}^2}{gW} = \frac{2.24D/W}{1+0.02D/W}, $

由分母上D/W前的系数为0.02可知, 当$ D/W\gg 50 $时, 可得到$ v_0\propto\sqrt{gW} $, 此时流量只与$ W $相关; 而当$ D/W\ll 50 $时, 可得$ v_0\propto\sqrt{gD} $, 此时流量只与$ D $相关, 满足Beverloo定律. 实际情况下, $ D/W\gg 50 $很难被满足, 因此, 在文献[22,23]中, 出口在底部的矩形筒仓卸载问题仍旧满足Beverloo定理, 没有出现流量与$ W $相关的情况.
2 -->

3.2.出口在侧面

图6显示了出口在侧面的情况下流量恒定期间某一时刻压强、速度以及无量纲常数$ I $的分布图. 与出口在底部的情况非常相似, 压强在出口处降低而速度显著增大, $ I $的值在出口处最大, 且最大值为$ 0.1 $左右. 接下来分析不同区域内的压强$ p^{\rm p} $和速度$ v^{\rm p} $值与容器厚度及出口尺寸的关系.
图 6 $ D = 0.40125L $以及$ W = 0.25L $情况下$ t/\sqrt{L/g} = 4 $时刻连续数值模拟结果, 从左至右: 容器内压强与其最大值之比; 总速度与其最大值之比, 其中黑色实线表示颗粒物流线; 无量纲常数$ I = d\sqrt{2}D_2/(\sqrt{p/\rho}) $
Figure6. Continuum simulation results with $ D = 0.40125L $ and $ W = 0.25L $ at time $ t/\sqrt{L/g} = 4 $, from the left to the right: pressure $ p^{\rm p} $ normalized by it's maximum value within the silo; the total velocity $ U^p $ normalized by it's maximum value within the silo, the black lines represent the streamlines; dimensionless number $ I = d\sqrt{2}D_2/(\sqrt{p/\rho}) $.

3 -->

3.2.1.压强结果分析

如图7所示, 研究沿着流线方向距离出口约为$ D $的区域, 可假设颗粒物几乎沿着竖直方向运动, 且流动区域沿$ x $方向的尺寸约为$ D $, 同出口在底部情况做相同的受力分析, 可以得到沿中心流线方向距离出口$ D $处的压强与$ D $和$ W $的关系为
图 7 出口在侧面情况下容器内不同区域受力图, 图中阴影部分代表颗粒物停滞区域
Figure7. Force diagram of different zones within the silo for the case with lateral orifice, the dashed area represents the stagnant zone.

$ \frac{p^{\rm p}(D) }{\rho g D} = \frac{ \gamma_1 }{1 + \gamma_2 D/W}, $

其中$ \gamma_1 = {1}/({2 \alpha_1 \mu_{\rm s}}) $, $ \gamma_2 = ({\alpha_2 \mu_{\rm w}})/({\alpha_1 \mu_{\rm s}}) $. 图8(a)给出了在$ W $不同的情况下距离出口$ D $处的颗粒物无量纲化压强$ p^{\rm p}(D)/(\rho g L) $随无量纲化出口尺寸$ D/L $的变化, 与出口在底部时一样, 当$ W $降低, 颗粒物在$ z = D $处的压强显著降低. 随着$ W $降低, 压强与$ D $的相关性则逐渐减弱. 图8(b)显示了无量纲化压强$ p^{\rm p}(D)/(\rho g D) $随无量纲化出口尺寸$ D/W $的变化, 图中黑色虚线是利用最小二乘法将模拟数据与方程(12)拟合得到的, 所得系数与出口在底部的系数非常接近(见图4(b)), 可见, 出口处的压强大小对于出口的位置并不敏感.
图 8 在$ W $不同的情况下距离出口$ D $处的颗粒物压强$ p^{\rm p}(D) $结果　(a) 无量纲化压强$ p^{\rm p}(D)/(\rho g L) $随无量纲化出口尺寸$ D/L $的变化; (b) 无量纲化压强$ p^{\rm p}(D)/(\rho g D) $随无量纲化出口尺寸$ D/W $的变化
Figure8. Results of granular pressure at the distance $ D $ from the orifice $ p^{\rm p}(D) $ for various $ W $: (a) Dimensionless pressure $ p^{\rm p}(D)/(\rho g L) $ vs. dimensionless orifice size $ D/L $; (b) dimensionless pressure $ p^{\rm p}(D)/(\rho g D) $ vs. dimensionless orifice size $ D/W $.

3 -->

3.2.2.出口速度结果分析

对于出口在侧面的情况, 颗粒物流动方向会发生偏离, 将出口处运动方向与竖直方向的夹角记为$ \theta $. 图9(a)为无量纲化总速度$ {U_0^p}^2/(2gL) $随无量纲化出口尺寸$ D/L $的变化, 可以看出, 出口处总速度大小$ U_0 = \sqrt{u_0^2+v_0^2} $几乎与$ W $无关($ u_0 $表示出口处水平方向的速度, $ v_0 $表示竖直方向的速度), 且$ U_0 \propto \sqrt{gD} $. 水平速度可被表示为一个与总速度$ U_0 $以及出口处运动倾斜角相关的函数($\theta =$$ {\rm {arctan}}(u_0/v_0) $):
图 9 在$ W $不同的情况下位于出口处中心流线上速度结果　(a) 无量纲化总速度$ {U_0^{\rm p}}^2/(2gL) $随无量纲化出口尺寸$ D/L $的变化; (b) 速度倾斜角$ {\rm {sin}} \theta $随无量纲化出口尺寸$ D/W $的变化
Figure9. Results of velocity on the central streamline at the orifice for various $ W $: (a) Dimensionless total velocity $ {U_0^{\rm p}}^2/(2gL) $ vs dimensionless orifice size $ D/L $; (b) angle of inclination $ {\rm {sin}} \theta $ vs dimensionless orifice size $ D/W $.

$ u_0 = U_0 \sin(\theta). $

图9(b)显示了速度倾斜角$ {\rm {sin}}\theta $随无量纲化出口尺寸$ D/W $的变化, 图中黑色虚线为将模拟数据拟合如下方程所得到, 可见$ \sin(\theta) $与$ W $呈现如下相关性:

$ \sin (\theta) = \sqrt{\frac{0.56}{1+0.25 D/W}}\;. $

上述结果表明, 出口处总速度大小只与$ D $有关, 但出口速度倾斜角$ \theta $与$ W $相关, 使得出口中心处法向速度(即水平速度$ u_0 $)受到$ W $影响, 出现了新的流量规律. 该现象可从N-S方程中得到解释, 当$ W $减小, 摩擦力增大时, 速度方向会与重力加速度方向更趋于一致^[18]. 而由方程(14)中分母上$ D/W $前的系数可知, 当$D/W \gg$4时, 可得到$ \sin (\theta) \propto W/D $, 则出口中心处法向速度满足$ u_0 \propto \sqrt{gW} $; 当$ D/W \ll 4 $时, 可得到$ \sin (\theta) $为常数, 则出口中心处法向速度满足$ u_0 \propto \sqrt{gD} $. 可以推测, 对于出口在底部的情况, 出口处速度方向基本只与$ D $相关, 除非$ D/W\gg 50 $. 由此可见, 侧面开口情况下, 出口中心处法向速度与$ W $和$ D $的两种相关性阶段的临界$ D/W $值要远远小于出口在底部的情况.

本文主要通过基于局部本构理论的连续数值模拟方法研究出口在底部和侧面的筒仓卸载问题. 深入研究了容器厚度$ W $(影响容器前后壁摩擦力大小)和出口尺寸$ D $对颗粒物内部物理量的影响, 着重分析了靠近出口处的压强及速度. 首先通过受力分析预测了不同位置的压强与$ W $以及$ D $的关系, 并利用模拟方法进行了验证. 理论分析和数值模拟结果都表明, 距离出口较近区域的压强存在两种相关性, 当$ D/W $足够小时, 压强只与$ D $相关, 当$ D/W $足够大时, 压强只与$ W $相关. 出口在底部和侧面两种情况均呈现上述规律. 随后, 分析了出口处法向速度与$ W $以及$ D $的相关性, 有趣的是, 出口中心处法向速度呈现出与压强不同的相关性. 对于出口在底部的情况, 对于模拟中任意$ D/W $值, 出口中心处法向速度仍然只与$ D $相关; 而出口在侧面时, 当$ D/W $足够小时, 出口中心处法向速度只与$ D $相关, 当$ D/W $足够大时, 出口中心处法向速度只与$ W $相关. 这个结果与一直以来人们所认为的压强是影响出口速度的主要因素明显相违背^[10], 而与近来Perge等^[12]以及Aguirre等^[13]所做实验得出的结论一致, 即筒仓出口处颗粒物流量与压强无关. 对于出口在侧面的情况, 总速度大小只与$ D $相关, 而出口中心处法向速度(水平速度)与$ W $相关, 其原因是出口处速度倾斜角与$ W $相关, 当$ W $降低, 出口处流线会更靠近竖直方向而总速度大小不变, 因此出口中心处法向速度(即水平方向速度)降低. 对于出口在底部的矩形筒仓卸载问题, 本文还揭示了当容器厚度$ W $与出口高度$ D $相比很小时, 流量仍然满足Beverloo定律的原因^[22,23]. 这是由于与出口在侧面相比, 出口在底部时, 造成相关性规律改变的$ D/W $临界值较大, 只有当$ D/W \gg 50 $时, 出口流速(竖直方向速度)才会受$ W $影响, 呈现不同的变化规律, 一般实际情况无法满足该条件.
综上所述, 本文利用基于局部本构理论的连续数值模拟得到了与文献中实验一致的结果^[18,22,23], 并深入剖析了颗粒物内部运动机制, 可见该模拟方法能够较好地预测颗粒流运动规律. 从前文分析中可看出, 停滞区的角度对流量影响较大, 而Kamrin和Koval^[24]所提出的非局部理论能够更好地描述该区域的运动情况, 未来可将局部本构模型修正为更为精确的非局部本构模型, 并比较二者在筒仓卸载问题中结果的区别.