不同飛行高度下超聲速來流/射流及其相互作用的數值模擬

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

(1．上海大學上海市應用數學和力學研究所，上海 200072;2. 中國科學院力學研究所高溫氣體動力學國家重點實驗室，北京 100190;3. 中國科學院大學工程科學學院，北京 100049)

后臺階流動、射流一直是流體力學關注的重要的基本流動現象[1-4]，而后臺階流動和射流的組合流動也普遍存在于流體力學的各個領域. 在火箭發射和飛行過程中，火箭處于不同的飛行高度，由于環境空氣密度隨著高度的增加而變小，火箭周圍的氣體靜壓也隨著高度的增加而變小，發動機尾噴管產生的射流形態也發生變化: 由地面附近的過膨脹形態轉變為高空時的欠膨脹形態[5-8]. 火箭噴管附近流場存在著復雜的超聲速來流和燃料射流相互作用情況，可通過兩步后臺階模型模擬超聲速來流與射流相互作用進行研究. 二者相互作用產生的回流區，會對射流形態產生影響，并影響到火箭發動機噴管的防熱[9-11]. 氣動噴嘴產生的燃料氣體與噴嘴周圍的空氣摻混[12]、超燃沖壓發動機中支板產生燃料射流和發動機內流相互作用[13]等現象也可視為兩步后臺階處超聲速來流和射流相互作用的流動現象. 因此超聲速來流與射流在兩步后臺階處相互作用的現象具有重要的科研價值和工程意義.

已有的工作主要運用后臺階研究流動產生流體分離和再附、后臺階回流區內渦結構等現象[14-16]: 超聲速來流在臺階拐角處由于外形的突變，產生流動分離現象; 氣體快速膨脹并在臺階下游位置產生一個剪切層，剪切層和后臺階下壁面之間形成回流區; 回流區內存在著復雜的渦結構. Biswas等[17]研究了膨脹比H/h= 1.942 3，2.5 和3.0 的后臺階流動，以及在流動方向上來流的空間演化和三維效應. 結果發現，在靠近內側壁處的剪切層存在向回流區噴射的現象. Hasan等[18]通過實驗研究了Re= 11 000 時層流分離的后臺階流動. 結果表明，與未受擾動的流動相比，擾動增加了剪切層的生長速率和湍流強度，并減少了再附長度. 擾動頻率及其諧波幅度的測量結果也表明，剪切層存在失穩. Reddeppa等[19]在高超聲速激波風洞中研究了Ma=7.6 時的后臺階流動. 結果發現，后臺階表面傳熱率會受到臺階高度的影響.

射流普遍存在于燃燒反應過程、發動機噴流等研究領域. Zapryagaev等[20]通過實驗確定，從軸對稱噴嘴排出的超聲速射流剪切層的初始部分具有三維結構，欠膨脹射流初始段是由膨脹扇、剪切層、攔截激波等基本結構構成的. Imamoglu等[21]運用大渦模擬和加權本質無振蕩(weighted essentially non-oscillatory，WENO)格式求解Navier-Stokes 方程組，發現模擬射流激波胞格結構良好地符合實驗結果. 劉昕等[22]研究了不同射流與環境靜壓比下，欠膨脹超聲速射流近場的流場結構. 結果發現，在欠膨脹超聲速射流情況下，欠膨脹射流方向出現的擬周期性攔截激波是欠膨脹射流穩定發展的特征，而這種擬周期性欠膨脹激波結構的消失是射流失穩開始的標志.

本工作分別研究了不同飛行高度下超聲速后臺階流動、射流以及超聲速來流/射流相互作用的流場結構和流動規律. 采用高精度格式離散Navier-Stokes 方程，通過改變不同來流條件參數，重點對流場的回流區渦結構、射流膨脹扇、兩個剪切層的相互作用等展開研究和討論.

1 物理模型和數值方法

1.1 物理模型

