射流與超聲速來流剪切層作用的數(shù)值模擬研究

鄧放 韓桂來 劉美寬 丁玨 翁培奮 姜宗林

摘要：本文采用高精度格式求解二維Navier-Stokes方程，分別應(yīng)用五階WEND格式、六階中心差分格式離散對流項和黏性項，時間推進(jìn)采用三階Runge-Kutta格式，并應(yīng)用MPI非阻塞式通信進(jìn)行并行計算，模擬了超聲速射流與同向超聲速流動的相互作用，目的在于研究流長產(chǎn)生的兩個剪切層的相互作用，發(fā)現(xiàn)了有一定夾角的剪切層作用產(chǎn)生周期性次生射流的現(xiàn)象。本文研究了射流/超聲速流剪切層相互作用、作用產(chǎn)生的次生射流具有脈沖性等方面，著重對次聲射流的周期性等展開分析和討論。

關(guān)鍵詞：射流;剪切層;超聲速;次聲射流;數(shù)值模擬

中圖分類號：V211 文獻(xiàn)標(biāo)識碼：A

射流和超聲速來流的組合流動是空氣動力學(xué)的重要流動現(xiàn)象，兩者之間相互作用非常普遍，如雙模態(tài)超燃沖壓發(fā)動機(jī)中，為增強(qiáng)氣流和燃料的摻混，燃燒室中使用燃料噴注支板，以達(dá)到燃料與超聲速氣流混合的目的[1]?；鸺l(fā)射過程中，噴管出口的超聲速射流與箭體周圍氣流相互作用，并隨著高度變化形成不同類型的射流特征[2-4];超聲速巡航導(dǎo)彈在飛行過程中，發(fā)動機(jī)尾部噴管形成的射流與彈體周圍氣流形成相互作用[5]。因此，射流、超聲速流動及其剪切層相互作用具有重要的工程實用背景和理論研究價值。

現(xiàn)階段剪切層的研究主要集中采用凹腔、后臺階等模型：Stalling[6]等根據(jù)剪切層的再附形式的不同，最早將凹腔流動的形式分為開口式流動和閉口式流動;Michael[7]等在試驗中發(fā)現(xiàn)不穩(wěn)定模態(tài)下有超聲速剪切層的大尺度卷吸渦結(jié)構(gòu)。Shen[8]等通過大量試驗測量及分析，將后臺階起始渦流的整個發(fā)展過程可以被歸結(jié)為4個階段：（1）第一加速階段;（2）第二加速階段;（3）第二誘導(dǎo)分離渦與第二分離渦階段;（4）破碎階段。Makoto[9]等使用數(shù)值模擬和試驗的方法，研究了后臺階高度對剪切層再附長度的影響。

目前的很多研究還局限于剪切層與激波膨脹波作用、剪切層與邊界層作用等領(lǐng)域，而沒有合適的模型來研究剪切層與剪切層相互作用：Suzuki[10]等使用幾何聲學(xué)理論和DNS這兩種方法研究了二維超聲速剪切層與可壓縮波弱激波膨脹波的相互作用問題，研究發(fā)現(xiàn)可壓縮波在剪切層渦之間的鞍點處發(fā)生泄漏;Kobayachi T[11]等對不同的臺階入口速度分布情況在后臺階流動的影響進(jìn)行了討論。

本文在后臺階處構(gòu)造了超聲速來流與射流的剪切層相互作用的模型。分別應(yīng)用五階WEND格式[12]、六階中心差分格式離散對流項和黏性項[13]，Steger-Warming流通量矢量分裂法[14]對特征值進(jìn)行分裂，時間推進(jìn)采用三階精度Runge-Kutta格式[15]，借助高精度數(shù)值模擬手段，開展流場基本結(jié)構(gòu)和流動規(guī)律研究。

1 物理和數(shù)學(xué)模型

1.1 射流與超聲速來流作用的物理模型

