H-1NF仿星器標(biāo)準(zhǔn)磁場位形分析與高能量離子運(yùn)動軌道模擬

向 東，夏彥博，宋澤宸，曹錦佳，杜 丹，管滬楠，龔學(xué)余

(1.南華大學(xué) 核科學(xué)技術(shù)學(xué)院，湖南 衡陽 421001；2.南華大學(xué) 數(shù)理學(xué)院，湖南 衡陽 421001)

仿星器和托卡馬克裝置是目前國際上研究最多的兩大磁約束可控?zé)岷司圩冄b置類型。由于建造仿星器的成本高、技術(shù)難度大等諸多原因，從20世紀(jì)70年代開始，大多數(shù)核聚變研究都聚焦于托卡馬克。目前中國擁有EAST、HL-2M、J-TEXT等托卡馬克裝置，還沒有仿星器裝置，因而在仿星器方面的理論研究相對較少。然而，由于仿星器內(nèi)等離子體全部依靠外部磁場約束，易于實現(xiàn)自然穩(wěn)態(tài)運(yùn)行，激發(fā)大破裂的可能性非常低，這些優(yōu)點(diǎn)在聚變界一直備受青睞[1-3]。特別是在2015年12月，目前世界上最大、設(shè)計最復(fù)雜、版本最新的德國W7-X仿星器成功實現(xiàn)高參數(shù)等離子體放電而轟動全球[4]。W7-X仿星器設(shè)計的遠(yuǎn)景目標(biāo)是實現(xiàn)1 800 s超長脈沖穩(wěn)態(tài)運(yùn)行[5]。目前，建造中國聚變工程實驗堆(CFETR)已被列為國家聚變能源發(fā)展的路線圖[6]。仿星器作為未來聚變反應(yīng)堆備選方案之一，其相關(guān)研究同樣也受到了國內(nèi)科學(xué)界的關(guān)注。西南交通大學(xué)和日本國立聚變研究所合作在中國設(shè)計籌建一個新型準(zhǔn)軸對稱仿星器CFQS[7-8]，浙江大學(xué)、中國科技大學(xué)等高校聚變研究團(tuán)隊也開始了關(guān)于仿星器的理論研究[9-10]。為了在國內(nèi)廣泛開展H-1NF仿星器的實驗和理論研究，本文對H-1NF仿星器的磁場線圈組成和分布特征進(jìn)行分析，研究其磁場位形特點(diǎn)及單粒子運(yùn)動軌道特征，并模擬計算高能量離子在一種標(biāo)準(zhǔn)磁場位形中的典型運(yùn)動軌道。

1 H-1NF仿星器磁場線圈分布

H-1NF仿星器的磁場沿大環(huán)方向三周期120°對稱分布[11-13]，如圖1所示，其磁場線圈系統(tǒng)共包括36餅環(huán)向場線圈(TFC)，1個中心環(huán)導(dǎo)體，即極向場線圈(PFC)，1條螺旋線圈(HCW)，1對內(nèi)垂直場線圈(IVFC)，1對外垂直場線圈(OVFC)。其中，環(huán)向場線圈、極向場線圈、螺旋線圈和內(nèi)垂直場線圈均安裝在仿星器真空腔體內(nèi)，外垂直場線圈安裝在真空腔體外。

a——三維視圖；b——俯視圖圖1 H-1NF仿星器磁場線圈分布Fig.1 Distribution of magnetic field coil of H-1NF stellarator

36餅環(huán)向場線圈的大小規(guī)格相等，每餅環(huán)向場線圈內(nèi)部繞有10圈銅導(dǎo)線，其平均等效半徑RTFC=0.38 m。這些環(huán)向場線圈的圓心不在同一平面上，若采用大柱坐標(biāo)(R，φ，Z)來描述各線圈的分布位置，則每個環(huán)向場線圈的圓心坐標(biāo)(Rc，φc，Zc)滿足以下關(guān)系式：

Rc=R0+ρscos(3φc)

(1)

Zc=ρssin(3φc)

(2)

