基于傅里葉級數(shù)展開的裂變位能曲面與碎片質(zhì)量分布計算

宿 陽，劉麗樂，陳永靜，葛智剛

(中國原子能科學(xué)研究院 核數(shù)據(jù)重點(diǎn)實(shí)驗(yàn)室，北京 102413)

自1938年原子核裂變現(xiàn)象被發(fā)現(xiàn)以來，原子核裂變理論研究一直是核物理研究的重要領(lǐng)域。正如1939年Meitner、Frisch、Bohr和Wheeler[1-2]的開創(chuàng)性文章中所討論的，原子核裂變可定性看作是原子核形狀從單一復(fù)合核演化到斷裂成兩個碎片的過程。目前原子核裂變研究在實(shí)驗(yàn)和理論兩個方面有了長足的進(jìn)步，但仍缺乏能對裂變過程嚴(yán)格和精確描寫的理論。

裂變碎片質(zhì)量分布是原子核裂變的重要觀測量，人們發(fā)展了各種方法用于計算裂變碎片質(zhì)量分布，如唯像方法[3-6]，這種方法通常是基于已有的物理圖像和基本概念，采用簡化的理論公式并引入一些模型參數(shù)，根據(jù)實(shí)驗(yàn)數(shù)據(jù)擬合確定參數(shù)再現(xiàn)裂變產(chǎn)額實(shí)驗(yàn)數(shù)據(jù)或其他可觀測量。唯像方法通常在確定其參數(shù)的核區(qū)和能區(qū)可再現(xiàn)實(shí)驗(yàn)數(shù)據(jù)，但外推到未知區(qū)域時有很大的不確定性。此外，唯像模型未考慮動力學(xué)效應(yīng)，無法描述裂變動力學(xué)過程，不能加深人們對裂變過程的理解。

另一種是基于能量密度泛函理論的原子核裂變微觀理論方法[7-8]。裂變微觀理論從核子-核子有效相互作用出發(fā)，原則上無需自由參數(shù)，具有較好的自洽性和預(yù)言能力。裂變微觀理論分絕熱和非絕熱兩類，絕熱的微觀動力學(xué)理論研究首先需計算原子核位能曲面，在位能曲面基礎(chǔ)上進(jìn)行動力學(xué)計算給出裂變碎片分布；非絕熱的裂變動力學(xué)研究不需計算位能曲面，但只適合研究裂變過程和單個裂變事件，只能給出碎片最可幾分布。需要強(qiáng)調(diào)的是，微觀理論方法需較大的計算機(jī)資源，非常耗時，目前也只對有限的偶偶核裂變系統(tǒng)進(jìn)行研究。

基于宏觀-微觀模型的裂變理論是目前發(fā)展較成熟的一種計算裂變碎片分布量的方法，這種方法首先要給定原子核的形狀，然后基于宏觀-微觀方法計算裂變位能曲面，在位能曲面基礎(chǔ)上進(jìn)行含時或不含時的動力學(xué)計算給出碎片質(zhì)量分布。這種方法既考慮了裂變核的靜態(tài)性質(zhì)，如位能曲面，又考慮了裂變動力學(xué)效應(yīng)，是理論基礎(chǔ)較好又具有一定預(yù)言能力的方法。相對裂變唯像模型，宏觀-微觀模型具有一定的物理基礎(chǔ)，相對微觀理論模型，宏觀-微觀模型無需大的計算資源。因此，宏觀-微觀模型是近年來有望成為替代唯像模型進(jìn)行裂變后數(shù)據(jù)計算和評價的工具，成為國際上裂變研究的熱點(diǎn)。

