喬磊， 白俊強(qiáng)， 邱亞松,*， 華俊,2， 張揚(yáng)

(1. 西北工業(yè)大學(xué) 航空學(xué)院, 西安 710072； 2. 中國航空工業(yè)集團(tuán)有限公司 中國航空研究院, 北京 100012；3. 西安交通大學(xué) 航天航空學(xué)院， 西安 710049)

在飛行器和風(fēng)力機(jī)械氣動特性的評估和優(yōu)化設(shè)計(jì)中，需要進(jìn)行大量的定常流動模擬。在非定常數(shù)值模擬中，需要用定常流動的解作為初始解。因此，提高定常流場數(shù)值模擬的效率，可以使計(jì)算流體力學(xué)方法在氣動外形設(shè)計(jì)應(yīng)用中發(fā)揮更有力的作用。目前，常用的定常流動問題的控制方程是定?？蓧嚎s雷諾平均納維-斯托克斯(Reynolds Averaged Navier-Stokes, RANS)方程，此方程是一個(gè)由對流項(xiàng)主導(dǎo)的非線性偏微分方程組，其求解分為顯式和隱式2種方法。由于定常問題不需要考慮時(shí)間精度，隱式解法得以充分發(fā)揮其穩(wěn)定性強(qiáng)、計(jì)算效率高的優(yōu)點(diǎn)，因而得到了廣泛的應(yīng)用。

隱式解法的實(shí)質(zhì)是非線性方程的牛頓迭代解法。牛頓迭代的優(yōu)點(diǎn)是可以實(shí)現(xiàn)快速的平方收斂，但是要求初值足夠接近方程的解才能避免發(fā)散。在實(shí)際應(yīng)用中，一般情況下在迭代開始前是無法獲得足夠接近真實(shí)解的初值的。在外流流場模擬中，定常計(jì)算的初值通常為按自由來流條件設(shè)定的一個(gè)均勻場，直接應(yīng)用牛頓解法往往會導(dǎo)致求解發(fā)散。因此，研究者們發(fā)展了多種類型的牛頓迭代全局化方法。

針對非線性問題的牛頓迭代全局化方法可以分為2種類型。第一類是不涉及非線性問題本身的方法。在這類方法中，比較常見的有網(wǎng)格序列法[1-2]和線性搜索法[3-4]。這類方法不是本文研究的重點(diǎn)，故不做深入介紹。因?yàn)檫@些方法獨(dú)立于非線性問題本身，所以可以與本文所討論的方法疊加使用，構(gòu)造進(jìn)一步提高求解效率的方法。另外一類方法是與非線性問題本身相關(guān)的。這類方法是通過修改所要求解的問題，保證迭代的收斂。針對定常流場模擬問題，這類方法的典型代表有邊界松弛(boundary condition relaxation)法和方程延拓(equation continuation)法[5]。在均勻初始解中，求解域內(nèi)部殘差全部為零。在物面邊界處，由于邊界條件(無滑移或無穿透條件)的存在，會產(chǎn)生一個(gè)階躍型非零殘差。針對這一特點(diǎn)，邊界松弛法通過逐步施加邊界條件，降低迭代的不穩(wěn)定性。Lyra[6]和Kuzmin[7-8]等對邊界松弛法有較多的研究和應(yīng)用。

在方程延拓法中，最為常見的就是經(jīng)典的偽時(shí)間推進(jìn)法(Pseudo Time Marching，PTM)。方程延拓法是指在控制方程中引入額外的項(xiàng)，得到性質(zhì)改善的近似方程。然后逐次求解拓項(xiàng)系數(shù)不斷減少的近似方程，逼近原方程的解。由于這種延拓處理不改變非線性迭代的起點(diǎn)和終點(diǎn)，所以此類延拓方法又被稱為同倫延拓。偽時(shí)間推進(jìn)法作為一種方程延拓法，在定?？刂品匠讨幸肓艘粋€(gè)時(shí)間項(xiàng)，隨著求解的收斂時(shí)間項(xiàng)的作用逐漸消失。對于隱式方法，偽時(shí)間項(xiàng)使控制方程雅可比矩陣的非奇異性和對角占優(yōu)得到增強(qiáng)，保證了隱式迭代穩(wěn)定性。作為一種經(jīng)典方法，偽時(shí)間推進(jìn)法得到了非常廣泛的應(yīng)用，也有很多研究和發(fā)展。Coffey等[9]提出了一種差分代數(shù)方程形式的時(shí)間推進(jìn)法，提高了可壓縮燃燒問題的收斂速度。Kelley等[10]提出了一種帶約束的偽時(shí)間推進(jìn)法，提高了迭代的穩(wěn)定性。Ceze和Fidkowski[11-12]提出了一種針對非物理解的罰函數(shù)法，以提高偽時(shí)間推進(jìn)法的魯棒性，降低迭代發(fā)散的幾率，從而提高計(jì)算效率。

