劉君,魏雁昕,陳潔
大連理工大學 航空航天學院,大連 116024
從2015年開始本文作者團隊致力于激波裝配法研究[1],采用非結構動網格技術較好地解決了這類算法遇到的復雜幾何拓撲難題,實現了間斷結構的自適應裝配[2-3]。計算中不斷進行流場結構辨識和裝配,從幾何角度看,收斂后的流場空間被分割為多個邊界幾何參數確定的相鄰區域,原則上可以采用結構網格的多塊拼接技術來提高效率和精度,近年來探索將激波裝配法和有限差分法高精度格式相結合,實現全場一致高精度計算。
與常規的多塊拼接網格技術不同,由于裝配法的激波邊界在計算過程中存在運動,應用前要解決網格變形引起的幾何守恒律問題。通過調研幾何守恒律文獻,認識到貼體坐標變換過程中采用的數學恒等式在離散條件下不再成立,因此曲線坐標系下的離散等價方程應該是帶有源項的。所謂的幾何守恒律就是構建與對流項離散格式相匹配的坐標變換系數的差商算法,以保證離散條件下源項為0。由于幾何守恒律與差分格式相關,目前還沒有普適性的算法,本文作者團隊從帶有源項的離散等價方程出發,采用在不滿足幾何守恒律條件下直接離散源項的新思路,提出差分格式應用于曲線坐標系的離散準則[4],理論上適合任意格式和差商的組合,目前采用一階迎風格式[5]和多種WENO(Weighted Essentially Non-Oscillatory)格式進行驗證[6],計算結果表明可以消除坐標變換引起的數值誤差。
如果流場中存在激波相交等現象,間斷裝配后經常在局部出現復雜的網格形狀。按照離散準則應用有限差分格式,在離散點僅用到當地坐標變換系數,利用這個特點建立了基于非結構網格的有限差分法,對于定常斜激波、激波相交和運動激波等空間二維流場,在僅有3條網格線連接離散點上,采用一階迎風格式計算獲得穩定收斂解,計算精度接近直角坐標系均勻結構網格的一階迎風格式計算結果[7]。這種非結構網格有限差分法可以有效解決多塊拼接網格技術應用于激波裝配法時遇到的難題[7]。
本文目的是將非結構網格有限差分法推廣到空間三維并在此基礎上拓展其應用范圍。
考慮笛卡爾坐標系(t,x,y,z)下三維守恒型Euler方程:
(1)
從(t,x,y,z)變換到計算空間(τ,ξ,η,ζ)上,得到包含源項的形式:
(2)
式中:
(3)
(4)
在坐標變換函數空間導數連續條件下S=0,因此CFD經典教科書給出方程形式:
(5)

對于S≠0引起數值誤差的現象很早就被注意到,1977年Steger發現在貼體坐標系(Body-Fitted Coordinates,BFC)中采用二階中心差分格式計算均勻流場會出現物理量的不守恒,數值結果表明非零源項S引起的誤差為二階精度小量[8-9],該源項由坐標變換產生,在流場中造成非物理解,嚴重時會引起數值振蕩,影響計算穩定性[10-11]。Thomas和Lombard在研究動網格時提出了目前這個研究領域常用的幾何守恒律概念,對于有限差分法,如果在離散條件下式(3)等于0稱為滿足體積守恒律,如果式(4)等于0稱為滿足面積守恒律,進一步論證滿足幾何守恒律不會出現坐標變換的數值誤差[12-13]。20世紀CFD領域以二階精度格式為主,基于Steger得出的結論,認為S≠0引起的誤差不嚴重,相關研究工作逐步停止,直到20多年后高精度格式進入到應用階段,幾何守恒律的研究才再次成為熱點[14]。近10年來,幾何守恒律研究取得很大進展,相關研究可以參考文獻[15]。
有限體積法中計算面積和體積幾何參數是代數表達式,如果注意到二者的相容性,不存在面積守恒律問題,只有應用變形網格才需要解決體積守恒律。本文作者團隊采用非結構動網格技術模擬多體分離問題,為了保證流固耦合界面傳遞的能量守恒,提出了一種新的體積守恒律,其主要特點是將不守恒量作為源項加入到對流項中[16-17]。解決有限差分法的動網格問題時繼續沿著這個思路,采用直接離散S≠0的源項算法。根據源項來自于坐標變換過程,2018年提出離散等價方程源項和對流項的相容性離散準則(以下簡稱離散準則)[4]。按照準則應用差分格式得到的源項和對流項的計算形式有相似之處,合并處理可以提高計算效率,先是通過算例測試表明可行性[4],然后進行理論證明合并處理不會引入過程誤差[5]。文獻[4-6]中離散準則和驗證算例是基于空間二維問題的,下面從空間三維方程出發進行介紹。

