浮升力对航空煤油传热流动影响的数值研究

2021-06-28姚长鑫

航空工程进展 2021年3期

姚长鑫

（重庆建筑工程职业学院城市轨道交通与机电工程系，重庆400072）

0 引言

超燃冲压发动机是助力飞行器实现超高声速飞行的核心部件。吸入发动机的高速空气被压缩后温度可达上千摄氏度，高温空气继而与燃料混合燃烧产生巨大推力，从而大幅提升飞行速度。2018年我国关于超燃冲压发动机的风洞实验取得重大突破，极大地推进了我国地空往返运输系统动力技术的革新。

超燃冲压发动机具有重量轻、载荷大的优势，但同时也存在难以攻克的技术难关。高温空气和燃烧后的高温气体同时造成燃烧室壁面温度大幅增加。当马赫数为6.5时，来流总温将超过1 800 K，燃烧后产生的气体温度甚至会高于2 800 K［1］。目前尚未找到合适的材料能在此高热载荷状态下长时间工作，因此对燃烧室的热防护不可或缺。再生冷却技术较为可靠高效，其基本原理是将发动机燃料作为冷却剂，在其流经冷却通道时将燃烧室壁面上过多的热量吸收，在冷却发动机壁面的同时起到预热燃料的作用。再生冷却技术优势显著，发展潜力大，近年来成为各国研究热点［2-7］。

再生冷却效果受到多种因素的影响，外因包括冷却结构、涂层材料、压力等［8-12］，超临界状态下物性参数骤变引起的冷却效果改变则属内因。目前对于超临界水、超临界CO2的研究较多，通过对超临界流体特性的研究［13-15］，一般认为航空煤油在超过超临界压力时物性突变，造成的近壁流体湍动能下降或浮升力引起的边界流动流化都可能引起传热恶化。近年来对于这一主题已有较多研究，Zhang Jingzhi等［15］通过数值模拟得出管壁温度超过拟临界温度或流体平均温度高于临界温度时管内超临界航空燃油会出现传热恶化的结论，可通过提升煤油压力实现对传热恶化的改善；王彦红等［16］对竖直管内的RP-3航空煤油的传热和流动进行了实验研究，认为当航空煤油边界层流体处于临界点附近，热物性剧变形成的浮升力和热加速效应造成了管内的传热恶化，当施加在管壁上的热流密度减小、进口压力提升或进口温度增加时，传热恶化现象得以缓解；李勋锋等［17］对圆管内航空煤油的流动与传热特性以及黏性底层内的传热机理进行了探讨，结果证明航空煤油在加热起始段的对流传热系数增长较快，但是当持续加热至壁温超过RP-3的拟临界温度时，对流传热系数呈现先小幅下降后显著升高的变化趋势。另外，RP-3本身的物性变化及湍流在临近壁面处湍流强度的大小都会对其传热特性产生影响。程泽源等［18］对竖直圆管内的燃油RP-3换热恶化机理进行了研究，并得到了管径和化热恶化起始点的联系。数值模拟结果证明，质量流量较小时，换热恶化程度与管径成正比，发生节点也会随着管径增加而提前，大质量流量下无传热恶化发生。

除了超临界流体物性参数的变化之外，浮升力作为影响燃油再生冷却效果的一大重要因素近年来也是诸多学者的研究重点。张斌等［19］对竖直圆管内超临界碳氢燃料的传热进行了实验研究，并发现在所有工况中管道入口处均出现传热恶化现象，且不论流向朝上还是朝下，浮升力都会产生显著影响；徐可可［20］在其竖直模拟研究中讨论了浮升力产生对管内航空煤油传热的影响，结果表明，浮升力影响下下管壁传热强化，温度低于上管壁，入口段处浮升力的存在能够削弱传热恶化；孙星等［21］通过数值模拟方法对水平方管中浮升力对超临界航空煤油的换热进行了讨论，同样得出浮升力造成下壁面换热增强的结论，同时指出，浮升力在加热段上游区域对传热影响大，但随着流速增加，其影响渐渐弱化；P.Forooghi等［22］用数值模拟方法研究了倾斜管道内浮升力对湍流对流的影响，结果表明浮升力引起的传热恶化在管道倾斜角度为85°或75°时更为显著，且传热系数的不均匀性在浮升力作用下更明显。

现有文献多数集中于超临界流体内流动和传热的研究，但其中大多研究都针对常重力条件。在实际的运行过程中，飞行器在加速过程中往往会遇到超重力条件。因此，本文建立RNGk-ε湍流三维模型，并对该模型进行对比验证；在常重力、零重力以及多倍重力条件下，对水平管内航空燃油RP-3的流动及传热进行三维数值模拟研究。

