CN108280273A - 一种基于非等距网格下的有限体积流场数值计算方法 - Google Patents

Abstract

一种基于非等距网格下的有限体积流场数值计算方法，涉及了一种线性权任取的五阶有限体积WENO格式，在不降低整体数值精度的前提下，只需任取三个加起来为等于1的正数作为线性权，然后进行后续计算即可。本发明构造的数值格式由于不需要计算非等距网格下跟网格尺度有关的复杂线性权，所以更易于推广到更高阶的WENO格式、移动网格或是自适应网格中的流场计算。此外，本发明同传统的WENO格式一样采用稳定性好，操作简单和易于编程的Runge‑Kutta时间离散方法进行求解。最后，本发明通过具体的数值算例与传统格式进行比较从而验证其在流场计算中的优越性。

Description

一种基于非等距网格下的有限体积流场数值计算方法

技术领域

本发明公开了一种基于非等距网格下的有限体积流场数值计算方法，可用于计算数学、空气动力学、磁流体力学、航空航天、导弹发射、汽车工业、土木工程、环境工程、船舶、气象工程等技术领域。

背景技术

近几十年，双曲守恒律方程的高阶数值求解格式一直是重点的研究内容，其应用领域包括计算流体力学、计算天文学、半导体模拟、交通流模型以及磁流体力学等。求解这些问题的方程主要困难在于其解的复杂性，即使初始值光滑也会经过一段时间的演变产生激波、接触间断和稀疏波等复杂解结构。另外，低阶数值格式会抹平这些间断，不能很好的分辨出这些复杂解结构，所以针对此类问题的高阶数值格式的研究非常有必要。1987年，Harten和Osher首次提出了ENO(Essentially Non-Oscillatory基本无振荡)格式并成功应用于相关领域。ENO格式的主要思想是在所有候选模板中选出最光滑的重构目标单元的值，从而实现高精度和基本无振荡的性质。而因为其存在模板选择过程从而使得计算过程中会造成计算结果的浪费，所以在此基础上构造出了WENO(Weighted ENO)格式。1994年，WENO格式首次被提出并构造出一维情形的三阶有限体积形式。1996年，Jiang和Shu构造了3阶和5阶有限差分WENO格式，并给出了WENO格式的一般计算框架，包括重构多项式的构造，线性权的计算，光滑指示器的计算以及非线性权和最终重构值的计算。基于Jiang和Shu的工作，后续很多学者参与到有限差分/有限体积WENO格式的研究中来，比如：Hu，Qiu等人构造和实现了非结构网格下的有限体积WENO格式和中心WENO(Central WENO)格式；Balsara和Shu等人构造出了7到11阶等更高阶的有限差分WENO格式；Qiu在WENO格式的基础上提出HWENO(Hermite WENO)格式等。此外，其应用领域也非常广泛，比如：Dietrich，Bernuzzi，Grimm-Strele，Kupka，Leung等人应用于天体物理学和地球物理学；Fierro，Pressel等应用于气候科学。相比于ENO格式，WENO格式将所有候选模版进行凸组合，从而避免了计算结果的浪费，并且保持了ENO格式具有的基本无振荡的性质以及提高了数值格式的整体精度。

需要注意的是，WENO格式的主要内容在于重构多项式的构造、线性权的选取、光滑指示器以及非线性权的计算。而这些计算过程在等距网格中相对比较简单，并且其具体的表达式也很容易得到。但是一旦推广到非等距网格中，其构造相对复杂，因为其具体的表达式不仅跟网格点上的物理量有关还与网格尺度有关。此外，不管是等距网格还是非等距网格，上述过程中的线性权必须是最优线性权才能使整体格式达到预期设计的精度，不然会使精度降低甚至容易出现负线性权导致格式不稳定。对于此，Zhu等设计出了一种线性权可以任取，并且能够保证格式整体精度不丢失的有限体积和有限差分WENO格式。

