馮選燕，燕振國(guó)，朱華君，馬燕凱，馮新龍

(1.新疆大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院，烏魯木齊 830000；2.空天飛行空氣動(dòng)力科學(xué)與技術(shù)全國(guó)重點(diǎn)實(shí)驗(yàn)室，綿陽(yáng) 621000)

0 引 言

非線性問(wèn)題廣泛存在于各類學(xué)科和實(shí)際問(wèn)題中，是科學(xué)研究中的難點(diǎn)之一[1]。其高效數(shù)值求解方法的研究不僅在計(jì)算數(shù)學(xué)基礎(chǔ)理論上具有重要意義，而且在科學(xué)與工程計(jì)算中處于重要地位，能夠促進(jìn)計(jì)算流體力學(xué)（Computational Fluid Dynamics, CFD）等科學(xué)研究和工程應(yīng)用的發(fā)展。

具體到CFD 領(lǐng)域，其控制方程（Navier-Stokes equations，即N-S 方程組）具有高度非線性特征，對(duì)于絕大部分實(shí)際流動(dòng)問(wèn)題，往往只能對(duì)控制方程進(jìn)行離散求解，為了獲得全離散方程，需要同時(shí)進(jìn)行空間離散和時(shí)間離散。通常的空間離散包括限差分法、有限體積法和間斷有限元（discontinuous Galerkin-finite element method, DG-FEM）等數(shù)值方法[2-4]，時(shí)間離散則通常被劃分為顯式和隱式兩類[5]。其中隱式時(shí)間推進(jìn)方法是解決CFD 剛性流動(dòng)模擬效率問(wèn)題的重要手段[6]，然而隱式時(shí)間離散后不可避免地需要求解非線性離散系統(tǒng)。發(fā)展這類大型非線性離散系統(tǒng)的基礎(chǔ)理論和高效求解方法是提高CFD 計(jì)算效率的關(guān)鍵手段[7]。

Newton 法是求解非線性方程的常用方法之一，其二階的收斂效率在高效CFD 計(jì)算中十分具有吸引力[8]。然而傳統(tǒng)的Newton 法需要計(jì)算非線性殘差方程的Jacobian（雅可比）矩陣來(lái)構(gòu)造線性迭代方程組，并且要求對(duì)線性方程組進(jìn)行精確求解，上述兩點(diǎn)在大規(guī)模CFD 計(jì)算中都是較難實(shí)現(xiàn)的。為了避免上述困難，通常會(huì)在上述兩個(gè)環(huán)節(jié)引入近似。例如，CFD 中常用的Jacobian-free Newton-Krylov（JFNK）方法就采用有限差分來(lái)近似Jacobian 矩陣，引入了Jacobian 矩陣計(jì)算誤差[9-11]；而線化后的方程組系統(tǒng)通常采用迭代方法進(jìn)行求解，迭代方法會(huì)引入線性迭代誤差[12-13]。所以CFD 等大型模擬中實(shí)際上采用的是非精確Newton 方法[14]。

