隨移動窗推進的帶電粒子束團長程傳輸模擬分析*

郝建紅，王 希，張 芳，趙 強，薛碧曦，范杰清，董志偉

(1. 華北電力大學 電氣與電子工程學院， 北京 102206； 2. 北京應用物理與計算數(shù)學研究所, 北京 100094)

通過向太空發(fā)射相對論粒子束，在地面站追蹤其傳輸過程中的參數(shù)變化，可用于追蹤地球磁層磁力線，獲取用于天文物理研究的相關參數(shù)及資料等[1]。此外，探空火箭亦是通過向太空發(fā)射高能粒子束進行長達上萬千米的長距離傳輸來獲取空間科學數(shù)據(jù)。因而，在天體物理、太空主動空間探測等研究領域，帶電粒子束團的長程傳輸是一個基礎研究問題。但以往相關研究大多關注粒子加速器等設備內的帶電粒子米量級及以下的短程傳輸，這種條件下，由于傳輸距離短且磁場的限制抑制了帶電粒子束團的徑向彌散，故常忽略束團尺寸上的變化[2-3]。然而研究帶電粒子束團長程傳輸時，即便在距地35 786 km的地球靜止軌道(Geostationary Earth Orbit，GEO)等準真空環(huán)境下，由于傳輸周期長，其自身初始發(fā)射度和自生空間電荷場導致的徑向尺寸上的膨脹不可忽略。

由于粒子束長程傳輸實驗受限，其研究方法主要采用理論和粒子模擬。對于帶電粒子束團傳輸?shù)难芯浚琀umphries給出了較為完整的帶電粒子束傳輸?shù)陌j方程模型，即K-V方程[4]，對片狀和圓柱形帶電粒子束團在外加電磁場、自生電磁場及初始發(fā)射度等影響下發(fā)生的徑向彌散分別進行了說明。Kochikov等在自由場漂移區(qū)中利用電子脈沖動力學，分別研究了低中高相對論效應的電子束團傳輸4 ns(即1.2 m)的短程傳輸過程，完成了與動力學模擬結果的擬合，但沒有涉及更長距離的傳輸計算[5]。相比帶電粒子束團的短程傳輸，帶電粒子束團的長程傳輸模擬由于影響因素多、計算困難大等原因而更難進行。國內關于帶電粒子束傳輸?shù)奈墨I相對較少，文獻[6-8]對帶電粒子束團在真空長程傳輸時自身空間電荷效應導致的徑向擴散進行了理論研究，但未綜合考慮發(fā)射度的影響。李等基于磁透鏡的強流束傳輸，編寫了適用于模擬計算強流束在磁透鏡中傳輸?shù)氖鱾鬏斈M程序，比較分析了不同束團流強條件下的模擬結果，結果表明：束流越強，自生空間電荷力越大，束流包絡膨脹越明顯[9]。因此，研究強流束團的長程傳輸問題時，必須同時考慮帶電粒子束團自生空間電荷效應及初始發(fā)射度對束團徑向尺寸膨脹的影響。

粒子云(Particle In Cell，PIC)模擬方法是研究帶電粒子運動行為的重要工具之一[10]。在研究帶電粒子束團的長程傳輸時，除初始發(fā)射度的彌散效應，其自身空間電荷效應的影響不可忽略。使用PIC方法可以模擬束團內部的自生電磁場，從而實現(xiàn)長程傳輸過程中束團自生空間電荷效應對其徑向尺寸的影響。利用PIC方法進行粒子模擬時，必須對整個計算空間進行離散，為了考察厘米甚至毫米量級的束團徑向膨脹，束團徑向包絡內至少需要3～5個網(wǎng)格分辨，而PIC方法本身要求離散空間網(wǎng)格的縱橫比不能很大(一般小于10)，故空間網(wǎng)格步長一定在亞毫米量級，時間步長則由有限差分或有限體積方法中滿足穩(wěn)定性及收斂性的收斂條件限定。同時為了給束團長程傳輸?shù)膹较蜃杂膳蛎浟粝伦銐虻目臻g，計算空間的徑向尺度也不能很小，故計算千米量級的長程傳輸時，內存占用量與計算量非常巨大。而由于包含實際有效的束團運動的計算空間占比非常小，通常慣用的空間配分的并行計算效率非常低，從而導致帶電粒子束團長程傳輸?shù)腜IC模擬難以進行。將移動窗技術和PIC模擬方法相結合可以有效解決這一問題，移動窗技術選用大小合理的計算窗口與粒子束團同步運動，且僅離散剖分粒子束團現(xiàn)刻的傳輸區(qū)域，計算保留粒子束團現(xiàn)刻的傳輸結果，通過將無數(shù)個短程傳輸窗口連接起來實現(xiàn)千米量級的長程傳輸，大量節(jié)省計算資源及時間，提高計算效率。

本文考慮GEO準真空環(huán)境，忽略地磁場與稀薄等離子體的碰撞，以PIC靜電模型為基礎，結合移動窗技術建立了一種2D3V(空間二維、速度三維)[11]的片狀帶電粒子束團的長程傳輸模型，有效緩解了制約帶電粒子束團長程傳輸模擬的計算資源難題。通過將模擬結果與解析模型K-V方程的計算結果進行對比分析，驗證所建模型的合理性，并利用所建模型實現(xiàn)了高能相對論電子束團長程傳輸(千米量級)過程的模擬與分析。

1 帶電粒子束團長程傳輸模型

基于PIC模擬技術中時間與空間步長的匹配關系[12]，建立帶電粒子束團長程傳輸模型時，移動窗前進速度、粒子速度、網(wǎng)格大小、時間步長之間需要相互滿足匹配條件。

真空環(huán)境下，帶電粒子束團長程傳輸過程中徑向的自由彌散由自生靜電場驅動，同時伴隨感生磁場的徑向箍縮效應。引入移動窗技術，設置計算窗口以接近光速沿傳輸方向前進時，諸多因素如電磁信號邊界反射、宏粒子數(shù)漲落、空間電荷加載和發(fā)射效應等都會產(chǎn)生噪聲以破壞電磁場環(huán)境，影響傳輸過程。為克服上述問題并降低噪聲干擾，利用PIC靜電模型中電、磁現(xiàn)象退耦的特點，將靜電模型與外加磁場相結合來綜合體現(xiàn)自生電磁場效應的影響，同時通過設置移動窗邊界條件來凈化計算空間的電磁場環(huán)境。

1.1 粒子束團建模

研究帶電粒子束團的長程傳輸過程時，若采用三維帶電粒子束團模型進行模擬，即使結合移動窗技術推進傳輸過程，仍對計算資源提出較高要求。為合理節(jié)省計算資源并提高計算效率，選用2D3V(空間二維、速度三維)坐標系，將z方向默認為單位長度的周期性邊界，建立空間二維片狀帶電粒子束團，如圖1所示。采用二維片狀束團模型，空間上簡化了z向的網(wǎng)格剖分，顯著提高了計算效率；且由于速度與場量均是三維分布，故仍可實現(xiàn)帶電粒子束團長程傳輸過程中z向感生磁場的徑向箍縮效應。

建立帶電粒子束團長程傳輸模型時，空間網(wǎng)格剖分需綜合考慮因果律的限制及移動窗前進速度與束團速度的匹配。傳輸過程中束團沿x軸向傳輸，其軸向(x)尺寸應遠大于徑向(y)尺寸，帶電粒子束團電荷在(x,y)空間中呈均勻分布。

圖1 粒子束團模型Fig.1 Particle beam model

1.2 移動窗技術

移動窗技術通常應用在基于恒定速度傳輸?shù)募す饷}沖或粒子束的相關研究領域。其允許模擬窗口以最高接近光速的速度沿所選方向移動，此功能每時刻僅對包含激光脈沖或粒子束實際運動的有效計算空間進行計算和記錄，最終通過連接無數(shù)分時刻來還原整體模擬，可以顯著減小計算空間尺寸且不影響整體物理傳輸過程。應用于長程傳輸可以節(jié)省大量用于計算空白相空間的計算資源，顯著提高計算效率，縮短計算所需時間。

設置移動窗與束團以相同速度沿軸向運動，故考慮粒子相對論效應時還需進行坐標和速度的洛倫茲變換[13]。移動窗計算空間的大小及邊界條件的設置均需保證帶電粒子束團長程傳輸?shù)目蓪崿F(xiàn)性。真空環(huán)境下，模擬帶電粒子束團的長程傳輸時，除初始發(fā)射度的影響外只需考慮自身空間電荷效應，但此環(huán)境下模擬計算易受邊值條件的影響，為消除和降低模擬過程中的這個影響，將移動窗傳輸方向兩邊界設置為第一類邊界條件，使其模擬環(huán)境更接近于無限大真空。同時在傳輸過程中，特別是長程傳輸后程，應盡可能使束團位于移動窗的中央位置，以保證上述邊界條件成立。

1.3 束團初始發(fā)射度

粒子束團初始發(fā)射度反映粒子束的質量特征，用來表征粒子束團自身匯聚或發(fā)散的品質，較小的橫向發(fā)射度對應于良好準直的小角度發(fā)散束。狹義上，發(fā)射度指帶電粒子在相空間中所占的體積(三維)或面積(二維)，如式(1)所示。

(1)

在模擬中，將宏粒子在二維空間(x,y)中按照邊緣發(fā)射度的定義來抽樣，即

(2)

1.4 束團自生電磁場建模

帶電粒子束團在傳輸過程中，自生電場的排斥力與感生磁場的箍縮力共同作用于束團，這會導致其徑向尺寸的變化[14]。

PIC靜電模型將束團內部宏粒子攜帶的電荷離散至空間網(wǎng)格節(jié)點上求得各節(jié)點處的電荷密度ρ，通過求解真空中的泊松方程，如式(3)所示，得到各頂點靜電勢Φ與網(wǎng)格各邊處的徑向電場分量Ey。

(3)

帶電粒子束團前進過程中形成的電流會在束團周圍感生出上下(y方向)對稱的z向磁場Bz，從而對束團徑向產(chǎn)生箍縮作用。束團傳輸過程中自生電場與感生磁場的數(shù)值關系由式(4)給出[4]。

(4)

模擬帶電粒子束團的長程傳輸時，采用PIC靜電模型求解出自生電場后利用式(4)計算出感生磁場的大小，以外加場的形式使其作用于束團傳輸模擬。此方法區(qū)別于電磁模型[15]，可以有效避免高頻噪聲、振蕩的激勵，節(jié)省此類高頻反應所需的計算空間，提高計算效率。

2 模擬結果驗證

K-V方程是求解帶電粒子束團傳輸?shù)睦碚摻馕龇椒ǎ竟?jié)通過將PIC模擬得到的帶電粒子束團長程傳輸徑向膨脹特性與K-V方程的解析計算結果相對比，以此驗證靜電粒子模擬與移動窗技術相結合的技術手段的可行性及所建模型的合理性。

2.1 K-V方程

Humphries在文獻[4]中提供了旁軸近似情況下的包絡方程，如式(5)所示，描述了片狀帶電粒子束團在傳輸過程中徑向尺寸均方根的變化。

(5)

(6)

其中，R為束團徑向尺寸的均方根，R″為R關于x坐標的二階偏導數(shù)，q為粒子帶電量，J為束團沿z方向的線電流密度，ε0為真空介電常數(shù)，m0為粒子靜止質量，γ和β為相對論因子。

式(5)等號右邊第一項表示非彈性碰撞和軔致輻射所引起的能量損失導致的束團徑向膨脹；第二、第三項表示外加電場和磁場所引起的徑向膨脹或箍縮效應；第四項與第五項則分別表示自生空間電荷效應與束團初始發(fā)射度對徑向尺寸的影響[2]。考慮GEO準真空環(huán)境下無外加電磁場、無粒子碰撞與能損的粒子束團傳輸問題，且認為束團整體始終沿x軸以恒定速度傳輸，故均無須考慮式(5)等號右邊第一、第二、第三項。本模型對應的K-V方程最終形式為：

(7)

2.2 模擬結果與解析結果對比分析

如1.1節(jié)所述，采用電荷均勻分布的空間二維片狀束團模型進行長程傳輸模擬。下面利用結合移動窗技術的PIC方法和K-V方程分別模擬和計算質子與電子兩種典型帶電粒子束團在真空環(huán)境下的長程傳輸過程，對比了兩種方法得到的束團徑向自由膨脹特性。

帶電粒子的相對論效應與其質量與能量均相關，如表1所示，質子質量較大，中高能情況下其相對論因子γ近似為1，屬于非相對論束團；但電子質量較小，故其相對論效應不可忽略。電子的長程傳輸模擬中，移動窗以接近光速前進，故結合移動窗技術研究電子束團傳輸特性時，需進行洛倫茲變換。

表1 典型帶電粒子參數(shù)對比

圖2是初始能量為1 MeV，電流分別為1 mA和10 mA的質子束團長程傳輸?shù)挠嬎憬Y果。圖2縱坐標R/R0為束團徑向尺寸相對初始尺寸的膨脹倍數(shù)。結果表明，不同電流強度下，質子束團的PIC模擬結果與K-V方程的計算結果都吻合得較好。圖2(a)中，由于電流強度較小，PIC模擬中，電荷均勻分布的質子束團在傳輸至約1.7 km時，自身空間電荷力的作用使其徑向膨脹至約初始的65倍(對于束流邊緣的粒子而言，這種相互作用本質上會逐漸發(fā)展為非線性)，且PIC模擬結果與K-V方程(假設束團z方向線電流密度恒定)的計算結果在徑向尺寸膨脹接近65倍之前幾乎重合。但由于式(7)默認式(6)中的線電流密度J恒定，而PIC傳輸模擬過程中J會由于帶電粒子束團各個尺寸上的膨脹而減小，故隨傳輸距離增加，K-V方程解析結果將比PIC模擬結果大。如圖2(b)所示，當束團其余參數(shù)均相同僅電流強度擴大10倍(I=10 mA)時，束團傳輸至1 km，徑向尺寸膨脹至初始尺寸約60倍后K-V方程解析結果大于PIC模擬結果。對比圖2(a)、圖2(b)并結合文獻[9]可知，相同能量下帶電粒子束團的初始電流強度越大，其傳輸過程中的徑向膨脹效應越明顯，式(6)中J的變化速度越快，解析結果的偏大效應就會更早出現(xiàn)且更明顯。

(a) E=1 MeV，I=1 mA

圖3(a)是E=1 MeV、I=100 mA的電子束團在兩種計算方法下的長程傳輸結果的對比，電子束團呈現(xiàn)與質子相似的演化規(guī)律。但傳輸過程中由于電子束團的相對論效應不可忽略，需要進行洛倫茲變換。故在進行坐標和速度變換過程中的某些近似條件會導致電子束團的PIC模擬結果與解析結果略有偏差，但兩計算結果的偏差相較于其整體彌散所占的比例非常小，仍可認為兩計算結果具有較好的擬合度。

圖3(b)給出了E=10 MeV、I=100 mA的中能電子束團在兩種方法下的長程傳輸?shù)挠嬎憬Y果對比，其膨脹特性的變化仍呈相似規(guī)律。對比圖3(a)、圖3(b)可知，相同電流強度時電子束團的初始能量越高，傳輸相同距離時其徑向的膨脹程度越弱，因此傳輸至更遠的距離才出現(xiàn)式(6)中線電流密度J恒定的假設所帶來的誤差，故此模型更適宜用于高能強相對論電子束團的長程傳輸過程模擬。

(a) E=1 MeV, I=100 mA

3 束團傳輸過程及結果

上文驗證了移動窗技術與PIC模擬相結合的方法的可行性及合理性，并指出此模型更適宜用于高能強相對論電子束團長程傳輸。故本節(jié)選取初始能量E=100 MeV的強相對論電子束團在真空環(huán)境下的長程傳輸過程進行分析。PIC模擬中參數(shù)的選擇主要由數(shù)值收斂條件限定，且隨束團的γ增大，靜止坐標系下經(jīng)洛倫茲變換的收斂限制更加嚴格，故強相對論電子束團傳輸模型的參數(shù)選擇需要綜合更多因素，即要同時兼顧收斂條件與傳輸模擬的計算效率，以保證千米量級以上傳輸?shù)目蓪崿F(xiàn)性。由式(6)和式(7)的分析可知，電子束團的初始能量越高，空間電荷效應的影響與發(fā)射度的影響相比就越小，為了更好地體現(xiàn)和分析兩者共同作用下對傳輸過程的影響，模擬選用束團電流為50 A。

傳輸過程中由于內部各粒子之間的排斥力，帶電粒子束團會在所有尺寸上變寬，圖4為上述高能電子束團傳輸?shù)? km處的移動窗口圖。對比圖1可以看出，束團徑向存在明顯的膨脹。其隨傳輸距離的變化如圖5所示，當傳輸距離在一定范圍內時，束團徑向尺寸的膨脹倍數(shù)與傳輸距離呈近似指數(shù)變化規(guī)律；隨傳輸距離進一步增大，膨脹倍數(shù)則與傳輸距離呈近似線性變化規(guī)律，此結果與理論K-V方程的計算結果吻合。

圖4 束團傳輸5 km處的彌散情況Fig.4 Dispersion at 5 km of cluster transmission

圖5 電子束團徑向膨脹隨傳輸距離的變化曲線Fig.5 Curve of radial expansion of the electron beam varies with the propagation distance

束團前進過程中，末端的電子被正向電子減速，而首端的電子被正向電子加速，導致內部的速度分布發(fā)生變化，引起束團軸向尺寸上的膨脹，如圖4所示。圖6(a)給出了上述高能電子束團傳輸至1 km時內部粒子軸向相對論速度λux的分布情況。且隨傳輸距離增加，這種速度變化逐漸累加，當束團傳輸至5 km時，其軸向相對論速度分布的展寬明顯增加，如圖6(b)所示。

(a) x=1 km

關于傳輸過程中自生電荷場的研究，1.4節(jié)中提到自生電場與磁場的大小存在固定線性關系，呈相同變化規(guī)律。由于真空環(huán)境下忽略外場與粒子碰撞，故束團傳輸過程中，自生場在軸向與徑向均高度對稱。但由于強相對論粒子束團的傳輸涉及洛倫茲變換，故變換后靜止坐標系下自生電場與磁場的大小在軸向呈輕微前沖分布，而徑向仍然高度對稱。圖7給出了上述高能電子束團傳輸至5 km時自生電場的分布情況。其中，圖7(a)是電場的分布云圖，顯示了軸向的輕微前沖趨勢；圖7(b)選擇束團軸向中心坐標處，即圖7(a)中黑線處做徑向剖面(線)記錄移動窗口內該線上各位置的自生電場值，表明了自生電場在徑向的高度對稱性。

4 結論

在準真空環(huán)境下帶電粒子束團長程傳輸?shù)倪^程中，束團由于受初始發(fā)射度和空間電荷效應的影響而發(fā)生各方向的膨脹。但受限于計算資源和計算能力、影響因素難以模擬等因素，帶電粒子束團長程傳輸?shù)哪M較難實現(xiàn)，其膨脹效應過程及變化規(guī)律也較難得到。本文將移動窗技術與PIC靜電模型相結合，建立二維片狀帶電粒子束團模型實現(xiàn)了長程傳輸模擬。取質子和電子兩種典型帶電粒子分別對比了不同電流強度、不同初始能量下PIC方法與K-V方程計算得到的束團徑向自由膨脹特性，兩結果具有較好的擬合度，驗證了本文模型的合理性與正確性。并結合K-V方程進行誤差分析，說明此模型更適宜用于模擬高能強相對論電子束團的長程傳輸。通過模擬初始能量為100 MeV的電子束團的長程傳輸過程，分析了其膨脹效應特性及傳輸過程中束流參數(shù)與場量的變化。為進一步研究真空帶電粒子束團長程傳輸效應與電流強度、束團初始能量、初始包絡半徑等因素之間的規(guī)律及傳輸過程中的電磁場特性、束團內部特性等奠定了有力的基礎。