CN111177965B - 基于定常问题求解的多重分辨weno格式定点快速扫描方法 - Google Patents

Abstract

基于定常问题求解的多重分辨WENO格式定点快速扫描方法，应用于模拟定常可压流场问题，是根据双曲型方程沿特征线方向传播的性质加上Gauss‑Seidel迭代构造出的一种时间离散方法，该方法可以明显提高迭代法的收敛速度。相对于传统的Runge‑Kutta时间离散方法，可以节约迭代时间一半左右。该方法不止能有效应用于矩形规则计算区域，在复杂的不规则计算区域也得到较好的计算结果，可使大量的定常问题的残差快速下降到机器零。本发明不但保持了多分辨加权基本无振荡格式的优点，还可以明显加快残差收敛到机器零的速度。最后，本发明采用该格式对多个经典的定常问题进行了数值模拟，充分验证了本发明的有效性和可靠性。

Description

基于定常问题求解的多重分辨WENO格式定点快速扫描方法

技术领域

本发明属于计算流体力学工程技术领域，具体涉及基于定常问题求解的多重分辨WENO格式定点快速扫描方法。

背景技术

双曲守恒律方程的稳态问题是流体力学领域中常见的数学问题，在工程应用中十分重要。因此，构造解决这类问题的高鲁棒性和高精度数值模拟方法至关重要，也吸引了很多研究者的兴趣。在计算大型定常问题时，虽然大型超级计算机的广泛应用可以有效节约时间，但在算法设计中程序运行效率依然至关重要。

1959年，Godunov为解流场问题提出了一阶精度的数值格式。一阶精度的数值方法在捕捉激波时能抑制非物理振荡但会把强间断过度抹平，而强间断对问题的后续研究有着重要意义，因此需构造高精度数值格式来准确捕捉强间断。为了提高格式的计算精度，模拟解的光滑结构以及准确地捕捉激波位置，Harten于1983年首次提出了全变差减少(TVD)格式，并在此基础上与Osher于1987年提出了基本无振荡格式(ENO)格式。ENO格式的主要思想是在逐次扩展的模板中选用最光滑的模板构造多项式求出单元边界处的值，进而在光滑区域达到高阶精度，同时在间断附近实现基本无振荡的效果。但在格式的构造过程中，ENO格式是选用所有候选模板中的最优模板，其他的模板全部浪费掉，且数值精度越高浪费得越多，严重影响了计算效率。为了提高模板信息的使用，Liu，Osher和Chan等于1994年提出了加权基本无振荡(WENO)格式，提高了计算效率和计算精度。1996年，Jiang和Shu进一步改善了WENO格式，使得数值精度能够提高到2r-1阶，并设计出了光滑因子和非线性权的构造框架。WENO格式的主要思想是通过低阶重构通量的线性凸组合获得高阶近似。但该经典WENO格式的实现过程中，线性权的计算复杂，且在很多定常问题中残差降不下去。因此，2018年Zhu和Shu提出了多重分辨WENO格式，线性权可以任取为和为1的正数，且在光滑区域格式的数值精度保持最优，可使许多经典的定常问题算例残差可以下降到接近机器零。但对于经典的WENO格式进行空间离散和TVD Runge-Kutta格式进行时间离散，迭代次数较多和迭代的CPU时间较长，迭代效率不够高。于是为了提高迭代效率，在时间离散上提出了新的离散方法，比如快速行进算法和快速扫描算法。快速行进算法是让解严格按照递增(递减)的顺序更新点值，因此需要两个基本的要素，迎风性和排序算法。快速行进算法的时间复杂度为O(NlogN)阶，其中N为网格点，logN部分来自排序算法。为了继续加快迭代效率，又提出了快速扫描算法。与快速行进法相比，快速扫描法是一种并行的算法。根据特征方向扫描方向分为有限数量的组，每个组按特定扫描方向进行Gauss-Seidel迭代，这样可以覆盖每个特征方向，因此不需要进行排序算法。快速扫描法最早被用来求解定常的Hamilton-Jacobi方程，其时间复杂度为O(N)。但时间复杂程度的系数与方程是息息相关的。2016年Wu和Zhang将该快速扫描算法应用到求解双曲守恒律方程，也可以明显地加快格式的迭代速度。

