基于氮气挤代的推进剂供应系统动态特性研究

王 谦,秦 侃,郝常乐,张安静,罗 凯,党建军

(西北工业大学 航海学院,陕西 西安,710072)

0 引言

随着作战需求的提升与科技水平的进步,世界各国都把增强鱼雷综合性能作为水下武器装备的重要发展方向之一[1-2]。能源供应系统作为鱼雷动力系统的重要组成部分,对鱼雷航速、航程和航深等性能都有较大影响。典型鱼雷的能源供应系统原理[3]如图1 所示。

图1 某型鱼雷能源供应系统原理图Fig.1 Principle of torpedo energy supply system

挤代式推进剂供应系统由高压气舱、电爆活门、管路以及推进剂舱组成。当鱼雷入水后,海水电池一方面激活点火器点燃固体药柱,产生大量燃气后使互锁阀打开;另一方面激活电爆活门,使高压气舱与推进剂舱连通。此时,高压气体快速流出高压气舱并向推进剂舱充气,挤压其中的燃料填充管路空隙并在燃料泵前建立一定压力,以确保燃料泵供出足量的燃料进入燃烧室进行稳定燃烧[4]。若在固体药柱燃尽之前,燃料泵前未建立起指定压力,会使动力系统偏离设计点工况运行,甚至导致启动失败等情况发生[5]。由此可见,挤代式推进剂供应系统[6-7]的动态特性与鱼雷动力系统的启动加速性能密切相关。图1 中的挤代式供应系统与泵供式供应系统采用级联的方式,确保燃料能够快速稳定地向燃烧室供应,但在挤代过程中,管内压力的剧烈变化会对气体产生强扰动,在流场中产生激波从而影响系统的动态稳定性。因此,有必要建立一种方法用于研究挤代式推进剂供应系统的动态特性。

国内针对水下能源供应系统的相关研究较少。罗凯等[8]研制了一种水下热动力能源供应系统,采用比例器对三组元燃料的动态供应进行精确控制,三组元供应比例最大误差不超1.5%;李代金等[9]建立了由燃烧室、雾化喷嘴以及旋转燃烧室组成的水下热动力能源供应系统的数学模型,讨论了柱塞燃料泵以及喷嘴与燃烧室的动态匹配关系,并提出了设计和选择时应当避免的问题。目前尚无针对挤代式推进剂供应系统动态性能的相关公开研究报道。

在管内激波流动研究方面,Roy[10]采用黎曼解算器对管内流体流动进行了理论研究与数值计算,并针对经典Sod 激波管问题,采用Python 代码对其密度、压力、马赫数和熵值进行了近似求解。Chen 等[11]采用高阶点隐式求解方法研究了激波管内激波与边界层的相互作用,得到了二维激波管中激波与边界层的非定常相互作用过程。Qiu等[12]提出了可压缩格子玻尔兹曼方法(lattice Boltzmann method,LBM)仿真激波/边界层在激波管中的相互作用,结果表明加权无波动无自由参数耗散差分格式(weighted non-oscillatory and non-freeparameter dissipation difference,WNND)-Z 方案和WNND-L 方案都能很好地捕捉涡和激波的结构。目前国内外针对激波管内流动的数值计算较为丰富,能够为挤代式推进剂供应系统内部的激波捕捉以及流动仿真提供理论依据。

考虑到国内外对水下挤代式推进剂供应系统动态特性的实验研究有限,且由于系统模型内产生的激波问题,以及管内产生的流动分离和涡流等多动态结构,加大了流体计算软件计算此类问题的复杂性[13-15]。针对此问题,开发了一种一维瞬态可压缩求解器用于解决挤代式推进剂供应系统的动态特性仿真问题。建立一维程序分析推进剂供应系统内部的非定常流动状态,采用计算流体力学(computational fluid dynamics,CFD)仿真和经典文献进行对比,验证了程序的可靠性;采用一维程序研究了不同活门直径、高压气舱压力和推进剂舱容积对系统模型的平衡时间和稳定压力的影响,分析了三者对系统动态特性产生影响的规律,并给出了相关因素对系统动态特性影响的设计参考。通过对系统动态特性的研究显示一维程序可以更加便捷高效地得到所需系统动态特性下的模型结构参数,具有较好的工程应用价值。

1 推进剂供应系统数值模型

推进剂供应系统由高压气舱、管路、活门和推进剂舱等部分组成,其基本结构如图2 所示。

图2 推进剂供应系统结构简图Fig.2 Structure of propellant supply system

