FVM in CFD 学习笔记_第6章_有限体积网格

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

划分网格是有限体积法的首要工作，即，把原本连续的计算域划分成一系列不重叠的单元集合，这些单元由一系列面构成（对于每个面来说，其分属于两个单元），同时还要标记边界面来定义物理边界。随后，还要计算网格相关的几何信息（比如单元体积、面积矢量等），以及拓扑关系（单元-面，面-角点的编号关系）。这些是如何处理的呢？本章便着重阐明有限体积法网格中的几何信息和拓扑信息。

1 区域离散

如图所示，对于（a）所示的计算域，可以用结构网格剖分（b），也可用非结构网格剖分（c），将原本控制方程在连续区域上的解转化为控制方程在离散区域上的解，最终获得的是离散点上的变量值（即，有限个离散点处的值）。

结构网格的拓扑关系是规整的，所以无需额外存储，其编程较为简单、计算效率较高，然而其适用的几何计算域较为苛刻，并不是任何计算域都能很容易地划分成结构网格的，当然，可以用分块结构网格的剖分方式来提高其几何适应性，但是工作量很大，且难度较大。

非结构网格的拓扑关系比较复杂，需要额外存储（单元-面，面-角点的关系表），且编程较为复杂（需要考虑网格的非正交修正），计算效率也比较低，然而非结构网格的优点也很明显，其对于复杂几何区域的适应性很强，其划分网格的过程非常方便快捷。

实际上，非结构网格是更一般的情形，而结构网格则是非结构网格的一种特例而已，因此非结构网格的所有计算是兼容结构网格的。

在OpenFOAM中，使用的是非结构网格，但是也能使用分块结构网格和分块嵌套结构网格。下面就先讲讲简单的结构网格，然后再讲讲复杂得非结构网格。

2 有限体积网格

以梯度计算为例，看看结构网格和非结构网格都需要存储什么几何量和拓扑关系。

梯度计算所需的网格信息

如图，变量一般存储在单元形心$C$处，若采用Green-Gauss方法来计算有限体积单元形心处的变量梯度值$\nabla {\varphi }_C$，有
$$\nabla {\varphi }_C=\frac{1}{V_C}\sum _{f=1}^{nb(C)}{\varphi }_f{\vec{S}}_f$$
即，用面形心处的变量值${\varphi }_f$与面积矢量${\vec{S}}_f$复合后，加和，再除以单元体积$V_C$，便可得到单元形心处的变量梯度值了。

而面形心处的变量值${\varphi }_f$则一般用该面两侧单元形心处变量值的加权平均来获取，即
$${\varphi }_f=g_F{\varphi }_F+g_C{\varphi }_C$$
一种计算权重系数的方法是
$$g_F=\frac{V_C}{V_C+V_F}{g}_C=\frac{V_F}{V_C+V_F}=1-g_F$$
这里只是做简要叙述，梯度的具体计算方法（含非正交修正方法），权系数的其它插值方法将在后续章节中详解。

可见，要存储的几何信息和拓扑关系还真不少呢，有

单元形心坐标值、面形心坐标值、面积矢量值、单元标识与组成单元的面标识的对应关系、面标识与面左右两侧单元标识的对应关系，等

3 结构网格

结构网格的几何信息和拓扑信息要简单明了得狠，以下图为例

在结构网格中，$i,j,k$方向分别划分$N_i,N_j,N_k$个单元，则单元总数为$N_i\ast N_j\ast N_k$，且直接可用标识$(i,j,k)$索引到特定的单元，而特定单元的邻近单元可由标识$(i,j,k)$加或减1来寻访到，至于邻近单元的数目，3维情况单元为六面体单元，每个单元有6个面、8个角点、有6个邻近单元，2维情况单元为四边形单元，每个单元有4个面、4个角点、有4个邻近单元。由于拓扑关系明晰，所以无需额外存储拓扑信息。

3.1 拓扑信息

