本文许可遵循 Creative Commons Attribution 4.0 International License.

Part I: 低雷诺数条件下的圆柱绕流

有关本问题的数学描述，请参阅 Schäfer and Turek (1996) 的文章。 目前的 OpenFOAM 模拟是 Darshan Thummar 和 Fabian Gabriel 创建的设置修改版本。

检查并运行测试用例模拟

每个 OpenFOAM 模拟都包含文件夹和文件的层次结构。 最上层的模拟文件夹始终包含以下三个子文件夹：

1. system - 用于求解器执行、方程离散化、线性求解器、网格划分等的控制字典
2. constant - 传输和热物理属性、网格（polyMesh 文件夹）和几何文件
3. a time folder; 物理场写入时间文件夹，其中文件夹名称反映拍摄快照的物理时间； 通常，有一个 0（零）文件夹来定义每个场的初始条件； 0.org 文件夹表示在开始模拟之前可以修改一个或多个物理场的初值，例如设置更复杂、不均匀的初始条件

本练习的模拟文件夹位于 test_cases/Column2D。 在检查和运行模拟之前，我们在 exercises 文件夹中创建一个副本。 如果您尚未创建练习文件夹，请运行：

# assuming you are at the repository's top level
mkdir exercises

对于本讲座中的所有 OpenFOAM 模拟，我们使用专用的 Apptainer 镜像。 已修改的常用 OpenFOAM 命令的逻辑隐藏在仓库顶层的 RunFunctions 脚本中。 该脚本包含“runApplication”和“runParallel”等 shell 函数，它们简化了 OpenFOAM 应用程序的执行并将常规输出以及错误消息写入日志文件。 这些函数又依赖于环境变量，这些变量定义了 Apptainer 镜像名称以及你的系统上讲座仓库的位置。 脚本 setup-env 定义所有必需的 shell 变量：

# run at the repository's top level
source setup-env

您可以通过键入 *$M{L}_{C}FD\ast 并按TAB键两次（快速）来检查新定义的变量。\ast \ast 请注意，每当打开新终端时，您都必须获取环境变量。\ast \ast 或者，您可以将上面的行添加到"$HOME/.bashrc"文件中。

现在让我们递归复制模拟文件夹（-r 选项），然后导航到新文件夹：

cp -r test_cases/cylinder2D/ exercises/
cd exercises/cylinder2D/

您可以使用 ls 命令来检查构成模拟案例的文件和文件夹。 测试用例的顶层有两个脚本。 Allrun 脚本执行所有必要的预处理步骤并使用 Apptainer 容器启动模拟。 Allclean 脚本将文件夹重置为其初始状态。

在开始模拟之前，请尝试浏览 Column2D 文件夹内的所有文件和文件夹来回答以下问题：

• 使用多少个 MPI 等级（多少个处理器核心）执行模拟？ 提示：大多数 CFD 代码中的并行化基于域分解。
• 对于哪些场我们需要初始条件和边界条件？ 提示：检查初始时间文件夹。
• 写下定义入口速度的函数。 提示：我们为哪个场定义此边界条件？
• 哪个实用程序用于域离散化？ 提示：查看官方用户指南中有关网格生成的部分。
• 基于平均入口速度和圆柱直径的雷诺数是多少？ 提示：尝试分四个步骤回答这个问题：
  • 测定运动粘度
  • 确定圆柱体的直径
  • 确定平均入口速度
  • 将这些数字代入 $Re$ 的定义中
• 数值求解器是可压缩的还是不可压缩的？ 提示：在扩展代码指南中搜索求解器名称。
• 使用什么线性求解器来求解压力泊松方程？ 提示：检查 fvSolution 字典。
• 动量方程的对流项是如何离散的？ 提示：检查 fvSchemes 字典。
• 采用哪种湍流模型？
• 您希望创建多少个时间文件夹？ 提示：检查controlDict。
• 哪个文件包含计算力系数的字典？

现在，让我们实际运行模拟：

./Allrun

模拟大约需要 20-30 分钟才能完成。 同时，尝试回答以下问题：

