基于格子玻尔兹曼方法的局部网格加密算法<br/>——粗细网格间的数据转换1)

基于格子玻尔兹曼方法的局部网格加密算法
——粗细网格间的数据转换1)

2023-12-16刘春友李作旭王连平

力学学报 2023年11期

关键词：黏性分辨率加密

刘春友李作旭王连平 ,†, ,2)

* (南方科技大学力学与航空航天工程系,广东深圳 518055)

† (南方科技大学复杂流动及软物质研究中心,广东深圳 518055)

** (广东省湍流基础研究与应用重点实验室,广东深圳 518055)

引言

在实际的工程应用中,各种各样的流体力学问题随处可见,而如何准确且高效的描述其对应的流动特征一直是工业界和科学界重点关注的问题.近半个世纪以来,计算机技术的迅猛发展使得大规模数值求解复杂问题成为可能,而数值模拟作为一种有效的研究手段也得到了更加广泛的关注和应用.不同于直接数值求解纳维-斯托克斯(Navier-Stokes,N-S)方程的传统计算流体力学方法,介观计算流体动力学方法是直接数值求解基于分子动力学理论的Boltzmann 方程,而格子玻尔兹曼方法(lattice Boltzmann method,LBM)则是近30 年来出现的一种较为高效的介观计算流体动力学方法.它较高的计算效率和灵活性使其可以适用于各种高度复杂的物理现象,如湍流[1]、燃烧[2]、多相流[3]、磁流体[4]和粒子悬浮[5-6]等问题.由于其数值耗散较低,故LBM也常被应用于气动声学模拟中[7-8].同时,LBM 也凭借演化过程中“碰撞-迁移”的局部性,使得其非常适合于大规模的并行计算[9-10].除此之外,许多学者为了充分利用LBM 的优点,还通过反设计的思想在原始的Boltzmann 方程中添加一些特定的源项,以使得修改后的LBM 可以满足特定需求[11-12],而这也使得LBM 方法的适用范围得到了极大的扩展.

流体的流动特点一般是随着Reynolds 数的增加,会产生越来越精细的动力学结构.而工业领域里涉及到的流动问题大多是高Reynolds 数的湍流,通过直接数值计算来解析所有的湍流结构往往需要较高的计算成本.因此,为了在有限的计算资源下获得较为准确的结果,往往需要提高局部分辨率来解析这些局部流域内包含的精细流动结构.但由于LBM算法在时空离散时是沿着特征线进行的离散,故在给定离散速度集后LBM 中的时间步长和空间步长便通过离散速度相互耦合,这一特点导致标准的LBM只能使用均匀的直角网格.

为了克服标准LBM 只能使用均匀网格这一缺陷,Filippova 等[13]提出了基于LBM 的局部网格加密(FH)方法.该方法首先需要使用均匀的直角粗网格(背景网格)覆盖整个流场,然后在需要进行加密的区域上对网格进行加密,不同网格分辨率区域下仍使用标准的LBM 演化关系且使用同一离散速度集.而粗细网格之间的耦合则通过流场物性参数以及分布函数间的转换实现,如可以通过保证粗细网格间运动黏度系数一致来得到粗细网格间格子松弛系数的关系,然后通过保证粗细网格间宏观量一致得到分布函数的转换关系.文中以二维圆柱扰流为测试算例,在使用局部网格加密技术后得到了较好的数值结果,这也验证了网格局部加密的有效性.局部网格加密算法的提出对于标准LBM 的意义是重大的,这意味着其可以通过局部加密来节省大量的计算成本,也使得LBM 被大规模应用于工业界成为了可能.

与FH 方法不同,Lin 等[14]建议不对粗细网格间分布函数的非平衡部分做任何缩放,而是直接交换粗细网格的分布函数.Dupuis 等[15]比较了FH 方法和Lin 的方法后发现,当考虑非平衡部分的转换时,可以得到更精确的数值结果.而Touil 等[16]更是直接指出,Lin 的方法实际上是混淆了连续分布函数和离散分布函数之间的区别,其直接在粗细网格间交换分布函数会导致该方法的精度下降.此外,Rohde等[17]提出了一种基于有限体积思想且可以适用于任何碰撞模型的局部网格加密算法,但其也犯了与Lin 相同的错误,即不对分布函数做任何转换,这种处理方式在模拟高Reynolds 数下的流动时,很容易在网格加密界面处出现间断和非物理震荡.而许多学者[18-21]的相关研究也表明,在局部网格加密算法中,粗细网格间非平衡态分布函数的转换是必要且不可或缺的.

众所周知,LBM 的灵活性很大程度上来自于可以通过Chapman-Enskog (CE)展开[22]灵活设计源项.因此,在考虑源项的前提下,建立粗细网格间分布函数的转换关系,对于扩展局部网格加密算法在LBM中的应用范围至关重要.Huang 等[23]基于CE 展开,提出了一种在多松弛碰撞模型(multiple relaxation time,MRT)[24]下利用多层网格技术计算被动标量的模型.该方法在考虑外力项存在的情况下,在矩空间下执行粗细网格间的数据转换.Liou 等[25]同样也基于CE 展开,提出了一种在BGK (Bhatnagar-Gross-Krook)模型[26]下,同时考虑外力和标量源项的加密方法.值得注意的是,这些方法都是在 CE 展开的基础上导出的,其推导过程较为复杂,且在推导中对分布函数的非平衡态使用了一阶CE 近似,而这有可能限制局部网格加密算法在高阶LBM[27]中的应用.本文的主要目标是,在考虑任意源项存在的前提下,提出一套规范且简洁的粗细网格在重合点处的数据转换关系推导方法.

