基于时域差分法的雷击回击通道周围空间电磁场研究

景 弘，王思华

（兰州交通大学自动化与电气工程学院，兰州730070）

基于时域差分法的雷击回击通道周围空间电磁场研究

景 弘，王思华

（兰州交通大学自动化与电气工程学院，兰州730070）

由于雷击回击过程发生了强烈的电荷中和，会在周围空间形成强烈的电磁脉冲，对各种电气设备产生干扰。为此，对回击通道周围的电磁场进行模拟，提出了基于时域有限差分法的电磁场数值计算方法。雷电流模型采用双指数函数与Heidler函数组合的函数模型，回击通道采用改进型的DU模型。将模拟的结果与已有文献的结果进行比较，得出时域有限差分法在模拟雷击回击通道周围电磁场上具有一定有效性。模拟结果表明：雷电回击通道周围空间的电磁场受到回击速度、大地相对介电常数、大地电导率等因素的影响。电场的初始峰值随回击速度的增加而减小，磁场随回击速度的增大变化不大，只是略微有所上升；电场随着大地相对介电常数的增加而增加，波形变越来越陡峭；电场的初始峰值随着电导率的减小而增大，波形上升变的越来越陡峭，这些模拟结果与实测雷电电磁场变化是一致的。

雷电流模型；回击通道；雷电电磁场；回击速度；大地相对介电常数；大地电导率

0 引言

雷电是一种常见的大气放电现象，主要分为云层内部，云层与云层之间，云层与大地间放电。其中，云层与大地间放电对人类生活危害最大。雷电放电时会在周围产生强大的电磁感应效应、高电压波入侵和电磁辐射效应，对附近的建筑物、人员和电子设备构成严重的威胁。因此，对雷电电磁脉冲进行研究显得尤为重要。传统的雷电电磁场理论计算有Quasi-Images算法（也称频域法）和偶极子方法，偶极子方法中，将有限电导大地假定为理想导体，忽略了大地对雷电电磁场的影响，当电磁场在有限电导的地面上传播时，它的高频分量会被有选择地衰减，用偶极子方法得到的电磁场将会与实际的电磁场有较大的差异；频域法考虑了大地电导率的影响，但是其公式中积分存在极点且收敛速度很慢，也很难用数值计算得到精确的结果。在过去10年，FDTD（时域有限差分法）作为麦克斯韦方程解的一种数值计算方法，赢得了极为广泛的应用，它是基于简单的公式迭代，不需要复杂的渐进逼近和格林函数，并且可以通过改变网格的大小来提高精度。

一直以来，人们对雷电电磁场进行大量的研究，文献[1-3]采用偶极子方法对Maxwell方程进行求解，得到不同高度不同水平距离处雷电电磁场的分布规律，而对于雷电电磁场微分方程如何求解并没有给出具体的方法。文献[4-6]分析了几种常用的闪电回击工程模型的有效性，得到MTLL回击模型具有相对更高的精确度和有效性，但是DU模型具有较为明确的物理意义，对DU模型加以改善将得到与实际模型相吻合的理论模型。文献[7-13]采用时域有限差分法分析了雷电电磁场，但对雷电流的选型、雷电模型的选型并没有给出具体的说明，仿真精度有待于提高。

笔者将基于DU模型研究雷电回击通道周围电磁场的分布规律，在求解电磁场微分方程时，采用时域差分法，对比偶极子方法，得到时域有限差分法的有效性，并分析了回击速度，大地相对介电常数，大地电导率对回击通道周围空间电磁场的影响，这对于雷电防护及预防雷电电磁干扰具有指导意义。

1 通道雷电流选择

目前回击通道底部电流通常采用双指数函数和Heidler函数来表示[14-15]，然而双指数函数在初始时刻电流的导数最大，不符合闪电物理机制，Heidler函数不可直接积分，不便于雷电磁场的计算。这里对标准雷电流进行拆分、拟合，得到下式雷电流函数：

式中前半部分表示击穿电流，后半部分表示电晕电流，η为击穿电流峰值的修正因子，τ1、τ2为时间常数，α为波头衰减常数，β为波尾的衰减常数，IBD为击穿电流的峰值，IC为电晕电流的峰值。对于击穿电流和电晕电流相关参数典型的取值为：IBD=28 kA，τ1=0.0612 μs，τ2=15.2312 μs，α=198.5 μs，β=6.4125 μs，IC=20.5 kA。

由该函数得到典型回击通道底部的雷电流及其导数波形如图1所示。

图1 回击通道底部电流及其导数波形Fig.1 The bottom current and its derivative waveform of the return stroke

