快譜小堆堆芯三維瞬態(tài)核-熱-力耦合計(jì)算方法研究

胡田亮,姜奪玉,李 達(dá),王立鵬,曹 璐,張信一,蘇春磊

(西北核技術(shù)研究所 強(qiáng)脈沖輻射環(huán)境模擬與效應(yīng)國(guó)家重點(diǎn)實(shí)驗(yàn)室,陜西 西安 710024)

快譜小堆堆芯小、能譜硬、溫度高、中子泄漏強(qiáng),堆芯核-熱-力之間存在著復(fù)雜的耦合關(guān)系。不同于熱譜反應(yīng)堆,溫度載荷作用下的力學(xué)效應(yīng)導(dǎo)致的幾何變形是快譜小堆堆芯反應(yīng)性反饋的重要來(lái)源[1],與傳統(tǒng)堆芯中圍繞工質(zhì)(如水)的臨界沸騰而延伸出來(lái)的傳統(tǒng)堆芯安全分析方法不同,快譜小堆中的高力學(xué)行為與耦合效應(yīng)也是快譜小堆安全分析中最為關(guān)鍵的內(nèi)容之一[2]。因此必須研究基于核-熱-力耦合的計(jì)算方法,為開(kāi)展快譜小堆的瞬態(tài)安全分析奠定基礎(chǔ)。

快中子脈沖反應(yīng)堆在爆發(fā)脈沖時(shí),可在瞬發(fā)超臨界狀態(tài)運(yùn)行,裂變速率會(huì)在幾毫秒內(nèi)增加數(shù)個(gè)數(shù)量級(jí),同時(shí)堆芯相對(duì)體積發(fā)生變化,具有熱膨脹自熄滅機(jī)制。脈沖過(guò)程涉及了反應(yīng)堆動(dòng)力學(xué)、非穩(wěn)態(tài)傳熱、力學(xué)分析及耦合過(guò)程。因此,開(kāi)展快中子脈沖堆的瞬態(tài)分析,不僅可以驗(yàn)證中子學(xué)計(jì)算的精度,而且可對(duì)核-熱-力耦合計(jì)算及力學(xué)變形反饋模型進(jìn)行有效的驗(yàn)證[3]。本文開(kāi)展快譜小堆堆芯三維瞬態(tài)核-熱-力耦合計(jì)算方法初步研究,計(jì)算快中子脈沖堆Godiva的瞬發(fā)超臨界瞬態(tài)過(guò)程。

1 多物理耦合計(jì)算模型

1.1 傳熱計(jì)算模型

對(duì)于傳熱計(jì)算,可用熱傳導(dǎo)方程表示:

(1)

式中:ρ為密度;cp為比定壓熱容;T為溫度;k為導(dǎo)熱系數(shù);qf為體積釋熱率。

其中熱源由中子學(xué)模型計(jì)算得到:

(2)

式中:G為中子能群數(shù);φg為第g群中子通量密度;κΣf,g為第g群中子的能量產(chǎn)生截面。

1.2 力學(xué)計(jì)算模型

(3)

式中:u為位移;σ為應(yīng)力;f為外力。

在本文中采用小變形假設(shè),幾何方程可表示為:

(4)

式中,ε為應(yīng)變。

溫度變化也能引起變形,由于溫度變化引起的膨脹或收縮受到約束產(chǎn)生的應(yīng)力就是熱應(yīng)力,當(dāng)考慮溫度變化時(shí),根據(jù)Duhamel-Neumann原理,實(shí)際應(yīng)力需在彈性應(yīng)力中加入溫度載荷產(chǎn)生的附加應(yīng)力。若初始狀態(tài)下,結(jié)構(gòu)內(nèi)部初始應(yīng)力為零,此時(shí)結(jié)構(gòu)的本構(gòu)關(guān)系可表示為:

σij=Cijklεkl+βTΔTδij

(5)

式中:βT為熱膨脹系數(shù);ΔT為溫度的變化。

1.3 中子動(dòng)力學(xué)計(jì)算模型

在快中子脈沖堆Godiva中,由于采用高濃度的235U金屬核燃料,中子能譜很硬,同時(shí)金屬系統(tǒng)的溫升相對(duì)較小,相關(guān)研究表明多普勒效應(yīng)對(duì)反應(yīng)性的影響可忽略[3-5]。主要的反饋效應(yīng)是由于體積變形導(dǎo)致的中子泄漏引起的,因此必須考慮變形對(duì)中子學(xué)的影響。在本文中動(dòng)力學(xué)計(jì)算采用小變形假設(shè),則計(jì)算網(wǎng)格的體積應(yīng)變可表示為:

(6)

式中:V為網(wǎng)格變形后的體積;V0為初始網(wǎng)格體積。

與其他類(lèi)型的龍蝦相比，淡水龍蝦出現(xiàn)病害的情況比較少見(jiàn)。但是，如果稻田的水質(zhì)較差或者飼養(yǎng)管理措施不合理，那么出現(xiàn)病蟲(chóng)害的可能性就會(huì)大大提高，為了最大限度地降低龍蝦病害發(fā)生率，一方面要加強(qiáng)日常管理力度，另一方面要緊跟時(shí)代發(fā)展潮流，對(duì)龍蝦病害問(wèn)題進(jìn)行預(yù)防，在施肥過(guò)程中盡量采用少量多次的方式，提高龍蝦的抵抗力，保障龍蝦健康。

假設(shè)在瞬態(tài)過(guò)程中,中子的能譜不變,則中子宏觀反應(yīng)截面的變化主要是由計(jì)算網(wǎng)格體積改變導(dǎo)致的原子核密度改變引起的,瞬態(tài)過(guò)程中中子截面更新計(jì)算方法如下式所示:

(7)

式中:Σ為網(wǎng)格變形后的中子宏觀反應(yīng)截面;Σ0為初始的中子宏觀反應(yīng)截面。

對(duì)于中子學(xué)計(jì)算,本文采用中子擴(kuò)散時(shí)空動(dòng)力學(xué)進(jìn)行模擬:

(8)

式中:vg為第g群中子平均速度;φg為第g群中子通量密度;χp,g為瞬發(fā)中子到第g群中子的概率;χd,i,g為第i組緩發(fā)中子先驅(qū)核衰變后的緩發(fā)中子出現(xiàn)在第g群中的概率;λi為第i組緩發(fā)中子先驅(qū)核的衰變常量;βi為第i組緩發(fā)中子所占的份額;I為緩發(fā)中子先驅(qū)核總的組數(shù)。

1.4 中子動(dòng)力學(xué)計(jì)算模型耦合計(jì)算模型

方程(1)、(3)和(8)及相應(yīng)的初始條件、邊界條件即構(gòu)成了核-熱-力耦合的方程組,聯(lián)立求解該方程組即可得到中子通量、功率、溫度、應(yīng)力及位移等物理量的分布。在本文中空間離散采用有限元方法,對(duì)于時(shí)間離散,中子擴(kuò)散時(shí)空動(dòng)力學(xué)方程和瞬態(tài)傳熱方程均采用全隱式向后差分進(jìn)行離散,動(dòng)力方程采用直接積分法中的Newmark隱式算法[6]。經(jīng)過(guò)空間和時(shí)間離散,則可將描述核-熱-力的偏微分方程組轉(zhuǎn)換為非線(xiàn)性方程組:

(9)

式中:φ為中子通量密度的節(jié)點(diǎn)參量;T為溫度的節(jié)點(diǎn)參量;u為位移的節(jié)點(diǎn)參量;Kij和P為與中子通量、溫度和位移相關(guān)的耦合系數(shù)。

方程(9)可直接聯(lián)立求解,采用Picard或者Newton迭代法處理方程的非線(xiàn)性問(wèn)題,但直接聯(lián)立求解形成的矩陣規(guī)模很大,對(duì)計(jì)算機(jī)硬件要求高。本文首先采用算子分裂法,將中子場(chǎng)的求解分離出來(lái),則可將核-熱-力耦合問(wèn)題轉(zhuǎn)換為順序求解熱彈性動(dòng)力學(xué)方程組(式(10))和中子動(dòng)力學(xué)方程(式(11))。

(10)

(11)

采用JFNK方法(Jacobian-Free Newton Krylov method)求解熱彈性動(dòng)力學(xué)方程組(10),對(duì)于非線(xiàn)性方程組(10)將右端項(xiàng)移到左邊,并寫(xiě)為算子形式:

F(X)=0 (X=[T,u])

(12)

采用牛頓法求解方程(式(12)),可將非線(xiàn)性問(wèn)題轉(zhuǎn)換為迭代求解如下的線(xiàn)性問(wèn)題:

J(Xk)ΔXk+1=-F(Xk)

Xk+1=Xk+ΔXk+1

(13)

JFNK方法在外層采用Newton-Raphson方法處理偏微分方程組的非線(xiàn)性問(wèn)題,內(nèi)層采用克雷洛夫子空間法求解線(xiàn)性方程[7]。JFNK方法區(qū)別與傳統(tǒng)Newton方法的是由于內(nèi)層采用克雷洛夫子空間法,因此不需要顯式的構(gòu)造Jacobi矩陣,只需要Jacobi矩陣和向量的乘積,可采用差分近似方法計(jì)算得到,如式(14)所示。JFNK方法的求解效率與預(yù)處理方法相關(guān),Krylov迭代中的預(yù)處理方法很多,包括應(yīng)用于一般矩陣的代數(shù)預(yù)處理子到利用特殊問(wèn)題類(lèi)之特征的問(wèn)題相關(guān)預(yù)處理子。本文中非線(xiàn)性方程的數(shù)值求解都是基于petsc程序包,計(jì)算中采用的預(yù)處理方法為代數(shù)多重網(wǎng)格方法[8-9]。

(14)

式中,ε為差分步長(zhǎng),計(jì)算方式如文獻(xiàn)[9]所示。

瞬態(tài)過(guò)程中每一時(shí)間步的計(jì)算流程如圖1所示。圖1中:1) 根據(jù)上一步中子學(xué)計(jì)算得到的中子通量密度進(jìn)行裂變熱源計(jì)算;2) 采用JFNK方法求解熱彈性動(dòng)力學(xué)方程組;3) 根據(jù)熱彈性動(dòng)力學(xué)方程組求得的節(jié)點(diǎn)位移,計(jì)算每個(gè)網(wǎng)格的體積應(yīng)變,根據(jù)式(7)更新計(jì)算中采用的宏觀截面,同時(shí)將計(jì)算得到的節(jié)點(diǎn)位移傳遞給中子動(dòng)力學(xué)計(jì)算的網(wǎng)格進(jìn)行變形;4) 根據(jù)更新的宏觀截面和網(wǎng)格進(jìn)行中子動(dòng)力學(xué)計(jì)算,得到通量分布。

2 基于MOOSE的多物理耦合程序開(kāi)發(fā)與驗(yàn)證

本文中,三維核-熱-力耦合求解程序是基于MOOSE多物理耦合框架開(kāi)發(fā)的。MOOSE框架以面向?qū)ο蟮木幊趟枷脒M(jìn)行開(kāi)發(fā),對(duì)物理場(chǎng)偏微分方程中的各數(shù)學(xué)算子進(jìn)行抽象封裝,通過(guò)平臺(tái)內(nèi)不同的程序模塊組合實(shí)現(xiàn)對(duì)某個(gè)物理場(chǎng)具體偏微分方程的描述。具備以下特征:1) 與維度無(wú)關(guān)的用戶(hù)代碼,一套代碼對(duì)一、二、三維問(wèn)題都適用;2) 基于連續(xù)或非連續(xù)Galerkin有限元的網(wǎng)格離散;3) 全隱式耦合;4) 允許采用非結(jié)構(gòu)化網(wǎng)格,包含多種形狀函數(shù);5) 很強(qiáng)的網(wǎng)格自適應(yīng)性;6) 內(nèi)置并行功能,用戶(hù)不用開(kāi)發(fā)額外的并行代碼;7) 內(nèi)置的后處理等[10]。

作者基于MOOSE開(kāi)發(fā)了基于中子擴(kuò)散理論和偶階輸運(yùn)理論的中子學(xué)計(jì)算模塊,該中子學(xué)求解模塊經(jīng)過(guò)一系列基準(zhǔn)題的驗(yàn)證[11]。傳熱和動(dòng)力學(xué)方程求解是基于MOOSE中內(nèi)置的HeatConduction模塊和Tensormechanics模塊,在本文中將這兩個(gè)模塊耦合用以求解熱彈性問(wèn)題。為驗(yàn)證傳熱和力學(xué)模塊的正確性,本文選取勞倫斯利弗莫爾國(guó)家實(shí)驗(yàn)室(LANL)構(gòu)建的金屬鈾-鉬球殼堆[12]進(jìn)行熱應(yīng)力分析計(jì)算,其中球殼的內(nèi)徑為5.08 cm,外徑為7.62 cm。文獻(xiàn)中給出的熱加載關(guān)系為:

(15)

式中:Tmax為最大溫升;b為脈沖半高寬;tpp為峰功率時(shí)刻。

在絕熱邊界條件下,計(jì)算脈沖半高寬為41 μs,堆芯平均溫升為420 ℃,堆芯內(nèi)外表面位移計(jì)算結(jié)果如圖2所示。本文計(jì)算得到的結(jié)果與文獻(xiàn)中給出的解析解吻合良好,說(shuō)明程序可準(zhǔn)確模擬脈沖堆熱-力耦合過(guò)程。

中子學(xué)求解和熱彈性問(wèn)題的耦合求解是通過(guò)MOOSE的MultiApp和Transfer功能實(shí)現(xiàn)的。在計(jì)算中熱力學(xué)求解程序作為MaterApp,中子學(xué)求解程序作為MultiApp,通過(guò)Transfer功能實(shí)現(xiàn)功率和位移在MaterApp和MultiApp之間的傳遞。

3 快中子脈沖堆Godiva-Ⅰ三維核-熱-力耦合瞬態(tài)計(jì)算研究

Godiva-Ⅰ是由LANL建造的世界上第1座快中子脈沖堆。Godiva-Ⅰ是LANL為驗(yàn)證熱膨脹自熄滅裂變脈沖的機(jī)理,將臨界裝置Lady Godiva改建成快中子脈沖堆。它是一個(gè)無(wú)反射層的高富集鈾金屬球形快中子脈沖堆。該裝置的堆芯是由235U富集度為93.5%、密度18.75 g/cm3的鑄造鈾金屬構(gòu)成,緩發(fā)臨界質(zhì)量約為52.04 kg,堆芯金屬球直徑約為17 cm[13]。本文中對(duì)Godiva-Ⅰ初始周期為16.2 μs的瞬態(tài)過(guò)程進(jìn)行了模擬計(jì)算。

3.1 計(jì)算模型

在本文的計(jì)算中,采用三維幾何進(jìn)行建模,網(wǎng)格如圖3所示。計(jì)算中采用的傳熱和力學(xué)參數(shù)列于表1。傳熱計(jì)算中,金屬球體表面設(shè)置為絕熱邊界條件,溫度場(chǎng)初始條件設(shè)置為均勻溫度場(chǎng),溫度為293.6 K。動(dòng)力學(xué)計(jì)算中,金屬球體表面設(shè)置為自由膨脹邊界條件,初始位移和速度都設(shè)置為0[14-15]。

表1 Godiva-Ⅰ基本參數(shù)Table 1 Basic parameter of Godiva-Ⅰ

圖3 Godiva-Ⅰ網(wǎng)格劃分模型Fig.3 Meshing model of Godiva-Ⅰ

多群中子截面和緩發(fā)中子先驅(qū)核參數(shù)是采用Serpent程序在293.6 K的溫度下計(jì)算制作得到的[16]。中子學(xué)計(jì)算采用的能群結(jié)構(gòu)列于表2[17]。中子學(xué)計(jì)算先求解穩(wěn)態(tài)特征值問(wèn)題,并將功率歸一化為500 W,作為瞬態(tài)計(jì)算的輸入條件。瞬態(tài)計(jì)算的持續(xù)時(shí)間為0.001 s,時(shí)間步長(zhǎng)為1×10-7s。

表2 能群結(jié)構(gòu)Table 2 Group structure

3.2 脈沖瞬態(tài)過(guò)程計(jì)算結(jié)果

使用編寫(xiě)的多物理耦合程序?qū)odiva-Ⅰ初始周期為16.2 μs的瞬發(fā)超臨界瞬態(tài)過(guò)程進(jìn)行了計(jì)算,計(jì)算結(jié)果如圖4所示。從圖4可看出,程序計(jì)算值和實(shí)驗(yàn)值吻合良好。

圖4 瞬態(tài)過(guò)程中裂變率計(jì)算值和實(shí)驗(yàn)值的比較Fig.4 Comparison between numerical and experimental results of fission rate

