梁子軒, 丁 玨, 翁培奮

(上海大學上海市應用數(shù)學和力學研究所,上海 200072)

二維翼型繞流計算中預處理和多重網(wǎng)格方法的應用

梁子軒, 丁 玨, 翁培奮

(上海大學上海市應用數(shù)學和力學研究所,上海 200072)

結合 Turkel矩陣預處理方法和多重網(wǎng)格方法,發(fā)展一種適合低速粘性流動計算的高效數(shù)值方法.通過對Navier-Stokes方程的時間導數(shù)項實施 Turkel矩陣預處理,使得可壓縮 Navier-Stokes方程在低速情況下的系統(tǒng)剛性得到改善.為進一步加速收斂,提高計算效率,采用多重網(wǎng)格的 3層 V循環(huán)方式,對 RAE2822超臨界翼型的低速粘性繞流流場進行數(shù)值模擬.計算結果表明,該預處理及多重網(wǎng)格方法能夠大幅度地提高低速粘性流場的收斂性能,具有較高的計算精度,能夠將可壓縮 Navier-Stokes方程的馬赫數(shù)計算區(qū)域擴展到低速不可壓縮區(qū)域,對于低速粘性流場的計算非常有效.

低速流動;多重網(wǎng)格;預處理;隱式求解

20世紀 80年代,Choi等[3]采用預處理方法成功求解了定常粘性流場.國內(nèi)也有諸多學者針對不同的問題對預處理方法進行了深入研究.韓忠華[4]運用預處理方法對旋翼進行了高效的數(shù)值計算;黃典貴[5]嘗試了一種新的預處理方法,成功地進行了可壓縮、不可壓縮、高馬赫數(shù)與低馬赫數(shù)流場的計算.

多重網(wǎng)格方法是一種非常有效的加速收斂技術,對于求解線性橢圓方程是一種最優(yōu)值方法,對求解 Euler方程也有較好的收斂效果.近年來,多重網(wǎng)格法被廣泛應用到加速收斂Navier-Stokes方程的求解中,并取得了顯著的成效.例如,楊愛明等[6]采用多重網(wǎng)格方法有效計算了旋翼跨聲速無粘流場.

基于以上討論,本研究將 Turkel[7]的預處理矩陣和多重網(wǎng)格方法相結合,發(fā)展了一種適合低速流動計算的高效數(shù)值方法——時間積分采用隱式 LUSGS(lower-upper symmetric Gauss-Seidel)方法[8],對流項通量采用 Jameson等[9]提出的有限體積法進行離散,同時對人工粘性項進行修正.通過對RAE2822翼型的低速繞流進行數(shù)值模擬,對比分析在預處理方法和多重網(wǎng)格法相結合時的計算效率,同時對計算所得的流場壓力分布和實驗值進行比較.

1 控制方程及計算網(wǎng)格

微分形式的二維Navier-Stokes方程可表示為

式中,

對式 (1)進行積分,并應用高斯定理將對流項體積分轉化為面積分,得到積分形式的 Navier-Stokes方程為

式中,

i,j分別為沿 x,y,z 3個坐標軸方向的單位向量.

計算網(wǎng)格采用 C形網(wǎng)格,如圖 1所示,網(wǎng)格的大小為 181(弦向)×61(法向),翼型表面占據(jù) 120個網(wǎng)格單元.在 3層多重網(wǎng)格計算中,最稀疏的網(wǎng)格為46×16,其中翼型表面占據(jù) 5個網(wǎng)格單元,第一層網(wǎng)格距離物面的距離約為 1.0×10-6倍弦長,附面層內(nèi)大約有 18～21層左右的網(wǎng)格單元,同時對翼尖和翼根的網(wǎng)格進行局部加密.

圖 1 計算網(wǎng)格Fig.1 C-shapemeshes for RAE2822

2 數(shù)值計算方法

2.1 預處理技術

通常情況下,為保持預處理形式相對簡單,更好地計算粘性問題,一般采用原始變量 Q=(P,u,v,T)T進行求解,即式 (1)可轉換為

