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

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

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

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

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

1 幾何模型

定位格架幾何結(jié)構(gòu)復(fù)雜，因彈簧、剛突與燃料棒間是線接觸或小弧面接觸，使局部網(wǎng)格劃分困難，網(wǎng)格質(zhì)量差，對CFD模擬結(jié)果影響較大。本文簡化方法是在彈簧、剛突與燃料棒間留出0.1 mm的縫隙，如圖1所示，同時(shí)刪除原條帶上的一些定位孔等不影響流動特性的微小結(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所示。在燃料棒與剛突、彈簧縫隙處采用線控制和面控制進(jìn)行網(wǎng)格細(xì)化，如圖3所示。在格架條帶和燃料棒束壁面設(shè)置附面層網(wǎng)格以提高網(wǎng)格質(zhì)量和效率，如圖4所示。

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

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

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

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

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

表2 不同湍流模型Y+值

3 計(jì)算模型

3.1 邊界條件

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

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

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

3.2 湍流模型

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

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

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

4 結(jié)果分析

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

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

4.1 沿程壓降

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

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

4.2 流場攪混特性

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

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

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

圖6 沿程壓降

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

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

湍流強(qiáng)度因子Tt為：

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

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

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

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

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

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

4.3 換熱分析

1) 沿程平均Nu

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

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

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

2) 棒束周向Nu

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

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

圖8 不同湍流模型Nu沿程分布與實(shí)驗(yàn)值

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

5 結(jié)論

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

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

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

參考文獻(xiàn)：

[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.