宏观-微观模型和朗之万方法在低能核裂变研究中的应用

刘丽乐，陈永静，吴锡真，李祝霞，葛智刚，沈彩万，宿 阳，黄小龙，舒能川

(1.中国原子能科学研究院 核数据重点实验室，中国核数据中心，北京 102413；2.中国原子能科学研究院 核物理研究所，北京 102413；3.湖州师范学院 理学院，浙江 湖州 313000)

核裂变研究在核天体物理及放射性核束物理等基础科学领域具有重要的学术意义，同时在核能开发及核技术应用领域具有一定的应用价值。近年来，随着计算技术和测量技术的进步，核裂变在理论及实验研究上均取得了一定的进展，有必要对核裂变机制进行深入研究，加深对裂变现象的认识。

目前，核裂变理论研究的主要途径有唯象模型[1-2]、断点模型[3-4]、基于宏观-微观模型的动力学方法[5-14]及裂变微观理论[15-18]。其中基于宏观-微观模型的裂变动力学方法(即布朗运动模型)，较有代表性的是多维朗之万方法，可在合理计算时间内对裂变核从基态形变至断点的动力学过程进行描述，更关键的是该方法具有对大范围锕系核甚至超锕系核裂变碎片分布量进行定量计算的能力。早期，朗之万方法主要被用于高激发能下核裂变和重离子融合反应的研究，这种情况下微观壳效应的影响可忽略。近年来，朗之万方法逐渐被用于低能核裂变研究中，如日本Chiba等[5-7]采用朗之万方法对锕系核及超锕系核的裂变碎片质量分布及总动能分布进行了研究；美国Sierk[8]将朗之万方法应用于五维形变空间；波兰Pomorski等[9]采用三维朗之万方法对超重核性质进行了研究；美国Randrup等[10]在Smoluchouski极限条件下，采用裂变核在五维位能曲面随机游走的方法描述低能裂变碎片分布量；北京大学樊铁栓课题组[11]基于五维位能曲面开展了锕系核裂变性质的系列研究；本课题组近年来采用朗之万模型对核裂变动力学过程进行了细致的研究，包括对朗之万路径初始位置、断点位置和能级密度参数等对裂变碎片质量分布的影响，裂变断点构型[12]，不同的朗之万方程输入量对计算结果的影响[13]，以及核黏滞性对裂变过程的影响[14]等的研究。

本工作在之前研究[12-14]的基础上，进行补充和延伸，主要研究拉长形变空间和壳衰减因子对裂变碎片质量分布及总动能分布的影响，以及断点处拉长形变与质量非对称度间关联的影响，以进一步加深对裂变过程和断点构型的认识。最后，给出14 MeV中子诱发233,235U裂变碎片发射中子前与发射中子后的质量分布计算结果，并与ENDF/B-Ⅷ.0评价数据[19]进行对比。

1 模型方法

1.1 双中心壳模型

裂变过程通常被看作裂变核形状演变的过程，核形状描述是裂变研究的第一步。本工作采用双中心壳模型来描述裂变核形状，该模型是Maruhn等[20]在Nilsson模型基础上发展起来的，较适用于核裂变和重离子融合反应等大形变核的描述。模型假设原子核在形变过程中体积保持不变，核表面为一等势面，不随形变而发生变化。因核表面等密度面正比于等势面，可用等势面描述原子核表面的形状，如图1a所示。双中心壳模型中独立粒子势表示为：

图1 基于双中心壳模型描述的裂变核形状(a)和对称轴z方向的真实势能及形变谐振子势(b)Fig.1 Nuclear shape described with two-center shell model (a) and corresponding actual potential and deformed oscillator potential along symmetry axis z (b)

(1)

其中，z′i=z-zi，i=1，2。沿对称轴z方向可将势中心分为z1和z2两部分(图1b)，通过z=0所在横截面分割开来。通过一定的约束条件，即z=0、z=z1和z=z23个位置的势能连续和势能对z的一阶导连续，以及原子核形变中体积守恒，则式(1)引入的12个参数可由5个自由的形状参数来表示：{Z0，δ1，δ2，η，ε}。其中，Z0为拉长形变参数，表示双中心之间的距离，即Z0=z2-z1，通常以无量纲的Z0/R0表示，R0为与形变核体积相同的球形核半径；δi为碎片变形参数，描述左、右两端碎片的变形程度，δi=(3βi-3)/(1+2βi)，βi=ai/bi，i=1，2；η为质量非对称度，η=(V2-V1)/(V2+V1)=(Af2-Af1)/(Af2+Af1)，V1和V2分别为左、右两部分的体积，Af1和Af2分别为相应的质量数；ε为颈部参数，表示在z=0处沿对称轴z方向的真实谐振子势与常规谐振子势比值，即ε=E/E′=1+cz′+dz′2，如图1b所示。

