石黎銘 吳雪科 萬迪 李會東? 樊群超? 王中天 馮灝 王占輝 馬杰

1)(西華大學(xué)理學(xué)院,高性能科學(xué)計算省高校重點實驗室,成都 610039)

2)(山西大學(xué)量子光學(xué)與光量子國家重點實驗室,太原 030006)

1 引 言

盡管托卡馬克被認(rèn)為是未來最具有實用價值的磁約束聚變裝置,但其聚變β值較低、經(jīng)濟效益不明顯,且實現(xiàn)可控核聚變反應(yīng)至今還有很多技術(shù)難題未得到很好的解決.具有較高β值的類磁鏡聚變反應(yīng)裝置在經(jīng)濟和效率上更具優(yōu)勢[1,2].串節(jié)磁鏡[3]、反場箍縮(reversed field pinch,RFP)[4]和場反位形(field reversed configuration,FRC)[5-7]等類磁鏡結(jié)構(gòu)聚變反應(yīng)裝置作為實現(xiàn)聚變能源商業(yè)化的一支潛在力量引起了世界上不少科學(xué)家的研究興趣[1,3,8].

近年來,串節(jié)磁鏡裝置KMAX、多勢阱裝置EMC2、RFP等類磁鏡裝置在約束高溫等離子體方面取得了很大進展[3-14].緊湊型聚變反應(yīng)裝置(compact fusion reactor,CFR)是由洛克希德馬丁(Lockheed Martin)公司根據(jù)“磁鏡約束”原理提出的一種“高β聚變反應(yīng)堆”,其體積較小,具有較大的潛在應(yīng)用價值.CFR裝置由多個線圈構(gòu)成(圖1),中性束由加熱裝置注入,約束磁場由多組同軸電流線圈產(chǎn)生,分別是一個中心線圈,一組內(nèi)部線圈,兩組封裝線圈和一組磁鏡線圈,各組線圈的半徑大小、位置、電流強度以及線圈數(shù)量均可根據(jù)需求適當(dāng)調(diào)節(jié).裝置內(nèi)部線圈的電流方向與其他幾組線圈電流方向相反,在裝置內(nèi)形成會切磁場位形,構(gòu)成了一個邊界附近磁場較強而芯部磁場較弱的磁阱結(jié)構(gòu),使高溫等離子體能夠較好地約束在裝置內(nèi)部[14].當(dāng)軸心附近的帶電粒子向外運動時,周圍的強磁場就會再把它推回來.盡管一開始這種作用力很小,但是粒子偏離軸心越遠(yuǎn),周圍磁場的推力就會越大[14].最新的T4B實驗結(jié)果表明,CFR裝置中的等離子體在高β值條件下具有穩(wěn)定的膨脹率；零維約束模型預(yù)測結(jié)果表明,在T4B裝置中可以通過中性He粒子束在加熱條件下所發(fā)生的電荷交換反應(yīng)向等離子體傳遞能量,加熱等離子體,使等離子體獲得毫秒量級的良好約束,同時裝置內(nèi)會產(chǎn)生高能電子和He離子[15].CFR裝置內(nèi)高能粒子的鞘層損失和軸向尖端損失對高溫等離子體的約束具有重要影響,研究這些高能粒子的約束性能具有重要的實際意義[16].

