小推力液體火箭發(fā)動機在軌熱分析

閆 波，張會強，王 兵

(1.清華大學 航天航空學院，北京 100084；2.中國運載火箭技術研究院 研究發(fā)展中心，北京 100076)

0 引言

空間探測是航天活動的一個熱點方向[1-2]，而空間探測離不開小推力發(fā)動機。小推力發(fā)動機的工作環(huán)境通常是太空，所以需要進一步考慮太空星體對發(fā)動機的影響，特別是太陽輻射對模型的非均勻性熱影響。當飛行器在空間預定軌道飛行時，由于飛行器相對于太陽和地球的位置和方向發(fā)生著不斷的變化，其外熱流和輻射交換也發(fā)生著改變，從而使得發(fā)動機遭受高低溫變化的影響。因此，對小推力液體火箭發(fā)動機進行空間在軌熱分析是十分有意義的。

本文利用ANSYS Workbench 有限元分析軟件，對小推力發(fā)動機推力室在地球同步軌道(GEO)的結構熱特性開展了分析研究，研究成果為小推力液體火箭發(fā)動機結構熱設計提供了依據。

1 GEO軌道外熱流計算方法

空間飛行器繞地球軌道運行時，受到太陽直接輻射、地球反照和地球的紅外輻射，這三部分通常稱為空間外熱流(簡稱外熱流)，屬于空間飛行器的外部熱環(huán)境[3-5]。

當空間飛行器在約200 km以上高度的低地球軌道(LEO)上飛行時，其接受的空間外熱流主要是太陽輻射、地球反照輻射和地球紅外輻射三種，對于GEO軌道，地球反射和紅外輻射可以忽略不計，空間飛行器接受的空間外熱流主要是太陽輻射[5]。以下只針對太陽輻射進行分析。

從LEO至GEO軌道的高度上，太陽光被認為是均勻的平行光束，太陽輻射強度為太陽常數S，目前規(guī)定S=(1 353±21)W/m2。

空間飛行器外表面任一微元面積dA上(如圖1)受到的太陽輻射外熱流為

dq1=ScosθdA

(1)

令

φ1=dq1/SdA

則

φ1=cosθ

式中φ1為太陽輻射角系數。任一面元的太陽輻射外熱流為

dq1=αsSφ1dA

(2)

式中αs為面元表面的太陽吸收率，太陽輻射角系數需要分析空間飛行器在空間的運動規(guī)律，以確定微元表面與地球、太陽的相對關系[5]。

圖1 太陽輻射熱流圖Fig.1 Heat flow diagram of solar radiation

綜上所述，本文對熱源的基本假定如下：

1)太陽輻射計算時采用的太陽常數為1 353 W/m2(熱流密度)；

2)只考慮太陽對發(fā)動機的直射和斜射，不考慮太陽的散射，太陽光為平行光；

3)不考慮地球的紅外輻射，地球反照；

4)忽略空間飛行器對發(fā)動機的輻射和導熱，不計其他行星的熱輻射；

5)外層空間是絕對黑體；

6)只考慮推力室自身遮擋光線，不考慮其他物體遮擋[6]。

2 有限元數值仿真模型

為了研究太陽輻射對推力室結構熱特性的影響，本文采用有限元法對深空工作環(huán)境下推力室穩(wěn)態(tài)工作時和不工作時分別進行了模擬計算[7-8]，兩種分析都采用三維穩(wěn)態(tài)熱分析模型。

2.1 推力室工作參數

本文主要研究對象是某小推力液體發(fā)動機，該型發(fā)動機相關參數如表1所示，發(fā)動機額定推力為445 N，N2O4為氧化劑，混肼-50或一甲基肼(MMH)為燃料，采用液膜加輻射組合的冷卻方式。發(fā)動機有較長的工作壽命，穩(wěn)態(tài)工作時間大于500 s。該發(fā)動機由推力室(噴注器、燃燒室、噴管)、各種閥門、調節(jié)器及機架等各部分組成[9]。

表1 發(fā)動機綜合參數Tab.1 Integrated parameters of engine

2.2 定義材料屬性及單元類型

推力室三維實體模型在Pro/E中創(chuàng)建，通過CAD和ANSYS Workbench的無縫連接將模型導入。建立三維實體模型后，需要定義推力室結構材料屬性，包括材料密度、比熱、導熱等參數。

