一种可压缩流体的数值模拟方法

摘要

本发明适用于计算流体力学领域，提供了一种可压缩流体的数值模拟方法，包括：S10：读取计算网格、时刻的物理变量分布、边界条件、试验时长；S20：在计算网格上，根据时刻的物理变量获得守恒变量和对流通量、、；S30：分别计算对流通量、、在物理空间上的数值通量、、；S40：将守恒变量和数值通量、、代入离散欧拉方程得到n+1时刻的守恒变量及其对应的物理变量；S50：n从1~k遍历，其中k为时刻试验时长的序号，将守恒变量及时刻的物理变量作为最终的试验结果。通过本发明的方法提高了计算效率、保证了数值计算时的本质无振荡特性、实现了稳定的计算。

权利要求

1.一种可压缩流体的数值模拟方法，其特征在于，包括如下步骤：步骤S10：读取计算网格、 时刻的物理变量分布、边界条件、试验时长；

步骤S20：在计算网格上，根据 时刻的物理变量获得守恒变量 和对流通量 、、 ，其中， 为笛卡尔坐标系下 方向上的对流通量、 为笛卡尔坐标系下方向上的对流通量； 为笛卡尔坐标系下 方向上的对流通量；

步骤S30：分别计算对流通量 、 、 在物理空间上的数值通量 、 、 ；

步骤S40：将守恒变量 和数值通量 、 、 代入离散欧拉方程得到 时刻的守恒变量 及其对应的物理变量；

步骤S50： 从1 k遍历，其中k为 时刻试验时长的序号，将守恒变量 及 时刻~

的物理变量作为最终的试验结果；

其中，步骤S30包括如下步骤：步骤S301：构建对流通量 、 、 在物理空间上的正负混合通量 、 、；

步骤S302：对物理空间上的正负混合通量 、 、 进行间断检测，若检测到正负混合通量 、 、 均未存在间断，将对流通量 、 、 进行分裂和逼近，得到物理空间上的数值通量 、 、 ；若检测到正负混合通量 、 、存在间断，将对流通量 、 、 投影至特征空间，并分裂为正负通量 、、 ，并根据特征空间上的正负通量 、 、 计算得到物理空间上的数值通量 、 、 ；

步骤S302中，根据特征空间上的正负通量 、 、 ，计算得到物理空间上的数值通量 、 、 的步骤如下：步骤S3021：对特征空间上的正负通量 、 、 进行间断检测，若未检测到正负通量 、 、 存在间断，则采用加权本质无振荡格式对应的高阶线性格式对正负通量 、 、 进行逼近，得到特征空间上的数值通量；若检测到正负通量 、 、 存在间断，则采用加权本质无振荡格式对正负通量 、 、进行逼近，得到特征空间上的数值通量；

步骤S3022：将特征空间上的数值通量投影回物理空间，得到物理空间上的数值通量、 、 ；

步骤S302中和步骤S3021中正负混合通量和正负通量的间断检测的方法如下：获取对流通量 、 、 在物理空间上的正负混合通量 、 、 和特征空间上的正负通量 、 、 中的任意一个通量的分量，并将所述分量作为目标分量 ；将 在计算网格的网格单元 中心处的离散值记为 ，遍历所有计算网格单元 ，对目标分量 在网格单元 附近依次进行间断检测；=1,2,...,N，N为网格单元的总个数；

获取目标分量 在网格单元 对应的加权本质无振荡格式大模板 中的最大斜率值 ；

获取目标分量 在网格单元 对应的加权本质无振荡格式阈值模板 中的斜率阈值 ，其中加权本质无振荡格式阈值模板 为加权本质无振荡格式大模板 的子模板中任一迎风子模板；

当 < 时，则检测到目标分量 在网格单元 附近不存在间断；当斜率阈值时，则检测到目标分量 在网格单元 附近存在间断。

