不同液固动量交换系数模型对颗粒流化行为数值模拟的适应性

孙铭阳，于传兵，吕东，韦鲁滨

(1.中国恩菲工程技术有限公司，北京，100038；2.中国矿业大学(北京) 化学与环境工程学院，北京，100083)

由于其优异的混合、传质和传热特性，液固流态化技术已被广泛用于冶金、选矿、化工和环保行业[1-4]。国内外学者针对液固流化床进行了大量研究，以揭示流化床内颗粒运动特点和床层膨胀等行为。胡新辉等[5]对液固流化床内颗粒循环运动机理进行了研究，提出了液固流化床内流体和颗粒径向速度分布的数学模型，并认为颗粒循环运动由流体径向速度分布不均导致；DOROODCHI等[6]借助CFD研究了倾斜板对颗粒在液固分选流化床内流化行为的影响，并比较了CFD模型和理论模型的优劣；张锴等[7-8]改进了双流体模型，从而较全面地考虑了离散颗粒对液固两相动量的影响，并用该改进的模型对入口速度突变时液固流化床内固体颗粒时空分布特性进行了详细研究；沈志恒等[9]借助CFD对倒置液固流化床内颗粒速度、浓度和床层孔隙率等分布特性进行了模拟，发现倒置业液固流化床内存在中间向下两侧向上的颗粒运动规律；张卫义等[10]对液固流化床内颗粒最小流化速度进行了研究，提出了一个Ergun型最小流化速度模型，并与其他模型进行了对比，在其所列举的试验数据下，该模型的预测结果比Wen-Yu模型的优。 液固动量交换系数是影响颗粒流化过程数值模拟准确性的重要因素。目前针对液固动量交换系数对液固流化床数值模拟准确性影响方面的研究较少。本文作者对不同粒度或密度的9种颗粒流化过程分别进行试验研究，并结合试验结果对比分析不同液固动量交换系数对颗粒流化行为数值模拟的准确性，本文试验方法和数值模拟结论可为相关液固流化床数值模拟研究提供参考。

1 颗粒流化试验与结果

1.1 流化颗粒与流化试验装置

流化颗粒分别选择不同粒度、密度的石英砂和煤颗粒，共包括如表1所示的9种颗粒。

搭建了如图1所示的颗粒流化试验系统。上升水流由平行管流体分布器射出后，沿柱体向上流动并在径向上逐渐均匀分布，经第1层孔径0.045 mm筛网后，流体轴向速度基本实现均匀分布。流化颗粒置于第2层孔径0.045 mm筛网上，筛网孔径远小于流化颗粒直径，可以保证底层颗粒受到均匀向上的流体曳力，且不会沿近壁面区漏下。流化床溢流给到过滤箱，过滤箱内亦布置1块孔径0.045 mm筛网，溢流经筛网过滤后返回至清水池，从而实现带出颗粒的回收和流化水的循环利用。

表1 流化颗粒物性Table 1 Properties of fluidized particles

图1 液固流化试验系统Fig.1 Experimental system for liquid-solid fluidization

1.2 流化试验结果

令H为上升水流速度为uF时颗粒床层高度，则此时液固流化床层内颗粒平均体积分数可表示为

床层孔隙率为

床层有效密度则可按下式计算：

式中：ρl和ρs分别为液相和颗粒相密度。

图2所示为不同颗粒所形成液固流化床层有效密度随上升水流变化情况。从图2可以看出：颗粒粒度越小，密度越小，所形成的液固流化床层密度对上升水流越敏感。当2种颗粒密度相差较大时，在相同流化速度下，不同颗粒形成的床层密度流化曲线基本不能相交。通过回归分析发现：床层密度ρb与上升水流速度u0呈现形式的自然指数函数关系，同时也可被形式的幂函数较好拟合。

图2 液固流化床层密度与上升水流关系Fig.2 Bed effective density at various fluidized velocities

2 颗粒流化过程数值模拟方法

2.1 液固两相流动控制方程

采用 Euler-Euler法模拟液固流化床层内固液两相流动。

液、固两相的质量守恒均满足下式：

式中：αp，up和ρp分别为p相体积分数、速度和密度；p和k分别表示固相或液相，当p为固相时则k为液相，反之亦然；mpk和mkp分别为由p相到k相的转移质量和由k相到p相的转移质量；Sp为p相的质量源相项。对于惰性颗粒流化过程，mpk，mkp和Sp可忽略，且液固流化床内流动可看作不可压缩，因此，式(4)可进一步简化为

液相与各颗粒相存在相互作用，液相动量方程为