單元類型選用ANSYS 軟件中的solid87，單元尺寸取為0.002 m，因為推力室中同時存在熱對流和熱輻射換熱方式，為了避免施加載荷的覆蓋情況出現，這里要在內表面和外表面各設定一個表面效應單元SURF152。在內表面把燃氣輻射熱流密度施加在實體單元上，把對流換熱施加在表面效應單元上，在外表面把太陽對推力室的輻射施加在表面效應單元上。

2.3 邊界條件與加載過程

對于小推力液體火箭發(fā)動機，由于受發(fā)動機結構本身及擠壓式供應系統(tǒng)供應壓力的限制，一般采用液膜冷卻結合輻射冷卻的方法。同時，其中包含的傳熱過程有燃氣對結構及結構和環(huán)境的輻射換熱、燃氣對結構及推進劑對結構的對流換熱和結構內部的熱傳導。小推力發(fā)動機的工作過程，涉及到三種基本傳熱方式，即熱傳導、熱對流和熱輻射。

加載載荷包括：模型初場溫度分布、工作環(huán)境溫度、燃氣輻射熱流密度和燃氣恢復溫度、流體與壁面間對流換熱系數等。本文全場發(fā)射率均設置為0.9，初場溫度分布均為298 K，工作環(huán)境溫度為4 K。利用巴茲公式和發(fā)動機的熱力計算數據，得到了推力室內的燃氣壁面表面對流換熱系數、燃氣恢復溫度和燃氣溫度的分布，如圖2和圖3中所示[10-11]。噴注面上的冷卻劑孔噴出冷卻液，在燃燒室前段形成冷卻液膜，液膜長度約為燃燒室的一半，液膜溫度約為430 K。因為液膜厚度很小，且其透射率較高，故忽略燃氣和壁面對液膜的輻射傳熱過程。噴注器前端面和噴管末端面絕熱處理。

圖2 燃氣(液膜)沿軸向溫度分布Fig.2 Distribution of gas(liquid film)temperature along the axial

圖3 燃氣(液膜)沿壁面對流換熱系數分布Fig.3 Distribution of gas(liquid film)convection heat transfer coefficient along the wall

外壁對太空的輻射采用軟件自帶輻射模型施加。在太空中，由于不同時刻太陽光線照射到推力室的位置不同，所以不同時刻施加到發(fā)動機表面的熱量也不同，直射時通過單元熱流密度為q=S，S為太陽常數，斜射時q=Scosθ，太陽輻射施加的載荷類型為熱流密度[6]。太陽光在坐標系中方向余弦為(cosα，cosβ，cosγ)，在ANSYS Workbench 下插入APDL 語言加載燃氣輻射、對流載荷和太陽輻射，將燃氣輻射熱流密度、燃氣恢復溫度和對流換熱系數設置為表格數組加載，把太陽光方向余弦設為參數變量，可根據發(fā)動機在太空中所處的方位設計不同的太陽輻射角度。本節(jié)假設太陽沿x軸垂直照射，即方向余弦為(-1，0，0)，如圖4所示。

圖4 太陽輻射方向Fig.4 Direction of solar radiation

噴注面有燃氣輻射，需要施加輻射邊界條件，本文計算模型中對噴注面輻射的燃氣溫度取為2 000 K是合適的。需要對噴注器管道施加對流換熱邊界條件，需要得到管道對流換熱系數和流體溫度對流換熱系數，如表2所示。

表2 噴注器管道對流換熱系數Tab.2 Convective heat transfer coefficient of injector pipe

推進劑入口溫度取290 K。但在施加邊界條件時，為簡便起見，各段流體溫度設置為平均溫度，氧化劑主噴孔、燃料主噴孔、冷卻劑孔內流體溫度取320 K[12]。噴注器的輻射和對流換熱采用軟件自帶模型施加。

推力室不工作工況只有推力室對太空的輻射和太陽對推力室的輻射，不計地球反照輻射、地球紅外輻射和其他行星的熱輻射。

3 數值仿真結果分析

3.1 發(fā)動機穩(wěn)態(tài)工作工況

經模擬計算得到有太陽輻射和沒有太陽輻射推力室室壁溫度分布規(guī)律，有太陽輻射推力室溫度分布云圖如圖5所示，圖6為沒有太陽輻射推力室穩(wěn)態(tài)工作時溫度分布。

