王亞輝，馬 宇

(中山大學(xué) 中法核工程與技術(shù)學(xué)院，廣東 珠海 519082)

核反應(yīng)堆堆芯是一個典型的多物理耦合系統(tǒng)，其中的核-熱-流耦合過程是反應(yīng)堆工程模擬中的一個主要研究點[1-2]，尤其是對于反應(yīng)堆設(shè)計及安全分析。中子輸運(yùn)過程決定了核燃料中中子裂變產(chǎn)生的能量，而熱工水力過程決定了溫度分布，并通過多普勒效應(yīng)影響宏觀中子截面。由于這一特點，核反應(yīng)堆行為的可靠分析需要涉及中子輸運(yùn)和熱工水力過程的多物理耦合[3-4]。從數(shù)學(xué)角度來看，這些多物理模擬的實現(xiàn)依賴于在不同物理場之間有效傳輸數(shù)據(jù)的能力。因此，在相同方法下解決耦合問題，可有效消除外部數(shù)據(jù)交換和插值。然而，由于不同物理場控制方程之間的差異性以及中子輸運(yùn)方程的復(fù)雜性，對其耦合求解存在一定的難度。當(dāng)前的研究工作普遍基于不同數(shù)值方法對不同物理場進(jìn)行求解，隨后在不同物理場之間增加外部數(shù)據(jù)接口和插值過程，這使得多物理耦合實現(xiàn)困難，并且可能會降低精度[5]。

格子Boltzmann方法(LBM)[6-9]是一種起源于流體動力學(xué)求解的數(shù)值方法，該方法采用介觀分布函數(shù)來求解宏觀變量的演化過程，可將不同的控制方程轉(zhuǎn)換為相似的線性LBM方程，采用相似的格式進(jìn)行求解，這使得LBM在多物理耦合模擬中具有突出優(yōu)勢。相比于傳統(tǒng)數(shù)值方法而言，LBM具有極為簡單的形式，只需上百行程序即可實現(xiàn)復(fù)雜物理過程的模擬；同時由于LBM方程只與當(dāng)前的節(jié)點以及上游的一個節(jié)點有關(guān)，該方法具有極強(qiáng)的并行特性?；谝陨蟽?yōu)勢，LBM在流動傳熱[10-13]、多相流[14-15]、流固耦合[16]、多孔介質(zhì)流[17-19]、磁流體[20-21]、化學(xué)反應(yīng)[22]、堆芯物理[5，23-26]等方面受到了大量的關(guān)注。

本文建立求解反應(yīng)堆堆芯核-熱-流耦合過程的LBM模型，將核-熱-流不同物理過程統(tǒng)一到相似的LBM格式下采用相似的數(shù)據(jù)結(jié)構(gòu)和離散格式進(jìn)行求解，消除不同物理場之間的外部數(shù)據(jù)傳輸和插值過程，降低多物理耦合求解的實現(xiàn)難度。

1 基于LBM的核-熱-流耦合模擬

反應(yīng)堆堆芯內(nèi)部存在著多種物理過程的相互耦合，工程中常?？紤]核-熱-流耦合過程。對于燃料區(qū)域，溫度分布通過改變中子截面進(jìn)而影響中子物理過程，同時裂變產(chǎn)熱又作為熱源影響傳熱過程。對于慢化劑區(qū)域，除以上效應(yīng)外，溫度分布還受到流動過程的影響?？傮w而言，核-熱-流耦合過程可以采用如下控制方程組進(jìn)行描述，包括中子輸運(yùn)方程、緩發(fā)中子先驅(qū)核守恒方程、大渦模擬(LES)方程以及能量守恒方程：

(1)

(2)

(3)

(4)

(5)

(6)

其中，Ef為1次裂變釋放的能量。式(2)左側(cè)第1項(對流項)主要存在于液體燃料堆芯中(如熔鹽堆)，是緩發(fā)中子隨液體燃料流動遷移而產(chǎn)生的。在核-熱-流耦合計算中，中子宏觀截面Σx與溫度T相關(guān)，即Σx=Σx(T)。

為實現(xiàn)堆芯核-熱-流耦合過程的高效準(zhǔn)確模擬，本文采用了具有多物理耦合特性的LBM進(jìn)行求解，該方法能將不同的物理過程轉(zhuǎn)化為相似的LBM形式并在相似的數(shù)據(jù)結(jié)構(gòu)以及離散格式下采用相似的求解格式實現(xiàn)耦合求解。對堆芯核-熱-流耦合過程，采用如下LBM模型：

(7)

二維條件下典型的離散速度模型有D2Q4(2維4速)和D2Q9(2維9速)模型，如圖1所示。

a——D2Q4；b——D2Q9

D2Q4離散速度定義為：

(8)

D2Q9離散速度定義為：

e=

(9)

其中，c為格子遷移速度，c=δx/δt，δx為格子寬度，即網(wǎng)格尺寸。

該模型的平衡態(tài)分布函數(shù)和無量綱弛豫時間可表示為：

(10)

(11)

其中：A、B、C為物理場參數(shù)；?為LBM格子方向權(quán)重；Ψ為廣義宏觀變量；ζ和U分別為廣義擴(kuò)散系數(shù)及廣義宏觀速度；I為單位對角矩陣；τf為為廣義無量綱弛豫時間；cs為LBM格子聲速。

在統(tǒng)一LBM框架下，所有物理場可采用相似的過程進(jìn)行演化計算，包括碰撞和遷移過程如下：

(12)

(13)

其中，H*為碰撞后分布函數(shù)。宏觀變量可通過對所有格子方向的分布函數(shù)求格子速度的0階矩得到：

(14)

通過對方程中的變量賦予不同的定義，以上的統(tǒng)一LBM模型可被用來求解不同的問題。本研究中，對中子輸運(yùn)問題的模擬按照精度需求的不同可選擇中子輸運(yùn)SN方程LBM模型、SP3方程LBM模型以及擴(kuò)散方程LBM模型。為了使該核-熱-流耦合LBM模型能用于固體燃料反應(yīng)堆及液體燃料反應(yīng)堆，本研究中對考慮流動效應(yīng)的緩發(fā)中子先驅(qū)核守恒方程采用有限Boltzmann形式的LBM模型。

基于中子擴(kuò)散方程和中子輸運(yùn)SP3方程的特性，二者求解的LBM模型均采用D2Q4格子速度模型。由于緩發(fā)中子方程和SN方程均具有強(qiáng)對流特性，這兩個LBM模型均采用D2Q9格子速度模型。流動速度場的LBM模型也采用D2Q9格子速度進(jìn)行模擬。對于傳熱方程LBM模型，可以與LES方程的LBM模型類似地采用D2Q9格子速度模型，但對于不考慮流動效應(yīng)的中子輸運(yùn)-傳熱耦合過程，傳熱方程LBM模型可采用更高效且穩(wěn)定的D2Q4格子速度模型進(jìn)行模擬。

本文給出核-熱-流耦合過程中幾種典型的邊界條件處理方式，包括中子輸運(yùn)過程的真空邊界、反射邊界、流動傳熱過程中的對稱邊界、絕熱邊界、無滑移邊界、自由出流邊界等。從LBM的邊界處理角度來說，這些邊界可統(tǒng)一用兩種方式處理，包括反彈(BB)邊界及非平衡外推邊界。

1)反彈邊界

反彈邊界的定義是在邊界上出射的分布函數(shù)經(jīng)過邊界反彈后，沿其鏡像對稱方向再次反彈回計算域中，這種邊界一般被用于反射邊界、對稱邊界及絕熱邊界。以右邊界為例，此時對稱邊界的定義如下：

H3(B)=H1(B) D2Q4