本文的其余部分按如下顺序展开.在第 1 节,梳理BGK 模型和MRT 模型下LBM 的时空离散过程,尤其是指出在时空离散过程中离散分布函数变量的引入,从而为后文局部网格加密算法中的数据转换进行必要的铺垫;同时,梳理BGK 模型和MRT 模型下非平衡态分布函数的初始化方法,为后文数值验证中分布函数的初始化提供依据.第 2 节则是凭借第 1 节中的理论基础,分别详细给出BGK 模型和MRT 模型下粗细网格间数据转换关系的推导过程.在第 3 节中,通过对强迫泰勒-格林涡流动和平板泊肃叶流中对流-扩散问题进行数值模拟,以验证第 2 节中转换关系对复杂源项的适应性;通过对顶盖驱动方腔流动进行数值模拟来展示局部网格加密技术在处理复杂流动问题方面的优势;通过对一维剪切波问题的模拟,来观察局部网格加密技术对LBM 在数值耗散方面的影响.最后,在第 4 节给出本文结论.

1 理论基础

速度离散后的Boltzmann 方程(discrete velocity Boltzmann equation,DVBE)可以表示如下[28]

这里,ξα表示第 α 个离散速度点对应的离散速度;fα(x,t) 代表t时刻在位置x处以速度 ξα移动的粒子对应的分布函数;Ωα代表碰撞项;Sα为可以代表包括体力项在内的任意源项.以二维问题中被广泛使用的D2Q9 速度集为例,如图1 所示,其离散速度点可以表示如下[28]

图1 D2Q9 离散速度模型示意图Fig.1 Schematic diagram of D2Q9 discrete velocity model

不同离散速度点对应的权重为

需要注意的是,本文使用的D2Q9 离散速度集长度未经缩放,即格子声速:Cs=1.

在标准的LBM 中,时空离散是沿着特征线进行的离散,这种离散方法导致了时间步长和空间步长通过离散的速度集进行了耦合.这就使得整个计算域都必须使用均匀的直角网格,以保证流体粒子在一个时间步长 ∆t后可以精确地移动到邻近的网格节点上,这也是LBM 需要使用 on-lattice 速度集的原因.而对于其他的介观方法,如:基于有限体积法对DVBE 经行离散的DUGKS[29]等方法,这些方法中时间步长和空间步长是解耦的,故其可以使用非均匀网格,且既可以使用on-lattice 的速度集也可以使用off-lattice 的速度集.接下来,我们将给出标准LBM 中对DVBE 的时空离散过程直接对DVBE 左右两侧沿特征线积分,可得如下结果

从上式不难看出,分布函数fα从t时刻到t+∆t时刻,时间上推进了 ∆t而空间上推进的长度为∆x=ξα∆t,这也正是LBM 中时空耦合的原因.同时在推导上式的过程中,我们对DVBE 等号右边项使用了梯形积分公式来近似,其中分别代表Ωα(x+ξα∆t,t+∆t),Sα(x+ξα∆t,t+∆t).不难发现式 (4)右边项中包含了t+∆t时刻的未知量,而引入新的变量gα(x,t) 可以消除这种隐式[30],该离散分布函数变量的表达式如下

现在可以使用离散分布函数变量gα对略去高阶小量的式 (4) 重新表示

在给定碰撞项 Ωα以及源项Sα的具体形式后,结合式 (5) 对式 (6) 简单整理后便可以得到对应碰撞模型下的LBM 演化式.下面以BGK 碰撞模型和MRT 碰撞模型为例,给出两种模型下LBM 演化式的详细推导过程.

1.1 BGK 碰撞模型

速度离散后的BGK 碰撞模型可以直接表示为

式中,q为离散速度集中离散速度点的个数.而H(n)(ξα) 代表以 ξα为自变量的第n阶Hermite 多项式[31],前两阶的Hermite 多项式可以表示如下

其中,式 (8) 中的密度 ρ 和速度u可以由fα以及离散速度 ξα表示如下

而根据式 (11)和式(12) 以及Hermite 多项式的具体表达式可知,密度 ρ、速度u以及黏性应力张量σ正好分别对应连续分布函数变量fα的第 0 阶、第1阶Hermite 矩以及非平衡态分布函数的第2阶Hermite 矩

将式 (7)代入到式(5)后,可以反解出

式中,τν代表格子松弛系数,其与网格分辨率有关

由于时空离散前后平衡态部分并没有发生任何改变,故为了统一符号,取将式 (14) 和式(15) 代入式 (7),可以得到BGK 碰撞模型下完全由新变量gα表示的碰撞项Ωα

将式 (16) 代入式 (6),便可以得到使用BGK 碰撞模型的标准LBM 演化方程

而通过Hermite 展开得到的力项表达式 (18),也与Guo 等[32]通过反设计得到的结果完全一致.

除此之外,由于式 (17) 的解gα≠fα,故由离散分布函数变量gα计算宏观量的表达式需要重新给出.通过式 (14),可以得到连续分布函数变量fα对应的各阶Hermite 矩a(n)与离散分布函数变量gα的各阶Hermite 矩间的关系

进而可以通过式 (13)、式(20) 以及式 (24) 得到当源项为体力项时,由离散分布函数变量gα表示的密度 ρ、速度u以及由离散非平衡态分布函数变量表示的黏性应力张量σ

最后,在经过合适的初始化后,式 (17) 便可以显式推进.而当初始流场中包含非平衡态时,正确的初始化分布函数对于一些非定常问题是至关重要的.其中平衡态分布函数可以直接根据式 (8) 通过初始场的密度以及速度给出,而非平衡态分布函数的初始化可以借由Hermite 展开给出[33].对于等温的流动问题,分布函数的非平衡态部分可以直接使用其前两阶Hermite 展开来近似,即

将式 (24) 代入到上式中,便可以得到初始流场在BGK 模型下且考虑源项为体力项时,非平衡态分布函数的具体表达形式

由于我们的目的是保证初始时刻分布函数所对应的密度、速度以及黏性应力张量与初始条件一致,故这里只是展到二阶.

以上便是从DVBE 到LBM 的详细推导过程.其中注意的是,在离散过程中引入了离散分布函数变量gα,而且核心演化关系式 (17) 实际求解的也是离散分布函数变量gα,而不是时空离散前的连续分布函数变量fα.此外,从式 (5) 不难看出,离散分布函数变量gα与时间步长 ∆t相关,进而也与网格分辨率有关.故当使用式 (17) 作为计算流域上的演化公式时,即使忽略时空离散误差,同一时空点下不同网格分辨率对应的gα仍然是不一致的,所以粗细网格重合点上的gα并不能直接交换,而这也是局部网格加密算法中不同网格分辨率在重合点处离散分布函数gα需要转换的原因[16].

