李治剛，安 萍，潘俊杰，趙文博，劉盧果，宮兆虎，蘆 韡,強勝龍，賀 濤

(中國核動力研究設計院 核反應堆系統(tǒng)設計技術重點實驗室，成都610213)

由于輕水堆的中子物理與熱工水力之間存在強烈的反饋現(xiàn)象，故在進行壓水堆的堆型設計與安全分析時必須考慮物理-熱工耦合現(xiàn)象[1]。為能準確模擬反應堆物理-熱工耦合現(xiàn)象，近年來，國內外針對輕水堆開展了大量的物理-熱工耦合研究，研發(fā)了反應堆多物理耦合計算方法和模擬軟件。例如，針對堆芯的物理-熱工耦合計算程序有PARCS[2]，DYN3D[3]，ANCK[4]，NODAL3[5]，NALSANMT/CORBA-IV[6]等；包含反應堆系統(tǒng)及堆芯物理-熱工耦合的計算程序有PARCS/TRACE[7]，CTSS[8]，DYN3D/ATHLET[3]等。在對存在明顯物理-熱工反饋且堆芯功率分布不均的事故過程分析中，必須采用3維物理-熱工耦合分析方法。彈棒事故是反應性事故中最嚴重的一種[9]，可能會威脅燃料棒包殼和壓力容器的完整性，在壓水堆安全分析中必須對其進行評價。傳統(tǒng)的彈棒事故評價方法是采用2維與1維相結合的方法模擬堆芯中子動力學。近年來，美國和法國發(fā)展了3維動力學分析方法用于彈棒事故分析，如法國AREVA公司開發(fā)的SMART 3D+FLICA III+THEMIS+COMBAT系統(tǒng)[10]、美國西屋公司開發(fā)的TWINKLE+FACTRAN+LOFTRAN系統(tǒng)[11]及美國核管理委員會開發(fā)的PARCS/TRACE/ROBIN系統(tǒng)[12]，已分別用于嶺澳核電站、大亞灣核電站和秦山核電站二期的彈棒事故分析。我國也開發(fā)了多套3維時空中子動力學計算軟件，如中國核動力研究設計院開發(fā)的CORCA-K[13]及西安交通大學開發(fā)的NECP-Bamboo[14]等，這些軟件為研發(fā)自主的彈棒事故3維動力學分析系統(tǒng)提供了基礎。

本文基于堆芯3維時空中子動力學計算軟件CORCA-K、堆芯熱工水力計算軟件CORTH[15]及截面更新計算軟件PACFAC，針對典型壓水堆的物理-熱工3維瞬態(tài)耦合計算問題，自主研發(fā)了堆芯多專業(yè)耦合計算軟件COMPCS，采用壓水堆彈棒基準題NEACRP[16]和MOX/UO2[17]，對COMPCS軟件進行了驗證，并對“華龍一號”堆芯首循環(huán)壽期初和壽期末的彈棒事故進行了計算與分析。

1 理論模型

1.1 3維時空中子動力學方程求解方法

3維時空中子動力學多群中子擴散方程可表達為[18]

(1)

(2)

其中，g為能群號；i為緩發(fā)中子先驅核組號；G為能群數(shù)；χp,g和χd,g分別為瞬發(fā)中子裂變譜和緩發(fā)中子裂變譜；Dg為第g群中子擴散系數(shù)；Σf,g′為第g′群裂變宏觀截面；φg為第g群中子注量率；υg′為第g′群平均裂變中子數(shù)；Σs,g′→g為第g′群到第g群的散射截面；Σr,g為第g群的移除截面；vg為第g群中子速度；βi和β分別為第i組緩發(fā)中子份額和緩發(fā)中子總份額；Ci為第i組緩發(fā)中子先驅核濃度；λi為第i組緩發(fā)中子先驅核衰變常數(shù)。

目前，求解式(1)的典型方法主要為有限差分法、節(jié)塊法和有限元方法。其中，節(jié)塊法具有較高的計算效率和精度，在商業(yè)軟件中應用廣泛。節(jié)塊法的特點是對式(1)進行橫向積分，將3維問題的求解變成聯(lián)立求解3個1維方程。在網(wǎng)格m內，對給定的坐標方向x，對式(1)沿與x方向垂直的另2個方向(y,z)進行積分，可得到1維方程為

Qg,x(x)-Lg,x(x)，g=1,…，G

(3)

其中，φg,x，Qg,x，Lg,x分別為x方向橫向積分中子注量率、源項(包括裂變源項和散射源項)和橫向泄露項。同理可得y和z方向相應的1維方程。