1.1 高压气舱和管路数值模型

1.1.1 计算模型

由于鱼雷空间的限制,高压气舱不易过大且需要有一定量的气体储备,因此设计为细长罐结构。同时,考虑到管内流动以轴向流动为主,高压气舱内的气体流动与管路中的流动类似,均可简化为一维流动,其一维连续性方程为

式中:ρ为密度;t为时间;ux为轴线方向的速度分量。

一维动量方程为

式中:τij为垂直于i轴的微元表面上受到的沿j方向的切应力;fx为轴线方向上的单位质量力。

一维能量方程为

式中:E为流体微团的总能,即内能、动能和势能之和;T为温度;λ为热传导系数;h为焓。

采用有限体积法将高压气舱和管路离散为储存流体信息的网格单元,每个网格单元都视为1 个控制体,网格单元之间的交界面视为控制面。图3为典型的储存流体信息的控制体(j),以及相邻控制体的交界面

图3 储存流体信息的控制体Fig.3 Control body for storing fluid information

每一个网格单元内的气体数据都被储存于数据结构之中,包括边界条件、时间步长、网格单元中心位置坐标、网格单元交界面位置坐标、网格单元内流体信息(密度、速度、温度、压力、内能和音速)变化的平均值以及状态通量(质量通量、动量通量和能量通量)的时间导数。为保证程序的2 阶精度,采用网格中心点数据来代表网格单元内气体的各种数据进行计算。

网格单元的交界面上储存着流体经过此壁面的质量通量、动量通量和能量通量,因此流体在高压气舱和管路中的流动就可转化为在多个控制体内流体信息的变化以及控制面上的通量交换。仅考虑理想气体的一维流动,忽略粘性项且假设管路区域中任何状态变化都是连续的,则根据方程(1)～(3)计算出的网格单元内的质量、动量和能量,可用来更新网格单元内的流体信息,即

式中:umass,j、umomentum,j和uenergy,j分别为j单元网格上的质量、动量和能量;vj为单元网格上流体速度;ej为网格单元上流体比内能。挤代系统中的氮气可被认为是热完全气体,则其比内能与温度成正比,表示为

式中,cv为气体的定容比热容。继而根据已求解出的密度与温度通过理想气体状态方程可求出压力

式中,R为气体常数。

声速为

式中:aj为网格单元上的音速;γ为比热比。通过假设每个网格单元内流体信息的线性变化,将网格单元内流体信息的平均值插值到控制体的交界面上,从而得到每个控制体左右两端的流体信息。采用AUSMDV(advection upstream splitting methoddifference vector)方法[16]计算交界面的状态通量。

1.1.2 边界条件

高压气舱与管路的边界条件统一采取在外部扩建2 层虚拟网格的方法,以保证流体的流动在空间分布和时间推进上都具有2 阶精度。设置的2 层网格虽在实际的气体流动物理网格之外,但却包含应用于边界条件的各种数据信息。图4为增加了虚拟网格后,每个储存气体各种数据信息的网格单元结构。

图4 根据实际网格流体信息复制的虚拟网格图Fig.4 Virtual grid copied from actual grid fluid information

若边界条件为壁面时,网格中的参数采用以单元壁面为对称轴的方式复制对称面网格单元中的流体信息,但要注意的是,此时速度的方向应与对称面单元的速度相反;若边界条件为入口时,单元中的质量流量和温度从上一级流入的网格单元中获取;若边界条件为出口时,单元中的温度与压力将传递到下一级的网格单元中。

1.2 活门计算模型

在高压气体向推进剂舱充气挤代推进剂的过程中,活门某一截面处气体流速会达到音速从而发生壅塞现象,因此,通过活门的质量流量和出口温度按照壅塞和非壅塞2 种状态来处理。

将高压气舱出口处的压力和温度赋值为活门入口处的压力和温度,分别记作PL,valve和TL,valve,将管路入口处的压力赋值为活门出口处的压力PR,valve,则活门出口处的温度为

式中:A为活门横截面积;Cd为流量系数,其受活门直径、压力和温度的影响,根据经验Cd的取值范围为0.71～0.93。

1.3 推进剂舱计算模型

考虑到推进剂舱为直径与高近似相等的圆柱状结构,可采用集总参数法忽略推进剂舱的实际体积尺寸,将其视为一点,用该点上流体信息随时间变化的一元函数关系来表征推进剂舱内部整体的变化情况,如图5 所示。

