基于两相流理论的稀疏和致密颗粒流统一模型

何 康，施华斌，余锡平

(1.清华大学水利水电工程系，北京 100084；2.澳门大学土木与环境工程系，澳门 999078；3.港澳海洋研究中心，香港 999077；4.南方科技大学海洋科学与工程系，深圳 518055)

颗粒体和流体作为基本的两相介质，广泛存在于河流输沙、海底滑坡等流动问题中[1−4]。在颗粒体和流体混合物的运动过程中，颗粒的体积分数和运动状态均可能发生变化。如在海底滑坡过程中，滑坡体上方由于发生剪切掺混会形成低浓度的浊流区，其底部则是高浓度的致密颗粒流动区域，两者均可能对深水工程和海底管道等设施造成重大破坏[5−8]。因此，研究多流动情形下颗粒体的运动过程具有重要意义。

数十年来，很多学者采用实验和数值手段对稀疏颗粒体的流动过程进行了研究。早在20 世纪，研究者们就通过室内实验对异重流等稀疏颗粒流的结构，及颗粒粒径对流动的影响进行了分析[9−10]，也有学者建立简化模型对颗粒体的运动过程进行预测[11]。随着数值模拟技术的发展，近些年有学者采用直接数值模拟[12]，大涡模拟[13]和雷诺平均方程模拟[14]对稀疏颗粒流问题进行了研究。其中对颗粒体和流体的瞬时质量、动量守恒方程进行平均后可获得两相平均方程，进一步对各应力项、相间作用项进行封闭可建立基于连续介质假设的两相流模型[15]。该方法在计算效率上更具优势且保证了一定的计算精度，更适用于解决实际问题，但以往的模型大多仅考虑了稀疏颗粒体间的碰撞作用，而忽略了致密堆积状态下颗粒体间的相互作用，因而还很难准确刻画近底高浓度颗粒体的流动过程[16]。

近些年来，致密颗粒流问题也吸引了很多研究者的兴趣。与稀疏颗粒流不同的是，致密颗粒体中颗粒体的体积分数达到自然堆积状态，此时颗粒与颗粒间的持续接触作用不可忽略[17]。此外，致密颗粒体和流体的相互作用，如拖曳力、孔隙水压等均对颗粒体的运动过程有着重要影响，需要进行细致刻画[18−19]。通过考虑颗粒体的摩擦和碰撞作用，耦合颗粒与流体的相间作用建立的两相流模型为研究颗粒体的流动过程提供了新的视角[20]。然而，目前还未探究该两相流模型在统一描述稀疏和致密颗粒流问题的可能性。此外，分析不同粒径颗粒体和流体混合物的运动过程也有助于理解实际的颗粒流问题。因此，本文采用统一的两相流模型同时对不同粒径的稀疏和致密颗粒体流动问题开展了研究，并对流体紊动的影响，颗粒体的运动过程及其与流体的相互作用进行了分析，以加深对多流动情形下颗粒体和流体混合物流动问题的理解。

1 数学模型

1.1 基本方程

本文采用文献[20]中的两相流模型，平均后的颗粒体和流体控制方程为：

式中：下标 s 和 f分别为颗粒和流体相；本文仅考虑立面二维流动问题，下标i和j分别为水平和垂直方向坐标，并遵循求和约定；α为体积分数；ρ为相密度；U为速度；p为压强；τ为剪切应力；I为相间作用力；g为重力加速度。

1.2 颗粒相本构关系

将颗粒相应力划分为摩擦与碰撞应力，可表示为：

式中，上标 f 和 c分别为颗粒应力的摩擦与碰撞部分。

摩擦剪切应力可写为：

摩擦区域内的颗粒相压强采用以下公式进行计算[23]：

颗粒相碰撞应力基于颗粒动理论模化，可表示为[25]：

式中：Θs为颗粒温度，可由颗粒温度控制方程求解，详见文献 [20]；η=(1+e)/2，e为颗粒的回弹系数；R=(2−αs)/[2(1−αs)3]为颗粒的径向分布函数[26]；和为碰撞粘性系数和碰撞体积粘性系数，有：

图1 给出了不同区域内颗粒相压强随颗粒体积分数的变化。从图1 可以看到颗粒体积分数较小时，稀疏颗粒体仅承担碰撞压强；当颗粒体积分数增大到可堆积状态时，致密颗粒体将同时承担碰撞和摩擦压强，且颗粒体总压强将逐渐由碰撞压强主导转变为由摩擦压强主导。因此，通过综合考虑颗粒体的摩擦与碰撞作用，本文实现了对不同浓度颗粒流的统一描述。

