目錄
第8 章 有限體積法.2
§8.1 有限體積法的基本概念.2
8.1.1 幾個術語2
8.1.2 控制體積的選擇.3
8.1.3 結構與非結構網格.5
§8.2 有限體積法構造7
8.9-1 方程離散及基本格式7
8.9-2 物理特性要求.11
8.9-3 迎風型通量格式14
8.9-3 TVD 格式18
§8.3 非結構網格上的有限體積法23
8.3.1 基本方程.23
8.3.2 離散基本思路24
8.3.3 數值通量近似.25
第9 章 在水力學問題中的應用.29
§9.1 滲流問題中的應用.29
9.1.1 飽和—非飽和地下水運動基本控制方程29
9.1.2 方程的離散31
9.1.3 算例【陳揚 碩士論文】33
§9.2 二維明渠非恒定流計算.38
9.2.1 水流基本方程38
9.2.2 控制方程的離散39
9.2.3 計算實例.53
§9.3 三維紊動分層流計算.64
9.3.1 紊動分層流基本方程64
9.3.2 紊流模型及控制方程離散.65
9.3.3 壓力校正法66
9.3.3 邊界條件69
9.3.4 鹽度引起的負浮力流動71

第8章 有限體積法
有限差分方法是從描述各種物理現象的基本微分方程出發構造離散方程的，前文已經對其作了翔實、周密的論述。該部分將從基礎算法入手分析介紹在計算流體力學界廣為應用的有限體積法。基于有限體積法的實用算法在計算流體力學、計算傳熱學等領域得到了飛速發展 [1-3]。在水力學諸多問題，如水流物質輸運模擬，水工水力學模擬以及潰壩洪水波演進等水流模擬中也得到了廣泛應用。
§8.1 有限體積法的基本概念
有限體積法，又稱為有限容積法，它正是從物理量守恒這一基本要求出發提出的。這也是其受計算流體力學界廣為稱道和喜歡之處。其以守恒型的方程為出發點，通過對流體運動的有限子區域的積分離散來構造離散方程。有限體積法有兩種導出方式，一是控制容積積分法，另一個是控制容積平衡法。不管采用哪種方式導出的離散化方程，都描述了有限各控制容積物理量的守恒性，所以有限體積法是守恒定律的一種最自然的表現形式。
該方法適用于任意類型的單元網格，便于應用來模擬具有復雜邊界形狀區域的流體運動；只要單元邊上相鄰單元估計的通量是一致的，就能保證方法的守恒性；有限體積法各項近似都含有明確的物理意義；同時，它可以吸收有限元分片近似的思想以及有限差分方法的思想來發展高精度算法。由于物理概念清晰，容易編程；有限體積法成為了工程界最流行的數值計算手段。
8.1.1 幾個術語
在進行數值計算時，要把計算區域劃分成一系列互不重疊的離散小區域，然后在該小區域上離散控制方程求解待求物理量。在有限差分法中只涉及到網格節點的概念，而有限體積法因為物理解釋需要，形成了以下幾個常用幾何要素的相關名詞。
控制體積（control volume）：圖中陰影部分，方程積分離散時的小體積單元（二維為面積單元）。
單元（cell）：控制體積的中心，常用形心來表示，為待求物理量的幾何位置。圖中用空心園來表示，如點P、W、E、N、S 等。常用單元來代表整個控制體積。以下如若不作特殊說明，則用“單元”來代表控制體積。
網格線（grid lines）：用來分割計算區域內各控制體積的交錯曲線簇，如圖中折線N1N2、N2N3、N3N4、N4N1等
網格節點（nodes）：網格線之間的交點， 圖中用黑圓點來表示，如N1、N2、N3、N4等單元界面（cell faces）：相鄰兩個控制體積間的公共面(二維則可以認為是公共邊，這樣看起來就和網格線一致了，但是要注意這不是同一個概念)，圖中用小寫字母e、n、w、s 表示。通常定義e、n、w、s 幾何位置位于交界面的形心點，二維則認為在公共邊的中心點。
8.1.2 控制體積的選擇
當你開始用有限體積法模擬流體流動時，而且劃分好網格后，你必須選定控制體積的形成方式。目前，常用的有兩種方法：單元中心方式（cell-centered）和頂點中心方式
（vertex-centered）。另外一些學者還發展了由兩種方式綜合形成的混合方式。根據問題的特點和要求，不同的變量可以采用不同的控制體積，因此又產生了交錯網格和同位網格的稱謂，這里不再深入介紹，讀者根據需要可以參考相關文獻[1-3]。
前面我們我們在討論有限體積法的術語時，已經看到了單元中心控制體積的形成方式。
這里再作一個說明并對兩者作簡單的比較。圖8-2 為兩種控制體積選擇方式示意圖。陰影部分表示單元的控制體積。空心圓點表示單元，實心圓點表示網格節點，在頂點中心方式中，單元和網格節點重合。

