垂直弯管段高压密相煤粉气力输送特性的数值模拟

裴 宇，熊源泉，周海军

（东南大学能源热转换及其过程测控教育部重点实验室，南京210096）

高压超浓相煤粉气力输送是大型气流床高压气化的关键技术之一，由于受到目前试验方法和条件的限制，还不能实时准确地测量管道中气固两相的流动参数，因此对高压密相煤粉气力输送过程进行数值模拟研究有利于输送系统的设计与优化.

类比稠密气体分子运动论，Savage等、Jenkins等、Shahinpour 等、Lun 等、Johnson 等和Syamlal发展了颗粒动力学理论，引入颗粒碰撞恢复系数的概念和稠密气体的Boltzmann输运方程，修正了经典的分子运动论对两颗粒碰撞的处理，推导出颗粒相应力的表达式和适合于描述颗粒流动的封闭方程组.颗粒动力学理论的引入是近几十年来双流体模型最重要的发展.在高压密相模拟中，水平管和垂直管的研究取得了一定进展［1-2］，但弯管方面的模拟研究鲜有报道.Mcglinchey等［3］在颗粒动力学基础上提出了常压下90°弯管的计算模型（以下简称Mcglinchey弯管模型），该模型中忽略了摩擦应力，而且固相壁面采用的是无滑移条件，模拟的压降结果误差在10%～90%内.Ma等［4］运用Mcglinchey弯管模型研究了粒径对水平转垂直向上弯管压降的影响，模拟得到的整体管道压降与试验较吻合，但就弯管段而言，管道压降误差高达50%.Vashisth等［5］采用镜面系数和颗粒与壁面碰撞恢复系数来表示壁面剪切和拟热流条件，对Mcglinchey弯管模型进行了修正，模拟获得的弯管固相体积分数和速度分布云图与实际较吻合，但没有对压降进行分析.

在本文高压密相煤粉气力输送管道内，离心加速度垂直指向管道壁面外侧，而且远大于重力加速度，是垂直弯管气力输送过程中影响气固流动的最主要因素.在输送过程中，离心力的作用导致弯管壁面外侧煤粉堆积形成堆积层，并以颗粒群的形式运动，而弯管壁面内侧颗粒呈现悬浮流动.悬浮流动时，固相应力主要来源于颗粒之间的瞬时碰撞和平移产生的动量传递，这在颗粒动力学理论中已阐述清楚；在颗粒堆积层内，颗粒间除了瞬时碰撞和平移外，由于离心力作用，沉积层内的颗粒相互挤压且长时间接触产生了摩擦应力，并通过摩擦应力进行颗粒间动量和能量的传递与耗散.由于离心力较大，堆积层顶部体积分数接近自然堆积体积分数，颗粒与颗粒以及颗粒与壁面都保持长时间接触，摩擦应力成为堆积层内固相应力和壁面剪切应力的重要部分，必须加以充分考虑.

以往文献在运用运动颗粒动力学解决弯管问题时，在颗粒内和颗粒与壁面处忽略了摩擦应力的影响，造成数值模拟产生较大的误差.为此，笔者在颗粒动力学的基础上考虑了摩擦应力的影响，通过一体化管道对垂直向上转水平弯管（其曲率半径为0.2m，管道内径为10mm，固相速度为5～8m／s）的高压密相煤粉气力输送进行模拟，并结合多组试验工况，验证了模型的准确性，同时分析了垂直弯管内气固流动情况以及补充风体积流量对输送特性的影响.

1 数学模型

1.1 控制方程

连续方程

动量方程

式中：φ为各相体积分数；ρ为密度；v为速度矢量；τ为应力张量；pg为气相压力；Fsg表示气固相间作用力；g为重力加速度；下标g、s分别表示气相、固相.

在密相煤粉气力输送过程中气固相间作用构成复杂，由于气固速度差而产生的曳力是相间作用的主导因素，所以本文模型中的气固相间作用只考虑曳力.

式中：β为气固相间的曳力系数，大量参考文献采用Gidaspow 提出的曳力模型［6］（以下简称Gidaspow模型）.

当φg≤0.8时，β采用Ergun提出的适用于密相的曳力模型：

式中：μg为气体黏性系数；ds为颗粒直径.

当φg＞0.8时，β采用Wen等提出的适用于稀相的模型：

式中：CD为流体阻力系数.

