基于时域有限差分的微放电仿真算法

张剑锋，游检卫，王洪广，李韵，崔万照，崔铁军,*

1.东南大学 毫米波国家重点实验室，南京 210096 2.西安交通大学 电子物理与器件研究所，西安 710049 3.中国空间技术研究院 西安分院 空间微波技术重点实验室，西安 710110

基于时域有限差分的微放电仿真算法

张剑锋1，游检卫1，王洪广2，李韵3，崔万照3，崔铁军1,*

1.东南大学 毫米波国家重点实验室，南京 210096 2.西安交通大学 电子物理与器件研究所，西安 710049 3.中国空间技术研究院 西安分院 空间微波技术重点实验室，西安 710110

为克服传统微放电阈值预测方法建模粗糙、精度低的缺点，提高阀值预测精度和效率，研究了基于时域有限差分的精确微放电阀值预测算法。基于时域有限差分算法和粒子追踪算法，通过时空网格自洽互耦实现复杂结构微波器件微放电阈值准确计算。其中，共形和并行算法是提高微放电仿真精度和效率的关键。文章基于空间蛙跳策略实现了基于时域有限差分的微放电仿真算法，利用算法分析了典型微波器件的微放电阈值，仿真与实验结果吻合良好，误差小于0.3 dB；同时，并行效率最高可达83%，验证了算法的准确性和高效性。

时域有限差分；粒子追踪；微放电；共形；并行

相比于静电放电[1]，卫星通信系统中的微放电现象危害更大，分析更困难。微放电仿真算法的早期研究主要针对简单的平板类电极和传输线，结合实验拟合得到阈值参考曲线，依此指导后续器件设计[2]。然而，现实器件的复杂性使得实际阈值曲线与参考曲线间存在较大误差。为提高设计成功率，一般采取过设计，但这将显著增加产品成本；更严重的是，一旦产品验证中发现微放电现象，必须重新设计，大大延长了产品研制周期，这种不确定性已成为星载有效载荷电性能设计的瓶颈。

为解决上述问题，各国科研机构先后开展了微放电仿真算法研究，其中ESA开发出的MEST[3]和FAST3D[4]是该领域的标志性研究成果。MEST和FAST3D都是基于计算电磁学相关算法，阈值预测精度较高，适用范围较广，为实际工程中微波器件的微放电阈值预测提供了有效的分析手段，缩短了产品设计周期，提高了设计成功率。

相比其他星载技术，微放电仿真算法研究起步较晚，不仅存在诸多待研问题，而且已有研究成果也具有较大的提升和完善空间。例如，MEST虽然能够对较复杂的微波器件进行微放电阈值快速预测，但它基于平板近似理论仅考虑了器件中发生微放电的局部结构，提高仿真速度的同时牺牲了计算精度，尤其对同轴类结构误差较大。基于积分方程的FAST3D虽然较差分算法精度高，但电磁粒子自洽互耦建模复杂，且对实际工程中普遍存在的填充了非均匀媒质的微波器件求解困难。

因此，一种能够克服现有算法缺点、在保证计算精度的前提下尽可能提高计算速度的微放电仿真算法成为当前的研究焦点。微放电是多物理场耦合问题，涉及电磁学和力学，数值仿真精度取决于电磁和粒子建模准确性。根据计算电磁学理论，时域有限差分算法(Finite Difference Time Domain，FDTD)[5]是求解时域非均匀电磁问题的精确方法，其网格形式与主流粒子建模方法——粒子追踪算法(Particle In Cell，PIC)[6]相同，都是直六面体，两者耦合建模精度高；更重要的是，FDTD具有天然并行性，计算效率可籍由高性能计算机硬件得到极大提升，是所有电磁算法中并行效率最高的算法。

基于此，本文研究了基于时域有限差分的微放电仿真算法，包括FDTD和PIC的耦合思想、实现方法以及关键技术。本文实现的算法不仅突破了传统微放电仿真算法的局限，改善了仿真精度，而且利用并行技术大幅提高了计算速度，算法更适于复杂结构微波器件微放电设计。

