用计算机模拟闪电形成的尝试课件.ppt

1、用计算机模拟闪电形成的尝试By By 金秀儒金秀儒物理三班物理三班PB05206218PB05206218Outline摘要引言总体设想设想流程图模块1模块2模块3计算机模拟程序模拟结果结果讨论絮话分形混沌内禀随机性遗憾&展望 致谢参考文献程序源代码摘要 本文提出的模型属于DBM模型的变体，用以模拟自然界中的闪电现象。模型首先通过在二维空间划分网格把实际问题离散化，然后通过适当的算法对闪电形成过程进行模拟，由于计算机时过长的原因，本文最终只能舍弃精确的算法，得到一个初步的甚至是不完整的模型进行计算机模拟，然后进行了一些定性的讨论。引言 Niemeyer,Pientronero和Wielsman

2、n于1983年研究了平行玻璃之间SF6中的电击现象，他们得到了多枝杈的径向放电图形（1s的30kV电脉冲），他们对图形进行了计算机模拟，提出的介电击穿模型(dielectric breakdown model，DBM)。闪电现象属于气体电介质击穿现象，用类似DBM模型的方法进行计算机模拟是比较自然的想法，本文提出的模型即属于DBM模型的变体，用以模拟自然界中的闪电现象。总体设想总体设想 在雷雨季节，大块的云层顶部带正电，而底部则有过剩的负电，于是在接近地面时，地面感应产生正电，地面为高电势，产生的电势差使得空气被击穿，产生闪电。我们以云端为起点，地面为终点模拟闪电的形成我们以云端为起点，地面为

3、终点模拟闪电的形成。方便起见，我们开辟一个矩形的二维空间，通过在这个二维空间划分网格把实际问题离散化，我们之后的讨论就限定在这个二维空间的网格内。在二维情况把原则性的问题解决了,推广到三维情况就只是技术问题了。推广到三维情况就只是技术问题了。由此，我们可以命云层的电势为云层的电势为=0，所有的击穿点电势也应为0,地面的电势为常值，为了避免复杂的量纲计算避免复杂的量纲计算，我们可以采用自然单位制，即不论实际电势差多少，命地面的电势为电势为=1，（实际程序中为方便起见命为100也可以），不失一般性，我们命上边界（网格第一行）的电势=0，以模拟云层，而下边界（网格最后一行）的电势=1，以模拟地面。击

4、穿点的电势为击穿点的电势为=0总体设想总体设想 忽略一些过于细节的影响，我们认为除云层和地面外,空间没有电荷分布。另一方面，为完为完全确定边界条件全确定边界条件,把所有击穿点视作变化的边界把所有击穿点视作变化的边界,又由于电击对两侧边界电势影响很小而忽略其影响，则可假定两侧的边界的电势自始不变，问题简化为Dirichlet 问题；联立描写静电场的两个基本方程可得到Poisson方程，如前所述，全边界给定，空间没有电荷分布，则Poisson方程就变为LapLace方程；于是，求解LapLace方程，即可得到全空间电势；我们把最初的击穿点设为第二排中间一格，重解LapLace方程，即可得到全空间电

5、势；再通过概率的方法决定下一个击穿点(下一个击穿点有可能来自所有已击穿点的相邻单元，下一个击穿点有可能来自所有已击穿点的相邻单元，否则不可能出现分叉闪电否则不可能出现分叉闪电)，重解LapLace方程，得到全空间电势；如是循环，即模拟闪电作为电介质击穿的过程，最后到达边界，结束程序。设想流程图：击穿第一个点击穿第一个点通过离散的通过离散的LaplaceLaplace方程方程重新计算全平面电势重新计算全平面电势计算所有相邻点击穿概率计算所有相邻点击穿概率取随机数取随机数决定下一个击穿点决定下一个击穿点击穿点是否到达边界？击穿点是否到达边界？NONOYESYES程序结束，输出图形程序结束，输出图形

