Boltzmann

Boltzmann（精选7篇）

Boltzmann 篇1

0引言

格子Boltzmann算法, 因其计算简单, 适合处理复杂边界和易于并行, 现已成功地应用于许多流体问题的模拟和建模方面[1]。虽然前人已对格子Boltzmann算法的并行性能进行了一些研究[2,3,4,5,6], 但由于测试性能时基于的流体问题不同, 计算机系统不同 (硬件体系结构、编程实现方式等) , 所得到的结果 (加速比、效率) 也有很大的差异。并且缺少定量化的理论分析, 其结果往往依赖于硬件环境, 不利于指导今后的大规模数值计算。

本文基于顶盖驱动的方腔流[7,8,9], 对格子Boltzmann算法的并行性能作了全面系统的测试与分析, 这为以后编写高效的并行程序提供了理论和实践依据, 能够较大地提高程序的效率, 节省计算时间, 充分利用并行计算机的计算资源。

1格子Boltzmann算法

格子Boltzmann算法是一种模拟流体流动的数值方法, 它不再基于连续介质假设, 而是把流体看成许多只有质量没有体积的微粒组成, 这些微粒可以向空间的若干方向任意流动。它的主要思想就是以简单规则的微观粒子运动代替复杂多变的宏观现象。粒子在每个时间步的运动由两个子步构成, 即:

(1) 碰撞 当多个粒子到达同一网格结点时, 它们按碰撞规则相互作用并改变各自的速度方向;

(2) 流动 每个粒子按其速度方向移动到最近的网格结点。

文献[10]中提出格子Boltzmann算法的单松弛模型, 即:LBGK (Lattice Bhatnagar-Gross-Krook) 模型。目前, LBGK模型是格子Boltzmann算法中最主要的, 也是应用最广泛的模型。LBGK模型抛弃了Fermi-Dirac分布, 采用Maxwell分布, 使得各向同性、Galilean不变性和压力独立于流速等问题都得到满足。Qian等将格子Boltzmann方程中的碰撞项完全线性化, 用单松弛时间逼近碰撞项, 其对应的格子Boltzmann演化方程为:

fα (x + eα, t + 1) = fα (x, t) + ω (fαeq (x, t) -fα (x, t) ) (1)

其中, ω为松弛因子, 为了达到良好的稳定性, 必须在0到2之间, 小于1为亚松弛, 大于1为超松弛; fα为单粒子分布函数, 表示在t时间上, x格点上沿运动有1个粒子的概率; eα称为粒子运动方向或离散速度方向;fαeq为平衡态分布函数, 是流体动力学量 (如密度、动量等) 的函数。

构造LBGK模型的关键是选择恰当的平衡态分布函数feqα。一旦选定了离散速度eα, 只须确定出权系数tp就可得到平衡态分布函数。DnQb系列模型的分布函数为:

$f{}_{\alpha }^{eq}(x,t)=t{}_p\rho \left\{1+\frac{e{}_{\alpha }\cdot u}{c{}_s^2}+\frac{(e{}_{\alpha }\cdot u){}^2}{2c{}_s^4}-\frac{u{}^2}{2c{}_s^2}\right\}$ (2)

其中cs为声速;tp对应粒子速度的权系数;u为宏观速度;ρ为宏观密度。

对于本文所选用的D2Q9和D3Q19模型, 权系数选取如下:

D2Q9 :$t{}_0=\frac{4}{9}t{}_1=\frac{1}{9}t{}_2=\frac{1}{36}c{}_s^2=\frac{1}{3}$

D3Q19:$t{}_0=\frac{12}{36}t{}_1=\frac{2}{36}t{}_2=\frac{1}{36}c{}_s^2=\frac{1}{3}$

2格子Boltzmann算法的并行化

对于格子Boltzmann算法, 首先就要将整个计算分解成一些小的任务, 尽量开拓并行执行的机会。计算网格的规则划分方式通常有三种, 分别是一维带状 (Slice) 、二维块状 (Block) 或棋盘 (Checkerboard) 、三维盒状 (Box) 划分, 如图1所示。

对于不同域分解方式执行效率的比较, 可以从各个方向上的网格点数的接近程度、各子域的通信量以及总通信量进行分析。

负载均衡对于流体计算来说, 就是每个计算子域包含的网格数尽量接近[11]。假设要求解问题的维度为dim (一般, dim≤3) , np为处理器数, 每维网格数size[i]的权重为:

$weights[i]=\frac{size[i]}{\underset{0\le i\le dim-1}{{\displaystyle \max }}size[i]}$ (3)

划分后每个方向上的处理器数p[i], 为了使得划分后各个方向上的网格点尽量接近, 即在约束条件:

${\displaystyle \prod _{0\le i\le dim-1}p}[i]=np$ (4)

下求:

$\min {\displaystyle \sum _{0\le i\le dim-1}p}[i]\times weights[i]$ (5)

特别地, 当dim=3时,

${\displaystyle \sum _{0\le i\le 2}p}[i]\times weights[i]\ge 3\sqrt[3]{np}\times weights[0]\times weights[1]\times weights[2](6)$

如果size[0]=size[1]=size[2], 即:

weights[0]=weights[1]=weights[2]=1

则当p[0]=p[1]=p[2]=np1/3, 式 (6) 取得最小值。

流场的划分方式同样也决定了通信量的大小[12]。三维流场可以采用三种方式来进行划分, 下面以三维方腔为例来介绍流场在不同划分方式下的通信量。假设流场的尺寸为M×N×P。

(1) 一维slice划分 将三维流场沿一个方向上划分成若干立方块。

设在x方向划分为np块, 则x方向上每块的长度为M/np, 每块的计算量为总计算量的1/np。对于非边界子块, 它需要与左右两个方向的相邻子块交换的共享数据, 其通信量为2N×P, 总的通信量为 (np-1) ×N×P。在y、z方向上的划分与x方向类似。对于不同的M、N、P的值, 首先在最大的维数上分割所产生的通信量最小, 而当M=N=P时, 在各个方向上划分无异。

(2) 二维block划分 将三维流场沿两个方向划分成若干立方块。

设在x方向划分为px块, y方向划分为py块, 每块的计算量为总计算量的1/ (px×py) 。对于非边界子块, 它需要与周围四个方向上的相邻子块交换边界数据, 其通信量为2 (M/px+N/py) ×P, 总的通信量为 (px-1) × (N×P/py) + (py-1) × (M×P/px) 。

(3) 三维box划分 将三维流场沿三个方向划分成若干立方块。

设在x、y、z各个方向上分别划分成px、py、pz块, 每块的计算量为总计算量的1/ (px×py×pz) 。对于非边界子块, 它需要与六个近邻交换边界数据, 其通信量为2[ (M/px) × (P/pz) + (M/px) × (N/px) + (N/px) × (P/pz) ], 总的通信量为:

(px-1) × (N×P) / (py×pz) + (py-1) × (M×P) / (px×pz) + (pz-1) × (M×N) / (px×py)

因此, 二维和三维划分在工作负载尽量均衡的条件下, 单个子块的通信量以及总的通信量均达到最优。

根据以上理论分析的结果, 考虑三维方腔流流场 (M=N=P≡W) , 在对流场进行一维、二维、三维划分时, 子域的计算量均为W3/np, 其中二维划分时pxblock×pyblock=np, 三维划分时pxbox×pybox ×pzbox =np。子域之间需要交换的信息量分别为2W2、2 (pxblock+ pyblock) W2/np、2 (pxbox+pybox +pzbox) W2/np。若处理器平均计算一个网格点数据所需时间为τc, 子域之间平均传递一个数据所需时间τm, 记$\kappa =\frac{\tau {}_c}{\tau {}_m}$。那么, 一维、二维、三维划分的计算通信比分别为:

$\beta {}_{slice}=\frac{W{}^3/np\times \tau {}_c}{2W{}^2\times \tau {}_m}=\frac{1}{2}\kappa \frac{W}{np}$ (7)

$\begin{array}{l}\beta {}_{block}=\frac{W{}^3/np\times \tau {}_c}{2(px{}_{block}+py{}_{block})W{}^2/np\times \tau {}_m}\\ =\frac{1}{2}\kappa \frac{W}{px{}_{block}+py{}_{block}}(8)\end{array}$

$\begin{array}{l}\beta {}_{box}=\frac{W{}^3/np\times \tau {}_c}{2(px{}_{box}+py{}_{box}+pz{}_{box})W{}^2/np\times \tau {}_m}\\ =\frac{1}{2}\kappa \frac{W}{px{}_{box}+py{}_{box}+pz{}_{box}}(9)\end{array}$

根据上述公式, 三种划分方式计算通信比之间的关系为:

$\beta {}_{block}=\frac{np}{px{}_{block}+py{}_{block}}\beta {}_{slice}$ (10)

$\beta {}_{box}=\frac{np}{px{}_{box}+py{}_{box}+pz{}_{box}}\beta {}_{slice}$ (11)

$\beta {}_{box}=\frac{px{}_{block}+py{}_{block}}{px{}_{box}+py{}_{box}+pz{}_{box}}\beta {}_{block}$ (12)

3性能分析

本文所有的测试都是使用GNU C调用MPI来实现的。在本节中, 通过分析实验的结果来验证理论分析得出的结论, 并和文献中的结果进行比较。

3.1不同域分解方式下的并行性能比较

(1) 各个方向上网格点数对并行性能的影响

以三维方腔流问题为例, 采用三维Box划分方法对计算区域进行分解, 选择问题规模为257×257×257, 模型为D3Q19, 迭代步数为1000步, 捆绑阻塞型发送和接收, 在至强3000上使用64个处理器进行测试。

图2中, 同等问题规模, 在相同数量的处理器上运行, 三维划分方式下处理器分配方式的不同导致各子区域网格点数不同。墙上时间随着各个方向上的网格点数的逐渐接近而逐步下降, 当各个方向上的网格点数相等, 当处理器划分为4×4×4, ${\displaystyle \sum _{i}p}[i]\times weights[i]=4\times 1+4\times 1+4\times 1=12$最小, 各个方向上的网格点数均为257/4, 即负载尽量均衡时, 执行速度最快。

(2) 三种划分方式对并行性能的影响

Satofuka和Nishioka在1999年对格子Boltzmann算法在二维问题中不同划分方式下 (水平划分、垂直划分和块划分) 的并行性能进行了分析, 采用Fortran90在拟向量机Hitachi SR2201上进行测试, 得出的结论是:垂直划分 (按列) 优于按行及二维块划分。这里我们将问题扩展到三维, 对格子Boltzmann算法在三种划分方式下在负载均衡时的并行加速比进行了分析 (如图3所示) 。

测试条件:流场规模129×129×129, 处理器数np<32, 捆绑阻塞式通信方式。

在处理器数np≤16三种划分方式的加速比差别不大, 尤其是二维和三维划分方式。公式 (7) - (9) 对三种划分方式计算通信比的理论分析在一定程度上能描述总体趋势, 而在实际的测试过程中, 并行加速比受多方面因素的影响, 如在程序运行过程中处理器和磁盘资源不是独占的, 在处理器多于8之后, 二维和三维划分所受影响要比一维划分的大, 从而导致测试结果出现波动。但随着处理器数的增多三维划分比其它两种划分的加速比要高, 并有持续升高的趋势, 这说明流场的多维划分更适合大规模处理器上的并行化计算。

(3) 正方形 (体) 流场下各个方向上划分的块数的排列对并行性能影响

不管采用哪种划分方式, 在满足各个方向上的网格点数尽量接近的条件下, 对于正方形 (体) 网格来说, 处理器个数的排列如二维划分px×py和py×px对性能基本没有影响。C语言是行优先方式存储数组, 如图4所示, 二维数组划分成2×2的块, 假设将其映射到2个处理器上。左图为二维数组分块后在水平划分方式“H”、垂直划分方式“V”下的处理器映射关系;右图是分块后各块在内存空间中的地址排列。

在共享存储模型下, 垂直划分后的数据在地址空间上不连续, 增加了访存开销, 同时在与相邻子域交换边界时需对不连续数据进行特殊处理;而在分布式存储的集群环境下, 数据在局部存储器中连续存储, 垂直划分相对水平划分不会增加额外的访存开销。

本文以二维正方形方腔1025×1025、D2Q9模型、捆绑阻塞通信为例, 比较了一维划分时水平划分和垂直划分在提高性能上的差别。从图5中可以看出, 水平划分和垂直划分的加速比极为接近, 这是因为在问题规模相同, 处理器数也相同时, 各个子域的网格点数均为10252/np, 计算通信比也相同。同时, 分布式存储环境下, 也消除了C语言行优先方式导致的访存开销上的差异。

3.2不同通信方式下的并行性能比较

在划分方式选定之后, 各子域之间产生的数据交换采用什么通信方式, 对程序的并行性能也会产生一定的影响。MPI中最常用的数据传送方式是标准的阻塞 (Blocking) 和非阻塞 (NonBlocking) 两种机制[13]。非阻塞通信方式在某些并行机上可以实现计算与通信重叠, 文献[5]均对格子Boltzmann算法在不同通信方式下的并行加速比作了分析和比较。这里, 以二维方腔流问题为例, 网格规模为1025×1025, D2Q9模型, 迭代步数1000, 从图6中的在至强3000上的实验结果来看, 阻塞通信 (捆绑式MPI_Sendrecv) 和非阻塞通信 (MPI_Isend/recv) 对格子Boltzmann算法在提高算法性能上基本没有什么差别, 原因在于格子Boltzmann算法是计算密集性算法, 通信时间只占用总执行时间的很少一部分, 实验的结果与文献[5]中的结论相吻合。

3.3算法的可扩放性分析

(1) 调整问题规模及处理器数时分析算法的性能

首先对二维流动问题进行分析, 增加问题规模从1292到10252;增加处理器数从2到32个处理器。如图7所示。

固定问题规模, 加速比随处理器数的增加基本呈线性的增长, 增加处理器数虽然提高了执行速度, 但随之而来的通信时间所占的比例也会大幅增加, 计算通信比降低。在129×129的网格上测试时, 当处理器数达到18以后, 加速比的变化幅度降低, 基本维持恒定。Amdahl加速比定律也说明了这一点, 加速比在处理器数增大时会有一个上限。

在增加处理器数量同时增大问题规模的情况下, 加速比会随着的处理器数量的增加而增大。

固定处理器数, 随着问题规模的增加, 加速比呈上升趋势。

(2) 不同体系结构集群上的算法性能

可扩放性是算法和体系结构的结合体。考察格子Boltzmann算法, 在不同体系结构的集群环境 (IWill 刀片服务器集群环境及至强3000 SMP高性能集群) 下的性能, 发现在格子Boltzmann算法同样具有良好的并行性能。图8给出了在IWill 刀片服务器集群环境下的格子Boltzmann算法的加速比。问题规模10252, 捆绑阻塞式通信, 迭代1000步。

4结论

本文以方腔流为例对格子Boltzmann算法的并行性能进行了系统分析, 并从理论上和实验数据上系统地进行了分析对照, 得出以下结论:

(1) 格子Boltzmann算法在不同体系结构的并行机上, 均具有良好的并行性及可扩放性;

(2) 在域分解时, 各个方向上的网格点数尽量接近, 可优化通信时间, 从而加快执行时间;

(3) 高维区域划分技术适合于大规模处理器上的并行计算, 在处理器个数较少时Block划分和Box划分效果相当;

(4) 负载均衡时, 划分时处理器排列方式不会影响格子Boltzmann算法并行性能;

(5) 至强3000, 非阻塞和阻塞两种通信方式下均能获得良好的加速比。

最后, 感谢上海高校网格技术E-研究院提供的实验环境。

Boltzmann 篇2

扩散过程在物理、化学、生物等许多领域中起着很重要的作用。热扩散方程是描述传热过程的一个重要方程, 但在复杂的边界和初始条件下, 解析求解是很困难的。许多学者利用格子Boltzmann方法求解扩散问题, 并取得了很多成果。刘慕仁等人给出了求解一维有源扩散方程的格子Boltzmann模型, 确定了局部平衡函数Chapman-Enskog展开的待定系数[5], 他们还利用格子Boltzmann方法求解了一维对流扩散方程, 确定了方法中的粘滞系数与对流系数的关系[6]。徐世英等人利用浓度分布的Chapman-Enskog展开及多尺度技术, 得出了Tyson反应扩散的反应物和催化剂随时间的浓度空间分布值[7]。本文在前人的研究基础上, 提出了格子Boltzmann方法求解一维热扩散方程的有效数值方法。

1 热扩散方程的Boltzmann模型

一维热扩散方程表达式

式中T———温度, 要求T≥0;

u———速度;

λ/ρcp———热扩散系数;

λ———导热系数;

ρ———密度;

cp———定压比热容。

1.1 平衡态分布函数

将一维空间离散成均匀的线性格子, 每个节点与相邻的节点相连。将速度ca=cea离散为三个方向, ea可取值-1, 0, 1, 分别对应粒子速度方向向左、不动和向右三种情况。c=Δx/Δt为粒子迁移速率, 3-速格子模型如图1所示。

用fi (x, t) 表示沿ea方向的粒子分布函数, 表示在时间t时, 在位置x处粒子出现的概率, 应有fi≥0, 宏观量T定义为

式中fieq———平衡态分布函数。

根据统计物理, 分布函数遵守以下格子Boltzmann方程

式中τ———弛豫时间。

式 (3) 也称LBGK (Lattice BGK, 简称LBGK) 方程, 为满足稳定性条件, 要求τ≥0.5。

由式 (3) 可知, 只要得到feq (x, t) 就可以得出以后时间的分布函数f (x, t) 。由LBGK方程恢复宏观方程的关键是选择适当的平衡态分布函数。由统计力学可知, 在局域平衡条件下, 分布函数可表示为局域平衡量的函数, 考虑到ea的对称性, 将平衡态分布函数表示为

式中a0、a1、a2———待定系数。

1.2 格子Boltzmann方程

对式 (3) 的左侧项进行空间和演化时间的Taylor级数展开, 取其到二级项可得

应用Chapman-Enskog方法对时间系数、空间系数和分布函数进行多尺度展开, 可得到

式中ε———小参量。

将式 (6) 代入式 (5) , 在ε一阶小量下有

在ε二阶小量下有

将式 (7) 和式 (8) 分别对下标求和得出

由此得到的式 (9) 、式 (10) 即为时间t1和t2尺度下的密度平衡方程。

1.3 宏观热扩散方程

将式 (4) 代入式 (2) , 可以得到

为了在t1尺度上恢复方程的对流项, 令

由式 (9) 、式 (12) 得到

由式 (7) 和式 (14) 得到

将式 (15) 代入式 (10) , 得到

将式 (14) 乘以ε, 式 (16) 乘以ε2, 两式相加便可回到宏观尺度下的方程 (1) 。

令, 则

由式 (11) 、式 (13) 和式 (17) 可得出待定系数

2 数值实验

为验证得到的热扩散方程的Boltzmann模型, 采用文献[6, 8]中的算例, 进行了验证性模拟实验。

2.1 验证算例1

算例1中热扩散方程的边界条件为第一类边界条件且为恒温, 具体的控制方程

算例1的解析解在文献[6]中已给出

算例1中热扩散方程的格子Boltzmann模型求解参数:τ=1.2, u=0.2, v=0.5, 时间步长Δt=0.01, 空间步长Δx=0.1, 粒子的迁移速率为c=10, 总长度l=40, t分别为280、290、300和310。

图2给出了不同时刻t下算例1的模拟值与解析解。从图中可以看出, t为280、290、300和310时, 模拟值与解析解的最大绝对误差分别为1.31E-04、8.70E-05、7.20E-05和7.00E-06, 从而说明数值解与解析解吻合良好。

网格数目对模拟结果精度有一定的影响。本文分析网格数目为80、200和400时, 对计算精度的影响情况。

表1为τ=1.2, u=0.2, v=0.5, 时间步长Δt=0.01, t=300, 总长度l=40, 空间步长Δx分别为0.5, 0.2和0.1时, 各个节点在不同网格数目下的解析解与数值解及其绝对误差。从表中可以看出, 网格数目分别为80、200和400时的最大绝对误差分别为3.02E-03、3.31E-04和7.00E-05。由此可见, 随着网格的细化, 模拟结果的绝对误差越来越小。

2.2 验证算例2

算例2中的热方程的边界条件是随着时间t变化的, 具体方程如下

