烟气绕流换热管流动特性的无网格法数值分析

2018-10-09闫顺林王皓轩

发电设备 2018年5期

闫顺林, 王皓轩, 韩韦

(华北电力大学能源动力与机械工程学院，河北保定 071003)

换热管作为一种换热元件，被广泛应用于锅炉水冷壁、省煤器、凝汽器等换热设备中。由于烟气携带大量的飞灰颗粒，且烟气流动具有很大的随机性，使得飞灰对换热管的冲击碰撞点也相对随机，这对于换热管的防磨带来一定的困难，换热管的磨损将影响到电厂设备运行的安全性和经济性[1-2]。

虽然国内外存在大量的换热管的磨损模拟研究，也提出很多相应的防磨措施，但是大量的模拟研究基本都是把烟气正向冲刷作为入口条件进行二维流动模拟，不能有效地反映出实际情况[3-7]。因此，研究三维下涡结构的发展及运动特性，更能反映出流动的真实情况[8-9]。由于烟气流动中具有很大的随机性，多数情况下烟气冲刷换热管都具有一定的攻击角度，即便安装梳行管、导流板来使烟气流动均匀，但仍然不能保证烟气正向冲刷换热管。

因此，笔者采用一种基于格子Boltzmann方法(LBM)的计算流体仿真软件——XFlow，在三维空间中研究不同攻击角度下烟气的流动情况、涡结构的演化规律，以及换热管的静态气动力特性等，以期为换热管的防磨措施提供一些指导性的建议。

1 物理模型和边界条件

图1为简化的烟气绕流换热管模型。为了使模拟的结果更接近实际情况，特选取一个16 m×8 m×8 m的长方体作为烟气通道，直径D=50 mm的单个换热管光管作为研究对象,光管位置选择管束中前排管子，因为前排管束受烟气冲刷的直接作用较强。

图1 烟气绕流换热管简化模型

左侧设定为烟气速度入口边界，v=9 m/s,且入口的湍流强度设定为5%，右侧为烟气自由出口边界，其余四面设置为无滑移壁面边界。烟气温度设为800 K，动力黏度设为3.328×10-5Pa·s，烟气质量流量定为1.5 kg/m3，不考虑重力因素的影响。

2 数值模型及计算方法

2.1 计算方法

XFlow最大特点就是具备无网格方法，即不需要对传统的计算域进行网格划分，突破了模拟实体结构表面的复杂性限制，避免了因网格质量不高而带来的误差，从而使模拟结果更符合实际情况。

LBM源于格子气自动机(LGA)方法[10-11]。作为一类新的流体模型和计算方法，LGA方法也存在不足之处：统计噪声、碰撞算子的指数复杂性，不满足伽利略不变性等。为了克服LGA的这些不足，LBM方法开始出现并发展。

(1)

(2)

(3)

(4)

(5)

然后在对其在时间和空间上进行离散，从而得到完全离散、含外力项的格子LBGK，即

(6)

式中：δt为时间步长；fa为连续分布函数f的离散形式(a=0,1,2,…,N)；faeq为离散速度空间的局部平衡态分布函数；Fa为离散速度空间的外力项；τ=τ0/δt为无量纲弛豫时间。由于在求解过程中沿特征线进行积分，式(6)具有Lagrange特征，使模型能够达到二阶数值计算精度。

2.2 格子模型

DdQm(d维空间，m个离散速度)系列模型是LBM的基本模型，常用的模型参数有D1Q3、D2Q7、D2Q9、D3Q19、D3Q27等。

格子Boltzmann方程(LBE)可以通过把Boltzmann方程的速度空间离散成有限数量的离散速度得到。虽然有许多关于三维流动模拟离散速度方案，但考虑湍流流动数值计算的稳定性和精度，笔者在模拟烟气绕流换热管的三维流动中，采用的是基于多松弛格子Boltzmann方法(MRT-LBE)[12]下的D3Q27模型[13-14]。

该模型分布函数的碰撞过程不是在速度空间完成，而是通过线性变换，在矩空间完成。这种碰撞形式能够提供一个自然、方便的途径去表示碰撞中不同时间尺度物理量的松弛过程，从而避免单松弛模型中所有物理量均是在同一时间尺度下向其平衡态弛豫，其粒子分布函数f的演化方程可表示为:

(7)

(8)

(9)

表1 D3Q27离散速度模型参数

该模型的变换矩阵为：

(10)

式中：s4=1.54;s5=s7=0.62;s10=1.50;s13=1.83;s16=1.40;s17=1.61;s18=s20=1.98;s23=s25=1.74。

运动黏度及体积黏度系数分别为：

(11)

(12)

2.3 计算采用的湍流模型

湍流模型的数值求法主要有直接数值模拟(DNS)、雷诺平均法(RANS)和大涡模拟(LES)。DNS方法直接采用瞬时的N-S方程对湍流进行数值求解计算，计算量大，对计算机本身的要求也很高，目前只可适用于低雷诺数范围(见图2)。RANS是一种基于整体全局观念模拟湍流的方法，也是应用最广的方法，RANS计算变量的时间平均，从而消除求解上时间独立性问题，虽然能够降低计算成本，但是需要在N-S方程上添加一些新的输运方程项。目前已经有很多的RANS模型，但都是针对特殊问题的模型，并且模型中的许多参数需要依靠人工经验设置。LES是解决湍流的一种折中方法，它是基于局部的数值方法，对微观小尺度进行模拟，更加接近真实物理。

