柳 建, 李樹民, 金 鋼

(1. 成都理工大學(xué) 工程技術(shù)學(xué)院， 四川 樂山 614007；2. 中國空氣動(dòng)力研究與發(fā)展中心， 四川 綿陽 621000)

完全可壓縮流和不可壓縮流在計(jì)算流體力學(xué)中有比較完整的研究和應(yīng)用[1-4].近20年來,出現(xiàn)了許多對高低速流采用統(tǒng)一格式求解的研究[5-6]，這種求解方式通常采用條件預(yù)處理或者通量限制器組合的方法.近年來,對低速流的計(jì)算方法研究主要集中在條件預(yù)處理和守恒變量方法上[7-8]，但對流體的可壓縮性研究較少.總的來說，對于弱可壓縮流的研究相對較少.通常在弱可壓縮流的計(jì)算中采用Bossinesq假設(shè)或攝動(dòng)法，這樣一來對低速流體密度變化的捕捉十分困難，而在一些計(jì)算中密度變化量又是很重要的量[9-11]，因而需要研究求解低速弱可壓縮流密度變化的方法.Wesseling等[12]對低馬赫數(shù)下的弱可壓縮流體的計(jì)算進(jìn)行了探討.基于他們的方法，采用全場音速為歸一化因子的壓力原變量方法，以消除在流馬赫數(shù)趨于零時(shí)由于歸一化引起的壓力項(xiàng)奇性和右端源項(xiàng)消失的困難.對這種方法做出完整的推導(dǎo)并給出離散形式，并用幾個(gè)算例進(jìn)行數(shù)值計(jì)算.計(jì)算結(jié)果顯示該方法能夠較好地統(tǒng)一處理從零馬赫數(shù)到較高馬赫數(shù)的流場計(jì)算問題，尤其是能較好處理弱可壓縮流問題.

1 壓力原變量控制方程

基于壓力原變量的全馬赫數(shù)流模型的控制方程，可從Euler方程組得到.一維Euler方程組的守恒形式[1-2]為

Ut+?f(U)/?x=Q，

(1)

首先,從(1)式中消去內(nèi)能項(xiàng)，把能量方程表示為溫度的形式.利用理想氣體狀態(tài)方程

(2)

其中,e為內(nèi)能，γ為比熱容.把連續(xù)方程代入能量方程中，就可以得到能量方程的溫度T表示形式

(3)

繼而利用密度、溫度和壓力的關(guān)系得出壓力的控制方程.由狀態(tài)方程，可以把密度寫成下面的形式

ρ=ρPP+ρTT，

(4)

代入連續(xù)方程得到壓力的控制方程

ρPPt+ρPTt+mx=0,

(5)

其中，ρP=1/(RT)，ρT=-ρ/T.

由于密度是壓力和溫度的函數(shù)，只要得到溫度和壓力的空間分布值就可以利用狀態(tài)方程

ρ=P/(RT)

(6)

求出密度.

上面的推導(dǎo)中動(dòng)量方程保持不變，仍然寫作

mt+(um+P)x=ρfb,

(7)

其中,速度可以通過其與密度和動(dòng)量的關(guān)系得到.到此，(3)、(5)、(6)和(7)式構(gòu)成以溫度、壓力和動(dòng)量為基本變量的封閉方程組.

2 壓力原變量方程組的離散和求解

下面給出壓力原變量方程組的離散求解過程.為避免出現(xiàn)壓力項(xiàng)的鋸齒波問題，因而采用交錯(cuò)網(wǎng)格[3-4].為表述簡便起見，盡量略去一些無關(guān)緊要項(xiàng).令λ為時(shí)間步長與空間步長的比值，對(3)式略去外力項(xiàng)做功后離散得

(8)

對(5)式離散得

(9)

對(7)式離散得

(10)

計(jì)算中，當(dāng)變量的取值位置與網(wǎng)格預(yù)定變量位置不相同時(shí),使用vanAlbada通量限制器進(jìn)行插值

其中

(11)

