CN106354925A - 一种基于电晕放电空间电势分布的模拟方法 - Google Patents

Abstract

本发明是一种基于电晕放电空间电势分布的模拟方法，其仿真模型以泊松方程及粒子连续性方程为基础，增加了能够描述粒子流动性的Navier‑stokes方程。原因在于电晕放电过程是一个粒子移动及碰撞的过程，粒子的快速移动会造成放电空间中气体的流动，对空间电势分布造成影响。利用能够描述气体流动的Navier‑stokes方程与泊松方程和粒子连续性方程进行耦合，通过调整仿真模型的参数得到样品表面的电势分布，再与样品表面电势测量结果进行对比，以验证该模型的有效性。能够快速计算出放电空间的电势分布，包括样品表面的电势分布，为研究电晕放电机理和微观放电过程以及材料极化效果模拟提供了重要的理论支持。

Description

一种基于电晕放电空间电势分布的模拟方法

技术领域

本发明属于空间电势计算领域，具体涉及一种基于电晕放电空间电势分布的模拟方法。

背景技术

电晕放电现象广泛存在于电力系统中。电力设备在制造和运行中的绝缘缺陷在直流高压下会引起电晕放电，同时，高压输电线路的导线附近，雷云静电场下的尖端物体附近，都存在大量的电晕放电现象。电晕放电在施加持续的直流电压时将会产生大量的空间电荷的漂移和积累，主要表现为其迁移区内单一极性离子的大量累积。根据气体放电理论，这种空间电荷的分布与运动方向可能会影响到空气间隙的绝缘性能以及绝缘材料表面的电荷集聚效应。因此，对电晕放电的空间电势分布进行高效、准确的仿真与模拟，将为深入研究气体放电机理提供重要参考，同时对于提高电气系统的绝缘性能及电力工业的安全性有着重要意义。

目前，为更好的了解和应用电晕放电技术，许多学者都在重点研究放电机理和仿真模拟。如Han Yin等人(Modeling of Trichel Pulses in the Negative Corona on aLine-to-Plane Geometry，IEEE Transactions on Magnetics 2014年第50卷)利用泊松方程及粒子连续性方程对针-板电晕放电粒子运动过程进行了模拟。然而电晕放电的过程是粒子快速移动及碰撞的过程，在该过程中会造成气体流动形成电晕风，影响放电过程中粒子的分布，显然利用泊松方程及粒子连续性方程不能够准确地刻画电晕放电过程。

发明内容

本发明是从电晕放电的物理本质出发，通过与Navier-stokes方程的耦合，提供一种科学合理，适用性强，效果佳，能够对电晕放电的空间电势分布进行准确且高效计算的基于电晕放电空间电势分布的模拟方法，其特征是，它包括的步骤有：

1)建立空间电势控制方程

由于电晕放电过程中产生的电子和正负离子改变了实验装置内的空间电势分布，因此实验空间内的电场分布情况由式(1)-式(3)三个方程进行模拟：

$\▿\•\stackrel{\→}{D}=(N_p-N_e-N_n)---\left(1\right)$

$\stackrel{\→}{D}={\mathrm{\ϵ}}_0\stackrel{\→}{E}---\left(2\right)$

$\stackrel{\→}{E}=-\▿\mathrm{\φ}---\left(3\right)$

其中，为静电位移，为电场强度，φ为电势，ε₀为气体介电常数，N_e、N_p、N_n分别为电子、阳离子、阴离子的密度。将式(2)和式(3)代入式(1)中，得到式(4)的泊松方程，即：

${\▿}^2\mathrm{\φ}=-\frac{(N_p-N_e-N_n)}{{\mathrm{\ϵ}}_0}---\left(4\right)$

电介质的表面电势可由电场强度的积分形式表示：

$V={\∫}_0^{d}E\·dx---\left(5\right)$

其中，V为电介质表面电势，d为电介质的厚度；

2)建立粒子连续性方程

为了描述电子和正、负离子的运动过程，引入传统的流体动力学模型，其控制方程由电子连续性方程、正负离子连续性方程和泊松方程组成，电子的密度和正负离子的密度由式(6)-式(8)粒子连续性方程求得：

