Palabos Tutorial 阅读笔记 | 2.5 边界条件 Boundary conditions

前言

本人是正在借助AI学习Lattice Boltzmann method 和 Palabos的小白。在此将理解的Palabos Tutorial的内容记录下来留做备忘，并与大家分享，希望得到各位大神指点。欢迎随时指出其中的错误，也欢迎有相同目标的同学一起交流。

2.5.0 介绍

在这里这里我们第一次通过数值模拟还原了一个著名物理情形：Poiseuille流。该流动发生在一个有不可滑动壁的通道内，并且流速与通道壁平行。该流动的解析解由沿通道壁平行的速度分量的抛物线状轮廓表示。在模拟中，解析解用于在域的入口和出口上采取 Dirichlet边界条件。作为初始条件，采用了零速度场，并开始模拟直至在域内的每个点上收敛以得到Poiseuille解。

• Dirichlet边界条件（固定边界条件）是一种边界条件，用于描述解决偏微分方程问题时，通过指定问题解的边界值来限制解的行为。在Dirichlet边界条件下，已知在边界上的解的值。例如，对于一个二维的热传导问题，在边界上指定温度值就是一个Dirichlet边界条件。这个条件通常表示为 u(x,y)=f(x,y)，其中 u 是问题的解，f(x,y) 是在边界上给定的函数。

请注意，模拟最初在初始条件中的零速度和边界条件中的有限速度之间存在不连续性。如果在调整参数时模拟变得数值不稳定，请尝试增加晶格分辨率或减小晶格单位中的速度。

2.5.1 tutorial_1_5.cpp 代码思路

• 定义流体问题的边界条件和初始条件
• 设置模拟参数
• 通过Lattice Boltzmann迭代求解
• 然后进行数据分析和输出结果。

代码中的主要函数包括：

• poiseuilleVelocity：计算泊肖流模型的速度分布。
• PoiseuilleVelocity类：一个仿函数，用于初始化速度和密度的边界条件。
• channelSetup：设置边界条件和初始条件。
• writeGifs：保存流体流动的动画。
• main函数：程序的入口，定义了模拟的总时间和图像保存的时间间隔，创建Lattice Boltzmann模拟器和边界条件，执行模拟和保存结果。

总体来说，这段代码演示了如何使用Lattice Boltzmann方法模拟Poiseuille流，并将结果输出为动画和数据文件。

2.5.2 tutorial_1_5.cpp 代码结构

这段C++代码是基于Palabos库的一个简单的二维泊肃叶流（Poiseuille flow）模拟。

2.5.2.1 包含头文件

包括C++标准库中的iomanip和iostream，以及Palabos库中的相关头文件。

2.5.2.2 命名空间

将全局命名空间简化为“plb”和“std”。

2.5.2.3 类型定义和描述符

定义数据类型T为double，并设置DESCRIPTOR为Palabos库中的D2Q9Descriptor。

2.5.2.4 函数和类定义

• poiseuilleVelocity：计算泊肃叶速度分布。

• PoiseuilleVelocity：一个模板类，用于初始化速度边界条件和初始条件。

• channelSetup：设置模拟通道（设置边界条件、指定边界速度、创建初始条件）。

• writeGifs：将模拟结果保存为GIF图像。

2.5.2.5 主函数 main

• 初始化Palabos库。
• 设置输出目录。
• 定义参数、模拟时间等。
• 创建一个多块晶格，并为其分配一个动态模型（BGK）。
• 创建一个边界条件对象。
• 调用channelSetup设置模拟通道。
• 在主循环中进行格子玻尔兹曼迭代，并在特定时间间隔保存GIF图像。
• 删除边界条件对象。

2.5.2.6 说明

在这个模拟中，通道的宽度和高度都是有限的。速度边界条件设置为泊肃叶速度分布，初始条件根据泊肃叶速度分布设置。模拟的主循环执行格子玻尔兹曼迭代，并在特定时间间隔保存速度场的GIF图像。

这个例子展示了如何使用Palabos库设置一个简单的二维泊肃叶流模拟，并保存速度场的可视化结果。

