求解欧拉方程的预处理隐式无网格算法

曹 骋，陈红全

(南京航空航天大学 航空学院，南京 210016)

随着现代计算机技术和数值分析的发展，计算流体力学(简称CFD)已取得了长足的进步，涌现出了可用于先进飞行器设计的各类CFD计算软件[1-5]. 这类软件大多选用网格单元进行计算区域的离散，因此，这类方法可统称为“网格方法”，通常在流场数值模拟前都要求先对计算区域进行网格生成. 由于网格拓扑的约束，要对包含飞机全机等复杂几何外形的计算区域生成合适的整体网格并非易事，尤其要求刻画出多体小间隙等几何细节时. 因此，如何局部或完全摆脱网格的约束问题，成为人们开始关注和研究的热点[6-8]，出现了一类意在完全打破网格拓扑限制的算法，无网格算法. 这类算法对于处置和刻画复杂几何外形细节将更具灵活性，其理由主要有：计算区域是用无网格点云点布点离散的，具有灵活性和兼容性，既可按要求直接布点，也可利用已有的结构或非结构网格的网格点；无网格点上的空间导数近似仅依赖于当地点云点上的信息，而点云点无须连成网格. 许多学者正是受这些重要特征的驱动，开始对其研究，各具特色的无网格算法相继涌现[9-14]. 就作者熟悉的空气动力学领域而言，最引人注目的是早期Batina[6]的工作. 他采用中心差分格式，发展了显式无网格算法，成功地求解了带有激波的可压缩绕流. 之后，Morinishi[11]采用虚拟中点(点云中心点与卫星点之间)迎风处理策略和加权最小二乘空间导数逼近，成功地发展了隐式无网格算法. 近年来，成功应用无网格算法的算例[15-18]已不难在文献中找到，借助于现代图形处理器(GPU)等并行加速技术，算法的收敛速度也进一步得到了强化[19-22]. 但同时可以注意到，上述无网格算法大多是基于可压缩流发展的，不能直接拓展到模拟几乎不可压的低马赫数流.

本文致力于发展出一种带预处理的隐式无网格算法，既能模拟可压跨音速流，又能模拟几乎不可压的低马赫数流. 众所周知，在低马赫数下，特征型欧拉方程的3个不同特征值代表的波速存在巨大差异，会导致原控制方程产生系统病态，即所谓的“刚性”问题[23]. 而这反过来也会导致迭代计算收敛困难，数值模拟解往往偏离实际物理解. 为此，本文先沿用网格算法中著名的Weiss-Smith型预处理矩阵[24]，对无网格半离散型欧拉方程进行矩阵预处理，以期缩减方程特征值在低马赫数条件下的差异，改善上述所谓的“刚性”问题；接着给出了无网格点上受预处理影响的人工耗散项及远场边界条件等实施细节，并通过无网格点排序与点云分割，实现了LU-SGS算法[25]的时间推进求解，形成了预处理隐式无网格算法. 最后，先选用典型翼型算例进行了数值模拟，并与文献或实验结果进行了验证比较；接着给出了机翼及翼身组合体等三维跨音速和低马赫数绕流算例. 算例表明所发展的隐式算法如预期比相应的显式算法收敛更快，已从单纯模拟可压缩流动拓展到模拟几乎不可压的低马赫数流.

1 控制方程

在直角坐标系(x,y,z)上，与时间t相关的三维欧拉方程的守恒形式为[26]

(1)

其中W为守恒矢变量，E，F和G为矢通量，可分别写为

W=[ρ,ρu,ρv,ρw,ρE]T,

(2)

E=[ρu,ρu2+p,ρuv,ρuw,ρuH]T,

(3)

F=[ρv,ρvu,ρv2+p,ρvw,ρvH]T,

(4)

G=[ρw,ρwu,ρwv,ρw2+p,ρvH]T.

(5)

式中：ρ，p，E和H分别代表了密度、压强、单位质量的总能和总焓；u，v和w为速度分量. 完全气体满足状态方程，有如下关系式成立：

(6)

ρH=ρE+p.

(7)

其中，γ为流体的比热比，对于空气而言，γ=1.4. 通过求解上述方程组寻求满足∂W/∂t=0的定常解.

2 无网格空间导数近似

采用无网格算法，计算域先需布点离散. 布点离散后，即可形成由每一点临近点构成的当地点云结构[6,8]. 记C(i)为第i点的点云，其构成参见图1. 图中，点云的中心点即为i点，点云的其余点被称作卫星点.

图1 三维点云C(i)示意图Fig.1 Three-dimensional cloud of points for C(i)

无网格空间导数近似是在计算域离散点上进行的. 基于图1所示的点云结构，任一可微函数f的空间导数可近似为[8,27]

