朱 帆，伍建輝，余呈剛，馬玉雯,*，陳金根,*，蔡翔舟

(1.中國科學院 上海應(yīng)用物理研究所，上海 201800；2.中國科學院 先進核能創(chuàng)新研究院，上海 201800；3.中國科學院大學，北京 100049)

考慮到材料的工作溫度、慢化吸收比、化學穩(wěn)定性和價格等因素，熔鹽堆一般采用石墨作為慢化劑，堆芯由六邊形或四邊形通道式石墨組件排布而成，但石墨的慢化能力較弱且輻照后需定期更換[1-4]。氫化鋯作為另一種較好的慢化劑材料，其慢化能力優(yōu)于石墨且具有較好的熱穩(wěn)定性和耐高溫抗輻照等特點，不過因其物理和材料特性，氫化鋯一般被制成棒狀插入熔鹽中[5-7]。因此，棒狀氫化鋯慢化釷基熔鹽堆(ZrH-MSR)在中子學和熱工水力學方面與傳統(tǒng)壓水反應(yīng)堆和石墨慢化熔鹽堆均有較大不同。與傳統(tǒng)壓水反應(yīng)堆相比，熔鹽既是裂變熱源也是冷卻劑，需依靠燃料鹽流動帶走熱量。而與石墨慢化熔鹽堆相比，ZrH-MSR為開式通道，相鄰通道間熔鹽存在橫向的質(zhì)量、動量和能量交換(橫向交混)，各熔鹽通道內(nèi)質(zhì)量流密度沿軸向可能不斷發(fā)生變化，使通道內(nèi)熔鹽溫度也產(chǎn)生相對應(yīng)的變化。

ZrH-MSR作為一種新型熔鹽堆，在核熱耦合模擬方面尚缺少適用的分析工具。本文基于蒙特卡羅粒子輸運程序MCNP與自主開發(fā)的子通道熱工水力程序SubTH，進行核熱耦合程序MCNP-SubTH的開發(fā)，可為ZrH-MSR堆芯設(shè)計提供分析和優(yōu)化手段。

1 MCNP-SubTH耦合程序計算方法

MCNP-SubTH的中子學程序由MCNP5計算程序、NJOY99截面加工程序和ENDF/B-Ⅶ核數(shù)據(jù)庫構(gòu)成，利用NJOY99生成不同溫度點下的截面數(shù)據(jù)并采用偽材料法插值，可得到較高的MCNP計算精度。為實現(xiàn)高精度的核熱耦合計算，需針對ZrH-MSR流動傳熱特性開發(fā)適用的熱工水力學程序。本文基于子通道模型開發(fā)了熱工水力學程序SubTH，該程序不僅充分考慮了ZrH-MSR相鄰通道間熔鹽存在的橫向交混問題，也考慮了燃料鹽既是內(nèi)熱源也是冷卻劑的特點。

MCNP-SubTH采用交錯網(wǎng)格的映射方式實現(xiàn)物理-熱工網(wǎng)格匹配，物理模型與熱工模型中慢化棒與慢化棒相互對應(yīng)，物理模型中熔鹽網(wǎng)格是由熱工模型中每個慢化棒周圍的4個熔鹽通道所組成，解決了穩(wěn)態(tài)核熱耦合程序因網(wǎng)格類型不同難以耦合的問題，可為ZrH-MSR提供更準確的堆芯核熱耦合計算。

1.1 中子學模型

MCNP是一款基于蒙特卡羅方法的粒子輸運計算程序，已廣泛運用于反應(yīng)堆臨界計算、輻射屏蔽設(shè)計等領(lǐng)域，其計算精度得到了國際認可[8]。因此采用MCNP對ZrH-MSR復雜的堆芯結(jié)構(gòu)建模是一合理選擇，它能為熱工程序提供準確的三維內(nèi)熱源分布。

