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

（1.北京石油化工學(xué)院機(jī)械工程學(xué)院，北京 102617；2.上海海事大學(xué)商船學(xué)院，上海 201306；3.中國(guó)石油大學(xué)（北京）機(jī)械與儲(chǔ)運(yùn)工程學(xué)院，北京 102249）

0 引言

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

1 貼體網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格的差異

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

圖1 貼體網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格

2 貼體網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格的有限容積法實(shí)施過程差異

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

表1 基于兩種網(wǎng)格系統(tǒng)的有限容積法整體差異對(duì)比

3 貼體網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格有限容積法控制方程的異同及變換要點(diǎn)

以下式所示的二維對(duì)流擴(kuò)散方程為例，對(duì)比介紹基于這兩種網(wǎng)格的有限容積法所采用的控制方程的差異及方程轉(zhuǎn)換的要點(diǎn)。

式中：ρ為流體密度；x、y為空間坐標(biāo)分量；t為時(shí)間坐標(biāo)分量；u、v分別為物理平面上x和y方向的速度分量；φ為通用變量；Γφ為擴(kuò)散系數(shù)；Sφ為源項(xiàng)。

3.1 貼體網(wǎng)格有限容積法控制方程及轉(zhuǎn)換要點(diǎn)

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

圖2 貼體網(wǎng)格控制方程的坐標(biāo)變換方法

式中：J為雅克比系數(shù)；（）x，（）ξ，（）η分別表示括號(hào)里的變量（或變量組合）對(duì)坐標(biāo)x、ξ和η的偏導(dǎo)數(shù)。

（2）邊界條件的變換。通過邊界條件的統(tǒng)一形式aφ+b（?φ/?n）=c來引出邊界條件也應(yīng)進(jìn)行坐標(biāo)變換，原因是其中的一階導(dǎo)數(shù)項(xiàng)跟坐標(biāo)有關(guān)。然后以圖3 所示的左邊界為例介紹一階導(dǎo)數(shù)的變換方法。

圖3 邊界條件的處理

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

圖4 左邊界外法向?qū)?shù)的推導(dǎo)思路

3.2 非結(jié)構(gòu)化網(wǎng)格有限容積法采用的控制方程及邊界條件

與貼體網(wǎng)格不同，非結(jié)構(gòu)化網(wǎng)格基于直角坐標(biāo)系，因此，進(jìn)行有限容積法離散時(shí)控制方程無須進(jìn)行坐標(biāo)變換。但是，應(yīng)向?qū)W生講清楚非結(jié)構(gòu)化網(wǎng)格有限容積法采用矢量形式的控制方程，其原因是：非結(jié)構(gòu)化網(wǎng)格不存在統(tǒng)一的網(wǎng)格節(jié)點(diǎn)拓?fù)浣Y(jié)構(gòu)，難以直接給出統(tǒng)一的偏導(dǎo)數(shù)的離散格式，因此采用如下式所示的矢量形式的控制方程［1-3］，即

式中，u為速度矢量。

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

4 兩種網(wǎng)格控制方程離散過程的差異

通過圖5 介紹兩種方法的計(jì)算網(wǎng)格系統(tǒng)及網(wǎng)格參數(shù)，讓學(xué)生了解貼體網(wǎng)格系統(tǒng)跟直角坐標(biāo)系結(jié)構(gòu)化網(wǎng)格系統(tǒng)無本質(zhì)區(qū)別，而非結(jié)構(gòu)化網(wǎng)格系統(tǒng)與前面所學(xué)的結(jié)構(gòu)化網(wǎng)格系統(tǒng)有本質(zhì)的區(qū)別。以此引出貼體網(wǎng)格離散方法與結(jié)構(gòu)化網(wǎng)格完全一致，而非結(jié)構(gòu)化網(wǎng)格的離散是一種全新的方法，本次課將重點(diǎn)學(xué)習(xí)。

圖5 貼體網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格以及計(jì)算節(jié)點(diǎn)布置

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

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

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

讓學(xué)生觀察式（6）所示的離散格式以及前面所講的邊界條件的表達(dá)式，讓他們發(fā)現(xiàn)梯度計(jì)算在非結(jié)構(gòu)化網(wǎng)格有限容積法離散中是核心。接著介紹求解梯度常采用的Green-Gauss 梯度法［11-13］及最小二乘法［14-15］，具體方法用圖6 進(jìn)行對(duì)照式講解。

圖6 兩種梯度求解法的思路對(duì)比

5 兩種網(wǎng)格有限容積法的適應(yīng)性

采用表2 對(duì)比介紹兩種網(wǎng)格有限容積法的差異，讓同學(xué)們從整體上把握二者的特征及適應(yīng)性。

表2 貼體網(wǎng)格與非結(jié)構(gòu)化網(wǎng)格有限容積法綜合對(duì)比

結(jié)合表2 的分析，可以給出以下建議：對(duì)一些幾何形狀相對(duì)規(guī)則、流動(dòng)結(jié)構(gòu)較為復(fù)雜又需要高精度的問題，建議用貼體網(wǎng)格法，而對(duì)于不規(guī)則性較大的一般性流動(dòng)與傳熱問題，建議使用非結(jié)構(gòu)化網(wǎng)格。

6 結(jié)語

不規(guī)則區(qū)域問題是流動(dòng)與傳熱數(shù)值計(jì)算普遍遇到的問題，基于貼體網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格有限容積法是解決此類問題普遍采用的方法。作為商業(yè)軟件的使用者，應(yīng)當(dāng)了解該數(shù)值計(jì)算理論，作為肩負(fù)歷史重任的未來潛在軟件開發(fā)者，更是要踏實(shí)掌握該根基性的理論。雖然該內(nèi)容一直是數(shù)值傳熱學(xué)教學(xué)的難點(diǎn)，但只要從本文所述的5 個(gè)方面進(jìn)行對(duì)比式教學(xué)，可以讓學(xué)生在對(duì)比學(xué)習(xí)中深入了解兩種方法的本質(zhì)，準(zhǔn)確掌握兩種方法的技術(shù)細(xì)節(jié)，為將來成為流動(dòng)傳熱仿真軟件開發(fā)人才、實(shí)現(xiàn)計(jì)算模擬軟件國(guó)產(chǎn)化、突破技術(shù)封鎖奠定基礎(chǔ)。

·名人名言·

大學(xué)的榮譽(yù)，不在它的校舍與人數(shù)；而在于它一代一代人的質(zhì)量。

——柯南特