本研究采用 Turkel[7]發(fā)展的一種適用于計算粘性流動的預處理矩陣,具體形式為

2.4 多重網(wǎng)格法

多重網(wǎng)格法[6]就是利用一系列粗細不同的網(wǎng)格來求解同一離散問題.將迭代誤差分為高頻波動誤差和低頻光滑誤差,在細網(wǎng)格上利用有效的迭代方法消除高頻波動誤差,利用粗網(wǎng)格消除低頻光滑誤差.通過插值算子在粗網(wǎng)格與細網(wǎng)格間插值,最終達到一致收斂.多重網(wǎng)格法的兩個主要過程就是限制和插值,限制是細網(wǎng)格向粗網(wǎng)格傳值的過程,插值是粗網(wǎng)格向細網(wǎng)格傳值的過程.本研究在限制過程中采用面積加權平均法,插值過程中采用雙線性插值法.有關多重網(wǎng)格法的詳細介紹,請參考文獻[11].

2.5 邊界條件

本研究采用的邊界條件包括物面邊界條件、遠場邊界條件和割縫邊界條件,其中物面邊界條件采用無滑移邊界條件.

黎曼遠場邊界條件理論是建立在方程特征系統(tǒng)的基礎上,適用于未經(jīng)預處理的系統(tǒng).由于經(jīng)過預處理之后,系統(tǒng)的特征值會發(fā)生變化,因此,相應的遠場邊界條件也相應地發(fā)生變化.本研究采用簡單的遠場邊界條件,計算中不涉及跨音速和超音速情況,因而只給出亞音速遠場邊界條件,即

(1)流入為亞音速 (-a＜qn＜0),則

ub=u∞,vb=v∞,Tb=T∞,pb=pint;

(2)流出 (0＜qn＜a),則

ub=uint,vb=vint,Tb=Tint,pb=p∞.

式中,下標 b表示邊界上的值,下標 int表示和邊界相鄰一層的值.

3 算例分析

本研究對 RAE2822翼型的粘性繞流進行數(shù)值模擬.來流馬赫數(shù)M a為 0.01,攻角α為 1.89o,雷諾數(shù) Re為 5.7×106,湍流模型采用BL(baldw in-lomax)湍流模型[12].分別對如下 4種情況進行對比分析:①未采用預處理和多重網(wǎng)格方法;②只采用多重網(wǎng)格方法;③只采用預處理方法;④預處理和多重網(wǎng)格方法均被采用.然后對預處理方法和多重網(wǎng)格方法的收斂性能和計算精度進行討論.結果如圖 2和圖 3所示.

圖 2(a)給出了只采用預處理法和只采用多重網(wǎng)格法計算的壓力曲線.可以看出,采用預處理技術計算的壓力曲線比較光滑而且和實驗值[13]吻合得比較好;而采用多重網(wǎng)格法,計算的壓力曲線不僅不光滑,還在翼尖發(fā)生大幅度的振蕩,并且與實驗值相比,翼型根部和翼型下表面的壓力值和實驗值差別較大,這說明預處理技術能夠有效地提高計算精度.圖 2(b)給出了預處理和多重網(wǎng)格技術均被采用時與兩種方法均未被采用時的壓力分布情況.可以看出,兩種方法均被采用時,計算得到的壓力曲線與實驗值吻合較好,而兩種方法均未被采用時計算的壓力曲線有明顯的振蕩,并且在翼型根部、翼型下表面和翼尖部計算的壓力曲線與實驗值差別較大.圖2(c)給出了上述 4種情況下的壓力曲線變化,可以看出,在預處理和多重網(wǎng)格法均被采用時與只采用預處理方法時計算的壓力曲線基本上是重合的,而且曲線比較光滑,與實驗值吻合很好;而只采用多重網(wǎng)格法和預處理與多重網(wǎng)格法均未被采用時計算的壓力曲線基本上也是重合的,但曲線有明顯的振蕩,與實驗值差別較大.通過對比分析可知,預處理技術能夠有效地改善可壓縮 Navier-Stokes方程的系統(tǒng)“剛性”,使得馬赫數(shù)計算范圍從可壓縮區(qū)域擴展到不可壓縮低速流動區(qū)域,而多重網(wǎng)格法卻不具備這一性能.

