利用開源CFD程序計(jì)算可壓縮流場(chǎng)

鄭群，李義進(jìn)，孫蘭昕

(哈爾濱工程大學(xué) 動(dòng)力與能源學(xué)院，黑龍江 哈爾濱 150001)

計(jì)算機(jī)計(jì)算能力的迅速增長(zhǎng)加速了CFD技術(shù)的發(fā)展，推進(jìn)了CFD技術(shù)在工程實(shí)際中的應(yīng)用.通過(guò)數(shù)值模擬，可以得到復(fù)雜流場(chǎng)的多種信息，從而為流場(chǎng)內(nèi)部流動(dòng)結(jié)構(gòu)的深入研究提供了有力的手段.

可壓縮全流場(chǎng)流動(dòng)(如葉柵內(nèi)流動(dòng))可以用密度方法或壓力方法進(jìn)行數(shù)值模擬.基于密度方法，求解馬赫數(shù)較高的流動(dòng)問(wèn)題已取得了很大的成功，但對(duì)于馬赫數(shù)很小的流動(dòng)計(jì)算，此法受到諸多限制，可用預(yù)處理方法來(lái)解決［1-2］.基于壓力方法求解流場(chǎng)問(wèn)題，對(duì)馬赫數(shù)的大小沒有限制，即任何馬赫數(shù)的流動(dòng)均是適用的.

1 流體運(yùn)動(dòng)的控制方程組

絕對(duì)坐標(biāo)系下的流體運(yùn)動(dòng)控制方程組在Favre平均(除了壓力和密度本身以外，都用密度加權(quán)平均)下，得到如下形式的Navier-Stokes方程組:

其中:

為了使上述流動(dòng)方程組封閉，需要引入下面的幾個(gè)補(bǔ)充方程.工質(zhì)可看作完全氣體，其狀態(tài)方程為

式中:Cp為定比熱常數(shù).工質(zhì)的分子粘性系數(shù)由薩瑟蘭公式確定:

式中:T0=273.16 K，μ0為一個(gè)大氣壓下0℃時(shí)氣體的動(dòng)力粘性系數(shù)，Ts為薩瑟蘭常數(shù)，與氣體性質(zhì)有關(guān).

2 邊界條件的處理

2.1 進(jìn)出口條件

在可壓縮流體計(jì)算中，進(jìn)口為亞音速時(shí)，一般給定氣體進(jìn)口總壓、總溫，同時(shí)給定氣體進(jìn)口氣流角.根據(jù)氣體動(dòng)力學(xué)理論，總壓力與總溫度的計(jì)算公式為

實(shí)施這一進(jìn)口條件的簡(jiǎn)單方法為:每一步迭代用上一步(第一步用初始化的數(shù)據(jù))的壓力和氣流角計(jì)算氣流進(jìn)口速度，再由總溫的公式計(jì)算溫度.

用k－ε模型計(jì)算湍流時(shí)，進(jìn)口的k及ε值給定.沒有具體的實(shí)驗(yàn)數(shù)據(jù)時(shí)，一般采用:

式中:L=0.07l，Cμ=0.09.

進(jìn)口條件為超音速時(shí)，給定所有相關(guān)量.用kε模型來(lái)計(jì)算湍流時(shí)，進(jìn)口k及ε值的給定與亞音速方法一樣.

2.2 壁面條件

在標(biāo)準(zhǔn)k－ε模型中，采用壁面函數(shù)法處理湍流粘性系數(shù).壁面臨近節(jié)點(diǎn)P的湍流粘性系數(shù)計(jì)算如下:

出口為亞音速時(shí)，一般給定氣體出口的靜壓，其他量的梯度為零.出口為超音速時(shí)，所有量的梯度為零，即

其中:

kp按k方程計(jì)算.k方程壁面條件為方程壁面為

3 Navier-Stokes控制方程組的離散

OpenFOAM［5］運(yùn)用張量方法和面向?qū)ο蠹夹g(shù)處理連續(xù)介質(zhì)問(wèn)題.只要能夠?qū)?wèn)題寫成張量形式的偏微分方程，就可以根據(jù)OpenFOAM提供的基本程序?qū)懸粋€(gè)新的求解器［6-8］，從而求解此問(wèn)題.在OpenFOAM中處理計(jì)算問(wèn)題的流程如圖1.

圖1 OpenFOAM計(jì)算流程Fig.1 Calculation flow process of OpenFOAM

一個(gè)OpenFOAM求解例子一般包括constant、system和0文件夾.constant文件夾包括網(wǎng)格文件、物性參數(shù)和湍流參數(shù)等;system文件夾包括時(shí)間控制、離散格式和求解方法;0文件夾包括初始化參數(shù)和邊界條件.計(jì)算結(jié)果存儲(chǔ)在1000文件夾中(如圖2).

