稠密柵元不同子通道內(nèi)湍流流動(dòng)的RANS和URANS模擬

于意奇，顧漢洋，楊燕華，程 旭

（上海交通大學(xué) 核科學(xué)與工程系，上海 200240）

稠密柵元的燃料組件廣泛應(yīng)用于先進(jìn)的反應(yīng)堆設(shè)計(jì)。因稠密的柵元結(jié)構(gòu)將減少水鈾比從而硬化了中子能譜，這將提高轉(zhuǎn)換比并增加燃料利用率［1－3］。

最早的棒束內(nèi)的湍流流動(dòng)實(shí)驗(yàn)研究開始于20世紀(jì)60年代。近年來，隨著實(shí)驗(yàn)測量手段的發(fā)展，越來越多的實(shí)驗(yàn)研究開始關(guān)注于棒束內(nèi)的湍流流動(dòng)［4－9］。

實(shí)驗(yàn)研究表明:棒束內(nèi)的湍流流動(dòng)和圓管內(nèi)的完全不同，在棒束間隙區(qū)有很強(qiáng)的交混。這一現(xiàn)象曾被歸結(jié)于二次流，不過，最近的實(shí)驗(yàn)研究表明這并不是影響棒束間強(qiáng)交混的主要因素。棒束間隙區(qū)的準(zhǔn)周期流動(dòng)振動(dòng)才是影響子通道間能量質(zhì)量交換的最主要因素。這種流動(dòng)振動(dòng)和子通道的幾何構(gòu)造與雷諾數(shù)密切相關(guān)，當(dāng)雷諾數(shù)低于一定的臨界值時(shí)，這種振動(dòng)便消失［10］。實(shí)驗(yàn)發(fā)現(xiàn)，壁面通道中，在1.015≤W／D≤1.250（W 為棒最遠(yuǎn)端距壁面的距離，D為棒徑）的范圍內(nèi)，流動(dòng)振動(dòng)的頻率隨間隙的減小而減小［11］。盡管對這一流動(dòng)振動(dòng)現(xiàn)象還未完全了解，但學(xué)術(shù)界普遍認(rèn)為這是由于流動(dòng)的不穩(wěn)定造成的。

如今，RANS、LES、URANS和 DNS方法均被應(yīng)用于模擬棒束內(nèi)的湍流流動(dòng)［12－16］。研究表明:采用各向同性的湍流模型的穩(wěn)態(tài)模擬（RANS）無法再現(xiàn)棒束流動(dòng)的強(qiáng)各向異性的特點(diǎn)。盡管采用各向異性的湍流模型的穩(wěn)態(tài)模擬可模擬出湍流驅(qū)動(dòng)的二次流現(xiàn)象，但在P／D（P為兩棒之間的中心距離）較小的稠密柵元子通道內(nèi)，二次流不是影響流動(dòng)的主要因素。所以，在這些子通道中，模擬的結(jié)果往往與實(shí)驗(yàn)有一定的差距。盡管數(shù)值模擬粗略高估了振動(dòng)的波長，URANS方法對稠密柵元棒束間隙的流動(dòng)振動(dòng)有較好的模擬結(jié)果。

在這些方法中，DNS方法是最理想的一種，因?yàn)閷τ谡w的流動(dòng)振動(dòng)現(xiàn)象還不是十分了解。LES方法幾乎可得到和DNS相同的結(jié)果。但由于LES和DNS的計(jì)算代價(jià)太大，URANS是目前比較實(shí)際的一種數(shù)值模擬方法。

1 數(shù)值方式

在本研究中，將選用不同P／D和通道形狀的子通道實(shí)驗(yàn)作為參考。詳細(xì)的實(shí)驗(yàn)參數(shù)和數(shù)值模擬信息列于表1［4，17－18］。實(shí)驗(yàn)截面圖示于圖1。

表1 本工作的實(shí)驗(yàn)和數(shù)值參數(shù)Table 1 Experimental and numerical parameters

圖1 實(shí)驗(yàn)截面圖Fig.1 Set－up cross－section

1.1 RANS模擬

由于子通道的對稱性，選用1／6和1／4子通道作為RANS模擬的計(jì)算域，如圖2所示。棒、通道和表面采用非滑移的邊界條件，進(jìn)口采用均流的邊界條件，出口設(shè)置定壓的邊界條件，其余為對稱的邊界條件。通過網(wǎng)格敏感性分析，分別得到25萬的中心子通道網(wǎng)格和50萬的壁面子通道網(wǎng)格。采用渦粘性和雷諾應(yīng)力模型來評估它們對于稠密柵元內(nèi)湍流模擬的適用性。渦粘性模型包括k－ε模型和SST模型，雷諾應(yīng)力模型包括ORS、SSG和BSL模型。所有模型均采用壁面函數(shù)模擬近壁面的流動(dòng)。所有近壁面網(wǎng)格保持y＋＝1。

圖2 RANS的計(jì)算域Fig.2 Computational domain for RANS

1.2 URANS模擬

URANS計(jì)算域包括兩個(gè)由間隙連接的子通道（圖3）。整個(gè)計(jì)算包括3組周期性邊界條件和非滑移的壁面邊界條件。進(jìn)出口的周期性邊界條件固定質(zhì)量流速，計(jì)算長度為0.6m（4倍于實(shí)驗(yàn)測量的平均波長λ）。

圖3 URANS的計(jì)算域Fig.3 Computational domain for URANS

網(wǎng)格總數(shù)為60萬，時(shí)間步長滿足:Δt?1／f，f為最小振動(dòng)頻率，Δt選為10－4s。幾乎所有應(yīng)用于RANS模擬的湍流模型均被應(yīng)用于URANS計(jì)算。另外，專為瞬態(tài)模擬設(shè)計(jì)的SAS湍流模型也將應(yīng)用于URANS計(jì)算。

2 結(jié)果和討論

本工作對所有工況均進(jìn)行了RANS模擬，對P／D＝1.06的三角形排列的中心通道和4棒束平行排列的子通道進(jìn)行了URANS計(jì)算，并將主流速度、壁面剪應(yīng)力、湍動(dòng)能和壁面溫度等參數(shù)與文獻(xiàn)［4］的實(shí)驗(yàn)結(jié)果進(jìn)行比較。

2.1 主流速度

對于P／D為＝1.17的三角形排列子通道，RANS模擬計(jì)算結(jié)果示于圖4。圖中，y為徑向上與壁面距離，v為主流速度。棒束內(nèi)的二次流雖尺度很小，但它將影響子通道內(nèi)的湍流分布。雷諾應(yīng)力模型將湍流的各向異性考慮在內(nèi)，所以，雷諾應(yīng)力湍流模型（BSL，ORS，SSG）較渦粘性湍流模型（k－ε，SST）得到的速度分布更接近實(shí)驗(yàn)測量結(jié)果。由圖3b可見，SSG湍流模型的模擬計(jì)算結(jié)果與實(shí)驗(yàn)測量結(jié)果吻合良好。

P／D＝1.06的三角形排列稠密柵子通道內(nèi)的速度分布示于圖5。圖中，Wb為主流速度，ymax為徑向上速度最大值與壁面的距離。在棒束間隙區(qū)（φ＝30°），RANS和 URANS的計(jì)算結(jié)果有很大差異，其原因在于棒束間隙區(qū)存在強(qiáng)烈的速度振動(dòng)。在URANS的計(jì)算結(jié)果中，各向異性的SSG湍流模型給出了對實(shí)驗(yàn)更好的預(yù)測。4棒平行排列稠密柵元子通道的速度分布示于圖6。由圖6同樣發(fā)現(xiàn)，在棒束間隙區(qū)，URANS的結(jié)果較RANS有很大改進(jìn)。

圖4 P／D＝1.17的三角形排列子通道主流速度分布Fig.4 Velocity distributions in triangular array with P／D＝1.17

圖5 P／D＝1.06的三角形排列稠密柵子通道內(nèi)速度分布Fig.5 Velocity distributions in triangular array with P／D＝1.06

圖6 P／D＝1.12、W／D＝1.06的4棒束平行排列子通道速度分布Fig.6 Velocity distributions in 4rod parallel array with P／D＝1.12，W／D＝1.06

2.2 壁面剪應(yīng)力

圖7示出P／D＝1.17的三角形排列子通道內(nèi)的壁面剪應(yīng)力tau的分布。圖中，taum為壁面平均剪應(yīng)力。圖7a顯示，雷諾應(yīng)力湍流模型（SSG、BSL、ORS）的模擬結(jié)果更為平坦，較渦粘性湍流模型（k－ε、SST）也更接近實(shí)驗(yàn)結(jié)果。這是由于雷諾應(yīng)力湍流模型對二次流的模擬均化了壁面剪應(yīng)力的分布。圖7b示出了雷諾數(shù)更高的情況，模擬結(jié)果同樣證實(shí)了上述結(jié)論。

圖8示出P／D＝1.06的三角形排列子通道內(nèi)的壁面剪應(yīng)力分布。由圖8可見，在URANS計(jì)算方式下，即便采用各向同性的k－ε湍流模型，得到的結(jié)果也比在RANS計(jì)算方式下采用SSG湍流模型得到的結(jié)果更接近實(shí)驗(yàn)測量結(jié)果。這也證明了在P／D較小的稠密柵元子通道中，二次流不再是主導(dǎo)的影響因素。