• 创建了哪些新文件夹？
• 处理器[0-1]文件夹包含什么？
• 网格由多少个单元组成？ 让以下步骤指导你：
  • 获取使用容器的函数：source $ML_CFD_BASE/RunFunctions
  • 运行 checkMesh 实用程序：runApplication checkMesh
  • 检查 log.checkMesh 文件
• 力系数写入文件coefficient.dat； 该文件位于哪里？ Tip: use find . -type f -name coefficient.dat
• coefficient.dat 的前四列包含什么（标题之后）？ Tip: use head -n 20 path/to/coefficient.dat and refer to the documentation for the acronyms

在我们开始检查和可视化流场之前，浏览求解器的日志文件始终是一个好主意。 打开 log.pimpleFoam 文件，最好使用 vim、nano 或 emacs 等命令行编辑器，并尝试回答以下问题：

• 完成模拟需要多少秒？ 提示：转到日志文件的最末尾
• 每个时间步长 p-U 耦合平均需要多少次迭代？
• 最大库朗数（大约）是多少？

检查日志文件应该始终是第一个操作，以防您的模拟未按预期执行（崩溃、意外结果、执行速度极慢）。 有一个非常有用的实用程序可以从名为 foamLog 的日志文件中提取有关残差和迭代解决方案的信息。 它可以按如下方式使用：

source $ML_CFD_BASE/RunFunctions
runApplication foamLog log.pimpleFoam

将创建一个新的 logs 文件夹。 为了快速查看一些残差和迭代计数，我们使用一个名为 Gnuplot 的小工具。 如果您的系统上未安装 Gnuplot，请运行：

sudo apt install gnuplot

然后，导航到 logs 文件夹，并绘制一些可用文件，例如：

cd logs
gnuplot
# now we are in a gnuplot shell
# continuity error after the first p-U-coupling iteration over time
gnuplot> plot 'contCumulative_0'
# sometimes, it's helpful to use a logarithmic scale on the ordinate
gnuplot> set logscale y
# final residual for p and Ux after fist p-U-coupling iteration
gnuplot> plot 'pFinalRes_0', 'UxFinalRes_0'
# iterations to solve for p in the first p-U-coupling iteration
gnuplot> unset logscale y
gnuplot> plot 'pIters_0'
# exit Gnuplot -> q + Enter

在 ParaView 中可视化流程

我们来看看流场。 要打开 ParaView，请运行：

# run this command at the simulation folder's top level
paraview post.foam

post.foam 是一个空的虚拟文件，用于指示 ParaView 我们正在提供 OpenFOAM 模拟。 **单击“应用”之前，**将“案例类型”设置为“分解案例”。 要检查网格，请选择 Surface With Edges 作为表示。 然后，返回到 Surface 表示并设置速度大小的动画：

• 在字段下拉菜单中选择 U
• 转到最后一个时间步骤，然后单击按钮 Rescale To Data Range（将鼠标指针悬停在按钮上时会显示文本）
• 返回第一个时间步并单击播放按钮
• enjoy the vortex street~

您还可以通过转到File->Save Animation来保存动画。

该域分为两个子域，每个子域对应一个 MPI 等级（核）。 要可视化分解，请关闭 ParaView，然后执行以下步骤：

touch processor0/post.foam
touch processor1/post.foam
# inspect the first processor domain
paraview processor0/post.foam
# inspect the second processor domain
paraview processor1/post.foam

评估升力和阻力系数

经过圆柱的流动的典型目标量是作用在圆柱上的阻力和升力。 作为本教程的最后一步，我们将在 Jupyter 笔记本中绘制随时间变化的阻力和升力系数。 打开一个新的 Jupyter 笔记本，为文件指定一个合理的名称，并在第一个单元格中添加一些描述，例如：

# Flow past a circular cylinder at low Reynolds number
## Drag and lift coefficients

