一維Fokker-Planck方程的有限體積法求解及其在中性束注入加熱中的應用

張 能，龔學余,，黃千紅，侯 偉，謝寶藝，余江妹，李景春

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

?

一維Fokker-Planck方程的有限體積法求解及其在中性束注入加熱中的應用

張 能1，龔學余1,2，黃千紅2,*，侯 偉2，謝寶藝1，余江妹1，李景春2

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

采用有限體積法數(shù)值求解一維Fokker-Planck方程，通過模擬電子自碰撞過程進行程序校驗。研究表明，有限體積法能高效求解Fokker-Planck方程，能確保分布函數(shù)的非負性和粒子數(shù)密度守恒，同時計算程序能有效克服傳統(tǒng)求解方法中出現(xiàn)的分布函數(shù)對麥克斯韋分布的過沖現(xiàn)象。模擬了HL-2A裝置在中性束注入加熱等離子體過程中，離子分布函數(shù)和溫度隨時間的演化情況。結果表明，隨著中性束的注入，離子分布函數(shù)出現(xiàn)非麥克斯韋化，離子溫度迅速增加后穩(wěn)定；計算結果和實驗結果符合較好。進一步討論束能量和功率的影響，隨束能量和功率的增加，等離子體離子溫度均升高，離子溫度隨束能量的增加升高的幅值較大，而隨束功率的增加升高的幅值較小。

有限體積法；Fokker-Planck方程；分布函數(shù)演化；中性束注入加熱

中性束注入(neutral beam injection，NBI)加熱的物理機制較明晰且具備加熱效率高等特點，已成為高溫等離子體輔助加熱最有效的方法之一，目前在國內(nèi)外的主要托卡馬克裝置中已得到廣泛應用[1-4]。在HL-2A托卡馬克裝置中，已進行了不同參數(shù)下的中性束注入加熱實驗，并取得很好的加熱效果。2009年，在中性束注入和電子回旋共振加熱(ECRH)條件下實現(xiàn)了國內(nèi)的首次H模放電[5]；此后，借助中性束注入加熱系統(tǒng)，開展了一系列重要的等離子體物理實驗研究，包括等離子體輸運和約束改善研究[6]、L-H轉換機理研究[7]、高能量粒子物理學和磁流體不穩(wěn)定性研究[8]等；2013年，在1.5 MW中性束注入加熱實驗中，單獨用0.6 MW的束功率同樣獲得了H模放電[9]。因此，對HL-2A裝置進行中性束注入加熱實驗的同時，有必要展開相關的數(shù)值模擬研究，為以后的實驗提供參考。

為更深入了解中性束注入后等離子體溫度的演化及其影響因素，需求解相應的Fokker-Planck方程。Fokker-Planck方程傳統(tǒng)的求解是采用有限差分法進行，計算結果對網(wǎng)格點步長和時間步長均有很嚴格的要求以確保分布函數(shù)的非負性和粒子數(shù)守恒[10-13]；文獻[14]利用有限元法數(shù)值求解一維Fokker-Planck方程，模擬了高斯分布的電子經(jīng)過自碰撞趨于麥克斯韋分布的過程，但已開發(fā)的程序與傳統(tǒng)上用差分方法的求解程序一樣，對麥克斯韋分布存在過沖現(xiàn)象[15]。

本文推導一維Fokker-Planck方程有限體積法計算格式，數(shù)值模擬電子自碰撞由超高斯分布趨于麥克斯韋分布的過程，并進行程序校驗。進一步模擬HL-2A裝置中性束注入加熱過程中，離子分布函數(shù)和溫度隨時間的演化情況，并討論束功率和能量的影響。

1 一維Fokker-Planck方程的有限體積法求解

等離子體中的帶電粒子與德拜球內(nèi)其他粒子相互作用中，會發(fā)生多重小角度碰撞，使帶電粒子速度不斷出現(xiàn)微小的隨機性變化，能準確描述這種粒子系統(tǒng)的動力學方程是Fokker-Planck方程：

(1)

(2)

(3)

(4)

