超聲速高度欠膨脹沖擊射流的大渦模擬

劉 海，姚朝暉

（清華大學 航天航空學院，北京100084）

0 引 言

欠膨脹超聲速沖擊射流在許多工程技術領域得到廣泛的應用，如STOVL（Short Take Off and Vertical Landing）飛行器起飛和著陸、火箭的發射，以及除塵、除水等工程問題；而且其流場結構復雜，激波和渦之間存在著強烈的相互作用。因此無論是工程應用還是基礎研究，沖擊射流都極具研究價值［1－3］。

Alvi［1］等通過實驗給出了不同工況下的流場陰影圖和壁面上的壓強分布曲線，并指出沖擊平板上的壓強脈動的主頻與沖擊單音的頻率相同。姚朝暉等［3］利用紋影和高速攝影技術發現在馬赫盤出現之前，沖擊射流的流場是整體振蕩的；而馬赫盤出現之后，流場振蕩局限于馬赫盤與沖擊平板之間。Henderson［4］討論了在高度欠膨脹工況下流場的沖擊單音的存在與馬赫盤出現位置的關系，觀察到了流場的環形激波的生成和消失，給出一些瞬時脈動流場的渦的信息，并對反饋環進行了新的論述。

圖1為高度欠膨脹（來流總壓p0跟環境壓強pa的比值NPR大于3.8）沖擊射流的流場示意圖。由于RANS只能得到流場的平均場的信息，而相比實驗LES可以獲得更豐富的流場信息，所以本文使用大渦模擬對高度欠膨脹的超聲速沖擊射流的流場進行了數值模擬，主要是對流場的振蕩、內外剪切層之間的環形激波的生成和消失過程以及內剪切層的渦結構的演化進行了仔細的研究。

圖1 高度欠膨脹沖擊射流流場示意圖Fig.1 A schematic of the highly underexpanded supersonic impinging jet

1 數值模擬方法

首先對柱坐標系下的三維可壓縮非定常NS方程進行基于密度的加權過濾［5］：

其中，x、θ、r分別表示的是軸向、周向和徑向坐標，各項參數如下：

式中τ表示指分子粘性項，而σ和Q則代表過濾產生的亞格子應力和亞格子熱通量項。

由于方程無量綱化參考的是環境流場和噴嘴出口直徑，所以氣體的狀態方程為：

粘性系數的計算采用Sutherland公式：

對于空氣，C＝110.4／Ta，Ta為環境溫度。

本文的亞格子應力模式采用的是Smagorinsky模型：

式中Prt為湍流普朗特數，而渦粘系數的計算公式為：

式中Cs為Smagorinsky系數，在壁面附近再進行如下的修正：

計算中對流項采用五階WENO格式進行離散，粘性通量項采用顯式的4階中心格式進行離散。時間積分采用四步三階具有強穩定保持性質的Rung－Kutta格式［6］，其CFL數等于2，所以計算效率較高。并對計算程序進行了并行化。

本文在計算取噴嘴出口為臨界聲速來流條件；出口和遠場邊界取與時間相關的無反射邊界條件［7］；沖擊壁面和噴嘴壁面取絕熱無滑移固壁邊界條件；中心軸在柱坐標系下為奇軸，為保證速度在軸線上單值，采用Lele基于級數展開得到的公式進行插值［8］。

2 程序驗證

本文程序驗證所用的算例與Alvi的實驗工況［1］相同：其中收縮噴嘴的出口直徑為d＝25.4mm，沖擊距離h／d＝2.0，來流總壓跟環境壓強的比值NPR＝2.5，沖擊平板是大平板。計算網格數：x×θ×r＝300×48×250。本文程序有效性的驗證主要是將平均場的數據結果和實驗進行對比。

觀察圖3、圖4可以發現本程序計算得到的流場的激波柵格結構的大小和位置以及沖擊壁面上的壓強系數（Cp＝（p－pa）／（p0－pa））分布與實驗結果吻合較好。從而驗證了程序可以有效地模擬超聲速沖擊射流的流場。

3 流場的振蕩

下面討論的算例是基于Henderson［4］的實驗工況：噴嘴直徑d＝25.4mm，沖擊距離h／d＝2.08，噴嘴唇厚0.1d，壓比NPR＝4.03，沖擊平板為大平板。計算網格數為：x×θ×r＝300×48×250。此種工況流場為對稱模態，所以后文分析主要是基于θ＝0的平面進行的，圖中r／d＝0對應于計算域的中心軸線。

由圖5的馬赫數分布可以看出馬赫盤及斜激波向流場下游運動的時候，斜激波后面的流場的亞聲速區域逐漸增大，同時馬赫盤后面的流場的速度也在減小，而馬赫盤往上游運動時斜激波和馬赫盤后的流場變化趨勢相反。原因是當馬赫盤向下游運動的時候，馬赫盤前面的流場膨脹區域變大，所以氣體會進一步膨脹，馬赫盤前的馬赫數提高，所以馬赫盤的激波強度變大，導致馬赫盤后面的流場速度減小。同時隨著馬赫盤往下游運動，斜激波的激波角度變大，所以斜激波強度也變大。馬赫盤在向下游運動時，可以發現馬赫盤的直徑有減小的趨勢。

4 環形激波的生成與消失

