時間譜方法中的高效GMRES算法

貢伊明,劉戰合,劉溢浪,張偉偉,*

1.西北工業大學 航空學院,西安 710072

2.鄭州航空工業管理學院航空工程學院,鄭州450046

時間譜方法中的高效GMRES算法

貢伊明1,劉戰合2,劉溢浪1,張偉偉1,*

1.西北工業大學 航空學院,西安 710072

2.鄭州航空工業管理學院航空工程學院,鄭州450046

研究了時間譜方法求解周期性非定常流場的計算效率,并對時間譜方法應用于周期性非定常流動的隱式求解方法進行探討。當采樣點數增加或減縮頻率增大時,時間譜方法對應的雅可比矩陣對角占優性質迅速惡化,導致很多傳統的迭代方法失效。為了解決上述問題,論文采用帶預處理的廣義極小殘差(GMRES)算法來提高雅可比系數矩陣的計算收斂性。使用時間譜方法對NACA0012翼型強迫振蕩算例進行計算,并與時域差分方法的計算效率和精度進行對比。研究表明在保證計算精度的同時,時間譜方法普遍可將計算效率提高一個量級左右。對于跨聲速周期性流動,廣義極小殘差算法不論是穩定性還是收斂性都優于對稱SGS迭代算法。

時間譜方法;廣義極小值殘差(GMRES)算法;周期性非定常流動;預處理;計算效率

近些年來,隨著計算機性能的提高,計算流體力學(CFD)在工程中得到廣泛的運用。在實際的工程應用中,很多問題涉及到的非定常流動都具有周期性特征。比如繞直升機旋翼的非定常流動和渦輪機葉片繞流等。在傳統周期性非定常流動的計算方法中,時間項均采用時域差分法(Backward Difference Format,BDF),為保證計算精確性,一個周期往往被分為幾十個時刻點并且要經過若干周期才能得到收斂的周期性非定常解。因此BDF在解決周期性非定常流動問題上計算效率偏低,沒有利用流動的周期性特征。

對于周期性非定常流動,Hall和Crawley[1]提出了一種線化頻域快速計算方法。在該方法中,時域流場變量分解為時均量與非定常一階諧波項。其計算效率很高,但只局限于非定常性較弱和非線性較低的非定常計算中。1998年Ning和He[2]提出了一種非線性計算方法,不同之處在于其通過額外的應力建立時均部分與非定常項之間的關系,通過這種方法能夠求解某些比較復雜和變化劇烈的非定常流動,但該方法仍然無法求解非簡諧周期性流動。隨后Hall等[3]提出了基于傅里葉展開的諧波平衡法。該方法不直接求解流場變量的諧波系數,其在時域內通過譜算子將方程進行時間方向的離散。Hall等將該方法成功應用于模擬壓氣機葉柵繞流。Mc Mullen等[4-5]發展了一種主要在頻域里求解流動控制方程的非線性頻域方法。該方法與諧波平衡法的思想一致只是構造方法有所不同且成功應用到模擬圓柱渦脫落運動和模擬翼型強迫俯仰振蕩等問題。但該方法需要在每步迭代中計算守恒變量和殘值的傅里葉變換與逆變換。類似的還有Dai等[6]提出的時間配置法。Gopinath等[7-8]推導了離散時間導數譜算子的顯式公式并進一步提出時間譜方法(Time Spectral Method,TSM)。該方法利用流動周期性特征,通過傅里葉變換與逆變換將時間導數項通過其他所有時刻點的插值進行離散,從而將非定常問題中各時刻的順序排列求解轉化為多個時刻耦合的定常求解問題。TSM在工程中得到廣泛應用,并取得了很好的實際效果,如直升機旋翼問題[9]、渦輪機葉片繞流[10]以及渦脫落[11]問題。國內對TSM研究相對較少,比較典型的是楊小權等[12]應用基于隱式LUSGS(Lower-Upper Symmetric Gauss-Seidel)迭代的TSM模擬俯仰翼型和機翼的強迫運動問題,謝立軍等[13]應用基于SGS(Symmetric Gauss-Seidel)迭代的TSM計算動導數問題以及詹磊和劉鋒[14]對基于龍格庫塔算法的TSM計算精度進行研究。