本工作采用兩步后臺階作為基本模型，研究了超聲速后臺階流動、射流和超聲速來流/射流組合流動這3 種流動形式. 超聲速后臺階流動的物理模型如圖1(a)所示: 從邊界AB 處產生Ma= 3 的超聲速來流，來流經過點C 時發生分離，產生來流剪切層. 射流的物理模型如圖1(b)所示: EF 處噴出Ma= 2 的高溫高壓射流，射流邊界附近形成射流剪切層. 超聲速來流/射流組合流動的物理模型如圖1(c)所示: 邊界AB 處產生Ma= 3 的超聲速來流，氣流從點C 處膨脹分離，形成一個來流剪切層; EF 處產生Ma=2 的射流. 因為射流出口處的靜壓顯著高于環境靜壓，因此此時射流的類型為欠膨脹射流，在射流邊界附近會出現射流剪切層. 來流剪切層和射流剪切層在第二個臺階上方相互作用，形成復雜的流場結構. 考慮到流場的對稱性，本工作流場的計算區域、CUP 和網格分布如圖2 所示，其中CUP 總數為72，網格總數為180 萬，采用矩形結構網格. 由于網格太密，圖2(c)為稀疏后(總數為1.8 萬)的網格分布. 壁面BC，CD，DE 為絕熱壁，HG 采用外推的無窮遠邊界條件，AH 采用超聲速無反射邊界條件，FG 采用對稱邊界條件，AB，EF 相應的邊界條件如表1 所示.

表1 邊界條件Table 1 Boundary conditions

圖1 物理模型Fig.1 Physical models

圖2 計算域、CUP 和網格分布Fig.2 Computational domain，CPU and grid distributions

當超聲速來流的Ma= 3 時，選取了3 種不同飛行高度的國際標準大氣(international standard atmosphere，ISA)條件作為來流的參數條件，選取的高度分別是20，30 和50 km，具體參數如表2 所示. 后臺階EF 處射流的Ma=2，參數設置如表3 所示，組分也為空氣.

表2 20～50 km 飛行高度的高空大氣標準參數Table 2 Atmospheric standard parameters at 20～50 km flight altitudes

表3 射流參數Table 3 Jet flow parameters

1.2 控制方程

在笛卡爾坐標系下，二維無量綱化Navier-Stokes 方程表示為

式中: 下標v表示粘性項; 各矢量形式為

式(2)～(6)中的ρ和p分別表示氣體的密度和壓力;u，v分別表示x，y方向上的速度分量; 各黏性應力分量為

各方向上的熱流分量為

其中黏性系數μ根據Sutherland 公式確定，Pr=0.72.

控制方程需要通過變換關系實現物理空間(x，y)向計算空間(ξ，η)的變換. 直角坐標系下的控制方程經過Jacobian 變換之后，計算空間中的Navier-Stokes 方程表示為

式中: 各矢量形式為

Jacobian 行列式為

取參考量L為特征長度，u0=a0，即來流速度為特征速度，溫度T0、分子黏性系數μ0、密度ρ0為特征變量. 對應的特征時間t0=L/u0，特征壓強，各無量綱量分別為

1.3 計算方法

本工作分別應用五階精度WENO 格式[24]、六階精度中心差分格式離散對流項和粘性項[25]，時間推進應用三階精度Runge-Kutta 格式[26]. 采用MPI 非阻塞并行模式，網格總數180 萬，使用72CPU 在廣州天河二號上運行.

Steger-Warming 流通量矢量分裂法通過對特征值進行分裂，表示為

五階精度WENO 格式對流項離散形式[24]如下:

式中，CONV 表示對流項.

矢通量構造為

式中:

其中IS 為光滑度量函數;n為一個整數. 本工作根據Jiang等[24]的建議，取n=2.

組合系數如下:

G通量的構造方式類似，替換F為G即可獲得表達式.

粘性項半離散逼近式[25]如下:

式中，六階精度中心差分格式導數項可以離散表示為

時間積分采用三階精度Runge-Kutta 方法[26]，即

1.4 數值驗證和網格無關化驗證

使用平板超聲速流動的算例來驗證計算方法的正確性和可靠性. 超聲速來流的Ma= 7，雷諾數Re=1.233×106. 如圖3 所示，在前緣激波之后，超聲速來流密度變小，并且邊界層沿著超聲速來流方向逐漸變厚. 這與平板超聲速流動的實驗情況一致.

