數值反應堆原型系統開發及示范應用研究進展

楊 文，胡長軍，劉天才，吳明宇，王昭順，姜 強

(1.中國原子能科學研究院，北京 102413；2.北京科技大學，北京 100083)

更安全高效的核反應堆對快速精準的堆芯多物理耦合理論計算、燃料和材料服役老化行為及剩余壽命預測評估等提出了迫切需求。基于先進耦合建模和大規模并行計算技術的數值反應堆(簡稱數值堆)已成為國內外核反應堆工程與技術領域前沿熱點。數值堆不僅有望使上述需求成為可能，更可能為先進反應堆的設計優化、不同工況運行模擬優化、嚴重事故序列演示預測及燃料和材料研發提供經濟、高效的數值試驗驗證平臺。美歐等核能先進國家高度重視數值堆技術研發與應用推廣[1]。

近些年，在CASL[2-3]、NEAMS[4-5]、NURESIM[6-7]等系列大型項目持續支持下，美歐在數值堆技術研發方面取得了明顯的研發進展及初步工業應用。美國于2010年部署立項了CASL、NEAMS等大型多部門聯合研發項目。CASL項目致力于發展先進模擬與仿真技術，以改善在役輕水堆電廠的性能和效率，包括延長壽命、提升功率和增加燃耗，著重于壓力容器內部的物理現象模擬。最終，CASL推出了一套反應堆分析程序包VERA，并于2020年取得商業應用許可證[8]。NEAMS項目則主要面向四代快堆、小型模塊堆和其他新型反應堆，開發反應堆和核燃料循環系統先進數值模擬技術，重點開發了燃料產品線FPL、反應堆產品線RPL、集成產品線IPL 3個產品系列[9]。近期，美國又立項支持了第3個數值堆國家級研發項目CESAR[10]，其核心目標是在CASL和NEAMS基礎上，升級開發可在正在開發的E級超算上部署運行的確定性中子輸運(UNIC)、Monte Carlo中子輸運(OpenMC)、熱工CFD(Nek5000) 3個軟件系統，最終形成耦合的下一代反應堆仿真工具包(TRIDENT)，能在E級平臺高效執行，用于快堆設計、許可和安全分析等。歐洲數值堆的研發主要體現在NUR系列項目上。NUR旨在建立一個供歐洲核反應堆仿真通用的參考平臺，提供高精度的物理模型、軟件工具和評價結果，包括NURESIM(2005—2008)、NURISP(2009—2012)[11]、NURESAFE(2013—2016)[7]及NURE-NEXT等階段，是一整套適用于現有輕水堆及未來堆型反應堆設計、安全分析及運行應用程序的開發和驗證。

我國三大核電集團(中國核工業集團有限公司、中國廣核集團有限公司、國家電力投資集團公司)及部分高校等緊跟國際前沿，近幾年積極開展數值堆技術研發，并取得一定進展。在科技部“十三五”國家重點研發計劃“數值反應堆原型系統開發及示范應用”重點專項支持下，中國原子能科學研究院聯合北京科技大學、中國科學院計算機網絡信息中心等積極開展了數值堆關鍵技術研發。目前初步完成了基于大規模并行的數值反應堆原型系統CVR1.0主要模塊的開發，正在開展示范應用推廣工作。本文介紹數值反應堆原型系統、核心軟件開發與驗證、反應堆物理-熱工-結構多物理耦合計算，以及數值反應堆原型系統初步示范應用等研究進展。

1 數值反應堆原型系統CVR1.0

CVR1.0是核反應堆堆芯典型穩態運行工況數值模擬計算系統，可實現穩態運行工況下反應堆中子物理、熱工水力、結構力學以及反應堆燃料和材料等典型物理過程和失效行為的多物理耦合高精細模擬計算和分析預測。其核心是一套基于E級超算的涉及反應堆中子物理、熱工水力、結構力學、燃料和材料5個專業的高性能模擬計算軟件，采用多物理耦合與軟件測試架構(表1)以及“反應堆物理、熱工、力學多尺度精細化耦合計算技術”“材料輻照脆化和輻照腫脹多尺度模擬計算技術”“基于E級計算機的超大規模并行求解技術”“多源數據、多模型、多物理裝置結合的模擬驗證技術”4項關鍵技術。

表1 數值核反應堆原型系統CVR1.0Table 1 Virtual reactor system CVR1.0

2 CVR1.0核心專業軟件開發與驗證

2.1 中子物理計算軟件研發

1) 中子輸運-特征線計算軟件ANT-MOC

