液氧/甲烷發(fā)動機再生冷卻和膜冷卻傳熱數(shù)值研究

蘇 展，高玉閃，張曉光，邢理想，張 航

(西安航天動力研究所，西安 710100)

1 引言

在低溫可重復使用液體火箭發(fā)動機工作過程中，推進劑燃燒會在推力室產(chǎn)生極高的熱流，釋放出大量的熱量，產(chǎn)生極高的溫度和壓力，使推力室壁面中產(chǎn)生巨大的溫度梯度，壁溫急劇升高，且推力室壁面不同位置的溫度不同。為了保證結(jié)構(gòu)安全，在設(shè)計推力室結(jié)構(gòu)時，必須設(shè)計有效的冷卻系統(tǒng)以降低推力室壁面溫度，防止壁面溫度過高造成結(jié)構(gòu)損壞。常見的液體火箭發(fā)動機推力室冷卻系統(tǒng)包括再生冷卻和膜冷卻，工程上常使用多種冷卻方式相結(jié)合的方法共同進行推力室熱防護。

國內(nèi)對再生冷卻進行了大量的流場仿真，建立了多種二維和三維仿真模型。楊成驍?shù)柔槍σ后w火箭發(fā)動機分別建立了二維和三維仿真模型，其中二維模型的仿真精度較低，但計算效率較高，適用于發(fā)動機方案篩選和優(yōu)化設(shè)計；胡江玉等通過仿真分析了不同構(gòu)型面板對發(fā)動機再生冷卻的影響，為液體火箭發(fā)動機再生冷卻通道的設(shè)計提供了基礎(chǔ)。但是，流場仿真模型較為復雜，且對于較為復雜的再生冷卻模型，計算效率較低。因此在發(fā)動機設(shè)計階段，有必要采用計算效率較高的模型計算冷卻套的溫升和壓降。

在發(fā)動機系統(tǒng)仿真中，對于再生冷卻傳熱，工程上通常采用巴茲半經(jīng)驗公式計算分析推力室軸向的熱流密度、冷卻劑溫升和壓降等關(guān)鍵參數(shù)。針對甲烷低溫推進劑發(fā)動機，冷卻劑臨界溫度低，傳熱情況復雜，物性參數(shù)變化較為劇烈，巴茲公式已不再適用。北京航天動力研究所對此進行了大量的研究，利用修正的巴茲公式計算了氫氧發(fā)動機和液氧甲烷發(fā)動機冷卻套溫升，計算結(jié)果與試車結(jié)果偏差不大。但是對于膜冷卻的推力室，其邊區(qū)傳熱情況更為復雜，且膜冷卻傳熱對冷卻劑物性參數(shù)影響很大，普通的再生冷卻模型已不能滿足推力室傳熱分析的計算精度。

為更好地計算再生冷卻和膜冷卻推力室冷卻套溫升和壓降，本文通過分析冷卻套傳熱過程，考慮液膜冷卻對邊區(qū)傳熱的影響，對低溫推進劑冷卻套傳熱模型進行修正，并對比使用該模型的計算結(jié)果與發(fā)動機試驗數(shù)據(jù)，驗證該模型對甲烷推進劑傳熱計算的準確性和適應(yīng)性。

2 推力室再生冷卻模型

推力室再生冷卻傳熱過程可以通過隔著多層隔層的兩股流體間的換熱描述，采用無限大平面模型，計算兩股流體間的對流換熱。由于推力室外壁與環(huán)境溫度相差不大，可忽略外壁與環(huán)境間的對流換熱，且外壁兩側(cè)溫度差較小，忽略外壁間的熱傳導。由于燃氣輻射換熱量較小，忽略燃氣輻射換熱和推力室外壁與環(huán)境的換熱對計算結(jié)果影響不大。因此，再生冷卻分段傳熱過程模型如圖1所示，每一段中高溫燃氣將熱量通過邊界層傳遞到燃燒室內(nèi)壁，通過內(nèi)壁的導熱將熱量傳遞到冷卻劑邊界層而后被流動著的冷卻劑帶走。其一維穩(wěn)態(tài)傳熱關(guān)系式如式(1)所示。