在计算的时候，最终还是要生成整体矩阵来做计算的，那么局部标识$(i,j)$与整体标识$n$是如何关联的呢？其实很简单，如果整体标识是按照一行一行来标记变量的话，有
$$n=i+(j-1)N_{i}1\le i\le N_{i}1\le j\le N_j$$
那么，单元$(i,j)$邻近单元的局部标识和整体标识为
$$\begin{gathered} (i,j)\to n \\ (i+1,j)\to n+1 \\ (i-1,j)\to n-1 \\ (i,j+1)\to n+N_j \\ (i,j-1)\to n-N_j \end{gathered}$$
如果是3维问题，则
$$n=i+(j-1)N_i+(k-1)N_i{N}_{j}1\le i\le N_{i}1\le j\le N_{j}1\le k\le N_k$$
单元$(i,j,k)$邻近单元的局部标识和整体标识为
$$\begin{gathered} (i,j,k)\to n \\ (i+1,j,k)\to n+1 \\ (i-1,j,k)\to n-1 \\ (i,j+1,k)\to n+N_j \\ (i,j-1,k)\to n-N_j \\ (i,j,k+1)\to n+N_i{N}_j \\ (i,j,k-1)\to n-N_i{N}_j \end{gathered}$$
非常方便，组装整体矩阵的时候也很快捷，因为在整体矩阵中的存储位置一目了然，不烦取寻找对应关系了。

3.2 几何信息

访问某特定单元的局部几何信息较为简便，对单元$(i,j)$来说，其构成面为$\overrightarrow{S1}(i,j)$，$\overrightarrow{S2}(i,j)$，$\overrightarrow{S1}(i+1,j)$，$\overrightarrow{S2}(i,j+1)$，由于构成单元的面矢量是向外为正的，所以有
$${\vec{S}}_{i-1/2,j}=-\overrightarrow{S1}(i,j){\vec{S}}_{i,j-1/2}=-\overrightarrow{S2}(i,j)$$
另外俩面的正方向与坐标轴正方向相同，故
$${\vec{S}}_{i+1/2,j}=\overrightarrow{S1}(i+1,j){\vec{S}}_{i,j+1/2}=\overrightarrow{S2}(i,j+1)$$

3.3 访问单元场

以计算梯度为例，面$(i+1/2,j)$处的值可由两侧单元值插值获取，即
$${\varphi }_{i+1/2,j}=g_{i+1/2,j}{\varphi }_{i+1,j}+(1-g_{i+1/2,j}){\varphi }_{i,j}$$
单元形心梯度值可以算得，即
$$\begin{gathered} \nabla {\varphi }_{i,j}=\frac{1}{V_{i,j}}\left({\varphi }_{i+1/2,j}{\vec{S}}_{i+1/2,j}+{\varphi }_{i-1/2,j}{\vec{S}}_{i-1/2,j}+{\varphi }_{i,j+1/2}{\vec{S}}_{i,j+1/2}+{\varphi }_{i,j-1/2}{\vec{S}}_{i,j-1/2}\right) \\ =\frac{1}{V_{i,j}}\left({\varphi }_{i+1/2,j}{\overrightarrow{S1}}_{i+1,j}-{\varphi }_{i-1/2,j}{\overrightarrow{S1}}_{i,j}+{\varphi }_{i,j+1/2}{\overrightarrow{S2}}_{i,j+1}-{\varphi }_{i,j-1/2}{\overrightarrow{S2}}_{i,j}\right) \end{gathered}$$
如果用上图（b）的整体标识，则
$$\nabla {\varphi }_n=\frac{1}{V_{i,j}}\left({\varphi }_{n+1/2}{\vec{S}}_{n+1/2}+{\varphi }_{n-1/2}{\vec{S}}_{n-1/2}+{\varphi }_{n+N_i/2}{\vec{S}}_{n+N_i/2}+{\varphi }_{n-N_i/2}{\vec{S}}_{n-N_i/2}\right)$$
不难发现，除了边界面之外，每个内部面都被两个单元享有，那么这个面的面积矢量对于其中一个是正的，而对于另一个是负的，那么计算的时候统一算好了每个内部面的面积矢量，然后根据标识关系去添加正负号就行了（对于单元(i,j)，面标识比(i,j)大的为正，面标识比(i,j)小的为负），而不用一把每个面的面积矢量都重复算两次了。