ANT-MOC是基于三維特征線法的堆芯中子輸運計算軟件，軟件采用直接三維特征線方法，攻克了“特征線法多堆型高精細建模與預處理技術”“面向E級架構的并行特征線法優化技術”等關鍵技術，具備穩態中子輸運方程求解能力[12]。

ANT-MOC通過了C5G7、Takeda等多個國際基準題驗證。表2和圖1分別為ANT-MOC求解C5G7基準例題[13]的有效增殖因數keff和中子通量密度分布計算結果，3種配置下keff與Monte Carlo程序計算值的相對偏差分別為0.140%、0.175%、0.270%，符合較好。目前，ANT-MOC已實現了中國實驗快堆(CEFR)活性區127組件模擬，正積極推廣其他應用。

圖1 C5G7基準測試中子通量密度分布Fig.1 Neutron flux density distribution of C5G7 benchmark test

表2 C5G7算例測試結果Table 2 Test result of C5G7

2) 中子輸運-Monte Carlo軟件ANT-RMC

ANT-RMC是基于Monte Carlo方法的堆芯中子輸運計算軟件，由項目團隊與清華大學合作開發完成。ANT-RMC能處理復雜幾何結構、采用連續能量點截面對復雜能譜和材料進行描述，并能根據實際問題的需要對臨界問題本征值和本征函數計算、精細核素鏈燃耗模擬、中子動態參數計算、在線核截面處理、核熱耦合等進行計算。支持燃耗區超百萬、計數超千億的大規模并行算法。

目前，ANT-RMC通過了多個國際基準題(VERA#3基準題、VERA#5基準題、VERA#8基準題、VENUS-Ⅱ基準題和深穿透基準題等)的驗證。表3為ANT-RMC求解VERA#5全堆基準題算例2[14]的keff計算結果，和MCNP的計算結果符合較好。

表3 VERA#5測試結果Table 3 Test result of VERA#5

2.2 熱工水力模擬軟件研發

1) 子通道分析軟件CVR-PASA

CVR-PASA是基于數據庫進行數據管理的高性能并行子通道分析軟件，采用MPI+OpenMP混合編程模型實現，使用區域分解方式進行并行任務劃分，在劃分時采用“適應于多幾何堆芯類型的自適應全堆芯子通道任務劃分與映射”關鍵技術劃分全堆芯子通道任意求解域個數。CVR-PASA可實現全堆芯、全通道熱工水力模擬，可獲得穩態流場、壓力場、溫度場以及空泡份額分布和臨界熱流密度。

CVR-PASA通過了OECD/NRC標準題等算例、壓水堆實堆數據的驗證，在天河2號超算上進行壓水堆實堆全堆芯157組件(控制體約600萬)的穩態模擬，在8 792核時模擬時間為354.66 s，同類型算例CTF需約23 min。表4和圖2是CVR-PASA模擬商用壓水堆實堆全堆芯的冷卻劑壓力等和燃料棒包殼外表面溫度分布計算結果，與實堆數據基本一致，能反映各參數在堆內的真實分布。

表4 CVR-PASA計算結果與實測數據對比Table 4 Comparison of CVR-PASA results with measured data

圖2 燃料棒包殼外表面溫度分布Fig.2 Temperature distribution on the outer surface of fuel rod cladding

2) 計算流體力學模擬軟件CVR-PACA

CVR-PACA是數值堆單相精細化計算流體力學熱工水力模擬軟件。軟件采用精確大渦模擬模型的高精度譜元方法(最高可達24階精度，常用商用軟件一般不超過5階精度)，支持異構架構，在“神威·太湖之光”超算上采用主從核異構并行方案。CVR-PACA可對反應堆的堆芯進行精確的數學物理建模，支持模擬的網格數量超百億，選用合適的湍流模型來閉合平均質量、動量和能量守恒方程，精確求解反應堆容器內部的壓力場、流場、溫度場[15-16]。CVR-PACA通過了Matis-H等多個基準題驗證。從Matis-H基準題[17]撕裂式攪混翼下游不同位置處各時均速度分量與實驗值的對比[15]可看出，CVR-PACA的LES模型及URANS SSTk-ω模型均與實驗值較接近，能較好預測由于撕裂式攪混翼引起的湍流流動狀態。

2.3 結構力學模擬軟件研發

HARSA是用于反應堆堆芯結構力學行為數值模擬的軟件，實現大規模撕裂有限元求解方法及負載均衡技術、性能分析方法，支持靜力學、流致振動、磨損等計算，具備全堆芯組件靜力變形、流致振動分析能力[18]。

