晁嫣萌，楊立新，張玉相，龐錚錚

(1.北京交通大學(xué) 機械與電子控制工程學(xué)院，北京 100044；2.中科華核電技術(shù)研究院有限公司，廣東 深圳 518026)

定位格架是反應(yīng)堆燃料組件的重要部件，影響著堆芯的熱工水力性能。開發(fā)自主知識產(chǎn)權(quán)的燃料組件需深入了解定位格架對熱工水力特性的影響，相比于周期長、費用高的試驗研究，CFD模擬已成為安全、快速的重要設(shè)計輔助手段之一，應(yīng)用CFD進行燃料組件內(nèi)定位格架對流動傳熱特性的影響分析具有重要工程價值和意義。

1997年，Lee等[1]運用非線性k-ε湍流模型在未作任何模型參數(shù)調(diào)整情況下開展了數(shù)值仿真，得出該模型很大程度上低估了湍流強度，說明湍流模型對CFD模擬結(jié)果有著重要的影響。2002年，Smith等[2]采用RNGk-ε模型對帶攪混翼格架燃料組件內(nèi)流場進行了研究，結(jié)果與實驗數(shù)據(jù)吻合較好。2003年，Yadigaroglu等[3]對棒束通道的數(shù)值計算結(jié)果表明，標準k-ε模型在窄縫空間不能很好地預(yù)測湍流流動。2012年，Liu等[4]采用Fluent軟件，利用簡化的燃料組件棒束通道模型，計算了6種不同湍流模型下Nu的沿程分布和典型位置的周向分布，其結(jié)果表明，Realizablek-ε湍流模型結(jié)合近壁面函數(shù)模擬結(jié)果與實驗結(jié)果吻合較好。Holloway等[5-6]對不同類型的格架進行了單相流動換熱實驗研究，對棒束通道內(nèi)換熱系數(shù)沿程及周向分布進行了測量，為后續(xù)開展CFD模擬的標定提供了重要實驗參考數(shù)據(jù)。

目前研究表明，采用不同湍流模型會使相同CFD模型的模擬結(jié)果不同。本文針對中國廣核集團研發(fā)的格架結(jié)構(gòu)，建立包含5×5典型格架結(jié)構(gòu)和棒束通道CFD分析模型，采用ANSYS CFX軟件，通過模擬獲得不同湍流模型下的燃料組件棒束沿程及周向換熱系數(shù)分布，與其他學(xué)者的相關(guān)實驗數(shù)據(jù)進行對比分析，通過對3個典型攪混效果評價因子的分析，探討子通道內(nèi)流動與換熱的內(nèi)在關(guān)系，同時對比不同湍流模型對結(jié)果的影響，以得到最適宜本文所研究定位格架及棒束通道內(nèi)流動傳熱特性計算的湍流模型。

1 幾何模型

定位格架幾何結(jié)構(gòu)復(fù)雜，因彈簧、剛突與燃料棒間是線接觸或小弧面接觸，使局部網(wǎng)格劃分困難，網(wǎng)格質(zhì)量差，對CFD模擬結(jié)果影響較大。本文簡化方法是在彈簧、剛突與燃料棒間留出0.1 mm的縫隙，如圖1所示，同時刪除原條帶上的一些定位孔等不影響流動特性的微小結(jié)構(gòu)。

圖1 接觸部位的簡化

2 網(wǎng)格生成

采用非結(jié)構(gòu)化四面體網(wǎng)格加拉伸網(wǎng)格的形式，應(yīng)用ICEM CFD完成模型網(wǎng)格劃分。格架部分由于結(jié)構(gòu)復(fù)雜使用非結(jié)構(gòu)化網(wǎng)格，流道部分幾何規(guī)則且長徑比較大，故采用拉伸網(wǎng)格，網(wǎng)格區(qū)域劃分如圖2所示。在燃料棒與剛突、彈簧縫隙處采用線控制和面控制進行網(wǎng)格細化，如圖3所示。在格架條帶和燃料棒束壁面設(shè)置附面層網(wǎng)格以提高網(wǎng)格質(zhì)量和效率，如圖4所示。

圖2 網(wǎng)格類型區(qū)域的劃分

圖3 格架局部表面網(wǎng)格示意圖

