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

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


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).

注意接下来的几行是如何调用类中实现的函数.timeName().C().Cf()。同样重要的是,Mesh.C().Cf()返回表示每个单元和内部单元面心的坐标矢量。因此,调用Mesh.C().size()方法可以得到单元的数量。注意,mesh.Cf()仅仅包含内部单元面的面心坐标矢量,不包含边界面!

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;

也可以采用标准的cpp循环形式对每一个网格进行遍历:

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

每一个单元都是由一组单元面构成的,这些面可以是由两个单元共享的内部面或者边界面(在OpenFOAM中边界面也被成为Patch)。每一个内部面都有一个单元被称为owner,而另一侧的单元则被称为neighbour,即:owner,neighbour和内部单元面一一对应,mesh.owner().size()mesh.neighbour().size()mesh.Cf().size()的返回值相同,都是内部单元面的个数。

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)的形式,第一个数为包含面的数量,括号中的数为包含的面的索引值。

除了内部面之外,如果我们还想获取边界网格的信息,可以通过boundaryMesh对象。每个边界面也包含在constant/polyMesh/faces文件中。但是,该文件首先定义内部面及其包含的顶点信息。构成边界面(patch)的单元面的索引及其起止范围则是在constant/polyMesh/boundary文件中定义。OpenFOAM中,所有的边界单元信息都可以从对mesh.boundaryMesh()的遍历中获取。

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();

可以使用Mesh.boundary()[patchI].Cf().size()Mesh.boundary()[patchI].start()方法来检查一个单元面是内部还是位于边界上。

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)的起止面编号。

    forAll(fcs,faceI)
        if (faceI%80==0)
        {
            if (faceI<mesh.Cf().size())
                Info << "Internal face ";
            else
            {
                forAll(mesh.boundary(),patchI)
                    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:";
            forAll(fcs[faceI],vertexI)
                // 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;
    }

边界面集合patch也可以用它们的名字从网格中检索出来。如果用户要从字典中引用一个特定的patch,这可能很有用(就像你在计算一个特定patch的受力情况时那样)。

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. 案例文件下载

点击下载案例文件

Everything not saved will be lost.
最后更新于 2023-05-07