基于統一耦合框架的堆芯物理熱工耦合程序的開發及驗證

姜 榮， 馮文培， 陳紅麗， 強勝龍， 李治剛， 潘俊杰， 張喜林， 羅 曉

(1. 中國科學技術大學核科學技術學院， 合肥 230027； 2. 中國核動力研究設計院， 成都 610213)

1 引 言

反應堆系統內中子物理、熱工水力、燃料力學等多物理過程間存在復雜的耦合關系. 為了實現反應堆多物理多尺度耦合，提高模擬計算的精確度，擴大計算規模，減少重復性工作[1]，基于統一耦合框架的多物理程序耦合成為研究熱點. 這些耦合框架或耦合平臺包括ICoCo、MOOSE平臺、OpenPALM、preCICE等，其中通過ICoCo理念封裝實現多物理程序的耦合方法在國外已開展較多研究.

法國CEA基于ICoCo理念實現了中子學程序CRONOS2、兩相流分析程序FLICA4、系統分析程序CATHARE程序間的耦合[2]. 德國學者在NURESAFE項目框架下實現了子通道程序CTF和中子程序DYN3D的耦合[3]. 為了更好地模擬反應堆壓力容器中的三維熱工水力現象，德國KIT使用ICoCo接口實現了系統程序TRACE和CFD程序TrioCFD之間的耦合，該耦合系統已集成于SALOME平臺[4]. 中國科學技術大學和德國KIT合作開發了CFD程序TrioCFD和子通道程序SubChanFlow之間基于ICoCo理念和MEDCoupling庫的統一耦合框架的耦合[5]. 為了更好地模擬鈉冷快堆熱工水力現象，CEA實現了在ASTRID項目框架下系統程序CATHARE、子通道程序TrioMC和CFD程序TrioCFD之間的耦合[6]. 國內對于ICoCo封裝理念的研究還較少，因此對其開展更多的研究從而為多物理多尺度耦合提供基礎是有必要的.

本文基于ICoCo封裝理念和MED庫數據傳遞模型，對中子時空動力學程序CORCA-K和子通道分析程序CORTH進行了基于耦合框架的物理-熱工耦合程序的開發. 采用NEACRP國際輕水堆彈棒基準[7]進行了程序測試，包括對封裝后程序和耦合程序的測試，并與基準結果以及CORCA-K原有耦合接口的計算結果進行對比.

2 理論依據

ICoCo接口是在NURISP項目范圍內[8]開發的面向對象的通用代碼耦合接口，首先應用于CATHARE和TrioCFD的多尺度熱工水力耦合[9]. ICoCo接口基于C++開發，通過定義Problem通用母類來定義函數，這些函數的功能包括初始化、時間步推進、保存和還原以及場數據傳遞，表1給出了ICoCo定義的重要接口及功能.

表1 重要的ICoCo接口及功能

基于ICoCo理念對程序進行封裝，即對程序開發上述接口. 其中，在進行場數據交換時采用功能強大的網格和區域操作庫MEDCoupling實現場數據的交互和映射. MEDCoupling場數據結構包含足夠多的信息，可以進行各種插值操作，也可以在不同進程間通過并行模式和分布式模式交換數據結構信息. 進行ICoCo封裝時需要對每個程序定義一套或多套MED網格，以進行輸入場數據和輸出場數據的傳遞和映射.

為了實現基于統一耦合框架的物理-熱工程序的耦合，采用C++ 編寫的supervisor程序實現對各子程序的調用. 圖1給出了C++ supervisor進行程序調用的示意圖.

圖1 supervisor調用程序示意圖

3 程序開發

3.1 CORCA-K的ICoCo封裝

CORCA-K程序是中國核動力研究設計院自主研發的三維瞬態中子學計算軟件，其核心功能是進行三維瞬態中子擴散方程的數值求解計算，空間離散采用第二類邊界條件格林函數節塊法，時間離散采用二級二階對角線隱式龍格庫塔格式[10]. CORTH程序是中國核動力研究設計院自主研發的子通道分析程序，可以進行反應堆熱工水力設計與安全分析，也可以用來進行新型燃料組件的設計與研發[11].

