負三角形變托卡馬克位形下高能量離子激發(fā)魚骨模的模擬研究*

任珍珍 申偉

1) (安徽大學物理與光電工程學院,合肥 230601)

2) (中國科學院等離子體物理研究所,合肥 230031)

負三角形變位形下的托卡馬克具有更低的湍流輸運和更好的能量約束,被認為是未來聚變堆一個更好的選擇.為了探索負三角形變位形下高能量粒子激發(fā)不穩(wěn)定性的特征,使用動理學-磁流體混合模型程序M3D-K 開展了此位形下高能量離子激發(fā)魚骨模的線性不穩(wěn)定性和非線性演化的模擬研究.基于類EAST 參數條件,模擬發(fā)現負三角形變解穩(wěn)理想內扭曲模不穩(wěn)定性,但會致穩(wěn)魚骨模不穩(wěn)定性.非線性模擬發(fā)現在沒有磁流體非線性效應時,負三角形變位形下的魚骨模更不容易飽和,可能的解釋是相比于正三角形變位形,在負三角形變位形下的高能量離子軌道更接近與芯部,因而更容易驅動魚骨模不穩(wěn)定性.這些結果表明考慮高能量粒子激發(fā)的魚骨模不穩(wěn)定性后,負三角形變位形相比于正三角形變位形并沒有明顯優(yōu)勢.

1 引言

托卡馬克是最有前景的磁約束聚變裝置之一,其中等離子體截面的形變對于其約束與穩(wěn)定性具有顯著的影響.因此,對于未來的聚變裝置,形變參數是關鍵的設計參數之一[1?4].近年在TCV 托卡馬克裝置上的研究發(fā)現相比于傳統(tǒng)的正三角形變,負三角形變可降低湍流進而改善等離子體約束[5,6].隨后,DIII-D 托卡馬克裝置上也采用了負三角形變位形在L 模下放電,發(fā)現達到了較高的歸一化比壓值(βN=2.7)以及相當于H 模的約束性能(H98y2=1.2),并且沒有邊界局域模發(fā)生[7].另外,負三角形變位形可以降低偏濾器的熱負載,這些優(yōu)點使其被認為是未來聚變堆的一個更好的選擇[8].

高能量粒子物理對于未來聚變堆中的燃燒等離子體來說是一個很重要的研究方向[9].這是由于在未來的磁約束核聚變堆中,實現自持燃燒的等離子體中的氘氚聚變反應會產生大量的快阿爾法粒子.此外,輔助加熱下如中性束注入、離子回旋加熱等也會產生大量的快粒子.這些快粒子可以激發(fā)出各種不穩(wěn)定性,這些不穩(wěn)定性反過來會引起很強的快粒子輸運并導致快粒子損失到裝置的第一壁,大量的快阿爾法粒子輸運會嚴重降低其對背景等離子體的加熱效率.另外,快粒子可能對磁流體不穩(wěn)定性如新經典撕裂模、內扭曲模、電阻壁模等有顯著的影響,而這些磁流體不穩(wěn)定性有可能嚴重降低等離子體約束性能甚至會導致托卡馬克破裂.

高能量粒子激發(fā)的不穩(wěn)定性包括阿爾芬本征模不穩(wěn)定性、快粒子模、魚骨模等.其中,魚骨模不穩(wěn)定性最早在PDX 裝置上有中性束注入的條件下被觀測到[10],是環(huán)向模數和極向模數都為 1 的一種芯部不穩(wěn)定性.魚骨模不穩(wěn)定性在實驗中被觀測到之后,有兩種理論被用來解釋魚骨模不穩(wěn)定性形成的物理機制.1984 年Chen等[11]認為魚骨模是捕獲高能量離子通過進動頻率共振激發(fā)的.另外,1986 年,Coppi 和Porcelli[12]在考慮熱離子逆磁效應后認為魚骨模是一種類磁流體模,其頻率與熱離子逆磁漂移頻率相當.隨后,魚骨模在多個裝置上被觀測到[13?20],因此魚骨模是一種常見的高能量粒子激發(fā)的不穩(wěn)定性.

然而,目前針對負三角形變位形的磁約束等離子體開展的高能量粒子物理相關的工作較少,DIII-D 裝置上開展了負三角形變位形下的高能量粒子激發(fā)阿爾芬本征模不穩(wěn)定性的實驗,發(fā)現相比于正三角形變位形阿爾芬本征模引起的高能量粒子輸運并沒有明顯減弱[21],但是相關的工作沒有理論模擬方面的分析,主要是實驗觀測結果.因此,本工作針對負三角形變位形等離子體,使用動理學-磁流體混合模型程序M3D-K 開展了高能量粒子激發(fā)的魚骨模不穩(wěn)定性的模擬研究,線性模擬結果發(fā)現負三角形變位形解穩(wěn)理想內扭曲模不穩(wěn)定性,但會致穩(wěn)魚骨模不穩(wěn)定性.然而,非線性模擬發(fā)現在沒有磁流體非線性效應時,負三角形變位形下的魚骨模更不容易飽和.這些結果表明考慮高能量粒子激發(fā)的魚骨模不穩(wěn)定性后,負三角形變位形相比于正三角形變位形并沒有明顯優(yōu)勢.

本文分為以下幾個部分,第2 節(jié)簡單介紹了M3D-K 程序所用的物理模型以及模擬所用的主要參數.第3 節(jié)展示了負三角形變下魚骨模不穩(wěn)定性的線性模擬結果.第4 節(jié)展示了負三角形變下魚骨模不穩(wěn)定性的非線性模擬結果.最后,第5 節(jié)給出了本工作的結論.

2 M3D-K 程序模型以及參數設置

2.1 M3D-K 程序模型

本文使用了動理學-磁流體混合模型程序M3D-K[22,23].該程序使用電阻磁流體模型描述熱等離子體,使用漂移動理學或者回旋動理學模型描述高能量粒子.磁流體方程組在極向截面采用有限元方法求解,在大環(huán)方向采用4 階有限差分或者擬譜法.由于M3D-K 代碼采用柱坐標系,在磁軸處不存在奇異點,因而適用于等離子體芯部磁流體不穩(wěn)定性的數值模擬.對于高能離子部分,為了降低粒子噪聲以及模擬所需粒子數量,程序采用了δf方法.M3D-K 程序通過在動量方程中包含高能離子壓強張量,將高能離子與背景磁流體自洽耦合.對于背景等離子體部分,程序還可以包含雙流體效應,以及包括熱離子逆磁漂移、電子壓強等物理效應.M3D-K 程序已經成功被用來模擬各種磁流體不穩(wěn)定性與高能離子的相互作用[24?35],是模擬高能離子不穩(wěn)定性的有力工具.

M3D-K 程序所用的方程組如下,高能量粒子效應通過壓強耦合進入動量方程:

連續(xù)性方程:

能量方程:

法拉第定律:

安培定律:

歐姆定律:

這里ρ,P,v,E,B,J,γ,μ0,κ,η,ν分別表示質量密度、壓強、流體速度、電場、磁感應強度、電流密度、絕熱系數、真空磁導率、熱傳導系數張量、電阻率、黏滯系數.d/dt=?/?t+v·?.

高能量粒子壓強張量在忽略非對角項的情況下使用Chew-Goldberger-Low (CGL)形式:

式中,I表示單位矩陣;b是磁場方向單位矢量;平行和垂直磁場方向的壓強利用高能量粒子分布函數F在回旋中心坐標 (X,v//,μ) 中計算,

其中ρh=v⊥×b/Ω是回旋半徑矢量,Ω ≡eB/m是回旋頻率,e是電子電荷量,(x,v) 是高能量粒子相空間坐標,μ是磁矩,θ是回旋角.

高能量粒子回旋中心分布函數F可以表示為

其中δ是狄拉克函數.

高能量粒子的運動用回旋動理學或者漂移動理學來描述:

其中,E是總的電場,B是總的磁場,b0是沿平衡磁場方向的單位矢量,B0和δB分別是平衡和擾動的磁場.〈·〉 表示回旋平均,如果是漂移動理學的話直接用粒子導心處的電磁場.變量B＊和B＊＊分別表示為

2.2 初始參數與平衡剖面設置

主要的平衡剖面和參數根據EAST 托卡馬克裝置的#71320 炮設置[36],主要參數如下: 小半徑a=0.45 m,大半徑R0=1.86 m,橢圓形變κ=1.63,三角形 變δ=-0.436,磁軸處磁場B0=1.615 T,磁軸處電子密度ne0=5.52×1019m-3,阿爾芬速度vA=B0/(μ0ρ0)1/2=3.35×106m/s,阿爾芬時間τA=R0/vA=5.55×10-7s,磁軸處包括熱等離子體和快離子的總壓強βtotal,0=4.61%.這里μ0是真空磁導率,ρ0是磁軸處的熱等離子體質量密度.平衡的安全因子剖面和壓強剖面如圖 1 中藍線所示,其中Ψ是歸一化極向磁通,其在磁軸處值為 0,在等離子體最外閉合磁面處為 1.q=1 面的位置如 圖 1 中紅色虛線所 示,在Ψ=0.171 處.在本工作中為了簡化分析,熱等離子體密度的徑向分布取為常數.

圖1 安全因子與總壓強平衡剖面Fig.1.Equilibrium profiles of safety factor and total pressure.

考慮中性束注入產生的高能量離子,其入射能量為E0=58 keV.高能量離子在速度方向以及在實空間的徑向分布表達式為

3 線性模擬結果

首先分析在負三角形變條件下的魚骨模的線性不穩(wěn)定性結果.圖2 對比了正三角形變δ=0.436,無三角形變δ=0,負三角形變δ=-0.436位形下快離子激發(fā)內扭曲模以及魚骨模的情況,圖中Phot,0/Ptotal,0是快離子壓強與總壓強比值,Phot,0是磁軸處快離子壓強值,Ptotal,0是磁軸處總壓強值.可以看出,當沒有快離子效應時,三種位形下的理想內扭曲模都是不穩(wěn)定的,并且負三角形變對理想內扭曲模不穩(wěn)定性起解穩(wěn)作用,正三角形變對理想內扭曲模起致穩(wěn)作用.當快離子壓強比值Phot,0/Ptotal,0增大時,快離子的動理學效應會致穩(wěn)內扭曲模[37?39].然而,當快離子壓強增大到超過一定閾值時,魚骨模不穩(wěn)定性被快離子共振激發(fā).另外,從圖 2 可以看出,無三角形變或者負三角形變下魚骨模很難激發(fā)(激發(fā)閾值值比較高),但一旦被激發(fā)則增長更快;正三角形變下激發(fā)魚骨模的閾值比較低,但被激發(fā)后卻增長得慢一些.對于三角形變影響內扭曲模的分析,之前Eriksson 和Wahlberg[40]的工作發(fā)現正的橢圓形變和正的三角形變對理想內扭曲模起穩(wěn)定作用.Martynov等[41]發(fā)現三角形變對理想內扭曲模的勢能 δW有重要的貢獻,考慮到三角形變后解析結果與KINX 程序結果一致.另外,他們使用KINX 磁流體本征程序模擬發(fā)現一定參數條件下負三角形變對理想內扭曲模起解穩(wěn)作用.這些理論分析結果與本文模擬的結果一致.

圖2 不同三角形變下模頻率和線性增長率與快離子壓強比值 Phot,0/Ptotal,0 的關系Fig.2.Mode frequency and linear growth rate as a function of the fast ion pressure fraction Phot,0/Ptotal,0.

圖3 給出了不同高能量離子壓強比值的模結構.在負三角形變位形下,當快離子壓強比值Phot,0/Ptotal,0=0 時,理想內扭曲模頻率為 0,且模結構是上下對稱的(圖 3(a));當快離子壓強比值Phot,0/Ptotal,0=0.15時,模結構變得扭曲并且有一定大小的模頻率(圖 3(b));當快離子壓強比值Phot,0/Ptotal,0=0.4時,模式變?yōu)轸~骨模不穩(wěn)定性,并且有更高的模頻率(圖 3(c)).作為對比,圖 3(d)給出了正三角形變位形下的魚骨模,可以看出模結構與負三角形變位形下的魚骨模很相似.

圖3 不同高能量離子壓強比值 Phot,0/Ptotal,0 下的流函數U(a) δ=-0.436,Phot,0/Ptotal,0=0 ;(b) δ=-0.436,Phot,0/Ptotal,0=0.15 ;(c) δ=-0.436,Phot,0/Ptotal,0=0.4 ;(d) δ=0.436,Phot,0/Ptotal,0=0.4Fig.3.Velocity stream function U at different fast ion pressure fraction Phot,0/Ptotal,0 : (a) δ=-0.436,Phot,0/Ptotal,0=0 ;(b) δ=-0.436,Phot,0/Ptotal,0=0.15 ;(c) δ=-0.436,Phot,0/Ptotal,0=0.4 ;(d) δ=0.436,Phot,0/Ptotal,0=0.4.

圖4 芯部壓強剖面平坦下模頻率和線性增長率與快離子壓強比值 Phot,0/Ptotal,0 的關系Fig.4.Mode frequency and linear growth rate as a function of the fast ion pressure fraction Phot,0/Ptotal,0 with flat pressure profile.

另外,由于未來聚變堆運行時需要先穩(wěn)定磁流體模式,選取了沒有高能量粒子影響下理想內扭曲模穩(wěn)定的算例開展研究.考慮到環(huán)效應對內扭曲模穩(wěn)定性的影響時,當q=1 面之內的壓強梯度足夠小時,內扭曲模會趨于穩(wěn)定[42].因此,選取了一個芯部壓強剖面平坦并且磁軸處壓強值保持不變的算例,重新計算了不同快粒子驅動時3 個三角形變參數下的內扭曲模的線性穩(wěn)定性,壓強剖面如圖 1的紅線所示.新壓強剖面下的魚骨模增長率以及模頻率如圖 4 所示.可以看出,在新壓強剖面下魚骨模的激發(fā)閾值更高,并且正三角形變下魚骨模的激發(fā)閾值仍然低于負三角形變.

4 非線性模擬結果

本部分考慮沒有磁流體非線性效應下的魚骨模非線性演化,這里將熱等離子體的磁流體響應通過只保留n=1 環(huán)向擾動而限制為線性的.由于線性部分 中Phot,0/Ptotal,0=0.4 條件下負三角形變位形的魚骨模比正三角形變位形的魚骨模增長率高.為了更好地對比,選取了一個更低總比壓值βtotal,0=3.91%的負三角形變位形下的魚骨模算例,并且快離子壓強比值仍保持為Phot,0/Ptotal,0=0.4,如圖 5 所示.可以看出,圖5(a) 中正三角形變魚骨模與負三角形變魚骨模的線性增長率基本一樣,增長率均為γτA=0.00906.但是正三角形變的魚骨模可以非線性飽和,而負三角形變的魚骨模不能飽和.圖5(b) 給出了兩種三角形變位形下魚骨模向下掃頻的特征,正三角形變位形下魚骨模頻率從ω=0.0413ωA向下掃頻到ω=0.00707ωA,負三角形變位形下魚骨模頻率從ω=0.0336ωA向下掃頻到 0.圖6為Phot,0/Ptotal,0=0.4,βtotal,0=4.61%,δ=-0.436 算例的動能和磁能的n=1 分量隨時間演化.可以看出,它們都隨時間線性增長隨后飽和,線性增長率一致.但是,由于本文的非線性模擬是考慮沒有磁流體非線性效應下的演化,在演化中磁能的n=0 分量被限制為不變而只演化n=1 分量.

圖5 沒有磁流體非線性效應的魚骨模非線性演化 (a) 動能的 n=1 分量的 演化;(b) 模頻率演化Fig.5.Time evolution of the fishbone without fluid nonlinearity: (a) n=1 kinetic energy;(b) mode frequency.

圖6 動能和磁能的 n=1 分量隨時間演化Fig.6.Time evolutions of n=1 kinetic energy and magnetic energy.

