花瓣形燃料元件棒束通道內過冷流動沸騰特性數值研究

杜利鵬，蔣澤平，崔 軍，張文超，侯延棟，蔡偉華,*

(1.東北電力大學 熱流科學與核工程實驗室，吉林 吉林 132012；2.東北電力大學 能源與動力工程學院，吉林 吉林 132012；3.中廣核研究院有限公司，廣東 深圳 518031)

過冷流動沸騰是指在主流溫度低于飽和溫度情況下，發生在加熱壁面附近的流動沸騰換熱現象，該現象的發生極大增強了壁面的換熱系數，因此，其廣泛存在于核反應堆、電子散熱系統等換熱設備中。對于核反應堆，為提高核反應堆功率和獲得更高的傳熱效率，高性能的燃料元件開發已成為重要手段之一。花瓣形燃料元件通過其自身螺旋結構使冷卻劑攪混，換熱性能增強，且無需定位格架，能夠較大地提升核反應堆熱工水力性能；與傳統壓水堆設計相比，在相同堆芯功率下，最大臨界熱流比提高58%[1]。因此，花瓣形燃料元件棒束通道內流動換熱特性的研究，得到了國內外學者的廣泛關注。

花瓣形燃料元件的螺旋結構會對流動產生重要影響。Shirvan和Kazimi[2]通過數值模擬方法對4×4花瓣形燃料元件棒束通道內的單相流動換熱特性進行研究，分析了花瓣形燃料元件的結構及螺旋節距對通道內橫向流動的影響。Xiao等[3]通過數值模擬，分析了通道內的橫向流動分布，得到了橫向攪混模型。蔡偉華等[4]研究發現，燃料元件螺旋產生的橫向流動主要集中在燃料元件的內凹弧區域。同時，花瓣形燃料元件具有的螺旋結構，使其區別于傳統的棒束通道內的換熱特性。Fang等[5]對花瓣形燃料元件和非花瓣形燃料元件棒束通道內的高溫氫氣流動換熱進行數值分析；從熱工和水力綜合性能來看，花瓣形燃料元件相較于十字非螺旋形和圓棒形元件分別提高了28%和6.1%。蔡偉華等[4]發現通道內產生的橫向流動使得加熱壁面溫度呈不均勻分布，內凹弧的壁面溫度明顯大于外凸弧。Zeng等[6]發現入口流速和花瓣燃料元件的截面幾何結構顯著影響冷卻劑換熱特性。Shirvan[7]和Cong等[8]對花瓣形燃料元件的偏離核態沸騰(DNB)性能進行了數值研究，均發現蒸汽主要集中在燃料元件的內凹區域。此外，Cong等[8]發現，發生沸騰危機的位置不是相鄰燃料棒間距最小時的截面位置。

在核反應堆高壓高溫運行工況下，通過實驗來研究燃料元件棒束通道內冷卻劑的過冷流動沸騰是極為困難的。隨著計算流體力學(CFD)方法的高速發展，尤其是歐拉(Eulerian-Eulerian)兩流體模型建立之后，進行燃料元件棒束通道內的兩相流動數值模擬成為了可能。Krepper等[9]通過實驗數據評估了歐拉兩流體和RPI(Rensselaer Polytechnic Institute)壁面沸騰模型對過冷沸騰現象的預測潛力，發現大部分模擬工況下的計算值都與實驗值符合較好，高壓工況下預測偏差較大。Colombo等[10]通過20組實驗數據評估了歐拉兩流體模型對過冷沸騰預測能力，認為封閉模型限制了數值模擬的普適性。Pan等[11]驗證了相間作用力模型對過冷沸騰的影響。Gu等[12]對適用于高壓下的RPI壁面沸騰模型的封閉模型進行評估，并給出相對最優的模型組合。對于復雜通道內的過冷沸騰，Yang等13]和張蕊等[14]分別都使用RPI壁面沸騰模型，對具有定位格架和攪混翼的燃料棒束通道內的過冷沸騰現象進行了數值研究。

