边界条件分为两种，一种是和网格对齐的，直线型的边界条件。一种则是弯曲型的边界条件。这两种边界条件，都可以在Palabos里面轻松实现。

1 网格对齐边界

1.1 OnLatticeBoundaryConditionXD类

有一些边界条件是局部的，这种情况用Dynamic实现，有些是非局部的，这种通常会产生复合Dynamic，通常要借助于数据处理器。OnLatticeBoundaryConditionXD类相当于提供了一个统一的接口，不管是局部的还是非局部的，或者是这两者都有，这取决于采用的边界条件。不过，不管是什么情况OnLatticeBoundaryConditionXD都有相应的解决方案。

比如可以用下面的语句来实例化一个Zou/He边界条件：

1
2
OnLatticeBoundaryCondition3D < T ，DESCRIPTOR > *  boundaryCondition  = 
        createZouHeBoundaryCondition3D < T ，DESCRIPTOR > （）;

创建了一个OnLatticeBoundaryCondition3D类的指针boundaryCondition，采用的是createZouHeBoundaryCondition3D方法。createZouHeBoundaryCondition3D函数不属于任何类，位于plb里面的一个函数。其定义为：

1
2
3
return new BoundaryConditionInstantiator3D <
 T, Descriptor,
 WrappedZouHeBoundaryManager3D > ();

可见是返回一个BoundaryConditionInstantiator3D对象。

而BoundaryConditionInstantiator3D对象又是OnLatticeBoundaryCondition3D的子类。继承关系如下图：

OnLatticeBoundaryCondition3D的继承关系
OnLatticeBoundaryCondition3D的继承关系

Zou/He是用的相当多的一种边界条件。处理速度、密度、压力边界条件的一种方法。属于非平衡态反弹格式。具有二阶精度和良好的数值稳定性。

目前通过OnLatticeBoundaryConditionXD可以实现的边界条件方法有四种：

BC Method
Regularized BC createLocalBoundaryConditionXD
Skordos BC createInterpBoundaryConditionXD
Zou/He BC createZouHeBoundaryConditionXD
Inamuro BC createInamuroBoundaryConditionXD

四种边界条件的简要介绍可以看官方文档。

利用这些边界条件，可以实现速度Dirichlet边界（所有速度分量具有强加值）和压力边界（施加压力，并且切向速度分量为零）。此外，提出了各种类型的Neumann边界条件。在这些情况下，通过从相邻单元复制该变量的值，以一阶精度施加给定变量的零梯度边界条件。

注：Dirichlet边界条件指函数本身在边界上给定，Neumann边界条件指函数的法向导数在边界上给定。

1.2 常规矩形域的边界

这种应该是最多的吧，求解域是一个矩形，或者是一个立方体。

所有OnLatticeBoundaryConditionXD对象都有两个方法：setVelocityConditionOnBlockBoundaries和setPressureConditionOnBlockBoundaries，用于设置位于矩形域上的边界。如果给定2D框边界框边界上的所有节点都实现了Dirichlet条件，请使用以下命令：

1
2
3
Box2D boundaryBox(x0,x1, y0,y1);
boundaryCondition->setVelocityConditionOnBlockBoundaries (
   lattice, boundaryBox, locationOfBoundaryNodes );

这里先是指定了一个boundaryBox，这个boundaryBox就是边界框了。boundaryCondition其实就是此前创建的那个OnLatticeBoundaryConditionXD对象。在这里调用了这个对象的方法setVelocityConditionOnBlockBoundaries。来看看这个方法对用的声明：

1
2
3
4
5
template<typename T, template<typename U> class Descriptor>
 void ‪OnLatticeBoundaryCondition2D::setVelocityConditionOnBlockBoundaries (
         ‪BlockLattice2D& lattice,
         ‪Box2D block, ‪Box2D applicationDomain,
         ‪boundary::BcType bcType )

结合声明可以知道，第一个参数lattice表示的是创建边界条件所在的block-lattice区域。第二个boundaryBox表示的边界所作用的是一个矩形区域。第三个参数locationOfBoundaryNodes实际上也是一个BoxXD类，它的作用你可以把他看成是一个选择框。如下图所示。locationOfBoundaryNodes 通常用的是一个边。