2009年Jameson和Shankaran[15]研究發現基于龍格庫塔算法的TSM計算撲翼時會出現不收斂的情況,計算效率不如BDF。另有研究表明,對于非線性強或非簡諧的周期非定常流動問題,TSM則需要大量的采樣點才能擬合,TSM計算效率也大幅下降[16]。上述2種情況的出現是因為隨著減縮頻率增大和采樣點數增加,雅可比矩陣對角占優性質惡化,使得計算效率下降甚至數值發散。

廣義極小殘值(Generalized Minimum RESidual,GMRES)算法[17-18]通過對迭代矩陣進行預處理,降低迭代系數矩陣的條件數,可有效提高高維病態線性方程組的求解效率和魯棒性。本文針對TSM出現的雅可比矩陣對角占優性質惡化的問題,通過引入GMRES算法來改善TSM求解的魯棒性和計算效率。

1 數值方法

1.1 時間譜方法

如果流動具有周期性,則在每一個計算單元內的守恒量U隨時間均呈周期性變化。在一個周期T內設置N 個等時間間距的采樣點,若N為奇數,對守恒量U進行正向傅里葉變換得

式中:Δt=T/N為時間間隔;k為諧波數;ω為角頻率;n代表第n時刻。對應的傅里葉反變換為

將式(2)對時間求導可得

將式(1)代入式(3)中,化簡可得[7]

式中:

若N為偶數,則類似進行計算可得

于是,可將含有N個采樣點在第n個時刻的Navier-Stokes方程的半離散形式寫為

式中:V為單元體的體積;Rn為殘差向量。采用中心格式計算黏性通量與無黏通量。

在式(5)中引入偽時間項,則該方程進一步寫為

式中:τn為偽時間。當顯式離散偽時間項時,為了保證計算穩定性,偽時間步長須限制為

式中:CFL為CFL數;λ為譜半徑;k′為最大諧波數:

當采樣點數增加或減縮頻率增大時,偽時間步長會明顯減小從而嚴重降低計算效率。采用隱式格式可以放寬偽時間步長限制,因此在計算中多采用隱式格式對偽時間項進行離散。對偽時間導數進行一階向后差分離散得到式(8)所示的全隱式格式:

式中:其中:Jn為第n個時刻的通量雅可比矩陣。當矩陣A非對角項為零時,視為進行顯式處理。

1.2 隱式求解穩定性

常用時間推進格式如LU-SGS方法,對稱SGS方法均要求系數矩陣A對角占優。但是TSM增加了非對角項元素會減弱矩陣A的對角占優特性,增加矩陣A條件數。從而導致常用時間推進格式的計算效率降低甚至會發生不收斂情況。非對角項中包含了線性化的空間偏導數項和時間偏導數項。從物理意義上分析,空間通量近似為

式中:S為網格單元表面積;un為當前網格點速度;a為聲速。采用TSM引入非對角項近似為

Atime～Nω (10)則非對角時間導數離散項與原始空間導數離散項之比近似為

從式(11)中可以看出,角頻率ω增大或采樣點數N增加,均導致系數矩陣主對角占優特性變差。若此時仍采用LU-SGS或對稱SGS等要求系數矩陣主對角嚴格占優的方法求解,計算穩定性和收斂性會嚴重降低。采用不同減縮頻率f流動計算的收斂歷史對比如圖1所示。采用不同采樣點數計算的收斂歷史對比如圖2所示。圖1和圖2中縱坐標lg(Residmax)為迭代殘差的對數,從圖中可以明顯看出,當減縮頻率f增大或采樣點數增加時,收斂速度明顯減慢。對于大減縮頻率或采樣點數多的情況,采用對稱SGS迭代方法已經不能保證TSM的高效,需要尋找一種在系數矩陣A相對病態的情況下仍能高效求解的算法來進行時間推進。

1.3 GMRES算法