H3,6,7(B)=H1,5,8(B) D2Q9

(15)

其中，B為邊界節(jié)點。

2)非平衡外推邊界

非平衡外推邊界基于這樣一個假設(shè)：節(jié)點上的分布函數(shù)Fi可表示為平衡態(tài)部分和非平衡態(tài)部分的和：

(16)

在非平衡外推邊界中，邊界節(jié)點的平衡態(tài)分布函數(shù)可由式(10)直接給出，而假設(shè)邊界節(jié)點的非平衡部分可由其內(nèi)一層節(jié)點的非平衡態(tài)部分外推得到：

(17)

其中，O為邊界向內(nèi)一層節(jié)點。

由于非平衡外推邊界采用外推方式，理論上該方法可適用于所有類型的邊界條件處理，包括前文中的反射邊界、對稱邊界及絕熱邊界等。

2 數(shù)值實驗

選取典型問題對本文建立的核-熱-流耦合LBM模型進(jìn)行測試與驗證。在流場計算中，Smagorinsky系數(shù)設(shè)為Cs=0.1；所有模擬的誤差限均取為10-6。由于中子輸運(yùn)LBM模型(包括SN方程、SP3方程以及擴(kuò)散方程)及核-熱耦合LBM模型均已在文獻(xiàn)[27]中進(jìn)行過驗證，本文不再重復(fù)，只著重對基于LES方程的流場流動LBM進(jìn)行驗證。

2.1 LBM模型驗證

為驗證核-熱-流耦合LBM在不同流態(tài)流場模擬的準(zhǔn)確性，對典型的方腔驅(qū)動流進(jìn)行了模擬研究[28-30]。為研究耦合LBM模型對不同Re的適應(yīng)性，分別對Re=100、1 000以及10 000的流場進(jìn)行了模擬。當(dāng)Re=100和1 000時，分別采用128×128網(wǎng)格，當(dāng)Re=10 000時，由于湍流流態(tài)的復(fù)雜性，采用256×256網(wǎng)格進(jìn)行模擬，格子速度模型統(tǒng)一采用D2Q9模型。方腔的高和寬分別取為H和L，流體驅(qū)動速度為U0。

沿水平中線的橫向相對速度分量以及沿垂直中線的縱向相對速度分量如圖2所示，并與直接數(shù)值模擬(DNS)方法進(jìn)行對比。由圖2可看出，本文核-熱-流耦合LBM的模擬結(jié)果在不同Re下與參考解均吻合得很好，證明當(dāng)前的耦合LBM模型對層流和湍流條件均有較強(qiáng)的適應(yīng)性及較高的精度。

a——水平中線橫向相對速度分量u/U0；b——垂直中線縱向相對速度分量v/U0

2.2 基于LBM模型的簡化熔鹽堆核-熱-流耦合分析

液體燃料熔鹽堆(MSR)是第4代反應(yīng)堆中唯一一種液體燃料堆芯，該堆芯能在高溫條件下保持較低的蒸汽壓，具有較高的安全性，同時可實現(xiàn)小型化應(yīng)用。不同于固體燃料堆芯，MSR的核燃料溶解在熔鹽慢化劑中，隨著熔融鹽在堆芯中流動而遷移，熔融鹽的流動效應(yīng)對緩發(fā)中子先驅(qū)核的遷移行為產(chǎn)生影響，同時緩發(fā)中子先驅(qū)核的遷移影響中子注量率和溫度分布，而溫度分布又通過中子截面和材料熱物性影響中子輸運(yùn)和流動過程。因此，對于MSR而言，中子輸運(yùn)-傳熱-流動過程是強(qiáng)烈耦合在一起的，必須對其進(jìn)行耦合研究。

