周海军，熊源泉

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

引 言

就目前形势来看，在未来相当长的时间内，煤炭仍将作为我国最主要的一次能源[1]。因此选择清洁高效的煤炭技术对其进行合理开发和利用势在必行。而煤气化作为煤炭清洁高效利用的核心技术，正受到越来越多的关注。尤其是大规模高效气流床煤气化技术，它是目前最具前景的煤气化技术，而高压密相气力输送正是其关键技术之一[2]。

与其他输送方式相比，高压密相气力输送具有能耗低、耗气量低、管道磨损小以及固相浓度高等优点，因此近年来被广泛应用于能源、化工以及冶金等行业。但由于它是一种非线性动态响应的气固两相流，而且具有固相浓度高、影响因素多以及流动形态复杂等特点[3-4]，导致目前仍未完全弄清其输送特性。

关于高压密相气力输送，通常可采用试验和数值模拟两种方法对其进行研究。但受其复杂的输送特性以及现有测量技术的限制，仅依靠试验研究已很难再有进一步的突破。然而通过数值模拟方法却可以捕捉计算区域内每一个单元格的详细流动信息，从而有助于揭示其输送特性[5]；而且还能节省试验成本。所以数值模拟方法是探讨高压密相气力输送非常必要且有效的研究手段。但是目前高压密相气力输送的相关研究主要集中在试验研究领域，数值模拟的研究仍然是非常有限的。其中，Yuan 等[4-5]采用欧拉-拉格朗日方法对高压密相气力输送做了一系列的模拟研究。但该方法不仅计算量大、周期长而且消耗计算资源多，所以并不适于模拟高压密相气力输送。而谢灼利等[6]则通过引入颗粒动理学理论在双欧拉方法的基础上对水平管气力输送进行了模拟研究，但是他们却忽略了摩擦应力的作用。所以蒲文灏等[7-8]在此基础上通过修正Johnson-Jackson 摩擦压应力模型并结合Syamlal摩擦黏度模型引入了摩擦应力模型，同时也因此而获得了较为完整的三维非稳态数理模型；而且他们采用该数理模型对水平管高压密相气力输送进行了模拟研究，并取得了一些具有价值的研究成果。但是该数理模型依然存在一些理论上的缺陷与不足。一方面，颗粒动理学理论是类比于稠密气体分子运动学理论而构建的一种固相应力封闭理论，该理论认为颗粒在气固两相流系统中呈均匀分布，而且颗粒间作用为随机的双颗粒瞬时碰撞[9-10]。所以它通常用于表征稀相流的颗粒间碰撞作用[11]。而摩擦应力模型则是基于不同屈服准则而构建的固相应力模型，主要用于表征颗粒处于持续接触时的挤压和摩擦特性。只有当固相体积浓度较高时，颗粒间才会形成挤压并产生摩擦应力。因此摩擦应力模型适宜描述密相流中的固相应力。同时当颗粒间作用由随机的双颗粒瞬时碰撞逐渐转变为颗粒持续接触引起的挤压和摩擦作用时，其间必然存在一个过渡过程；而且很多学者[12-14]已然证实了这一点：在稀相流和密相流之间存在一个过渡流。在过渡流中，颗粒间作用主要表现为相关的多颗粒碰撞[15]。但蒲文灏等认为摩擦应力的临界固相体积浓度为0.1，导致其提前引入了摩擦应力，从而高估了摩擦应力的作用。另一方面，Johnson-Jackson 摩擦压应力模型是基于粒径为1 mm 的玻璃珠的物性而获得的，虽然蒲文灏等根据模拟结果对该模型中的常数系数做了一定的修正，但却未给出相关的理论依据。