HARSA通過了IAEA基準算例1、單根組件受限熱彎曲試驗算例、IAEA基準算例6等算例驗證，在天河2號上，模擬網格量達137億。表5和圖3是使用HARSA按照IAEA基準算例1的要求建立單盒組件自由變形模型，開展有限元靜力計算的計算結果。與IAEA基準算例1給出的理論解對比，HARSA計算結果的偏差在3.8%以內，計算準確性較好。

表5 組件軸向伸長量計算結果對比Table 5 Comparison of fuel assembly axial elongation deformation displacement

圖3 組件軸向伸長形變位移模擬Fig.3 Simulation of fuel assembly axial elongation deformation displacement

2.4 燃料性能分析軟件研發

ATHENA是用于對壓水堆燃料元件進行單棒和多棒性能分析的軟件，軟件開發了燃料溫度分布、變形計算、裂變氣體釋放及內壓等模型，結合燃料元件熱工-力學多物理耦合計算分析耦合方案，基于先進并行計算方法，建立了燃料并行分析能力。支持計算燃料元件溫度分布、變形、裂變氣體釋放、燃料棒內壓、包殼腐蝕等參數，具備高效、精細化全堆芯燃料性能分析的能力[19]。

ATHENA通過國際標準基準例題、同類程序計算結果、典型商用壓水堆核電站數據進行了驗證，獲得了燃料溫度、燃料變形、燃料內壓等參數計算對比結果，分析表明ATHENA軟件計算結果合理可靠。圖4為ATHENA和對標軟件METEOR基于典型商用壓水堆核電站燃料數據計算得到的燃料芯塊中心溫度隨時間變化的趨勢圖，計算結果基本符合。

圖4 峰值溫度隨時間變化Fig.4 Peak temperature change with time

2.5 材料多尺度模擬軟件研發

MISA是項目團隊自主研發的系列反應堆用結構材料輻照損傷多尺度模擬軟件集，包括微觀尺度的分子動力學模擬軟件MISA-MD、空位躍遷機制的動力學Monte Carlo模擬軟件MISA-KMC、空位和間隙雙躍遷機制的動力學Monte Carlo模擬軟件MISA-AKMC[20]、介觀尺度的相場模擬軟件MISA-PF、隨機團簇動力學模擬軟件MISA-SCD。軟件集可對核材料輻照損傷進行多尺度模擬，獲得用于RPV鋼輻照脆化、堆內構件輻照腫脹預測所需的微觀信息。

MISA系列軟件主要通過國外軟件模擬算例及結果進行驗證，模擬結果與對比軟件較吻合。圖5分別為對LAMMPS和MISA-MD模擬結果作統計分析得到的Frenkel缺陷對數量Nfp隨時間變化趨勢圖，結果較吻合。

圖5 級聯碰撞過程中Frenkel缺陷對數量隨時間變化圖Fig.5 Evolution of number of displacement cascades’ Frenkel defect with time

3 多物理耦合系統與數值模擬軟件測試框架研發

3.1 多物理耦合系統研發

數值堆多物理耦合模擬平臺SPIDER是一款能實現核反應堆高精細耦合模擬的集成軟件平臺，其整體架構如圖6所示，支持各物理模塊之間的不斷交互和統一協作，并針對各模塊設計建立數據庫，用于促進核反應多物理耦合過程中的數據交互和模擬數據的統一管理，具備數據支撐能力。

圖6 SPIDER系統架構Fig.6 System structure of SPIDER

SPIDER平臺中包括核熱耦合和流固耦合兩大核心模塊。核熱耦合過程十分復雜，存在著強烈的相互反饋，一方面，通過物理計算可得出裂變反應率和中子通量分布，從而計算出功率分布，并進一步求得反應堆的傳熱特性和熱工水力特性；另一方面，熱工水力的計算會對物理過程進行反饋，影響中子通量分布。SPIDER平臺中核熱耦合模塊針對自主研發的并行子通道模擬軟件CVR-PASA、特征線法中子輸運模擬軟件ANT-MOC設計實現了基于文件的數據交互耦合流程，為反應堆物理-熱工耦合的數值模擬提供了更為準確的堆內參數變化。流固耦合模塊針對自主研發的計算流體力學軟件CVR-PACA和結構力學軟件HARSA同樣設計并實現了基于文件的數據交互耦合流程[21]。

SPIDER通過了C5G7、中國實驗快堆燃料組件、Hoogenboom-Martin等多個算例驗證[22]。圖7為采用OECD的C5G7-Rodded B例題[13]進行核熱耦合功能驗證的計算結果，燃料棒表面溫度分布合理，與文獻符合較好。采用中國實驗快堆燃料組件算例進行流固耦合功能驗證的計算，燃料棒束最大位移出現在中心燃料棒束的頂端，為27.83 mm[21]。