為研究耦合LBM模型對MSR模擬的適應(yīng)性及MSR內(nèi)部的多物理耦合行為，基于本文建立的耦合LBM模型對簡化MSR內(nèi)部的核-熱-流耦合過程進(jìn)行模擬分析。基于相關(guān)文獻(xiàn)[31]建立了正方形簡化MSR模型，其內(nèi)部熔鹽流動被作為自由流動處理。反應(yīng)堆核-熱耦合過程的準(zhǔn)確性已在之前的研究中進(jìn)行了充分驗證，同時本文建立的LBM模型對于流動復(fù)雜過程模擬的準(zhǔn)確性也已進(jìn)行了驗證，為此本文不再進(jìn)行重復(fù)驗證。簡化堆芯結(jié)構(gòu)為2 m×2 m的方形區(qū)域，其內(nèi)部存在自由流動的液體熔鹽，核燃料被溶解在液體熔鹽中，隨液體熔鹽流動而流動。該問題考慮雙群6組緩發(fā)中子，模擬過程中考慮中子物性及材料熱物性的溫度反饋。

對該問題考慮中子輸運(yùn)-傳熱-流動耦合過程模擬，不考慮熔融鹽壓縮效應(yīng)，并且使用Boussinesq方程近似熔鹽流動過程中的浮升力隨溫度的變化。該堆芯是一個均勻液體反應(yīng)堆，四周邊界均為真空無入射中子邊界，緩發(fā)中子先驅(qū)核在壁面處均被反射回堆芯，沒有泄漏。熔鹽在壁面處的流動邊界均為無滑移邊界，同時左右壁面及下壁面均為靜止壁面，上壁面以Ud的速度驅(qū)動堆芯內(nèi)熔鹽流動。堆芯內(nèi)熱量不通過壁面?zhèn)鬟f出堆芯，即四壁面均為絕熱邊界。熔鹽堆芯熱量導(dǎo)出通過堆芯內(nèi)的散熱器實現(xiàn)，散熱器散熱過程近似可表示為：

q?=γ(Tref-T)

(18)

其中：Tref為參考溫度，Tref=900 K；γ為均勻體積換熱系數(shù)，γ=106W/(m3·K)；q?為散熱器換熱功率密度。

分析了不同條件下的堆芯內(nèi)耦合行為，分析條件分別為：中子截面無溫度反饋(Ud=0.5 m/s)、中子截面有溫度反饋(Ud=0.5 m/s，5 m/s)。圖3示出不同條件下瞬發(fā)中子注量率φ1、緩發(fā)中子先驅(qū)核C1以及功率P的分布，圖4示出堆芯流線圖。圖3中flagT=0表示無溫度反饋，flagT=1表示有溫度反饋。由圖3可看出，溫度反饋強(qiáng)烈地影響中子注量率分布和功率分布，同時對于緩發(fā)中子也有較為明顯的影響。

當(dāng)不考慮溫度反饋條件時，由于堆芯內(nèi)部中子截面分布均勻，此時中子注量率分布從堆芯中心向堆芯邊緣逐漸降低，中子注量率最大的區(qū)域集中在中心區(qū)域附近，因此此時堆芯功率的最大值也集中在堆芯中部，如圖3a、c實線所示。由于低速下熔鹽內(nèi)部浮升力起主要作用，緩發(fā)中子將隨著熔鹽向堆芯上部遷移，緩發(fā)中子密度在堆芯上部較大，如圖3b實線所示。