宏觀-微觀模型方法的基礎(chǔ)是裂變過程中原子核形狀描述，給定原子核一個形狀，求出與這個形狀相關(guān)的勢能，因此可把原子核裂變過程描述為原子核形狀在勢能面上的演化。各種宏觀-微觀模型方法的不同主要在于采用的核形狀描述方法不同。如美國洛斯阿拉莫斯國家實(shí)驗(yàn)室的Nix、M?ller、Randrup等采用連接二次曲面描述核形狀[9-10]，德國法蘭克福大學(xué)Maruhn和Greiner提出雙中心殼模型[11]，中國原子能科學(xué)研究院Liu等、日本東京工業(yè)大學(xué)Chiba等采用雙中心殼模型描述核形狀[12-15]來研究裂變產(chǎn)物，北京大學(xué)樊鐵栓課題組則采用廣義Lawrence形狀描述法[16-17]。近年來，波蘭Pomorski等[18]發(fā)展了基于傅里葉級數(shù)展開的核形狀描述方法，這種核形狀描述方法優(yōu)勢在于用少量的形變參數(shù)就可描述原子核裂變過程中的各種形狀，且計算速度快，無需迭代計算。

235U是最重要的燃料核，其中子誘發(fā)的裂變產(chǎn)額數(shù)據(jù)是核反應(yīng)堆設(shè)計與運(yùn)行維護(hù)、核裝置設(shè)計與測試、軍控核查等的重要數(shù)據(jù)，同時也是裂變物理研究的基礎(chǔ)數(shù)據(jù)。本工作以中子誘發(fā)235U形成的復(fù)合核236U為例，對236U裂變過程進(jìn)行研究，分析236U位能曲面隨不同形變參數(shù)的變化，計算236U裂變碎片質(zhì)量分布，分析模型參數(shù)以及核溫度對碎片質(zhì)量分布的影響。

1 傅里葉級數(shù)展開的核形狀描述方法

用盡量少的形變參量精確描述原子核形狀是一項(xiàng)艱巨的任務(wù)，特別是與裂變有關(guān)的核形狀，人們發(fā)展了一些有效的核形狀描述方法[9-11，16，18]，本文采用文獻(xiàn)[18]提出的傅里葉級數(shù)展開的核形狀描述方法描述裂變過程中的核形狀。圖1示出了在拉長形變、質(zhì)量不對稱度以及頸部寬度3個自由度下從鞍點(diǎn)到斷點(diǎn)過程中典型的原子核形狀示意圖。

圖1 柱坐標(biāo)下原子核形狀示意圖Fig.1 Schematic visualization in cylindrical coordinates of nuclear shape

在柱坐標(biāo)(ρ,z,φ)下，軸對稱核形狀方程可用傅里葉級數(shù)展開如下：

(1)

(2)

假定在給定的z處，原子核橫向截面的形狀(圖1中紅色虛線)用橢圓的形式來處理，a(z)和b(z)分別為橢圓的半短軸和半長軸，也可得到裂變過程中非軸對稱的原子核形狀，具體形式如下：

(3)

這種傅里葉級數(shù)表示的形狀參量可將核形狀展開到任意階數(shù)，從而達(dá)到想要的精度。與常見的使用球諧函數(shù)描寫的核形狀相比，傅里葉形狀參量能快速擬合核形狀，在6階內(nèi)即可得到很好擬合結(jié)果，而球諧函數(shù)通常需展開到12階以上。圖2所示為使用傅里葉級數(shù)描述的核形狀，傅里葉級數(shù)截斷到6階。

圖2 傅里葉級數(shù)展開各階系數(shù)對核形狀的貢獻(xiàn)Fig.2 Contribution of Fourier series coefficients to shape

雖然傅里葉級數(shù)展開系數(shù)an與核形狀有直接聯(lián)系，前3階a2、a3、a4分別對應(yīng)四極、八極和十六極形變，對裂變而言，這3種形變分別對應(yīng)裂變核拉長量、質(zhì)量不對稱度和頸部寬度，但它們之間的對應(yīng)關(guān)系較復(fù)雜。如式(2)中原子核拉長量c與展開系數(shù)中偶數(shù)項(xiàng)a2n的關(guān)系，首先，a2對c占主要影響，但c與所有的偶數(shù)項(xiàng)系數(shù)有關(guān)。其次，c與a2n呈負(fù)相關(guān)的關(guān)系。這導(dǎo)致通過系數(shù)a2直接反映拉長量時，對一些關(guān)鍵的裂變過程進(jìn)行分析會變得復(fù)雜。因此，可用展開系數(shù)an組合1組新的集體坐標(biāo){qn}來表示形變參量，這樣得到的集體形變參量可體現(xiàn)裂變過程的特征。新定義集體坐標(biāo){qn}如下：

