Walsh,函数有限体积法的多重网格特征研究1)

王刚 干源 任炯

(西北工业大学航空学院,西安 710072)

为了满足不断增长的工程应用需求,计算流体力学工作者一直在努力发展各种新型的数值算法,提升算法和软件的精准度、计算效率和鲁棒性,追求以更少的计算代价高效地获得更加准确和精细的流场模拟结果[1-2].为此,人们发展了各种高分辨率、高精度的有限差分方法[3]、有限体积方法[4]和有限元方法[5-6]来提高流场模拟的精准度,其中包括近些年衍生出的重构间断有限元方法(reconstructed discontinuous Galerkin(RDG) method)[7-8]、谱体积方法(spectral volume(SV) method)[9-10]和谱差分方法(spectral difference(SD) method)[11-12]等.目前的高精度方法依然存在两个方面的缺憾:其一,在网格规模固定的情况下,通过提高格式精度提升流场模拟精细程度的同时,会带来求解变量数目(自由度)的大幅增加,进而导致计算量剧增、收敛速度下降和鲁棒性变差等问题;其二,这些方法大都假设流动变量在网格内部连续变化,仅允许流场间断出现在网格单元的交界面处[13].即使增加了格式精度,也不能从根本上改善间断捕捉的分辨率.

为了破除流动参数在网格内部连续变化这个预设条件对分辨率的限制.任炯等[14-15]提出了一种可以在网格内部捕捉间断的Walsh 函数有限体积方法(finite volume method with Walsh basis functions,FVM-WBF).FVM-WBF 方法利用具备间断特征的Walsh 正交基函数[16]将计算网格单元内部的解变量表达为级数展开形式的离散逼近[17-18],从而在网格单元内部刻画出流动变量的分片连续区域,针对每个分片连续区域进行流动控制方程的有限体积离散,最终建立可以在网格单元内部捕捉间断的高分辨率算法.这种方法虽然能够实现网格内部的间断捕捉,但其流动模拟分辨率的高低依然取决于Walsh 基函数级数的规模,亦即参与离散的Walsh 基函数的数目.与DG和SV 等高精度方法类似,FVM-WBF 方法在增加基函数数目提高求解分辨率的同时,其计算量增加和收敛速度下降的问题同样会不可避免地显现[14].在提高FVM-WBF 方法流场模拟精细程度的同时,发展合适的加速收敛算法,提升流场求解的效率显得非常必要.

当前,比较实用的流场加速收敛算法主要是当地时间步长[19]和多重网格方法[20-21].当地时间步长方法的主要原理是:在定常流动问题的求解过程中,在每个计算网格上采用满足稳定性约束的最大时间步长来进行时间推进,从而提高流场求解的收敛速度.多重网格方法根据实现策略的不同划分为几何多重网格(geometrical multigrid,GMG)[22]和代数多重网格(algebraic multigrid,AMG)[23]两类.GMG 方法通过对细网格的连续“粗化”形成单元尺度由细到粗的多重网格系列,引入满足守恒性条件的插值、限制算子传递不同粗细网格上的求解信息[13],将细网格上难以消除的低频误差以粗网格上高频误差的形式快速消除.GMG 方法实现的难点在于网格的粗化过程和插值/限制算子的合理构造[24-25],特别是对于非结构网格和大规模分布式并行计算情景.GMG方法中比较常见的是聚合多重网格[26-27]的策略,即通过网格的不断聚合形成一系列不同几何尺度的粗网格,这会增加网格前处理的难度,在复杂的几何拓扑下,它们还可能存在鲁棒性问题.AMG 主要用于隐式求解控制方程,其不依赖于空间几何信息,也不需要在粗网格上计算通量和源项[28],仅通过“粗化”系数矩阵即可完成粗细网格之间的信息传递[29].AMG 方法凭借上述优点,得到了广泛的关注[30],但在非线性效应显著的求解问题中,其加速收敛效果不及GMG 方法[28].

在前文所提及的FVM-WBF 方法中,不同级别的基函数级数项所影响的局部均值区域尺度是有“粗细”之分的,在迭代求解流场的过程中,不同级别基函数的系数之间相互影响,天然隐含着类似几何多重网格的不同空间尺度上的流场信息传递,具备自动引入多重网格策略的潜力.本文将依据FVMWBF 方法内在的多重网格特征,提出一种结合多重网格策略的新型FVM-WBF 方法,以期在流场计算中获得显著的加速收敛效果.

