欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页  >  IT编程

LBM中的straight boundary及部分代码(以D2Q9为例)

程序员文章站 2022-04-02 11:15:39
本文将从何雅玲老师的《格子Boltzmann方法的理论和应用》、Yuanxun Bill&Justin Meskas的《Lattice Boltzmann Method for Fluid SImulations》、Timm Kruger的《The Lattice Boltzmann Method》以及默罕默德的《格子玻尔兹曼方法——基础与工程应用(附计算机代码)几本书中挑选一些重要的或常用到的直边界进行介绍。值得一提的是,Timm Kruger里把边界分成两类:link-wise和wet-.....

本文将从何雅玲老师的《格子Boltzmann方法的理论和应用》、Yuanxun Bill&Justin Meskas 的《Lattice Boltzmann Method for Fluid SImulations》、Timm Kruger 的《The Lattice Boltzmann Method》以及默罕默德的《格子玻尔兹曼方法——基础与工程应用(附计算机代码)几本书中挑选一些重要的或常用到的直边界进行介绍。

值得一提的是,Timm Kruger里把边界分成两类:link-wise和wet-node。他俩的区别可见下图。Link-wise包括反弹边界条件,wet-node包括平衡格式,非平衡插值格式,非平衡反弹格式等。而何雅玲老师和*的《格子玻尔兹曼方法——基础与工程应用(附计算机代码)》等里并没有这两个概念。

LBM中的straight boundary及部分代码(以D2Q9为例)

在边界条件中应注意未知的是哪几个方向的分布函数,左边界是f1 f5 f8,右边界是f3 f6 f7,上边界是f4 f7 f8,下边界是f2 f5 f6。当然,在每个边界上也可以把个方向的分布函数都给定,这并没有什么影响(这段话主要是针对我学习时搞不懂为什么有的人给了3个,有的给了4个,有的给了8个。。。。后来才发现多给了是不要紧的)。

角点处需要特殊处理,有时往往是角点处理不好导致计算发散或结果不理想。角点处会专门介绍。

 

  1. 周期性边界

(1)标准周期性边界

适用条件:流场在空间呈现周期性变化或在某个方向无穷大(注意:是流场呈周期而不是几何呈周期)

表达式:周期性边界常常需要在左右两侧增加虚拟节点x0和xN+1

LBM中的straight boundary及部分代码(以D2Q9为例)

代码:注意每种语言的下角标不同

for (j = 1; j < NY; j++)

{

    f[0]][j][1] = f[NX+1][j][1];

    f[0]][j][5] = f[NX + 1][j][5];

    f[0]][j][8] = f[NX + 1][j][8];

}

尽管Timm里也提到了一种不需要虚拟节点的周期性边界,但仍建议使用虚拟节点。因此在此不介绍无虚拟节点的周期性边界。

(2)广义周期性边界条件

适用条件:速度呈现周期性但压力场/速度场非周期

表达式:分布函数分解成两部分——平衡态和非平衡态

LBM中的straight boundary及部分代码(以D2Q9为例)

 

平衡态部分:

LBM中的straight boundary及部分代码(以D2Q9为例)

非平衡态部分:

LBM中的straight boundary及部分代码(以D2Q9为例)

 

 

因此可得:

LBM中的straight boundary及部分代码(以D2Q9为例)

2、反弹格式

适用条件:固体壁面,无滑移。其主要思想是:分布函数在迁移过程中撞到固体壁面后原路弹回。

(1)标准反弹格式

2个时间步完成迁移-反弹-迁移的过程

(2)半步长反弹

1个时间步内完成迁移-反弹-碰撞的过程

LBM中的straight boundary及部分代码(以D2Q9为例)

3)静止壁面

以下边界为例:需要确定的是f2 f5 f6

LBM中的straight boundary及部分代码(以D2Q9为例)

代码:

for (i = 0; i <= NX; i++)

         {

                   f[i][0][2] = f[i][0][4];//下边界

                   f[i][0][5] = f[i][0][7];

                   f[i][0][6] = f[i][0][8];

}

4)移动壁面

LBM中的straight boundary及部分代码(以D2Q9为例)

rhow是流体在固体壁面处的速度,与固体的实际物密度无关。可以是局部流体密度,或者系统平均密度。

 

总结:

  1. 对于稳态流动,标准反弹与半步长反弹结果接近;但对于非稳态流动,半步长反弹更准确;通常认为标准反弹是一阶精度,半步长反弹是二阶精度;半步长反弹的另一个优点是不需要固体node,遇到边界后直接弹回
  2. 建议使用半步长反弹
  3. 反弹格式优点:稳定,对于静止壁面严格保证质量守恒
  4. 缺点:当反弹格式与BGK模型一起使用时,无滑移边界实际位置与粘度有关(?)

 

3、非平衡态反弹格式(Zou-He边界)

以速度边界为例。

联立质量守恒、x方向动量守恒、y方向动量守恒三个关系式,并假设与边界垂直的方向上,反弹格式对非平衡部分仍然成立。可得:

西边界:

LBM中的straight boundary及部分代码(以D2Q9为例)

LBM中的straight boundary及部分代码(以D2Q9为例)

LBM中的straight boundary及部分代码(以D2Q9为例)

LBM中的straight boundary及部分代码(以D2Q9为例)

代码:

for (j = 1; j < NY; j++)

         {

                   u[0][j][0] = U;

                   u[0][j][1] = 0;

                   rho[0][j] = (f[0][j][0] + f[0][j][2] + f[0][j][4] + 2 * (f[0][j][3] + f[0][j][6] + f[0][j][7])) / (1 - U);

                   f[0][j][1] = f[0][j][3] + 2 * rho[0][j] * U / 3;

                   f[0][j][5] = f[0][j][7] + rho[0][j] * U / 6;

                   f[0][j][8] = f[0][j][6] + rho[0][j] * U / 6;

}

4、非平衡态外推格式

将分布函数分成平衡态和非平衡态两部分。

平衡态部分:使用边界节点上宏观物理量求得,如果边界节点上的值未知,则使用内部相邻节点的值替代。

非平衡态部分:用内部相邻节点的非平衡部分替代

LBM中的straight boundary及部分代码(以D2Q9为例)

代码:

左右边界:

for (j = 1; j < NY; j++)

                            for (k = 0; k < Q; k++)

                            {

                                     rho[NX][j] = rho[NX - 1][j];   //右边界

                                     u[NX][j][0] = u[NX - 1][j][0];

                                     u[NX][j][1] = u[NX - 1][j][1];

                                     f[NX][j][k] = feq(k, rho[NX][j], u[NX][j]) + f[NX - 1][j][k] - feq(k, rho[NX - 1][j], u[NX - 1][j]);

 

                                     rho[0][j] = rho[1][j];         //左边界

                                     u[0][j][0] = U;

                                     u[0][j][1] = 0;

                                     f[0][j][k] = feq(k, rho[0][j], u[0][j]) + f[1][j][k] - feq(k, rho[1][j], u[1][j]);

                            }

上下边界:

for (i = 0; i <= NX; i++)

                    for (k = 0; k < Q; k++)

                    {

                             u[i][0][0] = 0;//下边界

                             u[i][0][1] = 0;

                             rho[i][0] = rho[i][1];

                             f[i][0][k] = feq(k, rho[i][0], u[i][0]) + f[i][1][k] - feq(k, rho[i][1], u[i][1]);

 

                             rho[i][NY] = rho[i][NY - 1]; //上边界

                             u[i][NY][0] = 0;

                             u[i][NY][1] = 0;

                             f[i][NY][k] = feq(k, rho[i][NY], u[i][NY]) + f[i][NY - 1][k] - feq(k, rho[i][NY - 1], u[i][NY - 1]);

                    }

5、开放边界条件

(1)出口速度未知,通常采用外插法来求未知的分布函数

如右边界:

LBM中的straight boundary及部分代码(以D2Q9为例)

 

 

 

但是,在某些条件下,上述边界条件可能会产生不稳定解,使用一阶内插法比二阶内插法具有更高的稳定性。

代码:

for (j = 1; j < NY; j++)

         {

                   u[NX][j][0] = u[NX - 1][j][0];

                   u[NX][j][1] = 0;

                   rho[NX][j] = rho[NX - 1][j];

                   f[NX][j][3] = 2 * f[NX - 1][j][3] - f[NX - 2][j][3];

                   f[NX][j][6] = 2 * f[NX - 1][j][6] - f[NX - 2][j][6];

                   f[NX][j][7] = 2 * f[NX - 1][j][7] - f[NX - 2][j][7];

         }

(2)另一种方法是假设边界上的压力已知,及密度rho已知。通常,边界上的密度需要设定为常数(如1)。

LBM中的straight boundary及部分代码(以D2Q9为例)

 

本文地址:https://blog.csdn.net/zmyue1314/article/details/107303018