軸對稱后向臺階不穩定流動及壓強振蕩數值研究①

張翔宇，何國強，劉佩進

(西北工業大學燃燒、熱結構與內流場重點實驗室，西安 710072)

0 引言

后向臺階流動具有幾何結構簡單、分離點固定和來流條件容易控制等特點，是研究湍流邊界層分離和再附流動的重要形式之一。目前，國內外對層流至湍流條件下的流動進行了大量的試驗和數值模擬［1－3］，試驗大多在風洞和水洞中進行，測量手段有熱線、熱膜探針和PIV等。目前，為提高戰術固體火箭發動機性能，大多采用高能復合推進劑及復雜裝藥構型，尤其是翼柱形裝藥的固體火箭發動機在翼槽燒完時，也會形成具有自身特點的后向臺階流動。這種臺階結構是軸對稱的，具有一定傾斜角度，并且由于側向加質的作用，近壁面速度垂直于軸線，流場中主流速度較大，一般為完全發展的湍流流動。前柱后翼形裝藥形成的后向臺階會使流動在臺階附近形成剪切層，進而產生轉角渦脫落現象，這是固體火箭發動機中3種典型的旋渦脫落現象之一［4－5］。當渦脫落頻率與聲場某階固有頻率接近時，就可能形成共振，這種聲渦耦合產生的共振是大長徑比固體火箭發動機燃燒室壓強振蕩的源頭。Flandro和Jacobs［6］最早指出，流動經過障礙物時，形成的旋渦脫落會將部分平均流能量轉化為聲能，從而產生聲能增益。Mason 和 Folkman［7］以及 Brown［8］等得出結論，旋渦脫落可誘發燃燒室的軸向基頻。Culick 和 Magiawala［9］、Dunlap［10］通過實驗研究表明，當旋渦脫落頻率與燃燒室聲場固有頻率接近時，會引起壓強振蕩現象。Vuillot［11］、Lupoglazoff［12］設計了一種轉角渦脫落致壓強振蕩的發動機，并通過試驗及數值模擬證明了轉角渦脫落產生的聲渦作用，可引起與聲場固有頻率相耦合的壓強振蕩現象。Shanbhogue［13］等通過冷流試驗，研究了轉角渦脫落引起的壓強振蕩問題，得到了St數隨流場結構參數變化的規律。

本文對戰術發動機內形成的軸對稱后向臺階流動進行模型化處理，設計了一種無需外加擾動的旋渦脫落致壓強振蕩的流場結構。運用穩定性理論得到聲渦耦合模態預測方法，并采用大渦模擬(LES)方法，針對寬來流條件進行數值模擬。通過對比計算結果，驗證了模態預測方法的準確性，分析了轉角渦脫落產生的一般過程及旋渦多尺度運動規律。通過分析流場關鍵位置的壓強振蕩，研究了氣動聲學特性及旋渦運動與壓強振蕩的關系。

1 計算構型設計

1.1 聲渦耦合模態預測方法

對于某一特定構型的發動機燃燒室，如果其平均流邊界層的不穩定頻率處于軸向聲模態范圍內，就有可能形成旋渦脫落引起的壓強振蕩現象。Lupoglazoff等通過求解不可壓Orr－Sommerfeld方程，得到了后向臺階結構燃燒室的徑向速度分布曲線。分析表明，具有拐點的雙曲正切形速度分布是剪切層發展為離散旋渦的重要條件。應用空間增長率理論，分析流體動力學穩定性問題，旋渦波動量為

其中，α =αr+iαi∈C，ω∈R，B 是積分常數，f(y)、φ(y)、π(y)是渦量波動分布。任意給定ω，計算求解本征問題的解為α=αr+iαi，并引入以下無量綱角頻r=2πfδ/ΔUr+iai。

通過處理可得到cr與實際頻率的關系曲線。分析表明，在cr的臨界值(cr)cr以上，邊界層是穩定的。此外，計算中還需定義以下參數:無量綱旋渦位移速度將邊界層厚度表示成燃燒室直徑D的函數δ=KDD;將剪切層速度差ΔUe表示成平均流速度U的函數ΔUe=KUU=KU(a0M)。最終形成后向臺階不穩定流動需同時滿足以下2個條件:

