孫 暢，焦守華，柴 翔，楊 云

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

液態(tài)鉛鉍合金具有良好的中子學(xué)性能和較高的中子經(jīng)濟(jì)性，使反應(yīng)堆可在常壓下運(yùn)行，同時(shí)實(shí)現(xiàn)高功率密度。這不僅大幅降低了由高壓控制引起的失水事故的發(fā)生概率，同時(shí)也大幅降低了一回路設(shè)備對(duì)承壓性能的要求，是我國(guó)未來(lái)鉛基反應(yīng)堆的重要材料之一[1]。鉛鉍冷卻快堆燃料組件為六邊形密集棒束結(jié)構(gòu)，其采用了具有強(qiáng)化換熱和減小機(jī)械震蕩效果[2]的金屬繞絲進(jìn)行固定。

由于液態(tài)鉛鉍合金的高密度、強(qiáng)腐蝕及實(shí)驗(yàn)溫度要求較高等因素，開(kāi)展液態(tài)鉛鉍合金在燃料組件內(nèi)的熱工水力實(shí)驗(yàn)難度較大。此外由于液態(tài)鉛鉍合金的普朗特?cái)?shù)Pr很低(約0.025)，若使用相似準(zhǔn)則法在常規(guī)介質(zhì)中進(jìn)行?；瘜?shí)驗(yàn)也存在一定的挑戰(zhàn)。綜合實(shí)驗(yàn)技術(shù)、設(shè)備及經(jīng)費(fèi)等原因，目前國(guó)外較大規(guī)模地開(kāi)展液態(tài)鉛鉍合金在燃料組件內(nèi)的熱工水力實(shí)驗(yàn)的科研機(jī)構(gòu)僅有德國(guó)的KIT和意大利的ENEA實(shí)驗(yàn)室。德國(guó)KIT的KALLA實(shí)驗(yàn)室[3]開(kāi)展了液態(tài)鉛鉍合金在燃料棒束間的熱工水力特性研究實(shí)驗(yàn)，以研究冷卻劑的流動(dòng)阻力特性和流動(dòng)傳熱特性。Litfin等[4]開(kāi)展了液態(tài)鉛鉍合金和水冷卻劑在19棒束結(jié)構(gòu)內(nèi)的流動(dòng)阻力特性實(shí)驗(yàn)，并采用歐拉相似準(zhǔn)則進(jìn)行了關(guān)聯(lián)比較。Pacio等[5]開(kāi)展了使用格架定位的19棒束燃料組件內(nèi)的液態(tài)鉛鉍合金流動(dòng)傳熱行為實(shí)驗(yàn)，通過(guò)實(shí)驗(yàn)測(cè)得的壓降和流動(dòng)速度估算鉛鉍合金在組件內(nèi)的流動(dòng)特性，最終得出“液態(tài)鉛鉍合金與水等常規(guī)流體介質(zhì)具有相似的流動(dòng)特性”這一結(jié)論。意大利的Brasimone研究中心[6-7]對(duì)液態(tài)重金屬作為反應(yīng)堆冷卻劑的可能性展開(kāi)了研究，針對(duì)液態(tài)重金屬堆型相關(guān)設(shè)備材料耐高溫、抗腐蝕性測(cè)量，堆芯組件內(nèi)的熱工水力等開(kāi)展了一系列實(shí)驗(yàn)。

針對(duì)Pacio等[8-9]設(shè)計(jì)的帶繞絲的19棒束鉛鉍組件，本文使用商用計(jì)算流體力學(xué)(CFD)軟件STAR-CCM+對(duì)其建立全尺寸模型，劃分網(wǎng)格后選取數(shù)值模型進(jìn)行模擬計(jì)算，將該組件中子通道內(nèi)的流動(dòng)傳熱與實(shí)驗(yàn)值進(jìn)行對(duì)比驗(yàn)證，之后對(duì)組件的進(jìn)口質(zhì)量流量和總熱功率進(jìn)行敏感性分析。

1 正常工況模擬

本文以Pacio等[8-9]設(shè)計(jì)的組件為研究對(duì)象，組件由截面為正六邊形的外套管與其中呈正三角形排列的19根被金屬繞絲纏繞的加熱棒組成。燃料組件幾何參數(shù)列于表1。

表1 燃料組件幾何參數(shù)Table 1 Geometric parameter of fuel assembly

1.1 幾何模型

圖1 燃料棒和子通道編號(hào)示意圖Fig.1 Diagram of number of rod and sub-channel

1.2 網(wǎng)格劃分

網(wǎng)格劃分使用STAR-CCM+中自動(dòng)生成非結(jié)構(gòu)化網(wǎng)格，該工具可自動(dòng)根據(jù)局部區(qū)域的幾何特征修改附近網(wǎng)格尺寸，根據(jù)壁面網(wǎng)格質(zhì)量選取不同的y+值。為盡可能提高模擬精度，開(kāi)展了網(wǎng)格敏感性分析。網(wǎng)格劃分時(shí)選取3種網(wǎng)格數(shù)量進(jìn)行敏感性分析，在不同網(wǎng)格數(shù)量下計(jì)算得到1號(hào)加熱棒平均溫度在軸向的發(fā)展趨勢(shì)，如圖2所示。計(jì)算發(fā)現(xiàn)，3種網(wǎng)格數(shù)量在相同的計(jì)算條件下溫度發(fā)展趨勢(shì)相同，溫度模擬值非常相近。為保證計(jì)算機(jī)可負(fù)荷的情況下盡可能提高網(wǎng)格質(zhì)量，采用4 000萬(wàn)網(wǎng)格尺寸進(jìn)行計(jì)算。各計(jì)算區(qū)域的網(wǎng)格細(xì)節(jié)如圖3所示。

圖2 網(wǎng)格敏感性分析Fig.2 Mesh sensitivity analysis

1.3 湍流數(shù)值模型

在各湍流數(shù)值模型中，k-ε模型能較好地模擬充分發(fā)展的湍流流動(dòng)，k-ω模型則能較好地模擬各壓力梯度下近壁面處的流動(dòng)情況。Menter[11]提出SSTk-ω模型，不僅在近壁面處保留原始的k-ω模型，也對(duì)充分發(fā)展的湍流流動(dòng)狀況運(yùn)用k-ω模型，k方程和ω方程為：

(1)

(2)

(3)

(4)

(5)

其中：ρ為密度；Gk為湍流動(dòng)能k的產(chǎn)生項(xiàng)；Gω為比耗散率ω的產(chǎn)生項(xiàng)；t為時(shí)間；x為三維坐標(biāo)系；μ為動(dòng)力黏度；μt=ρkT，T為湍流時(shí)間尺度；σk、σω2、β*和β可通過(guò)混合函數(shù)求解得到；γ為混合系數(shù)；S為流體平均應(yīng)力張量的模；U為流體速度；y為距壁面的距離；ν為流體的運(yùn)動(dòng)黏度；CDkw為與流體的橫向擴(kuò)散項(xiàng)相關(guān)的項(xiàng)。

a——整體網(wǎng)格；b——局部網(wǎng)格圖3 流體區(qū)域和固體區(qū)域的網(wǎng)格Fig.3 Mesh around fluid region and solid region

由于鉛鉍合金普朗特?cái)?shù)較低、熱導(dǎo)率較高，與水冷卻劑相比，鉛鉍合金在包殼表面形成的熱邊界層薄于流動(dòng)邊界場(chǎng)，故該模擬計(jì)算中采用SSTk-ω湍流模型。

湍流普朗特?cái)?shù)Prt表現(xiàn)了湍流流動(dòng)的一種屬性，其數(shù)值為湍流動(dòng)量擴(kuò)散率與湍流熱擴(kuò)散率的比值，本文選取了不同學(xué)者關(guān)于Prt的研究結(jié)果[12-14]。