2.如权利要求1所述的可压缩流体的数值模拟方法，其特征在于，步骤S302中，若检测到正负混合通量 、 、 均未存在间断，将对流通量 、 、 进行分裂和逼近，得到物理空间上的数值通量 、 、 ，获得物理空间上数值通量 、 、的方法如下：

采用通量分裂格式将对流通量 、 、 分裂为正负通量 、 、 ；对正负通量 、 、 采用加权本质无振荡格式对应的高阶线性格式进行逼近，获得物理空间中的数值通量 、 、 。

3.如权利要求1所述的可压缩流体的数值模拟方法，其特征在于，步骤S302，将对流通量 、 、 投影至特征空间，并分裂为正负通量 、 、 的方法如下：将对流通量 、 、 由物理空间投影至特征空间，得到特征空间上对流通量、 、 ；采用通量分裂格式对特征空间上的对流通量 、 、 分裂得到正负通量 、 、 。

4.如权利要求1所述的可压缩流体的数值模拟方法，其特征在于，获取目标分量 在网格单元 对应的加权本质无振荡格式大模板 中的最大斜率值 的计算方法如下：

在网格单元 对应的加权本质无振荡格式的大模板 上对目标分量 采用加权本质无振荡格式对应的高阶线性格式进行逼近，得到逼近多项式函数 ；

计算逼近多项式函数 在大模板 上的最大斜率值 。

5.如权利要求1所述的可压缩流体的数值模拟方法，其特征在于，获取目标分量 在网格单元 对应的加权本质无振荡格式阈值模板 中的斜率阈值 的计算方法如下：

选取网格单元 对应的加权本质无振荡格式中的阈值模板 ；

在加权本质无振荡格式阈值模板 中计算基于网格最小可分辨波长的目标分量 的斜率阈值 ，其中，网格最小可分辨波长为网格单元长度的两倍。

6.如权利要求4所述的可压缩流体的数值模拟方法，其特征在于，计算逼近多项式函数在大模板 上的最大斜率值 的公式如下：其中， 为逼近多项式函数 在大模板 上的导数。

7.如权利要求1‑5之一所述的可压缩流体的数值模拟方法，其特征在于，计算斜率阈值的公式如下：

其中， 为网格单元长度， 和 为加权本质无振荡格式中阈值模板 对应的整数，m为网格单元编号， 和 为网格单元的中心位置处目标分量 的离散值， 为正数小量。

说明书

一种可压缩流体的数值模拟方法

技术领域

[0001] 本发明涉及计算流体力学领域 ，尤其是涉及一种可压缩流体的数值模拟方法。

背景技术

[0002] 研究关于飞行器三维复杂流动的数值模拟问题是目前计算流体力学的前沿热点问题，同时也是工程应用中的重要需求，计算流体力学为低成本地研究真实流动现象提供
了巨大支持。和传统的风洞试验相比较，计算流体力学技术具有计算成本低、飞行大气环境
仿真度高的优点，可以以较低计算代价获得较为精确的气动特性，为后续飞行器的相关设
计提供快速和准确的指导。

[0003] 对流扩散方程是粘性流体力学及其他实际问题中常见的一类基本控制方程。对于这类方程，除了极少数简单情形，大部分问题目前还无法求得精确解，只能利用数值方法进
行数值模拟。

[0004] 流体控制方程组是指基于连续介质假设由质量守恒、动量守恒、能量守恒定律推导得到的连续方程、动量方程和能量方程组成的方程组，是常规流动问题的核心数学描述。
对于不同流动问题，流动控制方程组可以简化省略部分项，或增加其他方程。特别地，对于
可压缩无粘流动，流动控制方程组简化为欧拉方程组。

[0005] 加权本质无振荡格式（Weighted Essentially Non‑oscillatory， 简称WENO）方法是近年来广泛流行的一种高分辨率数值方法，用于解决以对流为主的对流扩散方程，特
别是双曲守恒律问题。加权本质无振荡格式具有在光滑区可以得到任意高阶精度，在间断
区保持稳定的、本质无振荡性质和陡峭的间断过渡等优点。

