太赫茲回旋脈塞的粒子模擬方法研究

劉大剛，王重陽，廖樹清

?

太赫茲回旋脈塞的粒子模擬方法研究

劉大剛1，王重陽1，廖樹清2

(1. 電子科技大學(xué)物理電子學(xué)院 成都 610054 ; 2. 中國工程物理研究院流體物理研究所 四川綿陽 621900)

分析了粒子模擬方法在模擬共形邊界，特別是帶此類邊界的太赫茲回脈塞時(shí)引入計(jì)算誤差的原因。給出了傳統(tǒng)的共形網(wǎng)格差分算法的穩(wěn)定性條件及兩種通用的共形處理方法，并對兩種方法的穩(wěn)定性及準(zhǔn)確性進(jìn)行了分析。提出了一種新型的自適應(yīng)共形網(wǎng)格剖分方法，該方法解決了共形網(wǎng)格的局限性問題。在有自主知識產(chǎn)權(quán)的全電磁粒子模擬軟件CHIPIC上實(shí)現(xiàn)了該方法，并以一個(gè)圓柱波導(dǎo)驗(yàn)證了該方法的正確性。采用該方法模擬了工作于TE26模式的頻率為0.42 THz的回旋脈塞及HE06模式0.22 THz準(zhǔn)光腔回旋脈塞。

共性網(wǎng)格; 截止頻率; 回旋脈塞; 粒子模擬; 太赫茲波

電子回旋脈塞(electron cyclotron resonate maser, ECRM)的研究和發(fā)展，是微波電子學(xué)中的一個(gè)重大突破，推動(dòng)了一大類新型高功率毫米波、亞毫米波器件——回旋管的誕生?；匦芸商钛a(bǔ)傳統(tǒng)微波器件和激光器在毫米波和亞毫米波段的缺口[1]，其波段也可發(fā)展到太赫茲波領(lǐng)域。當(dāng)工作波段位于太赫茲波段時(shí)，微小的工作波長將導(dǎo)致回旋脈塞結(jié)構(gòu)尺寸的進(jìn)一步縮小。因此，太赫茲回旋脈塞的工作性能對結(jié)構(gòu)參數(shù)非常敏感，這對數(shù)值模擬的模型及算法精度提出了很高的要求。

另一方面，粒子模擬方法(PIC)是一種使用計(jì)算機(jī)來模擬電磁場與粒子相互作用問題的數(shù)值模擬方法。從電磁場算法上來看，該方法采用時(shí)域有限差分算法[2](finite difference time domain, FDTD)，F(xiàn)DTD直接從概括電磁場普遍規(guī)律的麥克斯韋旋度方程出發(fā)，將其轉(zhuǎn)換為差分方程組，是在一定體積內(nèi)和一段時(shí)間上對連續(xù)電磁場的數(shù)據(jù)取樣。因此，它是對電磁場問題的最原始、最本質(zhì)、最完備的數(shù)值模擬，具有較廣泛的適用性。從粒子算法上來看，PIC方法從洛倫茲、動(dòng)量、沖量等定律出發(fā)來描述粒子在電磁場中的運(yùn)動(dòng)，其所計(jì)算的粒子運(yùn)動(dòng)行為是全面的，繼而更能反映出實(shí)際的粒子運(yùn)動(dòng)規(guī)律。也正是因?yàn)镻IC方法的以上特點(diǎn)，目前許多電真空微波源的研究者都將其作為數(shù)值模擬的首選方法。然而PIC方法中FDTD算法采用的是YEE網(wǎng)格模型，其在復(fù)雜邊界建模上與實(shí)際結(jié)構(gòu)的差異將可能降低模擬精度(主要體現(xiàn)在階梯近似誤差和網(wǎng)格數(shù)值色散誤差)，這在對結(jié)構(gòu)參數(shù)較敏感的太赫茲回旋脈塞的模擬中是致命的。雖然階梯近似誤差和網(wǎng)格數(shù)值色散誤差隨著網(wǎng)格尺寸的減小而減小，但是網(wǎng)格尺寸的減小會(huì)使在模擬一定物理空間所需的總網(wǎng)格數(shù)急劇增加，以及保證算法穩(wěn)定性所需的時(shí)間步長急劇減小，從而導(dǎo)致計(jì)算機(jī)所需存儲(chǔ)空間和CPU計(jì)算時(shí)間急劇增大。

