Palabos User Guide中文解读 | 第九章 | 定义边界条件

最新推荐文章于 2023-04-08 12:24:45 发布

方雨岚(Yulan Fang)

最新推荐文章于 2023-04-08 12:24:45 发布

阅读量2.8k

点赞数 4

分类专栏： Palabos入门教程&初学文章标签：流体力学

本文链接：https://blog.csdn.net/qq_40259141/article/details/104436838

版权

Palabos入门教程&初学专栏收录该内容

52 篇文章 192 订阅

订阅专栏

作者的话：本人在学习palabos时，发现国内中文资料甚少，恰好网上可以直接搜到palabos user guide这种英文资料，加之时间充裕，便打算开始翻译，翻了一节后发现这可能算侵权，就比较伤脑筋，突然想到自己写中文解读即可，便有了下面的博客。

Palabos User Guide
Release 1.0 r1
Copyright © 2019 University of Geneva
Jul 05, 2019

Chapter Nine

定义边界条件

9.1 概述

大致呢，palabos库提供了一些不同的边界条件，可供平直对齐的或者弯曲的边界使用。

重要的一点：弯曲边界只能在1.0版本中使用，包括并行的体素化和I/O扩展。目前弯曲边界的记录不是很详细。
可在这里查找：showCases/aneurysm

9.2 Grid-aligned boundaries平直对齐的边界

9.2.1 OnLatticeBoundaryConditionXD 类

有些边界条件是局部的，所以是dynamics类的形式被运行，有些不是局部的，则被数据处理器运行，有些则是半局部半非局部的。

Palabos里有OnLatticeBoundaryConditionXD，用于初始化dynamics对象，也可添加数据处理器，看你选哪个边界条件。如下例子是选择Zou/He边界条件。

OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition =
createZouHeBoundaryCondition3D<T,DESCRIPTOR>();

除此以外还有4个其他选项：（palabos.org-bc 这个网址里有更多细节）

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

这些边界都可以使用velocity Dirichlet boundaries和压力边界。

9.2.2 规则域（ regular block domain）的边界

2D情况下，Palabos需要知道给定的边界是平的还是角落，3D情况下Palabos需要知道边界是平的墙还是角落的边界。

为了简化上述的需求，所有的OnLatticeBoundaryConditionXD对象都含有一个方法：setVelocityConditionOnBlockBoundaries和setPressureConditionOnBlockBoundaries，
来设置四边形域的边界。

如果boundaryBox的2D域边界上所有点都应用了Dirichlet条件，就用下面的代码：

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

第一行表示边界位置，第二行规定了哪些点上初始化速度。

如果边界上所有点都在格子的反弹域之外，上面的代码需要变成：

boundaryCondition->setVelocityConditionOnBlockBoundaries(lattice, locationOfBoundaryNodes);

如果说，所有反弹域之外的节点都需要设置速度条件，即需：

boundaryCondition->setVelocityConditionOnBlockBoundaries(lattice);

好的，到这里我们假定有一个长nx，宽ny的流域，使用自由滑移条件，左侧边界使用Dirichlet条件，右侧边界为出口，代码如下：

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

需要牢记的是，如果边界不确定在lattice之外，必须要使用一个额外的 Box3D 语句来表示边界的位置。毕竟你需要告诉Palabos流体的方向和边界的情况。

9.2.3 零梯度条件（Zero-gradient conditions）

在boundary::outflow这块地方，可以修改成其他参数，如下面所示。

boundary::dirichlet (默认值)
Dirichlet condition: imposed velocity value.

boundary::outflow 或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.

下面是压力边界。

boundary::dirichlet(默认值)
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速度节点，如果你重新定义它为压力边界点，是没有用的，由此在定义边界的时候你一定要当心。
当然边界节点的数值是可以按你需求随便修改的。

9.2.4 设置Dirichlet边界的速度或压强

setBoundaryVelocity 函数可以赋值一个常数速度。下面这行代码的意思是2D域下将所有之前已经设为Dirichlet速度节点的点的速度设置为0

setBoundaryVelocity(lattice, lattice.getBoundingBox(), Array<T,2>(0.,0.) );

上行代码对所有别的节点无影响，若是设置压力值（等同来讲，密度值），使用 _setBoundaryDensity_命令。

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

A non-constant, space dependent velocity resp. pressure profile can be defined by replacing the velocity resp. pressure argument by either a function or by a function object (an instance of a class which overrides the function call operator)which yields a value of the velocity resp. density as a function of space position. An example is provided in the file examples/showCases/poiseuille/poiseuille.cpp.