φc=φ′-0.009 7sin(3φ′)

(3)

(4)

其中：R0為大半徑，R0=1.00 m；ρs為擺動半徑，ρs=0.22 m。

中心極向場線圈位于Z=0.00 m的赤道面上，其內(nèi)部繞有36圈銅導(dǎo)線，其線圈等效半徑RPFC=1.00 m。1對內(nèi)垂直場線圈分別對稱分布于赤道平面Z=±1.07 m的平面上，每個內(nèi)垂直場線圈內(nèi)部繞有16圈銅導(dǎo)線，其線圈等效半徑RIVFC=0.72 m。1對外垂直場線圈分別對稱分布于赤道平面Z=±0.70 m的平面上，其內(nèi)部繞有12圈(通過外部引線可設(shè)置成4、8、12圈3種，其標(biāo)準(zhǔn)配置為8圈)銅導(dǎo)線，其線圈等效半徑ROVFC=2.13 m。螺旋線圈內(nèi)部并排4條銅導(dǎo)線，它們沿環(huán)向呈3周期纏繞中心極向場線圈，纏繞的等效小半徑rHCW=0.10 m。

H-1NF仿星器線圈標(biāo)準(zhǔn)接線方式是將環(huán)向場線圈、極向場線圈、內(nèi)垂直場線圈、外垂直場線圈中的銅導(dǎo)線全部串聯(lián)，由ABB電源的主回路供電，而螺旋線圈由ABB電源的從回路供電。通過設(shè)置ABB電源主、從兩路輸出的電流強(qiáng)度，可獲得不同的磁場位形。H-1NF仿星器標(biāo)準(zhǔn)運(yùn)行模式是外垂直場線圈設(shè)置為8圈銅導(dǎo)線，且螺旋線圈不供電。在標(biāo)準(zhǔn)運(yùn)行模式下，當(dāng)主路輸入的銅導(dǎo)線電流強(qiáng)度為13 889 A時，若將中心極向場線圈等效看作是一圈圓電流，則其電流強(qiáng)度約為5.0×105A，此時磁軸上的磁感應(yīng)強(qiáng)度約為1 T。本文以此種H-1NF仿星器標(biāo)準(zhǔn)運(yùn)行模式下磁場位形為例，分析其位形特點(diǎn)，并模擬計算研究高能量離子在該標(biāo)準(zhǔn)磁場位形中的典型運(yùn)動軌道特征。

2 標(biāo)準(zhǔn)配置下的磁場位形

H-1NF仿星器運(yùn)行時，磁面將圍繞中心極向場線圈沿環(huán)向呈三周期麻花狀分布，令最外層磁面為等離子體邊界，繪制出其等離子體位形如圖2所示?；贖-1NF仿星器標(biāo)準(zhǔn)運(yùn)行模式參數(shù)，采用HELIAC程序追蹤計算各磁面的磁場線[13]，分別繪制出不同環(huán)向角φ=0°、15°、30°、45°、60°、75°、90°、105°、120°等真空磁場的龐加萊圖，如圖3所示。由圖3可看出，仿星器真空磁面沿環(huán)向螺旋狀分布，這與托卡馬克磁面為環(huán)向?qū)ΨQ分布有很大不同，且仿星器的磁軸不在一平面上。另外，圖3a和圖3i顯示出φ=0°和120°磁面的極向截面完全相同，這說明H-1NF仿星器的磁場分布是沿大環(huán)方向3周期120°對稱。

圖2 H-1NF仿星器內(nèi)等離子體位形Fig.2 Plasma shape in H-1NF stellarator

