(19)国家知识产权局 (12)发明 专利 (10)授权公告 号 (45)授权公告日 (21)申请 号 202111413401.6 (22)申请日 2021.11.25 (65)同一申请的已公布的文献号 申请公布号 CN 113987691 A (43)申请公布日 2022.01.28 (73)专利权人 中国人民解 放军国防科技大 学 地址 410073 湖南省长 沙市开福区德雅路 109号 (72)发明人 谢文佳　任伟杰　田正雨　张烨　 于航　杨帆　 (74)专利代理 机构 长沙国科天河知识产权代理 有限公司 432 25 专利代理师 段盼姣 (51)Int.Cl. G06F 30/15(2020.01)G06F 30/28(2020.01) G06F 113/08(2020.01) G06F 119/08(2020.01) G06F 119/14(2020.01) 审查员 沈晴 (54)发明名称 激波失稳的高精度混合计算方法、 装置、 设 备和存储介质 (57)摘要 本申请涉及一种上述激波失稳的高精度混 合计算方法、 装置、 设备和存储介质， 所述方法包 括： 对流场单元信息进行初始化， 并将当前时刻 的重构方式和数据格式设置为预设方式， 根据流 场信息以及当前时刻的重构方式， 对流场单元信 息进行重构， 根据重构结果， 采用当前时刻的数 据格式计算流场界面通量， 根据网格界面的通量 得到稳定性矩阵， 通过矩阵稳定性分析方法， 分 析稳定性矩阵的特征向量， 进而得到激波失稳 发 生的位置， 在失稳位置使用高鲁棒性的数值格式 进行计算， 其他位置使用低耗散数值格式计算。 通过在激波失稳位置施加针对性的处理， 实现了 稳定地捕捉激波。 该方法能够稳定捕捉激波， 进 而得到了高超声速飞行器表面热流的准确分布。 权利要求书4页 说明书10页 附图4页 CN 113987691 B 2022.08.23 CN 113987691 B 1.一种激波失稳的高精度混合计算方法， 其特 征在于， 所述方法包括： 对流场单 元信息进行初始化， 并将当前时刻的重构方式和数据格式设置为预设方式； 根据流场信息以及当前时刻的重构方式， 对流场单 元信息进行重构； 根据重构结果， 采用当前时刻的数据格式计算 流场界面 通量； 根据流场界面 通量， 对流场单 元信息进行 更新； 根据所述网格界面的通 量得到稳定性矩阵； 对稳定性矩阵的特 征向量进行分析， 确定 激波失稳位置； 在所述激波失稳位置和其他网格位置选择不同的重构方式和数值格式作为下一时刻 的重构方式和数值格式， 并更新当前时刻的重构方式和数据格式； 根据预设判断标准， 判断计算结果是否满足所述预设判断标准， 如果满足， 则计算完成 输出计算结果。 2.根据权利要求1所述的方法， 其特征在于， 根据所述网格界面的通量得到稳定性矩 阵， 包括： 根据所述网格界面的通量和矩阵稳定性分析的思路， 得到整个流场扰动随时间的传播 方程； 所述整个流场扰动随时间的传播方程 为： 其中， δQa为第a个网格状态量在当地的扰动， a＝1、 2、 ......、 M， M为网格界面中网格的 数量， S为稳定性矩阵， t为时间； 根据传播方程和预设公式进行计算， 得到稳定性矩阵； 所述预设公式为： 其中： Qm与Qk分别为第m个网格和第k个网格中的状态量在当地坐标系下的矢量， m和k为 大于等于1且小于等于M的整数， δQm与δQk分别为第m个网格和第k个网格的状态量在 当地的 扰动， Ωm为第m个网格的体积， 为第m个网格的边界， Amk为第m个网格与第k个网格边界 的面积， Fmk为第m个网格与第k个网格边界上的数值通量， Zm,Zk分别为第m个网格和第k个网 格中的状态量在全局坐标系下的矢量， Mmk为坐标转换矩阵。 3.根据权利要求1所述的方法， 其特 征在于， 激波失稳位置是指发生激波失稳的网格； 对稳定性矩阵的特 征向量进行分析， 确定 激波失稳位置， 包括： 求解稳定性矩阵的最大特征值， 并计算所述最大特征值对应的特征向量； 所述特征向 量包含四个守恒量对应的分量； 采用绝对值最大值的方式对所述 四个守恒量对应的分量进行归一化处理， 并将归一化 结果取算 术平均得到 激波失稳判断依据； 根据所述激波失稳判断依据和预设条件进行判断， 满足预设条件时， 则对应的网格发 生激波失稳。 4.根据权利要求1所述的方法， 其特征在于， 根据所述激波失稳位置， 确定网络中不同权　利　要　求　书 1/4 页 2 CN 113987691 B 2位置在下一时刻的重构方式和数值格式， 并更新当前时刻的重构方式和当前时刻的数据格 式， 包括： 在所述激波失稳位置在下一时刻选择MUSCL重构和HL L格式进行计算； 在所述激波失稳位置以外的其他网格在下一时刻选择5阶WENO格 式重构和全速域HLLC 格式计算； 根据不同位置在下一时刻选择选择的重构方式和数据格式对对应位置的当前时刻的 重构方式和数据格式进行 更新。 5.根据权利要求4所述的方法， 其特征在于， 在所述激波失稳位置在下一时刻选择 MUSCL重构和HL L格式进行计算， 步骤中， 所述激波失稳位置 MUSCL重构流场状态的表达式为： 其中： QL和QR分别是间断左侧和右侧的状态量， Qi‑1、 Qi、 Qi+1、 Qi+2分别为第i ‑1个网格、 第 i个网格、 第i+1个网格和第i+2个网格中的状态量在当地坐标系下的矢量， Δl为网格中心 格点到网格 边界的距离， ψL、 ψR为限制器 函数， 所述限制器 函数采用Van  Albada限制器， 所述 限制器函数的表达式为： 其中： T＝ L、 R， 所述激波失稳位置的HL L数值通量的表达式为： 其中， QL和QR分别为间断左侧和右侧的状态量； SL和SR指黎曼解的两个波速， SL为间断左 侧声波速度， SR为间断右侧声波速度； FR和FL分别为间断右侧和左侧的数值通量， 可以通过 公式计算： 其中， ρ 为密度， n＝(nx,ny)为法向矢量， p为压力， q＝unx+vny为逆变速度， u、 v分别为x方 向和y方向的速度， 下 标T代表间断左侧或右侧； H表示总焓。 6.根据权利要求4所述的方法， 其特征在于， 在所述激波失稳位置以外的其他网格在下 一时刻选择五阶WENO格式重构和全速域HL LC格式计算， 步骤中： 五阶WENO格式重构流场信息左值的公式为： QL＝ω0Q(0)+ω1Q(1)+ω2Q(2)权　利　要　求　书 2/4 页 3 CN 113987691 B 3