学习自F. Moukalled, L. Mangani, M. Darwish所著The Finite Volume Method in Computational Fluid Dynamics - An Advanced Introduction with OpenFOAM and Matlab

Chapter 7 The Finite Volume Mesh in OpenFOAM and uFVMhtml

OpenFOAM是强大高效的开源代码，而uFVM则侧重教育学习(便于理解却丧失效率)，本章着重讲解OpenFOAM的网格文件格式，以及uFVM的网格数据结构是如何架构的，可见uFVM与OpenFOAM的实现细节是很是相似的，能够做为学习OpenFOAM的先导。node

1 uFVM

1.1 OpenFOAM测试算例

uFVM代码是直接读入OpenFOAM的算例配置文件的，因此先要理解OpenFOAM的算例是如何设置的，下图展现了一个名为cavity的算例，其设置文件所有放在cavity文件夹下面。c++

在cavity文件夹下，有三个文件夹：web

其中文件夹0下面存放的是初始物理场(含内部量与边界量)的信息，即初始速度场U和初始压力场p存放在0/P和0/U文件中。一旦计算开始后，若设置了输出间隔，则会生成对应时刻标号的文件夹，其里面存放着对应时刻的流场变量。算法

在文件夹system中存放的是与FVM算法相关的三个设置文件：controlDict是计算的控制参数，如模拟的开始和结束时间、使用的时间步长、数据输出间隔等信息；fvSchemes是定义离散格式的，好比梯度格式、插值格式等；fvSolution则定义的是求解算法(求解Ax=b方程的法子，如共轭梯度、多重网格等)、松弛因子、收敛指标等信息。数组

在const文件夹中，transportProperties文件存放的是相关物理特性，如粘性系数；polyMesh中则存放着描述网格信息的文件points、faces、owner、neighbour、boundary共5个文件。有人可能会很好奇，那个blockMeshDict是干啥的？它实际上是用来剖分分块结构网格的配置文件，里面大概是写每一个块的角点以及划分单元数目和单元长度增加比率的信息，写好后，用blockMesh命令就能生成较为简单的分块结构网格了，也就是前面提到的那5个文件信息了。实际上，若是网格不是用blockMesh工具来划分的，而是由别的格式的非结构网格转化而来的，那么在polyMesh文件夹下就见不到这个blockMeshDict文件了。数据结构

下面我们看看polyMesh文件夹中这5个文件是如何给出网格信息的。架构

1.2 polyMesh文件夹

如上图所示(来自OpenFOAM开发者Hrvoje Jasak大神的博士论文)，OpenFOAM中的单元是任意多面体，而每一个面能够由任意多边形构成，因此很是灵活通用，固然结构网格是这种多面体网格的一种特例。app

OpenFOAM中只处理3维网格，若是是2维问题，把它沿着展向拉伸一层网格让它变成3维网格，同时把展向两端的边界面定义成empty就能够了。svg

points

points文件存放的是角点的坐标(

x

i

,

y

i

,

z

i

)

(x_i,y_i,z_i)(xi​,yi​,zi​)列表，注意，这些角点，既是单元的角点，也是面的角点。标识为0的角点(第1个角点)坐标存放在第1行，标识为1的角点(第2个角点)坐标存放在第2行。注意：因为C++中的数组标识是从0开始而非从1开始的，因此第1个量的标识是0，而若是有N个量，则第N个量的标识是N-1(有点反人类，跟天然数究竟是从0开始仍是从1开始有点像哈)。这个points文件的格式以下

#number of points

