基于熱解過程的變熱物性碳/酚醛能量擴散數值研究①

薛海峰，陳 雄，周長省

(南京理工大學 機械工程學院，南京 210094)

?

基于熱解過程的變熱物性碳/酚醛能量擴散數值研究①

薛海峰，陳 雄，周長省

(南京理工大學 機械工程學院，南京 210094)

為了研究碳/酚醛材料內部能量傳遞以及體積燒蝕過程，基于熱解動力學模型，提出了碳/酚醛復合材料熱物性隨溫度和時間的變化模型。通過熱解過程中材料本身不斷發生變化的密度來反推酚醛樹脂、炭纖維、樹脂碳以及材料孔隙的體積比，以此來推斷材料的瞬態物性參數。在該前提下，常用的碳/酚醛三層模型中的分層結構可在程序內部通過對密度的判定來獲取，實現了傳熱燒蝕的耦合計算。研究結果表明，在受熱初期，熱解層厚度及材料質量損失速率迅速增高；隨著時間的推進，能量逐漸向材料內部進入，并在進入過程中同樣由于熱解吸熱、氣體逸出以及對外界熱輻射在逐漸衰減，使得能量滲透速度減緩；仿真結果與氮氣氛圍下的激光燒蝕試驗結果吻合較好。

碳/酚醛；變熱物性；能量擴散；體積燒蝕

0 引言

相對精確的傳熱與燒蝕數值模擬是建立在準確的材料熱物性描述基礎上的。熱解炭化類材料熱物性的實驗室獲取基本針對完全炭化狀態，或者原始材料狀態，不同狀態下的參數很難由實驗獲取[1]。目前，對熱解炭化類材料傳熱燒蝕研究基本建立在常物性[2]或者分層常物性前提下的[3-5]，研究者預設了炭化層、熱解層及原始材料層，通過分別求解各層的能量方程來獲取溫度分布；還有一部分研究建立在熱物性隨溫度變化而變化的基礎之上[6-7]。事實上，碳/酚醛材料達到熱解臨界溫度會發生熱解，導致材料本身發生變化，故而溫度是材料熱物性變化因素之一；作為化學反應，熱解需要時間的作用，那么材料熱物性不僅與溫度有關，還與高溫環境下的作用時間有關。

本文研究提出了碳/酚醛復合材料密度、熱導率和比熱容等熱物性隨溫度和時間變化模型，并以此模型對第二類邊界條件下碳/酚醛材料一維溫度擴散進行了求解；在該前提下，常用的碳/酚醛三層模型中的分層結構可通過對密度的判定來獲取，并不需要對傳熱過程進行分層求解。通過對氮氣氛圍下激光燒蝕試樣的分析，證實了相關假設與模型的準確性。研究了材料內部熱物性變化以及熱解燒蝕過程隨溫度及時間的變化規律，為相關仿真設計提供了參考。

1 物理模型和計算方法

1.1 基本假設

(1)熱解氣體是化學惰性，熱解氣體密度與固體材料相比為小量，并與固相材料處于熱平衡狀態，且假設熱解氣體生成后，在當前存在于材料孔隙中，下一個時刻被新生成的熱解氣體替代并逸出固體材料區域。

(2)由于熱解導致的體積變化及熱膨脹忽略不計，不考慮材料表面退移現象。

(3)在工作過程中，炭纖維不發生化學反應，且不存在升華過程。

(4)假設材料的熱導率為各向同性。

1.2 物理模型

碳/酚醛復合材料的耐高溫機制已經研究的非常明確，按照燒蝕機理將其劃分為熱解炭化類材料。在高溫邊界作用下，熱量通過熱傳導作用進入碳/酚醛材料內部；盡管材料表面迅速升溫，但在一開始溫度并未達到臨界熱解溫度。在很短的時間內，材料表面及靠近表面的內部區域溫度達到了臨界熱解溫度，開始出現熱解現象。隨著工作時間的推移，材料表面附近持續熱解直至炭化。在這個過程中，當地密度逐漸降低，直至炭化層密度。至此，材料表面形成炭化層。

