控制方程中较为关键的一项即源项,源项中包含了质量守恒方程的源项、动量守恒方程的源项和能量守恒方程的源项。其中质量守恒方程的源项代表液态推进剂蒸发成气态,从而给气相控制方程带来的加质项。其数学模型可以用蒸发和掺混的子模型来描述。动量守恒方程的源项代表了气体与壁面的摩擦损失。能量守恒方程的源项代表了气相推进剂之间化学反应的放热项以及气体与壁面热交换带来的加热或放热项。如何通过各个子模型准确地描述上述物理过程是燃烧室数学建模的关键。

加质项可表示为^[13-14]

(6)

其物理意义为假设液态推进剂的蒸发沿燃烧室的分布形态为正弦信号的前半个周期,并且蒸发的范围限定于区间[l_s,l_f]内,式(6)满足,即保证了质量添加的守恒性^[15]

(7)

释热项主要是求反应焓H_s的问题,通过热力计算可以得到推进剂化学反应的平衡温度,即

T_eq=F_m(K_c,p_c)(8)

式中:K_c 为混合比; p_c 为燃烧室压力。因此反应焓可以表示为

Hs=∫TeqT0CpdT(9)

值得注意的是,混合比K_c 以入口的氧化剂和燃料的质量流量的比值来表示[见式(10)],而不是以燃烧室内总的氧化剂和燃料的质量积存量来表示,因此燃烧室的热力参数可以敏感地对入口流量的波动进行响应。这样就解决了一大关键问题,即流量的高频扰动可以实时地被模型所反映,进而反馈到燃烧室的其他热力参数的计算上。故而该模型能够本质地适应高频相应的计算。

K_c=(q_m,o)/(q_m,f)(10)

空间离散采用了ROE格式,具体格式见文献[16]。由于存在源项和动量方程中面积变化带来的附加项,直接采用ROE格式计算会导致在面积变化处产生流量不守恒的现象^[17],需对式(1)的特征向量进行修正,具体形式为

(13)

式中符号的含义见文献[16-17]。

对于火箭发动机而言,燃烧室内蒸发和化学反应的时间尺度非常小,而流动和传热的时间尺度较大,即意味着其数学模型存在很大的刚性,显式的欧拉法或龙格-库塔法求解刚性微分方程较为困难^[18],因此时间离散采用了隐式变步长的时间推进Dassl方法^[19]。

以线性分布热源的问题为例,模型示意图如图1所示。气体工质为空气。流通通道为等截面,忽略气体的摩擦力,对于理想无黏流动,初始时刻加热量为零,则各个截面上气流参数相同,入口静压p=1×10⁵ Pa,静温T =300 K,单位横截面积上的气流流量q_m=0.1 kg/s,管道总长度L=0.1 m,管道截面积A=0.01 m²,气流通道内的热源分布为

(14)

图1 线性热源问题示意图
Fig.1 Schematic of the flow problem with linear heat source

取C_p=1 004.5 J/(kg·K),ΔT=1 900 K。计算结果如图2所示。结果表明,在平直段没有加热,参数不变; 在加热段,当亚音速气流被加热时,气流温度和流速增加,密度和压力减小。对于密度、速度和温度的计算,数值仿真结果和精确解吻合较好。对压力的计算,数值仿真结果在加热段的起点和终点处与精确解存在一定的误差,100个分段单元时,最大误差仅为0.02,增加分段单元数时,误差会逐渐减小。

图2 线性热源问题计算结果
Fig.2 Calculation result of the flow problem with linear heat source

以文献[20]中的试验结果为例,对含熵波和声波的气路算例的求解进行验证。如图3所示,该试验装置是由一段长约1.5 m和直径0.12 m的圆柱形流动管路组成,流路出口设置短喷管使流动达到临界音速,冷空气和热空气经过临界气喷嘴送入该试验装置流路。冷空气流路喷嘴前设置谐波干扰,使冷空气产生流量振荡,进而引起掺混气体的温度变化,以形成熵波,增大冷空气的流量会导致混合气温度降低,因此将冷空气等效为燃料,热空气等效为氧化剂。该方案可等效于富燃燃烧室,可以模拟富燃燃气发生器燃料流路中的干扰。高低温气体在头部处通过很多气喷嘴直接进行掺混,混合距离可以忽略不计,认为在燃烧室入口处已经完成了混合。

图3 试验装置示意图
Fig.3 Schematic of the testing apparatus

给定入口边界条件,热空气温度T_h=780 K,流量q_m,h0=6.84 kg/s,冷空气温度T_c=300 K,流量q_m,c为周期性变化,表达式为

q_m,c=q_m,c0+Δq_m,csin(2πft)(15)

式中q_m,c0=10 kg/s。