此外，由于ZrH-MSR堆芯同時存在熱反饋和反應(yīng)性反饋，其熔鹽溫度和密度、慢化劑溫度、中子截面等會產(chǎn)生相互耦合變化。材料密度變化會影響各柵元材料的密度設(shè)定，而溫度變化則影響計算所需的中子截面數(shù)據(jù)，從而迭代計算時需更新對應(yīng)的密度和截面數(shù)據(jù)供MCNP使用。目前，針對中子截面數(shù)據(jù)更新提出了3種方法：NJOY在線加工法、多普勒展寬法和偽材料法[8]。與偽材料法相比，NJOY在線加工法需在每次迭代過程中對所有溫度點運行NJOY計算所對應(yīng)的中子截面，這種方法計算結(jié)果最為精確但所需計算時間也最多。多普勒展寬法涉及中子截面數(shù)據(jù)的高精度擬合，其擬合參數(shù)取決于中子能量和材料溫度，但該方法尚未應(yīng)用于MNCP5中。因此，綜合考慮中子截面數(shù)據(jù)精度與程序計算效率，本文選擇偽材料法，用于與溫度相關(guān)截面數(shù)據(jù)的更新。

偽材料法認為，由于共振能量附近溫度所對應(yīng)的截面近似服從平方根規(guī)律，則核素在某一溫度下的截面可用相鄰的較低溫度和較高溫度所對應(yīng)截面的加權(quán)平均值近似：

Σ(T)=wlowΣ(Tlow)+whighΣ(Thigh)

wlow=1-whigh

(1)

式中：T為目標溫度；Σ(T)為核素在溫度T下所對應(yīng)的宏觀截面；Tlow和Thigh分別為溫度T所在區(qū)間的較低和較高溫度，其對應(yīng)的宏觀截面分別為Σ(Tlow)和Σ(Thigh)；wlow和whigh分別為溫度Tlow和Thigh處的權(quán)重系數(shù)。

1.2 熱工水力學模型

針對ZrH-MSR堆芯中燃料鹽既是內(nèi)熱源也是冷卻劑，且相鄰通道間熔鹽存在橫向交混的特點，本文基于子通道模型[9]編寫了熱工水力分析程序SubTH，其計算流程如圖1所示。首先，讀入堆芯功率分布、子通道幾何參數(shù)等文件，并假設(shè)各子通道內(nèi)初始軸向壓降損失相同，計算軸向動量守恒方程得到燃料鹽軸向流動速度，同時修正各子通道內(nèi)軸向壓降，直至總質(zhì)量流量收斂。其次，進行燃料鹽質(zhì)量守恒方程和橫向動量守恒方程計算，得出各子通道內(nèi)燃料鹽的橫向流速和壓降。然后，進行能量守恒方程計算，得出每個網(wǎng)格內(nèi)的溫度分布。最后，完成所有軸向網(wǎng)格計算并判斷燃料鹽溫度收斂后，輸出燃料鹽軸向速度、密度和溫度分布等參數(shù)。

燃料鹽質(zhì)量守恒方程為：

(2)

式中：ρi為子通道i內(nèi)燃料鹽的密度；z為軸向高度；ui為子通道i內(nèi)燃料鹽的流速；Ai為子通道i的流通面積；ρ′為子通道i與j內(nèi)燃料鹽的平均密度；vij為子通道i的燃料鹽向相鄰子通道j的橫流速度；Sij為子通道i與j間的連接長度。式(2)左邊第1項為子通道i內(nèi)軸向質(zhì)量流量的變化，第2項為子通道i內(nèi)流出和流進的橫向流量的變化。

圖1 SubTH計算流程Fig.1 Calculation flow chart of SubTH

燃料鹽能量守恒方程為：

(3)

式中：hi、hj分別為子通道i與j內(nèi)燃料鹽的焓值；Ti、Tj分別為子通道i與j內(nèi)燃料鹽的溫度；h′為子通道i與j內(nèi)燃料鹽的平均焓值；qfuel-salt為燃料鹽線功率；k為導熱系數(shù)；lij為子通道i與j間的質(zhì)心距；wij為湍流交混量，wij=β·G·Sij，其中湍流交混系數(shù)β=a·Reb，系數(shù)a、b均由經(jīng)驗關(guān)系式給出[10-11]。式(3)左邊第1項為子通道i內(nèi)軸向能量的變化，左邊第2項為子通道i內(nèi)橫向能量的變化；右邊第1項為燃料鹽的裂變功率，第2～4項分別為軸向熱傳導帶走的熱量、熔鹽由于橫向溫差引起的橫流帶走的熱量、相鄰通道間熔鹽湍流交混帶走的熱量。

燃料鹽軸向動量守恒方程為：

