鱼雷燃烧室气-液-固三相耦合换热数值仿真

陈 肯,冯要飞,宗 潇,伊进宝,肖炎彬,史小锋

(中国船舶集团有限公司 第705 研究所,陕西 西安,710077)

0 引言

燃烧室是鱼雷热动力装置中最重要的部件之一[1],推进剂在燃烧室内燃烧,将化学能转变为热能,并推动发动机做功[2]。为保证燃烧室壳体在高温高压的工作环境下具有足够的热强度,需通过冷却水对燃烧室进行冷却换热[3]。在冷却过程中,燃烧室因其结构复杂会出现局部区域过热导致壳体强度失效的问题。为指导燃烧室的结构改进,需要对燃烧室的冷却换热过程进行研究。刘学等[4]利用数值仿真计算的方式,研究了燃烧室多孔介质几何结构对其冷却效果的影响;Chen 等[5]基于双向热流固耦合的方法,采用大涡模拟(large eddy simulation,LES)对非涡流环形燃烧室中的自激旋转模式方位不稳定性进行了研究,并在此基础上构建了频率预测匹配模型;Laurent 等[6]建立了液体火箭发动机的湍流喷雾模型,并根据建立的模型对发动机进行了数值分析;Zhukov 等[7]采用了扩展的涡耗散模型,用CFX 软件对多孔喷嘴头的火箭燃烧室进行了仿真。上述研究聚焦于陆用及航空航天用燃烧室的冷却换热过程,对水下涡轮发动机燃烧室的冷却换热过程不适用。

赵卫兵等[8-9]利用理论解析的方式计算了鱼雷燃烧室在低频不稳定燃烧状态下的温度特性,并针对三组元推进剂的燃烧室流场进行了基于Fluent 的数值研究。Qin 等[10]采用双向热流固耦合的方式,比较了部分进气式涡轮机燃烧室的质量流量和效率的试验与仿真结果,分析了5 kW 功率下涡轮机的输出特性,并建立了涡轮机输出参数与输入边界条件间的一维设计计算模型。上述工作对水下燃烧室的工作状况及冷却特性进行了数值计算研究,但上述研究并未同时涉及气-液-固三相的耦合换热过程。

文中以燃烧室2 种不同冷却水流道结构作为仿真对象,利用Fluent 软件首次开展燃烧室内高氯酸羟胺(HAP)三组元推进剂燃烧,燃气、冷却水流动换热,以及壳体参与导热的气-液-固三相传热换热过程的热流固耦合数值仿真,分析相同工况及冷却水流量下燃烧室不同结构的冷却换热特性差异及造成局部区域过热的原因,为进一步研究和改进燃烧室的结构提供理论参考。

1 冷却换热理论分析

燃烧室的冷却换热过程涉及到的物理模型包括: 燃料的雾化、燃烧模型,燃气及冷却水的流动模型,燃气、冷却水的对流换热模型,固体域的导热模型,以及燃气、固体域的辐射传热模型。

1.1 流体控制方程

采用N-S 方程描述流体域的内流场。

连续性方程

式中: ρ为流体密度;t为时间;u为流体速度;Sm为流体控制体在单位时间质量因源项产生的变化量。

动量方程

式中:P为流体控制体表面的压力;µ为流体的动力粘度;g为重力加速度;F为作用于流体控制体的质量力合力。

能量方程

式中:e为流体控制体的能量;κ为导热系数;Q为控制体内单位时间热源项提供的热量。

式中:x代表一维换热过程中的位置变化;i、j指代节点坐标。

文中,重力对仿真计算的影响较小,故不考虑其影响;冷却水域无源项,故Sm、F、Q均等于0。在燃烧域中,Sm为注入的燃料在单位时间蒸发的流量;F为燃气连续相受到的除重力外的体积力,F=0;Q为液态燃料蒸发吸收的热量以及燃烧产生的化学能之和。

1.2 燃料雾化及燃烧控制方程

1.2.1 燃料雾化

液滴直径取中间质量的定义方式,根据液滴雾化的累计分布R-R 函数,可得液滴直径在R-R 函数分布情况下,按中间质量方式定义的直径表达式为[2]

式中:dm为中间质量直径;d*为液滴破碎后的标称直径;n为破碎系数。

1.2.2 液滴运动及蒸发

液滴经破碎雾化进入燃烧室内的运动控制方程为[11]

式中:up为颗粒速度;ρp为颗粒密度;FX为其他作用力;dp为颗粒直径;Re为颗粒雷诺数;CD为阻力系数。

未达到蒸发温度的温度控制方程为[12]

达到蒸发温度后的液滴质量控制方程为

