孫梓源 王元震2) 劉悅?

1)(大連理工大學(xué)物理學(xué)院,三束材料改性教育部重點(diǎn)實(shí)驗(yàn)室,大連 116024)

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

基于HL-2A 實(shí)驗(yàn)參數(shù),利用TOQ 程序構(gòu)建了具有不同臺(tái)基結(jié)構(gòu)的平衡,在BOUT++三場模塊下對(duì)臺(tái)基磁流體力學(xué)(magnetohydrodynamics,MHD)穩(wěn)定性進(jìn)行數(shù)值模擬研究.線性模擬表明,減小臺(tái)基高度、增大臺(tái)基寬度、減小臺(tái)基電流能夠提高臺(tái)基MHD 穩(wěn)定性,利用色散關(guān)系理論,對(duì)上述現(xiàn)象進(jìn)行解釋.在MHD穩(wěn)定性的前提下,預(yù)測了不同臺(tái)基寬度對(duì)應(yīng)的最高臺(tái)基高度,對(duì)數(shù)據(jù)進(jìn)行擬合,得到可以預(yù)測臨界臺(tái)基高度的公式,并在此基礎(chǔ)上結(jié)合動(dòng)理學(xué)氣球模(kinetic ballooning mode,KBM)理論,同時(shí)預(yù)測了臺(tái)基高度和寬度.本文也研究了臺(tái)基結(jié)構(gòu)對(duì)MHD 不穩(wěn)定模式的影響,線性模擬表明,臺(tái)基高度能夠微弱影響不穩(wěn)定模式的徑向模展寬;非線性模擬表明,不穩(wěn)定模式的前期增長主要受單一主導(dǎo)模的影響,模式增長到一定大小會(huì)發(fā)生臺(tái)基坍塌,爆發(fā)邊緣局域模(edge localized mode,ELM),ELM 尺寸的演化與主導(dǎo)模幅值的演化同步,總體來說具有較大線性增長率的平衡在非線性模擬中具有更大的ELM 尺寸和更廣范圍的臺(tái)基坍塌.

1 引言

高約束模式(H 模),作為托卡馬克中一種較高參數(shù)的運(yùn)行模式,從1982 年德國ASDEX 裝置的實(shí)驗(yàn)中首次發(fā)現(xiàn)以來,已被廣泛研究,因其具有較好的約束性能,被ITER 確定為標(biāo)準(zhǔn)放電方案[1?3].H 模最典型的特征是在等離子體邊界區(qū)域存在邊界輸運(yùn)壘(edge transport barrier,ETB),即在該區(qū)域存在較高的溫度和密度梯度,形成一個(gè)類似于臺(tái)基結(jié)構(gòu),被稱為臺(tái)基區(qū)[4].臺(tái)基區(qū)較大的等離子體梯度存儲(chǔ)大量的自由能,容易激發(fā)各種MHD 不穩(wěn)定性,包括剝離模、氣球模、剝離-氣球模和電阻氣球模等[5?7],其中最重要的是電流驅(qū)動(dòng)的剝離模與壓強(qiáng)梯度驅(qū)動(dòng)的氣球模耦合而成的剝離-氣球模,這種不穩(wěn)定性是大型ELM 爆發(fā)的主要原因[8].爆發(fā)ELM時(shí),臺(tái)基坍塌,同時(shí)向外噴發(fā)粒子和熱流,這些粒子和熱流有可能會(huì)使面向等離子體的材料和組件嚴(yán)重?fù)p壞[9].根據(jù)爆發(fā)頻率的高低、加熱功率的大小、能量損失的大小、模結(jié)構(gòu)以及是否出現(xiàn)前兆震蕩等可將ELM 分為Type I,Type II,Type III,Type IV和grassy ELM等[10].

理論研究發(fā)現(xiàn),s-α模型和j-α模型可以判斷平衡的剝離-氣球模不穩(wěn)定性[5].當(dāng)壓強(qiáng)梯度一定時(shí),增大電流容易造成剝離模不穩(wěn)定;而壓強(qiáng)梯度較大時(shí),容易造成氣球模不穩(wěn)定;如果壓強(qiáng)梯度和電流都超過穩(wěn)定性閾值,兩種模式會(huì)耦合成剝離-氣球模,可能會(huì)爆發(fā)具有較大能量損失的Type I ELM[11];因此臺(tái)基結(jié)構(gòu)和臺(tái)基電流與邊界MHD穩(wěn)定性密切相關(guān).目前有很多程序能夠計(jì)算邊界MHD 穩(wěn)定性,也有些模型能預(yù)測臺(tái)基結(jié)構(gòu),例如ELITE[12],GATO[13],MISHKA[14]和BOUT++[15]等程序能計(jì)算邊界MHD 穩(wěn)定性,EPED 等模型能預(yù)測臺(tái)基結(jié)構(gòu)[16].使用ELITE 程序能夠計(jì)算得到剝離-氣球模的MHD 穩(wěn)定性邊界,比如在DIII-D等裝置上發(fā)現(xiàn)壓強(qiáng)臺(tái)基高度與臺(tái)基寬度接近線性相關(guān)[17],在確定臺(tái)基寬度的前提下預(yù)測了臺(tái)基高度[12].EPED 模型在使用ELITE 計(jì)算剝離-氣球模的前提下,考慮了動(dòng)理學(xué)氣球模,同時(shí)考慮這兩種不穩(wěn)定模式時(shí)能夠同時(shí)預(yù)測臺(tái)基高度和寬度[18],李凱[19]在EPED 模型的基礎(chǔ)上發(fā)展了REPED 模型,預(yù)測了EAST 裝置真實(shí)位形下H 模放電的臺(tái)基結(jié)構(gòu).本文使用的BOUT++是一個(gè)大型流體模擬程序框架,已經(jīng)發(fā)展了三場[20]、四場[21]、五場[22]、六場[23]等模塊,這些模塊既能判斷理想MHD 穩(wěn)定性,又加入了抗磁效應(yīng)、E×B漂移、電阻、高階電阻等非理想效應(yīng),可以進(jìn)行線性和非線性的模擬研究[24].

