纤维增强复合材料激光烧蚀效应的数值模拟

彭国良，闫 辉，刘 峰，2，杜太焦，2，王玉恒，2

（1.西北核技术研究所，陕西西安710024；2.激光与物质相互作用国家重点实验室，陕西西安710024）

纤维增强复合材料激光烧蚀效应的数值模拟

彭国良1，闫 辉1，刘 峰1，2，杜太焦1，2，王玉恒1，2

（1.西北核技术研究所，陕西西安710024；2.激光与物质相互作用国家重点实验室，陕西西安710024）

考虑材料的热解、氧化、相变及辐射和内外对流换热等物理过程，给出了激光烧蚀纤维增强复合材料的物理模型及数学模型。以碳纤维/环氧树脂复合材料为例，编程计算了材料的激光烧蚀过程，计算结果与实验结果符合得较好。计算结果表明：考虑复合材料的内对流时得到的结果更准确；较强功率密度激光辐照时，氧化对烧蚀的贡献可以忽略；功率密度一定时，烧蚀质量随时间近似为线性变化，功率密度越高，烧蚀效率越高。以辐照结束时背表面温度及烧蚀质量为目标物理量，对烧蚀过程做了参数敏感性分析，结果表明：热容及热导率对背表面温度的影响较大；树脂含量对烧蚀质量的影响较大，但其相对敏感度随激光功率密度增加而下降；激光功率密度超过1 kW/cm2时，辐射系数对烧蚀质量影响较大，但其相对敏感度随激光功率密度增加而下降。

纤维增强复合材料；激光烧蚀；氧化；内对流；数值模拟

1 引 言

纤维增强复合材料（如K/E、C/E等）具有高比强度、高比模量等特点，在航天航空领域得到了广泛的应用，目前正逐步取代黑色金属、有色金属等传统材料，成为轻质化结构的主要材料，如飞行器的壳体结构（发动机、过渡段、控制舱等）逐步由纤维缠绕复合材料取代钢质材料已成为主流方向。近年来，国外对复合材料的热化学烧蚀进行了大量研究[1-6]，但对于复合材料的激光烧蚀效应研究，因其特殊的应用背景而少有详细报道；国内对复合材料激光烧蚀效应的研究起步较晚，在研究方法上以试验研究为主[7-8]。对复合材料激光烧蚀效应的数值模拟研究以传统的热传导方程为主[9-12]，对复合材料特有的热化学反应等考虑较少。

本文应用自行开发的激光效应模拟程序，详细考虑了纤维增强复合材料激光烧蚀过程中的热解、氧化、升华、辐射及内外对流等效应，建立了纤维增强复合材料激光烧蚀效应的有限元模型，完成了对烧蚀过程的定量描述。以碳纤维/环氧树脂复合材料为例，将模拟结果与实验结果进行了比对，结果符合较好。

2 计算模型

2.1 控制方程

为了得到复合材料烧蚀的控制方程，假设烧蚀过程满足以下条件：

（1）热解气体为理想不可压气体；

（2）气体与固体处于热平衡状态，即材料内部任意位置气体与固体温度相等。

根据上述假设，考虑热解气体流动导致的内对流效果，纤维增强复合材料的激光烧蚀过程可由以下3个控制方程来描述。

质量守恒方程：

热分解方程：

能量守恒方程：

2.2 边界条件

当激光辐照复合材料时，激光能量的吸收存在体吸收及面吸收两种方式。树脂基复合材料受激光辐照时，表面会很快碳化。中红外激光在碳材料表面的趋肤深度很小，可以近似为表面吸收。考虑辐射及外对流，边界条件为：

式中：n代表法向，α为激光耦合系数，I为激光功率密度，ε为材料表面发射率，σ为Stefan-Bo1tz-man常数，T∞为环境温度，h为对流换热系数。对碳纤维/环氧树脂复合材料，激光耦合系数为0.90[4]；受激光辐照时，表面会很快碳化，因此材料表面发射率可取用碳的发射率（0.92）；对流换热系数由经验公式给出：

式中：Nu、Pr、Re分别为努赛尔数、普朗特数、雷诺数，λair为气体热导率，x为沿表面的吹风方向的长度。

综合式（1）～（4）即可模拟烧蚀过程，当某个单元温度超过材料升华温度且单元的焓值超过相变焓时，利用生死单元技术删除该单元。

2.3 氧化过程的处理

材料在空气中会发生氧化反应。高温下，氧化反应速率取决于氧气扩散速率，故可利用传质系数的经验公式来估计氧化反应的速率。

固体材料界面上的气体传质速度为：

式中，α为传质系数（单位为m/s），由经验公式给出：

得到氧气的传质系数以后，氧化对表面烧蚀速率的贡献为：

