RELAP5/FLUENT耦合程序的開發

何 帆，蔡翔舟，郭 威，何 龍，崔 蕾，趙 恒

(1.中國科學院 上海應用物理研究所，上海 201800；2.中國科學院大學，北京 100049； 3.中國科學院 先進核能創新研究院，上海 201800)

RELAP5是美國愛德華國家工程實驗室為美國核管會開發的輕水堆冷卻系統事故工況的瞬態行為最佳估算程序[1-3]，已被廣泛應用于反應堆的瞬態事故分析和安全評審等方面，并被核工業接受成為安全分析工具。在RELAP5/MOD4.0版本中，添加了鉛鉍合金、熔鹽等多種計算流體，其可應用于熔鹽堆等第4代核能系統的瞬態事故分析[4-5]。同時，隨著流體動力學方法(如商業化程序FLUENT)和計算機性能的高速發展，有關反應堆的三維數值模擬越來越多[6-8]，堆芯上下腔室、冷卻劑通道等局部構件計算流體動力學(CFD)模擬已被報道。以FLUENT為代表的CFD程序將會成為開發第4代核反應堆的一個強大工具。為綜合利用系統代碼和CFD程序的優點，國際上有學者已開始系統程序與CFD程序間的耦合研究，開發了一些系統/CFD耦合程序，并進行了初步驗證，如Aumiller等[9]首先利用并行虛擬機技術(PVM)實現了RELAP5-3D和CFD代碼的耦合，從原理上證明了RELAP5-3D/CFD顯式耦合是可行的；Papukchiev等[10]將ATHLET程序和CFX軟件耦合，模擬了流體在管道流動、傳熱等物理過程；Anderson等[11]針對超高溫氣冷堆冷卻劑流出堆芯進入出口腔室的熱混合問題，利用PVM耦合了RELAP5和CFD代碼，在出口腔室建立了三維流動效應模型。

本文討論顯式耦合RELAP5/MOD4.0和FLUENT方法，通過RELAP5/MOD4.0源代碼的二次開發和FLUENT用戶自定義函數(UDF)功能，RELAP5和FLUENT會分別在每個時間步開始時讀入耦合計算所需要的變量進行計算，并在該時間步結束時輸出耦合計算所需要的變量，然后再進行下一時間步的計算。本文為測試耦合程序的正確性，分別使用RELAP5/FLUENT耦合程序和獨立的RELAP5或FLUENT對水平圓管進行模擬分析。最后，利用耦合程序對2 MW熔鹽堆進行穩態工況模擬和功率突變的瞬態分析。

1 耦合方法

程序間的數據信息交換是耦合程序的本質[12]。在已有的一維系統程序和三維CFD程序耦合的研究成果中，主要是基于PVM來實現。在PVM中，額外編寫接口控制程序，利用PVM內部所提供的Send和Recv函數在每個時間步的計算中調用系統程序和CFD程序的邊界參數，并將這些參數值分別賦予相應的程序作為邊界條件以進行下一時間步的計算。但對RELAP5/FLUENT顯式耦合而言，這種方法相對較為復雜，且需較高的編程能力[13-14]。在本文中直接對RELAP5/MOD4.0源代碼進行二次開發輸入輸出接口模型，同時FLUENT基于UDF功能直接讀入和輸出邊界參數，這種方法相對簡單，也較易實現。

圖1為RELAP5和FLUENT在耦合計算時需實現交換的相關參數。FLUENT模擬區域的耦合出口邊界的輸出參數主要包括流體壓力p、空泡份額α、流體溫度T和質量流量W，將作為下游RELAP5計算區域的入口邊界參數；同時RELAP5模擬區域也需反饋流體壓力p、流體回流溫度T及空泡份額α來作為FLUENT壓力出口邊界條件。對FLUENT的速度入口邊界而言，需輸入上游RELAP5流體壓力p、空泡份額α、流體溫度T和流體流速v等參數，同時反饋給上游RELAP5包括流體壓力p、空泡份額α、流體溫度T等參數來作為RELAP5計算的邊界條件。

圖1 耦合程序交換參數示意圖Fig.1 Exchanged parameter of coupled program

(1)

(2)

其中：Ai為FLUENT邊界單元格的面積；n為時間步數。輸出的平均溫度采用質量加權的方式：

(3)

RELAP5整體結構如圖2所示，其主要分為3部分：輸入部分(INPUT)，瞬態/穩態計算部分(TRNCTL)，輸出部分(OUTPUT)。瞬態/穩態計算部分負責將輸入數據按照變量所屬類型進行相應的熱工水力、控制邏輯等方面的計算，并將計算結果保存在輸出文件中，它由TRNSET、TRAN和TRNFIN等子程序組成。其中，TRAN控制瞬態求解的步進，進行矩陣運算，是程序最重要的部分，也是程序計算過程中耗費時間最長的部分。在RELAP5瞬態計算過程中，TRAN子程序下的Dtstep子程序模塊控制時間步長的大小與計算結果輸出、繪圖編輯的頻率。在瞬態/穩態計算部分中，Dtstep子程序模塊中的Majout子程序控制大編輯輸出(圖3)，其特點如下：1) 該模塊全程參與瞬態/穩態部分的計算；2) 按照輸入卡設定將計算過程中相關變量的數據輸出到輸出文件中；3) 程序中含有大量的變量信息，并可按照一定格式將其輸出到輸出文件中；4) 該子程序模塊具備一定的擴展性，可添加變量，并進行賦值、邏輯判斷等操作。

圖2 RELAP5主程序3大部分Fig.2 Three parts of RELAP5 main program

圖3 瞬態/穩態計算部分的模塊結構Fig.3 Modular structure of transient/steady state calculation part

基于RELAP5的上述特點，本文將程序所需要的輸入輸出模塊放置在該子程序模塊中，即從外部(如FLUENT)讀取數據，并賦值給程序變量；或將程序變量的值輸出到外部。在熱工參數輸出功能上，RELAP5提供了小編輯輸出功能，通過設置輸出變量，可將數據輸出到輸出文件中，本文耦合程序動態輸出功能則借鑒了該特點。通過修改小編輯卡相關模塊，禁止小編輯卡輸出，確保小編輯卡不會對耦合程序的輸出產生干擾。同時，將小編輯卡所含變量的數據傳遞給輸出變量；從而完成小編輯卡的輸出變量向動態輸出變量的轉換。

如原RELAP5輸入卡301 mflowj 110010000命令，將會輸出編號為110的管道的質量流到RELAP5輸出文件中，程序修改后的該命令將直接把質量流輸出到外部。為了與原輸入卡小編輯輸出區分，同時擴展動態輸出數據的數量，本文進一步修改RELAP5小編輯模塊，在rmiedt子程序模塊中將輸入卡編號范圍從301～399修改為5 001～5 999；在Majout子程序中添加小編輯輸出代碼。程序修改后可通過輸入卡命令5001 mflowj 110010000直接將結果輸出到外部。

在RELAP5中，邊界條件主要是由時間相關變量來設置，包括時間相關控制體(TDV)、時間相關接管(TDJ)和隨時間變化的功率等。這些邊界條件參數主要包括液體/氣體的流速或流量，液體/氣體的溫度、壓力以及流體的空泡率等，其在輸入卡中以表格的形式進行輸入，表值在RELAP5讀入輸入文件后保存在對應的變量中，在計算過程中通過插值函數來獲得不同時刻的變量值。因此，可通過修改表值來實現修改邊界條件的目的。以RELAP5輸入卡110TDJ質量流為例：1100201 0.0 0.0，1100201 1.0 2.0，表示在0.0～1.0 s期間，110TDJ線性引入了2.0 kg/s的質量流量。若在計算中動態地將該表值0.0和2.0修改為3.0，則表示當前1.0 s時，110TDJ的質量流量為3.0 kg/s。其他邊界條件參數如流體的溫度、壓力等參數也可以通過這種方式輸入。通過這種動態修改表值的方式可實現將FLUENT邊界條件傳遞給RELAP5。

耦合程序計算流程示意圖如圖4所示。在耦合程序中，RELAP5和FLUENT分別讀入各自的輸入文件，并進行初始化。鑒于RELAP5是系統級程序而FLUENT多應用于局部構件的分析，RELAP5首先進行第1個時間步的計算更為適宜。RELAP5在第1個時間步計算結束后，將FLUENT所需要的邊界參數傳遞給FLUENT；FLUENT在進行該時間步的計算后將邊界參數傳遞給RELAP5，以便RELAP5進行下一時間步的計算。通過這種方式在每個時間步不停地更新邊界參數，直到最后1個時間步完成整個瞬態計算過程。