綜上所述，要對太赫茲回旋脈塞進(jìn)行準(zhǔn)確高效的粒子模擬，必須對PIC方法中的FDTD算法有所改進(jìn)。目前，國內(nèi)外研究者采取的改進(jìn)方法主要有3種途徑：1) 加入并行算法；2) 采用交替隱式差分算法；3) 采用共形方法。并行算法通過多線程協(xié)作的方式來提高運(yùn)算速度，但其缺點(diǎn)是隨著線程數(shù)的增加，并行加速比將趨于飽和[3]。交替隱式差分方法通過放大時(shí)間步長的方法來提高運(yùn)算速度，但時(shí)間步長增長將導(dǎo)致粒子運(yùn)算的不穩(wěn)定性[4]。共形方法[5]能在較大網(wǎng)格尺寸情況下，保證建模的準(zhǔn)確性，從而以較快的運(yùn)算速度得到較準(zhǔn)確的模擬結(jié)果。目前，該方法在計(jì)算電磁散射等問題中得到廣泛應(yīng)用，但在粒子模擬領(lǐng)域使用較少。

1 自適應(yīng)共形網(wǎng)格剖分方法

1.1 矩形網(wǎng)格模型缺陷

圖1所示為矩形網(wǎng)格及共形網(wǎng)格對一斜漸變金屬導(dǎo)體結(jié)構(gòu)的擬合比較，其中黑影部分為導(dǎo)體區(qū)域，從圖可看出當(dāng)網(wǎng)格尺寸較大時(shí)，矩形網(wǎng)格將不能準(zhǔn)確地反映該結(jié)構(gòu)的斜漸變特性。圖2所示為用矩形網(wǎng)格對一回旋脈塞建模的示意圖。實(shí)線是設(shè)計(jì)尺寸，是建模后的互作用區(qū)長度，￠是實(shí)際設(shè)計(jì)的互作用區(qū)長度。很明顯經(jīng)過矩形網(wǎng)格剖分后，回旋脈塞互作用區(qū)長度變長，從而導(dǎo)致模擬出現(xiàn)較大誤差。

a. 矩形網(wǎng)格擬合?????b. 共形網(wǎng)格擬合

圖2 回旋脈塞矩形建模剖面圖

1.2 傳統(tǒng)共形差分公式

如圖3a所示，設(shè)一共形元胞其中黑色部分表示導(dǎo)體所占的元胞區(qū)域，白色表示真空所占元胞區(qū)域。在此元胞中，電場的遞推公式保持不變，仍按照一般FDTD進(jìn)行計(jì)算；而對于其磁場，則需按下式進(jìn)行迭代[5]：

1) 共形元胞在導(dǎo)體外部的面積應(yīng)該大于元胞總面積的5%，即該元胞相對面積(不含導(dǎo)體部分面積與總面積之比)要大于1/20；

2) 共形元胞中的最長相對邊長(不含導(dǎo)體部分)與該元胞相對面積的數(shù)值比應(yīng)小于12。

共形元胞的局限性主要體現(xiàn)在圖3b和圖3c所示的兩種元胞模型中。對于此類元胞，要么不滿足共形網(wǎng)格特性的制約條件1)，要么不滿足制約條件2)，當(dāng)用式(1)計(jì)算此類元胞時(shí)將出現(xiàn)畸變場，從而影響計(jì)算的穩(wěn)定性。處理此類元胞的方法有兩種[3]，即后向加權(quán)平均的修正方法(MCFDTD)及相對面積修正方法(SC-FDTD)。MCFDTD方法的實(shí)質(zhì)是將用式(1)計(jì)算出的當(dāng)前時(shí)間步的電磁場與它上一個(gè)時(shí)間步的電磁場再作平均，以達(dá)到消除畸變場的目的。該方法在一定程度上會(huì)改善計(jì)算的穩(wěn)定性，但需要在迭代時(shí)再對電磁場進(jìn)行額外處理，增加了計(jì)算時(shí)間。另一方面，當(dāng)畸變場過大時(shí)該方法還是不能滿足計(jì)算的穩(wěn)定性要求。SC-FDTD方法是將模擬區(qū)域內(nèi)所有相對面積小于1/6的共形元胞的相對面積都近似為1/6。該方法在一定程度上也會(huì)改善計(jì)算的穩(wěn)定性，但其影響了計(jì)算的準(zhǔn)確性。因?yàn)樵S多滿足共形限定條件的相對面積在1/20～1/6之間的共形元胞的相對面積都被當(dāng)作1/6來處理，從而降低了模擬的準(zhǔn)確性。

1.3 自適應(yīng)共形元胞剖分

圖4 共形邊界擬合方法示意

圖4所示為解決共形元胞限定問題的一種新型的共形元胞剖分方法。圖4a所示為按實(shí)際共形邊界剖分后所得共形元胞示意圖，斜虛線是實(shí)際的共形邊界，C1、C2、C3是按此邊界剖分出的共形元胞，C1元胞將不滿足共形限定條件。圖4b所示為擬合后的效果示意，該方法的思路是在保證共形邊界基本形狀不變的前提下，適當(dāng)增加不滿足共形限定條件的元胞的邊長，從而增加其共形面積以達(dá)到使其滿足共形限定條件的目的。具體作法是：首先判斷共形邊界所在的元胞C1是否滿足共形限定條件1)，如果不滿足，則按比率增加其共形邊長直至滿足共形條件1)。然后判斷是否滿足共形條件2)，如果不滿足則增加短邊共形邊長，直至滿足共形條件2)。從圖4b可看出，經(jīng)過此處理后雖然也會(huì)引入一定的誤差，但該方法是以共形限定條件作為依據(jù)來處理的，避免了SC-FDTD方法所引入的誤差。另外，該方法處理過程主要在剖分元胞階段進(jìn)行，無需在迭代時(shí)增加額外的處理，也就彌補(bǔ)了MCFDTD的缺陷。

1.4 自適應(yīng)共形元胞的模擬驗(yàn)證

本文采用直角坐標(biāo)系，用自適應(yīng)共形元胞及非共形元胞分別對圓柱波導(dǎo)中的TE11模進(jìn)行模擬。圓柱波導(dǎo)半徑10 mm，長20 cm，理論截止頻率c= 8.797 654 GHz，如圖5所示。模擬時(shí)，元胞尺寸大小為1 mm，在波導(dǎo)激勵(lì)源[8]注入峰值功率為5 kW，模式為TE11的電磁波，如圖5c所示。

圖5 圓柱波導(dǎo)結(jié)構(gòu)

圖6a和圖6b為注入頻率為截止頻率的1.1倍時(shí)的共形與非共形輸出功率比較。從圖可看出，由于輸入頻率大于截止頻率，電磁場反射較小，輸出功率接近于注入功率5 kW。

圖6 自適應(yīng)共形與非共形輸出比較

從圖6c和圖6d可看出，當(dāng)注入頻率接近fc時(shí)，自適應(yīng)共形計(jì)算輸出平均功率接近于0，而非共形還有一定輸出。這是因?yàn)榉枪残嗡惴ㄔ跀M合圓柱波導(dǎo)時(shí)是用矩形元胞擬合的，如圖5b所示。這種擬合與自適應(yīng)共形元胞擬合比較，擬合出的圓柱波導(dǎo)半徑比實(shí)際半徑偏差更大，且在圓周橫截面上有更多毛刺所致。這就證明了在相同元胞步長情況下，自適應(yīng)共形算法計(jì)算準(zhǔn)確性比非共形高。

2 自適應(yīng)共形網(wǎng)格模擬驗(yàn)證

2.1 TE06模式0.42 THz回旋脈塞的共形模擬

回旋脈塞的結(jié)構(gòu)如圖7a所示[7]，其中為9 mm，為11 mm，為9 mm，為3°，為2°，為2.2 mm，電子注引導(dǎo)中心半徑為1.15 mm，電子注橫縱速度比為1.5，工作電流1.5 A，電壓20 kV，互作用區(qū)磁場8.11 T。

模擬時(shí)為了保證邊界條件封閉在脈塞左端、右端加入吸收邊界條件，同時(shí)為了引出電子在脈塞左端設(shè)置一環(huán)狀金屬導(dǎo)體作為發(fā)射體，如圖7b所示。

在圖7b中，由于該脈塞工作于太赫茲頻段，其結(jié)構(gòu)尺寸偏小(徑向總長不到3 mm)，所以其徑向及軸向的元胞剖分都很小。另外，為了降低磁場，該器件被設(shè)計(jì)成工作于二次諧波TE26模式上，這不僅要求必須是3維建模，而且要求要有更高的建模準(zhǔn)確性，因?yàn)榻r(shí)誤差過大將會(huì)導(dǎo)致該器件工作點(diǎn)的偏移[9-10]。

圖8a和圖8b所示為模擬所得模式分布圖，圖8c和圖8d所示為模擬所得頻譜及功率圖。從模式分布圖可看出，其工作模式在角向有兩個(gè)周期，徑向有6個(gè)半周期，其模式是TE26模式，從頻譜及功率圖可看出，其輸出平均功率為7.8 KW左右，工作頻率為0.419 4 THz，這與參考文獻(xiàn)[7]中的結(jié)果圖9基本吻合。

圖9 各個(gè)模式下功率圖

2.2 HE06模式0.22 THz準(zhǔn)光腔回旋脈塞的模擬

該器件被設(shè)計(jì)成工作于HE06模式[11]，工作頻率為0.22 THz，輸入功率60 W，放大輸出功率10 KW。腔體結(jié)構(gòu)由上下兩個(gè)半徑為4.52 mm的反射鏡組成[12]，如圖10b所示。反射鏡側(cè)面有兩個(gè)側(cè)弧用于衍射返波，腔體總長57 mm，如圖10a所示。

圖10b和圖10c分別是準(zhǔn)光腔非共形與共形建模結(jié)果比較，由圖可知矩形網(wǎng)格擬合時(shí)會(huì)引入較大的誤差。由于此管工作頻率較高，且屬于放大類器件，此誤差會(huì)導(dǎo)致工作點(diǎn)的偏移，使其工作于振蕩模式，從而干擾放大器的正常工作。

圖11是共形輸出功率與頻譜。由圖11a與圖11b可知，自適應(yīng)共形擬合較好的保證了建模精度，從而使器件穩(wěn)定工作于放大狀態(tài)。

圖12是非共形輸出功率與頻譜，矩形擬合時(shí)，由于誤差的引入使得模式為HE05的返波振蕩出現(xiàn)(圖12a中頻率191.2 GHz處)，從而破壞了放大工作原理，最終導(dǎo)致輸出功率的急劇下降，如圖12b所示。

3 結(jié)束語

粒子模擬是研究粒子與電磁場相互作用的一種重要手段，它可以幫助研究人員認(rèn)識和理解太赫茲回旋脈塞中的一些復(fù)雜的物理問題，也可以幫助研究人員優(yōu)化相應(yīng)器件的設(shè)計(jì)，從而使器件工作性能達(dá)到最優(yōu)。然而，要對帶有共形邊界的太赫茲回旋脈塞進(jìn)行電磁粒子模擬，如果采用標(biāo)準(zhǔn)的粒子模擬方法，要么保證不了必要的運(yùn)算精度，要么使運(yùn)算時(shí)間過長從而導(dǎo)致模擬無法完成。共形網(wǎng)格算法能較好的解決這一類問題，它能在較大的元胞尺寸下保證運(yùn)算的精度，從而在較短的時(shí)間內(nèi)模擬出正確結(jié)果。然而，共形算法亦有其局限性，本文在分析了標(biāo)準(zhǔn)時(shí)域有限差分法對共形邊界模擬所引入計(jì)算誤差的原因，以及兩種通用共形處理方法的局限性基礎(chǔ)上，提出了自適應(yīng)共形剖分方法。該方法能在盡量減少誤差的前提下，保證共形邊界得以較準(zhǔn)確地?cái)M合，從而確保模擬的準(zhǔn)確性。給出了工作于0.42 THz模式為TE26及工作于0.22 THz模式HE06的回旋脈塞算例，算例結(jié)果與文獻(xiàn)給出的實(shí)驗(yàn)結(jié)果及理論推導(dǎo)結(jié)果基本穩(wěn)合，進(jìn)而證明了該方法的有效性。

[1] CHU K R.The electron cycltron maser[J]. Rev Mod Phys, 2004,76: 489-540.

[2] 王秉中. 計(jì)算電磁學(xué)[M]. 北京: 科學(xué)出版社, 2002.

全面標(biāo)準(zhǔn)化管理與服務(wù)業(yè)從業(yè)人員素質(zhì)提升相輔相成，通過實(shí)施與宣貫標(biāo)準(zhǔn)體系，進(jìn)行團(tuán)隊(duì)整體培訓(xùn)和管理經(jīng)驗(yàn)的積累，提高從業(yè)人員的職業(yè)道德、知識水平、服務(wù)技能。多渠道開展服務(wù)業(yè)管理人員和從業(yè)人員的標(biāo)準(zhǔn)化知識培訓(xùn)，鼓勵(lì)高校及科研院所與企業(yè)的聯(lián)合共建工作，在服務(wù)業(yè)示范點(diǎn)設(shè)置標(biāo)準(zhǔn)化創(chuàng)新研究基地，培養(yǎng)既有標(biāo)準(zhǔn)化理論功底又有豐富標(biāo)準(zhǔn)化工作經(jīng)驗(yàn)的人才。

WANG Bing-zhong. Electromagnetics calculation[M]. Beijing: Science Press, 2002.

[3] 廖臣, 劉大剛, 劉盛綱. 三維電磁粒子模擬并行計(jì)算的研究[J]. 物理學(xué)報(bào), 2009, 58(10): 6709-6718.

LIAO Chen, LIU Da-gang, LIU Sheng-gang. Three- dimensional electromagnetic particle-in-cell simulation by parallel computing[J]. Acta Phys Sin, 2009, 58(10): 6709- 6718.

[4] LIU Shao-bin, LIU San-qiu. One-step alternating direction implicit FDTD algorithm[J]. Chin Phys B, 2004,13(11): 01892.

[5] 葛德彪, 閻玉波. 電磁波時(shí)域有限差分法[M]. 西安: 西安電子科技大學(xué)出版社, 2011.

GE De-biao, YAN Yu-bo. Finite-difference time-domain method for electromagnetic waves[M]. Xi￠an: Xidian University Press, 2011.

[7] 袁雪松, 鄢揚(yáng), 傅文杰, 等. 高次諧波太赫茲回旋管的多模工作研究[J]. 紅外與毫米波學(xué)報(bào), 2012, 31(4): 342-347.

YUAN Xue-song, YAN Yang, FU Wen-jie, et al. Multi- model high harmonic operation in a terahertz gyrotron[J]. Infrared Millim Journal of Infrared and Millimeter Waves, 2012, 31(4): 342-347.

[8] 周俊, 祝大軍, 劉大剛, 等. 粒子模擬中波導(dǎo)激勵(lì)源的設(shè)計(jì)與實(shí)現(xiàn)[J]. 強(qiáng)激光與粒子束, 2006, 18(12): 2019-2024.

ZHOU Jun, ZHU Da-jun, LIU Da-gang, et al. Deign and realization of waveguide excitation source in particle-in-cell simulation[J]. High Power Laser and Particle Beam, 2006, 18(12): 2019-2024.

[9] HAN Y, YUAN X S. Study of a grotron oscilltor withcorrugated interaction cavity[J]. Acta Phys Sin, 2012, 61(6): 64102-064102.

[10] 謝鴻全, 劉濮鯤. 等離子體填充波紋波導(dǎo)中低頻模式特性分析[J]. 物理學(xué)報(bào), 2004, 53(9): 3114-3118.

XIE Hong-quan, LIU Pu-kun. Analysis of the characteristics of low-frequency models in a corrugated waveguide filled with plasma[J]. Acta Phys Sin, 2004, 53(9): 3114-3118.

[11] SIRIGIRI J R, SHAPIRO M A, TEMKIN R J. High-power 140-GHz quasioptical gyrotron traveling-wave amplifier [J]. PhysRevtt, 2003, 258302.

[12] 王平. 太赫茲準(zhǔn)光回旋器件的三維數(shù)值模擬仿真[D]. 成都: 電子科技大學(xué), 2014.

WANG Ping. Three-dimensional numerical simulation of terahertz quasi-optical devices[D]. Chengdu: University of Eelectronic Science and Technology of China, 2014.

編 輯 黃 莘

Research of Terahertz Cyclotron Oscillator Particle Simulation

LIU Da-gang1, WANG Chong-yang1, and LIAOShu-qing2

(1. School of Physical Electronics, University of Electronic Science and Technology of China Chengdu 610054; 2. Institute of Fluid Physics, China Academy of Engineering Physics Mianyang Sichuan 621900)

The calculation errors of particle-in-cell(PIC) method at the conformal boundary are analyzed, especially for the conformal boundary at the THz cyclotron oscillator. The condition of stability and two general treatment methods for the conformal grid are put forward. And the stability and accuracy of these two conformal treatment methods are analyzed. Based on these analysis, an adaptive conformal cell subdivision method is developed. With this method, the problems of the limitation of the conformal grid are solved. The method is implemented on the electromagnetic particle simulation software CHIPIC and is verified by simulating a cylindrical waveguide. Finally, using this method, the THz cyclotron maser with TE26 and the quasi optical cavity cyclotron oscillator with HE06 are simulated successfully.

conformal grid; cut-off frequency; cyclotron oscillator; particle-in-cell simulation; terahertz waves

TN12

A

10.3969/j.issn.1001-0548.2017.02.009

2016-01-30；

2016-07-05

中國工程物理研究院科學(xué)技術(shù)發(fā)展基金(2014B0402057)

劉大剛(1973-)，男，博士，副教授，主要從事高功率微波源粒子模擬和太赫茲源模擬等方面的研究.