[OpenFOAM] 编程基础 04 Understanding the numerical mesh in OpenFOAM

发布于 2023-05-07  212 次阅读

Please refresh the page if equations are not rendered correctly.

1. 代码分解

求解器的主函数总是以如下形式开始,代码的第7-8行分别会创建runTimemesh 两个实例(instance of object/class)。如果你不熟悉什么是类或对象,强烈建议你访问这个网站,只有当你读完关于类、继承和多态的所有内容后再回来:http://www.cplusplus.com/doc/tutorial/classes/。

#include "fvCFD.H"

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

    // These two create the time system (instance called runTime) 
    #include "createTime.H"
    #include "createMesh.H" // and fvMesh (instance called mesh).


Note how the next lines call functions .timeName(), .C() and .Cf() implemented in the objects. It is also important to realize that mesh.C() and .Cf() return vector fields denoting centers of each cell and internal face. Calling the mesh.C().size() method therefore yields the total size of the mesh.

    Info << "The most recent time folder is " << runTime.timeName() << nl
         << "The mesh has " << mesh.C().size() << " cells and " << mesh.Cf().size()
         << " internal faces in it!" << nl << endl;


It's possible to iterate over every cell in a standard cpp for loop:

for (label cellI = 0; cellI < mesh.C().size(); cellI++)
    if (cellI%20 == 0) // only show every twentieth cell not to spam the screen too much
        Info << "Cell " << cellI << " with centre at " << mesh.C()[cellI] << endl;
Info << endl; // spacer


Each cell is constructed of faces - these may either be internal or constitute a boundary, or a patch in OpenFOAM terms; internal faces have an owner cell and a neighbour.

for (label faceI = 0; faceI < mesh.owner().size(); faceI++)
    if (faceI%40 == 0)
        Info << "Internal face " << faceI << " with centre at " << mesh.Cf()[faceI]
             << " with owner cell " << mesh.owner()[faceI]
             << " and neighbour " << mesh.neighbour()[faceI] << endl;
Info << endl; // spacer


Info << mesh.cells()[CellI] << endl;

其中,CellI是单元的索引值,输出结果为6(49 50 941 1341 11 47)的形式,第一个数为包含面的数量,括号中的数为包含的面的索引值。


Boundary conditions may be accessed through the boundaryMesh object. In reality, each boundary face is also included in the constant/polyMesh/faces description. But, in that file, the internal faces are defined first. In addition, the constant/polyMesh/boundary file defines the starting faceI indices from which boundary face definitions start. OpenFOAM also provides a macro definition for for loops over all entries in a field or a list, which saves up on the amount of typing.

forAll(mesh.boundaryMesh(), patchI)
    Info << "Patch " << patchI << ": " << mesh.boundary()[patchI].name() << " with "
         << mesh.boundary()[patchI].Cf().size() << " faces. Starts at total face "
         << mesh.boundary()[patchI].start() << endl;
Info << endl;


Faces adjacent to boundaries may be accessed as follows. Also, a useful thing to know about a face is its normal vector and face area.

    label patchFaceI(0);  // initialize patchFaceI
    forAll(mesh.boundaryMesh(), patchI)
        Info << "Patch " << patchI << " has its face " << patchFaceI << " adjacent to cell "
             << mesh.boundary()[patchI].patch().faceCells()  //[patchFaceI]
             << ". It has normal vector " << mesh.boundary()[patchI].Sf()[patchFaceI]
             << " and surface area " << mag(mesh.boundary()[patchI].Sf()[patchFaceI])
             << endl;
    Info << endl;

获取内部单元面的法向量可以直接调用mesh.Sf()方法。另一个调取所有面(内部面和边界面)的面积的方法是mesh.magSf()。对于内部单元面,面的法向量从owner指向neighbour,owner的单元格索引比neighbour小。对于边界面,法线总是指向求解域外(它们有不存在的 "虚 "邻居)。

For internal faces, method .Sf() can be called directly on the mesh instance. Moreover, there is a shorthand method .magSf() which returns the surface area as a scalar. For internal faces, the normal vector points from the owner to the neighbour and the owner has a smaller cell index than the neighbour. For boundary faces, the normals always point outside of the domain (they have "imaginary" neighbors which do not exist).


It is possible to look at the points making up each face in more detail. First, we define a few shorthands by getting references to the respective objects in the mesh. These are defined as constants since we do not aim to alter the mesh in any way. NOTE: these lists refer to the physical definition of the mesh and thus include boundary faces.

    const faceList& fcs = mesh.faces();
    const List<point>& pts = mesh.points();
    const List<point>& cents = mesh.faceCentres();


Use can be made of the mesh.boundary()[patchI].Cf().size() and mesh.boundary()[patchI].start() methods to check whether the face is internal or lies on a boundary.

在第8行的forAll循环中,mesh.boundary()[patchI].start()<= faceI) && (faceI < mesh.boundary()[patchI].start()+mesh.boundary()[patchI].Cf().size() 判断条件分别给出的每一个边界面集合(patch)的起止面编号。

        if (faceI%80==0)
            if (faceI<mesh.Cf().size())
                Info << "Internal face ";
                    if ((mesh.boundary()[patchI].start()<= faceI) &&
                        (faceI < mesh.boundary()[patchI].start()+mesh.boundary()[patchI].Cf().size()))
                        Info << "Face on patch " << patchI << ", faceI ";
                        break; // exit the forAll loop prematurely

            Info << faceI << " with centre at " << cents[faceI]
                 << " has " << fcs[faceI].size() << " vertices:";
                // Note how fcs[faceI] holds the indices of points whose coordinates
                // are stored in the pts list.
                Info << " " << pts[fcs[faceI][vertexI]];
            Info << endl;
    Info << endl;

在原来的Cavity教程中,frontAndBack边界被定义为 "empty "类型。这是一个特殊的边界条件案例,可能会导致意外的行为,因为它的.Cf()字段的大小为0。如果遇到类型为empty的patch且存在很大的风险,可以用以下代码预先判断patch的类型,以避免遇到这个问题。

In the original cavity tutorial, on which the test case is based, the frontAndBack boundary is defined as and "empty" type. This is a special BC case which may cause unexpected behavior as its .Cf() field has size of 0. Type of a patch may be checked to avoid running into this problem if there is a substantial risk that an empty patch type will appear.

    label patchID(0);
    const polyPatch& pp = mesh.boundaryMesh()[patchID];
    if (isA<emptyPolyPatch>(pp))
        // patch patchID is of type "empty".
        Info << "You will not see this." << endl;


Patches may also be retrieved from the mesh using their name. This could be useful if the user were to refer to a particular patch from a dictionary (like when you do when calculating forces on a particular patch).

    word patchName("movingWall");
    patchID = mesh.boundaryMesh().findPatchID(patchName);
    Info << "Retrieved patch " << patchName << " at index " << patchID << " using its name only." << nl << endl;

    Info<< "End\n" << endl;

    return 0;

2. 案例文件下载


最后更新于 2023-05-07