气固两相流穿越液池过程颗粒运动及分布特性

吴晅，梁盼龙，王丽芳，李海广，武文斐

（1内蒙古科技大学内蒙古自治区白云鄂博矿多金属资源综合利用重点实验室，内蒙古 包头014010；2内蒙古科技大学能源与环境学院，内蒙古 包头 014010）

引言

在激冷式水煤浆气化炉洗涤冷却室内，存在携带固态灰渣的合成气穿越液池以完成洗涤的工艺过程[1]。不仅如此，在浸没燃烧装置中[2-3]，在冲击水浴除尘器[4]中都广泛存在气固两相流穿越液池过程。与其他干法或湿法洗涤方式不同，该过程涉及复杂的气液固三相复杂湍流流动。同时，冲击式洗涤过程、泡沫式洗涤过程和淋浴式洗涤过程3种洗涤过程[5]贯穿其中，共同构成了气固两相流穿越液池的复杂气固分离过程。气固分离过程与颗粒的运动行为密切相关，而颗粒的分布特性又是表征颗粒运动行为的方式之一。为此，了解颗粒运动及其分布特性对提高气体的洗涤除渣效果以及完善设备结构都具有重要的指导意义。在实验方面，付碧华等[6]对激冷式水煤浆气化炉洗涤冷却室内的气液两相流动特性进行了研究，但其忽略了固体颗粒的存在。为此，贺必云[7]搭建了三相流动冷态实验平台，实现了洗涤冷却室内固相颗粒流动的实验和数值模拟研究并获得了颗粒的浓度分布特性，但其没有对不同粒径颗粒的运动进行对比，而且其在数值模拟方面还无法实现对颗粒微观运动过程的直观描述。在气液固三相数值模拟算法方面， Xu等[8]建立了CFD- VOF-DPM三相流动模型，在不同操作条件、液体性质和管口形状等因素下仅对单一气泡在液固混合液中形成、上升的运动过程进行了模拟，但该模型忽略了多气泡间相互聚并、破碎等作用的影响，从而与实际情况存在较大差距。Baltussen等[9]则采用直接模拟的方法，对气液固三相浆态反应器中的气液固三相流动进行数值模拟，但其模型将气泡和颗粒均视为离散相，并逐一跟踪，因此该算法计算量较大，只能针对小尺度物理模型进行计算。目前对气固两相流穿越液池过程中的三相流动数值模拟还少有相关文献可以借鉴。

为此，鉴于前人研究中存在的不足和现状，本文将采用实验和数值模拟的方法，分别对气固两相流在穿越液池过程中的不同粒径固体颗粒运动过程及其分布规律进行研究。以期加深对该洗涤净化工艺及其中气液固三相流动分布特征的认识。

1 实验平台的建立

1.1 实验装置

图1 实验系统Fig.1 Schematic diagram of experiment

实验系统如图1所示。液池内径480 mm、高1700 mm，下降管内径80 mm、外径为92 mm，下降管长度为1200 mm。为便于摄像，实验台主体采用钢化玻璃制作。实验台其他部分采用有机玻璃制作。在液池壁面上布置多个取样孔进行取样分析。实验中的气液固三相系统采用空气-纯净水-实心玻璃微珠三相系统。

实验系统主体主要由下降管和液池组成，液池内装有液态纯净水，下降管出口顶部浸没在水中。气体由空气压缩机送入储气罐，经稳压后通过流量计调节流量，随后经三通分流阀分成主、副两个通道。主通道为主气体通道，气体通过主通道进入气固混合室。同时，副通道中的气体进入颗粒发生器，颗粒发生器所产生的气固两相流进入混合室与来自主通道的气体混合。混合后的气固两相流进入下降管，经下降管出口进入液池，从而在气固两相流穿越液池的过程中，实现气固分离。

1.2 实验数据处理方法