(4)

表1 球形核對應(yīng)的傅里葉級數(shù)展開系數(shù)(奇數(shù)項(xiàng)為零)Table 1 Values of Fourier expansion coefficients for spherical shape (Odd coefficients are zero)

集體坐標(biāo)中，q2、q3和q4分別對應(yīng)原子核拉長量、質(zhì)量不對稱度和頸部寬度。根據(jù)文獻(xiàn)[19]，q5和q6對小形變區(qū)液滴能的影響可忽略，在大形變區(qū)(剪裂線附近)也僅不超過0.5 MeV的影響。所以在之后的計算中不考慮q5和q6。上述集體坐標(biāo)與原子核不對稱形變自由度η一起構(gòu)成了四維形變空間。

基于上述的核形狀描述，采用宏觀-微觀方法計算裂變核的裂變位能曲面。對于一給定的形變的原子核，總能量Etot(η,q2,q3,q4)由式(5)計算。

Etot=Emac+Emic

(5)

其中，宏觀能Emac通常在液滴模型框架內(nèi)計算。文獻(xiàn)[20]的Lublin-Strasbourg Drop (LSD)模型很好再現(xiàn)了原子核質(zhì)量和裂變位壘高度，本文用LSD模型計算給出宏觀能。微觀修正能包括殼修正能和對修正能，分別采用標(biāo)準(zhǔn)的Strutinsky方法[21]和BCS方法[22]計算，其中用到的單粒子能級在給定形變參數(shù)(η,q2,q3,q4)下由Yukawa-Folded勢進(jìn)行計算。

Etot=Emac+Emicφ(T)

(6)

溫度修正因子隨溫度的變化趨勢如圖3所示，當(dāng)溫度為0時，φ(T)接近1，當(dāng)溫度增大到3 MeV時，φ(T)接近0，說明此時微觀修正能對位能曲面已無貢獻(xiàn)，只剩下宏觀能的貢獻(xiàn)。

圖3 溫度函數(shù)φ(T)形式Fig.3 Form of temperature function φ(T)

2 集體模型描述裂變過程

裂變動力學(xué)過程的基本思想是原子核在裂變過程拉長方向(q2方向)的運(yùn)動是相對緩慢的，同時伴隨著垂直于裂變方向(q3，q4方向)上的快速振動，在Born-Oppenheimer近似(BOA)下(BOA近似認(rèn)為當(dāng)1個多自由度體系不同自由度的運(yùn)動快慢相差懸殊時，可近似把這些自由度分開處理)，總能量為E的裂變核體系波函數(shù)形式如下：

ΨnE(q2,q3,q4)=unE(q2)φn(q3,q4;q2)

(7)

其中：unE(q2)為裂變方向運(yùn)動的本征函數(shù)，主要與原子核拉長形變q2有關(guān)，可通過WKB近似得到[24]；φn(q3,q4;q2)模擬給定拉長形變q2后q3、q4方向上的n聲子的快速集體振動。對于低能裂變的情況，在給定拉長形變q2時，(q3±dq3,q4±dq4)區(qū)域內(nèi)的概率密度W(q3,q4;q2)可通過如下方式得到：

(8)

模型的進(jìn)一步簡化是使用下述形式的Wigner函數(shù)代替式(8)中概率密度。

(9)

其中：Vmin(q2)為給定拉長形變q2下最小的位能；T*為廣義溫度[25]，考慮了裂變核的熱激發(fā)和集體零點(diǎn)能E0。

T*=E0/tanh(E0/T)

(10)

