基于分形理論的雷電先導(dǎo)三維建模與仿真

1、基于分形理論的雷電先導(dǎo)三維建模與仿真魏錐，王之猛 （聊城供電公司，聊城，252000）摘要：雷電是一種自然放電現(xiàn)象，其放電的隨機(jī)性以及雷電先導(dǎo)在產(chǎn)生和發(fā)展過程中易受周圍環(huán)境影響的特點導(dǎo)致雷電活動對電力系統(tǒng)產(chǎn)生巨大危害。如何描述這種特性在雷電先導(dǎo)發(fā)展過程中所起的作用以及規(guī)律是特高壓輸電系統(tǒng)外絕緣優(yōu)化配置中需要解決的關(guān)鍵問題之一。為此采用分形模型模擬了先導(dǎo)的發(fā)展過程，在三維中描述雷電放電通道的分支和扭曲。本課題基于對雷電發(fā)展過程準(zhǔn)靜態(tài)場的假設(shè)，采用有限差分法作為電場計算的方法，將偏微分方程的問題化為線性方程組問題，并用迭代法進(jìn)行線性方程的數(shù)值解。結(jié)合以上計算出的空間電場，采用分形方法完成了三維雷電

2、先導(dǎo)放電通道的建模與仿真，并實現(xiàn)了雷電先導(dǎo)發(fā)展過程及相應(yīng)的空間電勢圖形靜態(tài)和動態(tài)顯示。在此基礎(chǔ)上，本課題還討論了避雷針對閃電軌跡的影響。關(guān)鍵詞：分形模型；三維建模；雷電先導(dǎo)；有限差分法；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 引言分形先導(dǎo)發(fā)展模型是一種基于分形理論的雷電先導(dǎo)模型，不僅模擬過程與雷電放電過程基本類似，而且結(jié)果也能夠體現(xiàn)通道的

8、分岔和彎曲。閃電現(xiàn)象屬于氣體電介質(zhì)擊穿現(xiàn)象，用類似DBM模型何金良，張薛巍，董林，等. 輸電線路雷擊過程分析的雷電通道分形模型J. 中國科學(xué)E 輯: 技術(shù)科學(xué)， 2009， 39(11): 1818-1823.的方法進(jìn)行計算機(jī)模擬是比較自然的想法，因此一些學(xué)者以此為依據(jù)來模擬自然界中的閃電現(xiàn)象。中國科學(xué)技術(shù)大學(xué)物理2班的王寧遠(yuǎn)，用MATLAB實現(xiàn)了有限差分法解靜電場邊值問題的算法王寧遠(yuǎn).有限差分法解靜電場的邊值問題的算法實現(xiàn)及相關(guān)問題討論.百度文庫.，將偏微分方程的問題化為線性方程組問題，并使用了迭代法進(jìn)行線性方程組的數(shù)值解。討論了從幾個角度去優(yōu)化迭代法的措施。并運(yùn)用這樣的方法解決了閃電的二維

9、靜態(tài)模擬問題。但作者由于一些原因未進(jìn)行三維及動態(tài)的模擬，因此，將分形模型從二維推廣到三維，從靜態(tài)推廣到動態(tài)有很大的理論價值和研究意義。 1 模型的建立與簡化關(guān)于理想模型的建立與文司馬文霞，李建標(biāo)，楊慶，等. 雷電先導(dǎo)分形特性及其在特高壓線路耐雷性能分析中的應(yīng)用J. 高電壓技術(shù)，2010，36(1): 86-91.基本相同，所以此處不再贅述。文2在二維模擬中對理想模型進(jìn)行了簡化，本文在其基礎(chǔ)上研究三維模型的建立與仿真方法。由于三維模型對計算機(jī)提出了更高的要求，普通計算機(jī)將難以完成100100100的網(wǎng)格計算。所以本文提出的三維模型只能降低先導(dǎo)路徑的精度要求，采用下列參數(shù)模擬。雷電先導(dǎo)從1000米

10、高空開始，計算步長為40米，分級前伸速度300000m/s，先導(dǎo)頭部對地電勢50MV，經(jīng)過查閱相關(guān)資料，以上假設(shè)有一定的科學(xué)依據(jù)。因此編程時，需開辟了252525空間，故每一單位坐標(biāo)代表40米，在計算先導(dǎo)打到地面或避雷針的總時間tt時已做了相應(yīng)變換，不影響最后結(jié)果。因為計算電勢需要完整的邊界條件，所以我們做一個近似的假定，我們認(rèn)為在整個過程中前后左右點保持最初的線性變化的分布。假設(shè)(i,j,k)為通道中的一個點，從該點出發(fā)考慮18個可能擊穿的點，6個單位坐標(biāo)距離的點分別為(i-1,j,k)、(i+1,j,k)、(i,j-1,k)、(i,j+1,k)、(i,j,k-1)、(i,j,k-1)，12

