堆芯燃料組件流阻特性模塊式與整體式三維數值模擬方案比較

楊紅義，周志偉，林 超，李 淞，王予燁

(中國原子能科學研究院 反應堆工程技術研究部，北京 102413)

快堆組件中相同功能組件按照其具體形狀及所承擔的具體功能任務分為不同類型的組件。每類組件對應的功率有所差異，故在反應堆運行時，不同類型組件對應入口流量不同。若燃料組件分配到的流量偏大，則堆芯出口處冷卻劑鈉的溫度偏低，進而影響汽輪機參數，使得快堆經濟性變差；若組件分配到的流量偏小，可能導致燃料組件包殼溫度升高甚至超過極限溫度，威脅快堆的安全。因此，須對燃料組件進行合理的設計和計算使其能匹配不同流量的冷卻劑。而組件的阻力特性是影響堆芯流量分配至關重要的因素，所以建立合理的阻力特性計算方案成為當下亟待研究的問題。

堆芯包含大量組件，每盒組件又包含數十根燃料棒，因此要對整個堆芯進行三維流動換熱計算模型的建模是極其不現實的。為方便示范快堆設計，中國原子能科學研究院堆芯及一回路項目組開發了適用于示范快堆熱工水力設計及優化的程序，該程序需獲得單個組件的流阻特性，以準確預估堆芯的流量分配是否合理。因此單個組件流阻特性的預測精度將很大程度上決定堆芯流量分配預測的準確性。

丁振鑫[1]對CEFR的Ⅰ-Ⅱ型小柵板聯箱的流動阻力水力特性進行了實驗。馮預恒等[2-4]則對 CEFR Ⅰ-Ⅱ型小柵板聯箱及其節流件進行了穩態模擬計算，給出了節流件內的壓力和速度分布，并與文獻[2]中的試驗結果進行對比，計算結果驗證了三維數值模擬結果的可靠性，證明其完全可在核反應堆設計(至少在熱工流體領域)中使用。王雨等[5]利用CFX對七根式小柵板聯箱中快堆燃料組件管腳段均勻開口和非均勻開口兩種模式下的速度場和壓力場分布進行了數值模擬。李淞等[6]對169棒束的快堆燃料組件稠密棒束進行了三維數值模擬，根據模擬結果可更準確地預測棒束通道內的流動情況，可為今后稠密棒束組件水力學設計和子通道內流量測量試驗提供參考。國外對池式鈉冷快堆的流動特性也進行了大量的實驗和數值模擬研究[7-8]。Padmakumar等[9]使用縮比模型對印度原型快堆(PFBR)進行實驗，給出了反應堆內冷卻劑流動分析、溫度分布和流動振動特性。Péniguel等[10]對鈉冷快堆冷卻劑的三維流動特性和耦合傳熱進行了分析，但同樣缺乏對局部壓力分布及整體壓降的模擬計算。燃料組件壓降主要由3部分構成：管腳段壓降、六角管棒束區壓降和操作頭段壓降。上述相關研究都是對燃料組件各部分進行單獨研究，目前尚未發現關于燃料組件全組件的流阻特性計算。戚展飛等[11]對帶孔洞的復雜套管在不同冷卻劑流動方向下的三維管段單相流動進行了數值模擬，研究表明，管道局部流動截面相對于流動方向的改變方式和改變順序，是產生不同局部壓降的重要原因。各組件部分之間過渡段也存在可觀局部壓降，分開計算并不能表征這部分的壓降，同時分開計算也忽略了各部分之間流動的互相作用。為得到更為精確的單組件流阻特性，對燃料組件模擬件進行全組件三維數值模擬研究至關重要。

本文選用燃料區某兩類組件作為研究對象，分別對其進行整體式和模塊式三維數值模擬，獲得其流阻特性并與試驗結果進行對比，通過分析兩種方法的精確度并結合計算資源耗費，綜合評估兩種分析方法的優勢。然后深入研究速度和溫度對兩類組件流阻特性的影響，并分析這些因素對兩種方案預測精度的影響。

