楊 文，胡長軍，劉天才，吳明宇，王昭順，姜 強(qiáng)

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

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

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

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

1 數(shù)值反應(yīng)堆原型系統(tǒng)CVR1.0

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

表1 數(shù)值核反應(yīng)堆原型系統(tǒng)CVR1.0Table 1 Virtual reactor system CVR1.0

2 CVR1.0核心專業(yè)軟件開發(fā)與驗證

2.1 中子物理計算軟件研發(fā)

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

ANT-MOC是基于三維特征線法的堆芯中子輸運(yùn)計算軟件，軟件采用直接三維特征線方法，攻克了“特征線法多堆型高精細(xì)建模與預(yù)處理技術(shù)”“面向E級架構(gòu)的并行特征線法優(yōu)化技術(shù)”等關(guān)鍵技術(shù)，具備穩(wěn)態(tài)中子輸運(yùn)方程求解能力[12]。

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

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

表2 C5G7算例測試結(jié)果Table 2 Test result of C5G7

2) 中子輸運(yùn)-Monte Carlo軟件ANT-RMC

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

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

表3 VERA#5測試結(jié)果Table 3 Test result of VERA#5

2.2 熱工水力模擬軟件研發(fā)

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

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

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

表4 CVR-PASA計算結(jié)果與實測數(shù)據(jù)對比Table 4 Comparison of CVR-PASA results with measured data

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

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

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

2.3 結(jié)構(gòu)力學(xué)模擬軟件研發(fā)

HARSA是用于反應(yīng)堆堆芯結(jié)構(gòu)力學(xué)行為數(shù)值模擬的軟件，實現(xiàn)大規(guī)模撕裂有限元求解方法及負(fù)載均衡技術(shù)、性能分析方法，支持靜力學(xué)、流致振動、磨損等計算，具備全堆芯組件靜力變形、流致振動分析能力[18]。

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

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

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

2.4 燃料性能分析軟件研發(fā)

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

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

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

2.5 材料多尺度模擬軟件研發(fā)

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

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

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

3 多物理耦合系統(tǒng)與數(shù)值模擬軟件測試框架研發(fā)

3.1 多物理耦合系統(tǒng)研發(fā)

數(shù)值堆多物理耦合模擬平臺SPIDER是一款能實現(xiàn)核反應(yīng)堆高精細(xì)耦合模擬的集成軟件平臺，其整體架構(gòu)如圖6所示，支持各物理模塊之間的不斷交互和統(tǒng)一協(xié)作，并針對各模塊設(shè)計建立數(shù)據(jù)庫，用于促進(jìn)核反應(yīng)多物理耦合過程中的數(shù)據(jù)交互和模擬數(shù)據(jù)的統(tǒng)一管理，具備數(shù)據(jù)支撐能力。

圖6 SPIDER系統(tǒng)架構(gòu)Fig.6 System structure of SPIDER

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

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

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

3.2 數(shù)值模擬軟件測試框架研發(fā)

數(shù)值模擬軟件測試框架(NuSTA)是為數(shù)值模擬程序提供一種可用的測試框架，可生成算例，進(jìn)行單元測試、蛻變測試、屬性測試、精度測試等?？蚣艿闹饕δ転椋航y(tǒng)一測試用例定義規(guī)范，支持多種不同輸入輸出格式的高性能數(shù)值模擬軟件；支持Verification & Validation流程中的主要活動，實現(xiàn)傳統(tǒng)測試、蛻變測試、精度階分析、回歸測試、誤差分析等功能接口的集成和封裝；完成蛻變測試、誤差分析工具(ScDebug)的開發(fā)；實現(xiàn)測試數(shù)據(jù)的管理和展示?？蚣艿募夹g(shù)特點在于：將蛻變測試和數(shù)值模擬程序相結(jié)合，有效緩解測試用例少且難構(gòu)造的問題；研發(fā)ScDebug精度測試工具，相比目前動態(tài)精度分析工具性能提高70%。NuSTA框架的總體結(jié)構(gòu)設(shè)計如圖8所示。

圖8 NuSTA框架總體結(jié)構(gòu)Fig.8 Overall structure of NuSTA framework

4 數(shù)值反應(yīng)堆原型系統(tǒng)示范應(yīng)用

完成CVR1.0的開發(fā)和驗證后，目前正開展相關(guān)示范應(yīng)用工作，重點在關(guān)鍵材料失效分析和反應(yīng)堆工程設(shè)計或設(shè)計驗證。

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

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

圖9 中子輻照下Fe0.3at.%Cu微觀結(jié)構(gòu)半徑和數(shù)密度隨損傷劑量演化關(guān)系Fig.9 Relationship of radius microstructure and number density of Fe0.3at.%Cu and damage dose under neutron irradiation

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

5 總結(jié)與展望

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

CVR1.0核心是基于E級超算的反應(yīng)堆中子物理、熱工水力、結(jié)構(gòu)力學(xué)、燃料和材料5個專業(yè)高性能模擬計算軟件，可實現(xiàn)反應(yīng)堆堆芯穩(wěn)態(tài)運(yùn)行工況下中子物理、熱工水力、結(jié)構(gòu)力學(xué)以及反應(yīng)堆燃料和材料等典型物理過程和失效行為的多物理耦合高精細(xì)模擬計算和分析預(yù)測。

數(shù)值堆后續(xù)研發(fā)重點是持續(xù)推進(jìn)CVR1.0應(yīng)用反饋和升級優(yōu)化的同時，進(jìn)一步拓展CVR1.0至全堆系統(tǒng)、全工況運(yùn)行模擬，并逐步實現(xiàn)反應(yīng)堆設(shè)計研發(fā)與運(yùn)行工程應(yīng)用。