羅一鳴 王占輝? 陳佳樂 吳雪科 付彩龍 何小雪 劉亮楊曾辰 李永高 高金明 杜華榮 昆侖集成模擬設計組

1) (核工業(yè)西南物理研究院,聚變科學技術所,成都 610041)

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

3) (西南交通大學物理科學與技術學院,成都 611756)

托卡馬克等離子體物理過程時空尺度跨度大,不同空間區(qū)域(如芯部、臺基區(qū)、刮削層、靶板區(qū))的主要物理過程不同,因此需要采用系統(tǒng)集成方法開展全域多時空尺度物理問題分析.為了更加深入地研究托卡馬克等離子體放電實驗的穩(wěn)態(tài)運行及爬升期間的輸運與約束過程,通常采用多種物理程序開展集成模擬研究,對放電實驗結果進行集成模擬對照,相互驗證并進一步開展物理分析.本文基于OMFIT 平臺,結合HL-2A裝置第37012 炮高比壓放電實驗結果完成了集成模擬驗證與分析,驗證了程序的可靠性與適用性.在該流程中,通過選取適當?shù)哪Ｐ?對實驗參數(shù)進行了校核與補充,經演化后模擬結果與實驗結果比較吻合.在此基礎上,本文進一步采用TGLF 模型開展了芯部靜電漂移波線性不穩(wěn)定性分析,結果顯示NBI 離軸加熱導致H 模約束改善的原因是,該實驗在NBI 功率沉積位置的ETG 不穩(wěn)定性處于被抑制的狀態(tài),輸運由ITG 不穩(wěn)定性占據主導,同時輸運水平降低至新經典水平.

1 引言

目前以托卡馬克為代表的可控磁約束核聚變是核聚變研究的一種主要方式.托卡馬克中等離子體具有復雜的物理現(xiàn)象,其時間尺度、空間尺度的跨度很大[1],且在不同區(qū)域的物理機制也不盡相同,因此無法使用單一理論模型或者物理程序進行準確且快速的模擬.對不同尺度的物理采用不同的模塊分別計算,將平衡、輸運、不穩(wěn)定性、加熱及電流驅動等程序通過適當?shù)牧鞒淘O計耦合起來,即可在盡可能快的時間內,得到盡可能精確的計算結果.

目前較為常見的集成模擬平臺有OMFIT[2],CRONOS[3],METIS[4]和IMAS[5]等,OMFIT 作為目前國際上最為成熟、最為全面的集成模擬平臺,集成了多個等離子物理計算程序,可以涵蓋輸運[6,7]、電流演化[8]、平衡[9]、不穩(wěn)定性分析[10]等物理.目前在其他裝置上,已采用OMFIT 開展了實驗結果的集成模擬分析與驗證[11-15].最早在DⅢ-D 的L 模放電[11]中進行了集成模擬分析與驗證,而后在H 模放電實驗中,實現(xiàn)了芯部與臺基的自洽耦合集成模擬,對實驗結果實現(xiàn)了很好的模擬驗證[12].在EAST 上已經完成了對穩(wěn)態(tài)長脈沖放電[13]、高比壓放電[14]和高自舉電流放電[15]實驗的模擬驗證.然而,在國內HL-2A 托卡馬克裝置上開展的實驗結果的集成模擬分析與驗證工作還比較少.

此外,集成模擬在新型托卡馬克裝置、未來聚變堆的放電運行模式設計中也扮演著極為重要的角色.在完成對裝置的零維設計后,采用集成模擬對裝置的運行模式進行設計,有利于校核設計參數(shù),以及為放電實驗提供重要的參考依據.目前已采用OMFIT 的集成模擬流程,完成了對ITER[16]和CFETR[17-19]等裝置的運行模式設計與優(yōu)化.此外,基于METIS和CRONOS 程序,在DEMO上完成了穩(wěn)態(tài)運行模式的設計工作[20].