1.2 MRT 碰撞模型

MRT 模型是d'Humieres[24]在1992 年提出的一种碰撞模型,该模型可以被视为广义的BGK 模型.相比于只有一个松弛参数的原始BGK 模型,MRT模型中则有着多个松弛参数用以描述不同物理量趋向平衡态的过程.除了少数几个有着明确物理含义的松弛参数外,MRT 模型还有若干个可调参数,相关文献[34-35]在对MRT 模型进行系统研究后指出,对这些参数进行合理的调节可以有效改善数值稳定性以及减小边界条件的误差.

MRT 模型对应的速度离散后的碰撞项可以表示如下[36]

其中 Λ 代表物理松弛系数矩阵,其与网格分辨率无关;而的表达式仍然使用式 (8).上式中使用了与D'Humieres 相同的记号,以“〈·| ”,“|·〉 ”来分别表示行向量以及列向量,即

这里上标“T ”代表转置操作符.将式 (28) 代入到式(5) 的矩阵形式中后,可以反解出由离散分布函数变量 |g〉所表示的连续分布函数变量 |f〉如下

其中I为单位矩阵;上标“-1 ”代表矩阵求逆操作符.进而可以得到MRT 模型下由离散分布函数变量|g〉表示的碰撞项

将式 (30) 代入到式 (6) 的矩阵形式中,可以得到

若记

则将上式代入到式 (31) 中后,便可以得到MRT 模型下关于离散分布函数变量 |g〉的演化方程

其中,M是一个q×q的正交变换矩阵;M-1为M的逆矩阵,且满足以下关系

式中,|m〉代表转换矩阵M下分布函数 |g〉对应的一组矩.对式 (33) 两边同时左乘转换矩阵M后便可得到MRT 模型下关于 |m〉的演化方程[23]

这里 |ms(x,t)〉=M|S(x,t)〉.而转换矩阵M与速度集的选取有关,以D2Q9 速度集为例,其对应的转换矩阵M为[34]

其中,sν与运动黏度系数 ν,以及se与体积黏度系数ζ的关系可以由CE 分析给出[37]

进而我们可以通过式 (13)、式(19)、式(39 a)以及附录A 中MH得到MRT 模型下当源项为体力项时,由离散分布函数变量gα表示的密度 ρ、速度u

上述结果也与Chen 等[37]通过CE 分析得到结果保持一致.

当初始流场中包含非平衡态时,MRT 模型下的分布函数也需要正确的初始化.同样可以使用Hermite 展开来初始化MRT 模型下的分布函数,其中平衡态分布函数的初始化与BGK 模型一致,而离散非平衡态分布函数的初始化,则需要将附录A 中 |gneq〉的前6 个矩代入到式 (26),进而可以得到MRT 模型下当源项为体力项时,初始流场中的具体表达式

上式中各参数的具体形式可以参见附录A.

此外,由于矩阵M和矩阵初始化后就已确定,且其在式 (33) 演化过程中不再改变,故可以在初始化时将计算后单独存储起来以避免在执行碰撞步骤时重复计算,从而节省计算量.且直接使用式 (33) 去执行分布函数的碰撞步骤相比于使用式(34) 和式 (35) 更加节省时间,经我们程序测试后发现CPU 消耗时间减少 25% 左右.

2 局部网格加密中的数据转换

如前所述,同一点处不同网格分辨率下的离散分布函数变量gα之间的转换是必要且不能忽略的.使用局部网格加密算法需要保证,当粗细网格均使用同一离散速度集且在忽略时空离散误差后,粗细网格区域在同一时空点下数值求解后得到的连续分布函数变量fα必须是一致的,而且不同分辨率下同一点处的物性参数如:物理松弛系数 τ 以及运动黏度系数 ν 等也必须保持一致.

LBM 在对式 (1) 做时空离散时为了消除隐式,引入了与分辨率有关的离散分布函数变量gα以及相应的格子松弛系数 τν(BGK 模型)和(MRT 模型)

根据上式可知,不同网格分辨率下的同一时空点上,格子松弛系数 τν和以及实际求解得到的离散分布函数变量gα和 |g〉是不一致的,故粗细网格在交换这些量之前需要经过合适的转换.与之相反的是,物理松弛系数 τ 以及 Λ-1是不随网格分辨率变化的;而当考虑粗细网格区域均使用同一离散速度集时,则数值求解得到的连续分布函数变量fα和 |f〉,在不考虑时空离散误差的前提下,也是不随网格分辨率变化的.故接下来将以此为桥梁,建立起BGK模型和MRT 模型下,不同网格分辨率间的数据转换关系.

下文中,与粗网格和细网格区域相关的任何量分别用上标“c”和“f”表示.设局部加密时粗细网格的空间步长比例为

局部网格加密算法中不同网格分辨率下的区域仍然使用的是同一离散速度集,故在考虑对流尺度下,粗细网格的时间步长也有同样的比例关系

故为了保证粗细网格上物理时间的统一推进,粗网格上运行一步则细网格上需要运行n步.故n常设为大于 1 的整数.

2.1 BGK 模型下的局部网格加密

为了将我们的方法与现有推导过程更好地对比,在给出新的方法前,先梳理一下Liou 等[25]在考虑源项时,分布函数转换关系的推导过程.

首先对式 (17) 使用泰勒展开

考虑CE 展开如下

其中 ϵ 为与克努森数Kn同阶的小量.将式 (47) 代入到式 (46) 中,便可以得到与 ϵ0,ϵ1,ϵ2同阶的方程

而式 (48b) 等号左边项均为宏观量与离散速度ξα的函数,故当粗细网格区域使用同一离散速度集时,其与网格分辨率无关,故如下关系式成立

而格子松弛系数 τν的转换则是通过保证动力黏度系数 µ=p(τν-0.5)∆t一致的基础上推导出来的,即

