一种适用于DSMC方法的高空羽流中相变模拟的实现方法

技术领域

本发明涉及飞行器空气动力试验技术领域，具体涉及适用于DSMC方法的高空羽流中相变模拟的实现方法。

背景技术

近年来，随着对临近空间战略意义认识的加深，各航空航天大国发展了众多长航时、大载荷、高超声速巡航的临近空间飞行器。这类飞行器在高空/在轨飞行时经常利用姿控/轨控发动机进行姿态调整、变轨等动作。由于飞行器所处环境压力低，燃气羽流会快速膨胀形成巨大的羽流流场。羽流的快速膨胀，流场中的温度、压力会急剧下降，导致其中H

在工程实践中，通常采用连续流区的计算流体力学方法模拟发动机内流，将发动机喷流出口宏观参数转换成DSMC方法微观模拟的输入参数，并特殊处理出口附近的高密度区，实现羽流污染/喷流干扰的DSMC方法模拟，具体来说，如《一种适用于DSMC方法的高空羽流中相变模拟的实现方法》就是这样的一种模拟方式：基于经典成核理论，采用DSMC方法计算得到的宏观参数，转化为适用于DSMC方法的仿真颗粒，对羽流流场中不同位置发生相变的颗粒参数进行计算，并得到在流场参数变化时颗粒参数的变化值。通过双向耦合算法，解耦计算每个仿真颗粒与周围气相的相互作用，得到固体颗粒温度和速度的变化值，以及气相流场参数的变化值。

发明内容

本发明的一个目的是解决至少上述问题和/或缺陷，并提供至少后面将说明的优点。

为了实现本发明的这些目的和其它优点，提供了一种适用于DSMC方法的高空羽流中相变模拟的实现方法，包括：

步骤一，对高空发动机喷管出口宏观参数、喷管出口燃气宏观参数进行转换，以得到DSMC计算所需的模拟输入参数；

步骤二，模拟分子从喷管出口进入流场，并模拟分子的运动和分子间的碰撞，通过DSMC计算得到的羽流流场中的当地宏观参数；

步骤三，基于经典成核理论，以及步骤二得到的当地宏观参数，转化为适用于DSMC方法的仿真颗粒数，对羽流流场中不同位置发生相变的颗粒参数进行计算，以得到在流场参数变化时颗粒参数的变化值Ⅰ；

步骤四，通过双向耦合算法，解耦计算步骤三中每个仿真颗粒与周围气相的相互作用，得到固体颗粒温度和速度的变化值Ⅱ，以及气相流场参数的变化值Ⅲ。

优选的是，在步骤三中，所述颗粒参数被配置为包括：颗粒核数、颗粒半径、颗粒质量、颗粒温度、颗粒速度。

优选的是，在步骤三中，流场参数变化时颗粒参数变化值Ⅰ的获取流程被配置为包括：

S31，确定每一个网格内产生的固体颗粒数，并基于以下公式获得对应的仿真颗粒数I：

I＝JΔtV/F

其中，J为凝结的成核率，V为网格体积，F

S32，在网格内的固体颗粒位置变换后，基于以下公式确定在新位置的气体参数条件下的颗粒参数变化值Ⅰ：

其中，λ

优选的是，在步骤四中，所述变化值Ⅱ是通过以下方式获得：

S410，设固体颗粒处于当地自由分子流的状态，且不考虑多原子气体的振动激发的情况下，在同一个网格里，每个DSMC气体仿真分子作用到一个固体颗粒上的力

其中：R

S411，则颗粒速度变化量能通过合力向量与因数Δt/M

