燃氣自由射流的正格式數值模擬

朱孫科，陳二云，馬大為，樂貴高

(1.南京理工大學機械工程學院，江蘇 南京 210094;2.上海理工大學動力學院，上海 200093)

0 引言

在火箭導彈、航空航天等軍事工程技術領域都會涉及到燃氣射流問題，由于射流實驗受測試設備限制且費用高昂，不能提供較詳細的流場細節，高分辨率的數值格式能夠捕捉到燃氣射流復雜流場的波系結構，數值模擬是研究燃氣射流的有效手段。

目前，模擬此類流動的數值格式主要有TVD、NND、ENO和間斷有限元等［1-4］。20世紀 90年代末，Liu等［5-6］提出了一種求解雙曲守恒律的時間和空間均具有二階精度的正格式，其數值通量由一個采用中心差分格式的數值通量和一個包含弱耗散和強耗散的上風格式數值通量組成，時間上利用Runge-Kutta方法進行離散。該格式構造過程簡單，計算經濟性好，是一種性能優良的計算格式。

本文在Liu等人工作的基礎上［5-6］，將雙曲守恒方程的正格式推廣到軸對稱Euler方程的求解，并將其應用于燃氣自由射流的數值模擬。計算結果與實驗照片反映的流動特征基本吻合，與高精度、高分辨率間斷有限元方法的計算結果相當，馬赫盤的位置與實驗測量結果誤差較小，表明該方法具有較強的激波捕捉能力。

1 控制方程和離散方法

1.1 控制方程

在忽略粘性和熱傳導效應的假設下，軸對稱流動的Euler方程組的守恒形式為:

其中:Ω∈R2，T 是時間變量;U=［ρ，ρu，ρv，E］T;F(U)=［ρu，ρu2+p，ρuv，(E+p)u］T;G(U)=［ρv，ρuv，ρv2+p，(E+p)v］T;S(U)=(- ρv/y)［1，u，v，(E+p)/ρ］T;ρ是密度;u和v分別對應于x和y方向的速度分量;p是壓強;E=ρe+ρ(u2+v2)/2，表示單位體積的總能量;e是單位質量的內能，對于理想氣體，狀態方程為p=(γ-1)ρe;γ 是比熱比。

1.2 空間離散

采用正格式方法對控制方程式(1)進行離散，得到如下軸對稱Euler方程組的半離散有限差分方程:

其中，Δx和Δy分別為x和y方向上的空間步長。以通量Fi+1/2，j為例，正格式數值通量的構造過程如下。

1.2.1 弱耗散格式構造

將式(2)中的數值通量Fi+1/2，j采用一個中心差分通量和一個上風通量組合而成，引入限制器L0=R diag(φ0(θk))R-1，限制器函數 φ0(θ)滿足［7］:

可以構造如下的數值通量:

這里|?A|是A=▽F的模，它等于R|Λ|R-1，其中|Λ|=diag(|λk|)。

這樣可以得到數值通量的如下弱耗散格式:

1.2.2 強耗散格式構造

若A=▽F的模|˙A|=R diag(μk)R-1，這里對角矩陣diag(μk)≥|Λ|。引入限制器 L1=R diag(φ1(θk))R-1，其中，φ1(θ)是 minmod 限制器函數滿足［8］:

可以進一步得到數值通量的如下強耗散格式:

1.2.3 正格式數值通量

將以上的弱耗散格式(4)和強耗散格式(6)合并，可以構造如下的正格式數值通量:

1.3 時間離散

為了與空間精度匹配，時間離散采用二階精度的Runge-Kutta 法［9-10］，即:

2 計算結果與分析

通過噴管射流任意子午面截得的一半區域如圖1所示，計算區域取為［-1∶11］×［0∶4］，軸向距離為12de，徑向取4de，其中de為噴口直徑，計算網格為Δx=1/20，Δy=1/40。噴流流動參數和外流動參數如表1所示，表中下標“j”表示噴管出口流動參數，“∞”表示外流流動參數。噴流中心軸AB、噴管固壁表面EF和FG采用鏡面反射邊界條件，噴管出口邊界GH采用表1所示的射流條件，下游邊界BC采用外推，上游邊界DE和外邊界CD采用靜止大氣條件，初場賦予靜止大氣條件，在t=0時刻噴管出口突然噴出射流。

圖1 計算區域Fig.1 Calculation domain

表1 噴流流動參數Table1 Parameters of fluid flow

圖2給出了工況1在某時刻的密度和壓力等值線分布，從圖中可以看到，射流的流場結構類似鉆石狀，馬赫盤略有呈現，入射激波、反射激波、馬赫盤等構成的波胞交替產生，由噴管噴出的燃氣射流射入壓力更低的大氣環境，形成的膨脹波使壓力能夠逐漸過渡到靜止大氣條件，膨脹波穿過射流傳播并在射流邊界處被反射回噴流，這種反射過程在射流下游會多次重復出現，形成準周期性的波胞結構。同時由于噴流與靜止大氣間的對流作用，流場中產生了接觸間斷面。另外，由于Euler方程可以看作是高Reynolds數下的近似，因此邊界層區域附近將出現Kelvin-Helmholtz不穩定性，越往下游不穩定性越明顯。該結果與三階間斷有限元方法的數值結果［4］吻合較好，與相同流動條件下的實驗紋影照片［11］所反映的流場結構基本吻合，如圖2(c)所示。

圖2 工況1等值線分布和紋影照片Fig.2 Contours and schlieren photograph in case 1

