薛海峰,陳 雄,周長省
(南京理工大學(xué) 機械工程學(xué)院,南京 210094)
?
基于熱解過程的變熱物性碳/酚醛能量擴散數(shù)值研究①
薛海峰,陳 雄,周長省
(南京理工大學(xué) 機械工程學(xué)院,南京 210094)
為了研究碳/酚醛材料內(nèi)部能量傳遞以及體積燒蝕過程,基于熱解動力學(xué)模型,提出了碳/酚醛復(fù)合材料熱物性隨溫度和時間的變化模型。通過熱解過程中材料本身不斷發(fā)生變化的密度來反推酚醛樹脂、炭纖維、樹脂碳以及材料孔隙的體積比,以此來推斷材料的瞬態(tài)物性參數(shù)。在該前提下,常用的碳/酚醛三層模型中的分層結(jié)構(gòu)可在程序內(nèi)部通過對密度的判定來獲取,實現(xiàn)了傳熱燒蝕的耦合計算。研究結(jié)果表明,在受熱初期,熱解層厚度及材料質(zhì)量損失速率迅速增高;隨著時間的推進,能量逐漸向材料內(nèi)部進入,并在進入過程中同樣由于熱解吸熱、氣體逸出以及對外界熱輻射在逐漸衰減,使得能量滲透速度減緩;仿真結(jié)果與氮氣氛圍下的激光燒蝕試驗結(jié)果吻合較好。
碳/酚醛;變熱物性;能量擴散;體積燒蝕
相對精確的傳熱與燒蝕數(shù)值模擬是建立在準(zhǔn)確的材料熱物性描述基礎(chǔ)上的。熱解炭化類材料熱物性的實驗室獲取基本針對完全炭化狀態(tài),或者原始材料狀態(tài),不同狀態(tài)下的參數(shù)很難由實驗獲取[1]。目前,對熱解炭化類材料傳熱燒蝕研究基本建立在常物性[2]或者分層常物性前提下的[3-5],研究者預(yù)設(shè)了炭化層、熱解層及原始材料層,通過分別求解各層的能量方程來獲取溫度分布;還有一部分研究建立在熱物性隨溫度變化而變化的基礎(chǔ)之上[6-7]。事實上,碳/酚醛材料達(dá)到熱解臨界溫度會發(fā)生熱解,導(dǎo)致材料本身發(fā)生變化,故而溫度是材料熱物性變化因素之一;作為化學(xué)反應(yīng),熱解需要時間的作用,那么材料熱物性不僅與溫度有關(guān),還與高溫環(huán)境下的作用時間有關(guān)。
本文研究提出了碳/酚醛復(fù)合材料密度、熱導(dǎo)率和比熱容等熱物性隨溫度和時間變化模型,并以此模型對第二類邊界條件下碳/酚醛材料一維溫度擴散進行了求解;在該前提下,常用的碳/酚醛三層模型中的分層結(jié)構(gòu)可通過對密度的判定來獲取,并不需要對傳熱過程進行分層求解。通過對氮氣氛圍下激光燒蝕試樣的分析,證實了相關(guān)假設(shè)與模型的準(zhǔn)確性。研究了材料內(nèi)部熱物性變化以及熱解燒蝕過程隨溫度及時間的變化規(guī)律,為相關(guān)仿真設(shè)計提供了參考。
1.1 基本假設(shè)
(1)熱解氣體是化學(xué)惰性,熱解氣體密度與固體材料相比為小量,并與固相材料處于熱平衡狀態(tài),且假設(shè)熱解氣體生成后,在當(dāng)前存在于材料孔隙中,下一個時刻被新生成的熱解氣體替代并逸出固體材料區(qū)域。
(2)由于熱解導(dǎo)致的體積變化及熱膨脹忽略不計,不考慮材料表面退移現(xiàn)象。
(3)在工作過程中,炭纖維不發(fā)生化學(xué)反應(yīng),且不存在升華過程。
(4)假設(shè)材料的熱導(dǎo)率為各向同性。
1.2 物理模型
碳/酚醛復(fù)合材料的耐高溫機制已經(jīng)研究的非常明確,按照燒蝕機理將其劃分為熱解炭化類材料。在高溫邊界作用下,熱量通過熱傳導(dǎo)作用進入碳/酚醛材料內(nèi)部;盡管材料表面迅速升溫,但在一開始溫度并未達(dá)到臨界熱解溫度。在很短的時間內(nèi),材料表面及靠近表面的內(nèi)部區(qū)域溫度達(dá)到了臨界熱解溫度,開始出現(xiàn)熱解現(xiàn)象。隨著工作時間的推移,材料表面附近持續(xù)熱解直至炭化。在這個過程中,當(dāng)?shù)孛芏戎饾u降低,直至炭化層密度。至此,材料表面形成炭化層。