圖 2 計算的壓力曲線與實驗值的對比Fig.2 Calculated pressure curves compared w ith exper imental data

圖 3為采用預處理及多重網(wǎng)格方法時的殘值隨工作單元的收斂曲線.表1則給出了 4種情況下收斂性能和計算精度的對比.由圖 3和表1可得如下實驗結果.

圖 3 殘值收斂曲線Fig.3 Residual convergence curve

(1)在未采用預處理或多重網(wǎng)格方法時,殘值振蕩幅度較大,收斂比較困難,在工作單元達到大約11 000個時,殘值下降到 10個量級.結合圖 2可知,這時計算的壓力曲線有明顯的振蕩,且與實驗值差別較大.

(2)在采用了預處理方法后,雖然殘值有較小幅度的振蕩,但在工作單元達到近 4 000個時,殘值下降到 10個量級;與未采用預處理或多重網(wǎng)格法時的計算結果相對比,預處理方法的工作單元加速比約為3,壓力曲線有明顯的改善,且與實驗值吻合很好.

(3)采用多重網(wǎng)格 3V-1-2-3循環(huán)方式時,在工作單元達到約 3 000個時,殘值下降 10個量級;與未采用預處理或多重網(wǎng)格法時的計算結果相比,工作單元加速比約為 4,但二者計算的壓力曲線是一致的,均有明顯的振蕩,且與實驗值差別較大;與只采用預處理方法計算的結果相對比,雖然多重網(wǎng)格方法的收斂性要比預處理方法好些,但計算的壓力曲線不及預處理計算的壓力曲線.

表1 預處理及多重網(wǎng)格法在翼型繞流中的收斂性能和計算精度的對比Table 1 Contrast of convergence capability and accuracy by using precond itioned m ethodsand mutigr id in airfoil flow

(4)預處理和多重網(wǎng)格法均被采用時,收斂性能有大幅度的改善,在工作單元達到約 110個時,殘值很快地降到 10個量級;只采用預處理法或只采用多重網(wǎng)格法或預處理多重網(wǎng)格法均被采用時的工作單元加速比分別約為 27.5,36.7和 100.0.計算的壓力曲線與只采用預處理法的曲線相同,都與實驗值吻合很好.

這些算例說明:①多重網(wǎng)格法能夠有效地提高低速粘性流場的收斂性能,但不能夠改善可壓縮方程在低速不可壓縮流動計算中的方程“剛性”問題,因此,在低速流動情況下,計算精度不及預處理技術;②預處理和多重網(wǎng)格法的結合,不僅能夠改善可壓縮方程的系統(tǒng)“剛性”,而且具有較高的計算精度,可將馬赫數(shù)計算區(qū)域擴展到低速不可壓縮區(qū)域,能夠明顯地提高低速粘性流場的收斂性能.

4 結 束 語

本研究提出的預處理多重網(wǎng)格法能夠改善可壓縮 Navier-Stokes方程的系統(tǒng)“剛性”,大幅度地提高低速粘性流場的收斂性能,并且具有較高的計算精度,能夠將可壓縮 Navier-Stokes方程的馬赫數(shù)計算區(qū)域擴展到低速不可壓縮區(qū)域,對于低速粘性流場的計算非常有效.算例表明,結合預處理與多重網(wǎng)格法,低速粘性流場的收斂性能大幅提高,與只采用預處理法或只采用多重網(wǎng)格法或預處理多重網(wǎng)格法均未被采用時計算的結果相對比,工作單元加速比分別約為36.7,27.5和 100.0,同時能夠顯著提高流場計算精度,對低速粘性流場的計算非常有效.

[1] ROGERS SE,KWAK D.An upwind differencing scheme for the incomp ressible Navier-Stokes equations[J].Applied Numerical Mathematics,1991,8(1):43-46.

[2] CHOI Y H,MERKLE C L. The application of p recondition in viscous flows [J]. Journal of Computational Physics,1993,105:207-223.

