楊月誠，吳朋朋，高雙武，錢治凱

(第二炮兵工程大學(xué) 二系201室，西安 710025)

快速升壓過程噴管側(cè)向載荷流固耦合分析①

楊月誠，吳朋朋，高雙武，錢治凱

(第二炮兵工程大學(xué) 二系201室，西安 710025)

針對固體火箭發(fā)動機(jī)大膨脹比噴管出現(xiàn)的側(cè)向載荷，采用三維數(shù)值仿真方法進(jìn)行分析。通過集成軟件平臺MpCCI，連接計算流體動力學(xué)軟件FLUENT和有限元軟件ABAQUS，結(jié)合FLUENT中UDF功能，對燃?xì)饬鲃优c噴管結(jié)構(gòu)運動變形進(jìn)行了耦合計算。計算結(jié)果發(fā)現(xiàn)，此大膨脹比噴管在快速升壓階段先后經(jīng)歷了自由激波分離和受限激波分離;分析得出了較強(qiáng)的側(cè)向載荷主要由不對稱的燃燒波、激波轉(zhuǎn)變和噴管出口部位的激波震動3種狀態(tài)產(chǎn)生;側(cè)向載荷的大小也與快速升壓的時間快慢有較大關(guān)系。采用流固耦合計算方法能體現(xiàn)噴管的結(jié)構(gòu)變形，從而更準(zhǔn)確地反映噴管與燃?xì)饬飨嗷ビ绊懙恼鎸嵀h(huán)境，為更深入研究大膨脹比噴管側(cè)向載荷提供了支撐。

固體火箭發(fā)動機(jī);大膨脹比噴管;快速升壓;側(cè)向載荷;流固耦合

0 引言

噴管是固體火箭發(fā)動機(jī)的基本部件之一，在許多情況下，它決定了發(fā)動機(jī)的外形和能量質(zhì)量完善程度。改進(jìn)噴管是提高火箭發(fā)動機(jī)性能重要途徑之一。氣動性能設(shè)計、結(jié)構(gòu)強(qiáng)度設(shè)計都是噴管設(shè)計的重要內(nèi)容。

隨著推進(jìn)技術(shù)的發(fā)展，運載火箭的助推級或第一級發(fā)動機(jī)正采用越來越大面積比的噴管，以提高高空性能。但大面積比噴管在地面試車以及發(fā)動機(jī)的啟動和關(guān)機(jī)過程中，經(jīng)常產(chǎn)生分離流動現(xiàn)象。在瞬變狀況下，分離流動經(jīng)常呈現(xiàn)非常復(fù)雜的非軸對稱性，造成嚴(yán)重的側(cè)向載荷。如美國的J2發(fā)動機(jī)和航天飛機(jī)主發(fā)動機(jī)(SSME)、俄羅斯的RD-0120發(fā)動機(jī)、歐洲的火神發(fā)動機(jī)和日本的LE-7A等發(fā)動機(jī)［1-2］在研制過程中，均遇到了嚴(yán)重的側(cè)向載荷。國內(nèi)外學(xué)者對于分離流已經(jīng)開展了一些研究，胡海峰等［3］采用流固耦合中的流動與換熱模型，實現(xiàn)了溫度場解算與流場解算的耦合數(shù)值分析。Francesco Nasuti等［4］研究了大面積比固體火箭發(fā)動機(jī)噴管的3種分離結(jié)構(gòu):自由激波分離(FSS)、受限激波分離(RSS)、含有無粘分離的自由激波分離。Ralf H Stark等［5］對噴管在啟動和關(guān)閉階段，由于分離流產(chǎn)生的側(cè)向載荷進(jìn)行了試驗研究。

嚴(yán)重的噴管側(cè)向載荷問題，會對噴管造成危害，導(dǎo)致噴管氣動性能下降，噴管側(cè)向載荷也是影響無控火箭射擊精度的一個主要因素。因此，對側(cè)向載荷的深入研究是非常重要的。在很大的空氣彈性變形情況下，由于流動和結(jié)構(gòu)相互作用，將會引起側(cè)向載荷巨大的增長［6］。以往的仿真計算大多單獨進(jìn)行流場計算，但它沒有反映噴管結(jié)構(gòu)與燃?xì)饬飨嗷ビ绊懙恼鎸嵀h(huán)境，因此缺乏足夠的準(zhǔn)確性。隨著數(shù)值算法的發(fā)展，為得到高質(zhì)量的數(shù)值仿真結(jié)果，采用流固耦合方法對其進(jìn)行耦合求解逐漸成為一種趨勢。該方法可揭示內(nèi)流場和結(jié)構(gòu)相互影響及規(guī)律，較真實地反映發(fā)動機(jī)工作狀態(tài)，提高固體火箭發(fā)動機(jī)設(shè)計水平［7］。

文中運用MpCCI(Mesh—based parallel Code Coupling Interface)耦合器，作為計算流體動力學(xué)軟件FLUENT和有限元分析軟件ABAQUS的數(shù)據(jù)交換平臺，對燃?xì)饬鲃优c噴管結(jié)構(gòu)運動變形進(jìn)行了耦合計算。

1 流固耦合數(shù)理模型

1.1 物理模型

基于主要研究噴管的側(cè)向載荷，因此把大面積比噴管和外場作為研究對象。文中以某固體火箭發(fā)動機(jī)噴管為研究背景。該噴管總長1.6 m，噴管面積比ε=55.2。

1.2 計算模型

計算區(qū)域為噴管和外場。流場區(qū)域由于噴管結(jié)構(gòu)變形，需采用動網(wǎng)格技術(shù)來模擬，所以接近噴管內(nèi)壁的流場局部區(qū)域采用了四面體非結(jié)構(gòu)化網(wǎng)格，為減小網(wǎng)格對流場非對稱性的影響，流場局部區(qū)域采用的四面體非結(jié)構(gòu)化網(wǎng)格采用預(yù)先劃分網(wǎng)格節(jié)點，以盡量減小網(wǎng)格在yz截面上的不對稱因素。流場其他區(qū)域均采用六面體結(jié)構(gòu)化網(wǎng)格(見圖1)。

有限元軟件ABAQUS用于構(gòu)造噴管結(jié)構(gòu)，噴管結(jié)構(gòu)共劃分26 326個網(wǎng)格，單元類型為C3D8R(見圖2)。

圖1 流場計算網(wǎng)格模型Fig.1 Layout of computational grid for flow field

1.2.1 結(jié)構(gòu)參數(shù)

文中所選噴管結(jié)構(gòu)的主要性能參數(shù)為密度ρ=7 850 kg/m3，泊松比γ=0.33，彈性模量E=2.0 ×1011Pa，屈服強(qiáng)度σs=1 000 MPa。

圖2 結(jié)構(gòu)計算網(wǎng)格模型Fig.2 Layout of computational grid for the nozzle

1.2.2 初始條件

文中對點火過程中的快速升壓階段進(jìn)行研究。參考典型的燃燒室壓強(qiáng)、溫度時間變化曲線［8-10］。以快速升壓階段起點為時間0點，分別對5 ms快速升壓(噴管入口壓強(qiáng)在5 ms時間段由1 MPa升至5 MPa，噴管入口溫度在4 ms時間段由500 K升至1 500 K);10 ms快速升壓(噴管入口壓強(qiáng)在10 ms時間段由1 MPa升至5 MPa，噴管入口溫度在8 ms時間段由500 K升至1 500 K);20 ms快速升壓(噴管入口壓強(qiáng)在20 ms時間段由1 MPa升至5 MPa，噴管入口溫度在16 ms時間段由500 K升至1 500 K)進(jìn)行研究。噴管結(jié)構(gòu)內(nèi)表面為耦合邊界。噴管流場取海平面的壓強(qiáng)與溫度作為流固耦合初始條件。噴管頭部壁面為絕熱固壁邊界。壓力遠(yuǎn)場邊界條件取海平面的壓強(qiáng)和溫度。

1.2.3 模型驗證

由于側(cè)向載荷較難算準(zhǔn)，這是當(dāng)前學(xué)術(shù)界公認(rèn)的問題，側(cè)向載荷對網(wǎng)格密度有較強(qiáng)的依賴性。因此，模型驗證工作比較重要。為選取合適的計算模型網(wǎng)格，文中對3種不同計算模型網(wǎng)格數(shù)，使用計算流體動力學(xué)軟件FLUENT對流場進(jìn)行網(wǎng)格無關(guān)性驗證，范圍從網(wǎng)格668 561(網(wǎng)格A)到1 570 975個網(wǎng)格(網(wǎng)格C)。文中主要研究噴管側(cè)向載荷現(xiàn)象，通過設(shè)定入口壓強(qiáng)恒為4 MPa時，對0～5 ms時間段內(nèi)，測得的Y軸方向上的側(cè)向力峰值及產(chǎn)生時間點進(jìn)行對比。計算模型網(wǎng)格無關(guān)性驗證在表1中進(jìn)行了顯示。

表1 燃燒室壓強(qiáng)為4 MPa時計算模型網(wǎng)格數(shù)Table 1 Grid distributions under 4 MPa combustion pressure

分析發(fā)現(xiàn)，網(wǎng)格B與網(wǎng)格C在此階段側(cè)向載荷峰值的差異為3.6%，時間點的差異為1.7%，網(wǎng)格數(shù)為1 158 635(網(wǎng)格B)的計算模型網(wǎng)格已可較好的用于文中研究內(nèi)容，故采用網(wǎng)格B作為發(fā)動機(jī)內(nèi)流場計算模型網(wǎng)格。

使用有限元軟件ABAQUS來模擬噴管工作時受力狀況，通過單獨使用軟件ABAQUS進(jìn)行計算。結(jié)果表明，該網(wǎng)格可較好地反映結(jié)構(gòu)變形與結(jié)構(gòu)應(yīng)力。

1.2.4 計算方法

流場計算中，控制方程采用N-S方程，湍流模型采用RNGk-ε模型，粘性按Sutherland定律給定。Sutherland定律采用三系數(shù)形式:

式中T0為參考溫度;μ0為T0時的參考粘性系數(shù);S為等效溫度。

對于動態(tài)網(wǎng)格的處理，具有運動邊界的任意有限大小的控制容積，積分形式的守恒方程為

式中ρ為流體密度;φ為待求變量為流動速度為網(wǎng)格移動速度;Γ為擴(kuò)散系數(shù);Sφ為源項;?V為有限容積的控制邊界。

在整個耦合過程中，采用12個處理器，時間步為10-5s。發(fā)動機(jī)內(nèi)流場由FLUENT軟件計算，采用耦合隱式求解方法，對流項和粘性項的離散采用二階迎風(fēng)格式，整個計算具有二階精度。噴管結(jié)構(gòu)運動變形計算采用ABAQUS隱式求解模塊——ABAQUS/Standard模塊。雙方在耦合區(qū)域部分的網(wǎng)格可不匹配，而網(wǎng)格數(shù)據(jù)之間的轉(zhuǎn)換是通過MpCCI的插值來實現(xiàn)，從而將FLUENT軟件和ABAQUS軟件每一個迭代步的計算結(jié)果進(jìn)行數(shù)據(jù)交換。

2 結(jié)果分析

在耦合過程中，通過計算流體動力學(xué)軟件FLUENT計算結(jié)果，顯示了噴管不同時刻的側(cè)向載荷與流場的變化。

2.1 快速升壓時間對側(cè)向載荷產(chǎn)生的影響

圖3為噴管yz方向上側(cè)向載荷數(shù)值Fyz隨時間變化曲線。由圖3可見，隨著氣流向下游傳播，側(cè)向載荷數(shù)值先后經(jīng)歷了由不對稱的燃燒波、激波轉(zhuǎn)變和噴管出口部位的激波震動3種狀態(tài)導(dǎo)致的較強(qiáng)側(cè)向載荷。5 ms快速升壓時間條件下產(chǎn)生的側(cè)向載荷峰值為1 520 N，10 ms快速升壓時間條件下產(chǎn)生的側(cè)向載荷峰值為2 453 N，20 ms快速升壓時間條件下產(chǎn)生的側(cè)向載荷峰值為6 000 N。對各種快速升壓時間條件下產(chǎn)生的側(cè)向載荷峰值進(jìn)行對比可見，5 ms快速升壓過程產(chǎn)生的側(cè)向載荷普遍比10 ms快速升壓過程產(chǎn)生的側(cè)向載荷要小，側(cè)向載荷峰值是10 ms快速升壓過程側(cè)向載荷峰值的3/5。20 ms快速升壓過程產(chǎn)生的側(cè)向載荷普遍比10 ms快速升壓過程產(chǎn)生的側(cè)向載荷要大，側(cè)向載荷峰值是10 ms快速升壓過程側(cè)向載荷峰值的2.5倍。由此可得出，對于此噴管快速升壓時間，將對噴管側(cè)向載荷造成較大的影響。一定時間范圍內(nèi)，較短的快速升壓時間將會明顯減小側(cè)向載荷峰值。

圖3 不同升壓時間條件下噴管側(cè)向載荷時間變化曲線Fig.3 Computed side forces for the nozzle

2.2 流場分析

為了進(jìn)一步分析快速升壓過程中側(cè)向載荷產(chǎn)生的原因，文中以10 ms快速升壓時間為例進(jìn)行流場分析。

在0～0.5 ms階段，噴管側(cè)向載荷數(shù)值很小，以圖3中10 ms快速升壓為例，在0.5 ms時，噴管側(cè)向載荷僅為26.729 2 N。0.5 ms之后，噴管側(cè)向載荷迅速躍升，在0.812 6 ms時達(dá)到極大值，此時噴管側(cè)向載荷為1 126 N。圖4為0.812 6 ms時噴管XY平面截圖，顯示了此時的噴管內(nèi)流場速度和溫度云圖。由圖4可看出，噴管內(nèi)流場出現(xiàn)較明顯不對稱現(xiàn)象。在0.812 6 ms時，噴管入口氣流正在經(jīng)歷壓強(qiáng)與溫度的迅速躍升，呈現(xiàn)快速膨脹的壓力波與快速升高的溫度。因此，在0.812 6 ms時出現(xiàn)的較大側(cè)向載荷主要看作是燃燒波導(dǎo)致。

0.812 6 ms之后，側(cè)向載荷處于波動狀態(tài)，燃燒波逐漸向約束激波轉(zhuǎn)變，在激波轉(zhuǎn)變的過程中又出現(xiàn)了較大的側(cè)向載荷的波動。之后，在1.1～8.34 ms階段主要是約束激波流動模態(tài)。此時間段內(nèi)，沒有出現(xiàn)較大的側(cè)向載荷波動。

8.34 ms之后，噴管側(cè)向載荷再一次迅速躍升，在9.05 ms時達(dá)到快速升壓階段的最大值，此時噴管側(cè)向載荷為2 453 N。

圖4 0.812 6 ms時噴管內(nèi)流場速度及溫度云圖Fig.4 Velocity and temperture contour of the nozzle at 0.812 6 ms

圖5(a)、(b)為9.05 ms時噴管XY平面截圖，顯示了此時的噴管內(nèi)流場速度和溫度云圖，可見內(nèi)流場出現(xiàn)明顯不對稱現(xiàn)象。圖5(c)為9.05 ms時，噴管內(nèi)流流場馬赫數(shù)云圖，清晰地顯示了此時的流場狀態(tài)。導(dǎo)致此階段出現(xiàn)側(cè)向載荷最大值的原因是噴管出口部位的激波震動，此時約束激波在噴管出口部位來回移動，而由于激波移動的不對稱性，出現(xiàn)了與噴管口相接觸的激波部分已經(jīng)在噴管口外面，另一部分激波卻還在噴管出口部位，這樣就必然導(dǎo)致了更為嚴(yán)重的側(cè)向載荷。此時，噴管流場還呈現(xiàn)馬赫盤激波模態(tài)［11］(Mach disk shock pattern)。在此階段也會發(fā)現(xiàn)馬赫盤來回移動的現(xiàn)象，這些現(xiàn)象都主要因為入口壓強(qiáng)不斷增加，流場區(qū)域不斷進(jìn)行調(diào)整，以適應(yīng)噴管擴(kuò)張段流動而導(dǎo)致。受噴管內(nèi)馬赫盤的強(qiáng)烈阻擋，氣流有繞開其流動的趨勢，從而使馬赫盤下游形成了低速區(qū)。低速區(qū)內(nèi)的氣流溫度較高，這在噴管內(nèi)流場溫度云圖上反映很明顯。此階段最終在12.49 ms之后趨于穩(wěn)定。原先在噴管出口部位的激波完全進(jìn)入外場，側(cè)向載荷回落至平穩(wěn)的狀態(tài)，此階段側(cè)向載荷數(shù)值很小，且噴管出口質(zhì)量流量平穩(wěn)。將此階段稱為噴管內(nèi)流場充分流動狀態(tài)，如圖6所示。

圖5 9.05 ms時噴管內(nèi)流場速度、溫度及馬赫數(shù)云圖Fig.5 Velocity，temperture and Mach contour of the nozzle at 9.05 ms

2.3 結(jié)構(gòu)分析

耦合過程中，通過有限元軟件ABAQUS計算結(jié)果，顯示了噴管結(jié)構(gòu)在不同時刻的應(yīng)力、應(yīng)變分布圖。

圖6 15 ms時噴管內(nèi)流場速度云圖Fig.6 Velocity contour of the nozzle at 15 ms

圖7為9 ms時噴管結(jié)構(gòu)Mises應(yīng)力云圖。有限元軟件ABAQUS計算得到的Mises等效應(yīng)力遵循材料力學(xué)的第四強(qiáng)度理論，用于判斷結(jié)構(gòu)是否發(fā)生塑性變形。由圖7可見，此時噴管結(jié)構(gòu)最大的Mises應(yīng)力約為20 MPa，遠(yuǎn)小于噴管結(jié)構(gòu)屈服強(qiáng)度1 000 MPa。因此，此時噴管不會發(fā)生破壞。

圖8 為噴管外壁中點(1.6，0.52，0)的應(yīng)力時間變化曲線，圖9為噴管出口外壁中點(1.6，0.52，0)的位移時間變化曲線。

圖7 9 ms時噴管結(jié)構(gòu)Mises應(yīng)力圖Fig.7 Mises contour of the nozzle at 9 ms

圖8 噴管外壁中點的應(yīng)力時間變化曲線Fig.8 Stress of the center of the outer wall

由圖8、圖9可見，曲線較好反映了流場參數(shù)變化對結(jié)構(gòu)造成的影響。對于此噴管厚度較大，結(jié)構(gòu)變形相應(yīng)較小。噴管結(jié)構(gòu)的微小變形對噴管壁壓也相當(dāng)敏感［12］。對于大型噴管如美國的J2發(fā)動機(jī)噴管，以及軍用固體火箭發(fā)動機(jī)噴管，需要追求輕質(zhì)殼體，均會導(dǎo)致結(jié)構(gòu)變形較大，結(jié)構(gòu)發(fā)生變化又必定對流場造成影響，這也是造成分離流不穩(wěn)定性的一個重要因素，可作為解釋在燃燒室壓力不變時分離位置不斷移動的一個原因。

圖9 噴管外壁中點的位移時間變化曲線Fig.9 Displacement of the center of the outer wall

3 結(jié)論

(1)文中對大膨脹比噴管流動采用了三維流固耦合數(shù)值模擬，展現(xiàn)了噴管流場與結(jié)構(gòu)之間的相互影響，并對噴管結(jié)構(gòu)的應(yīng)力、位移進(jìn)行分析。通過對FLUENT軟件二次開發(fā)及數(shù)據(jù)后處理，展現(xiàn)了側(cè)向載荷隨時間變化曲線。經(jīng)分析表明，較強(qiáng)的側(cè)向載荷主要由不對稱的燃燒波、激波轉(zhuǎn)變和噴管出口部位的激波震動3種狀態(tài)產(chǎn)生，而且側(cè)向載荷的大小也與快速升壓的時間快慢有較大關(guān)系。在一定時間范圍內(nèi)，隨著快速升壓時間的增長，側(cè)向載荷峰值將會急劇增大。

(2)MpCCI軟件能很好地將計算流體動力學(xué)軟件FLUENT和有限元軟件ABAQUS聯(lián)合進(jìn)行流固耦合計算，同步對噴管流場和噴管結(jié)構(gòu)進(jìn)行分析。由于其較好的功能，這種方法可被用來構(gòu)造多種工況下發(fā)動機(jī)的失效機(jī)理，尤其是準(zhǔn)確呈現(xiàn)發(fā)動機(jī)的工作過程，為合理設(shè)計發(fā)動機(jī)和延壽工作提供準(zhǔn)確的技術(shù)支持。

(3)分離流狀態(tài)下的流固耦合分析非常復(fù)雜，文中方法只重點考慮了在一定條件下噴管流場與結(jié)構(gòu)的相互影響，后續(xù)還應(yīng)通過大量試驗并進(jìn)行對比，不斷修改完善分析模型，為工程實際應(yīng)用打下基礎(chǔ)。

［1］Hadjadj A，Onofri M.Nozzle flow separation［J］.Shock Waves，2009，19(3):163-169.

［2］Yasuhide W，Norio S.LE-7A engine nozzle flow separation phenomenon and the possibility of RSS suppression by the step inside the nozzle［R］.AIAA 2004-4014.

［3］胡海峰，鮑福廷，王藝杰.噴管分離流流動-熱-結(jié)構(gòu)順序耦合數(shù)值模擬及試驗研究［J］.宇航學(xué)報，2011，32(7):1534-1541.

［4］Francesco Nasuti，Marcello Onofri.Numerical Analysis of Flow Separation Structures in Rocket Nozzles［R］.AIAA 2007-5473.

［5］Ralf H Stark，Chloe Genin.Experimental study on rocket Nozzle side load reduction［R］.AIAA 2011-389.

［6］Zhang S J，F(xiàn)uchiwaki T.Aeroelastic coupling and side loads in rocket Nozzles［R］.AIAA 2008-4064.

［7］Dennis C W，Belfield J.Towards identifying rocket motor failure modes using coumputational fluid dynamic and finite element codes［R］.AIAA 2008-5209.

［8］陳洪章，賈慶英.小型固體火箭發(fā)動機(jī)壓力-時間曲線預(yù)估技術(shù)的研究［J］.彈箭與制導(dǎo)學(xué)報，1981(1).

［9］楊月誠.火箭發(fā)動機(jī)理論基礎(chǔ)［M］.西安:西北工業(yè)大學(xué)出版社，2010:95-133.

［10］武曉松，陳軍，王棟，等.固體火箭發(fā)動機(jī)工作過程數(shù)值仿真［M］.北京:高等教育出版社，2006:257-266.

［11］Frey M，Hagemann G.Status of flow separation prediction in rocket nozzles［C］//34th AIAA/ASME/SAE/ASEE Joint Propulsion Conference.

［12］Nave L H，Coffey G A.Sea level side loads in high-area-ratio rocket engines［R］.AIAA 73-1284.

Rapid pressurization side load fluid-structure coupled analysis in SRM nozzle

YANG Yue-cheng，WU Peng-peng，GAO Shuang-wu，QIAN Zhi-kai
(Second Artillery Engineering University No201 Staff Room，Xi＇an 710025，China)

A side load in the overexpanded nozzle was studied by means of a three-dimensional numerical simulation method.The MpCCI software was used to link the FLUENT CFD code whose UDF is applied with the ABAQUS FE code to analyze gas flow and nozzle deformation.Results show that three types of asymmetric shock physics incur strong side loads:the generation of combustion wave，shock transitions，and shock pulsations across the nozzle lip.Moreover，the ramp time results show that the ramp rate has the potential effect on the nozzle side loads.The fluid-structure coupled numerical simulation support the accuracy of the method of numerical simulation.The simulation provides the base for further study.

solid rocket motor;over-expanded nozzle;rapid pressurization;side load;fluid-structure coupled analysis

V435

A

1006-2793(2012)04-0463-05

2012-04-13;

2012-06-08。

楊月誠(1954—)，男，教授/博導(dǎo)，研究方向為航空宇航推進(jìn)理論與工程。E-mail:wupeng_ch@163.com

(編輯:崔賢彬)