垂直管内固-液两相流全耦合CFD-DEM模型研究

2022-02-08张德胜赵睿杰施卫东

农业机械学报 2022年12期

张德胜周游赵睿杰施卫东

(1.江苏大学国家水泵及系统工程技术研究中心，镇江 212013； 2.南通大学机械工程学院，南通 226019)

0 引言

随着我国城市化发展越来越快，基础设施中排污设施的建设也越来越多[1-3]。污水管网作为城镇、工业和农业区的重要排污设施，对当地生活和经济发展起着重要的保障作用。

生活污水和工业废水中往往含有纤维杂质、纸屑和颗粒等固体物，其中较轻的固体漂浮在水面并随污水漂流，较重的固体分布在水流断面上并呈悬浮状态流动，造成固体物之间发生缠绕、摩擦和碰撞，或者与管道壁面发生摩擦碰撞，导致管道的磨损和堵塞，使得整个排污系统无法正常工作[4-6]。因此，对污水污物的两相流动展开研究，可对排污系统的设计、安全可靠的运行以及后期的维护起到重要的指导作用。

固-液两相输送系统的流态极其复杂，受各工况参数影响[7-10]，其中流体的流动特性、固体物的运动特性以及流体-固体之间复杂的耦合机制受到广泛关注。目前在固-液两相流领域的主要研究方法包括实验研究和数值模拟。过去，许多学者基于实验数据将压降或者其他过程参数与操作参数(如颗粒体积分数、输送速度和颗粒直径)关联起来[11-13]，但这种方法只在特定条件下才具有适用性。此外，这种分析结果一般为流体的宏观信息，无法通过其深入了解固-液两相流中的微观信息。随着计算机技术的发展，数值模拟已经成为研究固-液两相流的另一个重要手段。数值模拟方法可以分为Euler-Euler方法[14-17]和Euler-Lagrange[8-9,18]方法。在前者中，液相和颗粒相均被视为连续介质，求解时均Navier-Stokes方程，但该方法很大程度上依赖于本构关系，不能获得足够的关于颗粒-颗粒/壁面碰撞的信息。而在Euler-Lagrange方法中，流体被视为连续相，直接求解时均Navier-Stokes方程，而颗粒被视为离散相，通过计算流场中每个颗粒的运动，从而在仿真过程中能够很好地捕捉颗粒的碰撞信息。流体和颗粒之间的相互作用通过在各个方程中添加源项实现。

在实验方面，文献[7]利用激光多普勒测速仪(LDA)测量了垂直管道内固-液两相的速度分布，并且发现重颗粒趋向于集中在管道中心，而轻颗粒趋向于集中在管壁附近。在数值模拟方面，文献[8]以NaCl颗粒和盐水为载体进行了参数化研究，考察操作条件对水平管中NaCl颗粒两相流动的影响。文献[9]采用CFD-DEM (Computational fluid dynamics-discrete element method)方法，对垂直管内不同参数对颗粒水力输送的影响进行研究，但未考虑颗粒对湍流的影响。文献[16]利用Euler-Euler方法对颗粒-液浆流动进行模拟，结果表明，随着颗粒尺寸增大，摩擦压力损失减小，颗粒体积分数的径向分布增大，且管道底部颗粒体积分数增高。文献[19]比较具有不同叶轮型式潜水排污泵的通过性能，发现双流道污水泵内颗粒速度更快，颗粒聚集程度更小，泵的通过性能强于双叶片污水泵。文献[20]对后掠式双叶片污水泵固液两相流场进行数值模拟，分析了固相颗粒的分布情况和泵内过流部件的磨损情况，并探讨了颗粒体积分数及颗粒大小对泵水力性能的影响规律。

目前大部分模型采用传统Saffman和Magnus升力模型，且忽略其他次要力以及颗粒-湍流的相互作用。本文提出一种全耦合CFD-DEM模型，考虑所有次要力和颗粒-湍流的相互作用，分析进口颗粒体积分数、输送速度和颗粒直径对固-液两相流动的影响。

1 数学模型

1.1 流体控制方程

流体被视为连续相，其控制方程是质量守恒和动量守恒方程，通过在动量方程中添加源项Fp-f来考虑颗粒对流体的影响，公式为