基于此思想，本发明构造了一种非等距网格下的线性权可以任取的有限体积基本加权无振荡格式，并结合3阶满足TVD性质的Runge-Kutta时间离散方法数值求解并得到某一区域内的流场解结构。本发明在随机网格、周期变化网格以及伸缩网格等非等距笛卡尔网格中，将基于Runge-Kutta时间离散的新的线性权可以任取的有限体积WENO格式求解流动问题并通过数值实验来验证方法的可行性。

发明内容

本发明针对现有技术中的不足，提供一种基于非等距网格下的有限体积流场数值计算方法，能够在一般的非等距笛卡尔网格下利用高精度计算方法对流场内的流动问题进行数值计算模拟。

为实现上述目的，本发明采用以下技术方案：

一种基于非等距网格下的有限体积流场数值计算方法，其特征在于，在非等距不均匀笛卡尔网格中，构造5阶精度的基本加权无振荡格式进行流场内的数值计算模拟，具体包括以下步骤：

步骤1、在流场数值计算区域内，生成非等距笛卡尔坐标系；

步骤2、在建立的非等距笛卡尔坐标系中构造有限体积WENO格式，其线性权可以任意取值且与网格尺度和类型无关，并且能保持格式的5阶精度，随后结合Runge-Kutta时间离散方法进而求解无粘流体的控制方程组，得到流场内所有非等距网格单元点的物理量的高精度数值近似值。

为优化上述技术方案，采取的具体措施还包括：

所述步骤1中，非等距笛卡尔坐标系具体为：

随机网格：

Δx_i＝0.5Δx+(1.5Δx-0.5Δx)·random_number(i)，

Δy_j＝0.5Δy+(1.5Δy-0.5Δy)·random_number(j)，

其中，下标i，j分别表示x和y方向的网格序号，Δx_i表示x方向上第i个网格对应的步长，Δy_j表示y方向上第j个网格对应的步长，Δx表示x方向上的平均步长，Δy表示y方向上的平均步长，random_number(i)和random_number(j)分别表示第i个和第j个网格对应的0到1之间的随机数；

或周期变化网格：

其中，L_x，L_y分别表示计算区域x方向和y方向的长度，N_x，N_y分别表示计算区域x方向和y方向的网格单元数；

或伸缩网格：

Δx_i＝Δx_i-1·α，

Δy_j＝Δy_j-1·β，

其中，α和β分别表示为x方向和y方向的伸缩比。

所述步骤2具体包括：

2.1对无粘流体的控制方程组两边同时在目标单元内积分，构造有限体积WENO格式中的空间导数项，进而得到半离散的有限体积格式；

2.2利用Runge-Kutta时间离散方法求解半离散有限体积格式，从而得到控制方程组的时空全离散有限体积格式，其中的一组非线性权可以任意取值，与网格尺度和类型无关，只需满足和为1的正数即可；

2.3时空全离散有限体积格式具体为时间方向上的一个迭代公式，给定流场初始状态的相关物理量，根据迭代公式不断循环，即可得到流场某一时刻或稳定状态时的物理量的高精度近似值。

所述步骤2中，无粘流体的控制方程组为：

其中，t表示时间变量，x，y表示空间变量，U＝(ρ，ρu，ρv，E)^T表示守恒变量，f(U)＝(ρu，ρu²+p，ρuv，u(E+p))^T，g(U)＝(ρv，ρuv，ρv²+p，v(E+p))^T，f(U)，g(U)表示通量，f(U)_x表示f(U)对x求导，g(U)_y表示g(U)对y求导，ρ，p，u，v，E分别表示流体密度、压强、水平方向速度、竖直方向速度以及能量，T表示转置，U₀表示初始状态值。

所述步骤2.1具体包括：

对控制方程两边同时在目标网格单元I_i，j内积分有：

其中，表示守恒变量U(x，y，t)在目标单元I_i，j内的平均值，表示对时间变量t求导，表示对应的被积函数Φ在对应积分区间[a，b]内积分；

构造有限体积WENO格式中的空间导数项，具体过程如下：

2.1.1利用三点Gauss求积公式求解上式右端中的积分项，即

其中，G_m表示三点Gauss求积公式对应的求积节点；

