Arnoldi方法簡述及其在流動(dòng)穩(wěn)定性中的應(yīng)用

李武庸, 涂國華, 陳 曦

(1. 云南大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院, 云南昆明 650500;2. 中國空氣動(dòng)力研究與發(fā)展中心空氣動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室， 四川綿陽 621000;3. 中國空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所, 四川綿陽 621000)

引 言

矩陣的特征值問題可寫為Ax=λx, 其中A∈Rn×n為矩陣,λ∈C為矩陣的特征值,x∈Cn為特征值λ所對(duì)應(yīng)的特征向量; 特征值問題在科學(xué)和工程中應(yīng)用非常廣泛, 如流動(dòng)穩(wěn)定性問題, 樓房、 橋梁及渦輪機(jī)的防震設(shè)計(jì), 電路網(wǎng)與動(dòng)力系統(tǒng)中分支模式分析等[1]。一個(gè)最直接的特征值問題求解方法就是求解方程det(A-λI)=0。 除了一些特殊的矩陣(如三角矩陣)外, 該方法對(duì)于一般矩陣面臨如下困難: 其一, 如果A的階數(shù)比較大, 則方程det(A-λI)=0的計(jì)算量非常大; 其二, det(A-λI)=0的根不穩(wěn)定?；谏鲜鲈? 研究人員通常尋求其他求解方法。目前, 求解矩陣特征值問題的方法大體上可以分為兩類: 變換法和迭代法。變換法是直接對(duì)原矩陣進(jìn)行處理, 通過一系列變換, 使之變成容易求解特征值的形式, 如Jacobi方法、 Givens方法、 QR方法[2]等。

上述基于Arnoldi方法的改進(jìn)都利用了投影矩陣的特征信息, Jia[29]證明了即使n階矩陣A的Ritz值收斂,A的Ritz向量不一定收斂。為了克服這一弊端, 1997年, Jia[15,30-31]提出了改進(jìn)的Arnoldi方法, 此方法通過最小化已收斂的Ritz值對(duì)應(yīng)的剩余范數(shù)來優(yōu)化Ritz向量。2006年曹陶桃[32]應(yīng)用了求解大型非對(duì)稱矩陣特征問題的精化Arnoldi-Chebyshev方法。然而采用基本的Arnoldi算法計(jì)算巨型矩陣的特征值和相應(yīng)的特征向量, 會(huì)出現(xiàn)數(shù)值計(jì)算存儲(chǔ)無限增長的情況, 限制了該方法的適用性。1992年, Sorensen[33]提出隱式重啟Arnoldi方法, 此方法將Arnoldi過程中的迭代數(shù)固定為預(yù)先給定的值, 再利用已更新的初始向量重啟Arnoldi過程。該迭代方案應(yīng)用隱式位移QR迭代法, 避免了計(jì)算初始向量存儲(chǔ)需求量大的問題。2015年Fazeli等[34]基于隱式重啟Arnoldi方法提出了并行隱式重啟Arnoldi方法, 此算法是隱式重啟Arnoldi方法[35]的一個(gè)變體, 它以嵌套的方式生成多個(gè)子空間, 在每個(gè)重啟周期內(nèi)動(dòng)態(tài)地選擇最佳的子空間。2016年Fender等[37]基于隱式重啟Arnoldi方法提出了并行混合Arnoldi算法, 并行性和計(jì)算機(jī)性能對(duì)計(jì)算收斂速度至關(guān)重要, 因此該方法使用了加速器并改進(jìn)文獻(xiàn)[38-39]中所提出的方法。

Arnoldi方法已在流向渦失穩(wěn)、 后掠翼橫流渦失穩(wěn)、 接觸線失穩(wěn)、 三維邊界層首次失穩(wěn)問題等全局流動(dòng)穩(wěn)定性問題求解中得到廣泛應(yīng)用, 相關(guān)研究內(nèi)容可參考綜述[40-41]以及博士論文[42-44]。本文將利用Arnoldi方法的不同變體并結(jié)合譜位移技術(shù)去求解不可壓Poiseuille流動(dòng)穩(wěn)定性問題、 超聲速邊界層穩(wěn)定性問題和高超聲速流向渦穩(wěn)定性問題, 結(jié)果表明結(jié)合譜位移技術(shù)的多重隱式重啟Arnoldi方法求解效率最高。