(1)

(2)

(3)

式中ρf——流体密度，kg/m3

t——时间，s

p——压力，Pa

μf——流体动力粘度，Pa·s

μt——流体湍流粘度，Pa·s

g——重力加速度，m/s2

Fp-f——颗粒对流体的力源项，N/m3

Ff-p——流体对颗粒的总作用力，N

Vcell——流体网格体积，m3

n——流体网格内颗粒数

αf——流体体积分数，%

颗粒被划分为64个部分，计算流体网格中的流体体积分数，公式为

(4)

式中k——与该网格重叠的颗粒数

Vi,part——颗粒i与该网格重叠的总体积，m3

采用标准k-ε湍流模型来封闭方程，通过在方程中添加湍动能源项和湍流耗散源项来考虑颗粒对湍流的影响[21]，公式为

(5)

(6)

(7)

(8)

(9)

式中Sk——湍动能源项，kg/(m·s3)

Sε——湍流耗散源项，kg/(m·s4)

vf——流体瞬时速度，m/s

vp——颗粒速度，m/s

dp——颗粒直径，m

νf——流体运动粘度，m2/s

fp——颗粒雷诺数有关的系数

Rep——颗粒雷诺数

C′ε3——颗粒雷诺数有关的系数

其中C′ε2=1.92,Cεp=0.085。

1.2 颗粒运动方程

颗粒运动采用离散元法(DEM)求解，粒子的运动由牛顿第二定律描述，公式为

(10)

(11)

式中mp——颗粒质量，kg

Fg——重力，N

Fdrag——Dif曳力[22]，N

Floth——Loth升力[23]，N

Fp——压力梯度力，N

Fvm——虚拟质量力，N

Fc——颗粒-颗粒/壁面间相互作用力，N

Ip——转动惯量，kg·m2

ωp——颗粒角速度，s-1

Tc——颗粒接触力矩，N·m

Tf——流体作用力矩，N·m

许多因素都会导致升力的产生，而流体涡度和颗粒旋转是产生升力的两个主要机制。为表示流体涡度和颗粒旋转对升力的影响，Loth升力模型将升力系数考虑为剪切和旋转引起的升力系数的组合，这扩大了升力模型的适用范围。根据文献[23]，升力表示为

(12)

式中Cl——Loth升力系数

ωf——流体涡度，s-1

为比较不同升力模型的效果，本文给出Saffman和Magnus升力的表达式。流体剪切运动引起的Saffman升力[24]Fsaf表示为

(13)

式中Cls——Saffman升力系数

根据文献[25]，升力系数Cls可表示为

(14)

(15)

颗粒旋转引起的Magnus升力计算公式[26]为

(16)

(17)

(18)

(19)

式中Rer——颗粒旋转雷诺数

ωf-p——相对角速度，s-1

当颗粒密度远大于流体密度时(ρf/ρp≪1)，虚拟质量力和压力梯度力并不重要(气-固流)。当ρf/ρp大于0.1时，虚拟质量力和压力梯度力变得显著，并且需要被考虑，即

(20)

(21)

式中Vp——颗粒体积，m3

Cvm——虚拟质量系数，默认值为0.5

1.3 离散随机游走模型

(22)

(23)

式中ξi、ξj、ξk——i、j、k方向的高斯分布随机数

u′、v′、w′ ——i、j、k方向的脉动速度分量，m/s

由于RANS模型的简单性，假设湍流为各向同性，计算脉动速度的均方根为

(24)

式中kf——流体湍动能，m2/s2

将上述流体的瞬时速度用于以上各种模型的计算中，所以在计算颗粒的运动时也考虑了湍流效应。

2 模型设置

2.1 数值格式

控制方程用有限体积法离散，采用标准k-ε湍流模型，采用SIMPLE算法处理速度-压力耦合。在CFD-DEM模拟中，颗粒时间步长Δtp通常为Rayleigh时间步长ΔtR的10%～40%，并且流体时间步长需为颗粒时间步长的10～100倍。综合考虑计算精度和时间成本，固体颗粒时间步长设为2×10-5s，流体时间步长设为2×10-4s。ΔtR计算式为

(25)

式中Ri——颗粒i半径，m

vi——颗粒i泊松比