图5 推进剂舱控制体Fig.5 Propellant camber control body

选取推进剂舱为控制体(control volume,CV),管路出口处的接触面为控制面(control surface,CS),显然推进剂舱内恒有质量流入而无质量流出,根据式(1)和式(3)可得

式中:V为体积;e为单位质量流体所具有的能量;Acs为管路出口处的横截面积。

又因

由此,可以获得时间推进当中t+∆t时刻下推进剂舱内部气体流动的质量与能量。根据式(20)～(22)即可得到推进剂舱内部能量、温度和压力随时间的分布变化情况

1.4 时间推进

对于每一个网格单元,其内部各状态参量从时间n推进到时间n+1,采用欧拉预测校正格式迭代计算出每个网格内的质量、动量和能量,再由单元体内已知信息推导出内部各状态参量,即

式中:x为离散控制体单元的空间步长;包含控制体j单元内的质量、动量和能量;包含控制体前后j±1/2单元壁面上的质量通量、动量通量和能量通量;是由上一次迭代计算出的控制体内部的状态参数得到的控制体内储存的质量、动量和能量。

为了保持稳定,时间推进步长限制为

式中: Δtallowed是对所有网格单元设定的最小值;α是收敛条件判断数,在仿真中通常限定α<0.5。由于流动过程为可压缩流动,在Δtsignal的计算上考虑音速a的影响,对于每个网格单元取值为

2 模型验证

2.1 管路流动模型验证

针对经典一维Sod 激波管问题,分别采用一维程序与L1d 程序[17]进行求解,进行交互对比来验证一维程序在流体流动方面的可靠性。

将激波管管径设为D=0.01 m,长度为1 m,其内填充氮气,分为高压和低压2 个部分。壁面边界条件应用于管路两端,交换边界条件应用于管路中间,忽略黏性效应,管路内部初始状态的设定如表1 所示。

表1 计算域模型初始状态Table 1 Initial state of computing domain model

取流动时间t=0.6 ms 时2 种方式下的温度、压力、密度和速度在轴向分布的对比情况,如图6所示。结果显示,一维程序在仿真流体流动方面与Sod 激波管的程序仿真相近,误差较小,通过对比分别采用一维程序、5 阶目标本质无振荡(targeted essentially non-oscillatory,TENO)格式、5 阶加权本质无振荡(weighted essentially non-oscillatory,WENO)格式、5 阶加权紧致非线性方案(weighted compact nonlinear schemes,WCNS)格式、5 阶Steager-Warming 格式和加入人工黏性项,计算经典Sod 激波管问题后得到的数值计算情况,证明了文中数值方法在仿真管内流动方面的精确性,因此所建立的一维程序可以用于管内流动状态的仿真计算。

图6 一维程序与L1D 代码的验证对比Fig.6 Validation comparison of one-dimensional program and L1D code

2.2 系统模型整体验证

文中采用Fluent 软件仿真计算的结果与一维程序计算结果进行交叉验证。

建立如图7 所示的推进剂供应系统模型,模型共分为高压气舱段、活门段、管路段和推进剂舱段4 个部分。给出如表2 所示的典型推进剂供应系统模型各部分尺寸及初始状态设定。

表2 系统模型尺寸及初始状态设定Table 2 System model size and initial state setting

图7 典型推进剂供应系统模型几何尺寸Fig.7 Geometry of typical propellant supply system model

为验证一维程序对模型在某一时刻数值仿真的准确性,如图8 是一维程序求解得到的0.5、1.0和1.5 ms 流动时刻下,模型轴线处的流体参数与CFD 仿真所得结果的对比情况。

图8 不同时间下系统压力轴线变化图Fig.8 System pressure axis changes under different time

从图中可以看出,各个时刻内高压气舱段压力的变化曲线与CFD 仿真结果一致。而推进剂舱中CFD 仿真的管路内出现了压力的震荡,其主要是由于管径变化[18-20]使活门后产生了如图9 所示的膨胀波导致管路内部压力的震荡,而一维程序虽然无法很好地捕捉到其震荡特性,但由图8(d)～(f)可以看出,其对于震荡减弱后的数值捕捉良好,并且其震荡数值均值与一维程序仿真结果的相对误差小于10%。

图9 速度云图Fig.9 Velocity cloud picture

总体来说,一维程序基本可以反映系统模型在轴线方向上流体信息的变化状态。