源項表示為
(6)
采用算子表示離散格式,例如δξ為對流項導數?/?ξ對應差分格式,Φ和[Γ,Κ,Π]呈線性函數,可得如下離散形式:

(7)
按照離散準則,源項也使用同樣的格式離散,得到

(8)
代入離散等價方程,得到半離散形式:
(9)


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

采用RHS表示空間離散格式計算結果,得到
(12)
可以看出,對于均勻流場RHS=0,不論網格如何變形,不會出現誤差,現有二維算例表明離散準則可以同時消除不滿足體積守恒律和面積守恒律引起的數值誤差[4-6]。
基于上述分析,對幾何守恒律研究進展進行調研分析[13],提出如下看法:鄧小剛等提出的SCMM(Symmetrical Conservative Metric Method)算法很好地解決了緊致類格式的幾何守恒律問題[18-19],近年來部分學者借鑒SCMM思路構建WENO 格式的幾何守恒律,做法是把格式分解為中心差分和數值耗散兩部分[20-22],其中具有迎風特性的數值耗散部分在滿足幾何守恒律條件下作為源項直接加入,從計算結果看,可以有效降低誤差,但是不能完全消除,給出的均勻流誤差曲線存在緩慢上升。
幾何守恒律和離散準則的目標是消除有限差分法坐標變換引起的數值誤差,出發點均為帶有S≠0源項的離散等價方程,但是技術路線完全不同,最終達到的效果也存在如下明顯差異:
1)坐標變換系數采用準確的解析表達式是不滿足幾何守恒律的。離散準則對坐標變換系數沒有限制,解析式或任何差商格式均可。
2)不同差分格式需要不同的幾何守恒律算法。所謂幾何守恒律就是根據網格幾何參數計算坐標變換系數的差商,常見的面積和體積計算表達式大多是相對于形心的中心差商,很少考慮迎風特性,原則上幾何守恒律僅適用于中心格式,至今未看到MUSCL(Monotone Upstream-centred Schemes for Conservation Laws)和帶有minmod限制器等包含迎風格式的幾何守恒律。從以上推導看,離散準則對差分格式沒有限制,給出的驗證算例不僅包含一階迎風格式[4-5],對于WENO格式,誤差曲線同樣保持平直[6]。

文獻調研未發現幾何守恒律在六面體結構網格以外的應用情況。下面將介紹利用離散準則的特點,構建基于非結構網格的有限差分法。
對比算子表達式式(9)和式(10),按照離散準則應用差分格式和經典CFD的主要差異在于離散點以外的相鄰網格上通量分裂的計算表達式。使用離散準則時明顯增加計算量,對于三維問題,一階迎風格式涉及到7個節點,離散準則需要進行6次通量分裂格式的計算,應用于TVD(Total Variation Diminishing)等二階格式涉及到13個節點,需要12次通量分裂;常規的有限差分法每個點上通量分裂格式只要1次,每個方向掃描1次,不管一階格式還是高階格式,通量分裂只需要計算3次即可。
為了消除坐標變換引起的數值誤差,離散準則確實付出了較高的計算效率為代價,不過利用差分格式中僅需要當地坐標變換系數的特性,可以在非結構網格點建立局部貼體坐標系從而采用有限差分格式進行計算。
為便于原理介紹,先針對二維情況進行討論。差分格式沿著網格線離散,從物理空間(x,y)映射到計算空間(ξ,η)至少需要離散點與周圍4個鄰點連成4條網格線,對于離散點0和周圍3個鄰點只能連成3條網格線,構成的非結構網格如圖1所示。

圖1 二維坐標變換對應關系
采用差商近似這些坐標變換系數,影響數值解的是網格坐標的相對位置。例如,采用二階中心差商求解,1-0-2三點連成曲線的斜率和1-2兩點直線段斜率一樣,因此把1-2線段平移到d-b不影響斜率。同樣處理1-0-3三點連成的曲線,得到對應的線段a-c。圍繞離散點0的a、b、c、d這4個點構成的網格可以得到坐標變換系數數值解,結果與1、2、3點構成的網格是一樣的。坐標變換系數的表達式為