(1)流場第p階固有模態fp位于剪切層不穩定范圍內，即無量綱頻率(cr)p小于臨界頻率(cr)cr;

本文采用半封閉聲腔模型fp=(2p－1)a0/(4L)。

(2)旋渦運動距離l應至少包含一個渦波長度λv，即 λv＜l。

代入KV、KU表達式，對于給定構型，可認為二者為常數。聯立以上二式，整理得:

式中的關鍵參數參考文獻［11］中的結果，代入數值計算得到:

式(4)就是后向臺階流動中，給定長徑比及臺階位置條件下主流速度與轉角渦脫落可能誘發的聲模態之間的關系。需指出的是，上述結論僅是發生不穩定的必要條件。盡管如此，經驗表明，輕微的不穩定趨勢，就有可能引發聲渦耦合現象。這也就是為什么當滿足必要條件時，就很容易發生不穩定，而制造完全穩定的發動機卻十分困難。

1.2 研究對象的幾何結構

為得到后向臺階流動渦脫落及氣動聲學的一般規律，需對計算模型進行適當簡化。忽略側向加質的影響，將后向臺階擴張角θ確定為90°。長徑比參考戰術發動機燃燒室特征，并進行適當縮比，整個流動區域為軸對稱結構，總長1 200 mm;將后向臺階布置在流場中央，臺階前后段L1、L2等長，且均為600 mm;單側后向臺階高度H為25 mm;臺階前后管路直徑D1、D2分別為50 mm和100 mm。流動區域示意圖見圖1。

2 數值研究方法

在后向臺階流動中，臺階后緣會產生類似平板剪切層的大渦序列，整個流場中的旋渦結構十分復雜。工程應用中，普遍采用的雷諾時均(RANS)方法難以模擬這種復雜的流動特性和旋渦結構。因此，本文采用大渦模擬(LES)方法，對多尺度旋渦的運動過程進行數值模擬。

2.1 控制方程離散及求解

LES在湍流數值模擬中只計算大尺度脈動，將小尺度脈動對大尺度運動的作用建立亞格子模型，過濾小尺度脈動，并得到亞格子應力項，亞格子模型的選擇至關重要。后向臺階流動中，剪切層與近壁面流動關系密切。因此，選用對近壁面區域進行修正的WALE(Wall－Adatping Local Eddy－viscosity)亞格子模型，從而避免了Smagorinsky模型難以準確預測近壁區流動及湍流轉捩現象的不足［14］。計算中，不考慮化學反應以及凝相粒子的影響，控制方程采用純氣相可壓N－S方程?？臻g參數離散采用具有二階精度的BCD(Bounded Central Differencing)格式，時間上采用二階隱式格式，時間步長為10－5s。為提高計算的收斂性，先采用空間二階離散格式的RANS方法得到流場的穩態解，然后采用上述LES方法求解非定常流動。

2.2 計算網格與邊界條件

由于LES受網格密度及人為劃分的影響顯著，須將可能出現大尺度回流區與渦脫落區的計算域進行高精度的空間網格剖分。在臺階前后、轉角及壁面附近進行加密。計算采用二維軸對稱方法，通過網格收斂驗證，總網格數約為55 000(圖2)。介質為常溫常壓可壓氣體，粘性模型采用Sutherland定律，入口邊界給出空氣質量流率、來流總壓和初始溫度，出口邊界給出靜壓和回流溫度，壁面采用無滑移絕熱壁面條件。該計算模型是一個進口為聲學反射邊界的半封閉聲腔。

2.3 計算方法驗證

后向臺階流動數值模擬的關鍵是對旋渦運動及壓強振蕩的準確模擬。陳曉龍等［15］通過搭建冷流試驗系統，在高壓空氣中加入示蹤粒子，對轉角渦脫落進行了可視化研究。

