传统的间断有限元方法是一种空间离散方法，为了求解控制方程经空间离散后得到的常微分方程组，需要合适的时间离散方法。以Runge-Kutta方法为代表的显式方法具有便于编程、易于达到高精度等优点，得到了广泛的应用[2-4]。尽管如此，显式方法的时间步长受到了稳定性的严格限制，整体上效率并不高。相反，尽管编程复杂、单步计算量大，隐式方法对时间步长的限制却远比显式方法松，部分隐式方法甚至不限制时间步长的大小[5]，整体上，相较显式方法，隐式方法所需的计算时间更少。隐式方法与间断有限元方法的结合正逐渐成为相关领域的趋势[6-8]。
采用隐式方法，每个时间步都要求解非线性方程组，一种经典、有效的方法是牛顿法。牛顿法的每个迭代步都要求解线性方程组，基于非结构网格、采用间断有限元方法进行空间离散时，该线性方程组的系数矩阵具有阶数大、稀疏性强、数值非对称等特点，基于此，通常采用广义极小残量(Generalized Minimal Residual，GMRES)方法求解该线性方程组。该方法具有内存占用少、计算能力节约的优点，是当下求解大型稀疏线性方程组的主流方法[9-10]。
矩阵带宽缩减技术常被用于以直接方法求解线性方程组[11-13]。基于以隐式间断有限元方法求解定常、二维、无黏流动问题，本文将该技术应用于以GMRES方法求解线性方程组。这样做能够加快GMRES方法的收敛速度，增大程序可承受的时间步长；另外，该技术的引入不会破坏原有的程序结构，易于编程实现。
1 基本技术与方法 1.1 矩阵带宽缩减技术 采用直接方法求解线性方程组
 | (1) |
时，需要将A分解成下三角矩阵Lo与上三角矩阵Up的乘积，然后求解2个三角形方程组：
 | (2) |
若A为稀疏矩阵，且存在远离对角线的非零元素，则Lo、Up中的非零元素将远比A中的多。过多的非零元素不仅占用大量的存储空间，还导致式(2)的求解非常缓慢。所谓矩阵带宽缩减，就是通过不断交换A的行、列，将其非零元素集中到对角线附近。设P为一系列初等置换矩阵的乘积，则调节矩阵非零元素分布的过程可以由PAPT表示。记
 | (3) |
由于
因为只涉及到矩阵元素的位置交换，实际应用中，用于缩减矩阵带宽的算法多是基于图论构建的，而不是真的构建出置换矩阵P。给定一个稀疏的矩阵，如果将其看作一个反映结点连接关系的邻接矩阵，就可以根据它绘制出一张由结点和连接线构成的网络图，而矩阵的带宽缩减问题就转化成基于这张图的编号排序问题。下面举例说明。
图 1(a)定性地显示了一个稀疏矩阵，其中的圆点代表了矩阵的非零元素，空位代表 0，则相应的网络图如图 1(b)所示。仅依据后者就能看出，矩阵的最大带宽为4，分别出现在矩阵第1行的第1、5元素，及矩阵第6行的第2、6元素之间。前述基于图论的算法，就是将图 1(b)的编号重新排列，使得任意两个由线连接的编号尽可能接近。
 |
| 图 1 稀疏矩阵和结点网络(带宽未经缩减) Fig. 1 Sparse matrix and corresponding node network(without bandwidth reduction) |
| 图选项 |
基于图论的结点重排方法有很多，比较典型的有CM(Cuthill-Mckee)方法、RCM(Reverse Cuthill-Mckee)方法、GPS方法、Burgess方法等。下面简单介绍一下当前应用非常广泛的RCM方法的操作步骤[11-13]。
步骤1 ??准备两个空队列Y、Z。
步骤2??从网络图中选择一个编号尚未出现于Y的、分支最少的结点，记为B。
步骤3??将B的编号置于Y的第一个空位。
步骤4??按照分支数目递增的顺序，将与B相连的结点的编号依次置于Z中。
步骤5??记Z中的第一个编号代表的结点为C。若C的编号尚未出现于Y，则将其置于Y的第1个空位，同时按照分支数目递增的顺序，将尚未出现于Y的、与C相连的所有结点的编号置于Z中。
步骤6??若Z非空，则回到步骤5。
步骤7??若图中有未被涉及的结点(比如孤立的结点)，则回到步骤2。
步骤8 ??将Y中的排列顺序倒过来。
按照上述步骤，由图 1(b)，可以得到队列：
 | (4) |