在控制方程組的求解過(guò)程中，離散代數(shù)方程的求解耗時(shí)最多.OpenFOAM用Krylov子空間法(PCG和PBiCG等)和幾何代數(shù)多重網(wǎng)格法(GAMG)進(jìn)行加速求解.

圖2 OpenFOAM計(jì)算目錄結(jié)構(gòu)Fig.2 Calculation directory structure of OpenFOAM

Krylov子空間法是基于Krylov空間正交投影的參數(shù)依賴、非靜態(tài)迭代方法，其思想是在Krylov子空間遍歷搜尋解.即給一個(gè)線性系統(tǒng)Ax=b，這里A是一個(gè)大的、稀疏的n×n非奇異矩陣，根據(jù)標(biāo)準(zhǔn)的Richardson迭代，有

在變換的Krylov子空間上產(chǎn)生近似解:

其中，對(duì)于給定初始矢量x0，有殘量r0=b－Ax0.為了保證收斂性和健壯性，預(yù)處理技術(shù)是Krylov子空間法有效實(shí)施的關(guān)鍵組成部分.本文采用ILU(0)預(yù)處理技術(shù)，將Krylov子空間法中幾種不同代數(shù)方程求解方法.

基于幾何和分析的代數(shù)多重網(wǎng)格法(geometricalgebraic multi-grid solvers，GAMG)是加速線性代數(shù)方程迭代收斂的有效方法，這是目前國(guó)際上代數(shù)多重網(wǎng)格法研究領(lǐng)域中新發(fā)展起來(lái)的方法.其基本思想是在不同層次的網(wǎng)格上循環(huán)迭代，結(jié)合幾何和代數(shù)2種途徑來(lái)濾掉不同波長(zhǎng)的誤差分量，由于各種頻率的誤差分量都可以得到比較均勻的衰減，因而可以加快了迭代收斂速度.

3.1 TVD格式

Harten首次構(gòu)造了二階精度的TVD格式［9］，并且建立了一整套構(gòu)造TVD格式的理論.TVD格式［10］如下:

式中:(φf(shuō))UD為一階的迎風(fēng)格式，(φf(shuō))HD為高階精度的格式，φ(r)是限制器，

3.2 壓力、速度和密度的耦合

3.2.1 離散動(dòng)量方程

將動(dòng)量方程中的壓力梯度項(xiàng)獨(dú)立出來(lái)，可得離散形式為

3.2.2 離散質(zhì)量方程

質(zhì)量方程對(duì)每個(gè)單元做積分，可得離散形式為

或

式中:“*”代表上一步的計(jì)算值，式(23)中右端第4項(xiàng)可忽略不計(jì).

在可壓縮流中，壓力通過(guò)動(dòng)量方程影響速度，通過(guò)狀態(tài)方程影響密度，從而影響質(zhì)量守恒，因此計(jì)算可壓縮流場(chǎng)時(shí)，需要耦合壓力、速度和密度.

3.2.3 壓力方程

由狀態(tài)方程(7)得

將式(24)、(25)代入式(22)，得面通量:將式(20)中的各項(xiàng)插值到單元界面上:

最終可得壓力方程如下:

3.2.4 壓力修正方程

壓力修正量影響密度和速度的修正量:

將式(28)、(29)代入式(23)，得面通量:

當(dāng)流場(chǎng)計(jì)算達(dá)到收斂時(shí)，壓力修正量趨近于0.用動(dòng)量插值法［13］求界面上的速度，以便克服壓力與速度的失耦.為了防止亞松弛因子的數(shù)值影響數(shù)值計(jì)算結(jié)果，動(dòng)量插值時(shí)，動(dòng)量方程不采用亞松弛.

3.3 SIMPLER算法

SIMPLER算法中壓力修正值只用來(lái)修正速度，壓力場(chǎng)是基于當(dāng)前速度場(chǎng)由壓力方程解得.與SIMPLE算法相比，SIMPLER算法多解一個(gè)壓力方程，使壓力場(chǎng)與速度場(chǎng)始終保持協(xié)調(diào)，故整體效率較高.SIMPLER算法計(jì)算步驟見文獻(xiàn)［11-12］

穩(wěn)態(tài)計(jì)算時(shí)，基于壓力的算法不包含時(shí)間導(dǎo)數(shù)項(xiàng)，這時(shí)可以通過(guò)變量松弛和方程松弛使程序計(jì)算順利收斂.變量松弛方程是顯式的:

式中:α是松弛因子.方程松弛是隱式的:

最終可得壓力修正方程:

為了不破壞質(zhì)量守恒，壓力方程和壓力修正方程不能采用方程松弛，而應(yīng)采用變量松弛.

