基于准一维模型的低温流体预冷充填管路仿真分析

2023-05-05任孝文周晨初陈宏玉刘亚洲

火箭推进 2023年2期

任孝文,周晨初,陈宏玉,李平,刘亚洲

(1.西安航天动力研究所液体火箭发动机技术国家重点实验室,陕西西安 710100;2.航天推进技术研究院,陕西西安 710100)

0 引言

在液氧/煤油发动机、液氧/液氢发动机以及液氧/甲烷发动机[1-4]等新一代运载火箭主动力中,均存在低温推进剂预冷充填的复杂相变过程。例如,在发射前的地面加注阶段需要对主/副阀前管路进行预冷,在起动阶段则需要对阀后管路及喷注器进行预冷充填。在这些过程中,低温推进剂吸收管壁的热量产生蒸汽,当蒸汽排出的速率小于蒸汽产生的速率时,蒸汽受到压缩会在管内产生瞬时的压力升高,同时液体与蒸汽的相互作用会在管路内形成一定的压力波动;随着管壁温度的逐渐降低,管内蒸发量减小,流动由气液两相流动逐渐转化为纯液相流动。

对于使用低温推进剂的上面级高空发动机而言,贮箱隔离阀至主阀前管路无法如一级/助推级发动机在地面进行预冷,需要在高空开展这一过程,较为准确地预测预冷时间以及推进剂消耗量,能够显著提高推进剂的利用率,最大限度提升有效载荷。同时掌握低温推进剂对阀后管路及容腔的预冷充填瞬态特性[5],能够为燃烧组件点火时序的制定提供重要参考。

目前国内外现有的管路预冷模型大多忽略充填的动态过程而将关注点集中于预冷消耗量和预冷时间方面。典型的研究工作有:李永兵等进行的输送系统循环预冷数值仿真,其通过一维两流体模型详细研究了各种参数对预冷效果的影响[6];王娇娇等采用准稳态方法建立了一维预冷计算模型,其采用一维均相法,考虑了辐射漏热的影响,结果表明存在一个最优质量流速使得预冷时间缩短的同时减少推进剂的浪费[7];Darr等假设流动速度恒定,进行了液氮预冷管路的一维计算[8]。

在考虑充填等瞬态流动的预冷过程方面,程谋森等针对氢/氧液体推进剂供应管路的充填过程进行了研究,其采用一维均相模型以统一的方程描述亚临界和超临界态的流体流动,使用特征线差分方法进行数值求解,计算结果近似反映了充填过程的动态特性[9];同样针对氢氧发动机供应管路,刘昆等使用有限元分割的方法采用一维均相模型进行了仿真计算,其中两相单元使用了等效流容方程[10]。此外,美国GFSSP[11]、欧洲ESPSS[12]、德国THESEUS[13]等程序在处理两相流动时均使用了一维均相模型,该模型在每个控制体内假设所有相态的压力、温度及速度是均匀平衡的,忽略气液相间的滑移速度。

在两流体模型方面,由于其极大的复杂度,至今未有成熟的、完整的两流体模型应用于充填过程的系统级仿真中。近年来,针对欧洲新一代低成本可重复使用液氧甲烷发动机的动态仿真要求[14-15],冯卡门研究中心(Von Karman Institute)在与罗马大学和西班牙EAI公司(Empresarios Agrupados Internacional)的紧密合作下基于ESPSS仿真模型库在两相流动领域进行了持续大量的投入,其研究表明得到预期稳定的非均匀两相流模型是存在一定难度的[16],最新发布的ESPSS 3.2.4测试版本加入了两流体六方程模型来处理两相流动,但目前通过的测试仅限于不含相变的理想空气与水的两相流动问题[17],能够处理相变情况下的两相流动仍在持续的研发完善中。