以上便是Liou 等在考虑源项存在时,给出的BGK 模型下粗细网格间离散非平衡态分布函数变量以及格子松弛系数 τν转换关系的推导过程,而转换关系的推导过程本质上是利用CE 展开寻找与网格分辨率无关的量.不难发现整个推导过程是较为复杂的,而且使用了来近似对于恢复N-S 方程而言,这样的近似是足够的,但对于需要捕捉稀薄效应的高阶LBM 方法[38]而言,一阶近似显然是不够的.下面将给出一种更加简洁且直观的推导方法.

如前所述,BGK 模型下的物理松弛系数 τ 和连续分布函数变量fα与网格分辨率无关,故由式 (15)和式 (14)可知如下关系是成立的

对上式简单整理后,便可以得到不同网格分辨率在重合点处的数据转换关系如下

将式 (52b) 二边减去平衡态并利用式 (52a),便可以得到与式 (50) 一致的转换关系式.

相比之下,我们的推导过程中除忽略时空离散误差外,并没有做任何额外的假设(没有对非平衡态分布函数做近似假设),故上述转换关系同样保证了粗细网格重合点处高阶非平衡态分布函数的一致.然而,Liou 等在使用的假设后,也得到了与我们相同的结果.而这主要是由于各阶非平衡态分布函数间的递推关系导致的(详细推导过程见附录 B)

上述推导过程不仅解释了我们的推导结果与前人结果一致的原因,同时也为之前学者对非平衡态分布函数做一阶近似提供了理论依据,而且上述推导过程也表明了局部网格加密算法同样可以被扩展到高阶LBM 中.

2.2 MRT 模型下的局部网格加密

下面将继续使用新方法给出MRT 模型下,不同网格分辨率在重合点处离散分布函数变量 |gc〉和|gf〉以及格子松弛系数矩阵间的转换关系.与上一节BGK 模型下的推导过程类似,不同分辨率下的物理松弛系数矩阵 Λ-1应保持一致,故由附录A 可得

而由CE 分析可知,sq和sϵ并不会直接出现在N-S 方程中,故可以被视为自由参数,但其对数值稳定性及边界条件误差有一定影响[34].Chen 等[37]认为这两个参数在不同网格分辨率下可以不需要转换,即但本文中为了保持转换关系的一致性,仍然采用式 (58) 作为粗细网格间sϵ和sq的转换关系式.

同样,直接从式 (29) 出发,基于连续分布函数变量 |f〉在不同网格分辨率下的重合点处保持一致,便可以得到MRT 模型中粗细网格在时空重合点处的离散分布函数变量 |g〉的转换关系

且该结果与Huang 等[23]的结果也完全一致.

综上,我们以保证不同网格分辨间连续分布函数变量fα,|f〉以及物理松弛系数 τ,Λ-1一致为基础,以更加直观且简洁的方法分别推导了BGK 模型和MRT 模型下,不同网格分辨率在时空重合点处离散分布函数变量gα,|g〉和格子松弛系数 τν,间的转换关系,其相比于之前学者[23,25]的推导过程更加简洁且直观.

3 局部网格加密算法的数值验证

本节数值模拟中将采用与Yu 等[18]相同的网格排布方式,其示意图如图2 所示.同时,为了避免局部网格加密算法中插值误差对数值模拟的影响,我们使用与Gendre 等[39]相同的高阶插值公式.而对于加密界面处缺失的离散分布函数,采用Chen 等[40]提出的只使用空间插值的方法.

图2 粗细网格排布示意图.灰色区域为粗细网格的重叠区域Fig.2 Schematic diagram of the coarse and fine grid arrangement.The grey area is the overlap between the coarse and fine grids

下面将通过使用局部网格加密技术的LBM来对强迫泰勒-格林涡流动、平板泊肃叶流中对流-扩散问题、顶盖驱动方腔流动和一维剪切波问题进行模拟,在前两个算例中均包含非均匀且非定常的源项,故以此来验证粗细网格间数据转换关系对复杂源项的适应性;而有着较为复杂流动现象的顶盖驱动方腔流动则可以较好地验证局部网格加密技术的有效性;一维剪切波问题的模拟则是用来检验局部网格加密技术对原本的LBM 算法在耗散方面的影响.

3.1 强迫泰勒-格林涡流动

在标准的 Taylor-Green vortex 流动中没有外力的参与,因此流动会随时间单调衰减.为了验证源项存在时,粗细网格间数据转换关系的正确性,我们以Min 等[41]使用的二维强迫泰勒-格林涡流动作为验证算例.该流动中不仅有体力的参与,且该体力非均匀且非定常,故可以较为充分地验证BGK 模型以及MRT 模型下转换关系的正确性.

该流动问题对应的控制方程为

流动的物理区域被设为:(x,y)∈(0:Lx,0:Ly),4 周为周期边界条件,考虑单位质量受到的体力为b=k2ν(1-Q)u时,则该流动对应的理论解如下

式中,U0代表初始速度的幅值;kx=2π/Lx和ky=2π/Ly分别代表x和y方向上的波数P0代表参考压力(可以为任意值).为了简单,我们采用kx=ky,Lx=Ly=L=1 以及P0=1.此外,Q代表衰减因子,可以通过调节参数Q来控制流动的衰减速度.

(1)Q=1 时,代表没有外力,各项宏观量按标准的泰勒-格林涡流动衰减;

(2)Q=0.5 时,代表各项宏观量的指数衰减率仅为体力不存在时的一半;

(3)Q=0 时,体力输入的能量与黏性耗散平衡,流场中各项宏观量与初始流场保持一致且不随时间变化;

(4)Q=-0.5 时,体力输入的能量大于黏性耗散,流场中的压力和速度等宏观量的绝对值将随着时间的增加而增加.

使用t=0 时的解作为流场各宏观量的初始值,并根据式 (8) 对平衡态分布函数进行初始化;且由于该流动问题对应的初始流场中包含非平衡态部分,故需要根据式 (27) 和式 (42) 分别对BGK 模型和MRT 模型下的离散非平衡态分布函数初始化;而两者之和便可以作为初始流场对应的离散分布函数变量gα.