图1 颗粒相压强随颗粒体积分数的变化Fig.1 Variation of particle pressures with the particle volume fraction

1.3 流体相本构关系

考虑紊流的影响，流体相剪切应力可模化为：

基于k−ε紊流模型，有：

式中：Cμ为常数，一般取为0.09；流体相紊动能kf和紊动能耗散率 εf通过求解考虑拖曳力修正后的k−ε紊流方程得到，详细介绍可参考文献[20]。

1.4 相间作用力

这里主要考虑颗粒与流体的相间拖曳力，表达式为：

式中：Kfs为流体与颗粒的相间拖曳力系数；右边第二项为颗粒由于流体紊动导致的漂移速度项，其中为流体紊动运动粘性系数；为施密特数，根据已有的研究结果取值1.0[27]。

这里采用经典的Gidaspow(1994)拖曳力公式[28]，Kfs可表示为：

式中，Cd经由以下公式确定：

式中，Res=αfρf|Uf−Us|ds/μf为颗粒雷诺数。

1.5 计算方法及边界条件

该两相流模型基于开源平台MFIX(multiphase flow with interphase exchanges)离散并求解。模型的控制方程通过矩形结构化网格，以有限体积法离散，并采用改进后的SIMPLE 算法求解压力相关方程[29]；对流项的离散采用了二阶Superbee 格式，并采用了延迟校正方法处理[30]；模型的计算在满足Courant-Friedrichs-Lewy(CFL)收敛条件下，采用自适应时间步长的方法在时空上推进。

本文以稀疏和致密颗粒流为研究对象进行数值模拟。其中稀疏颗粒流算例以文献[10]中的实验进行设置，如图2(a)所示；致密颗粒流算例采用了文献[31]中的实验数据，其设置如图2(b)所示。实验中，初始时刻稀疏和致密颗粒体均放置在水槽一侧，随着闸门开启，颗粒体在重力和密度差作用下向前发展，并最终形成稳定的沉积体和堆积体。模型应用时将所有固壁设置为无滑移边界，对顶部采用刚盖假定，忽略流体自由表面的波动影响。

图2 稀疏和致密颗粒流数值模拟设置 /mFig.2 Numerical setup for dilute and dense granular flows

2 稀疏颗粒流研究

2.1 模型验证

本文选取文献[10]中不同粒径的稀疏颗粒流实验进行验证。其中细、粗颗粒粒径分别为ds=0.025 mm和ds=0.069 mm ，密 度 均 为ρs=3217 kg/m3，初始颗粒体积分数为=3.49×10−3，实验设置详见图2(a)。在模拟设置中，颗粒摩擦角ϕs=30◦，忽略μ0s的影响，回弹系数e=0.8；实验中流体为水，密度取为 ρf=1000 kg/m3，粘性系数μf=0.001 kg/(m·s)。经网格敏感性分析，模拟选取5 mm×2.5 mm的均匀网格进行计算。

模拟与实测的稀疏细、粗颗粒体的前端运动距离变化结果如图3 所示，其中以 αs/=1等值线的前端位置作为稀疏颗粒体的前端运动距离。从图3 可以看到，该模型可以很好地刻画稀疏细颗粒体与粗颗粒体的运动过程，证明了该模型在稀疏颗粒流问题中的有效性。此外，模拟与实验结果均表明初始运动阶段稀疏细、粗颗粒体的前端运动距离基本一致，但在t≈20 s后，稀疏粗颗粒体的前端运动距离开始明显落后于稀疏细颗粒体的，说明不同粒径稀疏颗粒体的运动过程存在较大的差异，需要开展进一步的分析。

图3 模拟与实验的稀疏颗粒体前端运动距离比较Fig.3 Comparisons of simulated and measured front distances in dilute granular flows

2.2 相间作用的影响

图4 给出了稀疏细、粗颗粒相对浓度 αs/的分布变化。由于颗粒和流体间存在密度和速度差异，可以发现，稀疏颗粒体在与流体的剪切运动过程中前端会呈现出近似椭圆形的头部结构，后方则是悬浮颗粒体的主体结构和狭长的尾部结构，并在其中能发现数个涡旋结构，颗粒体和流体交界面整体呈现出一系列的开尔文－亥姆霍兹不稳定性结构。此外，也能观察到稀疏细颗粒能在较长时间内保持向前运动，并维持明显的头部结构；而t=15 s和t=20 s时保持悬浮状态的粗颗粒浓度则明显更少，其头部结构范围也逐渐缩小，尽管此时细、粗颗粒体的前端运动距离差异较小。