(8)

其中fik为第k个卫星点与中心点i连线中点处的近似值. 显然，系数αik，βik和γik满足

(9)

上述空间导数近似式中的系数，可用加权的最小二乘曲面逼近来求得[8,19].

3 预处理隐式无网格算法

3.1 通量计算

应用(8)式，欧拉方程的通量可近似写为

(10)

其中i和k连线中点上的数值通量Lik具体可写为

(11)

这里Uik与逆变速度定义类似，可定义为

Uik=αikuik+βikvik+γikwik.

(12)

假定中心点i和卫星点k连线的中点处存在一个虚拟的交界面，那么交界面上的通量Lik可以通过选择合适的格式，从i点和k点处的守恒变量计算确定. 采用中心格式，交界面处的守恒变量可简单地写为

(13)

显然，交界面上Lik可由Wik计算得到.

必须指出，上述中差计算形成的格式与网格算法中采用的中心差分格式类似，迭代计算会出现不稳定. 一般在激波这样的强间断附近会产生解的非物理振荡，这种振荡也可能因所谓的奇偶失联而发生在任何无网格点上. 因此，本文沿用Jameson[28]在网格算法中的做法，在离散的方程中添加人工耗散项Disi以保证算法稳定. 基于无网格点云结构，人工耗散项可写为

(14)

其中ε(2)和ε(4)为可调系数，λik则为雅克比矩阵Aik=∂Lik/∂Wik的谱半径，其形式为

(15)

应用(10)式，添加上人工耗散项后，得到方程(1)的半离散形式，具体可简写为

(16)

3.2 预处理

可以注意到，上述交界面上欧拉方程的3个不同的特征值为

(17)

在低马赫数下，如前述，这3个特征值代表的波速会存在巨大差异，雅克比矩阵Aik的条件数过大，会导致迭代计算收敛困难. 为此，本文考虑对方程(16)进行矩阵预处理以减小矩阵条件数. 从网格算法中可知，适用于低马赫数流动的预处理矩阵[24,29-32]可有多种选择，其中较为著名的是Weiss-Smith型预处理矩阵[24]. 对守恒型方程，该矩阵具体可写为

Γ=

(18)

其中：

(19)

(20)

(21)

这里的Ma和Ma分别代表了当地和来流马赫数.

应用(18)式进行矩阵预处理，方程(16)可改写为

(22)

(23)

因此，预处理以后人工耗散项可改写为

(24)

本文将采用LU-SGS隐式算法，对预处理以后的半离散方程(22)进行时间离散求解.

3.3 时间离散

依据Jameson和Yoon的做法[25]，对预处理以后的数值通量引入如下线性化处理:

(25)

这里通量雅克比矩阵Bik按文献[25]中的分裂方式定义为

(26)

其中I为单位矩阵. 那么方程(22)的时间后差隐式格式可以写为

(27)

其中Δti为点云中心点i的时间步长，R为空间导数用点云近似的残值，ΔW定义为

ΔW=Wn+1-Wn+1

(28)

因为系数αik，βik和γik具有式(9)的特性，雅克比矩阵Bik也应该具有如下特性:

(29)

因而，隐式格式(27)可退化为

(30)

式(30)用LU-SGS算法[25]求解，格式可分为如下两步:

(31)

(32)

其中，L(i)和U(i)为点云C(i)卫星点的LU分割.Di定义为

(33)

实际计算中，沿用Sharov和Nakahashi[33]的做法，涉及的矩阵运算Aik(Wk)ΔWk由通量的增量ΔLik(Wk)来简化计算，减少工作量.

必须指出的是，L(i)和U(i)如果分割不当会导致算法退化为运算较慢的Jacobi迭代[33]. 为避免该问题，文献[33]依据网格间的连接关系先进行网格分层编号，再进行LU分割. 层号比当前网格层号小的相邻网格单元被归入L(i)，反之归入U(i). 对于无网格算法，点云点间不存在网格拓扑这样的刚性连接关系，为此，本文利用卫星点与中心点的虚拟连接关系来构造类似的分层结构. 作者曾提出一种多层点云的重排算法，优化了GPU的内存访问模式， 实现了GPU并行无网格算法的进一步加速[22]. 本文沿用这一方法对无网格点进行类似的分层编号重排. 分层后，对点云C(i)的卫星点进行LU分割. 将层号比中心点i层号小的卫星点分割纳入L(i)，反之则纳入U(i). 本文将这一分割做法用于了后续的算例运算.

3.4 边界条件