本文将首先对FVM-WBF 方法的数值原理进行简要介绍,接着通过理论分析,推导FVM-WBF 方法中实现多重网格加速收敛的主要原理和算法,进而提出一种结合多重网格策略的FVM-WBF 方法.最后以圆柱和NACA0012 翼型两种外形的定常无黏绕流为数值算例,对上述改进型FVM-WBF 方法的多重网格特征进行算例验证,通过与不引入多重网格策略的情形进行对比,考察所发展的结合多重网格策略的FVM-WBF 方法的计算效率.

1.1 二维Walsh 基函数级数逼近

不同于传统的有限体积方法,FVM-WBF 方法利用具有间断性质的Walsh 基函数,对网格单元内部的流动守恒变量Q(x,y,t) 进行级数逼近

其下标(m,n) 对应两组一维Walsh 基函数的张量积顺序.cm,n(t) 是二维Walsh 基函数gm,n(x,y) 对应的基系数.一维Walsh 基函数gm(x)和gn(x) 可由简单方波经过一系列二分演化形成,具体过程可参考文献[14].延拓一维Walsh 基函数的层级划分方式[15],二维Walsh 基函数gm,n(x,y) 所属的层级p可由下式确定

其中,p1和p2分别表示两个一维Walsh 基函数gm(x)和gn(y) 的所属级别.

按照式(3),定义在单位正方形区域上的前三级二维Walsh 基函数的取值如图1 所示.图中蓝、红色部分分别对应基函数正、负值区域.可以看出二维Walsh 基函数具有与一维Walsh 基函数类似的性质:(1)具有等幅方波的属性;(2)增加第p级基函数,网格单元内原有的分片连续子区域在x和y方向上分别递进二分,形成新的 1 /2p尺度的子区域;(3)p≥1时,第p(p≥1) 级的基函数在前p级基函数形成的各子区域上的面积分始终为零.

图1 单位正方形上二维Walsh 基函数示意图Fig.1 Schematic diagram of the two-dimensional Walsh basis function on the unit square

根据上述性质,Walsh 基函数级数式(1)表达的是一个尺度由大到小的逐级逼近过程,第p级基函数对应着 1 /2p尺度的流动变量描述,且基函数的级数越高,对流动变量的描述尺度就越小.在前p+1级Walsh 基函数级数式的基础上继续延伸一级基函数,新增加的第p+1 级基函数的系数只决定新产生的12p+1网格尺度上的流动变量描述,既有尺度上的流动变量描述不受影响.

为了直观理解上述Walsh 基函数级数数值逼近特性,图2 以单位正方形区域上的解析函数Q(x,y)=-(x-0.4)2-(y-0.6)2+0.1xsin(2y)+1为例,具体演示了Walsh 基函数级数式的逐级数值逼近过程.图2 中左上方给出了级数式数值逼近时的前三级基系数,三维立体图中的蓝、红色分别表示基函数项取值的正负,1×1,2×2和4×4 分别代表采用第一级、前两级、前三级Walsh 基函数时形成的分片连续子区域数目.

图2 FVM-WBF 方法对单位正方形区域上的变量Q(x,y)=-(x-0.4)2-(y-0.6)2+0.1xsin(2y)+1拟合数值结果Fig.2 The numerical results of two-dimensional variable Q(x,y)=-(x-0.4)2-(y-0.6)2+0.1xsin(2y)+1simulated by the FVM-WBF method on unit square area

从图中可以看到,当采用第一级Walsh 基函数级数式对变量进行逼近时,整个控制体内变量分布由基函数g1,1和对应的系数c1,1确定,表示为常数分布.在此基础上,累加第二级基函数组合项c1,2g1,2+c2,1g2,1+c2,2g2,2,并不会改变第一级基函数系数c1,1,仅涉及每个1/2 网格尺度区域上的均值调整,实际效果是将控制体内的变量描述成了2×2 个分片连续的分布.继续引入第三级Walsh 基函数,新增加的12 个基函数组合项并不会改变前两级的基函数系数,仅对每个1/4 网格尺度区域上的均值进行调整,效果是将控制体内的变量描述成了4×4 个分片连续的分布.以此类推,在前p级基函数级数式的基础上,增加第p级基函数组合项,其效果只是调整每个 1 /2p网格尺度区域上的均值,对前p级基函数对应的系数不产生影响.显然,当p增加到足够大时,控制体内变量分布趋近解析分布.