圖6對計算和實驗的瞬時流場進行了簡單的對比，可以看出流場的馬赫盤，斜激波和環形激波的形狀基本一致。在Henderson的實驗陰影圖中，由于沖擊壁面附近的圖像不清晰，在內外剪切層之間只觀察到一道環形激波，而計算中觀察到兩道環形激波的生成與消失的周期過程。

圖6 數值計算的密度云圖和實驗陰影圖Fig.6 Computed density contour and measured shadowgraph photograph

圖7為密度等值線和旋轉強度云圖，由圖7（a）可見在馬赫盤后面有兩道環形激波a和b存在于流場的內外剪切層之間，在剪切層上存在著大尺度渦結構，此時馬赫盤往下游運動；圖7（b）中渦結構A、D追上環形激波a一起向下游運動，而此時激波a和激波b強度變弱；圖7（c）中可看出激波b已經消失，激波a下行過程中在靠近它的地方開始有新環形激波bb生成，環形激波a進一步變弱；在圖7（d）中環形激波a消失，在斜激波后面形成新的環形激波aa，而激波bb隨大尺度渦結構A、D在往下游運動的過程中逐漸變強。

斜激波后第一道環形激波形成后在內外剪切層的大尺度渦結構追上它之前強度逐漸變強，而兩個渦結構追上它一起下行過程中則逐漸變弱，消失前在靠近它的附近生成一道新的環形激波，新的環形激波會隨大尺度渦下行過程中先增強然后變弱消失。同時可以發現環形激波生成與消失的周期、大尺度渦結構周期和馬赫盤的振蕩周期一致。

5 內剪切層的渦

Henderson在實驗結果中給出了流場的脈動場，此處本文對流場的內剪切層的渦結構的演化進行進一步的討論。

由于流場運動的周期性，所以首先看圖8的最后兩個圖8（d）和（e）然后接著再看圖（a），馬赫盤由向右運動變成了向左運動，相比其它時刻的圖，可見在馬赫盤與斜激波的交點處，圖8（e）中一個相對較大的渦結構開始卷起。即在馬赫盤與斜激波相交的地方，當馬赫盤運動到最下游時，內剪切層上一個相對較大的渦結構開始卷起，然后在下行過程中發展成大尺度渦結構，如圖8中的渦B。

圖8 內剪切層的渦結構演化（dt＝0.225）Fig.8 The evolution of the vortex（dt＝0.225）

對于圖8中的渦結構A，與第二道環形激波一起下行，在到達沖擊壁面時與壁面作用然后破碎，所有這些破碎的渦結構我們用A′表示，其中破碎的渦結構大部分隨壁射流沿徑向流出，小部分向回流區運動。由于渦A破碎成A′時，此時內外剪切層間的第二道環形激波靠近沖擊壁面，在環形激波和沖擊壁面之間形成高壓區，使得破碎的渦結構沿射流壁面的徑向運動受阻，會在原地打轉一段時間，其中一部分破碎渦跟隨回流往中心區移動。隨著時間的發展，靠近壁面的環形激波會慢慢變弱消失，高壓區消失，這些破碎的渦結構大部分將跟隨著壁射流沿徑向運動。

在圖8的壁射流區，可以觀察到在內外剪切層的作用下形成了壁射流區內外交錯的渦結構序列。

6 結 論

基于柱坐標系下的NS方程，利用大渦模擬的方法創建了求解超聲速沖擊射流的三維并行計算程序，驗證了程序的有效性。

對所模擬的高度欠膨脹沖擊射流算例，在斜激波、馬赫盤及渦結構的共同作用下，內外剪切層之間的環形激波存在著周期性的生成和消失過程。環形激波的周期過程與大尺度渦運動、馬赫盤振蕩周期一致。

內剪切層上的大尺度渦結構的形成與馬赫盤和斜激波的振蕩相關；環形激波與沖擊壁面之間的高壓會對內剪切層的渦沿徑向運動形成一定的阻礙，內外剪切層的作用形成了壁射流區內外交錯的渦結構。

［1］ALVI F S，LADD J A.Experimental and numerical investigation of supersonic impinging jets［R］.AIAA Paper 2000－2224，2000.

［2］ALVI F S，IYER K G.Mean ＆unsteady flow properties of supersonic impinging jets with lift plates［R］.AIAA Paper 99－1829，1999.

［3］姚朝暉，何楓，韓標，等.高速沖擊射流流場特性與噪聲機理的研究［J］.流體力學實驗與測量，2003，17（2）：84－87.

［4］HENDERSON B，BRIDGES J，WERNET M.Experimental study of the oscillatory flow structure of tone producing supersonic impinging jets［J］.J.Fluid Mech.，2005，542：115－137.

［5］張兆順，崔桂香，許春曉.湍流大渦數值模擬的理論與應用［M］.北京：清華大學出版社，2009：120－134.

［6］RUUTH S J，SPITERI R J.High－order strong－stability－preserving Runge－Kutta methods with downwind－biased spatial discretizations［J］.SIAMJournalofNumericalAnalysis，2004，42（3）：974－996.

［7］張涵信，沈孟育.計算流體力學－差分方法的原理和應用［M］.北京：國防工業出版社，2003：313－341.

［8］MORINISHI Y，VASILYEV OV，OGI T.Fully conservative finite difference scheme in cylindrical coordinates for incompressible flow simulations［J］.Journal ofComputationalPhysics，2004，197：686－71.