基于速調管2.5維粒子模型的非線性注波互作用的研究

高冬平 丁耀根 劉濮鯤 杜朝海 趙 鼎

①(中國科學院電子學研究所中國科學院高功率微波源與技術重點實驗室 北京 100190)②(中國科學院研究生院 北京 100039)

1 引言

速調管自1939年問世以來，作為一種高功率、高增益和高效率的微波放大器件，已經在科學研究、國防建設和工業生產等諸多領域得到了廣泛的應用。速調管工作穩定可靠，使用壽命長，在大功率微波應用領域有著不可替代的作用。但是，速調管的工作機理很復雜，實際研究過程中通常依靠多次制作實驗樣管來修正設計，不僅研制周期長，而且耗資巨大。為更加準確地設計速調管，國內外的許多研究機構開發出一系列模擬軟件，比如基于1維(1D)模型的KLY6，JPNDSK；基于2.5維(2.5D)模型的FCI[1]，Arsenal. MSU，中國科學院電子所的LSP[2]，TESLA[3]，和Magic2D[4]；以及基于3維(3D)模型的Magic3D[4,5]。1維模型通常采用圓盤模型。該模型忽略了電子的徑向運動和電子間的徑向作用力，認為每個圓盤中包含的電子具有全同的運動狀態。1D模型能較好地模擬聚焦良好且電子注波動小的速調管。但是對于導流系數較高，電子注波動較大的速調管，電子的徑向運動已經不可忽略，故此時得到的模擬結果與實測值偏離較大。3維粒子模擬程序(PIC：Particle-In-Cell simulation)可以完全考慮粒子在3個方向的運動信息，通過求解全3維的Maxwell方程來推進粒子的運動。理論上講，該方法最接近速調管的實際物理情形。然而，3維粒子模型由于是全時域粒子模擬代碼，需要巨大的計算資源，這嚴重影響該模型的實用性。

常規的速調管具有軸對稱結構，可以忽略電子的角向運動對互作用系統的影響，從而將全3維模型簡化為2.5D(2維位移，3維速度)模型。此外對于導流系數較大、電子注徑向波動較劇烈的情形，2.5D粒子模型能夠準確地考慮電子的徑向運動，很好地處理漂移管壁對電子的截獲，全面準確地分析輸出腔入口和出口位置處電子注的群聚狀態和電子注的運動電流。這些問題都是1維模型所不能處理的。現有的大多數2.5D模型都是基于端口近似的頻域穩態模型，不能研究系統的瞬態變化過程；此外Magic2D采用全數值的時域模型，雖然能夠較全面地模擬軸對稱結構的互作用系統，但是消耗資源較大，運算時間太長，缺乏實用性。本文改進的2.5D模型是一個將端口近似和PIC相結合的時域粒子模型。該模型可以模擬互作用系統隨時間的演化過程，從而有助于分析諧振腔間隙電壓從開始建立至達到穩定的復雜過程；同時，由于采用端口近似模擬高頻諧振腔，減少了高頻場的計算時間，比全數值的Magic2D高效，更具有實用性。

本文將著重介紹用于研究速調管中復雜注波互作用過程的2.5D粒子模型。該模型是首次將端口近似方法和PIC結合起來的時域粒子模型。大量的研究表明端口近似方法能夠比較準確地模擬諧振腔對互作用系統的影響[1?4]，所以該模型具有高效、準確和節省計算資源等優點。在每一個時間步中，該模型可以選擇采用面積權重或體積權重將粒子電荷恰當地分配到空間網格各點上，再根據電流連續性方程將電流分配到網格各邊上；采用時域有限差分方法求解麥克斯韋方程組以精確解出電子注在漂移管中產生的空間電荷場；采用端口近似來模擬高頻諧振腔場，諧振腔場分布可以采用解析場(漂移頭間隙電場分布采用雙曲余弦函數近似表示)或數值場(通過Superfish或CST等軟件仿真求得)兩種形式；最后，采用Boris旋動格式求解洛倫茲力方程來更新粒子的速度和位置，并在柱坐標系下對粒子的速度推進作了修正。C++語言特有的類設計也使得KLY2D具有很強的擴展性，它提供了多種模擬邊界，有金屬邊界、對稱軸邊界和漂移管右端的吸收邊界。本文將在最后給出對具體速調管實例的計算結果與相關討論。

2 物理模型

2.1 運動粒子電荷與電流的PIC分配

