openFOAM post代码结构（另附mesh类）

PS: 这是一个对师兄写的后处理代码的阅读笔记。没有普遍性。

找到一个类似的代码：读取OpenFOAM数据文件并进行后处理-物联网技术文章-傲云电气网 (weixiuzhan.cn)

运行步骤就C++一样：

1. 要检查参数、公式和输出路径（主要读post.cpp，foampost.cpp文件）

2. 编译（./make.sh）

3. 运行（./run.sh）。

Ps=010这个代码里面有三个后处理的文件夹：

post

foampost-pg-derived

foampost-vof

接下来对每一个进行我自己的解读：

post

这个没有编译，因为variables.h里面明显是二维的参数，post.cpp这些是三维的。估计要修改。

• post.cpp：主函数，调用了需要的函数进行各种步骤
• read_para.hpp：从文件读各种parameter
• stat_info.hpp：这个也是二维的，各种accumulate和collect的命令，输出文件有result_ver_es.dat，result_hor_es.dat，result_hor_usx.dat，result_hor_usy.dat，result_hor_us.dat，result_gs_ins.dat，result_gs_avg.dat
• manage_mem.hpp：给各个变量分配内存
• post.h：头文件
• variables.h：post的namespace定义，2维的

 编译不动，不知道咋啦。

foampost-pg-derived

• foampost.cpp：后处理数据，输出axial_es.dat，gs_avg.dat文件
• readFields.H：读取p和Us

通常在实验中利用，取两点测两点之间压降，在通过压降公式求两点之间空隙率。这和我们说的轴向分布不一样，轴向分布是同一高度下先时间平均再空间平均，时空平均值随轴向高度的变化（也就是后面的foampost-vof文件）。

outletUsAvg:

这个已经编译好了，直接运行就行

./run.sh
或者
./run.sh >&log&

可能输出三个文件(runTime.value()>10)：

log：日记

gs_avg.dat：outlet的固相平均Gs（因为openFoam从0开始，N-1也就是最后一个patch了其实）

axial_es.dat：epsilon_g，这里的epsilon是根据压力求的，这种方法认为重力=压差力，通过压差计算两点之间的平均空隙率。

 关于polyMesh和fvMesh 

 凡是int函数内的调用头文件，大概率是插入一段代码。这也是OF的套路之一。

 在文件中使用

#include "createMesh.H"

等效于
    fvMesh mesh
    (
        IOobject
        (
            fvMesh::defaultRegion,
            runTime.timeName(),
            runTime,
            IOobject::MUST_READ
        )
    ); 
这段代码

1. polyMesh和fvMesh都能调用网格信息，但是网格信息fvMesh给的更多更细。
2. 一般计算过程都是调用fvMesh。
3. 通常不会显式使用fvMesh，而是使用#include "createMesh.H"来调用。OF中有很多这种玩法要注意。

mesh类

OpenFOAM_C++学习 - 知乎 (zhihu.com)

mesh.C();这个返回的是向量场，代表每个单元的中心。简称:体心 center
mesh.Cf();这个返回的是向量场，代表内部面的中心。简称:面心 face center

mesh.C().size();这个应该是整型，代表所有的单元个数。
mesh.C()[cellI];这种用法要理解。

mesh.owner().size();内部面的个数。

mag()是大小，Sf()是surfaceVectorField，mag是surfaceScalarField

内部面

1. 在类实例中可以直接使用.Sf()方法计算。这个是面积法向量
2. 可以直接使用.magSf()计算面积，这是个标量
3. 对于内部面法向量，明确规定由ower指向neighbour
4. owner的单元编号小于neighbour的单元编号
5. 对于边界面，法向永远指向域外，这个域外是不存在的，是想象出来的。

OpenFOAM中计算某一个边界的面积，代码为：

计算OpenFOAM中某一个边界的面积_黄岛主的博客-CSDN博客

// 找到inlet这个patch的ID
label patchID = mesh.boundaryMesh().findPatchID("inlet"); 
 
// 创建inlet的polyPatch类型的引用
const polyPatch& cPatch = mesh.boundaryMesh()[patchID]; 
 
// 创建面标量场magSf，代表的是网格单元某一个面的面积
const surfaceScalarField& magSf = mesh.magSf(); 
 
// 定义标量，并初始化为0.0
scalar patchArea = 0.0; 
 
// 遍历inlet上的所有单元面，然后将单元面的面积加起来
forAll(cPatch,faceI) 
{ 
    patchArea += magSf.boundaryField()[patchID][faceI]; 
} 
 
// 输出inlet这个边界面的面积
Info << "Patch Area " << patchArea << endl; 

 或者

label patchID = mesh.boundaryMesh().findPatchID("inlet");
 
 
const polyPatch& cPatch = mesh.boundaryMesh()[patchID];
 
 
const surfaceScalarField& magSf = mesh.magSf();
 
 
// 直接使用gSum函数来求inlet的面积
scalar AREA = gSum(magSf.boundaryField()[patchID]);
 
 
Info << "AREA = " << nl << AREA << endl;

foampost-vof

• foampost.cpp：后处理数据，输出axial_es.dat，gs.dat，gs_avg.dat文件
• readFields.H：读取voidfraction和Us

这个已经编译好了，直接运行就行

./run.sh
或者
./run.sh >&log&

可能输出4个文件(runTime.value()>10)：通常用越长时间平均越好，但是也得看计算量

log：日记

gs.dat：不同运行时间（>10）空间（>LOUT=2.6）的Gs，在>2.6范围的cell的气体速度的体积求和然后除以总体积求平均。算是outlet了吗？

axial_es.dat：纵向设置了50个level，算的是每个level的总和，求平均epsilon，然后取每个level的中心为坐标。轴向分布是同一高度下先时间平均再空间平均，时空平均值随轴向高度的变化。

gs_avg.dat：时间（>10）空间（>LOUT=2.6）平均Gs，也就是gs.dat的求和平均

孔隙率求解方法

 可以通过测量压降来计算，这是实验上的方法。