瓦斯保护是电力变压器的主保护之一，气体取样分析能有效反映变压器内部是否存在过热或放电故障。传统的以气体含量为主要判据的轻瓦斯监测只发告警信号，通过人工带电检测的方式对瓦斯继电器中气体进行取样分析，但在产生气体体积大于气体继电器整定值前已经发生多次放电气体集聚，从故障预警角度分析该方法有一定滞后性^[1-2]。智能集气盒产品能在轻瓦斯动作后自动取气和分析，但气体成分仍以分析可燃气体为主^[3]，存在放电类型不明的问题。因此，为了及时准确地了解变压器放电情况，要对气体组分的产生机理进行研究，为轻瓦斯报警判据提供理论支撑。

现有研究主要对矿物绝缘油过热或油中放电情况单独进行分析。过热产气对应的热性故障研究相对发展较慢，分析方法主要采用热力学和反应分子动力学。文献[4-5]使用反应分子动力学模拟得出了矿物绝缘油热解的主要反应路径以及小分子气体的生成路径。文献[6-7]利用宏观热力学方法计算出矿物绝缘油热故障下裂解产气过程，总结出矿物绝缘油热裂解的两个阶段。文献[8]提出了基于能量加权的溶解气体分析方法。油中放电对应的电性故障包括局部放电、火花放电和电弧放电，故障严重程度依次升级。目前对于电性故障的研究主要通过3种典型电极放电实验以及分子动力学模拟进行。文献[9]研究了3种典型电极油中放电气体含量变化规律。文献[10]对比了两种电压类型作用下矿物绝缘油产气特性。文献[11]将特征气体的体积分数与放电能量进行关联性分析，通过特征气体的单位焦耳产气量、C₂H₂/H₂比值以及CO的体积分数来判断电性故障的严重程度。文献[12]总结出了矿物绝缘油在相同电应力、不同温度下的裂解产物以及小分子的生成路径。

目前对于电热联合作用下矿物绝缘油裂解产生气体的研究较少，因此本文利用ReaxFF力场模拟电热单独作用下矿物绝缘油的裂解过程，分析了矿物绝缘油的产气特性。在传统的运动方程式中加入电磁理论，分析了热场和电场对矿物绝缘油裂解过程的影响，同时对比了电热联合作用与热场单独作用下特征气体组分比例。最后，通过柱板电极放电实验与仿真结果进行对比，验证仿真结果的合理性，为矿物绝缘油电热作用下的裂解研究提供路径。

矿物绝缘油主要由烷烃（C_nH_2n+2，质量分数约为60%）、环烷烃（C_nH_2n，质量分数为25%~40%）和芳香烃（C_nH_2n+6，质量分数为0~15%）3类化合物混合而成，平均碳原子数为20。矿物绝缘油的裂解过程主要包括各类烃的断链、脱氢、开环和缩合等化学反应。

油浸式变压器发生热性故障或电性故障时，会释放出大量的能量，故障点局部温度迅速升高，附近矿物绝缘油长碳链分子中C—H键和C—C键断裂，同时生成少量的H⁺氢游离基和小分子烃类自由基（如CH₃·，CH₂·和CH·等）。处于高能状态的氢游离基和烃类自由基之间会互相作用，生成氢气（H₂）、甲烷（CH₄）、乙炔（C₂H₂）、乙烯（C₂H₄）、乙烷（C₂H₆）等^[5]。

ReaxFF力场由VAN DUIN A C T等人在2001年提出^[13]。其原理是根据原子间不同的距离产生的势能确定键级，根据键级描述化学体系与成键相关的相互作用，可以平滑地过渡成键和断键之间的体系能量。ReaxFF分子模拟计算的势能曲线与密度泛函理论计算得到的结果吻合较好，且计算时间少了一个量级。相对于量子力学，ReaxFF分子模拟可以完全通过体系驱动确定反应路径。

ReaxFF力场中，势能函数定义如下：

