冷卻劑與含格架燃料棒束多場耦合分析

馬玉琢，黃 美，黃繼緣，張沛健

(華北電力大學 核科學與工程學院，北京 102206)

壓水堆核電廠通過冷卻劑流經燃料組件進行熱量交換將核裂變產生的能量帶出堆芯，冷卻劑的流動會對燃料組件產生沖擊與振動，使組件產生位移和形變，與此同時組件結構形狀和位置的改變會反過來對冷卻劑的流動產生顯著影響。冷卻劑與燃料組件間的這種相互作用直接關系到反應堆的安全運行與核電廠的運營效率。

在構成燃料組件的燃料棒束結構中，定位格架是用于支撐堆芯燃料棒束的彈性結構構件，對燃料組件的熱工和結構特性影響較大。利用計算流體力學(CFD)工具對燃料組件內單相和兩相流動現象的研究已十分豐富[1-3]，對組件結構與冷卻劑相互作用進行的流固耦合研究也日漸增多[4-6]，但目前的研究大多僅限于單向流固耦合，對涉及到結構場、溫度場和流場的雙向多場耦合情況研究較少。

本文采用數值模擬計算中的多場耦合分析方法對冷卻劑與含格架燃料棒束進行多場耦合模擬，分析研究冷卻劑與燃料組件間的熱工水力特性和結構形變特性，進而為堆芯的安全分析和燃料組件的優化設計提供參考。

1 計算模型

1.1 幾何模型

本文通過對實際燃料組件的適當簡化，忽略對流場影響輕微的固定件，保留定位格架中的條帶彈簧、剛性突起和攪混翼片，建立了總長1 050 mm的5×5燃料棒束計算域。其中：入口段長度500 mm，定位格架段長度50 mm，出口段長度500 mm。燃料棒束長1 050 mm，外徑9.5 mm，柵距12.6 mm；定位格架尺寸63.4 mm×63.4 mm×38.38 mm，厚度0.425 mm，定位格架上的條帶彈簧高1.4 mm，剛性突起高1.18 mm，攪混翼片高5.38 mm，彎折角度30°[7]。燃料棒束與格架的幾何模型如圖1所示。

1.2 網格劃分

基于上述幾何模型，本文在入口段與出口段使用結構化網格，格架段由于結構較為復雜選用非結構化網格；格架段與入口段和出口段間通過交界面連接，3部分共同構成混合網格，如圖2所示。對網格無關性進行了驗證，網格參數列于表1[8]。當單元數達到485萬時，各項參數滿足計算需要，綜合計算能力的考慮本文選用2號網格。

圖1 棒束與格架的幾何模型Fig.1 Rod and grid geometric model

圖2 混合網格結構Fig.2 Hybrid mesh structure

表1 網格參數統計Table 1 Mesh parameter statistics

1.3 計算條件

在帶有定位格架棒束通道幾何條件下，冷卻劑流過定位格架附近時受壁面黏性作用，因而湍流模型對壁面處理方法的不同對模擬結果有較大影響。參考現有研究成果[9-10]，對比k-ε、RNGk-ε和SSTk-ω3種湍流模型，考慮了旋流效應的SSTk-ω模型計算結果在預測壓降和溫度分布時更為準確，因此本文在流體計算部分采用SSTk-ω模型。

系統初始溫度292 ℃，壓力15.5 MPa，入口邊界初速度4.8 m/s，出口為壓力邊界條件[11]；燃料棒為粗糙壁面，根據堆芯功率密度分布特性燃料棒線熱流密度沿軸向呈余弦分布[12]；其他結構為無滑移、絕熱固定壁面，燃料棒兩端固定約束邊界，定位格架為彈性約束條件。

2 數值模擬結果與分析

根據流道所在位置的不同可將冷卻劑流動的子通道分為中心通道、邊通道和角通道，并將燃料棒按位置分為中心棒、中層棒和外層棒。同時為便于模型對比和參數分析，分別沿燃料棒束軸向和徑向選取4個截面，并在中心流道內選取一點，截面和點的位置及編號如圖3所示。

2.1 數值模擬可行性驗證

為驗證數值模擬的可行性，選取Holloway實驗布置作為模擬參數輸入，模擬與實驗對照如圖4所示[13]，結果表明歸一化努塞爾數隨冷卻劑流動的變化趨勢一致，模型能較好地模擬實際工況。

圖4 數值模擬可行性驗證Fig.4 Feasibility verification of numerical simulation

2.2 冷卻劑流動與傳熱的數值模擬