在方程延拓法中，還有一類基于黏性或人工耗散的方法。這種方法在計(jì)算流體力學(xué)發(fā)展的早期就被用于定常流場的計(jì)算。Young等提出了用于加速全速勢方程收斂的黏性松弛法[13]。Hicken等對比了黏性延拓和偽時(shí)間推進(jìn)法在定常RANS方程求解中的效率和穩(wěn)定性，認(rèn)為“黏性延拓法”具有較高的魯棒性和效率，是偽時(shí)間推進(jìn)法的一種可能的替代方法”[14-15]。在計(jì)算流體力學(xué)領(lǐng)域以外，Pollock發(fā)表了一種針對線性對流輸運(yùn)方程的，通過拉普拉斯算子解決由于解中不光滑成分而導(dǎo)致病態(tài)雅可比的問題[16]。黏性延拓法中還有一種雷諾數(shù)延拓法[17]。這種方法非常便于實(shí)現(xiàn)，需要做的就是從一個(gè)較低的雷諾數(shù)開始定常迭代，這樣較大的物理黏性會使問題具有足夠的耗散從而保持穩(wěn)定。根據(jù)Hicken和Zingg的研究結(jié)論[15]，這種方法存在一個(gè)不足，就是Navier-Stokes方程的物理黏性不涉及連續(xù)方程，所以增穩(wěn)作用比較有限。

通常外流流場定常模擬的迭代過程就是把繞流物體的壁面邊界對流動的影響傳播到流場內(nèi)部的過程。從這方面考慮，黏性延拓法有助于改善非線性求解部分的效率。然而，在牛頓迭代中，延拓項(xiàng)會對雅可比矩陣的性質(zhì)產(chǎn)生直接影響，從這個(gè)角度看，偽時(shí)間法帶來的較強(qiáng)的對角占優(yōu)特性又相對更具有優(yōu)勢。為綜合利用兩者的優(yōu)勢，本文提出了一種拉普拉斯增穩(wěn)偽時(shí)間推進(jìn)(Laplacian stabilized Pseudo Time Marching, LPTM)法。通過引入拉普拉斯算子增加方程的穩(wěn)定性，提高可用的CFL數(shù)，達(dá)到加速收斂、節(jié)約計(jì)算時(shí)間的目的。

本文首先介紹了控制方程和牛頓迭代法等工作基礎(chǔ)，然后詳細(xì)闡述了拉普拉斯增穩(wěn)的迭代思路，分析其與經(jīng)典偽時(shí)間推進(jìn)法的優(yōu)缺點(diǎn)，并給出了LTPM法的實(shí)現(xiàn)細(xì)節(jié)，最后，通過無黏NACA0012翼型、湍流RAE2822翼型和三維ONERA M6機(jī)翼3個(gè)算例，對LPTM法的收斂加速效果進(jìn)行了驗(yàn)證。

1 控制方程及其隱式解法

本文計(jì)算采用格心格式有限體積法求解可壓縮RANS方程，無黏通量通過三階MUSCL重構(gòu)的Roe格式計(jì)算[18]，無黏通量采用二階中心格式離散[18],隱式時(shí)間推進(jìn)中采用對稱高斯-賽德爾(Symmetric Gauss-Seidel, SGS)迭代預(yù)處理的廣義最小殘差法(Generalized Minimal RESidual method，GMRES)求解線性子問題，并使用多重網(wǎng)格技術(shù)加速收斂。程序通過基于MPI的分布式并行策略提高計(jì)算速度。本文在計(jì)算中使用的湍流模型為SA(Spalart-Allmaras)一方程湍流模型[19]。

1.1 控制方程

本文所研究問題的控制方程是無量綱的定常可壓縮RANS方程，表示為

(1)

式中：R(w)為殘差，流體狀態(tài)向量w=(ρ,U,p)，ρ、U和p分別為密度、速度和壓力的流場自變量；F(w)和Fv(w)分別為無黏和黏性通量；Ω為求解域或控制體；v為積分變量?？刂品匠讨械耐牧黟ば皂?xiàng)需要通過求解湍流模型得到。在本文的討論中，湍流模型的求解是和流動控制方程獨(dú)立的。因此，本文不涉及關(guān)于湍流模型求解的問題。