图4 细颗粒和粗颗粒相对浓度分布变化Fig.4 Evolutions of the relative concentration fieldfor fine particles and coarse particles

从相间作用来看，不同粒径颗粒体与流体的运动差异主要由其受到的相间拖曳力决定。这里给出t=5 s时细颗粒和粗颗粒与流体的相对速度Us−Uf和相对拖曳力K(Uf−Us)/ρfg的分布，如图5所示。从图5 可以看到，Us−Uf的方向基本垂直向下，表明颗粒体相对流体主要做垂向沉积运动，且粗颗粒与流体的相对速度要明显大于细颗粒与流体的相对速度，即粗颗粒的沉速更大。从相间拖曳力的分布来看，其方向基本垂直向上，即对颗粒体的沉降起阻碍作用。值得注意的是，壁面处颗粒受到的拖曳力明显更大，这是壁面附近沉积颗粒浓度较大的缘故。此外还能发现，涡旋结构内粗颗粒受到的拖曳力小于细颗粒的，而壁面附近粗颗粒受到的拖曳力作用则明显大于细颗粒的，即粗颗粒受到的沉积作用更强。

图5 t =5 s 时稀疏细、粗颗粒流中流体与颗粒相间速度 U s −Uf 和相对拖曳力 K(Uf −Us)/ρfg 比较(黄色和白色虚线分别为=1 和 =0.1的等值线)Fig.5 Comparisons of interphase velocity Uas −nUdf relative draKg(U ff −oUrs)c/ρefg between the fluid and solid phase in dilute fine and coarse particle flows at t(=th5e s yellow and white dash lines are contour lini =th1 and =0.1)

图6 展示了t=20 s时细、粗颗粒体中流体和颗粒速度的比较结果。从图中可以看到，细颗粒中流体和颗粒体的速度更大，这是因为细颗粒的重力沉降作用偏弱，更容易在运动中保持悬浮状态，因而悬浮颗粒体及周围流体所能维持的运动强度更大。粗颗粒的沉降速度明显大于细颗粒，垂向沉积作用更强，更容易丢失运动能量，因而颗粒体及其周围流体速度要更小。故而在颗粒体的向前发展过程中能发现，t≈20 s后细颗粒体的前端运动距离明显大于粗颗粒体的，且在运动过程中处于悬浮状态的颗粒更多。

图6 t =20 s 时稀疏细、粗颗粒流中流体速度 Uf 与颗粒速度 Us 比较(黄色和白色虚线分别为=1 和 =0.1的等值线)Fig.6 Comparisons of fluisd velocity Uf and solid velocity Us in dilute fine and coarse particle flows at t=20 s(the yellow and white dash lines are contour lines with and =0.1)

2.3 流体紊动的影响

在稀疏颗粒体运动过程中，流体的紊动掺混可能对颗粒体和流体的不稳定结构产生影响。图7给出了稀疏细颗粒流中流体速度 |Uf|和相对紊动粘性系数的分布变化。从t=5 s到t=20 s，可以看到流体速度 |Uf|逐渐减小，其较大值主要分布在颗粒体前端和中间区域。值得注意的是，该阶段流体的紊动粘性系数仍然呈现增大的趋势，的值甚至可达100 以上，且大都分布在颗粒体的后上方，前端和下方的值则相对较小，并在底部形成了长条形的低紊动粘性系数区域。

图7 稀疏细颗粒流中流体速度 | Uf|和 相对紊动粘性系数的分布变化(黄色和白色虚线分别为=1 和 =0.1的等值线)Fig.7 Evolutions of the fluid velocity |aUnfd| the relative turbulent viscin dilute fine particle flows (the yellow and white dash lines are contour lines with and =0.1)

为了进一步分析紊流对流体运动的影响，图8给出了是否考虑紊流模型时细颗粒中流体涡量|ωf|=1/2|∂Uf,j/∂xi−∂Uf,i/∂xj|的比较结果。从图8可以明显地看到，不考虑紊流模型时的流体涡量要明显大于考虑紊流模型时的，即不考虑紊流模型时流体的涡旋运动更为强烈，这是因为此时流体的运动没有受到紊流的影响，因而在运动过程中耗散的能量更少。

