劉君，魏雁昕，陳潔

大連理工大學(xué) 航空航天學(xué)院，大連 116024

從2015年開始本文作者團隊致力于激波裝配法研究[1]，采用非結(jié)構(gòu)動網(wǎng)格技術(shù)較好地解決了這類算法遇到的復(fù)雜幾何拓?fù)潆y題，實現(xiàn)了間斷結(jié)構(gòu)的自適應(yīng)裝配[2-3]。計算中不斷進行流場結(jié)構(gòu)辨識和裝配，從幾何角度看，收斂后的流場空間被分割為多個邊界幾何參數(shù)確定的相鄰區(qū)域，原則上可以采用結(jié)構(gòu)網(wǎng)格的多塊拼接技術(shù)來提高效率和精度，近年來探索將激波裝配法和有限差分法高精度格式相結(jié)合，實現(xiàn)全場一致高精度計算。

與常規(guī)的多塊拼接網(wǎng)格技術(shù)不同，由于裝配法的激波邊界在計算過程中存在運動，應(yīng)用前要解決網(wǎng)格變形引起的幾何守恒律問題。通過調(diào)研幾何守恒律文獻，認(rèn)識到貼體坐標(biāo)變換過程中采用的數(shù)學(xué)恒等式在離散條件下不再成立，因此曲線坐標(biāo)系下的離散等價方程應(yīng)該是帶有源項的。所謂的幾何守恒律就是構(gòu)建與對流項離散格式相匹配的坐標(biāo)變換系數(shù)的差商算法，以保證離散條件下源項為0。由于幾何守恒律與差分格式相關(guān)，目前還沒有普適性的算法，本文作者團隊從帶有源項的離散等價方程出發(fā)，采用在不滿足幾何守恒律條件下直接離散源項的新思路，提出差分格式應(yīng)用于曲線坐標(biāo)系的離散準(zhǔn)則[4]，理論上適合任意格式和差商的組合，目前采用一階迎風(fēng)格式[5]和多種WENO(Weighted Essentially Non-Oscillatory)格式進行驗證[6]，計算結(jié)果表明可以消除坐標(biāo)變換引起的數(shù)值誤差。

如果流場中存在激波相交等現(xiàn)象，間斷裝配后經(jīng)常在局部出現(xiàn)復(fù)雜的網(wǎng)格形狀。按照離散準(zhǔn)則應(yīng)用有限差分格式，在離散點僅用到當(dāng)?shù)刈鴺?biāo)變換系數(shù)，利用這個特點建立了基于非結(jié)構(gòu)網(wǎng)格的有限差分法，對于定常斜激波、激波相交和運動激波等空間二維流場，在僅有3條網(wǎng)格線連接離散點上，采用一階迎風(fēng)格式計算獲得穩(wěn)定收斂解，計算精度接近直角坐標(biāo)系均勻結(jié)構(gòu)網(wǎng)格的一階迎風(fēng)格式計算結(jié)果[7]。這種非結(jié)構(gòu)網(wǎng)格有限差分法可以有效解決多塊拼接網(wǎng)格技術(shù)應(yīng)用于激波裝配法時遇到的難題[7]。

本文目的是將非結(jié)構(gòu)網(wǎng)格有限差分法推廣到空間三維并在此基礎(chǔ)上拓展其應(yīng)用范圍。

1 離散等價方程及其離散準(zhǔn)則

考慮笛卡爾坐標(biāo)系(t,x,y,z)下三維守恒型Euler方程：

(1)

從(t,x,y,z)變換到計算空間(τ,ξ,η,ζ)上，得到包含源項的形式：

(2)

式中：

(3)

(4)

在坐標(biāo)變換函數(shù)空間導(dǎo)數(shù)連續(xù)條件下S=0，因此CFD經(jīng)典教科書給出方程形式：

(5)

