棒束子通道湍流交混特性的數(shù)值模擬研究

張 銳，馬在勇，陸定晟，何清澈，張盧騰，孫 皖，潘良明

(重慶大學 低品位能源利用技術及系統(tǒng)教育部重點實驗室，重慶 400044)

反應堆堆芯棒束燃料組件內(nèi)部流動傳熱特性影響著堆芯相關熱工參數(shù)，子通道間湍流交混效應使得通道內(nèi)臨界熱流密度(CHF)提高。研究反應堆堆芯棒束子通道間湍流交混效應有助于提高反應堆的安全性和經(jīng)濟性[1]。

子通道間的湍流交混效應是一種由多因素(通道間壓差、幾何結構和機械部件等)導致的一種綜合性交混效應[2]。國內(nèi)外學者針對子通道間湍流交混效應展開了大量研究。王正杰等[3]通過無因次分析方法得出了描述子通道間單相湍流交混率的無因次方程式，并通過示蹤法實驗結果計算了近似方形排列的子通道單相湍流交混率以及無因次方程式的系數(shù)和指數(shù)。康森厚等[4]推導出了對任何形式排列棒束子通道皆適用的單相湍流交混無因次通用方程。曹念等[5]采用帶格架的5×5全長棒束開展子通道單相交混系數(shù)實驗研究，使用子通道分析程序配合實驗數(shù)據(jù)間接獲得了單相交混系數(shù)。謝峰等[6]對19根棒束繞絲結構燃料組件進行實驗，并結合子通道程序計算得到湍流交混系數(shù)。張琦等[7]對5×5螺旋十字型棒束(HCF)組件進行熱工水力實驗，基于能量平衡法對HCF組件的交混特性進行了分析，發(fā)現(xiàn)HCF組件的等效交混系數(shù)不隨雷諾數(shù)的增加而明顯變化，其均值為0.019。Kawahara等[8]提出了一種用于預測核燃料棒束子通道間氣液兩相湍流交混率的彈狀-攪混流模型，并用湍流交混率數(shù)據(jù)證實了該模型的有效性。Sharabi等[9]研究了幾種湍流模型在非圓截面通道內(nèi)的傳熱預測能力，并通過實驗比較了湍流模型用于預測觀測現(xiàn)象的性能。Wang等[10]提出了計算流體力學方法在反應堆工程中面臨的主要挑戰(zhàn)和發(fā)展方向，有利于推動其進一步應用。Ju等[11]提出了一種結合高通濾波大渦模擬(LES)模型的高可擴展高性能譜元方法，用于模擬單相工況下通過平行壁面通道和矩形排列棒束通道中的湍流流動，并通過實驗和理論驗證了湍流交混系數(shù)依賴于雷諾數(shù)的模擬結果，證明了譜元法結合LES模型的可行性。于意奇等[12]運用CFD方法對稠密柵元內(nèi)典型子通道湍流流動進行RANS和URANS模擬分析，分析比較不同湍流計算模型特性。李小暢等[13]對無格架和帶交混翼格架的5×5棒束進行數(shù)值模擬分析，通過計算在不同工況下的努賽爾數(shù)及流體湍流強度，對不同湍流模型的影響及近壁面網(wǎng)格處理的適用范圍進行了比較分析。王雅峰[14]基于CFD數(shù)值計算理論對不同節(jié)徑比和不同雷諾數(shù)的流體湍流工況狀態(tài)進行了數(shù)值模擬，從而得到橫向交混率的計算公式。對于不同子通道結構，裸棒束與帶混合葉片的棒束也有相關研究，In等[15]通過一系列CFD模擬，分析了帶攪混翼格架的全加熱棒束的強化傳熱，與已有的實驗數(shù)據(jù)進行了比較，并進行了相關分析。Ikeno等[16]對裸棒束內(nèi)湍流流動進行了大渦模擬，分析了方形排列棒束內(nèi)的湍流流動，并表明CFD的重要作用之一是利用湍流模型模擬流動中的物理現(xiàn)象，以及預測平均流動。Ju等[17]采用LES和非定常RANS(URANS)方法的譜元法研究了反應堆堆芯中帶有定位格架和攪混翼格架的棒束中的湍流，通過將湍動能等值線與橫向速度矢量相結合，研究了攪混翼格架形成的渦流行為，揭示了渦流行為。