对于物面边界条件，无粘条件下满足的无穿透边界条件[6]在预处理前后都一样. 而远场基于特征分析的无反射边界条件，预处理后由于特征值的改变，必须作相应的修正. Turkel[34]早期给出了精确的基于特征分析的预处理远场边界条件，之后又对其进行了简化，用于处理几乎不可压流动. 本文沿用Turkel的简化处理[32]，预处理远场边界条件是结合无网格布点离散灵活的特点实施的. 为此，本文把远场边界面处理成临近内点的镜面，每一离镜面最近的内点对应一个虚拟的镜像点，纳入该点的点云中. 那么，虚拟点上的物理量可直接根据Turkel的方法确定.

以远场边界入流为例，虚拟点处的原始变量可取为

ρg=ρ,ug=u,vg=v,wg=w,pg=pi.

(34)

这里的下标g为虚拟点，而下标i为内部流场点，下标代表了自由流. 注意，这里g是i的镜像点.

需要指出的是，位于远场边界外引入的镜像虚拟点是在计算域布点离散过程中生成的. 这样做，相当于远场边界条件位于虚拟点和对应内部点的中点处，而所有位于计算域内与远场边界条件相关的临近边界点都可以按照内部点的方式处理. 可以看出，这样处理的远场边界条件实施更为简单方便.

4 算例与分析

本文用上述预处理隐式无网格算法，先对NACA 0012对称翼型和AGARD三段翼型进行了低马赫数绕流的数值模拟，通过与文献或实验结果的对比分析，验证所提算法的计算效率和准确性，展示隐式无网格算法的快速收敛特性；接着给出了机翼及翼身组合体等跨音速和低马赫数两类绕流的算例，展示算法通过预处理对这两类流动模拟的兼容能力.

4.1 NACA 0012翼型低马赫数绕流

本文选用NACA 0012翼型几乎不可压的低马赫数绕流对发展的算法进行了考核计算验证. 如图2所示，计算域用3343个无网格点进行离散. 先用马赫数Ma=0.001，攻角为0°的绕流模拟展示了预处理对收敛历程的影响(见图3). 图中可看出，发展的预处理隐式算法，在如此小的马赫数下仍能快速收敛，图中同时给出了不加预处理的收敛困难的情形供比较. 对应的翼型表面压强系数分布(见图4)取得了与文献[35]和实验[36]一致的结果. 图中横坐标X/C代表离前缘的相对距离，C为对应翼型的弦长，纵坐标Cp则为压强系数. 图5则给出了对应的马赫数等值线和云图. 图中反映的流场对称性与零攻角对称翼型绕流的物理特性吻合. 接着，通过不同马赫数的数值模拟展示了马赫数对算法收敛历程的影响. 从图6中可看出，模拟的4个来流马赫数(0.3，0.1，0.05，0.01)都能很好收敛. 图中同时给出了显式预处理算法[23]的收敛历程供比较. 必须指出，在几乎不可压的低马赫数下，很有必要使算法具有良好收敛特性. 这一点从捕捉的流场解可以看出. 图7给出了隐式无网格算法带与不带预处理计算获得的流场密度等值线(Ma=0.05)比较. 图中可看出，预处理后流场解更加合理，表现为等值线更加光滑.

图2 NACA 0012翼型周围点分布Fig.2 Points around NACA 0012 airfoil

图3 预处理对收敛历程的影响Fig.3 Effects of preconditioning on convergence histories

图4 翼型表面压强系数分布(Ma=0.001)Fig.4 Pressure coefficient distribution(Ma=0.001)

图5 马赫数等值线及云图(Ma=0.001)Fig.5 Mach number contours (Ma=0.001)

(a) 预处理隐式

(b) 预处理显式图6 马赫数对收敛历程的影响Fig.6 Effects of Mach numbers on convergence histories

(a) 带预处理

(b) 不带预处理图7 隐式无网格算法密度等值线比较(Ma=0.05)

Fig.7 Density contours of implicit gridless methods(Ma=0.05)

4.2 AGARD三段翼型低马赫数绕流

这里给出外形更为复杂的AGARD三段翼型[37]几乎不可压的低马赫数绕流的计算结果. 几何外形可从图8中看出. 计算域用5865个无网格点进行离散. 来流马赫数按实验条件[37]取为0.197，迎角取为4.01°. 如图9收敛历程所示，发展的隐式预处理无网格算法再次呈现出快速收敛特性，图中同时给出了隐式不带预处理、显式带与不带预处理无网格算法的收敛历程供比较. 对应的翼面压强系数分布(见图10)与实验[37]吻合的较好，这一点可从主翼上捕捉的激波位置和强度看出. 图11则给出了对应的马赫数等值线和云图，展示出了多段翼型的流动特征.

图8 多段翼型周围点分布 Fig.8 Points around multi-element airfoil