基于圆柱直径 <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.6944em;"></span><span class="mord mathnormal">d</span></span></span></span>、平均入口速度 <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.8201em;"></span><span class="mord accent"><span class="vlist-t"><span class="vlist-r"><span class="vlist" style="height:0.8201em;"><span style="top:-3em;"><span class="pstrut" style="height:3em;"></span><span class="mord mathnormal" style="margin-right:0.10903em;">U</span></span><span style="top:-3.2523em;"><span class="pstrut" style="height:3em;"></span><span class="accent-body" style="left:-0.2222em;"><span class="mord">ˉ</span></span></span></span></span></span></span></span></span></span> 和运动粘度 <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.4306em;"></span><span class="mord mathnormal" style="margin-right:0.06366em;">ν</span></span></span></span> 的雷诺数为 <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.6833em;"></span><span class="mord mathnormal" style="margin-right:0.00773em;">R</span><span class="mord mathnormal">e</span><span class="mspace" style="margin-right:0.2778em;"></span><span class="mrel">=</span><span class="mspace" style="margin-right:0.2778em;"></span></span><span class="base"><span class="strut" style="height:1.0701em;vertical-align:-0.25em;"></span><span class="mord mathnormal">d</span><span class="mord accent"><span class="vlist-t"><span class="vlist-r"><span class="vlist" style="height:0.8201em;"><span style="top:-3em;"><span class="pstrut" style="height:3em;"></span><span class="mord mathnormal" style="margin-right:0.10903em;">U</span></span><span style="top:-3.2523em;"><span class="pstrut" style="height:3em;"></span><span class="accent-body" style="left:-0.2222em;"><span class="mord">ˉ</span></span></span></span></span></span></span><span class="mord">/</span><span class="mord mathnormal" style="margin-right:0.06366em;">ν</span><span class="mspace" style="margin-right:0.2778em;"></span><span class="mrel">=</span><span class="mspace" style="margin-right:0.2778em;"></span></span><span class="base"><span class="strut" style="height:0.1056em;"></span><span class="mord">...</span></span></span></span>。

Pandas 模块提供了强大的功能来读取逗号分隔值（CSV）文件。 使用 Matplotlib 可以将系数可视化为简单的线图。

import matplotlib.pyplot as plt
from pandas import read_csv

# increase plot resolution
plt.rcParams["figure.dpi"] = 160

在下一个单元格中，使用以下要点作为起点，并通过检查 coefficient.dat 文件并查阅 read_csv 文档。

path = "path/to/coefficient.dat"
names = ["t", "cd", "cl"]
data = read_csv(path, skiprows=???, header=0, sep=???, usecols=[0, 1, 3], names=names)
data.head()

最后，我们创建一个具有两个轴的图形，并将阻力系数和升力系数绘制为线条。 替换绘图函数的指定参数。

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 4), sharex=True)

