白俊華, 胡春波, 李佳明
(西北工業(yè)大學(xué) 燃燒、熱結(jié)構(gòu)與內(nèi)流場重點(diǎn)實(shí)驗(yàn)室, 陜西 西安 710072)
隨著高能固體推進(jìn)劑的廣泛應(yīng)用,燃燒室燃燒產(chǎn)物的溫度逐漸提高,在發(fā)動機(jī)工作過程中,噴管喉襯需要暴露在最嚴(yán)酷的傳熱、熱應(yīng)力和高溫環(huán)境中,噴管喉部的熱流密度高達(dá)12 MW/m2以上。急劇的熱交換造成了噴管喉部嚴(yán)重的燒蝕現(xiàn)象,所以發(fā)動機(jī)噴管喉襯熱結(jié)構(gòu)以及換熱規(guī)律的研究對噴管結(jié)構(gòu)設(shè)計優(yōu)化具有重大的意義。
研究噴管的換熱情況以及溫度分布是噴管熱分析與設(shè)計的主要途徑。而發(fā)動機(jī)噴管的傳熱是復(fù)雜的傳導(dǎo)、對流、輻射的耦合換熱過程,實(shí)際噴管的溫度分布不僅與其工況、粘性邊界層狀態(tài)等有關(guān),而且與噴管本身的結(jié)構(gòu)、材料等密切相關(guān),所以必須作為耦合傳熱問題求解。國內(nèi)外學(xué)者針對噴管喉部的熱交換規(guī)律及影響因素開展了一定的理論研究和實(shí)驗(yàn)分析。熊永亮等人[1]建立了軸對稱的有限元計算模型,通過給定熱邊界條件對喉襯組件和擴(kuò)張段結(jié)構(gòu)的瞬態(tài)溫度分布情況進(jìn)行了研究。Thakre等人[2]建立了碳/碳以及石墨喉襯熱化學(xué)燒蝕的理論模型,該模型考慮了推進(jìn)劑的化學(xué)反應(yīng)、燃?xì)鉄嵛镄宰兓庀嗟幕瘜W(xué)動力學(xué)和噴管表面的異相反應(yīng)等問題。陳博等人[3]利用酒精/氧氣燃?xì)獍l(fā)生器模擬火箭發(fā)動機(jī)富氧燃?xì)猸h(huán)境,從而分析了燃?xì)鈪?shù)對燒蝕性能的影響。李佳明等人[4]以鍛壓鎢為材料設(shè)計噴管,通過在噴管喉襯徑向的一定位置填埋熱電偶來測量噴管結(jié)構(gòu)溫度,從而剝離了熱化學(xué)燒蝕對噴管喉部的影響,獲得噴管內(nèi)部結(jié)構(gòu)的溫度場分布。雖然上述研究工作取得了一定成果,但是仍然缺乏精準(zhǔn)的計算模型來獲取噴管喉部流固交界面的相關(guān)參數(shù),從而使噴管的研究工作受阻。
本文建立了噴管喉襯流固耦合換熱模型,通過對參考文獻(xiàn)[4]實(shí)驗(yàn)工況的計算,驗(yàn)證了計算模型的準(zhǔn)確性;應(yīng)用該數(shù)值仿真模型對堵多因素對噴管喉襯熱交換規(guī)律的影響進(jìn)行了研究,為優(yōu)化噴管設(shè)計提供理論參考。
本文選用參考文獻(xiàn)[4]的實(shí)驗(yàn)構(gòu)型尺寸,考慮噴管構(gòu)型的軸對稱性,為了減少計算量,建立二維軸對稱數(shù)值模型。計算區(qū)域型如圖1所示,計算區(qū)域中,參照文獻(xiàn)[4]取石墨組件軸向長25 mm,徑向厚度10 mm;噴管套層采用45#鋼制作,徑向厚度為20 mm,噴管喉襯的結(jié)構(gòu)尺寸如圖2a)所示。在噴管壁面處劃分邊界層網(wǎng)格,底層網(wǎng)格厚度為0.01 mm,模型總網(wǎng)格數(shù)約為57萬,單個網(wǎng)格體積約為5.2×10-9m3。采用SIMPLE算法,二階迎風(fēng)格式。
圖1 計算區(qū)域
圖2 噴管結(jié)構(gòu)及測溫點(diǎn)位置示意圖[4]
本文采用的是粘性可壓非穩(wěn)態(tài)二維N-S控制方程。噴管喉部換熱模型的能量守恒方程包含流體換熱、流固耦合換熱以及固體內(nèi)部的熱傳導(dǎo)。
氣相能量守恒方程如下:
(1)
式中:keff=k+ki為有效導(dǎo)熱率,ki為湍流引起的導(dǎo)熱率,由采用的湍流模型以及壁面函數(shù)確定。等號右邊的前3項(xiàng)分別表示由于熱傳導(dǎo)、組分?jǐn)U散、黏性耗散引起的能量轉(zhuǎn)移,源項(xiàng)Sh是由化學(xué)反應(yīng)、輻射等其他因素引起的熱源。
在固體區(qū)域僅考慮導(dǎo)熱,如(2)式所示,建立熱傳導(dǎo)方程:
(2)
式中:速度v表示固體區(qū)域由于旋轉(zhuǎn)或平移等運(yùn)動的速度,·(kT)表示熱傳導(dǎo)引起的熱流,源項(xiàng)Sh表示固體區(qū)域的內(nèi)熱源。
為了獲取更精準(zhǔn)的噴管喉部換熱數(shù)值仿真模型,針對時均化過程中產(chǎn)生的脈動量乘積建立Reynolds應(yīng)力微分方程式進(jìn)行求解,具體如下:
(3)
式中:DT,ij表示湍流擴(kuò)散項(xiàng);DL,ij是分子擴(kuò)散項(xiàng);Pij為湍流切應(yīng)力在主流速度方向做的功;φij表示壓力應(yīng)變項(xiàng)。
在噴管喉襯內(nèi)壁面處,對低Re數(shù)的湍流耗散率方程ω加入粗糙度的修正。
(4)
(5)
(6)
式中:Ks表示壁面當(dāng)量粗糙度。
在固體壁面附近,由于分子黏性的作用,湍流脈動受到阻尼,因此對壁面附近分子黏性起作用的區(qū)域要做特別的處理。在劃分網(wǎng)格時,把第一個內(nèi)節(jié)點(diǎn)P布置到對數(shù)分布律成立的范圍內(nèi),即配置到旺盛湍流區(qū)域。第一個內(nèi)節(jié)點(diǎn)與壁面之間區(qū)域的當(dāng)量導(dǎo)熱系數(shù)λt按如下公式確定:
(7)
式中:qw由對數(shù)分布定律所規(guī)定,Tw為壁面上的溫度,根據(jù)此式可以導(dǎo)得第一個內(nèi)節(jié)點(diǎn)上當(dāng)量導(dǎo)熱系數(shù)λt的計算式。
溫度對數(shù)分布律的表達(dá)式為:
(κ′=0.465,Ec=4.75)
(8)
由(7)式、(8)式可導(dǎo)出壁面上當(dāng)量導(dǎo)熱系數(shù)λt的表達(dá)式。
燃?xì)廨椛淠P筒捎没覛怏w加權(quán)平均模型。該模型的基本假設(shè)是對于一定厚度的氣體吸收量,其發(fā)射率為:
(9)
式中:aε,i為第i組組分的發(fā)射率加權(quán)系數(shù);κi為第i組組分的吸收系數(shù);p為所有吸收性氣體的分壓總和;s為輻射的行程長度。
計算區(qū)域如圖1所示,噴管入口根據(jù)實(shí)驗(yàn)過程中采集到的燃燒室平均壓強(qiáng)給定壓力入口邊界;噴管出口建立外推邊界條件,忽略大氣環(huán)境對噴管內(nèi)流動的影響;石墨及45#鋼左側(cè)邊界設(shè)置為絕熱;45#鋼其余邊界均設(shè)定與大氣環(huán)境發(fā)生自然對流換熱;計算區(qū)域其他不同組件的交界面均設(shè)置為耦合邊界條件。
在數(shù)值仿真計算區(qū)域相應(yīng)給定如參考文獻(xiàn)[4]中3個測點(diǎn)的位置,其坐標(biāo)位置如圖2b)所示,提取發(fā)動機(jī)穩(wěn)定工作時間段5~9 s時測點(diǎn)(A、B、C)的溫度隨時間的變化,可以得到數(shù)值仿真結(jié)果與實(shí)驗(yàn)結(jié)果的對比,如表1所示。
表1 計算與實(shí)驗(yàn)結(jié)果對比
由表1計算與實(shí)驗(yàn)結(jié)果對比可知,在發(fā)動機(jī)穩(wěn)定工作段(5~9 s)時,數(shù)值仿真結(jié)果與實(shí)驗(yàn)結(jié)果吻合良好,最大誤差不超過10%。從而驗(yàn)證了流固耦合換熱數(shù)值模型的準(zhǔn)確性。
針對噴管壁面粗糙度、燃?xì)庵休椛錃怏w組分、燃燒室壓強(qiáng)及推進(jìn)劑燃溫等因素對噴管喉襯的換熱影響開展數(shù)值仿真研究,計算工況如表2所示。
表2 噴管喉部換熱數(shù)值仿真工況
由于噴管喉襯惡劣的工作環(huán)境,在噴管喉襯內(nèi)壁面發(fā)生復(fù)雜且嚴(yán)重的燒蝕現(xiàn)象,導(dǎo)致噴管喉襯內(nèi)壁面是極其粗糙的。本文以當(dāng)量粗糙度的概念表征壁面粗糙程度。對噴管喉襯壁面粗糙度為0.1 mm、0.3 mm、0.5 mm時與光滑壁面即0 mm時的工況(1~4)進(jìn)行仿真。不同壁面粗糙度下噴管喉襯內(nèi)壁面溫度以及熱交換的熱流密度如圖3、圖4所示。
圖3 粗糙度對壁面溫度的影響
圖4 粗糙度對壁面熱流密度的影響
仿真結(jié)果表明,隨著噴管喉襯內(nèi)壁面粗糙度的增加,噴管喉襯內(nèi)壁面溫度上升,熱交換加劇,熱交換的熱流密度增加。在發(fā)動機(jī)噴管喉襯擴(kuò)張段,由于燃?xì)饬魉僭龃?燃?xì)鈩幽苻D(zhuǎn)化為勢能的部分增大,所以燃?xì)獗诿娲植诙葘α鲃拥挠绊懜鼮閺?qiáng)烈,熱流密度分布出現(xiàn)一個爬升峰值。同時,由于噴管喉襯擴(kuò)張段主流溫度急劇下降,因此擴(kuò)張段壁面溫度主要呈現(xiàn)下降趨勢,僅出現(xiàn)一個較小的峰值。粗糙度為0~0.1 mm時溫度與熱流密度相差很大,平均溫度相差280℃左右,熱流密度相差約2 MW/m2。但0.1~0.3 mm工況下的差別相對較小。
不同燃?xì)饨M分下(工況5、6、7),噴管喉襯內(nèi)壁面溫度以及輻射換熱的熱流密度如圖5、圖6所示。
圖5 燃?xì)饨M分對壁面溫度的影響
圖6 燃?xì)饨M分對壁面輻射熱流密度的影響
工況5和工況7對比結(jié)果表明,收斂段的溫度相差60℃左右,熱流密度約差0.5 MW/m2,但在直段和擴(kuò)張段,二者相差很小。在噴管喉襯收斂段,由于燃?xì)鉁囟容^高,輻射換熱較為強(qiáng)烈,因此當(dāng)燃?xì)饨M分中強(qiáng)輻射氣體的質(zhì)量分?jǐn)?shù)增加時,噴管喉襯收斂段的壁面溫度隨著強(qiáng)輻射氣體質(zhì)量分?jǐn)?shù)的增加而升高。而在噴管喉襯直段及其擴(kuò)張段,由于燃?xì)鉁囟缺容^低,輻射換熱相對較弱,因此燃?xì)庵袕?qiáng)輻射氣體質(zhì)量分?jǐn)?shù)的增加對噴管喉襯壁面溫度以及輻射換熱強(qiáng)度影響較小。
不同燃燒室壓強(qiáng)下即工況1、8、9噴管喉襯內(nèi)壁面溫度以及熱交換的熱流密度如圖7、圖8所示。仿真結(jié)果表明,燃燒室壓強(qiáng)的增加導(dǎo)致噴管喉襯內(nèi)壁面換熱加劇、壁面溫度升高。噴管喉襯內(nèi)壁面溫度和總熱流密度在噴管喉部直段的前端達(dá)到最大。在發(fā)動機(jī)噴管喉襯直段,溫度變化較為平緩,在噴管喉襯的擴(kuò)張段,沿噴管喉部的距離逐漸增加,總熱流密度逐漸降低。
圖7 燃燒室壓強(qiáng)對壁面溫度的影響
圖8 燃燒室壓強(qiáng)對壁面熱流密度的影響
設(shè)定燃燒室壓強(qiáng)為7.5 MPa,針對燃?xì)鉁囟确謩e為2 000 K、2 500 K和3 000 K的工況(1、5、10)開展數(shù)值仿真,提取發(fā)動機(jī)工作4 s時噴管喉襯內(nèi)壁面沿軸向的溫度分布及總熱流密度如圖9、圖10所示。仿真結(jié)果表明,燃?xì)鉁囟葘姽芎硪r的熱交換以及溫度分布影響很大,燃溫越高則壁面溫度越高、熱流密度越大。初始燃?xì)鉁囟认嗖?00 K時,內(nèi)壁面溫度約差200 K,總熱流密度增加1 MW/m2左右。當(dāng)燃?xì)鉁囟葹? 000 K時,噴管喉部熱交換總的熱流密度高達(dá)8 MW/m2以上。
圖9 推進(jìn)劑燃溫對壁面溫度的影響
圖10 推進(jìn)劑燃溫對壁面熱流密度的影響
本文基于二維軸對稱計算構(gòu)型建立了噴管喉襯流固耦合的換熱模型,通過與噴管喉襯熱結(jié)構(gòu)測溫結(jié)果[4]做對比分析,驗(yàn)證了模型的準(zhǔn)確性。并且針對噴管壁面粗糙度、燃?xì)饨M分等因素對噴管喉部的換熱規(guī)律進(jìn)行了分析研究。仿真結(jié)果表明:
1) 噴管喉襯內(nèi)壁面溫度和總熱流密度在噴管喉部直段的前端達(dá)到最大;
2) 隨著粗糙度增大,噴管喉襯內(nèi)壁面的溫度和總熱流密度增大;在發(fā)動機(jī)噴管喉襯擴(kuò)張段,壁面粗糙度對熱流密度的影響更為強(qiáng)烈;
3) 隨著燃?xì)饨M分中強(qiáng)輻射氣體質(zhì)量分?jǐn)?shù)的增加,在噴管喉襯收斂段,輻射熱流密度逐漸增加,在噴管喉襯直段以及擴(kuò)張段,輻射熱流密度變化不大;
4) 隨著燃燒室壓強(qiáng)以及推進(jìn)劑燃溫的升高,噴管喉襯內(nèi)壁面溫度和總熱流密度增大。
參考文獻(xiàn):
[1] 熊永亮,郜冶. 噴管溫度與應(yīng)力場的數(shù)值研究[J]. 哈爾濱工程大學(xué)學(xué)報, 2007, 28(8): 852-855
Xiong Yongliang, Gao Ye. Numerical Study of Temperature and Stress Field of Nozzles[J]. Journal of Harbin Engineering University, 2007, 28(8): 852-855 (in Chinese)
[2] Thakre P, Yang V. Chemical, Erosion of Carbon-Carbon/Graphite Nozzles in Solid-Propellant Rocket Motors[J]. Journal of Propulsion and Power, 2008, 24(4): 40-50
[3] 陳博,張立同. 燃?xì)獍l(fā)生器條件下穿刺復(fù)合材料噴管的燒蝕性能研究[J]. 無機(jī)材料學(xué)報,2008,23(6):1159-1164
Chen Bo, Zhang Litong. Ablation Characteristic of the Pierced C/C Composite Nozzle in a Combustion Chamber Generator[J]. Journal of Inorganic Material, 2008,23(6):1159-1164 (in Chinese)
[4] 李佳明,胡春波. 固體發(fā)動機(jī)噴管喉襯溫度場測量與分析[J]. 實(shí)驗(yàn)流體力學(xué),2012,26(5): 57-60
Li Jiaming, Hu Chunbo. The Measurement and Analysis of Temperature Field in a Solid Motor Nozzle Throat Insert[J]. Journal of Experiments in Fluid Mechanics, 2012, 26(5): 57-60 (in Chinese)
[5] Evans B, Kuo K K, Ferrara P J. Characterization of Nozzle Erosion Phenomena in a Solid-Propellant Rocket Motor Simulator[C]∥44th AIAA/ASME/SAE/ASEE Joint Propulsion Conference. Hartford, CT, 2008
[6] Yu D, Toru S. Evaluation of Ablation and Longitudinal Vortices in Solid Rocket Motor by Computational Fluid Dynamics[C]∥42th AIAA/ASME/SAE/ASEE Joint Propulsion Conference. California, 2006