架空电力线路杆塔处雷电感应电压计算方法

摘要

本发明请求保护一种架空电力线路杆塔处雷电感应电压的计算方法。主要内容：1.在考虑杆塔影响情况下的空间雷电电磁场计算。选择采用基于亚网格技术的时域有限差分法，对空气空间进行粗网格剖分，对杆塔进行细网格剖分。2.将雷电电磁场二维计算结果应用于架空电力线路与雷电流回击通道三维空间结构的场量分解。雷电回击电流激发形成的电磁场。3.架空电力线路雷电感应电压计算及杆塔位置雷电感应电压提取。采用Agrawal模型和FDTD法计算架空线路上的雷电感应电压；通过网格查询定位于杆塔位置，提取杆塔处的雷电感应电压。本发明能够有效且更加准确地评估架空电力线路和杆塔处绝缘子、避雷器等所承受的雷电感应电压。

说明书

技术领域

本发明涉及架空电力线路防雷、绝缘配合计算领域，具体为一种架空电力线路杆塔处雷电感应电压计算方法。

背景技术

雷击大地时，雷电释放出的脉冲电流将在数千米的空间范围内产生强烈的辐射电磁场。架空输配电线路距离长，裸露于空气当中，极易与雷电电磁场相耦合，形成线路过电压，严重的会造成线路对地或相间闪络、损坏端接变压器及开关设备等。特别是对于35kV及以下的架空配电线路而言，由于绝缘水平较低，雷电感应过电压已成为线路绝缘闪络、跳闸增多的主要原因。正确评估架空电力线路，特别是位于杆塔处的绝缘子、避雷器等装置所承受的雷电感应电压的强度和时变特性，对于线路防雷、确保其安全稳定运行具有重要意义。

已有文献均在忽略杆塔影响的情况下线路雷电感应过电压的计算，所得计算结果是存在缺陷的。因为根据电磁波传播的基本原理可知，在电磁波传播的路径上出现两种媒质时，由于电磁参数发生突变，会出现反射、折射等现象。因此，当雷电电磁场触及杆塔时，会发生畸变，这样的电磁场耦合到线路上所形成的感应电压也随之发生变化。

若要更准确地计算架空电力线路以及安装在杆塔上的绝缘子、避雷器等所承受的雷电感应电压，从而更合理地进行绝缘配合和防雷设计，就有必要考虑杆塔的影响。

发明内容

针对现有技术中的问题，本发明提供了一种考虑杆塔影响的情况下更加准确地计算架空电力线路所承受的雷电感应电压，为线路导线、绝缘子选型和线路绝缘配合设计提供理论依据的架空电力线路杆塔处雷电感应电压计算方法，本发明的技术方案如下：一种架空电力线路杆塔处雷电感应电压计算方法，其包括以下步骤：

