再生冷却推力室准二维传热数值计算

2023-05-05吴有亮丁煜朔刘阳旻

火箭推进 2023年2期

吴有亮,丁煜朔,刘潇,刘阳旻,田原

(北京航天动力研究所,北京 100076)

0 引言

液体火箭发动机燃烧室内充满着高温、高压以及高速流动的燃气,燃气将对壁面产生巨大的对流和辐射热流,若不采取有效的冷却措施,壁温将高到现有工程材料无法承受的程度,进而产生灾难性的后果[1-2]。因此,火箭发动机传热和热防护是发动机设计中的一个突出问题。现代液体火箭发动机中,再生冷却是应用最为广泛的冷却技术。

再生冷却传热过程涉及燃烧、流动和换热等过程之间的耦合相互作用,导致传热计算相当复杂。许多学者利用商用CFD软件对推力室传热过程进行三维仿真计算[3-7]。吴峰等应用经验公式计算燃气侧对流、辐射换热,采用气—固耦合算法进行冷却通道—冷却剂三维耦合传热计算[8]。康玉东等考虑冷却剂温度分层,应用经验公式计算燃气侧对流、辐射换热,采用气—固耦合算法进行冷却通道—冷却剂三维耦合传热计算[9]。Pizzarelli对高深宽比再生冷却通道进行了三维流固耦合传热数值仿真,得到了较为精确的计算结果[10]。Divalentin则利用Fluent研究了通道曲率变化引起的二次流对传热过程的影响[11]。CFD计算虽然能够获得更为精确的计算结果,但是计算过程复杂,收敛性差,耗时长且对计算机配置要求高,不利于推力室冷却结构的快速优化设计。

目前工程上传热和冷却剂流阻计算主要利用一维的传热计算模型[12-14],将传热过程简化为燃气与内壁面之间的对流和辐射换热、通过推力室内壁的热传导以及再生冷却剂与推力室壁之间的对流换热3个过程。一维传热计算模型相对实际传热过程做了大量简化,忽略了冷却通道内非对称加热导致的热分层现象,随着冷却通道深宽比的增大,径向温度分层愈发明显,远离燃气侧壁面的冷却剂保持着较低温度,仍可有效地吸收热量。采用一维传热模型,无法准确描述远端冷却剂的冷却能力,导致计算结果与试验测量结果存在较大误差。燃气侧对流换热采用经典的Bartz公式,未考虑近壁面边界层内的导热过程。

为了提高传热计算精度,本文参考文献[15]的冷却剂传热方程,建立了再生冷却推力室的准二维传热数值模型,考虑了冷却剂层间导热导致的温度分层效应,燃气侧增加了直接求解边界层控制方程得到热流密度的方法。最终基于MATLAB开发完成了通用的再生冷却推力室准二维传热程序,较为准确地计算了某氢氧发动机再生冷却推力室的传热情况,并与一维和三维传热计算结果以及热试验结果进行了对比,验证了模型的准确性。

1 计算方法

再生冷却通道的剖面结构如图1所示,准二维传热模型中传热过程主要包括燃气与内壁面之间的对流和辐射换热、内壁导热、冷却剂沿径向的层间导热、冷却剂与内壁的对流换热。

图1 通道结构示意图Fig.1 Channel cross-section

1.1 燃气侧传热

针对燃气侧对流换热,本文研究了两种计算方法,一种是传统的Bartz公式,另一种是直接求解边界层控制方程得到热流密度。

1.1.1 传统Bartz公式

qgc=hg(Taw-Twg)

(1)

式中:hg为燃气侧的换热系数;Taw是燃气绝热壁温;Twg是燃气侧的壁温。

燃气侧对流换热系数用Bartz[16]提出的管内充分发展的湍流换热准则方程表示,即

(2)

式中:μ、cp、Pr均以总温(Tc)ns为定性温度;dt为燃烧室喉部直径;d为沿燃烧室轴线计算截面的直径;(pc)ns为燃烧室压力;c*为特征速度;rwt为喉部处喷管外形曲率半径。

通过附面层时气体性质变化的修正系数为

(3)

