基于格子Boltzmann方法的反應堆核-熱-流耦合模擬

王亞輝，馬 宇

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

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

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

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

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

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

(1)

(2)

(3)

(4)

(5)

(6)

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

為實現堆芯核-熱-流耦合過程的高效準確模擬，本文采用了具有多物理耦合特性的LBM進行求解，該方法能將不同的物理過程轉化為相似的LBM形式并在相似的數據結構以及離散格式下采用相似的求解格式實現耦合求解。對堆芯核-熱-流耦合過程，采用如下LBM模型：

(7)

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

a——D2Q4；b——D2Q9

D2Q4離散速度定義為：

(8)

D2Q9離散速度定義為：

e=

(9)

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

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

(10)

(11)

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

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

(12)

(13)

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

(14)

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

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

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

1)反彈邊界

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

H3(B)=H1(B) D2Q4

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

(15)

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

2)非平衡外推邊界

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

(16)

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

(17)

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

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

2 數值實驗

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

2.1 LBM模型驗證

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

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

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

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

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

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

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

q?=γ(Tref-T)

(18)

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

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

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

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

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

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

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

3 結論

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

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