流体力学论文

第3章 渗流微分方程的数值模拟方法

4 嵌入边界条件 如果要求解，必须将边界条件加入。但由于这是局部元素的求解，通解不必将边界条

件嵌入，只作为总体元求解的基础，在总体元求解时边界条件的处理必不可少。

（iv）总体有限元求解

Poisson方程：

L?u??p??u?Q?0

（3.1.2.20）

u?x,2??1,u??0,y??1,ux边界条件：

x?2?0,uyy?0?0

其中

Q?Qw??x?1??(y?1)

（3.1.2.21）

?是Dirac函数。

1 作图并划分元素，注明标号

由边界条件u3?u6?u9?u1?u2?1,只要求出u4,u5,u7,u8即可。

基本公式（3.1.2.13）（3.1.2.15）和（3.1.2.16）。

其中

?2?3a14??1???a24???6a34??1???3a44??1???6?162??131?AiNM?a11?a21???a31??a41a12a22a32a42a13a23a33a4331?31?36231?61?6?1???3? 1??6?2??3?33 --

这是已知局部元的结果，即（3.1.2.18）式，

FMi?1???1?Qw??

?1????1?即（3.1.2.19）式。

2 求Ajk(9?9)

?2?3????????????????14?16?104对3?1?1?1?18根据（3.1.2.15）式（其中I=4）,最后求得：

60?126333330?1?10?14336000?1?102称363000?1?1?1?1433360??0??0????1?3?0??0??1?6?2?3?33Ajk3

33Fk????Qw??x?1???y?1??kd???Qw?kx?1,y?1??Qw?5

3 加入边界条件

?4?3???对????18?1?12称?1?3??1?3??16?4?3??加入本质边界条件后由（3.1.2.13）可得：

3631??u4??2????5?u5???Qw?3?，其中Q?1

w??u??07????1??u8????2??T3

3解得?u4u5u7u8???0.783T0.5250.655T0.783?

0.7190.756而根据解析解可得：?u4u5u7u8???0.7560.127?T，比较可知，u5误

差较大，因为对复杂的温度分布，用线性分布Lagrange插值函数作为基函数不太合适，用

二次插值函数作为基函数可使解更合理，但计算复杂。

3.2 Lattice Boltzmann (LB) 方法模型的建立

从16 世纪以来, New ton、Poisson、Stokes、Euler 和Navier 等物理学家将流体

第3章 渗流微分方程的数值模拟方法

视为不间断的整体,基于连续介质模型假设,采用了不同的研究对象和研究体系(如Lagrange 法和Euler 法),用微积分方法计算流体运动的参数,用一些关于流体物理量,如密度、速度、压力等在时间和空间上的连续函数的偏微分方程(如N avier2Stokes 方程简称N-S 方程) 来描述流体运动,由于N-S 方程是高度非线性化的偏微分方程,仅仅一些具有简单边界或者有比较严格物理限制的现象才能够得到理论分析解，一般情况下严格求解是非常困难。计算流体力学的任务是流体力学的数值模拟,最广泛的数值解法是从偏微分方程出发,采用有限差分、有限元、有限体积等离散方法,得到网格单元上的代数方程组, 而求得微分方程的数值解,是自上而下的求解过程。而Boltzmann 另辟蹊径,从分子运动论的观点出发,以微观粒子运动为基础, 建立离散的速度模型,在满足质量、动量和能量守恒的条件下,得出粒子分布函数,然后对粒子分布函数进行统计计算,得到压力、流速等宏观变量。方法先给出微观粒子团在网格节点上运动与碰撞的演化模型(研究的关键),通过多尺度分析技术,使之再现流场控制方程(描述流体运动的偏微分方程),则建模成功,就可以

用该演化模型进行数值求解,这是一种全新的方法,是自下而上的求解过程。

3.2.1 Lattice Boltzmann 方法的主要思想

将流体存在的区域划分为均匀网格,将流体想象成许多只有质量没有体积的微小粒子组成,在同一时刻同一网格节点上,粒子向九个方向运动(如图9)，移动到最近的网格节点。其中每个节点上允许一个静止粒子存在,加上与其相邻的有8个节点,所以称为二维九点格子模型,记为D2Q9模型。到目前为止已建立的LBM模型有：D1Q3, D2Q9, D2Q7, D2Q13, D3Q15,

这里以二维中常用的D2Q9 模型为例,来介绍Lattice Boltzmann 方法。

D3Q18, D3Q27等(D指维数,Q指粒子运动方向的总数)。

图（9） 二维九点格子模型速度方向以及网格划分

其演化过程主要分两个步骤：(a)迁移,粒子从一个节点在一个时间步长内,以恒定的速度运动到相邻节点；(b)碰撞,在一个节点上从相邻节点运动来的粒子发生碰撞, 根据质量、动量和能量守恒规则改变粒子的速度, 然后各个粒子又以改变后的速度迁移。这两个

35 --

步骤交替循环,直到流场达到收敛。

设?X,t?代表在时刻t, 位置X??x,y?处的节点,流体的密度为???(X,t),流体的速度为u?u?X,t?, 时间步长为?t,t?0,?t,2?t,?,m?t，格子步长为?x,粒子迁移速率c??x?t,

粒子离散速度：

c0??0,0?ci?ci??cos?i,sin?i?c????i?1??22??4i?1,?4?

2?cos?i,sin?i?c??i??i?5??i?5,?8?

粒子分布函数fi??X,t?表示节点?X,t?处运动速度为ci的粒子数量, i?0,1,?,8，则粒子

分布函数演化过程如下：

fiX?ci?t,t??t?fi?X,t????f?X,t??,i?0,1,?,8

iger 等在1989 年作了一个非常重??f?X,t??称为碰撞算子,表示碰撞引起的变化。H

要的简化,假定流场的分布接近于局部平衡状态,并将碰撞操作过程??f?X,t??进行线性化处理,用单一松弛时间使流场逐渐达到局部平衡状态并且线性稳定；后来, Qian、Chen等学者提出,如果松弛形式采用Bhatnagar-Gross-Krook (BGK) 操作就能再现宏观N 一S 方

??程,这就是所谓的Lattice BGK(LBGK)模型,它所对应的LB方程(LBE)是:

fiX?ci?t,t??t?fi?X,t?????1??f?X,t??ifi?eq??X,t??,i?0,1,?,8 (3.2.1.1)

在低Mach(马赫)数的假设下u??cs,其中粒子平衡态分布函数为

?c?uc?u?eq?fi?X,t?????1?i2?i4cs2cs??????????u?u??? （3.2.1.2）

22cs??2且 cs?c3,?0?49,?1??2??3??4?19;?5??6??7??8?136;

其中?是松弛时间尺度,控制达到平衡的速度(可根据需要进行设置),由于稳定性的原

因, 经过实际测算?必须大于12。

事实上不同的网格剖分有着不同的平衡分布函数,LBM 建立模型的核心问题就是根据

不同的网格确定对应的平衡分布函数,对D2Q9模型,我们就取上面的平衡分布函数。

????X,t??在节点(X , t)上根据质量和动量守恒规则,流体的宏观密度、速度、压强定义如下:

?f?X,t???iiiiiifi?eq??X,t? （3.2.1.3）

fi?eq? ?u??u?X,t???cf?X,t???cii?X,t? （3.2.1.4）

p?cs2??X,t? （3.2.1.5）