[0001] 本发明涉及发动机真空环境气固两相流的数值模拟技术领域，且具体地涉及一种CFD/DSMC耦合DEM的羽流月尘模拟技术。

[0002] 月尘是一直困扰月面着陆/上升探测器性能及航天员健康等突出问题的污染来源，已造成多次着陆月球的精密仪器和太阳能电池帆板的失效以及航天员的肺部损伤。月尘对探测器敏感仪器、热控系统、密封系统、机械结构和材料表面磨损等方面都具有很大的破坏性。通过对月尘在着陆器发动机高速羽流冲击下的动力学特性进行研究，模拟在探测器羽流高速冲击作用下月尘的卷起、扩散和下降过程，可以配合探测器找到有利的登月地点，有效降低月尘对探测器、月球车及宇航员的危害；进一步提出可靠有效的防尘和除尘方案，为载人登月和建立月球空间科研站重大战略任务中的人机长时驻月提供支撑，具有重要应用价值。

[0003] 月球表面无大气且低重力的环境使得发动机羽流扬尘特性研究非常复杂，国内外研究学者先后提出了地面实验、着陆图像分析和仿真研究三种方法。地面实验研究过程中无法长时间保持高速羽流的真空、低重力环境条件，无法在短时间内对不同任务的整个羽流扬尘过程进行定量实验。着陆图像虽然可以客观记录扬尘过程中的月面侵蚀特性，但有限的探月任务无法获得大量的数据，难以支撑羽流扬尘特性的全阶段研究。仿真研究不受限于实验条件和数据数量，可以通过图像数据分析结果，持续完善羽流扬尘模型，深入分析羽流扬尘特性，成为了此类研究的主要手段。

[0004] Zhang等研究了“嫦娥五号”低推力羽流表面侵蚀过程和侵蚀特性，建立了嫦娥五号着陆器发动机羽流仿真模型，采用Roberts侵蚀模型，使用拉格朗日框架下的离散相模型来研究月球尘埃粒子的动力学特性，但是该研究中流场的数值模拟只使用了CFD方法，没有考虑流场稀薄时N‑S方程失效的问题，也没有考虑月尘颗粒对羽流流场的反向影响。

[0005] Li等将模拟粒子运动的动力学方法与模拟稀薄羽流场的直接模拟蒙特卡罗方法相结合，建立了模拟羽流与月球尘埃相互作用的宏观气粒双向耦合方法，考虑了宏观气动力、羽流与月球尘埃之间的对流换热以及月尘粒子对羽流的反作用，并采用微观气粒双向耦合方法与宏观方法进行比较，但是该研究中没有建立月尘污染影响区模型，只考虑了流场中的月尘颗粒对羽流流场的影响而没有考虑月面侵蚀坑对羽流流场的影响。

[0006] Morris等采用DSMC方法对发动机喷管内的羽流特性及月面真空环境羽流扩散特性进行了研究，同样采用气固双向耦合的方法，但也没考虑到月面侵蚀坑对羽流流场的影响，并且没有对月尘颗粒的侵蚀建立物理模型。综上所述，国内外对于羽流冲击月尘运动的数值模拟，缺乏一种建立月尘颗粒与流场双向耦合模型并描述月尘详细侵蚀过程以及侵蚀坑对羽流反向影响的模拟方法，无法反映真实的羽流冲击月尘的物理过程。

[0007] 综上可知，现有技术中并没有能够同时模拟月尘与羽流相互作用、月面侵蚀坑与羽流相互作用，从而模拟出真实的羽流冲击月尘的物理过程的方法。

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

[0025] 基于N‑S方程的计算流体力学CFD(Compute Fluid Dynamics)是一种对流场进行数值模拟的方法。直接模拟蒙特卡罗方法DSMC(Direct Simulation Monte Carlo)是1969年由Bird G A提出的，通过用大量仿真分子模拟真实气体分子的运动与碰撞，并通过统计的方法获得流场各点的宏观流动参数，从而对流动进行模拟。离散元法DEM(Discrete Element Method)是由Cundall P.A.在1971年提出的，最早应用于岩土力学分析的一种基于分子动力学的分析方法

