鄒云鵬 陳錫熊 陳偉

1) (核工業(yè)西南物理研究院,成都 610041)

2) (通用原子能,加州 CA 92186-5608)

3) (中國科學(xué)技術(shù)大學(xué)核科學(xué)與技術(shù)學(xué)院,合肥 230026)

基于臨界梯度模型,使用TGLFEP 和EPtran 兩個程序可以模擬分析阿爾芬本征模引起的高能量粒子輸運問題.本文在原有模型的基礎(chǔ)上,加入了兩點改進使模擬結(jié)果更接近實驗.其一,考慮閾值剖面的演化過程.判斷閾值的物理量由密度梯度(dn/dr)改為歸一化的密度梯度((dn/dr)/(n/a)),并且使用TGLFEP 模擬證明閾值與高能量粒子密度成反比例關(guān)系,也就是說,當(dāng)密度降低時,閾值會增大.第二,考慮有限軌道寬度效應(yīng).使用OBRIT 程序計算高能量粒子的損失錐,并輸入到EPtran 程序中,從而增加了一種高能量粒子的損失通道.利用DIII-D#142111 和#153071 進行實驗驗證,結(jié)果表示改進后的模型更接近實驗.最后,利用神經(jīng)網(wǎng)絡(luò)代替TGLFEP 計算臨界梯度,并實現(xiàn)EPtran 的并行計算以縮短運行時間.以此建立一個EP 模塊并加入到OMFIT 集成模擬中,模擬結(jié)果表示當(dāng)阿爾芬本征模驅(qū)動高能量粒子輸運,會導(dǎo)致芯部的壓強和電流降低,從而提升當(dāng)?shù)氐陌踩蜃?改變平衡位形.

1 引言

中國聚變工程實驗堆(CFETR)[1,2]的主要目的之一是實現(xiàn)自持燃燒,因此高能量粒子(EP)分布的預(yù)測成為了一個重要課題.在過去的研究中,通常使用經(jīng)典慢化分布(classical slowing down)[3,4]模型,但是DIII-D 的實驗表明在不穩(wěn)定的阿爾芬本征模(AE)存在時,診斷得到的EP 剖面與經(jīng)典慢化分布預(yù)測的剖面相差甚遠[5,6].類似的現(xiàn)象在其他實驗裝置上[7-9]也同樣被觀測到.在理論上,EP 在實空間的梯度會驅(qū)動AE,而不穩(wěn)定的AE會反過來增強EP 的徑向輸運,從而導(dǎo)致EP 密度剖面變平.當(dāng)考慮多個不穩(wěn)定的AE 時,在實空間會產(chǎn)生島狀結(jié)構(gòu)的重疊,這會導(dǎo)致粒子的隨機擴散[10],當(dāng)然這種現(xiàn)象需要AE 的振幅足夠大時才會發(fā)生,因此本文重點討論的是單模影響EP 分布的情況.此外,AE 還會改變EP 的運動軌跡從而增加軌道損失,當(dāng)AE 驅(qū)動EP 向外輸運時會將粒子推進損失錐中,從而增大EP 的損失份額.因此,在實際托卡馬克中,EP 分布和AE 穩(wěn)定性構(gòu)成了一個復(fù)雜的負反饋系統(tǒng).

為了解釋經(jīng)典慢化模型與實驗結(jié)果的區(qū)別,并使模擬結(jié)果更接近實驗,國際上研發(fā)(或改進)了,例如ORBIT[11],MEGA[12]和kick model[13]等模型.然而自洽的數(shù)值模擬需要大量的計算資源,因此本文選擇一個簡化的模型來預(yù)測EP 輸運和再分布,即臨界梯度模型.該模型的主要思路為AE的不穩(wěn)定性主要由EP 壓強梯度驅(qū)動,那么EP 的壓強梯度應(yīng)該存在一個閾值使AE 達到臨界穩(wěn)定(這個閾值的存在已經(jīng)在理論[14]和實驗中被證實[15]),當(dāng)EP 壓強梯度超過這個閾值時,AE 的振幅快速上升,EP 輸運系數(shù)快速增大(stiff transport),直到EP 壓強梯度降回到閾值.基于這個理論,本文使用TGLFEP 和EPtran 的組合來計算EP 的輸運和再分布,其中TGLFEP 用來計算臨界梯度剖面,而EPtran 用來計算EP 的再分布.并在此基礎(chǔ)上,增加兩個改進點: 第一,優(yōu)化臨界閾值的計算方法.EP 梯度的閾值不再是一個固定值,而是與EP 密度成反比例關(guān)系,當(dāng)EP 密度降低時,閾值會隨之增大;第二,考慮有限軌道寬度效應(yīng)(FOW)造成的EP 損失.舊模型僅考慮EP 在最外層閉合磁面(LCFS)處,由于擴散而造成的損失,本文使用ORBIT 程序計算EP 的運動軌跡并確定其相空間的損失錐.

