姜 榮, 馮文培, 陳紅麗, 強勝龍, 李治剛, 潘俊杰, 張喜林, 羅 曉
(1. 中國科學(xué)技術(shù)大學(xué)核科學(xué)技術(shù)學(xué)院, 合肥 230027; 2. 中國核動力研究設(shè)計院, 成都 610213)
反應(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é)果進行對比.
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)用程序示意圖
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)鏈接庫.
在得到封裝后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)格
選擇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ù)的描述.
通過編寫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ù)
使用基于耦合框架的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最高燃料溫度對比
(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) 下一步可以進行耦合程序的敏感性測試,進一步確定其誤差范圍.