李 煜，程廣益，李宗陽，竇怡彬，蔣君庭

（上海機電工程研究所，上海 201109）

0 引 言

噴管是固體火箭發(fā)動機的基本部件之一。噴管部件在高溫燃氣的作用下溫度急劇上升，在進行火箭發(fā)動機各部件結(jié)構(gòu)設(shè)計時，需要考慮各部件受熱狀態(tài)下的實際強度以保證發(fā)動機穩(wěn)定工作。此外，燃氣傳熱會引起各部件內(nèi)部溫度分布不均，從而引入熱應(yīng)力載荷的作用。因此，要進行合理的火箭發(fā)動機噴管結(jié)構(gòu)設(shè)計就必須考慮傳熱的影響。在兩相流動時，顆粒在噴管內(nèi)會不斷撞擊噴管內(nèi)壁，對噴管的熱防護、變形量以及應(yīng)力的大小和分布造成影響［1］。因此，研究噴管在兩相流動下的傳熱特性及熱固耦合有重要意義。

本文考慮噴管結(jié)構(gòu)（包括燒蝕層、喉襯、絕熱層以及殼體）的傳熱特性、噴管外壁與大氣之間的傳熱因素，分析噴管結(jié)構(gòu)對內(nèi)部溫度場的影響；利用顆粒軌道模型［2］，研究不同顆粒直徑D和不同負載比η時兩相流動對噴管結(jié)構(gòu)內(nèi)部的溫度場的影響；進行熱固耦合計算，研究純氣相和兩相流動時噴管的熱變形分布，分析兩相流動對熱變形的影響。

1 數(shù)值計算方法

計算時進行以下簡化：不考慮化學反應(yīng)；不考慮顆粒相的燃燒、蒸發(fā)、破碎過程，且各個工況均為穩(wěn)態(tài)計算；不考慮內(nèi)壁面的燒蝕；不考慮顆粒接觸傳熱；不考慮各層材料之間的接觸熱阻；燃氣為理想氣體。

1．1 計算模型

考慮到計算量較大，本文采用了二維軸對稱模型。噴管結(jié)構(gòu)如圖1所示。各結(jié)構(gòu)材料分別為：耐燒蝕層為石墨布；喉襯為碳碳復(fù)合材料；絕熱層為碳酚醛；殼體為45號鋼。各材料物性參數(shù)如表1所示［3－7］。

圖1 噴管模型Fig．1 The model of nozzle

表1 材料物性參數(shù)Tab．1 The physical property parameters of the materials

1．2 計算方法

對于高速可壓縮流場計算，通常使用密度基求解［8］，本文采用隱式Reo格式求解方法。目前對固體火箭發(fā)動機流場進行模擬時，湍流模型大多采用標準k－ε模型，可取得較高的精度。因此，本文主要采用標準k－ε模型來進行計算。

描述氣固兩相流動時，采用顆粒軌道模型，顆粒相的瞬時速度vp由式（1）積分確定。

式中：mp為顆粒相的質(zhì)量；vp為顆粒相速度；t為時間；Fp為拖曳力，表達式為

式中：ρ為氣相密度；dp為顆粒直徑；CD為拖曳系數(shù)，v為氣相速度。

顆粒的位置可通過式（3）計算。

式中：Xp是t＋Δt時刻顆粒的位置；為t時刻顆粒的位置，對式（1）和式（3）進行數(shù)值積分可以獲得顆粒的軌跡。

1．3 邊界條件

根據(jù)實際情況給定適當?shù)倪吔鐥l件。假定壁面無滑移，流體與噴管內(nèi)壁為耦合傳熱。噴管入口設(shè)為壓力入口，給定壓力入口總壓7 MPa，總溫3 000 K；噴管出口為壓力出口；出口和入口均給定湍流強度和水力直徑。噴管外壁面、殼體尾部端面、絕熱層尾部端面、燒蝕層尾部端面的對流換熱系數(shù)根據(jù)大空間自然對流換熱計算公式得到［9］，為10 W／（m2·K），發(fā)射率為0．8。殼體前部端面、絕熱層前部端面、燒蝕層前部端面均為絕熱面，相鄰固體材料之間為內(nèi)部壁面。外部自然來流溫度為300 K。顆粒相采用三氧化二鋁，密度為2 970 kg／m3，比熱容為897 J／（kg·K），導熱率為202．4 W／（m·K），計算選擇的顆粒直徑分別為1!m、10!m、20!m 和50!m，負載比分別為0．01、0．05、0．1、0．2 和0．3。

