復合噴管結構溫度場流固耦合仿真分析①

傅學金，鄧 恒，寇軍強，宋 軍，陳 剛

(1.第二炮兵工程大學，西安 710025;2.西安現(xiàn)代控制技術研究所，西安 710065;3.第二炮兵裝備部駐航天科技集團型號辦事處，北京 100076)

復合噴管結構溫度場流固耦合仿真分析①

傅學金1，2，鄧 恒2，寇軍強2，宋 軍2，陳 剛3

(1.第二炮兵工程大學，西安 710025;2.西安現(xiàn)代控制技術研究所，西安 710065;3.第二炮兵裝備部駐航天科技集團型號辦事處，北京 100076)

為了評估復合噴管熱防護性能以及獲取噴管燒蝕和結構應力分析的工況條件，運用Fluent流體動力學軟件，對復合噴管的結構溫度場進行了數(shù)值仿真。分析中，采用了兩方程RNG k-ω湍流模型和增強型壁面函數(shù)，利用流固耦合的計算方法，獲得了噴管結構瞬態(tài)溫度場的計算結果，重點分析了結構溫度場最終分布狀態(tài)和初期傳播特點，以及喉襯溫度隨時間的變化規(guī)律，估算了喉襯的燒蝕。分析結果表明，噴管結構熱防護性能滿足要求，溫度最高區(qū)域位于噴管收斂段中后部，喉襯線燒蝕量約為2.1 mm，為噴管結構進一步優(yōu)化設計提供了重要參考依據(jù)。

流固耦合;溫度場;復合噴管;數(shù)值仿真

0 引言

在高壓力、高燃溫、長時間工作的中大型發(fā)動機設計中，多選用復合材料噴管的方案。復合噴管的結構設計是噴管設計中十分重要的內(nèi)容，包括了耐燒蝕層、隔熱層和喉襯的結構設計和各部件之間的接縫匹配設計。耐燒蝕層、隔熱層和喉襯的結構設計主要是以對噴管的燒蝕和傳熱的評估作為依據(jù)，而各部件之間的接縫匹配設計則是以噴管的熱應力分析作為依據(jù)的。無論是噴管的燒蝕、傳熱和熱應力的分析，都要以噴管溫度場的分析為前提條件。因此，在噴管結構設計中，噴管結構溫度場分析是十分重要的工作。

目前，對噴管結構溫度場的計算有3種計算方法:(1)對噴管內(nèi)壁面施加溫度邊界條件，計算噴管結構的溫度場［1］;(2)考慮燃氣與壁面的對流換熱作用，利用巴茲公式計算燃氣與壁面之間的對流換熱系數(shù)，以此作為邊界條件計算結構的溫度場［2－8］;(3)將噴管內(nèi)流場與結構的交界面作為流體域和固體域的熱耦合邊界，進行流場和傳熱的耦合計算［9－11］。在第1種方式中，對噴管內(nèi)壁面施加的是溫度邊界條件，用熱傳導的傳熱方式替代了噴管內(nèi)壁面與燃氣之間的對流換熱方式，這改變了邊界條件的性質(zhì)。在第2種方式中，雖考慮到了噴管內(nèi)壁面與燃氣的對流換熱邊界條件，但巴茲公式所需的壁面溫度、壓力和馬赫數(shù)都是從絕熱壁面邊界條件的噴管流場計算結果中取得的，沒有考慮結構對溫度場結果的耦合影響作用。第3種計算方式就將流場和溫度場耦合起來進行計算，這符合噴管的實際狀態(tài)，計算結果更為準確可信。文獻［10］就明確指出，噴管內(nèi)流場和結構溫度場的解耦計算是不合適的;文獻［11］則利用試驗測試的噴管內(nèi)壁面的溫度數(shù)據(jù)驗證了耦合計算的正確性，與解耦計算結果的對比分析表明，耦合計算比解耦計算結果更為準確。綜上分析，本文將采用流場與結構耦合傳熱的方式計算某型噴管結構的溫度場，分析初始階段和工作結束時溫度場，以及喉襯溫度隨時間的變化規(guī)律，為噴管結構優(yōu)化設計提供重要參考依據(jù)。

1 計算方法