$\frac{{\∂N}_e}{\∂t}+\▿(N_e{\stackrel{\→}{W}}_e-D_e{\▿N}_e)=(k_i-k_a)N_e+k_d{N}_n-\frac{{\mathrm{\β}}_{\mathrm{ep}}{N}_e{N}_p}{e_0}---\left(6\right)$

$\frac{{\∂N}_p}{\∂t}+\▿\left(N_p{\stackrel{\→}{W}}_p\right)=k_i{N}_e-\frac{{\mathrm{\β}}_{\mathrm{ep}}{N}_e{N}_p}{e_0}---\left(7\right)$

$\frac{{\∂N}_n}{\∂t}+\▿\left(N_n{\stackrel{\→}{W}}_n\right)=k_a{N}_e-k_d{N}_n---\left(8\right)$

其中，分别为电子、阳离子、阴离子的漂移速率，表示为 式中μ_e＝0.05m²/vs、μ_p＝2.24×10^-4m²/vs、μ_n＝2.16×10^-4m²/vs分别为电子、阳离子、阴离子的迁移速率；k_a＝6×10⁶(1/s)、k_d＝0分别为电离系数、附着系数、解离系数，其中EN＝E×10²¹/NO₂，NO₂的空气分子密度约为2.46×10²⁵(1/m³)；β_ep＝5×10^-13m³/s为电子与阳离子的复合系数；为电子扩散系数，其中玻尔兹曼常数k_B＝1.38065×10^-23m²kg·s^-2k^-1，T为温度，e₀＝1.602×10^-19C；

在电晕放电过程中，由于光电离产生的电子所引发二次电子崩影响了空间电荷分布，为了能够准确地描述电晕放电微观粒子的发展过程，在电子连续性控制方程中增加了光电离项，根据实验装置，针极与接地电极近似成相距为3cm的两个单元dV₁和dV₂，dV₁单元辐射的光子被dV₂单元吸收后电离出的电子数由式(9)计算：

$S_{ph}=\frac{1}{4\mathrm{\π}}\·\frac{P_q}{P+P_q}\·{\mathrm{w\μ}}_e{N}_e\∫\frac{f\left(r\right)}{r^2}{\mathrm{dV}}_1{\mathrm{dV}}_2---\left(9\right)$

其中，P和P_q分别为大气压和激发态氮原子的衰减压强，P_q＝3997Pa，w为辐射光子概率；为光子被dV₂单元吸收的概率，f(r)的计算为式(10)：

$f\left(r\right)=\frac{\exp (-k_1{P}_{O_2}r)-\exp (-k_2{P}_{O_2}r)}{rlg(k_2/k_1)}---\left(10\right)$

其中，k₁＝2.63×10^-4cm^-1Pa^-1和k₂＝0.015cm^-1Pa^-1分别为氧气对光波的最小和最大吸收系数，为氧气的压强；

3)建立空气流动方程

在电晕放电过程中粒子会发生碰撞和扩散，在电极附近引起空气流动从而影响粒子分布，为了描述气体的流动性，用Navier-Stokes方程对其进行描述，其推导过程为式(11)和式(12)：

$\frac{\mathrm{Dv}}{\mathrm{Dt}}=-\frac{\▿p}{\mathrm{\ρ}}---\left(11\right)$

是流体梯度算子，得一般形式为：

$\mathrm{\ρ}(\frac{\∂v}{\∂t}+(v\•\▿)v)=-\▿p+f---\left(12\right)$

利用Navier-Stokes方程来描述空气流动对整个放电空间粒子分布的影响，由于所研究的是稳态过程即所以应用的计算为式(13)：

$\mathrm{\ρ}{\stackrel{\→}{W}}_e\•\▿{\stackrel{\→}{W}}_e=-\▿p+\mathrm{\μ}{\▿}^2{\stackrel{\→}{W}}_e-q\▿\mathrm{\φ}---\left(13\right)$

其中，ρ为空气密度，p为空气压力，q为空间电荷密度，为动力粘度；

将泊松方程、粒子连续性方程以及Navier-Stokes方程进行耦合，即得到电晕放电空间电势分布的整体模型为式(14)：