$ \begin{align*} E_{\text {system }}= & E_{\text {bond }}+E_{\text {over }}+E_{\text {under }}+E_{\text {val }}+E_{\text {pen }}+ \\ & E_{\text {tor }}+E_{\text {conj }}+E_{\text {vdWaals }}+E_{\text {Coulomb }} \end{align*} $

(1)

式中：E_system——系统总势能；

E_bond——键能；

E_over，E_under——过配位和欠配位的能量矫正项；

E_val，E_pen——价角能和惩罚能；

E_tor，E_conj——扭转角能和共轭键能；

E_vdWaals，E_Coulomb——范德华能和库仑能。

根据原子间势能计算键级，并判断化学键的形成和断裂，进而模拟化学反应过程。

反应分子动力学模拟过程中运动方程求解有多种方法，其中Velocity Verlet算法^[14]的精度和稳定性较好，但模拟时间较长。为了同时模拟热场和电场对矿物绝缘油裂解过程，需要改进原来的Velocity Verlet算法，在传统的运动方程式中加入电磁理论。反应分子动力学模拟理论框架如图 1所示。

通过加入一个基于牛顿第二定律的扰动项，模拟带电粒子在外加电场下的运动过程。其运动方程可以表示为

$ \boldsymbol{p}_{\mathit{\boldsymbol{i}}}=\boldsymbol{F}_{\boldsymbol{i}}+q_{i} \boldsymbol{E}(\boldsymbol{t})+q_{i} \boldsymbol{v}_{\mathit{\boldsymbol{i}}} \boldsymbol{B}(\boldsymbol{t}) $

(2)

式中：p_i——第i个带电粒子动量；

F_i，v_i——第i个带电粒子受到的力和速度；

q_i——第i个带电粒子的电荷量；

E(t)，B(t)——电场强度、磁场强度的瞬时值。

假定电磁场均匀并在z平面偏振，即电场分量在z轴方向作用，磁场分量在y轴方向作用，电场强度和磁场强度的瞬时值可表示为

$ \boldsymbol{E}(\boldsymbol{t})=E_{\max } \cos (\omega t)(0 \mathbf{i}+0 \mathbf{j}+1 \mathbf{k}) $

(3)

$ \boldsymbol{B}(\boldsymbol{t})=B_{\max } \cos (\omega t)(0 \mathbf{i}+1 \mathbf{j}+0 \mathbf{k}) $

(4)

式中：E_max，B_max——电场强度、磁场强度最大值；

ω——角频率。

q_iE(t)作用方向为z轴，q_iv_iB(t)作用方向为x轴和z轴。进而可以推导出Velocity Verlet算法外加电磁场后改进方程为

$\begin{array}{l} v_{i, x}(t+\Delta t)= \\ \frac{v_{i, x}(t+\Delta t / 2)+\left(\Delta t / 2 m_{i}\right) f_{i, x}(t+\Delta t)}{P}-\frac{(\Delta t / 2) Q\left\{v_{i, z}(t+\Delta t / 2)+\left(\Delta t / 2 m_{i}\right)\left[f_{i, z}(t+\Delta t)+q_{i}(t+\Delta t) E(t+\Delta t)\right]\right\}}{P} \end{array} $

(5)

$ \begin{array}{l} v_{i, z}(t+\Delta t)= \\ \frac{v_{i, z}(t+\Delta t / 2)+\left(\Delta t / 2 m_{i}\right)\left[f_{i, z}(t+\Delta t)+q_{i}(t+\Delta t) E(t+\Delta t)\right]}{P}+\frac{(\Delta t / 2) Q\left[v_{i, x}(t+\Delta t / 2)+\left(\Delta t / 2 m_{i}\right) f_{i, x}(t+\Delta t)\right]}{P} \end{array} $

(6)

其中：

$ P=1+\left[\frac{1}{2} \Delta t q_{i}(t+\Delta t) B(t+\Delta t)\right]^{2} $