圖5 有太陽輻射穩(wěn)態(tài)工作推力室壁溫分布Fig.5 Wall temperature distribution of thrust chamber with solar radiation in steady working

圖6 無太陽輻射穩(wěn)態(tài)工作推力室壁溫分布Fig.6 Wall temperature distribution of thrust chamber without solar radiation in steady working

有太陽輻射時選取如圖5中兩條曲線上的溫度進行分析，經過分析比較發(fā)現，曲線1溫度略高于曲線2溫度，最大溫差出現在噴管出口處，溫度相差10 K左右，其余位置溫度相差在3 K以下。對有太陽輻射和沒有太陽輻射推力室室壁溫度分布云圖比較可知，溫度分布和大小幾乎一致，噴管喉部溫差在1 K以下，噴管出口處相差最大，最大差值10 K左右。由此可以得出結論，太陽輻射對推力室穩(wěn)態(tài)工作時溫度有一定影響，但影響不大。

3.2 發(fā)動機不工作工況

推力室發(fā)動機不工作工況只有推力室對太空的輻射和太陽對推力室的輻射，模擬計算得到推力室溫度分布云圖如圖7所示。推力室最高溫度在太陽直射面噴管尾部為338.81 K，由于有熱傳導作用，溫度最小值不是太空溫度4 K而是267.27 K，最小值在推力室背光面的噴管尾部，推力室溫度沿負x軸方向逐漸遞減。

如圖4所示在推力室身部內壁面兩條母線上的不同位置選取了6個采樣點，分別位于x1=-0.074 m，x2=-0.022 m，x3=0 m，x4=0.070 m，x5=0.123 m，x6=0.186 m處。各處溫度變化曲線如圖8所示。從圖8中可以看出曲線1上的溫度變化趨勢為溫度先上升后有小幅度下降，過了喉部再次上升，噴管出口處達到最大值；曲線2上的溫度一直保持下降變化趨勢，噴管出口處達到最小值。從兩曲線溫度差值可以看出，曲線1和曲線2溫度差值在燃燒室和噴注器連接處最小，只相差4.15 K，然后溫差開始逐漸增加直到燃燒室結束，而后溫度差又開始減少，喉部溫度差為8.95 K，過了喉部溫度差開始逐漸增加，到達噴管尾部達到最大溫度差值71.41 K。

圖7 有太陽輻射發(fā)動機不工作時推力室壁溫分布Fig.7 Wall temperature distribution of thrust chamber with solar radiation under not working condition

圖8 推力室壁溫沿軸向變化曲線Fig.8 Changing curve of thrust chamber wall temperature along the axial

由此可以得出結論，當發(fā)動機穩(wěn)態(tài)工作時太陽輻射對推力室溫度有一定影響，但影響不大。而當發(fā)動機不工作時，由于小推力發(fā)動機工作環(huán)境是太空，發(fā)動機機體一半接受太陽輻射，溫度較高，一半面對深空冷環(huán)境，溫度較低，太陽輻射對模型的非均勻性影響較大。

4 結論

對考慮太陽輻射的深空工作環(huán)境下推力室穩(wěn)態(tài)工作時和發(fā)動機不工作時分別進行了模擬計算，太陽輻射沿x軸方向垂直照射推力室，兩種分析都采用三維穩(wěn)態(tài)熱分析模型。

1)當發(fā)動機穩(wěn)態(tài)工作時太陽輻射對推力室溫度有一定影響，但影響不大。有太陽輻射時推力室最大溫差出現在直射母線和背光母線噴管出口處，溫度相差10 K左右，其余位置溫度相差在3 K以下。有太陽輻射和無輻射室壁溫度分布和大小幾乎一致，噴管喉部溫差在1 K以下，噴管出口處相差最大，最大差值10 K左右。

2)發(fā)動機不工作工況只有推力室對太空的輻射和太陽對推力室的輻射，最高溫度在太陽直射面噴管尾部為338.81 K，由于有熱傳導作用，溫度最小值不是太空溫度4 K而是267.27 K，最小值在推力室背光面的噴管尾部，推力室溫度沿負x軸方向逐漸遞減。當發(fā)動機不工作時，由于小推力發(fā)動機工作環(huán)境是太空，發(fā)動機機體一半接受太陽輻射，溫度較高，一半面對深空冷環(huán)境，溫度較低，太陽輻射對模型溫度分布均勻性有一定影響。