(13)
二維算例表明流場沒有明顯不對稱現象,且網格的相鄰線段間最大夾角由120°變化到180°過程中,計算結果無明顯變化[7]。
三維空間中存在3條坐標軸,結構網格除了離散點外還需要6個相鄰點。由于三維幾何特性比較復雜,當4個相鄰點之間滿足如下條件:離散點0落在其他4個點構成的四面體內部,那么可以采用以上二維空間類似的方法進行映射,即在相鄰4點基礎上構造非結構網格并且進行有限差分法計算。圖2為0-1線段進行3次映射的示意圖,把曲線段1-0-4映射成計算空間ξ軸、0-1為正向,曲線1-0-3映射成計算空間η軸、0-3為正向,曲線1-0-2映射成計算空間ζ軸、0-1為正向。以變換矩陣行列式為例來說明坐標變換系數的計算方法:

圖2 三維非結構網格示意圖
(14)
在各點坐標ri=(xi,yi,zi)已知的情況下,采用中心差商計算式(14)中偏導數,例如:
xξ=0.5(x1-x4)
xη=0.5(x3-x1)
xζ=0.5(x1-x2)
這是對0-1線段多重映射形成的一種局部坐標系,考慮到映射在計算空間坐標軸方向特性,還應該有5種組合方式,總計需要6次計算。在非結構網格中4條連線等權重,組合以后存在24種局部坐標系,大大增加計算量,定性分析沒有實用價值。為了改善這種非結構網格有限差分法的計算效率,提出在直角坐標系下結構網格基礎上局部應用非結構有限差分法的扎染算法。
以二維半模球柱組合體為例,對扎染算法進行原理介紹。在包含物體在內的整個計算區域上分布直角坐標系下的均勻結構化網格,如圖3(a)所示,網格生成只需要一個控制參數Δx。根據網格點位置和物體幾何參數之間的相對關系,將已經生成的空間網格點分為3類:用黃色實心圓標識的是完全位于物體內部或物體表面附近點,定義為內部節點;用紅色矩形標識的是受物體影響的外部點,定義為非結構節點;沒有標識的是不受物體影響的外部點,定義為外部節點。標識內/外部節點需要設定判斷準則,本文采用與物體間距d<0.5Δx作為準則,結果如圖3(b)所示。

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

圖4 非結構節點連接關系
首先,在圖3(a)的均勻網格上進行計算,因為不需要Δx以外的其他網格相關信息,無需進行坐標變換,通量分裂計算表達式非常簡單,對于確定的差分格式,這種結構網格可以實現最高的計算效率。目前,本文算法只能采用時間顯式格式。定義全場標量cij,在圖3(b)示中非結構網格點和內部網格點上cij=0,其余點cij=1。由于cij=0點不影響物體外部流場,原則上可以不用計算。為避免跳過這些點需要的邏輯判斷,時間推進采用如下統一計算式:
(15)
下文中對是否進行跳點邏輯判斷進行了效率對比分析。
其次,僅在圖3(c)所示的混合網格區域內(紅色和藍色節點)進行計算,求解式(15)中cij=0的非結構節點和物面邊界節點參數,物體邊界對流場的影響在這一步體現。從圖3(c)可以看出,即使是非結構網格,離散點3條連線中只有1條特殊線,根據其所處象限確定唯一的映射關系或局部坐標系,只需計算1次,不像第2節中進行12次計算后算術平均。
混合網格區域內非結構節點與物面節點一一對應,很多工程應用中,物面節點數在總網格中所占的比例較小,所以大部分流場計算是在均勻正交網格下完成,計算效率高。物體內部對流場無貢獻,始終保持初始流場參數,如果云圖顯示整個流場,總有部分區域同色,如圖3(d)所示。本文的方法是一種直角坐標系均勻網格和局部曲線坐標系非結構網格的拼接算法,計算過程、流場圖像和中國傳統的扎染工藝有相似之處:根據設計圖形把部分白布包扎起來,直接投入染缸,操作簡單,從染缸出來后,為了精細化表現主題,保持變成藍底色部分不變,在其他區域進行局部人工彩繪,很多成品中還保留部分沒有上色的白布,所以稱其為扎染算法。
扎染算法應用于三維空間,點的分類過程完全相同:每個方向上根據距離物體d>0.5Δx準則進行分類,標識出內部節點,然后把相鄰點包含有內部節點的離散點標記為非結構點。在每個點的局部坐標系下差分離散是按照3個坐標軸6個節點進行的,只是離散點之間的連接關系按照非結構網格處理,因此,二維算法可以直接推廣至三維空間中。盡管三維問題增加了幾何復雜性,但是如果在每個坐標軸方向不多于1個內部節點,按照上面的算法也可以迅速實現:① 非結構節點連接1個內部點,延伸坐標軸到物面找點;② 非結構節點連接2個內部點,對沒有交點的坐標軸不需要特殊處理,退化為包含內部點的二維問題,根據最小距離準則尋找物面點,進行2次映射;③ 非 結構節點連接3個內部點,延伸坐標軸到物面3點,構成三角形平面,離散點和三角形面心或垂點的連線進行3次映射。對于連接內部點大于3個或同一坐標軸上有2個內部點的情況,根據工程經驗,這些情況下幾何構型非常復雜,要么需要加密網格來刻畫,要么需要修改簡化模型,總體來說需要人工干預,本文暫時不考慮這些特殊情況。
扎染算法計算流程簡述如下:
Step1讀取待計算物體模型(描述物體的散點集合或函數),給定計算區域邊界,在此基礎上對流場中所有節點屬性進行標識,確定結構節點、非結構節點、物面節點和內部節點。
Step2在包含物體的整個區域采用直角坐標系下的均勻網格進行計算,直接使用常規的有限差分方法進行求解,該過程中無需引入坐標變換,計算效率高。
Step3由于內部節點對流場計算不產生貢獻,且Step 2計算中沒有考慮物體的影響,因此將Step 2中內部節點和非結構節點的計算結果舍去。
Step4將非結構節點和物面節點一一對應,建立非結構網格“帶”,在此區域內使用非結構有限差分法進行求解,得到物面節點和非結構節點的流場參數。
Step5將Step 2中計算得到的結構網格節點與Step 4中得到的非結構網格節點及物面節點進行疊加,得到該時刻完整的流場信息。
Step6重復Step 1~Step 5進行時間推進求解,達到收斂條件,輸出流場結果。
選取二維/三維超聲速流動問題驗證本文算法。其中二維算例用于說明該方法的適用性及收斂穩定性,三維算例通過對比及理論分析展示扎染算法計算效率高的特點。
使用扎染算法對二維圓柱繞流問題進行數值模擬,驗證該方法在求解超聲速流動時的適用性。算例來流馬赫數Ma∞=3.0,無量綱流動參數(ρ,u,v,p)=(1.4,3.0,0,1.0),流場初始化為來流,時間離散采用Rung-Kutta法,空間離散使用一階迎風格式。無量綱計算域為x∈[0,0.1],y∈[0,0.4]的矩形區域,圓柱圓心位于(0.1,0.2),半徑R=0.045,分別選取5組不同Δx的網格進行數值模擬。定義5組中節點間距值最大的算例為Δx0=0.002 5,其余4組網格控制參數分別為Δx0/2、Δx0/4、Δx0/8以及Δx0/16。
5種網格條件下的收斂曲線如圖5所示,殘差單調下降,沒有出現異常波動,表明扎染算法可以得到穩定的收斂解,圖中t為無量綱時間。圖6(a)為5種網格條件下計算得到的對稱軸線上壓力P隨x變化曲線,物體表面壓力分布曲線如圖6(b)所示。從圖6(a)看出,隨著網格的不斷加密,對激波的分辨能力逐漸增強,對稱軸線上壓力間斷更加陡峭,峰值也相應增大。當進一步加密網格為Δx0/8和Δx0/16時,藍色曲線與綠色曲線重合,后續即使繼續加密網格,壓力曲線也不再發生明顯變化。從圖6(b)物面壓力分布曲線同樣可以看出上述結論,隨著網格加密,駐點處壓力峰值逐漸增大,當網格加密至一定程度后,壓力不再隨網格發生變化。在網格收斂性測試中通常采用沿滯止流線處駐點壓力來評估計算準確性,其理論值由Rayleign-pitot公式求得

圖5 殘差收斂曲線

