基于分形理论的雷电先导三维建模与仿真

1、基于分形理论的雷电先导三维建模与仿真魏锥，王之猛 （聊城供电公司，聊城，252000）摘要：雷电是一种自然放电现象，其放电的随机性以及雷电先导在产生和发展过程中易受周围环境影响的特点导致雷电活动对电力系统产生巨大危害。如何描述这种特性在雷电先导发展过程中所起的作用以及规律是特高压输电系统外绝缘优化配置中需要解决的关键问题之一。为此采用分形模型模拟了先导的发展过程，在三维中描述雷电放电通道的分支和扭曲。本课题基于对雷电发展过程准静态场的假设，采用有限差分法作为电场计算的方法，将偏微分方程的问题化为线性方程组问题，并用迭代法进行线性方程的数值解。结合以上计算出的空间电场，采用分形方法完成了三维雷电

2、先导放电通道的建模与仿真，并实现了雷电先导发展过程及相应的空间电势图形静态和动态显示。在此基础上，本课题还讨论了避雷针对闪电轨迹的影响。关键词：分形模型；三维建模；雷电先导；有限差分法；MATLABModeling and Simulation of Lightning Leader in 3-D based on the Fractal TheoryWEI Zhui, WANG Zhi-meng(Liaocheng Power Supply Company, Liaocheng 252000, China)Abstract: Lightning is a natural discharge

3、phenomenon. It causes great damage to power system due to the stochastic nature of lightning discharges and the vulnerability of lightning leaders to various environments in the process of inception and propagation. How to describe the function and regular pattern of this factor in the progressing p

4、rocess of lightning leader is one of the key problems to be solved for the optimal allocation of external insulation in uhv transmission network. Consequently, we introduced the fractal model to simulate the developing process of the lightning leader, reproducing branches and tortuosity of lightning

5、 channels in 3-D simulations. This topic is based on the assumption of quasi-static field in lightning development process, using the finite difference method to calculate the electric field, changing partial differential equations into linear equations, using the iterative method to get the numeric

6、al solution of linear equations. Combining with the above electric field, this topic uses the fractal theory to model and simulate the three-dimensional lightning leader discharge channel, achieving development process of lightning leader, at the same time, the potential changes are shown statically

7、 and dynamically on the graph. And on this basis, the subject also discussed the influence to the lightning trajectory caused by air rod. Keywords: Fractal model; Three-dimensional modeling; Lightning leader; Finite difference method; MATLAB0 引言分形先导发展模型是一种基于分形理论的雷电先导模型，不仅模拟过程与雷电放电过程基本类似，而且结果也能够体现通道的

8、分岔和弯曲。闪电现象属于气体电介质击穿现象，用类似DBM模型何金良，张薛巍，董林，等. 输电线路雷击过程分析的雷电通道分形模型J. 中国科学E 辑: 技术科学， 2009， 39(11): 1818-1823.的方法进行计算机模拟是比较自然的想法，因此一些学者以此为依据来模拟自然界中的闪电现象。中国科学技术大学物理2班的王宁远，用MATLAB实现了有限差分法解静电场边值问题的算法王宁远.有限差分法解静电场的边值问题的算法实现及相关问题讨论.百度文库.，将偏微分方程的问题化为线性方程组问题，并使用了迭代法进行线性方程组的数值解。讨论了从几个角度去优化迭代法的措施。并运用这样的方法解决了闪电的二维

9、静态模拟问题。但作者由于一些原因未进行三维及动态的模拟，因此，将分形模型从二维推广到三维，从静态推广到动态有很大的理论价值和研究意义。 1 模型的建立与简化关于理想模型的建立与文司马文霞，李建标，杨庆，等. 雷电先导分形特性及其在特高压线路耐雷性能分析中的应用J. 高电压技术，2010，36(1): 86-91.基本相同，所以此处不再赘述。文2在二维模拟中对理想模型进行了简化，本文在其基础上研究三维模型的建立与仿真方法。由于三维模型对计算机提出了更高的要求，普通计算机将难以完成100×100×100的网格计算。所以本文提出的三维模型只能降低先导路径的精度要求，采用下列参数模

10、拟。雷电先导从1000米高空开始，计算步长为40米，分级前伸速度300000m/s，先导头部对地电势50MV，经过查阅相关资料，以上假设有一定的科学依据。因此编程时，需开辟了25×25×25空间，故每一单位坐标代表40米，在计算先导打到地面或避雷针的总时间tt时已做了相应变换，不影响最后结果。因为计算电势需要完整的边界条件，所以我们做一个近似的假定，我们认为在整个过程中前后左右点保持最初的线性变化的分布。假设(i,j,k)为通道中的一个点，从该点出发考虑18个可能击穿的点，6个单位坐标距离的点分别为(i-1,j,k)、(i+1,j,k)、(i,j-1,k)、(i,j+1,k

11、)、(i,j,k-1)、(i,j,k-1)，12个坐标距离的点分别为(i-1,j-1,k)、(i-1,j+1,k)、(i+1,j-1,k)、(i+1,j+1,k)、(i-1,j,k-1)、(i-1,j,k+1)、(i+1,j,k-1)、(i+1,j,k+1)、(i,j-1,k-1)、(i,j-1,k+1)、(i,j+1,k-1)、(i,j+1,k+1)。2 程序流程根据以上的理论和假设可以做出程序流程图，如图2-1。显然，该算法包括以下3大模块：模块1：初始情况模块2：计算电势模块3：计算概率，决定下一个击穿点图2-1 程序流程图2.1 模块1 初始情况相关程序代码如下：clc %清除命令窗口