雖然仿星器磁面沿環(huán)向螺旋狀分布不具有環(huán)向?qū)ΨQ性，但是不同環(huán)向角的磁面極向截面除了沿磁軸旋轉(zhuǎn)一定角度外，其形貌特征相近。因此，基于磁軸的位置和磁面沿環(huán)向角的變化規(guī)律，以磁軸為旋轉(zhuǎn)軸，將圖3所示的不同環(huán)向角的極向截面的磁場龐加萊圖在原有截面坐標(biāo)(R，Z)內(nèi)，圍繞磁軸旋轉(zhuǎn)一定的角度，并將坐標(biāo)原點(diǎn)平移至磁軸的位置，定義旋轉(zhuǎn)后的極向截面坐標(biāo)為(R′，Z′)，此截面稱為旋轉(zhuǎn)坐標(biāo)系下的等效標(biāo)定極向截面。因此得到了一種旋轉(zhuǎn)坐標(biāo)系下等效標(biāo)定極向截面上的磁場龐加萊圖，圖3所示的各種極向截面的磁場龐加萊圖經(jīng)過旋轉(zhuǎn)變換后對應(yīng)的等效標(biāo)定極向截面的磁場龐加萊圖如圖4所示。

在H-1NF仿星器標(biāo)準(zhǔn)運(yùn)行模式下，螺旋線圈不通電，而環(huán)向場線圈、極向場線圈、內(nèi)垂直場線圈和外垂直場線圈均可看作是一個等效的單匝圓形線圈導(dǎo)體。圓形導(dǎo)體產(chǎn)生的磁場可用第二類橢圓積分來表示，而這種積分可采用Hastings型多項式來計算[14]。因此，基于圓形導(dǎo)體產(chǎn)生磁場的計算方法和場強(qiáng)疊加原理，可實現(xiàn)H-1NF仿星器內(nèi)任一點(diǎn)的磁感應(yīng)強(qiáng)度計算。根據(jù)這一原理，HELIAC程序計算磁場以實現(xiàn)任一磁場線的追蹤。本文采用此方法，利用Matlab編程計算出環(huán)向角φ=0°時極向截面的磁場分布特征(圖5)。圖5a為旋轉(zhuǎn)坐標(biāo)系下環(huán)向角φ=0°時等效標(biāo)定極向截面的磁場龐加萊圖；圖5b為圖5a中示意磁面(極向截面為一個環(huán))隨極向角變化的磁感應(yīng)強(qiáng)度分布；圖5c為旋轉(zhuǎn)坐標(biāo)系下Z′=0直線上沿R′方向的磁感應(yīng)強(qiáng)度分布。

a——φ=0°；b——φ=15°；c——φ=30°；d——φ=45°；e——φ=60°； f——φ=75°；g——φ=90°；h——φ=105°；i——φ=120°圖3 不同環(huán)向角的極向截面的磁場龐加萊圖Fig.3 Poincaré map of magnetic field in poloidal cross-section with different toroidal angles

3 高能量離子的典型運(yùn)動軌道模擬

忽略重力場等次要因素，離子在真空磁場B中的運(yùn)動遵循動力學(xué)方程：

(5)

其中：m、q和v分別為離子的質(zhì)量、電荷和運(yùn)動速度；t為運(yùn)動時間。

計算運(yùn)動軌道的數(shù)值算法采用四階龍格庫塔算法[15]。另外，數(shù)值精度還與每個離散時間步Δt所對應(yīng)的空間尺寸有關(guān)，因此計算前必須估算離子在計算區(qū)域內(nèi)磁場中的拉莫爾半徑(rL=mv⊥/qB)，以保證離子在磁場中運(yùn)動1個回旋周期的數(shù)值模擬經(jīng)過了足夠多的時間步數(shù)演算。在磁場中運(yùn)動離子的漂移速度與離子的能量、荷質(zhì)比、投擲角和背景磁場均有關(guān)系。若離子的初始位置不同，其所在位置的磁場就可能不同，且在運(yùn)動中由于離子位置處于不停地變化中，因此其軌道上每一點(diǎn)的磁場會不停地變化，投擲角也會不斷地改變。對于低能量離子，其漂移速度很小，回旋運(yùn)動的導(dǎo)心基本會沿磁力線走，因此，以下主要研究高能量離子的運(yùn)動軌道特征。

a——φ=0°；b——φ=15°； c——φ=30°；d——φ=45°；e——φ=60°； f——φ=75°；g——φ=90°；h——φ=105°；i——φ=120°圖4 旋轉(zhuǎn)坐標(biāo)系下不同環(huán)向角的等效標(biāo)定極向截面的磁場龐加萊圖Fig.4 Poincaré map of magnetic field in equivalent calibrated poloidal cross-section with different toroidal angles in rotating coordinate system

圖5 環(huán)向角φ=0°時極向截面的磁場分布特征圖5 Distribution characteristic of magnetic field in poloidal cross-section with toroidal angle φ=0°

不失一般性，本文令測試離子為質(zhì)子，能量為10 keV，假設(shè)其初始位置的導(dǎo)心在環(huán)向角為φ=0°時圖5a所示磁面上的點(diǎn)A，其極向截面坐標(biāo)為(R=1.30 m，Z=0.00 m)，相應(yīng)旋轉(zhuǎn)坐標(biāo)系中等效標(biāo)定極向截面坐標(biāo)為(R′=0.06 m，Z′=0.00 m)。令離子的初始投擲角為η，則cosη=v∥/v，v∥為離子速度平行磁場方向的分量。通行粒子軌道和捕獲粒子軌道間的臨界投擲角ηc滿足：

(6)

其中：Bs為離子初始位置的磁感應(yīng)強(qiáng)度；Bm為離子運(yùn)動軌道上的最大磁感應(yīng)強(qiáng)度。在0°～90°范圍內(nèi)，當(dāng)離子的初始投擲角η<ηc時為通行粒子，反之為捕獲粒子。根據(jù)圖5b的磁場數(shù)據(jù)可初步理論估計ηc，但因為仿星器磁場分布不具有環(huán)向?qū)ΨQ性，在不同環(huán)向角對應(yīng)的旋轉(zhuǎn)坐標(biāo)系中標(biāo)定極向截面上的磁感應(yīng)強(qiáng)度沿極向角分布會有一定的差異，所以實際的ηc與根據(jù)圖5b簡單理論的估計值不同。經(jīng)過精細(xì)計算，得到ηc=57°。若令離子的初始投擲角η=30°，則模擬計算出其運(yùn)動全軌道如圖6所示，顯然這是一種典型的通行粒子軌道。若將離子的初始投擲角增加到η=60°，則可得到一種典型的捕獲粒子軌道，其模擬計算出的全軌道如圖7所示。

由于仿星器磁面不具有環(huán)向?qū)ΨQ性，因而無法像研究托卡馬克中的粒子軌道一樣將三維粒子軌道直接投影到二維極向截面上觀察。從圖6、7中無法直接觀察離子運(yùn)動軌道相對于磁面的位置，也不能直接判斷離子是否飛出最外層閉合磁面。然而，H-1NF仿星器內(nèi)不同環(huán)向角的極向截面沿磁軸旋轉(zhuǎn)一定角度后形貌特征相近(圖4)。因此，基于磁軸位置和和磁面沿環(huán)向角的旋轉(zhuǎn)變化規(guī)律，可將高能量離子的三維運(yùn)動軌道投影到類似圖4所示的旋轉(zhuǎn)坐標(biāo)系中不同環(huán)向角的等效標(biāo)定極向截面上，從而可更加直觀清晰地顯示出高能量離子在仿星器磁場位形中的運(yùn)動軌道特征。雖然旋轉(zhuǎn)坐標(biāo)系中不同環(huán)向角的等效標(biāo)定極向磁面不完全相同(圖4a～i)，但是差別不大。為了簡化，統(tǒng)一取φ=0°的標(biāo)定極向磁面作為背景，分別繪制出圖6、7所示的三維粒子軌道對應(yīng)在旋轉(zhuǎn)坐標(biāo)系中等效標(biāo)定極向截面的二維軌道投影(圖8、9)。從圖8可看出，仿星器中的通行粒子軌道與托卡馬克中的通行粒子軌道特征相似，但仿星器中的通行粒子在旋轉(zhuǎn)等效標(biāo)定極向截面上繞幾圈后才閉合。從圖9可看出，仿星器中的捕獲粒子與托卡馬克中的捕獲粒子類似地形成了香蕉軌道，但在旋轉(zhuǎn)等效標(biāo)定極向截面上仿星器中香蕉軌道并未閉合，且軌道逐漸向外側(cè)漂移，最終有可能飛出最外層閉合磁面導(dǎo)致粒子損失。

a——三維視圖；b——俯視圖圖6 高能量通行粒子全軌道模擬圖Fig.6 Full-orbit simulation of energetic passing particle

a——三維視圖；b——俯視圖圖7 高能量捕獲粒子全軌道模擬圖Fig.7 Full-orbit simulation of energetic trapped particle

圖8 高能量通行粒子軌道在旋轉(zhuǎn)坐標(biāo)系下的 等效標(biāo)定極向截面上的軌道投影Fig.8 Orbit projection of energetic passing particle on equivalent calibrated poloidal cross-section in rotating coordinate system

采用上述方法對不同初始位置的測試離子進(jìn)行大量試算，以進(jìn)一步證明此方法的有效性。例如，設(shè)置測試離子的初始位置的導(dǎo)心是環(huán)向角為φ=0°時圖5a所示磁面上的點(diǎn)B，其旋轉(zhuǎn)坐標(biāo)系中的標(biāo)定極向截面坐標(biāo)為(R′=0.00 m，Z′=-0.12 m)，初始投擲角為η=60°，則該測試離子在旋轉(zhuǎn)標(biāo)定極向截面的軌道投影如圖10所示。雖然該測試離子的初始投擲角仍為60°，但其所在初始位置的磁感應(yīng)強(qiáng)度變大，則其臨界投擲角ηc相應(yīng)變大，從而滿足條件η<ηc，此時該測試離子就變?yōu)橥ㄐ辛Ｗ印?/p>

圖9 高能量捕獲粒子軌道在旋轉(zhuǎn)坐標(biāo)系下的 等效標(biāo)定極向截面上的軌道投影Fig.9 Orbit projection of energetic trapped particle on equivalent calibrated poloidal cross-section in rotating coordinate system

圖10 高能量粒子軌道在旋轉(zhuǎn)坐標(biāo)系下的 等效標(biāo)定極向截面上的軌道投影Fig.10 Orbit projection of energetic particle on equivalent calibrated poloidal cross-section in rotating coordinate system

4 結(jié)論

本文采用HELIAC程序計算了H-1NF仿星器在標(biāo)準(zhǔn)運(yùn)行模式下不同環(huán)向角的真空磁場龐加萊圖，研究其磁場位形特點(diǎn)，并模擬計算出高能量離子在該標(biāo)準(zhǔn)磁場位形中典型的通行粒子和捕獲粒子三維運(yùn)動全軌道?；贖-1NF仿星器標(biāo)準(zhǔn)磁場位形的磁軸位置和磁面沿環(huán)向角的旋轉(zhuǎn)變化規(guī)律，以磁軸為旋轉(zhuǎn)軸，按照旋轉(zhuǎn)規(guī)律將不同環(huán)向角的極向截面旋轉(zhuǎn)后得到旋轉(zhuǎn)坐標(biāo)系下的等效標(biāo)定極向截面，并將高能量離子的三維運(yùn)動軌道投影到此等效標(biāo)定極向截面上，從而能更加清晰地顯示出高能量離子在該磁場位形中的運(yùn)動軌道特征。結(jié)果表明，H-1NF仿星器中的通行粒子軌道特征與一般托卡馬克中的通行粒子軌道特征相似，不同之處是仿星器中的通行粒子軌道在等效標(biāo)定極向截面上繞幾圈后才閉合；H-1NF仿星器中的捕獲粒子類似于一般托卡馬克中的捕獲粒子形成的香蕉軌道特征，但H-1NF仿星器中香蕉軌道并未閉合，且軌道逐漸向外側(cè)漂移，最終可能導(dǎo)致粒子損失。

感謝澳大利亞國立大學(xué)Boyd Blackwell博士和Clive Michael博士為本工作提供的HELIAC程序代碼及建議。