6、初始情况初始情况（STARTSTART！）！）显然，该算法包括以下3大模块：模块1:初始情况：模块2计算概率，决定下一个击穿点：模块3计算电势：模块模块1(初始情况初始情况)如前所述，上边界（网格第一行）的电势=0，以模拟云层，而下边界（网格最后一行）的电势=1，以模拟地面。初始情况尤如一个平行板电容器，电势随所在行数，呈线性变化；然后假定两侧的边界（第一列和最后一列）的电势自此之后不再变化，全边界被完全确定，问题简化为Dirichlet 问题;模块模块2(计算概率，决定下一个击穿点计算概率，决定下一个击穿点)根据介电击穿模型(dielectric breakdown model，DBM)，我

7、们可以利用一个现成的结果，即击穿点（i，j）的某一最近邻点（i+1，j）的击穿概率为：式中分母中的求和号遍及所有最近邻点，且（i+1，j）可以替换为（i-1，j）（i，j+1）（i，j-1）；M=1时，模型满足Laplace方程，由击穿点（i，j）=0，即我们所要用的公式为：式中分母中的求和号遍及所有最近邻点；由此式计算所有最近邻点的击穿概率；再由计算机产生的伪随机数决定下一个击穿点；以下为由计算机产生伪随机数决定下一个击穿点的具体方法：我们以列表的形式给出可能的击穿点的概率分布：从中显而易见，作为概率函数的pi满足，由是，我们命F(x)作为概率分布函数有，由此，我们便可以由计算机实现对下一个

8、击穿点的决定：模块模块2(计算概率，决定下一个击穿点计算概率，决定下一个击穿点)由此，我们便可以由计算机实现对下一个击穿点的决定：1.把所有最邻近点取出，共n个点，任意命名为a1，a2，a3aian;2.由计算机生成一个伪随机数，这个数必然在这个数必然在0,1区间中随机选取区间中随机选取；3.这个数必然对应于唯一一个，由此即可判定x所在的区间，进而判定pi以及所对应的ai；4.这样，我们就确定了下一个击穿点，5.命这个击穿点=0；模块模块3(计算电势计算电势)右图表示的是数值计算中的一个基本的结果:Laplace方程的离散形式方程的离散形式它表明了一个貌似简单的结果，即某点的电势等于与它相连的

9、四点电势的平均值；本模块要实现的是计算各离散点的电势，因此，通过Laplace方程的离散形式以及之前完全确定的边界条件（击穿点视为边界的一部分），就可以完全确定网格中所有点的电势；最自然的想法是直接通过递归的方法用Average函数求解，但实现过程中发现，这个方法是完全行不通的，因为它会产生所谓循环引用的问题，即若干步后所求变量计算时调用自己的值，使程序进入死循环。这个问题困扰作者许久，直到想到以下这个貌似“十分显然”的方法：命xi为第i点电势，把Laplace方程的离散形式写成则对全部n个点，可以列n个方程：方程组中，第i个方程中xj前的系数aj如此给定：当j=i时，aj=-4；当j为i邻近

10、的四个点之一时，aj=1；其他情况，aj=0，即有即有 n个方程个方程,每个方程中的每个方程中的n个元素中个元素中,除了这个点和与它最邻近的除了这个点和与它最邻近的4个点外个点外,都为都为0元元;把给定的m个边界条件（包括击穿点）的电势值代入方程组，剩余的剩余的n-m独立方程，恰能独立方程，恰能给出剩余的给出剩余的n-m个未知网格的电势个未知网格的电势；这样就把求各点电势问题转化为求解线性方程组把求各点电势问题转化为求解线性方程组的问题了；0.0.0.0.0.0.1)1(2211)1(1)1)(1()1(22)1(11)1()2(1)1)(2()2(22)2(11)2(31)1(3323213

11、121)1(2222212111)1(11212111nnnnnnininnnnnnnniinnnnnnnnniinnnnnnniinnnniinnnniixaxaxaxaxaxaxaxaxaxaxaxaxaxaxaxaxaxaxaxaxaxaxaxaxaxaxaxaxaxa模块模块3(计算电势计算电势)04 xxxxxdownuprightleft模块模块3(计算电势计算电势)接下来就是该线性方程组的求解问题，这个问题最终没有解决，原因不是理论上的，而是来自运行过程中计算机时过长运行过程中计算机时过长，以致无法求解，说明如下：如果直接用Cramer法则求解，即使是100阶的方程（相当于10X

12、10网格）乘除法次数约为(101100!99),即使每秒计算即使每秒计算1033也需要近也需要近10120年年，必须优化解法；改用直接法计算（如高斯消元法），一般只能算200阶以下的方程（相当于10X20网格），加上需反复计算多次以及模块2所用时间仍然无法接受，也需要优化解法；由于方程中多为0元，为稀疏方程组，故可以用雅可比迭代法计算，但由于要判断迭代是要判断迭代是否收敛以及收敛速度否收敛以及收敛速度的问题，实现时可能需要用并行算法，皆超出作者能力范围，而且计算机的性能上也不符合条件，需要进一步优化解法，需要加强计算机功能；综上，由于缺乏可行的实施方法，以及时间和精力的限制，作者最终只好放弃进

13、一步的修正，程序未涉及电势变化，只是初步的甚至是不完整的程序，故本文题名为“尝试”，遗憾终究是不免的；计算机模拟程序 通过上述分析,以下用计算机模拟闪电的形成;这只是一个初步的结果 程序未涉及电势变化 程序未推广到三维情况模拟结果(10 X 10像素)模拟结果(10 X 10像素)模拟结果(10 X 10像素)模拟结果(10 X 10像素)模拟结果(20 X 20像素)模拟结果(20 X 20像素)模拟结果(20 X 20像素)模拟结果(20 X 20像素)模拟结果(50 X 50像素像素,100 X 100像素像素)由于电脑功能尚不够强大，运行容易导致系统崩溃，无法实现。实际上,程序运行时很

14、可能出现这种情况,即程序运行一段时间后击穿区域先抵达两恻的边界,导致结果完全不可信,唯一的解决方法是横向扩展区域,这样做就对计算机提出了更高的要求。结果讨论结果讨论从以上模拟结果可以看到，10 X 10像素以及20 X 20像素的模拟结果基本令人满意;作者原计划是把像素调到几十万，这样肉眼就无法分辨出方格了，结果或许可以和真实的闪电媲美，但现在看来，要求这样性能的计算机时，是无法实现的；事实上，由于在本次模拟中尚未考虑电势的影响，是一个不完善的模拟；以上这些都有待进一步的工作。虽然我们得到的只是一个不完善的模拟，但仍可以拿来做一些定性的分析如下：往下击穿的概率大于往两边击穿的概率，往两边击穿的

15、概率大于往上击穿的概率；越是接近击穿尖端的点被击穿的概率越大；由于图形自上而下都是从最高点击穿然后分出许多“枝杈”，与真实闪电形状近似，故由此可以反过来验证地面为高电势（云层底部带负电荷，地面感应出正电荷）的假设；遗憾&展望 由于缺乏计算电势模块，程序是不完整的，相信运用雅可比迭代或其他优化的算法，加上性能强大的计算机，或许可以做出较完整的精确的程序;奈何造化之工,终究不是人力所及!或许类比到其他方面后别有应用之处。比如模拟法求解微分方程。絮话(分形)闪电作为分形的一种,具有某些分形的特征.具有精细结构:分形集都具有任意小尺度下的比例细节.具有自相似形式:部分与整体以某种形式相似.这两点也正是

16、本文进行模拟的重要基础.我们不能说天地之间只有400个区域，甚至也不能说天地之间有几十万个区域，但我们可以认为在开始的一刹那，云端的一个宏观足够小，微观足够大（不至于引起量子效应）的区域，那个过程确实是这样的，而天地之间的过程与那个过程是自相似的。絮话(混沌)闪电作为混沌的一种表象,具有某些分形的特征.对于初始值的极端敏感依赖性只要作为实验出发点的初始值有一个微不足道的差异，在混沌状态下同一种过程将导致截然不同的图像。不可预言性由于不可能以无限的精度测量初始值，因此不可能预言任何混沌系统的最后结果。吸引子尽管混沌看起来是杂乱无章的，但仍然具有某种条理性。这种内在有序性的源泉是一种被数学家称之为

17、吸引子的东西，它因具有倾向于把一个系统或一个方程吸引到某个终态而得名。确定性系统中产生了随机性，而微观的随机性又在宏观上表现出某种相似性，这正是本文的精彩之处。絮话(内禀随机性)通过以上的模拟以及现实的经验,我们是否会有这样的感觉:闪电的形成在微观上类似于一种随机过程,而宏观上又表现出某种规律性,描述确定性的系统,却不可避免的引入了随机性,这种随机性往往被认为是确定性系统中的一种内禀随机性。这种随机性是来自于我们理论体系描述方式的先天不足,还是客观现象本就如此?That is the question致谢致谢 在本文写作过程中，得到了很多的帮助，韦威同学在编程上提供了大力支持，张景拓、董鹤飞、

18、熊俊、李翔等同学参与了讨论，提出了有益的意见，在此表示感谢！指导老师：程福臻、林宣滨参考文献DBM模型(PRL):Niemeyer,L.Pietronero L and Wiesmann H J.Phys Rev Lett,1984,52:1033一些教材:电磁学 胡友秋 程福臻 刘之景 高等教育出版社1992分形原理及其应用 孙霞 吴自勤 黄(田匀)中国科技大学出版社2003概率论与数理统计 陈希孺 中国科技大学出版社1992计算物理学 马文琻 科学出版社2002数值计算方法和算法 张韵华 奚梅成 陈效群 科学出版社 2000相关文献:用带电悬浮粒子云模拟雷电放电的研究 I.维列夏金(莫斯科

19、动力学院)吴维韩(清华大学,北京)工科物理 1998年副刊介电击穿模型(DBM)的相变问题研究:重正化群方法 常福宣 李后强 刘德 林理彬 大自然探索 1998年 第1期超薄膜外延生长的计算机模拟 吴锋民 朱启鹏 吴自勤(物理学报 第47卷第9期1998年9月薄膜生长的随机模型 刘祖黎 魏合林 王汉文 王均震 物理学报 第48卷第7期1999年7月射频溅射制备的BST薄膜介电击穿研究 卢肖 吴传贵 张万里 李言荣 物理学报 第55卷第5期2006年5月有限步扩散反应置限分形聚集 吴锋民 朱启鹏 施建青 物理学报 第47卷第4期1998年4月非均质基底表面上团簇生长的Monte Carlo模拟

20、高国良 钱昌吉 钟瑞 罗孟波 叶高翔 物理学报 第55卷第9期2006年9月程序源代码：（程序源代码：（Generated by Microsoft Visual Basic 6.0）边界电势For i=1 To m e(i,1)=e1 e(i,n)=e2Next i初始化矩阵0For i=1 To m For j=2 To n-1 e(i,j)=e(i,j-1)+(e2-e1)/(n-1)Next jNext i初始化占据矩阵For i=1 To m For j=1 To n o(i,j)=False Next jNext iFor i=1 To m o(i,1)=TrueNext i初始化

21、电击For i=1 To m For j=1 To n b(i,j)=False Next jNext iPicture1.ScaleHeight=Picture1.HeightPicture1.ScaleWidth=Picture1.Widthc(1)=vbBlackc(3)=RGB(255,0,0)c(2)=RGB(0,0,255)兰c(4)=vbGreenRandomizeEnd SubPublic Sub Addnode(i,j As Integer)s(i,j)=e(i,j)TotalE=TotalE+e(i,j)iii=iii+1If iii=m+1 Then jjj=jjj+1

22、iii=1End If Nod(iii,jjj).i=i Nod(iii,jjj).j=jEnd SubPublic Sub CountG(i,j As Integer)Dim ii,jj As Integerii=ijj=j-1If jj 0 Then If o(ii,jj)=fasle And b(ii,jj)=False And s(ii,jj)=-1 Then Addnode ii,jj End IfEnd Ifii=i-1jj=jIf ii 0 Then If o(ii,jj)=False And b(ii,jj)=False And s(ii,jj)=-1 Then Addnode

23、 ii,jj End IfEnd Ifii=ijj=j+1Const maxN=1000Const Accu=0.01Dim m,n As Integer 边长Dim c(1 To 4)As Long 颜色Dim e(1 To maxN,1 To maxN)As Double 电势矩阵Dim o(1 To maxN,1 To maxN)As Boolean 是否占据Dim b(1 To maxN,1 To maxN)As Boolean 是否为电击Dim s(1 To maxN,1 To maxN)As Double 概率矩阵Dim e1,e2 As Double 上下界电势Dim Total

24、E As DoubleDim Nod(1 To maxN,1 To maxN)As node 链表Dim iii,jjj As Integer 链表指针Dim Continue As BooleanDim countloop As IntegerDim debugt As BooleanDim sta As IntegerPrivate Sub Command1_Click()PreworkShowgridDim k As Integerk=m 2Continue=TrueBreakzero k,1Showgrid 显示一下While Continue MainworkWendEnd SubP

25、ublic Sub Prework()Dim i As IntegerDim j As Integercountloop=0e1=0程序源代码：（程序源代码：（Generated by Microsoft Visual Basic 6.0）Expand 扩展节点算概率Prepg=0If jjj 1 Then k=mElse k=iiiEnd IfFor j=1 To jjj For i=1 To k s(Nod(i,j).i,Nod(i,j).j)=s(Nod(i,j).i,Nod(i,j).j)/TotalE s(Nod(i,j).i,Nod(i,j).j)=s(Nod(i,j).i,Nod

26、(i,j).j)+Prepg Prepg=s(Nod(i,j).i,Nod(i,j).j)If debugt Then Text1.Text=Text1.Text+CStr(countloop)+CStr(s(Nod(i,j).i,Nod(i,j).j)+vbCrLf End If Next i If(jjj-j)=1 Then k=iii End IfNext j抽随机数G=Rnd()flag=FalseIf jjj 1 Then k=mElse k=iiiEnd IfIf j=m Then If o(ii,jj)=False And b(ii,jj)=False And s(ii,jj)=

27、-1 Then Addnode ii,jj End IfEnd Ifii=i+1jj=jIf ii=Accu And(G-s(Nod(i,j).i,Nod(i,j).j)=Sx)And(Int(newy)=Sy)And(Int(newx)=(Sx+Dx-1)And(Int(newy)=100 ThenPicture1.Height=Form1.ScaleHeight-100End IfIf Form1.ScaleWidth=60 ThenPicture1.Width=Form1.ScaleWidth-60End IfEnd SubLabel4.Caption=电势为:+CStr(e(nowx,

28、nowy)a1=Fix(e(nowx,nowy)-e1)/maxe*255)b1=0 c1=255-Fix(e(nowx,nowy)-e1)/maxe*255)cor=RGB(a1,b1,c1)Label4.ForeColor=corEnd IfEnd SubPrivate Sub Text1_Change()e2=Val(Text1.Text)End SubPrivate Sub Timer1_Timer()sta=sta+10While Form1.Height=8085 Form1.Height=Form1.Height+staWendEnd SubThats all.Thank you！