（1）單元中心方式
單元位于控制體積的中心，即將單一的網格單元作為控制體積，網格單元互不重疊。目前用非結構網格做流動計算多使用單元中心格式,它有利于處理陸地和給定流量邊界條件。
（2）單元頂點方式
以網格節點為中心，它的形成有多種方式。其常用的構成方式是由以該節點為頂點的格子的形心以及各共頂點的網格線中心點的一系列連接線段所構成一個多邊形控制體積。也可以由環繞該節點的一組格子組成控制體。
因此計算同樣多的變量，單元中心方式變量布置簡單直觀，易于處理邊界條件和保持離散的守恒性，而且需要的網格數要比單元頂點方式少得多，可節省計算時間。
8.1.3 結構與非結構網格
1 結構網格
結構網格，指的是一類具有一定的分布特征，可以用相應的行列關系來順序描述的網格。
常用的矩形網格、曲線網格以及塊結構網格等。
圖 8-3 結構網格
矩形網格是FDM 或FVM 最為常用的網格系統。它不對原始的控制方程（直角坐標）作任何坐標變換，在邊界處采用簡單的階梯形網格近似復雜的邊界。或者，用規則矩形網格覆蓋整個包括陸邊界在內的計算區域，對陸邊界及陸域處的網格作特殊的技術處理，如所謂的凍結法，即令某些控制體積“凍結”，或者說，“堵塞”某些控制體積，則剩下的“活動”的控制體積可形成所需要的不規則區域[1]，“凍結”亦即令該控制體積中的因變量在運算過程中始終保持不變。常采用大數值源項和大粘性系數法實現凍結區內流速恒為零。盡管這類對規則網格作特殊的技術處理的方法對邊界概化過于嚴重，邊界計算值過于粗糙，陸邊界附近會形成虛假的曲折水流。但是，其網格生成方便，計算方法簡單等優點使該方法得到了大量的應用。
在計算天然水域時，利用規則矩形網格階梯近似邊界不僅改變了原始物理邊界附近的流態，而且增加了邊界條件設置的復雜性；計算精度因此也大受影響，尤其在近邊界高梯度區域。自然地，研究人員引入廣義曲線坐標變換的思想。即在物理計算區域內構筑曲線網格，使得網格曲線與所求解區域的邊界相重合，而后利用坐標變換將復雜的物理域變換到規則的計算空間內；在規則計算區域上離散求解變換后的控制方程，將得到的解投影至原物理區域得到數值解；或者在曲線網格上直接應用原始因變量求解。除將復雜的物理空間變換成為易于求解的規則網格之外，曲線坐標變換方法還能式計算網格點與物理域的活動網格節點相對應，從而可適應非恒定流動邊界的模擬。目前常用的生成曲線網格的方法有微分方法和代數方法。微分方法常用求解橢圓型方程或雙曲型方程實現。其中，尤以Thompson 為代表提出
的利用Poisson 方程生成貼體曲線網格的方法最為著名和流行。該類型方法可通過調整源項控制網格的分布，在梯度大的區域人工或自動加密網格。因此，該方法提出后，得到了廣泛的發展和應用。

