定制后处理应用程序

当OpenFOAM中不存在执行所需计算的类似应用程序或函数对象时，可编程定制后处理应用程序。

本节介绍了一个计算上升气泡面积和速度的示例后处理应用程序的编程。在这个例子中，上升的气泡是用一个用于跟踪气泡的场的水平集来近似的。用interIsoFoam求解器对上升气泡进行了模拟[4,5]。interIsoFoam求解器实现了非结构分段线性界面计算(PLIC)流体体积方法[3]的一种变体（参考[2]最近的评论）。

在interIsoFoam求解器中，由符号距离重建等值面（ $ψ_{R D F}$ ）[5]和体积分数 $ψ_{α}$ 之间存在一个选择。由于重构气泡中的差异可以忽略不计，因此在本例中，我们使用 $ψ_{α} \equiv α$ 作为水平集场。因此，气泡是一个水平集： $Σ (t) = {x : α_{c} (x) = 0.5}, x \in R^{3} (8.1)$ ${α_{c}}_{c} \in C$ 体积分数与网格中心相关联，在 $x$ 点处近似为 $α_{c} := α_{c} (x)$ 以重建等值面。 OpenFOAM中的等值面重建算法实现了正则化推进四面体算法[6]，该算法在与 $\sum (t)$ 相交的网格边对 $α_{c} (x)$ 使用线性逼近。利用反向距离加权(IDW)插值方法，从以网格为中心的体积分数 ${α_{c}}_{c} \in C$ 求取网格角点体积分数值。

等值面算法的任务是将气泡表面重建为三角形网格，从场场等值面体积分数场（本例）。在这个例子中，气泡的面积是三角形等值面网格的面积，并且气泡质心同样是网格等值面的中心根据气泡质心的时间演化计算气泡速度。

本章依赖于来自OFPrimer存储库的Cases/Chapter04/RisingBubble2D案例。

示例应用程序名为isoSurfaceBubbleCalculator，它也可以在源代码库中找到，位于文件夹：

applications/utilities/postProcessing/

isoSurfaceBubbleCalculator使用的库的名称为bubbleCalc，其源代码位于文件夹：

src/bubbleCalc

该示例从实现bubbleCalc库开始，然后是isoSurfaceBubbleCalculator后处理应用程序。

// 程序清单20
isoSurfacePoint
(
    const volScalarField& cellValues,
    const scalarField& pointValues,
    const scalar iso,
    const isoSurfaceParams& params = isoSurfaceParams(),
    const bitSet& ignoreCells = bitSet()
);

使用等值面对上升的气泡建模是向现有类添加功能的原型示例。对于等值面，此类称为isoSurfacePoint，其源文件位于文件夹：

$FOAM_SRC/sampling/surface/isoSurface

当然，需要扩展的类通常是未知的，但可以通过在扩展代码指南中搜索具有类似功能的类来找到它，

isoSurfacePoint类的接口显示，类构造函数负责重构等值面三角形网格，如程序清单20所示。

isoSurfacePoint构造函数对上升气泡规定了以下要求：

应在单元格角点（pointValues）处计算附加物理场
成员函数尚不可用于面积、质心或上升速度计算。

上述项目为isoSurfacePoint类提供了附加功能。在后处理应用程序中使用全局函数和数据实现新计算是可能的，但这使得无法在其他设置中重复使用它们：函数对象或另一个应用程序。

因此，等值面计算将被封装到一个名为isoBubble的新类中，并编译到一个可动态加载的共享库BubbleCalc中。将后处理代码编译到共享库中，可以在不同的应用程序和其他可能需要的库之间共享功能。一个例子可以是计算两个气泡之间的符号距离。

isoSurfacePoint需要一个与单元角点关联的附加pointField来计算等值面，而气泡重建仅依赖于以单元为中心的volScalarField，因此pointField的计算被封装到一个小类中，如清单21所示。此示例中使用的是反向距离加权插值，也可以使用其他插值方法。