圖3是應用上述數值方法對流場進行的二維大渦模擬結果與冷流試驗中的高速攝影結果。在激光照射下，混合了粒子的流向氣體呈亮白色，而照片中的“黑洞”現象是由于旋渦中高速運動氣體將粒子甩出造成的?？煽闯觯疚牟捎玫臄抵捣椒ǖ玫降男郎u位置和結構與試驗結果較為吻合。

此外，還對VKI實驗室針對Ariane 5 P230助推器縮比模型的冷流試驗進行了數值模擬校驗。該實驗器為軸對稱結構［16］(圖4)，通過高頻壓強傳感器測量了流場不同位置的壓強振蕩。圖5是計算域網格劃分，與后向臺階流動數值模擬采用的網格剖分尺寸一致。

發動機頭部壓強振蕩頻率的試驗結果和數值模擬結果見表1。由表1可看出，采用的二維軸對稱LES數值方法得到的壓強振蕩頻率與試驗結果較為接近，前四階頻率與實驗結果間的最大誤差小于6%，可較好地對湍流場氣動聲學問題進行研究。

表1 LES結果與試驗結果對比Table 1 Results of LES and cold experiments Hz

3 后向臺階流動的旋渦結構

3.1 轉角渦脫落的形成過程

在后向臺階流動中，一般以雷諾數作為特征參數進行研究，定義Re=UH/ν。其中，U為來流速度;H為臺階高度;ν為介質的運動粘度。為方便分析流場特征，需將流動區域進行劃分。

以Re=11 288(Ma=0.02)的算例為例，處理穩態計算結果繪制流函數分布圖(圖6)，在緊靠臺階后緣突擴區域形成一個低速角渦區，在其下游是一個狹長的近壁面回流區，接著流動再附著并繼續發展。通過上述分析，可得如圖7所示的后向臺階湍流分離流動流場的基本劃分方法。

連續介質經過后向臺階時，會在臺階下游形成典型的流動剪切層，隨著剪切層向下游發展，形成流向渦，這種大尺度旋渦運動對流動穩定性及氣動聲學問題都具有重要影響。

LES可獲得旋渦結構的運動演化過程，通過整個計算區域內的等渦量線圖可觀察剪切層的運動特征，包括旋渦的產生、脫落和合并過程。

圖8描述了流動穩定后轉角渦脫落形成的典型過程。圖8(a)是一個旋渦脫落后新旋渦產生的初始過程，隨著剪切層在后向臺階下游持續發展，由于流動自身的不穩定性，剪切層沿徑向擺動，同時受到近壁面回流區的卷吸作用，沿流向運動了一定距離后，剪切層向壁面方向發展(X=0.67 m)。

遠離臺階后，流動速度逐漸降低，向壁面移動的部分剪切層又受到后續發展的剪切層的擠壓作用，此時新的渦核基本形成(X=0.68 m)，如圖8(b)所示。

此后，隨著旋渦向下游繼續發展，新的渦量注入渦核，旋渦體積增加。當旋渦增長到一定體積時，完整的渦結構基本形成(X=0.69 m)，其受剪切層的影響也逐漸減小(圖8(c))。

此后，逐漸形成了穩定的渦脫落(X=0.70 m)，渦量呈離散式分布(圖8(d))。下文分析中，采用的旋渦脫落頻率按上述過程進行判定。

3.2 轉角渦脫落的形成過程

后向臺階流動的軸向速度分布特點是臺階下游流速低，且距離臺階越遠，流速越低。由于流道結構突擴及壁面的粘性作用，管流中心流速大于近壁面處流速。由于湍流的多尺度性和隨機性，旋渦運動也呈現類似特點。流場中旋渦運動的典型過程如圖9所示。

臺階后緣由剪切層發展出離散的小尺寸旋渦(X=0.64 m)，其在向下游發展過程中，運動速度逐漸降低，被后續形成的旋渦追趕，并發生初始合并過程(X=0.74 m)，隨著壁面渦量及離散渦量的注入(X=0.9 m)，進而形成大尺度旋渦向下游運動，直至離開流場(X=1.0 m)。