對(duì)于非精確Newton 方法，國(guó)內(nèi)外開(kāi)展了近似方法的系列研究工作。Vanden 和Orkwis[15]討論了精確和數(shù)值Jacobian 矩陣的求解過(guò)程及其對(duì)計(jì)算的影響，數(shù)值測(cè)試結(jié)果表明：在這兩種情況下獲得了較一致的迭代收斂次數(shù)及相近的收斂殘差，并指出對(duì)于較為復(fù)雜的非線性系統(tǒng)，數(shù)值Jacobian 矩陣可能更具有優(yōu)勢(shì)。Ezertas[16]等在此基礎(chǔ)上進(jìn)一步研究了數(shù)值Jacobian 矩陣計(jì)算中誤差的來(lái)源，探究并給出了使得有限差分近似誤差最小化的有限差分?jǐn)_動(dòng)幅值。測(cè)試發(fā)現(xiàn)：在該條件下，采用數(shù)值Jacobian 矩陣表現(xiàn)出了與精確Jacobian 矩陣相當(dāng)?shù)氖諗啃阅堋；诖耍瑪?shù)值Jacobian 矩陣具有實(shí)現(xiàn)簡(jiǎn)單等優(yōu)勢(shì)，被廣泛應(yīng)用于CFD 大規(guī)模數(shù)值模擬中。然而，Yang 等[17]的數(shù)值測(cè)試結(jié)果表明：采用無(wú)矩陣（Jacobian-free, JF）方法的數(shù)值Jacobian 矩陣并不能獲得與精確Jacobian 矩陣相當(dāng)?shù)氖諗刻匦裕邢薏罘终`差的引入會(huì)對(duì)迭代的收斂過(guò)程產(chǎn)生明顯的影響。

另外，非精確Newton 方法在每一步的內(nèi)迭代中通過(guò)線性迭代收斂判據(jù)來(lái)控制求解線性方程組時(shí)解的精確程度，在此過(guò)程中會(huì)產(chǎn)生迭代的收斂誤差。對(duì)于該問(wèn)題，文獻(xiàn)[18]將其稱為Newton 迭代的強(qiáng)制項(xiàng)，其中強(qiáng)制項(xiàng)控制系數(shù) ηk的取值大小對(duì)收斂性有很大的影響。一些學(xué)者[19]結(jié)合理論分析和數(shù)值測(cè)試給出了該參數(shù)的不同選取方式，并給出了強(qiáng)制項(xiàng)能夠保證Newton 方法二階收斂速度的理論分析。然而，在考慮Jacobian 矩陣誤差的情況下，線性迭代收斂判據(jù)的選取方式尚未得到深入研究，需要進(jìn)一步考慮Jacobian矩陣誤差的影響并展開(kāi)對(duì)比研究。

因此，本文重點(diǎn)對(duì)存在Jacobian 矩陣誤差情況下不同線性迭代收斂判據(jù)的表現(xiàn)進(jìn)行了研究，分析了不同線性迭代判據(jù)所引起的過(guò)度求解問(wèn)題。結(jié)合兩種類型的判據(jù)發(fā)展了一種新的線性迭代收斂判據(jù)，并利用數(shù)值測(cè)試對(duì)新提出線性迭代收斂判據(jù)的優(yōu)勢(shì)進(jìn)行了驗(yàn)證。

對(duì)于本文所關(guān)心的CFD 問(wèn)題，其求解的Euler（歐拉）和N-S 方程，可以描述為以下形式：

1 控制方程和數(shù)值方法

1.1 隱式時(shí)間推進(jìn)方法

對(duì)于Euler 或N-S 方程，式（3）是非線性方程組。采用Newton 方法對(duì)該方程進(jìn)行線性化后迭代求解，迭代式為：

其中，J(Q)是 殘 差 向量N(Q)的 Jacobian 矩 陣，ΔQn=Qn+1 ?Qn。

為了增強(qiáng)計(jì)算的穩(wěn)定性，通常會(huì)在式（4）左端添加偽時(shí)間項(xiàng)，此時(shí)不考慮時(shí)間項(xiàng)的計(jì)算精度，時(shí)間離散通常采用實(shí)現(xiàn)簡(jiǎn)單、存儲(chǔ)量小的Euler 隱式格式，此時(shí)迭代方程變?yōu)椋?/p>

式（5）中時(shí)間步長(zhǎng) Δt→∞時(shí)，退化為式（4）的形式。

1.2 Newton-Krylov 方法

Newton-Krylov 方法[20]目前是CFD 領(lǐng)域常用的大型非線性方程組求解算法，基于Krylov 空間，其具體形式為：

該類方法在Krylov 空間內(nèi)構(gòu)造正交的基向量，并在該空間中尋找使得線性方程組誤差最小的數(shù)值解，因此有殘差單調(diào)下降等優(yōu)勢(shì)。其收斂速度依賴于線性方程組系數(shù)矩陣的條件數(shù)（由矩陣的特征值分布決定），當(dāng)矩陣的條件數(shù)較小時(shí)，Krylov 方法可以在幾次搜索內(nèi)將殘差降到很小的量級(jí)。這里為了求解Euler 和N-S 方程產(chǎn)生的非對(duì)稱系統(tǒng)，采用廣義最小殘差法（generalized minimum residual，GMRES）[21-22]，具體如下：

其中：M代表最大內(nèi)迭代次數(shù)， σ為GMRES 迭代的收斂判據(jù)。Jqi的計(jì)算可采用有限差分近似得到，具體形式在本文1.3 節(jié)中詳細(xì)介紹。N代表GMRES 的最大重啟次數(shù)，重啟可以有效避免存儲(chǔ)過(guò)多Krylov 空間向量，減少正交化操作，是大規(guī)模計(jì)算中常用的技術(shù)。

而采用Krylov 方法求解線性方程組時(shí)帶來(lái)的線性迭代誤差，主要依賴于線性迭代收斂判據(jù)：

其 中， ΔQk=Qk+1?Qk； σ可控制線性迭 代 的 求 解 精度，當(dāng)其取值過(guò)大時(shí)，較大的線性迭代誤差可能對(duì)Newton 迭代的收斂造成較大影響，而取值過(guò)小又容易導(dǎo)致過(guò)度求解，造成計(jì)算浪費(fèi)，因此 σ的選擇對(duì)于發(fā)展高效的數(shù)值方法來(lái)說(shuō)是十分重要的。

1.3 線性化方法及其誤差

在每一步Newton 迭代的過(guò)程中，需要采取直接法或迭代法求解由Jacobian 矩陣J組成的線性方程組。在標(biāo)準(zhǔn)的Newton 方法中需要采用精確方法計(jì)算J，然而對(duì)于大規(guī)模CFD 計(jì)算而言，計(jì)算精確的雅可比矩陣所產(chǎn)生的計(jì)算量和存儲(chǔ)量都是較大的。因此，歷史上提出了多種近似雅可比矩陣方法，其中一種是對(duì)雅可比矩陣的表達(dá)式進(jìn)行近似，例如在經(jīng)典的LUSGS 方法中對(duì)流通量和黏性通量的Jacobain 矩陣分別采用一階精度近似和黏性譜半徑近似[5]；另外一種方法則是采用數(shù)值方法計(jì)算Jacobian 矩陣，例如CFD中常用的JF 方法就是采用有限差分方法近似Jacobian 矩陣；此外，還存在上述兩種近似的混合方法[11]。本文主要討論采用JF 方法的誤差及其影響。

在構(gòu)造Krylov 子空間（式6）過(guò)程中，令q=Jj?1r(k)，則Jq的計(jì)算可采用有限差分近似計(jì)算Jacobian 矩陣與向量的乘積：

2 誤差影響研究

本文考慮的非精確Newton 迭代法在求解大型非線性方程組時(shí)相對(duì)精確Newton 方法存在兩類誤差：Jacobian 矩陣誤差和線性迭代誤差。這里首先給出存在兩種誤差情況下非精確Newton 方法的迭代式，其次給出了兩種誤差對(duì)Newton 迭代式的影響，最后分析了Jacobian 矩陣誤差來(lái)源驗(yàn)證及其影響。

2.1 誤差及其分析

非精確Newton 迭代的過(guò)程中，Jacobian 矩陣誤差和線性迭代誤差都對(duì)整個(gè)非線性方程組求解過(guò)程中的收斂性造成一定影響。對(duì)于經(jīng)典Newton 法，解的迭代式為：

作為示例，這里給出Euler 方程采用DG 方法離散后的具體離散形式和Jacobian 矩陣的形式，詳細(xì)矩陣系數(shù)和N-S 方程的離散形式參照文獻(xiàn)[11]。具體殘差形式為：

這里d為空間維數(shù)，M=BTΛ(wJ)B為 質(zhì)量矩陣， Λ為對(duì)角矩陣，D j為第j個(gè)方向的微分矩陣，B為從求解點(diǎn)或者模態(tài)系數(shù)到通量點(diǎn)的變換矩陣、BΓ為到單元邊界上的通量點(diǎn)的變換矩陣，w和J分別代表數(shù)值積分權(quán)重和網(wǎng)格變換Jacobian，Mc為 單元邊界基函數(shù) ?Γ和基函數(shù) ?之間的映射矩陣，F(xiàn)^n為對(duì)流數(shù)值通量。

其中，空間殘差項(xiàng)N(Q) 關(guān) 于守恒變量Q求導(dǎo)的Jacobian 矩陣J(Qn)形式：

?L(Qn)/?Q具 有 塊 結(jié) 構(gòu)，例 如 ?Le1(Qn)/?Qe2代 表 第e1行和第e2列的塊Jacobian 矩陣元素，其含義為第e1個(gè)單元的 L 對(duì)第e2個(gè)單元的Q的Jacobian 矩陣塊。其具體形式為：

2.2 Jacobian 矩陣誤差來(lái)源驗(yàn)證及其影響

本小節(jié)首先簡(jiǎn)單驗(yàn)證了采用有限差分近似的Jacobian 矩陣的相對(duì)誤差變化規(guī)律，驗(yàn)證了Jacobian矩陣誤差主要來(lái)源于差分格式帶來(lái)的截?cái)嗾`差和計(jì)算機(jī)的舍入誤差[16]，而且表明在特定的數(shù)據(jù)精度和計(jì)算格式下，都存在一個(gè)使得誤差最小的ε，即最優(yōu)ε，如圖1 所示。圖中FD 表示一階前向差分，CD 表示二階中心差分。

圖1 數(shù)值Jacobian 矩陣的相對(duì)誤差變化規(guī)律Fig.1 Relative error variations of the numerical Jacobian matrices

為進(jìn)一步測(cè)試Jacobian 矩陣誤差對(duì)收斂速度的影響，這里利用周期計(jì)算域內(nèi)的一維Euler 系統(tǒng)所搭建的簡(jiǎn)單模型問(wèn)題進(jìn)行研究，此處選取120 個(gè)網(wǎng)格點(diǎn)，通過(guò)控制網(wǎng)格拉伸比使得網(wǎng)格向計(jì)算域兩側(cè)匯聚，并通過(guò)控制計(jì)算網(wǎng)格的大小比值來(lái)控制問(wèn)題的剛性，此處網(wǎng)格拉伸比為1.05?？臻g離散采用一階迎風(fēng)格式，該模型問(wèn)題Jacobian 矩陣的條件數(shù)約為1 ×108。圖2 中，橫坐標(biāo)（Non-Iter）代表非線性迭代次數(shù)，結(jié)果驗(yàn)證了當(dāng)數(shù)值Jacobian 矩陣近似誤差較大時(shí)會(huì)引起效率降低甚至無(wú)法收斂的問(wèn)題[16-17]，然而當(dāng)Jacobian矩陣誤差較小時(shí)則可能避免上述問(wèn)題。