(7)

$ Q=q_{i}(t+\Delta t) B(t+\Delta t) $

(8)

式中：v_i，x(t+△t)，v_i，z(t+△t)——t+△t时刻第i个粒子在x轴和z轴方向的速度；

△t——仿真设置的步长；

m_i——第i个粒子的质量；

f_i，x(t+△t)，f_i，z(t+△t)——t+△t时刻第i个粒子在x轴和z轴方向受到的力；

q_i(t+△t)——t+△t时刻第i个粒子的电荷量；

E(t+△t)，B(t+△t)——t+△t时刻电场强度和磁场强度的值。

由于矿物绝缘油体系为非磁性物质，故在分子动力学模拟的数值求解过程中可忽略磁场的作用，在计算过程中可设B(t+△t)为零。

选取矿物绝缘油3种典型的成分，即链烷烃（C₂₀H₄₂）、双环烷烃（C₂₀H₃₈）和双环芳香烃（C₂₀H₂₆）进行模拟，具体建模过程如下。

首先，使用Material Studio软件分别搭建3种烃类的3维模型，并进行结构优化，使分子处于能量最稳定的状态。3种烃类单体模型如图 2所示。

将3种烃类单体建立不定形液态系统，共包含30个分子。依据矿物绝缘油的成分组成比例，烷烃、环烷烃、芳香烃的比例为6∶3∶1，分子个数分别为18、9、3。采用周期性边界条件作为约束条件，以避免不必要的重叠。

然后，对所构建的不定形液态系统进行动力学弛豫，使系统达到稳定状态。计算的截断半径均设为1.25 nm，并选用Smart算法。其具体步骤如下。

（1）退火处理  选择正则[简写为NVT，表示具有确定的粒子数（N）、体积（V）、温度（T）]系综，对模型进行5次循环退火处理。

（2）加压和解压处理  利用等温等压[简写为NPT，表示具有确定的粒子数（N）、压强（P）、温度（T）]系综，对模型进行0.1~10 MPa的加压和解压处理。

（3）能量最小化  采用Forcite模块，对整体模型进行能量最小化。

经过处理后模型内部应力减小，处于稳定状态，液态系统的体系密度约为0.85 g/cm³，接近矿物绝缘油的实际密度，这也说明动力学弛豫的过程是有效的。稳定状态下，液态系统的边长为2.56 nm。

最后，采用Forcite模块的Dynamic功能，系综选项设置为NVT，让系统达到温度为300 K的平衡状态。优化前后的矿物绝缘油分子模型如图 3所示。

由于Material Studio软件不支持电场仿真，需将模型以pdb格式导出，转换成data格式后导入Lammps仿真软件，并进行电热作用仿真。

对矿物绝缘油分子模型进行反应分子动力学模拟，力场参数选择ReaxFF C/H/O^[15]，采用NVT系综，设定模拟时间为100 ps，步长为0.1 ps。根据简谐过渡态理论^[16]和矿物绝缘油裂解原理，温度只会影响反应速率而不会改变反应路径。为减少模拟时间，将温度设定在2 000~3 000 K，并采用PuReMD并行计算方法^[17]。由于电场强度过低不足以导致矿物绝缘油裂解，过高会发生热溢溃现象^[18]，所以电场强度设定在1~10 MV/m。

设定电场强度为零，在温度分别为2 000、2 200、2 400、2 600、2 800 K的情况下，得到的矿物绝缘油分子模型模拟结果如图 4所示。

由图 4可知：2 000 K温度下，矿物绝缘油主要以大分子碳链形式存在，仅产生少量小分子气体；随着温度的升高，大分子碳链数量减少，当温度达到2 800 K时，大分子碳链全部裂解为小分子气体，产生大量H₂、CH₄、C₂H₂、C₂H₄、C₂H₆等特征气体。

统计出不同温度下特征气体分子数如表 1所示。

由表 1可知：H₂、CH₄、C₂H₂分子数均随着温度上升而增加，其中C₂H₂分子数在温度达到2 400 K之前几乎为零，到达2 600 K后快速上升；C₂H₄和C₂H₆分子数随着温度上升先增加后减少，是由于达到一定温度后，气体进一步裂解会导致分子数下降。

将矿物绝缘油分子模型在2 000、2 200、2 400、2 600、2 800 K温度下重复进行多次模拟，结果显示不同温度下的微观裂解过程相似。以C₂₀H₄₂为例，热裂解断链、脱氢反应如图 5所示。

在反应早期阶段主要是碳链的断裂，碳链从靠近中间位置开始断链的比例远远高于从两端开始断链。随着碳链的断裂产生大量自由基，促使脱氢反应加速，生成大量小分子气体，如H₂、CH₄、C₂H₆、C₃H₈。随着温度升高，C₂H₆和C₃H₈会进一步脱氢产生C₂H₄和C₂H₂。环烷烃与芳香烃相较于烷烃更稳定，反应温度更高，在高温下会发生开环反应，热解过程与烷烃类似。此结论与文献[18]中所述热裂解过程吻合。

研究同一温度下不同电场强度对矿物绝缘油裂解的影响，以温度为2 400 K时为例，在电场强度分别为0、2、4、6 MV/m的情况下，得到的矿物绝缘油分子模型模拟结果如图 6所示。

由图 6可知，高强度电场下碳链裂解更加彻底，说明电场强度增大对矿物绝缘油裂解具有促进作用。

不同电场强度下特征气体分子数如表 2所示。

由表 2可知：H₂、CH₄、C₂H₂分子数均随着电场强度增大而增加；C₂H₄和C₂H₆分子数随着电场强度增大先增加后减少，是由于达到一定温度后进一步裂解导致分子数减少；电场强度在0~2 MV/m时，特征气体分子数总量缓慢增加；电场强度超过2 MV/m后，特征气体分子数总量快速增加。

矿物油模型在电场的作用下极化，链内与链间的范德华力降低，使主链的卷曲程度降低。同时也会影响系统整体动能，电场强度增大，系统动能也会随之增大，分子间碰撞加剧，导致系统温度升高，进一步促进矿物绝缘油裂解。

裂解产生的小分子和极性基团加速系统温度上升，在这种正反馈的作用下矿物绝缘油裂解更加剧烈。因此，电场强度过大会导致热溢溃现象。在0~2 MV/m的电场强度下，分子数缓慢增加，矿物绝缘油裂解并不剧烈。当电场强度超过一定值，电场力的作用超过范德华力削弱的极限，模型中分子被极化拉伸以至于断裂，矿物绝缘油裂解反应急剧加速，相对于无外加电场情况，裂解的起始温度更低。

在不同温度和不同电场强度的情况下，电热联合作用与热场单独作用的特征气体分子数总量对比如图 7所示。

由图 7可知：在电场强度小于4 MV/m时，电热联合作用下特征气体分子数总量相对于热场单独作用下少量增加；在电场强度大于4 MV/m时，电场对系统的作用明显，分子数总量快速增加；相同电场强度下，温度越高，分子数总量越多。

以电场强度为4 MV/m为例，分别模拟2 000、2 200、2 400、2 600、2 800 K温度下电热联合作用与热场单独作用的矿物绝缘油裂解过程，裂解过程产生气态小分子C₁-C₃（C₁、C₂、C₃）以及H⁺氢游离基，电热联合作用与热场单独作用的C₁-C₃和H⁺分子数对比如图 8所示。

由图 8可知：电热联合作用下产生的C₁-C₃和H⁺分子数在不同温度下均比热场单独作用下要多，这说明电场对矿物绝缘油的裂解具有一定的促进作用，电场加速矿物绝缘油裂解的作用机理与电场的极化作用和增加系统总动能有关；H⁺分子数在2 200、2 400、2 600 K温度下有大幅增加，这说明电场力的作用使模型中分子被极化拉伸以至于断裂，矿物绝缘油裂解反应急剧加速；C₁-C₃和H⁺分子数在2 800 K温度下小幅增加，这表明系统中大分子碳链已经全部裂解成小分子烃类。

