竖直环形通道内的湍流混合对流数值分析

胡相杰,温济铭,高璞珍,田瑞峰,*

(1.哈尔滨工程大学 核科学与技术学院,黑龙江省核动力装置性能与设备重点实验室,黑龙江 哈尔滨 150001;2.中国原子能科学研究院 核安全研究所,北京 102413)

以空气为传热介质、大气环境为热阱的反应堆舱室冷却系统具有系统简单、结构紧凑、不依赖冷却水、可靠性高等特点,可以用于小型反应堆的余热导出[1]。压力容器与舱室壁构成环形间隙,间隙内空气的流动受到风机强迫驱动以及密度差自然循环驱动的共同作用,处于混合对流状态,具有特殊的传热特性,对系统设备的安全运行有重要影响[2]。

Jackson等[3]是最早系统研究混合对流现象的学者之一,指出自下而上的湍流混合对流,在特定的工况范围内可能出现对流传热被浮力效应抑制或强化的现象,并从边界层理论出发推导出了混合对流的传热关联式。关于该现象的发生机理,多位学者通过实验和数值模拟等方法给出了解释。Carr等[4]认为浮力效应使得流场结构发生明显畸变,同时湍流切应力减小,黏性底层厚度增加,导致传热被抑制;Aicher等[5]指出黏性底层边缘和主流之间的平均速度梯度降低导致湍流的产生率减小;Shome[6]认为浮力效应抑制了边界层内的湍流产生,导致湍流层化;He等[7]将浮力效应对传热的影响方式归纳为直接效应和间接效应两种机制。

直接数值模拟(DNS)方法可以精确模拟混合对流传热问题[8],但介于其高昂的计算成本,更多学者选择了基于雷诺时均方法开展研究。混合对流情况下近壁面流场、切应力和湍动能发生显著畸变,Cotton等[9]指出高雷诺数k-ε方法不适用于求解混合对流问题,因此采用了低雷诺数LS(Launder and Sharma)模型,该模型是在标准k-ε模型基础上增加阻尼函数发展而来,对近壁面流体层有更高的求解精度。Mikielewicz等[10]对比分析了多种湍流模型,其中低雷诺数LS模型与经验式的预测结果吻合最好。

现有研究中混合对流换热系数关联式的精度大多不高,难以在工程上推广应用,对于混合对流传热特性产生机制的研究尚不够完善。本文建立环形竖直通道内空气混合对流的稳态数值计算模型,从局部换热系数、速度场、湍动能以及无量纲数等角度开展混合对流传热特性分析。

1 几何模型

本文研究对象如图1所示,环形加热通道由两个竖直圆筒壁嵌套组成,间距为d。空气自下而上流经加热通道。以底部的圆心为坐标原点建立坐标系,竖直向上为x方向。通道内侧壁面的中部(L/5≤x≤6L/5,L为加热段长度)采用均匀壁温(Tw)方式加热,模拟压力容器外壁面,其余壁面均为绝热。加热面及与之相对壁面的热辐射发射率(ε)均为0.9,在图1中由加粗线表示,其余壁面发射率为0.1。通道两端分别增加了拓展计算域,边界条件为环境温度和压力(T0=300 K,p0=101.325 kPa),该计算域的作用是维持通道进出口的自由流动状态[11]。

给定通道入口正压头为pin,使得进出口形成压差(Δp=pin-p0)。压差作为强迫对流的外部驱动力,传热温差(ΔT=Tw-Tf,Tf为通道内气体平均温度)作为自然对流的内在驱动力,两者的强度分别由无量纲数Re和Gr表示:

(1)

(2)

式中:u为气体平均流速,m/s;D为当量直径,取两倍通道宽度d,m;μ为动力黏度,Pa·s;重力加速度g=9.8 m/s2;β为膨胀系数,1/K;空气密度ρ由理想气体状态方程计算:

(3)

式中,空气的气体常数Rair=287 J/(kg·K)。

混合对流和纯强迫对流的无量纲换热系数分别定义为NuM和NuF:

(4)

式中:h为换热系数,W/(m2·K);q为加热面热流密度,W/m2;λ为热导率,W/(m·K)。

2 数学模型

2.1 基本控制方程

稳态单相流动传热基本控制方程为:

(5)

(6)

(7)

式中:U为速度矢量,m/s;T为温度,K;cp为比定压热容,J/(kg·K);物性参数cp、λ、μ均由关于温度的多项式函数定义;S为应变率张量;ρ0为环境温度T0和压力p0对应的空气密度,kg/m3。

2.2 湍流模型

本文采用SSTk-ω湍流模型[12],该模型是可求解高、低雷诺数问题的普适性模型,结合了k-ε和k-ω模型两者的优势。

(8)

(9)

式中:k为湍流动能;ω为比耗散率;Гk、Гω分别为k和ω的有效扩散系数;Gb为浮力产生的湍动能;Gω为ω产生的湍动能;Yk、Yω分别为k和ω的湍流耗散;Sk、Sω为源项;Dω为阻尼交叉扩散项;Gk为平均速度梯度产生的湍动能。