粗细网格上均采用D2Q9 离散速度集.初始速度幅值U0对应的Mach 数为 0.01;Reynolds 数Re=U0L/ν=100;网格数:Nx=Ny=100.加密区域范围为 (x/L,y/L)∈(0.4:0.6,0.1:0.4),加密比例为2.MRT 模型下各松弛系数如下

在细网格中,MRT 模型下的各松弛系数直接使用式 (57) 及式 (58) 转换而来.此外还需注意的是原来在计算速度u时,有如下公式(考虑源项为体力项)

但本算例中b=k2ν(1-Q)u,其大小与u直接相关,故需要将b代入到式 (63) 中来显式得到u的计算公式

为了充分验证粗细网格间数据转换关系的正确性,首先分别绘制了BGK 模型和MRT 模型在Q=0.5 以及tU0/L=1 时,相对压力P-P0、水平方向速度ux以及黏性应力 σxx的分布云图.数值计算结果如图3 所示,图中的黑色方框为网格加密的界面,其所包围的区域为加密区域.从图中可以看出,各宏观量在加密界面两侧均光滑无间断.

图3 强迫泰勒-格林涡各宏观量的分布云图(Q=0.5, tU0/L=1,黑色矩形框内的区域为加密区域)Fig.3 The contours of each macroscopic quantity of the forced Taylor-Green vortex (Q=0.5, tU0/L=1,the area inside the black rectangular box is the refinement region)

此外,还分别计算了当Q为 0.5,0.0 和 -0.5 时,在 (x/L,y/L)=(0.4,0.25) 点(粗细网格重合点)处,相对压力P-P0、水平方向速度ux以及黏性应力σxx随时间的变化,其数值结果如图4 所示.在每张子图中,绘制了在Q取不同的值时,BGK 模型以及MRT 模型下数值结果与理论解果的对比.其中实线代表由式 (61) 给出的理论解,圆圈和倒三角分别代表BGK 模型和MRT 模型下的数值结果,而不同的颜色代表Q取不同的值.观察图4 后发现,两种碰撞模型的数值结果都与理论解吻合得非常好.

图4 强迫泰勒-格林涡中在 (x/L,y/L)=(0.4,0.25) 点处各宏观量随时间的变化.(a)～(c)分别为BGK 模型和MRT 模型下的相对压力P-P0、水平方向速度 ux 和黏性应力 σxx 与理论结果的对比Fig.4 The variation of each macroscopic quantity with time at the point (x/L,y/L)=(0.4,0.25) in the forced Taylor-Green vortex.(a)～(c) are comparisons between the relative pressure P-P0,velocity in the x direction ux and viscous stress σxx under BGK model and MRT model,respectively,and the theoretical results

最后,为了进一步观察粗细网格界面两侧宏观量的连续性,我们还绘制了当Q为 0.5,0.0 和-0.5以及tU0/L=1 时,在y/L=0.25 处BGK 模型和MRT模型下各宏观量随x/L变化的剖线图,同时也与理论解进行了对比.其数值结果如图5 所示,图中两条黑色虚直线间的区域为加密区域.从图中可以看出,两种碰撞模型都取得了非常好的数值结果,且各宏观量在加密界面两侧光滑连续.

图5 强迫泰勒-格林涡中各宏观量的剖线图(tU0/L=1, y/L=0.25,两条黑色虚线间的区域为加密区域).(a)～(c)分别为BGK 模型和MRT 模型下相对压力 P-P0、水平方向速度 ux 和黏性应力 σxx 与理论结果的对比Fig.5 Profile diagram of each macroscopic quantity of forced Taylor-Green vortex (tU0/L=1, y/L=0.25,the area between the two black dashed lines is refined).(a)～(c) are comparisons between the relative pressure P-P0,velocity in the x direction ux and viscous stress σxx under BGK model and MRT model,respectively,and the theoretical results

BGK 模型和MRT 模型下良好的数值结果也印证了第 2 节中粗细网格间数据转换关系的正确性.接下来将以平板泊肃叶流中对流-扩散问题为验证算例,进一步测试BGK 模型下粗细网格间数据转换关系对复杂源项的适应性.

3.2 平板泊肃叶流中对流-扩散问题

如前所述,合理设计LBM 中的源项S可以极大地扩展该方法的适用范围,许多学者通过CE 分析反设计源项S的形式,以使得修改后的LBM 可以满足特定要求[11].甚至在很多情况下,LBM 被当作一款求解一些特殊偏微分方程的高效通用求解器[42].其中,对流-扩散方程作为一种可以描述由于对流和扩散过程而引起的传递热量、质量或其他物理量输运现象[43]的特殊偏微分方程,被广泛应用于求解各种对流-扩散问题.该方程可以被表示如下

式中,ϕ(x,t) 代表一个标量,u(x,t) 是由N-S 方程控制的速度,D为扩散系数.

对于如何使用LBM 求解对流-扩散方程,之前的学者们利用CE 分析给出了许多不同的方案[44-48].这里使用Chai 等[48]提出的通过CE 分析反设计源项进而恢复对流-扩散方程的方案,该方法的标量场中分布函数的演化方程如下

这里 ϕα对应标量场中的分布函数,而Rα则是为了在恢复式 (65) 时减小额外计算误差所添加的源项,且和Rα分别表示如下

这里 ∇ϕ 可以直接使用如下公式计算

接下来将用上述方案去求解平板泊肃叶流中对流-扩散问题.在计算这个问题时,需要两套分布函数:gα用于计算速度场中的密度 ρ 以及速度u;ϕα用于计算标量场中的 ϕ.

该问题对应的计算简图如图6 所示.其中上下板处的灰色区域为加密区域,加密范围为:y/H∈{(0.0:0.1)∪(0.9:1.0)},加密比例为 2.其他参数为:定常时中心速度U0对应的Mach 数为 0.01;下板标量:ϕb=0.0,上板标量:ϕt=0.1;Reynolds 数Re=HU0/ν=10,H为槽道宽度;Peclet 数Pe=HU0/D=10;网格数Nx=6,Ny=20,粗细网格均采用D2Q9 速度集以及BGK 碰撞模型.该流动采用水平方向的体力驱动,单位质量所受的体力为:bx=8νU0/H2.左右为周期边界条件,上下边界均使用非平衡态外推边界条件[49],其中边界点处密度的计算使用非平衡态反弹边界条件[50]中密度的计算方法.

