大膨脹比火箭發動機噴管流動分離與氣動彈性分析①

胡海峰，鮑福廷，蔡 強，劉 旸

(西北工業大學航天學院，西安 710072)

0 引言

目前的運載助推級動力系統，為了發揮發動機的性能潛力，采用的噴管膨脹比越來越大。大膨脹比的噴管在地面試車及發動機的起動和關機過程中，都會產生氣流分離現象。氣流分離模態的轉涙，導致在瞬變狀態下噴管內三維氣流的分離流動常呈現復雜的非軸對稱性，可能引發噴管側向載荷。嚴重的噴管側向載荷會造成噴管震動、發動機使用壽命縮短、噴管結構破壞等后果。

美國的J2S發動機、航天飛機主發動機(SSME)、俄羅斯的RD-120發動機、歐洲的Vulcain發動機和日本的LE-7等發動機的研制過程中，均遇到了較嚴重的噴管側向載荷問題。國外對大膨脹比的噴管性能分析計算開展了大量研究，并取得了很多有用結論。文獻［1］對噴管分離流試驗結果進行分析歸納，給出了判別分離流的經驗公式;文獻［2-3］對分離流產生的側向載荷進行了分析仿真;文獻［4-6］對影響分離流各個分離模式進行了數值仿真及實驗對比分析研究;文獻［7］對分離流下的噴管熱性能進行了探討;文獻［8］開展了分離狀態下噴管氣動響應的理論探索;文獻［9］對某型噴管開展了分離模式下結構氣動響應的試驗研究;文獻［10-12］開展了分離條件下噴管流固耦合數值模擬。

本文以某大膨脹比噴管二維軸對稱簡化模型為分析對象，采用數值模擬方法研究其分離模式，分析在瞬變狀態下噴管結構的氣動彈性響應，討論了噴管瞬態下徑向位移分布。側向載荷在三維狀態下能夠預示，為二維簡化三維的分析提供基礎，并為后續三維深入分析側向載荷及噴管結構氣動彈性響應提供支持。

1 研究對象與方法

1.1 研究對象

本研究選取VAC公司的Volvo S1噴管為研究對象［6］。該噴管為最大推力噴管，其擴張段內型面滿足Rao-Shmyglevsky方程:

幾何構型見圖1，幾何參數見表1。

圖1 噴管基本構型簡圖Fig.1 Basic nozzle geometry sketch

表1 計算幾何結構基本參數Table 1 calculation geometry data

1.2 流場數值計算方法

本文采用SA湍流模型，求解直角坐標系下二維穩態雷諾時均N-S方程。這是由于該模型相對于一般的兩方程湍流模型，對逆壓梯度下的分離流動有著更為準確的模擬。該湍流模型對壁面網格節點的要求更為嚴格。因此，采用結構化網格，壁面網格局部加密，保證第一層網格的y+接近于1，并保證在有粘性影響的近壁面區域(Rev＜200)內，至少有10個單元網格，網格均勻過渡。對N-S湍流時均方程，采用耦合隱式求解;對各參數的離散，采用二階精度的迎風格式。為準確描述理想氣體各物性，數值計算中根據經驗多項式進行定壓比熱容計算，氣體粘性系數由三系數的Sutherland定律給定。Sutherland定律具體形式為

式中 T0為參考溫度，取為273.11 K;S為等效溫度，取為110.56 K;μ0為 T0時的參考粘性系數，取為1.716×10-5kg/(m·s)。

1.3 網格無關性驗證

分離流對網格密度有很強的依賴性，本文在計算獲得穩態解的基礎上，對網格進行加密對比計算結果。當網格加密后，計算結果不發生變化，即認為計算得到的收斂解與網格無關。本文在保證壁面y+接近1的前提下，計算得到收斂解，然后對網格加密，適當調整網格間距。為區別比較，本文在流向及垂直流向方向(徑向)加密。本文列舉燃燒室總壓與環境背壓比值(NPR)為14、16兩種工況進行對比分析。表2所示為具體的噴管區域給出網格的分布及總網格數，外流場網格數據未列舉。

