WENO差分方法的一种边界处理技术

摘要

本发明提供了一种利用WENO差分数值方法进行流场计算时的一种边界处理技术，能够使得进行流场计算时达到一致高精度，并且能够处理复杂边界，结构网格上该发明的主要技术点和步骤为：(1)网格划分，将网格划分为内部区域和边界区域两个区域，边界区域只有一层网格，从而尽量减少边界区域DG的复杂计算。(2)空间离散，对内部区域利用WENO有限差分方法进行空间离散。(3)边界处理，对边界区域利用DG方法进行空间离散。(4)耦合处理，将边界处理部分与内部区域进行耦合，边界区域通过DG多项式逼近提供导数信息。(5)时间离散，通过TVD-Runge-Kutta方法进行显示时间离散。(6)后处理，通过Tecplot或者Paraview等商业软件对流场结果进行可视化模拟。

说明书

技术领域

本发明涉及计算流体力学数值方法领域，特别是涉及一种求解双曲守恒律方程的高精度有限差分方法的边界处理技术。 

背景技术

飞行器三维复杂流动的数值模拟和相关多目标优化问题是目前计算流体力学中的前沿热点问题，同时也是一个面向工程实际需求的应用问题。然而在目前的计算机规模和求解能力的条件下，目前主流的数值方法在计算效率上还不能满足这一工程实际应用问题的需要，解决这个问题的关键之一是提高流场解算器的效率。 

目前已经出现了多种为了充分利用流行的高精度数值方法的优点的耦合方法，例如WENO限制器间断有限元方法，高阶的SV方法，耦合的DG/FV方法，多区域耦合的DG和WENO方法，其中多区域耦合的DG和WENO方法具有高效、高精度且易处理复杂边界的优点，但是目前的所有的高阶方法中效率最高的仍然为有限差分方法，但是该方法一般都需要在结构网格上进行计算，且不能够处理复杂边界，处理边界时也没有紧致性，因而缺少一种能够帮助有限差分方法紧致处理复杂边界的高阶数值方法。 

发明内容

为了克服现有技术的不足，我们提出了一种能够处理复杂边界的WNEO有限差分方法作为处理WENO差分方法的一种边界处理技术，该方法发展了多区域耦合的DG和WENO-FD方法，更加充分的利用了有限差分方法的优势和DG的优势，几乎整个区域都利用WENO-FD方法计算，从而效率更高，仅仅在边界的一层网格处使用DG方法，边界处理利用DG方法处理，从而能够紧致处理复杂边界。 

本发明的技术方案是： 

将所求解的问题进行结构网格剖分，然后将所求解区域分为两种单元区域，一种是内部单元区域，另一种是边界单元区域，而且边界单元区域只有一层计算网格。 

利用传统的WENO-FD方法对内部单元区域进行空间离散，利用传统的DG方法对边界单元区域进行空间离散，并且处理边界。 

利用多区域耦合DG和WENO方法中的耦合方式实现WENO-FD方法和DG方法的耦合处理。耦合过程中需要构造内部单元区域边界处的数值通量，WENO数值通量的构造方法使用特殊的HWENO构造过程实现。 

边界单元处DG方法计算得到该单元上解的n阶多项式逼近结果，通过这个多项式逼近结果可以求出该单元中心点处的值以及不超过n阶的导数值。 

在空间离散结束之后，利用高阶的TVD Runge-Kutta方法求解得到下一时刻的数值结果。 

利用Tecplot等可视化软件对二维标量双曲守恒律方程和二维Euler方程组等验证算例的结果进行可视化处理。 

本发明的有益效果是： 

本发明充分利用了高阶WENO有限差分方法高效简单的优点，结合了DG方法能够处理复杂边界并且具有局部性和紧致性的特点，发展了多区域耦合DG和WENO方法，提出了一种基于DG处理边界的WNEO-FD方法，该方法仅仅在边界一层网格处使用DG方法，在处理大规模问题是效率几乎等价于WENO-FD方法，而且也具有处理复杂边界的能力，并且在处理复杂边界时不再需要构造多个虚拟网格点，减小了WENO-FD边界处理时的模板大小从而处理时具有紧致性和局部性。 

附图说明

图1为一维算例网格划分示意图 