VCdIbD.png
VCdIbD.png

第四个参数boundary::BcType省略掉了，这是因为有默认参数：dirichlet。所以，除非不是dirichlet条件，可以不需要明指出具体的BcType。BcType是一个枚举类型：

1
2
3
4
enum  	BcType { 
  dirichlet, neumann, freeslip, density, 
  outflow, normalOutflow 
}

如果lattice的外框中的某条边作为边界条件，那么代码可以简化如下：

1
boundaryCondition - > setVelocityConditionOnBlockBoundaries （lattice ， locationOfBoundaryNodes ）;

这里的locationOfBoundaryNodes 就是一条边。而如果lattice的外框都是边界条件，并且是同样的，那么代码可以写成：

1
boundaryCondition - > setVelocityConditionOnBlockBoundaries （lattice ）;

现在，假设长为nx，宽慰ny的block-lattice中的上壁和下壁（包括角）实现自由滑动条件（切向速度分量为零梯度），左壁为Dirichlet条件，右壁为流出条件（所有速度分量的零梯度）。则代码实现为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
Box2D inlet(0, 0, 2, ny-2);//入口
Box2D outlet(nx-1, nx-1, 2, ny-2);//出口
Box2D bottomWall(0, nx-1, 0, 0);//上壁面
Box2D topWall(0, nx-1, ny-1, ny-1);//下壁面

//默认Dirichlet条件
boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, inlet );

boundaryCondition->setVelocityConditionOnBlockBoundaries (
                               lattice, outlet, boundary::outflow );
boundaryCondition->setVelocityConditionOnBlockBoundaries (
                               lattice, bottomWall, boundary::freeslip );
boundaryCondition->setVelocityConditionOnBlockBoundaries (
                               lattice, topWall, boundary::freeslip );

如果边界条件不在lattice的外边，那么需要增加一个BoxXD来进行限定。就像那种最复杂的情况。

1.3 BcType

前面提到，BcType是一个命名空间plb:boundary下面的枚举类型。下面是各种不同边界条件的物理含义。

velocity boundaries Note
boundary::dirichlet (default value) Dirichlet condition: imposed velocity value.
boundary::outflow or boundary::neumann Zero-gradient for all velocity components.
boundary::freeslip Zero-gradient for tangential velocity components, zero value for the others.
boundary::normalOutflow Zero-gradient for normal velocity components, zero value for the others.
pressure boundaries Note
boundary::dirichlet (default value) Dirichlet condition: imposed pressure value, zero value for the tangential velocity components.
boundary::neumann Zero-gradient for the pressure, zero value for the tangential velocity components.

值得注意的是：边界条件不能被覆盖

无法覆盖边界的类型。一旦边界条件被定义为是Dirichlet速度边界条件，它就会保持Dirichlet速度边界条件。将其重新定义为别的边界条件是没有用的。因此，必须仔细定义边界条件，并且分段定义，如上例所示。

不过定义的边界条件的数值是可以随时改变的。

1.4 在Dirichlet边界上设置速度或压力值

使用函数setBoundaryVelocity施加恒定速度值。以下命令在2D域的所有节点上将速度值设置为零，这些节点之前已被定义为Dirichlet速度节点：

1
setBoundaryVelocity(lattice, lattice.getBoundingBox(), Array2>(0.,0.) );

然而在其他的节点，这个命令没有作用。

如果是要设定压力边界条件，可以采用以下命令：

1
setBoundaryDensity(lattice, lattice.getBoundingBox(), 1. );

在LBM中压力和密度是直接挂钩的，因此设置压力其实就是设置密度。

时变的，与空间相关的速度和压力边界条件也可以实现。具体参考案例：examples/showCases/poiseuille/poiseuille.cpp. 即伯肃叶流动。

值得注意的是，setBoundaryVelocity属于plb下定义的公共函数，不属于任何类。

1.5 内部的和外部的边界条件

有时，边界不位于规则的块状区域的表面上。在这种情况下，setVelocityConditionOnBlockBoundaries等函数是无用的，并且必须手动构造边界。下面是一个例子，介绍如何手动构建沿顶壁的自由滑动条件，包括角节点：

