鐘翊君，龔學(xué)余,*，黃千紅，譚清懿，龔 蕾，杜 丹

(1.南華大學(xué)，湖南 衡陽 421001；2.福清核電站，福建 福清 350300)

高溫等離子體的約束和輔助加熱是實(shí)現(xiàn)受控?zé)岷朔磻?yīng)的基本問題。自從ASDEX[1]首次在托卡馬克上發(fā)現(xiàn)改善約束性能的H模后，大量實(shí)驗(yàn)和理論研究表明，由負(fù)徑向電場引起的等離子體極向旋轉(zhuǎn)速度的急劇增加是誘發(fā)L-H模轉(zhuǎn)化和改善約束性能的重要原因[2-6]。在托卡馬克磁場位形下，相同能量的電子和離子，由于離子在磁場分布不均勻引起的梯度漂移和曲率漂移作用下產(chǎn)生的跨越磁面的運(yùn)動比電子的大，從而形成電子和離子相對分離，有些邊界離子還會損失出等離子體區(qū)域，在邊緣出現(xiàn)凈電子，凈電子的空間分布會在邊緣產(chǎn)生一個負(fù)的徑向電場。與此同時，由于負(fù)徑向電場的產(chǎn)生，帶電粒子運(yùn)動過程中增加了一個較大的電漂移，反過來影響電子和離子的運(yùn)動軌跡和約束性能。因此，關(guān)于負(fù)徑向電場存在時，粒子運(yùn)動情況的研究是磁約束聚變等離子體研究領(lǐng)域的重要問題。另一方面，電子俘獲率在電子回旋波電流驅(qū)動和功率沉積的研究中有著重要意義，俘獲電子會吸收更多的波功率，使波功率的沉積更加集中，并使電流驅(qū)動效率明顯降低[7-9]。同時，在中性束注入下快離子產(chǎn)生的自舉電流[10-11]及低雜波電流驅(qū)動和功率沉積等的研究中，也需先知快離子的俘獲率。然而，對于負(fù)徑向電場存在時俘獲率的計算尚未見相關(guān)報道。

本文采用有限元法和高斯迭代法[12-13]，結(jié)合實(shí)驗(yàn)測得的中心邊緣密度、溫度，等離子體總電流及中心電流密度，數(shù)值求解平衡方程[14]，得到等離子體磁場，從而模擬氘離子在負(fù)徑向電場存在時的運(yùn)動軌跡，然后采用蒙特卡羅方法抽取滿足麥克斯韋速度分布的氘離子速度樣本，統(tǒng)計氘離子的俘獲率，研究負(fù)徑向電場對氘離子俘獲率的影響。

1 物理模型

1.1 平衡位形的求解

結(jié)合EAST裝置實(shí)驗(yàn)平衡數(shù)據(jù)，采用有限元法和高斯迭代法求解Grad-Shafranov方程[14]，計算軸對稱等離子體平衡位形下的磁通分布：

(1)

p(ψ)和g(ψ)2均可表示為關(guān)于角向磁通的多項(xiàng)式形式：

(2)

(3)

其中：n為電子離子密度；ti、te分別為離子、電子溫度；nc與ne分別為電子和離子的中心密度與邊緣密度；ψc為中心角向磁通；tec、tic分別為電子和離子的中心溫度；tee、tie分別為電子和離子的邊緣處溫度；b1和b2為常數(shù)。

為了解式(1)，式(2)、(3)中的參數(shù)必須被確定。先假定一個中心角向磁通ψc，然后假定b1和b2的初值，最后用有限元法進(jìn)行網(wǎng)格剖分，用程序求出的總電流、中心電流密度和實(shí)驗(yàn)可間接測得的總電流Ip、中心電流密度j0相等以及程序求出的中心角向磁通和假定的中心角向磁通相等做迭代條件，反演出與實(shí)驗(yàn)相符合的b1、b2和ψc。

由上述方法可得到各網(wǎng)格點(diǎn)的環(huán)向電流密度和極向磁通等物理量，間接求出平衡位形下各點(diǎn)的磁場：

(4)

其中，B為磁場。

1.2 帶電粒子運(yùn)動軌跡的求解

粒子在托卡馬克位形中運(yùn)動，滿足洛倫茲力方程：

(5)

其中：m為質(zhì)量；v為速度；t為時間；q為帶電量；E為電場；B為磁場，通過式(4)得到，直角坐標(biāo)系中，磁場采用(Bx,By,Bz)表示，具體對應(yīng)關(guān)系示于圖1。

圖1 磁場分布Fig.1 Magnetic field distribution

托卡馬克中粒子運(yùn)動方程可表示為：

(6)