[0026] 本发明提供了一种CFD/DSMC耦合DEM的羽流月尘模拟技术应用方法，在二维轴对称的坐标体系下(也称为在月面着陆/上升器发动机类喷管燃气流动的轴对称坐标体系下)，将羽流冲击月尘的计算域划分为连续流区和稀薄流区，使用CFD‑DEM双向耦合模型和DSMC‑DEM双向耦合模型进行数值模拟，以连续流区的出口边界作为稀薄流区的入口边界来实现两个区域的耦合，模拟着陆器发动机羽流冲击月尘的物理过程，计算月尘在羽流冲击下的运动距离、滞空时间、数密度分布和速度分布，实现月尘在着陆器发动机羽流冲击下的动力学特性的仿真；

[0027] 其中，连续流区包括喷管及其附近的羽流核心区流场和流场下方的月尘污染区，分别使用CFD(computational fluid dynamics)或CFD/DSMC(direct simulation Monte Carlo)和DEM(discrete element method)方法进行数值模拟。稀薄流区包含距离喷管较远的羽流流场和相应流场月尘污染区，分别使用DSMC和DEM方法进行数值模拟，所述两个计算域之间由入口和出口边界条件连接，以连续流区的出口边界作为稀薄流区的入口边界来实现两个区域的耦合。

[0028] 基于DEM方法建立月尘粒子可计算模型，结合发动机羽流具体数据进行发动机接近月表工作状态下的月尘动力学仿真。假设月尘粒子按照一定的分布规律覆盖在月球表面，初速度为零，在每个时间步根据弹簧‑阻尼器模型计算月尘粒子之间的接触力，并计算羽流和电场对月尘粒子的作用力，更新粒子运动状态。

[0029] 在喷管出口附近的羽流致密区依然满足流体的连续性假设，通过求解N‑S方程对这部分羽流进行CFD模拟，避免粒子数密度过高时DSMC计算负担大的问题。使用DEM离散元方法模拟月尘颗粒，对月尘颗粒进行受力分析，包括流场的曳力、重力、电场力、颗粒之间的接触力、范德华力，从而根据牛顿第二定律更新月尘颗粒的运动状态，实现月尘颗粒侵蚀过程和空中运动轨迹的模拟。

[0030] 在距离喷管出口较远的羽流稀薄区不再满足流体的连续性假设，使用DSMC直接模拟蒙特卡洛方法进行流场的仿真，DSMC计算域的入口边界为CFD计算域的出口边界。DSMC中的模拟分子进行气体分子之间的概率碰撞和气体分子‑固体颗粒之间的碰撞，前者实现了稀薄流中对月尘颗粒的曳力效果，后者实现了月尘颗粒对气相的反向作用，从而实现了羽流低密度区的气‑固双向耦合。

[0031] 实施例：

[0032] 在二维轴对称的坐标体系下，将羽流冲击月尘的计算域划分为连续流区和稀薄流区。其中，连续流区包括喷管附近的核心区羽流流场和流场下方的月尘污染区，分别使用CFD和DEM方法进行数值模拟。稀薄流区包含距离喷管较远的羽流流场和相应流场下方的月尘污染区，分别使用DSMC和DEM方法进行数值模拟。所述两个计算域之间由入口和出口边界条件连接，以连续流区的出口边界作为稀薄流区的入口边界来实现两个区域的耦合。

[0033] 参照图1，为本发明的CFD/DSMC耦合DEM模拟方法的流程图。

[0034] 所述的连续流区的流体模型，使用CFD方法进行数值计算，控制方程为：

[0035]

[0036] 上式中，Q为守恒变量，F1、F2为无粘通量的两个分量，G1、G2为粘性通量的两个分量，t为时间，x为水平方向位置，y为竖直方向位置；

[0037] 其中，守恒变量Q为：

[0038]

[0039] 上式中，ρ为气体密度，u为x方向速度分量，v为y方向速度分量，et为总能量密度；

[0040] 无粘通量的两个分量F1、F2为：

[0041]

[0042] 上式中，ht为总焓密度；

[0043] 粘性通量的两个分量G1、G2为：

[0044]

[0045] 上式中，τxy为xy平面的剪切应力，qx为在x方向上的热流，qy为在y方向上的热流，τxx为x方向的正应力，τyy为y方向的正应力；

[0046] 其中，

[0047]

[0048]

[0049] 上式中，τij为粘性应力张量，p为压强，Cv为定容比热容，T为温度，μ为粘性系数，μt为湍流粘度， 为速度矢量，为梯度算子，Cp为定压比热容，Pr为普朗特数，i为，j为，i、j为剪切应力的指标。

