一种高效的隐式间断Galerkin方法研究

郭永恒， 杨 永， 张 强

（西北工业大学 翼型叶栅空气动力学国防科技重点实验室，陕西 西安710072）

0 引 言

在计算流体力学领域，为了高精度地求解非线性偏微分方程组，研究更加复杂的流动现象，间断Galerkin方法已经引起人们越来越多的关注。特别是经过Cuckburn和Shu的长期探索，一种具有TVD性质的显式Runge－Kutta间断Galerkin（RKDG）格式得以逐步完善［1］，被广泛应用于双曲守恒律问题的数值求解，取得了大量令人满意的结果，显示了间断Galerkin方法的优越性。然而，美中不足的是，随着逼近精度的提高，RKDG格式对应的稳定性条件将越来越严格，时间步长受到明显的限制，从而导致更多CPU时间的消耗。对于定常流场的计算问题，尽管可以把当地时间步长技术与RKDG格式相结合，在一定程度上加速收敛过程，但是即便如此，最大的时间步长依然受到当地稳定性条件的限制［2］。缓慢的收敛速度在很大程度上制约着RKDG方法在工程中的应用。基于以上分析，本文建立了一种隐式间断Galerkin（Implicit Discontinuous Galerkin，IMDG）求解器，并通过对翼型亚声速和跨声速流场的模拟，检验了该求解器的计算效率。

1 间断Galerkin方程

在二维区域D上，Euler方程可以写成如下守恒形式

Q和F＝（F，G）T分别代表守恒型流动变量和通量向量。把D划分成Ne个互不重叠的子区域Di，并且设解函数空间为

P（Di）是定义在Di上的多项式空间。设Di上的测试函数集合为，它符合如下条件

为了控制数值离散产生的小量误差，使用加权残量法［3］，在每个单元上作如下内积运算

运用Green公式，我们由方程（4）得到间断Galerkin方程

∂Dij和n分别表示当前单元的边界和相关的单位外法向量。在式（5）左边第三项的积分中，我们使用迎风格式对数值通量进行计算，实现相邻单元之间的信息传递［4］：

在基函数的构造过程中，我们应用Gram－Schmidt方法对多项式序列

进行规范正交化，所得的结果即为数值格式中使用的基函数，其具体形式为：

三角形单元经过坐标变换后，在计算区域D＝｛（ξ，η）｜0≤η≤1－ξ，0≤ξ≤1｝上，基函数满足如下条件：

这样，式（5）第一项对应的质量矩阵就成为对角矩阵，使离散格式得到明显简化。当采用前三项时，离散格式为二阶精度，当采用前六项时，离散格式为三阶精度。在本文的数值实验中，显式格式和隐式格式一律设置为二阶精度用来比较收敛速度。

2 隐式时间离散格式

我们在时间方向上运用Euler向后差分，并把第n和第n＋1时间层之间物理量X的增量记为δXn，即

那么，在第i号单元上，间断Galerkin方程（5）就变形为

其中，Ri，l（Qn）为残差向量，即

设整个流场中守恒型变量为

Q对应的有限元系数向量为

这样，Di上的有限元解就可以表示为

对式（10）左侧第二、第三项的被积表达式做线性化处理，得

A为系数矩阵，R为残差向量

R的分量记为

对比式（10）和式（18）可以发现，对于定常问题，IMDG格式退化成Newton迭代法，它具有二次收敛阶［5］。因为在每一个时间步上，建立系数矩阵A需要进行大量的数值积分，所以受文献［5］的启发，我们对式（18）进行一个重要修正，得

其中，

“mod（n－1，k）＝0”表示n－1被k（k≥2）整除，这样，每进行k次迭代，残值向量R更新k次，而系数矩阵A只更新一次，与式（18）相比，计算量大大降低，迭代过程也从整体上得到了进一步的优化。

3 线性系统的求解

通过类比，我们把有限体积法中的LU－SGS方法［6］进行推广，用来求解大型稀疏线性系统（22）。首先，把矩阵A进行如下分解

L、D和U分别为严格的分块下三角矩阵、分块对角矩阵和分块上三角矩阵。令矩阵

我们可以根据如下两个步骤对线性系统（25）进行快速求解

4 计算结果分析

为了测试IMDG方法的计算效率，我们分别求解了NACA0012翼型对应的亚声速和跨声速流场，并与RKDG格式得到的结果作比较。首先，在NACA0012翼型周围生成非结构网格，如图1、图2所示。

对于亚声速情形，设定计算状态为：Ma∞＝0．63，α＝2°。在IMDG格式中，逐步加大CFL数，可以发现该格式是无条件稳定的，这说明本文的隐式格式具有优良的稳定性。图3是压强系数分布曲线，可以看出，RKDG和IMDG对应的结果几乎完全一致，这说明二者的计算精度是相同的。图4是残值随迭代步数的变化曲线，图5是残值随CPU时间的变化曲线，可以看出IMDG对应的残值不仅在大范围内单调下降，而且当下降到相同量级时，IMDG使用的迭代步数和CPU时间分别比RKDG节省了90%和85%以上，收敛速度几乎提高了一个数量级。

为了进一步检测IMDG求解器的计算效率，我们求解了NACA0012翼型的跨声速流场，设定计算状态为：Ma∞＝0．8，α＝1．25°。在此算例中，本文没有附加任何限制器，依然能够得到收敛的结果，如图6～图8所示。

可以看出，IMDG能够捕捉到位置和RKDG完全一致的激波；同时，在高效计算方面，它再一次展示了自身的优越性。

5 结 论

为了提高间断Galerkin方法求解定常流场问题的效率，本文建立了与之相关的隐式离散格式，并在一定程度上对迭代过程进行了优化。数值实验表明该格式是无条件稳定的，这非常有利于计算效率的大幅提高。今后，可以在本文基础上展开更深层次的探索，特别是简化隐式格式的建立过程，提高相关的大型稀疏线性系统的求解精度，从而使IMDG的计算效率进一步增强，为间断Galerkin方法在工程计算中的广泛应用打下坚实的基础。

［1］ COCKBURN B，SHU C W．TVD Runge－Kutta local projection discontinuous Galerkin finite element method for scalar conservation laws II：General Framework［J］．Math．Comp，1989，（52）：411－435．

［2］ PATRICK RASETARINERA，HUSSAINI M Y．An efficient implicit discontinuous spectral Galerkin method［J］．Journal of Computational Physics，2001（172）：718－738．

［3］ 王烈衡，许学军．有限元方法的数学基础［M］．北京：科学出版社，2004．（WANG L H，XU X J．The mathematical foundations of the finite element method［M］．Beijing：Science Press，2004．）

［4］ ROE P L．Approximate Riemann solver，parameter vectors and different schemes［J］．Journal of Computational Physics，1981，43：357－372．

［5］ 奥特加J M，莱因博尔特 W C．多元非线性方程组迭代解法［M］．北京：科学出版社，1983．（ORTEGA J M，RHEINBOLDT W C．Iterative solution of nonlinear equations in several variables［M］．Beijing：Science Press，1983．）

［6］ JAMESON A，TURKEL E．Implicit scheme and LU－decompositions［J］．Math．Comput．，1981，37：385－397．