一种多重分辨率WENO格式结合ILW边界处理的定点快速扫描方法

摘要

本发明公开了一种多重分辨率WENO格式结合ILW边界处理的定点快速扫描方法，包括：把定常双曲守恒律问题转化为依赖时间的双曲守恒律问题，用新型的五阶ILW边界处理方法处理物面边界；将双曲守恒律方程空间部分用有限差分多重分辨加权基本无振荡格式离散；对控制方程中的时间部分用三阶龙格库塔方法和定点快速扫描法离散；根据时空全离散方法得到下一时间层每一个点的近似值，得到计算区域内守恒变量的残差在趋于稳定时的数值模拟结果。该方法能保持物面边界附近计算格式高阶精度不降低的优点，有效减少空间模板的使用数量和减少需要人为调整的参数，实施更加简单且易于实现。

说明书

技术领域

本发明属于计算流体力学工程技术领域，具体涉及一种多重分辨率WENO格式结合ILW边界处理的定点快速扫描方法。

背景技术

双曲守恒律方程的稳态问题是流体力学领域中经常遇到的数学问题，在工程应用中也占据重要的位置。因此，构造解决此类问题的高鲁棒性和高精度数值模拟方法变得很重要，很多学者也在这一领域攻坚克难。由于实际的问题当中计算区域多种多样，导致其物理量残差很难收敛到机器零。因此设计一种可以使残差下降到机器零的物面边界的高精度算法显得格外重要。在计算大型定常问题时，虽然大型超级计算机的广泛应用可以使计算对时间没那么敏感，但在算法设计中，程序运行的效率依然非常重要。

为捕捉流场内的激波，很多激波捕捉格式被设计出来。Godunov最先提出了一阶精度的激波捕捉格式。一阶精度的数值方法能抑制非物理振荡但会把大梯度处过度抹平，而强间断对物理现象的描述十分重要，因此需构造高阶精度的数值格式来更精确捕捉强间断。Harten于1983年首次提出了全变差减少(TVD)格式，之后Osher接着提出了本质无振荡格式(ENO)格式。ENO格式的主要思想是在逐次扩展的模板中选用最光滑的模板构造多项式求出单元半点的值，进而在光滑区域达到高阶精度，同时在间断附近实现基本无振荡的效果。但通过格式的构造过程可以看到，ENO格式是选用所有候选模板中的最优模板，其他的模板全部浪费掉，且数值精度越高浪费得越多，严重影响了计算效率。为了提高模板更有效的使用，Liu，Osher和Chan等于1994年提出了加权本质无振荡(WENO)格式，提高了一阶的计算精度。1996年，Jiang和Shu进一步改善了WENO格式，使得数值精度能够提高到2r-1阶，并设计出了光滑因子和非线性权的构造框架。WENO格式的主要思想是通过低阶重构通量的线性凸组合获得高阶近似。但该经典WENO格式的实现过程中，线性权的计算复杂，且在很多定常问题中残差无法下降到机器零。因此，2018年Zhu和Shu提出了多重分辨WENO格式，通过设计了不等距模板，让线性权可以任取为和为1的正数，且在光滑区域格式的数值精度保持最优，可使许多经典的定常问题算例残差可以下降到接近机器零。

但是对于复杂的计算区域的数值模拟，有限差分方法最大的困难就是物面边界的处理。目前非贴体网格的主要方法有浸入边界法、嵌入边界法、SharpInterface方法、虚拟单元法。但这些方法都有各自的局限性，而且有一个共同的缺点，数值精度达不到高阶。于是2010年Sirui Tan提出了ILW边界处理方法，通过泰勒展开构造虚拟流场点上的值，耦合高阶WENO外推成功模拟了静止和运动物体的流动问题，数值结果显示对于静止问题可以达到五阶精度，对于运动问题可以达到三阶精度。之后ILW边界处理方法又被推广到了求带源项的欧拉方程，使得求解带粘性的NS方程也成为可能。后来又被改进为了简化的ILW边界处理方法，即高阶导数可直接由多项式外推得到，不用ILW过程依然不影响精度。

