基于CFD方法的鉛鉍冷卻燃料棒束的熱工水力特性分析

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

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

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

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

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

1 正常工況模擬

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

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

1.1 幾何模型

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

1.2 網格劃分

網格劃分使用STAR-CCM+中自動生成非結構化網格，該工具可自動根據局部區域的幾何特征修改附近網格尺寸，根據壁面網格質量選取不同的y+值。為盡可能提高模擬精度，開展了網格敏感性分析。網格劃分時選取3種網格數量進行敏感性分析，在不同網格數量下計算得到1號加熱棒平均溫度在軸向的發展趨勢，如圖2所示。計算發現，3種網格數量在相同的計算條件下溫度發展趨勢相同，溫度模擬值非常相近。為保證計算機可負荷的情況下盡可能提高網格質量，采用4 000萬網格尺寸進行計算。各計算區域的網格細節如圖3所示。

圖2 網格敏感性分析Fig.2 Mesh sensitivity analysis

1.3 湍流數值模型

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

(1)

(2)

(3)

(4)

(5)

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

a——整體網格；b——局部網格圖3 流體區域和固體區域的網格Fig.3 Mesh around fluid region and solid region

由于鉛鉍合金普朗特數較低、熱導率較高，與水冷卻劑相比，鉛鉍合金在包殼表面形成的熱邊界層薄于流動邊界場，故該模擬計算中采用SSTk-ω湍流模型。

湍流普朗特數Prt表現了湍流流動的一種屬性，其數值為湍流動量擴散率與湍流熱擴散率的比值，本文選取了不同學者關于Prt的研究結果[12-14]。

根據不同Prt關系式計算得到的1號加熱棒平均溫度T在軸向上的發展狀況，如圖4所示。計算發現，在相同的計算條件下的溫度發展趨勢相同，溫度模擬值差異不大。由于Cheng等[14]關于Prt的研究更適用于恒定熱流密度情況下的鉛鉍冷卻劑流動計算，故選擇其研究的Prt關系式。

圖4 湍流普朗特數敏感性分析Fig.4 Turbulent Prandtl number sensitivity analysis

1.4 邊界條件

在計算過程中，燃料組件的總功率為197 kW。

經歷過17世紀“黃金時代”的阿姆斯特丹，如今的倚重科技創新，重塑城市生活方式，兼顧城市發展與居民生活質量，成為歐洲乃至全球市場中最具競爭力綠色城市之一，是歐洲最安全的城市。

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

(6)

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

2 模型驗證

在Pacio等實驗中，若干數量的熱電偶布置在沿軸向發展的3個不同位置處，其對應高度分別為繞絲螺距H的1/6、11/6和15/6倍。將各溫度測量點的模擬值與實驗結果進行對比，以驗證幾何建模與數值模型的適用性，結果如圖5所示。圖5中，Rod1表示1號加熱棒，SCH-1表示1號子通道。

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

時，Pe與Nu的關系圖5 不同高度處加熱棒溫度和子通道溫度的實驗結果與模擬值對比Fig.5 Comparison of temperature at different heights between simulation value and experimental data

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

3.1 分析算例

本文旨在研究繞絲棒束關于進口質量流量和功率的敏感性，擬定的質量流量和功率工況列于表2。當功率為197 kW時，調整組件入口的質量流量，分別為Pacio等實驗參數的0.25、0.5、1和2倍。當組件入口質量流量為19 kg/s時，將組件功率設為Pacio等實驗參數的0.6、0.8、1和1.2倍。

3.2 結果分析

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

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

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

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可看出，邊緣子通道左上部區域中冷卻劑的流動強度明顯低于右下部區域，圖8b中邊緣子通道左上部區域中冷卻劑的流動強度明顯高于右上區域，這是因為繞絲在加熱棒上沿軸向方向呈逆時針向上旋轉，不同高度的截面上橫向流動受到方位角的影響，使速度的方向隨繞絲的轉動而發生變化，最大速度位置也發生周期性旋轉。在同一平面中，鉛鉍合金橫向流速度在中心通道的分布較為平緩，而在邊角位置的變化較劇烈。圖9示出冷卻劑流動過程中的流線圖。

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

考慮到棒束的對稱性，如圖10所示，選取棒1、4、11、12研究溫度隨軸向高度的變化。對于處在棒束中心區域的燃料棒，如棒1、4，其表面平均溫度沿軸向的發展呈穩定的周期性震蕩，這是因為繞絲的存在使中心通道間的鉛鉍合金呈穩態交叉流動，從而使溫度呈現穩定的周期性震蕩。對于處在棒束邊緣區域的加熱棒，如棒11、12，邊角子通道內的二次流形態較為復雜，故未觀察到加熱棒溫度沿軸向的周期性變化。

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

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

(7)

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

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

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

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

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

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

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

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 結論

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