對于S≠0引起數(shù)值誤差的現(xiàn)象很早就被注意到，1977年Steger發(fā)現(xiàn)在貼體坐標(biāo)系(Body-Fitted Coordinates，BFC)中采用二階中心差分格式計算均勻流場會出現(xiàn)物理量的不守恒，數(shù)值結(jié)果表明非零源項S引起的誤差為二階精度小量[8-9]，該源項由坐標(biāo)變換產(chǎn)生，在流場中造成非物理解，嚴(yán)重時會引起數(shù)值振蕩，影響計算穩(wěn)定性[10-11]。Thomas和Lombard在研究動網(wǎng)格時提出了目前這個研究領(lǐng)域常用的幾何守恒律概念，對于有限差分法，如果在離散條件下式(3)等于0稱為滿足體積守恒律，如果式(4)等于0稱為滿足面積守恒律，進一步論證滿足幾何守恒律不會出現(xiàn)坐標(biāo)變換的數(shù)值誤差[12-13]。20世紀(jì)CFD領(lǐng)域以二階精度格式為主，基于Steger得出的結(jié)論，認(rèn)為S≠0引起的誤差不嚴(yán)重，相關(guān)研究工作逐步停止，直到20多年后高精度格式進入到應(yīng)用階段，幾何守恒律的研究才再次成為熱點[14]。近10年來，幾何守恒律研究取得很大進展，相關(guān)研究可以參考文獻[15]。

有限體積法中計算面積和體積幾何參數(shù)是代數(shù)表達式，如果注意到二者的相容性，不存在面積守恒律問題，只有應(yīng)用變形網(wǎng)格才需要解決體積守恒律。本文作者團隊采用非結(jié)構(gòu)動網(wǎng)格技術(shù)模擬多體分離問題，為了保證流固耦合界面?zhèn)鬟f的能量守恒，提出了一種新的體積守恒律，其主要特點是將不守恒量作為源項加入到對流項中[16-17]。解決有限差分法的動網(wǎng)格問題時繼續(xù)沿著這個思路，采用直接離散S≠0的源項算法。根據(jù)源項來自于坐標(biāo)變換過程，2018年提出離散等價方程源項和對流項的相容性離散準(zhǔn)則(以下簡稱離散準(zhǔn)則)[4]。按照準(zhǔn)則應(yīng)用差分格式得到的源項和對流項的計算形式有相似之處，合并處理可以提高計算效率，先是通過算例測試表明可行性[4]，然后進行理論證明合并處理不會引入過程誤差[5]。文獻[4-6]中離散準(zhǔn)則和驗證算例是基于空間二維問題的，下面從空間三維方程出發(fā)進行介紹。

源項表示為

(6)

采用算子表示離散格式，例如δξ為對流項導(dǎo)數(shù)?/?ξ對應(yīng)差分格式，Φ和[Γ,Κ,Π]呈線性函數(shù)，可得如下離散形式：

(7)

按照離散準(zhǔn)則，源項也使用同樣的格式離散，得到

(8)

代入離散等價方程，得到半離散形式：

(9)

半離散形式表示為

(10)

采用同樣的差分格式離散齊次方程式(5)，得到的算子形式為

(11)

通過對比可以看出，如果齊次方程離散形式式(11)采用“凍結(jié)坐標(biāo)變換系數(shù)”就得到式(10)。從經(jīng)典CFD教科書角度看，前者似乎降低了精度，但是從推導(dǎo)過程看，按照離散準(zhǔn)則應(yīng)用差分格式，計算過程包含源項、自動消除坐標(biāo)變換引起的誤差，式(10)就是離散算子δ的精度。根據(jù)Steger的結(jié)論[8-9]，源項S≠0是二階精度小量，即在離散條件下齊次方程(5)本身就是原始出發(fā)方程(1)的二階精度的近似，不管算子δ精度如何，均無法消除坐標(biāo)變換源項S≠0引入的誤差。

采用RHS表示空間離散格式計算結(jié)果，得到

(12)

可以看出，對于均勻流場RHS=0，不論網(wǎng)格如何變形，不會出現(xiàn)誤差，現(xiàn)有二維算例表明離散準(zhǔn)則可以同時消除不滿足體積守恒律和面積守恒律引起的數(shù)值誤差[4-6]。

基于上述分析，對幾何守恒律研究進展進行調(diào)研分析[13]，提出如下看法：鄧小剛等提出的SCMM(Symmetrical Conservative Metric Method)算法很好地解決了緊致類格式的幾何守恒律問題[18-19]，近年來部分學(xué)者借鑒SCMM思路構(gòu)建WENO 格式的幾何守恒律，做法是把格式分解為中心差分和數(shù)值耗散兩部分[20-22]，其中具有迎風(fēng)特性的數(shù)值耗散部分在滿足幾何守恒律條件下作為源項直接加入，從計算結(jié)果看，可以有效降低誤差，但是不能完全消除，給出的均勻流誤差曲線存在緩慢上升。

