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

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

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

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

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

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

針對固體火箭發動機大膨脹比噴管出現的側向載荷，采用三維數值仿真方法進行分析。通過集成軟件平臺MpCCI，連接計算流體動力學軟件FLUENT和有限元軟件ABAQUS，結合FLUENT中UDF功能，對燃氣流動與噴管結構運動變形進行了耦合計算。計算結果發現，此大膨脹比噴管在快速升壓階段先后經歷了自由激波分離和受限激波分離;分析得出了較強的側向載荷主要由不對稱的燃燒波、激波轉變和噴管出口部位的激波震動3種狀態產生;側向載荷的大小也與快速升壓的時間快慢有較大關系。采用流固耦合計算方法能體現噴管的結構變形，從而更準確地反映噴管與燃氣流相互影響的真實環境，為更深入研究大膨脹比噴管側向載荷提供了支撐。

固體火箭發動機;大膨脹比噴管;快速升壓;側向載荷;流固耦合

0 引言

噴管是固體火箭發動機的基本部件之一，在許多情況下，它決定了發動機的外形和能量質量完善程度。改進噴管是提高火箭發動機性能重要途徑之一。氣動性能設計、結構強度設計都是噴管設計的重要內容。

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

嚴重的噴管側向載荷問題，會對噴管造成危害，導致噴管氣動性能下降，噴管側向載荷也是影響無控火箭射擊精度的一個主要因素。因此，對側向載荷的深入研究是非常重要的。在很大的空氣彈性變形情況下，由于流動和結構相互作用，將會引起側向載荷巨大的增長［6］。以往的仿真計算大多單獨進行流場計算，但它沒有反映噴管結構與燃氣流相互影響的真實環境，因此缺乏足夠的準確性。隨著數值算法的發展，為得到高質量的數值仿真結果，采用流固耦合方法對其進行耦合求解逐漸成為一種趨勢。該方法可揭示內流場和結構相互影響及規律，較真實地反映發動機工作狀態，提高固體火箭發動機設計水平［7］。

文中運用MpCCI(Mesh—based parallel Code Coupling Interface)耦合器，作為計算流體動力學軟件FLUENT和有限元分析軟件ABAQUS的數據交換平臺，對燃氣流動與噴管結構運動變形進行了耦合計算。

1 流固耦合數理模型

1.1 物理模型

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

1.2 計算模型

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

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

圖1 流場計算網格模型Fig.1 Layout of computational grid for flow field

1.2.1 結構參數

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

圖2 結構計算網格模型Fig.2 Layout of computational grid for the nozzle

1.2.2 初始條件

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

1.2.3 模型驗證

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

表1 燃燒室壓強為4 MPa時計算模型網格數Table 1 Grid distributions under 4 MPa combustion pressure

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

使用有限元軟件ABAQUS來模擬噴管工作時受力狀況，通過單獨使用軟件ABAQUS進行計算。結果表明，該網格可較好地反映結構變形與結構應力。

1.2.4 計算方法

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

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

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

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

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

2 結果分析

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

2.1 快速升壓時間對側向載荷產生的影響

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

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

2.2 流場分析

為了進一步分析快速升壓過程中側向載荷產生的原因，文中以10 ms快速升壓時間為例進行流場分析。

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

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

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

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

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

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

2.3 結構分析

耦合過程中，通過有限元軟件ABAQUS計算結果，顯示了噴管結構在不同時刻的應力、應變分布圖。

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

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

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

圖7 9 ms時噴管結構Mises應力圖Fig.7 Mises contour of the nozzle at 9 ms

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

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

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

3 結論

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

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

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

［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］胡海峰，鮑福廷，王藝杰.噴管分離流流動-熱-結構順序耦合數值模擬及試驗研究［J］.宇航學報，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，Fuchiwaki 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］陳洪章，賈慶英.小型固體火箭發動機壓力-時間曲線預估技術的研究［J］.彈箭與制導學報，1981(1).

［9］楊月誠.火箭發動機理論基礎［M］.西安:西北工業大學出版社，2010:95-133.

［10］武曉松，陳軍，王棟，等.固體火箭發動機工作過程數值仿真［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—)，男，教授/博導，研究方向為航空宇航推進理論與工程。E-mail:wupeng_ch@163.com

(編輯:崔賢彬)