分别采用集中参数模型和本文的分布参数模型对该试验装置的动态特性进行仿真计算,并将时域解变换至频域,结果如图4所示,结果表明,在ω=0～600 s^-1的低频范围内,气路的幅值响应有逐渐降低的趋势,即对低频幅值响应更敏感,相位差滞后随频率的增加有增加的趋势,在100 s^-1以上的频率范围,在幅值和相位特性上均产生了熵波响应。

图4 试验装置末端的频率特性
Fig.4 Frequency characteristics of the testing apparatus

从幅值上看,本文的分布参数模型和集中参数模型均和试验结果吻合较好,区别在于,集中参数模型不能捕捉到100 s^-1以上的熵波现象,而分布参数模型完全能够反映出熵波现象。再从相位上看,集中参数模型在大于100 s^-1的范围内,完全和试验结果相背离,而分布参数模型与试验结果能够吻合,再次验证了集中参数模型只能适用于低频范围内,ω=100 s^-1 对应的频率为100/2π,即不到20 Hz。这与文献中提到的频率范围相一致。而本文提出的分布参数模型所适用的频率范围远大于集中参数模型。

进一步地,将研究的频率范围扩充至中高频范围,取无量纲角频率0<ωL/a<8,计算气路的频率特性,并与试验数据比较,图5为气路试验系统管路末端的压力频率响应特性,图6为气路试验系统管路中点位置的压力频率响应特性。

图5 气路试验系统管路末端压力频率响应特性
Fig.5 Pressure frequency characteristics of the testing apparatus at the end point of the gas tube

结果表明,在更高的频率范围内,采用集中参数模型计算的压力响应幅值趋向于0,说明了集中参数模型无法反映高频的响应特性。而分布参数模型和试验结果吻合较好,特别是在ωL/a=π的位置,即一阶纵向声学频率处,气路的末端应处于压力振荡的波腹位置,振幅最大,气路中点处应处于压力振荡的波节位置,振幅最小。本文的模型准确地捕捉到了气路末端的压力响应峰和气路中间的压力波节,峰值和频率均与试验结果吻合。甚至到了二阶声学频率处(ωL/a=2π位置),该模型也有较好的计算结果,二阶声学频率下气路的末端和中点均为压力振荡的波腹,振幅最大,模型的计算结果和试验结果也能够吻合。从相频特性曲线上看,当入口流量的扰动频率增加时,气路中点和末端与入口的相位角滞后是逐渐增大的,相对而言,气路末端的相位角滞后相对更大,可以看到,分布参数模型对相位角的计算与试验结果也能较好吻合。

图6 气路试验系统管路中点压力频率响应特性
Fig.6 Pressure frequency characteristics of the testing apparatus at the middle point of the gas tube

在低频范围内试验数据表明了气路中存在明显的熵波波动特征。但在无量纲角频率ωL/a>1.5后,试验数据沿着比较平滑的曲线移动,这一点在相频曲线中体现得更加明显。这是由于熵波是一种物质波,其传播速度等于流速,文献[20]提到,在通常湍流度水平的圆柱形流路中,当频率大于400 Hz时,熵波会产生耗散效应。原因在于,熵波的产生机理是:当入口推进剂流量产生波动时,在燃烧室内会周期性地生成温度不同的燃气团,当扰动频率较低时,前一时刻生成的燃气团有足够的时间向下游流动,后一时刻在燃烧室头部生成新的燃气团时,前一时刻的燃气团已经流向了下游,各个温度不同的燃气团存在位置空间上的分布,气团随流体整体向下游流动,从而造成了燃烧室内的压力波动。当入口扰动频率很高时,上一个时刻生成的燃气团还没来得及流动到下游,又在同一个位置生成了另一个燃气团,高低温的燃气团几乎在同一位置生成,进而迅速达到热平衡,使得燃气温度在空间分布趋于平均效应,因此无法形成压力波动。在表象上熵波表现出了高频耗散效应,即在低频下,燃烧室的压力响应曲线上会叠加熵波响应,而在高频下,压力响应曲线不会叠加熵波响应,压力幅频响应曲线应趋于平滑。

与文献[8]中的燃烧室线性化频域模型相比较,见图5和图6。线性化频域模型在高频下仍存在叠加的熵波,在相频特性曲线上尤为明显,与上述物理过程不相符合。而本文所建立的时域非线性分布参数模型计算得到的响应曲线在低频处存在熵波叠加,在高频处趋于平滑,能够准确地捕捉到高频下的熵波耗散现象。原因在于,在时域的控制方程中,已经考虑了速度、压力和温度等参数之间的耦合关系,其传播和耗散特性已经在方程中得以体现。因此,只要差分方程和微分方程满足相容性关系,数值积分计算的结果就能够反映原有的物理特性。