Yi——颗粒i杨氏模量，Pa

2.2 模型设置及边界条件

建立与文献[7]实验相同的管道，并且进行网格划分，如图1所示。为更好捕捉壁面处流体信息，对边界层网格进行了细化，壁面第1层网格厚度为0.26 mm，且满足30

表1为模拟中的颗粒和流体属性[9]。表2为本文所有的模拟情况，研究了进口颗粒体积分数Cv、输送速度Uf和颗粒直径dp对固-液两相流动的影响。

表1 颗粒和流体属性Tab.1 Particle and fluid properties

表2 本文模拟中使用的参数Tab.2 Parameters used in this simulation

3 结果与讨论

3.1 模型验证

为检验模型的正确性，对直径为2.32 mm的球形颗粒进行数值模拟，并将结果与文献[7]的实验进行比较，图2(图中r表示测量点到管道中心的径向距离，R表示管道半径，并进行无量纲处理)所示为全耦合模型的验证。由图2可知，无升力的模型不能很好地预测管道中心处颗粒的分布，且管道壁面附近的颗粒体积分数略高于实验。传统的Saffman和Magnus升力模型预测的颗粒体积分数在管道中心和壁面处误差很大。而考虑Loth升力的模型能够较好地预测管道内颗粒体积分数的分布，同时各相轴向速度和流体均方根速度的模拟结果与实验吻合较好，表明考虑Loth升力和湍流调制的全耦合模型能够较好地预测颗粒的运动和流体的关键流动特性。

图2 模型验证Fig.2 Model validation

本文建立了网格数分别为6×104(粗网格)、9×104(中等网格)和1.1×105(细网格)的3套模型，对比3套网格预测的颗粒体积分数及各相轴向速度分布，发现粗网格和细网格与参考网格的误差均小于3.0%。通过对比不同模拟时间(4、5、6 s)下的计算结果，发现颗粒体积分数及各相轴向速度的误差均小于2.5%。所以，综合考虑计算精度和成本，本文采用中等网格和总模拟时间为4 s的方案来进行数值计算。

3.2 进口颗粒体积分数对水力输送的影响

进口颗粒体积分数Cv是固-液两相排污系统中最重要的参数之一，其直接影响颗粒的运动以及流体的流动特性。研究Cv为1.0%、2.5%和5.0%的影响，分别对应表2中序号1、2和3。

图3 不同Cv下颗粒分布Fig.3 Distribution of particles at different Cv values

不同Cv下颗粒的分布情况如图3所示。从图3中看出，当Cv为1.0%时，颗粒集中在管道中心处，而管道半径方向的中部(r/R=0.5)和壁面附近(r/R=0.9)颗粒存在较少，甚至没有。当Cv增加到5.0%时，颗粒同样主要集中在管道中心，但半径方向的中部和壁面附近均存在一定的颗粒。这是因为随着进口颗粒体积分数的增加，管道内颗粒的数量急剧增加，颗粒更加分散，同时更多的颗粒由于升力作用向管道中心聚集。

图4 不同Cv下碰撞次数Fig.4 Number of collisions at different Cv values

碰撞次数可分为颗粒间和颗粒-壁面间碰撞次数，在z=1.7～2.2 m区域内提取碰撞次数。不同Cv下，碰撞次数随时间变化如图4所示。可以看出，颗粒间碰撞从t=0.37 s左右才开始，这是因为颗粒需要经过一定时间运动到所研究区域内，而在t=0.5 s以后，颗粒间碰撞次数和颗粒-壁面间碰撞次数均已稳定，即颗粒两相流已经充分发展，考虑到稳定计算和其他工况，对2～4 s之间的数据进行平均处理(虚线所示)。从图4a、4b可以看出，颗粒间和颗粒-壁面间碰撞次数都随着Cv增加而增加，这是因为颗粒数量的增加使得更多颗粒向管道中心聚集，管道内发生更加频繁的碰撞。同时，管壁附近颗粒体积分数有所增加，该处的颗粒开始与壁面发生一定次数的碰撞。但颗粒-壁面间碰撞次数远低于颗粒间碰撞次数，且颗粒-壁面间碰撞次数的波动幅度随Cv的增加而增加。