通過分析渦量云圖，將不同來流條件下的旋渦多尺度運動頻率列于表2。由于旋渦運動的隨機性，取旋渦運動3～5個周期的平均值?？煽闯觯S著來流速度增加，各尺度的旋渦運動頻率均相應增加。

表2 多尺度旋渦運動頻率Table 2 Frequency of the multi-scale vortex motion

4 壓強振蕩及氣動聲學特性

后向臺階流動中的旋渦脫落和旋渦合并過程，會引發流動的不穩定性。Powell方程指出［17］，在流速很低的條件下，僅存在渦的地方才會產生聲。當旋渦運動與聲場固有頻率接近時，就可能發生聲渦耦合現象，造成流場中持續的壓強振蕩。LES方法不但可模擬不同來流條件下的旋渦發展過程，還可得到流場各處壓強隨時間的變化規律。通過FFT分析，獲得壓強波動的特征頻率和振幅。計算中，記錄流場中關鍵部位的壓強值，壓強監測點分布如圖10所示。

圖10中，1位于流場進口壁面上，試驗中高頻振蕩壓強傳感器一般布置在燃燒室頭部，是研究燃燒不穩定問題的重要布測點;2位于臺階轉角處;3表示布置在角渦區和回流區的監測點;4表示布置在再發展區的監測點。

對于不規則的聲腔結構，ANSYS軟件的模態分析，可確定固有頻率及振型。試驗證實，該方法對頻率預測準確性很高。ANSYS分析聲場前5階軸向模態分別為 99.5、182.6、381.5、464.7、663.6 Hz。對各工況下數值模擬記錄的監測點1壓強值進行FFT分析，得到結果如圖11所示，橫坐標表示壓強振蕩頻率，縱坐標表示振幅(波動量與平均壓強之比)。將流場的結構及各工況來流參數代入式(4)進行計算，得到誘發聲場前5階軸向振型所需的來流Ma范圍(表3)。并將上述結果列于表4。

表3 聲渦耦合模態預測Table 3 Prediction of acoustic-flow coupling mode

表4 數值模擬與聲渦耦合模態預測結果對比Table 4 Results of the mode prediction and LES

對比ANSYS模態分析結果與數值模擬結果可知，本文設計的后向臺階流動出現了較為明顯的壓強振蕩現象。在設計來流范圍下，流動最多可誘發聲場的前5階軸向頻率。隨著來流速度增加，各階軸向振型相繼被激發，同時壓強振蕩振幅也逐漸增大，在流場頭部前兩階軸向頻率的振幅最大，最大振幅接近平均壓強的5‰。聲渦耦合模態預測方法對不同速度下可能誘發的軸向聲模態預測結果與數值模擬結果十分吻合，當來流速度接近激發某階固有頻率所需的速度下限時，該頻率逐漸被激發。分析可知，激發軸向基頻所需的來流速度最低。工況1對應的來流Ma=0.01，此時沒有發生聲渦耦合現象。因此，通過該方法，也可對后向臺階穩定流動的速度上限進行預估。

對比相同工況下流場中不同壓強監測點的FFT分析結果發現，監測點2、3與1的壓強振蕩頻率一致，但各階頻率的振幅更低，流場頭部壓強振蕩最為強烈。處于再發展區的監測點4各階頻率的振幅更小，并出現了小于軸向基頻的低振幅的寬峰振蕩，此頻率隨來流速度增加而逐漸增加。

工況2～5監測點4的FFT分析結果如圖12所示。結合表2給出的多尺度旋渦運動結果分析發現，低于軸向基頻的頻率與大尺度旋渦運動頻率接近。在來流速度較低的情況下(工況2)，再發展區中大尺度旋渦運動頻率對壓強振動的影響較大。隨來流速度增加，聲渦耦合現象逐漸增強，當聲場某階固有頻率被激發后，壓強振蕩頻率與聲場固有頻率相耦合，流場中的壓強振蕩以聲場軸向各階頻率為主。當大尺度旋渦運動頻率與基頻接近時，會造成基頻的多峰現象(工況4、5)。表2的結果也表明，當近臺階小尺寸旋渦形成頻率接近聲場某階固有頻率時，該頻率及更低階頻率就有可能被激發。這說明后向臺階流動中，轉角處旋渦脫落及旋渦運動是誘發聲渦耦合現象的源頭。