水平管高压密相气力输送本质上是一种典型的稠密气固两相流，在气体-颗粒、颗粒-颗粒以及颗粒-壁面间的相互作用下，通常处于非线性非平衡状态。所以在流动过程中颗粒存在多尺度结构：微尺度结构（分散的单颗粒）和介尺度结构（颗粒团），导致流动呈现强烈的非均匀结构特征。其中，介尺度结构的存在是其呈现非均匀结构特征的根本原因，同时也会导致气固滑移速度增大，曳力下降[16]。因此在水平管高压密相气力输送的数值模拟中需要考虑介尺度结构效应。然而先前的数理模型采用均是Gidaspow 曳力模型[6,8,11]，该曳力模型是基于均匀化的气固两相流而构建的[17]，未考虑介尺度结构引起的曳力下降，所以不能准确地描述水平管高压密相气力输送中的非均匀结构输送特性。

补充风是高压密相气力输送试验中一个非常重要的操作参数，它是由缓冲罐直接引入到输送管道中的一股风量，所以并不会影响煤粉在发送罐内的流化状态。但是由于缓冲罐的压力基本上是整个输送系统中最高的，所以补充风势必会改变输送管道中的压力分布，进而影响煤粉的出料驱动力，导致煤粉的出料量/质量流量发生改变。同时，补充风也会增加输送管道中的输送风量，从而提高输送表观气速，调节输送固气比[18]，进而改变输送管道中的流动形态，导致输送机理也发生相应的变化。所以考察和研究补充风对水平管高压密相气力输送的影响具有十分重要的现实意义。

本文以水平管高压密相气力输送为模拟研究对象，针对其数理模型中存在的缺陷与不足，做进一步的改进和完善。并采用改进后的数理模型对一组仅改变补充风的高压密相气力输送试验工况进行模拟计算和分析，考察补充风对水平管高压密相气力输送的影响机制，从而揭示其输送特性，为高压密相气力输送的合理设计与全面优化提供更可靠的理论基础和经验指导。

1 输送试验及试验装置

图1为东南大学自主研发的上出料式高压密相气力输送试验装置，整个输送试验装置主要由供气系统、储料罐系统、输送管道及其附件、测量传感器以及数据采集与控制系统五部分构成。在输送试验中，两个储料罐分别用作发送罐（保持其压力为3.0 MPa）和接收罐（保持其压力为2.5 MPa）。高压气瓶中的氮气流经缓冲罐后被分成充压风、流化风和补充风三路风量。其中，流化风（保持其流量为0.4 m3/h）由发送罐底部布风板引入用于流化罐内煤粉，充压风由发送罐上部引入用于维持罐内压力恒定，补充风则由输送管道前端引入用于调节输送管道内气固比。在两罐压差驱动下，被流化的煤粉由发送罐流出，经管径为10 mm 的输送管道输送后流入接收罐。在接收罐内含煤粉的气流经布袋除尘器除尘后，采用电动排气阀放空以维持接收罐的压力。煤粉输送完毕后，切换输送试验装置中的相关阀门，使发送罐和接收罐互换，进入下一组输送实验。本文的试验目的是考察补充风对水平管高压密相气力输送的影响，因此在输送试验中仅改变补充风流量，试验工况详细参数见表1。输送物料为内蒙褐煤煤粉，其主要物性参数见表2。关于输送试验更详尽的介绍可见文献[19-20]。

图1 高压密相气力输送试验装置Fig.1 Experimental facility schematic diagram of dense-phase pneumatic conveying under high pressure

表1 输送试验工况参数Table 1 Conveying experiment parameters

表2 内蒙褐煤煤粉的主要物性参数Table 2 Main physical properties of pulverized lignite

2 数理模型

2.1 控制方程

连续性方程

式中，αg、αs分别为气固两相体积浓度；ρg、ρs分别为气固两相密度；vg、vs分别为气固两相速度。

动量方程

式中，g为重力加速度，m/s2；Pg为气相压力，Pa；σg、σs分别是气固两相应力张量,Pa；Fsg为气固两相间曳力，N；β为曳力系数。

2.2 固相应力模型

如引言所述，水平管高压密相气力输送实际上是由三种流动形态稀相流、过渡流以及密相流共同构成的非均匀结构稠密气固两相流。对于稀相流，可用颗粒动理学理论对其固相应力进行描述，该理论公式见文献[7]。但是考虑到颗粒动理学理论在描述过渡流和密相流时存在一定的局限性，所以本文引入Savage 径向分布函数g0,ss[式（6）]对其进行了修正，从而使其能更精准地描述密相流中的碰撞作用[21]。而且Lee 等[22]也验证了该方法的可靠性和适用性。