以上主要是针对正方形网格的描述,借助双线性映射[31],可以将式(2) 所定义的Walsh 基函数gm,n(x,y)由单位正方形网格转换到二维任意凸四边形网格,并构建出任意四边形网格上的Walsh 基函数级数式.图3 以前三级Walsh 基函数级数式为例,给出了其在任意四边形网格上形成的4×4 个分片连续离散子区域示意,其中子区域 Ωi,j的下标(i,j) 由两组对应一维Walsh 基函数级数形成的子区间下标i和j按照张量积顺序给定.

图3 四边形网格的Walsh 基函数级数子区域划分示意图Fig.3 Subdivision of a quadrilateral grid according to Walsh basis function series

1.2 二维Euler 方程的FVM-WBF 离散

下面对求解二维Euler 方程的FVM-WBF 方法进行简要介绍.考虑二维守恒型Euler 方程,其积分形式如下

对于单个守恒变量,上式中共包含M×N个待求解的系数.显然,仅靠方程(5)不能提供足够的定解条件.由于Walsh 基函数本身具有间断特征,在控制体内不能保证处处可微,因此无法利用类似有限元的变分过程[6]建立定解条件.为了获得满足定解条件的方程数目,FVM-WBF 方法将方程(5)在Walsh 基函数级数式形成的各子区域上进行有限体积的空间离散,得到如下关于基函数系数时间导数的半离散方程

由于矩阵A可逆[14],因此可将方程组(7)改写为

式中,残值项R=A-1F.以Walsh 基函数系数c(t)为解变量的方程组(8)可以采用经典(传统)有限体积方法中各类成熟的时间离散格式进行推进求解.本文使用的是满足TVD 性质的显式三步龙格-库塔格式[32],具体步骤如下

式中上标n和n+1 代表时间层级.cn和cn+1分别为当前时间层和下一时间层的基函数系数,残值项R(k)由第k步基函数系数c(k)所确定的流动变量经过运算得到,是时间步长.在定常问题的求解中,传统有限体积方法一般采用当地时间步长方法进行计算,具体计算方法如下

这里 σ为CFL 数,| Ω|和c分别为网格单元的面积和当地声速,ΔS1和ΔS3,ΔS2和ΔS4分别是网格单元Ω两两相对的边界边长,n1,n2,n3和n4分别是对应边界上的外法线矢量.在早前的FVM-WBF 方法中[14],考虑到第p级Walsh 基函数系数的时间推进应满足1/2p网格尺度区域上的稳定性条件约束,如果对Ω上的全部待解系数使用相同的时间步长可将其近似取为Δt/2p.

1.3 采用多重网格策略的FVM-WBF 方法

Walsh 基函数级数的数值特性决定了方程组(8)中各级基函数系数cm,n在网格单元内影响的局部均值区域尺度不同.方程组(8)中的方程根据影响的均值区域尺度划分为一系列连续的级别,在时间推进过程中隐含着不同尺度的流场信息传播,形成了多尺度均值耦合求解的封闭系统,该系统提供了天然的多重网格平台.显然,对于不同级别区域尺度的均值,其相应的方程求解对应的当地稳定性条件约束也不一样.早前的FVM-WBF 方法按照控制体Ω中最小的 1 /2p网格尺度区域所约束的稳定性条件,将方程组(8)中各个待解系数的时间步长统一设置为Δt/2p.这种处理方式过于保守,会导致大于1/2p网格尺度区域的流场信息以偏小的时间步长传播和演化,导致低频误差的消除速度减慢,给计算效率带来不利影响.