本文在常温推进剂充填管路的水击和流动振荡研究的基础上[18-22],采用非因果建模的Modelica语言[23]基于MWorks仿真平台建立了准一维有限体积的预冷充填模型。本文模型与公开文献已有模型的不同之处在于,除基本控制方程外新增了充填率对流方程,并采用一维形式求解该方程及能量方程,质量方程及动量方程则使用集总参数形式。通过定义相同的接口模块,本文建立的预冷充填管路组件模型可与液体火箭发动机瞬态特性模块化通用仿真模型库[24]中的其他组件模型进行数据交换,拓展了该模型库的仿真能力。最后,为获得低温流体预冷充填管路的主要特性,利用本文模型分别研究了不同结构初温及不同低温介质对预冷充填过程的影响。

1 计算模型

1.1 数学模型

低温流体预冷充填管路的典型流动物理过程如图1所示,根据流型图谱可细分为泡状流、反环状流、弥散流等[25]。

图1 预冷充填典型物理过程Fig.1 Physical process of pre-cooling filling

为简化模型并提升系统仿真效率,假设蒸汽与低温流体之间具有明确的分界面且垂直于流动方向,分界面的移动速度即为液体充填的速度,且气液之间不存在速度滑移,即气液等速;忽略流体之间的导热,管路横截面沿流向为定值;流动过程使用集总参数的方法描述,换热过程和充填率的计算采用一维分布参数的形式描述,如此得到的计算模型如图2所示。

图2 预冷充填计算模型(i为网格索引)Fig.2 Simulation model of pre-cooling filling (i is index of node)

当液体未充满管路,即管路内总的充填率αt< 1时,根据质量守恒定律,液体部分的压力计算式为

(1)

根据质量守恒定律,气体的积存量变化由液体汽化率和气体排出的流量决定,即

(2)

结合气体状态方程,气体部分的压力计算式为

(3)

式中:pout为管路出口压力;L为管路总长度;R为蒸汽气体常数;T为流体温度。

液体部分流出流量的计算式为

(4)

式中:d为管路直径;λ为达西摩擦系数;ρL为液体密度。

由于气液等速的假设,气液界面的移动速度usf计算式为

(5)

气液界面的位置,即已充满液体的长度l计算式为

(6)

管路总的充填率αt为

(7)

当液体充满管路,即管路内总的充填率αt=1时,上述方程均切换为如下液体的质量与动量方程。

液体的压力为

(8)

(9)

式中pin和pout分别为管路充满液体后入口和出口的压力。

液体的流量为

(10)

一维分布参数形式的能量方程为

(11)

式中:q为壁面传给流体的热流密度;表观比热容cp由饱和状态的气相比热容cpg和液相比热容cpL计算得到,即

cp=αcpL+(1-α)cpg

(12)

式中α为单个网格内的液体充填率,表观密度ρ具有相同的形式。

一维形式的充填率对流方程为

(13)

这样在求解过程中,随着充填过程的进行,每个网格中的充填率α将从0到1变化,每个网格内能量方程中的密度与比热容即会从气相逐渐转变为液相。

(14)

低温流体与过热壁面的沸腾换热过程可根据壁面温度Tw与再浸润温度TMHF、临界热流温度TCHF及流体饱和温度Tsat的大小划分为膜态沸腾、过渡沸腾以及核态沸腾3个区域:

1)Tw>TMHF,膜态沸腾阶段;

2)TCHF

3)Tsat

具体计算方法可参见文献[7]。

不考虑管壁之间的导热以及环境的热泄漏,每个网格内的管壁释热方程为

(15)

式中:Tw为管壁温度;mw为管壁质量;cw为管壁比热容。

1.2 计算方法

本文的数值计算基于Modelica语言的面向对象建模的MWorks仿真平台,由于该平台能够使用时间步长自适应的DASSL算法[26]求解代数方程和常微分方程组,因此只需要将能量方程及充填率对流方程的空间偏导部分进行数值离散,将分布参数的偏微分方程组转化为网格单元内的常微分方程组。

由于充填锋面垂直于管路轴向的假设,网格内充填率只有0和1两种状态,因此随着充填状态由未充填到充满的变化,单个网格中的充填率是一个从0至1的阶跃过程。这就要求空间偏导的离散具有迎风性好、耗散性小的特点,对此,本文采用TVD格式[27]进行空间离散。

