CN110705183A - 一种带缓冲区的多层网格lbm演化方法 - Google Patents

Abstract

本发明公开了一种带缓冲区的多层网格LBM演化方法，解决了现有网格加密计算效率低、复杂度高计算量大的弊端，其技术方案要点是在粗网格和细网格的交界处记录两列已经加细的粗网格点和加细后生成的细网格点，形成设定的缓冲区的网格点；对缓冲区的多层网格进行初始化，并分别对粗网格及细网格的网格点进行LBM演化，至多层网格的粗网格和细网格在时间上一致；分别对粗网格及细网格根据设定的插值公式进行空间插值计算，获得粗网格及细网格的分布函数；进行流场信息量的计算，本发明的一种带缓冲区的多层网格LBM演化方法，在不同层次网格的交接处设置缓冲区，降低了计算方法和算法的复杂度，提高了方法的并行性及计算结果精度，缩短了计算时间。

Description

一种带缓冲区的多层网格LBM演化方法

技术领域

本发明涉及流体动力学计算领域，特别涉及一种带缓冲区的多层网格LBM演化方法。

背景技术

格子Boltzmann方法源于格子气自动机方法(LGCA)，而LGCA是更广泛的元胞自动机(CA)在流体力学中的应用。20世纪50年代Von Neumann率先提出在时间和空间都离散的CA数学模型，并将其应用于模拟生物的自复制功能。此后，Broadwell等在提出了用于研究流体激波结构的离散速度CA模型，但时间和空间还是保持连续的。由于CA可以从微观的角度以粗粒度的尺度数值求解非线性偏微分方程描述的宏观物理体，所以为了研究流体的运输性质，法国的Hardy、Pomeau和Pazzis首次提出基于CA的完全离散HPP模型。在该模型中，流体被离散成一系列粒子，且在时间和空间也被离散到规则的正方形格子上。HPP模型满足质量守恒和动量守恒定律，但是由于正方形格子缺乏足够的对称性，HPP模型中的应力张量不能满足各向同性，且存在伪随机守恒量使得HPP模型不能恢复Navier-Stokes方程的非线性项和耗散项，因此HPP模型并没有被广泛应用于流体流动的研究。1986年，法国的Frisch、Pomeau和美国的Hasslacher针对HPP模型存在的问题，提出了具有足够对称性的二维正六边形模型(FHP模型)。同年，法国的d'Humières、Lallemand和Frisch提出了四维面心立方(FCHC)模型及其在三维空间中的投影。FHP模型和FCHC模型都克服了对称性不足的缺点，能分别恢复二维和三维的不可压缩Navier-Stokes方程。因此，LGCA得到了极大关注，并提出了各种更为复杂的LGCA模型，如颜色模型等，这些模型的发展促进了LGCA在各领域的应用。在LGCA中，流体粒子存在于离散的格子节点上，并沿着格线迁移，所有粒子按照一定的碰撞规则同步地相互碰撞与迁移。

1988年McNamara和Zanetti提出了最早的格子Boltzmann模型。该格子Boltzmann模型继承了LGCA的诸多优点，并消除了格子气模型中的大部分统计噪音。但是该模型仍然采用LGCA碰撞方式，存在不满足伽利略不变量、指数复杂性和粒子状态方程与速度有关等缺点。Higuera和Jimenez在1989年对上述模型进一步简化，提出线性化碰撞算子模型，该模型引入平衡态分布函数将碰撞算子线性化，使计算复杂性大大下降，但该方法的数值稳定性较差。同年，Higuera等进一步提出强化算子格子Boltzmann模型，使得平衡态分布函数和碰撞算子的选择不再依赖原有的LGCA模型，而是由所描述的宏观方程确定，碰撞矩阵是一个对称的循环矩阵，且其元素仅仅与离散速度的夹角有关。该模型增强了计算稳定性，并满足伽利略不变性。为了进一步简化碰撞项，Chen和Qian等分别独立提出了单松弛(SingleRelaxation Time，SRT)或BGK(Bhatnagar-Gross-Krook)模型。该模型使用单一松弛系数来代替碰撞项的矩阵模式，用以控制不同方向的粒子靠近各自平衡态的快慢。单松弛模型与矩阵模型类似，但与前面两种模型不同的是当粒子种类数增加时，碰撞算子不会随着粒子数增加而变得复杂。此外，该模型中不再使用Fermi-Dirac分布作为平衡态分布函数，而是使用气体中的Maxwell-Boltzmann平衡态分布函数，克服了格子气模型在推导Navier-Stokes方程对流项时不满足伽利略不变性的问题。该模型可以通过Chapman-Enskog展开恢复到Navier-Stokes方程，且计算效率得到了极大提高。在单松弛模型提出不久，法国学者d'Humières在1992年的一次会议上提出了一种广义格子Boltzmann模型。该模型曾长期被忽视，直到2000年Lallemand和Luo对这一模型做了细致的理论分析，表明该模型在物理原理、参数选取和数值稳定性方面都有很大的优势，该模型也称为多松弛(Multiple-Relaxation-Time，MRT)模型。