HL-2A 托卡馬克在2009 年首次實(shí)現(xiàn)了偏濾器位形下ELMy H 模放電運(yùn)行[25],針對(duì)HL-2A裝置,已有很多實(shí)驗(yàn)和模擬工作開展了對(duì)臺(tái)基MHD 不穩(wěn)定性的研究.在實(shí)驗(yàn)工作上,發(fā)展了多個(gè)緩解ELM的方法,如低雜波驅(qū)動(dòng)、超聲分子束注入、共振磁擾動(dòng)和雜質(zhì)注入等[26?28];在模擬工作上,Tang等[29]利用BOUT++六場模塊解釋了密度梯度可以觸發(fā)ELM 期間的準(zhǔn)相干模,Wu等[30]利用剝離-氣球模理論對(duì)等離子體垂直擺動(dòng)過程中的ELM 特性進(jìn)行了研究.然而對(duì)于HL-2A 裝置,臺(tái)基結(jié)構(gòu)影響MHD 穩(wěn)定性的研究仍然較少,因此針對(duì)實(shí)驗(yàn)和模擬的需求,本工作將基于HL-2A的實(shí)驗(yàn)參數(shù),給出大量的具有不同臺(tái)基結(jié)構(gòu)的臺(tái)基剖面,利用BOUT++中的MHD 程序做穩(wěn)定性的線性掃描,得到了在不同臺(tái)基寬度下臺(tái)基高度的極限值,并對(duì)數(shù)據(jù)進(jìn)行擬合,得到了臺(tái)基寬度和臨界臺(tái)基高度的關(guān)系式,最后在此基礎(chǔ)上結(jié)合了動(dòng)理學(xué)氣球模理論,同時(shí)預(yù)測臺(tái)基寬度和高度;針對(duì)觸發(fā)MHD 不穩(wěn)定性的臺(tái)基結(jié)構(gòu),也進(jìn)行了非線性模擬研究,對(duì)產(chǎn)生ELM的演化過程進(jìn)行了分析.另外,本工作模擬得到的數(shù)據(jù)將會(huì)整合到HL-2A的集成模擬平臺(tái)中,促進(jìn)HL-2A 集成模擬工作的開展.

本文首先介紹了本工作采用的模型及參數(shù)設(shè)置介紹,對(duì)具有不同臺(tái)基結(jié)構(gòu)的平衡進(jìn)行了MHD穩(wěn)定性分析,并結(jié)合動(dòng)理學(xué)氣球模理論得到了可以預(yù)測HL-2A 臺(tái)基高度和寬度的方法,對(duì)于MHD 不穩(wěn)定的平衡,也對(duì)ELM的演化進(jìn)行了分析和總結(jié).

2 模型

選取了HL-2A 托卡馬克H 模放電的一般實(shí)驗(yàn)參數(shù)為基礎(chǔ)進(jìn)行研究.主要參數(shù)包括:大半徑R0=1.65 m,小半徑a=0.4 m,等離子體中心的環(huán)向場B0=1.27 T.為了便于模擬,忽略x點(diǎn),圖1是忽略了x點(diǎn)后的等離子體邊界區(qū)域,即本文的模擬區(qū)域,并且在模擬中,分形面以外的區(qū)域被近似為真空,等離子體電流近似為零.

圖1 模擬區(qū)域內(nèi)的歸一化極向磁通分布Fig.1.Normalized poloidal magnetic flux in the simulated region.

基于實(shí)驗(yàn)參數(shù),首先利用TOQ 程序根據(jù)構(gòu)建了大量的不同臺(tái)基結(jié)構(gòu)和臺(tái)基電流的平衡,它們具有不同的臺(tái)基高度Pped、臺(tái)基寬度Δped和邊界安全因子q剖面,臺(tái)基部分的壓強(qiáng)剖面被設(shè)置為

其中,ae用于確定臺(tái)基的高度,ωe用于確定臺(tái)基的寬度,常數(shù)const 用于確定邊界處的壓強(qiáng),ψ是歸一化的極向磁通,ψe用于控制臺(tái)基的位置.

利用BOUT++程序中的三場模塊對(duì)平衡的MHD 穩(wěn)定性進(jìn)行了模擬,三場模型包括渦度方程(2)、能量方程(3)和歐姆方程(4)[15,24,31]:

在方程組(2)—(4)中,帶有上標(biāo)“～”的物理量是擾動(dòng)量,帶有下標(biāo)“0”的物理量是平衡量,其中

對(duì)于物理量F,?//F=B?//(F/B),而[f,g]=b0·(?f×?g)是泊松括號(hào)項(xiàng),此外還有:

(5)式中的右端第2 項(xiàng)為抗磁項(xiàng),考慮此項(xiàng)時(shí)需要保持等離子體平衡流為0,即:

因此:

在理想MHD 模擬中,不考慮抗磁效應(yīng),此時(shí):

同時(shí)在理想MHD 模擬中,電阻η、高階電阻ηH、熱擴(kuò)散系數(shù)χ//、離子平行黏度μi,//和泊松括號(hào)項(xiàng)都為零.在非線性模擬中,加入了抗磁效應(yīng)、電阻(倫奎斯特?cái)?shù)S= μ0R0vA/η=108)、高階電阻(高階倫奎斯特?cái)?shù)和離子平行黏度其中vA為阿爾芬速度,阿爾芬時(shí)間τA=R0/vA≈0.41 μs,這些非理想效應(yīng)的加入,有利于數(shù)值收斂.為了提高模擬效率,在線性模擬中,當(dāng)環(huán)向模數(shù)為n時(shí),1/n的環(huán)向范圍被模擬,初始擾動(dòng)的環(huán)向模數(shù)也被設(shè)置為n;在非線性模擬中,1/5的環(huán)向范圍被模擬,初始擾動(dòng)的環(huán)向模為混合模.線性模擬的網(wǎng)格點(diǎn)數(shù)為260× 64× 17,非線性模擬的為260× 64× 65.

模擬中使用了沿磁力線的坐標(biāo)系(x,y,z)[32],x,y,z分別為徑向、極向、環(huán)向坐標(biāo),環(huán)向和極向?yàn)檠h(huán)邊界,徑向的邊界條件為

徑向內(nèi)邊界:

徑向外邊界:

3 模擬結(jié)果

3.1 線性模擬結(jié)果

首先固定安全因子q剖面不變,改變臺(tái)基的高度和寬度.臺(tái)基高度參數(shù)規(guī)定為臺(tái)基環(huán)向比壓其中Pped表示臺(tái)基頂部的壓強(qiáng),臺(tái)基寬度參數(shù)規(guī)定為在ψ坐標(biāo)下的歸一化臺(tái)基寬度?ψ=?ped/a,其中a為等離子體小半徑.

上文提到,線性模擬的網(wǎng)格點(diǎn)數(shù)為260× 64×17,為證明選取的分辨率對(duì)線性模擬足夠精確.對(duì)于臺(tái)基寬度和高度分別為?ψ=0.088,βt=6.0×10?3的平衡,不同分辨率的網(wǎng)格被生成,該平衡的壓強(qiáng)、平行電流以及安全因子剖面如圖2所示.將徑向nx,極向ny和環(huán)向nz的網(wǎng)格點(diǎn)數(shù)分別設(shè)定為:36,68,132,260,516;32,64,128,256,512 以及17,33,65,129,257.改變n x時(shí),n y和n z固定為64和17,改變n y時(shí),n x和n z固定為260和17,改變n z時(shí),n x和n y固定為260和64,計(jì)算環(huán)向模數(shù)n=15 時(shí)的線性增長率,結(jié)果如圖3 所示,其中γ/ωA代表歸一化線性增長率.從圖3 可以看出,相對(duì)于每組最大分辨率,選擇分辨率為260× 64×17 計(jì)算得到的結(jié)果誤差均小于0.3%,所以該分辨率對(duì)于線性模擬足夠精確.

圖2 臺(tái)基高度 βt=6.0×10?3,臺(tái)基寬度?ψ=0.088 時(shí)的壓強(qiáng)(a)、平行電流和安全因子(b)剖面Fig.2.Pressure(a),parallel current density and safety factor profiles(b)with pedestal height βt=6.0×10?3 and pedestal width ?ψ=0.088.

圖3 利用不同分辨率的網(wǎng)格得出的線性增長率Fig.3.Linear growth rates corresponding to different resolutions.

圖4是臺(tái)基高度和寬度不同時(shí)對(duì)應(yīng)的臺(tái)基區(qū)壓強(qiáng)、平行電流以及安全因子剖面.臺(tái)基高度βt的取值為2.0×10?3,3.0×10?3和4.0×10?3;臺(tái)基寬度?ψ取值為0.066,0.077,0.088,0.099和0.110,在圖4(b)和(c)中可以看到,隨著臺(tái)基高度的增大或者臺(tái)基寬度的減小,對(duì)應(yīng)自洽的平行電流剖面升高.

圖4 臺(tái)基高度 βt為2 .0×10?3,3.0×10?3和4 .0×10?3,臺(tái)基寬度 ?ψ為0.066,0.077,0.088,0.099和0.110 時(shí)的壓強(qiáng)(a)、平行電流和安全因子(b)(c)剖面Fig.4.Pressure(a),parallel current density and safety factor profiles(b)(c)with pedestal heights βt of 2.0×10?3,3.0×10?3 and 4 .0×10?3,and pedestal widths ?ψ of 0.066,0.077,0.088,0.099 and 0.110.

計(jì)算具有不同臺(tái)基結(jié)構(gòu)平衡的理想MHD 穩(wěn)定性,結(jié)果如圖5 所示,發(fā)現(xiàn)隨著臺(tái)基高度的增大或?qū)挾鹊臏p小,線性增長率增大,這是因?yàn)楦淖兣_(tái)基的高度和寬度相當(dāng)于改變了臺(tái)基梯度,在s-α模型中,壓強(qiáng)梯度增大會(huì)增大氣球模的不穩(wěn)定性[33],并且在壓強(qiáng)梯度變大的同時(shí),等離子體電流也隨之增大,這些都會(huì)導(dǎo)致臺(tái)基區(qū)等離子體的自由能增大,從而使增長率變大.為解釋該現(xiàn)象,在Huuang等[34]研究的基礎(chǔ)上加入電流驅(qū)動(dòng)項(xiàng),推導(dǎo)得到了理想剝離-氣球模的色散關(guān)系,推導(dǎo)過程見附錄A.其中,n是環(huán)向模數(shù),k//是平行波矢,k⊥是垂直波矢,對(duì)于本文所研究的剝離-氣球模,k//?k⊥.由(16)式可知,與? P0和? J//0成正相關(guān).增大臺(tái)基高度,減小臺(tái)基寬度時(shí),不僅? P0增大,而且對(duì)應(yīng)的平行電流梯度? J//0也增大,? P0驅(qū)動(dòng)氣球模,而? J//0驅(qū)動(dòng)剝離模,所以改變臺(tái)基結(jié)構(gòu)對(duì)剝離模成分和氣球模成分都有影響.(16)式也可以說明,氣球模成分與環(huán)向模數(shù)n成正比,而剝離模成分與環(huán)向模數(shù)n成反比,圖5 中總的剝離-氣球模的線性增長率γpb(n)∝n,這說明氣球模和剝離-氣球模主導(dǎo)的中高n模增長率較大,而剝離模主導(dǎo)的低n模增長率較小.