1.2 宏观-微观模型

裂变位能曲面是朗之万方程的重要输入量之一，其对于理解裂变机制及裂变产物的各种性质具有决定性的意义。本工作采用宏观-微观模型计算位能曲面，该模型中裂变核的位能由描述平滑变化趋势的宏观液滴能和局部涨落效应的微观修正能两部分构成。由于裂变过程中最重要的是形变位能的影响，宏观液滴能主要包含与形变相关的表面能和库仑能，采用有限程液滴模型[21]计算。微观修正能包含壳修正能和对修正能，分别采用Strutinsky方法[22-23]和BCS方法[24]计算，其中的单粒子能级基于双中心壳模型得到。以236U为例，宏观液滴能和总形变能在拉长形变-质量非对称度二维空间的投影如图2所示。可见，宏观液滴能随着拉长形变的增加，呈现先上升后下降的趋势，这是表面能和库仑能相互竞争的结果，同时可看到，当仅考虑宏观液滴能时，最优裂变路径为对称裂变道(图2a中虚线)。当考虑微观修正能时，第1个鞍点位于{Z0/R0=0.5，η=0.0}附近，随着拉长形变的增加，对称裂变道第2个鞍点的位能较非对称裂变道鞍点位能高3～5 MeV，此时最优裂变路径为非对称裂变道(图2b中虚线)。

位能曲面不仅与原子核的形状有关，还与核温度有关。一般随着核温度的升高，微观壳效应和对效应会明显减弱，相应地微观修正能会减小。而核温度不够高时，宏观液滴能随温度变化较小，可忽略其温度相关性。因此，含温度的位能曲面可表示为：

ε=0.35；δ=0.18图2 236U的宏观液滴能(a)及总形变能(b)在Z0/R0-η空间的投影Fig.2 Macroscopic liquid drop energy (a) and total potential energy (b) of 236U projected on Z0/R0-η plane

Etot(q,T)=Emac(q)+δEmic(q,T=0)Φ(T)

Φ(T)=exp(-aT2/Ed)

(2)

式中：Φ(T)为温度相关性因子，采用Ignatyuk提出的形式；a为能级密度参数，由于裂变碎片分布量计算对a不敏感，本工作将其设为常数，取a=Acn/10 MeV-1，Acn为裂变复合核质量数；Ed为壳衰减因子，本工作Ed取值60 MeV，与Randrup等[25]的取值一致。

1.3 朗之万方程

裂变动力学过程采用朗之万方程描述，如式(3)所示，其描述了裂变过程中广义坐标qi和广义动量pi随时间变化的规律。

(3)

式中：n为动力学计算中所考虑的形变空间的维数，即可自由变化的形状参数的数量；V为形变位能；mij和γij为质量张量和粘滞张量，分别采用Werner-Wheeler近似方法[26]和墙加窗一体模型[27]计算；gij和Γj分别为随机力强度和归一化的随机力，将归一化的随机力假设为白噪声，则满足如下关系：

〈Γi(t)〉=0

〈Γi(t1)Γj(t2)〉=2δijδ(t1-t2)

(4)

可通过高斯随机数来模拟归一化的随机力Γj。随机力强度gij根据涨落-耗散原理计算：

gikgjk=γijT*

(5)

其中，T*为有效核温度，考虑了较低激发能时量子效应的影响，其与常规核温度T的关系[28]为：

(6)

由式(6)可看到，当核温度T=0 MeV时，T*=1.0 MeV；当T较高时，T*逐渐与T相等。有效核温度的引入主要是为了合理描述低激发能时由零点能带来的涨落。根据费米气体模型，核温度T与核体系内能Eint的关系为Eint=aT2。Eint由式(7)计算：

(7)

式中，E*为裂变核的总激发能，约为入射粒子动能与相应结合能之和。在朗之万方程数值计算中，每一步相应的核形状和集体动能均不同，相应的核体系内能不同，于是核温度不同，因此每一步均需计算相应的有效核温度。

本工作应用三维朗之万方法，以14 MeV中子诱发235U(14 MeV n+235U)为例描述裂变过程。在双中心壳模型形状参数描述中，固定颈部参数ε=0.35，并令左、右两碎片变形程度相同，即δ1=δ2=δ，最终三维形变空间的广义坐标为{Z0/R0，δ，η}。在朗之万动力学计算中，每条到达断点的朗之万路径即为1个裂变事件，为节省计算时间，设定朗之万计算从第1鞍点态{Z0/R0=0.5，δ=0.2，η=0.0}出发，可得到与基态出发相同的结果。本工作中断点根据裂变核的颈部半径来定义，通常断点处的颈部半径近似于核子尺度，假设当裂变核颈部半径小于0.5 fm即到达断点。

2 计算结果

为节省计算时间，在动力学计算前需建立位能曲面、质量张量和粘滞张量在三维形变空间的网格点，朗之万计算中采用三维抛物线插值的方法得到相应取值。在三维形变空间中，双中心距离Z0/R0、变形参数δ和质量非对称度η的变化范围及相应步长如下：

Z0/R0=-0.3 (0.1) 3.5

δ=-0.45 (0.03) 0.81

η=-0.62 (0.04) 0.62

其中，括号中的数值为各形状参数相应的步长。对于每个裂变系统，朗之万计算从第1个鞍点态出发，共模拟25万条朗之万路径可得到稳定的碎片分布量。

2.1 拉长形变空间对裂变碎片质量分布、总动能分布及断点构型的影响

图3 拉长形变空间对14 MeV n+235U裂变碎片质量分布和总动能分布的影响Fig.3 Effect of elongated deformation space on fragment mass distribution and total kinetic energy distribution of 14 MeV n+235U fission

Z0max分别取3.0R0、3.5R0和4.0R0时断点处Z0/R0与η的关联情况如图4所示。从图4可看到，当Z0max分别取3.5R0和4.0R0时，Z0/R0与η关联曲线基本一致，而当Z0max取3.0R0时，对称裂变区域核拉长明显偏低，这与图3b相应的裂变碎片总动能在对称裂变区域偏低的情况相符。同时可看到，在对称裂变区域核拉长最大，表明对称裂变道对应于超长形变。随着质量非对称度η的增大，核拉长逐渐减小，当η≈0.175(重碎片质量数AH=138)时，核拉长最小，这是源自Z=50、N=82壳效应的影响，导致满壳附近的碎片核形状接近球形。当η>0.175时，随着质量非对称度的增大，核拉长逐渐变大，这是由于碎片核逐渐远离满壳，相应的碎片核形变增大所致。可见，断点处核拉长与质量非对称度的关联也说明图3b中碎片总动能随质量非对称度的变化规律主要源自相应核拉长的贡献。通过研究核拉长形变空间对裂变碎片分布量计算的影响，确定了在朗之万动力学计算中，核拉长边界值Z0max至少应取3.5R0，为节省计算时间，本工作取Z0max=3.5R0。

图4 核拉长形变空间对断点处核拉长与质量非对称度关联的影响 Fig.4 Effect of elongated deformation space on correlation between elongation and mass asymmetry at scission point

2.2 壳衰减因子对裂变碎片质量分布、总动能分布及断点构型的影响

为描述一定激发能下的裂变碎片分布量，位能曲面的温度相关性通过式(2)引入，式中壳衰减因子Ed影响壳效应随激发能升高而减弱的快慢程度，本工作研究了壳衰减因子Ed对裂变碎片质量分布、总动能分布及断点构型的影响，结果如图5、6所示。

图5 壳衰减因子对14 MeV n+235U裂变碎片质量分布和碎片总动能分布的影响Fig.5 Effect of shell damping parameter on fragment mass distribution and total kinetic energy distribution of 14 MeV n+235U fission