式中:Ma为计算截面的马赫数;k为燃气比热比。

1.1.2 直接求解边界层控制方程

连续方程为

(4)

动量方程为

(5)

能量方程为

(6)

式中:ρ为密度;u为流向速度;v为纵向速度;r为横曲率半径;H为总焓;T为温度;p为压力;μ为动力黏度;cp为定压比热容;上标“—”为雷诺平均;上标“′”为湍流波动量;N为流动指标数,N=0时是二维流动,N=1时是轴对称流动。

针对上述边界层控制方程进行柱坐标转换后,使用有限差分方法求解,并完成TDK(二维化学反应喷管流场)程序开发。具体求解过程可参考文献[17-18]。利用自编TDK程序求解燃气侧边界层控制方程后,直接求解不同恒壁温条件下的热流密度沿轴向位置分布,得到了一个燃气侧热流密度插值表。在传热计算模块中,这张热流密度表作为输入,通过对壁温和轴向位置插值计算不同位置和不同壁温下的热流。

1.2 冷却剂侧传热

冷却剂侧考虑传热和摩擦的影响,基于质量、动量和能量的稳态守恒定律,建立了冷却通道内冷却剂流动模型。其中质量和动量方程以一维形式描述,而由于非对称加热的冷却通道中,在径向会产生显著的温度分层现象,因此流体能量方程以轴向和径向二维形式描述。由于冷却通道内的传热模型一维和二维形式并存,因此本模型称为准二维传热模型。

由于动量方程为一维形式,因此冷却剂压力沿径向相同,而能量方程为二维形式,所以温度、流速沿轴向(x方向)、径向(y方向)不同。冷却剂物性参数如密度、焓值、黏性等沿轴向、径向不同。

p=p(x),T=T(x,y),u=u(x,y)

其中速度可表示为

u(x,y)=umh(x)F(y)

(7)

式中:umh(x)表示通道中间位置速度沿轴向分布;F(y)表示速度剖面的形状。

根据一维质量守恒定律

(8)

可得到中间位置速度

(9)

式中h和b分别为冷却通道的高度和宽度。

速度剖面可表示为

(10)

式中:B为速度的不对称率;n为与雷诺数有关的修正系数。冷却通道内热流传递过程如图2所示。

图2 通道内热流传递示意图Fig.2 Heat fluxes in a slice of cooling channel

二维能量方程为

(11)

式中:ρ为冷却剂密度;H0为冷却剂总焓;qc(y)为冷却剂层间导热;qw为侧面肋与冷却剂的对流换热。

冷却剂层间热流为

(12)

式中kt为沿径向湍流导热系数,与冷却剂雷诺数Re和导热系数λ有关,这是文献[19]在光壁条件下提出的,kt=λ0.008Re0.9。

侧面肋与冷却剂的对流换热热流密度为

qw=hl(Twl-Tl)

(13)

式中:hl为冷却剂侧的对流换热系数;qw为冷却剂从壁面吸收的热流密度;Twl为液壁温;Tl为冷却剂温度。

对于液氢,其对流换热系数可根据Hess-Kunz公式表示为[20]

(14)

式中:Nu为努塞尔数;Re为雷诺数;Pr为普朗特数;ν为运动黏度。

冷却剂压力根据一维动量方程可知

pn=pn-1-Δpn-1,f-Δpn-1,m

(15)

式中:pn为计算截面冷却剂压力;pn-1为前一计算截面冷却剂压力;Δpn-1,f为黏性压力损失;Δpn-1,m为动量损失。

冷却剂状态方程为

ρ=ρ(p,T),H=H(p,T),k=k(p,T)

(16)

1.3 内壁导热

由于冷却通道内考虑了冷却剂温度分层,因此侧面肋上需要考虑径向的壁面温度分布。燃气侧热流通过内壁到达冷却剂和肋,进入侧肋的热流等于侧肋传向冷却剂的热流,其热流可根据稳态传热平衡方程表示为

(17)

式中:kw为内壁导热系数;tw为侧肋宽度。

进入侧肋的热流等于燃气侧对流换热热流密度,可表示为