图6 平板泊肃叶流中对流-扩散问题示意图Fig.6 Schematic diagram of diffusion problem in a planar Poiseuille flow

该流动问题中速度ux(y,t) 所满足的控制方程为

对应的解析解可以通过分离变量法求得,其级数形式可以表示为

式中k=2n+1.此外,标量 ϕ (y,t) 所满足的控制方程可以表示如下

其对应的解析解也同样可以通过分离变量法给出,且级数形式表示如下

其中J=-D∇ϕ+ϕu代表扩散通量与对流通量之和,而Jx,Jy分别为x,y方向上的分量.

我们使用局部网格加密技术,分别计算了在t∗=tν/H2=0.04,0.1,1.0 时,水平方向的速度ux、标量 ϕ、通量Jx和Jy并与通过式 (71)和式(73) 得到的理论结果进行了对比.数值结果如图7 所示,图中的黑色虚线代表粗细网格的界面,实线代表理论解,圆圈代表数值解,不同的颜色代表不同的时刻.从图中可以看出,不同时刻下的数值结果都与理论解吻合得非常好,而且在加密界面的两侧各宏观量均光滑且连续.良好的数值结果也进一步验证了第 2节中粗细网格间数据转换关系的正确性,同时也验证了该转换关系对复杂源项的良好适应性.

图7 平板泊肃叶流中对流-扩散问题(黑色虚线为加密界面).(a)～(d)分别为BGK 模型下水平方向速度 ux、标量 ϕ、通量 Jx 和 Jy 与理论结果的对比ux,scalar ϕ,flux Jy and Jx compared with theoretical results under the BGK modelFig.7 Diffusion problem in a planar Poiseuille flow (the black dashed lines are the refinement interface).(a)～(d) represent the velocity in x direction

3.3 顶盖驱动方腔流动

顶盖驱动方腔流作为经典的基准算例已被广泛用作数值方法的测试案例.为了对数值结果进行合理的评价,使用Tamer 等[51]的数值结果作为参考数据并与之进行对比.使用LBM 算法去计算顶盖驱动方腔流时,如果使用的边界条件较为粗糙且网格分辨率较低时,顶盖左右两个角点处往往会产生明显的压力“噪声”.

下面使用局部加密算法对该流动进行数值模拟,用以进一步展示局部加密算法的优势.本次模拟使用BGK 模型以及D2Q9 速度集,雷诺数Re=UwL/ν=1000,其中L代表顶盖长度,Uw为顶盖速度其对应的马赫数为 0.1,ν 为运动黏度系数.边界条件:左右壁面以及下壁面均使用修正反弹格式,顶盖处使用运动边界对应的反弹格式[52],即

图8 顶盖驱动方腔流4 个角点示意图,其中黑色粗实线代表壁面Fig.8 Diagram of four corner points of the square cavity flow driven by the top lid,where the thick black line represents the wall surface

其中 (xb,yb) 代表右上角点对应的坐标,而 2,3,6 方向均使用静止边界的反弹格式.左上角点也是同样处理.

对整个顶盖附近实施加密,加密范围为:y/L∈(0.7:1.0),加密比例为 2.不同分辨率下压力分布云图如图9 所示.其中,不仅对比了使用局部加密(UCG-L)和使用均匀粗网格(分辨率与局部加密中的粗网格分辨率一致,UCG)以及均匀细网格(分辨率与局部加密中的细网格分辨率一致,UFG)下的结果,同时还添加了与UCG-L 达到稳态时理论计算耗时相同的均匀网格(EQ)下的结果,用来对比使用局部网格加密算法的非均匀网格与相同计算量下的均匀网格两者的数值结果.从图9 中不难看出,采用局部加密后原来的“噪声”得到了极大的消除,而且其对“噪声”的抑制效果比相同计算耗时下的均匀网格更好.这也侧面验证了局部网格加密方法在消除局部“噪声”的有效性.值得一提的是,Dong 等[35]通过详细的分析指出,这里的压力噪声主要是边界条件隐藏误差导致的,且可以通过利用双松弛时间碰撞模型优化加以有效消除,其文章中的数值结果也证明了这一点.

图9 Re=1000 时顶盖驱动方腔流的压力分布云图对比,(a)～(d)分别为UCG,UFG,EQ 和UCG-L 对应的压力分布云图 ((d)中的黑色直线为加密界面,界面以上为加密区域)Fig.9 When Re=1000,the contours of pressure of the lid-driven cavity flow are compared.(a)～(d) are the contours of pressure corresponding to UCG,UFG,EQ and UCG-L,respectively (the black line in (d) is the grid refinement interface,and the area above the black line is the refinement area)

图10 中绘制了使用不同分辨率时沿垂直方向穿过几何中心的水平速度ux随y/L的变化关系以及沿水平方向穿过几何中心的垂直速度uy随x/L的变化关系.而且从图10(b) 以及图10(d) 中可以很清楚地看到即使是远离加密区域UCG-L 的数值结果仍要好于与其所需计算量相同的均匀网格下的结果,而且在使用更少计算资源的情况下还取得了与UFG 非常相近的数值结果,这就表明在局部流域执行网格加密可以很好地减小局部数值误差对全流域的影响.

图10 Re=1000 时顶盖驱动方腔流中无量纲速度的剖线图,UCG,UFG,EQ 和UCG-L 与参考数据的对比.(a) 沿垂直直线穿过几何中心的ux/Uw,(b)为其局部放大,(c) 沿水平直线穿过几何中心的 uy/Uw,(d)为其局部放大Fig.10 Profile diagram of dimensionless velocity in a lid-driven cavity flow when Re=1000,comparison of UCG,UFG,EQ and UCG-L with reference data.(a) ux/Uw in a vertical straight line through the center of the geometry,(b) local amplification of (a),(c) uy/Uw in a horizontal straight line through the center of the geometry,(d) local amplification of (c)