本文开发的模型将一阶、二阶迎风格式及7种TVD格式全部集成在了该组件模型中,并可以在用户界面下拉菜单进行选择,如图3所示,大大降低了模型应用难度。

图3 参数设置用户界面Fig.3 User interface of parameter setting

2 模型验证

2.1 液氮预冷充填管路

美国弗洛里达大学Darr等在NASA马歇尔空间飞行中心的资助下,开展了液氮预冷管路的试验研究[8]。试验的目的在于获得低温液氮充填垂直不锈钢管路的预冷曲线,为了能够对预冷沸腾中3种沸腾机制提供足够多的数据支撑,试验工况涵盖了较为宽广的质量流量和压力范围。试验中选择液氮作为低温介质,不仅是出于安全性的考虑,还因为液氮的表面张力、沸点和汽化潜热等物性与液氧较为相似。试验中,流动方向与重力方向相同,对于所有的工况,试验段管路的初始温度保持在室温293 K。

图4所示为针对Darr试验段搭建的仿真模型,其中预冷充填管路分为20段。

图4 针对Darr试验的仿真模型Fig.4 Simulation model for Darr experiment

图5所示为6种不同工况下管路中部壁面温度随时间的变化曲线,其中G为质量通量,ReL为流动雷诺数,pin为入口压力。从对比结果来看,本文仿真模型能够较好地模拟各个工况下的预冷过程中壁面温度的变化,而Darr等学者的模型与试验数据有一定的差距,原因在于其模型在计算过程中假设流动过程的流量及压力不随时间变化,而仅仅考虑了温度的动态变化过程;而本文模型在考虑低温流体对管路充填这一动态过程的同时计算了流体与管路之间的热交换。

图5 仿真结果与试验数据对比Fig.5 Simulation results versus experimental data

对于管路预冷过程而言，在工程实践中通常比较关心的是预冷时间，即管壁温度需要多长时间才会降至低温流体的温度。因此，应以预冷时间为指标进行仿真模型的误差对比，误差计算公式为

(16)

式中:tsim为预冷时间仿真值;texp为预冷时间实验值。

表1给出了本文仿真模型与Darr等仿真模型对预冷时间的计算误差。

表1 预冷时间仿真误差

从表1可以看到相比于实验数据,两种模型的计算误差均不超过13%,而本文模型误差相对更小(小于6%),这表明本文模型对预冷流动过程的描述更为准确,且使用的对流换热经验关联式也较为合理。

2.2 液氧预冷充填管路

为进一步确认本文建立的管路预冷充填模型的有效性,对Jin等学者的试验[28]进行了仿真。试验中流体介质为低温液氧,试验系统如图6所示,管路全长7 m,内径10.2 mm,管路外包覆有厚度21 mm的绝热层,以减少周围环境向管路的热泄漏,在管路外壁面如图A、B、C位置布置有温度传感器。仿真模型见图4,管路分段数为100。

图6 Jin等试验系统示意图Fig.6 Diagram of experimental system of Jin et al

图7及图8给出了两种工况下的实验数据与仿真结果对比,本文仿真模型计算得到的壁面温度变化基本与实验趋势一致。由于管路较长,且质量通量G较小,因此管路预冷时间较长,在实验工况范围内最长预冷时间达到了400多秒。且从不同管路位置的壁面温度曲线来看,管路上游(A位置)的冷却速度快于管路下游(B位置)的冷却速度,管路中部(C位置)的冷却速度介于两者之间。

图7 仿真结果与试验数据对比[G=24.5 kg/(m2·s), Re=2 071, pin=557 kPa]Fig.7 Simulation results versus experimental data[G=24.5 kg/(m2·s), Re=2 071, pin=557 kPa]

图8 仿真结果与试验数据对比[G=41.6 kg/(m2·s), Re=4 010, pin=817 kPa]Fig.8 Simulation results versus experimental data[G=41.6 kg/(m2·s), Re=4 010, pin=817 kPa]