式中，αs,max是颗粒堆积时的最大固相体积浓度。

对于密相流，可用摩擦应力模型描述其摩擦应力，并结合颗粒动理学理论表征其固相应力[23]。然而Johnson-Jackson 摩擦压应力模型引入了三个经验常数来表征材料的摩擦特性，却又未提供确定这三个经验常数的方法；而且该模型未考虑任何相关的物性参数[24]。所以其经验性特别强，而其准确性以及适用性却相对较差。从力学机制上来看，摩擦应力与颗粒刚度kn以及颗粒平均粒径ds是相关的，所以Berzi 等[25]引入这两个物性参数构建了Berzi 摩擦压应力模型[式（7）]。同时本文又通过耦合Pitman-Schaeffer-Gray-Stiles 屈服准则构建了与Berzi 摩擦压应力模型相对应的摩擦黏度模型[式（8）][26]，与Berzi 摩擦压应力模型共同构建了一个新的摩擦应力模型。该摩擦应力模型不仅考虑了材料的相关物性，而且同时兼顾了密相流中颗粒相的膨胀性和压缩性，从而克服了Johnson-Jackson 摩擦压应力模型和Syamlal摩擦黏度模型存在的缺陷。

式中，μf、pf分别是摩擦黏度和摩擦压应力；a是表征材料摩擦特性的经验常数；ϕi是内摩擦角，（°）；Θs是颗粒拟温度，m2/s2；αs,min是摩擦应力的临界固相体积浓度。考虑到在自然堆积状态下颗粒自重引起的静摩擦作用，所以本文设αs,min=0.4。而且Schneiderbauer 等[27]也认为当固相体积浓度大于0.4时，颗粒间开始产生挤压和摩擦作用。

2.3 曳力模型

气力输送实际上就是以气体为载体并借助其压能和动能将粉体颗粒输送至指定位置的过程。在这个过程中，固相颗粒通过与气相间的相互作用从而获得能量并跟随气体运动至指定位置。因此气固两相间的作用力是其进行动量交换和能量传递的枢纽，也是气力输送中最重要的作用力之一，它主要包括曳力、升力、浮力、Basset 力以Saffman 力等。但在水平管高压密相气力输送中，其他相间作用相对较小，而曳力却是最主要的相间作用力，所以通常只需考虑曳力[8]。

图2 为当Qs=0.8 m3/h 时采用电容层析成像技术（ECT）测得的水平管纵截面固相体积浓度分布。该图进一步证实了水平管高压密相气力输送是多种流动形态（稀相流、过渡流以及密相流）并存的非均匀结构气固两相流。在稀相流中曳力占主导地位，固相颗粒通常以单颗粒的形式处于悬浮状态，所以流动呈均匀分布，介尺度结构一般较难形成。而在密相流中固相应力占主导地位，从而限制了颗粒团的形成，所以介尺度结构难以维持和稳定，流动依然呈均匀分布[16]。但在过渡流中，曳力、固相应力以及重力呈协调竞争关系，所以在过渡流中始终进行着激烈的动量、质量以及能量交换[7]，同时也必然引起颗粒团聚的形成，从而使过渡流呈现强烈的非均匀结构特征。因此介尺度结构通常主要存在于过渡流中。

图2 水平管纵截面固相体积浓度分布Fig.2 Solids volume fraction distribution at longitudinal section of horizontal pipe