4 計(jì)算結(jié)果

4.1 圓弧凸包流動(dòng)

凸包流動(dòng)是一個(gè)內(nèi)流問(wèn)題.下壁面有凸包的管道無(wú)粘流動(dòng)，長(zhǎng)期以來(lái)被作為檢驗(yàn)程序和計(jì)算方法的經(jīng)典算例.圓弧凸包的弦長(zhǎng)為1.0，高度為0.1(亞音速和跨音速計(jì)算)或0.04(超音速計(jì)算).

4.1.1 亞音速流動(dòng)

給定圓弧凸包通道進(jìn)口處的總溫、總壓及速度方向(馬赫數(shù)為0.5)、出口處的靜壓.壁面采用無(wú)滲透滑移邊界條件.圖3表示計(jì)算區(qū)域中的馬赫數(shù)分布.由圖4可見，馬赫數(shù)關(guān)于圓弧凸包通道的弦線中點(diǎn)具有良好的對(duì)稱性.

圖3 亞音速流動(dòng)計(jì)算區(qū)域中的馬赫數(shù)分布Fig.3 Mach number distribution of subsonic flow calculation region

圖4 亞音速流動(dòng)時(shí)的壁面馬赫數(shù)分布Fig.4 Surface Mach number distribution for subsonic flow

4.1.2 跨音速流動(dòng)

進(jìn)口邊界條件設(shè)定中隱含著馬赫數(shù)為0.675，其余邊界條件的設(shè)定同亞音速流動(dòng).

如圖5(a)所示，激波捕捉不明顯，這是因?yàn)橐浑A上風(fēng)格式將激波摸平.如圖5(b)所示，用TVD格式捕捉激波，其效果比圖5(a)好.而在文獻(xiàn)［12］的計(jì)算結(jié)果中，激波基本被抹平了.這體現(xiàn)了本文采用的TVD格式對(duì)流場(chǎng)中的間斷具有良好的捕捉能力.圖6計(jì)算結(jié)果表明流場(chǎng)中出現(xiàn)了超音速區(qū)，馬赫數(shù)關(guān)于圓弧凸包的弦線中點(diǎn)不再具有良好的對(duì)稱性，最大馬赫數(shù)點(diǎn)向下游方向偏移.

圖5 跨音速流動(dòng)計(jì)算區(qū)域中的馬赫數(shù)分布Fig.5 Mach number distribution of transonic flow calculation region

圖6 跨音速流動(dòng)時(shí)的壁面馬赫數(shù)分布Fig.6 Surface Mach number distribution for transonic flow

4.1.3 超音速流動(dòng)

給定進(jìn)口邊界處的靜溫、靜壓和進(jìn)口速度;出口邊界，外插所有流場(chǎng)變量;上下壁面采用無(wú)滲透滑移邊界條件.以進(jìn)口的條件初始化流場(chǎng).如圖7所示，運(yùn)用上風(fēng)格式計(jì)算所得的結(jié)果與文獻(xiàn)［12］相似，而運(yùn)用TVD格式計(jì)算所得的結(jié)果表明激波捕捉比較清晰，圓弧凸包前沿處有一個(gè)前緣壓縮波系，它同上壁面相交并反射后，又與后緣壓縮波系相交.TVD格式對(duì)馬赫數(shù)等參數(shù)分布梯度很大的流動(dòng)現(xiàn)象具有良好的捕捉特性.如圖8所示的壁面馬赫數(shù)分布中，壓縮波系從端壁發(fā)出，被端壁反射，從而形成陡峭參數(shù)分布.

圖7 超音速流動(dòng)計(jì)算區(qū)域中的馬赫數(shù)分布Fig.7 Mach number distribution of supersonic flow calculation region

圖8 超音速流動(dòng)時(shí)的壁面馬赫數(shù)分布Fig.8 Surface Mach number distribution for supersonic flow

4.2 平面葉柵

本文計(jì)算了一個(gè)平面葉柵內(nèi)部流場(chǎng).計(jì)算網(wǎng)格如圖9.給定進(jìn)口處的總壓169.6 kPa、總溫293 K和速度方向(垂直于進(jìn)口方向)，出口處的靜壓140 kPa.計(jì)算采用高雷諾數(shù)k-ε湍流模型.圖10給出了平面葉柵的馬赫數(shù)分布.

圖9 平面葉柵網(wǎng)格Fig.9 Plane cascade grid

圖10 平面葉柵的馬赫數(shù)分布Fig.10 Mach numer distribution of plane cascade

5 結(jié)論

基于OpenFOAM程序，用SIMPLER算法模擬亞音速、跨音速和超音速條件下的圓弧凸包流場(chǎng)，然后模擬平面葉柵流場(chǎng).結(jié)果表明:

1)文中采用TVD格式能夠很好捕捉激波，并且具有很高的模擬精度;

2)針對(duì)不同馬赫數(shù)進(jìn)口求解，說(shuō)明經(jīng)過(guò)可壓縮推廣的SIMPLER算法適用于任何馬赫數(shù)流場(chǎng)計(jì)算;

3)平板葉柵計(jì)算與其他CFD軟件計(jì)算比較，可利用OpenFOAM程序開發(fā)的求解程序來(lái)替代標(biāo)準(zhǔn)商業(yè)CFD求解軟件包.

本文只對(duì)二維的算例進(jìn)行了考核，今后工作中將增加對(duì)三維計(jì)算的算例的考核，使程序滿足實(shí)際計(jì)算的需要.

［1］黃典貴.基于Roe格式的可壓與不可壓流的統(tǒng)一計(jì)算方法［J］.應(yīng)用數(shù)學(xué)和力學(xué)，2006，27(6):669-675.

HUANG Diangui.Unified computation of flow with compressible and incompressible fluid based on Roe＇s scheme［J］.Applied Mathematics and Mechanics，2006，27(6):669-675.

［2］嚴(yán)明，趙秋穎，梁磊.預(yù)處理方法在葉輪機(jī)械三維數(shù)值模擬中的應(yīng)用［J］.航空動(dòng)力學(xué)報(bào)，2007，22(1):41-47.

YAN Ming，ZHAO Qiuying，LIANG Lei.Application of preconditioning method in 3D numerical simulation of turbomachinery［J］.Journal of Aerospace Power，2007，22(1): 41-47.

［3］周紅梅，蘇莫明，任雄.基于可壓縮SIMPLE算法的葉柵通道湍流流場(chǎng)的數(shù)值模擬［J］.汽輪機(jī)技術(shù)，2007，49 (4):266-268.

ZHOU Hongmei，SU Moming，REN Xiong.Numerical investigation of turbulent flow for cascade passage based on compressible SIMPLE arithmetic［J］.Turbine Technology，2007，49(4):266-268.

［4］伏曉艷，高歌.求解可壓縮流動(dòng)的同位網(wǎng)格SIMPLE方法研究［J］.航空動(dòng)力學(xué)報(bào)，2007，22(10):1674-1677.

FU Xiaoyan，GAO Ge.Study of SIMPLE algorithm for compressible flows on non-staggered grids［J］.Journal of Aerospace Power，2007，22(10):1674-1677.

［5］WELLER H，TABOR G，JASAK H.A tensorial approach to computational continuum mechanics using object-oriented techniques［J］.Computers in Physics，1998，12(6):620-631.

［6］MANGANI L.Development and validation of a C++object oriented CFD code for heat tansfer analysis［C］// ASME-JSME Thermal Engineering and Summer Heat Transfer Conference.Vancouver，2007:1-16.

［7］MANGANI L，ANDREINI A.Application of an object-oriented CFD code to heat transfer analysis［C］//ASME Turbo Expo.Berlin，2008:1-13.

［8］李孫偉.OpenFOAM模擬風(fēng)洞邊界層風(fēng)場(chǎng)［J］.山西建筑，2008，34(23):70-71.

LI Sunwei.Simulating boundary layer in wind tunnel using OpenFOAM［J］.Shanxi Architecture，2008，34(23):70-71.

［9］HARTEN A.High resolution schemes for hyperbolic conservation laws［J］.Journal of Computational Physics，1983，49(3):357-393.

［10］JASAK H.Error analysis and estimation for the finite volume method with applications to fluid flows［D］.London: Imperial College of Science，Technology and Medicine，1996:97-100.

［11］陶文銓.計(jì)算傳熱學(xué)的近代進(jìn)展［M］.北京:科學(xué)出版社，2000:159-200.

TAO Wenquan.Modern development of computational heat tansfer［M］.Beijing:Science Press，2000:159-200.

［12］KARKI K C，PATANKAR S V.Pressure based calculation procedure for viscous flows at all speeds in arbitrary configurations［J］.AIAA Journal，1989，27(9):1167-1174.

［13］RHIE C M，CHOW W L.Numerical study of the turbulent flow past an airfoil with trailing edge separation［J］.AIAA Journal，1983，21(11):1525-1532.

［14］TRAORE P，AHIPO Y M，LOUSTE C.A robust and efficient finite volume scheme for the discretization of diffusive flux on extremely skewed meshes in complex geometries［J］.Journal of Computational Physics，2009，228(14): 5148-5159.

［15］CHEN Hailong，DAI Shaoshi，LI Jia，YAO Xiongliang.Three-dimensional numeriacal simulaition of the flow past a circular cylinder based on LES method［J］.Journal of Marine Science and Application，2009，8(2):110-117.