发明内容

本发明针对现有技术中的不足，提供一种基于定常问题求解的多重分辨WENO格式定点快速扫描方法，能针对各种可压定常流场问题进行高精度数值模拟。本发明给出了该方法的具体的构造过程。相比于经典WENO空间离散和TVD Runge-kutta时间离散格式，可以明显提高迭代效率，节约一半左右的迭代CPU时间。该快速扫描法是通过在所有方向上依次扫描，总有一次可以与计算节点的迎风方向相同，再结合Gauss-Seidel迭代的特点，有新值用新值没有新值用旧值，可以提高欧拉向前时间离散格式的迭代效率，在保证格式收敛的情况下CFL数可以取较大值。Runge-kutta时间离散格式虽然也可以取到很大的CFL数，但多计算了两个虚拟时间层，导致了计算效率的低下。对于空间离散，经典五阶有限差分WENO格式的残差无法下降到接近机器零，而多重分辨WENO格式能较好地使残差快速下降到机器零附近。此格式的构造极为简单，数值精度较高，易于推广到多维情形。

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

基于定常问题求解的多重分辨WENO格式定点快速扫描方法，其特征在于，包括如下步骤：

步骤1：把双曲守恒律方程离散为空间半离散有限差分格式，用多重分辨WENO格式重构通量的近似值；

步骤2：将重构通量代入含有时间导数项的半离散有限差分格式，得到关于时间导数的常微分方程，对方程中的时间导数使用快速扫描公式将半离散有限差分格式离散成时空全离散有限差分格式；

步骤3：根据时空全离散有限差分格式得到下一时间层上的近似值，依次迭代，得到计算区域残差稳定时流场的数值结果。

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

进一步地，所述步骤1具体如下：

考虑一维双曲守恒律方程：

其半离散格式的形式为：

其中，U＝(ρ,ρu,E)^T表示守恒变量，f(U)＝(ρu,ρu²+p,u(E+p))^T表示通量，U_t表示U对t求导，f_x(U)表示f(U)对x求导，ρ，u，p，E分别表示流体密度，速度，压强，能量，T表示转置，U₀表示初始状态值，L(U)表示-f_x(U)的空间离散形式，t表示时间变量，x表示空间变量；

把空间离散成统一长度的网格单元

单元长度

单元中心为

其中i为坐标序号，有

其中，

和

分别表示通量f(U)在目标单元I_i的边界

和

处的五阶近似的数值通量，U_i(t)表示U在网格单元I_i内点x_i处的值U(x_i,t)。

进一步地，求通量f(U)在目标单元I_i的边界

和

处的五阶近似值

和

具体如下：

1)用Lax-Friedrichs分裂把通量分裂为

其中

2)将目标单元I_i以及其周围共5个单元组成一个大模板T₃＝[I_i-2,I_i-1,I_i,I_i+1,I_i+2]，从大模板中选择一个包含一个单元的小模板T₁＝[I_i]和一个包含三个单元的小模板T₂＝[I_i-1,I_i,I_i+1]，其中I_i为对应序号的网格单元；

3)在每个模板上分别重构代数多项式q₁(x)、q₂(x)和q₃(x)，使得它们在单元边界有五阶精度，其具体过程如下：在三个模板T₁、T₂和T₃上分别构造代数多项式q₁(x)，q₂(x)和q₃(x)，使其满足：

4)任意取线性权为：γ₁₂＝1/11，γ₂₂＝10/11，γ₁₃＝1/111，γ₂₃＝10/111，γ₃₃＝100/111；重新构造出p₁(x)，p₂(x)和p₃(x)，满足：

p₁(x)＝q₁(x) (7)

5)计算光滑指示器β_l，用于衡量重构多项式p_l(x)在目标单元上的光滑度，计算公式为：