基于气固两相流中存在多种流动形态的考虑，一些研究者通过引入多段式曳力模型对各种流动形态的流动特性进行了模拟研究。其中，lyu 等[28]通过引入Mckeen 和Gidaspow 曳力模型构建了三段式曳力模型，并采用该曳力模型对稠密气固两相流中稀相区、过渡区以及密相区的流动特性进行了模拟研究，从而验证了三段式曳力模型的可靠性。其中，Mckeen 曳力模型采用一个小于1 的尺度因子对Gibilaro 曳力模型进行修正以考虑介尺度结构引起的曳力下降[29]，所以该曳力模型适用于过渡流的模拟。而Gidaspow 曳力模型则是基于平均化方法而获得的典型曳力模型，它可用于稀相流以及密相流的模拟。因此本文对lyu 等构建的三段式曳力模型进行了适当的修正，从而构建了一个能同时兼顾水平管高压密相气力输送中三种流动形态的曳力模型。

修正的三段式曳力模型为

式中，Ψ1和Ψ2均是反正切转换函数，用于连接各段曳力系数，使其在断点处光滑过渡，从而改善模拟计算的收敛性和稳定性；μg是气体黏度；CD是阻力系数；Res是颗粒Reynolds数。

2.4 气固湍流模型

为考察水平管高压密相气力输送中的气固两相湍流效应，本文引入了realizablek-ε模型分别对气固两相湍流（realizablek-εmodel for per phase）进行模拟和分析。

3 数理模型

3.1 气固湍流模型

（1）进口边界条件 进口边界条件设为速度进口边界条件。对于气相,设径向及切向气速为零，轴向气速分布vg,inlet(r)为

式中，r为到进口截面圆心的距离,m；D是管道直径,m。

对于固相,根据经验设置相同的气固两相速度分布模拟更易收敛，故其轴向速度分布vs,inlet(r)为

其中，式（20）是用于计算us,inlet的经验公式[30]。

（2）出口边界条件 出口边界条件设为压力出口边界条件，出口压力设为Pout。

（3）湍流设置 设水力直径为管道直径，气固两相湍动强度分别为10%和5%。

（4）壁面边界条件 对于气相，采用无滑移壁面边界条件。对于固相，引入了考虑摩擦应力项的Johnson-Jackson壁面边界模型[31]

式中，τsw和qw分别是颗粒-壁面间剪切应力以及拟热流；ϕ为镜面系数；vsw是颗粒-壁面间的滑移速度,m/s；esw是颗粒-壁面间碰撞恢复系数；μw是颗粒-壁面间的摩擦系数。

3.2 模拟设置

在三维双精度算法的基础上，采用本文改进后的水平管高压密相气力输送数理模型对表1中的输送试验工况进行模拟计算。在计算中，设置时间步长为5×10-5s，单个时间步长内迭代次数设为40，收敛残差设为10-4。其他相关的数理模型参数设置见表3。

需要明确指出的是本文所用颗粒径向分布函数、摩擦应力模型、速度分布函数及曳力模型等均采用FLUENT 提供的UDF 进行编译，并加载到软件ANSYS FLUENT 15中进行模拟计算。

表3 数理模型参数Table 3 Mathematical model parameters

4 结果与讨论

4.1 网格划分及其无关性分析

模拟对象为长2.4 m 的水平管，该水平管位于图1 中的A 点位置。采用GAMBIT 软件对水平管进行网格划分，如图3所示。为了保障网格质量，进口端面采用铜钱型网格且呈非均匀分布，其中壁面附近采用相对较密的网格以确保其计算效果。

图3 水平管的网格划分Fig.3 Meshes generation of horizontal pipe

为了尽量提高数值模拟的计算精度和计算效率，采用表4中Mesh A、B、C、D四种不同规格的网格对补充风Qs=0.8 m3/h 的典型试验工况进行模拟计算，考察网格无关的界限并确定合适的网格尺寸。

由表4可以看出，随着网格尺寸的减小，模拟预测的水平管压降不断接近试验值，尤其是当管道总网格数大于46.08 万个（Mesh C 和Mesh D）时模拟预测的水平管压降几乎不再变化。图4为不同网格尺寸下模拟预测的气相速度沿高度方向的分布，由该图发现：随着网格尺寸的减小，气相速度分布逐渐趋于一致。尤其在Mesh C 时已基本与网格尺寸无关，因此本文的数值模拟计算选用Mesh C。

4.2 压降的对比和分析