1 模型构建

1.1 物理模型

本文采用的计算模型如图1所示，水平圆管内径为1.8 mm，壁厚为0.2 mm，管长为300 mm；管壁上施加均匀热流q；管内流动的工质为航空燃油RP-3；管道内入口流量m=2.5 g/s；管道入口温度为373 K，压力为5 MPa。

图1 水平圆管结构示意图Fig.1 Sketch map of horizontal tube

1.2 控制方程及计算方法

由于管内超临界RP-3航空燃油的物性会随温度发生变化，故选用精度较高的RNGk-ε湍流模型结合增强壁面处理方法。控制方程如式（1）～式（5）所示［17］。

质量守恒方程：

动量守恒方程：

能量守恒方程：

湍动能方程：

湍动能消耗率方程：

式中：μe为有效黏度；λe为有效导热系数；Gk、Gb为湍动能产生项，但二者产生原因不同，前者为层流速度梯度，后者为浮升力；ak、aε为湍流普朗特数；C1ε、C2ε、C3ε为模型常量，取值分别为1.42、1.68和0.084 5［23］。

研究过程中运用有限容积法进行数值计算，采用商业软件Ansys 15.0进行计算，双精度分离求解器进行数值模拟。入口处给定恒定温度和压力条件，湍流程度为0.1。出口边界为压力边界，气液界面为无滑移边界条件，且无渗透。加热段为恒定热流密度边界，加热段长度为300 mm。

本文共选取4套网格进行网格无关性验证，网格数量分别为588 863，881 020，1 452 546，2 014 659，选用q=300 kW/m2，压力为5 MPa，进口质量流量为3 g/s的工况进行数值结果验证。在4套网格的计算结果中，选取水平管内流体平均速度和流体平均温度作为比较参数，如表1所示，可以看出：采用网格数为1 452 546的网格计算得到的流体平均温度和流体平均速度与其他网格计算结果的偏差在1%以内。综合计算速度和计算准确度两个因素考虑之后选取第3套网格进行本文中所有工况的计算。

表1 不同网格的计算结果比较Table 1 Comparison of numerical results from different grids

1.3 模型和数值方法验证

本文计算过程中所使用的物性参数均由文献［24］中的物性参数拟合而来。燃油密度ρ随温度的分段变化为

航空燃油的动力黏度η与温度之间的变化关系为

比热Cp的变化为

导热系数k与温度间的关系为

航空燃油RP-3的临界温度近似为660 K，温度达到该值附近时，物性参数发生剧烈变化。

将本文所计算出的结果与文献［25］进行对比以验证本文数值方法的有效性，比较过程选取上管壁壁温作为比较参数，如图2所示，可以看出：本文计算结果和文献［25］中的数值模拟结果及实验数据均相近，误差不超过5%，因此本文的计算结果是可靠的。

图2 上管壁温度的计算值与实验值的比较Fig.2 Comparison of the wall temperature between the computational and experimental data

2 计算结果分析

2.1 基本工况下温度场及湍动能分析

模拟过程中对常重力条件下，入口流量为2.5 g/s，壁面热流密度q=400 k W/m2时管内航空燃油RP-3的传热与流动进行数值计算。管内流体温度分布云图和沿管长的变化线图如图3所示，水平管内航空燃油在流动方向上逐渐被加热，因此在水平方向上流体温度呈现不断升高的变化趋势，这一趋势在图3（b）中体现得更加明显。靠近水平管中心线的流体温度相近，这一区域的范围大概为z=-0.75 mm到z=0.75 mm。在靠近管壁处流体温度具有较大梯度。

图3 水平管道中流体的温度分布Fig.3 The distribution of bulk temperature in the horizontal pipe

根据管内湍动能强度（对于湍动能E的定义为的分布，如图4所示，可以看出：管内湍流强度沿流向并非线性变化，而是先轻微下降后大幅增加；相比于管中心位置处，越靠近管壁处的湍动能拐点越延后出现。

图4 水平管道中湍动能分布Fig.4 The distribution of turbulent kinetic energy in the horizontal pipe

2.2 浮升力对流动换热的影响

