OpenFOAM 网格拓扑与计算表面向量和 (2)

OpenFOAM 网格拓扑与计算表面向量和 (2)

在上一篇文章中我们讨论了OpenFOAM的网格拓扑，以及通量的计算方式。在这篇文章中我们来对OpenFOAM的网格与面做一系列的处理。

在 The Finite Volume Method in Computational Fluid Dynamics^[1]书的第六章 Exercise中，有一道练习题 (Exercise 2)，计算每个网格的表面向量和是否为0：

Write a program that reads an OpenFOAM® mesh and checks that for each element the sum of the surface vectors is zero.

本文尝试作如下解答：
首先，新建文件夹Exercise_6_8_2，在其中创建 Exercise6_8_2.C 文件：

#include "fvCFD.H"

int main(int argc, char *argv[])
{
    #include "setRootCase.H"

    // runTime类 与 mesh类的实例化
    #include "createTime.H"
    #include "createMesh.H"
    // 在循环中，按照mesh的数量来遍历所有的网格
    for (label meshID = 0; meshID < mesh.C().size(); meshID++){
        // 初始化 temp 向量
        vector temp = vector::zero;
        // 子循环遍历网格的所有face：
        for (label faceID = 0; faceID < mesh.faces().size(); faceID++){
        	// 如果此face的本网格(owner) 序号等于外层循环的网格序号
            if (mesh.owner()[faceID] == meshID){
                Info << "the face vector is: " << mesh.Sf()[faceID] << endl;
                // 对面法向量求和
                temp += mesh.Sf()[faceID];
            }
         }
		// 再遍历所有的内部面 (InternalFaces)
        for (label neighbourID = 0; neighbourID < mesh.neighbour().size(); neighbourID++){			// 如果内部面中有相邻网格(neighbour) 的序号等于外层循环的网格序号
            if (mesh.neighbour()[neighbourID] == meshID){
                    Info << "the neighbour vector is: " << mesh.Sf()[neighbourID] << endl;
                    // 减去此面法向量
                    temp -= mesh.Sf()[neighbourID];
            }
        }
		// 判定 temp 向量的模值是否为0
        if (mag(temp) != 0){
        	// 如果不为0的话，则输出错误信息 和 网格序号，并跳出循环
            Info << "error and break" << " " << meshID << " " << temp << endl;
            break;
        }
    }
    Info<< "End\n" << endl;
    return 0;

接下来是大家所熟悉的编译环节, 编辑Exercise_6_8_2/Make/files:

Exercise6_8_2.C

EXE = $(FOAM_USER_APPBIN)/Exercise6_8_2

并且 wmake 编译。

接下来，在测试算例 testCase 中执行 blockMesh，画出$4\times 4\times 1$ 的网格:

Writing polyMesh
----------------
Mesh Information
----------------
  boundingBox: (0 0 0) (1 1 0.1)
  nPoints: 50
  nCells: 16
  nFaces: 72
  nInternalFaces: 24
----------------
Patches
----------------
  patch 0 (start: 24 size: 4) name: movingWall
  patch 1 (start: 28 size: 12) name: fixedWalls
  patch 2 (start: 40 size: 32) name: frontAndBack

End

接下来，执行 Exercise_6_8_2，可以看到如下信息：

可以看到，输出了每个 $cell$ 的各个 $face$ 的面法向量，并且把它们做了 $owner$ 与 $neighbour$ 的区分。

程序并没有输出错误信息，可见各个网格的面法向量之和为0。