反應堆核-熱-燃耗多物理耦合框架研究與應用

吳明宇，朱 迎，盧 旭，苗 雪，吳宗蕓，李 龍，胡 赟，趙民富

(1.中國原子能科學研究院 反應堆工程技術研究所，北京 102413；2.北京科技大學 計算機與通信工程學院，北京 100083)

反應堆中子學、熱學、流體力學是相互耦合的復雜系統，每個物理過程由偏微分方程或方程組描述。例如：中子物理過程由中子輸運方程描述，燃耗過程由燃耗方程組描述，瞬態過程由中子動力學方程描述，冷卻劑的換熱流動過程由傳熱方程和流體力學方程組描述。上述方程中求解的物理量往往是互相影響、相互耦合的。例如，中子輸運計算得到反應堆內的中子通量分布情況，是進行燃耗計算的必需數據，也為熱工水力計算提供了相對功率分布；由熱工水力輸運計算得到的堆內冷卻劑的密度和燃料溫度分布，可通過冷卻劑核素密度的變化及多普勒效應改變材料的宏觀截面，從而對中子的輸運過程產生影響；燃耗計算得到的新的燃料核素及核素密度，不但會通過改變燃料總截面影響中子輸運過程，也給熱工水力計算提供了總功率及各燃耗區的實際功率。實際反應堆的運行過程是中子輸運、燃耗、熱工水力學等多物理過程相互作用、相互影響的結果。

隨著計算機硬件計算能力的不斷提升，可視化先進建模、多物理、多尺度、精細化模擬獲得了快速發展。美國制定并執行了CASL計劃，構建了反應堆應用虛擬環境(VERA)，主要是充分利用現有軟件建立松耦合環境，產生一個相比傳統工具具有更高分辨率的耦合計算系統，對反應堆整個壽期運行歷史進行預測。進一步的NEAMS計劃是基于美國愛達荷國家實驗室開發的面向對象多物理模擬(MOOSE)框架[1]，開發多物理耦合集成平臺，實現系統-物理-熱工-結構-材料及燃料的高精細耦合模擬軟件開發。MOOSE框架采用了模塊化的方法，允許科學家和工程師創建新的完全耦合的多物理應用程序。法國電力公司的數值模擬軟件平臺Salome能夠用于數值模擬的通用預處理、后處理和工作流管理操作，能便捷地進行軟件間的相互耦合運算和處理。近些年，國內核能領域也逐漸開始重視數值堆相關研究開發工作，并取得了不同程度的進展。中國工程物理研究院高性能數值模擬軟件中心分析了數值反應堆多物理耦合關鍵技術[2]，提出采用中間件把許多共性的計算方法、領域并行編程的復雜問題，通過區域分解、數據分解、多級并行得到化解。

本文根據反應堆堆芯多物理耦合分析需求，研究多物理耦合算法，并完成耦合程序開發，實現中子物理、燃耗、熱工子通道的多物理耦合計算。

1 多物理耦合求解方法

堆芯多物理耦合示意圖如圖1所示。耦合系統求解方法分為兩種，如圖2所示，一種是采用等效近似來進行物理量的解耦，然后通過參數傳遞和迭代計算來得到滿意的結果；另外一種方法就是構建微分方程組直接聯立求解，所有物理量同時更新，如近幾年研究的JFNK算法[3-4]。

圖1 堆芯多物理耦合示意圖Fig.1 Schematic diagram of reactor core multi-physical coupling

圖2 多物理耦合計算方法示意圖Fig.2 Schematic diagram of multi-physical coupling calculation method

根據求解方法的不同，多物理系統設計中采用的耦合方式主要有以下3種：1) 松耦合，通常是在現有計算分析程序的基礎上開發接口程序，程序間的數據傳遞基于雙方的輸入輸出文件展開，這種耦合方式對程序的改動小，能充分利用已有程序，但程序間的相容性不好，容易存在收斂性的問題；2) 緊耦合，將耦合系統中描述不同物理過程的偏微分方程或方程組進行聯立，在同一時間步長內，同時求解耦合系統中所有的非線性微分方程并進行變量更新，在計算過程中全部或部分使用內存進行數據的傳遞，這類程序涉及大量源代碼的編寫、重構與修改，需要解決相容性和通用性，是目前研發的熱點；3) 耦合框架，搭建統一的I/O平臺進行數據的處理和傳遞，各專業程序分別作為獨立求解器內嵌于平臺中，數據可通過文件或內存傳遞，是前兩種耦合方式的集成，代表程序有美國愛達荷國家實驗室開發的MOOSE系統[5-6]，ANASYS和COMSOL等大型商用軟件也具備耦合框架平臺的特點。