圖6 不同網格尺度下壓力曲線
(16)

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

圖7 二維球頭繞流:一階迎風格式與WENO格式密度和壓力云圖對比
4.2.1 跳點法與覆蓋法效率對比
以三維球頭繞流為例,對扎染算法計算效率進行初步分析討論。計算模型如圖8所示,計算區域無量綱長度x=0.08、y=0.28、z=0.28,球頭半徑R=0.045。流場初始化及來流參數與4.1節二維算例相同。
圖8(b)為初始流場xy平面切片,圖中白色半圓內為物體區域,扎染算法求解過程中物體內部對流場計算無貢獻。對內部節點分別采用兩種處理方法,方法1表述為每一個時間步內對全場結構網格進行循環求解,后將對流場計算無貢獻的內部節點物理量賦為來流值,稱之為覆蓋法;方法1表述為循環過程中對結構網格的節點屬性進行判定,當節點位于物體內部時跳過該節點進行下一行列循環,內部節點不進行求解,稱之為跳點法。分別采用兩種處理方法求解三維球頭繞流,網格控制參數Δx=Δy=Δz=0.002,節點個數Nx=41,Ny=Nz=141。CFL數取值0.5,時間離散采用Runge-Kutta法,空間離散采用一階迎風格式,通量分裂使用Vanleer求解器。網格節點屬性如圖9所示,綠色標識為結構化節點,紅色標識為非結構節點,藍色標識為物面邊界節點,物體內部節點未顯示。

圖8 扎染算法計算模型

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

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

圖11 扎染算法CPU耗時對比

表1 節點屬性分布
4.2.2 常規差分算法與扎染算法計算效率分析
使用常規差分算法處理三維曲邊界問題時,在計算域中的所有網格點均需要由物理坐標系轉到貼體坐標系下進行格式離散。以本文采用的VanLeer格式為例,求解通量形式。定義ξ方向馬赫數:

(17)
式中:
(18)

(19)
式中:
(20)
式中:ρ、u、v、w、p分別表示密度、x方向速度(直角坐標系下)、y方向速度(直角坐標系下)、z方向速度(直角坐標系下)、壓力。
當|Maξ|<1(局部亞聲速)時,
(21)
式中:
(22)
(23)

(24)
其中:

(25)
從式(17)~式(25)中可知對于常規結構網格的有限差分法,通量求解過程中包含坐標變換系數,在增大計算量的同時內存占用也相應增加,計算效率較低。扎染算法求解過程中在遠離物面的空間計算區域內為均勻且正交的結構化網格,在此類節點求解時無需使用坐標變換系數,通量求解形式如下,以x方向超聲速流動條件為例,定義方向馬赫數:
(26)
當|Max|≥1時,
(27)
式中:
(28)
由式(26)~式(28)可知,均勻正交網格條件下,通量求解公式簡化,節省計算時間的同時占用內存更小,相較于常規結構網格有限差分算法計算效率更高,更加適合進行大規模并行計算。
本文以非結構網格有限差分方法為基礎,提出一種局部多重分塊網格的新算法。在均勻正交網格區域使用無坐標變換的差分算法,計算效率高;物面附近使用非結構有限差分法求解,可以有效處理復雜的拓撲結構。從上述算例演示可知,扎染算法適用于多種計算格式且滿足穩定性要求,通過對求解公式的對比分析表明扎染算法數據結構簡單,計算效率高,有利于開展大規模網格并行計算。
通過本文理論證明及算例驗證,定性分析了扎染算法的優勢及應用前景:
1)計算中大部分區域為均勻正交結構化網格,通量分裂中無需引入坐標變換系數,理論上可以達到最高的計算效率。
2)扎染算法應用中主要計算區域采用無需坐標變換的直角結構網格,使用差分方法進行求解,對網格分區無要求,適應于任意構架的高性能計算平臺。
3)由于物面節點是由空間非結構節點依據最小距離判據確定,因此僅需可以描繪外形的散點集合即可開展數值模擬,無需人工進行(表面/空間)網格生成,實現高度自動化的快速流場計算。
當前研究中針對具有解析表達式的簡單三維外形已經實現扎染算法計算,后續計劃與相關學者開展合作研究,將CAD等通用建模軟件與本文所提算法相結合,最終實現輸入實體輸出流場結果的快速計算,節省人力成本,提高計算效率。