基于統一耦合框架的堆芯物理程序與熱工程序的耦合，需要對CORCA-K和CORTH進行ICoCo封裝. 首先，對CORCA-K的源程序進行功能模塊劃分，如劃分為啟動模塊、初始化模塊、時間步計算模塊、中止模塊等. 然后根據功能模塊開發表1所示的ICoCo接口.

在編寫場數據交換與映射相關接口時，根據CORCA-K程序的網格特征，建立三維結構化MEDCoupling網格來進行數據交換. 同時，設置CORCA-K的輸入物理場為有效燃料溫度、慢化劑溫度、冷卻劑密度和燃料芯塊溫度，CORCA-K的輸出物理場為功率分布.

完成CORCA-K的ICoCo封裝后，編譯得到CORCA-K的動態鏈接庫，編寫supervisor調用動態鏈接庫，即可實現CORCA-K的計算功能. 按照同樣的方法對CORTH進行封裝并編譯得到CORTH的動態鏈接庫.

3.2 耦合實現

在得到封裝后CORCA-K和CORTH的動態鏈接庫后，編寫supervisor程序，編譯得到基于統一耦合框架的CORCA-K和CORTH的耦合程序，圖2為該耦合程序進行計算的流程圖.

圖2 CORCA-K與CORTH耦合計算流程圖Fig.2 The coupling calculation flowchart of CORCA-K and CORTH

圖3表示的是CORCA-K和CORTH的MED網格. CORCA-K創建的是三維結構化網格，CORTH創建的是三維非結構化網格. 在進行數據傳遞時，CORCA-K計算得到堆芯三維功率分布，通過MED網格映射傳給CORTH，CORTH將新的功率分布作為輸入，重新計算得到平均多普勒燃料溫度、慢化劑平均溫度、最高燃料溫度和冷卻劑密度，再通過MED網格映射傳給CORCA-K，CORCA-K重新計算得到功率分布，這樣循環下去，實現二者的耦合.

CORTH非結構化網格

4 數值驗證

4.1 測試算例

選擇NEACRP輕水堆彈棒基準中的C1和C2為測試算例，對封裝后的CORCA-K和基于耦合框架的耦合程序進行測試. 由OECD/NEA發布的輕水堆彈棒基準，即在熱態零功率和熱態滿功率下，堆芯中一組控制棒快速彈出造成正反應性引入. NEACRP基準中堆芯由157組燃料組件和64組反射組件組成，其中49組燃料組件中布置有控制棒，各組件的徑向尺寸是21.606 cm×21.606 cm，軸向總長度為427.3 cm，如圖4所示. 瞬態分析計算5 s，0.1 s時發生彈棒事故，0～1 s內，時間步長為0.005 s，1～5 s內，時間步長為0.05 s. 表2給出了NEACRP彈棒基準中C1和C2運行數據的描述.

4.2 封裝后CORCA-K測試

通過編寫supervisor程序檢驗每個ICoCo接口是否可以正常調用和實現其功能，從而對封裝后的CORCA-K進行測試. 以C1和C2為研究對象，分別采用不同的計算模式如臨界計算、彈棒計算等進行計算，全面驗證封裝后CORCA-K的功能. 經過多次驗證計算，對于相同的算例和相同的計算模式，封裝后CORCA-K與封裝前CORCA-K計算得到的輸出卡內容完全一致. 由此可見，CORCA-K的封裝是正確的. 由于封裝前后輸出卡完全一致，因此這里不列出具體的數據.

圖4 輕水堆基準C1和C2徑向堆芯布置圖

表2 C1和C2的堆芯運行數據

4.3 耦合程序測試

