基于OpenFOAM的超临界CO2流动与传热求解器的开发与验证

2021-02-02曾剑锐张尧立

厦门大学学报（自然科学版） 2021年1期

曾剑锐,张尧立,2*,徐宏,叶楷，许松，洪钢,2

(1.厦门大学能源学院,福建厦门361102；2.福建省核能工程技术研究中心，福建厦门361102)

超临界CO2(sCO2)具备良好的传热特性、较易达到的临界温度与压力、无毒、容易获得等特点，可作为良好的循环工质应用于sCO2布雷顿循环系统，降低压缩机功率，提高系统循环效率，并大大减小设备的体积[1]，在未来有望应用于第4代先进核能系统.

由于sCO2在临界点附近物性变化非常剧烈，这对其流动和传热影响很大，sCO2的对流传热问题在理论和应用上存在许多问题有待解决.Jiang等[2-3]对sCO2流过竖直加热圆管的对流传热问题进行了系列实验研究，并用ANSYS Fluent商业软件进行数值模拟，分析热流密度和浮升力等因素对流动传热的影响.Lei等[4]利用Fluent软件对低质量流量下sCO2的竖直流动进行数值模拟，研究流动参数和热物性参数的规律，探索传热机理.

然而，目前sCO2的数值模拟大多是通过商业软件来进行的[5-6]，而商业软件不能够准确地模拟sCO2的对流与传热.在我国大力发展自主知识产权软件的大背景下，OpenFOAM的开源特性有助于打破软件的版权壁垒，自主地开发研究工具.此外，sCO2有3种对流传热模式：正常传热、传热恶化和传热强化.使用OpenFOAM可以定制化开发求解器，对于研究3种传热模式具有很强的拓展性.研究人员曾利用OpenFOAM对sCO2的对流传热展开模拟研究.Xiong等[7]对湍流模型进行修改并植入OpenFOAM中，将新模型用于sCO2流过竖直加热圆管的模拟计算，重点研究新模型对sCO2流动传热规律的预测表现，但他们对如何将sCO2的热物性参数导入OpenFOAM中进行计算更新的介绍甚少.由于sCO2特殊的热物性对计算结果影响很大，所以，本研究具体阐述了sCO2的热物性参数导入OpenFOAM 5.0中并实现精准更新的方法，同时对新开发的热物理模型进行严格的验证，确保不会因为对sCO2热物性处理不当而导致计算结果失准，此部分工作也补齐了适用于OpenFOAM 5.0并且以基于密度的基本热力学类为框架进行开发的热物理模型的缺口；此外，本研究还开发了可应用于sCO2流动与传热的sCO2Foam求解器，并使用该求解器进行模拟，将模拟结果分别与实验数据和/或经验关联式进行对比，以此验证sCO2Foam在sCO2流动传热以及流动阻力方面的适用性.

1 sCO2Foam求解器的开发

图1为8.8 MPa下CO2的热物性随温度变化的情况，可以看出在拟临界点附近，CO2热物性变化非常剧烈，需精准获取每个状态点的热物性参数，才能保证模拟的准确性.

Cp、μ、k、ρ、h和Tpc分别为定压比热容、动力黏度、热导率、密度、焓和拟临界温度.图1 8.8 MPa下CO2的热物性随温度的变化Fig.1Variations of the thermophysical properties of CO2 with temperature at 8.8 MPa

OpenFOAM原有的热物理派生类，例如：更新动力黏度的transport模型、比热容与焓的thermo模型和密度的equationOfState模型，均无法满足sCO2的热物性变化规律.因此，开发OpenFOAM新的热物理模型来解决这一问题，开发流程分为热物性表的制备、原热物理派生类的改写以及热物性库的编译三部分.

1.1 热物性表的制备与读取方式

sCO2的热物性库来自于REFPROP NIST，将其改写成OpenFOAM可读的热物性表.表中的压力为定值p0，虽然计算过程中压力随着迭代求解会有变动，但热物性参数仍以设定的压力值p0为基础进行查表.

温度表由焓表反向插值得到：T=T(p0,h)，迭代过程中流体温度由经过计算得到的焓值来查表更新，温度值更新后，可以根据热物性表查得其他对应的热物性参数值.因此，热物性参数表格形式如下：

ρ=ρ(p0,T);h=h(p0,T);k=k(p0,T);