采用Fluent軟件，對噴管的流場和流場對固體的傳熱進行整體分析。Fluent軟件運用了有限體積法離散流體動力學微分控制方程，將其變換為一組代數(shù)方程來求解。流體動力學控制方程主要包括了質(zhì)量守恒方程、動量守恒方程和能量守恒方程，可寫為如下的通用形式:

式中 φ為通用變量，可代表速度u、v、w和溫度T等變量;Γ為廣義擴散系數(shù);s為廣義源項，它們在不同的控制方程中具有不同的表達式。

對于涉及到傳熱求解的問題，必須使用能量守恒方程來對溫度場進行求解，此處特將能量方程展開如下:

式中 T為溫度;k為傳熱系數(shù);cp為比熱容;sT為粘性耗散項。

從式(2)可看出，在流體域內(nèi)，溫度的變化與密度、速度和溫度梯度都有關系，在固體域內(nèi)，介質(zhì)的速度為零，密度為常數(shù)，因而溫度的變化僅與溫度的梯度有關系，在流體域和固體域的耦合界面上，溫度在要滿足對流體域內(nèi)溫度求解的同時，還要滿足對固體域的求解，即符合連續(xù)性的要求。從以上分析可看出，要獲得噴管結構準確的溫度場，必須將流場域的流動分析與結構的固體傳熱分析耦合起來。因此，設置流體域與固體域直接的邊界為耦合邊界。

在工業(yè)化和城市化快速推進過程中，歐美國家也經(jīng)歷了鄉(xiāng)村發(fā)展滯后于城市、農(nóng)村居民收入增長滯后于城市居民收入增長的不協(xié)調(diào)階段。中共中央黨校(國家行政學院)《行政管理改革》編輯部編輯張世貴梳理了發(fā)達國家破解城鄉(xiāng)發(fā)展差距問題的經(jīng)驗。

由于噴管的結構溫度場取決于噴管內(nèi)的流動，所以對噴管內(nèi)流場的準確分析顯得十分必要。湍流模型對噴管內(nèi)流場的計算具有很大的影響，在本分析中，選擇RNG k-ε湍流模型。該模型應用的范圍比標準k-ε模型更廣泛，對低雷諾數(shù)和高雷諾數(shù)的流動都能進行很好的分析，能更好地反映出高切應變率和流線曲率對流動影響的效應。固體火箭發(fā)動機噴管內(nèi)的燃氣流動經(jīng)歷了收斂段的亞音速、喉部的音速和擴張段的超音速，雷諾數(shù)可能從低雷諾數(shù)變化到高雷諾數(shù)，并且先收斂、后擴張的構型，使得靠近噴管壁面的流線存在較大的曲率。因此，RNG k-ε湍流模型適合對噴管內(nèi)流場的計算。

對噴管壁面函數(shù)的處理直接關系到流固界面的傳熱，從而直接影響到噴管結構的溫度場，因而選擇合適的壁面函數(shù)也是噴管結構溫度場分析的一個重要方面。本分析中，將選擇增強的壁面處理(Enhanced Wall Treatment)，該壁面函數(shù)結合了壁面的二層模型。當邊界層網(wǎng)格細密到可滿足直接求解粘性邊界層時為二層模型，當臨近邊界層網(wǎng)格較稀疏時，則采用增強的壁面函數(shù)求解邊界層。

2 復合噴管結構溫度場分析

2.1 物理問題

為了給噴管熱噴管的耐燒蝕層和背壁厚度設計提供依據(jù)，以及將來為喉襯熱的應力分析和接縫設計提供輸入條件，分析噴管結構的溫度場分布。某發(fā)動機復合噴管結構如圖1所示，收斂段和擴張段的耐燒蝕層為碳/酚醛材料，喉襯為炭/炭材料，噴管背壁為高硅氧/酚醛材料，最外層為鋼殼體。發(fā)動機燃燒室平均壓力為8 MPa，工作時間為35 s。

圖1 噴管結構簡圖Fig.1 Nozzle structural schematic

對于該噴管的結構溫度場分析，可忽略一些次要影響因素，對模型進行必要的簡化，作以下假設和簡化:

(1)燃氣簡化為理想氣體

(2)不考慮噴管各部件之間縫隙中的粘接劑和填充物，各部件無縫隙接觸;

(3)忽略壁面的炭化、燒蝕和輻射效應，只考慮燃氣與噴管內(nèi)壁面的對流換熱效應;