但φg＞0.8（即φs≤0.2）时并不完全是稀相流动，而且Gidaspow 模型在φg＝0.8 时不具有连续性，因此Lu等［7］对Gidaspow 模型进行了修正，提出Huilin-Gidaspow 模型：

1.2 颗粒动力学理论

Savage等认为颗粒相内的应力可以看成是动力学应力与摩擦应力之和，即，其中为颗粒动力学应力张量，为摩擦应力.

式中：为颗粒动力学应力张量中的固相压力；为颗粒相黏性系数；λks为颗粒相体积黏性系数；Ss为颗粒相速度变形率张量；I为单位张量.

固相压力由颗粒运动引起的动力效应和颗粒碰撞引起的碰撞效应组成：

颗粒相黏性系数包括颗粒相碰撞黏度和动力黏度2部分：

式中：Θs为颗粒拟温度；ess为颗粒间碰撞恢复系数；g0，ss为颗粒的径向分布函数.

颗粒相体积黏性系数表示颗粒相对压缩的悬浮阻力，其表达式为：

径向分布函数为：

颗粒拟温度的输运方程如下：

式中：qs为拟热流密度，在颗粒相内，表示沿颗粒拟温度Θs梯度方向的拟热流密度，其中kΘs为扩散系数；表示颗粒相内由剪切应力产生的颗粒脉动；γΘs表示颗粒与颗粒之间碰撞导致的颗粒脉动的耗散；φsg表示气体与颗粒间脉动的交换［1-2，7］.

Johnson等的摩擦应力模型认为固相体积分数达到0.5后摩擦应力才起作用，但高压密相煤粉气力输送中，固相体积分数低于0.5时也存在摩擦应力.因此，笔者认为固相体积分数达到0.1后摩擦应力开始起作用，并根据实际情况修正经验参数，具体表达式如下：

式中：pf为摩擦应力；φ为颗粒间的摩擦角；经验参数F＝0.1，r＝2，s＝5，φs，min＝0.1，φs，max＝0.55.

在气固模拟中气固湍流输运方程较统一，具体可参考文献［1］、文献［7］和文献［8］.

1.3 边界条件

1.3.1 入口和出口条件

对于气相，入口处径向速度分布为光滑管充分发展湍流流动，出口处设置为压力出口.

式中：Ug为表观气速；φs，in为入口处固相体积分数；D为管道半径.

对于固相，入口处径向速度为均匀分布，出口处设置为压力出口.

其中，qm，s为固相质量流量.

1.3.2 壁面条件的修正

对于气相，壁面采用无滑移边界条件.对于固相，壁面剪切应力仅考虑颗粒与壁面通过碰撞产生的颗粒动力学应力.

式中：usw为颗粒与壁面的滑移速度，下标w 表示壁面；φ′为镜面系数，与管道材料、颗粒类型和壁面加工条件等因素有关，本文弯管段的镜面系数设定为0.05，水平段和垂直段的镜面系数设定为1×10-5.

式中：esw为颗粒与壁面碰撞恢复系数，取值为0.6.

Jenkins等［9］认为颗粒流动分为2种区域：快速剪切区域和缓慢剪切区域.快速剪切区域中固相体积分数不高，固相应力主要来源于颗粒碰撞与平移导致的动量传递；缓慢剪切区域中固相体积分数高，单个颗粒与周围接触面紧密接触产生相互作用，固相应力主要来源于颗粒与周围接触面长时间接触而产生的静摩擦应力.离心加速度是影响流动的最主要因素，本文中垂直弯管的离心加速度远大于重力加速度，这使得颗粒与壁面的接触压力较大，此时摩擦应力对壁面剪切应力的影响不容忽视.Johnson等［10］在流化床模型中将摩擦应力表示为·tanδ，其中pf为摩擦正应力即摩擦压力，δ为颗粒与壁面的摩擦角.引用该表达式，可得到考虑壁面摩擦应力的壁面剪切应力表达式：

此时，壁面拟热流密度可采用式（24）进行计算.

2 模拟条件

如图1所示，模拟中采用一体化管道，管道内径为10mm，其中垂直管长2m、水平管长4m、垂直弯管曲率半径为0.2m，管道截面网格中轴向采用六面体cooper，轴间距为2mm，共计56.8万个网格.