//程序清单21
class isoPointFieldCalc
{
    // Computed point iso-field.
    scalarField field_;
    public:
        isoPointFieldCalc(const volScalarField& isoField)
        {
            this-> calcPointField(isoField);
        }
        virtual void calcPointField(const volScalarField& isoField)
        {
            volPointInterpolation pointInterpolation (isoField.mesh());
            field_ = pointInterpolation.interpolate(isoField)();
        }
        virtual ~isoPointFieldCalc() = default;
        const scalarField& field() const
        {
            return field_;
        }
};

清单22显示了isoBubble的私有类接口。

class isoBubble:public regIOobject
{
    isoPointFieldCalc isoPointField_;
    // Bubble geometry described with an iso-surface mesh.
    isoSurfacePoint bubbleSurf_;
    // Number of zeros that are prepended to the timeIndex()
    // for written VTK files.
    label timeIndexPad_;
    // Output format.
    word outputFormat_;
    fileName paddWithZeros(const fileName& input) const;

从类声明开始，注意isoBubble继承自regIOobject类。OpenFOAM将对象的IO操作封装在regIOobject类中，该类注册到对象注册表。对象注册表是一个类，当在应用程序代码中执行写操作时，它调用所讨论对象的写成员函数。它也是程序中的全局对象，向订阅它的所有对象发出写入请求。这种设计代表了一种称为“观察者模式”（[1]）的常见面向对象设计（OOD）模式，并且由于以下原因（其中），它非常适合CFD求解器：

一个面向对象的CFD应用程序可能会使用许多对象
对象不必在每个时间步长被写出。写入频率遵循从字典文件读取并由注册表应用的用户定义的规则。

调用runtime.write()使Foam::Time注册表循环访问一系列指向已注册对象（在观察者模式中为“主题”）的指针，并在每个时间步骤结束时将write调用转发给每个已注册对象。只有当有问题的时间步是为输出指定的时间步时，才会发生实际的写入。例如，用户可以规定输出每N个时间步进行一次。

为了使用OpenFOAM中可用的写入控制（例如，每5个时间步长或每0.01秒）写入isoBubble对象，isoBubble类从regIOobject类公开继承。isoBubble的padWithZeros私有成员函数扩展文件名，以便paraView可视化应用程序将其作为时间序列进行读取。然后，指定VTK格式的实际输出将从isoBubble委托给isoSurface类。这些文件被写入实例目录中，其名称是作为IOobject构造函数参数的子参数给出的。

isoBubble类的私有属性为：

isoPointFieldCalc，用于根据以网格为中心的等值面场值计算网格角点等值面场值
isoSurfacePoint执行等值面重构并存储曲面网格数据
等值面格式化输出所需的数据和功能

在OpenFOAM模拟中，点场是默认不可用的，OpenFOAM在数学模型中使用网格中心物理场作为因变量。 isoSurfacePoint不存在的空构造函数和以网格为中心(vol)场作为参数的构造函数排除了从isoSurfacePoint继承的可能性。取而代之的是使用组合，isoPointCalc对象计算isoSurfacePoint构造函数所需的点场。

通过将计算出的isoPointField和所有必要的参数传递给构造函数，isoSurfacePoint类负责等值面网格的实际重构，如清单23所示。

// 程序清单23
isoBubble::isoBubble
(
    const IOobject& io,
    const volScalarField& isoField,
    scalar isoValue,
    bool isTime ,
    label timeIndexPad,
    word outputFormat,
    bool regularize
)
:
    regIOobject(io, isTime),
    isoPointField_(isoField),
    bubbleSurf_(isoField, isoPointField_.field(), isoValue),
    timeIndexPad_(timeIndexPad),
    outputFormat_(outputFormat)
{}

isoSurfacePoint的成员函数可用于执行气泡的中心和面积计算。我们将气泡的面积计算为等值面面积法向矢量的大小之和： $A_{b} = f \sum ∥ S_{f} ∥ (8.2)$ 气泡中心计算为所有等值面网格点的代数平均值： $C_{b} = \frac{1}{N} N \sum x_{N} (8.2)$ 其中N是等值面网格的点数。计算由isoBubble::area和isoBubble::center成员函数执行，如清单24所示。

// 程序清单24
scalar isoBubble::area() const
{
    scalar area = 0;
    const pointField& points = bubblePtr_-> points();
    const List<labelledTri> & faces = bubblePtr_-> localFaces();
    forAll (faces, I)
    {
        area += mag(faces[I].normal(points));
    }
    return area;
}

vector isoBubble::center() const
{
    const pointField& points = bubblePtr_-> points();
    vector bubbleCenter (0,0,0);
    forAll(points, I)
    {
        bubbleCenter += points[I];
    }
    return bubbleCenter / bubblePtr_-> nPoints();
}

在重用了isoSurface类并封装了IO操作以及气泡区域和中心的计算之后，现在可以轻松地编写应用程序并实例化任意多的气泡对象。 isoBubble类是为可供选择的计算而准备的，并且可以与其他与气泡相关的计算一起扩展。

为了方便地与其他人（人或客户端程序）共享isoBubble的实现，isoBubble实现被设置为某个库（即bubbleCalc库）的一部分。列出 $PR I ME R_{E} X A MP L E S_{S} RC / b u bb l e C a l c 目录会显示库文件： ‘‘‘ ba s h ? > l s i so B u bb l e . C i so B u bb l e . d e p i so B u bb l e . H l n I n c l u d e M ak e ‘‘‘ 库构建规范在 mak e 文件夹中的 o pt i o n s 和 f i l es 文件给出。 ‘ f i l es ‘ 文件列出了编译到库目标代码中的文件，链接器随后使用该代码绑定到在包含的头文件中声明并在应用程序代码中调用的函数调用。该文件还指定了编译后的库的安装位置，在本例中，按照用户定义库的建议，库安装在文件夹$ foam_user_libbin中，以避免污染OpenFOAM系统的库安装目录。

isoBubble.C
LIB = $(FOAM_USER_LIBBIN)/libbubbleCalc

options文件列出了所有必要的目录路径，这些目录包含isoBubble中使用的头文件，以及安装库代码的目录路径。预编译库代码允许我们扩展现有代码，而不必总是重新编译所有代码。通过在bubbleCalc目录中发出命令wmake libso，使用OpenFOAM构建系统wmake编译该库。

EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/triSurface/lnInclude \
    -I$(LIB_SRC)/surfMesh/lnInclude \
    -I$(LIB_SRC)/surfMesh/MeshedSurfaceProxy/ \
    -I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
    -lmeshTools \
    -lfiniteVolume \
    -ltriSurface \
    -lsampling

通过编译的库，可以在示例后处理实用程序isoSurfaceBubbleCalculator中使用isoBubble类。

示例后处理实用程序isoSurfaceBubbleCalculator可在代码存储库的applications/utilities/ postProcessing文件夹中找到。

首先，必须包含isoBubble类的声明，才能使库在应用程序中使用。此外，若要处理既有模拟结果，请使用时间步长提取器。相应的代码片段如清单27所示。

// 程序清单27
#include "fvCFD.H"
#include "timeSelector.H"
#include "isoBubble.H"
using namespace bookExamples;

接下来，为应用程序分配一个命令行选项，然后初始化根用例、模拟时间和网格，如清单28所示。

// 程序清单28
int main(int argc, char *argv[])
{
    argList::addOption
    (
        "isoField",
        "Name of the field whose level set is the bubble surface."
    );
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    if (!args.found("isoField"))
    {
        FatalErrorInFunction
            << "Provide the name of the iso-surface field."
            << "Use '-help' for more information."
            << abort(FatalError);
    }

负责实际计算的应用程序部分如清单29所示。

Foam::instantList timeDirs = Foam::timeSelector::select0(runTime, args);
OFstream bubbleFile("bubble.csv");
bubbleFile << "TIME,AREA,CENTER_X,CENTER_Y,CENTER_Z\n";
forAll(timeDirs, timeI)
{
    runTime.setTime(timeDirs[timeI], timeI);
    Info<< "Time = " << runTime.timeName() << endl;

    volScalarField isoField
    (
        IOobject
        (
            args.get<word> ("isoField"),
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        ),
        mesh
    );

    isoBubble bubble
    (
        IOobject
        (
            "bubble",
            "bubble",
            runTime,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        isoField
    );

    const auto area = bubble.area();
    const auto center = bubble.center();

    Info << "Bubble area = " << bubble.area() << endl;
    Info << "Bubble center = " << bubble.center() << endl;
    bubbleFile << runTime.timeOutputValue() << ","
        << area << ","
        << center[0] << "," << center[1] << "," << center[2] << "\n";
    bubble.write();
}

从regIOobject继承会强制客户端代码使用IOobject初始化isoBubble对象。IOobject构造函数的参数定义如下：

“bubble”。对象的名称（在本例中，是用于保存等值面数据的文件名）
“bubble”。实例，或存储对象的目录（在本例中与文件名相同）
runTime 。对象注册到的对象注册表-我们希望isoBubble的IO操作由模拟时间来调节，
IOobject::。为IOobject定义运行时wirte/read操作的标记（枚举）
inputField。跟踪的iso-field。

然后使用timeSelector类在模拟案例目录的文件和文件夹中查找时间步长（迭代）目录。在每次迭代中，读取等值线场，重建等值线气泡对象，计算气泡的面积和中心，并将气泡等值线曲面网格写入磁盘。从这个意义上说，isoBubble类的行为方式与任何其他OpenFOAM类的行为方式相同。

通过调用Allrun脚本来运行模拟，初始化气泡，运行interIsoFoam解算器并调用以计算气泡等参曲面。

>  isoSurfaceBubbleCalculator -isoField alpha.air

气泡曲面作为一系列“vtk”文件存储在气泡目录中。物理时间、气泡面积和气泡中心坐标存储在bubble.csv中。

图8.3中等值面的可视化显示了OpenFOAM中的2D等值面算法与ParaView中的2D等值面算法之间的一些差异（在isoSurface.H文件中提到）。

图8.3 α=0.5时OpenFOAM和ParaView等值面，最终时间步长

在Jupyter笔记本iso-bubble.ipynb中，使用Pandas Python数据分析库处理buble.csv格式的isoSurfaceBubbleCalculator的CSV输出。清单30列出了Python pandas代码。这个使用pandas的小例子只是为了向读者介绍pandas库：pandas可以进行比这个例子中演示的更强大的数据处理。不过，有两个优点可以提高对数据处理代码的理解：阅读CSV文件很简单，列通过它们的名称访问，导数的计算非常简单。除了这些简单的优点外，panda在CFD环境中最大的好处是它能够轻松地存储和处理使用panda.MultiIndex with pandas. DataFrame进行参数研究的结果。

import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams["figure.dpi"] = 200
rcParams["text.usetex"] = True
bubble_csv = pd.read_csv("bubble.csv")
plt.plot(bubble_csv["TIME"], bubble_csv["AREA"])
plt.xlabel("Time in $s$")
plt.ylabel("Area in $m^2$")
plt.savefig("iso-bubble-area.png")
plt.plot(bubble_csv["TIME"], bubble_csv["CENTER_Y"])
plt.ylabel("Y-coordinate of the bubble center")
plt.xlabel("Time in $s$")
plt.savefig("iso-bubble-center-y.png")
y_velocity = bubble_csv["CENTER_Y"].diff() / bubble_csv["TIME"].diff()
plt.plot(bubble_csv["TIME"], y_velocity)
plt.ylabel("Bubble y-velocity in $m/s$")
plt.xlabel("Time in $s$")
plt.savefig("iso-bubble-center-y-velocity.png")

气泡面积随时间的分布如图8.4所示。当气泡上升到溶液域的顶壁并在其上扩展时，面积几乎减半。图8.5给出了用式（8.3）计算的气泡质心的Y坐标分布。该分布在肉眼看来足够平滑。然而，用有限差分法从式（8.3）中气泡质心的Y坐标计算气泡上升速度的Y−分量容易出错。图8.6证实了这种精确度的灾难性损失。以下建议的练习改进了示例后处理计算中气泡上升速度的评估。

图8.4 气泡面积随时间的分布

图8.5 气泡质心的Y坐标

图8.6 气泡Y向速度

练习：

如果用更精确的时间有限差分格式计算上升速度，会发生什么？扩展isoBubble以计算气泡速度，作为在气泡等值面的三角形中心处计算的平均速度？提示：使用interpolationCellPoint在三角形中心插入速度。

编写另一个后处理应用程序，仅使用alpha.air和U场计算气泡速度。计算速度的差异是什么？这些差异来自哪里？