郭永恒， 楊 永， 張 強

（西北工業(yè)大學 翼型葉柵空氣動力學國防科技重點實驗室，陜西 西安710072）

0 引 言

在計算流體力學領域，為了高精度地求解非線性偏微分方程組，研究更加復雜的流動現(xiàn)象，間斷Galerkin方法已經(jīng)引起人們越來越多的關注。特別是經(jīng)過Cuckburn和Shu的長期探索，一種具有TVD性質的顯式Runge－Kutta間斷Galerkin（RKDG）格式得以逐步完善［1］，被廣泛應用于雙曲守恒律問題的數(shù)值求解，取得了大量令人滿意的結果，顯示了間斷Galerkin方法的優(yōu)越性。然而，美中不足的是，隨著逼近精度的提高，RKDG格式對應的穩(wěn)定性條件將越來越嚴格，時間步長受到明顯的限制，從而導致更多CPU時間的消耗。對于定常流場的計算問題，盡管可以把當?shù)貢r間步長技術與RKDG格式相結合，在一定程度上加速收斂過程，但是即便如此，最大的時間步長依然受到當?shù)胤€(wěn)定性條件的限制［2］。緩慢的收斂速度在很大程度上制約著RKDG方法在工程中的應用。基于以上分析，本文建立了一種隱式間斷Galerkin（Implicit Discontinuous Galerkin，IMDG）求解器，并通過對翼型亞聲速和跨聲速流場的模擬，檢驗了該求解器的計算效率。

1 間斷Galerkin方程

在二維區(qū)域D上，Euler方程可以寫成如下守恒形式

Q和F＝（F，G）T分別代表守恒型流動變量和通量向量。把D劃分成Ne個互不重疊的子區(qū)域Di，并且設解函數(shù)空間為

P（Di）是定義在Di上的多項式空間。設Di上的測試函數(shù)集合為，它符合如下條件

為了控制數(shù)值離散產(chǎn)生的小量誤差，使用加權殘量法［3］，在每個單元上作如下內(nèi)積運算

運用Green公式，我們由方程（4）得到間斷Galerkin方程

?Dij和n分別表示當前單元的邊界和相關的單位外法向量。在式（5）左邊第三項的積分中，我們使用迎風格式對數(shù)值通量進行計算，實現(xiàn)相鄰單元之間的信息傳遞［4］：

在基函數(shù)的構造過程中，我們應用Gram－Schmidt方法對多項式序列

進行規(guī)范正交化，所得的結果即為數(shù)值格式中使用的基函數(shù)，其具體形式為：

三角形單元經(jīng)過坐標變換后，在計算區(qū)域D＝｛（ξ，η）｜0≤η≤1－ξ，0≤ξ≤1｝上，基函數(shù)滿足如下條件：

這樣，式（5）第一項對應的質量矩陣就成為對角矩陣，使離散格式得到明顯簡化。當采用前三項時，離散格式為二階精度，當采用前六項時，離散格式為三階精度。在本文的數(shù)值實驗中，顯式格式和隱式格式一律設置為二階精度用來比較收斂速度。

2 隱式時間離散格式

我們在時間方向上運用Euler向后差分，并把第n和第n＋1時間層之間物理量X的增量記為δXn，即

那么，在第i號單元上，間斷Galerkin方程（5）就變形為

其中，Ri，l（Qn）為殘差向量，即

設整個流場中守恒型變量為

Q對應的有限元系數(shù)向量為

這樣，Di上的有限元解就可以表示為

對式（10）左側第二、第三項的被積表達式做線性化處理，得

A為系數(shù)矩陣，R為殘差向量

R的分量記為

對比式（10）和式（18）可以發(fā)現(xiàn)，對于定常問題，IMDG格式退化成Newton迭代法，它具有二次收斂階［5］。因為在每一個時間步上，建立系數(shù)矩陣A需要進行大量的數(shù)值積分，所以受文獻［5］的啟發(fā)，我們對式（18）進行一個重要修正，得

其中，