2.5.3 分行代码详解

2.5.3.1 包含头文件

这段C++代码包含了四个头文件的引用：

• <iomanip>：C++ 标准库头文件，包含了一些与格式化 I/O 相关的函数和类型，例如 setprecision() 和 setw()。
• <iostream>：C++ 标准库头文件，包含了输入输出流的定义以及与其相关的类、函数、常量等。
• <vector>：C++ 标准库头文件，包含了 STL 中的 vector 容器的定义以及与其相关的类、函数、常量等。
• "palabos2D.h" 和 "palabos2D.hh"：这两个头文件是 Palabos 库中用于二维流体模拟的头文件，包含了对应的类和函数的定义。

其中，"palabos2D.h" 包含了一些较为基础的类的定义，例如 MultiBlockLattice2D 和 OnLatticeBoundaryCondition2D；而 "palabos2D.hh" 包含了一些更高级的功能的定义，例如图像的读写、流场数据的可视化等。

通过引用这些头文件，可以使用相应的函数和类，方便进行二维流体模拟。

2.5.3.2 命名空间

这一段C++代码使用了C++的命名空间（namespace）功能，用来简化代码中的名称冲突问题。使用命名空间后，可以直接使用命名空间中的成员函数、变量等，而不需要再使用作用域运算符（::）来指明所属的命名空间。 具体来说，这段代码使用了plb和std两个命名空间。

• plb命名空间是Palabos库的命名空间，包含了库中的各种类和函数；
• std命名空间是C++标准库的命名空间，包含了C++标准库中的各种类和函数。

因此，在代码中使用了plb命名空间和std命名空间后，就可以直接使用这些命名空间中的类和函数，而不需要再使用作用域运算符来指明所属的命名空间。

2.5.3.3 类型定义和描述符

这段代码定义了两个常量 T 和 DESCRIPTOR。

• T 是一个类型别名（typedef），它表示使用双精度浮点数来存储数值。
• DESCRIPTOR 是一个宏定义，它被定义为 plb::descriptors::D2Q9Descriptor，是一个类型别名，表示使用 Lattice Boltzmann 方法中的 D2Q9 模型进行流场模拟。在 Palabos 库中，描述子（descriptor）定义了模型的性质，例如模型中的速度集合、权重和速度方向等。

2.5.3.4 函数和类定义

这段代码定义了一个名为poiseuilleVelocity的函数，它的作用是计算在泊肃叶流中，某个位置的速度值。

具体来说，这个函数接受两个参数：

• iY：当前位置在网格中的纵坐标。
• parameters：一个引用，指向一个 IncomprFlowParam<T> 类型的对象，该对象包含了泊肃叶流的相关参数，如参考速度、通道高度等等。

第20行：将纵坐标 iY 转换为物理空间中的坐标 y，这里用到了 parameters 对象中的 getResolution() 方法，它返回了一个标准单位在网格尺度中的长度，从而将纵坐标转换为了物理空间中的坐标。

注意：该函数计算了给定高度的通道中的抛物线泊肃叶流速度剖面，并给出通道中心处的最大速度。对象**IncomprFlowParam**将在下面进行讨论：它存储了模拟的各种参数。

第21行：4. * parameters.getLatticeU() * (y - y * y) 计算了当前位置的速度值。其中，parameters.getLatticeU() 返回了参考速度，在泊肃叶流中这是流体在通道中心的最大速度；(y - y * y) 对应了流体速度沿通道宽度方向的变化，这是泊肃叶流的一种典型速度分布，形状是一个抛物线。因此这个函数的计算方式也被称为“抛物线型泊肃叶流速度分布”。

最后，该函数返回了当前位置的速度值。

这是所谓的函数对象或函数式的例子。它是一个类，重载了函数调用运算符operator()。该类的实例就像普通函数一样，但该对象比一个纯函数更智能。例如，它可以在构造函数中接受参数，并拥有可配置的行为。本函数对象封装了函数**poiseuilleVelocity，并配置了一个IncomprFlowParam**类型的对象。这是一个有用的技巧，可以避免将模拟参数声明为全局变量并使它们对所有人可访问，就像在代码1.2中所做的那样。

