劉永超,龔學余,*,楊文軍,向 東,曹錦佳

(1.南華大學 電氣工程學院,湖南 衡陽 421001;2.南華大學 核科學技術學院,湖南 衡陽 421001)

0 引 言

能源危機問題是當今最主要的問題，可控核聚變是徹底解決這場危機的有效途徑之一。但可控核聚變需要將等離子體加熱到大約1億攝氏度，常規(guī)的容器不能承受如此高溫，所以需要很好的約束反應物質(zhì)。常用的約束手段有兩種，慣性約束和磁約束。慣性約束是通過強激光電子束對聚變材料進行作用，在極短的時間內(nèi)讓物質(zhì)達到高溫高密度狀態(tài)，從而發(fā)生聚變反應。磁約束聚變是利用強磁場來約束高溫等離子體態(tài)的聚變核燃料，持續(xù)可控的發(fā)生核反應釋放能量，與瞬間釋放能量的聚變方式相比該方案更有利于商用發(fā)電[1]。受控磁約束聚變最具代表性的位形有：托卡馬克和仿星器。而托卡馬克位形被認為是最具前景的。目前國際上在建的ITER[2](international thermonuclear experimental reactor)以及國內(nèi)正在立項設計的CFETR[3](China fusion engineering test reactor)，正是采用的托卡馬克位形。

由于托卡馬克裝置是軸對稱系統(tǒng)，為了計算方便，利用其環(huán)向?qū)ΨQ性，L.L.LAO等人采用的二維模型[4]。袁保山等人對鐵芯結(jié)構(gòu)的托卡馬克裝置，從研究等離子體區(qū)域極向場的角度考慮，近似把中心柱看作無限長，研究了鐵芯的存在對極向場的影響[5]。E.R.Solano等人建立了卷軸模型，計算了鐵芯在非飽和情況下對極向場的影響，給出了極向場的解析表達式[6]。杜世俊等人在研究鐵芯對平衡計算的影響時，把鐵芯對平衡的影響看作是外部極向場線圈鏡像電流的作用，巧妙地將問題轉(zhuǎn)化為空氣鐵芯的磁流體力學(magnetohydrodynamics，MHD)方程的求解[7]。以上結(jié)構(gòu)都忽略環(huán)向不對稱性，視為環(huán)形對稱，皆可看成二維結(jié)構(gòu)。程際等人建立帶鐵芯的極向磁場線圈三維數(shù)值模型，計算并研究鐵芯托卡馬克的三維磁場的極向分量在環(huán)向上的不對稱性[8]。

理想的托卡馬克磁場位形是二維結(jié)構(gòu)，主要包含環(huán)向磁場和極向磁場。但是，電流母線接頭、環(huán)向場線圈和極向場線圈的安裝偏差以及鐵磁性質(zhì)的物質(zhì)等等都會使磁場不再對稱，產(chǎn)生徑向磁場(誤差場)，這會使得磁場位形變成三維結(jié)構(gòu)[9-10]。誤差場的存在會對托卡馬克裝置產(chǎn)生極大影響，例如導致撕裂模鎖模，甚至引發(fā)破裂[11]。故而研究磁場位形的三維結(jié)構(gòu)顯得尤其重要。

J-TEXT(joint texas experimental tokamak)[12]是典型的圓截面鐵芯中型托卡馬克裝置，其大半徑R=1.05 m，等離子體小環(huán)半徑約0.35 m，縱場最大磁感應強度Bt=3.0 T，等離子體最大電流I=400 kA，等離子體中心密度為5.8×1019m-3，中心電子溫度0.85 keV，中心離子溫度0.45 keV。

本論文主要是基于J-TEXT平衡參數(shù)，利用VMEC程序求解得到J-TEXT裝置的三維磁場位形。第1節(jié)介紹了VMEC程序的工作原理，給出了等離子體在磁面坐標下的Grad-Shafranov方程。第2節(jié)中通過VMEC程序仿真得到結(jié)果，論證VEMC程序嵌入J-TEXT裝置的可行性。第3節(jié)分析了等離子體壓強、等離子體電流對三維磁場的影響。

1 物理模型

1.1 VMEC程序簡介

1.2 平衡方程

容易得到，大半徑坐標系(R,φ,Z)下軸對稱模型的等離子體平衡Grad-Shafrov方程[14-15]：

(1)

為了獲得通量坐標系下等離子體平衡Grad-Shafrov方程，人們通常用標準方法求解網(wǎng)格上的(R,φ,Z)的Grad-Shafrov方程，然后映射得到(ψ,θ,ζ)坐標下的Grad-Shafrov方程。

