浮升力和流动加速对超临界CO2管内流动传热影响*

张海松 朱鑫杰 朱兵国 徐进良 刘欢

(华北电力大学能源动力与机械工程学院, 低品位能源多相流与传热北京市重点实验室, 北京102206)

(2019 年 10 月 8日收到; 2019 年 12 月 5日收到修改稿)

基于单相流体的概念, 超临界流体的异常传热行为已经被研究很多年了, 但是关于其流动传热机理仍没有统一的认识. 本文通过理论分析和实验研究了超临界二氧化碳在竖直管内向上流动过程中, 浮升力和流动加速效应对其流动结构和传热过程的影响. 结果表明, 没有确凿的实验证据表明超临界流体的异常传热行为是浮升力和流动加速直接导致的, 存在的估计浮升力和流动加速效应准则均是在常物性流体的基础上, 做了大量假设得出的, 不同的研究者采用浮升力和流动加速准则分析超临界流体的传热恶化得出的结论不一致.最后, 基于拟沸腾理论分析超临界流体的传热恶化过程, 提出超临界沸腾数区分了超临界流体正常传热与恶化传热的转换边界, 为超临界流体流动传热研究提供新思路, 超临界沸腾数对建立用于不同技术的超临界流体动力循环的最佳运行条件具有重要意义.

1 引言

以超临界流体为工作媒介的能量循环, 在热力循环效率、能量容量、安全性和经济性以及使用性能等方面具有很大的提高潜力, 因此, 被广泛地应用于工业领域中[1]. 最近, Xu等[2]提出了一个1000 MWe S-CO2的发电设计方案, 主蒸汽参数为 620 ℃/30 MPa, 功率效率可达到 48.37%, 这优于以超临界水为工作媒介的水蒸气朗肯循环. 毫无疑问, 这些新系统的设计和优化都需要深入地了解S-CO2独特的热工水力特性, 尤其是一定条件下发生的传热恶化现象[3].

较早发现超临界流体传热危机现象的是Shiralkar和Griffith[4], 他们在竖直光管中发现热流密度高于一定值后, 加热壁面壁温会突然升高至 863 K, 并且热流密度越高, 壁温峰值越大. 此后, 各国学者对超临界流体的传热恶化现象展开了研究, 主要包括 H2O, CO2及有机工质等流体, 并提出了诸多传热恶化的发生机理和预测传热恶化发生模型[5−16]. Bourke等[5]通过对超临界二氧化碳传热实验研究发现, 高热流密度下, 传热恶化仅发生在向上流动过程中, 认为传热恶化是由浮升力引起. Bae[6]通过实验发现向下和向上流动加热过程中均存在传热恶化现象, 向下流动恶化程度要比向上恶化程度低. Brassington和Cairns[7]对超临界氦进行了研究, 发现对于向上流动, 可以用浮升力诱导传热恶化很好地解释, 但是一些壁面峰值没有出现在临界点附近, 认为浮升力并不是唯一的影响因素. Hall和 Jackson[8]根据常物性的假设, 结合自然对流与强制对流解释了超临界流体局部传热恶化现象, 认为浮升力是导致恶化的原因, 并提出了一个可以预测传热恶化发生的理论模型.Mceligot等[9]提出了一个无量纲参数, 描述了由于加热而产生的流动加速度的影响. Liu等[10]通过理论分析和实验研究了浮力和加速度对混合对流的影响, 并发展了估计浮升力和流动加速效应的新准则. Huang等[11]回顾了公开文献中的浮升力和流动加速准则, 发现存在的区分强制对流和混合对流的准则估算与实验数据有较大的差异, 不同研究者采用相同的准则得到的临界值不一致. 除此之外,有研究者发现浮升力准则临界值还与流体流动方向、工质种类有关[12]. 从公开的文献来看, 缺少确凿的实验证据, 证明超临界流体传热恶化发生时,浮升力或流动加速是直接原因, 研究者仅仅通过流动方向和常物性假设去解释超临界流体传热异常行为, 认为浮升力或流动加速是导致超临界流体传热恶化的原因, 得出的结论仅在各自的实验范围内适用性较好, 这些无量纲准则不具有通用性. 除此之外, 研究者进行实验和理论研究时, 可能忽略了一些因素, 导致不同研究者得出了不同的结论[13−17].近年来, 超临界流体的单相概念不断受到物理学家的质疑, 2010年Simeoni等[18]通过非弹性X射线散射和超临界氩分子动力学证明了Widom线(WL)可以将超临界区域划分为两个不同的动力学机制:类汽和类液区域. 2013年, Gorelli等[19]通过分子动力学模拟表明, 超临界状态下的水可以通过WL划分为类汽和类液两个区域. 2015年, Banuti[20]采用拟沸腾概念解释了超临界流体跨越WL, 从类汽到类液转变这一过程, 并给出了临界点以上汽液共存的拟沸腾线, 这个转变过程的WL和临界点的汽液共存线类似.

本文研究了超临界CO2在垂直向上流动过程中的流动传热特性, 检查了浮升力和流动加速效应对局部传热恶化的影响, 工况参数范围为: 管子内径din= 8 mm, 质量流速G= 100—1500 kg/(m2·s),热流密度qw= 20—400 kW/m2, 进口压力Pin=8.220—20.822 MPa. 最后, 根据实验结果, 依据拟沸腾理论对超临界流体的恶化进行解释和研究, 为超临界布雷顿循环设计及优化提供理论依据和参考.

2 实验系统和数据处理

实验系统装置如图1所示, 实验回路主要包括主循环回路系统、冷却水循环系统以及二氧化碳充液回收系统, 实验系统设计压力和温度分别为30 MPa和 500 ℃. 实验所使用的工质为纯度高于99.0%的工业CO2, 工质从高压泵出来分为两路:一路流过缓冲罐, 用以平衡系统的流量和压力脉动, 然后进入预热器并加热到实验所需要的入口温度, 最后工质流入实验段; 另一路在实验过程中作为旁路调节, 可以结合高压泵对实验流量进行精确的调节, 压力通过背压阀进行调节, 其中电加热系统采用低压大电流通过不锈钢管本身的电阻直接发热, 通过连续调节电压, 实现电加热功率连续改变, 整个系统的最大加热功率为 160 kW. 实验中涉及到的主要仪器及相关参数的不确定度如表1所列.

表1 测量仪器的精度和范围Table 1. Accuracies and ranges of measuring instruments.

式中Pe为电加热功率;h为热平衡效率, 热平衡效率的定义为实验段达到稳态热平衡时, 流过整个加热段流体的焓升与电加热功率的比值,

其中m为质量流量;ib,out和ib,in分别为流体的出口和进口温度对应的焓值, 这两个值分别通过NIST物性软件查得;U和I分别为电压和电流.给定入口温度和质量流速, 假设主流焓升沿加热段轴向线性增加, 根据能量守恒, 可计算任意轴向位置的流体主流焓值:

图2 实验段Fig. 2. Vertically positioned test tube.

垂直光滑圆管内壁温可简化为具有内热源的一维稳态导热问题, 圆柱坐标下的导热微分方程描述为

式中j为内热源,lw为壁面导热系数. 给定边界条件:r=R,t=tw,o, dt/dr= 0, 求解得

式中din和dout分别为管子的内、外径;qw为表面热流密度;tw,i和tw,o分别为管子的内、外壁温. 实验段内表面热流密度可由下式计算:

3 浮升力和加速度理论分析

3.1 浮力对流动结构和湍流输运的理论分析

考虑一半径为R的充分发展圆管管道内流动,二维坐标系下的简单流动几何模型如图3所示. 对动量方程从壁面到通道中心积分有

式中右边第一项是流体受热膨胀引起的加速项, 第二项是壁面产生的剪切力项, 第三项是流体剪切力项, 最后一项是由于密度差引起的浮升力项. 因此,对于竖直管道中向上流动的单相流体, 任意位置y处的剪切力可以表达为

其中dBu和dAc分别为与浮升力和流动加速有关的厚度, 从(9)式可以看出, 浮升力和流动加速总是使剪切力趋于减小. 为了得到方程(9)的解, 将任意位置的剪应力积分:边界条件为时则可得圆管内的圆柱形流体元的力平衡:

图3 圆管中简单的二维流动几何模型Fig. 3. Vertically positioned test tube.

代入 (10)式, 可得

其中

Bu为浮升力参数, 表征了浮升力与壁面剪切应力的相对大小, 决定了右侧第二项的重要性. 如果Bu, 浮力的影响可以忽略不计; 当Bu足够大时, 剪切力分布会发生改变, 在极端情况下, 剪切力会变为负值.