本文共分為五個部分.第一部分引言介紹OMFIT集成模擬對實驗結果驗證的研究概況.第二部分介紹本文所采用的OMFIT 集成模擬工作流程及主要程序模塊的物理模型.第三部分介紹HL-2A 裝置實驗放電參數(shù)和初始等離子體溫度和密度剖面,以及集成模擬計算結果.第四部分進一步分析芯部不穩(wěn)定性,探討了芯部約束改善的物理機制.第五部分為本文的總結.

2 OMFIT 集成模擬工作流及物理模型

2.1 OMFIT 集成模擬工作流

本文所采用的集成模擬工作流如圖1 所示.該工作流按照任務目標分成兩部分:芯部實驗結果的集成模擬驗證和芯部不穩(wěn)定性分析.在實驗結果驗證的集成模擬流程中,包含了平衡程序EFIT[21]、電流演化及源項程序ONETWO[22]和剖面演化程序TGYRO[6,7].ONETWO 同時耦合了輔助加熱計算模塊,如計算中性束的NUBEAM[23]、計算低雜波的GENRAY[24]、計算ECRH 的TORAY[25]程序模塊.本工作流程重點在于芯部實驗結果的集成模擬與不穩(wěn)定性分析,所以臺基剖面將固定為實驗剖面.

圖1 OMFIT 芯部等離子體剖面集成模擬流程圖Fig.1.The integrated simulation workflow of core plasma with OMFIT.

考慮到電流擴散時間尺度遠大于輸運時間尺度,為了數(shù)值上的準確性,集成模擬中采用輸運時間尺度作為迭代的時間步長.一個完整的電流演化過程需要非常多的時間迭代步數(shù)才能實現(xiàn),因此在單一程序內同時對電流分布和等離子體動理學剖面如電子密度ne、電子溫度Te、離子溫度Ti等進行高效地演化是不現(xiàn)實的.這里將計算電流演化和等離子體動理學剖面演化的過程分開計算,由ONETWO進行各種源項計算和電流演化,而等離子體動理學剖面的演化則由包含了芯部湍性和新經典輸運的TGYRO 程序進行計算.

將電流和輸運分開計算的基本假設是,在特定條件下的平衡態(tài)物理計算過程中,電流和輸運可以分開計算多次迭代集成達到最終平衡態(tài),等離子體可以在一個遠小于電阻擴散時間尺度的時間內,輸運達到熱平衡狀態(tài).這樣我們可以按照在電流演化的時間尺度上取相應的步長,電流每演化一個步長,直接計算新的電流分布下輸運能夠達到的穩(wěn)定狀態(tài),從而避免了按時間演化輸運所需要的大量計算成本,同時保證了所關心的平衡狀態(tài)演化的計算精度.該流程無需對小尺度輸運過多關注,完成對電流分布和動理學剖面演化后,即可在此基礎上引入新模塊對更為細致的物理問題進行研究.

具體的集成模擬工作流程如下:

1) 從HL-2A 放電平頂段的某一時刻t0出發(fā),將該時刻處理后的實驗數(shù)據導入OMFIT 集成平臺,包括電子密度ne、電子溫度Te、離子溫度Ti、雜質離子密度nimp或有效電荷Zeff、輻射功率Qrad、環(huán)向旋轉速度Ω的剖面分布,以及等離子體電流Ip、環(huán)電壓V、環(huán)向磁場強度Bt和平衡位形等.

2) 保持等離子體動理學剖面不變,在一個很短的輸運時間步長τ內使用ONETWO 計算電流分布,同時由ONETWO 調用輔助加熱計算程序來計算熱源項.由于本文所分析的實驗只有中性束注入(neutral beam injection,NBI)加熱,因此只需調用NUBEAM 程序來計算驅動電流以及快離子分布.在HL-2A 的實驗中,由于能量約束時間通常在10 ms 的量級,為了能夠根據剖面及時更新加熱、電流驅動,因此將時間步長設置為1 ms.

3) 把ONETWO 計算得到的壓強p和電流分布FF′剖面?zhèn)鬟f給EFIT,由EFIT 進行平衡計算,得到新的電流分布下的等離子體平衡.

4) 將2)得到的熱源項和3)得到的等離子體平衡交給剖面演化TGYRO 程序,由TGYRO 計算芯部區(qū)域的等離子體動理學剖面,如電子密度ne、電子溫度Te、離子溫度Ti.

5) 將實驗上的臺基區(qū)等離子體剖面與TGYRO計算得到的芯部等離子體剖面相結合,便得到歸一化小半徑從0到1 上的完整等離子體剖面.

6) 基于更新后的剖面和平衡,重新循環(huán)進行第2)—5)步的操作,在迭代N個循環(huán)后,等離子體剖面計算達到收斂,最終得到一個穩(wěn)態(tài)的解.

2.2 物理模型

EFIT 是托卡馬克平衡重建最常用的程序,其原理是計算磁面坐標下的Grad-Shafranov 方程:

其中p是壓強;F是電流通量;ψ是約化后的極向磁通.在給定初始值p(ψ)和F(ψ) 后,即可由EFIT程序計算得到二維的平衡位形.

ONETWO和TGYRO 都是計算輸運的程序,但原理上不盡相同.ONETWO 采用GLF23[26]和一些簡化的新經典輸運模型分別計算湍流輸運系數(shù)和新經典輸運系數(shù),TGYRO 則是采用TGLF和NEO 程序來分別計算湍流輸運和新經典輸運.在得到輸運系數(shù)后,ONETWO 求解輸運方程來獲得隨時間演化的溫度、密度、電流剖面,TGYRO[6,7]則是通過通量匹配的方法來得到演化終態(tài)的溫度、密度、環(huán)向旋轉剖面.如2.1 節(jié)所述,采用TGYRO來模擬溫度、密度剖面,用ONETWO 來模擬電流剖面.下面詳細介紹TGYRO 相關的物理模型與求解方法.

TGLF 既可以作為線性程序[10]計算不穩(wěn)定性波模的本征函數(shù)、頻率和增長率,也可以作為準線性輸運模型[27]計算湍性輸運通量,包括粒子、能量和環(huán)向動量通量等.在集成模擬工作流程中,TGLF以準線性模型得到湍性輸運通量,傳遞給TGYRO進行剖面演化計算.目前TGLF 湍性通量計算僅適用于捕獲離子模(trapped ion mode,TIM)、捕獲電子模(trapped electron mode,TEM)、離子溫度梯度模(ion temperature gradient,ITG)、電子溫度梯度模(electron temperature gradient,ETG)的漂移波湍流,而無法計算由動理學氣球模(kinetic ballooning mode,KBM)、環(huán)向阿爾芬本征模(toroidal Alfvén eigenmode,TAE)以及高能量粒子模(energetic particle mode,EPM)等引起的湍性輸運.TGLF 的飽和準則通過與GYRO 或CGYRO非線性動理學模擬結果校核而確定的,發(fā)展出SAT0[27]和SAT1[28]兩種不同的飽和準則.通常認為SAT0更適用于低q95的情形(如q95≤6),后者適用于高q95的情形[29].本文研究的實驗q95～4,因此選擇使用SAT0.

新經典物理程序NEO[30-32]是直接求解漂移動理學方程的第一性原理程序,基于線性Fokker-Planck 碰撞算子,通過直接求解一系列由DKE 推導出來的方程來得到等離子體的分布方程,進而得到相應的輸運系數(shù).它可以比較精確地求解新經典效應主導的任何物理量,比如新經典輸運通量、自舉電流、等離子體旋轉等.在計算新經典輸運通量方面,NEO 與流體程序NCLASS和解析公式Chang-Hinton 都有很好的校核,但NEO 比后兩者的適用范圍都廣.

TGYRO 可以快速地將輸運計算到穩(wěn)態(tài),忽略中間的演化過程.TGYRO 通過調用TGLF和NEO分別計算湍性輸運通量和新經典輸運通量,實現(xiàn)總通量Qσ(σi,e)與源項計算得來的目標通量相匹配[32].

