¹四川大学物理学院，四川 成都

²中国核动力设计研究设计院，四川 成都

收稿日期：2021年3月27日；录用日期：2021年4月15日；发布日期：2021年4月29日

摘要

超临界二氧化碳在拟临界点附近具有良好的导热性能和易压缩性，采用超临界二氧化碳作为换热工质的布雷顿循环系统拥有可观的系统热效率，但物性的剧烈变化可能会导致流动不稳定性问题。采用三维数值模拟方法对超临界二氧化碳双通道流动不稳定性问题进行了研究，与实验进行了对比，分析了不稳定起始至结束整个过程中双通道间的状态变化。已有的数值结果表明超临界二氧化碳流动不稳定性分为三个阶段，随着加热功率的线性增加，流动状态经历波动起始、持续波动和波动结束三个过程。质量流量、流体密度和主流温度均会呈现正余弦波动状态，三者存在一定的相位延迟。壁面温度和表面换热系数会随着流体的波动而出现相应的变化，壁面出现最大温度点和最大表面系数换热点。

关键词

超临界二氧化碳，双通道，流动不稳定性

Three-Dimensional CFD Numerical Study of Supercritical Carbon Dioxide Two-Channel Flow Instability

Jiajian Huang¹, Yuan Zhou^1*, Yanping Huang², Guangxu Liu²

¹College of Physical, Sichuan University, Chengdu Sichuan

²Nuclear Power Institute of China, Chengdu Sichuan

Received: Mar. 27^th, 2021; accepted: Apr. 15^th, 2021; published: Apr. 29^th, 2021

ABSTRACT

Supercritical carbon dioxide has outstanding thermal conductivity and compressibility near the pseudo-critical point. The Brayton cycle system which uses supercritical carbon dioxide as heat exchanger has considerable thermal efficiency. However, the drastic change of physical properties may lead to the problem of flow instability. Three-dimensional numerical simulation is used to study the supercritical carbon dioxide two-channel flow. The calculated results are compared with the experimental results and the state changes between the two channels are analyzed during the whole process. The numerical results show that flow instability of the supercritical carbon dioxide is divided into three stages. With the linear increase of heating power, the flow state goes through three processes: the beginning of fluctuation, the continuous and the end of fluctuation. The mass flow rate, fluid density and flux temperature show sine and cosine fluctuation state, and there is a certain phase delay. The wall temperature and surface heat transfer coefficient will change with the fluctuation of the fluid, and the maximum temperature and the maximum surface heat transfer coefficient will occur on the wall.

Keywords:Supercritical Carbon Dioxide, Two-Channel, Flow Instability

Copyright © 2021 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

随着环境污染和化石能源短缺问题的日益加剧，亟需寻找先进的更加效率的能量转换线系统 [1]。采用超临界二氧化碳布雷顿热力循环方案的能量转化系统在近二十年来备受关注，超临界二氧化碳工质无相变、密度大、拟临界附近易压缩以及布雷顿热力循环压比低等使得超临界二氧化碳系统拥有结构紧凑、体积小、热效率和安全性高的特点 [2]。超临界二氧化碳简单布雷顿循环系统中工质再汽轮机做功后会经历回热、冷却和再压缩过程，为降低压缩功耗，将冷端温度控制在工质拟临界温度附近，在拟临界附近工质的物理性质会剧烈变化，这种剧烈变化可能会导致传热恶化和流动不稳定性现象的发生。

对于超临界流体不稳定性研究，起初，国内外许多学者是在第四代反应堆超临界水堆背景下，开展了有关超临界水流动不稳定性的研究。大部分关于超临界流体流动不稳定性的研究都是基于简单自然循环单通道回路。Ting Xiong [3] 等人采用实验研究的方法研究了超临界水在细长双通道实验段内的流动不稳定性现象，给出了不稳定性起始功率的判断方式，得到了不稳定功率阈值与系统压力和工质进口温度之间的关系。Xixi [4] 等人通过实验和数值模拟对超临界水双通道流动不稳定性进行了深入研究，实验采用尺寸为19 × 6.5 mm的厚不锈钢管，得到了关于不同进口温度、不同运行压力、不同质量流量和不同加热功率的不稳定边界图，但厚壁通道不会加强系统的不稳定性；数值模拟了双通道间流量异相振动的不稳定性现象，同时分析了热工参数对振幅和周期的影响。Wenyu Wang [5] 等基于超临界水循环锅炉采用简单回路对超高压超临界水不稳定性进行了实验研究，实验系统的包含回热部分，实验得到了表征不稳定性的无量纲参数图，将不稳定性划分为两个区域：一类不稳定性和二类不稳定性。陈宇宙 [7] 等人采用自然循环实验回路，研究了超临界水双通道流动不稳定性，在有限工况内观察到流动不稳定性现象。随着二十世纪末，高温涡轮系统、高压紧凑型换热器及高温材料等领域的快速突破，核动力技术发达国家开始启动对超临界二氧化能量转换技术的研发计划，越来越多的学者开始对超临界二氧化碳不稳定流动问题进行研究。Sharma [6] 等人建立了包含水平竖直加热和水平竖直冷却等多种加热冷却方式的简单自然循环回路，研究了超临界二氧化碳在单通道自然循环条件的换热及流动不稳定性，实验只在水平加热和水平冷却的布置方式下，且加热功率0.3~1.0 kW，进口温度接近临界温度工况，发现了流动不稳定性。Vijay Chatoorgoon [7] 等人采用自然循环回路，对超临界二氧化碳双通道不稳定性进行研究，实验研究了在一定的运行压力下不稳定发生的功率起始点与流量、进口工质温度的关系，得到了不稳定流动的起始功率和流量的对应关系。对于超临界二氧化碳流动不稳定的数值研究很多，如Lomperski [8]，Guangxu Liu [9]，Lei Zhang [10] 等均是对单通道回路的自然循环不稳定性进行研究。

