交替LU分裂算法及其在CFD中的應(yīng)用

相 倩 吳頌平

(北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京100191)

徐 悅

(中國航空研究院 航空數(shù)值模擬技術(shù)研究應(yīng)用中心，北京100012)

隨著CFD(Computational Fluid Dynamics)模擬能力的提高，計算中需要求解的方程組的規(guī)模也急劇增加.快速求解大型方程組已經(jīng)成為CFD計算的瓶頸之一［1］.迭代算法是求解大型方程組的有效方法，其收斂速度依賴于方程組系數(shù)矩陣的性質(zhì).針對超聲速流動離散方程組系數(shù)矩陣往往是本質(zhì)正定的這一特點，本文提出了一種迭代算法，稱為交替LU分裂 (ALUS，Alternating Lower-Upper Splitting)算法.該算法計算量小，易于實現(xiàn)，且具有很好的魯棒性，可大大節(jié)省計算時間.

1 迭代算法

常用的迭代算法分為兩類.一類是基于系數(shù)矩陣分裂的古典迭代方法，例如Jacobi方法，Gauss-Seidel方法，松弛法 (SOR，Successive O-ver Relaxation)等.其中交替方向算法因其簡單容易實現(xiàn)，在實際計算中有著廣泛的應(yīng)用［2］.近些年逐步發(fā)展起來的HSS(Hermitian and skew-Hermitian Splitting)方法［3］和 PSS(Positive definite and skew-Hermitian Splitting)方法［4］及其擴展算法也可以看作是廣義的交替方向算法，具體參見文獻［5］.另一類算法是基于變分極值原理的迭代算法，其中最具代表性的有適用于對稱正定系統(tǒng)的共軛梯度法 (CG，Conjugate Gradient)以及適用于非對稱正定的廣義極小殘量法 (GMRES，Generalized Minimal Residual).

1.1 交替方向算法

設(shè)矩陣A∈Cn×n，考慮線性方程組

若A可以分解成A=Q+S，則迭代算法

稱為交替方向算法.該算法亦可以寫成

并且A=M－F.由式 (3)可知，矩陣M可以看成是一個預(yù)處理算子.一個好的預(yù)處理算子M應(yīng)該構(gòu)造簡單易于求逆，并且譜半徑ρ(M－1F)應(yīng)該盡可能小.

對于交替方向算法式 (2)，有以下定理成立.

定理1 對任意的正實數(shù)α＞0，若αI+Q和αI+S都是可逆矩陣，迭代過程式 (2)收斂的充要條件是:ρ［(αI－ Q)(αI+Q)－1(αI－ S)·(αI+S)－1］＜1，并且倘若迭代過程式 (2)收斂，則必收斂到方程組 (1)的解.

1.2 HSS算法和PSS算法

若D、L和U分別代表矩陣A的對角線、嚴格下三角和嚴格上三角部分，則A還可以分解成下三角矩陣和反Hermite矩陣之和，A=P+，其中，P=D+L+UH，=U － UH.

引理1 設(shè)矩陣 P∈Cn×n.如果 P+PH正定，則對任意的實數(shù)α＞0，當αI+P可逆時，‖ (αI－P)(αI+P)－1‖2＜1.

引理2 設(shè)矩陣S∈Cn×n是反Hermite矩陣，則‖ (αI－S)(αI+S)－1‖2=1，?α ＞0.

定理2 若H+HH正定，則對任意實數(shù)α＞0，迭代算法

是收斂的，則稱為HSS迭代算法.

定理3 若P+PH正定，則對任意實數(shù)α＞0，迭代算法

是收斂的，稱其為PSS迭代算法.

共軛斜量法是求解對稱正定線性方程組的有效方法，消去法是求解下三角線性方程組的有效方法.文獻［6］充分利用反對稱矩陣S的特點，給出了快速求解帶位移的反對稱問題 (αI+S)y=的SSS(Solving a Shifted Skew-Symmetric Systems)算法.將這些求解方法結(jié)合起來，使得HSS迭代算法和PSS迭代算法成為有效的廣義交替方向迭代算法［7－8］.盡管如此，HSS算法和PSS算法仍然存在很多問題，例如:當矩陣A的階數(shù)為奇數(shù)時，反對稱矩陣S是奇異的，而參數(shù)α的取值通常比較小，這將導(dǎo)致計算不穩(wěn)定;無論是HSS還是PSS方法都需要求解帶位移的反對稱方程組，但上述SSS算法比較復(fù)雜計算緩慢.為此，本文提出了交替LU分裂算法.

1.3 ALUS算法

將對角矩陣D分解成D1，D2兩部分，則矩陣A可分裂成下三角矩陣=D1+L與上三角矩陣之和，于是由引理1，可得到以下結(jié)論.

是收斂的，本文稱為交替LU分裂算法.

按照式 (3)，記

則有

2 數(shù)值實驗

計算效率以及迭代矩陣的譜半徑是衡量迭代算法優(yōu)劣的重要標準.通過以下算例對本文的ALUS算法進行評估.

2.1 線性問題

在區(qū)域D={0≤x，y≤1}上考慮邊值問題

其中，p＜0;f(x，y)=(3p－5)e2x+y+4xp+p－4，該問題的解析解為u=e2x+y+2x2+y+1.

取自然數(shù)N，令h=1/N，采用均勻網(wǎng)格Δx=Δy=h進行離散區(qū)域D，考慮式 (7)的差分近似