幾何守恒律和離散準(zhǔn)則的目標(biāo)是消除有限差分法坐標(biāo)變換引起的數(shù)值誤差，出發(fā)點均為帶有S≠0源項的離散等價方程，但是技術(shù)路線完全不同，最終達到的效果也存在如下明顯差異：

1)坐標(biāo)變換系數(shù)采用準(zhǔn)確的解析表達式是不滿足幾何守恒律的。離散準(zhǔn)則對坐標(biāo)變換系數(shù)沒有限制，解析式或任何差商格式均可。

2)不同差分格式需要不同的幾何守恒律算法。所謂幾何守恒律就是根據(jù)網(wǎng)格幾何參數(shù)計算坐標(biāo)變換系數(shù)的差商，常見的面積和體積計算表達式大多是相對于形心的中心差商，很少考慮迎風(fēng)特性，原則上幾何守恒律僅適用于中心格式，至今未看到MUSCL(Monotone Upstream-centred Schemes for Conservation Laws)和帶有minmod限制器等包含迎風(fēng)格式的幾何守恒律。從以上推導(dǎo)看，離散準(zhǔn)則對差分格式?jīng)]有限制，給出的驗證算例不僅包含一階迎風(fēng)格式[4-5]，對于WENO格式，誤差曲線同樣保持平直[6]。

文獻調(diào)研未發(fā)現(xiàn)幾何守恒律在六面體結(jié)構(gòu)網(wǎng)格以外的應(yīng)用情況。下面將介紹利用離散準(zhǔn)則的特點，構(gòu)建基于非結(jié)構(gòu)網(wǎng)格的有限差分法。

2 非結(jié)構(gòu)網(wǎng)格的有限差分法

對比算子表達式式(9)和式(10)，按照離散準(zhǔn)則應(yīng)用差分格式和經(jīng)典CFD的主要差異在于離散點以外的相鄰網(wǎng)格上通量分裂的計算表達式。使用離散準(zhǔn)則時明顯增加計算量，對于三維問題，一階迎風(fēng)格式涉及到7個節(jié)點，離散準(zhǔn)則需要進行6次通量分裂格式的計算，應(yīng)用于TVD(Total Variation Diminishing)等二階格式涉及到13個節(jié)點，需要12次通量分裂；常規(guī)的有限差分法每個點上通量分裂格式只要1次，每個方向掃描1次，不管一階格式還是高階格式，通量分裂只需要計算3次即可。

為了消除坐標(biāo)變換引起的數(shù)值誤差，離散準(zhǔn)則確實付出了較高的計算效率為代價，不過利用差分格式中僅需要當(dāng)?shù)刈鴺?biāo)變換系數(shù)的特性，可以在非結(jié)構(gòu)網(wǎng)格點建立局部貼體坐標(biāo)系從而采用有限差分格式進行計算。

為便于原理介紹，先針對二維情況進行討論。差分格式沿著網(wǎng)格線離散，從物理空間(x,y)映射到計算空間(ξ,η)至少需要離散點與周圍4個鄰點連成4條網(wǎng)格線，對于離散點0和周圍3個鄰點只能連成3條網(wǎng)格線，構(gòu)成的非結(jié)構(gòu)網(wǎng)格如圖1所示。

圖1 二維坐標(biāo)變換對應(yīng)關(guān)系

采用差商近似這些坐標(biāo)變換系數(shù)，影響數(shù)值解的是網(wǎng)格坐標(biāo)的相對位置。例如，采用二階中心差商求解，1-0-2三點連成曲線的斜率和1-2兩點直線段斜率一樣，因此把1-2線段平移到d-b不影響斜率。同樣處理1-0-3三點連成的曲線，得到對應(yīng)的線段a-c。圍繞離散點0的a、b、c、d這4個點構(gòu)成的網(wǎng)格可以得到坐標(biāo)變換系數(shù)數(shù)值解，結(jié)果與1、2、3點構(gòu)成的網(wǎng)格是一樣的。坐標(biāo)變換系數(shù)的表達式為

(13)

二維算例表明流場沒有明顯不對稱現(xiàn)象，且網(wǎng)格的相鄰線段間最大夾角由120°變化到180°過程中，計算結(jié)果無明顯變化[7]。

