王時佳，王少杰

(1.中國電子科技集團公司第20研究所 雷達事業(yè)部, 陜西 西安 710068;2.中國科學(xué)技術(shù)大學(xué) 近代物理系，安徽 合肥 230026)

?

托卡馬克聚變等離子體參數(shù)演化的數(shù)值模擬

王時佳1，王少杰2

(1.中國電子科技集團公司第20研究所 雷達事業(yè)部, 陜西 西安 710068;2.中國科學(xué)技術(shù)大學(xué) 近代物理系，安徽 合肥 230026)

基于一維粒子-能量輸運方程組，采用差分方法編寫了模擬托卡馬克聚變等離子體的密度和溫度剖面演化的數(shù)值求解程序。在其中使用預(yù)估-校正方法處理了方程中的非線性項并利用能量約束時間定標律求得輸運系數(shù)。模擬得到的聚變性能與之前文獻采用其他方法計算的結(jié)果符合較好。該程序可用于研究托卡馬克反應(yīng)堆的聚變性能以及溫度和密度的剖面形狀。

數(shù)值模擬；預(yù)估校正算法；非線性方程；定標律

國際熱核聚變實驗反應(yīng)堆(ITER)將成為人類歷史上第一個能長時間維持聚變增益因子Q(聚變釋放功率與裝置消耗的輔助加熱功率之比)大于10的聚變實驗裝置[1-3]。為了達到這樣的高性能狀態(tài)，ITER中的反應(yīng)燃料必須成為燃燒等離子體狀態(tài)，即裝置中的等離子體主要由聚變反應(yīng)所釋放的能量進行加熱。研究表明聚變反應(yīng)的發(fā)生率與燃料的密度呈非線性正相關(guān)關(guān)系，而與溫度呈非線性非單調(diào)關(guān)系[4]。因此在燃燒等離子體中，溫度的演化將是高度非線性的。

作為一種托卡馬克類型的磁約束聚變實驗裝置，ITER的主要設(shè)計參數(shù)是由現(xiàn)有的托卡馬克裝置實驗運行參數(shù)進行統(tǒng)計和外推而得到的[2,5]。ITER 要實現(xiàn)設(shè)計目標中的高參數(shù)運行仍具有較大的不確定性和挑戰(zhàn)性。因此，在實驗裝置建造之前進行可靠的數(shù)值模擬驗證工作仍具有重要的價值。目前主流的關(guān)于ITER 聚變性能計算的模擬工作主要使用粒子輸運和能量輸運相互獨立的模型來進行[3]，對于密度和溫度演化的自洽性仍有欠缺。

本文從一維柱位形下的等離子體流體輸運方程出發(fā)，使用能量約束定標律建立了一種約化的等離子體輸運模型。其中利用預(yù)估-校正算法處理了能量輸運方程中的非線性源項，模擬了等離子體密度剖面和溫度剖面的自洽演化，計算了反應(yīng)堆的聚變性能。

1 輸運方程及其性質(zhì)

基于電中性等離子體條件以及等溫等離子體假設(shè)(電子和離子具有相同的溫度T)，燃燒等離子體的輸運過程可由以下一維柱坐標輸運方程組描述[6]

(1)

(2)

(3)

(4)

其中，n是電子密度；qr是熱傳導(dǎo)導(dǎo)致的徑向熱流密度；Γr是徑向粒子流密度；ur是粒子的徑向?qū)α魉俣?；D和χ分別是粒子與熱量的輸運系數(shù)；Sp代表粒子源，Sα和Saux分別是氘氚聚變反應(yīng)和輔助加熱設(shè)備帶來的能量源項；Srad代表因輻射帶走能量而導(dǎo)致的能量損失項；r代表等離子體小半徑方向上的位置。

與傳統(tǒng)的能量輸運方程不同[7-8]，燃燒等離子體中的物理過程是高度非線性的：一方面等離子體主要由核聚變反應(yīng)生成的阿爾法粒子碰撞加熱，而局域內(nèi)的阿爾法粒子生成率是密度以及溫度的非線性函數(shù)；另一方面在等離子體能量損失中占據(jù)相當比重的輻射損失的單位體積發(fā)射強度也是溫度和密度的非線性函數(shù)。具體表現(xiàn)在輸運方程中，有

Sα=kn2<σv>Eα

(5)

Srad=SBrm+SSyn+SL

(6)

SBrm=CBn2T1/2

(7)

(8)

sl=clNt

(9)

其中，k是由氘氚燃料密度與電子密度相對比例決定的系數(shù)；<σv>是氘氚反應(yīng)的碰撞截面，根據(jù)文獻[4]，<σv>=3.68×10-12T-2/3exp(-19.94T-1/3)cm3sec-1(T的單位為keV)；Eα=3.5 MeV是聚變產(chǎn)生的阿爾法粒子攜帶的能量；SBrm，SSyn和SL分別代表韌致輻射，回旋輻射和線輻射帶來的能量損失；CB，CS，CL和A是由裝置參數(shù)決定的常數(shù)[9]。

2 數(shù)值方法

3 能量定標律與輸運系數(shù)

本文使用了一種約化方法來求出輸運系數(shù)：在任意時刻根據(jù)ITER能量約束時間定標律，由系統(tǒng)的各宏觀物理量計算出此時的能量約束時間并據(jù)此計算出下一時刻的系統(tǒng)中等離子體的總能量；固定輸運系數(shù)的剖面形狀，而對于其絕對大小數(shù)值做猜測值并帶入輸運方程組，計算出相應(yīng)的下一時刻等離子體總能量，將其同根據(jù)能量約束時間計算出的數(shù)值進行比較并據(jù)此修正輸運系數(shù)大小的猜測值，再將新的輸運系數(shù)大小猜測值帶入輸運方程組并迭代計算猜測值，直到輸運系數(shù)對應(yīng)的下一時刻等離子體總能量與根據(jù)定標律計算得到的結(jié)果相符，即得到欲求得的輸運系數(shù)。具體來說，假設(shè)熱輸運系數(shù)有以下形式

χ=χ0g(r)

(10)

其中，g(r)為給定的輸運系數(shù)剖面形狀因子；χ0為待求數(shù)值。χ0在tm時刻到tm+1時刻演化時所用的數(shù)值通過迭代求解數(shù)值方程來確定，符合要求的χ0可使得tm時刻和tm+1時刻之間等離子體總能量W的變化滿足

(11)

(12)

本文設(shè)定粒子輸運系數(shù)D與χ有以下關(guān)系

D=0.6χ

(13)

這是一個被廣泛采用的近似，并被許多實驗驗證過[10-11]。

4 計算實例

根據(jù)上述數(shù)值方法，這里給出在ITER典型長脈沖放電參數(shù)下，對于燃燒等離子體演化過程的模擬過程和結(jié)果。

圖1 輸運系數(shù)形狀

對于輸運系數(shù)的剖面形狀，這里假設(shè)(10)式中的輸運系數(shù)形狀因子g(r)有如圖1中所示的剖面，注意到圖中的輸運系數(shù)形狀具有典型的外部輸運壘(ETB)特征，輸運壘的寬度為0.05 個等離子體小半徑。

邊界條件設(shè)置為

(14)

(15)

在設(shè)置了一系列初值條件測試后我們發(fā)現(xiàn)，最后的穩(wěn)態(tài)剖面對于初值并不敏感，這與預(yù)期的由燃燒等離子體自主演化決定剖面的性質(zhì)相符。這里設(shè)置溫度初始分布T0(r)和密度分布初值n0(r)如圖2所示。

圖2 初始溫度和密度剖面

對于能量定標律中涉及到的裝置宏觀物理量，這里取值與文獻[5]相同，其具體取值參見表1。對于徑向粒子流密度中包含的粒子對流速度ur，假設(shè)為粒子箍縮速度，根據(jù)文獻[13],其大小和分布為

(16)

其中，q=q(r)為托卡馬克的安全因子，這里取q為ITER感應(yīng)放電運行下的典型形狀[12]，其大小和形狀如圖6所示。

表1 能量約束定標律所需參數(shù)

由以上條件，計算得到了穩(wěn)定的等離子體各參數(shù)剖面，其形狀如圖3～圖7所示。由穩(wěn)定剖面得到的聚變增益因子Q=10.84，與之前文獻中通過其他模型計算得到的ITER性能基本相同[5,12,14-15]。

圖3 穩(wěn)態(tài)密度剖面

圖4 穩(wěn)態(tài)溫度剖面

圖5 穩(wěn)態(tài)壓強剖面

圖6 穩(wěn)態(tài)箍縮速度及安全因子剖面

圖7 穩(wěn)態(tài)輻射功率密度剖面

5 結(jié)束語

對于托卡馬克聚變裝置來說，其聚變性能主要依賴于密度剖面和溫度剖面的維持。本文采用數(shù)值求解等離子體粒子和能量輸運方程組的方法模擬了ITER參數(shù)下托卡馬克聚變反應(yīng)堆的聚變性能。其中考慮到聚變反應(yīng)加熱源項的非線性，對其采用了預(yù)估-校正方法處理。模擬中需要使用的輸運系數(shù)通過能量約束定標律進行了計算。模擬得到的結(jié)果與之前文獻用其他方法得到的結(jié)果符合較好，可以用來評價聚變反應(yīng)堆性能。此外，該結(jié)果還可以用于進一步研究加料和輔助加熱對于聚變反應(yīng)堆性能的影響。

[1] ITER Physics Expert Group.Chapter 1: Overview and summary[J].Nucl Fusion,1999,39(12):2137-2174.

[2] Shimada M,Campbell D J,Mukhovatov V,et al.Chapter 1:Overview and summary[J].Nucl Fusion, 2007,47(6):S1-S17.

[3] Sips A C C,Giruzzi G,Ide S,et al.Progress in preparing scenarios for operation of the international thermonuclear experimental reactor[J].Physics Plasmas,2015,22(2):021804-021827.

[4] Huba J D.NRL plasma formulary[M].NA,USA:Naval Research Laboratory,2009.

[5] Doyle E J,Houlberg W A,Kamada Y,et al.Chapter 2:Plasma confinement and transport[J].Nucl Fusion,2007,47(6):S18-S127.

[6] Cardozo N J L.Perturbative transport studies in fusion plasmas[J].Plasma Physics Controlled Fusion,1995,37(8):799-852.

[7] 熊波,楊威.油罐溫度場的數(shù)值計算[J].電子科技,2010,23(6):91-93.

[8] 程超. ITER屏蔽塊的熱流體和熱應(yīng)力計算[D].成都:電子科技大學(xué),2011.

[9] Uckan N A.ITER physics design guidelines[R].Vienna:IAEA,1990.

[10] Becker G.Study of anomalous inward drift in tokamaks by transport analysis and simulations [J].Nucl Fusion,2004,44(9):933-944.

[11] Becker G,Kardaun O.Anomalous particle pinch and scaling of vin/d based on transport analysis and multiple regression[J].Nucl Fusion,2007,47(1):33-43.

[12] Uckan N A.ITER technical basis [R].Vienna:IAEA,2002.

[13] Baker D R,Rosenbluth M N.Density profile consistency and its relation to the transport of trapped versus passing electrons in tokamaks[J].Physics Plasmas,1998,5(8):2936-2941.

[14] Pereverzev G V,Angioni C,Peeters A G,et al.Theoretical predictions of the density profile in a tokamak reactor[J].Nucl Fusion,2005,45(4):221-225.

[15] 馮開明.可控核聚變與ITER計劃[J].現(xiàn)代電力,2006(5):82-88.

Numerical Simulation on the Evolution of Plasma Parameters in Tokamak Fusion Reactor

WANG Shijia1,WANG Shaojie2

(1. Radar Department,20th Research Institute of China Electronics Technology Group,Xi’an 710068,China;2. Department of Modern Physics, University of Science and Technology of China,Hefei 230026,China)

The evolution of plasma temperature and density profiles is simulated with a code which solving the one-dimensional transport formulas with difference method. The nonlinear terms in the formulas are treated with a predictor-corrector method and the transport coefficient is evaluated by matching the scaling law of energy confinement time. The fusion performance calculated with the modeling results accord well with the results calculated with other models in the references. The code in this paper can be used to study the fusion performance and the profiles of temperature and density in tokamak fusion reactors.

numerical simulation; predictor-corrector method; non-linear formula; scaling law

2017- 03- 20

國家自然科學(xué)基金(11175178)

王時佳(1986-)，男，博士，助理工程師。研究方向：數(shù)值分析及數(shù)據(jù)處理。

10.16180/j.cnki.issn1007-7820.2017.06.002

TN01;O241.82

A

1007-7820(2017)06-005-04