用ALUS算法求解方程組 (8)，初始近似取零向量，迭代的收斂標準是殘差小于10－6.

參數(shù)p決定了線性方程組的系數(shù)矩陣性質(zhì).|p|越小，系數(shù)矩陣越接近對稱矩陣;|p|越大，系數(shù)矩陣的非對稱性越強.參數(shù)N決定了系數(shù)矩陣的規(guī)模.參數(shù)α決定了迭代矩陣的譜分布，影響了迭代的收斂速度，記能使迭代步數(shù)最少的α值為α*，通過數(shù)值實驗，表1給出了不同p和N下α*的取值.

表1 α 的取值

盡管不能給出α*的精確表達式，但從表1中可以看出，當|p|≥10時，α*與|p|近似成正比，與N成反比，可用經(jīng)驗公式表示為α*≈2|p|/N，如圖1所示.

圖1 α*與p、N關(guān)系圖

表2和表3給出了不同p、N和α下，迭代矩陣的譜半徑ρ(M－1ALUSFALUS).由表2和表3可以看到，參數(shù)α可以在很大范圍內(nèi)取值，這表明ALUS算法具有很好的魯棒性.表4給出了p=－1，N=100時，GMRES算法，HSS算法、PSS算法以及ALUS算法收斂所需的CPU時間及迭代步數(shù).由表4可知，本文的新算法，計算所需CPU時間較小，收斂快是一種高效迭代算法.

表2 p=－10，N=32時，迭代矩陣譜半徑

表3 p=－103，N=32時，迭代矩陣譜半徑

表4 不同迭代算法的比較

2.2 二維超聲速圓柱繞流

將ALUS算法應(yīng)用于超音速圓柱繞流的CFD計算中.算例中，來流馬赫數(shù)Ma=8，雷諾數(shù)Re=1.59 ×105.

在計算坐標系下，二維Navier-Stokes方程的形式可以表示為

式中的粘性通量采用中心差分，無粘通量采用Harten-Yee的TVD(Total Variation Diminishing)格式，時間推進則采用了新的ALUS算法.計算網(wǎng)格為結(jié)構(gòu)網(wǎng)格，大小為100×80.迭代的收斂標準是殘差小于10－8.

表5給出了不同參數(shù)α下，ALUS算法的迭代步數(shù)以及CPU計算時間，并與經(jīng)典的LU-SGS(Lower-Upper Symmeffic Gauss-seidel)算法進行對比.可以看出，本文的ALUS算法較經(jīng)典的LU-SGS算法收斂的更快，計算時間甚至可以節(jié)約20%以上;而且α在很廣的范圍內(nèi)，迭代都是收斂的.

為了進一步分析LU分裂算法的計算效率，圖2給出了 CFL值分別為1.3以及2.4時，ALUS算法以及LU-SGS算法的收斂史.

圖2 不同CFL數(shù)下的收斂史

表5ALUS算法與LU-SGS算法的比較

數(shù)值實驗的結(jié)果表明，本文給出的ALUS算法是一種效率高、魯棒性好的迭代算法，可用于CFD的計算.

3 結(jié)論

對于求解大型方程組的交替方向算法，提出了新的ALUS算法.該算法在每步迭代中只需要求解兩個三角方程組，計算量很小.與其他同類迭代算法相比，收斂得更快，因此是一種高效的迭代算法.數(shù)值實驗還表明，算法中的參數(shù)α在很廣的范圍內(nèi)取值，迭代均能收斂，可見新算法魯棒性很好.二維超聲速圓柱繞流的算例則表明，ALUS算法可用于CFD計算.

References)

［1］李良.大型線性方程組求解技術(shù)及在計算電磁學(xué)中的應(yīng)用研究［D］.成都:電子科技大學(xué)數(shù)學(xué)科學(xué)學(xué)院，2009

Li Liang.Study of solutions to large linear systems with applications in computational electromagnetic［D］.Chengdu:School of Mathematical，University of Electronic Science and Technology，Jr of China，2009(in Chinese)

［2］Peaceman D W，Rachford H H，Jr.The numerical solution of parabolic and elliptic differential equations［J］.J Soc Indust and Appl Math，1955，3(1):28－41

［3］Bai Z Z，Golub G H，Nq M K.Hermitian andskew-Hermitian splitting methods for non-Hermitian positive definite linear systems［J］.J Matrix Anal and Appl，2003，24(3):603 －626

［4］Bai Z Z，Golub G H，Lu L Z，et al.Block triangular and skew-Hermitian splitting methods forpositive-definite linear systems［J］.J Sci Comput，2006，26(3):844 －863

［5］蔣爾雄.矩陣計算［M］.北京:科學(xué)出版社，2008:104－120

Jiang Erxiong.Matrix computations［M］.Beijing:Science Press，2008:104－120(in Chinesec)

［6］Jiang Erxiong.Algorithm for solving a shifted skew-symmetric linear system［J］.Front Math China，2007，2(2):227 －242

［7］Benzi M，Guo X P.A dimensional split preconditioner for Stokes and linearized Navier-Stokes equations［J］.Applied Numerical Mathematics，2011，61(1):66 －76

［8］Bai Z Z，Benzi M，Chen F.On preconditioned MHSS iteration methods forcomplex symmetric linear systems［J］.Numer Algorithms，2011，56(2):297 －317