段凡, 贺煊, 周强丨中等雷诺数双分散悬浮系统流固曳力以及固固曳力

文摘 2024-04-15 18:07 北京

点击上方蓝字关注我们

中等雷诺数双分散悬浮系统流固曳力以及固固曳力

段凡 ¹ 贺煊 ¹周强 ^1,²

1. 西安交通大学化学工程与技术学院，陕西西安 710049
2. 西安交通大学动力工程多相流国家重点实验室，陕西西安 710049

DOI：10.12034/j.issn.1009-606X.223212

摘要采用一套全解析数值方法模拟了颗粒可自由移动的双分散悬浮系统，并对文献中已有的流固以及固固曳力公式的准确性进行了检验。模拟的参数范围为整体固含率0.1, 0.2, 0.3，粒径比1.5和2，小颗粒固含率占比0.1, 0.3, 0.5，颗粒-流体密度比10, 100, 500, 1000，整体颗粒雷诺数10, 20, 50。结果显示，对于流固曳力的预测，文献中已有的三类公式中静态均匀系统公式最准，动态悬浮系统公式次之，单分散扩展公式最差。进一步分析发现双分散悬浮系统流固曳力受局部固含率、颗粒相间滑移速度、颗粒拟温度、颗粒Stokes数以及颗粒微结构的影响。在这些因素的共同作用下，流固曳力随颗粒-流体密度比变化不明显，小颗粒相与大颗粒相的流固曳力差异小于静态均匀系统。对于固固曳力的预测，当颗粒-流体密度比等于10或100时，润滑力的作用会使碰撞次数在不同颗粒对之间分布不均，导致分子混乱假设不成立，基于颗粒动力学理论的公式远远高估固固曳力。

关键词流固曳力；固固曳力；双分散悬浮系统；直接数值模拟；颗粒动力学理论；分子混乱假设

1 前言

流体-颗粒两相流动无处不在，自然界中典型的例子有沙尘暴、火山喷发和河流泥沙沉积，工业过程中有颗粒物流化，如煤炭燃烧、生物质气化以及石油催化裂化裂解。在这些流动中，颗粒物尺寸往往不是单一的，而是具有一定分布，甚至存在截然不同的两种尺寸，如生物质和用来辅助生物质流化的沙子。不同尺寸的颗粒终端速度不同，产生不同颗粒相分离的现象，最终对流化床的混合或反应造成显著影响。实验和数值模拟是研究流体-颗粒两相流动现象的两种主要手段。相比于实验，数值模拟的优势在于一方面成本较低，有望应用于工业尺度，另一方面，能够获得更为丰富的信息，允许对流动现象进行更为深入的分析以及更为精密的调控，数值模拟最大的问题在于其结果的准确性。常见的用来模拟流体-颗粒两相流动的数值方法有解析颗粒到表面的直接数值模拟(Particle-resolved Direct Numerical Simulation, PR-DNS)，计算流体力学-离散单元法(Computational Fluid Dynamics-Discrete Element Method, CFD-DEM)和多流体方法(Multi-fluid Method, MFM)。这三种方法对流动的解析度依次降低，计算量相应地依次减小，随着近年来计算机技术的迅猛发展，MFM已经能够直接服务于工业应用^[1-3]，而PR-DNS和CFD-DEM仍局限于基础研究。PR-DNS理论上是完全准确的，因为其从最基本的物理方程出发，未引入任何假设或经验性本构关系^[4]，而CFD-DEM和MFM的准确性主要依赖其所采用的本构关系。

分离现象的发生表明不同颗粒相速度不同，而不同颗粒相之间的滑移速度大小取决于流固曳力、固固曳力以及外力三者的平衡。某一颗粒相的固固曳力定义为其所受到的来自其他颗粒相的平均碰撞力。对于颗粒非常稀少的系统，固固曳力可以忽略，由流固曳力等于外力可得颗粒相终端速度近似正比于粒径的平方，这意味分离现象极易发生。然而，对于相对稠密的流化床，固固曳力会显著降低不同颗粒相之间的滑移速度，从而抑制不同颗粒相分离。由此可知，准确计算流固曳力和固固曳力对于预测分离现象是否发生尤为重要^[5]。CFD-DEM和MFM的计算网格通常远大于颗粒尺寸，无法解析颗粒周围流体的速度分布，因此需要额外的本构关系来计算流固曳力。另外，MFM将颗粒相看作连续介质，无法解析颗粒碰撞，因此还需要额外的固固曳力本构关系。本构关系的研究可在微尺度和介尺度^[6-8]两个层面上进行，后者基于前者的研究成果进一步考虑由不稳定性导致的非均匀结构所带来的影响，本研究仅考虑微尺度层面。

在微尺度颗粒群流固曳力研究中，研究者们习惯于用Stokes力来对其无量纲化，并默认流固曳力指代无量纲后的流固曳力^[9-16]。Stokes力为低雷诺数下单个颗粒在无限大空间、均匀来流中受到的力。对于多颗粒系统，由于颗粒之间通过流体相互影响，流固曳力一般随固含率增大而增大，包含固含率这一参数的经典曳力本构关系有来自实验的Wen等^[17]、Gidaspow^[18]公式以及基于PR-DNS的Beetstra等^[10]、Tenneti等^[13]、Zhou等^[19,20]的公式。这些公式假设所有颗粒尺寸相同，因此被称作单分散曳力公式。对于存在两种以及两种以上颗粒尺寸的系统(双分散系统以及多分散系统)，某一颗粒相所受的流固曳力不仅与整体固含率相关，还受粒径比和固含率占比(该颗粒相的固含率与整体固含率之比)影响。van der Hoef等^[9]首次对低雷诺数双分散系统进行了直接数值模拟，并构建了包含粒径比和固含率占比的曳力公式，而后，Beetstra等^[10]将其扩展至中等雷诺数，该曳力公式如下：

(1)