質(zhì)量數(shù)為A的原子核的溫度T是通過唯像的關(guān)系(E*=aT2)從核的熱激發(fā)能(E*)計算得到的，其中a=A/(10 MeV)。當(dāng)T很小時，廣義溫度T*近似等于零點(diǎn)能；當(dāng)T遠(yuǎn)大于E0時，T*近似等于T。零點(diǎn)能E0是本文動力學(xué)計算過程中的一個可調(diào)參量。

為得到給定拉長形變q2處的裂變碎片質(zhì)量分布，需對概率密度W(q3,q4;q2)中的q4參量進(jìn)行積分：

(11)

裂變概率很大程度上依賴于頸部的粗細(xì)，即頸部半徑Rneck。假設(shè)頸部斷裂概率P等于：

(12)

其中：Pneck為與頸部粗細(xì)呈正比的頸部斷裂概率；k0/k為裂變方向較大的集體速度,在最終計算碎片質(zhì)量分布時會直接消除。本文頸部斷裂概率Pneck使用Gauss函數(shù)的形式來描述：

Pneck(Rneck)=exp[-lg 2(Rneck/d)2]

(13)

其中，d描述的是頸部斷裂概率的半寬度，是本文動力學(xué)計算過程中的第2個可調(diào)參量。

引入頸部斷裂概率Pneck(q3,q4;q2)后，概率密度分布(式(11))中對q4的積分可寫為：

(14)

在上述近似中，對于一個給定的q3，裂變可能會以不同的概率發(fā)生在q2的一段區(qū)域內(nèi)。若認(rèn)為一部分原子核在q′2處發(fā)生裂變，則在下一步q2處就不能再考慮。因此，為得到在給定q2處真實(shí)的裂變概率密度w′(q3;q2)，需排除q′2

(15)

對q2進(jìn)行積分后再進(jìn)行歸一處理就得到裂變碎片質(zhì)量分布：

(16)

由于q3和質(zhì)量不對稱度存在一一對應(yīng)的關(guān)系，由式(16)即可得到與實(shí)驗(yàn)數(shù)據(jù)可直接比較的裂變碎片質(zhì)量分布。

3 計算結(jié)果及其分析

3.1 各集體坐標(biāo)下的位能曲面分析

基于LSD模型和Yukawa-Folded勢的宏觀-微觀模型計算了n+235U形成的復(fù)合核236U在非軸對稱形變、拉長、左右碎片質(zhì)量不對稱以及頸部寬度4個集體自由度下的裂變位能曲面。計算時形變參數(shù)空間(η,q2,q3,q4)網(wǎng)格點(diǎn)的取值如下：

η=0.00(0.03)0.21

q2=-0.60(0.05)2.35

q3=0.00(0.03)0.21

q4=-0.21(0.03)0.21

(17)

基于上述形變參數(shù)網(wǎng)格點(diǎn)，由式(6)可計算不同溫度下每個網(wǎng)格點(diǎn)上的原子核的勢能。

圖4 (q2,η)平面上236U位能曲面 (在q3、q4方向上能量取極小)Fig.4 Potential energy surface in (q2,η) plane for 236U (as obtain after minization with respect to q3 and q4)

圖4為(q2,η)平面上236U位能曲面(零溫)，對應(yīng)非對稱裂變。可看出，從q2>0.1開始，η=0的方向，能量均是最小的，即不存在使能量減小的其他η。其次，η=0方向存在明顯的內(nèi)外壘結(jié)構(gòu)，內(nèi)壘位于0.5

圖5示出了236U在頸部寬度q4取最小時位能曲面在(q2,q3)平面的投影。從圖5可看出，236U具有雙峰裂變位壘結(jié)構(gòu)，即對稱裂變內(nèi)壘和不對稱裂變外壘。基態(tài)位置處于q2=0.35、q3=0附近，隨著原子核繼續(xù)拉伸，在q2=0.6、q3=0附近出現(xiàn)第1個鞍點(diǎn)，為對稱裂變內(nèi)壘；越過內(nèi)壘原子核繼續(xù)拉伸在對稱裂變方向出現(xiàn)第2個極小(圖中b點(diǎn))，也即同質(zhì)異能態(tài)；從第2極小出發(fā)，在非對稱裂變方向q2=1.0、q3=0.09處出現(xiàn)第2個鞍點(diǎn)(圖中c點(diǎn))，為不對稱裂變外壘?？汕宄吹?，q3>0.06以后，在1.2

