CFETR參數(shù)下α粒子慢化過程的數(shù)值模擬*

吳相鳳 王豐 林展宏 陳羅玉 于召客 吳凱邦 王正洶

(大連理工大學(xué)物理學(xué)院,大連 116024)

氘氚聚變產(chǎn)生的高能量α 粒子是維持未來托卡馬克反應(yīng)堆等離子體高溫的主要加熱源,良好的α 粒子約束對于維持穩(wěn)態(tài)燃燒等離子體至關(guān)重要.在持續(xù)發(fā)生聚變反應(yīng)的系統(tǒng)中,α 粒子遠(yuǎn)離熱平衡,呈現(xiàn)非麥克斯韋分布.如果忽略軌道效應(yīng),基于局域庫侖碰撞的假設(shè)可以得到α 粒子的經(jīng)典慢化分布,然而由于α 粒子存在較大的漂移軌道寬度,空間輸運不容忽視,為得到更為準(zhǔn)確的α 粒子分布函數(shù),需要開展相關(guān)的數(shù)值計算.本文使用模擬程序PTC (particle tracer code)在中國聚變工程試驗堆(CFETR)不同的放電模式下,采用粒子軌道跟蹤和蒙特卡羅碰撞方法,對α 粒子慢化過程進(jìn)行了數(shù)值模擬,獲得了更為真實的α 粒子分布函數(shù),并將其與經(jīng)典慢化分布進(jìn)行了對比.結(jié)果顯示分布函數(shù)在中等能量附近和經(jīng)典慢化分布存在較大差異.進(jìn)一步的分析表明,這是由于中等能量下α 粒子的較強(qiáng)的徑向輸運引起的.本文的研究結(jié)果對準(zhǔn)確評估α 粒子加熱背景等離子體的能力具有重要參考價值.

1 引言

氘氚聚變反應(yīng)是目前最有希望實現(xiàn)商業(yè)化應(yīng)用的聚變反應(yīng)之一[1,2].隨著國際熱核聚變實驗堆(ITER)工程的持續(xù)推進(jìn),人們開始關(guān)注在ITER 之后的聚變裝置的設(shè)計和建造問題,以便盡早地實現(xiàn)聚變能源商業(yè)化[3,4],中國聚變工程試驗堆(CFETR)將作為ITER 和聚變電站商用示范堆(DEMO)之間的關(guān)鍵銜接點,以推進(jìn)中國聚變能源發(fā)展,推動聚變能源商業(yè)化進(jìn)程,以及實現(xiàn)可持續(xù)地清潔能源方面發(fā)揮重要作用[5,6].

在CFETR 燃燒等離子體中,因聚變生成的攜帶有 200 MW 聚變功率的α 粒子是加熱等離子體的主要方式之一,然而,這些高能量的α 粒子只有被有效慢化才能起到維持等離子體自持燃燒的作用[7?10].在氘氚等離子體中關(guān)于3.5 MeV 的α 粒子在等離子體中的慢化分布或能量沉積的研究已經(jīng)進(jìn)行了半個多世紀(jì)[11?15].α 粒子慢化過程中,在非均勻磁場中漂移并受到電子和離子的摩擦以及與背景離子的角度散射,基于局域庫侖碰撞假設(shè)求解福克-普朗克方程的穩(wěn)態(tài)解可以得到經(jīng)典慢化分布[16,17],它被廣泛應(yīng)用于包括α 粒子、中性束注入粒子等在內(nèi)的快離子相關(guān)物理研究,如加熱、電流驅(qū)動、快離子與Alfvén 本征模的相互作用及快離子輸運等[18?21].然而這種慢化分布具有一定局限性,在真實情況下,軌道效應(yīng)和空間輸運是不可忽略的,同時,粒子漂移軌道寬度也會影響α 粒子在等離子體中的輸運[20?22].這種簡化模型無法準(zhǔn)確描述α 粒子的行為,有必要采用更豐富的物理模型描述非局域的α 粒子慢化過程[23,24].

為了探索α 粒子的慢化過程,本文回顧了經(jīng)典慢化分布理論模型并采用數(shù)值模擬方法進(jìn)行對比分析.在數(shù)值模擬方面使用粒子軌道追蹤耦合蒙特卡羅模擬程序PTC(particle tracer code)[25]在CFETR 穩(wěn)態(tài)運行模式和混雜運行模式下對α 粒子的慢化過程進(jìn)行了計算,模型考慮了粒子的有限軌道效應(yīng)以及新經(jīng)典輸運,采用蒙特卡羅方法計算α 粒子與背景等離子體之間的庫侖碰撞,得到了等離子體平衡條件下α 粒子的慢化分布函數(shù).同時,統(tǒng)計了α 粒子熱化和損失的比例,計算了α粒子對背景等離子體的加熱功率,并得到了α 粒子的總體粒子約束時間和能量約束時間.本研究通過在不同運行模式下對α 粒子整個慢化過程中的α 粒子的產(chǎn)生、約束、對背景等離子體的加熱以及整個分布演化的數(shù)值模擬分析,還有理論和程序模擬的對比分析,為CFETR 的物理設(shè)計提供了有價值的參考依據(jù).

本文第2 節(jié)描述了經(jīng)典慢化分布理論模型、CFETR 穩(wěn)態(tài)和混雜運行模式下的等離子體參數(shù)、PTC 程序物理模型.第3 節(jié)是兩種運行模式下α 粒子慢化過程的模擬結(jié)果,包括α 粒子的粒子數(shù)變化、對背景等離子體的加熱功率、α 粒子的能量變化、α 粒子的慢化分布函數(shù)等.第4 節(jié)是討論分析,包括兩種運行模式下慢化過程的分析、約束時間的分析、理論與程序模擬的慢化分布函數(shù)的分析.第5 節(jié)是總結(jié).

2 經(jīng)典慢化分布與PTC 模型

2.1 經(jīng)典慢化分布

經(jīng)典慢化分布是指高能粒子與熱背景等離子體碰撞導(dǎo)致高能粒子的經(jīng)典減速形成的速度或能量分布,它可以作為描述α 粒子慢化過程的一種理論模型.經(jīng)典慢化分布的推導(dǎo)基于忽略速度擴(kuò)散項和軌道效應(yīng)的假設(shè),根據(jù)這個假設(shè),可以簡化福克-普朗克方程,得到如下式子[16]:

其中fα為α 粒子速度分布函數(shù),Sα為α 粒子源項.此外,假設(shè)vi?vα ?ve,其中vi和ve為背景離子和背景電子的熱速度,vα為典型的α 粒子速度,并假設(shè)背景電子溫度與背景離子溫度相等,即Ti=Te,有高能α 粒子的慢化碰撞算子:

其中聚變反應(yīng)速率S0=nDnT〈σv〉,nD為氘的粒子數(shù)密度,nT為氚的粒子數(shù)密度,σ為反應(yīng)截面,v為相對速度大小.對于能量為3.5 MeV 的α 粒子,初始速度vb=1.3×107m/s.將(2)式和(3)式代入(1)式并進(jìn)行計算,得到穩(wěn)態(tài)時三維的α 粒子速度慢化分布函數(shù):

其中H(vb-v) 為階躍函數(shù),將(4)式乘以 4πv2,得到一維的α 粒子速度慢化分布:

根據(jù)能量E與速度v之間的變換關(guān)系以及階躍函數(shù)的定義,得到一維的α粒子能量慢化分布:

其中,Eb為α 粒子初始能量,Ec為臨界能量,根據(jù)(6)式,可得到CFETR 穩(wěn)態(tài)運行模式下的局域穩(wěn)態(tài)能量慢化分布如圖1 所示.

圖1 電子溫度分別為27.78,14.4 和6.7 keV,對應(yīng)電子密度分別為1.14×1020,9.34×1019 和7.47×1019 m?3 參數(shù)下得到的經(jīng)典能量慢化分布f1, f2, f3Fig.1.Classical energy slowing-down distributions f1,f2 and f3 obtained for the electron temperatures of 27.78,14.4 and 6.7 keV,and their corresponding electron densities of 1.14×1020,9.34×1019 and 7.47×1019 m-3.