本文結(jié)構(gòu)如下: 第1部分簡述基本的Arnoldi算法; 第2部分簡述基于Arnoldi算法的幾類變體及譜位移技術(shù); 第3部分通過求解具體的流動(dòng)穩(wěn)定性問題, 比較各方法的有效性。

1 基本的Arnoldi算法

展開最后一項(xiàng)得

span{v1，v2，…，vk}=span(v1，Av1，…，Ak-1v1)

上述過程稱為k步Arnoldi分解[20,33], 可歸納為如下方程

AVk=VkHk+rkekT

其中,Vk=(v1，v2，…，vk)，ek=Ik(:，k)

此Arnoldi分解過程可視為n×n階矩陣A簡化為k×k階上Hessenberg矩陣H的過程, 適用于計(jì)算n階巨型非對(duì)稱矩陣A的近似特征對(duì)(λi，xi)i∈{1，2，…，n}， 將給定的巨型矩陣A逐步化簡為上Hessenberg矩陣H。通過計(jì)算m(m<

2 Arnoldi算法的變體

Arnoldi算法數(shù)值比較穩(wěn)定, 但計(jì)算的每個(gè)基向量需要與前面計(jì)算出的所有基向量進(jìn)行正交化操作, 隨著計(jì)算的進(jìn)行, 運(yùn)算量和存儲(chǔ)量都會(huì)增加。Arnoldi方法還需計(jì)算一個(gè)m階Hessenberg矩陣的特征值和特征向量, 此方法的運(yùn)算量為O(m3)。 因此, 為了節(jié)省存儲(chǔ)和計(jì)算量, 重啟通常是必需的。在實(shí)際應(yīng)用中常采用重啟Arnoldi算法及其變體, 即選取一個(gè)合適的m值, 用新的初始向量v+重新啟動(dòng)Arnoldi過程; 其中v+是所需特征向量對(duì)應(yīng)Ritz向量的線性組合[22,24]。重啟算法的思想是期望新的初始向量v+包含關(guān)于所需特征向量的更多信息。即, 將擴(kuò)大在所需特征向量方向上的系數(shù), 同時(shí)縮小在不需要的特征向量方向上的系數(shù)。為了提高重啟Arnoldi算法的效率, 文獻(xiàn)[24,27,33,36，46]研究了基于顯式重啟Arnoldi的一些加速技術(shù)。

2.1 顯式重啟Arnoldit算法

考慮已運(yùn)行m步的Arnoldi過程, 從Arnoldi向量v1，v2，…，vm中選擇其線性組合為新的初始向量v+， 重新開始運(yùn)行Arnoldi過程, 即:v1=v+。 由于Krylov空間的性質(zhì)

span{v1，v2，…，vk}=span(v1，Av1，…，Ak-1v1)

因此v+的形式可為:v+=p(A)v1, 其中p(A)為次數(shù)不超過m-1的多項(xiàng)式。

假設(shè)A∈Cn×n可對(duì)角化, 且滿足Axi=λixi，i=1∶n。若v1有特征向量展開的形式v1=a1x1+a2x2+…+anxn， 則v+為

x=a1p(λ1)x1+a2p(λ2)x2+…+anp(λn)xn

的系數(shù)倍。若p(λα)>>p(λβ)， 則v+在xα方向上比在xβ方向上占比多。通過選擇p(λ)， 可得相應(yīng)特征方向上的v+， 即它在某些特征向量方向上的分量被加強(qiáng), 而在其他特征向量方向上的分量被弱化[47]。

例如, 若選取

p(λ)=c(λ-μ1)(λ-μ2)…(λ-μp)

(1)

其中，c是常數(shù),v+是沿a1p(λ1)x1+a2p(λ2)x2+…+anp(λn)xn方向上的向量。因此, 如果λβ接近濾波值μ1，μ2， …，μp中的某個(gè), 而λα不接近; 則相對(duì)于xα而言,xβ方向不再被加強(qiáng)。因此, 從K(A，v1)中選取一個(gè)較好的重啟向量v+就是選取一些多項(xiàng)式濾波來調(diào)出矩陣譜中不需要的部分[36,49]。

該算法流程如圖1所示。

圖1 顯式重啟Arnoldi算法流程圖Fig. 1 Flow chart of explicitly restarted Arnoldi algorithm

2.2 隱式重啟Arnoldi算法

隱式重啟Arnoldi算法是1992年Sorensen基于

Arnoldi[33]重啟過程提出的算法, 使用帶位移的QR迭代法[50]隱式確定重啟向量v+[33,35,51]。假設(shè)H∈Rm×m是上Hessenberg矩陣, 選取濾波值μ1，μ2， …，μp為位移QR法的位移, 則位移QR迭代算法如下

H0=H

fori=1∶p

qr(Hi-1-μiI)=ViRi

Hi=RiVi+μiI

end
H+=Hp

(2)

每個(gè)Hi，i∈(1，2，…，p)是上Hessenberg矩陣; 此外, 令

則H+=VTHV， 及

VR=(H-μ1I)…(H-μPI)

(3)

則上述結(jié)果表明, 多項(xiàng)式濾波(1)等價(jià)于(2)。注意到, (3)中的矩陣R是上三角矩陣, 因此V(:，1)=p(H)e1。其中p(λ)是多項(xiàng)式濾波(1)且c=1/R(1，1)[2]。

若已用初始向量v1運(yùn)行了c步Arnoldi迭代, 有

(4)

AVce1=VcHe1

其中,V+=VcV。令v+為矩陣V+的第1列

v+=V+(:，1)=VcV(:,1)
=c·Vc(H-μ1I)…(H-μpI)e1

式AVce1=VcHe1，表明對(duì)于任意的μ有

(A-μ·I)Vce1=Vc(H-μ·I)e1

因此

v+=c·(A-μ1I)…(A-μpI)Vce1
=p(A)v1=Vc·V(:,1)

該算法流程如圖2所示。

圖2 隱式重啟Arnoldi算法流程圖Fig. 2 Flow chart of implicitly restarted Arnoldi algorithm

2.3 多重隱式重啟Arnoldi算法

在隱式重啟Arnoldi方法中子空間的大小憑經(jīng)驗(yàn)選擇, 若選擇不當(dāng)可能會(huì)導(dǎo)致結(jié)果誤差偏大, 方法不收斂。多重隱式重啟Arnoldi方法便是一種改進(jìn)的隱式重啟Arnoldi方法, 改進(jìn)了隱式重啟Arnoldi方法中憑經(jīng)驗(yàn)選擇子空間的大小這一弊端。此方法也稱為多重隱式重啟嵌套子空間Arnoldi方法, 它是基于巨型矩陣A投影在多個(gè)嵌套子空間上而不是單個(gè)子空間上[53], 再利用Arnoldi算法計(jì)算巨型矩陣A在一組嵌套的Krylov子空間上的Ritz元素, 其嵌套子空間為K(mi，v)， (i=1，2，…，l)，K(mi，v)?K(mi+1，v)。 若在這些子空間中計(jì)算的Ritz元素的精度不理想, 則選擇這些子空間中Ritz元素精度“最佳”的所對(duì)應(yīng)的子空間大小記為mbest。若殘差R滿足R(m1)

算法流程如圖3所示。

圖3 多重隱式重啟Arnoldi算法流程圖Fig. 3 Flow chart of multiple implicitly restarted Arnoldi algorithm

2.4 譜位移技術(shù)

在譜空間中位于邊緣處特征值對(duì)應(yīng)的特征向量和其他特征向量有很大差異, 故邊緣處的特征值容易求解。而聚集在一塊的特征值對(duì)應(yīng)的特征向量的方向很接近, 需要大量的迭代步數(shù)才能準(zhǔn)確提取各個(gè)特征向量對(duì)應(yīng)的特征方向。為有效地將所關(guān)注的特征值從譜中提取出來, 需要引入譜位移技術(shù)。譜位移技術(shù)就是在不改變原問題解的情況下, 將原始解對(duì)應(yīng)的譜空間進(jìn)行一定的變換, 把所需的特征值移到譜的邊緣, 以便于求解。推導(dǎo)如下

譜位移技術(shù)的優(yōu)勢大致歸納為以下幾點(diǎn):

(1)可以直接處理一些奇異矩陣問題。如求解廣義特征值問題Ax=λBx時(shí), 若B是個(gè)奇異矩陣。直接求解B-1A不可能, 若引入譜位移σ后, 則可以直接求解

(A-σB)-1B

(2)可以計(jì)算譜內(nèi)部的特征值, 不僅是那些特殊的特征值(模最大、 最小等), 還包括任意區(qū)間上的特征值。

(3)可以加速收斂, 在迭代法求解特征值問題的過程中, 特征值的收斂速度取決于特征值之間的距離, 通過譜位移技術(shù)可以直接將需要的特征值從聚集的特征值中提取出來, 從而加快迭代收斂速度。

下一部分我們利用上述所提3種方法結(jié)合譜位移技術(shù)求解流動(dòng)穩(wěn)定性問題, 經(jīng)比較得出多重隱式重啟Arnoldi算法更容易算出分布在譜邊緣的特征值。Tezuka等[55]較詳細(xì)地介紹了采用基本的Arnoldi方法求解特征值問題的過程, 并在三維不可壓縮流動(dòng)的穩(wěn)定性分析中取得了成功。

3 數(shù)值實(shí)驗(yàn)

利用上面介紹的Arnoldi方法及其變體求解流動(dòng)穩(wěn)定性的中3個(gè)典型問題, 即不可壓Poiseuille流動(dòng)穩(wěn)定性問題、 超聲速邊界層穩(wěn)定性問題和高超聲速流向渦穩(wěn)定性問題[56-57], 以展示各方法的優(yōu)缺點(diǎn)。首先, 為不失一般性, 以不可壓流動(dòng)為例詳細(xì)介紹流動(dòng)穩(wěn)定性方程的推導(dǎo)過程, 可壓流動(dòng)穩(wěn)定性方程的推導(dǎo)過程見附錄A。

只考慮不可壓縮二維流動(dòng)情形(Squire證明不可壓縮流動(dòng)的三維流動(dòng)情形可歸結(jié)為二維情況)[48]

(5)

設(shè)主流方向沿x軸, 引入小擾動(dòng):u=U+u′，v=v′，p=P+p′代入式(5)得

(6)

φ=φ(y)e[iα(x-ct)]

(7)

p′=π(y)e[iα(x-ct)]

(8)

其中, 在時(shí)間模態(tài)下, 波數(shù)α是實(shí)數(shù), 波速c=cr+ici是復(fù)數(shù), 若ci<0， 則為穩(wěn)定模態(tài), 若ci>0， 則為不穩(wěn)定模態(tài), 若ci=0， 則為中性穩(wěn)定模態(tài)。將式(8)代入小擾動(dòng)方程(5)～(7), 得到

-αcφ′+iαUφ′-iαU′φ=

(9)

(10)

(D2-α2)2φ=iαR[(U-c)(φ″-α2φ)-U″φ]

(11)

即Aφ=cBφ， 其中

A=(D2-α2)2-iαR(UD2-Uα2-U″)

B=iαR(α2-D2)

3.1 Poiseuille流擾動(dòng)特征值問題

首先對(duì)Poiseuille流,U=1-y2， Orr-Sommer-

feld方程(11)系數(shù)是關(guān)于x軸的偶函數(shù), 而方程中微分的階數(shù)是偶數(shù), 故特征函數(shù)一定可以表示成奇模態(tài)(奇函數(shù))和偶模態(tài)(偶函數(shù))的線性組合, 其中每個(gè)奇、 偶模態(tài)都是該問題的特解。研究表明奇模態(tài)都是穩(wěn)定的; 對(duì)偶模態(tài), 邊界條件為

φ?(0)=φ′(0)=φ(1)=φ′(1)=0

在對(duì)稱面處, 1階和3階導(dǎo)數(shù)等于零是偶函數(shù)自帶的條件。將特征函數(shù)用N項(xiàng)Chebyshev 多項(xiàng)式來逼近

其中

T2n(y)=cos(2narccos(y))，y∈[-1，1]

是第n階第1類 Chebyshev多項(xiàng)式。為確定N項(xiàng)未知系數(shù), 考慮到有兩個(gè)邊界條件選定N-2個(gè)配置點(diǎn)

yj=cos(jπ/(2N-4))，j=1，2，…，N-2

設(shè)在N-2個(gè)配置點(diǎn)精確滿足 Orr-Sommerfeld 方程, 考慮邊界條件

整理后得矩陣特征值問題

Aφ=cBφ

計(jì)算條件為:Re=10 000.0; 流向波數(shù)α=1.0; Chebyshev配置點(diǎn)數(shù)N=100; 基本流分布U(y)=1-y2。 初始向量均為隨機(jī)向量, 收斂誤差均為10-10, 計(jì)算機(jī)設(shè)備條件為Intel(R) Core(TM) i5-5250U CPU @ 1.60GHz。數(shù)值結(jié)果見圖4, 表1。

(a) RA

(b) IRA

(c) MIRA圖4 給定譜位移參數(shù)σ=0.3情形下3種不同的Arnoldi方法求解的Poiseuille流擾動(dòng)特征值問題的特征譜Fig. 4 Characteristic spectra of three different Arnoldi methods for solving the perturbation eigenvalue problem of poiseuille flow at the spectral displacement parameter σ=0.3

表1 3種方法求解不可壓縮二維流動(dòng)比較Table 1 Comparison of three methods for solving two-dimensional incompressible flow

在圖4中， 子圖(a)RA表示重啟Arnoldi, 子圖(b)IRA表示隱式重啟Arnoldi, 子圖(c)MIRA表示多重隱式重啟Arnoldi。 黑色+為matlab中eig求解的特征譜。 在表1中“*”表示未滿足誤差要求。 結(jié)合圖4與表1可知, 采用重啟Arnoldi算法只有一個(gè)最不穩(wěn)定的模態(tài)收斂, 其他模態(tài)并未收斂; 采用隱式重啟Arnoldi算法5個(gè)模態(tài)都收斂, 但空間維數(shù)m憑經(jīng)驗(yàn)選取為25耗時(shí)0.445 065 s, 耗時(shí)較長; 而采用多重隱式重啟Arnoldi算法5個(gè)模態(tài)都收斂且只需一次迭代, 耗時(shí)最短。

由圖5可知, 重啟Arnoldi算法的最大殘差的模一直在10-1～10-3附近徘徊, 直至280步左右下降至大約10-6處, 但達(dá)到迭代上限一直沒有收斂; 隱式重啟Arnoldi算法最大殘差的模一開始就比較小, 迭代35步就降至10-10處; 而多重隱式重啟Arnoldi算法的最大殘差的模一開始就位于10-10處, 迭代3～5步同樣在10-10周圍浮動(dòng)。

(a) RA

(b) IRA

(c) MIRA圖5 3種不同的Arnoldi方法求解的Poiseuill流擾動(dòng)特征值問題所對(duì)應(yīng)的殘差Fig. 5 Residual corresponding to the perturbation eigenvalue problem of poiseuille flow solved by three different Arnoldi methods

譜位移參數(shù)是特征值的初始猜測值, 由圖6可見, 多重隱式重啟Arnoldi方法在不同位移參數(shù)下, 皆能收斂到有效特征值; 特別地, 在從0.1～0.6的范圍內(nèi), 多重隱式重啟Arnoldi方法都能得到譜位移參數(shù)附近的模態(tài), 說明該方法魯棒性較好。

(a) σ=0.1

(b) σ=0.3

(c) σ=0.5

(d) σ=0.6

(e) σ=0.8

(f) σ=0.9圖6 不同譜位移參數(shù)σ=0.1～0.9情形下多重隱式重啟Arnoldi方法得到的特征譜Fig. 6 Characteristic spectra obtained by multiple implicitly restarted Arnoldi method under different spectral displacement parameters σ=0.1～0.9

3.2 Ma=2.5超聲速平板邊界層的特征值問題

下面考慮Ma=2.5超聲速平板邊界層的特征值問題[52]。采用與第1個(gè)算例相似的方法, 將流場分解為基本流和擾動(dòng)場, 代入Navier-Stokes方程, 去掉基本流部分和擾動(dòng)非線性項(xiàng), 用Chebyshev配置點(diǎn)法離散后同樣可得到特征值方程如下

A(y)q=wB(y)q，q=(u′，v′，p′，T′，w′)

(a) RA

(b) IRA

(c) MIRA圖7 3種方法分別求解Ma=2.5高超聲速邊界層穩(wěn)定性與函數(shù)eig的比較Fig. 7 Comparison of three methods for solving hypersonic boundary layer stability at Ma=2.5 with function eig

表2 3種方法分別求解來流Ma=2.5的平板邊界層流動(dòng)比較Table 2 Comparison of three methods for solving the boundary layer on a flat plate at incoming Mach number Ma=2.5

3.3 Ma=6高超聲速6°攻角圓錐背風(fēng)流向渦穩(wěn)定性分析

最后考慮計(jì)算量更大的二維特征值問題, 以Ma=6高超聲速6°攻角圓錐背風(fēng)流向渦穩(wěn)定性分析為例, 如圖8所示, 等值線是背風(fēng)流向渦的基本流; 云圖是流向擾動(dòng)速度特征函數(shù)實(shí)部分布。該流動(dòng)情況下的邊界層穩(wěn)定性問題歸結(jié)于求解如下二維特征值問題

圖8 Ma=6高超聲速6°攻角圓錐邊界層背風(fēng)流向渦不穩(wěn)定模態(tài)特征函數(shù)分布Fig. 8 Characteristic function distribution of unstable modes of vortices in the leeward side of conical boundary layer at 6° angle of attack and hypersonic speed Ma=6

A2d(y，z)q=αB2d(y，z)q

其中,y，z分別為法向與展向坐標(biāo)

q=(u′，v′，p′，T′，w′，αu′，αv′，αp′，αT′，αw′)， 流向波數(shù)α為待求特征值; 法向邊界條件為:u′=v′=T′=w′=0，展向?yàn)橹芷谶吔鐥l件; 法向和展向分別用4階中心差分離散, 得到巨型稀疏矩陣特征值問題, 分別用上述3種算法求解頻率為 150 kHz 的不穩(wěn)定模態(tài), 結(jié)果如表3所示, 可見3種方法中多重隱式重啟Arnoldi耗時(shí)最短且精度較高(達(dá)到收斂條件時(shí)剩余范數(shù)最小)。

