梁子軒, 丁 玨, 翁培奮

(上海大學(xué)上海市應(yīng)用數(shù)學(xué)和力學(xué)研究所,上海 200072)

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

梁子軒, 丁 玨, 翁培奮

(上海大學(xué)上海市應(yīng)用數(shù)學(xué)和力學(xué)研究所,上海 200072)

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

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

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

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

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

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

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

式中,

對(duì)式 (1)進(jìn)行積分,并應(yīng)用高斯定理將對(duì)流項(xiàng)體積分轉(zhuǎn)化為面積分,得到積分形式的 Navier-Stokes方程為

式中,

i,j分別為沿 x,y,z 3個(gè)坐標(biāo)軸方向的單位向量.

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

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

2 數(shù)值計(jì)算方法

2.1 預(yù)處理技術(shù)

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

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

2.4 多重網(wǎng)格法

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

2.5 邊界條件

本研究采用的邊界條件包括物面邊界條件、遠(yuǎn)場(chǎng)邊界條件和割縫邊界條件,其中物面邊界條件采用無(wú)滑移邊界條件.

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

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

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

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

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

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

3 算例分析

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

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

圖 2 計(jì)算的壓力曲線與實(shí)驗(yàn)值的對(duì)比Fig.2 Calculated pressure curves compared w ith exper imental data

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

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

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

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

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

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

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

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

4 結(jié) 束 語(yǔ)

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

[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ù)值計(jì)算方法及主動(dòng)流動(dòng)控制研究[D].西安:西北工業(yè)大學(xué),2007.

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

[6] 楊愛明,翁培奮,喬志德.用多重網(wǎng)格方法計(jì)算旋翼跨聲速無(wú)粘流場(chǎng)[J].空氣動(dòng)力學(xué)報(bào),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] 楊小權(quán).高效數(shù)值算法在懸停旋翼繞流計(jì)算中的應(yīng)用研究 [D].上海:上海大學(xué),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

國(guó)家自然科學(xué)基金資助項(xiàng)目(10802046)

丁 玨 (1973～),女,副研究員,研究方向?yàn)榄h(huán)境流體力學(xué)、空氣動(dòng)力學(xué).E-mail:dingjue_lu@shu.edu.cn

(編輯:孟慶勛)

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