图2 湍流模型示意图

LES分析研究重点在小尺度上的湍流模型，不需要任何主观参数来刻画湍流现象，XFlow就是基于LES方法来模拟湍流流动，在LES模型里，引入一个额外的黏度，称为湍流涡流黏度υt，模拟湍流流动，即

υeffective=υt+υturbulent

(13)

式中：υeffective为有效湍流黏度；υturbulent为湍流黏度。

另外，对于边界层的模拟，XFlow采用的是Wall-Modeled Large Eddy Simulation (WMLES)方法[15]，不需要在尾涡外的剪切区添加人工湍流黏度，Wall模型实施如下：

(14)

(15)

3 模拟结果及分析

笔者分别模拟计算了入口烟气攻击角度θ=0°、10°、20°、45°、60°五个工况，θ设定见图3。分析讨论了五种工况下换热管静气动力特性的变化情况，计算域涡结构的演化规律及其影响。

图3 烟气冲刷管壁示意图

3.1 计算域拓扑结构尺寸的设定

为了使模拟结果更加接近实际的流动情况，传统的CFD软件往往会在网格结构和数量上进行加密和修改，通过模拟结果与以往文献或是实验结果相比较，再加上计算时长的因素，来选取一个相对较好的网格结构和数目，但仍然避免不了在网格的画取上，人为因素所导致较大误差的产生。而在XFlow中，因为它是基于LBM的空间域离散化的域结构，可以根据设定的计算尺寸自动生成，一定程度上减少了误差的产生；同时在计算运行中又具备格子自适应细化功能，可以根据速度场的速度梯度变化程度，自动对速度梯度变化剧烈的位置进行格子加密，确保准确有效地计算出流场的分布情况。

分别取计算运行中0.05 s和0.2 s两个时刻点对应的域结构形式，从图4中可以看到：随着烟气绕流换热管尾流区的增大，计算域内的域结构发生了很大变化，尤其在换热管后面速度梯度变化较为剧烈的区域，格子密度和数量明显增加，并且在设置该计算尺寸的情况下，运行过程中的稳定性参数也基本维持在0.3处(稳定性参数<1,就说明计算域结构满足计算需求)，说明了计算设置的格子尺寸满足了计算的需要，保证了模拟结果的正确性。

图4 空间离散化的域结构图

3.2 气动力特性的变化情况

笔者就烟气流动冲刷换热管角度的不确定性，模拟计算在不同的θ下烟气绕流换热管的气动力特性随时间的变化情况(θ为五个值情况下的三分力系数及阻力矩变化见图5)。

从图5可以看到：几种情况下的换热管的阻力系数波动趋势具有一致性，且存在短时间内阻力系数为负值的情况，但是时间平均的阻力系数不会出现明显的负值；另外，较大攻击角度的负值出现得较早，而且其升力系数也相对较高，扭转系数波动也相对剧烈，说明大冲击角度下烟气流动分离加剧，非定常效应加强，导致某时刻出现阻力减小，甚至为负，提供推力的现象发生，而0°冲击角下的升力系数接近于0，扭转系数波动较为缓慢，阻力矩的绝对值也是最小，接近0，这都有利于流动均匀。

3.3 涡结构的演化规律及影响

为了更清晰地看到不同攻击角度下流体域内涡结构的发展情况，就攻击角度为0°、45°、60°，图6给出了这三个情况下不同时刻对应的涡结构3D云图。

图6 流动空间涡结构发展3D云图

从图6可以看出：流场的发展速度之快，形成多个在三维空间内相互交错的发夹涡等涡结构，体现出了湍流的混浊性。从涡结构的演化中看到，随着θ的增大，湍流发展越迅速，涡结构发生倾斜运动的时间也越早。θ=0°时换热管尾迹区域内湍流发展相对较慢，涡运动也基本保持在水平位置，较长距离的涡运动途径也表现出了涡能量的耗散过程；θ=60°时发生湍流转捩的时间较短，而且后期的湍流扰动也最剧烈。

另外，从图6中也可以看到，随着尾流域内离换热管距离的增大，涡结构的能量逐渐耗散，甚至消失。图7表示部分攻击角度下的瞬时流动涡矢量图，发现低攻击角度下在尾流区远场区域，出现两个或是三个涡的涡对现象，在两个涡的涡对中，两个涡的旋度相反，涡量大小却基本相等。三个涡的涡对中，其中一个涡与另外两个涡的涡量大小大致相等，但方向相反，而较大攻击角度下却很难形成远场二次涡对，但却在计算域内形成大范围的涡结构，其主要是侧壁面存在和换热管后形成的涡街共同作用的结果，侧壁面的存在又阻碍了气流的旋转，制约了涡量的输运，由于壁面的剪切作用也消耗掉部分旋转动量，使涡量损失较大，阻碍气流流动。另外，大范围的涡结构的存在，会出现暂时性的气流逆向冲刷换热管表面的现象，这也正解释了图5中换热管阻力系数短时间内出现负值的原因。