求解歐拉方程的預處理隱式無網格算法

曹 騁，陳紅全

(南京航空航天大學 航空學院，南京 210016)

隨著現代計算機技術和數值分析的發展，計算流體力學(簡稱CFD)已取得了長足的進步，涌現出了可用于先進飛行器設計的各類CFD計算軟件[1-5]. 這類軟件大多選用網格單元進行計算區域的離散，因此，這類方法可統稱為“網格方法”，通常在流場數值模擬前都要求先對計算區域進行網格生成. 由于網格拓撲的約束，要對包含飛機全機等復雜幾何外形的計算區域生成合適的整體網格并非易事，尤其要求刻畫出多體小間隙等幾何細節時. 因此，如何局部或完全擺脫網格的約束問題，成為人們開始關注和研究的熱點[6-8]，出現了一類意在完全打破網格拓撲限制的算法，無網格算法. 這類算法對于處置和刻畫復雜幾何外形細節將更具靈活性，其理由主要有：計算區域是用無網格點云點布點離散的，具有靈活性和兼容性，既可按要求直接布點，也可利用已有的結構或非結構網格的網格點；無網格點上的空間導數近似僅依賴于當地點云點上的信息，而點云點無須連成網格. 許多學者正是受這些重要特征的驅動，開始對其研究，各具特色的無網格算法相繼涌現[9-14]. 就作者熟悉的空氣動力學領域而言，最引人注目的是早期Batina[6]的工作. 他采用中心差分格式，發展了顯式無網格算法，成功地求解了帶有激波的可壓縮繞流. 之后，Morinishi[11]采用虛擬中點(點云中心點與衛星點之間)迎風處理策略和加權最小二乘空間導數逼近，成功地發展了隱式無網格算法. 近年來，成功應用無網格算法的算例[15-18]已不難在文獻中找到，借助于現代圖形處理器(GPU)等并行加速技術，算法的收斂速度也進一步得到了強化[19-22]. 但同時可以注意到，上述無網格算法大多是基于可壓縮流發展的，不能直接拓展到模擬幾乎不可壓的低馬赫數流.

本文致力于發展出一種帶預處理的隱式無網格算法，既能模擬可壓跨音速流，又能模擬幾乎不可壓的低馬赫數流. 眾所周知，在低馬赫數下，特征型歐拉方程的3個不同特征值代表的波速存在巨大差異，會導致原控制方程產生系統病態，即所謂的“剛性”問題[23]. 而這反過來也會導致迭代計算收斂困難，數值模擬解往往偏離實際物理解. 為此，本文先沿用網格算法中著名的Weiss-Smith型預處理矩陣[24]，對無網格半離散型歐拉方程進行矩陣預處理，以期縮減方程特征值在低馬赫數條件下的差異，改善上述所謂的“剛性”問題；接著給出了無網格點上受預處理影響的人工耗散項及遠場邊界條件等實施細節，并通過無網格點排序與點云分割，實現了LU-SGS算法[25]的時間推進求解，形成了預處理隱式無網格算法. 最后，先選用典型翼型算例進行了數值模擬，并與文獻或實驗結果進行了驗證比較；接著給出了機翼及翼身組合體等三維跨音速和低馬赫數繞流算例. 算例表明所發展的隱式算法如預期比相應的顯式算法收斂更快，已從單純模擬可壓縮流動拓展到模擬幾乎不可壓的低馬赫數流.

1 控制方程

在直角坐標系(x,y,z)上，與時間t相關的三維歐拉方程的守恒形式為[26]

(1)

其中W為守恒矢變量，E，F和G為矢通量，可分別寫為

W=[ρ,ρu,ρv,ρw,ρE]T,

(2)

E=[ρu,ρu2+p,ρuv,ρuw,ρuH]T,

(3)

F=[ρv,ρvu,ρv2+p,ρvw,ρvH]T,

(4)

G=[ρw,ρwu,ρwv,ρw2+p,ρvH]T.