3.2 流动加速度对流动结构和湍流输运的理论分析

当流体接收到壁面的热量后, 由于热膨胀会使流体产生加速效应, 考虑物性剧烈变化的单相流体沿x轴竖直向上流动, 应用能量和质量守恒及与温度相关的流体性质, 假设流体的密度变化主要是由轴向温度变化引起的, 密度沿流动方向的变化可以表示为

定义无量纲量r*=r/rpc,X=x/din,rpc为拟临界温度对应的密度, 则方程(15)变为

式中的b为膨胀系数,cP为定压比热. 管内流动过程中的任意横截面均满足质量守恒, 这里r和u存在如下关系:

定义无量纲速度为

其中du/dx的值可以通过流体的膨胀率估计, 应用理想气体假设和能量守恒、质量守恒定律, 可以近似得到

因此, 在边界层区域, 由于低密度层的存在,壁面附近流体的速度比核心区流体的速度低, 因此, 主流流体加速度产生一个压力梯度, 该压力梯度均匀作用于边界层, 由于密度梯度发生变化, 导致速度梯度发生变化, 最终改变了近壁面剪切力分布. 为了估计流动加速对剪切力的改变, 忽略浮升力的影响, 因此, 可以得到类似上述浮升力的无量纲表达式

这里

(20)式右侧第二项的重要性将由Ac决定,Ac为流动加速参数, 表示流动加速度和壁面切应力的相对大小, 当Ac时, 流动加速效应可以被忽略. 实际上, 由于热膨胀使流体加速导致湍流边界层的“分层”现象, 对超临界流体的传热过程有着重要意义, 但是没有引起足够的重视.

公开文献中用来估计浮升力和流动加速对剪切力改变的准则均做了大量假设, 大都是分别考虑了浮升力和加速效应, 实际上, 浮升力和加速效应总是同时存在的. 对 (10) 式变形, 并无量纲化, 可以得到一个同时包含浮生力项和流动加速项的无量纲方程:

从以上推导来看, 得到的准则均是无量纲准则, 这些准则仅在某些情况下近似满足假设条件,并且都假定了壁面剪应力与切应力之间的一般关系. 事实上, 当传热恶化发生时, 并没有确凿的实验证据表明浮升力和流动加速效应是直接原因, 这些假设与实际情况偏离较大. 从公开的文献来看,不同的研究者得出的结论也不一致. 除此之外, 这些模型均在“两区”模型的假设基础上进行理论分析, 这需要一个合理的特征温度, 但是“两区”模型并不是总是成立的, 模型的有效性很大程度上取决于壁面热流密度的水平. 本文选取两个典型准则即浮升力准则[8]和流动加速准则[9], 通过实验数据检验存在的经典准则的适用性, 即:

式中的Gr为格拉晓夫数,Re为雷诺数.

4 实验结果及讨论

4.1 局部壁温Tw,in, Bu, Ac随主流焓值ib的分布

不同质量流速下的Tw,in,Bu以及Ac随主流焓值分布, 如图4—图6所示, 所有传热恶化均发生在主流温度Tb

当质量流速为700 kg/(m2·s)的某一热流密度下时(如图5(a)所示), 恶化与Bu和Ac准则不存在一定的对应关系. 随着质量流速进一步增大到1000 kg/(m2·s), 发现Ac数的最大值与壁面温度最值发生了对应关系, 浮升力完全不能预测恶化起始点的发生, 如图5(b) 所示. 因此, 有研究者认为流动加速是导致传热恶化的主要原因[15], 实际上,这个准则大小仅与流体的主流温度有关, 对于大质量流速, 当流体主流温度在拟临界点附近时,Ac数通常就会达到最大值.

图4 局部壁温 Tw,in, Bu, Ac 随主流焓值 ib的分布关系(a) P = 8.220 MPa, G = 200 kg/(m2·s), qw = 60 kW/m2,(b) P = 8.220 MPa, G = 520 kg/(m2·s), qw = 42 kW/m2Fig. 4. Local inner wall (Tw, in), Bu, Ac distributions with bulk fluid enthalpy (ib): (a) P = 8.220 MPa,G =200 kg/(m2·s), qw = 60 kW/m2, (b) P = 8.220 MPa, G =520 kg/(m2·s), qw = 42 kW/m2.