[0050] 直接模拟蒙特卡罗方法(DSMC)通过用大量仿真分子模拟真实气体分子，对仿真分子记录其位置坐标、运动速度和内能，然后跟随着模拟分子的运动，更新其位置坐标、运动速度和内能，最后通过统计平均的方法获得流场各点的宏观参数，从而对流动进行模拟的方法。

[0051] 连续介质的宏观特性有：气体密度ρ，压强p和温度T，速度c0是矢量，其在X,Y,Z方向的分量分别为u0,v0,w0(需要说明的是，下标0代表宏观速度)。DSMC程序中采用了气体分子的微观属性。假设气体为混合气体，ns代表分子数密度，即单位体积内的s种类的分子数目。那么有：

[0052]

[0053] 上式中，ms为s种类气体的分子的质量， 为气体分子的平均质量，n为气体分子的平均数密度；

[0054] 宏观流动速度c0与分子速度c相关，故有：

[0055]

[0056] 上式中， 为s种类气体的分子的平均速度， 为气体分子的平均质量，为气体分子的平均速度；

[0057] 热运动速度c'的计算方法如下：

[0058] c'＝c‑c0

[0059] 在N‑S方程中，假设压力遵循帕斯卡准则，而且是一个标量。压强的标量可以看作是三个正应力法向值的平均值，即：

[0060]

[0061] 上式中，u'为热运动速度在X方向的分量，v'为热运动速度在Y方向的分量，w'为热运动速度在Z方向的分量， 为热运动速度平方的平均值(需要说明的是，上标'代表热运动速度)；

[0062] 在分子模型中，温度其实是分子所具有的平动能和内能的反映。体现平动能部分的温度和分子的热运动速度有关。能量均分原理指出：对每一个自由度来说,与其相关联的能量大小为kT/2，由于分子平动有三个方向的自由度，因此有：

[0063]

[0064] 上式中，k为玻尔兹曼常数，且k＝1.38×10‑23J/K，Ttr为平动温度；

[0065] 结合以上两式，故有：

[0066] p＝nkTtr

[0067] 由于在DSMC程序中存储的是每个分子的总速度，而非仅仅其热运动速度，因此可将上式写为：

[0068]

[0069] 至此，可以得到：与平动温度相关联的能量就等于分子总动能减去宏观流动的动能。

[0070] 分子内能是转动能、振动能和电子能之和，即：

[0071] εint＝εrot+εvlb+εel

[0072] 上式中，εint为内能，εrot为转动能，εvlb为振动能，εel为电子能；

[0073] 与内能相关联的温度是：

[0074]

[0075] 上式中，Tint为与内能相关联的温度， 为平均内能；

[0076] 其中，ξ为分子自由度，这样，温度即可以表示为：

[0077]

[0078] 离散元法(DEM)对离散的每一个颗粒计算其在每个时间步内的受力以及位移，根据计算的结果更新所有颗粒的位置。

[0079]

[0080] Fne、Fnv为法向的弹性力和粘滞力，nij为颗粒接触法向单位向量， 为法向接触力。

[0081] 线性接触模型中颗粒间的法向弹性力和粘滞力可采用以下方式计算：

[0082]

[0083]

[0084] 式中， 为两颗粒的法向相对速度； 为两颗粒的法向重叠长度，kn为等效刚度系数，cn为等效阻尼系数，且kn和cn可分别表示为：

[0085]

[0086] 式中， 为单元i、j的刚度系数； 为单元i、j的质量；为无量纲阻尼系数，ep为颗粒的回弹系数。

[0087] 本发明在探月着陆器发动机羽流冲击月尘的计算域内建立了时间推进的瞬态气固两相流模型，分为连续流区的CFD‑DEM双向耦合模型和稀薄流区的DSMC‑DEM双向耦合模型，实现了月尘颗粒堆积、月尘颗粒与流场相互作用的模拟，解决了传统模型难以描述月尘详细侵蚀过程以及侵蚀坑对羽流反向影响的难题，能够反映真实的羽流冲击月尘的物理过程，通过连续流区和稀薄流区进出口边界的耦合，实现了全计算域的统一时间推进，通过计算月尘在羽流冲击下的运动距离、滞空时间、数密度分布和速度分布，实现了月尘在着陆器发动机羽流冲击下的动力学特性的仿真。至此，本发明的CFD/DSMC耦合DEM的羽流月尘模拟技术过程构造完成。

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

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