GMRES以Galerkin原理為基礎,先建立Krylov子空間1組標準正交基,然后在Krylov子空間中利用最小二乘法得到解向量。其中,可以采用Arnoldi方法或修正的Gram-Schmidt方法(Modified Gram-Schmidt,MGS)計算Krylov子空間Km(A,v)中的標準正交基。帶預處理的GMRES算法具體過程為[19]

1)給定初值x(0),然后計算r0=b-Ax(0),對r0預處理,即r0=ML-1r0,記θ=(r0, r0)1/2,得到第一個正交基v1=r0/θ。

2)進入GMRES(m)迭代,j=1,2,…,kmax。

① 計算wj=Avj,并對wj進行預處理,wj=ML-1wj。

② 計算hij=(wj,vi),i=1,2,…,j。

j

③計算

⑤ 可通過hj+1,j＜δ來判定是否已經收斂。δ為給定內迭代收斂標準,其中在沒有滿足此條件下記錄所求得的殘值向量作為初始殘值向量繼續返回步驟1)重新迭代,直至滿足hj+1,j＜δ,停止該時間步的迭代,即為重啟型的GMRES算法。

3)求解方程組的近似解方程組新解可表示為

式中:ym=arg minθe1-Hˉmz2為待定的m維實z型列矢量,θ為殘差范數,Hˉm為Hessenberg矩陣,e1為單位向量,使 θe1-Hˉmz2達到最小的z向量即為所求的ym。具體求解方法見文獻[19]。

在本文中需要對式(8)進行迭代求解。在每個偽時間步將其簡化為線性方程組Ax=b形式,其中x對應式(8)中的ΔU,A對應式(8)中的系數矩陣A,b對應式(8)中等式的右邊項。

對于需要求解的式(8),當左端系數矩陣條件數cond(A)=A A-1較大時,方程組是病態的,此時傳統GMRES算法不具有高速收斂特性,如果要保證計算的高效性,就必須對系數矩陣進行有效的預處理。

所謂預處理,就是將方程組轉換成另外一個同解方程組,而后者比前者更容易求解。預處理的實質就是在不改變方程解的情況下,將矩陣特征值集中在某個值附近,從而改善矩陣的性質,加快迭代法的收斂速度。預處理將線性方程組Ax=b轉換為同解的方程組M-1Ax=M-1b(左預處理)或AM-1Mx=b(右預處理)[20]。

本文采用的是左預處理方法并通過LU-SGS方法對預處理矩陣的逆矩陣和列向量的乘積進行處理。

2 數值模擬及分析

本文對強迫俯仰振蕩的NACA0012翼型進行數值模擬。計算網格如圖3所示,網格單元總數為8 568,節點總數為6 004,翼型表面節點為180,附面層第一層高度為1×10-5,附面層總高度為0.02。用有限體積法求解,空間離散采用中心格式,時間推進采用雙時間步法。

2.1 跨聲速算例

通過跨聲速標模AGARD1對TSM進行驗證并研究不同算法TSM的計算效率。AGARD1算例在迎角較大時會出現跨聲速。

翼型俯仰運動規律為

式中:αm為平均迎角;α0為俯仰運動振幅;ω為俯仰振蕩的角速度。角速度ω與無量綱減縮頻率f的關系為

式中:c為翼型根部弦長;U∞為自由來流速度。記xca為質心位置,表1列出了主要的計算條件。

物理時間離散分別采用TSM和BDF,收斂標準均為5.0×10-7。采用BDF計算時,每個時間周期分為120個時間等份,每一步子迭代最大為1 000步,共進行4個周期的計算。在后面計算中收斂標準以及BDF的子迭代條件保持不變。對于AGARD1算例,圖4給出了TSM與BDF計算得到的升力系數CL和力矩系數Cm并與實驗值[21]對比。從圖中可以看出對于該算例的升力系數CL,TSM僅需要5個采樣點就可以達到與BDF同樣精度。而對于力矩系數,TSM則需要7個采樣點才可以達到與BDF同樣的精度。

表1 跨聲速計算條件Table 1 Transonic computational conditions

