超临界二氧化碳矩形回路自然循环稳态数值研究

Steady State Numerical Study on the Natural Cyclic of Supercritical Carbon Dioxide in a Rectangle Loop

摘要: 超临界二氧化碳在超临界区具有良好的热力学性能，作为布雷顿循环的工质，相比于传统的朗肯循环有着更高的效率，自然循环是提升系统非能动安全性的技术基础，超临界二氧化碳循环系统的稳态特性是设计的关键。为研究超临界二氧化碳自然循环稳态特性，采用计算流体力学(CFD)方法，对超临界二氧化碳在小型矩形回路内自然循环换热特性进行三维稳态数值模拟，冷热段高差1.0 m，压力范围7~11 MPa，加热壁面温度320~370 K。模拟结果显示，当流体温度达到拟临界温度附近时，质量流量、换热功率和表面换热系数显著升高。由于在拟临界点附近剧烈的物性变化，受浮升力和重力影响，水平管内局部温度呈现底部向顶部逐渐升高。基于模拟结果拟合了稳态雷诺数Re和修正格拉晓夫数Grm的关系式。为超临界二氧化碳非能动安全系统的设计和分析提供数据基础。

Abstract: Due to the admirable thermal characteristics in the supercritical region, carbon dioxide is considered for use as working fluid of Brayton cycle, which could achieve higher efficiency compared with steam Rankine cycle. Natural circulation is the technical basis for improving passive safety of systems, and the steady-state characteristics are essential for the design of supercritical carbon dioxide power conversion system. In order to study the steady-state characteristics of the natural circulation of supercritical carbon dioxide, a three-dimensional steady-state simulation of the heat transfer characteristics of the natural circulation of supercritical carbon dioxide in a small rectangular loop was carried out by using the CFD method. The simulation parameters are as follows: the height difference of cold and hot section is 1.0 m, the pressure range is 7 - 11 MPa, and the heating wall temperature range is 320 - 370 K. The simulation results show that the local temperature in the horizontal tube increases gradually from the bottom to the top due to the drastic thermal properties changes near the quasi-critical point and the influence of buoyancy and gravity. When the fluid temperature reaches near the quasi-critical temperature, the mass flow rate, heat transfer power and surface heat transfer coefficient increase significantly. Based on the simulation results, the relationship between steady-state Reynolds number Re and modified Grashof number Grm is fitted. The simulation results can provide reference for the design and analysis of passive safety system of supercritical carbon dioxide.

ABSTRACT

Due to the admirable thermal characteristics in the supercritical region, carbon dioxide is considered for use as working fluid of Brayton cycle, which could achieve higher efficiency compared with steam Rankine cycle. Natural circulation is the technical basis for improving passive safety of systems, and the steady-state characteristics are essential for the design of supercritical carbon dioxide power conversion system. In order to study the steady-state characteristics of the natural circulation of supercritical carbon dioxide, a three-dimensional steady-state simulation of the heat transfer characteristics of the natural circulation of supercritical carbon dioxide in a small rectangular loop was carried out by using the CFD method. The simulation parameters are as follows: the height difference of cold and hot section is 1.0 m, the pressure range is 7 - 11 MPa, and the heating wall temperature range is 320 - 370 K. The simulation results show that the local temperature in the horizontal tube increases gradually from the bottom to the top due to the drastic thermal properties changes near the quasi-critical point and the influence of buoyancy and gravity. When the fluid temperature reaches near the quasi-critical temperature, the mass flow rate, heat transfer power and surface heat transfer coefficient increase significantly. Based on the simulation results, the relationship between steady-state Reynolds number Re and modified Grashof number Gr_m is fitted. The simulation results can provide reference for the design and analysis of passive safety system of supercritical carbon dioxide.

Keywords:Supercritical Carbon Dioxide, Steady-State, Natural Circulation, CFD, Heat Transfer

Copyright © 2020 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]、地热系统 [3]、废热回收系统 [4] 以及先进的核能反应堆系统 [5]。安全性是上述系统可持续发展的前提和基础，作为提升非能动安全性的重要途径，对系统的设计和运行至关重要，自然循环成为了持续关注的重点。自然循环是一种非能动循环方式，无需机械泵等驱动装置，仅通过回路中热端和冷端的流体密度差在重力作用下产生驱动压头推动循环进行，将热量从热源带出。自然循环系统结构简单、可靠性好，采用自然循环可以提高系统的固有安全性。因此，研究超临界二氧化碳工质在自然循环中的换热能力，以及影响其流动换热的因素对系统非能动安全性的设计与分析尤为重要。

在超临界二氧化碳自然循环特性数值研究方面，研究者主要采用二维和三维CFD模型分析换热性能。Zhang等人 [6] 采用二维CFD模型研究了冷热源为定壁面温度边界条件下简单超临界二氧化碳自然循环回路内的温度场、速度场和压力分布。Chen等人 [7] 通过二维CFD模型分析了加热段放置位置对超临界二氧化碳流动传热特性的影响。研究了在固定工作压力和冷却段壁面温度下，分析了加热段位于底部水平中段、底部水平左侧、回路左侧竖直中段和左侧底部四种位置下超临界CO₂流动的传热效应和不稳定性。Z. Du等人 [8] 研究了超临界二氧化碳在拟临界点附近的换热行为，分析得到浮升力的影响对超临界二氧化碳自然循环不可忽略，尤其是工况运行在临界点附近。Yadav等人 [9] 以地热能利用为背景，针对加热段和冷却段均为套管式换热器的自然循环回路，采用三维CFD模型，研究了在同一回路内单相、两相和超临界态下循环流量的影响，提出了超临界二氧化碳自然循环稳态特性的表征式，并与Vijayan [10] 实验得到的关系式进行比较。Ajay等人 [11] 采用三维模型研究了摩擦因子对超临界二氧化碳自然循环换热的影响，与Vijayan的Re和Gr_m关系式不同，提出了摩擦因子和Re的关系式。

目前，关于超临界二氧化碳自然循环的研究较少。为研究超临界二氧化碳在自然循环条件下的稳态运行特性，包括温度温度梯度分布、表面换热强弱以及与换热功率密切相关的质量流量。本文采用三维模型，模拟超临界二氧化碳在小型矩形回路内的自然循环流动换热特性，工况覆盖拟临界区和超临界区，使得二氧化碳可以在超临界状态下不发生相变，类似于单相。经模拟分析，发现了二氧化碳在临界点附近物性的剧烈变化，会导致质量流量和表面换热系数的变化，浮升力对水平管道内的温度和速度分布有显著影响。对于往后的系统设计和运行工况的选择需要考虑以上问题。

2. 计算模型建立

本文研究主要通过ICEM17.0软件对既定的几何结构进行网格划分，再通过Fluent17.0软件进行数值计算，物性软件查询软件采用REFPROP。

2.1. 几何模型建立

模拟的矩形回路由加热段、冷却段和两根竖直管段组成，如图1所示。管道内CO₂由底部的加热段(HHX)进行加热，由顶部的冷却段(CHX)进行冷却，除加热段和冷却段外，剩余管段视为绝热壁面。由于加热段和冷却段温度不同，矩形回路内的工质密度随回路位置的变化存有差异，在密度差引起的浮力以及重力的作用下，在矩形回路内形成循环流动，从而进行换热。

Figure 1. Schematic diagram of the geometry of the simulation loop

图1. 模拟回路的几何结构示意图

回路的具体参数如下：

回路直径：d = 0.008 m；CHX长度：L₁ = 0.6 m；HHX长度：L₂ = 0.6 m；水平管道总长度：L₀ = 1.0 m；管道高度：H₀ = 1.0 m；上端绝缘管长度：L₃ = 0.2 m；下端绝缘管长度：L₄ = 0.2 m。

本模拟建立在以下假设条件下：

a) CO₂在矩形回路中以单相状态流动；b) 系统运行在稳态条件下；c) 除CHX和HHX以外，剩余管段为绝热壁面；d) 加热段和冷却段初始壁面为恒温。

2.2. 数值计算方法

本数值模拟采用FLUENT17.0进行模拟，采用ICEM进行三维网格的划分。质量、能量和动量控制方程采用控制体积方法进行离散。动量和能量项以及湍流项均采用二阶迎风格式。壁面设置为无滑移壁面。湍流模型选择RNG k-ε模型，此模型已成功模拟和预测过有关CO₂的研究 [12]。计算过程中同时对温度、压力和速度进行监测，所有监测参数的残差值保持恒定时，计算收敛。不同的参数的残差收敛条件不同，当连续项残差低于10⁻³，能量项的残差低于10⁻⁶时，认为模拟收敛，完成一次计算。