(

(#x #y #z)

......

)

看看points文件的范例，共有1074个角点：

1074

(

(32 16 0.9377383239)

(33.9429245 16.11834526 0.9377383239)

(35.84160614 16.46798134 0.9377383239)

(37.67648315 17.04080009 0.9377383239)

(39.42799377 17.82870483 0.9377383239)

(41.07658768 18.82359314 0.9377383239)

(...)

...

)

faces

也是个列表，存放的是构成面的角点标识列表，一样，第1行是0号标识面的全部角点标识，第2行是1号表示面的全部角点标识。faces文件的格式以下

#number of faces

(

#number of points for face 1 (#p1 #p2 #p3 ...... )

#number of points for face 2 (#p1 #p2 #p3 ...... )

......

)

faces文件的范例以下，有3290个面(全部面=内部面+边界面)，这里每一个面由4个角点组成

3290

(

4(36 573 589 52)

4(41 578 634 97)

4(44 81 618 581)

4(30 82 619 567)

4(121 50 587 658)

4(39 120 657 576)

......

)

owners

owners文件存储着面的所属(owner)单元标识，即，标识为0的第1个面的owner单元标识存放在第一行，标识为1的第2个面的owner单元标识存放在第2行，以此类推。注意：owner单元的数目等于面的总数目，即内部面数目+边界面数目。

单元的数目等于owner的最大值+1(由于标识从0开始的)。

owners文件的格式以下

#number of owners

(

#owner of face1

#owner of face2

......

)

owners文件范例以下

3290

(

0

1

2

3

4

5

6

......

)

单元总数也能够从owners文件的头部信息行中获取，nCells后面跟的就是单元数目

note "nPoints:1074 nCells:918 nFaces:3290 nInternalFaces:1300";

neighbours

neighbours文件存放的是面的邻居单元(neighbour)标识列表，注意邻居单元的数目等于内部面的数目，由于边界面只有一个单元，就是own单元，边界面是没有neighbour单元的。

neighbours文件的格式为

#number of neighbour

(

#neighbour of face1

#neighbour of face2

......

)

neighbours文件的范例为

1300

(

22

68

29

96

31

34

......

)

boundary

boundary文件存放的是计算域的边界列表，每一个边界类型所含的面将做为一个patch存在，并赋予名字。每一个边界patch的类型将用其所含面的总数(nFaces)和起始面标识(startFace)来指定。换言之，对于面而言，首先标识的是内部面，而后再标识外部边界面，而每一个boundary patch上的面将被连续标识，以便在boundary中指定名字与类型。

boundary patch的格式为

#boundary patch name

{

type #patchtype;

nFaces #number of face in patch set;

startFace #starting face index for patch;

}

boundary patch的范例以下，即从1300号到1399号面命名为wall-4，其边界类型为壁面。

wall-4

{

type wall;

nFaces 100;

startFace 1300;

}

有人可能会很好奇？为啥没有单元的面标识列表呢？由于彻底不须要，owners和neighbours就彻底标清楚了面和单元的编码关系，由这些关系就足以获得element-faces的标识关系，不必再画蛇添足地整上一个单元所含面标识的列表文件出来。并且，FVM的计算中大可能是作面循环，而较少会用到单元循环，因此以面为基础的数据架构是很实用的。

1.3 uFVM网格

uFVM读入的是OpenFOAM的网格文件，用的是cfdOpenFoamMesh(2018年V1.5版本中是cfdReadPolyMesh)函数，这个函数依次去读取points、faces、owner、neighbour、boundary文件信息，并后处理来得到额外的elements、points的拓扑信息。

网格读取并重构拓扑关系成功后，将存储在一个结构体当中，以elbow算例为例，完成cfdOpenFoamMesh后，mesh的信息以下(注意，在V1.5版本的uFVM中是把这些几何信息和拓扑信息所有一股脑地存放在global变量Region的mesh中，而在老版本中(书上那样子)是分开来分别存放在mesh下面的nodes、faces、elements和boundaries里面，这里仍是以书本为例把，毕竟分开存放更便于理解数据架构，虽然寻访的时候要两层剥离稍微麻烦点)。mesh的详情以下

m = cfdReadOpenFoamMesh('elbow')

m =

nodes: [1x1074 struct]% 1074个节点的信息列表

numberOfNodes: 1074 % 节点总数为1074

caseDirectory: 'elbow'% 算例名字elblow

numberOfFaces: 3290% 面总数为3290

numberOfElements: 918% 单元总数为918

faces: [1x3290 struct]% 3290个面的信息列表

numberOfInteriorFaces: 1300% 内部面总数为1300

boundaries: [1x6 struct]% 6个boundary patch的信息列表

numberOfBoundaries: 6% 边界的数目为6(这个量好像跟下面的重复了！)

numberOfPatches: 6% 边界patch的数目为6

elements: [1x918 struct]% 918个单元的信息列表

numberOfBElements: 1990% 边界单元数目为1990

numberOfBFaces: 1990% 边界面的数目为1990，这个与边界单元的数目是一致的

这个网格能够用cfdPlotMesh函数将其可视化

在nodes列表中，每一个元素为一个结构体，表明着第i个节点的信息(注意不一样于c++，matlab列表的下标是从1开始算的)，那么第1个节点的信息以下

n1= m.nodes(1)

n1 =

centroid: [3x1 double]% 节点形心坐标，即节点的坐标x,y,z

index: 1% 节点总体标识(编号，编码)

iFaces: [172 328 1355 1386 1677 1891 1893]% 节点被哪些面所享有，这些面的标识列表

iElements: [112 219 220]% 节点被哪些单元所享有，这些单元的标识列表

在faces列表中，每一个元素为一个结构体，表明着第i个面的信息，以第3个面为例

m.faces(3)

ans =

iNodes: [45 82 619 582]% 构成该面的节点列表

index: 3% 该面的标识

iOwner: 3% 该面owner单元标识

iNeighbour: 30% 该面neighbour单元标识(若为边界面，则该标识为-1)

centroid: [3x1 double] % 该面形心坐标x,y,z

Sf: [3x1 double]% 该面的面积矢量Sx,Sy,Sz

area: 5.3046% 该面的面积S

CN: [3x1 double]% 该面所属单元形心到该面形心的距离矢量Cf

geoDiff: 4.5940% ？面几何扩散系数 gDiff_f = Ef / CF，见第8章

T: [3x1 double]% 面所属单元和邻近单元形心之间的距离矢量CF？(感受更像是Tf=Sf-CF为非正交修正矢量，见第8章内容)

gf: 0.4226% 面插值中的几何权重系数gf

walldist: 0% 面所属单元形心到壁面的垂直距离(某些湍流模型中会用到)

iOwnerNeighbourCoef: 1% ？

iNeighbourOwnerCoef: 1% ？

在elements列表中，每一个元素为一个结构体，存放着第i个单元的信息，以第20个单元为例

m.elements(20)

ans =

index: 20% 该单元标识

iNeighbours: [100 103]% 该单元的邻近单元(与该单元共享面的那些单元)标识列表

iFaces: [33 34 1317 1493 1494] % 构成该单元的面标识列表

iNodes: [168 79 616 705 617 80] % 构成该单元的节点标识列表

volume: 3.2484% 该单元体积

faceSign: [1 1 1 1 1]% 该单元的构成面是否为其owner面(==1)或neighbour(==-1)

numberOfNeighbours: 2% 该单元邻近单元数目

centroid: [3x1 double]% 该单元的形心坐标x,y,z

注意，单元标识和面标识是统一的，以保证二者间在逻辑上的属从关系，边界面是在内部面以后才开始编号的。能够用cfdPlotElements来标明特定的单元，如

cfdPlotElements([20 300])

m.elements(300)

ans =

index: 300

iNeighbours: [278 302 590]

iFaces: [407 435 436 2053 2054]

iNodes: [283 820 679 142 290 827]

volume: 1.9083

faceSign: [-1 1 1 1 1]

numberOfNeighbours: 3

centroid: [3x1 double]

20单元位于左上角蓝色的，300单元位于右下角红色的表示。20单元确实有两个邻近面，有5个构成面(2维网格展向拉伸成3维后，展向先后还有俩面)；300单元有3个邻近面，有5个构成面，其体积是小于20单元的；与前面给出的信息是一致的。

最后给出的是boundary patches的信息，存放在boundaries列表中，每一个元素为一个结构体，存放着第i个边界片的信息，以第1个边界片为例

>> m. boundaries(1)

ans =

userName: 'wall-4'% 该边界片的名字(用户随便起的名字)

index: 1% 该边界片的标识

type: 'wall'% 该边界片的类型(物理类型，这里是壁面)

numberOfBFaces: 100% 该边界片的面总数为100个

startFace: 1301% 该边界片的起始面标识为1301

也就是说，第1-1300都是内部面，而从1301开始的后面那些面都是边界面(再啰嗦一句，C++下标从0开始，matlab下标从1开始，因此这里uFVM中是1301，而OpenFOAM中是1300)，那么我们看下1301号面这个边界面的信息是什么

>> m.faces(1301)

ans =

iNodes: [38 53 590 575]

index: 1301

iOwner: 1

iNeighbour: -1

centroid: [3x1 double]

Sf: [3x1 double]

area: 3.7510

CN: [3x1 double]

geoDiff: 5.6264

T: [3x1 double]

gf: 1

walldist: 0.6667

iOwnerNeighbourCoef: []

iNeighbourOwnerCoef: []

这些信息中能够发现边界面与内部面的不一样之处，即边界面的neighbour单元是没有的(标识为-1)，边界面的几何权重系数gf为1。

那么，若是要对对某个边界patch作循环的话，该如何处理呢？对这个边界patch的起始面和终止面循环就行了，好比要对boundary patch的第2个patch的面作循环，能够这样子

theMesh = cfdGetMesh;

iPatch = 2;

iBFaces = cfdGetFaceIndicesForBoundaryIndex(iPatch)

for iBFace = iBFaces

theBFace = theMesh.faces(iBFace);

disp(theBFace) %display theBFace internal fields

end

cfdGetFaceIndicesForBoundaryIndex是这样定义的

theIndices = cfdGetFaceIndicesForBoundaryIndex(theBoundaryIndex)

%

theBoundary = cfdGetBoundary(theBoundaryIndex);

theNumberOfBFaces = theBoundary.numberOfBFaces; % 该片边界所含面的数目

theStartFace = theBoundary.startFace;% 起始面标识

theIndices = [theStartFace:theStartFace+theNumberOfBFaces-1];% 起始面标识 到 起始面标识+边界片的面总数-1

%

end

最后，再看下新版本(2018 V1.5版本)uFVM中的mesh信息，其实和上面是同样的，只是糅一块儿了。

global Region

>> Region.mesh

ans =

nodeCentroids: [1074x3 double]

numberOfNodes: 1074

faceNodes: {3290x1 cell}

numberOfFaces: 3290

owners: [3290x1 double]

numberOfInteriorFaces: 1300

numberOfBFaces: 1990

neighbours: [1300x1 double]

numberOfElements: 918

numberOfBElements: 1990

cfdBoundaryPatchesArray: {6x1 cell}

numberOfBoundaryPatches: 6

closed: 0

elementNeighbours: {918x1 cell}

elementFaces: {918x1 cell}

elementNodes: {918x1 cell}

upperAnbCoeffIndex: [1300x1 double]

lowerAnbCoeffIndex: [1300x1 double]

nodeElements: {1074x1 cell}

nodeFaces: {1074x1 cell}

elementCentroids: [918x3 double]

elementVolumes: [918x1 double]

faceCentroids: [3290x3 double]

faceSf: [3290x3 double]

faceAreas: [3290x1 double]

faceWeights: [3290x1 double]

faceCF: [3290x3 double]

faceCf: [3290x3 double]

faceFf: [3290x3 double]

wallDist: [3290x1 double]

wallDistLimited: [3290x1 double]

1.4 uFVM Geometric Field(不一样几何位置上的变量场)

在FVM的求解的过程当中，已知变量、待求变量、中间变量的信息常常被存放在不一样的地方，如单元形心、面形心、角点(节点)处，而变量的类型又有标量、向量、矢量之分，所以便有了位于单元、面、节点处的标量、向量、矢量场这么众多不一样类型的场。

1.4.1 Element Fields(存储在单元上的场)

单元场能够用以下函数来定义

cfdSetupMeshField(theUserName, theLocale, theType, theTimeStep)

其中theUserName为该场的名字，theLocale为该场在几何意义上的存储位置(Elements， Faces， Nodes)，即存放在单元、面仍是节点上，theType为存放变量的类型(Scalar或Vector)，标量仍是矢量，最后的theTimeStep代表该场的不一样时刻(Step0， Step1， 等)，例如

>> UField = cfdSetupMeshField('U:water','Elements','Vector','Step0')

UField =

userName: 'U:water'

name: 'U_fluid01'

type: 'Vector'

locale: 'Elements'

phi: [2908x3 double]

便定义了一个存储在单元形心上的矢量场，其存储的时间步是Step0，即，当前时间步。

如上图，该变量数组的长度是NumberOfElements + NumberOfBoundaryFaces，即，单元数目 + 边界面数目。也就是说，虽然是存储在单元上的场，可是把边界面也当成一种特殊类型的单元，因此单元场上存的变量值，既有每一个单元上的值，**也有边界面上的值，边界面的值是做为单元场的边界条件存在的！**它们是按照先单元后边界面的顺序排列的。

那么，若是要把UField在patch1上的边界值设成[1 0 0]该如何作呢？以下

% get the mesh

theMesh = cfdGetMesh;% 提取网格

% get information about the boundary patch% 提取boundary patch信息

theBoundary = theMesh.boundaries(iPatch);% 拿出第iPatch个boundary patch的信息

numberOfElements = theMesh.numberOfElements;% 总体单元数目

numberOfInteriorFaces = theMesh.numberOfInteriorFaces;% 总体内部面数目

numberOfBFaces = theBoundary.numberOfBFaces;% 总体边界面数目

% Starting face

iFaceStart = theBoundary.startFace;% 第iPatch个boundary patch的起始面标识

% get information about starting and ending elements

% 第iPatch个boundary patch在Element Fields中，起始单元标识的计算为

% 起始单元标识 = 总体单元数目 + (第iPatch个边界patch的开始面标识 - 内部面数目)

iElementStart = numberOfElements+iFaceStart-numberOfInteriorFaces;

% 第iPatch个boundary patch在Element Fields中，结束单元标识的计算为

% 结束单元标识 = 起始单元标识 + (第iPatch个边界patch的面数目 - 1)

iElementEnd = iElementStart+numberOfBFaces-1;

% define the indices as an index array

iBElements = iElementStart:iElementEnd;% 起始单元标识 到 结束单元标识

>> UField.phi(iBElements,:) =

cfdComputeFormulaAtLocale('[1;0;0]','BPatch1','Vector')

ans =

1 0 0

1 0 0

1 0 0

...

...

其中的

cfdComputeFormulaAtLocale(theFormula,theLocale,theType)

计算在theLocale处的theFormla的值，而且返回特定长度的theType(Scalar 或 Vector)数组。

1.4.2 Face Fields(存储在面上的场)

面场是存储位置在面上的场，即theLocale设置为’Faces’的场

cfdSetupMeshField(theUserName, theLocale, theType, theTimeStep)

这里，数组的长度就是总体面的数目了，即内部面数目 + 外部边界面数目，即numberOfFaces = numberOfInteriorFaces + numberOfBoundaryFaces

对于某个边界片boundary patch的边界面的访问方法以下

theMesh = cfdGetMesh; % 提取网格

numberOfElements = theMesh.numberOfElements; % 总体单元数目

numberOfInteriorFaces = theMesh.numberOfInteriorFaces; % 内部面数目

theBoundary = theMesh.boundaries(iPatch); % 提取出第i个边界片boundary patch

numberOfBFaces = theBoundary.numberOfBFaces; % 第i个boundary patch(该片边界)的边界面数目

%

iFaceStart = theBoundary.startFace;% 该片边界的开始边界面标识

iFaceEnd = iFaceStart+numberOfBFaces-1; % 该片边界的结束边界面标识 = 起始面标识 + 边界面总数 - 1

iBFaces = iFaceStart:iFaceEnd;% 该片边界面标识范围

%

iElementStart = numberOfElements+iFaceStart-numberOfInteriorFaces; % 该片边界的起始单元标识

iElementEnd = iElementStart+numberOfBFaces-1;% 该片边界的结束单元标识

iBElements = iElementStart:iElementEnd;% 该片边界的单元范围

thBFaces = theMesh.faces(iBFaces)% 取出该片边界的边界面信息(几何与拓扑信息)

若是要取出该片边界boundary patch的边界单元值，则可用(phi是单元场)

phi_b = phi[iBElements];

1.4.3 Node Fields (存储在节点上的场)

这个比较简单，就是在节点上存放的变量场，数组的长度就是节点的总体数目，节点也不用区分什么内部边界外部边界的，因此没啥好赘述的。

1.5 uFVM网格的对应操做

在离散和求解过程当中，最多见的操做莫过于对单元、内部面、边界面、边界单元、或边界片(boundary patch)的循环遍历操做了。那么它们要如何实现呢？

1.5.1 单元循环遍历

作单元标识从1到numberOfElements(0到numberOfElements-1)的循环就行了，不论是对单元几何拓扑信息的提取，仍是对单元场变量的提取，都是同样的，以下

for iElement=1:numberOfElements

theElement = theMesh.elements(iElement) % 取出第iElement个单元的拓扑几何信息

phi(iElement) % this is field phi at centroid of element iElement

% 取出第iElement个单元形心上存放的物理量phi值

....

end

若是要对边界单元作循环处理，以下

% 对于边界单元(只有单元场才有边界单元的概念)作循环

% 循环范围从 总体单元数目+1 到 总体单元数目 + 边界面数目(边界面数目等于边界单元的数目)

for iBElement = numberOfElements + 1: numberOfElements + numberOfBFaces

phi(iBElement) = 0;% 将边界单元上的phi值设为0

end

1.5.2 面循环遍历

面的编码是先作内部面编码，而后再依次对边界面作编码，所以，若对内部面作循环，只须要从1到numberOfInteriorFaces循环就行了，即

for iFace=1:numberOfInteriorFaces

theFace = theMesh.faces(iFace)% 取出第iFace个面的几何拓扑信息

end

若是要对边界面作循环，则

for iBFace= numberOfInteriorFaces+1:numberOfFaces

theBFace = theMesh.faces(iBFace)

end

若是要对特定某片的boundary patch的边界面作循环，则

startFace = theMesh.boundaries(n).startFace% 第n个boundary patch的起始面标识

nFaces = theMesh.boundaries(n).nFaces% 第n个boundary patch的边界面数目

for iBFace = startFace: startFace+nFaces-1% 循环范围从起始面 到 起始面+该片边界面总数目-1 便可

...

end

固然，能够用函数cfdGetFaceIndicesForBoundaryIndex直接得到循环范围startFace: startFace+nFaces-1。

1.6 计算Gauss Gradient

在uFVM中，对于单元形心处梯度值的计算是在函数CFDComputerGradientGauss0中进行的，0表示没有非正交修正操做，具体细节能够参考函数代码，因为这部份内容和第9章的梯度计算是重复的，而第9章中除了说起单元中心梯度，还讲了面中心梯度，节点梯度的算法，因此直接参考第9章的代码讲解就行了。

2 OpenFOAM

to be continued…