1 基本理论

如图1所示，微放电是空间微波器件中电子在外加电磁场作用下发生的谐振现象，其中初始电子来自宇宙射线和热发射等，它们在外加电场作用下加速碰撞金属腔下壁(见图1(a))并产生二次电子(见图1(b))，若此时电场恰好反向，这些二次电子将在外加电场作用下加速运动直至碰撞金属上壁(假设碰撞前电场未反向)，此时将产生更多的二次电子(见图1(c))，电场此时若再次反向，上壁新产生的二次电子将重复刚刚的过程，即加速碰撞下壁进而产生更多的二次电子，该过程将一直持续下去，直至发生雪崩效应→微放电或随时间增加二次电子数目又逐渐减少→不发生微放电。

理论分析可知，微放电发生与否主要取决于器件结构、传输信号波形和输入功率。微放电阈值是指在确定的激励条件下能够诱发微放电的最小输入功率。显然，微波器件内部电磁场的准确计算是微放电阈值准确预测的关键。实际工程中，虽然微放电主要发生在器件中的类平板电容处，但器件的整体复杂性使这些位置的电磁场与平板波导差别很大，利用平板近似得到的电磁场进行微放电阈值预测误差过大。因此，基于时域麦克斯韦方程的时域有限差分算法成为电磁粒子耦合算法中电磁场求解的优先选择。

2 电磁粒子耦合算法

2.1 时域有限差分算法

FDTD算法是由K.S.Yee提出的一种求解时域麦克斯韦方程组的差分算法[7]。相比于计算电磁学其他主流精确方法，FDTD最大的优点是公式简单、物理概念清晰且无需求解矩阵方程，其核心包括如图2所示的Yee网格和差分公式[8]。

如图2所示，Yee网格中电场和磁场分别位于主网格和对偶网格的棱边中心，并分别由4个磁场和电场包围，是麦克斯韦两个旋度方程的直接体现。Yee网格将求解空间离散化，蛙跳策略则实现时间离散。通过差分近似，FDTD实现了连续算子空间到离散算子空间的转化，进而可计算得到整个求解区域中任意位置上电磁场的时间分布。因此，基于FDTD的电磁粒子耦合算法有助于理解微放电原理并对整个演变过程有更为直观清晰的物理认识。下面讨论如何将带电粒子运动效应耦合入电磁计算，根据电磁理论可知：

(1)

式(1)是全电流定律的微分形式，其中J代表电流源，可具体表示为：

(2)

式中：Jc=σE为传导电流；Jv=ρv为运流电流；Ji为外加电流源。实际工程中，绝大多数问题不涉及自由电荷在真空中运动产生的运流电流，因此传统FDTD算法中并不包含该项贡献。根据第1节基本理论可知，微放电由真空中自由电子谐振运动产生，要描述电磁场与电子的自洽互耦，必须在FDTD迭代计算中加入运流电流贡献，即需将传统差分公式改写为：

(3)

2.2 自洽互耦算法

如上所述，运流电流是实现电磁和粒子自洽互耦的关键，它由微波器件中自由电子运动产生。基于此，本文将采用图3所示的策略进行电磁粒子耦合建模，具体步骤描述为：

1)引入初始电子(采用高斯随机分布)。

2)利用FDTD计算给定微波器件中的电磁场分布，利用插值算法得到器件内所有自由电子所在位置处的电磁场值。

3)利用PIC算法追踪电子运动轨迹并判断其在本Δt时间内是否穿过了器件内金属壁。如果穿过，则借助二次电子发射模型计算新产生的二次电子状态(位置和初速度)，否则进入步骤4。

4)籍由PIC算法计算电子在Δt时间段内的运动轨迹，并将其产生的运流电流加权分配到其在Δt时间段内穿过的所有网格的棱边中心。

5)重复步骤2～4直至达到最大迭代时间。

分析可知，本文采用的电磁粒子耦合建模方法理论上未进行任何近似，实现了电磁和粒子底层自洽互耦，因此能够准确模拟微放电整个过程。需要强调的是，耦合建模的关键在于利用FDTD算法的蛙跳格式，通过时间离散采样实现电子与电磁场相互作用，真实模拟空间中自由电子在微波器件中的微放电过程。