[0006] 传统的计算流体力学在求解含激波、接触间断、相界面等间断流场结构的流动问题时，一阶格式会引起间断结构的过度抹平，高阶线性格式则会在间断附近产生非物理的
数值振荡，均难以实现有效模拟。加权本质无振荡格式方法对子模板上低阶格式进行非线
性加权，既可以实现对变量光滑区域的高精度计算，又可以保证间断附近的本质无振荡特
性，捕捉到陡峭的间断过渡，被广泛应用于含激波等间断结构的数值计算中。但经典的加权
本质无振荡格式在整个计算区域进行非线性加权，计算量较大；可采用一定的判定准则检
测间断，仅在变量间断附近进行加权本质无振荡格式非线性加权计算，提高格式效率。

[0007] 采用加权本质无振荡格式求解欧拉方程组时，可以采用分量型加权本质无振荡格式，即在物理空间对分量形式方程进行求解，但这样仍会产生数值振荡；也可以采用特征型
加权本质无振荡格式，即在特征空间对特征形式方程进行求解，这样可更好地抑制振荡，但
涉及特征变换等矩阵计算，计算量过大。

发明内容

[0008] 本发明实施例的目的是提供一种可压缩流体的数值模拟方法，目的在于解决现有技术中存在的技术问题，以欧拉方程组为例展开说明，包括如下步骤：

[0009] 步骤S10：读取计算网格、 时刻的物理变量分布、边界条件、试验时长；

[0010] 步骤S20：在计算网格上，根据 时刻的物理变量获得守恒变量 和对流通量、 、 ，其中， 为笛卡尔坐标系下x方向上的对流通量、 为笛卡尔坐标系下y方
向上的对流通量； 为笛卡尔坐标系下z方向上的对流通量；

[0011] 步骤S30：分别计算对流通量 、 、 在物理空间上的数值通量 、 、 ；

[0012] 步骤S40：将守恒变量 和数值通量 、 、 代入离散欧拉方程得到n+1时刻的守恒变量 及其对应的物理变量；

[0013] 步骤S50：n从1 k遍历，其中k为 时刻试验时长的序号，将守恒变量 及 时刻~
的物理变量作为最终的试验结果。

[0014] 进一步地，步骤S30包括如下步骤：

[0015] 步骤S301：构建对流通量 、 、 在物理空间上的正负混合通量 、 、

[0016] 步骤S302：对物理空间上的正负混合通量 、 、 进行间断检测，若检测到正负混合通量 、 、 均未存在间断，将对流通量 、 、 进行分裂和逼近，
得到物理空间上的数值通量 、 、 ；若检测到正负混合通量 、 、 存在间
断，将对流通量 、 、 投影至特征空间，并分裂为正负通量 、 、 ，并根据特
征空间上的正负通量 、 、 计算得到物理空间上的数值通量 、 、 。

[0017] 进一步地，步骤S302中，根据特征空间上的正负通量 、 、 ，计算得到物理空间上的数值通量 、 、 的步骤如下：

[0018] 步骤S3021：对特征空间上的正负通量 、 、 进行间断检测，若未检测到正负通量 、 、 存在间断，则采用加权本质无振荡格式对应的高阶线性格式对正
负通量 、 、 进行逼近，得到特征空间上的数值通量；若检测到正负通量 、
、 存在间断，则采用加权本质无振荡格式对正负通量 、 、 进行逼近，得
到特征空间上的数值通量；

[0019] 步骤S3022：将特征空间上的数值通量投影回物理空间，得到物理空间上的数值通量 、 、 。

[0020] 进一步地，步骤S302中，若检测到正负混合通量 、 、 均未存在间断，将对流通量 、 、进行分裂和逼近，得到物理空间上的数值通量 、 、 ，获得物理空
间上数值通量 、 、 的方法如下：

[0021] 采用通量分裂格式将对流通量 、 、 分裂为正负通量 、 、 ；对正负通量 、 、 采用加权本质无振荡格式对应的高阶线性格式进行逼近，获得
物理空间中的数值通量 、 、 。