2 堆芯多物理耦合系統架構設計

含有多個偏微分方程的耦合系統直接求解的難度很大，即便是精細化求解單一物理過程的偏微分方程也存在諸多困難，因此本文工作首先采用基于參數傳遞的松耦合方法，基于單一物理過程的蒙特卡羅程序、特征線程序、燃耗程序、子通道程序，采用Master-Slave主從架構、Python語言編程，構建中子物理、燃耗、熱工子通道的多物理耦合系統，如圖3所示。在耦合架構設計中采用分層與封裝策略，按照功能需求分為系統控制層、驅動層、數據層、執行層。同時基于網格類型和模擬參數，建立網格數據和參數的表示規范以及統一接口，節省存儲資源、提高訪問效率；基于不同空間離散形式的高效定位數據索引與差值算法，實現參數的高效轉換與傳遞，提高耦合效率；基于具體耦合計算模型，建立耦合參數傳遞函數的映射關系，提高耦合參數傳遞的準確性；基于單級與多級時間步長控制方法，平衡模擬精度與計算時間。針對不同模擬模型，研究收斂控制方式的選擇對程序穩定性和收斂性的影響。

圖3 堆芯多物理耦合系統架構設計Fig.3 Architecture design of core multi-physical coupling system

3 耦合系統初步測試

3.1 輸運-燃耗耦合計算測試

輸運-燃耗耦合計算中，堆芯的中子通量密度、能譜分布及同位素核密度的參數傳遞與迭代計算原理如圖4所示。蒙特卡羅程序進行臨界計算，得到中子通量、單群反應截面等數據，傳遞給點燃耗程序；點燃耗程序進行燃耗計算，得到新的核素密度，傳遞給蒙特卡羅程序。通過數據的往返傳遞，從而處理燃耗計算的全過程。本次工作中子輸運計算程序采用清華大學反應堆工程計算分析實驗室(REAL)研發的三維粒子輸運蒙特卡羅程序ANT-RMC[7-8]，該程序包含燃耗計算模塊DEPTH[9]，采用程序內耦合模式開發，能夠處理含1 500余種核素的精細燃耗鏈。程序首先通過臨界計算(連續能量)模塊得到中子通量、單群反應截面等數據，傳遞給點燃耗模塊DEPTH，通過燃耗計算得到新的核素密度，傳遞給臨界計算模塊。

圖4 輸運-燃耗耦合計算流程Fig.4 Coupled calculation process of transport-burnup

ANT-RMC程序通過負載平衡和區域分解能實現大規模輸運燃耗耦合計算，本文采用ANT-RMC程序對Hoogenboom-Martin基準題進行驗證。Hoogenboom-Martin基準題[10]是一個含有燃耗計算的全堆芯基準計算模型，堆芯共包含241個燃料組件，每個燃料組件包含17×17個柵元，其中包括264個燃料棒，25個控制棒，燃料棒的高度為366 cm，沿軸向劃分成24等份，堆芯結構與區域分解如圖5所示。

圖5 Hoogenboom-Martin基準題堆芯結構與區域分解Fig.5 Core structure and domain decomposition of Hoogenboom-Martin benchmark

對Hoogenboom-Martin基準題模型進行精細化pin-by-pin描述，采用96核并行規模，實現了百萬燃耗區全堆三維精細燃耗計算。時間步長為250 d時的徑向功率分布與有效增殖因數keff變化如圖6所示。

圖6 徑向功率分布與keff的變化Fig.6 Change of radial power distribution and keff

3.2 核熱耦合計算測試