图1 管道示意图和截面网格Fig.1 Piping schematic and the cross-section grid

3 试验装置

高压密相煤粉气力输送试验装置见图2.高压氮气经缓冲罐分为流化风、充压风和补充风3路，发料罐内的煤粉经流化风、充压风驱动进入输送管道，在发料罐出口引入补充风调节输送过程中的固气比（即固相与气相质量比），充压风用来维持发料罐内的压力.进入受料罐的含粉气流经布袋除尘器分离煤粉后放空，受料罐压力由电动调节阀调节.在输送管的水平段、垂直段和垂直弯管段均设有差压变送器来测量压差.在高压发料罐上安装电子秤，测量煤粉质量流量.试验装置上压力、流量、压差和煤粉质量流率数据均采用计算机数据采集系统采集.输送管道内径为10mm，输送距离为53m.

图2 试验系统和发料罐示意图Fig.2 Schematic diagram of the experimental system and supply tank

4 结果与讨论

选取内蒙褐煤的输送试验进行模拟，试验中保持背压2.5 MPa、总输送压差0.5 MPa和流化风体积流量0.4m3／h不变，通过改变补充风体积流量来改变输送固气比；模拟中验证了不同固气比下模型的准确度，其中褐煤全水分质量分数为5.45%，密度为1 496kg／m3，平均粒径为204.3μm，选取管程中间段的6m 为几何模型，具体参数设置见表1.

表1 模型参数设置Tab.1 Model parameter settings

4.1 考虑壁面摩擦应力对模拟结果的影响

分别采用未考虑摩擦应力和考虑摩擦应力的壁面条件对4组工况进行模拟对比，压降模拟结果见表2.

表2 考虑摩擦应力对垂直弯管压降的影响Tab.2 Effect of friction stress on pressure drop of the vertical elbow

注：1）表示未考虑摩擦应力时模拟得到的垂直弯管压降与试验值的误差；2）表示考虑摩擦应力时模拟得到的垂直弯管压降与试验值的误差.

由表2可知，4组工况下采用未考虑摩擦应力的壁面条件模拟得到的压降与试验值间的误差均较大，远超出合理范围，采用考虑摩擦应力的壁面条件模拟得到的压降均大于采用未考虑摩擦应力的壁面条件模拟得到的压降，计算准确性得到较大改善，基本达到了模拟精度的要求.

针对是否考虑摩擦应力对壁面剪切应力的影响较大，以补充风体积流量为0.6m3／h的工况为例，图3给出了45°截面壁面剪切应力.从图3 可以看出，未修正前的壁面剪切应力只考虑颗粒与壁面的碰撞，最大值为180Pa；而修正后的壁面剪切应力引入了摩擦应力的影响，最大值增大到240Pa，壁面的流动阻力增大，相应的垂直弯管压降从9.03 kPa增大到12.25kPa.

图3 考虑摩擦应力对45°截面壁面剪切应力的影响Fig.3 Effect of friction stress on the wall shear stress of 45°crosssection

图4给出了45°截面沿径向分布的固相速度和固相湍动能.从图4（a）可以看出，由于修正后的壁面条件考虑了颗粒与壁面长时间接触而产生的摩擦应力，从而增大了壁面剪切应力，增大了壁面流动阻力，造成颗粒密集区域速度减小.对于气相而言，颗粒密集区域流动阻力增大，速度减小，而在颗粒稀疏区域加速带动固相速度增大.从图4（b）可以看出，由于壁面是固相湍动能产生和耗散的重要部分，考虑摩擦应力后，颗粒在壁面处的剪切应力增大，壁面拟热流密度中的产生项随之增大，壁面拟热流密度也随之增大，壁面产生的湍动能经颗粒间传递造成截面颗粒密集区域的湍动能增大，而在颗粒稀疏区域，由于气固速度增大，颗粒间曳力增大，气固湍动能交换更频繁，固相湍动能随之增大，因此考虑壁面摩擦应力后弯管内的固相湍动能均有所增大.

图4 考虑摩擦应力对45°截面沿径向分布的固相速度和固相湍动能的影响Fig.4 Effect of friction stress on the radial solid velocity distribution and turbulent kinetic energy of 45°cross-section

4.2 弯管整体流动过程的分析

