本文的研究对象包括推力室热燃气区域与再生冷却区域两部分,图1为计算所采用的推力室示意图,为了简化计算,采用非平衡流模型,不考虑推进剂的喷注与燃烧过程,已充分燃烧的热燃气从推力室入口进入推力室,随后从喷管出口流出。冷却剂从喷管扩张段处的冷却通道入口进入,流经通道对壁面进行冷却之后从头部附近的出口流出。其中推力室部分的尺寸为,推力室圆柱段长度L₁=100 mm,收敛段长度L₂=40 mm,扩张段长度L₃=50 mm。推力室圆柱段直径D_c=52 mm,喉部直径D_t=26 mm,喷管出口直径D_o=56 mm。图2给出了再生冷却区域几何结构以及壁面纵向肋示意图,其中冷却通道部分尺寸为矩形通道高5 mm,宽3 mm,内壁厚1 mm,在内壁面靠近燃气侧设置均匀分布的矩形纵向肋。

图1 仿真所采用的推力室示意图
Fig.1 Schematic diagram of thrust chamber used for simulation

图2 再生冷却区域几何结构及壁面纵向肋示意图
Fig.2 Schematic diagram of regeneration coolingregion and wall longitudinal ribs

首先,在推力室圆柱段热燃气侧壁面设置均匀分布的纵向肋,其中沿周向肋总数目为120条,肋高、宽、间距均为0.68mm。此外,在冷却通道内部底面设置均匀分布的人工粗糙度,其凸台高、宽均为0.3mm,凸台间距为10 mm。为了比较,本文共计算4种工况如图3所示,其中Case 1为参考工况,无冷却通道底部的人工粗糙度与燃气侧壁面上的纵向肋,Case 2与Case 3工况分别设置有冷却通道底部人工粗糙度与燃气侧壁面上的纵向肋,而Case 4工况则同时设置有人工粗糙度与纵向肋。

图3 不同工况下几何模型设置
Fig.3 Different geometry setting under four cases

由于几何体的对称性,只取原几何的1/40即9^o作为计算域,图4给出了再生冷却通道区域的计算域与网格划分,整体网格划分全部采用结构化网格由于其高计算精度与计算效率。在凸台与纵向肋附近网格加密以捕捉复杂的流动特征。Case 1工况下热燃气区域与冷却区域网格总数分别为21万与33万,其余工况下由于在肋与凸台附近加密因此网格总数略高于Case 1。

图4 再生冷却区域计算域与网格
Fig.4 Calculation domain and mesh generation of the regenerative cooling region

本文计算区域划分为热燃气区域与再生冷却区域两个独立的部分,并采用文献^[5]中的耦合方式来处理这两部分之间的耦合传热,耦合方式如下:

1)在热燃气区域壁面施加一个初始温度分布并计算直至收敛,在此基础上,得到了热燃气侧壁面的热流密度分布;

2)以第1步得到的热流密度分布作为边界条件,对再生冷却区域进行了传热计算,得到一个新的热燃气侧壁温度分布;

3)以第2步得到的壁面温度分布作为边界条件,再次对热燃气区域进行了数值模拟,得到了一个新的壁面热流密度分布;

4)重复步骤2和步骤3,直至热气侧壁面热流密度和壁面温度在两次迭代前后之间的误差小于1%即可认为计算已经收敛。

本文采用商业软件Fluent 17.2作为求解器进行计算,湍流模型采用重整化群(RNG)k-SymboleA@湍流模型,该模型中包含了涡流对湍流的影响,因此能够处理冷却通道中超临界流体复杂的湍流流动^[16,17],近壁区域处理采用标准壁面函数。燃气区采用有限速率14组分20步化学反应模型来模拟燃气的非平衡流动^[18],其中采用涡耗散概念(EDC)模型来考虑化学反应模型与湍流模型之间的相互作用,它能够包含比较完整的化学反应的相关机理^[19]。

入口边界取质量流量入口,其中热燃气入口质量流量为2 kg/s,温度为3498 K,其中温度与各组分摩尔分数由一维化学平衡计算软件CEA得到。冷却剂入口质量流量为5 kg/s,温度为120 K。

出口边界取压力出口,其中喷管出口压力设置为101 325 Pa,冷却通道出口压力设置为8 MPa。计算域对称面取对称边界条件,流固耦合壁面取无滑移条件,壁面粗糙度取3.2 μm,除加热面和流固耦合壁面以外其余壁面均取绝热边界。