3 电磁粒子共形算法

工程中，星载微波器件发生微放电的位置主要是类平板结构和类同轴传输线结构。对于前者，Yee网格阶梯近似仅影响电磁计算精度，且误差相对较小；对于后者，后续分析可知，阶梯近似对电磁计算和粒子模拟都会造成较大误差，尤其对粒子模拟，加密网格并不能减小误差。因此，研究适用于电磁粒子耦合算法的共形算法是提高复杂结构微波器件微放电仿真精度的关键。

图4(a)是典型多区域问题，其中区域1填充金属，边界为曲面。通过判断网格中心位于哪一区域，曲面边界将近似为图中粗红线所示的折线边界，此即传统阶梯剖分。阶梯近似剖分算法简单，但基于此的电磁计算误差较大。为提高电磁计算精度，在综合考虑计算效率和稳定性后，本文选择了由DeyS.和MittraR.提出的共形算法[10]，即对图4(b)所示的共形网格，磁场计算公式修正为：

(4)

如前所述，即使不引入共形算法，FDTD的计算精度也可通过加密网格得以提高；然而，

非共形剖分将对PIC算法造成很大误差。根据二次电子发射理论[11]，电子入射到金属表面时，金属内部将激发出二次电子，其出射状态取决于入射角度θinc和入射能量Vi，而θinc的确定涉及入射位置外法向。如图5所示，当电子由上向下入射时，对图5(a)所示的阶梯边界，入射角均为0°；而对图5(b)的共形边界，数值入射角与真实入射角可近似认为一致。根据相关理论，θinc的计算误差将严重影响新生二次电子的出射状态，进而大大降低微放电仿真精度。显然，阶梯近似对粒子追踪精度的影响远大于对电磁计算精度的影响，且无法通过加密网格减小误差。因此，粒子共形是耦合算法的必要技术，引入共形边界可提高θinc的计算精度，修正二次电子出射状态，提高后续微放电数值仿真精度。

综上所述，共形算法是高精度微放电阈值预测的关键，而曲面边界共形剖分是算法实现的前提条件。图6是球体目标直六面体网格剖分结果，由图6(b)(d)可知，曲面边界阶梯剖分误差大，虽然增大剖分密度可减小误差，但计算量将急剧增加；同时，如前所述，该方法无法减小由此导致的粒子追踪误差。为此，本文研究并实现了共形剖分算法，其中的关键技术是奇异点处理和曲面边界与直六面体网格求交投影算法[9]。图6(c)(e)是利用共形剖分算法得到的剖分结果。由图可见，共形剖分几乎能够完全重构原始曲面边界，不仅几何剖分精度远高于阶梯剖分，而且边界信息的重建为实现电磁和粒子共形算法提供了必要条件。

4 电磁粒子并行算法

实际工程中，典型微波器件尺寸可达几十个波长，网格数目将达上亿甚至百亿量级；同时，微放电数值仿真中后期粒子数将达百万量级，此时电磁和粒子自洽互耦计算极端耗时，严重影响微放电设计效率，研究加速算法是解决该问题的有效途径。当前，随着64位计算机系统和多核CPU的普及，基于OpenMP的多核并行算法以其高性价比越来越受到研究人员的重视。基于此，本文研究并实现了基于OpenMP的电磁粒子耦合并行加速算法，大大提高了微放电分析速度。

如图7所示，蓝色粗实线代表FDTD或PIC的0代主进程，程序运行至位置①时，进入嵌套并行区域，0代主进程分解为由绿色点划线表示的多个1代子进程，每个1代子进程由CPU的独立核心执行计算任务；所有1代子进程进入位置②后，再次分解成黑色虚线表示的2代子进程；运行到位置③，2代子进程进行归并操作，此时程序由所有1代子进程组成；进入位置④后，1代子进程归并为0代子进程，跳出嵌套并行区域；运行到位置⑤进入并行区，分解生成1代子进程，并行运算至位置⑥，重新归并为0代子进程。基于FDTD和PIC的算法公式可知，上述并行策略可用于FDTD差分迭代、S参数计算和粒子轨迹追踪。需要注意的是：由于网格和粒子数远大于CPU核心数，并行算法将涉及嵌套并行，进程会派生多个线程，因此，同时占用所有CPU核心并不能达到最高并行效率。