曲線網格應用遺留眾多需要解決的問題。首先，大多數網格的生成方法只提供了離散點的變換，而不給出解析函數形式的變換關系。使用不光滑的網格時，對變換關系的差分近似會造成了很大的數值誤差，甚至會導致不切實際的值。其次，如果網格嚴重偏離正交性，就會極大損壞原有的迭代方法的收斂速率。再次，因變量的選擇也須謹慎考慮。在曲線網格中，可取原始笛卡爾坐標系變量或曲線坐標系中沿網格方向的協變量兩種作為因變量。
滿足正交性的網格具有優良的特性。在正交曲線坐標系下，變換后的某些項消失；控制方程更為簡潔；數值穩定性和解的精度高；計算量和收斂速度也有了很大的改善。
但是，正交曲線網格極大的受制于邊界的幾何形狀，對于天然水域，由于邊界極其復雜，
很難生成完全正交的曲線網格，特別對于三維問題，正交曲線網格幾乎是不可能實現的。因此，在實際工程水流計算上較少采用。
正是由于復雜邊界下生成正交坐標系的困難，而非正交曲線網格不僅相對較容易生成，而且可隨意調節網格的疏密，在高梯度區域人工或自動加密網格。因此，研究人員退而求其次，放棄完全正交性的要求，轉而尋求非正交曲線網格生成和數值處理。如生成邊界正交曲
線網格，或生成盡量接近正交的曲線網格，以及探求非正交的光滑或不光滑或網格的求解方法[2]等。為充分利用網格正交性的優點，常希望盡可能生成接近正交的網格，至少在接近
邊界高梯度區域，以恢復數值近似損失的精度。為此，Sorenson[11]開發了一種能控制網格疏密和角點的邊界正交曲線網格生成的微分方法。
2 非結構網格
八十年代以來，為了更好地擬合自然邊界，非正交網格上的數值計算開發和應用得到了飛速的發展。
為了更好的說明非結構網格，可以先看圖8-4。
圖 8-4 非結構網格
從圖中可以看出，非結構網格中單元格分布不再規則一致，其位置很難再憑借行列索引關系確定。非結構網格可以采用任意形狀的單元格，單元邊的數目也無限制，彌補了結構化網格不能夠解決任意形狀和任意連通區域的網格剖分的缺欠。實用上，為簡化編程以及貼合邊界要求，一般應用四面體、六面體（二維為三角形、四邊形）等已經足夠擬合天然水域邊界。非結構網格最重要的一個特征是控制方程離散得到的代數方程的系數矩陣不再是結構網格下有規律的對角結構；若用對角形式存儲，其帶寬只能通過適當的布置單元編號順序來減少。非結構網格原則上可應用于任何類型的數值方法，包括FDM，但FEM 和非結構網格上的FVM 算法更成熟，應用更廣。

非結構網格最早用于FEM，但流體流動是高度非線性問題，而且FEM 計算量較大，這些問題使得基于FEM 的非結構網格技術未能在對流問題為主的地面水流（如淺水流動，水波運動等）計算上得到重視。八十年代以來，基于FVM 的非結構網格技術在空氣動力學得到了廣泛的發展和應用。九十年代開始一些專家學者根據淺水流動特征，將這些算法引入到計算淺水動力學中，并在模擬涌潮，潰壩等水力計算難題上取得了成功[5]。
與此同時，粘性流動的非結構網格FVM 模擬也開始出現。并在20 世紀90 年代中后期掀起了研究高潮[6]。作為全球計算流體力學軟件供應商和技術服務商的Fluent 公司已經將最新的非結構網格研究成果集成，實現了研究成果的商業化。
非結構網格能很好地模擬自然邊界及水下地形，利于邊界調節的實現；便于控制網格密
度，易作修改和適應性調整；網格生成有眾多富有成效的方法和自適應技術[9]，比曲線網格更易得到高質量的單元格。但是非結構網格應用也帶來一系列需要解決的問題：單元格排列不規則，須建立相應的數據結構存儲單元格信息；控制方程離散得到的代數方程的系數矩陣是高度稀疏的非對角型矩陣，需要尋求合適的存儲方式及解法；隱格式較難實現，粘性項處理困難；數值解后處理工作量大；二階非結構FVM 較易實現，若要擴展到高階格式，則
需花費較大的代價。
§8.2 有限體積法構造
8.9-1 方程離散及基本格式
1 基本思路
與有限差分法不同，有限體積法是基于守恒型的控制方程。這里以一維對流方程為例，
說明有限體積法離散的基本思路。

為推導離散化的方程，我們必須事先對計算區域進行網格劃分，然后在各個單元上進行積分平均。假設我們對一維計算域劃分得到一個網格系統，如圖 所示的為該網格系統中某單元鄰接關系示意圖。重點考察單元P，單元E 和W 分別為它的兩個相鄰單元。虛線為單元界面，圖中用小寫字母表示，界面e、w 包圍形成的陰影部分即為單元P 的控制體積。控
制體積實際是三維體的的概念，只不過在考慮一維控制體積時，實際上把y、z 兩個坐標方向假設成單位厚度，這樣在一維問題中看起來成為了線的概念。

2 狀態變量分布近似
從前面可以知道，用有限體積法推導離散化方程時，必須確定物理量的局部分布，這是離散過程極為重要的一步，不僅關系到守恒性是否保持，而且關系到算法精度，這是由于不同的狀態變量分布對應了不同的離散格式。

