固體火箭發(fā)動(dòng)機(jī)噴管兩相流動(dòng)下的熱固耦合研究

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

（上海機(jī)電工程研究所，上海 201109）

0 引 言

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

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

1 數(shù)值計(jì)算方法

計(jì)算時(shí)進(jìn)行以下簡(jiǎn)化：不考慮化學(xué)反應(yīng)；不考慮顆粒相的燃燒、蒸發(fā)、破碎過程，且各個(gè)工況均為穩(wěn)態(tài)計(jì)算；不考慮內(nèi)壁面的燒蝕；不考慮顆粒接觸傳熱；不考慮各層材料之間的接觸熱阻；燃?xì)鉃槔硐霘怏w。

1．1 計(jì)算模型

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

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

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

1．2 計(jì)算方法

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

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

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

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

顆粒的位置可通過式（3）計(jì)算。

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

1．3 邊界條件

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

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

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

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

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

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

2 流場(chǎng)計(jì)算結(jié)果與分析

2．1 純氣相流場(chǎng)分析

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

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

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

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

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

2．2 不同顆粒直徑對(duì)流場(chǎng)的影響

計(jì)算不同顆粒直徑的兩相流動(dòng)時(shí)，選擇的顆粒直徑D分別為1!m、10!m、20!m和50!m，負(fù)載比取0．3。

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

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

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

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

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

2．3 不同負(fù)載比對(duì)流場(chǎng)的影響

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

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

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

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

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

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

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

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

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

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

圖12所示為純氣相時(shí)噴管的溫度分布，可以看出，溫度的最高點(diǎn)出現(xiàn)在噴管喉部上游附近處，并且沿軸線方向先增大后減小，與流場(chǎng)計(jì)算結(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)生的變形只能向自由端傳播，變形量沿軸線不斷增大，在尾部達(dá)到最大；同理，軸向變形的變形量分布與總變形量分布相同。徑向變形的最大值也出現(xiàn)在噴管的尾部端面，變形量為2．056 6 mm；最小值出現(xiàn)在燒蝕層前部，變形量為0．352 7 mm。由于殼體前端面固定，殼體在該處的變形量為零，而燒蝕層在殼體內(nèi)側(cè)，其徑向變形指向內(nèi)側(cè)，即指向圓心，因而變形量在數(shù)值上為負(fù)值。噴管外部的徑向變形量的分布與噴管內(nèi)壁溫度的變化趨勢(shì)一致，先增大后減小再增大，在噴管喉部上游附近出現(xiàn)高溫。

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

3．2 兩相流動(dòng)對(duì)噴管變形影響分析

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

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

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

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

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

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

圖16 不同負(fù)載比時(shí)的變形量Fig．16 Maximum deformation at different load ratios

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

4 結(jié)束語

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

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

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

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

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