(18)

式中δw为内壁厚度。

外壁视为绝热壁面,因此侧肋顶端导热为0,即

(19)

1.4 辐射换热

最大辐射热流可表示为

(20)

式中:qr为辐射热流密度;Tg为燃气静温;εw,ef为壁面有效黑度;σ为斯忒藩—波尔兹曼常数;εg为燃气黑度;αw为壁面吸收率;c为边区混合比修正系数,有边区混合比时取值0.65。

辐射热流沿推力室轴线变化,根据最大热流插值近似计算得到。

1.5 计算过程

根据上述的传热计算模型,使用MATLAB编写了再生冷却准二维方程程序,具体迭代求解过程如下:

1)将冷却通道沿轴向划分为N个截面,沿径向划分为M层,沿轴向依次计算每个截面上的参数分布,计算时物性参数直接调用Nist;

2)根据初始值或前一截面计算结果,给定当前截面的中间位置速度umh、冷却剂压力p和燃气侧内壁温度Twg;

3)主要根据能量方程(11)、侧肋导热方程(17)、状态方程(16)以及相关方程,由已知的umh、p、Twg,可求出冷却剂温度Tl、壁面温度Tw;

4)根据质量守恒方程(8)、冷却剂压力方程(15)、内壁导热方程(18),由上一步求出的Tl和Tw,可求出新的umh、p、Twg;

5)重复进行第3)、4)步的过程,迭代计算直到umh、p、Twg满足收敛要求,即前后两次计算差值小于0.001。

2 计算结果

2.1 计算模型

利用此程序计算了某型号氢氧发动机再生冷却通道内的传热情况,并与一维传热计算结果和Fluent三维传热计算结果进行了对比。冷却剂参数与推力室主要工作参数见表1,采用逆流冷却方式。

表1 推力室主要工作参数

三维计算中考虑到推力室的结构对称性,周向选择60°作为计算区域,网格划分如图3所示。冷却通道壁面粗糙度给定6.3 μm。计算过程为了提高收敛速度,采用燃气区域和再生冷却通道分区计算、互给边界条件的迭代方法。求解过程中,给定燃气壁面温度,求解燃气区控制方程直至收敛,得出燃气向壁面的热流密度,以此热流密度为边界条件,求解冷却剂与冷却通道的耦合流动换热控制方程直至收敛,得出气壁面温度。如此循环几次,如果前后两次气壁面热流密度和温度相差满足一定的精度,则认为计算收敛。为了使收敛加速,可以先用一维经验公式得出气壁面热流密度和温度的初步分布。

湍流模型采用标准k-ε模型,应用标准壁面函数法处理壁面物理量与近壁面区域物理量之间的相互联系,无量纲距离y+=30～300。冷却剂热物理性质和输运性质随压力和温度的不同而发生变化,采用自定义函数编写。

图3 网格模型Fig.3 Computational grids of CFD models

2.2 结果对比

图4为燃气侧壁面热流沿程分布,图中CFD_3D和1D分别为三维CFD和一维传热程序计算结果,TDK_2D和Bartz_2D均为准二维传热程序计算结果,只是燃气侧热流密度分别采用TDK程序和Bartz公式计算。不同计算方法中,一维传热得到的最大热流最小,准二维Bartz公式算得热流次之,二者得到的热流密度曲线变化趋势基本一致,在冷却通道的突扩突缩处,壁面热流有小幅度的波动。TDK程序和三维计算得到的最大热流较为接近,圆柱段以后二者曲线吻合得非常好。但是在靠近喷注面区域,三维计算得到的热流呈不断增大趋势,而TDK算得热流反而呈减小趋势。主要是因为三维计算考虑了燃烧过程,在喷注面附近为雾化混合区,燃烧未充分进行,热流密度较小。此外,喉部以后区域,一维和准二维Bartz公式算得的热流比三维和TDK程序计算结果大,可能是因为喉部以后区域内沿轴向马赫数不断增大,边界层逐渐增厚,冷却剂与内壁面之间的换热热阻增大,影响了传热过程,而Bartz公式未考虑附面层效应。