圖1 碳/酚醛復(fù)合材料三層模型示意圖Fig.1 Schematic of three layer model about carbon-phenolic
根據(jù)碳/酚醛復(fù)合材料的傳熱燒蝕機理,通用的三層模型如圖1所示,包括Ⅰ炭化層、Ⅱ熱解層和Ⅲ原始材料層。羅永康[8]的研究表明,在2 800 K溫度以下,炭化形成的酚醛樹脂碳主要成分為無定形碳(Free Carbon)。結(jié)合本文的假設(shè):炭化層中的主要構(gòu)成為無定形碳、炭纖維以及酚醛樹脂熱解形成的孔隙;熱解層中的主要構(gòu)成為酚醛樹脂、部分酚醛樹脂熱解形成的無定形碳、炭纖維以及熱解形成的孔隙;原始材料層中主要成分為酚醛樹脂和炭纖維。
1.3 控制方程
本文的研究對象為圓柱體試樣在自然對流環(huán)境下的傳熱燒蝕過程。已知加熱壁面的熱流密度,故熱壁面采用第二類邊界條件,冷壁面定為絕熱壁。由于研究對象的工作環(huán)境為自然對流環(huán)境,其換熱系數(shù)量級為5~25 W/(m2·K)。對柱體而言Bi數(shù)為
式中h為表面換熱系數(shù);r為半徑;k為材料熱導(dǎo)率。
本研究對象的Bi量級約為0.009~0.044,所以試樣同一軸向剖面處各點過余溫度偏差小于5%[9]。前文中假設(shè)碳/酚醛材料在工作過程中熱導(dǎo)率為各向同性,故而該研究最終簡化為一維非定常熱擴散問題:
(1)
其中,源項S包含了酚醛樹脂熱解潛熱、熱解氣體逸出所攜帶的能量以及試樣側(cè)表面對外界的熱輻射,材料的當(dāng)?shù)孛芏圈?、比定壓熱容cp和熱導(dǎo)率k均為當(dāng)?shù)販囟萒和工作時間t的函數(shù)。
邊界條件:
x=0,q=qhot
x=L,q=0
1.4 熱物性及能量源項處理
1.4.1 密度模型
基于“多組分模型[10]”給出碳/酚醛原始材料密度表達(dá)式:
ρv=ΓρΑ+(1-Γ)ρB
(2)
式中Γ為酚醛數(shù)值在原始材料中所占體積比;下標(biāo)v、A、B分別代表原始材料、酚醛樹脂和炭纖維。
在熱解過程中,以阿累尼烏斯方程來表征酚醛樹脂的熱解速率[11],材料密度的變化過程為
且ρ≥ρch
(3)
式中T0為碳/酚醛熱解臨界溫度;A0為熱解反應(yīng)指前因子;E0為熱解反應(yīng)活化能;下標(biāo)ch代表炭化層。
在當(dāng)前時間步內(nèi)積分式(3),且在該時間步內(nèi),視溫度T為常量,有
(4)
式中 下標(biāo)1和2分別代表當(dāng)前時間步和下個時間步。
當(dāng)材料密度值達(dá)到炭化層密度ρch時,熱解過程結(jié)束,當(dāng)?shù)夭牧限D(zhuǎn)化為炭化層。文中炭化層密度ρch取決于材料初始狀態(tài)的配方以及酚醛樹脂的成碳率,由式(4)給出:
ρch=εcharΓρA+(1-Γ)ρB
(5)
式中εchar為酚醛樹脂的成碳率。
1.4.2 比熱容模型
酚醛樹脂在較高溫度下分解成為小分子氣體逸出并留下殘?zhí)糩11]。本文假設(shè)將酚醛樹脂的熱解過程等效為部分樹脂瞬間熱解成無定形碳和小分子氣體,其余保持原始材料結(jié)構(gòu)。在該假設(shè)條件下,熱解層中含有的材料主要為酚醛樹脂、炭纖維、由部分酚醛樹脂轉(zhuǎn)變成的無定形碳和材料孔隙。其中,存在于材料孔隙中的熱解氣體處于流動狀態(tài),其比熱容在能量源項中考慮。設(shè)熱解區(qū)域中轉(zhuǎn)變?yōu)闊o定形碳的酚醛樹脂占當(dāng)?shù)禺?dāng)時原始酚醛樹脂的質(zhì)量分?jǐn)?shù)為τA→C,在該假設(shè)前提下,有
ρ=τA→CεcharΓρA+(1-τA→C)ΓρA+(1-Γ)ρB
結(jié)合式(2),即
ρ=ρV-(1-εchar)τA→CΓρA
(6)
則熱解層中當(dāng)?shù)胤尤渲D(zhuǎn)化率為
(7)
若材料處于初始狀態(tài),有τA→C=0;材料處于完全炭化狀態(tài),有τA→C=1。
該狀態(tài)下,樹脂碳、酚醛樹脂和炭纖維占當(dāng)?shù)責(zé)峤鈱又械馁|(zhì)量分?jǐn)?shù)分別為
(8)
同樣,根據(jù)“多組分模型”可給出碳/酚醛材料在整個工作過程中的比定壓熱容:
cp=τAcpA+τBcpB+τCcpC
(9)
式中 下標(biāo)C代表樹脂碳。
1.4.3 熱導(dǎo)率模型
由于炭化層和熱解層中存在大量的孔隙,孔隙產(chǎn)生有兩個來源:一是酚醛樹脂熱解,二是炭纖維和酚醛樹脂由于高溫導(dǎo)致的粘合面剝離;而酚醛樹脂熱解產(chǎn)生的孔隙是炭化層和熱解層中孔隙的主要來源,假設(shè)酚醛樹脂完全熱解產(chǎn)生的無定形碳孔隙率為φC,有
(10)
基于多孔介質(zhì)傳熱理論[12]中“有效導(dǎo)熱系數(shù)法”,給出碳/酚醛在不同狀態(tài)下的熱導(dǎo)率表達(dá)方程:
k=ΓAkA+ΓBkB+ΓCkC+φkg
(11)
(12)
式中 下標(biāo)g代表熱解氣體。
1.4.4 能量源項
能量源項S(W/m3)包含了酚醛樹脂熱解潛熱、熱解氣體逸出所攜帶的能量以及試樣側(cè)表面對外界的熱輻射的等效源項:
S=Spy+Sen+Srad
(13)
其中
式中Spy為酚醛樹脂熱解能量源項;Sen為熱解氣體逸出攜帶的能量源項;Srad為側(cè)表面對外界熱輻射的等效源項;ε為材料表面發(fā)射率;σ為斯特潘·波爾茲曼常數(shù);hA為酚醛樹脂熱解潛熱;hg為熱解氣體顯焓。
1.4.5 熱解氣體熱物性模型
碳/酚醛材料在工作過程中,溫度變化較大,這會導(dǎo)致前文推導(dǎo)過程中用到的熱解氣體導(dǎo)熱系數(shù)以及比熱容變化較大,為保證數(shù)值計算的準(zhǔn)確性,需要對其進行估算。
相關(guān)研究表明,酚醛樹脂熱解產(chǎn)物主要成分為:甲烷(CH4)、乙烯(C2H4)、乙炔(C2H2)和苯蒸氣(C6H6)。表1給出了酚醛樹脂高溫?zé)峤猱a(chǎn)物分布[11]。

表1 酚醛樹脂熱解產(chǎn)物摩爾分?jǐn)?shù)Table 1 Mole fraction of pyrolysis products
本文忽略壓強變化對熱解氣體熱物性的影響。查閱JANAF表,可獲得表1中各組分在不同溫度下比定壓熱容(J·kg-1·K-1),進行分段多項式擬合如下形式:
cp=a+bT+cT2+dT3+eT4
具體結(jié)果見表2。根據(jù)多原子氣體的EuckenA關(guān)系式[13]及表2中定壓比熱函數(shù)可得出各熱解組分熱導(dǎo)率隨溫度變化關(guān)系如下形式:
式中T、M、σ、Ωk、cp和R0分別為溫度(K)、摩爾質(zhì)量(g/mol)、碰撞直徑(?)、倫納德-瓊斯參數(shù)、比定壓熱容(J/(kg·K))、標(biāo)準(zhǔn)氣體常數(shù)(J/(mol·K))。
關(guān)于各組分的L-J參數(shù)及碰撞積分可查表獲取。在獲取了各組分的定壓比熱和熱導(dǎo)率之后,分別通過質(zhì)量平均和多組分氣體混合物熱導(dǎo)率公式,求得混合氣體的定壓比熱及熱導(dǎo)率。
通過以上的熱物性模型,可將炭化層、熱解層以及原始材料層的傳熱方程進行統(tǒng)一,不需要對方程進行分層計算,簡化了數(shù)學(xué)模型。

表2 各組分比定壓熱容多段擬合系數(shù)表Table 2 Coefficient of polynomial fitting about specific heat of each component
本文使用上述方法及參數(shù)對壁面熱流密度為1.402 MW/m2圓柱體碳/酚醛試樣的熱解及熱擴散過程進行了數(shù)值模擬。試樣直徑為6.84 mm,長度為15 mm。碳/酚醛材料的相關(guān)參數(shù)由合作單位給出:酚醛樹脂密度為1 186 kg/m3,炭纖維密度為1 800 kg/m3,酚醛樹脂的體積分?jǐn)?shù)為0.6,酚醛樹脂的成碳率為0.5,酚醛樹脂完全炭化后的真實密度為1 500 kg/m3;熱解反應(yīng)臨界溫度為573 K,反應(yīng)指前因子為185 000 s-1,活化能為100 810 J/mol[11],熱解潛熱為420 000 J/kg[14];材料初始溫度與環(huán)境溫度保持一致為300 K,材料表面發(fā)射率為0.85[6]。作為對比,采用了一種常用的依據(jù)溫度作線性插值獲取熱解層內(nèi)材料物性參數(shù)[14]方法來計算材料內(nèi)部能量傳遞過程,并對比了2種方法對體積燒蝕過程的影響。
圖2為材料受熱端面溫度隨時間變化曲線。材料表面在受熱后短時間內(nèi)迅速上升,在5 ms左右材料表面即達(dá)到了熱解臨界溫度,此刻表面開始發(fā)生熱解反應(yīng);隨著時間的推移,材料表面溫度在1 s時已經(jīng)達(dá)到了1 284 K,約為最高溫度的65%;在以后的時間內(nèi),材料表面溫度上升緩慢,直至10 s材料表面基本達(dá)到該工況下的最高溫度1 959 K。
圖3為材料內(nèi)部溫度在不同時刻點的分布曲線。熱流加載初期,材料表面很快達(dá)到了最高溫度,這在圖2中已經(jīng)有所體現(xiàn)。溫度在材料內(nèi)部的傳遞隨著時間的推進而逐漸變緩,這主要是由于酚醛樹脂的熱解吸熱、熱解氣體逸出過程中攜帶了部分熱量、材料受熱溫度升高對外界的熱輻射隨著溫度的升高呈幾何級數(shù)增長以及材料本身的熱容所造成的。在材料受熱過程后期,炭化層內(nèi)的溫度分布趨于一致,即達(dá)到了一個近似定常的能量傳遞過程。

圖2 材料表面溫度Fig.2 Temperature of surface

圖3 材料內(nèi)部溫度分布Fig.3 Temperature distribution inside material
因為材料表面最高溫度約為1 595 K,并未達(dá)到炭纖維以及無定形碳的升華點;且研究對象處于非氧化氛圍,故而材料表面并不存在熱化學(xué)燒蝕過程。
圖4為材料內(nèi)部在不同時刻點密度分布曲線,在本研究的物理模型中,材料密度是表征材料分層狀態(tài)的一個特征參數(shù)。圖5、圖6分別更為詳細(xì)地給出了炭化層厚度變化和熱解層厚度變化曲線。在熱流加載初期,材料表面很快發(fā)生熱解,且熱解層厚度隨時間近乎線性增大,且增速較高;在1 s之前材料表面一直處于熱解狀態(tài),并未形成炭化層。1 s時,材料表面開始出現(xiàn)炭化層,此刻熱解層厚度為0.92 mm。由此往后,熱解層厚度增速逐漸降低,而炭化層厚度以較高增速開始增加。隨著熱流作用時間的推進,炭化層以及熱解層厚度基本以一個較低的增速均勻增加。這主要是因為能量在往材料內(nèi)部傳遞過程中,由于材料表面對外界的熱輻射以及內(nèi)部熱解吸熱等行為,導(dǎo)致能量衰減,故而炭化層厚度的增速在降低。

圖4 材料內(nèi)部密度分布Fig.4 Density distribution inside material

圖5 炭化層厚度Fig.5 Thickness of charring layer

圖6 熱解層厚度Fig.6 Thickness of pyrolysis layer
圖7為材料內(nèi)部熱解氣體生成率分布,材料內(nèi)部熱解氣體生成率是當(dāng)?shù)責(zé)峤夥磻?yīng)速率的直接表征。隨著時間推進,熱解層逐漸向材料內(nèi)部移動,結(jié)合圖3可發(fā)現(xiàn),熱解層內(nèi)當(dāng)?shù)販囟戎饾u降低,這直接導(dǎo)致了熱解氣體生成率峰值的下降。

圖7 熱解氣體生成率分布Fig.7 Local generation rate of pyrolysis gas
根據(jù)阿累尼烏斯方程,本研究中的熱解速率與當(dāng)?shù)夭牧厦芏群彤?dāng)?shù)販囟认嚓P(guān)。圖8給出了t=30 s時,材料熱解層內(nèi)溫度、密度以及熱解氣體生成速率分布曲線。熱解氣體生成速率在熱解層內(nèi)近乎呈正態(tài)分布,但仔細(xì)觀察可發(fā)現(xiàn),熱解層在靠近受熱面處氣體生成率變化較快,而在遠(yuǎn)離受熱面處氣體生成率變化較為平緩。造成這個現(xiàn)象的原因首先是由于靠近受熱面處材料溫度較高,導(dǎo)致化學(xué)反應(yīng)速度較快;其次是由于密度對熱解氣體生成速率的貢獻為線性的,而溫度的貢獻卻是指數(shù)形式的。這說明在熱解過程中盡管密度與溫度對熱解氣體生成速率的影響是綜合的,但溫度仍舊占主導(dǎo)地位。

圖8 溫度、密度和熱解氣體生成率分布(t=30 s)Fig.8 Distribution of temperature,density and local generation rate of pyrolysis gas (t=30 s)
圖9為熱解氣體生成總通量變化曲線。顯然,在受熱初期,熱解氣體生成總通量快速增長,在1.06 s時達(dá)到了0.1467 kg/(m2·s),隨后即很快下降,10 s之后,進入緩慢下降區(qū)。
為了對本研究提出的模型以及程序進行驗證,在氮氣氛圍中對碳/酚醛試樣進行激光燒蝕試驗。試驗一共分為6組,燒蝕時間從10~60 s,以10 s為一個時間間隔;對燒蝕前后的試樣進行稱重分析,獲取了每組試驗試樣的質(zhì)量損失。圖10給出了兩種模型試樣質(zhì)量損失仿真結(jié)果與試驗結(jié)果的對比。由于文獻[15]中根據(jù)溫度線性插值材料物性模型,并未考慮熱解反應(yīng)速率的影響,導(dǎo)致仿真結(jié)果與試驗結(jié)果偏差較大;而本文提出的模型仿真與試驗結(jié)果有較好的一致性,說明了模型的正確性。

圖9 熱解氣體生成總通量變化Fig.9 Total mass flux of pyrolysis gas

圖10 結(jié)構(gòu)質(zhì)量損失Fig.10 Mass loss of samples
(1)在受熱初期,材料表面溫度以及熱解氣體生成速率迅速升高;而后的過程中,由于材料熱解吸熱以及氣體逸出攜帶的能量,材料表面的溫升速率開始降低。
(2) 隨著時間的推進,能量逐漸向材料內(nèi)部進入,并在進入過程中同樣由于熱解吸熱、氣體逸出,以及對外界熱輻射在逐漸衰減,使得能量滲透速度減緩,該過程符合材料的熱防護需求。
(3) 溫度對熱解過程的貢獻呈指數(shù)形式,在熱解層內(nèi),溫度仍占熱解主導(dǎo)地位。
(4) 仿真結(jié)果與氮氣氛圍下的激光燒蝕試驗結(jié)果吻合較好,對后續(xù)研究有一定的指導(dǎo)意義。
[1] 張杰,孫冰.基于熱解動力學(xué)的絕熱材料燒蝕研究[J].固體火箭技術(shù),2010,33(4):454-458.
[2] 王作良.噴管熱防護的有限元數(shù)值分析[D].哈爾濱:哈爾濱工程大學(xué),2005.
[3] 陳蘭,余澤楚.燒蝕材料三維瞬態(tài)熱響應(yīng)數(shù)值模擬[J].工程熱物理學(xué)報,1996,17(12):103-106.
[4] 王春光,田維平,楊德敏,等.脈沖發(fā)動機中隔層傳熱炭化模型[J].推進技術(shù),2012,33(2):259-262.
[5] 張曉光,劉宇,王長輝.雙脈沖固體發(fā)動機噴管傳熱燒蝕特性[J].航空動力學(xué)報,2012,27(6):1391-1397.
[6] Robert L Potts.Application of integral methods to ablation charring erosion,a review[J].Journal of Spacecraft and Rockets,1995,32(2):200-209.
[7] Shinn-Shyong Tzeng,Ya-Ga Chr.Evolution of microstructure and properties of phenolic resin-based carbon/carbon composites during pyrolysis[J].Materials Chemistry and Physics,2002,73(2-3):162-169.
[8] 羅永康,彭維舟,王為民.燒蝕復(fù)合材料用酚醛樹脂研究[J].宇航材料工藝,1988(5).
[9] 楊世銘,陶文銓.傳熱學(xué)[M].北京:高等教育出版社,1998.
[10] Chen Y K,Milos F S.Two-dimensional implicit thermal response and ablation program for charring materials[J].Journal of Spacecraft and Rockets,2001,38(4):473-481.
[11] 馬偉.酚醛樹脂的熱解研究[D].重慶:重慶大學(xué),2007.
[12] 林瑞泰.多孔介質(zhì)傳熱傳質(zhì)引論[M].北京:科學(xué)出版社,1995.
[13] 陳軍.火箭發(fā)動機燃燒基礎(chǔ)[M].南京:南京理工大學(xué),2011.
[14] 徐曉亮.熱防護機理與燒蝕鈍體繞流的渦方法研究[D].北京交通大學(xué),2011.
[15] 梁華.沖壓發(fā)動機補燃室內(nèi)絕熱層傳熱燒蝕特性研究[D].南京理工大學(xué),2009.
(編輯:薛永利)
Numerical research on energy diffusion of carbon-phenolic with variable thermal properties based on pyrolysis process
XUE Hai-feng,CHEN Xiong,ZHOU Chang-sheng
(School of Mechanical Engineering,Nanjing University of Science and Technology,Nanjing 210094,China)
A one-dimensional heat transfer and volumetric ablation code was developed to study energy diffusion process in carbon-phenolic.Based on pyrolysis kinetics model,a model on variable thermal properties with local temperature and working time was proposed.Transient physical property parameter of material can be deduced using variable densing of material during pyrolysis process to calculate the volume ratio among phenolic resin,carbon fiber,resin carbon and porosity of material.Under the premise,the commonly used carbon-phenolic hierarchical structure of the three layer model can be obtained by the judgement of density within the material.The simulated results show that at the beginning of heating,thickness of pyrolysis layer and mass loss increase rapidly.After the formation of charring,thickness of pyrolysis layer grows slowly and mass loss rate decreases.As heating time goes on,the internal temperature of charring gradually converges,and pyrolysis layer gradually enters into material with a steady speed.
carbon-phenolic;variable thermal properties;energy diffusion;volumetric ablation
2014-09-19;
:2014-10-27。
薛海峰(1986—),男,博士生,研究領(lǐng)域為固體火箭發(fā)動機熱防護。E-mail:liangwangongli@163.com
V258
A
1006-2793(2015)01-0130-06
10.7673/j.issn.1006-2793.2015.01.025