图5为不同补充风下模拟预测的水平管压降及其试验值的对比关系，通过该图计算发现：与水平管压降的试验数据相比，模拟结果的相对误差为-2.7%～+4.1%；所以模拟结果与试验数据是相当吻合的。而且模拟还合理地预测了水平管压降随补充风的变化规律：补充风增加，水平管压降先减小后增加。从而验证了本文所用数理模型的可靠性及适用性。

表4 不同网格尺寸下模拟预测的水平管压降Table 4 Predicted pressure drop of horizontal pipe with different grid scale

图4 不同网格尺寸下模拟预测的气相速度沿高度方向的分布Fig.4 Predicted gas velocities vs dimensionless height with different grid scale

图5 不同补充风下模拟预测的水平管压降与其试验值的对比Fig.5 Comparison of predicted pressure drop of horizontal pipe with its experimental data at different supplementary gas flow rates

图6 水平管中不同横截面气固两相速度沿高度方向的分布Fig.6 Velocities of gas and solids phase vs dimensionless height at different cross sections of horizontal pipe

4.3 固相体积浓度分布的对比和分析

图6 为当Qs=0.8 m3/h 时水平管不同横截面的气固两相速度沿高度方向（竖直方向）的分布。由该图可以发现：由水平管进口沿流动方向（轴向方向），气固两相速度很快便达到稳定状态，在距离进口L=120D的截面后，气固两相速度已经几乎不再随流动距离的变化而变化，说明水平管中的气固两相流已经进入了充分发展段。所以本文选择距离进口L=180D截面处的力学特性参数以及输送特性参数进行详细讨论和分析，探究水平管高压密相气力输送特性。

图7 水平管横截面固相体积浓度分布云图Fig.7 Solids volume fraction distribution at cross section of horizontal pipe

图8 内蒙褐煤煤粉颗粒粒径分布Fig.8 Particle size distribution of pulverized lignite

图7 为当Qs=0.8 m3/h 时水平管横截面固相体积浓度分布云图，由该图可以看出：无论是根据ECT技术测得的横截面固相体积浓度分布云图（为方便描述下文一律简称为ECT 图）还是根据模拟预测的横截面固相体积浓度分布云图，水平管中的固相体积浓度分布均可划分为三个区域：对应于稀相流的上部悬浮区、对应于过渡流的中部过渡区以及对应于密相流的底部沉积区。而且各区域间的分界线并不是十分明显，因为过渡区与其他两个区域始终进行着激烈的质量、动量以及能量交换。另外，对比图7（a）、（b）两图还可以发现：ECT 图中的固相体积浓度随高度方向的变化更为平缓，这是由于在数值模拟中未考虑颗粒粒径分布（图8）而引起的。

图9为不同补充风下模拟预测的固相体积浓度分布云图与ECT 图的对比关系，通过对比该图中9（a）、（b）两图可以发现：模拟预测的固相体积浓度分布云图与ECT 图总体上是相吻合的。而且随着补充风的增加，模拟结果与ECT 图均获得了相同的变化规律：水平管上部悬浮区不断增大，而其底部沉积区却不断减小。从而进一步证实了数理模型的可靠性及适用性。

4.4 补充风对水平管高压密相气力输送特性以及力学特性的影响

图10 为不同补充风下模拟预测的输送特性参数沿高度方向的分布。由表1 可知：表观气速随补充风的增加而增大，因此模拟预测的气相速度亦随补充风的增加而增大[图10（a）]，从而导致气固两相间动量交换和能量传递变得更加激烈，其表现为曳力随着补充风的增加而增大（图11），所以固相速度也会随补充风的增加而增大[图10（b）]。同时曳力增大也会导致气固两相脉动速度以及相关的扰动作用增强，所以气固两相湍动能以及颗粒拟温度均随着补充风的增加而增大[图10（d）、（e）、（f）]。

图9 不同补充风下模拟预测的固相体积浓度分布云图与ECT图的对比Fig.9 Comparison of predicted solids volume fraction distribution with ECT image at different supplementary gas flow rates