如圖3中虛線所示，當(dāng)考慮溫度反饋時，堆芯中心處中子注量率更大更加集中，同時中子注量率分布和功率分布集中在堆芯偏下區(qū)域。這一現(xiàn)象是由緩發(fā)中子先驅(qū)核和溫度分布共同直接作用影響的。一方面，由溫度反饋引起的宏觀截面的變化直接影響了中子輸運(yùn)過程，間接導(dǎo)致堆芯下部中子注量率分布集中。當(dāng)Ud=0.5 m/s時，浮升力和對流傳熱同時起作用。在相對較強(qiáng)的浮升力作用下，由于密度降低，溫度較高的燃料隨熔鹽材料向上流動，如圖4a所示。因此，高溫流體集中在堆芯上部，使得上部區(qū)域中的宏觀截面小于下部區(qū)域的宏觀截面。由于燃料在下部區(qū)域具有更強(qiáng)的慢化和裂變效應(yīng)，因此中子分布更集中在該區(qū)域，同時相應(yīng)地功率分布也更集中于該區(qū)域。中子在下部的積聚導(dǎo)致了緩發(fā)中子的積聚，同時由于燃料與熔融鹽一起流動，所以對流效應(yīng)使得緩發(fā)中子先驅(qū)核向堆芯的下部蓄積，這增加了堆芯下部的緩發(fā)中子源，也相應(yīng)地增加了瞬發(fā)中子在下部的積聚。

圖3 不同條件下快群中子(a)、緩發(fā)中子先驅(qū)核(b)及功率(c)分布對比

圖4 低速(a)和高速(b)條件下堆芯流線圖

圖3中點劃線展示了流速對于MSR堆芯工況的影響，即對流傳熱強(qiáng)度的影響。相較于虛線情況而言，點劃線所帶有的條件具有較高的驅(qū)動速度，并因此具有較高的對流傳熱效果。在研究溫度和流動條件時，可發(fā)現(xiàn)當(dāng)Ud=5 m/s時浮升力對流體流動的影響相對較小，此時流體呈現(xiàn)強(qiáng)迫對流狀態(tài)，如圖4b所示。在這種情況下，由于浮升力導(dǎo)致的左側(cè)渦結(jié)構(gòu)變小，而由于驅(qū)動力導(dǎo)致的右側(cè)渦結(jié)構(gòu)主導(dǎo)了流場。在這種條件下，不同溫度的流體通過流體流動相互混合變得均勻，并且由浮升力引起的高溫流體的上浮過程被減弱。主流區(qū)域中相對均勻的溫度分布導(dǎo)致該區(qū)域中子截面更均勻。因此，中子分布更集中在堆芯的中心區(qū)域，堆芯功率也集中在靠近中心區(qū)域。

3 結(jié)論

本文針對反應(yīng)堆堆芯內(nèi)中子輸運(yùn)-傳熱-流動多物理耦合過程，構(gòu)建了核-熱-流耦合過程的多物理場耦合LBM模型，其中中子輸運(yùn)過程涵蓋了SN方程、SP3方程以及擴(kuò)散方程3種不同近似；緩發(fā)中子先驅(qū)核守恒方程考慮了MSR堆芯中燃料流動效應(yīng)，使得該模型對于不同尺度、不同特點的堆芯均具有較強(qiáng)的適應(yīng)性?；贚BM的實現(xiàn)簡單性以及多物理耦合特性，降低了多物理耦合模擬的難度，同時消除了傳統(tǒng)外耦合方法的數(shù)據(jù)傳輸與插值過程。數(shù)值驗證結(jié)果表明，本文建立的耦合LBM模型對不同流態(tài)的流動過程均具有較強(qiáng)的適應(yīng)性與較高的精度?；隈詈螸BM模型對MSR內(nèi)部核-熱-流耦合行為的研究表明，溫度反饋對于高溫堆芯有較為明顯的作用，不可忽略；同時提高液體熔鹽流速能夠有效展平功率和溫度分布，改善傳熱。耦合過程的模擬結(jié)果表明本文的LBM模型對核-熱-流耦合過程具有良好的適應(yīng)性，同時能夠?qū)σ簯B(tài)燃料堆芯條件給出合理的研究結(jié)果。

本文研究是LBM在堆芯多物理耦合模擬領(lǐng)域的一個初步嘗試，在將來的工作中，將進(jìn)一步擴(kuò)展所建立的耦合LBM模型，包括多相流動過程、耦合燃耗過程、力學(xué)過程等，進(jìn)一步提高模型適用性以及物理場耦合的全面性。