表2 NPR=14、16時計算模型網格數Table 2 Grid disributions when NPR=14，16

圖2、圖3分別為NPR=14、16時，計算2組不同疏密條件下得到的穩態收斂解中噴管壁面壓力分布與試驗數據的對比曲線。由圖2、圖3可直觀看出，計算結果和試驗結果吻合較好。在網格加密后，壁面壓強分布和之前相對較稀疏網格比較基本沒有變化。說明在計算過程中，得到的穩態收斂解和網格無關，基本揭示了真實的噴管內部流動。

圖2 NPR=14下不同網格類型下噴管壁面壓力分布線Fig.2 Nozzle wall pressure when NPR=14

圖3 NPR=16下不同網格類型下噴管壁面壓力分布線Fig.3 Nozzle wall pressure when NPR=16

1.4 流場結構分析

該數值模擬部分對 NPR 分別為 8、10、12、14、16.4、21、25、30、35、40、45、50 等 12 種工況，按照上述的數值方法進行仿真分析。每個工況均以上一計算結果為初值。結果表明，入口壓強由高到低不同工況下，噴管流場先后出現2種不同的激波分離模態:自由激波和受限激波分離，如圖4所示。2種不同激波分離模態下，噴管內流場及壓強分布情況等呈現不同特點。圖4中，a為馬赫盤;b為斜激波;c為受限激波;d為內激波。

圖4 自由激波和受限激波簡圖Fig.4 Schematic illustrations of FSS and RSS

1.4.1 自由激波(FSS)

圖4所示自由激波結構簡圖中，噴管內的氣流由壁面分離后，向外延伸流出噴管，由于噴管膨脹比較大，噴管內氣流處于過膨脹狀態，燃燒室總壓與環境背壓之比較小，導致在噴管內部出現激波。但由圖4可看到，氣流經過正激波a及斜激波b后，在內部形成馬赫盤，壁面的氣體分離后，未附著在壁面上。同時，隨燃燒室總壓的提高，噴管內部的激波往噴管出口移動，但在激波前，流場的參數基本不變。

圖5為燃燒室在自由激波模態下噴管擴張段的壁面壓強分布曲線。由圖5可看出，在自由激波模式下，壁面出現分離后，壓強經過激波迅速上升到環境壓強大小。

圖5 自由激波階段噴管壁面壓強分布Fig.5 Profile of the pressure distributionsalong the nozzle wall

1.4.2 受限激波(RSS)

圖4所示受限激波結構簡圖中，分離后的激波附著到了噴管壁面，形成一個附著回流區。

圖6所示為燃燒室在受限激波模態下，噴管擴張段的壁面壓強分布。

圖6 受限激波階段噴管壁面壓強分布Fig.6 Distributions of the pressure along nozzle wall during the RSS model

由圖6可看出，壁面壓強同樣首先隨氣流膨脹逐漸降低，在分離點由于分離斜激波的作用，壁面壓強迅速升高到一個相對穩定壓強(圖6中所示平臺區域)。與自由激波分離不同的是短暫分離后，氣流再次附著到壁面。由于再附著激波的強烈作用，在再附著點壁面壓強急劇升高到超過環境壓強。此后，壁面壓強通過一系列逐漸減弱的波動，最終回落到略低于環境壓強，直到噴管出口。這種逐漸減弱的壓強波動是由于激波的反射特性造成的。再附著激波有反射離開壁面的趨勢，這種趨勢的強弱由激波的強弱決定，由此造成氣流在遠離和靠近壁面之間反復，在受限區域內形成了菱形激波串，進而引起壁面壓強的波動。

上述兩種不同的分離模態下，每種分離激波模態噴管內部都會產生一道內激波(圖4中d曲線)。這道內激波源于噴管拋物線型面起始點。最大推力噴管的設計方法，決定了噴管型面在這一點為二階不連續點，是一個流動不穩定點，容易產生激波。圖4流場結構示意圖表明，這道內激波將在很大程度上決定中心區激波的形態。中心區流場正激波后的渦流具有推動流動向壁面發展的趨勢。當流體朝向壁面流動的動量達到一定程度時，便產生了分離再附著現象。

1.5 流固耦合分析計算

本文基于計算流體力學(CFD)和計算結構動力學(CSD)的時域耦合仿真方法，分析噴管流固耦合氣動彈性問題。通過非定常計算流體力學求解器，直接計算噴管殼體在任意時刻的非定常載荷，包括熱和氣動力，在時域內推進結構運動方程，給出噴管殼體結構的詳細時間響應歷程。

根據流場和結構兩個獨立域求解的時間同步推進技術，文獻［13］將CFD/CSD耦合分為全耦合、緊耦合、松耦合3種方式。本文采用松耦合，其流體計算與結構計算在時間上耦合基本邏輯見圖7。

在流場計算方面，通過前述方法計算N-S方程;在結構場計算方面，求解結構靜力學平衡方程。

流固耦合計算關鍵在于流固交界面的處理［14］，流固交界面上數據需雙向交互。流場對固體區域的作用，包括力的作用等載荷。計算時，在一個時間步內，首先計算流場，然后將流場的作用通過流固交界面加載到固體區域，計算出固體的溫度場和應力應變，如果變形較大，則改變流固界面。基于無限插值(transfinite interpolation)的變形網格技術方法，重新生成計算區域的網格，進行下一個時間步的計算。具體如圖8所示。

圖7 CFD/CSD松耦合結構示意圖Fig.7 CFD/CSD loose coupled structure

圖8 流固耦合表面數據交互Fig.8 Data exchanged of the FSI interface

為了提高計算效率，采用并行化的MDICE［15］并行計算環境，遵照MPI(Message Passing Interface)規范，實現計算過程中各不同模塊數據在不同進程間的數據傳遞。CFD/CSD松耦合的具體流程見圖9。

圖9 CFD/CSD松耦合程序流程示意圖Fig.9 Flow chart of the CFD/CSD loosely coupled

計算噴管殼體應力分布時，首先求解在外力作用 下應變的分布，再通過應力和應變的關系，求解應力的分布。外力和應變之間的關系，通過虛功原理建立平衡方程。假設是作用在噴管殼體計算單元上的外力，通過虛功原理建立壁面的力學平衡方程。節點虛位移為則外力所做的功為

2 仿真結果分析

計算中，取噴管殼體壁厚為等厚度，具體材料屬性及厚度如表3所示。模擬噴管入口壓力為3 MPa的瞬態起動工況，分析該過程中噴管殼體的氣動響應問題。

圖10所示為發動機工作起動階段不同時刻流場速度云圖和噴管殼體徑向位移變形云圖顯示。由圖10可看出，在噴管起動過程中，殼體在內壓作用下沿徑向出現變形微小變形，在1×10-6m量級水平。分析可能的原因為燃燒室內3 MPa壓強偏低，噴管壁厚11 mm較厚，同時2 ms作用時間極短，導致徑向形變較小。同時，不同時刻殼體徑向位移最大的也不是在殼體的末端。隨噴管內部壓強的建立，噴管不同部分的變形不盡相同。由云圖顯示部分可看出，在噴管收斂段部分，由于其壓強較高，導致該段變形較大。

表3 噴管殼體材料數據Table 3 Property of the nozzle case material

圖10 不同時刻噴管內部流場和噴管殼體結構變形云圖Fig.10 Contour of the flow inside nozzle and deformation of nozzle case at different times

圖11所示為噴管殼體末端的徑向位移隨時間的變化曲線?？煽闯?，在噴管起動階段，末端變形位移最大，隨噴管內部流場建立，在端點的位移逐漸退化為震蕩變形。但由圖示幅值可看到，在模擬時間內，位移均為負值，模擬仿真的時間較短，仿真時間段內未形成周期變化。圖示端點位移振蕩并非是周期性，這是因為在燃燒室室壓3 MPa條件下，噴管內分離狀態為受限激波，分離條件下的氣流再附著到壁面，再分離再附著循環，直至氣流流出噴管，在附著區域內壓強出現變化，導致其壁面受力復雜，最終表現為結構振蕩。