式中：αl和ul分别为液相体积分数和速度；us为颗粒相的速度；τl为液相应力张量；g为重力加速度；Ksl(ul-us)为液相与颗粒相的交换动量，Ksl为液相与固相间的动量交换系数；最后一项为液相因施加给离散固体颗粒相升力、压力梯度力和湍流分散力等而受到的相应反作用。

颗粒相s的动量方程为

式中：αs为颗粒相体积分数；τs为颗粒相应力张量；Flift,s为颗粒相所受升力；Fvm,s为颗粒相所受附加质量力；Ftd,s为湍流中颗粒相所受湍流分散力；Kls(ul-us)为颗粒相所受流体阻力。

式(7)中τs由应力张量可求颗粒相微元发生形变时所受法向和切向应力。通过颗粒相应力张量得到由速度和压力表示的颗粒相运动方程，进而通过求解该方程得到颗粒相流场。颗粒相应力张量具有以下形式：

式中：μs为颗粒相剪切黏度；I为单位张量；Ps为颗粒相压力；λs为颗粒相体积黏度；这些变量为颗粒脉动能或法向恢复系数的函数。

类比气体分子动理论，可以得到颗粒温度的概念。颗粒温度实际表征了颗粒相中固体颗粒随机脉动速度的强度，因此又称颗粒脉动能，其表达式如下：

式中：us,i为颗粒i方向的脉动速度。

根据颗粒动理学理论，颗粒脉动能的输运方程为

2.2 计算域离散和求解方法

数值模拟中液固流化床几何结构如图3所示，将液固流化床层计算域划分为结构网格。速度入口给入均匀稳定的上升水流，上升水流速度根据试验值进行设定。数值模拟选择二维双精度瞬态求解器。虽然表观上升水流速度在0.1 m/s以下，但流化床内的流动是工程上典型湍流流动，选择RNGk-ε湍流模型来预测流化床内液固两相运动[11-12]。采用有限体积法(FVM)对控制方程进行离散；选择 QUICK插值方法，将控制单元节点处各物理量值表示成相邻控制单元面上值；选择 Phase coupled SIMPLE算法来进行压力与速度的耦合求解。当液固两相体积分数残差小于10-4，其余物理量的残差小于10-3时，认为该时间步内的迭代达到收敛[13]。

图3 数值模拟中液固流化床几何模型Fig.3 Geometry and mesh of fluidized bed

3 不同液固动量交换系数模型模拟结果

3.1 液固动量交换系数模型

式(6)和式(7)中Kls=Ksl，为颗粒相与液相动量交换系数。

式中：τs为颗粒弛豫时间；f为与颗粒相体积分数和阻力系数等相关的函数。

由式(11)可知：动量交换系数是曳力系数、颗粒相体积分数和液固相对速度等物理量的函数。目前，液固两相动量交换系数的数学模型主要有Syamlal-O’Brien[14]、Wen-Yu[15]、Gidaspow[16]、Gibilaro[17]和Huilin-Gidaspow[18]。

Syamlal-O’Brien模型如下：

阻力系数CD采用DALLA VALLE表达式[19]求得：

Res为颗粒相雷诺数，有

式中：ds为s相颗粒直径。

式(12)中ur,s为s相颗粒沉降末速修正，可按下式求得：

Wen-Yu模型数学形式如下：

Wen-Yu模型中曳力系数表达式为

Gidaspow模型为Wen-Yu模型和Ergun方程的组合，当液相体积分数α1＞0.8时，表达式为Wen-Yu模型；当液相体积分数α1≤0.8时，表达式为Ergun方程，具体形式如下：

Gibilaro模型具有以下形式：

Huilin-Gidaspow模型也是Wen-Yu模型和Ergun方程的组合，当床层空隙率的变化跨越临界点(0.8)时，利用一个平滑过渡函数来实现Wen-Yu模型和Ergun方程间的转换，其形式如下：

式中：Ψ为平滑过渡函数，Ψ=0.5+arctant(262.5×(αs-0.2))/π。

3.2 不同液固动量交换系数模拟床层高度对比

利用以上 5个液固动量交换系数模型模拟0.25～0.35 mm石英砂在12.3 mm/s上升水流速度下流化过程。床层运动稳定后，对比5个模型模拟所得床层内固体颗粒体积分数分布如图4所示。从图4可知：预测床层体积分数最大的为Gibilaro模型，最小的为Syamlal-O’Brien模型；Gidaspow 和 Huilin-Gidaspow模型预测的床层内颗粒体积分数分布完全一致；Wen-Yu模型与Gidaspow 和Huilin-Gidaspow模型相比，预测的颗粒体积分数略高。

图4 12.3 mm/s流化速度下各模型预测0.25～0.35 mm石英砂颗粒体积分数分布Fig.4 Distribution of 0.25-0.35 mm quartz particles at 12.3 mm/s fluidization rate