整个回路的网格划分采用O型网格的形式，使单个网格更均匀。截面上最大网格尺寸0.1 mm，最小网格尺寸0.04 mm，网格质量均高于0.4。网格在径向壁面附近进行加密处理。截面网格和弯管处网格如图2(a)和图2(b)所示。经网格敏感性分析，最终网格数量为291,336。

(a) 截面网格 (b) 弯管处网格

Figure 2. Mesh generation of the loop

图2. 回路网格划分

后续结果分析和处理中需要利用以下公式：

截面上的质量流量：

$m={\displaystyle {\int }_0^A\rho udA}$

稳态雷诺数：

$\mathrm{Re}=\frac{4m}{\pi du}$

修正的修正格拉晓夫数：

$G{r}_m=\frac{g\beta d^3{\rho }^{2}Q{H}_0}{A{u}^3{C}_p}$

其中，

$Q=m{C}_p\Delta T$

在计算Re和Gr_m时，对于稳态导热问题，相关的CO₂物性的特征温度取进出口主流的平均温度，由于竖直管道为绝热壁面，加热段出口主流温度等于冷却段入口主流温度，冷却段出口主流温度也等于加热段入口主流温度，即：

$T_{avg}=\frac{T_{CHX}+T_{HHX}}{2}$

修正的格拉晓夫数Gr_m 与经典的定义的格拉晓夫数有所不同，H.H.Bau [13] 等人称之为无量纲传热速率，式中QH₀/Aμ³C_p和温温具有相同的量纲。

2.3. 2拟临界物性分析

CO₂拟临界点附近的物性变化剧烈，图3显示的是在给定压力条件下比热容随温度的变化。从图中可以看出在达到拟临界温度时，比热容瞬间达到最大值，同时随着压力的升高，拟临界温度也逐渐升高。在拟临界点附近物性的剧烈变化对CO₂作为换热工质的影响程度如何值得关注，更重要的是在模拟过程中对CO₂物性的处理上，利用理想气体方程不能准确反映物性剧烈变化的，应采用真实物性。

Figure 3. The changes of C_p in different pressure

图3. 比热容在不同压力下随温度的变化

2.4. 物性处理

由于CO₂在拟临界点附近的物性变化剧烈，因此在模拟过程中需要采用真实的气体物性。计算过程中涉及到的物性有在给定压力条件下的密度、比热容、导热系数和粘度系数在所需温度区间内的对应数值。以比热容为例，如图4所示，本次模拟通过标准的NIST物性查询软件REFPROP，在给定压力下，以1 K为间隔进行所需温度区间内CO₂物性的查询，在物性变化剧烈的区域适当的缩小取值温度的间隔，通过分段线性的方法将物性曲线分成若干段输入到FLUENT中的物性模块中，在计算过程中采用插值法得到确定温度点的CO₂物性。

3. 结果分析

在本文的模拟过程中，加热段和冷却段均采用恒温壁面的方式进行设置，更接近于加热段和冷却段均为套管式换热器的情况。冷却段的壁面温度在整个模拟过程中一直保持305 K，加热段的壁面温度以10 K为步长从320 K变化370 K。回路压力以1 MPa为步长从7 MPa变化到11 MPa。共30组工况，工况如表1所示。

Figure 4. Real properties of CO₂ (C_p as example, the data comes from REFPROP)

图4. CO₂真实物性输入(以C_p为例，物性数据来源于REFPROP软件)

Table 1. Text matrix used for the current investigation

表1. 本文模拟工况范围

3.1. 回路温度分布

图5显示的是在7 MPa回路压力下，管内CO₂主流温度沿回路位置的温度分布随加热壁面温度的变化关系。CO₂在加热段中迅速升温，当流过竖直管段时流体的温度损失很小，随后在冷却段迅速被冷却，再通过绝热的竖直管段进去加热段，如此往复形成自然循环流动。加热壁面温度为320 K、350 K和360 K时的流体流动方向与加热壁面温度为330 K、340 K和370 K时的流动方向不同。回路内流体流动方向的随机流动行为是因为回路呈现热对称和几何对称分布。Chen等人 [14] 也观察到类似的流动方向随机性现象。

Figure 5. Distribution of bulk fluid temperature along the loop

图5. 管内主流温度沿回路位置的分布情况

3.2. 回路横截面温度分布

图6显示的是在压力为7 MPa、加热壁面温度为340 K的条件下，加热段和冷却段不同位置的截面温度和速度沿径向的变化。分别记录了加热段和冷却段的进口、出口和中心的温度沿直径的变化情况。在进口处，管道中心主流温度分布近乎均匀，由于加热段和冷却段均为水平放置。在浮升力的作用下，温度和速度在水平管道横截面的竖直方向上并不是对称分布的，温度在截面上的分布呈现为从底部逐渐向顶部升高(−4.0 mm为底部，4.0 mm为顶部，以0.8 mm作为间隔)，在截面的水平方向上，温度和速度分布是对称分布的。冷热段进口来流为竖直管段，浮升力不会对温度和速度在径向产生影响。

(a) 加热段温度 (b) 冷却段温度 (c) 加热段流体速度 (d) 冷却段流体速度

Figure 6. Distribution of temperature and velocity in cross section

图6. 回路横截面温度和速度分布

3.3. 速度和温度分布

图7显示的是回路加热段和冷却段的不同位置处的温度分布云图，从冷热段中心和出口的温度云图可以直观的看出，温度分布呈现一个底部往顶部升高的不对称现象，由于采用恒定壁面加热和冷却的模拟设置，加热段和冷却段近壁面温度恒定，浮升力影响不重要。如上述相同，水平放置的加热段和冷却段的顶部温度高于底部温度是因为浮升力的影响，Liu G等人 [15] 通过实验研究了浮升力对自然循环挂内的影响。从冷热段进口的温度云图可以看出，在进口位置温度几乎是对称分布的，说明在非加热管段内，浮升力对温度的影响是很小的。

Figure 7. Temperature contour in cross section

图7. 回路横截面温度分布云图

图8显示的是回路整体的速度分布及弯管处的速度分布情况，从图中可以看出回路内CO₂呈现逆时针的流动方向。在弯管处，流体的内侧速度明显高于外部速度，但是经过弯管后，外侧的速度高于内侧的速度，这是由于弯头弯曲下离心力的作用使得原本在弯头内侧的流体运动到了外侧。加之浮升力的影响，在回路上部的冷却段也因此形成了管内顶部速度高于底部速度的现象。

从温度和速度在回路局部的截面上的分布结果可以看出，在浮升力和重力影响下，管内温度和速度分布沿径向不是对称分布的，温度和速度的不均匀分布对自然循环的稳定运行具有一定的影响。因此简单的二维模拟是不能准确得到这些因素的影响结果的。

Figure 8. Velocity distribution of the whole loop

图8. 回路速度分布图

3.4. 质量流量及换热功率

质量流量和循环换热功率是自然循环研究最关注的物理量，图9和图10为不同加热段壁面温度下，质量流量和换热功率随着压力的变化关系。结果表明，在相同的加热条件下，质量流量随着压力的增加而增加，在到达拟临界点后，质量流量有所下降。拟临界点附近的质量流量变化显著，而换热功率在拟临界点附近有所上升，但相比于质量流量的变化，则不那么剧烈。质量流量和换热功率受密度和比热容的影响较大，而密度的拟临界温度在相同压力下总是低于比热容，传热功率在较小温度下的传热功率会有所降低。而在模拟的过程中只能设置冷却壁面和加热壁面在拟临界点的前后，才能确保主流温度经过拟临界点。对于壁面温度为320 K和330 K在压力分别为9 MPa和10 MPa时出现拐点的情况，是壁面温度在这两个压力的临界温度附近且高于临界温度，致使主流温度接近临界温度，密度等物性变化剧烈，从而发生较大的质量流量增加。

Figure 9. Effect of heating wall temperature on mass flow rate

图9. 加热壁面温度对质量流量的影响(HHX)

在较高的加热壁面温度下，加热段的CO₂密度较小，从而质量流量较小，但是换热功率却有所增加，这是因为粘度系数随着温度的升高而减小，粘度系数的减小有助于提高增强换热的雷诺数，从而达到换热增强的效果。