$\left\{\begin{array}{c}{\▿}^2\mathrm{\φ}=-\frac{(N_p-N_e-N_n)}{{\mathrm{\ϵ}}_0}\\ \frac{{\∂N}_e}{\∂t}+\▿(N_e{\stackrel{\→}{W}}_e-D_e{\▿N}_e)=(k_i-k_a)N_e+k_d{N}_n-\frac{{\mathrm{\β}}_{\mathrm{ep}}{N}_e{N}_p}{e_0}+S_{\mathrm{ph}}\\ \mathrm{\ρ}{\stackrel{\→}{W}}_e\•\▿{\stackrel{\→}{W}}_e=-\▿p+\mathrm{\μ}{\▿}^2{\stackrel{\→}{W}}_e-q\▿\mathrm{\φ}\\ \frac{{\∂N}_p}{\∂t}+\▿\left(N_p{\stackrel{\→}{W}}_p\right)=k_i{N}_e-\frac{{\mathrm{\β}}_{\mathrm{ep}}{N}_e{N}_p}{e_0}\\ \frac{{\∂N}_n}{\∂t}+\▿\left(N_n{\stackrel{\→}{W}}_n\right)=k_a{N}_e-k_d{N}_n\end{array}---\left(14\right)\right.$

4)模型求解

采用有限元算法(FEM)对方程组式(14)进行求解，把求解域离散成有限个互不重叠的单元，在每个单元内选择一些合适的节点作为插值点，然后把待求的偏微分方程改写成插值函数组成的线性方程组，对方程组式(14)求解得到所需的解。

对于电晕放电达到稳定状态后的聚合物表面电势衰减特性，电子、正负离子密度的变化率为常数，泊松方程及Navier-Stokes方程直接通过有限元算法进行求解，而粒子连续性方程需要化成偏微分方程的弱解形式进行求解，求解过程为式(16)，