图5 不同Cv下各物理量分布Fig.5 Distribution of various physical quantities at different Cv values

颗粒体积分数及各相的速度分布如图5所示。由图5可知，颗粒趋向于聚集在管道中心，颗粒体积分数在管道中心出现最大值，并沿管径向壁面逐渐降低，在壁面减为零。流体和颗粒的轴向速度在管道中心存在最大值，并且沿管径向壁面降低，而流体均方根速度趋势相反，其在管道中心处最小，并沿管径向壁面逐渐增大，到壁面附近处增加到最大值。

由图5a可见，随着Cv增加，整个管道内颗粒体积分数都增大，且管道中心处颗粒体积分数增加最为明显。这是因为Cv增大直接导致颗粒数量增多，更多的颗粒在升力作用下向管道中心聚集。而颗粒间碰撞频率随着颗粒数量的增加而显著增加，更多的碰撞使得颗粒更加分散，这也就能解释管道半径方向的中部和壁面附近颗粒体积分数的增加。

由图5b、5c可见，随着Cv增加，流体轴向速度有所降低，这种变化在管道中心区域更明显，在管壁附近忽略不计。这是因为随着Cv增加，流动状态从以流体为主导的流动模式转变为流体-颗粒作用和颗粒-颗粒作用为主导的流动模式，更多数量的颗粒使得流体-颗粒间的动量交换加剧，导致流体轴向速度降低，而靠近壁面区域颗粒较少，两相间动量交换较少使得靠近壁面区域的流体轴向速度几乎不变。随着Cv增加，颗粒轴向速度也出现了一定的降低，这是由于颗粒的动能被更加频繁的颗粒间相互作用所耗散。由图5d可以看出，颗粒对湍流的调制作用导致流体均方根速度有所降低，且调制作用随Cv增加更加明显。

图6 Cv对无量纲压降的影响Fig.6 Effect of Cv on dimensionless pressure drop

3.3 输送速度对水力输送的影响

输送速度Uf与输送系统(管道、过滤系统和泵系统等)的流量和效率密切相关，是两相流动系统中另一重要的参数。本节研究了Uf为1.5、3.0、4.5 m/s对水力输送特性的影响，分别对应表2中序号4、5和6。

不同Uf下颗粒分布情况如图7所示。从图7看出，在Uf=1.5 m/s时，颗粒在整个管道内分布最均匀，并且壁面附近处存在较多的颗粒。当Uf增加到3.0 m/s时，颗粒的分布不均匀，在管道中心处较为集中，且管道壁面附近的颗粒数减少。当Uf继续增加到4.5 m/s时，颗粒的分布更加不均匀，大部分颗粒集中在管道中心，壁面附近几乎没有颗粒。因此，随着Uf增加，管道壁面附近的颗粒逐渐减少，更多的颗粒向管道中心聚集，颗粒在管道内的分布逐渐不均匀。

图7 不同Uf下颗粒分布Fig.7 Distribution of particles at different Uf values

不同Uf下，碰撞次数随时间的变化如图8所示。可以看出，在不同Uf下，开始出现颗粒间碰撞和颗粒流达到稳定的时间不同。在较低的Uf下，开始出现颗粒间碰撞的时间和颗粒流达到稳定的时间都较长，分别为1.1 s和1.5 s。这是因为颗粒从管道入口运动到出口的时间由颗粒速度决定，较低的Uf导致颗粒速度较低，颗粒到达管道出口的时间就较长。从图8还可以看出，颗粒间碰撞次数和颗粒-壁面间碰撞次数都随Uf的增加而增加，这是因为随着Uf的增加更多的颗粒向管道中心聚集，中心区域内发生更为频繁和剧烈的碰撞，导致少部分颗粒在碰撞的作用下向管壁运动，并与管壁发生碰撞。同样，颗粒间碰撞次数远大于颗粒-壁面间碰撞次数,且颗粒-壁面间碰撞次数的波动幅度随Uf增加而增加。

图8 不同Uf下碰撞次数Fig.8 Number of collisions at different Uf values

图9 不同Uf下各物理量分布Fig.9 Distribution of various physical quantities at different Uf values