射流與超聲速來流作用的物理模型如圖1所示。超聲速來流經(jīng)過后臺階時，在后臺階拐角處形成剪切層，稱為來流剪切層。射流的上方也產(chǎn)生一剪切層，稱為射流剪切層。由于射流剪切層和來流剪切層的相互作用，射流剪切層與超聲速來流剪切層交匯處一些物質(zhì)會從射流剪切層和來流剪切層之間射出，稱為次生射流。來流剪切層和射流剪切層之間的區(qū)域為回流區(qū)。本文模擬了超聲速后臺流動/超聲速射流相互作用下的流場結(jié)構(gòu)，相應(yīng)的邊界條件見表1。射流出口馬赫數(shù)為3，參數(shù)設(shè)置見表2。超聲速來流馬赫數(shù)為3，具體參數(shù)見表3。

1.2 控制方程

采用可壓縮Navier-Stokes（N-S）方程作為流動控制方程，N-S方程采用來流參數(shù)進(jìn)行無量綱化：

各矢量形式如下：式中：ρ和p分別為氣體的密度和壓力，u、v和w分別為x、y和z方向上的速度分量，為單位體積總能，采用理想氣體狀態(tài)方程p=ρRT，其中R₀=8.314J/（mol·K），Mair為氣體的摩爾質(zhì)量。

直角坐標(biāo)系下的控制方程經(jīng)過Jacobian變換之后，計算空間中的Navier-Stokes方程形式為：

各矢量形式如下：

式中，Jacobian行列式為：

各黏性應(yīng)力分量如下：

各方向上的熱流分量如下：

1.3 計算方法

對流項離散形式如下：

式中的矢通量構(gòu)造如下：

其中：

式中：ω_xk⁺采用Jiang和Shu[16]建議的最佳權(quán)重系數(shù)為：式中：n為一個整數(shù)，通常取1或2，我們根據(jù)Jiang和Shu的建議，取n為2。而組合系數(shù)分別為：

G通量類似，輪換下標(biāo)即可獲得相應(yīng)表達(dá)式。黏性項半離散逼近式為：

其中導(dǎo)數(shù)項可以離散表示為：

本文在時間積分上采用了三階精度Runge-Kutta方法，計算公式為：

2 計算結(jié)果及分析

2.1 基本流動現(xiàn)象

在第一個臺階BC的上方，是馬赫數(shù)Ma=3的超聲速來流。回流區(qū)氣流流動速度較小，回流再附到臺階下游表面的卷挾，來流方向的壁面邊界層延伸形成超聲速來流剪切層。

射流從EF處噴出，首先經(jīng)過三角形狀的等速核心區(qū)，其中溫度、密度等都為常數(shù)，等速核心區(qū)外是射流邊界層。氣體的膨脹一方面使得氣體加速并沿流向壓力降低，另外一方面氣體的過膨脹使得射流壓力低于環(huán)境壓力，導(dǎo)致射流在剪切層邊界被壓縮，形成反射激波或稱攔截激波，射流邊界和攔截激波之間為射流剪切層。

在超聲速來流剪切層和射流剪切層的交匯處，會形成一股新的射流射入回流區(qū)，我們稱為次生射流。研究發(fā)現(xiàn)，次生射流具有周期性，與剪切層的波動有關(guān)。當(dāng)超聲速來流剪切層的水平傾角達(dá)到一定高度時，新的次生射流射出;如果超聲速來流剪切層的水平低于一定高度，次生射流將停止噴射。

在超聲速來流剪切層和射流剪切層相互作用后，兩個剪切層匯聚到一起，形成一股新的剪切層，稱為混合剪切層。由于剪切層兩側(cè)的速度是不匹配的，因此它受到Kelvin-Helmholtz不穩(wěn)定性的影響壓力的波動對剪切層產(chǎn)生新的擾動，在混合剪切層兩側(cè)出現(xiàn)了微弱的激波群結(jié)構(gòu)。這些微弱的激波被包裹在由兩個剪切層相互作用形成的兩個激波之間。

次生射流是由超聲速來流剪切層和射流剪切層之間的相互作用產(chǎn)生的結(jié)構(gòu)現(xiàn)象，同時，由于次生射流的射出，次生射流的前端會產(chǎn)生一道微弱的激波，我們稱之為次生激波。如果次生激波向右移動，則次生激波的波前的超聲速來流剪切層部分會被抬高。另一方面，如果次生激波向左移動，則次生激波的波前的超聲速來流剪切層部分將被壓下。在超聲速來流剪切層的另一側(cè)，二次沖擊波也產(chǎn)生其自身的透射激波。如圖2～圖5所示。