負(fù)徑向電場Er表示為：

其中：E0為常數(shù)；r為徑向位置。用四階龍格庫塔方法求解式(6)，計算平衡位形下各點(diǎn)的帶電粒子運(yùn)動軌跡。

2 數(shù)值模擬結(jié)果與分析

2.1 氘離子在平衡位形中的運(yùn)動軌跡

根據(jù)上述理論，首先采用EAST裝置(炮號33068)的放電數(shù)據(jù)(R=1.8 m，a=0.45 m，Ip=0.5 MA，j0=1.5 MA·m-2，中心處的環(huán)向磁場Bt=2.0 T，nc=3×1019m-3，ne=1×1018m-3，tec=3 keV，tic=2 keV，tie=tee=0.1 keV)求解式(1)，得到平衡位形。然后根據(jù)式(4)計算截面各網(wǎng)格點(diǎn)的磁場，最后采用龍格庫塔方法求解式(6)，模擬氘離子在托卡馬克電磁場中的運(yùn)動軌跡。帶電粒子在托卡馬克電磁場中的運(yùn)動是以速度vm沿螺旋磁力線的移動加上由磁場不均勻性引起的曲率漂移和梯度漂移合成速度vDB以及電場引起的電漂移速度vE所帶來的運(yùn)動。其中vDB在該平衡位形中離子的漂移方向是向上的，vE在該平衡位形中離子的漂移方向是沿極向磁場反方向的。圖2示出不同負(fù)徑向電場下氘離子的運(yùn)動軌跡在小環(huán)方向上的投影，氘離子初始位置(x0，y0，z0)=(2.227 5 m，0，0)，初始速度(vx，vy，vz)=(0，1.65×105m/s，0)。圖2a中軌跡1為通行粒子擴(kuò)張軌跡，軌跡2為外翻香蕉粒子軌跡，軌跡3為內(nèi)翻香蕉粒子軌跡，軌跡4為通行粒子擴(kuò)張軌跡。由圖2a可看出，隨負(fù)徑向電場的增大，氘離子軌跡由擴(kuò)張通行軌跡向外翻香蕉軌跡再向內(nèi)翻香蕉軌跡，最后向壓縮通行軌跡演變。圖2b中速度vmE表示離子沿磁力線方向的運(yùn)動速度vm與電場引起的電漂移速度vE的合成速度，與圖2a所示軌跡相互對應(yīng)。軌跡1、2所代表的氘離子由于洛倫茲力從初始點(diǎn)開始往赤道面以上運(yùn)動，加上向上的vDB，離子的運(yùn)動軌跡在初始磁面外出現(xiàn)擴(kuò)張軌跡。軌跡3和4所代表的粒子從初始點(diǎn)開始往赤道面以下運(yùn)動，加上向上的vDB，離子的運(yùn)動軌跡在初始磁面內(nèi)出現(xiàn)壓縮軌跡。軌跡2的合成速度小于軌跡1的，于是氘離子軌跡由通行軌跡向香蕉軌跡演變。

1——E0=0 V/m；2——E0=15 000 V/m；3——E0=30 000 V/m；4——E0=80 000 V/m圖2 初始速度(vx，vy，vz)=(0，1.65×105 m/s，0)時氘離子運(yùn)動軌跡在小環(huán)方向上的投影(a)及速度合成(b)Fig.2 Projection of poloidal cross-section of deuterium ion orbit (a) and velocity compose (b) under initial velocity (vx, vy, vz)=(0, 1.65×105 m/s, 0)

1——E0=0 V/m；2——E0=15 000 V/m；3——E0=30 000 V/m；4——E0=80 000 V/m圖3 初始速度(vx，vy，vz)=(1.34×105 m/s，6.7×104 m/s，6.7×104 m/s)時氘離子運(yùn)動軌跡在小環(huán)方向上的投影(a)及速度合成(b)Fig.3 Projection of poloidal cross-section of deuterium ion orbit (a) and velocity compose (b) under initial velocity (vx, vy, vz)=(1.34×105 m/s, 6.7×104 m/s, 6.7×104 m/s)

圖4示出初始位置(x0，y0，z0)=(2.227 5 m，0，0)、初始速度(vx，vy，vz)=(-1.34×105m/s，-6.7×104m/s，-6.7×104m/s)的氘離子在不同負(fù)徑向電場下氘離子的運(yùn)動軌跡在小環(huán)方向上的投影及速度合成。圖4中軌跡1為無負(fù)徑向電場下的氘離子軌跡，初始入射角大，從初始位置開始沿磁力線負(fù)方向做回旋運(yùn)動，受向上的漂移速度vDB的影響，形成內(nèi)翻香蕉粒子軌跡。負(fù)徑向電場為15 000 V/m時，粒子運(yùn)動速度增加，形成壓縮的通行軌跡。當(dāng)負(fù)徑向電場進(jìn)一步增加時，依然為壓縮的通行軌跡。

