格子Boltzmann方法在GPU平台下对多孔介质流动的模拟

D01:10.19392/j.cnki.l671-7341.201622037

格子Boltzmann方法在GPU平台下对多孔介质

流动的模拟

顾超

北京计算科学研究中心北京100193

：GPU有着强大的并行计算能力，现在越来越多的被人应用到科学计算领域了。格子Boltzmann是一种模拟不可压缩流动的计算方法，

由于其具有天然的并行性，因此本文利用GPU的并行计算能力，在GPU平台上，利用格子Boltzmann模拟周期性的绕方块多孔介质流动，在保证

摘要

程序结果的计算精度的情况下，极大的提高了程序计算效率&同时将计算结果同有限元的结果进行比较，也得到了较好的计算结果。

关键词：GPU;格子Boltzmann;多孔介质；并行计算

20年来，格子Boltzmann方法在计算流体流动的模拟中，起到了很大的作用6但是，与此同时，格子BoltzmamULBM)方法，在计算过程中，有着汁算fi大，计算时间长，内存消耗大的特点。但是由于LBM方

在近

法是一种显示计算格式，在计算的过程中，只耑要利用邻近网格点的信息，因此该方法有着良好的并行性，尤其是在高的并行效率s

GPU平台上，可以达到很

GPU是一个大规模并行计算架构，被广泛的用于图像和非图像计算领域，机器学习，图像，语音等领域都冇着丨、V:用。GPU的主要优势在于单位时间的浮点汁算能力远远的高于CPU的浮点计算能力m。

对于多孔介质流动的数值模拟，多数情况下，多孔介质流动，都属于缓慢流动，而在模拟缓慢流动的数值方法中,LBM方法足其中一种 相当广泛的方法。关于LBM对于多孔介质的流动，在很多领域中都有 着里要的I、V:用:复合材料、流变学、地质学、统计物理、生物科学等诸多 领域都有着广泛的应用3

1

-440M =0

000_ 0

\"

其中式子(1)中f形式所示：

1-1-21-200101

—1-1-122-2-2-2110-101-10201-110-111-20211-11-1000001-1

1

1

1

1

112211

-1-1

11-1-1

一 1

-1

001-1

\\* MERGEFORMAT (4)

t为矩阵QzlT1 x S的分量，这里面矩阵S为如卞

GPU上运行的程序，为

了减少数据从CPU内存上到GPU内存上拷贝的消耗，我们将LBM算 法中的每一步计算都放在GPU上面来实现。

本文的程序是遥于nVIDIA公司CUDA C语言实现的。一共分为三 个部分,第一部分首先介绍了 LBM方法，第二部分介绍了 LBM方法在 GPU平台上的实现。第三部分，我们利用实现的程序对多孔介质流动进

行模拟，并进行相关的计算，将计算得到的结果他文献结果进行比 较，都吻合的很好。

在本文中，我们实现了一个通用的LBM在

S = diag{^sv,sv^sq^,scrsv,sv) \\* MERGEFORMAT(5)

在这个模型中,矩空间{%(x,仏=1，2,…,9}有三个守fr:变物质 的密度r,动M的两个分量，这M动ill为j=ni。

因为LBM模型模拟的是不可压缩多孔介质流动，所以在本文中，

P ~ P〇 +

将会采用如下近似来封傳问题：

P〇 ~\\^P — ^fj

M

9

\\* MERGEFORMAT (6)

1格子Boltzmann方法

j = AU = P〇(M，v), p〇u =

其中，对于矩空间的平衡态分布函数; '的具体形式由如下式子给出：

Boltzmann模咽（MRT),对于MRT模型的方程如下所；qK:

+ + = Qy {mp-mj) + Fi

\\* MERGEFORMAT(l)

本文采用的足二维情况下的9个离散速度方向的多松弛格子

mlq = 8p,mlq = -2Sp +3p0u2

m^1 — Sp — 3p0u2 ,me4q — p^u