2.1.2利用Lax-Friedrichs数值通量近似代替Gauss求积公式中每个求积节点的流通量，即

其中，表示对应点右侧高阶近似值，表示对应点左侧高阶近似值，表示对应点上侧高阶近似值，表示对应点上侧高阶近似值，α表示通量f(U)或g(U)对守恒变量U求导的最大值；

2.1.3在非等距网格下高阶近似的值，其它点的高阶近似值类似求得，具体如下：

2.1.3.1固定空间变量y，即假设y＝y_k，利用网格平均值得到的高阶近似值，过程如下：

a)选择一个母模板进而得到一个四次重构多项式p₁(x，y_k)，其满足以下条件：

得到其具体表达式为：

其系数满足线性方程组U＝C×A，具体形式如下：

A＝(a₁，a₂，a₃，a₄，a₅)^T，

其中U，A表示向量，C表示系数矩阵，T表示对向量或矩阵转置；

b)选择两个小模板从而得到两个线性重构多项式p₂(x，y_k)，p₃(x，y_k)，其分别满足：

其具体表达式为：

c)任意选择三个小于1且和为1的正数作为线性权γ_n(n＝1，2，3)；

d)计算每个模板对应的光滑指示器β_n，k(n＝1，2，3)，其计算公式如下：

其中，n表示对应多项式的序号，α表示求和指标，r表示对应多项式的次数，表示多项式p_n(x，y_k)对自变量x求α次偏导数；

对应光滑指示器的具体表达式为：

e)根据线性权γ_n和光滑指示器β_n，k求非线性权ω_n，k，公式如下：

其中，为计算过程中的过渡值；

f)得到的重构近似表达式：

2.1.3.2分别令y＝y_k-2，y_k-1，y_k+1，y_k+2，重复步骤2.1.3.1，从而得到

的值；

2.1.3.3固定空间变量x，即分别令x＝x_i+1/2和x＝x_i-1/2，利用步骤2.1.3.1和步骤2.1.3.2得到的值，得到的最终高阶近似表达式：

2.1.3.4重复步骤2.1.3.1到步骤2.1.3.3，类似得到的值；

2.1.4将步骤2.1.3的值代入步骤2.1.2，再代入步骤2.1.1，即可得到流体控制方程组的半离散有限体积格式，并表示为：

所述步骤2.2中，假设时间步长为Δt，tⁿ表示第n时间层，表示的值，在tⁿ时间层进行步骤2.1的过程即可得到然后利用满足TVD性质的Runge-Kutta时间离散方法离散半离散有限体积格式从而得到时空全离散有限体积格式，如下：

其中，和为计算中间的过渡值。

所述步骤2.3中，空间方向和时间方向的离散完成后得到的时空全离散有限体积格式为关于时间层的迭代公式，已知流场初始状态的相关物理量的值，根据迭代公式不断循环，即可得到流场某一时刻或稳定状态时的物理量的高精度近似值。

本发明的有益效果是：在随机网格、周期变化网格以及伸缩网格等非等距笛卡尔网格中构造新的具有5阶精度的线性权可以任取的有限体积WENO格式，结合Runge-Kutta时间离散方法用来数值模拟流场内包含激波、接触间断以及稀疏波等复杂解结构的流动问题，其优势在于：1)在保证格式整体精度的前提下，构造过程中的线性权可以任意取值，与网格尺度无关，只需满足大于0且求和等于1即可；2)避免处理传统有限体积WENO格式中出现的负非线性权；3)能够在各类非等距网格中处理流体流动问题；4)能够更简单的推广到高维问题或是移动网格下的问题。

传统的有限体积WENO数值方法在非等距网格中的线性权必须数最优线性权才能达到预期设计精度，并且其与网格类型和网格尺度有关，计算相对复杂，尤其是推广到高维情形则需要计算每个Gauss求积节点对应的线性权，显得更加复杂繁琐，且极易出现负线性权从而导致格式不稳定，但通过本发明的方法不需要计算非等距网格下跟网格尺度有关的线性权以及每个Gauss求积节点对应的线性权，可以任意取值，从而能够高效的准确的在非等距网格中数值模拟流动问题，从而高分辨率的计算出激波、接触间断以及稀疏波等解结构并同时在解的光滑区域能够保持时空一致高阶精度。