圖5 理 想MHD 線性增長率(a)臺(tái)基高度 βt 分別為2.0×10?3,3.0×10?3和4.0× 10–3,臺(tái)基寬度 ?ψ=0.110;(b)臺(tái)基高度 βt=3.0×10?3,臺(tái)基寬度 ?ψ分別為0.066,0.077,0.088,0.099和0.110Fig.5.Linear growth rates with different pedestals:(a)Pedestal heights βt are 2 .0×10?3,3.0×10?3 and 4 .0×10?3,while pedestal width ?ψ is 0.110;(b)pedestal height βt is 3.0×10?3 while pedestal widths ?ψ are 0.066,0.077,0.088,0.099 and 0.110.

中等環(huán)向模數(shù)(n≈ 15)的模同時(shí)具有剝離模成分和氣球模成分,兩者耦合為剝離-氣球模,可能會(huì)造成較大能量損失的Type I ELM[11],為了預(yù)防這種ELM的產(chǎn)生,定義了相對(duì)于剝離-氣球模穩(wěn)定的臨界平衡,即在n≤15時(shí)γpb=0的平衡.如圖6所示,臨界平衡的臺(tái)基高度隨臺(tái)基寬度的增大近似線性增大,這是因?yàn)樵龃笈_(tái)基寬度,臺(tái)基不穩(wěn)定性減小,臨界臺(tái)基高度增大.

圖6 臨界臺(tái)基高度 βt 隨臺(tái)基寬度 ?ψ的變化Fig.6.The critical pedestal heights βt corresponding to different pedestal widths ?ψ.

選取具有不同臺(tái)基高度和寬度的平衡,計(jì)算n=15模式的理想MHD 不穩(wěn)定性,得到了不穩(wěn)定性模式對(duì)應(yīng)的模結(jié)構(gòu)(擾動(dòng)壓強(qiáng)),包括徑向一維模結(jié)構(gòu)和二維模結(jié)構(gòu),如圖7 所示.圖7(a)是固定臺(tái)基寬度?ψ=0.066,臺(tái)基高度不同時(shí)的歸一化徑向模結(jié)構(gòu);圖7(b)是固定臺(tái)基高度βt=4.0×10?3,臺(tái)基寬度不同時(shí)的歸一化徑向模結(jié)構(gòu).從圖7(a)中可以看到,壓強(qiáng)擾動(dòng)的均方根峰值位于壓強(qiáng)梯度最大的位置,隨著臺(tái)基高度的增大,線性增長率增大,同時(shí)模的徑向展寬也有微弱增大的趨勢;當(dāng)臺(tái)基寬度增大時(shí),壓強(qiáng)梯度最大值的位置左移,從圖7(b)中觀察到,壓強(qiáng)擾動(dòng)均方根峰值也隨之左移.在圖7(c)中,主要發(fā)生不穩(wěn)定性的位置位于模擬域外側(cè)壞曲率磁場處,同時(shí)模結(jié)構(gòu)比較細(xì)密,呈現(xiàn)剝離-氣球模的典型結(jié)構(gòu).

固定臺(tái)基寬度?ψ=0.066、臺(tái)基高度βt=4.0×10?3,當(dāng)q剖面不同時(shí),平衡的剖面分布如圖8 所示.在圖8中,q1—q5五個(gè)安全因子剖面對(duì)應(yīng)的q95值分別為2.15,2.38,2.62,2.86和3.10,圖2—圖7中模擬選擇的平衡的安全因子剖面為圖8(a)中的q3.重復(fù)之前尋找臨界平衡的過程,得到了基于剝離-氣球模理論的臨界臺(tái)基高度隨臺(tái)基寬度的變化,如圖9(a)(PBM(q1—q5))所示.對(duì)于具有不同的q剖面的平衡,計(jì)算得到的臺(tái)基高度幾乎都與臺(tái)基寬度線性相關(guān),因此將數(shù)據(jù)其擬合成βt=k?ψ+b的形式,發(fā)現(xiàn)斜率k與安全因子值q95成負(fù)相關(guān),而截距b的數(shù)值相對(duì)較小,將其作為與q95二次相關(guān)的擬合修正量,如圖9(b)所示.擬合結(jié)果如(17)式、(18)式所示,給定臺(tái)基寬度?ψ和邊界安全因子值q95,可以預(yù)測臨界臺(tái)基高度βt:

圖7 (a)臺(tái)基高度 βt 分別為3.0×10?3,4.0×10?3,5.0×10?3和6 .0×10?3,臺(tái)基寬度 ?ψ為0.066,n=15的歸一化徑向模結(jié)構(gòu);(b)臺(tái)基高度 βt=4.0×10?3,臺(tái)基寬度?ψ 分別為0.066,0.077,0.088和0.099,n=15的歸一化徑向模結(jié)構(gòu);(c)臺(tái)基高度 βt=4.0× 10?3,臺(tái)基寬度?ψ=0.066,n=15的歸一化二維模結(jié)構(gòu)Fig.7.(a)The normalized radial mode structure of n=15 when pedestal heights βt are 3.0×10?3,4.0×10?3,5.0×10?3and 6 .0×10?3,while pedestal width ?ψ is 0.088;(b)the normalized radial mode structure of n=15when pedestal height βt is 4 .0×10?3,pedestal widths ?ψ are 0.066,0.077,0.088 and 0.099 respectively;(c)the normalized two-dimensional mode structure of n=15 when pedestal height βt is 4 .0×10?3 and pedestal width ?ψ is 0 .066.

圖8 臺(tái)基高度 βt=4.0×10?3、臺(tái)基寬度 ?ψ=0.066 時(shí)的壓強(qiáng)、安全因子(a)和平行電流(b)剖面Fig.8.Pressure,safety factor(a)and parallel current density profiles(b)with pedestal height βt=4.0×10?3 and pedestal width ?ψ=0.066.

當(dāng)邊界區(qū)域的安全因子有明顯增大時(shí),會(huì)導(dǎo)致邊界區(qū)域的磁剪切增大,等離子體電流也會(huì)隨之增大,如圖8(b)所示.在圖9(a)中看到,隨著邊界安全因子剖面的升高或邊界電流密度的增大,不穩(wěn)定性增強(qiáng),線性增長率變大,臨界臺(tái)基高度下降.

圖9 (a)不同安全因子的條件下,剝離-氣球模(PBM)穩(wěn)定性限制下的臺(tái)基高度和寬度的關(guān)系,以及動(dòng)理學(xué)氣球模(KBM)穩(wěn)定性限制下的臺(tái)基高度和寬度的關(guān)系,兩者的交點(diǎn)即為臺(tái)基高度和寬度的預(yù)測值;(b)剝離-氣球模(PBM)穩(wěn)定性約束擬合直線的斜率 k和截距b隨q95的變化Fig.9.(a)Relationship between the pedestal height and width under the stability constraints of peeling-ballooning mode(PBM),and relationship between the pedestal height and width under the stability of constraint of kinetic ballooning mode(KBM),the insterection of the two is the predicted value of pedestal height and width;(b)the slope k and intercept b of the fitting lines with the stability constraints of peeling-ballooning model(PBM)vary with q95.

根據(jù)剝離-氣球模的穩(wěn)定性理論,只能在確定臺(tái)基寬度的前提下預(yù)測臺(tái)基高度,為了同時(shí)預(yù)測臺(tái)基高度與寬度,需要另外引入臺(tái)基寬度的約束條件,動(dòng)理學(xué) 氣球模(kinetic ballooning mode,KBM)理論目前被認(rèn)為是可靠的約束條件[18],該理論認(rèn)為,歸一化臺(tái)基寬度?ψ與臺(tái)基極向比壓的 1/2次方即線性相關(guān)[35],即:

其中G是比例系數(shù),與碰撞率v?,環(huán)徑比ε以及其他參數(shù)存在弱相關(guān)關(guān)系,根據(jù)理論分析,設(shè)置此系數(shù)為0.1[35].經(jīng)過推導(dǎo),得到基于KBM 理論的約束邊界表達(dá)式(20),推導(dǎo)過程見附錄B.

根據(jù)HL-2A 托卡馬克H 模放電的一般實(shí)驗(yàn)結(jié)果,將等離子體電流Ip設(shè)定為160 kA,同時(shí)將其他一般性參數(shù),即第2 節(jié)所敘述的參數(shù)代入(20)式,得到基于HL-2A 托卡馬克實(shí)驗(yàn)參數(shù)的KBM 約束邊界表達(dá)式:

在圖9(a)中,動(dòng)理學(xué)氣球模約束曲線(利用KBM 表示)與剝離-氣球模約束擬合直線的交點(diǎn)(利用星號(hào)表示),即為本文預(yù)測的具有不同安全因子剖面的臺(tái)基結(jié)構(gòu).由于時(shí)間所限,本文在對(duì)剝離-氣球模進(jìn)行MHD 穩(wěn)定性計(jì)算時(shí)保持了環(huán)向磁場B0不變,而在實(shí)際托卡馬克H 模放電時(shí)小半徑a,環(huán)向場B0以及等離子體電流Ip都會(huì)在本文選取的參數(shù)附近微小變化,如果想對(duì)HL-2A的H 模放電臺(tái)基結(jié)構(gòu)進(jìn)行更加精準(zhǔn)的預(yù)測,需要在本文工作的基礎(chǔ)上微弱改變參數(shù),進(jìn)行大量的計(jì)算,這也是未來工作的重點(diǎn).

計(jì)算平衡的理想MHD 穩(wěn)定性,固定臺(tái)基寬度?ψ=0.066,臺(tái)基高度,安全因子剖面為圖8(a)中的q3,在考慮抗磁效應(yīng)的條件下,設(shè)定離子密度分別為1.0× 1019,2.0× 1019,3.0× 1019,3.0× 1019,3.0× 1019,4.0× 1019,以及5.0× 1019m–3,計(jì)算了MHD 穩(wěn)定性,結(jié)果如圖10 所示.可以看到,密度越小,抗磁效應(yīng)對(duì)不穩(wěn)定性的抑制作用越明顯,特別是環(huán)向模數(shù)較大時(shí),這是因?yàn)殡x子抗磁頻率ωi=?b0×?Pi0·k⊥/n0eB0∝n/n0,即ωi與環(huán)向模數(shù)n成正比,與等離子體密度n0成反比.在(16)