圖2 不同量級(jí)Jacobian 矩陣誤差對(duì)非精確Newton 迭代的影響Fig.2 Influence of the Jacobian matrix error order on the inexact Newton iteration

3 線性迭代收斂判據(jù)研究

對(duì)于Newton 類非線性系統(tǒng)求解方法，線性迭代誤差雖然不影響計(jì)算最終收斂后的精度，但是對(duì)非線性系統(tǒng)能否收斂以及收斂時(shí)的計(jì)算效率有較大的影響。其中，相關(guān)研究[18]給出的線性系統(tǒng)迭代收斂判據(jù)的選取準(zhǔn)則要求，在保證Newton 方法二階收斂效率的前提下盡量避免線性系統(tǒng)的過(guò)度求解。所謂過(guò)度求解就是當(dāng)線性系統(tǒng)的求解精度達(dá)到一定程度以后，繼續(xù)提升線性系統(tǒng)的求解精度將無(wú)法顯著提升非線性系統(tǒng)的收斂效率，此時(shí)額外增加的線性系統(tǒng)迭代步數(shù)就造成了線性系統(tǒng)的過(guò)度求解。從提高計(jì)算效率角度，避免過(guò)度求解顯然是必須考慮的問(wèn)題。

常用的迭代誤差收斂判據(jù)有兩種，一種是絕對(duì)量型收斂判據(jù)，即選取某個(gè)絕對(duì)的數(shù)值作為迭代收斂的判據(jù)，為了使得該策略中的參數(shù)更具有普適性，本文采用如下的形式：

另外一種是相對(duì)型收斂判據(jù)，即將線性系統(tǒng)的迭代收斂判據(jù)與非線性系統(tǒng)的殘差進(jìn)行關(guān)聯(lián)，其形式如下：