(10)

(11)

式中,μt为湍流黏度,Pa·s。

3 网格划分及敏感性分析

采用ICEM进行网格划分,控制近壁面第1层网格的y+约等于1,以满足SSTk-ω模型对近壁面网格密度的要求。数值求解采用Coupled压力速度耦合算法,自然对流为体积力驱动的物理现象,因此压力离散方案选择body force weighted格式,能量残差收敛到10-6,连续性、湍流等残差收敛到10-4。以79万、115万、150万、203万和316万的网格进行网格无关性验证,提取不同网格的流量和加热面平均对流换热系数两个计算结果,如图2所示。结果表明,203万与316万网格计算得到的流量和平均对流换热系数相对误差分别为

图2 网格敏感性分析Fig.2 Grid sensitivity analysis

0.66%和0.34%。因此在同时考虑计算准确性和效率的基础上,本文选择203万的网格进行后续研究。

4 计算结果与讨论

4.1 传热抑制与强化

图3示出无量纲换热系数(NuM)在加热面上沿流动方向x的分布,对应工况A1～A5、B1～B5的主要参数列于表1,通道纵横比L/D=25。B组工况的外部驱动压头小于A组,因此Re较低。从工况1到工况5,逐步提高加热面温度,使得Gr增加,同时保持外部驱动压头不变。

表1 L/D=25时计算工况Table 1 Calculation case at L/D=25

图3 换热系数沿程分布Fig.3 Local heat transfer coefficient profile

湍流对流换热系数在沿程上的分布一般包括下降段和水平段两部分,分别对应热入口段和充分发展段,如图3所示。热入口段的长度与Re正相关,因此A组的热入口段长度显著大于B组。

在高Re情况下,从工况A1～A5,随Gr的增加NuM减小,混合对流传热受到抑制,如图3a所示;在低Re情况下,从工况B1～B5,随Gr的增加NuM先减小,在B2工况达到极小值55.5,随后再增加,混合对流传热逐渐恢复甚至增强,这是由于充分发展段的自然对流湍流局部增强,使得传热强化,如图3b所示。

在纯强迫对流的理想情况下,即排除重力场引起的浮力效应后,随着加热面温度升高,强迫对流换热系数NuF在近壁面黏性效应的作用下降低[14]。前文分析中工况A1～A5传热抑制的成因包括了黏性效应以及特殊的混合对流传热机制,因此引入无量纲换热系数比NuM/NuF以评估混合对流换热系数相对于纯强迫对流的抑制程度。NuF由强迫对流的经典关联式Gnielinski公式[14]计算得到,包含了黏性效应修正:

(12)

式中:f为管内湍流的Darcy阻力系数;Pr为0.7;Re与对应NuM的Re保持一致。

浮力数Bo表征Re和Gr的相对大小:

(13)

考虑30种运行工况,包络系统所有可能的运行状态,并计算分析各个工况Bo和NuM/NuF的关系,如图4所示。当Bo趋于0时,Re极大,此时混合对流传热受强迫对流主导,自然对流作用可以忽略,NuM/NuF趋于1。随着Bo增加,自然对流的相对强度增加,NuM/NuF减小,在Bo=10-5左右达到最小值,此时传热抑制作用最强烈,NuM/NuF约为0.9,即浮力效应使得混合对流换热系数相比纯强迫对流减少了10%左右。随着Bo进一步增加,传热抑制作用开始减弱,传热恢复,NuM/NuF上升,当Bo达到3×10-5左右,NuM/NuF=1。随后Bo超过1,浮力效应开始在混合对流中起传热强化的作用。

图4 各工况NuM/NuF与Bo的关系Fig.4 Relationship between NuM/NuF and Bo

现有的NuM/NuF与Bo关联式拟合多采用Jackson等的形式[3],如式(14)所示,但该式在NuM/NuF的极小值附近拟合精度并不理想。

(14)

式中,C=104,n=0.46。

本文采用Symolon等的形式[15]对图4中的数据进行最小二乘法拟合,获得了NuM/NuF和Bo的关联式,如式(15)所示,R2=0.910 2,适用于8.89×10-7≤Bo≤5.99×10-5。

(15)

将研究计算拓展到其他的通道纵横比,并对结果进行拟合。当L/D=16.7、7.72×10-7≤Bo≤1.02×10-4时,拟合结果如式(16)所示,R2=0.967 4。

(16)

当L/D=50、1.08×10-6≤Bo≤2.17×10-5时,拟合结果如式(17)所示,R2=0.739 8。

(17)

不同纵横比的传热关联式汇总如图5所示。

图5 传热关联式汇总Fig.5 Summary of heat transfer correlation