圖11 噴管端點徑向變形隨時間變化曲線Fig.11 Point at the end of the case radial direction deformation vs the time

3 結論

(1)對于本文計算的給定噴管構型，在燃燒室與噴管出口壓強之比(NPR)變化過程中，噴管內出現了自由激波和受限激波2種分離模態。不同的分離模態，導致噴管壁面承受不同的壓力分布。

(2)在對流動的分析基礎上，結合CFD/CSD耦合分析方法，分析了分離情況下噴管殼體的氣動彈性問題。結果表明，在起動過程中，隨噴管內部壓強建立，噴管不同部分出現不同的結構變形。在噴管擴張段，其承受的壓力較高、變形較大，噴管擴張段端點部分變形隨時間并非嚴格意義上的周期變化。這主要是因為在受限分離下，壁面的壓力出現波動。本文為后續詳細分析側向載荷提供基礎。

值得指出的是噴管內部的流動較為復雜，分離狀況下產生的側向載荷通常表現為非對稱性。因此，為詳細揭示分離條件下結構響應，需開展三維分析。

［1］Frey M，Hagemann G.Status of flow separation pre-diction in rocket nozzles［C］//34th AIAA/ASME/SAE/ASEE Joint Propulsion Conference ＆ Exhibit，Cleve-land，Oh，USA，July 13-15，1998.

［2］Kurt B Smalley，Andrew M Brown and Joseph Ruf.Flow separation side-loads excitation of rocket nozzle FEM［C］//48th AIAA/ASME/ASCE/AHS/ASC Structures，Structural Dynamics，and Materials Con.，Honolulu，Hawaii，USA.A-pril 23-26，2007.

［3］劉亞冰，王長輝，許曉勇.噴管分離流動與側向載荷定常數值模擬［J］.航空動力學報，2008，23(11):2115-2118.

［4］Emanuele Martelli，Francesco Nasuti，Marcello Onofri.Numerical calculation of FSS/RSS transition in highly overexpanded rocket nozzle flows［J］.Shock Waves，2010，20:139-146.

［5］Vincent Lijo，Heuy Dong Kim，Toshiaki Setoguchi，et al.Numerical simulation of transient flows in a rocket propulsion nozzle［J］.International Journal of Heat and Fluid Flow，2010，31:409-417.

［6］Jan ?stlund.Flow processes in rocket engine nozzles with focus on flow separation and side-loads［R］.TRITA-MEK Technical Report，2002:9.

［7］王藝杰，鮑褔廷，杜佳佳.固體火箭發動機噴管分離流流動數值模擬及實驗研究［J］.固體火箭技術，2010，33(4):406-408.

［8］Pekkari L O.Aeroelastic analysis of side load in supersonic nozzles with separated flow［C］//30th AIAA/ASME/SAE/ASEE Joint Propulsion Conference June 27-29，1994/Indianapolis，IN.AIAA 94-3377.

［9］Dr Andrew M Brown，Joseph Ruf，Darren Reed，et al.Characterization of side load phenomena using measurement of fluid/structure interaction［R］.AIAA 2002-3999.

［10］Lefrancois E.Numerical validation of a stability model for a flexible over-expanded rocket nozzle［J］.Int.J.Numer.Meth.Fluids.2005，49:349-369.

［11］Zhang S J ，Fuchiwaki T.Aeroelastic coupling and side Loads in rocket nozzles［C］//38th Fluid Dynamics Conference and Exhibit 23-26 June 2008，Seattle，Washington，AIAA 2008-4064.

［12］Luciano Garelli.Fluid-structure interaction study of the start-up of a rocket engine nozzle［J］.Computers ＆ Fluids，2010，39:1208-1218.

［13］Hurka J，Ballman J.Elastic panels in transonic flow［R］.AIAA 2001-2722.

［14］Lv X，Zhao Y，Huang X Y，et al.A matrix-free implicit unstructured multigrid finite volume method for simulating structural dynamics and fluid-structure interaction［J］.Journal of Computational Physics，2007，225:120-144.

［15］CFD-F ASTRAN V2009.0 user manual［M］.ESI CFD inc.