(4)

式中：u′為子通道i與j內(nèi)燃料鹽的平均流速；pi為子通道i的軸向壓力；f為摩擦阻力系數(shù)；Dh為當量水力直徑；K為形阻系數(shù)；Δz為子通道i的軸向高度；g為重力加速度；θ為子通道i的軸線與水平面間的夾角；CT為燃料鹽橫向交混因子；w′為有效湍流交混量。式(4)左邊第1、2項之和為子通道i內(nèi)軸向壓降的變化；右邊第1項為子通道i的軸向壓力變化，第2～5項分別為子通道i的摩擦壓降、形阻壓降、提升壓降及子通道i與子通道j內(nèi)耦合引起的流阻壓降。

燃料鹽橫向動量守恒方程為：

(5)

式中：vk為子通道i的燃料鹽向相鄰的第k個子通道的橫流速度；βk為子通道i與k之間的夾角；pj為子通道j的軸向壓力；kg為橫向阻力系數(shù)。式(5)左邊第1、2項之和為子通道i內(nèi)橫向壓降的變化，右邊第1項為相鄰通道間的橫向壓力變化，右邊第2項為通道間耦合引起的橫向摩擦壓降。

1.3 耦合方案

MCNP-SubTH采用外耦合的方式，實現(xiàn)MCNP與SubTH之間功率、密度和溫度等數(shù)據(jù)的交換，其計算流程如圖2所示。

MCNP-SubTH的計算步驟如下：1) 利用NJOY程序借助ENDF/B-Ⅶ數(shù)據(jù)庫產(chǎn)生不同溫度點的截面數(shù)據(jù)庫；2) 耦合計算前劃分MCNP輸入卡中熔鹽和慢化劑區(qū)域與SubTH程序中控制體相對應(yīng)，并將每個熔鹽和慢化棒的控制體進行編號，其網(wǎng)格映射關(guān)系如圖3所示，迭代過程中可根據(jù)中子或熱工程序的模擬需求，提供各網(wǎng)格節(jié)點對應(yīng)的數(shù)據(jù)(包括功率、密度和溫度等)；3) 給定一初始溫度和密度的MCNP輸入文件，調(diào)用MCNP進行計算，將F7卡統(tǒng)計的堆芯裂變沉積能量分布存入指定文件；4) 根據(jù)物理-熱工網(wǎng)格對應(yīng)關(guān)系，將存入的裂變沉積能量分布數(shù)據(jù)進行處理，并傳遞給子通道程序輸入文件以更新輸入?yún)?shù)，隨后調(diào)用SubTH程序執(zhí)行計算，計算結(jié)束后將各節(jié)點的溫度及密度等參數(shù)分別存儲；5) 根據(jù)網(wǎng)格間的映射關(guān)系，將SubTH計算得到的熱工數(shù)據(jù)進行處理，以更新MCNP輸入卡；6) 將熔鹽溫度是否收斂作為耦合程序結(jié)束判定標準，如式(6)所示。

max(ern)≤0.001n∈(1,N)

(6)

圖2 MCNP-SubTH耦合程序計算流程Fig.2 Calculation flow chart of MCNP-SubTH coupling code

a——物理模型；b——熱工模型；c——耦合模型圖3 MCNP-SubTH網(wǎng)格映射關(guān)系Fig.3 Mesh mapping relationship in MCNP-SubTH

2 MCNP-SubTH耦合程序驗證

目前還沒有可用的ZrH-MSR耦合基準題或是已經(jīng)驗證過的ZrH-MSR穩(wěn)態(tài)核熱耦合程序，因此采用分模塊驗證方式，分別對MCNP-SubTH耦合程序的物理模塊、熱工模塊和數(shù)據(jù)交換模塊的準確性進行校驗。

2.1 物理模塊

MCNP計算偏差主要來源為偽材料法產(chǎn)生的截面數(shù)據(jù)的精確度。因此，本文選取了ZrH-MSR六邊形單柵元模型(慢化棒直徑為2.63 cm，柵元邊長為2.31 cm，燃料鹽為77.5 mol%LiF-22.5 mol%(HM)F4)來進行偽材料法的精確性驗證，分別對比了采用NJOY在線加工法和偽材料法在800、900、1 000和1 100 K溫度下的計算結(jié)果。表1列出有效增殖因數(shù)、裂變反應(yīng)率和俘獲反應(yīng)率計算結(jié)果的對比。表1中，CASE A為采用NJOY直接加工得到某一溫度點的截面數(shù)據(jù)(基準值)，CASE B為采用偽材料法以50 K為間隔進行插值得到某一溫度點的截面數(shù)據(jù)。結(jié)果表明，兩種方法得到的計算結(jié)果符合較好。