对比雷电流回击通道中的工程模型[16-17]，发现DU模型具有较为明确的物理意义，其电磁场的计算结果也能反映实测波形的大部分特征。DU模型的电流表达式如下：

式中，τD为电流放电时间常数，v为雷电流的传播速度，一般取（1/3～1/2）c，c为光速。

将（1）式得到的回击通道底部雷电流代入式（2），得到回击通道的雷电流，其波形如图2所示。

图2 DU回击通道电流波形Fig.2 The DU current waveform of the return stroke

2 空间电磁场计算

从Maxwell方程组出发，根据偶极子理论，可以推导出空间任意一点P（r，φ，z）的电磁场：

式中，r，φ，z分别为P点在柱坐标系中的径向坐标，方位角和轴向坐标；R为电流原点到场点的距离，其中R=（（z-z’）2+r2）1/2；ε0，μ0和c分别为真空电导率、磁导率和光速；Er，Ez和Hφ分别为P点水平电场、垂直电场和磁场。（3）和（4）式的第一项为静电场分量，第二项为感应电场分量，第三项为辐射电场分量；（5）式第一项为感应磁场分量，第二项为辐射磁场分量。图3是基于偶极子理论的雷电回击通道电磁场计算模型。

图3 雷电回击通道的天线模型Fig.3 The antenna model of lightning return stroke channel

通过上面3个式子可知，空间任意一点电磁场即有与通道电流和电流导数有关的项，又有与电流对时间的积分有关的项，很难直接积分得到解析公式。因此，必须采用数值方法进行计算。

3 雷电电磁场时域有限差分计算

20世纪60年代K.S.Yee提出了时域有限差分法，并首先运用到解决电磁散射问题。时域有限差分法（FDTD）从Maxwell方程组出发，将时间和空间离散化，在时间和空间上构建网格，用过去时刻的电场和磁场分量来计算下一时间步的电场磁场分量，进而得到整个空间随时间变化的电磁场。图4是基于FDTD的计算模型。

图4 雷电电磁场的计算模型Fig.4 The calculation model of lightning electromagnetic field

对于图4雷电电磁场的计算模型，一般取雷电回击通道高度H=7.5 km，回击速度v=1.3×108，吸收边界条件采用一阶Mur吸收边界条件。在图5采用的二维坐标系下，相应的麦克斯韦方程组为：

图5 二维雷电电磁场的FDTD模型Fig.5 Two dimensional FDTD model of lightning electromagnetic field

将上式的3个方程中心离散可以转换为如下3个方程：

由于计算机内存有限，时域差分问题的空间也是有限的，需要用特殊的边界条件来截断，这里采用一阶Mur吸收边界条件，其对应的表达式如下。

在z轴方向：

在r轴方向：

在轴线上的差分方程：

（1）如果为无源区域：

（2）如果为有源区域：

其中I（0，j+1/2）为在回击通道高度为Δz（j+0.5）处的电流元。

在进行时域有限差分法求解雷电流周围空间的电磁场时，雷电流回击通道采用DU模型，通道基部电流采用双指数函数和Heidler函数结合的函数表示，坐标轴上空间间隔Δr=Δz=0.5 m，地面的电导率σ=2.5×10-4S/m，相对介电常数εr=10，计算距离雷电回击通道r=15 m处的垂直电场，垂直电场的波形图如图6所示。

通过对比图6可以发现，（a）、（b）两图的波形几乎一致，这就说明时域有限差分法在计算雷电空间电磁场的有效性。

为了进一步说明FDTD有限差分法在计算雷电空间电磁场的正确性，这里分别用FDTD有限差分法和偶极子方法计算地面50 m的垂直电场和磁场，计算结果如图7所示。

图6 FDTD方法计算和来自文献[13]实测的r=15 m处的电场Fig.6 Calculation of the electric field at r=15 m by the FDTD methodandTheelectricfieldfromther=15mmeasuredbythe[13]

图7 地面50 m处切向磁场和垂直电场Fig.7 Tangential magnetic field on the ground 50 m and Vertical electric field at the ground 50 m

通过图7可以看出，两种方法的计算结果基本一致，这也进一步说明时域有限差分法在计算雷电周围空间电磁场的正确性。

4 基于FDTD的计算结果

从式（9）-（13）可以看出雷电周围的电磁场与雷电流模型、通道高度、回击速度、大地介电常数、大地电导率有着密切的关系。下面采用FDTD数值方法分析上述因素对雷电周围空间电磁场结果的影响效果。图8（a）、（b）是r=15 m处，回击速度v=1.5×108，1.9×108，2.1×108m/s时，电磁场的变化情况。

