非均匀矩形网格的局部网格细化LBM 算法研究1)

2023-11-16孟欣雨杨双骏桑为民

力学学报 2023年10期

关键词：物面插值计算结果

安博孟欣雨杨双骏桑为民,2)

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

† (中国空气动力研究与发展中心结冰与防除冰重点实验室,四川绵阳 621000)

** (翼型、叶栅空气动力学重点实验室,西安 710072)

引言

格子玻尔兹曼方法[1-3](LBM)在经历几十年的快速发展后已逐渐成为流体力学研究中的重要分支和热点话题.同时在计算流体力学(CFD)领域中得到广泛应用,比如不可压缩流动[4]、可压缩流动[5]、多孔介质流[6]、气动声学[7]、多相流[8]、燃烧[9]、渗流[10]等.由于其经典的碰撞迁移理论完美地契合均匀正方形网格的网格结构,因此在LBM 相关数值模拟中应用非常广泛.虽然基于均匀正方形网格的单松弛LBM (传统LBM)在数值模拟方面表现不俗,但由于其鲁棒性和数值稳定性较差,尤其当松弛因子 τ=0.5+3UL/(Re∆x)接近0.5 时,数值稳定性逐渐丧失[1-3],因此其应用范围基本仅限于低雷诺数Re的数值模拟.为了应对逐渐增大的雷诺数,就必须减小网格间距 ∆x,同时还能满足物面处以及流动变化剧烈区域流动细节的捕捉.而均匀正方形网格的细化必然导致网格总量的骤然增加,使得计算效率大幅下降,且对计算设备的要求也急剧提高.极大地限制传统LBM 的应用与发展.

作者在之前的研究中提出了基于树网格的LBM算法[11-12],可以在保证计算精度,满足局部网格细化要求的前提下,大幅提高传统LBM 的计算效率.此外,众多学者也开展了传统LBM 局部网格细化的算法优化研究[13-27].比如有限差分LBM[13-14]、插值补充LBM[15]、多重网格LBM[16-18]、多区块LBM[19-21]、树网格LBM[22-25]以及均匀矩形网格LBM[26-27]等.Mei 等[13]以及Guo 等[14]借鉴传统的求解N-S 方程的手段,将有限差分方法与LBM结合.由于控制方程,即格子玻尔兹曼方程不再通过碰撞迁移演化求解,转而依靠差分格式逼近.所以网格结构的使用也变得较为灵活,随着贴体网格的引入,在网格整体数量不变的情况下解决了网格物面处的局部加密的问题.但是贴体网格的使用会带来复杂的坐标变换以及在求解过程中需要构造高精度差分格式,总体来讲未能有效提升计算效率.同样地,基于贴体网格,He 等[15]使用的插值补充格式,依然保留了LBM的精髓,即分布函数的碰撞迁移过程,然后通过插值来近似求解碰撞迁移后虚拟点的分布函数,但计算过程比较复杂且插值格式精度较低.Marvriplis[16]、Patil 等[17]和王兴勇等[18]将传统CFD 中的多重网格技术应用于LBM 的数值模拟中.其核心思想是通过粗细网格的信息交换,使得低频误差在粗网格快速降低,高频误差在细网格得到光滑处理,从而提升计算效率.总体来讲,多区块LBM[19-21]和树网格LBM[22-25]算法基本是殊途同归,虽然每个学者的实现手段不尽相同,但是核心思想是一致的.即对流场内部的不同区域采取局部加密.这些区域往往是流动变化较为剧烈的区域或者是物面流动区.Zhou[26-27]提出了基于均匀矩形网格的LBM 算法,不同于本文非均匀矩形网格LBM 算法,Zhou 的计算模型没有插值过程,用重新构建的矩形离散速度模型以及对应的平衡态分布函数取代传统的正方形离散速度模型(见图1),可以保证分布函数以矩形方式完成格点之间的碰撞迁移.虽然该方法可以在一定程度上减少总体网格数量,但局部网格细化效果不佳且不具备普遍的适用性.以上国内外的相关研究对本文工作予以很大的启发,尤其是插值补充LBM 算法[15]中虚拟点分布函数的插值求解过程,对本文基于非均匀矩形网格25 点拉格朗日插值LBM 算法的构建提供了思路.

图1 LBGK-D2Q9 离散速度模型Fig.1 Discrete velocity model of LBGK-D2Q9