在表1 中,我们比较了3 种情况下每 1000 个迭代步的 CPU 实际消耗时间.如果假设在 UCG 上运行到稳态所需的 CPU 时间是t,若考虑对流时间尺度以及二维情况下则UFG 下运行到相同状态下,理论上所需的 CPU 时间是 8t,而上述加密的UCG-L只需要 3.4t,表1 中的测量结果也基本支持上述论断.值得注意的是,在本算例中UCG-L 所需的计算时间还不及 UFG 的一半,但其数值结果却与 UFG非常相近.有理由相信在某些特殊算例中通过合理布置加密区域,局部网格加密可以实现更高的加速比.

表1 不同分辨率下程序演化1000 步时CPU 消耗时间对比Table 1 Comparison of CPU consumption time when the program evolves 1000 steps at different resolutions

3.4 一维剪切波问题

为了观察局部网格加密技术对LBM 算法在数值耗散方面的影响,我们选择了一维剪切波问题[53]作为研究算例,其示意图如图11 所示.该问题对应的控制方程及初始化流场如下

图11 一维剪切波问题示意图Fig.11 Schematic diagram of one-dimensional shear wave problem

式中,u0代表初速流场中速度的幅值.由于物理黏性的存在,速度ux会随着时间逐渐衰减,其对应的理论解为

针对该问题数值模拟的各参数和设置如下:上下使用周期边界条件;粗细网格上均采用BGK 碰撞模型和D2Q9 速度集;非平衡态分布函数的初始化仍然使用式 (27).初始速度幅值u0对应的Mach数为 0.1;Reynolds 数Re=Hu0/ν=10,H为槽道宽度;网格数Nx=4,Ny=20.分别计算了加密区域范围为:(x/L,y/L)∈(0.0:1.0,0.1:0.25),(x/L,y/L)∈(0.0:1.0,0.1:0.4),(x/L,y/L)∈(0.0:1.0,0.1:0.6) 和(x/L,y/L)∈(0.0:1.0,0.4:0.6),且加密比例均为2 的情况下,UCG-L 和UCG 下的数值耗散 νN与物理黏性ν之比,其中数值黏性可以通过反解式 (78) 求得

分别绘制出UCG-L 和UCG 下的 νN/ν 随着y/H变化的剖线图,其数值结果如图12 所示.图12 展示了在不同区域使用局部网格加密后,在给定时刻(t∗=tν/H2=4π2/ln2) 下的数值黏性与物理黏性之比(红色线)以及均匀粗网格下的数值黏性与物理黏性之比(蓝色线),图中的黑色虚线代表粗细网格的分界线.从图中可以发现,不同加密区域下,加密界面两侧的数值黏性相比于均匀网格下,表现出了不一样的结果,既有增大的情况也有降低的情况.而其中数值黏性小于物理黏性现象的出现可能与加密界面附近的非物理震荡[19,54-55]有着密切的关系,故还需要一系列更加系统的理论研究将局部网格加密技术引起的数值耗散和数值色散与加密区域附近的非物理震荡联系起来,而这也将作为我们后续工作的重点.

图12 一维剪切波问题在不同加密区域下均匀网格和非均匀网格对应的数值黏性和物理黏性之比的分布情况(黑色虚线为加密界面).(a)～(d)分别为不同加密区域下UCG-L 和UCG 对应的数值黏性和物理黏性之比随 y/H 变化的剖线图Fig.12 For one-dimensional shear wave problem,the distribution of the ratio of numerical viscosity and physical viscosity corresponding to uniform and non-uniform grids in different refinement regions (the black dashed line is the refinement interface).(a)～(d) are the cross-section plots of the ratio of numerical viscosity and physical viscosity corresponding to UCG-L and UCG with y/H in different refinement regions,respectively

图12 一维剪切波问题在不同加密区域下均匀网格和非均匀网格对应的数值黏性和物理黏性之比的分布情况(黑色虚线为加密界面).(a)～(d)分别为不同加密区域下UCG-L 和UCG 对应的数值黏性和物理黏性之比随 y/H 变化的剖线图 (续)Fig.12 For one-dimensional shear wave problem,the distribution of the ratio of numerical viscosity and physical viscosity corresponding to uniform and non-uniform grids in different refinement regions (the black dashed line is the refinement interface).(a)～(d) are the cross-section plots of the ratio of numerical viscosity and physical viscosity corresponding to UCG-L and UCG with y/H in different refinement regions,respectively (continued)

4 总结

本文在考虑任意源项存在的前提下,提出了一套规范且简洁的粗细网格在重合点处的数据转换关系推导方法,该方法既可以用于BGK 模型也可以适用于MRT 模型,同时也容易被推广到其他碰撞模型下的数据转换.该推导方法从保证粗细网格重合点处对应的连续分布函数一致出发,除忽略离散误差外,并不依赖于CE 展开以及非平衡态近似.这也说明,局部网格加密算法不仅适用于连续流,而且也适用于高阶非平衡态流动问题.此外,我们还从理论上证明了保证粗细网格间非平衡态部分的一阶CE 近似一致,便可以保证整个非平衡态部分的一致,这也给之前学者对非平衡态部分做一阶CE 近似提供了理论依据.同时,我们还通过引入新的Hermite 转换矩阵MH直接构建了MRT 模型下,时空离散前后分布函数对应的各阶Hermite 矩间的关系,从而将BGK 模型下基于Hermite 展开的初始化方法推广到了MRT 模型下.

