計算全馬赫數流場的溫度和壓力基本變量方法

柳 建, 李樹民, 金 鋼

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

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

1 壓力原變量控制方程

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

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

(1)

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

(2)

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

(3)

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

ρ=ρPP+ρTT，

(4)

代入連續方程得到壓力的控制方程

ρPPt+ρPTt+mx=0,

(5)

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

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

ρ=P/(RT)

(6)

求出密度.

上面的推導中動量方程保持不變，仍然寫作

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

(7)

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

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

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

(8)

對(5)式離散得

(9)

對(7)式離散得

(10)

計算中，當變量的取值位置與網格預定變量位置不相同時,使用vanAlbada通量限制器進行插值

其中

(11)

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

且

(12)

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

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

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

(13)

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

3 算例結果及分析

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

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.

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

圖 1 Sod激波管算例計算結果

圖 2 馬赫3算例計算結果

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

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

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

圖 4 在t=0.5 s時的流線圖

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

表 1 各種吹速下流場內密度參數值

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

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

4 結論

從守恒格式的Euler方程組推導出了基于壓力原變量的流體控制方程.通過在交錯網格上離散，給出了該方程的完整離散格式和數值求解過程，并通過4個算例驗證了其可以很好地應用于可壓縮流和弱可壓縮流的計算.求解過程簡單，易于推廣到熱、流耦合計算中去.

[1] 朱自強,吳子牛,李津,等. 應用計算流體力學[M]. 北京：北京航空航天大學出版社，1998：1-20.

[2] 張君. 基于笛卡爾網格的預處理方法及應用研究[D]. 南京：南京航天航空大學,2016.

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

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

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

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

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

[8] 韋安陽,樊建人,羅坤,等. 低馬赫數方法在自然對流數值模擬的應用[J]. 工程熱物理學報,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] 王世慶,柳建,張翔,等. 光路交疊對強光傳輸的影響[J]. 光子學報，2008,37(1):173-176.

[11] 柳建,金鋼,王世慶,等. 研究連續激光內通道傳輸的弱可壓縮流模型[J]. 強激光與粒子束，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] 柳建,李樹民,趙杰,等. 鏡面熱變形及吹氣流場對光束的聯合影響[J]. 光學精密工程，2014,22(8):56-62.

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