基于CFD程序的物理熱工耦合計算不確定性分析

于 濤，雷洲陽，趙鵬程,*，錢冠華，李 捷

(1.南華大學 核科學技術學院，湖南 衡陽 421001；2.南華大學 湖南省數字化反應堆工程技術研究中心，湖南 衡陽 421001)

物理熱工耦合方法中的不確定性量化是保障物理熱工耦合計算結果可靠性的重要基礎，是物理熱工耦合的重要研究方向之一。不確定性分析主要采用輸入不確定性傳播方法，研究輸入參數對模型帶來的不確定性[1]。將不確定性分析方法應用到物理熱工耦合計算，能更精確地預估模型的參數范圍，為核反應堆設計、安全分析和評審提供可靠的計算置信區間，具有重要的學術意義和實用價值[2-3]。

國際上，由美國橡樹嶺國家實驗室(ORNL)牽頭的輕水堆先進仿真聯盟(CASL)于2011年起采用DAKOTA軟件開展對多物理耦合的不確定性分析，取得了較為豐富的成果[4-9]；美國能源部(DOE)核能辦公室成立的核能先進建模與仿真項目(NEAMS)于2017—2019年采用PROTEUS-Nek5000耦合工具箱SHARP開展了鈉冷快堆和鉛冷快堆安全特性的不確定性分析研究[10]；歐盟在“歐洲核能技術平臺可持續發展(SNETP)”戰略規劃中開展了NURESIM系列項目，在該項目中基于SALOME平臺進行了多物理、多尺度耦合相關研究，并采用URANIE工具箱開展了相關不確定性分析工作[11-13]。近年來，國內也進行了一系列不確定性研究：沈昊等[14]采用敏感性系數法對棒束單相流動傳熱計算流體力學(CFD)進行了不確定性量化分析；潘昕懌等[15]采用混合法和隨機抽樣法研究了多群核數據不確定性對堆芯物理計算的影響；高新力等[16]基于SNAP-DAKOTA-RELAP5程序開展了壓水堆大破口事故的不確定性分析；張梟羽等[17]采用多項式混沌法進行了壓水堆堆芯冷卻劑通道計算的不確定性分析；趙陽[18]基于DAKOTA/SALib程序開展了釷基熔鹽堆系統事故不確定性分析；曹志偉等[19]采用確定論統計法(GSM)對CPR1000核電廠大破口事故進行了不確定性量化分析；郭炯等[20]基于VSOP-UAM程序開展了高溫氣冷堆的不確定性分析；劉勇等[21]采用兩步法的敏感性計算策略，針對快堆BN-600基準例題的有效增殖因數進行了不確定性量化分析；鄧程程等[22]基于RELAP5最佳估算程序對先進熱工水力試驗(ACME)臺架開展了小破口失水事故的不確定性和敏感性分析。一方面，國內在堆物理、系統程序分析不確定性量化方面進行了大量工作，但有關物理熱工耦合不確定性分析工作開展較少；另一方面，第4代反應堆中鈉冷快堆以及鉛基快堆存在極為復雜的物理熱工耦合現象，應用不確定性分析可更精確地評估反應堆安全特性，因此有必要進行深入研究。

本文基于CFD程序FLUENT的用戶自定義函數(UDF)，耦合中子動力學計算模型、燃料棒熱傳導計算模型、不確定性分析程序SIMLAB，開發物理熱工耦合計算不確定性分析平臺CFD/PFS。進而基于CFD/PFS開展小型自然循環鉛冷快堆SNCLFR-10的無保護超功率(UTOP)事故不確定性分析，研究燃料多普勒系數等不確定性輸入參數對功率峰值等目標參數帶來的不確定性，并評估反應性反饋系數對目標參數的影響程度。

1 不確定性分析CFD/PFS平臺開發

1.1 點堆中子動力學模型

點堆中子動力學模型(PKM)主要用于計算反應堆內中子密度的動態變化，進而更新堆芯功率，為冷卻劑和慢化劑熱工水力計算提供熱源[23]。本工作采用6組緩發中子的點堆中子動力學模型：

(1)

(2)

