求解歐拉方程的預(yù)處理隱式無網(wǎng)格算法

曹 騁，陳紅全

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

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

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

1 控制方程

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

(1)

其中W為守恒矢變量，E，F(xiàn)和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分別代表了密度、壓強(qiáng)、單位質(zhì)量的總能和總焓；u，v和w為速度分量. 完全氣體滿足狀態(tài)方程，有如下關(guān)系式成立：

(6)

ρH=ρE+p.

(7)

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

2 無網(wǎng)格空間導(dǎo)數(shù)近似

采用無網(wǎng)格算法，計(jì)算域先需布點(diǎn)離散. 布點(diǎn)離散后，即可形成由每一點(diǎn)臨近點(diǎn)構(gòu)成的當(dāng)?shù)攸c(diǎn)云結(jié)構(gòu)[6,8]. 記C(i)為第i點(diǎn)的點(diǎn)云，其構(gòu)成參見圖1. 圖中，點(diǎn)云的中心點(diǎn)即為i點(diǎn)，點(diǎn)云的其余點(diǎn)被稱作衛(wèi)星點(diǎn).

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

無網(wǎng)格空間導(dǎo)數(shù)近似是在計(jì)算域離散點(diǎn)上進(jìn)行的. 基于圖1所示的點(diǎn)云結(jié)構(gòu)，任一可微函數(shù)f的空間導(dǎo)數(shù)可近似為[8,27]

(8)

其中fik為第k個(gè)衛(wèi)星點(diǎn)與中心點(diǎn)i連線中點(diǎn)處的近似值. 顯然，系數(shù)αik，βik和γik滿足

(9)

上述空間導(dǎo)數(shù)近似式中的系數(shù)，可用加權(quán)的最小二乘曲面逼近來求得[8,19].

3 預(yù)處理隱式無網(wǎng)格算法

3.1 通量計(jì)算

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

(10)

其中i和k連線中點(diǎn)上的數(shù)值通量Lik具體可寫為

(11)

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

Uik=αikuik+βikvik+γikwik.

(12)

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

(13)

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

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

(14)

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

(15)

應(yīng)用(10)式，添加上人工耗散項(xiàng)后，得到方程(1)的半離散形式，具體可簡(jiǎn)寫為

(16)

3.2 預(yù)處理

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

(17)

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

Γ=

(18)

其中：

(19)

(20)

(21)

這里的Ma和Ma分別代表了當(dāng)?shù)睾蛠砹黢R赫數(shù).

應(yīng)用(18)式進(jìn)行矩陣預(yù)處理，方程(16)可改寫為

(22)

(23)

因此，預(yù)處理以后人工耗散項(xiàng)可改寫為

(24)

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

3.3 時(shí)間離散

依據(jù)Jameson和Yoon的做法[25]，對(duì)預(yù)處理以后的數(shù)值通量引入如下線性化處理:

(25)

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

(26)

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

(27)

其中Δti為點(diǎn)云中心點(diǎn)i的時(shí)間步長(zhǎng)，R為空間導(dǎo)數(shù)用點(diǎn)云近似的殘值，ΔW定義為

ΔW=Wn+1-Wn+1

(28)

因?yàn)橄禂?shù)αik，βik和γik具有式(9)的特性，雅克比矩陣Bik也應(yīng)該具有如下特性:

(29)

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

(30)

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

(31)

(32)

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

(33)

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

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

3.4 邊界條件

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

以遠(yuǎn)場(chǎng)邊界入流為例，虛擬點(diǎn)處的原始變量可取為

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

(34)

這里的下標(biāo)g為虛擬點(diǎn)，而下標(biāo)i為內(nèi)部流場(chǎng)點(diǎn)，下標(biāo)代表了自由流. 注意，這里g是i的鏡像點(diǎn).

需要指出的是，位于遠(yuǎn)場(chǎng)邊界外引入的鏡像虛擬點(diǎn)是在計(jì)算域布點(diǎn)離散過程中生成的. 這樣做，相當(dāng)于遠(yuǎn)場(chǎng)邊界條件位于虛擬點(diǎn)和對(duì)應(yīng)內(nèi)部點(diǎn)的中點(diǎn)處，而所有位于計(jì)算域內(nèi)與遠(yuǎn)場(chǎng)邊界條件相關(guān)的臨近邊界點(diǎn)都可以按照內(nèi)部點(diǎn)的方式處理. 可以看出，這樣處理的遠(yuǎn)場(chǎng)邊界條件實(shí)施更為簡(jiǎn)單方便.

4 算例與分析

本文用上述預(yù)處理隱式無網(wǎng)格算法，先對(duì)NACA 0012對(duì)稱翼型和AGARD三段翼型進(jìn)行了低馬赫數(shù)繞流的數(shù)值模擬，通過與文獻(xiàn)或?qū)嶒?yàn)結(jié)果的對(duì)比分析，驗(yàn)證所提算法的計(jì)算效率和準(zhǔn)確性，展示隱式無網(wǎng)格算法的快速收斂特性；接著給出了機(jī)翼及翼身組合體等跨音速和低馬赫數(shù)兩類繞流的算例，展示算法通過預(yù)處理對(duì)這兩類流動(dòng)模擬的兼容能力.

4.1 NACA 0012翼型低馬赫數(shù)繞流

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

圖2 NACA 0012翼型周圍點(diǎn)分布Fig.2 Points around NACA 0012 airfoil

圖3 預(yù)處理對(duì)收斂歷程的影響Fig.3 Effects of preconditioning on convergence histories

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

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

(a) 預(yù)處理隱式

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

(a) 帶預(yù)處理

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

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

4.2 AGARD三段翼型低馬赫數(shù)繞流

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

圖8 多段翼型周圍點(diǎn)分布 Fig.8 Points around multi-element airfoil

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

圖10 表面壓強(qiáng)系數(shù)分布Fig.10 Pressure coefficient distribution

圖11 馬赫數(shù)等值線和云圖Fig.11 Mach number contours

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

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

圖12 M6 機(jī)翼周圍點(diǎn)分布Fig.12 Points around ONERA M6 wing

(a)η=0.65

(b)η=0.9圖13 機(jī)翼表面壓強(qiáng)系數(shù)分布 (Ma=0.84)Fig.13 Pressure coefficient distribution(Ma=0.84)

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

圖14 機(jī)翼表面壓強(qiáng)系數(shù)分布(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 機(jī)翼上表面壓強(qiáng)系數(shù)等值線及云圖(Ma=0.01)

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

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

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

圖17 DLR-F4 翼身組合體周圍點(diǎn)分布Fig.17 Points around DLR-F4

圖18 收斂歷程Fig.18 Convergence histories

圖19 組合體表面壓強(qiáng)系數(shù)等值線云圖Fig.19 Pressure coefficient contours

5 結(jié) 論

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