防雷/防热涂层的雷击热解损伤的跨尺度模拟

董 琪，赵致远，李 婷，贾玉玺，张如良

（1.山东科技大学，青岛，266590；2.航天材料及工艺研究所，北京，100076；3.山东大学，济南，250061）

0 引言

运载火箭、导弹等航天器飞行将面临超高温热流环境，其铝合金或碳纤维复合材料壳体不能直接裸露在外，需涂覆硅橡胶类烧蚀防热涂层［1］，以阻止高温破坏壳体并向内传递。而目前防热涂层通常采用绝缘材料，在雷电环境中飞行会面临雷击破坏，致使航天器在后续飞行中失去防热保护，导致任务失败。因此，开展防热涂层雷击损伤研究十分必要。

雷电作用下防雷/防热涂层产生的直接效应具有显著的电-热-化学-力多物理场耦合特征，涉及涂层材料的绝缘击穿、电阻生热、电弧烧蚀、树脂热解及陶瓷化以及力学破坏等。然而，雷击试验的成本高、危险性高、复杂瞬态过程难以表征，因此，使用计算机进行模拟计算成为目前雷击损伤研究的重要手段。其中，以有限元模拟为代表，能够考虑雷击损伤复杂的电-热-化学-力多物理场耦合特征，量化不同耦合场作用下不同类型的损伤情况［2］，可以有效阐明防雷/防热涂层雷击损伤的机理，为雷击损伤预测提供理论基础。Dong 等［3］建立了碳纤维增强复合材料（Carbon Fiber Reinforced Plastic，CFRP）雷击损伤过程的电-热-化学分析模型并加以模拟分析；Alexandros 等［4］建立了CFRP 层合板雷击损伤过程的电-热耦合模型，对比分析了层合板由螺栓连接和粘贴连接时的雷击损伤情况；Zhu 等［5］则利用CFRP 复合材料雷电损伤有限元模型讨论了其几何形状、纤维排布顺序等因素对损伤程度的影响。

但以有限元为代表的多物理场耦合技术由于时间尺度的限制，在模拟涂层在雷击热效应下发生热解这一过程时尚存在不足，采用的电-热-化学耦合有限元模型计算热解损伤时缺乏可靠的热解反应动力学方程。这是由于热解反应动力学方程的获取通常来源于对材料热重试验曲线的热分析，而热重试验条件远达不到实际雷击所造成的极端升温速率和极端高温情况，所得的动力学方程不能很好地与实际吻合。因此，为更合理地研究防雷/防热涂层雷击损伤机理，进而准确预测其雷击烧蚀损伤，探明涂层材料在雷击极端温度条件下的热解动力学机制及行为具有重要意义。

目前针对此类问题的研究较为少见，然而，针对这一问题，分子模拟手段能够从原子-分子的微观尺度模拟聚合物分子在极端升温速率和极端高温下的热解行为。Zhang等［6］利用分子模拟方法对PVC和巴塞木进行了热解模拟，并分析了模拟结果与雷击损伤的关系；Paajanen等［7］模拟了高温下纤维素的热解并分析了其反应机理与动力学。

本文采用反应分子动力学模拟方法研究甲基苯基硅橡胶在雷电极端升温速率下的热解反应机制，进而通过雷击过程电-热-化学耦合有限元模拟，跨尺度计算甲基苯基硅橡胶基防雷/防热涂层材料体系的雷击热解损伤，并将模拟结果与雷击试验损伤结果作对比，旨在构建能更准确反映涂层雷击实际过程的跨尺度模拟方法，为涂层雷击热解损伤预测以及防热层防雷材料开发提供理论支持。

1 基础理论与模拟方法

1.1 涂层树脂热解反应分子动力学模拟

选取甲基苯基硅橡胶的交联网络状结构作为研究体系，其分子的通用结构式如图1所示。

苯基为硅橡胶提供卓越的耐热、耐低温及耐辐照性能，本文以30%苯基含量构建甲基苯基硅橡胶分子链模型作为交联前的单体模型，如图2所示。

