李红文， 张 涛

（天津大学 电气与自动化工程学院，天津 300072）

0 引 言

在以能源动力及化工为代表的许多工业领域，气固两相流有着非常广泛的应用，其中稀相气固两相流的管道流动中，流量测量非常常见，如燃煤电站煤粉输送管道的在线测量、锅炉系统中含尘烟气的流量测量、设备中含有粉尘气体的管道测量等等。

文献［1］采用在煤粉吹送管道中加装标准文丘里管等节流元件，通过差压法测量煤粉浓度；文献［2－3］采用文丘里管测量煤粉气固两相流流量。由于节流件的作用，管道内部流体流动情况多变，流场复杂，在这种情况下，固相颗粒的存在对两相流的测量结果产生影响。通常节流件往往用于单相流体测量，在用于两相流测量时，需要进行优化设计，此时固相颗粒的影响是需要深入研究的问题，而两相流测量中所有的问题都通过实流实验来研究，显然不方便。

随着计算机软硬件及仿真应用技术的发展，计算流体力学（computational fluid dynamics，简称CFD）仿真模拟技术在流体力学领域中得到很大发展。它突破了实物条件的限制，使用较方便，应用较广泛。其中Fluent就是最具代表性的仿真软件之一，已经在工程实践中得到较广泛的应用［4－5］，效果较为满意。Fluent中开发了专门的离散相模型（discrete phase model，简称DPM），是研究气固两相流离散相最有力的工具。它能针对上面提到的这些问题进行离线仿真模拟研究，但是在文丘里管的仿真模拟中发现，其计算的准确性与上述工况的实流实验及文献中的数据结果存在一定的偏差，结果不尽如人意。

本文针对两相流中颗粒受力情况，结合气相流场的特点，对离散相颗粒行为进行了细化分析，对颗粒所受的各个力进行了合理的取舍，借助用户自定义函数（user－defined function，简称 UDF）这一有力工具，对文丘里管内低浓度气固两相流进行DPM模型算法的改进与优化，达到了提高仿真模拟准确性的目的。

1 建模仿真与算法原理

1.1 Gambit网格模型与气相仿真流场描述

在Fluent中，DPM是求解气固两相流问题的一种常用方法，其中气体为连续相，固体颗粒为离散相。气相流场的准确计算与模拟，是准确研究离散相颗粒行为的先决条件。依据试验用管道实体情况建立Gambit网格模型，文丘里管及其前后直管段尺寸采用文献［6］中数据。

采用Gambit对文丘里管模型进行网格划分时，在喉部适当加密，既保证网格质量，又不至于使网格数目过多。

仿真计算中气相控制方程采用标准k－ε黏性模型，采用三维单精度求解器，残差设为0.000 1，离散格式为一阶迎风格式，管道入口为速度入口模式，出口为自由流出口。为提高仿真的准确性，采用了网格自适应adapt功能［7］。

气相仿真中文丘里管附近轴向剖面速度云图，如图1所示。

图1 文丘里管附近轴向剖面仿真速度云图

从图1中可见，由于具有前置直管段，文丘里管前1D处轴向剖面速度场形成规则的管内湍流速度分布，文丘里管后2D范围以外流场又开始逐渐恢复成流入文丘里管之前的湍流状态；而在文丘里管喉部前后，形成高速的速度核心区，此处流场复杂，速度与压力梯度大。关于文丘里管内节流流场的详细分析见文献［3］。

1.2 固体颗粒在气相流场中的受力分析

仿真中离散相固体颗粒定义为煤粉颗粒，形状为球形，离散相颗粒在气相流场中除受到重力与浮力外，在运动过程中受到多种力的作用，不同力对颗粒运动所起的作用不同［8－9］。

颗粒处在运动状态时具有的惯性力为：

其中，up为颗粒的速度；ρp为颗粒的密度；dp为颗粒直径。

气相流场中，颗粒对气体具有阻力，其反作用力是气体对颗粒的曳力，表达式为：

其中，u为气相流体（空气）的速度；ρ为空气的密度；rp为颗粒半径；CD为曳力系数或阻力系数。工程应用中，（2）式可有不同的表达形式。

压力梯度力是气相流场中压力梯度对颗粒引起的作用力，表达式为：

其中，Vp为颗粒体积。

由表观质量效应产生虚假质量力，表达式为：

Magnus升力是由于颗粒的旋转而产生的升力，一般与重力有相同的数量级，表达式为：

其中，ω为颗粒自转角速度。

Saffman升力是由于颗粒在具有速度梯度的流场中运动，颗粒上下侧流体的速度差对颗粒产生的升力，表达式为：