该问题的解析解在文献[8]中已给出

图3给出了不同时刻t下算例2的模拟值与解析解。从图中可以看出, t分别为1、4、7和10时的最大绝对误差分别为6.21E-04、4.62E-04、4.98E-04和5.09E-04, 从而说明数值解与解析解吻合良好。

3 结语

本文基于格子Boltzmann方法, 采用ChapmanEnskog展开和多尺度技术, 建立了求解一维热扩散方程的D1Q3模型, 模型的数值解与文献的解析解吻合良好, 其两者的误差随网格细化而大幅度减小, 从而说明了本文模型可用于求解一维热扩散方程。

参考文献

[1]Raabe D.Overview of the lattice Boltzmann method for nano-and microscale fluid dynamics in materials science and engineering[J].Modeling and Simulation in Materials Science Engineering, 2004, 12 (6) :R13-R46.

[2]Tang G H, Tao W Q, He Y L.Gas slippage effect on microscale porous flow using the lattice Boltzmann method[J].Physical Review E, 2005, 72 (5) :056301.

[3]Grunau D, Chen S, Eggert K.A lattice Boltzmann model for multiphase fluid flows[J].Physics of Fluids, 1993, 5 (10) :2557-2562.

[4]Miller W, Succi S, Mansutti D.Lattice Boltzmann model for anisotropic liquid-solid phase transition[J].Physical review letters, 2001, 86 (16) :3578-3581.

[5]邓敏艺, 刘慕仁, 李珏, 等.一维有源扩散方程的格子Boltzmann方法[J].广西师范大学学报:自然科学版, 2000, 18 (1) :9-12.

[6]刘慕仁, 何云, 陈若航, 等.一维对流扩散方程的格子Boltzmann方法[J].广西科学, 1999, 6 (3) :168-169.

[7]徐世英, 卫玉敏, 吴春光, 等.一维Tyson反应扩散系统的格子Boltzmann方法模拟[J].计算机与应用化学.2008, 25 (5) :537-540.

Boltzmann 篇3

顶盖驱动流场的数值模拟是流体力学中的一个经典问题, 其数值分析和实验结果都取得了一些成果, 本文采用BGK碰撞模型, D2Q9格子划分来对顶盖驱动流进行数值模拟计算。

1 二维问题LB模型

Boltzmann方程是统计力学中用来描述非平衡态分布函数演化规律的方程, 其方程式如下[3]:

$\frac{\partial f}{\partial t}+\begin{array}{l}\to \\ u\end{array}\cdot \nabla {}_{\begin{array}{l}\to \\ x\end{array}}f+\frac{\begin{array}{l}\to \\ F\end{array}}{m}\cdot \nabla {}_{\begin{array}{l}\to \\ u\end{array}}f=\text{Ω}(f)(1)$

其中, $f(\begin{array}{l}\to \\ x\end{array}‚\begin{array}{l}\to \\ u\end{array}‚t)$表示t时刻, 空间$\begin{array}{l}\to \\ x\end{array}$处, dxdydz体积内, 速度在$\begin{array}{l}\to \\ u\end{array}+\begin{array}{l}\\ d\end{array}\begin{array}{l}\to \\ u\end{array}$范围内的平均分子数;$\begin{array}{l}\to \\ F\end{array}$表示外力, m为分子量, Ω (f) 为碰撞函数。对于计算繁杂的碰撞函数, 采用BGK碰撞模型来近似如下:

$\text{Ω}(f)=-\frac{1}{\tau }(f-f{}^{eq})(2)$

式中:τ——无量纲松弛时间;

feq——平衡态分布函数;

对于D2Q9模型, 如图1所示, 格点上的速度为[4]:

$\begin{array}{l}\to \\ e{}_a\end{array}=\{\begin{array}{l}0,a=0\\ \left(\cos \frac{a-1}{2}\pi ‚\sin \frac{a-1}{2}\pi \right)‚a=1‚2‚3‚4\\ \sqrt{2}\left(\cos \left(\frac{a-5}{2}\pi +\frac{\pi }{4}\right)‚\sin \left(\frac{a-5}{2}\pi +\frac{\pi }{4}\right)\right)‚a=5‚6‚7‚8\end{array}(3)$

Boltzmann方程的离散化速度模型为:

$f{}_a(\begin{array}{l}\to \\ x\end{array}+\begin{array}{l}\to \\ e{}_a\end{array}\Delta t‚t+\Delta t)$

$=f{}_a(\begin{array}{l}\to \\ x\end{array}‚t)-\frac{f{}_a(\begin{array}{l}\to \\ x\end{array},t)-f{}_a^{eq}(\begin{array}{l}\to \\ x\end{array},t)}{\tau }a=0,1,\cdots ‚8(4)$

其中fa和f${}_a^{eq}$分别为与第a个方向上的速度$\begin{array}{l}\to \\ e{}_a\end{array}$相对应的粒子速度分布函数和平衡态分布函数。

$\begin{array}{l}f{}_a^{eq}(\begin{array}{l}\to \\ x\end{array})=w{}_a\rho (\begin{array}{l}\to \\ x\end{array})\left[1+3\frac{\begin{array}{l}\to \\ e{}_a\end{array}\cdot \begin{array}{l}\to \\ u\end{array}}{c{}^2}+\frac{9}{2}\frac{(\begin{array}{l}\to \\ e{}_a\end{array}\cdot \begin{array}{l}\to \\ u\end{array}){}^3}{c{}^4}-\frac{3}{2}\frac{\begin{array}{l}\to \\ u\end{array}\cdot \begin{array}{l}\to \\ u\end{array}}{c{}^2}\right]\\ a=0,1\cdots ‚8(5)\end{array}$

式 (5) 中wa为粒子分布函数的加权系数, 其在D2Q9模型上取值如下[4,5,6]:

$w{}_a=\{\begin{array}{l}4/9‚a‚0\\ 1/9‚a=1\sim 4\\ 1/36‚a=5\sim 8\end{array}(6)$

方程的计算可以分为两步:

(1) 碰撞:$f{}_a^{*}(\begin{array}{l}\to \\ x\end{array},t)=f{}_a(\begin{array}{l}\to \\ x\end{array},t)-\frac{f{}_a(\begin{array}{l}\to \\ x\end{array},t)-f{}_a^{eq}(\begin{array}{l}\to \\ x\end{array},t)}{\tau }$

(2) 迁移:$f{}_a(\begin{array}{l}\to \\ x\end{array}+\begin{array}{l}\to \\ e{}_a\end{array}\Delta t,t+\Delta t)=f{}_a^{*}(\begin{array}{l}\to \\ x\end{array},t)$

其中宏观量与微观量的关系式如下:

$\rho ={\displaystyle \sum _{a=0}^{8}f}{}_a(7)$

$\begin{array}{l}\to \\ u\end{array}=\frac{1}{\rho }{\displaystyle \sum _{a=0}^{8}f}{}_a\begin{array}{l}\to \\ e{}_a\end{array}(8)$

$v=\frac{2\tau -1}{6}\frac{\delta x}{\delta t}(9)$

2 数值计算步骤

对于给定的问题, 首先要确定特征尺度, 将物理量无量纲化, 然后根据Re数, 以及确定的格子常数 (格子空间长度和时间步长) 计算无量纲运动粘度, 再由式计算出松弛时间τ。开始迭代计算时, 首先计算宏观变量$\rho ‚\begin{array}{l}\to \\ u\end{array}$, 其次计算顶盖边界各个格点的微观变量fa, 再次各点进行碰撞计算, 对其余三个边界采用bounceback边界处理, 最后对各格点进行迁移计算, 重复该过程最后得到整个流场的稳定解。

其顶部运动边界如下:

$\{\begin{array}{l}\rho =\frac{f{}_0+f{}_1+f{}_3+2(f{}_2+f{}_5+f{}_6)}{1+v{}_0}\\ f{}_4=f{}_2-\frac{2}{3}\rho v{}_0\\ f{}_7=f{}_5+\frac{1}{2}(f{}_1-f{}_3)-\frac{1}{2}\rho u{}_0-\frac{1}{6}\rho v{}_0\\ f{}_8=f{}_6-\frac{1}{2}(f{}_1-f{}_3)+\frac{1}{2}\rho u{}_0-\frac{1}{6}\rho v{}_0\end{array}(10)$

3 流场数值模拟

本文对顶盖驱动流采用D2Q9模型进行模拟, 初始条件设置为最顶部节点无量纲速度为u=0.1, v=0, 其余各节点速度为0, 网格划分为128×128。顶部边界无量纲边界。

各种雷诺数条件下方腔内计算结果的流线如图1所示。

各种雷诺数条件下方腔内计算结果的等压线如图2所示。

4 结 论

本文采用的D2Q9模型处理方腔内顶盖驱动流场数值模拟, 对各种雷诺数条件下的流场进行了模拟, 从计算结果可以看出, 在低雷诺数, 和较高雷诺数的条件下, Lattice Boltzmann方法都能对方腔内, 顶盖驱动流场作出准确的模拟, 但在高雷诺数下, 数值计算需要较多时间步骤才能收敛。

LBM方法是与时间过程相关, 非常适合处理非稳定流动问题。同时其又是基于分子统计力学方程, 因此更能从介观层面上去描述物理现象, 对于化工领域中许多复杂的流动传热现象能够做出较为精准的模拟。LBM方法边界条件易于实施, 因此对于处理复杂几何形状内的流动以及处理多孔介质内的流动与传热问题有很大的优势。作为一种新的模拟和建模方法, 其理论思想简单, 编程容易, 应用非常广泛, 逐渐成为一种有效的流场模拟手段, 甚至弥补传统计算方法的不足, 比如传统计算方法难以模拟的多孔介质问题, 而Lattice Bohzmann方法却能从容处理。

目前Lattice Bohzmann方法还处于不断发展之中, 作为一种新生事物, 人们对它还缺乏了解和比较陌生, 需要大量的论证和实践, 挖掘其潜力, 开拓新的应用领域。而且这种打破传统的数学建模观念, 自底向上建模, 不需建立和求解复杂的偏微分方程, 在针对化工领域中复杂传热流动及相态变化的模拟中, 具有广阔应用前景的一种高效建模数值算法。

摘要：Lattice Boltzmann方法是近十年来发展起来的基于介观层次上的新算法, 是一个崭新的求解流体系统偏微分方程的行之有效且效率比较高的数值计算方法, 以其为物理现象建模的手段, 在很多领域得到了应用, 尤其在化工领域中具有复杂边界和复杂流场的模拟中取得巨大成功, 如湍流、多孔介质流、化学反应流、燃烧模拟、结晶过程、磁流体等等。采用LB方法中的D2Q9模型, 采用BGK碰撞模型及BounceBack边界条件数值模拟流体在顶盖驱动中的流动现象, 给出流场的流线和等压线图, 计算结果表明该方法运用到顶盖驱动流动计算中, 所得结果与基准解吻合良好。

关键词：Lattice Boltzmann方法,数值模拟,顶盖驱动,BGK碰撞模型,D2Q9模型

参考文献

[1]李元香, 康立山.格子气自动机[M].南宁:广西出版社, 1994.

[2]G M, G Z.Use of the Boltzmann Equation to Simulate Lattice Gas Au-tomata[J].Physical Review Letters, 1988, 61:2332-2335.

[3]Chen H, Chen S, Matthaeus W H.Lattice Boltzmann Model for Simu-lation of Magnethydrodynamics.Phs.Rev.A, 1992, 45:5339.

[4]谭玲燕, 陈明杰, 施卫平, 等.用Lattice Boltzmann方法计算矩形柱的绕流问题[J].吉林大学学报, 2009, 47 (2) :273-276.

[5]郭照立, 郑楚光.格子Boltzmann方法的原理及应用[M].北京:科学出版社, 2008.

Boltzmann 篇4

格子Boltzmann方法 (LBM) [1,2]是近几十年发展起来的一种新型的数值方法, 它不同于传统数值方法的流体计算和建模方法。它的微观粒子背景使得它具有许多其他数值方法所没有的独特优点。同时格子Boltzmann方法 (尤其是LBGK模型) 的演化过程非常简单清晰, 程序比较简洁。所以从格子Boltzmann方法诞生之日起就受到国内外专家的广泛关注。

格子Boltzmann方法除了能够用于等温流动外, 还可以用于有温度变化的热力学问题。包括两种模型:多速度模型和双分布函数模型[3]。其中的双分布函数模型由于数值稳定性较好, 在应用中具有一定的优势。双分布函数模型的基本思想是, 如果粘性热耗散和压力所作的功可以忽略, 温度就可以看作是一个跟随流体运动的被动量, 并且满足一个简单的对流扩散方程, 因此可以使用一个独立于密度分布的新分布函数来模拟温度场。所以双分布模型使用两类分布函数:密度分布函数和温度分布函数。其中密度分布函数用于模拟速度场, 而温度分布函数则用来模拟温度场[4]。

采用这种新的热格子Boltzmann模型模拟封闭方腔中的自然对流方式对在不同Ra值下的流动涡运动和温度扩散进行了比较。

1 物理模型

图1所示为二维封闭方腔的物理模型, 方腔的左壁和右壁温度恒定, 且左壁的温度较高, 右壁的温度较低, 而上下两壁是绝热的。坐标系的原点取为方腔的左下角, 并以水平方向为x方向, 重力方向为y方向。流动的初始条件为:u=v=0, T=T0;边界条件是:

左壁:u=v=0, T=Th

右壁:u=v=0, T=Tc

上壁和下壁:u=v=0, ∂T/∂y=0

其中, u和v分别为速度在x方向和y方向上的分量;Th>Tc, T0= (Th+Tc) /2。

自然对流的两个最基本的无量纲参数是Rayleigh和Prandtl数, 分别定义为:

Ra=|g|βΔTH3/ (vD) , Pr=v/D

式中:ΔT是高温壁和低温壁的温度差;H为方腔的高度;β为体积膨胀系数;v为运动粘性系数;D为热扩散率。

2 热格子Boltzmann方法

2.1 模拟速度场的格子Boltzmann方程

模拟速度场的格子Boltzmann方程采用9个离散速度。在连续的玻尔兹曼方程的基础上, 通过对时间和空间的离散并采取单一松弛时间方法, 得到格子BGK模型:

$\begin{array}{l}g{}_i(x+c{}_i\text{Δ}t,t+\text{Δ}t)-g{}_i(x,t)=\\ -\frac{1}{\tau {}_u}[g{}_i(x,t)-g{}_i^{(eq)}(x,t)],(1)\end{array}$

式中:i=0～8;gi为单粒子分布函数;g${}_i^{(eq)}$为相应的平衡态单粒子分布函数。

定义为:

式中:$s{}_i(u)=\omega {}_i[\frac{c{}_i\cdot u}{c{}_s^2}+\frac{(c{}_i\cdot u){}^2}{2c{}_s^4}+\frac{u{}^2}{2c{}_s^2}];$

$\{\begin{array}{l}\lambda +\gamma =\sigma \\ \lambda +2\gamma =\frac{1}{2}。\end{array}$

通常, $\sigma =\frac{5}{12},\lambda =\frac{1}{3},\gamma =\frac{1}{12}$;

$c{}_i=\{\begin{array}{l}(0,0)i=0\\ (\cos [(i-1)\text{π}/2],\sin [(i-1)\text{π}/2])ci=1~4\\ \sqrt{2}(\cos [(i-5)\text{π}/2+\text{π}/4],\sin [(i-5)\text{π}/2+\text{π}/4])ci=5~8\end{array}$

2.2 模拟温度场的格子Boltzmann方程

模拟温度场的格子Boltzmann方程使用5个离散速度。模型的演化方程为:

式中:i=0～4, Ti为温度分布函数;T${}_i^{(eq)}$为其平衡态分布函数。

定义为:

${\rm T}{}_i^{(eq)}=\frac{{\rm T}}{5}[1+\frac{5}{2}\frac{c{}_{i}u}{c{}^2}],i=0,1,2,3,4$ (4)

2.3 模型的耦合

许多流动问题 (如自然对流) 可以采用Boussinesq近似。采用这种近似的流体运动方程组称为Boussinesq方程组:

此方程组可以使用由式 (1) 和式 (3) 耦合起来的复合模型模拟。这种耦合是通过在密度演化方程的的右端增加一个力项得到的:

$f{}_i=-\frac{1}{2c}\text{Δ}t\alpha {}_{i}e{}_i\times g\beta ({\rm T}-{\rm T}{}_0)$ (6)

当i=2, 4时, αi=1;否则αi=0。

这样, 密度演化方程变为:

宏观物理量:

$u=\underset{i=1}{{\displaystyle {\displaystyle \sum }^8}}c{}_{i}g{}_i,p=\frac{c{}^2}{4\sigma }[\underset{i=1}{{\displaystyle {\displaystyle \sum }^8}}g{}_i+s{}_0(u)],{\rm T}=\underset{i=1}{{\displaystyle {\displaystyle \sum }^4}}{\rm T}{}_i$

输运系数:

2.4 边界条件的处理

边界条件的处理方法在格子Boltzmann方法中起重要的作用, 对格子Boltzmann模型的精度和稳定性都有很大的影响。本文采用非平衡态外推法[4], 其基本思想是, 将边界节点上的分布函数分解为平衡态和非平衡态两部分, 并根据具体的边界条件定义新的平衡态分布来近似平衡态部分, 而非平衡态部分则使用非平衡态外推确定。

3 数值结果

图2, 3, 4是采用上述的双分布函数模型动态模拟封闭方腔自然对流, 并对不同Ra值下流动和传热机制进行分析。

a) 从流线分布图可以看出, 当Ra较小时, 流动的典型特征是在方腔中央出现一个大涡。随着Ra的增大, 涡逐渐变成椭圆形, 并在Ra达到106时, 流线图变得更加复杂, 椭圆形涡分裂成两个涡并且分别向左壁和右壁移动, 方腔中间出现第三涡, 流体呈现湍流特征。

2) 从等温线分布图可以看出传热机制随Ra增加的变化情况。当Ra较小时, 传热主要是由热壁和冷壁之间的热传导引起的, 在这种情况下等温线几乎是垂直的。随着Ra的增大, 传热逐渐由热传导占统治地位变为对流占统治地位, 因此等温线在方腔中央变得水平, 并且只在热壁和冷壁附近的薄边界层内保持垂直。

3) 从等压线分布图可以看出, 随着Ra的增大等压线逐渐变为水平线, 且压力呈现中间低, 上、下底面高的对称分布。

以上观察到的这些现象与文献[5,6,7,8]中都是一致的, 为了定量比较, 我们分别计算了Ra为103, 104, 105, 106情况下壁面的平均努塞尔数, $\stackrel{-}{{\displaystyle {\rm N}u}}=\frac{1}{{\rm H}}{\displaystyle {\int }_0^{{\rm H}}-}\frac{{\rm H}}{\Delta {\rm T}}\frac{\partial {\rm T}}{\partial y}|{}_{\text{w}\text{a}\text{l}\text{l}}\text{d}y$。结果如表1所示。

从表1中可以看出, 采用本文所介绍的双分布函数TD2G9模型, 所得到的结果与文献[5,6,7,8]是相吻合的, 说明本文提出的数值方法是有一定的可靠性。

4 结论

本文采用一种新的不可压格子Boltzmann模型, 在原有格子Boltzmann离散速度的基础上, 引入一类新的分布函数, 并结合温度分布函数组成一个新的不可压的双分布函数TD2G9模型。因此以此模型动态模拟了二维方腔自然对流的形成和演化过程, 数值结果表明:1) 该模型在一定程度上降低了格子Boltzmann模型由可压缩效应引起的误差, 提高了数值计算结果的精确度;2) 和以往的其它格子Boltzmann有很大的不同, TD2G9模型中, 两类分布函数可以使用不同的格子和离散速度, 使得计算更加灵活, 适应性更强;3) TD2G9除了具有一般双分布模型的特点外, 没有增加任何额外的计算量, 仍然保持了一般双分布模型计算量小和数值稳定的优点;4) 数值结果表明Ra越大, 自然对流的作用越强, 流场的演化越复杂, 并呈现湍流特征, 壁面附近的换热越强, 同时等压线逐渐出现方腔中心低压, 上、下壁面高压的轴对称分布。