附图说明

图1是三种非等距网格示意图。

图2是实施例二中，在随机网格和周期网格中计算双马赫问题的压强等值线图。

图3是实施例二中，在随机网格和周期网格中计算双马赫问题的密度等值线图。

图4是实施例三中，在随机网格和周期网格中计算对称Riemann问题的密度等值线图。

图5是实施例四中，计算Rayleigh-Taylor不稳定问题的网格示意图。

图6是实施例四中，任取三种线性权计算Rayleigh-Taylor不稳定问题的密度等值线图。

具体实施方式

现在结合附图对本发明作进一步详细的说明。

首先，构造新的非等距网格下有限体积WENO格式。

考虑二维流场控制方程：

其中t表示时间变量，x，y表示空间变量，U＝(ρ，ρu，ρv，E)^T表示守恒变量，f(U)＝(ρu，ρu²+p，ρuv，u(E+p))^T，g(U)＝(ρv，ρuv，ρv²+p，v(E+p))^T，f(U)，g(U)表示通量，f(U)_x表示f(U)对x求导，g(U)_y表示g(U)对y求导，ρ，p，u，v，E分别表示流体密度、压强、水平方向速度、竖直方向速度以及能量，T表示转置，U₀表示初始状态值。值得注意的是，上述所有的物理量都是关于时间变量和空间变量的函数，即随时间和空间位置的变化而变化。

令x_i+1/2-x_i-1/2＝Δx_i，y_j+1/2-y_j-1/2＝Δy_j，网格中心为网格单元为I_i，j＝[x_i-1/2，x_i+1/2]×[y_j-1/2，y_j+1/2]，其中下标i，j为坐标序号。为了方便说明，本发明取随机网格、周期变化网格和伸缩网格三种非等距网格，如图1所示。

定义表示U在网格单元I_i，j内的平均值，即tⁿ表示第n时间层，Δt表示时间步长，则tⁿ⁺¹＝tⁿ+Δt。表示在第n时间层的值，下面的过程全在第n层时间进行，所以暂时不考虑时间变量t，即将简写为

步骤1，在非等距网格下高阶近似的值，G_m表示Gauss节点，表示对应点右侧高阶近似值，表示对应点左侧高阶近似值，具体如下：

步骤1.1，固定空间变量y，假设y＝y_k，利用网格平均值得到的高阶近似值，过程如下：

步骤1.1.1，选择一个母模板进而得到一个四次重构多项式p₁(x，y_k)，其满足以下条件：

得到其具体表达式为

其中a₁，a₂，a₃，a₄，a₅满足线性方程组U＝C×A，具体形式如下：

A＝(a₁，a₂，a₃，a₄，a₅)^T， (5)

步骤1.1.2，选择两个小模板从而得到两个线性重构多项式p₂(x，y_k)，p₃(x，y_k)，其分别满足：

其具体表达式为

步骤1.1.3，任意选择三个小于1的正数(和为1)作为线性权γ_n(n＝1，2，3)，传统的5阶有限体积WENO格式这里必须是最优线性权，且和网格尺度以及对应的目标点有关，计算相对复杂，而本发明在进行数值计算时为了简便全部取为同一组数。

步骤1.1.4，计算每个模板对应的光滑指示器β_n，k(n＝1，2，3)，其计算公式如下：

其中，n表示对应多项式的序号，α表示求和指标，r表示对应多项式的次数，表示多项式p_n(x，y_k)对自变量x求α次偏导数。

对应光滑指示器的具体表达式为

步骤1.1.5，根据线性权γ_n和光滑指示器β_n，k求非线性权ω_n，k，公式如下：

其中，为计算过程中的过渡值，ε为很小的数，防止分母为0，本发明取值10^-6。

步骤1.1.6，最后，得到的重构近似表达式：