圖5 236U在(q2,q3)平面的位能曲面Fig.5 Potential energy surface in (q2,q3) plane for 236U

更詳細(xì)的位能曲面結(jié)果如圖6所示，分別給出了不同拉長形變下的位能曲面分布圖，從圖可很清楚看出，在q2=1.05之前(圖6a、b)，(q3,q4)平面僅存在對稱裂變方向上(q3=0)的1個能量最小點(diǎn)，對應(yīng)圖5中a、b兩點(diǎn)的形狀。當(dāng)原子核繼續(xù)拉伸到q2=1.05，即第2鞍點(diǎn)處，即圖6c，原子核勢能處于一較高水平，此時各種形狀的核能量均較高，能量最低點(diǎn)出現(xiàn)在非對稱裂變方向(圖6c)，說明了不對稱裂變外壘存在的真實(shí)性(圖5中c點(diǎn))。越過外壘后原子核到達(dá)裂變谷所在的區(qū)域，從q2=1.2的位能曲面(圖6d)可看出，由于對稱裂變路徑上較高外壘的存在，非對稱裂變谷先于對稱裂變谷出現(xiàn)(圖5中d點(diǎn))，圖6d中僅1個能量最小點(diǎn)，對應(yīng)的左右不對稱自由度q3=0.1，說明非對稱裂變在兩種裂變模式的競爭中處于領(lǐng)先位置。當(dāng)q2=1.6時(圖6e)，從圖6e可看出，對稱裂變谷已出現(xiàn)，即圖6e中的g點(diǎn)，對應(yīng)的q3=0，而非對稱裂變谷(q3=0.1)的深度更深。隨著原子核繼續(xù)拉長，在q2=2.0后，非對稱裂變谷較對稱裂變谷深得多，說明在236U低能裂變過程中，非對稱裂變在裂變過程中擁有優(yōu)勢，這與前面的位能曲面分析以及碎片質(zhì)量分布的實(shí)驗(yàn)結(jié)果一致。從圖6a到圖6f位能曲面的變化可看出，原子核從最初的1條裂變路徑演化為2條裂變路徑的過程中頸部自由度q4的演化起重要作用。

圖6 給定q2值的(q3,q4)平面的位能曲面Fig.6 (q3,q4) potential energy surfaces for selected values of q2

圖7為考慮了溫度對微觀修正能的影響后，計算得到的不同溫度236U位能曲面在(q2,q3)平面上的投影。從圖7可看出，隨著溫度升高，位能曲面對稱裂變路徑上的內(nèi)壘(q2=0.6，q3=0)以及外壘(q2=1.2，q3=0)的高度均降低。其次，可看出，隨著溫度的升高，位能曲面的結(jié)構(gòu)發(fā)生明顯的變化，這是因?yàn)槲⒂^修正能的影響隨溫度增加而減弱。當(dāng)溫度T為2.0 MeV時，可清楚看出位能曲面趨近于僅宏觀能的貢獻(xiàn)(等勢面比較規(guī)則)。從圖中還可看出，隨著溫度的升高，q3約為0.10對應(yīng)的非對稱裂變道逐漸消失，到T=2.0 MeV只存在對稱裂變道(q3=0)。

3.2 裂變產(chǎn)額計算結(jié)果分析

基于前面計算的位能曲面和第2章介紹的計算裂變碎片質(zhì)量分布的簡單集體模型，本文計算了零溫236U裂變碎片質(zhì)量分布(放中子前)，并與熱中子誘發(fā)235U裂變碎片質(zhì)量分布進(jìn)行了比較，如圖8所示，計算中用到的兩個參數(shù)E0和d分別取3.0 MeV和1.0 fm。可看出，理論計算結(jié)果與實(shí)驗(yàn)[27]符合較好，特別是輕重碎片質(zhì)量分布的峰位。