圖1 碳/酚醛復合材料三層模型示意圖Fig.1 Schematic of three layer model about carbon-phenolic

根據碳/酚醛復合材料的傳熱燒蝕機理，通用的三層模型如圖1所示，包括Ⅰ炭化層、Ⅱ熱解層和Ⅲ原始材料層。羅永康[8]的研究表明，在2 800 K溫度以下，炭化形成的酚醛樹脂碳主要成分為無定形碳(Free Carbon)。結合本文的假設：炭化層中的主要構成為無定形碳、炭纖維以及酚醛樹脂熱解形成的孔隙；熱解層中的主要構成為酚醛樹脂、部分酚醛樹脂熱解形成的無定形碳、炭纖維以及熱解形成的孔隙；原始材料層中主要成分為酚醛樹脂和炭纖維。

1.3 控制方程

本文的研究對象為圓柱體試樣在自然對流環境下的傳熱燒蝕過程。已知加熱壁面的熱流密度，故熱壁面采用第二類邊界條件，冷壁面定為絕熱壁。由于研究對象的工作環境為自然對流環境，其換熱系數量級為5～25 W/(m2·K)。對柱體而言Bi數為

式中h為表面換熱系數；r為半徑；k為材料熱導率。

本研究對象的Bi量級約為0.009～0.044，所以試樣同一軸向剖面處各點過余溫度偏差小于5%[9]。前文中假設碳/酚醛材料在工作過程中熱導率為各向同性，故而該研究最終簡化為一維非定常熱擴散問題：

(1)

其中，源項S包含了酚醛樹脂熱解潛熱、熱解氣體逸出所攜帶的能量以及試樣側表面對外界的熱輻射，材料的當地密度ρ、比定壓熱容cp和熱導率k均為當地溫度T和工作時間t的函數。

邊界條件：

x=0，q=qhot

x=L，q=0

1.4 熱物性及能量源項處理

1.4.1 密度模型

基于“多組分模型[10]”給出碳/酚醛原始材料密度表達式：

ρv=ΓρΑ+(1-Γ)ρB

(2)

式中Γ為酚醛數值在原始材料中所占體積比；下標v、A、B分別代表原始材料、酚醛樹脂和炭纖維。

在熱解過程中，以阿累尼烏斯方程來表征酚醛樹脂的熱解速率[11]，材料密度的變化過程為

且ρ≥ρch

(3)

式中T0為碳/酚醛熱解臨界溫度；A0為熱解反應指前因子；E0為熱解反應活化能；下標ch代表炭化層。

在當前時間步內積分式(3)，且在該時間步內，視溫度T為常量，有

(4)

式中 下標1和2分別代表當前時間步和下個時間步。

當材料密度值達到炭化層密度ρch時，熱解過程結束，當地材料轉化為炭化層。文中炭化層密度ρch取決于材料初始狀態的配方以及酚醛樹脂的成碳率，由式(4)給出：

ρch=εcharΓρA+(1-Γ)ρB

(5)

式中εchar為酚醛樹脂的成碳率。

1.4.2 比熱容模型

酚醛樹脂在較高溫度下分解成為小分子氣體逸出并留下殘碳[11]。本文假設將酚醛樹脂的熱解過程等效為部分樹脂瞬間熱解成無定形碳和小分子氣體，其余保持原始材料結構。在該假設條件下，熱解層中含有的材料主要為酚醛樹脂、炭纖維、由部分酚醛樹脂轉變成的無定形碳和材料孔隙。其中，存在于材料孔隙中的熱解氣體處于流動狀態，其比熱容在能量源項中考慮。設熱解區域中轉變為無定形碳的酚醛樹脂占當地當時原始酚醛樹脂的質量分數為τA→C，在該假設前提下，有

ρ=τA→CεcharΓρA+(1-τA→C)ΓρA+(1-Γ)ρB

