李孝偉,王正裕

(上海大學,上海市應(yīng)用數(shù)學和力學研究所,上海200072)

0 引 言

計算流體力學的發(fā)展過程是一個流體動力學模型的精確化和細致化的過程,同時流場數(shù)值模擬的收斂性和計算效率也成為了計算流體力學中的關(guān)鍵問題。為了提高計算效率,廣大計算流體力學工作者發(fā)展了多種加速收斂的措施,典型的如采用當?shù)貢r間步長和隱式殘值光順技術(shù)等。但人們并不滿足,因為現(xiàn)實期待一些更為有效的加速收斂方法的出現(xiàn)。

1993年,G.M.Shroff與H.B.Keller[1]提出一種遞歸投影算法(the Recursive Projection Method,RPM),目的是為了穩(wěn)定和加速非線性分叉和混沌的數(shù)值研究中的迭代過程。對于離散的偏微分方程,該方法對迭代變量的殘差構(gòu)造Krylov子空間,并從中分離出影響迭代收斂的強勢特征子空間,進而通過使用牛頓迭代法消除其對迭代收斂進程帶來的負面影響。

1999年,Philip Love[2]在三維的Kolmogorov流動和泰勒渦流的分叉研究中應(yīng)用了 RPM方法來加快定位流場的分叉點,以探討湍流的起源過程。在該研究中,作者證明了RPM方法可以有效地幫助加快發(fā)現(xiàn)低雷諾數(shù)下 Kolmogorov流動的分叉現(xiàn)象。2001年,為了提高RPM方法對于不變子空間的近似精度,V.Janovsky和Liberda.Praha[3]提出一種改進的RPM方法。該方法使用Cayley變形處理雅可比矩陣,以修正超出穩(wěn)定邊界特征值所對應(yīng)的子空間的正交基；在得到的修正的正交基上,提出一種耦合的雙重迭代的 RPM方法,最后證明該修正的RPM方法的能夠發(fā)現(xiàn)一些原始RPM方法發(fā)現(xiàn)不了的Hopf分叉點。2005年,M.Dorobantu和J.M?ller將RPM用于計算流體力學的NSMB程序塊中作為加速數(shù)值收斂進程的外殼[4-7],計算了二維超音速噴管流動和二維亞音速翼型繞流,并證明RPM方法在二維定常流場的計算中可以達到約5倍的加速效果。

從以上研究者們的結(jié)果可以看出,RPM方法的加速收斂效果和作用非常顯著,這激發(fā)了本文作者對它的研究興趣。在研究中發(fā)現(xiàn),原有RPM方法在實施過程中不易保證正交基的精度,針對此,本文提出了一種改型的RPM方法。計算實踐表明,該改型的RPM方法,不僅可以更好地保證正交基的精度,還可以節(jié)省相關(guān)的計算時間。文中一方面還研究了幾個重要參數(shù)對RPM方法加速收斂性能的影響規(guī)律,另一方面還研究了在流場數(shù)值模擬中采用了常規(guī)的加速收斂措施以后改型的RPM方法的加速收斂性能。

1 RPM方法

對于存在定常解的偏微分方程(組),數(shù)值求解時一般是通過空間離散,構(gòu)造一個包含時間導(dǎo)數(shù)項的常微分方程(組)來進行求解,可以表示為

在時間推進上,可以使用一個迭代格式來表述方程(1)的求解過程,

令Fy為由離散的方程組所定義的雅克比矩陣。假設(shè)對于迭代格式(2)式,存在解y*=F(y*),則雅克比矩陣可表達為

從而,(2)式的收斂特性就由雅克比矩陣J的特征值的譜分布來決定。如果該矩陣的所有的特征值的模都小于1,那么(2)式將漸進地收斂到y(tǒng)*；當J有一小部分特征值的模小于1但非常地接近于1時,迭代式(2)的收斂進程會變得非常的緩慢；如果J存在任意一個特征值的模大于1,那么(2)式的迭代格式就會發(fā)散。對于模接近1或大于1的J的特征值,由于它們是影響離散方程(2)式收斂的強勢因子,所以稱之為J的強勢特征值。