三維空間中存在3條坐標(biāo)軸，結(jié)構(gòu)網(wǎng)格除了離散點外還需要6個相鄰點。由于三維幾何特性比較復(fù)雜，當(dāng)4個相鄰點之間滿足如下條件：離散點0落在其他4個點構(gòu)成的四面體內(nèi)部，那么可以采用以上二維空間類似的方法進行映射，即在相鄰4點基礎(chǔ)上構(gòu)造非結(jié)構(gòu)網(wǎng)格并且進行有限差分法計算。圖2為0-1線段進行3次映射的示意圖，把曲線段1-0-4映射成計算空間ξ軸、0-1為正向，曲線1-0-3映射成計算空間η軸、0-3為正向，曲線1-0-2映射成計算空間ζ軸、0-1為正向。以變換矩陣行列式為例來說明坐標(biāo)變換系數(shù)的計算方法：

圖2 三維非結(jié)構(gòu)網(wǎng)格示意圖

(14)

在各點坐標(biāo)ri=(xi,yi,zi)已知的情況下，采用中心差商計算式(14)中偏導(dǎo)數(shù)，例如：

xξ=0.5(x1-x4)

xη=0.5(x3-x1)

xζ=0.5(x1-x2)

這是對0-1線段多重映射形成的一種局部坐標(biāo)系，考慮到映射在計算空間坐標(biāo)軸方向特性，還應(yīng)該有5種組合方式，總計需要6次計算。在非結(jié)構(gòu)網(wǎng)格中4條連線等權(quán)重，組合以后存在24種局部坐標(biāo)系，大大增加計算量，定性分析沒有實用價值。為了改善這種非結(jié)構(gòu)網(wǎng)格有限差分法的計算效率，提出在直角坐標(biāo)系下結(jié)構(gòu)網(wǎng)格基礎(chǔ)上局部應(yīng)用非結(jié)構(gòu)有限差分法的扎染算法。

3 扎染算法

以二維半模球柱組合體為例，對扎染算法進行原理介紹。在包含物體在內(nèi)的整個計算區(qū)域上分布直角坐標(biāo)系下的均勻結(jié)構(gòu)化網(wǎng)格，如圖3(a)所示，網(wǎng)格生成只需要一個控制參數(shù)Δx。根據(jù)網(wǎng)格點位置和物體幾何參數(shù)之間的相對關(guān)系，將已經(jīng)生成的空間網(wǎng)格點分為3類：用黃色實心圓標(biāo)識的是完全位于物體內(nèi)部或物體表面附近點，定義為內(nèi)部節(jié)點；用紅色矩形標(biāo)識的是受物體影響的外部點，定義為非結(jié)構(gòu)節(jié)點；沒有標(biāo)識的是不受物體影響的外部點，定義為外部節(jié)點。標(biāo)識內(nèi)/外部節(jié)點需要設(shè)定判斷準(zhǔn)則，本文采用與物體間距d<0.5Δx作為準(zhǔn)則，結(jié)果如圖3(b)所示。

圖3 扎染算法原理示意圖

由于非結(jié)構(gòu)節(jié)點的4個相鄰網(wǎng)格點中存在物體內(nèi)部節(jié)點，該點不能直接采用直角均勻網(wǎng)格計算。在非結(jié)構(gòu)點集和物面之間的狹小空間范圍內(nèi)使用非結(jié)構(gòu)有限差分法進行求解。上文介紹的非結(jié)構(gòu)網(wǎng)格有限差分法，本質(zhì)是通過多重映射來滿足曲線局部坐標(biāo)系下建立結(jié)構(gòu)網(wǎng)格的坐標(biāo)軸需求，按照這個思路本文提出一種不需要進行物面剖分網(wǎng)格的新算法：① 如果非結(jié)構(gòu)點只連接1個內(nèi)部節(jié)點，在兩點之間按照最小距離準(zhǔn)則尋找與其對應(yīng)的的物面點，由4條線可以構(gòu)成結(jié)構(gòu)網(wǎng)格，如圖4(a)所示；② 如果非結(jié)構(gòu)點連接2個不同坐標(biāo)軸上的內(nèi)部節(jié)點，那么物面必然位于該非結(jié)構(gòu)點作為角點的4個直角網(wǎng)格內(nèi)，同樣按照垂直投影或者最小距離準(zhǔn)則尋找物面點，由3條線可以構(gòu)成非結(jié)構(gòu)網(wǎng)格，如圖4(b)所示；③ 如果非結(jié)構(gòu)點連接2個同坐標(biāo)軸上的內(nèi)部節(jié)點，或者連接3個內(nèi)部點，表明物體局部幾何構(gòu)型非常復(fù)雜，需要進行特殊處理，圖3所示的球柱組合體外形沒有這種情況。最后建立的非結(jié)構(gòu)點集和物面之間的網(wǎng)格，如圖3(c)所示，藍色節(jié)點為按照最小距離準(zhǔn)則尋找的物面節(jié)點。