图9 给出了不考虑紊流模型时细颗粒和粗颗粒的相对浓度分布。与图3 中考虑紊流模型的计算结果相比，不考虑紊流模型后悬浮颗粒体的前端会呈现明显的“抬起”现象，即颗粒体的头部会被水流抬起，且前端运动距离要更小，与实验结果不符。此外，不考虑紊流模型后悬浮颗粒体的后端结构发展得更为无序，可见紊流对颗粒体和流体交界面处的结构有着重要影响。

图9 不考虑紊流模型时细颗粒和粗颗粒浓度分布演变Fig.9 Evolutions of the concentration field for fine particles and coarse particles without turbulence model

为了进一步分析流体紊动对流体和颗粒体的影响，图10 给出了不考虑紊流t=20 s时细、粗颗粒体前端流体速度Uf和颗粒速度Us。与图5 相比，不考虑紊流模型时涡旋结构处流体和颗粒的速度均要比考虑紊流时的更大，这是因为在运动过程中没有受到流体紊动粘性的抑制作用。从运动方向上来看，不考虑紊流模型时流体和颗粒的运动有明显的上下起伏变化，而考虑紊流后流体和颗粒的流动方向则基本一致向前。数值模拟结果表明，颗粒体与流体在剪切运动过程中会在颗粒体后上方产生较大的紊动粘性系数，其分布结构抑制了流体涡旋结构的发展。在流体紊动耗散作用下，流体的运动方向倾向于向前发展，进而也促进了颗粒体的运动。因而紊流尽管使得流体的涡量和速度减少了，但在运动方向上有利于流体和颗粒体的发展。

图10 不考虑紊流模型 t =20 s 时刻稀疏细、粗颗粒流中流体速度 Uf 和颗粒速度 Us 比较(黄色和白色虚线分别为 =1 和 =0.1的等值线)Fig.10 Comparisons of the fluid velocity Uf and the solid velocity Us in dilute fine and coarse particle flows at t=20 swithout the turbulence model (the yellow and white dash lines are contour lines with and =0.1)

图11 模拟与实验中的颗粒堆积表面比较Fig.11 Comparisons of simulated and measured granular deposit surfaces

3 致密颗粒流研究

3.1 模型验证

本文同时对不同粒径颗粒堆积体的坍塌过程进行了模拟，以对该模型在致密颗粒流问题中的有效性进行验证。其中细、粗颗粒粒径分别为ds=0.84 mm和ds=6.02 mm，颗粒摩擦角分别为ϕs= 2 6°和 ϕs= 3 1°，密度为 ρs=3.6×103kg/m3，初始体积分数为=0.59，实验设置详见图2(b)。在模拟中，细、粗颗粒体中的取值分别为2.5 kg/(m·s)和0.5 kg/(m·s) ，回弹系数均取为e=0.8；实验中流体为水，密度取为 ρf=1000 kg/m3，粘性系数μf=0.001 kg/(m·s)。经网格敏感性分析，模拟选取5 mm×5 mm的均匀网格进行计算。

3.2 拖曳力的影响

为了对细、粗颗粒体的坍塌过程开展进一步分析，图12 给出了不同时刻细、粗颗粒与流体的相对速度Us−Uf分布，其中虚线表示颗粒体的堆积表面。从中可以看到，细颗粒坍塌体中颗粒与流体的相对速度较小，且最大值仅集中于颗粒体表面；粗颗粒坍塌体中颗粒与流体的相对速度则较大，且在颗粒体内部也会存在较大的相对速度。

图12 颗粒体与流体相对速度 U s −Uf 的变化Fig.12 Evolutions of the relative velocity between solid and fluid phases Us −Uf

进一步地，对流体作用在颗粒体上的相对拖曳力K(Uf−Us)/ρfg进行分析，所得结果如图13所示。在t˜=1时刻，流体拖曳力的较大值均分布在颗粒体堆积表面附近，其中细颗粒堆积表面附近的流体拖曳力朝向外侧，有利于颗粒体的向外发展；而粗颗粒体受到的流体拖曳力则朝向左上方，与颗粒体坍塌方向相反；细、粗颗粒体在内部受到的流体拖曳力则均朝向左上方，对颗粒体的运动起阻碍作用。在t˜=3 和t˜=5时刻，细颗粒堆积表面附近的仍存在较大的流体拖曳力，其方向仍然朝向外侧；而粗颗粒堆积表面附近的流体拖曳力则明显更小，其中=3时流体拖曳力对颗粒体前端仍然起阻碍作用；细、粗颗粒体内部的拖曳力方向则均呈现较为复杂的局部变化，既可能促进颗粒体流动也可能阻碍颗粒体的运动。