1.2 偽時(shí)間推進(jìn)法的基本作用

方程式(1)是關(guān)于w的非線性方程。在非線性迭代法中，牛頓迭代由于其二階收斂能力而得到廣泛應(yīng)用。關(guān)于方程式(1)的牛頓迭代可以表示為

(2)

式中：n為迭代步數(shù)；R′為非線性算子的雅可比矩陣；δw為解向量w的增量。

但是牛頓迭代的收斂對方程的性質(zhì)和初始解有較苛刻的要求：必須滿足雅可比矩陣非奇異，初始解足夠接近方程解。所謂“足夠接近”，具體是指對于任意第n步牛頓迭代，須滿足

(3)

式中：R″為非線性算子的海森矩陣；ε為當(dāng)前解的誤差。因此，增強(qiáng)牛頓迭代收斂性的一般思路是降低問題雅可比矩陣的奇異性，以及增強(qiáng)問題的線性。作為一種常見的特例，經(jīng)典的全局化方法是偽時(shí)間推進(jìn)法，通過引入偽時(shí)間項(xiàng)，所求問題轉(zhuǎn)化為

(4)

式中：RPT(w)為包含偽時(shí)間項(xiàng)的總殘差；Q(w)為守恒變量；τ為偽時(shí)間。

在實(shí)際計(jì)算中，Q(w)可以取為Navier-Stokes方程守恒變量，也可以簡單取w本身。出于簡化分析的考慮，本文取w。這樣，牛頓迭代轉(zhuǎn)化為

(5)

式中：dτ為偽時(shí)間步長。由于對角陣I的引入，所解問題的線性增強(qiáng)，式(4)中的雅可比矩陣在對角占優(yōu)和條件數(shù)方面相對于式(2)有所改善，使牛頓迭代得以收斂。

2 拉普拉斯延拓

2.1 拉普拉斯算子的引入

1.2節(jié)提到，要保證牛頓迭代的收斂，基本策略就是降低雅可比矩陣的奇異性，增加問題的線性，使近似解盡快靠近問題真實(shí)解。要滿足這些要求，偽時(shí)間推進(jìn)法并不是唯一的選擇。實(shí)際上，任意非奇異的線性算子，都能起到減少海森矩陣范數(shù)、增加雅可比矩陣范數(shù)的作用。因此，本文考慮引入另外一種線性算子——拉普拉斯算子作為延拓項(xiàng)。引入拉普拉斯項(xiàng)的控制方程如下：

(6)

式中：RLPT(w)為由流體狀態(tài)向量w計(jì)算得到的殘差；cLP為拉普拉斯項(xiàng)的縮放系數(shù)；Δ為拉普拉斯算子。這樣，牛頓迭代公式變?yōu)?/p>

(7)

需要特別說明的一點(diǎn)是，在具體實(shí)現(xiàn)中，式(5)和式(7)的殘差矢量都是由原始方程得到的，并不包含由延拓項(xiàng)產(chǎn)生的殘差。根據(jù)作者的經(jīng)驗(yàn)，這種選擇得到的迭代格式比嚴(yán)格用牛頓迭代求解延拓后的方程具有更高的收斂效率。

另外，在黏性計(jì)算中，速度場在無滑移邊界附近會形成邊界層。此時(shí)，若對動量方程增加拉普拉斯算子，則會較大程度地?cái)U(kuò)散由無滑移邊界條件引起的低速區(qū)，反而不利于收斂。因此，對于存在無滑移邊界的問題，本文只對密度和壓力等受諾依曼邊界條件約束的分量施加拉普拉斯項(xiàng)。

2.2 拉普拉斯延拓作用

偽時(shí)間項(xiàng)向雅可比矩陣中加入了一個(gè)對角陣，拉普拉斯項(xiàng)在雅可比矩陣中引入了一個(gè)拉普拉斯算子。兩者都是對稱正定的線性算子，在改善控制方程正則性方面具有類似的作用。不過拉普拉斯逆算子與偽時(shí)間項(xiàng)逆算子相比，一個(gè)重要的特點(diǎn)是全局性。在實(shí)際計(jì)算中，數(shù)值解的誤差以及雅可比矩陣和海森矩陣的范數(shù)都難以得到，延拓算子與方程本身復(fù)合后的逆算子的具體形式更是難以求出。因此，此處僅對拉普拉斯算子和偽時(shí)間項(xiàng)的區(qū)別做定性分析。