[3] CHOI D,MERKLE C L.App lication of time-iterative schemes to incompressible flow[J].A IAA Journal,1985,23(17):1518-1524.

[4] 韓忠華.旋翼繞流的高效數(shù)值計算方法及主動流動控制研究[D].西安:西北工業(yè)大學,2007.

[5] 黃典貴.一種新的流體力學預處理方法[J].工程熱物理學報,2005,26(4):593-595.

[6] 楊愛明,翁培奮,喬志德.用多重網(wǎng)格方法計算旋翼跨聲速無粘流場[J].空氣動力學報,2004,22(3):313-318.

[7] TURKEL E. Preconditioned methods for solving the incompressible and low speed compressible equations[J].Journal of Computational Physics,1987,72:277-298.

[8] VAN LEER B,LEEW T,ROE PL.Characteristic timestepping or local preconditioning of the Euler equations[C]∥A IAA Paper.1991:260-282.

[9] JAMESON A,SCHM IDT W,TURKEL E.Numerical solution of the Euler equations by finite volume methods with Runge-Kutta time stepping schemes[C]∥A IAA Paper.1981:81-1259.

[10] NORBERT K,RAMESH K J.Solution of two-dimensional Euler equations:experience with a finite volume code[R].DLR-FB,1987:87-41.

[11] 楊小權.高效數(shù)值算法在懸停旋翼繞流計算中的應用研究 [D].上海:上海大學,2009:48-50.

[12] BALDWIN B S,LOMAX H.Thin-layer app roximation and algebraic model for separated turbulent flows[C]∥A IAA Paper.1978:78-257.

[13] COOK PH,MCDONALD M A,FIRMAN M C P.Aerofoil RAE 2822-p ressure distributions,and boundary layer and wake measurements[R].Experimental Data Base for Computer Program Assessment,AGARD Report AR-138,1979.

Application of Precond ition ing and M ulti-gr id Techn ique to Two-D im ensional Flow Calculation

L IANG Zi-xuan, D INGJue, WENG Pei-fen
(Shanghai Institute of Applied Mathematics and Mechanics,Shanghai University,Shanghai200072,China)

Turkel preconditioning matrix and multi-grid technique are combined to develop an efficient method suitable for numerical computation of low-speed viscous flow. Rigidity of the system of comp ressible Navier-Stokes equations is improved using the Turkel preconditioning matrix,thus the comp ressible range of Mach of two-dimensional flow is extended for low-speed flows.By using the LUSGS(lower-upper symmetric Gauss-Seidel)implicit method and introducing preconditioning methods together with amultigrid scheme,an RAE2822 airfoil is simulated at low Reynolds number and a small angle of attack.Numerical results show that the preconditioning methods and multi-grid technique can greatly imp rove convergence of viscous flow s.Besides,the above method can extend the computational domain of Mach number of compressible Navier-Stokes equations to a low-speed incompressible area,which is effective in calculating low-speed viscous flows.

low speed flow;multi-grid;preconditioning method;implicitmethod

O 355

A

1007-2861(2011)02-0158-06

10.3969/j.issn.1007-2861.2011.02.010

2009-11-04

國家自然科學基金資助項目(10802046)

丁 玨 (1973～),女,副研究員,研究方向為環(huán)境流體力學、空氣動力學.E-mail:dingjue_lu@shu.edu.cn

(編輯:孟慶勛)

低速不可壓縮流動的數(shù)值模擬一直是計算流體動力學 (computational fluid dynamics,CFD)研究中的活躍領域之一.在低速流動情況下,流場近似為不可壓縮,如果采用可壓縮控制方程來求解,往往會出現(xiàn)收斂緩慢、計算精度低等問題.其原因在于:①此時速度與音速相比很小,條件數(shù)太大;②時間相關的流體力學方程組的對流項系數(shù)矩陣的幾個特征值不處于同一數(shù)量級,方程組奇異.為了使流場的計算能夠快速收斂到定常解,同時使得時間推進方法可以有效地模擬不可壓縮流動,目前主要采用 2種方法:擬壓縮性方法[1]和預處理方法[2].