圖1 冷卻套傳熱示意圖Fig.1 Schematic diagram of one-dimensional heat transfer model of cooling jacket

其中，為燃氣溫度，為燃氣側(cè)壁溫，為冷卻劑側(cè)壁溫，為冷卻劑溫度，為燃氣側(cè)對流換熱系數(shù)，為室壁熱導率，為室壁厚度，為冷卻劑側(cè)對流換熱系數(shù)。

2.1 燃氣側(cè)對流換熱

燃氣與推力室內(nèi)壁之間的對流傳熱屬于強迫對流傳熱。高溫燃氣高速流動會在內(nèi)壁面形成湍流邊界層，推力室燃燒情況較為復雜，邊界層內(nèi)部燃氣組分和流動狀態(tài)復雜。而推力室燃氣側(cè)對流換熱情況與邊區(qū)燃氣的狀態(tài)有關(guān)，因此，推力室湍流邊界層內(nèi)部傳熱情況較為復雜，通過理論計算的方式得到對流換熱系數(shù)的方法不再適用，一般使用半經(jīng)驗公式進行計算。最常使用的是巴茲公式，如式(2)～(3)所示。

燃氣的動力粘度和普朗特數(shù)根據(jù)燃氣的其他熱力參數(shù)近似計算，在喉部附近換熱系數(shù)可以采用喉部曲率半徑修正。

2.2 推力室內(nèi)壁間的熱傳導

在內(nèi)壁中，熱能是通過熱傳導的方式進行傳遞，根據(jù)傅里葉導熱定律，通過內(nèi)壁的熱流為式(4):

式中，為推力室壁熱傳導熱流密度，為推力室壁熱導率，為冷卻劑側(cè)壁溫。

2.3 冷卻劑側(cè)對流傳熱

由于甲烷冷卻劑的臨界溫度低，在冷卻通道中，隨著溫度的升高，冷卻劑由亞臨界狀態(tài)逐漸向超臨界狀態(tài)變化，熱量經(jīng)由冷卻劑蒸氣膜邊界層傳遞，基本關(guān)系可由式(5)表示:

式中，為冷卻劑側(cè)對流換熱熱流密度，為肋條傳熱修正系數(shù)，為冷卻劑側(cè)對流換熱系數(shù)，為冷卻劑溫度。

肋條傳熱修正系數(shù)計算方法如式(6)、(7)所示。

其中，為冷卻道槽寬，為肋高。

對于甲烷，冷卻劑側(cè)對流換熱系數(shù)采用經(jīng)驗公式(8)，該公式適用范圍約為90～500 K:

2.4 壓降

壓降分為沿程摩擦損失和動量損失，其中，延程摩擦損失常用式(9)、(10)計算:

式中，為沿程損失系數(shù)，為通道表面積，為通道直徑，為冷卻劑密度，為冷卻劑流速。

2.5 物性和數(shù)值計算方法

燃氣的物性由液氧甲烷熱力學計算獲取。采用自定義函數(shù)處理冷卻劑甲烷物性變化與溫度和壓力的關(guān)系，甲烷物性由REFPROP軟件計算獲得，能夠較為準確地計算出甲烷在跨臨界狀態(tài)時，其物性參數(shù)產(chǎn)生的突變。

推力室冷卻套一維計算模型采用集中參數(shù)法，但由于冷卻套型面較為復雜，且推進劑物性參數(shù)沿著推力室軸向變化很大，無法對整個推力室冷卻套采用集中參數(shù)法，因此采用分段計算的形式。根據(jù)推力室的幾何形狀和流動特點，對推力室進行分段，在段內(nèi)采用集中參數(shù)法，將上一段出口參數(shù)傳遞給下一段作為其入口參數(shù)。

3 液膜、氣膜冷卻修正

冷卻劑進入推力室，與高溫燃氣換熱，形成液膜冷卻區(qū)。隨著冷卻劑溫度的升高，冷卻劑逐漸氣化，形成氣膜冷卻區(qū)。液膜和氣膜冷卻區(qū)的換熱主要包括:燃氣與液膜(氣膜)之間的對流換熱、液膜(氣膜)與壁面之間的對流換熱、壁面內(nèi)的軸向?qū)岷捅诿媾c冷卻劑的對流換熱，如圖2所示。