由于偽時(shí)間項(xiàng)是一個(gè)對角陣，對角陣的逆算子仍然是對角陣，顯然它對殘差矢量的響應(yīng)是當(dāng)?shù)氐?、局部的。而拉普拉斯算子則不同。在三維歐氏空間中，通過格林函數(shù)法[20]可知，對于殘差向量R(w(x))，在拉普拉斯逆算子作用下的響應(yīng)為

(8)

這意味著任意一點(diǎn)的殘差都會對解向量產(chǎn)生全局性的影響；同時(shí)任意一點(diǎn)的解向量變化，都包含了整個(gè)求解域中的殘差信息。這正是拉普拉斯算子橢圓性的體現(xiàn)?？紤]到定常流場求解通常初始化為均勻來流狀態(tài)，初始?xì)埐顑H在壁面處不為零，則在計(jì)算初始階段，拉普拉斯算子向流場內(nèi)部傳遞邊界信息的效果要遠(yuǎn)遠(yuǎn)強(qiáng)于偽時(shí)間項(xiàng)，解的殘差也就可以更快地降低。這樣，由式(3)可知，牛頓迭代的穩(wěn)定性也就更容易保證，對CFL數(shù)的要求也就可以放得更寬，最終可以通過使用較大的CFL數(shù)，得到更快的收斂速度。

但是，拉普拉斯延拓相對于偽時(shí)間推進(jìn)法有2個(gè)方面的不足。第一，偽時(shí)間項(xiàng)在解逐步收斂時(shí)會自動消失，而拉普拉斯算子的作用在不加特殊處理的情況下會一直存在。這樣，即使選用的較小的CFL數(shù)，其負(fù)面影響也只是局限于使收斂速度變慢，只要進(jìn)行充分的迭代，仍能得到原始問題的解。這一特點(diǎn)給實(shí)際應(yīng)用帶來極大便利。而對于拉普拉斯延拓，必須特別設(shè)計(jì)一種機(jī)制，使得拉普拉斯項(xiàng)隨迭代推進(jìn)而消失。不過，這是個(gè)很容易克服的缺點(diǎn)。并且，在偽時(shí)間推進(jìn)法中，為追求較高的收斂效率，通常也不會使用單一CFL數(shù)進(jìn)行計(jì)算，同樣涉及到在迭代過程中改變CFL數(shù)的問題。第二，正如式(8)所示，拉普拉斯算子的逆是具有全局影響的，反映在離散系統(tǒng)中，就意味其值對應(yīng)的矩陣不再具有稀疏性。這樣，基于近似LU分解的線性求解方法，包括ADI、LU-SGS以及本文使用的SGS方法，效率就會有所下降，從而導(dǎo)致線性問題求解的代價(jià)增加。在實(shí)際應(yīng)用中，需要綜合權(quán)衡拉普拉斯算子的正負(fù)作用，才能得到整體效率更優(yōu)的求解策略。

2.3 迭代策略

本文的求解策略包含3層迭代。最外層是延拓迭代，隨著迭代的進(jìn)行，要保證延拓項(xiàng)不斷降低。第2層是牛頓迭代，每個(gè)延拓方程仍然是一個(gè)非線性問題，本文通過牛頓迭代法進(jìn)行求解。由于延拓問題不是最終關(guān)心的問題，所以牛頓迭代并不需要精確進(jìn)行。本文針對每個(gè)延拓問題，只進(jìn)行一步牛頓迭代。第3層迭代是針對牛頓迭代產(chǎn)生的線性問題，本文使用SGS預(yù)處理的GMRES方法進(jìn)行求解。同樣，由于牛頓迭代只是近似求解，線性求解同樣也不需要嚴(yán)格進(jìn)行。本文的策略是線性殘差降低2個(gè)數(shù)量級時(shí)認(rèn)為GMRES迭代收斂。

延拓項(xiàng)的推進(jìn)策略使用成熟的CFL數(shù)遞增策略(Switched Evolution Relaxation, SER)[21-22]。本文對cLP和CFL數(shù)的導(dǎo)數(shù)采用同樣的遞減模式，如式(9)所示。

(9)

(10)