通量坐標系下的等離子體平衡Grad-Shafrov方程[16-17]：

(2)

2 VMEC程序仿真平衡參數(shù)與結(jié)果分析

2.1 平衡參數(shù)

VMEC是三維平衡求解器，在提供初始等離子體邊界、等離子體壓強剖面、旋轉(zhuǎn)變換剖面和外部電流的情況下，程序?qū)⒃诤侠淼臅r間內(nèi)計算整個環(huán)形等離子體的平衡位形。

本文選取J-TEXT裝置第1059184炮平衡文件數(shù)據(jù)，基于平衡文件中壓強剖面數(shù)據(jù)作為VMEC程序的輸入?yún)?shù)，如圖1所示。等離子體壓強在磁軸位置最大，向邊界遞減，在邊界處趨于零。

由于旋轉(zhuǎn)變換τ與安全因子q的關系為：τ～1/q。旋轉(zhuǎn)變換τ可以通過安全因子q來計算求得。同壓強剖面，基于平衡文件中安全因子剖面數(shù)據(jù)作為VMEC程序的輸入?yún)?shù)，如圖1所示。安全因子q的徑向剖面通常在磁軸處或靠近磁軸處有最小值，并向外增大。在大縱橫比和圓形截面的情況下，q簡單地由環(huán)向電流密度剖面j(r)決定。因旋轉(zhuǎn)變換τ與安全因子q成反比例關系，則旋轉(zhuǎn)變換的徑向剖面通常在磁軸處或靠近磁軸處有最大值，并向外遞減。

J-TEXT裝置是一個圓形截面的托卡馬克裝置，根據(jù)平衡文件，大半經(jīng)R=1.05 m，小半徑a=0.35 m。VMEC程序利用一組傅里葉邊界系數(shù)確定其邊界[18]：

式中，RBC，ZBC，RBS，ZBS為傅里葉邊界系數(shù)；m，n分別為極向模數(shù)和環(huán)向模數(shù)；u，v分別為極向角和環(huán)向角。

外部電流的設置適用于自由邊界問題，本文討論的固定邊界問題，故將外部電流設置為固定值。

2.2 結(jié)果分析

圖2(a)為三維磁場邊界輪廓分布模型。從圖2中可以看出，圓形截面，螺旋結(jié)構(gòu)磁場位形，磁場強度由內(nèi)向外(靠近和遠離圓環(huán)中心)逐級遞減。最大磁場約為2 T，磁軸磁場約為1.5 T，基本接近J-TEXT第1059184炮磁場大小。圖2(b)是徑向磁場BR在小截面上的分布。徑向磁場BR是極向磁場Bθ的徑向分量。從圖中可以看出，BR關于Z=0軸對稱，同時滿足在Z=0時等于零。在極向方向0～2π，BR先反向增大后遞減為零，再正向增大后減小為零。圖2(d)是軸向磁場BZ在小截面上的分布。軸向磁場BZ是極向磁場Bθ的軸向分布，與徑向磁場BR共同組成極向磁場。從圖2中可以看出，BZ關于Z=0軸對稱。在Z=0軸上，從大環(huán)內(nèi)側(cè)向大環(huán)外側(cè)逐級遞增。在極向方向0～2π，BZ先減小再增大。圖2(c)是環(huán)向磁場Bφ在小截面上的分布。環(huán)向磁場的基本位形可由安培環(huán)路定律得到，它與大半經(jīng)R呈反比。從圖中可以看出，Bφ關于Z=0軸對稱。在Z=0軸上，從大環(huán)內(nèi)側(cè)向大環(huán)外側(cè)逐級遞減，與J-TEXT第1059184炮實驗數(shù)據(jù)基本吻合。在極向方向0～2π，BZ先增大再減小。

圖2 基于J-TEXT第1 059 184炮平衡數(shù)據(jù)，三維磁場位形分布Fig.2 Based on the equilibrium data of J-TEXT No.1,059,184 gun, the three-dimensional magnetic configuration distribution

圖3主要是徑向磁場BR，軸向磁場BZ和環(huán)向磁場Bφ在小截面的分析。為分析磁場的環(huán)向分布，取磁面r=1(磁軸附近)和r=49(大環(huán)外側(cè)邊界附近)，在極向角θ=0°，計算BR、BZ、Bφ隨環(huán)向角φ的變化。