[0022] 进一步地，步骤S302，将对流通量 、 、 投影至特征空间，并分裂为正负通量 、 、 的方法如下：

[0023] 将对流通量 、 、 由物理空间投影至特征空间，得到特征空间上对流通量、 、 ；采用通量分裂格式对特征空间上的对流通量 、 、 分裂得到正负通量
、 、 。

[0024] 进一步地，步骤S302中和步骤S3021中的通量间断检测的方法如下：

[0025] 获取对流通量正负混合通量 、 、 和特征空间上的正负通量 、 、中的任意一个通量的分量，并将所述分量作为目标分量g；将g在计算网格的网格单元i
中心处的离散值记为gi，遍历所有计算网格单元i，对目标分量g在网格单元i附近依次进行
间断检测；i=1,2,...,N，N为网格单元的总个数；

[0026] 获取目标分量g在网格单元i对应的加权无本质振荡格式大模板Swhole中的最大斜率值Si,max；

[0027] 获取目标分量g在网格单元i对应的加权本质无振荡格式阈值模板Sthreshold中的斜率阈值Si,threshold，其中加权无本质振荡格式阈值模板Sthreshold为加权无本质振荡格式大模
板Swhole的子模板中任一迎风子模板；

[0028] 当Si,max

[0029] 进一步地，获取目标分量g在网格单元i对应的加权无本质振荡格式大模板Swhole中的最大斜率值Si,max的计算方法如下：

[0030] 在网格单元i对应的加权无本质振荡格式的大模板Swhole上对目标分量g采用加权本质无振荡格式对应的高阶线性格式进行逼近，得到逼近多项式函数 ；

[0031] 计算逼近多项式函数 在大模板Swhole上的斜率的绝对值的最大值Si,max。

[0032] 进一步地，获取目标分量g在网格单元i对应的加权本质无振荡格式阈值模板Sthreshold中的斜率阈值Si,threshold的计算方法如下：

[0033] 选取网格单元i对应的加权本质无振荡格式中的阈值模板Sthreshold；

[0034] 在加权本质无振荡格式阈值模板Sthreshold中计算基于网格最小可分辨波长的目标分量g的斜率阈值Si,threshold，其中，网格最小可分辨波长为网格单元长度的两倍。

[0035] 进一步地，计算逼近多项式函数 在大模板Swhole上的斜率绝对值的最大值Si,max的公式如下：

[0036]

[0037] 其中， 为逼近多项式函数 在大模板Swhole上的导数。

[0038] 进一步地，计算斜率阈值Si,threshold的公式如下：

[0039]

[0040] 其中， 为网格单元长度，r和b为加权本质无振荡格式中阈值模板Sthreshold对应的整数，m为网格单元编号，gm+1和gm为网格单元的中心位置处目标分量g的离散值， 为正
数小量。

[0041] 本发明产生的技术效果至少具有以下方面：

[0042] 1）通过本发明提供的可压缩流体的数值模拟方法，该数值模拟方法在求解欧拉方程组的过程中进行了两次间断检测，在未检测到间断的区域采用加权本质无振荡格式对应
的高阶线性格式进行计算，简化了非必须的加权本质无振荡格式非线性加权计算，实现了
加权本质无振荡格式整体计算效率的提高；在检测到间断的区域采用加权本质无振荡格式
计算，保证数值计算时的本质无振荡特性，实现稳定的计算。需要说明的是本发明的方法主
要适用于有限差分加权本质无振荡格式，但易于扩展至有限体积加权本质无振荡格式。

[0043] 2）本发明的可压缩流体的数值模拟方法中，在第一次对对流通量在物理空间上进行间断检测时，对于检测到的正负混合通量均未存在间断时，直接在物理空间中对对流通
量进行分裂和逼近计算，不对对流通量进行特征投影等耗时巨大的矩阵运算，减小了试验
时的计算量，提高了试验效率。

[0044] 3）本发明的可压缩流体的数值模拟方法中，在物理空间和特征空间进行了两次间断检测后的计算过程中，仅针对在特征空间中检测出间断的通量进行加权本质无振荡格式
的逼近，大大减小了加权本质无振荡非线性加权计算的计算量。

[0045] 4）本发明的可压缩流体的数值模拟方法中，在物理空间进行间断检测时，构建了正负混合通量，进行间断检测时仅是对正负混合通量进行间断检测，而不是对原对流通量
中的各个分量分别进行间断检测，减少了多次间断检测的计算量。

[0046] 5）本发明的可压缩流体的数值模拟方法中，在进行通量间断检测的方法中将斜率阈值作为间断检测的标准，除了正数小量外，不包含任何经验参数，具有更好的普适性，同
时，斜率阈值与所研究的网格单元所处的阈值模板内通量差值的绝对值正相关，对于不同
的网格单元取值不同，具有更好的自适应性。

附图说明

[0047] 图1是本发明可压缩流体的数值模拟方法的总体流程图 ；

[0048] 图2是本发明中通量进行间断检测的流程图；

[0049] 图3是本发明中 5阶加权本质无振荡格式大模板与高阶线性格式逼近多项式函数示意图；

[0050] 图4是本发明5阶加权本质无振荡格式各子模板与斜率阈值示意图。

具体实施方式

[0051] 以下的说明提供了许多不同的实施例、或是例子，用来实施本发明的不同特征。以下特定例子所描述的元件和排列方式，仅用来精简的表达本发明，其仅作为例子，而并非用
以限制本发明。

[0052] 如图1‑4所示，本发明实施例的目的是提供一种可压缩流体的数值模拟方法，目的在于解决现有技术中存在的技术问题，包括如下步骤：

[0053] 步骤S10：读取计算网格、时 刻的物理变量分布、边界条件、试验时长；

[0054] 步骤S20：在计算网格上，根据 时刻的物理变量获得守恒变量 和对流通量、 、 ，其中， 为笛卡尔坐标系下x方向上的对流通量、 为笛卡尔坐标系下y方
向上的对流通量； 为笛卡尔坐标系下z方向上的对流通量；

[0055] 步骤S30：分别计算对流通量 、 、 在物理空间上的数值通量 、 、 ；

[0056] 步骤S40：将守恒变量 和数值通量 、 、 代入离散欧拉方程得到n+1时刻的守恒变量 及其对应的物理变量；

[0057] 步骤S50：n从1 k遍历，其中k为 时刻试验时长的序号，将守恒变量 及 时刻~
的物理变量作为最终的试验结果。

[0058] 对流扩散方程是粘性流体力学及其他实际问题中常见的一类基本控制方程。对于这类方程，除了极少数简单情形，大部分问题目前还无法求得精确解，只能利用数值方法进
行数值模拟。

[0059] 流体控制方程组是常规流动问题的核心数学描述，具体是指基于连续介质假设由质量守恒、动量守恒、能量守恒定律推导得到的连续方程、动量方程和能量方程组成的方程
组。对于不同流动问题，流动控制方程组可以简化省略部分项，或者增加其他方程。特别地，
对于可压缩无粘流动，流动控制方程组简化为欧拉方程组，欧拉方程组的形式如下：

[0060]

[0061] ， ， ， ，

[0062] 其中：U为守恒变量；F、G、H分别为笛卡尔坐标系下x、y、z三个方向上的对流通量；S为源项； 为密度；u、v、w分别为笛卡尔坐标系下x、y、z方向上的速度分量；E为单位质量流
体总能；p为压力

[0063] 在均匀结构网格上，欧拉方程组的半离散形式如下：

[0064]

[0065] 其中，i、j、k为网格单元编号； 、 、 为x、y、z方向上对流通量对应的离散后的数值通量， 、 、 为x、y、z方向的网格单元长度，Si,j,k为网格单元中
心源项值。

[0066] 计算网格具体为：根据设计时的任务要求，构建需要设计对象的模型，所述模型为不同种类飞行器的模型，并建立复杂流场计算所需要的网格，由给定的网格得到需要的网
格单元信息，网格单元信息包括网格尺度、节点坐标等。

[0067] 边界条件指的是计算区域边界的物理变量的分布和取值。

[0068] 试验时长指的是本次试验需要的时间。

[0069] 具体地，物理变量指的是 、u、v、w、p，物理变量通过简单的计算后可以获得守恒变量 和对流通量 、 、 。

[0070] 上述方案中，根据物理空间上的对流通量 、 、 计算出物理空间上的数值数值通量 、 、 ，然后将守恒变量 和数值通量 、 、 代入离散欧拉方程得
到n+1时刻的守恒变量 及其对应的物理变量；n从1 k遍历，其中k为 时刻试验时长的
~
序号，将守恒变量 及 时刻的物理变量作为最终的试验结果。

[0071] 物理空间上的对流通量 、 、 计算出物理空间上的数值数值通量 、 、的过程中，采用了两次间断检测，包括物理空间间断检测和特征空间间断检测，在物理
空间间断检测时，首先构建对流通量 、 、 在物理空间上的正负混合通量 、
、 ，然后对正负混合通量 、 、 进行间断检测，若检测到正负混合通量 、
、 均未存在间断，则直接对 、 、 进行分裂和逼近，得到物理空间上的数值通
量 、 、 ，若检测到正负混合通量 、 、 存在间断，将对流通量 、 、
投影至特征空间，并分裂为正负通量 、 、 ，并根据特征空间上的正负通量
、 、 ，对特征空间上的正负通量 、 、 进行间断检测，计算得到物理空
间上的数值通量 、 、 。

[0072] 因此，通过本发明提供的可压缩流体的数值模拟方法，该数值模拟方法在求解欧拉方程组的过程中进行了两次间断检测，在未检测到间断的区域采用加权本质无振荡格式
对应的高阶线性格式进行计算，简化了非必须的加权本质无振荡非线性加权计算，实现了
加权本质无振荡格式整体计算效率的提高；在检测到间断的区域采用加权本质无振荡格式
计算，保证数值计算时的本质无振荡特性，实现稳定的计算。需要说明的是本发明的方法主
要适用于有限差分加权本质无振荡格式，但易于扩展至有限体积加权本质无振荡格式。

[0073] 进一步地，步骤S30包括如下步骤：

[0074] 步骤S301：构建对流通量 、 、 在物理空间上的正负混合通量 、 、；

[0075] 步骤S302：对物理空间上的正负混合通量 、 、 进行间断检测，若检测到正负混合通量 、 、 均未存在间断，将对流通量 、 、 进行分裂和逼近，
得到物理空间上的数值通量 、 、 ；若检测到正负混合通量 、 、 存在间
断，将对流通量 、 、 投影至特征空间，并分裂为正负通量 、 、 ，并根据特
征空间上的正负通量 、 、 计算得到物理空间上的数值通量 、 、 。

[0076] 上述方案中，特征空间上的正负混合通量 的计算方式如下：

[0077]

[0078] 其中， ，ci为网格单元i中心处音速。

[0079] 正负混合通量 、 的计算方式与正负混合通量 的计算方式相似。

[0080] 在间断检测时，检测到正负混合通量 、 、 均未存在间断指的是正混合通量 、 、 和负混合通量 、 、 两者都没有检测到间断，此时直
接对对流通量 、 、 进行分裂和逼近，得到物理空间上的数值通量 、 、 ；当
正混合通量 、 、 和负混合通量 、 、 中任何一个检测到间断，
则说明正负混合通量 、 、 存在间断，此时对流通量 、 、 投影至特征空
间进行进一步的间断检测。

[0081] 上述方案中，在第一次对对流通量 、 、 在物理空间上进行间断检测时，对于检测到的正负混合通量 、 、 均未存在间断时，直接在物理空间中对对流通量
、 、 进行分裂和逼近计算，不用将对流通量 、 、 进行特征投影等耗时巨大
的矩阵运算，减小了试验时的计算量，提高了试验效率。

[0082] 在物理空间对对流通量 、 、 进行间断检测时，构建了正负混合通量 、、 ，进行间断检测时仅是对正负混合通量 、 、 进行间断检测，而不是针
对原对流通量 、 、 中的各个分量分别进行间断检测，减少了多次间断检测的计算
量。

[0083] 进一步地，步骤S302中，根据特征空间上的正负通量 、 、 ，计算得到物理空间上的数值通量 、 、 的步骤如下：

[0084] 步骤S3021：对特征空间上的正负通量 、 、 进行间断检测，若未检测到正负通量 、 、 存在间断，则采用加权本质无振荡格式对应的高阶线性格式对正
负通量 、 、 进行逼近，得到特征空间上的数值通量；若检测到正负通量 、
、 存在间断，则采用加权本质无振荡格式对正负通量 、 、 进行逼近，得
到特征空间上的数值通量；

[0085] 步骤S3022：将特征空间上的数值通量投影回物理空间，得到物理空间上的数值通量 、 、 。

[0086] 上述方案中，在对特征空间上的正负通量 、 、 进行间断检测时，当未检测到间断时采用加权本质无振荡格式对应的高阶线性格式对正负通量 、 、 进行
逼近，当检测到间断时，则采用加权本质无振荡格式对正负通量 、 、 进行逼近，
在此过程中，仅针对在特征空间上检测出间断的通量进行加权本质无振荡格式的逼近，大
大减小了加权本质无振荡非线性加权计算的计算量。

[0087] 进一步地，步骤S302中，若检测到正负混合通量 、 、 均未存在间断，将对流通量 、 、 进行分裂和逼近，得到物理空间上的数值通量 、 、 ，获得物
理空间上数值通量 、 、 的方法如下：

[0088] 采用通量分裂格式将对流通量 、 、 分裂为正负通量 、 、 ；对正负通量 、 、 采用加权本质无振荡格式对应的高阶线性格式进行逼近，获得
物理空间中的数值通量 、 、 。

[0089] 上述方案中，通量分裂格式优选为Lax–Friedrichs格式，也可以选用其他的通量分裂格式，在此不做限制。

[0090] 进一步地，步骤S302，将对流通量 、 、 投影至特征空间，并分裂为正负通量 、 、 的方法如下：

[0091] 将对流通量 、 、 由物理空间投影至特征空间，得到特征空间上对流通量、 、 ；采用通量分裂格式对特征空间上的对流通量 、 、 分裂得到正负通量
、 、 。

[0092] 上述方案中，通量分裂格式优选为Lax–Friedrichs格式，也可以选用其他的通量分裂格式，在此不做限制。

[0093] 进一步地，步骤S302中和步骤S3021中的通量间断检测的方法如下：

[0094] 获取对流通量正负混合通量 、 、 和特征空间上的正负通量 、 、中的任意一个通量的分量，并将所述分量作为目标分量g；将g在计算网格的网格单元i
中心处的离散值记为gi，遍历所有计算网格单元i，对目标分量g在网格单元i附近依次进行
间断检测；i=1,2,...,N，N为网格单元的总个数；

[0095] 获取目标分量g在网格单元i对应的加权无本质振荡格式大模板Swhole中的最大斜率值Si,max；

[0096] 获取目标分量g在网格单元i对应的加权本质无振荡格式阈值模板Sthreshold中的斜率阈值Si,threshold，其中加权无本质振荡格式阈值模板Sthreshold为加权无本质振荡格式大模
板Swhole的子模板中任一迎风子模板；

[0097] 当Si,max

[0098] 进一步地，获取目标分量g在网格单元i对应的加权无本质振荡格式大模板Swhole中的最大斜率值Si,max的计算方法如下：

[0099] 在网格单元i对应的加权无本质振荡格式的大模板Swhole上对目标分量g采用加权本质无振荡格式对应的高阶线性格式进行逼近，得到逼近多项式函数 ；