图8 不同的回击速度时的电场和磁场Fig.8 Electric field at different speed of return stroke

从图8（a）中可以看出回击速度对电场的影响较大，随着回击速度的增大，电场的初始峰值在减小，陡度在增大，说明电场的高频分量在增大；从图8（b）中可以看出回击速度对磁场的影响较小，随着回击速度的增大，磁场初始峰值有略微的上升，说明磁场的高频分量在下降。这与实际观察雷电电磁场变化是一致的。

为了分析大地相对介电常数，大地电导率对电场的影响效果，这里分别取大地相对介电常数εr=3，5，10、大地电导率σ=1×10-2，1×10-3，2.5×10-4S/m，运用FDTD分析观测点离场点r=15 m处电场的变化情况。图8（c）、（d）是电场的变化情况。

从图8（c）中可以看出相对介电常数对电场有一定的影响，随着相对介电常数的增加，电场的初始峰值在增加，波形上升越来越陡峭，且上升的时间在下降；随着大地电导率的增加，电场的初始峰值在下降，波形变得越来越陡峭，说明电场的高频分量在增加。

5 结论

1）采用双指数函数和Heidler函数组合表示的回击通道基部电流模型能够很好的反映自然条件下雷电机理，方便计算回击通道周围空间的电磁场；采用改进型的DU回击通道模型，能够很好的再现雷电回击测量场的4个特征，得出与实测结果相吻合的模型。

2）与传统的偶极子方法相比，FDTD数值计算方法能够很容易计算回击通道周围空间的电磁场，并且随着空间步长的缩小，计算结果的精度可以进一步的提高。用FDTD计算的回击通道周围的电磁场能够反映回击速度、大地相对介电常数、大地电导率等因数对电磁场的影响效果。

3）随着回击速度的增加，r=15 m处电场初始峰值在下降，电场的高频分量在增加；磁场随回击速度的增加变化不大，初始峰值只是略微上升；随着相对介电常数的增加，电场的初始峰值逐渐增加，波形上升越来越陡峭；随着大地电导率的增加，电场的初始峰值逐渐下降，波形也变得越来越陡峭，这与实测结果基本上是一致的。

[1]曹芙蓉，贺昌辉，杨河林.雷电放电空间的电磁场分析[J].船舶电子工程，2004，24（6）:86-89.CAO Furong，HE Changhui，YANG Helin.Electromag⁃netic field lightning space analysis of[J].Ship Electronic Engineering，2004，24（6）:86-89.

[2]王晓嘉，陈亚洲，万浩江，等.LEMP空间表达式求解及分布规律[J].微波学报，2012，28（6）:16-21.WANG Xiaojia，CHEN Yazhou，WAN Haojiang，et al.LEMP spatial expression and distribution law[J].Journal of microwaves，2012，28（6）:16-21.

[3]刘荣美，汪友华，张岩.雷电回击电磁模型[J].电瓷避雷器，2014，137（2）:81-89.LIU Rongmei，WANG Youhua，ZHANG Yan.Lightning electromagnetic model of[J].Insulators and surge arrest⁃ers，2014，137（02）:81-89.

[4]李韦霖，杨琳，李慧，等.雷电通道模型研究与应用[J].电子技术，2010（5）:174-178.LI Weilin，YANG Lin，LI Hui，et al.Research and appli⁃cation of lightning channel model[J].Electronic technolo⁃gy，2010（5）:174-178.

[5]陈亚洲，王晓嘉，万浩江.闪电回击工程模型的有效分析[J].高电压技术，2012，38（10）:2683-2690.CHEN Yazhou，WANG Xiaojia，WAN Haojiang.Effec⁃tive analysis of lightning engineering model[J].High volt⁃age engineering，2012，38（10）:2683-2690.

[6]MEGUMU Miki，VLADIMIR A.RAKOV，KEITH Rambo，GEORGE H.Schnetzer，and Martin A.Uman.Electric fields near triggered lightning channels measured with Pockels sensors[J].Joural of Geophysical Research，2002，107（16）:1-11.

[7]AKIYOSHI Tatematsu，TAKU Noda.Three-Dimensional FDTD calculation of Lightning Induced voltages on a mul⁃tiphase distribution line with the lightning arresters and an overhead shielding wire[J].IEEE Transactions on Electro⁃magnetic Compatibility，2013，PP（99）:1-9.

