托卡马克聚变等离子体参数演化的数值模拟

王时佳，王少杰

(1.中国电子科技集团公司第20研究所 雷达事业部, 陕西 西安 710068;2.中国科学技术大学 近代物理系，安徽 合肥 230026)



托卡马克聚变等离子体参数演化的数值模拟

王时佳1，王少杰2

(1.中国电子科技集团公司第20研究所 雷达事业部, 陕西 西安 710068;2.中国科学技术大学 近代物理系，安徽 合肥 230026)

基于一维粒子-能量输运方程组，采用差分方法编写了模拟托卡马克聚变等离子体的密度和温度剖面演化的数值求解程序。在其中使用预估-校正方法处理了方程中的非线性项并利用能量约束时间定标律求得输运系数。模拟得到的聚变性能与之前文献采用其他方法计算的结果符合较好。该程序可用于研究托卡马克反应堆的聚变性能以及温度和密度的剖面形状。

数值模拟；预估校正算法；非线性方程；定标律

国际热核聚变实验反应堆(ITER)将成为人类历史上第一个能长时间维持聚变增益因子Q(聚变释放功率与装置消耗的辅助加热功率之比)大于10的聚变实验装置[1-3]。为了达到这样的高性能状态，ITER中的反应燃料必须成为燃烧等离子体状态，即装置中的等离子体主要由聚变反应所释放的能量进行加热。研究表明聚变反应的发生率与燃料的密度呈非线性正相关关系，而与温度呈非线性非单调关系[4]。因此在燃烧等离子体中，温度的演化将是高度非线性的。

作为一种托卡马克类型的磁约束聚变实验装置，ITER的主要设计参数是由现有的托卡马克装置实验运行参数进行统计和外推而得到的[2,5]。ITER 要实现设计目标中的高参数运行仍具有较大的不确定性和挑战性。因此，在实验装置建造之前进行可靠的数值模拟验证工作仍具有重要的价值。目前主流的关于ITER 聚变性能计算的模拟工作主要使用粒子输运和能量输运相互独立的模型来进行[3]，对于密度和温度演化的自洽性仍有欠缺。

本文从一维柱位形下的等离子体流体输运方程出发，使用能量约束定标律建立了一种约化的等离子体输运模型。其中利用预估-校正算法处理了能量输运方程中的非线性源项，模拟了等离子体密度剖面和温度剖面的自洽演化，计算了反应堆的聚变性能。

1 输运方程及其性质

基于电中性等离子体条件以及等温等离子体假设(电子和离子具有相同的温度T)，燃烧等离子体的输运过程可由以下一维柱坐标输运方程组描述[6]

(1)

(2)

(3)

(4)

其中，n是电子密度；qr是热传导导致的径向热流密度；Γr是径向粒子流密度；ur是粒子的径向对流速度；D和χ分别是粒子与热量的输运系数；Sp代表粒子源，Sα和Saux分别是氘氚聚变反应和辅助加热设备带来的能量源项；Srad代表因辐射带走能量而导致的能量损失项；r代表等离子体小半径方向上的位置。

与传统的能量输运方程不同[7-8]，燃烧等离子体中的物理过程是高度非线性的：一方面等离子体主要由核聚变反应生成的阿尔法粒子碰撞加热，而局域内的阿尔法粒子生成率是密度以及温度的非线性函数；另一方面在等离子体能量损失中占据相当比重的辐射损失的单位体积发射强度也是温度和密度的非线性函数。具体表现在输运方程中，有

Sα=kn2<σv>Eα

(5)

Srad=SBrm+SSyn+SL

(6)

SBrm=CBn2T1/2

(7)

(8)

sl=clNt

(9)

其中，k是由氘氚燃料密度与电子密度相对比例决定的系数；<σv>是氘氚反应的碰撞截面，根据文献[4]，<σv>=3.68×10-12T-2/3exp(-19.94T-1/3)cm3sec-1(T的单位为keV)；Eα=3.5 MeV是聚变产生的阿尔法粒子携带的能量；SBrm，SSyn和SL分别代表韧致辐射，回旋辐射和线辐射带来的能量损失；CB，CS，CL和A是由装置参数决定的常数[9]。

2 数值方法

3 能量定标律与输运系数

本文使用了一种约化方法来求出输运系数：在任意时刻根据ITER能量约束时间定标律，由系统的各宏观物理量计算出此时的能量约束时间并据此计算出下一时刻的系统中等离子体的总能量；固定输运系数的剖面形状，而对于其绝对大小数值做猜测值并带入输运方程组，计算出相应的下一时刻等离子体总能量，将其同根据能量约束时间计算出的数值进行比较并据此修正输运系数大小的猜测值，再将新的输运系数大小猜测值带入输运方程组并迭代计算猜测值，直到输运系数对应的下一时刻等离子体总能量与根据定标律计算得到的结果相符，即得到欲求得的输运系数。具体来说，假设热输运系数有以下形式

χ=χ0g(r)

(10)

其中，g(r)为给定的输运系数剖面形状因子；χ0为待求数值。χ0在tm时刻到tm+1时刻演化时所用的数值通过迭代求解数值方程来确定，符合要求的χ0可使得tm时刻和tm+1时刻之间等离子体总能量W的变化满足

(11)

(12)

本文设定粒子输运系数D与χ有以下关系

D=0.6χ

(13)

这是一个被广泛采用的近似，并被许多实验验证过[10-11]。

4 计算实例

根据上述数值方法，这里给出在ITER典型长脉冲放电参数下，对于燃烧等离子体演化过程的模拟过程和结果。

图1 输运系数形状

对于输运系数的剖面形状，这里假设(10)式中的输运系数形状因子g(r)有如图1中所示的剖面，注意到图中的输运系数形状具有典型的外部输运垒(ETB)特征，输运垒的宽度为0.05 个等离子体小半径。

边界条件设置为

(14)

(15)

在设置了一系列初值条件测试后我们发现，最后的稳态剖面对于初值并不敏感，这与预期的由燃烧等离子体自主演化决定剖面的性质相符。这里设置温度初始分布T0(r)和密度分布初值n0(r)如图2所示。

图2 初始温度和密度剖面

对于能量定标律中涉及到的装置宏观物理量，这里取值与文献[5]相同，其具体取值参见表1。对于径向粒子流密度中包含的粒子对流速度ur，假设为粒子箍缩速度，根据文献[13],其大小和分布为

(16)

其中，q=q(r)为托卡马克的安全因子，这里取q为ITER感应放电运行下的典型形状[12]，其大小和形状如图6所示。

表1 能量约束定标律所需参数

由以上条件，计算得到了稳定的等离子体各参数剖面，其形状如图3～图7所示。由稳定剖面得到的聚变增益因子Q=10.84，与之前文献中通过其他模型计算得到的ITER性能基本相同[5,12,14-15]。

图3 稳态密度剖面

图4 稳态温度剖面

图5 稳态压强剖面

图6 稳态箍缩速度及安全因子剖面

图7 稳态辐射功率密度剖面

5 结束语

对于托卡马克聚变装置来说，其聚变性能主要依赖于密度剖面和温度剖面的维持。本文采用数值求解等离子体粒子和能量输运方程组的方法模拟了ITER参数下托卡马克聚变反应堆的聚变性能。其中考虑到聚变反应加热源项的非线性，对其采用了预估-校正方法处理。模拟中需要使用的输运系数通过能量约束定标律进行了计算。模拟得到的结果与之前文献用其他方法得到的结果符合较好，可以用来评价聚变反应堆性能。此外，该结果还可以用于进一步研究加料和辅助加热对于聚变反应堆性能的影响。

[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] 熊波,杨威.油罐温度场的数值计算[J].电子科技,2010,23(6):91-93.

[8] 程超. ITER屏蔽块的热流体和热应力计算[D].成都:电子科技大学,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].现代电力,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

国家自然科学基金(11175178)

王时佳(1986-)，男，博士，助理工程师。研究方向：数值分析及数据处理。

10.16180/j.cnki.issn1007-7820.2017.06.002

TN01;O241.82

A

1007-7820(2017)06-005-04