圖4 非結(jié)構(gòu)節(jié)點連接關(guān)系

首先，在圖3(a)的均勻網(wǎng)格上進行計算，因為不需要Δx以外的其他網(wǎng)格相關(guān)信息，無需進行坐標(biāo)變換，通量分裂計算表達式非常簡單，對于確定的差分格式，這種結(jié)構(gòu)網(wǎng)格可以實現(xiàn)最高的計算效率。目前，本文算法只能采用時間顯式格式。定義全場標(biāo)量cij，在圖3(b)示中非結(jié)構(gòu)網(wǎng)格點和內(nèi)部網(wǎng)格點上cij=0，其余點cij=1。由于cij=0點不影響物體外部流場，原則上可以不用計算。為避免跳過這些點需要的邏輯判斷，時間推進采用如下統(tǒng)一計算式：

(15)

下文中對是否進行跳點邏輯判斷進行了效率對比分析。

其次，僅在圖3(c)所示的混合網(wǎng)格區(qū)域內(nèi)(紅色和藍色節(jié)點)進行計算，求解式(15)中cij=0的非結(jié)構(gòu)節(jié)點和物面邊界節(jié)點參數(shù)，物體邊界對流場的影響在這一步體現(xiàn)。從圖3(c)可以看出，即使是非結(jié)構(gòu)網(wǎng)格，離散點3條連線中只有1條特殊線，根據(jù)其所處象限確定唯一的映射關(guān)系或局部坐標(biāo)系，只需計算1次，不像第2節(jié)中進行12次計算后算術(shù)平均。

混合網(wǎng)格區(qū)域內(nèi)非結(jié)構(gòu)節(jié)點與物面節(jié)點一一對應(yīng)，很多工程應(yīng)用中，物面節(jié)點數(shù)在總網(wǎng)格中所占的比例較小，所以大部分流場計算是在均勻正交網(wǎng)格下完成，計算效率高。物體內(nèi)部對流場無貢獻，始終保持初始流場參數(shù)，如果云圖顯示整個流場，總有部分區(qū)域同色，如圖3(d)所示。本文的方法是一種直角坐標(biāo)系均勻網(wǎng)格和局部曲線坐標(biāo)系非結(jié)構(gòu)網(wǎng)格的拼接算法，計算過程、流場圖像和中國傳統(tǒng)的扎染工藝有相似之處：根據(jù)設(shè)計圖形把部分白布包扎起來，直接投入染缸，操作簡單，從染缸出來后，為了精細(xì)化表現(xiàn)主題，保持變成藍底色部分不變，在其他區(qū)域進行局部人工彩繪，很多成品中還保留部分沒有上色的白布，所以稱其為扎染算法。

扎染算法應(yīng)用于三維空間，點的分類過程完全相同：每個方向上根據(jù)距離物體d>0.5Δx準(zhǔn)則進行分類，標(biāo)識出內(nèi)部節(jié)點，然后把相鄰點包含有內(nèi)部節(jié)點的離散點標(biāo)記為非結(jié)構(gòu)點。在每個點的局部坐標(biāo)系下差分離散是按照3個坐標(biāo)軸6個節(jié)點進行的，只是離散點之間的連接關(guān)系按照非結(jié)構(gòu)網(wǎng)格處理，因此，二維算法可以直接推廣至三維空間中。盡管三維問題增加了幾何復(fù)雜性，但是如果在每個坐標(biāo)軸方向不多于1個內(nèi)部節(jié)點，按照上面的算法也可以迅速實現(xiàn)：① 非結(jié)構(gòu)節(jié)點連接1個內(nèi)部點，延伸坐標(biāo)軸到物面找點；② 非結(jié)構(gòu)節(jié)點連接2個內(nèi)部點，對沒有交點的坐標(biāo)軸不需要特殊處理，退化為包含內(nèi)部點的二維問題，根據(jù)最小距離準(zhǔn)則尋找物面點，進行2次映射；③ 非 結(jié)構(gòu)節(jié)點連接3個內(nèi)部點，延伸坐標(biāo)軸到物面3點，構(gòu)成三角形平面，離散點和三角形面心或垂點的連線進行3次映射。對于連接內(nèi)部點大于3個或同一坐標(biāo)軸上有2個內(nèi)部點的情況，根據(jù)工程經(jīng)驗，這些情況下幾何構(gòu)型非常復(fù)雜，要么需要加密網(wǎng)格來刻畫，要么需要修改簡化模型，總體來說需要人工干預(yù)，本文暫時不考慮這些特殊情況。