為了穩(wěn)定形如(2)式的收斂緩慢或不穩(wěn)定的迭代格式,G.M.Shroff和H.B.Keller[1]1993年提出了RPM方法。RPM方法的主要思想是,通過分離出雅可比矩陣J的強勢特征值所決定的不變子空間,進而使用牛頓迭代法在該子空間中消除強勢特征值對收斂進程所帶來的負面影響。

假設(shè)式(2)的雅可比矩陣J存在小部分的特征值位于區(qū)域Kδ={(|z|≤1-δ} (δ＞0)之外,即有如下形式,|λ1|≥|λ2|≥… ≥|λp|＞1-δ≥|λp+1|≥…≥|λN|,(其中p?N)。對應(yīng)于以上特征值的分布,定義子空間P和Q如下:

P≡對應(yīng)于特征值{λ1,…,λp}的 F*y的最大不變子空間；

Q≡RN-P,為P的正交補。

對于所定義的子空間,通過對特征值{λ1,…,λp}的分離和提取,可以構(gòu)造一組對應(yīng)的不變子空間P的正交基Vp,進而在Vp的基礎(chǔ)上構(gòu)造正交投影算子P和Q。

通過向兩個正交子空間的投影,方程的解也就分解成為了如下兩部分,

其中,對應(yīng)于構(gòu)造的不變子空間的正交基Vp,投影算子可表達為,

將(5)式應(yīng)用到(4)式有

由于不變子空間P是分離出來的影響收斂進程的強勢特征子空間,對應(yīng)于該空間的投影分量pn是收斂緩慢或發(fā)散的解分量。為了對其收斂特性做相應(yīng)的改善,把牛頓方法

應(yīng)用到(6a)式,從而得到如下的RPM 方法的迭代形式,

若設(shè)Δqn=qn-q*,將解 qn+1在 y*=p*+q*處泰勒展開可得,

可以證明[4],上式中g(shù)*q=QJ*Q的所有特征值位于Kδ的鄰域內(nèi),并且QJ*P=0。因此,在以上分解的迭代關(guān)系中,RPM方法首先通過遞歸地對投影分量qn的殘值Δqn進行監(jiān)控而得到補空間的迭代特性變化,并從中構(gòu)造和抽取Jn的強勢特征值所對應(yīng)的不變子空間的正交基Vp,進而相應(yīng)地對迭代格式進行投影改進。這也是“遞歸”和“投影”(Recursive&Projection)所表述的意義。

對于(8a)式中出現(xiàn)的Jn,我們并不需要顯式地直接求取,在所構(gòu)造的正交基Vp的基礎(chǔ)上,通過對(3)式構(gòu)造差分格式,可以直接得到

其中,vj為Vp的列向量；ε為大于0的小量,一般可以取為10-8～10-12。通過差分格式非顯式地求取離散方程的雅可比矩陣也是RPM方法的能夠?qū)崿F(xiàn)以小的計算代價而獲得良好的加速效果的重要原因之一。

至于求取投影算子P和Q的遞歸算法,包括基于累積方式的正交基的計算、正交基的精度的維持等,請參見文獻[1],在此不再重述。下面僅敘述一下RPM 方法對“正交基維數(shù)的減小”的處理

RPM方法在遞歸求取不變子空間的正交基的過程中,檢測的是n時刻的雅可比矩陣的變化。隨著計算進程的發(fā)展,Jn是與n時刻計算所得到的流場變量yn相關(guān)的,因此Jn的強勢特征值也應(yīng)該是隨著計算進程不斷變化的。所以,前面所積累的不變子空間P也就可能無法近似為當前時刻的雅可比矩陣的優(yōu)勢特征子空間,其正交基會變得不準確。這進而會導(dǎo)致整個計算的發(fā)散和失敗。對此,G.M.Shroff和H.B.Keller在文獻[1]中提出需要在適當?shù)臅r刻減小不變子空間的基,其思想如下。