对于欧拉方程经典的三阶TVD Runge-Kutta时间离散，迭代次数较多和迭代的CPU时间较长，迭代效率有些低。为了提高迭代效率，于是在时间离散上提出了很多新的离散方法，比如快速行进算法和快速扫描算法。快速行进算法是让解严格按照递增或递减的顺序更新点值，这又会增加排序的时间复杂度。为了继续加快迭代效率就提出了快速扫描算法。与快速行进法相比，快速扫描法是一种并行的算法。由迎风性扫描覆盖了每一个特征方向，这样就省去了排序的时间，再结合Gauss-Seidel迭代，可以很大程度加快迭代的收敛。快速扫描法最开始使被用来求解静态的Hamilton-Jacobi方程。2016年Wu和Zhang将该快速扫描算法应用到求解双曲守恒律方程，也可以明显地加快格式的迭代速度。2020年Li和Zhu又将该算法与多重分辨WENO结合求解简单计算区域的定常流问题，达到了快速的完全收敛。

发明内容

本发明所要解决的技术问题是针对上述现有技术的不足，提供一种多重分辨率WENO格式结合ILW边界处理的定点快速扫描方法。

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

一种多重分辨率WENO格式结合ILW边界处理的定点快速扫描方法，其特征在于，所述方法用于针对多种复杂计算区域的可压定常流场问题进行高精度数值模拟，包括：

步骤1.在笛卡尔坐标系下，把定常双曲守恒律问题转化为依赖时间的双曲守恒律问题，用新型的五阶ILW边界处理方法处理物面边界；

步骤2.将双曲守恒律方程空间部分用有限差分多重分辨加权基本无振荡格式进行离散；

步骤3.对控制方程中的时间部分用三阶龙格库塔方法和定点快速扫描法离散成全离散的有限差分格式；

步骤4.根据时空全离散方法得到下一时间层每一个点的近似值，依次迭代，得到计算区域内守恒变量的残差在趋于稳定时的数值模拟结果。

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

上述的步骤1中，对于一维双曲守恒律方程：

其半离散格式的形式为：

其中，U＝(ρ,ρu,E)

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

上述的步骤1中，用新型的五阶ILW边界处理方法处理物面边界，物面边界外的虚拟点用ILW边界处理,虚拟点由物面边界点经泰勒展开求出，包括：

步骤1.1虚拟点用泰勒展开公式求出，U的第m个分量在点x

然后通过已知的边界条件确定该边界是出流还是入流，如果U

入流边界用时间导数转换成空间导数求，出流边界用多项式外推求得；

