黃瑀琦，蔡杰進

(華南理工大學 電力學院，廣東 廣州 510640)

反應堆燃料組件是堆芯的主要部件，是堆芯熱工水力性能的重要影響因素之一。與研究周期長、費用高的試驗研究相比，CFD模擬已成為安全、快速研究手段之一，因此利用CFD進行燃料組件子通道的流動分析研究具有重要意義。燃料組件棒束子通道內湍流流動模擬問題是核熱工水力領域長期存在的問題之一，而決定模擬質量的核心便是其內部的算法代碼。

MATIS-H基準題是韓國原子能中心主導啟動的第二個檢驗CFD代碼的國際基準題試驗計劃(IBE-2)，其設計目的是檢驗CFD軟件代碼在三維子通道單相湍流中的模擬質量。MATIS-H基準題的格架是標準壓水堆燃料組件的5×5格架，試驗使用了兩種攪混翼，分別為分離式和旋流式。

Chang等[1]總結了MATIS-H的試驗數據，至2010年12月，共有25款商用和開源的CFD軟件代碼進行盲測，參與模擬測試的網格數從330萬至1.44億不等，近避面y+值從10-2到150不等。根據測試結果得出一個針對單相流動模擬的最佳預測指導(BPG)[2]，本文中的模擬設置參考了該最佳預測指導。Cinosi等[3]利用商業(yè)CFD軟件STAR-CCM+計算MATIS-H基準題，指出基于RANS模型求解的速度和壓力在峰值的預測上有偏差，低估了格架下游的湍流強度。Miku?等[4]在2012年新歐洲核能會議上提出OpenFOAM在反應堆熱工水力領域，尤其是子通道CFD模型代碼研究上具有廣泛的應用前景。其團隊針對不同工況下棒束子通道，在小尺度幾何模型下，同時進行WALE-LES模型和DNS模型的模擬。利用DNS模型數據驗證植入的WALE-LES模型，證明OpenFOAM具有模型植入能力[5]。商用CFD軟件為便于工程使用以及保護自身知識產權，其絕大多數模型的核心算法代碼均被封裝。用戶并不能閱讀和修改模型代碼。當商用軟件對特定問題進行模擬時(如單相棒束流動)，對部分物理量(如速度、壓力)的預測值較為滿意，但部分物理量(如湍流強度)預測值與試驗值偏差較大，且這種偏差并非來自于網格和收斂精度的影響。Chen等[6]使用商業(yè)軟件CFX v14.5研究了基于RANS方法的4種流動模型對子通道棒束流動的CFD模擬，在速度和壓力上能得到與試驗一致性較好的結果，但在湍流強度上卻和試驗數據相差2～5倍。用戶只能通過更換模型才能解決這種問題，但更換模型時需要重新考慮新模型的適用范圍和重新設置各種參數，甚至新模型在計算資源占用和耗時上是舊模型的幾倍甚至幾十倍，因此開源CFD軟件在近年來備受研究者的青睞。本文借助開源軟件OpenFOAM v8.0，對MATIS-H基準題的幾何模型進行全尺寸建模，模擬得到速度分布和壓力分布，并研究OpenFOAM對子通道CFD模擬的有效性。

1 5×5燃料組件建模

1.1 幾何建模

以MATIS-H基準題為標準建立幾何模型，攪混翼為分離式葉片，各圍板尺寸均為寬33 mm、長63.40 mm、厚度0.4 mm。定位格架的幾何結構復雜，因彈簧、剛突與燃料棒間是線接觸或弧面接觸，導致局部網格劃分難度大，網格質量差，目前的研究大多對定位格架進行簡化。本文的簡化方法是將彈簧、剛突與燃料棒間改為面接觸，并保留定位孔等結構。簡化后的格架如圖1所示。

圖1 簡化后的格架Fig.1 Simplified grid

1.2 網格劃分

通過OpenFOAM中的blockMesh工具實現外流域的劃分，snappyHexMesh工具進行細化網格。利用blockMesh工具生成1個65 mm×65 mm×520 mm長方體外流域。在前處理環(huán)節(jié)上，OpenFOAM允許用戶使用其他商用軟件進行網格劃分等前處理操作，本文為驗證OpenFOAM自身的模擬有效性，利用OpenFOAM具有的snappyHexMesh工具，僅在格架處設置細化區(qū)域，并將劃分精度等級設置為1級、2級、3級、4級，以便進行網格敏感性分析。其中，3級精度網格劃分如圖2所示。將不同精度網格所計算出的軸向平均速度與MATIS-H的試驗數據進行對比。

網格敏感性分析結果列于表1，可見，4級精度網格數是3級精度網格的2.1倍，但計算精度僅提升了1.33%，計算時長是其3.6倍，3級精度網格已能滿足計算要求，誤差小于5%。本文將采用3級精度網格進行模擬。

圖2 3級精度網格劃分Fig.2 Level 3 mesh

表1 網格敏感性分析Table 1 Sensitivity analysis of mesh

1.3 模擬條件設置

