薛海峰，陈 雄，周长省

(南京理工大学 机械工程学院，南京 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