为进一步验证一维程序对模型在某一截面处流体参数随时间变化的准确性,图10 是采用数值方法求解得到的高压气舱前10 mm 位置和推进剂舱中部处的流体参数与一维程序仿真所得数据结果的对比情况。

图10 高压气舱与推进剂舱流体参数对比Fig.10 Comparison of fluid parameters between high-pressure chamber and propellant chamber

从图10(a)～(c)中可以看出,高压气舱段一维程序仿真得到的结果与CFD 仿真结果基本一致,能够反映高压气舱的泄气过程。而从图10(e)和图10(f)中可以看出,在推进剂舱段的仿真曲线中两者仿真的压力和质量流量效果也接近一致,而图10(d)中一维仿真程序仿真的温度与CFD 仿真结果有些许偏差,其主要原因是CFD 仿真了流体在管路中流动与壁面因黏性摩擦而产生的热量[21],激波较强时温度的真实值较之以热完全气体计算出的相应数值相差较大。因此,CFD 仿真相比于一维程序有更高的温度,两者总体趋势接近,相对误差不超过7%。总体来说,一维程序可以仿真系统模型各个截面流体信息随时间的变化情况。

综上所述,一维程序在仿真系统模型流动方面有较高的准确性,可以节省CFD 数值计算因网格数过大所需的过长时间,因此可以采用一维程序来求解系统模型内流体参数的变化。

3 推进剂供应系统动态特性研究

采用一维程序选取表2 中不同活门直径(D=4、6、8 mm)下的系统模型进行数值计算,计算过程中通过改变高压气舱压力以及推进剂舱容积,获得如图11 所示的系统模型内流动达到稳定时,系统的平衡时间和稳定压力云图。

图11 不同活门直径下的系统动态特性分析Fig.11 Analysis of system dynamic characteristics under different valve diameters

通过不同活门直径条件下,高压气舱压力对系统平衡时间的影响可以看出: 由图11(a)所示,在D=4 mm 时,平衡时间随着高压气舱压力的升高而缓慢增大,高压气舱的压力每升高2 MPa,最终平衡的时间增加0.006～0.01 s,而当活门直径增大到D=8 mm 时如图11(c)所示,平衡时间基本不会再随着高压气舱压力的升高而变化。这可能是由于在活门直径较小时,系统内部不平衡的压差导致通过活门产生的入射激波与反射激波相互干扰而导致系统达到平衡的时间延长,而增大活门直径后,活门处产生的入射激波大幅减小[22],对平衡时间的干扰也逐渐减弱。

类似的,不同活门直径条件下,推进剂舱容积变化对系统平衡时间也有相应的影响。从图11(a)～(c)中可以看出,在活门直径为D=4 mm 时,推进剂舱的容积每增加2 L,最终平衡时间增加0.06 s,随后随着活门直径的增大,推进剂舱容积的改变量对最终平衡时间的影响逐渐减弱,在活门直径变为原来的2 倍后,推进剂舱容积的改变量对最终平衡时间的影响减弱为原来的1/5。由式(12)～(14)可以看出,通过活门的质量流量也会随着活门直径的增大而增加,从而导致平衡时间可以在一定的容积范围内更快地达到平衡。

由图11(d)～(f)中不同的活门直径条件下,高压气舱压力及推进剂舱容积变化对系统稳定压力的影响可以看出: 系统稳定压力随着高压气舱压力的升高而增加,但增加的幅度又随着推进剂舱容积的增大而减小,这也符合理想气体状态方程中的描述。随着推进剂舱容积的增大,最终平衡的压力从一开始的高压气舱压力每增加2 MPa,平衡压力增加0.12～0.16 MPa,减小到高压气舱压力每增加2 MPa,平衡压力增加0.1 MPa;随着高压气舱压力的升高,最终平衡的压力从一开始的推进剂舱容积每增加2 L,平衡压力减小0.04 MPa,增大到推进剂舱容积每增加2 L,平衡压力减小0.08MPa。

不同活门直径、高压气舱压力以及推进剂舱容积条件下系统的平衡时间与稳定压力各有不同,说明不同因素的变化情况对不同条件下系统的动态特性均有影响。下文通过一维程序计算得到的不同位置处的流体信息,进一步分析三者对系统动态特性的影响。

3.1 高压气舱压力对系统动态特性的影响

改变表2 中高压气舱压力,分别取P=10、15、20 MPa 后进行一维程序仿真,得到高压气舱末端与推进剂舱处的流体压力特性曲线如图12 所示。