其中：N(t)為中子密度；ρ(t)為反應性；Λ為中子代時間；U為緩發中子總份額；t為時間；Ui為第i組緩發中子份額；λi為第i組緩發中子先驅核衰變常量；Ci(t)為第i組緩發中子先驅核濃度；q為外加中子源源強。采用二階泰勒多項式積分法求解點堆中子動力學方程，將中子密度N(t)使用全隱式二階泰勒近似展開，以迭代的方式求解不同時刻的中子密度。

點堆中子動力學模型中，中子密度主要由總反應性控制，總反應性包括外部引入反應性和內部反應性反饋[24]。反饋反應性主要包括由燃料溫度引起的燃料多普勒反應性、材料密度變化引起的軸向膨脹反應性和徑向膨脹反應性以及冷卻劑密度反應性。總反應性為：

αc,r(Tc(t)-Tc(0))+αc(Tc(t)-Tc(0))

(3)

其中：KD為多普勒系數；αc,a為堆芯軸向膨脹反饋系數；αc,r為堆芯徑向膨脹反饋系數；αc為冷卻劑溫度反饋系數；Tf(t)和Tc(t)分別為t時刻的燃料、冷卻劑平均溫度。

1.2 燃料棒熱傳導模型

燃料棒熱傳導模型(FTM)包括穩態和瞬態熱傳導模型，用于求解燃料棒芯塊、氣隙、包殼的溫度，將包殼表面熱流密度傳遞給CFD程序，并反饋PKM燃料芯塊和包殼的溫度。

1) 穩態熱傳導模型

芯塊熱傳導：

(4)

氣隙熱傳導：

(5)

包殼內部熱傳導：

(6)

包殼與冷卻劑間的對流換熱：

Tcu-Tw=ql/(2π·Rcu·h)

(7)

換熱系數h的計算關系式：

(8)

積分熱導率[25]：

(9)

其中：Ru為芯塊半徑；Rci為包殼內半徑；Rcu為包殼外半徑；T0和Tu分別為芯塊中心、外表面溫度；Tci和Tcu分別為包殼內、外表面溫度；Tw為冷卻劑溫度；ku、kg、kz、kc依次為芯塊、間隙、包殼、冷卻劑的熱導率；h為包殼與冷卻劑間的換熱系數；ql為線功率密度；Qv為體積釋熱率；Nu為努塞爾數；De為當量直徑；k(T)為對應材料的熱導率公式；k為材料的平均熱導率；T為材料溫度；Ta、Tb為相鄰兩節點溫度。采用積分熱導率法求解燃料棒穩態溫度，根據邊界條件由外向內迭代求解燃料包殼、芯塊的溫度分布。

2) 瞬態熱傳導模型

根據能量守恒，燃料棒瞬態熱傳導模型的控制方程[26]為：

(10)

忽略燃料棒軸向傳熱，則：

(11)

(12)

其中：下標m代表材料，包括芯塊、氣隙、包殼材料；ρm為材料密度；cp,m為材料比定壓熱容；r為半徑；k為材料熱導率；i為徑向所在的節點數。采用有限差分方法求解燃料棒瞬態導熱模型，沿徑向將芯塊和包殼劃分多個節點，迭代求解燃料棒瞬態溫度分布。

1.3 計算流體力學模型

CFD模型用于計算冷卻劑溫度、速度等參數，將冷卻劑溫度、速度傳遞給FTM，并反饋PKM冷卻劑的溫度。CFD模型中的控制方程包括質量守恒方程、動量守恒方程和能量守恒方程，表達式如下。

質量守恒方程：

(13)

動量守恒方程：

(14)

能量守恒方程：

(15)

實際反應堆的幾何結構非常復雜，為降低建模難度、節省計算時間，采用多孔介質模型分析堆芯和熱交換器模型，將冷卻劑流經這些區域的壓降等效為動量守恒方程中的附加動量源項[27]，表達式為：

(16)

其中：u為速度矢量；t為時間；ρ為液體密度；ui為速度場分量；μ為湍動能黏度；p為液體壓力；xi為方向坐標；Si為外加動力源項；k為流體的熱導率；cp為比定壓熱容；ST為體積熱源項；vi為局部流速；μ(T)為隨溫度變化的流體動力黏度；ρ(T)為隨溫度變化的流體密度；Di為黏性壓降系數；Ci為慣性壓降系數。采用有限體積法(FVM)建立離散控制方程組，求解離散方程時選取SIMPLE算法計算流體的溫度與速度[28]。