為得到更高的時(shí)間精度和穩(wěn)定性，使用四階龍格-庫塔法分步求解(8)和(10)式.為簡化計(jì)算，假設(shè)

且

(12)

用Pn代替(10)式中的Pn+1，然后用四階龍格-庫塔法對(8)和(10)式聯(lián)立求解得到Tn+1和m*.

把(12)式代入(9)式得到如下關(guān)于δP的方程

α1δPj-1+α2δPj+α3δPj+1=Fp.

(13)

這是一個(gè)三對角方程組，可以通過高斯-賽德爾迭代法求解.最后把δP代入(12)式得到mn+1.

3 算例結(jié)果及分析

采用壓力原變量方法，對Sod[13]提出的激波管問題和Arora等[14]提出的馬赫3問題進(jìn)行了計(jì)算，初始條件為：

Sod激波管問題case1：t=0,u1=0，P1=1，ρ1=1，u4=0，P4=0.1，ρ4=0.125；

馬赫3問題case2：t=0,u1=0.92，P1=10.333，ρ1=3.857，u4=0，P4=1，ρ4=1.

以上問題中的時(shí)間、壓力、密度都為無量綱的歸一化單位.將計(jì)算值和理論值進(jìn)行比較，分別得到圖1和圖2的結(jié)果.

圖 1 Sod激波管算例計(jì)算結(jié)果

圖 2 馬赫3算例計(jì)算結(jié)果

圖1是case 1在t=0.15時(shí)的結(jié)果，圖2是case 2在t=0.087 5時(shí)的結(jié)果，實(shí)線是理論值，離散點(diǎn)是計(jì)算值.從圖中可以明顯看出，計(jì)算值與理論值的間斷位置符合較好，誤差小于3%，間斷兩側(cè)的密度符合也較好，普遍小于10%;但在小幾何區(qū)間的連續(xù)間斷點(diǎn)處，符合不太好，密度誤差達(dá)到5%，在個(gè)別間斷點(diǎn)超過10%.考慮到其他以密度為原變量的守恒方法的計(jì)算結(jié)果與理論解也存在較大的誤差，因而仍然可以認(rèn)為本方法具有較好的符合度.

對低速弱可壓流，計(jì)算了封閉方腔中的微加熱例子.初始條件為：ρ=1.226kg/m3,T=300K,q=[47 015·exp(-r2/0.036)]2·α,R=297,P=ρRT,v=0m/s，其中，ρ是密度，T是溫度，P是壓強(qiáng)，R是氣體常數(shù)，v是速度，q是方腔內(nèi)中心對稱呈高斯分布的熱源，α=6.5×10-5/m是氣體吸收系數(shù).方腔的邊長是 0.6m，g=9.8m/s2是重力加速度.圖3和圖4給出了在t=0.5 s時(shí)的腔內(nèi)密度等值線分布和流線情況.

圖 3 在t=0.5 s時(shí)的密度等值線分布

圖 4 在t=0.5 s時(shí)的流線圖

在強(qiáng)光控制中，有時(shí)為抑制熱暈效應(yīng)，會(huì)采用吹入或充入惰性氣體的方式[15].其中，在吹氣速度較大時(shí)，流場密度也會(huì)變得不均勻，屬于典型的弱可壓縮問題.對這樣的情況也采用壓力原變量方法進(jìn)行了計(jì)算，表1給出了5、10、20和30 m/s吹氣速度下流場密度的最大、最小值和均方根值(采用初始密度進(jìn)行了歸一化),圖5給出了3種吹速下的密度分布云圖.

表 1 各種吹速下流場內(nèi)密度參數(shù)值

圖 5 吹氣速度為5、10、30 m/s時(shí)的流場密度云圖