使用基于耦合框架的CORCA-K和CORTH耦合程序計算C1和C2，并使用CORCA-K已開發的與CORTH直接耦合的接口CORTHV2計算相同的算例作為對比. 由于CORTHV2耦合接口計算的結果與基準參考解已經存在偏差，基于耦合框架的CORCA-K與CORTH的耦合程序是在CORCA-K和CORTH兩程序的基礎上開發，因此為驗證耦合框架的耦合方式的準確性，同時排除CORCA-K和CORTH本身的計算誤差，主要對比基于耦合框架的CORCA-K與CORTH的耦合程序與CORTHV2的計算結果. NEACRP基準的PANTHER修正后的計算結果[12]作為參考.

表3、表4給出了使用兩種耦合方式計算C1和C2得到的穩態和瞬態計算結果，下文以“ICoCo”代表基于耦合框架的耦合程序，“CORTHV2”代表CORCA-K原來的耦合接口. 從表3中可知，C1使用兩種方式進行穩態計算得到的臨界硼濃度均為1 128.79 ppm，與基準結果僅差0.49 ppm，C2使用兩種方式計算得到的臨界硼濃度均為1 154.73 pcm，與基準結果僅差1.87 pcm. 一方面可以驗證ICoCo和CORTHV2兩種耦合方式進行穩態計算的準確性，另一方面也可以驗證封裝后CORCA-K進行臨界計算的準確性.

C1中，兩種耦合方式最大功率出現的時間僅差0.002 5 s，最大功率和5 s時的功率相差0.228 9，最大功率相差0.23. 考慮到CORCA-K與CORTH之間的數據傳遞需要經過MED網格的交互與映射，因此可能會造成一定誤差，之后可以進行網格敏感性研究，進一步確定誤差范圍. 對于5 s時的功率，兩種耦合方式計算的結果相同，均為0.16. 對于5 s時的平均多普勒燃料溫度，二者的計算結果也非常接近，最大不超過0.5 ℃，燃料最高溫度相差不超過3 ℃. 同時，兩種耦合方式的結果與基準結果也基本符合. C2中，兩種耦合方式的計算結果更加一致. 二者最大功率到達時間僅差0.002 5 s，最大功率及5 s時的功率完全一致(1.07和1.03)，5 s時的平均多普勒燃料溫度和燃料最高溫度相差均不超過1 ℃.

表3 不同耦合方式計算得到的C1彈棒基準的穩態和瞬態結果

圖5～圖10分別給出了使用兩種耦合方式計算得到的相對功率變化、平均多普勒燃料溫度、最高燃料溫度的變化趨勢. 由圖中可見，最大功率、最大功率出現的時間、5 s時的功率、平均多普勒燃料溫度和燃料最高溫度的變化趨勢二者皆吻合得很好，與PANTHER的結果略有差距，但相差很小.

結合以上圖表和分析可知，基于統一耦合框架的堆芯物理-熱工耦合程序的計算結果是準確的.

表4 不同耦合方式計算得到的C2彈棒基準的穩態和瞬態結果

圖5 ICoCo和CORTHV2耦合模式計算得到的C1功率變化對比

圖6 ICoCo和CORTHV2耦合模式計算得到的C1平均多普勒燃料溫度對比

圖7 ICoCo和CORTHV2耦合模式計算得到的C1最高燃料溫度對比

圖9 ICoCo和CORTHV2耦合模式計算得到的C2平均多普勒燃料溫度對比

圖10 ICoCo和CORTHV2耦合模式計算得到的C2最高燃料溫度對比

5 結 論

(1) 本文基于ICoCo統一封裝理念和MED庫統一數據傳遞模型，分別對中子時空動力學程序CORCA-K和子通道分析程序CORTH進行了ICoCo封裝. 然后基于統一耦合框架實現了CORCA-K和CORTH的耦合.

(2) 以NEACRP輕水堆彈棒基準中C1和C2為研究對象進行耦合程序的測試，并與CORCA-K原有的耦合接口CORTHV2計算的結果進行對比，以基準結果作為參照. 經過對比分析可知，基于統一耦合框架的物理熱工耦合程序計算的結果與CORTHV2計算的結果符合地很好，各參數的變化趨勢完全一致，由此證明基于統一耦合框架的物理-熱工耦合程序的計算是正確的.

(3) 下一步可以進行耦合程序的敏感性測試，進一步確定其誤差范圍.