非均匀加热管内超临界CO2传热特性研究

2019-12-04吴家荣李红智韩万龙

热力发电 2019年11期

吴家荣，李红智，杨玉，韩万龙

（西安热工研究院有限公司，陕西西安 710054）

超临界CO2布雷顿循环(SCO2-BC)因布局简单灵活，循环效率高，设备紧凑，工质无毒、廉价等优势在核电、煤电、余热回收和舰船推进等多个领域得到了越来越多的研究[1-4]。作为SCO2-BC 发电系统的关键环节，超临界CO2在锅炉气冷壁内的传热流动规律关系着整个系统的安全和效率。不同于亚临界流体，超临界CO2的物性随着压力和温度的变化呈现强烈的非线性变化，尤其在拟临界点附近，流体的比定压热容和导热系数会发生明显突变，如图1所示。物性变化及其引起的自然对流使流动状态由单一的强制对流发展为复杂的混合对流，伴随而来的是传热恶化和强化的发生，这给锅炉气冷壁的热力设计造成了困难。

目前，比定压热容和导热系数的增大被认为是传热强化的主要原因；而关于传热恶化，一种观点是Ackerman 等人[5]根据超临界水传热提出的拟膜态沸腾理论，另一种则是Jackson 等人[6-8]提出的传热恶化与浮力效应和加速效应对湍流的产生和传热有关。

图1 CO2 物性Fig.1 Physical properties of CO2

针对上述观点，各国学者开展了大量超临界流体的传热研究。但是，这些研究[9-12]绝大多数都是在均匀热流条件下开展的，这与气冷壁在炉膛内单侧受热的实际情况不符。因此，也有学者开始重点研究超临界流体在非均匀热流条件下的流动传热。Bai 等人[13]比较了超临界水在均匀和非均匀加热圆管内不同母线处的壁温、温差、对流传热系数随主流焓的变化，模拟结果表明，非均匀加热时，圆管周向壁温的分布是不均匀的，传热强化也仅发生在局部区域，某些母线处甚至出现流体温度高于壁温的情况。Fan 等人[14]利用数值模拟对高质量流速下超临界CO2在竖直非均匀加热圆管中的传热特性进行了研究，结果表明，热流密度很大的顶部母线处浮力和加速效应对传热的影响远大于热流密度较小的底部母线处。强化传热的主要原因是拟临界区比定压热容的增大；而传热恶化除受浮力和加速效应的影响外，还源于管内黏性底层的增厚。Gu 等人[15]比较了超临界水在内螺纹管内均匀和非均匀受热的情况，数值模拟结果显示，非均匀受热时，因管内旋流的影响，流体混合更好，使得管内最高内壁温比均匀受热时的最高内壁温低了10 K，且最高内壁温并不总出现在受热面的中点。上述对超临界流体在非均匀热流条件下的研究都集中在内径＞ 10 mm、质量流量较大的情况，对小管径情况下超临界CO2传热特性的研究还存在不足。

综上所述，本文开展了超临界CO2在非均匀热流密度下圆管内传热的数值模拟研究。所选取的参数为：质量流速100 kg/(m2·s)、热流密度28～50 kW/m2、压力7.6～8.6 MPa。通过在加热段一侧施加热流以模拟气冷壁受热情况，比较了均匀和非均匀加热的差别，分析了非均匀加热情况下热流密度和入口压力对对流传热系数、浮力效应和加速效应的影响，总结并检验了4 种常用传热关联式对本文数据的适用性，最后提出并检验了新的传热关联式。

1 物理模型和控制方程

图2为加热圆管的几何模型。图中竖直圆管外径为8 mm，内径为5 mm，加热段长660 mm，进出口段长110 mm。在圆管外壁的一侧施加热流，由于管壁周向和轴向导热的存在，整个圆管的热流密度是非均匀的。

图2 加热圆管几何模型Fig.2 Geometry model of the heated circular pipe

笛卡尔坐标系下可压缩流体的连续性方程为

式中：ρ为密度，kg/m3；ui为速度分量，m/s。

动量方程为

式中：gi为重力加速度分量，m/s2；为黏性系数，kg/(m·s)；uj、uk为速度分量，m/s；u′为速度脉动量，m/s。

能量方程为

式中：λeff为有效导热系数，W/(m·K)；h为焓，J/kg；cp为比热容，J/(kg·K)。

湍动能方程为

式中：μt为湍流黏性系数，kg/(m·s)；k为湍动能，J/kg；σk、Gk、Yk为模型系数。

耗散率方程为

式中：ω为耗散率；σω、Gω、Yω、Dω为模型系数。

式(1)—(5)中，上标“—”表示时间平均项，上标“～”表示Favre 平均项，各项具体的表达式可参考文献[16]。

2 边界条件和网格无关性检验

在ICEM CFD软件中划分六面体结构网格并加密流固耦合面的网格。流体域第1 层网格高度设为0.000 2，以保证第1 层网格质心到内壁面的无量纲距离y+小于1。数值模拟借助ANSYS Fluent 软件进行，固体域材料为316L 不锈钢，流体域材料选取ANSYS Fluent 软件内嵌套的NIST 真实气体物性库中的实际CO2，开启Z方向的重力加速度。入口边界条件选取可压缩流体常用的质量入口，出口边界条件选取压力出口。湍流模型选取SSTk-ω两方程模型，采用压力速度耦合的SIMPLIC 算法，动量、能量、湍动能和耗散率方程的离散均采用二阶迎风格式，松弛因子保持默认。连续性方程、动量方程、能量方程、湍动能和耗散率方程的收敛残差设置为1.0×10-6，当残差曲线和进出口平均温度不发生明显变化，且进出口质量流量的误差小于1%时，可认为计算收敛。

选取质量流速G=100 kg/(m2·s)、热流密度qwi=50 kW/m2、入口压力p=7.6 MPa、入口温度Tin=295.15 K 的工况，分别采用4 套网格计算所给几何模型的加热段压降，结果见表1。

表1 网格数与加热段压降Tab.1 The grids number and the pressure drop of heating section

由表1可知，加热段压降随网格数量增大而增大，但当网格数增大到848 194 时，压降的变化很小。考虑到计算成本和求解精度，本文选取网格数848 194 进行数值模拟。

3 模型验证

为验证数值方法的可靠性，需进行模型验证。因为现存文献中鲜有非均匀热流条件下超临界CO2的传热试验，所以选取Kim 等人[12]在均匀热流条件下进行的试验中的1 组数据(G=238 kg/(m2·s)，qwi=52 kW/m2，p=7.621 MPa)进行模型验证，得到内壁温和主流温度沿管长的分布如图3所示。结果表明，试验和模拟所得内壁温和主流温度的误差均在3%以内，认为模型可靠。

图3 内壁温和主流温度沿管长的分布Fig.3 The distribution of inner wall temperature and bulk temperature along the tube length

4 结果与讨论

4.1 均匀与非均匀加热

与均匀加热的情况不同，非均匀加热时，管壁加热面和未加热面的局部热流密度不同，壁温和对流传热系数也不同。图4给出了p=7.6 MPa、G=100 kg/(m2·s)、qwi=50 kW/m2工况下非均匀加热和均匀加热的对流传热系数和内壁温随流体主流温度的变化曲线。

为便于比较，对于均匀加热，热流密度即实际施加于加热面的热流密度；对于非均匀加热，热流密度为整个圆周面的平均热流密度，实际每处的热流密度并不相同。其中，加热面对流传热系数定义如下：

未加热面对流传热系数定义如下：

总对流传热系数定义如下：

式中，hh、huh、h为加热面对流传热系数、未加热面对流传热系数、总对流传热系数，kW/(m2·K)；Twi、Tb为内壁面温度和主流温度，K。

由图4a)可知，2 种加热情况对应的对流传热系数变化趋势相同：拟临界点前，都随比定压热容的增大而增大；拟临界点后，随比定压热容的减小而减小。非均匀加热时，未加热面的对流传热系数最大，加热面的对流传热系数最小；均匀加热时，加热面的对流传热系数则介于两者之间。

图4b)的内壁温变化反映着图4a)中对流传热系数的变化：非均匀加热时，加热面的局部热流密度最大，对流传热系数最小，内壁温最高。因此，较之均匀加热，非均匀加热受热面壁温更容易超温。