12、中显示的所有输入和输出tic; %启动秒表计时器，输出tic使用后的流逝时间，与toc结合记录程序运行时间lx=25;ly=25;lz=25; %定义三维矩阵的最大行、列、层数 v1=zeros(lz,ly,lx); %三维全零矩阵，用来存放电位值v=v1; %v用来确定边界r=3; %发展概率指数c=3e5; %分级前伸速度T=0; %时间矩阵tt=0; %返回雷电先导历时temp=0;%- 初始电位，从上到下线性增加-for i=1:lz for j=1:ly for k=1:lx v1(i,j,k)=(lz-i)*50/(lz-1); end endendo=sqrt(2);v(1,:,

13、:)=1;v(lz,:,:)=1;v(:,1,:)=1;v(:,ly,:)=1;v(:,:,1)=1;v(:,:,lx)=1; %确定边界b=0; %记录电位的迭代次数i=lz-1; %第一个点的高度，即z坐标j=fix(ly/2); %第一个点的y坐标k=fix(lx/2); %第一个点的x坐标axis(0,25,0,25,0,25) %设置x轴与y轴的限度xlabel('X-axis'),ylabel('Y-axis'),zlabel('Z-axis')%为x,y,z坐标轴命名title('三维模拟得到的先导通道电势示例')

14、%为图形添加标题line(10,10,10,10,0,7,'Color','r','linewidth',3);%画一条红色的线代表避雷针2.2 模块2 计算电势命xi为第i点电势，Laplace方程的离散形式可以写成： (4-1)则对全部n个点，可以列n个方程，同方程组(3-4)。方程组中，第i个方程中xj前的系数aj如此给定：当j=i时，aj=-6；当j为i邻近的六个点之一时，aj=1；其他情况，aj=0。即有 n个方程,每个方程中的n个元素中,除了这个点和与它最邻近的6个点外,都为0元；把给定的m个边界条件(包括击穿点)的电势值代入方程组，

15、剩余的n-m独立方程，恰能给出剩余的n-m个未知网格的电势；这样就把求各点电势问题转化为求解线性方程组的问题了。迭代程序如下：while(maxt>1e-5) b=b+1; maxt=0; for l=2:lz-1 for m=2:ly-1 for n=2:lx-1 if v(l,m,n)=0 v2(l,m,n)=(v1(l,m,n+1)+v1(l,m+1,n)+v1(l+1,m,n)+v1(l,m-1,n)+v1(l,m,n-1)+v1(l-1,m,n)/6; t=abs(v2(l,m,n)-v1(l,m,n); if(t>maxt) maxt=t; end end end en

16、d end v1=v2; end为了加快迭代速度，可以提前使用新值进行迭代，因为在上述程序中，我们的扫描赋值方式是一行一行扫描，在对v2(i，j)赋值时，它周围四个点其中有2个已经被扫描过了，即已经获得了新的数值，这个数值应该更优，所以我们尽量使用新算出的数值进行迭代，因此，只要把上述代码中：v2(l,m,n)=(v1(l,m,n+1)+v1(l,m+1,n)+v1(l+1,m,n)+v1(l,m-1,n)+v1(l,m,n-1)+v1(l-1,m,n)/6;改为：v2(l,m,n)=(v1(l,m,n+1)+v1(l,m+1,n)+v1(l+1,m,n)+v2(l,m-1,n)+v2(l,m

17、,n-1)+v2(l-1,m,n)/6; 就可以将迭代次数减少。2.3 模块3 计算概率，决定下一个击穿点同样是基于分布函数，再由rand(1,1,1)生成一个随机数，这个数必然在0,1区间中随机选取,而且这个数必然对应于唯一一个可能击穿的点，由此即可判定x所在的区间，进而判定pi以及所对应的ai，这样，我们就确定了下一个击穿点，命这个击穿点=0。相关代码与二维类似2，此处不再赘述。3 动画效果和电势的变化3.1 动画效果三维动画采用程序动画，相关代码如下： head=line(x0,k,y0,j,z0,i,'linestyle','-','erasem

18、ode','none','linewidth',2);%连接击穿点与引发此点击穿的点drawnow;对于三维图形，为了更直观地观察图形，可以采用view函数来改变视角 MATLAB 动画制作及举例应用EB/OL.百度文库. 。与z=0平面所成的方向角叫仰角，与x=0所成的夹角叫方位角。这样，默认的三维视角方向仰角为30度，方位角为-37.5度，在MATLAB中，view函数改变所有类型的二维和三维图形的图形视角。view(az,el)和view(az,el)将视角改变到所指定的方位角az和仰角el中。因此view(0,0)可显示x-z轴视图，view(9

19、0,0)可显示y-z轴视图。图3-1为模拟得到的三维雷电先导路径的最终结果，图3-2为x-z轴视图，图3-3为y-z视图。其实程序运行时同样是动态变化的。3.2 电势的变化为了更好地理解雷电先导的特点，我们可以通过以下程序实现在图形上显示电势的动态变化。x,y,z=meshgrid(0:25,0:25,0:25);for ii=1:lz for jj=1:ly for kk=1:lx v3(ii,jj,kk)=v1(kk,jj,ii); end endend slice(v3,10,10,0) colorbar('vert') %显示表征颜色刻度的块，为当前轴添加垂直块 dra

20、wnow; %更新图形窗口head=line(x0,k,y0,j,z0,i,'linestyle','-','erasemode','none','linewidth',2);%连接击穿点与引发此点击穿的点drawnow; 三维电势图，可以采用slice命令实现 MATLAB颜色的使用EB/OL.收稿日期:2013年9月30日作者简介：魏锥（1990），女，本科，主要研究方向为继电保护，740676027王之猛（1986），男，本科，从事继电保护工作该命令显示通过立体图形的矩形切片图。slice(V,sx,sy,sz)显示三元函数V=V(X,Y,Z)确定的超立体形在