延拓參數(shù)的選擇對計(jì)算的穩(wěn)定性和效率有至關(guān)重要的影響。較大的延拓參數(shù)(對于時(shí)間推進(jìn)法對應(yīng)較小的CFL數(shù))會使求解穩(wěn)定但是收斂較慢，反之亦然。由于控制方程的非線性性質(zhì)，對延拓參數(shù)數(shù)值的選擇進(jìn)行先驗(yàn)的理論分析較為困難，因此參數(shù)數(shù)值的選擇具有一定的經(jīng)驗(yàn)性。但是，這一問題在實(shí)際應(yīng)用中是比較容易克服的。在氣動外形優(yōu)化工作中，通常需要對相似的計(jì)算狀態(tài)進(jìn)行數(shù)十、數(shù)百乃至上千次重復(fù)計(jì)算。在這種情況下，通過對目標(biāo)狀態(tài)進(jìn)行數(shù)次試算，選取相對高效的CFL數(shù)和拉普拉斯項(xiàng)參數(shù)是完全可行的。

3 算例驗(yàn)證

3.1 計(jì)算效率的比較方法

本文通過計(jì)算收斂所消耗的CPU時(shí)間t衡量計(jì)算方法的效率。因此，首先需要保證計(jì)算硬件環(huán)境的一致。具體地，二維算例是在1個(gè)Intel Xeon E5-2 620 v3 2.40 GHz CPU上通過單進(jìn)程計(jì)算，三維M6算例是在2個(gè)同樣的CPU上通過10個(gè)MPI進(jìn)程并行計(jì)算。其次要統(tǒng)一收斂判斷標(biāo)準(zhǔn)。由于本文牛頓迭代所用的是近似雅可比矩陣，所以這里不追求控制方程非線性殘差收斂到數(shù)值極限，而是以氣動設(shè)計(jì)中常用的升力和阻力系數(shù)的收斂為標(biāo)準(zhǔn)。文獻(xiàn)[23]指出一般工程問題對氣動力的要求為：升力系數(shù)CL精確到0.001，阻力系數(shù)CD精確到0.000 1。本文采用上述容差的一半作為氣動力收斂判斷標(biāo)準(zhǔn)。具體地，針對最后10步非線性迭代，如果升力系數(shù)的變化范圍小于0.5×10-3，并且阻力系數(shù)的變化范圍小于0.5×10-4，則認(rèn)為計(jì)算收斂。

(11)

需要說明的是，選取最后10步迭代作為觀察范圍是不具有一般性的，特別是在CFL比較小的計(jì)算中，很容易得到虛假收斂判定。因此，出于嚴(yán)謹(jǐn)性的考慮，對于本文所涉及的算例，在達(dá)到這一收斂判據(jù)后僅記錄當(dāng)時(shí)的CPU時(shí)間，計(jì)算迭代仍繼續(xù)進(jìn)行，后續(xù)的計(jì)算結(jié)果可以作為對收斂性進(jìn)行人工輔助判斷的依據(jù)。如果后續(xù)計(jì)算沒有出現(xiàn)CL和CD的明顯漂移或波動，則認(rèn)為此收斂判定有效。實(shí)際結(jié)果表明，這一收斂判據(jù)對本文涉及的算例來說是合理準(zhǔn)確的，以此為基礎(chǔ)的計(jì)算效率對比也是有效的。

3.2 無黏NACA0012翼型

本算例考察了LPTM法對二維無黏問題的收斂加速作用。計(jì)算狀態(tài)為來流馬赫數(shù)0.8，迎角1.25°。計(jì)算網(wǎng)格總量為44 640單元，物面第一層網(wǎng)格高度為1.0×10-3，法向增長率為1.2。網(wǎng)格的大致結(jié)構(gòu)和疏密分布如圖1所示，Y/C和X/C分別為以弦長C無量綱化的縱、橫坐標(biāo)。

圖1 無黏NACA0012翼型算例的計(jì)算網(wǎng)格

算 例CFL0c0LPnt/s相對時(shí)間節(jié)約/%PTM1615238.79-28.2PTM2810930.250LPTM1205×10-26725.7714.8LPTM2205×10-36622.0727.0LPTM3201×10-36422.5025.6LPTM4205×10-46221.3729.4LPTM5105×10-59225.5615.5

圖2 無黏NACA0012翼型算例升力系數(shù)、阻力系數(shù)及殘差收斂曲線