2 驗證及應用

2.1 flibe熔鹽在水平圓形管道流動問題

1) 圓管流動問題描述

首先利用flibe熔鹽在水平圓管流動問題來測試耦合程序。假設一水平圓管長1.0 m，流通面積為3.14 159×10-4m2，初始壓力為0.101 325 MPa，初始流速為1.0 kg/s，初始流體溫度為873.15 K，管壁粗糙度為10-6m且絕熱。整個瞬態模擬時間為10 s，水平圓管入口處flibe熔鹽的流速初始時刻為1.0 kg/s，在2～4 s期間線性增加至1.5 kg/s，6～8 s期間線性減少至1.0 kg/s并保持不變直到10 s結束。

圖4 耦合程序計算流程示意圖Fig.4 Schematic of coupled program calculation process

水平圓管入口處flibe熔鹽的溫度初始時刻為873.15 K，4～6 s期間溫度線性增加到893.15 K并保持不變直到10 s結束。為了保持一致性，FLUENT采用RELAP5中的flibe熔鹽物性。圖5為flibe熔鹽圓管流動的模型示意圖，分別采用RELAP5、RELAP5/FLUENT耦合程序及FLUENT對圓管進行模擬。

2) 結果與討論

RELAP5、RELAP5/FLUENT耦合程序及FLUENT在圓管中耦合界面1和耦合界面2處的質量流量和溫度變化分別如圖6、7所示。在整個瞬態模擬分析中，RELAP5/FLUENT耦合程序質量流量先保持不變，然后線性增加到1.5 kg/s后保持約2 s不變，最后線性減少到1.0 kg/s并保持不變。其變化趨勢符合預期，也與RELAP5和FLUENT單獨計算的結果保持了較好的一致性。同時，耦合程序的流體溫度先保持873.15 K不變；當圓管入口處溫度線性

圖5 flibe熔鹽圓管流動的模型示意圖Fig.5 Model schematic of flibe salt flow in round tube

圖6 耦合界面1(a)和耦合界面2(b)的質量流量Fig.6 Mass flow rate at coupled-boundary 1 (a) and coupled-boundary 2 (b)

圖7 耦合界面1(a)和耦合界面2(b)的流體溫度Fig.7 Fluid temperature at coupled-boundary 1 (a) and coupled-boundary 2 (b)

增加后，由于溫度熱傳導和流體流動需一定時間，流體溫度在耦合界面1和耦合界面2位置處會在t=4.1 s時刻溫度開始線性增加到893.15 K，并保持893.15 K不變。圖8為在耦合界面1和耦合界面2處RELAP5/FLUENT耦合程序與RELAP5、FLUNET的溫度差值隨時間變化的情況，其中，C-R表示RELAP5/FLUENT耦合程序和RELAP5的溫度差值；C-F表示RELAP5/FLUENT耦合程序和FLUENT的溫度差值。RELAP5/FLUENT耦合程序與RELAP5的溫度差值均小于0.01 K，與FLUENT的溫度差值也均小于0.30 K，相對誤差較小。其變化趨勢符合預期，與RELAP5和FLUENT單獨計算的結果保持了很好的一致性。

圖8 耦合界面處耦合程序與RELAP5、 FLUNET的溫度差值Fig.8 Temperature difference between coupled program and RELAP5, FLUNET at coupled-boundary 1 and coupled-boundary 2

RELAP5、RELAP5/FLUENT耦合程序及FLUENT在耦合界面1和耦合界面2位置處的管道中流體壓力變化分別如圖9所示。在t=2～4 s過程中，由于flibe熔鹽流速增加，管道中流體壓降變大，在出口壓力不變的情況下，管道中流體的壓力會逐漸變大；在t=4～6 s過程中flibe熔鹽溫度增加，黏性降低導致壓降變小，在耦合界面1和耦合界面2處中流體的壓力會降低；在t=6～8 s過程中，由于flibe熔鹽流速降低，管道中流體壓降變小，管道中流體的壓力會逐漸變小。RELAP5、RELAP5/FLUENT耦合程序和FLUENT在不同時刻耦合界面1和耦合界面2之間的壓降列于表1。RELAP5/FLUENT耦合程序與FLUENT計算的壓降比較相符，在t=4 s時最大相對誤差為2.0%，在可接受的范圍內。RELAP5/FLUENT耦合程序和FLUENT計算的壓降，與RELAP5計算的壓降差距較大，這是由于RELAP5基于流體充分發展情況下計算流體壓降而RELAP5/FLUENT耦合程序和FLUENT采用低Re湍流模型，在進出口效應的影響下，耦合界面1處的流體尚未充分發展，導致壓降計算相對偏大。總體而言，RELAP5/FLUENT耦合程序與RELAP5、FLUENT計算的壓力變化比較一致。