图2 甲基苯基硅橡胶的单体分子模型Fig.2 Molecular model of methyl phenyl silicone rubber monomer

将上述16 个硅橡胶单体分子链装入大小为30.4 Å×30.4 Å×30.4 Å 的具有用周期性边界条件的立方盒子中，构建交联体系，初始建模密度为1 g/cm3。根据氧化脱氢的硫化反应法，使用交联脚本程序［8］代替催化剂的作用以简化交联过程，将硅橡胶链上甲基的C 原子标记为反应位点直接进行交联，形成C—C 键的同时删除甲基上的一个氢原子。模拟交联过程采用COMPASS Ⅱ力场，目标交联度为100%。为增加分子的移动性，设置温度为500 K，最终得到交联度为95.1%的硅橡胶交联结构模型，如图3所示。

图3 甲基苯基硅橡胶交联分子模型Fig.3 Molecular model of cross-linked methyl phenyl silicone rubber

基于上述硅橡胶分子模型，本文通过Lammps 2020 软件，采用Adri van Duin 等［9］开发的ReaxFF 反应力场模拟硅橡胶的热解过程。

在实际的雷击过程中，复合材料雷击作用点处升温速率可达1011K/min以上，极短时间内可达3 000 K以上的极高温度［10］。为使热解反应模拟的条件尽可能接近雷击的实际情况，同时兼顾计算成本和计算机运行效能，本文进行了升温速率位于1012～1014K/s 区间范围内的多次热解反应分子模拟，且将模拟的终止温度设置在7 000 K以上，以保证分子模型完全热解，具体算例的温度设置条件如表1所示。此外，模拟在NPT系综下进行，设置步长为0.25 fs。

表1 不同极端升温速率下硅橡胶热解的分子模拟条件Tab.1 Molecular simulation conditions for the pyrolysis of silicone rubber under different extreme heating rates

根据数据处理需要，本文设置每100步输出一次当前时刻模拟系统内产物信息，该信息包括当前时刻下系统中所有存在的产物的分子式及其相应数量。

1.2 树脂热解反应分子模拟的动力学方程

基于上述硅橡胶热解反应分子模拟输出的不同时刻产物信息，本文根据式（1）计算不同时刻下固相产物所占的质量比：

式中W为质量保留率；x为当前体系中存在的化合物种类总数量；N0为反应开始时体系的微粒数；M0为反应开始时体系的相对分子质量；Nti为t时刻下第i种剩余固相产物的微粒数；Mti为t时刻下第i种剩余固相产物的相对分子质量。

由此可得到硅橡胶热解反应的虚拟热解度曲线。通过不同升温速率下的虚拟热解度曲线，进一步计算动力学参数（活化能Ea、指前因子A、反应级数n）。本文采用的热解动力学参数计算方法如下：

a）Starink法。

Starink［11］对比了热分析动力学方法中经典的Kissinger 方程、Ozawa 方程和Boswell 方程的计算精度，并提出了较高计算精度的Starink方程，其具体形式如下：

式中β为升温速率；T为热力学温度；Rg为普适气体常量；G(α)为反应机理函数的积分式。

根据该方程进行拟合可计算出活化能Ea，且在确定机理函数的前提下，可进一步得到指前因子A。

b）Malek法。

硅橡胶的热解过程是聚合物受热解聚、无规则断链的结构转变，是分子量变小的反应，符合n级化学反应机理，Malek法［12］能够在得到活化能Ea的基础上计算反应级数n。对于n级化学反应机理，Malek法计算反应级数n的具体形式如下：

式中αp为热解速率最大时的热解度；up表示Ea/RgTp；π(up)为Luke近似式［13］中引入的有理函数。

将αp、Tp、Ea、Rg代入式（3）中迭代可求出反应级数n。

在得到活化能Ea、指前因子A、反应级数n后，代入式（4）中，即得到热解反应分子模拟动力学方程。

1.3 防雷/防热涂层雷击热解损伤的有限元模拟

首先，构建与雷击试验中试样一致的涂层平板几何模型，如图4所示。模型由三层材料构成，表层为130 μm的防雷层、次表层为4 mm的防热层、底层为3 mm 的铝基板。为提高计算效率，构建的几何模型为1/4的平板模型，平面尺寸为250 mm×250 mm。

图4 雷击损伤电-热-化学耦合有限元分析1/4涂层平板几何模型Fig.4 Quarter geometric model for electrical-thermal-chemical coupling analysis of lightning damage of the coating

a）有限元边界条件与载荷。

在模型的两个非对称侧面设置零电势，并在模型上表面施加对流和辐射热边界；将雷电弧的电流密度和热流密度定义为面载荷，施加在模型上表面。雷电流密度分布方程定义如下：

式中J为电流密度；r为径向坐标；c为常数；Q为热流密度；R(t)为雷电通道半径，见式（7）。

式中Ip为峰值电流强度；t为时间。

电流载荷遵循人工模拟雷电试验SAE APR5412［14］标准，该标准依据自然雷电现象中的电流变化定义了用于试验研究的电流波形。为使模拟中加载的雷电流波形与试验电流波形保持一致，提取文献［1］中雷击试验所用雷电流波形曲线上测得的数据点进行拟合，见图5。

图5 雷电流波形Fig.5 Lightning current waveform

模型施加电流载荷分两个阶段，第1 阶段施加500 μs 的雷电分量A，第2 阶段不施加雷电流，散热至60 s。其中，第1阶段电流载荷为图5b中拟合得到的红色曲线，其具体函数形式如下：

b）控制方程。

涂层中雷电流传导过程由如下方程控制：

式中V为控制单元面积；S为控制单元表面；rc为单位体积内体积电流源；m表示S的法向量；σ为电导率矩阵；φ为电势。

雷电流在涂层中传导过程中电能会以焦耳热的形式部分转换为热能，相应的温度场控制方程为

式中ρ为密度；cp为比热；k为热导率矩阵；θ为内热源。

其中，内热源又分为焦耳热θe和热解反应热Qp，两者分别由式（11）与式（12）表示。

式中η为能量转换系数；Hs为聚合物热解反应的焓。

式（12）的热解反应热可由1.2 节中热解反应分子模拟动力学方程式（4）中定义的热解速率代入计算得到。

上述涂层树脂热解反应在有限元热-化学反应模块的定义是通过用户自定义子程序HETVAL 和USDFLD实现的。

本文利用子程序USDFLD计算并定义树脂材料的热解反应和热解程度，从而控制雷击过程中材料的电学、热学参数随温度和热分解反应程度的演变；另一方面利用子程序HETVAL定义了材料的升华吸热、热分解吸热等热通量。在仿真计算的每一个增量步开始时，Abaqus求解器会将当前的温度、时间和状态变量等数据传递给子程序USDFLD，随后USDFLD根据这些数据和定义的算法等计算增量步结束时状态变量的值，随后，这些变量的值将被传递给Abaqus 求解器和HETVAL，由HETVAL计算当前增量步所产生的热通量，并将计算结果返回给Abaqus 求解器。Abaqus求解器首先根据USDFLD返回的状态变量确定当前的热导率、电导率、比热容等材料参数，而后再结合HETVAL返回的热通量计算增量步结束时的温度。如此不断迭代，直至仿真结束。

c）材料性能参数。

防热层、防雷层和铝基板的材料属性设置见表2至表4。其中，防雷层的电导率随温度变化的规律如图6所示。

表2 防热层材料性能参数Tab.2 Material performance parameters of the thermal insulation layer

图6 防雷层电导率随温度变化的曲线Fig.6 Curve of electrical conductivity of the lightning protection layer as a function of temperature

表2中防热层材料性能参数取自参考文献［15］。由于防雷层是在与防热层相同的树脂基体中掺杂了金属铜粉制备得到，因此表3中给出的防雷层材料性能是基于表2中硅橡胶的性能参数，根据铜、硅橡胶的比例，按照混合法则计算得到。图6中防雷层电导率随温度变化的数据为试验测试得到。由于超过300 ℃后，防雷层材料的树脂基体将开始发生热分解反应，继续升温过程（高于300 ℃）中的电导率测试难以进行，而在实际仿真中，雷击区域的温度将达到上千摄氏度，电导率的实测数据显然不满足仿真要求。因此，考虑到防雷层中树脂基体的热分解会使其中导电金属粒子接触电阻降低，本文通过唯象的方法，在防雷材料温度超过300 ℃后，将其电导率设置为热解度的线性函数，假设防雷材料完全热解时的电导率为纯铜的电导率（5.714×107S/m）；而在达到铜的沸点（2 562 ℃）时，防雷材料不再导电，电导率退化至1×10-14S/m。

表3 防雷层材料性能参数Tab.3 Material performance parameters of the lightning protective layer

由于雷电流并未击穿防热层到达铝基板，即铝基板并不会参与导电，所以其温度变化不明显，表4中铝基板的材料参数无需考虑温度的影响。

表4 铝基板材料性能参数Tab.4 Material performance parameters of the aluminum substrate

2 结果与分析

2.1 涂层树脂热解反应分子模拟的动力学参数计算

涂层树脂热解过程分子模拟结果可视化处理如图7所示，可以发现，在温度上升的过程中交联甲基苯基硅橡胶分子不断解聚，热分解产生的小分子体现出气体扩散现象，系统体积也随之扩大。

图7 1×1013 K/s升温速率下硅橡胶热解反应分子模拟可视化结果Fig.7 Visualization results of molecular simulation of silicone rubber pyrolysis reaction at a heating rate of 1×1013 K/s

基于硅橡胶热解过程产物数据，根据式（1）计算表1中所列算例的剩余固相产物，进而得到硅橡胶在不同升温速率下的热解度随温度的变化曲线，部分升温速率下所得的虚拟热解度曲线如图8所示。

图8 在不同极端升温速率下硅橡胶热解反应分子模拟所得虚拟热解度曲线Fig.8 Virtual pyrolysis degree curves obtained from molecular pyrolysis simulations of silicone rubber at different extreme heating rates

取图8 中每条热解度曲线上α为0.1、0.2、…、0.9时的数据，依据式（2）的Starink方程，分别计算ln（β/T1.92）和-1.000 8/RgT，绘制散点图并拟合，如图9 所示。图9 中每组拟合直线的斜率为对应不同热解度下的活化能Ea，其数值与拟合相关系数Rc如表5所示。

表5 Starink法所得硅橡胶热解反应分子模拟对应不同热解度下的活化能Tab.5 Activation energies at different pyrolysis degrees obtained from molecular pyrolysis simulation of silicone rubber by Starink method

图9 Starink法所得硅橡胶热解分子模拟对应不同热解度的活化能拟合结果Fig.9 Fitting results of activation energies for different pyrolysis degrees obtained from molecular pyrolysis simulation of silicone rubber by Starink method

根据式（3）的Malek法，将表5计算得到的活化能Ea分别代入不同升温速率下的αp，Tp，Ea，Rg，迭代计算得到反应级数n如表6 所示，取算数平均值，得到的反应级数n为1.041。

表6 Malek法所得硅橡胶热解反应分子模拟的反应级数nTab.6 Reaction order n obtained from the molecular pyrolysis simulation of silicone rubber by Malek method

将活化能Ea、n级化学反应机理函数的积分式G(α)代入式（2）可得到指前因子A。其中，n级化学反应机理函数的积分式G(α)代入反应级数n后形式如下：

将该式与活化能Ea代入式（2），得到指前因子A如表7所示。

表7 Starink法所得硅橡胶热解反应分子模拟对应不同热解度的指前因子ATab.7 Pre-exponential factors A at different pyrolysisi degrees obtained from molecular pyrolysis simulations of silicone rubber by Starink method

将活化能Ea、指前因子A、反应级数n代入式（4），就得到了涂层热解分子模拟的反应动力学方程。

2.2 有限元模拟结果及试验验证

根据上述分子模拟计算得到的反应动力学方程编写用户自定义子程序HETVAL 和USDFLD，在Abaqus软件中模拟防雷/防热涂层的雷击热解损伤。

图10 为雷电流分量A 结束时刻和雷击结束后自然冷却60 s时的温度分布云图。结果表明，雷击结束时，雷击中心处温度接近2 000 ℃，而散热60 s 后导电层与防热层温度均降至30 ℃左右，远低于硅橡胶的热解温度，因此认为雷击结束后60 s时防护涂层的热解损伤过程已经结束。

图10 标准雷电波形A作用下1/4涂层平板模型温度分布Fig.10 Temperature distribution of the quarter flat model under the standard lightning waveform A

图11为标准雷电波形A作用下涂层结构在不同时刻的热解云图。结果显示，标准雷电波形A作用下涂层的热解损伤形貌呈圆形。表面导电防护层热解损伤面积在雷击作用开始的前100 μs 快速扩大，100～500 μs变化速率减慢，而在雷击作用结束后的散热过程中几乎不再扩大。在厚度方向上，热解损伤在雷击作用阶段及其结束后的前1 ms内主要集中在防雷层，此后才开始向防热层扩展。热解损伤深度在雷击结束后约5 ms时达到最大值。

图11 标准雷电波形A作用下1/4涂层平板模型在不同时刻下的热解度云图Fig.11 Pyrolysis degree contours of quarter plat model at different times under the standard lightning waveform A

这是由于防雷层的热解损伤主要由雷电流在材料中传导产生的焦耳热导致，雷击过程的前100 μs雷电流较大，在材料中产生了大量的焦耳热，导致热解区域快速扩张。在100 μs 后，随着雷电流强度的减弱，热解损伤的扩展速率降低。而防热层不导电，其升温热源为防热层与防雷层之间的热传递，传递效率受两者间温差及材料导热效率的影响。而由于防热层具有较好的隔热性能，因此需要其与防雷层温差较大时温度才会显著升高并导致材料的热解损伤，具有一定的滞后性。

为了验证仿真模型与结果的合理性，对比雷击热解损伤有限元模拟结果与硅橡胶防护涂层人工雷击试验的损伤形貌，如图12 所示。由图12 可以发现，试验结果与雷击热解损伤模拟结果均体现出圆形损伤形貌，但试验得到的损伤形貌由靠近中心的完全热解区和外侧的放射状未完全热解区组成。考虑到试验中试样的接地状态，雷电流的传导难以形成未完全热解区的放射状损伤形貌，因此，推测放射状的损伤由雷电电弧扫略并烧蚀形成。而本文所建立的有限元模型仅考虑了雷电流在模型中传导产生的焦耳热及其导致的烧蚀损伤，暂未考虑电弧扫略注入的电弧热等，因此，仅取试验结果中的完全热解区与仿真结果进行对比。试验结果中完全热解区的直径约为220 mm，仿真结果得到的热解区域直径约为195 mm，两者相对误差约为11.4%。以上结果表明，仿真结果与试验结果较为吻合，所建立的有限元仿真模型较为可靠。

图12 防雷/防热涂层雷击损伤模拟与试验结果对比Fig.12 Comparison of lightning damage simulation and experimental results for lightning protection/thermal insulation coatings

3 结束语

针对甲基苯基硅橡胶基防雷/防热涂层，采用反应分子动力学模拟与有限元模拟两种手段，建立了在雷击极端温度环境下涂层热解损伤的跨尺度模拟方法。模拟所得的涂层完全热解形貌及面积与人工模拟雷击试验结果吻合较好。但由于有限元模型未考虑电弧扫略造成的烧蚀作用，仿真结果未能体现试验结果中涂层未完全热解区的放射状形貌。综上，文中构建的跨尺度模拟方法能较可靠地预测硅橡胶基防雷/防热涂层的雷击热解损伤，具有较好的工程应用价值。