值得指出的是，在用有限體積法計算時，和有限差分法一樣，方程的解是用單元節點上離散點值構成的，而不關心單元間的狀態變量是怎么變化的，也就是不關心解的分布。這和有限單元法中一旦選定了分布曲線，就確定了狀態變量的分布函數是不同的。盡管我們在離散方程的時候要假定單元分布曲線，但是這只是為了推導公式時計算積分近似而采用的一些輔助關系式。一旦離散化方程推導出來了，就可以不用再管這些分布曲線近似。這種觀點使我們在采用何種分布曲線近似方法有完全的自由。在積分離散時，根據數值模擬的需要，對控制方程中的每一項都可以采用不同的分布曲線來近似單元界面上的狀態變量，而不必要追求近似假設的一致性。例如，我們用一階迎風格式離散對流擴散方程時，對流項中因變量是用階梯型分布近似，而擴散項中也取為階梯分布的話，則根本導不出離散式，至少要有二階的分段線性分布才行。上面已經提到，分布曲線近似關系到算法的精度，實際上，有限體積法中，每一種不同的狀態變量分布近似，體現了不同的離散格式的幾何構造方法。如對流問題中，一階迎風格式為階梯型分布，中心格式采用的分段線性分布，更為高階的格式則需要更多的單元加入，形成復雜的狀態變量分布近似關系，例如對流項二次迎風插值QUICK，分段拋物插值PPM 等。
3 基本格式
利用這些交界面表達式可以得到具體的離散化的代數方程，這里不在敘述。在這里，我們還可以看出，有限體積法離散的思路和有限差分法的有著本質的區別，有限差分法是應用Taylor 展開式近似微分控制方程中的導數，而有限體積法則是通過積分將導數項化成交界面的狀態變量的表達式，然后，代入根據局部分布近似導出的交界面值得到完整的離散化的代數方程。
對于中心格式和一階迎風格式的特點和有限差分法的一致。中心格式精度較高，具有二階精度，但是不能反映對流傳輸機制。對非線性問題，常會出現非線性不穩定，為了消除該影響，需要添加人工粘性。一階迎風格式穩定性好，但是數值耗散大。值得強調的是，一階迎風格式離散思想中蘊涵的物理本質，即波動傳波的信息主要來自上游方向，是我們構造高分辨率格式的基礎。

8.9-2 物理特性要求
1 守恒性
如果對一個離散方程在定義域的任一有限空間內作求和的運算（相當于連續問題中對微分方程作積分），所得的表達式滿足該區域上物理量守恒的關系時，則稱該離散格式具有守恒特征。
我們知道，根據流體運動的不同特征，可以采用拉格朗日（Lagrange）法和歐拉（Euler）法來描述流體的運動。其中歐拉法著眼于研究固定空間位置上各流體質點運動要素的空間變化情況，在推導流體運動的控制方程時，歐拉法被廣為采用。讓我們回顧一下連續方程的推導過程，首先在流場中劃定一個固定的空間區域作為流體運動的控制體積，考察流體流入、流出該控制體積的情況，利用質量守恒原理，可導出流體運動的連續性方程。其他描述流體中物質（能量）輸運的控制方程也可根據各自的守恒律推導出來。在數值計算中，物理量的守恒性也必須得到保證，否則會導致非物理解的產生。因此守恒性成了數值求解過程中非常重要的一個方面。
有限體積法正是從物理量守恒這一基本要求出發提出的。它有兩種導出方式[3]，一是控制容積積分法，另一個是控制容積平衡法。不管采用哪種方式導出的離散化方程，都描述了有限各控制容積物理量的守恒性。就像微分方程是表示關于無窮小控制體積內物理量的守恒一樣，帕坦卡[1]（S.V. Patankar）非常生動的對此作了注解：如果把我們想象成處于微積分以前的時期，控制容積方程就會成為我們表達守恒原理的唯一途徑。有限體積法的離散化方程滿足了單個控制體積的平衡，當然在整個計算區域內，諸如質量、動量等物理量的積分守恒也就都能精確得到滿足。無論在數值計算中采用巨大數目的細網格和少數的粗網格，數值解也照樣顯示準確的積分平衡。有限體積法的離散思想自動滿足守恒定律，如質量守恒，動量守恒，能量守恒等等。所以有限體積法是守恒定律的一種最自然的表現形式。
為了得到守恒型的離散方程，需要滿足一定的條件。這里仍以前文導出的全離散的守恒型顯格式離散方程（8-12）為例。一般對l + k +1個單元（稱為節點模板），守恒型數值通量

