杨海龙，张庆河，谢 琳

(天津大学 水利工程仿真与安全国家重点实验室，天津 300072)

传统泥沙运动力学主要是通过各种推移质、悬移质输沙公式或方程来描述泥沙输移和地形演变[1-2]，对局部冲刷等较为剧烈的泥沙运动进行数值模拟研究，大多也是将床面剪切应力作为泥沙起动的主要动力因素，并将其耦合到泥沙冲淤模型中，无法揭示水流泥沙运动的内在机制[3]。近年来，随着计算机计算能力的提高，从细观尺度研究泥沙运动规律的工作越来越多，例如，以DPM(Discrete Particle Method)为代表的离散颗粒模型，能够通过对单个颗粒运动进行模拟来分析水流泥沙相互作用机理[4]。DPM方法可以描述颗粒运动过程，模拟精度较高，但存在计算量巨大，目前只能模拟颗粒数为几万或者几十万级的问题，而对于百万甚至千万级的大规模颗粒运动还很难满足计算要求，将其应用于模拟实际的工程问题尚存在困难。

MP-PIC(MultiPhase Particle-In-Cell)称为多相流质点网格模型[5]，是一种新型的离散颗粒模拟方法，近年来在化工等领域得到了广泛应用。与传统的DPM方法不同，该方法将离散的颗粒碰撞力视为一个连续的各向同性的颗粒相应力梯度力场，该梯度应力的大小与颗粒体积分数密切相关，在实际计算过程中不需要进行颗粒邻域搜索，无须计算颗粒间的接触力。另外，该模型可以采用一个“代表颗粒”来代表具有相同动力学特性的颗粒群，因此能够极大地节省计算量，提高颗粒模拟的计算效率，适用于模拟百万甚至千万级的大规模的颗粒系统[6]，有望应用于实际工程问题的模拟。为此，本文将基于MP-PIC方法，采用欧拉方法描述液相，用拉格朗日法描述颗粒相运动，建立描述泥沙静水沉降的三维数值模型，并采用该模型对泥沙静水沉降进行模拟研究。

1 模型的建立与计算方法

1.1 流体控制方程及拖曳力模型

1.1.1 流体运动控制方程

MP-PIC模型中，考虑颗粒组分影响的三维流体运动的控制方程由以下两相流连续性方程和时均Navier-Stokes方程组成[7]

(1)

(2)

(3)

式中：f为计算颗粒运动函数；mp为颗粒质量；up为颗粒运动速度；uf为颗粒位置处的流体速度；Dp为拖曳力系数；P为颗粒位置处的流体压力梯度；ρp为颗粒密度。

1.1.2 紊流模型

为了使式(3)在求解时封闭，需要引入描述紊动粘性系数μt的紊流模型。MP-PIC模型采用了大涡模拟(LES)，它是介于直接数值模拟(DNS)与雷诺平均法(RANS)之间的一种紊流数学模型。此方法通过滤波函数将紊流运动中小尺度的涡动滤掉，而小尺度涡对大涡的影响则通过向方程中加入亚格子尺度模型(SGS)进行描述。本文采用Samgorinsky[8]提出的SGS模型，其μt系数如下

(4)

Cs=Cs0(1-ey+/A+)

(5)

△=(δxδyδz)1/3

(6)

式中：Cs为Samgorinsky常数；△为亚格子过滤尺寸；Cs0=0.1为Van Driest常数[9]；y+是到壁面无量纲化的最近距离；A+=25.0是半经验化的常数；δi表示沿i轴方向的网格尺度。

1.1.3 拖曳力模型

本文采用由Gidaspow[10]提出的Wen-Yu[11]与Ergun[12]混合拖曳力模型，该模型由各分段函数组成，公式如下

(7)

本文中，αcp是堆积状态下的颗粒体积分数，一般取0.6，αp为沉降过程中的颗粒体积分数，D1为Wen-Yu拖曳力模型，D2为Ergun模型，公式如下

(8)

(9)

式中：C1=180；C2=2；Cd为Wen-Yu拖曳力系数；rp为颗粒半径；uf为颗粒位置处流体的流速。

(10)

(11)

1.2 颗粒运动方程及数值求解

1.2.1 颗粒运动方程

在MP-PIC方法中，固相采用离散颗粒方法，考虑颗粒的拖曳力、重力、颗粒相互碰撞的接触力等。颗粒相的运动通过求解关于颗粒分布函数(PDF)f的Liouville输运方程[13]得到。由于颗粒不存在化学反应，故在这里并不考虑颗粒的温度和材料的变化。其中，f为颗粒空间位置xp、速度up、质量mp以及时间t的函数