5 实验验证

为验证电磁粒子耦合算法对复杂结构微波器件微放电阈值预测的准确性和高效性，对图8所示的脊波导滤波器进行微放电数值仿真和实验测量。脊波导滤波器尺寸为：20mm×10mm×109mm，波导壁镀银处理，工作频率11.5GHz，TE10模激励。

为确保实验结果的可靠性，中国空间技术研究院西安分院空间科学与技术实验室对该器件进行了多次微放电实验，取实验结果波动范围不超过10%的3组数据平均值作为最终结果，得到微放电阈值为525W。

如图9所示，利用本文实现的电磁粒子耦合算法仿真得到了不同输入功率下电子数随时间变化曲线，其中二次电子发射采用修正Vaughan模型[11]，参数为：Vmax=165eV，δmax=2.22，V0=12.5 eV，ks=1。仿真时间35ns，初始宏粒子数100个(每宏粒子包含100个真实电子)。由图可知，微放电阈值约为490W，与实验结果相比，仿真误差小于0.3dB，验证了算法的准确性。

为验证OpenMP并行技术对微放电仿真速度的提升效果，定义并行效率η，即

(5)

式中：Ts为串行算法运行时间；Tp为并行算法使用P个CPU物理核心时的运行时间。

测试平台采用8核CPU、主频3.0GHz的惠普工作站，仿真器件采用图8所示的脊波导。表1给出了不同网格数目下串行(单核)和并行(4核)的计算时间和并行效率。由表可知，相比串行算法，并行算法最多可节省70%的计算时间，并行效率随未知数增多由64%提高至83%。显然，基于OpenMP的并行算法对微放电仿真效率改善明显，是有效的加速方法。需要说明的是，并行效率不仅取决于网格数目和进程数，而且受器件类型影响，器件越复杂并行效率越高。同时，并行计算使用的CPU核心数为实际核心数的一半时性价比最高。

表1 并行效率随未知数变化表

6 结束语

本文研究了时域有限差分算法和粒子追踪算法的自洽互耦建模方法，实现了能够对复杂结构微波器件进行微放电仿真的耦合算法。为提高微放电仿真精度，研究了电磁共形和粒子共形算法，通过修正共形网格磁场迭代公式提高了电磁计算精度；同时，消除了阶梯近似导致的边界法向计算误差，使二次电子出射状态更接近真实情况，提高了粒子追踪精度。为提高微放电仿真效率，研究了基于OpenMP的并行算法，计算速度显著提高。最后，利用本文实现的算法对脊波导滤波器进行了微放电数值仿真，仿真结果与实验结果吻合良好，误差小于0.3dB，验证了算法的正确性；同时，与串行算法相比，并行算法计算时间最大可节省70%，验证了算法的高效性。

References)

[1] 王立，秦晓刚，李凯，等. 空间静电放电传导干扰分析方法研究[J]. 中国空间科学技术，2004，24(5)：51-55.

WANG L，QIN X G，LI K，et al. Research on the analysis method of conductor interference induced by space electrostatic discharges[J]. Chinese Space Science and Technology，2004，24(5)：51-55(in Chinese).

[2] KLEMAN K J. Improved RF system for Aladdin[C]//Proceedings of the Particle Accelerator Conference，Washington, D.C.: IEEE，1993：1235-1237.

[3] DE L J, PÉREZ F, ALFONSECA M, et al. Multipactor prediction for on-board spacecraft RF equipment with

theMEST software tool[J]. IEEE Transactions on Plasma Science, 2006, 34(2):476-484.