Figure 10. Effect of heating wall temperature on heat transfer power (HHX)

图10. 加热壁面温度对换热功率的影响(HHX)

3.5. 表面换热系数

图11显示的是给定压力条件下，冷却段表面换热系数随加热段壁面温度的变化关系。在较低的压力下(7 MPa、8 MPa和9 MPa)，表面换热系数随着加热壁面温度的升高而逐渐减小，但在较高压力下(10 MPa和11 MPa)，表面换热系数随着温度的升高而升高，到达峰值后逐渐减小。随着压力的升高，峰值在更高的温度出现，这是由于拟临界温度随着压力升高而增大。因此，拟临界点附近具有很高的表面传热系数，即具有很强的换热能力。如对照CO₂物性表，9 MPa，320 K加热壁面温度条件下，也是表面换热系数的峰值所在，因此表面换热系数较大。

Figure 11. Effect of system pressure on surface heat transfer coefficient

图11. 系统压力对表面换热系数的影响

3.6. Re和Gr_m关系式的拟合

为了检验本模拟的正确性，实验结果与前人Vijayan [10] 进行比较，Vijayan实验研究了水的自然循环，得到了体现循环能力的无量纲参数Re和修正格拉晓夫数Gr_m之间的关系式。修正的格拉晓夫数是利用流体的平均温度作为特征温度而计算得到的，因此对整体的循环能力具有概括性的体现，但是不能表征局部换热能力。

图12显示的是模拟数据及拟合曲线同Vijayan关系式的对比，虽然在趋势上都是吻合的，但是数值上模拟得到的数据均未落在该关系式上，且数值均低于经实验得到的关系式。造成这样偏差的原因可能有两个，一是加热方式，Vijayan在实验过程中采用的是恒定热流密度的方式进行加热，模拟采用的是恒温壁面的加热方式；二是工质不同，Vijayan实验中采用的工质是水，模拟中采用的是二氧化碳，且处于拟临界和超临界状态下。

为了更好的对比，与Ajay [7] 的研究方式一样，通过拟合得到的关系式和Vijayan得到的关系式在数值上进行对比，结果如下：

Vijayan关系式 [10]：

$\mathrm{Re}=1.96{\left(G{r}_{m}d/L_t\right)}^{\frac{1}{2.75}}$

模拟结果拟合关系式：

$\mathrm{Re}=1.786{\left(G{r}_{m}d/L_t\right)}^{\frac{1}{2.74}}\left(30000\le \mathrm{Re}\le 90000\right)$

Figure 12. Fitting formula compare with Vijiayan

图12. 拟合公式与Vijiayan公式对比

从关系式的形式上一致，只在系数有较小的差别。上述关于CO₂自然循环模拟得到的关系式正确性得到证实。该关系式在给定的雷诺数范围内，对超临界二氧化自然循环的设计和分析具有一定的参考价值。

4. 结论

本研究通过CFD方法对超临界二氧化碳在简单矩形回路内自然循环换热特性进行了模拟和分析。采用三维建模的方式进行了稳态数值计算，得出了以下结论：

1) 在加热段和冷却段为理想的几何对称和热对称时，自然循环的流动方向是具有随机性的，因此，需要自然循环沿着既定的方向流动，需要按照需求布置不对称性，更利于系统的稳定运行；

2) 在拟临界附近，由于物性的剧烈变化及浮升力的作用，温度和速度在三维空间中的分布不是对称均匀的，呈现管内上部温度高于下部温度，水平管段尤为明显，竖直管段温度分布均匀；

3) 在7~11 MPa范围内，质量流量和换热功率随着加热壁面温度的升高而升高，且在拟临界点处达到最大，随后开始降低；

4) 基于模拟结果，拟合了新的关于超临界二氧化碳流动换热的Re和修正Gr_m关系式，对超临界二氧化碳自然循环设计和分析具有一定的参考价值。

本文初步进行了超临界二氧化碳自然循环稳态特性数值模拟，数值计算结果可以用于以超临界二氧化碳为工质的换热系统的设计和安全运行的分析。后续还要进行有关超临界二氧化碳自然循环瞬态不稳定性研究，以及不稳定性理论研究。