1
2
3
4
5
6
7
8
9
10
// The automatic creation of a top-wall ...
Box2D topWall(0, nx-1, ny-1, ny-1);
boundaryCondition->setVelocityConditionOnBlockBoundaries (
                               topWall, lattice, boundary::freeslip );

// can be replaced by a manual construction as follows:
Box2D topLid(1, nx-2, ny-1, ny-1);
boundaryCondition->addVelocityBoundary1P( topLid, lattice, boundary::freeslip );
boundaryCondition->addExternalVelocityCornerNP( 0, ny-1, lattice, boundary::freeslip );
boundaryCondition->addExternalVelocityCornerPP( nx-1, ny-1, lattice, boundary::freeslip );

第一组代码，是自己定义了一个topWall。可以看到这其实就是一个2Dbox的顶线。

第二组代码，用的不是setVelocityConditionOnBlockBoundaries这个方法，而是里面对顶线和角点分别进行边界条件设定的方法。topLid也是一个边，只是不包含两边的角点。

根据角落是凸面还是凹面，区分外角和内角。如下图所示：

内角点和外角点
内角点和外角点

对于addVelocityBoundary1P后面的1P，其实是代表了边的方向，你可以理解成边的法向量。这样的后缀还有0P，0N，1N。加上1P一共是四个。

  • 代码1P表示：“墙法线指向正y方向”
  • 代码0N标记为“负x方向”
  • 代码1N表示正x方向
  • 代码0P表示y负方向
  • 一般边界的法向量都是指向流场的外面。这和边界条件设定为入口还是出口没有任何关系，这些后缀代表的就是纯几何上的法向量。

在3D中也有类似的东西。

在后面的文章中我会对这些做一个更为细致的介绍。

1.6 总结

总结一下设置速度和压力边界条件，以及outflow这种边界条件的基本方法。

  • 首先，要定义一个接口，也就是定义一个OnLatticeBoundaryConditionXD类。采用不同的边界类型算法，采用不同的初始化方法，比如，要让边界采用Zou/H算法来定义计算边界条件，就用createZouHeBoundaryConditionXD来进行初始化；
  • 得到了这个接口之后，调用创建好的这个接口下的成员函数：setVelocityConditionOnBlockBoundaries或者setPressureConditionOnBlockBoundaries，来指定某个边界的类型。其中，指定的边界可能需要某种方式去进行构建；
  • 采用plb:setBoundaryVelocity()这样的方法来设定边界上的具体数值。

至此，完成了边界条件的设定。

2 周期性边界

这部分因为没有详细了解过，所以先挖坑吧。

文档

3 反弹边界条件

首先，反弹边界条件主要用于对wall实现不可滑移边界条件的实现。

Palabos里面的反弹边界条件的执行方法分为2种：Bounce-Back 1和Bounce-Back 2。

从文档的解读可以知道：Bounce-Back 1实际上对应的是 halfway bounce-back。Bounce-Back 2对应的是 fullway bounce-back。

需要注意的是，不管是Bounce-Back 1还是Bounce-Back 2，都模拟的是边界位于固体节点和流体节点的中间位置。如下图所示：

Bounce-Back 1 VS Bounce-Back 2
Bounce-Back 1 VS Bounce-Back 2

其中图a表示的是全步长反弹，也就是Bounce-Back 2，图b表示的是半步长反弹，也就是Bounce-Back 1。下面来看看这两者到底有什么区别。

注意：图片来自于文献[1]。

3.1 Bounce-Back 1 VS Bounce-Back 2

  • Bounce-Back 1是一种半步长反弹（HBC）。

官方文档里指出，这种算法是在边界节点邻近的流体节点进行下一次的碰撞之前，该点的分布函数取上一次碰撞结束后的相反方向的分布函数值，用公式来表达，就是：

$f_{3}\left(x_{N}, t+\Delta t\right)=f_{1}^{\star}\left(x_{N}, t\right)$

其中带星号的分布函数是代表着进行碰撞后得到的值，而不是迁移后得到的值。此外，f1和f3的方向是相反的。