1.4 不確定性分析程序SIMLAB

SIMLAB是由歐洲聯合研究中心(JRC)研發的一款不確定性開源程序，具有可視化操作界面，包括統計預處理模塊、模型執行模塊、統計后處理模塊，主要用于不確定性分析和全局敏感性分析[29]。SIMLAB程序在管理經濟可行性評估以及模糊推理系統等領域已經開展了不確定性和敏感性分析[30-31]。

SIMLAB程序的統計預處理模塊提供參數分布類型及抽樣方法，用戶根據模型要求選擇不同的參數分布及抽樣方法。SIMLAB程序中參數分布類型包括正態分布、均勻分布、韋伯分布、指數分布、γ分布、β分布等，也可以自行創建函數關系，且能可視化參數分布。抽樣方法主要包括FAST法、Sobol法、拉丁超立方方法、隨機抽樣法等。

使用非參數統計法(Wilks公式)確定抽樣次數，Wilks公式根據分布區間的置信度b及輸出參數的概率a確定所需最小抽樣數目n[32]。Wilks公式為：

1-an≥b

(17)

(1-an)-n(1-a)an-1≥b

(18)

其中式(17)適用于單側限值分布，式(18)適用于雙側限值分布。

SIMLAB程序的統計后處理模塊用于不確定性分析及敏感性評估。不確定性分析能可視化參數的概率分布、累積分布、逆累積分布，給出參數的均值、方差、標準差、上下限值等統計結果；敏感性評估給出輸入參數與輸出參數的相關系數或敏感性指數等。

1.5 物理熱工耦合不確定性分析平臺建立

熱工計算包括計算流體力學模型和燃料棒熱傳導模型，物理計算為點堆中子動力學模型，3個模型的耦合方法如圖1所示。PKM進行堆芯中子學計算，對FTM提供堆芯功率；FTM進行瞬態熱傳導計算，對CFD提供熱流密度，對PKM提供燃料溫度；CFD進行流動換熱計算，對FTM提供冷卻劑的溫度和速度。

圖1 物理熱工耦合方法

物理熱工耦合程序計算流程如圖2所示，首先讀取初始、邊界條件，調用CFD模塊計算穩態的冷卻劑溫度、速度，將CFD模塊計算結果傳遞給FTM計算燃料棒溫度分布，當前后兩次計算的各節點溫度分布小于0.1 K時，輸出穩態計算值；讀取穩態計算結果作為瞬態初始參數，將初始化參數傳遞給PKM計算堆芯功率；FTM通過堆芯功率計算燃料棒溫度分布和表面熱流密度，CFD模塊根據表面熱流密度計算冷卻劑溫度、速度分布；將燃料棒、冷卻劑溫度分布反饋給PKM用于計算反應性反饋，并將冷卻劑溫度、速度反饋給FTM。當計算時間達到耦合時間步長，物理與熱工模塊進行數據傳遞，之后進行下一時間步的計算，直至FTM計算出的燃料、包殼溫度幾乎不變，則認為瞬態事故區域穩定，停止計算。

圖2 物理熱工耦合程序計算流程

通過SIMLAB程序中模型執行模塊的外部接口，實現物理熱工耦合程序與SIMLAB程序的耦合計算，開發了物理熱工耦合不確定性分析平臺CFD/PFS。耦合方法如圖3所示：通過SIMLAB程序對耦合程序中的不確定性輸入參數進行隨機抽樣，并生成樣本文件；將樣本文件傳遞給物理熱工耦合程序，經計算獲得輸出文件；通過執行模塊外部接口，將輸出文件傳遞回SIMLAB程序，并應用SIMLAB程序統計分析樣本參數與輸出參數的相關系數，對參數進行不確定性量化分析。

基于CFD/PFS平臺的不確定性計算流程如圖4所示，計算流程如下：1) 識別重要的不確定性輸入參數；2) 在SIMLAB程序中定義重要輸入參數的取值范圍及概率分布，選取抽樣方法，獲得不同的樣本組合；3) 將每一組樣本組合傳遞給物理熱工耦合程序進行計算獲得響應值；4) 將響應值傳遞給SIMLAB程序，利用SIMLAB程序對結果進行不確定性量化，計算輸入參數與響應值的相關系數。

