基于蒙特卡罗-离散纵标方法的氘氚激光等离子体聚变反应率数值模拟*

陈忠 赵子甲 吕中良 李俊汉 潘冬梅

1) (西南科技大学国防科技学院,绵阳 621010)

2) (国防科技大学文理学院,长沙 410073)

3) (中国科学技术大学,合肥 230027)

惯性约束聚变(ICF)是实现受控热核聚变可能途径之一.聚变中子源项是氘氚激光等离子体物理设计与分析的重要参数之一,其准确性直接影响分析结果的可靠性.目前国内外对于ICF氘氚聚变反应产生的中子源项研究主要基于解析公式法,在温度和反应类型等方面适用范围有限.本文采用粒子云概念对氘、氚粒子云团开展了随机抽样与时空网格划分,然后基于麦克斯韦速率分布律对氘氚聚变反应开展了多普勒能量展宽效应分析与微分截面温度修正工作,耦合蒙特卡罗方法和离散纵标方法,开展了激光等离子体中D−T粒子云团聚变反应率的数值模拟工作.研究结果显示,与原核数据库截面相比,D−T,D−D,T−D截面经修正后多普勒温度效应显著.在20－100 keV的等离子体温度范围内,相较传统的解析公式法,本文模拟结果更符合最新的ENDF核数据库的氘氚反应截面数据,且与解析公式法结果在低能区存在较大误差,可能是计算方法不同与低温截面差异过大原因导致.

1 引 言

惯性约束聚变(ICF)是实现受控热核聚变可能途径之一[1,2].点火是研究惯性约束核聚变的关键,所获得的聚变中子源项是激光等离子体物理设计与分析的重要参数之一,其准确性直接影响分析结果的可靠性[3−14].目前国内外对于ICF氘氚聚变反应产生的中子源项研究主要基于解析公式法,如对于氘氚各占50%密度的情况,其产生率计算公式为

其中,S(ni,Ti)为中子产生率,ni为氘或者氚粒子密度,Ti为等离子体温度,σ为反应截面,vDT为麦克斯韦分布下的氘氚粒子相对速度,表示两者的加权平均[15].

目前大量使用的解析公式法主要针对氘氚粒子密度相同且温度不大于100 keV,在温度和反应类型等方面适用范围有限.本文通过结合蒙特卡罗(MC)方法和离散纵标(SN)方法,开展了激光等离子体中D−T粒子云团聚变反应率的数值模拟算法程序设计.

2 算法设计

2.1 物理模型与碰撞物理过程

对于所给定的氘氚等离子体集团,在进行模拟时,由于数据量过于庞大,不可能把每个粒子都进行聚变反应的模拟,因此本文首先引入了粒子云概念,即多个位置、速度相同的同类粒子的集合,且认为粒子云所发生的作用均为同一反应.其次使用了网格模型,即将所模拟的时空进行了网格划分,其中整个等离子体空间经网格化后构成一个三维的空间坐标,一个空间网格中可包含多个粒子云.

氘氚等离子体为大量粒子组成,满足麦克斯韦速率分布律,即

其中,m为粒子质量,k为玻尔兹曼常数,T为等离子团温度.

本文选取了ENDF/B−VI和ENDF/B−VII的氘氚聚变反应数据库,其所给出的粒子能量为系统的质心系能量,因此需要将实验室系下入射粒子与靶粒子碰撞时的对应能量转换为质心系能量.等离子体中,氘、氚粒子由于速度和质量不同,用动量来进行矢量合成.根据动量守恒

图1 动量矢量分解Fig.1.Momentum vector decomposition.

将(3)式写成标量形式：

其中,ϕ为两矢量夹角,Px和Py分别为质心系动量的x轴和y轴分量.ϕ角即为需要离散化处理的角度分布.