图4 热流沿轴向分布Fig.4 Variations of wall heat flux

造成热流密度差异的主要原因是三维CFD与TDK程序考虑了湍流和附面层的影响,能够较为精确地计算壁面温度梯度,并且在计算燃气物性中对壁面热流计算影响较大的密度和黏度时,CFD与TDK程序均使用了物理上较为精确的公式,而巴兹公式直接使用以经验公式表征的燃气密度与黏度。NASA也曾使用TDK耦合RTE进行推力室传热计算,其热流计算结果也与推力室CFD仿真结果吻合较好,且计算时间大大减少[21]。NASA的计算方法与国内使用的一维热平衡传热程序并无本质区别,主要区别在于燃气侧热流密度使用附面层程序求解而非巴兹公式。因此从原理上看,基于边界层控制方程直接求解燃气侧热流密度,对于提高工程上推力室燃气侧传热计算精度和效率是非常重要的。

图5为冷却通道中冷却剂平均温度沿程分布。在不同计算方法中,冷却剂入口温度均为35 K,在流动过程中平均温度不断增大。由于冷却剂温升主要受燃气侧总热流影响,在喉部以后区域,采用Bartz公式的一维和准二维程序算得冷却剂温度偏高。在喉部以前区域,由于热流相对较小,冷却剂温度增长放慢,但总的温升与TDK_2D较为接近,具体温升见表2。

图5 冷却剂温度沿轴向分布Fig.5 Variations of coolant temperature

表2 不同方法计算结果

图6为采用准二维传热程序计算的冷却剂温度分层情况,计算过程中冷却通道内沿径向划分了18层单元,第1层靠近内壁,第18层靠近外壁。

图6 冷却剂分层温度沿轴向分布Fig.6 Streamwise variations of coolant temperature between different layers

可以看到冷却剂在径向上存在明显的温度分层现象,靠近内壁的冷却剂温度最大,出口处可达233 K,靠近外壁的冷却剂温度最小,出口处为76 K。因此冷却剂的热力学性质在同一通道截面上也存在较大差异,本文在计算过程中直接调用Nist求解每层冷却剂的热力学性质,避免了不均匀物性场对换热的影响。

图7为三维CFD和TDK_2D计算的出口截面冷却剂温度分布云图,三维计算的温度分层更为明显,最大温度可达350 K。

图7 冷却剂出口温度分布Fig.7 Distribution of outlet temperature of coolant

图8为冷却剂平均压力沿轴向分布情况,各种方法计算得到的流阻在表2中给出。在经过突扩突缩处,冷却剂的总压发生不同程度的损失。因为通道的横截面积突然发生变化,冷却剂的流速大小发生急剧改变,造成了较大的局部损失,冷却剂的总压下降。

图8 冷却剂压力沿轴向分布Fig.8 Variations of coolant pressure

图9为CFD算得壁面温度分布云图,图10为燃气侧壁面温度沿程分布。各种方法计算得到最大气壁温在表2中给出,最大温度分别在-0.03、0、-0.07、-0.07 m处出现。由图中可以看出,TDK_2D和三维计算得到的壁温曲线在喉部以后较为接近,喉部以前TDK_2D计算壁温偏高。一维和Bartz_2D计算壁温曲线总体趋势较为一致,但是由于圆柱段以后温度分层影响更为显著,此时壁温比未考虑温度分层时高50 K左右。

图9 气壁温三维仿真结果Fig.9 Distribution of gas wall temperature for CFD_3D

图10 气壁温不同方法计算结果Fig.10 Variations of gas wall temperature

将计算得到的一些关键参数进行了对比,见表2。 4种方法计算的压降均较试验值偏小,主要是因为计算未考虑入口和出口处局部流阻损失。采用准二维传热程序计算的温升与试验值吻合得非常好,误差在10%以内,优于一维传热结果。CFD算得热流峰值最大,TDK_2D算得气壁温峰值最大,由于缺乏试验值比较,计算误差无法进行判断。本程序虽然计算时长相对一维计算偏大,但是远小于三维计算时间,且不用考虑前处理过程,计算方便。