圖3 SIMLAB與物理熱工耦合程序的耦合方法

圖4 不確定性方法的計算流程

2 物理耦合程序驗證

2.1 分析模型

小型自然循環鉛基快堆(SNCLFR-10)是由中國科學院設計的功率為10 MW的鉛鉍合金冷卻的池式反應堆，燃料采用富集度為19.75%的UO2，燃料包殼采用316L不銹鋼，一回路冷卻劑采用液態鉛鉍合金共金體[33]。SNCLFR-10的結構示意圖示于圖5，堆芯主要設計參數參見文獻[34]。采用國際著名快堆多物理耦合分析程序SIMMER-Ⅲ，基于SNCLFR-10無保護超功率(UTOP)事故對CFD/PFS平臺物理熱工耦合程序開展Code-to-Code驗證。

圖5 SNCLFR-10的結構

2.2 耦合程序驗證

瞬態超功率事故是反應堆設計必須考慮的典型事故瞬態工況，是指向堆內突然引入一個外部反應性，導致功率急劇上升的事故。本文將反應堆初始運行條件設為穩態，在2 s內線性引入1 $反應性，整個瞬態過程無停堆動作。

UTOP事故中反應性的計算值和功率變化如圖6所示。反應性在2 s內激增，但由于燃料和冷卻劑的負反饋作用，反應性只能提高到0.88 $左右，之后由于燃料和冷卻劑的負反饋增大，反應性逐漸降低，反應性最后穩定在臨界狀態。反應堆功率在2 s內升高30倍左右，由于負反饋作用，功率逐漸降低，最后穩定在58 MW左右。SIMMER-Ⅲ的功率最后穩定在60 MW左右，耦合程序與SIMMER-Ⅲ的相對誤差為3.33%。

圖7示出燃料最高溫度和包殼最高溫度隨時間的變化。隨著功率的增大，燃料最高溫度和包殼最高溫度逐漸升高，之后隨著負反應性的增大，堆芯功率降低，燃料和包殼的最高溫度會隨之降低。相對于功率的變化，燃料和包殼的最高溫度變化幅度較低，其變化在時間上也比功率延遲。耦合程序和SIMMER-Ⅲ計算的燃料最高溫度最后分別穩定在2 500 K和2 450 K左右，兩者的相對誤差為2.04%；包殼最高溫度分別穩定在1 098 K和1 080 K左右，兩者的相對誤差為1.7%。

圖6 反應性和堆芯功率隨時間的變化

圖7 燃料最高溫度和包殼最高溫度隨時間的變化

對比物理熱工耦合程序與SIMMER-Ⅲ的計算值，兩者的誤差較小、符合度較高，說明CFD/PFS平臺計算結果具有良好的準確性和可靠性，可用于開展反應堆物理熱工耦合不確定性量化分析研究。

3 UTOP事故不確定性分析

3.1 不確定性參數的選取

反應堆中，熱工水力反饋參數會通過影響相關核素的核密度或微觀截面，從而直接影響組件等效均勻化宏觀截面參數，因此將燃料多普勒系數、堆芯軸向膨脹系數、包殼徑向膨脹系數和冷卻劑溫度系數作為不確定性輸入參數[35]。為減少抽樣取值的偶然性，假設所有輸入參數為均勻分布；將輸入參數的相對不確定度假設為5%。輸入參數的不確定性分布列于表1。

表1 輸入參數的不確定性分布

在SNCFLR-10的UTOP事故分析中，將表征反應堆運行狀態及安全特性的參數，包括總反應性峰值、功率峰值、燃料及包殼最高溫度，選取為不確定性輸出參數。

3.2 不確定性量化

選取拉丁超立方方法對輸入參數抽樣，根據式(18)，滿足兩個95%雙側限值分布的最小樣本量為93。抽取100組工況作為輸入樣本，為直觀檢測輸入參數的樣本量，將輸入參數進行歸一化處理[36]。假設第i個參數的上限值為Uij，下限值為Lij，第j個樣本的參數值為Kij，則歸一化結果Xij為：

(19)

Xij的取值在0～1之間，靠近0的樣本接近下限值，靠近1的樣本接近上限值。