9.2.5 内部和外部边界

setVelocityConditionOnBlockBoundaries 仅仅在设置规则域最外表面边界条件时用得到。

下面的例子是如何手动建立一个自由滑移条件的顶层边界（包括拐角点）。

// The automatic creation of a top-wall …自动生成从(0, ny-1)到(nx-1, ny-1)的自动滑移条件顶层边界
Box2D topWall(0, nx-1, ny-1, ny-1);
boundaryCondition->setVelocityConditionOnBlockBoundaries (
lattice, topWall, boundary::freeslip );

// … can be replaced by a manual construction as follows:手动生成自动滑移条件顶层边界
Box2D topLid(1, nx-2, ny-1, ny-1); //从(1, ny-1)到(nx-2, ny-1)，左右各空了一个(0, ny-1)，(nx-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 );

其中1P, NP, PP表示边界面法线的方向。值得注意的是，边界法线仅仅是几何因素，而不会取决于边界是进出口，在任何情况下，避免法线都会指向流体之外，朝向非流体区域。

对于平直的边界，1P表示y轴正向，1N表示y轴负向，0P表示x轴正向，0N表示x轴负向。
拐角处，NP表示X负向与Y正向（即法线指向左上）。

若想判断拐角点为外拐角点还是内拐角点，只需要看这个直角的尖头指内还是指外，下图有一个指内的内拐角点，五个指外的外拐角点。

3D情况下，你可以使用下面代码：

addVelocityBoundaryDO

D的取值有0,1，2，分别对应x,y,z。
O取值有N和P，对应负和正。

而在面面交界处的例子如下：

addInternalVelocityEdge0NP
addExternalVelocityEdge1PN

末尾三个数对应的xyz是周期性的，取决于数字012对应xyz开始。
0NP：x面，负向y值，正向z值。
1PN：y面，正向z值，负向x值，

3D拐角点例子：

addInternalVelocityCornerPNP
addExternalVelocityCornerNNN

9.3 循环边界

每当Palabos添加或者执行数据处理器的时候，都需要知道边界是否循环，由此每当你创建一个新的区域之后，尽快定义这块区域是否边界循环。

循环边界与之前讨论的边界条件不同。它只能作用于对立的，atomic-block或multi-block的最外边界。只要对立边界也是流体节点，它就可以实现循环边界。

在atomic-block情况下，它只是一个流动的操作，不影响标量场和张量场。
但在multi-block情况下，它就影响了，在这种情况下设置循环边界则意味着你定义了一个非局部的数据处理器，这个数据处理器会访问周边最近的节点，这些邻居节点的值会由multi-block的最外层边界的周期性决定。

默认情况下，所有的block都是不循环的，默认情况下最外边界节点的行为是一个版本的反弹算法。
若想打开沿着x轴的循环功能，添加如下代码：

lattice.periodicity().toggle(0, true); // Periodic block-lattice.
scalarField.periodicity().toggle(0, true); // Periodic scalar-field.

你也可以用一行代码让所有边界都循环：

lattice.periodicity().toggleAll(true);

最后，Neumann类型的边界（出口，自由滑移，绝热边界等）万万不可被设为循环边界，如果你非要这么做，程序会崩。

9.4 反弹边界

在Palabos里，所有默认的边界都是反弹边界。可添加dynamics对象BounceBack<T,Descriptor>到选定的格子中。

在算法实现中，我们假定t-1时刻未知粒子群获得对立粒子群碰撞后的值，这种版本我们称之为反弹-1.
若是从t-2时刻中获得值，则成为反弹-2。这两种反弹的版本都会在两个节点之间生成一个无滑移边界，在反弹1版本，与反弹流体格子只有半个格子长的宽度，相比反弹2版本，这只有一半的路程。这是种牺牲计算精度的方法，反弹2无需考虑边界朝向问题，也无需知道拐角点，所以在设置流体问题的时候，最好先用反弹-2的方式先作尝试。
在特别复杂的流场环境中，如多孔介质，反弹边界往往比其他方法更好。

9.4.1 反弹区域的分析解释

Palabos中的第零规则是：不要写终端用户处写循环。

将反弹dynamics添加到lattice格子里，一般用数据处理器，或者defineDynamics函数，在例子：examples/codesByTopic/bounceBack/instantiateCylinder中可以发现。

首先，iX和iY作为反弹点的位置坐标的函数，被写到一个类里来定义2D圆柱域。

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