圖7 燃料棒表面溫度分布Fig.7 Temperature distribution of fuel rod surface

3.2 數值模擬軟件測試框架研發

數值模擬軟件測試框架(NuSTA)是為數值模擬程序提供一種可用的測試框架，可生成算例，進行單元測試、蛻變測試、屬性測試、精度測試等。框架的主要功能為：統一測試用例定義規范，支持多種不同輸入輸出格式的高性能數值模擬軟件；支持Verification & Validation流程中的主要活動，實現傳統測試、蛻變測試、精度階分析、回歸測試、誤差分析等功能接口的集成和封裝；完成蛻變測試、誤差分析工具(ScDebug)的開發；實現測試數據的管理和展示。框架的技術特點在于：將蛻變測試和數值模擬程序相結合，有效緩解測試用例少且難構造的問題；研發ScDebug精度測試工具，相比目前動態精度分析工具性能提高70%。NuSTA框架的總體結構設計如圖8所示。

圖8 NuSTA框架總體結構Fig.8 Overall structure of NuSTA framework

4 數值反應堆原型系統示范應用

完成CVR1.0的開發和驗證后，目前正開展相關示范應用工作，重點在關鍵材料失效分析和反應堆工程設計或設計驗證。

RPV鋼輻照脆化多尺度模擬研究，主要利用自主開發軟件針對RPV鋼及其模型合金的輻照脆化機理開展研究。項目團隊利用MISA-KMC軟件研究了合金元素Mn對富Cu團簇形核的影響，分析了Mn含量變化對團簇析出的影響[23]；利用位錯動力學方法結合分子動力學和分子靜力學計算獲得的缺陷釘扎力，研究了FeCu模型合金中Cu析出物導致硬化的機理，該工作將微觀結構演化模擬與宏觀性能預測建立橋梁，對深入研究RPV鋼輻照硬化機理以及預測輻照脆化趨勢具有重要意義[24]。

基于MISA-MD與MISA-SCD耦合模擬了Fe0.3at.%Cu合金中微觀結構演化及其導致的輻照硬化、脆化行為，圖9為Fe0.3at.%Cu合金經中子輻照后微觀結構的數密度和平均半徑隨損傷劑量(dpa)的變化，輻照后產生的微觀結構主要為純Cu團簇、Cu_Vac復合團簇、孔洞和位錯環。位錯環的直徑小于1 nm，接近0.75 nm，可認為位錯環對位錯運動阻礙作用較弱；孔洞的直徑約為2 nm，但孔洞在演化后期才出現，且數密度很低，可忽略其對位錯運動阻礙作用；Cu_Vac類型團簇尺寸較大，但輻照前期Cu_Vac數密度較低，可不考慮其對位錯運動阻礙作用。當僅考慮純Cu團簇對輻照硬化、脆化的影響時，基于Orowan硬化模型及RPV鋼輻照后硬化與脆化關系，可獲得純Cu團簇引起材料韌脆轉變溫度增量隨損傷劑量變化關系。圖10為純Cu團簇引起材料韌脆轉變溫度增量(DBTT)與實驗結果對比，由圖可知，模擬結果與實驗結果[25]較吻合，驗證了結果的可靠性。

圖9 中子輻照下Fe0.3at.%Cu微觀結構半徑和數密度隨損傷劑量演化關系Fig.9 Relationship of radius microstructure and number density of Fe0.3at.%Cu and damage dose under neutron irradiation

圖10 韌脆轉變溫度增量隨損傷劑量變化關系Fig.10 Predicted increasement of DBTT as a function of damage dose

5 總結與展望

在國家重點研發計劃“數值反應堆原型系統開發及示范應用”重點專項支持下，初步完成了基于大規模并行的數值反應堆原型系統CVR1.0主要模塊的開發，并在反應堆壓力容器輻照脆化機理認知與輻照脆化趨勢預測、快堆設計分析等方面實現了初步示范應用。

CVR1.0核心是基于E級超算的反應堆中子物理、熱工水力、結構力學、燃料和材料5個專業高性能模擬計算軟件，可實現反應堆堆芯穩態運行工況下中子物理、熱工水力、結構力學以及反應堆燃料和材料等典型物理過程和失效行為的多物理耦合高精細模擬計算和分析預測。

數值堆后續研發重點是持續推進CVR1.0應用反饋和升級優化的同時，進一步拓展CVR1.0至全堆系統、全工況運行模擬，并逐步實現反應堆設計研發與運行工程應用。