表1 有效增殖因數(shù)、裂變反應(yīng)率和俘獲反應(yīng)率計算結(jié)果的對比Table 1 Comparison of calculation result of keff, fission reaction rate and capture reaction rate

此外，為驗證偽材料法對于功率分布的影響，本文分別對比了采用NJOY在線加工法和偽材料法在800、900、1 000和1 100 K溫度下各網(wǎng)格內(nèi)的功率(軸向劃分200個控制體，每個控制體軸向高度為0.5 cm)。結(jié)果表明，在不同溫度下，偽材料法與NJOY在線加工法計算得到的各網(wǎng)格功率數(shù)據(jù)吻合較好，最大相對偏差為0.01%，證明了偽材料法可適用于ZrH-MSR核熱耦合計算。

2.2 熱工模塊

熱工模塊采用Fluent計算結(jié)果作為驗證基準。選取了3種典型的堆芯子通道劃分方案(4棒束矩形通道、7棒束六邊形通道和7棒束圓形通道)進行SubTH程序準確性驗證。使用Gambit軟件建模，并進行網(wǎng)格劃分。將Gambit劃分完成的網(wǎng)格文件導入Fluent中，設(shè)定進、出口等邊界條件，加載熔鹽(77.5 mol%LiF-22.5 mol%(HM)F4)、慢化劑(ZrH1.6)等物性參數(shù)，最后進行模擬計算。由于ZrH-MSR中熔鹽的雷諾數(shù)較高，F(xiàn)luent計算中采用標準k-ε湍流模型，使用Simple算法，收斂殘差設(shè)為10-6。

矩形通道的子通道劃分方案如圖4a所示，橫向劃分了9個子通道，軸向劃分了60個控制體。其中，矩形通道邊長為8 cm，高為150 cm，慢化棒直徑為2.63 cm，入口流速為0.5 m/s，入口質(zhì)量流量為8.7 kg/s，入口溫度為923.15 K，功率均勻分布，功率密度為120 MW/m3。網(wǎng)格劃分方案如圖4b所示，整個模型網(wǎng)格數(shù)量達到72萬。

圖4 矩形通道的子通道(a)和網(wǎng)格(b)劃分方案Fig.4 Division scheme of sub-channel (a) and mesh (b) for rectangle channel

圖5示出Fluent計算的矩形通道溫度及熔鹽軸向溫度分布。從圖5b可見，SubTH與Fluent計算結(jié)果吻合較好，最大相對偏差為0.45%。

圖5 矩形通道溫度(a)及熔鹽軸向溫度(b)分布Fig.5 Distribution of rectangle channel temperature (a) and axial temperature of molten salt (b)

2) 7棒束六邊形通道

六邊形通道的子通道劃分方案如圖6a所示，橫向劃分18個子通道，軸向劃分60個控制體。其中，六邊形通道邊長為6 cm，高為150 cm，慢化棒直徑為2.63 cm，入口流速為0.5 m/s，入口質(zhì)量流量為11.4 kg/s，堆芯入口溫度為923.15 K，功率均勻分布，功率密度為120 MW/m3。網(wǎng)格劃分方案如圖6b所示，整個模型網(wǎng)格數(shù)量達到216萬。

圖6 六邊形通道的子通道(a)和網(wǎng)格(b)劃分方案Fig.6 Division scheme of sub-channel (a) and mesh (b) for hexagon channel

圖7示出Fluent計算的六邊形通道溫度及熔鹽軸向溫度分布。由圖7b可見，1號子通道的溫度最低，8號子通道的溫度最高，SubTH與Fluent計算結(jié)果吻合較好，最大相對偏差為0.40%。

3) 7棒束圓形通道