其中l＝2,3表示对应模板序号，

表示多项式p_l(x)对x的α阶导数，r＝2；β₁的独立定义为：

其中，

γ_1,1＝1-γ_0,1，

6)通过线性权γ_l和光滑指示器β_l计算非线性权ω_l，其计算公式为：

其中l＝1,2,3表示对应模板序号，

τ为计算过程中的过渡值，β_l为光滑指示器，ε＝10^-6防止分母为零；

7)求出数值通量分裂f⁺(U)在点

处的多重分辨WENO重构值：

相应地，求出数值通量分裂f^-(U)在点

处的多重分辨WENO重构值

进一步地，步骤2中，时间离散的定点快速扫描公式为：

其中*代表有n+1层的新值就用n+1层的新值，没有新值就用第n层的旧值；而扫描的顺序为：i＝i₁,…,i_N＝1,…,N和i＝i₁,…,i_N＝N,…,1，交替扫描，N表示网格数量；

使用上述快速扫描公式将半离散有限差分格式离散成时空全离散有限差分格式。

进一步地，步骤3中，时空全离散有限差分格式为关于时间层的迭代公式，初始状态值已知，通过迭代公式求出下一时间层的近似值，依次得到残差稳定时计算区域内的数值模拟值；残差ResA的定义如下：

其中，

Δt_n表示第n层的时间步长。

本发明的有益效果是：相比于经典WENO格式，该多重分辨WENO格式通过采用一系列不等距中心模板的信息使定常问题的残差下降得更快且其值能接近于机器零。该格式能精确捕捉激波，且在解的光滑区域能保持最优数值精度。能任意选择线性权的取值，在减少计算量的同时不降低格式在解的光滑区域的的数值精度。相比于经典的TVD Runge-Kutta时间离散和欧拉向前时间离散，快速扫描法可以取较大的CFL数，可极大减少格式的迭代次数和节省大量的CPU时间，并且该方法易于推广到高维情形。

附图说明

图1a～1c是例一的残差下降图，空间离散用多重分辨WENO离散，时间离散分别用欧拉向前离散、定点快速扫描离散和三阶龙格库塔离散，CFL数分别用0.1、1.0和1.0。图中不同的线是测试的不同网格时的残差下降曲线。

图2a～2c是例二的残差下降图，空间离散用多重分辨WENO离散，时间离散分别用欧拉向前离散、定点快速扫描离散和三阶龙格库塔离散，CFL数分别用0.1、1.0和1.0。图中不同的线是测试的不同网格时的残差下降曲线。

图3a～3c是例三的残差下降图，空间离散用多重分辨WENO离散，时间离散分别用欧拉向前离散、定点快速扫描离散和三阶龙格库塔离散。图中不同的曲线是测试取不同的CFL数时的残差下降曲线。

图4a～4b是例四中的采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝0.5时的密度等值线和不同CFL数时残差的下降曲线。

图4c～4d是例四中的采用多重分辨WENO空间离散，三阶龙格库塔空间离散，得到的CFL＝0.5时的密度等值线和不同CFL数时残差的下降曲线。

图5a～5b是例五中的采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝0.6时的密度等值线和不同CFL数时残差的下降曲线。

图5c～5d是例五中的采用多重分辨WENO空间离散，三阶龙格库塔空间离散，得到的CFL＝0.6时的密度等值线和不同CFL数时残差的下降曲线。

图6a～6b是例六中的采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝1.4时的密度等值线和不同CFL数时残差的下降曲线。

图6c～6d是例六中的采用多重分辨WENO空间离散，三阶龙格库塔空间离散，得到的CFL＝1.2时的密度等值线和不同CFL数时残差的下降曲线。

图7a～7b是例七中的采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝1.3时的密度等值线和不同CFL数时残差的下降曲线。

图7c～7d是例七中的采用多重分辨WENO空间离散，三阶龙格库塔空间离散，得到的CFL＝1.2时的密度等值线和不同CFL数时残差的下降曲线。