圖3 工況2不同時刻密度等值線分布Fig.3 Density contours at different times in case 2

圖3給出了工況2下自由射流在不同時刻(時間為無量綱時間)的密度等值線發展過程。觀察圖3(a)～圖3(e)可知，在射流噴出噴管的初期，形成類似球狀的不斷向外擴張的初始膨脹波，膨脹波在傳播過程中遇到射流邊界后發生反射，形成的入射激波遇到馬赫盤后發生反射，產生反射激波，于是在馬赫盤的邊緣與入射激波和反射激波交匯，并形成三叉激波，出現三波交點和滑移線。隨著時間的推移，流場繼續發展，馬赫盤上游的流場結構相對穩定，波胞結構并未隨時間而改變(如圖3(f)所示)，流場趨于穩定，馬赫盤下游的流場被滑移線分成內部射流以及由射流邊界和滑移線包圍的外部射流，并出現明顯的Kelvin-Helmholtz不穩定性。與工況1的噴流流場結構相比，由于壓力比增加較大，流場中只形成了一個波胞。本文的計算結果第一馬赫盤距噴管出口平面的距離為4.07de，與理論預估值［12］(3.972de)相對誤差為2.4%，與間斷有限元計算值(4.05de)的誤差為0.5%，吻合較好。另外，從圖3(f)可以看到，第一馬赫盤的過渡區間很窄，且三波點清晰可見，表明正格式方法對激波有較強的捕捉能力，得到的馬赫盤位置較精確。

圖4 工況3等值線分布和紋影照片Fig.4 Contours and schlieren photograph in case 3

圖4給出了工況3在某時刻的密度和壓力等值線分布，此時流場中形成了比較穩定的馬赫盤結構，出現了含膨脹波、桶形激波、反射激波、馬赫盤、滑移線、射流邊界等非常復雜的流場波系結構，三波點清晰可見，該復雜波系結構與相同流動條件下的實驗照片［13］所反映的流動特征符合較好。從圖4中可以看到，馬赫盤下游區域出現了波動，已有的數值結果和實驗均表明，自由射流的流場波動現象是普遍存在的［14］。與工況2對比發現，在保持噴管和馬赫數不變的情況下，由于壓力的進一步增大，馬赫盤出現在射流下游的更遠處，并且直徑更大。第一馬赫盤距噴管出口平面的距離為4.68de，與實驗測量值［13］(4.56de)，相對誤差為2.6%，與間斷有限元計算值(4.6de)的誤差為1.7%，吻合較好。

3 結論

將二維正格式發展到軸對稱Euler方程組的數值求解，并對三種欠膨脹壓力比下的燃氣自由射流進行了數值模擬。計算結果與實驗照片所反映的流動特征基本吻合，雖然該方法只有二階精度，但其計算結果可以與三階精度的間斷有限元相比較。利用該方法捕捉到的馬赫盤位置，與實驗測量值和理論預估值的相對誤差小于3%，且得到的燃氣射流復雜流場波系結構較清晰，這表明正格式具有較強的激波捕捉能力和較高的分辨率，在燃氣自由射流的數值模擬中具有一定的應用前景，可以進一步深入研究。

［1］鄭敏，張涵信.無波動、無自由參數的耗散差分格式(NND)在噴流計算中的應用［J］.空氣動力學學報，1989，7(3):273-281.

［2］樂貴高，馬大為，藏國才.火箭底部流動的TVD數值模擬［J］.彈道學報，1995，7(1):35-40.

［3］徐文燦，黃振宇.高精度ENO格式在射流數值模擬中的應用［J］.空氣動力學學報，2001，19(4):401-406.

［4］陳二云，馬大為，趙改平，等.燃氣自由射流的高精度間斷有限元數值模擬［J］.彈道學報，2009，21(1):79-82.

［5］LIU X D，LAX P D.Positive schemes for solving multi-dimensional hyperbolic systems of conservation laws［J］.CFD Journal，1996，5(2):133-156.

［6］LIU X D，LAX P D.Positive schemes for solving multi-dimensional hyperbolic systems of conservation laws II［J］.Journal of Computational Physics，2003，187(2):428-440.

［7］SWEBY P K.High resolution schemes using flux limiters for hyperbolic conservation laws［J］.SIAM Journal on Numerical Analysis，1984，21(5):995-1011.

［8］KURGANOV A，TADMOR E.New high-resolution central schemes for nonlinear conservation laws and convection-dif-fusion equations［J］.Journal of Computational Physics，2000，160(1):241-282.

［9］SHU C W.Total-variation-diminishing time discretizations［J］.SIAM Journal on Scientific and Statistical Computing，1988，9(6):1073-1084.

［10］SHU C W，OSHER S.Efficient implementation of essentially non-oscillatory shock-capturing schemes，II［J］.Journal of Computational Physics，1989，83:32-78.

［11］MATSUDA T，UMEDA Y，ISHII R，et al.Numerical and experimental studies on chocked underexpanded jets［R］.AIAA 87-1378，1987.

［12］張福祥.火箭燃氣射流動力學［M］.北京:國防工業出版社，1988.

［13］TESHIMA K，MORIYA T，MORI T.Visualization of a free jet by a laser induced fluoresense method［J］.Journalof the Japan Society for Aeronautical and Space Sciences，1984，32(5):61-64.

［14］ISHII R，UMEDA Y，YUHI M.Numerical analysis of gasparticle two-phase flows［J］.Journal of Fluid Mechanics，1989，203:475-515.