图9与图10所示为对应于两种工况下,本文模型计算得到的管路内压力及充填率变化曲线。由于预冷过程中,液氧吸收常温壁面的热量发生了汽化,因此在管路充满前,管内处于气液两相状态,管内压力主要由汽化而来的氧气决定,且从图中看到随着充填的进行,管内压力出现了一定的波动。而在液氧充满的瞬间,由于图6中阀门3流阻的作用,产生了明显的水击压力。

观察图7及图9,在pin=557 kPa的工况下,图9中充填率曲线表明管路内液氧在300 s左右即显示充满,而在图7中,300 s时刻管路中部(C位置)仍未完全冷却至流体温度;在pin=817 kPa的工况下,图10中充填率曲线表明管路内液氧在220 s左右显示充满,在图8中,此时管路末端壁面温度未冷却至流体温度。这一不合理之处是由于本文模型的假设造成,其假设充填过程中流动具有明确的气液界面,且垂直于管路内流体流动方向,该气液界面左侧为纯液相,右侧为纯气相,如图2所示。

图9 管路压力及充填率[G=24.5 kg/(m2·s), Re=2 071, pin=557 kPa]Fig.9 Pressure and filling rate of pipeline [G=24.5 kg/(m2·s), Re=2 071, pin=557 kPa]

图10 管路压力及充填率[G=41.6 kg/(m2·s), Re=4 010, pin=817 kPa]Fig.10 Pressure and filling rate in pipeline[G=41.6 kg/(m2·s), Re=4 010, pin=817 kPa]

该假设导致准一维管路网格内的充填率非0即1,并没有介于0～1之间的状态,而实际上在壁面温度还未完全冷却至流体温度的情况下,流体中依然含有较少的蒸汽。即本文建立的准一维管路预冷充填模型中,充填率实际并不是严格的纯液相的充填率,而是液体含量较大的两相流体与纯气相之间的分界面。这也是为什么在液氧质量通量较大时,图8中管路末端温度降低至液体温度的时刻与图10中充填率达到1的时刻较为接近,分别为270 s与220 s左右;而在液氧质量通量较小时,图7与图9中的相应时刻分别为400 s与300 s左右。

这一点是该模型的局限性,但是与壁面温度的实验数据对比可知,该模型在预测管路预冷时间方面是较为准确的。同时由于该模型计算得到充填率为1时,从图7和图8中也可看到,膜态沸腾阶段已经基本结束,大流量的低温流体已经充填至管路末端,而燃烧组件点火时序的制定依据也是大流量推进剂的进入时刻,因此后续将该模型应用于发动机起动过程的动态仿真分析是合理的,是有一定的工程意义的。

3 仿真分析

3.1 壁面初始温度对预冷充填管路的影响

使用图4搭建的管路预冷系统模型,研究壁面初始温度对预冷充填管路的影响。低温介质为液氧,管长2 m,内径10 mm,壁厚1.25 mm,贮箱供应压力1.0 MPa,该参数下对应的稳态流动质量通量G=12 755 kg/(m2·s),0 s时刻阀门开启,充填开始进行。在上述参数一定的前提下,分别计算管路初始温度300 K、400 K、500 K、600 K时的预冷充填过程。

图11所示为计算得到的管路不同位置的壁面温度变化曲线,图12所示为不同初始壁温下管路内压力变化曲线的对比,图13所示为不同初始壁温下管路内充填率的对比。从以上计算结果可以看到,随着管路初始温度的增高,液氧对管路的预冷时间逐渐增长:管路末端位置预冷完全的时间从管壁初始温度300 K工况下的3 s逐渐增加至管壁初始温度600 K情况下的6 s。

当管路初始壁温升高时,管路内充填开始时由蒸发产生的氧气量也逐渐增多,导致图13中的充填时间逐渐增长,同时导致了在图12中蒸发气体的建压时间缩短,气体建压的压力峰值逐渐增大,管路内压力波动也更加明显;并且由于气体的缓冲作用,在低含气率的液氧充满管路时管路末端节流元件位置的水击压力逐渐减小,各项具体数据见表 2。

图12 不同初始壁温下管路出口压力变化Fig.12 Change of pipeline outlet pressure under different initial wall temperatures

