CN103116705A - 一种补燃循环火箭发动机故障仿真分析方法 - Google Patents

Abstract

本发明涉及一种补燃循环火箭发动机故障仿真分析方法，包括以下步骤：1】针对补燃循环发动机系统组成及其组件特性，运用流体力学、气体动力学、热力学理论建立模型；2】在模型中嵌入故障特征方程；3】采用线性化的牛顿迭代法求解上述步骤所生成的方程组，得出仿真结果；4】根据上述热力参数特征，通过事先获得的统计参数对发动机热试车或飞行状态可能出现的故障进行试前故障预判，或者通过与发动机试车或飞行后实际测量参数变化趋势进行比对，进行试后故障定位。本发明克服现有技术的不足之处，本发明方法使故障分析的效率、准确性大大提高。

Description

一种补燃循环火箭发动机故障仿真分析方法

技术领域

本发明涉及一种补燃循环火箭发动机故障仿真分析方法。

背景技术

补燃循环发动机是我国新一代运载火箭主要动力装置，其系统较现役常规型号发动机更为复杂，组件种类和数量更多，试车过程检测参数达数十个，氧系统若出现故障后果往往是产品结构严重破坏。以往故障分析工作完全靠设计人员的经验、产品分解检查情况完成，工作量很大，而且受设计人员的个人经验的影响大。本发明可大大提高故障定位的效率和准确性。

发明内容

本发明解决的技术问题是：克服现有技术的不足之处，提供一种补燃循环火箭发动机故障仿真分析方法，该方法使故障分析的效率、准确性大大提高。

本发明的技术解决方案：

一种补燃循环火箭发动机故障仿真分析方法，其特殊之处在于包括下列步骤：

1】针对补燃循环发动机系统组成及其组件特性，运用流体力学、气体动力学、热力学理论建立以下模型：

A 泵静态数学模型：

（1）扬程

$\mathrm{\Δp}=a\·n^2\·\mathrm{\ρ}+b\·n\·q_m+c\·\frac{{q_m}^2}{\mathrm{\ρ}}---\left(1\right)$

式中：Δp-扬程；

a、b、c-通过水力试验获得的扬程特性参数；

n、q_m、ρ-转速、介质流量和密度；

（2）功率

$P=A\·n^3\·\mathrm{\ρ}+B\·n^2\·q_m+C\·n\·\frac{{q_m}^2}{\mathrm{\ρ}}---\left(2\right)$

式中：P-功率；

A、B、C-功率特性参数；

B涡轮静态数学模型：

（1）气涡轮

a涡轮绝热速度

$C_{\mathrm{ad}}=\sqrt{\frac{2k}{k-1}\·R_i\·T_i\·[1-{\left(\frac{p_e}{p_i}\right)}^{\frac{k-1}{k}}]}---\left(3\right)$

式中：C_ad-涡轮绝热速度；

T_i、p_i、p_e-涡轮入口的温度、压力和出口压力；

R_i、k-气体常数、绝热指数；

b涡轮功率

$P={\mathrm{\η}}_t\·q_m\·\frac{k}{k-1}\·R_i\·T_i\·[1-{\left(\frac{p_e}{p_i}\right)}^{\frac{k-1}{k}}]---\left(4\right)$

式中：P-涡轮功率；

T_i、p_i、p_e-涡轮入口的温度、压力和出口压力；

R_i、k-气体常数、绝热指数；

η_t、q_m-涡轮效率、燃气流量；

c涡轮出口温度

$T_e=T_i-T_i\·{\mathrm{\η}}_t\·[1-{\left(\frac{p_e}{p_i}\right)}^{\frac{k-1}{k}}]---\left(5\right)$

式中：T_i、p_i、p_e-涡轮入口的温度、压力和出口压力；

k、T_e-气体绝热指数、涡轮出口温度；

η_t、q_m-涡轮效率、燃气流量；

d涡轮流量

当流动状态为临界或超临界状态时，

即 $\frac{p_e}{p_i}\≤{\left(\frac{2}{k+1}\right)}^{k/(k-1)}$ 时，

$q_m=\frac{\mathrm{\μA}\·p_i}{\sqrt{R_i\·T_i}}\·\sqrt{k\·{\left(\frac{2}{k+1}\right)}^{(k+1)/(k-1)}}---\left(6a\right)$

当流动状态为亚临界时，即 $\frac{p_e}{p_i}>{\left(\frac{2}{k+1}\right)}^{k/(k-1)}$ 时，

$q_m=\frac{\mathrm{\μA}\·p_i}{\sqrt{R_i\·T_i}}\·\sqrt{\frac{2k}{k-1}[{\left(\frac{p_e}{p_i}\right)}^{2/k}-{\left(\frac{p_e}{p_i}\right)}^{(k+1)/k}}---\left(6b\right)$

式中：T_i、p_i、p_e-涡轮入口的温度、压力和出口压力；

k、R_i-气体绝热指数、气体常数；

q_m-燃气流量；

μA-涡轮最小流通面积，可通过涡轮吹风试验获得；

e涡轮效率

${\mathrm{\η}}_t=a_2\·{\left(\frac{n}{C_{\mathrm{ad}}}\right)}^2+a_1\·\left(\frac{n}{C_{\mathrm{ad}}}\right)+a_0---\left(7\right)$

式中：η_t、C_ad-涡轮效率、绝热速度；

a₂,a₁,a₀,n-涡轮的效率特性系数和转速；

（2）液涡轮

a静叶喷嘴出口速度

$C=\sqrt{2\·\frac{(p_i-p_e)}{{\mathrm{\ρ}}_i}}---\left(8\right)$

式中：C-静叶喷嘴出口速度；

p_i、ρ_i-喷嘴入口压力和密度；

p_e-喷嘴出口静压，在数值上等于冲击式涡轮的出口压力，近似等于泵出口压力；

b涡轮流量

    q_m=C·ρ_i·μA              (9)

式中：C-静叶喷嘴出口速度；

ρ_i、μA-流体密度、涡轮静叶喷嘴流通面积；

c涡轮功率

P=η_t·q_m·(p_i-p_e)/ρ_i               (10)

式中：P、η_t、p_i、p_e-涡轮功率、效率、入口压力、出口压力；

ρ_i-流体密度；

C供应管路静态数学模型，包括管道、过滤器、孔板及阀门：

（1）液路元件流阻方程

$\mathrm{\Δp}=\mathrm{\ξ}\·\frac{{q_m}^2}{{\mathrm{\ρ}}_i}---\left(11\right)$

式中：Δp、ρ_i、q_m、ξ-液路元件压降、平均密度、流量和流阻系数；

（2）气路元件流阻方程

${p_i}^2-{p_e}^2=\mathrm{Cg}\·\mathrm{RT}\·{q_m}^2---\left(12\right)$

式中：Cg-气路元件流阻系数；

p_i、p_e、q_m-气路元件入口、出口压力和流量；

R、T-气体介质的气体常数和温度；

D发生器静态数学模型

（1）理论热力特性

根据已有的热力计算软件计算出理论数据表，再使用最近点线性插值法计算实际数据点上的发生器燃气热力特性参数；

$[T,k,R,\mathrm{Cp}]={[T,k,R,\mathrm{Cp}]}_{\mathrm{fore}}+\frac{{[T,k,R,\mathrm{Cp}]}_{\mathrm{back}}-{[T,k,R,\mathrm{Cp}]}_{\mathrm{fore}}}{K_{\mathrm{back}}-K_{\mathrm{fore}}}\·(K-K_{\mathrm{fore}})---\left(13\right)$

式中：T、k、R、Cp-燃气温度、绝热指数、气体常数、定压比热容；

[ ]_fore、[ ]_back-插值前点、插值后点数值；

（2）考虑流体焓值对发生器温度影响

a流体在泵中的温升

$\mathrm{\ΔT}=\frac{\mathrm{\Δp}}{\mathrm{\ρ}\·\mathrm{Cp}}(\frac{1-\mathrm{\η}}{\mathrm{\η}}-\frac{T}{\mathrm{\ρ}}\·\frac{\∂\mathrm{\ρ}}{\∂T})---\left(14\right)$

式中：ΔT、Δp-流体温升、压降；

T、ρ、η、Cp-流体温度、密度、泵效率、流体定压比热容；

-流体密度对温度的偏导数；

b预压泵后流体混合温度

$T=\frac{q_1\·T_1+q_2\·T_2}{q_1+q_2}---\left(15\right)$

式中：T、T₁、T₂、q₁、q₂-流体混合后温度、流体温度1、温度2、流量1、流量2；

c流体焓值对发生器温度的影响

$T_{\mathrm{gg}}=T_{\mathrm{ggth}}+\frac{K_{\mathrm{gg}}\·{\mathrm{Cp}}_o\·(T_o-T_{o0})+{\mathrm{Cp}}_f\·(T_f-T_{f0})}{(K_{\mathrm{gg}}+1)\·{\mathrm{Cp}}_{\mathrm{gg}}}---\left(16\right)$

式中：T_o0、T_f0、T_ggth-理论热力计算时流体温度和对应的发生器理论温度；

Cp_o、Cp_f、T_o、T_f-氧化剂和燃料的定压比热和实际温度；

K_gg、T_gg-发生器混合比、温度；

E推力室静态数学模型

（1）基本方程

a流量

$q_m=q_{\mathrm{mo}}+q_{\mathrm{mf}}=\frac{p_c\·A_t}{C^{*}}---\left(17\right)$

式中：q_m、q_mo、q_mf-推力室燃气流量、氧化剂流量、燃料流量；

p_c、A_t、C^*-推力室压力、喉部流通面积、特征速度；

b混合比

$K_c=\frac{q_{\mathrm{mo}}}{q_{\mathrm{mf}}}---\left(18\right)$

式中：K_c、q_mo、q_mf-推力室混合比、氧化剂流量、燃料流量；

c真空推力

    F_sv＝η_c·η_n·q_m·I_sv          （19）

式中：q_m、F_sv、η_c、I_sv、η_n-推力室燃气流量、真空推力、燃烧效率、真空比冲、喷管效率；

d地面推力

      F=F_sv-p_a·A_e            （20）

式中：p_a、A_e-环境压力和推力室喷口面积；

F、F_sv-推力室地面推力、真空推力；

（2）热力数据处理

根据已有的热力计算软件计算出理论数据表，再使用最近点线性插值法计算实际数据点上的推力室燃气热力特性参数；

$T(\mathrm{Kc},\mathrm{Pc})=T({\mathrm{Tc}}_0,{\mathrm{Pc}}_0)+\frac{\mathrm{\Δ}{T}_{\mathrm{Pc}}}{\mathrm{\ΔPc}}\·(\mathrm{Pc}-{\mathrm{Pc}}_0)+\frac{{\mathrm{\ΔT}}_{\mathrm{Kc}}}{\mathrm{\ΔKc}}\·(\mathrm{Kc}-{\mathrm{Kc}}_0)---\left(21\right)$

式中：Kc0、Pc0、Kc、Pc-插值前点与后点的混合比、室压；-线性插值斜率，由插值前后点数值决定；

F节流组件静态数学模型

（1）节流阀

   Δp=a₀'+a1'·α+a₂′·α²       （22）

式中：Δp、α-节流阀压降、转角；

a0′、a1′、a2′-特性系数，通过液流试验获得；

（2）流量调节器

$q_m=\sqrt{\mathrm{\ρ}}\·(b_0+b_1\·\mathrm{\α})---\left(23\right)$

式中：q_m、ρ、α-调节器流量、介质密度、流量调节器转角；

b0、b1-特性系数，通过液流试验获得；

G平衡方程式

（1）流量平衡方程

∑q_mi=∑q_me          （24）

式中：q_mi、q_me-入口流量、出口流量；

（2）功率平衡方程

P_t=∑P_p           （25）

式中：P_t、P_p-涡轮功率、泵功率；

（3）压力平衡方程

     p_in＝p_out+Δp         （26）

式中：p_in、p_out、Δp-入口压力、出口压力、压降；

H推进剂物理性质

（1）密度

ρ=(a₀₀+a₀₁·T+a₀₂·T²)+(α₁₀+α₁₁·T+α₁₂·T²)·p+(a₂₀+a₂₁·T+a₂₂·T²)·p²             （27）

式中：ρ、p、T-密度、压力、温度；

a₀₀、a₀₁、a₀₂、a₁₀、a₁₁、a₁₂、a₂₀、a₂₁、a₂₂-密度拟合系数；

（2）动力粘度

μ=(b₀₀+b₀₁·T+b₀₂·T²)+(b₁₀+b₁₁·T+b₁₂·T²)·p  （28）

式中：μ、p、T-动力粘度、压力、温度；

b₀₀、b₀₁、b₀₂、b₁₀、b₁₁、b₁₂-粘度拟合系数；

（3）定压比热

Cp=c₀+c₁·T+c₂·T²           （29）

式中：C_p、T-定压比热、温度；

c₀、c₁、c₂-温度拟合系数；

2】嵌入故障特征方程：

所述故障特性方程如下：

A液体或气体泄漏故障

当发动机某个组件出现液体或气体泄漏时，在相应组件模型中补充流量方程；

$\stackrel{\·}{m}={\stackrel{\·}{m}}_{\mathrm{xl}}+{\stackrel{\·}{m}}_{\mathrm{out}}---\left(30\right)$

式中：

入口流量、出口流量、泄漏流量；

B管道、过滤器与孔板堵塞故障

在供应管路模型中增加下述方程，可改变供应管路流阻系数；’

ξ’＝c·ξ                （31）

式中：ξ′、ξ、c-故障状态流阻系数、正常状态流阻系数、常数大于1；

C阀门故障

在供应管路模型中增加下述方程，可改变供应管路流阻系数；

   ξ’=e^α·ξ或 ${\mathrm{\ξ}}^{,}=f\left(\frac{1}{{\sin }^4\mathrm{\α}}\right)\·\mathrm{\ξ}---\left(32\right)$

式中：ξ′、ξ、α-故障状态流阻系数、正常状态流阻系数、阀门打开角度；

D涡轮故障

涡轮故障包括喷嘴烧蚀和效率下降，可在涡轮模型中增加下述方程进行模拟：

喷嘴烧蚀：

   A’=c·A          （33）

式中：A′、A、c-故障状态流通面积、正常状态流通面积、常数大于1，效率下降：

      η’＝c·η            （34）

式中：η′、η、c-故障状态效率、正常状态效率、常数大于0且小于1；E燃气发生器、推力室故障，包括喷嘴堵塞与烧蚀；

通过下述两个方程进行模拟：

喷嘴烧蚀：

    A’=c·A            （35）

式中：A′、A、c-分别为故障状态流通面积、正常状态流通面积、常数大于1；

喷嘴堵塞：

      ξ’=c·ξ         （36）

式中：ξ′、ξ、c-分别为故障状态流阻系数、正常状态流阻系数、常数大于1；

F泵故障

泵故障主要包括转子损坏、碰磨、汽蚀或泵叶片断裂，均通过改变扬程或效率来模拟；

η’＝c·η或ΔP’=c·ΔP         （37）

式中：η′、η、ΔP′、ΔP、c-分别为故障状态效率、正常状态效率、故障状态扬程、正常状态扬程、常数大于0且小于1；

3】采用线性化的牛顿迭代法求解上述步骤所生成的方程组，得出仿真结果；所述仿真结果为发动机在不同故障模式下的热力参数特征；

4】根据上述热力参数特征，通过事先获得的统计参数对发动机热试车或飞行状态可能出现的故障进行试前故障预判，或者通过与发动机试车或飞行后实际测量参数变化趋势进行比对，进行试后故障定位。

本发明与现有技术相比的有益效果是：

1、可在热试车或飞行前进行故障预测，采取相应的规避措施，降低风险。

2、本发明能够实现发动机热试车或飞行后故障分析与定位，提高故障分析效率与准确性，减少人力、物力成本，降低研制费用。

3、本发明能较为全面地模拟补燃发动机典型故障模式，解决了无法通过试验确定故障特性的难题。

4、本发明为准确评估发动机工作状态、分析故障发生的部位和原因提供更为客观、便捷的方法，大大提高了分析的科学性和效率。

附图说明

图1为本发明发动机系统组成结构示意图；

图2为本发明涡轮效率下降故障模式下的仿真计算得到的热力参数特征与热试车实测的参数特征对比示意图；

其中附图标记为：1-涡轮；2-氧化剂泵；3-推力室；4-推力室燃料阀；5-节流阀；6-燃料泵；7-流量调节器；8-发生器；9-发生器氧化剂阀；10-发生器燃料阀。

具体实施方式

具体方法如下：

（1）依据发动机具体系统组成和组件特性，运用上述建模理论建立发动机稳态工作数学模型；

（2）根据需要模拟的故障模式，在发动机工作模型中相应部位嵌入故障特征方程；

（3）运用特定的求解方法求解上述发动机工作模型和故障特性方程所组成的非线性方程组；

（4）上述非线性方程组的求解结果即为特定故障模式下的发动机热力参数特征。该参数特征可在一定程度上表征热试或飞行前发动机的故障影响，另外，在热试或飞行参数出现异常时，可通过模拟参数变化特征，对故障进行初步定位。

实施例1：

1、依据发动机具体系统组成（见图1）和组件特性，建立发动机稳态工作数学模型。

A泵静态数学模型：

（1）扬程

$\mathrm{\Δp}=a\·n^2\·\mathrm{\ρ}+b\·n\·q_m+c\·\frac{{q_m}^2}{\mathrm{\ρ}}---\left(1\right)$

式中：Δp—扬程；

a、b、c—扬程特性参数（通过水力试验获得）；

n、q_m、ρ—转速、介质流量和密度。

（2）功率

$P=A\·n^3\·\mathrm{\ρ}+B\·n^2\·q_m+C\·n\·\frac{{q_m}^2}{\mathrm{\ρ}}---\left(2\right)$

式中：P—功率；

A、B、C——功率特性参数。

B涡轮静态数学模型：

a．涡轮绝热速度

$C_{\mathrm{ad}}=\sqrt{\frac{2k}{k-1}\·R_i\·T_i\·[1-{\left(\frac{p_e}{p_i}\right)}^{\frac{k-1}{k}}]}---\left(3\right)$

式中：C_ad—涡轮绝热速度；

T_i、p_i、p_e—涡轮入口的温度、压力和出口压力；

R_i、k—气体常数、绝热指数。

b．涡轮功率

$P={\mathrm{\η}}_t\·q_m\·\frac{k}{k-1}\·R_i\·T_i\·[1-{\left(\frac{p_e}{p_i}\right)}^{\frac{k-1}{k}}]---\left(4\right)$

式中：P—涡轮功率；

T_i、p_i、p_e—涡轮入口的温度、压力和出口压力；

R_i、k—气体常数、绝热指数；

η_t、q_m—涡轮效率、燃气流量。

c．涡轮出口温度

$T_e=T_i-T_i\·{\mathrm{\η}}_t\·[1-{\left(\frac{p_e}{p_i}\right)}^{\frac{k-1}{k}}]---\left(5\right)$

式中：T_i、p_i、p_e—涡轮入口的温度、压力和出口压力；

k、T_e—气体绝热指数、涡轮出口温度；

η_t、q_m—涡轮效率、燃气流量。

d．涡轮流量

当流动状态为临界或超临界状态时，即

时，

$q_m=\frac{\mathrm{\μA}\·p_i}{\sqrt{R_i\·T_i}}\·\sqrt{k\·{\left(\frac{2}{k+1}\right)}^{(k+1)/(k-1)}}---\left(6\right)$

当流动状态为亚临界时，即

时，

$q_m=\frac{\mathrm{\μA}\·p_i}{\sqrt{R_i\·T_i}}\·\sqrt{\frac{2k}{k-1}[{\left(\frac{p_e}{p_i}\right)}^{2/k}-{\left(\frac{p_e}{p_i}\right)}^{(k+1)/k}}---\left(7\right)$

式中：T_i、p_i、p_e—涡轮入口的温度、压力和出口压力；

k、R_i—气体绝热指数、气体常数；

q_m—燃气流量；

μA——涡轮最小流通面积，可通过涡轮吹风试验获得。

C供应管路静态数学模型，包括管道与阀门。

（1）液路元件流阻方程

$\mathrm{\Δp}=\mathrm{\ξ}\·\frac{{q_m}^2}{{\mathrm{\ρ}}_i}---\left(8\right)$

式中：Δp、ρ_i、q_m、ξ—液路元件压降，平均密度（用入口密度近似代替）、流量和流阻系数。

（2）气路元件流阻方程

${p_i}^2-{p_e}^2=\mathrm{Cg}\·\mathrm{RT}\·{q_m}^2---\left(9\right)$

式中：Cg—气路元件流阻系数；

p_i、p_e、q_m—气路元件入口、出口压力和流量；

R、T—气体介质的气体常数和温度（用入口参数）。

D发生器静态数学模型

（1）理论热力特性

根据已有的热力计算软件计算出理论数据表，再使用最近点线性插值法计算实际数据点上的发生器燃气热力特性参数。

$[T,k,R,\mathrm{Cp}]={[T,k,R,\mathrm{Cp}]}_{\mathrm{fore}}+\frac{{[T,k,R,\mathrm{Cp}]}_{\mathrm{back}}-{[T,k,R,\mathrm{Cp}]}_{\mathrm{fore}}}{K_{\mathrm{back}}-K_{\mathrm{fore}}}\·(K-K_{\mathrm{fore}})---\left(10\right)$

式中：T、k、R、Cp—燃气温度、绝热指数、气体常数、定压比热容；

[ ]_fore、[ ]_back——插值前点、插值后点数值。

（2）考虑流体焓值对发生器温度影响

a．流体在泵中的温升

$\mathrm{\ΔT}=\frac{\mathrm{\Δp}}{\mathrm{\ρ}\·\mathrm{Cp}}(\frac{1-\mathrm{\η}}{\mathrm{\η}}-\frac{T}{\mathrm{\ρ}}\·\frac{\∂\mathrm{\ρ}}{\∂T})---\left(11\right)$

式中：ΔT、Δp—流体温升、压降；

T、ρ、η、Cp—流体温度、密度、泵效率、流体定压比热容；

—流体密度对温度的偏导数。

c．流体焓值对发生器温度的影响

$T_{\mathrm{gg}}=T_{\mathrm{ggth}}+\frac{K_{\mathrm{gg}}\·{\mathrm{Cp}}_o\·(T_o-T_{o0})+{\mathrm{Cp}}_f\·(T_f-T_{f0})}{(K_{\mathrm{gg}}+1)\·{\mathrm{Cp}}_{\mathrm{gg}}}---\left(12\right)$

式中：T_o0、T_f0、T_ggth—理论热力计算时流体温度和对应的发生器理论温度；

Cp_o、Cp_f、T_o、T_f—氧化剂和燃料的定压比热和实际温度；

K_gg、T_gg—发生器混合比、温度。

E推力室静态数学模型

（1）基本方程

a.流量

$q_m=q_{\mathrm{mo}}+q_{\mathrm{mf}}=\frac{p_c\·A_t}{C^{*}}---\left(13\right)$

式中：q_m、q_mo、q_mf—推力室燃气流量、氧化剂流量、燃料流量；

p_c、A_t、C^*—推力室压力、喉部流通面积、特征速度。

b．混合比

$K_c=\frac{q_{\mathrm{mo}}}{q_{\mathrm{mf}}}---\left(14\right)$

式中：K_c、q_mo、q_mf—推力室混合比、氧化剂流量、燃料流量；

c．真空推力

    F_sv＝η_c·η_n·q_m·I_sv              （15）

式中：q_m、F_sv、η_c、I_sv、η_n—推力室燃气流量、真空推力、燃烧效率、真空比冲、喷管效率。

d．地面推力

F=F_sv-p_a·A_e             （16）

式中：p_a、A_e—环境压力和推力室喷口面积；

F、F_sv—推力室地面推力、真空推力。

（2）热力数据处理

根据已有的热力计算软件计算出理论数据表，再使用最近点线性插值法计算实际数据点上的推力室燃气热力特性参数。

$T(\mathrm{Kc},\mathrm{Pc})=T({\mathrm{Kc}}_0,{\mathrm{Pc}}_0)+\frac{{\mathrm{\ΔT}}_{\mathrm{Pc}}}{\mathrm{\ΔPc}}\·(\mathrm{Pc}-{\mathrm{Pc}}_0)+\frac{{\mathrm{\ΔT}}_{\mathrm{Kc}}}{\mathrm{\ΔKc}}\·(\mathrm{Kc}-{\mathrm{Kc}}_0)---\left(17\right)$

式中：Kc0、Pc0、Kc、Pc—插值前点与后点的混合比、室压；

—线性插值斜率，由插值前后点数值决定。

F节流组件静态数学模型

（1）节流阀

   Δp＝a₀'+a₁'·α+a₂'·α²             （18）

式中：Δp、α—节流阀压降、转角；

a0′、a1′、a2′-特性系数，通过液流试验获得。

（2）流量调节器

$q_m=\sqrt{\mathrm{\ρ}}\·(b_0+b_1\·\mathrm{\α})---\left(19\right)$

式中：q_m、ρ、α——调节器流量、介质密度、流量调节器转角；

b0、b1—特性系数，通过液流试验获得。

G平衡方程式

（1）流量平衡方程

    ∑q_mi=∑q_me          （20）

式中：q_mi、q_me—入口流量、出口流量。

（2）功率平衡方程

   P_t=∑P_p             （21）

式中：P_t、P_p—涡轮功率、泵功率。

（3）压力平衡方程

    p_in＝p_out+Δ_p          （22）

式中：p_in、p_out、Δp—入口压力、出口压力、压降。

H推进剂物理性质

（1）密度

ρ=(a₀₀+a₀₁·T+a₀₂·T²)+(a₁₀+a₁₁·T+a₁₂·T²)·p+(a₂₀+a₂₁·T+a₂₂·T²)·p²                 （23）

式中：ρ、p、T—密度、压力、温度；

a₀₀、a₀₁、a₀₂、a₁₀、a₁₁、a₁₂、a₂₀、a₂₁、a₂₂—拟合系数。

（2）动力粘度

μ=(b₀₀+b₀₁·T+b₀₂·T²)+(b₁₀+b₁₁·T+b₁₂·T²)·p    （24）

式中：μ、p、T—动力粘度、压力、温度；

b₀₀、b₀₁、b₀₂、b₁₀、b₁₁、b₁₂—拟合系数。

（3）定压比热

Cp＝c₀+c₁·T+c₂·T²            （25）

式中：Cp、T-定压比热、温度；

c₀、c₁、c₂—拟合系数。

2、计算涡轮效率下降与推力室堵塞为例：

根据需要模拟的故障模式在发动机稳态工作数学模型的相应部位嵌入故障特征方程，仿真计算涡轮效率下降与推力室堵塞两种故障模式。

（1）涡轮效率下降

       η’＝c·η

式中：η′、η、c—故障状态效率积、正常状态效率、常数（大于0且小于1）。

（2）推力室堵塞

          ξ’=c·ξ

式中：ξ′、ξ、c—故障状态流阻系数、正常状态流阻系数、常数（大于1）。

3、涡轮效率下降故障模式下的仿真计算得到的热力参数特征与热试车实测的参数特征对比情况如图2所示。其中，纵轴参数代表故障状态参数相对正常参数的相对变化量，斜杠线柱图代表热试车实测参数，点状柱图代表仿真计算参数。其中：A-发动机推力；B-发动机氧化剂流量；C-发动机燃料流量；D-氧化剂泵出口压力；E-燃料一级泵出口压力；F-燃料二级泵出口压力；G-发生器压力；H-涡轮出口压力；I-主涡轮泵转速。

4、上述参数相对变化趋势对比结果表明，本次试车参数异常与涡轮效率下降相关，可初步定位为涡轮效率下降故障。

实施例2：

典型故障模式仿真、效应特性分析

1、依据发动机具体系统组成（见图1）和组件特性，建立发动机稳态工作数学模型。

2、对如下几种典型故障模式进行仿真、效应特性分析得出故障模式效应表，详见表1。

表1故障模式效应表

参数的故障敏感性定义为：

$f\left(x\right)\left\{\begin{array}{cc}1& \mathrm{\&Delta;a}\%>1\%\\ 0& -1\%\&le;\mathrm{\&Delta;a}\%\&le;1\%\\ -1& \mathrm{\&Delta;a}\%<-1\%\end{array}\right.$

其中Δa%为参数相对变化量，上述定义对应堵塞压降增加5倍，流量泄漏、面积增加、效率下降50%的情形。

如果由上述十个参数形成故障效应，可表示为向量（发动机推力、发动机氧化剂流量、发动机燃料流量、氧化剂泵出口压力、燃料一级泵出口压力、燃料泵二级泵出口压力、发生器压力、涡轮出口压力、主涡轮泵转速），根据参数的故障敏感性定义，正影响记为1，负影响记为-1，无影响或弱影响记为0。

3、通过上述故障模式仿真与效应特性分析，可以建立特定类型发动机的故障模式效应表，形成效应特性对比分析数据集。对故障预测与确定故障特性起到了一定的积极作用。

Claims (2)

1.一种补燃循环火箭发动机故障仿真分析方法，其特征在于包括下列步骤：

1】针对补燃循环发动机系统组成及其组件特性，运用流体力学、气体动力学、热力学理论建立以下模型：

A泵静态数学模型：

（1）扬程

$\mathrm{\&Delta;p}=a\&CenterDot;n^2\&CenterDot;\mathrm{\&rho;}+b\&CenterDot;n\&CenterDot;q_m+c\&CenterDot;\frac{{q_m}^2}{\mathrm{\&rho;}}---\left(1\right)$

式中：Δp-扬程；

a、b、c-通过水力试验获得的扬程特性参数；

n、q_m、ρ-转速、介质流量和密度；

（2）功率

$P=A\&CenterDot;n^3\&CenterDot;\mathrm{\&rho;}+B\&CenterDot;n^2\&CenterDot;q_m+C\&CenterDot;n\&CenterDot;\frac{{q_m}^2}{\mathrm{\&rho;}}---\left(2\right)$

式中：P-功率；

A、B、C-功率特性参数；

B涡轮静态数学模型：

（1）气涡轮

a涡轮绝热速度

$C_{\mathrm{ad}}=\sqrt{\frac{2k}{k-1}\&CenterDot;R_i\&CenterDot;T_i\&CenterDot;[1-{\left(\frac{p_e}{p_i}\right)}^{\frac{k-1}{k}}]}---\left(3\right)$

式中：C_ad-涡轮绝热速度；

T_i、p_i、p_e-涡轮入口的温度、压力和出口压力；

R_i、k-气体常数、绝热指数；

b涡轮功率

$P={\mathrm{\&eta;}}_t\&CenterDot;q_m\&CenterDot;\frac{k}{k-1}\&CenterDot;R_i\&CenterDot;T_i\&CenterDot;[1-{\left(\frac{p_e}{p_i}\right)}^{\frac{k-1}{k}}]---\left(4\right)$

式中：P-涡轮功率；

T_i、p_i、p_e-涡轮入口的温度、压力和出口压力；

R_i、k-气体常数、绝热指数；

η_t、q_m-涡轮效率、燃气流量；

c涡轮出口温度

$T_e=T_i-T_i\&CenterDot;{\mathrm{\&eta;}}_t\&CenterDot;[1-{\left(\frac{p_e}{p_i}\right)}^{\frac{k-1}{k}}]---\left(5\right)$

式中：T_i、p_i、p_e-涡轮入口的温度、压力和出口压力；

k、T_e-气体绝热指数、涡轮出口温度；

η_t、q_m-涡轮效率、燃气流量；

d涡轮流量

当流动状态为临界或超临界状态时，

即 $\frac{p_e}{p_i}\&le;{\left(\frac{2}{k+1}\right)}^{k/(k-1)}$ 时，

$q_m=\frac{\mathrm{\&mu;A}\&CenterDot;p_i}{\sqrt{R_i\&CenterDot;T_i}}\&CenterDot;\sqrt{k\&CenterDot;{\left(\frac{2}{k+1}\right)}^{(k+1)/(k-1)}}---\left(6a\right)$

当流动状态为亚临界时，即 $\frac{p_e}{p_i}>{\left(\frac{2}{k+1}\right)}^{k/(k-1)}$ 时，

$q_m=\frac{\mathrm{\&mu;A}\&CenterDot;p_i}{\sqrt{R_i\&CenterDot;T_i}}\&CenterDot;\sqrt{\frac{2k}{k-1}[{\left(\frac{p_e}{p_i}\right)}^{2/k}-{\left(\frac{p_e}{p_i}\right)}^{(k+1)/k}}---\left(6b\right)$

式中：T_i、p_i、p_e-涡轮入口的温度、压力和出口压力；

k、R_i-气体绝热指数、气体常数；

q_m-燃气流量；

μA-涡轮最小流通面积，可通过涡轮吹风试验获得；

e涡轮效率

${\mathrm{\&eta;}}_t=a_2\&CenterDot;{\left(\frac{n}{C_{\mathrm{ad}}}\right)}^2+a_1\&CenterDot;\left(\frac{n}{C_{\mathrm{ad}}}\right)+a_0---\left(7\right)$

式中：η_t、C_ad-涡轮效率、绝热速度；

a₂,a₁,a₀,n-涡轮的效率特性系数和转速；

（2）液涡轮

a静叶喷嘴出口速度

$C=\sqrt{2\&CenterDot;\frac{(p_i-p_e)}{{\mathrm{\&rho;}}_i}}---\left(8\right)$

式中：C-静叶喷嘴出口速度；

p_i、ρ_i-喷嘴入口压力和密度；

p_e-喷嘴出口静压，在数值上等于冲击式涡轮的出口压力，近似等于泵出口压力；

b涡轮流量

    q_m=C·ρ_i·μA          (9)

式中：C-静叶喷嘴出口速度；

ρ_i、μA-流体密度、涡轮静叶喷嘴流通面积；

c涡轮功率

   P＝η_t·q_m·(p_i-p_e)/ρ_i       (10)

式中：P、η_t、p_i、p_e-涡轮功率、效率、入口压力、出口压力；

ρ_i-流体密度；

C供应管路静态数学模型，包括管道、过滤器、孔板及阀门：

（1）液路元件流阻方程

$\mathrm{\&Delta;p}=\mathrm{\&xi;}\&CenterDot;\frac{{q_m}^2}{{\mathrm{\&rho;}}_i}---\left(11\right)$

式中：Δp、ρ_i、q_m、ξ-液路元件压降、平均密度、流量和流阻系数；

（2）气路元件流阻方程

${p_i}^2-{p_e}^2=\mathrm{Cg}\&CenterDot;\mathrm{RT}\&CenterDot;{q_m}^2---\left(12\right)$

式中：Cg-气路元件流阻系数；

p_i、p_e、q_m-气路元件入口、出口压力和流量；

R、T-气体介质的气体常数和温度；

D发生器静态数学模型

（1）理论热力特性

根据已有的热力计算软件计算出理论数据表，再使用最近点线性插值法计算实际数据点上的发生器燃气热力特性参数；

$[T,k,R,\mathrm{Cp}]={[T,k,R,\mathrm{Cp}]}_{\mathrm{fore}}+\frac{{[T,k,R,\mathrm{Cp}]}_{\mathrm{back}}-{[T,k,R,\mathrm{Cp}]}_{\mathrm{fore}}}{K_{\mathrm{back}}-K_{\mathrm{fore}}}\&CenterDot;(K-K_{\mathrm{fore}})---\left(13\right)$

式中：T、k、R、Cp-燃气温度、绝热指数、气体常数、定压比热容；

[ ]_fore、[ ]_back-插值前点、插值后点数值；

（2）考虑流体焓值对发生器温度影响

a流体在泵中的温升

$\mathrm{\&Delta;T}=\frac{\mathrm{\&Delta;p}}{\mathrm{\&rho;}\&CenterDot;\mathrm{Cp}}(\frac{1-\mathrm{\&eta;}}{\mathrm{\&eta;}}-\frac{T}{\mathrm{\&rho;}}\&CenterDot;\frac{\&PartialD;\mathrm{\&rho;}}{\&PartialD;T})---\left(14\right)$

式中：ΔT、Δp-流体温升、压降；

T、ρ、η、Cp-流体温度、密度、泵效率、流体定压比热容；

-流体密度对温度的偏导数；

b预压泵后流体混合温度

$T=\frac{q_1\&CenterDot;T_1+q_2\&CenterDot;T_2}{q_1+q_2}---\left(15\right)$

式中：T、T₁、T₂、q₁、q₂-流体混合后温度、流体温度1、温度2、流量1、流量2；

c流体焓值对发生器温度的影响

$T_{\mathrm{gg}}=T_{\mathrm{ggth}}+\frac{K_{\mathrm{gg}}\&CenterDot;{\mathrm{Cp}}_o\&CenterDot;(T_o-T_{o0})+{\mathrm{Cp}}_f\&CenterDot;(T_f-T_{f0})}{(K_{\mathrm{gg}}+1)\&CenterDot;{\mathrm{Cp}}_{\mathrm{gg}}}---\left(16\right)$

式中：T_o0、T_f0、T_ggth-理论热力计算时流体温度和对应的发生器理论温度；

Cp_o、Cp_f、T_o、T_f-氧化剂和燃料的定压比热和实际温度；

K_gg、T_gg-发生器混合比、温度；

E推力室静态数学模型

（1）基本方程

a流量

$q_m=q_{\mathrm{mo}}+q_{\mathrm{mf}}=\frac{p_c\&CenterDot;A_t}{C^{*}}---\left(17\right)$

式中：q_m、q_mo、q_mf-推力室燃气流量、氧化剂流量、燃料流量；

p_c、A_t、C^*-推力室压力、喉部流通面积、特征速度；

b混合比

$K_c=\frac{q_{\mathrm{mo}}}{q_{\mathrm{mf}}}---\left(18\right)$

式中：K_c、q_mo、q_mf-推力室混合比、氧化剂流量、燃料流量；

c真空推力

    F_sv＝η_c·η_n·q_m·I_sv      （19）

式中：q_m、F_sv、η_c、I_sv、η_n-推力室燃气流量、真空推力、燃烧效率、真空比冲、喷管效率；

d地面推力

     F=F_sv-p_a·A_e       （20）

式中：p_a、A_e-环境压力和推力室喷口面积；

F、F_sv-推力室地面推力、真空推力；

（2）热力数据处理

根据已有的热力计算软件计算出理论数据表，再使用最近点线性插值法计算实际数据点上的推力室燃气热力特性参数；

$T(\mathrm{Kc},\mathrm{Pc})=T({\mathrm{Tc}}_0,{\mathrm{Pc}}_0)+\frac{\mathrm{\&Delta;}{T}_{\mathrm{Pc}}}{\mathrm{\&Delta;Pc}}\&CenterDot;(\mathrm{Pc}-{\mathrm{Pc}}_0)+\frac{{\mathrm{\&Delta;T}}_{\mathrm{Kc}}}{\mathrm{\&Delta;Kc}}\&CenterDot;(\mathrm{Kc}-{\mathrm{Kc}}_0)---\left(21\right)$

式中：Kc0、Pc0、Kc、Pc-插值前点与后点的混合比、室压；

-线性插值斜率，由插值前后点数值决定；

F节流组件静态数学模型

（1）节流阀

   Δp=a₀'+a₁'·α+a₂'·α²         （22）

式中：Δp、α-节流阀压降、转角；

a0′、a1′、a2′-特性系数，通过液流试验获得；

（2）流量调节器

$q_m=\sqrt{\mathrm{\&rho;}}\&CenterDot;(b_0+b_1\&CenterDot;\mathrm{\&alpha;})---\left(23\right)$

式中：q_m、ρ、α-调节器流量、介质密度、流量调节器转角；

b0、b1-特性系数，通过液流试验获得；

G平衡方程式

（1）流量平衡方程

    ∑q_mi=∑q_me        （24）

式中：q_mi、q_me-入口流量、出口流量；

（2）功率平衡方程

P_t=∑P_p            （25）

式中：P_t、P_p-涡轮功率、泵功率；

（3）压力平衡方程

    p_in＝p_out+Δp         （26）

式中：p_in、p_out、Δp-入口压力、出口压力、压降；

H推进剂物理性质

（1）密度

ρ＝(a₀₀+a₀₁·T+a₀₂·T²)+(a₁₀+a₁₁·T+a₁₂·T²)·p+(a₂₀+a₂₁·T+a₂₂·T²)·p²         （27）

式中：ρ、p、T-密度、压力、温度；

a₀₀、a₀₁、a₀₂、a₁₀、a₁₁、a₁₂、a₂₀、a₂₁、a₂₂-密度拟合系数；

（2）动力粘度

μ=(b₀₀+b₀₁·T+b₀₂·T²)+(b₁₀+b₁₁·T+b₁₂·T²)·p     （28）

式中：μ、p、T-动力粘度、压力、温度；

b₀₀、b₀₁、b₀₂、b₁₀、b₁₁、b₁₂-粘度拟合系数；

（3）定压比热

Cp=c₀+c₁·T+c₂·T²         （29）

式中：Cp、T-定压比热、温度；

c₀、c₁、c₂-温度拟合系数；

2】嵌入故障特征方程：

所述故障特性方程如下：

A液体或气体泄漏故障

当发动机某个组件出现液体或气体泄漏时，在相应组件模型中补充流量方程；

$\stackrel{\&CenterDot;}{m}={\stackrel{\&CenterDot;}{m}}_{\mathrm{xl}}+{\stackrel{\&CenterDot;}{m}}_{\mathrm{out}}---\left(30\right)$

式中：

-入口流量、出口流量、泄漏流量；

B管道、过滤器与孔板堵塞故障

在供应管路模型中增加下述方程，可改变供应管路流阻系数；

ξ’＝c·ξ        （31）

式中：ξ′、ξ、c-故障状态流阻系数、正常状态流阻系数、常数大于1；

C阀门故障

在供应管路模型中增加下述方程，可改变供应管路流阻系数；

ξ’=e^α·ξ或 ${\mathrm{\&xi;}}^{,}=f\left(\frac{1}{{\sin }^4\mathrm{\&alpha;}}\right)\&CenterDot;\mathrm{\&xi;}---\left(32\right)$

式中：ξ′、ξ、α-故障状态流阻系数、正常状态流阻系数、阀门打开角度；

D涡轮故障

涡轮故障包括喷嘴烧蚀和效率下降，可在涡轮模型中增加下述方程进行模拟：

喷嘴烧蚀：

     A’=c·A     （33）

式中：A′、A、c-故障状态流通面积、正常状态流通面积、常数大于1，效率下降：

     η’＝c·η       （34）

式中：η′、η、c-故障状态效率、正常状态效率、常数大于0且小于1；

E燃气发生器、推力室故障，包括喷嘴堵塞与烧蚀；

通过下述两个方程进行模拟：

喷嘴烧蚀：

    A’=c·A       （35）

式中：A′、A、c-分别为故障状态流通面积、正常状态流通面积、常数大于1；

喷嘴堵塞：

   ξ’＝c·ξ         （36）

式中：ξ′、ξ、c-分别为故障状态流阻系数、正常状态流阻系数、常数大于1；

F泵故障

泵故障主要包括转子损坏、碰磨、汽蚀或泵叶片断裂，均通过改变扬程或效率来模拟；

    η’＝c·η或ΔP’=c·ΔP         （37）

式中：η′、η、ΔP′、ΔP、c-分别为故障状态效率、正常状态效率、故障状态扬程、正常状态扬程、常数大于0且小于1；

3】采用线性化的牛顿迭代法求解上述步骤所生成的方程组，得出仿真结果；所述仿真结果为发动机在不同故障模式下的热力参数特征；

4】根据上述热力参数特征，通过事先获得的统计参数对发动机热试车或飞行状态可能出现的故障进行试前故障预判，或者通过与发动机试车或飞行后实际测量参数变化趋势进行比对，进行试后故障定位。

2.根据权利要求1所述的补燃循环火箭发动机故障仿真分析方法，其特征在于：所述步骤4】的具体实现方式为：统计实际测量参数变化，后与故障模式效应表进行比对，进行试前故障预判或试后故障定位。

Images