基于已有研究成果分析，目前对超临界流体流动不稳定性的研究主要以超临界水为主，对于超临界二氧化碳的流动不稳定性研究较少，尤其是平行通道间的流动不稳定性研究。本文运用计算流体力学软件Fluent，对超临界二氧化碳在双通道内流动不稳定性进行三维数值研究，采用真实二氧化碳物理性质，研究不同进口质量流量、进口温度和不同运行压力对不稳定性起始功率点的影响，采用表征超临界流体不稳定性的无量纲参数绘制不稳定性边界。将为超临界二氧化碳双通道流动不稳定性提供理论基础

2. 研究对象

本文主要研究对象超临界二氧化碳在双通道结构内的流动不稳定性，通道间的流动不稳定现象属于局部振荡，因此，在选取物理模型时，只关注双通道结构之间发生异相振荡即可。采用CFD方法(计算流体力学)，构建三维全尺寸模型，对超临界二氧化碳流动不稳定现象进行数值模拟，采用Fluent17.0作为计算分析工具。

2.1. 物理模型

图1显示的是双通道几何结构，该结构尺寸为Vijay Chatoorgoon [7] 自然循环双通道实验回路的实际管道尺寸。图2显示的是双通道结构横截面的网格结构，为消除圆形管道边缘网格畸变，采用O型网格划分方式。

图1. 双通道几何结构

在图1的双通道结构中，包含进出口管段、预混器和双通道管段，所有管道均为圆管。通道1和通道2采用相同几何尺寸，管道内径为16.51 mm，长度为2100 mm，与预混器连接的进出口管段长度为300 mm，有效加热部分为1500 mm，进出口管道内径20 mm，长度100 mm，双通道与进口管道和出口管道几何对称分布。计算考虑重力方向与流动方向相反。

2.2. 网格划分及敏感性分析

在进行数值计算之前，利用ICEM对双通道几何结构进行网格划分。图2显示的是双通道结构横截面的几何网格划分，进出口管段、预混器和双通道结构均采用结构化网格划分，采用T型管道O型网格划分方式。采用结构化O型网格划分，消除了靠近壁面处的网格畸变，有利于计算精度的提高。

图2. 双通道结构网格划分截面图

通过恒定加热功率的计算方法，通过稳态计算，得到不同网格数量下管内主流温度沿加热管道的分布情况，从图3可以看出，当网格数量为1701512和2352563的时候，管内主流温度沿管道分布情况基本一致，同时考虑计算时间，网格数量为1701512已经满足计算要求，以下计算均采用该网格数量划分方式。

图3. 网格敏感性分析

3. 计算模型

3.1. 湍流模型

由于超临界二氧化碳在拟临界处的物性变化剧烈且导热系数不高，在加热过程中，近壁面处的一层流体由于加热而密度急剧减小，该层流体的流速会略大于主流速度，产生近似的分层流动，因此湍流模型采用SST模型，Xixi [4] 和Li Jingjing [11] 等人已经证实在超临界流体计算中采用SST湍流模型更加准确。对于湍流模型，本文考虑采用SST模型，图4显示的是SST模型和标准K-Epsilon模型计算结果对比，SST模型计算结算结果流量的双通道对称性更合理，同时不稳定性起始功率值更接近实验值6.6 kW，因此本文采用SST模型作为湍流模型。对于SST模型，湍动能方程式(1)和湍流耗散方程式(2)为：

$\frac{\partial \left(\rho k\right)}{\partial t}+\frac{\partial \left(\rho k{u}_i\right)}{\partial x_i}=\frac{\partial }{\partial x_j}\left[{\tau }_k\frac{\partial k}{\partial x_j}\right]+G_k+G_b+S_k-Y_M$ (1)