μ=μ(p0,T);Cp=Cp(p0,T);T=T(p0,h).

(1)

此外，热物性表的数据点采用非均匀布置，在拟临界点附近局部加密，使其在确保精度的前提下减少了计算成本.

1.2 OpenFOAM新热物理派生类

specie类是热物性参数类的基类，它所衍生出的派生类可以看作是更新各个热物性参数的子模型，每个子模型具有特定的功能.以新的子模型tabularEOS为例，状态方程类tabularEOS的构造函数主要定义了从算例配置文件夹constant中读取密度表densityTable的机制，并通过成员函数定义了密度场的计算方式，即返回密度表中的密度值.同理，新子模型tabularTransport定义了动力黏度与热导率的计算方式；新子模型hTabularThermo定义了比热容与焓值的计算方式.这些子模型通过联合编译，组合成同时拥有上述子模型功能的新库，供后续使用.

OpenFOAM中basicThermo基类下的热物理派生类通过引用specie基类下的热物性参数相关类，并利用makeThermo宏函数的调用组合，生成可选的热物理模型.OpenFOAM热物理类的部分继承派生关系如图2所示，其中实线表示OpenFOAM原热物理类继承派生关系，虚线框中为修改后的新派生类，虚线箭头为新生成的类的继承派生关系.可以看出，heRhoThermo继承自heThermo，heThermo又继承自rhoThermo.另外，pureMixture类的模板参数为constTransport，而输运方程类constTransport正好是specie类中的派生类(在新的热物理子模型中，输运方程类为tabularTransport)，这样就将两部分的类联系起来.因此，还需对OpenFOAM中basicThermo基类衍生出的热物理派生类进行改写，并引用specie基类下的新热物理子模型，才能联合编译生成另一个新库.

图2 OpenFOAM部分原热物理类与修改后的热物理类的继承派生关系Fig.2Inherited derivation relation of the original and modified thermo-physical classes of OpenFOAM

由于本研究所用的求解器的热物理类接口为基于密度的基本热力学类rhoThermo，所以热物理派生类需根据rhoThermo类的框架进行改写.首先，新派生类rhoTabularThermo定义了密度、动力黏度和压缩性等，rhoTabularThermo类中的成员函数根据查表法的机制返回了从热物性表中读取的密度值.接着，基于heRhoThermo类开发新派生类heRhoTabularThermo，其构造函数中调用了calculate函数，实现了各热物性参数的更新.温度的更新机制需进一步加以分析，calculate函数中的forAll循环实现了每个网格单元的热物性参数的遍历更新.原本heRhoThermo类中定义的温度更新是由三参数的THE函数经过一系列的调用，得到温度和焓值相关的牛顿迭代递推公式，由焓值计算温度.但是sCO2拟临界区附近的焓值发生剧变，容易导致在这附近的温度求解发散.而新派生类中定义了温度直接由焓表的反向插值来更新，避免了上述问题，同时保证了热物性更新的准确性.

1.3 热物性库与求解器的编译

完成新热物理类开发后，在rhoTabularThermos.C中调用makeThermo宏函数将上述的各个新派生类组合形成的新热物理模型添加到hashTable中，并编译生成热物性库.求解器做相应修改，植入编译好的两个库，最终形成求解器sCO2Foam.同时，在算例的配置文件中写入模型对应的关键词，如下所示：

热物性库通过lookupThermo函数将关键词组合起来，从hashTable中找到相应的元素，在算例运行时就可成功调用热物性库.

若用户需要计算的工质发生改变，只需将constant文件夹中的热物性表替换成新工质的热物性表，无需对求解器进行修改，即可调用新工质的热物性参数.相比于多项式拟合法，热物性库更加具有便捷性和广泛适用性.

2 sCO2Foam求解器的适用性验证

2.1 新模型处理热物性参数的适用性验证

为了验证新开发的热物理模型能够正确地处理sCO2的热物性参数，现将经sCO2Foam读取、迭代运算和输出得到的热物性参数与美国国家标准与技术研究院(National Institute of Standards and Technology,NIST)数据库相应数据进行对比，如图3所示，经sCO2Foam读取、运算并输出的sCO2热物性参数随温度的变化关系与NIST库的数据完全吻合.因此，植入新热物理模型的sCO2Foam可以准确读取并正确处理sCO2的热物性参数，这为后续的模拟计算奠定基础.