CORCA-K 求解式(3)時，采用第二類邊界格林函數(shù)節(jié)塊法進行空間離散，采用二階對角線隱式龍格庫塔格式(diagonal implicit Runge Kutta，DIRK)進行時間離散[19]，詳細的處理過程可參見文獻[13]。

1.2 熱工水力計算方法

COMPCS軟件中除采用閉式單通道計算模型和圓柱瞬態(tài)導熱模型計算冷卻劑溫度和燃料溫度外，還建立了與堆芯熱工水力計算軟件CORTH對應的接口，從而實現(xiàn)對堆芯熱工參數(shù)子通道的模擬。CORTH軟件基于具有滑速比的四方程模型，適用于反應堆堆芯或棒束實驗熱工水力分析。本文主要介紹典型單通道模型的計算方法，子通道模型可參見文獻[15]。

1.2.1閉式單通道計算模型

冷卻劑單相單通道的瞬態(tài)守恒方程為

(4)

(5)

其中，ρ為冷卻劑密度， kg·m-3；h為冷卻劑質量焓， J·kg-1；u為冷卻劑流動速度， m·s-1；q為體積熱流密度， W·m-3；t為時間， s。

1.2.2圓柱瞬態(tài)導熱模型

圓柱瞬態(tài)導熱模型的守恒方程為[20]

(6)

其中，cp為材料的質量定壓熱容， J·(kg·K)-1；ρf為材料密度， kg·m-3；λ為熱導率， W·(m·K)-1；T為熱力學溫度， K。

在徑向上，將燃料芯體劃分為n個網(wǎng)格，氣隙為1個網(wǎng)格，將包殼劃分為2個網(wǎng)格。在空間上，采用有限體積法或有限差分法對式(6)進行離散。在時間上，采用全隱式格式。采用高斯-賽德爾迭代計算得到燃料元件的徑向溫度分布，采用羅蘭公式計算燃料有效溫度[21]。

1.3 物理-熱工耦合算法

COMPCS軟件的耦合計算流程，如圖1所示。采用CORCA-K控制時間迭代流程，基于CORCA-3D軟件提供的穩(wěn)態(tài)燃耗狀態(tài)庫[22]，在每個迭代步內調用CORTH或閉式單通道模型更新熱工水力參數(shù)，如冷卻劑溫度和密度，調用圓柱瞬態(tài)導熱模型計算燃料有效溫度，PACFAC軟件可基于新的狀態(tài)參數(shù)更新材料截面。在穩(wěn)態(tài)時，采用Picard迭代策略[23]，保證每個迭代步內中子物理和熱工水力均收斂后才進入下一個迭代步。

2 數(shù)值驗證

本文采用壓水堆彈棒基準題NEACRP和MOX/UO2對COMPCS軟件的物理-熱工耦合計算精度進行驗證。

圖1 COMPCS軟件的耦合計算流程Fig.1 Flow chart of coupling calculation process of COMPCS

2.1 NEACRP基準題

NEACRP基準題是1991年由Finnemann等建立的，主要用于輕水堆堆芯3維物理-熱工耦合軟件的驗證。該基準題參考了典型壓水堆的幾何尺寸和運行狀態(tài)，堆芯布置157個燃料組件和1圈徑向反射層，軸向分為18層。根據(jù)初始功率水平和控制棒彈出位置的不同，可將堆芯分為A1，A2，B1，B2，C1，C26種工況。其中，A1，B1，C1為熱態(tài)零功率(hot zero power, HZP)，A2, B2, C2為熱態(tài)滿功率(hot full power，HFP)，堆芯功率為2 775 MW，共有6組緩發(fā)中子，總緩發(fā)中子份額為0.007 6。NEACRP基準題的幾何尺寸劃分、材料截面及控制棒移動信息可參見文獻[16]。

本文利用COMPCS軟件計算了6個工況的主要參數(shù)，并與NEACRP基準題的參考結果進行比對。該參考結果是1997年發(fā)布的修訂版數(shù)據(jù)，由PANTHER程序采用精細時空網(wǎng)格計算得到。

表1列出了穩(wěn)態(tài)及瞬態(tài)主要參數(shù)的對比情況。由表1可見，臨界硼質量分數(shù)與參考結果的偏差在5×10-6以內，相對偏差均小于0.5%；堆芯相對功率峰值的最大偏差出現(xiàn)在B1和C1工況中，相對偏差分別為2.37%和2.35%，其余工況的相對偏差在1%以內；瞬態(tài)結束時堆芯相對功率最大偏差出現(xiàn)在零功率工況中， A1，B1，C1的相對偏差分別為1.57%，1.35%，1.37%；瞬態(tài)結束時燃料多普勒溫度與參考結果的相對偏差均在1%以內。

圖2 為NEACRP基準題相對功率隨時間變化情況的對比。為進行比較，圖2還給出了PARCS軟件的計算結果。由圖2可見，在滿功率工況中，COMPCS軟件計算結果與PARCS結果基本一致，堆芯相對功率峰值與參考結果的相對偏差均在1%以內；在零功率工況中，COMPCS軟件計算結果與PARCS結果趨勢一致，堆芯相對功率峰值更接近參考結果。

表1 NEACRP基準題穩(wěn)態(tài)和瞬態(tài)主要參數(shù)與參考結果的對比Tab.1Comparison of steady and transient main parameters with reference results in NEACRP benchmark

(a)A1

(b)A2

(c)B1

(d)B2

(e)C1

(f)C2圖2NEACRP基準題相對功率隨時間變化情況的對比Fig.2 NEACRP benchmark relative power vs. time

2.2 MOX/UO2基準題

由于武器級MOX燃料的中子動力學參數(shù)與典型壓水堆有顯著差別，2006年由OECD/NEA組織建立了壓水堆彈棒MOX/UO2基準題。該基準題基于西屋四環(huán)路全尺寸的壓水堆核電站模型，堆芯額定功率為3 565 MW，共有193個燃料組件(UO2燃料組件141個，MOX燃料組件52個)，可模擬平衡循環(huán)中一組控制棒的彈出過程。MOX/UO2基準題包含4個部分：第1部分為2維固定物理-熱工工況，給出了特征值、棒價值、組件和柵元精細功率；第2部分為3維滿功率工況，給出了臨界硼質量分數(shù)、組件和柵元功率；第3部分為3維零功率工況，給出了臨界硼質量分數(shù)、組件和柵元功率；第4部分為基于第3部分的瞬態(tài)工況，給出了控制棒彈棒瞬態(tài)響應參數(shù)，包括堆芯功率和反應性等。采用組件程序建模可得到MOX/UO2基準題組件的少群參數(shù)。 MOX/UO2基準題堆芯、組件幾何尺寸及控制棒移動信息可參見文獻[17]。

圖3為MOX/UO2基準題瞬態(tài)工況相對功率隨時間變化情況的對比。選取PARCS 2G軟件的計算結果作為比對參考結果，表2列出了MOX/UO2基準題穩(wěn)態(tài)工況計算結果的對比情況。由于第1部分為固定截面計算，不涉及物理-熱工耦合迭代，所以本文未對第1部分進行模擬計算。

由圖3可見， COMPCS計算的相對功率峰值為1.644，該值與PARCS 2G軟件、PARCS 8G 軟件及CORETRAN軟件計算結果的相對偏差分別為14.49%，-4.6%，1.2%。由表2可見，COMPCS軟件計算的穩(wěn)態(tài)關鍵參數(shù)中，第3部分的Fxy與參考結果的相對偏差為1.04%，其他參數(shù)的相對偏差均在1%以內。

驗證結果表明，COMPCS軟件對NEACRP基準題和MOX/UO2基準題的計算精度與國際同類軟件相當，計算精度較高。

圖3 MOX/UO2基準題瞬態(tài)工況堆芯相對功率隨時間的變化Fig.3 Relative power vs. time in MOX/UO2benchmark transient condition

表2 MOX/UO2基準題穩(wěn)態(tài)工況的計算結果對比Tab.2Comparison of results of MOX/UO2 benchmark in steady condition calculated by COMPCS and PARCS 2G

3 “華龍一號”堆芯彈棒事故模擬

“華龍一號”(HPR1000)反應堆是我國具有自主知識產(chǎn)權的第三代核電壓水堆[24-25]，堆芯由177個AFA3G燃料組件構成，每個組件包含呈17×17方形排列的264根燃料棒、24根導向管和1根測量管。堆芯控制棒組件按功能可分為控制棒組和停堆棒組，控制棒布置可參見文獻[24]?？刂瓢艚M由功率補償棒(G1，G2，N1，N2)和溫度調節(jié)棒(R)構成。停堆棒組(SA，SB，SC)的功能是提供反應堆停堆所必需的負反應性?？刂瓢艚M的提出價值由大到小依次為G1，G2，N1，N2，R，但在實際運行中，功率補償棒采用疊步方式插入或提出，R棒只能在插入限之內運行，G1棒全部插入堆芯的概率最高。因此，本文將模擬G1棒處于全插、其他控制棒處于全提時，G1棒在首循環(huán)壽期初(beginning of life, BOL)和壽期末(end of life, EOL)從堆芯彈出后堆芯相對功率及溫度等參數(shù)隨時間的變化情況。計算中設置4個功率水平，分別為25%FP，50%FP，75%FP，100%FP，“華龍一號”首循環(huán)燃耗由CORCA-3D軟件計算得到[22]。