$\left\{\begin{array}{c}{\∫}_{\mathrm{\Ω}}({\mathrm{\μ}}_2{N}_e\▿u\▿v-f_2\▿v-D_{\mathrm{\ϵ}}{\▿N}_{e}v)=0\\ {\∫}_{\mathrm{\Ω}}({\mathrm{\μ}}_1{N}_p\▿u\▿v-f_1\▿v)=0\\ {\∫}_{\mathrm{\Ω}}({\mathrm{\μ}}_3{N}_n\▿u\▿v-f_3\▿v)=0\end{array}\right.---\left(16\right)$

根据电晕放电的实验，需要对仿真模型设置边界条件，其中针极的边界条件为：

$\left\{\begin{array}{c}{W}_e{N}_e-{\stackrel{\→}{D}}_e{\▿N}_e=N_p{W}_p\\ -{\stackrel{\→}{D}}_p{\▿N}_p=0\\ N_n=0\\ V=-10000\left(V\right)\end{array}\right.-$

地电极的边界条件为：

$\left\{\begin{array}{c}-{\stackrel{\→}{D}}_e{\▿N}_e=0\\ N_p=0\\ -{\stackrel{\→}{D}}_n{\▿N}_n=0\\ V=0\left(V\right)\end{array}\right.$

绝缘外壳的边界条件为：

$\left\{\begin{array}{c}-{\stackrel{\→}{D}}_e\▿N_e=0\\ -{\stackrel{\→}{D}}_p\▿N_p=0\\ -{\stackrel{\→}{D}}_n\▿N_n=0\\ {\mathrm{\ϵ}}_0{\mathrm{\ϵ}}_r\stackrel{\→}{E}=0\end{array}\right.$

电介质的边界条件为：

$\left\{\begin{array}{c}-{\stackrel{\→}{D}}_e\▿N_e=0\\ W_p{N}_p-{\stackrel{\→}{D}}_p\▿N_p=0\\ -{\stackrel{\→}{D}}_n\▿N_n=0\\ {\stackrel{\→}{D}}_p-{\stackrel{\→}{D}}_e-{\stackrel{\→}{D}}_n=-{\mathrm{\ϵ}}_r(N_p-N_e-N_n)\end{array}\right.$

利用步骤1)-4)能够实现对电晕放电空间电势分布的模拟，通过实验结果验证模型的有效性，并优化电晕放电参数，提高电晕放电模型空间电势的计算精度。

本发明的一种基于电晕放电空间电势分布的模拟方法从电晕放电的物理本质出发，通过与Navier-stokes方程的耦合，能够对电晕放电的空间电势分布进行准确且高效的计算，有利于研究电晕放电粒子运动的微观过程。具有科学合理，适用性强，效果佳等优点。

附图说明

图1是电晕放电空间电势图；

图2是纯硅橡胶样品表面电势分布曲线图；

图3是掺杂氧化铁硅橡胶样品表面电势分布曲线图；

图4是样品表面电势衰减曲线图。

具体实施方式

下面结合附图和具体实施例对本发明的一种基于电晕放电空间电势分布的模拟方法作进一步描述。

本发明的一种基于电晕放电空间电势分布的模拟方法，包括的步骤有：

1)建立空间电势控制方程

由于电晕放电过程中产生的电子和正负离子改变了实验装置内的空间电势分布，因此实验空间内的电场分布情况由式(1)-式(3)三个方程进行模拟：

$\▿\•\stackrel{\→}{D}=(N_p-N_e-N_n)---\left(1\right)$

$\stackrel{\→}{D}={\mathrm{\ϵ}}_0\stackrel{\→}{E}---\left(2\right)$

$\stackrel{\→}{E}=-\▿\mathrm{\φ}---\left(3\right)$

其中，为静电位移，为电场强度，φ为电势，ε₀为气体介电常数，N_e、N_p、N_n分别为电子、阳离子、阴离子的密度。将式(2)和式(3)代入式(1)中，得到式(4)的泊松方程，即：

${\▿}^2\mathrm{\φ}=-\frac{(N_p-N_e-N_n)}{{\mathrm{\ϵ}}_0}---\left(4\right)$

电介质的表面电势可由电场强度的积分形式表示：

$V={\∫}_0^{d}E\·dx---\left(5\right)$

其中，V为电介质表面电势，d为电介质的厚度；

2)建立粒子连续性方程

为了描述电子和正、负离子的运动过程，引入传统的流体动力学模型，其控制方程由电子连续性方程、正负离子连续性方程(亦称粒子连续性方程)和泊松方程组成，电子的密度和正负离子的密度由式(6)-式(8)粒子连续性方程求得：

$\frac{\∂N_e}{\∂t}+\▿(N_e{\stackrel{\→}{W}}_e-D_e\▿N_e)=(k_i-k_a)N_e+k_d{N}_n-\frac{{\mathrm{\β}}_{\mathrm{ep}}{N}_e{N}_p}{e_0}---\left(6\right)$

$\frac{\∂N_p}{\∂t}+\▿\left(N_p{\stackrel{\→}{W}}_p\right)=k_i{N}_e-\frac{{\mathrm{\β}}_{\mathrm{ep}}{N}_e{N}_p}{e_0}---\left(7\right)$

$\frac{\∂N_n}{\∂t}+\▿\left(N_n{\stackrel{\→}{W}}_n\right)=k_a{N}_e-k_d{N}_n---\left(8\right)$

其中，分别为电子、阳离子、阴离子的漂移速率，表示为 式中μ_e＝0.05m²/vs、μ_p＝2.24×10^-4m²/vs、μ_n＝2.16×10^-4m²/vs分别为电子、阳离子、阴离子的迁移速率；k_a＝6×10⁶(1/s)、k_d＝0分别为电离系数、附着系数、解离系数，其中EN＝E×10²¹/NO₂，NO₂的空气分子密度约为2.46×10²⁵(1/m³)；β_ep＝5×10^-13m³/s为电子与阳离子的复合系数；为电子扩散系数，其中玻尔兹曼常数k_B＝1.38065×10^-23m²kg·s^-2k^-1，T为温度，e₀＝1.602×10^-19C；

在电晕放电过程中由于光电离产生的电子所引发二次电子崩影响了空间电荷分布，为了能够准确地描述电晕放电微观粒子的发展过程，在电子连续性控制方程中增加了光电离项，根据实验装置，针极与接地电极近似成相距为3cm的两个单元dV₁和dV₂，dV₁单元辐射的光子被dV₂单元吸收后电离出的电子数由式(9)计算：

$S_{ph}=\frac{1}{4\mathrm{\π}}\·\frac{P_q}{P+P_q}\·{\mathrm{w\μ}}_e{N}_e\∫\frac{f\left(r\right)}{r^2}{\mathrm{dV}}_1{\mathrm{dV}}_2---\left(9\right)$

其中，P和P_q分别为大气压和激发态氮原子的衰减压强，P_q＝3997Pa，w为辐射光子概率；为光子被dV₂单元吸收的概率，f(r)的计算为式(10)：

$f\left(r\right)=\frac{\exp (-k_1{P}_{O_2}r)-\exp (-k_2{P}_{O_2}r)}{rlg(k_2/k_1)}---\left(10\right)$

其中，k₁＝2.63×10^-4cm^-1Pa^-1和k₂＝0.015cm^-1Pa^-1分别为氧气对光波的最小和最大吸收系数，为氧气的压强；

3)建立空气流动方程