图13 流体与颗粒体相对拖曳力 K(Uf −Us)/ρfg 的变化Fig.13 Evolutions of the relative drag force between solid and fluid phases K(Uf −Us)/ρfg

3.3 流体动压的影响

如图14 所示，受到颗粒坍塌体的驱动，周边流体也会产生运动。在=1时刻，颗粒体和流体的最大速度均集中于右前方，并且可以看到流体在运动过程中发展出了明显的涡旋结构，即在颗粒体右上方形成了流动中心，该处流体速度较小而周围流体速度较大。在t˜=3时刻，颗粒和流体的流动速度明显变大，其最大值仍集中于右前方，且在颗粒堆积表面上方形成了2 个涡旋区域。在t˜=5时刻，颗粒体的速度明显减小，但颗粒体表面上方的涡旋区域仍维持较大的运动强度，可继续影响颗粒体的堆积形态。

图14 细颗粒体坍塌过程中颗粒速度 Us 和流体速度 Uf的变化Fig.14 Evolutions of the particle velocity Us and the fluid velocity Uf during the collapse process of fine particles

考虑到流体压强的变化可能对颗粒体的运动产生影响，这里去除流体静压的影响，定义流体相 的 动 压 强=pf−ρfg(Hw−y)。图15 显 示 了细、粗颗粒体坍塌过程中流体动压及相对动压梯度力的变化。从图15 可以看到，在t˜=1时流体正压大多集中于颗粒体前端，负压大多集中于颗粒体的上方，且细颗粒中的动压绝对值和动压梯度力明显更大，这是流体受到的颗粒阻力越大，产生的动压越不容易消散的缘故。t˜=3 和t˜=5时颗粒体内部的正压逐渐消散，但在颗粒堆积表面上方则发展出了2 个负压区域，相应的相对动压梯度力大小可达1 以上，即与静水压梯度力相比不可忽略。值得注意的是，t˜=1时坍塌体前端正压区域内动压梯度力促进了颗粒体的向前发展；t˜=3 和t˜=5时发展出的负压区域则主要对颗粒表面形态的塑造起作用，对颗粒堆积体前端的影响较小。

图15 坍塌过程中流体动压 (云图)和动压梯度力 − ∇/ρfg(箭头)的变化Fig.15 Evolutions of the dynamic fluid pressure (cloud charts) and the gradient of dynamic fluid pressure−∇/ρfg(arrows) during the collapse process

4 结论

本文应用综合考虑不同流态下颗粒体本构关系、流体与颗粒体相间作用以及流体紊动影响的统一两相流模型，对稀疏和致密颗粒流问题开展了数值模拟研究，揭示了不同粒径颗粒体与流体的相互作用机制。得到以下结论：

(1) 该两相流模型可以很好地模拟稀疏和致密颗粒体的运动过程，对稀疏颗粒体的前端运动距离和致密颗粒体的堆积表面均能很好地刻画，证明了该模型在稀疏和致密颗粒流问题中的有效性。

(2) 对稀疏颗粒流的数值研究结果表明，不同粒径颗粒体与流体的相间作用有明显的差异。其中，细颗粒体与流体的相间速度明显更小，受到的沉积作用更弱，更容易在流体中保持悬浮状态，因而稀疏细颗粒体及周围流体所能维持的运动强度要更大，颗粒体的前端运动距离也要更长。

(3) 紊流在稀疏颗粒体的运动过程中起着重要作用。模拟结果表明：紊动粘性系数与流体自身粘性系数的比值可达100 以上，考虑紊流后流体的运动明显受到紊动粘性系数的耗散作用，流体涡旋结构的发展受到了抑制，流体则倾向于向前运动，进而促进了稀疏颗粒体的发展。

(4) 流体与颗粒的相间拖曳力在不同粒径致密颗粒体中有着不同的分布。其中细颗粒体堆积表面在较长时间内均有着较大的拖曳力，其方向朝向堆积表面外侧，促进了颗粒体的向外发展；粗颗粒体堆积表面的拖曳力则在初始时刻较大，阻碍了颗粒体的运动，但在后续坍塌过程中逐渐减小。

(5) 流体动压对致密颗粒体的运动起着重要作用，且在细颗粒体中更为明显。颗粒坍塌过程中产生的流体动压梯度力与静水压梯度力相当，其中，正压在初始时刻促进了颗粒体前端的发展，负压则主要影响了颗粒体堆积表面的塑造。