一种基于空间推进算法的旋转爆震发动机流场参数分布预测方法与流程

1.本发明属于旋转爆震发动机流场数值计算领域，具体涉及一种基于空间推进算法的旋转爆震发动机流场参数分布预测方法。


背景技术：

2.旋转爆震发动机作为一种新型推进装置，如图2所示，常见的结构为同心环柱，环柱底端开有可燃混合气喷注槽口。环腔内存在单个或多个爆震波朝喷注平面倾斜，且沿周向传播。紧邻爆震波后的气体经过迅速释热，温度和压力得到显著提升，而后一系列膨胀波使气流膨胀加速。当靠近喷注平面的燃烧产物的压力低于可燃混合气的喷注总压时，可燃混合气得以注入环腔。爆震波和可燃混合气/燃烧产物间断面，围成一个充满可燃混合气三角区，确保爆震燃烧的可持续性。在流场结构上，爆震波、剪切层以及诱导激波相交于三角区上顶点。大部分气流沿着旋转爆震发动机出口轴向排出，但流动参数沿周向的变化仍显著。
3.相较于燃气涡轮发动机和冲压发动机，旋转爆震发动机在工作中熵增较小，热效率更高，受到各研究机构的广泛关注。为了获取旋转爆震发动机流场的完整信息，高精度数值仿真技术必不可少。爆震燃烧过程中运动激波与化学放热反应强烈耦合，导致传统的基于时间推进的数值模拟网格尺度和离散时间步过小。为了获取旋转爆震发动机的气动性能，往往耗费较多的时间和计算资源。


技术实现要素：

