（等离子体物理专业论文）潘宁阱中电子等离子体的数值模拟.pdf

中田科掌技术大掌司e 士掌位论文擅墨 摘要 本论文的主要工作是成功地编写了一套针对电子束潘宁阱非中性等离子体 装置( e b p t ) 的粒子模拟程序。并且对潘宁阱内电子聚心状态的定标关系进行 了研究。 我们通过全面的考虑和设计，成功地编写了整套程序。进行了一系列的调 试和实验运行后，表明该程序能够正常运行，进行非中性等离子体物理方面的 数值实验工作。目前，我们已经利用该程序完成了对电子聚心约束定标关系的 研究。程序的主要特点有：在旋转坐标系里研究粒子的运动；采用非均匀网格 划分来贴近边界；p i c 法分配电荷和求粒子受力：迭代法求解泊松方程；蛙跳 格式推进粒子。 我们利用该程序对e b p t 的电子等离子体进行了数值实验研究，主要的实 验结论有： 1在e b p t 中心附近可以形成密度可以突破布里渊极限的稳 态电子等离子体。模拟结果得到稳态时潘宁阱中心的电子 密度可以达到l o ”。左右。 2 稳态时，电子的势能最小值位置不是潘宁阱中心而是内 部平均密度为布罩渊密度的径向位置( r = 1 5 m m ) ，电子 能谱在低能端出现截断。 3 潘宁阱中心会出现一定程度的电子势垒，把低能( 截止能 量以上) 电子约束在阱中心附近。 4 数值结果表明了电磁场匹配条件的满足是达到较高聚心的 必要条件，与e b p t 的实验结果一致。 5 阱中的电子总数对稳态中心密度的影响不大，电子较多 时，相应的近热化的约束电子增多；电子运动范围增大， 即能量展宽增大，聚心密度越大；聚心密度与有效势的大 小成正比，与p f x 和e b p t 实验结果一致。 6 达到稳定的时阃与初始电子数成反比，同碰撞驰豫时间与 电子密度成反比的关系相吻合；达到稳定的时间与潘宁阱 的有效势有近似线性的正比关系。 以上对潘宁阱中的等离子体定标关系的研究，为进步研究潘宁阱中的非 中性等离子体性质提供了必要的、可靠的数据。) ， ! ! 苎竺苎查查竺! 主竺苎丝奎 一一生翌墨叁! ! a b s t r a c t b e c a u s eo f t h e p o o ra c c e s s i b i l i t ya n d t h ei n c a p a b i l i t yo f e m i s s i o n ，i ti sd i f f i c u l tt o d e t e c td i r e c t l yt h ep a r a m e t e r so fc o n f i n e de l e c t r o n si nt h ee l e c t r o nb e a mp e n n i n g t r a p ( e b p t ) b e s i d e s t h e e x p e r i m e n t s o nt h ee b p t , n u m e r i c a ls i m u l a t i o ni s a p o w e r f u lw a y t os t u d ys u c hp l a s m a si nd e t a i l i nt h i st h e s i s ，ap a r t i c l e 。i n - c e l l ( p i c ) c o d eh a sd e v e l o p e dt oa n a l y z et h e s t a t eo fe l e c t r o np l a s m a si nt h ee b p t t h e c h a r a c t e r i s t i cf e a t u r e so ft h ec o d ea r e ：r o t a t i o n a lc o o r d i n a t e sa r eu s e d t oe l i m i n a t et h e s t a t i cm a g n e t i cf i e l d ；n o n - u n i f o r ms p a t i a lg r i d sa r es e t t om a k et h ec o m p u t a t i o n a l b o u n d a r yb e a sc l o s et ot h es u r f a c eo ft h ee l e c t r o d e sa sp o s s i b l e ；p i cm e t h o di su s e d t o a s s i g np a r t i c l ec h a r g et on e i g h b o r i n gg r i dp o i n t s a n dc a l c u l a t et h ef o r c ea tt h e p a r t i c l e s ；t h ep o i s s o ne q u a t i o ni ss o l v e d w i t hi t e r a t i v ea l g o r i t h ma n dt h ep a r t i c l e sa r e a d v a n c e db yl e a p - f r o gm e t h o d t h es t e a d ys t a t eo fe l e c t r o np l a s m a si ne b p ta n dt h ea s s o c i a t et i m ee v o l u t i o n p r o c e s s a r es t u d i e db yt h i sp i cc o d e t h em a i nc o n c l u s i o n sa r ea sf o l l o w s ： 1 t h es t e a d ys t a t eo fe l e c t r o nf o c u s i n gi sf o u n di nt h ee b p tw i t h t h ec e n t r a ld e n s i t yo f1 0 y b o rm o r e 2 t h es e l f - f i e l d o fc o n f i n e de l e c t r o n sd e c r e a s e st h ee l e c t r i c p o t e n t i a li nt h et r a p t h em a x i m u mp o t e n t i a li sn ol o n g e ri nt h e c e n t e r o f p e n n i n g t r a pb u t i n t h e p o s i t i o n o f r 2 1 5 m m 3 t h e r ee x i s t sab a r r i e rp o t e n t i a li nt h ec e n t e rf o re l e c t r o n s ，i n w h i c hs o m e l o w - e n e r g y e l e c t r o n sa r e t r a p p e d t h et r a p p e d e l e c t r o n sa r ea l m o s tt h e r m a l i z e d 4 i ti sa l s oi d e n t i f i e dt h a tt h ee l e c t r i cf i e l dm a t c h e st h em a g n e t i c f i e l di st h en e c e s s a r yc o n d i t i o nf o rs t r o n gf o c u s i n go fe l e c t r o n s i n e b p t 5 t h ed e n s i t yo ft h ef o c u si sa l m o s ti n d e p e n d e n to fa m o u n to f t h e i n i t i a l l yi n j e c t e de l e c t r o n s b u ti t i sp r o p o r t i o n a lt ot h ew i d t ho f t h ee n e r g yd i s t r i b u t i o na n dt h ee f f e c t i v ep o t e n t i a l 6 t h et i m et or e a c ht h es t e a d ys t a t ei si n v e r s e l yp r o p o r t i o n a lt ot h e i n i t i a l l yi n j e c t e de l e c t r o n s ，a n dp r o p o r t i o n a l t ot h ee f f e c t i v e p o t e n t i a l i i 中田科学技术大学硕士学位论文致谢 致谢 在本论文完成之际，感谢所有指导和帮助过我的人们。 首先对我的导师刘万东教授和杨维宏副教授致以深深的谢 意。论文工作的每一步进展都离不开他们对我的悉心指导。他 们渊博的学识、严谨的治学态度使我受益非浅。 感谢郑坚博士和远在法国的庄革博士几年来给我的热心鼓 励和直接指导。 感谢俞昌旋教授、胡希伟教授和马锦秀教授在百忙中抽出 时间给我的指导和支持。 感谢本实验室的谢锦林博士、赵凯、于治、梁小平、迟霁 对我工作和生活上的协助。 付成江 中啊科掌童诛大掌习e 士掌位论文第一幸5 i t 和i 遗 i 1引言 第一章引言和综述 1 9 3 6 年潘宁( f m p e n n i n g ) 提出利用外加的静电场和静磁场构成的潘 宁阱装置( p e n n i n gt r a p ) 来约束带电粒子，此后潘宁阱主要应用在原子分子物 理、粒子物理实验领域中，约束少量带电粒子。近些年来，潘宁阱对非中性等 离子体所具有的良好约束位形受到人们的日益重视，用于研究纯电子( 或离子) 、 正电子或反质子等离子体的平衡、稳定和输运性质。l t u r n e r 和d c b a r n e s 等人提出在潘宁阱中可以得到高度聚心的电子等离子体状态，以用于受控热核 聚变研究【2 】。我们实验组建立了一套电子束潘宁阱系统( e b p t ) 来进行这方面 的相关实验i “1 。由于潘宁阱的形态比较封闭，电子等离予体不具有自发或受激 辐射的特性，对其参数的诊断，特别是等离子体参数的空间分布比较困难。本 文建立了一套完整的粒子模拟程序来描述阱中的电子等离子体状态，研究其中 电子等离子体的运动规律及参数定标关系，以达到与实验相互印证和互补的目 的。作为引言，我们首先介绍一下惯性静电约束的原理，然后介绍利用潘宁阱 进行聚心实验的原理和相关的实验结果。 1 2 惯性静电约束原理及进展 受控热核聚变是高温等离子体物理学的主要研究领域，半个世纪以来一直 受到各国科学家的关注。惯性静电约束( i n e r t i a le l e c t r o s t a t i cc o n f i n e m e n t 或 i e c ) n 是最早尝试应用于受控热核聚变研究的方式之一，并且获得过一定通量 的热核聚变中子。它的原理是在一外部建立的球对称静电势阱中，注入( 或由 背景中性气体电离) 离子使其加速会聚至球心。当外加势场足够大、中心处的 离子密度足够高时，可使中心达到热核聚变所需的温度和密度，产生热核聚变 反应。但由于非中性等离子体m i 中自电场作用的影响，中心的密度受到布里渊 ( b r i l l o u i n ) 密度极限的制约，从而限制了整个体系的聚变反应率。一般而言， 布里渊密度远低于聚变反应所需的密度。密度的限制以及i e c 装置内电极对注 入离子的吸收几和过大的热负载使得静电惯性约束热核聚变的研究在七十年代 以后一段较长时间内停滞不前。 1 9 9 3 年，美国的洛斯阿拉莫斯国家实验室( l o s a l a m o s n a t i o n a l l a b o r a t o r y 或l a n l ) 的l t u r n e r 和d b a r n e s 等人提出在潘宁阱中，利用低能、小正则 中圈丰h 扛技术大学硕士攀位论文第一窜引富和壤篷 角动量的电子束注入方式，可以得到高度聚心的电子等离子体状态”, 9 1 。他们注 意到在非中性等离子体的情况下，布里渊密度极限仅在体积平均上才有意义。 通过适当的方式可以在约束体积范围之内产生个比较大的空间密度梯度，从 而在约束的局部区域内密度可以突破布里渊极限。此时的非中性电子等离子体 的自电场可以形成类似于i e c 中由外电场构造的球对称离子势阱。聚心电子团 构成了虚阴极，取代了i e c 装置中的内电极，消除了使用实物电极的障碍。在 保证虚阴极存在的条件下，注入的离子可以加速进入到密度极高的非中性( 电 子为主) 等离子体中心，获得聚变反应所需的能量产生聚变。 1 3 潘宁阱中球对称有效势的产生 典型的潘宁阱的结构如图1 1 ，主要由内表面为旋转对称的两个端极和一 个环极组成。三个电极封闭在真空室内，根据约束粒子的不同，在环极端极之 间外加正电压或负电压，可以在阱中产生电四极势，对阱中带电粒子提供轴向 约束：由平行于对称轴方向的均匀外磁场对带电粒子提供径向约束。 图1 t 潘宁阱原理图 一般而言，端极、环极是会切的旋转双曲面。在图1 1 所示的以阱中心为 原点柱坐标歹“ 上、下端极的剖面双曲线为2 2 2 一，2 = 2 z ：，而环极的为 r 2 - - 2 2 2 = 疗。这里z 为轴向坐标，r 为径向坐标，沿柱对称轴的上、下端极间 距为2 ，沿径向环极的最小半径为，0 。此时，在环极上相对两端极施加一电压 圪，则在柱坐标系( ，z ) 中由电极应用电压产生的静电四极场( 设坐标原点的 中豳科掌崮l 术大掣嘎士学位惋。文第一幸引t 和缘遗 电位为零) - - f 以表示为3 盹加去( ，2 _ 2 2 2 ) ( 1 1 ) 在以角频率为拉莫频率国。的旋转参考系中，由外加均匀磁场b o 和加在端 极和环极之间的电压k 的形成的潘宁阱电磁场为一轴对称的静电场，可以表为 柱坐标系p ，2 ) 中的空间等效电势蝣( 设坐标原点的电位为零) = 焘( ，2 _ 2 z2 ) 一器一 z ， 等效空间电势其实就是潘宁阱的电磁场在一特定角频率的旋转参考系中的重新 表述。( 1 2 ) 式右边的第一项为电场的贡献，第二项为磁场的贡献。 通过调整外加的磁场和四极静电场，潘宁阱的电磁场可以提供我们需要的 有效球形势阱。当电场和轴向磁场满足 = 。 ；e b 0 ( t 0 2 + 2 z 2 0 ) ( 1 3 ) 2 4 m 。 形成球对称等效电场球。我们称( 1 3 ) 式为潘宁阱电磁场匹配条件。在电磁场 匹配情况下，其等效电势为 驴一盟1 2 m , o2 ) = 一矗p 2 。， 这里，p ：；“：+ z 2 o 由于等效电势为球对称形，在潘宁阱中运动电子的受力 将指向球心。正则角动量很小的电子在它逃逸前将反复穿越势阱中心附近。 1 4 利用潘宁阱进行聚变实验的具体实现 如果潘宁阱中约束的非中性等离子体中存在较大的密度梯度，那么它的密 度在某一局部如芯部是能够超越布里渊密度极限的。实验上采用零角动量的低 能电子束注入潘宁阱的方式来实现超越布里渊密度的稠密电子芯。 在潘宁阱的约束位形中，粒子的运动保持能量守恒与正则角动量守恒。正 则角动量守恒是指 易2 ，( p 日一m r o ) 。2 ) m r 2 p c o j = c o n s t( 1 5 ) 也就是说在电子的运动过程中，b 将保持不变，其中( ，口) 为柱坐标。 中国科学技术大掌司e 士掌位论文】疗一 粤f 曹- 和崩i 述 在初始注入电子的正则角动量很小( bzo ) 的情况下，b 守恒就意味着电 子在运动的过程中，r 很小( r * 0 ) 抑或毋z 国，2 。所以电子在每一次运动周期 中，将穿过原点附近一次。如果在潘宁阱内整个空间中分布着很多这样运动的 电子，则在潘宁阱内可以实现较大的密度梯度，阱中的电子密度将与力戎反比， 从而在原点附近将会出现较高的密度分布。对于这一过程，在潘宁阱电磁场匹 配的条件下，也可以在旋转参考系中考虑。 实验中为了保证注入的电子具有小正则角动量，将电子沿下端极中一t 2 , 轴的 , f - 于= l 注入到潘宁阱之中( 即，z0 ) 。实际的操作时，在上端极上加一接近于地 的负偏压，电子束沿中心轴从下端极注入到阱中，然后在上端极附近沿轴向折 回。入射束和折回束在阱的中心附近会碰撞而发生偏转，把注入电子的部分轴 向动能转化为径向的动能，从而使得注入电子有了径向运动的自由度。当少量 的偏转电子被潘宁阱的势场约束，在阱中旋转时，将使中心区域的等效空阳j 电 荷密度增加，从而使得继续注入的电子束更容易通过碰撞发生偏转，被潘宁阱 的势阱所约束，进一步增加中心空间电荷密度。上述过程是一个典型的f 反馈 过程，如果能够持续进行在中心附近所造成的等效电荷密度就可能会超越布 罩渊密度极限。 为了保证偏转了的约束电子能够重新折回并经过原点附近，有必要提供一 个合适的约束势场，这就是潘宁阱中在旋转坐标系下的球形势阱。在这个等效 球形势阱的作用之下，阱中零角动量的电子轨迹在获得足够的能量逃逸前会反 复经过阱的中心。大量的电子运动在这类径迹上，所造成芯部的密度就可能会 超越布里渊密度极限，实现和建立起电子聚心状态。 1 5潘宁阱上的聚心电子实验结果 到目前为止，已经有两个小组利用这种原理进行了相关实验。我们下面简 单介绍一下两个相关的实验结果。 1 5 1p f x 上的实验结果 1 9 9 7 年美国l a n l 的d b a r n e s 领导的小组在p h y s i c a lr e v i e wl e t t e r s 上报 告实验结果，宣布在潘宁聚变实验装置( p e n n i n gf u s i o ne x p e r i m e n t 或p f x ) 内 实现电子聚心状态，芯部的电子密度高达布里渊密度的3 5 倍 1 0 - 1 2 l 。p f x 主要的 实验结果如下： 中曲科掌技术犬掌硬士掌位论文第一章引言和壤避 l 、电子聚心状态的形成：通过观察环电极收集电流与中断约束m c p 收集电子 与电磁场状态的关系，可以判断电子聚心状态已经形成。 2 、注入电子流阈值及其与电磁场的关系：在潘宁阱形成等效的球形势阱情况 下，实验发现出现电子聚心状态对注入的电子束流有阈值要求，且此阂值与 施加的电场成正比。 3 、聚心电子密度：p f x 用基于f o k k e r p l a n c k 方程的模型，对所得的实验结果 进行数值模拟计算，结果表明在电子聚心时芯部密度可以达到布里渊密度的 3 5 倍。 1 5 2e b p t 上的实验结果 我们实验室设计和建造了一台电子束潘宁阱( e l e c t r o nb e a mp e n n i n gt r a p ， e b p t ) 非中性等离子体装置，于1 9 9 8 年正式运行。主要用于进行非中性等离 子体物理方面的实验，目前己完成电子聚心约束实验。e b p t 主要的实验结论 为： 1 、通过对有能量增量的逃逸电子的测量，证实了在e b p t 装置中，获得了有较 强聚心的非中性电子等离子体。 2 、与p f x 实验相类似，证实了潘宁阱电磁场匹配是实现较强的电子聚心状态 的必要条件：电子聚心状态的形成和建立对注入电子的束流有阙值要求；束 流阈值与环极上的应用电压成线性关系，表明在此装置中，约束电子的约 束时间与磁场强度平方值成正比，符合经典的碰撞输运系数的定标关系。 3 、利用电子束方波开关调制注入方式，着重研究了电子聚心状态的时间演化特 征，发现在此装置中电子聚心状态是不稳定的。我们建立了一个自洽的物理 模型，定性和半定量地解释了这种时间演化行为。上升过程为捕获电子的累 积和电子聚心状态建立的过程，下降过程则为电子的热化过程，热化导致了 电子的能量和正则角动量展宽，芯部密度下降和有效的约束体积减小，从而 破坏了电子聚心状态。 4 、在一定的磁场范围( 磁场较小) 内，潘宁阱的聚心电子密度与磁场强度的平 方值，或布里渊密度成正比，即聚心电子密度与布里渊密度之比与不随磁场 变化，表明理想的潘宁阱约束电子的定标关系成立。但在较大的磁场时，这 一比值随磁场增大而减小，理想的定标关系遭破坏。 5 、根据实验数据和建立的自洽模型，可以估算电子聚心状态时芯部的电子密度 约为布里渊密度极限的l o 倍。 中田科掌执术大掌习e 士掌位 文第一幸引t 和壤遗 1 6 本文的研究目的、内容以及意义 e b p t 潘宁阱系统的三个电极安装在个不锈钢的圆筒真空室内，其可近 性很差，无法对阱中的等离子体参数进行直接的诊断。现有的各种诊断方法都 只是间接诊断，例如通过法拉第筒或微道板测量逃逸电子的能量、速度和空间 分布来间接测量阱内的等离子体状态。这些测量的准确性很难得到保证。本论 文的目的是通过数值模拟的方法，分析e b p t 中电子聚心状态，了解其中的物 理过程，了解阱中的等离子体参数，以期与实验相互补充和相互印证。论文内 容主要是针对e b p t 的位形中，对阱中的电子等离子体进行粒子模拟。首先在 调研的基础上编写了整套粒予模拟源程序。程序采用非均匀的网格划分，迭代 方法求解泊松方程，获得粒子与场的瞬时空间分布；用蛙跳法推进粒子运动， 实现系统的时间动力学演化。通过数值实验，分析研究了形成电子聚心的条件， 研究了电子在阱中的分布及演化，得到了潘宁阱中的等离子体参数与各种运行 参数之间的定标关系。从而为进一步研究潘宁阱中的非中性等离子体性质提供 必要的、可靠的数据。 中田科掌技术大掌疆士掌位论文第= 章潘宁阱中非中性等离子体的截值模拟 第二章潘宁阱中非中性等离子体的数值模拟 我们在第一章中对潘宁阱中实现电子聚心状态的实验原理及相关的实验结 果作了简单地介绍。这一章中，我们将介绍潘宁阱中非中性等离子体的数值模 拟，首先是粒子模拟的基本原理和一些有关的概念，然后是针对潘宁阱位形的 粒子模拟程序的具体实现方式，最后介绍程序的流程。 2 1 粒子模拟的概念 随着计算速度的提高、内存的增大以及广泛的普及，计算机模拟在各种领 域中越来越重要。等离子体物理研究领域，计算机模拟已经成为与实验、理论 相比拟的研究方法，针对不同应用目的和背景人们已经建立了许多程序，对 复杂的等离子体体系的行为进行研究，并能得到具有一定可靠性的结论。计算 机模拟可分成数值模拟与粒子模拟描述两种，数值模拟通过对已建立的理论模 型的数值求解，得到解析计算中无法实现结果。粒子模拟则以等效的微观粒子 作为模拟的主体，模拟粒子在各种力场的作用下按照动力学规律进行运动的情 况。由于计算机内存及速度的关系，直接对微观粒子的模拟是不可能的，原因 是体系的微观粒子数太大，同时直接计算粒子之间的两两相互作用计算量太大。 所有数值模拟通常是用一个模拟粒子代表了一团相同的粒子，同时对粒子间的 相互作用的处理采取粒子产生场、场对粒子进行作用的间接作用的方法。 为了消除点电荷对自身作用的奇点通常对带电粒子用粒子云来处理，即 模拟粒子代表了一团微观粒子云，同时对外界及粒子云本身产生空间场在空间 进行适当的网格化，这即是在等离子体粒子模拟计算中采用粒子云一网格法 1 “。 这里我们以静电模式为例，粒子云一网格法的要点概述如下： l 、粒子由其位置和速度( x 。，v 。) 设定，坐标x 。在模拟域中连续变化，表示第p 个粒子的中心位置。 2 、每个粒子代表一个电荷为有限分布的粒子云。粒子云的电荷g 。按定的规 则分配到x 。邻近的网格点一，上，统计所有粒子对网格点的贡献，计算出 网格点上的电荷密度n ，。t ，表示第( j j 力个网格点的位置。 3 、由离散的电场方程( 在静电模式下) v e = p 6( 2 1 ) 计算网格点上电场e 。 中囡科掌技术大攀硬士掌位论支第；章潘宁阱中非中性等离子体的戴值模拟 4 、按一定规则计算与第p 个粒子有关( 与屯邻近) 的网格点电场对它的作用 力( ) 。 5 、由力t 把粒子p 按牛顿方程( 非相对论情况) m p 每吨 出。 d t p ( 2 2 ) 推进一步。同样的运算扫遍所有粒子，完成整个系统的一步时间演化。 以上做法的关键在于：第一，不直接计算粒子与粒子之间的静电作用，而 是用宏观的电场来传递这种作用：第二，粒子对电荷密度的贡献并不在点粒子 的所在位置统计，而是在个或若干个相邻的网格点上统计。这样就把点粒子 对电场的贡献变为有限体积的粒子( 粒子云) 对电场的贡献，从而消除了点粒子 电场的近距离发散。同样，用若干网格点电场对粒子的作用力，代替粒子所在 位置电场对粒子的作用，也把点粒子相互作用变为有限体积的粒予云的相互作 用。 粒子云网格法的物理图象为，很多个体积不为零的粒子云，在计算域中游 动，并可以互相重叠、穿透。它们之间的相互作用则通过它们各自贡献于网格 点上的电磁场，再由此电场、磁场作用于每一个粒子云来实现。粒子云与粒子 云之间不直接相互作用。 与通常的统计物理及磁流体力学方法相比较，粒子云网格法是一种本质上 全然不同的描述等离子体的方法。物理上，由于它舍弃了统计的概念与方法， 力图回归到对个别粒子的描述，就更能表现出个别粒子在等离子体集体行为中 的作用；数学上，由于对粒子的处理舍弃了流体处理中的解析方法和偏微分方 程组，也就绕过了由微分到差分所带来的收敛性、稳定性等可能的数学困难。 可以说，等离子体粒子模拟是至今最彻底的模拟。 2 2 坐标系的选取 我们的模拟工作是在p c 机上进行的。由于p c 机性能的限制，只能进行二 维的模拟，我们选择二维速度和二维空间。用于约束电子等离子体的潘宁阱的 构造原理如图1 1 。如1 4 节所述，在旋转坐标系中，有效势是球对称的。根 据这种对称性，我们选择柱坐标系进行模拟。如图1 1 所示，选择，、z 二维空 间和速度为我们的模拟空间。以磁场的方向为z 轴，垂直于磁场的方向为r 轴， 潘宁阱的中心点为坐标原点。在这种情况下，程序中的模拟粒子实际上相当于 三维空间里以z 轴为对称轴的电子环，粒子在极向的运动不予考虑，电子环在 中田科掌拉术大掌硬士掌位论文第= 章潘宁阱中非中性等育子体的数值模拟 球对称的等效静电场中运动。在球对称情况下，用球坐标的径向坐标作为变量 更好，为了区分，本文中用r 代表柱坐标系下的径向坐标，胄代表球坐标系下 的径向坐标。 2 3 边界条件的处理和网格的划分 2 3 ，1 等效电场的边界条件与网格的划分 我们实验所用的潘宁阱的上下端极是会切的旋转双曲面，这使得泊松方程 的求解非常困难。一般的快速解法都要求解域比较规则。而我们的问题不仅是 双曲线边界，而且解域不是封闭的，在端极和环极的交会点是奇点。经过推导， 发现很难找到合适的坐标变换或保角变换。在没有更好的方法之前，我们采用 迭代法来求解泊松方程。虽然迭代法的运算效率比较低，但它能够灵活地处理 边界条件并控制求解的精度。 因为e b p t 的电极电压由外加电压控制，把泊松方程的边界取在电极表面， 可以得到较为准确的数值。网格的划分则围绕着如何贴近边界这个问题来进行。 为得到足够的精度，我们采用均匀网格结合非均匀网格的划分方法。在阱中心 为均匀网格，电极附近则为非均匀网格，使大多数的边界点取在电极表面上。 由于潘宁阱的端极和环极剖面是两条无限接近的双曲线，我们不能也没有必要 取太多的边界点在它们交会点附近的曲线上。在这些地方，我们把一些边界点 取在阱中的空间内，用阱中没有电荷时的电四极势解析解来作为边界值。由于 在阱中运动的电子很少能到达交会点附近，电子产生的自恰势对其电势的影响 很小，所以由此而引起的误差是我们的计算能够接受的。下面我们来详细介绍 一下网格的划分方式。 潘宁阱的上、下端极的剖面双曲线方程为2 2 2 一r 2 = 2 z ：，而环极的为 r2 - - 2 z2 = 求，这里z 为轴向坐标，r 为径向坐标，沿柱对称轴的上、下端极间 距为2 2 0 ，沿径向环极的最小半径为。在我们的实验装置中，2 为2 1 2 r a m ，h 为1 5 m m 。则四条直线，= 0 、，= r 0 、z = 一面和z = z o 围成一个矩形。在这 个矩形的区域内，将r 方向均匀划分3 0 等分，z 方向均匀划分4 0 等分构成矩形 网格，r = 0 5 哪，a z = 0 5 3 眦。这些网格的边线的延长线与双曲线的边 界相交，得到与端、环极相交的交点。从与，轴平行的已有直线与边界的交点 作z 轴的平行线；类似的，从与z 轴平行的已有直线与边界的交点作，轴的平 行线。这些平行线在均匀分割的矩形外面相交，与双曲线边界组成非均匀网格。 注意在渐近线附近取最外侧的直线为边界。最终得到的网格中心均匀部分为 3 1 x 4 1 个网格，总体为1 0 1 x 5 1 个的非均匀网格。这样的划分方法可以保证边界 最大限度的贴近电极的表面，从而由电极电压得到准确的边界值；中心的均匀 中田科掌技术大掌习e 士学位论文第：；章潘宁阱中葺e 中性等，r q 蕾的数值期拟 部分又简化了粒子运动的计算。 2 3 2 粒子运动的边界条件 粒子运动的边界包括电极表面和柱坐标系的中心对称轴。在程序运行一定 时间后，每步都可能有部分粒子到达电极表面。我们假设电极是完全吸收的介 质，并且电极的电压不变，二次电予发射也被忽略。当粒子出了到达边界后， 表明它已经获得了足够多的能量从潘宁阱中逃逸，记录它的位置与速度后把它 从计算循环中去掉。由于我们选择的旋转坐标系的轴对称性，= 0 的对称轴形 成了计算边界，此边界相当于一弹性反射界面，当粒子穿越中心轴，= 0 时，要 将它的r 轴速度和坐标进行反转。 2 , 4 粒子初始条件的设置 一般给出的初始条件是粒子的统计分布函数f ，v ) ，我们需要根据分布按 统计意义给出每一个粒子的初始位置x 和初始速度p 。通常把位置空间设定与 速度空间设定分开考虑。设粒子总数为氘我们以粒子在半径r ，的球内均匀、r ， 外按l r - 分布的情况为例，说明如何根据统计分布设置粒子的初始条件。这种 分布是潘宁阱聚心常采用的一种空间分布1 6 ”，即 i 啊足) = ”。r r 沁= 簪耻州。 他3 r 为球坐标系下的径向坐标，r 。为分布的最大半径，n 。为常数。在区间( 0 ，r ) 之 间的粒子数为 n ( r ) = 4 x 【”( p ) p 二d p ( 2 4 ) 定义积分分布函数月尺) ： 删= 意当 眨s ， 注意只0 ) 2 o j f ( r 一) = l 。因为( r ) 0 ，所以，( r ) 为单调增函数，把目r 。) 分 成k - 1 等份，求解以下方程： ，( r ) 2 而i - i ，( 7 = 1 ，2 ，k ) ( 2 6 ) 就可得到满足粒子初始分布的每个粒子所处的球半径r 。 因为粒子的初始分布是球对称的，故可以认为在二维坐标下，粒子与坐标 原点的连线和r 轴的夹角口在0 - n 之间是均匀分配。用上面类似的做法即可求 出只功，进而求出每个粒子的方位角6 1 。 要注意粒子的球半径r ，和方位角6 ：都是递增的，产生了不符合初始分布的 1 0 中田科掌技术大掌硕士掌位论文 第= 章潘宁阱中非中性等离子体的数值概报 相关性，因此要采用随机方法把它们的次序打乱。有了粒子的球半径r ，和方位 角岛，就可以唯一确定粒子的位置。使用类似的方法就可以由粒子速度分布函 数确定每个粒子的初始速度。 2 5电荷分配与粒子受力 2 5 1 电荷分配 电荷分配有许多种方法，考虑到精度和计算效率，我们采用网格中粒予云 法( p a r t i c l e s i n c e l l ，简称p i c ) 。把二维空间( ，z ) 划分为矩形网格后，格距为 a r ，z 。每一格点都占据 一块以该格点为中心的面 积为r 止的空间，即网 格点( ，) 所在的室。 把粒子p 看成与网格形状 相同的粒子云，边长甜和 业，电荷均匀分布。一般 来说，每个粒子云可以跨 越4 个邻近室( 均匀网格 情况下) ，p i c 法在分配电 荷时，按电子云落入各室 的电荷量，累加到各室对 应的格点上。 这种分配方法，相当 于把粒子云落入各室的面 i 一1ii + l - i - 一一一1 一 - i i ii i ? l 一卜- ( r 。 z 。) 1 - + i ；o i i ； i i i i iii iil 图2 1 粒子云电荷分配示意图 积作为分配于该网格的权重。用图2 1 示意的方法，更容易计算这些面积。图 中实线代表网格线，网格线的交点为网格点，点划线构成的矩形是以中心为网 格点的网室，粒子p 所在的点为( r p ，z p ) ，相应的粒子云由虚线所示的矩形所代 表。显然，粒子云落入( f ，j ) 、( f + 1 ，) 、( f j + 1 ) 、( f + l ，尹1 ) 四个网格点所 在的网室中。于是电荷分布的权重函数可写出如下： 彬，j ( 0 z ，) = ( 0 一0 1 ) ( z ，一2 川) ，业， 现麓：巍刊麓若( x - - x p - - x , + i y y y ) 眩，彬p l ( 0 z p ) = ( o 一+ 1 ) ( z p z j ) 以k ， 一 l 彤+ 1 + 1 ( r p z p ) = ( r p 一) ( z p 一：，) a r a z ， 对其它格点( f ，- ，) ，彬，( 0 ，：，) = 0 a 中固科攀扭拳夫攀硬士学位论文第= 幸潘宁阱中非中性等离子俸的致仉模拟 2 5 2 粒子的受力 得到网格点上的电场后，可以用插值法求出粒子所在位置的等效电场强度 值，从而简单地求得粒子的受力。这等价于把粒子当作粒子云，用( 2 7 ) 式的 权重函数分配到邻近四个格点上，分别求出每格点上所分配电荷受到的电场 作用力，这些力的合力就是粒子受到的电场力，它与粒子的电荷之比即是粒子 位置的等效电场e 。我们在程序中所用到的权重函数与上述分配电荷的权重函 数相同，这样可以可以避免自身力的作用( 即粒子加速它自己) 1 5 l o 2 6 泊松方程的求解与有效场强的计算 虽然迭代法求解泊松方程的速度比较慢，但它的边界处理灵活，精度也可 以有足够的控制。在柱坐标系下，泊松方程v 2 u = f 的表达式如下： 乓“+ ! 塑+ 三“：， ( 2 8 ) 研2r 甜a z 我们采用五点法进行差分，如图2 2 所示，在( i ，) 点的差分值涉及周围 ( i + 1 ，) 、( i ，p 1 ) 、( i - l ，) 、( i ，j - 1 ) 四个网格点，为方便计，我们令这五点的 i 1ii j r l i _ 】广_ l h 2 i + t 图2 2 两维差分格式图 j + 1 1 - l 标号分别为0 、1 、2 、3 、4 、5 。设，向网格间距为h ，：向间距为也，当，0 时泊松方程的差分式为忡i ； 陪+ 爿等+ 陪一爿等+ 斧+ 斧飞降+ 考 - ，o = 。c 2 m 中豳科攀技术大学硬士掌位论文第= 章潘宁僻中非中性等离子体的藏值攘拟 当，= 。时，即节点位于对称轴上时，需另外处理。这时( 就- 0 _ 。，由洛必 达法则可知 啡考) = 旺 所以对称轴上的泊松方程为 2 塑o r2 + 参= ( 2 出2 。 其差分式为 寿“，+ 吉“! + 古蝴一l 毒+ 寺卜一f o o - q 挖， 可q + 虿如+ 可蝴一l 可+ 虿严一_ o _ 1 2 非均匀网格部分与此类似，只是间距不同。按2 3 1 划分网格后，就知道 了各网格间的间距，由( 2 9 ) 或( 2 1 2 ) 就可以得到任何点的差分系数( 即 方程中“o 、“l 、u 2 、u 3 、前的系数，记为a 、b 、c 、d 、e ) 。 给定粒子初始位置后，就可以得到空间电荷分布( 2 7 介绍电荷密度的 计算方法) 。由已知的边界条件，就可以迭代求解泊松方程了。第一次求解时以 拉普拉斯方程的解析解作为迭代的初值，反复求解时以老时间步的结果为新时 间步的初值。根据( 2 9 ) 或( 2 1 2 ) 式，用周围四点的前次迭代值求出中一心 点的新值，所有的新值求出后，即完成一次迭代。 在迭代中，通过控制残差决定所需要的运算步数，进而控制精度。每步的 残差r e s i d 定义为： r e s i d = a u o + b “l + c “2 + d “3 + e “4 一f o ( 2 1 3 ) 于是新的可表示为： “：“= u o “一国掣 ( 2 1 4 ) 式中为松弛系数。迭代前，首先计算所有网格点上f 值的和作为残差的初始 值。每完成一步迭代后计算所有网格点r e s i d 值之和，与初始残差比较，当减 少到e p s 倍时即认为精度达到要求。终止迭代。我们通过权衡精度和运算时间， 在程序中设置e p s 为1 0 。 得到新的电势后，根据事例中磁场的强度，根据( 1 2 ) 式计算得出有效势 的数值。在我们的计算中，粒子很少到达非均匀分割的网格区域。为了简化计 算，只求解均匀区域的有效场强。使用中心差分方法可以得到二极精度。在 1 0 1 5 1 的网格上，对潘宁阱的边界情况v 2 庐= 0 的拉普拉斯方程的迭代解与解 析解的绝对误差分布如图2 3 所示，此时电极间的电压为l o k v ，中心部分的误 中囡科掌拭术大掌硬士掌位论文 第= 章潘宁阱中非中性等离子体的双值模拟 差很小。奇点( 渐近线) 附近的误差大一些，但也只有0 2 v ，只有整个极间电 口1 0 菩o1 毫 - 0 2 - 03 2 0 1 5 图2 3 拉普拉斯方程的绝对误差 压的o 0 2 。经过单粒子简谐运动检验，在经过1 0 0 0 个周期的运行后，粒子运 动轨迹与解析的轨迹的偏差在半个网格尺度内，精度满足我们模拟的要求。 2 7 粒子推进的格式 当粒子的所在空间位置的等效电场e 确定后，粒子的运动简单地由牛顿方 程所确定。对第p 个粒子： 一d v p ：生生 d t m p 出。 2 d t 2 ”n ( 2 9 ) 常微分方程组有很多种很好的数值解法。对我们的粒子模拟来说，计算量 成为突出的问题，最常用的简单而精度又满足要求的积分格式是蛙跳格式。把 时间离散化为 f n = 聍r1 ( 2 1 0 ) ，= + l 2 ) a t j 按数值模拟文献的习惯写法，物理量右上角的上标均表示该量对应的时刻。例 如式中，”、r ”“2 分别表示第h 个离散时刻和第”个与第n + l 离散时刻的中间时 刻，或称第n 个半时刻：a t 为时间步长。 在蛙跳格式中，时间离散化的特点表现在引入半时刻。所有的粒子位置及 中圈科掌擅术大掌硕士掌位论文第；章潘宁阱中非中性等离子俸的敦值模拟 依赖于粒子位置的物理量都在整时刻取值。而所有的粒子速度及依赖于速度的 物理量都在半时刻取值。对我们的近似静电情形，推进格式可写成： v ”n ：v ”l ，2 + q , x t e ”i m，( 2 1 1 ) x n + 12x ”+ a t v “2i 蛙跳格式兼有精度高、计算简单、不需要中间存储的优点。蛙跳格式中，要求 速度始终比位置坐标落后半个时刻，因此设定了t = 0 时刻的粒子位置和速度后， 还要把速度的初始值后退半个时刻。蛙跳格式的相对误差与时间步长的二次方 成正比，这也是由蛙跳的中心差分带来的优点。但若，取得太小，则所需要的 计算步太多。在本论文计算中，e b p t 中小正则角动量电子的运动周期为1 0 。9 秒 量级，故取a t = 5 0 1 0 。秒。 2 8 程序中参数的处理与数据的整理输出 程序中使用的原始物理量包括粒子的速度和位置，空间网格点上的电势和 电场强度。用到的数组相当大，通过直接认读这些数据来判断结果是不可能的， 我们使用m a t l a b 输出图形来显示数据，以得到简洁直观的认读。 在程序中使用的物理量不一定是我们最关心的，例如粒子模拟用每个粒子 的坐标和速度为原始变量，而我们关心的是粒子的能谱、空间的密度等，因此 建立专门的程序来实现这些量的转换。根据每次运行事例的需要，我们有选择 的输出空间电势分布、空间粒子密度和粒子能谱，并检查程序运行时的能量守 恒性。当事例运行时间较长时，我们在程序结束时输出全部粒子的坐标和速度， 以便需要时可以继续运算。 阱中的电子密度是我们很关心的一个参数。因为我们计算中的粒子云实际 在三维空间中是一个圆环，每个圆环带有一定的电荷。所以统计密度时是先统 计单位面积中的圆环个数，然后再换算成单位体积中的电子数。我们把( cz ) 平面划分成均匀方格，a r = a z = o 2 5 m m 。计算密度也就是计算截面积为方形 的圆环中的电子数，计算中取的单位面积是a r a z = o 2 5 0 2 5 m m2 ：o 0 6 2 5 m m2 。 在迭代法中需要用到的电荷密度的计算与此类似。因为电荷是分配到网格 点上的，故计算时的单位面积是以网格点为中心的方格，间距与计算网格相同， r = a z = 0 5 r a m 。 计算粒子的能谱需要计算粒子的动能和势能。由粒子的二维速度可以直接 得到粒子的动能，计算粒子的势能时除了需要知道粒子的位置外，还需要知道 势场的大小。在我们的计算中，使用的是包括粒子的自恰场与外场的有效势。 中田科掌技术大掌硬士掌脚潍文第二潘宁阱中爿e 中性等瞳卜俸的羲值棚瓤 包括自恰场是因为在中心处粒子电荷产生的自恰场很强，将会强烈的影响粒予 的能量。 2 9 程序流程图 初始化各变量，设定程序参数 ( 网格划分、电磁场、粒子数等) 根据特定粒子的初始分布，给出各 粒子的位置、速度( 或由数据文