“mod（n－1，k）＝0”表示n－1被k（k≥2）整除，這樣，每進行k次迭代，殘值向量R更新k次，而系數(shù)矩陣A只更新一次，與式（18）相比，計算量大大降低，迭代過程也從整體上得到了進一步的優(yōu)化。

3 線性系統(tǒng)的求解

通過類比，我們把有限體積法中的LU－SGS方法［6］進行推廣，用來求解大型稀疏線性系統(tǒng)（22）。首先，把矩陣A進行如下分解

L、D和U分別為嚴格的分塊下三角矩陣、分塊對角矩陣和分塊上三角矩陣。令矩陣

我們可以根據(jù)如下兩個步驟對線性系統(tǒng)（25）進行快速求解

4 計算結果分析

為了測試IMDG方法的計算效率，我們分別求解了NACA0012翼型對應的亞聲速和跨聲速流場，并與RKDG格式得到的結果作比較。首先，在NACA0012翼型周圍生成非結構網(wǎng)格，如圖1、圖2所示。

對于亞聲速情形，設定計算狀態(tài)為：Ma∞＝0．63，α＝2°。在IMDG格式中，逐步加大CFL數(shù)，可以發(fā)現(xiàn)該格式是無條件穩(wěn)定的，這說明本文的隱式格式具有優(yōu)良的穩(wěn)定性。圖3是壓強系數(shù)分布曲線，可以看出，RKDG和IMDG對應的結果幾乎完全一致，這說明二者的計算精度是相同的。圖4是殘值隨迭代步數(shù)的變化曲線，圖5是殘值隨CPU時間的變化曲線，可以看出IMDG對應的殘值不僅在大范圍內(nèi)單調下降，而且當下降到相同量級時，IMDG使用的迭代步數(shù)和CPU時間分別比RKDG節(jié)省了90%和85%以上，收斂速度幾乎提高了一個數(shù)量級。

為了進一步檢測IMDG求解器的計算效率，我們求解了NACA0012翼型的跨聲速流場，設定計算狀態(tài)為：Ma∞＝0．8，α＝1．25°。在此算例中，本文沒有附加任何限制器，依然能夠得到收斂的結果，如圖6～圖8所示。

可以看出，IMDG能夠捕捉到位置和RKDG完全一致的激波；同時，在高效計算方面，它再一次展示了自身的優(yōu)越性。

5 結 論

為了提高間斷Galerkin方法求解定常流場問題的效率，本文建立了與之相關的隱式離散格式，并在一定程度上對迭代過程進行了優(yōu)化。數(shù)值實驗表明該格式是無條件穩(wěn)定的，這非常有利于計算效率的大幅提高。今后，可以在本文基礎上展開更深層次的探索，特別是簡化隱式格式的建立過程，提高相關的大型稀疏線性系統(tǒng)的求解精度，從而使IMDG的計算效率進一步增強，為間斷Galerkin方法在工程計算中的廣泛應用打下堅實的基礎。

［1］ COCKBURN B，SHU C W．TVD Runge－Kutta local projection discontinuous Galerkin finite element method for scalar conservation laws II：General Framework［J］．Math．Comp，1989，（52）：411－435．

［2］ PATRICK RASETARINERA，HUSSAINI M Y．An efficient implicit discontinuous spectral Galerkin method［J］．Journal of Computational Physics，2001（172）：718－738．

［3］ 王烈衡，許學軍．有限元方法的數(shù)學基礎［M］．北京：科學出版社，2004．（WANG L H，XU X J．The mathematical foundations of the finite element method［M］．Beijing：Science Press，2004．）

［4］ ROE P L．Approximate Riemann solver，parameter vectors and different schemes［J］．Journal of Computational Physics，1981，43：357－372．

［5］ 奧特加J M，萊因博爾特 W C．多元非線性方程組迭代解法［M］．北京：科學出版社，1983．（ORTEGA J M，RHEINBOLDT W C．Iterative solution of nonlinear equations in several variables［M］．Beijing：Science Press，1983．）

［6］ JAMESON A，TURKEL E．Implicit scheme and LU－decompositions［J］．Math．Comput．，1981，37：385－397．