上述研究表明，對于子通道間湍流交混效應的研究是非常有必要的。目前對子通道間湍流交混效應的實驗研究主要是通過化學示蹤劑方法來完成，這種方法只能得到局部點或面的參數(shù)，無法得到全局參數(shù)，同時進行實驗需要消耗巨大的人力、物力和財力。為此，本文擬采用CFD方法對子通道間的湍流交混效應進行數(shù)值模擬研究。

1 數(shù)值計算模型

1.1 實驗本體建模

本文以方形棒束排列的相鄰兩個子通道作為研究對象，采用DesignModeler軟件建立節(jié)徑比(P/D)為1.326的幾何模型。燃料棒外直徑(D)為9.5 mm，節(jié)距(P)為12.6 mm，軸向流動方向長度為800 mm。流動橫截面如圖1a所示，幾何模型在橫截面上網(wǎng)格劃分情況如圖1b所示。

圖1 流動橫截面(a)和網(wǎng)格劃分(b)

本文選擇單相液態(tài)水作為工質(zhì)，設定兩子通道入口水溫分別為300 K、320 K，子通道溫度平均值對應物性參數(shù)，即密度ρ=992.5 kg/m3，動力黏度μ=7.02×10-4Pa·s，導熱率λ=0.6 W/(m·K)，求解單相工況下對應的湍流交混系數(shù)。

1.2 網(wǎng)格敏感性驗證

采用ICEM軟件進行幾何模型的網(wǎng)格劃分。保持z方向長度4 mm不變，流動橫截面網(wǎng)格采用相同劃分方法，通過不斷加密z方向節(jié)點數(shù)量得到4套網(wǎng)格。表1列出4套網(wǎng)格劃分情況，采用M1、M2、M3、M4 4套網(wǎng)格分別進行完全相同工況下的數(shù)值計算，通過分析計算結果實現(xiàn)網(wǎng)格敏感性驗證。

表1 網(wǎng)格數(shù)量對比

采用入口速度u=1.5 m/s，其他條件不變，計算對應工況流體速度分布及壁面切應力分布。為保證流動充分發(fā)展，取z=0.75 m處流動橫截面上直線L1處的速度分布以及中心位置燃料棒與流體交界的圓弧R1上的應力分布，L1與R1位置如圖2所示。圖3示出M1、M2、M3、M4對應L1位置速度分布與R1位置壁面切應力分布。

圖2 R1與L1位置示意圖

圖3 L1位置速度分布(a)與R1位置壁面切應力分布(b)

由圖3可知，4種網(wǎng)格計算所得L1位置處的速度分布變化趨勢一致且數(shù)值上差別不大，而R1位置處的壁面切應力分布變化趨勢基本一致，但數(shù)值上有所差別。M1與M4最大相差近5%，M2與M4最大相差近4%，而M3與M4之間數(shù)值上非常接近，相差程度都在1%以內(nèi)。因此通過R1位置處的壁面切應力分布比較，可得出當網(wǎng)格數(shù)量在154萬以上時計算結果已與網(wǎng)格數(shù)量無關，故選擇M3網(wǎng)格進行數(shù)值模擬計算，不影響計算結果準確性且能節(jié)約計算資源，降低計算成本。

1.3 湍流模型選擇

選擇M3網(wǎng)格進行數(shù)值模擬計算，采用SSTk-ω模型、標準k-ω模型、標準k-ε模型、RNGk-ε模型和RSM模型分別在入口速度u=0.5、1、1.5、2 m/s對應工況進行數(shù)值模擬計算。為保證湍流流動充分發(fā)展，選取計算區(qū)間為z=0.25 m流動橫截面與z=0.75 m流動橫截面之間的流體壓降Δp。圖4示出不同湍流模型的摩擦阻力系數(shù)f與Blasius公式計算的理論值。

由圖4可知，5種不同湍流模型對應的f與理論值在雷諾數(shù)較大時都具有較好吻合性。其中SSTk-ω模型和標準k-ω模型的f與理論值在雷諾數(shù)較小時相差較小，而標準k-ε模型、RNGk-ε模型以及RSM模型的f在雷諾數(shù)較小時則相差很大。SSTk-ω模型和標準k-ω模型之間，前者相比于后者在低雷諾數(shù)下的f更加接近于理論值，且其隨雷諾數(shù)變化趨勢更加接近于理論值變化趨勢。在5種湍流模型中，SSTk-ω模型對于子通道湍流流動數(shù)值模擬最為合理精確，因此選擇SSTk-ω模型作為本文數(shù)值模擬的湍流模型。

圖4 不同湍流模型摩擦阻力系數(shù)與理論值

2 數(shù)值模擬結果分析

2.1 間隙湍流熱流法

保持兩個相鄰子通道具有相同的幾何形狀和相等的入口質(zhì)量流速。兩個相鄰子通道內(nèi)流體保持一定入口溫差，其中通道i中流體溫度較高，通道j中流體溫度較低。保證沒有其他的熱流輸入，兩個相鄰子通道間的傳熱只能夠通過熱傳導和湍流交混兩個過程來實現(xiàn)。

由湍流交混和通過間隙傳導而產(chǎn)生的熱量(在穩(wěn)態(tài)流動下)等于：

Qtotal=Qin-Qout=q″totalSijΔx

(1)

通過間隙的熱流為：

(2)

其中：ρ為子通道中流體密度，kg/m3；U為子通道中流速，m/s；Aflow為子通道面積，m2；cp為比定壓熱容，J/(kg·K)；Sij為子通道之間的間隙，m；Δx為作用的長度，m；Qtotal為由湍流交混和通過間隙傳導而產(chǎn)生的熱量，kJ；Qin為流體進口熱量，kJ；Qout為流體出口熱量，kJ；q″total為通過間隙的總熱流，kW；Tin為入口流體溫度，K；Tout為出口流體溫度，K；i、j為子通道。

湍流交混系數(shù)β定義為湍流所導致的熱流和軸向熱流的比值，即：

(3)

因為水的湍流導熱率遠大于其分子導熱率，式(3)可改寫為：

(4)

從理論上來講，式(4)是將湍流和分子擴散的交混系數(shù)一并計算出來的，同時在這里應強調(diào)基于子通道熱平衡的熱流將包括對流和擴散傳輸?shù)呢暙I。真正的差別只在于局部(近間隙)的熱平衡。由此可得：

(5)

設置入口速度u=1.5 m/s，保持其他條件不變進行數(shù)值模擬，計算軸向z=0.05～0.75 m每間隔0.1 m處兩子通道橫截面溫度，再采用間隙湍流熱流法計算湍流交混系數(shù)軸向變化，如圖5所示。

圖5 間隙湍流熱流法計算的湍流交混系數(shù)軸向變化

由于存在入口段效應，在軸向z=0.35 m前所算得的湍流交混系數(shù)較大且沿著軸向方向迅速減小。在軸向z=0.35 m后湍流交混系數(shù)趨于穩(wěn)定近似等于2×10-3。將對應雷諾數(shù)Re以及動力黏度μ代入Galbraith計算關系式：

(6)

其中，W′ij為湍流交混率。

由單相湍流交混系數(shù)計算式可得：

(7)

對應單相湍流交混系數(shù)β計算結果為2.79×10-3，與間隙湍流熱流法所得結果保持數(shù)量級10-3相吻合，且數(shù)值差異不大。因此，采用間隙湍流熱流法能較準確地計算湍流交混系數(shù)，且沿流體軸向流動方向能反映局部湍流交混系數(shù)的變化。

2.2 類比濃度計算法

通過在子通道冷卻劑中加入一定濃度示蹤劑，再測定示蹤劑濃度在子通道中軸向變化情況，從而推導計算湍流交混系數(shù)。通常將濃度沿軸向變化數(shù)據(jù)擬合如式(8)來確定湍流交混率W′ij：

(8)

其中：Gi和Ai分別為子通道i的質(zhì)量流速和橫截面積；Ci(z)為子通道i軸向z處的示蹤劑濃度；Δz為兩個軸向位置之間的距離。

由于傳質(zhì)方程與傳熱方程相類比，濃度沿軸向的變化來推算湍流交混率，可類比為溫度沿軸向的變化來推算湍流交混率，從而通過軸向溫度變化得到湍流交混率：

(9)

其中，Ti(z)為子通道i在軸向z處的流體溫度。

本文以2005～2017年長江中游城市群31個地級市為研究對象，研究產(chǎn)業(yè)集聚、技術創(chuàng)新對經(jīng)濟增長的影響。專利申請量來源于各城市每年的國民經(jīng)濟和社會發(fā)展統(tǒng)計公報，消費價格指數(shù)來源于湖北、湖南、江西3省的統(tǒng)計年鑒，2005～2016年其余數(shù)據(jù)來源于2006～2017《中國城市統(tǒng)計年鑒》，2017年數(shù)據(jù)來源于各市2017年國民經(jīng)濟和社會發(fā)展統(tǒng)計公報。

設置入口速度u=1.5 m/s，保持其他條件不變，計算軸向z=0.05～0.75 m每間隔0.1 m處兩子通道橫截面溫度，再采用類比濃度計算法計算湍流交混系數(shù)的軸向變化，如圖6所示。

圖6 類比濃度計算法得到的湍流交混系數(shù)軸向變化

類比濃度計算法所計算結果同樣存在入口段效應，在軸向z=0.35 m前所算得湍流交混系數(shù)較大且沿著軸向方向迅速減小。在z=0.35 m后湍流交混系數(shù)趨于穩(wěn)定近似等于2×10-3。湍流交混系數(shù)數(shù)量級在10-3，同樣與Galbraith計算關系式計算結果相吻合，且在數(shù)值上差異不大。

2.3 兩種計算方法比較

對比間隙湍流熱流法與類比濃度計算法所得到的湍流交混系數(shù)軸向變化，可看出兩者幾乎一致，相對偏差不超過1.5%，并且兩者具有相同的軸向變化趨勢，能很好地體現(xiàn)出流體流動入口段效應，對于湍流交混系數(shù)計算具有較高精度。

2.4 單相湍流交混系數(shù)擬合關系式

本文在子通道水力直徑、燃料棒直徑以及間隙寬度都一定的情況下，對相鄰兩子通道湍流交混進行數(shù)值模擬，計算得到不同雷諾數(shù)下單相湍流交混系數(shù)。其他條件一定下，分析雷諾數(shù)對單相湍流交混系數(shù)的影響，并擬合雷諾數(shù)與單相湍流交混系數(shù)的關系式。雷諾數(shù)范圍為4 500～34 000。

假設湍流交混系數(shù)與雷諾數(shù)關系式為冪指數(shù)關系：

β=bRea

(10)

其中，a、b為待定常數(shù)。

lgβ=lgb+algRe

(11)

對其進行擬合可得如圖7所示的結果。

由圖7可知，類比濃度計算法與間隙湍流熱流法擬合直線斜率分別為a1=0.092 7、a2=0.091 5，與y軸截距分別為b1=-3.108 2、b2=-3.113 2，對應關系式中a1=0.092 7、a2=0.091 5，b1=0.000 78、b2=0.000 77，故單相湍流交混系數(shù)擬合關系式為：

圖7 冪指數(shù)對數(shù)化擬合圖

β1=0.000 78Re0.092 7

(12)

β2=0.000 77Re0.091 5

(13)

假設湍流交混系數(shù)與雷諾數(shù)關系式為多項式關系：

β=aRe2+bRe+c

(14)

對其進行擬合，結果如圖8所示。

圖8 多項式擬合曲線

由圖8可知，擬合曲線a1=-1.377×10-12、a2=-1.309×10-12，b1=5.461×10-8、b2=5.209×10-8，c1=0.001 43、c2=0.001 41，故單相湍流交混相系數(shù)擬合關系式為：

β1=-1.377×10-12Re2+

5.461×10-8Re+0.001 43

(15)

β2=-1.309×10-12Re2+

5.209×10-8Re+0.001 41

(16)

上述冪指數(shù)、多項式擬合關系式都是通過最小二乘法進行曲線擬合，各類關系式對應相關系數(shù)R2列于表2。

表2 各類擬合曲線相關系數(shù)

在雷諾數(shù)小于15 000時，湍流交混系數(shù)隨雷諾數(shù)的增大而增大，但增長幅度較小。而在雷諾數(shù)大于15 000時，隨著雷諾數(shù)的變化湍流交混系數(shù)幾乎不變。因此在經(jīng)過充分發(fā)展后，可認為湍流交混系數(shù)β為一常數(shù)，擬合湍流交混系數(shù)常數(shù)如圖9所示。

圖9 湍流交混系數(shù)常數(shù)擬合圖

由圖9可知，擬合湍流交混系數(shù)β常數(shù)分別為β1=0.001 9、β2=0.001 85，在雷諾數(shù)小于15 000時β與擬合值相差最大近10%，而在雷諾數(shù)大于15 000時β與擬合值相差最大近3%，類比濃度計算法與間隙湍流熱流法所對應擬合常數(shù)相差僅為2.6%。

綜合上述幾種擬合方式，多項式擬合雖相關系數(shù)R2較高，但由于多項式擬合關系式結構復雜且一次項與二次項系數(shù)太小，因而其擬合關系式不太適用較廣范圍。雖然冪指數(shù)擬合關系式相關系數(shù)R2不是最高，但其相關性能滿足要求，且擬合單相湍流交混系數(shù)關系式結構簡單以及系數(shù)大小合適，較為簡單清晰便于廣泛適用。

使用間隙湍流熱流法計算湍流交混系數(shù)時忽略了分子導熱，進而導致使用類比濃度計算法所得的湍流交混系數(shù)略大于間隙湍流熱流法的結果。同時隨著雷諾數(shù)增大，一方面湍流交混強度增強，另一方面間隙處渦旋的影響范圍較小，因此由湍流引起的能量轉(zhuǎn)移會減小，兩方面的綜合作用會使得湍流交混系數(shù)處于變化波動較小的狀態(tài)。

但總的來說，采用類比濃度計算法與間隙湍流熱流法能準確計算子通道湍流交混系數(shù)，兩者計算湍流交混系數(shù)的軸向分布沒有顯著差異。

3 結論

本文使用CFD數(shù)值模擬的方法對單相工況下反應堆堆芯棒束子通道間湍流交混特性進行研究，研究發(fā)現(xiàn)：采用類比濃度計算法與間隙湍流熱流法能夠準確計算子通道湍流交混系數(shù)，兩者計算湍流交混系數(shù)的軸向分布沒有顯著差異；同時發(fā)現(xiàn)湍流交混系數(shù)隨雷諾數(shù)的增大而增大，但雷諾數(shù)足夠大時保持基本不變。基于CFD計算結果擬合得到單相湍流交混系數(shù)的計算關系式。本文對于子通道間湍流交混效應的研究有一定的幫助。