圖3 平板上超聲速流動密度分布圖Fig.3 Density distributions of the supersonic flow on a plate

采用可壓縮邊界層的自相似解用于數值比較. 可壓縮邊界層的自相似解是半解析解. 本工作采用五階Runge-Kutta 方法結合打靶法求解自相似解. 將半解析解與數值模擬解進行比較，來驗證數值模擬的代碼和代數方法. 如圖4 所示，數值模擬曲線與半解析解曲線吻合良好，表明該算法是可靠的.

圖4 平板上超聲速流動邊界層上的速度分布Fig.4 Velocity distributions of the supersonic flow boundary layer on a plate

采用飛行高度為20 km 的后臺階流動和射流相互作用流場進行網格無關化驗證，以確定流場合適的計算網格. 選擇y= 0 處的數據來驗證網格無關性. 選擇t= 0.015 s 對稱線位置的靜壓曲線，使用了3 種不同總數的網格(180 萬，360 萬和540 萬). 圖5 顯示了不同網格的壓力分布. 結果發現，3 種網格的結果吻合良好. 本工作選擇總數為180 萬的網格進行數值模擬.

圖5 3 種網格的壓力曲線變化Fig.5 Changes of the pressure curves under three kinds of grids

2 數值結果與分析

2.1 不同飛行高度下的超聲速后臺階流動

在超聲速后臺階流動中，超聲速來流從兩步后臺階兩側流過(見圖6). 可以看出: 超聲速來流經過后臺階上方壁面時，由于壁面附近氣體速度梯度急劇變化，沿來流方向壁面起始處產生了微弱的前緣激波; 氣體經過拐點時，幾何外形發生突變產生分離形成剪切層，剪切層和臺階之間形成一個回流區; 回流區內氣流速度小，動能轉化為內能，溫度較其他區域明顯增大; 氣體在經過臺階拐角時發生膨脹，氣體靜壓急劇變小且流速大小和方向均發生改變，拐角附近形成膨脹扇區. 由于剪切層/激波作用，回流區內形成了渦結構，剪切層也因此擾動形成了激波(見圖7). 圖6 中3 種飛行高度下的流場結構類似，但當飛行高度為50 km 時，因為飛行高度的增加，密度和來流靜壓會有量級的降低，Re也會隨之有量級的減小，尾渦有明顯拉長的現象，剪切層的水平傾角也變小(見圖6(c)).

圖6 不同飛行高度下的后臺階流動密度梯度Fig.6 Supersonic flow acts on the double backward-facing steps at different flight altitudes

圖7 20 km 飛行高度下的后臺階流動底部渦結構Fig.7 Vortex structure at the bottom of the double backward-facing step at 20 km flight altitude

回流區內氣體流動速度較慢，并且有渦的相互作用，來流動能在回流區內轉化為內能，壁面溫度和來流溫度產生了明顯的差異，壁面溫度相對于來流溫度要高. 3 種飛行高度下的壁面溫度接近，后臺階壁面BC，CD，DE 的溫度為特征溫度的2.3～2.7 倍(見表4).

表4 不同飛行高度下的后臺階流動的壁面溫度Table 4 Wall temperatures of the double backward-facing step at different flight altitudes

2.2 不同飛行高度下的射流

射流從后臺階底部處產生后，首先經過一個三角形的等速核心區. 等速核心區內的溫度、密度等都為常數且不會改變. 等速核心區外產生膨脹扇區，膨脹扇區沿著射流的流動方向逐漸變寬. 膨脹扇區的發展造成了射流寬度的增大和三角形等速核心區的線性收縮，如圖8 所示.射流氣體的膨脹使得氣體加速并且速度方向發生改變. 同時，由于氣體的膨脹引起射流邊界附近的靜壓低于環境靜壓，使得射流氣體在剪切層邊界被壓縮，形成反射激波和射流剪切層. 在不同飛行高度下，流場基本結構類似.