这屮K E (是袼子〈上面的一个格点，{cf:i=l, 2,…,9}是介观粒子的 离散速度集合，其中各个方向的粒子速度分别为：

m7 = ~pu->mlq = pv

〇〇\\* MERGEFORMAT (7)

Cj = (0,0), c2 = -c4 = c(l, 0)

c3 =-c5 =c(0,l)

\\* MERGEFORMAT(2)

c6 =-c8 = c(l,l)C7 =-c9 =c(-l,l)其中C 表斯时间步长，时间t为当前时间步；

{/：(^出=1,2,-9}是以速度^运动的单个粒子的分布函数,其中

,2,…,9丨为矩空间分布函数，对于是由如下变换得

me9q - p0uvs _

m? =-p〇v^mlq = a

(^2_v2)

在本文中》对角矩阵S中的两个参数、&分别为如下所示的式子3 V_6v +汾’V

\\* MERGEFORMAT (8)

(8_〜）

那么在本文中，该方法由如下几个步骤来实现：1) 2)

对/进行流动演化；

通过,来计算守恒变ill办和U;

到的。

%=2X

户：i

,/

3) 利用u^u+aG^/2计算uV其中a=F/r〇,这F代表外力；

\\* MERGEFORMAT (3)

4) 利用计算出来的和办计算矩S间的平衡态分布函数， 5 )利用,来计算

6)计算 u**=u*+acfe/2=u+ac^;

7 )将映射到分布函数空间得到分布函数/：;

其中竭为矩阵M的分量，矩阵M的具体形式如下所示s

38 •

科技风2016年11月下8)计算碰撞后的分布函数。

通过该方法实质上是求解如下所所的Navier-Stokes方程，该方程 的形式为如下所示：

5.u + u-Vu =-—p〇

Vz?+

+ F

\\* MERGEFORMAT(9)

V.u = 0

2模拟多孔介质流动的GPU实现在格子

Boltzmann方法中，介观粒子的移动是在如下图所序：的格子

图中的圆点代表每个粒子，从图中可以科出，粒子的分布函数，只 会迁移到周围相邻的几个点〇在实际利用格子

Boltzmann计算中，需要

计算的介观粒子往往很多，几且在模拟多孔介质的流动过程中，往忭需 要迭代很多的时间步，才能得到稳态流动的数值解。因此一个二维问题 的串行程序，忭忭耑要计算儿天甚至一个月的时间才能得到计算结果。

在GPU

上、最基本处理单兀是SP( streaming processor)

,计算过程屮，凡

体的指令和任务都是在SP上进行处理的。GPU进行并行计算,也就是 很多个SP同时做处理。

本文程序的编写，采用的是CUDA平台。在CUDA平台上，最基本 的处理单元是thread,—千CUDA的并行程序会有多个threads来执行。 一定个数的threads会组成一个block,同一个block中的threads可以同 步，也可以通过共享内存(shared memory)通信。同时，多个blocks则会再 构成grid。

具体的示意图如下所示：

图 2 CUDA 平台 Thread,Block, Grid 关系图

电子信息4

在

GPU中，一个GPU所拥有的Thread忭忭有很多个在本文所使

用的Nvidia-K40的显卡中，采用的计算架构是Tesla Kepler,该架构T 面，所能够调用的Thread多达两千万个《

因此，为了提升计算速度，本文在编写多孔介质流动模拟程序的时 候，采用的是一个Thread处理一个格子上面的计算。同时，因为计算过 程屮，：要经常用到格子的离散速度⑷=1,2,…,9},矩阵M,矩阵M的 逆矩阵，所以为了优化程序对这些M的访问速度，本文将这些M放入常 量内存。

同时,在每一步计算过程中，都会用到矩空间的平衡态函数，为了 提高函数的调用效率，本文对于矩空间的平衡态函数采用宏来进行预 定义

