gulie的个人博客分享 http://blog.sciencenet.cn/u/gulie

博文

[转载]LBM:PalaBos 几何与边界的处理

已有 6376 次阅读 2011-12-1 21:56 |个人分类:知识的海洋|系统分类:科研笔记|文章来源:转载

在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>());



http://blog.sciencenet.cn/blog-250582-513842.html

上一篇:在Word2007中添加公式编辑器图标
下一篇:[转载]LBM:PalaBos进行程序编制

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备14006957 )

GMT+8, 2019-11-14 08:36

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部