圓形通道的子通道劃分方案如圖8a所示，橫向劃分12個子通道，軸向劃分60個控制體。其中，圓形通道直徑為12 cm、高為150 cm，慢化棒直徑為2.63 cm，入口流速為0.5 m/s，入口質(zhì)量流量為15.4 kg/s，堆芯入口溫度為923.15 K，功率均勻分布，功率密度為120 MW/m3。網(wǎng)格劃分方案如圖8b所示，整個模型網(wǎng)格數(shù)量達到223萬。

(3)課后鞏固。課后鞏固環(huán)節(jié)很重要。在這一環(huán)節(jié)要求學生進行知識點的自我總結(jié)，完成相應(yīng)的測試，進行自我檢查。開展師生交流，學生可以通過發(fā)帖提出自己的疑惑，學生之間的交流不但解決了問題，還開拓了視野。通過這些環(huán)節(jié)的實施，將課程的學習融入整個過程中，強化了過程學習，培養(yǎng)學生終生學習的品質(zhì)。

圖7 六邊形通道溫度(a)及熔鹽軸向溫度(b)分布Fig.7 Distribution of hexagon channel temperature (a) and axial temperature of molten salt (b)

圖8 圓形通道子通道(a)及網(wǎng)格(b)劃分方案Fig.8 Division scheme of sub-channel (a) and mesh (b) for circle channel

圖9示出Fluent計算的圓形通道熔鹽溫度及熔鹽軸向溫度分布。圓形通道的子通道只有內(nèi)、外兩層，對應(yīng)的熔鹽軸向溫度曲線也為兩條。從圖9b可見，SubTH與Fluent計算結(jié)果吻合較好，最大相對偏差為0.42%。

基于Fluent對SubTH計算的不同棒束結(jié)構(gòu)模型結(jié)果進行了驗證，熔鹽溫度相對偏差小于0.45%，證明了SubTH在ZrH-MSR熱工水力分析中的正確性和可靠性。

圖9 圓形通道溫度(a)及熔鹽軸向溫度(b)分布Fig.9 Distribution of circle channel temperature (a) and axial temperature of molten salt (b)

圖10 不同子通道熔鹽軸向溫度和功率的分布Fig.10 Distribution of axial temperature and power of molten salt in different sub-channels

2.3 數(shù)據(jù)交換模塊

由于耦合程序采用外耦合方式，因此需驗證迭代計算時物理與熱工模塊間數(shù)據(jù)處理及傳遞的正確性，從而確保耦合程序的可靠性。本文分別采用手動輸入和程序輸入兩種方式對上述3種典型的堆芯子通道模型進行計算。計算流程如下：首先將MCNP計算得到的功率分布手動輸入到SubTH中，計算燃料鹽溫度和密度場分布；其次將計算得到的燃料鹽溫度和密度場分布再手動輸入到MCNP中，重新計算功率分布；最后將采用手動輸入方式迭代計算后的熔鹽溫度和功率分布與直接使用耦合程序的計算結(jié)果進行對比，來驗證耦合程序的可靠性。以4棒束矩形通道為例(熱功率為1 MW，功率余弦分布)，圖10示出不同子通道熔鹽溫度和功率的分布。結(jié)果表明，兩者的計算結(jié)果符合較好，證明了耦合程序的準確性和可靠性。

3 ZrH-MSR堆芯燃料組件穩(wěn)態(tài)核熱耦合分析

基于優(yōu)化完成的ZrH-MSR堆芯模型[12]，選取六邊形燃料組件作為算例，其子通道劃分方案如圖11所示，橫向劃分78個子通道。此外，考慮到熔鹽堆在線添料會造成熔鹽軸向分布不均勻效應(yīng)，將軸向劃分640個控制體，近似認為熔鹽成分在每個細網(wǎng)格內(nèi)不變。該組件邊長為12.68 cm，高為320 cm，由37根直徑為2.63 cm ZrH1.6慢化棒組成，棒與棒之間區(qū)域為77.5 mol%LiF-22.5 mol%(HM)F4的燃料鹽，入口流速為1.5 m/s，入口質(zhì)量流量為141.0 kg/s，入口溫度為923.15 K，熱功率為20 MW。此外，基于ORNL針對TAP反應(yīng)堆的研究結(jié)果，本文假定ZrH慢化棒中存在1.25%的裂變沉積能量[13]。