本文应用Soave-Redlich- Kwong(SRK)状态方程来求解热燃气区域的流体密度,流体的各项输运性质,包括粘度和导热系数等利用Chung方法计算^[20]。冷却通道区域内壁材料取铜,外壁材料取镍,其各项物理性质均取常数。冷却剂为甲烷,综合考虑计算精度与计算成本两方面的因素,其各项热物理性质采用文献^[5]中的方法来处理。甲烷在各个压力下的密度、定压比热、粘性和热导率等随温度的变化数据由National Institute of Standards and Technology(NIST)得到,再利用二维插值的方法使用用户自定义函数(UDF)写入到Fluent中进行计算。在本研究中所取压力范围为6～20 MPa,温度为120～600 K。

表1给出了推力室热燃气区域、冷却剂、内壁与外壁4个区域3种网格密度下x× y × z方向的节点数设置,图5为三种网格数目下Case 1热燃气侧壁面温度沿轴向变化。由图可知,粗网格计算结果与其他两种网格下的结果相比存在较大区别,而后两种网格数目所计算得到的结果差别在10 K以内。因此综合考虑计算精度与计算成本,选择中间数目的网格来进行后续的研究。

表1 网格无关性验证所用不同网格设置
Tab.1 Various grid arrangements for grid independence test

图5 三种网格密度下Case 1壁面温度沿流向变化

为了研究在圆柱段添加纵向肋之后对推力室内部流动与传热的影响,在z=0.05 m处取一个截面,图6和图7分别给出了有无从纵向肋即Case 1与Case 3两种工况下该截面上热燃气区域与冷却通道区域的温度分布。当在壁面设置纵向肋之后,肋顶处贴近热燃气部分的温度分布与无肋工况下基本一致,而两条肋之间区域处的燃气温度明显较低。因此,此处的传热会受到影响。而由图7可知,虽然肋间燃气温度较低,但是由于增大了换热面积,因此有肋工况下壁面温度高于无肋工况,且整体壁面温度梯度也有所增大。

图6 2种工况下推力室z=0.05 m截面上温度分布

图7 2种工况下冷却通道z=0.05 m截面上温度分布

图8和图9分别给出了4种工况下z=0.05 m截面处壁面热流密度与温度沿周向变化,通过比较Case 1与Case 3工况下壁面热流密度与温度变化可知,正如前文所言,肋顶面处的热流密度与无肋工况下的基本一致,而肋侧面与肋底面处的热流密度由于该处燃气温度较低因而明显偏低。而通过比较两种工况下壁面温度变化可知,带肋工况下各个面上的温度均高于无肋工况,这是由于虽然肋侧面与肋底面热流密度较低,但由于整体传热的提高使得靠近通道底面处的冷却剂温度升高,削弱了其冷却效果,导致壁面温度显著升高。

图8 4种工况下z = 0.05 m处壁面热流密度沿周向变化

图9 4种工况下z = 0.05 m处壁面温度沿周向变化

图 10给出了Case 1与Case 2两种工况下圆柱段冷却通道中心截面上近壁区域流线图,由图可知,当通道底部无人工粗糙度时,流线光滑而均匀,而当设置凸台之后,靠近凸台处的流动发生了扰动,在凸台下游处产生了漩涡。

图 10 2种工况下通道底部近壁区域流线图

为了定量分析人工粗糙度所产生的涡流强度的变化规律,在这里给出螺旋度的定义为:

(1)

式中:H为螺旋度; V为速度矢量。涡流的强度可以通过螺旋度的绝对值来测量,图 11给出了2种工况下沿流动方向螺旋度变化。无人工粗糙度工况下,由于流线光滑而均匀,因此螺旋度基本为0,而设置人工粗糙度之后使得螺旋度有了显著的提升。

图 11 2种情况下螺旋度沿流动方向变化

螺旋度的提升会使得流体之间的传热得到改善,图 12给出了Case 1和Case 2冷却通道z=0.05 m截面上温度分布,设置凸台能够使得靠近通道底面凸台附近冷却剂温度梯度减小,促进了流体之间的传热,进而使得壁面温度也有了明显的降低。而由图8和图9可知,无论壁面是否有纵向肋,添加人工粗糙度均可使壁面热流密度升高而壁面温度显著降低。