在电晕放电过程中粒子会发生碰撞和扩散，在电极附近引起空气流动从而影响粒子分布，为了描述气体的流动性，用Navier-Stokes方程对其进行描述，其推导过程为式(11)和式(12)：

$\frac{\mathrm{Dv}}{\mathrm{Dt}}=-\frac{\▿p}{\mathrm{\ρ}}---\left(11\right)$

是流体梯度算子，得一般形式为：

$\mathrm{\ρ}(\frac{\∂v}{\∂t}+(v\•\▿)v)=-\▿p+f---\left(12\right)$

利用Navier-Stokes方程来描述空气流动对整个放电空间粒子分布的影响，由于所研究的是稳态过程即所以应用的计算为式(13)：

$\mathrm{\ρ}{\stackrel{\→}{W}}_e\•\▿{\stackrel{\→}{W}}_e=-\▿p+\mathrm{\μ}{\▿}^2{\stackrel{\→}{W}}_e-q\▿\mathrm{\φ}---\left(13\right)$

其中，ρ为空气密度，p为空气压力，q为空间电荷密度，为动力粘度；

将泊松方程、粒子连续性方程以及Navier-Stokes方程进行耦合，即得到电晕放电空间电势分布的整体模型为式(14)：

$\left\{\begin{array}{c}{\▿}^2\mathrm{\φ}=-\frac{(N_p-N_e-N_n)}{{\mathrm{\ϵ}}_0}\\ \frac{{\∂N}_e}{\∂t}+\▿(N_e{\stackrel{\→}{W}}_e-D_e{\▿N}_e)=(k_i-k_a)N_e+k_d{N}_n-\frac{{\mathrm{\β}}_{\mathrm{ep}}{N}_e{N}_p}{e_0}+S_{\mathrm{ph}}\\ \mathrm{\ρ}{\stackrel{\→}{W}}_e\•\▿{\stackrel{\→}{W}}_e=-\▿p+\mathrm{\μ}{\▿}^2{\stackrel{\→}{W}}_e-q\▿\mathrm{\φ}\\ \frac{{\∂N}_p}{\∂t}+\▿\left(N_p{\stackrel{\→}{W}}_p\right)=k_i{N}_e-\frac{{\mathrm{\β}}_{\mathrm{ep}}{N}_e{N}_p}{e_0}\\ \frac{{\∂N}_n}{\∂t}+\▿\left(N_n{\stackrel{\→}{W}}_n\right)=k_a{N}_e-k_d{N}_n\end{array}---\left(14\right)\right.$

4)模型求解

采用有限元算法(FEM)对方程组式(14)进行求解，把求解域离散成有限个互不重叠的单元，在每个单元内选择一些合适的节点作为插值点，然后把待求的偏微分方程改写成插值函数组成的线性方程组，对方程组式(14)求解得到所需的解。

对于电晕放电达到稳定状态后的聚合物表面电势衰减特性，电子、正负离子密度的变化率为常数，泊松方程及Navier-Stokes方程直接通过有限元算法进行求解，而粒子连续性方程需要化成偏微分方程的弱解形式进行求解，求解过程为式(16)，