将图 1(b)中的编号i置于编号Y(i)的位置，得到的新编号如图 2(b)所示。分别记编号调整前、后的稀疏矩阵为
 |
| 图 2 稀疏矩阵和结点网络(经RCM缩减带宽) Fig. 2 Sparse matrix and corresponding node network(with RCM bandwidth reduction) |
| 图选项 |
 | (5) |
图 2(a)描绘了
1.2 隐式间断有限元方法
1.2.1 空间离散方法 以二维标量守恒律方程
 | (6) |
为例，简单介绍间断有限元方法。
式中：u=u(x, t)为守恒变量，x=(x, y)为空间坐标向量，t为时间；f和g为通量向量的分量；u(x, 0)为变量u在t=0时刻的取值；u0(x)为初值函数，其决定了u在t=0时刻的取值；下标t、x、y分别表示变量对t、x、y的偏导数。
用Th代表计算区域Ω的一套网格划分，则式(6)的间断有限元解所在的空间为
 | (7) |
式中：P(Ωk)为定义在单元Ωk上、阶数不大于p的全体多项式的集合。本文选择定义于Ωk内均布结点上的Lagrange插值多项式lik(x)(i=1, 2, …, N)作为P(Ωk)的基函数[14]。在单元Ωk内，以
 | (8) |
逼近u，将式(8)代入式(6)，然后用测试函数lik(x)与之相乘，利用分部积分，得到
 | (9) |
式中：uhknb为单元Ωk的邻居单元Ωknb上的数值解；nk为?Ωk上的外单位法向量；H^为数值通量，可理解为对(f, g)·nk的近似，有多种类型可供选择[15]; dS为沿单元边界进行线积分的长度微元。依照式(9)，对全部K个单元的全体测试函数列式，可以得到下述常微分方程组：
 | (10) |
式中：
 | (11) |
对于
 | (12) |
Rk为与Uk维数相同的向量，其第i个分量为测试函数为lik(x)时，式(9)等号右边的表达式。
1.2.2 时间离散方法 对于式(10)，选择时间步长可变的向后差分公式(BDF)作为时间离散方法。该方法属于隐式方法，可通过增加级数达到高阶精度，且每个时间步内的计算量不随级数的增加而增加。其形式为
 | (13) |
式中：Δt=tn－tn－1；Un－l代表t=tn－l时刻的数值解；q为方法的级数，当q=3时，方法具有3阶精度; 系数αq, l的计算方法如下[16]：
 | (14) |
1.2.3 弯曲边界处理 采用高阶间断有限元方法求解含有弯曲物面的问题时，通常需要在物面附近布置曲边单元，目的是防止边界表达不当对格式精度的破坏[17]。考虑到常见的网格生成器只能生成直边网格，首先，本文在局部以圆弧逼近弯曲边界[18]，从而获得用于定义曲边单元的结点坐标。接下来，利用等参变换[19]，基于弯曲边界上的结点坐标，本文建立起参考单元和曲边单元的映射关系。至于利用该映射关系，在参考单元上计算基于弯曲单元的体积分、面积分的公式，可参见文献[16]第8.6.2节。
2 方程组求解 2.1 Jacobi矩阵的构造 采用牛顿法解非线性方程组(13)时，需要计算等号左边的表达式对解向量U的Jacobi矩阵。由式(10)、式(13)，主要的工作在于得到R对U的Jacobi矩阵。
对于标量方程(6)，记上述R对U的Jacobi矩阵为J，则J含有(N×K)2个元素。将J按行、列划分为K2个子阵，则序号为(m, n)的子阵记为Jm, n，其对应于Rm对Un的偏导数。由式(9)，除n=m、单元Ωn与Ωm相邻2种情况以外，Jm, n内的元素均为0。
对于Jk, k，其第(i, j)个元素为
 | (15) |
至于Jk, knb，其第(i, j)个元素为
 | (16) |