根據(jù)不同Prt關(guān)系式計(jì)算得到的1號(hào)加熱棒平均溫度T在軸向上的發(fā)展?fàn)顩r，如圖4所示。計(jì)算發(fā)現(xiàn)，在相同的計(jì)算條件下的溫度發(fā)展趨勢(shì)相同，溫度模擬值差異不大。由于Cheng等[14]關(guān)于Prt的研究更適用于恒定熱流密度情況下的鉛鉍冷卻劑流動(dòng)計(jì)算，故選擇其研究的Prt關(guān)系式。

圖4 湍流普朗特?cái)?shù)敏感性分析Fig.4 Turbulent Prandtl number sensitivity analysis

1.4 邊界條件

在計(jì)算過(guò)程中，燃料組件的總功率為197 kW。

經(jīng)歷過(guò)17世紀(jì)“黃金時(shí)代”的阿姆斯特丹，如今的倚重科技創(chuàng)新，重塑城市生活方式，兼顧城市發(fā)展與居民生活質(zhì)量，成為歐洲乃至全球市場(chǎng)中最具競(jìng)爭(zhēng)力綠色城市之一，是歐洲最安全的城市。

入口邊界條件為：質(zhì)量流量為19 kg/s，冷卻劑進(jìn)口溫度為473.15 K。組件外圍壁面和上、下兩個(gè)端面設(shè)為絕熱邊界條件以忽略鉛鉍合金對(duì)壁面間的熱傳遞，另外將其余面設(shè)為交界面以實(shí)現(xiàn)不同計(jì)算區(qū)域間的熱傳遞。壓力設(shè)為出口邊界條件，在加熱棒內(nèi)壁面處給定熱流密度，熱流密度qw由下式計(jì)算：

(6)

其中：N為加熱棒數(shù)量；Q為總功率。

2 模型驗(yàn)證

在Pacio等實(shí)驗(yàn)中，若干數(shù)量的熱電偶布置在沿軸向發(fā)展的3個(gè)不同位置處，其對(duì)應(yīng)高度分別為繞絲螺距H的1/6、11/6和15/6倍。將各溫度測(cè)量點(diǎn)的模擬值與實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比，以驗(yàn)證幾何建模與數(shù)值模型的適用性，結(jié)果如圖5所示。圖5中，Rod1表示1號(hào)加熱棒，SCH-1表示1號(hào)子通道。

由圖5可見(jiàn)，隨著冷卻劑充分發(fā)展，溫度模擬值與實(shí)驗(yàn)結(jié)果的差距減小。在眾多模擬值與實(shí)驗(yàn)結(jié)果的對(duì)比中，子通道溫度的最大相對(duì)誤差為4.12%，加熱棒溫度的最大相對(duì)誤差為2.84%。由于增加考慮固體域內(nèi)的導(dǎo)熱，本文結(jié)果與Chai等[15]的結(jié)果相比，尤其在子通道處，模擬值與實(shí)驗(yàn)結(jié)果更為接近。

時(shí)，Pe與Nu的關(guān)系圖5 不同高度處加熱棒溫度和子通道溫度的實(shí)驗(yàn)結(jié)果與模擬值對(duì)比Fig.5 Comparison of temperature at different heights between simulation value and experimental data

3 燃料棒束熱工水力特性分析

3.1 分析算例

本文旨在研究繞絲棒束關(guān)于進(jìn)口質(zhì)量流量和功率的敏感性，擬定的質(zhì)量流量和功率工況列于表2。當(dāng)功率為197 kW時(shí)，調(diào)整組件入口的質(zhì)量流量，分別為Pacio等實(shí)驗(yàn)參數(shù)的0.25、0.5、1和2倍。當(dāng)組件入口質(zhì)量流量為19 kg/s時(shí)，將組件功率設(shè)為Pacio等實(shí)驗(yàn)參數(shù)的0.6、0.8、1和1.2倍。