圖5示出初始位置(x0，y0，z0)=(2.227 5 m，0，0)，初始速度(vx，vy，vz)=(0，-1.65×105m/s，0)的氘離子在不同負(fù)徑向電場下氘離子的運(yùn)動軌跡在小環(huán)方向上的投影及速度合成。圖5中軌跡1為無負(fù)徑向電場下的氘離子軌跡，初始入射角小，為壓縮的通行粒子軌跡。加入負(fù)徑向電場時，依然為壓縮的通行軌跡。

2.2 氘離子俘獲率計算

托卡馬克中，帶電粒子運(yùn)動速度整體上滿足麥克斯韋分布，其分布函數(shù)形式為：

1——E0=0 V/m；2——E0=15 000 V/m；3——E0=30 000 V/m圖4 初始速度(vx，vy，vz)=(-1.34×105 m/s，-6.7×104 m/s，-6.7×104 m/s)時氘離子運(yùn)動軌跡在小環(huán)方向上的投影(a)及速度合成(b)Fig.4 Projection of poloidal cross-section of deuterium ion orbit (a) and velocity compose (b) under initial velocity (vx, vy, vz)=(-1.34×105 m/s, -6.7×104 m/s, -6.7×104 m/s)

1——E0=0 V/m；2——E0=15 000 V/m圖5 初始速度(vx，vy，vz)=(0，-1.65×105 m/s，0)時氘離子運(yùn)動軌跡在小環(huán)方向上的投影(a)及速度合成(b)Fig.5 Projection of poloidal cross-section of deuterium ion orbit (a) and velocity compose (b) under initial velocity (vx, vy, vz)=(0, -1.65×105 m/s, 0)

(10)

其中：f(v)為單位速率區(qū)間內(nèi)粒子出現(xiàn)的概率；k為玻耳茲曼常數(shù)；T為氘離子對應(yīng)位置處平均等離子體的溫度。

用蒙特卡羅方法對粒子速度進(jìn)行抽樣，獲得整體上滿足麥克斯韋分布的100 000個粒子速度樣本，然后跟蹤每個粒子的運(yùn)動軌跡，統(tǒng)計俘獲粒子的個數(shù)，得到粒子俘獲率。

圖6示出不同初始速度及不同負(fù)徑向電場下的氘離子俘獲率。圖6中初始位置(x0，y0，z0)=(2.227 5 m，0，0)。對比相同初始速度、不同負(fù)徑向電場下的氘離子俘獲率可看出：隨負(fù)徑向電場的增大，離子的俘獲率逐漸降低，當(dāng)負(fù)徑向電場很大時，俘獲率很低，基本都是通行軌跡。對比不同初始速度下負(fù)徑向電場引起的氘離子俘獲率可看出：初始速度小時，負(fù)徑向電場增加引起的氘離子俘獲率的變化大。這是因?yàn)槌跏妓俣刃r，沿磁力線方向的速度vm小，在負(fù)徑向電場不變的情況下，電漂移速度不變，合成速度的相對變化大，所以俘獲率變化大。另外，對比負(fù)徑向電場為0處的俘獲率可看出，平均溫度為0.46 keV的氘離子的俘獲率略低于0.3 keV的氘離子。這是因?yàn)槟芰吭礁?，粒子在邊界損失率越高。另一方面，損失粒子占的份額小，所以在圖6中顯示出略低的趨勢。

圖6 不同初始速度、不同負(fù)徑向電場下的氘離子俘獲率Fig.6 Capture rate of deuterium ion in different initial velocities and different negative radial electric fields

3 結(jié)論

根據(jù)給定的總電流、中心電流密度、中心密度溫度和邊緣密度溫度等實(shí)驗(yàn)參數(shù)求解Grad-Shafranov方程，得到了EAST典型放電實(shí)驗(yàn)(炮號33068)的平衡位形。結(jié)合粒子運(yùn)動方程得到了考慮負(fù)徑向電場情況下氘離子在該平衡位形中的運(yùn)動軌跡，統(tǒng)計了不同初始速度、不同負(fù)徑向電場的俘獲率。模擬結(jié)果表明：隨負(fù)徑向電場的增大，氘離子軌跡由擴(kuò)張軌跡向外翻香蕉軌跡再向內(nèi)翻香蕉軌跡，最后向壓縮軌跡演變；氘離子俘獲率隨負(fù)徑向電場的增大而減小，氘離子初始速度越小，其變化越大。