这是一个定义了一个Poiseuille边界条件的类PoiseuilleVelocity。下面是每一行代码的解释和作用：

• 第25行：template <typename T> 使用泛型编程的方式定义了一个类型T，这样PoiseuilleVelocity类可以用于任何类型T。
• 第26行：class PoiseuilleVelocity { 定义了一个名为PoiseuilleVelocity的类。
• 第27行：public: 这个类的公共元素，这意味着这些元素可以从类的外部访问。
• 第28-30行： 定义了一个构造函数，用于初始化PoiseuilleVelocity类的参数。
• 第33行： 定义了一个函数操作符，将Poiseuille边界条件的速度分量赋值给u。
• 第34行：将Poiseuille边界条件的速度分量赋值给u的第一个元素。
• 第35行：将Poiseuille边界条件的速度分量赋值给u的第二个元素为0。
• 第39行：定义了一个函数操作符，将Poiseuille边界条件的速度和密度分别赋值给u和rho。
• 第40行：将Poiseuille边界条件的速度分量赋值给u的第一个元素。
• 第41行：将Poiseuille边界条件的速度分量赋值给u的第二个元素为0。
• 第42行：将Poiseuille边界条件的密度赋值为1。
• 第44行：private:: 这个类的私有元素，这意味着这些元素只能在类的内部访问。
• 第45行：这个类的一个私有元素，用于存储传入的参数。

这段代码定义了一个名为channelSetup的函数，该函数包含三个参数：

• 一个名为lattice的MultiBlockLattice2D对象
• 一个名为parameters的IncomprFlowParam对象
• 一个名为boundaryCondition的OnLatticeBoundaryCondition2D对象。

该函数的主要目的是为lattice对象设置流体通道的初始状态，具体步骤如下：

• 第54行：使用boundaryCondition对象的setVelocityConditionOnBlockBoundaries函数，创建速度边界条件。
• 第55行：指定域的所有外部边界实现速度的Dirichlet 边界条件。此时，边界上的速度值仍未定义。
• 第57-59行：使用setBoundaryVelocity函数，指定速度边界条件。该函数使用PoiseuilleVelocity<T>(parameters)对象作为边界速度函数。
• 注意：第58行开始定义了外部边界上速度的Dirichlet边界条件。此时，边界上的速度值尚未定义。然后，将速度边界节点的值从分析的Poiseuille剖面中定义。请注意，尽管此函数应用于整个域，但它仅对先前定义为Dirichlet边界节点的节点产生影响。
• 第62-63行：使用initializeAtEquilibrium函数，创建初始条件。该函数使用PoiseuilleVelocity<T>(parameters)对象作为速度函数。
• 第65行：调用lattice.initialize()函数，初始化多块晶格对象。

这段代码定义了一个名为writeGifs的函数，它的作用是将给定的多块格子数据结构lattice中的速度场数据以GIF格式写入到磁盘上，同时指定图片的尺寸为imSize。函数的第二个参数iter表示当前的迭代步数，用于命名输出文件的后缀。函数的具体实现包括以下几行：

• 第70行：定义一个常量imSize，表示输出GIF图片的尺寸为600x600像素。
• 第71行：创建一个ImageWriter对象，用于将数据写入到磁盘上的文件。
• 第72-74行：调用imageWriter对象的writeScaledGif函数，将速度场数据以GIF格式写入到磁盘上。其中第一个参数指定了输出文件名的前缀，第二个参数是一个指向计算速度范数的函数的指针，第三个和第四个参数分别指定了输出GIF图片的宽度和高度。

2.5.3.5 主函数 main

这段C++代码是整个程序的主函数，其作用是运行整个流体模拟的主要逻辑。

• 第77行：定义了一个返回整型值的主函数，它带有两个参数argc和argv，分别是命令行参数的数量和指向参数字符串的指针数组。
• 第78行：Palabos库初始化，处理输入的命令行参数。调用了Palabos库中的plbInit()函数，用于初始化Palabos库和MPI。
• 第80行：设置输出文件夹路径，该路径是相对于当前目录的路径，所有的输出文件都将保存在该目录下，以便进行后续的分析和可视化。

这段C++代码定义了一个类型为IncomprFlowParam<T>的变量parameters，用于在不可压缩流体情况下从无量纲变量转换为格子单位。变量parameters通过一个类构造函数进行初始化，该构造函数的参数列表包括：

• (T)1e-2：无量纲最大速度，它是Poiseuille流中的最大速度，并以格子单位表示。
• (T)100.：雷诺数。
• 100：参考长度的分辨率（通道高度）。
• 2.：无量纲变量中的通道长度。
• 1.：无量纲变量中的通道高度。

注意：类型为 IncomprFlowParam 的对象存储模拟参数（雷诺数、格子分辨率、晶格速度中的参考速度等），并执行单位转换。例如，它用于从雷诺数自动计算弛豫参数omega。

这段代码定义了两个常量imSave和maxT，以及调用了一个函数writeLogFile。

• 第93行：定义了一个常量imSave，表示保存GIF图片的时间间隔，以无量纲时间为单位，即每隔0.1个无量纲时间就保存一次图片。
• 第95行：定义了一个常量maxT，表示总的模拟时间，以无量纲时间为单位，即模拟3.1个无量纲时间的流动。
• 第98行：调用了一个函数writeLogFile，将流体的参数和流动名称写入日志文件，方便后期查看和分析。其中，参数parameters表示流体参数，"Poiseuille flow"是流动名称。

这段代码定义了一个名为lattice的对象，这个对象是一个二维的多块格子（MultiBlockLattice2D），并使用模板参数T和DESCRIPTOR来指定其数据类型和描述符类型。这个对象被初始化为一个 BGK 动力学（BGKdynamics），这是一种处理流体动力学的算法，它使用了一个叫做 BGK 模型的经典的格子 Boltzmann 方法来模拟液体的流动。在这里，BGKdynamics 用于描述流体的宏观运动和它的本质相互作用之间的联系。这个对象被初始化为一个二维网格，其大小由parameters.getNx()和parameters.getNy()指定。它的初始状态将根据指定的初始条件来设置，这个条件将在下面的代码中被定义。

这段代码的作用是创建边界条件，并将边界条件应用到模拟的通道中。

• 第104-105行：创建一个指向OnLatticeBoundaryCondition2D类对象的指针，并使用createLocalBoundaryCondition2D()函数将其初始化为一个新的局部边界条件对象。这个指针将用于将边界条件应用于模拟通道中。
• 注意：选择用于实现边界条件的算法。类型 LocalBoundaryCondition 实现了规则化的边界条件，完全是局部的（不需要访问相邻节点）。类型 InterpBoundaryCondition**** 使用有限差分方案计算壁面应变率张量，因此是非局部的。
• 第107行：调用名为channelSetup()的函数，该函数的参数为模拟通道、IncomprFlowParam对象和OnLatticeBoundaryCondition2D对象指针。在channelSetup()函数中，边界条件被设置，用于给通道设置边界速度，并初始化通道的初始条件。最后调用lattice.initialize()对通道进行初始化。

这段代码实现了主要的计算过程。在这里，我们用一个for循环执行时间步长的迭代。每次迭代，程序都会检查是否需要将当前状态写入输出文件中，然后通过执行Lattice Boltzmann迭代来更新流场的状态。当完成所有迭代时，我们会释放创建的边界条件内存。

• 在第110行，使用for循环来迭代时间步长。这个循环将一直执行，直到总的时间达到预设的最大值maxT。
• 在第111-113行，通过检查当前时间是否是保存GIF图像的时间点，来决定是否将当前状态写入GIF文件。如果是，则调用writeGifs函数生成文件。
• 在第116行，执行Lattice Boltzmann迭代来更新流场状态。这里调用lattice的collideAndStream()函数实现。
• 在119行，释放之前创建的边界条件内存，防止内存泄漏。

2.5.4 结果输出

2.5.5 附录

2.5.5.1 tutorial_1_5.cpp 源码

2.5.5.2 开头注释