1．4 網(wǎng)格無關(guān)性驗證

不同數(shù)量的網(wǎng)格可能會對計算結(jié)果產(chǎn)生一定的影響，因而需要進行網(wǎng)格無關(guān)性驗證，確保網(wǎng)格對計算結(jié)果不會產(chǎn)生明顯影響。整體網(wǎng)格及壁面邊界層附近的網(wǎng)格如圖2所示。本文選取的3種網(wǎng)格數(shù)量及相關(guān)的純氣相計算結(jié)果如表2所示。

圖2模型網(wǎng)格Fig．2 The mesh of model

表2 網(wǎng)格數(shù)量及計算結(jié)果Tab．2 Grid number of meshes and the calculation results

經(jīng)對比發(fā)現(xiàn)，3種網(wǎng)格對計算結(jié)果的影響甚微。為了節(jié)省計算資源，本文選擇Mesh1作為計算網(wǎng)格。

2 流場計算結(jié)果與分析

2．1 純氣相流場分析

圖3所示為流固交界面、耐燒蝕層和絕熱層交界面、絕熱層和殼體交界面以及殼體外壁面的溫度分布。圖4為文獻［3］中流固交界面上的溫度分布。對比圖3和圖4可以發(fā)現(xiàn)，由于文獻中的燃面不斷推移，因而在圖中左側(cè)有較長的平直段，同時文獻中的噴管擴張段較短，因而溫度下降較為明顯，但在噴管的收縮段及喉部的溫度分布和文獻中的噴管內(nèi)壁面溫度分布一致，從而驗證了本文計算結(jié)果的準確性。

由圖3可以明顯地看出：絕熱層碳酚醛的導熱率較小，對殼體起到了很好的隔熱作用；燒蝕層的石墨布、喉襯的碳碳復(fù)合材料和殼體的鋼導熱率較大，因而在燒蝕層和殼體內(nèi)部溫度變化較小，交界面之間的溫度變化也較小。

各界面溫度的最高點出現(xiàn)在噴管喉部的上游附近，沿噴管軸線方向，溫度先增大后減小。噴管內(nèi)流體溫度沿軸線方向逐漸減小，由于喉襯的導熱率較大、傳熱快，大量熱量通過交界面迅速地傳遞給燒蝕層，因而在上游出現(xiàn)最高溫。下游的流體溫度相對較低，溫度呈現(xiàn)下降趨勢。噴管尾部溫度下降的梯度較大，這是由于尾部端面的對流換熱及輻射換熱導致了大量熱量散失到外部自然來流。

圖3 各交界面的溫度分布Fig．3 Temperature distribution at each interface

圖4 文獻［3］中流固交接面的溫度分布Fig．4 Temperature distribution at the fluid－solid interface in the literature［3］

2．2 不同顆粒直徑對流場的影響

計算不同顆粒直徑的兩相流動時，選擇的顆粒直徑D分別為1!m、10!m、20!m和50!m，負載比取0．3。

圖5給出了純氣相及不同直徑顆粒下兩相流場中噴管溫度沿軸向的分布。圖6為軸線氣相的馬赫數(shù)分布。對比文獻［10］發(fā)現(xiàn)，本文的計算結(jié)果在軸向位置0．05 m之前的曲線分布趨勢同文獻計算結(jié)果保持一致。這是由于兩相流動時顆粒相的慣性以及顆粒相與氣相之間的黏性，軸線附近的流場受到顆粒相的傳熱作用，使得兩相流動時軸線溫度較純氣相情況下的軸線溫度更高。但在軸向位置0．05 m之后，曲線的分布有較大差異。這是由于本文的計算模型和文獻［10］所采用的噴管模型不相同，本文所選計算模型的擴張段更長，在擴張段尾部，直徑較大顆粒的隨流性較差，形成的無顆粒區(qū)域較大，大量粒子在靠近軸線附近運動，因而軸線附近的溫度就較高。這也說明，在噴管喉部下游的一小段區(qū)域內(nèi)，軸線馬赫數(shù)和溫度分布具有一定的規(guī)律，但是在下游一定區(qū)域之后，其趨勢并不呈現(xiàn)出明顯的規(guī)律。由此可知，不同顆粒直徑下兩相流場沿軸線的溫度分布與模型擴張段有重要關(guān)系。

