郝 雯，馬聰慧，邵躍躍

(中國空空導(dǎo)彈研究院，河南 洛陽 471009)

0 引 言

在固體火箭發(fā)動機(jī)中，通常會發(fā)生熱傳遞和質(zhì)量傳遞的綜合過程，該過程伴隨有吸熱或放熱的化學(xué)反應(yīng)。在固體火箭發(fā)動機(jī)噴管壁的熱交換過程中，同時存在著導(dǎo)熱、熱對流和熱輻射這三種基本熱量傳遞方式。火箭推進(jìn)劑在燃燒時，燃?xì)馔ㄟ^熱對流、熱輻射和凝聚相微粒的直接接觸導(dǎo)熱，將熱量傳遞給噴管壁內(nèi)表面;被加熱的噴管壁以熱傳導(dǎo)的方式，將熱量由內(nèi)表面經(jīng)過噴管壁向外表面?zhèn)鬟f;通過熱對流和熱輻射，熱量由外表面向周圍空間散失［1－4］。

固體火箭發(fā)動機(jī)噴管的傳熱過程十分復(fù)雜，包括熱傳導(dǎo)、熱輻射、熱對流。實際上，發(fā)動機(jī)噴管的熱狀態(tài)除與噴管初始條件相關(guān)之外，還與發(fā)動機(jī)噴管的結(jié)構(gòu)、材料物性參數(shù)等也有關(guān)系［5－6］。實際計算過程中，由于無法提前給出噴管壁面的熱邊界條件，同時還要考慮換熱中的輻射效應(yīng)，所以，噴管傳熱計算必須當(dāng)作是耦合的傳熱問題來進(jìn)行求解。

高速旋轉(zhuǎn)的工況會縮短固體火箭發(fā)動機(jī)的燃燒時間，增加推力與壓強(qiáng)，對發(fā)動機(jī)的傳熱造成嚴(yán)重的影響［7］，此外高速旋轉(zhuǎn)所致的強(qiáng)旋流動現(xiàn)象也對發(fā)動機(jī)結(jié)構(gòu)的熱防護(hù)帶來負(fù)面效應(yīng)。因此，計算高速旋轉(zhuǎn)工況下的固體火箭發(fā)動機(jī)噴管受熱狀態(tài)是十分有意義的。

1 噴管喉部換熱數(shù)值模型

本文建立了噴管流固熱耦合的換熱模型:先進(jìn)行發(fā)動機(jī)噴管內(nèi)流場數(shù)值模擬，之后將高精度的計算結(jié)果與噴管固壁的溫度場計算相耦合，其中考慮固壁中的熱傳導(dǎo)、內(nèi)流場與噴管壁面的湍流換熱以及包含吸收－發(fā)射性氣體介質(zhì)的熱輻射交換，并給予噴管構(gòu)型將該換熱模型應(yīng)用與軸對稱噴管的流場與溫度場的解算?？紤]噴管構(gòu)型的軸對稱性，為了減少數(shù)值仿真過程中的計算量，建立二維軸對稱數(shù)值模型。

二維軸對稱［8－9］質(zhì)量守恒方程如下:

上述方程是二維軸對稱的質(zhì)量守恒方程，該方程適用范圍包括:可壓和不可壓流動。式(1)中x是軸向坐標(biāo)，r 是徑向坐標(biāo)，u 和v 分別是軸向和徑向的速度分量，源項Sm是稀疏相增加到連續(xù)相中的質(zhì)量或其他加質(zhì)流動中增加的質(zhì)量。

對于二維軸對稱流場，建立考慮旋流在內(nèi)的軸向和徑向的動量守恒方程，如下:

噴管喉部換熱模型的能量守恒方程包含流體換熱、流固耦合換熱以及固體內(nèi)部的熱傳導(dǎo)等，流固耦合界面處的輻射換熱以及湍流邊界層內(nèi)的粘性換熱對固體內(nèi)部的熱結(jié)構(gòu)分布影響較大［10］，必須加以建模分析。能量守恒方程如下:

式中:keff= k + ki為有效熱導(dǎo)率，ki為湍流引起的導(dǎo)熱率，由采用的湍流模型以及壁面函數(shù)確定［11］。式(4)右端的前三項分別表示由于熱傳導(dǎo)、組分?jǐn)U散、粘性耗散引起的能量轉(zhuǎn)移，源項是由于化學(xué)反應(yīng)、輻射等其他因素引起的熱源。

在固體區(qū)域，能量控制方程如下:

式中，速度v 表示固體區(qū)域由于旋轉(zhuǎn)或平移等運(yùn)動的速度。式(5)右端第一項表示熱傳導(dǎo)引起的熱流，第二項即源項表示固體區(qū)域的內(nèi)熱源。