由于在计算过程中，会涉及到介观粒子的迁移，所以对于介观粒子 的分布函数，只能存放在GPU的全局内存4{面(global memory)。但足木 文模拟的是多孔介质的流动，该问题的计算区域示意图如下图所示：

图3多孔介质示意图，图(a)为多孔介质，图(b)为该多孔介质的最小体元如上图中，图(a)为整个多孔介质流动的结构，黑色部分为固体，空

白区域为液体。由于多孔介质流动是属于缓慢流动，所以在本文中，可 以将问题的计算区域抽象成如图(b )所示的计算区域。边界条件，对于 固体边界采用的是无滑移边界条件，对于上下左右的流体边界，采用的 是周期性边界条件。

3数值结果3.1数值解的收敛性

从文献[2]可以知道，格子Boltzmann程序的精度是二阶精度，也就 是说，多孔介质的程序计算结果足正确的话，那么该程序的收敛阶必须 是二阶收敛。

下而给出+同间体体积分数的情况下的收敛阶。在本文中s固体体 积分数由如下所示的式子表示：

固体面积\\* MERGEFORMAT(IO)

计算区域的总面积

如下表所示，给出了 a=0.9025,0.01这两种情况下的收敛阶。

表1固体体积分数『0.9025,0.01情况下的收敛阶

0.01

误差

1.656409e-001

6.797929e-0028.561572e-0032.1370

从上面表可以看出，不管a=0.9025是『0.01还是的情况下，计算 出来的收敛阶都接近2,这也就说明了本文所编写的程序收敛阶为2, 所以由此n丨以狞出此多孔介质流动的模拟程序的结果与理论相符合。

3.2渗透率的计算

同时，由于多孔介质流动多为缓慢流动,对于缓慢流动，存在达西 定律，达西定律的公式为如下所示：

U

=—Vp V MERGEFORMAT(l.l)

其中U为多孔介质流动的平均速度,K为多孔介质的渗透

本文将格子

Boltzmann方法计算出来的渗透率同文献[3]中的结果 进行比较，结果图如下所示S

* 39

电子信息科技风2016年11月下

元计算得到的结果也是一样的，与此同时，本文的程序是基于GPU平台 下的大规模几行程序，在模拟同样的问题，同样的网格尺度的情况下， 本文所用的时间是普通中行程序的几十倍,并且随着网格的增多，由于

GPU里而依然是一个thread处塵一个格子，伹是中行程序却足一个线

间却不会有太大的变化。

参考文献：

[1] Nit a C, Itu L M, Suciu C. GPU accelerated blood flow computa­

程处理所有的格子,相比之下，在迭代相同的步数情况下中行程序所消 耗的时间会越来越多，但是本文的GPU平台下的并行程序所消耗的时

图4有限元方法与本文LBM方法的计算结果比较

从上图可以看出，本文所采用的数值方法的GPU程序，4文献I ”所 采用的有限元方法的计算结果吻合的相当的好u 4结论

从前面的章节可以看出，本文基于GPU平台的多孔介质数值模拟， 所得到的结果与理论吻合的很好，，与此同时，计算出来的渗透率同有限

tion using the Lattice Boltzmann Method [CJ.High Performance Extreme Computing Conference (HPEC),2013 IEEE. IEEE,2013:1-6.

[2] He X,Luo L S.Theory of the lattice Boltzmann method:From the Boltzmann equation to the lattice Boltzmann equation|J].Physical Review E, 1997,56(6):68U.

[3] Yazdchi K, Srivastava S,Luding S. Microstructural effects on the permeability of periodic fibrous porous media |J].International Journal of Multiphase Flow，2011,37(8):956_966.

作者简介：

顾超(1991 -)，男，汉族，湖北仙桃，硕士，研究方向：计算流体力学。

(上接第31页）

d 麵

1 口

三、结论

利用Simulink构建了严格遵守客观物理定律的模拟仿真系统，该系 统不受实验环境限制，能够使物理实验以一个比较客观的角度展现出 来，实验结果较准确，与

参考文献：[1]

冯鉴，郭世伟.《基于Simulmk的机械系统可视化建》模仿真分析

FLASH方法模拟实验相比更科学、也更准确地

表达物理量的关系和数值变化。

图6无阻尼时动能势能机械能仿真曲线

[J1.煤矿机械，2002(6):24-6,

[2] http://www.mathworks.com/[3]

吴迪，孙洪毅，等，《基于Matlab Simulink的物理实验-简谐振动

仿真研究》.大连大学，116622.

图7有阻尼时动能势能机械能仿真曲线

40