瞬態(tài)過(guò)程中平均溫升隨時(shí)間的變化如圖5所示。由圖5可看出,引入正反應(yīng)性后,在開(kāi)始階段,由于系統(tǒng)功率上升較慢,平均溫度上升也較小,隨后堆芯功率急劇上升,累積的熱量導(dǎo)致金屬鈾的溫度升高,熱膨脹導(dǎo)致球體發(fā)生膨脹,產(chǎn)生位移;隨后由于堆芯膨脹,金屬密度減少,中子泄漏增強(qiáng),引入負(fù)的反應(yīng)性,堆芯處于次臨界狀態(tài),導(dǎo)致堆芯功率迅速降低,堆芯平均溫度變化逐漸減緩。圖6為外表面位移和速度。由圖6可發(fā)現(xiàn),由于裂變功率的迅速升高和降低,堆芯震蕩,導(dǎo)致堆芯相繼處于壓縮和膨脹的狀態(tài),產(chǎn)生的熱慣性效應(yīng),導(dǎo)致堆芯表面位移、速度和加速度產(chǎn)生震蕩,同時(shí)由于變形導(dǎo)致堆芯反應(yīng)性的改變,堆芯功率在500 μs后也產(chǎn)生周期性的震蕩,如圖7所示。在瞬態(tài)計(jì)算結(jié)果中還觀察到堆芯功率分布形狀變化較小,因此傳統(tǒng)的基于點(diǎn)堆動(dòng)力學(xué)的中子學(xué)分析方法也是一種誤差較小的近似方法。

圖5 平均溫升隨時(shí)間的變化Fig.5 Time evolution of average temperature increase

圖6 外表面位移和速度隨時(shí)間的變化Fig.6 Time evolution of outer surface displacement and velocity

圖7 裂變率隨時(shí)間的變化Fig.7 Time evolution of fission rate

在圖1中可看出,本文建立的多物理耦合模型,力學(xué)變形效應(yīng)的反饋?zhàn)饔冒▋蓚€(gè)方面:1) 熱膨脹導(dǎo)致堆芯體積增大,核子密度減小,引起宏觀反應(yīng)截面的變化,中子泄漏增強(qiáng),這是一種負(fù)反饋效應(yīng);2) 堆芯體積增大,導(dǎo)致堆芯尺寸的變化,這是一種正的反饋?zhàn)饔谩Ｔ趥鹘y(tǒng)的快堆分析中,一般認(rèn)為第1種反饋是主導(dǎo)因素,第2種是二階效應(yīng),在小變形條件下,一般忽略堆芯尺寸變化引起的反應(yīng)性反饋效應(yīng)。本文分析了Godiva-Ⅰ在瞬發(fā)超臨界瞬態(tài)條件下,變形導(dǎo)致的密度變化和幾何形狀變化對(duì)瞬態(tài)計(jì)算結(jié)果的影響。圖8中的兩個(gè)算例都考慮了密度負(fù)反饋效應(yīng),區(qū)別是是否考慮堆芯尺寸變化的反饋效應(yīng)。從圖8可看出,采用動(dòng)網(wǎng)格對(duì)幾何進(jìn)行變形,幾何尺寸增加引入的是正的反應(yīng)性,會(huì)使峰值功率增加,而只考慮密度變化反饋效應(yīng)得到的計(jì)算值偏低。

總之,在Godiva瞬發(fā)超臨界瞬態(tài)計(jì)算過(guò)程中,由于多普勒反饋效應(yīng)是可忽略的,熱載荷導(dǎo)致的力學(xué)變形效應(yīng)是計(jì)算中唯一考慮的反應(yīng)性反饋機(jī)制。從計(jì)算結(jié)果可看出,本文的方法可準(zhǔn)確考慮由于熱膨脹導(dǎo)致的反應(yīng)性反饋效應(yīng)。

4 結(jié)論

本文基于MOOSE多物理耦合框架,開(kāi)展了三維核-熱-力耦合計(jì)算研究,開(kāi)發(fā)了多物理耦合計(jì)算程序,并應(yīng)用該程序計(jì)算了快中子脈沖堆Godiva-Ⅰ的瞬發(fā)超臨界瞬態(tài),實(shí)現(xiàn)了瞬態(tài)過(guò)程中,裂變率、溫度、位移等物理量的準(zhǔn)確模擬,計(jì)算結(jié)果與實(shí)驗(yàn)值符合較好。計(jì)算結(jié)果表明,該計(jì)算模型可準(zhǔn)確模擬熱膨脹引起的反應(yīng)性反饋效應(yīng),為進(jìn)一步開(kāi)展快譜小堆堆芯多物理耦合安全分析奠定了基礎(chǔ)。