柳 建, 李樹(shù)民, 金 鋼

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

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

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)方程，可以把密度寫(xiě)成下面的形式

ρ=ρPP+ρTT，

(4)

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

ρPPt+ρPTt+mx=0,

(5)

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

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

ρ=P/(RT)

(6)

求出密度.

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

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

(7)

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

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

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

(8)

對(duì)(5)式離散得

(9)

對(duì)(7)式離散得

(10)

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

其中

(11)

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

且

(12)

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

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

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

(13)

這是一個(gè)三對(duì)角方程組，可以通過(guò)高斯-賽德?tīng)柕ㄇ蠼?最后把δP代入(12)式得到mn+1.

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

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

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

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

以上問(wèn)題中的時(shí)間、壓力、密度都為無(wú)量綱的歸一化單位.將計(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í)線(xiàn)是理論值，離散點(diǎn)是計(jì)算值.從圖中可以明顯看出，計(jì)算值與理論值的間斷位置符合較好，誤差小于3%，間斷兩側(cè)的密度符合也較好，普遍小于10%;但在小幾何區(qū)間的連續(xù)間斷點(diǎn)處，符合不太好，密度誤差達(dá)到5%，在個(gè)別間斷點(diǎn)超過(guò)10%.考慮到其他以密度為原變量的守恒方法的計(jì)算結(jié)果與理論解也存在較大的誤差，因而仍然可以認(rèn)為本方法具有較好的符合度.

對(duì)低速弱可壓流，計(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)中心對(duì)稱(chēng)呈高斯分布的熱源，α=6.5×10-5/m是氣體吸收系數(shù).方腔的邊長(zhǎng)是 0.6m，g=9.8m/s2是重力加速度.圖3和圖4給出了在t=0.5 s時(shí)的腔內(nèi)密度等值線(xiàn)分布和流線(xiàn)情況.

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

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

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

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

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

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

4 結(jié)論

從守恒格式的Euler方程組推導(dǎo)出了基于壓力原變量的流體控制方程.通過(guò)在交錯(cuò)網(wǎng)格上離散，給出了該方程的完整離散格式和數(shù)值求解過(guò)程，并通過(guò)4個(gè)算例驗(yàn)證了其可以很好地應(yīng)用于可壓縮流和弱可壓縮流的計(jì)算.求解過(guò)程簡(jiǎn)單，易于推廣到熱、流耦合計(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] 李國(guó)良,尚慶,袁湘江. 基于全局預(yù)處理方法的全馬赫數(shù)算法研究[C]//第2屆空間工程及信息技術(shù)國(guó)際國(guó)際會(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] 韋安陽(yáng),樊建人,羅坤,等. 低馬赫數(shù)方法在自然對(duì)流數(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] 王世慶,柳建,張翔,等. 光路交疊對(duì)強(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] 柳建,李樹(shù)民,趙杰,等. 鏡面熱變形及吹氣流場(chǎng)對(duì)光束的聯(lián)合影響[J]. 光學(xué)精密工程，2014,22(8):56-62.

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