圖4 附面層局部網(wǎng)格細節(jié)圖

本文對格架四面體網(wǎng)格大小、拉伸網(wǎng)格數(shù)量以及附面層網(wǎng)格參數(shù)進行了敏感性分析，最終確定的計算網(wǎng)格模型參數(shù)列于表1。根據(jù)計算結(jié)果對燃料棒表面的Y+(第1層網(wǎng)格中心到壁面的無量綱距離)取平均值，如表2所列。不同湍流模型Y+值均在3～5范圍內(nèi)，說明計算模型壁面網(wǎng)格能滿足各湍流模型要求。

表1 網(wǎng)格數(shù)量統(tǒng)計

表2 不同湍流模型Y+值

3 計算模型

3.1 邊界條件

計算邊界條件與Holloway等[5-6]的實驗工況保持一致。實驗工況為常壓，冷卻劑水入口溫度為20 ℃，加熱棒表面熱流密度為10 kW/m2，Re分別為28 000和42 000，對應(yīng)的CFD計算模型具體邊界條件設(shè)置列于表3，表中未提到的邊界參數(shù)均按軟件默認設(shè)置。

表3 邊界條件設(shè)置

計算采用不可壓縮穩(wěn)態(tài)流動傳熱模型，選取6種湍流模型進行了對比分析，水工質(zhì)物性取CFX軟件IAPSW數(shù)據(jù)庫參數(shù)，計算求解采用高階精度差分格式，收斂標準為各物理量殘差小于5×10-5，同時出口流量和平均溫度不再發(fā)生變化。

3.2 湍流模型

本文進行了渦粘和雷諾應(yīng)力兩類湍流模型計算，其中包括渦粘模型中的標準k-ε、RNGk-ε和SST 3種模型以及雷諾應(yīng)力模型中LRR-IP、LRR-QI和SSG 3種模型。渦粘模型是假設(shè)雷諾應(yīng)力與平均速度梯度呈正比，雷諾應(yīng)力模型未應(yīng)用渦粘假設(shè)，而是求解流體雷諾應(yīng)力輸運方程。

標準k-ε模型目前在科學(xué)研究及工程實際中得到最為廣泛的檢驗和應(yīng)用，但其用于強旋流動和彎曲壁面流動時仍會產(chǎn)生一定的失真，RNGk-ε對湍流黏度進行了修正，SST模型考慮了湍流剪切應(yīng)力的輸運，對在逆壓力梯度下流動的分離起點和發(fā)展能預(yù)測得很準。與渦粘性模式相比，雷諾應(yīng)力模型考慮了流動中的旋轉(zhuǎn)及旋轉(zhuǎn)流動的影響，包含以下因素：流線曲率、變形率突變、二次流及浮力，更適用于應(yīng)變場復(fù)雜的各向異性流場，尤其是存在大流線曲率或渦旋的流場。3種雷諾應(yīng)力模型由不同學(xué)者發(fā)展，有各自不同的模型常數(shù)。

6種湍流模型計算時，除SST外，均采用了默認的Scalable壁面函數(shù)，SST模型使用自動壁面函數(shù)。

4 結(jié)果分析

本文以Holloway實驗結(jié)果為參考，對比分析了不同湍流模型CFD計算結(jié)果。Holloway實驗格架為不帶彈簧剛突的撕裂式攪混翼結(jié)構(gòu)，本文CFD模型格架為帶彈簧剛突的分開式攪混翼結(jié)構(gòu)，圖5為兩種格架結(jié)構(gòu)俯視圖。通過數(shù)值結(jié)果與實驗結(jié)果的對比，分析了格架不同攪混翼類型對流動的影響，同時在模型結(jié)構(gòu)完全一致的棒束通道區(qū)域進一步驗證了不同湍流模型的計算結(jié)果。下面從沿程壓降、流動攪混性能評價因子以及Nu分布對比分析不同湍流模型的計算結(jié)果。

圖5 不同攪混翼格架俯視圖

4.1 沿程壓降

圖6為兩種雷諾數(shù)工況下格架和棒束通道計算壓降與實驗結(jié)果對比。圖6a示出了格架段的壓降和棒束通道的壓降測點位置，圖6b、c分別為兩種雷諾數(shù)工況下的各湍流模型計算壓降與實驗結(jié)果比較情況，圖中兩條水平直線位置分別代表實驗的兩段壓降測量值。