1 25 点插值LBM 计算模型

本文研究使用最常用的LBGK-D2Q9 模型[28].其平衡态分布函数的构造如下式所示

其中,ωα是离散速度模型中不同离散方向的权系数;ρ为流体密度;u为流体速度矢量;eα是离散速度,其大小和方向构造如下所示

其中,c=∆x/∆t=1 为格子速度;cs为格子声速;∆x为网格步长;∆t为时间步长.图1 介绍了LBGK-D2Q9计算模型的离散速度模型.详细说明了分布函数fα的碰撞迁移方式.

式(3)定义了LBGK 计算模型中碰撞迁移理论的控制方程,其中r为空间向量;ξ 为速度向量;τ 为松弛时间;分别为分布函数和平衡态的分布函数.

下式定义了流场宏观物理量,速度u和密度 ρ 与微观粒子分布函数fα之间的联系.每一次演化过后,流场的速度和密度都会通过下式得到更新

传统的LBM 算法基于均匀正方形网格,根据其碰撞迁移理论,分布函数的迁移过程必须保证在网格格点上.因为在传统的LBM 算法中,格子速度 ψ.如图2 所示,点A和O为网格格点.点O处的流体粒子,其2 方向的分布函数经过一个时间步长 ψ 的演化会迁移至A点,迁移距离为网格间距 ω,即.其他方向的分布函数也遵循同样的规则迁移至O点周围的各点.由此可见均匀正方形网格天然匹配LBM 的演化机理.通常情况下,为了得到较为理想的数值模拟结果,尤其是捕捉物面处以及流动剧烈变化区的流动细节,往往需要极为细腻的网格.这就使得均匀正方形网格的网格总量急剧增加,极大地降低了传统LBM 的计算效率和鲁棒性,限制了LBM 的应用.因此,本文采用了非均匀矩形网格用于流场局部加密既,保证了物面处和流动变化剧烈区的流场细节捕捉,同时相较于均匀正方形网格,在保证计算精度的前提下又能提升LBM 的计算效率.

图2 均匀直角网格的迁移演示Fig.2 Streaming on uniform Cartesian grid

由于本文的计算网格是非均匀的矩形网格,每一层网格的网格间距都不相同,传统LBM 中基于格点和格点之间迁移条件无法满足.

如图3 所示,点A和O均为网格格点.由于点A所在的网格层和点O所在的网格层之间的网格间距不同,∆xOA<∆xOC.所以当经过一个时间步长∆t的演化后,针对O点,只有2 方向的分布函数f2可以正常迁移到同为网格格点的A点其他方向的分布函数均无法完成正常迁移,比如,7 方向的分布函数只能迁移至B'点,而此点为虚拟点,且我们并不知道此点的信息.因此本文数值模拟,除了LBM 原有的碰撞迁移计算还需要额外构造一个插值格式近似求解B'点的信息.在插值LBM 算法[15]的基础上,本文提出了基于25 点拉格朗日插值格式的非均匀矩形网格LBM 算法.

图3 非均匀矩形网格的迁移演示Fig.3 Streaming on non-uniform rectangular grid

图4,以A33点5 方向的分布函数f5为例,介绍了本文25 点拉格朗日插值格式的构造过程.由图可见点Aij(i=1,2,···,5;j=1,2,···,5)均为计算域的实际网格点,而点aij(i=1,2,···,5;j=1,2,···,5)为沿5 方向迁移后的虚拟点,并且包括点A33以及其周围的24 个点都参与了插值格式的构建.网格沿横(ξ1,ξ2,ξ3和 ξ4)、纵(η1,η2,η3和 η4)向均有4 层网格且网格间距均不相同.∆=∆t·c为流体粒子单位时间步长 ∆t内移动的距离(c=1.0),同时也是最细网格间距(∆=∆xfinest).本文沿横、纵向靠近物面的第一层最细网格为 ∆xfinest=1/1024,见图5(b)中计算域四角黑色区域.以5 方向的分布函数为例,下式给出了本文25 点拉格朗日插值格式

图4 非均匀矩形网格的迁移演示(α=5)Fig.4 Streaming on discrete direction (α=5)

图5 网格对比Fig.5 Comparison between different mesh

2 数值模拟

2.1 计算域构建及边界条件