还有一种写成上图（c）所示的东西南北的局部标识形式
$$\nabla {\varphi }_C=\frac{1}{V_C}\left({\varphi }_e{\vec{S}}_e+{\varphi }_w{\vec{S}}_w+{\varphi }_n{\vec{S}}_n+{\varphi }_s{\vec{S}}_s\right)$$

4 非结构网格

非结构网格更加灵活，然而也更加复杂。单元、面和角点都是顺序编号的，每个单元的邻近单元数目并不相等，每个单元的面数目也并不相等，甚至每个面的角点数目也并不相等，所以并不能从单元、面、角点的标识直接得到其邻近单元、面、角点的标识（与结构网格迥然不同），也就是需要存储这些额外的拓扑关系了。比如上图中，单元9的邻近单元为1、2、6、8、10、11，这是无法从单元9的标识中直接靠加减增量来获得的。

上图（a）为局部标识，（b）为整体标识，可见，单元C的整体标识为9，其邻近6个单元的整体标识为10、11、8、6、1、2，其所含6个面的整体标识为16、22、23、25、11、10，其所含6个角点的整体标识为21、22、21、14、13、12（角点标号图上并未标明，这里显然有俩21角点是不对的，仅仅作为示意不用太纠结了）。

那么梯度的计算，用局部标识表示出来便是
$$\nabla {\varphi }_C=\frac{1}{V_C}\left({\varphi }_{f_1}{\vec{S}}_{f_1}+{\varphi }_{f_2}{\vec{S}}_{f_2}+{\varphi }_{f_3}{\vec{S}}_{f_3}+{\varphi }_{f_4}{\vec{S}}_{f_4}+{\varphi }_{f_5}{\vec{S}}_{f_5}+{\varphi }_{f_6}{\vec{S}}_{f_6}\right)$$
用整体标识表示出来便是
$$\nabla {\varphi }_9=\frac{1}{V_9}\left({\varphi }_{f_{16}}{\vec{S}}_{f_{16}}+{\varphi }_{f_{22}}{\vec{S}}_{f_{22}}-{\varphi }_{f_{23}}{\vec{S}}_{f_{23}}-{\varphi }_{f_{15}}{\vec{S}}_{f_{15}}-{\varphi }_{f_{11}}{\vec{S}}_{f_{11}}-{\varphi }_{f_{10}}{\vec{S}}_{f_{10}}\right)$$
由于每个面的面积矢量仅计算了一次，所以对某个特定单元来说，需要判断该面积矢量是正还是负，一般的规定是，计算面的面积矢量时，其正方向为面左侧单元指向右侧单元，左单元的整体标识要小于右单元的整体标识，也常称左单元为owner单元，右单元为neighbour单元。所以，上式中出现了正负号，对于某个标识为k的单元，其梯度算式可写成如下形式
$$\nabla {\varphi }_k=\frac{1}{V_k}\left(-\sum _{n\leftarrow (f-nb(k))<k}{\varphi }_n{\vec{S}}_n+\sum _{n\leftarrow (f-nb(k))>k}{\varphi }_n{\vec{S}}_n\right)$$
其中$n\leftarrow (f-nb(k))<k$表示$k$单元的面$f$的另一侧单元的标识小于$k$，而$n\leftarrow (f-nb(k))>k$则标识大于k。

在OpenFOAM中，把内部面的面积矢量方向的正方向定义为从单元1到单元2，分别称为owner和neighbor，即所属单元和邻近单元，如下图所示

具体计算梯度的时候，如果是一个单元一个单元做循环就太麻烦了，所以做面循环比较划算，即：

1. 定义整场所有单元的梯度数组并让其初始值为0
2. 对于内部面做循环
   计算面通量flux_f=phi_f * Sf
   累加flux_f到面的owner单元梯度上，并累加-flux_f到面的neighbor单元梯度上
3. 对边界面做循环
   计算面通量flux_f=phi_f * Sf
   累加flux_f到面的owner单元梯度上
4. 对单元做循环
   单元梯度除以单元体积

5 几何信息

非结构网格除了上述拓扑信息外，还需要用到很多几何信息，比如单元体积、面面积、单元形心、面形心、面积矢量与面owner和neighbor单元形心连线的夹角，等。

5.1 单元类型