图3 经sCO2Foam读取、迭代运算和输出得到的sCO2热物性参数与NIST库相应数据的对比(p0=8.8 MPa)Fig.3The comparison of the thermophysical properties of sCO2 in NIST database with those obtained after being read, processed and outputted by sCO2Foam (p0=8.8 MPa)

2.2 sCO2的流动传热

sCO2流动传热部分适用性验证以Li等[8]的部分实验为基准，初始条件同实验工况，入口温度T0为298.15 K，热力学压力p0为8.8 MPa，入口雷诺数Re0≈9 000，质量流速为3.68 kg/h.模拟基于笛卡尔坐标下的Navier-Stokes方程，能量方程中忽略对低流速影响较小的黏性耗散项、重力项和单位质量动能项.对流项采用二阶linearUpwind格式，压力速度耦合采用SIMPLE算法，湍流模型采用k-ω进行模拟.具体的物理模型如图4所示.头尾两段均为绝热段，中间段施加恒定热流密度，总管长为500 mm，sCO2向上流过竖直圆管.

图4 物理模型与坐标系统Fig.4Physical model and coordinate system

分别划分475万,576万与700万的结构化网格进行网格无关性分析，经比较，三者的壁面温度相对误差不超过0.05%.考虑到资源消耗成本，本研究选用475万网格进行模拟，并在近壁面处划分相对精密的网格，对于所有计算工况，满足y+在0.5左右.

选取3个sCO2传热经验关联式：Dittus-Boelter[9]、Gnielinski[10]和Jackson[11]，与本研究模拟的沿圆管轴向的壁温结果一同和实验数据进行对比.图5给出了不同热流密度(qw)下，三者的数据比较.可以看出，相较于3个sCO2传热经验关联式给出的沿圆管轴向的壁温分布，本研究模拟结果与实验值更为贴近.壁温沿轴向不断上升，浮升力和流动加速效应对管内对流传热影响不大，此时处在sCO2传热正常的区域.在热流密度分别为6 498和13 626 W/m2条件下，本研究模拟的壁温值与实验值的最大误差分别为0.13%和0.22%.可见对于浮升力和流动加速效应影响不显著的情况下，sCO2Foam采用k-ω湍流模型能较好地模拟出与实验相近的结果，sCO2Foam可用于sCO2流动传热问题的研究.

图5 模拟结果和经验关联式计算的壁温与实验数据的比较Fig.5The comparison between the experimental data and the results obtained both by simulation and the empirical correlations

2.3 sCO2的流动阻力

通常，流体在水平圆管中流动的总压降可以用下式计算[12]：

Δp=Δpf+Δpac，

(2)

(3)

(4)

其中，Δp为总压降，Δpf为摩擦压降，Δpac为由流动加速引起的压降，f为摩擦阻力系数，L为圆管长度，D为圆管直径，u为流速，G为质量流速，ρ为密度，下标out和in分别表示出口处和入口处.因此Δpf=Δp-Δpac，摩擦阻力系数f可从经验关联式中获得.

对水平圆管的加热段进行流动阻力分析，其物理模型和基本工况与上述一致.将模拟计算的摩擦压降与摩擦阻力系数经验关联式进行联合对比，选取的3个摩擦阻力系数经验关联式(Itaya[13]、Blasius[14]和Filonenko[15])的表达式如下：

(5)

(6)

fF=(1.82lgRe-1.64)-2.

(7)

图6为不同入口流速下，sCO2Foam计算求得的sCO2摩擦压降与通过经验关联式推算的摩擦压降的对比.整体上看，模拟计算出的摩擦压降与经验关联式估计的摩擦压降均较为接近，其中模拟结果与式(5)的压降预测最为接近.3种流速从低到高排列，模拟结果与3个经验关联式(式(5)～(7))计算出的摩擦压降的误差范围分别为4.87%～6.73%，4.84%～6.98%和4.99%～7.27%.可见sCO2Foam可用于sCO2的流动阻力问题的研究.

图6 不同入口流速下，模拟计算的摩擦压降与3个经验关联式推算的摩擦压降的对比(qw =13 631 W/m2)Fig.6Comparison of the frictional pressure drops obtained by simulation with those by the empirical correlations at different inlet velocities (qw =13 631 W/m2)