床层高度模拟值确定方法如图5所示，用不同液固动量交换系数模拟所得不同上升水流速度下床层高度与试验值对比如图6所示。从图6可知：Wen-Yu，Gidaspow 和Huilin-Gidaspow对床层高度的预测值与试验值最吻合，且Gidaspow 和Huilin-Gidaspow床层高度预测值完全相同；Wen-Yu，Gidaspow，Syamlal-O’Brien，Gibilaro和Huilin-Gidaspow这5个模型所得床层高度预测值与试验值间的均方根误差分别为5.67，5.01，47.14，24.42和5.01，表明Gidaspow和Huilin-Gidaspow模型准确度略高于Wen-Yu模型的准确度。

图5 预测的床层高度h确定方法示意图Fig.5 Method of determining height of predicted fluidized bed

图6 不同液固动量交换系数模型模拟结果对比Fig.6 Predicted bed height by using different fluid-solid exchange coefficients

图7所示为图5中对应的石英砂颗粒相速度矢量图。从图7可知：Gibilaro和Gidaspow模型得到的石英砂相只有在液固界面附近存在较大的速度变化，床层的其余位置处石英砂颗粒的速度很小，可认为准静止地悬浮于上升水流中；Wen-Yu，Syamlal-O’Brien和 Huilin-Gidaspow模型得到的石英砂相存在由中心向上再流向壁面、沿壁面下降的循环流动，这与试验中观察到的床层流动状态吻合。因此，从床层高度预测值和床层内液固两相流态 2方面来讲，Huilin-Gidaspow模型比其他4个模型更合适。

图7 石英砂颗粒相速度矢量分布Fig.7 Velocity vector distribution of quartz particles phase

图8 0.25～0.35 mm 石英砂在12.3 mm/s流化速度下流化过程Fig.8 Fluidization process of 0.25-0.35 mm quartz particles at 12.3 mm/s fluidization rate

图8所示为0.25～0.35 mm石英砂在12.3 mm/s上升水流速度下流化过程。从图8可知：随着流动时间的增长，床层膨胀度逐渐增大；40 s时床层高度基本达到稳定，但此时床层内液固两相流动还未达到平衡，由石英砂相的流动速度可以看出，0～100 s之间颗粒相运动比较混乱，没有明显规律；到达120 s后，床层内液固两相流动基本实现稳定，颗粒相表现出由中心向上，到达床层顶端后，向两侧壁面流动，最后沿壁面低速区下降的循环流动趋势，与试验中观察到的现象吻合。

图9 4种颗粒流化床层颗粒体积分数模拟值与试验值Fig.9 Comparison between simulated and experimental volume fraction of particles phase

进一步选择Hui-Gidaspow模型模拟0.15～0.20和0.50～0.63 mm石英砂颗粒以及 0.25～0.35和 0.90～1.00 mm精煤颗粒在不同上升水流速度下流化过程。图9 所示为4种颗粒在不同上升水流速度下床层颗粒体积分数试验值和预测值。从图9可知：床层颗粒体积分数预测相对误差基本都在 5%以内，说明在用Euler-Euler法模拟颗粒流化过程中，选择 Hui-Gidaspow模型得到的模拟结果具有较高的准确度。

4 结论

1) 5种液固动量交换系数模型对粒度为 0.25～0.35 mm的石英砂颗粒流化床层高度预测准确度由高到低依此为：Gidaspow模型和Huilin-Gidaspow模型，Wen-Yu模型，Gibilaro模型，Syamlal-O’Brien模型，其中Huilin- Gidaspow和Gidaspow模型均方根误差都为5.01，Syamlal-O’Brien模型均方根误差为47.14。

2) 利用Huilin-Gidaspow模型模拟得到的液固流化床层内颗粒相存在由中心向上，到达床层顶端后，再向两侧壁面流动，最后沿壁面低速区下降的循环流动过程，与试验中观察到的现象相符，而 Gidaspow模型模拟结果中颗粒相不存在循环流动。从床层高度预测值和颗粒相运动状态 2方面考虑，5个模型中Huilin-Gidaspow模型最适合液固流化床内颗粒流化行为数值模拟。

3) 选择 Huilin-Gidaspow模型对 0.15～0.20和0.50～0.63 mm石英砂颗粒以及 0.25～0.35和 0.90～1.00 mm精煤颗粒在不同上升水流速度下流化过程进行了模拟，床层内颗粒相体积分数预测值相对误差基本都处于5%以内，说明利用Euler-Euler法，同时选择 Hui-Gidaspow液固动量交换系数模型模拟颗粒流化过程具有较高准确度。