1、用计算机模拟闪电形成的尝试,By 金秀儒 物理三班 PB05206218,Outline,摘要 引言 总体设想 设想流程图 模块1 模块2 模块3 计算机模拟程序 模拟结果 结果讨论 絮话 分形 混沌 内禀随机性 遗憾 & 展望 致谢 参考文献 程序源代码,摘要,本文提出的模型属于DBM模型的变体，用以模拟自然界中的闪电现象。模型首先通过在二维空间划分网格把实际问题离散化，然后通过适当的算法对闪电形成过程进行模拟，由于计算机时过长的原因，本文最终只能舍弃精确的算法，得到一个初步的甚至是不完整的模型进行计算机模拟，然后进行了一些定性的讨论。,引言,Niemeyer,Pientronero和Wie
2、lsmann于1983年研究了平行玻璃之间SF6中的电击现象，他们得到了多枝杈的径向放电图形（1s的30kV电脉冲），他们对图形进行了计算机模拟，提出的介电击穿模型(dielectric breakdown model，DBM)。闪电现象属于气体电介质击穿现象，用类似DBM模型的方法进行计算机模拟是比较自然的想法，本文提出的模型即属于DBM模型的变体，用以模拟自然界中的闪电现象。,总体设想,在雷雨季节，大块的云层顶部带正电，而底部则有过剩的负电，于是在接近地面时，地面感应产生正电，地面为高电势，产生的电势差使得空气被击穿，产生闪电。我们以云端为起点，地面为终点模拟闪电的形成。方便起见，我们开辟
3、一个矩形的二维空间，通过在这个二维空间划分网格把实际问题离散化，我们之后的讨论就限定在这个二维空间的网格内。在二维情况把原则性的问题解决了,推广到三维情况就只是技术问题了。由此，我们可以命云层的电势为=0，所有的击穿点电势也应为0,地面的电势为常值，为了避免复杂的量纲计算，我们可以采用自然单位制，即不论实际电势差多少，命地面的电势为=1，（实际程序中为方便起见命为100也可以），不失一般性，我们命上边界（网格第一行）的电势=0，以模拟云层，而下边界（网格最后一行）的电势=1，以模拟地面。击穿点的电势为=0,总体设想,忽略一些过于细节的影响，我们认为除云层和地面外,空间没有电荷分布。另一方面，为
4、完全确定边界条件,把所有击穿点视作变化的边界,又由于电击对两侧边界电势影响很小而忽略其影响，则可假定两侧的边界的电势自始不变，问题简化为Dirichlet 问题；联立描写静电场的两个基本方程可得到Poisson方程，如前所述，全边界给定，空间没有电荷分布，则Poisson方程 就变为LapLace方程 ；于是，求解LapLace方程，即可得到全空间电势；我们把最初的击穿点设为第二排中间一格，重解LapLace方程，即可得到全空间电势；再通过概率的方法决定下一个击穿点(下一个击穿点有可能来自所有已击穿点的相邻单元，否则不可能出现分叉闪电)，重解LapLace方程，得到全空间电势；如是循环，即模拟
5、闪电作为电介质击穿的过程，最后到达边界，结束程序。,设想流程图：,显然，该算法包括以下3大模块： 模块1: 初始情况： 模块2 计算概率，决定下一个击穿点： 模块3 计算电势：,模块1(初始情况),如前所述，上边界（网格第一行）的电势=0，以模拟云层，而下边界（网格最后一行）的电势=1，以模拟地面。初始情况尤如一个平行板电容器，电势随所在行数，呈线性变化；然后假定两侧的边界（第一列和最后一列）的电势自此之后不再变化，全边界被完全确定，问题简化为Dirichlet 问题;,模块2(计算概率，决定下一个击穿点 ),根据介电击穿模型(dielectric breakdown model，DBM)，我
6、们可以利用一个现成的结果，即击穿点（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)作为概率分布函数有，由此，我们便可以由计算机实现对
7、下一个击穿点的决定：,模块2(计算概率，决定下一个击穿点 ),由此，我们便可以由计算机实现对下一个击穿点的决定：把所有最邻近点取出，共n个点，任意命名为a1，a2，a3aian;由计算机生成一个伪随机数，这个数必然在0,1区间中随机选取；这个数必然对应于唯一一个，由此即可判定x所在的区间，进而判定pi以及所对应的ai；这样，我们就确定了下一个击穿点，命这个击穿点=0；,模块3(计算电势 ),右图表示的是数值计算中的一个基本的结果: Laplace方程的离散形式它表明了一个貌似简单的结果，即某点的电势等于与它相连的四点电势的平均值；本模块要实现的是计算各离散点的电势，因此，通过Laplace方程
8、的离散形式以及之前完全确定的边界条件（击穿点视为边界的一部分），就可以完全确定网格中所有点的电势；最自然的想法是直接通过递归的方法用Average函数求解，但实现过程中发现，这个方法是完全行不通的，因为它会产生所谓循环引用的问题，即若干步后所求变量计算时调用自己的值，使程序进入死循环。这个问题困扰作者许久，直到想到以下这个貌似“十分显然”的方法：,命xi为第i点电势，把Laplace方程的离散形式写成则对全部n个点，可以列n个方程：方程组中，第i个方程中xj前的系数aj如此给定：当j=i时，aj= -4；当j为i邻近的四个点之一时，aj=1；其他情况，aj=0，即有 n个方程,每个方程中的n个
9、元素中,除了这个点和与它最邻近的4个点外,都为0元;把给定的m个边界条件（包括击穿点）的电势值代入方程组，剩余的n-m独立方程，恰能给出剩余的n-m个未知网格的电势；这样就把求各点电势问题转化为求解线性方程组的问题了；,模块3(计算电势 ),模块3(计算电势 ),接下来就是该线性方程组的求解问题， 这个问题最终没有解决，原因不是理论上的，而是来自运行过程中计算机时过长，以致无法求解，说明如下：如果直接用Cramer法则求解，即使是100阶的方程（相当于10X10网格）乘除法次数约为(101100!99),即使每秒计算1033也需要近10120年，必须优化解法；改用直接法计算（如高斯消元法），一
10、般只能算200阶以下的方程（相当于10X20网格），加上需反复计算多次以及模块2所用时间仍然无法接受，也需要优化解法；由于方程中多为0元，为稀疏方程组，故可以用雅可比迭代法计算，但由于要判断迭代是否收敛以及收敛速度的问题，实现时可能需要用并行算法，皆超出作者能力范围，而且计算机的性能上也不符合条件，需要进一步优化解法，需要加强计算机功能；综上，由于缺乏可行的实施方法，以及时间和精力的限制，作者最终只好放弃进一步的修正，程序未涉及电势变化，只是初步的甚至是不完整的程序，故本文题名为“尝试”，遗憾终究是不免的；,计算机模拟程序,通过上述分析,以下用计算机模拟闪电的形成;这只是一个初步的结果程序未涉
11、及电势变化 程序未推广到三维情况,模拟结果 (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像素),由于电脑功能尚不够强大，运行容易导致系统崩溃，无法实现。实际上,程序运行时很可能出现这种情况,即程序运行一段时间后击穿区域先抵达两恻的边界,导致结果完全不可信,唯一的解决方法是横向扩展区域,这样做就对计算机提出了更高的要
12、求 。,结果讨论,从以上模拟结果可以看到，10 X 10像素以及20 X 20像素的模拟结果基本令人满意;作者原计划是把像素调到几十万，这样肉眼就无法分辨出方格了，结果或许可以和真实的闪电媲美，但现在看来，要求这样性能的计算机时，是无法实现的；事实上，由于在本次模拟中尚未考虑电势的影响，是一个不完善的模拟；以上这些都有待进一步的工作。虽然我们得到的只是一个不完善的模拟，但仍可以拿来做一些定性的分析如下：往下击穿的概率大于往两边击穿的概率，往两边击穿的概率大于往上击穿的概率； 越是接近击穿尖端的点被击穿的概率越大； 由于图形自上而下都是从最高点击穿然后分出许多 “枝杈”，与真实闪电形状近似，故由
13、此可以反过来验证地面为高电势（云层底部带负电荷，地面感应出正电荷）的假设；,遗憾 & 展望,由于缺乏计算电势模块，程序是不完整的，相信运用雅可比迭代或其他优化的算法，加上性能强大的计算机，或许可以做出较完整的精确的程序;奈何造化之工,终究不是人力所及!或许类比到其他方面后别有应用之处。比如模拟法求解微分方程。,絮话(分形),闪电作为分形的一种,具有某些分形的特征.具有精细结构: 分形集都具有任意小尺度下的比例细节 .具有自相似形式 : 部分与整体以某种形式相似.这两点也正是本文进行模拟的重要基础.我们不能说天地之间只有400个区域，甚至也不能说天地之间有几十万个区域，但我们可以认为在开始的一刹
14、那，云端的一个宏观足够小，微观足够大（不至于引起量子效应）的区域，那个过程确实是这样的，而天地之间的过程与那个过程是自相似的。,絮话(混沌),闪电作为混沌的一种表象,具有某些分形的特征.对于初始值的极端敏感依赖性 只要作为实验出发点的初始值有一个微不足道的差异，在混沌状态下同一种过程将导致截然不同的图像。 不可预言性 由于不可能以无限的精度测量初始值，因此不可能预言任何混沌系统的最后结果。 吸引子 尽管混沌看起来是杂乱无章的，但仍然具有某种条理性。这种内在有序性的源泉是一种被数学家称之为吸引子的东西，它因具有倾向于把一个系统或一个方程吸引到某个终态而得名。确定性系统中产生了随机性，而微观的随机
15、性又在宏观上表现出某种相似性，这正是本文的精彩之处。,絮话(内禀随机性),通过以上的模拟以及现实的经验,我们是否会有这样的感觉:闪电的形成在微观上类似于一种随机过程,而宏观上又表现出某种规律性,描述确定性的系统,却不可避免的引入了随机性,这种随机性往往被认为是确定性系统中的一种内禀随机性。 这种随机性是来自于我们理论体系描述方式的先天不足,还是客观现象本就如此?That is the question,致谢,在本文写作过程中，得到了很多的帮助，韦威同学在编程上提供了大力支持，张景拓、董鹤飞、熊俊、李翔等同学参与了讨论，提出了有益的意见，在此表示感谢！指导老师：程福臻、林宣滨,参考文献,DBM模
16、型(PRL):Niemeyer,L. Pietronero L and Wiesmann H J. Phys Rev Lett,1984,52:1033一些教材:电磁学 胡友秋 程福臻 刘之景 高等教育出版社1992分形原理及其应用 孙霞 吴自勤 黄(田匀) 中国科技大学出版社2003概率论与数理统计 陈希孺 中国科技大学出版社1992计算物理学 马文琻 科学出版社2002数值计算方法和算法 张韵华 奚梅成 陈效群 科学出版社 2000相关文献:用带电悬浮粒子云模拟雷电放电的研究 I.维列夏金(莫斯科动力学院) 吴维韩(清华大学,北京) 工科物理 1998年副刊介电击穿模型(DBM)的相变问题
17、研究:重正化群方法 常福宣 李后强 刘德 林理彬 大自然探索 1998年 第1期超薄膜外延生长的计算机模拟 吴锋民 朱启鹏 吴自勤 (物理学报 第47卷第9期1998年9月薄膜生长的随机模型 刘祖黎 魏合林 王汉文 王均震 物理学报 第48卷第7期1999年7月射频溅射制备的BST薄膜介电击穿研究 卢肖 吴传贵 张万里 李言荣 物理学报 第55卷第5期2006年5月有限步扩散反应置限分形聚集 吴锋民 朱启鹏 施建青 物理学报 第47卷第4期1998年4月非均质基底表面上团簇生长的Monte Carlo模拟 高国良 钱昌吉 钟瑞 罗孟波 叶高翔 物理学报 第55卷第9期2006年9月,程序源代码
18、：（Generated by Microsoft Visual Basic 6.0）,边界电势 For i = 1 To me(i, 1) = e1e(i, n) = e2 Next i 初始化矩阵0 For i = 1 To mFor j = 2 To n - 1e(i, j) = e(i, j - 1) + (e2 - e1) / (n - 1)Next j Next i 初始化占据矩阵 For i = 1 To mFor j = 1 To no(i, j) = FalseNext j Next i For i = 1 To mo(i, 1) = True Next i 初始化电击 For
19、 i = 1 To mFor j = 1 To nb(i, j) = FalseNext j Next i Picture1.ScaleHeight = Picture1.Height Picture1.ScaleWidth = Picture1.Widthc(1) = vbBlack c(3) = RGB(255, 0, 0) c(2) = RGB(0, 0, 255) 兰 c(4) = vbGreen Randomize End Sub,Public Sub Addnode(i, j As Integer) s(i, j) = e(i, j) TotalE = TotalE + e(i,
20、j) iii = iii + 1 If iii = m + 1 Thenjjj = jjj + 1iii = 1 End IfNod(iii, jjj).i = iNod(iii, jjj).j = j End Sub Public Sub CountG(i, j As Integer) Dim ii, jj As Integer ii = i jj = j - 1 If jj 0 ThenIf o(ii, jj) = fasle And b(ii, jj) = False And s(ii, jj) = -1 ThenAddnode ii, jjEnd If End If ii = i -
21、1 jj = j If ii 0 ThenIf o(ii, jj) = False And b(ii, jj) = False And s(ii, jj) = -1 ThenAddnode ii, jjEnd If End If ii = i jj = j + 1,Const maxN = 1000 Const Accu = 0.01 Dim 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 是
22、否占据 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 TotalE As Double Dim Nod(1 To maxN, 1 To maxN) As node 链表 Dim iii, jjj As Integer 链表指针 Dim Continue As Boolean Dim countloop As Integer Dim debugt As Boolean Dim sta As Integer
23、Private Sub Command1_Click() Prework Showgrid Dim k As Integer k = m 2 Continue = True Breakzero k, 1 Showgrid 显示一下 While ContinueMainwork Wend End Sub Public Sub Prework() Dim i As Integer Dim j As Integer countloop = 0 e1 = 0,程序源代码：（Generated by Microsoft Visual Basic 6.0）,Expand 扩展节点 算概率 Prepg =
24、0 If jjj 1 Thenk = m Elsek = iii End If For j = 1 To jjjFor i = 1 To ks(Nod(i, j).i, Nod(i, j).j) = s(Nod(i, j).i, Nod(i, j).j) / TotalEs(Nod(i, j).i, Nod(i, j).j) = s(Nod(i, j).i, Nod(i, j).j) + PrepgPrepg = s(Nod(i, j).i, Nod(i, j).j) If debugt Then Text1.Text = Text1.Text + CStr(countloop) + “ “
25、+ CStr(s(Nod(i, j).i, Nod(i, j).j) + vbCrLf End IfNext iIf (jjj - j) = 1 Thenk = iiiEnd If Next j 抽随机数 G = Rnd() flag = False If jjj 1 Thenk = m Elsek = iii End If,If j = m ThenIf o(ii, jj) = False And b(ii, jj) = False And s(ii, jj) = -1 ThenAddnode ii, jjEnd If End If ii = i + 1 jj = j If ii = m T
26、henIf o(ii, jj) = fasle And b(ii, jj) = False And s(ii, jj) = -1 ThenAddnode ii, jjEnd If End If End Sub Public Sub Expand() Dim i, j As Integer Dim flag As Boolean For j = 1 To nflag = FalseFor i = 1 To m 从行开始扫描，节省时间If b(i, j) = True ThenCountG i, jflag = TrueEnd IfNext iIf flag = False ThenExit Su
27、bEnd If Next j End Sub,Public Sub Breakzero(ByVal x As Integer, ByVal y As Integer) e(x, y) = 0 o(x, y) = True b(x, y) = True If y = m ThenContinue = False End If If debugt Then Text1.Text = Text1.Text + “decide“ + CStr(x) + “,“ + CStr(y) + vbCrL End If countloop = countloop + 1 End Sub Public Sub M
28、ainwork() Dim Prepg As Double Dim G As Double Dim ia, ja As Integer Dim flag As Boolean Dim k As Integer prework Dim i, j As Integer improvable For i = 1 To mFor j = 1 To ms(i, j) = -1Next j Next i TotalE = 0 iii = 0 jjj = 1 ,程序源代码：（Generated by Microsoft Visual Basic 6.0）,If b(i, j) = False Thencor
29、 = RGB(a1, b1, c1)Else: cor = c(4)End IfPicture1.Line (Sx, Sy)-(Sx + Dx - 1, Sy + Dy - 1), cor, BFNext j Next i End Sub Private Sub Form_Load() e2 = 100 Timer1.Enabled = True m = 20 预处理 n = 20 debugt = False Prework 预处理 Setscalenow 划坐标 Updatelabel 更新标签 Showgrid 显示 End Sub Private Sub Form_MouseMove(
30、Button As Integer, Shift As Integer, x As Single, y As Single) Label4.Caption = “不在画图区域内“ Label4.ForeColor = vbBlack End Sub Private Sub HScroll1_Change() m = HScroll1.Value n = HScroll1.Value debugt = False Prework 预处理 Setscalenow 划坐标 Updatelabel 更新标签 Showgrid 显示 Label5.Caption = “阶数:“ + CStr(HScro
31、ll1.Value) End Sub,If iii = 1 And jjj = 1 Thenia = Nod(1, 1).ija = Nod(1, 1).jflag = True ElsePrepg = 0For j = 1 To jjjFor i = 1 To kIf (G - Prepg) = Accu And (G - s(Nod(i, j).i, Nod(i, j).j) = Accu Thenia = Nod(i, j).ija = Nod(i, j).jflag = TrueExit ForEnd IfPrepg = s(Nod(i, j).i, Nod(i, j).j)Next
32、iIf (jjj - j) = 1 Thenk = iiiEnd IfNext jIf Not flag Then Text1.Text = Text1.Text + “!WQEQEWQ“ & vbCrLfia = Nod(iii, jjj).ija = Nod(iii, jjj).jEnd If End If Breakzero ia, ja Showgrid End Sub Public Sub Setscalenow() Dim i As Integer Picture1.Scale (0, 0)-(Picture1.ScaleWidth, Picture1.ScaleHeight),F
33、or i = 0 To m - 1Picture1.Line (i * Picture1.ScaleWidth / m, 0)-(i * Picture1.ScaleWidth / m, Picture1.ScaleHeight), c(1)Picture1.Line (0, i * Picture1.ScaleHeight / n)-(Picture1.ScaleWidth, i * Picture1.ScaleHeight / n), c(1) Next i End SubPublic Sub Updatelabel() Label2.ForeColor = c(2) Label3.For
34、eColor = c(3) End Sub) Public Sub Showgrid() Dim maxe As Double 最大改变量 Dim i As Integer Dim j As Integer Dim Dx As Double Dim Dy As Double Dim Sx As Double Dim Sy As Double Dim cor As Long 颜色 Dim a1 As Long Dim b1 As Double Dim c1 As Double maxe = -e1 + e2 For i = 1 To mFor j = 1 To nDx = Picture1.Sc
35、aleWidth / m - 1Dy = Picture1.ScaleHeight / n - 1Sx = 1 * i + (i - 1) * DxSy = 1 * j + (j - 1) * Dya1 = Fix(e(i, j) - e1) / maxe * 255)b1 = 0c1 = 255 - Fix(e(i, j) - e1) / maxe * 255),程序源代码：（Generated by Microsoft Visual Basic 6.0）,Private Sub Picture1_MouseMove(Button As Integer, Shift As Integer,
36、x As Single, y As Single) Dim i As Integer Dim Dx, Dy, Sx, Sy As Integer Dim nowx, nowy As Integer Dim a1, b1, c1 As Integer Dim cor As Long Dim maxe As Double Dim newx, newy As Single newx = x newy = y maxe = e2 - e1 Label1.Caption = newx Label2.Caption = newy If (Int(newx) Mod (Picture1.ScaleWidth
37、 / m) = 0) Or (Int(newy) Mod (Picture1.ScaleHeight / m) = 0) ThenLabel4.Caption = “边框“Label4.ForeColor = vbblck ElseFor i = 1 To mFor j = 1 To nDx = Picture1.ScaleWidth / m - 1Dy = Picture1.ScaleHeight / n - 1Sx = 1 * i + (i - 1) * DxSy = 1 * j + (j - 1) * DyIf (Int(newx) = Sx) And (Int(newy) = Sy)
38、And (Int(newx) = (Sx + Dx - 1) And (Int(newy) = (Sy + Dy - 1) Thennowx = inowy = jExit ForEnd IfNext jNext i,Private Sub Label10_Click() m = 200 n = 200 debugt = False Prework 预处理 Setscalenow 划坐标 Updatelabel 更新标签 Showgrid 显示 End Sub Private Sub Label11_Click() m = 500 n = 500 debugt = False Prework
39、预处理 Setscalenow 划坐标 Updatelabel 更新标签 Showgrid 显示 End Sub Private Sub Label12_Click() m = 1000 n = 1000 debugt = False Prework 预处理 Setscalenow 划坐标 Updatelabel 更新标签 Showgrid 显示 End Sub Private Sub Label6_Click() m = 10 n = 10 debugt = False Prework 预处理 Setscalenow 划坐标 Updatelabel 更新标签 Showgrid 显示 End
40、Sub,Private Sub Label7_Click() m = 20 n = 20 debugt = False Prework 预处理 Setscalenow 划坐标 Updatelabel 更新标签 Showgrid 显示 End Sub Private Sub Label8_Click() m = 50 n = 50 debugt = False Prework 预处理 Setscalenow 划坐标 Updatelabel 更新标签 Showgrid 显示 End Sub Private Sub Label9_Click() m = 100 n = 100 debugt = Fa
41、lse Prework 预处理 Setscalenow 划坐标 Updatelabel 更新标签 Showgrid 显示 End Sub Private Sub Form_Resize() Picture1.Left = 30 Picture1.Top = 20 If Form1.ScaleHeight = 100 Then Picture1.Height = Form1.ScaleHeight - 100 End If If Form1.ScaleWidth = 60 Then Picture1.Width = Form1.ScaleWidth - 60 End If End Sub,Lab
42、el4.Caption = “电势为:“ + CStr(e(nowx, nowy)a1 = Fix(e(nowx, nowy) - e1) / maxe * 255)b1 = 0c1 = 255 - Fix(e(nowx, nowy) - e1) / maxe * 255)cor = RGB(a1, b1, c1)Label4.ForeColor = cor End If End Sub Private Sub Text1_Change() e2 = Val(Text1.Text) End Sub Private Sub Timer1_Timer() sta = sta + 10 While Form1.Height = 8085Form1.Height = Form1.Height + sta Wend End Sub,Thats all . Thank you！,