图 12 Case 1和Case 2冷却通道z=0.05 m截面上温度分布

图 13给出了4种工况下壁面平均热流密度与温度沿轴向变化,由于计算省去了推进剂的喷注与燃烧过程,因此使得圆柱段靠近喷注面板附近的温度高于实际情况下的壁面温度。由图8可知,由于肋侧面与肋底面热流密度均低于无肋工况,因此有纵向肋的两种工况下圆柱段壁面平均热流密度均低于无肋工况。此外,通过在通道底部设置人工粗糙度,使得该处的热流密度有所提升。四种工况下收敛段与扩张段由于几何构型完全一致,因此该区域的热流密度基本一致。而通过比较壁面温度变化可知,带有纵向肋工况下虽然平均热流密度较低,但是由于换热面积大大增加,因此使得圆柱段温度急剧升高,几乎达到了无肋工况下的两倍。通过设置人工粗糙度可以使壁面温度整体下降约100 K,而在实际应用中可以通过优化凸台尺寸与布局使得该效果进一步增强。

图 13 4种工况下平均热流密度与壁面温度随流向变化
Fig.13 Streamwise variation of average wall heat flux and temperature under four cases

由前文分析可知,由于在推力室圆柱段设置纵向肋导致换热面积增加,使得平均热流密度与壁面温度分布规律产生差异,表明此时单纯使用沿周向平均热流密度无法准确的描述发动机壁面的实际换热情况。在这里引入等效平均热流密度的概念,其中等效热流密度q_eq,ave的计算公式为:

q_eq,ave=(Q_rib)/(S_smooth)(2)

式中,Q_rib、S_soomth分别代表轴向该位置处带肋壁面沿周向总热流量与光滑壁面表面积。

图 14给出了4种工况下壁面等效平均热流密度沿轴向变化,通过对比图 13(b)与图 14可知,有无纵向肋工况下壁面热流密度变化规律与壁面温度保持一致,表明使用该参数能够描述带肋发动机推力室壁面的实际换热特征。此外,圆柱段带肋工况下在z=0.1m处热流密度与温度均存在一个突降过程,这是由于该位置处为纵向肋的末端,即推力室圆柱段与收敛段的交界处。在此位置处由于肋的消失使得在扩张段换热面积突然减小,因此壁面温度也随之降低。

图 14 4种工况下等效平均热流密度沿流向变化
Fig.14 Streamwise variation of equivalent-average wall heat flux under four cases

为了定量研究添加纵向肋与人工粗糙度对发动机流动与传热的影响,表2给出了4种工况下冷却剂压降,温升及圆柱段壁面平均温度比较,由表2可知,相较于Case 1,Case 2壁面温度降低了85.4 K,压降增大了0.11 MPa,同样,Case 4相较于Case 3,壁面温度降低了66.6 K,而压降增大了0.11 MPa。这表明,在通道底面添加人工粗糙度可以有效降低壁面温度,但同时会使压降产生明显的升高。此外,而相比于Case 1,Case 3的冷却剂温升提高了24.2 K,但同时圆柱段壁面平均温度升高了276.4 K。这表明在推力室燃气侧壁面添加纵向肋可以有效提高冷却剂温升,但同时也会使该处的壁面温度急剧升高,因此需要额外注意该处热防护工作。

表2 4种工况下冷却压降、温升及圆柱段壁面平均温度比较
Tab.2 Comparison of coolant pressure drop, temperature rise and average wall temperature of cylinder section under four cases

值得一提的时,Case 2和Case 4相较于Case 1和Case 3,冷却剂温升分别升高了1.3 K和降低了4.7K。这表明,通过在通道内底部设置凸台,虽然由于扰流与增大换热面积使得靠近底面处的传热增加,但并不一定能会使得冷却剂整体温升提高。图 15给出了Case 1和Case 2工况下冷却通道出口处冷却剂温度分布,由图可知,Case 2工况下由于凸台的存在,通过扰流与增大换热面积,使得由于非对称加热导致的温度分层有所削弱。但同时由于壁面温度较低,使得靠近通道底部处的流体温度明显较低,因此两种工况下冷却剂整体温升比较接近。

图 15 Case 1和Case 2冷却通道出口处冷却剂温度分布
Fig.15 Coolant temperature distribution at cooling channel outlet under Case 1 and Case 2