圖10 在考慮抗磁效應(yīng)時(shí),具有不同離子密度的剝離-氣球模增長率,臺(tái)基高度βt=4.0×10?3,臺(tái)基寬度?ψ=0.066Fig.10.When the diamagnetic drift effect is included,linear growth rates of peeling-ballooning modes with different ion densities,while pedestal height βt=4.0×10?3 and pedestal width ?ψ=0.066.

式中,等式右邊的前兩項(xiàng)是壓強(qiáng)驅(qū)動(dòng)項(xiàng),即氣球模的驅(qū)動(dòng)項(xiàng),可以看出,氣球模成分的線性增長率(γb)正比于,減小離子密度和增大環(huán)向模數(shù)時(shí),γb增大的幅度沒有離子抗磁頻率ωi增大的幅度大,抗磁致穩(wěn)的相對(duì)效果更強(qiáng),模式被抑制.由(A20)式知,不只是氣球模,如果抗磁效應(yīng)足夠強(qiáng),當(dāng)γpb=ωi/2時(shí),剝離-氣球模的增長率降為零[32].

3.2 非線性模擬結(jié)果

在非線性模擬中,考慮了離子的抗磁漂移,根據(jù)圖10的結(jié)果,為了避免抗磁效應(yīng)過大,將離子密度設(shè)置為ni=5×1019m–3.邊緣局域模(ELM)尺寸被定義為芯部能量的相對(duì)損失[20]:

其中Rin和Rout分別是內(nèi)模擬邊界和最大壓強(qiáng)梯度的徑向位置.當(dāng)臺(tái)基結(jié)構(gòu)和安全因子剖面不同時(shí),ELM尺寸如圖11 所示,其中橫坐標(biāo)表示歸一化的時(shí)間.

圖11 不同臺(tái)基結(jié)構(gòu)和安全因子剖面的ELM 尺寸(a)臺(tái)基高度 βt 分別為2 .0×10?3、2.5 × 10?3和3 .0×10?3,臺(tái)基寬度?ψ分別為0.066,0.077和0.088,安全因子剖面為圖8(a)中的 q3;(b)臺(tái)基高度 βt=3.0×10?3,臺(tái)基寬度 ?ψ=0.088,安全因子剖面分別為圖8(a)中的 q3、q4和q5Fig.11.ELM size corresponding to different pedestal structures and safety factor profiles:(a)βt are 2 .0×10?3,2 .5×10?3 and 3.0×10?3,?ψ=0.066,0.077,0.088,while safety factor profile is the q3 profile in Fig.8(a);(b)βt is 3 .0×10?3,?ψ is 0 .088,while safety factor profiles are the q3,q4 and q5 profiles in Fig.8(a).

從圖11 中可以看到,ELM 尺寸與前面的線性不穩(wěn)定性研究相互對(duì)應(yīng),即臺(tái)基寬度越窄,臺(tái)基高度越高,或者邊界安全因子剖面越高,線性增長率越大,對(duì)應(yīng)的ELM 尺寸越大;并且可以將ELM的演化分為3個(gè)階段,即線性增長階段、臺(tái)基初始坍塌階段和湍流階段.圖12是擾動(dòng)壓強(qiáng)在徑向和環(huán)向上的二維分布隨時(shí)間的演化,結(jié)合圖11 可以看到,在前期的線性增長階段,擾動(dòng)范圍較小,擾動(dòng)量持續(xù)增大,擾動(dòng)壓強(qiáng)的分布如圖12(a)所示;當(dāng)擾動(dòng)量增大到一定幅值,大約t=100τA之后,擾動(dòng)壓強(qiáng)的正值部分向外側(cè)移動(dòng),而負(fù)值部分向芯部移動(dòng),如圖12(b)和(c)所示,此時(shí)臺(tái)基剖面出現(xiàn)坍塌,ELM 尺寸迅速增大;之后進(jìn)入湍流輸運(yùn)階段,模式互相耦合,正負(fù)擾動(dòng)繼續(xù)向兩側(cè)發(fā)展,ELM 尺寸增長的速度變緩,如圖12(d)—(f)所示.從圖11也可以看出,ELM 尺寸增長的越快,越早進(jìn)入湍流階段.在圖13 中將臺(tái)基高度βt固定為3.0×10?3,顯示了最外中平面上歸一化擾動(dòng)壓強(qiáng)均方根隨時(shí)間演化的徑向分布;從圖13(a)和(b)可以看到,在安全因子剖面相同時(shí),臺(tái)基寬度較窄時(shí)擾動(dòng)壓強(qiáng)消散的更快;而在圖13(b)—(d)中看到,在臺(tái)基寬度相同時(shí),邊界安全因子剖面越高,擾動(dòng)壓強(qiáng)消散的越快;同時(shí)對(duì)于消散更快的算例,擾動(dòng)相對(duì)平衡位置向兩側(cè)擴(kuò)張的范圍也有增大的趨勢;對(duì)比圖11 發(fā)現(xiàn)擾動(dòng)大約在剛進(jìn)入湍流階段時(shí)刻開始消散,所以對(duì)于ELM 尺寸增長越快的算例,擾動(dòng)消散的也越快.

圖12 當(dāng)臺(tái)基高度 βt=3.0×10?3,臺(tái)基寬度 ?ψ=0.088,安全因子剖面為圖8(a)中的 q4時(shí),歸一化的擾動(dòng)壓強(qiáng)在徑向 ψ和環(huán)向 ζ 平面隨時(shí)間的演化Fig.12.Evolution of pressure perturbation in the frame of normalized poloidal flux ψ and toroidal angle ζ with pedestal height βt=3.0×10?3 and pedestal width ?ψ=0.088,while safety factor profile is the q4 profile in figure 8(a).