常见的3维单元有四面体、六面体、三棱柱、多面体等，构成单元的面有四边形、三角形和多边形，显然，由任意多边形构成的任意多面体是一种最通有的情形，而其他单元则是该单元的特例。

5.2 计算面的面积和面的形心

首先计算形心（centroid），如果是三角形，其形心和几何中心（geometric center）是重合的，很好计算，用其三个角点平均即可
$${\vec{x}}_{CE}={\vec{x}}_G=\frac{1}{3}\sum _{i=1}^3{\vec{x}}_i$$

然而，对于普遍意义上的多边形来说，其形心的计算就没那么简单了，需要先用一个点将该多边形剖分为许多小三角形，然后再用小三角形的形心来加权计算该多边形的形心，这个剖分点一般选取为几何中心即可，对于$k$边形来说，用其角点坐标平均可得其几何中心，即
$${\vec{x}}_G=\frac{1}{k}\sum _{i=1}^k{\vec{x}}_i$$
加权系数用小三角形的面积即$S_t$可，则多边形的形心为
$${\vec{x}}_{CE}=\frac{\sum _{t=1}^{subTriangles(C)}({\vec{x}}_{CE}{)}_t{S}_t}{\sum _{t=1}^{subTriangles(C)}{S}_t}$$
嗯，你说啥？三角形面积值$S_t$该咋算呢？

5.2.1 三角形的面积

用两条边矢量的叉乘除以2就好了，面积矢量为
$$\vec{S}=\frac{1}{2}({\vec{r}}_2-{\vec{r}}_1)\times ({\vec{r}}_3-{\vec{r}}_1)=\frac{1}{2}\left|\begin{array}{ccc}\mathbf{i}& \mathbf{j}& \mathbf{k}\\ x_2-x_1& y_2-y_1& z_2-z_1\\ x_3-x_1& y_3-y_1& z_3-z_1\end{array}\right|=S_x\mathbf{i}+S_y\mathbf{j}+S_z\mathbf{k}$$
面积幅值为
$$S=\sqrt{S_x^2+S_y^2+S_z^2}$$
那么如何判断面积矢量的方向是向外还是向内的呢？计算连接单元形心CE和面形心ce的矢量，然后再与$\vec{S}$做点积，如果大于0，则面积矢量是朝外的，如果小于0，则面积矢量是朝内的。

对于2维网格，面的面积相当于三角形拉伸了单位长度的控制体的体积，即
$$\begin{gathered} V=\frac{1}{2}|({\vec{r}}_2-{\vec{r}}_1)\times ({\vec{r}}_3-{\vec{r}}_1)|=\frac{1}{2}\left|\begin{array}{ccc}1& 1& 1\\ x_1& x_2& x_3\\ y_1& y_2& y_3\end{array}\right| \\ =\frac{1}{2}[x_1(y_2-y_3)+x_2(y_3-y_1)+x_3(y_1-y_2)] \end{gathered}$$
注意，若角点1、2、3是逆时针排序，则体积（面积）为正值，如果是顺时针排序则为负值，即右手准则，这也就是为什么2维问题的角点总是逆时针排序的。

5.2.2 单元体积和形心

对于多边形组成的多面体，其形心和体积的计算和前面的多边形很相似，即，先用几何中心点把该多面体分成很多小多棱锥体，然后计算每个棱锥体的体积（=底面积*高/3），而每个棱锥体的形心位于其底面形心和顶部角点连线的1/4处，将所有小多棱锥体的形心做体积加权平均，便得到了多面体的形心坐标。

用下面几何中心剖分多面体为很多小多棱锥体
$${\vec{x}}_G=\frac{1}{k}\sum _{i=1}^k{\vec{x}}_i$$
计算多棱锥体的形心
$$({\vec{x}}_{CE}{)}_{pyramid}=0.75({\vec{x}}_{CE}{)}_f+0.25({\vec{x}}_G{)}_{pyramid}$$
计算多棱锥体的体积
$$V_{pyramid}=\frac{1}{3}{\vec{d}}_{Gf}\cdot {\vec{S}}_f$$
加权计算多面体的形心
$$({\vec{x}}_{CE}{)}_C={\vec{x}}_C=\frac{\sum ^{subPyramids(C)}({\vec{x}}_{CE}{)}_{pyramid}{V}_{pyramid}}{\sum ^{subPyramids(C)}{V}_{pyramid}}$$

5.2.3 面加权系数

对于1维情况，已知面$f$两侧单元$C、F$的变量值，可插值算出面上的变量值，有
$${\varphi }_f=g_f{\varphi }_F+(1-g_f){\varphi }_C$$
其中
$$g_f=\frac{d_{Cf}}{d_{Cf}+d_{fF}}$$
对于高维系统就没这么简单了，几何加权系数的定义并不唯一，可以选择对应的体积作为加权系数，即
$$g_f=\frac{V_C}{V_C+V_F}$$
这种方式在有些情况下是不靠谱的，比如下图

显然，$fC=fF$，有$g_f=0.5$，然而若用体积加权，由于$V_C<V_F$，可得$g_f<0.5$，使得${\varphi }_f$偏于${\varphi }_C$了。

还有更加麻烦的事情，当$C,f,F$不共线的时候，情况就复杂得多了，这时候需要用到后续章节的非正交修正方法了，我们后面再谈。

如上图，一种处理方式是依据到面的法向距离来计算加权系数，即
$$g_f=\frac{{\vec{d}}_{Cf}\cdot {\vec{e}}_f}{{\vec{d}}_{Cf}\cdot {\vec{e}}_f+{\vec{d}}_{fF}\cdot {\vec{e}}_f}$$
其中${\vec{e}}_f$为面的单位法矢量
$${\vec{e}}_f=\frac{{\vec{S}}_f}{||{\vec{S}}_f||}$$

6 代码讲解

uFVM是仿照OpenFOAM，用matlab编写的CFD学习代码，其重在实现理论和后续深入理解OpenFOAM，而计算效率不及OpenFOAM（matlab的运行效率一向偏低，肯定比不过C++，但matlab的语言书写简便，便于理解），主要是便于学习如何实现理论，毕竟OpenFOAM理解起来确实太难了。

uFVM中，对于几何信息和拓扑信息的处理是在由cfdProcessOpenFoamMesh程序（2018年的V1.5版本中是cfdProcessGeometry函数）中进行的，该程序在cfdReadOpenFoamMesh程序读入OpenFOAM网格之后执行，读入的网格是用OpenFOAM的格式给出的，即网格由points，faces，owners，neighbours，boundaries文件给出。

cfdProcessOpenFoamMesh程序中，先计算face面的几何信息（面形心、面积、面积矢量），如下：

function cfdProcessGeometry
%--------------------------------------------------------------------------
%
%  Written by the CFD Group @ AUB, Fall 2018
%  Contact us at: cfd@aub.edu.lb
%==========================================================================
% Routine Description:
%   This function reads poly mesh files from "constant/polyMesh" directory
%   and stores them in the database ~domain
%--------------------------------------------------------------------------

% Get info
theNumberOfElements = cfdGetNumberOfElements;			% 单元数目
theNumberOfFaces = cfdGetNumberOfFaces;					% 面数目
theNumberOfInteriorFaces = cfdGetNumberOfInteriorFaces;	% 内部面数目
theFaceNodesIndices = cfdGetFaceNodeIndices;			% 面的角点标识（所有面的）
theNodeCentroids = cfdGetNodeCentroids;					% 角点形心（其实就是角点）坐标（所有角点的）
theElementFaceIndices = cfdGetElementFaceIndices;		% 单元的面标识（所有单元的）
owners = cfdGetOwners;									% 面的owner单元标识（所有面的）
neighbours = cfdGetNeighbours;							% 面的neighour单元标识（所有面的）