為了進一步確認三角形變對于魚骨模非線性演化的影響,分析了不同三角形變參數下魚骨模的非線性演化特征.圖7(a)給出了在Phot,0/Ptotal,0=0.4 以及βtotal,0=4.61% 條件下三角形變參數 從δ=0 增大到δ=0.436 時魚骨模的非線性演化.可以看出,當沒有三角形變時,魚骨模非線性并不能飽和;當三角形變從δ=0 增大到δ=0.3 時,魚骨模的線性增長率逐步降低,但是魚骨模在非線性階段仍然不能飽和;當三角形變增大到δ=0.436時,魚骨模的增長率進一步減小,其非線性演化到飽和狀態(tài).作為對比,圖7(b)給出了Phot,0/Ptotal,0=0.4 以及βtotal,0=3.91% 條件下三角形變參數從δ=0 減小到δ=-0.436 時魚骨模的非線性演化.當三角形變δ=0 時,魚骨模的非線性演化到飽和狀態(tài),與圖7(a)中δ=0 的算例對比說明總比壓值降低有利于魚骨模非線性演化趨于飽和.當三角形變從δ=0 減小到δ=-0.2 時,魚骨模的非線性演化仍然是飽和狀態(tài),當三角形變進一步減小至δ=-0.3以下時,魚骨模的非線性演化趨于無法飽和.值得注意的是,三角形變參數從δ=-0.2 降低到δ=-0.436 時,魚骨模的線性增長率增幅不大,但是非線性演化卻變化非常明顯.

圖7 不同三角形變參數下的魚骨模動能的 n=1 分量的演化(a) βtotal,0=4.61% ;(b)βtotal,0=3.91%Fig.7.The n=1 kinetic energy evolution of the fishbone with different triangularity: (a) βtotal,0=4.61% ;(b)βtotal,0=3.91%.

為了解釋負三角形變參數下魚骨模非線性演化不容易飽和的原因,分別在圖 5 的兩個算例中取了兩組快離子對比,每組快離子的能量、拋射角以及所在的平均徑向位置都基本相同,如圖8所示.圖中紅色的粒子軌道對應粒子能量E=57.7 keV,Λ=1.05,Pφ=-0.27,其中Pφ是快離子的環(huán)向角動量,對應于快離子在徑向的位置,Pφ=ZeΨ+mv//RBφ/B.另外,黑色的粒子軌道對應粒子能量E=49.8 keV,Λ=0.985,Pφ=-0.33.可以看出,兩組高能量離子都是捕獲粒子且位于等離子體芯部,與魚骨模結構的位置基本相同.研究發(fā)現,在負三角形變位形下的高能量離子軌道更接近芯部,因此更有利于快離子與魚骨模之間的能量交換,這可能是在負三角形變位形下魚骨模更不容易飽和的原因.

圖8 不同三角形變位形下的捕獲高能量離子軌道 (a)δ=-0.436,βtotal,0=3.91% ;(b) δ=0.436,βtotal,0=4.61%Fig.8.Orbits of trapped fast ions with different triangularity: (a) δ=-0.436,βtotal,0=3.91% ;(b) δ=0.436,βtotal,0=4.61%.

5 結論

本文使用動理學-磁流體混合模型程序M3D-K研究了負三角形變位形下高能量離子激發(fā)魚骨模的線性特征以及非線性演化.基于類EAST 參數條件,線性模擬結果發(fā)現負三角形變對內扭曲模起解穩(wěn)作用,正三角形變對內扭曲模起致穩(wěn)作用.當快離子壓強增大到超過一定閾值時,不同三角形變位形下的魚骨模不穩(wěn)定性都會被快離子共振激發(fā),并且正三角位形下魚骨模的快離子壓強激發(fā)閾值更低.另外,正三角形變會致穩(wěn)魚骨模不穩(wěn)定性,而負三角形變對魚骨模起輕微的致穩(wěn)作用.非線性模擬研究結果發(fā)現降低總比壓值有利于魚骨模非線性演化趨于飽和.另外,發(fā)現負三角形變下魚骨模非線性演化不容易飽和,其可能的原因在于負三角形變位形下的高能量離子軌道更接近芯部,因此更有利于高能量離子驅動魚骨模時將能量傳輸給魚骨模不穩(wěn)定性.

本研究沒有考慮熱等離子的磁流體非線性,因而在魚骨模的非線性演化中沒有考慮模模耦合以及n=0 帶狀流,這些效應應該有助于魚骨模的非線性演化趨于飽和,這將作為以后工作考慮的物理效應.另外,本文采用的是類EAST 參數條件,未來可以對ITER 裝置下做進一步的負三角形變魚骨模不穩(wěn)定性的理論模擬分析.

感謝浙江大學傅國勇教授的討論.本文的數值計算是在中國科學院等離子體物理研究所的“神馬”大型高性能計算集群上開展的.