Numerical simulation of ablation and carbonization of GFRP irradiated by laser induced microwave transmission decay

玻璃纤维增强树脂基复合材料（简称玻璃钢，英文缩写GFRP）是一种使用广泛的功能复合材料，具有比质量轻、比强度高，易于成形制造等优点。玻璃钢中的增强纤维为玻璃纤维，树脂基体是一种有机高分子聚合物，均具有优良的介电特性，对微波传输的衰减作用很小，因而，玻璃钢常用于制造微波器件的容器或外形蒙皮，如雷达罩等。

当玻璃钢受到激光辐照、气动加热等外界热流载荷作用时，树脂基体受热发生热降解，聚合物分子链被打断、缩聚，形成高度无序的热解碳产物^[1-3]。热解碳是无定形碳，当其温度达到一定程度后，其微观结构产生变化，局部向石墨层状结构转变，形成石墨微晶^[4-5]。当石墨微晶层数和尺寸增加时，石墨层状结构形成的π键电子增加，π键电子在微波作用下形成局域电流，从而将部分微波能量转换为电路的焦耳热，使微波传输产生损耗。在这个过程中，由于无定形碳结构向石墨结构转变存在阈值温度^[6-9]，热解碳的温度演变过程对石墨微晶产生与转化有重要影响。高瑞等^[6]研究了碳化温度对PAN基多孔无定形碳介电性能和吸波性能的影响，发现热处理温度达到800 ℃及以上时，损耗角正切出现了明显增大。V.Zubkova^[10]、M.P.Monsef^[11]、Ge^[12]、Peng^[13]等人分别对碳材料热解、碳化等过程中介电性能变化进行研究，发现当温度小于500 °C时，介电性能没有明显变化，当温度更高时，介电常数和损耗角正切都出现明显变化。V.V.Semak等^[14]根据玻璃钢热解过程，利用三级热解反应方程对激光烧蚀碳产物产生和移除过程进行了建模，并将热解碳体积分数关联到碳颗粒半径进而关联到射频波衰减表达式，但其模型复杂，参数众多，模型中也未考虑碳产物微观结构的影响。

文中针对玻璃钢激光烧蚀碳化后树脂碳微观结构变化导致的微波传输衰减效应开展了数值模拟研究，基于玻璃钢温度场时间演变过程，对影响树脂碳微观结构变化过程的相关量进行分析，提取出微波传输衰减的模型表征量，对数值计算模型进行了验证。

针对玻璃纤维/环氧树脂（E-51环氧树脂）复合材料层合板开展了玻璃钢激光烧蚀碳化致微波传输衰减效应的数值建模研究。

为了模拟玻璃钢激光烧蚀碳化致微波传输衰减效应，将激光辐照-微波传输衰减效应分解为激光辐照耦合、材料热响应、提取表征量、微波传输衰减分析等过程，将激光辐照过程解耦，基于材料温度响应开展分析。

数值建模分析的总体思路如图1所示，首先，设计实验获取材料对激光的耦合特性，作为温度响应计算模型的输入，构建激光辐照玻璃钢温度场的计算模型；其次，通过实验获得特定功率密度激光辐照下材料的表面温度响应，对计算模型进行校验；随后，采用模型计算得到不同激光功率密度参数下试样的温度响应，获得材料的体温度场，实现从实验测量可获得的点/面温度等有限信息向体温度信息的升维处理；根据材料体温度场的时间演变过程，分析提取出模型计算结果中可用于表征微波传输衰减效应的表征量，对该表征量的合理性进行分析。具体方法是，从玻璃钢烧蚀碳化微波传输衰减的实验结果中选取一个数据点（通常选择绝对值较大的数据点），对表征量进行标定，然后根据标定系数对计算获得的表征量变化曲线进行整体校准，分析校准后的表征量变化曲线与微波传输衰减测量结果间的匹配一致性，判断表征量是否具有合理性。如果校核后的表征量与实验结果一致，则说明该表征量可用于模拟微波传输衰减效应。

Figure 1.  Idea scheme of numerical modeling analysis

采用积分球对玻璃钢样品激光反射率、透过率进行测试，得到了材料耦合特性变化规律。实验布局如图2所示，激光经过扩束和光束整形，通过光阑、反射率测量积分球后辐照在试样上，试样后放置积分球测量透过的激光。

Figure 2.  Experimental setup of laser coupling characteristic measurement of GFRP

测量得到试样反射率、透过率随时间变化的曲线后，移除积分球，采用相同实验条件通过热像仪测量试样表面温度，从而得到材料在不同温度下的激光耦合特性曲线，如图3所示。从图3可知，材料对激光的吸收率在459 ℃之前有先下降后上升再下降的变化过程，这和玻璃钢在激光加热作用下的表面状态、树脂整体均匀度变化等有关；在459 ℃后，材料表面热解产生的热解碳产物使材料对激光的吸收由体吸收模式转变为面吸收模式，激光能量主要沉积在材料表面。

Figure 3.  Experiment results of laser coupling characteristics test of GFRP

根据实验现象和材料耦合特性的测试结果，基于商业有限元软件建立了玻璃钢层合板激光辐照温度场有限元计算模型，模型基于傅里叶传热方程，涉及热传导、辐射散热、对流散热过程，激光辐照等效为材料体热源或表面热流。其控制方程组为：

式中：ρ为密度；c为比热容；T为材料温度；T₀为环境温度；k为热导率；$\dddot q$为体生热率；n为表面法向单位矢量；q_conv为对流热流；q_rad为热辐射热流；q_flux为激光等效面热流；h_f为对流换热系数；ε为材料发射率；σ为Stefan-Boltzmann常数，其值为5.67×10⁻⁸ W/(m²·K⁴)；α为材料对激光的吸收率；I₀为入射激光功率密度。

层合板厚度为2 mm，单层玻璃钢厚度为0.1 mm，层合层数为10层。实验试样长宽尺寸为22.8 mm×10.2 mm。激光光斑为均匀分布的矩形光斑，光斑尺寸为16 mm×10.2 mm，即在材料宽度方向为全覆盖。由于玻璃钢层合板热导率较低，激光辐照时间也较短，热量在层合板面内方向的扩散长度较小，故选取材料宽度方向的中截面为特征面，构建了二维计算模型，材料宽度方向任一截面上的温度分布可认为与中截面上的温度分布一致。计算模型根据玻璃钢的层压层数分别构建各层的几何形态，各层间通过接触换热系数（thermal contact conductance coefficient, TCC）定义层间热接触，TCC随材料热解状态变化。模型计算流程如图4所示。在模型中，采用热解温度判据（459 ℃）判断材料对激光的吸收是否出现体吸收模式到面吸收模式的转变。若未出现转变，模型中材料对激光的吸收以体热源形式进行模拟；否则，材料对激光的吸收以面热流形式进行模拟。

Figure 4.  Flow diagram of temperature field calculation model for laser irradiating GFRP laminated plate

模型中，材料的比热容取1044 J/(kg·K)，密度为2 000 kg/m³，热导率1.18 W/(m·K)。

设计实验对计算模型进行了验证。验证实验采用图2所示的实验布局，用功率密度为75 W/cm²的近红外激光辐照在玻璃钢试样上，采用热像仪测量试样表面的最高温度。模型计算结果与实验结果对比如图5所示。从图中可得，模型计算结果与实验结果基本一致。另外，图中实验曲线在300 ℃时截止，是由于热像仪在选择高温档位时，其温度测量区间下限为300 ℃。

Figure 5.  Comparison between calculation results and experiment results of temperature field calculation model for laser irradiating GFRP laminated plate

采用不同功率密度的激光对玻璃钢试样进行烧蚀，获得不同烧蚀碳化程度的测试试样。典型的玻璃钢烧蚀碳化形貌如图6所示。

Figure 6.  Typical ablation and carbonization morphology of GFRP

采用矢量网络分析仪和波导管对试样烧蚀碳化后在10 GHz微波下的介电特性进行测试，得到试样介电常数ε和损耗角正切tanδ，结果如图7所示。

Figure 7.  Experiment results of permittivity ε and tangent of loss angle tanδ of GFRP samples

根据公式（1）可得到试样的微波传输衰减量^[15-16]：

式中：P₀为进入样品前的微波功率；P为经样品传输后的功率；γ为微波的传播系数，其实部与虚部分别为衰减系数α（单位N·p/m，其中1 N·p/m＝8.686 dB/m）与相位系数β（单位rad/m）；l为微波在材料中传输距离；ε_r为介电常数实部；λ₀为微波波长。计算得到的试样微波传输衰减量变化趋势如图8所示。可见，微波传输衰减量变化趋势与损耗角正切tanδ的变化相似。这是由于，在公式（1）中，微波传输衰减量与介电常数实部的方根成正比，与损耗角正切近似成正比，而图7中介电常数实部方根的相对变化率较小（以75 W/cm²处的数据为比较基准，最大相对变化率为15%），损耗角正切tanδ相对变化率较大（最大相对变化率为60%），故微波传输衰减量曲线形态主要与损耗角正切曲线形态相关。

Figure 8.  Tendency of microwave transmission decay varying with laser power density

在环氧树脂热解碳化过程中，温度是影响树脂碳微观结构转变的重要因素。在这个过程中，存在温度阈值T_c，只有材料温度超过该温度阈值时，树脂碳微观结构才发生转变。树脂碳温度越高，微观结构转变速率越快。另外，树脂碳微观结构转变过程是一个与持续时间相关的过程，当树脂碳温度超过T_c后，其微观结构以一定速率从对微波无吸收的无定形碳结构转变为对微波存在吸收的石墨微晶结构。这个过程持续的时间越长，最终转变得到的石墨微晶结构含量越多，对微波传输衰减作用越强。

根据上述分析，为模拟树脂烧蚀碳化所致微波传输衰减变化规律，预设了2个模型表征量S_t,Tc、S_Tt,Tc来提取整合玻璃钢层合板激光辐照温度场计算模型得到的材料体温度场信息，研究这2个表征量是否可以用于描述微波传输衰减变化趋势。其中，S_t,Tc表示网格单元温度超过T_c的持续时间加权和：

式中：N为网格单元总数；t₀为计算总时间；T_i为网格单元i的温度；t_i为第i个网格单元温度超过T_c的总持续时间；n为温度超过T_c的网格单元数。S_Tt,Tc表示网格单元温度超过T_c的温度与持续时间乘积加权和：

式中：T₀为所有计算结果中的最高温度；T_ij为分析时间步j（其时间步长为$\Delta {t_j}$）中网格单元i的温度。N、t₀、T₀可视为表征量中的权重因子。

假定S_t,Tc和S_Tt,Tc与微波传输衰减量A间存在比例关系，即

式中：c_t、c_Tt为比例系数。取一个实验数据点对系数c_t、c_Tt定标，观察A_t、A_Tt与其他实验数据点的匹配程度。若A_t或A_Tt与实验数据A匹配程度高，则说明该表征量具有合理性。

以玻璃钢烧蚀碳化试样介电特性测试实验中的激光参数作为参考，采用玻璃钢层合板激光辐照温度场计算模型进行温度场计算，处理温度场计算结果时在数据处理程序中取不同的温度阈值T_c对温度场的演进过程数据进行处理，得到S_t,Tc、S_Tt,Tc。其中，温度阈值T_c取值时，依据文献研究结果取值，并覆盖较大的区间，以进行进一步的分析。以微波传输衰减实验数据中的最大值为标定点，分别对S_t,Tc、S_Tt,Tc进行标定，可得到标定系数，从而得到模型计算的材料微波传输衰减量变化曲线。典型的模型计算结果与实验结果对比如图9所示。由图可见，温度阈值为873 ℃时，计算结果与实验结果的匹配程度良好，而温度阈值为800、1000 ℃时，计算结果与实验结果在激光功率密度较低时匹配程度较差，产生微波传输衰减作用的功率密度阈值存在减小/增大的现象。分析出现这种现象的原因认为，若模型中温度阈值T_c取值高于树脂碳实际温度阈值T_c0，则当材料温度$T \in [{T_{{\rm{c}}0}},{T_{\rm{c}}})$时，实验中树脂碳已经产生石墨微晶，对微波传输产生衰减作用，但计算模型中提取的S_t,Tc、S_Tt,Tc却因未满足条件未能反映这种作用，表现为计算结果出现微波传输衰减的功率密度阈值增大，如图9(b)所示。反之，若模型中温度阈值T_c取值低于树脂碳实际温度阈值T_c0，则提取的S_t,Tc、S_Tt,Tc提前反映出微波传输衰减作用，表现为计算结果出现微波传输衰减的功率密度阈值减小，如图9(c)所示。另外，图9(b)、(c)中计算结果曲线出现折线的原因尚需进一步分析研究。

Figure 9.  Typical matching results between model calculating results and experiment data

为了量化计算结果与实验结果的匹配一致性，选择决定系数R²为表征量。决定系数R²常用于曲线拟合或回归分析，用于判断拟合或回归模型的优劣，其值越大代表拟合/回归模型越接近实际，其计算式为：

式中：y_i为实测值；${\hat y_i}$为拟合值；$\bar y$为实测值的平均值。对取不同温度阈值时的R²进行比对，得到如图10所示的结果。

Figure 10.  Results matching R² curves of simulation and experiment varying with temperature threshold

可见，计算模型中采用表征量S_t,Tc时R²最大为0.995 6，采用S_Tt,Tc时R²最大为0.9933，故采用表征量S_t,Tc更匹配实验结果。此时，对应的温度阈值为873 ℃。另外，采用表征量S_t,Tc得到的R²最大为0.9956，表明此时计算结果与实验结果间的匹配偏差很小。这说明，玻璃钢激光烧蚀碳化致微波传输衰减效应的数值计算模型是合理有效的。

（1）构建了玻璃钢激光烧蚀碳化致微波传输衰减数值模型分析总体方案，建立了激光辐照玻璃钢温度场计算模型，设计实验获得了材料的温度相关性耦合特性，并对模型进行了校验，结果表明温度场计算模型与实验结果基本一致。

（2）采用不同功率密度激光辐照玻璃钢样品，获得了不同烧蚀碳化程度玻璃钢的介电常数实部、损耗角正切的变化曲线，从而计算得到了微波传输衰减量随激光功率密度的变化规律。

（3）基于玻璃钢温度场演变过程，从数值模型计算结果中提取了网格单元温度超过阈值温度的持续时间加权和、网格单元温度超过阈值温度的持续时间与温度乘积加权和两个表征量进行计算结果与实验数据匹配性分析。结果表明，前者更适合作为烧蚀碳化致微波传输衰减效应的模型表征量，且碳化温度阈值为873 ℃时，可以很好地计算预测玻璃钢烧蚀碳化后的微波传输衰减量。