其中，F_fp,_i为第i种颗粒相的(无量纲)流固曳力，F_fp(ϕ, Re₃₂)代表单分散系统曳力公式，ϕ为整体固含率，即整个颗粒相的体积占比，Re₃₂为基于索特平均直径的颗粒雷诺数，y_i为第i种颗粒相直径与索特平均直径之比。对于常见的流态化系统，颗粒雷诺数的量级不超过10²，研究者们将量级小于1的颗粒雷诺数称为低颗粒雷诺数，量级在1到10²之间的颗粒雷诺数称为中等颗粒雷诺数^{[5,10,13,20-23]}。在低颗粒雷诺数范围内，流固曳力公式可借助低雷诺数相关理论来构建，而中等雷诺数则无成熟的理论，往往需要通过实验或数值模拟构建经验性流固曳力公式。由于对动态双分散系统进行直接数值模拟需要庞大的计算资源以及格子玻尔兹曼方法(Lattice Boltzmann Method, LBM)的弱可压性，本研究仅模拟10, 20, 50这三种相对较低的整体颗粒雷诺数。在van der Hoef等^[9]和Beetstra等^[10]的模拟中，为了简化问题以及节省计算量，颗粒随机均匀排布且速度恒为0。颗粒速度恒定为0意味着颗粒Stokes数趋于无穷大，颗粒速度脉动和颗粒相间滑移速度均为0。采用相同的做法，Yin等^[11]、Cello等^[12]以及Rong等^[14]提出了各自的双分散流固曳力公式，将这类曳力公式统称为静态均匀系统曳力公式。需要注意的是，Beetstra等^[10]、Yin等^[11]以及Cello等^[12]的整体固含率范围为0.1~0.6，而Rong等^[14]为0.62~0.726。静态均匀系统曳力公式都会计算出一个较小的小颗粒相曳力和较大的大颗粒相曳力，如式(1)中，对于小颗粒相，y_i<1，对于大颗粒相，y_i>1，这主要是因为在颗粒随机均匀排布的系统中，大颗粒相实际感受到的固含率高于小颗粒相^[6,11]。在van der Hoef等^[9]的突破性工作之前，研究者们将单分散曳力公式中的颗粒雷诺数替换为每种颗粒相自己的雷诺数(注意固含率为整体固含率)，从而将其扩展至双分散甚至多分散系统，称这类公式为单分散扩展公式。单分散扩展公式包含了不同颗粒相之间滑移速度对曳力的影响，但未考虑局部固含率的差异。低颗粒雷诺数下，单分散扩展公式预测小颗粒相的曳力与大颗粒相相同，这样的结果显然是错误的。尽管如此，Rotondi等^[24]发现在液固沉降系统中，单分散扩展公式的表现远优于静态均匀系统公式，Duan等^[25]认为这是由静态均匀系统公式未考虑颗粒相间滑移速度导致。Rong等^[14]建议将静态均匀系统公式中Re₃₂替换为Re_32,_i来考虑颗粒相之间的滑移速度，Re_32,_i为基于第i种颗粒相自身相滑移速度以及索特平均直径的颗粒雷诺数，称这类公式为动态悬浮系统公式。然而，动态悬浮系统公式是否优于静态均匀系统公式仍无验证。本研究的第一个目标是基于PR-DNS数据对各种静态均匀系统曳力公式、单分散曳力扩展公式以及动态悬浮系统曳力公式进行评价与改进。

动态悬浮系统流固曳力与静态均匀系统的差别以及导致这种差别的因素，尤其对于双分散动态系统，目前仍不清晰。Rubinstein等^[15]进行了低雷诺数(颗粒拟温度非常小)单分散动态悬浮系统的直接数值模拟，发现动态悬浮系统曳力随Stokes数增大而增大，当Stokes数趋于无穷大时，动态悬浮系统曳力回归到静态均匀系统曳力公式。Huang等^[16]通过给位置固定的颗粒施加一个随机速度脉动证明流固曳力随颗粒拟温度升高而增大。由于Huang等施加的速度脉动不随时间变化，因此模拟的是Stokes无穷大的系统。Koch等^[26]采用Stokes动力学相关方法对单分散悬浮系统进行了模拟，发现当Stokes数较大时，颗粒排布结构几乎等同于硬球流体结构，随着Stokes数降低，润滑力的作用使颗粒排布结构逐渐偏离硬球流体结构甚至产生聚团。需要特别强调的是，介尺度研究中的非均匀结构是由不稳定性引起的，与由润滑力导致的非均匀结构不同。Rubinstein等^[15,27]和Duan等^[25]也观察到由润滑力导致的非均匀结构，并发现这种结构明显降低了流固曳力。Stokes数和颗粒拟温度均与颗粒-流体密度比相关，Stokes与颗粒-流体密度比成正比，而颗粒拟温度与颗粒-流体密度比的2/3次方成反比^[28,29]，当颗粒-流体密度比趋于无穷大，Stokes数趋于无穷大，颗粒拟温度趋于0，颗粒排布结构趋于均匀，动态悬浮系统转变为静态均匀系统。这表明颗粒-流体密度比是研究动态悬浮系统流固曳力与静态均匀系统差别的关键。Tavanashad等^[30]通过直接数值模拟研究了颗粒-流体密度比对单分散动态悬浮系统流固曳力的影响，发现随着颗粒-流体密度比降低，Stokes数降低，颗粒拟温度升高，颗粒排布的非均匀性增强，在这三个因素共同作用下，流固曳力随颗粒-流体密度比变化不明显。对于双分散系统，Duan等^[25]发现低雷诺数下，动态悬浮系统小颗粒相与大颗粒相的曳力差异小于静态均匀系统，但并未给出一个清晰的解释。对于中等雷诺数双分散动态悬浮系统，目前还尚未有相关文献报道。本研究的第二个目标是基于第一个目标的评价结果探究颗粒-流体密度比等因素对双分散动态悬浮系统流固曳力的影响。

固固曳力公式一般基于颗粒动力学理论(Kinetic Theory of Granular Flow, KTGF)推导获得^[31-34]。不同固固曳力公式的区别来源于推导过程中采用的数学方法不同。Syamlal^[31]忽略了颗粒速度脉动，Iddir等^[32]以及Chao等^[33]都对速度分布函数进行了泰勒展开，并部分地忽略了颗粒相之间的滑移速度，Zhao等^[34]未采用任何数学简化手段，因此他们所提出的公式理论上最为准确。Duan等^[35]发现在颗粒相对稠密(>0.01)的系统中，Zhao等的公式高估固固曳力，因为分子混乱假设不成立。分子混乱假设是指两个颗粒碰撞前的速度无任何关联性，然而，在稠密系统中，颗粒平均自由程非常短，碰撞后的颗粒很容易再次发生碰撞，导致碰撞前速度存在一定关联性。Duan等研究的局限于颗粒-流体密度比较大的系统，当颗粒密度与流体较为接近时，润滑力对固固曳力的影响不可忽略。一方面，两个颗粒在靠近或分离时，润滑力使其相对速度降低，导致固固曳力减小，当润滑力作用非常强时，颗粒将无法接触，于是固固曳力为0；另一方面，润滑力会影响颗粒排布结构，从而间接地影响固固曳力。固固曳力随颗粒-流体密度比的变化规律是本研究的第三个目标。

2 模拟方法

本研究采用具有二阶精度的格子玻尔兹曼-浸没边界法(Lattice Boltzmann-Immersion Boundary Method, LB-IBM)求解流场^[36]以及流固耦合，离散单元法追踪颗粒运动，软球模型模拟颗粒碰撞，该方法的准确性已在文献^{[19,20,25,37]}中得到验证。当相互靠近(分离)的两个颗粒的间隙远小于一个颗粒直径时，由于颗粒间的流体受到挤压(拉伸)，这两个颗粒将受到一对大小相等且方向相反的斥力(吸引力)，该力被称为润滑力^[38]。由于PR-DNS通常无法解析如此狭小间隙内的流动，因此为了准确计算润滑力，采用如下理论公式对PR-DNS计算值进行修正：

(2)

其中，F_l,_ηγ代表颗粒η与颗粒γ之间颗粒η受到的润滑力，μ为流体动力学黏度，d_η和d_γ为颗粒η和γ的直径，c_η和c_γ为颗粒η和γ的瞬时速度，k_ηγ为由颗粒η中心指向颗粒γ中心的单位向量，又称作碰撞的法向向量，r_ηγ为颗粒η和γ之间的距离，d_ηγ=(d_η+d_γ)/2。由该式可知，随着颗粒间隙尺寸r_ηγ-d_ηγ趋于0时，润滑力趋于无穷大。实际上，由于连续介质假设失效以及颗粒表面存在一定粗糙度，润滑力并不能趋于无穷大，而是在一定的间隙尺寸达到饱和。式(2)在PR-DNS中的具体施加方法以及准确性验证请参考文献[39-43]，根据Brändle de Motta等^[41]的研究结果，取润滑力起始修正位置为0.0625d_s (d_s为小颗粒直径，在本研究中等于16个LBM网格)，即一个LBM网格，润滑力饱和位置0.0005d_s。润滑力饱和位置是一个可调节的参数，与流体平均自由程和颗粒表面粗糙度有关，假如d_s=200 μm，那么润滑力饱和位置为0.1 μm，其对应常压常温下空气的平均自由程。