結合式(2)，即

ρ=ρV-(1-εchar)τA→CΓρA

(6)

則熱解層中當地酚醛樹脂轉化率為

(7)

若材料處于初始狀態，有τA→C=0；材料處于完全炭化狀態，有τA→C=1。

該狀態下，樹脂碳、酚醛樹脂和炭纖維占當地熱解層中的質量分數分別為

(8)

同樣，根據“多組分模型”可給出碳/酚醛材料在整個工作過程中的比定壓熱容：

cp=τAcpA+τBcpB+τCcpC

(9)

式中 下標C代表樹脂碳。

1.4.3 熱導率模型

由于炭化層和熱解層中存在大量的孔隙，孔隙產生有兩個來源：一是酚醛樹脂熱解，二是炭纖維和酚醛樹脂由于高溫導致的粘合面剝離；而酚醛樹脂熱解產生的孔隙是炭化層和熱解層中孔隙的主要來源，假設酚醛樹脂完全熱解產生的無定形碳孔隙率為φC，有

(10)

基于多孔介質傳熱理論[12]中“有效導熱系數法”，給出碳/酚醛在不同狀態下的熱導率表達方程：

k=ΓAkA+ΓBkB+ΓCkC+φkg

(11)

(12)

式中 下標g代表熱解氣體。

1.4.4 能量源項

能量源項S(W/m3)包含了酚醛樹脂熱解潛熱、熱解氣體逸出所攜帶的能量以及試樣側表面對外界的熱輻射的等效源項：

S=Spy+Sen+Srad

(13)

其中

式中Spy為酚醛樹脂熱解能量源項；Sen為熱解氣體逸出攜帶的能量源項；Srad為側表面對外界熱輻射的等效源項；ε為材料表面發射率；σ為斯特潘·波爾茲曼常數；hA為酚醛樹脂熱解潛熱；hg為熱解氣體顯焓。

1.4.5 熱解氣體熱物性模型

碳/酚醛材料在工作過程中，溫度變化較大，這會導致前文推導過程中用到的熱解氣體導熱系數以及比熱容變化較大，為保證數值計算的準確性，需要對其進行估算。

相關研究表明，酚醛樹脂熱解產物主要成分為：甲烷(CH4)、乙烯(C2H4)、乙炔(C2H2)和苯蒸氣(C6H6)。表1給出了酚醛樹脂高溫熱解產物分布[11]。

表1 酚醛樹脂熱解產物摩爾分數Table 1 Mole fraction of pyrolysis products

本文忽略壓強變化對熱解氣體熱物性的影響。查閱JANAF表，可獲得表1中各組分在不同溫度下比定壓熱容(J·kg-1·K-1)，進行分段多項式擬合如下形式：

cp=a+bT+cT2+dT3+eT4

具體結果見表2。根據多原子氣體的EuckenA關系式[13]及表2中定壓比熱函數可得出各熱解組分熱導率隨溫度變化關系如下形式：

式中T、M、σ、Ωk、cp和R0分別為溫度(K)、摩爾質量(g/mol)、碰撞直徑(?)、倫納德-瓊斯參數、比定壓熱容(J/(kg·K))、標準氣體常數(J/(mol·K))。

關于各組分的L-J參數及碰撞積分可查表獲取。在獲取了各組分的定壓比熱和熱導率之后，分別通過質量平均和多組分氣體混合物熱導率公式，求得混合氣體的定壓比熱及熱導率。

通過以上的熱物性模型，可將炭化層、熱解層以及原始材料層的傳熱方程進行統一，不需要對方程進行分層計算，簡化了數學模型。

表2 各組分比定壓熱容多段擬合系數表Table 2 Coefficient of polynomial fitting about specific heat of each component

2 計算結果及分析

