“數值傳熱學”中貼體網格與非結構化網格上的有限容積法對比教學

宇 波， 禹國軍， 韓東旭， 李敬法， 王 藝， 孫東亮

（1.北京石油化工學院機械工程學院，北京 102617；2.上海海事大學商船學院，上海 201306；3.中國石油大學（北京）機械與儲運工程學院，北京 102249）

0 引言

流動與傳熱數值計算中多涉及不規則的計算區域，因此不規則區域的流動與傳熱數值計算理論是流動傳熱仿真軟件開發的重要基礎理論，其中貼體網格和非結構化網格有限容積法是核心內容。然而，教學經驗表明，初學者不易準確區分貼體網格和非結構化網格，認為不規則的網格都是非結構化網格，對兩者的有限容積法計算理論（如控制方程的離散與求解等）更是困惑。此外，由于兩種網格的定義、內涵、有限容積理論等知識點在教材及教學過程中又較為分散，學生難以對非規則區域有限容積法體系進行整體把握，因此在分別完成非結構化有限容積法各教學內容之后，進行一次完整理論體系涉及的各知識要點的綜合辨析教學將會是一個有效的手段。綜合辨析教學可以讓學生更容易厘清兩者的異同、掌握兩者對應的控制方程離散方法、理解兩者對物理問題的適應性，實現對該理論體系理解和應用能力的升華，從而為將來進行流動傳熱數值仿真軟件的自主研發奠定基礎。

1 貼體網格和非結構化網格的差異

推薦采用圖1 所示的兩套網格系統來區分貼體網格和非結構化網格。為了講清楚貼體網格的定義“貼體網格是一種通過建立物理域邊界點與計算域邊界點的映射關系，從而得到物理域與計算域內部節點間的對應關系，將不規則的物理域轉化成規則的計算域，然后由計算域上規則的結構化網格根據對應關系轉換到原物理域上的網格”［1］。以圖1（a）來進行解釋：為了生成圖1（a1）所示的不規則六邊形區域的網格，將該不規則區域的邊界映射到圖1（a2）所示的正方形邊界上，然后生成計算域上的規則的結構化網格，最后根據邊界的映射關系將計算域上的網格轉換到原物理域上，得到圖1（a1）所示的貼體網格。通過該例，可以讓學生更容易理解為什么貼體網格本質上是結構化網格：貼體網格是結構化網格的一個數學變換，其具有統一的拓撲結構。對于非結構化網格，建議從網格生成的過程來理解其定義：即非結構化網格是一種直接生成適應物理域幾何形狀的不規則的網格，網格點生成的順序不存在幾何上的連續性，即相鄰的網格點的生成順序不一定是連續的。以圖1（b）所示的非結構化三角形網格為例進行說明：對于圖1（b）所示的某個三角形單元，無法直接給出其鄰點的關系，因為網格點的生成順序是無規律的，因此鄰點的關系須在計算時進行判斷。通過這樣的講述，學生就很容易理解非結構化網格的本質特征“網格不具有統一的拓撲結構”。通過以上兩套網格的對比給出結論：貼體網格和非結構化網格的本質差異就在于結構化網格節點具有統一的拓撲結構，而非結構化網格節點拓撲結構雜亂。

圖1 貼體網格和非結構化網格

2 貼體網格和非結構化網格的有限容積法實施過程差異

課堂教學中可采用表1 來對比介紹基于貼體網格和非結構化網格的有限容積法實施過程的差異。通過該表，讓學生從宏觀上理解兩種方法的差異及造成差異的原因，即：由于網格坐標系的差異，兩種方法所采用的控制方程的形式、控制方程及邊界條件的離散均存在較大的差異?；谫N體網格的有限容積法在與物理域不同的計算域上進行控制方程的離散求解，因此控制方程及邊界條件需要轉化到計算域對應的坐標系，而非結構化網格則直接在物理域上進行離散求解，無須進行坐標變換，但由于非結構化網格不具有固定的鄰點關系，不易直接計算偏導數，因而采用矢量形式的控制方程［1-3］。

表1 基于兩種網格系統的有限容積法整體差異對比

3 貼體網格和非結構化網格有限容積法控制方程的異同及變換要點

以下式所示的二維對流擴散方程為例，對比介紹基于這兩種網格的有限容積法所采用的控制方程的差異及方程轉換的要點。

式中：ρ為流體密度；x、y為空間坐標分量；t為時間坐標分量；u、v分別為物理平面上x和y方向的速度分量；φ為通用變量；Γφ為擴散系數；Sφ為源項。

3.1 貼體網格有限容積法控制方程及轉換要點

（1）控制方程變換。貼體網格對應的計算域所基于的坐標系不同于直角坐標系，因此須對式（1）所示的適用于直角坐標系的控制方程進行坐標變換后方能使用。初學者對于貼體坐標系下的控制方程的形式較為陌生，建議以圖2 所示的轉換過程進行教學。通過該圖，讓學生掌握坐標變換的本質，即通過一定的坐標變換法則將物理域上的控制方程轉換為計算域上的控制方程，從而在規則的計算域上進行離散求解。而該轉換法則［1］就是方程轉換的關鍵。接著舉例介紹轉換法則的使用方法：圖2 中ψ為待求導的物理量或物理量的組合，例如式（1）中的ρuφ和Γφ（?φ/?x）等均可作為ψ代入圖中所示的相應坐標變換法則中進行轉換。二階偏導數的轉換等價于作2 次一階偏導數的轉換，例如將Γφ（?φ/?x）作為ψ值進行第1 次坐標變換后得到了下式所示的表達式，其中含有φx（即φ對x的一階偏導），再將其進行一次坐標變換即可［4-6］。此外，通過此圖也讓學生了解為何只有跟空間坐標有關的項需要進行坐標變換，非穩態項無須進行變換，源項中只需變換跟坐標有關的部分。