對冷卻劑流動與傳熱進行數值模擬計算，燃料棒線熱流密度沿軸向呈余弦分布，平均熱流密度為8.473×105W/m2，冷卻劑沿燃料棒束軸向各截面(Plane1～4)上的二次流速度及溫度分布云圖如圖5所示。由圖5可見：冷卻劑自下而上流動的過程中，從入口段到定位格架前幾乎沒有橫向流動的產生，棒束與冷卻劑間的換熱較弱，溫度梯度較大；冷卻劑流入格架段后由于格架壁面、條帶彈簧和剛性突起等結構的存在，開始出現流道間的橫向流動，繼續流動至格架出口位置，由于攪混翼片對流動的明顯擾動使橫向流動進一步顯著增強，中心位置的二次流速度可達3.1 m/s，橫向流動的攪混作用產生的渦流促進了不同位置和不同流道之間的熱量交換，使溫度分布趨于均勻，避免近棒束壁面溫度過高，并對定位格架下游出口段的流動和傳熱產生影響。橫向二次流一方面增強了冷卻劑與燃料棒束間的熱量交換，另一方面也會對燃料棒束結構產生沖擊，這種影響將在下述多場耦合模擬中進行分析。

圖5 軸向截面冷卻劑二次流速度及溫度分布云圖Fig.5 Contour of secondary flow velocity and temperature distributions at axial-intersection for coolant

圖6示出燃料棒束徑向各截面(Plane5～8)二次流速度分布云圖及局部子通道速度矢量圖。由圖6可見：由于定位格架內部條帶彈簧和剛性突起對冷卻劑流動的擾動作用并不強烈，冷卻劑在定位格架出口(Plane5，z=33 mm)之前雖出現了橫向流動，但并不明顯；冷卻劑從定位格架流出(Plane6，z=39 mm)后受到攪混翼片的擾動出現強烈的二次流，最大流速達到約3.8 m/s，對應Plane6局部速度矢量圖可見形成明顯的流動旋渦；隨流動繼續發展，多個流動旋渦逐漸融合，環繞燃料棒束形成對角線方向更大的繞流(Plane7，z=60 mm)，對應Plane7局部速度矢量圖可見旋渦的跨棒束分布，直至超過定位格架出口約20倍水力直徑(Dh)處繞流的影響才逐步消失(Plane8，z=120 mm)。

2.3 多場耦合條件下冷卻劑流動的數值模擬

上述冷卻劑流動和傳熱的數值模擬分析是在固體結構部件完全不變的條件下進行的。實際上由于冷卻劑流動產生的沖擊和振動，燃料棒束和定位格架會產生相應的形變和移動，特別是定位格架結構中的條帶彈簧、剛性突起和攪混翼片形狀與位置的改變會反過來對冷卻劑的流動特性產生顯著影響，進而影響到冷卻劑與燃料組件間的傳熱與溫度場的分布。因此需通過多場耦合分析流體與固體間的相互影響。

本文采用ANSYS WORKBENCH軟件中的多場求解器MFX進行雙向多場耦合分析，建立固體計算域和流體計算域分別使用ANSYS和CFX同時迭代求解。固體計算域內通過功率分布設定將燃料棒束熱流量傳遞給流體域的冷卻劑，冷卻劑在棒束流道內流動的過程中冷卻棒束交換熱量，并對燃料組件結構施加力的影響，組件由于流體力和熱應力產生形變和位移反過來影響流動和溫度的分布。在每個時間步長內通過交錯循環的方式實現結構場、溫度場、流場三者間的雙向數據傳遞，以此保證每個耦合的時間步都是真實的收斂解，求解流程及多場間的數據傳遞如圖7所示。

圖6 徑向截面二次流速度云圖及局部速度矢量圖Fig.6 Contour of secondary flow velocity and local velocity vector at radial-intersection

圖7 求解流程(a)及數據傳遞(b)Fig.7 Solving procedure (a) and data transfer (b)

在多場耦合的條件下對冷卻劑流動和傳熱進行數值模擬計算，并與單流場模擬結果進行對比，結果如圖8所示。由圖8可見，兩種條件下二次流速度和溫度分布隨冷卻劑流動的變化趨勢是一致的，均是在冷卻劑流經定位格架后產生更大的二次流，隨流動對格架下游產生影響并逐漸減弱，溫度從入口到出口逐漸升高。從二次流速度和溫度分布的變化規律來看，兩種條件下二次流速度和溫度均符合中心通道大于邊通道大于角通道的分布規律，考慮到管內流動中壁面對流動的客觀影響，這種分布規律符合在相同功率分布的條件下，角通道的流速較低、熱量交換效果較差的事實，需注意角通道可能出現局部熱點的情況。

圖8 單流場與多場耦合條件下二次流速度與溫度的對比Fig.8 Comparison of secondary flow velocity and temperature between single flow field and multi-field coupled cases

另一方面，單流場條件下二次流的峰值速度和平均速度均明顯高于多場耦合條件下的分布，冷卻劑的溫度波動和平均溫度均高于多場耦合條件下的分布。這是由于多場耦合時定位格架內結構與冷卻劑相互作用使結構形變，特別是攪混翼片受冷卻劑沖擊產生沿流速方向的變形使冷卻劑在流經攪混翼片時受到的阻力降低，較小的流動壓降產生較低的橫向二次流速度，致使冷卻劑與棒束間的換熱較弱，冷卻劑溫度更低。多場耦合條件更貼近于冷卻劑流動的實際情況，這一結果對比可避免單流場條件下數值模擬對流動和增強換熱過于樂觀的估計，可提升安全分析過程的可靠性。以Plane7為例，橫向截面的二次流速度分布云圖和局部子通道二次流速度矢量圖對比(圖9)同樣可見相同趨勢。

2.4 多場耦合條件下燃料棒束的結構特性分析

實際堆芯中冷卻劑與燃料棒束是相互作用的：一方面冷卻劑受燃料棒束形狀結構的影響改變了流場和溫度分布，另一方面燃料棒束受到冷卻劑傳熱和流動沖擊產生形變和振動改變了結構，是一個雙向多場耦合的過程。為分析燃料棒束的結構特性，采用多場耦合條件模擬全長燃料組件7跨中的1跨。考慮到結構和材料特性對棒束兩端施加固定約束，定位格架附加彈性約束，格架彈簧基礎剛度為2.9 N/mm3 [14]。

a——單流場條件；b——多場耦合條件圖9 Plane7的二次流速度云圖及局部速度矢量圖對比Fig.9 Comparison of secondary flow velocity contour and local velocity vector of Plane7

從整個燃料棒束結構來看，定位格架攪混翼片自身的功能特性決定了最大形變出現在攪混翼片頂端，最大變形量為2.057 μm(圖10)。燃料棒束的形變分布如圖11所示。從對安全影響更為關鍵的燃料棒束來看，最大形變出現在定位格架下游外層棒束的中部，最大變形量為0.653 μm(圖11a)。整個棒束結構的形變均為μm量級，形變對結構安全的影響很小。

進一步的模擬中將定位格架從燃料棒束模型中去除，結果顯示燃料棒束的形變小于含定位格架的模型，最大形變出現在外側棒的中部，最大值僅為0.320 μm，燃料棒的形變由中間到兩邊呈對稱分布(圖11b)。

圖10 定位格架的形變分布Fig.10 Deformation distribution of grid spacer

a——隱藏定位格架；b——無定位格架圖11 燃料棒束的形變分布Fig.11 Deformation distribution of fuel rod

圖12 不同位置燃料棒束的形變分布Fig.12 Deformation distribution of fuel rod in different positions

中心棒、中層棒和外層棒沿冷卻劑流動方向的形變分布如圖12所示。圖12中實線所示為有格架情況下的形變，從軸向沿程分布看，由于兩端的固定約束和定位格架的彈性約束使燃料棒形變均符合先增長再減小再明顯增長最后減小的過程，且定位格架下游棒束形變量明顯大于格架上游棒束形變量，這是由于定位格架的攪混作用和燃料棒功率分布特性致使處于不同位置的燃料棒形變分布呈現出不均勻的特點所致。從形變的最大值和平均值來看，外層棒的形變大于中層棒，二者均大于中心棒。結合圖9b Plane7的二次流分布可知，形變是由于流動分布不均勻導致的，這種大范圍的環繞棒束繞流在棒的兩側分布越不均勻，對棒束橫向產生的應力就越大。外層棒由于單側流動最顯著出現了所有棒最大的形變位移，而中心棒附近雖然出現了二次流速度最大值，但兩側流動分布較為均勻，棒的形變是3個位置最小的，中層棒的分布介于兩者之間。

圖12中虛線所示為無格架情況下的形變，通過對比有無定位格架燃料棒束的形變可知，定位格架對冷卻劑流動的攪混作用是明顯的，冷卻劑的橫向流速提升對燃料棒束的沖擊增大，使棒束產生更大的徑向形變。但冷卻劑對棒束的沖擊作用有限，燃料棒束最大形變小于1 μm，棒束可保持較好的結構穩定性。

3 結論

本文通過建立冷卻劑與含格架燃料棒束的幾何模型，采用數值模擬方法獲得了冷卻劑在5×5含定位格架燃料棒束通道內流動的分布，對比了單流場條件與多場耦合條件下冷卻劑的流動和溫度分布特性，并開展了冷卻劑與燃料棒束間的雙向多場耦合分析，得出如下結論。

定位格架通過擾動冷卻劑使棒束通道間形成橫向流動，并在定位格架下游形成沿對角線流動的繞流，顯著增強了冷卻劑與燃料棒間的對流換熱。多場耦合條件下二次流的峰值速度和平均速度及冷卻劑的溫度波動和平均溫度均小于單流場條件下的分布，前者更貼近于冷卻劑流動的實際情況，采用多場耦合條件分析可避免對流動和增強換熱過于樂觀的估計。在多場耦合條件下冷卻劑的橫向流動與燃料棒束的熱應力使棒束發生軸向和徑向的形變，燃料棒束的最大形變量小于1 μm，流動分布不均勻導致燃料棒形變分布呈現出不均勻的特點。有無定位格架的模擬對比表明定位格架的存在使冷卻劑的攪混流動明顯，冷卻劑對燃料棒的沖擊也更強，但總地來說兩種情況下棒束產生的形變均很小，仍能保持原本結構的穩定。