圖2 膜冷卻傳熱示意圖Fig.2 Schematic diagram of heat transfer calculation model for film cooling

3.1 液膜冷卻區(qū)

液膜長度與液體冷卻劑的熱容和蒸發(fā)潛熱成正比。在本文中，由于室壓已經(jīng)遠遠高于甲烷的臨界壓力(4.64 MPa)，此時甲烷不存在氣液兩相共存的沸騰狀態(tài)，當被加熱到臨界溫度后直接過渡到氣態(tài)，即不存在蒸發(fā)潛熱，因此液膜區(qū)只存在液膜升溫段，不存在液膜蒸發(fā)段。認為甲烷液膜升溫到臨界溫度(190 K)即為氣膜冷卻的起始點。

液膜冷卻效率是液膜雷諾數(shù)的函數(shù)，如圖3所示。

圖3 液膜冷卻效率與雷諾數(shù)的關(guān)系Fig.3 Relationship between liquid membrane cooling efficiency and liquid membrane Reynolds number

液膜熱平衡方程為式(11)～(12):

3.2 氣膜冷卻區(qū)

氣膜效率計算方法如式(13)所示:

絕熱壁溫為式(14):

式中，為燃氣靜溫，為氣膜初始溫度(甲烷臨界溫度)。

壁面熱平衡方程為式(15):

式中，為氣膜與壁面之間的對流換熱熱流密度。

4 結(jié)果和討論

4.1 算例驗證

本文采用半經(jīng)驗公式和數(shù)值計算相結(jié)合的方法，具有一定的近似性。為了檢驗該方法的正確性和合理性，對10 t液氧甲烷發(fā)動機開展了額定工況和80%工況下的推力室冷卻傳熱計算。10 t液氧甲烷發(fā)動機冷卻液流向示意圖如圖4所示。發(fā)動機冷卻液從噴管前端流入冷卻套，流至推力室收擴段后折返，一小部分甲烷在推力室喉部前端流入推力室形成膜冷卻，剩余的甲烷流至推力室頭部冷卻套出口后進入推力室燃燒。

圖4 冷卻液流向示意圖Fig.4 Diagram of coolant flow

運用上述仿真模型，分別計算得到了100%工況和80%工況下的冷卻套溫升和壓降，計算結(jié)果如表1所示。由表可知，該計算模型仿真計算的不同工況下的冷卻套溫升比試車數(shù)據(jù)偏小，最大偏差為-3.7%。冷卻套壓降計算結(jié)果與試車數(shù)據(jù)相比，偏差較大，最大偏差為3.4%。其計算精度可靠，滿足工程運用需求。

表1 仿真結(jié)果與試驗結(jié)果對比Table 1 Comparison between simulation results and test data

圖5為2個工況下的熱流密度、氣壁溫、液壁溫、甲烷溫度和壓力的軸向分布。由圖可知，推力室冷卻套熱流密度在推力室頭部逐漸上升。在推力室收縮段，由于推力室直徑逐漸減小，換熱面積逐漸減小，熱流密度逐漸增大，因此，需要在推力室喉部前設(shè)置液膜冷卻，降低推力室的氣壁溫和液壁溫。液膜冷卻起始處，熱流密度迅速降低，隨著液膜溫度升高，液膜冷卻區(qū)逐漸向氣膜冷卻區(qū)過渡，熱流密度逐漸增大，到推力室喉部熱流密度最大，推力室氣壁溫和液壁溫達到最大。至噴管擴張段時熱流密度急劇下降，氣壁溫和液壁溫逐漸降低。冷卻通道的壓降受冷卻劑流動狀態(tài)影響，在冷卻套入口段，溫度較低，甲烷為液態(tài)，密度較高，流速較低，流阻較小。在推力室喉部區(qū)域，甲烷溫度較高，進入超臨界狀態(tài)，密度變化較大，且喉部流速最大，流阻最大，壓力迅速降低。隨著溫度的升高，甲烷密度逐漸減小，流速逐漸上升，流阻逐漸增大，壓力快速降低，在出口處壓力最低。