此處 ηk為線性系統(tǒng)收斂判據(jù)的主要控制參數(shù)，其含義為線性系統(tǒng)迭代誤差要相對(duì)當(dāng)前的非線性系統(tǒng)殘差小 ηk倍以上。文獻(xiàn)[19]對(duì)式（20）的線性迭代收斂判據(jù)進(jìn)行了較系統(tǒng)的研究，并從能夠保持Newton 方法二階收斂效率的角度給出了相應(yīng)的數(shù)學(xué)證明和數(shù)值測(cè)試；然而并未對(duì)式（19）的絕對(duì)量型收斂判據(jù)和存在Jacobian 矩陣誤差情況下式（20）的迭代收斂誤差的影響開(kāi)展研究。

3.1 不同 ηk的對(duì)比研究

這里采用風(fēng)雷高精度軟件平臺(tái)[24]開(kāi)展測(cè)試研究?？刂品匠滩捎肗-S 方程，算例為黏性NACA 0012翼型流動(dòng)問(wèn)題，來(lái)流Ma=0.1，雷諾數(shù)取100，攻角為0。，翼型表面采用無(wú)滑移邊界條件，計(jì)算區(qū)域采用20 倍翼型長(zhǎng)度，在外邊界采用遠(yuǎn)場(chǎng)邊界條件。采用2842個(gè)四邊形網(wǎng)格對(duì)計(jì)算域進(jìn)行離散，如圖3 所示。在當(dāng)前的定常計(jì)算中，添加了偽時(shí)間項(xiàng)來(lái)控制線性系統(tǒng)的剛性以增強(qiáng)計(jì)算的穩(wěn)定性，計(jì)算結(jié)果如圖4 所示。本計(jì)算主要關(guān)注殘差的收斂過(guò)程。

圖3 NACA 0012 算例的計(jì)算網(wǎng)格分布Fig.3 Computational grid distribution for the case of NACA 0012

圖4 NACA 0012 算例計(jì)算結(jié)果密度分布Fig.4 Density distribution for the flow around NACA 0012

如圖5 所示，首先對(duì)比測(cè)試了 ηk分別取1×10?3、1×10?4和1×10?6情況下采用式（8）的數(shù)值Jacobian 矩陣（JF1）的迭代收斂情況，作為對(duì)比，本文同時(shí)給出了采用精確Jacobian 的結(jié)果（JF0）。從線性迭代步數(shù)（Iter）上來(lái)講，顯然采用更小的 ηk所需的線性迭代步數(shù)顯著增加，線性方程組的求解精度相應(yīng)提升。然而除了 ηk=1×10?3以外，所有非線性系統(tǒng)殘差基本完全重疊；當(dāng) ηk=1×10?3時(shí)，收斂曲線在收斂初段和最終階段已經(jīng)出現(xiàn)了較明顯的偏差。上述結(jié)果表明，當(dāng)ηk≤1 ×10?4時(shí)線性迭代已經(jīng)足夠精確，繼續(xù)降低 ηk則會(huì)造成過(guò)度求解問(wèn)題，使計(jì)算效率降低。

圖5 NACA 0012 算例在不同Jacobian 矩陣下非線性迭代殘差和GMRES 迭代步數(shù)對(duì)比Fig.5 Comparison of the nonlinear iterative residual and GMRES iteration steps for different Jacobian matrices in the case of NACA 0012

另外需要強(qiáng)調(diào)的是，在迭代接近收斂時(shí)，基于ηk的迭代收斂判據(jù)都會(huì)造成過(guò)度求解問(wèn)題。本文進(jìn)一步通過(guò)增加GMRES 的迭代重啟次數(shù)來(lái)增加允許的最大線性系統(tǒng)迭代總步數(shù)，從而進(jìn)一步研究造成上述問(wèn)題的原因，結(jié)果如圖6 所示。在迭代的最終收斂階段，無(wú)論將允許的最大線性系統(tǒng)迭代步數(shù)增大到多大（在本文的測(cè)試范圍內(nèi)），線性方程組均無(wú)法收斂到 ηk=1×10?6，相應(yīng)的線性方程迭代的總步數(shù)達(dá)到允許的最大值。與之相反，采用精確Jacobian 矩陣的計(jì)算結(jié)果，雖然在收斂的最終階段迭代步數(shù)同樣增加，但是并不存在上述現(xiàn)象，說(shuō)明上述問(wèn)題是當(dāng)所要求的線性系統(tǒng)迭代誤差太小時(shí)，Jacobian 矩陣誤差會(huì)造成線性系統(tǒng)無(wú)法收斂引起的。然而從非線性系統(tǒng)的殘差角度來(lái)看，各條曲線是基本重合的，進(jìn)一步說(shuō)明采用近似Jacobian 所增加的線性系統(tǒng)迭代是由于收斂判據(jù)過(guò)嚴(yán)引起的過(guò)度求解，需要發(fā)展更加合理的迭代收斂判據(jù)避免上述問(wèn)題。

圖6 不同GMRES 允許的最大迭代步數(shù)下的線性迭代次數(shù)Fig.6 Linear iterative number under different allowable maximum GMRES iteration steps

3.2 不同 ζk的對(duì)比研究

類似3.1 節(jié)，這里對(duì)式（19）的收斂判據(jù)進(jìn)行測(cè)試，結(jié)果如圖7 所示。本文分別給出了 ζk取1×10?6、1×10?10、1×10?14的測(cè)試結(jié)果，為了對(duì)比，圖中同時(shí)給出了 ηk=1×10?6情況下的收斂曲線?？梢钥吹剑捎忙苉類 型的迭代收斂判據(jù)，只有在 ζk取很小的值時(shí)才能保證迭代收斂。但是，當(dāng) ζk取值很小時(shí)（如1×10?14），在迭代的初始階段會(huì)造成較嚴(yán)重的過(guò)度求解現(xiàn)象。而在迭代接近收斂的階段，采用 ζk類型的迭代收斂判據(jù)卻可以較好地避免過(guò)度求解的問(wèn)題。

圖7 NACA 0012 算例中不同ζ 取值情況下非線性迭代殘差和GMRES 迭代步數(shù)對(duì)比Fig.7 Comparison of the nonlinear iterative residual and GMRES iteration steps for different ζ in the case of NACA 0012

3.3 一種新的線性迭代收斂判據(jù)

基于以上研究結(jié)果，本文發(fā)展了一種新的迭代收斂判據(jù)：

其中 max函數(shù)內(nèi)兩項(xiàng)分別為式（19）和式（20）對(duì)應(yīng)的絕對(duì)收斂判據(jù)和相對(duì)收斂判據(jù)，這里取兩項(xiàng)中的較大值。 在非線性迭代初期||N(Qk)||2較大，因此ηk||N(Qk)||2的量級(jí)會(huì)大于ζk||Q0||2，從而避免了后者在非線性迭代初期所引起的嚴(yán)重過(guò)度求解問(wèn)題；隨著迭代逐漸接近收斂， ||N(Qk)||2逐漸減小，只要取合理的參數(shù)，可以使得 max函 數(shù)在接近收斂時(shí)取 ζk||Q0||2，此時(shí)可以有效避免圖5 中的過(guò)度求解問(wèn)題。在此基礎(chǔ)上，本文還在 max函 數(shù)外添加了相應(yīng)的 min函數(shù)，其作用是保證每個(gè)線性迭代至少要下降一定量級(jí)（本文中取ηmin=0.01），從而避免類似圖7 中迭代不能完全收斂的問(wèn)題。本文中 ζk取1×10?13，主要從計(jì)算的舍入誤差角度進(jìn)行選取，由于最終線性迭代誤差會(huì)加到Qk上， ζk=1×10?13的參數(shù)選擇可以避免線性方程組收斂到比Newton 迭代步本身舍入誤差更小的量級(jí)，從而避免類似圖5 中的過(guò)度求解問(wèn)題。