摘要：采用耦合温度的方法构建一类新的热格子Boltzmann模型, 并对二维封闭方腔中流体的流动状态和传热进行模拟。数值结果表明, 该模型克服传统模型的可压缩效应, 在没有添加限制条件的情况下, 节省了计算时间和提高了计算精度, 并且以此模型动态模拟二维封闭方腔自然对流的结果与已有文献的结果吻合良好。

关键词：热格子Boltzmann,双分布函数,自然对流

参考文献

[1]Dieter A, Wolf-Gladrow.Lattice-gas cellular automata and latticeBoltzmann models.Springer-Verlag Berlin Heidelberg, 2000.

[2]Sauro Succi.The lattice boltzmann equation for fluid dynamics andbeyond.Numerical Mathematics and Scientific Computation (2001) .

[3]Xiaoyi He, Shiyi Chen, Gary D.et al.Anovel thermal model forthe lattice Boltzmann method in incompressible limit.journal ofcomputational physics., 1998 (146) :282-300.

[4]郭照立, 郑楚光, 李青, 等.流体动力学的格子Boltzmann方法[M].武汉:湖北科学技术出版社, 2002.

[5]D Orazio A, Corcione M, Celata G P.Application to natural con-vection enclosed flows of a lattice Boltzmann BGK model coupledwith a general purpose thermal boundary condition[J].Interna-tional Journal of Thermal Sciences, 2004 (43) :575-586.

[6]de Vahl Davis G.Natural convection of air in a square cavity:abench mark numerical solution[J].Int.J.Numer.Meth.Flu-ids, 1983 (3) :249-264.

[7]Barakos G, Mitsoulis E, Assimacopoulos D.Natural convectionflow in a square cavity revisited:Laminar and turbulent modelswith wall functions[J].Int.J.Numer.Meth.Fluids, 1994 (18) :695-719.

Boltzmann 篇5

路基在荷载作用下, 沉降将随时间发展, 其发展规律主要是通过以下两大类方法来加以描述。第一类是以固结理论及各种土的本构模型为基础的计算沉降量的各种有限元法;第二类是根据现场实测资料来推算沉降量与时间关系的预测方法[1]。但是由于固结理论的假设条件以及确定计算指标在试验技术上的一些问题, 因而在一般工程设计中难以采用, 使得实测路基沉降过程数据在某种意义上较理论计算更为重要和实用。

目前常见的沉降预测方法主要有双曲线法、指数曲线法、抛物线法、三点法、沉降速率法、星野法、Asoaka法、泊松曲线法、灰色理论和人工神经网络等。由于每种沉降预测方法都有一定的适用范围, 需要结合具体的工程项目中沉降变形的特点, 选择合适的预测方法进行沉降预测分析。针对这个问题, 本文结合高速铁路路基工程实例, 采用Boltzmann模型, 对路基后期及最终沉降进行预测。实例研究表明该预测模型具有较好的精度和适应性。

2 路基沉降变形的“S”形特征

路基的沉降按其发展可以分为瞬时沉降Sd (t) 、固结沉降Sc (t) 和次固结沉降Ss (t) 三部分。合理的沉降监测数据曲线应表现为“S”形[2], 即随着加载过程路基沉降可分为四个阶段:

1) 发生阶段:路基刚填筑时, 土体尚处于弹性状态, 路基沉降量随荷载的增加呈近似线性增减。

2) 发展阶段:随着路基填筑, 荷载的不断加大, 使其逐步进入塑性状态。随着塑性区不断开展, 路基沉降速率也不断增加, 直到荷载不再增加为止。

3) 成熟阶段:当路基填筑完成, 部分尚未完成的固结和土体的流变导致沉降随着时间的推移而继续, 但沉降速率递减。

4) 到达极限:随着时间的不断延长, 沉降速率快速减小并趋于稳定。

3 Boltzmann预测模型的建立与分析

3.1 模型的建立

在时间预测序列中, Boltzmann模型的表达式为:

其中, St为t时刻对应的预测值, 其单位为长度单位;t为时间;B为待定参数且为正;A1为最小值;A2为最大值;t0为当St= (A1+A2) /2时所对应的时间, B为无量纲, A1, A2单位为长度单位 (见图1) 。

3.2 模型的特点

Boltzmann预测模型具有以下3个特点, 分别为:

1) 单调递增性, 随着时间t的增长, St也不断增长, 即:

2) 有界性, 当时间t趋于无穷大时, St趋近于A2, 即:

3) 呈“S”形, 存在拐点, 该曲线对时间t呈“S”形, 即:

显然t=t0时, ;当t<t0时, , 曲线开口凹向上方;当t>t0时, , 曲线开口凹向下方, 所以呈“S”形。

4) 良好的适应性, 确定参数A1, A2之后, 通过调节t0和B的值可以模拟相当大范围的曲线, 如图2所示。

3.3 模型参数的求解

Boltzmann预测模型含有3个未知参数A1, A2, B, t0可以通过A1和A2求得。本文采用高斯—牛顿迭代算法求得各参数的最优估计, 通过此方法获得Boltzmann预测模型参数的最小二乘无偏估计[6]。

4 工程实例

4.1 实例介绍

时速350 km/h的沪昆客运专线杭长湖南段线路通过的地区多为岩溶地基、黄土地基、软土地基等不良地层和特殊地层。全段路基工程15.33 km, 占线路总长度的16%, 线路均为无砟轨道, 地基的处理方式有换填、打入桩和CFG桩等。由于高速行车要求路基提供一个高度平顺和稳定的轨下基础, 所以控制路基的沉降变形成为路基设计的关键。因此在路基线下工程施工过程中, 必须加强现场沉降观测和实验分析, 掌握工程特性的变化规律, 及时验证并修正理论计算结果, 确保线下工程沉降变形满足无砟轨道铺设条件[7]。

4.2 预测分析

选取DK876+301断面处的沉降观测值, 以2011年10月~2012年2月的沉降实测数据为基础, 采用高斯—牛顿迭代算法, 利用MATLAB编写程序, 对Boltzmann预测模型各参数进行最优估计, 然后再对指定时间点沉降及最终沉降进行预测[3]。参数的计算结果为:A1=-0.52, A2=6.06, t0=37, B=15.94, 代入式 (1) 得:

应用式 (2) 拟合预测各时间点的沉降量, 路基沉降实测值与时间关系见表1, 路基沉降预测值与实测值对比如表2, 图3所示。

预测稳定后的250 d, 500 d的沉降量分别为6.06 mm, 6.08 mm, 说明该路基沉降已经稳定。预测数据的绝对误差小于0.5 mm, 相关系数达到0.98, 说明Boltzmann预测模型具有较高的精度。

5 结语

1) Boltzmann预测模型所预测的路基沉降量与时间的关系曲线和路基实测沉降变化规律相一致, 都呈“S”形, 符合全过程的沉降量与时间的关系, 说明该模型能够反映路基沉降量与时间的关系[4]。2) 实测沉降值与预测沉降值较吻合, 表明Boltzmann预测模型能够较准确地对路基最终沉降进行预测, 能够满足工程精度要求, 可以为高速铁路路基沉降预测提供有效地参考[8]。

摘要：通过研究高速铁路路基沉降量与时间的关系, 并结合目前工程实际中常用的沉降预测方法, 建立了Boltzmann预测模型, 将该模型应用于具体的工程, 证明了该模型应用于高速铁路路基沉降预测的合理性。

关键词：Boltzmann模型,高速铁路,沉降预测

参考文献

[1]宰金珉, 梅国雄.全过程的沉降量方法研究[J].岩土力学, 2000, 21 (4) :322-325.

[2]梅国雄, 宰金珉, 殷宗泽, 等.沉降时间曲线呈S形的证明——从一维固结理论角度[J].岩土力学, 2004, 25 (1) :20-22.

[3]潘林有, 谢新宇.用曲线拟合的方法预测软土地基沉降[J].岩土力学, 2004, 25 (7) :1053-1058.

[4]余闯, 刘松玉.路堤沉降预测的Gompertz模型应用研究[J].岩土力学, 2005, 26 (1) :82-86.

[5]雷长顺, 肖世伟.高速铁路预应力管桩地基沉降理论分析与计算[J].山西建筑, 2012, 38 (3) :146-147.

[6]夏元友, 刘鹏, 莫介臻.高速公路软基沉降预测系统及其应用研究[J].公路, 2005 (8) :275-279.

[7]王炳龙.高速铁路软土路基工后沉降的预测与控制[D].上海:同济大学博士学位论文, 2003.

Boltzmann 篇6

在过去的十年中,CPU设计者受功率的限制不再在单一CPU上添加更多的晶体管,把主要精力放在发展多核CPU体系结构上。由于图形处理本身就是一个并行任务,因而采用流处理体系结构的GPU更适合于并行计算,更容易采用多核策略。当代GPU比相同标准的CPU能提供更多数量级的处理器和浮点计算能力,NVIDIA公司的GPU相对Intel的CPU,在浮点计算能力和存储器带宽方面有大约10倍的优势[1]。

1 GPGPU和CUDA的相关研究

一直以来,用GPU进行多领域的模拟计算是一个复杂的问题,需要研究者有深厚的图形学领域的编程功底。随着近来多种编程模型的出现,GPU编程的门槛一步步降低。特别是NVIDIA公司针对自己的GPU产品发布的编程模型:计算统一设备架构(ComputeUnifiedDeviceArchitecture,CUDA)[1]。CUDA是一个更为通用的编程模型,它通过基于C语言的编程接口调用显存,建立线程块和网格,实现数据的并行计算。

在计算流体力学领域,关于GPU的数值模拟运算已日趋成熟。周季夫等人[2]在舍弃pixclbuffer离屏渲染的同时,采用最新的帧缓存对象,多重纹理、多通道渲染和乒乓技术来设计一套基于方腔的LMB数值模拟程序,最终使GPU的计算时间缩短到CPU计算时间的六分之一。CUDA模型出现后,Tolke和Krafczyk[3]利用CUDA实现了三维Lattice-Boltzman多孔介质流动模型。JulienC.Thibault和InancSenocak[4]通过CUDA采用GPU集群模拟Navier-Stokes方程。

2 相关理论

2.1 CUDA模型

GPU是单指令、多数据的多处理器,专为密集型、高度并行化的计算而设计,更多的晶体管被用来进行数据处理。CUDA编程模型中,一个任务被分解为可独立处理的粗粒度子问题,再细分成细粒度的片段,以便通过协作的方法并行解决。在CUDA中GPU是由多组多处理器所构成,其中每组多处理器又包含多个处理器。GPU的SIMD也就意味着每个处理器执行同一指令,操作不同的数据来形成并行的概念。