11、個坐標(biāo)距離的點分別為(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 程序流程根據(jù)以上的理論和假設(shè)可以做出程序流程圖，如圖2-1。顯然，該算法包括以下3大模塊：模塊1：初始情況模塊2：計算電勢模塊3：計算概率，決定下一個擊穿點圖2-1 程序流程圖2.1 模塊1 初始情況相關(guān)程序代碼如下：clc %清除命令窗口中顯示的所有輸入和輸出tic; %啟動秒表計時器

12、，輸出tic使用后的流逝時間，與toc結(jié)合記錄程序運(yùn)行時間lx=25;ly=25;lz=25; %定義三維矩陣的最大行、列、層數(shù) v1=zeros(lz,ly,lx); %三維全零矩陣，用來存放電位值v=v1; %v用來確定邊界r=3; %發(fā)展概率指數(shù)c=3e5; %分級前伸速度T=0; %時間矩陣tt=0; %返回雷電先導(dǎo)歷時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,:,:)=1;v(lz,:,:)=1;v(:,1,:

13、)=1;v(:,ly,:)=1;v(:,:,1)=1;v(:,:,lx)=1; %確定邊界b=0; %記錄電位的迭代次數(shù)i=lz-1; %第一個點的高度，即z坐標(biāo)j=fix(ly/2); %第一個點的y坐標(biāo)k=fix(lx/2); %第一個點的x坐標(biāo)axis(0,25,0,25,0,25) %設(shè)置x軸與y軸的限度xlabel(X-axis),ylabel(Y-axis),zlabel(Z-axis)%為x,y,z坐標(biāo)軸命名title(三維模擬得到的先導(dǎo)通道電勢示例) %為圖形添加標(biāo)題line(10,10,10,10,0,7,Color,r,linewidth,3);%畫一條紅色的線代表避雷針2

14、.2 模塊2 計算電勢命xi為第i點電勢，Laplace方程的離散形式可以寫成： (4-1)則對全部n個點，可以列n個方程，同方程組(3-4)。方程組中，第i個方程中xj前的系數(shù)aj如此給定：當(dāng)j=i時，aj=-6；當(dāng)j為i鄰近的六個點之一時，aj=1；其他情況，aj=0。即有 n個方程,每個方程中的n個元素中,除了這個點和與它最鄰近的6個點外,都為0元；把給定的m個邊界條件(包括擊穿點)的電勢值代入方程組，剩余的n-m獨(dú)立方程，恰能給出剩余的n-m個未知網(wǎng)格的電勢；這樣就把求各點電勢問題轉(zhuǎn)化為求解線性方程組的問題了。迭代程序如下：while(maxt1e-5) b=b+1; maxt=0;

15、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(tmaxt) maxt=t; end end end end end v1=v2; end為了加快迭代速度，可以提前使用新值進(jìn)行迭代，因為在上述程序中，我們的掃描賦值方式是一行一行掃描，在對v2(i，j)賦值時，它周圍四個點其中有2個已經(jīng)被掃描過了，即已經(jīng)獲得了

16、新的數(shù)值，這個數(shù)值應(yīng)該更優(yōu)，所以我們盡量使用新算出的數(shù)值進(jìn)行迭代，因此，只要把上述代碼中：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,n-1)+v2(l-1,m,n)/6; 就可以將迭代次數(shù)減少。2.3 模塊3 計算概率，決定下一個擊穿點同樣是基于分布函數(shù)，再由rand(1,1,1)生成一個隨機(jī)數(shù)，這個數(shù)必然在0,1區(qū)間中隨機(jī)選取,而

17、且這個數(shù)必然對應(yīng)于唯一一個可能擊穿的點，由此即可判定x所在的區(qū)間，進(jìn)而判定pi以及所對應(yīng)的ai，這樣，我們就確定了下一個擊穿點，命這個擊穿點=0。相關(guān)代碼與二維類似2，此處不再贅述。3 動畫效果和電勢的變化3.1 動畫效果三維動畫采用程序動畫，相關(guān)代碼如下： head=line(x0,k,y0,j,z0,i,linestyle,-,erasemode,none,linewidth,2);%連接擊穿點與引發(fā)此點擊穿的點drawnow;對于三維圖形，為了更直觀地觀察圖形，可以采用view函數(shù)來改變視角 MATLAB 動畫制作及舉例應(yīng)用EB/OL.百度文庫. 。與z=0平面所成的方向角叫仰角，與x=

18、0所成的夾角叫方位角。這樣，默認(rèn)的三維視角方向仰角為30度，方位角為-37.5度，在MATLAB中，view函數(shù)改變所有類型的二維和三維圖形的圖形視角。view(az,el)和view(az,el)將視角改變到所指定的方位角az和仰角el中。因此view(0,0)可顯示x-z軸視圖，view(90,0)可顯示y-z軸視圖。圖3-1為模擬得到的三維雷電先導(dǎo)路徑的最終結(jié)果，圖3-2為x-z軸視圖，圖3-3為y-z視圖。其實程序運(yùn)行時同樣是動態(tài)變化的。3.2 電勢的變化為了更好地理解雷電先導(dǎo)的特點，我們可以通過以下程序?qū)崿F(xiàn)在圖形上顯示電勢的動態(tài)變化。x,y,z=meshgrid(0:25,0:25,

19、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) %顯示表征顏色刻度的塊，為當(dāng)前軸添加垂直塊 drawnow; %更新圖形窗口head=line(x0,k,y0,j,z0,i,linestyle,-,erasemode,none,linewidth,2);%連接擊穿點與引發(fā)此點擊穿的點drawnow; 三維電勢圖，可以采用slice命令實現(xiàn) MATLAB顏色的使用EB/OL.收稿日期:2013年9月30日作者簡介：魏錐（1990），女，本科，主要研究方向為繼電保護(hù)，740676027王之猛（1986），男，本科，從事繼電保護(hù)工作該命令顯示通過立體圖形的矩形切片圖。slice(V,sx,sy,sz)顯示三元函數(shù)V=V(X,Y,Z)確定的超立體形在x軸、y軸與z軸方向上的若干點(對應(yīng)若干平