3.2 結(jié)果分析

垂直于軸向、距進(jìn)口2.25倍螺距(z=738 mm)和2.5倍螺距(z=820 mm)時(shí)，工況3下截面橫向流速度分布如圖6所示。由圖6可見(jiàn)，隨繞絲沿軸向高度的發(fā)展，冷卻劑在子通道內(nèi)發(fā)生交混，使得截面上的速度峰值從左上區(qū)域向右上區(qū)域轉(zhuǎn)移。由于邊角子通道的當(dāng)量直徑較大，因此摩擦阻力較小，該區(qū)域的冷卻劑速度則相對(duì)較大。繞絲的存在使得子通道間隙寬度減小，鉛鉍合金的流通面積減小，則其流通速度增大。

表2 擬定的質(zhì)量流量和功率工況Table 2 Selected case of mass flow rate and power

圖7為與圖6對(duì)應(yīng)高度的溫度截面，由于速度的不均勻性分布使得橫截面上的溫度分布呈不均勻狀態(tài)。在邊角通道周圍加熱功率較少，且鉛鉍合金速度較大，鉛鉍合金得不到充分加熱而呈現(xiàn)出較低的溫度分布。而在速度較低的子通道區(qū)域，由于較低的鉛鉍合金速度分布，鉛鉍合金可得到充分加熱，故圖7中中心區(qū)域的溫度較高。

a——z=738 mm；b——z=820 mm圖6 工況3下的截面橫向流速度分布Fig.6 Velocity distribution of cross section transverse flow in case 3

a——z=738 mm；b——z=820 mm圖7 工況3下的截面橫向流溫度分布Fig.7 Temperature distribution of cross section transverse flow in case 3

圖8示出工況3下截面橫向流速度矢量。由圖8a可看出，邊緣子通道左上部區(qū)域中冷卻劑的流動(dòng)強(qiáng)度明顯低于右下部區(qū)域，圖8b中邊緣子通道左上部區(qū)域中冷卻劑的流動(dòng)強(qiáng)度明顯高于右上區(qū)域，這是因?yàn)槔@絲在加熱棒上沿軸向方向呈逆時(shí)針向上旋轉(zhuǎn)，不同高度的截面上橫向流動(dòng)受到方位角的影響，使速度的方向隨繞絲的轉(zhuǎn)動(dòng)而發(fā)生變化，最大速度位置也發(fā)生周期性旋轉(zhuǎn)。在同一平面中，鉛鉍合金橫向流速度在中心通道的分布較為平緩，而在邊角位置的變化較劇烈。圖9示出冷卻劑流動(dòng)過(guò)程中的流線圖。

a——z=738 mm；b——z=820 mm圖8 工況3下的截面橫向流速度矢量Fig.8 Velocity vector of cross section transverse flow under case 3

圖9 工況3下的速度流線圖Fig.9 Streamline plot in case 3

考慮到棒束的對(duì)稱性，如圖10所示，選取棒1、4、11、12研究溫度隨軸向高度的變化。對(duì)于處在棒束中心區(qū)域的燃料棒，如棒1、4，其表面平均溫度沿軸向的發(fā)展呈穩(wěn)定的周期性震蕩，這是因?yàn)槔@絲的存在使中心通道間的鉛鉍合金呈穩(wěn)態(tài)交叉流動(dòng)，從而使溫度呈現(xiàn)穩(wěn)定的周期性震蕩。對(duì)于處在棒束邊緣區(qū)域的加熱棒，如棒11、12，邊角子通道內(nèi)的二次流形態(tài)較為復(fù)雜，故未觀察到加熱棒溫度沿軸向的周期性變化。

圖10 工況3下不同加熱棒的溫度沿軸向的變化Fig.10 Temperature of different heating rods vs. axial position in case 3