守恒格式可以自動算出激波并給出其正確的位置，因此這種數值方法稱為激波捕獲法（shock-capturing method）。如果使用非守恒格式，則需要在有關位置使用局部Rankine—Hugoniot 條件，以獲得準確的激波位置，這種方法稱為激波擬合法(shock-fitting method)。
2 迎風性
在有限差分法部分中的基本理論部分，我們已經深入的討論了雙曲型方程（對流）的特征線概念。特征值體現了微分方程的解（或者說擾動、信息）的傳播方向。它既有很強的數學意義，更反映了問題的物理機制，表達了波動、能量等的傳播方向。這啟示我們在理論分
析和數值方法的設計上，應該充分考慮這一物理特性，與之相適應，而不可反其道而行之。根據特征（信息）傳播方向，設計得到的數值方法可統稱為特征型或迎風型格式。例如前文
提到的傳統的一階迎風格式，盡管非常簡單，卻是蘊含了最原始的的迎風設計思想，即在計算單元界面處的數值通量時，應根據特征傳播的方向性，對于向右傳播的分量，應該使用左
邊的單元值來計算(因為影響來自左邊)；對于向左傳播的分量，應該使用右邊的值來計算(因為影響來自右邊)。
當前，這種根據信息傳播方向設計數值格式的迎風思想已經成為了流體數值計算發展的重要基石，它不斷地啟發研究人員構造出各種各樣的離散格式，如特征迎風格式，矢通量分裂格式，通量差分裂格式等。并在這些基于迎風思想的格式基礎上，發展了分門別類的高階、高分辨率格式，如MUSCL 格式、QUICK 格式、ENO（WENO）格式等等。
3 TVD 性

既要利用格式的高階精度，又能使格式具有TVD 性質，從而求得高精度的物理上有意義的解，是研究人員的目標，為此，過去的二十幾年中，各種二階（甚至更高階）的TVD 格式得到了大量的發展。這類高階TVD 格式開發應用成為了計算流體力學界最前沿的研究。
8.9-3 迎風型通量格式
對于雙曲型問題，不同特征分量傳播的方向不同。考慮傳播方向的影響從物理上分析更為合理。我們對傳統的一階迎風格式已經作了介紹。這里我們將介紹20 世紀80 年代發展起來的通量向量分裂格式（Flux Vector Splitting），通量差分裂格式（Flux Difference Splitting）以及Godunov 法和近似黎曼解算子法（Approximate Riemann Solvers）。
1 矩陣分裂
實際流動的控制方程是一個非線性系統，各個變量之間相互耦合。我們這里考慮一維的雙曲型方程組
2 通量向量分裂格式（FVS）
Steger 和Warming 根據矩陣分裂方法，于1981 年首次提出了一種通量向量分裂格式。
最早該格式應用于空氣動力學計算上，20 世紀90 年代，開始擴展到計算水力學領域。

在把通量向量格式應用到水力學中時，仍需要注意方程組的特性。因為水流運動的控制方程——圣維南（Saint Venant）方程組的通量項不具有齊次性，需要補充一個冗余的能量方程得到具有齊次通量的方程組，然后就可以按照前面的方法建立離散化方程。
3 通量差分裂格式（FDS）
通量向量分裂格式是將通量函數按不同特征方向進行分裂，然后積分離散。與通量向量分裂格式不同，通量差分裂格式是對任意區間上通量函數的差進行分裂的。
4 Godunov 格式和近似黎曼解法

8.9-3 TVD 格式
在前文討論中，我們已經知道，滿足TVD 特性的格式求出的數值解才是真正符合實際的物理解。然而，一階迎風格式盡管具有TVD 性，但是數值耗散過大，將本來應該很陡峭的激波給抹得很平．為了在解光滑的區域獲得較高精度的解，并且在激波等間斷區域獲得沒有數值振蕩或者是數值振蕩可以令人接受的較陡峭的數值解，人們開始研究高分辨格式（High Resolution schemes）。1983 年Harten 在他的論文中提出了高分辨率和總變差減少的概念，首先證明了計算格式具有TVD 性質的充分條件，具體構造了修正通量的高分辨率（二階）TVD 格式。他的成果在數值算法的研究上具有十分重要的意義，，為構造和分析高分辨率格式提供了有力的工具。隨后，各種各樣的TVD 格式如雨后春筍般的涌現出來，因其具有對間斷的自動高分辯識別能力，而得到了廣泛的應用。這種格式具有兩個特點：（1）對標量非線性方程及常系數雙曲型方程組，格式的解是總變差遞減的（Total Variation Diminishing）；（2）與守恒律和熵不等式是相容的。TVD 格式提高了對激波的捕捉能力，但其不足之處在極值點上格式只有一階逼近精度。

§8.3 非結構網格上的有限體積法
前面主要對有限體積法基本概念和離散格式作了介紹。在這節中，我們將介紹二維非結構網格上的有限體積法，以便于應用它來模擬自然中復雜區域內的流動及物質輸運現象。本節只對算法的空間離散進行討論，因為時間的離散和有限差分法一致，因此，不在專門介紹。

標簽： 點擊： 評論: