李治剛，安 萍，潘俊杰，盧 川，蘆 韡,*，楊洪潤(rùn)

(1.中國(guó)核動(dòng)力研究設(shè)計(jì)院，四川 成都 610213； 2.核反應(yīng)堆系統(tǒng)設(shè)計(jì)技術(shù)重點(diǎn)實(shí)驗(yàn)室，四川 成都 610213)

在壓水堆中燃料溫度、慢化劑密度、冷卻劑溫度等熱工參數(shù)會(huì)影響中子的截面，進(jìn)而影響堆芯中子通量和功率分布，而功率又反過(guò)來(lái)影響堆芯燃料溫度、冷卻劑溫度等熱工參數(shù)，即壓水堆中子物理與熱工之間存在強(qiáng)烈的反饋現(xiàn)象[1]，因此目前在壓水堆的設(shè)計(jì)和安全分析中均須考慮物理-熱工耦合計(jì)算，以便更加真實(shí)地模擬堆芯穩(wěn)態(tài)和瞬態(tài)過(guò)程。自20世紀(jì)80年代以來(lái)，國(guó)內(nèi)外針對(duì)壓水堆物理-熱工耦合現(xiàn)象開(kāi)展了大量研究，提出了較多有效的耦合計(jì)算方法，如Picard迭代法[2]、JFNK迭代法[3]，并研發(fā)了大量三維物理-熱工耦合計(jì)算軟件，如CRONOS/FLICA4、PARCS/TRACE[4]、CSSS[5]、DYN3D/ATHLET[6]、NALSANMT/CORBA-Ⅳ[7]等。

目前典型的堆芯物理-熱工耦合計(jì)算程序常采用外耦合方式實(shí)現(xiàn)，且常采用1種中子物理或熱工水力計(jì)算方法求解守恒方程。為研究不同中子物理、熱工水力計(jì)算方法對(duì)壓水堆物理-熱工耦合計(jì)算的影響，本文通過(guò)內(nèi)耦合方式實(shí)現(xiàn)物理-熱工耦合計(jì)算，研制典型壓水堆堆芯物理-熱工耦合計(jì)算軟件ARMcc，其中采用節(jié)塊展開(kāi)法(NEM)[8]和格林函數(shù)節(jié)塊法(NGFM)[9]求解中子擴(kuò)散方程，分別采用單通道計(jì)算模型[10]和一維圓柱導(dǎo)熱計(jì)算[11](采用有限差分法或有限體積法求解)冷卻劑溫度和燃料溫度，并采用典型壓水堆耦合基準(zhǔn)題NEACRP[12-13]彈棒初始參數(shù)對(duì)軟件穩(wěn)態(tài)計(jì)算功能進(jìn)行驗(yàn)證。

1 理論模型

1.1 穩(wěn)態(tài)中子擴(kuò)散方程求解方法

多群穩(wěn)態(tài)中子擴(kuò)散方程[14-15]可表達(dá)為：

(1)

式中：χg為裂變譜；Dg為第g群中子擴(kuò)散系數(shù)；Σfg′為第g′群裂變宏觀(guān)截面；φg為第g群中子注量率；υg′為第g′群平均裂變中子數(shù)；Σsg′→g為第g′群到第g群的散射截面；ΣRg為第g群的移除截面；k為特征值或有效增殖因數(shù)。

目前式(1)的典型求解方法包括有限差分法、節(jié)塊法、有限元方法等，其中節(jié)塊法由于具有較高的效率和精度，在商業(yè)軟件中被廣泛應(yīng)用。節(jié)塊法的重要特點(diǎn)是對(duì)式(1)進(jìn)行橫向積分，將三維問(wèn)題的求解變成聯(lián)立求解3個(gè)一維問(wèn)題。在網(wǎng)格m內(nèi)，對(duì)給定的坐標(biāo)方法u(交替等于x,y,z)，對(duì)式(1)沿與u方向垂直的另兩個(gè)方向v和w進(jìn)行積分，得到3個(gè)一維方程(為便于描述，全文略去節(jié)塊編號(hào)m)：

g=1,…,G;u=x,y,z;u≠v≠w

(2)

式中，φgu、Qgu和Lgu分別為橫向積分通量、源項(xiàng)(包括裂變?cè)错?xiàng)和散射源項(xiàng))和橫向泄漏項(xiàng)。