將沿加熱棒周向的平均溫度作為包殼溫度，研究質(zhì)量流量對(duì)加熱棒溫度分布的影響。選取棒1、4、11和12為研究對(duì)象，采用式(7)對(duì)每根加熱棒的包殼溫度進(jìn)行無(wú)量綱化：

(7)

其中：T′為無(wú)量綱的加熱棒包殼溫度；T為加熱棒包殼溫度；Tc、Th分別為加熱棒包殼溫度的初始值和最大值。

在總功率相同的條件下，不同質(zhì)量流量下的溫度分布如圖11所示。由圖11可見(jiàn)，質(zhì)量流量的增加使得進(jìn)口位置處加熱棒的溫度變化更加劇烈，隨鉛鉍合金在軸向的發(fā)展，質(zhì)量流量對(duì)加熱棒中后段的溫度變化幅度影響較小。

為研究燃料組件內(nèi)冷卻劑的流動(dòng)特性，將子通道i與子通道j交界面中心位置處的橫向速度定義為ωi，j(i

由圖12可見(jiàn)，由于絕熱外壁的存在和同一截面中繞絲相對(duì)交界面的方向變化，ω2，3和ω27，28分別在軸向呈現(xiàn)出不同的周期性變化。在劃分網(wǎng)格過(guò)程中，由于繞絲向加熱棒圓心處移動(dòng)0.1 mm，鉛鉍合金流經(jīng)該區(qū)域時(shí)，其流通面積急速降低，冷卻劑橫向交混速度則迅速增加，形成異常峰值。另外，隨鉛鉍合金進(jìn)口質(zhì)量流量的增大，其交混程度變得更劇烈。

a——棒1；b——棒4；c——棒11；d——棒12圖11 加熱棒溫度的質(zhì)量流量敏感性分析Fig.11 Sensitivity analysis of mass flow rate for heating rod temperature

圖12 ωi，j沿軸向高度隨質(zhì)量流量的變化Fig.12 ωi，j change along axial height with mass flow rate

在入口質(zhì)量流量均為19 kg/s的條件下，研究加熱棒溫度的功率敏感性，如圖13所示。圖13結(jié)果與圖10的相同，功率大小只影響加熱棒溫度的大小，但功率變化并不影響加熱棒平均溫度隨軸向的變化。圖14示出ωi，j沿軸向高度隨功率的變化。由圖14可見(jiàn)，ω2，3和ω27，28在不同功率下均呈現(xiàn)出相同的周期性變化，由此可見(jiàn)功率大小對(duì)冷卻劑的橫向交混速度影響較小。

a——棒1；b——棒11圖13 加熱棒溫度的功率敏感性分析Fig.13 Power sensitivity analysis of heating rod temperature

圖14 ωi，j沿軸向高度隨功率的變化Fig.14 ωi，j change along axial height with power

4 結(jié)論

本文基于商用計(jì)算流體力學(xué)軟件STAR-CCM+，使用流固耦合的方法對(duì)帶有繞絲結(jié)構(gòu)的19棒束鉛鉍組件建立幾何模型，選取適當(dāng)?shù)木W(wǎng)格數(shù)量，確定合適的湍流模型。模擬計(jì)算結(jié)果表明，由于邊角子通道與中心通道間的摩擦阻力存在差異，使得鉛鉍合金在橫截面上的速度呈不均勻性分布，進(jìn)而導(dǎo)致溫度分布的不均勻性。在進(jìn)口質(zhì)量流量為19 kg/s和總功率為197 kW工況下，中心區(qū)域和邊緣區(qū)域的加熱棒溫度分別呈現(xiàn)穩(wěn)態(tài)周期性震蕩和非穩(wěn)態(tài)震蕩。另外，隨質(zhì)量流量的增加，子通道間橫向交混增大。功率變化對(duì)通道間的橫向交混速度的影響較小，對(duì)冷卻劑溫度的橫向分布也無(wú)明顯差異。