式(1)～(5)是二維軸對稱噴管喉部流固耦合換熱數(shù)值仿真模型的基本控制方程，通過求解方程組可以獲取噴管喉部計算域內(nèi)的溫度場、壓力、速度場等分布，為噴管喉部換熱研究提供數(shù)據(jù)支持。

2 發(fā)動機(jī)噴管受熱狀態(tài)數(shù)值計算

對噴管在靜態(tài)和旋轉(zhuǎn)條件下的受熱狀態(tài)做計算研究，對轉(zhuǎn)速為0 r/min，4 000 r/min，8 000 r/min 時進(jìn)行流固熱耦合數(shù)值模擬。

2.1 轉(zhuǎn)速為0 r/min 時噴管受熱狀態(tài)計算

取燃燒室平均壓強(qiáng)7.19 MPa 為入口條件，45#鋼給定密度7. 85 ×103 kg/m3，比熱給定450 J/(kg·K)，導(dǎo)熱系數(shù)給定35 W/(m·K)，進(jìn)行非穩(wěn)態(tài)傳熱計算，時間步長取0.1 s，計算5 s 內(nèi)噴管中燃?xì)獾牧鲃优c換熱。

不同時刻噴管溫度云圖如圖1 所示。

圖1 0 r/min 時不同時刻溫度云圖

隨著發(fā)動機(jī)工作時間的延續(xù)，噴管壁面的溫度處于持續(xù)上升狀態(tài)，噴管壁面的溫度傳遞呈二維特性，靠近噴管內(nèi)壁面溫度相對較高。

表1 是轉(zhuǎn)速為0 r/min 時，發(fā)動機(jī)在1.8 s，3.8 s，4.8 s 時不同位置的熱流密度值。分別取三個觀測點，噴管收斂段(x =40 mm)、噴管喉部(x=85 mm)、噴管擴(kuò)張段(x=100 mm)處，得到該點處熱流密度的值。

表1 0 r/min 時 不同時刻壁面總熱流密度表(W/m2)

由上表可知，噴管喉部區(qū)域總熱流密度最大，表示噴管喉部換熱最為強(qiáng)烈。因此，喉部前端部分區(qū)域溫度最高，收斂段溫度逐漸降低。而不同時刻t=1.8 s，t=3.8 s，t =4.8 s 時，燃?xì)鈧?cè)對噴管壁面熱交換的總熱流密度逐漸降低，發(fā)動機(jī)噴管喉部位置的總熱流密度最大，換熱最為強(qiáng)烈，當(dāng)發(fā)動機(jī)工作1.8 s 時，可達(dá)8 400 000 W/m2以上。由于噴管喉部燃?xì)饬鲃幼兓瘎×遥虼丝偀崃髅芏茸兓?，在噴管的收斂段和擴(kuò)張段，熱流密度均比喉部熱流密度要低。

由計算區(qū)域溫度云圖可知，噴管壁面的溫度傳遞呈二維特性，靠近噴管內(nèi)壁面溫度相對較高，噴管喉部區(qū)域換熱最為強(qiáng)烈，這也與不同時刻壁面總熱流密度示意圖相對應(yīng)。喉部區(qū)域熱流密度最大。因此，喉部前端部分區(qū)域溫度最高，收斂段溫度逐漸降低。

2.2 轉(zhuǎn)速為4 000 r/min 時噴管受熱狀態(tài)計算

取燃燒室平均壓強(qiáng)25. 54 MPa 為入口條件，其他條件與2.1 節(jié)同，計算5 s 內(nèi)噴管中燃?xì)獾牧鲃雍蛽Q熱。

不同時刻噴管溫度云圖如圖2 所示。

表2 是轉(zhuǎn)速為4 000 r/min 時，發(fā)動機(jī)在0.5 s，1.5 s，3.0 s 時不同位置的熱流密度值。分別取噴管收斂段(x=40 mm)、噴管喉部(x=85 mm)、噴管擴(kuò)張段(x=100 mm)三個點計算熱流密度的值。

由于在高速旋轉(zhuǎn)狀態(tài)下，發(fā)動機(jī)噴管內(nèi)壓強(qiáng)增加，導(dǎo)致噴管內(nèi)壁面的熱流密度增加，壁面溫度升高。熱流密度的變化規(guī)律與0 r/min 時大致相同，在噴管喉部達(dá)到峰值。在噴管的收斂段和擴(kuò)張段，隨著到噴管喉部距離逐漸增加，熱流密度逐漸減小。

圖2 4 000 r/min 時不同時刻溫度云圖

表2 4 000 r/min 時不同時刻壁面總熱流密度表(W/m2)

2.3 轉(zhuǎn)速為8 000 r/min 時噴管受熱狀態(tài)計算