由于次生激波和第二個臺階上壁面的作用，導(dǎo)致界面上產(chǎn)生Richtmyer-Meshkov不穩(wěn)定性。如圖6所示。我們可以發(fā)現(xiàn)在CD和DE的交界處存在明顯的蘑菇結(jié)構(gòu)，這些蘑菇結(jié)構(gòu)是次生激波加速卷人壁面的射流剪切層的物質(zhì)。

2.2 過膨脹射流和欠膨脹射流

我們研究了來流剪切層和射流剪切層的壓比為0.6情況下的流場結(jié)構(gòu)特點。射流壓力低于環(huán)境壓力，因此，此時的射流是過膨脹射流。過膨脹射流和超聲速來流作用于第二個后臺階的位置，從而形成了復(fù)雜的流場結(jié)構(gòu)。在0.1s的時間內(nèi)，射流剪切層出現(xiàn)了振蕩現(xiàn)象，在噴流剪切層和流入剪切層之間也不會產(chǎn)生次生噴流。這是因為如果射流是典型的過度膨脹射流，超聲速來流剪切層提前失穩(wěn)，超聲速來流剪切層和射流剪切層相互作用微弱，此時則不會產(chǎn)生二次射流。

我們還研究了壓比為1.2/2.5/10/30情況，在這5種情況下，射流的壓力高于周圍的壓力，這表明此時的射流是典型的欠膨脹射流。首先，我們以壓力比2.5為例。從0.01s到0.10s，如圖7～圖10所示，共有4個周期，每個周期約為0.021s。我們選擇了第一個周期：從0.015～0.040s。利用組分標(biāo)定的方法，第一個臺階上方的超聲速流動被標(biāo)注為藍(lán)色，第二個臺階左端射出的射流標(biāo)記為紅色。如果射流剪切層和超聲速來流剪切層之間物質(zhì)的顏色變成金色，則意味著在回流區(qū)內(nèi)噴入了超聲速來流和射流的混合物質(zhì)。這種方法可以很好地直觀描述次生射流的產(chǎn)生過程。次生射流注入回流區(qū)后，射流水平傾角變小，射流剪切層壓低，次生射流開始停止射入，隨后射流剪切層再次抬高直到新的次生射流射入，新的周期產(chǎn)生，從而循環(huán)往復(fù)。

2.3 次生射流的周期性

為了研究不同壓力比對次生射流的影響，我們選擇不同的壓力比為2.5，10，30。由于次生射流對回流區(qū)壓力變化有影響，我們選擇圖10所示的壁面5個壓力點來觀察次生射流的周期性。

如圖11所示，這5個點的壓力變化幾乎是相同的頻率，因此我們可以從數(shù)據(jù)角度準(zhǔn)確觀察次生射流。

如圖12（a）所示，我們可以準(zhǔn)確觀察4個周期的間隔。它們的周期大概在0.025s，周期頻率在50Hz左右。這意味著當(dāng)超聲速來流和射流壓比變大時，整個周期波動的頻率幾乎沒有變化。由于射流剪切層局部也會震動，我們可以觀察到主周期也有其他細(xì)微的波動，當(dāng)壓力比較小時，射流剪切層整體周期性波動的影響處于主導(dǎo)地位，同時，射流剪切層局部振動影響相對較小但存在。

然而，如圖12（b）所示，當(dāng)壓力比為30時。除了主循環(huán)外，我們可以注意到有很大的波動性。壓比變大時，射流剪切層局部振動影響變大，甚至超聲速來流剪切層出現(xiàn)失穩(wěn)現(xiàn)象。每個主循環(huán)大約為0.025s。即使壓力比變化，壓力點的周期不會改變。射流剪切層整體的震動形成了主要周期和次生射流。由于射流剪切層的局部振動和超聲速來流剪切層可能會提前失穩(wěn)，也會導(dǎo)致來流區(qū)壁面附近的壓力變化。但是，我們?nèi)匀豢梢杂^察到完整的循環(huán)。盡管壓力比變化，主循環(huán)仍然約為0.025s，周期頻率約為50Hz。

次生射流的周期性與壓力比變化之間存在一定的關(guān)系，但壓比變大會使得剪切層失穩(wěn)更加顯著，這就導(dǎo)致了射流剪切層的局部振動和超聲速來流剪切層的提前失穩(wěn)。

