波紋場和新經(jīng)典撕裂模對托卡馬克中快離子輸運(yùn)影響的數(shù)值模擬研究

劉 錦 陳嘉盈 常悅悅 詹 研 虞立敏

（華東理工大學(xué)理學(xué)院上海200237）

化石能源的枯竭是人類必須面對的問題，可控核聚變?yōu)槿祟惤鉀Q能源危機(jī)提供了一個(gè)可能的思路。托卡馬克是當(dāng)前熱門的可控核聚變裝置［1-2］。受控磁約束核聚變領(lǐng)域一個(gè)重要的問題是如何更好地約束快離子，通常擾動(dòng)場（包括裝置本身固有的波紋場和不穩(wěn)定性產(chǎn)生的擾動(dòng)場）引起快離子再分布和損失，而快離子損失不僅會(huì)降低等離子體溫度，降低聚變反應(yīng)率，而且也會(huì)損傷裝置器壁，影響裝置的壽命。因此研究快離子的輸運(yùn)和損失的物理機(jī)制對提升核聚變裝置中等離子體性能和延長裝置壽命有重要的意義［3-5］。

波紋擾動(dòng)和新經(jīng)典撕裂模（Neoclassical Tearing Mode，NTM）對快粒子輸運(yùn)乃至對裝置運(yùn)行造成重要影響，因而得到廣泛的研究。波紋擾動(dòng)指由于環(huán)向場線圈有限，有線圈的地方磁場大、無線圈的地方磁場小導(dǎo)致的類似波紋的環(huán)向場［3，6-7］。Wu等［8］使用ORBIT和NUBEAM/TRANSP代碼計(jì)算了不同中性束注入方式下的束離子瞬態(tài)損失和波紋損失，提出波紋損失的主要過程是捕獲粒子拐點(diǎn)徑向漂移引起的波紋隨機(jī)擴(kuò)散，以及提高等離子體電流可以有效減少離子損失；李吉波［3］利用ORBIT模擬了不同的等離子體電流、快離子螺距角分布等因素對波紋損失的影響，得出等離子體電流較小和快離子螺距角較大的時(shí)候，損失份額較大；Nave等［9］研究了離子回旋射頻（Ion-Cyclotron Radio Frequency，ICRF）加熱下的等離子體中波紋場振幅由0.08%～1.5%的情況，研究顯示：由于熱離子和快離子的輸運(yùn)產(chǎn)生的扭矩會(huì)引起等離子體反向旋轉(zhuǎn)。

NTM是在經(jīng)典撕裂模（TM）上發(fā)展而來。TM考慮了q=m/n磁面以內(nèi)的區(qū)域有限電阻效應(yīng)和慣性效應(yīng)，NTM額外考慮了自舉電流［10］。Strumberger等［11］計(jì)算了中性束注入下由NTM和波紋場引起的快離子損失，發(fā)現(xiàn)NTM造成的通行粒子損失集中在tloss=40～200μs，造 成 的 捕 獲 粒 子 損 失 集 中 在tloss≥200μs，波紋場造成的捕獲粒子損失集中在tloss≥35μs，而波紋場幾乎不影響通行粒子；Gobbin等［12］利用ORBIT模擬還原了ASDEX Upgrade中的NTM引起的磁島效應(yīng)，提出NTM不僅增加了損失粒子總量，并且使損失粒子通量集中在一個(gè)更小的區(qū)域；Zhang等［7］在對NTM和波紋場的研究中得出粒子在NTM和波紋的共同作用下會(huì)產(chǎn)生更大的軌道漂移，尤其是當(dāng)粒子更多或更長時(shí)間地通過磁島中心時(shí)；Bardóczi等［13］利 用TRANSP-‘Kick’輸 運(yùn) 模 型 對NTM進(jìn)行的模擬研究，發(fā)現(xiàn)相空間中重疊共振會(huì)導(dǎo)致快離子在重疊島鏈間的混沌運(yùn)動(dòng)，重疊共振在磁島邊緣產(chǎn)生且造成快離子損失。

本文主要使用數(shù)值模擬的方法，依靠導(dǎo)心運(yùn)動(dòng)代碼ORBIT［14］，研究波紋擾動(dòng)和新經(jīng)典撕裂模（NTM）對快粒子輸運(yùn)的影響。我們計(jì)算帶1個(gè)正電荷，質(zhì)量為2個(gè)質(zhì)子質(zhì)量的氘離子，采用NTM的模數(shù)為m/n=2/1（其中：m為極向模數(shù)；n為環(huán)向模數(shù)）和全超導(dǎo)托卡馬克核聚變實(shí)驗(yàn)裝置（Experimental and Advanced Superconducting Tokamak，EAST）中的波紋場模型，對不同參數(shù)情形下粒子損失和環(huán)向角動(dòng)量擴(kuò)散展開了研究。

1 計(jì)算方法

1.1 ORBIT程序理論模型

根據(jù)漂移動(dòng)力學(xué)哈密頓理論，可以得到ORBIT代碼采用的微分方程模型為［14-15］：

式中：H為帶電粒子的動(dòng)能和靜電勢能之和；Φ為靜電勢能；ρ//為歸一化的回旋半徑；μ為磁矩；g正比于ψp外的等離子體極向電流；I正比于ψp內(nèi)環(huán)向總電流。在磁面坐標(biāo)系下，θ是極向角，ζ是環(huán)向角，ψp為極向磁通，ψ為環(huán)向磁通。平衡磁場表示為

1.2 NTM擾動(dòng)模型

ORBIT中可以處理的擾動(dòng)形式為：

式中：ω是擾動(dòng)的角頻率，根據(jù)低頻NTM的情形［16-17］，本文分別對頻率f為0 kHz、3 kHz、6 kHz、9 kHz的情況進(jìn)行了討論；m和n分別代表了擾動(dòng)在極向和環(huán)向的模數(shù)，下面計(jì)算中我們?nèi)/n=2/1。

模的徑向結(jié)構(gòu)使用解析表達(dá)式［18］：

式中：ψw是等離子體邊界的極向磁通；ψp是r處的磁通；α0是擾動(dòng)的幅度。擾動(dòng)的形狀類似鐘形曲線。擾動(dòng)的最大值等于α0，擾動(dòng)磁場δBr最大值出現(xiàn)位置受安全因子q(r)的控制，出現(xiàn)在q(r)=m/n=2處［10］。

擾動(dòng)幅度α0的選取參考了文獻(xiàn)［12］中根據(jù)實(shí)驗(yàn)磁島寬度選取的方法。EAST上測量的磁島寬度為Δr/a=0.15［19］，本文中a=40 cm，則磁島寬度應(yīng)為6 cm。因此取α0=5×10-5，使得磁島寬度與實(shí)驗(yàn)值相近。

1.3 波紋場模型

波紋擾動(dòng)指由于環(huán)向場線圈有限，有線圈的地方磁場大，無線圈的地方磁場小導(dǎo)致的環(huán)線磁場不均勻性。波紋場可以表述為：

典型圓截面中，?。?/p>

式中：R0是波紋中心的大半徑；br是托卡馬克截面的橢圓率，對于本文的圓截面托卡馬克，br=π/4。李吉波［3］對EAST中快粒子波紋損失做了詳細(xì)的研究。為了貼合實(shí)際情況，式（5）的其余參數(shù)均選取文獻(xiàn)［3］中EAST波紋同樣的值：

1.4 平衡磁面

為了簡化計(jì)算，本文采用解析數(shù)值平衡下圓截面托卡馬克，大徑R0=185cm，截面半徑a=40 cm。在這種情況下，將平衡磁場分解為環(huán)向BT和極向Bp，則可以近似表示為：

本文取BT0=1.8 T。式（7）中q=dψ/dψp，稱為安全因子。平衡磁面如圖1所示。

圖1 平衡磁面Fig.1 Balanced magnetic surface

而安全因子采用三次曲線擬合的表達(dá)式，如下：

滿足：

取α0=5×10-5時(shí)，安全因子q和δBr/B隨歸一化磁面的關(guān)系如圖2。

圖2 安全因子、模幅度和擾動(dòng)磁場大小(a)安全因子和模幅度，(b)安全因子和擾動(dòng)磁場Fig.2 Safety factor,mode amplitude,and magnitude of the disturbance magnetic field(a)Safety factor and mode amplitude,(b)Safety factor and disturbance magnetic field

1.5 其他條件

本文中計(jì)算采用的粒子是帶一個(gè)正電荷、質(zhì)量為兩個(gè)質(zhì)子質(zhì)量的氘離子，不考慮粒子的碰撞和散射。

2 計(jì)算結(jié)果與討論

2.1 擾動(dòng)對粒子損失和擴(kuò)散的影響

為了探究擾動(dòng)對粒子損失數(shù)的影響，我們分為以下4類情況進(jìn)行討論，分別是無徑向磁擾動(dòng)、添加波紋擾動(dòng)、添加NTM擾動(dòng)以及同時(shí)添加波紋和NTM擾動(dòng)。令α0=5×10-5，粒子能量為60 keV。初始粒子均勻分布在同一個(gè)磁面上，而投射角λ（v///v）相同。為了使計(jì)算結(jié)果收斂并節(jié)省計(jì)算時(shí)間，我們采用的粒子數(shù)目為3 000個(gè)。掃描的參數(shù)范圍為：歸一化的磁面從0.1變化到0.9、步長為0.1，投射角從0變化到0.9、步長為0.1。采用插值的方法繪制熱力圖來統(tǒng)計(jì)粒子的損失和擴(kuò)散情況。其中，損失是指一旦粒子到達(dá)最外層磁面就認(rèn)為發(fā)生損失，沒有考慮粒子再回到等離子體的情形；而擴(kuò)散則通過歸一化的環(huán)向角動(dòng)量（即Pζ/ψw，為了簡潔記為Pζ）的方差來表示，設(shè)第i個(gè)粒子t時(shí)刻的環(huán)向角動(dòng)量為Pζi(t)，所有粒子t時(shí)刻的環(huán)向角動(dòng)量的平均值為-Pζ(t)，則有如式（10）［20］：

由于方差每個(gè)時(shí)刻都在變化，而最值得關(guān)注的是方差的最大值，因?yàn)樗砹肆Ｗ訑U(kuò)散程度最大的情況，因此取其所有時(shí)刻的最大值。

圖3 中，橫坐標(biāo)為初始粒子的歸一化磁面，縱坐標(biāo)是初始粒子的投射角。顏色棒表明初始粒子在該磁面和投射角下最終損失粒子數(shù)，數(shù)值模擬使用了3 000個(gè)粒子，因此顏色棒最大值為3 000。圖3（a）是沒有加擾動(dòng)情況下粒子損失情況，屬于第一軌道損失。為了清楚擾動(dòng)造成的粒子損失，在圖3（b）、（c）、（d）中去除了第一軌道損失粒子。從圖3（a）中可以得出，粒子靠外磁面、低投射角分布的情況下?lián)p失粒子數(shù)多，這和預(yù)期是一致的。低投射角意味著粒子平行速度分量小，有大的概率成為捕獲粒子。磁面高意味著粒子初始位置接近托卡馬克內(nèi)壁，容易造成第一軌道損失。與沒有添加任何擾動(dòng)相比，添加波紋增加了低投射角、內(nèi)磁面的粒子的損失并有跳變增加現(xiàn)象，添加NTM顯著增加了磁面0.4～0.9下投射角大于0.4的粒子損失，磁面為0.4處劇烈的粒子損失是由于該磁面處于δBr極大值附近。而NTM和波紋都添加的損失情況近似于它們單獨(dú)添加時(shí)造成損失之和。

環(huán)向角動(dòng)量的擴(kuò)散同樣是我們關(guān)注的一個(gè)物理量，可以間接反映離子在擾動(dòng)作用下的擴(kuò)散情況。我們用粒子群運(yùn)行所有時(shí)刻中全部粒子環(huán)向角動(dòng)量的方差的最大值來描述。統(tǒng)計(jì)結(jié)果如圖4所示。

圖3 粒子損失統(tǒng)計(jì)圖(a)沒有徑向磁擾動(dòng)，(b)添加NTM，(c)添加波紋，(d)添加波紋和NTM（(b)、(c)、(d)圖中已去除第一軌道損失粒子）Fig.3 Particle loss statistics(a)Without radial magnetic perturbation,(b)With NTM,(c)With ripple field,(d)With both of ripplefield and NTM(First orbit losses have been removed in figure(b),(c),(d))

圖4 ?統(tǒng)計(jì)圖(a)沒有徑向磁擾動(dòng)，(b)添加NTM，(c)添加波紋，(d)添加波紋和NTM Fig.4?statistics(a)Without radial magnetic perturbation,(b)With NTM,(c)With ripple field,(d)With both of ripple field and NTM

由上述討論可知，波紋和NTM擾動(dòng)并不是對所有粒子都有影響，而是分別對捕獲粒子和通行粒子造成影響。我們分別計(jì)算了捕獲粒子和通行粒子在這兩種不同條件下的軌道。設(shè)初始粒子能量40 keV，捕獲粒子初始磁面0.7，投射角0.3；通行粒子初始磁面0.7，投射角0.8。

圖5 中顯示了不同擾動(dòng)條件下相同初始條件的粒子軌跡。在波紋影響下，捕獲粒子軌道出現(xiàn)劇烈偏移，但是通行粒子軌道并沒有發(fā)生明顯變化。在NTM的影響下捕獲粒子軌道同樣出現(xiàn)偏移，但是程度不如波紋影響的劇烈。通行粒子的軌道線條同樣變寬了，體現(xiàn)了NTM對通行粒子造成影響。當(dāng)兩者都疊加時(shí)，捕獲粒子的軌道漂移更加劇烈，說明波紋和NTM的協(xié)同作用將造成更大的環(huán)向角動(dòng)量擴(kuò)散。單粒子運(yùn)行的結(jié)果印證了波紋擾動(dòng)對捕獲粒子作用強(qiáng)烈，是捕獲粒子環(huán)向角動(dòng)量擴(kuò)散增加的主要原因。這解釋了圖4中波紋場影響下環(huán)向角動(dòng)量擴(kuò)散在λ＜0.5的區(qū)域特別強(qiáng)烈的現(xiàn)象。

2.2 擾動(dòng)對不同相空間粒子的影響

為研究擾動(dòng)對不同種類的粒子影響，我們分析了相空間中均勻分布的粒子的變化。高能粒子可以按其損失與否分為損失（Lost）和磁約束（Confined）。按照運(yùn)動(dòng)狀態(tài)分為捕獲（Trapped）、正向通行（Copassing）、逆 向 通 行（Counter-passing）和 停 滯（Stagnation）。用軌道分類圖來表述這6種粒子的分布，橫坐標(biāo)為歸一化環(huán)向角動(dòng)量Pζ/ψw，縱坐標(biāo)為歸一化磁矩，如圖6所示。

我們設(shè)置粒子能量為40 keV，選取初始粒子在相空間中均勻分布，調(diào)整NTM和波紋擾動(dòng)的有無進(jìn)行計(jì)算。從計(jì)算結(jié)果中，統(tǒng)計(jì)損失粒子和未損失粒子的初始相空間中的分布，得到如下的結(jié)果。

圖5 通行粒子和捕獲粒子軌道圖(a)沒有徑向磁擾動(dòng)，(b)添加NTM，(c)添加波紋，(d)添加波紋和NTMFig.5 Trajectories of passing particles and trapped particles(a)Without radial magnetic perturbation,(b)With NTM,(c)With ripple field,(d)With both

圖6 不同類型的粒子在相空間分布Fig.6 Distribution of different types of particles in phase space

圖7 中顯示了某損失粒子和某未損失的粒子的初始相位。對比發(fā)現(xiàn)，在不加擾動(dòng)的情況下，損失和未損失粒子按相位圖的邊界分布。加了NTM的情況下，逆向通行損失粒子和正向通行損失粒子開始入侵未損失粒子的邊界。我們感興趣的是隨著擾動(dòng)幅度α0增加，相位圖會(huì)發(fā)生怎樣的改變。于是在其他條件都相同的前提下，令粒子能量為40 keV，以便觀察到寬的Trapped Confined區(qū)域。α0的變化區(qū)間從3×10-5變化到9×10-5，不添加波紋，觀察粒子相位圖的變化。如圖8所示，隨著擾動(dòng)幅度的增加，位于Co-passing Confined和Counter-passing Confined區(qū)域的粒子逐漸損失，Trapped Confined區(qū)域也有同樣的轉(zhuǎn)變，但是比例明顯小于其他區(qū)域。

同時(shí)，我們也發(fā)現(xiàn)波紋場的主要影響造成Trapped Confined區(qū)域內(nèi)粒子大量損失。隨著粒子能量的升高，Trapped Confined區(qū)域會(huì)越來越小，不難預(yù)見波紋場造成的影響會(huì)越來越小。為了證實(shí)這一觀點(diǎn)，我們在運(yùn)行中只添加了波紋，并改變粒子能量從40～120 keV，觀察相位圖變化。我們得到的結(jié)果如圖9所示。

從圖9（a）～（c）三幅圖的變化規(guī)律中可以看到，隨著粒子能量增加，波紋影響的區(qū)域越來越小，最后消失不見，說明粒子能量越高，波紋擾動(dòng)的影響就越小。

2.3 加入NTM頻率對損失和擴(kuò)散的影響

上面的研究都是基于靜態(tài)NTM擾動(dòng)，但實(shí)際上大多數(shù)粒子都是在毫秒量級內(nèi)損失，這和NTM的周期是同一個(gè)數(shù)量級［12］。因此需要考慮NTM的頻率帶來的影響。假設(shè)NTM以恒定頻率旋轉(zhuǎn)，我們設(shè)置頻率為0 kHz、3 kHz、6 kHz、9 kHz，擾動(dòng)幅度為α0=5×10-5，粒子能量為60 keV，不加入波紋場。同樣做出粒子的損失和環(huán)向角動(dòng)量擴(kuò)散圖。

圖7 損失粒子和未損失粒子的初始相位(a)沒有徑向磁擾動(dòng)，(b)添加NTM，(c)添加波紋，(d)添加波紋和NTMFig.7 The initial phase of the lost and confined particles with(a)Without radial magnetic perturbation,(b)With NTM,(c)With ripple field,(d)With both.

圖8 改變擾動(dòng)幅度對相位圖的影響(a)α0=3×10-5，(b)α0=6×10-5，(c)α0=9×10-5Fig.8 The effect of changing the amplitude of the disturbance on the phase map(a)α0=3×10-5,(b)α0=6×10-5,(c)α0=9×10-5

圖9 只添加波紋，改變粒子能量，觀察相位圖變化(a)E=40 keV，(b)E=80 keV，(c)E=120 keVFig.9 The change of phase diagram caused by various particle energies when only ripple field is applied(a)E=40 keV,(b)E=80 keV,(c)E=120 keV

圖10 NTM不同頻率的粒子損失統(tǒng)計(jì)圖(a)f=0 kHz，(b)f=3kHz，(c)f=6 kHz，(d)f=9 kHzFig.10 Statistics of particle loss at different frequencies of NTM(a)f=0 kHz,(b)f=3 kHz,(c)f=6kHz,(d)f=9 kHz

圖11 NTM不同頻率的統(tǒng)計(jì)圖(a)f=0 kHz，(b)f=3 kHz，(c)f=6 kHz，(d)f=9 kHzFig.11 Statistics of at different frequencies of NTM(a)f=0 kHz,(b)f=3 kHz,(c)f=6 kHz,(d)f=9 kHz

從圖10、圖11中我們發(fā)現(xiàn)頻率加入后，粒子損失情況保持了原有的特征，對初始磁面0.7～0.8、投射角大于0.4的范圍，可以看到隨著頻率的提高，損失數(shù)量逐步減少。而對于環(huán)向角動(dòng)量擴(kuò)散情況，不變之處在于磁面從0.7～0.9、投射角0～0.2的范圍有一個(gè)明顯的高峰，這可能與該處粒子非常容易損失有關(guān)。對投射角大于0.4的區(qū)域，隨著頻率加入發(fā)生了明顯的變化：原先的條帶狀分布被破壞，磁面為0.4和0.6的亮紋產(chǎn)生的混合、0.8處的亮紋逐漸減弱；淺藍(lán)色區(qū)域增加，尤其是f=9 kHz的情形。這些跡象都表明NTM頻率增加使得環(huán)向角動(dòng)量擴(kuò)散變得更加均勻。

3 結(jié)語

本文使用導(dǎo)心運(yùn)動(dòng)代碼ORBIT，研究了波紋場和新經(jīng)典撕裂模對圓截面托卡馬克粒子輸運(yùn)的影響。計(jì)算發(fā)現(xiàn)，擾動(dòng)能增加離軸、低投射角的快粒子的損失率。具體地，與沒有添加任何擾動(dòng)相比，添加波紋增加了低投射角粒子的損失，添加NTM顯著增加了磁面0.4～0.9下投射角大于0.4的粒子損失，而兩個(gè)都添加的損失情況近似于單獨(dú)添加時(shí)造成損失之和。添加波紋會(huì)顯著增加環(huán)向角動(dòng)量擴(kuò)散，而在NTM擾動(dòng)最強(qiáng)的區(qū)域環(huán)向角動(dòng)量擴(kuò)散也相對顯著，波紋擾動(dòng)和NTM的對環(huán)向角動(dòng)量擴(kuò)散的協(xié)同作用遠(yuǎn)大于單獨(dú)作用之和。波紋擾動(dòng)對λ＜0.5的區(qū)域影響?yīng)q為劇烈，這部分粒子大多屬于捕獲粒子。從單粒子軌道圖上，在波紋場影響下捕獲粒子軌道發(fā)生偏移，其程度遠(yuǎn)大于NTM單獨(dú)作用導(dǎo)致的偏移。證明了波紋擾動(dòng)是造成粒子環(huán)向角動(dòng)量擴(kuò)散的主要原因。

在相位圖上，位于同向通行約束和反向通行約束區(qū)域的粒子會(huì)受到NTM影響而損失，波紋場會(huì)造成磁捕獲約束區(qū)域內(nèi)粒子大量損失。隨著粒子能量升高，磁捕獲約束區(qū)域變小，波紋擾動(dòng)的影響就越小。

當(dāng)考慮NTM頻率時(shí)，粒子損失情況基本保持原始特征，但是對初始磁面0.7～0.8、投射角大于0.4的范圍，隨著頻率的提高，損失粒子數(shù)量逐步減少。同時(shí)我們還發(fā)現(xiàn)，NTM頻率增加使得環(huán)向角動(dòng)量擴(kuò)散變得更加均勻。

我們的工作研究了波紋場和2/1 NTM對不同初始狀態(tài)粒子的作用情況，對托卡馬克裝置的設(shè)計(jì)、高能粒子的約束都有參考價(jià)值。后續(xù)還將從托卡馬克截面的設(shè)計(jì)、安全因子曲線的選取等方面做進(jìn)一步研究。