为了验证包含复杂源项的粗细网格间数据转换关系,考虑了两个具体算例.第1 个算例是外力作用下的不可压强迫泰勒-格林涡流动,这个问题涉及非均匀非定常流场,并且由源项表示的外力也为非均匀场.初始流场对应的分布函数采用Hermite 展开来准确构建.流场内部的部分区域采用细网格.计算结果显示粗细网格边界处压力、速度、应力场均光滑衔接,这些宏观量随时间和空间的演化剖面都与解析解完全吻合,并且BGK 碰撞模型和MRT 碰撞模型下的计算结果也基本完全一致,证明了我们推导的基于两种碰撞模型转换关系的正确性.第2 个算例是平板泊肃叶流驱动的非定常对流-扩散问题,这里标量场对应的Boltzmann 方程需要引进源项,才能准确恢复宏观非定常对流-扩散方程.在上下两个壁面附近,我们采用细网格局部加密.本问题的标量场在流向有对流通量,在垂直壁面方向有扩散通量,并且它们分布都不均匀,且随时间演化.通过分离变量法,给出了标量场的非定常解析解.计算结果显示速度场、标量场和标量通量的两个分量在粗细网格界面都光滑过渡,并与解析解完全吻合,进一步验证了一般源项下转换关系的正确性.此外,我们还使用局部网格加密技术对经典的顶盖驱动方腔流动问题进行了数值模拟,数值结果表明局部网格加密技术在处理局部压力噪声方面有着明显的优势,而且通过对比使用局部加密的非均匀网格和不使用局部加密的均匀网格的计算耗时以及宏观速度的剖线图,可以较为清楚地看到局部网格加密技术在处理复杂流动问题时的优势.最后,为了观察局部网格加密技术对LBM 算法在数值耗散方面的影响,我们通过在一维剪切波中的不同区域使用局部网格加密,并将其数值黏性与均匀粗网格下的数值黏性进行了对比,发现由局部加密引起的数值黏性与加密区域的选取有很大的关系,且在部分加密区域下还出现了数值黏性小于物理黏性的情况,一系列更加精细的算例和理论分析需要被考虑来阐述该现象的发生.

在上述的数值验证中,我们的算例均为二维算例,并没有涉及更加复杂的三维流动问题,这主要是由于LBM 中的局部网格加密技术经过20 多年的发展已经趋于成熟,相关的文献中[25,37,56]均已报道了大量使用局部网格加密技术计算复杂三维流动问题的测试结果,相关的数值结果均展示了局部网格加密技术在三维复杂流动问题中的优势及适应性.本文的主要创新点在于从粗细网格间分布函数需要转换的本质出发,构建了一套直观且简洁的包含任意源项下分布函数转换关系的推导过程,且我们的推导结果与前人的结果完全一致,这一点在 2.1 中也给出了解释.而如上所述,LBM 中的局部网格加密技术已经被广泛应用到各种复杂流动问题的计算中,且其使用的分布函数转换关系与我们的也完全一致,故本文中并没有对更加复杂的三维流动进行模拟,而是借用二维情况下的强迫泰勒-格林涡流动和平板泊肃叶流中的对流-扩散问题,着重考察了分布函数转换关系对复杂源项的适应性,以及通过对顶盖驱动方腔流动的模拟初步展示局部网格加密技术在处理复杂流动问题方面的优势.

同时需要注意的是,粗细网格重合点处的分布函数转换关系式 (48)和式(49) 中的粗细网格分辨率之比n理论上可以采用任意正整数,但本文粗细网格间分布函数的转换关系是在忽略时空离散误差的前提下得出的,而实际情况下由于离散误差的存在,不同网格分辨率下得到的数值结果中所包含的离散误差也是不一样的,故粗细网格界面处的数值结果间往往会产生一定的“数值间断”,且该“数值间断”随着加密比例的增加而增加,对于一些较为复杂的流动,这往往会在粗细网格界面产生非物理的震荡,故上面的算例中局部网格的加密倍数均为 2,而且绝大多数文献的网格加密比例在计算中往往也被限定为 2.

必须指出的是,本文的分析并没有考虑粗细网格间不同数值离散误差有可能在网格加密界面处产生的数值震荡.这个问题在已有文献中有一些讨论[19,57],但需要进一步做理论分析并找到消除数值震荡的有效方法.其中,不同碰撞模型对粗细网格界面的数值稳定性有着不同的影响[19],但也需要从理论上进一步做分析.如果网格加密边界与物理壁面相交,那么物理壁面的反弹格式局部误差[35]与粗细网格界面误差交汇,此时往往会在壁面附近产生一定的数值震荡[54],该问题的理论分析会更具挑战性,但对实际应用问题却非常重要.这些都是需要进一步研究的问题.

此外,测试和分析局部网格加密算法在数值色散、数值耗散方面的表现对于该技术的实际应用有着非常重要的意义,尤其是有助于加密界面处非物理震荡问题[19,54-55]的解决.在LBM 的局部网格加密算法中,粗细网格上仍然使用标准的LBM 算法,而不同的是需要在粗细网格交界面处进行合理的数据转换.其中,数据转换时往往涉及未知分布函数的构建,目前构建方式可以分为有限差分方法和有限体积方法两类,而不同的构建方式对局部网格加密算法的数值耗散和数值色散都有不同的影响.故研究局部网格加密算法对LBM 在数值色散和数值耗散方面的影响,不仅需要详细对比不同的构建格式在数值耗散和数值色散方面的表现,还需要从理论上分析不同构建格式间结果差异的原因,以及数值耗散和数值色散对加密界面处非物理震荡的影响,而这也将是我们下一步的工作重点.考虑到该问题的复杂性,故我们将在后续文章中尝试对该问题做较为完整和详细的介绍.

附录A MRT 模型下时空离散前后的分布函数对应的各阶Hermite 间的关系

首先,需要对式(29)两边同时左乘Hermite多项式矩阵H

为了让上式可以直接反映连续分布函数变量 |f〉的各阶Hermite 矩与离散分布函数变量 |g〉各阶Hermite 矩间的关系,上式中的Hermite 多项式矩阵H可以直接由Hermite 多项式组成,以D2Q9 速度集为例

上式中各阶Hermite 多项式的具体形式可以表示如下[58]

此时,式 (A1) 中的H|f〉可以表示如下

进而式 (A1)可以表示为

上式直接给出了连续分布函数变量 |f〉和离散分布函数变量 |g〉对应的Hermite 矩间的关系.在对式 (32) 两边取逆后,可以得到

若记Hermite 转换矩阵MH为如下形式

结合式 (A7) 与式 (32) 后,式 (A6) 可以被进一步简化如下

D2Q9 速度集下的Hermite 转换矩阵MH可以显式表示为

上述Hermite 转换矩阵MH中各参数表示的具体表达式如下

当考虑源项为体力项时,|gneq〉的前6 个矩的具体形式如下