步骤1.2，随后，分别令y＝y_k-2，y_k-1，y_k+1，y_k+2，重复步骤1.1，只需将公式(2)-(16)中的下标k分别换成k-2，k-1，k+1，k+2即可，从而类似得到

的值。

步骤1.3，固定空间变量x，分别令x＝x_i+1/2和x＝x_i-1/2，利用步骤1.1和步骤1.2得到的值，重复步骤1.1，只需将这五个值替换步骤1.1中对应的值即可，从而得到的最终高阶近似表达式：

步骤2，在非等距网格下高阶近似的值，G_m表示Gauss节点，表示对应点上侧高阶近似值，表示对应点上侧高阶近似值。类似于步骤1，先固定空间变量x＝x_l，利用网格平均值得到的高阶近似值；其次，分别令x＝x_l-2，x_l-1，x_l+1，x_l+2分别得到 的值；进而固定空间变量y，分别令y＝y_j-1/2和y＝y_j+1/2利用上述五个值重复上述步骤即可得到的最终高阶近似表达式：

步骤3，将步骤1和步骤2得到的和值代入以下Lax-Friedrichs数值通量近似公式代替Gauss求积公式中每个求积节点对应的真实流通量，即

其中，α表示通量f(U)或g(U)对守恒变量U求导的最大值。

步骤4，进一步通过三点Gauss求积公式得到通量在对应区间的积分(即为对应边的流通量)，即将式(17)和(18)代入下式：

其中，G_m表示三点Gauss求积公式对应的求积节点，分别为0；α_m表示三点Gauss求积公式每个节点对应的求积系数，分别为

步骤5，对控制方程(1)两边同时在目标网格单元I_i，j内积分有：

步骤6，将式(19)和(20)代入式(21)从而得到(21)式的右端项的高阶近似值，并将其表示为则得到包含时间导数项的半离散有限体积格式：

步骤7，(24)式为关于时间变量t的常微分方程，利用三阶具有TVD性质的Runge-Kutta方法求解，如下：

其中，和为计算中间的过渡值。

步骤8，(24)和(25)式联合得到了时空全离散格式，即为关于时间t的迭代公式，假设流场物理量的初始状态已知，则可以循环迭代从而得到某一时刻或是稳定流场各物理量的值，达到数值模拟流场流动的目的。

下面给出四个算例作为本发明所公开方法的具体实施例。

实施例一、精度测试问题。考虑方程组(1)，初始条件为：ρ(x，y，0)＝1+0.2sin(π(x+y))，u(x，y，0)＝0.7，v(x，y，0)＝0.3，p(x，y，0)＝1。精确解为：ρ(x，y，t)＝1+0.2sin(π(x+y-t))，u(x，y，t)＝0.7，v(x，y，t)＝0.3，p(x，y，t)＝1。计算终止时刻为t＝2。计算区域为(x，y)∈[0，2]×[0，2]。为了说明本发明数值计算方法中的线性权可以任意取值，本发明只取五组不同规律的数作为线性权：(1)γ₁＝0.98，γ₂＝0.01，γ₃＝0.01；(2)γ₁＝0.495，γ₂＝0.495，γ₃＝0.01；(3)γ₁＝1/3，γ₂＝1/3，γ₃＝1/3；(4)γ₁＝0.001，γ₂＝0.001，γ₃＝0.998；(5)γ₁＝0.3，γ₂＝0.4，γ₃＝0.3。

表1任取五种线性权的数值精度与传统WENO格式的比较

实施例二、双马赫反射问题。如图2、图3所示，在距左边界1/6处有一马赫数为10的强激波，以与x轴夹角为60°射入，左右状态值分别为：(ρ_L，u_L，v_L，E_L)＝(8，7.145，-4.125，563.544)，(ρ_R，u_R，v_R，E_R)＝(1.4，0，0，2.5)。CFL数取0.6，计算网格300×100，终止时刻t＝0.2。

实施例三、径向对称Riemann问题。如图4所示，考虑流体力学方程组(1)，取初始条件为：