固体颗粒温度能通过总的热传递率和Δt/(c

其中，c

颗粒的速度和温度通过下式得到：

T

其中，

在固体颗粒速度

优选的是，在步骤四中，所述变化值Ⅲ是通过以下方式获得：

S420，对每个固体颗粒，基于以下公式计算可能与其发生碰撞的仿真分子数n

其中，N

S422，对所有n

在整体坐标系下，

而碰撞后气体分子的绝对速度

其中，ε为方位角，δ为偏折角，镜面反射碰撞的偏转角分布为

f(δ)＝1/2sinδ，δ∈[0,π]

漫反射碰撞的偏转角分布为：

f(δ)＝0.02042δ

通过“取舍”法得到δ角。

本发明至少包括以下有益效果：

(1)本发明将固体颗粒转化为仿真颗粒，与DSMC仿真中的仿真分子一致，方便气相和固相的相间相互作用的计算。

(2)本发明将相间相互作用解耦计算，既考虑气相对固相的作用，又考虑固相对气相的作用，简化了相间相互作用的模拟过程。

(3)本发明将相变过程模拟为随机过程，与DSMC方法的随机过程一致，方便与DSMC方法的耦合

本发明的其它优点、目标和特征将部分通过下面的说明体现，部分还将通过对本发明的研究和实践而为本领域的技术人员所理解。

附图说明

图1为本发明一种适用于DSMC方法的高空羽流中相变模拟的实现方法的流程示意图。

图2为本发明中颗粒半径变化示意图。

具体实施方式

下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

如图1，本发明的一种适用于DSMC方法的高空羽流中相变模拟的实现方法，涉及高空/在轨飞行器羽流污染的数值模拟，具体针对多组分燃气等的羽流/喷流中容易发生相变的H

步骤一、基于连续流区的计算流体力学方法模拟发动机内流，其包括：

S11、确定发动机喷管燃烧室总温、总压、出口燃气组分，其过程为：假定发动机燃烧室总温为3150K、总压3.5MPa，出口燃气组分为N

S12、利用轴对称NS方程求解发动机喷管内流获得喷管出口燃气宏观参数，其过程为：在发动机型面喷管坐标系中，采用常规计算流体力学方法求解纳维-斯托克斯方程模拟喷管内部流动，获得喷管出口燃气宏观参数，即组分数密度、温度、速度；

步骤二、基于统计力学原理，将步骤一得到的发动机喷管出口宏观参数转化到DSMC方法模拟所需的分子微观量，将喷管出口燃气宏观参数根据麦克斯韦速度分布函数和接受-拒绝取舍法转换成当地模拟分子的数密度、速度、位置坐标等，作为DSMC方法模拟的输入参数，其包括：

S21、由于喷管出口流动向外，因此出口分子的速度分布函数为单边麦克斯韦速度分布函数，积分后得到穿过平面的分子数密度通量，其过程为：首先，由于DSMC方法的计算域边界由大量平面元组成。对于入流曲面，需要在每个模拟时间步指定模拟分子通过每个平面元的通量及其特性。考虑一个面积为A的平面单元，其单位法向量

与面元单位法向量之间的点积用于定义面元法向上的速度比(s

式中

其中n是数密度，erf(*)为误差函数。F

其中A为所考虑的面元的面积，Δt为DSMC时间步长，W

(2.08974617E+23 5.70999585E+23 3.06078424E+226.92001673E+211.03002865E+22)

S22、利用两个独立的随机数从完全麦克斯韦-玻尔兹曼分布中取样生成分子y和z方向的速度分量；从不完全的、有偏向的麦克斯韦-玻尔兹曼分布取样得到喷管出口外法向即x方向的速度分量，其过程为：通过面元的通量所需产生的分子数为：

N

注意，这里的

这样，在

因而有：

其中：

这里，

(1)计算值y＝-3+6R

(2)生成第二个随机数，0≤R

(3)计算f(y)/M，并检查

R

(a)如果上述表达式成立，则通过设置

(b)如果上述表达式不成立，则拒绝值y，返回到(1)，并生成另一个值y。

总之，对于N

C

C

C

由于转换前速度矢量垂直于喷管出口的y轴，z方向速度默认是0，可得C

S23、利用随机数在喷管出口计算网格单元上给出分子的位置坐标，其过程为:所有模拟分子的位置坐标通过随机数在面元上随机分布。

步骤三、DSMC计算，模拟分子从喷管出口进入流场，模拟碰撞分子的运动和碰撞过程。喷口附近高密度区域进行快速处理以适于DSMC方法模拟，其包括：

S31、确定每一碰撞网格和每一时间步长内每一模拟分子的碰撞次数，当平均碰撞次数达到十次时认为碰撞网格中的状态已达到平衡；

S32、在保证质量、动量以及能量守恒的前提下，采用平衡分布函数取样碰撞网格中的宏观参数。

在本实施例中，发动机喷管出口附近气体密度很高，若在高密度区使用DSMC方法，大量的时间花费在碰撞计算上。针对DSMC方法的这一缺陷，在发动机喷管出口附近采用平衡粒子模拟方法(EPSM)，即在模拟分子的平均碰撞次数大于十次以上的网格中，在保证质量守恒、动量守恒以及能量守恒的前提下，无需任何碰撞计算直接给出模拟分子的速度分布与能量分布。

步骤四、基于DSMC方法计算得到的羽流流场中的当地宏观参数，根据经典成核理论，计算羽流流场中不同位置发生相变的颗粒核数、颗粒半径、颗粒质量、颗粒温度、颗粒速度，并转化为适用于DSMC方法的仿真颗粒数，计算在流场参数发生变化后颗粒参数的变化，其包括：

S41、确定每一个网格内产生的固体颗粒数，给出相应的仿真颗粒数。确定每个仿真颗粒的半径、质量、温度、速度等参数；

在一个网格中，用过饱和度S判断当地是否发生相变，

其中：P为实际蒸气压；P

如果蒸气压力大于饱和蒸汽压，过饱和度S>1，蒸气可能出现凝结。凝结的成核率J指的是单位时间单位体积内所形成的达到临界半径的凝结核心的数量。

其中：N

ΔΦ

其中：r

转化为DSMC仿真凝结核的数量为：

I＝JΔtV/F

其中：V是网格体积；F

当过饱和度S的值增大到一定程度时，气态水分子会逐步凝聚在一起，形成含有若干水分子的核心，这些凝结核心是以r为半径的球状液滴。临界半径即为初步形成的凝结核心的半径，用r

/>

其中：ρ

确定网格内所有I个颗粒的参数：

颗粒半径均等于r

M

颗粒位置在网格内随机分布；

颗粒速度等于网格内气体宏观速度加上随机热运动速度，热运动速度与步骤2中的气体分子热运动速度一致。

颗粒温度T

T

其中：ΔT为相变产生的汽化潜热导致的颗粒温度增量。

S42、确定每个固体颗粒位置变化后，确定在新位置的气体参数条件下的颗粒参数。

颗粒在运动到一个新网格后(颗粒的运动将在步骤五介绍)，气相的参数发生变化，颗粒周围气体分子会吸附到颗粒表面，颗粒本身的分子会蒸发，导致颗粒也会发生变化(如图2，其中图2中A表示气体分子，箭头指向颗粒表示气体分子被吸附，即凝结，箭头指向外，表示颗粒表面分子的蒸发)。颗粒的尺寸变化为：

其中：λ

此处的ΔT为温度差Ts-Tv，Ts为压力P对应的饱和温度。

颗粒尺寸变化后，质量相应发生变化。

步骤五、采用双向耦合算法，解耦计算每个颗粒与周围气相的相互作用，得到固体颗粒温度和速度的变化，以及气相流场参数的变化，其包括：

S51、计算每个固体颗粒在运动过程中与同一网格内气体分子作用后，速度、温度的变化，确定新的速度、温度和位置。

考虑气相对固体颗粒的作用。假设固体颗粒处于当地自由分子流的状态，不考虑多原子气体的振动激发，在同一个网格里每个DSMC气体仿真分子作用到一个固体颗粒上的力和能量分别为：

其中：R

颗粒速度通过合力向量与因数Δt/M

固体颗粒温度通过总的热传递率和Δt/(c

其中，c

可得颗粒的速度和温度：

T

其中，

固体颗粒速度

S52、计算每个固体颗粒与可能的气体分子的碰撞，确定每个与之碰撞的气体分子碰撞后的参数。

考虑固体颗粒对周围气体的影响，首先要确定在每个时间步长内哪个仿真分子将与颗粒进行碰撞。对Bird的非时间计数方法进行修正，来确定与所选的颗粒可能发生碰撞的计算分子数n

其中，N

对所有n

整体坐标系下

碰撞后气体分子的绝对速度为

ε为方位角，在360°范围内随机给出；

δ为偏折角，镜面反射碰撞的偏转角分布为：

f(δ)＝1/2sinδ，δ∈[0,π]

漫反射碰撞的偏转角分布为：

f(δ)＝0.02042δ

可以通过“取舍”法得到δ角。

以上方案只是一种较佳实例的说明，但并不局限于此。在实施本发明时，可以根据使用者需求进行适当的替换和/或修改。

这里说明的设备数量和处理规模是用来简化本发明的说明的。对本发明的应用、修改和变化对本领域的技术人员来说是显而易见的。

尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用。它完全可以被适用于各种适合本发明的领域。对于熟悉本领域的人员而言，可容易地实现另外的修改。因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。