根据图5中NuM/NuF和Bo的关系,可将混合对流现象划分为3个区间。当Bo≤1.3×10-6时,NuM/NuF接近1,该区间为强迫对流主导,自然对流作用可以忽略,因此称为强迫对流区;当1.3×10-63.3×10-5时,NuM/NuF大于1,热浮力起到强化传热作用,该区间称为自然对流传热强化区。Bo在定义上与几何参数无关,如式(12),因此混合对流传热特性的分区结论不仅限于本文的3个几何参数,也可推广到其他宽度的流道结构。

Symolon等[15]进行了竖直圆管内的混合对流实验,Shin等[16]进行了竖直窄矩形通道内的混合对流实验,两者的通道纵横比分别为90.3和100,分别拟合得到了混合对流传热关联式,在图5中与本文的拟合结果进行了对比。NuM/NuF的极小值代表了传热抑制的最恶劣工况。纵横比为16.7、25和50时,NuM/NuF最小值分别为0.927、0.882和0.782。当纵横比达到100,NuM/NuF最小值达到0.5左右。随着纵横比的增加,通道宽度变窄,曲线在图5中向下偏移,即传热抑制作用增强。另外可以观察到,NuM/NuF达到极小值的工况均出现在Bo约为1×10-5时,此时传热抑制作用最强烈。

4.2 速度场与湍流分布

图6～8分别为气体流速、湍动能k和雷诺应力τt在垂直于壁面方向上的分布。图6的横坐标为无量纲y+,纵坐标为无量纲流速u+,分别由式(18)、(19)定义:

图6 半对数坐标下的无量纲速度分布Fig.6 Dimensionless velocity profile on semi-logarithmic coordinate

(18)

(19)

式中:uτ为近壁面摩擦速度,m/s;y为与壁面的距离,m;τω为壁面切应力,Pa。

图7、8横坐标为无量纲径向距离R,R=0代表加热面,R=0.5代表主流中心。

图7 湍动能分布Fig.7 Turbulent kinetic energy profile

沿垂直于壁面方向可将流体层依次划分为黏性底层、过渡层、对数率层和主流区。y+≤5为黏性力主导的黏性底层。过渡层为边界层内湍动能的主要产生区域,由图7、8可见,湍动能和雷诺应力的峰值大致位于R=0.02处,对应y+=40,实际离壁距离4 mm,因此将5

如图6a所示,工况A1由强迫对流主导,速度分布形状呈半倒U型。从工况A1～A5可知,随着Gr增加,近壁面气体温度升高,温度的不均匀使得边界层内的气体浮升力大于主流区,形成非均匀的体积力,使得边界层气体加速。过渡层内的速度梯度仍较大,而对数率层和主流区的速度分布先逐渐趋于平坦,而后出现较小的速度峰,向自然对流主导的半M型过渡。与此同时,局部速度梯度的减小使得雷诺应力降低,阻碍了湍动能的产生,如图7a、图8a所示,导致湍动能显著降低和湍流层化,混合对流传热相比纯强迫对流受到抑制,NuM/NuF从0.99减小到0.88。

图8 雷诺应力分布Fig.8 Reynolds stress profile

从工况B1～B5可知,随着Gr增加,速度峰进一步升高,速度分布逐渐趋向于自然对流主导的半M型分布,如图6b所示。边界层外缘和主流中心之间的负速度梯度导致主流区出现了负雷诺应力,如图8b,负雷诺应力提供了额外的湍动能,在主流区形成了另一个湍动能峰,与位于过渡层外缘的湍动能峰共同构成了双峰结构,如图7b所示。He等[7]也在这一过程中观察到了湍动能的双峰结构。主流区的负雷诺应力绝对值远大于边界层内的正应力,相比前述湍流层化的状态,通道内的湍动能得到恢复甚至增强,NuM/NuF从0.90增加到1.15,混合对流传热相比纯强迫对流得到加强,浮升力起到了强化传热的作用。

浮升力通过间接效应影响混合对流传热,即首先改变平均速度场,使得雷诺应力变化,影响了局部湍动能的产生,进而导致湍流的层化或恢复,影响混合对流传热。相比较宽的通道,窄通道内的Re更小,因此在发生湍流层化时,局部湍动能更低,传热抑制更为恶劣,如图5所示。

5 结论

本文基于CFD方法对竖直环形通道内的混合对流传热问题开展了数值模拟研究,主要结论如下。

1) 拟合得到了Symolon等形式的混合对流传热关联式。定义Bo=Gr/Re3,基于Bo将混合对流现象划分为3个区域:Bo≤1.3×10-6为强迫对流区;1.3×10-63.3×10-5为自然对流传热强化区。该分区方法对于不同纵横比的竖直通道具有一定的推广意义,可用于验证混合对流换热系统的运行可靠性。

2) 窄通道内发生湍流层化时的传热抑制作用更强,该结论有助于工程上合理选择传热通道的几何尺寸。

3) 浮升力对竖直通道内混合对流传热的影响机制可以归纳为间接效应,即浮升力使平均速度场发生畸变,使得雷诺应力变化,进而影响湍动能的产生,导致湍流的层化或恢复,以及传热被抑制或强化。