圖11 ZrH-MSR燃料組件子通道劃分方案Fig.11 Sub-channel division scheme of ZrH-MSR fuel assembly

3.1 初始溫度無關(guān)性驗證

采用800、923.15和1 100 K的初始溫度對MCNP-SubTH的初始無關(guān)性進行了驗證，如圖12所示。由圖12可見，雖然初始設(shè)定溫度值不同導致計算得到的初始keff差別較大，但是經(jīng)過1次耦合迭代后，它們的keff基本達到重合，且隨迭代次數(shù)的增加，keff趨近于同一值。因此，耦合程序具有初始無關(guān)性。

3.2 溫度分布

圖13示出穩(wěn)態(tài)情況下不同子通道熔鹽軸向溫度分布。由圖13可見，隨軸向高度的增加，各子通道內(nèi)熔鹽溫度均不斷增大，但在進出口附近熔鹽溫度增加緩慢。這是由于隨軸向高度的增加，軸向功率不斷增大，從而熔鹽溫度也逐漸增加，但因組件軸向功率呈余弦分布，導致進出口處附近功率增加緩慢，從而進出口處熔鹽溫度變化緩慢。

圖12 keff隨耦合迭代次數(shù)的變化Fig.12 Change of keff with iteration

圖13 不同子通道熔鹽軸向溫度分布Fig.13 Axial temperature distribution of molten salt in different sub-channels

此外，在相同軸向高度下不同子通道間熔鹽溫度差距較大。這主要是由于組件內(nèi)功率分布不均勻效應(yīng)造成的，以1號子通道和71號子通道為例，盡管1號子通道的流通面積和濕潤周長等參數(shù)較小應(yīng)不利于傳熱，但71號子通道位于組件中心區(qū)域功率較高，而1號子通道位于組件角落功率較低(圖14)，與流通面積和濕潤周長等熱工參數(shù)相比，此時功率分布不均勻效應(yīng)占主導地位，從而導致越靠近組件中心區(qū)域的熔鹽溫度越高。同時，計算可得該六邊形燃料組件徑向、軸向功率峰因子分別為1.43和1.25，總功率峰因子為1.79，高于傳統(tǒng)的壓水反應(yīng)堆總功率峰因子(約1.40)。當進行ZrH-MSR堆芯詳細設(shè)計時，需進一步優(yōu)化組件大小、慢化棒數(shù)目及布置位置等展平堆芯功率分布，以減小各通道出口溫差。

圖14 不同子通道熔鹽軸向功率分布Fig.14 Axial power distribution of molten salt in different sub-channels

圖15示出穩(wěn)態(tài)情況下最熱熔鹽通道(71號子通道)、ZrH慢化棒外表面(19號慢化棒)軸向溫度分布和ZrH慢化棒徑向溫度分布(19號慢化棒，軸向高度為236 cm)。由圖15可見，ZrH慢化棒軸向溫度比熔鹽溫度高，這是由于一定比例的裂變能量沉積在ZrH慢化棒內(nèi)，造成ZrH慢化棒溫度會高于熔鹽溫度，通道內(nèi)流動的熔鹽對ZrH慢化棒起到冷卻作用。此外，ZrH慢化棒軸向、徑向最高溫度分別為1 077.98 K和1 078.51 K，均低于ZrH最高破壞溫度1 200 K[5]。

圖15 熔鹽、ZrH慢化棒軸向溫度分布和ZrH慢化棒徑向溫度分布Fig.15 Axial temperature distribution of molten salt and radial temperature distribution of ZrH moderator

4 結(jié)論

本文針對ZrH-MSR中熔鹽既是裂變熱源也是冷卻劑，和相鄰通道間熔鹽存在橫向交混的問題，基于子通道方法的熱工程序SubTH，耦合粒子輸運程序MCNP開發(fā)了一套用于ZrH-MSR穩(wěn)態(tài)核熱耦合分析程序MCNP-SubTH，并與Fluent計算結(jié)果進行對比，驗證了耦合程序的準確性。使用該程序?qū)α呅蝂rH-MSR燃料組件進行模擬，計算其功率分布、熔鹽和ZrH慢化棒溫度分布等數(shù)據(jù)，驗證了MCNP-SubTH耦合程序的可行性及有效性。關(guān)于材料的相容性問題，需進一步實驗研究。