基于已有的強勢特征子空間的正交基Vp,有

可以證明矩陣H-的特征值為的特征值的一個子集,并為不變子空間P部分對應(yīng)的特征值。理論上RPM方法只需要對大小為m×m的矩陣的特征值進行重選,淘汰置于Kδ區(qū)域內(nèi)的特征值,并求得剩余特征值所對應(yīng)的基矢量就可以了。但是我們在具體實施的過程中發(fā)現(xiàn),對矩陣的特征值進行重選不難,但是對重選后得到的復(fù)數(shù)特征值對所對應(yīng)的不變子空間的正交基的構(gòu)造存在困難。對于重選后的矩陣的復(fù)數(shù)特征值對,其所對應(yīng)的特征向量也是復(fù)數(shù)形式。在此,問題的困難在于如何結(jié)合已有的正交基組和復(fù)數(shù)特征向量組構(gòu)造新的正交基組。也就是說,我們在將一個酉空間的基矢量轉(zhuǎn)化成單個的實空間的基矢量時找不到相關(guān)的理論依據(jù)。H.B.Keller提出,將矩陣的復(fù)特征向量分解為由復(fù)數(shù)的實部和虛部分別構(gòu)造的兩個實向量來處理。但是,對于這種處理方法,在計算時如果淘汰的特征值少于復(fù)特征值的數(shù)量時,我們原意是要減小和篩選原不變子空間的正交基,就有可能變?yōu)樵黾釉蛔冏涌臻g的正交基。在數(shù)值試驗中,我們發(fā)現(xiàn),這種處理方法很難保證正交基V p的計算精度和合理性,并會導(dǎo)致計算的失敗。為此,我們發(fā)展了一種更容易保持正交基V p精度的改型的RPM方法。

2 改型的RPM方法

在原RPM過程中,Vp的精度難以保持的最根本原因在于這種累積計算正交基的方式會使得早被計入的一些正交基過時,因而積累的正交基不能實時地反映當前雅可比矩陣Jn的變化。如果要保持正交基Vp的精度,就需要在某些時候減小正交基Vp的維數(shù),也就是需要剔除掉一些過時的正交基。但是,這樣的處理方式又會使算法變得非常復(fù)雜而不易于實現(xiàn)。為此,我們在流場計算的過程中考慮放棄這種累積計算正交基的Vp方式,而采用即時覆蓋的方式計算正交基Vp。

鑒于對所采用的計算正交基Vp的方式的改變,我們對解的投影分量qn的收斂速度的監(jiān)視也應(yīng)該相應(yīng)地改變?yōu)閷獾恼w收斂歷程的監(jiān)視:通過監(jiān)視Δyn+1=yn+1-yn的變化,實時地求取全局的收斂緩慢的強勢特征子空間的基,以覆蓋原有的基Vp,即在每次RPM迭代之后,我們得到新的yn+1,計算

并將得到的Δyn+1累積到Krylov向量組中,

對Krylov向量組采用QR分解并提取滿足Ka接受比條件的基Vp,覆蓋原有的正交基Vp,進而進入新一輪的RPM數(shù)值迭代。

本文在采用覆蓋的方式計算正交基Vp時,首先使用已有的正交基Vp,按(8)式對迭代進行加速計算,當累積了 ks組的 Δyn后,對由 Δyn構(gòu)造的 Krylov向量組求取新的正交基Vp,并覆蓋已有的正交基Vp,再重復(fù)(8)式的數(shù)值迭代。這種以覆蓋方式計算正交基Vp的RPM方法,不僅可以更好地保證正交基Vp的精度,還可以節(jié)省原RPM方法中維持基的精度和減小基的維數(shù)的相關(guān)計算。