圖13 固定臺(tái)基高度 βt=3.0×10?3,當(dāng)臺(tái)基寬度和安全因子剖面不同時(shí),外中平面處歸一化的擾動(dòng)壓強(qiáng)均方根隨時(shí)間的演化Fig.13.Time evolution of the root-mean-square of pressure perturbation at the outer mid-plane with βt=3.0×10?3,different pedestal widths and safety factor profiles.

為了觀察各個(gè)環(huán)向模在非線性模擬中的演化,我們?cè)诃h(huán)向進(jìn)行了傅里葉分解,結(jié)果如圖14所示.結(jié)合圖11 看到,在ELM 爆發(fā)前,以n=15為首的中等環(huán)向模數(shù)的模首先迅速增長,而超過30的高n和低于10的低n的環(huán)向模幅值幾乎為0,我們把前期增長最快速的模,即圖14中n=15的模稱為主導(dǎo)模.通過觀察可以發(fā)現(xiàn),圖14(a),(c),(d)和(b)的主導(dǎo)模的峰值逐漸增加,而在圖11中對(duì)應(yīng)的ELM 尺寸也逐漸增加,因此主導(dǎo)模的峰值基本決定了ELM的尺寸.ELM 爆發(fā)后,主導(dǎo)模及中等環(huán)向模數(shù)的模的幅值下降,模式之間互相耦合,逐漸進(jìn)入到湍流階段,其中n=0的帶狀流的發(fā)展明顯抑制了以主導(dǎo)模為首的其他環(huán)向模的增長,可以看到,圖14(a)中的抑制效果最明顯,圖14(d)中的抑制效果最不明顯,環(huán)向傅里葉分解結(jié)果與吳[6]等的研究結(jié)果相似.

圖14 不同臺(tái)基結(jié)構(gòu)和安全因子剖面時(shí)擾動(dòng)的環(huán)向傅里葉分解圖Fig.14.Toroidal Fourier analysis of the perturbation with different pedestal structures and safety factor profiles.

圖15是圖11對(duì)應(yīng)算例的主導(dǎo)模的演化圖,即n=15的模隨時(shí)間的演化,其中圖15(a)中的算例具有相同的邊界安全因子剖面,即圖8(a)中的q3剖面,以及不同的臺(tái)基結(jié)構(gòu);圖15(b)中的算例具有相同的臺(tái)基高度βt=3.0×10?3和臺(tái)基寬度?ψ=0.088,以及不同的邊界安全因子剖面,即圖8(a)中的q3—q5剖面;與圖11 比較,ELM 尺寸的演化與主導(dǎo)模幅值的演化相一致,主導(dǎo)模的峰值大約出現(xiàn)在臺(tái)基初始坍塌階段的中期,同時(shí),主導(dǎo)模在線性增長階段增長的越快,峰值越高,ELM 爆發(fā)的越快,尺寸也越大,因此可以通過主導(dǎo)模式的增長率判斷ELM爆發(fā)的快慢和尺寸的大小.

圖15 不同臺(tái)基結(jié)構(gòu)和安全因子剖面的等離子體歸一化的主導(dǎo)模隨時(shí)間的演化Fig.15.Time evolution of normalized dominant toroidal modes with different pedestal structures and safety factor profiles.

根據(jù)圖11 中ELM 尺寸以及圖12 中擾動(dòng)的隨時(shí)間演化情況,選取了具有不同臺(tái)基結(jié)構(gòu)和安全因子剖面的算例,圖16 所示為t=195τA時(shí)的歸一化擾動(dòng)壓強(qiáng)的二維分布,其中圖16(a)和(b)的臺(tái)基高度βt由2.0×10?3升高到 3.0×10?3,圖16(b)—(d)的臺(tái)基寬度?ψ由0.066 增大到0.088,圖16(d)—(f)的邊界安全因子剖面由圖8(a)中的q3升高到q5.可以看到臺(tái)基高度越高,臺(tái)基寬度越窄,邊界安全因子剖面越高,擾動(dòng)壓強(qiáng)越有向平衡位置兩側(cè)微弱擴(kuò)張的趨勢,對(duì)應(yīng)的臺(tái)基不穩(wěn)定性增強(qiáng),在圖11中的體現(xiàn)為ELM的尺寸也增大,這些說明圖16與圖13 得到了相似的結(jié)論.圖17 所示為磁面平均的總壓強(qiáng)(平衡壓強(qiáng)與擾動(dòng)壓強(qiáng)之和)同樣在t=195τA時(shí)的分布剖面,3 個(gè)算例具有相同的初始?jí)簭?qiáng)剖面,根據(jù)圖11,它們均處于湍流階段,可以看到,當(dāng)邊界安全因子剖面越高時(shí),不穩(wěn)定性越強(qiáng),臺(tái)基坍塌的范圍越大,代表著擾動(dòng)壓強(qiáng)向兩側(cè)擴(kuò)張的范圍也越大.

圖16 當(dāng)臺(tái)基結(jié)構(gòu)和安全因子剖面不同時(shí),t=195τA 時(shí)刻歸一化擾動(dòng)壓強(qiáng)在徑向 ψ和環(huán)向 ζ 平面上的分布Fig.16.Pressure perturbation in the frame of normalized poloidal flux ψ and toroidal angle ζ with different pedestal structures and safety factor profiles at t=195τA.

圖17 固定臺(tái)基高度 βt=3.0×10?3和臺(tái)基寬度?ψ=0.088,當(dāng)安全因子剖面不同時(shí),非線性模擬中 t=195τA 時(shí)刻磁面平均的總壓強(qiáng)分布Fig.17.Fixed βt=3.0×10?3,?ψ=0.088,when safety factor profiles are different,the surface-averaged pressure profiles distribution at t=195τA in the nonlinear simulations.

4 總結(jié)

本文基于HL-2A 托卡馬克裝置的實(shí)驗(yàn)參數(shù),生成具有不同臺(tái)基結(jié)構(gòu)和安全因子剖面的平衡,利用BOUT++三場雙流體模塊進(jìn)行了線性和非線性模擬,結(jié)果總結(jié)如下:

1)在線性模擬中,增大臺(tái)基高度與減小臺(tái)基寬度都會(huì)增大臺(tái)基梯度,同時(shí)自洽的平行電流梯度也增大,這些都會(huì)導(dǎo)致剝離-氣球模的增長率變大;安全因子q剖面的升高會(huì)導(dǎo)致平行電流密度及密度梯度的增大,同樣會(huì)導(dǎo)致不穩(wěn)定性增長率變大.

2)n=15 左右的中等環(huán)向模數(shù)的模在模擬中被發(fā)現(xiàn)具有主導(dǎo)意義,這是因?yàn)橹械拳h(huán)向模數(shù)的模耦合了剝離模與氣球模成分,對(duì)臺(tái)基不穩(wěn)定性起到了至關(guān)重要的作用,非線性模擬發(fā)現(xiàn),當(dāng)這些主導(dǎo)模式越不穩(wěn)定時(shí),ELM 尺寸越大.

3)在線性模擬中定義了處于MHD 穩(wěn)定性邊界的臨界平衡,通過計(jì)算具有不同臺(tái)基結(jié)構(gòu)平衡的MHD 穩(wěn)定性,找到了對(duì)應(yīng)不同臺(tái)基寬度的臨界臺(tái)基高度,發(fā)現(xiàn)臺(tái)基高度與寬度大致線性相關(guān),經(jīng)過擬合得到了預(yù)測臨界臺(tái)基高度的經(jīng)驗(yàn)公式,在此基礎(chǔ)上結(jié)合動(dòng)理學(xué)氣球模理論,同時(shí)預(yù)測了HL-2A 托卡馬克的臺(tái)基高度和寬度.

4)線性和非線性模擬都發(fā)現(xiàn)模式不穩(wěn)定性的增強(qiáng)會(huì)使擾動(dòng)的分布范圍更大:在線性模擬中發(fā)現(xiàn),隨著臺(tái)基高度的增大,模的展寬會(huì)微弱的增大,在非線性模擬中發(fā)現(xiàn),不穩(wěn)定性越強(qiáng),擾動(dòng)向兩側(cè)擴(kuò)張的范圍越大,臺(tái)基坍塌范圍也越大.

未來將會(huì)在本文工作的基礎(chǔ)上微弱改變參數(shù),例如小半徑a、環(huán)向場B0、等離子體電流Ip以及邊界安全因子值q95,利用機(jī)器學(xué)習(xí)等方法更準(zhǔn)確地預(yù)測臺(tái)基結(jié)構(gòu),并且會(huì)將結(jié)果集成到HL-2A 集成平臺(tái)下,進(jìn)而為HL-2A 托卡馬克H 模實(shí)驗(yàn)和集成模擬提供參考.

感謝美國勞倫斯利弗莫爾國家實(shí)驗(yàn)室的徐學(xué)橋研究員、英國約克大學(xué)的Dudson 博士以及中國科學(xué)院等離子體物理研究所的夏天陽研究員對(duì)BOUT++代碼的貢獻(xiàn),以及與衡陽師范學(xué)院黃艷清博士的有益討論.

附錄A 剝離-氣球模色散關(guān)系的推導(dǎo)

將(A8)—(A11)式代入(A7)式,得到:

垂直波矢表示為

采用局域化近似,將k//看成常數(shù),且在BOUT++三場模型的渦度表達(dá)式(A6)中,電子動(dòng)量的貢獻(xiàn)被忽略[34].將(A13)式代入(A12)式,整理后得到理想剝離-氣球模的色散關(guān)系:

考慮抗磁項(xiàng)時(shí),渦度變?yōu)?/p>

同樣對(duì)(A15)式進(jìn)行傅里葉變換,得到:

離子抗磁頻率的表達(dá)式為

考慮抗磁項(xiàng)時(shí),增長率用γ表示,將渦度方程傅里葉變換表達(dá)式(A12)中的γpb用γ替代,即為

將(A14),(A16),(A17)式代入(A18)式,同樣采用局域化近似,整理后得到:

(A19)式與Snyder等[8]等得到的一元二次方程形式相同.解方程(A19),取正根,得到:

由(A20)式可以看出,當(dāng)γpb=ωi/2時(shí),增長率降為零.

附錄B KBM 約束邊界表達(dá)式的推導(dǎo)

KBM湍流理論認(rèn)為,歸一化臺(tái)基寬度?ψ與臺(tái)基極向比壓βp,ped之間滿足[19]:

其中G(v?,ε,···)是與碰撞率v?,環(huán)徑比ε以及其他參數(shù)存在弱相關(guān)關(guān)系的系數(shù),根據(jù)理論分析,將其設(shè)置為0.1[35].

臺(tái)基極向比壓βp,ped定義為[19]

其中Pped是臺(tái)基頂部的壓強(qiáng),〈Bp〉是平均極向磁場,定義為

其中Ip是等離子體電流,l是最外磁面的周長,則:

其中a為等離子體小半徑,將(B1)式代入(B3)式:

將(B5)式代入(B2)式:

將(B6)式代入(B1)式:

臺(tái)基環(huán)向比壓βt定義為

其中B0為環(huán)向場,聯(lián)立(B7)式和(B8)式可得: