具体实施方式

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

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

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

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

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

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

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

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

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

，，，，

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

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

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

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

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

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

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

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

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

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

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

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

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

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

其中，，c _i 为网格单元i中心处音速。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

当S _i,max<S _i , _threshold时，则检测到目标分量g在网格单元i附近不存在间断；当S _i,max S _i , _threshold时，则检测到目标分量g在网格单元i附近存在间断。

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

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

计算逼近多项式函数在大模板S _whole上的斜率的绝对值的最大值S _i,max。

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

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

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

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

其中，为逼近多项式函数在大模板S _whole上的导数。

进一步地，计算斜率阈值S _i , _threshold的公式如下：

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

上述方案中，如图3-4所示，本发明实施例中以5阶的加权本质无振荡格式为例进行说明，其中5阶加权本质无振荡格式构造模板包括大模板和三个子模板S ₀ 、S ₁ 、S ₂ ，则可取为大模板S _whole，S ₁ 为阈值模板；S _threshold其中x为横坐标，x _i 为网格单元的中心，g _i 为目标分量g在网格单元的中心x _i 处的值。

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

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

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

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

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