实验主要是获得液池内的液体中所含固体颗粒的局部平均质量浓度。当整个三相流动过程处于准平衡状态时，利用针筒式取样管，通过取样孔抽取液池中的液固混合溶液。选取不同长度的取样针管，以确保获得不同径向位置处的液固混合溶液样品。取样针管外径非常小，其对液池内三相流动的影响可以忽略。采用精度为百万分之一的分析天平获得空管取样管的质量（mtube）及取样后取样管和样品悬浮液的总质量mtotal。读取取样管刻度获取悬浮液体积信息，该悬浮液体积由水的体积和悬浮液中的颗粒体积共同组成。得出与该悬浮液相同体积的水的质量（ml），再根据水与颗粒的密度由式(1)计算得出所取悬浮液样品中颗粒的质量（mp）。依据式(2)由悬浮液样品中颗粒的质量（mp）除以悬浮液体积（V）即可获得单位容积液池液体中所含颗粒质量浓度（C）。

式中，ρp为颗粒密度，ρl为液体密度。依据式(1)、式(2)，通过统计每次进入液池的颗粒总质量，获得输入单位质量颗粒总量下的液池液体中所含局部平均颗粒浓度。为保证测量精度，在同一高度取样时，沿液池径向选取多个取样点取样并求算术平均，从而获得液池内该高度处的局部平均颗粒浓度。

2 数值模拟平台建立

在准确判断颗粒在气体和液体中受力及颗粒间相互碰撞时，需要界定气液界面。同时Lagrange框架下的确定性颗粒轨道模型能形象、直观地描述出不同粒径离散颗粒在液池中沉降、悬浮和穿越等一系列的细微运动过程。为此，在本文数值模拟平台的建立上，采用 Euler框架下的气液两相流模型对气液流动进行模拟，并引入气液界面跟踪模型实现对气液界面的清晰界定；采用Lagrange框架下的确定性颗粒轨道模型对离散颗粒相的受力和运动进行跟踪计算。

模型简化和假设如下：视气体中的固体颗粒为球形颗粒；在气固两相流穿越液池过程中，各相间已达到热质平衡，无传热传质现象；忽略离散颗粒相对气液连续相的影响[10-11]。

2.1 气液连续相数学模型

气液两相连续方程为

动量方程为

式中，ui为i方向上的速度；p为压力；μ为黏度；为Reynolds应力项；Fv为体积力源项。

湍流模型采用 RNG k-ε湍流方程。湍动能 k方程为

湍动能耗散ε方程为

式中，Gk为湍动能产生项； C1ε= 1 .42，

引入 VOF气液界面追踪模型来描述气液两相间的运动界面，其流体体积函数方程为

采用连续表面张力模型[12]模拟相界面上的表面张力并以体积力的形式加入动量方程的源项中。

式中，σ为表面张力系数，κ为曲率。

2.2 离散相数学模型

离散相颗粒运动的计算采用确定性颗粒轨道模型。每个颗粒的运动过程可分解为在流体作用下运动的悬浮过程及在其他颗粒对其作用下运动的碰撞过程[13]。

单个颗粒在流体作用下的运动通过牛顿第二定律来控制，其运动方程如下

式中，up为颗粒的线速度，m·s-1；mp为颗粒质量，kg；F为颗粒受到的合力，N；ωp为颗粒的角速度，rad·s-1；Mp为颗粒所受的合力矩，N·m；Ip为颗粒的转动惯量，kg·m-2。

在经过一个时间步长Δt后的颗粒速度、角速度及位移可分别表示为

式中，上角标0表示上一个时间步长；Sp为颗粒的位移，m。

由于液池内水的黏性不大，故忽略Basset力对颗粒的影响。因此在本文中主要考虑颗粒受到的力有重力FG、曳力FD、Margnus力FM、浮力FB及附加质量力FA。

重力FG为

式中，dp为颗粒粒径，m；ρp为颗粒密度，kg·m-3。

单颗粒在流场中所受曳力FD由式(12)计算[14-15]

式中，rp为颗粒半径，m；ρc为气相或液相的密度，kg·m-3，计算时是选取气相还是液相的属性，则根据颗粒所处相来决定；uc为气液相的速度，m·s-1；up为颗粒线速度，m·s-1；CD为流体曳力系数，其是颗粒Reynolds数Rep的函数，可由式(13)计算得到[16]