为了消除上述缺陷,本文提出一种采用多重网格策略的FVM-WBF 方法.在定常问题求解中,放弃对方程组(8)采用统一的时间步长Δt进行时间推进,而是根据不同级别系数方程所关联的网格尺度,依据对应的当地稳定性约束条件,对每个方程取尽可能大的时间步长进行迭代,在各级网格尺度上以较快的速度同时消除不同频率的误差,从而达到提高整体流场收敛速度的目的.图4 具体给出了本文对不同基函数系数方程的时间步长取定方案.第0级方程是对Walsh 基函数系数c1,1的迭代求解,对应的是整体网格单元尺度的流场信息传播,因此其时间步长应该依据整体网格的几何尺度获得,即由式(10)计算得到的Δt;第1 级方程是对Walsh 基函数系数c1,2,c2,1和c2,2的迭代求解,对应的是1/2 网格单元尺度的流场信息传播,其时间步长应该基于1/2 网格尺度获得,即Δt/2 ;依次类推,第p级方程的时间步长依据 1 /2p网格尺度获得,即Δt/2p.

根据以上分析和图4 中的方案,考虑将方程(9)中的时?间步长改写为如下块对角矩阵

图4 不同级别的基函数所对应的时间步长Fig.4 Correspondence between different levels of basis functions and time steps

其中,Δtp是1 /2p网格尺度级别的方程对应时间步长形成的对角方阵子块,其阶数L由第p级Walsh 基函数的数目决定.对于采用多重网格策略的FVMWBF 方法,矩阵子块Δtp的取值如下

式中Δt由式(10)计算获得.不同层级的基函数系数所能调整的均值区域尺度对应着不同波长的流动信息传播,对于不同层级的方程采用式(11)设置不同的时间步长进行推进求解,这与几何多重网格方法中不同粗细的网格上采用不同时间步长的原理基本类似,其效果都是快速地消除不同波长(频率)的误差,进而提高流场计算的收敛速度.

为了方便描述,采用符号“MG”将采用了多网格策略的FVM-WBF 方法简写为FVM-WBF_MG 方法.由于FVM-WBF 方法提供了天然的多重网格平台,因此新的多重网格策略FVM-WBF 方法不存在传统有限体积方法实施几何多重网格算法时所需的插值和限制算子,相较于早前的FVM-WBF 方法,其求解的控制方程的数目不变,每一步的迭代过程并不会额外增加计算时间.

选取圆柱和NACA0012 翼型定常无黏绕流为算例测试FVM-WBF_MG 方法的加速收敛效果.在本文算例中,为了对FVM-WBF 方法中使用前两级、前三级和前四级Walsh 基函数的情形进行区别,采用每种情形所形成的分片连续子区域数目“2×2”,“4×4”和“8×8”对其进行标识.所有计算的空间离散均采用二阶精度的熵相容格式[33-34].对能量方程的最大残差

收敛历程和力系数渐进误差

的收敛历程进行考察,衡量和对照FVM-WBF 方法使用多重网格策略前后的收敛速度.其中,式(12)中的ΔEi,j为网格单元内子区域Ωi,j上的能量方程残差,式(13)中的C和Cconv分别代表当前迭代步和充分收敛后获得的升力或阻力系数值.

2.1 圆柱绕流

圆柱不可压无黏绕流作为存在解析解的经典的流体力学问题之一,经常被用以检验数值方法的合理性和有效性.本算例选取马赫数Ma=0.1 的计算状态,对FVM-WBF 方法的多重网格性能进行测试.时间推进的CFL 数取为0.5.图5 给出的是计算网格分布情况,网格数目为80×60.

图5 二维圆柱的计算网格Fig.5 The computational grid of the two-dimensional cylinder

图6 对比了采用与未采用多重网格策略的FVM-WBF 方法所获得的流场马赫数云图,可以看到两者完全相同.图7 给出了使用不同Walsh 基函数级数的情况下,采用与不采用多重网格策略的FVM-WBF 方法计算得到的壁面压力系数分布.与势流理论的解析解进行对比,不难发现,随着基函数级数的延伸,压力系数逐渐逼近解析解.图6和图7的结果显示:在基函数级数相同的情况下,使用本文提出的多重网格策略并不会影响流场最终的收敛结果,从实践角度验证了该策略的合理性.图8 给出了采用与不采用多重网格策略的FVM-WBF 方法所获得的能量方程残差和阻力系数渐进误差关于迭代次数(iteration)的收敛历程对比.可以看到,使用多重网格策略后的残差和力系数收敛所需迭代次数大幅减少.根据1.3 节的推导,本文多重网格策略不会增加FVM-WBF 方法在每个迭代步的计算时间,因此迭代步数的减少可直接反映CPU 时间的降低.对比图8 中相同基函数级数的计算结果,可以发现,多重网格策略具有显著的加速收敛效果,验证了多重网格策略的有效性.进一步对比不同基函数级数的结果,可以发现,不同基函数级数使用多重网格策略后的加速程度不同,随着Walsh 基函数级数式的延伸,其加速效果越明显.