图9为不同Uf下各物理量分布。由图9a可见，随着Uf增加，管壁附近颗粒体积分数逐渐降低，颗粒体积分数峰值逐渐向管道中心靠近且逐渐增加，当Uf为4.5 m/s时颗粒体积分数在中心存在最大值。这是因为随着Uf增加，壁面附近滑移速度增加，颗粒受到的升力变大，推动更多的颗粒向管道中心聚集。

由图9b、9c可见，随着Uf增加，流体和颗粒的轴向速度也随之增加，其径向梯度也随输送速度的增加变得更加陡峭，这在管壁附近更加明显。输送速度对流体均方根速度的影响如图9d所示。流体均方根速度随Uf增加而增加，且陡峭程度逐渐增加。这是因为随着Uf增大，管道内雷诺数增大，从而湍流强度总体变大。在壁面附近，流体均方根速度的增加更加明显。

图10 Uf对无量纲压降的影响Fig.10 Effect of Uf on dimensionless pressure drop

由图10可见，随着Uf增加，管内压降也快速增加，这与文献[29]中的实验一致。这是因为随着Uf增加，流体-壁面的摩擦损耗加剧，以及更加频繁和剧烈的粒子间相互作用，导致更多的能量被耗散，从而造成更大的压降。

3.4 颗粒直径对水力输送的影响

在排污设施中，输送的固体污物往往大小各异。因此，研究不同粒径的颗粒对两相流动的影响具有重要的意义。对相对较粗的颗粒(dp为1.0、1.5、2.0 mm)进行数值计算，研究颗粒直径对固-液两相流动的影响。

不同dp下颗粒的分布情况如图11所示。从图11 看出，当dp较小(dp为1.0 mm)时，颗粒在管道内分布非常均匀，并且壁面附近也存在一定量的颗粒。当增大颗粒直径到1.5 mm时，颗粒较为集中，管道壁面附近的颗粒减少。当dp继续增加到2.0 mm时，颗粒分布极不均匀，大部分颗粒集中在管道中心，壁面附近几乎没有颗粒。因此，从颗粒的分布可以看出，小颗粒在管道内分布均匀，而大颗粒在管道中心聚集，且壁面附近几乎不存在颗粒，表现出极不均匀性。

图11 不同dp下颗粒分布Fig.11 Distribution of particles at different dp values

在不同dp下，碰撞次数随时间的变化如图12所示。从图12a看出，颗粒间碰撞次数随dp增加而减少。随着dp增加，相同浓度下颗粒数量的急剧减少会降低颗粒间碰撞次数，而大颗粒在管道中心更加聚集会增加中心处颗粒间碰撞次数，由此可见，颗粒数量的急剧减少对颗粒间碰撞次数的影响大于颗粒向中心聚集带来的影响。从图12b中看出，在颗粒直径较小时，颗粒分布较为分散，颗粒-壁面间碰撞次数比颗粒间碰撞次数少，但未到可忽略的程度。在颗粒直径较大时，颗粒主要集中在管道中心处，颗粒-壁面间碰撞次数远小于颗粒间碰撞次数。

图12 不同dp下碰撞次数Fig.12 Number of collisions at different dp values

图13为不同dp下各物理量的分布。由图13a可知,同浓度下，颗粒直径越小，在管道内的分布越分散。这是因为小颗粒在管道内的颗粒数量较大，管道中心处颗粒间的接触次数增加，使得颗粒更分散，导致近壁面处存在一定颗粒。当输送速度不变时，流体轴向速度径向梯度基本一致，但是大颗粒与流体滑移速度更大，颗粒受到的升力就更大，管道壁面处颗粒向中心聚集，这也能解释大颗粒更加向管道聚集的现象。

由图13b、13c可知，流体轴向速度随着dp增加呈现略微降低的趋势，这种趋势在管道中心较为明显。同时发现，随着dp增加，颗粒速度也存在一定量的降低。这是因为随着dp增加，流体对大颗粒的夹带作用变小，即大颗粒的随动性较低，因此颗粒速度较低。由图13d可知,dp对流体均方根速度的影响很小。

由图14可知，随着dp增加，压降的变化不大，这在文献[30]的实验中也有相同的报道。因此，输送速度对压降的影响大于进口颗粒体积分数，而颗粒直径的影响最小。