圖9 為距離射流噴口0.05 處y方向半截面上的速度分布. 可以看出，3 種飛行高度的情況都屬于高度欠膨脹射流形態. 流體首先經過等速核心區，在等速核心區內流體速度恒定，均約為無量綱速度2. 隨后流體開始快速膨脹，形成對稱的膨脹扇區. 射流在經過膨脹激波時加速并且速度方向發生改變，此時速度開始增加. 膨脹扇區外因為流體靜壓與周圍環境靜壓不匹配，產生了攔截激波和剪切層，經過高強度的攔截激波和寬度較大的射流剪切層時，氣體速度急劇下降.

圖9 不同飛行高度下距射流噴口0.05 處半截面上的速度分布Fig.9 Velocity distributions on half section of jet flow at different flight altitudes(0.05 away from the nozzle)

不同飛行高度下膨脹扇區的大小有很大差異，尤其是在飛行高度為50 km 時氣體的密度、靜壓與其他兩個飛行高度相比有量級的降低，射流和環境氣體的靜壓比也明顯增加，射流邊界張角增大，膨脹扇區大小也隨之變大. 因此當飛行高度為50 km 時，速度在無量綱范圍0.05＜y ＜0.3 之間持續增加，而其他兩種高度下的速度增加范圍約為0.05＜y ＜0.1(見圖8). 該范圍反映了膨脹扇區的大小. 3 種飛行高度下氣體所能達到的最大速度均為射流初始速度的1.25～1.3 倍左右，且隨著飛行高度的增加，最大速度也增加. 飛行高度為50 km 時對應的最大速度明顯大于其他兩種飛行高度.

圖8 不同飛行高度下的射流密度梯度Fig.8 Density steps of jet flow at different flight altitudes

隨著飛行高度的增加，由于射流噴口靜壓與環境靜壓的比值顯著增大，射流的強度顯著增強，射流邊界附近的剪切層水平傾角也隨之增大，各飛行高度下的剪切層水平傾角如表5 所示，其中飛行高度為50 km 時的剪切層水平傾角最大，可達53.2°左右.

表5 不同飛行高度下的射流剪切層角度Table 5 Angles of jet flow shear layer at different flight altitudes

2.3 超聲速來流/射流在后臺階的相互作用

以飛行高度為20 km 時超聲速后臺階流動和射流相互作用的流場結構為例. 后臺階拐角處氣體發生分離產生來流剪切層，來流剪切層和后臺階之間形成回流區. 回流區內復雜的渦結構造成了來流剪切層的擾動，并在剪切層上形成小激波(如圖10(a)). 在射流起始的三角核心區外氣體發生膨脹形成膨脹扇區，氣體在膨脹扇區內的速度大小和方向均發生改變. 來流剪切層和射流剪切層交匯后發生相互作用，兩個剪切層匯合為一道混合剪切層. 混合剪切層兩側形成了兩個新的激波，兩側速度大小不同引發不穩定性，使得混合剪切層在后期失穩. 混合剪切層的失穩導致剪切層上產生了一些小激波結構，這些小激波結構被包裹在混合剪切層兩側的激波內. 如圖10 所示，不同飛行高度下流場結構類似，射流的形態有差異，飛行高度越高，射流剪切層的水平傾角越大，膨脹扇區面積越大.

圖10 不同飛行高度下超聲速來流/射流相互作用的流場密度梯度Fig.10 Density steps of supersonic flow/jet flow interaction at different flight altitudes

兩步后臺階兩側的超聲速來流會對射流本身形態產生影響，使得射流提前出現失穩現象:超聲速來流剪切層和射流剪切層發生相互作用產生混合剪切層，混合剪切層由于受到兩側來流和射流的雙重影響容易產生失穩. 以往研究表明，射流軸線壓力分布經過了典型的膨脹-壓縮-膨脹的循環過程[22]. 以飛行高度為20 km 為例，將只有射流的工況和超聲速來流/射流相互作用的工況下的射流軸線壓力進行對比發現: 在超聲速來流/射流相互作用的工況下，因為有超聲速來流的影響，流場結構提前出現了膨脹-壓縮-膨脹的射流結構，軸線壓力出現了明顯的波動現象，符合膨脹-壓縮-膨脹的現象(見圖11). 而在只有射流的工況下，由于計算域較小，射流并未出現膨脹-壓縮-膨脹的現象.