f=f(xp,up,mp,t)

(12)

f的输运方程源自于气体动力学理论中对碰撞项松弛时间有所改进的Boltzmann-BGK方程[14-15]

(13)

(14)

式中：Ap为颗粒加速度；up为颗粒速度散度；为碰撞源项，由碰撞阻尼力和恢复各向同性力两部分组成；fD和fG分别是速度为质量平均态和高斯分布态的颗粒分布函数；τD和τG为达到这两种局部平衡态的颗粒碰撞松弛时间，与颗粒的碰撞恢复系数有关。由此得出，颗粒运动方程为

(15)

(16)

式中：Xp为颗粒碰撞力，其表达式为

(17)

(18)

αp+αf=1

(19)

式(16)中右边第一项为运动状态下初始碰撞力，第二项为颗粒堆积状态下的平均碰撞力，只有当颗粒出现堆积状态时，该项才会生效。τp为颗粒碰撞正应力，是颗粒周围的多个颗粒对其的共同作用力，由欧拉网格计算得到的空间梯度应力来表征，而梯度应力通过插值到网格单元的颗粒体积分数得到

(20)

式中：Ps为常数，是模型的参数，具有压力的量纲；αp是颗粒的体积分数；αcp为颗粒堆积临界体积分数；β和ε为无量纲常数，一般建议取2≤β≤5，ε取10-7量级的小数。由于颗粒的堆积状态与颗粒的尺寸、形状及排列方式有关，该模型允许颗粒的体积分数达到或偶然超过临界体积分数，所以引入ε用来消除奇点，以保证颗粒碰撞应力的求解不受影响。

1.2.2 颗粒与壁面间的相互作用

在MP-PIC方法中，颗粒与壁面间的相互作用并不需要对颗粒的受力进行分析，而是引入2个系数[16]，分别为颗粒与壁面碰撞的正向恢复系数en,w和切向恢复系数et,w，由此得出颗粒与壁面碰撞后的运动轨迹。

uc,i,n,post=-en,w×uc,i,n,pre

(21)

uc,i,t,post=-et,w×uc,i,t,pre

(22)

式中：uc,i,n,pre和uc,i,n,post是碰撞前后的法向速度；uc,i,t,pre和uc,i,t,post是碰撞前后的切向速度，下标i表示一个计算颗粒。

1.2.3 颗粒运动的数值求解

在数值计算中，流体的三维连续控制方程采用有限体积法求解，颗粒的运动并不是直接求解Liouville输运方程，而是将颗粒的性质插值到计算网格单元中心，或者通过交错的标量与动量节点将流场的性质插值到颗粒位置处，使得流体与颗粒的控制方程通过动量交换隐性地耦合到一起，进而通过有限差分的方法求得代表颗粒的位置和速度并不断更新[17]

(23)

(24)

图1 颗粒与流体运动的耦合过程Fig.1 Procedure for coupling of particles and fluid

1.3 颗粒与流体运动的耦合

本文采用交错网格，并通过插值算子将颗粒的属性插值到欧拉网格当中，首先求解颗粒与流体之间的动量交换项，然后流体运动采用OpenFOAM中基于PISO算法的瞬态不可压缩流体求解器进行求解，再次应用插值算法将流场的特性插回到颗粒位置处，最后求解颗粒的运动方程，从而实现颗粒与流体运动的耦合。具体计算过程如图1所示。

2 颗粒沉降模拟的网格敏感性分析

在应用MP-PIC方法进行颗粒静水沉降模拟之前，本节首先通过模拟单个颗粒泥沙在静止水体中自由沉降过程，将颗粒沉速结果与经验公式对比并进行分析，讨论模型对计算单元网格的敏感性，并对网格尺度的选取范围给出合理建议。

2.1 几何模型及计算参数

表1 模型中颗粒与流体参数Tab.1 Particle and fluid parameters used for the present simulation