本文使用上述方法及參數對壁面熱流密度為1.402 MW/m2圓柱體碳/酚醛試樣的熱解及熱擴散過程進行了數值模擬。試樣直徑為6.84 mm，長度為15 mm。碳/酚醛材料的相關參數由合作單位給出：酚醛樹脂密度為1 186 kg/m3，炭纖維密度為1 800 kg/m3，酚醛樹脂的體積分數為0.6，酚醛樹脂的成碳率為0.5，酚醛樹脂完全炭化后的真實密度為1 500 kg/m3；熱解反應臨界溫度為573 K，反應指前因子為185 000 s-1，活化能為100 810 J/mol[11]，熱解潛熱為420 000 J/kg[14]；材料初始溫度與環境溫度保持一致為300 K，材料表面發射率為0.85[6]。作為對比，采用了一種常用的依據溫度作線性插值獲取熱解層內材料物性參數[14]方法來計算材料內部能量傳遞過程，并對比了2種方法對體積燒蝕過程的影響。

圖2為材料受熱端面溫度隨時間變化曲線。材料表面在受熱后短時間內迅速上升，在5 ms左右材料表面即達到了熱解臨界溫度，此刻表面開始發生熱解反應；隨著時間的推移，材料表面溫度在1 s時已經達到了1 284 K，約為最高溫度的65%；在以后的時間內，材料表面溫度上升緩慢，直至10 s材料表面基本達到該工況下的最高溫度1 959 K。

圖3為材料內部溫度在不同時刻點的分布曲線。熱流加載初期,材料表面很快達到了最高溫度，這在圖2中已經有所體現。溫度在材料內部的傳遞隨著時間的推進而逐漸變緩,這主要是由于酚醛樹脂的熱解吸熱、熱解氣體逸出過程中攜帶了部分熱量、材料受熱溫度升高對外界的熱輻射隨著溫度的升高呈幾何級數增長以及材料本身的熱容所造成的。在材料受熱過程后期,炭化層內的溫度分布趨于一致,即達到了一個近似定常的能量傳遞過程。

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

圖3 材料內部溫度分布Fig.3 Temperature distribution inside material

因為材料表面最高溫度約為1 595 K，并未達到炭纖維以及無定形碳的升華點；且研究對象處于非氧化氛圍，故而材料表面并不存在熱化學燒蝕過程。

圖4為材料內部在不同時刻點密度分布曲線，在本研究的物理模型中，材料密度是表征材料分層狀態的一個特征參數。圖5、圖6分別更為詳細地給出了炭化層厚度變化和熱解層厚度變化曲線。在熱流加載初期，材料表面很快發生熱解，且熱解層厚度隨時間近乎線性增大，且增速較高；在1 s之前材料表面一直處于熱解狀態，并未形成炭化層。1 s時，材料表面開始出現炭化層，此刻熱解層厚度為0.92 mm。由此往后，熱解層厚度增速逐漸降低，而炭化層厚度以較高增速開始增加。隨著熱流作用時間的推進，炭化層以及熱解層厚度基本以一個較低的增速均勻增加。這主要是因為能量在往材料內部傳遞過程中，由于材料表面對外界的熱輻射以及內部熱解吸熱等行為，導致能量衰減，故而炭化層厚度的增速在降低。

圖4 材料內部密度分布Fig.4 Density distribution inside material

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

圖6 熱解層厚度Fig.6 Thickness of pyrolysis layer

圖7為材料內部熱解氣體生成率分布，材料內部熱解氣體生成率是當地熱解反應速率的直接表征。隨著時間推進，熱解層逐漸向材料內部移動，結合圖3可發現，熱解層內當地溫度逐漸降低，這直接導致了熱解氣體生成率峰值的下降。

圖7 熱解氣體生成率分布Fig.7 Local generation rate of pyrolysis gas

