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

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

?

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

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

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

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

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

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

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

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

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

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

(1)

(2)

(3)

(4)

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

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

(5)

其中：

(6)

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

(7)

(8)

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

(9)

對(duì)速度項(xiàng)積分，得到：

(10)

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

(11)

(12)

(13)

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

(14)

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

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

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

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

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

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

(15)

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

(16)

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

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

(17)

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

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

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

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

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

圖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 結(jié)論

本文通過(guò)有限體積法數(shù)值求解一維Fokker-Planck方程，數(shù)值模擬一維速度空間超高斯分布電子經(jīng)自碰撞趨于麥克斯韋分布的過(guò)程，證明有限體積法能高效求解Fokker-Planck方程。進(jìn)一步模擬了HL-2A裝置中性束注入加熱的情況，得到了離子分布函數(shù)隨中性束注入出現(xiàn)非麥克斯韋化。計(jì)算了離子溫度的演化曲線，與實(shí)驗(yàn)結(jié)果符合較好。討論了中性束的功率和能量的影響，結(jié)果表明，隨著束功率和能量的升高，離子溫度升高，相較而言，隨束能量的增大離子溫度升高的幅值較大。由于數(shù)值模擬中性束注入加熱過(guò)程中采用了一維速度空間近似，所以現(xiàn)有計(jì)算程序及模擬結(jié)果存在不足。在中性束注入加熱過(guò)程中，需考慮中性束注入角度對(duì)加熱效率、功率沉積等的影響，以便研究更多的物理問(wèn)題，這也是下一步工作要考慮的問(wè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] 嚴(yán)龍文，段旭如，丁玄同，等. HL-2A裝置物理實(shí)驗(yàn)研究的進(jìn)展[C]∥第十六屆全國(guó)等離子體科學(xué)技術(shù)會(huì)議. 上海：復(fù)旦大學(xué)，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]. 計(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]王建國(guó)，李有宜，李建剛. HT-6M中性束注入加熱理論分析[J]. 物理學(xué)報(bào)，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

國(guó)際熱核聚變實(shí)驗(yàn)堆計(jì)劃專項(xiàng)資助項(xiàng)目(2014GB108002)；國(guó)家自然科學(xué)基金資助項(xiàng)目(11375085)；湖南省核聚變與等離子體物理創(chuàng)新團(tuán)隊(duì)建設(shè)資助項(xiàng)目(NHXTD03)；衡陽(yáng)市科技局資助項(xiàng)目(2013KJ24)；南華大學(xué)研究生科研創(chuàng)新資助項(xià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