其中

在液固悬浮液中，颗粒所受到的曳力与颗粒所处位置的液体容积份额有关，故在计算液体中颗粒所受曳力时，其曳力系数DC′可修正为[17]

式中，lα为液体容积份额。

Margnus力FM表达式为

颗粒所受的附加质量力FA可表示为

式中，ρl为液体密度，kg·m-3。

颗粒在液体中将受到液体对它的浮力FB作用

式中， Vp为颗粒体积，m3。

在液体中，颗粒碰撞前由于液膜的存在，颗粒的碰撞速度会有较大的衰减。因此颗粒在液体中应考虑液体剪切力对颗粒的作用[18]。这里采用近距离相互作用模型（close distance interaction model，CDIM）来考察液体剪切力对颗粒在液体中相互碰撞的影响。该模型认为由于颗粒间液膜的存在对颗粒的碰撞接触有十分重要的阻尼作用。为此需要确定颗粒碰撞前的垂直接触速度。其可描述为[17]

式中，h为两颗粒之间距离的中点到颗粒中心的距离；rp为颗粒半径；f、φ为修正系数，其计算如下

应用 Runge-Kutta法对式(18)进行计算可求得碰撞前的颗粒接触速度。

颗粒间的相互碰撞作用采用的是直接模拟Monte Carlo（DSMC）方法。在某一颗粒计算网格内，取样颗粒i和取样颗粒j所代表的真实颗粒发生碰撞的概率Pij与颗粒i以相对速度Gij在Δt时间内运动所掠过的体积和该体积内颗粒j的数浓度有关，可表达为[19]

式中，wj为取样颗粒j的数目权重；Vi为取样颗粒i所在网格的体积，m3；dp,i为颗粒i的直径，m；dp,j为颗粒j的直径，m；Gij为颗粒i和j的相对速度，m·s-1。

取样颗粒i和同一网格内其他所有颗粒的碰撞概率Pi为[19]

式中，N为取样颗粒i所在网格的取样颗粒总数。

根据式(20)，可运用修正的 Nanbu方法[20]（modified Nanbu method）确定与颗粒i发生碰撞的颗粒。即在产生[0,1]之间的均匀随机数R后，用式(21)选择与颗粒i可能发生碰撞的备选颗粒j

随后，利用随机数R(0＜R＜1)来判断颗粒i是否与颗粒j发生碰撞，如

满足，则颗粒i和颗粒 j在该时间段内将发生碰撞。颗粒i和颗粒j的线速度和角速度采用硬球模型计算得到，即

式中，f为摩擦系数；e为回复系数；ω为颗粒角速度。碰撞前的相对速度 G0=-；切向方向上单位矢量 t =/[]；相对速度切向分量

在数学模型的计算上，利用FLUENT软件为平台，在 Euler框架下对气液连续相进行模拟并逐一存储各时间段的气液连续相流场数据及网格数据。在Lagrange框架下运用VC++语言自编离散相颗粒运动程序，完成颗粒计算网格的划分、气液连续相数据的提取和离散相颗粒的受力及跟踪计算，从而获得每个颗粒的运动轨迹。在颗粒的受力计算上，依据所划分的颗粒计算网格中气相或液相所占容积份额的比重，来判定颗粒是处于气相中还是液相中，从而相应算出颗粒所受到的作用力。

3 实验与数值模拟结果及分析

3.1 气固两相流穿越液池的三相流动过程

图2 5～15 μm粒径颗粒的运动过程Fig.2 Movement of particles of 5—15 μm in diameter

图3 80～110 μm粒径颗粒运动过程Fig.3 Movement of particles of 80—110 μm in diameter

本文选取 5～15 μm、80～110 μm 和 380～530 μm 3种粒径范围的颗粒进行实验以观察不同粒径颗粒在随气体穿越液池过程中的运动规律。利用高清摄像机拍下固体颗粒随气体穿越液池的三相流动过程。并对3种情况下的颗粒运动过程进行对比，如图2～图4所示。图2为直径5～15 μm粒径颗粒的运动变化过程。从图中所示颗粒运动过程可以看出，颗粒随气体沿下降管进入液池后，一部分颗粒依靠惯性向下冲击进入液体并逐渐沉降到液池底部。此外，由于细小颗粒随气体运动的跟随性较好，大部分更为细小的颗粒被气泡所挟裹，反折向上运动。在向上穿越液池的过程中，部分颗粒被液体捕获，并在液面以下区域悬浮，从而该区域呈现较高浓度的颗粒分布。随后，由于液体的回流作用，在液面悬浮的颗粒逐渐向液池下部区域扩散，并最终充满整个液池空间。从图2中还可以看到，有一部分颗粒随气体穿越了液池，进入液池上部空间，因此在液池自由液面以上空间，可以看到明显的颗粒分布。

图3 和图4分别为80～110 μm和380～530 μm粒径颗粒随气体进入液池的运动变化过程。将图2～图4进行对比可以看出，随着颗粒粒径的增大，颗粒随气体的跟随性降低，而颗粒的惯性作用却增强，这使得颗粒更容易沉降到液池底部，从而较大粒径颗粒会在液池底部形成较高的浓度分布。此外，从图4所示的380～530 μm粒径颗粒运动过程可以看出，由于其粒径较大，惯性较强，在随气体一同冲击进入液池后，就直接进入液体中，并快速沉降到液池底部。同时从图4还可以看到，该粒径的颗粒几乎已不能随气体一起穿越液池并进入液面以上空间，这也说明该粒径颗粒基本都被液体所捕获。

3.2 液池内颗粒浓度分布的实验测量

图5 不同粒径颗粒沿液池高度方向浓度分布Fig.5 Concentration distribution of different size particles along height direction of liquid path

图5给出了气体流量为10 m3·h-1下，不同粒径颗粒沿液池高度方向上的浓度分布。本文以液池底部为零点，图中0.6 m处即为液池静态液面高度，0.36 m处即为下降管出口所处高度。从图中所示曲线分布可以看出，在液池内，颗粒浓度分布呈现多峰的分布特征。随着颗粒粒径的增大，颗粒的沉降作用增强，液池底部的颗粒浓度增大。其中，粒径为5～15 μm颗粒与气体的跟随性较好，因此在沿液池高度方向上，该粒径颗粒呈现相对较为平缓的分布趋势，只有在靠近下降管出口以下区域（即0.28 m高度），其浓度出现最大峰值。相比之下，80～110 μm 粒径颗粒与气体的跟随性较差，因此其浓度在液面以下区域（即0.52 m高度）出现较高峰值。对于380～530 μm粒径范围的颗粒，由于其粒径较大，惯性作用更强，沉降效果更好，因此其最大浓度峰值出现在液池底部。除此之外，由于液池内液体向上浮力及逆压力梯度[21-22]的作用下，该粒径颗粒在下降管出口处出现聚集并形成一个浓度峰值。从图中还可以看到，粒径较小颗粒在液池内部的浓度相对较高。而粒径较大的颗粒，因其沉降性较好，受气体和液体的影响较弱，故其浓度相对较低，整体波动相对较剧烈。

图6为气体流量对粒径为80～110 μm颗粒在液池高度方向上分布的影响。从图中所示曲线可以发现，一方面，随着气体流量的增加，使得该粒径颗粒更容易被气体所挟裹到液面附近悬浮；另一方面，气体流量的提高，使得颗粒冲击动量增加，惯性力增强，从而导致液池底部的颗粒浓度提高。因此，在气体流量为18 m3·h-1时，液池内沿高度方向上呈现在液池底部和液面附近颗粒浓度较高，而液池中间部位颗粒浓度较低，即两头高，中间低的分布趋势。

图6 不同气体流量下颗粒沿液池高度方向浓度分布Fig.6 Concentration distribution of particles along height direction of liquid path with different flow rate of gas(particle size：80—110 μm)

图7给出了气体流量为15 m3·h-1条件下，下降管出口静态浸没深度（即液体处于静止状态时，下降管出口到自由液面的垂直距离）对 5～15 μm粒径颗粒在液池高度方向上分布的影响。当下降管出口静态浸没深度为24 cm时，液池中的颗粒分布总体上较为平缓，只是在液面处（即高度h=0.6 m）的颗粒浓度出现最高峰值。随后由于液体回流作用，在位于液位和下降管出口之间（即高度h=0.44 m）再次出现峰值。在下降管出口下方（即高度h=0.28 m），由于颗粒的惯性作用，在此出现聚集，从而再次出现峰值。相比之下，当下降管出口静态浸没深度为2 cm时，液池内颗粒浓度分布波动较大。结合作者前期对该过程气液两相流动特性的研究[5]可知，当下降管出口静态浸没深度变小时，气体穿越液池所经历的行程缩短，气体对液池内液体的扰动会相应减弱，这使得液体内的颗粒浓度分布波动较大。而浸没深度的增加提高了液体的扰动程度，促进了颗粒运动，从而使得液池内的颗粒浓度分布相对平缓。

图8给出了80～110 μm粒径颗粒在液池不同高度上的径向分布。图中由于h=0.44 m高度以上正好有下降管，故其以上高度的径向测点是从0.04 m处开始，即下降管外壁面处。从图中所示曲线可以看出，由于h=0.44 m高度以上正好位于下降管出口以上区域。在该区域，径向颗粒浓度在靠近下降管外壁面区域出现峰值，并且在液面高度处（即h=0.6 m）呈现最大值。这是由于气体在穿越液池时，气泡主要在下降管外壁面附近聚集并向上流动，为此颗粒在气体的带动下在此区域出现较高浓度分布。h=0.36 m高度正好位于下降管出口下部区域，颗粒随气体从下降管出口处冲击进入液池。由于颗粒的惯性作用，在下降管出口下部区域出现较大的颗粒浓度值。此外，由于液体回流作用，随着液池高度的降低，径向颗粒浓度逐渐上升并且其波动幅度相对减弱。

3.3 气固两相流穿越液池运动过程的数值模拟

参考上述实验所选用的气液固三相系统以及所采用的颗粒粒径范围，表1给出了数值模拟所用到的主要计算参数和计算条件。数值模拟所选用的颗粒粒径及其加入级配如表2所示。

表1 计算参数Table 1 Calculation parameters

表2 入口颗粒级配Table 2 Particles size grading

图9显示了数值模拟得到的多种粒径颗粒共存时的气固两相流在穿越液池过程中不同粒径颗粒的运动过程及其轨迹。比较前面实验获得的不同粒径运动图像可知，数值模拟得到的不同粒径颗粒运动轨迹与实验结果较为接近，这说明数值模拟能够清晰地预测出气固两相流在穿越液池的过程中固体颗粒的运动演变过程。而所采用的确定性颗粒轨道模型，可以更加直观地观测到不同粒径颗粒细微的运动过程。

图9 液池内颗粒运动过程的数值模拟Fig.9 Numerical simulation results of movement of particles in liquid bath

图 10给出了数值模拟计算得到的某一时刻气固两相流穿越液池的气液固三相流动过程。由图可以清晰直观地看到不同粒径颗粒的运动过程。气固两相流在穿越液池的过程中，气体进入液池，在液池内形成大量气泡，而细小粒径颗粒被气泡所挟裹，并随气体一同钻出液面。同时，液面处气泡破碎所激发的液滴又将一部分已钻出液面的细小颗粒捕获，并在重力作用下一同回落至液池。此外，较大粒径的颗粒则在惯性作用下直接沉降到液池底部。还有一部分小粒径颗粒则在液池中悬浮并随着液体的回流而流动。可见，不同粒径颗粒在随气体穿越液池过程中分别经历了颗粒被气泡挟裹、颗粒在液体中悬浮和沉降、颗粒被液滴捕获和颗粒完成穿越液池等多个运动分离过程。最后只有少数细小粒径颗粒能随气体一同穿越液池并被气体携带进入液面以上空间而排出。

3.4 液池内颗粒浓度分布的数值模拟

图10 气固两相流穿越液池气液固三相流动Fig.10 Gas-liquid-solid three-phase flow in process of gas-solid flow passing through liquid bath

图11 不同粒径颗粒沿液池高度方向上的浓度分布的数值模拟结果Fig.11 Numerical simulation results of distribution of different size particles along height direction of liquid path

图 11给出了数值模拟得到的不同粒径颗粒在沿液池高度方向上的浓度分布。从图中可以看出，颗粒在液池内同样呈现多峰的波动分布规律。其中，10 μm粒径颗粒总体分布较为平缓。300 μm和90 μm 粒径颗粒浓度在液池中的波动较大。由于大粒径颗粒的沉降性较好，故其浓度的最大峰值出现在液池底部。而在液面以上的空间（即0.6 m高度以上），主要是10 μm颗粒占主，这说明该粒径颗粒存在随气体一道穿越液池的逃逸现象。

此外，对比数值模拟与实验结果，两者还存在一定的差异。分析其原因，在于不论是所加入的颗粒总量上，还是入口颗粒粒径分布上，由于计算时间限制，数值模拟还无法与实验实现精确一致。但是从数值模拟结果可以看出，所得到的颗粒运动分离过程以及颗粒浓度分布均能够反映出与实验一致的一般规律，即颗粒分布的波动性，其浓度峰值出现的位置也基本一致。这在一定程度上证明了本文所建立数学模型进行数值预测的合理性。

4 结 论

（1）基于Enler-Lagrange框架下的气液固三相数值模拟结果与实验结果均能反映出一致的气液固三相流动过程和固体颗粒的运动及其分布特征。

（2）数值模拟结果表明，颗粒在随气体一同穿越液池的过程中，不同粒径颗粒分别经历了被气泡挟裹、在液体中悬浮、在液体中沉降、被液滴捕获和穿越液池等多个运动分离过程。

（3）液池内不同粒径颗粒浓度沿液池高度方向和径向均呈现多峰的波动分布特征。

（4）5～15 μm粒径细小颗粒在整个液池内的浓度分布较为均匀。而380～530 μm粒径颗粒主要集中分布在液池底部。

（5）随着下降管出口静态浸没深度的增加，液池内颗粒浓度分布趋于平缓。

（6）随着气体流量的增大，液池内颗粒浓度呈现液池底部和液面处增高，而液池中部减低的分布态势。

符号说明

d——颗粒直径，m

Gij——颗粒i和j的相对速度，m·s-1

Ip——颗粒的转动惯量，kg·m-2

Mp——颗粒所受的合力矩，N·m

m——质量，kg

mtube——取样管空管质量，kg

u——速度，m·s-1

α——气相或液相含率

ρ——密度，kg·m-3

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

下角标

g ——气相

l——液相

p ——颗粒相

[1] Wang Fuchen (王辅臣), Yu Guangsuo (于广锁), Gong Xin (龚欣),et al. Research and development of large-scale coal gasification technology [J].Chemical Industry and Engineering Progress(化工进展), 2009, 28 (2)：173-180

[2] Jiang Han (江瀚), Liu Zhongliang (刘中良), Gong Xiaolong (宫小龙). Thermal economics and application of pressurized submerged combustion evaporators [J].CIESC Journal(化工学报), 2011, 62(12)：3498-3502

[3] Vaibhav K A, Ashwani K G. Jet characteristics from a submerged combustion system [J].Applied Energy, 2012, 89：246-253

[4] Wang Rongdong (王荣东), Du Hai’ou (杜海鸥), Gao Yaopeng (高耀鹏),et al. The experimental study on sodium aerosol removal efficiency of water bath scrubber [J].Nuclear Science and Engineering(核科学与工程), 2012, 32 (4)：319-325