图6 圆柱绕流马赫数云图Fig.6 Mach number contours of the low speed flow over cylinder

图7 圆柱表面压力系数对比图Fig.7 Comparison of pressure coefficients on cylinder surface

图8 FVM-WBF_MG 方法在圆柱低速绕流计算中的加速效果Fig.8 Effect of FVM-WBF_MG method in the calculation of the low speed flow over cylinder

2.2 NACA0012 翼型绕流

NACA0012 翼型作为航空工程应用中的典型模型,其绕流问题经常被用以检验数值方法的合理性和有效性.选取NACA0012 翼型低速、跨声速和超声速三个流动状态,进一步对FVM-WBF 方法的多重网格加速收敛性能进行测试.三个状态的计算中,时间推进的CFL 数均取0.5.图9 给出的是计算网格分布情况.网格数目为60×20.下面分别对三个计算状态的数值结果进行描述和分析.

图9 NACA0012 翼型计算网格Fig.9 The computational grid of the NACA0012 airfoil

(1)低速绕流

NACA0012 翼型低速绕流的计算状态为:马赫数Ma=0.15,迎角 α=10°.

图10给出了采用与未采用多重网格策略的FVM-WBF 方法所获得的流场马赫数云图对比,可以看到两者完全相同,说明了本文多重网格策略在NACA0012 翼型低速绕流计算中的合理性.图11 给出了使用不同Walsh 基函数级数的FVM-WBF 方法采用与未采用多重网格策略获得的能量方程残差和升力系数渐进误差关于迭代次数的收敛历程对比.可以发现,收敛到相同的误差水平时,使用多重网格策略的残差和力系数渐进误差所需迭代步数均大幅降低.对比相同和不同基函数级数的结果,均具有与上述圆柱低速绕流算例类似的表现,表明了多重网格策略在NACA0012 翼型低速计算中的有效性.

图10 NACA0012 翼型低速绕流马赫数云图Fig.10 Mach number contours of the low speed flow over NACA0012 airfoil

图11 FVM-WBF_MG 方法在NACA0012 翼型低速绕流计算中的加速效果Fig.11 Effect of FVM-WBF_MG method in the calculation of the low speed flow over NACA0012 airfoil

(2)跨声速绕流

NACA0012 翼型跨声速绕流的计算状态为:马赫数Ma=0.8,迎角 α=1.25°.

图12 给出了采用与未采用多重网格策略的FVM-WBF 方法所获得的流场密度云图,其中基函数级数选取前两级和前四级.可以发现:多重网格策略并不影响最终收敛后的流场结果;在FVM-WBF方法中,使用的基函数Walsh 级数越高,流场中激波结构的捕捉就更为精细.图13 以3 阶精度DG 方法[35]的计算结果作为参考值,对比了使用不同Walsh 基函数级数的FVM-WBF 方法采用与不采用多重网格策略计算得到的压力系数.流场和压力系数结果表明,基函数级数的级别越高,捕捉到的激波越锐利,其压力系数值越逼近3 阶精度参考值,然而相同基函数级数的FVM-WBF 方法采用与不采用多重网格策略的压力系数分布完全相同,说明本文多重网格策略并不会影响NACA0012 翼型跨声速流场的最终收敛结果,再次验证了该策略的合理性.图14 分别给出的是跨声速状态下的三种不同级数的FVM-WBF 方法采用与未采用多重网格策略所获得的能量方程残差和升力系数渐进误差关于迭代次数的收敛历程对比.不难发现,在NACA0012 翼型的跨声速计算中,多重网格策略能够以很少的迭代步数收敛到与不使用多重网格策略的计算结果相同的误差水平,说明本文多重网格策略在跨声速状态下也能够显著提升收敛速度.