图5 局部壁温 Tw,in, Bu, Ac 随主流焓值 ib 的分布关系(a) P = 8.220 MPa, G = 700 kg/(m2·s), qw = 245 kW/m2,(b) P = 8.220 MPa, G = 1000 kg/(m2·s), qw = 245 kW/m2Fig. 5. Local inner wall (Tw, in), Bu, Ac distributions with bulk fluid enthalpy (ib): (a) P = 8.220 MPa,G =700 kg/(m2·s), qw = 245 kW/m2, (b) P = 8.220 MPa, G =1000 kg/(m2·s), qw = 245 kW/m2.

图6 给出了不同质量流速下的壁温、Bu和Ac随焓值分布, 其中▲和□分别对应质量流速为1500和 1250 kg/(m2·s). 图6(a)给出了两者在相同压力不同热流密度下的壁温分布, 质量流速为1500 kg/(m2·s)在靠近拟临界区发生了局部恶化,壁面温度出现峰值, 然后传热恢复, 其中a为壁温恶化起始点, b为壁温恶化完全恢复点; 质量流速为 1250 kg/(m2·s)为非恶化工况, 壁面温度光滑、连续地上升, a'点、b'点分别为与 1500 kg/(m2·s)质量流速相同主流温度下对应的特征壁温. 两个质量流速下的浮升力准则值沿主流温度分布如图6(b)所示, 相同的主流温度下正常传热的Bu值始终大于恶化传热, 在恶化起始点处, 大质量流速的Bu值是小质量流速的二倍, 根据Hall和Jackson[8]的理论, 此时浮升力对大质量流速的剪切力改变要大于小质量流速剪切力的改变, 传热恶化应该更严重, 但是结果恰恰相反. 流动加速的影响如图6(c)所示,Ac在靠近拟临界焓值附近同时达到了最大值, 整个过程二者在数值上并没有大的差距, 尤其是质量流速为1500 kg/(m2·s)恶化起始点对应的主流温度下, 二者的Ac分别为 2.26 × 10–11和2.15 × 10–11, 质量流速为 1500 kg/(m2·s)稍大一点, 但是不足以构成不同寻常传热行为的决定因素.

Kurganov和Kaptilnyi[17]在相对较大的质量流速的实验研究过程中也发现了Hall和Jackson[8]的浮升力模型不适用, 但是他们认为该浮升力模型假定在y+= 30时,t/tw< 0.9, 浮升力和加速度的影响是明显的, 这过高地估计了壁面热流密度,作者认为t/tw< 0.97 是合理的. 实际上, 在大质量流速下, 管道横截面依然会出现较大的密度梯度. 图7给出了所有数据的Gr和Re2.7的分布关系, □代表恶化传热, ▲代表正常传热. Hall和Jackson[8]认为当Gr>Re2.7时, 浮升力会减小剪切力可能会导致传热发生恶化, 并从理论推导出一个临界值 1 × 10–5. 从图7 可以看出, 这个临界值并不适用于所有数据, 仅凭Gr和Re2.7的相对大小不能区分出正常和恶化. 除此之外, 正常传热和恶化传热也难以区分, 这在其他的文献也被报道过[21].

图6 不同质量流速下 (a)局部壁温 Tw,in, (b) Bu, (c) Ac随主流焓值 ib的分布 (NHT, 正常传热; HTD, 恶化传热)Fig. 6. (a) Local inner wall Tw, in, (b) Bu, (c) Ac distributions with bulk fluid enthalpy ib (NHT, normal heat transfer; HTD, heat transfer deterioration).

图7 不同质量流速下的 Gr 和 Re2.7 分布关系Fig. 7. Gr and Re2.7 distribution at different mass flow rates.

4.2 拟沸腾理论和拟沸腾数决定传热过程

将超临界传热类比为亚临界压力下的过冷沸腾, 对于亚临界沸腾过程, 气泡在壁面上成核并生长, 这些气泡脱离壁面与核心区的过冷液体融合并凝结, 因此, 蒸汽不会在管道的核心区中积聚. 相反, 如果气泡不离开壁面, 而是合并形成一个蒸汽层, 结果会导致相应的壁面温度急剧升高.Kandlikar[22]认为蒸汽层的形成与一个重要的无量纲参数K1有关, 对于气泡受到的蒸发动量与惯性力的大小由K1表征:

其中Bo为沸腾数,Bo=qw/(Gifg),ifg为潜热,rf和rg分布为液相密度和汽相密度. 对于给定工质的密度比, 沸腾数对确定传热机理具有重要意义. 较大的K1值表示蒸汽充分膨胀, 使气泡附着在壁上, 当气泡的结合形成了一个蒸汽层, 导致壁面温度的飞升, 称为传热恶化. 现在我们类比亚临界沸腾过程, 来处理超临界流体的异常行为, 如图8(a)和图8(b)所示, 当壁面的温度达到拟临界温度, 类汽流体首先在壁面附近产生, 并且类液态流体逐渐向类汽态流体转变. 此时,K1数变成

对于亚临界压力下的沸腾,qw/ifg表示汽液界面上的蒸发质量通量, 这个蒸发过程满足汽、液平衡条件, 对于特定的工质, 潜热ifg仅与压力有关.但是, 对于超临界流体的类液与类汽转化过程, 亚临界下的潜热的概念不再适用, 从类液到类汽的相变焓由∆i量化, ∆i被认为与拟临界焓值ipc有关,因此, 在WL上的由类液态向类汽态转变的蒸发通量可以由qw/kipc来衡量, 为了方便计算,K1数最后变成

K1值的大小代表壁面蒸汽层的膨胀力与惯性力对比, 其中qw/Gipc为超临界沸腾数(SBO)数或拟沸腾数,ipc为给定压力下的拟临界焓值. 根据Banuti[20]的拟沸腾理论, 超临界流体跨过“相变点”吸收的能量, 不仅用于克服分子间的引力, 而且还要提供流体本身温度的升高, 拟沸腾温度随着压力升高而增大, 而ipc是流体本身的物理属性,随着压力增大而增大, 也就是拟沸腾相变焓随着压力增大而增大, 这与Banuti[20]的拟沸腾理论一致.

图8 基于拟沸腾理论的超临界流体径向膨胀模型Fig. 8. Radial expansion model of supercritical fluids based on pseudo-boiling.

图9 超临界沸腾数区分两个传热机制Fig. 9. Supercritical boiling number distinguishes the two regimes of heat transfer.

为了检验SBO数的有效性, 通过实验获得了8 mm管径下大量S-CO2流动传热实验数据, 共89个工况, 参数范围涵盖100—1500 kg/(m2·s),qw= 20—400 kW/m2,pin= 8.220—20.822 MPa.这些数据对应的壁温恶化和SBO的分布关系如图9所示, 其中红色代表恶化工况, 黑色代表正常工况. 几乎所有正常和恶化工况均可以通过5.204 ×10–4这一临界值区分开. 图中 A 点工况为 8 mm管径下正常传热外壁温随ib的分布, 没有明显的壁面温度飞升, 壁面温度光滑分布, 其中SBO =3.685 × 0–4, 这个值小于临界值, 而 B 点对应的工况明显壁面温度飞升, 为恶化传热, 其中∆T= 103 ℃为壁面温度飞升值, 其中 SBO = 1.012 × 10–3, 这个值大于临界值.

5 结论

采用传统的单相流体概念对超临界流体的流动传热过程进行研究已经持续了近80年了, 这些研究大都是在常物性的基础上, 做了大量假设, 考虑浮升力和流动加速效应对其流动传热的影响, 得出了大量的无量纲准则来判断其发生传热恶化的起始点, 但是这些准则不具有通用性, 单相的理想流体概念很难处理超临界流体的异常流动传热过程. 关于超临界流体有关密度的相变过程一直都是争议的焦点, 本文采用拟沸腾理论处理超临界流体的异常流动传热过程, 通过类比超临界压力与亚临界压力之间的传热, 提出超临界沸腾数SBO, 对于8 mm管径超临界二氧化碳从正常传热向恶化传热转变的临界值为 5.204 × 10–4, 在一定条件下,这个无量纲数表示低密度层膨胀引起的动量力与惯性力对比, 成功地采用拟沸腾数区分了超临界二氧化碳正常传热和恶化传热的边界, 为理解超临界流体的流动传热机理开辟了新的途径. 除此之外,我们后续的工作会继续探讨管径对超临界二氧化碳传热的影响, 关注超临界状态下类汽和类液相互转变的微观和宏观过程的关联, 尤其是这种转变过程中的相变焓如何确定, 我们将继续采用分子动力学模拟对其进行研究.