图13 不同初始壁温下管路内充填率Fig.13 Filling rate under different initial wall temperatures

表2 不同初始壁温下管路充填过程的特征数据

3.2 不同低温介质对预冷充填过程的影响

在3.1节液氧预冷充填管路的基础上,更换流体介质为液甲烷和液氮进行不同低温介质预冷充填过程的对比研究,仿真过程参数设置与3.1节相同,其中管壁初始温度保持在300 K,且入口流体温度均设置为100 K。

图14给出了3种低温流体预冷充填过程中壁面温度的变化曲线,可以看到,在给定管路入口压力温度及管路流阻的情况下,液甲烷的预冷速度最快,液氮预冷速度最慢,液氧介于两者之间。但是由于3种低温介质的密度不同(3种介质中液甲烷密度最小,液氧密度最大),因此在入口压力及管路流阻相同的条件下,3种流体的稳态流速是不相同的,液氧、液甲烷、液氧流速分别为11.87、18.83、13.69 m/s;而由于流体的流速对对流传热过程的影响较大,因此若要分析介质的不同(主要是物性参数的不同)对预冷充填过程的影响,则需要将流速的影响剔除。

图14 预冷充填过程中管路出口位置壁面温度变化Fig.14 Change of wall temperature at outlet of pipe

基于上述分析,针对不同介质的预冷充填过程,保持管路结构参数一定,调整管路流阻。由于充填过程为非稳态过程,充填过程的流速是随时间变化的,而在非稳态过程中保持3种介质的流速相同是很困难的,因此此处采用调整管路流阻以保证3种低温介质的稳态流速相同的办法来尽可能剔除流速的影响,稳态流速均为11.87 m/s。

图15所示为相同稳态流速下的液氧、液甲烷及液氮预冷充填管路的壁面温度计算结果,该温度为管路出口位置的壁面温度。

图15 相同流速下管路出口位置的壁面温度Fig.15 Change of wall temperatures under same flow velocity

可以看到液氧预冷时间最短,用时3 s,液甲烷的预冷时间次之,用时4 s,其中液氧的预冷速度先慢后快,液甲烷的预冷速度先快后慢,而液氮的预冷速度最慢,所需要的预冷时间最长,为5.5 s。

为分析出现以上差别的原因,图16～图18分别给出了3种低温流体充填过程的壁面温度变化以及再浸润温度(temperature of minimum heat flux)与临界热流密度温度(temperature of critical heat flux)对沸腾换热区域的划分。

图16 液氧预冷充填管路的沸腾传热Fig.16 Boiling heat transfer of liquid oxygen pre-cooling filling pipe

图17 液甲烷预冷充填管路的沸腾传热Fig.17 Boiling heat transfer of liquid methane pre-cooling filling pipe

图18 液氮预冷充填管路的壁面温度及沸腾传热分区Fig.18 Boiling heat transfer of liquid nitrogen pre-cooling filling pipe

可以看到,对于图16中的液氧而言,预冷过程中膜态沸腾、过渡沸腾以及核态沸腾3种沸腾机制均较为明显;而对于图17中的液甲烷,预冷过程中过渡沸腾与核态沸腾占据主导,膜态沸腾对预冷过程的影响较小。由于核态沸腾阶段的对流换热系数较大,而膜态沸腾阶段的对流换热系数较小,因此在大部分时间内(液氧进入核态沸腾前),液甲烷的预冷速度快于液氧;而由于液氧与液甲烷物性的不同,在核态沸腾阶段,相比于液甲烷,液氧与壁面的对流换热系数更高,因此液氧预冷所用时间更短。

图18中的液氮预冷过程中,在大部分时间内膜态沸腾占据主导,因此液氮的预冷速度最慢。造成上述3种低温介质预冷换热分区不同的原因主要在于,相同压力下,液氮的饱和温度最低,液甲烷的饱和温度最高,液氧的饱和温度介于二者之间,因此充填过程中,图16～图18中液氮的膜态沸腾区域占比最大,液甲烷的膜态沸腾区域占比最小,液氧的膜态沸腾区域介于两者之间。