图8a～8b是例八中的采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝1.1时的密度等值线和不同CFL数时残差的下降曲线。

图8c～8d是例八中的采用多重分辨WENO空间离散，三阶龙格库塔空间离散，得到的CFL＝1.2时的密度等值线和不同CFL数时残差的下降曲线。

图9a～9b是例九中的采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝1.3时的密度等值线和不同CFL数时残差的下降曲线。

图9c～9d是例九中的采用多重分辨WENO空间离散，三阶龙格库塔空间离散，得到的CFL＝1.4时的密度等值线和不同CFL数时残差的下降曲线。

图10a～10b是例十中的采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝0.9时的密度等值线和不同CFL数时残差的下降曲线。

图10c～10d是例十中的采用多重分辨WENO空间离散，三阶龙格库塔空间离散，得到的CFL＝0.9时的密度等值线和不同CFL数时残差的下降曲线。

具体实施方式

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

基于定常问题求解的多重分辨WENO格式的快速扫描算法的构造。

考虑一维双曲守恒律方程：

其半离散格式的形式为：

其中，U＝(ρ,ρu,E)^T表示守恒变量，f(U)＝(ρu,ρu²+p,u(E+p))^T表示通量，U_t表示U对t求导，f(U)_x表示f(U)对x求导，ρ，u，p，E分别表示流体密度，速度，压强，能量，T表示转置，U₀表示初始状态值，L(U)表示-f_x(U)的空间离散形式，t表示时间变量，x表示空间变量。

把空间离散成统一长度的网格单元

单元长度

单元中心为

其中i为坐标序号，有：

其中，

和

分别表示通量f(U)在目标单元I_i的边界

和

处的五阶近似的数值通量，U_i(t)表示U在网格单元I_i内点x_i处的值U(x_i,t)。

求通量f(U)在目标单元I_i的边界

和

处的五阶近似值

和

具体步骤如下：

1)用最简单的Lax-Friedrichs分裂把通量分裂为

其中

为了简单起见，本发明只描述f⁺(U)在点

处的重构过程并将其定义为

2)将目标单元I_i以及其周围共5个单元组成一个大模板T₃＝[I_i-2,I_i-1,I_i,I_i+1,I_i+2]，从大模板中选择一个包含一个单元的小模板T₁＝[I_i]和一个包含三个单元的小模板T₂＝[I_i-1,I_i,I_i+1]，其中I_i为对应序号的网格单元。

3)在每个模板上分别重构代数多项式q₁(x)、q₂(x)和q₃(x)，使得它们在单元边界有五阶精度。其具体过程如下：在三个模板T₁、T₂和T₃上分别构造代数多项式q₁(x)，q₂(x)和q₃(x)，使其满足:

4)任意取线性权为：γ₁₂＝1/11，γ₂₂＝10/11，γ₁₃＝1/111，γ₂₃＝10/111，γ₃₃＝100/111。重新构造出p₁(x)，p₂(x)和p₃(x)，满足：

p₁(x)＝q₁(x) (7)

5)计算光滑指示器β_l，用于衡量重构多项式p_l(x)在目标单元上的光滑度，计算公式为：

其中l＝2,3表示对应模板序号，

表示多项式p_l(x)对^x的α阶导数，r＝2。但是β₁比较特殊，要独立定义：

其中，

γ_1,1＝1-γ_0,1，

6)通过线性权γ_l和光滑指示器β_l计算非线性权ω_l，其计算公式为：

其中l＝1,2,3表示对应模板序号，

τ为计算过程中的过渡值，β_l为光滑指示器，ε＝10^-6防止分母为零。

7)求出数值通量分裂f⁺(U)在点

处的多重分辨WENO重构值：

其次，将计算结果代入含有时间导数项的半离散有限差分格式，得到关于时间导数的常微分方程。

最后，给出了三种时间离散方法：

1)时间离散的三阶TVD Runge-Kutta时间离散公式：