在堆芯核熱耦合計算中堆芯熱工計算采用子通道分析方法。子通道分析方法將堆芯內復雜的流通面積劃分子通道，通過子通道間質量守恒、能量守恒、動量守恒方程求解堆芯組件內的流場和溫度場。本次計算采用的子通道計算程序為自主研發的子通道模擬程序CVR-PASA，CVR-PASA程序具備壓水堆子通道計算功能[11]，同時也開發了面向六角形組件的液態金屬快堆子通道計算功能，采用多層級索引實現從堆芯-組件-棒束-子通道多層級映射描述[12]，對鈉冷快堆進行精細建模。快堆組件子通道映射示意圖如圖7所示。

圖7 快堆組件子通道映射示意圖Fig.7 Mapping diagram of fast reactor assembly subchannel

中子輸運計算程序采用三維蒙特卡羅程序ANT-RMC，擁有中子光子輸運、全能區在線截面處理、平衡氙、燃耗計算和衰變計算、全堆倒換料、MPI與OpenMP混合并行、接續計算等功能，具備在線截面處理功能，并能夠在物理熱工燃耗耦合中綜合應用。ANT-RMC支持包括外耦合(由第三方腳本控制)、混合耦合(ANT-RMC系統調用控制熱工計算)、內耦合(兩者融為一個程序)在內的多種核熱耦合方式，適用于不同的反應堆分析場景和不同計算平臺的要求，具有分析實際反應堆多循環、多物理場(輸運、燃耗、熱工水力)耦合問題的能力。

核熱耦合系統聯立公式如圖8所示，直接求解帶熱工反饋的輸運方程是非常困難的，目前較常用的方式是Picard迭代方法：首先設定一個初始的溫度和密度場的分布，然后開展輸運計算，功率分布根據臨界計算結果不斷地予以更新，新的功率分布作為熱工水力計算的條件，通過流動、傳熱及燃料棒導熱計算，得到新的冷卻劑溫度、密度及燃料棒溫度參數，存儲在中子物理幾何模塊內的對應柵元下，迭代過程直到滿足收斂條件為止。

圖8 核熱子通道耦合計算示意圖Fig.8 Schematic diagram of nuclear thermal coupling calculation

核熱耦合計算模型為商用示范快堆堆芯與組件[13]，采用ANT-RMC開展堆芯pin-by-pin精細化計算，燃料活性區軸向等分為20段，通過輸運計算給出組件內部徑向與軸向功率分布，然后采用CVR-PASA開展子通道計算，在穩態工況下，經過迭代反饋，燃料組件出口溫度的計算結果與文獻[14]中的組件子通道出口溫度分布相近。耦合計算的溫度與冷卻劑密度分布如圖9所示。通過核熱耦合計算能夠更為精確地描述軸向與徑向功率分布對冷卻劑的影響。核熱耦合計算后續完善可用于鈉冷快堆堆芯熱工水力計算及流量優化分析。

圖9 快堆組件的溫度分布(a)和冷卻劑密度分布(b)Fig.9 Temperature distribution (a) coolant density distribution (b) of fast reactor assembly

4 結論與展望

單一的物理計算難以全面模擬堆芯物理過程，多物理耦合計算成為研究多工況下堆芯行為演化規律的重要手段。基于先進耦合建模和大規模并行計算技術的數值反應堆已成為國內外領域前沿熱點。數值反應堆不僅使上述需求成為可能，更為先進反應堆的設計優化、不同工況運行模擬優化、嚴重事故序列演示預測及燃料和材料研發提供一個經濟高效的試驗平臺。本文構建了基于中子輸運、燃耗、熱工子通道的堆芯多物理耦合系統，通過輸運-燃耗耦合計算測試和核熱耦合計算測試初步驗證了耦合系統的功能。

在后續工作中，將繼續拓展多物理耦合系統，形成較為完備的全堆芯中子學、熱工分析、燃料性能分析多物理耦合計算系統，同時開展瞬態工況下的緊耦合計算方法研究。通過燃料元件-組件-全堆芯的多尺度建模和物理-熱工-燃料性能多物理場耦合計算，給出全堆芯精細化中子場分布、流場分布、溫度場分布以及燃料組件性能在整個壽期內的精細化、量化分析結果，可以改變傳統的以局部最大功率或最熱組件為限值的保守性包絡計算方法，為建立堆芯行為演化規律提供重要支撐。