所模拟的物理过程为： 某一随机位置处,具有某一随机速率的入射粒子云沿某一随机方向发射,则根据核反应率的计算原理,沿途路径上的所有网格,均认为入射粒子云与其中所有粒子发生反应,且按照截面大小发生核反应.由此,在一个时间网格内,氘氚聚变反应的中子产生率为

其中,S为中子产生率;ni为入射粒子的数密度,nt为靶粒子的数密度,当氘氚粒子各占50%时,ni=nt;σ′为修正后的氘氚聚变反应截面;t为模拟的时间步长;l为入射粒子云在一个时间网格内走过的路程.

2.2 时空网格划分与误差控制

氘氚等离子体发生聚变,并实现可持续惯性热核聚变燃烧,必须满足以下3个基本条件：

1)劳森判据条件niτ≥1014s/cm3;

2)燃料等离子体温度条件Th≥5×107K ;

3)ρr乘积条件ρmr≥3g/cm2;

其中,ni为热核燃料等离子体密度,τ为由惯性维持该离子数密度不变的时间间隔,ρm为被压缩的低温热核燃料质密度,r为预压缩到高质密度热核燃料小球的半径[16].

为减少氘氚聚变反应率计算误差及计算量,激光等离子体的时空网格划分极为重要.从误差产生源来看,假设入射粒子云所经过的网格,均要计算其内部的粒子云与入射粒子云的反应,但如果入射粒子云只经过网格边缘部分,则其内部大部分粒子云并没有与入射粒子云发生碰撞,由此导致误差.以二维网格为例,如图2所示.

图2 二维网格误差分析Fig.2.Two−dimensional grid error analysis.

图2中的灰色网格代表着仅其边缘部分被入射粒子云穿过,可以看出,若将该类网格均计入氘氚聚变反应,模拟结果将由此偏大.同样,入射粒子云的初始点和终点位置也会造成误差,如图3所示.

图3 二维网格大小造成的误差Fig.3.Error caused by two−dimensional mesh size.

图3给出了两个入射粒子云,由于终点位置的不同,导致实际路径长度差别接近一个网格.因此,网格划分越细,则误差越小.而网格大小最主要的判定条件取决于入射粒子所走的距离长短,若网格过大,每个时间网格(步长)所通过的路程则有可能在一两个网格内,由此会导致巨大的误差.根据上述的误差来源,需要考虑两个条件： 一是时间网格大小即时间步长,二是氘氚等离子体热运动的速率.根据公式

计算所走路程L,以此推出合适的网格大小设置.

对于惯性约束核聚变装置,如美国国家点火装置,其实际发生核反应的时间大约在10－10－10－11s.根据流体力学中判断计算收敛的柯朗−弗里德里奇−列维条件,计算并设定时间步长为t=1 ns.

对于氘氚等离子体热运动速率,应以模拟中的最小速率为下限,选择模拟所设定的最小温度的最概然速率作为估算网格速率的参考速率值,则有

其中,vm为T温度下麦克斯韦速率分布的最概然速率.根据(9)式计算,等离子体温度为20 keV时,氘粒子入射,其最概然速率为1.9575 × 106m/s;氚粒子入射,其最概然速率为1.1325 × 106m/s.估算时,只考虑这个速率的数量级106,由(8)式可得

即,入射粒子云所走的路径长度在10－3m量级以上,为保证10%以下的误差,网格大小应至少在10－4m以下,故设定网格大小为10－4m.

2.3 不同反应类型截面修正

本文所用数据来自ENDF/B−VI和ENDF/B−VII的氘氚聚变反应数据库,其数据由美国洛斯阿拉莫斯国家实验室(LANL)利用EDA−R矩阵码对氘氚聚变反应进行研究所生成的评价中子数据.其中,D−T的反应数据是ENDF/B−VI和ENDF/B−VII的官方原数据,D−D的反应数据是洛斯阿拉莫斯国家实验室的一个初步结果,有ENDF格式、HTML格式和PDF格式三种[17].

本文模拟的是氘氚等离子体聚变反应的中子产生率,只需要考虑3个反应,分别是

对于给定温度,氘氚等离子体在满足麦克斯韦速率分布的条件下,在任何速度、任何角度的位置上都存在粒子,因此在进行截面修正时,应同时考虑这两个物理量对微分截面的影响.本文修正截面的目的在于解决大数据量模拟计算时计算空间不足的问题,将所有满足同一麦克斯韦速率分布的粒子的截面归一化为一个截面值,即当粒子数量足够多时,无论入射粒子与靶粒子的夹角有多大,靶粒子的速率大小有多大,可使用同一微分截面值.

按照麦克斯韦速率分布函数的概率做离散化处理,即

由于不需要考虑出射后的动量方向,所以只需要将入射粒子和靶粒子的动量矢量夹角进行归一化,即仰角归一化.同时,当仰角确定时,整个方向角的圆周上,都应有粒子存在,所以每一仰角所占概率,应为仰角所对应的球带占整个球表面积的比例,如图4所示.

其概率为

图4 仰角对应球带(球坐标)Fig.4.Elevation corresponding to the ball belt (ball co−ordinate).

其中,ϕ为仰角,∆ϕ为离散化的角度间隔.则归一化的仰角为

对于D−D,D−T,T−D反应截面,首先计算出某一温度下粒子能够具有的速率区间[vh1,vh2],计算入射粒子动量P1和靶粒子动量P2,计算速率所对应的麦克斯韦速率分布概率f.以仰角为循环变量,从0到π,计算碰撞后的动量P,以之计算质心系能量,并计算仰角对应的概率 ΔS/S.得到质心系能量后,从数据库中读取对应截面σ,利用公式

计算出修正后的截面值,获得对应不同入射粒子能量的矩阵.

2.4 路径函数与随机函数设计

本文所模拟的等离子体只包含氘离子和氚粒子,且认为这些粒子均匀分布.因为氘氚聚变反应率需要通过路径长度计算,因此计算入射粒子云的运动路径(路径函数)较为关键.这里对路径函数做一简要说明： 首先输入入射粒子云的初始位置(x,y,z)和入射粒子云速率的三个坐标分量(vx,vy,vz).然后建立一个n× 3的矩阵,其中n代表着入射粒子云走过的网格数,3即为三个坐标值.为避免漏算网格导致误差,将原时间步长1 ns细分成时间网格10－12s,计算每一个时间网格内入射粒子云走过的网格数,进而计算路径长度.

本文中,除等离子体温度和氘氚粒子数密度已知外,对于氘氚粒子云的位置、速率、方向等均随机产生.本文采用了全部随机方法,即无论是入射粒子云的种类、位置,或是速率,还是方向,均按照一定概率进行随机抽样.由于入射粒子云不同,对应的麦克斯韦速率分布并不完全相同,所以要先判断是哪种粒子,分两种情况进行后续的部分.对于入射氘粒子云,先求出氘粒子云在此温度下的最概然速率vm,并求出两个半高宽的速率范围[vh1,vh2],然后求解麦克斯韦速率分布函数,离散化后利用函数依概率进行随机,得到一个速率值.对极角也做同样的离散化处理.对于入射氚粒子云,与氚粒子云入射算法相同,更改质量参数即可.

2.5 算法流程图

总的算法流程见图5.

3 氘氚聚变反应模拟分析

3.1 初始边界条件

为与解析公式(1)进行比较,本文中,ni=1014cm－3,氘氚粒子各占50%,并满足劳森判据.惯性约束等离子体时间步长通常在纳秒量级[18],因此本文设置时间步长为1 ns.根据2.2节的工作,设定网格大小为10－4m.模拟的温度取值为20,40,60,80,100 keV.为确保精确度,对于每一个所模拟的温度,粒子云速率离散为1000个网格,极角离散为180个网格,方位角离散为360个网格.采用MC方法对105个网格单元进行采样,采样网格单元作为入射粒子云网格来控制统计误差.

3.2 模拟结果与讨论

三种聚变反应截面的温度修正结果如图6所示(1 barn=10－24cm2).

从图6(b)和图6(c)可明显地看到多普勒效应.这是由于氘氚等离子体中,靶粒子服从麦克斯韦速率分布,因此在截面温度修正后,能量将有所展宽,且温度越高,麦克斯韦速率分布范围越宽,由此靶粒子展宽越大,同时峰值截面也逐渐减小.经截面修正,将所有满足同一麦克斯韦速率分布的粒子的截面归一化为一个截面值,解决了大量靶粒子热运动所带来的大数据量模拟计算时计算空间不足的问题,大大节约了计算时间.

氘氚等离子体聚变反应模拟结果如图7所示.

图5 氘氚等离子体聚变反应模拟流程图Fig.5.Simulation flow of deuterium−tritium plasma fusion reaction.

图6 氘氚等离子体聚变反应修正 (a) D−D反应截面修正;(b) D−T反应截面修正;(c) T−D反应截面修正Fig.6.Correction of fusion reaction of deuterium−tritium plasma： (a) D−D fusion cross section correction;(b) D−T fusion cross sec−tion correction;(c) T−D fusion cross section correction.

从图7可以看出,根据解析公式法,随等离子体温度的上升,中子产生率首先增加,在60 keV左右达到峰值,而后开始降低;而数值模拟结果则是中子产生率随等离子体温度上升而一直增加.根据三种聚变反应截面的温度修正即图6,分析可知氘−氘反应、氘−氚反应、氚−氘反应的截面在20－100 keV的范围内呈递增趋势,在此范围内,本文采用的数值模拟法相较传统解析公式法,更符合ENDF/B−VI和ENDF/B−VII的氘氚聚变反应数据库的截面数据趋势.

图7 氘氚等离子体聚变反应中子产生率Fig.7.Neutron production rate of deuterium−tritium plasma fusion reaction.

下面分析等离子体处于低温时数值模拟结果与解析公式法结果误差过大的原因.首先,对于氘氚聚变反应率,解析公式法和本文采用了不同的计算方法,具体可参见本文第2节所述.另一方面,解析方法所用截面来自核物理分析方法,本文则采用了最新的ENDF核数据库.根据图7,解析方法所用的D−T聚变截面(聚变反应最重要的过程)在60 keV左右达到峰值,所获得的中子产生率也在60 keV左右达到峰值.而本文所用数值模拟方法所获得的中子产生率和ENDF/B−VII数据库的聚变截面则在20－100 keV的范围内呈递增趋势.两种方法相比,在约20 keV,本文方法给出的聚变反应率远低于解析方法,差异约70%.图8对比分析了两种方法所对应的聚变反应截面.根据该图,D−T聚变截面在低于60 keV的范围内出现明显差异.就20 keV而言,解析法使用截面为0.4077 barn,本文方法使用截面为0.0597 barn,差异约85%.这可能是导致这两种方法之间出现显著差异的根本原因.

图8 两种方法聚变反应截面对比[19,20]Fig.8.Comparison of fusion cross sections of the two meth−ods[19,20].

4 结 论

本文针对惯性约束的激光等离子体,采用MC方法和SN方法,开展了D−T粒子云团聚变反应率的数值模拟工作.得到以下结论： 1)与原核数据库截面相比,D−T,D−D,T−D截面经修正后多普勒温度效应显著;2)在等离子体温度20－100 keV范围内,本文数值模拟结果较解析公式法更符合ENDF/B−VI和ENDF/B−VII的氘氚聚变反应数据库的截面数据;3)本文数值模拟结果与与解析公式法结果在低能区存在较大误差,可能是计算方法不同与低温截面差异过大原因导致.