圖2(a)、(b)給出了無黏NACA0012翼型算例升力系數(shù)和阻力系數(shù)對計(jì)算CPU時(shí)間的收斂曲線。可以看出LPTM法的力系數(shù)收斂有明顯的提前，并且振蕩的幅度也有所降低。圖2(c)給出了系統(tǒng)殘差隨迭代步數(shù)的收斂過程。由圖2可見，收斂曲線依CFL數(shù)分為4簇，顯示了CFL數(shù)對收斂效率具有直接決定性作用。而LPTM算例之所以能使用更高的CFL數(shù)，原因是拉普拉斯算子的引入。在最初的幾步迭代中，拉普拉斯算子的引入可以使殘差更快收斂，驗(yàn)證了2.2節(jié)中對拉普拉斯項(xiàng)作用的分析。

圖3給出了本文計(jì)算得到的翼型壓力分布與Vassberg和Jameson[24]通過4 096×4 096規(guī)模的網(wǎng)格得到的壓力分布的對比，Cp為壓力系數(shù)。由于各計(jì)算結(jié)果收斂較為一致，為清晰起見，圖中只給出了LPTM法與PTM法各一個(gè)結(jié)果，并未逐一對比。從圖中可見，LPTM法與PTM法得到的計(jì)算結(jié)果有較高的一致性，并且與參考結(jié)果符合較好。這表明，本文LPTM法并沒有因?yàn)樘岣哂?jì)算效率而犧牲精度。

圖3 無黏NACA0012翼型算例的表面壓力系數(shù)分布

3.3 湍流RAE2822翼型

本算例考察了LPTM法對二維湍流問題的收斂加速作用。計(jì)算算例采用AGARD-AR-138報(bào)告中湍流RAE2822翼型實(shí)驗(yàn)[25]的Case10。實(shí)驗(yàn)來流馬赫數(shù)0.75，基于弦長的雷諾數(shù)6.2×106，實(shí)驗(yàn)迎角3.19°，升力系數(shù)0.743。計(jì)算時(shí)為匹配實(shí)驗(yàn)升力系數(shù)，將迎角修正為3.09°。計(jì)算網(wǎng)格總量為58 720單元，物面第一層網(wǎng)格高度為1.0×10-5，法向增長率為1.2，平均y+為0.9。網(wǎng)格的大致結(jié)構(gòu)和疏密分布如圖 4所示。

圖5(a)、(b)給出了本算例升力系數(shù)和阻力系數(shù)的收斂過程。比較值得注意的一點(diǎn)是，從收斂歷程的形態(tài)上看，在cLP較小時(shí)，收斂過程似乎只是由于CFL數(shù)較大而被壓縮了。而當(dāng)cLP較大時(shí)，收斂過程的形態(tài)也發(fā)生了變化，這意味著拉普拉斯算子的存在改變了迭代過程在解空間的路徑。圖5(c)給出了本算例系統(tǒng)殘差隨迭代步數(shù)變化的過程。在迭代初始階段，CFL數(shù)相同的LPTM1～4算例，較大cLP的算例殘差收斂速度較快，直接體現(xiàn)了拉普拉斯算子的作用。而整體收斂歷史同樣以CFL數(shù)分為4簇，體現(xiàn)了CFL數(shù)對收斂效率的重要影響，以及通過拉普拉斯延拓增加可用CFL數(shù)的意義。

圖4 湍流RAE2822翼型算例的計(jì)算網(wǎng)格

算 例CFL0c0LPnt/s相對時(shí)間節(jié)約/%PTM139451.00-20.4PTM247842.370LPTM185×10-45536.8613.0LPTM285×10-55231.4325.8LPTM381×10-55633.2521.5LPTM485×10-65732.8022.6LPTM565×10-76536.5213.8

圖5 湍流RAE2822翼型算例升力系數(shù)、阻力系數(shù)及殘差收斂曲線

圖6給出了本算例得到的表面壓力系數(shù)分布。可以看出，計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果符合較好，并且LPTM法與經(jīng)典的PTM法得到的壓力分布有較高的一致性。

3.4 三維ONERA M6機(jī)翼

本算例通過ONERA M6[26]算例來測試LPTM法在三維黏性計(jì)算中的表現(xiàn)。計(jì)算狀態(tài)為來流馬赫數(shù)0.839 5，基于參考弦長Cref=0.646 1 m的雷諾數(shù)為1.172×107，迎角3.06°。相對于前2個(gè)算例，本算例雷諾數(shù)更高，因此物面第一層網(wǎng)格尺度更小，數(shù)值離散后的控制方程穩(wěn)定性問題更為突出。本算例計(jì)算所用湍流模型仍為SA模型。計(jì)算網(wǎng)格5 997 680單元，物面第一層網(wǎng)格法向高度為2.0×10-6，法向增長率為1.2，平均y+為1.2。網(wǎng)格的大致結(jié)構(gòu)和疏密分布如圖7所示。

圖6 湍流RAE2822翼型算例的表面壓力系數(shù)分布

圖7 三維ONERA M6機(jī)翼算例的計(jì)算網(wǎng)格

算 例CFL0c0LPnt/s相對時(shí)間節(jié)約/%PTM13808864.73-12.5PTM24667877.230LPTM1105×10-4487167.509.0LPTM2105×10-5426231.5920.9LPTM3101×10-5436074.6422.9LPTM4105×10-6456154.6521.9LPTM555×10-7587190.528.7

圖8 三維ONERA M6機(jī)翼算例升力系數(shù)、阻力系數(shù)及殘差收斂曲線

圖8(a)、(b)給出了本算例升力系數(shù)和阻力系數(shù)的收斂曲線，可見不同cLP值對力系數(shù)收斂的加速效果。圖8(c)給出了本算例殘差隨迭代步數(shù)的收斂歷史，收斂曲線同樣依CFL數(shù)分為4簇。在迭代初期，引入拉普拉斯延拓的算例同樣顯示了較快的收斂速度。圖9給出了本文計(jì)算得到的表面壓力分布與實(shí)驗(yàn)結(jié)果的對比，y/b為截面展向站位?？梢钥闯觯?jì)算結(jié)果與實(shí)驗(yàn)結(jié)果符合較好，并且LPTM法與經(jīng)典的PTM法得到的壓力分布有較高的一致性。

圖9 三維ONERA M6機(jī)翼算例表面壓力系數(shù)分布

4 結(jié) 論

為提高速定常流動模擬的計(jì)算效率，構(gòu)造了一種在經(jīng)典偽時(shí)間推進(jìn)法基礎(chǔ)上附加一個(gè)拉普拉斯項(xiàng)的混合延拓方法。

1) 在常規(guī)的偽時(shí)間推進(jìn)法中，偽時(shí)間項(xiàng)作為延拓項(xiàng)，通過增強(qiáng)控制方程的線性以及改善雅可比矩陣的對角占優(yōu)實(shí)現(xiàn)了牛頓迭代的穩(wěn)定化。拉普拉斯算子是一個(gè)對稱正定的線性橢圓算子，不僅與偽時(shí)間項(xiàng)同樣有增強(qiáng)方程正則性的作用，還有提供耗散、加速邊界信息向流場內(nèi)部傳播的作用。綜合利用拉普拉斯算子和偽時(shí)間項(xiàng)各自的優(yōu)勢可以構(gòu)造效率更高的計(jì)算方法。

2) 拉普拉斯延拓法對收斂效率的實(shí)際影響與迭代參數(shù)的選擇有密切關(guān)系。本文通過數(shù)值實(shí)驗(yàn)，總結(jié)得到了計(jì)算效率較優(yōu)的一個(gè)拉普拉斯項(xiàng)參數(shù)平臺區(qū)，可以作為在實(shí)際問題中應(yīng)用本文方法的一個(gè)參考。具體地，通過無黏NACA0012翼型、湍流RAE2822翼型和三維ONERA M6機(jī)翼等3個(gè)算例，驗(yàn)證了拉普拉斯延拓方法分別在有黏/無黏、二維/三維對定常流動求解的加速效果。計(jì)算結(jié)果表明，拉普拉斯延拓法可以節(jié)省20%以上的計(jì)算時(shí)間，在迭代過程中引入拉普拉斯算子作為加速收斂措施，是有實(shí)際應(yīng)用價(jià)值的。此外，即使選定相對次優(yōu)的拉普拉斯參數(shù)時(shí)，本文方法仍能起到部分加速作用，而不至于完全失效或起到負(fù)效果，證明本文方法在實(shí)際應(yīng)用中足夠的魯棒性。

參考文獻(xiàn)(References)

[1] GEUZAINE P.Newton-Krylov strategy for compressible turbulent flows on unstructured meshes[J].AIAA Journal,2001,39(3):528-531.

[2] WONG J S,DARMOFAL D L,PERAIRE J.The solution of the compressible Euler equations at low mach numbers using a stabilized finite element algorithm[J].Computer Methods in Applied Mechanics and Engineering,2001,190(43-44):5719-5737.

[3] DENNIS J,SCHNABEL R.Numerical methods for unconstrained optimization and nonlinear equations[M].Philadelphia:Society for Industrial and Applied Mathematics,1996：116-126.

[4] MICHALAK C,OLLIVIER-GOOCH C.Globalized matrix-explicit Newton-GMRES for the high-order accurate solution of the euler equations[J].Computers & Fluids,2010,39(7):1156-1167.

[5] BROWN D A,ZINGG D W.Advances in homotopy continuation methods in computational fluid dynamics:AIAA-2013-2944[R].Reston:AIAA,2013.

[6] LYRA P R M,MORGAN K.A review and comparative study of upwind biased schemes for compressible flow computation.PartⅢ:Multidimensional extension on unstructured grids[J].Archives of Computational Methods in Engineering,2002,9(3):207-256.

[7] KUZMIN D,LOHNER R,TUREK S.Flux-corrected transport:Principles,algorithms,and applications[M].2nd ed.Berlin:Springer,2012:193-238.

[8] KUZMIN D,LOHNER R,TUREK S.Flux-corrected transport:Principles,algorithms,and applications[M].Berlin:Springer,2005:207-250.

[9] COFFEY T S,KELLEY C T,KEYES D E.Pseudotransient continuation and differential-algebraic equations[J].SIAM Journal on Scientific Computing,2003,25(2):553-569.

[10] KELLEY C T,LIAO L Z,QI L,et al.Projected pseudotransient continuation[J].SIAM Journal on Numerical Analysis,2008,46(6):3071-3083.

[11] CEZE M,FIDKOWSKI K J.A robust adaptive solution strategy for high-order implicit CFD solvers:AIAA-2011-3696[R].Reston:AIAA,2011.

[12] CEZE M,FIDKOWSKI K J.Constrained pseudo-transient continuation[J].International Journal for Numerical Methods in Engineering,2015,102(11):1683-1703.

[13] YOUNG D P,MELVIN R G,BIETERMAN M B,et al.Global convergence of inexact Newton methods for transonic flow[J].International Journal for Numerical Methods in Fluids,1990,11(8):1075-1095.

[14] HICKEN J,BUCKLEY H,OSUSKY M,et al.Dissipation-based continuation:A globalization for inexact-Newton solvers:AIAA-2011-3237[R].Reston:AIAA,2011.

[15] HICKEN J,ZINGG D.Globalization strategies for inexact-Newton solvers:AIAA-2009-4139[R].Reston:AIAA,2009.

[16] POLLOCK S.A regularized Newton-like method for nonlinear PDE[J].Numerical Functional Analysis and Optimization,2015,36(11):1493-1511.

[17] KNOLL D A,KEYES D E.Jacobian-free Newton-Krylov methods:A survey of approaches and applications[J].Journal of Computational Physics,2004,193(2):357-397.

[18] BLAZEK J.Computational fluid dynamics:Principles and applications[M].2nd ed. Oxford:Elsevier Science,2005:227-270.

[19] SPALART P R,ALLMARAS S R.A one-equatlon turbulence model for aerodynamic flows:AIAA-1992-0439[R].Reston:AIAA,1992.

[20] POLYANIN A D,NAZAIKINSKII V E.Handbook of linear partial differential equations for engineers and scientists[M].2nd ed.Boca Raton:Chapman and Hall/CRC,2015:1199-1231.

[21] CEZE M,FIDKOWSKI K.Pseudo-transient continuation,solution update methods,and CFL strategies for DG discretizations of the RANS-SA equations:AIAA-2013-2686[R].Reston:AIAA,2013.

[22] MULDER W A,LEER B V.Experiments with implicit upwind methods for the euler equations[J].Journal of Computational Physics,1985,59(2):232-246.

[23] 張涵信.關(guān)于CFD高精度保真的數(shù)值模擬研究[J].空氣動力學(xué)學(xué)報(bào),2016,34(1):1-4.

ZHANG H X.Investigation on fidelity of high order accuate numerical simulation for computational fluid dynamics[J].Acta Aerodynamica Sinica,2016,34(1):1-4 (in Chinese).

[24] VASSBERG J,JAMESON A.In pursuit of grid convergence,Part I:Two-dimensional Euler solutions:AIAA-2009-4114[R].Reston:AIAA,2009.

[25] COOK P H,MCDONALD M A,FIRMIN M C P.Aerofoil RAE2822—Pressure distribution and boundary layer and wake measurements:AGARD-AR-138[R].Reston:AGARD,1979.

[26] SCHMITT V,CHARPIN F.Pressure distributions on the ONERA-M6-wing at transonic mach numbers:AGARD-AR-138[R].Reston:AGARD,1979.