(4)忽略各部件之間的接觸熱阻;

(5)真實的高硅氧/酚醛背壁是分段的，在該模型中合并為一段。

2.2 計算模型

針對噴管的軸對稱特點，建立噴管的軸對稱計算模型，將噴管的內(nèi)流場區(qū)域設為流體域，將噴管的各部件設為固體域，如圖2所示。對各區(qū)域分別采用四邊形網(wǎng)格進行劃分，如圖3所示。噴管內(nèi)流場區(qū)域與內(nèi)壁面之間建立耦合邊界條件，噴管外壁面邊界設為與環(huán)境的對流換熱邊界，對流換熱系數(shù)設為5 J/m2。噴管的入口設為壓力入口邊界，壓力為燃燒室平均壓力8 MPa，燃氣溫度3 200 K，出口為壓力邊界，為一個大氣壓，溫度為300 K。材料參數(shù)如表1所示。

圖2 噴管軸對稱計算模型Fig.2 Nozzle axisymmetric model

圖3 噴管網(wǎng)格模型Fig.3 Nozzle computational mesh

表1 材料參數(shù)Table 1 Parameters of nozzle materials

在復合噴管中，流場的發(fā)展速率遠大于結構的熱傳導速率，因而在短時的點火增壓過程，可忽略結構的熱傳導效應，噴管結構的熱傳導計算就可從發(fā)動機進入平穩(wěn)工作階段時開始。如此，在噴管結構的熱傳導計算過程中，噴管流場都是穩(wěn)定狀態(tài)。因此，在本文分析中，首先對噴管進行穩(wěn)態(tài)內(nèi)流場分析;然后，在噴管穩(wěn)態(tài)內(nèi)流場的結果基礎上，對噴管的熱傳導進行瞬態(tài)分析。具體的分析過程如圖4所示。

圖4 噴管結構分析過程Fig.4 Nozzle fluid structrue analysis process schematic

2.3 計算結果及分析

2.3.1 工作35 s時溫度場分布

噴管在35 s時溫度場分布云見圖5。從圖5中可看出，熱沿噴管內(nèi)壁面近似平行向外壁面?zhèn)鬟f，噴管內(nèi)壁面附近區(qū)域為高溫區(qū)，逐漸向外壁降低，溫度的傳遞僅通過了由炭/酚醛收斂段、喉襯和炭/酚醛擴張段組成的第一層結構，高硅氧/酚醛背壁還處于300～500 K的溫度范圍，鋼殼體外壁溫度為300 K。由此可見，該噴管的結構是滿足熱防護的要求的，可能還有較大的余量。

從圖5中還可明顯看出，最高溫度位于噴管內(nèi)壁在從收斂段中部到喉部的一段區(qū)間，在2 680～2 860 K之間，見圖6。炭/酚醛收斂段與喉襯的接縫就位于該處，因此要注意此處燒蝕和接縫設計。

圖5 35 s噴管溫度場分布Fig.5 Nozzle temperature distribution at 35 s

根據(jù)噴管結構分布特點，分別在收斂段、喉部和擴張段取截面(見圖7)，考察溫度各部位沿徑向截面的分布(見圖8～圖10)。由曲線圖8～圖10可看出，溫度沿截面從內(nèi)壁向外壁逐漸降低，在與高硅氧/酚醛背壁的交界面上曲線發(fā)生了轉折，表明界面兩側材料參數(shù)的不同，導致了溫度分布的一階導數(shù)不連續(xù)，圖中橫坐標作了歸一化處理。高硅氧/酚醛背壁的溫度都在300～500 K范圍內(nèi)，其外壁面溫度為300 K。顯然，鋼殼體的溫度還屬于常溫300 K，沒有發(fā)生溫升。收斂段和擴張段截面的溫度呈線性下降，而喉部的截面的溫度下降逐漸趨緩，為非線性。喉部3個截面的溫度分布趨勢相同。其中，截面1溫度最高，截面2次之，截面3最低。表明喉襯在入口部位是整個喉襯的高溫區(qū)域，燒蝕可能最為嚴重。

圖6 噴管最高溫度區(qū)Fig.6 The highest temperature zone of nozzle

圖7 提取溫度的截面位置Fig.7 The positions where temperature is extracted

圖8 35 s收斂段截面溫度分布Fig.8 Temperature change with radial direction distance at nozzle convergent cross section at 35 s

