張雪貝， 王 馳， 陳紅麗

(中國(guó)科學(xué)技術(shù)大學(xué)物理學(xué)院， 合肥 230027)

1 引 言

傳統(tǒng)上，反應(yīng)堆設(shè)計(jì)和反應(yīng)堆安全分析一般都使用最佳估算程序. 隨著計(jì)算機(jī)性能的提高和并行計(jì)算的發(fā)展[1-2]，反應(yīng)堆的高置信度模擬在反應(yīng)堆設(shè)計(jì)、方案優(yōu)化和安全分析的研究中開(kāi)始受到重視，只有考慮了反應(yīng)堆模擬中的多物理反饋，才能實(shí)現(xiàn)高置信度模擬，而中子擴(kuò)散和熱工水力耦合計(jì)算是多物理耦合計(jì)算的重要組成部分[3-8]. 反應(yīng)堆實(shí)際的運(yùn)行就是中子學(xué)和熱工水力相互反饋的一個(gè)過(guò)程. 反應(yīng)性的溫度系數(shù)(燃料溫度系數(shù)和慢化劑溫度系數(shù))是反應(yīng)堆正常運(yùn)行時(shí)實(shí)現(xiàn)反應(yīng)性控制的重要因素[9]. 要實(shí)現(xiàn)反應(yīng)堆運(yùn)行和瞬態(tài)工況的準(zhǔn)確的計(jì)算，必須要實(shí)時(shí)考慮燃料溫度、慢化劑溫度和密度對(duì)局部中子通量和系統(tǒng)反應(yīng)性的影響.

計(jì)算流體力學(xué)(Computational Fluid Dynamics, CFD)程序FLUENT能夠通過(guò)對(duì)質(zhì)量連續(xù)方程，動(dòng)量方程，能量守恒方程的耦合求解實(shí)現(xiàn)反應(yīng)堆堆芯和燃料組件的精細(xì)模擬. 因此將CFD模型與中子動(dòng)力學(xué)模型耦合來(lái)進(jìn)行反應(yīng)堆安全分析的方法受到很大的關(guān)注[10]. FLUENT中的UDS(User Defined Scalar)能夠?qū)σ活?lèi)擴(kuò)散方程運(yùn)用Fluent內(nèi)的求解器進(jìn)行求解，其在多物理場(chǎng)耦合計(jì)算已得到廣泛應(yīng)用. 西安交通大學(xué)基于FLUENT的UDF和UDS功能開(kāi)發(fā)TASNAM程序[11]，主要用于熔鹽堆的中子擴(kuò)散數(shù)值計(jì)算. 海軍工程大學(xué)基于商用軟件CFX的用戶接口編程功能，增加三維時(shí)空中子動(dòng)力學(xué)模型，并與CFD熱工水力學(xué)進(jìn)行耦合，對(duì)壓水堆穩(wěn)態(tài)下局部三維流動(dòng)行為和三維物理特性進(jìn)行了數(shù)值模擬[12].

本文基于FLUENT的UDF和UDS功能，對(duì)中子擴(kuò)散方程進(jìn)行定義，利用Gambit進(jìn)行精細(xì)建模并利用FLUENT內(nèi)的求解器在同一套網(wǎng)格下實(shí)現(xiàn)中子擴(kuò)散與熱工水力的耦合計(jì)算. 迭代計(jì)算流程如圖1所示. 為驗(yàn)證耦合計(jì)算方法的可行性和數(shù)據(jù)傳遞的正確性，對(duì)5x5壓水堆組件模型[13]進(jìn)行建模和耦合計(jì)算，并將計(jì)算結(jié)果與其它耦合程序計(jì)算結(jié)果進(jìn)行對(duì)比. 然后通過(guò)對(duì)M2LFR-1000[14]熱組件的建模和耦合計(jì)算，得到組件穩(wěn)態(tài)條件下的中子通量分布和溫度分布，將熱工水力參數(shù)(燃料最高溫度，包殼外表面最高溫度)和子通道程序(KMC-SUB)[15]的計(jì)算結(jié)果進(jìn)行比較.

圖1 耦合計(jì)算流程圖Fig.1 The flow chart of coupling calculation

2 計(jì)算方法和數(shù)學(xué)模型

2.1 FLUENT的UDS模塊

FLUENT的UDS模塊能夠利用有限容積法對(duì)一類(lèi)方程進(jìn)行定義和運(yùn)用其內(nèi)在的求解器進(jìn)行離散和求解，形式如式(1)所示：

(1)

FLUENT中方程各項(xiàng)定義如表1所示.

瞬態(tài)中子擴(kuò)散方程與UDS定義中各項(xiàng)對(duì)應(yīng)關(guān)系如式(2)所示，式(2)方程左端第一項(xiàng)為非穩(wěn)態(tài)項(xiàng)，左端第二項(xiàng)為擴(kuò)散項(xiàng)，方程右端為源項(xiàng).

對(duì)于穩(wěn)態(tài)計(jì)算，略去方程中的非穩(wěn)態(tài)項(xiàng).

方程如式(3)所示.

(2)

(3)

式(3)中φg(r)代表中子通量，單位為cm-2s-1，D代表中子擴(kuò)散系數(shù)，單位為cm，∑f為宏觀裂變截面，單位為cm-1，∑g′→g代表從g′群到g群的宏觀轉(zhuǎn)移截面，單位為cm-1，∑r為宏觀消失截面，單位為cm-1，ν為每次裂變釋放的中子數(shù) ，χg為中子的裂變譜.

有效增殖系數(shù)由式(4)計(jì)算：

(4)

表1 Fluent中的UDS對(duì)方程的定義

表2 5×5組件材料熱物性

2.2 熱工水力計(jì)算模型

冷卻劑區(qū)域能量方程如式(5)所示：

(ρcoolantcPcoolantU·T)=βcoolantU·P+

(5)

燃料區(qū)域?qū)岱匠倘缡?6)所示：

(6)

(7)

γ為每次裂變釋放的能量.

包殼和氣隙導(dǎo)熱方程如式(8)，式(9)所示：

(8)

(9)

2.3 材料溫度截面庫(kù)

通過(guò)DRAGON程序在已知材料核子密度條件下，同時(shí)考慮溫度對(duì)材料密度的影響，通過(guò)能群歸并，計(jì)算得到燃料，氣隙，包殼，冷卻劑在400, 500, 600, 700, 800,900, 1 000, 1 100, 1 200, 1 300, 1 400, 1 500, 1 600 K處的宏觀截面，通過(guò)插值得到材料在400～1 600 K范圍內(nèi)的連續(xù)溫度截面分布[16]，建立材料宏觀截面和溫度的函數(shù)關(guān)系，簡(jiǎn)寫(xiě)為式(10)，然后將宏觀截面溫度分布函數(shù)通過(guò)UDF添加到FLUENT的求解器中，通過(guò)讀取每次迭代后每個(gè)網(wǎng)格的新的溫度分布更新每個(gè)網(wǎng)格的材料宏觀截面.

∑=f(T)

(10)

3 耦合計(jì)算方法驗(yàn)證

3.1 模型介紹

本文運(yùn)用基于FLUENT的UDF和UDS功

能得到中子擴(kuò)散和熱工水力耦合計(jì)算方法對(duì)5×5壓水堆組件模型進(jìn)行計(jì)算，材料熱物性由表2給出，模型結(jié)構(gòu)和參數(shù)由圖2和表3給出，冷卻劑密度和比熱隨溫度變化由式(11)，式(12)給出.

圖2 5×5組件徑向和軸向結(jié)構(gòu)Fig. 2 Radial and axial geometric structure of the assembly

表3組件尺寸和材料參數(shù)

Tab.3Sizeandmaterialparametersoftheassembly

參數(shù)半徑燃料棒半徑4.1 mm包殼內(nèi)徑4.2 mm包殼外徑4.8 mm燃料棒間距12.5 mm燃料高度3 m底部反射層高度0.2 m頂部反射層高度0.2 m燃料UOX(2%, 4% enrichment)冷卻劑Water, 1 000 ppm boron氣隙Helium, 0.1 MPa包殼Zircaloy-2功率12.5 MW能群數(shù)2能群分界能0.625 eV

3.2 建模與網(wǎng)格劃分

使用Gambit對(duì)5×5壓水堆燃料組件模型進(jìn)行建模和網(wǎng)格劃分，燃料棒徑向網(wǎng)格劃分如圖3所示，組件徑向網(wǎng)格劃分如圖4所示，中子擴(kuò)散和熱工水力計(jì)算采用同一套網(wǎng)格，其中軸向網(wǎng)格尺寸為0.1 m，模型網(wǎng)格總量3 100 000，利用Gambit的網(wǎng)格質(zhì)量檢驗(yàn)工具對(duì)所畫(huà)網(wǎng)格進(jìn)行檢測(cè)，EquiSize Skew在0～0.4的網(wǎng)格數(shù)占99.32%，湍流計(jì)算采用k-epsilon模型.

