OpenFOAM直接调用Python:pybFoam案例

免费CFD计算器（欢迎使用！！！）

本文直接基于cavity算例把读取、修改、计算和写回走一遍。算例只有400个网格，结果本身没什么好研究的，主要看这个接口到底能不能用。

安装

我这里的环境为OpenFOAM v2606、Python 3.10.12和pybFoam 0.5.3。

先加载OpenFOAM环境，再安装pybFoam：

source /home/jml/OpenFOAM/OpenFOAM-v2606/etc/bashrcpython3 -m pip install --user pybFoam

pybFoam会针对当前OpenFOAM版本编译扩展，并不是装完一个纯Python包就结束了。我这里第一次编译花了几分钟。

运行Python脚本前也要加载OpenFOAM环境，否则会找不到libfiniteVolume.so等动态库。我第一次直接运行时就报了下面这个错误：

ImportError: libfiniteVolume.so: cannot open shared object file

算例设置

使用OpenFOAM v2606官方cavity算例：

tutorials/incompressible/icoFoam/cavity/cavity

计算域为0.1 m × 0.1 m，网格为20 × 20，共400个单元。上壁面以1 m/s向右运动，其余壁面固定，运动黏度为0.01 m²/s，使用icoFoam计算到0.5 s。

这里没有先运行blockMesh，而是在Python中读取blockMeshDict并生成网格：

from pybFoam import Time, argList, dictionary, fvMeshfrom pybFoam.meshing import generate_blockmeshrun_time = Time(argList([str(case), "-case", str(case)]))generate_blockmesh(　run_time,　dictionary.read(str(case / "system" / "blockMeshDict")),)mesh = fvMesh(run_time)print(mesh.nCells())

输出的网格数为400，和blockMeshDict中的设置一致。

NumPy修改OpenFOAM场变量

网格生成以后，直接读取0/U：

from pybFoam import volVectorFieldimport numpy as npU = volVectorField.read_field(mesh, "U")U_view = np.asarray(U["internalField"])

此时U仍然OpenFOAM中的volVectorField，U_view是它的NumPy视图，数组形状为：

(400, 3)

400对应网格数量，3对应三个速度分量。这里不是把U复制一份再交给NumPy，修改U_view时，OpenFOAM管理的那块内存也会随之改变。

为了看得明显一些，我按网格中心坐标写入一个二维旋涡初场：

Ux = sin(πx/L) cos(πy/L)Uy = -cos(πx/L) sin(πy/L)

其中L = 0.1 m，代码如下：

centres = np.asarray(mesh.C()["internalField"])x = centres[:, 0]y = centres[:, 1]L = 0.1U_view[:, 0] = np.sin(np.pi*x/L)*np.cos(np.pi*y/L)U_view[:, 1] = -np.cos(np.pi*x/L)*np.sin(np.pi*y/L)U_view[:, 2] = 0.0U.correctBoundaryConditions()write(U)

write(U)执行以后，0/U中的internalField由原来的uniform (0 0 0)变成了包含400个矢量的nonuniform List<vector>。

随后直接运行icoFoam：

icoFoam | tee log.icoFoam

求解器可以正常读取Python写回的U场，并计算到0.5 s。

在Python中调用fvc

pybFoam不只是把OpenFOAM文件读成NumPy。求解完成以后，还可以在Python中直接调用OpenFOAM的有限体积算子。

下面先切换到最后一个时间步，再读取U：

instants = list(selectTimes(run_time, ["postProcess"]))positive = [t for t in instants if float(str(t)) > 0.0]latest = positive[-1]run_time.setTime(latest, len(positive))U = volVectorField.read_field(mesh, "U")

速度梯度可以直接用fvc.grad计算：

grad_U = fvc.grad(U)()write(grad_U)

计算得到OpenFOAM中的volTensorField，而不是Python自己用相邻点差分出来的数组。write(grad_U)执行后，最终时间目录中会生成：

0.5/grad(U)

如果还要在Python中继续处理，同样可以把这个张量场映射为NumPy数组：

grad_values = np.asarray(grad_U["internalField"])omega_z = grad_values[:, 1] - grad_values[:, 3]

OpenFOAM张量分量的排列为xx、xy、xz、yx、yy、yz、zx、zy、zz，上面的计算对应二维流动的面外涡量：

ωz = ∂Uy/∂x - ∂Ux/∂y

最后再从U中提取接近x = 0.05 m的竖直中心线。20 × 20网格的单元中心并不正好落在x = 0.05 m上，所以程序自动选择距离最近的一列，即x = 0.0475 m，共20个采样点。

x_line = np.unique(x)[np.argmin(np.abs(np.unique(x) - 0.05))]mask = np.isclose(x, x_line)order = np.argsort(y[mask])y_line = y[mask][order]Ux_line = values[mask, 0][order]np.savetxt("vertical_midline_Ux.csv",　np.column_stack([y_line, Ux_line]),　delimiter=",",)

0.5 s时，计算得到的最大速度模为0.8527 m/s，平均速度模为0.1866 m/s。这个网格太粗，本文也没有与Ghia等经典数据对比，所以这些数值只用于确认读取和后处理流程，不作为cavity算例的精度验证。

pybFoam和PyFoam不一样

这两个名字很像，但不是一个东西。

PyFoam主要在OpenFOAM外部工作，用来修改字典、运行求解器和分析日志。pybFoam则把OpenFOAM类直接绑定到Python，Python拿到的可以是fvMesh、volVectorField以及fvc、fvm算子。

这次实际跑下来，下面这条数据通路已经可以工作：

它比较适合做批量场处理、自定义初始化、在线后处理和数据集生成。以后如果要把机器学习模型接进计算过程，也不用只靠反复读写csv传数据。

目前这个接口还比较新，安装时要和具体OpenFOAM版本一起编译，能绑定的类和算子也没有覆盖OpenFOAM全部功能。不过相比只在外面调用命令，它已经可以直接碰到OpenFOAM的网格、场和离散算子了，这一点还是挺有意思的。