CN114996866A - 一种用于射流式脉冲发动机内部的数学模型和离散方法 - Google Patents

一种用于射流式脉冲发动机内部的数学模型和离散方法 Download PDF

Info

Publication number
CN114996866A
CN114996866A CN202210512362.3A CN202210512362A CN114996866A CN 114996866 A CN114996866 A CN 114996866A CN 202210512362 A CN202210512362 A CN 202210512362A CN 114996866 A CN114996866 A CN 114996866A
Authority
CN
China
Prior art keywords
flow
equation
grid
formula
phi
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202210512362.3A
Other languages
English (en)
Inventor
张国庆
王楠楠
徐勇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Institute of Technology BIT
Original Assignee
Beijing Institute of Technology BIT
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beijing Institute of Technology BIT filed Critical Beijing Institute of Technology BIT
Priority to CN202210512362.3A priority Critical patent/CN114996866A/zh
Publication of CN114996866A publication Critical patent/CN114996866A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/08—Thermal analysis or thermal optimisation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computing Systems (AREA)
  • Fluid Mechanics (AREA)
  • Algebra (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Aerodynamic Tests, Hydrodynamic Tests, Wind Tunnels, And Water Tanks (AREA)

Abstract

本发明提供一种用于射流式脉冲发动机内部的数学模型和离散方法，包括如下步骤：步骤一：建立数学模型，以RANS方程为控制方程，步骤二：对步骤一建立的数学模型及流体计算域进行网格划分，对步骤一中得到的水下固体火箭发动机喷管几何模型建立流体计算域；步骤三：控制方程的离散和求解；步骤四：基于步骤一至步骤三开展水下固体火箭发动机喷管摆动高速射流流动特性的数值计算分析，该设计解决了用单域网格生成方法无法解决或效果不佳的复杂构型的流场计算问题；分区网格生成以后，不同区域允许采用不同的网格***或用不同的数值计算格式求解。

Description

一种用于射流式脉冲发动机内部的数学模型和离散方法
技术领域
本发明是一种用于射流式脉冲发动机内部的数学模型和离散方法，属于射流式脉冲发动机技术领域。
背景技术
超音速射流元件结构复杂，难以生成高质量的单域结构计算网格，而内部流场又是复杂的粘性流场，采用非结构网格对边界层的计算精度影响较大，现有技术中采用非结构网格，但计算结果不太理想，不能很好地反映控制口附近的边界层分离现象和激波与边界层的相互干扰现象，现在急需一种用于射流式脉冲发动机内部的数学模型和离散方法来解决上述出现的问题。
发明内容
针对现有技术存在的不足，本发明目的是提供一种用于射流式脉冲发动机内部的数学模型和离散方法，以解决上述背景技术中提出的问题，本发明结构合理，解决用单域网格生成方法无法解决或效果不佳的复杂构型的流场计算问题。
为了实现上述目的，本发明是通过如下的技术方案来实现：一种用于射流式脉冲发动机内部的数学模型和离散方法，包括如下步骤：
步骤一：建立数学模型，以RANS方程为控制方程，结合Realizable k-epsilon湍流模式、近壁区处理、封闭性方程、边界条件和初始条件，作为超音速射流元件数值计算的数学模型；
步骤二：对步骤一建立的数学模型及流体计算域进行网格划分，对步骤一中得到的水下固体火箭发动机喷管几何模型建立流体计算域，并对两流体计算域交界面进行几何优化，在计算域之间进行数据传递，实现对喷管摆动工况的模拟，分块对流域进行网格划分，并进行网格无关性验证，得到最优网格数；
步骤三：控制方程的离散和求解；
步骤四：基于步骤一至步骤三开展水下固体火箭发动机喷管摆动高速射流流动特性的数值计算分析。
进一步地，在步骤一中，建立数学模型包括控制方程、湍流模式、近壁区的湍流模式以及改进壁面函数，控制方程根据超音速射流元件内部的流动特点，忽略体积力和热源，并假定壁面为绝热壁面，在以下各式中，以ui取代其时均值
Figure BDA0003639843150000021
u'i为脉动值，在笛卡儿坐标系中有：
(1)连续方程
Figure BDA0003639843150000022
(2)动量方程
Figure BDA0003639843150000023
式中，雷诺应力
Figure BDA0003639843150000024
将在湍流模型中给出其表达式，τij是粘性应力张量，可表示为：
Figure BDA0003639843150000025
流体运动的应变率—
Figure BDA0003639843150000026
质点的体积膨胀率—
Figure BDA0003639843150000027
δij是Kronecker符号；μ'是“容积粘性系数”或“第二粘性系数”，它反映由体积变化引起流体偏离热力学压强的粘性应力，除高温和高频声波等极端情况外，一般的气体运动中可近似认为μ'＝0；
(3)能量方程
Figure BDA0003639843150000031
式中：
Jj'是混合气体组元j'的扩散量
Figure BDA0003639843150000032
其中，
Figure BDA0003639843150000033
cp,j'是混合气体组元j'的定压比热，取Tref＝298.15K；
Figure BDA0003639843150000034
λ是导热系数，cp是定压比热，Prt＝0.85
Figure BDA0003639843150000035
有效湍流涡旋粘性系数μeff和湍流涡旋粘性系数μt将在湍流模式中给出其表达式；
湍流模式中在超音速射流元件的流动边界层存在强的逆压梯度和流动分离，并且在一定条件下还有较大回流出现，采用Realizable k-epsilon模式来计算湍流，湍动能k和湍流耗散率ε的定义如下：
Figure BDA0003639843150000036
(1)湍动能k的输运方程
Figure BDA0003639843150000037
式中：
Gk＝μtS2，平均应变率张量系数
Figure BDA0003639843150000038
YM＝2ρεMt 2，湍流马赫数
Figure BDA0003639843150000039
a是当地声速
(2)湍流耗散率ε的输运方程
Figure BDA0003639843150000041
式中：
Figure BDA0003639843150000042
且η＝Sk/ε
(3)湍流涡旋粘性系数μt
在高雷诺数时，湍流涡旋粘性系数μt表示为
Figure BDA00036398431500000413
式中：
Figure BDA0003639843150000043
而
Figure BDA0003639843150000044
Figure BDA0003639843150000045
与旋转参考面有关，在本计算中，不予考略；
A0＝4.04，
Figure BDA0003639843150000046
Figure BDA0003639843150000047
在低雷诺数时，有效湍流涡旋粘性系数μeff可由下式求得
Figure BDA0003639843150000048
式中：
Figure BDA0003639843150000049
Cυ≈100
(4)雷诺应力
Figure BDA00036398431500000410
由Boussinesq假设，雷诺应力
Figure BDA00036398431500000411
和时均速度梯度的关系为
Figure BDA00036398431500000412
(5)模型常数
C2＝1.9，σk＝1.0，σε＝1.2；
在近壁区的湍流模式中采用一种将二层壁面模型与改进壁面函数有机结合的壁面处理方法：在二层近壁模型中，按湍流雷诺数的大小将整个流动区域分成粘性影响区和完全湍流区；湍流雷诺数Rey的定义如下：
Figure BDA0003639843150000051
y是计算网格单元中心到壁面的法向距离；
在完全湍流区
Figure BDA0003639843150000052
用Realizable k-epsilon湍流模式；
在受粘性影响的近壁区
Figure BDA0003639843150000053
用Wolfstein一方程模型[89]，在该模型中，湍动能k的输运方程与Realizable k-epsilon湍流模式的相同，湍流涡旋粘性系数μt,2layer和湍流耗散率εt,2layer的计算如下：
Figure BDA0003639843150000054
εt,2layer＝k1.5/lε
上式中：
μt是高雷诺数时的湍流涡旋粘性系数，其计算见式
lμ和lε是长度尺度，计算如下：
lμ＝cly[1-exp(-Rey/Aμ)]
lε＝cly[1-exp(-Rey/Aε)]
上式中的常数
Figure BDA0003639843150000055
Aμ＝70，Aε＝2cl
为了使两层模型定义的μt,2layer能够与外层高雷诺数定义的μt平稳过渡，定义了一个改进的湍流涡旋粘性系数μt,enh，其定义如下
μt,enh＝λεμt+(1-λε)μt,2layer
式中λε是折衷函数，其作用是：当外层的k—ε解与两层模型计算的解不匹配时，防止解的收敛过程受阻，它的定义如下：
Figure BDA0003639843150000056
上式中，常数A确定折衷函数λε的宽度；若给定变量ΔRey，要使λε的值不超过其远场值的1％，A的值应为
Figure BDA0003639843150000061
而ΔRey的值可在
Figure BDA0003639843150000062
值的5％～20％之间给定；
改进壁面函数为一种混合函数，该函数将线性壁面律和对数壁面律有机结合在一起；对于平均速度，其改进壁面函数可表示为
Figure BDA0003639843150000063
上式中，
Figure BDA0003639843150000064
其中a＝0.01c，b＝5/c，而
Figure BDA0003639843150000065
对式求导得
Figure BDA0003639843150000066
改进湍流壁面律可以表示为
Figure BDA0003639843150000067
其中
Figure BDA0003639843150000068
Figure BDA0003639843150000069
为对数律斜率保持恒定的位置，在计算中取值为60；
Figure BDA00036398431500000610
Figure BDA00036398431500000611
Figure BDA00036398431500000612
从上述α、β和γ的表示式可知，改进后的湍流壁面律考虑了压强梯度和热的影响；
改进层流壁面律可以表示为
Figure BDA0003639843150000071
改进温度壁面函数表示为
Figure BDA0003639843150000072
其中，
Figure BDA0003639843150000073
Pr为分子普朗特数；
封闭性方程
为了使上述方程封闭，须加如下封闭性方程：
状态方程：f(p,ρ,T)＝0
本文利用完全气体状态方程计算流场的密度，即：ρ＝p/RT；
定压比热：
Figure BDA0003639843150000074
R是通用气体常数，M是气体的分子量，f是分子运动的自由度数，对多原子分子，取f＝7；
导热系数：
Figure BDA0003639843150000075
分子粘性系数μ，当温度不太高也不特别低时采用Sutherland公式计算，即
Figure BDA0003639843150000076
μ0是参考温度T0时的粘性系数，S为Sutherland常数；
对空气(210K～1900K)，μ0＝1.716×10-5Pa﹒s，T0＝273K,S＝111K；
当温度较低时，采用Keyes公式计算，即
Figure BDA0003639843150000077
超音速射流元件由三个入口、两个出口和壁面组成，在三个入口中，有一个是主气源入口，另外两个为控制流入口，根据超音速射流元件的流动特点和FLUENT软件中对边界条件的定义方法，将主气源入口边界条件定义为压力入口边界条件，两控制流入口边界条件定义为质量流入边界条件，两输出口边界条件定义为压力出口边界条件，对于定常或非定常流动，当采用时间推进算法求解时，除了应给出边界条件外，还需要给出初始条件，以下为计算超音速射流元件内部的定常或非定常流动时各边界条件和初始条件的提法及相应的数学处理；
在压力入口边界条件给出总压(P0)、总温(T0)、流动方向、湍流强度和水力直径；此外，还给出静压(Ps)，用来初始化流场，由理想气体的等熵关系式可得
Figure BDA0003639843150000081
上式中
Figure BDA0003639843150000082
其中，c是音速，k＝cp/cv是等熵指数，
密度可由理想气体的状态方程求得，即
Figure BDA0003639843150000083
在式上式中的静温Ts可由下式求出
Figure BDA0003639843150000084
可得到压力入口处的密度、静温和速度标量值，再根据给定的流动方向得到压力入口处在各方向上的速度分量值；
质量流入边界条件
给出质量流率
Figure BDA0003639843150000085
总温(T0)、静压(Ps)、流动方向、湍流强度和
水力直径；则质量流量按下式计算
Figure BDA0003639843150000091
上式中，vn为法向速度，A为入口面积
入口处的密度由理想气体状态方程求得，即
Figure BDA0003639843150000092
上式中，如果入口处的流动为超音速流动，Ps采用给定的静压值，入口处的流动为亚音速流动，则Ps由入口处计算网格的格心值外推得到；
式中的静温可由下式求出
Figure BDA0003639843150000093
在质量流入边界条件中，通过法向速度vn和变量的入口值，可得到入口处相应变量的通量，如：质量的通量为ρvn，湍动能的通量为ρkvn。
进一步地，网格分为超音速射流元件网格分区和边界层网格，整个计算区域分成10个子域，每个子域网格单独生成，其中，子域1采用非结构网格，子域8采用结构网格和非结构网格共存的混合网格，余下的8个子域采用结构网格，在对各子域生成计算网格时，尽量使各子域连接处的网格完全连续，不需要进行插值计算，以便减少区域边界的信息交换量，使边界的处理变得简单，这样，在射流元件内部流场的跨区域激波捕捉过程中，因为区域边界的变量传递守恒，就可以保证激波在跨区域处连续，根据超音速射流元件内各壁面附近的流动特点，通过给定适当的h1、m和R值，就可生成疏密程度与流场参数的变化梯度大体一致的边界层网格。
进一步地，在步骤三中，控制方程的离散和求解，超音速射流元件内部的流动为跨音速粘性内流，且存在扭曲区域，因此，采用有限体积法在物理空间曲线网格上将控制方程离散为能进行数值计算的代数方程；
有限体积法是用积分形式的控制方程进行描述的，故对任意二维控制体S，将控制方程改成积分形式的矢量方程为：
Figure BDA0003639843150000101
式中，
Figure BDA0003639843150000102
和
Figure BDA0003639843150000103
分别定义为：
Figure BDA0003639843150000104
其中，l是控制体S的环线长度，
Figure BDA0003639843150000105
是速度矢量，u和v是笛卡尔坐标系中的速度分量，E是单位质量的总能，τ是粘性应力张量，q是热流量；
空间离散及线性化
计算域所生成的网格中，对任一网格单元S，某一标量φ的稳态守恒输运方程可表示为：
Figure BDA0003639843150000106
式中
Figure BDA0003639843150000107
—环线长度矢量；
Γφ—φ的耗散系数；
Figure BDA0003639843150000108
—φ的梯度；
Sφ—单位面积φ的源项；
应用有限体积法对方程进行离散，则该方程在网格单元上的离散形式为：
Figure BDA0003639843150000109
式中
Nlines—包围网格单元S的边数；
φf—通过边f输运的φ值；
Figure BDA0003639843150000111
—通过边f的质量流量；
Figure BDA0003639843150000112
—边f的长度；
Figure BDA0003639843150000113
在边f法向上的分量；
S—网格单元的面积；
通过对网格单元中心值进行Taylor级数展开，可得到各个网格单元交界线上的流动参数φl的二阶迎风格式表达式，即：
Figure BDA0003639843150000114
式中
φ—网格单元的中心值；
Figure BDA0003639843150000115
—上游网格单元质心到边f中心的位移矢量；
Figure BDA0003639843150000116
—φ的梯度，定义如下：
Figure BDA0003639843150000117
Figure BDA0003639843150000118
是与边f相邻的两网格单元值的平均值，耗散项采用二阶精度的中心差分格式离散；
用有限体积法对湍动能k和湍动能耗散率ε的输运方程进行离散，离散方程是含有未知变量φ及其相邻网格单元未知量的非线性方程，须进行线性化，其线性化形式为：
Figure BDA0003639843150000119
下标nb代表相邻的网格单元，ap和anb分别为φ和φnb的线性化系数；
非定常方程空间方向的离散和定常方程的空间离散相同，时间方向的离散如下：
若某一标量φ随时间变化的历程用微分形式表示为
Figure BDA00036398431500001110
式中，F表示空间方向的离散项，则一阶精度的离散形式为
Figure BDA0003639843150000121
二阶精度的离散形式为
Figure BDA0003639843150000122
时间导数离散后，对F(φ)可以采取不同时间层上的φ值进行估算，根据估算F(φ)所采用φ值的时间层不同，时间方向的离散方法分隐式和显式两种算法；
(1)隐格式
在隐格式中，采用未知时间层上的φ值估算F(φ)，改写为
Figure BDA0003639843150000123
则一阶精度的隐格式为
φn+1＝φn+ΔtF(φn+1)
非线性方程，须先围绕φn+1进行线性化，然后采用迭代法求解，其迭代格式为
φi＝φn+ΔtF(φi)
类似地，可得二阶精度的隐格式为
Figure BDA0003639843150000124
迭代格式为
Figure BDA0003639843150000125
(2)显格式
在显格式中，采用未知时间层上的φ值估算F(φ)，改写为
Figure BDA0003639843150000126
第n+1时间层的φn+1值可以根据已存在的第n时间层的φn计算，即
φn+1＝φn+ΔtF(φn)。
进一步地，附壁流动工作状态的推力验证和切换过程中的切换时间验证是通过单通道冷气测试***完成的，该测试***由高压大流量调压***、数据采集与处理***、射流元件试验台动架和控制台四部分组成，高压大流量调压***包含高压压气机和高压气瓶组，主要为射流元件提供高压主气源和控制气流；试验台动架用于固定射流元件，在射流元件试验台动架上安装着两个高频响高精度的力传感器，用来检测射流元件在附壁流动工作状态时的合推力大小和切换过程中合推力的变化，数据采集与处理***用于采集力传感器传来的电压信号，并对其进行去噪处理；截止阀的通断、控制阀的转换、数据的采集与处理由控制台控制。
本发明的有益效果：本发明的一种用于射流式脉冲发动机内部的数学模型和离散方法，该设计解决了用单域网格生成方法无法解决或效果不佳的复杂构型的流场计算问题；分区网格生成以后，不同区域允许采用不同的网格***或用不同的数值计算格式求解；分区网格生成后的求解过程中，任何时刻，只有一个子域的数据存于计算机缓冲区，因此放宽了对计算机存储器的要求；采用分区网格计算，适合于多处理机上的并行计算，具有良好的守恒性，适合于内流和带激波的流场计算；能方便地描述复杂的计算域，可以采用笛卡儿坐标系处理扭曲区域，同时保持精度较高的近似解，第二类边界条件和第一类条件一样易于处理，适应于跨音速无粘流动和粘性流动计算。
附图说明
通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：
图1为本发明一种用于射流式脉冲发动机内部的数学模型和离散方法中射流元件测试***示意图；
图2为本发明一种用于射流式脉冲发动机内部的数学模型和离散方法中Ps为80ba推力变化波形图；
图3为本发明一种用于射流式脉冲发动机内部的数学模型和离散方法中Ps为100ba推力变化波形图；
图中：1-空气压缩机、2-过滤器、3-气瓶、4-减压阀、5-压力表、6-力传感器、7-超音速射流元件、8-三位三通电磁方向控制阀、9-可调节流阀、10-截止阀。
具体实施方式
为使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施方式，进一步阐述本发明。
请参阅图1，本发明提供一种技术方案：一种用于射流式脉冲发动机内部的数学模型和离散方法，一种用于射流式脉冲发动机内部的数学模型和离散方法，包括如下步骤：
步骤一：建立数学模型，以RANS方程为控制方程，结合Realizable k-epsilon湍流模式、近壁区处理、封闭性方程、边界条件和初始条件，作为超音速射流元件数值计算的数学模型；
步骤二：对步骤一建立的数学模型及流体计算域进行网格划分，对步骤一中得到的水下固体火箭发动机喷管几何模型建立流体计算域，并对两流体计算域交界面进行几何优化，在计算域之间进行数据传递，实现对喷管摆动工况的模拟，分块对流域进行网格划分，并进行网格无关性验证，得到最优网格数；
步骤三：控制方程的离散和求解；
步骤四：基于步骤一至步骤三开展水下固体火箭发动机喷管摆动高速射流流动特性的数值计算分析。
进一步地，在步骤一中，建立数学模型包括控制方程、湍流模式、近壁区的湍流模式以及改进壁面函数，控制方程根据超音速射流元件内部的流动特点，忽略体积力和热源，并假定壁面为绝热壁面，在以下各式中，以ui取代其时均值
Figure BDA0003639843150000151
u'i为脉动值，在笛卡儿坐标系中有：
(1)连续方程
Figure BDA0003639843150000152
(2)动量方程
Figure BDA0003639843150000153
式中，雷诺应力
Figure BDA0003639843150000154
将在湍流模型中给出其表达式，τij是粘性应力张量，可表示为：
Figure BDA0003639843150000155
流体运动的应变率—
Figure BDA0003639843150000156
质点的体积膨胀率—
Figure BDA0003639843150000157
δij是Kronecker符号；μ'是“容积粘性系数”或“第二粘性系数”，它反映由体积变化引起流体偏离热力学压强的粘性应力，除高温和高频声波等极端情况外，一般的气体运动中可近似认为μ'＝0；
(3)能量方程
Figure BDA0003639843150000158
式中：
Jj'是混合气体组元j'的扩散量
Figure BDA0003639843150000159
其中，
Figure BDA00036398431500001510
cp,j'是混合气体组元j'的定压比热，取Tref＝298.15K；
Figure BDA00036398431500001511
λ是导热系数，cp是定压比热，Prt＝0.85
Figure BDA0003639843150000161
有效湍流涡旋粘性系数μeff和湍流涡旋粘性系数μt将在湍流模式中给出其表达式；
湍流模式中在超音速射流元件的流动边界层存在强的逆压梯度和流动分离，并且在一定条件下还有较大回流出现，采用Realizable k-epsilon模式来计算湍流，湍动能k和湍流耗散率ε的定义如下：
Figure BDA0003639843150000162
(1)湍动能k的输运方程
Figure BDA0003639843150000163
式中：
Gk＝μtS2，平均应变率张量系数
Figure BDA0003639843150000164
YM＝2ρεMt 2，湍流马赫数
Figure BDA0003639843150000165
a是当地声速
(2)湍流耗散率ε的输运方程
Figure BDA0003639843150000166
式中：
Figure BDA0003639843150000167
且η＝Sk/ε
(3)湍流涡旋粘性系数μt
在高雷诺数时，湍流涡旋粘性系数μt表示为
Figure BDA0003639843150000168
式中：
Figure BDA0003639843150000171
而
Figure BDA0003639843150000172
Figure BDA0003639843150000173
与旋转参考面有关，在本计算中，不予考略；
A0＝4.04，
Figure BDA0003639843150000174
Figure BDA0003639843150000175
在低雷诺数时，有效湍流涡旋粘性系数μeff可由下式求得
Figure BDA0003639843150000176
式中：
Figure BDA0003639843150000177
Cυ≈100
(4)雷诺应力
Figure BDA0003639843150000178
由Boussinesq假设，雷诺应力
Figure BDA0003639843150000179
和时均速度梯度的关系为
Figure BDA00036398431500001710
(5)模型常数
C2＝1.9，σk＝1.0，σε＝1.2；
在近壁区的湍流模式中采用一种将二层壁面模型与改进壁面函数有机结合的壁面处理方法：在二层近壁模型中，按湍流雷诺数的大小将整个流动区域分成粘性影响区和完全湍流区；湍流雷诺数Rey的定义如下：
Figure BDA00036398431500001711
y是计算网格单元中心到壁面的法向距离；
在完全湍流区
Figure BDA00036398431500001712
用Realizable k-epsilon湍流模式；
在受粘性影响的近壁区
Figure BDA00036398431500001713
用Wolfstein一方程模型[89]，在该模型中，湍动能k的输运方程与Realizable k-epsilon湍流模式的相同，湍流涡旋粘性系数μt,2layer和湍流耗散率εt,2layer的计算如下：
Figure BDA0003639843150000181
εt,2layer＝k1.5/lε
上式中：
μt是高雷诺数时的湍流涡旋粘性系数，其计算见式
lμ和lε是长度尺度，计算如下：
lμ＝cly[1-exp(-Rey/Aμ)]
lε＝cly[1-exp(-Rey/Aε)]
上式中的常数
Figure BDA0003639843150000182
Aμ＝70，Aε＝2cl
为了使两层模型定义的μt,2layer能够与外层高雷诺数定义的μt平稳过渡，定义了一个改进的湍流涡旋粘性系数μt,enh，其定义如下
μt,enh＝λεμt+(1-λε)μt,2layer
式中λε是折衷函数，其作用是：当外层的k—ε解与两层模型计算的解不匹配时，防止解的收敛过程受阻，它的定义如下：
Figure BDA0003639843150000183
上式中，常数A确定折衷函数λε的宽度；若给定变量ΔRey，要使λε的值不超过其远场值的1％，A的值应为
Figure BDA0003639843150000184
而ΔRey的值可在
Figure BDA0003639843150000185
值的5％～20％之间给定；
改进壁面函数为一种混合函数，该函数将线性壁面律和对数壁面律有机结合在一起；对于平均速度，其改进壁面函数可表示为
Figure BDA0003639843150000186
上式中，
Figure BDA0003639843150000187
其中a＝0.01c，b＝5/c，而
Figure BDA0003639843150000188
对式求导得
Figure BDA0003639843150000191
改进湍流壁面律可以表示为
Figure BDA0003639843150000192
其中
Figure BDA0003639843150000193
Figure BDA0003639843150000194
为对数律斜率保持恒定的位置，在计算中取值为60；
Figure BDA0003639843150000195
Figure BDA0003639843150000196
Figure BDA0003639843150000197
从上述α、β和γ的表示式可知，改进后的湍流壁面律考虑了压强梯度和热的影响；
改进层流壁面律可以表示为
Figure BDA0003639843150000198
改进温度壁面函数表示为
Figure BDA0003639843150000199
其中，
Figure BDA00036398431500001910
Pr为分子普朗特数；
封闭性方程
为了使上述方程封闭，须加如下封闭性方程：
状态方程：f(p,ρ,T)＝0
本文利用完全气体状态方程计算流场的密度，即：ρ＝p/RT；
定压比热：
Figure BDA0003639843150000201
R是通用气体常数，M是气体的分子量，f是分子运动的自由度数，对多原子分子，取f＝7；
导热系数：
Figure BDA0003639843150000202
分子粘性系数μ，当温度不太高也不特别低时采用Sutherland公式计算，即
Figure BDA0003639843150000203
μ0是参考温度T0时的粘性系数，S为Sutherland常数；
对空气(210K～1900K)，μ0＝1.716×10-5Pa﹒s，T0＝273K,S＝111K；
当温度较低时，采用Keyes公式计算，即
Figure BDA0003639843150000204
超音速射流元件由三个入口、两个出口和壁面组成，在三个入口中，有一个是主气源入口，另外两个为控制流入口，根据超音速射流元件的流动特点和FLUENT软件中对边界条件的定义方法，将主气源入口边界条件定义为压力入口边界条件，两控制流入口边界条件定义为质量流入边界条件，两输出口边界条件定义为压力出口边界条件，对于定常或非定常流动，当采用时间推进算法求解时，除了应给出边界条件外，还需要给出初始条件，以下为计算超音速射流元件内部的定常或非定常流动时各边界条件和初始条件的提法及相应的数学处理；
在压力入口边界条件给出总压(P0)、总温(T0)、流动方向、湍流强度和水力直径；此外，还给出静压(Ps)，用来初始化流场，由理想气体的等熵关系式可得
Figure BDA0003639843150000211
上式中
Figure BDA0003639843150000212
其中，c是音速，k＝cp/cv是等熵指数，
密度可由理想气体的状态方程求得，即
Figure BDA0003639843150000213
在式上式中的静温Ts可由下式求出
Figure BDA0003639843150000214
可得到压力入口处的密度、静温和速度标量值，再根据给定的流动方向得到压力入口处在各方向上的速度分量值；
质量流入边界条件
给出质量流率
Figure BDA0003639843150000218
总温(T0)、静压(Ps)、流动方向、湍流强度和
水力直径；则质量流量按下式计算
Figure BDA0003639843150000215
上式中，vn为法向速度，A为入口面积
入口处的密度由理想气体状态方程求得，即
Figure BDA0003639843150000216
上式中，如果入口处的流动为超音速流动，Ps采用给定的静压值，入口处的流动为亚音速流动，则Ps由入口处计算网格的格心值外推得到；
式中的静温可由下式求出
Figure BDA0003639843150000217
在质量流入边界条件中，通过法向速度vn和变量的入口值，可得到入口处相应变量的通量，如：质量的通量为ρvn，湍动能的通量为ρkvn。
进一步地，网格分为超音速射流元件网格分区和边界层网格，整个计算区域分成10个子域，每个子域网格单独生成，其中，子域1采用非结构网格，子域8采用结构网格和非结构网格共存的混合网格，余下的8个子域采用结构网格，在对各子域生成计算网格时，尽量使各子域连接处的网格完全连续，不需要进行插值计算，以便减少区域边界的信息交换量，使边界的处理变得简单，这样，在射流元件内部流场的跨区域激波捕捉过程中，因为区域边界的变量传递守恒，就可以保证激波在跨区域处连续，根据超音速射流元件内各壁面附近的流动特点，通过给定适当的h1、m和R值，就可生成疏密程度与流场参数的变化梯度大体一致的边界层网格。
进一步地，在步骤三中，控制方程的离散和求解，超音速射流元件内部的流动为跨音速粘性内流，且存在扭曲区域，因此，采用有限体积法在物理空间曲线网格上将控制方程离散为能进行数值计算的代数方程；
有限体积法是用积分形式的控制方程进行描述的，故对任意二维控制体S，将控制方程改成积分形式的矢量方程为：
Figure BDA0003639843150000221
式中，
Figure BDA0003639843150000222
和
Figure BDA0003639843150000223
分别定义为：
Figure BDA0003639843150000224
其中，l是控制体S的环线长度，
Figure BDA0003639843150000225
是速度矢量，u和v是笛卡尔坐标系中的速度分量，E是单位质量的总能，τ是粘性应力张量，q是热流量；
空间离散及线性化
计算域所生成的网格中，对任一网格单元S，某一标量φ的稳态守恒输运方程可表示为：
Figure BDA0003639843150000231
式中
Figure BDA0003639843150000232
—环线长度矢量；
Γφ—φ的耗散系数；
Figure BDA0003639843150000233
—φ的梯度；
Sφ—单位面积φ的源项；
应用有限体积法对方程进行离散，则该方程在网格单元上的离散形式为：
Figure BDA0003639843150000234
式中
Nlines—包围网格单元S的边数；
φf—通过边f输运的φ值；
Figure BDA0003639843150000235
—通过边f的质量流量；
Figure BDA0003639843150000236
—边f的长度；
Figure BDA0003639843150000237
在边f法向上的分量；
S—网格单元的面积；
通过对网格单元中心值进行Taylor级数展开，可得到各个网格单元交界线上的流动参数φl的二阶迎风格式表达式，即：
Figure BDA0003639843150000238
式中
φ—网格单元的中心值；
Figure BDA0003639843150000239
—上游网格单元质心到边f中心的位移矢量；
Figure BDA0003639843150000241
—φ的梯度，定义如下：
Figure BDA0003639843150000242
Figure BDA0003639843150000243
是与边f相邻的两网格单元值的平均值，耗散项采用二阶精度的中心差分格式离散；
用有限体积法对湍动能k和湍动能耗散率ε的输运方程进行离散，离散方程是含有未知变量φ及其相邻网格单元未知量的非线性方程，须进行线性化，其线性化形式为：
Figure BDA0003639843150000244
下标nb代表相邻的网格单元，ap和anb分别为φ和φnb的线性化系数；
非定常方程空间方向的离散和定常方程的空间离散相同，时间方向的离散如下：
若某一标量φ随时间变化的历程用微分形式表示为
Figure BDA0003639843150000245
式中，F表示空间方向的离散项，则一阶精度的离散形式为
Figure BDA0003639843150000246
二阶精度的离散形式为
Figure BDA0003639843150000247
时间导数离散后，对F(φ)可以采取不同时间层上的φ值进行估算，根据估算F(φ)所采用φ值的时间层不同，时间方向的离散方法分隐式和显式两种算法；
(1)隐格式
在隐格式中，采用未知时间层上的φ值估算F(φ)，改写为
Figure BDA0003639843150000248
则一阶精度的隐格式为
φn+1＝φn+ΔtF(φn+1)
非线性方程，须先围绕φn+1进行线性化，然后采用迭代法求解，其迭代格式为
φi＝φn+ΔtF(φi)
类似地，可得二阶精度的隐格式为
Figure BDA0003639843150000251
迭代格式为
Figure BDA0003639843150000252
(2)显格式
在显格式中，采用未知时间层上的φ值估算F(φ)，改写为
Figure BDA0003639843150000253
第n+1时间层的φn+1值可以根据已存在的第n时间层的φn计算，即
φn+1＝φn+ΔtF(φn)。
附壁流动工作状态的推力验证和切换过程中的切换时间验证是通过单通道冷气测试***完成的，该测试***由高压大流量调压***、数据采集与处理***、射流元件试验台动架和控制台四部分组成，高压大流量调压***包含高压压气机和高压气瓶组，主要为射流元件提供高压主气源和控制气流。试验台动架用于固定射流元件，在射流元件试验台动架上安装着两个高频响高精度的力传感器，用来检测射流元件在附壁流动工作状态时的合推力大小和切换过程中合推力的变化，数据采集与处理***用于采集力传感器传来的电压信号，并对其进行去噪处理。截止阀的通断、控制阀的转换、数据的采集与处理由控制台控制。
附壁流动工作状态中输出推力的大小是衡量超音速射流元件性能好坏的重要指标之一，因此，本课题首先计算出多种主气源压强产生的推力值，然后在同等条件下通过实验获得了各主气源压强对应的推力值。表1列出了不同主气源压强Ps对应的推力值，其中，test为射流元件在X方向上产生的合推力的测试值。fl、fr分别表示射流元件两输出口在X方向上的计算推力值，F为fl和fr的合力。在计算Ps＝80bar、90bar和100bar时，fl、fr和F呈周期性波动，表中的计算值是时均值。Error表示计算值F与测试值Test的相对误差。
表1为推力的计算值与实验值
Ps(bar) Test(N) f<sub>r</sub>(N) f<sub>l</sub>(N) F(N) Error(％)
20 44.0 49.65 0.42 49.23 11.89
30 79.9 84.34 1.95 82.40 3.13
40 106.0 108.91 4.43 104.48 -1.43
45 126.3 142.09 4.58 137.51 8.87
50 142.6 164.59 5.66 158.93 11.45
55 161.4 170.27 5.46 164.81 2.11
60 179.3 186.87 6.43 180.44 0.64
70 202.8 220.68 8.55 212.12 4.6
80 239.9 263.05 18.04 245.01 2.1
90 275.4 296.36 21.78 274.58 -0.3
100 313.2 346.18 27.89 318.29 1.63
从表1中可看出，计算值与实验值吻合较好，尤其是当主气源压强超过50bar后，相对误差的绝对值小于5％。
图2为相同条件下通过实验获得的推力随时间变化的波形图，其中，上方的方波为输入的控制信号，下方的曲线为两个力传感器测得的合推力输出，曲线为消除噪音并光滑后的合推力输出，图2Ps＝80bar，图3Ps＝100bar，图2中，推力上升时间是1.82ms，推力下降时间是1.28ms；图3中，推力上升时间是1.18ms，推力下降时间是0.85ms。
由计算结果与实验结果的对比可知，计算获得的推力上升时间和推力下降时间与实验测得的数据比较接近
以上显示和描述了本发明的基本原理和主要特征和本发明的优点，对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。
此外，应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。

Claims (5)

1.一种用于射流式脉冲发动机内部的数学模型和离散方法，其特征在于：包括如下步骤：
步骤一：建立数学模型，以RANS方程为控制方程，结合Realizable k-epsilon湍流模式、近壁区处理、封闭性方程、边界条件和初始条件，作为超音速射流元件数值计算的数学模型；
步骤二：对步骤一建立的数学模型及流体计算域进行网格划分，对步骤一中得到的水下固体火箭发动机喷管几何模型建立流体计算域，并对两流体计算域交界面进行几何优化，在计算域之间进行数据传递，实现对喷管摆动工况的模拟，分块对流域进行网格划分，并进行网格无关性验证，得到最优网格数；
步骤三：控制方程的离散和求解；
步骤四：基于步骤一至步骤三开展水下固体火箭发动机喷管摆动高速射流流动特性的数值计算分析。
2.根据权利要求1所述的一种用于射流式脉冲发动机内部的数学模型和离散方法，其特征在于：在步骤一中，建立数学模型包括控制方程、湍流模式、近壁区的湍流模式以及改进壁面函数，控制方程根据超音速射流元件内部的流动特点，忽略体积力和热源，并假定壁面为绝热壁面，在以下各式中，以ui取代其时均值
Figure FDA0003639843140000011
u'i为脉动值，在笛卡儿坐标系中有：
(1)连续方程
Figure FDA0003639843140000012
(2)动量方程
Figure FDA0003639843140000013
式中，雷诺应力
Figure FDA0003639843140000021
将在湍流模型中给出其表达式，τij是粘性应力张量，可表示为：
Figure FDA0003639843140000022
流体运动的应变率—
Figure FDA0003639843140000023
质点的体积膨胀率—
Figure FDA0003639843140000024
δij是Kronecker符号；μ'是“容积粘性系数”或“第二粘性系数”，它反映由体积变化引起流体偏离热力学压强的粘性应力，除高温和高频声波等极端情况外，一般的气体运动中可近似认为μ'＝0；
(3)能量方程
Figure FDA0003639843140000025
式中：
Jj'是混合气体组元j'的扩散量
Figure FDA0003639843140000026
其中，
Figure FDA0003639843140000027
cp,j'是混合气体组元j'的定压比热，取Tref＝298.15K；
Figure FDA0003639843140000028
λ是导热系数，cp是定压比热，Prt＝0.85
Figure FDA0003639843140000029
有效湍流涡旋粘性系数μeff和湍流涡旋粘性系数μt将在湍流模式中给出其表达式；
湍流模式中在超音速射流元件的流动边界层存在强的逆压梯度和流动分离，并且在一定条件下还有较大回流出现，采用Realizable k-epsilon模式来计算湍流，湍动能k和湍流耗散率ε的定义如下：
Figure FDA0003639843140000031
(1)湍动能k的输运方程
Figure FDA0003639843140000032
式中：
Gk＝μtS2，平均应变率张量系数
Figure FDA0003639843140000033
YM＝2ρεMt 2，湍流马赫数
Figure FDA0003639843140000034
a是当地声速
(2)湍流耗散率ε的输运方程
Figure FDA0003639843140000035
式中：
Figure FDA0003639843140000036
且η＝Sk/ε
(3)湍流涡旋粘性系数μt
在高雷诺数时，湍流涡旋粘性系数μt表示为
Figure FDA0003639843140000037
式中：
Figure FDA0003639843140000038
而
Figure FDA0003639843140000039
Figure FDA00036398431400000310
与旋转参考面有关，在本计算中，不予考略；
A0＝4.04，
Figure FDA00036398431400000311
Figure FDA00036398431400000312
在低雷诺数时，有效湍流涡旋粘性系数μeff可由下式求得
Figure FDA0003639843140000041
式中：
Figure FDA0003639843140000042
Cυ≈100
(4)雷诺应力
Figure FDA0003639843140000043
由Boussinesq假设，雷诺应力
Figure FDA0003639843140000044
和时均速度梯度的关系为
Figure FDA0003639843140000045
(5)模型常数
C2＝1.9，σk＝1.0，σε＝1.2；
在近壁区的湍流模式中采用一种将二层壁面模型与改进壁面函数有机结合的壁面处理方法：在二层近壁模型中，按湍流雷诺数的大小将整个流动区域分成粘性影响区和完全湍流区；湍流雷诺数Rey的定义如下：
Figure FDA0003639843140000046
y是计算网格单元中心到壁面的法向距离；
在完全湍流区
Figure FDA0003639843140000047
用Realizable k-epsilon湍流模式；
在受粘性影响的近壁区
Figure FDA0003639843140000048
用Wolfstein一方程模型[89]，在该模型中，湍动能k的输运方程与Realizable k-epsilon湍流模式的相同，湍流涡旋粘性系数μt,2layer和湍流耗散率εt,2layer的计算如下：
Figure FDA0003639843140000049
εt,2layer＝k1.5/lε
上式中：
μt是高雷诺数时的湍流涡旋粘性系数，其计算见式
lμ和lε是长度尺度，计算如下：
lμ＝cly[1-exp(-Rey/Aμ)]
lε＝cly[1-exp(-Rey/Aε)]
上式中的常数
Figure FDA0003639843140000051
Aμ＝70，Aε＝2cl
为了使两层模型定义的μt,2layer能够与外层高雷诺数定义的μt平稳过渡，定义了一个改进的湍流涡旋粘性系数μt,enh，其定义如下
μt,enh＝λεμt+(1-λε)μt,2layer
式中λε是折衷函数，其作用是：当外层的k—ε解与两层模型计算的解不匹配时，防止解的收敛过程受阻，它的定义如下：
Figure FDA0003639843140000052
上式中，常数A确定折衷函数λε的宽度；若给定变量ΔRey，要使λε的值不超过其远场值的1％，A的值应为
Figure FDA0003639843140000053
而ΔRey的值可在
Figure FDA0003639843140000054
值的5％～20％之间给定；
改进壁面函数为一种混合函数，该函数将线性壁面律和对数壁面律有机结合在一起；对于平均速度，其改进壁面函数可表示为
Figure FDA0003639843140000055
上式中，
Figure FDA0003639843140000056
其中a＝0.01c，b＝5/c，而
Figure FDA0003639843140000057
对式求导得
Figure FDA0003639843140000058
改进湍流壁面律可以表示为
Figure FDA0003639843140000059
其中
Figure FDA0003639843140000061
Figure FDA0003639843140000062
为对数律斜率保持恒定的位置，在计算中取值为60；
Figure FDA0003639843140000063
Figure FDA0003639843140000064
Figure FDA0003639843140000065
从上述α、β和γ的表示式可知，改进后的湍流壁面律考虑了压强梯度和热的影响；
改进层流壁面律可以表示为
Figure FDA0003639843140000066
改进温度壁面函数表示为
Figure FDA0003639843140000067
其中，
Figure FDA0003639843140000068
Pr为分子普朗特数；
封闭性方程
为了使上述方程封闭，须加如下封闭性方程：
状态方程：f(p,ρ,T)＝0
本文利用完全气体状态方程计算流场的密度，即：ρ＝p/RT；
定压比热：
Figure FDA0003639843140000069
R是通用气体常数，M是气体的分子量，f是分子运动的自由度数，对多原子分子，取f＝7；
导热系数：
Figure FDA0003639843140000071
分子粘性系数μ，当温度不太高也不特别低时采用Sutherland公式计算，即
Figure FDA0003639843140000072
μ0是参考温度T0时的粘性系数，S为Sutherland常数；
对空气(210K～1900K)，μ0＝1.716×10-5Pa﹒s，T0＝273K,S＝111K；
当温度较低时，采用Keyes公式计算，即
Figure FDA0003639843140000073
超音速射流元件由三个入口、两个出口和壁面组成，在三个入口中，有一个是主气源入口，另外两个为控制流入口，根据超音速射流元件的流动特点和FLUENT软件中对边界条件的定义方法，将主气源入口边界条件定义为压力入口边界条件，两控制流入口边界条件定义为质量流入边界条件，两输出口边界条件定义为压力出口边界条件，对于定常或非定常流动，当采用时间推进算法求解时，除了应给出边界条件外，还需要给出初始条件，以下为计算超音速射流元件内部的定常或非定常流动时各边界条件和初始条件的提法及相应的数学处理；
在压力入口边界条件给出总压(P0)、总温(T0)、流动方向、湍流强度和水力直径；此外，还给出静压(Ps)，用来初始化流场，由理想气体的等熵关系式可得
Figure FDA0003639843140000074
上式中
Figure FDA0003639843140000081
其中，c是音速，k＝cp/cv是等熵指数，
密度可由理想气体的状态方程求得，即
Figure FDA0003639843140000082
在式上式中的静温Ts可由下式求出
Figure FDA0003639843140000083
可得到压力入口处的密度、静温和速度标量值，再根据给定的流动方向得到压力入口处在各方向上的速度分量值；
质量流入边界条件
给出质量流率
Figure FDA0003639843140000084
总温(T0)、静压(Ps)、流动方向、湍流强度和水力直径；则质量流量按下式计算
Figure FDA0003639843140000085
上式中，vn为法向速度，A为入口面积
入口处的密度由理想气体状态方程求得，即
Figure FDA0003639843140000086
上式中，如果入口处的流动为超音速流动，Ps采用给定的静压值，入口处的流动为亚音速流动，则Ps由入口处计算网格的格心值外推得到；
式中的静温可由下式求出
Figure FDA0003639843140000087
在质量流入边界条件中，通过法向速度vn和变量的入口值，可得到入口处相应变量的通量，如：质量的通量为ρvn，湍动能的通量为ρkvn。
3.根据权利要求1所述的一种用于射流式脉冲发动机内部的数学模型和离散方法，其特征在于：网格分为超音速射流元件网格分区和边界层网格，整个计算区域分成10个子域，每个子域网格单独生成，其中，子域1采用非结构网格，子域8采用结构网格和非结构网格共存的混合网格，余下的8个子域采用结构网格，在对各子域生成计算网格时，尽量使各子域连接处的网格完全连续，不需要进行插值计算，以便减少区域边界的信息交换量，使边界的处理变得简单，这样，在射流元件内部流场的跨区域激波捕捉过程中，因为区域边界的变量传递守恒，就可以保证激波在跨区域处连续，根据超音速射流元件内各壁面附近的流动特点，通过给定适当的h1、m和R值，就可生成疏密程度与流场参数的变化梯度大体一致的边界层网格。
4.根据权利要求1所述的一种用于射流式脉冲发动机内部的数学模型和离散方法，其特征在于：在步骤三中，控制方程的离散和求解，超音速射流元件内部的流动为跨音速粘性内流，且存在扭曲区域，因此，采用有限体积法在物理空间曲线网格上将控制方程离散为能进行数值计算的代数方程；
有限体积法是用积分形式的控制方程进行描述的，故对任意二维控制体S，将控制方程改成积分形式的矢量方程为：
Figure FDA0003639843140000091
式中，
Figure FDA0003639843140000092
和
Figure FDA0003639843140000093
分别定义为：
Figure FDA0003639843140000094
其中，l是控制体S的环线长度，
Figure FDA0003639843140000095
是速度矢量，u和v是笛卡尔坐标系中的速度分量，E是单位质量的总能，τ是粘性应力张量，q是热流量；
空间离散及线性化
计算域所生成的网格中，对任一网格单元S，某一标量φ的稳态守恒输运方程可表示为：
Figure FDA0003639843140000101
式中
Figure FDA0003639843140000102
—环线长度矢量；
Γφ—φ的耗散系数；
▽φ—φ的梯度；
Sφ—单位面积φ的源项；
应用有限体积法对方程进行离散，则该方程在网格单元上的离散形式为：
Figure FDA0003639843140000103
式中
Nlines—包围网格单元S的边数；
φf—通过边f输运的φ值；
Figure FDA0003639843140000104
—通过边f的质量流量；
Figure FDA0003639843140000105
—边f的长度；
Figure FDA0003639843140000106
—
Figure FDA0003639843140000107
在边f法向上的分量；
S—网格单元的面积；
通过对网格单元中心值进行Taylor级数展开，可得到各个网格单元交界线上的流动参数φl的二阶迎风格式表达式，即：
Figure FDA0003639843140000108
式中
φ—网格单元的中心值；
Figure FDA0003639843140000109
—上游网格单元质心到边f中心的位移矢量；
Figure FDA0003639843140000111
—φ的梯度，定义如下：
Figure FDA0003639843140000112
Figure FDA0003639843140000113
是与边f相邻的两网格单元值的平均值，耗散项采用二阶精度的中心差分格式离散；
用有限体积法对湍动能k和湍动能耗散率ε的输运方程进行离散，离散方程是含有未知变量φ及其相邻网格单元未知量的非线性方程，须进行线性化，其线性化形式为：
Figure FDA0003639843140000114
下标nb代表相邻的网格单元，ap和anb分别为φ和φnb的线性化系数；
非定常方程空间方向的离散和定常方程的空间离散相同，时间方向的离散如下：
若某一标量φ随时间变化的历程用微分形式表示为
Figure FDA0003639843140000115
式中，F表示空间方向的离散项，则一阶精度的离散形式为
Figure FDA0003639843140000116
二阶精度的离散形式为
Figure FDA0003639843140000117
时间导数离散后，对F(φ)可以采取不同时间层上的φ值进行估算，根据估算F(φ)所采用φ值的时间层不同，时间方向的离散方法分隐式和显式两种算法；
(1)隐格式
在隐格式中，采用未知时间层上的φ值估算F(φ)，改写为
Figure FDA0003639843140000118
则一阶精度的隐格式为
φn+1＝φn+ΔtF(φn+1)
非线性方程，须先围绕φn+1进行线性化，然后采用迭代法求解，其迭代格式为
φi＝φn+ΔtF(φi)
类似地，可得二阶精度的隐格式为
Figure FDA0003639843140000121
迭代格式为
Figure FDA0003639843140000122
(2)显格式
在显格式中，采用未知时间层上的φ值估算F(φ)，改写为
Figure FDA0003639843140000123
第n+1时间层的φn+1值可以根据已存在的第n时间层的φn计算，即
φn+1＝φn+ΔtF(φn)。
5.根据权利要求1所述的一种用于射流式脉冲发动机内部的数学模型和离散方法，其特征在于：附壁流动工作状态的推力验证和切换过程中的切换时间验证是通过单通道冷气测试***完成的，该测试***由高压大流量调压***、数据采集与处理***、射流元件试验台动架和控制台四部分组成，高压大流量调压***包含高压压气机和高压气瓶组，主要为射流元件提供高压主气源和控制气流；试验台动架用于固定射流元件，在射流元件试验台动架上安装着两个高频响高精度的力传感器，用来检测射流元件在附壁流动工作状态时的合推力大小和切换过程中合推力的变化，数据采集与处理***用于采集力传感器传来的电压信号，并对其进行去噪处理；截止阀的通断、控制阀的转换、数据的采集与处理由控制台控制。
CN202210512362.3A 2022-05-12 2022-05-12 一种用于射流式脉冲发动机内部的数学模型和离散方法 Pending CN114996866A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210512362.3A CN114996866A (zh) 2022-05-12 2022-05-12 一种用于射流式脉冲发动机内部的数学模型和离散方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210512362.3A CN114996866A (zh) 2022-05-12 2022-05-12 一种用于射流式脉冲发动机内部的数学模型和离散方法

Publications (1)

Publication Number Publication Date
CN114996866A true CN114996866A (zh) 2022-09-02

Family

ID=83027785

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210512362.3A Pending CN114996866A (zh) 2022-05-12 2022-05-12 一种用于射流式脉冲发动机内部的数学模型和离散方法

Country Status (1)

Country Link
CN (1) CN114996866A (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116522740A (zh) * 2023-06-30 2023-08-01 中国空气动力研究与发展中心计算空气动力研究所 一种发动机喷嘴的油气界面捕捉方法、装置、设备及介质
CN116579259A (zh) * 2023-04-24 2023-08-11 中国人民解放军陆军装甲兵学院 弹道三维瞬态流场建模及多物理场数值计算方法及装置
CN116579271A (zh) * 2023-07-13 2023-08-11 中国空气动力研究与发展中心计算空气动力研究所 基于笛卡尔网格下对壁面函数的修正方法及装置

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116579259A (zh) * 2023-04-24 2023-08-11 中国人民解放军陆军装甲兵学院 弹道三维瞬态流场建模及多物理场数值计算方法及装置
CN116579259B (zh) * 2023-04-24 2024-02-09 中国人民解放军陆军装甲兵学院 弹道三维瞬态流场建模及多物理场数值计算方法及装置
CN116522740A (zh) * 2023-06-30 2023-08-01 中国空气动力研究与发展中心计算空气动力研究所 一种发动机喷嘴的油气界面捕捉方法、装置、设备及介质
CN116522740B (zh) * 2023-06-30 2023-09-05 中国空气动力研究与发展中心计算空气动力研究所 一种发动机喷嘴的油气界面捕捉方法、装置、设备及介质
CN116579271A (zh) * 2023-07-13 2023-08-11 中国空气动力研究与发展中心计算空气动力研究所 基于笛卡尔网格下对壁面函数的修正方法及装置
CN116579271B (zh) * 2023-07-13 2024-03-29 中国空气动力研究与发展中心计算空气动力研究所 基于笛卡尔网格下对壁面函数的修正方法及装置

Similar Documents

Publication Publication Date Title
CN114996866A (zh) 一种用于射流式脉冲发动机内部的数学模型和离散方法
Alkishriwi et al. A large-eddy simulation method for low Mach number flows using preconditioning and multigrid
Deiwert Numerical simulation of high Reynolds number transonic flows
Arunajatesan et al. Hybrid RANS-LES modeling for cavity aeroacoutics predictions
Meister Comparison of different Krylov subspace methods embedded in an implicit finite volume scheme for the computation of viscous and inviscid flow fields on unstructured grids
Duraisamy et al. Adjoint based techniques for uncertainty quantification in turbulent flows with combustion
Haider et al. Mathematical analysis of flow passing through a rectangular nozzle
Chen et al. A new robust carbuncle-free roe scheme for strong shock
Brehm et al. Towards a viscous wall model for immersed boundary methods
Chen et al. Cartesian grid method for gas kinetic scheme on irregular geometries
Kawai et al. A dynamic wall model for large-eddy simulation of high Reynolds number compressible flows
Hinman Laminar near wake of hypersonic blunt bodies
Subrahmanyam et al. A universal velocity profile for near-wall flows
Zhan et al. Linear discrete velocity model-based lattice Boltzmann flux solver for simulating acoustic propagation in fluids
Gopalakrishnan et al. Computation of high-subsonic and transonic flows by a lattice Boltzmann method
Isaev et al. Simulation of a turbulent supersonic underexpanded jet flowing into a submerged space with the help of a shear stress transfer model
Hamdi et al. Link-Wise Artificial Compressibility Method for attached and separated flows past airfoils
Hiejima Development of linear unstable modes in supersonic streamwise vortices using a weighted compact nonlinear scheme
Volkov Preconditioning of the Euler and Navier-Stokes equations in low-velocity flow simulation on unstructured grids
Zhao et al. Application of Gas-Kinetic Scheme for Continuum and Near-Continuum Flow on Unstructured Mesh
Suzen et al. Application of several turbulence models for high speed shear layer flows
Raverdy et al. Large-Eddy Simulation of the Flow Around a Low Pressure Turbine Blade
Jemcov et al. A High-Resolution Primitive Variable Solver for Compressible Flow Simulation
Titarev et al. Rarefied gas flow into vacuum through a long circular pipe composed of two sections of different radii
Hsu et al. Numerical study of vortex-dominated flows for wings at high incidence and sideslip

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination