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

王時(shí)佳，王少杰

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

?

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

王時(shí)佳1，王少杰2

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

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

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

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

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

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

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

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

(1)

(2)

(3)

(4)

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

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

Sα=kn2<σv>Eα

(5)

Srad=SBrm+SSyn+SL

(6)

SBrm=CBn2T1/2

(7)

(8)

sl=clNt

(9)

其中，k是由氘氚燃料密度與電子密度相對(duì)比例決定的系數(shù)；<σv>是氘氚反應(yīng)的碰撞截面，根據(jù)文獻(xiàn)[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分別代表韌致輻射，回旋輻射和線輻射帶來(lái)的能量損失；CB，CS，CL和A是由裝置參數(shù)決定的常數(shù)[9]。

2 數(shù)值方法

3 能量定標(biāo)律與輸運(yùn)系數(shù)

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

χ=χ0g(r)

(10)

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

(11)

(12)

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

D=0.6χ

(13)

這是一個(gè)被廣泛采用的近似，并被許多實(shí)驗(yàn)驗(yàn)證過(guò)[10-11]。

4 計(jì)算實(shí)例

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

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

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

邊界條件設(shè)置為

(14)

(15)

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

圖2 初始溫度和密度剖面

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

(16)

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

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

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

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

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

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

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

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

5 結(jié)束語(yǔ)

對(duì)于托卡馬克聚變裝置來(lái)說(shuō)，其聚變性能主要依賴于密度剖面和溫度剖面的維持。本文采用數(shù)值求解等離子體粒子和能量輸運(yùn)方程組的方法模擬了ITER參數(shù)下托卡馬克聚變反應(yīng)堆的聚變性能。其中考慮到聚變反應(yīng)加熱源項(xiàng)的非線性，對(duì)其采用了預(yù)估-校正方法處理。模擬中需要使用的輸運(yùn)系數(shù)通過(guò)能量約束定標(biāo)律進(jìn)行了計(jì)算。模擬得到的結(jié)果與之前文獻(xiàn)用其他方法得到的結(jié)果符合較好，可以用來(lái)評(píng)價(jià)聚變反應(yīng)堆性能。此外，該結(jié)果還可以用于進(jìn)一步研究加料和輔助加熱對(duì)于聚變反應(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] 熊波,楊威.油罐溫度場(chǎng)的數(shù)值計(jì)算[J].電子科技,2010,23(6):91-93.

[8] 程超. ITER屏蔽塊的熱流體和熱應(yīng)力計(jì)算[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] 馮開(kāi)明.可控核聚變與ITER計(jì)劃[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

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

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

10.16180/j.cnki.issn1007-7820.2017.06.002

TN01;O241.82

A

1007-7820(2017)06-005-04