扎染算法計算流程簡述如下：

Step1讀取待計算物體模型(描述物體的散點集合或函數(shù))，給定計算區(qū)域邊界，在此基礎(chǔ)上對流場中所有節(jié)點屬性進行標(biāo)識，確定結(jié)構(gòu)節(jié)點、非結(jié)構(gòu)節(jié)點、物面節(jié)點和內(nèi)部節(jié)點。

Step2在包含物體的整個區(qū)域采用直角坐標(biāo)系下的均勻網(wǎng)格進行計算，直接使用常規(guī)的有限差分方法進行求解，該過程中無需引入坐標(biāo)變換，計算效率高。

Step3由于內(nèi)部節(jié)點對流場計算不產(chǎn)生貢獻，且Step 2計算中沒有考慮物體的影響，因此將Step 2中內(nèi)部節(jié)點和非結(jié)構(gòu)節(jié)點的計算結(jié)果舍去。

Step4將非結(jié)構(gòu)節(jié)點和物面節(jié)點一一對應(yīng)，建立非結(jié)構(gòu)網(wǎng)格“帶”，在此區(qū)域內(nèi)使用非結(jié)構(gòu)有限差分法進行求解，得到物面節(jié)點和非結(jié)構(gòu)節(jié)點的流場參數(shù)。

Step5將Step 2中計算得到的結(jié)構(gòu)網(wǎng)格節(jié)點與Step 4中得到的非結(jié)構(gòu)網(wǎng)格節(jié)點及物面節(jié)點進行疊加，得到該時刻完整的流場信息。

Step6重復(fù)Step 1～Step 5進行時間推進求解，達到收斂條件，輸出流場結(jié)果。

4 演示算例

選取二維/三維超聲速流動問題驗證本文算法。其中二維算例用于說明該方法的適用性及收斂穩(wěn)定性，三維算例通過對比及理論分析展示扎染算法計算效率高的特點。

4.1 二維圓柱繞流

使用扎染算法對二維圓柱繞流問題進行數(shù)值模擬，驗證該方法在求解超聲速流動時的適用性。算例來流馬赫數(shù)Ma∞=3.0，無量綱流動參數(shù)(ρ,u,v,p)=(1.4,3.0,0,1.0)，流場初始化為來流，時間離散采用Rung-Kutta法，空間離散使用一階迎風(fēng)格式。無量綱計算域為x∈[0,0.1],y∈[0,0.4]的矩形區(qū)域，圓柱圓心位于(0.1,0.2)，半徑R=0.045，分別選取5組不同Δx的網(wǎng)格進行數(shù)值模擬。定義5組中節(jié)點間距值最大的算例為Δx0=0.002 5，其余4組網(wǎng)格控制參數(shù)分別為Δx0/2、Δx0/4、Δx0/8以及Δx0/16。

5種網(wǎng)格條件下的收斂曲線如圖5所示，殘差單調(diào)下降，沒有出現(xiàn)異常波動，表明扎染算法可以得到穩(wěn)定的收斂解，圖中t為無量綱時間。圖6(a)為5種網(wǎng)格條件下計算得到的對稱軸線上壓力P隨x變化曲線，物體表面壓力分布曲線如圖6(b)所示。從圖6(a)看出，隨著網(wǎng)格的不斷加密，對激波的分辨能力逐漸增強，對稱軸線上壓力間斷更加陡峭，峰值也相應(yīng)增大。當(dāng)進一步加密網(wǎng)格為Δx0/8和Δx0/16時，藍色曲線與綠色曲線重合，后續(xù)即使繼續(xù)加密網(wǎng)格，壓力曲線也不再發(fā)生明顯變化。從圖6(b)物面壓力分布曲線同樣可以看出上述結(jié)論，隨著網(wǎng)格加密，駐點處壓力峰值逐漸增大，當(dāng)網(wǎng)格加密至一定程度后，壓力不再隨網(wǎng)格發(fā)生變化。在網(wǎng)格收斂性測試中通常采用沿滯止流線處駐點壓力來評估計算準(zhǔn)確性，其理論值由Rayleign-pitot公式求得