$\frac{\partial \left(\rho \omega \right)}{\partial t}+\frac{\partial \left(\rho \omega u_i\right)}{\partial x_i}=\frac{\partial }{\partial x_j}\left[{\tau }_{\omega }\frac{\partial \omega }{\partial x_j}\right]+G_{\omega }+D_{\omega }+S_{\omega }-Y_{\omega }$ (2)

图4. SST模型和标准K-Epsilon模型的计算结果比较

3.2. 超临界二氧化碳物性

在拟临界点附近，超临界二氧化碳的物理性质会发生剧烈变化。图5显示的是超临界二氧化碳在不同压力下随温度的变化。从图5中可以看出，超临界二氧化碳的比热容、密度、导热系数和动力粘度在拟临界点附近剧烈变化。在FLUENT软件中的二氧化碳物性不具备上述特征，计算需要采用真实的二氧化碳物性。计算采用的二氧化碳的物性通过NIST(National Institute of Standards and Technology)物性查询软件REFPROP获得。

4. 计算结果与讨论

4.1. 与实验对比

为验证计算方法的可行性，从Vijay Chatoorgoon [9] 的超临界二氧化碳自然循环双通道流动不稳定性研究的工况中选择3组工况进行模拟，通过比较计算得到的不稳定出现的功率起始点与实验测得的不稳定出的功率起始点进行比较。计算结果如图6所示，计算结果总体偏小，误差在15%以内，对于超临界流体的计算，认为该方法可以用于超临界二氧化碳流动不稳定性的模拟。

图5. 超临界二氧化碳的物理性质

图6. 数值计算结果与实验对比

4.2. 不稳定流动现象

4.2.1. 不稳定流动过程分析

图7显示的是超临界二氧化碳双通道发生流动不稳定现象的典型工况。压力7.4 MPa，加热进口温度288.16 K，质量流量0.0938 kg/s，壁面热流密度上升速率700 W/m²s。通过随着流动时间线性增加功率的方式，寻找不稳定发生功率起始点。随着功率的增加，双通道内的流体温度同步增加，质量流量基本保持恒定，当加热功率增加到5.851 kW时，通道内的质量流量和加热段出口温度开始出现轻微的抖动，随之双通道内的流体温度开始出现剧烈的振荡，两通道内的质量流量和出口流体温度呈现此消彼长的异相振动形式，温度的不稳定发生功率要比质量流量发生不稳定振荡的加热功率稍有滞后。随后逐渐增加功率值，振动幅值逐渐增大，当功率达到8.704 kW时，振动幅值达到最大，最后增加功率，振动幅值开始减小，即双通道的不稳定异相振荡逐渐消失，再升高功率，双通道之间的流量一直保持稳定流动状态。

4.2.2. 不稳定功率判断

对于不稳定功率起始点，按照双通道间瞬态流量出现明显的异相振动，且在该波动后双通道间流量仍能做持续性振动，该流量点对应的加热功率即为不稳定功率起始点。图8为压力7.4 MPa，进口温度为288.98 K工况下流动不稳定功率起始点的取值方法，可见稳态时稳态流量为0.0469 kg/s，在点A时双通道间的瞬态流量出现明显的异相振动，在功率点A以后，双通道的瞬态流量逐渐增大，通道间流量呈现此消彼长的状态。

图7. 不稳定性流动现象

图8. 不稳定起始点和不稳定功率区间

4.3. 流动不稳定过程分析

4.3.1. 波动情况分析

对于超临界二氧化碳双通道流动不稳定现象，流体的流量、密度和温度在波动期间相互关联，形成相互反馈。通过典型工况，对并联双通道流动不稳定现象进行瞬态模拟计算，观察波动期间流量、密度及温度随流动时间的变化规律，选取工况的压力为7.4 MPa，进口温度为280 K的欠临界状态，通道的总质量流量为0.1 kg/s，加热速率700 W/m²s。沿通道以150 mm为间隔布置11个监测节点，分别对流体密度、质量流量和主流温度在波动期间进行监测。

1) 流体密度沿轴向位置波动

流体密度沿轴线不同位置处，随流动时间的波动情况如图9所示，流体密度在轴向不同位置处均呈现周期性波动，振幅在加热段靠近出口位置处较小，在加热段的中部较大，这是由于处于实验段进口处于欠临界状态，中部的流体温度接近拟临界温度，密度在拟临界点附近波动较大，因此波动幅值也越大。

图9. 密度沿加热段不同位置的波动情况

2) 质量流量随轴向位置波动