线程模型是CUDA编程结构中最核心的部分,也是CUDA程序的基本构成。如图1所示,一个CUDA应用程序需要分割为几个在设备上执行的函数,这些函数会被编译成设备指令集,也就是图1中的kernel程序,每个kernel都对应着设备上一个由线程块组成的批处理线程,在kernel程序下载到设备上后就要按照批处理线程块来执行。

设备上的每个线程都只能使用设备的DRAM和On-Chip内存,主要的内存空间种类分为寄存器、本地内存、共享内存、全局内存、常量内存、纹理内存。按照时钟周期和数据读写速度,寄存器最快,本地内存其次,共享内存慢于本地内存,之后是纹理内存与常量内存,最后是全局内存。

2.2 LBM介绍

LBM模型是一种流体运动的介观模型。LBM方法是对连续Boltzmann方程的一个特殊的离散格式的求解。求解过程是时间推进式的,并且求解过程具有良好的区域性,所以很适合以并行的形式来求解。

实际上,对方程的求解在LBM方法中被转化为一系列的规则,它们被称为迁移规则和碰撞规则,这些规则作用于粒子的分布函数。对流场的求解因此转化为对受迁移规则和碰撞规则制约的分布函数的求解。

2.2.1 碰撞规则

$f{}_i^{*}(x,t)=f(x,t)-\frac{1}{\tau }[f{}_i(x,t)-f{}^{(eq)}(x,t)](1)$

2.2.2 迁移规则

fi(x,t+Δt)=f*i(x-ciΔt,t) (2)

其中,i表示粒子离散速度的不同方向,x表示离散节点的位置,fi(x,t)为离散速度方向i的粒子分布函数,fσi(x,t)为粒子的平衡分布函数,Δt为时间步长。ci为离散方向上粒子的速度。f*i(x,t)为粒子碰撞更新后的平衡分布函数,τ表示松弛时间,计算公式为:

$\tau =\frac{3U{}_{0}L}{Rec{}^2\Delta t}+0.5(3)$

式(3)中U0为特征速度,L为特征长度。Re为雷诺数,为流体力学中惯性力与粘性力的比值,在算例中为给定的值。

在图2所示的Lattice-Boltzmann的D2Q9模型中,平衡分布函数具体定义下:

$f{}_i^{(eq)}(x,t)=w{}_i\rho \left[1+\frac{3c{}_{i,\alpha }u{}_{\alpha }}{c{}^2}+\frac{9(c{}_{i,\alpha }u{}_{\alpha }){}^2}{2c{}^4}-\frac{3u{}^2}{2c{}^2}\right](4)$

式(4)中w0=4/9,w1=…w4=1/9,w5=…w8=1/36,在每个格点中都有9个不同方向的平衡分布函数。

格点上的宏观变量如密度、压强、速度等都可以由分布函数得到,具体计算式为:

$\rho ={\displaystyle \sum _{i}f}{}_i,p=\frac{\rho }{3}c{}^2=\rho c{}_s^2,u=\frac{1}{\rho }{\displaystyle \sum _{i}f}{}_{i}c{}_i(5)$

式(5)中$c{}_s=\frac{c}{\sqrt{3}}$为声速。

2.2.3 边界条件的设置

边界条件处理在LBM的应用中占有重要地位,它不但影响计算结果的精度,还影响计算的稳定性。一般的边界条件的格式有反弹格式、无滑移格式、外推格式。

LBM方法的演化方程可看作连续Boltzmann方程的一种特殊有限差分格式。基于这种观点,提出了边界处理的外推格式。外推法具有二阶精度,意味着比其他方法更快的收敛速度,但外推格式的稳定性较差,因此Guo等人[5]提出了一种非平衡态外推格式,并证明这种格式同样具有二阶精度。

其基本思想是将边界格点上的未知分布函数分解为平衡态和非平衡态两部分,然后利用一阶精度的外推格式得到非平衡部分。而平衡部分由物理边界条件求得。该格式可以表示为:

f${}_i^0$-f${}_i^{0.(eq)}$=f${}_i^1$-f${}_i^{1.(eq)}$ (6)

式(6)中f${}_i^0$,f${}_i^1$分别表示真实边界上的格点和紧邻的流体格点外的分布函数,f${}_i^1$,f${}_i^{1.(eq)}$分别为相应的平衡态分布。本文即采用非平衡外推的边界条件设置。

3 计算流程

在方腔模型算例中,主要变量有:9个分布函数,X和Y两个方向的速度,每个格点粒子的密度等。计算的基本步骤为:先将分布函数、速度、位置都以二维数组存储,在CPU端(Host端)和GPU端(Device端)分配存储空间,初始化数据后通过cudaMemcpy2D()函数把Host端数据传递到Device端。在Device端按照划分好的线程块进行并行计算,迭代数次后再将最终结果传递回Host端,按常规方法输出即可得到结果。在Device端的内核程序设计中,按照方腔算例的过程设计了4个内核函数,分别处理粒子碰撞过程、粒子迁移过程、左右下边界处理、上边界处理。

依次用以下4个内核来模拟方腔的计算过程。

(1) 粒子碰撞过程。

在这个过程中按照公式(4)对非固壁粒子计算平衡分布函数,然后重新计算分布函数。并设计三个一维数组记录固壁邻边粒子的平衡分布函数和分布函数,以供边界处理的内核函数调用。

(2) 粒子迁移过程。

把纹理绑定到cuda数组后,将上一步内核中计算的分布函数赋给cuda数组,通过公式(2)将分布函数进行交换,即将图2中1和3,2和4,5和7,6和8这几个方向的分布函数互换。

(3) 左右下边界处理。

运用非平衡外推的边界处理公式(6),调用粒子碰撞过程中记录的固壁邻边的数据,对三个边界的分布函数重新进行计算。完成后需要对左下角和右下角两个点进行单独处理,因为这两个点的邻边也是固壁,所以需要选择这两个点斜边上的邻点作为计算对象。

(4) 上边界处理。

通过公式重新给上边界的各个粒子给予初速度,完成下一步的粒子碰撞。

在device端的内存分配细节时,考虑到计算的是二维数组,因此没有采用cudaMalloc()和cudaMemcpy()函数分配和复制内存,而采用cudaMallocPitch()和cudaMemcpy2D()函数。cudaMallocPitch()能确保合理填充已分配的存储器,满足对齐要求,从而确保访问行地址或者执行二维数组与设备存储器的其它区域之间的复制时获得最优性能。用cudaMemcpy()函数复制内存时会将每个行都分开复制,导致调用很多次cudaMemcpy()函数,降低计算效率。而cudaMemcpy2D函数复制二维数组时可以指定数组的pitch,从而通过一次函数调用就可将数组全部复制。

本程序中使用全局内存和纹理内存,没有选取共享内存。共享内存虽然存取数据时比全局内存快,但需要先将数据从全局内存复制到共享内存,计算结束后再将数据复制到全局内存。而在本算例中对每个网格点上的数据仅进行了一次运算,并没有反复进行多次存取,所以选用共享内存对运算速度并没有提升。纹理内存是只读内存,可以将数组绑定到其上,通过纹理缓存加速读取速度,速度比全局内存快。本程序中离子迁移过程仅将分布函数的值进行交换,更适合用纹理内存。

4 实验结果和效率对比

本文的实验采用Intel®CoreTM 2CPU6300@1.86 GHz,2GB内存,NVIDIAGeForce9600GT显卡,512 M显存,显卡的核心频率为650 MHz,操作系统为WindowsXP平台。

本实验用方腔中的二维黏性流动作为算例,分析和验证LBM的计算能力。在方腔中,不可压流体被一个正方形围栏所限制,并且流动被顶部匀速的平移所驱动,在方腔里所产生的流体运动是一个封闭的流线问题。本文根据所用显卡的性能,为每个线程块划分16×16个线程,分别测试了网格数为雷诺数为1 000,网格数为128×128、256×256、512×512、768×768、1 024×1 024的方腔模型迭代1 000、5 000、10 000、50 000、100 000次时GPU相对CPU的加速比。不但顺利完成了基于CUDA的LBM方程数值计算,而且提高了GPU的运算速度,效率比CPU计算有大幅提升。

图3所示为256×256网格迭代100 000次时方腔内部流体流线模拟结果。

从图4可以看出,当网格数为128×128时GPU的加速比已经达到40左右,随着网格数的增加,加速比逐渐增大,256×256网格加速比为53左右,512×512网格加速比为60左右,768×768网格加速比达到65左右,这说明随着网格数量的增加,CUDA为GPU分配的线程数也增加,GPU的多线程并行的特性发挥得更加明显,相对CPU的加速比也就更大。但当1 024×1 024网格时加速比降为59左右,说明网格数增大到一定程度时由于硬件的限制不能再分配更多的线程,导致加速比下降,但相对CPU来说,速度的提升仍然相当可观。

5 结论

本文利用CUDA加速求解顶盖驱动方腔内流体的运动控制方程,实现了流体流动模拟,并设计了一种提高GPU迭代计算LBM加速比的算法。通过实验可以看出,加速比可以达到50倍,说明利用CUDA进行LBM计算不仅可行,而且相对CPU来说速度性能非常可观。

摘要：对Lattice Boltzmann方法(LBM)在CUDA下的建模和算法进行了研究,使得该方法在GPU下的计算速度得到提升,大大缩短了计算过程的时间消耗。利用非平衡外推边界条件处理,以LBM方法模拟了D2Q9模型的方腔顶盖驱动流动,采用全局内存和纹理内存存储数据,将模型中9个分布函数存储为二维网格,每个网格分配一个线程,每个线程块包括256个线程,多条线程并行计算。在普通个人计算机上,采用NVIDIA GeForce 9600 GT显卡和CUDA,实现了LBM模拟方腔流动,将计算速度提高到CPU的50倍。

关键词：计算统一设备架构(Compute Unified Device Architoctune,CUDA),GPU,Lattice Boltzmann方法,非平衡外推,边界处理,纹理内存,多线程,并行计算

参考文献

[1] NVIDIA Corporation NVIDIA CUDA Computer Unified Device Ar-chitecture Programming Guide.2nd ed.NVIDIA,2008