入口速度設為1 m/s，密度設為994 kg/m3，與MATIS-H試驗設置一致。出口選擇壓力法向梯度為0，壁面采用無滑移邊界處理，平均靜壓為0 MPa，殘差控制收斂精度設置為10-6，雷諾數約為1.10×105，參數設置列于表2。

求解器選用simpleFOAM。simpleFOAM是一種基于SIMPLE算法的求解器，SIMPLE算法由Patankar等[7]于1972年提出，并迅速成為計算不可壓縮流體的主流方法。

表2 參數設置Table 2 Parameter setting

2 結果分析

2.1 SST k-ω模型試算

在OpenFOAM下，設置如圖3所示4個截面和截面上的通道A作為數據獲取源。采用SSTk-ω模型在通道A上進行軸向速度場試算，將試算結果與MATIS-H基準題中的試驗數據進行對比。

截面3、4的通道A軸向速度如圖4所示，其中V為入口初速度，w為實際速度。截面3的通道A處于攪混翼下游，出現較為明顯的速度波動，OpenFOAM中實現的SSTk-ω模型低估了軸向速度波動幅度。截面4位于遠離攪混翼的下游，SSTk-ω模型的軸向速度預測效果較好，總體試算結果較為吻合，達到了商業(yè)CFD分析軟件的應用水平[8]。

2.2 SST k-ω模型模擬結果與討論

截面1位于格架下方15 mm處，截面2位于格架中部位置，距離格架入口15 mm處，截面3位于格架上方攪混翼下游出口，截面4位于遠離格架出口60 mm位置。各截面軸向速度v示于圖5。

圖3 橫截面選取和通道A示意圖Fig.3 Selection of plane and channel A

圖4 通道A軸向速度對比Fig.4 Comparison of axial speed of channel A

圖5 各截面軸向速度Fig.5 Axial velocity of each plane

由圖5可見，位于格架中部的截面2，由于單彈簧和雙彈簧實體的阻礙作用，彈簧附近的軸向流動速度較小，形成角通道流動速度大、棒束邊界流動速度小，受實體擠壓導致局部流速達到1.80 m/s。截面3位于攪混翼下游，軸向速度在攪混翼的作用下出現規(guī)則對稱的擾流，如圖6所示。在攪混翼下游出現對稱的軸向速度分布，中心通道的流速較角通道流速大，從流線圖可看出，在攪混翼下游出現渦流，這種渦流結構之間相互干涉，在攪混翼分布密集的中心通道，渦流之間干涉最為明顯。

各截面橫向流動速度u示于圖7。因位于攪混翼下游，受攪混翼的擾動影響，橫向流動速度開始增強，并規(guī)則地出現在攪混翼附近，形成1組大小相等、方向相反的橫向流動速度，且在攪混翼片兩側的流動速度最大，最大橫向流動速度達到入口速度的53%。位于遠離格架的截面4上，攪混翼的攪混效果減弱，最大橫向流動速度減小到入口速度的20%，但由攪混翼形成的橫向流動區(qū)域不斷擴大，在格架中心通道的橫向流場范圍最廣，這是由于中心通道的攪混翼片數目密度最大所致。

圖6 截面3的速度矢量圖和流線圖Fig.6 Vector and streamlinediagrams of plane 3

圖7 各截面橫向流動速度Fig.7 Cross flow velocity of each plane

各截面壓力分布云圖示于圖8。在格架中部壓力呈現中心通道高、四周通道較低，這是由于格架產生的束流效應使冷卻劑向中心流道靠攏，壓力升高。位于攪混翼下游的截面3上，經過攪混翼的攪混作用后，橫截面上的流體壓力分布發(fā)生很大變化。在靠近攪混翼處，流體壓力較小，在遠離攪混翼處，流體壓力較大。這是由于攪混翼加強冷卻劑擾流作用，使高壓區(qū)和低壓區(qū)的分布明顯。位于格架下游的橫截面4，經過一段規(guī)則的流動段，攪混翼對流體的擾動逐漸減弱，高壓區(qū)和低壓區(qū)相互擴散，右上角壓力與左下角壓力呈中心對稱，略高于中心通道壓力。

圖8 各截面壓力分布云圖Fig.8 Pressure distribution of each plane

3 結論

本文以MATIS-H試驗和文獻[8-9]中的數值模擬為事實基礎，利用OpenFOAM進行數值模擬。在攪混翼的作用下，攪混翼下游截面3上軸向流速最大達到1.28 m/s，相比于入口速度提升了28%，最大橫向流動速度達到入口速度的53%，中心通道冷卻劑攪混效果明顯。經數據對比得出，速度和壓力的模擬結果相比試驗有10%～20%的波動，屬于模型所產生的正常誤差范圍。壓力和溫度模擬趨勢和物理現象與試驗吻合度較高。OpenFOAM具備一定的單相棒束流動數值模擬能力，并能在標準求解器進行控制方程的植入，具有一定的研究前景。

本文僅是初步利用OpenFOAM進行的單相棒束流動模擬，未來還需進行多相流動的棒束流動研究，深入細致流場結構的研究和分析。