计算区域为[0，1]×[0，1]，CFL数取0.6，计算网格100×100，计算终止时刻为t＝0.13。

实施例四、Rayleigh-Taylor不稳定问题。如图5、图6所示，此问题考虑的是两种轻重流体在重力或惯性力的作用下，交界面不稳定，从而发生湍流混合，产生非常复杂的流场结构。考虑带重力源项的方程(1)，假设重力方向向上，取计算区域为[0，0.25]×[0，1.0]。初始界面位于y＝0.5处，重流体位于界面以下，初始状态为：ρ＝2.0，u＝0.0，v＝-0.025C·cos(8πx)，p＝2y+1；轻流体位于界面以上，初始状态为：ρ＝1.0，u＝0.0，v＝-0.025C·cos(8πx)，其中C为音速，γ＝5/3。左右边界为反射边界条件，上下边界设为无反射边界条件，上边界参数设为ρ＝1.0，u＝0.0，v＝0.0，p＝2.5，下边界参数设为ρ＝2.0，u＝0.0，v＝0.0，p＝1.0。计算终止时刻为t＝1.95。

以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。

Claims (7)

1.一种基于非等距网格下的有限体积流场数值计算方法，其特征在于，在非等距不均匀笛卡尔网格中，构造5阶精度的基本加权无振荡格式进行流场内的数值计算模拟，具体包括以下步骤：

步骤1、在流场数值计算区域内，生成非等距笛卡尔坐标系；

步骤2、在建立的非等距笛卡尔坐标系中构造有限体积WENO格式，其线性权可以任意取值且与网格尺度和类型无关，并且能保持格式的五阶精度，随后结合Runge-Kutta时间离散方法进而求解无粘流体的控制方程组，得到流场内所有非等距网格单元点的物理量的高精度数值近似值。

2.如权利要求1所述的一种基于非等距网格下的有限体积流场数值计算方法，其特征在于：所述步骤1中，非等距笛卡尔坐标系具体为：

随机网格：

Δx_i＝0.5Δx+(1.5Δx-0.5Δx)·random_number(i)，

Δy_j，＝0.5Δy+(1.5Δy-0.5Δy)·random_number(j)，

其中，下标t，j分别表示x和y方向的网格序号，Δx_i表示x方向上第i个网格对应的步长，Δy_j表示y方向上第j个网格对应的步长，Δx表示x方向上的平均步长，Δy表示y方向上的平均步长，random_number(i)和random_number(j)分别表示第i个和第j个网格对应的0到1之间的随机数；

或周期变化网格：

其中，L_x，L_y分别表示计算区域x方向和y方向的长度，N_x，N_y分别表示计算区域x方向和y方向的网格单元数；

或伸缩网格：

Δx_i＝Δx_i-1·α，

Δy_j＝Δy_j-1·β，

其中，α和β分别表示为x方向和y方向的伸缩比。

3.如权利要求2所述的一种基于非等距网格下的有限体积流场数值计算方法，其特征在于：所述步骤2具体包括：

2.1对无粘流体的控制方程组两边同时在目标单元内积分，构造有限体积WENO格式中的空间导数项，进而得到半离散的有限体积格式；

2.2利用Runge-Kutta时间离散方法求解半离散有限体积格式，从而得到控制方程组的时空全离散有限体积格式，其中的一组非线性权可以任意取值，与网格尺度和类型无关，只需满足和为1的正数即可；

2.3时空全离散有限体积格式具体为时间方向上的一个迭代公式，给定流场初始状态的相关物理量，根据迭代公式不断循环，即可得到流场某一时刻或稳定状态时的物理量的高精度近似值。

4.如权利要求3所述的一种基于非等距网格下的有限体积流场数值计算方法，其特征在于：所述步骤2中，无粘流体的控制方程组为：

