一种基于流场通量阶跃的高精度间断Galerkin人工粘性激波捕捉方法与流程

本文涉及计算流体力学技术流域，具体涉及一种基于流场通量阶跃的高精度间断Galerkin人工粘性激波捕捉方法。

背景技术：

高精度DG方法目前在数值耗散和色散方面的优良特性，理论上非常适合解决流体计算中的复杂多尺度问题，高精度DG方法在低速不可压缩流体的计算中取得了一定的成果并得到了广泛应用。但是当高精度DG方法应用到可压缩流体的计算时，流场中会产生激波，根据Godunov原理，高精度DG方法在激波间断附近会产生Gibbs现象，引起非物理解的产生，导致计算中断。激波捕捉已经成为阻碍高精度DG方法在可压缩流体计算发展的主要瓶颈。目前，高精度DG激波捕捉方法都不够成熟，以往的激波捕捉方法主要包括限制器和重构。这些方法都是采用后处理的方式对激波进行捕捉，当采用这些方式时，整个求解过程中会引入不可控的非线性项，由于非线性项的引入，物理过程与控制方程之间会存在不匹配的情况，这就会直接导致激波区域出现严重降解、收敛性差、鲁棒性差、容易出现非物理震荡，激波捕捉不够“锋利”等问题。

技术实现要素：

本发明的目的是提出一种高精度间断Galerkin人工粘性激波捕捉方法，该方法根据激波的特性和量纲分析原理，采用流场中单元界面上的通量阶跃，构造了一种新的人工粘性激波捕捉方法，主要用于解决高精度DG方法在激波捕捉过程中遇到的收敛性差、鲁棒性差、容易出现非物理震荡，激波捕捉不够“锋利”等问题。

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

一种基于流场通量阶跃的高精度间断Galerkin人工粘性激波捕捉方法，包括以下步骤：

步骤一：建立高精度DG框架，包括网格剖分、Euler控制方程、有限元方法中的基函数、测试函数、Gauss积分点等信息；

步骤二：以单元交界面上的通量阶跃为基础构造单元内的的人工粘性系数；

步骤三：将人工粘性系数带入到Euler控制方程，求解得到仿真计算结果。

本发明采用非结构网格对计算区域进行剖分，控制方程采用Euler方程，建立以基函数、测试函数、Gauss积分点为代表的DG高精度框架。同时在方程中以单元交界面上守恒变量的阶跃为基础构造新的人工粘性项，方程的对流项采用HLL格式离散求解。在有效捕捉激波的情况下，保证鲁棒性和计算精度。

与现有的手段比较，本发明的特点在于：

1、即使在不进行激波捕捉的情况下，单元交界面处的通量也是方程求解必须的中间变量，选择流场单元交界面处的通量阶跃来构造人工粘性，相比其他方法可以减少计算量，从而节省计算的时间。

2、采用通量阶跃人工粘性方法进行激波捕捉，压力分布曲线在激波附近更“锋利”，同时收敛性优于其他方法。

附图说明

本发明将通过例子并参照附图的方式说明，其中：

图1为本发明方法的流程图；

图2为采用本发明方法计算的一维sod问题；

图3为采用本发明方法对NACA0012翼型计算得到的收敛曲线与其他方法的对比。

具体实施方式

本说明书中公开的所有特征，或公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合。

本发明公开的一种基于流场通量阶跃的高精度间断Galerkin人工粘性激波捕捉方法，包括三个部分，如图1所示。

第一部分：建立高精度DG框架，包括网格剖分、Euler控制方程、有限元方法中的基函数、测试函数、Gauss积分点等信息的计算和存储。包括以下步骤：

步骤101、采用非结构网格对计算区域进行网格剖分，对于二维计算域，剖分的网格类型包括三角形和四边形，对于三维计算域，网格类型包含四面体、六面体、三棱柱和金字塔外形。

步骤102、构建微分形式下的Euler方程

其中U代表流场中的守恒量，U是一个矢量，U＝(ρ,ρu,ρv,ρw,ρE)^T，代表守恒量相对于时间t的偏导数，F^c代表守恒通量，代表守恒通量的散度。

步骤103、选择Taylor基作为基函数和测试函数，流场中的守恒量采用基函数的线性组合表示计算不同类型网格下的体积分Gauss积分点和面积分Gauss积分点，并在内存中进行存储备用。

步骤104、将守恒量的线性组合带入Euler控制方程(1)，并对方程进行积分，然后在方程两边同时乘以基函数φ，利用格林高斯公式，得到弱形势下的DG求解方程

其中，M＝∫Ωφiφj代表质量矩阵，u代表守恒变量用基函数线性组合时采用的系数又称为自由度，代表自由度对时间的导数。为单元边界面，为单元边界面的外法线方向，代表基函数的梯度。

第二部分：以单元交界面上的通量阶跃为基础构造单元内的的人工粘性系数，包括以下步骤：

步骤201、选择拉普拉斯人工粘性模型，带入公式(1)并重复步骤103和104得到包含人工粘性项的DG求解方程

其中，ε为人工粘性系数，决定单元内人工粘性的数值，代表自由度的梯度。

步骤202、对人工粘性系数ε进行重新构造，选择单元交界面处通量的阶跃(和代表单元交界面左右两侧的通量值)以及守恒变量的平均值(和代表单元交界面左右两侧的守恒变量值)进行线性组合，构造交界面处的中间阶跃Ujump

步骤203、将步骤202中得到的中间阶跃量Ujump在单元交界面处进行积分，然后积分量除以单元的总面积S，得到单元体内阶跃量的分布Ucell

步骤204、利用单元体内的阶跃分布Ucell和单元体心处的压力梯度经验参数C∈以及本地网格单元的参考尺度h，三者相乘，然后除以单元体心处的压力p，从而构造出单元内的人工粘性系数ε

第三部分：将人工粘性系数带入到控制方程，求解得到仿真计算结果；

步骤301、将步骤204中得到的单元内的人工粘性系数ε带入到包含人工粘性项的DG求解方程方程(3)中。

步骤302、对DG求解方程进行离散，通过迭代计算求解控制方程，得到仿真的气动结果和流场。

如图2所示，为采用本发明方法计算的一维sod问题，可以看到在激波附近曲线分布更“锋利”，对于激波区域，其压力或者密度曲线在理论上应该是一条垂直的线，如果计算结果越“锋利”，其结果就越接近于理论值，从图2中可以明显看出，新方法在“锋利”程度上比其他方法更优越，结果更接近理论值。对于图2中新方法的计算结果，图3所示为采用本发明方法对NACA0012翼型计算得到的收敛曲线与其他方法的对比，可以明显看到该方法在计算收敛时间上明显少于其他方法。

本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中披露的新特征或任何新的组合，以及披露的任一新的方法或过程的步骤或任何新的组合。