[2]周季夫,钟诚文,尹世群,等.基于GPGPU的Lattice-Boltzmann数值模拟算法.计算机辅助设计与图形学学报,2008;20(7):912—918

[3] Tolke J,Krafczyk M.TeraFLOP.computing on a desktop PCwith GPUs for 3D CFD.International Journal of ComputationalFluid Dyanmics,2008;22(7):443—456

[4] Thibault J C,Senocak I.CUDA implementation of a Navier-Stokessolver on MultiGPU desktop platforms for incompressible flows.47th AIAA Aerospace Sciences Meeting Including The New Hori-zons Forum and Aerospace Exposition.Florida 2009:1—15

Boltzmann 篇7

在宏观条件下,经典的流体力学理论认为,黏性流体在固-液界面上均表现为无滑移边界,这种无滑移边界假设的唯一依据是,在连续假设的条件下,理论与实验观察结果相一致.在处理宏观流动问题时,无滑移边界已被广泛用于各种工程和实验中[1]但当流场特征尺寸为微米或亚微米量级时,一些宏观流动中可以被忽略的因素如表面张力、润湿性、相对粗糙度等可能在微流动中占据主导地位,从而导致各种奇异的微尺度流动现象,例如边界滑移效应.

针对不同的问题,流体的运动可以从宏观或者微观的角度来研究.从16世纪以来,Newton,Poisson,Stokes,Euler和Navier等物理学家将流体视为连续的整体,流体的速度、密度、压力和温度等宏观物理量是空间和时间的连续函数,通过质量、动量和能量守恒,可以推导出一组描述流体运动的非线性偏微分方程.用微积分方法计算流体运动的参数,尽管采用了不同的研究对象和研究体系(如Lagrange法和Euler法),但他们所建立的流体运动平衡都是基于共同的连续介质模型,而忽略了离散过程中某些物理量的守恒性.微观尺度的流动模拟是把流体视为分子的集合,在此基础上产生了基于分子的模型.这些模型分为确定性方法和统计性方法两大类[2].在微纳尺度流动中,克努森数(Kn)是一个重要参数,Kn被定义为流体分子的平均自由程与流场的特征尺度之比,被认为是流体流动是否满足连续性条件的标准.根据Kn取值范围确定适用的数值模拟方法,如图1所示.需要强调的是,在微尺度流动中由于特征尺度很小,流动阻力较大,流动速度一般都很低,虽然也可以按照连续性假设条件求解,但是由于比表面积大,表面润湿性和表面滑移等在宏观问题中可以忽略的问题,在微观条件下必须被重新认识,固体表面的浸润吸附特性和静电力作用也必须考虑[3].

通过分子动力学等数值模拟方法也可以证明滑移现象的存在[5,6],但是分子动力学方法计算尺度小的问题[7]以及低流速与分子热运动的非线性耦合问题[8]限制了它的应用.格子Boltzmann方法继承了格子气自动机的主要原理并对其进行了改进,建立了许多开创性的思想,特别是实现了从模拟流体运动的连续介质模型向离散模型的转变[9].格子Boltzmann方法把计算区域离散为许多格子,格子的尺度远比分子平均自由行程要大,但又比有限差分的步长或有限体积法中的控制单元要小,从而将流体的流动抽象成大量的微观粒子按既定的规则碰撞和运移的过程,并可采用多种格式来处理各种形式的复杂边界,包括非滑移和完全滑移边界.这些粒子既比分子级别要大,但其质量又远比有限体积法中的控制容积质量小.粒子在格子之间运动,满足微观的动量守恒条件.流体宏观层次上的密度、速度等参数是通过对这些粒子的有关特性值进行平均获得,因而格子Boltzmann方法在微观流动或多孔介质的渗流规律模拟方面有较大的优势,并已成为研究流体在多孔介质内流动规律的新手段.

1 微流动格子Boltzmann方法模拟的相关理论

1.1 格子BGK模型

Boltzmann方程与流体力学基本方程有密切的联系,可以通过数值求解Boltzmann方程来模拟流体的宏观运动[10].求解Boltzmann的难点在于计算碰撞项,为了简化求解过程,1954年由Bhatnagar,Gross和Krook提出一种碰撞函数模型方法[11],用一个简单的算子来代替碰撞项,认为碰撞的效应是改变分布函数并使其趋于Maxwell平衡态分布.格子Boltzmann方程是Boltzman-BGK方程的一种特殊的离散形式.这一离散包括速度离散,时间离散和空间的离散.离散的时间和空间通过粒子的离散速度联系起来,这一特征把粒子在物理空间的运动分为迁移和碰撞,即粒子在两个时间步之间由一个节点恰好运动到对应的相邻节点,并在该节点上与其他粒子碰撞,这使得格子Boltzmann方法具备了很好的并行计算特性和较强复杂边界的处理能力[10].单松弛时间的格子Boltzman演化方程为

式中,τ是松弛因子,Ci是离散速度,fi是粒子分布函数,fieq是平衡态分布函数.

1992年,Qian等人提出DnQb模型(n代表空间维数,b是离散速度数),其平衡态分布函数统一写为[12]

式中,cs为格子声速,ωi是权系数,这两个参数决定了格子取值依赖于速度离散模型ci的选取.对于D2Q9速度离散模型

权函数与格子声速为

其中,c△x/△t为格子速度,△x和△t分别为格子步长和时间步长.

流体的宏观密度、宏观动量分别为

流体黏性系数和模型松弛因子之间的关系如下式所示

1.2 Shan-Chen模型

格子Boltzmann方法反映的是微观粒子的守恒特性,所以可以方便地描述不同相之间的相互作用,是研究多相流系统的一种有效途径.用格子Boltzmann方法模拟两相流主要有两种方法:一种是从系统自由能角度求得分子间作用力,即自由能方法[13,14],该方法能够满足局部质量和动量守恒,但是在密度变化较大的区域如两相界面处,不满足伽利略不变性,会导致一些非物理现象出现[15].另一种是以Shan-Chen模型[16,17,18]为代表的引入微观分子间作用力的方法,这种方法能够反应出流体动力学的本质,且易于耦合微流动问题中占主导作用的微观作用力.根据研究问题的不同,可分为单组分模型和多组分模型.

(1)单组分模型

在单组分Shan-Chen模型中,微观分子间作用力为

流体与固体之间的作用力为

式中ωi是权函数,其取值如下

式(9)中G是两相液体相互作用强度,式(10)中Gads是吸附参数.当x+ci△t是流体所在位置时,s为0;当x+ci△t是固体所在位置时,s为1.可以通过调整G来改变流体的表面张力.在确定了参数G的基础上可以通过调整参数Gads来调整液体在固体表面的接触角大小.Ψ(x,t)是相互作用势能,其它形式的相互作用势[16,17]也经常被采用.但是这些作用势不满足严格的等温过程,为了满足等温过程,可采用如下形式的相互作用势能[17]

格子Boltzmann方法中引入力项有多种方法,原始的Shan-Chen模型通过改变平衡态分布函数中的平衡态速度来反映粒子间相互作用力的影响.平衡态速度为

流体的真实宏观速度为

此时状态方程为[19]

(2)多组分模型

单组分多相流模型较适合于近似模拟气液两相系统,例如水、气系统,而对于多组分系统例如油、气、水系统需要用到多组分的Shan-Chen多相流模型.在多组分Shan-Chen模型中,将第σ组分的密度设为,从而得到平衡态速度

这里u'被定义为

第σ组分的速度按以下格式计算

作用在第σ组分的黏附力为[20]

式中,σ和表示两种不同的组分,Gc是用来控制两组分液体之间黏性力大小的参数.

作用在第σ组分的表面力为[20]

Gads,σ代表第σ种组分与壁面的作用强度,以前研究结果[20,21,22,23]表明对于非润湿流体Gads,σ取正值,而对于润湿流体Gads,σ取负值[20,21,22,23].

2 基于Shan-Chen模型的格子Boltz-mann方法在微流动研究中的应用

Shan-Chen模型直接对微观相互作用力进行描述,能够反映多相或多组分流体动力学的物理本质,因而得到了比较广泛的应用.这里限于篇幅,重点分析该方法在接触角、微观流、多孔介质微流动方面的应用,以便为石油储层微流动研究提供参考.

2.1 接触角模拟

接触角是指在气、液、固三相交点处所作的气-液界面的切线穿过液体与固-液交界线之间的夹角,是液体润湿固体程度的量度,Benzi率先对单组分Shan-Chen模型进行了理论分析[24],给出了适合Shan-Chen模型的接触角理论计算公式.接触角滞后的研究也具有重要的理论和应用价值,因为它不仅是判断一个表面是否具有超疏水性的重要标准[25],还可以看作液体在表面黏附性能的一个量度[26],Hyväluoma采用Shan-Chen模型模拟了液滴在倾斜表面的形状[27],研究了粗糙表面对接触角滞后的影响,结果表明增加粗糙壁面的疏水性,接触角滞后将减小,同时也指出接触角滞后是研究壁面润湿性的重要参考依据.Huang[21]利用多组分的Shan-Chen模型,模拟研究了三相接触角与参数之间的关系,得到经验公式

在粗糙表面,接触角的理论模型是Wenzel模型[28]和C-B模型[29].Patankar[30]和Murmur[31]分别发表文章,根据Gibbs自由能最小理论分析了两个模型的适用条件.Bico等[32]从实验发现从C-B模型向Wenzel模型转化需要额外作功的现象,Patankar则就实验现象从理论上进行了分析[30],Dupuis通过自由能模型的格子Boltzmann方法模拟了从C-B模型向Wenzel模型转化[33].

接触角是用来表征表面润湿性的关键参数,接触角的大小反映了表面润湿性的强弱,这就提供了一种新的模拟思路:当采用Shan-Chen多相流模型模拟工程问题时,首先通过接触角模拟可以确定吸附参数与接触角的关系,从而得到满足任意润湿性的壁面,从模拟角度确定了吸附参数就相当于确定了满足条件的材料属性,进而可以模拟流体在由这种材料属性构成的系统中的运动特性,这是数值模拟的优点所在.

2.2 微管道流模拟

随着科技的发展,研究对象的微型化是近20年来自然科学和工程技术发展的一个重要趋势,微纳米尺度通道内的流动问题也相应地越来越受到关注.微尺度和宏观尺度上的流体力学的主要区别,大致可分为4个方面:非连续效应、表面主导效应、低雷诺效应、多尺度和多物理场耦合效应等[34].流体相态不同,影响微尺度流动的主要因素就不同.对于气体而言,流体的奇特性主要是由气体的稀薄性引起的;而微尺度液态流体的特性则主要受表面张力和分子间力的影响[15].由于尺度问题,无滑移的宏观流体边界条件受到了挑战,在微观条件下,滑移效应不可以忽略,越来越多的实验研究和模拟结果也都表明了固/液界面上滑移现象的存在[7,35].基于滑移效应的人工超疏水表面已经在工业生产和人们的日常生活中得到越来越广泛的应用[36].例如,石油储层微通道的表面改性及注水作业中的降压减阻技术[37,38]、汽车挡风玻璃、建筑物的门窗、天线表面的防雪防雨[39]等.

对于这种滑移现象,Navier[40]早在1823年就提出了滑移速度边界条件,通过在固-液边界设置一个滑移速度,来解决工程实际问题.Succi等[41]把这种滑移边界条件引入到格子Boltzmann方法中,但只能适合直边界情况,随后有学者[42]对其进行了改进,使之适合于任意曲面边界.因为滑移速度的大小很难精确确定,所以计算结果的精确程度受到很大的影响.基于固体力学的塑性流动的启示,吴承伟等[43]发现界面滑移与界面剪贴强度密切相关,提出了应力控制边界滑移模型,这种模型的理论分析基础仍然是建立在宏观经典流体力学的连续介质理论.陈艳燕等[44]通过采用Shan-Chen单组分模型,对Couette流进行了模拟研究,指出滑移长度的大小与固体壁面和流体之间的作用力有直接关系,并给出流体与固体相互作用强度与滑移长度的关系,但是这种流体与固体之间的相互作用难于和实际物理参数相比较,所以很难应用于工程实际.

运用Shan-Chen单组分模型,Raiskinmaki等人研究了毛细现象[45],通过改变毛细管的直径,研究由于润湿性的原因,毛细管内液面上升的高度与毛细管直径的关系,模拟结果与理论结果表现一致.此外狄勤丰课题组研究了液滴在二维微管道内的运动[46,47],研究结果表明接触角越大,流动速度越快,减阻越显著;通过对液滴在粗糙壁面上运动的模拟,发现基底的润湿性对流动特性影响显著,同时也表明亲水材料构成的壁面,表面粗糙度使壁面更亲水,而由疏水材料构成的壁面表面粗糙度使壁面更疏水,模拟结果与实验结果相一致[48].

Kang等[23]采用Shan-Chen多组分模型模拟了两种不相溶液体在重力作用下的动力学过程,主要研究了接触角、Bo数(重力与表面张力的比)、黏度比、密度比对流动特性的影响.2006年,Harting等[49]用Shan-Chen多组分模型模拟了疏水管道壁面和流体间的相互作用,得到了与小Knudsen数时的实验结果[50]和其他数值模拟结果[51,52,53]相同的结论,即随着流固相互疏水作用的增强滑移效果增大,滑移长度只与壁面-流体之间的作用强度有关,而与流体速度无关.Kunert和Harting[54,55,56]研究了微通道流中粗糙度引起边界滑移的现象,引进了等效无滑移面的概念,模拟了粗糙疏水微管中的流体流动,如图2所示,并发现表面粗糙度和表面疏水性会影响等效无滑移面的位置,壁面粗糙度与滑移长度间存在非线性的关系,在高斯分布粗糙度情况下等效无滑移面的位置与分布的宽度成线性关系,给出了平均覆盖率与滑移长度的关系,如图3所示.文献[56]还研究了黏度对滑移长度的影响,如图4所示.

2.3 多孔介质微流动模拟

多孔介质是指固体材料中含有孔隙、微裂缝等各种类型毛细管体系的介质,流体在多孔介质中的流动称之为渗流.格子Boltzmann方法本身局部计算的特性,天然适合于并行计算,并且具有对边界处理简单的特点[57,58],使得格子Boltzmann方法在渗流微观机理研究方面具有很大的优势,并且随着其有效性得到不断的验证,逐渐被人们认为是可与实际实验相媲美的一种方法.这类数值研究与微观物理模拟使人们可以从微观的角度来解释宏观现象,从微观层次加深人们对渗流机理的认识[59].姚军课题组在油藏多孔介质渗流领域采用数值模拟方法,基于模拟退火算法[60]和择多算子的随机搜索方法[61]构建数字岩心,然后通过格子Boltzmann方法对构建的数字岩心进行评价[62],分析了岩心内流场的分布并给出了岩心的渗透率值,充分体现了数值模拟技术的高效性.钱吉裕等[63]通过格子Boltzmann方法确定了多孔介质流动参数,验证了Darcy定律.许友生[64]根据格子Boltzmann方法及相关理论,建立了一个新的模拟渗流运动的数值模型[65],该模型没有在边界上采取相应平均措施,同时避免了一些非物理副产品的出现,实例计算数值结果与精确解符合较好.阎光武等[66]对多孔介质中的流体流动进行了模拟,通过研究流量和压力梯度的关系,给出了格子方法模拟的Darcy定律,最后通过调节孔隙率,给出黏性系数与渗透率的关系曲线.这些模拟都是建立在无滑移边界条件上的,没有考虑壁面的润湿性.考虑液-液之间和固-液之间的相互作用,Martys等[20]采用多组分Shan-Chen模型模拟了两种不相溶液体在多孔介质中的驱替过程,研究了不相溶液体在多孔介质中的体积比与相对渗透率的关系,模拟结果与实验结果相一致.还有学者利用Shan-Chen模型的优势对复杂多孔介质内气液两相流的流动[67]进行了研究,分析了多孔介质中不同气液体积比时液体的输运现象.此外,有学者运用单组分Shan-Chen模型研究了液体在纸张中的渗透[68]和在简单多孔介质中模拟了两相驱替的过程[69,70],均得到了很好的结果.文献[71]采用ShanChen单组分多相模型从毛细数、黏度比和接触角3个方面研究了多孔介质中不相容两相液体的黏度耦合对流动特性的影响.结果表明当润湿相比非润湿相黏度小时,非润湿相的相对渗透率要大.文献[22]采用Shan-Chen多组分模型,通过Cα数(黏性力与表面张力的比)、M数(两相黏度比)、Bo数(重力与表面张力的比)3个无量纲数根据相似原理研究了理想非均质多孔介质的毛细压力与饱和度的关系.

2.4 其他方面的应用

随着研究的深入,Shan-Chen模型越来越得到发展.最初的Shan-Chen模型由于伪速度的存在,并不能模拟大密度比多相流体,但是Yuan等提出可以通过选择适当的状态方程来提高液气密度比[72].有人还提出通过精确的数值方法与双松弛时间方法[73],来减少伪速度,使密度比变大.原始的ShanChen单组分模型的表面张力是温度的函数[17],这就意味着不论采取多大的特征尺度,只要代表温度的参数G—定,气液密度比就一定,表面张力就一定.这种模型符合等温过程的物理意义,但该模型由于气体状态方程和表面张力不互相独立,又限制了它的应用.在求解两相流体之间作用力方程(9)时,通过采用不同离散方法[74,75],引入另外一个或多个参数使表面张力独立于温度参数,虽然这些方法备受争议[76],但是对于满足相似原理的流体力学问题无疑是更方便和能说明问题的.为了对非等温过程的模拟,Zhang等[77]用传统的离散方法把能量方程加入到修改后的单组分Shan-Chen模型中,模拟了液体沸腾现象.对于耦合电场的作用效果问题,有学者采用Shan-Chen模型研究了表面带电液滴在外加电场作用下的变形和动运[78].另外,耦合了悬浮颗粒模型[79],有学者对颗粒在气液自由面的动力学过程进行了模拟.总之,在微纳尺度下,液体的表面张力,静电力,范德华力,Casimir力等占主导作用,而Shan-Chen模型易于耦合这些力的作用,并且相对于其他双分布函数模型更易于编程实现,使得Shan-Chen模型得到了广泛应用.

3 结语

从上面的分析和叙述中可以看出,基于Shan-Chen模型的格子Boltzmann方法在微流动模拟研究中的应用越来越受到人们的关注和重视,并正发挥着愈来愈重要的作用.

首先,基于Shan-Chen模型格子Boltzmann方法是一个适合于模拟微流动问题的有效工具.对于微流动的模拟,其特征尺度一般在微米级范围,其Mach数和Reynolds数一般都不大,属于低Mach数层流状态.而格子Boltzmann方法是一种有别于传统计算流体力学方法的新方法,它基于统计物理,以极其简单的形式描述粒子的微观行为,来反映流体在宏观上的运动.由于该方法高效并行性、可以方便地对边界进行处理、程序编制简单、计算准确等特点,越来越被人们所关注并接受,它比宏观N-S方程更接近微观层次,更能反映流动的细节.Shan-Chen模型易于耦合微观条件下占主导作用微观力,能够反映多相/多组分流体动力学的物理本质,从而很好地拓宽了格子Boltzmann方法在模拟微流动问题方面的有效性.

其次,基于Shan-Chen模型的格子Boltzmann方法在微结构表面的滑移效应模拟方面具有很好的应用空间.微结构表面的滑移效应具有非常广泛的应用前景,从本质上来说,这种滑移现象是由流体在壁面的润湿性决定的,而Shan-Chen模型适合具有润湿性边界条件的滑移模拟.对于这一类问题,在满足相似理论(例如在达西渗流研究过程中,当研究黏性的影响时一般参考无量纲Re数,而当研究表面张力的影响时则主要参考无量纲Ca数)条件下,这种方法就具有很好的模拟效果.

摘要：对格子Boltzmann方法的本质及Shan-Chen模型的核心机制进行了全面阐述,并从应用实例角度对基于Shan-Chen模型的格子Boltzmann方法在微流动模拟方面的有效性、适应性进行了详细分析.结果表明,Shan-Chen模型易于耦合微观条件下占主导作用的微观力,拓宽了格子Boltzmann方法在微流动模拟方面的应用.同时,Shan-Chen模型在润湿性边界条件表征方面的优势,使得这种方法在微结构表面的滑移效应模拟方面具有很好的应用前景.