式中:mp,cp,Tp,εp,Ap分别为颗粒的质量、比热、温度、发射率和表面积;h为对流传热系数;T∞为气相温度;σ为波尔兹曼常数;θR为辐射温度;hfg为潜热。

1.2.3 气体燃烧

组分输运模型可以模拟燃料蒸发产生气体的燃烧过程。组分输运体积反应模型涡耗散(eddydissipation)模型的反应速率Ri,k的控制方程为[6]

式中:mq、mR分别为所有产物组分和某一特定反应物的质量比重;Mw,i、Mw,R分别为第i、R种物质的分子量;分别是反应物和生成物的化学当量系数;A和B是经验常数,分别取4 和0.5。

HAP 三组元推进剂完全燃烧的化学反应方程式为

2 模型前处理及网格划分

2.1 计算域简化及提取

进行燃烧室热流固耦合仿真计算,对燃烧室的几何结构进行适当简化,可以降低网格划分的难度,降低网格数量,提高网格质量,帮助仿真模型更快更好地收敛。

利用提取功能处理简化完成的燃烧室固体域,得到研究对象的全体计算域,如图1 所示。

图1 计算域提取图Fig.1 Extraction graph of computational domain

燃烧室在工作过程中,固体之间的接触部位会产生缝隙,燃气通过缝隙泄漏进入内衬与壳体之间的保温层,扰动保温层内的气体流动,保温层区域由死腔变为流动区域,主要换热形式由自然对流变为强制对流。将固体域变形产生的缝隙修正到提取的计算域模型中,保温层域与燃烧域通过0.25 mm 的缝隙连通。

文中燃烧室采用2 种不同冷却水流道结构,差异主要集中在冷却水域的碗状底部以及最后一段的冷却水出口流道上,如图2 所示。不同结构参数对比见表1。

表1 不同结构参数对比Table 1 Parameters comparison for different structures

图2 不同结构冷却水域流道差异Fig.2 Channel differences of the cooling water area for different structures

2.2 计算域网格划分

采用ICEM 软件完成冷却水域、燃烧域与保温层域的网格划分。对于冷却水域的螺旋流道部分,使用结构化网格生成其中1 条流道的网格,通过周期阵列得到8 条螺旋流道的结构化网格,对冷却水域的碗状底部、出口流道及燃烧域全域,利用八叉树算法生成四面体非结构化网格。

采用Mesh 模块完成内衬与壳体域的网格划分。采用Mechanical 算法生成内衬域及壳体域的四面体非结构化网格。

设置接触面网格尺寸,进行网格的局部加密,保证不同计算域在交界面处的网格单元节点个数与单元尺寸大小大致相等。将划分后的网格导入Fluent 中,利用Make Polyhedral 功能将四面体网格转为多面体网格,显著降低网格单元与节点数量,提高仿真计算的收敛精度与速度。网格划分结果如图3 所示(其中,图(b)为图(a)圈内放大部分)。

图3 计算域网格划分结果Fig.3 Grid division result of computational domain

2.3 网格无关性验证

通过改变网格的大小得到不同的网格划分方案,进行网格无关性验证,最终用于网格无关性验证的4 种网格划分方案见表2,网格无关性验证结果见图4。

表2 网格无关性验证网格划分方案Table 2 Grid division scheme of grid independence verification

图4 网格无关性验证结果Fig.4 Verification results of grid independence

图4 以计算域冷却水域出口处温度分布作为衡量网格数量对计算结果精度影响的指标。从图中可以看出,冷却水域的可接受网格尺寸为0.4 mm,燃气域的可接受网格尺寸为3 mm,在可接受的网格尺寸基础上,继续加密网格的数量不会对仿真计算的结构造成影响,方案2 的网格划分可满足对计算精度的要求。

3 仿真模型选择及验证

3.1 湍流模型选择

流体机械中雷诺数大于105认为处于高雷诺数湍流状态,针对湍流流动状态,仿真模型适合选用两方程模型中的k-ε模型进行流体运动的仿真。k-ε模型包括标准k-ε、RNG(renormailization group)k-ε和可实现k-ε等3 种模型。

以结构1 为仿真的计算域,采用网格划分方案2 的方式划分网格。表3 为进行湍流模型验证时设置的计算域边界条件。

表3 湍流模型选择边界条件设置Table 3 The boundary condition setting of turbulence model selection

监测3 种不同湍流模型的燃气域出口温度和冷却水域出口温度,计算冷却水温升,并分别与同工况下的试验结果对比,求出燃气出口温度与冷却水温升的相对误差,如表4 所示。

表4 不同湍流模型仿真与结果Table 4 Simulation and experiment results of different turbulence models