圖8 同樣采用NACA 0012 算例對(duì)不同迭代收斂判據(jù)進(jìn)行了對(duì)比，直至迭代第19 步采用式（21）的新迭代收斂判據(jù)（NewTol）得到了與 ηk=1×10?6完全相同的計(jì)算結(jié)果。此后采用 ηk=1×10?6的線性迭代步數(shù)迅速增加到允許的最大值，線性迭代無(wú)法收斂，而采用式（21）的線性迭代步數(shù)則開(kāi)始下降，說(shuō)明如設(shè)計(jì)的那樣，式（21）中 ζk||Q0||2部分開(kāi)始發(fā)揮作用從而大幅減少了線性迭代的計(jì)算量。直至第23 步（殘差收斂至1×10?13）時(shí)采用新判據(jù)的非線性系統(tǒng)殘差依然與ηk=1×10?6重疊，隨后采用式（21）的線性迭代步數(shù)開(kāi)始顯著減小，非線性殘差曲線不再重疊。表1 給出了殘差收斂到不同量級(jí)時(shí)所需要的線性迭代總步數(shù)，該指標(biāo)能夠很好地反映出所消耗的總計(jì)算量。當(dāng)殘差收斂到1×10?10以下時(shí)，JF0- ηk=1×10?6（精確Jacobian、ηk=1×10?6）、JF1- ηk=1×10?6（近似Jacobian、 ηk=1×10?6）和JF1-NewTol（近似Jacobian、新判據(jù)）所需線性迭代步數(shù)分別為13113、19670 和12692 步；當(dāng)殘差收斂到1×10?13以下時(shí)，所需線性迭代步數(shù)分別為16964、29670和13870 步；當(dāng)殘差收斂到1×10?16以下時(shí)，所需線性迭代步數(shù)分別為21473、44670 和15453 步?？梢钥吹?，新發(fā)展的線性迭代收斂判據(jù)能夠有效避免過(guò)度求解的問(wèn)題，顯著降低所需要的線性系統(tǒng)迭代步數(shù)。線性迭代步數(shù)最多降低3 倍，從而提高計(jì)算效率。

表1 不同判據(jù)情況下迭代收斂所需線性迭代步數(shù)對(duì)比Table 1 Comparison of the linear iteration steps under different convergence criteria

圖8 新線性迭代收斂判據(jù)的殘差和迭代步數(shù)Fig.8 Nonlinear iterative residual and GMRES iteration steps for the new linear iteration convergence criterion

為了進(jìn)一步驗(yàn)證上述迭代收斂策略在不同計(jì)算狀態(tài)中的通用性，本文通過(guò)改變離散格式的計(jì)算精度來(lái)改變非線性系統(tǒng)的特性和規(guī)模。將DG 格式的計(jì)算精度提高到三階，此時(shí)總共求解的網(wǎng)格自由度變?yōu)?7052 個(gè)。圖9 給出了三階精度的計(jì)算結(jié)果，可以看到，采用新的線性迭代收斂判據(jù)同樣能夠很好地避免接近收斂狀態(tài)下的過(guò)度求解問(wèn)題，從而顯著提高計(jì)算效率。

圖9 三階DG 格式計(jì)算中新線性迭代收斂判據(jù)的殘差和迭代步數(shù)Fig.9 Nonlinear iterative residual and GMRES iteration steps for the new linear iteration convergence criterion with third-order DG scheme

4 結(jié) 論

針對(duì)CFD 中常見(jiàn)的基于Newton 類方法的非線性系統(tǒng)求解問(wèn)題，本文通過(guò)數(shù)值測(cè)試，對(duì)存在Jacobian矩陣誤差情況下每個(gè)Newton 迭代步中線性迭代收斂判據(jù)進(jìn)行了研究，最終結(jié)合兩種類型的迭代收斂判據(jù)，發(fā)展了一種新的判據(jù)。得到的主要結(jié)論如下：

1）存在Jacobian 矩陣誤差情況下，線性迭代收斂判據(jù)的影響會(huì)更加顯著；

2）基于絕對(duì)量的線性迭代收斂判據(jù)容易在非線性迭代初期引起過(guò)度求解問(wèn)題；相對(duì)非線性殘差的線性迭代收斂判據(jù)容易在接近最終收斂時(shí)引起過(guò)度求解問(wèn)題；

3）新型線性迭代收斂判據(jù)能夠有效地降低過(guò)度求解問(wèn)題，顯著降低收斂所需要的線性迭代步數(shù)，從而提高計(jì)算效率。

下一步，我們將進(jìn)一步對(duì)Jacobian 矩陣誤差對(duì)非線性系統(tǒng)迭代收斂效率的影響開(kāi)展研究，并將發(fā)展的線性迭代收斂誤差應(yīng)用到更加復(fù)雜的流動(dòng)模擬中。