本节采用与苏东升[18]基于DPM方法进行颗粒沉降模拟一样的算例设置，模拟区域为长方体，坐标原点位于模拟区域底边界，z轴为高度方向，计算网格采用三维正方体结构化网格。模拟初始时刻，颗粒位于模拟区域中轴线上，距底面一定距离从静止状态释放，流体初始时刻保持静止，流速为零。开始阶段，泥沙颗粒在重力作用下加速下沉，随着颗粒速度不断增大，颗粒所受拖曳力不断增大，颗粒做加速度不断减小的变加速运动，当颗粒所受的拖曳力与重力保持平衡时，达到稳定沉降速度并保持匀速下沉。在模拟中要保证颗粒达到沉降速度时与计算底面仍有足够的高度，与x轴和y轴方向壁面拥有足够的宽度，排除模型壁面对颗粒沉速的影响。模拟中，颗粒形状简化为球体，本节选取颗粒直径为0.1 mm和2 mm两组算例，每一组算例分别设置计算网格长度Δx与颗粒直径D的比例为2、5、10、20、40五种情况，共计10个算例。其中，粒径为0.1 mm的算例采用层流模型，粒径为2 mm的算例采用LES紊流模型。当Δx/D小于40时，计算域长度和宽度均取为颗粒粒径的100倍，高度取为粒径的500倍，当Δx/D等于40时，计算域长度和宽度均取为颗粒粒径的120倍，高度取为粒径的520倍，以保证泥沙颗粒在沉降过程中不受壁面的干扰。模拟中采用的流体为清水，颗粒为泥沙，所有算例采用的颗粒与流体参数均相同，具体参数如表1所示。

2.2 模拟结果分析

Cheng[19]基于Brown等[20]收集的大量实验资料提出了计算泥沙颗粒无量纲沉速经验公式

(25)

式中：w*=((ρp-ρf)/ρf·gv)-1/3，w是无量纲沉速；dp*=((ρp-ρf)g/(ρfv2))1/3，dp是无量纲粒径。这里，将所有算例模拟结果与Cheng经验公式计算值的比较结果绘制于图2中。

由图2可以看出，大部分模拟结果比经验公式结果偏小，且随着网格尺度与粒径之比的增大，模拟结果偏小幅度增大，分析原因可能与本文选取的拖曳力模型有关，导致拖曳力系数选取稍大，而且随着网格越来越粗，计算精度降低，导致结果相差越来越大。

但总体而言，模拟沉速值与经验公式计算结果相差不大，对于粒径较小的0.1 mm颗粒，比较结果随网格尺度与颗粒直径之比Δx/D增大而有所增大，但在Δx/D小于40的条件下，相差不超过3.5%。对于粒径较大的2 mm颗粒，模拟沉速值与经验公式计算结果的比较结果同样随Δx/D增大而增大，但是当Δx/D≥5时，相差变化不明显，控制在1.7%以内，网格敏感性大幅降低。综合考虑网格对粗细颗粒的模拟精度以及计算量的问题，同时大涡模拟要求计算网格不能过粗，推荐采用2～5倍泥沙粒径模拟泥沙颗粒运动。

3 颗粒静水沉降模拟

为说明采用MP-PIC方法能够对水流泥沙相互作用机理进行分析并能够将其用于模拟和研究泥沙沉降的过程和现象，本节主要对不同粒径单颗粒沉速和群体颗粒制约沉降进行模拟研究。

图3 模拟结果与实验值和经验公式的比较Fig.3 Comparison among simulated results， experimental data and empirical formula

3.1 不同粒径单颗粒沉降模拟

模拟中，选取泥沙粒径在0.05～10 mm之间的14个算例，覆盖了层流、紊流过渡及紊流3个沉降区，每个算例的模拟区域长宽为颗粒粒径的100倍，高度为粒径的500倍，以保证泥沙颗粒在沉降过程中不受壁面的干扰。初始时刻，流体流速为零，颗粒位于计算域的中轴线，并由静止释放自由沉降。颗粒和流体物理参数与2.1节所述相同，网格尺度为2倍颗粒粒径。

将模拟结果与Brown等人实验值与Cheng经验公式一起绘于图3中。由图可知，不同粒径单颗粒泥沙沉速模拟结果与已有实验结果及经验公式值较为吻合，证明该模型可以合理模拟较大范围的单颗粒沉降。

将模拟结果与Cheng经验公式的比较结果列于表2中。由表中可得，模拟结果与经验公式的结果相差基本控制在6%以内，模拟结果较为理想，满足一般的模拟精度要求。

表2 模拟值与经验公式计算值的比较结果Tab.2 Comparisons between simulation results and empirical formula calculations

3.2 群体颗粒沉降模拟

3.2.1 模型及参数设置