(5)

式中：ρ，p，E和H分別代表了密度、壓強、單位質量的總能和總焓；u，v和w為速度分量. 完全氣體滿足狀態方程，有如下關系式成立：

(6)

ρH=ρE+p.

(7)

其中，γ為流體的比熱比，對于空氣而言，γ=1.4. 通過求解上述方程組尋求滿足?W/?t=0的定常解.

2 無網格空間導數近似

采用無網格算法，計算域先需布點離散. 布點離散后，即可形成由每一點臨近點構成的當地點云結構[6,8]. 記C(i)為第i點的點云，其構成參見圖1. 圖中，點云的中心點即為i點，點云的其余點被稱作衛星點.

圖1 三維點云C(i)示意圖Fig.1 Three-dimensional cloud of points for C(i)

無網格空間導數近似是在計算域離散點上進行的. 基于圖1所示的點云結構，任一可微函數f的空間導數可近似為[8,27]

(8)

其中fik為第k個衛星點與中心點i連線中點處的近似值. 顯然，系數αik，βik和γik滿足

(9)

上述空間導數近似式中的系數，可用加權的最小二乘曲面逼近來求得[8,19].

3 預處理隱式無網格算法

3.1 通量計算

應用(8)式，歐拉方程的通量可近似寫為

(10)

其中i和k連線中點上的數值通量Lik具體可寫為

(11)

這里Uik與逆變速度定義類似，可定義為

Uik=αikuik+βikvik+γikwik.

(12)

假定中心點i和衛星點k連線的中點處存在一個虛擬的交界面，那么交界面上的通量Lik可以通過選擇合適的格式，從i點和k點處的守恒變量計算確定. 采用中心格式，交界面處的守恒變量可簡單地寫為

(13)

顯然，交界面上Lik可由Wik計算得到.

必須指出，上述中差計算形成的格式與網格算法中采用的中心差分格式類似，迭代計算會出現不穩定. 一般在激波這樣的強間斷附近會產生解的非物理振蕩，這種振蕩也可能因所謂的奇偶失聯而發生在任何無網格點上. 因此，本文沿用Jameson[28]在網格算法中的做法，在離散的方程中添加人工耗散項Disi以保證算法穩定. 基于無網格點云結構，人工耗散項可寫為

(14)

其中ε(2)和ε(4)為可調系數，λik則為雅克比矩陣Aik=?Lik/?Wik的譜半徑，其形式為

(15)

應用(10)式，添加上人工耗散項后，得到方程(1)的半離散形式，具體可簡寫為

(16)

3.2 預處理

可以注意到，上述交界面上歐拉方程的3個不同的特征值為

(17)

在低馬赫數下，如前述，這3個特征值代表的波速會存在巨大差異，雅克比矩陣Aik的條件數過大，會導致迭代計算收斂困難. 為此，本文考慮對方程(16)進行矩陣預處理以減小矩陣條件數. 從網格算法中可知，適用于低馬赫數流動的預處理矩陣[24,29-32]可有多種選擇，其中較為著名的是Weiss-Smith型預處理矩陣[24]. 對守恒型方程，該矩陣具體可寫為

Γ=

(18)

其中：

(19)

(20)

(21)

這里的Ma和Ma分別代表了當地和來流馬赫數.

應用(18)式進行矩陣預處理，方程(16)可改寫為

(22)

(23)

因此，預處理以后人工耗散項可改寫為

(24)

本文將采用LU-SGS隱式算法，對預處理以后的半離散方程(22)進行時間離散求解.

3.3 時間離散

依據Jameson和Yoon的做法[25]，對預處理以后的數值通量引入如下線性化處理:

(25)

這里通量雅克比矩陣Bik按文獻[25]中的分裂方式定義為

(26)

其中I為單位矩陣. 那么方程(22)的時間后差隱式格式可以寫為

(27)