圖2 貼體網格控制方程的坐標變換方法

式中：J為雅克比系數；（）x，（）ξ，（）η分別表示括號里的變量（或變量組合）對坐標x、ξ和η的偏導數。

（2）邊界條件的變換。通過邊界條件的統一形式aφ+b（?φ/?n）=c來引出邊界條件也應進行坐標變換，原因是其中的一階導數項跟坐標有關。然后以圖3 所示的左邊界為例介紹一階導數的變換方法。

圖3 邊界條件的處理

由圖3 可知，對于物理域左邊界上任意一點P，其切線方向即為η坐標在當地的方向，按照基矢量變換原理［1］可得：gη=xηex+yηey。向學生解釋清楚這一點后，按圖4 所示的流程介紹如何由該切向基矢量推導得到左邊界的外法向導數。接著讓學生舉一反三，推導得到右邊界、上邊界和下邊界的外法向導數的表達式。由此，學生就比較容易掌握不同邊界的邊界條件在貼體坐標系下的表達式的來龍去脈。

圖4 左邊界外法向導數的推導思路

3.2 非結構化網格有限容積法采用的控制方程及邊界條件

與貼體網格不同，非結構化網格基于直角坐標系，因此，進行有限容積法離散時控制方程無須進行坐標變換。但是，應向學生講清楚非結構化網格有限容積法采用矢量形式的控制方程，其原因是：非結構化網格不存在統一的網格節點拓撲結構，難以直接給出統一的偏導數的離散格式，因此采用如下式所示的矢量形式的控制方程［1-3］，即

式中，u為速度矢量。

接著介紹邊界條件處理的特點：邊界條件通用表達式aφ+b（?φ/?n）=c中的一階偏導數通過與邊界相鄰的單元p0的梯度計算得到，即（?φ/?n）=?φp0·n（其中n為由p0指向相鄰節點的方向向量）。

4 兩種網格控制方程離散過程的差異

通過圖5 介紹兩種方法的計算網格系統及網格參數，讓學生了解貼體網格系統跟直角坐標系結構化網格系統無本質區別，而非結構化網格系統與前面所學的結構化網格系統有本質的區別。以此引出貼體網格離散方法與結構化網格完全一致，而非結構化網格的離散是一種全新的方法，本次課將重點學習。

圖5 貼體網格和非結構化網格以及計算節點布置

通過下式所示的對流項的ξ方向分量在貼體網格的離散讓學生理解貼體網格有限容積法方程離散方法跟結構化網格有限容積法完全一樣，即控制容積上的分部積分；不同之處在于貼體網格有限容積法的方程離散是在與物理域不同的計算域上進行［4-6］。其公式為

通過以下兩個步驟介紹非結構化網格有限容積法控制方程的離散方法：①結合高斯散度定理［7-8］，將式（3）在控制容積V上進行積分可得式（5）所示的積分形式的控制方程；②以擴散項為例，按式（6）進行離散。同理可寫出其他項的離散形式［9-10］。

式中：A和A分別為控制容積界面的面積和面積矢量；j為界面的編號。

讓學生觀察式（6）所示的離散格式以及前面所講的邊界條件的表達式，讓他們發現梯度計算在非結構化網格有限容積法離散中是核心。接著介紹求解梯度常采用的Green-Gauss 梯度法［11-13］及最小二乘法［14-15］，具體方法用圖6 進行對照式講解。

圖6 兩種梯度求解法的思路對比

5 兩種網格有限容積法的適應性

采用表2 對比介紹兩種網格有限容積法的差異，讓同學們從整體上把握二者的特征及適應性。

表2 貼體網格與非結構化網格有限容積法綜合對比

結合表2 的分析，可以給出以下建議：對一些幾何形狀相對規則、流動結構較為復雜又需要高精度的問題，建議用貼體網格法，而對于不規則性較大的一般性流動與傳熱問題，建議使用非結構化網格。

6 結語

不規則區域問題是流動與傳熱數值計算普遍遇到的問題，基于貼體網格和非結構化網格有限容積法是解決此類問題普遍采用的方法。作為商業軟件的使用者，應當了解該數值計算理論，作為肩負歷史重任的未來潛在軟件開發者，更是要踏實掌握該根基性的理論。雖然該內容一直是數值傳熱學教學的難點，但只要從本文所述的5 個方面進行對比式教學，可以讓學生在對比學習中深入了解兩種方法的本質，準確掌握兩種方法的技術細節，為將來成為流動傳熱仿真軟件開發人才、實現計算模擬軟件國產化、突破技術封鎖奠定基礎。

·名人名言·

大學的榮譽，不在它的校舍與人數；而在于它一代一代人的質量。

——柯南特