劉志遠(yuǎn)， 李國強(qiáng)， 徐立清

（中國科學(xué)院 等離子體物理研究所，安徽 合肥 230031）

基于TEQ的固定邊界托卡馬克等離子體平衡模擬

劉志遠(yuǎn)， 李國強(qiáng)， 徐立清

（中國科學(xué)院 等離子體物理研究所，安徽 合肥 230031）

文章設(shè)計(jì)了新的TEQ計(jì)算流程，實(shí)現(xiàn)了由磁面位置及磁面上磁通值計(jì)算網(wǎng)格點(diǎn)上磁通的算法，修正了TEQ求解固定邊界托卡馬克等離子體平衡時(shí)輸出g－file的程序，并以EAST裝置為基礎(chǔ)，用剪切和反剪切2種情況的安全因子、壓強(qiáng)分布進(jìn)行了檢驗(yàn)，發(fā)現(xiàn)TEQ的結(jié)果與EFIT符合得很好，最后用TEQ構(gòu)建了6組反剪切固定邊界平衡，簡單分析了壓強(qiáng)對磁軸水平位置的影響。

TEQ程序；g－file程序；EAST裝置；固定邊界等離子體平衡；數(shù)值模擬

0 引 言

NTCC（National Transport Code Collaboration）［1］中的 TEQ 程 序 是 一個(gè)通過 數(shù) 值 求 解Grad－Shafranov方程來計(jì)算托卡馬克磁流體力學(xué)平衡的程序，它是從Corsica［2］中獨(dú)立出的，可以求解自由邊界平衡，也可以求解固定邊界平衡，是一個(gè)很強(qiáng)大的程序庫［3－4］。

在托卡馬克研究領(lǐng)域，最常用的平衡數(shù)據(jù)文件格式是 EFIT 程序［5］輸出的 g－file，g－file幾乎已經(jīng)成為行業(yè)標(biāo)準(zhǔn)。但NTCC中的TEQ程序不能讀取g－file，因此，本文對TEQ這方面的功能進(jìn)行了補(bǔ)充并重新設(shè)計(jì)了相應(yīng)的程序模塊，對讀入g－file的程序重新編寫，修正輸出g－file的程序中的不足，實(shí)現(xiàn)了由磁面位置及磁面上磁通值計(jì)算網(wǎng)格點(diǎn)上磁通的算法，并設(shè)計(jì)了新的程序流程，調(diào)整了部分參數(shù)。通過數(shù)據(jù)驗(yàn)證，發(fā)現(xiàn)TEQ所得的計(jì)算結(jié)果與EFIT計(jì)算所得的結(jié)果吻合得比較好。用TEQ可以構(gòu)建多種平衡，計(jì)算所得的數(shù)據(jù)可以用來分析壓強(qiáng)對磁軸位置的影響，還可以進(jìn)一步分析平衡的穩(wěn)定性，也可以為構(gòu)建新的平衡提供一種依據(jù)并為設(shè)計(jì)新的裝置提供一種參考。本文所有的計(jì)算都是以EAST（Experimental Advanced Superconducting Tokamak）裝置為基礎(chǔ)進(jìn)行的。

1 托卡馬克等離子體平衡方程

由單流體MHD方程描述的等離子體平衡由（1）式表示，即

其中，p為等離子體壓強(qiáng)；B為磁感應(yīng)強(qiáng)度；j為電流密度。在軸對稱情況下，由（1）式可以得出角向磁通ψ所滿足的Grad－Shafranov方程：

其中，（r，φ，z）表示柱坐標(biāo)；jφ（r，z）為大環(huán)方向電流密度。

方程（1）所表示的平衡中，磁面就是ψ＝const表示的面，磁面也是平衡位形的等壓面，且j只有在磁面上的分量［6］，由此可得：

其中，Ψ和F分別為相對于對稱軸r＝0的角向磁通和電流通量。

軸對稱等離子體的MHD平衡由準(zhǔn)線性橢圓方程（2）和方程（3）來決定。一般情況下，給定P（ψ）和f（ψ）的具體函數(shù)形式，加上合適的邊界條件，就可以解出ψ＝ψ（R，Z），已知ψ后，就可以給出平衡位形的所有物理量［7］。

在給定P（ψ）和q（ψ）求解平衡問題時(shí)，考慮角向磁通的守恒條件為：

其中，ψm為磁軸處的角向磁通；ψb為邊界處的角向磁通。

在磁面上安全因子的分布為：

其中，Φ為穿過ψ＝const截面的環(huán)向磁通。

由f（ψ）和q（ψ）的定義，可得積分區(qū)域?yàn)棣祝絚onst所形成的環(huán)面，即

由（2）式、（3）式、（6）式可見，對于所考慮的平衡問題，只要知道P（ψ）和q（ψ）就可以了。引入對相鄰磁面間體積V（ψ）的平均：

可以得到：

將（2）式寫為非線性積分微分方程的形式，即

對（9）式取平均后，可得：

在給定了（4）式中的后，就可以得到方程（9）的唯一解［8］。

在程序中，一般用表示歸一化的ψ，即

其中，ψm表示磁軸處的角向磁通；ψb表示邊界處的角向磁通。為方便表示，后面也將ψn簡記為ψ。另外，也用（12）式表示（3）式中的f［9］，即

2 TEQ讀取g－file并進(jìn)行計(jì)算

TEQ本身可以使用自由邊界和固定邊界2種方法求解等離子體平衡方程，但是由于這2種方法的計(jì)算原理不同，它們使用的變量也不相同。本文首先使用自由邊界開始計(jì)算，需要注意的是，使用自由邊界只是使TEQ程序開始運(yùn)行的一種手段。用自由邊界設(shè)定計(jì)算區(qū)域后，將自由邊界轉(zhuǎn)變?yōu)楣潭ㄟ吔?，再逐個(gè)改變需要改變的物理量，包括邊界的位置、安全因子分布、壓強(qiáng)分布以及控制程序運(yùn)行的參數(shù)等，輸出數(shù)據(jù)后再與EFIT所得結(jié)果進(jìn)行對比?；镜牧鞒倘鐖D1所示。

在使用固定邊界求解Grad－Shafranov方程時(shí)，選取常用的q和P的位形作為計(jì)算平衡的輸入數(shù)據(jù)，這些數(shù)據(jù)都從由EFIT計(jì)算出的g－file中獲得。g－file中包含許多信息，但并不全部需要，一般只需要知道計(jì)算區(qū)域、固定邊界位置、真空時(shí)磁軸處的磁感應(yīng)強(qiáng)度、真空時(shí)磁軸的橫坐標(biāo)、等離子體電流、安全因子分布和壓強(qiáng)分布即可。

圖1 新的TEQ程序處理流程

通過EFIT－viewer可以得到平行方向電流的分布，它是衡量TEQ計(jì)算結(jié)果好壞的重要標(biāo)準(zhǔn)。此外，還需要對比磁面的分布，由于EFIT求解得到的是網(wǎng)格點(diǎn)上的磁通值，進(jìn)而可以得到其等值線，而TEQ使用固定邊界求解得到的是磁面的位置，因此在對比的時(shí)候只要兩者不出現(xiàn)相交的情況，即認(rèn)為計(jì)算是合理的。

在計(jì)算的過程中，需要注意的是TEQ使用的單位與EFIT使用的單位是不相同的，需要對g－file中的數(shù)據(jù)進(jìn)行相應(yīng)處理后才能進(jìn)行計(jì)算。

3 TEQ 輸出g－file

TEQ本身也可以產(chǎn)生g－file，但直接利用TEQ自帶的程序輸出的結(jié)果并不理想，輸出的F、磁面位形與EFIT的計(jì)算結(jié)果有差別。

經(jīng)過檢查程序，發(fā)現(xiàn)TEQ默認(rèn)以計(jì)算區(qū)域的中心的橫坐標(biāo)為真空時(shí)磁軸的橫坐標(biāo)，而一般情況下，這兩者可能并不相等。

例如采用的EFIT的計(jì)算區(qū)域，即左邊界為120cm，右邊界為260cm，則TEQ默認(rèn)的磁軸位置在190cm處，而實(shí)際上真空時(shí)磁軸位置（EAST）在185.8cm 處。根據(jù)（12）式對F進(jìn)行修正后發(fā)現(xiàn)，它可以與EFIT的計(jì)算結(jié)果符合較好，因此需要在程序中人為地指定真空時(shí)磁軸的橫坐標(biāo)，然后再進(jìn)行計(jì)算。

對于磁面位形的差別，經(jīng)過對TEQ程序的分析，發(fā)現(xiàn)這主要是由于自由邊界和固定邊界的計(jì)算方法不同引起的。

自由邊界計(jì)算時(shí)，會對整個(gè)計(jì)算區(qū)域進(jìn)行計(jì)算，當(dāng)把自由邊界轉(zhuǎn)化為固定邊界時(shí)，程序會根據(jù)自由邊界的計(jì)算結(jié)果來確定一組磁面的位置，而之后再以固定邊界計(jì)算時(shí)，程序都不會再計(jì)算最外閉合磁面以外的區(qū)域，即TEQ自帶的輸出g－file的程序輸出的其實(shí)是最后一次自由邊界計(jì)算時(shí)的磁面位形。固定邊界是對磁面位置進(jìn)行計(jì)算，而自由邊界是對網(wǎng)格點(diǎn)上的磁通進(jìn)行計(jì)算。

因此，固定邊界的計(jì)算結(jié)果沒有辦法以網(wǎng)格的形式直接輸出。但用固定邊界可以直接輸出磁面的位置以及磁面上的磁通值，這相當(dāng)于一組等值線，而輸出g－file需要的則是這樣一種網(wǎng)格，利用網(wǎng)格上的點(diǎn)可以畫出與固定邊界所給出的磁面相同的磁面位形，因此，設(shè)計(jì)了用磁面位置及磁面上磁通值來給出相應(yīng)網(wǎng)格點(diǎn)上磁通的程序。

用固定邊界求解時(shí)，TEQ輸出的磁面上點(diǎn)的坐標(biāo)是很有規(guī)律的，當(dāng)把這些點(diǎn)的坐標(biāo)變換到以磁軸為極點(diǎn)的極坐標(biāo)時(shí)，發(fā)現(xiàn)它們對應(yīng)的極角的分布都是固定的，相同極角的點(diǎn)從磁軸處向外輻射至邊界處，這樣的規(guī)律給計(jì)算帶來了較大的方便。

可以先把磁面上點(diǎn)的坐標(biāo)從直角坐標(biāo)轉(zhuǎn)變?yōu)闃O坐標(biāo)，如果要計(jì)算一個(gè)網(wǎng)格點(diǎn)（X，Y）上的磁通值，先把它變成極坐標(biāo)（ρ，θ），然后找出該點(diǎn)鄰近的2條輻射線L1和L2，其中L1上的點(diǎn)對應(yīng)的角為θ1，L2上的點(diǎn)所對應(yīng)的角為θ2，θ1和θ2滿足θ1≤θ＜θ2。

在L1上找出P1（ρ1，θ1）、P2（ρ2，θ1），使ρ1≤ρ＜ρ2，再在L2上找出P3（ρ3，θ2）、P4（ρ4，θ2），使ρ3≤ρ＜ρ4，依次在L1和L2上進(jìn)行樣條插值得出ρ處的磁通值，然后再根據(jù)該網(wǎng)格點(diǎn)到L1和L2的距離進(jìn)行對距離的加權(quán)平均，這樣就可以得到其磁通值了。對于邊界處的點(diǎn)，先判斷該點(diǎn)在邊界內(nèi)還是邊界外，在邊界內(nèi)就用上面的方法來計(jì)算。邊界外的點(diǎn)，由于不能確定其磁通值具體是多少，而且在計(jì)算中也不需要知道邊界外的網(wǎng)格點(diǎn)的磁通值具體是多少，因此，統(tǒng)一設(shè)為比邊界處的磁通值略小的值即可。

算法流程如圖2所示，用此方法算出的gfile畫最外閉合磁面時(shí)，會出現(xiàn)不光滑的現(xiàn)象，但對內(nèi)部的閉合磁面基本沒有影響，可以通過增加網(wǎng)格點(diǎn)的數(shù)量來降低最外磁面不光滑的程度。

圖2 修正磁面位形的程序流程

4 計(jì)算結(jié)果對比

本文以EFIT計(jì)算出的剪切及反剪切情況的g－file作為數(shù)據(jù)來源。以反剪切情況為例，安全因子分布如圖3中虛＊線所示，其中，橫坐標(biāo)為歸一化的ψ。

圖3 反剪切情況下的安全因子分布

壓強(qiáng)分布如圖4中虛＊線所示，此種情況下，用EFIT－viewer畫出的平行于磁場方向的電流，如圖5中虛＊線所示，等離子體電流為1.05MA。在TEQ中用這樣的安全因子和壓強(qiáng)分布進(jìn)行計(jì)算，得到的安全因子分布如圖3中實(shí)線所示，壓強(qiáng)分布如圖4中實(shí)線所示，平行電流分布如圖5中實(shí)線所示，等離子體電流為1.05MA，磁面位形如圖6所示，其中虛＊線表示求解固定邊界時(shí)設(shè)定的邊界。由于TEQ采用固定邊界求解，在邊界以外的區(qū)域沒有進(jìn)行求解，因此，也就無法直接得到邊界以外的磁面位形?？梢钥吹椒醇羟星闆r下，由EFIT和TEQ得到的安全因子、壓強(qiáng)分布符合較好，平行電流分布也基本相同。通過比較，磁面位形也相同。另外，對剪切情況進(jìn)行計(jì)算，可以得出相同的結(jié)論。

圖4 反剪切情況下的壓強(qiáng)分布

圖5 反剪切情況下平行電流分布

圖6 反剪切情況下TEQ計(jì)算出的磁面位形

5 構(gòu)建反剪切平衡位形

TEQ可以直接利用安全因子和壓強(qiáng)來進(jìn)行平衡的構(gòu)建，這對設(shè)計(jì)新的托卡馬克具有較大意義，可以方便地調(diào)節(jié)安全因子分布和壓強(qiáng)分布，從而獲得想要的磁面位形，比用EFIT方便。由于反剪切可以使等離子體性能大大改善，因此對先進(jìn)托卡馬克而言，反剪切是一種吸引人的選擇，所以用圖3中的安全因子分布、圖4中壓強(qiáng)分布分別乘以因子 0.50、0.75、1.00、1.25、1.50、1.75后，得到的一組壓強(qiáng)分布作為輸入來構(gòu)建固定邊界平衡，得到6組不同平衡參數(shù)，見表1所列［10－11］。以這種方式計(jì)算大量數(shù)據(jù)后，可以為構(gòu)建需要的平衡提供一種參考依據(jù)。

表1 6組反剪切平衡參數(shù)

由表1可以看出，隨著壓強(qiáng)的變大，磁軸橫坐標(biāo)變大，表明壓強(qiáng)變大時(shí)，磁軸向外側(cè)移動，而且壓強(qiáng)越大，磁軸移動得越多。

這可以定性解釋為：在托卡馬克中，磁力線的平衡由指向內(nèi)側(cè)的磁張力和指向外側(cè)的磁壓力、向外的壓力（由氣體壓強(qiáng)梯度產(chǎn)生）等來提供。在計(jì)算中，未改變磁場的大小，因此，磁壓力的影響可以忽略，而當(dāng)壓強(qiáng)等比例增大時(shí)，壓強(qiáng)梯度也變大，這樣向外的壓力就會增大，這只有依靠變大的磁張力來產(chǎn)生，而這就要求磁力線發(fā)生更大的形變來提供，于是便產(chǎn)生了上面的結(jié)果。

6 結(jié)束語

本文對TEQ輸出g－file的程序進(jìn)行了修正，使得計(jì)算出的安全因子、壓強(qiáng)、F、平行方向電流的分布以及磁面位形都能與EFIT符合得非常好。通過對剪切和反剪切2種情況的檢驗(yàn)，可以認(rèn)為在TEQ基礎(chǔ)上增加的讀取g－file，并進(jìn)行固定邊界求解的功能得到了比較好的實(shí)現(xiàn)。

利用TEQ，可以方便地在知道安全因子分布和壓強(qiáng)分布的情況下，進(jìn)行平衡的構(gòu)建，這對于設(shè)計(jì)期望的磁面位形具有重要意義。另外，從研究壓強(qiáng)對磁軸水平位置的影響發(fā)現(xiàn)，當(dāng)壓強(qiáng)變大時(shí)，磁軸向外移動，并且壓強(qiáng)變得越大，磁軸移動得越多。

［1］Kritz A H，Bateman G，Kinsey J，et al.The national transport code collaboration module library［J］.Computer Physics Communications，2004，164（1／3）：108－113.

［2］Crotinger J A，Lo Destro L，Pearlstein L D，et al.Corsica：a comprehensive simulation of toroidal magnetic－fusion devices［EB／OL］.（1997－06－12）［2011－10－19］.http：／／wormhole.ucllnl.org／software／corsica／finalreport.pdf.

［3］Carlsson J，Tech－X Corporation.TEQ library user manual［EB／OL］.（2007－07－31）［2010－10－15］.http：／／w3.pppl.gov／ntcc／Teq／TEQmanual.pdf.

［4］Carlsson J，Tech－X Corporation.Separatrix library user manual［EB／OL］.（2006－01－14）［2010－10－18］.http：／／w3.pppl.gov／ntcc／Separatrix／SeparatrixManual.pdf.

［5］Ferron J R，Walker M L，Lao L L，et al.Real time equilibrium reconstruction for tokamak discharge control［J］.Nuclear Fusion，1998，38（7）：1055－1066.

［6］胡希偉.等離子體理論基礎(chǔ) ［M］.北京：北京大學(xué)出版社，2006：23－69.

［7］查學(xué)軍，朱思錚，虞清泉.托卡馬克等離子體平衡的數(shù)值研究 ［J］.計(jì)算物理，2002，19（5）：413－418.

［8］Degtyarev L M，Drozdov V V.An inverse variable technique in the MHD－equilibrium problem ［J］.Computer Physics Reports，1985，2（7）：341－387.

［9］Bulmer R H，Pearlstein L D.DIII－D equilibrium and stability modeling with Caltrans［EB／OL］.（2003－01－28）［2010－10－15］.http：／／wormhole.ucllnl.org／software／corsica／d3d－demo.pdf.

［10］邱勵儉.聚變能及其應(yīng)用 ［M］.北京：科學(xué)出版社，2008：152－153.

［11］程曹軍，肖炳甲，劉成岳.基于Corsica的EAST等離子體平衡計(jì)算 ［J］.合肥工業(yè)大學(xué)學(xué)報(bào)：自然科學(xué)版，2011，34（4）：613－616.

Simulation of fixed－boundary tokamak plasma equilibrium with TEQ

LIU Zhi－yuan， LI Guo－qiang， XU Li－qing
（Institute of Plasma Physics，Chinese Academy of Sciences，Hefei 230031，China）

A new calculation process of TEQ is designed.An algorithm to calculate magnetic fluxes on grid points with positions of magnetic surfaces and corresponding values of the magnetic fluxes is developed.And the program of outputting g－file is improved for solving the fixed－boundary tokamak plasma equilibrium with TEQ.Based on the EAST device，two groups of safety factor and pressure profile which denote shear and reverse shear respectively are used to test the program.It is shown that the result of TEQ coincides with that of EFIT very well.Finally，six reverse shear fixed－boundary equilibriums are constructed and the effect of pressure on the horizontal coordinate of magnetic axis is analyzed.

TEQ program；g－file program；EAST device；fixed－boundary plasma equilibrium；numerical simulation

O539

A

1003－5060（2012）03－0412－05

10.3969／j.issn.1003－5060.2012.03.028

2011－09－14；

2011－10－25

國家自然科學(xué)基金資助項(xiàng)目（10975161）；科技部ITER專項(xiàng)課題資助項(xiàng)目（2009GB01001）

劉志遠(yuǎn)（1988－），男，山東濟(jì)寧人，中國科學(xué)院等離子體物理研究所碩士生；

李國強(qiáng)（1977－），男，河南杞縣人，博士，中國科學(xué)院等離子體物理研究所副研究員，碩士生導(dǎo)師.

（責(zé)任編輯 呂 杰）