基于国际电工委员会发布的IEC 60243规范，构建大型变压器内部放电故障模拟实验环境，包含高强度的封闭容器（模拟油箱）、电极（极间距可调，用以调整放电能量）、集气装置等，实现对内部放电故障的准确模拟。实验信号测量系统如图 9所示。

高强度的封闭容器中为25^#矿物绝缘油，柱板电极尺寸如图 9（b）所示。实验基本流程如下：首先，预处理矿物绝缘油，去除水分；其次，密封罐体并静置去除气泡减小外界干扰；然后，通过柱板电极进行放电故障模拟，控制电极之间的距离以及施加的电压来模拟故障的严重程度，采集示波器波形数据；接着，放电故障模拟后静置15 min，测定油中气体组分和含量；最后，更换矿物绝缘油并重复进行3次实验，特征参数取3次实验平均值。

实验信号测量系统中，在初始放电电压为23.4 kV时，电极间歇性放电，随着电场强度增大，相同时间内电极间歇性放电的次数增加，同时电弧亮度逐渐升高。柱板电极施加不同初始电压放电时特征气体的体积分数如表 3所示。总烃是CH₄、C₂H₆、C₂H₄ 3种气体的总和。

由表 3可知：当放电电压为23.4 kV时，H₂的体积分数为50.81%，占总气体体积的一半左右，总烃和C₂H₂各占约1/4；随着放电电压的升高，单位时间放电能量升高，H₂和总烃的体积分数逐渐下降，C₂H₂的体积分数不断上升达到37.78%。

通过改进Velocity Verlet算法求解反应分子动力学模拟，得到不同电压强度下特征气体的体积分数如表 4所示。由于反应分子动力学模拟中为了加速反应的速度，温度参数设置范围为2 000~3 000 K，远高于实际电弧温度，因此仅对比实验和模拟结果中体积分数的变化趋势。

由表 4可知：随着电场强度的增大，H₂和总烃的体积分数逐渐下降，分别降到46.76% 和17.62%，C₂H₂的体积分数不断上升到35.62%；C₂H₂与H₂的体积分数比值随着电场强度的增大从46.08%上升到76.18%。

将表 3中柱板电极放电实验的放电电压除以电极间距3 mm得到的电场强度分别为7.8、8.0、8.2、8.4 MV/m，将对应的特征气体的体积分数与表 4中的结果对比可知，在实验和模拟结果中，H₂、总烃和C₂H₂的体积分数变化趋势一致，随着电场强度的增大，H₂和总烃的体积分数逐渐下降，C₂H₂的体积分数不断上升。

本文从微观角度出发，构建了矿物绝缘油分子模型，通过反应分子动力学模拟，分析了热场和电场对矿物绝缘油裂解过程的影响，对比了电热联合作用与热场单独作用的差异，并通过柱板电极放电实验验证了仿真结果的合理性，得到结论如下。

（1）热场是矿物绝缘油裂解的主要影响因素，同时电场对矿物油裂解也起到促进作用。电场通过极化作用和增加系统总动能来加速矿物绝缘油的裂解。相同电场强度下，温度上升时，H⁺分子数急剧变化的原因是电场力作用使模型中分子被极化拉伸以至于断裂。

（2）随着电场强度的增大，H₂的体积分数逐渐下降，C₂H₂的体积分数不断上升，产生特征气体的分子数总量显著增加，但对起始放电电场的影响需要进一步研究。

（3）反应分子动力学模拟结果与实验结果中特征气体的体积分数变化规律一致。这说明通过反应分子动力学研究矿物绝缘油电热联合作用裂解是一种有效的途径。