1 模型介紹

本文以堆芯燃料組件為研究對象，選取其中兩類燃料組件進行流阻特性研究。燃料組件結構主要包括操作頭、六角管棒束段和管腳。燃料組件的棒束區與操作頭結構一致，通過改變管腳區域的結構來改變整個組件的流阻特性。兩類燃料組件的具體結構如圖1所示，其中二類燃料組件在管腳(二類管腳)區域設置節流件來增加管腳區域的流動損失。

圖1 燃料組件結構Fig.1 Structure diagram of fuel assembly

2 試驗方法

試驗模擬件與組件原型1∶1幾何相似，且加工精度也與設計施工圖紙要求保持一致。對于燃料組件，要求確保試驗組件在水作為工質時的雷諾數與鈉作為工質時的雷諾數相等，這意味著試驗組件在水中及液態鈉中應保持流速與運動黏度相同，進而充分模擬示范快堆實際工況中鈉的流動狀況與阻力特性，使得水力試驗結果可應用于示范快堆的各項設計和驗證。

沿燃料組件模擬件長度方向上冷卻劑鈉平均溫度TNa約為420 ℃。由于水的運動黏度ν水與該溫度下鈉的運動黏度ν鈉相等，可得到試驗組件對應水溫為93.2 ℃。由于兩種工質的雷諾數和運動黏度均相等，可得到二者具有相

同的流速，進一步可得通過試驗組件所要求的水工質體積流量Q水應等于組件平均溫度下鈉工質體積流量Q鈉，即：

lNa=l水，νNa=ν水

VNa=V水

Q水=AV水=AVNa=QNa

同樣，對于試驗組件，進口鈉溫TNa為358 ℃時，可得到試驗組件對應水溫為84 ℃。

試驗條件下采用恒定溫度的去離子水作為試驗工質，無法完全模擬實際的變溫過程，因此需要在兩個溫度下分別進行試驗，并進行溫度折算。首先，在84 ℃與93.2 ℃進行組件無管腳試驗，可在消除棒束區溫差所導致測量誤差的同時獲得最貼近實際情況的棒束段溫度；然后，在84 ℃下開展帶管腳的全組件試驗，獲得相應的全組件壓降；最后，84 ℃對應的管腳壓降與93.2 ℃對應的棒束段壓降構成組件的整體壓降。

3 計算方法

3.1 網格模型

根據結構的特性，將燃料組件分成4個區域進行網格劃分：管腳(兩類)、棒束及操作頭。不同類型燃料組件的棒束和操作頭結構相同，只有管腳不同。對各區域網格以壓降為標準，進行網格無關性驗證。其中棒束區由于網格規模較大，棒束根據繞絲單個繞距為單位具有可復制性，因此本文取1個繞距進行網格無關性驗證，結果如圖2所示。

根據網格無關性驗證結果，最終各區域具體網格數列于表1，具體網格劃分如圖3所示。

圖2 網格無關性驗證Fig.2 Grid independence verification

表1 各區域網格數Table 1 Grid type and number of each assembly

3.2 數學模型

燃料區組件管腳位于燃料棒下方，并無加熱源，考慮到液體介質為鈉液，可認為是動力黏度μ為常數的絕熱不可壓縮流體。因此，流動計算的控制方程可簡化為以下兩類方程[12]。

連續方程：

圖3 各區域網格Fig.3 Grid diagram of each assembly

N-S方程：

目前，對各流體微團的紊亂運動進行平均處理加湍流模型求解方式仍是全三維湍流N-S方程求解的主要方法。在湍流模式中雷諾平均法在工程中的應用非常廣泛，其基本思想是以雷諾假設為基礎，將湍流的物理量用平均值和脈動值的和來表示，湍流場被視為平均運動場和脈動運動場的疊加[13]。本文選用目前工程中應用最為廣泛的兩方程模型作為湍流模型，其中兩方程模型中RNGk-ε模型來源于嚴格的統計技術，它和標準k-ε模型很相似，但做了部分改進，使得RNGk-ε模型較標準k-ε模型在更廣泛的流動中有更高的可信度和精度。本文選用RNGk-ε雙方程湍流模型使方程組封閉。RNGk-ε模型具體方程[14]如下：

采用有限體積法對控制方程在空間上進行離散，選用SIMPLE算法對流場的離散方程進行求解。

3.3 邊界條件選取

1) 整體式與模塊式方案對比

為保證模擬條件和試驗一致，本文選取全組件試驗條件進行仿真計算，以去除數據處理過程產生的誤差。該條件下介質為84 ℃的去離子水，此時水的密度為969.4 kg/m3，動力黏度為3.393×10-4Pa·s。選取進口流量為53.8 m3/h的工況進行計算和詳細對比，進口為質量流量入口，出口采用壓力出口邊界條件。

2) 速度對壓降的影響

為保證計算結果的可靠性，需用試驗結果進行對比。本文選用整體式方案，并參照一類組件的全組件試驗條件，選取介質為84 ℃的去離子水，進口為質量流量入口，出口采用壓力出口邊界條件，研究速度對壓降的影響。速度與相應的體積流量與實驗中部分值對應(表2)。

3) 溫度對壓降的影響

為保證計算結果的可靠性，需用試驗結果進行對比。計算參照無管腳試驗，計算值包含棒束及操作頭區域，選用模塊式方案進行計算。選取進口流量為53.8 m3/h的工況為基準，分別選取84 ℃和93.2 ℃進行計算和詳細對比，進口為質量流量入口，出口采用壓力出口邊界條件。84 ℃去離子水的密度為969.4 kg/m3，動力黏度為3.393×10-4Pa·s；93.2 ℃去離子水的密度為963.4 kg/m3，動力黏度為3.03×10-4Pa·s。

表2 邊界條件Table 2 Boundary condition

3.4 計算方案

二類燃料組件只有管腳結構不同，其余部件相同。因此計算區域可分為4個模塊：一類管腳、二類管腳、棒束及操作頭。本文對兩種方案進行對比，并用試驗結果進行驗證。

1) 模塊式方案

模塊式方案對每個模塊進行單獨計算，獲得每個區域的壓降，通過各部分壓力損失疊加獲得燃料組件整體壓降，如圖4所示。

圖4 模塊式方案Fig.4 Modular calculation scheme

模塊式計算的計算量相對較小、周期短、靈活性高，針對不同類型燃料組件，棒束和操作頭的結構相同，具有可復制性，不用重復計算，只需對不同類型管腳進行計算即可。但模塊式計算忽略了各模塊間流動的相互作用的影響。

2) 整體式方案

整體式方案將4個模塊組合成2組，將各區域通過interface進行連接，最終合成整體進行計算(圖5)。

圖5 整體式方案Fig.5 Integral calculation scheme

整體式計算考慮各模塊之間流動的相互作用的影響，能更為精確地捕捉流動損失，獲得的整體壓降也更為精確。但整體式計算的計算量大、周期長、靈活性較低。

4 計算結果分析

4.1 模塊式與整體式三維數值模擬方案對比

1) 流阻特性

對模塊式和整體式兩種方案的整體壓降進行統計，并與試驗結果進行對比，結果列于表3。從表3可看出，整體式計算的結果更為接近試驗值。其中，一類組件整體式相對誤差較模塊式高4.9%，二類組件高1.5%。

對各子區域的壓降進行對比，可獲得不同計算方法影響最大的區域(表4)。可看出，絕對誤差最大的區域為棒束區，相對誤差影響最大的為操作頭區域，而影響最小的為管腳區域。

表3 總壓降對比Table 3 Comparison of total pressure drop

表4 各區域差值Table 4 Difference of pressure drop in each region

注：1) 絕對誤差為各區域壓降差值占組件整體壓降試驗值的百分比

2) 相對誤差為相對各區域壓降差值占各自區域壓降整體式計算的百分比

由此可推測，整體式計算法對上游區域的流態較對下游的影響大，而下游對上游幾乎無影響，所以位于最上游的管腳區域的整體式和模塊式計算的壓降幾乎一致。而位于最下游的操作頭區域雖然絕對誤差不大，但相較其自身壓降的相對誤差接近30%。

對存在壓降差值較大的區域進行分析，圖6為距棒束區進口不同高度截面的軸向速度分布。可看出，模塊式進口由于給定了流量進口，因此軸向速度分布均勻，而整體式由于上游管腳區域的影響，軸向速度存在高低速區。這種趨勢直至截面距離進口高度為0.05 m時，兩種計算模式的速度分布才逐漸趨于一致。因此在z<0.05 m區域兩部分流動損失不一致，整體式由于速度分布不均引起的損失更大。

圖7為操作頭中心截面流線圖，可看出整體式由于上游影響，入口流量存在偏離流動，在操作頭中間部分存在旋渦。但由于操作頭本身壓降占總體壓降的比例較小，因此該部分損失占總壓降不到1%，可忽略不計。

2) 計算資源對比

使用計算資源為32核、256內存計算機，最多可計算2億量級網格量，計算收斂時間約為120 h。

不同燃料組件只有管腳類型不同，在相同流量情況下，模塊式計算法可將2個方案分解為4個模塊，將4個模塊的壓降進行組合，即可得到整體壓降。4個模塊的網格量總和小于2億，可同時計算，計算時間約為120 h。

整體式計算時每個計算模型網格大于1.7億，2個計算模型總網格量為3.4億，需分2次進行，計算時間約為240 h。

圖6 棒束區z截面軸向速度分布Fig.6 z-section axial velocity distribution in bundle region

圖7 操作頭x截面流線圖Fig.7 Streamline diagram in x-section of operating head

針對本文研究對象，模塊式計算較整體式計算可節省一半時間。考慮到后續研究，堆芯共6種類型燃料組件，模塊式計算所需時間為整體式計算的1/6。

4.2 速度對壓降的影響

用模塊式和整體式分別對一類組件和二類組件進行計算，并與試驗結果進行對比，獲得速度對壓降的影響(表5)。可看出，隨著速度的上升，壓降上升；速度越大，整體式及模塊式計算結果與試驗值間的誤差越大。

表5 速度對壓降的影響Table 5 Effect of velocity to pressure drop

4.3 溫度對壓降的影響

由于全組件試驗只對84 ℃去離子水進行測試，而無管腳試驗對84 ℃和93.2 ℃分別進行測試，為保證計算結果的可靠近，選取無管腳試驗所包含區域為研究對象，即計算模型只包括棒束區和操作頭兩部分。模塊式和整體式計算所得的壓降與試驗結果的對比列于表6。從表6可看出，由于溫度越高去離子水的密度越低，在進口體積流量一定的情況下，溫度越高壓降越小。整體式預測壓降的相對誤差約為3%，模塊式預測壓降的相對誤差約為9%。溫度升高，預測精度小幅提升。

表6 不同溫度介質對壓降的影響Table 6 Effect of various temperature medium to pressure drop

注：進口流量為53.9 m3/h

5 結論

針對燃料區組件流阻特性預測精度問題，本文選用兩種計算方案進行計算預測，并與試驗結果進行對比，為后續研究提供依據。結果表明：

1) 整體式計算結果較模塊式更為精確，整體式較模塊式精度提高了4.9%及1.5%。

2) 針對本文研究對象，模塊式計算較整體式計算可節省一半時間，堆芯共6種類型燃料組件，模塊式計算時間為整體式計算的1/6。

3) 隨著速度的上升，壓降明顯升高，且整體式和模塊式預測精度隨速度的上升略有降低。

4) 溫度越高壓降越低，數值計算預測精度越高。