在使用标准LBM进行流体计算时，往往采用具有规则的几何形状的计算网格，但在计算中，虽然这类网格简单但缺乏灵活性，很难用于复杂流场。而对于均匀网格，如果需要获得流场中局部区域的精确解，则必须将整个计算区域的网格加密，这将导致所需的硬件资源和计算时间大幅提高，而在流场中经常存在一些区域的物理量变化比较剧烈，空间和时间梯度都较大，所以在这种情况下，经常需要使用局部加细网格的处理方法。

通常的网格加密格子Boltzmann方法往往需要在每个时间步上都进行时间和空间上的两次插值，这将增加计算量，同时会影响整个流场区域的求解精度。

现有的网格加密格子Boltzmann方法，往往采用递归的方式进行计算。导致大量内存的消耗，且调用时传递参数，申请空间，返回时恢复现场，都会增加时间的消耗。降低了整个算法的计算效率，也降低了方法的并行性。

发明内容

本发明的目的是提供一种带缓冲区的多层网格LBM演化方法，通过在不同层次网格的交接处设置缓冲区，解决了网格加密格子玻尔兹曼方法需要时间和空间上进行两次插值的难题，降低了计算方法和算法实现的复杂度，提高了方法的并行性及计算结果的精度，缩短了计算时间。

本发明的上述技术目的是通过以下技术方案得以实现的：

一种带缓冲区的多层网格LBM演化方法，包括有以下步骤：

在粗网格和细网格的交界处记录两列已经加细的粗网格点和加细后生成的细网格点，形成设定的缓冲区的网格点；

对缓冲区的多层网格进行初始化，并分别对粗网格及细网格的网格点进行LBM演化，至缓冲区粗网格和细网格在时间上一致；

分别对粗网格及细网格根据设定的插值公式进行空间插值，计算获得同一时间缓冲区上粗网格及细网格的分布函数；

进行流场信息量的计算。

作为优选，二层网格中，粗网格的分布函数通过中心差分格式获得，具体公式为：

式中，i为二层网格模型的9个离散方向，C_i(i＝0，1，2，…，8)表示粗网格点，R_i(i＝0，1，2，3)表示细网格。

作为优选，二层网格中，细网格R₀，R₁，R₂和R₃的分布函数由粗网格C₀，C₁，C₂，C₃，C₄，C₅，C₆，C₇和C₈这8个点上的分布函数插值获得，具体格式见如下：

i是二层网格模型的9个方向，取值范围为0，1，…，8，

表示从粗网格插值到细网格上的分布函数,

表示粗网格C_n上的f_i分布函数。

作为优选，LBM演化，演化方程如下：

f_i(x+e_iδ_t,t+δ_t)-f_i(x,t)＝Ω_i,

其中f为密度分布函数，x为空间位置，t为时刻，δ_t表示时间间隔，i为离散速度模型的速度方向，τ为格子无量纲松弛时间，Ω_i是碰撞项；

演化方程可分为碰撞及迁移两个过程，分别为：

碰撞过程：f_i ⁺(x,t)＝f_i(x,t)+Ω_i；

迁移过程：f_i(x+e_iδ_t,t+δ_t)＝f_i ⁺(x,t)。

作为优选，粗网格中的松弛时间τ^c和细网格中的松弛时间τ^f需满足如下关系：

其中

表示细网格的格子时间步，

表示粗网格的格子时间步，格子时间步即为演化一次所需要的时间。

为粗网格格子尺度，

为细网格格子尺度；

缓冲区粗网格和细网格具有不同的格子尺度，粗网格上进行一次LBM演化时刻t^c＝1，t^c为粗网格上的时间；细网格上进行2^n-1次LBM演化，时刻值为t^f＝1，保持细网格和粗网格上的时刻值一致。

作为优选，在粗网格和细网格上具体的演化计算如下：

在时刻t粗网格上演化计算：

碰撞：

迁移：

在时刻上细网格进行演化计算：

碰撞：

迁移：

作为优选，三层网格中：

粗网格的分布函数通过中心差分格式获得，具体公式为：

其中，i为三层网格模型的19个方向，

表示从细网格插值到粗网格C₀上的分布函数，

表示细网格R_j的分布函数；

细网格的分布函数由粗网格C_i(i＝0，1，2，…，26)这27个点上的分布函数插值获得，具体格式见如下：

i取值为三层网格模型的19个方向，

表示从粗网格插值到细网格上的分布函数，

表示粗网格C_n上的f_i分布函数。

综上所述，本发明具有以下有益效果：

在多层网格中，粗细网格的交界面采用缓冲区的，避免了使用递归的方法进行计算，提高了方法的并行性，有效的减少在空间上的插值，减少了粗细网格交界面上数据传递次数，缩短了计算时间；

通过设定的高精度插值公式计算出细网格的未知量，同样粗网格上的网格点可以由细网格上的已知量插值得到，既降低了方法的复杂度，也提高了计算结果的精度。

附图说明

图1为交界区域网格示意图；

图2为粗细网格的缓冲区；

图3为缓冲区中粗、细网格的计算情况分析图；

图4为D2Q9模型的二层网格示意图；

图5为D3Q19模型；

图6为D3Q19模型的三层网格示意图；

图7为LBM算法流程图；

图8为圆柱绕流多层网格划分示意图；

图9为缓冲区中粗、细网格的数据传递；

图10为圆柱绕流的涡量图。

具体实施方式

以下结合附图对本发明作进一步详细说明。

本实施例公开的一种带缓冲区的多层网格LBM演化方法，包括有以下步骤：

在粗网格和细网格的交界处记录两列已经加细的粗网格点和加细后生成的细网格点，形成设定的缓冲区的网格点；

对缓冲区的多层网格进行初始化，并分别对粗网格及细网格的网格点进行LBM演化，至缓冲区粗网格和细网格在时间上一致；

分别对粗网格及细网格根据设定的插值公式进行空间插值，计算获得同一时间缓冲区上粗网格及细网格的分布函数；

进行流场信息量的计算。

具体的，采用缓冲区的设置，去除时间上的插值和简化不同尺度网格交界面上的数据传递，便于算法的程序实现。参考图1，我们引入如下的缓冲区设置，见图2。

图2中的黑色实心圆点是粗网格中心点，空心方块点是细网格中心点，而虚线中的网格点即为缓冲区的网格点，空心圆点为缓冲区中粗网格的中心点，黑色实心方块点为缓冲区中细网格的中心点。根据图2中的显示，缓冲区中的网格点即是在粗细网格交界处的地方记录两列已加细的粗网格点和加细后生成的细网格点。

为了保证流体粘性ν和雷诺数Re在整个流场中的一致性，粗网格中的松弛时间τ^c和细网格中的松弛时间τ^f需要满足如下关系：

其中

表示细网格的格子时间步，

表示粗网格的格子时间步，格子时间步即为演化一次所需要的时间。

为粗网格格子尺度，为细网格格子尺度。

根据式(1)可知，粗细网格具有不同的尺度，假定在粗网格上进行一次LBM的碰撞和迁移过程，时刻t^c＝1，t^c为粗网格上的时间，而由式(1)可知，细网格的尺度是粗网格的1/n倍，那么细网格上进行一次LBM的碰撞和迁移的时刻值变为t^f＝1/n，t^f为细网格上的时间，为了保证粗细网格在计算时的时刻一致，那么细网格需要进行2^n-1次LBM的碰撞和迁移过程后，时刻值才为t^f＝1。为方便阐述，这里考虑二层网格，即n＝2的情况，这样粗网格进行一次碰撞迁移过程，细网格需要进行两次碰撞迁移过程才能使得两层网格的时刻值保持一致。