[8]SHOZO Sekioka，ICHIRO Matsubara，SHIGERU Yokoya⁃ma.Return-stroke model segmentation and its Application to lightning-induced surges calculation[J].IEEE Transac⁃tions on Electromagnetic Compatibility，2010，53（1）:122-130.

[9]杨栋新.基于时域有限差分法的雷电辐射电磁场的分析研究[D].保定:华北电力大学，2008.YANG Dongxin.Analysis and Research on electromagnet⁃ic field of lightning electromagnetic field based on finite difference time domain method[D].Baoding:North China Electric Power University，2008.

[10]宋庭新.基于网格计算的雷电电磁环境仿真[D].武汉:华中科技大学，2005.SONG Tingxin.Lightning electromagnetic environment simulation based on grid computing[D].Wuhan：Hua⁃zhong University of Science and Technology，2005.

[11]王立.时域有限差分法仿真二维电磁波传播[J].舰船科学技术，2011，33（2）:49-51.WANG Li.Finite difference time domain finite difference method simulation of two-dimensional electromagnetic wave propagation[J].Ship science and technology，2011，33（2）:49-51.

[12]杨华.雷电电磁场FDTD分析[D].青岛:中国海洋大学，2008.YANG Hua.FDTD analysis of lightning electromagnetic field[D].Qingdao:Ocean University of China，2008.

[13]孙建虎，刘颖芳，尹平.雷电流数学模型分析[J].后勤工程学院学报，2012，28（2）:86-91.SUN Jianhu，LIU Yingfang，YIN Ping.Mathematical mod⁃el of lightning current mathematical model analysis[J].Journal of Logistics Engineering College，2012，28（2）:86-91.

[14]陈娜娜.雷电流数学模型的对比[J].电气开关，2010，（3）:82-84.CHEN Nana.Comparison of mathematical model of light⁃ning current[J].Electric switch，2010，（3）:82-84.

[15]赵玉龙，刘光斌，余志勇.雷电流数学模型MATLAB仿真分析[J].电磁仿真，2012，（1）:53-55.ZHAO Yulong，LIU Guangbin，YU Zhiyong.Mathemati⁃cal model of lightning current mathematical model MAT⁃LAB simulation analysis[J].Electromagnetic simulation，2012，（1）:53-55.

[16]魏明.雷电电磁脉冲及其防护[M].北京:国防工业出版社，2010.WEI Ming.Thunder and lightning electromagnetic pulse and its protection[M].Beijing:national defence industry press，2010.

[17]陈则煌，张云峰，周中山，等.加速电荷电磁场求解闪电通道辐射场[J].电瓷避雷器，2015，138（4）:38-42.CHEN Zehuang，ZHANG Yunfeng，ZHOU Zhongshan，accelerated charge of solving the electromagnetic field of the lightning channel radiation field[J].Insulators and surge arresters，2015，138（4）:38-42.

Research of Electromagnetic Field around the Lightning Return-Stroke Channel based on FDTD Method

JING Hong，WANG Sihua
（School of Automatic and Electrical Engineering,Lanzhou Jiaotong University,Lanzhou 730070,China）

The electromagnetic pulse,caused by lightning,will be produced around the return stroke space during the process of Lightning discharge which disturbs electrical equipment.To this end,the elec⁃tromagnetic field around the channel is simulated,and the numerical method of electromagnetic field based on finite difference time domain method(FDTD)is proposed.The lightning current model adopts the function model of the combination of the double exponential function and the Heidler function,and the improved DU model is adopted for the lightning return-stroke channel.The simulation result is com⁃pared with existing literature,and obtained that the effectiveness of the FDTD method is effectiveness on simulating the electromagnetic field around the lightning return-stroke channel.The simulation result shows:The electromagnetic field around the lightning strikes channel is affected by the factors such as the speed of the return strike,the relative permittivity of the earth,the earth's conductivity and other fac⁃tors.The initial peak of the electric field decreases with the increase of the return strike speed,and the magnetic field increases little with the increase of the return strike speed.The electric field increases with the increase of the relative permittivity of the earth,and the wave shape becomes more and more steep.The initial peak of the electric field increases with the decrease of the conductivity,and the waveform will become steeper and steeper.The simulation results are consistent with the measured lightning electromag⁃netic field.

lightning current model；lightning return-stroke channel；lightning electromagnetic field；return stroke speed；relative permittivity of the earth；conductivity of the earth

10.16188/j.isa.1003-8337.2017.05.012

2016-06-02

景弘（1988—），男，硕士，研究方向：高电压与绝缘技术。

国家自然科学基金（编号:51567014）；中国铁路总公司科技研究开发计划项目（编号:2016J010－C）。