采用兩次安全因子剖面不同DIII-D 放電實驗來驗證改進后的臨界梯度模型,其中#142111 為反磁剪切[16],#153071 為單調(diào)遞增[17].EP 的初始分布使用TRANSP[18]無輸運的模型計算得到.使用臨界梯度模型模擬環(huán)向阿爾芬本征模(toroidal Alfvén eigenmode,TAE)引起的EP 輸運,并將輸運后的損失份額和沉積剖面與實驗進行對比.最后,利用神經(jīng)網(wǎng)絡(luò)模型替代TGLFEP 計算臨界梯度,并配合EPtran 程序的并行化,縮短計算時間.基于OMFIT 集成模擬框架,搭建一個高能量粒子模塊(module of EP,MOE),從而實現(xiàn)集成模擬中高能量粒子壓強自洽的迭代.

2 臨界梯度模型

2.1 臨界梯度計算

TGLF 程序基于回旋朗道流體模型,主要用來計算湍流物理,但是程序本身采用局域的流管(flux-tube)模型,因此同樣可以用于AE 的計算.由于TGLF 使用準(zhǔn)線性模型,因此其計算速度要遠大于回旋動理學(xué)程序.基于DIII-D 的H 模放電實驗,TGLF 程序已經(jīng)用GYRO 程序線性和非線性模擬校準(zhǔn).TGLF 程序主要求解6 個矩方程,即密度、平行速度、平行壓強、總壓強、平行能量流和總能量流.本征方程采用系列厄密函數(shù)表示.過去研究結(jié)果表示,TGLF 的模擬結(jié)果對基函數(shù)中的高斯寬度θw的設(shè)置非常敏感[19,20],通常選擇對應(yīng)增長率最大的θw.因為程序自動尋找該值的方法僅對湍流適用,并不適用于AE 的計算,所以改為掃描一定區(qū)間內(nèi)的增長率并選擇對應(yīng)最大值的θw.以#142111 為例,掃描高斯寬度0.4—2.0 區(qū)間.另外還有一個原則是模的頻率和增長率在一定區(qū)間內(nèi)有良好的魯棒性.圖1 展示了3 個環(huán)向模數(shù)n=3的模的增長率和頻率隨高斯寬度的變化.因為橙色和紫色曲線均在一定區(qū)間內(nèi)表現(xiàn)出良好的-魯棒性,所以被認為是TAE.而藍色曲線的波動過大,且不滿足輸出邊界條件(outgoing boundary condition)[21],因此被認為是一種數(shù)值噪聲.TGLF EP 是TGLF 的并行的封裝程序,可以自動求解多個磁面的高斯寬度并計算AE 的頻率和增長率.

圖1 (a)增長率和(b)頻率與高斯寬度的關(guān)系,其中橙色和紫色的曲線表示TAE,藍色線被認為是一種數(shù)值噪聲,作為參考,黑色實線為MEGA 的模擬結(jié)果Fig.1.(a) Growth rate and (b) frequency as a function of Gaussian width.The orange and purple curves represent acceptable TAE calculated using TGLFEP.Because the blue curve has a strong perturbation,the mode is treated as polluted by numerical noise.For comparison,the MEGA results are depicted by the black lines.

2.2 高能量粒子輸運

在臨界梯度剖面確定下來之后,就可以使用EPtran 程序來計算EP 的再分布.EP 在(r,E,λ,σ)空間的輸運方程如下:

