下载文档原格式
第39卷第2期原子能科学技术Vol.39,No.2 2005年3月Atomic Energy Science and TechnologyMar.2005100Me V 紧凑型回旋加速器主磁铁的三维磁场有限元分析技术钟俊晴,张天爵,杨建俊,储诚节(中国原子能科学研究院核技术应用研究所,北京 102413)摘要:100MeV 紧凑型回旋加速器主磁铁的几何结构十分复杂,但为了形成加速器束流动力学所要求的磁场分布,本文对初步设计的磁铁进行必要的简化。
综合采用各种适当的三维有限元网格剖分技术,对该磁铁的磁场进行数值分析,计算精度满足加速器物理设计的要求。
关键词:紧凑型回旋加速器;有限元法;剖分;双标量位法中图分类号:TL542 文献标识码:A 文章编号:100026931(2005)022*******3D Finite Element Analysis Method for Main Magnet Design of 100Me V Compact CyclotronZHON G J un 2qing ,ZHAN G Tian 2jue ,YAN G Jian 2jun ,C HU Cheng 2jie(China I nstitute of A tomic Energ y ,P.O.B ox 27523,B ei j ing 102413,China )Abstract : To get t he magnetic field distribution for beam dynamics calculation in t he 100MeV co mpact cyclot ron feat ured wit h a complex geomet ric magnet ,we simplified t he magnet struct ure and used various 3D finite element meshing met hods in t he initial design stage.After t he result s comparison by different meshing met hods ,we analysed numerically t he t hree dimentional magnetic field and found t hat t his 3D finite element calculation p rocess could meet t he requirement of t he p hysics design for 100MeV com 2pact cyclot ron.K ey w ords :compact cyclot ron ;finite element met hod ;meshing ;double scalar potential收稿日期:2004211225基金项目:国家自然科学基金资助项目(10125218)作者简介:钟俊晴(1978—),男,江西瑞金人,实习研究员,加速器专业 串列加速器升级工程主要由以下几部分组成:100MeV 回旋加速器、在线同位素分离器、超导直线增能器以及现有HI 213串列加速器注入器的改造。
5.3 三维静磁场的有限元分析5.3.1 边值问题以标量磁位m ϕ表示的无源区磁场的边值问题与电位的拉普拉斯边值问题的数学表达形式完全一样，可以如前节所述的有限元分析。
在此，考虑有电流存在以矢量磁位A 作为待求变量的有限元分析。
设在线性媒质中，磁场满足的边界条件：边界1S 面上有0A A =，在边界S 2面上取某种形式的对称面作为第二类齐次边界，在该面上磁场强度H 的切向分量为零：()0=⨯⨯∇=⨯n m n e A e H γ，有如下边值问题：()()⎪⎩⎪⎨⎧∈=⨯⨯∇∈=∈=⨯∇⨯∇210s s V n mm e A A A J A γγ5.3.2 场域剖分与插值对于求解场域V ，根据其形状和场的定性分布，选择合适的单元(例如四面体单元)，进行场域剖分，得到0Z 个单元、0N 个节点。
在单元e 内，对位函数A 进行插值。
若采用四面体单元：∑==41j je j N A A ~式中ej N 是单元形状函数，分量形式z j zje j y j yj ej x j xj ej zz y y x x A N A N A N A A A e e e e ~e ~e ~A ~⎪⎪⎭⎫⎝⎛+⎪⎪⎭⎫⎝⎛+⎪⎪⎭⎫ ⎝⎛=++=∑∑∑===414141 以矩阵表示磁矢量位A 在单元节点上的各分量[][]()z y x l A A A A A Tl l l l e l ,,==4321),,(][][~41z y x l A N A N A el T e j lj e j l ===∑=于是，磁矢量位z z y y x x A A A e ~e ~e ~A ~++=的矩阵表达式∴[][][][][][][][][][][][]e Te z x Te T eT ee z T ee y T ee x Tez y x A A A N N N A N A N A N A A A A 041000000N ~~~~=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=式中：Tz z y y x x x e A A A A A A A A ][][4141421 =[][][][]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=e ee e N N N 000000N 5.3.3 单元分析将近似函数A ~代入方程有余量()J A ~R -⨯∇⨯∇=m γ 取矢量权函数为W i ，有加权余量方程0>=<R ,W i单元剖分后余量加权的和式为零()[]{}∑∑⎰===∑=-⨯∇⨯∇⋅>=<0110Z e Z e i VemieidV εγJ A W R,W由矢量恒等式()()()B A A B B A ⨯∇⋅-⨯∇⋅=⨯⋅∇上式分解为：d()()()()()()()⎰⎰⎰⎰⎰⎰⎰⎰⋅-⋅⨯∇⨯-⨯∇⋅⨯∇=⋅-⨯∇⨯⋅∇-⨯∇⋅⨯∇=⋅-⨯∇⨯∇⋅eeeeeeeeV V i S m i m i V V i V m i m i V i m V i VV V V V VV d d d d d d d d J W S A W A W J W A W A W J W A W γγγγγ分析上式第二项，若有一部分处在外边界上，则由边界条件2s 上：()⇒=⨯⨯∇0n m e A γ()[]()[]0=⨯⨯∇⋅=⋅⨯∇⨯n m i n m i e A W e A W γγ被积函数为零。
若e S 全为内部边界，则相邻单元的交界面上，该面积分的合成将为零，所以该面积分不必计算，可放在i ε中，所以：()()i V V i mieeV V εγ=⋅-⨯∇⋅⨯∇⎰⎰d d J W A W在直角坐标系下，取伽辽金法，权函数为zi yi xi i W W W W ++=z i y i x i z zi y yi x xi N N N W W W e e e e e e ++=++=对于ix Wziy i xizy x xi yW z W W z y x e e e e e W ∂∂-∂∂=∂∂∂∂∂∂=⨯∇00 矩阵表示[][][]xi i i xi W y W z W ⨯∇=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡∂∂∂∂=⨯∇0W 其中[]⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡∂∂∂∂-∂∂-∂∂∂∂∂∂-=⨯∇000xy x z y z ， []()n i N W i xi ,,, 2100=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=按同样的道理，以矩阵形式表示yi W 和zi W[][]T iyiN W 00=， [][]Ti zi N W 00=则i W 可表示为：[][]i i i iN N N N =⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=00000W [][][][][]i N ⨯∇=⨯∇=⨯∇W W单元分析式设()n i z y x l lii ,,,,,W W 21===，则该式的矩阵形式为：[][][][][][][][]()⎰⎰==-⨯∇⨯∇Ve li Tli e T m TyVe li z y x l dV J W dV A W ,,N εγ 而取z zi y yi x xi i W W W e e e W ++= 时，应综合为下式[][][][][][][][][]()n i dV J dV A i Ve T i e T e m TV i e,,,N N N 21==-⨯∇⨯∇⎰⎰εγ式中：[][]Tz y xJ J J =J上式可以分解为3个标量方程式，即：xiVex i nj zj ji yj j i xj j i j i Ve m dV J N dV A x N z N A x N y N A y N y N z N z N εγ+=⎥⎦⎤⎢⎣⎡∂∂∂∂-∂∂∂∂-⎪⎪⎭⎫ ⎝⎛∂∂∂∂+∂∂∂∂⎰∑⎰=1yiVey i nj zj j i yj j i j i xj j i Ve m dV J N dV A y N z N A z N z N x N x N A y N x N εγ+=⎥⎦⎤⎢⎣⎡∂∂∂∂-⎪⎪⎭⎫ ⎝⎛∂∂∂∂+∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂∂∂-⎰∑⎰=1ziVez i nj zj j i j i yj ji xj j i m Ve dV J N dV A x N x N yN y N A z N y N A z N x N εγ+=⎥⎦⎤⎢⎣⎡⎪⎪⎭⎫ ⎝⎛∂∂∂∂+∂∂∂∂+∂∂∂∂-∂∂∂∂⎰∑⎰=1()n i ,,, 21=i 取不同值时应有不同的方程，一阶四面体单元，4,3,2,1=i，每取一个i 值，就有上面一组方程，总的单元分析将出现4组如上方程。
可整理为)4(=n 取：∑∑∑===+=++414141j j xixi zj xz j yj xy xj xxP A C A C A Cε∑∑∑===+=++414141j j yiyi zj yz j yj yy xj yxP A C A C A Cε∑∑∑===+=++414141j j zizi zj zz j yjzy xjzx P A CA CA Cε()4321,,,=i对比以上两组表达式，可得到xx C 、yy C 、zz C 、xy yx C C =、xz zx C C =、yz zy C C =等参数相应的计算式。
进一步整理单元分析的结果，可得矩阵形式的表达式e e e e p A k ][][][][ε+=其中： [][][]()43214321,,,,,,,===l k i k kkli e⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡==zz zx zx yz yy yxxz xy xx i C C C C C C C C C k ][子矩阵 T z y x z y x z y x z y x e p p p p p p p p p p p p p ][][444333222111=Tz y x z y x z y x z y x e A A A A A A A A A A A A A ][][444333222111=由上分析可知：(1) 按有限元法分析磁场的旋度旋度方程，单元中每一代数方程的变量数为3×n 。
(2) 单元分析的刚度矩阵为)3()3(n n ⨯⨯⨯矩阵，系数矩阵中的子矩阵对应为标量位分析中的系数矩阵元素。
4.3.4 综合集成将单元分析中各矩阵[]e k 、[]e A 、[]e p 等扩展为全域范围下的表示形式，综合集成得旋度旋度方程边值问题对应的有限元方程： e p A k ][][]][[ε∑+=∵0][⇒∑e ε得有限元方程：][]][[p A k = 5.3.5 ()02≠⨯⨯∇S n m e A γ时的处理 ()022≠⨯⎥⎦⎤⎢⎣⎡⎪⎪⎭⎫ ⎝⎛∂∂-∂∂+⎪⎭⎫ ⎝⎛∂∂-∂∂+⎪⎪⎭⎫ ⎝⎛∂∂-∂∂=⨯⨯∇S nz x y y z x x y z m S ny A x A x A z A z A y A e e e e e A γγ很显然，一般情况下上述边界条件的处理是十分困难的。
考虑到磁场分布的对称性，适当地设立坐标系或坐标的方位，可以使上述边界条件得以简化。
5.3.6 关于A 解答的唯一性在求解磁场的过程中，有两个问题是需要认真对待的。
其一是由于媒质的非线性问题，它表现为()m γμ的非线性，导致控制方程()J A =⨯∇⨯∇m γ为非线性方程，()m m m γγγ⇒=A 将是坐标的函数，这类问题需要专门讨论，可以参看“工程电磁场数值分析”书。
其二是控制方程中只考虑的A ⨯∇，而未确定A ⋅∇，这使该式的解析解对A 而言不是唯一的。
但由于计算机字长有限，按数值解而言可以得出计算结果，但不能保证0=⋅∇A ，按不同的解法，A 的数值可能不同，即体现了A 的解答不唯一，但A ⨯∇是唯一的，即是说B 的解答是唯一的，这正是我们所需要的。
若按以下方式推导()()[]J A A A =∇-⋅∇∇=⨯∇⨯∇2m m γγ将0=⋅∇A 代入有J A -=∇2m γ就有定解问题：()⎪⎩⎪⎨⎧∈=⨯⨯∇=∈-=∇2201S V n m S m e A A A J A γγ在直角坐标系下，控制方程可以写为()z y x u J A uu m ,,=-=∇2γ即三个独立的标量方程，如是可以按三个标量方程来分别求解，不是就可以免去上面很复杂的方程组求解吗？实际上并非如此。
首先，上述定解问题中三个标量方程彼此并不独立，相对独立的前提是0=⋅∇A ，即0=∂∂+∂∂+∂∂zA y A x A zy x ，可见A 的三个分量彼此相关，不能分开求解。
第二. ()0=⨯⨯∇n m e A γ 边界条件说明x A 、y A 、z A 三个分量在边界2S 处应满足的条件，三个分量不能独立，所以定解问题应为()⎪⎪⎪⎩⎪⎪⎪⎨⎧∈=⨯⨯∇=∈⎭⎬⎫=⋅∇-=∇21200s V n m S m e A A A A J A γγ这才与前面所讨论的定解问题是等价的。