圖11 20 km 飛行高度下只有射流和超聲速來流/射流相互作用時的工況軸線速度分布Fig.11 Axial velocity distributions of only jet flow and supersonic flow/jet flow interaction at 20 km flight altitudes

當飛行高度為20，30 和50 km 時，選取超聲速來流/射流相互作用和只有射流時距離射流噴口0.05 處半截面上的速度分布，結果如圖12 所示. 首先射流氣體經過三角形等速核心區，且在該區域內氣體速度恒定. 隨后流體開始快速膨脹，氣體速度開始增加，當飛行高度為50 km 時，兩種工況下的膨脹扇區大小都明顯大于其他飛行高度，因此速度增加時的y范圍明顯大于其他飛行高度. 但由于超聲速來流/射流的相互作用，超聲速來流壓縮射流，導致膨脹扇區小于只有射流的工況，尤其當飛行高度為50 km 時，超聲速來流/射流的靜壓明顯小于其他飛行高度. 來流對射流的影響顯著，膨脹扇區受到來流壓縮的影響明顯變小，超聲速來流/射流相互作用時速度增大的y范圍明顯小于只有射流時的速度增大范圍.膨脹扇區外存在攔截激波和射流剪切層，氣流經過攔截激波后速度急劇下降. 可以觀察到，在飛行高度為50 km 的超聲速來流/射流相互作用工況下，在穿過剪切層后速度會有波動. 這是因為已經到達了回流區，剪切層擾動導致回流區內產生大量的渦結構，從而造成速度的波動現象.

如圖12 所示，超聲速來流對射流的壓縮，不僅會改變膨脹扇區的大小，也會對膨脹扇區內的最大速度產生影響(見圖12). 在相同的飛行高度下，膨脹扇區內所能達到的最大速度，只有在射流的工況下略大. 在超聲速來流/射流相互作用和只有射流兩種工況下，同一飛行高度下的速度在等速階段以后變大，兩種工況下速度曲線基本吻合. 這說明膨脹扇區內氣體膨脹的速率基本一致，有無超聲速來流對于射流膨脹扇區的膨脹速率沒有太大影響.

圖12 不同飛行高度下只有射流和超聲速來流/射流相互作用情況時距射流噴口0.05 處半截面上的速度分布Fig.12 Velocity distributions on half section of only jet flow and supersonic flow/jet flow interaction at different flight altitudes (0.05 away from the nozzle)

超聲速來流對射流的影響主要有三個方面: ①超聲速來流會壓低射流邊界和射流剪切層，使得射流膨脹扇區縮小，膨脹區內速度最大值也會變小，但對射流膨脹速率影響甚微; ②超聲速來流和射流交匯后，兩個剪切層會發生相互作用并匯合為混合剪切層，因混合剪切層兩側的速度不同而產生不穩定性，造成混合剪切層的失穩，形成一些小激波結構，改變了射流本身結構; ③混合剪切層兩側由于氣體的壓縮產生兩個激波，兩個激波的出現，使得射流內部的膨脹-壓縮-膨脹現象提前出現.

射流剪切層的角度變化可以反映超聲速來流對射流的壓縮程度. 高度越高，來流靜壓越小，射流的影響相對變強，射流剪切層的水平傾角隨著飛行高度的增加而增大. 如表6 所示，同一飛行高度下，超聲速來流/射流相互作用時的射流剪切層水平傾角小于只有射流的情況. 當飛行高度為20 km 時，射流靜壓約為來流靜壓的3 倍左右，此時超聲速來流的影響較大. 在超聲速來流/射流相互作用的工況下，射流剪切層的水平傾角約為7°左右，而該飛行高度下只有射流工況下的射流剪切層水平傾角高達40°左右. 隨著飛行高度增加，來流對射流的壓縮效應變弱，剪切層的水平傾角也隨之增大. 尤其當飛行高度為50 km 時，射流靜壓約為來流靜壓的250 倍左右，此時來流影響較小，超聲速來流/靜流相互作用后射流剪切層的水平傾角為41°，相比于飛行高度為20 km 時明顯增大.