为了验证本文基于25 点拉格朗日插值的非均匀矩形网格LBM 算法,选取经典的顶盖驱动方腔内流为基准算例(Ma=0.173 2 且Re<10 000).本文数值模拟最细网格步长(∆xfinest=1/1024)的选取借鉴了作者之前研究工作[29-32]中针对腔体内流的网格独立性验证.基于此网格尺度,发现当Re<10 000 时能够保证y+<0.1,数值模拟结果可信且能够捕捉边界层附近的流场信息.如图5 所示,展示了顶盖驱动方腔内流的网格结构,图5(a)和图5(b)分别对应均匀正方形网格和非均匀矩形网格.流动信息采集点位于腔体中心P(x=L/2,y=L/2)处.

为了便于研究,本文非均匀矩形网格的建立借鉴了等比数列的控制函数,将计算域沿横向和纵向划分为4 个对称的子区间.每个子区间网格间距由下式控制

其中,N为网格点个数,σ 是等比控制函数的公比,在本文中当N=500 时,σ 为0.005 213 64.∆xfinest=∆yfinest定义了贴近物面处第一层网格的网格间距(最细网格).

2.2 边界条件处理

本文涉及的边界均为平直边界(方腔4 个边).其中顶盖为驱动边界,其余三边均为物面边界.由于非平衡态外推格式[33]易于代码编写且非常适用于平直边界点的整体处理.本文延续之前研究中已成熟应用的非平衡态外推格式,将边界点上的分布函数分为平衡态和非平衡态两部分.其中平衡态部分由平衡态分布函数的定义近似获得,而非平衡态部分则用非平衡态外推法求解.

如图6 所示点A,B和C为物面边界点,D点为流场点.根据LBM 的演化原理可知在每次演化之前需要知道每个点的分布函数,对于B点,其分布函数可分为平衡态和非平衡态两部分

图6 平直物面边界Fig.6 Straight wall boundary

以物面边界点B为例,该点的平衡态分布函数可用该点的宏观物理量来构造,如果该点的宏观物理量未知,则由D点的相应值代替.而非平衡态分布函数则由D点的非平衡态分布函数来近似代替

3 计算结果与分析讨论

3.1 定常流动数值模拟

本文以顶盖驱动方腔内流为例,针对本文RLBML25 算法在定常(Re=1000 和Re=5000)和非定常状态(Re=8800 和Re=9500)下开展数值模拟分析验证研究.

表1 介绍了当Re=1000 时本文的计算结果和其他文献[34-35]的数据对比.第1 列P 代表流场中涡漩的位置,LB 表示左下侧次级涡;RB 表示右下侧次级涡;PC 代表中心主涡(见图7).对应每个涡结构给出了4 个流动参数,分别为涡心横坐标(X)、纵坐标(Y)、流函数(ψ)和涡量(ω).由表1 可知,本文提出的RLBM-L25 算法其结果与其他文献的数据吻合较好,可以输出可信的数值模拟结果.

表1 不同计算结果对比 (Re=1000)Table 1 Different results comparison (Re=1000)

类似表1,表2 列出了当Re=5000 时的数据对比.数据结构和内容与表1 一致.再次证明RLBML25 算法对于定常流动数值模拟表现良好.

表2 不同计算结果对比 (Re=5000)Table 2 Different results comparison (Re=5000)

如图7(a)和图7(b)所示,展示了定常状态下的流场拓扑结构.可以明显观察到经典方腔驱动内流的典型流动特征,与其他文献[34-39]的数据完全一致.中心主涡覆盖大部分计算域从而主导整个流场,随着雷诺数的增加,次级涡逐渐增大并伴随边角处出现细碎涡的现象.图7(c)展示了Re=1000 时流场的局部放大图,左侧是基于均匀正方形网格的传统LBM 算法(SLBM)的计算结果,右侧是本文RLBML25 的计算结果.通过对比可以看出,本文基于非均匀矩形网格的网格加密方法可以保证物面处和流动突变区的流动细节捕捉.图7(d)和图7(e)分别展示了Re=1000 时基于均匀正方形网格的经典单松弛LBM(SLBM)且网格分辨率为1024 和基于非均匀矩形网格插值LBM(RLBM-500)且网格分辨率为500 时的涡量图对比,如图所示,涡量图吻合得较好.