(11)

(12)

3.3 計(jì)算邊界條件

耦合計(jì)算邊界條件如表4所示.

表4 耦合計(jì)算邊界條件

3.4 計(jì)算結(jié)果及分析

為驗(yàn)證網(wǎng)格無(wú)關(guān)性，本文分別計(jì)算了930萬(wàn)網(wǎng)格和310萬(wàn)網(wǎng)格兩套模型，得到310萬(wàn)網(wǎng)格條件下能滿足計(jì)算精度要求，組件有效增殖系數(shù)參考值為1.171 09[13]，計(jì)算得到的組件有效增殖系數(shù)為1.171. 圖5給出中子擴(kuò)散和熱工水力計(jì)算都達(dá)到收斂時(shí)2%和4%富集度燃料棒中心的軸向功率密度分布，藍(lán)色點(diǎn)和黑色點(diǎn)為其它耦合程序結(jié)果，紅色線和綠色線為本文的計(jì)算結(jié)果，4%富集度燃料棒中心最大功率密度計(jì)算值為4.25×108W/m3，2%富集度燃料棒中心最大功率密度計(jì)算值2.50×108W/m3，4%富集度燃料棒中心最大功率密度參考值為4.178×108W/m3，2%富集度燃料棒中心最大功率密度參考值為2.45×108W/m3.

由圖5可以看出，計(jì)算偏差主要在組件模型的進(jìn)口與出口處，功率密度參考值在組件進(jìn)出口處有微小的上升，這主要是由于上下反射層的影響，使得部分中子被反射進(jìn)入燃料區(qū)，進(jìn)而引起功率密度的微小上升，本文計(jì)算時(shí)由于考慮方便運(yùn)用Gambit建模和網(wǎng)格劃分，忽略了上下反射層的影響，所以燃料棒軸向功率密度在進(jìn)出口處未出現(xiàn)微小的上升. 圖6為燃料棒外徑處的溫度分布圖，圖7為燃料包殼內(nèi)徑處的溫度分布圖，圖8為燃料棒包殼外徑處的溫度分布圖，圖9為燃料組件進(jìn)出口處的冷卻劑和燃料棒的溫度分布圖. 圖10為相鄰燃料棒中心和對(duì)角燃料棒中心沿軸向的冷卻劑溫度分布. 圖11給出z=0.0 m,y=0.0 m處沿x方向的溫度分布(z=0為對(duì)稱(chēng)軸)，計(jì)算得到的4%富集度燃料中心最高溫度為1 506.97 K，2%富集度燃料中心最高溫度為1 066.42 K，4%富集度燃料中心最高溫度參考值為1 502.22 K，2%富集度燃料中心最高溫度參考值為1 047.60 K. 通過(guò)對(duì)壓水堆5×5組件模型的耦合計(jì)算驗(yàn)證了利用Fluent 的UDS和UDF功能實(shí)現(xiàn)中子擴(kuò)散和熱工水力耦合計(jì)算的可行性和數(shù)據(jù)傳遞的正確性.

圖3 燃料棒徑向網(wǎng)格劃分圖

圖4 組件徑向網(wǎng)格劃分圖Fig.4 Radial mesh of the assembly

圖5 燃料棒中心的軸向功率密度分布Fig.5 Axial power density distribution of the fuel rod center

Fig.6 Temperature distribution of the fuel rod outer diameter

Fig.7 Temperature distribution of the fuel rod internal diameter

圖8 燃料包殼外徑處溫度分布

Fig.8 Temperature distribution of the fuel cladding external diameter

圖9 組件進(jìn)出口溫度分布

Fig.9 Temperature distribution of the inlet and outlet

圖10 相鄰燃料棒中心處和對(duì)角燃料棒中心處冷卻劑溫度軸向分布Fig.10 Axial coolant temperature distribution of the contiguous fuel rod center and diagonal fuel rod center

圖11 z=0.0 m, y=0.0m處沿x方向的溫度分布Fig.11 Temperature distribution in the x axis direction (z=0.0 m,y=0.0 m)

4 耦合方法應(yīng)用

通過(guò)對(duì)5×5壓水堆組件模型進(jìn)行耦合計(jì)算，驗(yàn)證了利用FLUENT的UDF和UDS功能實(shí)現(xiàn)中子擴(kuò)散與熱工水力耦合計(jì)算是可行的，現(xiàn)利用該耦合方法計(jì)算M2LFR-1000反應(yīng)堆熱組件模型.

4.1 模型描述

M2LFR-1000是作者所在課題組設(shè)計(jì)的模塊化鉛冷快堆，組件和燃料棒結(jié)構(gòu)由圖12和圖13給出.堆芯燃料組件中的燃料棒采用正三角矩陣排布，構(gòu)成的棒束為正六邊形，被厚度為4 mm的組件盒包裹在其中，燃料棒的間距為14 mm，每個(gè)燃料組件包含有169根燃料棒. 燃料組件盒的內(nèi)對(duì)邊距為 185 mm，外對(duì)邊距為193 mm，組件中心距為198 mm. 燃料芯塊中心有直徑為 1.9 mm 的中心孔，可以在燃料棒相同線功率密度條件下降低芯塊中心溫度，提高堆芯安全裕度. 燃料芯塊外徑為 8.6 mm， 采用添加少量 MA 或不添加 MA 的 MOX燃料.

燃料芯塊和包殼之間有 0.15 mm 的間隙，其內(nèi)填充有～0.5 MPa 的 He，填充氣體壓力高于一回路的運(yùn)行壓力，這樣做一方面可以改善燃料芯塊和包殼之間的間隙導(dǎo)熱并提供惰性環(huán)境；另一方面可防止包殼因外壓與自身蠕變坍塌造成與燃料芯塊的接觸；此外，還可以檢驗(yàn)包殼是否破損. 燃料包殼為厚度0.55 mm的FMS T91，整個(gè)燃料棒的徑向直徑為10 mm，燃料棒活性區(qū)長(zhǎng)度為1 000 mm.

選擇綜合性能比較好的 FMS T91 作為其堆芯結(jié)構(gòu)和包殼的材料，F(xiàn)MS T91 的核素組成如表5所示[17]，冷卻劑為Pb.

堆芯入口冷卻劑溫度設(shè)為673.15 K，堆芯出口冷卻劑溫度設(shè)為753.15 K，在包括設(shè)計(jì)基準(zhǔn)事故(DBA)在內(nèi)的各種工況下，燃料芯塊最高溫度應(yīng)低于2 946.15 K[18]. 正常情況下，包殼最高溫度應(yīng)低于823.15 K[19]，并有足夠的安全裕度.

圖12 組件截面圖(單位:mm)

圖13 燃料棒截面圖Fig.13 Fuel rod cross section

4.2 建模和網(wǎng)格劃分

選取組件的1/6進(jìn)行建模計(jì)算，熱組件功率為3.6 MW，使用Gambit對(duì)M2LFR-1000組件模型進(jìn)行建模和網(wǎng)格劃分，其中軸向網(wǎng)格尺寸為0.05 m，模型網(wǎng)格總量為3 400 000，燃料棒徑向網(wǎng)格如圖14所示，組件徑向網(wǎng)格由圖15所示，利用Gambit的網(wǎng)格質(zhì)量檢驗(yàn)工具對(duì)所畫(huà)網(wǎng)格進(jìn)行檢測(cè)，Equi Size Skew在0～0.4的網(wǎng)格數(shù)占96.55%，說(shuō)明網(wǎng)格劃分質(zhì)量很好，湍流計(jì)算采用k-epsilon模型.

表5 FMS T91核子密度

圖14 M2LFR-1000燃料棒徑向網(wǎng)格劃分圖

Fig.14 Radial mesh of the M2LFR-1000 fuel rod

4.3 計(jì)算邊界條件和物性參數(shù)[18]

計(jì)算邊界條件由表6給出.

MOX熱導(dǎo)率:

kMOX=(1.0/2.85×0.03+0.035+

T×(0.286-0.715×0.03))+(6400/T2.5)

e(-16.35/T)

(13)

MOX密度:

(14)

圖15 M2LFR-1000組件徑向網(wǎng)格劃分圖

表6M2LFR-1000組件耦合計(jì)算邊界條件

Tab.6BoundaryconditionsofthecouplingcalculationoftheM2LFR-1000assembly

物理量邊界邊界條件類(lèi)型邊界值溫度(T)組件入口定值673 K中子通量(?)組件入口外推邊界邊界處梯度組件出口外推邊界邊界處梯度壓力(P)組件出口定值101 KPa速度(U)組件入口定值(0,0,1.66) m/s

T91熱導(dǎo)率

kT91=32.196-6.534e(-T/546.792)

(15)

T91比熱

CPT91=404.336+T(0.881+T(-2.04×

10-3+T×2.6755×10-6))

(16)

Pb熱導(dǎo)率

kPb=9.2+0.011T

(17)

Pb密度

densityPb=11441-1.2795T

(18)

Pb比熱

CPPb=176.2-4.923×10-2T+1.544×

10-5×T×T-1.524×10-6/(T×T)

(19)

Pb粘度

viscosityPb=4.55×10-4e1069/T

(20)

4.4 計(jì)算結(jié)果分析

圖16 和圖17 為組件出口處的快中子和熱中子通量分布(未歸一化)，由于燃料區(qū)的裂變，快中子集中在燃料區(qū)域，熱中子集中在冷卻劑區(qū)域；并且由于燃料區(qū)域溫度變化較大，造成此區(qū)域材料宏觀截面參數(shù)發(fā)生明顯變化，使得快中子和熱中子在燃料區(qū)有明顯變化；但是在冷卻劑區(qū)域的快中子和熱中子通量無(wú)明顯變化，主要是由于堆芯冷卻劑溫度673～753 K左右，冷卻劑Pb在此溫度范圍內(nèi)材料宏觀截面變化不大，所以對(duì)冷卻劑區(qū)域中子通量分布無(wú)明顯影響. 圖18為組件外邊界處的溫度分布；圖19為燃料芯塊中心溫度的軸向分布，由圖19可以看出燃料芯塊中心最高溫度出現(xiàn)在z=0.6 m處；圖20為z=0.6 m處，沿y軸方向上的組件溫度分布，由圖20可以看出由于中心孔的存在，使得燃料棒中心溫度分布變得平坦，減低燃料芯塊的最高溫度，提高堆芯的安全裕度. 圖21為軸向包殼外表面溫度分布. 耦合計(jì)算得到的包殼外表面最高溫度為774.6 K，燃料芯塊最高溫度為1 645.92 K，子通道程序KMC-SUB計(jì)算得到的包殼外表面最高溫度為777.28 K，燃料芯塊最高溫度為1 647.17 K.

圖16 出口快中子通量分布

圖17 出口熱中子通量分布

圖18 外邊界處溫度分布

圖19 軸向燃料中心溫度分布Fig.19 The fuel pellet centerline temperature distribution along Z axis direction

圖20 組件y方向溫度分布(z=0.6 m, x=0.0 m)Fig.20 Temperature distribution along the y axis direction (z=0.6 m, x=0.0 m)

圖21 軸向包殼外表面溫度分布Fig.21 The cladding outer surface temperature distribution along z axis direction

5 結(jié) 論

本文基于FLUENT的UDF和UDS功能，對(duì)中子擴(kuò)散方程進(jìn)行定義，對(duì)反應(yīng)堆組件模型進(jìn)行精細(xì)建模，利用FLUENT內(nèi)基于有限容積法的求解器對(duì)中子擴(kuò)散方程進(jìn)行迭代求解，同時(shí)進(jìn)行熱工水力的計(jì)算，實(shí)現(xiàn)了在同一個(gè)求解器和同一套網(wǎng)格下的中子擴(kuò)散和熱工水力的耦合計(jì)算.

通過(guò)對(duì)5×5壓水堆燃料組件模型的建模與耦合計(jì)算，得到燃料棒不同邊界處(燃料外徑，氣隙外徑，包殼外徑)處的溫度分布，以及組件的進(jìn)出口溫度. 穩(wěn)態(tài)條件下的組件有效增殖系數(shù)(1.171)與參考值(1.171 09)符合較好；將燃料棒中心軸向功率密度分布，相鄰燃料棒和對(duì)角燃料棒中心的軸向冷卻劑溫度分布，x方向(z=0.0 m，y=0.0 m)溫度分布與其他耦合程序計(jì)算結(jié)果進(jìn)行對(duì)比，驗(yàn)證了耦合計(jì)算方法的可行性和數(shù)據(jù)傳遞的正確性.

將該耦合方法應(yīng)用到M2LFR-1000中，對(duì)M2LFR-1000的熱組件進(jìn)行精細(xì)建模和耦合計(jì)算，得到組件內(nèi)的中子通量分布和溫度分布，并與子通道程序KMC-SUB計(jì)算結(jié)果進(jìn)行對(duì)比，燃料芯塊最高溫度偏差為1.25 K，包殼外表面最高溫度偏差為2.68 K，得到其熱工水力學(xué)參數(shù)都在設(shè)計(jì)限值范圍內(nèi).

下一步將考慮添加反射層的建模，然后進(jìn)行瞬態(tài)的中子擴(kuò)散與熱工水力耦合計(jì)算.