驻点流中的颗粒动力学

2021-06-24陈坚强袁先旭

空气动力学学报 2021年3期

李青，陈坚强，毕林，袁先旭，*

(1. 中国空气动力研究与发展中心空气动力学国家重点实验室，绵阳 621000；2. 中国空气动力研究与发展中心计算空气动力研究所，绵阳 621000)

0 引言

颗粒流体两相流广泛存在于化工及航空航天等领域，如流态化反应器、航空发动机和大推力火箭燃烧室等。在这些实际问题中，对于弥散相，一般颗粒数目NP从 O(1 02) 到 O(1 09)，颗粒相对携带流体的密度比从 O (1) 到 O(103)，如化工管道中的颗粒悬浮物或者是河流泥沙悬浮颗粒如气固流化床、航空航天中的颗粒湍流在某些极端条件下，如高温颗粒壁湍流，密度比可以达到O(1 04)，颗粒雷诺数体的密度，μ 是流体的动力黏度，Uf是流体在颗粒重心位置的欧拉速度，Vp是颗粒的拉格朗日速度，a是其中 ρf是流颗粒半径。从 O(1 0-3) 到 (1 03)，从 O(1 0-3)到 O(1 03)…… 特征参数覆盖的范围相当广泛，问题复杂，物理内容丰富。不同的工业背景下，我们遇到的问题也不同，因此提炼的科学问题以及解决的方法也不同。比如化工领域关心传热传质的宏观量、质量和热量增量，用于指导优化过程工程设计；航空航天领域则关心阻力、压力分布和热流密度分布，用于优化各种飞行器的性能；生物医学领域，比如心脑血管系统则关心颗粒沉积规律或者惯性聚集规律，用来预测血栓或从血液里过滤癌细胞，从而更好预防疾病或者设计相关医疗器械等。一般地，开发颗粒流体欧拉-欧拉两相连续模型是学术界的目标，从而为解决工业界问题提供可靠的技术手段。欧拉模型把弥散相当成连续介质考虑，然后根据具体问题进行各种简化假设，建立封闭关系式，然后用颗粒解析或点颗粒的方法把这些关系式耦合到欧拉框架下（Li[1], Morris[2]）。为了解决这个问题，我们需要理解颗粒两相流在不同几何边界下的流体动力学机制。比如Morris[2]只考虑平行剪切流的情况，驻点流情况就无法适用。驻点流的特点是，远场速度大，然后随着流体质点运动到近壁面，流体质点减速，动能转化为压力。因此，（以流体质点沿着对称轴法向运动为例），沿着法向的靠近壁面，流体质点速度不断减小，压力不断增大，当流体质点运动到驻点的时候将满足壁面无滑移和无穿透条件，滞止在驻点，此时驻点压力最大，流体质点动能为零；另一方面，驻点的巨大压力推动流体质点远离驻点，压力不断转化为流体质点动能，随着流体质点不断远离驻点，流体质点速度不断增加，压力减小。通过分析得知，驻点流的速度场和压力场是全场非均匀的。

在工程应用中，如页岩气开采，通过把带颗粒的非牛顿流体垂直向下注入地表层中，液体压力会使得岩石破裂，而悬浮颗粒会占据岩石裂缝的空间，使得岩石无法闭合。这样，水压力就可以进一步破裂岩石，从而为开采地下的天然气提供条件。因此近年来基于工程需求的推动，驻点流中的颗粒动力学蓬勃发展。例如 Vigolo[3]用点颗粒模型耦合直接数值模拟和物理实验研究T管流中的大密度比颗粒群运动（图1）；Li[4-5]用颗粒解析的直接数值模拟PR-DNS研究一个中性颗粒在三维轴对称驻点流对称轴上的运动；Manoorkar[6-7]通过物理实验和直接数值模拟，研究T管流中有限尺寸中性颗粒群的输运特征；Rallabandi[8]用理论解研究在低雷诺数条件下，一个颗粒在驻点流中的运动。Vigolo[3]发现，颗粒反弹系数和颗粒惯性高度相关，而且由于三维涡效应，颗粒反弹系数分布是弥散状的。Manoorkar[6-7]也同样发现T管流的三维涡效应，还发现颗粒输运惯性分选的特征。Wang[9]通过物理实验研究了驻点流的三维涡效应，他们通过在远离壁面的距离施加不同的扰动，得到了一个类似于卡门涡的脱离涡对，这个涡对在靠近壁面时始终不撞击对称中心（图2）。工业问题中的驻点流非常容易出现三维涡效应，抑或是流场本身出现三维涡结构（Vigolo[3]），抑或是来流扰动导致（Wang[9]）。鉴于此，Li[4-5]采用三维轴对称的数值控制域，考虑中性颗粒只在对称轴上运动的情况，移除了驻点流三维不对称效应带来的复杂性，得出了相似结论—颗粒惯性主导驻点流中的颗粒动力学，并且得出了更多的进一步细节。