其中Δti為點云中心點i的時間步長，R為空間導數用點云近似的殘值，ΔW定義為

ΔW=Wn+1-Wn+1

(28)

因為系數αik，βik和γik具有式(9)的特性，雅克比矩陣Bik也應該具有如下特性:

(29)

因而，隱式格式(27)可退化為

(30)

式(30)用LU-SGS算法[25]求解，格式可分為如下兩步:

(31)

(32)

其中，L(i)和U(i)為點云C(i)衛星點的LU分割.Di定義為

(33)

實際計算中，沿用Sharov和Nakahashi[33]的做法，涉及的矩陣運算Aik(Wk)ΔWk由通量的增量ΔLik(Wk)來簡化計算，減少工作量.

必須指出的是，L(i)和U(i)如果分割不當會導致算法退化為運算較慢的Jacobi迭代[33]. 為避免該問題，文獻[33]依據網格間的連接關系先進行網格分層編號，再進行LU分割. 層號比當前網格層號小的相鄰網格單元被歸入L(i)，反之歸入U(i). 對于無網格算法，點云點間不存在網格拓撲這樣的剛性連接關系，為此，本文利用衛星點與中心點的虛擬連接關系來構造類似的分層結構. 作者曾提出一種多層點云的重排算法，優化了GPU的內存訪問模式， 實現了GPU并行無網格算法的進一步加速[22]. 本文沿用這一方法對無網格點進行類似的分層編號重排. 分層后，對點云C(i)的衛星點進行LU分割. 將層號比中心點i層號小的衛星點分割納入L(i)，反之則納入U(i). 本文將這一分割做法用于了后續的算例運算.

3.4 邊界條件

對于物面邊界條件，無粘條件下滿足的無穿透邊界條件[6]在預處理前后都一樣. 而遠場基于特征分析的無反射邊界條件，預處理后由于特征值的改變，必須作相應的修正. Turkel[34]早期給出了精確的基于特征分析的預處理遠場邊界條件，之后又對其進行了簡化，用于處理幾乎不可壓流動. 本文沿用Turkel的簡化處理[32]，預處理遠場邊界條件是結合無網格布點離散靈活的特點實施的. 為此，本文把遠場邊界面處理成臨近內點的鏡面，每一離鏡面最近的內點對應一個虛擬的鏡像點，納入該點的點云中. 那么，虛擬點上的物理量可直接根據Turkel的方法確定.

以遠場邊界入流為例，虛擬點處的原始變量可取為

ρg=ρ,ug=u,vg=v,wg=w,pg=pi.

(34)

這里的下標g為虛擬點，而下標i為內部流場點，下標代表了自由流. 注意，這里g是i的鏡像點.

需要指出的是，位于遠場邊界外引入的鏡像虛擬點是在計算域布點離散過程中生成的. 這樣做，相當于遠場邊界條件位于虛擬點和對應內部點的中點處，而所有位于計算域內與遠場邊界條件相關的臨近邊界點都可以按照內部點的方式處理. 可以看出，這樣處理的遠場邊界條件實施更為簡單方便.

4 算例與分析

本文用上述預處理隱式無網格算法，先對NACA 0012對稱翼型和AGARD三段翼型進行了低馬赫數繞流的數值模擬，通過與文獻或實驗結果的對比分析，驗證所提算法的計算效率和準確性，展示隱式無網格算法的快速收斂特性；接著給出了機翼及翼身組合體等跨音速和低馬赫數兩類繞流的算例，展示算法通過預處理對這兩類流動模擬的兼容能力.

4.1 NACA 0012翼型低馬赫數繞流

本文選用NACA 0012翼型幾乎不可壓的低馬赫數繞流對發展的算法進行了考核計算驗證. 如圖2所示，計算域用3343個無網格點進行離散. 先用馬赫數Ma=0.001，攻角為0°的繞流模擬展示了預處理對收斂歷程的影響(見圖3). 圖中可看出，發展的預處理隱式算法，在如此小的馬赫數下仍能快速收斂，圖中同時給出了不加預處理的收斂困難的情形供比較. 對應的翼型表面壓強系數分布(見圖4)取得了與文獻[35]和實驗[36]一致的結果. 圖中橫坐標X/C代表離前緣的相對距離，C為對應翼型的弦長，縱坐標Cp則為壓強系數. 圖5則給出了對應的馬赫數等值線和云圖. 圖中反映的流場對稱性與零攻角對稱翼型繞流的物理特性吻合. 接著，通過不同馬赫數的數值模擬展示了馬赫數對算法收斂歷程的影響. 從圖6中可看出，模擬的4個來流馬赫數(0.3，0.1，0.05，0.01)都能很好收斂. 圖中同時給出了顯式預處理算法[23]的收斂歷程供比較. 必須指出，在幾乎不可壓的低馬赫數下，很有必要使算法具有良好收斂特性. 這一點從捕捉的流場解可以看出. 圖7給出了隱式無網格算法帶與不帶預處理計算獲得的流場密度等值線(Ma=0.05)比較. 圖中可看出，預處理后流場解更加合理，表現為等值線更加光滑.

圖2 NACA 0012翼型周圍點分布Fig.2 Points around NACA 0012 airfoil

圖3 預處理對收斂歷程的影響Fig.3 Effects of preconditioning on convergence histories

圖4 翼型表面壓強系數分布(Ma=0.001)Fig.4 Pressure coefficient distribution(Ma=0.001)

圖5 馬赫數等值線及云圖(Ma=0.001)Fig.5 Mach number contours (Ma=0.001)

(a) 預處理隱式

(b) 預處理顯式圖6 馬赫數對收斂歷程的影響Fig.6 Effects of Mach numbers on convergence histories

(a) 帶預處理

(b) 不帶預處理圖7 隱式無網格算法密度等值線比較(Ma=0.05)

Fig.7 Density contours of implicit gridless methods(Ma=0.05)

4.2 AGARD三段翼型低馬赫數繞流

這里給出外形更為復雜的AGARD三段翼型[37]幾乎不可壓的低馬赫數繞流的計算結果. 幾何外形可從圖8中看出. 計算域用5865個無網格點進行離散. 來流馬赫數按實驗條件[37]取為0.197，迎角取為4.01°. 如圖9收斂歷程所示，發展的隱式預處理無網格算法再次呈現出快速收斂特性，圖中同時給出了隱式不帶預處理、顯式帶與不帶預處理無網格算法的收斂歷程供比較. 對應的翼面壓強系數分布(見圖10)與實驗[37]吻合的較好，這一點可從主翼上捕捉的激波位置和強度看出. 圖11則給出了對應的馬赫數等值線和云圖，展示出了多段翼型的流動特征.

圖8 多段翼型周圍點分布 Fig.8 Points around multi-element airfoil

圖9 不同無網格算法收斂歷程比較Fig. 9 Convergence histories of different methods

圖10 表面壓強系數分布Fig.10 Pressure coefficient distribution

圖11 馬赫數等值線和云圖Fig.11 Mach number contours

4.3 M6機翼的跨音速和幾乎不可壓繞流

M6機翼[6,38]經常被用來考核發展的算法. 本文先進行了繞M6機翼跨音速流動的數值模擬. 計算采用了54320個無網格點(見圖12)，來流馬赫數為0.84，攻角為3.06°. 典型機翼翼剖面的壓強系數分布見圖13. 圖中η代表機翼翼剖面展向相對位置，翼尖η為1，翼根η為0. 計算結果與文獻[6]和實驗[39]進行了比較. 如預期，帶預處理的隱式無網格算法仍能用于該跨音速繞流計算，表現為計算結果與不帶預處理的結果十分吻合.

圖12 M6 機翼周圍點分布Fig.12 Points around ONERA M6 wing

(a)η=0.65

(b)η=0.9圖13 機翼表面壓強系數分布 (Ma=0.84)Fig.13 Pressure coefficient distribution(Ma=0.84)

接著，沿用文獻[38]的做法，進行了幾乎不可壓馬赫數繞流(Ma=0.01)的考核運算. 圖14給出了計算獲得的η=0.8處的翼剖面壓強系數分布，從圖中可以看出，發展的預處理隱式無網格算法與預處理顯式算法一樣，計算結果都與文獻[38]結果十分吻合. 比較而言，隱式預處理無網格算法較顯式收斂特性更優(見圖15). 圖16則給出了對應的上翼面壓強系數等值線云圖，可以看出捕捉的等值線光滑有序，一定程度上反映出M6機翼上翼面幾乎不可壓的低速流動特征.

圖14 機翼表面壓強系數分布(Ma=0.01，η=0.8)Fig.14 Pressure coefficient distribution (Ma=0.01, η=0.8)

圖15 收斂歷程(Ma=0.01)Fig.15 Convergence histories(Ma=0.01)

圖16 機翼上表面壓強系數等值線及云圖(Ma=0.01)

Fig.16 Pressure coefficient contours on the upper surface of the wing (Ma=0.01)

4.4 DLR-F4翼身組合體的低馬赫數繞流

面向實際應用，這里給出DLR-F4組合體低馬赫數繞流的演示算例. 眾所周知，馬赫數0.3在工程中常作為判斷一個流動是可壓與不可壓的參考馬赫數. 必須指出，按照可壓和不可壓的這一區分標準，對于來流馬赫數為0.3的翼身組合體繞流，其流場存在駐點附近等局部低馬赫數區和機翼上翼面等局部高馬赫數區，一定程度上可理解為可壓和不可壓并存的混合流場. 不難看出，隨著來流馬赫數的降低，流場中局部低馬赫數的不可壓區在整個流場中的占比將逐漸增大，這會一定程度上影響到迭代格式的收斂特性. 本文選用這一特別的來流馬赫數對發展的算法進行了收斂歷程測試. 計算采用了75060個無網格點(見圖17). 圖18則給出了對應的收斂歷程，圖中可以看出，加預處理可以明顯改善低馬赫數來流的收斂特性，一定程度上有助于提高解的捕捉精度(參見表面壓強系數等值線云圖19). 圖18中同時給出了來流馬赫數0.2的收斂歷程供比較，也可以看出，隨著來流馬赫數的減小，流場中低馬赫數的不可壓區占比擴大，如預期，加預處理的必要性更加顯現，表現為加與不加預處理，收斂差異更大.

圖17 DLR-F4 翼身組合體周圍點分布Fig.17 Points around DLR-F4

圖18 收斂歷程Fig.18 Convergence histories

圖19 組合體表面壓強系數等值線云圖Fig.19 Pressure coefficient contours

5 結 論

本文通過無網格點重排與點云分割，成功地結合LU-SGS算法，發展出一種基于Weiss-Smith預處理矩陣的預處理隱式無網格算法，用于求解三維歐拉方程. 借助于無網格布點，遠場邊界虛擬鏡像點的引入，簡化了預處理遠場邊界的實施. 發展的算法通過了二維和三維多個算例的計算考核. 算例表明，相較于顯式預處理算法，發展的預處理隱式算法收斂特性更優；算法的數值模擬能力已從單純模擬可壓縮跨音速等流動拓展到模擬幾乎不可壓的低馬赫數流. 此外，算法整體是基于無網格技術發展的，計算域只需布點離散，有助于靈活處理三維復雜氣動外形. 當然，實際工程繞流問題大多需考慮粘性的影響，這就要求發展的算法從求解歐拉方程拓展到求解Navier-Stokes方程，會涉及反映粘性繞流特征的異性點云結構的構造、合適湍流模型的嵌入、基于點云結構的粘性預處理譜半徑的確定等問題，將是未來研究的重點.