图5a为Ed分别取0、30、40、50、60 MeV时相应质量分布计算结果的比较，其中Ed=0对应于位能曲面上微观修正能消失的情况。从图5a可看到，当仅考虑宏观液滴能(Ed=0)时，质量分布为高斯形状的单峰分布，表明仅有对称裂变道的贡献。随着Ed的增大，壳效应强度逐渐增加，相应的质量分布峰区产额升高，谷区产额下降，说明对称裂变道的贡献降低而非对称裂变道的贡献增大，可见，Ed对质量分布峰谷比有较大的影响。在本工作中，通过对14 MeV中子诱发235U裂变碎片质量分布的计算，确定Ed的取值为60 MeV，在其他裂变系统计算中，不再调节Ed值。应说明的是，基于微观理论可实现能量相关的裂变位垒的较好描述[30]，原则上，Ed可通过微观计算获得。图5b为不同Ed下相应裂变碎片总动能分布计算结果的比较，可看到，当Ed=30、40、50、60 MeV时，计算结果相差不大，虽然在碎片质量数138附近，随着Ed的增大碎片总动能出现轻微升高，是由于Ed增大相应的壳效应增强引起的，但该影响极小。同时可看到，当仅考虑宏观能(Ed=0)时，碎片总动能在对称裂变区域最大，随着质量非对称度的增大而逐渐减小。

为进一步分析Ed对碎片总动能分布影响的原因，本工作研究了不同Ed下断点处核拉长Z0/R0随质量非对称度η的变化规律，如图6所示。可看到，当仅考虑宏观能时，核拉长随质量非对称度的变化较小，因此图5b显示的Ed=0时碎片总动能随着质量非对称度增大而减小的变化规律主要源自两碎片相应电荷项Z1Z2=(ZCN/ACN)2A1A2的贡献，当两碎片质量数相等时Z1Z2最大，相应碎片总动能也最大，随着质量非对称度的增加，两碎片相应的Z1Z2减小，于是总动能逐渐降低。当Ed=30、40、50、60 MeV时，对称裂变区域核拉长随Ed的增大略有增大，而在η=0.17附近核拉长随Ed的增大而略有降低，这与图5b结果基本相符，这是由于Ed增大使得壳效应增强所致。

图6 壳衰减因子对断点处核拉长与质量非对称度关联的影响Fig.6 Effect of shell damping parameter on correlation between elongation and mass asymmetry at scission point

2.3 14 MeV n+233,235U裂变碎片质量分布计算结果

基于三维朗之万模型，计算了14 MeV中子诱发233,235U裂变碎片发射中子前(pre-n)的质量分布，并与ENDF/B-Ⅷ.0评价数据进行比较，结果如图7中蓝线所示。计算的质量分布中，峰高、峰谷及峰宽均与评价数据符合，而轻、重峰位较评价数据向右偏离2～3个核子数，这是由于当前朗之万计算结果为初始碎片质量分布，而评价数据为碎片发射中子后的质量分布，因此峰位的轻微偏离是合理的。

图7 14 MeV n+233,235U裂变碎片发射中子前后质量分布计算结果与ENDF/B-Ⅷ.0评价数据的比较Fig.7 Calculated fragment mass distribution (pre- and post-neutron emission) in 14 MeV n+233,235U fission compared with evaluated data from ENDF/B-Ⅷ.0

3 结论

本工作基于三维朗之万模型对低能中子诱发裂变碎片质量分布、总动能分布和断点构型进行了研究，其中位能曲面采用基于双中心壳模型的宏观-微观模型计算得到。首先以14 MeV中子诱发235U裂变为例，研究了拉长形变空间对裂变碎片分布量和断点构型的影响，研究发现拉长形变空间对对称裂变道的影响较大，这是因为在轻锕系核区，对称裂变道对应于超长形变道，而较小的拉长形变空间会限制朗之万路径中拉长较大的核形状的出现。通过本研究，确定了模型计算中拉长形变空间边界至少应为3.5R0。其次，研究了温度相关的微观修正能公式中壳衰减因子对裂变碎片分布量和断点构型的影响，发现壳衰减因子对质量分布的影响较大，随着其数值的增大，质量分布峰谷比逐渐变大，这是由壳效应增强所引起，而其对碎片总动能分布及断点处核拉长影响较小，说明一定程度上壳效应的强弱对各质量非对称度相应的平均核拉长影响较小。通过计算，确定了本工作中壳衰减因子的合理取值为60 MeV。最后，基于该模型计算得到了14 MeV中子诱发233,235U裂变碎片质量分布，与ENDF/B-Ⅷ.0评价数据符合较好，说明该模型具有定量计算裂变碎片质量分布的能力。