图12 NACA0012 翼型跨声速绕流密度云图Fig.12 Density contours of the transonic flow over NACA0012 airfoil

图13 NACA0012 翼型表面压力系数对比图Fig.13 Comparison of pressure coefficients on NACA0012 airfoil surface

图14 FVM-WBF_MG 方法在NACA0012 翼型跨声速绕流计算中的加速效果Fig.14 Effect of FVM-WBF_MG method in the calculation of the transonic flow over NACA0012 airfoil

(3)超声速绕流

NACA0012 翼型超声速绕流的计算状态为:马赫数Ma=2.0,迎角 α=3°.

图15 中对比了采用与不采用多重网格策略的FVM-WBF 方法所获得的流场密度云图,其中基函数级数选取前两级和前四级.上述结果再次证明:采用多重网格策略不会改变收敛后的流场解;增加Walsh 基函数的级数可以有效提高流场求解的分辨率,在较稀疏的网格上精确地捕捉激波.图16 给出了使用与未使用多重网格策略的FVM-WBF 方法所获得的能量方程残差和升力系数渐进误差关于迭代次数的收敛历程对比.图中结果显示,本文多重网格策略在NACA0012 翼型的超声速计算中再次表现出良好的加速收敛效果.

图15 NACA0012 翼型超声速绕流密度云图Fig.15 Density contours of the supersonic flow over NACA0012 airfoil

图16 FVM-WBF_MG 方法在NACA0012 翼型超声速绕流计算中的加速效果Fig.16 Effect of FVM-WBF_MG in the calculation of the supersonic flow over NACA0012 airfoil

综合以上三种计算状态的数值结果比较,表明前述1.3 节中提出的多重网格策略在各种状态的定常流场数值模拟中均表现出显著的加速收敛效果.这与传统几何多重网格方法通过粗化增加网格层数进一步提高收敛速度的原理类似,Walsh 基函数级数式中包含的基函数级别越多,FVM-WBF 方法求解系统中隐含的不同级别的网格尺度也越多,相当于多重网格方法使用的网格重数越多,不同频率误差的消除速度均能得到提升,加速效果就越明显.

本文通过理论推导揭示了FVM-WBF 方法天然具备多重网格的算法特征.利用这一特征发展了一种结合多重网格策略的FVM-WBF 方法.通过二维圆柱和低速、跨声速和超声速下NACA0012 翼型绕流数值模拟,对该方法求解定常流场的加速收敛效果进行了测试和验证.上述工作可以得出以下主要结论:

(1)新发展的多重网格策略实现过程简单,不必进行网格的粗化聚合和构建插值、限制算子进行信息传递,在不引入额外计算量的同时,保证参数传递的守恒性,在同类多重网格方法中具有独特的优势.

(2)结合多重网格策略的FVM-WBF 方法不仅拥有良好的加速收敛效果,其稳定性也可得到很好的保证,能够适应低速、跨声速和超声速不同速度范围的流场数值模拟.

本文工作主要是以求解Euler 方程验证算法理论为主要目的.FVM-WBF 方法的多重网格特征在更复杂的黏性和三维流动问题中的应用问题将在后续研究工作中进一步开展.

猜你喜欢 级数步长流场 车门关闭过程的流场分析汽车实用技术(2022年19期)2022-10-19中心差商公式变步长算法的计算终止条件肇庆学院学报(2022年5期)2022-09-29基于Armijo搜索步长的BFGS与DFP拟牛顿法的比较研究成都信息工程大学学报(2021年5期)2021-12-30求收敛的数项级数“和”的若干典型方法数学学习与研究(2020年11期)2020-09-11基于随机森林回归的智能手机用步长估计模型中国惯性技术学报(2020年2期)2020-07-24无穷级数的柯西和与切萨罗和数码世界(2019年8期)2019-08-15一个非终止7F6-级数求和公式的q-模拟华东师范大学学报（自然科学版）(2019年3期)2019-06-24基于CFD新型喷射泵内流场数值分析舰船科学技术(2016年1期)2016-02-27基于动态步长的无人机三维实时航迹规划北京航空航天大学学报(2016年12期)2016-02-27天窗开启状态流场分析汽车实用技术(2015年8期)2015-12-26

相关热词搜索：网格 体积 函数