圖5給出了NACA0012翼型作俯仰振蕩時用TSM計算的收斂歷史。通過圖5可以看出,采用GMRES算法收斂速度明顯快于對稱SGS迭代。且隨著采樣點數增多,GMRES算法的收斂速度并沒有太大的變化。取采樣點數為4～9,表2對比了采用對稱SGS以及GMRES算法推進的計算效率。從表2中不難看出,隨著采樣點數增加,系數矩陣的條件數增加,對角占優特性減弱,GMRES算法的優勢就體現的更為明顯,甚至采用GMRES算法的計算效率可以達到對稱SGS迭代的6倍。當N=9時,采用對稱SGS迭代計算會發散,對稱SGS迭代已經失效。

表2 不同時間推進格式效率對比Table 2 Comparison of efficiency of different time marching methods

通過對比可以看出,采用基于GMRES算法的TSM相比于BDF方法計算效率提高了9倍左右,因此該方法有很好的應用前景。后面的算例沒有特殊說明,TSM均采用GMRES算法進行偽時間推進。

2.2 二維亞聲速俯仰振蕩

翼型的俯仰運動規律以及計算網格均與跨聲速算例相同,表3列出了亞聲速算例的計算條件。物理時間離散分別采用TSM和BDF。

圖6給出了NACA0012翼型作俯仰振蕩時用TSM計算的收斂歷史,通過圖6可以看出該算例在采樣點數N≤7時其收斂速度并非隨著采樣點數的增加而單調下降,這是因為在采樣點數較少時,雖然Atime有所增加,但Aspatial的影響仍占主導,其收斂速度主要與采樣點定常流場求解特性有關。當N＞7時采用TSM引入的非對角時間導數項Atime占主導作用,具體體現在隨著采樣點數增加,收斂速度呈現明顯的單調下降趨勢。圖7給出了不同采樣點數下的TSM計算效率。從圖中可以明顯看出,對于此算例而言,本文的TSM計算效率遠遠高于BDF。當N在4～11范圍內時,隨著采樣點數增加,TSM計算時間近似呈線性變化。當N=4時TSM的計算效率約是BDF的10倍,當N=12時TSM的計算效率約是BDF的2倍。

表3 亞聲速計算條件Table 3 Subsonic computational conditions

圖8給出了TSM與BDF計算得到的升力系數CL和力矩系數Cm對比。可以看出對于該亞聲速俯仰振蕩算例,當TSM采樣點數取N=5時就可以達到與BDF相當的精度,而計算效率提高了約10倍。因此,對于本算例而言,TSM的計算效率比BDF提高了一個量級。

2.3 非簡諧周期運動算例

對于翼型強迫非簡諧俯仰振蕩算例,本文借鑒文獻[16]中計算條件,通過兩個諧波的周期性運動為例來研究TSM在非簡諧周期運動中的計算效率。采用Euler方程求解。翼型俯仰運動規律為式中:ω1與ω2分別為翼型作強迫俯仰振蕩的兩個頻率;α1與α2為對應的俯仰運動振幅。計算狀態如表4所示。

圖9給出了BDF與TSM計算出的升力和阻力隨時間變化曲線。從圖10中可以看出采用TSM計算時需要至少7個點升力系數才與BDF的計算結果吻合,而要保證阻力系數與BDF計算結果吻合則需要13個點。表5對比了采用13個點的TSM與BDF的計算效率,從表中可以看出,雖然TSM需要13個點,但采用TSM相比于傳統BDF計算效率仍能提高6倍。

表4 非簡諧周期運動計算條件Table 4 Non harmonic periodic motion computational conditions

表5 非簡諧周期運動計算時間對比Table 5 Comparison of computational time in non harmonic periodic motion

因此,對于非簡諧周期性運動,雖然需要較多的采樣點數,但帶預處理的GMRES方法仍然可以保證TSM有較高的計算效率。

3 結 論

本文將GMRES算法應用到時間譜方法中,通過LU-SGS預處理降低系數矩陣條件數,解決了當采樣點數增加或減縮頻率增大所導致的普通迭代法失效的問題。本文采用了跨聲速AGARD1翼型俯仰振蕩、亞聲速翼型簡諧和非簡諧周期運動3個標準算例,計算結果顯示:

1)采用GMRES算法進行求解比對稱SGS迭代計算效率有很大提高,且隨著采樣點數增加,采用GMRES算法效率優勢更明顯。

2)一定范圍內隨著采樣點數增加,采用GMRES算法求解時間譜方法的收斂速度有所減慢但變化不大,計算時間近似呈線性變化。

3)對于線性或非線性較弱的算例,無論是簡諧或非簡諧周期性運動,采用GMRES算法的TSM相比于時域差分法具有明顯優勢,可將計算效率提高一個量級。

[1] HALL K C,CRAWLEY E F.Calculation of unsteady flows in turbomachinery using the linearized Euler equations[J].AIAA Journal,1989,27(6):777-787.

[2] NING W,HE L.Computation of unsteady flows around oscillation blades using linear and nonlinear harmonic Euler methods[J].Journal of Turbomachinery,1998,120(3):508-514.

[3] HALL K C,THOMAS J P,CLARK W S.Computation of unsteady nonlinear flows in cascades using a harmonic balance technique[J].AIAA Journal,2002,40(5):879-886.

[4] MCMULLEN M,JAMESON A,ALONSO J J.Acceleration of convergence to a periodic steady state in turbomachinery flows:AIAA-2001-0152[R].Reston:AIAA,2001.

[5] MCMULLEN M,JAMESON A,ALONSO J J.Application of a nonlinear frequency domain solver to the Euler and Navier-Stokes equations:AIAA-2002-0120[R].Reston:AIAA,2002.

[6] DAI H H,YUE X K,YUAN J P,et al.A time domain collocation method for studying the aeroelasticity of a two dimensional airfoil with a structural nonlinearity[J].Journal of Computational Physics,2014,270:214-237.

[7] GOPINATH A K,JAMESON A.Time spectral method for periodic unsteady computations over two-and threedimensional bodies:AIAA-2005-1220[R].Reston:AIAA,2005.

[8] VAN DER WEIDE E,GOPINATH A K,JAMESON A.Turbomachinery applications with the time spectral method:AIAA-2005-4905[R].Reston:AIAA,2005.

[9] CHOI S,POTSDAM M,LEE K,et al.Helicopter rotor design using a time-spectral and adjoint-based method:AIAA-2008-5810[R].Reston:AIAA,2008.

[10] GOMAR A,BONVY Q,SICOT F,et al.Convergence of Fourier-based time methods for turbomachinery wake passing problems[J].Journal of Computational Physics,2014,278(C):229-256.

[11] GOPINATH A K,JAMESON A.Application of the time spectral method to periodic unsteady vortex shedding:AIAA-2006-0449[R].Reston:AIAA,2006.

[12] 楊小權,程蘇堃,楊愛明,等.基于時間譜方法的振蕩翼型和機翼非定常黏性繞流數值模擬[J].航空學報,2013,34(4):787-797.YANG X Q,CHENG S K,YANG A M,et al.Time spectral method for numerical simulation of unsteady viscous flow over oscillating airfoil and wing[J].Acta Aeronautica et Astronautica Sinica,2013,34(4):787-797(in Chinese).

[13] 謝立軍,楊云軍,劉周,等.基于時間譜方法的飛行器動導數高效計算技術[J].航空學報,2015,36(6):2016-2026.XIE L J,YANG Y J,LIU Z,et al.A high efficient method for computing dynamic derivatives of aircraft based on time spectral method[J].Acta Aeronautica et Astronautica Sinica,2015,36(6):2016-2026(in Chinese).

[14] 詹磊,劉鋒.應用傅里葉時間譜方法求解二維跨音速流動問題解的精度研究[J].航空工程進展,2015,6(4):395-404.ZHAN L,LIU F.Study on accuracy of the solutions using the Fourier time spectral method for two-dimensional transonic flows[J].Advances in Aeronautical Science and Engineering,2015,6(4):395-404(in Chinese).

[15] JAMESON A,SHANKARAN S.An assessment of dualtime stepping,time spectral and artificial compressibility based numerical algorithms for unsteady flow with applications to flapping wings[C]//19th AIAA Computational Fluid Dynamics.Reston:AIAA,2009.