由这个式子及上门的图b可以看到，HBC里面，边界节点是不进行碰撞和迁移的，你可以认为它不对计算起任何作用，也就是说边界节点不用存储任何变量或数值。迁移只在流体节点进行，标准碰撞也只发生在流体节点。此外，靠近边界的这个流体节点迁移不到边界，而是直接反弹了回来。也就是说本来应该有边界迁移过来出来的分布函数，直接在该流体节点反弹了回来。这个分布函数作为下一次碰撞的初始条件。

也可以说，这种算法，反弹是发生在流体节点和固体节点中间的位置。改变了邻近边界节点的流体节点的迁移步骤

  • Bounce-Back 2是一种全步长反弹（FBC）。

文档中写到，如果将这种反弹的发生地点指定在固体节点中（文档是说的新的节点，实际上就是固体节点），这样，这些节点亦不发生碰撞，而是将碰撞过程代替为反弹过程。注意，这和HBC有一个很大的区别，要在固体边界上执行反弹。这些节点被视为纯算法节点，不具备真实的物理意义，也就是根据这些固体节点计算得到的宏观量（密度，动量）是不准确的。这是由于固体节点只起到一个反弹边界的作用。

用公式表达那么就是如下：

  1. 碰撞完成之后，正常迁移至边界节点：

$f_{1}\left(x_{N+1}, t\right)=f_{1}^{\star}\left(x_{N}, t-\Delta t\right)$

  1. 边界节点执行反弹规则，实际上是用反弹替换了碰撞：

$f_{3}^{\star}\left(x_{N+1}, t\right)=f_{1}\left(x_{N+1}, t\right)$

  1. 反弹完成之后，将该分布函数迁移至流体，作为下一次碰撞的初始值：

$f_{3}\left(x_{N}, t+\Delta t\right)=f_{3}^{\star}\left(x_{N+1}, t\right)$

从这些式子，并结合上图可以看到，边界节点不执行碰撞，而是利用反弹来对其进行替换。也就是说，反弹改变了边界节点的碰撞环节，但是并不改变它的迁移环节。

总结一下：

  • Bounce-Back 2是全步长反弹，其对时间精度为1阶。Bounce-Back 1是半步长反弹，对时间精度为2阶；
  • Bounce-Back 2需要固体边界节点，用于存储执行边界反弹和迁移步骤所需要的分布函数。而Bounce-Back 1不需要边界节点存储任何变量；
  • Bounce-Back 2在边界节点上执行正常的迁移，但是碰撞环节被反弹环节所取代，至于靠近边界节点的流体节点，不改变任何运作方式；而Bounce-Back 1在靠近边界节点的流体节点上执行正常的碰撞，但是迁移环节有所改变，本应由边界传递而来的分布函数在这里取值为与其相反方向的碰撞后的分布函数。注意，不是进过迁移后的分布函数。这也就意味着，在执行碰撞之后需要对靠近边界节点的流体节点在那些本应该从边界获得的分布函数进行先赋值，以免迁移过程中发生数值覆盖。

事实上，关于这两种碰撞方法的问题，困扰了我将近五天。

这个问题的复杂在于并不知道Bounce-Back 1和Bounce-Back 2到底代表的什么算法。实际上反弹规则还有很多种，分类方法也不一样，除了全步长和半步长反弹外，根据边界位于节点的哪个位置，有内推移法和外推移法，还有边界就位于边界节点的方法。此外，根据边界节点是否执行标准碰撞过程，还分为标准反弹和修正反弹格式。最终在Krüger, T.K.H.K.等于2017年的出版的书籍上获得了关于半步长和全步长的最为清晰的解答。也为此确定了Bounce-Back 1和Bounce-Back 2的运行机理。

关于这本书，大家可以在Google找到并且看到，由于版权原因，这里就不进行分享了，如果找不到可以取置顶文章里找的本博客的群，在群里询问群主。这本书写的很好，很有深度。广度也不错，我认为是这么多年来看到写的最好的一本关于LBM的教材，教材里面的内容也比较新，能跟得上时代，推荐初学者。

3.2 根据域类定义边界域