由圖可知，對于格架段壓降，6種湍流模型結(jié)果均高于實驗值，雷諾應(yīng)力模型和實驗值較吻合，渦粘模型結(jié)果偏高；對于棒束區(qū)壓降，除SST模型外，其余數(shù)值預(yù)測結(jié)果均低于實驗值，雷諾應(yīng)力模型平均低20%。僅就棒束區(qū)域壓降計算結(jié)果看，標準k-ε和SST模型與實驗最接近。對比CFD模型和實驗?zāi)Ｐ透窦芙Y(jié)構(gòu)，可知格架段壓降計算值偏高顯然是由于彈簧/剛突和攪混翼結(jié)構(gòu)不同而引起的，因格架段壓降中也包含一完整棒束區(qū)壓降，考慮模擬結(jié)果在棒束區(qū)的偏低量，可推算出由于格架結(jié)構(gòu)不同所引起的壓降約占總壓降的15%。

4.2 流場攪混特性

大量文獻[7-10]報道了格架攪混翼能在下游一定距離的橫截面產(chǎn)生并保持一定強度的渦流和交叉流。存在于子通道的渦流引起的離心力能帶走燃料棒壁面形成的汽泡，存在于子通道間的交叉流可均衡子通道間的焓升。不同形式的評價因子被定義來研究攪混翼所引起的流動特性，雖然其具體形式上略有差別，但所描述的流動本質(zhì)相同，這些因子可概括為渦流攪混因子、交叉流攪混因子及湍流強度因子，這3個是影響燃料棒表面換熱特性及臨界熱流密度的重要流動參數(shù)。

渦流攪混因子Svortex可用來衡量渦流對流體的攪混作用，定義為：

其中：R為子通道中心到燃料棒表面的距離；r為通道截面上流域內(nèi)各點到中心的距離；vt為橫截面上引起渦流的橫向速度；ua為橫截面上各點的軸向速度。

圖6 沿程壓降

交叉流攪混因子Fcross用來衡量子通道間交叉攪混的強度，定義為：

其中：s為燃料棒表面間距(柵距-棒直徑)；Vcross為穿過棒間截面的交叉流速度分量；Ubulk為截面平均軸向流速。

湍流強度因子Tt為：

其中：k為湍動能；U為子通道內(nèi)軸向平均流速。

取中心棒束左上角的子通道對上述3個因子進行積分計算，圖7示出了6種湍流模型3個因子在兩種雷諾數(shù)工況下的沿程變化曲線。圖中橫坐標采用了軸向距離與燃料棒直徑比值的無量綱參數(shù)，Dh為燃料棒直徑，零點位置為格架出口位置。由圖7可知，兩種Re工況下的攪混和交叉流因子在數(shù)量級及沿程分布趨勢上均一致，湍流強度因子分布趨勢一致，在高Re工況下該值略低。

圖7a、b表明，渦流攪混因子受攪混翼影響在剛出格架時值最大，在0～3Dh區(qū)間迅速衰減，3Dh后又開始增加并在10Dh左右達到第2個峰值，這區(qū)間受攪混翼的影響顯著，10Dh～25Dh區(qū)間攪混翼的影響逐漸減弱，其值波動衰減，25Dh后該因子已降至較低水平，并繼續(xù)沿程線性下降，攪混翼的影響基本消失。各湍流模型結(jié)果趨勢一致，除SST模型值略小外，其余曲線基本重合。

圖7c、d表明，交叉流攪混因子剛出格架時較小，在攪混翼作用下迅速上升，并在3Dh處達到最大值，3Dh～10Dh間波動下降，該區(qū)間交叉攪混因子受攪混翼的影響顯著，10Dh～25Dh區(qū)間攪混翼的影響逐漸減弱，該因子呈波動衰減，25Dh后該因子也降至較低水平，并沿程線性下降。雷諾應(yīng)力3個湍流模型結(jié)果基本重合，渦粘模型中SST模型與其他結(jié)果偏差略大。結(jié)合圖7a～d，可看到在渦流攪混因子較高的位置交叉流攪混因子較低(如格架出口及10Dh處)，同時當渦流攪混因子下降時交叉流攪混因子則有升高的趨勢(10Dh～25Dh區(qū)間)，可見攪混翼引起的繞子通道中心渦流強度與子通道間交叉流強度相互耦合。

圖7e、f表明，湍流強度因子出格架后迅速下降，在3Dh～10Dh間略有波動，5Dh～25Dh區(qū)間基本維持不變，在下游后段又呈現(xiàn)上升趨勢。SST和標準k-ε模型湍流強度因子明顯大于其他湍流模型的，這也間接說明了這兩個湍流模型計算壓降較大的原因。Re大的計算工況湍流因子量級略有減小，原因是Re增大導(dǎo)致軸向速度U增大，而湍動能k并不完全與U的二次方呈正比。

圖7 不同湍流模型3種評價因子沿程變化

4.3 換熱分析

1) 沿程平均Nu

取5×5格架中間燃料棒束表面，沿流動方向截取系列截面，對每個截面棒束圓周的Nu進行算術(shù)平均，繪制沿程平均Nu沿程變化曲線，如圖8所示。圖中分別示出了兩種雷諾數(shù)工況條件下6種湍流模型計算結(jié)果與Holloway實驗結(jié)果的對比曲線。

由圖8可知，兩種雷諾數(shù)工況下，隨Re的增大棒束表面平均Nu也增大，但沿程變化趨勢相同。6種湍流模型中，SST模型的計算結(jié)果明顯高于其他模型和實驗結(jié)果，這是由于SST模型采用的自動壁面函數(shù)對壁面網(wǎng)格更敏感導(dǎo)致的。在棒束通道區(qū)，渦粘類模型的Nu均高于雷諾應(yīng)力模型的結(jié)果，在高雷諾數(shù)情況下該趨勢更明顯。雷諾應(yīng)力類的3種湍流模型Nu基本重合，且在受格架影響較弱的棒束通道內(nèi)其值與實驗結(jié)果較為接近。

由圖8還可知，6種湍流模型計算得到的平均Nu均是在格架出口位置最大，然后迅速減小，在3Dh～10Dh區(qū)間出現(xiàn)劇烈波動，在10Dh～25Dh區(qū)間呈現(xiàn)波動下降趨勢，25Dh后Nu趨于穩(wěn)定值。沿程Nu變化區(qū)間的特征位置與上節(jié)中分析的3個流動因子變化區(qū)間位置是相對應(yīng)的，說明攪混翼引起的攪混流動直接影響著棒束表面的換熱特性。Holloway實驗中Nu曲線在離開格架位置后的0～5Dh區(qū)間迅速降低，然后有小幅上升，并在25Dh后Nu分別穩(wěn)定在210和280。對比可知，在25Dh位置前，計算結(jié)果與實驗值有較大差異，這是由于本文CFD模型采用了分開式攪混翼格架，攪混翼長度、寬度等形狀參數(shù)與實驗所用的撕裂式攪混翼格架不同(圖5)，計算結(jié)果表明分開式攪混翼會引起更強的渦流和交叉流，從而提高了格架下游子通道的換熱能力。經(jīng)過25Dh后，流體進入棒束通道，此時格架和攪混翼對流動的影響已大幅減弱，而棒束通道實驗段和CFD模型完全相同，所以這一區(qū)間計算與實驗Nu的對比能用來更直接地驗證模型的準確性。綜合分析不同Re下計算工況的結(jié)果，雷諾應(yīng)力類的3種湍流模型計算的Nu與實驗值吻合得更好。

2) 棒束周向Nu

Holloway實驗給出了2.2Dh、6.5Dh和35.7Dh位置的Nu周向分布，取計算模型的中間燃料棒束對應(yīng)位置截取3個截面，圖9左側(cè)示出了這3個位置棒束表面不同湍流模型計算Nu的分布，右側(cè)示出了對應(yīng)位置截面上的橫向流速度分布矢量圖。

由圖9左側(cè)的Nu曲線對比可知，除SST模型外，其他湍流模型預(yù)測值均與實驗值較接近。2.2Dh和6.5Dh位置的Nu周向分布實驗值與計算值在分布的波形上相位略有區(qū)別，這是由于計算模型彈簧剛突的存在以及攪混翼結(jié)構(gòu)的不同所引起的。Nu分布明顯受到攪混翼所引起的橫向速度分布的影響，對照右側(cè)的橫向速度分布矢量圖可發(fā)現(xiàn)：Nu較大的位置橫向速度較大(如棒束圓周105°、345°等位置)，而Nu處于波谷位置的橫向速度較小或局部有渦流出現(xiàn)(如棒束圓周45°、145°等位置)。在35.7Dh位置5種湍流模型的預(yù)測值與實驗值已基本重合，攪混翼引起的橫向流動影響也基本消失。

圖8 不同湍流模型Nu沿程分布與實驗值

圖9 不同湍流模型Nu周向分布

5 結(jié)論

1) 渦粘模型中標準k-ε和SST模型對壓降預(yù)測較好，而對換熱系數(shù)的預(yù)測結(jié)果與其他湍流模型和實驗值偏差較大，其中SST模型對Nu的過高預(yù)測與CFX中該模型采用的自動壁面函數(shù)有關(guān)；

2) 格架及攪混翼結(jié)構(gòu)給流動子通道內(nèi)帶來強烈的渦流和交叉流，雷諾應(yīng)力類湍流模型更適用于這種各向異性的流場模擬，3種雷諾應(yīng)力模型結(jié)果基本一致，Nu與實驗結(jié)果吻合較好，是后續(xù)進行這類格架詳細性能分析及設(shè)計優(yōu)化的首選湍流模型；

3) 攪混翼引起的沿程渦流攪混因子和交叉流攪混因子具有強烈的耦合關(guān)系，并共同影響著沿程平均Nu分布，這種影響趨勢在格架下游10Dh后逐漸減弱，在25Dh后基本消失。

參考文獻：

[1] LEE K B, JANG H G. A numerical prediction on the turbulent flow in closely spaced bare rod arrays by a nonlineark-εmodel[J]. Nuclear Engineering and Design, 1997, 172: 351-357.

[2] SMITH L D, CONNER M E, LIU B, et al. Benchmarking computational fluid dynamics for application to PWR fuel[C]∥Proceedings of the 10th International Conference on Nuclear Engineering. USA: ASME, 2002: 823-830.

[3] YADIGAROGLU G, ANDERANI M, DREIER J, et al. Trends and needs in experimentation and numerical simulation for LWR safety[J]. Nuclear Engineering and Design, 2003, 221: 205-223.

[4] LIU C C, FERNG Y M, SHIH C K. CFD evaluation of turbulence models for flow simulation of the fuel rod bundle with a spacer assembly[J]. Applied Thermal Engineering, 2012(40): 389-396.

[5] HOLLOWAY M V, McCLUSKY H L, BEASLEY D E, et al. The effect of support grid features on local, single-phase heat transfer measurements in rod bundles[J]. ASME Journal of Heat Transfer, 2004, 126: 43-53.

[6] HOLLOWAY M V, CONOVER T A, McCLUSKY H L, et al. The effect of support grid design on azimuthal variation in heat transfer coefficient for rod bundles[J]. ASME Journal of Heat Transfer, 2005, 127: 598-605.

[7] CUI X Z, KIM K Y. Three-dimensional analysis of turbulent heat transfer and flow through mixing vane in a subchannel of nuclear reactor[J]. Journal of Nuclear Science and Technology, 2003, 40(10): 719-724.

[8] LEE C M, CHOI Y D. Comparison of thermo-hydraulic performances of large scale vortex flow (LSVF) and small scale vortex flow (SSVF) mixing vanes in 17×17 nuclear rod bundle[J]. Nuclear Engineering and Design, 2007, 237(24): 2 322-2 331.

[9] IN W K, CHUN T H, OH D S, et al. CFD analysis of turbulent flows in rod bundles for nuclear fuel spacer design[C]∥Transactions of the 15th International Conference on Structural Mechanics in Reactor Technology. Seoul: SMiRT, 1999: 365-372.

[10] IN W K, CHUN T H, SHIN C H, et al. Numerical computation of heat transfer enhancement of a PWR rod bundle with mixing vane spacers[J]. Nuclear Technology, 2008, 161(1): 69-79.