第五章 场和电子轨迹求解

1、电子光学 n1.绪论 n2. 第一章 几何光学基础 n3.第二章 电子在均匀场中的运动 n4. 第三章 电子光学系统中的场 n5.第四章 电子轨迹方程 n6.第五章第五章 场和电子轨迹的求解场和电子轨迹的求解 n7.第六章 强流电子光学 教师：刘迎辉 电子科技大学物电学院 第五章 场和电子轨迹求解方法 本章组织 n5.1 有限差分法求解轴对称电场 n5.2 有限元法求解轴对称电场 n5.3 电荷密度法求解轴对称电场 n5.4 电子轨迹的数值求解方法 n5.5 有限差分法求解轴对称磁场 n5.6 测量磁场的实验方法。 电子光学系统 将一个给定的场（包括电场、磁 场）看做一个。 ：根据电子在场中的

2、运动，确 定电子光学系统的光学性质和光学参量。 （常用的物理方法：引力场、太阳系、黑洞） 确定电场和磁场的具体分布成为研究、确定电场和磁场的具体分布成为研究、 了解和设计电子光学系统必不可少的重了解和设计电子光学系统必不可少的重 要步骤。要步骤。 静电场和恒定电场的边值问题（物理），可归结静电场和恒定电场的边值问题（物理），可归结 为在给定边界条件下求解拉普拉斯方程或泊松方程为在给定边界条件下求解拉普拉斯方程或泊松方程 （数学）（数学） 。 常用的方法常用的方法 直接法直接法 间接法间接法 解析法解析法 数值法数值法 有限差分法有限差分法（FDM） 有限元方法有限元方法（FEM） 有限积分法有

3、限积分法（FIM） 实验法实验法 场的求解方法 边界积分法边界积分法（BIM） 场的求解方法 u1、实验法求解精度低，如电解槽法、电阻网法 电解槽法利用导电介质中的电流线模拟真空中的 电力线达到求解场分布 u2、解析法求解精度最高，适用简单、规则边界 结构 利用边界条件直接求解偏微分方程 u3、数值法求解精度较高，计算量大，适用范围 广 如有限差分法，将偏微分方程利用差分近似，通 过数值计算方法求得一系列关于场的离散值。 n电真空器件的研制涉及到电子光学、磁学、阴极电子学、微波电子学、电 磁场理论、材料学、机械与热分析诸多学科，工艺过程十分复杂。计算机 技术的发展与应用，极大的促进了微波管技术

4、的进步。它对提高微波管的 设计能力，缩短开发周期，减少整管硬件实验，改善微波管性能，固化已 有经验上发挥着越来越重要的作用，计算机辅助设计CAD技术已经成为微 波管设计的重要手段 软件计算方法应用领域 MAGICFDM三维静态场，高频电磁场，注波互作用 MAFIAFIM三维静态场，涡流场，高频电磁场，注波互作用 ANSYSFEM三维静态场，涡流场 HFSSFEM三维高频电磁场 SOSFDM三维静态场，高频电磁场，注波互作用 ARGUSFDM三维静态场，高频电磁场，注波互作用 MASKFDM三维静态场，高频电磁场，注波互作用 MEGAFEM三维静态场，涡流场，高频电磁场 GUN3DFDM三维电子

5、光学 国际部分常用大型电磁分析软件 nCAD技术已经成为研制新型微波管不可或缺的手段，其总体目标总体目标是：通过模拟， 一次装管成功。 n微波电子学数值模拟的实质是实质是在给定边界条件和初始条件下，对对Maxwell方程组方程组 和Lorentz方程方程进行求解。求解求解Maxwell方程组方程组的关键问题是在含有任意实际结 构形状、任意媒质分布的二维和三维空间内寻找方程的数值计算方法。该算法适 用于各种电磁现象，解决这些问题的方法是直接求解电磁场或电磁通量密度，或 直接引入矢量位、标量位等中间函数，通过有限差分法（FDM）、有限元法（ FEM）、有限积分法（FIM）、边界积分法（BIM）或其

6、他数值计算方法得到各类 电磁问题的数值近似解。 n对于带电粒子在电磁场中的运动，采用粒子模拟法（PIC）能获得精确的结果。由 于粒子模拟技术不再采用近似等效方法，而是根据微波器件的边界条件边界条件和初始条初始条 件件直接求解有源Maxwell方程组和Lorentz运动方程，因此可以获得精确的计算结 果。 n目前粒子模拟技术用于注波互作用计算多出现在大型通用电磁分析软件中，如 MAGIC,MAFIA，ARGUS等。如前表所示，目前国际上已有多个大型电磁分析软件 。除表中列出的大型电磁分析软件外，具有权威性的大信号模拟程序还有美国的 MAGY模拟软件和针对回旋速调管计算的专业软件MAGYKL，但目

7、前这两款程序都 对我国禁运禁运。 国际部分常用大型电磁分析软件 数值法的比较 n1）、有限差分法从电磁场方程的微分方程出发，在整个边界包 含的面积（二维）或体积（三维）区域里划分网格，用差分方程代替 微分方程，形成一个线性方程组。 n2）、有限元法从电磁场的变分原理出发，在整个区域内剖分 ， 形成一个线性方程组。 n3）、边界元法从库仑定律出发，只对边界进行离散，最适于求 解开放性边界问题，而且可以降低方程维数，使问题简单化。 5.1 有限差分法求解轴对称电场 n真空中，不考虑空间电荷，给定封闭边界上电位值，在封闭边 界包围的区域内，电位分布满足： 22 22 1 0 rrrz （5-1） 边

8、值问题边值问题 v存在边界面的电磁场问题（物理）。存在边界面的电磁场问题（物理）。 v根据给定边界条件对边值问题分类：根据给定边界条件对边值问题分类： 第一类边值问题狄里赫利(Dirichlet)问题： 已知电位函数整个边界面上的已知电位函数整个边界面上的 分布值分布值。 第二类边值问题纽曼（ Neumann ）问题： 已知函数在整个边界面上的已知函数在整个边界面上的法向导数法向导数 。 S f S f n 2 2 S f n 第三类边值问题（混合边值问题）： 已知已知一部分一部分边界面上的边界面上的函数值函数值，和，和另一部分另一部分边界面边界面 上函数的上函数的法向导数法向导数。 1 1

9、S f 12 SSS 边值问题边值问题 5.1 有限差分法求解轴对称电场 n 采用一定的网格分割方式离散化场域场域。 n 进行差分离散化处理偏微分方程偏微分方程。用离散的、只含 有限个未知数的差分方程组，来近似代替场域内具有 连续变量的偏微分方程以及边界上的边界条件（也包 括场域内不同媒质分界面上的衔接条件）。 n 结合选定的代数方程组的数值算法数值算法，编制计算机程 序，求解由上面所得对应于待求边值问题的差分方程 组，所得解答即为该边值问题的数值解。 主要求解步骤主要求解步骤: 有限差分法求解场 n一、边界内部处理五点差分 n二、轴上电位处理 n三、边界处理 0 y x 0 2 4 13 D

10、 L h h 注意（关键点）：注意（关键点）： 离散方法： 数值计算中主要的离散方法是泰勒级数法，即用差分来代替微商， 忽略高次项，把微分方程离散成差分方程。 一般采用正方形或矩形网格，网格形状规则简单，宜于求解边界 比较规则的电磁场问题。 5.1.1 轴对称电场 n一、边界内部不等距五点差分： U1 U0 U2 U3 U4 Z0 R0 h1h2 h3 h4 边界内五点差分法 nU0为待求电位，设其余的电位为已知值。将U1、 U2、 U3、U4，各点按泰勒级数展开并精确到二阶偏微分。 22 1 1100 2 22 2 2200 2 22 3 3300 2 22 4 4400 2 ()() 2

11、()()2 2 ()()3 2 ()()4 2 o o o o hUU UUh zz hUU UUh zz hUU UUh rr hUU UUh rr 1 五点差分法 n解得： 2 0120 2 11221212 2 0340 2 33443434 3344 0340 003340434034 222 () ()() 222 () ()() 1 () ()() U UUU h hhh hhh hz U UUU h hhh hhh hr hhhhU UUU rrr h hhr h hhr h h 五点差分法 n将以上方程组代入柱坐标系下的拉普拉斯方程 得五点差分方程 0011223344 C U

12、CUC UC UC U )( 2 211 1 hhh C )( 2 212 2 hhh C )( 2 4330 40 3 hhhr hr C 03 4 0434 2 () rh C r h hh 430 43 4321 0 22 hhr hh hhhh C 对称轴上 n轴上电位处r=0 2 2 0 1 lim r U rr U r 得对称轴上的拉普拉斯方程为： 02 2 2 2 2 z U r U 且，U3=U4，h3=h4，得对称轴上的拉普拉斯差分方程的系数为： 1 112 1 () C h hh 2 212 1 () C h hh 4 0C 3 2 3 2 C h 0 2 123 12 C

13、 h hh h 2 h 3 h 4 U 1 U 0 U 2 U 3 U 4 Z 0 R 0 h 1 0011223344 C UCUC UC UC U z 等距五点差分 n采用等间距的五点差分法，其差分方程的系数为： 12 1CC 0 0 3 5 . 0 r hr C 0 0 4 5 . 0 r hr C 4 0 C 在轴上的差分系数： 12 1CC0 3 C4 4 C6 0 C 0124 64UUUU 0123443 0 4()() 2 h UUUUUUU r 得内部和轴上差分拉普拉斯方程： 边界处理 n边界封闭问题：拉普拉斯方程和泊松方程必须在封闭边界内求解，否则可能得到不 稳定的解，而在

14、实际电子光学系统中，有些边界是敞开的，这些开放式边界处的函 数值是不知道。因此必须人为封闭这些边界。 )( )( )( ao a aoa zz zz UUUU 常用的方法有：线性插值和对数插值 )/ln( )/ln( )( ca c cac rr rr UUUU 线性插值对数插值 差分方程求解 n如图电极系统，求网格m*(n-1)个节点的电位 电极系统电极系统 04)( 2 )( 034 0 4321 UUU r h UUUU 064 0421 UUUU ,1,11,1, 111 (1)(1) 422 i ji ji jijij UUUUU ii 4 6 1 , 1 1 ,1,jijijiji

15、 UUUU 差分方程求解 n迭代法： n通过上述的处理，对于计算区域内的每一个网格点，都可以建立一个差分方程。如 果计算区域内有N个网格点，就形成了N个线性方程构成的大型方程组。解这样的方 程组采用直接法求解是比较困难的，又由于线性方程组的系数构成一个大型的稀疏 阵，因而可以采用迭代法进行电位求解。 n目前，解偏微分方程的迭代法常用如下的四种方法： n一：同步迭代法； n二：逐次超松弛迭代法（SOR）； n三：Chebyshev加速超松弛迭代法（SCA）； n四：交替方向隐式迭代法（ADI）; 其中，第一种方法需要先求出所有网格点上的k次近似值，代入差分方程再求出所 有网格点上的第k+1次近似

16、值，因而计算所需的计算机的内存量大，同时收敛速度 慢。对于最后一种迭代法而言，它比超松弛迭代法收敛快，也更有效，但是其程序 编制较为困难且有时不收敛，因而，在编制程序时可以采用了Chebyshev加速超松 弛迭代法。 同步迭代法 n首先任意给定节点（i，j）上电位的第0次数值作为解的零次近似，然后依次将近 似值带入方程右端，获得点（i，j）上第一次近似解，重复这样的过程，当计算 次数趋向于无穷的时候，就可以无限靠近所考察的微分方程的真实解。 ( ) ,1,11, (1)( )( )( ) 1, 111 (1)(1) 422 k i ji ji jij kkkk ij UUUUU ii ,1,1

17、1, (1)( )( )( ) 1 4 6 i ji ji jij kkkk UUUU )1()( max k n k n UU 2000个节点，误差为1*e-5，迭代次数20003000次 异步迭代法 n也称为赛德尔迭代法，是在计算第k+1次近似值时，位于此点左方和下 方的点一般已经计算出了第k+1次值，则将此近似值代入到方程式右端， 这种计算方式可以比同步迭代法节省一半的迭代次数。 (1) ,1,11, (1)(1)( )( ) 1, 111 (1)(1) 422 k i ji ji jij kkkk ij UUUUU ii ,1,11, (1)(1)( )( ) 1 4 6 i ji j

18、i jij kkkk UUUU 超松弛迭代法 n其中，w是一个介于1和2之间的常数，称为超松弛因子，当w取1的时候就是异步迭代方式， 但是要注意到w的选取也是非常重要的工作，一般会有各种近似公式或者经验公式选取w , (1)( )(1)( ) () i ji ji ji j kkkk UUUU 4)/*Jacobi(1 1 1)(k2 k 2 )Zmax/cos()axR/mcos( Jacobi 可以将迭代次数20003000次减小到100次左右 (1) ,1,11, (1)(1)( )( ) 1, 111 (1)(1) 422 k i ji ji jij kkkk ij UUUUU ii

19、,1,11, (1)(1)( )( ) 1 4 6 i ji ji jij kkkk UUUU 误差分析 n差分法的误差主要有两种： n（1）截断误差，由差分方程代替微分方程所 引起，两者的解之间的差别称为截断误差。截 断误差来源于采用五点差分方程代替微分方程 时舍去的h的三次方以上的项。显然，截断误 差和网格间距大小以及电极结构本身有关。 n （2）迭代误差，决定于迭代计算中误差控制 值，必须选择适当的节点及合适的迭代误差控 制值，以保证工程精度的要求。 计算框图 启动 给定边值 填写场域内 的初始值 叠代次数计数n:=0 n:=n+1 按超松弛迭代法进行 一次迭代，求 ) 1( ),( n

20、 ji 所有内点相邻二次迭代值 的相对误差是否小于W？ 打印迭代次数n，待求 数值解 )j , i ( 停机 是 否 5.1.2 平面对称电场 n平面对称系统：如静电偏转板，不考虑边缘场效应，将 y方向看做无限长。选择xoz面为对称面，则场满足微分 方程： 0 2 2 2 2 z U x U n及相应的边界条件。 5.1.2 平面对称电场 0011223344 C UCUC UC UC U )( 1 211 1 hhh C )( 1 212 2 hhh C )( 1 433 3 hhh C )( 1 434 4 hhh C 4321 0 11 hhhh C n通过划分网格，利用泰勒级数，在区域

21、内任选一个节点O和它相 邻的四个节点，截去高于三阶的项，获得五点差分格式： 5.1.3 尖端发射场 n导体尖端的电荷特别密集，尖端附近的电场特 别强，就会发生尖端放电尖端放电 5.1.3 尖端发射场 n场致发射阴极，在研究这种阴极区的电场分布时，一般 将钨尖看成是球状或者球锥状，采用球坐标。 222 2222222 2111 0 tansin UUUUU rrrrrr n由于系统轴对称，在方向角方向的场是均匀的，故而： 0 1 tan 12 2 2 222 2 U r U rr U rr U 尖端场区网格划分 n利用泰勒级数将相邻四个节点的电位在O点展开，去掉高次 项，得到方程的五点差分格式：

22、 0 0044332211 UCUCUCUCUC )( )1 (2 212 1 2 hhh r h C o )( )1 (2 211 2 1 hhh r h C o )( tan 2 433 4 3 hhh r h C oo 43 34 4321 12 21 0 tan 2)(22 hhr hh hhhhr hh hh C ooo )( tan 2 434 3 4 hhh r h C oo 13 o rh 24 o rh 尖端场区轴上 n当中心节点在z轴上时： 2 2 0 tan 1 lim UU r n并且，该球对称场相对于z轴而言，也可以看成轴对称场，利用近轴 场区特点亦可得出上式。 n则

23、z轴上方程： 0 22 2 2 22 2 U rr U rr U 轴上差分格式 0 00332211 UCUCUCUC 1 2 212 2(1) () o h r C h hh 2 1 112 2(1) () o h r C h hh 3 3 2 4 C h 3 21 0 2 1212 2()24 o hh C h hr h hh o rhhh 432 21 o r h h 1 2 1 简化： 均为函数 注意 n一般电子光学系统里，会有电子注通 过，特别是在强流电子光学系统里， 必须考虑电子注的影响，此时必须求 解泊松方程。 5.2 有限元法 基本思想： 在有限元方法中，场域被分割成许多很小的

24、子区域，通常称为“单元”或 “有限元”。 对所有子区域进行独立的处理和运算，便于对一个整体问题进行局部化处理。 通过选取恰当的尝试函数，使每个单元的计算都变得非常简单，经过对每个 单元重复而简单的计算，再将其结果总和起来，便可以得到用整体矩阵表达的 整个区域的解。 这一整体矩阵又常常是稀疏矩阵，可以更进一步简化和加快求解过程。由于 计算机非常适合于重复性的计算和处理过程，所以整体矩阵的形成过程很容易 使用计算机来实现。 5.2.1 变分原理 n轴对称静电场的能量积分： v dVEJ 2 2 1 v rdrdz z U r U J)()( 2 1 22 积分取极值的条件是变分为积分取极值的条件是

25、变分为0 0J 拉氏方程 n经过变换： ()() v UU JrrUdrdz rrzz 在邻域内电位变分取任何值的情况下，使得能量变分在邻域内电位变分取任何值的情况下，使得能量变分 为为0的条件是：的条件是： ()()0 UU rr rrzz 5.2.1 变分原理 5.2.2 有限元剖分 n应用变分原理，构造函数U（r，z）以使得能量积 分取极值。 n但是在整个域内构造电位函数太过困难，因而采 用剖分方法。 n步骤： n区域划分 n构造函数 n单元分析（能量积分） 一个静电场例子，讲二维有限元：一个静电场例子，讲二维有限元： n二维有限元法 同轴传输线，两个同芯长方形导体之间充满线性介质， 两

26、导体间加有直流电压10v，导体间贮有电荷，传输线的 长度远远大于其截面，可认为电场在传输线各个截面上 的分布都相同，只需求解电场在某个截面的分布。 场域剖分场域剖分 n二维有限元法 原则上讲，二维有限元可以取为各种多边形， 如三角形、四边形等等。与其它多边形相比， 三角形具有以下两个优点： (1)描述二维三角形的多项式有3项，该数 目与三角形的顶点数以及节点上未知量的个 数恰好相同，因而使得多项式形函数的利用 率最高 。 (2)三角形形状简单，能十分便利地表示复 杂的几何结构。 把两个要求解的量联系 起来，有限元中令待定 系数就是节点电位，当 然尝试函数要重新确定 bfK 任意单元内节点尝试函

27、数的选取：任意单元内节点尝试函数的选取： 二维有限元法 任意三角形单元由节点i，j，k构成，每个节点的尝试函数（三角形平面） 选择的规律一样. 注意：注意： 不能将一个三角形的顶点取为另一个 相邻三角形边的内点 每个单元和节点都要按逆时针编号。 节点尝试函数的表达式：节点尝试函数的表达式： n3. 二维有限元法 代入平面方程确定平面参数：（用矩阵形式作规范化求解） rz 平面方程(线性插值)： ( ,)( ,)( ,) iijjkk i r zj r zk r z 平面过三个点： ； iii jjj kkk rzi rzj rzk 1 1 1 iii jjj kkk rz rz rz 应用克莱

28、姆法则求节点尝试函数（平面方程）表达式：应用克莱姆法则求节点尝试函数（平面方程）表达式： n二维有限元法 S为三角形单元面积，为使其为正，i j k要 逆时针编号。平面方程系数有严格的规律， 且只与节点坐标有关 二元有限元 n令： 1 () jkkj ar zr z)( 2kiik zrzra)( 3ijji zrzra 1jk bzz ik zzb 2 ji zzb 3 jk rrc 1 ij rrc 3 ki rrc 2 123 () 2 ijk m aaa S 123 () 2 ijk m bbb S 123 () 2 ijk m ccc S n得 123123123 3 1 ()()(

29、) ( , ) 2 1 () 2 ijijijkkk m pppp pm aaabbbrcccz r z S ab rc z S n我们构造了电位函数，它由节点参数（坐标， 电位值）来描述，这个步骤称为构造函数 n二维有限元法 在小单元内进行能量积分： 22 1 ()() ) 2 m v Jrdrdz rz 22 123123 2 1 ()() 2 4 mm ijkijk s r bbbcccdrdz S 3 2 ,1 1 () 2 4 mm pqpqpq p q s r b bc cdrdz S n令 m m s m rdrdz S 2 4 1 则： )( 2 1 3 1, qp qp qp

30、qpmm ccbbJ 单元数共有单元数共有N个，则总的能量积分：个，则总的能量积分： 3 1,1 1 () 2 N mpqpqpq mp q Jb bc c n单元分析时，我们仅用了m单元 的节点i，j，k，但是m单元是和 四周有联系的，每个单元和相邻 单元有两个节点共有，而且单元 节点编号是按顺序排列。 n对于每个节点，他们的b，c系数 不同，但是电位积是相同的，因 此当某节点与其它几个单元相联 系，则有几组系数相加作为它的 系数。 n故： ji n ji ji J 1, , 2 1 注 意 n此时，我们将能量积分用多元函数来描述， 能量积分取极值的问题转化成多元函数取极 值的问题。 n多元

31、函数取极值即是对多元函数微分等于0： 0 i J 0 1, , N ji jji 同轴矩形电场分布同轴矩形电场分布 在大力推广在大力推广CAD技术的今天，从自行车到航天飞机，技术的今天，从自行车到航天飞机， 所有的设计制造都离不开有限元分析计算，所有的设计制造都离不开有限元分析计算，FEM在工在工 程设计和分析中将得到越来越广泛的重视。程设计和分析中将得到越来越广泛的重视。 国际上早在国际上早在20世纪世纪50年代末、年代末、60年代初就投入大量年代初就投入大量 的人力和物力开发具有强大功能的有限元分析程序。的人力和物力开发具有强大功能的有限元分析程序。 其中最为著名的是由美国国家宇航局（其中

32、最为著名的是由美国国家宇航局（NASA）在）在 1965年委托美国计算科学公司和贝尔航空系统公司年委托美国计算科学公司和贝尔航空系统公司 开发的开发的NASTRAN有限元分析系统。该系统发展至今有限元分析系统。该系统发展至今 已有几十个版本，是目前世界上规模最大、功能最强已有几十个版本，是目前世界上规模最大、功能最强 的有限元分析系统。的有限元分析系统。 v目前，世界各地的研究机构和大学发展了一批规模目前，世界各地的研究机构和大学发展了一批规模 较小但使用灵活、价格较低的专用或通用有限元分析较小但使用灵活、价格较低的专用或通用有限元分析 软件软件: v主要有德国的主要有德国的ASKA、英国的、

33、英国的PAFEC、法国的、法国的 SYSTUS、美国的、美国的ABQUS、ADINA、ANSYS、 BERSAFE、BOSOR、COSMOS、ELAS、MARC和和 STARDYNE等公司的产品。等公司的产品。 随着数值分析方法的逐步完善，尤其是计算机运随着数值分析方法的逐步完善，尤其是计算机运 算速度的飞速发展，整个计算系统用于求解运算的时算速度的飞速发展，整个计算系统用于求解运算的时 间越来越少，而数据准备和运算结果的表现问题却日间越来越少，而数据准备和运算结果的表现问题却日 益突出。益突出。 在现在的工程工作站上，求解一个包含在现在的工程工作站上，求解一个包含10万个方万个方 程的有限元

34、模型只需要用几十分钟。工程师在分析计程的有限元模型只需要用几十分钟。工程师在分析计 算一个工程问题时有算一个工程问题时有80%以上的精力都花在数据准备以上的精力都花在数据准备 和结果分析上。和结果分析上。 增强可视化的前置建模和后置数据处理功能增强可视化的前置建模和后置数据处理功能 v增强可视化的前置建模和后置数据处理功能增强可视化的前置建模和后置数据处理功能 目前几乎所有的商业化有限元程序系统都有功目前几乎所有的商业化有限元程序系统都有功 能很强的前置建模和后置数据处理模块。使用能很强的前置建模和后置数据处理模块。使用 户能以可视图形方式直观快速地进行网格自动户能以可视图形方式直观快速地进行

35、网格自动 划分，生成有限元分析所需数据，并按要求将划分，生成有限元分析所需数据，并按要求将 大量的计算结果整理成变形图、等值分布云图，大量的计算结果整理成变形图、等值分布云图， 便于极值搜索和所需数据的列表输出。便于极值搜索和所需数据的列表输出。 v与与CAD软件的无缝集成软件的无缝集成 当今有限元分析系统的另一个特点是与通用当今有限元分析系统的另一个特点是与通用CAD软件的集软件的集 成使用，即：在用成使用，即：在用CAD软件完成部件和零件的造型设计后，软件完成部件和零件的造型设计后， 自动生成有限元网格并进行计算，如果分析的结果不符合自动生成有限元网格并进行计算，如果分析的结果不符合 设计

36、要求则重新进行造型和计算，直到满意为止，从而极设计要求则重新进行造型和计算，直到满意为止，从而极 大地提高了设计水平和效率。大地提高了设计水平和效率。 当今所有的商业化有限元系统商都开发了和著名的当今所有的商业化有限元系统商都开发了和著名的CAD软软 件（例如件（例如Pro/ENGINEER、Unigraphics、SolidEdge、 SolidWorks、IDEAS、Bentley和和AutoCAD等）的接口。等）的接口。 v复合材料加工传热媒质问题复合材料加工传热媒质问题 v 不同加热灯丝不同加热灯丝位置位置情况下阴极的温度分布云图情况下阴极的温度分布云图 微观尺度材料设计微观尺度材料设

37、计 有限元方法 v半导体芯片温度场的数值仿真半导体芯片温度场的数值仿真 宏观尺度材料设计宏观尺度材料设计 有限元方法 v水轮机叶轮的受力分析模拟水轮机叶轮的受力分析模拟 5.3 电荷密度法求解轴对称电场 n许多物理问题可通过不同的途径归结为不同的数学模型许多物理问题可通过不同的途径归结为不同的数学模型 （大多数没有解析解）：（大多数没有解析解）： n a 偏微分方程的边值问题有限差分法偏微分方程的边值问题有限差分法 n b 区域上的变分问题有限元法区域上的变分问题有限元法 n C 边界上的积分问题边界元法边界上的积分问题边界元法 方法对比 n有限差分与有限元法异同点： n原理上： 有限差分法从

38、电磁场的微分方程出发，在整个区域里划分网格， 用差分方程代替微分方程，形成线性方程组 有限元法从电磁场的变分原理出发，在整个区域里剖分，形成 线性方程组 n共同的弱点： 1。都需要一个较为规则的封闭边界。 2。不论感兴趣的区域大小，都必须在整个区域内求解电场分布。 n边界元法从库仑定律出发，只对边界进行离散，最适 于求解开放性边界问题，而且可以降低方程维数，使问题简单化。 边界元法 nboundary element method又称边界积分方程-边界元法。 n边界元法是一种继有限元法之后发展起来的新型数值方法,与有限 元法在连续体域内划分单元的基本思想不同,边界元法仅在定义域 的边界上划分单

39、元,用满足控制函数去逼近边界条件. n它以定义在边界上的边界积分方程为控制方程，通过对边界分元 插值离散，化为代数方程组求解。 n它与基于偏微分方程的区域解法相比，由于降低了问题的维数， 而显著降低了自由度数，边界的离散也比区域的离散方便得多， 可用较简单的单元准确地模拟边界形状，最终得到阶数较低的线 性代数方程组。 边界元法 n又由于它利用微分算子的解析的基本解作为边界积分方程的核 函数 ，而具有解析与数值相结合的特点，通常具有较高的精度。 n特别是对于边界变量变化梯度较大的问题 ，如应力集中问题 ， 或边界变量出现奇异性的裂纹问题，边界元法被公认为比有限 元法更加精确高效。 n由于边界元法

40、所利用的微分算子基本解能自动满足无限远处的 条件，因而边界元法特别便于处理无限域以及半无限域问题。 边界元法 n边界元法的基础：边界元法是基于控制微分方程的基本解来建立 相应的边界积分方程，再结合边界的剖分而得到的离散算式。 nJaswon和Symm于1963年用间接边界元法求解了位势问题；Rizzo于 1967年用直接边界元法求解了二维线弹性问题；Cruse于1969年将 此法推广到三维弹性力学问题。1978年，Brebbia用加权余量法推 导出了边界积分方程，他指出加权余量法是最普遍的数值方法， 如果以Kelvin解作为加权函数，从加权余量法中导出的将是边界 积分方程边界元法，从而初步形成

41、了边界元法的理论体系， 标志着边界元法进入系统性研究时期。 边界元法的发展 n经过近40年的研究和发展，边界元法已经成为一种精确高效的工程 数值分析方法。在数学方面，不仅在一定程度上克服了由于积分奇异性 造成的困难，同时又对收敛性、误差分析以及各种不同的边界元法形式 进行了统一的数学分析，为边界元法的可行性和可靠性提供了理论基础。 n在方法与应用方面，现在，边界元法已应用到工程和科学的很多领 域，对线性问题，边界元法的应用已经规范化；对非线性问题，其方法 亦趋于成熟。 n在软件应用方面，边界元法应用软件已由原来的解决单一问题的计 算程序向具有前后处理功能、可以解决多种问题的边界元法程序包发展。

42、 n我国约在1978年开始进行边界元法的研究，目前，我国的学者在求 解各种问题的边界元法的研究方面做了很多的工作，并且发展了相应的 计算软件，有些已经应用于工程实际问题，并收到了良好的效果。 n边界元法是将边界元法是将区域内微分方程区域内微分方程 n通过积分定理变为通过积分定理变为边界上的积分方程边界上的积分方程 n再将积分方程在边界上离散为再将积分方程在边界上离散为代数方程代数方程。 电荷密度法 n从库仑定律出发，最适于求解开发性 边界问题 n对边界进行离散化处理 n电荷密度法也是一种边界元法 点电荷点电荷 可以简化为点电荷的条件可以简化为点电荷的条件: Q 1 r d dr 观察点 P 库

43、仑定律库仑定律 库仑定律：库仑定律：在真空中，两个静止点电荷之间相互在真空中，两个静止点电荷之间相互 作用力与这两个点电荷的电荷量作用力与这两个点电荷的电荷量q1和和q2的乘积成的乘积成 正比，而与这两个点电荷之间的距离正比，而与这两个点电荷之间的距离r12（或（或r21） 的平方成反比，作用力的方向沿着这两个点电荷的平方成反比，作用力的方向沿着这两个点电荷 的连线，的连线，同电相斥同电相斥，异电相异电相吸。吸。 i i i i r r qq kF 03 0 0 1785年，法国库仑（年，法国库仑（C.A.Coulomb) 适用于点电荷适用于点电荷 叠加性叠加性 q0 q1 q2 r02 F2

44、 r01 F1 F 库仑库仑 12 12 12 3 21 r r qq kF n将无穷远处看作电位零点，则点电荷在空间 中任意一点产生的电位为： 4 1 )( rr q r 库仑定律不仅对点电荷适用，对线电荷，面电荷和体电荷同样适用。 l l rr dl r 4 1 )( 4 1 )( s l rr ds r 4 1 )( rr d r l n在空间中，同时存在N个充满电荷的源，则他们在空中任意一点P 产生的电位满足电位叠加定理： N i i rr 1 )()( 在有限的区域内，忽略自由电荷，则区域内的电场是由加以一定电 压的电极形成的。整个区域的电位分布由电极表面电荷产生 如何求电极上的表面

45、电荷分布？ 充分利用电极电位已知的条件，将场点设在电极面上，则每个电极 的电位都满足： 4 1 )( s l rr ds r 4 1 )( s l rr ds r 本式表明本式表明: 电极表面的电荷分布确定了各个电极的电位，还确定了区域内的电位分布； 同时也表示为了使电极表面有这样的电荷分布，必须使各个电极的电位为给 定值。 电荷密度法步骤：电荷密度法步骤： 利用上式左端电位已知，采用数值法求解方程，求出电极表面电荷密度分布; 然后将电荷密度代入上式，求解空间任意一点的电位，则给定区域内的电位 分布就唯一确定了 5.3.2积分方程离散化 n上式的关键就是确定空间电荷分布，解析法难以完成任务，

46、而采用数值法，就必须进行离散化处理。 n将电荷存在的区域分成小区域，当区域足够小的时候，电荷 密度分布可以解析解，如最简单的认为：电荷分布是均匀的。 1 4 4 1 )( s N i i s l rr ds rr ds r N i ipip Cr 1 , )( Ci,p与小区域形状有关，与P点到小区域距离有关，但与 小区域的电位和电荷无关 简写为： 全部电荷在P点产生的电位： n将电极表面分成许多小区域，如N块，就可以 写成N个线性方程： N i ijij Cr 1 , )( 1,2,3,.,jN 5.3.3 用电荷密度法求轴对称电场 n真空中2个加上电压的半径不等，长短不 一，忽略厚度的圆筒

47、电极组成的轴对称电 极系统，求其内的电场分布。 2 , 0222 1 4 (coscos)(sinsin)() ii i ii zD siii i j zD iijjiijjij R ddz RRRRzz 第i个环带在P点产生的电位为： 上式中Di是第i个环带的半宽度。令 2 i 22 2 )()( 4 jiji ji zzRR RR k n故： 2 0 22 22 , ) 2 (sin1 ( )()( 4 j Dz Dz jiji ii s ji k d zzRR dzRii ii i i sji C , 其中，Cij表示第i个环上所有的面电荷在第j个环上P点处产生电位 的系数。 2 0 2

48、2 ) 2 (sin1 ( i k d 是第一类椭圆积分 当 i 从 1 到N 表示 N 个小环在第 j 个环上 P 点处产生的电位 ； 当 j 从 1 到N 表示第 i 个环上所有的面电荷在每个环上产生的电位； 合起来就形成了一个方程组或一个矩阵。 实质上只要求得Cij系数 , 就可以求出我们所要的电位分布了。 n检视上面的公式，需要注意： 当半径相等，角度相同，Z向坐标也一致的时 候（场点和源点重合），出现被积函数分母为 零，积分值为无穷大的现象。而实际上，该点 的电位和表面电荷都是有限的。这样的点称为 奇异点。 2 , 022 22 4 ()() (1sin () 2 ii i ii z

49、D s ii i j zD j ijij R dzd RRzz k 边界元法的主要缺点 n边界元法的主要缺点是它的应用范围以存在相应微分算子的基本 解为前提，对于非均匀介质等问题难以应用，故其适用范围远不 如有限元法广泛，而且通常由它建立的求解代数方程组的系数阵 是非对称满阵，对解题规模产生较大限制。对一般的非线性问题， 由于在方程中会出现域内积分项，从而部分抵消了边界元法只要 离散边界的优点。 n即：边界元法与有限元相比具有单元的未知数少,数据准备简单等 优点.但用边界元法解非线性问题时,遇到同非线性项相对应的区 域积分,这种积分在奇异点附近有强烈的奇异性,使求解遇到困难。 5.3.4 奇异

50、点的处理 n第一类：被积函数出现断点（eg: 当场点和源点重合，被积函数值 无穷大，而物理上说，该点的电位和电荷密度应当是有限的。） 处理方法： 1. 将被积函数区域分成奇异区和非奇异区，通过改变积分参数消除 奇异性； 2. 将被积函数作数学处理，使它变成可积函数； 3. 采用高斯型积分公式。 n另一类：奇异点出现在电极的两端，即曲率半径小的地方。由电磁 学可知：在一个处于电场中的导体，内部的E0，电位是常数，电 荷分布在其表面，曲率半径大的地方电荷分布少，曲率半径小的地 方，电荷分布多 - - -极端表现：尖端放电 5.3.4 奇异点的处理 n在上节分析中，将电极表面分成一个个小区域，每个小

51、区域的电荷密 度被认为是均匀的，可以看作是一个常数，这个数值往往是这个区域 电荷密度的平均值，这样的假设和许多实际情况很不相符。为了更精 确求出电极上电荷分布情况，对这一类奇异点也作特殊考虑。 n例如：采用不均匀划分子区域法对数形式划分，或者采用二进制 形式划分，设每个子区域上的电荷密度是常数，电极最边缘的子区域 较小，且数目较多，可以使电荷密度分布接近实际情况。 n另一方法是：均匀划分区域，但是区域内的电荷密度看成连续变化的 函数如： n进一步提高精度，可以采用二次插值函数，甚至采用切比雪夫多项式。 taatf os1 )( 5.3.5 系数矩阵的解法 n电荷密度法里，求解电极上的电荷分布和

52、空间中任意一点 的电位，最后都归结到求系数矩阵C和它的逆矩阵C-1。 n即需要确定C中的每一个元素Cij，需要进行大量的积分计 算，对系数矩阵中的积分常采用以下几种方法: 1.辛普生积分 2. 高斯积分法 3. 第一类、第二类完全椭圆积分的近似计算 辛普生积分 n将定积分区间a,b划分成2m个等分，得到2m个小 区间，区间长度为： m ab h 2 )( hiaxihax ii ) 12(2 122 在每对小区间上，用通过三点（x2i，x2i+1，x2i+2）的二 次抛物线来近似被积函数 （线性插值仅仅利用了两、三个节点的数据信息，因此逼近度自然不高） b a m i m i ii xfxfb

53、faf h dxxf )(4)(2)()( 3 )( 1 1 1 1 122 2 高斯积分法 ii t abab x 22 b a m i ii xfdxxf 1 )()( 上式中ti是勒让德多项式Lm（t）的第i个零点： 22 )()1 ( 2 imi i tLt 3.第一类、第二类完全椭圆积分的近似计算 n第一类椭圆积分： .ln.)( 4 410 4 410 mbmbbmmamaakK 2/ 022 )sin1 ( )( k d kK n第二类椭圆积分： 2/ 0 22 )sin1 ()( dkkE 两者的近似计算式： .ln.1 )( 4 41 4 41 mbmbmmamakE 2 1

54、 km 5.4 电子轨迹的数值求解方法 n电子在场中运动描述：1. 电子轨迹方程 2. 电子运动方程 n只需给出场的分布函数或者轴上电磁位的分布函数和 初始条件后，就可以求出电子运动解。 n困难在于：1. 场的分布常常不能用简单的解析式给出， 甚至无法用函数形式给出。 2. 即使给出了场的分布解析式，也无法得 出解的解析式。 给出离散值（通过数值计算或者实验） n电子轨迹方程的数学描述是一个二阶非线性非齐次微 分方程 n电子运动方程的数学描述是一个二阶齐次微分方程 n求解方法：1、图解法，橡皮膜法，自动轨迹仪等 2、数值法 一、高阶泰勒法一、高阶泰勒法 .),(t) (1) )( ),( 足够

55、光滑及的解ytfy ay btaytf dt dy 假设初值问题 ., )!1( )( ! )( ! 2 )( )()()( ,)( 1 1 )1()( 2 1 1 iii n i n n i n i iii ii tt h n y h n ty h ty htytyty ntty 其中 得阶泰勒展开处作在将 5.4.1 龙格库塔法 2 1 (1)( ) 1 ( )( , ( ),( , ) ( , ( ) () ( )( , ( ) 2! ( , ( )( , ( ) !(1)! 再将代入 得由已知 iii ii iiii nn nniiii dy y tf t y tf t yatb dt ft y t y ty tf t y thh ft y tft y t hh nn 称上式为n阶泰勒法 0 (1) 2 1 ( ,)( ,) ( ,) 2! n niiii iiii h y ft yft y yyf t y hhh n 当 充分小时 二、二阶龙格二、二阶龙格库塔法库塔法 ).()()( (3) ),( ),( )( , (2) )( ),( 3 1 12 1 22111 0 21 hOtyty hkyhtfk ytfk kckchyy