一種跨聲速定常流場求解加速方法

喬磊， 白俊強， 邱亞松,*， 華俊,2， 張揚

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

在飛行器和風力機械氣動特性的評估和優化設計中，需要進行大量的定常流動模擬。在非定常數值模擬中，需要用定常流動的解作為初始解。因此，提高定常流場數值模擬的效率，可以使計算流體力學方法在氣動外形設計應用中發揮更有力的作用。目前，常用的定常流動問題的控制方程是定常可壓縮雷諾平均納維-斯托克斯(Reynolds Averaged Navier-Stokes, RANS)方程，此方程是一個由對流項主導的非線性偏微分方程組，其求解分為顯式和隱式2種方法。由于定常問題不需要考慮時間精度，隱式解法得以充分發揮其穩定性強、計算效率高的優點，因而得到了廣泛的應用。

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

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

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

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

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

本文首先介紹了控制方程和牛頓迭代法等工作基礎，然后詳細闡述了拉普拉斯增穩的迭代思路，分析其與經典偽時間推進法的優缺點，并給出了LTPM法的實現細節，最后，通過無黏NACA0012翼型、湍流RAE2822翼型和三維ONERA M6機翼3個算例，對LPTM法的收斂加速效果進行了驗證。

1 控制方程及其隱式解法

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

1.1 控制方程

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

(1)

式中：R(w)為殘差，流體狀態向量w=(ρ,U,p)，ρ、U和p分別為密度、速度和壓力的流場自變量；F(w)和Fv(w)分別為無黏和黏性通量；Ω為求解域或控制體；v為積分變量。控制方程中的湍流黏性項需要通過求解湍流模型得到。在本文的討論中，湍流模型的求解是和流動控制方程獨立的。因此，本文不涉及關于湍流模型求解的問題。

1.2 偽時間推進法的基本作用

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

(2)

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

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

(3)

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

(4)

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

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

(5)

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

2 拉普拉斯延拓

2.1 拉普拉斯算子的引入

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

(6)

式中：RLPT(w)為由流體狀態向量w計算得到的殘差；cLP為拉普拉斯項的縮放系數；Δ為拉普拉斯算子。這樣，牛頓迭代公式變為

(7)

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

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

2.2 拉普拉斯延拓作用

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

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

(8)

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

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

2.3 迭代策略

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

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

(9)

(10)

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

3 算例驗證

3.1 計算效率的比較方法

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

(11)

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

3.2 無黏NACA0012翼型

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

圖1 無黏NACA0012翼型算例的計算網格

算 例CFL0c0LPnt/s相對時間節約/%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翼型算例升力系數、阻力系數及殘差收斂曲線

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

圖3給出了本文計算得到的翼型壓力分布與Vassberg和Jameson[24]通過4 096×4 096規模的網格得到的壓力分布的對比，Cp為壓力系數。由于各計算結果收斂較為一致，為清晰起見，圖中只給出了LPTM法與PTM法各一個結果，并未逐一對比。從圖中可見，LPTM法與PTM法得到的計算結果有較高的一致性，并且與參考結果符合較好。這表明，本文LPTM法并沒有因為提高計算效率而犧牲精度。

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

3.3 湍流RAE2822翼型

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

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

圖4 湍流RAE2822翼型算例的計算網格

算 例CFL0c0LPnt/s相對時間節約/%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翼型算例升力系數、阻力系數及殘差收斂曲線

圖6給出了本算例得到的表面壓力系數分布。可以看出，計算結果與實驗結果符合較好，并且LPTM法與經典的PTM法得到的壓力分布有較高的一致性。

3.4 三維ONERA M6機翼

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

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

圖7 三維ONERA M6機翼算例的計算網格

算 例CFL0c0LPnt/s相對時間節約/%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機翼算例升力系數、阻力系數及殘差收斂曲線

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

圖9 三維ONERA M6機翼算例表面壓力系數分布

4 結 論

為提高速定常流動模擬的計算效率，構造了一種在經典偽時間推進法基礎上附加一個拉普拉斯項的混合延拓方法。

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

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

參考文獻(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] 張涵信.關于CFD高精度保真的數值模擬研究[J].空氣動力學學報,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.