空間坐標(biāo)r,E,λ,σ 分別為徑向位置、能量、投擲角和運動方向(沿著磁場方向為正,逆著磁場方向為負).方程(1)等號右邊包括慢化項、源項和螺旋角空間的散射項[22].其中fEP表示高能量粒子分布函數(shù),νd表示螺旋角散射率,S0表示EP 源,E0為初始能量(α 粒子為3.5 MeV,中性束則為注入能量),τs為慢化時間,Ec為臨界能量,TEP為等效麥克斯韋溫度,a是小半徑.V '=?rV,V E==v///v,λ=μB/E,這里v表示粒子的速度,v//表示平行磁場方向的速度,μ為磁矩,B為磁感應(yīng)強度.在外中平面,0 ≤ λ ＜ λTP表示通行粒子,λTP≤ λ ≤ 1 表示捕獲粒子,其中λTP(r)=B(0)/B(π)[23].

其中,當(dāng)x＜ 0 時,[x]＞0=0,為局域的臨界梯度;nEP為高能量粒子密度;DITG/TEM為湍流輸運系數(shù);DAE并不是從物理推導(dǎo)的,通常將其設(shè)置得足夠大以驅(qū)動EP 快速的輸運.在過去的研究中,ITER 模擬設(shè)置為0.3[24],DIII-D 模擬設(shè)置為10[25].實際上,再分布剖面對DAE的設(shè)置并不敏感.因此根據(jù)已經(jīng)發(fā)表的文獻 [24,25],在計算DIII-D 的實驗時,設(shè)置DAE=10,在計算CFETR 時,設(shè)置DAE=0.3.

圖2 能量相關(guān)(橙)/不相關(guān)(藍) Angioni 模型對比Fig.2.Comparison between energy dependent (orange) and energy independent (blue) Angioni diffusion coefficients.

根據(jù)AE 增長率表達式[27],方程(4)改寫為

其中,LnEP為EP 密度的特征長度,將在3.1 節(jié)詳細說明此改動.

利用方程(1)計算EP 的分布后,對其積分可以得到EP 的密度和等效麥克斯韋溫度

EP 的初始分布采用經(jīng)典慢化分布:

根據(jù)方程(6)可知,EP 密度剖面為nsd(r)=S0τsI2(vc/v0),溫度為Tsd=(2I4/3I2)E0.其中積分函數(shù)被定義為.對于α 粒子,S0由背景DT 密度和相應(yīng)的反應(yīng)截面計算;對中性束注入來說,S0由NUBEAM 程序[28]計算.

3 臨界梯度模型的優(yōu)化

3.1 閾值剖面演化

在已經(jīng)發(fā)表的研究中,Waltz等[23-25]討論了多種臨界穩(wěn)定的判據(jù),本文使用判據(jù)γAE+ITG/TEM=γITG/TEM來判斷是否達到臨界穩(wěn)定(γAE+ITG/TEM指考慮背景等離子體時AE 的增長率,γITG/TEM指湍流的增長率).在舊模型中,EP密度的特征長度(LnEP)保持不變,不斷降低EP 密度(密度梯度成比例降低),直到達到臨界穩(wěn)定的判據(jù),此時的密度梯度則被視為臨界梯度.在這個假設(shè)下,密度剖面的形狀是不變的,顯然這不符合實際情況.因此使用TGLFEP 計算一個二維矩陣(a/LnEP,nEP)來確定臨界的歸一化梯度與EP 密度的對應(yīng)關(guān)系,利用這個矩陣,可以在EP 密度下降時,使a/LnEP上升.為了說明我們的改進,選擇一個解析的AE 增長率的表達式[27]:

其中,ωi為增長率,k//為平行波矢,vA為阿爾芬速度,n為環(huán)向模數(shù),βc和βα分別為背景粒子和高能量粒子的比壓.G和H是與粒子速度相關(guān)的擬合函數(shù),具體表達式見文獻 [27].等式右邊的第1項為背景離子和電子的朗道阻尼,當(dāng)假設(shè)平衡不變時,其為常數(shù).第2 項包括了EP 對AE 的阻尼和驅(qū)動,其中βα和δα(∝dpα/(pαdr),這里pα為高能量粒子壓強)兩個量會在EP 輸運時改變.由于EP 的驅(qū)動中同時包含壓強和壓強梯度,因此應(yīng)隨著EP 剖面變化而變化.令方程(8)等于0,則可以得到如下形式:

其中k1包含了背景粒子的朗道阻尼,k2包含了EP 本身的阻尼和溫度梯度驅(qū)動,在忽略溫度變化時,這兩個系數(shù)在每個磁面上都是固定的.到此為止,我們證明了與nEP成反比例關(guān)系.圖3(a)中,紅色虛線表示對應(yīng)初始EP 密度的a/LnEP,黑色實線表示經(jīng)典慢化分布的a/LnEP,黑色實線高于紅色虛線的區(qū)域表示AE 是不穩(wěn)定的.針對這個不穩(wěn)定區(qū)域,使用TGLFEP 掃描a/LnEP與nEP,圖3(b)展示了3 個不同位置,TGLFEP 的模擬結(jié)果也驗證了這種反比例關(guān)系.當(dāng)EP 輸運時,芯部的EP 密度是降低的,那么閾值則會相應(yīng)的上升.以ρ=0.4 為例(圖3(c)),紅色的叉表示初始時刻,黑色和紅色箭頭分別表示舊模型和改進后的模型下隨EP 密度的變化軌跡,這可以明顯看到兩個模型的差別.

圖3 (a)輸運前的AE 不穩(wěn)定區(qū)間;(b)在ρ=0.4,0.5,0.6 處,a/LnEP與nEP 的反比例關(guān)系;(c) ρ=0.4 處,舊模型和改進模型模擬的密度演化軌跡Fig.3.(a) Unstable AE region of n=3 without transport.Critical a/LnEP is depicted by dash red curve,and a/LnEP of classical slowing down distribution is depicted by solid black curve;(b) inverse proportional function between critical a/LnEP and nEP at ρ=0.4,0.5,0.6;(c) density evolution trajectories of previous CGM (black) and improved CGM (red) at ρ=0.4.

3.2 有限軌道寬度(FOW)效應(yīng)

磁場的不均勻性導(dǎo)致粒子導(dǎo)心軌道的漂移被稱為有限軌道寬度(FOW)效應(yīng),在等離子體邊界附近,能量越高的粒子就越容易在FOW 效應(yīng)的作用下從LCFS 損失掉,這種現(xiàn)象在實驗中是普遍存在的[16].然而舊模型并沒有考慮這種損失機制,因此只有很少的EP 會損失,而大量的粒子都沉積在等離子體邊界附近.本文使用ORBIT 程序計算EP在相空間中的損失錐,并將其加入到輸運模型中.

在ORBIT 計算中,初始EP 分布選擇在(ρ,E,ξ)空間的均勻分布,且僅考慮外中平面的粒子.由于ORBIT 程序采用二維空間坐標(biāo),因此這里僅考慮外中平面.當(dāng)然,由于粒子更容易在低場側(cè)損失,因此這會導(dǎo)致過度的估計EP 的損失份額.此外,在AE 擾動下,一部分原本可以約束的粒子會轉(zhuǎn)變?yōu)閾p失粒子,即AE 會增大EP 的損失錐.以#142111 為例,圖4(a)和圖4(b)分別展示了捕獲粒子和通行粒子在無擾動(藍色)和有AE 擾動(紅色)下的運動軌跡.其中通行粒子在高場側(cè)的軌跡被放大到圖4(c)中.在每幅圖中,黑色叉表示粒子的初始位置,紅色叉表示粒子從LCFS 逃逸的位置,黑色箭頭表示運動方向.

圖4 (a)捕獲粒子和(b)通行粒子在AE 擾動下的運動軌跡示例;(c)圖(b)在高場側(cè)的放大圖Fig.4.Representative trajectories of (a) trapped and (b) counter passing particles by including AE perturbation;the high magnetic field side of (b) is enlarged in (c) to reveal trajectory variations.

圖5(a)展示了無AE 擾動時不同能量的EP的損失邊界,在該邊界的右側(cè)表示粒子的損失區(qū)域.由于初始粒子均選擇在外中平面上,因此所有的損失粒子都具有負的螺旋角.損失錐的面積隨粒子能量增大而增大.對于能量為注入能量80 keV的粒子,其損失錐的峰靠近磁軸.圖5(b)展示了AE 擾動下?lián)p失錐的進一步增大.根據(jù)文獻 [5,16],AE 的擾動幅度設(shè)為dB/B≈O(10-4).將得到的損失錐輸入到EPtran 程序中,令損失錐內(nèi)的EP 分布函數(shù)為零,從而等效的考慮EP 由于FOW 效應(yīng)的損失.