2.2 2 MW熔鹽堆穩態模擬

1) 2 MW熔鹽堆簡介

液態燃料熔鹽堆基于在線干法處理技術,可實現釷鈾增殖,是一種重要的研究堆型[15]。本文以一種2 MW液態燃料熔鹽堆的初步概念設計為研究對象。該液態熔鹽堆設計熱功率為2 MW，堆芯采用石墨作為慢化劑，熱工水力設計堆本體熔鹽的進口溫度為873.15 K，出口溫度為893.15 K。一回路采用含有高富集度7Li的LiF-BeF2-ThF4-UF4熔鹽作為燃料鹽，二回路采用FLiNaK作為載熱劑，一回路質量流量為55.4 kg/s。在利用RELAP5/FLUENT耦合程序模擬分析時，堆芯部分選取了1/4模型進行CFD建模，堆芯外的系統回路采用RELAP5建模，如圖10a所示。下腔室設計采用喇叭狀結構而流量分配裝置處于下腔室正中心，其厚度為18 mm，其結構如圖10b所示。堆芯熔鹽通道1/4模型由圖11a所示的堆芯石墨組件構成192個熔鹽通道，熔鹽通道的分布如圖11b所示。由于實驗堆堆芯幾何結構較為復雜，對實驗堆幾何模型進行了簡化處理，CFD建模計算時省略了流量分配裝置支撐連接結構。熔鹽堆部分材料的物性參數列于表2。在CFD模型中，對幾何結構較為規則的堆芯活性區，采用六面體結構化網格進行劃分；對于堆芯上下熔鹽腔室為橢球形結構，采用非結構化四面體網格進行劃分。網格質量最差為0.45，約85%的網格質量控制在0.8以上。網格無關性分析表明，當計算模型網格數量到達3 000萬及以上時，可滿足CFD精確計算要求。

圖9 耦合界面1(a)和耦合界面2(b)的流體壓力Fig.9 Fluid pressure at coupled-boundary 1 (a) and coupled-boundary 2 (b)

表1 熔鹽在耦合界面1和耦合界面2之間的壓降Table 1 Salt pressure drop between coupled-boundary 1 and coupled-boundary 2

2) 結果分析

利用RELAP5/FLUENT耦合程序對2 MW熔鹽堆系統進行分析，相對單獨的RELAP5系統分析結果而言，耦合程序可通過FLUENT模型計算獲得堆芯更加詳細的溫度分布和流場分布。堆芯的溫度分布如圖12所示，熔鹽在流經堆芯的過程中，溫度逐漸上升并在堆芯活化區出口處達到最大值約905 K，隨后熔鹽在上腔室進行混合，在出口處溫度約893 K。

圖10 熔鹽堆模型示意圖(a)與流量分配裝置結構(b)Fig.10 Schematic of molten salt reactor (a) and structure of flow distribution device (b)

圖11 石墨組件(a)和熔鹽通道分布示意圖(b)Fig.11 Schematic of graphite assembly (a) and salt channel distribution (b)

表2 物性參數Table 2 Physical parameter

圖12 堆芯溫度分布示意圖Fig.12 Temperature distribution of reactor core

熔鹽堆下腔室流線示意圖如圖13所示，熔鹽由較窄的進口管進入較寬闊的下腔室，熔鹽經歷一突然擴展的流態，流線形狀近似于喇叭狀結構，有利于熔鹽的流動，避免了熔鹽在下腔室的渦流現象。下腔室設置流量分配板，有效抑制了熔鹽在下腔室中心區域的流速，使得熔鹽流體不是直接由進口沖擊中心位置及其附近熔鹽通道，而是流經分配板和其孔道向四周擴散流動。因而下腔室內由渦流產生的流動死區基本得到消除，在堆芯徑向流量分配相對均勻。