圖7 P／D＝1.17的三角形排列子通道內(nèi)的壁面剪應(yīng)力分布Fig.7 Wall shear stress distributions in triangular array with P／D＝1.17

圖8 P／D＝1.06的三角形排列子通道內(nèi)壁面剪應(yīng)力分布Fig.8 Wall shear stress distributions in triangular array with P／D＝1.06

同樣，在URANS計(jì)算方式下，雷諾應(yīng)力湍流模型（SSG、ORS）的模擬結(jié)果較渦粘性湍流模型（k－ε、SST）更接近實(shí)驗(yàn)測量結(jié)果。

URANS同樣得到更為均化的壁面通道壁面剪應(yīng)力的分布（圖9），并與實(shí)驗(yàn)吻和良好。在RANS計(jì)算中，只有各向異性的湍流模型可模擬出二次流。而在URANS計(jì)算中，流動(dòng)的對稱性被打破，在短時(shí)間平均下無法觀察到二次流，但在足夠長的時(shí)間下采用雷諾平均，即使采用各向同性的湍流模型也能重現(xiàn)二次流。

2.3 湍動(dòng)能

稠密柵間隙區(qū)的流動(dòng)振動(dòng)是一種相干結(jié)構(gòu)，在比較湍動(dòng)能時(shí)，有必要把這種流動(dòng)振動(dòng)貢獻(xiàn)考慮在內(nèi)。通常，相對湍動(dòng)能的定義如下:

式中:u、v、w分別為軸向速度波動(dòng)、徑向速度波動(dòng)和周向速度波動(dòng)；uτ為壁面剪應(yīng)力；k＋nc為湍動(dòng)能的非相干部分。由速度振動(dòng)引起的湍動(dòng)能稱為湍動(dòng)能的相干部分k＋c，其定義如下:

總湍動(dòng)能為兩部分之和:

圖10示出三角形排列和4棒束平行排列的子通道在不同周向角度上的湍動(dòng)能分布。由圖10可見，URANS的計(jì)算精確性高于RANS。專門為非穩(wěn)態(tài)計(jì)算設(shè)計(jì)的SAS湍流模型較其他湍流模型有著更高的模擬精度。

圖9 P／D＝1.12、W／D＝1.06的4棒束平行排列子通道壁面剪應(yīng)力分布Fig.9 Wall shear stress distributions in 4rod parallel array with P／D＝1.12，W／D＝1.06

圖10 稠密柵元子通道湍動(dòng)能分布Fig.10 Turbulent intensity distributions in tight lattice

2.4 壁面溫度

圖11示出P／D＝1.06的三角形排列子通道壁面溫度分布，溫度通過平均溫度來進(jìn)行無量綱化處理。由圖11可見，URANS的模擬結(jié)果更為均化并與實(shí)驗(yàn)測量結(jié)果吻合良好。因此，稠密柵元內(nèi)的流動(dòng)振動(dòng)均化了壁面溫度，對于反應(yīng)堆堆芯的工業(yè)設(shè)計(jì)是有利的。

圖11 P／D＝1.06的三角形排列子通道壁面溫度分布Fig.11 Wall temperature distributions in triangular array with P／D＝1.06

2.5 流動(dòng)振動(dòng)

稠密柵元子通道內(nèi)的流動(dòng)振動(dòng)緩慢發(fā)展，最后發(fā)展為穩(wěn)定的波動(dòng)（波幅達(dá)到常值或準(zhǔn)常值）。圖12示出三角形排列子通道窄縫區(qū)的橫向速度隨時(shí)間的變化和三維振動(dòng)快照。振動(dòng)幅度為2m／s，與文獻(xiàn)［4］的實(shí)驗(yàn)測量結(jié)果相近。然而，計(jì)算所得振動(dòng)波長比實(shí)驗(yàn)測量結(jié)果大。

總體而言，URANS仍無法模擬稠密柵元內(nèi)相干結(jié)構(gòu)的所有特性，但對于正確預(yù)測主流速度、壁面剪應(yīng)力和壁面溫度等工業(yè)關(guān)心的統(tǒng)計(jì)變量十分實(shí)用。

3 結(jié)論