從以上2個(gè)完全可壓縮流和2個(gè)低速弱可壓縮流算例的計(jì)算結(jié)果可以看出，壓力原變量方法能夠較好地解出激波層的密度、壓力和速度分布，盡管在激變處與理論值還有一定的偏差，但考慮到模型方程在求解過程中僅使用了二階中心差分格式，并且沒有計(jì)入擴(kuò)散項(xiàng)和湍流脈動(dòng)，這個(gè)結(jié)果還是可以接受的；而在處理低速弱可壓縮問題上，可以獲得細(xì)微的密度變化和流速變化，十分利于求解一些關(guān)注流場密度分布的問題[16].由此可見，該方法可作為一種統(tǒng)一的處理方法來處理工程計(jì)算中的強(qiáng)弱壓縮流問題.

4 結(jié)論

從守恒格式的Euler方程組推導(dǎo)出了基于壓力原變量的流體控制方程.通過在交錯(cuò)網(wǎng)格上離散，給出了該方程的完整離散格式和數(shù)值求解過程，并通過4個(gè)算例驗(yàn)證了其可以很好地應(yīng)用于可壓縮流和弱可壓縮流的計(jì)算.求解過程簡單，易于推廣到熱、流耦合計(jì)算中去.

[1] 朱自強(qiáng),吳子牛,李津,等. 應(yīng)用計(jì)算流體力學(xué)[M]. 北京：北京航空航天大學(xué)出版社，1998：1-20.

[2] 張君. 基于笛卡爾網(wǎng)格的預(yù)處理方法及應(yīng)用研究[D]. 南京：南京航天航空大學(xué),2016.

[3] 陶文銓. 數(shù)值傳熱學(xué)[M]. 2版. 西安：西安交通大學(xué)出版社,2001:1-26.

[4] PATANKAR S V. Numerical Heat Transfer and Fluid Flow[M]. New York:McGraw-Hill,1980:1-31.

[5] 李雪松,顧春偉,徐建中. 新型高低速流統(tǒng)一算法及其葉輪機(jī)應(yīng)用[J]. 航空學(xué)報(bào),2007,28(6):1334-1338.

[6] 李國良,尚慶,袁湘江. 基于全局預(yù)處理方法的全馬赫數(shù)算法研究[C]//第2屆空間工程及信息技術(shù)國際國際會(huì)議(AEIT 2012). HongKong:Hong Kong Education Societ，2012：128-135.

[7] 白曉征,郭正,劉君. 粘性極低速流動(dòng)的預(yù)處理方法研究[J]. 空氣動(dòng)力學(xué)學(xué)報(bào),2009,27(4):451-455.

[8] 韋安陽,樊建人,羅坤,等. 低馬赫數(shù)方法在自然對流數(shù)值模擬的應(yīng)用[J]. 工程熱物理學(xué)報(bào),2011,32(9):1549-1552.

[9] LIU J, WANG S Q. Weakly compressible fluid model to study thermal effects on laser propagating in closed tube[J]. Applied Materials Research,2012,1478(354):165-169.

[10] 王世慶,柳建,張翔,等. 光路交疊對強(qiáng)光傳輸?shù)挠绊慬J]. 光子學(xué)報(bào)，2008,37(1):173-176.

[11] 柳建,金鋼,王世慶,等. 研究連續(xù)激光內(nèi)通道傳輸?shù)娜蹩蓧嚎s流模型[J]. 強(qiáng)激光與粒子束，2006,18(12):1983-1986.

[12] WESSELING P. Principles of Computational Fluid Dynamics[M]. Berlin：Springer-Verlag,2000:578-601.

[13] SOD G A. A survey of several finite difference methods for systems of nonlinear conservation laws[J]. J Comput Phys,1978,27(1):1-31.

[14] ARORA M, ROE P L. A well-behaved TVD limiter for high-resolution calculations of unsteady flow[J]. J Comput Phys,1997,132(1):3-11.

[15] 柳建,李樹民,趙杰,等. 鏡面熱變形及吹氣流場對光束的聯(lián)合影響[J]. 光學(xué)精密工程，2014,22(8):56-62.

[16] 柳建. 近封閉環(huán)境下激光傳輸特性仿真分析平臺(tái)研究[D]. 北京：中國科學(xué)院研究生院，2005.