5 結論

(1)后向臺階下游形成剪切層，由于流動不穩定性和近壁面回流區的作用，剪切層會發展為渦量離散分布的旋渦結構。隨來流速度增加，旋渦脫落頻率增大，脫落位置也更靠近后向臺階。小尺寸旋渦被后續形成的旋渦追趕，并發生初始合并過程，隨壁面渦量及離散渦量的注入，最終形成大尺度旋渦，直至離開流場。

(2)從穩定性理論推導出后向臺階流動的聲渦耦合模態預測方法，可對不同速度下可能激發的聲場頻率進行合理預測。當來流速度增加至激發某階固有頻率所需的速度下限時，該頻率逐漸被激發。通過該方法，還可對后向臺階穩定流動的速度上限進行預估。

(3)流場頭部壓強波動振幅最大，頻率與聲場軸向模態相耦合。再發展區出現低于軸向基頻的寬峰振蕩受大尺度旋渦運動的影響，該頻率與基頻接近時，會造成基頻的多峰現象。當近壁面小尺寸旋渦形成頻率接近聲場某階固有頻率時，該頻率及更低階頻率就有可能被激發。

［1］Armaly B F，Durst F，Pereira J C F，et al.Experimental and theoretical investigation of backward－facing step flow［J］.Journal of Fluid Mechanics，1983，127:473－496.

［2］Lee T，Mateescu D.Experimental and numerical investigation of 2D backward facing step flow［J］.Journal of Fluids and Structures，1998，12:703－716.

［3］劉應征，樸英守，成亨鎮.后臺階分離再附湍流流動的實驗研究［J］.上海交通大學學報，2005，39(5):810－817.

［4］Fabignon Y，Dupays J.Instability and pressure oscillations in solid rocket motors［J］.Aerospace Science and Technology，2003，7:191－200.

［5］Telaro M，Avio Spa，Stella F，et al.Pressure oscillations in solid rocket motors［R］.AIAA 2006－4423.

［6］Flandro G A，Jacobs H R.Vortex－generated sound in cavities［R］.AIAA 1973－1014.

［7］Mason D R，Folkan S L，Behring M A.Thrust oscillations of the space shuttle solid rocket booster motor during static tests［R］.AIAA 79－1138.

［8］Brown R S，et al.Vortex shedding studies［R］.AIAA 80－1092.

［9］Culick F E C，Magiawala K.Excitation of acoustic modes in a chamber by vortex shedding［J］.Journal of Sound and Vibration，1979，64(3):455－475.

［10］Dunlap R，Brown R S.Exploratory experiment on acoustic oscillations driven by periodic vortex shedding［J］.AIAA Journal，1981，19(3):408－409.

［11］Vuillot F.Vortex－shedding phenomena in solid rocket motors［J］.Journal of Propulsion and Power，1995，11(4):626－639.

［12］Lupoglazoff N ，Vuillot F.Numerical simulation of vortex shedding phenomenon in 2D test case solid rocket motors［R］.AIAA 1993－3006.

［13］Shanbhogue S J，et al.Aeroacoustic instability in a solid rocket motor with finocyl grain［R］.AIAA 2003－4632－4636.

［14］張兆順，崔桂香，許春曉.湍流理論與模擬［M］.北京:清華大學出版社，2005.

［15］Chen Xiao－long，He Guo－qiang，Liu Pei－jin.Experimental investigation on pressure oscillation induced by flow instability in SRM［C］//60th IAC，Daejeon，Korea，2009.

［16］Anthoine J.Experimental and numerical investigations of nozzle geometry effect on the instabilities in solid－propellant boosters［R］.AIAA 2000－3560.

［17］Powell A.Theory of vortex sound［J］.J.Acoust.Soc.Am.，1964，36:177－195.