圖5 殘差收斂曲線

圖6 不同網(wǎng)格尺度下壓力曲線

(16)

由于計算域中結(jié)構(gòu)化網(wǎng)格占比大，在該類型網(wǎng)格點處可以直接使用各種成熟的計算格式，圖7 為一階迎風(fēng)格式和5階WENO格式在相同網(wǎng)格條件下的收斂流場。上半部分為一階迎風(fēng)格式計算結(jié)果，下半部分則是WENO格式的計算結(jié)果，通過對比密度和壓力流場可以看出兩種格式結(jié)果基本一致，僅在激波位置處存在一定的差異。圖7的結(jié)果表明扎染算法可以在各種離散格式下使用，均可得到穩(wěn)定的收斂解。

圖7 二維球頭繞流：一階迎風(fēng)格式與WENO格式密度和壓力云圖對比

4.2 三維球頭繞流

4.2.1 跳點法與覆蓋法效率對比

以三維球頭繞流為例，對扎染算法計算效率進行初步分析討論。計算模型如圖8所示，計算區(qū)域無量綱長度x=0.08、y=0.28、z=0.28，球頭半徑R=0.045。流場初始化及來流參數(shù)與4.1節(jié)二維算例相同。

圖8(b)為初始流場xy平面切片，圖中白色半圓內(nèi)為物體區(qū)域，扎染算法求解過程中物體內(nèi)部對流場計算無貢獻。對內(nèi)部節(jié)點分別采用兩種處理方法，方法1表述為每一個時間步內(nèi)對全場結(jié)構(gòu)網(wǎng)格進行循環(huán)求解，后將對流場計算無貢獻的內(nèi)部節(jié)點物理量賦為來流值，稱之為覆蓋法；方法1表述為循環(huán)過程中對結(jié)構(gòu)網(wǎng)格的節(jié)點屬性進行判定，當(dāng)節(jié)點位于物體內(nèi)部時跳過該節(jié)點進行下一行列循環(huán)，內(nèi)部節(jié)點不進行求解，稱之為跳點法。分別采用兩種處理方法求解三維球頭繞流，網(wǎng)格控制參數(shù)Δx=Δy=Δz=0.002,節(jié)點個數(shù)Nx=41,Ny=Nz=141。CFL數(shù)取值0.5，時間離散采用Runge-Kutta法，空間離散采用一階迎風(fēng)格式，通量分裂使用Vanleer求解器。網(wǎng)格節(jié)點屬性如圖9所示，綠色標(biāo)識為結(jié)構(gòu)化節(jié)點，紅色標(biāo)識為非結(jié)構(gòu)節(jié)點，藍色標(biāo)識為物面邊界節(jié)點，物體內(nèi)部節(jié)點未顯示。

圖8 扎染算法計算模型

圖9 xy切面流場節(jié)點屬性

由于內(nèi)部節(jié)點求解與流場無關(guān)，故兩種處理方法計算得到的流場結(jié)果完全一致，如圖10所示。為對比兩種方法的計算效率，分別在同一臺式計算機(主頻3.0 GHz，內(nèi)存16 G)上完成無量綱時間t=0.8時長計算，此時流場完全收斂。圖11 給出了兩種方法CPU耗時隨無量綱時間的變化曲線，從圖中可知覆蓋法CPU計算時間小于跳點法，計算時長約為后者的80%。由于內(nèi)部節(jié)點對非結(jié)構(gòu)點求解無影響，因此兩種方法計算中非結(jié)構(gòu)節(jié)點計算耗時相同，即圖11中紅色曲線所示。從表1中分析可知，計算區(qū)域節(jié)點總數(shù)818 002，內(nèi)部節(jié)點占總數(shù)的3.21%，非結(jié)構(gòu)點和物面邊界點均占總數(shù)的0.35%。覆蓋法計算過程中雖然對內(nèi)部節(jié)點進行求解，但是無需對節(jié)點是否位于物體內(nèi)部進行邏輯判斷，單位時間步計算時間0.301 s。跳點法求解時，雖然求解的節(jié)點個數(shù)減少3.21%，但是在計算中需要對節(jié)點類型判定，在此過程中邏輯運算帶來的時間增量遠大于內(nèi)部節(jié)點的求解時間，單位時間步計算時間0.378 s。綜上所述，扎染算法中由于物體內(nèi)部節(jié)點個數(shù)占節(jié)點總數(shù)比例較小，選取減少邏輯運算的覆蓋法可以有效提高計算效率，本文算例中均采用覆蓋法對內(nèi)部節(jié)點進行處理。