为了描述一个反弹域，一般是需要确定这个反弹域对应的节点。推荐采用继承一个DomainFunctionalXD类来进行。DomainFunctionalXD是一个矩形或者是立方体的域，其内包含了私有成员变量_iX_和_iX_（或者还有_iX_）。文档中的例子是继承DomainFunctional2D类，创建了一个圆形类，并且重写了operator方法，这个operator方法的作用是获得边界所在的边界节点。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
template<typename T>
class CylinderShapeDomain2D : public DomainFunctional2D {
public:
    CylinderShapeDomain2D(plint cx_, plint cy_, plint radius)
        : cx(cx_),
          cy(cy_),
          radiusSqr(util::sqr(radius))
    { }
    // The function-call operator is overridden to specify the location
    //   of bounce-back nodes.
    virtual bool operator() (plint iX, plint iY) const {
        return util::sqr(iX-cx) + util::sqr(iY-cy) <= radiussqr;< span>
    }
    virtual CylinderShapeDomain2D* clone() const {
        return new CylinderShapeDomain2D(*this);
    }
private:
    plint cx, cy;
    plint radiusSqr;
};

对这样的一个区域定义动态的语句为：

1
2
3
defineDynamics( lattice, lattice.getBoundingBox(),
                new CylinderShapeDomain2D(cx, cy, radius),
                new BounceBack );

这里调用了一个函数：defineDynamics。属于plb。第一个参数是各自域，第二个是这个各自域的所构成的Box2D，第三个是一个域类，也就是上面定义的CylinderShapeDomain2D类，第四个是BounceBack动态。反弹节点位于这个lattice.getBoundingBox()和定义的CylinderShapeDomain2D相交的单元格。

如果边界域被指定为更简单域集合的并集，则可以在每个更简单的域上迭代使用defineDynamics()，文档中写到了一共有15种形式的defineDynamics()。如果域是更简单域的交集或任何其他几何操作，则需要在域功能的定义内使用布尔运算符。

3.3 基于布尔运算符的反弹域的定义

假设域的几何形状是从外部源规定的。例如，当模拟多孔介质并拥有来自所讨论的多孔材料的扫描数据时。这种最简单的方法是将此数据存储为ASCII文件，该文件将流体节点与实体节点区分为0和1，用空格，制表符或换行符分隔。该文件表示常规数组的内容，仅存储原始数据，不存储矩阵维度的信息。数据布局必须与Palabos中的相同：z坐标的增量表示文件中的连续进度（在3D中），然后是y坐标，然后是x坐标。只需将该数据从文件刷新到标量字段即可。以下代码取自示例examples / codesByTopic / io / loadGeometry.cpp：

1
2
3
MultiScalarField2D boolMask(nx, ny);
plb_ifstream ifile("geometry.dat");
ifile >> boolMask;

下一步，使用函数defineDynamics的一个版本实例化反弹节点：

1
defineDynamics （lattice ， boolMask ， new  BounceBack < T ，DESCRIPTOR > ， true ）;

对于上述的操作，并没有自动错误检测。所有使用这个功能，使用者必须要确保标量字段的尺寸与几何文件中的数据大小相对应。

4 问题

  • 问题1：

这里提到了利用反弹边界条件来实现wall的定义，通常wall的条件表现为边界速度为0（静止的wall），那么是否能通过设定Dirichlet边界条件，并使其速度值=0来实现固态的边界条件呢？

与此同时，我还看到，在LBM中，运动的wall，通常还是用反弹方式来实现：

并没有用Dirichlet。

关于这一点，我认为还是可以通过利用Dirichlet来实现。因为本质上，不管是运动的wall还是静止的wall还是Dirichlet条件。

可以在案例中进行计算试试。

  • 问题2：

在定义反弹边界条件的时候，直接指定了反弹边界节点处应用BounceBack 这种动态行为。而在前面的边界条件的实施的时候，并没有看到指定动态行为的代码。是否这个工作包含在内部的代码以内，还是说通过数据处理器来实现了?

这个目前没有答案，需要进一步去看代码。

话说回来，文章开始写是5月22号，到现在才发布。足足写了5天，还遗留下很多疑惑。果然边界条件还是不好弄啊。

对上述两个问题有清楚的大佬，或者对本文有疑问的同学，均可以在下方留言。

另外，欢迎大家加群一起讨论。

[1]: Krüger, T.K.H.K., The Lattice Boltzmann Method: Principles and Practice. 2017.