學習自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 uFVM

OpenFOAM是強大高效的開源代碼，而uFVM則側重教育學習（便于理解卻喪失效率），本章著重講解OpenFOAM的網格文件格式，以及uFVM的網格數據結構是如何架構的，可見uFVM與OpenFOAM的實現細節是非常類似的，可以作為學習OpenFOAM的先導。

1 uFVM

1.1 OpenFOAM測試算例

uFVM代碼是直接讀入OpenFOAM的算例配置文件的，所以先要理解OpenFOAM的算例是如何設置的，下圖展示了一個名為cavity的算例，其設置文件全部放在cavity文件夾下面。

在cavity文件夾下，有三個文件夾：

其中文件夾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中的單元是任意多面體，而每個面可以由任意多邊形構成，所以非常靈活通用，當然結構網格是這種多面體網格的一種特例。

OpenFOAM中只處理3維網格，如果是2維問題，把它沿著展向拉伸一層網格讓它變成3維網格，同時把展向兩端的邊界面定義成empty就可以了。

points

points文件存放的是角點的坐標$(x_i,y_i,z_i)$列表，注意，這些角點，既是單元的角點，也是面的角點。標識為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 % 節點總數為1074caseDirectory: 'elbow' % 算例名字elblownumberOfFaces: 3290 % 面總數為3290numberOfElements: 918 % 單元總數為918 faces: [1x3290 struct] % 3290個面的信息列表numberOfInteriorFaces: 1300 % 內部面總數為1300boundaries: [1x6 struct] % 6個boundary patch的信息列表numberOfBoundaries: 6 % 邊界的數目為6（這個量好像跟下面的重復了！）numberOfPatches: 6 % 邊界patch的數目為6elements: [1x918 struct] % 918個單元的信息列表numberOfBElements: 1990 % 邊界單元數目為1990numberOfBFaces: 1990 % 邊界面的數目為1990，這個與邊界單元的數目是一致的

這個網格可以用cfdPlotMesh函數將其可視化

在nodes列表中，每個元素為一個結構體，代表著第i個節點的信息（注意不同于c++，matlab列表的下標是從1開始算的），那么第1個節點的信息如下

n1= m.nodes(1) n1 =centroid: [3x1 double] % 節點形心坐標，即節點的坐標x,y,zindex: 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,zSf: [3x1 double] % 該面的面積矢量Sx,Sy,Szarea: 5.3046 % 該面的面積SCN: [3x1 double] % 該面所屬單元形心到該面形心的距離矢量CfgeoDiff: 4.5940 % ？面幾何擴散系數 gDiff_f = Ef / CF，見第8章T: [3x1 double] % 面所屬單元和鄰近單元形心之間的距離矢量CF？（感覺更像是Tf=Sf-CF為非正交修正矢量，見第8章內容）gf: 0.4226 % 面插值中的幾何權重系數gfwalldist: 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: 300iNeighbours: [278 302 590]iFaces: [407 435 436 2053 2054]iNodes: [283 820 679 142 290 827]volume: 1.9083faceSign: [-1 1 1 1 1]numberOfNeighbours: 3centroid: [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: 1301iOwner: 1iNeighbour: -1centroid: [3x1 double]Sf: [3x1 double]area: 3.7510CN: [3x1 double]geoDiff: 5.6264T: [3x1 double]gf: 1walldist: 0.6667iOwnerNeighbourCoef: []iNeighbourOwnerCoef: []

這些信息中可以發現邊界面與內部面的不同之處，即邊界面的neighbour單元是沒有的（標識為-1），邊界面的幾何權重系數gf為1。

那么，如果要對對某個邊界patch做循環的話，該如何處理呢？對這個邊界patch的起始面和終止面循環就好了，比如要對boundary patch的第2個patch的面做循環，可以這樣子

theMesh = cfdGetMesh; iPatch = 2; iBFaces = cfdGetFaceIndicesForBoundaryIndex(iPatch) for iBFace = iBFacestheBFace = 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.meshans = nodeCentroids: [1074x3 double]numberOfNodes: 1074faceNodes: {3290x1 cell}numberOfFaces: 3290owners: [3290x1 double]numberOfInteriorFaces: 1300numberOfBFaces: 1990neighbours: [1300x1 double]numberOfElements: 918numberOfBElements: 1990cfdBoundaryPatchesArray: {6x1 cell}numberOfBoundaryPatches: 6closed: 0elementNeighbours: {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 01 0 01 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:numberOfElementstheElement = theMesh.elements(iElement) % 取出第iElement個單元的拓撲幾何信息phi(iElement) % this is field phi at centroid of element iElement % 取出第iElement個單元形心上存放的物理量phi值.... end

如果要對邊界單元做循環處理，如下

% 對于邊界單元（只有單元場才有邊界單元的概念）做循環 % 循環范圍從 整體單元數目+1 到 整體單元數目 + 邊界面數目（邊界面數目等于邊界單元的數目） for iBElement = numberOfElements + 1: numberOfElements + numberOfBFacesphi(iBElement) = 0; % 將邊界單元上的phi值設為0 end

1.5.2 面循環遍歷

面的編碼是先做內部面編碼，然后再依次對邊界面做編碼，因此，若對內部面做循環，只需要從1到numberOfInteriorFaces循環就好了，即

for iFace=1:numberOfInteriorFacestheFace = theMesh.faces(iFace) % 取出第iFace個面的幾何拓撲信息 end

如果要對邊界面做循環，則

for iBFace= numberOfInteriorFaces+1:numberOfFacestheBFace = 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…

總結

以上是生活随笔為你收集整理的FVM in CFD 学习笔记_第7章_OpenFOAM和uFVM中的有限体积网格的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。