对于二维Euler方程组：
 | (17) |
式中：
 | (18) |
状态方程为
 | (19) |
式中：E为动能和内能之和；μ和v分别为x和y方向速度；p为压强；ρ为密度；γ为比热容比，本文取值为1.4。
经过间断有限元空间离散后，相应的Jacobi矩阵含有(N×K×4)2个元素，可以被分为4×4个子阵，其中序号为(r, s)的子阵对应于fr、gr对ws的偏导数；该子阵的计算方法与前述标量方程Jacobi矩阵的计算方法相同。
见式(15)、式(16)，计算Jacobi矩阵涉及到数值通量对守恒变量的导数。大多数数值通量因含有取绝对值、取最大值等函数，本质上不可微。实际应用中，经常采用线性化手段，给出数值通量名义上的“导数”。本文选择Vijayasundaram型通量[20]
 | (20) |
作为区域内部的单元界面上的数值通量。在计算Jacobi矩阵时，令P+、P－分别作为数值通量对本单元Ωk、邻居单元Ωknb的守恒变量的导数。
2.2 解线性方程组 本文采用牛顿法求解源于隐式时间离散的非线性方程组。每个牛顿迭代步内，待求解的线性方程组为
 | (21) |
式中：
 | (22) |
类似于式(13)，We, n－l为二维Euler方程组内第e个方程在t=tn－l时刻的数值解向量；Qe为相应的间断有限元空间离散算子；其余参数的含义与式(13)、式(14)相同。
式(21)中，系数矩阵HW的表达式为
 | (23) |
式中：Jall为由4×4个子阵组成的Jacobi矩阵；I为单位矩阵；Mall为分块矩阵，除4个对角块均为式(11)中的M外，其余的元素均为0。HW阶数高(N×K×4)、稀疏性强，使得式(21)适于采用GMRES方法求解。
GMRES方法是一种迭代方法，理论上已经证明，线性方程组系数矩阵的条件数越小，迭代的收敛速度越快[10]；鉴于此，通常采用预处理技术加快方程组的求解速度。以式(1)为例，设A为稀疏矩阵，实际应用中，并非直接以GMRES方法求解式(1)，而是先构造一个接近于A、且便于求逆的矩阵C，然后采用GMRES方法求解下述线性方程组：
 | (24) |
C越接近A，C－1A的条件数就越小，GMRES迭代的收敛速度就越快。
GMRES方法并不需要精确的A－1。如果能够以较少的计算花费得到系数矩阵A的一个近似的LU分解：
 | (25) |
利用矩阵Li、Ui为三角形、便于求逆的特点，就可以快速构建出预处理矩阵：
 | (26) |
上述方法称为不完全LU分解(Incomplete LU decomposition，ILU)，是一种常用的构建GMRES预处理矩阵的方法。通过调整算法参数，可以使Li、Ui具有与A同等程度的稀疏性(称为无填充ILU)[21]。
不完全LU分解克服了精确LU分解占用存储空间大、计算缓慢的问题，但代价是预处理矩阵对系数矩阵的偏离。从GMRES方法的角度考虑，Ui-1Li-1越接近A-1，GMRES达到收敛的迭代步数就越少；事实上，如果用精确的LU分解构造预处理矩阵，只需迭代1步就能达到收敛。反之，预处理矩阵对系数矩阵偏离得越远，所需的GMRES迭代步数就越多。以计算资源作为评判标准，预处理矩阵的计算和GMRES迭代之间似乎存在一种矛盾。
回顾矩阵带宽缩减在以直接法求解线性方程组中的应用，同样是精确的LU分解，分解原始矩阵的速度明显慢于分解带宽缩减矩阵的速度，且后者的分解结果具有更强的稀疏性，占用的存储空间更小。本文将其应用于以GMRES迭代方法求解线性方程组。具体地，对于式(21)，首先按照RCM方法，依据HW的稀疏性构造队列Y，然后用带预处理矩阵的GMRES方法(预处理矩阵采用无填充ILU方法计算)求解线性方程组：
 | (27) |
式中：
 | (28) |
