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

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

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

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

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

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

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

1 數(shù)值方式

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

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

圖1 實驗截面圖Fig.1 Set－up cross－section

1.1 RANS模擬

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

圖2 RANS的計算域Fig.2 Computational domain for RANS

1.2 URANS模擬

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

圖3 URANS的計算域Fig.3 Computational domain for URANS

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

2 結(jié)果和討論

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

2.1 主流速度

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

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

圖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）也更接近實驗結(jié)果。這是由于雷諾應(yīng)力湍流模型對二次流的模擬均化了壁面剪應(yīng)力的分布。圖7b示出了雷諾數(shù)更高的情況，模擬結(jié)果同樣證實了上述結(jié)論。

圖8示出P／D＝1.06的三角形排列子通道內(nèi)的壁面剪應(yīng)力分布。由圖8可見，在URANS計算方式下，即便采用各向同性的k－ε湍流模型，得到的結(jié)果也比在RANS計算方式下采用SSG湍流模型得到的結(jié)果更接近實驗測量結(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計算方式下，雷諾應(yīng)力湍流模型（SSG、ORS）的模擬結(jié)果較渦粘性湍流模型（k－ε、SST）更接近實驗測量結(jié)果。

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

2.3 湍動能

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

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

總湍動能為兩部分之和:

圖10示出三角形排列和4棒束平行排列的子通道在不同周向角度上的湍動能分布。由圖10可見，URANS的計算精確性高于RANS。專門為非穩(wěn)態(tài)計算設(shè)計的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 稠密柵元子通道湍動能分布Fig.10 Turbulent intensity distributions in tight lattice

2.4 壁面溫度

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

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

2.5 流動振動

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

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

3 結(jié)論

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

圖12 P／D＝1.06的三角形排列子通道窄縫區(qū)橫向速度波動（a）和三維流動振動快照（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.