目前，針對花瓣形燃料元件棒束通道內冷卻劑的流動與換熱數值研究主要集中在單相流動與換熱、橫向攪混及偏離核態沸騰等方面。但花瓣形燃料元件棒束通道內的過冷流動沸騰特性的研究對核反應堆的安全設計具有重要的指導意義，而該方面的研究仍較為缺乏。因此，本文采用CFD方法，基于歐拉兩流體模型和RPI壁面沸騰模型，對2×2花瓣形燃料元件棒束通道的過冷流動沸騰開展相關數值研究。

1 兩相沸騰模型

過冷沸騰模型是基于歐拉兩流體模型——將汽液流體均看成連續介質，汽液之間可以進行質量、能量、動量交換，交換量的大小在數學模型上體現在氣液各相的質量、動量、能量方程的源項上。RPI壁面沸騰模型是在歐拉兩流體模型的基礎上，充分考慮了汽泡產生、生長及脫離現象對流動與換熱的影響，建立的壁面沸騰模型，能較為準確地反映出過冷流動沸騰特性。

1.1 歐拉兩流體模型

質量守恒方程：

(1)

式中，k、α、ρ、u、Δmk分別為氣相v或液相l、相體積份額、密度(kg/m3)、速度(m/s)、兩相之間質量傳遞(kg/(m3·s))。

動量守恒方程：

(2)

式中，p、g、Flv分別為壓力(Pa)、重力加速度(m/s2)、相間作用力(N)。

汽泡脫離加熱壁面后，進入主流區運動，會受到多種相間作用力的影響。相間作用力Flv包括曳力FD、升力FL、壁面潤滑力FWL、湍流耗散力FTD。

Flv=FD+FL+FWL+FTD

(3)

汽泡和液體間的曳力FD為：

(4)

式中：Ai為相界面面積，m2；μl為液相黏度，kg/(m·s)；Re為汽泡雷諾數；dB為汽泡平均直徑，m；CD為曳力系數，采用Ishii&Zuber[15]模型進行計算。

當連續相不均勻或旋轉時，汽泡會受到垂直于速度方向的升力FL：

(5)

式中，CL為升力系數，采用Tomiyama[16]模型計算。

汽泡從壁面產生，壁面潤滑力會使汽泡從壁面附近進入主流區域，壁面潤滑力FWL模型如下：

(6)

式中：n為壁面方向矢量；CWL為壁面潤滑力系數，采用Antal[17]給定的關系式計算。

由于主流區液相的湍流脈動產生湍流耗散力，使得汽泡從高濃度區向低濃度區運動。該力采用Lopez-de-Bertodano[18]模型進行計算。

(7)

式中：CTD為湍流耗散力系數，取值0.1；k為湍動能，m2/s2。

同時采用Sato[19]模型修正汽泡黏度，以此方法來考慮汽泡對液相湍流的影響。

能量守恒方程：

(8)

式中：h為比焓，J/kg；q為兩相中的能量傳遞，W，由牛頓冷卻公式計算：

q=hlAi(Tsat-Tl)

(9)

式中：hl為對流換熱系數，W/(m2·K)，由Ranz Marshall[20]關系式計算；Tsat、Tl分別為飽和溫度和液相溫度，K。

1.2 RPI壁面沸騰模型

壁面沸騰模型通過對壁面熱流密度分配來體現過冷沸騰中的傳熱現象，進而與汽泡的產生、生長、脫離等行為特性形成相互影響。由Kural和Podowski[21]提出的RPI壁面沸騰模型將壁面熱流q″w分成3部分：單相熱流密度q″c、淬火熱流密度q″q及蒸發熱流密度q″e:

q″w=q″c+q″q+q″e

(10)

q″c=hc(Tw-Tl)(1-Ab)

(11)

(12)

(13)

式中，hc、tw、Ab、f、λl、ddep、Nw、hlv分別為單相對流換熱系數(W/(m2·K))、汽泡等待時間(s)、汽泡影響區域面積(m2)、汽泡脫離頻率(1/s)、液相導熱系數(W/(m·K))、汽泡脫離直徑(m)、汽泡核化密度(m-2)、汽化潛熱(J/kg)。

