基于統(tǒng)一耦合框架的堆芯物理熱工耦合程序的開發(fā)及驗證

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

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

1 引 言

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

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

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

2 理論依據(jù)

ICoCo接口是在NURISP項目范圍內(nèi)[8]開發(fā)的面向?qū)ο蟮耐ㄓ么a耦合接口，首先應(yīng)用于CATHARE和TrioCFD的多尺度熱工水力耦合[9]. ICoCo接口基于C++開發(fā)，通過定義Problem通用母類來定義函數(shù)，這些函數(shù)的功能包括初始化、時間步推進、保存和還原以及場數(shù)據(jù)傳遞，表1給出了ICoCo定義的重要接口及功能.

表1 重要的ICoCo接口及功能

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

為了實現(xiàn)基于統(tǒng)一耦合框架的物理-熱工程序的耦合，采用C++ 編寫的supervisor程序?qū)崿F(xiàn)對各子程序的調(diào)用. 圖1給出了C++ supervisor進行程序調(diào)用的示意圖.

圖1 supervisor調(diào)用程序示意圖

3 程序開發(fā)

3.1 CORCA-K的ICoCo封裝

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

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

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

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

3.2 耦合實現(xiàn)

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

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

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

CORTH非結(jié)構(gòu)化網(wǎng)格

4 數(shù)值驗證

4.1 測試算例

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

4.2 封裝后CORCA-K測試

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

圖4 輕水堆基準(zhǔn)C1和C2徑向堆芯布置圖

表2 C1和C2的堆芯運行數(shù)據(jù)

4.3 耦合程序測試

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

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

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

表3 不同耦合方式計算得到的C1彈棒基準(zhǔn)的穩(wěn)態(tài)和瞬態(tài)結(jié)果

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

結(jié)合以上圖表和分析可知，基于統(tǒng)一耦合框架的堆芯物理-熱工耦合程序的計算結(jié)果是準(zhǔn)確的.

表4 不同耦合方式計算得到的C2彈棒基準(zhǔn)的穩(wěn)態(tài)和瞬態(tài)結(jié)果

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

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

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

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

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

5 結(jié) 論

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

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

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