典型的節(jié)塊法有NEM、解析節(jié)塊法(ANM)和NGFM等，不同節(jié)塊法的差異在于φgu、Qgu和Lgu等近似處理方式。本文將采用四階NEM和第二類(lèi)邊界條件NGFM進(jìn)行多群穩(wěn)態(tài)中子擴(kuò)散方程的求解，具體的理論詳見(jiàn)文獻(xiàn)[8-9]。

1.2 熱工水力計(jì)算方法

1) 單通道計(jì)算模型

冷卻劑單相單通道的守恒方程為：

(3)

式中：ρ為冷卻劑密度，kg/m3；h為冷卻劑焓，J/kg；u為冷卻劑流動(dòng)速度，m/s；q為體積熱流密度，W/m3。

對(duì)式(3)在網(wǎng)格m進(jìn)行z方向的積分得到：

(4)

根據(jù)入口邊界條件，采用式(4)可依次計(jì)算得到每個(gè)網(wǎng)格出口焓hout，則網(wǎng)格m的平均焓為:

(5)

2) 一維圓柱導(dǎo)熱模型

一維圓柱導(dǎo)熱穩(wěn)態(tài)模型的守恒方程為：

(6)

式中：λ為導(dǎo)熱系數(shù)，W/(cm·℃-1)；T為溫度，℃。

典型的壓水堆燃料元件幾何形式如圖1所示。在徑向上燃料芯體劃分為n個(gè)網(wǎng)格，氣隙為1個(gè)網(wǎng)格，包殼劃分為2個(gè)網(wǎng)格。

圖1 典型壓水堆燃料元件幾何形式Fig.1 Geometry of fuel element in typical PWR

(1) 有限差分法

采用有限差分法(DIF)對(duì)式(6)在網(wǎng)格i處進(jìn)行離散，得到網(wǎng)格i的溫度計(jì)算關(guān)系式：

(7)

式中：λ為材料的導(dǎo)熱系數(shù)；ri為第i個(gè)網(wǎng)格的半徑；Δri為第i個(gè)網(wǎng)格的間距。關(guān)于燃料中心、芯體外表面、包殼內(nèi)外表面的處理方式參見(jiàn)文獻(xiàn)[16]。

(2) 有限體積法

采用有限體積法對(duì)式(6)進(jìn)行離散，表達(dá)形式如下：

(λi-1/2ri-1/2+λi+1/2ri+1/2)·

(8)

式中，ri+1/2和ri-1/2分別為第i個(gè)網(wǎng)格左右邊界的半徑。

式(7)、(8)均為三對(duì)角矩陣，采用高斯-賽德?tīng)柕?jì)算得到燃料元件的徑向溫度分布，采用羅蘭公式加權(quán)計(jì)算燃料有效溫度。

2 驗(yàn)證與分析

基于上述中子物理和熱工水力計(jì)算理論，采用C/C++語(yǔ)言開(kāi)發(fā)了反應(yīng)堆堆芯多物理耦合計(jì)算軟件(ARMcc)。本文分別采用三維LWR基準(zhǔn)題[17]和NEACRP-L-335基準(zhǔn)題[12-13,18]對(duì)程序中子擴(kuò)散方程求解和物理-熱工耦合計(jì)算進(jìn)行驗(yàn)證。

2.1 中子擴(kuò)散方程求解驗(yàn)證

三維LWR基準(zhǔn)題的幾何布置及材料截面參數(shù)見(jiàn)文獻(xiàn)[17]。采用NEM和NGFM對(duì)三維LWR基準(zhǔn)題進(jìn)行模擬，堆芯有效增殖因數(shù)keff結(jié)果及偏差列于表1，三維LWR基準(zhǔn)題相對(duì)功率分布及偏差如圖2所示。

表1 keff結(jié)果及偏差Table 1 Result and deviation of keff

由表1和圖2可知，本軟件采用的NEM和NGFM方法模擬三維LWR基準(zhǔn)題均具有較高的精度。keff與參考結(jié)果相比，NEM和NGFM方法計(jì)算的keff偏差均在10 pcm以?xún)?nèi)；堆芯徑向相對(duì)功率與參考結(jié)果相比，采用NEM方法的最大偏差為0.553%，采用NGFM方法的最大偏差為1.03%。

圖2 三維LWR基準(zhǔn)題相對(duì)功率分布及偏差Fig.2 Relative power distribution and deviation of 3D LWR benchmark

2.2 物理-熱工耦合計(jì)算驗(yàn)證

圖3 PWR NEACRP彈棒基準(zhǔn)題堆芯模型Fig.3 Core model of PWR NEACRP rod ejection benchmark