计算域为边长6d_l的立方体，d_l代表大颗粒直径，边界条件为全周期。在这样一个微小区域内，未出现任何由不稳定性导致的非均匀结构^[35,44]。颗粒初始随机均匀地分布在计算域中，然后在重力以及压力梯度的作用下，流体与颗粒之间的滑移速度从0开始逐渐增大。模拟时间足够长(约为10到50倍Stokes松弛时间)，以确保有充足的统计稳定状态数据进行时间平均，在进行时间平均时，初始非稳定阶段的数据被剔除。不同整体颗粒雷诺数通过改变重力获得，压力梯度在模拟中可被动态调整以保持总质量通量恒定，即：

(3)

其中，u_i为第i种颗粒的平均速度或第i种颗粒相的速度，u_f为流体相的速度，ϕ_i和ϕ分别为第i种颗粒相的固含率和整体固含率，N_c为颗粒种类数，对于本研究关注的双分散系统，N_c=2。上式中i以及后面公式中任何出现多次的下标均为非哑指标。本研究所模拟的系统代表沉降(constant=0)或流态化(constant≠0)系统中远离壁面的微小区域。模拟参数范围如表1所示，此外，颗粒光滑且完全弹性。在LBM模拟中，各种量的单位为格子单位，格子单位可通过相似性准则转化为物理单位^[45]。对于表1中的参数，假设重力加速度为9.81 m/s²，流体为空气，那么小颗粒直径范围为162~2819 μm (A和B类颗粒)，大颗粒直径为280~4653 μm (C和D类颗粒)。该直径范围是本研究所有算例的小颗粒和大颗粒直径范围，对于任意一个算例，小颗粒和大颗粒直径均为单一值，直径比为1.5或2，如表1所示。

表1 模拟参数设置Table 1 Setting of simulation parameter

Note: When x_s=0.1 or 0.5, only the cases of Re₃₂ = 20 are simulated.

第i种颗粒相的固含率占比定义为第i种颗粒相的固含率与整体固含率的比值，即：

(4)

索特平均直径为所有颗粒直径的三次方之和比二次方之和，因此记作d₃₂，其可写作：

(5)

那么y_i=d_i/d₃₂。第i种颗粒相的雷诺数定义为：

(6)

其中，ρ_f为流体密度。整个颗粒相的雷诺数(简称为整体颗粒雷诺数)基于索特平均直径以及颗粒质量加权平均速度u_p定义，如下所示：

(7)

对于本研究考虑的所有颗粒密度相同的情况，质量加权平均等于体积加权平均，因此，

(8)

在动态悬浮系统曳力公式中，式(7)中u_p被u_i替换，即：

(9)

类似地，颗粒脉动雷诺数以及颗粒脉动Stokes数定义为：

(10)

(11)

其中，ρ_i为颗粒密度，T_i为颗粒拟温度，反映颗粒速度脉动的剧烈程度，定义如下所示：

(12)

其中，c_i代表某个第i种颗粒的速度，代表系综平均。上述各颗粒雷诺数以及Stokes数的物理意义分别为：Re_i反映第i种颗粒附近流体惯性力与黏性力的相对大小，Re₃₂整体地反映颗粒附近流体惯性力与黏性力的相对大小，Re_T,i用来衡量颗粒速度脉动的强弱，而Re_32,_i相比Re_i和Re₃₂并无明确的物理意义，引入Re_32,_i仅仅是为了在静态均匀系统曳力公式的基础上进一步考虑颗粒相间滑移速度的影响，St_T_,_i为颗粒脉动速度松弛时间与流体特征时间之比，反映颗粒惯性的大小。

在本研究中，上述各种变量均先在整个计算域空间平均，然后在统计稳定阶段时间平均。图1展示了ϕ=0.2, d_l/d_s=2, x_s=0.3, ρ_p/ρ_f=1000, Re₃₂=20下，第i种颗粒相雷诺数随时间的变化以及tν/d=0, 10, 80时刻(分别位于初始阶段、过渡阶段、统计稳定阶段)的颗粒排布结构，其中被绿色阴影覆盖的区间即对应统计稳定阶段。在任意一个时刻，u_f为计算域所有流体网格点的平均速度，颗粒的u_i, T_i分别为计算域所有第i种颗粒的平均速度、三个方向速度方差的平均。第i种颗粒相所受到的瞬时流固曳力和固固曳力分别按照如下公式统计。

(13)

(14)

图1 颗粒雷诺数随时间的变化以及tν/d=0, 10, 80时刻(分别对应初始阶段、过渡阶段、统计稳定阶段)的颗粒排布结构(被绿色阴影覆盖的区间即对应用来进行时间平均的统计稳定阶段，ϕ=0.2, d_l/d_s=2, x_s=0.3, ρ_p/ρ_f=1000, Re₃₂=20)Fig.1 Variations of the particle Reynolds number with time and the particle configurations at tν/d=0, 10, 80 (corresponding to the initial stage, transition stage, and statistical steady stage, the interval covered by the green shade corresponds to the statistical steady stage used for time averaging, ϕ=0.2, d_l/d_s=2, x_s=0.3,ρ_p/ρ_f=1000, Re₃₂=20)

式(13)和(14)中点积符号前面的式子分别代表计算域内所有第i种颗粒的平均流固曳力(不包含所施加的压力梯度力)和平均碰撞力，即固固曳力，点积符号后面的式子为流向单位向量，F_c_,ηγ为编号为η的第i种颗粒受到的来自编号为γ的第j种颗粒的碰撞力，该碰撞力根据弹簧-阻尼模型计算，更多信息请参考文献[46]，N_i和N_j分别为第i和j种颗粒相的颗粒数目。瞬时流固以及固固曳力除以Stokes力得到无量纲流固以及固固曳力，即：

(15)

(16)

其中，上标#代表无量纲，通常为了书写简便，省略该符号，分母为基于表观滑移速度(1-ϕ)|u_f-u_i|定义的Stokes力。最后，对这些变量进行时间平均得到将要展示的结果。值得注意的是，由于本研究所模拟的系统在垂直于流向的方向是统计对称的，因此流固和固固曳力的垂直分量(通常称作升力)非常接近0。由公式(15)可知，F为流固曳力与Stokes力之比，反映固含率、颗粒雷诺数等其他因素对流固曳力的影响，其与阻力系数C_D_,i以及曳力系数β_i的关系为：

(17)

(18)