在缓冲区中粗细网格点的LBM计算过程具体如下，假设计算网格如图2所示，那么在网格初始化后，进行第一次的碰撞迁移，进行粗网格点上的演化，缓冲区的粗网格点也进行演化，由于在缓冲区的第二列粗网格点右端没有粗网格点，只有细网格点，那么这一列的迁移过程可不进行，此刻粗网格点上的时刻已变为t+Δt。紧接着对细网格点进行第一次演化，这里缓冲区中细网格点也需要进行演化，此时细网格点的时刻值变为t+Δt/2。类似的，缓冲区中细网格点的第一列不进行迁移，然后再对细网格点进行第二次演化，含缓冲区，这时候缓冲区中细网格点的第二列虽然可以进行迁移，但由于第一列在第一次演化时计算结果已经不正确了，此时第二列的细网格点上的分布函数值也不正确，但其余细网格点上的分布函数值是准确的，这时候细网格点上的时刻也变为了t+Δt。结合图3，可知粗网格中c1列结果准确，c2列结果不准确，细网格中f1和f2列结果不准确，f3和f4列结果准确。

由于目前粗、细网格的时间是一致的，可以通过空间插值计算出粗网格中c2列中的准确数值和f1、f2列中的准确分布函数值。

图1中给出了交界区域网格示意图，黑色实心圆点表示粗网格，黑色实心方块点表示细网格，在虚线框内是交界区域网格。在交界区域网格中，粗网格上的分布函数需要通过插值格式分配到细网格上，同样，细网格上的分布函数也需要通过中心插值格式到粗网格上，下面分别给出二维和三维情况下的粗细网格分布函数间的插值格式。

二维D2Q9离散速度模型的粗细网格插值格式：

根据图4中的D2Q9模型的二层网格示意图，其中左图C_i(i＝0，1，2，…，8)表示粗网格点，R_i(i＝0，1，2，3)表示细网格，右图表示模型的9个离散方向。粗网格C₀的分布函数由细网格R₀，R₁，R₂和R₃通过中心差分格式获得，即

其中i是D2Q9离散速度模型的9个方向，如图4所示，

表示从细网格插值到粗网格上的分布函数，

表示细网格的分布函数。

细网格R₀，R₁，R₂和R₃的分布函数由粗网格C₀，C₁，C₂，C₃，C₄，C₅，C₆，C₇和C₈这8个点上的分布函数插值获得，具体格式见如下：

其中i是D2Q9离散速度模型的9个方向，如图5所示，取值范围为0，1，…，8，

表示从粗网格插值到细网格上的分布函数,

表示粗网格C_n上的f_i分布函数。

三维D3Q19离散速度模型的粗细网格插值格式：

根据图5中的D3Q19模型，其中i是D3Q19离散速度模型的19个方向，取值范围为0，1，…，18，图6中的三维二层网格示意图，粗网格点为C_i(i＝0，1，2，…，26)，粗网格点C₀加密后的细网格点为R_i(i＝0，1，2，…，7)。其中粗网格点C₀的分布函数由细网格R₀，R₁，R₂，R₃，R₄，R₅，R₆，R₇和R₈，通过中心差分格式获得，即

其中i是D3Q19离散速度模型的19个方向，如图5，其中

表示从细网格插值到粗网格C₀上的分布函数，

表示细网格R_j的分布函数。

细网格R₀，R₁，R₂，R₃，R₄，R₅，R₆，R₇和R₈的分布函数由粗网格C_i(i＝0，1，2，…，26)这27个点上的分布函数插值获得，具体格式见如下：

其中，i取值为D3Q19离散速度模型的19个方向，如图5，

表示从粗网格插值到细网格上的分布函数，

表示粗网格C_n上的f_i分布函数。

计算的插值公式如二维情况下的式(3)-式(7)或三维情况下的式(8)-式(16)，具体的计算示意图见图7。这样粗、细网格的全部网格点都在时刻t+Δt，且所有网格点上的分布函数都是准确的，之后可以进行流场信息量的计算，然后再开始下一个时间步的计算。