图9 不同无网格算法收敛历程比较Fig. 9 Convergence histories of different methods

图10 表面压强系数分布Fig.10 Pressure coefficient distribution

图11 马赫数等值线和云图Fig.11 Mach number contours

4.3 M6机翼的跨音速和几乎不可压绕流

M6机翼[6,38]经常被用来考核发展的算法. 本文先进行了绕M6机翼跨音速流动的数值模拟. 计算采用了54320个无网格点(见图12)，来流马赫数为0.84，攻角为3.06°. 典型机翼翼剖面的压强系数分布见图13. 图中η代表机翼翼剖面展向相对位置，翼尖η为1，翼根η为0. 计算结果与文献[6]和实验[39]进行了比较. 如预期，带预处理的隐式无网格算法仍能用于该跨音速绕流计算，表现为计算结果与不带预处理的结果十分吻合.

图12 M6 机翼周围点分布Fig.12 Points around ONERA M6 wing

(a)η=0.65

(b)η=0.9图13 机翼表面压强系数分布 (Ma=0.84)Fig.13 Pressure coefficient distribution(Ma=0.84)

接着，沿用文献[38]的做法，进行了几乎不可压马赫数绕流(Ma=0.01)的考核运算. 图14给出了计算获得的η=0.8处的翼剖面压强系数分布，从图中可以看出，发展的预处理隐式无网格算法与预处理显式算法一样，计算结果都与文献[38]结果十分吻合. 比较而言，隐式预处理无网格算法较显式收敛特性更优(见图15). 图16则给出了对应的上翼面压强系数等值线云图，可以看出捕捉的等值线光滑有序，一定程度上反映出M6机翼上翼面几乎不可压的低速流动特征.

图14 机翼表面压强系数分布(Ma=0.01，η=0.8)Fig.14 Pressure coefficient distribution (Ma=0.01, η=0.8)

图15 收敛历程(Ma=0.01)Fig.15 Convergence histories(Ma=0.01)

图16 机翼上表面压强系数等值线及云图(Ma=0.01)

Fig.16 Pressure coefficient contours on the upper surface of the wing (Ma=0.01)

4.4 DLR-F4翼身组合体的低马赫数绕流

面向实际应用，这里给出DLR-F4组合体低马赫数绕流的演示算例. 众所周知，马赫数0.3在工程中常作为判断一个流动是可压与不可压的参考马赫数. 必须指出，按照可压和不可压的这一区分标准，对于来流马赫数为0.3的翼身组合体绕流，其流场存在驻点附近等局部低马赫数区和机翼上翼面等局部高马赫数区，一定程度上可理解为可压和不可压并存的混合流场. 不难看出，随着来流马赫数的降低，流场中局部低马赫数的不可压区在整个流场中的占比将逐渐增大，这会一定程度上影响到迭代格式的收敛特性. 本文选用这一特别的来流马赫数对发展的算法进行了收敛历程测试. 计算采用了75060个无网格点(见图17). 图18则给出了对应的收敛历程，图中可以看出，加预处理可以明显改善低马赫数来流的收敛特性，一定程度上有助于提高解的捕捉精度(参见表面压强系数等值线云图19). 图18中同时给出了来流马赫数0.2的收敛历程供比较，也可以看出，随着来流马赫数的减小，流场中低马赫数的不可压区占比扩大，如预期，加预处理的必要性更加显现，表现为加与不加预处理，收敛差异更大.

图17 DLR-F4 翼身组合体周围点分布Fig.17 Points around DLR-F4

图18 收敛历程Fig.18 Convergence histories

图19 组合体表面压强系数等值线云图Fig.19 Pressure coefficient contours

5 结 论

本文通过无网格点重排与点云分割，成功地结合LU-SGS算法，发展出一种基于Weiss-Smith预处理矩阵的预处理隐式无网格算法，用于求解三维欧拉方程. 借助于无网格布点，远场边界虚拟镜像点的引入，简化了预处理远场边界的实施. 发展的算法通过了二维和三维多个算例的计算考核. 算例表明，相较于显式预处理算法，发展的预处理隐式算法收敛特性更优；算法的数值模拟能力已从单纯模拟可压缩跨音速等流动拓展到模拟几乎不可压的低马赫数流. 此外，算法整体是基于无网格技术发展的，计算域只需布点离散，有助于灵活处理三维复杂气动外形. 当然，实际工程绕流问题大多需考虑粘性的影响，这就要求发展的算法从求解欧拉方程拓展到求解Navier-Stokes方程，会涉及反映粘性绕流特征的异性点云结构的构造、合适湍流模型的嵌入、基于点云结构的粘性预处理谱半径的确定等问题，将是未来研究的重点.