[5] Wu Xuan (吴晅), Gao Huijie (高慧杰), Li Haiguang (李海广), Wu Wenfei (武文斐). Bubble characteristics of syngas during passing through cistern in gasifier [J].Journal of Chinese Society of Power Engineering(动力工程学报), 2013, 33 (12)：972-968

[6] Fu Bihua (付碧华), Wang Yifei (王亦飞), Lin Lan (林岚), Xu Jianliang (许建良),et al. Penetration depth and characteristics of interface fluctuation at quenching pipe outlet [J].CIESC Journal(化工学报), 2011, 62 (7)：1817-1823

[7] He Biyun (贺必云). Study on multiphase flow characteristics of the new type of scrubbing-cooling chamber [D]. Shanghai：East China University of Science & Technology, 2005：11-19

[8] Xu Yonggui, Liu Mingyan, Tang Can. Three-dimensional CFD-VOF-DPM simulations of effects of low-holdup particles on single-nozzle bubbling behavior in gas-liquid-solid systems [J].Chemical Engineering Journal, 2013, 222：292-306

[9] Baltussen M W, Seelen L J H, Kuipers J A M, Deen N G. Direct numerical simulations of gas-liquid-solid three phase flows [J].Chemical Engineering Science, 2013, 100：293-299

[10] Chen Min (陈敏), Yuan Zhulin (袁竹林), Gao Shiwang (郜时旺).Direct numerical simulation and experimental contrast of the moving granular bed dust filter [J].Proceedings of the CSEE(中国电机工程学报), 2003, 23 (11)：195-199

[11] Cai Guiying (蔡桂英), Yuan Zhulin (袁竹林). A study on the filtration performance of ceramic filter by numerical simulation [J].Proceedings of the CSEE(中国电机工程学报), 2003, 23 (12)：203-207

[12] Brackbill J U. A continuum method for modeling surface tension [J].Journal of Computational Physics, 1992, 100 (3)：335-354

[13] Ouyang J, Li J H. Particle-motion-resolved discrete model for simulating gas-solid fluidization [J].Chemical Engineering Science,1999, 54 (13)：2077-2083

[14] Xiong Y Q, Zhang M Y, Yuan Z L. Three-dimensional numerical simulation method for gas-solid injector [J].Powder Technology,2005, 160 (2)：180-189

[15] Kafui K D, Thornton C, Adams M J. Discrete particle-continuum fluid modeling of gas-solid fluidized beds [J].Chemical Engineering Science, 2002, 57 (13)：2395-2410

[16] Zou L M, Guo Y C, Chan C K. Cluster-based drag coefficient model for simulating gas-solid flow in a fast-fluidized bed [J].Chemical Engineering Science, 2008, 63 (4)：1052-1061

[17] Zhang Jianping, Fan L S, Chou Zhu,et al. Dynamic behavior of collision of elastic spheres in viscous fluids [J].Powder Technology,1999, 106 (1/2)：98-109

[18] Zhang Jianping, Li Yong, Fan Liang-Shih. Discrete phase simulation of gas-liquid-solid fluidization systems：single bubble rising behavior[J].Powder Technology, 2000, 113 (3)：310-326

[19] Tsuji Y, Tanaka T. Cluster patterns in circulating fluidized beds predicted by numerical simulation (discrete particle modelversustwo-fluid model) [J].Powder Technology, 1998, 95 (3)：254-264

[20] Li Jinghai (李静海), Ouyang Jie (欧阳洁), Gao Shiqiu (高士秋),et al.Multi-scale Simulation of Complex System of Particle Fluid (颗粒流体复杂系统的多尺度模拟) [M]. Beijing：Science Press, 2005：221

[21] Rouhani S Z. Two-phase flow patterns：a review of research results[J].Progress in Nuclear Energy, 1983, 11 (3)：219-259

[22] Wu Xuan (吴晅), Wang Lifang (王丽芳), Gao Huijie (高慧杰),et al.Numerical simulation of the initial impact behavioral characteristics during top-submerged gas injection [J].Boiler Technology(锅炉技术), 2013, 44 (3)：1-5