$\left\{\begin{array}{c}{\∫}_{\mathrm{\Ω}}({\mathrm{\μ}}_2{N}_e\▿u\▿v-f_2\▿v-D_{\mathrm{\ϵ}}{\▿N}_{e}v)=0\\ {\∫}_{\mathrm{\Ω}}({\mathrm{\μ}}_1{N}_p\▿u\▿v-f_1\▿v)=0\\ {\∫}_{\mathrm{\Ω}}({\mathrm{\μ}}_3{N}_n\▿u\▿v-f_3\▿v)=0\end{array}\right.---\left(16\right)$

根据电晕放电的实验，需要对仿真模型设置边界条件，其中针极的边界条件为：

$\left\{\begin{array}{c}{W}_e{N}_e-{\stackrel{\→}{D}}_e{\▿N}_e=N_p{W}_p\\ -{\stackrel{\→}{D}}_p{\▿N}_p=0\\ N_n=0\\ V=-10000\left(V\right)\end{array}\right.-$

地电极的边界条件为：

$\left\{\begin{array}{c}-{\stackrel{\→}{D}}_e{\▿N}_e=0\\ N_p=0\\ -{\stackrel{\→}{D}}_n{\▿N}_n=0\\ V=0\left(V\right)\end{array}\right.$

绝缘外壳的边界条件为：

$\left\{\begin{array}{c}-{\stackrel{\→}{D}}_e\▿N_e=0\\ -{\stackrel{\→}{D}}_p\▿N_p=0\\ -{\stackrel{\→}{D}}_n\▿N_n=0\\ {\mathrm{\ϵ}}_0{\mathrm{\ϵ}}_r\stackrel{\→}{E}=0\end{array}\right.$

电介质的边界条件为：

$\left\{\begin{array}{c}-{\stackrel{\→}{D}}_e\▿N_e=0\\ W_p{N}_p-{\stackrel{\→}{D}}_p\▿N_p=0\\ -{\stackrel{\→}{D}}_n\▿N_n=0\\ {\stackrel{\→}{D}}_p-{\stackrel{\→}{D}}_e-{\stackrel{\→}{D}}_n=-{\mathrm{\ϵ}}_r(N_p-N_e-N_n)\end{array}\right.$

利用步骤1)-4)能够实现对电晕放电空间电势分布的模拟，通过实验结果验证模型的有效性，并优化电晕放电参数，提高电晕放电模型空间电势的计算精度。

针对上述模型采用有限元方法进行空间电势的计算，以整个放电空间作为求解域，并满足上述边界条件即可得到整个放电空间的电势分布，包括电介质材料的表面电势分布，将仿真结果与实验测得的材料表面电势相比，可以验证该模型能够有效模拟电晕放电过程。

以硅橡胶薄膜及掺杂氧化铁的硅橡胶薄膜为试验样品，仿真后整个放电空间电势分布如图1所示，该电势分布图能够清晰地刻画放电空间内粒子的分布情况，可以更好地显示电晕放电的微观过程。而图2和图3分别显示了两种样品表面几何中心所在直线上的电势分布，其峰值为几何中心处的电势，可以看出硅橡胶薄膜及掺杂氧化铁硅橡胶薄膜的仿真结果分别为-1.17kV及-1.19kV。图4则是样品经电晕放电实验处理后表面电势衰减曲线，其衰减速率是描述样品存储电荷性能的指标，可以看出两种材料实验结果的初始值分别为-0.92kV及-0.94kV，实验结果与仿真结果存在较小的差异，原因是在实验过程中将电晕放电处理过的实验样品，移至表面电势衰减测量系统的过程中表面电荷快速衰减，造成实验结果略低于仿真结果。由于误差很小，因此该数值模型能够有效地模拟电晕放电的微观过程。

以上所述仅为本发明的优选实施方式而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种基于电晕放电空间电势分布的模拟方法，其特征是，它包括的步骤有：

1)建立空间电势控制方程

由于电晕放电过程中产生的电子和正负离子改变了实验装置内的空间电势分布，因此实验空间内的电场分布情况由式(1)-式(3)三个方程进行模拟：

$\▿\·\stackrel{\→}{D}=(N_p-N_e-N_n)---\left(1\right)$

$\stackrel{\→}{D}={\mathrm{\ϵ}}_0\stackrel{\→}{E}---\left(2\right)$

$\stackrel{\→}{E}=-\▿\mathrm{\φ}---\left(3\right)$

其中，为静电位移，为电场强度，φ为电势，ε₀为气体介电常数，N_e、N_p、N_n分别为电子、阳离子、阴离子的密度。将式(2)和式(3)代入式(1)中，得到式(4)的泊松方程，即：

${\▿}^2\mathrm{\φ}=-\frac{(N_p-N_e-N_n)}{{\mathrm{\ϵ}}_0}---\left(4\right)$

电介质的表面电势可由电场强度的积分形式表示：

$V={\∫}_0^{d}E\·dx---\left(5\right)$

其中，V为电介质表面电势，d为电介质的厚度；

2)建立粒子连续性方程

为了描述电子和正、负离子的运动过程，引入传统的流体动力学模型，其控制方程由电子连续性方程、正负离子连续性方程和泊松方程组成，电子的密度和正负离子的密度由式(6)-式(8)粒子连续性方程求得：

$\frac{\∂N_e}{\∂t}+\▿(N_e{\stackrel{\→}{W}}_e-D_e\▿N_e)=(k_i-k_a)N_e+k_d{N}_n-\frac{{\mathrm{\β}}_{ep}{N}_e{N}_p}{e_0}---\left(6\right)$

$\frac{\∂N_p}{\∂t}+\▿\left(N_p{\stackrel{\→}{W}}_p\right)=k_i{N}_e-\frac{{\mathrm{\β}}_{ep}{N}_e{N}_p}{e_0}---\left(7\right)$

$\frac{\∂N_n}{\∂t}+\▿\left(N_n{\stackrel{\→}{W}}_n\right)=k_a{N}_e-k_d{N}_n---\left(8\right)$

其中，分别为电子、阳离子、阴离子的漂移速率，表示为 式中μ_e＝0.05m²/vs、μ_p＝2.24×10^-4m²/vs、μ_n＝2.16×10^-4m²/vs分别为电子、阳离子、阴离子的迁移速率；k_a＝6×10⁶(1/s)、k_d＝0分别为电离系数、附着系数、解离系数，其中EN＝E×10²¹/NO₂，NO₂的空气分子密度约为2.46×10²⁵(1/m³)；β_ep＝5×10^-13m³/s为电子与阳离子的复合系数；为电子扩散系数，其中玻尔兹曼常数k_B＝1.38065×10^-23m²kg·s^-2k^-1，T为温度，e₀＝1.602×10^-19C；

在电晕放电过程中，由于光电离产生的电子所引发二次电子崩影响了空间电荷分布，为了能够准确地描述电晕放电微观粒子的发展过程，在电子连续性控制方程中增加了光电离项，根据实验装置，针极与接地电极近似成相距为3cm的两个单元dV₁和dV₂，dV₁单元辐射的光子被dV₂单元吸收后电离出的电子数由式(9)计算：

$S_{ph}=\frac{1}{4\mathrm{\π}}\·\frac{P_q}{P+P_q}\·{\mathrm{w\μ}}_e{N}_e\∫\frac{f\left(r\right)}{r^2}{\mathrm{dV}}_1{\mathrm{dV}}_2---\left(9\right)$

其中，P和P_q分别为大气压和激发态氮原子的衰减压强，P_q＝3997Pa，w为辐射光子概率；为光子被dV₂单元吸收的概率，f(r)的计算为式(10)：

$f\left(r\right)=\frac{\exp (-k_1{P}_{O_2}r)-\exp (-k_2{P}_{O_2}r)}{rlg(k_2/k_1)}---\left(10\right)$

其中，k₁＝2.63×10^-4cm^-1Pa^-1和k₂＝0.015cm^-1Pa^-1分别为氧气对光波的最小和最大吸收系数，为氧气的压强；

3)建立空气流动方程