不同重力条件下水平管内燃油RP-3的数值模拟结果表明，水平管内流体的流动在0g～0.5g范围内出现变化［23］。0g和0.5g时同一管截面上的流线图和相应的速度等值线图分别如图5所示，可以看出：零重力条件下管截面上的流动只有朝向管道中心的汇流，此时的速度等值线分布呈现为圆环状的分布；0.5g时管内流动发生变化，在管截面上观察到了明显的二次流动。二次流动产生的原因是管内温度分布的不均匀导致流体出现密度差，从而产生浮升力对流。靠近管壁处流体具有较高的温度，密度较小，管中心附近的流体由于距离加热管壁较远，温度低而密度具有较大值。在重力作用下，管中心的流体朝下流动，较高温度的流体则会向上流动，最终形成双涡卷的流线分布图如图5所示，可以看出：在重力作用下燃油流动速度等值线的分布更加复杂，速度值较大。

图5 不同重力条件下相同管截面上（y=150 mm）流线图及二次流速度等值线分布图Fig.5 The streamline and the secondary flow velocity contours at same tube cross-section（y=150 mm）under different gravitational conditions

不同重力下管壁的温度分布如图6～图7所示。

图6 q=400 k W/m2时下管壁温度沿着管长的分布Fig.6 The distribution of bottom wall temperature along the pipe at q=400 k W/m2

图7 q=500 kW/m2时下管壁温度沿着管长的分布Fig.7 The distribution of bottom wall temperature along the pipe at q=500 k W/m2

从图6～图7可以看出：q=400 k W/m2时管壁温度沿着管长呈现出逐渐增长的变化趋势，且0g和1g情况下管壁温度相近，而4g时管壁温度则低于前两种情况。造成这一差异的原因是浮升力对流的增强，重力加强时，浮升力对流也随之增强，使得更多的管中心处的低温流体流向管壁从而造成管壁温度的下降；当q=500 kW/m2时，在管道后半段中，由于此时流体温度已超过临界温度，物性发生变化，造成重力存在时温度更高的现象。

为了进一步说明重力改变时对流传热强度的变化，定义对流传热系数：

式中：q为表面热流密度；tw为壁面温度；tb为流体温度。

q=400 k W/m2时水平管内的对流传热系数变化如图8所示，可以看出：在水平管道进口处对流传热系数具有较大值，此时边界层较薄，对流强度较大，即“入口段效应”；随着流动的发展，流动边界层厚度逐渐增加而造成了对流传热系数下降，流动持续发展，对流传热系数出现第一个拐点；此后，对流传热系数呈现先增长后减小的变化趋势。0g情况下对流传热系数较小，随着重力的增加，浮升力对流的增强使得传热系数增加，但是传热系数的拐点位置并未发生改变。

图8 q=400 kW/m2时下对流传热系数沿管长的变化Fig.8 The distribution of heat transfer coefficient along the pipe at q=400 k W/m2

当维持其他边界条件不变，增大管壁热流密度至500 k W/m2时，管内对流传热系数随着流动的变化如图9所示，可以看出：对流传热系数在管道前段部分的变化与热流密度为400 kW/m2时的对流传热系数的变化趋势相近，在入口处仍然存在入口段效应，强度稍弱于热流密度为400 kW/m2的工况；对流传热系数的变化在略微增长后出现减小的变化趋势，这一变化相比于热流密度为400 kW/m2时的传热系数变化更加平缓；但在管道末端对流系数的变化则明显和热流密度为400 kW/m2时的变化存在差异，出现了对流传热系数小幅增长的变化；在接近管道出口处，温度已经超过了燃油RP-3的临界温度，因此物性发生突变，从而造成表面对流传热系数在接近管道出口处出现小幅增大的现象。

图9 q=500 kW/m2时下对流传热系数沿管长的变化Fig.9 The distribution of heat transfer coefficient along the pipe along the pipe at q=500 kW/m2

不同重力条件下同一管截面上湍动能E在重力方向上的分布如图10所示。

图10 q=400 kW/m2时湍动能在重力方向上的变化（x=150 mm）Fig.10 The variation of turbulent kinetic energy on the direction of gravity at q=400 kW/m2（x=150 mm）

从图10可以看出：由于湍动能分布是对称的，并且从管中心到管径湍动能的变化趋势是先增加后减小，最终在管壁处减小到0，这一变化与文献［14］中的变化相同。在本文中，为了体现重力的作用，只截取了z=0.25 mm到z=0.9 mm的部分。从图10可以看出：0g和1g下湍动能的变化几乎相同，重力增加至4g时由于浮升力作用的增强，湍动能在z=0.25 mm到z=0.75 mm之间明显增强，但是湍动能开始降低的拐点也相比于其他两种情况下的拐点离水平管壁更远。

q=500 k W/m2时管内湍流强度的变化如图11所示，取y=150 mm处管截面上的湍动能分布，为了直观体现湍动能受重力变化的影响，本文截取了部分变化。