[16] MAVRIPLIS D J,YANG Z.Time spectral method for periodic and quasi-periodic unsteady computations on unstructured meshes[J].Mathematical Modelling of Natural Phenomena,2011,6(3):213-236.

[17] SAAD Y,SCHULTZ M H.GMRES:A generalized minimal residual algorithm for solvingnonsymmetric linear systems[J].SIAM Journal on Scientific and Statistical Computing,1986,7(3):856-869.

[18] JAMESON A,YOON S.Lower-upper implicit schemes with multiple grids for the Euler equations[J].AIAA Journal,1987,25(7):929-935.

[19] 李春娜,葉正寅,王剛.基于二維非結構網格的GMRES隱式算法[J].西北工業大學學報,2007,25(5):630-635.LI C N,YE Z Y,WANG G.GMRES implicit algorithm based on 2D unstructured meshes for solving Euler equations[J].Journal of Northwestern Polytechnical University,2007,25(5):630-635(in Chinese).

[20] 劉巍,張理論,王勇獻,等.計算空氣動力學并行編程基礎[M].北京:國防工業出版社,2013.LIU W,ZHANG L L,WANG Y X,et al.Computational aerodynamics parallel programming basis[M].Beijing:National Defence Industry Press,2013(in Chinese).

[21] LANDON R H.NACA0012 oscillatory and transient pitching:AGARD-R-702[R].Paris:AGARD,1982.

Efficient GMRES algorithm in time spectral method

GONG Yiming1,LlU Zhanhe2,LlU Yilang1,ZHANG Weiwei1,*

1.School of Aeronautics,Northwestern Polytechnical University,Xi'an 710072,China
2.School of Aeronautic Engineering,Zhengzhou University of Aeronautics,Zhengzhou 450046,China

ln this paper,the computational efficiency of the time-spectral method for solving the periodic unsteady flow field is studied,and the implicit method of time spectral method for solving the periodic unsteady flow is discussed.When the number of sampling points increases or the reduced frequency magnifies,the diagonal dominant property of the Jacobian matrix corresponding to the time spectral method deteriorates rapidly,resulting in the failure of many traditional iterative methods.ln order to solve the problems above,the generalized minimum residual(GMRES)algorithm with preprocessing is used to improve the computational convergence of the Jacobian matrix.The time spectral method is used to compute the NACA0012 airfoil forced oscillation,and the computational efficiency and accuracy is compared with that of the time-domain difference method.The results show that the time spectral method can generally improve the computational efficiency an order of magnitude with saturated computational accuracy.For the transonic periodic flow,the GMRES algorithm is superior to SGS iterative algorithm both in stability and computational convergence.

time spectral method;generalized minimum residual(GMRES)algorithm;periodic unsteady flow;preprocess;computational efficiency

2016-10-27;Revised:2016-12-04;Accepted:2016-12-21;Published online:2016-12-26 15:17

URL:www.cnki.net/kcms/detail/11.1929.V.20161226.1517.004.html

s:National Natural Science Foundation of China for Excellent Young Scholars(11622220);Programme of lntroducing Talents of Discipline to Universities(B17037)

V211.3

A

1000-6893(2017)07-120894-09

10.7527/S1000-6893.2016.120894

2016-10-27;退修日期:2016-12-04;錄用日期:2016-12-21;網絡出版時間:2016-12-26 15:17

www.cnki.net/kcms/detail/11.1929.V.20161226.1517.004.html

國家自然科學基金優秀青年基金(11622220);高等學校創新引智計劃(B17037)

*通訊作者.E-mail:aeroelastic@nwpu.edu.cn

貢伊明,劉戰合,劉溢浪,等.時間譜方法中的高效GMRES算法[J].航空學報,2017,38(7):120894.GONG Y M,LlU Z H,LlU Y L,et al.Efficient GMRES algorithm in time spectral method[J].Acta Aeronautica et Astronautica Sinica,2017,38(7):120894.

(責任編輯:李明敏)

*Corresponding author.E-mail:aeroelastic@nwpu.edu.cn