4.技术目的：为解决上述技术问题，本发明提供了一种基于空间推进算法的旋转爆震发动机流场参数分布预测方法，通过耦合空间推进算法和单一γ的cj爆震波模型，并利用空间推进算法所对应的流线网格的随流特性，确定激波坐标系下爆震波后流动参数以及爆震波远下游流动参数，并通过坐标转换获得实验室坐标系下的流动参数，有效降低旋转爆震发动机流场解算耗时，实现旋转爆震发动机流动参数分布的快速预测。
5.技术方案：为实现上述目的，本发明采用的技术方案为：
6.一种基于空间推进算法的旋转爆震发动机流场参数分布预测方法，其特征在于：包括以下步骤：
7.(1)、预先估计喷注面进口平均速度，根据可压缩流动关系式，确定爆震波波前流动参数；
8.(2)、将爆震波前流动参数代入单一比热比的cj爆震模型中，确定爆震波后流动参数；
9.(3)、将爆震波后流动参数以及远下游流动参数赋值给初值线；第一次迭代时爆震波远下游流动参数由爆震波后压力膨胀至爆震波前压力获得，此后迭代中的爆震波远下游的流动参数则根据前一次计算结果；
10.(4)、将初值线代入空间推进算法中，确定激波坐标系下的旋转爆震发动机流场结构及流动参数分布；
11.(5)、将喷注面解点的压力值与喷注总压相比较，当解点压力值小于喷注总压时，对喷注流场和爆震波下游流场进行耦合求解，当解点压力值大于喷注总压时，仅求解爆震波下游流场；
12.(6)、确定爆震波远下游位置的爆震波前的流动参数的面积平均值，返回步骤(2)-步骤(5)进行迭代，直至远下游的平均气流角度为0，跳出循环，执行步骤(7)；
13.(7)、根据爆震波倾角大小，确定激波坐标系的旋转角度，并根据速度三角形确定实验室坐标系下的旋转爆震发动机流动参数分布。
14.进一步的，所述步骤(1)包括如下步骤：
15.(11)、根据旋转爆震发动机飞行工况，确定旋转爆震进口可燃混合气的总压pt和总温tt，并估计旋转爆震发动机进口平均速度
16.(12)、根据可压缩流关系式，确定爆震波前平均静温平均静压以及平均密度
17.进一步的，所述步骤(2)包括如下步骤：
18.(21)、构造单一比热比的cj爆震模型，模型方程为如下形式：
[0019][0020][0021]
其中，m
cj
表示爆震波马赫数，h表示无量纲释热量，q表示化学反应热，rg表示气体常数，γ为可燃混合气和燃烧产物的比热比；
[0022]
(22)、根据跨爆震波前后质量守恒定律、动量守恒定律以及能量守恒定律，确定爆震波后流动参数，包括爆震波后静压爆震波后密度和爆震波后温度
[0023]
进一步的，所述步骤(4)包括如下步骤：
[0024]
(41)、通过主流方向无粘通量e确定中间变量χ，进而获得初值线上的状态量q，中间变量χ的表达式为如下形式：
[0025][0026]
(42)、根据初值线上速度分布，确定下游解点坐标建立流线网格，表示x方向上第j个下游解点的x坐标，表示x方向上第j个初值线点的x坐标；
[0027]
(43)、根据本质无振荡插值方法和限制函数，获得步进半步的沿x方向的和的数值解，进而获取步进长度为1/2位置的状态量和
[0028]
(44)、根据激波极曲线理论，通过空间x方向上步进长度为1/2设定值的相邻两点，确定空间推进算法所对应的步进长度为1/2位置的解点压力和气流角，进而获得黎曼自相似解，相邻两点状态曲线的表达式为如下形式：
[0029][0030][0031]
其中，φb表示初值线下侧点的状态曲线，θb表示初值线下侧点的气流角，f(mb,αb)表示初值线下侧点的兰金许贡关系式，αb表示下游空间中任一点压力与初值线下侧点压力之比，ν(mb)表示初值线的下侧点的普朗特-迈耶函数，φ
t
表示初值线上侧点的状态曲线，θ
t
表示初值线上侧点的气流角，f(m
t
,α
t
)表示初值线上侧点的兰金许贡关系式，α
t
表示下游空间任一点压力与初值线上侧点压力之比，ν(m
t
)表示初值线的上侧点的普朗特-迈耶函数；
[0032]
(45)、通过对可压缩欧拉方程进行差分离散，获得单个空间递进步长下的流动参数，下游沿流动x方向的差分方程为如下形式：
[0033][0034]
其中，表示初值线相邻两点的y坐标差值，表示下游相邻两解点的y坐标差值，表示空间步进1/2长度位置的上侧格点的沿y方向的无粘通量，表示空间步进1/2长度位置的上侧格边的斜率，表示空间步进1/2长度位置的上侧格点的沿x方向上的无粘通量，表示空间步进1/2长度位置的下侧格点的沿y方向上的无粘通量，表示空间步进1/2长度位置的下侧格边的斜率，表示空间步进1/2长度位置的下侧格点的沿x方向上的无粘通量。
[0035]
进一步的，所述步骤(6)包括如下步骤：
[0036]
(61)、根据远下游沿y方向的流动参数分布，识别爆震波前可燃混合气所在区间；
[0037]
(62)、确定爆震波前可燃混气所在区间的流动参数的面积平均值，面积平均公式为如下形式：
[0038][0039][0040][0041][0042]
其中，表示爆震波前平均速度的x分量，表示爆震波前平均速度的y分量，n
x
表示爆震波面法矢的x分量，ny表示爆震波面法矢的y分量，表示单位质量所含内能，根据牛
顿-拉弗森多元迭代法，联立求解面积平均公式的四个方程，获得爆震波前流动参数以及
[0043]
进一步的，所述步骤(7)包括如下步骤：
[0044]
(71)、根据爆震波前气流角，修正爆震波相对于喷注平面的倾角，修正公式为如下形式：
[0045][0046]
其中，表示经过修正的爆震波相对于喷注平面的倾角，表示未经过修正的爆震波相对于喷注平面的倾角；
[0047]
(72)、根据多次迭代修正后的爆震波相对于喷注平面的倾角，确定激波坐标系的旋转角度，通过二维旋转矩阵，将出口方向转为水平，旋转角度以及二维旋转矩阵为如下形式：
[0048][0049][0050][0051]
其中，表示激波坐标系旋转角度，x
new
表示激波坐标系旋转后的流场中任意一点的x坐标，y
new
表示激波坐标系旋转后的流场中任意一点的y坐标，x
old
表示激波坐标系未旋转的流场中任意一点的x坐标，y
old
表示激波坐标系未旋转的流场中任意一点的y坐标，u
new
表示激波坐标系旋转后的流场中任意一点的速度x分量，v
new
表示激波坐标系旋转后的流场中任意一点的速度y分量，u
old
表示激波坐标系未旋转的流场中任意一点的速度x分量，v
old
表示激波坐标系未旋转的流场中任意一点的速度y分量；
[0052]
(73)、根据矢量加减原则，通过速度三角形法，将激波坐标系内的速度转换为实验室坐标系内的速度，转换关系式为如下形式：
[0053][0054]vlab
＝v
new
[0055]
其中，u
lab
表示实验室坐标系下任意一点速度的x分量，v
lab
表示实验室坐标系下任意一点速度的y分量。
[0056]
进一步的，所述步骤(6)中爆震波前可燃混气所在区间的上侧区间的气流角，在新一轮迭代的激波坐标系里修正为0。
[0057]
有益效果：由于采用了上述技术方案，本发明具有如下技术效果：
[0058]
本发明提供的一种基于空间推进算法的旋转爆震发动机流场参数分布预测方法，通过耦合空间推进算法和单一γ的cj爆震波模型，利用空间推进算法所对应的流线网格的随流特性，确定激波坐标系下爆震波后流动参数以及爆震波远下游流动参数，并通过坐标
转换获得实验室坐标系下的流动参数，可以快速解算旋转爆震发动机流场，快速评估旋转爆震发动机的气动性能，同时精确捕捉流场中的爆震波、剪切层以及诱导激波等流动结构，极大地缩短了旋转爆震发动机构型设计周期。
附图说明
[0059]
图1为本发明基于空间推进算法的旋转爆震发动机流场参数分布预测方法中的步骤(1)到步骤(6)确定的激波坐标系下的计算域；
[0060]
图2为旋转爆震发动机的实物结构示意图；
[0061]
图3为初值线上相邻两点的激波极曲线的相交过程；
[0062]
图4为空间推进算法从初值线到解点的步进过程，e-p表示格边点，c-p表示格心点，const x表示初值线以及解点所在直线，p-m表示普朗特-迈耶膨胀波，sh表示激波，sl表示无粘滑移线；
[0063]
图5为旋转爆震发动机流场中爆震波参数收敛历程；
[0064]
图6为激波坐标系下的旋转爆震发动机流线网格；
[0065]
图7为空间推进算法所获得的旋转爆震发动机温度场；
[0066]
图8为空间推进算法所获得的实验室坐标系下的旋转爆震发动机速度场；
[0067]
图中，11-上游爆震波对应初值线，12-爆震波远下游流场对应初值线，13-爆震波远下游流场冗余度段对应初值线，14-倾角与冗余段气流折转角一致的假想无粘壁面，15-倾角等于爆震波倾角的喷注平面阻塞段，16-可燃混合气模化进气段，17-爆震波远下游流场对应解点所在直线，18-下游爆震波对应解点所在位置，19-倾角与可燃混合气模化进气段的气流折转角一致的假想无粘壁面。
具体实施方式
[0068]
下面结合附图及实施例对本发明作更进一步的说明。
[0069]
本发明公开了一种基于空间推进算法的旋转爆震发动机流场参数分布预测方法，包括步骤：预先估计爆震波波前流动参数；根据单一比热比的cj爆震模型，确定爆震波后流动参数；将爆震波后流动参数和远下游流动参数赋值给初值线；通过空间推进算法确定激波坐标系下的流动参数；确定可燃混合气喷注起始位置和终了位置；迭代收敛旋转爆震发动机流场参数，并通过坐标转换获得实验室坐标系下的流动参数。
[0070]
结合图1所示，11为上游爆震波对应初值线，12为爆震波远下游流场对应初值线，13对应爆震波远下游流场冗余度段对应初值线。14对应倾角与冗余段气流折转角一致的假想无粘壁面，15对应倾角等于爆震波倾角的喷注平面阻塞段，16对应可燃混合气模化进气段，17对应爆震波远下游流场对应解点所在直线，18对应下游爆震波对应解点所在位置，19对应倾角与可燃混合气模化进气段的气流折转角一致的假想无粘壁面。
[0071]
激波坐标系是假设运动的爆震波线段ab保持不动，以线段ab方向作为y坐标方向，以与线段ab垂直的方向作为x坐标方向，形成的坐标系。实验室坐标系是以喷注平面的线段am方向组为x坐标方向，以与线段am垂直的方向作为y坐标方向，形成的坐标系。
[0072]
θ
ave
角为平均气流角，位于线段ij稍上游的位置，θ角为喷注平面(激波坐标系时倾斜的线段aj)和爆震波ab夹角的余角。
[0073]
预测方法具体包括以下步骤：
[0074]
步骤(1)，预先估计喷注面进口平均速度，根据可压缩流动关系式，确定爆震波波前流动参数；
[0075]
步骤(2)，将爆震波前流动参数代入单一比热比γ的cj爆震模型中，确定爆震波后流动参数；
[0076]
步骤(3)，第一次迭代中的爆震波远下游的流动参数由爆震波后压力膨胀至爆震波前压力获得，此后迭代中的爆震波远下游的流动参数则根据前一次计算结果，将爆震波后流动参数以及远下游流动参数赋值给初值线；如图1所示，初值线由线段ab和bd共同构成，而爆震波ab后的流动参数与前一次计算的线段ij爆震波后的流动参数一致，同时爆震波上方的线段bd的流动参数分布与前一次计算的线段ie的流动参数分布一致，具体数值则通过线性插值获得。
[0077]
步骤(4)，将初值线代入空间推进算法中，确定激波坐标系下的旋转爆震发动机流场结构及流动参数分布；
[0078]
步骤(5)，将喷注面解点的压力值与喷注总压相比较，当解点压力值小于喷注总压时，对喷注流场和爆震波下游流场进行耦合求解；若解点压力值大于等于总压时，喷注流场不存在，此时仅进行爆震波下游流场的求解。
[0079]
步骤(6)，确定爆震波远下游位置的爆震波前的流动参数的面积平均值，返回步骤(2)进行迭代，直至远下游的平均气流角度为0，跳出循环；
[0080]
步骤(7)，根据爆震波倾角大小，确定激波坐标系的旋转角度，并根据速度三角形确定实验室坐标系下的旋转爆震发动机流动参数分布；
[0081]
步骤(8)，将实验室坐标系下的旋转爆震发动机流动参数分布输出成数据文件。
[0082]
结合图3，本发明中，相邻两点的状态曲线的交点参数通过牛顿-拉弗森迭代方法获得。结合图4，本发明中，由初值线向下游解点推进过程中，通过直接计算激波和膨胀波之间的相互作用获得流动参数，不需额外构造单元点过程，通过自适应的改变格边的斜率，形成流线网格。
[0083]
本发明步骤(1)，具体包括如下步骤：
[0084]
(11)、根据旋转爆震发动机飞行工况，确定旋转爆震进口可燃混合气的总压pt和总温tt，并估计旋转爆震发动机进口平均速度
[0085]
(12)、根据可压缩流关系式，确定爆震波前流动参数，包括平均静温平均静压以及平均密度
[0086]
存在以下关系：
[0087][0088][0089][0090]
其中，所述c
p
和γ分别为可燃混合气和燃烧产物的定压比热和比热比，rg表示气体
常数。
[0091]
步骤(2)具体包括如下步骤：
[0092]
(21)、根据双γ的cj爆震模型，简化可得单一比热比γ的cj爆震模型：
[0093][0094][0095]
(22)、根据跨爆震波前后质量守恒定律、动量守恒定律以及能量守恒定律，确定爆震波后流动参数的表达式为如下形式：
[0096][0097][0098][0099]
其中，m
cj
表示爆震波马赫数，h表示无量纲释热量，q表示化学反应热，rg表示气体常数，表示爆震波前平均静温，表示爆震波前平均密度，表示爆震波前平均静压，表示爆震波后静压，表示爆震波后密度，表示爆震波后温度。
[0100]
步骤(4)具体包括如下步骤：
[0101]
(41)、通过主流方向无粘通量e确定中间变量χ，进而获得初值线上的状态量q，中间变量χ的表达式为如下形式：
[0102][0103]
e1、e2、e3、e4、是指主流方向无粘通量；
[0104]
(42)、根据初值线上速度分布，确定下游解点坐标，建立流线网格，下游解点的表达式为如下形式：
[0105][0106][0107]
表示x方向上第j个下游解点的x坐标，表示x方向上第j个初值线点的x坐标，δx表示网格沿x方向步进长度，表示y方向上第j个下游解点的y坐标，表示y方向上第j个初值线点的x坐标，表示y方向上第j初值线点的y向速度分量，表示y方向上第j初值线点的x方向速度分量；
[0108]
(43)、根据本质无振荡插值方法和限制函数，获得步进半步的沿x方向的和的数值解，进而获取步进长度为1/2位置的状态量和
[0109]
(44)、根据激波极曲线理论，通过空间x方向上步进长度为1/2设定值的相邻两点，确定空间推进算法所对应的步进长度为1/2位置的解点压力和气流角，进而获得黎曼自相似解，采用的相邻两点的状态曲线包括如下激波极曲线方程：
[0110][0111][0112]
其中，φb表示初值线下侧点的状态曲线，θb表示初值线下侧点的气流角，f(mb,αb)表示初值线下侧点的兰金许贡关系式，αb表示下游空间中任一点压力与初值线下侧点压力之比，ν(mb)表示初值线的下侧点的普朗特-迈耶函数，φ
t
表示初值线上侧点的状态曲线，θ
t
表示初值线上侧点的气流角，f(m
t
,α
t
)表示初值线上侧点的兰金许贡关系式，α
t
表示下游空间任一点压力与初值线上侧点压力之比，ν(m
t
)表示初值线的上侧点的普朗特-迈耶函数。
[0113]
(45)、通过对可压缩欧拉方程进行差分离散，获得单个空间递进步长下的流动参数，下游沿流动x方向的差分方程为如下形式：
[0114][0115]
其中，表示初值线相邻两点的y坐标差值，表示下游相邻两解点的y坐标差值，表示空间步进1/2长度位置的上侧格点的沿y方向的无粘通量，表示空间步进1/2长度位置的上侧格边的斜率，表示空间步进1/2长度位置的上侧格点的沿x方向上的无粘通量，表示空间步进1/2长度位置的下侧格点的沿y方向上的无粘通量，表示空间步进1/2长度位置的下侧格边的斜率，表示空间步进1/2长度位置的下侧格点的沿x方向上的无粘通量。
[0116]
步骤(5)中，将喷注面解点的压力值与喷注总压相比较，解点为流场区域中第n个constx线上的上游点通过步骤(44)和(45)获得的第n 1个constx线的下游点。其中constx代表图1中与爆震波ab平行的直线。
[0117]
所述步骤(6)包括如下步骤：
[0118]
步骤(61)、根据远下游沿y方向的流动参数分布，识别爆震波前可燃混合气所在区间；
[0119]
步骤(62)、确定爆震波前可燃混气所在区间的流动参数的面积平均值，面积平均公式为如下形式：
[0120]
[0121][0122][0123][0124]
其中，表示爆震波前平均速度的x分量，表示爆震波前平均速度的y分量，n
x
表示爆震波面法矢的x分量，ny表示爆震波面法矢的y分量，表示单位质量所含内能，a表示爆震波所占面积。根据牛顿-拉弗森多元迭代法，联立求解面积平均公式的四个方程，获得爆震波前流动参数以及
[0125]
迭代时，在步骤(3)赋值时用到的前一次的计算结果，是指根据步骤(62)得到的爆震波前流动参数进一步通过步骤(2)计算得到的爆震波后流动参数。
[0126]
所述步骤(7)包括如下步骤：
[0127]
(71)、根据爆震波前气流角，修正爆震波相对于喷注平面的倾角，修正公式为如下形式：
[0128][0129]
其中，表示经过修正的爆震波相对于喷注平面的倾角，表示未经过修正的爆震波相对于喷注平面的倾角；
[0130]
(72)、根据多次迭代修正后的爆震波相对于喷注平面的倾角，确定激波坐标系的旋转角度，通过二维旋转矩阵，将出口方向转为水平，旋转角度以及二维旋转矩阵为如下形式：
[0131][0132][0133][0134]
其中，表示激波坐标系旋转角度，x
new
表示激波坐标系旋转后的流场中任意一点的x坐标，y
new
表示激波坐标系旋转后的流场中任意一点的y坐标，x
old
表示激波坐标系未旋转的流场中任意一点的x坐标，y
old
表示激波坐标系未旋转的流场中任意一点的y坐标，u
new
表示激波坐标系旋转后的流场中任意一点的速度x分量，v
new
表示激波坐标系旋转后的流场中任意一点的速度y分量，u
old
表示激波坐标系未旋转的流场中任意一点的速度x分量，v
old
表示激波坐标系未旋转的流场中任意一点的速度y分量；
[0135]
(73)、根据矢量加减原则，通过速度三角形法，将激波坐标系内的速度转换为实验
室坐标系内的速度，转换关系式为如下形式：
[0136][0137]vlab
＝v
new
[0138]
其中，u
lab
表示实验室坐标系下任意一点速度的x分量，v
lab
表示实验室坐标系下任意一点速度的y分量。速度三角形如图1所示，表示激波坐标系下amlked区域里的流场速度，速度之间的换算通过矢量加减获得。由于激波坐标系中将激波的空间位置固定不动，而实际流动(实验室坐标系)中激波等流动结构沿着喷注平面切向高速运动，需要根据相对运动的变换规则，因此通过坐标转换获得实际旋转爆震发动机的流场。
[0139]
进一步的，所述步骤(6)中爆震波前可燃混气所在区间的上侧区间的气流角，在新一轮迭代的激波坐标系里修正为0。
[0140]
本发明的方法可用于快速评估旋转爆震发动机气动性能，即可以通过建立爆震过程和其下游流场结构的数学模型，避免传统时间推进算法耗费时间和计算资源的特点。
[0141]
为了更好的说明本发明，便于理解本发明的技术方案，本发明的典型但非限制性的实施例如下：
[0142]
旋转爆震发动机进口可燃混合气的喷注总压为5
×
105pa，喷注总温为300k，进口平均速度的初始值取450m/s，可燃预混气比热比为1.3961，可燃预混气气体常数为395.75j
·
kg-1
·
k-1
，燃烧产物比热比为1.1653，燃烧产物气体常数为346.20j
·
kg-1
·
k-1
，单位质量释热为5.4704
×
106j
·
kg-1
，可燃混合气活化温度15100k,单步化学反应指前因子1.0
×
109s-1
。旋转爆震发动机中环展开长度314mm，轴向长度75mm。
[0143]
图5是本发明快速解算迭代时爆震波各参数的收敛历程，θ
inc
表示爆震波倾角，hd表示爆震波长度，p1/p2表示爆震波前波后静压比，u
lab
表示爆震波速。图6是本发明快速解算时在激波坐标系下的流线网格。图7是本发明所得旋转爆震发动机的温度场，准确的捕捉到激波增压过程和爆震波后气流膨胀加速过程。图8是本发明所得旋转爆震发动机在实验室坐标系下的速度场，准确的捕捉到旋转爆震发动机的气流流向变化过程，实验室坐标系的定义是以喷注平面的线段am方向组为x坐标方向，以与线段am垂直的方向作为y坐标方向，形成的坐标系。
[0144]
表1是本发明快速解算所得的爆震波参数与传统时间推进算法解算所得爆震波参数的对比数据。
[0145]
表1
[0146]
[0147]
与传统时间推进算法所得旋转爆震发动机的爆震波参数相比，本发明所得爆震波高度相对偏差为1.38％，爆震波倾角相对偏差为1.35％，爆震波前后静压相对偏差为1.05％，爆震波波速相对偏差为0.97％，可见空间推进算法具有时间推进算法解算旋转爆震流场同等精度。
[0148]
表2是在4核8线程i7cpu环境下本发明解算旋转爆震发动机流场耗时与传统时间推进算法解算耗时的对比数据。
[0149]
表2
[0150]
算法cpu数目及型号耗时空间推进算法4核8线程，i75h传统时间推进算法4核8线程，i748h
[0151]
本发明是一种基于空间推进算法的旋转爆震发动机流场参数分布预测方法，使用4核8线程i7的cpu进行求解运算，总耗时5h。而传统的时间推进算法使用4核8线程进行求解运算，总耗时48h。本发明比传统时间推进算法求解速度快8.6倍，能有效的降低旋转爆震发动机构型设计和选取周期。
[0152]
以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。