从表中可知,对于燃烧域的出口温度、冷却水域的出口温升,可实现k-ε模型的计算精度最高,因此湍流模型选择可实现k-ε模型。

3.2 其他仿真模型选择

燃烧室实际工作过程中,燃料液滴在很短的时间内燃烧完全。文中选用组分输运模型作为燃烧模型,该模型的计算原理为气相反应物接触后瞬间完成体积反应,能很好地适用燃烧室的工作状况。

燃烧室的燃烧过程中注入的液滴颗粒会参与到辐射换热过程。由于数据/领域对象(data/domain object,DO)模型适用于所有光学深度的辐射过程及液滴颗粒燃烧过程的辐射模拟,故被选用进行计算。

3.3 仿真模型验证

以燃烧室结构1 为验证仿真模型准确性的计算域,以其在3 种不同工况下的工作环境作为仿真边界条件,进行热流固耦合仿真计算。因冷却水的温度变化数值能直观反映整个系统的冷却换热功率,故选用冷却水温升作为衡量模型准确性的指标。分别对比3 种工况下的仿真结果与试验在结果的误差,以验证模型的准确性。

仿真结果、试验结果以及计算的模型误差大小见表5。

表5 不同工况下仿真与试验结果对比Table 5 Results comparation between the simulation and experiment on the different working conditions

从表中可以看出3 种不同工况下,仿真模型的误差大小均控制在10%以内,建立的仿真模型准确。

4 燃烧室不同冷却水流道结构冷却换热仿真分析

燃烧室2 种冷却水流道结构冷却工况对比时,除结构差异外其余边界条件参数设置见表6。

表6 除结构差异外其余边界条件设置Table 6 The other boundary condition setting except the structural difference

表7 给出了表6 边界条件下针对不同结构的仿真结果。可知,结构2 相较于结构1 冷却水的出口温度更高,燃气温度更低,冷却水温升更高,壳体的内表面均温及最高温度均更低。结构2 的换热功率为101.3 kW,结构1 的换热功率为93.3 kW。结构2 相较于结构1 换热功率提高了8.5%。

表7 不同结构仿真结果对比Table 7 The comparison of simulation results for different structures

探究结构2 冷却效果提升的原因,具体分析2 种结构温度场与流场的差异。

4.1 燃烧域冷却特性

图5 是燃烧室2 种冷却水流道结构的燃烧域温度分布云图,可知2 种结构的燃烧域温度分布几乎一致。2 种结构燃烧域出口段的温度随位置的变化曲线如图6 所示,燃气温度随位置的变化规律大致相同,温度在数值上的最大差距为10 K,相较于燃烧产生的燃气1 400 K 左右的温度,10 K 温差对于传热的影响可以忽略。

图5 不同结构燃烧域温度分布云图Fig.5 Temperature contours of combustion area for different structures

图6 不同结构燃烧域出口端温度分布曲线Fig.6 Temperature curves of combustion area in the outlet for different structures

图7 是燃烧室不同结构燃烧域的速度流线图及在中截面处的速度矢量图,可知2 种结构的流场状态几乎相同,燃气最大流速均约61 m/s,且均在燃烧域下腔与出口流道交汇处达到最大。燃烧域的流场状况相似,可以佐证温度场的分布状况相似。

图7 不同结构燃烧域速度流线图及矢量图Fig.7 Velocity streamlines and vectors of combustion area for different structures

4.2 冷却水域冷却特性

图8 是燃烧室2 种结构冷却水域温度分布云图,可知冷却水在螺旋流道处的温度分布状况基本一致。图9 为燃烧室2 种结构在螺旋流道段的温度分布图,螺旋流道终端温度均约295 K,且温度随位置的变化规律大致相同。

图8 不同结构冷却水域温度分布云图Fig.8 Temperature contours of the cooling water area for different structures

图9 不同结构冷却水域螺旋流道温度分布曲线Fig.9 Cooling water temperature curves in the spiral channel for different structures

由图8 可知,2 种结构的碗状底部及出口流道段,冷却水的温度分布存在较大差别。结构1 碗状底部的冷却水平均温度小于结构2 的温度。结构1 碗状底部的6 条流道中,5 条流道的左侧存在一处大面积的高温区域,右侧则存在一处较为细长的条状低温区域,在碗状流道的最低端交汇后迅速消失。结构2 碗状底部的18 条流道中,同样存在部分流道内形成局部的高温区域。高温区域的大小相较于旧结构更小,温度更高,分布更为随机。但整体上分布于流道上半部分较宽的区域。与出口流道呈90°夹角的2 条流道的整体温度相较于其他流道更高。2 种结构的冷却水域出口流道的温度分布整体上呈现越靠近冷却水域出口处温度越高的趋势,这是冷却水流动换热的原因,符合传热学规律。结构1 相较于结构2 在出口流道段会形成许多局部高温或者低温区域。

