基于区域分解的耦合高精度复杂外形流场快速算法

摘要

本发明提供了一种基于区域分解的耦合高精度DG和WENO方法求解双曲守恒律方程和Euler方程组的快速计算方法，该方法首先对原问题进行区域分解，即在物理边界附近区域采用结构或非结构的间断有限元方法(DG)，在其余规则区域使用结构网格下的有限体差分型加权本质无振荡(WENO)格式。在处理区域耦合界面过程中，有两种处理方法，一种是守恒的耦合处理方法；另一种是非守恒的耦合处理方法。在实际计算过程中，我们使用坏单元指示器来判断界面附近解是否充分光滑，若解充分光滑，则界面处采用非守恒的耦合方法，否则，采用守恒的耦合方法。

说明书

技术领域：

本发明涉及计算流体力学数值方法领域，特别是涉及一种求解双曲守恒律方程的高精度数值方法。

背景技术：

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

目前流行的高精度数值方法主要包括间断有限元方法(DG)，高精度有限体积方法，如(k-exact)有 限体积方法以及高精度有限差分型方法，如有限差分型的加权本质无振荡格式(WENO)。间断有限元方 法具有高精度和易处理复杂边界这些优点，但其计算量大，计算效率低，不能满足工程实际需求；高精度 有限体积方法有传统的有限体积法推广而来，具有处理复杂边界的能力，但其重构模板一般不是紧致的， 这给该方法应用到实际三维问题上带来了一定的困难，另外此方法的计算量也比较大，计算效率较低；有 限差分型方法具有高精度和计算量小，计算效率高的优点，但有限差分方法一般只能在结构网格上应用， 很难处理复杂物形和边界。故针对三维飞行器复杂流动数值模拟和相关多目标优化这一问题，目前缺少一 种高精度、高效且易处理复杂边界的数值方法。

发明内容：

为了克服上述现有技术的不足，本发明提出一种基于区域分解的耦合DG和WENO方法，即将间断 有限元(DG)方法和有限差分型的加权本质无振荡(WENO)格式以区域分解的方式耦合在一起，在复 杂物形边界附近区域使用结构或非结构网格下的间断有限元方法处理计算区域的边界，同时在流场其余规 则区域使用有限差分型WENO方法以大规模的提高计算效率。多区域耦合DG和WENO方法具有高阶精 度，易处理复杂边界且计算量小，计算效率高等优势。

本发明的技术方案是：

对求解的实际问题首先进行区域分解，将整体计算区域划分为物形边界附近区域和其余规则计算区 域；对物形边界附近区域使用结构或非结构网格进行区域剖分，对其余规则计算区域采用结构网格进行区 域剖分；

对物形边界附近区域使用结构或非结构网格下的间断有限元方法进行初始化，对其余规则计算区域采 用结构网格下的有限差分型WENO进行初始化；

在耦合界面处分别构造两类数值通量，一类是耦合界面处DG数值通量，一类是耦合界面处有限差分 型WENO数值通量；

使用坏单元指示子判断耦合界面两侧邻居单元是否为坏单元，若存在坏单元，说明解在界面附近可能 存在间断，故在界面处使用是守恒的耦合方式，若不存在坏单元，说明解在界面附近充分光滑，故在界面 处使用非守恒的耦合方式；

确定界面处耦合方式之后，可对计算区域内各个子区域完成相应的空间离散，得到半离散个是的方程 (组)，此半离散方程(组)可使用三阶TVD Runge-Kutta方法求解。

本发明的有益效果是：

本发明结合了目前主流高精度算法的优点，采用区域分解的方式，耦合两种不同种类的高精度算法， 从而达到易于适应和处理各种复杂边界的同时，大幅度提高计算效率的目的。耦合的DG和WENO方法 可以方便的在混合网格上使用，在处理复杂物形边界附近区域时使用非结构网格的DG方法，在处理远场 规则计算区域时，使用结构网格的有限差分型WENO方法。相比于传统的DG方法，耦合方法可以幅度 降低计算量，提高计算效率；相比于传统的有限差分型格式，耦合方法在处理复杂边界时更加灵活，从而 满足工程实际的需要。耦合方法作为计算流体力学中求解双曲守恒律的一类数值方法，在工程实际应用中 有广阔的前景和应用价值。

附图说明：

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面对实施例的附图做简要介绍。

图1：NACA0012翼型示意图

图2：耦合算法针对计算NACA0012机翼绕流问题的区域分解示意图

图3：耦合算法混合网格下耦合界面局部示意图

图4：构造耦合界面处WENO-FD数值通量示意图

图5：构造耦合界面处DG数值通量示意图

图6：耦合算法计算亚音速NACA0012机翼绕流计算结果密度等值线示意图

图7：耦合算法整体流程示意图

具体实施方式：

下面将结合本发明实施例中的附图，对本发明实施例中的方法进行清楚、完整的描述。显然，所描述 的实例仅仅是本发明的一个应用实例。基于本发明中的实例，本领域技术人员在没有做出创造性工作前提 下所获得的所有其他实例，都属于本发明保护的范围。

针对图一物形，确定计算区域。在此问题中，需要我们进行数值模拟求解的问题时亚音速NACA0012 机翼绕流问题，流动初始条件为马赫数Ma＝0.4，攻角AoA＝5.0°。我们选定的计算区域为一规则的矩形区 域[-15.0,15.0]×[-15.0,15.0]；