圖5 不同顆粒直徑時軸線溫度分布Fig．5 Temperature distribution on the axis with different particle diameters

圖6 不同顆粒直徑時軸線氣相馬赫數(shù)分布Fig．6 Mach number distribution of axial gas phase with different particle diameters

圖7給出了純氣相及不同直徑顆粒下流固交界面、絕熱層和燒蝕層交界面、殼體與絕熱層交界面以及殼體外壁的溫度分布。由圖7可知，兩相流場對噴管結(jié)構(gòu)內(nèi)部溫度場的影響很大。與純氣相情況相比，兩相流動下，噴管各交界面在擴張段的溫度明顯偏高，而在噴管收縮段溫度無明顯差異，并且顆粒直徑越大，擴張段的溫度越低。

圖7 不同顆粒直徑時噴管各交界面溫度分布Fig．7 Temperature distribution at each interface of the nozzle with different particle diameters

2．3 不同負載比對流場的影響

針對不同負載比，計算兩相流動時的流場參數(shù)，選取的負載比分別為0．01、0．05、0．1、0．2和0．3，顆粒直徑為1!m。

圖8和圖9分別為不同負載比時沿軸線的溫度和馬赫數(shù)分布。

圖8 不同負載比時沿軸線的溫度分布Fig．8 Temperature distribution of the axis at different load ratios

圖9 不同負載比時沿軸線氣相的馬赫數(shù)分布Fig．9 Mach number distribution of axial gas phase at different load ratios

圖10所示分別為沿流固交界面、絕熱層和燒蝕層交界面、殼體和絕熱層交界面以及殼體外壁的溫度分布。由圖10可以看出，隨著負載比的減小，各交界面的溫度在噴管喉部及擴張段呈現(xiàn)下降的趨勢，而在收縮段溫度大小和變化趨勢基本一致，說明負載比對噴管喉部及擴張段結(jié)構(gòu)的溫度有較大的影響。

圖10 不同負載比時噴管各交界面的溫度分布Fig．10 Temperature distribution at each interface of the nozzle at different load ratios

3 熱固耦合計算結(jié)果與分析

通過對流場參數(shù)的計算獲得流固耦合交界面上的溫度分布后，可以進一步計算固體域的熱變形。為了對比分析兩相流動對噴管熱變形的影響，本文主要計算純氣相、不同顆粒直徑以及不同負載比時兩相流動的噴管熱固耦合，工況均與第2章所述工況相同。表3給出了噴管各結(jié)構(gòu)材料的力學參數(shù)。

噴管的網(wǎng)格劃分采用多區(qū)域法，并且在前部和尾部端面以及螺栓和殼體連接處加密?？偩W(wǎng)格數(shù)量為27．3萬。圖11給出了噴管的整體網(wǎng)格。

設(shè)置噴管外壁面及尾部端面為對流換熱和輻射換熱，對流換熱系數(shù)和發(fā)射率以及外部溫度的設(shè)置均與流場計算時采用的設(shè)置相同。設(shè)定螺栓的3個面為固定面，相鄰結(jié)構(gòu)之間采用綁定接觸，并且采用對稱接觸。

圖12所示為純氣相時噴管的溫度分布，可以看出，溫度的最高點出現(xiàn)在噴管喉部上游附近處，并且沿軸線方向先增大后減小，與流場計算結(jié)果的分布規(guī)律以及數(shù)值大小基本一致。

圖11 噴管整體網(wǎng)格Fig．11 The grids of the nozzle

圖12 噴管溫度分布Fig．12 Temperature distribution of nozzle

3．1 純氣相熱變形分析

圖13所示為純氣相情況下噴管的變形量分布。噴管的最大變形出現(xiàn)在噴管尾部，變形量為7．836 6 mm，最小值在螺栓處。由于螺栓的固定作用，熱膨脹產(chǎn)生的變形只能向自由端傳播，變形量沿軸線不斷增大，在尾部達到最大；同理，軸向變形的變形量分布與總變形量分布相同。徑向變形的最大值也出現(xiàn)在噴管的尾部端面，變形量為2．056 6 mm；最小值出現(xiàn)在燒蝕層前部，變形量為0．352 7 mm。由于殼體前端面固定，殼體在該處的變形量為零，而燒蝕層在殼體內(nèi)側(cè)，其徑向變形指向內(nèi)側(cè)，即指向圓心，因而變形量在數(shù)值上為負值。噴管外部的徑向變形量的分布與噴管內(nèi)壁溫度的變化趨勢一致，先增大后減小再增大，在噴管喉部上游附近出現(xiàn)高溫。

