李武庸, 涂國華, 陳 曦

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

引 言

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

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

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

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

1 基本的Arnoldi算法

展開最后一項得

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的過程, 適用于計算n階巨型非對稱矩陣A的近似特征對(λi，xi)i∈{1，2，…，n}， 將給定的巨型矩陣A逐步化簡為上Hessenberg矩陣H。通過計算m(m<

2 Arnoldi算法的變體

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

2.1 顯式重啟Arnoldit算法

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

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

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

假設(shè)A∈Cn×n可對角化, 且滿足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+， 即它在某些特征向量方向上的分量被加強, 而在其他特征向量方向上的分量被弱化[47]。

例如, 若選取

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

(1)

其中，c是常數(shù),v+是沿a1p(λ1)x1+a2p(λ2)x2+…+anp(λn)xn方向上的向量。因此, 如果λβ接近濾波值μ1，μ2， …，μp中的某個, 而λα不接近; 則相對于xα而言,xβ方向不再被加強。因此, 從K(A，v1)中選取一個較好的重啟向量v+就是選取一些多項式濾波來調(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)

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

則H+=VTHV， 及

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

(3)

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

若已用初始向量v1運行了c步Arnoldi迭代, 有

(4)

AVce1=VcHe1

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

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

式AVce1=VcHe1，表明對于任意的μ有

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

算法流程如圖3所示。

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

2.4 譜位移技術(shù)

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

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

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

(A-σB)-1B

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

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

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

3 數(shù)值實驗

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

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

(5)

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

(6)

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

(7)

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

(8)

其中, 在時間模態(tài)下, 波數(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)代入小擾動方程(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流擾動特征值問題

首先對Poiseuille流,U=1-y2， Orr-Sommer-

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

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

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

其中

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

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

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

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

整理后得矩陣特征值問題

Aφ=cBφ

計算條件為:Re=10 000.0; 流向波數(shù)α=1.0; Chebyshev配置點數(shù)N=100; 基本流分布U(y)=1-y2。 初始向量均為隨機向量, 收斂誤差均為10-10, 計算機設(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流擾動特征值問題的特征譜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種方法求解不可壓縮二維流動比較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算法只有一個最不穩(wěn)定的模態(tài)收斂, 其他模態(tài)并未收斂; 采用隱式重啟Arnoldi算法5個模態(tài)都收斂, 但空間維數(shù)m憑經(jīng)驗選取為25耗時0.445 065 s, 耗時較長; 而采用多重隱式重啟Arnoldi算法5個模態(tài)都收斂且只需一次迭代, 耗時最短。

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

(a) RA

(b) IRA

(c) MIRA圖5 3種不同的Arnoldi方法求解的Poiseuill流擾動特征值問題所對應(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個算例相似的方法, 將流場分解為基本流和擾動場, 代入Navier-Stokes方程, 去掉基本流部分和擾動非線性項, 用Chebyshev配置點法離散后同樣可得到特征值方程如下

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的平板邊界層流動比較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)定性分析

最后考慮計算量更大的二維特征值問題, 以Ma=6高超聲速6°攻角圓錐背風(fēng)流向渦穩(wěn)定性分析為例, 如圖8所示, 等值線是背風(fēng)流向渦的基本流; 云圖是流向擾動速度特征函數(shù)實部分布。該流動情況下的邊界層穩(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分別為法向與展向坐標

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

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

4 結(jié)論

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