圖10 三維球頭繞流：密度/壓力云圖

圖11 扎染算法CPU耗時對比

表1 節(jié)點屬性分布

4.2.2 常規(guī)差分算法與扎染算法計算效率分析

使用常規(guī)差分算法處理三維曲邊界問題時，在計算域中的所有網(wǎng)格點均需要由物理坐標(biāo)系轉(zhuǎn)到貼體坐標(biāo)系下進行格式離散。以本文采用的VanLeer格式為例，求解通量形式。定義ξ方向馬赫數(shù)：

(17)

式中：

(18)

(19)

式中：

(20)

式中：ρ、u、v、w、p分別表示密度、x方向速度(直角坐標(biāo)系下)、y方向速度(直角坐標(biāo)系下)、z方向速度(直角坐標(biāo)系下)、壓力。

當(dāng)|Maξ|<1(局部亞聲速)時，

(21)

式中：

(22)

(23)

(24)

其中：

(25)

從式(17)～式(25)中可知對于常規(guī)結(jié)構(gòu)網(wǎng)格的有限差分法，通量求解過程中包含坐標(biāo)變換系數(shù)，在增大計算量的同時內(nèi)存占用也相應(yīng)增加，計算效率較低。扎染算法求解過程中在遠離物面的空間計算區(qū)域內(nèi)為均勻且正交的結(jié)構(gòu)化網(wǎng)格，在此類節(jié)點求解時無需使用坐標(biāo)變換系數(shù)，通量求解形式如下，以x方向超聲速流動條件為例，定義方向馬赫數(shù)：

(26)

當(dāng)|Max|≥1時，

(27)

式中：

(28)

由式(26)～式(28)可知，均勻正交網(wǎng)格條件下，通量求解公式簡化，節(jié)省計算時間的同時占用內(nèi)存更小，相較于常規(guī)結(jié)構(gòu)網(wǎng)格有限差分算法計算效率更高，更加適合進行大規(guī)模并行計算。

5 結(jié)論與展望

本文以非結(jié)構(gòu)網(wǎng)格有限差分方法為基礎(chǔ)，提出一種局部多重分塊網(wǎng)格的新算法。在均勻正交網(wǎng)格區(qū)域使用無坐標(biāo)變換的差分算法，計算效率高；物面附近使用非結(jié)構(gòu)有限差分法求解，可以有效處理復(fù)雜的拓?fù)浣Y(jié)構(gòu)。從上述算例演示可知，扎染算法適用于多種計算格式且滿足穩(wěn)定性要求，通過對求解公式的對比分析表明扎染算法數(shù)據(jù)結(jié)構(gòu)簡單，計算效率高，有利于開展大規(guī)模網(wǎng)格并行計算。

通過本文理論證明及算例驗證，定性分析了扎染算法的優(yōu)勢及應(yīng)用前景：

1)計算中大部分區(qū)域為均勻正交結(jié)構(gòu)化網(wǎng)格，通量分裂中無需引入坐標(biāo)變換系數(shù)，理論上可以達到最高的計算效率。

2)扎染算法應(yīng)用中主要計算區(qū)域采用無需坐標(biāo)變換的直角結(jié)構(gòu)網(wǎng)格，使用差分方法進行求解，對網(wǎng)格分區(qū)無要求，適應(yīng)于任意構(gòu)架的高性能計算平臺。

3)由于物面節(jié)點是由空間非結(jié)構(gòu)節(jié)點依據(jù)最小距離判據(jù)確定，因此僅需可以描繪外形的散點集合即可開展數(shù)值模擬，無需人工進行(表面/空間)網(wǎng)格生成，實現(xiàn)高度自動化的快速流場計算。

當(dāng)前研究中針對具有解析表達式的簡單三維外形已經(jīng)實現(xiàn)扎染算法計算，后續(xù)計劃與相關(guān)學(xué)者開展合作研究，將CAD等通用建模軟件與本文所提算法相結(jié)合，最終實現(xiàn)輸入實體輸出流場結(jié)果的快速計算，節(jié)省人力成本，提高計算效率。