取燃燒室平均壓強(qiáng)41.92 MPa 為入口條件，其他條件與2.1 節(jié)同，計算5 s 內(nèi)噴管中燃?xì)獾牧鲃优c換熱。

不同時刻噴管溫度云圖如圖3 所示。

圖3 8 000 r/min 時不同時刻溫度云圖

表3 是8 000 r/min 時，發(fā)動機(jī)在0.5 s，1.5 s，3.5 s 時不同位置的熱流密度值。分別取噴管收斂段(x=40 mm)、噴管喉部(x =85 mm)、噴管擴(kuò)張段(x=100 mm)三個觀測點，計算熱流密度的值。

表3 8 000 r/min 時 不同時刻壁面總熱流密度表(W/m2)

在8 000 r/min 的轉(zhuǎn)速下，由于壓強(qiáng)，溫度等綜合作用，溫度傳遞與熱流密度的變化與前兩節(jié)變化相同。

2.4 對比分析

不同時刻的溫度對比云圖于前幾節(jié)所示，噴管壁面收斂段、直線段以及擴(kuò)張段不同位置的溫度均逐漸上升。不同時刻下噴管內(nèi)壁面溫度均在噴管喉部前端達(dá)到最大值，在噴管收斂段以及擴(kuò)張段壁面溫度均呈現(xiàn)下降趨勢。

對比分析同一工況下溫度場分布圖可知，噴管喉部流場流動區(qū)域的溫度在不同時刻變化不大，而噴管內(nèi)壁面的溫度隨著發(fā)動機(jī)的熱交換程度增大以及發(fā)動機(jī)工作時間的延續(xù)而逐漸增加，因此發(fā)動機(jī)喉部燃?xì)鈧?cè)對噴管喉部內(nèi)壁面的熱交換程度隨時間逐漸改變。

由表1 ～3 可知，發(fā)動機(jī)隨著工作時間的延續(xù)，燃?xì)鈧?cè)對噴管內(nèi)壁面熱交換的總熱流密度逐漸降低，發(fā)動機(jī)噴管喉部前端位置總熱流密度最大，換熱最為強(qiáng)烈。

燃?xì)鈧?cè)對噴管壁面熱交換的對流換熱系數(shù)主要由燃?xì)馕镄詤?shù)影響［12］，基本不隨發(fā)動機(jī)工作時間改變而改變，在噴管喉部前端處，分子活動最為劇烈，導(dǎo)致對流換熱系數(shù)最大，對流換熱最為強(qiáng)烈。在噴管擴(kuò)張段，隨著發(fā)動機(jī)工作時間的延續(xù)，流場流動逐漸發(fā)展，對流換熱系數(shù)略有上升。但與噴管喉部換熱系數(shù)相比要低于喉部數(shù)值。

高速旋轉(zhuǎn)對發(fā)動機(jī)壁面換熱的影響總趨勢是隨著旋轉(zhuǎn)過載的加大而增加。發(fā)動機(jī)燃?xì)馇邢蛩俣蕊@著增加，單位時間內(nèi)流過壁面的燃?xì)夥肿訑?shù)增多，對流換熱系數(shù)增大，導(dǎo)致對流換熱加劇。在燃燒室及噴管區(qū)域，均產(chǎn)生了燃?xì)鉁u旋，燃?xì)馑哂械牟糠謩幽茉跍u旋中逐漸耗散轉(zhuǎn)變?yōu)闊崮?，壁面的傳熱加?qiáng)。

圖4 為文獻(xiàn)［13］中的高速旋轉(zhuǎn)發(fā)動機(jī)試驗結(jié)果，由文獻(xiàn)可知，在高速旋轉(zhuǎn)條件下，噴管內(nèi)表面型面有明顯的燒蝕，特別是喉襯鑲嵌接縫處。噴管內(nèi)燃?xì)饬魉俑?、單位截面積的質(zhì)量流率大，使得噴管壁所受到的高溫高壓燃?xì)獾募訜嶙饔檬謬?yán)重，從而造成噴管材料的嚴(yán)重?zé)g。燒蝕現(xiàn)象隨著轉(zhuǎn)速的提高而愈加嚴(yán)重，與數(shù)值計算結(jié)果相對應(yīng)。

圖4 高速旋轉(zhuǎn)實驗發(fā)動機(jī)喉部照片

結(jié)合發(fā)動機(jī)高速旋轉(zhuǎn)時的內(nèi)流場流線示意圖5分析。由于實際固體推進(jìn)劑在燃燒過程中，將產(chǎn)生一定量的凝聚相微粒。這些凝聚相微粒的存在，對傳熱產(chǎn)生影響。高溫凝聚相微粒的導(dǎo)熱性較氣相高，對壁面通過撞擊直接接觸進(jìn)行熱傳導(dǎo)。

圖5 10 000 r/min 時發(fā)動機(jī)流場流線圖

同時，回旋渦流加劇了高溫燃?xì)庵心嗔Ｗ拥某练e和積聚，這些粒子凝固后的積屑和積瘤導(dǎo)致發(fā)動機(jī)內(nèi)部噴管熱傳遞效應(yīng)加劇，凝相粒子與高溫高壓燃?xì)猱a(chǎn)生的渦流共同沖刷該區(qū)域。

文獻(xiàn)［13］中噴管座螺紋連接處、噴管進(jìn)口端和出口端均有明顯的金屬流動痕跡和熔化。文獻(xiàn)［14］中發(fā)動機(jī)前封頭處被燒穿，也驗證了凝相粒子對壁面的傳熱影響，與數(shù)值模擬結(jié)果相符。

這些現(xiàn)象的揭示對固體火箭發(fā)動機(jī)的熱防護(hù)設(shè)計提出了更高的要求。

3 結(jié) 論

本文建立了二維非定常噴管壁面流固熱耦合模型，計算噴管受熱狀態(tài)、熱流密度等參數(shù)。計算結(jié)果表明，隨著發(fā)動機(jī)工作時間的延續(xù)，高溫燃?xì)鈱姽軆?nèi)壁面熱交換的總熱流密度逐漸降低，發(fā)動機(jī)噴管喉部前端位置總熱流密度最大，換熱最為強(qiáng)烈。

高速旋轉(zhuǎn)對發(fā)動機(jī)噴管壁面?zhèn)鳠嵊绊懙内厔菔请S著旋轉(zhuǎn)過載的增大而增大。發(fā)動機(jī)燃?xì)馇邢蛩俣蕊@著增加，單位時間內(nèi)流過壁面的燃?xì)夥肿訑?shù)增多，使對流換熱系數(shù)增大，導(dǎo)致對流換熱加劇。在噴管區(qū)域，產(chǎn)生了燃?xì)鉁u旋，燃?xì)馑哂械牟糠謩幽茉跍u旋中逐漸耗散轉(zhuǎn)變?yōu)闊崮?，壁面的傳熱加?qiáng)。高速旋轉(zhuǎn)對噴喉結(jié)構(gòu)強(qiáng)度與熱防護(hù)帶來嚴(yán)重影響。

［1］Cortopassi A C，Boyer E，Acharya R，et al. Design of a Solid Rocket Motor for Characterization of Submerged Nozzle Erosion［C］// 44th AIAA/ASME/SAE/ASEE Joint Propulsion Conference ＆ Exhibit，AIAA 2008－4889，2008.

［2］Klager K.The Interaction of the Efflux of Solid Propellants with Nozzle Materials［J］. Propellants and Explosives，1997，2(3):55－63.

［3］劉玉磊.燃?xì)舛媪鞴恬詈蟼鳠釘?shù)值分析［J］.航空兵器，2013 (3):41－43.

［4］王偉，王德升.噴管擴(kuò)張段絕熱層的燒蝕計算［J］.固體火箭技術(shù)，1999，22(3):16－19.

［5］Wilcox D C.Turbulence Modeling for CFD［M］.La Canada，California，DCW，1998. Industries，Inc，1998.

［6］Lien F S ，Leschziner M A. Assessment of Turbluent Transport Models Including Non－Linear RNG Eddy－Viscosity Formulation and Second－Moment Closure［J］.Computers and Fluids，1994，23(8):983－1004.

［7］Caveny L H. Extension to Analysis Ignition Transients of Segmented Solid Rocket Booster［R］. NASA－CR－150632，1978.

［8］陶文銓.數(shù)值傳熱學(xué)［M］.西安:西安交通大學(xué)出版社，1988.

［9］Patankar S V.傳熱與流體流動的數(shù)值計算［M］. 張政，譯.北京:科學(xué)出版社，1984.

［10］Golasfhani M. Computation of Two－Phase Viscous Flow in Solid Rocket Motors Using a Flux－Split Eulerian－Lagrangian Technique［R］.AIAA 1989－2178，1989.

［11］王福軍.計算流體動力學(xué)分析——CFD 軟件原理與應(yīng)用［M］.北京:清華大學(xué)出版社，2004.

［12］鄭亞，陳軍，鞠玉濤，等. 固體火箭發(fā)動機(jī)傳熱學(xué)［M］. 北京:北京航空航天大學(xué)出版社，2006.

［13］王棟，余陵，武曉松. 固體火箭發(fā)動機(jī)高速旋轉(zhuǎn)試驗研究［J］. 彈道學(xué)報，2004(4):87－91.

［14］邵愛民.大型固體發(fā)動機(jī)旋轉(zhuǎn)試車頭部熱防護(hù)工程分析［J］.固體火箭技術(shù)，1998(3):7－12.