圖5 不同工況下冷卻套性質(zhì)延軸向分布圖Fig.5 Axial distribution of cooling jacket system parameters under different operating conditions

在2種不同的工況下，熱流密度分布趨勢基本一致，推力室工況越高，燃氣溫度越大，熱流密度越大，導致冷卻套氣壁溫和液壁溫急劇增大，對換熱的要求也會更高。2種工況下的甲烷溫度、壓降分布趨勢也基本一致，推力室工況越低，冷卻劑溫升越小，密度越大，冷卻套壓降越小。因此，該模型計算結(jié)果可以較好地模擬出推力室冷卻套的基本換熱情況，滿足發(fā)動機在變工況下的使用。

4.2 冷卻劑流量的影響

圖6是甲烷流量改變時推力室壁面溫度和熱流密度的軸向分布。由圖可得，隨著冷卻劑流量的增加，熱流密度增大，推力室換熱效果增強。推力室氣壁溫和液壁溫均降低，推力室壁面溫度最高點出現(xiàn)在推力室喉部，這是因為在該點熱流密度最大。當冷卻劑流量為6 kg/s時，推力室喉部氣壁溫為751.3 K，當甲烷流量增大至7 kg/s時，推力室喉部氣壁溫為728.6 K，推力室冷卻效果增強。

圖6 壁面溫度在不同甲烷流量下延軸向的分布Fig.6 Wall temperature distribution under different coolant flow

由此可知，在系統(tǒng)允許的條件下，應(yīng)盡可能增大冷卻劑流量，從而增大推力室冷卻效果，降低推力室壁溫，避免推力室壁面因為局部高溫環(huán)境產(chǎn)生結(jié)構(gòu)性損傷。

4.3 推力室室壓的影響

圖7是推力室室壓改變時推力室壁面溫度和熱流密度的軸向分布。由圖可得，隨著推力室室壓增大，推力室壁面溫度增大。推力室喉部為冷卻套壁溫最高點。這是由于推力室室壓越高，發(fā)動機工況越高，冷卻套熱流密度越高。因此，當發(fā)動機推力室室壓過高時，應(yīng)采取多種冷卻措施降低推力室壁溫，以防壁溫過高造成結(jié)構(gòu)損壞。

圖7 不同推力室室壓下的計算結(jié)果Fig.7 Calculation results under different chamber pressures

4.4 冷卻套壁面厚度的影響

不同壁面厚度對推力室壁溫的影響見圖8。所選取的冷卻套壁厚分別為1，2，3 mm。由圖可得，當推力室冷卻套壁厚增加時，推力室氣壁溫逐漸升高，而推力室液壁溫逐漸降低，降低幅度小于氣壁溫升高幅度。當壁厚由1 mm增大到3 mm時，推力室喉部氣壁溫由747 K上升到877 K，溫度上升劇烈。若壁厚繼續(xù)增大，則氣壁溫過高，有燒穿推力室壁的風險。

圖8 不同壁厚下的計算結(jié)果Fig.8 Calculation results of different thrust chamber wall thickness

由此可知，在設(shè)計范圍內(nèi)，推力室冷卻套壁面厚度應(yīng)越薄越好，這樣能夠降低推力室燃氣壁面溫度，防止造成推力室壁面結(jié)構(gòu)性損壞。

5 結(jié)論

1)本文模型計算結(jié)果能夠與試車數(shù)據(jù)較好地吻合，冷卻套溫升和壓降與試車數(shù)據(jù)最大偏差不超過4%，滿足工程計算精度要求，反應(yīng)了實際情況冷卻套的換熱規(guī)律。

2)冷卻劑流量越大，熱流密度越大，換熱效果越好，同時推力室壁面溫度越低，冷卻套冷卻效果越好；在設(shè)計允許范圍內(nèi)，推力室壁面厚度越薄，冷卻效果越好。

3)推力室室壓越大，推力室壁面溫度越高，當發(fā)動機由低工況向高工況轉(zhuǎn)變時，應(yīng)采取適當手段增強冷卻套的冷卻效果。