其中，U⁽¹⁾,U⁽²⁾为中间过渡值，Δt为时间步长，上标n表示第n时间层，L(Uⁿ)、L(U⁽¹⁾)、L(U⁽²⁾)为-f_x(U)的高阶空间离散形式的近似值。

2)时间离散的欧拉向前格式为：

3)时间离散的定点快速扫描法格式为：

其中*代表有n+1层的新值就用n+1层的新值，没有新值就用第n层的旧值。而扫描的顺序为：i＝1,…,N和i＝N,…,1，交替扫描。

这样就得到时空全离散有限差分格式，时空全离散有限差分格式为关于时间层的迭代公式，初始状态值已知，通过迭代公式求出下一时间层的近似值，依次得到残差稳定时计算区域内的数值模拟值。对于二维问题，逐维用上面的重构过程。残差ResA的定义如下：

其中，

对于我们的例子欧拉方程，

其中

Δt_n表示第n层的时间步长。

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

例一：带源项的二维欧拉方程：

迭代初值为：ρ(x,y,0)＝1+0.2sin(x+y)，u(x,y,0)＝1，v(x,y,0)＝1，p(x,y,0)＝1+0.2sin(x+y)。计算区域为：(x,y)∈[0,2π]×[0,2π]。其数值结果如图1a～1c所示。

例二：二维欧拉方程：

它的迭代初值为：ρ(x,y,0)＝1+0.2sin(x-y),u(x,y,0)＝1,v(x,y,0)＝1,p(x,y,0)＝1。

计算区域为：x∈[-1,1]。其数值结果如图2a～2c所示。

例三.二维欧拉方程：

也可以写成：U_t+f(U)_x+g(U)_y＝0，迭代初值为

其中：

马赫数M＝2，计算区域为x∈[-1,1]，其数值结果如图3a～3c所示。

例四：经典的二维稳定斜激波算例。计算区域为0≤x≤4，0≤y≤2，激波与x轴有135度角。其数值结果如图4a～4d所示。

例五：经典的正规激波反射算例。下边界为反射边界，左边界和上边界为狄利克雷边界，右边界为超音速出口，计算区域为0≤x≤4，0≤y≤1，其数值结果如图5a～5d所示。

例六：超音速流通过短板算例。超音速流以与短板α＝10度角经过。马赫数M＝3，γ＝1.4。初值为：

ρ＝1，u＝cos(α)，v＝sin(α)，0≤x≤10，-5≤y≤5。短板位于x∈[1,2]，y＝0，其数值结果如图6a～6d所示。

例七：超音速流通过两个短板。超音速流以与短板α＝10度角经过。马赫数M＝3，γ＝1.4。初值为：

ρ＝1，u＝cos(α)，v＝sin(α)，0≤x≤10，-5≤y≤5。短板位于x∈[1,2]，y＝-2和y＝2，其数值结果如图7a～7d所示。

例八：超音速流通过三个短板。超音速流以与短板α＝10度角经过。马赫数M＝3，γ＝1.4。初值为：

ρ＝1，u＝cos(α)，v＝sin(α)，0≤x≤10，-5≤y≤5。短板位于x∈[1,2],y＝0和x∈[2,3],y＝-2和x∈[2,3],y＝2。其数值结果如图8a～8d所示。

例九：超音速流通过一个长板。超音速流以与短板α＝10度角经过。马赫数M＝3，γ＝1.4。初值为：

ρ＝1，u＝cos(α)，v＝sin(α)，0≤x≤7，-5≤y≤5。短板位于x∈[2,7],y＝0，其数值结果如图9a～9d所示。

例十：超音速流通过三个长板。超音速流以与短板α＝10度角经过。马赫数M＝3，γ＝1.4。初值为：

ρ＝1，u＝cos(α)，v＝sin(α)，0≤x≤5，-5≤y≤5。短板位于x∈[2,5]，y＝0，y＝-2和y＝2。其数值结果如图10a～10d所示。

需要注意的是，发明中所引用的如“上”、“下”、“左”、“右”、“前”、“后”等的用语，亦仅为便于叙述的明了，而非用以限定本发明可实施的范围，其相对关系的改变或调整，在无实质变更技术内容下，当亦视为本发明可实施的范畴。

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