图10图显示的是两个通道轴向不同位置处质量流量随流动时间的波动情况，质量流量波动的剧烈程度沿着进口到出口逐渐减小，这是由于密度沿轴向逐渐降低，且在局部位置剧烈减小，流量的可波动值随之下降。通道间各位置处的质量流量波动均呈现严格的异相振荡。

图10. 质量流量沿加热段不同位置的波动情况

3) 主流温度随轴向位置波动

图11显示的不稳定流动期间主流温度沿轴向的波动情况，温度的波动与密度的波动不同，虽然在非常靠近进口位置处的波动较小，但在离进口稍远的位置和靠近出口的位置温度的波动较为剧烈，而在管道中段温度的波动性较小，这是由于工质的温度在拟临界附近时Cp值较大，相应的载热能力较强，温度的波动幅值更小。可见，超临界二氧化碳并联通道流动不稳定性受其拟临界区物性的变化规律影响较大。温度与质量流量的波动相位与密度的波动相位具有一定的延迟性。

图11. 主流温度沿加热段不同位置的加热情况

4.3.2. 单个周期内不稳定过程分析

双通道异相振荡过程中，管内流体状态如密度和管道的壁面状态随着流动时间是周期性变化的，对单个周期的流体密度分布、管壁的温度和表面换热系数进行记录分析，均呈现周期性变化且通道间的状态变化趋势相反。

1) 单个周期内密度沿轴向分布

由于进口处于欠拟临界状态，密度分布会在轴向位置的某处剧烈降低，加热段出口位置处的密度可降低为进口密度的1/4。随着双通道异相振荡的进行，密度剧烈下降的区域会随流动时间沿轴向位置移动如图12所示，且两个通道间的密度分布变化趋势完全相反，密度的剧烈变化会导致管内局部位置的压力分布不同，通道间密度的交替变化使得不稳定振荡过程持续进行。

图12. 密度沿着轴向位置随时间的分布情况

2) 壁面分析

由于管内流体的物性在跨临界温度时发生突变，流动传热也呈现剧烈的变化，过程变得复杂，从图和图可以看出单个周期内通道一和通道二的表面换热系数沿轴线分布是极不均匀的。两个通道在同一时刻的表面换热系数的分布是不同的，但变化的形式完全相同，都是周期性的增加或减小，只是在时间上相差了半个周期的时间。

表面换热系数的剧烈变化导致壁面温度分布极不均匀，双通道的表面换热系数分布分别如图13和图14所示。在局部表面换热系数较大的区域，壁面温度非常低，这里指的表面温度近似管内壁温度，由于局部的表面换热系数极大，且超临界二氧化碳的比热容在该区域极大，形成瞬时传热，此刻工质的携带热量的能力很强，最终形成了局部温度很低的现象，壁面温度的分布分别如图15和图16所示。随后工质的温度升高，逐渐大于拟临界温度，密度、比热容和导热系数等开始急剧减小，换热能力降低但不会消失，管壁温的最高温度基本出现在加热段的出口位置处，最终呈现出进口壁面温度与出口温度相差较大的现象。

图13. 通道1的表面换热系数在单个周期内的变化情况

图14. 通道2的表面换热系数在单个周期内的变化情况

图15. 通道1的壁面温度在单个周期内的变化情况

图16. 通道2的壁面温度在单个周期内的变化情况

5. 结论

本文通过三维数值计算的方法，成功地实现了对简单的双通道结构内流动不稳定现象的模拟，与实验中不稳定起始功率点的对比，误差均在15%以内，分析了不稳定波动期间超临界二氧化碳工质的密度、质量流量和主流温度的变化，以及加热段管壁表面传热系数和壁面温度的分布情况。在既定的双通道结构和模拟工况下，得到以下结论：

1) 超临界二氧化碳在并联双通道结构下发生不稳定流动现象时，随着加热功率的增加，通道间质量流量、密度和温度经历不稳定起始、不稳定振幅增加和不稳定性现象消失的过程；

2) 发生不稳定流动现象时，双通道间流体的主流温度、密度和质量流量均呈现严格的异相异相振荡，流体主流温度的波动受拟临界温度的影响较大，较质量流量的波动在流动时间上存在明显的延迟，密度的波动趋势与质量的波动趋势相反；

3) 超临界二氧化碳的不稳定流动对管壁面状态的影响明显，实验段壁面温度和表面换热系数呈现周期性变化，单个周期内壁面温度沿轴向分布存在巨大差异。

符号说明

文章引用

黄家坚,周 源,黄彦平,刘光旭. 超临界二氧化碳双通道流动不稳定性三维CFD数值研究
Three-Dimensional CFD Numerical Study of Supercritical Carbon Dioxide Two-Channel Flow Instability[J]. 核科学与技术, 2021, 09(02): 113-125. https://doi.org/10.12677/NST.2021.92013

参考文献