表6 不同飛行高度下的射流剪切層角度Table 6 Angles of jet shear layer at different flight altitudes

射流剪切層水平傾角變化對于流場結構產生的影響如圖13 所示. 可以看出: 飛行高度增加，超聲速來流的影響變小，射流影響相對變大，射流剪切層水平傾角變大，底部渦區會有明顯的拉升. 對于不同飛行高度，由于射流剪切層和超聲速來流剪切層水平傾角的變化，導致回流區大小和形態發生改變，進而影響回流區大尺度渦的數量. 在飛行高度為20 km 時，超聲速來流靜壓要大于其他飛行高度，對流場結構的影響較大，此時超聲速來流剪切層向下發展. 當飛行高度為30 km 時，超聲速來流靜壓相對于20 km 時有所減小，對射流的影響變小，此時超聲速來流剪切層抬升，基本呈水平狀態. 上述兩個飛行高度中回流區內大尺度渦結構豐富. 在飛行高度為50 km 時，超聲速來流靜壓相對小于其他高度，射流在超聲速來流/射流相互作用中占主導地位，射流剪切層水平傾角變大，進而抬高了來流剪切層，回流區明顯變大，此時回流區內只觀察到兩個大尺度渦.

圖13 不同飛行高度下超聲速來流/射流相互作用時的底部渦結構Fig.13 Vortex structures of supersonic flow/jet flow interaction at different flight altitudes

超聲速來流/射流相互作用會對回流區壁面溫度產生影響. 與只有后臺階流動的工況相比，射流的存在使得回流區內渦結構相互作用更加劇烈，因而壁面溫度也有明顯提高. 選取兩步后臺階回流區內BC，CD 兩個壁面. BC 段溫度如表7 所示，可見不同飛行高度下壁溫約為來流溫度的2.8～4.3 倍左右，且高度越高、溫度越高，加熱的效果越明顯. CD 段溫度如表8 所示，約為來流溫度的4.0～5.3 倍左右. 與只有后臺階流動時的溫度(見表5)相比，超聲速來流/射流相互作用下BC，CD 的溫度比只有后臺階流動時大得多.

表7 不同飛行高度下超聲速來流/射流相互作用時的BC 段溫度Table 7 Temperatures of BC section of supersonic flow/jet flow interaction at different flight altitudes

表8 不同飛行高度下超聲速來流/射流共同作用時的CD 段溫度Table 8 Temperatures of CD section of supersonic flow/jet flow interaction at different flight altitudes

3 結束語

本工作分別應用五階精度WENO 格式、六階精度中心差分格式離散對流項和粘性項，時間積分采用三階精度Runge-Kutta 格式，研究了超聲速后臺階流動、射流的基本流場結構，并進一步研究了超聲速來流/射流在后臺階的耦合流動，得到如下研究結果.

(1) 在后臺階流動的工況中，氣體經過后臺階拐角形成一個剪切層，剪切層后臺階之間的回流區渦造成了剪切層的擾動，使剪切層上方形成了激波. 當飛行高度為50 km 時，剪切層明顯拉長.

(2) 在只有射流的工況中，射流首先經過三角形等速核心區，等速核心區外氣體膨脹形成膨脹扇區，膨脹扇區外形成攔截激波和射流剪切層. 射流剪切層的水平傾角隨著飛行高度的增加有小幅增大.

(3) 在超聲速來流/射流相互作用的工況中，射流剪切層和超聲速來流剪切層交匯并發生相互作用，形成一個混合剪切層. 混合剪切層兩側形成兩道激波，混合剪切層上由于失穩產生小激波結構. 隨著飛行高度的增加，射流剪切層明顯增大.

(4) 在超聲速來流/射流相互作用的工況中，有可能將射流氣體卷入了兩步后臺階的回流區中，對壁面形成顯著的加熱. 對兩個剪切層相互作用的流場機理尚不成熟，需要在以后的工作中進一步研究.