为了确保模拟结果达到网格无关，本研究采用d_s/12, d_s/16和d_s/20三种网格尺寸对算例ϕ=0.3, Re₃₂=100进行模拟，并采用Richardson外推获得无穷密网格的结果。需要注意的是，网格要求随固含率以及颗粒雷诺数减小而降低^[9,25,37]，因此能够满足算例ϕ=0.3, Re₃₂=100的网格尺寸必然能满足表1中的算例。计算结果发现，对于大颗粒相(小颗粒相)流固曳力，网格d_s/12, d_s/16, d_s/20相对于无穷密网格的偏差分别为4.5% (7.0%), 1.8% (4.2%), 0.82% (1.03%)，而固固曳力对网格尺寸不敏感，网格d_s/12, d_s/16和d_s/20之间的差距小于1%。为了节省计算量，本研究所选用的网格尺寸为d_s/16。

3 结果与讨论

3.1　流固曳力

3.1.1 颗粒-流体密度比对流固曳力的影响以及流固曳力公式构建

动态悬浮系统与静态均匀系统的本质区别在于颗粒-流体密度比，因此，在对比不同双分散曳力公式以及构建新公式之前，本研究首先探讨颗粒-流体密度比对流固曳力的影响。图2展示了不同整体颗粒雷诺数、不同整体固含率下流固曳力随颗粒-流体密度比的变化。由图可知，在目前参数范围内，流固曳力受颗粒-流体密度比影响较小，与密度比1000相比，曳力平均偏差最小仅为3.2%，最大为25%。对于单分散悬浮系统，Tavanashad等^[30]的计算结果同样表明流固曳力对颗粒-流体密度比不敏感，这种不敏感是颗粒速度脉动、颗粒Stokes数和颗粒微结构三个因素共同作用的结果。更深入的机理分析将在下节进行，本节重点对比各种双分散公式的准确性。

图2 不同整体颗粒雷诺数、不同整体固含率下流固曳力随颗粒-流体密度比的变化(x_s=0.3, d_l/d_s=2.0)Fig.2 Fluid-particle drag vs. particle-to-fluid density ratio at various overall particle Reynolds numbers and overall solid volume fractions (x_s=0.3, d_l/d_s=2.0)

动态悬浮系统颗粒-流体密度比为有限值，而静态均匀系统为无穷大，随着颗粒-流体密度比趋于无穷大，动态悬浮系统逐渐转变为静态均匀系统，同时又因为流固曳力随颗粒-流体密度比变化较小，所以可以推断出，对于目前的参数范围，静态均匀系统曳力公式能够较为准确地预测动态悬浮系统曳力。表2汇总了本研究将要检验的单分散以及双分散曳力公式，其中Re代表单分散系统的颗粒雷诺数，即Re=d(1-ϕ)|u_f-u_p|/μ。选择Beetstra等^[10]、Yin等^[11]以及Cello等^[12]的双分散曳力公式的原因是他们的参数范围与当前研究最为接近，更能突显出动态悬浮系统曳力与静态均匀系统的差别。图3(a)~3(c)分别展示了ϕ=0.1, 0.2, 0.3下表2中各单分散曳力公式计算值随颗粒雷诺数的变化，图3(d)展示了ϕ=0.2, x_s=0.3, Re₃₂=20下表2中各双分散曳力公式计算值随第i种颗粒直径与索特平均直径之比y_i的变化。在实际应用中，许多研究者常常将Gidaspow的单分散公式^[18]与Beetstra等的双分散公式相结合，即将BVK_B中的BVK_M替换为Gid_M，简称为Gid_M-BVK_B。由图3(a)~3(c)可知，单分散系统曳力随固含率和颗粒雷诺数增大而增大；由图3(d)可知，双分散系统中每一颗粒相曳力随y_i增大而增大。因为y_i=d_i/d₃₂，所以对于小颗粒相y_i<1，对于大颗粒相y_i>1，这表明小颗粒相曳力小于大颗粒相，并且随着粒径差异增大小颗粒相与大颗粒相曳力差异增大。又因为y_i=x_i+d_i/d_jx_j=(1-d_i/d_j)x_i+d_i/d_j，假如第i种颗粒相为小颗粒相，那么1-d_i/d_j>0，所以d_i/d_j固定情况下y_i随x_i增大而增大，这表明小颗粒相曳力随其固含率占比增大而增大；假如第i种颗粒相为大颗粒相，那么1-d_i/d_j<0，所以d_i/d_j固定情况下y_i随x_i增大而减小，即y_i随1-x_i=x_j增大而增大，这表明大颗粒相曳力也随小颗粒相固含率占比增大而增大，总之，小颗粒相曳力和大颗粒相曳力均随小颗粒相占比增大而增大。

表2 各种单分散以及双分散流固曳力公式Table 2 Various monodisperse and bidisperse fluid-particle drag relations

Note: Suffix _M and _B respectively represent monodisperse and bidisperse fluid-particle drag relations.

图3 ϕ=0.1 (a), 0.2 (b), 0.3 (c)下各单分散曳力公式计算值随颗粒雷诺数的变化；(d) ϕ=0.2, x_s=0.3, Re₃₂=20下各双分散曳力公式计算值随y_i的变化(由于y_i>x_i，公式CDD_B的横坐标范围为0.3到10)Fig.3 Predicted value of various monodisperse drag relation vs. the particle Reynolds number at ϕ=0.1 (a), 0.2 (b), 0.3 (c); (d) the predicted value of various bidisperse drag relation vs. y_i with ϕ=0.2, x_s=0.3, Re₃₂=20 (the range of the bidisperse drag relation of CDD_B is from 0.3 to 10 due to y_i>x_i)

图4将各种静态均匀系统曳力公式的预测值与PR-DNS结果进行了对比，横坐标为PR-DNS结果，纵坐标为预测值，分别用上标SIM和PRED表示，点越靠近对角线意味着公式预测值越准确。图中，紧跟在公式缩写名后面的是计算时所代入的颗粒雷诺数，颗粒雷诺数的定义见式(6), (7), (9)。子图是小颗粒相曳力与大颗粒相曳力比值(F_fp,s/F_fp,l)的对比。从图4可知，对于小颗粒相以及大颗粒相曳力，BVK_B, YS_B和CDD_B这三个公式整体表现良好，Gid_M-BVK_B公式预测值整体偏低，对于小颗粒相曳力与大颗粒相曳力之比，所有公式预测值均小于模拟值，即高估小颗粒相与大颗粒相之间曳力差异。由于公式BVK_B预测曳力比为粒径之比，因此图4(a)和4(d)中所有点排成两条水平线，反过来由图4(c)可知，公式CDD_B右端第二项远小于第一项。公式BVK_B, YS_B, CDD_B和Gid_M-BVK_B平均(最大)相对偏差分别为9.0% (37%), 7.6% (27%), 9.5% (34%)和19% (49%)。综合曳力以及曳力比看，YS_B公式最为准确，BVK_B和CDD_B公式次之，Gid_M-BVK_B公式最差。

图4 不同静态均匀系统曳力公式预测值与PR-DNS结果的对比(子图为小颗粒相曳力与大颗粒相曳力之比，R²代表拟合优度)Fig.4 Predictions by various drag relations of static homogeneous systems vs. PR-DNS results (the embedded subgraph shows the ratio of small-particle-phase drag to large-particle-phase drag, R² represents the goodness of fit)