解得
3 数值实验 本节采用隐式间断有限元方法求解3个典型的空气动力学问题。求解时，为了检验矩阵带宽缩减技术的效果，对该技术应用前、后的GMRES迭代过程进行了比较。
全部测试问题均为二维、无黏、定常流动，控制方程为式(17)~式(19)；采用时间相关法，以物理时间充分大时的解作为定常问题的解；空间离散均选用分片3阶多项式，时间离散均选择3级3阶向后差分；时间离散的步长按式(29)计算[2, 3, 16]：
 | (29) |
式中：CFL为量纲为1的参数，用于调节时间步长的大小，对于显式格式，CFL≤1；|Ωk|为单元Ωk的面积；λ|Γ为最大特征速度λ在单元边Γ上的取值；|Γ|为单元的边的长度；v=(u, v)为流体微团的速度矢量；c为声速。
当计算区域边界包含物体表面时，当地的数值解必须满足无穿透要求。见式(18)，对于与物面重合的单元界面，本文将条件v·nk=0代入(fe, ge)·nk，e=1~4，以修改通量函数的方式满足该要求。计算Jacobi矩阵时，直接求修改后的通量对守恒变量的偏导数即可。
3.1 圆柱绕流 圆柱直径为1，计算区域外边界到圆心距离为20；计算区域被划分为512个三角形网格，壁面附近的网格分布如图 3所示；外边界上采用特征边界条件，来流马赫数Ma∞=0.38，迎角α=0°；采用来流条件对问题进行初始化。
 |
| 图 3 圆柱附近的网格 Fig. 3 Grids around a cylinder |
| 图选项 |
图 4展示了原始的Jacobi矩阵HW，图 5则是经RCM方法缩减带宽的Jacobi矩阵
 |
| 图 4 带宽未经缩减的Jacobi矩阵 Fig. 4 Jacobian matrix without bandwidth reduction |
| 图选项 |
 |
| 图 5 带宽经过缩减的Jacobi矩阵 Fig. 5 Jacobian matrix with bandwidth reduction |
| 图选项 |
为了考察矩阵带宽缩减技术对GMRES迭代过程的影响，进行了对比实验：将迭代的收敛标准设定为绝对误差小于1×10-5，迭代步数上限定为100，以无填充ILU方法构造预处理矩阵；以残差达到1.6×10－4(对应物理时间tn=25.25)后的第1个牛顿迭代步作为观察区间，记录带宽缩减前、后的GMRES迭代过程，结果如表 1和表 2所示。
表 1 GMRES的迭代步数(圆柱绕流) Table 1 Iteration steps of GMRES(flow around a cylinder)
| CFL | 带宽未经缩减 | 带宽经过缩减 |
| 5 | 6 | 4 |
| 10 | 11 | 5 |
| 20 | 15 | 6 |
表选项
表 2 预处理矩阵与系数矩阵之差的Frobenius范数(圆柱绕流) Table 2 Frobenius norm for difference between preconditioner and coefficient matrix (flow around a cylinder)
| CFL | 带宽未经缩减 | 带宽经过缩减 |
| 5 | 40 | 34 |
| 10 | 205 | 85 |
| 20 | 40 157 | 203 |
表选项
表 1中，对于相同的时间步长，若线性方程组系数矩阵的带宽未经过缩减，则相应的GMRES收敛步数明显多于带宽经过缩减的情况。
表 2考察了带宽缩减前、后预处理矩阵与线性方程组系数矩阵的接近程度，考察标准为
 | (30) |
式中：A为线性方程组的系数矩阵；LiUi为经由无填充ILU方法得到的预处理矩阵；下标F表示矩阵的Frobenius范数。通过比较可以发现，带宽缩减后，预处理矩阵与系数矩阵的接近程度显著优于带宽缩减前的情况，这就解释了表 1中缩减带宽的矩阵所需的收敛步数更少的现象。
圆柱绕流问题壁面附近的马赫数等值线如图 6所示，计算结果具有良好的对称性，与文献[18]中的结果符合得很好。
 |
| 图 6 圆柱周围的马赫数等值线(Ma∞=0.38, α=0°, ΔMa=0.05) Fig. 6 Mach number contours around a cylinder (Ma∞=0.38, α=0°, ΔMa=0.05) |
| 图选项 |
3.2 NACA0012翼型绕流 翼型弦长为1，计算区域外边界到翼型表面的距离在翼型弦长的15~20倍之间；计算区域被划分为1 313个三角形网格，壁面附近的网格分布如图 7所示；外边界上采用特征边界条件，来流马赫数Ma∞=0.63，迎角α=2°，采用来流条件对问题进行初始化。
 |