圖5 (a)無擾動時,不同能量EP 在(ρ,v///v)空間的損失錐;(b)有無AE 擾動下,80 keV 的EP 損失錐對比圖Fig.5.(a) Loss boundary in (ρ,v///v) space with different energies;the loss region is on the right side of the curve.(b) loss cone without/with AE perturbation.The blue area represents the loss cone without AE,and the additional loss by including AE is highlighted by the red area.

4 臨界梯度模型的實驗驗證

4.1 #142111 放電實驗

實驗中利用ECE 測量到525 ms 時刻頻率為70—90 kHz 的低nTAE,這些模在空間上重疊,徑向位置均在ρ ≈ 0.4.MEGA 程序的模擬表示,n=1—5 的TAE 的能量演化如圖6(a)所示,n=3的TAE 增長率最大且最先達到飽和.該TAE 在線性階段((a)中紅色虛線對應(yīng)的位置)的模結(jié)構(gòu)如圖6(b)所示.圖7(a)展示了TGLFEP 程序計算的n=1—5 的模的增長率在徑向上的分布情況,可以看到從ρ ≈ 0.1 到0.6 這個范圍內(nèi),n=1—4 分別占主導(dǎo),因此不同徑向位置的臨界梯度均考慮當(dāng)?shù)卦鲩L率最大的模(圖7(b)).

圖6 MEGA 模擬結(jié)果 (a) n=1—5 的TAEs 的能量演化過程;(b) n=3 的TAE 的模結(jié)構(gòu),vr,cos 為徑向擾動速度的cos 分量,不同的極向諧波用不同的顏色區(qū)分Fig.6.(a)Evolutionof energy with n=1-5 TAEs by MEGA;(b) cosine part of radial velocity for the most unstable n=3 TAE.

圖7 (a)n=1—5 的TAEs在不同磁面的增長率;(b)單n和多n計算的臨界梯度剖面的對比Fig.7.(a) Growt hratein eachfluxsurfaceofn=1-5TAEsbyTGLFEP;(b)comparisonbetweencritical a/LnEP profiles with n=1-4 and n=3.

輸運后的EP 徑向密度剖面如圖8(a)中的紫色曲線所示,為了可以更清晰比較,圖8 中還加入了經(jīng)典慢化分布,即輸運前的剖面(黑色)、實驗結(jié)果(紅色)、舊臨界梯度模型的模擬結(jié)果(藍色)、改進后的臨界梯度模型但不考慮AE 對損失錐影響的模擬結(jié)果(綠色)以及MEGA 的結(jié)果[29].從圖8(a)可以看出,改進后的臨界梯度模型的模擬結(jié)果比舊模型更接近實驗結(jié)果.這里定義EP 的損失份額為損失的粒子數(shù)與經(jīng)典慢化分布預(yù)測的約束的總粒子數(shù)之比,那么舊模型的損失份額為15%,而改進后模型的模擬結(jié)果超過40%,這與實驗的約50%的損失份額非常接近.而且AE 區(qū)域梯度過小和EP 堆積在邊界附近的問題也都有所改善.在臨界梯度模型中,EP 的輸運過程可以理解為: 不穩(wěn)定的AE 不僅增強EP 的徑向輸運,而且會輕微增大損失錐,EP 向外輸運進入損失錐中后,由于FOW 效應(yīng)而直接從LCFS 逃逸出去.最后,由于考慮了負螺旋角的分布,模型可以計算EP 在螺旋角空間的分布(圖8(b)),黑色曲線對應(yīng)60 keV 的EP 的分布,紅色曲線則為0—80 keV全部EP 的分布.

圖8 (a) EP 密度剖面對比圖,經(jīng)典慢化分布(黑色)、實驗(紅色)、舊臨界梯度模型(藍色)、考慮(紫色)/不考慮(綠色)AE 擾動對損失錐影響的改進的臨界梯度模型、MEGA(黃色);(b) EP 在螺旋角空間的分布Fig.8.(a) Density profile comparison: Black curve represents classical slowing down;the red curve is inferred from experiment data;blue curve represents original CGM without loss cone effect;purple/green curve is improved CGM with loss cone from AE perturbed/unperturbed orbits;yellow curve represents MEGA results.(b) EP redistribution in pitch angle space.