图5对比了各种单分散曳力扩展公式的预测值与PR-DNS结果。对比图4和5可发现单分散曳力扩展公式表现明显不如静态均匀系统曳力公式。与图4所示情况不同的是，单分散曳力扩展公式倾向高估小颗粒相与大颗粒相曳力之比，即低估小颗粒相与大颗粒相的曳力差异。BVK_M, CDD_M, Gid_M和TGS_M公式平均(最大)相对偏差分别为10% (40%), 11% (40%), 17% (41%)和13% (39%)。虽然静态均匀系统曳力公式和单分散曳力扩展公式都会计算出一个较小的小颗粒相曳力和一个较大的大颗粒相曳力，但两者所考虑的机制不同，前者的曳力差异来自不同的局部固含率以及颗粒雷诺数，而后者则仅由颗粒雷诺数不同造成。需要注意的是，静态均匀系统曳力公式所考虑的不同颗粒相之间的雷诺数差异仅来源于颗粒粒径，而在动态悬浮系统中，不同颗粒相速度不同也会导致颗粒雷诺数不同。

图5 不同单分散曳力扩展公式预测值与PR-DNS结果的对比(子图为小颗粒相曳力与大颗粒相曳力之比，R²代表拟合优度)Fig.5 Predictions by various monodisperse drag extended relations vs. PR-DNS results (the embedded subgraph shows the ratio of small-particle-phase drag to large-particle-phase drag, R² represents the goodness of fit)

为了考虑由颗粒相间滑移速度造成的颗粒雷诺数差异，Rong等^[14]建议将静态均匀系统曳力公式中的质量加权平均颗粒相速度替换为每种颗粒相各自的速度，即将Re₃₂替换为Re₃₂_,i。图6展示了不同动态悬浮系统曳力公式的预测结果，可以发现，由于包含了滑移速度引起的颗粒雷诺数差异，动态悬浮系统曳力公式相比于静态系统均匀曳力公式进一步高估小颗粒相与大颗粒相的曳力差异。这表明存在除了局部固含率以及颗粒雷诺数之外的其他影响流固曳力的因素。BVK_B, YS_B, CDD_B和Gid_M-BVK_B公式平均(最大)相对偏差分别为11% (50%), 10% (39%), 12% (47%)和20% (61%)。

图6 不同动态悬浮系统曳力公式预测值与PR-DNS结果的对比(子图为小颗粒相曳力与大颗粒相曳力之比，R²代表拟合优度)Fig.6 Predictions by various drag relations of dynamic suspensions vs. PR-DNS results (the embedded subgraph shows the ratio of small-particle-phase drag to large-particle-phase drag, R² represents the goodness of fit)

总之，上述三种流固曳力公式的准确性由好到差依次为静态均匀系统曳力公式、动态悬浮系统曳力公式以及单分散曳力扩展公式。接下来分析这三种曳力公式在一些物理极限条件下是否成立。Yin等^[11]提出低雷诺数条件下，双分散曳力公式应满足以下三个条件：(1) 当所有颗粒粒径相同时，双分散曳力公式应能回归到单分散曳力公式，即当y_i=1，F_fp,_i=F_fp；(2) 当整体固含率趋近于0时，每种颗粒相的曳力都应接近1，即当ϕ→0时，F_fp,_i→1；(3)当小颗粒相的固含率占比以及粒径远小于大颗粒相时，所有小颗粒都可被视为点颗粒，因此每个小颗粒受到的曳力都是基于局部滑移速度的Stokes力，小颗粒平均曳力为基于流体平均速度以及自身平均速度的Stokes力，即当y_s→0时，F_fp,s→1/(1-ϕ) [式(15)]。中等雷诺数条件下，物理极限条件(1)依然成立，条件(2) 应修改为“当固含率趋近于0时，每种颗粒相的曳力都应回归到基于自身颗粒雷诺数的单颗粒曳力，也就是当ϕ→0时，F_fp,_i→F_iso(Re_i)”，条件(3)不再成立，这是因为此时曳力与滑移速度不再是线性关系。单分散曳力扩展公式能够同时满足物理极限条件(1)和(2)，而静态均匀系统曳力公式以及动态悬浮系统曳力公式均无法满足条件(2)，这是因为这些公式所使用的颗粒雷诺数(Re₃₂和Re₃₂_,i)无法随着整体固含率趋于0回归到每种颗粒相各自的雷诺数[式(7)和(9)]。本研究将静态均匀系统曳力公式和动态悬浮系统曳力公式中考虑局部固含率的方法(即修正因子y_i)引入单分散扩展公式，从而提出如下双分散曳力公式：

(19)

当ϕ→0时，F_fp,_i=F，而此时F与经典的单颗粒曳力公式基本一致，详情见文献[10]中图1，故上式同时满足物理极限条件(1)和(2)。图7为该公式预测值与PR-DNS结果对比，平均相对偏差、最大相对误差分别为5.2%, 31%。

图7 新公式(19)预测值与PR-DNS结果的对比(子图为小颗粒相曳力与大颗粒相曳力之比，R²代表拟合优度)Fig.7 Predictions of new relation (19) vs. PR-DNS results (the embedded subgraph shows the ratio of small-particle-phase drag to large-particle-phase drag, R² represents the goodness of fit)

3.1.2 动态悬浮系统流固曳力的影响因素分析

在颗粒排布随机均匀的静态系统中，大颗粒相局部固含率高于小颗粒相，因此大颗粒相曳力比小颗粒相大。这一规律在动态悬浮系统中是否依然成立？为此，本研究对动态悬浮系统中各个时刻的颗粒排布结构进行了Voronoi剖分，某一时刻剖分结果如图8所示。某一个颗粒的局部固含率定义为该颗粒体积比上其对应的多面体体积，表示为ϕ_local,_k，而小颗粒相或大颗粒相的局部固含率为：

(20)

其中，下标t代表时间平均。图9展示了局部固含率与整体固含率之比，横坐标代表动态悬浮系统，纵坐标代表静态均匀系统，分别用DYN和STAT表示。可发现，动态悬浮系统中，大颗粒相局部固含率依然高于小颗粒相，但相比于静态均匀系统，小颗粒相与大颗粒相之间的局部固含率差异缩小。由于流固曳力与固含率正相关，因此动态悬浮系统中小颗粒相曳力更接近大颗粒相，从而导致静态均匀系统公式(图4)以及动态悬浮系统公式(图6)高估小颗粒相与大颗粒相的曳力差异。

图8 Voronoi剖分Fig.8 Voronoi tessellation

图9 动态悬浮系统ϕ_local,_i/ϕ与静态均匀系统对比Fig.9 Comparison of ϕ_local,_i/ϕ between dynamic suspensions and static homogeneous systems