給定rr*處的邊界條件Tσ(r*)后,即可積分反演得到整個剖面:

3 實驗結果驗證

3.1 初始實驗等離子體溫度和密度剖面

HL-2A 托卡馬克是一個具有上、下封閉式偏濾器的中型聚變研究裝置,其大半徑為R1.65 m,小半徑為a0.4 m .通常情況下,裝置運行在下單零偏濾器位形.典型的等離子體放電參數(shù)為:等離子體電流IP180—200 kA,線平均電子密度ne(1.5-4)×1019m-3,環(huán)向磁場Bt1.2—1.6 T .此外,在HL-2A 裝置上具有1.5 MW 的中性束注入加熱系統(tǒng)(NBI),5 MW 的電子回旋加熱系統(tǒng)(ECRH)和1 MW 的低雜波電流驅動系統(tǒng)(LHCD).基于HL-2A 的第37012 炮高比壓高約束模放電進行模擬驗證,加熱方式只有NBI 加熱,無射頻波加熱,診斷數(shù)據較為齊全,其放電參數(shù)如圖2 所示.根據放電參數(shù),選取穩(wěn)態(tài)運行階段的某一時刻作為模擬的初始時刻t0,這里t0選取的是進入H 模后的1020 ms,該時刻具體的參數(shù)見表1.第37012 炮是常規(guī)偏濾器下單零位形放電,采用兩束NBI 進行加熱,等離子體參數(shù)相對較高.

表1 第37012 炮在1020 ms 時的參數(shù)Table 1.The parameters of the shot #37012 at the 1020 ms.

圖2 第37012 炮放電參數(shù)(a) 等離子體電流 Ip ;(b) 等離子體儲能 WE;(c) 歸一化比壓 βN和極向比壓 βp ;(d) 線平均電子密度 ;(e) NBI 加熱功率;(f)DαFig.2.The dischargement parameters of the shot #37012:(a) Plasma current Ip;(b) stored energy WE ;(c) normal-izedbeta βNand poloidal beta βp;(d)line-averaged elec-tron density;(e) NBI heating power PNBI;(f)Dα.

實驗診斷獲得的動理學剖面主要有電子密度ne、電子溫度Te、離子溫度Ti、離子旋轉速度ω.電子密度由FIR和微波反射兩種診斷方式進行測量,電子溫度由ECE 診斷進行測量,離子溫度和旋轉速度則由CXRS 系統(tǒng)給出.CXRS 測量 C6+的光譜,通過假定氘離子與碳離子的溫度與旋轉速度相等來獲取氘離子的信息.由于芯部電子密度較高,ECE診斷中出現(xiàn)了截止,因此芯部電子溫度也將結合離子溫度以及等離子體總儲能來近似構建.對第37012 炮1020 ms 時刻的離子溫度和旋轉速度數(shù)據處理如圖3 所示.

圖3 第37012 炮在1020 ms 時的(a)離子溫度和(b)旋轉速度剖面Fig.3.The ion temperature profile (a) and rotation profile(b) of the shot #37012 at the 1020 ms.

電子密度有兩種診斷方式,在邊界通過微波反射進行測量,精度較高,可以在排除異常數(shù)據點后直接使用,而靠近磁軸的電子密度是通過激光干涉進行測量,數(shù)據點較少,可以提供弦積分密度信息.對電子密度的處理則是以微波反射測量值為準,同時參考激光干涉弦積分后的芯部剖面,在實驗誤差允許的范圍內,二次重建一個同時滿足反磁測量得到的線平均密度和微波反射診斷數(shù)據的電子密度剖面,如圖4 所示.

圖4 第37012 炮在1020 ms 時的電子密度剖面處理Fig.4.The treatment of electron density profile of the shot#37012 at the 1020 ms.

在診斷數(shù)據不夠齊全且不自洽時,采納精度高、數(shù)據全的診斷數(shù)據,通過構造的方式得到與宏觀參量(線平均密度、儲能、βN)相匹配的且診斷未給出的電子溫度剖面.本文參考離子溫度剖面形狀,經過計算,最終得到與其他測量結果自洽的電子溫度剖面,如圖5 所示.在3.2 節(jié)中,使用圖3 擬合后的離子溫度和旋轉實驗剖面、圖4 二次優(yōu)化后的電子密度實驗剖面,以及圖5 構造的電子溫度剖面,一起作為OMFIT 集成模擬驗證的初始輸入.

圖5 第37012 炮在1020 ms 時得到的電子溫度剖面Fig.5.The profile of electron temperature of the shot#37012 at the 1020 ms.

3.2 集成模擬與實驗結果對照與分析

按照2.1 節(jié)中介紹的計算流程,在HL-2A 裝置上對第37012 炮進行了模擬驗證.在ONETWO計算中,采用耦合的NUBEAM 程序對中性束粒子進行追蹤計算.模擬時間段是1020—1030 ms,在這個時間段內等離子體處于平頂段穩(wěn)定運行階段,各個等離子體參數(shù)幾乎沒有變化.

由于臺基區(qū)物理的復雜性以及診斷上的誤差,目前對臺基區(qū)域的計算和對照不是很可靠.在對第37012 炮的模擬中,芯部 (ρ＜0.8) 的湍流輸運和新經典輸運采用TGYRO 進行計算,而在邊界處 (0.8＜ρ＜1) 固定剖面,與實驗保持一致.在第37012 炮的實驗診斷里,缺少對雜質密度或有效電荷的診斷,這里按照其他實驗炮的診斷,將有效電荷設置為全空間恒定的常數(shù)2.6.

經過多輪迭代后,各項宏觀參數(shù)及動理學剖面演化趨于收斂.從圖6—圖8 中可以看出,溫度、密度以及壓強略有下降,磁軸處溫度最大誤差也控制在10%以內,但整體上與實驗上吻合得很好,可以證明以上的集成模擬流程可信.對各種成分的電流剖面計算結果如圖7 所示,可以得到自舉電流的份額為52%,NBI 驅動電流份額在30%左右,歐姆電流為18%.

圖6 集成模擬計算中各物理量的多次迭代收斂性 (a1),(b1),(c1) 迭代前后對比;(a2),(b2),(c2) TGYRO 計算點的收斂過程Fig.6.The astringency of each physical quantity in the integrated simulation:(a1),(b1),(c1) The comparison between before and after the iteration;(a2),(b2),(c2) the convergence process of the TGYRO calculating points.

圖7 第37012 炮在1020 ms 時刻的各成分電流剖面Fig.7.The current profiles of each composition of the shot#37012 at the 1020 ms.

圖8 第37012 炮在1020 ms 時刻剖面的模擬結果與實驗結果對照 (a) 壓強剖面;(b) 電子密度剖面;(c) 離子溫度剖面;(d) 安全因子 q 剖面Fig.8.The experiment and simulation profiles comparation of the shot #37012 at the 1020 ms:(a) Pressure;(b) electron density;(c) ion temperature;(d) safety factor q .

在第37012 炮中,由兩束同向的NBI 對主等離子體進行加熱,總功率為1.4 MW.經過NUBEAM程序計算,可以從圖9 中看出,本炮放電采取的NBI 離軸注入,能量沉積位置主要在歸一化小半徑為0.2 處附近,且對離子成分的加熱效果更加顯著.

圖9 NBI 能量密度沉積分布Fig.9.The distribution of NBI deposed energy density.

4 芯部輸運分析

在計算達到穩(wěn)定狀態(tài)后,對TGYRO 最終結果進行分析,從圖10 可以看出,TGYRO 計算的總通量與目標通量可以相匹配,滿足(2)式—(4)式.電子能量輸運通道是湍性輸運(即TGLF 計算出的湍性通量)占據主導,而離子能量輸運通道上,湍性和新經典輸運(即TGLF 計算出的新經典通量)所占份額幾乎相當,在ρ0.2-0.3 附近的ITB(internal transport barrier)區(qū)域氘離子能量輸運由新經典輸運占主導,湍性輸運得到抑制.

圖10 第37012 炮放電在1020 ms 時刻模擬后得到 (a) 離子能量通量;(b) 電子能量通量Fig.10.The ion energy flux (a) and electron energy flux(b) of the shot #37012 at the 1020 ms.

使用TGLF-scan 模塊對上述結果進行靜電漂移波不穩(wěn)定性分析,圖10 為0.2—0.9 區(qū)間內兩支最不穩(wěn)定的本征模式掃描結果,橫坐標為歸一化小半徑,縱坐標為歸一化波數(shù)kθρs,顏色深淺代表不穩(wěn)定性增長率的高低,而紅色代表的是電子抗磁漂移方向 (ω ＞0),藍色代表的是離子抗磁漂移方向(ω＜0).通過圖11 的波數(shù)分析進行模式識別,在0.2—0.3區(qū)域內,高k的不穩(wěn)定性受到抑制,主要是ITG 模主導的湍流,在0.7 以外,有TEM和ITG 得到增長.

圖11 0.2—0.8 區(qū)域內兩支最不穩(wěn)定的本征模式的頻譜Fig.11.The spectrum of two most unstable eigenmode in the 0.2—0.8 region.

在經過TGLF 對不同位置(ρ0.3,0.5,0.8)的微觀不穩(wěn)定進行計算后,給出相應的增長率與波數(shù)的關系,如圖12所示.圖中橫縱坐標分別為歸一化的波數(shù)kθρs和γ(cs/a),其中cscs/Ωs,ΩseB/mic,a為裝置小半徑.按照線性分析的結果來看,高k模式在ρ0.3 附近受到抑制,湍流由紅色小圓點代表的ITG 主導.在偏離ITB 區(qū)域(ρ0.5)內,ETG 逐漸增長,成為湍流的主要因素.而在更靠近邊界處(ρ0.8),TEM和ITG 都得以增長,但主導輸運的仍是ETG 不穩(wěn)定性.

圖12 ρ=0.3,0.5,0.8 處線性不穩(wěn)定性的增長率與波數(shù)的關系(藍色為電子抗磁漂移方向,紅色為離子抗磁漂移方向)Fig.12.The relationship between the growth-rate and wavenumber of the linear instabilities in the ρ=0.3,0.5,0.8(the blue points represent the electron diamagnetic drift direction and the red points represent the ion diamagnetic drift direction).

5 結論

在OMFIT 平臺上,通過集成EFIT,ONETWO/NUBEAM,TGYRO/TGLF/NEO 得到了可以實現(xiàn)對HL-2A 裝置第37012 炮高比壓放電實驗結果的集成模擬驗證流程.在該流程中,驗證了HL-2A 裝置的程序適用性和可靠性,最終芯部集成模擬結果整體上可以很好地與實驗診斷數(shù)據相吻合,磁軸處溫度最大誤差仍然控制在10%以內.在第37012 炮的模擬結果基礎上,進一步采用TGLF 模型進行線性不穩(wěn)定性分析,結果顯示NBI 沉積位置附近的高k微觀ETG 不穩(wěn)定性受到抑制,輸運由低k的ITG 占據主導,輸運水平降低至接近新經典的水平,從而實現(xiàn)了NBI 離軸加熱導致的H模約束改善.

目前的結果尚未考慮臺基和邊界區(qū)物理,無法對邊界雜質輸運影響進行更多細致分析.同時,這里需要更為全面詳細的診斷數(shù)據作為支撐,未來將在本文工作基礎上進一步完善對HL-2A 實驗的集成模擬驗證流程.

感謝GA 方面提供的OMFIT 平臺支持和GACODE,感謝Orso Meneghini、毛瑞、李航、吳木泉對本文工作的指導、建議和幫助.