Claims (2)

1.基于定常问题求解的多重分辨WENO格式定点快速扫描方法，用于捕捉激波位置，其特征在于，包括如下步骤：

步骤1：把双曲守恒律方程离散为空间半离散有限差分格式，用多重分辨WENO格式重构通量的近似值；所述步骤1具体如下：

考虑一维双曲守恒律方程：

其半离散格式的形式为：

其中，U＝(ρ,ρu,E)^T表示守恒变量，f(U)＝(ρu,ρu²+p,u(E+p))^T表示通量，U_t表示U对t求导，f_x(U)表示f(U)对x求导，ρ，u，p，E分别表示流体密度，速度，压强，能量，T表示转置，U₀表示初始状态值，L(U)表示-f_x(U)的空间离散形式，t表示时间变量，x表示空间变量；

把空间离散成统一长度的网格单元

单元长度

单元中心为

其中i为坐标序号，有

其中，

和

分别表示通量f(U)在目标单元I_i的边界

和

处的五阶近似的数值通量，U_i(t)表示U在网格单元I_i内点x_i处的值U(x_i,t)；

步骤2：将重构通量代入含有时间导数项的半离散有限差分格式，得到关于时间导数的常微分方程，对方程中的时间导数使用快速扫描公式将半离散有限差分格式离散成时空全离散有限差分格式；所述步骤2中，时间离散的定点快速扫描公式为：

其中*代表有n+1层的新值就用n+1层的新值，没有新值就用第n层的旧值；而扫描的顺序为：i＝i₁,…,i_N＝1,…,N和i＝i₁,…,i_N＝N,…,1，交替扫描，N表示网格数量；

使用上述快速扫描公式将半离散有限差分格式离散成时空全离散有限差分格式；

步骤3：根据时空全离散有限差分格式得到下一时间层上的近似值，依次迭代，得到计算区域残差稳定时流场的数值结果。

2.如权利要求1所述的基于定常问题求解的多重分辨WENO格式定点快速扫描方法，其特征在于：求通量f(U)在目标单元I_i的边界

和

处的五阶近似值

和

具体如下：

1)用Lax-Friedrichs分裂把通量分裂为

其中

2)将目标单元I_i以及其周围共5个单元组成一个大模板T₃＝[I_i-2,I_i-1,I_i,I_i+1,I_i+2]，从大模板中选择一个包含一个单元的小模板T₁＝[I_i]和一个包含三个单元的小模板T₂＝[I_i-1,I_i,I_i+1]，其中I_i为对应序号的网格单元；

3)在每个模板上分别重构代数多项式q₁(x)、q₂(x)和q₃(x)，使得它们在单元边界有五阶精度，其具体过程如下：在三个模板T₁、T₂和T₃上分别构造代数多项式q₁(x)，q₂(x)和q₃(x)，使其满足：

4)任意取线性权为：γ₁₂＝1/11，γ₂₂＝10/11，γ₁₃＝1/111，γ₂₃＝10/111，γ₃₃＝100/111；重新构造出p₁(x)，p₂(x)和p₃(x)，满足：

p₁(x)＝q₁(x) (7)

5)计算光滑指示器β_l，用于衡量重构多项式p_l(x)在目标单元上的光滑度，计算公式为：

其中l＝2,3表示对应模板序号，

表示多项式p_l(x)对x的α阶导数，r＝2；β₁的独立定义为：

其中，

γ_1,1＝1-γ_0,1，

6)通过线性权γ_l和光滑指示器β_l计算非线性权ω_l，其计算公式为：

其中l＝1,2,3表示对应模板序号，

τ为计算过程中的过渡值，β_l为光滑指示器，ε＝10^-6防止分母为零；

7)求出数值通量分裂f⁺(U)在点

处的多重分辨WENO重构值：

相应地，求出数值通量分裂f^-(U)在点

处的多重分辨WENO重构值

Images