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

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

(西北工業(yè)大學(xué)航天學(xué)院，西安 710072)

0 引言

目前的運載助推級動力系統(tǒng)，為了發(fā)揮發(fā)動機的性能潛力，采用的噴管膨脹比越來越大。大膨脹比的噴管在地面試車及發(fā)動機的起動和關(guān)機過程中，都會產(chǎn)生氣流分離現(xiàn)象。氣流分離模態(tài)的轉(zhuǎn)涙，導(dǎo)致在瞬變狀態(tài)下噴管內(nèi)三維氣流的分離流動常呈現(xiàn)復(fù)雜的非軸對稱性，可能引發(fā)噴管側(cè)向載荷。嚴(yán)重的噴管側(cè)向載荷會造成噴管震動、發(fā)動機使用壽命縮短、噴管結(jié)構(gòu)破壞等后果。

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

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

1 研究對象與方法

1.1 研究對象

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

幾何構(gòu)型見圖1，幾何參數(shù)見表1。

圖1 噴管基本構(gòu)型簡圖Fig.1 Basic nozzle geometry sketch

表1 計算幾何結(jié)構(gòu)基本參數(shù)Table 1 calculation geometry data

1.2 流場數(shù)值計算方法

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

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

1.3 網(wǎng)格無關(guān)性驗證

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

表2 NPR=14、16時計算模型網(wǎng)格數(shù)Table 2 Grid disributions when NPR=14，16

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

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

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

1.4 流場結(jié)構(gòu)分析

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

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

1.4.1 自由激波(FSS)

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

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

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

1.4.2 受限激波(RSS)

圖4所示受限激波結(jié)構(gòu)簡圖中，分離后的激波附著到了噴管壁面，形成一個附著回流區(qū)。

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

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

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

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

1.5 流固耦合分析計算

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

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

在流場計算方面，通過前述方法計算N-S方程;在結(jié)構(gòu)場計算方面，求解結(jié)構(gòu)靜力學(xué)平衡方程。

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

圖7 CFD/CSD松耦合結(jié)構(gòu)示意圖Fig.7 CFD/CSD loose coupled structure

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

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

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

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

2 仿真結(jié)果分析

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

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

表3 噴管殼體材料數(shù)據(jù)Table 3 Property of the nozzle case material

圖10 不同時刻噴管內(nèi)部流場和噴管殼體結(jié)構(gòu)變形云圖Fig.10 Contour of the flow inside nozzle and deformation of nozzle case at different times

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

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

3 結(jié)論

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

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

值得指出的是噴管內(nèi)部的流動較為復(fù)雜，分離狀況下產(chǎn)生的側(cè)向載荷通常表現(xiàn)為非對稱性。因此，為詳細(xì)揭示分離條件下結(jié)構(gòu)響應(yīng)，需開展三維分析。

［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］劉亞冰，王長輝，許曉勇.噴管分離流動與側(cè)向載荷定常數(shù)值模擬［J］.航空動力學(xué)報，2008，23(11):2115-2118.

［4］Emanuele Martelli，F(xiàn)rancesco 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］王藝杰，鮑褔廷，杜佳佳.固體火箭發(fā)動機噴管分離流流動數(shù)值模擬及實驗研究［J］.固體火箭技術(shù)，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 ，F(xiàn)uchiwaki 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.