宏觀-微觀模型和朗之萬方法在低能核裂變研究中的應用

劉麗樂，陳永靜，吳錫真，李祝霞，葛智剛，沈彩萬，宿 陽，黃小龍，舒能川

(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時，隨著質量非對稱度的增大，核拉長逐漸變大，這是由于碎片核逐漸遠離滿殼，相應的碎片核形變增大所致?？梢姡瑪帱c處核拉長與質量非對稱度的關聯也說明圖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評價數據符合較好，說明該模型具有定量計算裂變碎片質量分布的能力。