基于時域有限差分的微放電仿真算法

張劍鋒，游檢衛，王洪廣，李韻，崔萬照，崔鐵軍,*

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