% Initialize mesh member arrays
elementCentroids = cfdVectorList(theNumberOfElements);	% 单元形心坐标（矢量列表，所有单元）
elementVolumes   = cfdScalarList(theNumberOfElements);	% 单元体积（标量列表，所有单元）
faceCentroids    = cfdVectorList(theNumberOfFaces);		% 面形心坐标（矢量列表，所有面）
faceSf           = cfdVectorList(theNumberOfFaces);		% 面面积矢量（矢量列表，所有面）
faceAreas        = cfdScalarList(theNumberOfFaces);		% 面面积（标量列表，所有面）
faceWeights      = cfdScalarList(theNumberOfFaces);		% 面权重系数（标量列表，所有面）
faceCF           = cfdVectorList(theNumberOfFaces);		% 面左右单元形心矢量（矢量列表，所有面）
faceCf           = cfdVectorList(theNumberOfFaces);		% 面左单元形心到面形心矢量（矢量列表，所有面）
faceFf           = cfdVectorList(theNumberOfFaces);		% 面右单元形心到面形心矢量（矢量列表，所有面）
wallDist         = cfdScalarList(theNumberOfFaces);		% 面owner单元形心到壁面的距离（标量列表，所有面），某些湍流模型里面用到壁面距离d
wallDistLimited  = cfdScalarList(theNumberOfFaces);		% 面owner单元形心到壁面距离的限制值（标量列表，所有面）

% Process basic face geometry
for iFace=1:theNumberOfFaces    						% 对所有面做循环（包含内部面 与 边界面）
    theNodesIndices = theFaceNodesIndices{iFace};    	% 该面（iFace号面）的角点标识
    theNumberOfFaceNodes = length(theNodesIndices);		% 该面总共有多少个角点
    
    % Compute a rough centre of the face    			% 计算面的大概中心（几何中心）
    local_centre = cfdVector(0,0,0);					% 所有角点坐标做平均处理
    for iNode=theNodesIndices
        local_centre = local_centre + theNodeCentroids(iNode,:);
    end
    local_centre = local_centre/theNumberOfFaceNodes;
    
    centroid = cfdVector(0,0,0);						% 初始化该面形心
    Sf = cfdVector(0,0,0);								% 初始化该面面矢量
    area = 0;											% 初始化该面面积
    
    % Using the centre compute the area and centoird of vitual triangles
    % based on the centre and the face nodes    
    for iTriangle=1:theNumberOfFaceNodes				% 面共划分的小三角形数目与角点数目相同
        point1 = local_centre;							% 三角形第1个点为面形心
        point2 = theNodeCentroids(theNodesIndices(iTriangle),:);	% 第2个点为第iTriangle个角点
        % 由于面角点轮换构成小三角形，即1-2,2-3,...,theNumberOfFaceNodes-1，所以最后一个点的第3个点的局部标识是1
        if(iTriangle<theNumberOfFaceNodes)	
            point3 = theNodeCentroids(theNodesIndices(iTriangle+1),:);
        else
            point3 = theNodeCentroids(theNodesIndices(1),:);
        end
        local_centroid = (point1+point2+point3)/3;			% 小三角形的形心等于重心等于几何中心
        local_Sf  = 0.5*cross(point2-point1,point3-point1);	% 叉乘计算小三角形的面积矢量
        local_area = norm(local_Sf);						% 取模得到小三角形的面积标量
        
        centroid = centroid + local_area*local_centroid;	% 小三角形形心与其面积相乘后累加
        Sf = Sf + local_Sf;									% 小三角形面积矢量的累加，获得整体面积矢量
        area = area + local_area;							% 小三角形的面积累加，获得面的面积
    end
    centroid = centroid/area;        % 除以面的整体面积，获得了面形心
    %
    faceCentroids(iFace,:) = centroid;	% 面形心存储到面形心列表中去
    faceSf(iFace,:)        = Sf;		% 面面积矢量存储到列表中去
    faceAreas(iFace)       = area;		% 面面积存储到列表中去
end

这个程序接下来是计算了单元的体积和形心，如下：