本文后面部分提到的RPM方法均指的是該改型的RPM方法。

3 幾個參數(shù)的影響

給出方程:

可以看出,以上方程的解析解為y*=1且與a的取值無關(guān)。對方程(14)構(gòu)造迭代格式

通過前面的描述可以知道,當存在某個a(i)的絕對值趨近于1時,它所對應(yīng)的y(i)值會很緩慢地收斂于y=1。取空間離散點數(shù)N=31,y0(i)=5.0,a(i)=0.1+0.89×i/N,i=1,2,…,N?？梢钥闯?有7個特征值大于0.80,有4個特征值大于0.90。

加入 RPM方法,程序中設(shè)定以下參數(shù):Ka:Krylov接受比；ε:計算 Jnvj時的小量 ；Nstep:使用迭代格式,每計算 Nstep步y(tǒng) 值,計算一次 Δyi；Ks:Krylov向量組的大小。

分別檢測了以上參數(shù)對RPM數(shù)值迭代效果的影響。收斂標準統(tǒng)一采用Ry＜10-10。計算結(jié)果見表1～表4。從表中可以看出:(1)Ka的取值大小對RPM的收斂效果存在明顯的影響:Ka越大,發(fā)現(xiàn)和提取第一組不變子空間正交基的時間越遲,該組正交基的維數(shù)也相對偏小；(2)ε對發(fā)現(xiàn)正交基的時間以及正交基的維數(shù)沒有顯著的影響；(3)在具有靜態(tài)的特征值數(shù)值試驗中,Nstep的取值極大地影響著RPM的加速收斂效果,高頻率地求取Δyi可以相對快速地發(fā)現(xiàn)強勢特征值所對應(yīng)的不變子空間,從而以相對較小的總迭代次數(shù)達到收斂標準；(4)Ks對RPM迭代的收斂效果不具有明顯的影響,對強勢特征值所對應(yīng)的不變子空間的正交基的添加數(shù)目也沒有明顯的影響。

以上結(jié)果說明,對于收斂緩慢的數(shù)值迭代格式,RPM方法具有非常明顯的加速收斂效果；而對RPM方法的效果影響較明顯的參數(shù)有Ka和Nstep。在所構(gòu)造的簡單模型中,RPM方法的加速最佳效果達到70倍。

表1 Ka對RPM數(shù)值過程的影響Table1 Effect of Ka on the process of RPM

表2 ε對RPM數(shù)值過程的影響Table 2 Effect ofεon the process of RPM

表3 Nstep對RPM數(shù)值過程的影響Table 3 Effect of Nstep on the process of RPM

表4 Ks對RPM數(shù)值過程的影響Table 4 Effect of Ks on the process of RPM

4 RPM方法流場數(shù)值模擬中的加速收斂效果

本節(jié)數(shù)值模擬了二維翼型繞流。流場控制方程采用Navier-Stokes方程。數(shù)值離散時,采用中心有限體積格式進行空間離散,在時間上采用混合五步Runge-Kutta格式加以求解。

本節(jié)的主要目的是想知道,在一般的流場解算器中加入了常規(guī)的加速收斂措施以后,再應(yīng)用RPM方法,流場計算的收斂性能是怎么樣的。所以,在流場解算器中我們首先加入了當?shù)貢r間步長和隱式殘值光順等加速收斂措施,然后再運用RPM方法。由于Nstep的取值與當?shù)貢r間步長的計算頻率相關(guān),此處選取最小取值5。Ks取值為10,ε取為1E-12。

以NACA0012翼型為研究對象,來流狀態(tài)為Ma=0.7584,α=0.0°,Re=9.0×106。

圖1展示了加入RPM方法之后數(shù)值模擬所得到的翼型表面壓力分布,并與試驗數(shù)據(jù)[8]進行了對比；圖2為有加入RPM 方法(考慮了多個Krylov接受比Ka取值)和沒有加入RPM方法分別計算得到的殘差收斂過程的對比。