本节参考Peysson等人[21]的物模实验，颗粒与流体参数均与物模实验保持不变，模拟不同浓度泥沙颗粒的制约沉降过程，并与物模实验进行对比验证。本模型计算域在xyz方向的坐标范围为[0，0.04]m×[0，0.04]m×[0，0.05]m，计算网格采用三维正方体结构化网格，模拟区域边界均为光滑壁面。初始时刻，泥沙颗粒按5%和10%的体积浓度随机分布在计算域一定高度处的正方体区域内，范围为[0，0.04]m×[0，0.04]m×[0，0.04]m，保证颗粒沉速分布达到稳定时与底面仍有足够高度，颗粒由静止自由沉降。颗粒及流体具体参数如表3所示。

表3 模型中颗粒与流体参数Tab.3 Particle and fluid parameters used for the present simulation

3.2.2 群体颗粒沉降模拟结果

按Peysson实验中的方法，用x=0.04n(m)(n=1,2,…,9)9组平面沿x轴将模拟区域平均分为10个体积相等的区间，待颗粒沉速达到稳定之后，统计每个区间内颗粒的平均沉速，从而得到颗粒沉速随x轴的分布。这里将体积浓度为5%和10%两组算例的的泥沙颗粒沉速结果与Peysson实验值绘制于图4中。

4-a Φ=5%4-b Φ=10%图4 泥沙颗粒沉速沿x轴的分布Fig.4 Distribution of settling velocity of sediment particle along x-axis

图4中实心点为每个区间内沿z轴的平均速度w‖(x)与整个模拟区域内沿z轴的平均速度w‖的比值；空心点为每个区间内垂直于z轴的平均速度w⊥(x)与整个模拟区域内沿z轴的平均速度w‖的比值。由图可知，模拟中间区域颗粒沉速较大，靠近边界壁面沉速较小，与群体颗粒沉降存在内部对流的物理现象较为符合。而且颗粒主要沿z轴方向运动，垂直于z轴方向速度较小。整体沉降速度分布趋势与Peysson物模实验大致相同。

由于沉降时颗粒之间会产生碰撞和干扰，群体颗粒沉降速度小于单颗粒沉速，并且随着颗粒浓度的增加，群体颗粒沉速越来越小。Richarson和Zaki[22]利用量纲分析法得到了群体沉速的计算公式为

wg/w0=(1-Φ)m

(26)

式中：wg与w0分别为群体颗粒及单颗粒的沉速；Φ为泥沙颗粒的体积浓度；m为幂指数，本文与Peysson相同取为5。将5%和10%的模拟结果与理论公式，实验值绘制于图6，发现MP-PIC方法模拟结果与Peysson实验值有一定偏离，但偏离均在实验结果的误差线内。模拟结果与理论公式的误差随体积浓度的增大而减小，当体积浓度Φ=5%时，二者相差3.49%，当体积浓度Φ=10%时，相差仅为0.17%。同时，模拟结果的颗粒沉速同理论公式一样具有随体积浓度增大而减小的趋势，能够较好地模拟群体颗粒沉降。

图5 对流强度随体积浓度的变化Fig.5 Variation of amplitude of the intrinsic convection with volume concentration图6 群体颗粒沉速与其他结果对比Fig.6 Comparison between population particle settling velocity and other results

4 结论

本文基于MP-PIC方法，采用欧拉方法描述液相，用拉格朗日方法描述颗粒相运动，建立了描述泥沙静水沉降运动的数值模型，分别对网格敏感性、单颗粒沉速和群体颗粒制约沉降等问题进行研究，可以得出：

(1)该模型具有对水流泥沙相互作用机理进行分析的能力，能够用于模拟和研究泥沙沉降的过程和现象。

(2)通过模拟单个泥沙颗粒在静止水体中自由沉降，讨论模型对计算单元网格的敏感性问题。为尽量缩小模拟结果与经验公式结果之间的差距，但是又不能使得计算量太大，同时考虑到大涡模拟对计算网格不能过粗的一般要求，推荐采用2～5倍泥沙粒径的网格尺度来模拟泥沙颗粒运动。

(3)将不同粒径单颗粒泥沙沉速与Brown等人实验值、Cheng经验公式进行对比，结果较为吻合，证明MP-PIC模型可以较好模拟不同粒径单颗粒泥沙运动。

(4)模拟了不同体积浓度泥沙颗粒的制约沉降过程，结果显示，模拟中间区域颗粒沉速较大，靠近边界处沉速较小，并且垂直于沉降方向的速度接近为零，与Peysson物模实验较吻合，充分体现了群体颗粒的沉降对流现象；同时，群体颗粒沉速模拟结果与理论公式相差不大，具有随体积浓度增大而减小的趋势，采用MP-PIC方法能够很好地对泥沙颗粒静水沉降进行三维数值模拟研究。