在电晕放电过程中粒子会发生碰撞和扩散，在电极附近引起空气流动从而影响粒子分布，为了描述气体的流动性，用Navier-Stokes方程对其进行描述，其推导过程为式(11)和式(12)：

$\frac{Dv}{Dt}=-\frac{\▿p}{\mathrm{\ρ}}---\left(11\right)$

是流体梯度算子，得一般形式为：

$\mathrm{\ρ}(\frac{\∂v}{\∂t}+(v\·\▿\left)v\right)=-\▿p+f---\left(12\right)$

利用Navier-Stokes方程来描述空气流动对整个放电空间粒子分布的影响，由于所研究的是稳态过程即所以应用的计算为式(13)：

$\mathrm{\ρ}{\stackrel{\→}{W}}_e\·\▿{\stackrel{\→}{W}}_e=-\▿p+\mathrm{\μ}{\▿}^2{\stackrel{\→}{W}}_e-q\▿\mathrm{\φ}---\left(13\right)$

其中，ρ为空气密度，p为空气压力，q为空间电荷密度，为动力粘度；

将泊松方程、粒子连续性方程以及Navier-Stokes方程进行耦合，即得到电晕放电空间电势分布的整体模型为式(14)：

$\left\{\begin{array}{c}{\▿}^2\mathrm{\φ}=-\frac{(N_p-N_e-N_n)}{{\mathrm{\ϵ}}_0}\\ \frac{\∂N_e}{\∂t}+\▿(N_e{\stackrel{\→}{W}}_e-D_e\▿N_e)=(k_i-k_a)N_e+k_d{N}_n-\frac{{\mathrm{\β}}_{ep}{N}_e{N}_p}{e_0}+S_{ph}\\ \mathrm{\ρ}{\stackrel{\→}{W}}_e\·\▿{\stackrel{\→}{W}}_e=-\▿p+\mathrm{\μ}{\▿}^2{\stackrel{\→}{W}}_e-q\▿\mathrm{\φ}\\ \frac{\∂N_p}{\∂t}+\▿\left(N_p{\stackrel{\→}{W}}_p\right)=k_i{N}_e-\frac{{\mathrm{\β}}_{ep}{N}_e{N}_p}{e_0}\\ \frac{\∂N_n}{\∂t}+\▿\left(N_n{\stackrel{\→}{W}}_n\right)=k_a{N}_e-k_d{N}_n\end{array}\right.---\left(14\right)$

4)模型求解

