1、粒子模拟课程设计1.物理模型本课程设计采用 PIC 模拟方法的静电模型，模拟电子枪阴极表面发射电子注，电子在静态场的作用下运动直至稳定的运动过程。由于电子注的环形特征以及圆柱系统具有良好的轴对称性，采用圆柱坐标系来进行理论分析和数值模拟。z 0图 1：多环形电子注在圆柱坐标系统中的运动2.理论分析理论分析包括：电子运动方程的求解、电荷源的求解以及泊松方程的求解三部分。2.1 电荷源的求解课程设计采用的是有限大小粒子模型，且形状、尺寸与所处的网格单元一致，位置由其中心点描述。假设代表宏电子位置为 ，距离其最近的最小网格点标识为 ，(,)rz(,)ijk然后将电荷平均分配到相邻的 8 个网格点上。
2、假设每个宏电子发射电流为 I，则宏粒子的电量为 。按照平均分配的原则，它对每个周围格点的贡献为 。粒子在有限qIt 18qt空间内均匀分布，则利用宏电子所在的网格单元的体积为：2.1.1 21*()2VSZri则每个格点分配的电荷密度大小为 2.1.224()(1)Itzri2.2 求解泊松方程已知空间各个格点上的电荷密度，求解泊松方程，获得各个格点上的点位。在国际单位制下泊松方程形式为：2.3.120在圆柱坐标系下，2.3.2221rrz（1 ） 考虑非轴线上格点的电势求解在 处将 2.3.2 式差分，并设定,ijk2.3.3rz（j=1,2.） 2.3.4(1)()jrj2.3.5124(
3、)cAjj得到：2.3.62, ,1,1,10,1,2()()(1) ()ijkcijkcijkijkijkcijkcijkijkzAAj 注意：差分表达式仅对于非轴向格点适用，即要求 。1（2 ）对于轴线上格点电位的求解2.3.72,1,1,21,11,0()636iiiiiz（3 ）左平面电位固定,为 100V:2.3.8(:,)0rZNV右平面电位为 0：2.3.9(1:,)r圆柱边界面上，除了左右边界点外，利用电势径向连续来求解电势。2.3.10,1,1rriNkik利用建立的电位离散方程组与边界条件，采用迭代法或者直接求解法，即可获得空间各个格点上的电位。2.3.112, ,1,1,
4、0,1,12()()() ()()r rrrr rriNkciNkciNkikcikcikiNkzAAii 其中， 。1r2.3 电子运动方程的求解粒子在场的作用下，运动状态发生变化，具有新的速度和位置。粒子的运动方程为2.4.1dvqEtm2.4.2rt由运动方程可知， （1）先差值求出电子所在位置的场， （ 2）后求解运动方程获得新时刻的速度与位置。（1 ）求粒子所在位置的电场：假设考察的电子位置为 ，距离电子最近的最小网格标识为 ，如图(,)zr(,)ijk3.3 所示，则电子所在位置处电场的分量为,1,1,1,1,1,1,14zijkijkijkijkijkijkijkijkE 2.4
5、.3,1,1,1,11,1,14rijkijijkijkijkijkijkijkEz 2.4.4,1,1,1,1,1,1,11( )( )4ijkijkijijkijkijkijkijkEr 2.4.5内部非轴线点: 1 1, 1, 1,2 222 1, , .204()4()() ()(1 4()ijk ijk ijkijkijk ijkijk ijkz 2.4.6轴线方向点: 2.4.721, 1.2,11,1,0()636kkkkkz（2 ）粒子运动状态更新采用静电模型中的蛙跳格式进行差分处理后,粒子速度在整个时间步长取值,粒子位置与电场力在半个时间步长取值。有：2.4.6112nnzz
6、zqtEm2.4.711 22()nn nrrrtt2.4.8112nn nrqttE3.程序编写此次课程设计，我采用了 c 语言进行编写3.1 主函数代码void main()double *p=(double*)U,*q=(double*)RO;gedianchushi(p,q,Nr,Nc,Nz);/格点初始化for(int t=0;t6000;t+)int n;n=t/10;if(t%10=0)zhurudianzi(lz,n);/第 n 次注入电子monibuchang(lz,n);/第 n 次注入的电子模拟半个时间步长dianhefenpei(lz,q,n,Nr,Nc,Nz); /1
7、.求前 n 次注入电子所处的网格，分配电荷密度dianshidiedai(p,q,Nc,Nz); /2.由电荷密度求电势sudugengxin(lz,q,n,Nc,Nz); /3.由电势求出格点处的电场,由电场更新粒子的速度weizhigengxin(lz,n); /4.由速度更新粒子位置3.2 程序流程图否是否是结束开始格点初始化：电荷密度 0，电势均匀增加， t=0t%10=0?？t注入 5 圈电子并模拟半个时间步长求出前 n 次注入的所有电子所处的网格，分配电荷密度由电荷密度求出电势由电势求出电场强度由电场强度更新粒子速度由粒子速度更新粒子位置t+t6000?3.3 初始化设置（1 ）物
8、理参数包含了各种物理结构形状、尺寸，以及各种物理条件。系统长：0.01m；系统半径：0.005m；入射的每个宏电子电流：-1e-3A最内圈宏电子注入半径：2e-3m，以后均匀递增 2e-4m。（2 ）计算参数轴向网格数 1:Nz,Nz=201;径向网格数 1:Nr,Nr=101;径向空间步长： 5.0*1rRmN3.3.1轴向空间步长： 5.1zL3.3.2角向空间步长： ， 02eN363.3.3时间步长： 13.*ts3.3.4每次注入宏电子 5 圈，共有 36*5=180 个宏电子。注入电子的时间间隔为 10 个时间步长。（3 ）粒子的初始化假设 5 圈电子以同样的速度从左平面入射进系统
9、，且横向分布满足轴对称性。3.3.571.0*/zms3.3.63.3.7rz=0 3.3.8r 和 的布置如图所示满足轴对称性。（4 ）电势初始化左平面固定电位 100V，3.3.9(,1:)10zrN右平面接地，电势为 0，3.3.10(,:)r两平面之间电位均匀变化，且横向等电位， (,1:)(1,:)10(,1:) *()zrrr z zNNi iiN3.3.11（5 ）电荷源初始化全部网格上的电荷源等于 0。（6 ）收敛流程这里只是为了讨论环形电子注在圆柱系统中的运动轨迹，因此对收敛没有要求，作为实例，只计算 6000 个时间步长即可。4.结果讨论4.1 注入五圈电子的截面图图 2：
10、matlab 中绘制的粒子截面图说明：每一圈 36 个，一次注入 5 圈。每圈粒子均匀分布于 360。4.2 粒子轨迹图下图为第一次注入的 5 圈电子运行 6000 个时间步长过程中的粒子轨迹图。图 3：粒子轨迹图说明：由于程序没有调试好，目前只得到这样一个结果。5.心得体会此次课程设计使我收获良多，首先清楚了 PIC 模拟方法的静电模型用于模拟多环电子注在圆柱系统中运动的基本流程，对于电子枪中电子的运动有了深刻清晰的认识。其次也锻炼了自己的编程能力，对 c 语言的使用也更加熟练。编程首先要画系统流程图，这样才能有全局观念，在此要感谢唐茂文同学的指导，是、使我更加清楚了系统流程。编程要做到条理
11、清晰，因此编程时我将其模块化，一个模块一个模块的进行编程，模块编好之后立即进行测试，看是否达到我想要的目的。编程之中的问题，要实现程序模块化，就牵扯到函数调用的问题，函数调用之中最大的问题就是形参和实参，要通过函数调用来改变实参的值就要用到指针，不使用指针就没有涉及编程的精髓，但是使用指针出现的错误较为隐蔽，很难发现，稍有不慎就会使整个程序崩溃。另外一个就是如何用指针来表示三维数组，表示二维数组没有什么难度，但三维数组我之前没有涉及过，是用指针数组还是用一维指针来表示，经过我的尝试，将三维数组强制转换成了一维指针，但是表示方法略微繁琐， （其他表示方法暂时没有尝试成功）致使程序的可读性略微下降。还有就是由于一开始的准备工作没有做充分，细节部分考虑的不够周全，比如轴线电荷密度和电势的处理，致使程序调试不甚理想。另外一点是，课本上的粒子下标是从 1 开始的，而我的程序是从 0 开始的，在公式代入的过程中也遇到了不小的麻烦。今后编程钱一定要充分做好准备，画好流程图，细节部分考虑好，做到磨刀不误砍柴工。