表3 3種方法求解不可壓縮二維流動(dòng)比較Table 3 Comparison of three methods for solving two-dimensional incompressible flow

4 結(jié)論

本文對(duì)重啟Arnoldi方法、 隱式重啟Arnoldi方法、 多重隱式重啟Arnoldi方法作了介紹; 指出重啟Arnoldi算法憑經(jīng)驗(yàn)選擇子空間的大小, 若子空間大小選擇不合適會(huì)導(dǎo)致結(jié)果誤差偏大, 所求的模態(tài)不收斂, 并且該算法還需顯式計(jì)算重啟向量, 隨著計(jì)算的進(jìn)行此算法的運(yùn)算量和存儲(chǔ)量都會(huì)增加。隱式重啟Arnoldi算法有效地改進(jìn)了重啟Arnoldi算法中顯式計(jì)算重啟向量時(shí)所需時(shí)間與空間這一弊端, 但子空間的大小仍憑經(jīng)驗(yàn)選擇。多重隱式重啟Arnoldi算法改進(jìn)了前兩種方法中憑經(jīng)驗(yàn)選擇子空間大小, 此方法迭代次數(shù)較少, 耗時(shí)較短。最后將上述方法應(yīng)用于3個(gè)典型的流動(dòng)穩(wěn)定性問題求解中, 結(jié)果表明結(jié)合譜位移技術(shù)的多重隱式重啟Arnoldi方法求解效率最高。