| 图 7 NACA0012翼型附近的网格 Fig. 7 Grids around NACA0012 airfoil |
| 图选项 |
与3.1节类似，对于本例，缩减带宽前，Jacobi矩阵内的非零元素集中在各个子块的对角线附近，经RCM方法缩减带宽后，全部非零元素都集中在Jacobi矩阵的主对角线附近。为了节约篇幅，本例未列出带宽缩减前、后的Jacobi矩阵。
本例同样进行了类似于3.1节的对比实验：GMRES迭代的收敛标准、迭代步数上限及预处理矩阵的构造方法与3.1节相同；以残差达到1.1×10-4(对应物理时间tn=6.82)后的第1个牛顿迭代步作为观察区间，记录带宽缩减前、后的GMRES迭代过程，结果见表 3、表 4。
表 3 GMRES的迭代步数(NACA0012翼型绕流) Table 3 Iteration steps of GMRES(flow around NACA0012 airfoil)
| CFL | 带宽未经缩减 | 带宽经过缩减 |
| 5 | 5 | 4 |
| 10 | 14 | 5 |
| 20 | 未收敛 | 6 |
表选项
表 4 预处理矩阵与系数矩阵之差的Frobenius范数(NACA0012翼型绕流) Table 4 Frobenius norm for difference between preconditioner and coefficient matrix(flow around NACA0012)
| CFL | 带宽未经缩减 | 带宽经过缩减 |
| 5 | 42 | 33 |
| 10 | 280 | 86 |
| 20 | 334 511 | 205 |
表选项
实验结论与3.1节类似。注意到在表 3中，当CFL=20时，未缩减Jacobi矩阵的带宽，会出现在指定迭代步数内无法收敛的情况。另一方面，实际计算本例时，如果应用了矩阵带宽缩减技术，即使令CFL=80，GMRES方法依然可以在指定迭代步数内达到收敛。
计算结果如图 8、图 9所示。其中图 8是翼型表面的压力系数Cp分布，图 9是翼型附近的马赫数等值线。可以看出，计算结果与文献[18]中的结果符合得很好。
 |
| 图 8 NACA0012翼型表面的压力系数(Ma∞=0.63, α=2°) Fig. 8 Pressure coefficient on surface of NACA0012 airfoil (Ma∞=0.63, α=2°) |
| 图选项 |
 |
| 图 9 NACA0012翼型周围的马赫数等值线(Ma∞=0.63, α=2°, ΔMa=0.03) Fig. 9 Mach number contours around NACA0012 airfoil (Ma∞=0.63, α=2°, ΔMa=0.03) |
| 图选项 |
3.3 RAE2822翼型绕流 翼型弦长为1，计算区域外边界到翼型表面的距离在翼型弦长的15~20倍之间；计算区域被划分为1 514个三角形网格，壁面附近的网格分布如图 10所示；外边界上采用特征边界条件，来流马赫数Ma∞=0.4，迎角α=2.79°；采用来流条件对问题进行初始化。
 |
| 图 10 RAE2822翼型附近的网格 Fig. 10 Grids around RAE2822 airfoil |
| 图选项 |
对于本例，带宽缩减前、后的Jacobi矩阵内的非零元素的分布情况与前面两例类似。
本例同样进行了类似于3.1节、3.2节的对比实验：GMRES迭代的收敛标准、迭代步数上限及预处理矩阵的构造方法与3.1节相同；以残差达到1.7×10－4(对应物理时间tn=2.21)后的第1个牛顿迭代步作为观察区间，记录带宽缩减前、后的GMRES迭代过程，结果如表 5、表 6所示。对比实验的结论与3.1节、3.2节类似。计算结果如图 11、图 12所示。其中图 11是翼型表面的压力系数分布，图 12是翼型附近的马赫数等值线。可以看出，计算结果与文献[22]中的结果较为吻合。
表 5 GMRES的迭代步数(RAE2822翼型绕流) Table 5 Iteration steps of GMRES (flow around RAE2822 airfoil)
| CFL | 带宽未经缩减 | 带宽经过缩减 |
| 10 | 7 | 5 |
| 20 | 10 | 5 |
| 40 | 33 | 6 |
表选项
表 6 预处理矩阵与系数矩阵之差的Frobenius范数(RAE2822翼型绕流) Table 6 Frobenius norm for difference between preconditioner and coefficient matrix (flow around RAE2822 airfoil)
| CFL | 带宽未经缩减 | 带宽经过缩减 |
| 10 | 70 | 63 |
| 20 | 231 | 155 |
| 40 | 349 645 | 369 |
表选项
 |