图12 不同高压气舱压力下系统动态特性对比Fig.12 Comparison of system dynamic characteristics under different pressures of high-pressure chamber

由图12(a)可以看出,高压气舱压力分别为P=10、15、20 MPa 时,系统平衡时间分别为0.85、0.86、0.86 s,可见高压气舱内压力的变化对高压气舱与推进剂舱达到平衡时间的影响较弱。根据小孔出流规律,更高的压差会增大活门出口处的质量流量,加速系统内整体压力的平衡,因此高压气舱压力的变化对系统平衡时间不敏感。由图12(b)可以看出,随着高压气舱内压力的上升,由于系统为密闭容器没有与外界联通,因此推进剂舱内最终稳定的压力也随之上升,在高压气舱压力分别为10、15、20 MPa 的条件下,系统最终稳定压力分别为0.64、0.95、1.28 MPa。

3.2 推进剂舱容积对系统动态特性的影响

改变表2 中推进剂舱容积,分别取V=22、24、26 L 进行一维程序仿真,得到高压气舱末端与推进剂舱处的流体压力特性曲线如图13 所示。

图13 不同推进剂舱容积下系统动态特性对比Fig.13 Comparison of system dynamic characteristics under different propellant volumes

由图13(a)可以看出,在推进剂舱容积分别为V=22、24、26 L 的条件下,系统最终平衡时间分别为0.79、0.82、0.86 s,随着推进剂舱容积的增大,系统最终平衡的压力也会随之降低,而通过活门的质量流量依旧不变,高压气舱需要更长的时间达到更低的压力,因此高压气舱与推进剂舱达到稳定的时间缓慢增大。由图13(b) 可以看出,随着推进剂舱容积的增大,根据理想气体状态方程可知,最终整个系统达到稳定时刻的稳定压力会随着推进剂舱容积的增大而减小。在推进剂舱容积分别为V=22、24、26 L 的条件下,系统最终稳定的压力分别为1.43、1.28、1.20 MPa。

3.3 活门直径对系统动态特性的影响

改变表2 中活门直径参数,分别取D=4、6、8 mm 后进行一维程序仿真得到高压气舱末端与推进剂舱处的流体压力特性曲线如图14 所示。

图14 不同活门直径下系统动态特性对比Fig.14 Comparison of system dynamic characteristics under different valve diameters

由图14(a)可以看出,在D=4、6、8 mm 的条件下计算得到的高压气舱最终平衡时间分别为2.01、0.82、0.45 s,可见随着活门直径的扩大,相同条件下活门出口流出的质量流量增大,系统达到稳定的时间前移。图14(b)中推进剂舱内稳定的压力会有不同,这是由于一维程序在流体仿真上无法反映流体在通过突扩管和突缩管内时流体的真实情况,如程序忽略了将活门看做整体进行集总参数的计算时会产生不同的压力损失和热损失导致最终的计算偏离理想气体的条件,以及流体在管内流动的局部损失等,从而会导致最终稳定的压力有些许偏差。

4 结束语

文中根据典型鱼雷能供系统建立了描述推进剂系统特性的一维数值仿真程序,通过与经典文献以及CFD 仿真的对比,证实了一维程序的可行性。利用一维程序对系统模型不同管径、高压气舱压力和推进剂舱容积的条件下展开仿真,得出以下结论:

1) 推进剂供应系统挤代过程中,活门管径的适当增大可以明显缩短系统达到平衡的时间。在活门管径由4 mm 变为8 mm 后,活门处质量流量增加,最终使得系统的平衡时间提前至原来的22.4%。

2) 高压气舱的压力对系统达到平衡时间的影响随着活门管径的增大而减弱。活门直径D=4 mm 时,高压气舱压力每升高2 MPa,平衡时间增加0.06 s,稳定压力增加0.128 MPa;在活门直径变为原来的2 倍后,高压气舱压力变化对平衡时间基本没有影响,其主要与不同管径下活门处产生的膨胀波和激波有关。

3) 推进剂舱容积的增大会使封闭系统内最终平衡的压力下降,导致高压气舱需要更长的泄气时间以达到更低的压力。在此条件下,将活门直径由4 mm 变为原来的2 倍,则会在出口处产生更大的质量流量,使得推进剂舱变化对系统的影响减弱为原来的1/5。

文中对于一维程序仿真挤代式供应系统精度的讨论还较为粗略,后续将结合仿真计算与试验分析,对挤代式系统开展深入的研究,以期获得更加适用于工程实践的挤代式供应系统仿真程序。