图8 描述了沿顶盖分布的剪切力fs,总体来讲,本文RLBM-L25 算法网格分辨率分别为500 和600 且雷诺数等于1000 时的计算结果与传统LBM算法网格分辨率为1024 时的结果较为吻合,但是在顶盖左右两侧附近仍存在一定量级的误差且不会因为插值格式的改变而消除,因此作者认为这种情况的出现首先跟顶盖上网格分布数量有关,其次跟网格局部细化的控制方程也有关.

图8 顶盖剪切力分布Fig.8 Shear force fs distribution along the lid

图9 速度均方根误差曲线Fig.9 History of root mean square error eRMS of velocityux

图10 扰动衰减率曲线Fig.10 History of perturbation decay rate ε of velocity ux

为了进一步讨论分析本文算法的表现,还开展了基于不同网格分辨率的研究.以传统算法SLBM且网格分辨率为1024 时的计算结果作为对比.

由表3 可见,第1 列给出了不同的网格分辨率(240～620).第2 列介绍了当Re=1000 时,信息采集点水平速度分量ux的收敛值.第3 列展示了每个网格分辨率相对于SLBM-1024 算法计算结果的相对误差 ∆ux(%).由表可见,随着网格分辨率的增加,相对误差逐渐减小,这与均匀正方形网格数值模拟所表现出来的特性一致.当网格分辨率从500 增加至620 时,相对误差从1.7%减小至0.99%.但是由于网格数量的增加会降低计算效率,同时1.7%的相对误差小于2.0%,是本文可接受的误差范围.综合考虑,本文所有的数值模拟都采用500 的网格分辨率.

表3 不同网格数量计算结果对比Table 3 Comparison of different resolutions

此外,本文还开展了不同插值格式之间的对比研究.同样地我们以传统算法SLBM 且网格分辨率为1024 时的计算结果作为对比.结合表4,可以看到第1 列列举了不同的插值格式,分别为拉格朗日9 点插值(RLBM-L9)、牛顿9 点插值(RLBM-N9)、埃尔米特9 点插值(RLBM-H9)以及拉格朗日25 点插值(RLBM-L25).第2 列介绍了当Re=1000 时,信息采集点水平速度分量ux的收敛值.第3 列展示了每个插值格式相对于SLBM-1024 计算结果的相对误差 ∆ux(%).由表可见,本文RLBM-L25 算法相对于其他插值算法表现优异.

4.2 非定常流动分数值模拟

为了进一步探讨本文算法(RLBM-L25)对于非定常流动的表现,本文开展了非定常周期性流动的数值模拟.选取两个非定常周期性雷诺数8800 和9500 作为数值模拟计算条件.

表5 展示了不同网格分辨率(RLBM-L25-500和RLBM-L25-620)雷诺数为9500 时计算结果,对SLBM-1024 的计算结果.第1 列数据(Data)分别列举了周期性速度曲线的最大值(uxmax)、最小值(uxmin)、平均值(uxmean)以及平均值相对SLBM-1024 计算结果的相对误差 ∆uxmean.通过表5 我们发现,对于非定常周期性流动,随着网格分辨率的增大,相对误差逐渐减小.相对于网格分辨率620,500 的分辨率产生较大的相对误差,但仍然在误差允许范围内(∆uxmean<5%).作者认为,这个网格分辨率可以保证非定常周期性流动的数值模拟精度.

表5 不同网格分辨率计算结果对比Table 5 Comparison of different resolutions

如图11 所示,介绍了流场在一个完整周期(T=1/f,f=0.432 6)内某一时刻的瞬时流线图.其表现出的流场典型特征与作者之前针对方腔内流的研究[41-42]完全一致.

图11 某一时刻周期性解的瞬时流线图Fig.11 Snapshots of periodic solutions at a given moment

图12 介绍了雷诺数为9500 时,信息采集点基于不同插值格式的计算结果对比.主画面展示了速度频谱图对比,可见不同的插值格式对于周期性流动的振荡频率(居中内嵌图)的预测基本一致,且本文RLBM-L25 算法表现优异.从右侧内嵌图可以看出,在预测周期性流动的速度振幅和平均值时,与SLBM-1024 的计算结果存在差异.相较于其他插值格式,RLBM-L25 算法的表现最好,即便如此,与SLBM-1024 的计算结果的相对误差为3.5%左右.作者认为导致这一误差的可能原因有两个.其一,与本文网格细化策略有关,比如网格分辨率,网格间距控制函数等.此外插值格式也是一个需要考虑的因素,作者认为更精确的插值格式能很大程度上消除或降低这一误差.