图4 对流传热系数和内壁温变化曲线Fig.4 Change curves of the convective heat transfer coefficient and inner wall temperature

4.2 非均匀加热热流密度的影响

图5为G=100 kg/(m2·s)、p=7.6 MPa 工况下，热流密度对非均匀传热的影响。

由图5a)可以看出，3 种热流密度条件下，总对流传热系数随主流温度的升高先增大后减小，在拟临界点附近出现峰值。拟临界点前，热流密度越小，对流传热系数越大；拟临界点后，3 种热流密度下的对流传热系数相差甚微。这是因为：拟临界点前流体因温度升高比定压热容不断增大，传热能力增强；而以上3 种工况对应的Gr/Re2.7均大于文献[17]所提到的临界值10‒5，传热受物性和浮力2 种因素的影响；由图5b)可知，热流密度越大，浮力越大，浮力效应对传热的抑制越严重。因此：28 kW/m2热流密度下的对流传热系数最大，39 kW/m2次之，50 kW/m2最小；拟临界点附近，Gr/Re2.7骤减，浮力对传热的抑制骤减，对流传热系数达到峰值；拟临界点后，流体比定压热容减小，传热系数减小；远离拟临界点的位置，流体的比定压热容随主流温度变化不大，3 种工况对应的Gr/Re2.7也随主流温度变化不大，且相差甚微。

图5 热流密度对非均匀传热的影响Fig.5 Effects of heat flux on non-uniform heat transfer

图5c)给出了3 种热流密度条件下管壁加热面和未加热面对流传热系数的变化。总体趋势和图5a)相同，但拟临界点前，未加热面的对流传热系数大于加热面的对流传热系数，拟临界点后，两者相差较小。这是因为，拟临界点前，加热面内壁温和主流温度相差较大而未加热面内壁温和主流温度相差较小，随着管壁周向导热和流体对流传热的不断进行，未加热和加热面的内壁面温差逐渐减小，因此，拟临界点后对流传热系数的差值也逐渐减小。

图5d)给出了3 种热流密度下内壁温随主流温度的变化。热流密度越大，内壁温越高。拟临界点前，内壁温变化平缓；拟临界点后，内壁温飞升；远离拟临界点的位置，温度变化又变得平缓。这与图5a)中对流传热系数的变化一致。

4.3 非均匀加热入口压力的影响

选取p=7.6、8.0、8.6 MPa 共3 组工况（G=100 kg/(m2·s)，qwi=50 kW/m2），分析入口压力对非均匀传热的影响，如图6所示。

由图6a)可以看出，不同于低热流密度下压力对对流传热系数的影响，在本文所给的高热流密度下，入口压力越大，对流传热系数越大。结合图6b)可知，在当前热流密度下，Gr/Re2.7远大于文献[17]提到的Gr/Re2.7临界值10‒5，浮力对传热的抑制作用十分明显。拟临界点前，流体因不断吸热，温度升高，比定压热容逐渐增大，传热能力增强；拟临界点处，流体比定压热容达到峰值，对流传热系数也达到峰值；拟临界点后，流体比定压热容骤减，对流传热系数减小；远离拟临界点的位置，流体比定压热容变化缓慢，对流传热系数随之变化缓慢。7.6 MPa 入口压力下流体的物性随主流温度变化的剧烈程度大于8.0 MPa 和8.6 MPa 压力下物性变化程度，管内近壁区流体与主流区流体密度差更大，导致浮力更大，对湍流传热的抑制作用更强。所以，入口压力较大的一组，对流传热系数更大。热流密度较低时，物性变化引起的浮力效应较弱，甚至可忽略其对湍流传热的抑制作用，此时，比定压热容的变化对湍流传热起主导作用。因此，在低热流密度下，压力越小，对流传热系数越大。

图6 入口压力对非均匀传热的影响Fig.6 Effects of inlet pressure on non-uniform heat transfer

图6c)对比了3 种入口压力下，管壁加热面和未加热面对流传热系数的变化。总趋势同图6a)。

图6d)给出了3 种入口压力下内壁温随主流温度的变化曲线。由图6d)可知，8.6 MPa 入口压力对应的内壁温较低，但与其他2 个条件下的内壁温相差不大。

4.4 热流密度和入口压力对加速效应的影响

除CO2物性变化和浮力影响外，加速效应也是一个影响超临界CO2在管内传热流动的重要因素。由前文可知，因流动过程的压降远小于入口压力，可忽略压降引起的加速效应。但是，随着管内流体温度的升高、密度的降低和热膨胀系数的变化，轴向的密度差也会使流体产生加速效应，进而使流体速度梯度发生变化，切应力随之变化，影响湍流的生成、扩散和传热能力，严重时易引发湍流的层流化。McEligot 等人[18]研究了湍流向层流转变过程中的传热特性，提出了无量纲准则数Kv的表达式（式(9)），以判别因密度差引起的加速效应对管内超临界流体传热的影响。

当Kv＜3×10‒6时，流动状态保持为湍流；当Kv≥3×10‒6时，湍流减弱，传热恶化。

图7给出了热流密度和入口压力对加速效应的影响。

图7 热流密度和入口压力对加速效应的影响Fig.7 Effects of heat flux and inlet pressure on acceleration effect

由图7a)可知，Kv随主流温度的升高先减小后增大，拟临界点为最小值。正如式(9)所示：随着温度的升高，流体的比定压热容增大，在拟临界点处达到最大值，此时Kv最小；拟临界点后，随着流体比定压热容减小，Kv逐渐增大。而热流密度越大，流体轴向密度场分布越不均匀，加速效应就越强。

由图7b)可知，Kv的变化与入口压力对浮力效应的影响类似：压力越小，Kv越大，加速效应就越强。原因是7.6 MPa 下CO2的物性变化较8.0 MPa和8.6 MPa 更为剧烈。比较热流密度和入口压力的影响可知，热流密度对加速效应的影响更大，这一点从Kv的表达式也可以看出。本文所有工况得到的Kv都远小于临界值，因此，加速效应的影响可忽略。

4.5 模拟结果与现有关联式的比较

对于亚临界单相流体，经典的Dittus-Boelter 关联式[19]对传热有着较好的预测，而超临界流体因其物性在拟临界点剧烈变化，以及引起的浮力和加速效应，使得Dittus-Boelter 关联式所预测的传热特性和试验结果存在很大偏差。针对这一问题，许多学者在各自试验的基础上，提出了许多超临界CO2在不同条件下的传热关联式。这些关联式大多把主流区流体和近壁区流体物性的比值作为修正项引入Dittus-Boelter 关联式。表2给出了4 种常见的传热关联式。

表2 4 种常见的传热关联式Tab.2 Four kinds of heat transfer correlations

图8为G=100 kg/(m2·s)、qwi=50 kW/m2、p=7.6 MPa 条件下由表2中传热关联式计算所得Nu和本文模拟值对比。由图8可见，不同传热关联式所得Nu相差很大，拟临界点附近更为明显。并且，4 种关联式预测的Nu与本文模拟值相差很大，主要原因在于本文所给的计算模型为非均匀热流密度下的竖直圆管内超临界CO2传热，而4 种关联式是针对均匀热流密度得到的，因此仅适用于各自的试验工况。

图8 4 种常见传热关联式计算Nu 与模拟值对比Fig.8 The Nu calculated by the above four common correlations and the simulated values

为较为准确地预测依据本文计算模型所得Nu，在Dittus-Boelter 关联式的基础上，结合浮力效应和物性变化，提出以下关联式：

式(10)中，浮力效应由近壁区流体和主流区流体的密度差产生，因此，用近壁区流体和主流区流体的密度之比作为浮力修正项。物性变化中，比定压热容的变化最明显，采用平均比定压热容与主流区流体比定压热容之比作为物性修正项。由4.4 节可知，加速效应的影响可忽略不计。利用多元线性回归的方法求得以下传热关联式：

任意选取3 组工况，比较新关联式计算的Nu和本文模拟值，结果如图9所示。由图9可见，新传热关联式与模拟值吻合较好。但应注意，虽然超临界CO2的传热规律具有一定的相似性，但CO2物性的剧烈非线性变化及其带来的浮力效应和加速效应给研究带来了很大的困难。目前，还不存在一个普遍适用且准确的传热关联式。因此，新关联式也只适用于本文中的数据。