3 結(jié)論

通過分析，可以得出以下結(jié)論：

（1）射流剪切層與超聲速來流剪切層相互作用，由于相互作用產(chǎn)生次生射流和混合剪切層。兩個弱的沖擊波位于混合剪切層的兩側(cè)。由于混合剪切層兩側(cè)的速度不一致，混合剪切層受到Kelvin-Helmholtz不穩(wěn)定性的影響，然后失穩(wěn)形成一些Kelvin-Helmholtz渦旋。

（2）次生射流具有周期性。周期性與射流剪切層的整體振蕩有關(guān)?；亓鲄^(qū)壁面附近的壓力變化受三個因素影響：射流剪切層整體振動、射流剪切層局部振動、超聲速來流剪切層提前失穩(wěn)。其中，射流剪切層整體振動是回流區(qū)壓力變化的主要因素。隨著壓比的增加，其他兩個因素的影響也增強(qiáng)。這意味著壓力的增加導(dǎo)致兩個剪切層之間相互作用的增強(qiáng)。

（3）現(xiàn)有工況下次生射流的周期約為0.025s，頻率約為50Hz。即使壓力比發(fā)生變化，射流剪切層波動的主要周期基本保持不變。因此，次生射流的周期并不隨壓力比的變化而變化。

參考文獻(xiàn)

[1]宋岡霖，田亮，江強(qiáng)，等.超燃沖壓發(fā)動機(jī)支板熱性能研究[J].推進(jìn)技術(shù)，2013，34（06）：809-816.

[2]Reardon J E，Nelson H F.Rocket plume base heatingmethodology[J].Journal of Thermophysics and Heat Transfer，1994，8（2）：216-221.

[3]Ebrahimi H B.Numerical investigation of jet interaction in asupersonic freestream[J].Journal of Spacecraft and Rockets，2008，45（1）：95.

[4]Troyes J，Dubois I，Borie V，et al.Multi-phase reactivenumerical simulations of a model solid rocket exhaust jet[C]//42^ndAIAA/ASME/SAE/ASEE Joint Propulsion Conference&Exhibit，2006： 4414.

[5]Michael R.Infrared signature generation of airborne targets[R].Proceedings of the SPIE（50277-786X），1993：114-122.

[6]Stallings R L，Wilcox F J.Experimental cavity pressuredistributions at supersonic speeds[R].NASA TP-2683，1987.

[7]Kegerise M A，Spina E F，Garg S，et al.Mode-switching andnonlinear effects in compressible flow over a cavity[J].Physicsof Fluids，2004，16（3）：678-687.

[8]Shen G X，Ma G Y.The investigation on the properties andstructures of starting vortex flow past a backward-facing step byWBIV technique[J].Experiments in Fluids，1996，21（1）：57-65.

[9]Oki M，Iwasawa T，Suehiro M，et al.Numerical simulationwithout turbulence model for backward-facing step flow[J].JSME International Journal Series B：Fluids and ThermalEngineering，1993，36（4）：577-583.

[10]Suzuki T，Lele S K.Shock leakage through an unsteady vortex-laden mixing layer：Application to jet screech[J].Journal ofFluid Mechanics，2003，490：139-167.

[11] Kobayashi T，Togashi S.Comparison of turbulence modelsapplied to backward-facing step flow[J].JSME InternationalJournal Series B：Fluids and Thermal Engineering，1996，39（3）：453-460.

[12] Jiang G S，Shu C W.Efficient implementation of weighted ENOschemes[J]·Journal of Computational Physics，1996，126（1）：202-228.

[13]傅德薰，馬延文，李新亮.可壓縮湍流直接數(shù)值模擬[M].北京：科學(xué)出版社，2010.

[14]Steger J L，Warming R F.Flux vector splitting of the inviscidgasdynamic equations with application to finite-differencemethods[J].Journal of Computational Physics，1981，40（2）：263-293.

[15]Dormand J R，Prince P J.A family of embedded Runge-Kutta formulae[J].Journal of Computational and AppliedMathematics，1980，6（1）：19-26.

[16]Jiang G S，Shu C W.Efficient implementation of weighted ENOschemes[J].Journal of Computational Physics，1996，126（1）：202-228.