% Compute volume and centroid of each element
for iElement=1:theNumberOfElements		% 对所有单元做循环
    
    theElementFaces = theElementFaceIndices{iElement};	% 该单元（iElement单元）的面的标识列表

    % Compute a rough centre of the element	% 计算该单元的几何中心（粗略中心），即，面形心的平均值
    local_centre = cfdVector(0,0,0);
    for iFace=1:length(theElementFaces)
        faceIndex = theElementFaces(iFace);
        local_centre = local_centre + faceCentroids(faceIndex,:);
    end
    local_centre = local_centre/length(theElementFaces);
    
    % Using the centre compute the area and centoird of vitual triangles
    % based on the centre and the face nodes
    %
    localVolumeCentroidSum = cfdVector(0,0,0);	% 小多棱锥的（形心*体积）加和
    localVolumeSum = 0;							% 小多棱锥的体积加和
    for iFace=1:length(theElementFaces)			% 单元的所有面做循环
        faceIndex = theElementFaces(iFace);     % 取出第iFace个面的整体标识faceIndex     
        
        Cf = faceCentroids(faceIndex,:) - local_centre;	% 面形心到体中心的矢量
                
        faceSign = -1;			% 判断该面面积矢量的正负号
        if iElement==owners(faceIndex)	% 如果该面的owner是该单元，则其面积矢量为正值，不然为负
            faceSign = 1;
        end
        
        local_Sf = faceSign*faceSf(faceIndex,:);	% 计算该面的面积矢量（向外为正）
        
        % calculate face-pyramid volume
        % 计算小多棱锥的体积：底面面积矢量 与 底面形心到顶点矢量 的点积 / 3
        localVolume = dot(local_Sf,Cf)/3;
        
        % Calculate face-pyramid centre % 计算小多棱锥的形心
        localCentroid = 0.75*faceCentroids(faceIndex,:) + 0.25*local_centre;
        
        %Accumulate volume-weighted face-pyramid centre	% (形心*体积) 加和
        localVolumeCentroidSum = localVolumeCentroidSum + localCentroid*localVolume;
        
        % Accumulate face-pyramid volume	% 小多棱锥体积加和，得到多面体的体积
        localVolumeSum = localVolumeSum + localVolume;
    end
    % 体积平均，得到多面体形心，并存储到单元形心列表中
    elementCentroids(iElement,:) = localVolumeCentroidSum/localVolumeSum;	
    elementVolumes(iElement)     = localVolumeSum;	% 多面体体积存储到单元体积列表中去
end

这个程序接下来做的事情是计算内部面形心f的own和neighbour单元形心C和F相关量，即矢量CF，矢量Cf，矢量Ff，几何权系数gC。以及边界面的矢量CF、矢量Cf，几何权系数gC=1，等

% Process secondary Face Geometry
for iFace=1:theNumberOfInteriorFaces   
    n = cfdUnit(faceSf(iFace,:));    
    %
    own = owners(iFace);
    nei = neighbours(iFace);        
    %
    faceCF(iFace,:)    = elementCentroids(nei,:) - elementCentroids(own,:);
    faceCf(iFace,:)    = faceCentroids(iFace,:)  - elementCentroids(own,:);
    faceFf(iFace,:)    = faceCentroids(iFace,:)  - elementCentroids(nei,:);
    faceWeights(iFace) = (faceCf(iFace,:)*n')/(faceCf(iFace,:)*n' - faceFf(iFace,:)*n'); % -(dot(faceFf(iFace,:),n))/(-dot(faceFf(iFace,:),n)+dot(faceCf(iFace,:),n))
end

for iBFace=theNumberOfInteriorFaces+1:theNumberOfFaces
    n = cfdUnit(faceSf(iBFace,:));    
    %
    own = owners(iBFace);       
    %
    faceCF(iBFace,:)        = faceCentroids(iBFace,:) - elementCentroids(own,:);
    faceCf(iBFace,:)        = faceCentroids(iBFace,:) - elementCentroids(own,:);
    faceWeights(iBFace)     = 1;
    wallDist(iBFace)        = max(dot(faceCf(iBFace,:),n), cfdSMALL);	% 边界单元中心到壁面的距离
    wallDistLimited(iBFace) = max(wallDist(iBFace), 0.05*cfdMag(faceCf(iBFace,:)));    
end

最后是把这些算好的几何信息存储到网格信息中去

% Save and Store
mesh = cfdGetMesh;
%
mesh.elementCentroids = elementCentroids;
mesh.elementVolumes   = elementVolumes;
mesh.faceCentroids    = faceCentroids;
mesh.faceSf           = faceSf;
mesh.faceAreas        = faceAreas;
mesh.faceWeights      = faceWeights;
mesh.faceCF           = faceCF;
mesh.faceCf           = faceCf;
mesh.faceFf           = faceFf;
mesh.wallDist         = wallDist;
mesh.wallDistLimited  = wallDistLimited;
%
cfdSetMesh(mesh);