圖13 純氣相情況下噴管的變形量分布Fig．13 Displacement distribution of nozzle in pure gas phase

3．2 兩相流動對噴管變形影響分析

圖14為顆粒直徑D＝1μm時噴管的各變形量分布。由圖14可以看出，噴管變形量的分布趨勢同純氣相情況下的分布趨勢一致，但兩相流動下的噴管的總變形量、軸向變形量以及徑向變形量比純氣相情況下的量值要大，這是因為兩相流動下的噴管結(jié)構(gòu)內(nèi)部溫度相對較高。

圖14 顆粒直徑1μm情況下噴管的變形量分布Fig．14 Displacement distribution of nozzle with 1μm partical diameter

圖15為不同顆粒直徑時噴管各變形量的最大值。比較發(fā)現(xiàn)，兩相流動時的總變形量、軸向變形量和徑向變形量的最大值比純氣相情況下的要大，并且隨著顆粒直徑的增大，各變形量的最大值逐漸減小。表4為不同顆粒直徑時各變形量最大值的變化率，相比于純氣相時的各變形量最大值，D＝1μm時的總變形量最大值增大了12．9%，D＝50μm時總變形量最大值增大了2%。隨著顆粒直徑的增大，各變形量最大值的變化率也不斷減小。

圖15 不同顆粒直徑時各變形量最大值Fig．15 Maximum deformation at different particle diameters

表4 不同顆粒直徑時各變形量最大值變化率Tab．4 The change rate of maximum deformation at different particle diameters

圖16為不同負載比時噴管各變形量的最大值。兩相流動情況下，噴管的各變形量最大值均大于純氣相情況下的值，并且隨著負載比的增大，各變形量的最大值逐漸增大。表5為不同負載比時各變形量最大值的變化率，相比于純氣相時的各變形量最大值，η＝0．01時總變形量最大值增加了0．4%，η＝0．3時總變形量最大值增加了12．9%。隨著負載比的增大，變化率不斷增大。

圖16 不同負載比時的變形量Fig．16 Maximum deformation at different load ratios

表5 不同負載比時各變形量最大值變化率Tab．5 The rate of change of maximum deformation at different load ratios

4 結(jié)束語

本文計算了純氣相條件下噴管的流場并利用顆粒軌道模型計算了兩相流場；分析了噴管結(jié)構(gòu)內(nèi)部溫度場的變化以及兩相流動對噴管內(nèi)部溫度的影響；通過熱固耦合計算，分析了變形量的分布以及兩相流動對噴管熱變形的影響，研究表明：

1）在考慮噴管結(jié)構(gòu)傳熱性能時，流固交界面及噴管內(nèi)部的溫度分布與絕熱壁面的溫度分布有較大差異；流固交界面及噴管的各交界面沿軸向溫度先增大后減小，在喉部上游附近達到最大值，在噴管尾部端面處溫度會發(fā)生驟降；噴管的最低溫度出現(xiàn)在絕熱層尾部端面徑向位置最大處，在相鄰結(jié)構(gòu)交界處溫度會急劇變化；絕熱層內(nèi)的溫度梯度很大，降溫幅度很大，而燒蝕層和殼體內(nèi)部的溫度梯度很小，降溫幅度很??；

2）不同顆粒直徑時的兩相流動對噴管軸線處的溫度影響僅在喉部下游的一小段區(qū)域內(nèi)呈現(xiàn)規(guī)律性的變化，在此段區(qū)域之外，顆粒直徑的影響不遵循相應(yīng)的規(guī)律；

3）在溫度場的作用下，噴管總變形量和軸向變形量沿軸線方向不斷增大，在尾部變形量最大，徑向變形量沿軸向先增大后減小再增大，在噴管尾部達到最大，最小值出現(xiàn)在燒蝕層前部；

4）兩相流動對噴管的熱變形有著較大影響，噴管的總變形量、軸向變形量及徑向變形量均增大。