| 图 11 RAE2822翼型表面的压力系数(Ma∞=0.4, α=2.79°) Fig. 11 Pressure coefficient on surface of RAE2822 airfoil (Ma∞=0.4, α=2.79°) |
| 图选项 |
 |
| 图 12 RAE2822翼型周围的马赫数等值线(Ma∞=0.4, α=2.79°, ΔMa=0.015) Fig. 12 Mach number contours around RAE2822 airfoil (Ma∞=0.4, α=2.79°, ΔMa=0.015) |
| 图选项 |
4 结论 本文将矩阵带宽缩减技术应用于以隐式间断有限元方法求解空气动力学问题当中，得到：
1) 针对基于非结构网格的间断有限元空间离散，构造了相应的Jacobi矩阵；该矩阵具有阶数高、稀疏性强、数值非对称的特点，相应的线性方程组适于采用GMRES方法迭代求解。
2) GMRES方法的迭代速度与预处理矩阵的品质有关；通过利用矩阵带宽缩减技术调整系数矩阵(与上述Jacobi矩阵密切相关)的非零元素分布，预处理矩阵与系数矩阵的接近程度显著提升，加速了GMRES方法的收敛，增大了内存有限的前提下、计算能够达到的最大时间步长。
参考文献
| [1] | 吕宏强, 张涛, 孙强, 等. 间断伽辽金法在可压缩流数值模拟中的应用研究综述[J]. 空气动力学学报, 2016, 35(4): 455-471. LYU H Q, ZHANG T, SUN Q, et al. Applications of discontinuous Galerkin method in numerical simulations of compressible flows:A review[J]. Acta Aerodynamica Sinica, 2016, 35(4): 455-471. (in Chinese) |
| [2] | COCKBURN B, SHU C W. The Runge-Kutta discontinuous Galerkin method for conservation laws V:Multidimensional systems[J]. Journal of Computational Physics, 1997, 141(2): 199-224. |
| [3] | COCKBURN B, SHU C W. Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[J]. Journal of Scientific Computing, 2001, 16(3): 173-261. DOI:10.1023/A:1012873910884 |
| [4] | FU G S, SHU C W. A new trouble-cell indicator for discontinuous Galerkin methods for hyperbolic conservation laws[J]. Journal of Computational Physics, 2017, 347: 305-327. DOI:10.1016/j.jcp.2017.06.046 |
| [5] | 颜庆津. 数值分析[M]. 4版. 北京: 北京航空航天大学出版社, 2012: 202-203. YAN Q J. Numerical analysis[M]. 4th ed. Beijing: Beihang University Press, 2012: 202-203. (in Chinese) |
| [6] | KANEVSKY A, CARPENTER M H, GOTTLIEB D, et al. Application of implicit-explicit high order Runge-Kutta methods to discontinuous-Galerkin schemes[J]. Journal of Computational Physics, 2007, 225(2): 1753-1781. DOI:10.1016/j.jcp.2007.02.021 |
| [7] | 郝海兵, 张强, 杨永, 等. 基于LU-SGS迭代的DGM隐式方法研究[J]. 西北工业大学学报, 2014, 32(3): 346-350. HAO H B, ZHANG Q, YANG Y, et al. An implicit scheme of discontinuous Galerkin method(DGM) based on LU-SGS iterative method[J]. Journal of Northwestern Polytechnical University, 2014, 32(3): 346-350. DOI:10.3969/j.issn.1000-2758.2014.03.003 (in Chinese) |
| [8] | BASSI F, BOTTI L, COLOMBO A, et al. Linearly implicit Rosenbrock-type Runge-Kutta schemes applied to the discontinuous Galerkin solution of compressible and incompressible unsteady flows[J]. Computers & Fluids, 2015, 118: 305-320. |
| [9] | SAAD Y, SCHULTZ M H. GMRES:A generalized minimal residual algorithm for solving nonsymmetric linear systems[J]. SIAM Journal on Scientific and Statistical Computing, 1986, 7(3): 856-869. DOI:10.1137/0907058 |
| [10] | 马昌凤, 柯艺芬, 唐嘉, 等. 数值线性代数与算法(MATLAB版)[M]. 北京: 国防工业出版社, 2017: 152-176. MA C F, KE Y F, TANG J, et al. Numerical linear algebra and algorithms(MATLAB edition)[M]. Beijing: National Defense Industry Press, 2017: 162-176. (in Chinese) |
| [11] | 华伯浩. 大型有限元程序系统的网格结点编码优化算法[J]. 应用数学和力学, 1981, 2(2): 207-218. HUA B H. The network node relabeling optimum algorithms for large finite element program system[J]. Applied Mathematics and Mechanics, 1981, 2(2): 207-218. (in Chinese) |
| [12] | 王家林. 有限元模型中自由度层次的带宽优化算法[J]. 西南交通大学学报, 2009, 44(2): 186-189. WANG J L. Bandwidth optimization algorithm of finite element models at level of degree of freedom[J]. Journal of Southwest Jiaotong University, 2009, 44(2): 186-189. DOI:10.3969/j.issn.0258-2724.2009.02.008 (in Chinese) |
| [13] | LIU W H, SHERMAN A H. Comparative analysis of the Cuthill-Mckee and the reverse Cuthill-Mckee ordering algorithms for sparse matrices[J]. SIAM Journal on Numerical Analysis, 1976, 13(2): 198-213. DOI:10.1137/0713020 |
| [14] | HESTHAVEN J S, WARBURTON T. Nodal discontinuous Galerkin methods:Algorithm, analysis, and applications[M]. New York: Springer Science and Business Media, 2008: 297-299. |
| [15] | QIU J X, KHOO B C, SHU C W. A numerical study for the performance of the Runge-Kutta discontinuous Galerkin method based on different numerical fluxes[J]. Journal of Computational Physics, 2006, 212: 540-565. DOI:10.1016/j.jcp.2005.07.011 |
| [16] | DOLEJSI V, FEISTAUER M. Discontinuous Galerkin method:Analysis and applications to compressible flow[M]. Berlin: Springer International Publishing, 2015: 423-424, 434-438. |
| [17] | BASSI F, REBAY S. High-order accurate discontinuous finite element solution of the 2D Euler equations[J]. Journal of Computational Physics, 1997, 138(2): 251-285. DOI:10.1006/jcph.1997.5454 |
| [18] | KRIVODONOVA L, BERGER M. High-order accurate implementation of solid wall boundary conditions in curved geometries[J]. Journal of Computational Physics, 2006, 211(2): 492-512. DOI:10.1016/j.jcp.2005.05.029 |
| [19] | 王勖成. 有限单元法[M]. 北京: 清华大学出版社, 2003: 131-132. WANG X C. Finite element method[M]. Beijing: Tsinghua University Press, 2003: 131-132. (in Chinese) |
| [20] | DOLEJSI V, FEISTAUER M. A semi-implicit discontinuous Galerkin finite element method for the numerical solution of inviscid compressible flow[J]. Journal of Computational Physics, 2004, 198(2): 727-746. DOI:10.1016/j.jcp.2004.01.023 |
| [21] | 方田.Krylov子空间法及预处理技术在CFD中的应用[D].武汉: 武汉理工大学, 2013: 45-52. FANG T.Krylov sub-space iterative method with preconditioning and its application in numerical simulation for CFD[D].Wuhan: Wuhan University of Technology, 2013: 45-52(in Chinese). |
| [22] | BURGESS N K, MARVRIPLIS D J. Robust computation of turbulent flow using discontinuous Galerkin method:AIAA-2012-0457[J]. Reston:AIAA, 2012, 25-27. |