NEACRP彈棒基準(zhǔn)題是1991年由Finnemann等[12-13]建立的，包含壓水堆和沸水堆兩種堆型，主要用于輕水堆堆芯三維物理-熱工耦合軟件的驗(yàn)證。其中壓水堆基準(zhǔn)題參考了典型壓水堆的幾何尺寸和運(yùn)行狀態(tài)，堆芯布置了157個(gè)燃料組件和1圈徑向反射層，軸向分為18層，如圖3所示[18]。根據(jù)初始功率水平和控制棒彈出位置的不同，共包含6種工況：A1，1/4堆芯、中心控制棒在熱態(tài)零功率(HZP)彈出；A2，1/4堆芯、中心控制棒在熱態(tài)滿(mǎn)功率(HFP)彈出；B1，1/4堆芯、外圍1組控制棒在熱態(tài)零功率彈出；B2，1/4堆芯、外圍1組控制棒在熱態(tài)滿(mǎn)功率彈出；C1，1/2堆芯、外圍1組控制棒在熱態(tài)零功率彈出；C2，1/2堆芯、外圍1組控制棒在熱態(tài)滿(mǎn)功率彈出。

關(guān)于NEACRP基準(zhǔn)題幾何尺寸劃分及材料截面參數(shù)可詳見(jiàn)文獻(xiàn)[12]，基準(zhǔn)題的參考結(jié)果由PANTHER程序采用精細(xì)時(shí)空網(wǎng)格計(jì)算，在1993年發(fā)布初版、1997年發(fā)布修訂版，計(jì)算結(jié)果包括穩(wěn)態(tài)和瞬態(tài)關(guān)鍵結(jié)果參數(shù)。本文將采用1997年修訂版的穩(wěn)態(tài)計(jì)算結(jié)果對(duì)本軟件物理-熱工穩(wěn)態(tài)耦合計(jì)算模塊進(jìn)行驗(yàn)證，而對(duì)于PANTHER未提供的結(jié)果參數(shù)(如堆芯組件徑向功率分布等)將采用PARCS軟件計(jì)算結(jié)果作為參考。2種物理計(jì)算方法和2種熱工計(jì)算方法組合后形成4種計(jì)算模式(NEM+VOL對(duì)應(yīng)節(jié)塊展開(kāi)法和有限體積法，NEM+DIF對(duì)應(yīng)節(jié)塊展開(kāi)法和有限差分法，NGFM+VOL對(duì)應(yīng)格林函數(shù)法和有限體積法，NGFM+DIF對(duì)應(yīng)格林函數(shù)法和有限差分法)。

1) 臨界硼濃度

采用ARMcc針對(duì)NEACRP基準(zhǔn)題計(jì)算的臨界硼濃度對(duì)比結(jié)果列于表2。與PANTHER(1997)相比，4種模式計(jì)算的6種工況的臨界硼濃度與參考結(jié)果符合較好，最大相對(duì)偏差在0.5%以?xún)?nèi)，其中在A(yíng)2和C2工況中，NEM方法計(jì)算臨界硼濃度的精度高于NGFM方法，而在剩余4個(gè)工況中，NGFM方法計(jì)算臨界硼濃度的精度高于NEM方法。

2) 堆芯組件徑向相對(duì)功率最大值

徑向相對(duì)功率最大值Fxy,max列于表3。與PARCS結(jié)果相比，ARMcc計(jì)算的堆芯組件徑向相對(duì)功率最大值與參考結(jié)果符合較好，相對(duì)偏差在0.5%以?xún)?nèi)。采用NEM方法計(jì)算的最大相對(duì)偏差-0.39%出現(xiàn)在A(yíng)2算例中，而采用NGFM方法計(jì)算的最大相對(duì)偏差0.27%出現(xiàn)在C2算例中。采用有限體積法或有限差分法對(duì)Fxy,max的影響較小。

B1算例堆芯組件徑向功率分布如圖4所示，其中在同一種物理方法中采用VOL和DIF的相對(duì)功率分布基本相同。從圖4可知，相對(duì)功率<0.9的組件功率最大偏差情況為：NEM方法為2.09%，NGFM方法為-0.14%，相對(duì)功率>0.9的組件功率最大相對(duì)偏差情況為：NEM方法為-1.48%，NGFM方法為-0.06%。

表2 不同工況下的臨界硼濃度Table 2 Critical boron concentration of different cases

表3 不同工況下的堆芯組件徑向相對(duì)功率最大值Table 3 Maximum radial relative power of core assembly of different cases

圖4 B1算例堆芯組件徑向功率分布(1/8堆芯)Fig.4 Radical power distribution of B1 case (1/8 core)