[4] VICENTE C, MATTES M, WOLK D, et al. FEST3D -a simulation tool for multipactor prediction[C]. International Workshop on Multipactor Corona & Passive Intermodulation in Space Rf Hardware, Noordwijk: ESTEC, 2005: 11-17.

[5] TAFLOVE A,HAGNESS S C.Computational electrodynamics[M].2nd ed.Boston:Artech House,2000.

[6] VAHEDI V, SURENDRA M. A Monte Carlo collision model for the particle-in-cell method: applications to argon and oxygen discharges[J]. Computer Physics Communications, 1995, 87(1): 179-198.

[7] YEE K S. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media[J]. IEEE Trans. Antennas Propag., 1966, 14(3): 302-307.

[8] 葛德彪. 电磁波时域有限差分方法[M]. 西安：西安电子科技大学出版社, 2011.

[9] YOU J W, WANG H G, ZHANG J F, et al. The conformal TDFIT-PIC method using a new extraction of conformal information (ECI) technique[J]. IEEE Transactions on Plasma Science, 2013, 41(11):3099-3108.

[10] DEY S, MITTRA R. A locally conformal finite difference time domain FDTD algorithm for modeling three-dimensional perfectly conducting objects[J]. IEEE Microwave & Guided Wave Letters, 1997, 7(9):273-275.

[11] 游检卫, 张剑锋, 李韵,等. Vaughan二次电子发射模型的深入研究与拓展[J]. 强激光与粒子束, 2013, 25(11):3035-3039.

YOU J W, ZHANG J F, LI Y, et al. Research and extension of Vaughan′s secondary electron emission[J]. High Power Laser and Particle Beams, 2013, 25(11):3035-3039(in Chinese).

(编辑：高珍)

Multipactor simulation algorithm based on finite difference in time domain

ZHANG Jianfeng1，YOU Jianwei1，WANG Hongguang2，LI Yun3，CUI Wanzhao3，CUI Tiejun1,*

1.StateKeyLaboratoryofMillimeterWaves,SoutheastUniversity,Nanjing210096，China2.KeyLaboratoryforPhysicalElectronicsandDevicesoftheMinistryofEducation,Xi′anJiaotongUniversity,Xi′an710049，China3.NationalKeyLaboratoryofScienceandTechnologyonSpaceMicrowave,ChinaAcademyofSpaceTechnology(Xi′an),Xi′an710110，China

Traditional multipactor prediction method based on plane approximation is fast but inaccurate, which cannot meet the electrical performance design requirements when devices become more and more complex. The finite difference time domain algorithm is one of the most popular methods in the computational electromagnetics. The particle tracing algorithm is a numerical method based on the Newton's force law.The two algorithms can be self-consistently combined by the discrete mesh in both time and space domains, which can be used to predict the multipactor breakdown threshold. In addition, the conformal and parallelizing technologies were well studied. Multipactor breakdown of classic microwave device was analyzed by the developed algorithm. The error is 0.3 dB compared with the experiments and the maximum parallel efficiency is 83%, which demonstrate the accuracy and efficiency of the algorithm.

finite difference time domain；particle tracing；multipactor；conformal；parallelization

10.16708/j.cnki.1000-758X.2017.0036

2016-08-31；

2017-03-07；录用日期：2017-03-17；

时间：2017-03-21 15:36:34

http://kns.cnki.net/kcms/detail/11.1859.V.20170321.1536.002.html

国家自然科学基金(61401096)；国家自然科学基金重点项目(U1537211)

张剑锋(1979-)，男，讲师，jfzhang@seu.edu.cn，研究方向为计算电磁学

*通讯作者：崔铁军(1965-)，男，教授，tjcui@seu.edu.cn，研究方向为计算电磁学、人工电磁材料

张剑锋,游检卫,王洪广,等.基于时域有限差分的微放电仿真算法[J].中国空间科学技术,2017,37(2)：89-95.ZHANGJF,YOUJW,WANGHG,etal.Multipactorsimulationalgorithmbasedonfinitedifferenceintimedomain[J].ChineseSpaceScienceandTechnology, 2017,37(2)：89-95(inChinese).

TN102

A

http://zgkj.cast.cn