除了上述局部固含率这一因素，流固曳力还受颗粒雷诺数、颗粒速度脉动、颗粒Stokes数以及颗粒微结构的影响。对于静态均匀系统，不同颗粒相的雷诺数差别仅在于颗粒粒径，颗粒无速度脉动，Stokes数为无穷大，颗粒随机均匀排布；而对于动态悬浮系统，不同颗粒相之间往往存在滑移速度，从而导致颗粒雷诺数差异增大，颗粒存在一定速度脉动，颗粒Stokes数为有限值，颗粒排布可能存在一定微结构。图10(a)~10(d)分别展示了不同整体颗粒雷诺数下小颗粒相和大颗粒相的雷诺数、脉动雷诺数、小颗粒相颗粒拟温度与大颗粒相颗粒拟温度之比以及脉动Stokes数[定义见式(6), (10), (12), (11)]随颗粒-流体密度比的变化，整体固含率、粒径比和小颗粒体积占比分别为0.2, 2和0.3。由图10(a)可知，随着密度比增大，小颗粒相雷诺数增加，大颗粒相雷诺数减小，小颗粒相与大颗粒相之间的滑移速度减小。由图10(b)和10(c)可知，随密度比增大，颗粒脉动雷诺数(即颗粒拟温度)减小，而且颗粒拟温度之比趋近于8，即粒径比2的3次方，这表明密度比越大，小颗粒相脉动能量和大颗粒相越接近均分。由图10(d)可知，脉动Stokes数随密度比减小而减小，而且在密度比较高时，小颗粒相脉动Stokes数大于大颗粒相，但在密度比较低时，小颗粒相脉动Stokes数小于大颗粒相，这是由于随着密度比降低，小颗粒相脉动能量逐渐低于大颗粒相。颗粒速度脉动和脉动Stokes数对流固曳力的影响机制可通过Balachandar的理论^[47]来解释，对于固含率极低的系统，每一个颗粒所受的曳力回归到单颗粒曳力，即：

(21)

图10 不同整体颗粒雷诺数下颗粒雷诺数(a)、脉动雷诺数(b)、小颗粒相颗粒拟温度与大颗粒相颗粒拟温度之比(c)以及脉动Stokes数(d)随颗粒-流体密度比的变化(ϕ=0.2, d_l/d_s=2, x_s=0.3)Fig.10 Particle Reynolds number (a), fluctuation Reynolds number (b), the ratio of small-particle-phase granular temperature to large-particle-phase granular temperature (c), and fluctuation Stokes number (d) at various overall particle Reynolds numbers (ϕ=0.2, d_l/d_s=2, x_s=0.3)

其中，F_iso代表单颗粒曳力公式，u_f,_k为颗粒k位置处的流体未扰动速度，u_k为颗粒k的速度，Re_k为基于u_f_,k-u_k的颗粒雷诺数。Re_k三个方向分量为Re_kx, Re_ky, Re_kz，x代表主流方向，y和z代表垂直于x的方向。基于该单颗粒曳力公式，平均无量纲曳力为：

(22)

其中，u_p和Re分别为平均颗粒速度以及平均颗粒雷诺数。令φ=F_iso(Re_k)Re_kx/Re，并对φ在(Re, 0, 0)处泰勒展开，得

(23)

其中，下标0代指(Re, 0, 0)。将上式代入式(22)，一阶导数项系综平均后为0，得

(24)

又因为

(25)

(26)

其中上撇代表相对于平均雷诺数的脉动值。将这两个式子代入到式(24)，最终可得

(27)

由该式可以看到，由于单颗粒曳力与颗粒雷诺数之间的非线性关系，平均曳力不仅与平均颗粒雷诺数相关，还与流体速度脉动、颗粒速度脉动以及流体-颗粒速度脉动关联矩有关。将Schiller等^[48]的单颗粒曳力公式(表2中F_iso)代入到φ中得φ的二阶导均大于0，如下式所示。

(28)

(29)

因此平均曳力随流体速度脉动以及颗粒拟温度增大而增大，随流体-颗粒速度脉动关联矩增大而减小，又因为流体-颗粒速度脉动关联矩随脉动Stokes数增大而减小，所以平均曳力随脉动Stokes数增大而增大。Huang等^[16]通过给静态均匀系统施加一个随机颗粒速度脉动证明曳力随颗粒拟温度增大而增大，然而在目前所研究的动态悬浮系统中，随着颗粒-流体密度比降低，颗粒拟温度升高，但曳力无明显变化。这是因为在Huang等的系统中流体-颗粒速度脉动关联矩为0，而在当前系统中，由于脉动Stokes数随密度比减小而减小，因此颗粒拟温度增大的同时流体速度脉动和颗粒速度脉动的关联性也增强，这两种因素的影响相互抵消，最终曳力随密度比变化不明显。另外，动态悬浮系统中小颗粒相颗粒拟温度明显高于大颗粒，而脉动Stokes数两者比较接近，这使得小颗粒相与大颗粒相的曳力差异低于静态均匀系统。

由于润滑力的作用，两个颗粒在相互靠近或相互分离时相对速度减小，该相对速度减小的比例与脉动Stokes数成反比，比例系数与流体分子自由程有关，对于气体，该比例系数约为O(1)^[35,49]。由图10(d)可知，随着颗粒-流体密度比减小，脉动Stokes数量级从10²减至1，两个颗粒的相对速度的减小比例从接近0%增至接近100%。这意味着当密度比很低时，颗粒对将难以靠近或分离，从而导致一定的颗粒微结构^[26]。此外，脉动Stokes数随颗粒雷诺数减小而减小[图10(d)]，因此当颗粒雷诺数很低时会产生颗粒微结构。Rubinstein等^[15,27]以及Duan等^[25]发现低雷诺数下，颗粒排布并非完全随机均匀，从而使流固曳力降低。本研究采用径向分布函数来衡量这种微结构的强弱，第i种颗粒相与第j种颗粒相之间的径向分布函数定义为：

(30)

其中，r为两个颗粒的距离，N_ij(r)代表距离在r到r+dr之间的颗粒对的数目，n_i为颗粒数密度，n_i=N_i/V，V为计算域体积。当i=j时，即一种颗粒相，在统计N_ii(r)时需注意一对颗粒要记两次，如对于编号为1和2的两个颗粒，12和21各计一次。图11(a)展示了整体固含率0.2、粒径比2、小颗粒占比0.3、整体颗粒雷诺数20时径向分布函数随r变化曲线，图11(b)为不同颗粒-流体密度比下径向分布函数在r=d_ij(接触点)处的值，横坐标代表PR-DNS结果，纵坐标为硬球流体理论值，分别用上标SIM和T来表示，该理论值可通过Lebowitz的公式^[50]来计算，即：

(31)

其中，d_ij=(d_i+d_j)/2。由图11(b)可知，PR-DNS结果与理论值在颗粒-流体密度等于1000时吻合良好，但随着密度比降低，模拟值逐渐高于理论值。当颗粒密度比等于100和10时，模拟值发散，即随着dr减小，模拟值趋于无穷大，这说明系统中存在两个或多个颗粒长时间接触的现象。总之，当颗粒-流体密度比很大时，动态悬浮系统颗粒排布结构为硬球流体结构，即颗粒排布结构是均匀的，随着密度比减小，两个或多个颗粒组成的聚团逐渐产生。值得注意的是，在目前参数范围内，这种由润滑力导致的颗粒聚团很难用肉眼直接观察。对于固含率且颗粒雷诺数极低的动态悬浮系统，Batchelor^[51]理论推导出曳力与径向分布函数接触点值的关系式，如下所示：

(32)

(33)

该式证明随着α增大，即由两个或多个颗粒组成的聚团越来越多，F_fp,_i逐渐减小，这与Rubinstein等^[15,27]以及Duan等^[25]的PR-DNS结果一致。

图11 (a)径向分布函数随颗粒间距离变化曲线(ϕ=0.2, d_l/d_s=2, x_s=0.3, Re₃₂=20)；(b) 径向分布函数接触点值的PR-DNS结果与理论值对比(d_l/d_s=2)Fig.11 (a) Variation of radial distribution function with the distance between particles (ϕ=0.2, d_l/d_s=2, x_s=0.3, Re₃₂=20); (b) comparison of the contact value of radial distribution function between PR-DNS results and theoretical value (d_l/d_s=2)