圖9 35 s喉襯截面上的溫度分布Fig.9 Temperature change with radial direction distance at nozzle throat cross section at 35 s

圖10 35 s擴張段截面溫度分布Fig.10 Temperature change with radial direction distance at nozzle divergent cross section at 35 s

2.3.2 工作初期溫度發(fā)展情況

考察了5.25 s內(nèi)喉襯以及與收斂段和擴張段接縫附近的溫度變化，見圖11。由圖11可看出，炭/炭材料喉襯的溫度傳導速度明顯高于兩側炭/酚醛材料的收斂段和擴張段，熱傳導開始時，接縫兩側溫差較大。隨時間推移，接縫處喉襯一側的溫度向高硅氧/酚醛的收斂段和擴張段一側傳導，接縫兩側的溫差逐漸縮小。

圖11 5.25 s內(nèi)溫度變化情況Fig.11 Temperature evolution within 5.25 s

接縫應力是接縫寬度設計的主要依據(jù)，其最大值往往發(fā)生在工作初期階段。接縫應力取決于接縫的接觸力和接觸面積。接觸力會隨溫度由內(nèi)壁向外傳導得越充分而變大，因為收斂段、喉襯和擴張段在溫度的作用下，將向軸向膨脹，相互之間發(fā)生擠壓，溫度越高，擠壓越嚴重，接觸力也就會越大。另一方面，溫度逐漸向外層傳導，具有一定寬度的接縫由開始不接觸到接觸，之后熱傳導越充分，接觸面積會越大，有降低接觸應力的趨勢。最大接觸應力發(fā)生在何時有多大與噴管的長度、材料的性質(zhì)、縫隙寬度、結構的厚度都有關系，在接縫設計時要考慮到這些因素。

2.3.3 喉襯溫度分布及變化

喉襯會因自身的熱應力產(chǎn)生破壞，因此有必要考察喉襯的溫度分布及變化情況。工作35 s結束時，喉襯的溫度分布云圖如圖12所示。整體上，溫度依據(jù)喉襯弧形內(nèi)型面由內(nèi)向外、由高到低的分布，喉襯入口段的溫度最高，達到了2 833 K，外壁中部區(qū)域溫度最低，在600～765 K之間，內(nèi)外溫差較大，到達了2 000 K左右。在這樣的溫度分布下，喉襯在兩個方向上可能存在發(fā)生斷裂的風險:一是軸向上，內(nèi)壁溫度高，外壁溫度低，喉襯產(chǎn)生由內(nèi)向外的彎曲變形，這使得喉襯有發(fā)生軸向斷裂的風險;二是周向上，外壁約束了內(nèi)壁周向膨脹，有發(fā)生周向斷裂的風險。是否發(fā)生破壞，還需要應力的分析。

觀察喉襯沿截面1(見圖7)由內(nèi)壁到外壁的溫度分布隨時間的變化，見圖13。

圖12 35 s喉襯溫度分布Fig.12 Throat temperatue distribution at 35 s

圖13 截面1溫度隨時間變化Fig.13 Temperatue change with time at the first cross section of throat in Fig.7

由圖13可知，隨著時間的推移，整個截面上的溫度在逐漸增高，同時溫度在逐步向外壁傳播，靠近內(nèi)壁的區(qū)域在開始階段升溫較快，然后趨緩，靠近外壁區(qū)域的升溫則較為勻速。由圖14可知，內(nèi)壁面溫度在5 s內(nèi)急劇上升，之后趨于平緩。

圖14 截面1內(nèi)壁面溫度隨時間的變化Fig.14 Inner wall temperature change with time at the first cross section of throat in Fig.7

3 喉襯燒蝕估算

由于喉襯壁面溫度較高，在湍流邊界層內(nèi)CO2和H2O與喉襯壁面的碳反應速率較快，大于CO2和H2O的向邊界層的擴散速率。因此，炭的質(zhì)量耗損計算采用擴散控制模型［12］，其炭的質(zhì)量耗損率表示式為

本計算中，CCO2e=0.049 3，CH2Oe=0.426 7，cv=3 208 J/(kg·K)，=0.027 3 kg/mol，hg從仿真計算的結果中提取出來，其隨時間的變化曲線如圖15所示。經(jīng)計算得到mw=3.766 kg/m2，折算為線燒蝕量r=2.1 mm(喉襯密度為 1 800 kg/m3)。

