年7月电子与信息学报

第32卷第7期电子与信息学报Vol.32No.72010年7月JournalofElectronics&InformationTechnologyJul.2010三维电子光学模拟器二阶有限元方法研究胡权黄桃杨中海李斌李建清(电子科技大学大功率微波电真空器件技术国防科技重点实验室成都610054)摘要：该文介绍了二阶有限元方法的基本理论，同时将其应用到3维微波管模拟器套装(MicrowaveTubeSimulatorSuite,MTSS)中的电子光学模拟器(ElectronOpticsSimulator,EOS)。比较了EOS分别采用一阶和二阶有限元方法的计算结果，得到了二阶具有更好的收敛性以及更快的收敛速度的结果。关键词：电子光学系统；电子枪；有限元；CAD中图分类号：O463.1;TN12文献标识码：A文章编号：1009-5896(2010)07-1726-05DOI:10.3724/SP.J.1146.2009.010493DSecond-OrderFiniteElementMethodforElectronOpticsSimulatorHuQuanHuangTaoYangZhong-haiLiBinLiJian-qing(VacuumElectronicsNationalLab,UniversityofElectronicScienceandTechnologyofChina,Chengdu610054,China)Abstract:Inthispaper,3Dsecond-orderfiniteelementmethodforelectronopticsystemispresented.AnditisimplementedinElectronOpticsSimulator(EOS).Theresultsarecomparedbetweenfirst-orderfiniteelementmethodandsecond-orderfiniteelementmethod.Second-orderfiniteelementmethodhasbetterconvergenceandfasterconvergencerate.Keywords:Electronopticsystem;Electrongun;Finiteelementmethod;CAD1引言随着计算机软硬件技术的不断提高，计算机模拟技术得到了飞速发展。可用于模拟电子光学系统的软件也由2维发展到了3维，一些专用和商用软件逐步出现，如MICHELLE[1,2]，OPERA[3]，Magic[4]，CST[5]等。电子科技大学研发的电子光学模拟器(ElectronOpticsSimulator,EOS)[6]，是3维微波管模拟器套装(MicrowaveTubeSimulatorSuite,MTSS)[7]的重要组成部分。EOS采用3维实体建模、3维有限元网格自动划分、静态轨迹迭代算法求解电子运动状态。可用于模拟设计轴对称电子枪、栅控电子枪、多注电子枪、带状注电子枪、轴对称及非轴对称多级降压收集极等。EOS采用的是常应变四面体单元，即一阶有限元方法。实际工程中的电场往往随坐标变化，但常应变四面体单元中的电场分量都是常量，难以适应急剧变化的电场。为了保证计算精度，必须采用密集的计算网格，这样节点数量将很多，方程组十分庞大。若采用高次插值函数，单元中的场是变化的，2009-07-29收到，2009-12-21改回国家自然科学基金(10476004,60601004)和国防科技重点实验室基金资助课题通信作者：胡权[e-mail]则可用较少的单元、较少的自由度而得到较高的计算精度，从而可降低方程组的规模。为此，我们在EOS中引入了二阶插值函数[8,9]。2有限元理论2.1体积坐标如图1所示四面体单元ijmp中，任意一点Q的位置可由下列4个比值来确定。图1四面体单元体积坐标/,/,/,/iijjmmppLVVLVVLVVLVV====(1)式中V为四面体ijmp的体积，Vi为四面体Qjmp的体积，Vj为四面体Qmpi的体积，Vm为四面体Qijp的体积，Vp为四面体Qijm的体积。第7期胡权等：三维电子光学模拟器二阶有限元方法研究1727111161iiijjjmmmpppxyzxyzVxyzxyz=(2)这4个比值称为Q点的体积坐标。由于Vi+Vj+Vm+Vp=V，因此Li+Lj+Lm+Lp=1。根据几何关系，直角坐标与体积坐标之间应符合下列关系：116iiiiijjjjjmmmmmpppppLabcdLabcdxyLabcdVzLabcd⎧⎫⎧⎫⎪⎪⎪⎪⎧⎫⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪−⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪=⋅⎨⎬⎨⎬⎨⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪−⎪⎪⎩⎭⎪⎪⎪⎪⎪⎪⎪⎪⎩⎭⎩⎭(3)式中ai，bi，ci，di(i，j，m，p)的计算公式如下1,1,1111,111jjjjjimmmimmpppppjjjjimmimmppppyzxyzaxyzbyzxyzyzxzxycxzdxyxzxy==−==−(4)为了使四面体的体积不致成负值，单元节点的编号i，j，m，p必须依照一定的顺序。在右手坐标系中，当按照i-j-m的方向旋转时，右手螺旋应向p的方向前进。求体积坐标的幂函数在四面体单元上的积分时，可应用式(5)!!!!ddd6(3)!abcdijmpVabcdLLLLxyzVabcd=++++∫∫∫(5)体积坐标的幂函数对直角坐标求逆时，可利用式(6)：161616ijmpijmpijmpijmpijmpijmpbbbbxVLLLLccccyVLLLLddddzVLLLL⎫⎛⎞⎪∂∂∂∂∂⎟⎪⎜⎟⎜⎪=−+−⎟⎜⎪⎟⎜⎟⎪∂∂∂∂∂⎝⎠⎪⎪⎪⎪⎛⎞∂∂∂∂∂⎪⎟⎜⎪⎟⎜=−+−⎬⎟⎜⎟⎪⎜⎟∂∂∂∂∂⎝⎠⎪⎪⎪⎪⎛⎞⎪∂∂∂∂∂⎟⎜⎪⎟⎜=−+−⎪⎟⎜⎟⎪⎜⎟∂∂∂∂∂⎪⎝⎠⎪⎭(6)2.2插值函数2.2.1一阶插值函数图2给出了一个四面体单元e的结构及结点编号。对于这样一个四面体单元，未知量Φ能够近似为(,,)eeeeexyzabxcydzΦ=+++(7)对应的一阶插值函数为()()1,,,,,,6eeeeekkkkkeNxyzabxcydzkijmpV=+++=(8)图210结点四面体单元且具有如下性质：()1,,,0,ektttktktNxyzktδ⎧=⎪⎪⎪==⎨⎪≠⎪⎪⎩(9)因此四面体e内任意一点(x，y，z)的电位()(),,,,peeekkkixyzNxyzΦΦ==∑(10)即四面体内任意一点的电位是四面体4个顶点电位的插值。对于常应变四面体单元，所采用的一阶插值函数可用体积坐标表示如下：,,,,kkNLkijmp==(11)2.2.2二阶插值函数一阶4结点四面体单元内的电场强度是常量，这种单元的精度不高。为了提高单元精度，可增加单元的自由度个数。除四面体4个顶点i，j，m，p外，在6条棱边的中点各增加一个结点，构成10结点2次四面体单元，如图2所示。每个节点有3个自由度，2次四面体单元有30个自由度，未知量Φ可近似为2次多项式12345678910222(,,)eeeeeeeeeeexyzaaxayazaxayazaxyayzazxΦ=+++++++++(12)与平面2次三角形单元类似，不难证明，2次四面体单元是完备单元。对于上式的待定系数a1~a10，若在整体坐标系下由节点未知量解出，其计算十分冗繁。为此，2次单元的未知量Φ可在体积坐标系下由插值函数表示：616iijjmmppkkkkkkiNNNNNNΦΦΦΦΦΦΦ===++++=∑∑(13)各节点处的插值函数可由体积坐标表示为14(21),,,,4,1,2,3;,,4,4,5,6;,,,iiijmipNLLijmpNLLijmNLLijmp⎫⎪=−⎪⎪⎪⎪=⎬⎪⎪⎪=⎪⎪⎭(14)与一阶插值函数一样，二阶插值函数同样具有式(9)的性质。2.3有限元采用有限元法计算空间的静电场分布，实际上1728电子与信息学报第32卷处理的是如下边值问题：xyzfxxyyzzΦΦΦααα⎛⎞⎛⎞⎛⎞∂∂∂∂∂∂⎟⎜⎟⎟⎜⎜−−⎟−=⎟⎟⎜⎜⎜⎟⎟⎟⎜⎜⎜⎟⎝⎠⎝⎠∂∂∂∂∂∂⎝⎠(15)边界条件为在上1,pSΦ=(16),xyzxyznqxyzΦΦΦαααγΦ⎛⎞∂∂∂⎟⎜++⎟⋅+=⎜⎟⎜⎟∂∂∂⎝⎠在上2S(17)上面的边值问题等价于如下变分问题：()在上10,FpSδΦΦ⎫=⎪⎪⎪⎬⎪=⎪⎪⎭(18)其中()222221d2dd(19)2xyzVSVFVxyzqSfVΦΦΦΦαααγΦΦΦ⎡⎤⎛⎞⎛⎞⎛⎞∂∂∂⎟⎢⎥⎜⎟⎟⎜⎜=+⎟+⎟⎟⎜⎜⎜⎢⎥⎟⎟⎟⎜⎜⎜⎟⎝⎠⎝⎠∂∂∂⎝⎠⎢⎥⎣⎦⎛⎞⎟⎜+−−⎟⎜⎟⎜⎝⎠∫∫∫∫∫∫∫∫整体的变分等于各个四面体单元的变分之和，()()1MeeeFFΦΦ==∑(20)e为网格单元编号，M为网格单元的总数。当γ=q=0时，()2221++d2d(21)eeeeeeexyzVeVFVxyzfVΦΦΦΦαααΦ⎡⎤⎛⎞⎛⎞⎛⎞∂∂∂⎟⎢⎥⎟⎟⎜⎜⎜=⎟⎟⎟⎜⎜⎜⎢⎥⎟⎟⎟⎜⎜⎜⎟⎝⎠⎝⎠∂∂∂⎝⎠⎢⎥⎣⎦−∫∫∫∫∫∫式(21)对eiΦ求偏导，得6/()++dd,,,,,1,,6eeeejeeiixVjieeeejjiiyzeiVNNFxxNNNNVxyxzfNVjijmpΦααα=⎡∂∂⎢∂∂=⎢∂∂⎣⎤∂∂∂∂⎥⎥∂∂∂∂⎦−=∑∫∫∫∫∫∫(22)式(22)可写成矩阵形式eeeeeFΦ∂=−∂KBΦ(23)当xα，yα，zα都等于1时，矩阵Ke和Be中的元素可表示为++d,,,,,,1,,6(24)elqeeeeeeqqqelllVNNNNNNKVxxyyzzlqijmp⎛⎞∂∂∂∂∂∂⎟⎜⎟=⎜⎟⎜⎟⎜∂∂∂∂∂∂⎝⎠=∫∫∫deellVBfNV=∫∫∫(25)一旦用有限元法求出了电位Φ，即可由下式求出电场强度EΦ=−∇(26)2.3.1一阶有限元对于一阶有限元，矩阵Ke和Be中的元素为1(),,,,,36eeeeeeklklklkleKbbccddklijmpV=++=(27),,,,4eeekVBfkijmp==(28)四面体单元e内任意一点(x,y,z)的电场强度为()16peeeeekkkkekiEbxcydzVΦ==−++∑(29)所以采用一阶有限元方法求解得到的同一四面体内的电场强度是不变的。2.3.2二阶有限元对于二阶有限元，矩阵Ke和Be中的元素相对比较复杂，可表示为()2221,,,,60eeeellllleKbcdlijmpV=++=(30a)()1180eeeeeeeijijijijeKbbccddV=++(30b)()()()11180eeeeeeeeeeiijmijmijmeKbbbccccccV⎡⎤=−+−+−⎢⎥⎣⎦(30c)()()()22221122122452(30d)eeeeeeeeemmjjmmjjeeeeemmjjKbbbbccccVcccc⎡=−++−+⎢⎣⎤+−+⎥⎦()()()21222124522eeeeeeeemimijjmeeeeeeeemimijjmeeeeeeemimijjmKbbbbbbbVcccccccddddddd⎡=+−−⎢⎣++−−⎤++−−⎥⎦(30e)1,,,,20eeelBVflijmp=−=(31a)1,1,,65eeelBVfl==(31b)四面体单元e内任意一点(x,y,z)的X向电场强度为161[(41)(41)64()+4()]xieeeeeeeiipppeeeeeeeeeeejmmjpmmpELbLbVLbLbLbLbΦΦΦΦ=−−++−+−+−(32)将上式中的b分别改为c，d便得到Ey，Ez。四面体单元e内的电场强度是体积坐标L的函数，而体积坐标L是x，y，z的函数，所以电场强度也是x，y，z的函数，即同一四面体内的电场强度是随坐标变化的。3模拟计算本文对同一把电子枪分别用