图10 是监测的2 种结构冷却水域出口流道一段长度上的温度分布随位置变化规律曲线,可以看出,结构2 冷却水温度越靠近出口温度越高,而结构1 的冷却水温度则在靠近出口处时出现了2 个局部高温。

图10 不同结构冷却水域出口段温度分布曲线Fig.10 Temperature curves of the cooling water in the outlet area for different structures

图11 是燃烧室2 种结构在冷却水域碗状底部差异处的温度云图。表8 为监测所得碗状底部流道不同结构差异处截面温度均值。

表8 不同结构差异处温度值Table 8 The temperatures at the differences for different structures

图11 不同结构差异处冷却水域温度分布云图Fig.11 Temperature contours of the cooling water area at the differences for different structures

从图11 可知,2 种结构的冷却水在流入碗状底部前,在碗状底部入口处的温度基本相等;在流出碗状底部进入出口流道入口时,结构2 的冷却水温度高于结构1。从表8 中可知,结构2 在碗状流道段与出口流道段的换热效率均大于结构1;且在碗状底部流道段的换热功率大于出口流道段,碗状底部的结构差异对于换热效率的影响贡献也大于出口流道段结构差异对于换热效率的影响贡献。

图12 是燃烧室2 种结构冷却水域的速度流线图,图13 是2 种结构冷却水域螺旋流道处水的流速随位置的变化规律。从图12 可看出在螺旋流道处,2 种结构的冷却水流动均近似于层流流动;从图13 可看出,在螺旋流道段2 种结构冷却水的流速分布规律相同且流速均为9 m/s 左右。

图12 不同结构冷却水速度流线图Fig.12 Velocity streamlines of the cooling water area for different structures

图13 不同结构螺旋流道水流速度分布曲线Fig.13 Cooling water velocity curves in the spiral channel for different structures

从图12 可知,结构1 碗状底部5 条流道的左侧区域均形成了较大的漩涡区域,右侧则存在狭长的近似于层流流动的区域。较大的漩涡区域导致冷却水滞止无法及时排出冷却水管道,造成较大面积的局部高温区。狭长的层流区域则通过对流换热的形式将热量带出,形成温度相对较低的区域。

结构1 的5 条流道的层流段在碗状底部最低处交汇,高速水流在底部对冲造成截流现象,速度突变并在第6 条流道以及冷却水域出口流道处形成大范围的湍流区与漩涡区,造成冷却水域出口流道段的冷却水温度分布不均,产生许多局部高温、低温区。

结构2 碗状底部各条流道的上半部分存在许多小的漩涡区域,且漩涡区的大小相较于旧结构更小。小的漩涡区导致许多局部的高温区域产生,其中与出口流道呈90°角的2 条流道的漩涡区域最大,几乎与流道宽度一致,致使通过此2 条流道流出的冷却水量相较于其他流道更少,水温整体较高。

结构2 冷却水域出口12 条流道的流道宽度较为狭窄,水的流动方式近似于层流,水温通过对流传热过程升高,传热距离越长水温越高。

5 结论

文中通过对比分析燃烧室2 种冷却水流道结构的冷却特性,分析其燃烧域及冷却水域的温度场及流场,得到以下结论:

1) 在燃烧域内部结构不变的前提下,壳体内冷却水域的结构变化对于燃烧室燃气的流动及温度影响不大。燃气出口温差仅为10 K,燃气流速大小均为61 m/s,分布基本相同。冷却水域结构的差异对于冷却水温度场及流场影响较大。结构2 的整体冷却效率相较于结构1 提升了8.5%。

2) 随着冷却水域结构底部流道数量的增加,冷却水域底部的湍流及漩涡状态减弱,冷却水域的湍流度降低,水的流速增大,存在的局部高温区域变小。增加底部流道数量,燃烧室换热功率提高,且影响效果明显。

3) 随着出口流道数量以及冷却水由底部流入横向出口流道的通道数量增加,出口处的湍流、漩涡及局部高温区域消失,碗状底部存在的截流现象消失。增加出口流道段的流道数,燃烧室的换热功率提高,但对换热功率的影响效果不如增加碗状底部流道处的流道数的影响效果明显。

文中研究证明了结构2 的冷却换热效果优于结构1,未来可针对结构2 在不同工况下的壳体热应力开展数值计算研究,论证结构2 的冷却水流道结构是否满足燃烧室的壳体强度使用需求,以确定燃烧室的冷却换热结构是否需要进一步优化改进。