101、获取多导架空线路雷电流波形特征参数，包括电流峰值I₀、波头时间τ₁和波尾时间τ₂、回击速度、衰减常数，建立基底雷电流模型i_s(0,t)和回击通道雷电流i_s(z',t)的数学模型，其中基底雷电流模型i_s(0,t)采用Heidler模型函数：

$>i_s(0,t)=\frac{I_0}{\eta }\frac{{(t/{\tau }_1)}^2}{1+{(t/{\tau }_1)}^2}\exp (-t/{\tau }_2)---\left(1\right)$>

式中：I₀用于限定电流峰值；为峰值电流修正系数；τ₁、τ₂分别表示波头、波尾时间。

回击通道雷电流i_s(z',t)的数学模型采用改进传输线模型：

i(z',t)＝e^-λz'i_s(0,t-z'/v_i) t≥z'/v_i    (2)

式中：z′为雷电通道高度；λ为衰减常数；v_i为雷电流回击速度；

102、获取空气、杆塔的介质参数，具体为在空气中的介电常数ε、磁导率μ和电导率σ；在杆塔中介电常数ε_f、磁导率μ_f和电导率σ_f。设置二维柱坐标系，对计算区域进行亚网格剖分，即对空气区域进行粗网格剖分，对杆塔进行细网格剖分，代入步骤101中的回击通道雷电流模型i_s(z',t)，根据亚网格时域有限差分法即Sub-grid FDTD，迭代计算雷电电磁场；

103、获取多导架空线路空间的位置参数，即在设定的计算坐标系中，获取多导架空导线所处的坐标位置，以及线路与雷击点的垂直距离d，根据架空电力线路与雷电流回击通道三维空间关系，将二维柱坐标系中计算所得的雷电场分量E_r沿架空线轴线方向进行分解转换E_x，此外，根据多导架空线路空间位置参数，计算多导架空线对应的单位长度电感系数矩阵L和单位长度电容系数矩阵C；

104、根据步骤103得到的雷电电场分量沿架空线路轴向的分解计算结果，及单位长度电感系数矩阵L和单位长度电容系数矩阵C，并输入架空线路始端和终端等值阻抗、架空线路长度，计算得出雷电电场对多导架空电力线路的电磁耦合的Agrawal模型，形成含有分布电源的传输线方程组；

105、将架空线路沿x方向以步长Δx剖分，时间离散步长取Δt，可得线路上雷电感应电压的迭代计算式

$>i_k^{n+1}=i_k^n-{\mathrm{\Delta tL}}^{-1}(\frac{u_{k+1}^{\mathrm{sn}}-u_k^{\mathrm{sn}}}{\mathrm{\Delta r}}-E_{\mathrm{xk}}^n)$>

$>u_k^{s(n+1)}=u_k^{\mathrm{sn}}-{\mathrm{\Delta tC}}^{-1}\left(\frac{i_{k+1}^{(n+1)}-i_k^{(n+1)}}{\mathrm{\Delta x}}\right)$>

式中：和分别表示第(n+1)Δt时刻线路长度为kΔx处的感应电流列相量和散射电压向量；

获取架空线路始端和终端等值阻抗，修正边界

$>u_1^{s(n+1)}={(\frac{C}{2\mathrm{\Delta t}}+\frac{1}{{2Z}_S\mathrm{\Delta t}})}^{-1}[(\frac{C}{2\mathrm{\Delta t}}-\frac{1}{{2Z}_S\mathrm{\Delta t}})u_1^{\mathrm{sn}}-\frac{i_1^n}{\mathrm{\Delta t}}+\frac{{\mathrm{hE}}_{\mathrm{zS}}^n}{Z_S\mathrm{\Delta x}}]$>

$>u_{k_{\max }+1}^{s(n+1)}={(\frac{C}{2\mathrm{\Delta t}}+\frac{1}{{2Z}_L\mathrm{\Delta t}})}^{-1}[(\frac{C}{2\mathrm{\Delta t}}-\frac{1}{{2Z}_L\mathrm{\Delta t}})u_{k_{\max }+1}^{\mathrm{sn}}-\frac{i_{k_{\max }}^n}{\mathrm{\Delta x}}+\frac{{\mathrm{hE}}_{\mathrm{zL}}^n}{Z_L\mathrm{\Delta x}}]$>

式中：和分别表示第(n+1)Δt时刻线路长度为kΔx处的感应电流列相量和散射电压向量；Z_S和Z_L分别为线路始端和终端端接阻抗向量；为导线处入射电场的水平分量向量；和分别为线路始端和终端位置入射电场的垂直分量向量；

106、求线路上的雷电感应电压，即散射电压与入射电压之和

$>u_k^{n+1}=u_k^{s(n+1)}-{\int }_0^h{E}_z^{n+1}\mathrm{dz}$>

然后根据架空线路起止位置、档距、杆塔与线路始端距离、网格剖分步长信息，确定杆塔所在位置对应的计算网格数，提取和输出杆塔处的雷电感应电压。

进一步的，步骤102中采用粗网格和细网格对二维雷电电磁场进行划分计算，的具体步骤为：

A、粗网格的空间步长和时间步长选择；Δl_max＜λ_min/10,

式中：Δl_max为网格的最大空间步长；λ_min为雷电电磁波的最小波长；

时间步长的选择依据Courant数值稳定条件

$>\mathrm{\Delta t}<\frac{1}{c_{\max }\sqrt{1/{\left(\mathrm{\Delta x}\right)}^2+1/{\left(\mathrm{\Delta y}\right)}^2+1/{\left(\mathrm{\Delta z}\right)}^2}}$>

B、细网格的空间步长和时间步长选择；

将一个粗网格作n_f等分，则细网格的空间步长为

$>{\mathrm{\Delta l}}_f=\frac{{\mathrm{\Delta l}}_c}{n_f}$>

细网格的时间步长为

$>{\mathrm{\Delta t}}_f=\frac{{\mathrm{\Delta t}}_c}{n_f}$>

C、采用中心差分离散Maxwell方程对粗网格上雷电电磁场分量的迭代计算；

D、基于亚网格FDTD方法的细网格上雷电电磁场分量的迭代计算。

进一步的，步骤C中采用中心差分离散Maxwell方程，得空间内雷电电磁场的FDTD迭代计算式

式中：i，j分别为r方向和z方向上的网格数，n为迭代时间步数；ε_c、μ_c和σ_c分别为空气的介电常数、磁导率和电导率。

而雷电流回击通道上的E_z需要做特别处理：

式中：表示(n+1)Δt时刻、高度为处的回击通道雷电流大小。

进一步的，步骤D中采用亚网格FDTD方法对细网格上雷电电磁场分量的迭代计算具体为：

d1、在杆塔所在位置进行细网格剖分，将原有一个粗网格作n_f等分，E_rf和E_zf表示细网格上的雷电电磁场在方向上的分量；

d2、对细网格从n→n+1时刻的推进分为n_f个小时间步进行，即细网格内部的场分量计算进行n_f次循环迭代计算，迭代计算公式为

式中：ii，jj分别为r方向和z方向上的细网格数，m为细网格内迭代时间步数；ε_f、μ_f和σ_f分别为杆塔介质材料的介电常数、磁导率和电导率。

进一步的，步骤d2中还包括通过插值法对粗细网格重合边界处的场分量的计算步骤，具体为利用n时刻和n+1时刻已知的场分量值以及细网格时间迭代步数进行插值，插值需满足波动方程

$>{▿}^{2}D-\frac{1}{v^2}\frac{{\partial }^{2}D}{{\partial t}^2}=0$>

式中：D表示E_r、E_z或其中之一；v为波速。

本发明的优点及有益效果如下：

本发明提出了基于亚网格技术的时域有限差分法(FDTD)，在考虑杆塔影响的情况下计算架空电力线路上形成的雷电感应电压的理论方法，这样既可以获得较准确的计算结果，又尽可能减小计算量、节省计算时间。

本发明还提出通过雷电场分量分解，将雷电电磁场二维计算结果应用于架空线路和雷击通道三维结构场-线耦合计算的技术方法，达到简化计算、提高效率的目的。上述要点弥补了现有理论研究的空白。

附图说明

图1是本发明优选实施例架空电力线路杆塔处雷电感应电压计算流程图；

图2是雷电电磁场计算粗网格剖分示意图；

图3是雷电电磁场计算细网格剖分示意图；

图4是亚网格FDTD计算程序实现流程图；

图5是将E_r分解为E_x的坐标及网格设置；

图6是雷电电磁场与多导架空线路耦合的Agrawal模型。

具体实施方式

以下结合附图，对本发明作进一步说明：

1.基于FDTD亚网格技术的雷电电磁场计算

说明：鉴于雷电电磁场的瞬变特性，选择时域有限差分法(FDTD)进行计算；进一步地，因为在大尺度空间电磁场的计算过程中需要考虑小尺度的杆塔结构，因此采用了FDTD的亚网格技术。

1.1获取雷电流波形特征参数，包括电流峰值、波头时间和波尾时间、回击速度、衰减常数等，建立基底雷电流模型i_s(0,t)和回击通道雷电流i_s(z',t)的数学模型。

基底雷电流推荐使用Heidler模型函数：

$>i_s(0,t)=\frac{I_0}{\eta }\frac{{(t/{\tau }_1)}^2}{1+{(t/{\tau }_1)}^2}\exp (-t/{\tau }_2)---\left(1\right)$>

式中：I₀用于限定电流峰值；为峰值电流修正系数；τ₁、τ₂用于限定波头、波尾时间。

回击通道雷电流推荐使用改进传输线模型：

i(z',t)＝e^-λz'i_s(0,t-z'/v_i) t≥z'/v_i    (2)

式中：z′为雷电通道高度；λ为衰减常数；v_i为雷电流回击速度。

1.2获取空气、杆塔的介质参数，具体为各自的介电常数ε、磁导率μ和电导率σ。

1.3粗网格的空间步长和时间步长选择

利用FDTD方法进行均匀空间雷电电磁场计算时，网格剖分的尺度在式(3)所示稳定性要求的前提下可适当取大一点，这类网格称为粗网格。网格过小会导致过大计算量，严重影响计算效率。

Δl_max＜λ_min/10    (3)

式中：Δl_max为网格的最大空间步长；λ_min为雷电电磁波的最小波长。

与此同时，时间步长的选择依据Courant数值稳定条件

$>\mathrm{\Delta t}<\frac{1}{c_{\max }\sqrt{1/{\left(\mathrm{\Delta x}\right)}^2+1/{\left(\mathrm{\Delta y}\right)}^2+1/{\left(\mathrm{\Delta z}\right)}^2}}---\left(4\right)$>

式中：c_max为最大波速；Δx、Δy、Δz为网格在x、y、z三个方向上的步长。

雷电电磁场在空气中传播，并且计算时取Δx＝Δy＝Δz＝Δl_c，则时间步长稳定性条件化简为

$>{\mathrm{\Delta t}}_c\le \frac{1}{\sqrt{3}c}{\mathrm{\Delta l}}_c---\left(5\right)$>

式中：c为光速。

1.4细网格的空间步长和时间步长选择

由于杆塔塔身直径相比粗网格的步长Δl_c小得多，如要体现杆塔对雷电电磁场的影响，就需要在杆塔位置处进行精细剖分。

假设将一个粗网格作n_f等分，则细网格的空间步长为

$>{\mathrm{\Delta l}}_f=\frac{{\mathrm{\Delta l}}_c}{n_f}---\left(6\right)$>

细网格的时间步长为

$>{\mathrm{\Delta t}}_f=\frac{{\mathrm{\Delta t}}_c}{n_f}---\left(7\right)$>

1.5粗网格上雷电电磁场分量的迭代计算

雷击大地或建筑物后雷电回击电流形成电磁场，雷电电磁场的空间分布以雷电流通道为轴对称分布，因此可以在二维柱坐标系内进行计算。

满足Δl_max、Δt_c的要求在二维柱坐标系中进行网格剖分，如图2所示，其中E_r、E_z和表示粗网格上的电磁场分量。

电场分量和磁场分量的取样分别为和采用中心差分离散Maxwell方程，可得空间内雷电电磁场的FDTD迭代计算式

式中：i，j分别为r方向和z方向上的网格数，n为迭代时间步数；ε_c、μ_c和σ_c分别为空气的介电常数、磁导率和电导率。

雷电流回击通道上的E_z需要做特别处理：

式中：表示(n+1)Δt时刻、高度为处的回击通道雷电流大小。

1.6细网格上雷电电磁场分量的迭代计算

在杆塔所在位置进行细网格剖分，如图3所示将原有一个粗网格作n_f等分，(以n_f＝3为例)。E_rf和E_zf表示细网格上的场分量。

在FDTD计算中，磁场分量的计算用到相邻4个电场分量，而电场分量的计算用到相邻2个磁场分量，这不仅适用于粗网格，同样适用于细网格。

进一步地，对细网格从n→n+1时刻的推进需分为n_f个小时间步进行，即细网格内部的场分量计算需进行n_f次循环迭代计算，迭代计算公式为

式中：ii，jj分别为r方向和z方向上的细网格数，m为细网格内迭代时间步数；ε_f、μ_f和σ_f分别为杆塔介质材料的介电常数、磁导率和电导率。

直接计算粗细网格重合边界处的场分量需要并不存在的细网格场分量，解决该问题的方法是利用n时刻和n+1时刻已知的场分量值以及细网格时间迭代步数，通过插值法进行求解。插值需满足波动方程

$>{▿}^{2}D-\frac{1}{v^2}\frac{{\partial }^{2}D}{{\partial t}^2}=0---\left(15\right)$>

式中：D表示E_r、E_z或其中之一；v为波速。

基于FDTD亚网格技术的雷电电磁场计算程序实现流程如图4所示。

2.电场分量在架空线路-雷击通道三维结构中的分解计算

(说明：通常情况下，架空电力线路与雷击通道的空间位置是三维结构，为了避免复杂而繁琐的三维雷电电磁场计算，提高计算效率，可将二维柱坐标系内计算的场分量进行分解转换。)

2.1获取多导架空线路空间位置参数，即在设定的计算坐标系中，多导架空导线所处的坐标位置，以及线路与雷击点的垂直距离d等。

2.2设架空线路轴向为x方向，对架空导线进行网格剖分，步长取Δx。

2.3如图5所示为将二维柱坐标系内计算所得E_r分解为E_x时所使用的坐标及网格设置。转换关系式满足

$>r_0=\sqrt{{\left(\mathrm{k\Delta x}\right)}^2+d^2}---\left(12\right)$>

$>E_{\mathrm{xk}}=\frac{E_r\left((i+1)\mathrm{\Delta r}\right)\frac{d}{(i+1)\mathrm{\Delta r}}+E_r\left(\mathrm{i\Delta r}\right)\frac{d}{\mathrm{i\Delta r}}}{2},\mathrm{i\Delta r}<r_0<(i+1)\mathrm{i\Delta r}---\left(13\right)$>

式中：i为r方向网格数；k为x方向网格数。

3.多导架空电力线路及绝缘子雷电感应电压计算

3.1根据多导架空线路空间位置参数，计算多导架空线对应的单位长度电感系数矩阵L和单位长度电容系数矩阵C。其中，C＝P^-1

P为多导架空线路对应电位系数矩阵，为n阶方阵，n为导线数。矩阵中的元素根据镜像法可求得：

$>p_{\mathrm{ii}}=\frac{1}{{2\mathrm{\pi ϵ}}_0}\ln \frac{{2h}_i}{R_i}---\left(14a\right)$>

$>p_{\mathrm{ij}}=\frac{1}{{2\mathrm{\pi ϵ}}_0}\ln \frac{L_{\mathrm{ij}}^{\prime }}{L_{\mathrm{ij}}},p_{\mathrm{ji}}=p_{\mathrm{ij}}(i\ne j)---\left(14b\right)$>

式中：ε₀为真空介电常数；h_i和R_i分别为i导线离地垂直距离和半径；L_ij和L'_ij分别为i导线及其镜像导线到j导线的距离。

架空电力线路近似可视为无损耗线，则

L＝μ₀ε₀P    (15)

式中：μ₀为真空磁导率。

3.2根据计算雷电电场场对多导架空电力线路的电磁耦合的Agrawal模型(如图6所示)，形成含有分布电源的传输线方程组：

$>\left(\begin{array}{c}\frac{\partial }{\partial x}{u}^{\mathrm{sca}}(x,t)+L\frac{\partial }{\partial x}i(x,t)=E_x(x,h,t)\\ \frac{\partial }{\partial x}i(x,t)+C\frac{\partial }{\partial x}{u}^{\mathrm{sca}}(x,t)=0\end{array}\right)---\left(16\right)$>

式中：u^sca(x,t)为散射电压；i(x,t)为线路中的感应电流；E_x(x,h,t)为入射电场在导线处的水平分量。

3.3将架空线路沿x方向以步长Δx剖分，时间离散步长取Δt，以中心差商近似计算式(14)中的偏导，可得线路上雷电感应电压的迭代计算式

$>i_k^{n+1}=i_k^n-{\mathrm{\Delta tL}}^{-1}(\frac{u_{k+1}^{\mathrm{sn}}-u_k^{\mathrm{sn}}}{\mathrm{\Delta r}}-E_{\mathrm{xk}}^n)---\left(17\right)$>

$>u_k^{s(n+1)}=u_k^{\mathrm{sn}}-{\mathrm{\Delta tC}}^{-1}\left(\frac{i_{k+1}^{(n+1)}-i_k^{(n+1)}}{\mathrm{\Delta x}}\right)---\left(18\right)$>

式中：和分别表示第(n+1)Δt时刻线路长度为kΔx处的感应电流列相量和散射电压向量。

获取架空线路始端和终端等值阻抗，修正边界

$>u_1^{s(n+1)}={(\frac{C}{2\mathrm{\Delta t}}+\frac{1}{{2Z}_S\mathrm{\Delta t}})}^{-1}[(\frac{C}{2\mathrm{\Delta t}}-\frac{1}{{2Z}_S\mathrm{\Delta t}})u_1^{\mathrm{sn}}-\frac{i_1^n}{\mathrm{\Delta t}}+\frac{{\mathrm{hE}}_{\mathrm{zS}}^n}{Z_S\mathrm{\Delta x}}]---\left(19\right)$>

$>u_{k_{\max }+1}^{s(n+1)}={(\frac{C}{2\mathrm{\Delta t}}+\frac{1}{{2Z}_L\mathrm{\Delta t}})}^{-1}[(\frac{C}{2\mathrm{\Delta t}}-\frac{1}{{2Z}_L\mathrm{\Delta t}})u_{k_{\max }+1}^{\mathrm{sn}}-\frac{i_{k_{\max }}^n}{\mathrm{\Delta x}}+\frac{{\mathrm{hE}}_{\mathrm{zL}}^n}{Z_L\mathrm{\Delta x}}]---\left(20\right)$>

式中：和分别表示第(n+1)Δt时刻线路长度为kΔx处的感应电流列相量和散射电压向量；Z_S和Z_L分别为线路始端和终端端接阻抗向量；为导线处入射电场的水平分量向量；和分别为线路始端和终端位置入射电场的垂直分量向量。

3.4线路上的雷电感应电压为散射电压与入射电压之和

$>u_k^{n+1}=u_k^{s(n+1)}-{\int }_0^h{E}_z^{n+1}\mathrm{dz}---\left(21\right)$>

3.5根据架空线路起止位置、档距、杆塔与线路始端距离、网格剖分步长等信息，确定杆塔所在位置对应的计算网格数，提取和输出杆塔处的雷电感应电压。

以上这些实施例应理解为仅用于说明本发明而不用于限制本发明的保护范围。在阅读了本发明的记载的内容之后，技术人员可以对本发明作各种改动或修改，这些等效变化和修饰同样落入本发明权利要求所限定的范围。