本工作對典型的子通道（中心通道和壁面通道）進(jìn)行URANS和RANS計(jì)算。計(jì)算結(jié)果表明，在較寬的稠密柵元通道（P／D＝1.17）內(nèi)，采用SSG湍流模型的RANS計(jì)算能很好地再現(xiàn)實(shí)驗(yàn)結(jié)果。對于P／D較?。≒／D＜1.1）的稠密柵元子通道，窄縫區(qū)內(nèi)大尺度的流動(dòng)振動(dòng)將顯著影響子通道的湍流流動(dòng)。URANS模擬雖無法再現(xiàn)稠密柵元子通道內(nèi)流動(dòng)振動(dòng)的所有動(dòng)力特性，但可很好地預(yù)測子通道內(nèi)的主流速度、壁面剪應(yīng)力、湍動(dòng)能分布。推薦采用SAS和SSG湍流模型進(jìn)行URANS計(jì)算。

圖12 P／D＝1.06的三角形排列子通道窄縫區(qū)橫向速度波動(dòng)（a）和三維流動(dòng)振動(dòng)快照（b）Fig.12 Flow pulsation（a）and dimension snapshot of oscillation（b）in narrow gap in triangular array subchannel with P／D＝1.06

［1］OLDEKOP W，BERGER H D，ZEGGEL W.General features of advanced pressurized water reactors with improved fuel utilization［J］.Nuclear Technology，1982，59:212－227.

［2］UCHIKAWA S，OKUBO T，KUGO T，et al.Investigations on innovative water reactor for flexible fuel cycle（LFWR）［C］∥Proceedings of GLOBAL.Tsukuba，Japan:［s.n.］，2005.

［3］CHENG X，LIU X J，YANG Y H.A mixed core for supercritical water－cooled reactors［J］.Nuclear Engineering and Technology，2008，40（1）:1－10.

［4］KRAUSS T，MEYER T.Experimental investigation of turbulent transport of momentum and energy in a heated rod bundle［J］.Nuclear Engineering and Design，1998，180:185－206.

［5］TRUPP A C，AZAD R S.The structure of turbulent flow in triangular array rod bundles［J］.Nuclear Engineering and Design，1975，32:47－84.

［6］TRIPPE G，WEINBERG D.Non－isotropic eddy viscosities in turbulent flow through rod bundles［M］∥Turbulent Forced Convection in Channels and Bundles:Vol.1.KAKAC S，SPALDING D B.Washington:Hemisphere Publishing Corporation，1979:505.

［7］SEALE W J.Turbulent diffusion of heat between connected flow passages［J］.Nuclear Engineering and Design，1979，54:183－195.

［8］REHME K.Simple method of predicting friction factors of turbulent flow in non－circular channels［J］.International Journal of Heat and Mass Transfer，1973，16:933－950.

［9］REHME K.The structure of turbulent flow through rod bundles［J］.Nuclear Engineering and Design，1987，99:141－154.

［10］MEYER L，REHME K.Large－scale turbulence phenomena in compound rectangular channels［J］.Experimental Thermal and Fluid Science，1994，8:286－304.

［11］BARATTO F，BAILEY S C C，TAVOULARIS S.Measurements of frequencies and spatial correlations of coherent structures in rod bundle flows［J］.Nuclear Engineering and Design，2006，236:1 830－1 837.

［12］BAGLIETTO E，NINOKATA H.Turbulence models evaluation for heat transfer simulation in tight lattice fuel bundles［C］∥Proceedings of the 10th International Topical Meeting on Nuclear Reactor Thermal Hydraulics （NURETH－10）.Seoul，Korea:［s.n.］，2003.

［13］CHANG D，TAVOULARIS S.Unsteady numerical simulations of turbulence and coherent structures in axial flow near a narrow gap［J］.ASME Journal of Fluids Engineering，2005，127（3）:458－466.

［14］CHANG D，TAVOULARIS S.Numerical simulation of turbulent flow in a 37－rod bundle［J］.Nuclear Engineering and Design，2007，237:575－590.

［15］MERZARI E，NINOKATA H，BAGLIETTO E.Large eddy simulation of the vortex street between rectangular channels［C］∥NTHAS 5.Jeju Island，Korea:［s.n.］，2006.

［16］NINOKATA H，MERZARI E，KHAKIM A.Analysis of low Reynolds number turbulent flow phenomena in nuclear fuel pin subassemblies of tight lattice configuration［J］.Nuclear Engineering and Design，2009，239（5）:855－866.

［17］KRAUSS T.Experimentelle untersuchung des turbulenten w?rme－und impulstransports in einem beheizten stabbündel ［R］.Germany:Karlsruhe University，1996.

［18］MANTLIK F，HEINA J，CHERVENKA J.Results of local measurements of hydraulic characteristic in triangular pin bundle，UJV－3778－R［R］.Rzez，Czech Republic:［s.n.］，1976.