式中：i、j為網(wǎng)格點；n為粒子種類數(shù)；Z為攜帶的電荷數(shù)；lnΛ為庫侖對數(shù)；Hn和Gn為Rosen-bluth勢函數(shù)。

假設等離子體各物理量在空間均勻分布且不受外力作用的影響，采用的粒子分布函數(shù)只與速度的大小和時間有關。Fokker-Planck方程可簡化成：

(5)

其中：

(6)

式中：a、b為粒子種類；上標a/b為粒子a碰撞b；nb為粒子b的密度；q為電荷數(shù)。則Rosenbluth勢函數(shù)表示為：

(7)

(8)

對式(5)使用有限體積法處理，方程兩邊同時積分：

(9)

對速度項積分，得到：

(10)

其中：ΔV為控制容積體積；A為控制容積邊界面積。采用Chang的方法[11]近似表示fi-1/2和fi+1/2：

(11)

(12)

(13)

為驗證計算格式的正確性，模擬速度空間一維情況下電子自碰撞由超高斯分布趨于麥克斯韋分布的過程，首先令電子速度分布函數(shù)為超高斯分布：

(14)

圖1a示出了用有限體積法計算的電子分布函數(shù)隨時間的演化，圖中標記為M的曲線為麥克斯韋分布函數(shù)。由圖可知，開始階段分布函數(shù)趨于麥克斯韋平衡分布的速度較快，經(jīng)很短時間后分布函數(shù)趨于麥克斯韋平衡分布的速度越來越慢，經(jīng)過100個時間步，分布函數(shù)和麥克斯韋分布函數(shù)基本重合，經(jīng)過1 000個時間步，分布函數(shù)幾乎無變化，說明此時基本達到平衡態(tài)。圖1b為有限元方法計算的結果[14]，由圖可知，低能區(qū)域的分布函數(shù)較麥克斯韋分布偏高，高能區(qū)域的分布函數(shù)偏低。這表明本文程序有效克服了計算過程中對麥克斯韋分布存在的過沖現(xiàn)象。圖2為電子數(shù)密度隨時間的演化關系，表明有限體積法能保證電子數(shù)守恒。

a——有限體積法；b——有限元方法

圖2 電子數(shù)密度的演化曲線Fig.2 Evolution curve of electron number density

通過有限體積法求解一維Fokker-Planck方程，數(shù)值模擬了一維情況下電子自碰撞由超高斯分布趨于麥克斯韋分布的過程，證明了計算格式和程序能確保電子分布函數(shù)的非負性和電子數(shù)密度守恒，表明了有限體積法對于求解Fokker-Planck的高效性。

2 HL-2A裝置上中性束注入加熱模擬

在上述高效求解Fokker-Planck方程的基礎上，進一步模擬在中性束注入加熱條件下的一維Fokker-Planck方程：

(15)

其中，S為中性束注入的源項，用高斯函數(shù)來模擬其束流：

(16)

計算采用HL-2A托卡馬克裝置中性束注入加熱實驗參數(shù)：大半徑R=1.65 m，小半徑r=0.4 m；等離子體電流Ip=480 kA，磁場強度Bt=2.8 T，等離子體中離子密度n0=1.5×1019m-3；在中性束注入前由歐姆加熱到的離子溫度(以能量表示)為Ti=400 eV；中性束注入功率為600 kW，能量為40 keV，計算得到的流強J=2.965×1018m-3·s-1。

由于離子與離子碰撞交換能量的時間較離子與電子碰撞短很多，因而中性束注入產(chǎn)生快離子的能量主要交給本底離子。所以在計算中，始終將電子的分布當作麥克斯韋分布來處理。計算離子的溫度隨中性束注入的變化：

(17)

圖3 離子分布函數(shù)的演化曲線Fig.3 Evolution curve of ion distribution function

圖3示出了中性束注入加熱過程中離子分布函數(shù)的演化曲線。結果表明，隨著中性束的注入，低能量的離子數(shù)減少，高能量的離子數(shù)增加，離子分布函數(shù)出現(xiàn)非麥克斯韋化。

圖4示出了中性束注入加熱過程中離子溫度的演化情況。由圖可知，中性束注入加熱等離子體初始階段離子溫度迅速上升，產(chǎn)生的快離子迅速將能量傳遞給本底離子；碰撞耗散和離子吸收能量平衡后，離子能量不再上升，在加熱的后階段，由于高能粒子能量的損失，離子溫度出現(xiàn)緩慢下降，這與實驗中中性束注入加熱等離子體時的離子溫度演化曲線相符[17]。

圖4 離子溫度的演化曲線Fig.4 Evolution curve of ion temperature

圖5為在其他參數(shù)不變的情況下，束功率分別取600、700、800 kW得到的溫度演化曲線。結果表明，提高中性束的注入功率，可提升離子溫度的爬升速度，使等離子體更快達更高溫度，但溫度幅值變化不大。圖6為在其他參數(shù)不變的情況下，束能量分別取40、50、60 keV得到的溫度演化曲線。結果表明，提高中性束的能量，能快速提升中性束的加熱溫度，且溫度幅值有較大的提升。

圖5 不同束功率下離子溫度的演化曲線Fig.5 Evolution curve of ion temperature for different beam powers

圖6 不同束能量下離子溫度的演化曲線Fig.6 Evolution curve of ion temperature for different beam energy

3 結論

本文通過有限體積法數(shù)值求解一維Fokker-Planck方程，數(shù)值模擬一維速度空間超高斯分布電子經(jīng)自碰撞趨于麥克斯韋分布的過程，證明有限體積法能高效求解Fokker-Planck方程。進一步模擬了HL-2A裝置中性束注入加熱的情況，得到了離子分布函數(shù)隨中性束注入出現(xiàn)非麥克斯韋化。計算了離子溫度的演化曲線，與實驗結果符合較好。討論了中性束的功率和能量的影響，結果表明，隨著束功率和能量的升高，離子溫度升高，相較而言，隨束能量的增大離子溫度升高的幅值較大。由于數(shù)值模擬中性束注入加熱過程中采用了一維速度空間近似，所以現(xiàn)有計算程序及模擬結果存在不足。在中性束注入加熱過程中，需考慮中性束注入角度對加熱效率、功率沉積等的影響，以便研究更多的物理問題，這也是下一步工作要考慮的問題。

[1] HU L Q, KURIYAMA M, AKINO N, et al. Beam divergence and power loading on the beamline components of the negative-ion based NBI system for JT-60U[J]. Fusion Engineering and Design, 2001, 54(2): 321-330.

[2] AKERS R J, APPEL L C, CAROLAN P G, et al. Neutral beam heating in the START spherical Tokamak[J]. Nucl Fusion, 2002, 42(2): 122-135.

[3] CIRIC D, BROWN D P D, CHALLIS C D, et al. Overview of the JET neutral beam enhancement project[J]. Fusion Engineering and Design, 2007, 82(5-14): 610-618.

[4] ZOHM H, ADAMEK J, ANGIONI C, et al. Overview of ASDEX upgrade results[J]. Nucl Fusion, 2009, 49(10): 104009.

[5] DUAN X R, DONG J Q, YAN L W, et al. Preliminary results of ELMy H-mode experiments on the HL-2A Tokamak[J]. Nucl Fusion, 2010, 50(9): 095011.

[6] XIAO W W, ZOU L N, DING X T, et al. Observation of a spontaneous particle-transport barrier in the HL-2A Tokamak[J]. Phys Rev Lett, 2010, 104: 215001.

[7] CHENG J, DONG J Q, ITOH K, et al. Dynamics of low-intermediate-high confinement transition in toroidal plasma[J]. Phys Rev Lett, 2013, 110: 265002.

[8] DUAN X R, DING X T, DONG J Q, et al. An overview of recent HL-2A experiments[J]. Nucl Fusion, 2013, 53(10): 104009.

[9] 嚴龍文，段旭如，丁玄同，等. HL-2A裝置物理實驗研究的進展[C]∥第十六屆全國等離子體科學技術會議. 上海：復旦大學，2013.