采用缓冲区的方法能够有效减少在空间上的插值，同时也减少了粗细网格交界面上的数据传递次数。上面的计算步骤给出的是

的情况，其中为粗网格格子尺度，

为细网格格子尺度。实际上对于

的情况以上的计算过程也是可行的。在这种情况下，一个粗网格进行加细后将变成n×n个细网格，在细网格计算的时候需要进行n次碰撞和迁移才能得到t＝1时刻的值，虽然n次碰撞和迁移完成后有f₁,f₂,…,f_n个标记的细网格点上的值不准确，但是可以通过本文创新发明的高精度插值公式计算出细网格的未知量，同样粗网格上标记为c2列的网格点可以由细网格上的已知量插值得到。这样，既降低了方法的复杂度，也提高了计算结果的精度。同时，粗细网格的交界面采用缓冲区的技术，避免了使用递归的方法进行计算，提高了方法的并行性，缩短了计算时间。

为表述清楚，现通过模拟圆柱绕流进行进一步说明：

LBM两种常用的碰撞模型：单松弛模型(LBGK)和多松弛模型(MRT)。这两种模型的LBM演化方程都可以表示为如下形式：

f_i(x+e_iδ_t,t+δ_t)-f_i(x,t)＝Ω_i, (17)

其中f为密度分布函数，x为空间位置，t为时刻，δ_t表示时间间隔，i为离散速度模型的速度方向，τ为格子无量纲松弛时间，Ω_i是碰撞项。

演化方程(17)可以分为以下两部分：

碰撞过程：f_i ⁺(x,t)＝f_i(x,t)+Ω_i, (18)

迁移过程：f_i(x+e_iδ_t,t+δ_t)＝f_i ⁺(x,t). (19)

从以上碰撞(18)和迁移(19)过程可以看出，碰撞过程仅在自身网格点上进行，不涉及其他网格点的信息。而迁移过程中涉及到时间和空间位置的变化，即在计算迁移过程时，需要使用到该网格点上沿着速度方向的相邻网格点信息。不管离散速度模型有多么复杂或者离散速度方向有多少，在迁移过程中的信息交换都仅是局部的，不涉及全局。在迁移计算后，需要更新网格点上的宏观密度、速度，

其中ρ为格点上的密度，u为格点上的速度，e_i为模型的离散速度方向。

从上式可知，宏观量的更新计算也仅在自身的网格点上进行即可，因此LBM方法具有天然的并行性。

针对LBGK和MRT模型，碰撞项的取值如下：

LBGK：

MRT：Ω_i＝[Ω]_i＝{-M^-1·S·[m(x,t)-m^eq(x,t)]}_i. (22)

其中，f_i ^eq为平衡态分布函数，M是变换矩阵，m^eq是矩空间的平衡态函数，S是对角矩阵diag(s₁,s₂,…,s_m)，m(x,t)与f(x,t)的关系如下：

LBM的具体算法流程如图7。

圆柱绕流的实验条件为：入口速度U＝0.1，计算域为-2.5D≤x≤21.5D和-3.5D≤y≤3.5D，其中D是圆柱的直径。雷诺数Re＝DU/ν，D是直径，ν是粘性系数，雷诺数的取值为Re＝100。多层网格划分的示意图见图8，表示对于进行了五层的网格划分。

多层网格LBM算法可以总结为：

1.对粗细网格点上的分布函数f_i ^k,k＝c,f(c表示粗网格，f表示细网格)进行初始化；

2.时刻t粗网格上演化计算：

i.碰撞：

ii.迁移：(缓冲区中与细网格的交界处网格点的分布函数暂不迁移)

3.在时刻

上细网格进行演化计算：

i.碰撞：

ii.迁移：

(缓冲区中与粗网格的交界处网格点的分布函数暂不迁移)；

4.在时刻t+1上，使用空间插值公式(3)-(7)，按照图9中的传递方式对步骤2、3中未进行迁移的网格点上的分布函数进行计算；

5.重复步骤2、3、4，直到满足设置终止条件，例如最大循环次数。