式中，a为化学计量数。对树脂基复合材料，热解剩余产物为焦碳，高温下碳的氧化产物为一氧化碳，a=2。由于热解气体燃烧会消耗部分氧气，因此上式给出的只是氧化对表面烧蚀速率贡献的上限。本文中氧化的贡献利用表面单元的密度变化来表征，即：

式中，Vce11为表面单元的体积。

2.4 材料参数

复合材料密度、热容及热导率均可用混合率表示为：

式中：下标v表示未分解材料，下标c表示完全分解材料。对于碳纤维/环氧树脂复合材料[2，9]，

沿纤维方向热导率kL=120（W/m·K），厚度方向热导率kh=0.173f+1.452（1-f）（W/m ·K），材料热解参数为A=1.0×1014，n=0.9，E=216.5 kJ/mo1，Q=996 J/g。

2.5 计算方法

由于能量方程由普通的热传导抛物方程变为式（3）的对流扩散方程，标准的有限元方法可能导致解不稳定。本文利用特征有限元方法求解式（3），标准有限元方法求解式（1）。特征有限元的求解步骤可参考文献[13]。

特征有限元方法离散格式为：

2.6 时间步长的选择

由于热化学反应过程是非线性很强的问题，为了得到正确的模拟结果，时间步长不能取得太大。对碳纤维/环氧树脂复合材料，图1给出了不同升温速率β下热化学反应模拟误差随时间步长的变化。

图1 不同升温速率下计算误差随时间步长的变化Fig.1 Variation of re1ative error as time step 1ength at different heating rates

由图1可看出，为保证模拟结果的正确性，时间步长要求小于0.01 s。

3 激光烧蚀模拟结果与讨论

为验证模型，以碳纤维/环氧树脂复合材料为例，模拟了复合材料的激光烧蚀。材料尺寸为5.6 cm×10 cm×2.54 mm，铺设角为[+45，-45，90，0，0，0-45，+45，0，0]。激光光斑直径为25 mm，表面气流速度为0.3 Ma。以材料中心点及背表面中心点为温度采样点，实验结果[9]及数值模拟结果如图2所示。从图中可以看到，本文的模型计算结果与实验结果符合得很好，而不考虑材料内部热解气体对流换热时的预测结果与实验结果具有一定偏差。

图3为功率密度为1810 W/cm2时烧蚀质量随时间的变化，可以看出烧蚀质量随时间的变化近似为线性变化。由于本文利用生死单元方法处理材料烧蚀，单元满足一定条件时将单元删除，因此烧蚀质量的模拟值表现为锯齿状，当单元越小时，曲线越光滑。考虑氧化的结果与不考虑氧化的结果相差极小，氧化的消耗可以忽略不计。实验结果略大于模拟结果，这是由于实际烧蚀过程中会存在力学剥蚀，即在气流及材料内部热应力联合作用下材料表面会有碎屑脱落，含力学剥蚀过程的模拟还有待进一步的研究。

图2 功率密度为1330W/cm2时采样点温度随时间的变化，实验结果来自文献[9]Fig.2 Variation of samp1ing point temperature as time at power density of1330W/cm2，the experiment resu1ts come from ref.[9]

图3 功率密度为1 810W/cm2时烧蚀质量随时间的变化，实验结果来自文献[9]Fig.3 Variation of ab1ation mass as 1oss time at power density of 1 810 W/cm2，the experiment resu1ts come from ref.[9]

图4为辐照时间5 s时烧蚀质量随功率密度的变化。从图中可以看出，烧蚀质量并不随输入能量线性变化，功率密度越高，单位烧蚀质量所需的能量越少。这主要是由于高温下碳材料表面极高的热辐射损失。由于碳的升华温度固定，辐射损失的功率密度也是一定的。激光功率密度越高，材料沉积能量的效率越高。当激光功率密度小于碳升华温度下的辐射功率密度时，材料表面不可能到达升华温度，质量损失只能来自基体热解及表面氧化。

图4 激光辐照5 s时烧蚀质量随功率密度的变化Fig.4 Variation of ab1ation mass 1oss as power density with 1aser irradiation time of 5 s

4 参数敏感性分析

对于激光防护来说，通常关心烧蚀时的材料背表面温度；而对激光烧蚀来说，更关心材料的烧蚀质量。本文取目标物理量分别为背表面温度及总烧蚀质量，利用上文给出的模型，计算了不同功率密度下辐射系数、对流换热系数、热解速率、热导率、树脂含量及热容等参数的相对敏感性。

图5 以烧蚀质量为目标物理量时各参数的相对敏感度随激光功率密度的变化Fig.5 Variation of re1ative sensitivity as 1aser power density when the object is ab1ation mass 1oss

图5给出了以总烧蚀质量为目标物理量时各参数的相对敏感度随激光功率密度的变化。从图中可以看出，树脂含量及辐射系数是较敏感的参数，其相对敏感度随激光功率密度变化较大。功率密度较低时，烧蚀质量主要来源于树脂的分解，因此树脂含量影响较大；功率密度较高时，碳纤维开始升华烧蚀，树脂分解在烧蚀质量中的比例开始下降，故树脂含量的相对敏感度下降。对辐射系数而言，激光功率密度较低时，表面温度较低，热辐射影响很小；随着激光功率密度增加，表面温度升高，热辐射的影响逐渐增大，但表面温度不可能超过碳的升华温度，因此当激光功率密度超过一定值后，碳开始升华，再增加激光功率密度时，辐射系数的相对敏感度又开始下降。因此，如果目标物理量是烧蚀质量，应对树脂含量及辐射系数等参数给予充分关注，尽量准确测量。

图6 以背表面温度为目标物理量时各参数的相对敏感度随激光功率密度的变化Fig.6 Variation of re1ative sensitivity as 1aser power density when the object is rear surface temperature

图6给出了以背表面温度为目标物理量时各参数的相对敏感度随激光功率密度的变化。从图中可以看出，热容及热导率是较敏感的参数，其相对敏感度随激光功率密度变化不大。其它参数的相对敏感度较低，对后表面温度影响不大。因此，如果目标物理量是背表面温度，应对热导率及热容等参数给予充分关注，尽量准确测量。

5 结 论

本文综合考虑了纤维增强复合材料激光烧蚀过程中的热解、氧化、辐射及相变等物理过程，给出了纤维增强复合材料激光烧蚀的物理模型和计算模型。以碳纤维/环氧树脂复合材料为例进行了计算，计算结果与实验结果符合得很好。对于碳纤维/环氧树脂复合材料，考虑内对流时计算结果比不考虑内对流时更准确；较强功率密度激光辐照时，氧化对烧蚀的贡献可以忽略；当激光功率密度一定时，烧蚀率趋于常数。由于存在表面辐射热损失，激光功率密度越高，烧蚀单位质量所需的能量越少；当功率密度小于一定值时，材料中只有基体热解及表面氧化造成的质量损失。通过敏感性分析，针对背表面温度，确定了热容及热导率2个高敏感度参数；对烧蚀质量，确定了树脂含量及辐射系数两个高敏感度参数。

[1] PARKAR Z，MANGUN CH.Ab1ation characteristics of an aromatic thermosetting copo1yester/carbon fiber composite[J]. J.Composite Mater.，2012，46：1819-1830.

[2] ZHANG ZH Y.Thermo-mechanica1behavior of po1ymer composites exposed fire[D].Virginia：Ba1cksburg，2010.

[3] BAHRAMIAN A R，KOKABIM.Ab1ation and therma1degradation behaviour of a composite based on reso1 type pheno1ic resin：processmode1ing and experimenta1[J].Polymer，2006，47：3661-3673.

[4] 陈敏孙，江厚满.切向空气气流对激光烧蚀碳纤维复合材料过程的影响[J].光学精密工程，2011，19（2）：482-486. CHEN M S，JIANG H M.Inf1uence of tangentia1airf1ows on process of 1aser ab1ating carbon-fiber composites[J].Opt. Precision Eng.，2011，19（2）：482-486.（in Chinese）

[5] 徐斌，伍晓宇，罗烽，等.0Cr18Ni9不锈钢箔的飞秒激光烧蚀[J].光学精密工程，2012，20（1）：45-51. XU B，WU X Y，LUO F，etal..Ab1ation of0Cr18Ni9 stain1ess stee1 fi1ms by femtosecond 1aser[J].Opt.Precision Eng.，2012，20（1）：45-51.（in Chinese）

[6] TARASENKO V F，PANCHENKO A N，BULGAKOVA N M，et al..Formation ofmicrostructure on 1iquid meta1 surface under nanosecond 1aser ab1ation[J].Chinese Optics，2011，4（1）：46-52.

[7] 陈林柱.蜂窝复合材料的激光烧蚀效应实验研究[M].西安：西北核技术研究所，2010. CHEN L ZH.Experiment Research of Laser Ablation Effect of Honeycomb Composite Material[M].Xi′an：Northwest Institute of Nuc1ear Techno1ogy，2010.（in Chinese）

[8] 常新龙，李正亮，胡波，等.碳纤维/环氧树脂层合板的激光烧蚀特性分析[J].红外与激光工程，2011，40（9）：1691-1695. CHANG X L，LIZH L，HU B，etal..Ana1ysis of carbon/epoxy 1aminates ab1ation subject to 1aser irradiation[J].Infrared and Laser Eng.，2011，40（9）：1691-1695.（in Chinese）

[9] 汤文辉，冉宪文，徐志宏，等.强激光对靶材烧蚀效应的数值模拟研究[J].航天器环境工程，2010，27（1）：32-34. TANGW H，RAN XW，XU ZH H，etal..Numerica1simu1ation of high intensity 1aser ab1ation ofmateria1s[J].Spacecraft Environment Eng.，2010，27（1）：32-34.（in Chinese）

[10] 李雅娣，张钢锤，吴平，等.炭纤维/环氧树脂复合材料层板连续激光烧蚀数值计算[J].固体火箭技术，2008，31（3）：262-265. LIY D，ZHANGGCH，WU P，etal..Numerica1ca1cu1ation for ab1ation of carbon fiber/epoxy resin 1aminated composites under continuous 1aser irradiation[J].J.Solid Rocket Technology，2008，31（3）：262-265.（in Chinese）

[11] SEMAK V V.Mode1ing of 1aser charring and materia1 remova1in fiberg1assmateria1[J].J.Directed Energy，2006，2：5-21.

[12] GRIFFISC A，MASUMURA R A.Therma1response of graphite epoxy composite subjected to rapid heating[J].J.Composite Materials，1981，15：427-442.

[13] ZIENKIEWICZO C，TAYLOR R L.The Finite ElementMethod：Fluid Fynamics[M].5th ed.Beijing：Tsinghua University Press，2007.

作者简介：

彭国良（1985—），男，湖北大冶人，硕士，助理研究员，2006年于清华大学获得学士学位，2009年于西北核技术研究所获得硕士学位，主要从事强激光效应的理论和数值模拟方面的研究。E-mai1：pg102@163.com

闫 辉（1987—），男，陕西宝鸡人，研究实习员，学士，2010年于南京大学获得学士学位，主要从事激光与物质相互作用的理论与数值模拟方面的研究。E-mai1：yanhui20133@126.com

刘 峰（1962—），男，陕西榆林人，博士，研究员，1983年于兰州大学获得学士学位，1991年于西北核技术研究所获得硕士学位，2006年于西安电子科技大学获得博士学位，主要从事强激光大气传输及目标效应的理论和数值模拟方面的研究。E-mai1：1iufeng-aser@163.com

杜太焦（1972—），男，山西襄汾人，博士，研究员，1994年于南开大学获得学士学位，2002年于中国科技大学获得博士学位，主要从事激光与物质相互作用的理论与数值模拟方面的研究。E-mai1：kewin-y1119@163.com

王玉恒（1977—），男，河南方城人，博士研究生，副研究员，1999年于西北工业大学获得学士学位，2002年于西北核技术研究所获得硕士学位，主要从事激光与物质相互作用的理论与数值模拟方面的研究。E-mai1：intree@126.com

Numerical simulation of laser ablation of fiber-reinforced compositematerials

PENG Guo-1iang1，YAN Hui1，LIU Feng1，2，DU Tai-jiao1，2，WANG Yu-heng1，2
（1.Northwest Institute of Nuclear Technology，Xi′an 710024，China；2.State Key Laboratory of Laser Interaction with Matter，Xi′an 710024，China）
*Corresponding author，E-mail：pgl02@163.com

The paper presents the physics and themath mode1s of 1aser ab1ation for fiber-reinforced composite materia1s in consideration of the physica1processes inc1uding pyrogenation，oxidation，radiation，phase transition，both inner and outer convections.By taking carbon/epoxy composites for examp1es，the 1aser ab1ation process is simu1ated，and the resu1ts agree we11with the experimenta1 resu1ts.The numerica1 resu1ts indicate that the inner convection p1ays an important ro1e in the ab1ation process；the oxidation effect can be ignored during high intense 1aser irradiation；the ab1ation mass is near1y 1inear to 1aser irradiation time when the 1aser power density is a constant and the ab1ation efficiency increaseswith increasing power density.The parameter sensitivity is ana1yzed with the 1ossmass and the rear surface temperature at the end of 1aser irradiation as the object.The ana1ysis resu1ts show that the therma1 capacity and therma1 conductivity havemore inf1uences onthe bear surface temperature；the epoxy vo1ume fraction hasmuch inf1uence on ab1ationmass，but its re1ative sensitivity decreaseswith increasing 1aser power density.Moreover，the radiation coefficient is an important factor when the 1aser power density ismore than 1 kW/cm2，but its re1ative sensitivity decreaseswith increasing 1aser power density.

fiber-reinforced compositemateria1；1aser ab1ation；oxidation；inner convection；numerica1 simu1ation

TG156.99

A

10.3788/CO.20130602.0216

1674-2915（2013）02-0216-07

2012-12-13；

2013-02-13

激光与物质相互作用国家重点实验室基金资助项目（No.SKLLIM1102）