[0100] 计算逼近多项式函数 在大模板Swhole上的斜率的绝对值的最大值Si,max。

[0101] 进一步地，获取目标分量g在网格单元i对应的加权本质无振荡格式阈值模板Sthreshold中的斜率阈值Si,threshold的计算方法如下：

[0102] 选取网格单元i对应的加权本质无振荡格式中的阈值模板Sthreshold；

[0103] 在加权本质无振荡格式阈值模板Sthreshold中计算基于网格最小可分辨波长的目标分量g的斜率阈值Si,threshold，其中，网格最小可分辨波长为网格单元长度的两倍。

[0104] 进一步地，计算逼近多项式函数 在大模板Swhole上的斜率绝对值的最大值Si,max的公式如下：

[0105]

[0106] 其中， 为逼近多项式函数 在大模板Swhole上的导数。

[0107] 进一步地，计算斜率阈值Si,threshold的公式如下：

[0108]

[0109] 其中， 为网格单元长度，r和b为加权本质无振荡格式中阈值模板Sthreshold对应的整数，m为网格单元编号，gm+1和gm为网格单元的中心位置处目标分量g的离散值， 为正
数小量。

[0110] 上述方案中，如图3‑4所示，本发明实施例中以5阶的加权本质无振荡格式为例进行说明，其中5阶加权本质无振荡格式构造模板包括大模板 和三个子模板S0、S1、S2，则可
取 为大模板Swhole，S1为阈值模板；Sthreshold其中x为横坐标，xi为网格单元的中心，gi为目
标分量g在网格单元的中心xi处的值。

[0111] 最大斜率阈值Si,max的计算方法：在网格单元i应的加权无本质振荡格式的大模板Swhole上对目标分量g采用加权本质无振荡格式对应的高阶线性格式进行逼近，得到逼近多
项式函数 ，在附3图中，虚线表示的为逼近多项式函数 ；然后通过求导计算逼近多
项式函数 的斜率值 ，取斜率值 的绝对值的最大值作为目标分量g在网格单
元i对应的加权无本质振荡格式大模板Swhole中的最大斜率值Si,max；

[0112] 在加权本质无振荡格式所有子模板中，取模板中心距界面 最近的迎风子模板为阈值模板Sthreshold，在Sthreshold上计算基于网格最小可分辨波长（即2倍的网格单元长
度）的斜率阈值Si,threshold。记阈值模板Sthreshold为 ，
为第i个网格单元，r和b为与所采用格式相关的整数，在附图4中，弯曲的
虚线的斜率的绝对值的最大值表示的是斜率阈值Si,threshold。

[0113] 本发明中间断检测方法的基本原理为：采用高阶线性格式对应的逼近多项式函数进行数值模拟时，在间断附近会产生非物理的数值振荡，这些振荡伴随有变量的大斜
率变化，这使得可以用 的斜率大小表征变量是否出现间断；而在构建的阈值模板
Sthreshold上的斜率阈值Si,threshold时，Si,threshold为波长为网格最小可分辨波长（即2倍的网格
单元长度）的正弦波在幅值为 时的斜率最大值，用以表示Sthreshold中允许
出现的最大斜率。这里幅值之所以取迎风的阈值模板Sthreshold中的相邻网格通量差值绝对
值的最小值 ，是为了减小Si,threshold的值，即减小模板内所允许的最大斜
率，使得计算更为稳定。最后，判断 对应的斜率值Si,max是否大于Si,threshold，当Si,max<
Si,threshold时，则检测到目标分量g在网格单元i附近不存在间断；当最大斜率值Si,max
Si,threshold斜率阈值时，则检测到g在网格单元i附近存在间断。

[0114] 在进行通量间断检测的方法中将斜率阈值作为间断检测的标准，除了正数小量外，不包含任何经验参数，具有更好的普适性，同时，斜率阈值与所研究的网格单元所处的
阈值模板内通量差值的绝对值正相关，对于不同的网格单元取值不同，具有更好的自适应
性。

[0115] 以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。