图2为内部单元处传统WENO模板示意图 

图3为内部单元边界处HWENO模板示意图 

图4为一维标量Burgers方程算例示意图 

图5为一维精度验证表格示意图 

图6为二维双马赫反射算例结构网格划分示意图 

图7为双马赫反射算例描述图 

图8为二维双马赫反射算例密度等值线示意图 

图9为算法实现过程流程图 

具体实施方式

首先我们针对一维标量Burgers方程说明我们的算法主要过程，如图1第一步的区域划分为两个部分，其中左边WENO-FD部分是内部单元区域，右边一层DG网格为一层边界单元区域。 

内部区域WENO-FD的具体实施过程如下: 

首先我们不考虑通量分裂，假设f'(u)＞0,那么在单元I_i守恒的有限差分格式形式为 

$\frac{d\left(u_i\left(t\right)\right)}{\mathrm{dt}}+\frac{1}{\mathrm{\Delta x}}({\hat{f}}_{i+\frac{1}{2}}-{\hat{f}}_{i-\frac{1}{2}})$

其中u_i(t)＝u(x_i,t)是单元I_i中点处的点值，为数值通量可以表示为 

$\hat{f}=f(u_{i-r},...,u_{i+s})$

如图2除了内部单元边界两个通量外，大部分内部单元通量我们选用图2的模板S＝(S₁,S₂,S₃)，通量 $\hat{f}=f\left(u_{i+\frac{1}{2}}\right),$其中$u_{i+\frac{1}{2}}=w_1{u}_{i+\frac{1}{2}}^1+w_2{u}_{i+\frac{1}{2}}^2+w_3{u}_{i+\frac{1}{2}}^3,$然后$u_{i+\frac{1}{2}}^s,s=\mathrm{1,2,3}$是由每一个小模板上包含单元的值通过拉格朗日插值得到的，对于5rd的WENO-FD小模板是通过3rd拉格朗日插值得出，最终的插值结果为: 

$\left(\begin{array}{c}{u}_{i+\frac{1}{2}}^1=-\frac{1}{6}\left(u_{i-1}\right)+\frac{5}{6}\left(u_i\right)+\frac{1}{3}\left(u_{i+1}\right)\\ u_{i+\frac{1}{2}}^2=\frac{1}{3}\left(u_i\right)+\frac{5}{6}\left(u_{i+1}\right)-\frac{1}{6}\left(u_{i+2}\right)\\ u_{i+\frac{1}{2}}^3=\frac{11}{6}\left(u_{i+1}\right)-\frac{7}{6}\left(u_{i+2}\right)+\frac{1}{3}\left(u_{i+3}\right)\end{array}\right)$

非线性权w₁,w₂,w₃的计算过程中满足w_r＝d_r+O(Δx³)： 

$w_r=\frac{{\alpha }_r}{\underset{s=1}{\overset{3}{\Sigma }}{\alpha }_s},{\alpha }_r=\frac{d_r}{{(ϵ+{\beta }_r)}^2},r=\mathrm{1,2},3$

其中β_r被称作“光滑因子”，这里选用光滑因子： 

${\beta }_r=\underset{l=1}{\overset{2}{\Sigma }}{\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\Delta x^{2l-1}{\left(\frac{{\partial }^l{p}_r\left(x\right)}{{\partial }_{l}x}\right)}^2\mathrm{dx}.$

非线性权的计算过程需要计算线性权：$u_{i+\frac{1}{2}}=d_1{u}_{i+\frac{1}{2}}^1+d_2{u}_{i+\frac{1}{2}}^2+d_3{u}_{i+\frac{1}{2}}^3=u\left(x_{i+\frac{1}{2}}\right)+O\left({\mathrm{\Delta x}}^5\right),$这样的线性系数$d_1=\frac{3}{10},d_2=\frac{3}{5},d_3=\frac{1}{10}.$

如图3，如果我们要计算内部单元边界的两个单元边界通量，我们的方法中使用更小的模板，图中的虚拟网格点我们不再需要，而是在这两个相对小的模板上采用特殊的Hermite插值得到我们要的每一个小模板的插值结果Hermite插值的结果为： 

$\left(\begin{array}{c}{u}_{N-\frac{3}{2}}^1=-\frac{1}{6}\left(u_{N-3}\right)+\frac{5}{6}\left(u_{N-2}\right)+\frac{1}{3}\left(u_{N-1}\right)\\ u_{N-\frac{3}{2}}^2=\frac{1}{3}\left(u_{N-2}\right)+\frac{5}{6}\left(u_{N-1}\right)-\frac{1}{6}\left(u_N\right)\\ u_{N-\frac{3}{2}}^3=\frac{13}{6}\left(u_{N-1}\right)-\frac{7}{6}\left(u_N\right)+\frac{2}{3}\mathrm{\Delta x}\left(u_N^{\prime }\right)\end{array}\right)$

和 

$\left(\begin{array}{c}{u}_{N-\frac{1}{2}}^1=-\frac{1}{6}\left(u_{N-2}\right)+\frac{5}{6}\left(u_{N-1}\right)+\frac{1}{3}\left(u_N\right)\\ u_{N-\frac{1}{2}}^2=\frac{1}{6}\left(u_{N-1}\right)+\frac{5}{6}\left(u_N\right)-\frac{1}{3}\left(u_N\right)\\ u_{N-\frac{1}{2}}^3=\left(u_N\right)-\frac{1}{2}\mathrm{\Delta x}\left(u_N^{\prime }\right)+\frac{1}{12}\Delta x^2\left(u_N^{\prime \prime }\right)\end{array}\right)$

然后再通过同样的方式得到我们需要的非线性权和光滑因子，从而可以求解方程得到完整的解。 

根据技术方案中的步骤我们还需要通过DG方法得到如图3边界单元I₁andI_N处单元中点处的值以及高阶导数值，例如I_N单元得到的多项式逼近解为: 

u(x)＝u₀p₀(x)+u₁p₁(x)+u₂p₂(x) 

系数u₀,u₁,u₂是由DG方法逼近得到的自由度。高阶导数信息也很容易通过多项式求导得出 

u′(x)＝u₀p′₀(x)+u₁p′₁(x)+u₂p′₂(x) 

u″(x)＝u₀p″₀(x)+u₁p″₁(x)+u₂p″₂(x) 

根据技术方案中的过程我们需要将两种方法的边界耦合到一起，DG和HWENO耦合过程如下： 

$\left(\begin{array}{c}\frac{d}{\mathrm{dt}}{u}_1-\frac{1}{\Delta x_1}({\tilde{f}}_{\frac{1}{2}}-{\tilde{f}}_{i+\frac{1}{2}})=0,\\ \frac{d}{\mathrm{dt}}{u}_i-\frac{1}{\Delta x_i}({\hat{f}}_{i-\frac{1}{2}}-{\hat{f}}_{i+\frac{1}{2}})=0,i=2,...,N-1,\\ \frac{d}{\mathrm{dt}}{u}_N-\frac{1}{\Delta x_N}({\tilde{f}}_{N-\frac{1}{2}}-{\tilde{f}}_{N+\frac{1}{2}})=0,\end{array}\right)$

其中代表通过通常的WENO重构得到的通量，代表DG方法得到的通量, 是通过上述HWENO重构的方式得到。 

为了验证上述算法，如图4我们针对一维标量的Burgers方程计算了强激波算例，并且我们针对线性的一维算例进行了算法的精度验证如图3。 

为了验证上述算法，我们针对二维双马赫反射算例进行了验证，双马赫反射问题的网格剖分如图6，其中上面的网格部分代表内部单元网格，最下面一层网格代表边界单元网格。 

如图7，假设初始时刻，马赫数Ma＝10的正激波向平板运动。激波与平板之间成60°夹角。通过上述描述的间断有限元处理边界的WENO有限差分方法计算流场的演化。 

如图3和图7，我们根据我们所述的间断有限元处理边界的WENO有限差分方法的计算结果，通过Tecplot可视化软件对我们的流场计算结果进行了可视化演示，图7中演示了双马赫反射问题中流场密度的变化，计算网格为320*80，初始马赫数为Ma＝10，图中演示了从ρ＝1.3965到ρ＝22.682总共30条密度等值线。 

本发明中没有详细说明的内容均属于本技术领域人员熟知的常识。