[10]MIRIN A A, McCOY M G, TOMACCHKE G P, et al. FPPAC88: A two-dimension multispecies nonlinear Fokker-Planck package[J]. Comput Phys Comm, 1988, 51(3): 373-380.

[11]SHOUCRI M, SHKAROFSKY I. A fast 2D Fokker-Planck solver with synergetic effects[J]. Comput Phys Commun, 1994, 82(2-3): 287-305.

[12]PEYSSON Y, SHOUCRI M. An approximate factorization procedure for solving nine-point elliptic difference equations: Application for a fast 2-D relativistic Fokker-Planck solver[J]. Comput Phys Commun, 1998, 109(1): 55-80.

[13]KARNEY C F F. Fokker-Planck and quasilinear codes[J]. Comput Phys Rep, 1986, 4(3-4): 183-244.

[14]翁蘇明，盛正明. 等離子體中Fokker-Planck方程的有限元模擬[J]. 計算物理，2007，24(2)：134-140.

WENG Suming, SHENG Zhengming. Finite element scheme for Fokker-Planck equation in plasmas[J]. Chinese Journal of Computational Physics, 2007, 24(2): 134-140(in Chinese).

[15]SHKRAROFSKY I P, JOHNSTON T W, BACKYNSKI M P. The particle kinetics of plasmas[M]. Netherlands: Addison-Wesley Publishing Company, 1966: 132-137.

[16]王建國，李有宜，李建剛. HT-6M中性束注入加熱理論分析[J]. 物理學報，1995，44(1)：92-97.

WANG Jianguo, LI Youyi, LI Jiangang. Theoretical study for neutron beam injection heating[J]. Acta Physica Sinica, 1995, 44(1): 92-97(in Chinese).

[17]PAN C H. Overview of the auxiliary heating systems for HL-2A Tokamak[M]∥Southwestern Institute of Physics Annual Report 2007. Chengdu: Southwestern Institute of Physics, 2007.

One-dimensional Fokker-Planck Equation Solving by Finite Volume Method and Its Application in Neutral Beam Injection Heating

ZHANG Neng1, GONG Xue-yu1,2, HUANG Qian-hong2,*, HOU Wei2,XIE Bao-yi1, YU Jiang-mei1, LI Jing-chun2

(1.SchoolofMathematicsandPhysics,UniversityofSouthChina,Hengyang421001,China; 2.SchoolofNuclearScienceandTechnology,UniversityofSouthChina,Hengyang421001,China)

The one-dimensional Fokker-Planck equation was solved by finite volume method. The reliability of program was verified by simulation of electron self-collision. The results show that the numerical scheme assures positive distribution function and particle number conservation. It is an effective method of solution for Fokker-Planck equation. And the overshoot of the Maxwell distribution that arouse from Fokker-Planck equation solving by traditional numerical method was eliminated. Furthermore, the evolutions of ion distribution function and temperature were simulated in neutral beam injection heating on the HL-2A. The results show that the ion distribution function tends to non-Maxwell distribution as the neutral beam injects, and the ion temperature rises rapidly and then reaches stability. The effects of beam power and energy on the ion temperature were investigated, and the results show that the ion temperature increases with the beam power and energy. The change is more obvious by increasing beam energy than power.

finite volume method; Fokker-Planck equation; evolution of distribution function; neutral beam injection heating

2014-12-02;

2015-01-26

國際熱核聚變實驗堆計劃專項資助項目(2014GB108002)；國家自然科學基金資助項目(11375085)；湖南省核聚變與等離子體物理創(chuàng)新團隊建設資助項目(NHXTD03)；衡陽市科技局資助項目(2013KJ24)；南華大學研究生科研創(chuàng)新資助項目(2014XCX08,2014XCX02)

張 能(1989—)，男，湖南桃江人，碩士研究生，從事核聚變與等離子體物理研究

*通信作者：黃千紅，E-mail: hqhhgy@163.com

O532.26

A

1000-6931(2015)06-0966-06

10.7538/yzk.2015.49.06.0966