其中，μ为空气的动力黏度。

Basset力属于瞬时流动阻力，在黏性流体中，它与流动的不稳定性有关，其大小与颗粒的运动经历有直接关系，其计算公式为：

1.3 通用DPM中仿真计算模型的算法

模型中假定离散相很稀疏，其体积分数小于10%，这样可以忽略颗粒之间的相互作用。模型具体解法如下：对连续相流体在欧拉框架下求解纳维斯托克斯方程，当计算结果收敛到适当程度后（通常残差设为0.01），将颗粒注入流场，以单个颗粒为对象在拉格朗日框架下求解颗粒轨道方程，同时颗粒的质量载荷会对气相流场产生影响，当计算再次收敛时（残差设为0.000 1），得到计算结果，结果通常以颗粒轨迹来表示。

Fluent中通过对颗粒作用力微分方程进行积分，来求解离散相颗粒的轨道。根据牛顿运动定律，其在直角坐标系下的形式（以x方向为例）为：

其中，m为单颗粒质量；Fx为x方向上的其他力，在1.2中均已提到；Fd为颗粒所受曳力，其计算式为：

其中，Rep为颗粒雷诺数或相对雷诺数，定义为：

CD有5种经验公式模型可选，一般采用如下的简化表达式：

其中，α1、α2、α3、η为 DPM 通用模型库中存储的常数，由相关经验公式得出。

对（8）式积分就得到了颗粒轨道上每一个位置上的颗粒速度，而颗粒轨道通过对（12）式进行积分，即沿着每个坐标方向求解方程就得到了离散相颗粒的轨迹。

2 通用DPM模型的分析与优化

2.1 通用DPM模型分析与优化策略

Fluent软件中通用的计算模型一般都为理想模型或者基本模型，不完全适用于气固两相流文丘里管流量测量这一特定工况。气相流场在文丘里管前后很复杂，速度和压力梯度大，固相颗粒受其作用，在文丘里管及其前后区域，颗粒与管壁、文丘里管的内壁发生多次撞击，同时颗粒自身还有快速自转的行为；而当气相速度低于一定临界值时，颗粒还会沉积在文丘里管前方的管壁下侧。因此，需要自定义其中的一些物理模型，以提高针对具体工况模拟的准确性。

2.1.1 颗粒所受各种力的量级与取舍

针对1.2中颗粒所受各种力，进行如下处理。通常，颗粒曳力最大最明显，压力梯度力与颗粒的曳力相比约小3个数量级，虚假质量力比颗粒的曳力小3～4个数量级，故后2种力可以忽略［9］。DPM通用模型面板中预留了颗粒的Saffman升力选项，通用模型中此力往往被忽略，但根据（6）式可知，Saffman升力在流动的主流区，在流体速度梯度小的情况下可以忽略，而在边界层中其作用很明显，因此要提高仿真模拟的准确性，在DPM模型面板中应选择激活此选项。

通用DPM模型中认为Magnus升力可以忽略，原因是将颗粒看作表面光滑的球形，所受外力矩为0。而事实上，不论其形状如何，颗粒在气相流场运动过程中始终有速度很高的自转行为，可达每秒1 000转的量级［9－10］。故根据（5）式，Magnus升力大小与重力具有相同数量级。而重力比颗粒的曳力小1个数量级，因在（8）式中已经考虑重力，所以Magnus升力不宜忽略。经计算其与颗粒重力大小接近、方向相反，两者大小差异在10%以内，故本文认为Magnus升力与重力是互相抵消的。

通用DPM模型忽略Basset力。针对本文所研究工况，气体密度与颗粒密度之比小于0.001，经计算，Basset小于重力3个数量级，故本文将其忽略。而对于两相密度接近的DPM工况，此力需要准确计算。

2.1.2 仿真模型气相速度入口的选择

从图1可见，设置模拟管道的气相速度入口时，入口截面上一般采用匀速的固定速度流入，而实际上装置实验段之前已经有足够长的直管段，在这段管道中气相已经过充分发展形成规则的管内湍流速度分布。通常仿真与实验中所设定入口速度仅是一个平均速度u0，于是当计算颗粒的轨道时，入口处颗粒运动状态的初值与实际工况有很大的差异，而随着迭代计算的进行，这个差异的影响会持续下去。根据（6）式，气相入口处速度流场的梯度对颗粒的Saffman升力也产生很明显的偏差，因此采用自定义的速度入口。

将文丘里管两相模型的气相入口处速度定义为：

其中，R为圆形管道半径；uin为管道入口处距轴线r处的速度；um为圆管道中心线上的速度；u0为已知的入口平均速度；λ＝um／u0。优化前后模型的气相速度入口形式，如图2所示。图2a是通用DPM模型速度入口形式，而图2b更符合真实情况，在优化模型中采用。当进行具体工况的仿真计算时，根据雷诺数所在的区间，（13）式中的n及λ值，按照表1选择。采用自定义速度入口uin，提高了仿真的准确性。

图2 优化前后模型速度入口形式

表1 充分发展的管内湍流流场分布系数

2.1.3 颗粒曳力模型的选择

通用DPM模型给出的5种曳力模型分别对应各种常见工况。在一定的颗粒雷诺数工况范围内，通常CD可以表示为Rep的一个全区域的经验公式函数。对于本文工况，在管内流动充分发展段，颗粒受气流作用，雷诺数变动范围较大；而在文丘里管前后，气体湍动剧烈，雷诺数的变动范围要超出充分发展段的数倍。所以，全区域的曳力系数模型具有一定误差，选择更加准确的曳力系数模型是模型优化手段之一。

在（9）式中，曳力系数CD采用文献［9］提供的“球形颗粒曳力系数关联式”，如（14）～（23）式，其颗粒雷诺数适用范围涵盖了从0到无穷大。在DPM仿真模型中，将（14）～（23）式代入曳力系数数据库中，以提高仿真模拟的准确性。

（14）～（23）式中，H＝lgRe，Re为颗粒雷诺数Rep。

2.1.4 壁面碰撞模型的偏差与改进

固相颗粒会与管道壁面及文丘里管内壁发生多次碰撞，颗粒的物理性质和运动状态、速度、入射角度等都影响颗粒与壁面的碰撞规律，目前还没有反映颗粒与壁面碰撞规律的通用模型。在通用DPM模型中，只能采用reflect（反射）、trap（捕获）、escape（逃逸）3种边界条件，而一般认为颗粒与壁面之间是弹性碰撞，不存在能量损失，或者简单地设定一定的常数值恢复系数。这2种方法都简化了颗粒与壁面碰撞的复杂过程，给模拟结果带来一些误差。因此，寻求合理的颗粒壁面碰撞模型，能较为准确地反映颗粒碰壁后的行为，同时可以描述更为准确的颗粒轨迹，对提高仿真准确性具有意义。

文献［10－12］对颗粒壁面碰撞模型进行了PIV及激光全息实验研究，得到一些经验公式，但每个经验公式都有其特定的工作条件限制，其观点都是认为恢复系数与入射角及入射速率有关。文献［11］所提供流体及颗粒的理化特性和工况与本文非常接近，本文中引用其经验公式如下：

其中，e为颗粒碰壁恢复系数；v1为碰撞前颗粒的速率；θ1、θ2分别为颗粒的入射角与出射角。

2.2 模型优化策略的实现

Fluent软件提供了用户自定义函数（UDF）功能，用户可以编写Fluent源代码来满足各自的特殊需要。UDF中可使用标准C语言的库函数，并使用Fluent提供的预定义宏，从Fluent求解器中获得数据。

本文使用编译型 UDF，针对（14）～（23）式及其他有关公式进行C语言编程，然后通过Fluent的几个自定义宏来实现通用DPM模型的优化，具体如下：使用DEFINE－PROFILE（自定义边界截面上的变量分布）来定义空气入口的湍流速度剖面；使用DEFINE－DPM－DRAG（自定义流体中颗粒的曳力系数）来定义颗粒曳力模型的表达式；使用DEFINE－DPM－BC（自定义颗粒到达边界后的状态）来实现自定义的颗粒碰壁规律。将包含上述宏语句的C语言代码，经反复调试运行成功后，添加到Fluent通用模型中，再进行DPM仿真模拟，这样，通用DPM模型适用范围得以拓宽，并提高了仿真精度，解决了对通用模型进行优化的实际需求。

3 DPM模型优化前后仿真对比

3.1 优化后DPM模型的验证

为验证优化后模型的正确性，针对文献［6］中的文丘里管进行工况对比仿真，管型几何参数如图3所示，管径D＝52.5mm，喉部管径d＝26.5mm，收缩段长度L1＝56.7mm，喉部长度L2＝12.2mm，渐扩段长度L3＝150.2mm，文丘里管的孔径比β＝d／D。颗粒的密度为1 340kg／m3，平均粒径为43μm，形状为球形，在管道入口通过面射流源注入流场。

仿真中为求准确，激活固相颗粒与气相相互作用选项，DPM模型优化前后仿真均采用ANSYS Fluent13版本软件。

图3 文丘里管示意图

计算结果如图4所示，纵坐标差压比是两相混合物、气体通过文丘里管时入口与喉部之间差压的比值，即ΔP1m／ΔP1g，这是文丘里管测量时的关键参数。

由图4可见，优化模型的模拟值与文献［6］的实验结果吻合得更好，验证了优化后模型的准确性更具优势。曲线趋势是，在仿真计算范围内，固气质量比越高，优化后模型仿真结果比优化前更接近文献［6］的实验数据。

图4 仿真优化前后与实验对比曲线

3.2 对标准文丘里管流出系数的验证

文丘里管属于标准节流装置，根据文献［13］的规定，对于机械加工收缩段经典文丘里管，其流出系数为0.995。

取β＝0.6，固气质量比为1.0，雷诺数Re＝（2.0～10.0）×105，采用2种模型仿真计算流出系数时，结果见表2所列。

由表2可见，优化后DPM模型仿真的准确性有明显提高。

取Re＝5.0×105，固气质量比为0.2～1.8，2种模型仿真计算流出系数对比见表3所列。

由表3可见，在固气质量比加大的情况下，仿真结果更接近于国家标准规定值，固气质量比越大，优化后模型的准确性越明显。

在CFD实践中，计算误差在2.0%之内，往往认为是合理的，仿真计算结果具有参考价值，可见优化后模型符合此指标。

表2 Re变化时2种模型流出系数对比

表3 固气质量比变化时2种模型流出系数对比

4 结束语

本文通过对通用DPM模型进行优化，得到的理论仿真结果与国家标准及文献［6］结果一致。仿真计算结果表明，对文丘里管气固两相流差压比值及流出系数进行模拟时，优化后模型仿真的准确性明显高于通用DPM模型。优化后模型有助于更准确分析两相流测量工况，给文丘里管的选型设计研究提供参考，以提高测量精度。研究过程中，将理论及经验公式通过编程并与宏命令灵活结合形成UDF模块，再嵌入Fluent仿真程序中的方法，是理论模型得以实现的关键所在。同时本文提出的方法具有通用性，除了可用于节流管道，也适用于弯头、变径、分支管、阀门等其他具有复杂气相流场特点的管道内气固两相流动模拟，以及流化床、反应釜等化学反应器中气固两相流研究之中，这一方法在工程上具有实用性。

［1］ 林宗虎.能源动力中多相流热物理基础理论与技术研究［M］.北京：中国电力出版社，2010：573－579.

［2］ 艾卫国，余战英，徐通模，等.煤粉浓度测量阻件性能比较的实验研究［J］.动力工程，2000，20（5）：892－895.

［3］ 吴占松，谢 菲.用于管道煤粉流量测量的文丘里管型设计及优化［J］.清华大学学报：自然科学版，2007，47（5）：666－669.

［4］ 杨汶雨，张 喆，袁景淇，等.气升式反应器结构参数优化设计CFD仿真研究［J］.合肥工业大学学报：自然科学版，2009，32（9）：1301－1303.

［5］ 尹梦晨，谷正气，容江磊.基于CFD分析的汽车空气动力学制动研究［J］.合肥工业大学学报：自然科学版，2011，34（3）：326－330.

［6］ Lee J，Crowe C T.Scaling law for metering the flow of gasparticle suspensions through venturis［J］.ASME Journal of Fluids Engineering，1982，104（1）：88－91.

［7］ 李鹏飞，徐敏义.精通CFD工程仿真与案例实战［M］.北京：人民邮电出版社，2011：57－66.

［8］ 车得福，李会雄.多相流及其应用［M］.西安：西安交通大学出版社，2007：80－85.

［9］ 岑可法.工程气固多相流动的理论和计算［M］.杭州：浙江大学出版社，1990：362－373.

［10］ 施学贵，徐旭常.颗粒在湍流气流中运动的受力分析［J］.工程热物理学报，1989，10（3）：320－325.

［11］ 饶 江，葛满初，徐建中，等.固体颗粒与通道壁面相互作用的 实验研 究 ［J］.工程 热物理 学 报，2003，24（1）：134－136.

［12］ Di Maio F P，Di Renzo A.Analytical solution for the problem of frictional－elastic collisions of spherical particles using the linear model［J］.Chemical Engineering Science，2004，59：3461－3475.

［13］ GB／T 2624.4－2006，用安装在圆形截面管道中的差压装置测量满管流体流量：第4部分 文丘里管［S］.