3) 堆芯軸向相對(duì)功率最大值

堆芯軸向相對(duì)功率峰值Fz,max列于表4。與PARCS結(jié)果相比，4種計(jì)算模式得到的堆芯軸向相對(duì)功率最大值與參考結(jié)果符合較好，最大相對(duì)偏差不超過(guò)0.1%。采用NEM+VOL模式計(jì)算時(shí)在滿(mǎn)功率工況A2、B2和C2中相對(duì)偏差達(dá)到最大，為0.1%。

4) 堆芯燃料最高溫度

表5列出不同工況下的堆芯燃料最高溫度。與PARCS結(jié)果相比，因零功率算例(A1、B1和C1)堆芯功率較小，燃料最高溫度與初始溫度相同，相對(duì)偏差為0.0%；在滿(mǎn)功率算例中，不同計(jì)算模式得到的燃料最高溫度存在明顯差異：1) 采用相同物理計(jì)算方法時(shí)，采用有限體積法得到的燃料最高溫度低于采用有限差分法；2) 采用相同熱工計(jì)算方法時(shí)，NEM方法計(jì)算的結(jié)果低于NGFM方法；3) NEM+VOL模式相對(duì)偏差最小(B2算例中為0.77%)，NGFM+DIF模式相對(duì)偏差最大(C2算例中為1.33%)。

表4 不同工況下的堆芯軸向相對(duì)功率峰值Table 4 Core axial relative power peak of different cases

表5 不同工況下的堆芯燃料最高溫度Table 5 Core maximum fuel temperature of different cases

5) 堆芯燃料多普勒溫度

堆芯燃料多普勒溫度影響中子截面參數(shù)，表6列出不同工況下的堆芯燃料多普勒溫度。與PARCS結(jié)果相比，零功率算例相對(duì)偏差為0.00%；在滿(mǎn)功率算例中，4種計(jì)算模式的最大相對(duì)偏差在0.5%以?xún)?nèi)，有限差分法較有限體積法更接近參考結(jié)果。

6) 堆芯出口冷卻劑最高溫度

堆芯出口溫度體現(xiàn)了熱量在軸向上的積分效果。堆芯出口冷卻劑最高溫度列于表7。與PARCS結(jié)果相比，零功率算例相對(duì)偏差為0.00%；在滿(mǎn)功率算例中，4種計(jì)算模式的最大相對(duì)偏差在0.2%以?xún)?nèi)，而NGFM方法的相對(duì)偏差小于NEM方法的，更接近參考結(jié)果。

表6 不同工況下的燃料多普勒溫度Table 6 Core fuel Doppler temperature of different cases

表7 不同工況下的堆芯出口冷卻劑最高溫度Table 7 Maximum temperature for core coolant outlet temperature of different cases

3 結(jié)論

本文針對(duì)壓水堆堆芯物理-熱工耦合計(jì)算現(xiàn)象，基于2種典型節(jié)塊法和一維圓柱導(dǎo)熱計(jì)算方法，研制了堆芯物理-熱工耦合穩(wěn)態(tài)計(jì)算程序，采用三維LWR基準(zhǔn)題對(duì)中子物理求解模塊進(jìn)行了驗(yàn)證，采用LWR彈棒基準(zhǔn)題NEACRP-L-335的穩(wěn)態(tài)結(jié)果對(duì)物理-熱工耦合穩(wěn)態(tài)計(jì)算模塊進(jìn)行了驗(yàn)證，分析了中子物理、熱工水力計(jì)算方法對(duì)物理-熱工耦合計(jì)算的影響，得到如下結(jié)論。

1) 采用NEM方法和NGFM方法在求解中子擴(kuò)散方程時(shí)具有良好的精度，對(duì)三維LWR基準(zhǔn)題模擬結(jié)果與參考結(jié)果符合較好。

2) 4種計(jì)算模式均能較為準(zhǔn)確地模擬堆芯物理-熱工耦合過(guò)程，但不同計(jì)算模式對(duì)關(guān)鍵參數(shù)的影響程度不同：對(duì)堆芯出口溫度和堆芯軸向相對(duì)功率最大值的影響較小，對(duì)堆芯臨界硼濃度、堆芯燃料最高溫度和堆芯多普勒溫度影響較大。

3) NGFM+DIF模式能更加準(zhǔn)確地模擬堆芯燃料多普勒溫度和堆芯功率分布；NGFM+VOL模式能更加準(zhǔn)確地模擬臨界硼濃度；NEM+VOL能更加準(zhǔn)確地模擬堆芯燃料最高溫度。