记L(U

记(V

3,4.用g

步骤1.2通过原方程时间和空间的转换得到方程：

通过(4)这个线性方程求出边界值的一阶导数，其余各阶导数也同样的方式求出；

然后代入(3)式求出物面边界附近虚拟点的值。

上述的V

步骤a：选取三个模板T

步骤b：计算光滑指示器β

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

步骤c：通过线性权γ

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

步骤d：通过非线性权和多项式函数p,可以得到V函数的各阶导数：

上述的步骤a具体过程如下：

在三个模板T

取线性权为:γ

重新构造出p

p

上述的步骤1中，对于二维双曲守恒律方程：

其半离散格式的形式为

其中，U＝(ρ,ρu,ρv,E)

计算过程中需要用到物面边界外的虚拟点的值，根据所求的虚拟点的不同变换方程，对于点U(i,j)，记该点为P点；

做P点与附近物面边界的外法线，外法线与物面边界交点记为P

原方程可变换成：

U

其中：

上述的步骤1中，用新型的五阶ILW边界处理方法处理物面边界，物面边界外的虚拟点用ILW边界处理,虚拟点由物面边界点经泰勒展开求出，包括：

步骤1.1点P的值可用泰勒展开公式求出，U的第m个分量在点P处的值为：

其中Δ为点P

步骤1.2通过原方程时间和空间的转换得到方程：

通过(29)这个线性方程求出边界值的一阶导数，其余各阶导数以同样的方式求出；

然后代入(28)式就可以求出物面边界附近虚拟点的值。

上述的步骤2中，一维情形时将方程写成守恒形式：

其中，

上述的步骤2中，求通量f(U)在目标单元I

步骤2.1用Lax-Friedrichs分裂把通量分裂为

步骤2.2将目标单元I

步骤2.3在每个模板上分别重构代数多项式q

步骤2.4任意取线性权为：γ

将多项式q

步骤2.5计算光滑指示器β

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

其中：

步骤2.6通过线性权γ

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

步骤2.7求出数值通量分裂f

从而得到空间离散的常微分方程。

上述的步骤3为：

使用三阶TVD Runge-Kutta时间离散公式和定点快速扫描公式将半离散有限差分格式离散成时空全离散有限差分格式，所述时空全离散有限差分格式为关于时间层的迭代公式；

所述步骤3中，时间离散的三阶TVD Runge-Kutta时间离散公式：

其中，U

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

其中*代表有n+1层的新值就用n+1层的新值，没有新值就用第n层的旧值，扫描的顺序为：

i＝1,…,N和i＝N,…,1交替扫描。

本发明具有以下有益效果：

相比于经典的ILW边界处理方法，新的WENO(Multi-resolution WeightedEssentially Non-Oscillatory，分多分辨加权基本无振荡格式)外推改进了线性权与网格尺度相关的问题，线性权可以任意取值。而且节省了使用的空间模板数量，减少了非线性权中参数的选取。使用起来更加方便简单，也有更广泛的适用性。

相比于经典WENO格式，该多重分辨WENO格式通过采用一系列不等距中心模版的信息使定常问题的残差下降得更快且其值能接近于机器零。该格式能精确捕捉激波，且在解的光滑区域能保持最优数值精度。能任意选择线性权的取值，在减少计算量的同时不降低格式在解的光滑区域的数值精度。

相比于经典的三阶Runge-Kutta时间离散，快速扫描法可以取较大的CFL数，可极大减少格式的迭代次数和节省大量的CPU时间，并且该方法易于推广到高维情形。

新型ILW(Inverse Lax-Wendroff)边界处理方法与原ILW边界处理方法相比，能在保持物面边界附近计算格式高阶精度不降低的优点，同时能有效减少空间模板的数量和减少需要人为调整的参数，在求解定常问题时物理量的残差可以快速收敛。此新型ILW边界处理方法在结合了定点快速扫描算法和多分辨加权基本无振荡格式后不但能保持原有高精度数值计算格式在激波捕捉时的各种优点，还可以明显加快物理量残差收敛到机器零的速度。最后，本发明采用上述方法对多个经典的二维定常可压流场问题进行了高精度数值模拟，充分验证了本发明的有效性和可靠性。

附图说明

图1为实施例一数值结果，其中图1a和图1c是实施例一中采用多重分辨WENO空间离散，龙格库塔时间离散，得到的CFL＝0.6时的压力等值线和残差的下降曲线；图1b和图1d是实施例一中采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝0.6时的压力等值线和残差的下降曲线。

图2为实施例二数值结果，其中图2a和图2c是实施例二中采用多重分辨WENO空间离散，龙格库塔时间离散，得到的CFL＝0.5时的压力等值线和残差的下降曲线；图2b和图2d是实施例二中采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝0.5时的压力等值线和残差的下降曲线。

图3为实施例三数值结果，其中图3a和图3c是实施例三中采用多重分辨WENO空间离散，龙格库塔时间离散，得到的CFL＝1.0时的压力等值线和残差的下降曲线；图3b和图3d是实施例三中采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝1.0时的压力等值线和残差的下降曲线。

图4为实施例四数值结果，其中图4a和图4c是实施例四中采用多重分辨WENO空间离散，龙格库塔时间离散，得到的CFL＝0.9时的压力等值线和残差的下降曲线；图4b和图4d是实施例四中采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝1.0时的压力等值线和残差的下降曲线。

图5为实施例五数值结果，其中图5a和图5c是实施例五中采用多重分辨WENO空间离散，龙格库塔时间离散，得到的CFL＝1.7时的压力等值线和残差的下降曲线；图5b和图5d是实施例五中采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝1.3时的压力等值线和残差的下降曲线。

图6为实施例六数值结果，其中图6a和图6c是实施例五中采用多重分辨WENO空间离散，龙格库塔时间离散，得到的CFL＝1.6时的压力等值线和残差的下降曲线；图6b和图6d是实施例五中采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝1.2时的压力等值线和残差的下降曲线。

图7为实施例七数值结果，其中图7a和图7c是实施例五中采用多重分辨WENO空间离散，龙格库塔时间离散，得到的CFL＝0.7时的压力等值线和残差的下降曲线；图7b和图7d是实施例五中采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝0.7时的压力等值线和残差的下降曲线。

图8是本发明方法流程图。

具体实施方式

以下结合附图对本发明的实施例作进一步详细描述。

参见图8，本发明的一种多重分辨率WENO格式结合ILW边界处理的定点快速扫描方法，用于针对多种复杂计算区域的可压定常流场问题进行高精度数值模拟，包括：

步骤1.在笛卡尔坐标系下，把定常双曲守恒律问题转化为依赖时间的双曲守恒律问题，用新型的五阶ILW边界处理方法处理物面边界；

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

其半离散格式的形式为

其中，U＝(ρ,ρu,E)

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

用新型的五阶ILW边界处理方法处理物面边界，物面边界外的虚拟点用ILW边界处理,虚拟点由物面边界点经泰勒展开求出，包括：

步骤1.1这些虚拟点用泰勒展开公式求出，U的第m个分量在点x

然后通过已知的边界条件确定该边界是出流还是入流，如果U

入流边界用时间导数转换成空间导数求，出流边界用多项式外推求得。记L(U

步骤1.2通过原方程时间和空间的转换可以得到方程：

通过(4)这个线性方程求出边界值的一阶导数，其余各阶导数也可以同样的方式求出。然后代入(3)式就可以求出物面边界附近虚拟点的值。

求

步骤a：选取三个模板T

其具体过程如下：

在三个模板T

任意取线性权为:γ

p

步骤b：计算光滑指示器β

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

步骤c：通过线性权γ

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

步骤d：通过非线性权和多项式函数p,可以得到V函数的各阶导数。

得到的各阶导数，代入式(3)中可以得到所需的所有虚拟点的点值。

再考虑二维双曲守恒律方程：

其半离散格式的形式为

其中，U＝(ρ,ρu,ρv,E)

计算过程中需要用到物面边界外的虚拟点的值，根据所求的虚拟点的不同变换方程。以求点U(i,j)为例，记该点为P点。做P点与附近物面边界的外法线，外法线与物面边界交点记为P

原方程可变换成：

U

其中：

步骤1中，用新型的五阶ILW边界处理方法处理物面边界，物面边界外的虚拟点用ILW边界处理,虚拟点由物面边界点经泰勒展开求出，包括：

步骤1.1点P的值可用泰勒展开公式求出，U的第m个分量在点P处的值为：

其中Δ为点P

步骤1.2通过原方程时间和空间的转换可以得到方程：

通过(29)这个线性方程求出边界值的一阶导数，其余各阶导数也可以同样的方式求出。然后代入(28)式就可以求出物面边界附近虚拟点的值。

求V

步骤a：选取三个模板T

任意取线性权为：γ

p

步骤b：计算光滑指示器β

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

步骤c：通过线性权γ

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

步骤d：通过非线性权和多项式函数p,得到V函数的各阶导数。

代入的式子(15)、(16)，可以得到所需的所有虚拟点的点值。

步骤2.将双曲守恒律方程空间部分用有限差分多重分辨加权基本无振荡格式进行离散；包括：

一维情形时将方程写成守恒形式

其中，

求通量f(U)在目标单元I

步骤2.1用Lax-Friedrichs分裂把通量分裂为

步骤2.2将目标单元I

步骤2.3在每个模板上分别重构代数多项式q

步骤2.4任意取线性权为：γ

步骤2.5计算光滑指示器β

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

其中：

步骤2.6通过线性权γ

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

步骤2.7求出数值通量分裂f

从而得到空间离散的常微分方程。对于二维和高维情况，逐维离散即可。

步骤3.对控制方程中的时间部分用三阶龙格库塔方法和定点快速扫描法离散成全离散的有限差分格式；

所述步骤3为：

使用三阶TVD Runge-Kutta时间离散公式和定点快速扫描公式将半离散有限差分格式离散成时空全离散有限差分格式，所述时空全离散有限差分格式为关于时间层的迭代公式。

所述步骤3，时间离散的三阶TVD Runge-Kutta时间离散公式：

其中，U

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

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

步骤4.根据时空全离散方法得到下一时间层每一个点的近似值，依次迭代，得到计算区域内守恒变量的残差在趋于稳定时的数值模拟结果。

所述步骤4为：

将所述时空全离散有限差分格式初始状态值设为精确解或者近似解，通过迭代公式求出下一时间层的近似值，依次迭代得到残差稳定时计算区域内守恒量的数值模拟值。

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

其中：

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

实施例一.经典的正规激波反射算例。下边界为反射边界，左边界和上边界为狄利克雷边界，右边界为超音速出口。计算区域为[0,4]×[0,1].其数值结果如图1所示。

实施例二.超音速圆柱绕流算例。计算区域为[-0.5,0]×[-0.5,0.5],圆柱半径0.05，圆心在原点，入流速度为3马赫，其数值结果如图2所示。

实施例三.超音速NACA0012翼型绕流算例。计算区域为[-1.5,1.5]×[-1.5,1.5]，入流速度为2马赫，冲击角度为1度。其数值结果如图3所示。

实施例四.超音速NACA0012翼型绕流算例。计算区域为[-1.0,2.0]×[-1.5,1.5]，入流速度为3马赫，冲击角度为5度。其数值结果如图4所示。

实施例五.亚音速NACA0012翼型绕流算例。计算区域为[-1.0,2.0]×[-1.5,1.5]，入流速度为0.8马赫，冲击角度为1度。其数值结果如图5所示。

实施例六.亚音速NACA0012翼型绕流算例。计算区域为[-1.0,2.0]×[-1.5,1.5]，入流速度为0.85马赫，冲击角度为1.25度。其数值结果如图6所示。

实施例七.亚音速圆柱绕流算例。计算区域为[-5,5]×[-5,5]。圆柱半径为0.5，圆心在原点。入流速度为0.38马赫。其数值结果如图7所示。

综上所述，本发明定点扫描算法推广到了计算复杂几何区域。ILW边界方法是高阶精度的边界处理方法，利用建立不等距模板的新型的WENO外推使得边界能够达到高阶精度且能有效避免振荡。快速扫描算法扫过物面边界时不需要再进行特殊处理，在网格点上有新点值就用新点值计算，没有新点值就用旧点值计算。通过在每个方向按坐标的拓扑结构顺序，上升和下降的顺序各扫描计算一次，可以极大加快迭代方法的迭代收敛速度。

本发明提出了一种新型ILW边界处理方法，并提出了一种将它与有限差分多重分辨加权基本无振荡格式结合的快速扫描法，能针对多种复杂计算区域的可压定常流场问题进行高精度数值模拟。本发明给出了该算法的具体的构造过程。相比于经典的ILW物面边界处理方法，通过设计不等距模板，利用多分辨率WENO的思想，使得线性权可以任取和为1的数。该方法不仅能保持物面边界附近计算格式高阶精度不降低的优点，同时能有效减少空间模板的使用数量和减少需要人为调整的参数，线性权也不需要求解可以任取，在求解定常问题时物理量的残差可以快速收敛，实施更加简单且易于实现。对于空间离散，多分辨率WENO离散相比于经典WENO空间离散通过设计不等距模板使得间断处精度可以降的更低，进而使得残差能够下降到机器零。快速扫描算法相比于三阶TVD Runge-Kutta时间离散格式，可以明显提高迭代效率，节约一半左右的迭代CPU时间。该快速扫描法是通过在所有方向上依次扫描，覆盖所有的迎风方向，再结合Gauss-Seidel迭代的特点，有新值用新值没有新值用旧值，可以提高欧拉向前时间离散格式的迭代效率，在保证格式收敛的情况下CFL数可以取较大值。Runge-Kutta时间离散格式由于多计算了两个虚拟时间层，相对计算效率比较低下。此格式的构造简单，数值精度较高，易于推广到多维情形。

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