從圖1、圖2中可以看出,在翼型繞流的數(shù)值模擬程序中加入RPM方法所得到的計算結(jié)果與試驗結(jié)果符合得很好,這也證明,RPM方法對數(shù)值過程沒有產(chǎn)生負面的影響；RPM方法對流場數(shù)值計算有明顯的加速收斂效果,即使在流場解算器中加入了當?shù)貢r間步長和隱式殘值光順來加速收斂措施的情況下,也能使總體收斂速度提高近1/3。但在殘差的收斂過程中,相對于原有的數(shù)值計算程序的收斂進程,采用RPM方法的殘差收斂曲線存在相對明顯的數(shù)值震蕩。這主要是由于在翼型繞流的計算時采用了當?shù)貢r間步長的緣故。不同網(wǎng)格點的時間步長是不一樣的,甚至同一網(wǎng)格點在不同迭代步其時間步長也是改變的,而RPM思想是通過監(jiān)視在前面的迭代步的殘差變化遞歸地發(fā)現(xiàn)和分離強勢特征子空間,因此,RPM方法遞歸的正交基提取方式不能很好地反映當前迭代步的時間步長所影響的離散雅可比矩陣的不變子空間。

圖1 翼型表面壓力Fig.1 Pressure distribution of the NACA0012 airfoil

圖2 殘值收斂曲線Fig.2 The curves of residual convergence

5 結(jié) 論

通過以上研究,可以得出如下幾點結(jié)論:

(1)改型的 RPM方法,不僅可以更好地保證正交基的精度,還可以節(jié)省相關(guān)的計算時間。

(2)靜態(tài)特征值的數(shù)值試驗表明,Ks和Nstep是影響RPM方法的加速收斂效果的兩個主要因素,其他因素的影響相對很小。Ka越大,發(fā)現(xiàn)和提取第一組不變子空間正交基的時間越遲,該組正交基的維數(shù)也相對偏小。Nstep值相對越小,可以相對快速地發(fā)現(xiàn)強勢特征值所對應(yīng)的不變子空間,從而以相對較小的總迭代次數(shù)達到收斂標準。

(3)在流場數(shù)值模擬中即使采用了常規(guī)的加速收斂措施以后,再使用RPM方法同樣可以使收斂速度提高近1/3,說明RPM方法在加速收斂方面具有很大的發(fā)展?jié)摿Α?/p>

[1]SHROFF G M,KELLER H B.Stablization of unstable procedures:the recursive projection method[J].Siam J.Numer.Annal,1993,30(4):1099-1120.

[2]PHILIP L.Bifurcations in Kolmogorov and Taylor-voertex flows[D].[Doctor Thesis].California:California Institute of Technology,1999.

[3]JANOVSKY V,LIBERDA O,PRAHA.Continueation of invariant subspacevia the recursive projection method[J].Ap plications of Mathematics,2003,48(4):241-255.

[4]DOROBANTU M,M?LLER J.RPM acceleration of the NSMB compressible flow code[R].TRITA-NA-0118,Royal Institute of Technology,Stockholm,Sweden,2001.

[5]M?LLER J.Aspects of the recursive projection method applied to flow calculations[D].[Doctor Thesis].Sweden Royal Institute of Technology,2005.

[6]STEFAN G,M?LLER J.Recursive projection method for efficient unsteady CFD simulations[A].In:24th-European Congress on Computational Methods in Applied Sciences and Engineering[C].2004.

[7]M?LLER J.An investigation of the recursive projection method on a scalar linear and non-linear partial differential equation[R].TRITA-NA-0118,Royal Institute of Technology,Stockholm,Sweden,2001.

[8]LADSON C L,HILL A S,JOHNSON W G.Pressure distributions from high Reynolds number transonic test of an NACA 0012 airfoil in the langley 0.3-meter transonic cryogenic tunnnel[R].NASA TM 100526,1987.