圖1 Lockheed Martin緊湊型聚變反應(yīng)磁鏡裝置 (a)構(gòu)造；(b)等離子體分布Fig.1.(a)The structure of CFR machine；(b)plasma distribution in CFR (https：//en.wikipedia.org/wiki/Lockheed_Martin_Compact_Fusion_Reactor).

徑向電場Er在磁約束聚變等離子體的研究中具有重要作用[17-23].1982年在ASDEX托卡馬克磁約束聚變中首次發(fā)現(xiàn)徑向電場在L模到H模轉(zhuǎn)化過程中的關(guān)鍵性作用,徑向電場的突然變化引起了L模到H模的轉(zhuǎn)化,從而抑制了湍流波動,由反常輸運過程控制的等離子體約束在H模下得到了改善[17-19].近年來,有關(guān)徑向電場的研究得到了迅速發(fā)展：理論研究表明徑向電場在形成輸運壘和L-H轉(zhuǎn)換機制中可能起著決定性作用[20],實驗上也證實了徑向電場和轉(zhuǎn)換機制的存在[21].此外,徑向電場對托卡馬克等離子體中粒子的“香蕉”軌道中心漂移、粒子軌道和損失錐邊界等都會產(chǎn)生影響,徑向電場對粒子軌道的改變會進一步影響等離子體輸運[22,23].在托卡馬克裝置中,可以采用在等離子體邊界區(qū)域安裝偏壓電極改善粒子和能量約束.1989年,最早在CCT托卡馬克裝置中實現(xiàn)了偏壓實驗,隨后在J-TEXT、ISTTOK托卡馬克等裝置上都進行了偏壓實驗,結(jié)果表明偏壓電極可以有效地提高徑向電場和等離子體約束性能[24-25].在類磁鏡裝置中徑向電場會影響粒子平行磁力線方向的運動速度,使得帶電粒子加速或減速,從而達到改善高能粒子約束性能的目的.因此,在CFR等線性裝置中可以通過末端加偏壓電極形成電場來改善等離子體的約束性能,這在KMAX等線性裝置中已得到有效驗證[26].

目前已有大量的數(shù)值研究工作討論了粒子約束的相關(guān)問題[27-29].與多粒子模擬、磁流體描述、統(tǒng)計描述等其他模擬方法相比,單粒子模擬能夠更加清晰地分析裝置中粒子運動的特點和軌道的拓?fù)浣Y(jié)構(gòu),為研究高能粒子、粒子輸運等提供理論基礎(chǔ)支撐[30].本文采用Boris算法[31-36]對單粒子進行模擬,研究在徑向電場作用下CFR裝置中粒子的運動軌道,并比較不同徑向電場作用對粒子約束性能的影響.第2節(jié)詳細(xì)介紹CFR磁場位形的計算和單粒子模擬方法,并研究粒子軌道特性；第3節(jié)研究徑向電場Er對CFR中高能α粒子損失率的影響；第4節(jié)總結(jié).

2 理論方法

2.1 CFR磁場位形的計算

在柱坐標(biāo)系(r,θ,z)中,對于以Z軸為中心的電流環(huán),由畢奧-薩伐爾定律可知徑向位置r和軸向位置z處產(chǎn)生的磁場分別為：

其中μ0為磁導(dǎo)率,a為電流線圈的半徑,I為電流線圈中的電流,K(k2)和E(k2)分別為第一類橢圓積分和第二類橢圓積分.

表1 CFR裝置線圈參數(shù)Table 1.Main parameters of coils in CFR.

本文研究的CFR裝置由9個半徑不同的共軸電流線圈組成,總磁場在徑向r和軸向z的分量表示為：

其中,a(n)為第n個電流線圈的半徑,I(n)為其中的電流,d(n)為第n個電流線圈所處的z軸位置,此時

CFR裝置中各線圈所處的軸向位置、線圈中半徑、線圈中所通電流的大小和方向如表1所列[14].

圖2 CFR裝置中R-Z平面上的磁場位形Fig.2.The magnetic flux of the R-Z plane in the CFR machine.

根據(jù)(3)式和(4)式以及表1中的數(shù)據(jù)編程計算出CFR裝置中磁場分布如圖2所示.從圖中可以看出,在裝置的邊緣區(qū)域磁場較強,而芯部磁場較弱,整個裝置的內(nèi)部形成了較好的磁鏡位形,在該裝置的內(nèi)部可以約束β值較高的高溫等離子體.

2.2 單粒子軌道演化方程

CFR裝置中,可以通過以下運動方程來演化高能帶電粒子的運動軌跡：

其中B為磁場強度,E為電場強度.在矢量方程(5)式和(6)式的數(shù)值求解過程中,針對因數(shù)值誤差的累積而造成的粒子能量衰減,從而導(dǎo)致軌道失真的問題,本文使用了目前比較流行的Boris算法.

Boris算法通過第k步相空間坐標(biāo)(xk,vk)求解第k+1步相空間坐標(biāo)(xk+1,vk+1).(5)式和(6)式可寫為如下離散格式：

其中h為時間步長,tk=kh,Boris算法將電場力和磁場力分開計算,通過(9)式和(10)式可從vk解析計算得到

從第k步遞推到第k+1步,電場力驅(qū)動的一半首先作用在vk上獲得vˉ,再通過(10)式計算得到v+,進一步將電場力驅(qū)動的一半加在v+上獲得vk+1.最終通過(11)式獲得vk+1.矩陣形式的Boris 算法可表示為

其中I是單位矩陣,Ωk的表達式為

圖3 Boris 算法和RK4算法模擬下第一個漂移周期和最后一個漂移周期的軌道(模擬時間1200 μs)(a)Boris 算法；(b)RK4算法Fig.3.Orbits of the first and the last drift motion periods by the Boris algorithm and the RK4 algorithm (simulation time is 1200 μs)：(a)The Boris algorithm；(b)the RK4 algorithm.

在CFR磁場位形中,由Boris算法計算得到的能量誤差和軌道失真度較小.在CFR裝置中,為驗證方法的優(yōu)越性,假定粒子位于只有z方向磁場的軸向z0=0 和徑向r0/a=0.15 處的位置上.給定粒子的初始速度與z軸垂直,φ0=0.5π.在驗證時間步長收斂性的前提下,取時間步長h為拉莫爾回旋周期的1/30.運行時間為1200 μs (約200個漂移周期),并設(shè)電場強度為2 kV/m,由原點沿半徑向外.在該磁場平面內(nèi),帶電粒子受到磁場力的作用進行拉莫爾回旋運動.由于電場、磁場梯度和磁場曲率的存在,帶電粒子會發(fā)生E×B漂移、磁場梯度漂移和磁場曲率漂移,漂移方向均沿著極向方向.圖3(a)和圖3(b)分別為Boris算法和四階龍格-庫塔RK4算法模擬得到的導(dǎo)心軌跡半徑約為0.15的粒子在第一個漂移周期和最后一個周期的軌道圖,圓環(huán)厚度近似為粒子拉莫爾回旋半徑長度.模擬結(jié)果表明,長時間迭代計算后Boris算法仍能保持軌道形狀不發(fā)生改變,與物理實際一致.RK4算法下的軌道形狀發(fā)生了改變,且完成最后一個周期需要更多的計算步數(shù),數(shù)值誤差較大.圖4中分別計算了Boris算法和四階龍格-庫塔RK4算法的相對能量誤差隨時間的變化情況.研究結(jié)果表明,Boris算法在無電場情況下能夠保證粒子的能量絕對守恒,而RK4算法的數(shù)值誤差隨模擬時間的積累不斷增大.

圖4 Boris算法和RK4算法數(shù)相對能量誤差隨時間變化對比Fig.4.Comparison of the relative energy errors between Boris method and RK4 method.

3 數(shù)值模擬結(jié)果及分析

3.1 徑向電場對α粒子軌道的影響

在徑向電場作用下,通過數(shù)值求解(5)式和(6)式模擬獲得了α粒子的運動軌跡.模擬采用T4B參數(shù)[15],裝置直徑為1 m,長2 m.功率為1 MW,在注入的中性束中,He粒子能量為25 keV.因此,本文模擬將粒子選定為能量為25 keV的α粒子.在模擬α粒子軌道的過程中,將粒子置于軸向位置z0=0 和徑向位置r0/a=0.4 處,此處磁場B0=5T,設(shè)定粒子的初始速度與Z軸的夾角φ0=0.5π,粒子的極向初始速度v?0=0.2vt,其中vt為粒子的熱速度,粒子的平行速度和垂直速度之比v‖/v⊥=1 .在CFR裝置中,將α粒子能量設(shè)定為25 keV,通過設(shè)置不同的徑向電場Er(—2 kV/m,0,+2 kV/m)研究α粒子軌道情況.如圖5所示,三種不同顏色的軌跡分別表示在無電場作用、負(fù)向徑向電場(方向指向軸心)和正向徑向電場作用下α粒子運動500 μs的軌道.研究表明,徑向電場Er會對粒子漂移速度產(chǎn)生影響,在正向徑向電場作用下,粒子漂移速度增大,α粒子反彈點向中心線圈位置移動；而在負(fù)向徑向電場作用下,粒子漂移速度減小,α粒子反彈點向兩端磁鏡線圈位置移動.

圖5 (a)CFR粒子軌道；(b)粒子軌道在X-Y截面的投影；(c)粒子軌道在R-Z截面的投影Fig.5.(a)The orbit of particles in CFR；(b)the projection of the particle orbital in the X-Y plane；(c)the projection of the particle orbital in the R-Z plane.

根據(jù)能量守恒,平行于磁場的速度v‖在徑向電場作用后可表示為

其中W為粒子總能量,m為粒子質(zhì)量,μM為磁矩,qα為α粒子電荷,Φ為電場電勢.(14)式中總能量W、電荷qα、質(zhì)量m和磁矩μM均為常量.從(14)式可知,當(dāng)粒子沿磁力線運動時,徑向電場會影響粒子平行磁力線方向的速度v‖,電勢越大,v‖越小,粒子沿Z向運動的距離越短,更不容易沿軸向尖端損失.

粒子在CFR中的漂移是由磁場的梯度漂移、磁場曲率漂移和電漂移共同產(chǎn)生的,粒子的速度可以表示為[37]

其中,ωc為回旋頻率,Rc為曲率矢量半徑,v‖,v⊥分別代表平行和垂直于磁場方向的速度,(15)式右邊第一項表示平行磁場的分量,第二項是磁場不均勻引起的總漂移,第三項是電漂移.

在正向電場作用下,電漂移引起的漂移速度和磁場不均勻引起的總漂移速度方向一致,導(dǎo)致粒子極向速度增大,平行于磁場方向的速度減小,有效地降低了粒子軸向尖端損失.而在負(fù)向電場作用下,電漂移引起的漂移速度和磁場不均勻引起的總漂移方向相反,此時徑向電場引起的電漂移和磁場梯度及曲率引起的總漂移存在競爭關(guān)系.在負(fù)向電場情況下,當(dāng)由磁場不均勻引起的漂移起主導(dǎo)作用時,粒子總漂移速度減小,平行于磁場方向的速度增加；但當(dāng)負(fù)向徑向電場增大到一定程度時,電漂移將超過磁場漂移逐漸起主導(dǎo)作用,粒子總漂移速度將增加.

3.2 絕熱區(qū)粒子模擬

在CFR裝置中心區(qū)域,由封裝的線會切以及兩個紡錘形的會切形成了一個磁勢阱結(jié)構(gòu),芯部區(qū)域的磁場趨于零,該結(jié)構(gòu)產(chǎn)生一個臨界磁面,將整個磁場位形分為絕熱區(qū)和非絕熱區(qū).在CFR裝置芯部附近磁場較弱,磁場梯度和曲率較大,類似于會切位形,該區(qū)域內(nèi)磁矩不守恒,可稱為非絕熱區(qū)；而其他區(qū)域磁矩守恒,稱為絕熱區(qū)[16].根據(jù)緊湊型先進磁鏡系統(tǒng)的磁場位形分布特點容易看出,絕熱區(qū)粒子更容易沿軸向尖端損失,當(dāng)高能粒子的初始速度與磁場平行時,粒子約束性最差,粒子極容易通過磁鏡的端口而損失.選取在絕熱區(qū)最容易損失的能量為25 keV的α粒子進行模擬.模擬過程中設(shè)定時間步長Δt= 1×10—12s,z0=0,φ0=0,粒子的極向初始速度v?0=0,在r0/a=0.2—0.5 范圍內(nèi),徑向每隔0.05的長度取一個徑向點.當(dāng)v‖/v=1,即初始速度與磁場方向平行時,粒子極容易沿著磁力線從端口逃離裝置.在當(dāng)前粒子初始條件下,α粒子僅僅運動了2.354 μs便全部直接損失(如圖6所示).

圖6 CFR中α粒子損失軌道Fig.6.The lost orbits of theαparticles in CFR.

為了研究在絕熱區(qū)不同徑向位置處徑向電場對能量為25 keV的α粒子最易損失方向(平行于Z軸方向)約束性能的影響,將粒子初始條件設(shè)置為z0=0,φ0=0,v?0=0 和v‖/v=1,在r0/a=0.2— 0.5范圍內(nèi),徑向每隔0.05的長度取一個徑向點,在同一個徑向位置,通過不斷調(diào)整正向徑向電場的大小研究粒子約束性能的變化.研究表明,在沒有徑向電場存在的情況下,能量為25 keV的α粒子只能在裝置中停留幾個微秒(如圖6).隨著正向徑向電場的增強,粒子在裝置內(nèi)的約束時間不斷增長,當(dāng)正向電場增大到某一臨界值時,在最容易損失位置處注入的α粒子軌道將會長時間停留在裝置內(nèi)部而不是逃離裝置,粒子能夠很好地約束在磁鏡內(nèi)Z處于—1.15—1.15之間(如圖7).在模擬中,選定不同的徑向初始位置,分析粒子在不同徑向電場作用下的運動軌跡得到了相同的結(jié)果.圖8中給出了在r0/a=0.3 處,不同徑向電場作用下α粒子的運動軌跡,由圖可見,當(dāng)正向電場由58620 V/m增大到58621 V/m時,粒子能夠被長時間地約束在CFR裝置內(nèi),說明了正向徑向電場的存在可以改善粒子約束性能.數(shù)值計算的結(jié)果表明,不同徑向位置處使得粒子能夠長時間約束的徑向電場臨界值與粒子的初始條件有關(guān)：粒子初始位置越靠近軸心位置(r= 0),正向的徑向電場臨界值越大.當(dāng)徑向電場大于臨界場值,在絕熱區(qū)運動的α粒子都能被很好地約束在CFR裝置內(nèi),區(qū)域內(nèi)粒子損失率將降至最低,甚至能將所有粒子約束.

圖7 粒子約束時間隨著正向電場強度的變化Fig.7.The particle confinement time changes with the intensity of the positive electric field.

圖8 初始位置 r0/a=0.3 處α粒子在不同正向徑向電場作用下的運動軌跡Fig.8.The orbits of α particle under the different positive radial electric fields at r0/a=0.3 .

同時研究了負(fù)向徑向電場對能量為25 keV的α粒子約束性能的影響,研究過程中選取粒子的初始狀態(tài)與正向電場情況下相同.研究表明,負(fù)向徑向電場的存在同樣可以改善粒子的約束性能.然而,負(fù)向徑向電場對粒子約束性能的影響隨著電場強度的增加,在所有位置的約束時間并未表現(xiàn)出都是單調(diào)遞增的趨勢：在負(fù)向電場強度比較小時,粒子的約束時間隨電場強度的增大而減小,表明此時電場使粒子更容易損失；而當(dāng)電場強度增大到一定程度時,在最容易損失位置注入的α粒子卻能夠很好地約束在CFR裝置內(nèi)(如圖9和圖10).

以上分析了在極限情況下,最容易損失的粒子通過外加徑向電場得到了很好的約束.為了更直觀地探究徑向電場對絕熱區(qū)無碰撞α粒子損失率的影響,選取初始條件為z0=0,φ0=0,v?0=0.2vt,使粒子連續(xù)分布在徑向初始位置r0/a= (0.2—0.5)之間,以不同的入射角(0-2π)發(fā)射500個α粒子,得到粒子運行500 μs后的損失率(如圖11所示).研究結(jié)果表明,對于絕熱區(qū)運動的高能粒子,隨著正向徑向電場增大,粒子損失率降低.而在負(fù)向徑向電場作用下,隨著電場強度增大,粒子損失率先增大后減小.盡管將高能α粒子全部約束在CFR裝置中需要提供極大的徑向電場,但在電極偏壓裝置可控的徑向電場下已經(jīng)能夠有效地改善粒子損失,提高粒子約束性.

圖9 粒子約束時間隨著負(fù)向徑向電場強度的變化Fig.9.The variations of particle confinement time with the intensity of the negative radial electric field.

圖10 初始位置 處α粒子在不同負(fù)向徑向電場作用下的軌道r0/a=0.3Fig.10.The orbits of α particle under the different negative radial electric fields at .r0/a=0.3

圖11 不同徑向電場作用下的α粒子損失率Fig.11.The loss rates of the α particles with different radial electric fields.

4 結(jié) 論

徑向電場的存在會改變裝置內(nèi)帶電粒子平行磁力線方向的速度,對裝置內(nèi)高能粒子的約束產(chǎn)生影響.本文從牛頓運動方程出發(fā),運用Boris算法模擬α粒子在CFR裝置中的運動軌道,探究了不同徑向電場作用下α粒子在CFR絕熱區(qū)約束性能的差異.模擬發(fā)現(xiàn)：在正向徑向電場作用下,α粒子運動電漂移速度與磁場梯度和曲率引起的漂移存在競爭,α粒子運動軌道反彈點沿磁力線向徑向電場方向移動.隨著電場強度增大,α粒子在CFR中約束性變好,當(dāng)徑向電場強度達到一定臨界值時,即使最容易損失的α粒子也能夠很好地約束在裝置內(nèi)；在負(fù)向徑向電場作用下,α粒子漂移速度增大.隨著電場強度的不斷增大,初始速度平行于Z軸方向的極容易損失的α粒子在CFR中約束時間先減少,但當(dāng)電場強度達到某一特定值時,粒子約束時間不斷增加,從而很好地約束在CFR裝置內(nèi)部.在CFR裝置中,當(dāng)所施加的電場強度增大到一定臨界值時,粒子能夠長時間地約束在裝置內(nèi)部,約束時間可以達到與當(dāng)前托卡馬克裝置粒子約束時間相當(dāng)?shù)暮撩肓考?由于在CFR裝置內(nèi)的磁場使用的是最優(yōu)化磁場位形,裝置內(nèi)部的高能粒子可以達到毫秒量級約束時間,然而該類線性裝置較托卡馬克裝置具有約束效率(β值)更高、尺寸更小和造價更低等諸多優(yōu)點,在未來的磁約束聚變領(lǐng)域具有較好的應(yīng)用前景.