LBM: PalaBos几何及边界的定义

在PalaBos中有很多中边界条件的定义方法，其可针对平直的边界、网格相关的边界（也就是我们常说的计算是在真是的物理边界上执行的）。 一般情形下，有这么两种选择：基于Bounce-Back的阶梯逼近；高阶插值的曲线或曲面边界逼近。目前， 在1.0版本中，曲线或者曲面边界已经实现，可参考动脉瘤的例子（aneurysm）。

网格相关的边界条件：

1. 所使用的类：OnLatticeBoundaryConditionXD （XD，中的X选取为2或3，针对于2维或3维）。

在边界的条件的执行中，有些是局部执行的，有些是非局部。局部执行方法在dynamics中直接实现，而非局部需要借助外部的数据处理方法来实现，比方：集成过程。PalaBos给我们提供一些接口，来帮助我们执行可能的方法。OnLatticeBoundaryConditionXD负责初始化dynamics对象，并添加数据处理机制。例如：Zou/He边界条件
OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition =
        createZouHeBoundaryCondition3D<T,DESCRIPTOR>();
这里给出四种常用的边界处理方法，并给出可能的优点或缺点的说明：
Regularized BC     createLocalBoundaryConditionXD
缺点：在处理outflow边界时，会引起边界不恰当的振荡。
Skordos BC     createInterpBoundaryConditionXD
优点：比较稳定
Zou/He BC     createZouHeBoundaryConditionXD
有点：处理入口和出口比好
Inamuro BC     createInamuroBoundaryConditionXD

借助上面的边界条件，在PalaBos中我们可以执行Dirichlet边界条件、压力边界条件。此外，也可以执行Neumann边界条件。

这需要提提边界条件处理的一些技巧，原因在于，在Palabos中，对象需要知道给定边是否是平直的或者是否是几何的角点（2D），平直的壁面或者边或者角点（3D），同时需要知道边界的方向。OnLatticeBoundaryConditionXD提供方法setVelocityConditionOnBlockBoundaries和setPressureConditionOnBlockBoundaries设定落在矩形块面上的边界。意思就是说，如果所有的边界点都落在定义的Box上，这些方法可以直接使用，如果不是，就需要重新定义。

2. 简单的边界设置方法给出如下：
Box2D boundaryBox(x0,x1, y0,y1);
boundaryCondition->setVelocityConditionOnBlockBoundaries (
                           lattice, boundaryBox, locationOfBoundaryNodes );

如果所有的边界点都落在外部的矩形或者立方体的边上，上述定义可简化为：
boundaryCondition->setVelocityConditionOnBlockBoundaries(lattice, locationOfBoundaryNodes);
如果所有的矩形或者立方体的格点都是边界点，进一步可以简化为
boundaryCondition->setVelocityConditionOnBlockBoundaries(lattice);

现在给出一个2D Channel的定义，首先定义边界相关的Box：
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);
设定边界类型
boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, inlet );
boundaryCondition->setVelocityConditionOnBlockBoundaries (
                               lattice, outlet, boundary::outflow );
boundaryCondition->setVelocityConditionOnBlockBoundaries (
                               lattice, bottomWall, boundary::freeslip );
boundaryCondition->setVelocityConditionOnBlockBoundaries (
                               lattice, topWall, boundary::freeslip );

boundary::outflow指0速度梯度边界条件；boundary::freeslip滑移边界条件。在上面的定义中，所有的边界点都落在外部的矩形上，所以不需要特殊的处理。


boundary::dirichlet (default value)     Dirichlet边界条件
boundary::outflow 或 boundary::neumann     指0速度梯度边界条件
boundary::freeslip                     速度在切向的投影为0
boundary::normalOutflow             法向0速度梯度边界条件

3. 压力边界条件安的选择
boundary::dirichlet (default value)     Dirichlet设定压力的取值，切向速度为0。
boundary::neumann                     0压力梯度，切向速度为0。

边界类型只能被定义一次，在运行过程中，不能再改变边界类型，如果进行重新定义，其并不起作用。这个在定义中一定要注意。尽管如此，定义的边界上的取值是可以改变的，比方可以定义时间依赖的边界取值。可以参见例子：examples/showCases/poiseuille/poiseuille.cpp。

外部和内部边界条件

有时，边界点没有落在正规的边界面上，在这种情形下，函数setVelocityConditionOnBlockBoundaries将不会起作用，此时边界必须手动重建。如何手动定义，这里以滑移边界条件为例，在上边界或者下边界定义滑移边界条件时:

1. 最直接的定义是：
Box2D topWall(0, nx-1, ny-1, ny-1);
boundaryCondition->setVelocityConditionOnBlockBoundaries (
                               topWall, lattice, boundary::freeslip );

也可以采用如下方法：
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 );

解释：
addVelocityBoundary1P 1表示沿着y，P表示沿着y的正向；同理，0P表示x的负方向；addExternalVelocityCornerNP 表示上边界右端点, NP便是方向：N 表示x的负方向；P表示y的正向。addExternalVelocityCornerPP的处理方法类似。 这都是针对2D情形。关凸或者凹的角点，这里就不再做过多的解释。

2. 对于3D平面采用使用如现的方法处理:

addVelocityBoundaryDO 这里 D可以是 0（x）, （1）y, （2）z； 后面字符O可以取P和N表示沿着D的正向还是负向。 对于内部和外部边的定义，示例如下:
addInternalVelocityEdge0NP （凹入）
addExternalVelocityEdge1PN （凸出）
第一个参数0或1暗示这条边平行于那个坐标轴,后面的两个N和P表示这跳边所在面的方向。0NP 表示 沿着x,其一个面方向为y的负向，另外一个面为z的正向。3D角点的处理方法：
addInternalVelocityCornerPNP  
addExternalVelocityCornerNNN
＊其解释于2D相同

周期边界条件：

1. 处理方式如下：

lattice.periodicity.toggle(0, true);      表示所有的量周期性沿着x
scalarField.periodicity.toggle(0, true);  表示标量沿着x周期

如下可以所有边界都是周期的：
lattice.periodicity.toggleAll(true);
对于周期性边界条件在定了块以后，就应该定义，这样PalaBos就可以添加执行相应的数据处理函数。

弹性式边界条件 Bounce-back （局部执行）

1. 在这类边界条件的处理中，边界点不再被看作是需要进行统计格点，仅仅就是算法层次，计算精度于网格尺度同阶。这种边界适合复杂几何体中流动的计算。下面给出一个示例：
    
定义几何体（此函数必须继承于DomainFunctionalXD）：
template<typename T>
class CylinderShapeDomain2D : public DomainFunctional2D {
public:
    CylinderShapeDomain2D(plint cx_, plint cy_, plint radius)
        : cx(cx_),
          cy(cy_),
          radiusSqr(util::sqr(radius))
    { }
    virtual bool operator() (plint iX, plint iY) const {
        return util::sqr(iX-cx) + util::sqr(iY-cy) <= radiusSqr;
    }
    virtual CylinderShapeDomain2D<T>* clone() const {
        return new CylinderShapeDomain2D<T>(*this);
    }
private:
    plint cx, cy;
    plint radiusSqr;
};

需要对建立的几何题定义dynamics:

defineDynamics( lattice, lattice.getBoundingBox(),
                new CylinderShapeDomain2D<T>(cx, cy, radius),
                new BounceBack<T,DESCRIPTOR> );

如果是多个几何体，有相交的情形，可以采用定义bool类型变量的方法，这方法很适合模拟多孔介质内的流动。示例如下：
MultiScalarField2D<T> boolMask(nx, ny);
plb_ifstream ifile("geometry.dat");
ifile >> boolMask;
defineDynamics(lattice, boolMask, new BounceBack<T,DESCRIPTOR>, true);

2. 当然采用Bounce-back定义边界条件也不失为一种好方法：如下：

Box2D top(   0,    nx-1, ny-1, ny-1);
Box2D bottom(0,    nx-1, 0,    0);
Box2D left(  0,    0,    1,    ny-2);
Box2D right( nx-1, nx-1, 1,    ny-2);
 
defineDynamics(lattice, top, new BounceBack<T,DESCRIPTOR>());
defineDynamics(lattice, bottom, new BounceBack<T,DESCRIPTOR>());
defineDynamics(lattice, right, new BounceBack<T,DESCRIPTOR>());
defineDynamics(lattice, left, new BounceBack<T,DESCRIPTOR>());

如有问题欢迎讨论……