電荷PIC分配采用面積加權[6]的算法來處理，具體過程如圖1所示。設該宏粒子所處的位置為xn，格點1所占的面積為a，格點2所占的面積為b，格點3所占的面積為c，格點4所占的面積為d，軸、徑向的網格步長分別為hz，hr，該網格總面積為S=h?h，另外再定義= (xn?)?z0，zr=(xn?xj,k)?r0，其中z0、r0分別為z、r方向的單位向量。

圖1 PIC電荷分配到網格點

根據面積權重分配原則，4個網格點的電荷密度為

在經過一個時間步長dt以后，宏粒子由初始位置A(xn)運動到了終點位置B(xn*)，位移為AB(如圖2所示)，下面僅討論該宏粒子位移處于一個網格之內的情形(該宏粒子位移穿過多個網格的情形可以分開成多段采用疊加原理進行處理)。

另Δw=wn*?wn，=(wn*+wn)/2，根據式(1)的電荷分配原則和電流連續性方程?Q/?t+∫○J?dS=0，對于圖2中節點1有ΔQj,k/Δt ++Ir,j,k+1/2=0，同理對于節點2、3和4可以列出類似的方程，聯立求解這組方程可以獲得相應網格邊電流密度，

圖2 PIC電流分配到網格各邊

這樣在每個時間步推進以后，都可以根據以上分析來確定空間各個網格點的電荷密度以及各網格邊上的電流密度，從而為下一步利用FDTD方法來求解麥克斯韋方程組作好準備。

2.2 粒子注的空間電荷場

本粒子模型在一個時間步中的推進流程如圖3所示：首先從右下角的方框開始，將電荷、電流分別加權分配于各網格點和各網格邊，然后結合一定的邊界條件采用時域有限差分算法來求解完整的麥克斯韋方程組，從而得到各個網格點的空間電荷場，即圖3中的(Esc,Bsc)，再結合高頻諧振腔場(Ecav,Bcav)得到施加于宏粒子的總電磁場，從而得到作用于宏粒子上的力，求解洛倫茲力方程，得到宏粒子新的速度與位置，接下來應用粒子在邊界上的發射、吸收條件(某個邊界作為粒子發射面使用，某個邊界作為吸收邊界使用，這些都可以在輸入文件中定義)。隨著時間步的推進，當一個高頻周期T到來時，再根據Ramo定理來計算諧振腔新的間隙電壓，一直迭代到所有的諧振腔間隙電壓都達到穩定為止。

粒子注在漂移管空間所產生的電磁場即空間電荷場，程序中采用時域有限差分方法[7,8]求解Maxwell旋度方程來獲得空間電荷場隨時間的演化。整個求解過程是在圓柱形漂移管中進行的，忽略了具體諧振腔的形狀，僅僅是一個漂移管模型，高頻諧振腔對粒子注的影響將在后面討論。

圖3 粒子的推進流程

2.3 高頻諧振腔場

速調管的高頻諧振腔采用端口近似(portapproximation)來處理，如圖4所示。

圖4 端口近似的示意圖

若已知諧振腔某模式(本程序中暫時先考慮TM010模)下的電磁場分布形式，如電場分布函數為fz(r, z)和fr(r, z)，磁場分布函數為fθ(r, z)，實際的諧振腔間隙電磁場值可以由諧振腔的間隙電壓求得。本程序中的場形分布函數可以采用解析場和數值場兩種形式來表征。解析場是根據諧振腔的特征結構參數由雙曲余弦級數求和表示，這樣得出來的諧振腔間隙場形與實際情形比較接近。數值場是先根據諧振腔的幾何參數利用第3方軟件，如CST微波工作室或Superfish來計算出諧振腔的諧振頻率，R/Q，Q值，以及場形分布函數，而具體網格點的電、磁場值再通過插值方法來得到。

運動粒子經過諧振腔時將會在諧振腔中激起感應電流，如果間隙之間的距離遠小于電子波長，此時可以略去滯后效應，從能量守恒的觀點出發可以導出感應電流定律(Ramo’s theorem)[9]如下：

其中iind(t)為諧振腔間隙的感應電流，Vgap(t)為間隙電壓，這里的Vgap(t)為縱向電場Ez沿諧振腔軸線的積分，即在式(3)右邊的積分符號內，J(r,t)為電子注的運動電流，Ecav(r,t)為諧振腔間隙電場，二者都是與位置r有關的矢量。

根據等效電路模型有

其中in為感應電流的第n次諧波分量，Zn為第n次諧波腔的間隙阻抗QL為諧振腔的有載品質因數，fr為諧振腔的諧振頻率)，若為基波腔則取n=1，若為2次諧波腔則取n=2。

綜上，先在每個時間步中計算各諧振腔的感應電流，當一個高頻周期結束后再將該周期內的感應電流作傅里葉分解以得到其基波分量和高次諧波分量，再通過式(4)求得新的間隙電壓，直至各腔間隙電壓達到穩定值為止。

2.4 粒子的推進

粒子的推進采用Boris[10]旋動格式實現，粒子的推進包括粒子的速度和位置推進，在Particle-In-Cell 算法中我們采用中心時間差分，粒子位置的更新是從(n?1)Δt到nΔt，粒子速度的更新是從(n?1/2)Δt 到(n+1/2)Δt ，粒子的位置和速度的更新是相差半個時間步長的，這即是傳統的蛙跳算法。

考慮相對論效應，采用粒子的歸一化動量u=γv，粒子的運動方程可寫為

其中m，q，x，v分別為粒子質量、電量、位置和速度；γ為相對論因子；E，B對應于當前粒子位置處的電場強度與磁感應強度，通過插值就可以由網格點上記錄的場值得到其它任意位置處的電磁場分量。

經過推導，粒子速度由時間步n?1/2的n?1/2u更新為n+1/2的n+1/2u推進流程為

式中2t=t? t。粒子的位置推進形式為

其中γn由u?根據來近似確定。

3 計算結果及相關討論

下面給出以中國科學院電子學研究所研制的S波段50 MW高峰值功率速調管為算例的研究結果。該速調管應用于電子直線加速器，工作頻率為2856 MHz，電子注直流電壓為309 kV，電子注電流為376 A，采用均勻磁場聚焦。該速調管由6個諧振腔組成，采用單間隙輸出腔，其余腔也均為單間隙諧振腔，輸入腔與輸出腔間隙中心間的距離為0.49 m。本程序計算時采用的輸入功率為275 W(與實驗時一致)，下面給出計算結果并對其進行分析討論。

圖5為最后一個時間步漂移管中的宏粒子分布，由于計算的需要本文給模擬區域加了一段必要的過渡段，即輸入腔的間隙中心的z軸位置為0.08 m，各間隙中心間的距離均保持不變，縱軸r為宏粒子的徑向位置對漂移管半徑R0的歸一化值，宏粒子是從發射面隨機發射出來的。從圖5中可以看到，在末前腔和輸出腔出現了明顯的電子注群聚塊，有很強的密度調制。

圖6表示在最后一個時間步漂移管空間中宏粒子歸一化軸向動量γβz=γvz/c 隨軸向距離z的分布，可知在輸出腔(其間隙中心對應圖6中的0.57 m處)中大部分電子由于受到間隙電場的減速作用，電子速度銳減，同時實現了從電子動能到微波能量的轉移。

圖7為在最后一個時間步漂移管空間中宏粒子歸一化徑向動量γβr=γvr/c 隨軸向距離z的分布。隨著運動距離的增加，宏粒子徑向動量的零散也就越大，特別是到了輸出腔的位置，由于有很強的徑向電場作用，宏粒子的徑向速度波動很大。

圖8為宏粒子的相軌跡，它反映了宏粒子在依次受到各個諧振腔的調制而逐漸出現群聚的整個變化過程，它可以作為速調管結構及各項參數選擇是否合理的重要參考之一。從圖8中可以很明顯地看出，在該參數結構下宏粒子出現了比較理想的群聚。

圖9為電子注歸一化運動電流的前4次諧波分量隨軸向距離z的變化曲線，隨著軸向距離的增大，電子注的群聚就越明顯，運動電流的幅值也就越大，其中基波分量幅值最大，其大小也決定了速調管效率的高低，二階及以上的諧波分量相對較小。

圖10為同一工作電壓(309.0 kV)下由KLY2D，LSP以及KLY6程序計算出的傳輸特性曲線，可見KLY2D與KLY6的差異比較大，由后者求出的最大輸出功率高于前者的計算結果，而且可以注意到隨著輸入功率的增加KLY6程序預測的輸出功率很快就達到了飽和，其原因是1維模型中電子動量都集中在軸向，高估了速調管中的互作用效率，而由KLY2D計算得到的輸出功率變化要相對平緩很多。從圖10中還可以看出，相同輸入功率下KLY2D的計算結果要略低于LSP的結果，其原因為輸出腔中電子注在軸向和徑向都會出現嚴重的超越和交叉，LSP圓環模型所依賴的層流性假設不再成立，而KLY2D則不存在這一問題。

圖5 r-z平面宏粒子的分布

圖6 宏粒子的軸向動量隨運動距離z的變化

圖7 宏粒子的徑向動量隨運動距離z的變化

圖8 宏粒子的相軌跡

圖9 電子注歸一化運動電流的諧波分量隨運動距離z的變化

圖10 由KLY2D、LSP和KLY6計算的傳輸特性曲線

當輸入功率為275 W時，由程序計算得到的輸出功率為51.5 MW，效率為44.3%，而對應樣管實驗測得的輸出功率為50.9 MW，由此可見，誤差是很小的。

4 結束語

本文所介紹的粒子模型首次將粒子模擬思想與等效電路端口近似相結合。與1維模型相比，它考慮了宏粒子的徑向運動，模擬了大量的宏粒子的運動，這樣更接近真實的物理情形。采用端口近似方法，能夠準確地模擬諧振腔高頻場對注波互作用系統的影響，同時比全數值的3D模型節約運算資源。基于該模型我們開發出了國內首個速調管專用粒子模擬程序KLY2D。其運算時間與通常的粒子模擬程序相比大為縮短，從而極大地提高了數值模擬的效率。此外，本文的2.5D模型又是一個時域的粒子模型，能夠記錄下所有粒子的瞬態信息，可以模擬互作用系統隨時間的演化過程，從而有助于分析諧振腔間隙電壓從開始建立至達到穩定的復雜過程，這是傳統的頻域代碼無法與之相比的。KLY2D運行穩定，計算結果可信度高，與本單位研發的實際管型實驗數據進行對比驗證，結果表明程序與實測值吻合得很好。

此模型還有很多需要完善的地方，接下來我們將考慮注波互作用區域入口處電子注的速度零散，使其更接近真實的物理情形。同時，我們將引入群聚腔為雙間隙耦合腔的模型，使其能夠用于擴展互作用速調管的模擬。

[1] Shintake T. Klystron simulation and design using the Field Charge Interaction (FCI) code[J]. Nuclear Instruments and Methods in Physics Research, 1995, Section A(363): 83-89.

[2] 趙鼎. 速調管非線性注波互作用程序的研究[D]. [博士論文],中國科學院電子學研究所，2007.Zhao D. Research on nonlinear beam-wave interaction program for klystrons[D]. [Ph.D.dissertation], Institute of Electronics, Chinese Academic of Sciences, 2007.

[3] Chernyavskiy I A, Cooke S J, and Vlasov A N, et al.. Parallel simulation of independent beam-tunnels in multiple-beam klystrons using TESLA[J]. IEEE Transactions on Plasma Science, 2008, 36(3): 670-681.

[4] Goplen B, Ludeking L, and Smithe D, et al..User-configurable MAGIC for electromagnetic PIC calculations[J]. Computer Physics Communications, 1995,87(1-2): 54-86.

[5] 鄒峰, 薛謙忠, 劉濮鯤. 大回旋電子注雙磁會切電子槍的數值模擬[J]. 電子與信息學報, 2008, 30(9): 2276-2278.Zou F, Xue Q Z, and Liu P K. Numerical simulation of large orbit gyrotron electron beam double CUSP gun[J]. Journal of Electronics & Information Technology, 2008, 30(9):2276-2278.

[6] Verboncoeur J P, Landon A B, and Glad N T. An object-oriented electromagnetic PIC code[J]. Computer Physics Communications, 1995, (87): 199-211.

[7] Huang Z Y and Pan G W. Universally applicable uniaxial perfect matched layer formulation for explicit and implicit finite difference time domain algorithms[J]. IEEE Transaction on Microwaves, Antennas & Propagation, 2008,2(7): 668-676.

[8] Wu D G, Chen J, and Liu C R. An efficient FDTD method for axially symmetric LWD environments[J]. IEEE Transactions on Geoscience and Remote Sensing, 2008, 46(6):1652-1656.

[9] 丁耀根. 大功率速調管的理論與計算模擬[M]. 北京：國防工業出版社, 2008: 33-43.Ding Y G. Theory and Computer Simulation of High Power Klystron[M]. Beijing: National Defense Industry Press, 2008:33-43.

[10] Birdsall C K and Langdon A B. Plasma Physics Via Computer Simulation[M]. Bristol: Adam Hilger, 1991: 58-63.