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

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

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

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

1 引言

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

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

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

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

2 經典慢化分布與PTC 模型

2.1 經典慢化分布

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

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

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

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

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

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

圖1 電子溫度分別為27.78,14.4 和6.7 keV,對應電子密度分別為1.14×1020,9.34×1019 和7.47×1019 m?3 參數下得到的經典能量慢化分布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,對分布函數(4)在三維空間積分,得到

根據粒子平衡方程和能量平衡方程:

接下來的模擬中,將使用PTC 程序對α 粒子的慢化過程進行模擬,并與理論模型進行比較分析.

2.2 CFETR 參數與PTC 物理模型

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

圖2 CFETR 中的背景等離子體參數 (a) 穩態運行模式下的密度、溫度和安全因子剖面;(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 穩態運行模式是一種基于等離子體穩定性和安全性的非感應電流驅動的運行模式,混雜運行模式是感應和非感應電流混雜驅動的脈沖式的運行模式.這兩種運行模式不是獨立的,它們可以相互轉換和組合,以實現更優秀的聚變反應效果和能量輸出.

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

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

3 模擬結果

3.1 α 粒子慢化過程

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

圖3 CFETR 穩態運行模式(實線)和混雜運行模式(虛線)下的各個物理量隨時間的變化 (a) α 粒子數量;(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 可以看到,在兩種運行模式下,被統計的在等離子體中的α 粒子數量、熱化的等離子體數量、損失的α 粒子數量變化趨勢幾乎一致,分別在0.68 s 和0.73 s 時,穩態運行模式和混雜運行模式下在等離子體中的α 粒子數量趨于平穩,達到穩態時的α 粒子數量在穩態運行模式和混雜運行模式下分別為 1.795×1020和1.709×1020,穩態運行模式下的數量略高于混雜運行模式.穩態運行模式下α 粒子損失率最后維持在 0.41%,混雜運行模式下α 粒子損失率最后維持在 0.38%.穩態運行模式和混雜運行模式下對背景等離子體的加熱功率分別在0.57 s 和0.59 s 時開始保持不變,穩態運行模式下α 粒子對背景電子的加熱功率達到110 MW,對背景氘離子和氚離子的加熱功率分別為45.7 MW 和31.1 MW,整體略高于混雜運行模式結果.穩態運行模式和混雜運行模式下α 粒子的平均能量大約經過0.64 s 和0.72 s 由初始的3.5 MW 趨于穩定的1.61 MW,即α 粒子經過慢化過程達到穩態.

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

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

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

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

3.2 約束時間的計算

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

4 討論

4.1 慢化分布對比分析

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

圖7 PTC 程序得到的能量慢化分布與理論能量慢化分布的對比 (a) 穩態運行模式下 ψ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.

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

經過修正的三維的穩態慢化分布函數為

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

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

4.2 約束時間的討論分析

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

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

5 結論

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