以补充风体积流量为0.6m3／h工况为例，对弯管整体流动过程进行分析.图5给出了沿弯管轴向的固相体积分数分布云图.图6给出了从弯管入口到中点每隔5°各截面固相体积分数云图.从图5和图6可以看出，颗粒相在入口段中心体积分数达到最大值，进入弯管后由于离心力的作用，颗粒相逐渐向弯管外壁移动形成高体积分数区，内壁处低体积分数区逐渐扩大.

在到达15°截面之前，靠近外壁面的部分颗粒开始与壁面发生碰撞，反弹后以一定的动量冲向未发生碰撞的颗粒群，此时颗粒群依旧保持垂直向上的速度且煤粉质量较大，反弹的颗粒无法突破颗粒群，会随着颗粒群再次向上运动，所以在15°截面之前外壁面的高体积分数区和内壁面的低体积分数区一直在扩大.在15°截面之后，颗粒群与外壁面发生碰撞并开始向内壁面反弹，由于低体积分数区中心气流速度较大，携带能力较强，反弹到中心的颗粒被气流迅速带走，而反弹到壁面的颗粒开始聚集并向低体积分数区中心移动；在15°截面到30°截面之间，从图6（d）～图6（g）可以清晰地看出，低体积分数区两侧的固相体积分数开始升高，颗粒向低体积分数区中心发展.当发展到45°截面时，高体积分数区颗粒与壁面的碰撞、颗粒与载气的曳力和离心力已经达到平衡，流形基本稳定，低体积分数区形成了倒“V”形的稳定区域.

图5 沿弯管轴向的固相体积分数分布云图Fig.5 Solid concentration contours along axial direction of elbow

图6 垂直弯管各截面固相体积分数云图Fig.6 Solid concentration cloud on various sections of vertical elbow

4.3 补充风体积流量对垂直弯管压降的影响

补充风体积流量增大导致输送管道补充风入口处的压力升高，减小了发料罐与补充风入口处的压差，可视为输送物料运动的能量减少，所以煤粉质量流量与体积通量均逐渐减小.另一方面，发料罐出煤量减少，用于填充煤粉流出后罐内形成空隙的冲压风体积流量减小，因此随补充风体积流量的增大，表观气速并非线性增大.补充风体积流量较小时，补充风体积流量增大对出煤量的影响不大，充压风体积流量变化不大，此时随补充风体积流量的增大，表观气速增大较快；补充风体积流量较大时，补充风体积流量增大对出煤量的影响较大，出煤量和补充风体积流量明显减小，此时随着补充风体积流量的增大，表观气速增大趋势减缓（见表3）.采用考虑摩擦应力的壁面条件模拟这4组工况，模拟中其他条件保持不变，改变入口处气相平均速度、固相平均速度和煤粉体积分数.

图7给出了垂直弯管压降随补充风体积流量的变化.从图7可以看出，修正后的垂直弯管段模拟所得压降与试验值的误差减小至20%以内，基本符合模拟要求.输送中垂直弯管段压降由气相压降和固相压降2部分组成，其中固相压降占主导地位，固相压降又分为动能压降、重力压降和摩擦压降，其表达式如下：

表3 不同补充风体积流量下对应的输送参数Tab.3 Conveying parameters at different supplementary air flows

式中：Gs为煤粉质量通量；H为高度变化量；f为固相摩擦压损系数；Us为固相表观速度.

图7 垂直弯管压降随补充风体积流量的变化Fig.7 Changes of vertical elbow pressure drop with supplementary air flow

4.4 补充风体积流量对垂直弯管流动特性的影响

选取垂直弯管段中点45°截面进行分析，各工况下的固相体积分数分布云图见图8.从图8可以看出，随着补充风体积流量的增大，截面煤粉质量通量减小，输送固气比降低，在云图中表现为体积分数分布变窄；补充风体积流量增大造成表观气速随之增大，从而增强了载气与颗粒间的曳力作用，煤粉更容易被携带，高体积分数区逐渐减小.

图8 不同补充风体积流量下垂直弯管段45°截面的固相体积分数云图Fig.8 Solid concentration cloud on 45°cross-section at different supplementary air flows

图9给出了各工况下垂直弯管段45°截面的固相湍动能和固相拟温度径向分布.从9（a）可以看出，垂直弯管段截面存在2个高固相湍动能区域：贴近壁面外侧的高体积分数区和贴近壁面内侧的低体积分数区.

随着补充风体积流量的增大，表观气速增大，截面煤粉质量通量减小，在高体积分数区，补充风体积流量不大时，贴近壁面处的固相体积分数最大值不变（见图8），此时固相速度增大导致颗粒与壁面的滑移速度增大，颗粒与壁面的剪切应力、壁面拟热流密度qs｜w和固相湍动能增大，当补充风体积流量达到1.0m3／h时，贴近壁面处的固相体积分数最大值开始减小，此时颗粒与壁面的碰撞概率和摩擦应力均减小，壁面拟热流密度qs｜w和固相湍动能也开始减小.

图9 各工况下垂直弯管段45°截面固相湍动能和固相拟温度径向分布Fig.9 Solid-phase turbulent kinetic energy and granular temperature distribution on 45°cross-section under different conditions

在低体积分数区，气固相间脉动传递强烈，表观气速增大对固相湍动能的增强效果更加明显，在补充风体积流量达到1.0m3／h时，截面煤粉质量通量急剧减小，造成内壁面低体积分数区煤粉体积分数进一步减小，此时该区域湍流强度有所减弱.固相拟温度曲线也存在着与固相湍动能相似的变化规律.

5 结 论

（1）考虑颗粒与壁面长时间接触而产生的摩擦应力后，模拟所得的垂直弯管压降与试验值的误差减小至20%以内，相比采用未考虑壁面摩擦应力时模拟所得的结果，计算精度得到较大改善.

（2）随着补充风体积流量的增大，垂直弯管压降先增大后减小，表观气速增大，垂直弯管壁面外侧煤粉堆积减少，低体积分数区范围增大，固相湍动能和固相拟温度均有所增大，当补充风体积流量达到1.0m3／h时，由于固相体积分数急剧下降，固相湍动能和固相拟温度均减小.

［1］蒲文灏，赵长遂.水平管加压密相煤粉气力输送沙丘流的数值模拟［J］.动力工程学报，2011，31（2）：120-126.

PU Wenhao，ZHAO Changsui.Numerical simulation on dune flow in pressurized dense phase pneumatic conveying of pulverized coal in horizontal pipe［J］.Journal of Chinese Society of Power Engineering，2011，31（2）：120-126.

［2］蒲文灏，赵长遂，熊源泉，等.垂直管密相输送的数值模拟［J］.动力工程，2008，28（1）：95-99.

PU Wenhao，ZHAO Changsui，XIONG Yuanquan，etal.Numerical simulation of the dense phase conveying in vertical pipe［J］.Journal of Power Engineering，2008，28（1）：95-99.

［3］MCGLINCHEY D，COWELL A，KNIGHT E A，et al.Bend pressure drop predictions using the Euler-Euler model in dense phase pneumatic conveying［J］.Particulate Science and Technology，2007，25（6）：495-506.

［4］MA A C，WILLIAMS K C，ZHOU J M，etal.Numerical study on pressure prediction and its main influence factors in pneumatic conveyors［J］.Chemical Engineering Science，2010，65（23）：6247-6258.

［5］VASHISTH S，GRACE J R.Simulation of granular transport of Geldart type-A，-B，and -D particles through a 90°elbow［J］.Industrial ＆ Engineering Chemistry Research，2012，51（4）：2030-2047.

［6］GIDASPOW D，JUNG J，SINGH R K.Hydrodynamics of fluidization using kinetic：an emerging paradigm 2002flour-daniel lecture［J］.Powder Technology，2004，148（2）：123-141.

［7］LU H L，GIDASPOW D，BOUILLARD J，etal.Hydrodynamic simulation of gas-solid flow in a riser using kinetic theory of granular flow［J］.Chemical Engineering Journal，2003，95（1／2／3）：1-13.

［8］PU W H，ZHAO C S，XIONG Y Q，etal.Threedimensional numerical simulation of dense pneumatic conveying of pulverized coal in a vertical pipe at high pressure［J］.Chemical Engineering ＆ Technology，2008，31（2）：215-223.

［9］JENKINS J T，SAVAGE S B.A theory for the rapid flow of identical，smooth，nearly，elastic spherical particles［J］.J Fluid Mech，1983，130（5）：187-202.

［10］JOHNSON P C，JACKSON R.Frictional-collisional constitutive relations for granular materials，with application to plane shearing［J］.J Fluid Mech，1987，176（9）：67-93.