ax1.plot(data.t, ???)
ax2.plot(data.t, ???)
ax2.set_xlabel(r"<span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.6151em;"></span><span class="mord mathnormal">t</span></span></span></span> in <span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.4306em;"></span><span class="mord mathnormal">s</span></span></span></span>")
ax2.set_ylabel(r"<span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.5806em;vertical-align:-0.15em;"></span><span class="mord"><span class="mord mathnormal">c</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.3361em;"><span style="top:-2.55em;margin-left:0em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mathnormal mtight" style="margin-right:0.01968em;">l</span></span></span></span><span class="vlist-s">​</span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span></span></span></span>")
ax1.set_ylabel(r"<span class="katex"><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.5806em;vertical-align:-0.15em;"></span><span class="mord"><span class="mord mathnormal">c</span><span class="msupsub"><span class="vlist-t vlist-t2"><span class="vlist-r"><span class="vlist" style="height:0.3361em;"><span style="top:-2.55em;margin-left:0em;margin-right:0.05em;"><span class="pstrut" style="height:2.7em;"></span><span class="sizing reset-size6 size3 mtight"><span class="mord mathnormal mtight">d</span></span></span></span><span class="vlist-s">​</span></span><span class="vlist-r"><span class="vlist" style="height:0.15em;"><span></span></span></span></span></span></span></span></span></span>")
ax1.set_ylim(2.8, 3.3)
plt.savefig("drag_and_lift.pdf", bbox_inches="tight")
plt.plot()

Part II: 轴对称气泡在静止液体中的上升

The bubble solver

Basilisk 是一个高度专业化的库，用于模拟空间尺度差异较大的多相流。 与 OpenFOAM 相比，Basilisk 中的模拟不受输入字典控制。 相反，定制的求解器是基于模块化库组件编译的。 test_cases 文件夹中的求解器对液槽中上升的轴对称气泡（围绕上升方向在方位角方向上对称）进行模拟。 求解器本身在 bubble.c 源文件中定义。 其他头文件以 VTK 格式写入输出文件，这对于 ParaView 中的后处理很有用。

Compiling the solver

编译和运行模拟需要 basilisk.sif Apptainer 容器。 我们可以使用类似于 runApplication 和 runParallel 的函数来自动运行容器内的命令。 然而，这一次，我们手动完成该过程，以便更好地理解与容器的交互。 位于源文件夹中的 Makefile 简化了求解器的编译。 从存储库的顶层运行以下命令：

# create a copy for the exercise
cp -r test_cases/basilisk_solver/ exercises/
# open a new shell in the Singularity container
apptainer shell basilisk.sif
# now we are operating inside the container
Apptainer> cd exercises/basilisk_solver/
Apptainer> make &> log.make

通过运行“cat log.make”检查日志文件以查看编译过程中是否遇到任何问题。 在同一文件夹中，您现在应该看到一个名为 bubble 的新二进制文件。

Performing a simulation

求解器采用 4 个参数来控制模拟：

1. 最大精细化级别(the maximum refinement level)； 求解器采用自适应网格细化； 在网格的某些区域使用更密集的网格以减少数值误差； 最大细化级别限制网格细化，使最小的单元不会变得太小（每次细化，时间步长都会下降，单元数会增加）
2. 结束时间(the end time)； 达到该物理时间后模拟停止
3. the Eötvös or Bond number
4. the Galilei number

我们将在第 5 讲中讨论无量纲数的定义和含义。现在，只要知道 $Eo=115$ 的 Eötvös 数与 $Ga=134$ 的伽利略数相结合会导致 所谓的球帽制度。 以下模拟大约在 10-15 分钟内完成。 如果您的机器有超过 2 个可用处理器核心，您可以通过向“mpirun”传递更高的值来减少运行时间，例如“mpirun -np 4”以在 4 个核心上运行模拟。

# create a dedicated folder for this test simulation
# in the exercises/basilisk_solver/ folder
Apptainer> mkdir run
Apptainer> cd run
Apptainer> mpirun -np 2 ../bubble 14 10 115 134 &> log.solver
# once the solver is done, leave the container
exit

您可以通过打开新终端并使用“tail -f log.solver”来跟踪求解器的进度。 要停止跟踪输出，请按 Ctrl+c。 请注意，此组合键仅停止 tail 命令，但不会停止求解器执行。

在 ParaView 中可视化流程

流量变量和气液界面的VTK文件存储在vtu文件夹中。 要打开快照，请导航至 exercises/basilisk_solver/run/vtu/ 并打开 ParaView。

# Note: you have to exit the Apptainer container before opening ParaView
# if you still see the 'Apptainer>' prompt, type 'exit' and press enter
cd exercises/basilisk_solver/run/vtu/
paraview

要打开整组快照，请单击文件->打开->bubble_..pvtu。 有几个有趣的选项可以对多相模拟进行后处理。 以下是一些建议：

• 通过选择 Surface With Edges 表示来可视化网格； 提前模拟时间并观察网格如何变化
• 选择体积分数字段 f； 气泡内的体积分数为 $f=0$，气泡外的 $f=1$，包含气液界面的单元中的 $0<f<1$
• 重建的界面作为线段存储在名为 plic_..pvtu 的文件中； 它们可以像现场快照一样打开； 线段有点难以看清； 以提高他们的可视度，
  • 选择线框表示
  • 将“属性”面板中的“线宽”值设置为 3
  • 在编辑颜色图下选择不同的颜色
• u.x 字段保存绝对速度场； 为了描绘潜在的再循环区域，更适合通过从速度场的 x 分量中减去上升速度来切换到移动参考系； 这个计算可以在 ParaView 中完成：
  • 10 个时间单位后（无量纲）上升速度约为 ${\tilde{U}}_b=0.735$； 请参阅求解器日志文件的第 6 列
  • 在 ParaView 中，应用计算器过滤器； 计算参考下图
  • 应用 Glyph 过滤器来可视化气泡尾流中的再循环区域

恭喜！ 第二次和第三次练习就这样完成了。