4.2 #153071 放電實驗

在#153071 放電實驗中,頻率在100—200 kHz的低nTAE 在靠外側(cè)的位置被激發(fā),但是忽略等離子體旋轉(zhuǎn),模擬得到了一個在ρ=0.2—0.6 范圍內(nèi)的一個較寬的模.增長率最大的n=4 的TAE如圖9 所示,圖中實線和虛線分別為擾動速度的余弦和正弦分量.FOW 效應(yīng)引起的損失錐如圖10所示,由于本次實驗AE 振幅較弱,且AE 的徑向?qū)挾扰c損失錐的重疊較小,因此在計算損失錐時沒有考慮AE 擾動的影響.

圖9 n=4 的TAE 的模結(jié)構(gòu),實線和虛線分別表示擾動速度的余弦和正弦分量Fig.9.Spatial profile of n=4 TAE,where cosine and sine part of radial velocity are depicted by solid and dash curve,respectively.

圖10 不同能量的EP 在(ρ,v///v)空間的損失錐Fig.10.Loss boundary in (ρ,v///v) space with different energies.

將模擬得到的閾值剖面和損失錐代入到EPtran 程序中后,計算得到的EP 再分布剖面如圖11 中的紅線所示,同樣,為了便于比較,圖中還給出了經(jīng)典慢化分布(黑色)、實驗分布(綠色)及舊模型的模擬結(jié)果(藍色).改進后的模擬結(jié)果呈類似高斯分布的形式,與實驗結(jié)果非常吻合,模擬曲線幾乎都在誤差棒內(nèi).

圖11 壓強剖面對比圖,經(jīng)典慢化分布為黑色曲線,舊臨界梯度模型為藍色曲線,改進的臨界梯度模型為紅色曲線,實驗結(jié)果為綠色三角并配有誤差棒Fig.11.Pressure profile of classical slowing down (black),previous (blue) and improved (red) CGM.For comparison,experimental data is depicted by green triangles with error bar.

4.3 基于神經(jīng)網(wǎng)絡(luò)構(gòu)建HL-3 的高能量粒子模擬模塊

目前OMFIT 集成模擬中,EP 的計算主要是用ONETWO (α 粒子)和NUBEAM (NBI)程序?qū)崿F(xiàn).但這兩個程序使用的是經(jīng)典慢化模型,模擬結(jié)果與實際分布相差很大.因此本章重點討論如何將臨界梯度模型嵌入集成模擬中.根據(jù)第2 節(jié)的模型,可知利用臨界梯度模型預(yù)測EP 分布分為兩步,第1 步是利用TGLFEP 計算臨界梯度,第2步是利用EPtran 程序計算EP 的再分布.目前EPtran 程序已經(jīng)改寫為并行程序,～5 min 就可以完成計算,因此,只有縮短第1 步的運行時間,才有可能將臨界梯度模型加入到集成模擬的迭代中,這里就選用神經(jīng)網(wǎng)絡(luò)的方法來實現(xiàn).

搭建神經(jīng)網(wǎng)絡(luò)第1 步需要合適的數(shù)據(jù)庫,樣本的輸入太多會導(dǎo)致模型預(yù)測不準(zhǔn)確,因此從物理出發(fā),盡量減少輸入變量.這里不考慮雜質(zhì)的影響,因此有效電荷數(shù)Zeff=1,而EP 的密度是掃描量,因此離子的密度可以通過準(zhǔn)中性條件來計算.為了方便,認為離子溫度與電子溫度相同.參考TGLFEP的輸入后,共有18 個輸入變量 (表1).模型的輸入則為3.1 節(jié)中方程(9)中的兩個系數(shù)k1和k2.根據(jù)輸入輸出,參照中國環(huán)流3 號(HL-3)的參數(shù)范圍,利用EFIT 生成150 個不同的平衡,每個平衡隨機選取8 個磁面的參數(shù)作為樣本,因此總共生成1200 個樣本.

表1 神經(jīng)網(wǎng)絡(luò)輸入變量Table 1.Variables for NN input.

神經(jīng)網(wǎng)絡(luò)模型采用最簡單的線性結(jié)構(gòu),主要有兩個功能: 1)預(yù)測是否存在不穩(wěn)定的AE;2)對于存在不穩(wěn)定AE 的磁面,計算k1和k2.流程圖如圖12 所示.

圖12 神經(jīng)網(wǎng)絡(luò)流程圖,第1 步判斷AE 是否被激發(fā),第2步計算系數(shù)k1和k2Fig.12.Flow chart of the neutral network (NN).The NN estimates if the AE can be excited at first,and calculates the two coefficients k1 and k2 for AE unstable location.

將1200 個樣本進行隨機洗牌,80%用來訓(xùn)練模型,剩下的20%用來對模型進行檢驗.第一階段是判斷是否存在不穩(wěn)定的AE,圖13 表示在模型學(xué)習(xí)200 次之后,均方差大約收斂為0.1,準(zhǔn)確率收斂到0.9 左右.繼續(xù)訓(xùn)練會發(fā)生過擬合.對于存在AE 的樣本,同樣選擇80%作為訓(xùn)練集,20%為驗證集.圖14 展示了預(yù)測值與目標(biāo)值之間的差距((a)k1和(b)k2),圓點越接近對角線(虛線)表示預(yù)測值越接近目標(biāo)值.從圖14 可以看到,迭代次數(shù)越高,預(yù)測值就越接近,但總體變化不大.

圖13 (a)訓(xùn)練集和驗證集的均方差隨學(xué)習(xí)次數(shù)的變化;(b)準(zhǔn)確度隨學(xué)習(xí)次數(shù)的變化Fig.13.(a) Loss (mean square error) and (b) accuracy for AE stability estimation versus training epoch.

圖14 預(yù)測值與目標(biāo)值的對比圖 (a) k1;(b) k2Fig.14.Predicted (a) k1 and (b) k2 compared with the targets.

此外,還額外制作了3 個HL-3 的平衡來驗證這個神經(jīng)網(wǎng)絡(luò)模型,3 個平衡具有不同的安全因子剖面,分別為單調(diào)遞增、芯部弱反磁剪切、強反磁剪切,這3 個平衡的壓強和安全因子剖面如圖15所示.此外,還根據(jù)高能量粒子臨界梯度與高能量粒子密度成反比這一物理,采用了一個自定義的損失函數(shù)作為對比.圖16為k1和k2的預(yù)測值和目標(biāo)值(TGLFEP 計算值)的徑向分布,可以看到預(yù)測值與目標(biāo)值相比誤差不大,判斷的AE 的范圍基本相同.將二者代入EPtran 中計算的高能量粒子剖面如圖17 所示,這說明訓(xùn)練的神經(jīng)網(wǎng)絡(luò)可以替代TGLFEP.

圖15 用來驗證神經(jīng)網(wǎng)絡(luò)的3 個平衡的(a)壓強剖面和(b)安全因子剖面Fig.15.(a) Pressure and (b) safety factor profile of three additional equilibria for NN validation.

圖16 HL-3 神經(jīng)網(wǎng)絡(luò)預(yù)測值(MSE 為藍色,自定義損失函數(shù)為紅色)與TGLFEP 計算值(黑色)對比圖 (a),(b) Case 1;(c),(d) Case 2;(e),(f) Case 3Fig.16.Coefficients of k1 and k2 predicted by NN with loss function of MSE (blue) and custom loss function (red): (a),(b) Case 1;(c),(d) Case 2;(e),(f) Case 3.For comparison,TGLFEP results are depicted by black curve.

圖17 具有(a)單調(diào)遞增、(b)芯部弱反磁剪切和(c)強反磁剪切安全因子剖面的平衡位形下,高能量粒子剖面對比圖.綠色曲線為經(jīng)典慢化模型計算的初始EP 剖面,黑色、藍色、紅色曲線分別為根據(jù)TGLFEP、MSE 和自定義損失函數(shù)得到的臨界梯度計算的EP 剖面Fig.17.EP profile comparison for the equilibrium with (a) monotonic,(b) weak and (c) strong shear q-profile.In each panel,green curve depicts initial EP profile with classical slowing down distribution and black curve depicts EP profile with the critical gradient calculated by TGLFEP.The blue and red curves depict EP profiles by NN with loss function of MSE and custom loss function,respectively.