图10 不同补充风下模拟预测的水平高压密相气力输送特性参数沿高度方向的分布Fig.10 Predicted conveying parameters of horizontal pipe vs dimensionless height at different supplementary gas flow rates

图11 不同补充风下模拟获得的曳力FsgFig.11 Predicted drag forces at cross section of horizontal pipe at different supplementary gas flow rates

图12 不同补充风下模拟获得的固相摩擦切应力分布云图Fig.12 Predicted solids frictional shear stress at cross section of horizontal pipe at different supplementary gas flow rates

由表1还可以发现：随着补充风的增加，煤粉的质量流量不断减小；而固相速度却不断增大[图10（b）]，所以固相体积浓度不断降低。而曳力，尤其是其竖直方向分量的增大[图11（a）]，会带动更多的煤粉颗粒处于悬浮状态。所以随着补充风的增加，悬浮区会不断增大，而沉积区却不断减小（图9）。从而导致水平管底部的固相摩擦切应力也随着补充风的增加而不断降低（图12），而且颗粒-壁面间的切应力也呈现出同样的变化规律[图13（b）]，因此颗粒-颗粒以及颗粒-壁面间的固相摩擦压损会随着补充风的增加而减小。相反，气体-壁面间的切应力却随着补充风的增加而增大[图13（a）]，这是因为气相压损主要取决于气速，而且通常与表观气速的2～3 次方呈正比[32]，而表观气速是随着补充风的增加而增大的（表1）。因此，气体-壁面的气相摩擦压损随着补充风的增加而增大。

对比图13（a）、（b）两图可以看出：气固两相与壁面的切应力始终处于竞争关系，当Qs＜0.8 m3/h时，颗粒-壁面间的切应力大于气体-壁面间的切应力并占主导地位。而当Qs≥0.8 m3/h 时，气体-壁面间的切应力开始反超颗粒-壁面间的切应力以及固相摩擦应力而占主导作用。所以水平管压降随补充风的增加而呈先减小后增加的规律（图5）。

图13 不同补充风下模拟获得的气固两相与壁面间的切应力分布云图Fig.13 Predicted frictional shear stresses between gas,solids and wall at different supplementary gas flow rates

5 结 论

本文以水平管高压密相气力输送为研究对象，在欧拉-欧拉方法的基础上，引入Savage 径向分布函数修正颗粒动理学理论以考虑稀颗粒间的碰撞作用，并采用基于Berzi 摩擦压应力模型构建的新摩擦应力模型以考虑摩擦作用，同时耦合修正的三段式曳力模型以考虑各种流动形态中的相间作用力，从而获得了一个改进的三维非稳态数理模型。并采用该数理模型考察了补充风对水平管高压密相气力输送的影响，获得如下主要结论。

（1）模拟精准地预测了水平管压降及其随补充风的变化规律，模拟结果的相对误差为-2.7%～+4.1%。

（2）模拟预测的水平管固相体积浓度分布与ECT图总体上是相吻合的。而且试验和模拟结果均表明：其固相体积浓度分布可划分为三个区域：对应于稀相流的上部悬浮区、对应于过渡流的中部过渡区以及对应于密相流的底部沉积区。随着补充风的增加，悬浮区不断增大，而沉积区却不断减小。

（3）随着补充风的增加，曳力以及气体-壁面间的切应力增大，而固相摩擦切应力以及颗粒-壁面间的切应力减小。

（4）随着补充风的增加，气固两相速度和湍动能以及颗粒拟温度增大，而固相体积浓度减小。

符 号 说 明

ess——颗粒-颗粒间碰撞恢复系数

h——颗粒在水平管中的高度位置，m

k——湍动能，J

L——与管道进口位置的距离，m

pf——固相摩擦压应力，Pa

qw——壁面拟热流，w

u——速度，m/s

v——速度矢量，m/s

α——体积分数

μ——动力黏度，Pa·s

μf——固相摩擦黏度，Pa·s

ρ——密度，kg/m3

τsw——颗粒-壁面间的切应力

下角标

g——气相

s——固相