针对求解问题物形特征，进行区域分解。对于本问题，由于处理物形的非规则区域集中于机翼附近区 域，因此我们将整体计算区域分为两大区域，区域一为机翼附近的非结构网格区域(图二中红色区域)， 该区域范围为[-0.4,1.4]×[-0.4,0.4]，该区域我们使用非结构网格上的DG方法计算；区域二为规则的远场 结构网格区域(图二中绿色区域)，该区域为规则区域，可以使用结构网格的有限差分型WENO-FD方法 计算。

确定网格尺度，进行网格剖分。针对本问题的特点和NACA0012机翼的物形特征长度，我们的网格剖 分参数如下：网格尺寸h＝0.05，非结构区域三角形网格单元数量N₁＝824，结构区域四边形网格单元数量 N₂＝57456，计算区域总体网格数量为N＝58520。

初始化流场区域，构造耦合界面处数值通量。在耦合界面处我们需要构造两类数值通量，分别是耦合 界面处WENO-FD的数值通量和DG的数值通量。

构造耦合界面处WENO-FD数值通量步骤如下：

寻找和获得耦合界面构造WENO-FD数值通量所需的虚拟节点的位置及节点上的函数值U_h，如图4 所示，设耦合界面的位置为在构造单元I_I+1，J处的WENO-FD数值通量时，需要虚拟节点 I_I，J，I_I-1，J，I_I-2，J，这三个节点上的函数值U_h由该虚拟节点对应的三角单元上的DG解函数多项式提供：

$U_h(x,y)=\underset{l=0}{\overset{K}{\Sigma }}{u}^{\left(l\right)}\left(t\right)v_{\left(l\right)}(x,y)$

其中u^(l)(t)为单元上DG自由度，v_(l)(x，y)相应的基函数；

在由虚拟节点I_I，J，I_I-1，J，I_I-2，J和WENO-FD区域内节点I_I+1，J，I_I+2，J做成的重构模板上，使用 WENO-FD数值通量构造方法构造单元I_I+1，J位于耦合界面处的WENO-FD数值通量

构造耦合界面处DG数值通量的步骤如下：

构造耦合界面处的DG数值通量需要提供界面左右两侧高斯积分节点上的和如图5所示，其中可以通过DG求解区域内相应单元上的解函数多项式获得，对于则需要使 用重构或插值手段获得，这里我们使用基于WENO的直接插值方法获得位于耦合界面处的对于 二维情况下，基于WENO的直接插值，我们逐维度进行插值的方法，即先沿y-轴方向进行WENO插值， 得到x-轴方向WENO插值所需节点值，再进行x-轴方向的WENO插值，从而得到耦合界面处位于高斯积 分节点上的

以耦合界面上点为例，给出沿x-轴方向进行WENO型插值的方法如下：

选取x-轴方向WENO插值的插值模板S＝{I_I-2，J，I_I-1，J，I_I，J，I_I+1，J，I_I+2，J}，将此模板分割为三 个小模板S₁＝{I_I-2，J，I_I-1，J，I_I，J}，S₂＝{I_I-1，JI_I，J，I_I+1，J}，S₃＝{I_I，J，I_I+1，J，I_I+2，J}；

在每个小模板上构造拉格朗日插值多项式P_l(x)，l＝1，2，3，对于每个模板S_l内的拉格朗日多项式，需 满足P_l(x_i，y_j)＝U_i，j，I_i，j∈S_l；

计算每个小模板上的拉格朗日多项式的线性权重d_l，l＝1，2，3和光滑因子β_l，l＝1，2，3，得到每个小 模板对应的非线性权，对于本例三个小模板S_l的线性权分别为：

$d_1=\frac{5}{16},d_2=\frac{10}{16},d_3=\frac{1}{16}$

每个小模板上的光滑因子β_l的计算方法如下：

${\beta }_l=\underset{k=1}{\overset{N}{\Sigma }}{\int }_{x_{I+\frac{1}{2}}}^{x_{I+\frac{3}{2}}}{\mathrm{\Delta x}}^{2k-1}{\left(\frac{{\partial }^k{P}_l\left(x\right)}{{\partial }^{k}x}\right)}^2\mathrm{dx},l=\mathrm{0,1,2}.$

其中N为小模板上拉格朗日插值多项式的次数。于是，可以得到每个小模板对应的非线性权如下：

$w_l=\frac{{\alpha }_l}{{\Sigma }_{s=0}^2{\alpha }_l},{\alpha }_l=\frac{d_l}{{(ϵ+{\beta }_l)}^2},l=\mathrm{0,1,2},$

其中ε＝1.0e^-6；

计算耦合界面处点的函数值，该值由上述模板中各个多项式的值与相应的非线 性权重加权组合得到

$U_{I+\frac{1}{2}.J+\frac{1}{2}}^{+}=\underset{l=0}{\overset{2}{\Sigma }}{w}_l{P}_l\left(x_{1+\frac{1}{2},J+\frac{1}{2}}\right)$

在完成耦合界面处WENO-FD数值通量和DG数值通量的构造后，在耦合界面位置使用坏单元指示 子，判断耦合界面两侧单元是否为坏单元：若耦合界面两侧单元存在坏单元，则耦合界面处使用守恒的耦 合方式，否则耦合界面处使用非守恒的耦合方式；对于守恒的耦合方式，即在耦合界面处采用唯一的数值 通量，这里我们可以选取WENO-FD的数值通量，或选取DG数值通量；对于非守恒耦合方式，耦合界面 处不同的子区域选择使用不同的数值通量，例如对于DG子区域，耦合界面处则使用DG数值通量，对于 WENO-FD子区域，耦合界面处使用WENO-FD数值通量。

最后，对于不同的子区域，使用相应的重构的数值通量完成空间离散，得到半离散形式的常微分方 程(组)，可以使用三阶TVD Runge-Kutta方法求解，本实例的计算结果如图6所示。