張劍鋒,游檢衛(wèi),王洪廣,李韻,崔萬照,崔鐵軍,*
1.東南大學(xué) 毫米波國家重點(diǎn)實(shí)驗室,南京 210096 2.西安交通大學(xué) 電子物理與器件研究所,西安 710049 3.中國空間技術(shù)研究院 西安分院 空間微波技術(shù)重點(diǎn)實(shí)驗室,西安 710110
基于時域有限差分的微放電仿真算法
張劍鋒1,游檢衛(wèi)1,王洪廣2,李韻3,崔萬照3,崔鐵軍1,*
1.東南大學(xué) 毫米波國家重點(diǎn)實(shí)驗室,南京 210096 2.西安交通大學(xué) 電子物理與器件研究所,西安 710049 3.中國空間技術(shù)研究院 西安分院 空間微波技術(shù)重點(diǎn)實(shí)驗室,西安 710110
為克服傳統(tǒng)微放電閾值預(yù)測方法建模粗糙、精度低的缺點(diǎn),提高閥值預(yù)測精度和效率,研究了基于時域有限差分的精確微放電閥值預(yù)測算法?;跁r域有限差分算法和粒子追蹤算法,通過時空網(wǎng)格自洽互耦實(shí)現(xiàn)復(fù)雜結(jié)構(gòu)微波器件微放電閾值準(zhǔn)確計算。其中,共形和并行算法是提高微放電仿真精度和效率的關(guān)鍵。文章基于空間蛙跳策略實(shí)現(xiàn)了基于時域有限差分的微放電仿真算法,利用算法分析了典型微波器件的微放電閾值,仿真與實(shí)驗結(jié)果吻合良好,誤差小于0.3 dB;同時,并行效率最高可達(dá)83%,驗證了算法的準(zhǔn)確性和高效性。
時域有限差分;粒子追蹤;微放電;共形;并行
相比于靜電放電[1],衛(wèi)星通信系統(tǒng)中的微放電現(xiàn)象危害更大,分析更困難。微放電仿真算法的早期研究主要針對簡單的平板類電極和傳輸線,結(jié)合實(shí)驗擬合得到閾值參考曲線,依此指導(dǎo)后續(xù)器件設(shè)計[2]。然而,現(xiàn)實(shí)器件的復(fù)雜性使得實(shí)際閾值曲線與參考曲線間存在較大誤差。為提高設(shè)計成功率,一般采取過設(shè)計,但這將顯著增加產(chǎn)品成本;更嚴(yán)重的是,一旦產(chǎn)品驗證中發(fā)現(xiàn)微放電現(xiàn)象,必須重新設(shè)計,大大延長了產(chǎn)品研制周期,這種不確定性已成為星載有效載荷電性能設(shè)計的瓶頸。
為解決上述問題,各國科研機(jī)構(gòu)先后開展了微放電仿真算法研究,其中ESA開發(fā)出的MEST[3]和FAST3D[4]是該領(lǐng)域的標(biāo)志性研究成果。MEST和FAST3D都是基于計算電磁學(xué)相關(guān)算法,閾值預(yù)測精度較高,適用范圍較廣,為實(shí)際工程中微波器件的微放電閾值預(yù)測提供了有效的分析手段,縮短了產(chǎn)品設(shè)計周期,提高了設(shè)計成功率。
相比其他星載技術(shù),微放電仿真算法研究起步較晚,不僅存在諸多待研問題,而且已有研究成果也具有較大的提升和完善空間。例如,MEST雖然能夠?qū)^復(fù)雜的微波器件進(jìn)行微放電閾值快速預(yù)測,但它基于平板近似理論僅考慮了器件中發(fā)生微放電的局部結(jié)構(gòu),提高仿真速度的同時犧牲了計算精度,尤其對同軸類結(jié)構(gòu)誤差較大?;诜e分方程的FAST3D雖然較差分算法精度高,但電磁粒子自洽互耦建模復(fù)雜,且對實(shí)際工程中普遍存在的填充了非均勻媒質(zhì)的微波器件求解困難。
因此,一種能夠克服現(xiàn)有算法缺點(diǎn)、在保證計算精度的前提下盡可能提高計算速度的微放電仿真算法成為當(dāng)前的研究焦點(diǎn)。微放電是多物理場耦合問題,涉及電磁學(xué)和力學(xué),數(shù)值仿真精度取決于電磁和粒子建模準(zhǔn)確性。根據(jù)計算電磁學(xué)理論,時域有限差分算法(Finite Difference Time Domain,F(xiàn)DTD)[5]是求解時域非均勻電磁問題的精確方法,其網(wǎng)格形式與主流粒子建模方法——粒子追蹤算法(Particle In Cell,PIC)[6]相同,都是直六面體,兩者耦合建模精度高;更重要的是,F(xiàn)DTD具有天然并行性,計算效率可籍由高性能計算機(jī)硬件得到極大提升,是所有電磁算法中并行效率最高的算法。
基于此,本文研究了基于時域有限差分的微放電仿真算法,包括FDTD和PIC的耦合思想、實(shí)現(xiàn)方法以及關(guān)鍵技術(shù)。本文實(shí)現(xiàn)的算法不僅突破了傳統(tǒng)微放電仿真算法的局限,改善了仿真精度,而且利用并行技術(shù)大幅提高了計算速度,算法更適于復(fù)雜結(jié)構(gòu)微波器件微放電設(shè)計。
如圖1所示,微放電是空間微波器件中電子在外加電磁場作用下發(fā)生的諧振現(xiàn)象,其中初始電子來自宇宙射線和熱發(fā)射等,它們在外加電場作用下加速碰撞金屬腔下壁(見圖1(a))并產(chǎn)生二次電子(見圖1(b)),若此時電場恰好反向,這些二次電子將在外加電場作用下加速運(yùn)動直至碰撞金屬上壁(假設(shè)碰撞前電場未反向),此時將產(chǎn)生更多的二次電子(見圖1(c)),電場此時若再次反向,上壁新產(chǎn)生的二次電子將重復(fù)剛剛的過程,即加速碰撞下壁進(jìn)而產(chǎn)生更多的二次電子,該過程將一直持續(xù)下去,直至發(fā)生雪崩效應(yīng)→微放電或隨時間增加二次電子數(shù)目又逐漸減少→不發(fā)生微放電。
理論分析可知,微放電發(fā)生與否主要取決于器件結(jié)構(gòu)、傳輸信號波形和輸入功率。微放電閾值是指在確定的激勵條件下能夠誘發(fā)微放電的最小輸入功率。顯然,微波器件內(nèi)部電磁場的準(zhǔn)確計算是微放電閾值準(zhǔn)確預(yù)測的關(guān)鍵。實(shí)際工程中,雖然微放電主要發(fā)生在器件中的類平板電容處,但器件的整體復(fù)雜性使這些位置的電磁場與平板波導(dǎo)差別很大,利用平板近似得到的電磁場進(jìn)行微放電閾值預(yù)測誤差過大。因此,基于時域麥克斯韋方程的時域有限差分算法成為電磁粒子耦合算法中電磁場求解的優(yōu)先選擇。
2.1 時域有限差分算法
FDTD算法是由K.S.Yee提出的一種求解時域麥克斯韋方程組的差分算法[7]。相比于計算電磁學(xué)其他主流精確方法,F(xiàn)DTD最大的優(yōu)點(diǎn)是公式簡單、物理概念清晰且無需求解矩陣方程,其核心包括如圖2所示的Yee網(wǎng)格和差分公式[8]。
如圖2所示,Yee網(wǎng)格中電場和磁場分別位于主網(wǎng)格和對偶網(wǎng)格的棱邊中心,并分別由4個磁場和電場包圍,是麥克斯韋兩個旋度方程的直接體現(xiàn)。Yee網(wǎng)格將求解空間離散化,蛙跳策略則實(shí)現(xiàn)時間離散。通過差分近似,F(xiàn)DTD實(shí)現(xiàn)了連續(xù)算子空間到離散算子空間的轉(zhuǎn)化,進(jìn)而可計算得到整個求解區(qū)域中任意位置上電磁場的時間分布。因此,基于FDTD的電磁粒子耦合算法有助于理解微放電原理并對整個演變過程有更為直觀清晰的物理認(rèn)識。下面討論如何將帶電粒子運(yùn)動效應(yīng)耦合入電磁計算,根據(jù)電磁理論可知:
(1)
式(1)是全電流定律的微分形式,其中J代表電流源,可具體表示為:
(2)
式中:Jc=σE為傳導(dǎo)電流;Jv=ρv為運(yùn)流電流;Ji為外加電流源。實(shí)際工程中,絕大多數(shù)問題不涉及自由電荷在真空中運(yùn)動產(chǎn)生的運(yùn)流電流,因此傳統(tǒng)FDTD算法中并不包含該項貢獻(xiàn)。根據(jù)第1節(jié)基本理論可知,微放電由真空中自由電子諧振運(yùn)動產(chǎn)生,要描述電磁場與電子的自洽互耦,必須在FDTD迭代計算中加入運(yùn)流電流貢獻(xiàn),即需將傳統(tǒng)差分公式改寫為:
(3)
2.2 自洽互耦算法
如上所述,運(yùn)流電流是實(shí)現(xiàn)電磁和粒子自洽互耦的關(guān)鍵,它由微波器件中自由電子運(yùn)動產(chǎn)生?;诖?,本文將采用圖3所示的策略進(jìn)行電磁粒子耦合建模,具體步驟描述為:
1)引入初始電子(采用高斯隨機(jī)分布)。
2)利用FDTD計算給定微波器件中的電磁場分布,利用插值算法得到器件內(nèi)所有自由電子所在位置處的電磁場值。
3)利用PIC算法追蹤電子運(yùn)動軌跡并判斷其在本Δt時間內(nèi)是否穿過了器件內(nèi)金屬壁。如果穿過,則借助二次電子發(fā)射模型計算新產(chǎn)生的二次電子狀態(tài)(位置和初速度),否則進(jìn)入步驟4。
4)籍由PIC算法計算電子在Δt時間段內(nèi)的運(yùn)動軌跡,并將其產(chǎn)生的運(yùn)流電流加權(quán)分配到其在Δt時間段內(nèi)穿過的所有網(wǎng)格的棱邊中心。
5)重復(fù)步驟2~4直至達(dá)到最大迭代時間。
分析可知,本文采用的電磁粒子耦合建模方法理論上未進(jìn)行任何近似,實(shí)現(xiàn)了電磁和粒子底層自洽互耦,因此能夠準(zhǔn)確模擬微放電整個過程。需要強(qiáng)調(diào)的是,耦合建模的關(guān)鍵在于利用FDTD算法的蛙跳格式,通過時間離散采樣實(shí)現(xiàn)電子與電磁場相互作用,真實(shí)模擬空間中自由電子在微波器件中的微放電過程。
工程中,星載微波器件發(fā)生微放電的位置主要是類平板結(jié)構(gòu)和類同軸傳輸線結(jié)構(gòu)。對于前者,Yee網(wǎng)格階梯近似僅影響電磁計算精度,且誤差相對較??;對于后者,后續(xù)分析可知,階梯近似對電磁計算和粒子模擬都會造成較大誤差,尤其對粒子模擬,加密網(wǎng)格并不能減小誤差。因此,研究適用于電磁粒子耦合算法的共形算法是提高復(fù)雜結(jié)構(gòu)微波器件微放電仿真精度的關(guān)鍵。
圖4(a)是典型多區(qū)域問題,其中區(qū)域1填充金屬,邊界為曲面。通過判斷網(wǎng)格中心位于哪一區(qū)域,曲面邊界將近似為圖中粗紅線所示的折線邊界,此即傳統(tǒng)階梯剖分。階梯近似剖分算法簡單,但基于此的電磁計算誤差較大。為提高電磁計算精度,在綜合考慮計算效率和穩(wěn)定性后,本文選擇了由DeyS.和MittraR.提出的共形算法[10],即對圖4(b)所示的共形網(wǎng)格,磁場計算公式修正為:
(4)
如前所述,即使不引入共形算法,F(xiàn)DTD的計算精度也可通過加密網(wǎng)格得以提高;然而,
非共形剖分將對PIC算法造成很大誤差。根據(jù)二次電子發(fā)射理論[11],電子入射到金屬表面時,金屬內(nèi)部將激發(fā)出二次電子,其出射狀態(tài)取決于入射角度θinc和入射能量Vi,而θinc的確定涉及入射位置外法向。如圖5所示,當(dāng)電子由上向下入射時,對圖5(a)所示的階梯邊界,入射角均為0°;而對圖5(b)的共形邊界,數(shù)值入射角與真實(shí)入射角可近似認(rèn)為一致。根據(jù)相關(guān)理論,θinc的計算誤差將嚴(yán)重影響新生二次電子的出射狀態(tài),進(jìn)而大大降低微放電仿真精度。顯然,階梯近似對粒子追蹤精度的影響遠(yuǎn)大于對電磁計算精度的影響,且無法通過加密網(wǎng)格減小誤差。因此,粒子共形是耦合算法的必要技術(shù),引入共形邊界可提高θinc的計算精度,修正二次電子出射狀態(tài),提高后續(xù)微放電數(shù)值仿真精度。
綜上所述,共形算法是高精度微放電閾值預(yù)測的關(guān)鍵,而曲面邊界共形剖分是算法實(shí)現(xiàn)的前提條件。圖6是球體目標(biāo)直六面體網(wǎng)格剖分結(jié)果,由圖6(b)(d)可知,曲面邊界階梯剖分誤差大,雖然增大剖分密度可減小誤差,但計算量將急劇增加;同時,如前所述,該方法無法減小由此導(dǎo)致的粒子追蹤誤差。為此,本文研究并實(shí)現(xiàn)了共形剖分算法,其中的關(guān)鍵技術(shù)是奇異點(diǎn)處理和曲面邊界與直六面體網(wǎng)格求交投影算法[9]。圖6(c)(e)是利用共形剖分算法得到的剖分結(jié)果。由圖可見,共形剖分幾乎能夠完全重構(gòu)原始曲面邊界,不僅幾何剖分精度遠(yuǎn)高于階梯剖分,而且邊界信息的重建為實(shí)現(xiàn)電磁和粒子共形算法提供了必要條件。
實(shí)際工程中,典型微波器件尺寸可達(dá)幾十個波長,網(wǎng)格數(shù)目將達(dá)上億甚至百億量級;同時,微放電數(shù)值仿真中后期粒子數(shù)將達(dá)百萬量級,此時電磁和粒子自洽互耦計算極端耗時,嚴(yán)重影響微放電設(shè)計效率,研究加速算法是解決該問題的有效途徑。當(dāng)前,隨著64位計算機(jī)系統(tǒng)和多核CPU的普及,基于OpenMP的多核并行算法以其高性價比越來越受到研究人員的重視?;诖耍疚难芯坎?shí)現(xiàn)了基于OpenMP的電磁粒子耦合并行加速算法,大大提高了微放電分析速度。
如圖7所示,藍(lán)色粗實(shí)線代表FDTD或PIC的0代主進(jìn)程,程序運(yùn)行至位置①時,進(jìn)入嵌套并行區(qū)域,0代主進(jìn)程分解為由綠色點(diǎn)劃線表示的多個1代子進(jìn)程,每個1代子進(jìn)程由CPU的獨(dú)立核心執(zhí)行計算任務(wù);所有1代子進(jìn)程進(jìn)入位置②后,再次分解成黑色虛線表示的2代子進(jìn)程;運(yùn)行到位置③,2代子進(jìn)程進(jìn)行歸并操作,此時程序由所有1代子進(jìn)程組成;進(jìn)入位置④后,1代子進(jìn)程歸并為0代子進(jìn)程,跳出嵌套并行區(qū)域;運(yùn)行到位置⑤進(jìn)入并行區(qū),分解生成1代子進(jìn)程,并行運(yùn)算至位置⑥,重新歸并為0代子進(jìn)程。基于FDTD和PIC的算法公式可知,上述并行策略可用于FDTD差分迭代、S參數(shù)計算和粒子軌跡追蹤。需要注意的是:由于網(wǎng)格和粒子數(shù)遠(yuǎn)大于CPU核心數(shù),并行算法將涉及嵌套并行,進(jìn)程會派生多個線程,因此,同時占用所有CPU核心并不能達(dá)到最高并行效率。
為驗證電磁粒子耦合算法對復(fù)雜結(jié)構(gòu)微波器件微放電閾值預(yù)測的準(zhǔn)確性和高效性,對圖8所示的脊波導(dǎo)濾波器進(jìn)行微放電數(shù)值仿真和實(shí)驗測量。脊波導(dǎo)濾波器尺寸為:20mm×10mm×109mm,波導(dǎo)壁鍍銀處理,工作頻率11.5GHz,TE10模激勵。
為確保實(shí)驗結(jié)果的可靠性,中國空間技術(shù)研究院西安分院空間科學(xué)與技術(shù)實(shí)驗室對該器件進(jìn)行了多次微放電實(shí)驗,取實(shí)驗結(jié)果波動范圍不超過10%的3組數(shù)據(jù)平均值作為最終結(jié)果,得到微放電閾值為525W。
如圖9所示,利用本文實(shí)現(xiàn)的電磁粒子耦合算法仿真得到了不同輸入功率下電子數(shù)隨時間變化曲線,其中二次電子發(fā)射采用修正Vaughan模型[11],參數(shù)為:Vmax=165eV,δmax=2.22,V0=12.5 eV,ks=1。仿真時間35ns,初始宏粒子數(shù)100個(每宏粒子包含100個真實(shí)電子)。由圖可知,微放電閾值約為490W,與實(shí)驗結(jié)果相比,仿真誤差小于0.3dB,驗證了算法的準(zhǔn)確性。
為驗證OpenMP并行技術(shù)對微放電仿真速度的提升效果,定義并行效率η,即
(5)
式中:Ts為串行算法運(yùn)行時間;Tp為并行算法使用P個CPU物理核心時的運(yùn)行時間。
測試平臺采用8核CPU、主頻3.0GHz的惠普工作站,仿真器件采用圖8所示的脊波導(dǎo)。表1給出了不同網(wǎng)格數(shù)目下串行(單核)和并行(4核)的計算時間和并行效率。由表可知,相比串行算法,并行算法最多可節(jié)省70%的計算時間,并行效率隨未知數(shù)增多由64%提高至83%。顯然,基于OpenMP的并行算法對微放電仿真效率改善明顯,是有效的加速方法。需要說明的是,并行效率不僅取決于網(wǎng)格數(shù)目和進(jìn)程數(shù),而且受器件類型影響,器件越復(fù)雜并行效率越高。同時,并行計算使用的CPU核心數(shù)為實(shí)際核心數(shù)的一半時性價比最高。
表1 并行效率隨未知數(shù)變化表
本文研究了時域有限差分算法和粒子追蹤算法的自洽互耦建模方法,實(shí)現(xiàn)了能夠?qū)?fù)雜結(jié)構(gòu)微波器件進(jìn)行微放電仿真的耦合算法。為提高微放電仿真精度,研究了電磁共形和粒子共形算法,通過修正共形網(wǎng)格磁場迭代公式提高了電磁計算精度;同時,消除了階梯近似導(dǎo)致的邊界法向計算誤差,使二次電子出射狀態(tài)更接近真實(shí)情況,提高了粒子追蹤精度。為提高微放電仿真效率,研究了基于OpenMP的并行算法,計算速度顯著提高。最后,利用本文實(shí)現(xiàn)的算法對脊波導(dǎo)濾波器進(jìn)行了微放電數(shù)值仿真,仿真結(jié)果與實(shí)驗結(jié)果吻合良好,誤差小于0.3dB,驗證了算法的正確性;同時,與串行算法相比,并行算法計算時間最大可節(jié)省70%,驗證了算法的高效性。
References)
[1] 王立,秦曉剛,李凱,等. 空間靜電放電傳導(dǎo)干擾分析方法研究[J]. 中國空間科學(xué)技術(shù),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]. 西安:西安電子科技大學(xué)出版社, 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] 游檢衛(wèi), 張劍鋒, 李韻,等. Vaughan二次電子發(fā)射模型的深入研究與拓展[J]. 強(qiáng)激光與粒子束, 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
國家自然科學(xué)基金(61401096);國家自然科學(xué)基金重點(diǎn)項目(U1537211)
張劍鋒(1979-),男,講師,jfzhang@seu.edu.cn,研究方向為計算電磁學(xué)
*通訊作者:崔鐵軍(1965-),男,教授,tjcui@seu.edu.cn,研究方向為計算電磁學(xué)、人工電磁材料
張劍鋒,游檢衛(wèi),王洪廣,等.基于時域有限差分的微放電仿真算法[J].中國空間科學(xué)技術(shù),2017,37(2):89-95.ZHANGJF,YOUJW,WANGHG,etal.Multipactorsimulationalgorithmbasedonfinitedifferenceintimedomain[J].ChineseSpaceScienceandTechnology, 2017,37(2):89-95(inChinese).
TN102
A
http://zgkj.cast.cn