图1 T型管流运动的示意图Fig. 1 Particles clustering in a T junction

除了Vigolo[3]的大密度比颗粒群课题，大多数研究的共同点都是：颗粒惯性与流体惯性相当St～O(1)，密度比背景流动是全场非均匀驻点流。这种情况下，一方面，颗粒惯性和环境流体惯性相当；另一方面，颗粒的有限尺度往往和流体微团最小水动力尺度相当，展现出高度的局部各向异性特征，换句话说，颗粒相与流体进行动量交换的方式是空间非均匀的[10]（图3）。因此，常规的点颗粒模型数值方法将不再适用，必须采用考虑有限尺寸效应的颗粒解析的数值算法：如有限元、格子玻尔兹曼法LBM[11]，贴体网格下的有限差分或有限体积[12]，虚拟点法[13]，浸没边界方法（Immersed Boundary Method，IBM）[14]。考虑颗粒的有限体积效应一直是颗粒流体两相流的挑战之一，尤其是在颗粒湍流中，有限体积的颗粒注入各向异性的动量扰动，会影响湍流的脉动统计量[15]。因此，为了保证数值实验的正确性，Li[4-5]采用了浸没式边界方法IBM耦合有限体积流体求解器的数值算法对驻点流的颗粒动力学进行分析。

图2 驻点流来流扰动涡对在近壁面的非对称特性[9]Fig. 2 The non-axisymmetric characteristic of a perturbed vortex pair in a stagnation point flow[9]

图3 有限尺寸悬浮颗粒尺度与流体微团最小运动尺度量级相当情况下的示意图Fig. 3 A schematic of a suspended particle with the same order of characteristic length scale as that of the flow

接下来，简要介绍和回顾Li[4-5]的核心工作。Li特别关注驻点流中的颗粒动力学机制问题，这个领域还罕有研究，但是却很重要。在实际工业问题中，颗粒悬浮的装置往往是复杂几何外形，包括驻点流情况。开发颗粒两相流的欧拉-欧拉连续模型则首先需要对复杂几何外形下的颗粒悬浮机理有清晰的认识。驻点流颗粒悬浮问题是一类特殊和典型的问题，比如常见的化工T型管道。

1 颗粒解析的直接数值模拟方法：浸没式边界方法IBM

方程(1)和方程(2)概括描述了颗粒解析的直接数值模拟方法IBM。式中 Ω 是颗粒体积，u是流场速度矢量，Vp是颗粒速度矢量，是满足流体与刚性颗粒边界无滑移条件的动量交换矢量,α是光滑函数以确保数值稳定性，是在颗粒表面拉格朗日点上流体相和颗粒相的滑移速度矢量，dt是流体求解器的时间步长[14]。一般来说，IBM方法可以计算颗粒滑移雷诺数Rep～O(100)的问题，IBM方法也可以用在复杂外形计算湍流问题[16]。

在湍流问题中，随着Re增加，对数值算法的网格解析精度要求增加，因为多尺度的湍流的最小脉动尺度K41随之减小[17]。值得注意的是，流场雷诺数和颗粒滑移雷诺数是两个不同的概念。当流场雷诺数Re很高的时候，颗粒滑移雷诺数Rep可以很小。一般来说，工业问题中的高雷诺数颗粒湍流，Rep～O(100)。所以，IBM 方法可以用来研究高雷诺数颗粒湍流。

2 单个中性颗粒在驻点流中的动力学