接下來,利用各向同性的α 粒子慢化分布計算出α 粒子的密度n和總能量nE,對分布函數(shù)(4)在三維空間積分,得到

根據(jù)粒子平衡方程和能量平衡方程:

接下來的模擬中,將使用PTC 程序?qū)Ζ?粒子的慢化過程進(jìn)行模擬,并與理論模型進(jìn)行比較分析.

2.2 CFETR 參數(shù)與PTC 物理模型

基于CFETR 穩(wěn)態(tài)運行模式和混雜運行模式(v201903)參數(shù)下進(jìn)行了α 粒子慢化過程的模擬,具體等離子體參數(shù)如圖2 所示.由圖2 可看到,在穩(wěn)態(tài)運行參數(shù)下,芯部區(qū)域和邊界區(qū)域等離子體參數(shù)相對比較均勻,混雜運行模式下的溫度略高于穩(wěn)態(tài)模式,密度量級與穩(wěn)態(tài)參數(shù)一致.

圖2 CFETR 中的背景等離子體參數(shù) (a) 穩(wěn)態(tài)運行模式下的密度、溫度和安全因子剖面;(b) 混雜運行模式下的密度、溫度和安全因子剖面Fig.2.Background plasma profiles in CFETR: (a) Density,temperature,and safety factor profiles in steady-state scenario;(b) density,temperature and safety factor profiles in hybrid scenario.

CFETR 穩(wěn)態(tài)運行模式是一種基于等離子體穩(wěn)定性和安全性的非感應(yīng)電流驅(qū)動的運行模式,混雜運行模式是感應(yīng)和非感應(yīng)電流混雜驅(qū)動的脈沖式的運行模式.這兩種運行模式不是獨立的,它們可以相互轉(zhuǎn)換和組合,以實現(xiàn)更優(yōu)秀的聚變反應(yīng)效果和能量輸出.

PTC 程序是一個擁有自主知識產(chǎn)權(quán)的粒子追蹤程序,可以模擬聚變反應(yīng)中產(chǎn)生的α 粒子在等離子體中的慢化過程[25].PTC 程序可以利用全軌道模型和漂移軌道模型來追蹤粒子運動,在α 粒子慢化過程模擬中,由于全軌道計算較為耗時,而且根據(jù)α 粒子的能量和位置估算出α 粒子的回旋半徑在0.5—5 cm,回旋半徑對慢化過程的影響是可以忽略不計的,因此,可以采用漂移軌道的導(dǎo)心運動方程追蹤α 粒子.

在模擬過程中,α 粒子源由單位體積的聚變速率來定義[25],粒子運動方程采用經(jīng)典的四階龍格庫塔方法求解,粒子與背景等離子的碰撞通過蒙特卡羅方法求解.PTC 程序?qū)⒙^程看作一系列彈性散射,在每個彈性散射中,PTC 先根據(jù)α 粒子的初始速度和位置,利用含時演化的軌道方程計算α 粒子的運動軌跡和所處的空間位置;然后,根據(jù)該位置區(qū)域內(nèi)的背景電子和離子的密度、溫度,計算出散射截面;接下來,使用概率分布函數(shù)來模擬碰撞過程中的能量和動量轉(zhuǎn)移.同時,在每個小時間段內(nèi)會隨機(jī)采樣一定數(shù)量的α 粒子進(jìn)行統(tǒng)計,記錄α 粒子的位置、能量等信息,通過大量的隨機(jī)采樣和彈性碰撞模擬,就可以得到α 粒子慢化過程中的分布函數(shù).相較于經(jīng)典慢化分布理論模型,PTC 程序在軌道運動的基礎(chǔ)上,考慮了粒子的新經(jīng)典輸運過程,包括考慮徑向擴(kuò)散和粒子漂移軌道寬度等因素,從而能夠更全面、準(zhǔn)確地模擬α 粒子的慢化過程.

3 模擬結(jié)果

3.1 α 粒子慢化過程

首先,在PTC 程序的輸入模塊中輸入CFETR穩(wěn)態(tài)運行模式下的平衡磁場和熱等離子體剖面、擾動電磁場以及α 粒子初始分布;然后在極向截面劃分三角形網(wǎng)格,將背景等離子體信息插值到網(wǎng)格,載入α 粒子源;接著,使用漂移軌道模型追蹤粒子,使用蒙特卡羅方法求解粒子與背景等離子體的彈性碰撞;在模擬過程中,程序會統(tǒng)計α 粒子數(shù)量、計算每個α 粒子的能量損失率、根據(jù)碰撞截面計算α 粒子對背景的加熱功率.同樣利用CFETR混雜運行模式下的平衡磁場、等離子體剖面等參數(shù)模擬α 粒子的慢化過程,得到兩種模式下慢化過程中的物理量變化對比如圖3 所示.

圖3 CFETR 穩(wěn)態(tài)運行模式(實線)和混雜運行模式(虛線)下的各個物理量隨時間的變化 (a) α 粒子數(shù)量;(b) α 粒子損失率;(c) α 粒子對背景等離子體的加熱功率;(d) α 粒子平均能量Fig.3.Time evolution of various physical quantities in CFETR steady-state scenario (solid lines) and hybrid scenario (dashed lines): (a) Number of α particles;(b) loss rate of α particles;(c) heating power of α particles to the background plasma;(d) average energy of α particles.

由圖3 可以看到,在兩種運行模式下,被統(tǒng)計的在等離子體中的α 粒子數(shù)量、熱化的等離子體數(shù)量、損失的α 粒子數(shù)量變化趨勢幾乎一致,分別在0.68 s 和0.73 s 時,穩(wěn)態(tài)運行模式和混雜運行模式下在等離子體中的α 粒子數(shù)量趨于平穩(wěn),達(dá)到穩(wěn)態(tài)時的α 粒子數(shù)量在穩(wěn)態(tài)運行模式和混雜運行模式下分別為 1.795×1020和1.709×1020,穩(wěn)態(tài)運行模式下的數(shù)量略高于混雜運行模式.穩(wěn)態(tài)運行模式下α 粒子損失率最后維持在 0.41%,混雜運行模式下α 粒子損失率最后維持在 0.38%.穩(wěn)態(tài)運行模式和混雜運行模式下對背景等離子體的加熱功率分別在0.57 s 和0.59 s 時開始保持不變,穩(wěn)態(tài)運行模式下α 粒子對背景電子的加熱功率達(dá)到110 MW,對背景氘離子和氚離子的加熱功率分別為45.7 MW 和31.1 MW,整體略高于混雜運行模式結(jié)果.穩(wěn)態(tài)運行模式和混雜運行模式下α 粒子的平均能量大約經(jīng)過0.64 s 和0.72 s 由初始的3.5 MW 趨于穩(wěn)定的1.61 MW,即α 粒子經(jīng)過慢化過程達(dá)到穩(wěn)態(tài).

此外,也在歸一化極向磁通ψ空間對α 粒子的加熱功率密度進(jìn)行了模擬計算,結(jié)果如圖4 所示,根據(jù)圖4 發(fā)現(xiàn),在兩種運行模式下α 粒子對背景電子和背景離子的加熱功率密度量級一致,在不同的ψ位置有所不同.根據(jù)粒子追蹤結(jié)果也得到了在兩種運行模式下達(dá)到穩(wěn)態(tài)時α 粒子密度分布和在能量空間、極向磁通空間的分布,如圖5 和圖6所示.

圖4 ψ 空間的加熱功率密度Fig.4.Heating power density in the ψ space.

圖5 穩(wěn)態(tài)時α 粒子的密度分布 (a) CFETR 穩(wěn)態(tài)運行模式;(b) CFETR 混雜運行模式Fig.5.The α particle density in steady-state: (a) CFETR steady-state scenario;(b) CFETR Hybrid scenario.

圖6 α 粒子分布函數(shù) (a) 能量空間;(b) 歸一化極向磁通空間Fig.6.The α particle distribution function: (a) Energy space;(b) normalized poloidal magnetic flux space.

3.2 約束時間的計算

根據(jù)模擬結(jié)果,由粒子平衡方程和能量平衡方程可以得到程序模擬得到總體的α 粒子約束時間和能量約束時間,α 粒子源項每秒鐘產(chǎn)生的總粒子數(shù)在穩(wěn)態(tài)運行模式和混雜運行模式下分別為3.5×1020和3.24×1020,穩(wěn)態(tài)時在等離子體中的粒子總數(shù)在穩(wěn)態(tài)運行模式和混雜運行模式下分別為 1.795×1020和1.709×1020,計算得到α 粒子總體粒子約束時間分別為 0.51 s和0.52 s.同樣,根據(jù)α 粒子初始總能量和穩(wěn)態(tài)時總能量計算得到α 粒子的總體能量約束時間在穩(wěn)態(tài)運行模式和混雜運行模式下分別為 0.24 s和0.22 s.

4 討論

4.1 慢化分布對比分析

當(dāng)α 粒子在非均勻磁場中運動時,粒子在軌道上的漂移會導(dǎo)致粒子的輸運方向發(fā)生變化,對于高能量α 粒子,由軌道效應(yīng)引起的輸運是可以忽略的,而對于中等能量下的α 粒子,軌道效應(yīng)引起的α 粒子速度擴(kuò)散較強(qiáng),輸運不可忽略.如圖7 所示,在中等能量E=0.2—2 MeV 附近,α 粒子輸運較強(qiáng),α 粒子在徑向的擴(kuò)散較為明顯,PTC 模擬結(jié)果較為平緩.在高能量尾部(＞3.5 MeV),在模擬中考慮了庫侖碰撞的熱化效應(yīng),存在一定的能量擴(kuò)散[26].在不同的空間位置,由于磁場的不同,軌道效應(yīng)對速度擴(kuò)散的影響也不同,導(dǎo)致在不同位置處的局部分布函數(shù)與理論分布函數(shù)符合程度不一致,在本文中,不同空間位置由不同的歸一化極向磁通(ψ)表示,ψa和ψb為穩(wěn)態(tài)運行參數(shù)下的兩個不同位置,ψc和ψd為混雜運行參數(shù)下的兩個不同位置.從圖7 可以看到,PTC 模擬結(jié)果與經(jīng)典慢化分布在穩(wěn)態(tài)運行模式ψb=0.5—0.6 的符合程度高于ψa=0.1—0.2,混雜運行模式下ψd=0.4—0.5 的符合程度高于ψc=0—0.1.

圖7 PTC 程序得到的能量慢化分布與理論能量慢化分布的對比 (a) 穩(wěn)態(tài)運行模式下 ψa=0.1—0.2和ψb=0.5—0.6 ;(b) 混雜運行模式下 ψc=0—0.1和ψd=0.4—0.5Fig.7.Comparison between the energy slowing-down distribution obtained by PTCcode and the classical energy slowing-down distribution:(a) In steady-state scenario at ψa=0.1 ?0.2 and ψb=0.5 ?0.6 ;(b)in hybrid scenario at ψc=0 ?0.1 and ψd=0.4?0.5.

為了進(jìn)一步驗證分布函數(shù)的不同是由輸運效應(yīng)引起的,分析并對比了Wilkie[27]使用回旋動理學(xué)方法考慮存在徑向輸運得到的修正慢化分布函數(shù).其中,利用高能量下以徑向擴(kuò)散項占主導(dǎo)的擴(kuò)散系數(shù)Drr,得到的動理學(xué)方程為

經(jīng)過修正的三維的穩(wěn)態(tài)慢化分布函數(shù)為

將其乘以 4πv2并根據(jù)能量E與速度v之間的變換關(guān)系以及階躍函數(shù)的定義,得到修正的一維的α 粒子能量慢化分布:

在其他一些研究中,也發(fā)現(xiàn)中等能量下的徑向輸運是不可忽略的.Hauff等[22]的研究發(fā)現(xiàn),經(jīng)典慢化分布中的局域慢化假設(shè)在中等能量下失效,需要考慮軌道效應(yīng)引起的徑向輸運,通過垂直去相關(guān)機(jī)制的研究發(fā)現(xiàn),漂移軌道和拉莫爾軌道導(dǎo)致粒子通量是個關(guān)于能量E的函數(shù),粒子擴(kuò)散系數(shù)隨著能量增大以1/E減小,中等能量下,粒子擴(kuò)散系數(shù)很大,徑向輸運很強(qiáng).Wilkie等[21]也發(fā)現(xiàn)經(jīng)典慢化分布中的通量面局域碰撞導(dǎo)致慢化的假設(shè)不成立,碰撞和輸運相關(guān)的時間尺度會相互競爭,在中等能量下,粒子徑向輸運時間小于碰撞特征時間,因此經(jīng)典慢化分布不再適用.同樣,在JET 的D-T 實驗中[28],也觀測到中等能量下的徑向通量較大,輸運較強(qiáng).因此準(zhǔn)確的α 粒子慢化分布函數(shù),需要考慮徑向擴(kuò)散,更真實的模擬結(jié)果有助于更好地預(yù)測α 粒子的輸運行為和能量沉積分布.

4.2 約束時間的討論分析

根據(jù)經(jīng)典慢化分布,可得到理論的α 粒子的粒子約束時間和能量約束時間的表達(dá)式.選取穩(wěn)態(tài)運行模式下ψ=0.1—0.2 內(nèi)的平均電子溫度和密度,得到理論α 粒子慢化時間為 0.93 s,粒子約束時間和能量約束時間分別為 0.76 s和0.31 s.ψ=0.5—0.6 內(nèi)α 粒子慢化時間為 0.42 s,粒子約束時間和能量約束時間分別為 0.48 s和0.17 s.對比PTC穩(wěn)態(tài)運行模式下的模擬結(jié)果,α 粒子總體慢化時間約為 0.68 s,總體粒子約束時間為 0.51 s,能量約束時間為 0.24 s,理論計算結(jié)果與程序模擬結(jié)果量級一致.

同樣局部選取混雜運行模式下ψ=0.1—0.2 內(nèi)的平均電子溫度和密度,計算得到理論α 粒子慢化時間為 1.01 s,粒子約束時間和能量約束時間分別為 0.79 s和0.33 s,ψ=0.5—0.6 內(nèi)α 粒子理論慢化時間為 0.43 s,粒子約束時間和能量約束時間分別為 0.48 s和0.17 s.穩(wěn)態(tài)運行模式下芯部的電子溫度和密度較為平緩,混雜運行模式下芯部的電子溫度和密度變化較快,使得芯部區(qū)域的約束性能有所不同,此外,穩(wěn)態(tài)運行模式下的聚變功率略大于混雜運行模式,單位時間內(nèi)產(chǎn)生的α 粒子更多,慢化時間也有所不同.

5 結(jié)論

本研究在CFETR 參數(shù)下用PTC 程序模擬研究了聚變產(chǎn)物α 粒子的慢化過程,對比分析了α粒子在慢化過程中不同運行模式下對背景等離子體的加熱功率、粒子數(shù)變化、約束時間、能量變化的區(qū)別,并得到了考慮軌道效應(yīng)和輸運之后更準(zhǔn)確的α 粒子慢化分布函數(shù).經(jīng)過PTC 模擬結(jié)果與經(jīng)典慢化分布的對比發(fā)現(xiàn),在中等能量下二者分布的峰值和能量變化存在不一致,徑向輸運會抹平中等能量下的能量慢化分布,使慢化分布峰值降低.本文進(jìn)一步對比了考慮輸運之后的修正慢化分布,發(fā)現(xiàn)考慮徑向擴(kuò)散之后,修正慢化分布只在中等能量下與慢化分布存在差異,且修正結(jié)果與PTC 模擬結(jié)果符合程度較好.研究表明,經(jīng)典慢化分布在中等能量下具有較大誤差,PTC 程序考慮了軌道效應(yīng)與輸運過程,可以得到更為準(zhǔn)確的α 粒子分布函數(shù),尤其在中等能量下的分布更符合真實的α 粒子的行為.本文的研究結(jié)果對于準(zhǔn)確預(yù)測α 粒子沉積分布和評估α 粒子加熱背景等離子體的能力具有重要參考價值.