正六边形电位分布的有限差分算法.doc

正六边形电位分布的有限差分算法

姓 名：
指导老师：
学号： 目 录
正六边形电位分布的有限差分算法 5
引言 5
1有限差分法的基本原理步骤 5
2正六边形二维场域电位的有限差分算法 6
2.1边界处理 6
2.2数学模型的建立 7
2.3正六边形电位分布的仿真程序 9
3讨论及分析 13

正六边形电位分布的有限差分算法
摘要：介绍了应用有限差分法求电位分布的一般步骤，针对静电场中轴对称情形下的正六边形场域的电位分布，建立了正三角形网格划分的有限差分法的计算模型，给出了Matlab仿真的程序设计流程图，并通过编程得到场域内的电位分布图形，对有限差分法的计算处理进行了讨论和分析。
关键词：正六边形；matlab；电位；有限差分法
引言
关于电磁场的数值计算，常用的方法有：有限差分法，时域有限差分法，有限单元法，矩量法，边界元素法等。由于这些方法只能获得近似解，因此，利用软件进行仿真或者求数值解就显得非常必要。本文基于Matlab，利用有限差分法求解静电场中的正六边形的电位分布。
1有限差分法的基本原理步骤
有限差分法的基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续的定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用求和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。然后利用差值方法便可以从离散解得到定解问题在整个区域上的近似解。有限差分法数值计算包括下列基本步骤：
1)区域的离散或子区域划分；
2)插值函数的选择；
3)方程组的建立；
4)求解方程组。
2正六边形二维场域电位的有限差分算法
轴对称电磁场问题是电工设备设计分析中常遇到的一大类问题。如图1所示，边长为b的正六边形二维场域内无电荷分布，6条边上的电位（V）依次为1，-1，1，-1，1，-1，求场域内的电位分布。
2.1边界处理
由对称性容易看出，正六边形外接圆的D 条直径EE’，FF’和GG’ 均为零电位线。因此，被这3条直径切割成的6个正三角形区域的电位函数不独立，而具有如下性质。
A点的电位与B点的电位满足
(1)
同理，有
(2)
(3)
其中，OA=OB=ρ。当然，即使在-30º≤φ<30º范围内，电位数据仍存在冗余现象。

所以，本题的正六边形二维场域电位分布的计算问题可以化为一个正三角形的电位分布问题，只要求出一个正三角形的电位分布，其它的就可以由（1）、（2）、（3）式的关系来确定。而一个正三角形场域中的电位计算，等价于下述拉普拉斯方程边值问题。
场域：φ=±30º和x=√3b/2≡x03条直线围成的等边三角形区域，如图1所示的△OGE’。
边界条件：ψ=0（当φ=±30º）
ψ=1（当x= x0）
但是，该问题是三角形场域，因此，如果用通常的正方形网格划分边界，那么边界就不能恰好地落在网格上，这样一方面给计算编程带来麻烦，其次会使计算产生边界取值的误差。所以，针对场域形状采用三角形网格划分是处理边界条件的好办法。一般在进行网格划分时采用对称性网格形式，这样既方便数学建模，也方便计算编程。

图1 正六边形场域的边值问题 图2 场域的正三角形网格划分

2.2数学模型的建立
二维场域的拉普拉斯方程可以用有限差分法进行近似计算。首先把求解的区域划分成网格，再把求解区域内连续的场分布用求网络节点上离散的数值解代替。网格必须划分得充分细，才能达到足够的精度。如图2所示，对于正三角形场域△OGE’，采用正三角形网格划分。其边界全部由网格点来划分，避免了边界取值的误差，也方便了计算编程。但域中任一点P的相邻点有6个，因此，用有限差分法计算编程需另建数学模型。
设每个正三角形网格边长为a（称为步长），网格节点（i，j）的电位为ψi,j，与其保持等距离的6个邻点的电位分别为ψi,j+1，ψi,j-1，ψi-1,j，ψi-1,j-1，ψi+1,j，ψi+1,j+1。在a充分小的情况下，可以ψi,j为基点进行泰勒级数展开。
（4）
（5）
（6）
（7）
（8）
（9）
其中：为沿l方向的方向导数，为沿方向的方向导数。由于方向导数可表为

所以，可得ψ的二次方向导数为
（10）
同理
（11）
把（4）-（9）式相加，得

把（10）、（11）式代入上式得
（12）
对于（12）式，由于拉普拉斯方程为

所以（12）式变为

表示a的4阶无穷小，可以略去不计，则有限差分的数学表达式为
（13）
2.3正六边形电位分布的仿真程序
对图1网格节点数设置为37×19=703，迭代精度为10-6。根据（13）式利用Matlab编制程序在计算机上运行，计算程序流程图如图3所示。程序运行得到计算结果的迭代次数为65。图4为程序计算结果的图示，它描述了正六边形二维区域内电位的等位线分布情况，其分布结果一目了然。

图3 计算程序流程图
1） 采用简单迭代法求解
简单迭代法的特点是用前一次迭代得到的网络点电位作为下一次迭代时的初值，迭代时计算公式为

程序运行得到计算结果的迭代次数为65，最后电位数值解收敛于某一固定值。表1、表2列出了上、下两个正三角形节点电位的差分运算结果。原则上，知道正六边形二维场域里的任意一个正三角形中的电位分布，由（1）-（3）式可得到另外的5个正三角形的电位分布，但正六边形二维场域里的6个正三角形的电位边值毕竟有两种：一种是外边界值为1，一种是外边界值为-1。 因此，只要把这两种边界的正三角形内的电位值算出，剩下的4个正三角形内电位分布就完全类似了，这样就能直观地分析正六边形中的电位分布情况。本文以上、下正三角形为两种不同边值的场域为例，列出计算结果比较在对角线上的电位值为0，相应的对称点的电位绝对值相等，符号相反，并沿着x轴越靠近中心，其电位绝对值越小，最后中心点O的电位也为零。

表1 上三角形节点电位的差分运算结果
0
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
0

0
-0.41
-0.65
-0.75
-0.79
-0.81
-0.82
-0.82
-0.81
-0.80
-0.78
-0.75
-0.70
-0.64
-0.54
-0.39
0

0
-0.28
-0.46
-0.57
-0.62
-0.65
-0.65
-0.64
-0.62
-0.58
-0.53
-0.46
-0.37
-0.23
0

0
-0.20
-0.35
-0.43
-0.48
-0.49
-0.49
-0.46
-0.42
-0.36
-0.27
-0.16
0

0
-0.15
-0.26
-0.32
-0.35
-0.35
-0.32
-0.28
-0.21
-0.12
0

0
-0.11
-0.18
-0.22
-0.23
-0.21
-0.16
-0.09
0

0
-0.08
-0.21
-0.13
-0.11
-0.07
0

0
-0.04
-0.06
-0.04
0

0
-1.00
0

0

表2 下三角形节点电位的差分运算结果

0

0
1.00
0

0
0.04
0.06
0.04
0

0
0.08
0.21
0.13
0.11
0.07
0

0
0.11
0.18
0.22
0.23
0.21
0.16
0.09
0

0
0.15
0.26
0.32
0.35
0.35
0.32
0.28
0.21
0.12
0

0
0.20
0.35
0.43
0.48
0.49
0.49
0.46
0.42
0.36
0.27
0.16
0

0
0.28
0.46
0.57
0.62
0.65
0.65
0.64
0.62
0.58
0.53
0.46
0.37
0.23
0

0
0.41
0.65
0.75
0.79
0.81
0.82
0.82
0.81
0.80
0.78
0.75
0.70
0.64
0.54
0.39
0

0
1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00
0

图4为程序计算结果的图示，它描述了正六边形二维区域内电位的等位线分布情况，其分布结果一目了然，与其物理分析的结果一致。源程序代码见附录。

图4 正六边形二维区域内的电位分布图

2） 采用超松弛法（SOR）求解
简单迭代法在解决问题时收敛速度比较慢，一般来说，实用价值不大。实际中常采用超松驰法， 相比之下它有两点重大的改进。 第一是计算每一网格点时，把刚才计算得到的邻近点的电位新值代入，即在计算点（i，j）的电位时，把点（i-1，j）、点（i，j-1）和点（i-1，j-1）的电位用刚才算得的新值代入，即

上式称为松驰法或赛德尔法（relaxation method）。由于提前使用了新值，使得收敛速度加快。
第二，再把上式写成增量形式

这时每次的增量（即上式右边的第二项）就是要求方程局部达到平衡时应补充的量。为了加快收敛，我们引进一个松驰因子w，将上式改写为

SOR法中w的取值对迭代公式的收敛速度影响很大，它的好坏直接影响到加速的快慢。为了保证迭代过程的收敛，必须要求0
w
迭代次数
1.7582
88
1.6500
57
1.5500
39
1.4500
29
1.4200
27
1.3900
27
1.3800
28
1.3700
29
1.3500
31
1.3000
35

由表3可以看到，当w为1.4左右时，收敛速度可以达到最快，迭代次数仅27次。
3） 在该正六边形场域内加矩形导体
为了与实际的电磁场问题结合起来，在该正六边形场域内加矩形导体，导体表面的电位（V）为1，计算此时的电位分布，其结果如图5所示。

图5 正六边形二维区域内加矩形导体的电位分布图

3讨论及分析
（1）传统的二维场域的电位数值差分计算，都是把整个场域网格划分，再进行计算，而很少考虑场域的对称性。而本文则是根据场域的对称性，尽量把场域的计算区域变小，这样可减少计算机的计算时间。对正六边形场域处理就是把整个区域归结为一个正三角形场域△OGE’的电位计算，其计算时间大大缩小。
（2）传统的二维场域网格划分一般是矩形或正方形。本文采用了正三角形网格划分，这是作者基于正六边形场域所作的特殊处理，因为，这样做可以使一些网格点落在边界上，边界网格点的值可准确确定，因此，可提高计算精度，由文计算结果可见，其正六边形场域电位计算的精度是很高的。由此说明，为提高计算精度，除网格划分变细之外，还可以采用特殊形状网格划分。
（3）场域网格划分的不同，其数值计算的数学模型也不同。传统的矩形或正方形网格划分便于数学建模及数值编程计算，本文采用正三角形网格划分，对数学建模及数值编程带来一定的麻烦，因此，在方向导数概念的基础上建立了某点相邻电位之和的1/6模型，这是与传统的四方形网格数值计算模型所不同的，并且其截断误差在4阶无穷小层次上，具有比较高的精度。
（4）迭代方法的不同，会带来收敛速度的不同。本文采用了超松弛迭代法，讨论了不同的松弛因子的选择对收敛速度的影响，并与普通的迭代方法进行了比较。当取到最优松弛因子时，迭代次数可以低至27次，收敛速度大大加快。
本文正是在如上所述的计算原则指导下，对正六边形二维场域边值问题进行了特殊处理并建立相应的数学建模，通过Matlab编程计算，只用了37×19=703个网格点及迭代27次就得到了精度为10-6的数值计算结果并大大地节省了机时，这是我们所预想和需要的。总之，通过对正六边形二维场域边值问题的数值计算处理，体现了数值计算中具体问题具体分析的原则（即特殊问题特殊处理的原则），只有这样，才能在实际的计算过程中有效地解决更多的问题。