驻点流的理论模型Hiemenz 流动可以在三维或者三维轴对称的条件下，通过自相似解求得[18]。Hiemenz 流动是全场非均匀的，远场速度大，近壁面速度小，压力则是远场小，近壁面大，这是因为流体动能转化为压力的缘故。远离壁面，壁面法向距离大于大约三个Hiemenz 特征长度δ，也叫黏性边界层厚度，Hiemenz理论解求出的涡量几乎为零[18]，流场是无黏的，是无黏挤压流动，挤压率可表示为B，根据具体问题不同，挤压率可以横跨好几个数量级B=O(1)～O(104)。在化工流动中，一般挤压率B～O(1)，以T 型管流为例，B=U/D，U是垂直管流的平均速度，D是等直径T型管的直径。而在航空航天飞行器的头部驻点流中，挤压率可以达到B～O(104)，B=U/D，U是脱体激波后的流体速度，D是脱体激波的距离。靠近壁面，由于壁面的无滑移滞止作用，黏性边界层会发展起来，Hiemenz理论解求出的涡量几乎为零[18]，Hiemenz流动中的黏性边界层厚度是常数v是流体的动力黏度，Bδ是流体的特征速度。因为Hiemenz流动是全场非均匀的，所以它没有特征流动雷诺数，只有局部雷诺数Rer=(r/δ)2，Rez=(z/δ)2，所以，通过局部最大雷诺数来控制数值计算域，防止数值发散[4]。防止数值发散至少有好几个方法，目前已知的一个是局部最大雷诺数不能太大，一个是远场采用非均匀粗网格。

图4表明，Li[4-5]的数值方法能很好地重构一个Hiemenz流场并与理论解吻合。图4 的横轴分别表示归一水平速度和法向速度，竖轴表示归一化的壁面距离，黑色虚线表示Hiemenz流动的理论解，从左到右不同彩色实线表示不同径向位置切面r/δ=1.56，3，8的DNS计算的速度剖面结果。Li[4]在此基础上研究单个有限尺寸中性悬浮颗粒在驻点流中的动力学（图5）。

2.1 与壁面发生弹性碰撞前

图6 显示了不同有限尺寸颗粒受到的无黏环境压力的理论解，Ba2是驻点流的特征水动力[19]，见式(3～5)。图6中横轴表示不同有限尺寸中性悬浮颗粒在驻点流中受到的归一化无黏环境压力的理论解，竖轴表示归一化的壁面距离，从左到右不同的彩色线表示从小到大的有限尺寸中性悬浮颗粒a/δ=0.8，1.6，2.4，3.2。可以看到，如果有限尺寸惯性颗粒只受到无黏环境压力作用，那么颗粒的运动导数和流体微团的随体导数是等价的，那这种情况可以认为颗粒有示踪粒子行为[20]。

图4 DNS流场归一化速度与驻点流理论解验证Fig. 4 Comparisons of dimensionless velocities between DNS and the theoretical solution

图5 一个有限尺寸中性悬浮颗粒在三维轴对称驻点流中的动力学的颗粒解析的直接数值模拟PR-DNS示意图[5]（ω/B表示归一化标量涡量）Fig. 5 A finite-size neutrally-buoyant particle in an axisymmetric stagnation point flow obtained by the particle-resolved DNS[5](ω/B represents thedimensionless vorticity intensity)

图6 不同有限尺寸颗粒受到的无黏环境压力的理论解[4]Fig. 6 The dimensionless theoretical ambient pressure force excerted on different finite size neutrally-buoyant particle[4]

图7（a）显示了不同有限尺寸颗粒沿法向对称轴受到的总水动力。远离壁面，总水动力线性递减，只受到无黏环境压力，呈示踪粒子行为。靠近壁面，由于壁面黏性边界层与颗粒相互作用，颗粒受到的总水动力不再是无黏环境压力，示踪粒子行为不再成立，见图7（b）。

图7 在不同壁面距离条件下，不同惯性颗粒受到的总水动力以及速度与流体微团速度的比较[4]Fig. 7 The hydrodynamic force on inertial particles at different gaps between the bottom of particles and the wall , and the velocities of inertial particles as a function of the distances between the center of particles and the wall [4]

图7（a）的横轴表示不同有限尺寸颗粒在驻点流中受到的归一化总水动力，竖轴表示归一化的壁面距离，从左到右不同的彩色线表示从小到大的有限尺寸中性悬浮颗粒a/δ=0.8，1.6，2.4，3.2。图7（b）的横轴表示不同有限尺寸颗粒在驻点流中的归一化速度，竖轴表示归一化的壁面距离，蓝色实线表示流体示踪粒子，其他的彩色线与图7（a）意义相同。图7上图表明靠近壁面，不同惯性颗粒受力趋势分为两种情况，一种情况颗粒受力趋于零，另一种趋于发散。这是因为近壁面的颗粒受力主要由润滑力控制，，当颗粒惯性小的时候，润滑力使得惯性颗粒动能耗散为零V～0，润滑力也收敛为零Flub～0；另一种情况，润滑力呈现发散趋势，因为惯性作用还不能在有限壁面距离内被流体的黏性作用完全耗散，颗粒以有限大小速度撞击壁面，而 ε～0，所以呈发散趋势，这里ε 是颗粒与壁面的无量纲壁面距离。这种情况下，弹性碰撞将会发生，纯粹描述流体动力学的N-S方程将不再能描述这个复杂物理问题，物理问题将变为流固耦合问题，必须考虑固体弹性碰撞的过程。

2.2 与壁面发生弹性碰撞后

当弹性碰撞发生后，方程（6）描述了颗粒动力学，方程（1）描述了有颗粒扰动的流体力学场。理论上是，求解一个考虑颗粒相动量注入的流体力学偏微分方程，同时求解一个颗粒动力学方程.之前的难点在于如何给出Fc的表达式。因为颗粒弹性碰撞力是物性材料和碰撞时间的函数，这里简要表示为Fc=f(Tc)，Tc是理论碰撞力表达式的理论碰撞时间。一般地，当颗粒惯性大的时候 St～O(1000)，颗粒反弹系数和流场特性无关，近似等效于颗粒本身的干碰撞系数[21-23]。但是当颗粒惯性小的时候 St～O(10)，反弹系数就和颗粒惯性相关，不再直接等于干碰撞系数，这是因为水动力效应黏性耗散了部分颗粒动能[21-22]。既然如此，当颗粒惯性小的时候，颗粒碰撞时间也没有理由是常数[24]。Li[5]总结了前人的实验和理论结果（ Zenit[23]和Birwa[24]）后，采用弹性力学解析解（Goldsmith[25]）来表达数值模拟中的Tc，获得了很好的结果。

图8表明，Li[5]的数值方法能很好地重构湿润碰撞的物理场景，准确预测了湿润碰撞的反弹系数，并且反应湿润碰撞的不同惯性颗粒的不同碰撞时间，湿润碰撞时间不取决于任意的数值光滑参数，数值上具有鲁棒性。图8（a）的横轴表示归一化的不同颗粒惯性，竖轴表示归一化的反弹速度，黑白符号代表实验数据（Joseph[21]和 Gondret[22]），彩色符号代表不同数值参数下的直接数值模拟DNS结果。图8（b）的横轴表示归一化的不同颗粒惯性，竖轴表示归一化的颗粒湿润碰撞时间，插图的竖轴表示归一化的干碰撞时间，彩色符号代表不同数值参数下的直接数值模拟DNS结果。Nc表示数值光滑参数，η表示无量纲归一化粗糙度，Stc表示不同粗糙度 η条件下的颗粒反弹发生的临界颗粒惯性。在化工领域，颗粒碰撞的典型惯性一般为是颗粒在静止环境流体中的定常沉降速度。所以Li[5]的数值方法是有针对性的。

图8 PR-DNS数值模拟的重力沉降问题的反弹系数和湿碰撞时间[5]Fig. 8 PR-DNS results of rebound coefficients and wet collision time for settling problem[5]

Li[5]给出了不同惯性的有限尺寸中性颗粒的反弹系数，数值实验表明，中性颗粒在驻点流中的动力学与沉降颗粒类似，都是依赖于颗粒惯性；数值方法上，湿润碰撞时间反应湿润碰撞的不同惯性颗粒的不同碰撞时间，湿润碰撞时间不取决于任意的数值光滑参数，数值上具有鲁棒性，见图9。这里Nc表示数值光滑参数，η表示无量纲归一化粗糙度，ac表示不同粗糙度 η条件下的颗粒反弹发生的临界颗粒半径。图9（a）的横轴表示归一化的不同颗粒惯性，竖轴表示归一化的反弹速度，彩色符号代表不同数值参数下的直接数值模拟DNS结果。图9（b）的横轴表示归一化的不同颗粒惯性，竖轴表示归一化的颗粒湿润碰撞时间，彩色符号代表不同数值参数下的直接数值模拟结果。我们发现，中性颗粒反弹的阈值、转捩点在a/δ～1处，这是合乎常理的。因为，一开始Li[4]把数值实验设置为中性颗粒，移除了两相密度比对颗粒惯性的影响，换句话说，颗粒惯性完全来自于颗粒本身的有限尺寸。如果颗粒有限尺寸与边界层厚度相当，那么颗粒惯性动能将会完全被黏性边界层耗散掉，见图7、图9和图10。如果颗粒有限尺寸比边界层厚度大很多，那么颗粒惯性将无法在有限截断尺度的壁面距离内被黏性边界层完全耗散（具体是润滑力主导这个过程），那么惯性颗粒将会以有限速度撞击壁面，见图7，弹性碰撞和反弹将会发生，见图9。总之，单个有限尺寸中性颗粒在驻点流中的动力学将完全由颗粒惯性a/δ控制[4]。

图10 不同有限尺寸中性颗粒与黏性边界层厚度对比示意图，a）和b）分别表示不同的有限尺寸中性悬浮颗粒在三维轴对称驻点流近壁面对流场的影响[4]Fig. 10 A schematic of surrounding flow fields around neutrallybuoyant particles with different sizes. Dashed lines represent the boundary layer edge [4]

3 一对中性颗粒在驻点流中的动力学

实际问题中，颗粒往往是一群的，它们聚集或者分散，因此研究一对中性颗粒的行为将为我们研究大量弥散颗粒在驻点流中的动力学打下研究基础（见图11（a））。在确保网格无关和时间步无关的前提下，我们测试了不同的惯性颗粒和不同的初始颗粒间距，确保图11（b）的结果是可重复的[5]。图11（b）的横轴表示归一化的时间，竖轴表示不同颗粒的归一化壁面距离，蓝色虚线代表单个中性颗粒，蓝色实线代表两倍于蓝色虚线的颗粒半径的中性颗粒，红黑线代表分别代表一对颗粒中的低位和高位颗粒，红色实线代表与蓝色虚线半径相同的低位颗粒，黑色实线代表与蓝色虚线半径相同的高位颗粒，绿色符号表示低位颗粒的无壁面接触的反弹点。我们发现，一对颗粒中的低位颗粒（红色轨迹）可以不接触壁面发生反弹，见图11（b）和插图里的绿色符号，而单个颗粒却没有这种行为（蓝色虚线轨迹）。这很反常，仅凭借水动力效应就能把一个惯性颗粒弹起吗？

图11 驻点流中的颗粒非接触碰撞反弹现象[5]Fig. 11 Contactless bouncing phenomenon of particle in a stagnation point flow[5]

图12解释了这种反常的反弹。图12的水平轴表示归一化时间。图12（a）的竖轴表示不同颗粒的归一化壁面距离，图12（b）的竖轴表示不同颗粒受到的总水动力，图12（c）的竖轴表示不同颗粒受到的无黏环境压力，图12（d）的竖轴表示不同颗粒受到的总水动力减去无黏环境压力的净效应。蓝色虚线代表单个中性颗粒，红黑线代表分别代表一对颗粒中的低位和高位颗粒，红色实线代表与蓝色虚线半径相同的低位颗粒，黑色实线代表与蓝色虚线半径相同的高位颗粒。简单来说，就是一对颗粒的情况下，高位颗粒替低位颗粒阻挡了来自上方的驻点流向下挤压效应[26],而低位颗粒没有物体替它阻挡。驻点流的非均匀性会产生一个向上托举的无黏环境压力[19]，见方程（3），见图6。图12也通过颗粒解析的直接数值模拟的结果表明，低位颗粒的总水动力是向上的，而单个颗粒则是向下的。无黏环境压力始终存在于驻点流，由于一对颗粒的高位颗粒的庇护作用，使得低位颗粒可以被无黏环境压力托举起来[5]。

Li[26]还发现惯性小的颗粒会抱团撞击壁面，惯性大的则无法抱团，见图13，水平轴表示归一化的颗粒惯性；左，右竖轴分别表示单个颗粒和一对颗粒中的低位颗粒的归一化反弹系数。Li[26]认为这是驻点流的挤压效应与颗粒惯性竞争的结果，当颗粒惯性小的时候，挤压效应大于颗粒惯性，颗粒间相互碰撞后无法弹开，被挤压流牢牢抱紧在一起；惯性大则不会。Li[26]归一化了抱团颗粒对的惯性后发现，抱团颗粒对的归一化反弹系数曲线和单个颗粒的反弹系数曲线是重合。这再一次印证了Li[4]的结论，驻点流中的颗粒动力学由颗粒惯性控制。

