摘要 为研究气力输送发送罐中煤粉出料稳定性及两相流动机理，采用基于Realizable k-ε双方程湍流模型并结合颗粒动力学理论的Two Fluid Model（TFM）对其进行数值模拟。根据实验数据对Gidaspow模型中的曳力系数进行修订，得到自定义的曳力模型，使之符合实际工况。对煤粉在出料口的质量流量、体积分数的波动特性进行分析，发现随着入口气速增大，煤粉出料稳定性呈现出先增加后降低的规律。最后对发送罐中空气-煤粉两相耦合运动流动机理进行分析，发现罐中空气主要以气流涡的形式进行运动，具体表现为多种尺度、相互碰撞、叠加的涡驱动煤粉沿流向扩散和展向迁移。
关 键 词 气力输送；TFM；曳力系数；数值模拟；流动机理；涡流
中图分类号 TQ022.4 文献标志码 A
0 引言
气力输送煤粉广泛应用于电力、冶金和环保等行业中，而发送罐是实现此类粉体输送的核心装置之一[1]。气体输送煤粉涉及到复杂的气-固耦合流动，国内外关于煤粉气力输送发送罐的研究针对不同操作条件下的输送规律较多[2]，对发送罐内气固两相流动机理进行分析还比较少，尚未有较为全面深入的研究。崔越等[3]对栓塞式气力输送灰罐进行了数值模拟，得到了物料发送过程中可分为起始、稳定和下降3个阶段，发现稳定输送阶段基本不受输送压力的影响，但是发送罐出料的稳定性会直接影响气力输送过程中管路内的流型[4-5]。徐贵玲，陈晓平等[6，7]在不同压力和物性条件的上出料式发送罐气力输送试验台上，以氮气为输送介质，研究煤粉气力输送特性，与高压输送相比，常压输送的固气质量比较高，气耗率及输送能耗较低。贺春辉，沈湘林等[8]以兖州烟煤煤粉作输送物料，发现上引式出料气力输送系统的煤粉质量流量、表观气速和固气质量比均略高于下引式出料系统，但二者均具有良好的输送稳定性。可见，不同工况下的研究会得出不同的结果，研究成果存在局限性。因此，对于罐内两相流动的机理研究显得尤为重要。
基于欧拉-欧拉框架的TFM模型已经广泛的应用于流化床[9]、喷动床[10]的数值模拟，其在模拟气固两相流的可行性也已经被验证[11-12]。颗粒动力学理论（KTGF）是基于稠密的气体分子动力学理论而建立的描述颗粒相运动的方程。该方程引入了颗粒温度和径向分布函数来描述颗粒相的速度脉动与颗粒碰撞等。该理论中[13]，颗粒相被处理成拟流体，将颗粒运动类比于气体分子的热运动，颗粒的随机运动动能被定义为“颗粒拟温度”，其正比于颗粒速度脉动的均方。本文主要分析发送罐内气固两相的流动特性，属于介观层次[14]的考察，选用TFM模型和KTGF相耦合的方法可满足模拟需求。
1 模型描述
1.1 TFM
在欧拉框架下构建的双流体模型属于比较复杂的多相流模型。它把每一相看做具有不同运动特性的连续性介质，而且各相之间存在相互渗透、耦合的现象，从而建立了一套完整的包含有各相的连续性方程与动量方程的方程组进行求解。方程把各界面间的交换系数与压力项耦合到一起，而耦合的方式则根据所含相的情况进行相应的选择，譬如颗粒流与非颗粒流的处理方式就有所不同。针对颗粒流，可以采用颗粒流体力学中的分子运动理论来得到流动特性。跟单流体模型比较而言，TFM优势在于考虑了气固两相间由于存在相互滑移而引起的阻力以及固相的湍流输运。因此在大多数情形下，其计算结果与实际情况更加吻合。由于每一相都有其单独的守恒方程，各相可以进行单独的计算，因此该模型适应性很强。当然，该方法也有其局限性。显然，方程数的增加必然导致计算量的增加。此外，顆粒相的导热系数，扩散系数等物性参数的确定依然存在一些问题，双流体模型在计算过程中也可能会发生伪扩散现象。
连续性方程（质量守恒方程）
[∂∂tαqρq+∇·αqρqυq=p=1nmpq]，（1）
式中：[αq]为q相的体积分数；[ρq]为q相的密度； [υq]为q相的速度； [mpq]为从第p相到q相的质量传递。
动量守恒方程
[∂∂tαqρqυq+∇·αqρqυqυq=-αq∇p+∇·τq+p=1nRpq+mpqυpq+αqρqFq+Flift，q+FVm，q]，（2）
式中：[τq]是第q相的压力应变张量；[Fq]是外部体积力；[Flift，q]是升力；[FVm，q]是虚拟质量力；[Rpq]是相间作用力；p为所有相共享的压力；[υpq]为相间速度。
1.2 颗粒动力学模型
引入颗粒流体力学来处理发送罐中颗粒相的相关问题。颗粒动力学理论源自于气体分子运动理论，借助速度分布函数的Boltzmann方程来推导颗粒相的物性参数及守恒方程，与热力学温度表征微观尺度分子的能量波动类似，颗粒相拟温度表达了宏观尺度颗粒随机运动的动能，通过颗粒相拟温度Θ来建立颗粒相宏观的黏度、压力、扩散系数等参数的本构关系：
[τs=[-ps+αsξs∇∙us]I-2αsμsSs]，（3）
[Ss=12∇us+∇usT-13∇∙usI]，（4）
[μs=10ρsdsΘπ0.596αs1+eg01+45g0αs1+e2+45αsρsdsg01+eΘπ0.5]，（5）
[ξs=43α2sρsdsg01+eΘπ0.5]，（6）
[ps=αsρsΘ+2ρs（1+e）α2sg0Θ]，（7）
式中：[ξs]是固相体积黏度，源于颗粒之间碰撞造成的动量交换；[μs]是颗粒相剪切黏度，由动力部分和碰撞部分组成；[ps]是颗粒相压力。
颗粒相拟温度[Θ]和颗粒随机运动的动能成正比，颗粒相脉动能守恒方程
[32∂∂tαsρsΘ+∇∙αsρsusΘ=τs：∇us+∇∙q-γ-3βΘ]，（8）
其中：
[q=k∇Q ，k=2rsa2sds1+eg0Qp0.5+150rsdsQp0.53841+eg01+65asg01+e2 ，g=31-e2as2rsg0Q4dsQp0.5 ，g0=1-asas，max13-1 ，]（9）
式（8）中：[Θ]是颗粒相拟温度，定义为[Θ=13u"su"s]；[u"s]为颗粒的脉动速度，表征颗粒之间因碰撞而产生的小尺度脉动强度；式右边第1项代表由于变形作用导致的能量生成；第2项是能量扩散；第3项是由于非弹性碰撞导致的能量耗散；最后1项是在2个相之间脉动能量的净转移率；[γ]是碰撞能量耗散；q是脉动能量通量；k是脉动能量传导率；e是颗粒碰撞的恢复系数；g0是径向分布函数；[αs，max]是对于任意填充的最大固体体积分数。
1.3 曳力模型
当黏性流体绕流物体时，物体会受到压力和切向应力的作用。在沿着物体横截面的流动平面中，物体所受的合力F可以分解为与来流方向同向的作用力FD和垂直于来流方向的升力FL。FD与流体和物体之间的相对速度有关，当物体运动速度大于流体流速时，FD即为物体在流体中受到的阻力。反之，当流体流速大于物体运动速度时，FD受到的就是流体的曳力。
本文对气固两相间的曳力系数β选用Gidaspow模型，该模型适用于密相气固流。当αg<0.8时系数[β]采用Ergun提出的公式，而当αg≥0.8时系数β采用Wen、Yu的公式：
[β=150α2sμgαgd2s+1.75ρgαsug-usds ，αg<0.8，β=34Cdαsαgρgug-usdsα-2.65g ，αg≥0.8，]（10）
[Cd=24Res1+0.15Res0.687 ，Res<1 000 ，Cd=0.44 ，Res≥1 000 ，]（11）
[Res=αsρgug-usdsμg]。（12）
2 数值计算说明
2.1 几何模型及网格划分
以上引式气力输送发送罐作为研究对象，几何模型尺寸如图1所示。模型整体高度为700 mm；填料高度为550 mm；上端直径为200 mm；底部为空气入口，直径为50 mm；出料管位于模型中心处，竖直向上，直径为10 mm。由图1所示可知，发送罐模型为轴对称模型，内部流场流动情况亦呈对称分布，因此取整体模型四分之一进行研究。利用ICEM软件对几何模型进行结构网格划分，如图2所示。
2.2 网格无关性检验
对于气固两相流流动，为使网格捕捉到所需流场的细节，要求网格的最大尺寸应低于最小的细节特征，即网格尺寸要小于颗粒直径的1/10[15]。如图3所示，分别划分10万、40万、80万、100万网格，发现80万网格满足求解需要，继续增加网格数对模拟结果影响较小。
2.3 数值模拟条件设置
容器底部进气口类型设置为速度入口，湍流强度和水力直径分别为5%，50 mm；出料口类型设置为出流条件；壁面条件选用无滑移壁面；初始物料填充高度设置为550 mm；采用SIMPLE算法即求解压力耦合方程组的半隐式方法进行求解，空间离散格式选用具有较高精度的二阶迎风格式和QUICK格式，防止一阶精度导致的数值上的扩散误差。湍流模型选用Realizable k-ε双方程模型，但该湍流模型适用于充分发展的湍流，不能应用于近壁面区域，因此要把第1层网格的计算节点放置在湍流核心区内。而过渡层与粘性底层区域中的物理量由壁面函数来计算完成。[Y+]的计算公式：[Y+=u*y/ν]；式中u*是近壁面摩擦速度；y为第1层网格节点与壁面的间距；[ν]为流体的运动粘度。在计算过程中设置边界层，第1层的高度为0.1 mm，增长率为1.1，计算得到[Y+]在30～300之间，满足壁面函数法的要求。
3 模拟的结果与讨论
3.1 曳力系数的修正和误差分析
首先将模拟数据与实验数据进行校正，由图4可以发现Gidaspow模型在速度较低时与实验数据的趋势一致，在数值上低于实验数据；但当速度增大到0.113 m/s时，需要对Gidaspow模型中的曳力系数进行修正。校正之后如图5所示，模拟数据与实验数据基本吻合。
从入口速度与出口质量流量的关系可以看出，随着入口气速的增大，出口的质量流量先增大后减小，在u=0.118 m/s时，出口质量流量达到最大值为0.021 kg/s，再增大入口气速，出口的质量流量逐渐降低。此外，衡量输送能力强弱还与输送总耗时有关。从图6可以看出，输送总耗时在速度较低的时候随着入口气速的增大而减小，当超过某一值（u=0.12 m/s）时，输送总耗时基本保持不变，最终维持在85 s左右。
这对于工程中实际的意义是获取最具性价比的操作条件。所以，以增大入口气体速度来提高输送量的方式前期作用效果明显，但当速度达到一定值之时，改善效果不再明显，只会造成成本显著提高。
3.2 出口体积分数和质量流量的波动
工程运用上，输送的稳定性一直是关注的热点。本节主要对输送过程中出口煤粉的质量流量和体积分数的波动进行分析。图7是不同入口气速下出口质量流量随着时间的变化，图8是不同入口气速下出口固体体积分数随时间的变化。从图中可以看出质量流量和体积分数都呈上下波动状，其波动程度可用标准差来表征。标准差值越大表示波动程度越大，输送过程越不稳定。波动程度剧烈会增大对管路的冲击，造成管路颤动，甚至是引发安全事故。
从图7可以看出，在入口气速较低的时候，出口的质量流量波动的整体振幅较小，随着入口气速的增大，波动的振幅也逐漸增大，其波动情况与入口气体速度不是单纯的线性关系。从图8中可以看出出口固体体积分数随气速变化并不明显，在速度为0.1 m/s时，波动程度相对较小。
从图9可以看出，出口固体体积分数的标准差随着入口气速的增大，先减小后增大；而出口固体体积分数的均方根则随着入口气速的增大，先增大后减小。通过对出口固体体积分数的均方根（αRMS）和标准差（αS）进行曲线拟合可得出方程
[αRMS=0.207+u-4.10u2]，（13）
[αS=0.15-0.50u+2.19u2。]（14）
通过计算可以得出以下结论：在入口气速u=0.10～0.12 m/s时，出口的质量流量最大、出口的固体体积分数的均方根最大并且出口固体体积分数的标准差最小。
3.3 发送罐内气固两相流流动机理分析
取三维模型1/2中心对称面作为观测面。图中X代表气力输送发送罐轴向尺寸，Y代表气力输送发送罐径向尺寸。如图10、图11所示，进入发送罐的空气主要沿着阻力较小的罐体壁面进行运动，从流线运动规律可以看出，在壁面处形成强烈漩涡，这些强烈漩涡对罐内平行于Z轴截面上的煤粉运动产生了两个方向的作用力：一个是展向的运动，即Y方向，该方向作用力使煤粉沿着展向进行移动，增强煤粉与空气的混合程度；另一个是流向的运动，即X方向，该方向的作用力是煤粉运动的主动力，力越大煤粉运动的越快，进而进入输送管的流量越大。
从空气流线图可以看出，由入口进入的空气可分为两部分，一部分穿过罐中物料沉积层，加速物料整体的流动，另一部分沿着壁面运动，在壁面附近与煤粉进行能量传递，涡流现象非常明显。从图中可以发现，在接近发送罐顶部的空气域存在两个旋转方向相反的空气漩涡，因为沿着管路运动的空气流首先达到罐体顶部，而沿着罐体壁面的空气流随后到达罐体顶部区域，这两股空气流在罐体顶部区域相遇，形成了外层逆时针方向和内层顺时针方向的空气涡。从煤粉流线图可以看出，首先煤粉受重力的作用整体往下运动，其次受到气流的曳力作用，煤粉的密度分布发生变化。随着输送时间的增加，罐中物料不断减少，顶部的空气涡尺度也随之越来越大。从图中可以看出煤粉在流向和展向上的运动速度都受到了空气涡强度和方向的影响。
对于流动产生大量的漩涡可以从两方面进行考滤。一方面，大量的漩涡会造成剧烈的湍流耗散，增大能量的消耗；另一方面，漩涡可以增大湍流强度，增强脉动程度，使气固之间的混合更加均匀，有利于两相之间的动量和能量传递。如图12所示，在轴向Z = 550 mm的截面观测煤粉沿径向的运动情况，可以看出煤粉的运动受到空气流的影响，沿Z轴方向翻滚，呈沸腾状，加速了气固两相间的动量和能量的交换。
4 结论
本文采用结合了颗粒动力学理论和Realizable k-ε的欧拉双流体模型，对上引式气力输送发送罐的流动机理和出料过程流动特性进行了数值模拟研究。对Gidaspow曳力模型进行了修正，使模拟结果与实验数据相吻合。对出料口固体质量流量和固体体积分数的波动特性进行分析，考察了不同入口气速下对出料稳定性的影响，在此基础上，对罐内气固两相的流动机理进行分析，得出以下结论：
1）随入口气速的增加，出料口煤粉的质量流量和体积分数都呈现先增大后减小的变化趋势，其稳定性呈现先增加后降低的变化趋势，入口气速在v=0.10～0.12 m/s范围内存在一个最优值，可使波动最小、煤粉体积分数和出料量最大。
2）初始阶段随着入口气速的增加，输送物料所需时间逐渐减小，但当入口气速超过0.12 m/s时，耗时基本保持恒定。
3）发送罐中的气体运动以各种类型的涡流运动为主，作用方式是各种尺度涡的叠加、碰撞和干涉；运动形式分解为流向和展向，驱动煤粉的展向扩散和流向移动，展向扩散加速煤粉和空气的混合程度，气流曳力驱动煤粉离开输送罐。
参考文献：
[1] 杨伦，谢一华. 气力输送工程 [M]. 北京：机械工业出版社，2006：172-179.
[2] 陆慧林，刘文铁，孙永立，等. 稠密气固两相湍流流动的实验和数值模拟 [J]. 力学学报，2000，32（4）：385-391.
[3] YUE C，HONG G，SUN J，et al. Numerical simulation of gas-solid flow in a conveying vessel [J]. Powder Technology， 2012，226（8）：34-42.
[4] JONES M G，WILLIAMS K C. Predicting the mode of flow in pneumatic conveying systems [J]. Particuology， 2008，6（5）：289-300.
[5] Setia G，Mallick S S. Modelling fluidized dense-phase pneumatic conveying of fly ash[J]. Powder Technology， 2015，270：39-45.
[6] 徐贵玲，陈晓平，梁财，等. 煤粉上出料式发送罐气力输送特性[J]. 化工学报，2012，63（6）：1709-1716.
[7] 徐贵玲，卢平，许盼，等. 煤粉外水含量对上出料式发送罐供料特性的影响[J]. 化工学报，2016，67（7）：2767-2776.
[8] 贺春辉，陈先梅，倪红亮，等. 发送罐出料方式对煤粉高压密相气力输送特性的影响[J]. 中国电机工程学报，2012，32（23）：32-39.
[9] ZHANG Y，LEI F，WANG S，et al. A numerical study of gas-solid flow hydrodynamics in a riser under dense suspension upflow regime[J]. Powder Technology，2015，280 （20）：227-238.
[10] FATTAHI M，HOSSEINI S H，AHMADI G. CFD simulation of transient gas to particle heat transfer for fluidized and spouted regimes[J]. Applied Thermal Engineering，2015，105（6）：385-396.
[11] 沈志恒，李巍，陆慧林，等. 不同曳力模型对提升管内气固流动特性的影响[J]. 石油机械，2011，39 （10）：32-38.
[12] CLOETE S，JOHANSEN S T，AMINI S. Investigation into the effect of simulating a 3D cylindrical fluidized bed reactor on a 2D plane [J]. Powder Technology，2013，239（17）：21-35.
[13] CHAN C K，GUO Y C，LAU K S. Numerical modeling of gas-particle flow using a comprehensive kinetic theory with turbulence modulation [J]. Powder Technology，2005，150（1）：42-55.
[14] TSUJI Y. Multi-scale modeling of dense phase gas-particle flow[J]. Chemical Engineering Science，2007，62（13）：3410-3418.
[15] AGRAWAL K，LOEZOS P N，SYAMLAL M，et al. The role of meso-scale structures in rapid gas-solid flows[J]. Journal of Fluid Mechanics，2001，445：151-185.
[責任编辑 田 丰]