圖15 喉徑處對流傳熱系數(shù)時間曲線Fig.15 Time curve of throat heat transfer coefficient

4 結論

(1)工作35 s結束時，噴管溫度僅傳導過第一層結構，噴管的結構滿足熱防護要求。

(2)溫度最高區(qū)域位于內(nèi)壁收斂段中后部，在2 680～2 860 K之間，此處可能是燒蝕最為嚴重的地方。

(3)喉襯溫度由內(nèi)到外從高到低平行分布，最高溫度為2 833 K，位于入口段，內(nèi)外溫差達到2 000 K。

(4)喉襯線燒蝕量估算為2.1 mm。

［1］沈勇.碳/碳固體火箭發(fā)動機柔性噴管熱及熱-結構耦合作用研究［D］.華中科技大學，2004.

［2］王天波，薛談順，周長省，等.復合結構噴管溫度場及應力場數(shù)值模擬［J］.彈道學報，2012，24(2):88-105.

［3］王作良.噴管熱防護的有限元數(shù)值分析［D］.哈爾濱工程大學，2006.

［4］藤亞莉，黃立威.復合材料噴管熱應力分析［C］//第七屆中國CAE工程分析技術年會論文集，2011:352-356.

［5］熊永亮，郜冶.噴管溫度與應力場的數(shù)值研究［J］.哈爾濱工程大學學報，2007，28(8):852-855.

［6］胡江華，孟松鶴，常新龍，等.錐形套式連接C/C噴管擴張段溫度場與應力場分析［J］.固體火箭技術，2012，35(1):64-68.

［7］田四朋，唐國金，李道奎，等.固體火箭發(fā)動機噴管結構縫隙設計［J］.推進技術，2005，26(5):448-451.

［8］田四朋，唐國金，李道奎，等.固體火箭發(fā)動機噴管結構完整性分析［J］.固體火箭技術，2005，28(3):180-183.

［9］胡海峰，鮑福廷，王藝杰，等.噴管分離流流動-熱-結構順序耦合數(shù)值模擬及試驗研究［J］.宇航學報，2011，32(7):1534-1542.

［10］鄭曉亞，陳鳳明，蔡飛超.復合噴管結構耦合計算的一種策略［J］.航空動力學報，2011，26(1):223-227.

［11］黃宏艷，王強.軸對稱噴管內(nèi)外流場與結構溫度場耦合計算［J］.推進技術，2008，29(2):195-200.

［12］陳汝訓.固體火箭發(fā)動機設計與研究(下)［M］.北京:宇航出版社，1993.

(編輯:崔賢彬)

Coupled numerical simulation on flow field and thermo-structure of composite nozzle

FU Xue-jin1，2，DENG Heng2，KOU Jun-qiang2，SONG Jun2，CHEN Gang3
(1.Second Artillery Engineering University，Xi＇an 710025，China;2.Xi＇an Modern Control Technology Research Institute，Xi＇an 710065，China;3.Equipment Representive Office，Second Artillery Corps，Stationad in China Academy of Launch Vehicle Technology，Beijing 100076，China)

In order to estimate the heat protection performance and obtain the temperature load for erosion calculation and structural stress analysis，temperature numerical simulation of composite nozzle thermo-structure was implemented by the Fluent CFD software using the fluid-structure coupled method.In the computation，the RNG k-ω turbulence model and the enhanced wall treatment method were selected.The temperature distributions at the end of working time and its evolution in initial period were analyzed，and the feature of nozzle throat temperature field with time-varying was investigated.The throat erosion was also estimated by diffusion control model.The simulation results show that the thermo-structure meets the heat protection requirement.The highest temperature zone is located in the middle and posterior portion of nozzle convergent section.The linear erosion thickness of nozzle throat is about 2.1 mm.The results also offer technical support for further optimizing the design of the nozzle thermo-structure.

coupled fluid-structure;temperature field;composite nozzle;numerical simulation

V438

A

1006-2793(2014)06-0781-06

10.7673/j.issn.1006-2793.2014.06.008

2014-03-21;

2014-06-29。

傅學金(1978—)，男，博士，研究方向為固體火箭發(fā)動機結構數(shù)值仿真。E-mail:fxej_epgc@sina.com