圖4為G1棒在0.1 s內從全插位彈出堆芯后，堆芯相對功率隨時間的變化情況。由圖4可見，在相同功率水平下，堆芯瞬態(tài)功率水平在壽期末和壽期初隨時間的變化趨勢是一致的，在壽期末4個初始功率水平下的堆芯瞬態(tài)功率峰值均高于壽期初，在初始功率水平為50%FP時，壽期末堆芯瞬態(tài)功率峰值達到壽期初的1.85倍。

圖5和圖6分別給出了G1棒在0.1 s內從全插位彈出堆芯后，堆芯燃料芯塊溫度峰值和包殼溫度峰值隨時間的變化情況。

(a)BOL

(b)EOL圖4堆芯相對功率隨時間的變化Fig.4 Core relative power vs. time

(a)BOL

(b)EOL圖5燃料芯塊溫度峰值隨時間的變化Fig.5 Fuel core peak temperature vs. time

(a)BOL

(b)EOL圖6包殼溫度峰值隨時間的變化Fig.6 Cladding peak temperature vs. time

由圖5和圖6可見，在壽期初，隨著時間的推進，燃料芯塊溫度峰值和包殼溫度峰值均呈現(xiàn)先增加后穩(wěn)定在一定水平上的變化趨勢；隨著功率水平的增加，包殼溫度峰值升高的速率也增加。在壽期末，功率水平在50%FP以下時，燃料芯塊及包殼溫度峰值的變化趨勢與壽期初相似，而當功率水平在50%FP以上時，燃料芯塊及包殼溫度峰值呈現(xiàn)先快速增加、到達峰值后又逐漸減小、最后穩(wěn)定在一定水平上的變化趨勢。

在壓水堆中，隨著燃耗的加深，UO2燃料芯塊會發(fā)生腫脹，并釋放裂變氣體，導致燃料熱導率隨著燃耗的加深呈下降趨勢。此外，隨著堆芯中鈾的消耗和钚的積累，堆芯中子動力學也會發(fā)生明顯變化。在相同功率水平下，壽期末由于緩發(fā)中子份額的下降，使堆芯即使在與壽期初相同的反應性下，中子增長速度也相對增加，這是圖4(b)中功率峰增長速度快于壽期初的主要原因。

為分析彈棒時間對彈棒過程的影響，針對堆芯功率有顯著變化的壽期末滿功率工況，設置彈棒時間te分別為0.1， 1.0， 3.0 s，計算了堆芯相對功率、燃料芯塊溫度峰值和包殼溫度峰值隨時間的變化情況，結果分別如圖7和圖8所示。

圖7 G1棒不同彈棒時間下壽期末滿功率工況堆芯相對功率隨時間的變化Fig.7 Core relative power vs. time after G1 ejected withdifferent ejection time at EOL HFP condition

(a)Fuel core peak temperature

(b)Cladding peak temperature圖8G1棒不同彈棒時間下壽期末滿功率工況燃料溫度峰值隨時間的變化Fig.8 Fuel temperature peak vs. time after G1 ejectedwith different ejection time at EOL HFP condition

由圖7和圖8可見，彈棒時間從0.1 s增加到3 s時，堆芯相對功率峰值從4.958減小到1.36，最終穩(wěn)定在1.172，燃料芯塊溫度峰值和包殼溫度峰值增加的速率呈現(xiàn)下降趨勢。這是因為彈棒引入的反應性是相同的，堆芯相對功率峰值、芯塊溫度峰值和包殼溫度峰值最終均趨于各自的穩(wěn)定值。

由于COMPCS軟件未與系統(tǒng)程序耦合，無法提供一回路壓力峰值，所以本文只選擇堆芯參數(shù)對彈棒過程進行評價。計算結果表明，堆芯燃料芯塊溫度峰值遠低于芯塊熔化溫度，包殼溫度峰值低于包殼熔化溫度，因此，在不同功率水平和不同彈棒時間下，G1棒彈棒事故不會威脅燃料包殼的完整性。

4 結論

本文研制的COMPCS軟件在模擬壓水堆堆芯物理-熱工耦合穩(wěn)態(tài)工況及瞬態(tài)彈棒過程時，計算結果與NEACRP和MOX/UO2基準題的參考結果符合良好，具有較高的精度；不同初始功率水平、燃耗狀態(tài)和彈棒時間對“華龍一號”彈棒過程有明顯影響，在壽期末、初始功率水平較高及彈棒時間較短時，堆芯相對功率上升的速率更快，會形成功率峰值及燃料溫度峰值；“華龍一號”首循環(huán)G1棒彈出的典型工況下，燃料芯塊溫度峰值和包殼溫度峰值均低于設計限值，不會威脅燃料包殼的完整性。