图12 单个和一对有限尺寸中性颗粒在驻点流中的受力以及各力分量比较示意图[5]Fig. 12 A illustration of the dynamics of a single and a pair of neutrally-buoyant particles in a stagnation point flow[5]

图13 单个（蓝色）和一对有限尺寸中性颗粒（红色）在驻点流中的反弹系数比较示意图[26]Fig. 13 The dimensionless rebound velocities of a single (blue)and a pair of (red) isolated neutrally-buoyant particles[26]

4 结论与展望

驻点流中的颗粒动力学与非驻点流颗粒动力学非常不同。研究通过选取有限尺寸中性悬浮颗粒在三维轴对称驻点流对称轴上运动这个基本模型，移除了真实来流扰动或者数值扰动的影响，移除了密度比对颗粒惯性的影响，这样颗粒惯性完全来自于颗粒有限尺寸，颗粒动力学将局限于法向，颗粒与环境流体和壁面的相互作用，不需要考虑来流扰动导致的非定常效应，不需要考虑颗粒不在对称轴上导致的剪切诱导升力。

远离壁面，中性颗粒会像示踪流体微团一样运动。通过数学推导，研究得出，颗粒动力方程的牛顿惯性项等于流体微团的随体导数。这在数学上解释了：为什么有限尺寸惯性颗粒与数学上无穷小的流体微团示踪粒子的运动轨迹重合。

靠近壁面，颗粒惯性的不同将导致颗粒偏移流体示踪颗粒轨迹的幅值不同: 惯性大的颗粒，将会诱导出大的润滑力，反之亦然。但是润滑力是数学理想模型。当颗粒惯性小，黏性润滑力足以耗散掉所有的颗粒动能，颗粒将会减速至零并停留在壁面；当颗粒惯性大，在撞击有限尺度的粗糙度壁面前，润滑力不足以耗散掉颗粒动能，颗粒将会以有限速度撞击壁面，弹性碰撞将会发生，这个有限速度与颗粒惯性成正比。当弹性碰撞发生后，流体力学问题变成了流固耦合问题。因此我们耦合一个固体弹性力学的解析表达式到碰撞模型里，与流体求解器耦合求解，与经典沉降实验数据吻合。

研究发现，单个中性颗粒在驻点流中的反弹系数与沉降颗粒类似，与颗粒惯性相关。但是当研究到一对中性颗粒的时候，反常的现象发生了。首先，根据颗粒惯性不同，颗粒对有可能会抱团，这是驻点流挤压效应与颗粒惯性相互竞争的结果；其次，当颗粒对以抱团撞击壁面的时候，颗粒的反弹系数与单个颗粒反弹系数是重合的；最后，当颗粒对不抱团的时候，低位颗粒可以不撞击壁面发生反弹。

研究阐明了驻点流中的有限尺寸中性颗粒动力学完全由颗粒惯性/有限尺寸控制的机理,为进一步研究驻点流中的颗粒动力学铺垫了基础。对于某些工业应用中的驻点流中的颗粒动力学问题，颗粒运动轨迹往往不是在对称轴上的，那么随之而来的平行加速流带来的剪切诱导升力就必须要考虑；此外，颗粒与携带流体的密度比也不总是1，颗粒的惯性将不再由有限尺寸单独控制，颗粒相对流体的密度比效应也需要考虑。未来我们将针对这些真实问题带来的复杂效应，进行更深入的基础研究。

致谢感谢国家重点研发计划（2019YFA0405200）经费的支持，感谢图卢兹大学Abbas 教授、Climent教授、Magnaudet教授和纽约城市大学Morris教授在作者博士工作期间给予的帮助和支持。

APPENDIX: NOMENTCLATURE

ρfcarried fluid density, kg/m3

ρpparticle density, kg/m3

Npparticle number, 1

aparticle radius, m

dparticle diameter, m

Ufunperturbed fluid velocity at particle center, m/s

Vpparticle velocity vector at particle center, m/s

Vnormal particle velocity at particle center, m/s

gnormal gravity constant, m2/s

qparticle flux over unit surface, m2/s