其中，t表示时间变量，x，y表示空间变量，U＝(ρ，ρu，ρv，E)^T表示守恒变量，f(U)＝(ρu，ρu²+p，ρuv，u(E+p))^T，g(U)＝(ρv，ρuv，ρv²+p，v(E+p))^T，f(U)，g(U)表示通量，f(U)_x表示f(U)对x求导，g(U)_y表示g(U)对y求导，ρ，p，u，v，E分别表示流体密度、压强、水平方向速度、竖直方向速度以及能量，T表示转置，U₀表示初始状态值。

5.如权利要求4所述的一种基于非等距网格下的有限体积流场数值计算方法，其特征在于：所述步骤2.1具体包括：

对控制方程两边同时在目标网格单元I_i，j内积分有：

其中，表示守恒变量U(x，y，t)在目标单元I_i，j内的平均值，表示对时间变量t求导，表示对应的被积函数Φ在对应积分区间[a，b]内积分；

构造有限体积WENO格式中的空间导数项，具体过程如下：

2.1.1利用三点Gauss求积公式求解上式右端中的积分项，即

其中，G_m表示三点Gauss求积公式对应的求积节点；

2.1.2利用Lax-Friedrichs数值通量近似代替Gauss求积公式中每个求积节点的流通量，即

其中，表示对应点右侧高阶近似值，表示对应点左侧高阶近似值，表示对应点上侧高阶近似值，表示对应点上侧高阶近似值，α表示通量f(U)或g(U)对守恒变量U求导的最大值；

2.1.3在非等距网格下高阶近似的值，其它点的高阶近似值类似求得，具体如下：

2.1.3.1固定空间变量y，即假设y＝y_k，利用网格平均值得到的高阶近似值，过程如下：

a)选择一个母模板进而得到一个四次重构多项式p₁(x，y_k)，其满足以下条件：

得到其具体表达式为：

其系数满足线性方程组U＝C×A，具体形式如下：

A＝(a₁，a₂，a₃，a₄，a₅)^T，

其中U，A表示向量，C表示系数矩阵，T表示对向量或矩阵转置；

b)选择两个小模板从而得到两个线性重构多项式p₂(x，y_k)，p₃(x，y_k)，其分别满足：

其具体表达式为：

c)任意选择三个小于1且和为1的正数作为线性权γ_n(n＝1，2，3)；

d)计算每个模板对应的光滑指示器β_n，k(n＝1，2，3)，其计算公式如下：

其中，n表示对应多项式的序号，α表示求和指标，r表示对应多项式的次数，表示多项式p_n(x，y_k)对自变量x求α次偏导数；

对应光滑指示器的具体表达式为：

e)根据线性权γ_n和光滑指示器β_n，k求非线性权ω_n，k，公式如下：

其中，为计算过程中的过渡值；

f)得到的重构近似表达式：

2.1.3.2分别令y＝y_k-2，y_k-1，y_k+1，y_k+2，重复步骤2.1.3.1，从而得到

的值；

2.1.3.3固定空间变量x，即分别令x＝x_i+1/2和x＝x_i-1/2，利用步骤2.1.3.1和步骤2.1.3.2得到的值，得到的最终高阶近似表达式：

2.1.3.4重复步骤2.1.3.1到步骤2.1.3.3，类似得到的值；

2.1.4将步骤2.1.3的值代入步骤2.1.2，再代入步骤2.1.1，即可得到流体控制方程组的半离散有限体积格式，并表示为：

6.如权利要求5所述的一种基于非等距网格下的有限体积流场数值计算方法，其特征在于：所述步骤2.2中，假设时间步长为Δt，tⁿ表示第n时间层，表示的值，在tⁿ时间层进行步骤2.1的过程即可得到然后利用满足TVD性质的Runge-Kutta时间离散方法离散半离散有限体积格式从而得到时空全离散有限体积格式，如下：

其中，和为计算中间的过渡值。

7.如权利要求6所述的一种基于非等距网格下的有限体积流场数值计算方法，其特征在于：所述步骤2.3中，空间方向和时间方向的离散完成后得到的时空全离散有限体积格式为关于时间层的迭代公式，已知流场初始状态的相关物理量的值，根据迭代公式不断循环，即可得到流场某一时刻或稳定状态时的物理量的高精度近似值。