采用有限元算法(FEM)对方程组式(14)进行求解，把求解域离散成有限个互不重叠的单元，在每个单元内选择一些合适的节点作为插值点，然后把待求的偏微分方程改写成插值函数组成的线性方程组，对方程组式(14)求解得到所需的解。

对于电晕放电达到稳定状态后的聚合物表面电势衰减特性，电子、正负离子密度的变化率为常数，泊松方程及Navier-Stokes方程直接通过有限元算法进行求解，而粒子连续性方程需要化成偏微分方程的弱解形式进行求解，求解过程为式(16)，

$\left\{\begin{array}{c}{\∫}_{\mathrm{\Ω}}({\mathrm{\μ}}_2{N}_e\▿u\▿v-f_2\▿v-D_{\mathrm{\ϵ}}\▿N_{e}v)=0\\ {\∫}_{\mathrm{\Ω}}({\mathrm{\μ}}_1{N}_p\▿u\▿v-f_1\▿v)=0\\ {\∫}_{\mathrm{\Ω}}({\mathrm{\μ}}_3{N}_n\▿u\▿v-f_3\▿v)=0\end{array}\right.---\left(16\right)$

根据电晕放电的实验，需要对仿真模型设置边界条件，其中针极的边界条件为：

$\left\{\begin{array}{c}{W}_e{N}_e-{\stackrel{\→}{D}}_e\▿N_e=N_p{W}_p\\ -{\stackrel{\→}{D}}_p\▿N_p=0\\ N_n=0\\ V=-10000\left(V\right)\end{array}\right.$

地电极的边界条件为：

$\left\{\begin{array}{c}-{\stackrel{\→}{D}}_e\▿N_e=0\\ N_p=0\\ -{\stackrel{\→}{D}}_n\▿N_n=0\\ V=0\left(V\right)\end{array}\right.$

绝缘外壳的边界条件为：

$\left\{\begin{array}{c}-{\stackrel{\→}{D}}_e\▿N_e=0\\ -{\stackrel{\→}{D}}_p\▿N_p=0\\ -{\stackrel{\→}{D}}_n\▿N_n=0\\ {\mathrm{\ϵ}}_0{\mathrm{\ϵ}}_r\stackrel{\→}{E}=0\end{array}\right.$

电介质的边界条件为：

$\left\{\begin{array}{c}-{\stackrel{\→}{D}}_e\▿N_e=0\\ W_p{N}_p-{\stackrel{\→}{D}}_p\▿N_p=0\\ -{\stackrel{\→}{D}}_n\▿N_n=0\\ {\stackrel{\→}{D}}_p-{\stackrel{\→}{D}}_e-{\stackrel{\→}{D}}_n=-{\mathrm{\ϵ}}_r(N_p-N_e-N_n)\end{array}\right.$

利用步骤1)-4)能够实现对电晕放电空间电势分布的模拟，通过实验结果验证模型的有效性，并优化电晕放电参数，提高电晕放电模型空间电势的计算精度。