用计算机模拟闪电形成的尝试 By 金秀儒 物理三班 PB05206218.

Presentation on theme: "用计算机模拟闪电形成的尝试 By 金秀儒 物理三班 PB05206218."— Presentation transcript:

1 用计算机模拟闪电形成的尝试 By 金秀儒 物理三班 PB

2 Outline 摘要 引言 总体设想 计算机模拟程序 遗憾 & 展望 致谢 参考文献 程序源代码 设想流程图 模拟结果 结果讨论 絮话
模块1 模块2 模块3 计算机模拟程序 模拟结果 结果讨论 絮话 分形 混沌 内禀随机性 遗憾 & 展望 致谢 参考文献 程序源代码

3 摘要 本文提出的模型属于DBM模型的变体，用以模拟自然界中的闪电现象。模型首先通过在二维空间划分网格把实际问题离散化，然后通过适当的算法对闪电形成过程进行模拟，由于计算机时过长的原因，本文最终只能舍弃精确的算法，得到一个初步的甚至是不完整的模型进行计算机模拟，然后进行了一些定性的讨论。

4 引言 Niemeyer,Pientronero和Wielsmann于1983年研究了平行玻璃之间SF6中的电击现象，他们得到了多枝杈的径向放电图形（1µs的30kV电脉冲），他们对图形进行了计算机模拟，提出的介电击穿模型(dielectric breakdown model，DBM)。闪电现象属于气体电介质击穿现象，用类似DBM模型的方法进行计算机模拟是比较自然的想法，本文提出的模型即属于DBM模型的变体，用以模拟自然界中的闪电现象。

5 总体设想 在雷雨季节，大块的云层顶部带正电，而底部则有过剩的负电，于是在接近地面时，地面感应产生正电，地面为高电势，产生的电势差使得空气被击穿，产生闪电。我们以云端为起点，地面为终点模拟闪电的形成。 方便起见，我们开辟一个矩形的二维空间，通过在这个二维空间划分网格把实际问题离散化，我们之后的讨论就限定在这个二维空间的网格内。在二维情况把原则性的问题解决了,推广到三维情况就只是技术问题了。 由此，我们可以命云层的电势为Φ=0，所有的击穿点电势也应为0,地面的电势为常值，为了避免复杂的量纲计算，我们可以采用自然单位制，即不论实际电势差多少，命地面的电势为Φ=1，（实际程序中为方便起见命为100也可以），不失一般性，我们命上边界（网格第一行）的电势Φ=0，以模拟云层，而下边界（网格最后一行）的电势Φ=1，以模拟地面。击穿点的电势为Φ=0

6 总体设想 忽略一些过于细节的影响，我们认为除云层和地面外,空间没有电荷分布。另一方面，为完全确定边界条件,把所有击穿点视作变化的边界,又由于电击对两侧边界电势影响很小而忽略其影响，则可假定两侧的边界的电势自始不变，问题简化为Dirichlet 问题； 联立描写静电场的两个基本方程可得到Poisson方程， 如前所述，全边界给定，空间没有电荷分布， 则Poisson方程 就变为LapLace方程 ； 于是，求解LapLace方程，即可得到全空间电势； 我们把最初的击穿点设为第二排中间一格，重解LapLace方程，即可得到全空间电势；再通过概率的方法决定下一个击穿点(下一个击穿点有可能来自所有已击穿点的相邻单元，否则不可能出现分叉闪电)，重解LapLace方程，得到全空间电势； 如是循环，即模拟闪电作为电介质击穿的过程，最后到达边界，结束程序。

7 设想流程图： 显然，该算法包括以下3大模块： 模块1: 模块2 模块3 初始情况： 计算概率，决定下一个击穿点： 计算电势： 击穿第一个点
通过离散的Laplace方程 重新计算全平面电势 计算所有相邻点击穿概率 取随机数 决定下一个击穿点 击穿点是否到达边界？ NO YES 程序结束，输出图形 初始情况 （START！）

8 模块1(初始情况) 如前所述，上边界（网格第一行）的电势Φ=0，以模拟云层，而下边界（网格最后一行）的电势Φ=1，以模拟地面。
初始情况尤如一个平行板电容器，电势随所在行数，呈线性变化；然后假定两侧的边界（第一列和最后一列）的电势自此之后不再变化，全边界被完全确定，问题简化为Dirichlet 问题;

9 模块2(计算概率，决定下一个击穿点 ) 根据介电击穿模型(dielectric breakdown model，DBM)，我们可以利用一个现成的结果，即击穿点（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)作为概率分布函数有， 由此，我们便可以由计算机实现对下一个击穿点的决定：

10 模块2(计算概率，决定下一个击穿点 ) 由此，我们便可以由计算机实现对下一个击穿点的决定：
把所有最邻近点取出，共n个点，任意命名为a1，a2，a3…ai…an; 由计算机生成一个伪随机数，这个数必然在[0,1]区间中随机选取； 这个数必然对应于唯一一个，由此即可判定x所在的区间，进而判定pi以及所对应的ai； 这样，我们就确定了下一个击穿点， 命这个击穿点Φ=0；

11 模块3(计算电势 ) 右图表示的是数值计算中的一个基本的结果: Laplace方程的离散形式
它表明了一个貌似简单的结果，即某点的电势等于与它相连的四点电势的平均值； 本模块要实现的是计算各离散点的电势，因此，通过Laplace方程的离散形式以及之前完全确定的边界条件（击穿点视为边界的一部分），就可以完全确定网格中所有点的电势； 最自然的想法是直接通过递归的方法用Average函数求解，但实现过程中发现，这个方法是完全行不通的，因为它会产生所谓循环引用的问题，即若干步后所求变量计算时调用自己的值，使程序进入死循环。 这个问题困扰作者许久，直到想到以下这个貌似“十分显然”的方法：

12 模块3(计算电势 ) 命xi为第i点电势，把Laplace方程的离散形式写成 则对全部n个点，可以列n个方程：
方程组中，第i个方程中xj前的系数aj如此给定： 当j=i时，aj= -4；当j为i邻近的四个点之一时，aj=1；其他情况，aj=0， 即有 n个方程,每个方程中的n个元素中,除了这个点和与它最邻近的4个点外,都为0元; 把给定的m个边界条件（包括击穿点）的电势值代入方程组，剩余的n-m独立方程，恰能给出剩余的n-m个未知网格的电势； 这样就把求各点电势问题转化为求解线性方程组的问题了；

13 模块3(计算电势 ) 接下来就是该线性方程组的求解问题， 这个问题最终没有解决，原因不是理论上的，而是来自运行过程中计算机时过长，以致无法求解，说明如下： 如果直接用Cramer法则求解，即使是100阶的方程（相当于10X10网格）乘除法次数约为(101∙100!∙99),即使每秒计算1033也需要近10120年，必须优化解法； 改用直接法计算（如高斯消元法），一般只能算200阶以下的方程（相当于10X20网格），加上需反复计算多次以及模块2所用时间仍然无法接受，也需要优化解法； 由于方程中多为0元，为稀疏方程组，故可以用雅可比迭代法计算，但由于要判断迭代是否收敛以及收敛速度的问题，实现时可能需要用并行算法，皆超出作者能力范围，而且计算机的性能上也不符合条件，需要进一步优化解法，需要加强计算机功能； 综上，由于缺乏可行的实施方法，以及时间和精力的限制，作者最终只好放弃进一步的修正，程序未涉及电势变化，只是初步的甚至是不完整的程序，故本文题名为“尝试”，遗憾终究是不免的；

14 计算机模拟程序 通过上述分析,以下用计算机模拟闪电的形成; 这只是一个初步的结果 程序未涉及电势变化 程序未推广到三维情况

15 模拟结果 (10 X 10像素 )

16 模拟结果 (10 X 10像素 )

17 模拟结果 (10 X 10像素 )

18 模拟结果 (10 X 10像素 )

19 模拟结果 (20 X 20像素 )

20 模拟结果 (20 X 20像素 )

21 模拟结果 (20 X 20像素 )

22 模拟结果 (20 X 20像素 )

23 模拟结果(50 X 50像素, 100 X 100像素) 由于电脑功能尚不够强大，运行容易导致系统崩溃，无法实现。
实际上,程序运行时很可能出现这种情况,即程序运行一段时间后击穿区域先抵达两恻的边界,导致结果完全不可信,唯一的解决方法是横向扩展区域,这样做就对计算机提出了更高的要求 。

24 结果讨论 从以上模拟结果可以看到，10 X 10像素以及20 X 20像素的模拟结果基本令人满意;
作者原计划是把像素调到几十万，这样肉眼就无法分辨出方格了，结果或许可以和真实的闪电媲美，但现在看来，要求这样性能的计算机时，是无法实现的； 事实上，由于在本次模拟中尚未考虑电势的影响，是一个不完善的模拟； 以上这些都有待进一步的工作。 虽然我们得到的只是一个不完善的模拟，但仍可以拿来做一些定性的分析如下： 往下击穿的概率大于往两边击穿的概率，往两边击穿的概率大于往上击穿的概率； 越是接近击穿尖端的点被击穿的概率越大； 由于图形自上而下都是从最高点击穿然后分出许多 “枝杈”，与真实闪电形状近似，故由此可以反过来验证地面为高电势（云层底部带负电荷，地面感应出正电荷）的假设；

25 遗憾 & 展望 由于缺乏计算电势模块，程序是不完整的，相信运用雅可比迭代或其他优化的算法，加上性能强大的计算机，或许可以做出较完整的精确的程序; 奈何造化之工,终究不是人力所及! 或许类比到其他方面后别有应用之处。 比如模拟法求解微分方程。

26 絮话(分形) 闪电作为分形的一种,具有某些分形的特征. 具有精细结构: 分形集都具有任意小尺度下的比例细节 . 具有自相似形式 :
部分与整体以某种形式相似. 这两点也正是本文进行模拟的重要基础. 我们不能说天地之间只有400个区域，甚至也不能说天地之间有几十万个区域，但我们可以认为在开始的一刹那，云端的一个宏观足够小，微观足够大（不至于引起量子效应）的区域，那个过程确实是这样的，而天地之间的过程与那个过程是自相似的。

27 絮话(混沌) 闪电作为混沌的一种表象,具有某些分形的特征. 对于初始值的极端敏感依赖性
只要作为实验出发点的初始值有一个微不足道的差异，在混沌状态下同一种过程将导致截然不同的图像。 不可预言性 由于不可能以无限的精度测量初始值，因此不可能预言任何混沌系统的最后结果。 吸引子 尽管混沌看起来是杂乱无章的，但仍然具有某种条理性。这种内在有序性的源泉是一种被数学家称之为吸引子的东西，它因具有倾向于把一个系统或一个方程吸引到某个终态而得名。 确定性系统中产生了随机性，而微观的随机性又在宏观上表现出某种相似性，这正是本文的精彩之处。

28 絮话(内禀随机性) 通过以上的模拟以及现实的经验,我们是否会有这样的感觉:闪电的形成在微观上类似于一种随机过程,而宏观上又表现出某种规律性,描述确定性的系统,却不可避免的引入了随机性,这种随机性往往被认为是确定性系统中的一种内禀随机性。 这种随机性是来自于我们理论体系描述方式的先天不足,还是客观现象本就如此? That is the question…

29 致谢 在本文写作过程中，得到了很多的帮助，韦威同学在编程上提供了大力支持，张景拓、董鹤飞、熊俊、李翔等同学参与了讨论，提出了有益的意见，在此表示感谢！ 指导老师：程福臻、林宣滨

30 参考文献 DBM模型(PRL): 一些教材: 相关文献:
Niemeyer,L. Pietronero L and Wiesmann H J. Phys Rev Lett,1984,52:1033 一些教材: 电磁学 胡友秋 程福臻 刘之景 高等教育出版社1992 分形原理及其应用 孙霞 吴自勤 黄(田匀) 中国科技大学出版社2003 概率论与数理统计 陈希孺 中国科技大学出版社1992 计算物理学 马文琻 科学出版社2002 数值计算方法和算法 张韵华 奚梅成 陈效群 科学出版社 2000 相关文献: 用带电悬浮粒子云模拟雷电放电的研究 I.维列夏金(莫斯科动力学院) 吴维韩(清华大学,北京) 工科物理 1998年副刊 介电击穿模型(DBM)的相变问题研究:重正化群方法 常福宣 李后强 刘德 林理彬 大自然探索 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月

31 程序源代码：（Generated by Microsoft Visual Basic 6.0）
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 '是否占据 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 Private Sub Command1_Click() Prework Showgrid Dim k As Integer k = m \ 2 Continue = True Breakzero k, 1 Showgrid '显示一下 While Continue Mainwork Wend End Sub Public Sub Prework() Dim i As Integer Dim j As Integer countloop = 0 e1 = 0 Public Sub Addnode(i, j As Integer) s(i, j) = e(i, j) TotalE = TotalE + e(i, j) iii = iii + 1 If iii = m + 1 Then jjj = jjj + 1 iii = 1 End If Nod(iii, jjj).i = i Nod(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 Then If o(ii, jj) = fasle And b(ii, jj) = False And s(ii, jj) = -1 Then Addnode ii, jj ii = i - 1 jj = j If ii <> 0 Then If o(ii, jj) = False And b(ii, jj) = False And s(ii, jj) = -1 Then jj = j + 1 '边界电势 For i = 1 To m e(i, 1) = e1 e(i, n) = e2 Next i '初始化矩阵0 For j = 2 To n - 1 e(i, j) = e(i, j - 1) + (e2 - e1) / (n - 1) Next j '初始化占据矩阵 For j = 1 To n o(i, j) = False o(i, 1) = True '初始化电击 b(i, j) = False Picture1.ScaleHeight = Picture1.Height Picture1.ScaleWidth = Picture1.Width c(1) = vbBlack c(3) = RGB(255, 0, 0) c(2) = RGB(0, 0, 255) '兰 c(4) = vbGreen Randomize End Sub

32 程序源代码：（Generated by Microsoft Visual Basic 6.0）
If j <= m Then If o(ii, jj) = False And b(ii, jj) = False And s(ii, jj) = -1 Then Addnode ii, jj End If ii = i + 1 jj = j If ii <= m Then If o(ii, jj) = fasle And b(ii, jj) = False And s(ii, jj) = -1 Then End Sub Public Sub Expand() Dim i, j As Integer Dim flag As Boolean For j = 1 To n flag = False For i = 1 To m '从行开始扫描，节省时间 If b(i, j) = True Then CountG i, j flag = True Next i If flag = False Then Exit Sub Next j 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 Then Continue = 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 Mainwork() 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 m For j = 1 To m s(i, j) = -1 Next j Next i TotalE = 0 iii = 0 jjj = 1 ' Expand '扩展节点 '算概率 Prepg = 0 If jjj > 1 Then k = m Else k = iii End If For 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(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 Next j '抽随机数 G = Rnd() flag = False

33 程序源代码：（Generated by Microsoft Visual Basic 6.0）
For i = 0 To m - 1 Picture1.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 Sub Public Sub Updatelabel() Label2.ForeColor = c(2) Label3.ForeColor = 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 m For j = 1 To n Dx = Picture1.ScaleWidth / m - 1 Dy = Picture1.ScaleHeight / n - 1 Sx = 1 * i + (i - 1) * Dx Sy = 1 * j + (j - 1) * Dy a1 = Fix((e(i, j) - e1) / maxe * 255) b1 = 0 c1 = Fix((e(i, j) - e1) / maxe * 255) If iii = 1 And jjj = 1 Then ia = Nod(1, 1).i ja = Nod(1, 1).j flag = True Else Prepg = 0 For j = 1 To jjj For i = 1 To k If (G - Prepg) >= Accu And (G - s(Nod(i, j).i, Nod(i, j).j)) <= Accu Then ia = Nod(i, j).i ja = Nod(i, j).j Exit For End If Prepg = s(Nod(i, j).i, Nod(i, j).j) Next i If (jjj - j) = 1 Then k = iii Next j If Not flag Then ' Text1.Text = Text1.Text + "!WQEQEWQ" & vbCrLf ia = Nod(iii, jjj).i ja = Nod(iii, jjj).j Breakzero ia, ja Showgrid End Sub Public Sub Setscalenow() Dim i As Integer Picture1.Scale (0, 0)-(Picture1.ScaleWidth, Picture1.ScaleHeight) If b(i, j) = False Then cor = RGB(a1, b1, c1) Else: cor = c(4) End If Picture1.Line (Sx, Sy)-(Sx + Dx - 1, Sy + Dy - 1), cor, BF Next j Next i End Sub Private Sub Form_Load() e2 = 100 Timer1.Enabled = True m = 20 '预处理 n = 20 debugt = False Prework '预处理 Setscalenow '划坐标 Updatelabel '更新标签 Showgrid '显示 Private Sub Form_MouseMove(Button As Integer, Shift As Integer, x As Single, y As Single) Label4.Caption = "不在画图区域内" Label4.ForeColor = vbBlack Private Sub HScroll1_Change() m = HScroll1.Value n = HScroll1.Value Label5.Caption = "阶数:" + CStr(HScroll1.Value)

34 程序源代码：（Generated by Microsoft Visual Basic 6.0）
Private Sub Label10_Click() m = 200 n = 200 debugt = False Prework '预处理 Setscalenow '划坐标 Updatelabel '更新标签 Showgrid '显示 End Sub Private Sub Label11_Click() m = 500 n = 500 Private Sub Label12_Click() m = 1000 n = 1000 Private Sub Label6_Click() m = 10 n = 10 Private Sub Label7_Click() m = 20 n = 20 debugt = False Prework '预处理 Setscalenow '划坐标 Updatelabel '更新标签 Showgrid '显示 End Sub Private Sub Label8_Click() m = 50 n = 50 Private Sub Label9_Click() m = 100 n = 100 '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 Sub Private Sub Picture1_MouseMove(Button As Integer, Shift As Integer, 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 / m) = 0) Or (Int(newy) Mod (Picture1.ScaleHeight / m) = 0) Then Label4.Caption = "边框" Label4.ForeColor = vbblck Else For i = 1 To m For j = 1 To n Dx = Picture1.ScaleWidth / m - 1 Dy = Picture1.ScaleHeight / n - 1 Sx = 1 * i + (i - 1) * Dx Sy = 1 * j + (j - 1) * Dy If (Int(newx) >= Sx) And (Int(newy) >= Sy) And (Int(newx) <= (Sx + Dx - 1)) And (Int(newy) <= (Sy + Dy - 1)) Then nowx = i nowy = j Exit For End If Next j Next i Label4.Caption = "电势为:" + CStr(e(nowx, nowy)) a1 = Fix((e(nowx, nowy) - e1) / maxe * 255) b1 = 0 c1 = 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) Private Sub Timer1_Timer() sta = sta + 10 While Form1.Height <= 8085 Form1.Height = Form1.Height + sta Wend

35 That’s all . Thank you！