综上所述，随着颗粒-流体密度比降低，颗粒拟温度升高，脉动Stokes数减小，两个或多个颗粒组成的聚团增多。颗粒拟温度升高增大曳力，而脉动Stokes数减小和聚团增多减小曳力，最终曳力随密度比无明显变化，如图2所示。对于小颗粒相与大颗粒相之间的曳力差异，局部固含率差异缩小，以及小颗粒相颗粒拟温度高于大颗粒相颗粒拟温度这两个因素使得曳力差异缩小，而颗粒相间滑移速度这一因素使得曳力差异增大，最终曳力差异减小，因此静态均匀系统公式以及动态悬浮系统公式都高估曳力差异，如图4和6所示。

3.2　固固曳力

当颗粒速度满足Maxwell分布时，Zhao等^[34]精确地推导出固固曳力表达式：

(34)

(35)

其中，f_pp,_ij代表单位体积区域内第i种颗粒相受到来自第j种颗粒相的固固曳力，单个第i种颗粒的平均固固曳力F_pp,_i=f_pp,_ij/n_i，m_i和m_j、e_ij分别为第i和j种颗粒的质量、第i和j种颗粒的恢复系数。无量纲相对速度为：

(36)

其中，g_ij=u_i-u_j。误差函数为：

(37)

图12展示了R(g)g随g的变化，两条虚线分别为R(g)g在g→0和g→+∞时的渐近线，其表达式如下所示

(38)

图12 R(g)g随g的变化(两条虚线分别为R(g)g在g→0和g→+∞时的渐近线)Fig.12 R(g)g vs. g (the two dashed lines are the asymptotic lines of R(g)g when g→0 and g→+∞, respectively)

由图12以及式(38)可知，R(g)g随g增大而增大，当g=0，R(g)g=0，此外，随着g增大，R(g)g由线性增长逐渐转变为平方增长。结合式(34)和(36)可推出，固固曳力随颗粒相间滑移速度以及颗粒速度脉动均呈线性到平方增长，注意式(34)中R(g)g与T_i+T_j相乘。

图13展示了不同整体固含率下大颗粒相所受到的固固曳力随颗粒-流体密度比的变化，实线(左轴)为PR-DNS结果，虚线(右轴)为理论预测值与模拟值之比。由于固固曳力为相互作用力，因此图13仅展示大颗粒相的固固曳力。由图13可知，固固曳力随整体固含率升高而升高，随颗粒-流体密度比降低而降低。这是因为整体固含率越高，碰撞概率越大，而密度比越低，即颗粒惯性越弱，润滑力作用越强，导致碰撞颗粒对的相对速度损失越来越严重。当润滑力作用非常强时，颗粒在接触之前相对速度完全损失，此时碰撞不会发生，于是固固曳力为0。这种润滑力对固固曳力的降低作用往往被研究者囊括在恢复系数之中^[25,52,53]，Sundararajakumar等^[54]理论推导出该恢复系数为：

(39)

其中，δ₀为润滑力开始明显起作用的位置比上润滑力达到饱和且不再增长的位置，在本研究模拟中，δ₀≈125，值得注意的是，恢复系数e对δ₀不敏感[式(39)中δ₀以lnδ₀形式出现]，因此无须知道δ₀精确值。当St_T≤lnδ₀-1.28时，两个颗粒碰撞后法向相对速度为0，如果无切向相对速度，该颗粒对接下来将以一个整体一起运动。图13中，ZW代表Zhao等^[34]的原始公式，该公式未包含润滑力作用，对于目前所研究的弹性颗粒，其中恢复系数e=1，ZW-SK代表包含润滑力作用的Zhao等的公式，其中e根据Sundararajakumar等的公式(39)来计算。当颗粒-流体密度比等于1000和500时，颗粒脉动Stokes数较大，由式(39)计算得到恢复系数接近1，因此公式ZW和ZW-SK预测值近似相同。与此同时，两个公式预测值均高于模拟值，且偏差随整体固含率增大而增大，这是由分子混乱假设失效造成的^[35]。对于整体固含率较高的系统(>0.01)，颗粒运动范围有限，两个已经发生碰撞的颗粒很容易再次发生碰撞，又因为不同颗粒相之间的碰撞降低其相对速度，从而导致不同颗粒相之间的实际碰撞前相对速度小于整体相对速度，因此固固曳力被高估。当密度比等于100和10时，公式ZW预测值远远大于模拟值，最大偏差可达419% (整体雷诺数等于10时可达809%)，而公式ZW-SK预测值明显低于ZW，更接近PR-DNS结果。

图13 不同整体固含率下固固曳力模拟值(左轴，实线)以及固固曳力预测值与模拟值之比(右轴，虚线)随颗粒-流体密度比的变化(d_l/d_s=2, x_s=0.3, Re₃₂=20)Fig.13 Variations of simulated particle-particle drag (left axis, solid line) and the ratio of predicted particle-particle drag to simulated particle-particle drag (right axis, dotted line) with particle-fluid density ratio at different overall solid volume fractions (d_l/d_s=2, x_s=0.3, Re₃₂=20)

即使包含了润滑力的作用，图13中公式ZW-SK的偏差依然随颗粒-流体密度比降低而升高，这同样是由分子混乱假设失效所导致。图14展示了不同颗粒-流体密度比下小颗粒相与大颗粒相之间的碰撞次数分布，总碰撞次数为3000次，横坐标为小颗粒编号，从1到198，纵坐标为大颗粒编号，从1到57，圆圈的大小以及颜色反映相应颗粒对的碰撞次数。当密度比等于10和100时[图14(a)和14(b)]，碰撞次数在不同颗粒对之间分布不均，存在明显“聚团”。该“聚团”意味着一些颗粒对由于润滑力阻碍其分离，因而长时间紧密相邻，反复碰撞[这与3.1.2节中图11(b)相互印证]。然而，这些颗粒对的碰撞对固固曳力几乎没有贡献，因为紧密相邻的两个颗粒的相对速度接近0。换言之，在密度比等于10和100的系统中，总是速度基本相同的颗粒发生碰撞，这表明分子混乱假设不成立，因此固固曳力大大降低。随着密度比增大，润滑力的束缚作用减弱，两个颗粒碰撞后能够轻易分离，再与其他颗粒发生碰撞，从而碰撞次数逐渐在不同颗粒对之间均匀分布[图14(c)和14(d)]。

图14 小颗粒相与大颗粒相之间碰撞次数分布[总碰撞次数为3000次，横坐标为小颗粒编号，从1到198，纵坐标为大颗粒编号，从1到57，(a)~(d) 颗粒-流体密度比分别为10, 100, 500, 1000，d_l/d_s=2, x_s=0.3, Re₃₂=20]Fig.14 Distributions of collision numbers between small and large particle phases [the total number of collisions is 3000 times, the abscissa represents the number of small particles which ranges from 1 to 198, the ordinate represents the number of large particles which ranges from 1 to 57, the particle-fluid densities of (a)~(d) are 10, 100, 500, and 1000, respectively, d_l/d_s=2, x_s=0.3, Re₃₂=20]