圖8示出拉丁超立方抽樣樣本歸一化分布。由圖8可見，歸一化后的參數服從均勻分布，因此輸入參數也服從均勻分布，驗證了輸入參數抽樣的合理性。

圖8 拉丁超立方抽樣樣本歸一化分布

圖9示出總反應性峰值、功率峰值、包殼最高溫度、燃料最高溫度的分布。由圖9可知：總反應性峰值的上限值為0.889 13 $，下限值為0.883 89 $，名義值為0.886 63 $，名義值與限值的最大相對偏差為3.09%；功率峰值的上限值為322.2 W，下限值為298.04 W，名義值為310.29 MW，名義值與限值的最大相對偏差為3.948%；燃料最高溫度的上限值為2 746.966 K，下限值為2 560.997 K，名義值為2 656.86 K，名義值與限值的最大相對偏差為3.608%；包殼最高溫度的上限值為1 221.189 K，下限值為1 169.641 K，名義值為1 196.409 K，名義值與限值的最大相對偏差為2.237%。各目標參數的名義值處于合理的結果不確定性范圍之內，表明選取的輸入參數分布是合理的；名義值與限值的相對偏差均不超過3.95%，表明計算的目標參數精度較高。

3.3 敏感性分析

輸入參數之間的相關性可能影響輸入和輸出的相關系數，利用偏相關系數(PCC)對計算結果進行分析。偏相關系數表示消除其他影響因素后輸入參數對響應值的影響程度，取值在-1～1之間。PCC為0時，輸入變量與響應值不相關；PCC為正呈正相關，為負呈負相關；PCC絕對值靠近1時，相關性最強。假設給定隨機變量X1和X2作為輸入，輸出變量Y，在消除X2所產生的影響后，PCC計算表達式為：

圖9 總反應性峰值、功率峰值、燃料最高溫度、包殼最高溫度的分布

(20)

其中：RX1Y|X2為偏相關系數；RX1Y為X1與Y的簡單相關系數；RX1X2為X1與X2的簡單相關系數變量間的相關系數；RX2Y為X2與Y的簡單相關系數。

圖10示出輸入參數與響應值的偏相關系數。由圖10可見，輸入參數與響應值均為正相關，其中，冷卻劑溫度系數與總反應性峰值、功率峰值、燃料最高溫度及包殼最高溫度的偏相關系數均小于0.16，屬于極弱相關；包殼徑向膨脹系數、堆芯軸向膨脹系數與各響應值的偏相關系數在0.4～0.8之間，相關性較強；燃料多普勒系數與各響應值的偏相關系數大于0.94，屬于極強相關。

圖10 輸入參數與響應值的偏相關系數

4 結論

針對物理熱工耦合計算進行不確定性分析，能得到更精確的關鍵參數分布信息，可進一步提高反應堆的安全性。本文介紹了基于CFD程序FLUENT開發的不確定性分析平臺CFD/PFS，并采用國際著名快堆多物理耦合分析程序SIMMER-Ⅲ對CFD/PFS平臺開展對比驗證，進而基于CFD/PFS平臺對SNCLFR-10的UTOP事故開展了不確定性分析，并研究了不確定性參數的敏感性，所得結論如下。

1) CFD/PFS平臺和SIMMER-Ⅲ程序的事故分析結果吻合良好，所開發的CFD/PFS平臺計算結果具有較好的準確性和可靠性，可用于物理熱工耦合不確定性量化分析。

2) 通過不確定性分析計算，得到了SNCLFR-10的UTOP事故中總反應性峰值、功率峰值、燃料最高溫度及包殼最高溫度等瞬態安全參數的95/95雙側容忍限值，且瞬態安全參數的名義值均處于雙側容忍限值內，表明對UTOP事故的不確定性分析具有合理性。量化分析中名義值與限值的相對偏差小于3.95%，表明獲得的瞬態安全參數不確定性較小，計算結果精度較高。

3) 對于SNCLFR-10的UTOP事故敏感性分析，通過對4組反應性反饋系數的偏相關系數比較得知，堆芯軸向膨脹系數、包殼徑向膨脹系數及燃料多普勒系數對瞬態安全參數的影響較為顯著。其中燃料多普勒系數對結果的影響極為強烈，對反應堆安全影響最大。