最后采用神經(jīng)網(wǎng)絡(luò)模型代替TGLFEP 程序,與并行的EPtran 程序組合,構(gòu)建MOE.為了提高運行速度,MOE 中使用解析公式[30]代替ORBIT程序計算EP 的損失錐.圖18 為OMFIT 工作流,右側(cè)的TGYRO,ONETWO 和EFIT 三個程序構(gòu)成了目前常用的集成模擬(不考慮臺基區(qū)).TGYRO計算各粒子的密度、溫度剖面;ONETWO 計算電流和粒子源;EFIT 構(gòu)建等離子體平衡.我們搭建的MOE 則讀取TGYRO 計算的粒子剖面信息、EFIT 生成的平衡的磁場信息,以及ONETWO 計算的中性束粒子源(α 粒子源則根據(jù)DT 的剖面進行計算),計算得到EP 的壓強和電流剖面?zhèn)鬟f給ONETWO,然后計算得到總壓強和總電流(P'為壓強的導(dǎo)數(shù),即壓強梯度,F是一個與極向電流相關(guān)的函數(shù),F'為F的導(dǎo)數(shù))再傳遞給EFIT 完成迭代.

圖18 OMFIT 集成模擬流程圖,藍色框里是舊OMFIT 迭代流程,紅色框中的是高能量粒子模塊Fig.18.OMFIT workflow with MOE.The typical iteration is in the blue border,and MOE is in the red border.

利用帶有MOE 的集成模擬計算的等離子體平衡如圖19 所示,其中綠線僅考慮EP 輸運對總壓強的影響,紅線考慮EP 輸運對壓強和電流的影響,藍線則為不包含MOE 的集成模擬的結(jié)果.可以看到EP 輸運導(dǎo)致芯部的總壓強和總電流下降,從而使得芯部的安全因子增大,在一定程度上有利于反磁剪切位形的形成.

圖19 利用帶有MOE 的集成模擬計算的(a)磁面、(b)總壓強、(c)安全因子剖面,其中綠線僅考慮了EP 輸運對總壓強的影響,紅線考慮了EP 輸運對壓強和電流的影響,藍線則為不包含MOE 的集成模擬的結(jié)果Fig.19.(a) Flux surface,(b) total pressure and (c) safety factor profile calculated by the OMFIT integrated simulation with MOE.Green curve only considers pressure modification,red curve considers both pressure and current modification,and blue curve is calculated without MOE.

5 結(jié)論

本文介紹了臨界梯度模型以及使用TGLFEP和EPtran 的組合模擬AE 導(dǎo)致EP 再分布的方法.兩點改進使模擬結(jié)果更接近實驗結(jié)果: 第一,根據(jù)解析理論,將判斷臨界穩(wěn)定的物理量從密度梯度改為歸一化的密度梯度,并修改相應(yīng)的擴散系數(shù)表達式.同時優(yōu)化閾值掃描的數(shù)值方法,令臨界的歸一化密度梯度成為EP 密度的反比例函數(shù),系數(shù)由TGLFEP 計算得到,從而使閾值可以隨著EP的變化而變化.第二,使用ORBIT 計算EP 在相空間的損失錐,并輸入到EPtran 程序中,從而增加一種由FOW 效應(yīng)導(dǎo)致的EP 損失通道.改進后的臨界梯度模型使用#142111 和#153071 兩次實驗進行驗證,模擬得到的EP 剖面與MEGA 結(jié)果類似,與舊模型相比更接近實驗結(jié)果.EP 的損失率也達到42%,這與實驗觀測的結(jié)果很接近.

在此基礎(chǔ)上,還利用神經(jīng)網(wǎng)絡(luò)代替TGLFEP,實現(xiàn)了EPtran 程序的并行計算,加快了整個流程的計算速度,并為HL-3 的集成模擬搭建了一個高能量粒子模塊(MOE).集成模擬的結(jié)果表示當(dāng)AE驅(qū)動EP 輸運,使得芯部壓強和電流減小,從而影響整個平衡.