根據阿累尼烏斯方程，本研究中的熱解速率與當地材料密度和當地溫度相關。圖8給出了t=30 s時，材料熱解層內溫度、密度以及熱解氣體生成速率分布曲線。熱解氣體生成速率在熱解層內近乎呈正態分布，但仔細觀察可發現，熱解層在靠近受熱面處氣體生成率變化較快，而在遠離受熱面處氣體生成率變化較為平緩。造成這個現象的原因首先是由于靠近受熱面處材料溫度較高，導致化學反應速度較快；其次是由于密度對熱解氣體生成速率的貢獻為線性的，而溫度的貢獻卻是指數形式的。這說明在熱解過程中盡管密度與溫度對熱解氣體生成速率的影響是綜合的，但溫度仍舊占主導地位。

圖8 溫度、密度和熱解氣體生成率分布(t=30 s)Fig.8 Distribution of temperature,density and local generation rate of pyrolysis gas (t=30 s)

圖9為熱解氣體生成總通量變化曲線。顯然，在受熱初期，熱解氣體生成總通量快速增長，在1.06 s時達到了0.1467 kg/(m2·s)，隨后即很快下降，10 s之后，進入緩慢下降區。

為了對本研究提出的模型以及程序進行驗證，在氮氣氛圍中對碳/酚醛試樣進行激光燒蝕試驗。試驗一共分為6組，燒蝕時間從10～60 s，以10 s為一個時間間隔；對燒蝕前后的試樣進行稱重分析，獲取了每組試驗試樣的質量損失。圖10給出了兩種模型試樣質量損失仿真結果與試驗結果的對比。由于文獻[15]中根據溫度線性插值材料物性模型，并未考慮熱解反應速率的影響，導致仿真結果與試驗結果偏差較大；而本文提出的模型仿真與試驗結果有較好的一致性，說明了模型的正確性。

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

圖10 結構質量損失Fig.10 Mass loss of samples

3 結論

(1)在受熱初期，材料表面溫度以及熱解氣體生成速率迅速升高；而后的過程中，由于材料熱解吸熱以及氣體逸出攜帶的能量，材料表面的溫升速率開始降低。

(2) 隨著時間的推進，能量逐漸向材料內部進入，并在進入過程中同樣由于熱解吸熱、氣體逸出，以及對外界熱輻射在逐漸衰減，使得能量滲透速度減緩，該過程符合材料的熱防護需求。

(3) 溫度對熱解過程的貢獻呈指數形式，在熱解層內，溫度仍占熱解主導地位。

(4) 仿真結果與氮氣氛圍下的激光燒蝕試驗結果吻合較好，對后續研究有一定的指導意義。

[1] 張杰,孫冰.基于熱解動力學的絕熱材料燒蝕研究[J].固體火箭技術,2010,33(4):454-458.

[2] 王作良.噴管熱防護的有限元數值分析[D].哈爾濱:哈爾濱工程大學,2005.

[3] 陳蘭,余澤楚.燒蝕材料三維瞬態熱響應數值模擬[J].工程熱物理學報,1996,17(12):103-106.

[4] 王春光,田維平,楊德敏,等.脈沖發動機中隔層傳熱炭化模型[J].推進技術,2012,33(2):259-262.

[5] 張曉光,劉宇,王長輝.雙脈沖固體發動機噴管傳熱燒蝕特性[J].航空動力學報,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] 羅永康,彭維舟,王為民.燒蝕復合材料用酚醛樹脂研究[J].宇航材料工藝,1988(5).

[9] 楊世銘,陶文銓.傳熱學[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].重慶:重慶大學,2007.

[12] 林瑞泰.多孔介質傳熱傳質引論[M].北京:科學出版社,1995.

[13] 陳軍.火箭發動機燃燒基礎[M].南京:南京理工大學,2011.

[14] 徐曉亮.熱防護機理與燒蝕鈍體繞流的渦方法研究[D].北京交通大學,2011.

[15] 梁華.沖壓發動機補燃室內絕熱層傳熱燒蝕特性研究[D].南京理工大學,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—)，男，博士生，研究領域為固體火箭發動機熱防護。E-mail：liangwangongli@163.com

V258

A

1006-2793(2015)01-0130-06

10.7673/j.issn.1006-2793.2015.01.025