总之，Zhao和Wang的理论公式高估固固曳力有两方面原因：一是润滑力导致相对速度减小；二是分子混乱假设失效。高固含率系统中颗粒运动范围受限以及低颗粒-流体密度比系统中润滑力阻碍颗粒对分离都会使分子混乱假设不成立。润滑力所造成的相对速度减小可通过Sundararajakumar和Koch的理论公式(39)来计算，然而目前暂无任何理论能够考虑两个颗粒碰撞前速度的相关性。一种可行的办法是进行不同润滑力饱和距离下(即不同δ₀)的直接数值模拟，我们将其作为未来的工作之一。

4 结论

本研究对颗粒可自由移动的双分散动态悬浮系统进行了直接数值模拟，整体固含率为0.1, 0.2, 0.3，粒径比为1.5, 2，小颗粒占比为0.1, 0.2, 0.3，颗粒-流体比为10, 100, 500, 1000，整体颗粒雷诺数为10, 20, 50。基于PR-DNS数据，本研究对文献中的三类流固曳力公式，即静态均匀系统公式、单分散扩展公式以及动态悬浮系统公式，进行了检验，并提出了新的公式；分析了局部固含率、颗粒相间滑移速度、颗粒拟温度、颗粒脉动Stokes数、颗粒微结构这几个因素对双分散动态悬浮系统流固曳力的影响；将固固曳力模拟值与颗粒动力学理论预测值进行了对比，并分析了两者之间产生偏差的原因，得到如下结论：

(1) 静态均匀系统流固曳力公式表现最好，动态悬浮系统流固曳力公式次之，单分散流固曳力扩展公式最差。

(2) 随着颗粒-流体密度比降低，颗粒拟温度升高，颗粒脉动Stokes数减小，两个或多个颗粒组成的聚团越来越多，其中颗粒拟温度升高增大流固曳力，而颗粒脉动Stokes数减小以及聚团增多降低流固曳力，最终导致流固曳力随颗粒-流体密度比变化不明显。

(3) 双分散动态悬浮系统中，小颗粒相与大颗粒相之间的局部固含率差异比静态均匀系统小，大颗粒相滑移速度比小颗粒相大，而小颗粒相颗粒拟温度比大颗粒相高，其中局部固含率差异缩小以及小颗粒相颗粒拟温度高于大颗粒相这两个因素降低小颗粒相与大颗粒相的流固曳力差异，而大颗粒相滑移速度比小颗粒相大使曳力差异增大，最终动态悬浮系统曳力差异小于静态均匀系统。

(4) 由颗粒动力学理论获得的公式远远高估固固曳力，尤其当颗粒-流体密度比较低时，这是由润滑力所导致的碰撞相对速度减小以及分子混乱假设失效两方面原因造成的。

(5) 当颗粒-流体密度比非常低时，由于润滑力的作用，碰撞次数在不同颗粒对之间分布不均，存在明显“聚团”，从而导致分子混乱假设不成立，固固曳力大幅度降低。

目前参数范围并不能覆盖流化床系统中所有流动情况，如缺少整体固含率大于0.3、粒径比大于2以及整体颗粒雷诺数大于50的算例。需要注意的是，由于颗粒数以及网格解析度的增加，模拟这些算例需要相当庞大的计算量。另外，本研究未给出新的固固曳力公式，这是因为固固曳力和润滑力饱和距离有关，而当前模拟中，该距离被固定为0.0005d_s。未来工作将进行不同润滑力饱和距离的直接数值模拟，并构建包含润滑力饱和距离的固固曳力公式。

略

Fluid-particle and particle-particle drag forces in moderate-Reynolds-number bidisperse suspensions

Fan DUAN ¹ Xuan HE ¹Qiang ZHOU ^1,²

1. School of Chemical Engineering and Technology, Xi'an Jiaotong University, Xi'an, Shaanxi 710049, China
2. State Key Laboratory of Multiphase Flow in Power Engineering, Xi'an Jiaotong University, Xi'an, Shaanxi 710049, China

Abstract: A set of fully resolved numerical methods are employed to simulate the bidisperse suspensions where the particles are free to translate and rotate according to the effects of the surrounding fluid, and the fluid-particle and particle-particle drag relations in the literature are examined. Three overall solid volume fractions of 0.1, 0.2, and 0.3, two diameter ratios of 1.5 and 2, three small-particle-phase fractions of 0.1, 0.3, and 0.5, four particle-to-fluid density ratios of 10, 100, 500, and 1000, and three overall particle Reynolds numbers of 10, 20, and 50 are chosen. Simulation results show that, among the fluid-particle drag relations available in the literature, in terms of the model accuracy, the relations obtained from static homogeneous systems are the best, the next are those of dynamic suspensions, and the monodisperse drag extended relations are the worst. Based on the simulation data, a new fluid-particle drag relation that meets all physical requirements is proposed. Further analysis reveals that the fluid-particle drag of bidisperse suspensions is influenced by five factors, that is the local solid volume fraction, the slip velocity between different particle phases, the granular temperature, the particle Stokes number, and the particle microstructure. Under the action of these factors, the change of the fluid-particle drag is not significant as the particle-fluid density ratio varies, and the difference of the fluid-particle drag between small and large particle phases is smaller than that in static homogeneous systems. For the particle-particle drag, when the particle-fluid density ratio equals 10 or 100, the collision numbers are unevenly distributed between different particle pairs because of the lubrication force. This uneven distribution of the collision numbers leads to the invalid of the assumption of molecular chaos, and for this reason, the particle-particle drag is highly overestimated by the relation derived from the kinetic theory of granular flow.

Keywords: fluid-particle drag；particle-particle drag；bidisperse suspensions；direct numerical simulation；kinetic theory of granular flow；assumption of molecular chaos

引用本文：段凡, 贺煊, 周强. 中等雷诺数双分散悬浮系统流固曳力以及固固曳力. 过程工程学报, 2024, 24(3): 297-314. (Duan F, He X, Zhou Q. Fluid-particle and particle-particle drag forces in moderate-Reynolds-number bidisperse suspensions (in Chinese). Chin. J. Process Eng., 2024, 24(3): 297-314, DOI: 10.12034/j.issn.1009‑606X.223212.)

作者简介：段凡，博士研究生，动力工程及工程热物理专业，E-mail: duanfan@stu.xjtu.edu.cn；

作者简介：周强，博士，教授，研究方向为多相流模拟和实验、颗粒流体介尺度结构、直接碳/煤燃料电池，E-mail: zhou.590@mail.xjtu.edu.cn

基金信息： ‍国家自然科学基金资助项目(编号：21978228)

中图分类号： O359

文章编号：1009-606X(2024)03-0297-18

文献标识码： A

收稿日期：2023-08-05

修回日期：2023-09-27

出版日期：2024-03-28

网刊发布日期：2024-04-03

http://mp.weixin.qq.com/s?__biz=MzI0MzM2MjgyOQ==&mid=2247528767&idx=2&sn=7a7499863d001cbb046c263ad2249f65

过程工程学报

《过程工程学报》(月刊)创刊于1976年，由中国科学院过程工程研究所主办、科学出版社出版。《过程工程学报》以过程工程科学为学科基础，重点刊登材料、化工、生物、能源、冶金、石油、食品、医药、资源及环境保护等领域中涉及过程工程的原创论文。

段凡, 贺煊, 周强丨中等雷诺数双分散悬浮系统流固曳力以及固固曳力

1 前 言