為使上述方程封閉，還需要參數——汽泡脫離直徑ddep、汽泡脫離頻率f、汽泡核化密度Nw及汽泡影響區域面積Ab的模型。上述參數計算關系式如下：

(14)

(15)

Nw=2101.085(Tw-Tsat)1.085

(16)

(17)

式中：a、b、φ可參考文獻[22]；Ja為過冷度雅克布數。

2 模型驗證

采用Bartolemei[23]的均勻加熱圓管高壓實驗及相關實驗數據進行模型驗證。幾何模型如圖1a所示，實驗工況如下：熱流密度為570 kW/m2、壓力為4.5 MPa、質量流速為900 kg/(m2·s)。模擬采用Standardk-ω結合增強壁面函數。通過對比實驗與模擬得到的軸向壁面溫度、流體溫度及空泡份額分布情況，分析二者誤差大小，驗證模型可靠性。圖1b為模擬獲得的結果與實驗數據對比，并通過式(18)計算相應的壁面溫度、流體溫度及空泡份額的均方根誤差，計算結果分別為1.61 K、2.07 K、0.044。可知，模擬結果和實驗數據吻合良好。

(18)

a——Bartolemei實驗幾何模型；b——計算值與實驗值對比圖1 模型驗證Fig.1 Model validation

3 幾何模型及網格劃分

3.1 幾何模型

為研究花瓣形燃料元件對棒束通道內過冷沸騰的流動和傳熱特性影響，本文選取2×2花瓣形燃料元件棒束通道為研究對象，其幾何結構形式如圖2～4所示。具體結構參數列于表1。實際應用中，在燃料元件旋轉了1/4的螺距處，燃料棒之間通過4個接觸點支撐；而網格劃分中，支撐點處網格難以劃分，所以在燃料元件間增加一定的距離。蔡偉華等[4]分析了不同距離對冷卻劑傳熱的影響，發現距離d為0.5 mm時，子通道內流體溫度相對偏差小于0.1%，所以本文采用上述參數建立模型。

圖2 花瓣形燃料元件Fig.2 Petal-shape fuel element

圖3 流體域Fig.3 Fluid domain

圖4 燃料元件排列方式Fig.4 Fuel element arrangement

3.2 網格劃分及邊界條件

使用ICEM軟件對2×2花瓣形燃料元件棒束通道進行結構化網格劃分。將整個流體域劃分成含燃料元件的內流體域和方形壁面的外流體域，中間通過interface面進行兩部分的插值計算。具體網格如圖5所示。入口和出口分別設置為速度進口和壓力出口，燃料元件壁面定義為無滑移壁面，采用均勻熱流條件，方形外壁面設置為周期性邊界條件。表2中列出了數值模擬工況具體參數(具體值參考壓水堆可能出現的運行工況)。

表1 幾何參數Table 1 Geometric parameter

圖5 計算域網格Fig.5 Calculation domain grid

表2 模擬工況Table 2 Simulation of working condition

3.3 湍流模型及網格無關性驗證

湍流模型在流動換熱模擬中起著至關重要的作用。壓水堆燃料組件通道內，由于組件產生旋流和近壁湍流各向異性導致通道內形成次級渦流，優先考慮選擇SSTk-ω、RNGk-ε、RSM湍流模型。如Liu等[24]對典型壓水堆燃料組件通道內的流動進行數值模擬，評估了各種湍流模型，發現SSTk-ω湍流模型對近壁面的處理性能更優。對于花瓣形燃料元件棒束通道，因其螺旋結構產生了復雜流動，蔡偉華等[4]和Fang等[5]對花瓣形燃料組件通道內單相流動評估了不同湍流模型的影響，根據相關實驗數據表明，SSTk-ω模型能獲得更好的結果。所以本文湍流模型選擇SSTk-ω模型。

根據上述模型設置，對研究對象進行網格無關性驗證。Fluent對SSTk-ω模型引入增強壁面函數，使得該模型對近壁面網格精度要求降低，允許1

圖6 網格無關性驗證Fig.6 Grid independent validation

4 結果與討論

4.1 流動特性分析

圖7為軸向不同截面的液相流動速度分布。由圖7可看出，冷卻劑以1.4 m/s的均勻速度進入通道內后，受到花瓣形燃料元件的螺旋導向作用，流場迅速發生變化，中心區域的速度最大達到2 m/s；燃料元件的內凹弧壁面附近，冷卻劑的速度低于截面平均速度，最低值為1.05 m/s；外凸弧處因其周向角度變化更大，使得流動阻力較小，近壁面附近的流體速度比內凹弧處的流速大，二者最大相差0.6 m/s。在內凹弧區域內，由于流體圍繞燃料元件流動，使得流場的分布特點也逐漸隨順時針方向旋轉，呈現不均勻性分布。花瓣形燃料元件和繞絲型燃料元件具有一定相似性，根據Song等[26]對繞絲型燃料元件的研究可知，在花瓣形燃料元件凸弧與主流相遇的一側稱為迎風側，對側為背風側。由于流動阻力，使得迎風側的壓力變得高于背風側的壓力，如圖8所示。結合圖7可看出，在外凸弧的兩側壁面附近，流速呈現非對稱分布，迎風側附近速度與主流速度更接近，背風側附近速度存在一定的梯度。

圖7 工況4軸向液相流動速度分布Fig.7 Axial liquid phase flow velocity distribution of case 4

紅色高，藍色低圖8 截面局部壓力分布Fig.8 Cross-sectional local pressure distribution

圖9 沿直線Line1和Line2液相流動速度Fig.9 Liquid phase flow velocity along Line1 and Line2

為量化分析凸弧兩側的流場不均勻性，同時鑒于燃料元件的對稱性，在內凹弧區域的迎風側(Line1)和背風側(Line2)取兩條直線，如圖7所示。圖9為直線Line1和Line2上液相流動速度變化情況。可看出，在近壁面附近，Line1處的液相速度梯度變化大于Line2處，且在一定距離內，Line2處的速度低于Line1處，此結果勢必使得Line1處的流動換熱強于Line2處。圖10為不同入口流速下，Line2處的液相流動速度隨徑向距離的變化趨勢。由圖10可知，入口流動速度從1.4 m/s增大到2.5 m/s過程中，直線上的局部流動速度達到截面平均速度的徑向距離逐漸增大，由此可知，入口流速增大，將增加內凹區域的流場不均勻性程度。

圖10 不同入口流速下Line2位置處液相流動速度Fig.10 Liquid phase flow velocity at Line2 position at different inlet flow velocities

圖11為工況3中不同軸向位置處，液相橫向流動速度分布。從圖11可見，橫向流動主要集中在花瓣形燃料元件的內凹區域，橫向流動方向與燃料元件扭轉方向一致，迎風側橫向速度大于背風側速度，其原因為在迎風側，主流流體與迎風側壁面相遇產生逆流，增大橫向流動速度，使得內凹弧區域橫向速度分布不均勻，最大值為0.045 m/s。從燃料元件相對入口位置旋轉了90°的位置處(z/H=0.25)起，通道的中間位置逐漸形成一個相對速度較小、與橫向流動方向相反的渦旋，使燃料元件內凹處與中間區域的流體攪混，起到強化換熱的效果。

圖11 工況3中不同軸向位置的液相橫向流動速度分布Fig.11 Liquid phase transverse flow velocity distribution at different axial positions in case 3

為體現花瓣形燃料元件產生的橫向流動能力，定義平均橫向流動強度F，定義如下：

(19)

式中:Vx、Vy、Vz分別為橫向速度分量和軸向速度，m/s；A、Ai分別為截面面積和網格面積，m2。

圖12為不同入口流速下橫向流動強度沿軸向變化曲線。花瓣形燃料元件由相對角度0°增到45°的過程是棒間隙增大的過程，之后相對角度由45°增到90°的過程是棒間隙減小的過程，因此燃料元件每1/4的螺距變化會形成1個變化周期。從圖12可看出，流體在經過充分發展后，橫向流動強度沿軸向呈波動變化，且其變化趨勢與花瓣形燃料元件棒束的間距周期變化一致，在間距增大時橫向流動強度減小，反之增大；由不同入口流速的模擬數據分析，入口流速從1.4 m/s增到2.0 m/s，橫向流動強度僅小幅增加，而在速度為2.0 m/s和2.5 m/s下，橫向強度變化基本一致，表明入口流動速度對橫向流動強度的影響基本可以忽略。

圖12 不同入口速度下橫向流動強度Fig.12 Transverse flow strength at different inlet flow velocities

圖13為不同熱流密度下橫向流動強度沿軸向變化趨勢。熱流密度為450 kW/m2時，通道剛發生過冷沸騰，隨熱流密度的增加，過冷沸騰程度增強，3個不同工況下的液相橫向流動強度變化趨于一致，表明氣相的產生和過冷沸騰的劇烈程度并不會對橫向流動產生影響，同時可注意到，氣相橫向流動強度大于液相橫向流動強度，這是由于氣相從壁面進入主流區時，除了受到花瓣形燃料元件的作用，還受到湍流耗散力和壁面潤滑力的影響。對于熱流密度600 kW/m2工況中的氣相橫向流動強度大于800 kW/m2工況，由式(4)、(5)可知，曳力與汽-液界面面積呈正比，升力大小與空泡份額呈正比，因此，熱流密度增大時，空泡份額和汽泡尺寸都在增加，曳力和升力同時增大，當汽液之間的曳力與升力的作用大于湍流耗散力等橫向力的作用時，氣相的橫向流動被削弱。

圖13 不同熱流密度下氣液兩相橫向流動強度Fig.13 Transverse flow intensity of gas-liquid two-phase at different heat fluxes

4.2 空泡份額分布特性

圖14為工況4中軸向不同位置處的空泡份額分布。由圖14可知，在z/H=0.375截面處，過冷沸騰最先在內凹弧區域內發生。這是由于在內凹弧區域處，流體的換熱面積更大，產生的熱量最多，且內凹區域處的流動速度更小，因此，壁面附近的流體更易達到飽和溫度，進而發生過冷沸騰。z/H=0.5位置后，發生過冷沸騰的區域增加，流體在花瓣形燃料元件的導向作用下，使得氣相擴散到不同子通道，且受中間通道內的旋流影響，在通道內形成相對均勻的分布；至通道出口位置處，內凹區域空泡份額最大為0.11，而外凸弧區域空泡份額為0.04，花瓣形燃料元件的內凹弧處的空泡份額仍大于其他位置。

圖14 工況4中不同軸向位置空泡份額分布Fig.14 Distribution of void fraction at different axial positions in case 4

圖15示出了燃料元件不同軸向位置處空泡份額沿周向的詳細分布。處于不同截面位置時，圖15a對燃料元件周向方向作出定義，其角度和圖15b、c相對應。由圖15b、c可知，z/H=0.5截面內凹區域內，空泡分額呈非對稱分布，燃料元件的背風側方向(15°、110°、195°、285°)附近的空泡份額大于其在外凸弧對應的迎風側處(70°、160°、250°、340°)的空泡份額，且在燃料元件相對z/H=0.5位置扭轉45°(z/H=0.875)后，上述現象依然存在，結合圖11的z/H=0.5截面橫向流速分布可知，在花瓣形燃料元件內凹區域中，橫向流速大的位置，其空泡份額小，反之亦然。造成這種現象的原因是較大的橫向流速，增強了冷熱流體攪混能力，導致過冷沸騰減弱，空泡份額相對較小。通過不同工況中空泡份額分布可知，在一定熱流密度下，對于壁面附近的過冷沸騰程度，內凹弧區域始終大于外凸弧區域。因此，對于壓水堆工況下，更應關注花瓣形燃料元件內凹弧區域的熱工水力特性。

a——周向角度定義；b——橫面z/H=0.5空泡份額；c——截面z/H=0.875空泡份額圖15 空泡份額沿周向的詳細分布Fig.15 Detailed distribution of void fraction along circumferential direction

4.3 換熱特性研究

圖16為工況4的不同軸向位置處，花瓣形燃料元件的壁面溫度沿周向分布情況。由圖16可知，燃料元件內凹弧區域的壁面溫度分布均勻，沿主流方向，在該區域內的過冷沸騰程度逐漸增加，壁面溫度整體呈下降趨勢；在燃料元件的外凸弧區域，存在劇烈的壁面溫度變化，結合圖7和圖14可知，z/H=0.25位置處，壁面與流體處于單相對流換熱階段，外凸弧壁面附近液體流動速度相對較大，換熱效果更好，因此，壁面溫度相較于內凹弧壁面顯著更低。在z/H=0.25截面后，外凸弧壁面溫度快速上升，在接近出口位置時，外凸弧壁面溫度大于內凹弧壁溫，外凸弧兩側的壁面溫度出現峰值。通過圖17中z/H=0.75和z/H=1位置處的單相熱流密度分布可看出，內凹弧和外凸弧交接處的單相熱流密度也存在峰值，分析其原因，主要是由于隨主流溫度的升高，單相對流換熱能力逐漸減弱，雖然該位置處存在單相對流換熱的峰值，但相變所需的汽化潛熱不足以抵消單相對流換熱的減弱，導致壁面溫度較高。

圖16 工況4中壁面溫度周向分布Fig.16 Circumferential distribution of wall surface temperature in case 4

圖17 單相熱流密度周向分布Fig.17 Circumferential distribution of single-phase heat flux

圖18為工況6單相熱流密度、淬火熱流密度、蒸發熱流密度沿軸向變化。其中蒸發熱流密度為在氣泡生長過程中，由于從液體到蒸汽的相變而發生的熱傳遞，淬火熱流密度是對氣泡脫離處的壁面與冷流體接觸時的熱傳導的量化值。在入口段，由于流體過冷度高，主要發生單相對流換熱，隨壁面溫度的增加，壁面上成核點增加，汽泡脫離壁面，蒸發熱流密度和淬火熱流密度都獲得增加，相應的單相熱流密度下降，但可看到其并未呈現線性的快速下降趨勢，還出現小幅上升，其原因是花瓣形燃料元件的螺旋結構加大了流體的攪混作用，使得單相熱流密度增強。同時可注意到，蒸發熱流密度快速增長，而淬火熱流密度維持在較低增長。根據Hsu[27]提出的核化理論，要使汽泡生長，汽泡頂部的液體溫度要達到所需的過熱度，因此對于附著在加熱壁上的小汽泡，需要一定的過熱層，使汽泡不斷變大后脫離。結合圖11可知，主要發生過冷沸騰的內凹區域內存在持續的橫向流動，使得加熱壁面上過熱層較薄，當汽泡的頂部到達過冷區域時，汽泡開始坍塌。大量汽泡生長，使得蒸發熱流密度在出口處達到最大，而汽泡脫離頻率低，淬火熱流增長慢。

圖18 工況6中熱流密度分配Fig.18 Heat flux partitioning in case 6

5 結論

采用RPI壁面沸騰模型對2×2花瓣形燃料元件棒束通道內過冷流動沸騰現象進行數值研究，分析了通道內的流速、空泡份額、壁面溫度、熱流分配等熱工水力參數的變化規律及影響因素，得出如下結論。

1) 棒束通道內，流體流動速度呈不均勻性分布，外凸弧處的流速大于內凹弧處的流速；隨入口流速增加，內凹弧區域內的流速不均勻性程度隨之增強。

2) 橫向流動強度隨相鄰燃料元件間距變化，而產生周期性波動；入口流速和熱流密度對截面平均橫向流動強度影響基本可忽略。

3) 在內凹弧區域過冷沸騰較為劇烈，空泡份額較大；橫向流動影響導致內凹區域空泡份額徑向上呈不對稱分布。

4) 內凹弧和外凸弧處，壁面溫度在周向上存在顯著差異，主要是由于流場不均勻性和換熱方式的不同；受橫向流動影響，蒸發熱流密度沿軸向逐漸升高，而淬火熱流變化不大。