这里

上标"+"表示碰撞后的值。

实验的模拟结果如图10，从图中可以看出卡门涡街。

本具体实施例仅仅是对本发明的解释，其并不是对本发明的限制，本领域技术人员在阅读完本说明书后可以根据需要对本实施例做出没有创造性贡献的修改，但只要在本发明的权利要求范围内都受到专利法的保护。

Claims (7)

1.一种带缓冲区的多层网格LBM演化方法，其特征是，包括有以下步骤：

在粗网格和细网格的交界处记录两列已经加细的粗网格点和加细后生成的细网格点，形成设定的缓冲区的网格点；

对缓冲区的多层网格进行初始化，并分别对粗网格及细网格的网格点进行LBM演化，至多层网格的粗网格和细网格在时间上一致；

分别对粗网格及细网格根据设定的插值公式进行空间插值，计算获得同一时间缓冲区上粗网格及细网格的分布函数；

进行流场信息量的计算。

2.根据权利要求1所述的带缓冲区的多层网格LBM演化方法，其特征是，二层网格中，粗网格的分布函数通过中心差分格式获得，具体公式为：

式中，i为二层网格模型的9个离散方向，C_i(i＝0，1，2，…，8)表示粗网格点，R_i(i＝0，1，2，3)表示细网格。

3.根据权利要求2所述的带缓冲区的多层网格LBM演化方法，其特征是，二层网格中，细网格R₀，R₁，R₂和R₃的分布函数由粗网格C₀，C₁，C₂，C₃，C₄，C₅，C₆，C₇和C₈这8个点上的分布函数插值获得，具体格式见如下：

i是二层网格模型的9个方向，取值范围为0，1，…，8，

(j＝0，1，2，3)表示从粗网格插值到细网格上的分布函数,(n＝0，1，2，…，8)表示粗网格C_n上的f_i分布函数。

4.根据权利要求3所述的带缓冲区的多层网格LBM演化方法，其特征是，LBM演化，演化方程如下：

f_i(x+e_iδ_t,t+δ_t)-f_i(x,t)＝Ω_i,

其中f为密度分布函数，x为空间位置，t为时刻，δ_t表示时间间隔，i为离散速度模型的速度方向，τ为格子无量纲松弛时间，Ω_i是碰撞项；

演化方程可分为碰撞及迁移两个过程，分别为：

碰撞过程：f_i ⁺(x,t)＝f_i(x,t)+Ω_i；

迁移过程：f_i(x+e_iδ_t,t+δ_t)＝f_i ⁺(x,t)。

5.根据权利要求4所述的带缓冲区的多层网格LBM演化方法，其特征是，粗网格中的松弛时间τ^c和细网格中的松弛时间τ^f需满足如下关系：

其中

表示细网格的格子时间步，

表示粗网格的格子时间步，格子时间步即为演化一次所需要的时间；

为粗网格格子尺度，

为细网格格子尺度；

缓冲区粗网格和细网格具有不同的格子尺度，粗网格上进行一次LBM演化时刻t^c＝1，t^c为粗网格上的时间；细网格上进行2^n-1次LBM演化，时刻值为t^f＝1，保持细网格和粗网格上的时刻值一致。

6.根据权利要求5所述的带缓冲区的多层网格LBM演化方法，其特征是，在粗网格和细网格上具体的演化计算如下：

在时刻t粗网格上演化计算：

碰撞：

迁移：

在时刻

上细网格进行演化计算：

碰撞：

迁移：

7.根据权利要求6所述的带缓冲区的多层网格LBM演化方法，其特征是，三层网格中：

粗网格的分布函数通过中心差分格式获得，具体公式为：

其中，i为三层网格模型的19个方向，

表示从细网格插值到粗网格C₀上的分布函数，

(j＝0，1，2，…，8)表示细网格R_j的分布函数；

细网格的分布函数由粗网格C_i(i＝0，1，2，…，26)这27个点上的分布函数插值获得，具体格式见如下：

i取值为三层网格模型的19个方向，

(j＝0，1，2，…，8)表示从粗网格插值到细网格上的分布函数，

(n＝0，1，2，…，26)表示粗网格C_n上的f_i分布函数。

Images