固體火箭發動機喉襯流場及熱結構耦合分析①

張曉光，王長輝，劉 宇，任軍學

(北京航空航天大學宇航學院，北京 100191)

0 引言

喉襯構成燃氣流動的最小通道，控制燃燒室壓強，決定發動機性能，是固體火箭發動機的關鍵部件;而其所處燃氣環境又最為惡劣。因此，喉襯也是發動機設計工作的薄弱環節。

流場及熱結構耦合分析是評估喉襯氣動特性、熱防護有效性及結構完整性的重要手段［1］。傳統的分析過程通常采用流場及熱結構順序計算的方法:先由一維等熵氣體動力學公式獲取燃氣溫度、壓強等流場參數，由半經驗公式(如Bartz公式)估算設定內壁溫度下的內壁對流換熱系數;再以燃氣溫度及內壁對流換熱系數為熱邊界條件，進行固相材料導熱計算，得到溫度分布;最后，以流場及傳熱分析所得喉襯內壁壓強與溫度場為邊界條件和載荷條件進行結構分析［2－5］。計算中一般未考慮內壁對流換熱系數隨內壁溫度的變化，未實現流動計算與傳熱計算的耦合，這與實際情況不符;另外，流場簡化計算模型及內壁對流換熱系數半經驗估算公式的應用，也對分析精度有較大影響。

針對上述問題，本文基于FLUENT流體計算軟件與ANSYS結構分析軟件，開展了喉襯流場及熱結構耦合分析。計算分為2部分:首先，由FLUENT完成流固耦合燒蝕傳熱計算;然后，由ANSYS基于FLUENT的流場及熱分析結果進行結構計算。計算實現了流場、熱的雙向耦合及流場、熱到結構的單向耦合，兼顧了分析精度與計算資源兩方面的因素，并充分發揮了2個軟件各自優勢，有助于準確預示喉襯在燃氣環境中的燒蝕傳熱及熱結構行為。

1 計算模型

本文的計算模型為Swaminathan［6］的傳熱試驗發動機石墨喉襯，如圖1所示。在喉部位置，4個熱電偶沿圓周方向間隔90°埋置于不同深度:T1(81，20)、T2(81，24)、T3(81，29)、T4(81，34)。發動機采用聚氨酯推進劑，燃燒室壓強3.43 MPa，燃燒溫度3 203 K，燃燒時間7 s。喉襯石墨材料密度1 810 kg/m3，熱導率37.7 W/(m·K)，比定壓熱容 711.3 J/(kg·K)。

圖1 喉襯結構及熱電偶布置Fig.1 Schematic of throat insert and thermocouple locations

2 流場與燒蝕傳熱耦合計算

由FLUENT流體計算軟件采用整場離散、整場求解的方法，模擬燃氣與固相材料間的非穩態流固耦合燒蝕傳熱。計算區域包括燃氣流動區與固相材料區，各區域采用通用的控制方程，區別僅在于廣義擴散系數及廣義源項的不同;流固耦合界面成為計算區域的內部邊界，其對流換熱系數由軟件本身的計算得出，而不再是分區求解時，需要人為設定的外邊界條件;控制方程由有限體積法離散，可保證流固耦合界面滿足溫度連續、熱流密度連續條件［7］。湍流模型采用Realizable k-ε兩方程模型，近壁處理采用增強型壁面函數。

喉襯的燒蝕計算則依據熱化學燒蝕理論，由FLUENT流體計算軟件中的Wall Surface Reaction模型完成。熱化學燒蝕理論認為，喉襯燒蝕由內壁面C與燃氣氧化性組分H2O、CO2、OH等發生異相化學反應所致，反應速率遵循Arrhenius定律，而由反應所致的壁面退移速率就是燒蝕率［8］。

由NASA-CEA［9］程序進行熱力計算，結果表明，燃氣組分濃度沿喉襯軸向變化較小。因此，可認為喉襯流動為凍結流。氧化性組分 H2O、CO2質量分數為0.176、0.09，OH因濃度較小而不予考慮。燒蝕反應方程式及其動力參數見表1［8］。

表1 反應模型及動力參數Table 1 Reaction schemes and kinetics parameters

建立喉襯的二維軸對稱計算模型，并劃分四邊形結構化網格，見圖2。在喉部附近沿軸向加密，在內壁附近沿徑向加密，且使流場近壁面第一層網格y+＜1，以滿足增強型壁面函數法的要求，并更好地識別湍流邊界層。由于喉襯入口就是噴管入口［6］，且為亞聲速，因此給定總壓3.43 MPa、總溫3 203 K及進口氣流角0°，靜壓則由內場外推;喉襯出口為超聲速，全部參數外推;喉襯入口、出口的湍流條件則通過給定湍流強度和水力直徑來確定。

圖2 喉襯網格劃分Fig.2 Computational grid of throat insert

本文喉襯非穩態流固耦合燒蝕傳熱計算只針對發動機穩態工作段，而不涉及達到穩態工作之前的壓強上升段。計算時，首先設置內壁為絕熱，外壁及端面為環境溫度，進行穩態計算獲得喉襯流場分布;然后，設置外壁及端面為絕熱，內壁為流固耦合壁面及燒蝕反應發生表面，以穩態流場計算結果為初場，對流場區域與固體區域作整場求解，進行流場與燒蝕傳熱的非穩態耦合計算，計算時間7 s，時間步長取為0.000 5 s。

各熱電偶測點溫度計算值與試驗值的對比見圖3。二者基本相符，但計算值普遍高于試驗值，且距離內壁越近，該現象越明顯。分析其原因:在計算方面，外壁及端面設置為絕熱，而實際是導熱的;在試驗方面，主要是熱電偶動態熱響應特性的影響。計算與試驗均表明，距離內壁越近，溫度及溫升速率也越高。

圖3 溫度計算值與試驗值對比Fig.3 Comparison between calculated and measured temperatures

圖4為t=1、4、7 s時刻喉襯流場與固相材料的溫度分布。由圖4可見，隨發動機工作進行，熱響應逐漸深入到固相材料內部，而燃氣與內壁之間強烈的對流換熱及固相材料內部的導熱，使固相材料形成明顯的徑向溫度梯度，且以喉部區域最為顯著。

圖4 不同時刻喉襯流場與固相材料的溫度分布Fig.4 Temperature distribution of flow field and solid material at various time

圖5為t=1、4、7 s時刻喉襯內壁溫度及燒蝕率的軸向分布。由圖5可知，喉襯內壁溫度及燒蝕率在喉部區域高，離開喉部區域則迅速下降，且峰值不在幾何喉部截面，而是位于喉部上游。這基本遵循內壁熱流密度的分布規律，見圖6。文獻［10－11］中，噴管燒蝕試驗結果也證實，燒蝕率峰值位于喉部上游，而不是在喉部截面。

圖5 不同時刻喉襯內壁溫度及燒蝕率軸向分布Fig.5 Axial distribution of temperature and erosion rate of throat inser wall at various time

圖6 不同時刻喉襯內壁熱流密度軸向分布Fig.6 Axial distribution of heat flux at various time

圖7為以燃氣總溫為定性溫度時的內壁瞬時對流換熱系數軸向分布。與圖7不同，由Bartz公式預示的內壁對流換熱系數，其峰值總在喉部截面。這是因為Bartz公式是基于管內旺盛湍流換熱試驗關聯式及一維等熵氣體動力學公式推導而來，未詳細考慮湍流邊界層發展及喉襯流道二維效應的影響［12－13］。實際上，從收斂段開始，邊界層將逐漸減薄，傳熱速率不斷增大;而在喉部上游某處，邊界層外緣將達到聲速，該處邊界層最薄，燃氣質量流率ρeue、熱流密度均達到峰值，材料的熱化學燒蝕率也最大;之后在擴張段，邊界層又將迅速增厚，傳熱速率也急劇下降［12－14］。

另外，由圖7也可知，隨傳熱過程的進行，燃氣與喉襯內壁的溫差逐漸減小，對流換熱逐漸減弱，對流換熱系數也逐漸降低，從中可見將流動與傳熱耦合計算的必要性。

圖7 不同時刻喉襯內壁對流換熱系數軸向分布Fig.7 Axial distribution of convection heat transfer coefficient at various time

3 結構分析

在發動機工作過程中，喉襯承受燃氣壓強及由傳熱燒蝕所致變溫載荷的作用，為了分析喉襯的結構完整性，需進行基于流場、熱分析結果的結構計算。由于結構位移、應變對流場、溫度場影響較小，因此忽略結構到流場、熱的耦合。

將FLUENT計算的網格及流場溫度場數據導入ANSYS結構分析軟件。刪去燃氣區域網格，將固相材料區域網格單元類型設置為PLANE 42;喉襯瞬態溫度場及內壁燃氣壓強作為載荷及邊界條件加載;施加位移約束條件，使喉襯外壁徑向位移為0，端面軸向位移為0。取彈性模量E=10.68 GPa，泊松比 μ =0.3，線膨脹系數 α =1.7×10－6/K。

由圖4所示溫度載荷及圖8所示壓強載荷計算了t=1、4、7 s時刻喉襯的應力分布。表2列出了各應力分量的最大、最小值，符號為正代表拉應力，符號為負代表壓應力。由于為空間軸對稱問題，因此僅有4個應力分量:徑向正應力σρ、軸向正應力σz、環向正應力σφ及作用在圓柱面而沿軸向的切應力τρz。由表2可知，各應力分量中最主要的是環向壓應力。

圖9為t=1、4、7 s時刻喉襯的環向應力分布。熱響應的深入，使喉襯溫度分布不均加劇，喉襯整體應力水平隨之升高，而喉部區域溫度梯度最大，相應地其應力值也最大。

圖8 喉襯內壁壓強載荷分布Fig.8 Pressure loading distribution on the inner wall

表2 不同時刻喉襯應力最大值Table 2 Maximum/minimum stress at various time

圖9 不同時刻喉襯環向應力分布Fig.9 Hoop stress distribution at various time

4 結論

(1)基于FLUENT流體計算軟件及ANSYS結構分析軟件，建立了固體火箭發動機喉襯流場及熱結構耦合分析模型，實現了流場與燒蝕傳熱的雙向耦合及流場、熱到結構的單向耦合。

(2)算例分析結果表明，喉襯溫度計算值與試驗值吻合較好;內壁對流換熱系數隨壁溫升高而降低;內壁熱流密度的峰值在喉部上游入口處，內壁溫度及燒蝕率也遵循相同的分布規律;喉部區域熱流密度大，燒蝕率大，溫度及其梯度高，應力水平也較高，環向壓應力是主要應力分量。

(3)該方法兼顧了分析精度與計算資源兩方面的因素，并充分發揮了2個軟件各自的優勢，有助于準確預示喉襯在燃氣環境中的燒蝕傳熱及熱結構行為。

［1］Lemoine L.Solid rocket nozzle thermostructural behavior［R］.AIAA 75-1339.

［2］吳景福.喉襯熱結構數值分析［J］.推進技術，1985，6(4):38-52.

［3］李耿，周為民，張鋼錘.噴管喉襯溫度場及熱應力計算分析［C］//中國航空學會航空動力分會火箭發動機專業委員會2006年學術年會文集，2006:183-188.

［4］Cozart A B，Shivakumar K N.Stress analysis of a 3-D braided composite ablative nozzle［R］.AIAA 99-1277.

［5］付鵬，蹇澤群，張鋼錘，等.發動機噴管喉襯燒蝕及熱結構工程計算［J］.固體火箭技術，2005，28(1):15-19.

［6］Swaminathan V，Setty S，Ramamadhav B S，et al.A method for the evaluation of temperature profiles and surface temperature in rocket motor nozzles with varying heat transfer coefficients［R］.AIAA 75-1347.

［7］陶文銓.數值傳熱學［M］.西安:西安交通大學出版社，2001:485-486.

［8］Thakre P，Yang V.Chemical erosion of carbon-carbon/graphite nozzles in solid-propellant rocket motors［J］.Journal of Propulsion and Power，2008，24(4):822-833.

［9］Gordon S，McBride B J.Computer program for calculation of complex chemical equilibrium compositions and applications［R］.NASA RP-1311，1994.

［10］陳博，張立同，成來飛，等.3D C/SiC復合材料噴管在小型固體火箭發動機中的燒蝕規律研究［J］.無機材料學報，2008，23(5):938-944.

［11］陳博，張立同，成來飛，等.燃氣發生器條件下穿刺C/C復合材料噴管的燒蝕性能研究［J］.無機材料學報，2008，23(6):1159-1164.

［12］Back L H，Massier P F，Gier H L.Convective heat transfer in a convergent-divergent nozzle［J］.International Journal of Heat and Mass Transfer，1964，7:549-568.

［13］Keizers H L J，Veraar R G.Analysis of heat transfer in nozzles［R］.AIAA 96-3290.

［14］Williams F A，Huang N C，Barrere M.固體推進劑火箭發動機的基本問題(上冊)［M］.京固群，譯.北京:國防工業出版社，1976:122-124.