圖7 不同溫度下(q2,q3)平面上236U位能曲面Fig.7 Potential energy surface in (q2,q3) plane for 236U with different temperatures

圖8 236U裂變碎片質(zhì)量分布計算結(jié)果 以及與實(shí)驗(yàn)值比較Fig.8 Fission fragment mass distribution of 236U and compared with experiment data

計算裂變碎片分布量的模型中有兩個可調(diào)參數(shù)E0(零點(diǎn)能修正)和頸部斷裂幾率中出現(xiàn)的半寬度參數(shù)d。圖9給出了不同E0(圖9a)和d(圖9b)的碎片質(zhì)量分布計算結(jié)果。可看出，E0主要影響裂變碎片分布的寬度，隨著E0的增大，碎片分布寬度也在增大，由于歸一的影響，分布寬度的變化影響了峰的高度。而與頸部斷裂相關(guān)的參數(shù)d對碎片質(zhì)量分布峰位以及對稱裂變和不對稱裂變產(chǎn)額的貢獻(xiàn)均有輕微影響。隨著d的增加，重峰的峰位向碎片質(zhì)量大的方向移動。需強(qiáng)調(diào)，低能裂變碎片質(zhì)量分布的主要特征是由裂變核的位能曲面特征決定的。

圖9 236U裂變碎片質(zhì)量分布計算結(jié)果 隨E0和d的變化Fig.9 Fragment mass distribution of 236U for different values of E0 and d

圖10示出了不同溫度236U裂變碎片質(zhì)量分布計算結(jié)果。從圖10可看出，隨著溫度的增加，236U裂變碎片質(zhì)量分布的不對稱裂變模式貢獻(xiàn)逐漸降低，對稱裂變模式貢獻(xiàn)逐漸增加。這一結(jié)果與圖7中位能曲面隨溫度變化趨勢一致。圖7中溫度為2.0 MeV時，計算的位能曲面中非對稱裂變道已消失，與此處溫度為2.0 MeV時裂變碎片質(zhì)量分布主要呈對稱裂變模式一致。

圖10 不同溫度情況下236U裂變碎片質(zhì)量分布Fig.10 Fragment mass distribution of 236U with different temperatures

4 結(jié)論

基于傅里葉展開級數(shù)的核形狀描述方法結(jié)合宏觀(Lublin-Strasbourg-Drop)-微觀(Yukawa-Folded)模型計算了中子誘發(fā)235U的復(fù)合核236U在拉長形變、質(zhì)量不對稱度以及頸部寬度自由度下的裂變位能曲面，分析了裂變過程中236U位能曲面隨不同集體自由度的變化特征。用基于Born-Oppenheimer近似的三維集體模型計算了236U裂變碎片質(zhì)量分布，分析了兩個主要參數(shù)零點(diǎn)能、頸部斷裂概率半寬度對裂變碎片質(zhì)量分布的影響，研究了溫度效應(yīng)對裂變位能曲面以及裂變碎片質(zhì)量分布計算結(jié)果的影響。研究結(jié)果表明傅里葉展開級數(shù)可合理描述原子核從基態(tài)演化到斷點(diǎn)的裂變過程中各種核形狀，LSD模型計算的宏觀能和基于Yukawa-Folded勢給出的微觀修正可清楚給出236U裂變位能曲面在裂變核演化過程中隨不同集體自由度的變化特征?；谏鲜鑫荒芮娴娜S集體模型再現(xiàn)了236U裂變碎片質(zhì)量分布的主要特征，不同溫度裂變碎片質(zhì)量分布隨溫度變化趨勢的計算結(jié)果與位能曲面分析結(jié)果一致，并且與現(xiàn)有碎片質(zhì)量分布實(shí)驗(yàn)測量結(jié)果的變化趨勢相符。