圖3(a),3(b)是BR在極向角θ=0°，磁軸和邊界的隨環(huán)向角φ變化的圖。從圖中可以看出，BR在磁軸和邊界處的變化規(guī)律相同，都是先遞減再遞增再遞減。在φ=90°和φ=270°附近取得最值，并且φ=180°截面對稱。圖3(c),3(d)是Bφ在極向角θ=0°，磁軸和邊界的隨環(huán)向角φ變化的圖。從圖中可以看出，Bφ在磁軸和邊界處的變化規(guī)律相同，都是先遞增再遞減。在φ=180°處取得最值，并且關于φ=180°截面對稱。圖3(e),(f)是BZ在極向角θ=0°，磁軸和邊界的隨環(huán)向角φ變化的圖。從圖3中可以看出，BZ在磁軸和邊界處的變化規(guī)律相同，都是先遞增再遞減。在φ=180°處取得最值，并且φ=180°關于截面對稱。

3 等離子體參數(shù)對三維磁場的影響

3.1 等離子體電流對三維磁場的影響

為分析等離子體電流對三維磁場的影響，選取等離子體電流Ip=-100 kA，Ip=-220 kA兩種電流在壓強剖面不變的情況下，分別計算徑向磁場BR，軸向磁場BZ和環(huán)向磁場Bφ在環(huán)向位置φ=0°的小截面分布，如圖4所示。

結(jié)合圖2的小截面分布圖，可以發(fā)現(xiàn)的極向方向，的徑向方向磁場變化較為明顯，而變化不明顯。為更加清楚的比較，分解得到在和極向方向分布圖(圖5)以及在和徑向分布圖(圖6)。

圖3 極向角度為零時，BR、Bφ、BZ在磁軸和邊界附近隨環(huán)向角的變化圖Fig.3 When the pole angle is zero, BR, Bφ and BZ near the magnetic axis and the boundary change with the ring angle

圖4 不同等離子體電流下BR、Bφ、BZ小截面分布Fig.4 Distribution of small cross-sections of BR, Bφ and BZ under different plasma currents

圖5 環(huán)向角φ=0°，不同等離子體電流下BR在中部以及邊界附近的極向分布；Fig.5 Circumferential angle φ=0°, the polar distribution of BR in the middle and near the boundary under different plasma currents

圖5(a)，(b)分別是BR在r=25(中間)附近以及在磁面r=49(邊界)，φ=0°，不同壓強下的極向分布圖。從圖中可以看出，在磁面r=25(中心)附近以及在磁面r=49(邊界)附近，Z軸上方，BR在Ip=-220 kA比Ip=-160 kA達到最值要慢，但是達到最值的值要大，而在Ip=-100 kA比Ip=-160 kA達到最值要快，但是最值的值較小。Z軸下方，BR在Ip=-220 kA比Ip=-160 kA達到最值要快，且達到最值的值也要大，而在Ip=-100 kA比Ip=-160 kA達到最值要慢，且最值的值也較小。由此可見，等離子電流越大，在等離子體中間靠邊界處，BR的極向方向變化會變緩，但是數(shù)值會變大。

圖6 環(huán)向角φ=0°時，在不同等離子體電流下BZ在極向角θ=0°及θ=180°的徑向分布Fig.6 Circumferential angle φ=0°, the radial distribution of BZ at the polar angle θ=0° and the polar angle θ=180° under different plasma currents

圖6(a)，(b)分別是BZ在θ=0°以及θ=180°，φ=0°，不同等離子體電流下的徑向分布圖。從圖中可以看出，在θ=0°附近，由中心向邊界沿著徑向方向，BZ在Ip=-220 kA比Ip=-160 kA遞增地更高，而在Ip=-100 kA比Ip=-160 kA遞增低一些。在θ=180°附近，由中心向邊界沿著徑向方向，BZ在Ip=-220 kA比Ip=-160 kA遞減地更快，而在Ip=-100 kA比Ip=-160 kA遞減慢一些。由此可見，在Z軸的左側(cè)，等離子體電流越大，遞增越大；在Z軸右側(cè)，等離子體電流越大，遞減越大。

等離子體電流產(chǎn)生極向磁場以及少部分的環(huán)向磁場。隨著Ip增大，極向磁場分量BR，BZ也隨之增大。BR在θ=0°附近極向方向變化緩慢，在θ=180°附近極向方向變化變快。BZ在磁軸左側(cè)變化幅度變小，磁軸右側(cè)變化幅度變大。而環(huán)向磁場Bφ基本保持不變。

3.2 等離子體壓強對磁場的影響

為分析等離子體壓強對三維磁場的影響，選取磁軸處的最大壓強p0=2 100 Pa和p0=8 400 Pa兩種壓強，如圖7(a)、7(b)所示。在等離子體電流Ip=-160 kA不變的情況下，分別計算徑向磁場BR，軸向磁場BZ和環(huán)向磁場Bφ在環(huán)向位置φ=0°的小截面分布，如圖8所示。

圖7 兩種不同的壓強剖面圖Fig.7 Two different pressure profiles

圖8 不同等離子體壓強下BR、Bφ、BZ小截面分布Fig.8 Distribution of small cross-sections of BR, Bφ and BZ under different plasma pressures

結(jié)合圖2的小截面分布圖，可以發(fā)現(xiàn)BR的極向方向，BZ的徑向方向磁場變化較為明顯，而Bφ變化不明顯，但磁軸向大環(huán)外側(cè)偏移。為更加清楚的比較，分解得到BR在r=25和r=49極向方向分布圖(圖9)以及BZ在θ=0°和θ=180°徑向分布圖(圖10)。

圖9(a)、圖9(b)分別是BR在r=25(中間)附近以及在磁面r=49(邊界)，φ=0°，不同壓強下的極向分布圖。從圖9可以看出，在磁面r=25(中心)附近以及在磁面r=49(邊界)附近，Z軸上方，BR在p0=8 400 Pa比p0=4 200 Pa更快達到最值，而在p0=2 100 Pa比p0=4 200 Pa要慢一些到達。Z軸下方，BR在p0=8 400 Pa比p0=4 200 Pa要慢，p0=2 100 Pa比p0=4 200 Pa更快達到最值。由此可見，壓強越大，在等離子體中間靠邊界處，BR的極向方向變化更陡峭。

圖9 環(huán)向角φ=0°，在不同等離子體壓強下BR在中部以及邊界附近的極向分布Fig.9 Circumferential angle φ=0°, the polar distribution of BR in the middle and near the boundary under different plasma pressures

圖10(a)，10(b)分別是BZ在θ=0°以及θ=180°，φ=0°，不同壓強剖面下的徑向分布圖。從圖中可以看出，在θ=0°附近，由中心向邊界沿著徑向方向，BZ在p0=8 400 Pa比p0=4 200 Pa遞增地更高，而在p0=2 100 Pa比p0=4 200 Pa遞增低一些。在θ=180°附近，由中心向邊界沿著徑向方向，BZ在p0=8 400 Pa比p0=4 200 Pa遞減地更緩，而在p0=2 100 Pa比p0=4 200 Pa遞減快一些。由此可見，在Z軸的左側(cè)，壓強越大，遞增越大；在Z軸右側(cè)，壓強越小，遞減越大。

圖10 環(huán)向角φ=0°，在不同等離子體壓強下BZ在極向角θ=0°以及極向角θ=180°的徑向分布Fig.10 Circumferential angle φ=0°, the radial distribution of BR at the polar angle θ=0° and the polar angle θ=180° under different plasma pressures

等離子體的一個重要參數(shù)β，它的定義為等離子體壓強與磁場壓強的當量比。隨著p0增大，等離子體趨于環(huán)向膨脹，致使β值變大，磁軸外移。BR在極向方向稍微變化增大，但大致保持原來的形狀。BZ在Z軸右側(cè)徑向方向變化更大，而Bφ磁軸向外偏移，其余基本保持原來形狀。

4 結(jié) 論

本文運用VMEC程序，基于J-TXET裝置第1059184炮平衡文件數(shù)據(jù)，對J-TXET裝置等離子體平衡位形進行計算。計算得到的圓形截面，螺旋結(jié)構(gòu)磁場位形，磁場強度由內(nèi)向外(靠近和遠離圓環(huán)中心)逐級遞減。結(jié)果表明，得到的磁場位形以及磁場強度基本符合J-TEXT裝置實驗數(shù)據(jù)，故可以初步論證VMEC程序可以滿足J-TEXT裝置三維磁場平衡位形的計算。

本文進一步分析了等離子體電流和壓強剖面對三維磁場的影響。等離子體電流直接影響極向磁場的強度，環(huán)向磁場影響不大。壓強對三維磁場的強度影響甚小，當?shù)入x子體壓強較大時，等離子體傾向于環(huán)狀膨脹，而當?shù)入x子體壓強降低時，等離子體則會環(huán)狀收縮。本文只討論了固定邊界，外部電流恒定的情況，下一步的工作對于自由邊界問題展開更深入的研究。