在堆芯活性區不同高度(z=0.0，0.4，0.8，1.1 m)下堆芯截面溫度分布如圖14所示。由于下腔室喇叭狀結構和流量分配板的作用，熔鹽在進入堆芯活性區(z=0.0 m處)各通道的流速和溫度相對分布均勻。熔鹽在流經堆芯活性區過程中，熔鹽溫度逐漸上升。在同一高度下，堆芯中心區域熔鹽溫度略高于堆芯外圍熔鹽溫度。堆芯石墨受到熔鹽加熱，溫度會略比通道中的熔鹽溫度高。在z=0.0 m處堆芯石墨溫度約在878～882 K，隨著堆芯熔鹽加熱和石墨自身功率加熱，堆芯石墨溫度逐漸上升，在z=1.1 m熔鹽出口處，石墨溫度達到894～906 K。

a——0.0 m；b——0.4 m；c——0.8 m，d——1.1 m圖14 在不同高度下堆芯活性區截面溫度分布Fig.14 Temperature distribution of active core at different heights

2.3 2 MW熔鹽堆功率突變瞬態分析

功率突變瞬態分析是在上述穩態分析的基礎上進行的。假定在t=0 s時堆芯功率突然增加10%并保持不變，熔鹽堆系統無其他安全設施引入。考慮到熔鹽堆系統的復雜性及計算機計算性能的限制，整個瞬態分析過程中時間步長設置為0.01 s，模擬時間為200 s。突然引入的功率變化會導致堆芯和系統的溫度分布發生變化。耦合程序模擬分析的系統質量流量和堆芯進出口溫度隨時間的變化如圖15所示。在瞬態計算中，一回路系統的質量流量為55.4 kg/s保持不變；由于熔鹽具有較大的比熱容及在堆系統中的滯留效應，堆芯進口溫度在前期基本不變，出口溫度逐漸上升。在t=20 s時刻，堆芯出口的平均溫度為896.05 K，堆芯進出口溫差為22.9 K。由于二回路換熱能力有限，流出堆芯的熔鹽溫度逐漸上升后，在流經換熱器冷卻后溫度會高于873.15 K，且使得堆芯入口溫度在約t=30 s后開始緩慢上升，在t=200 s時堆芯進口溫度達到875.0 K。

圖16為瞬態計算過程中，在t=5、10、15、20 s時刻堆芯的溫度分布。在初始時刻穩態情況下，堆芯內部最高溫度約906 K。由于堆芯功率的增加，堆芯內部的最高溫度也發生變化。在t=5、10、15、20 s時刻，堆芯內部最高溫度分別達到909.0、910.4、910.7和911.0 K。初始時刻，堆芯下腔室熔鹽平均溫度在874 K，局部區域平均溫度會達到875 K；在t=20 s時刻，由于堆芯功率增加，下腔室熔鹽平均溫度增加，局部區域溫度可達879 K；堆芯內石墨溫度也相對升高，靠近堆芯上腔室附近的石墨平均溫度可達903 K。

圖15 系統質量流量和堆芯進出口溫度的變化Fig.15 Mass flow rate of system and temperature of reactor core inlet and outlet

a——5 s；b——10 s；c——15 s；d——20 s圖16 不同時刻堆芯截面溫度分布Fig.16 Temperature distribution of reactor core at different time

3 結論

以RELAP5與FLUENT程序為基礎，利用對RELAP5源代碼的二次開發和FLUENT的用戶自定義函數進行編程，采用顯式耦合的方式，開發了RELAP5/FLUENT耦合程序，實現兩個代碼之間數據交換。RELAP5在每個時間步進行瞬時計算后，返回FLUENT所需的邊界條件參數，FLUENT在該時間步結束后返回RELAP5下一步計算所需的邊界條件參數，重復此過程，直到達到預定問題的計算時間為止。利用RELAP5/FLUENT耦合程序對2 MW熔鹽堆系統進行數值模擬分析，在進行系統分析的同時可獲得堆芯詳細的三維熱工水力效應。該耦合程序可適用于包括flibe熔鹽等多種流體的管道流動情況下的熱工水力分析；在進行熔鹽堆系統分析中，堆芯及上下腔室等存在三維復雜結構，局部區域存在明顯的三維流動或需要獲得三維溫度分布，也可基于RELAP5/FLUENT耦合程序進行熔鹽堆的熱工水力分析。