田曉艷，陳 森，楊 寧，朱 磊，李華琪，馬騰躍，胡 攀，康小亞

(西北核技術(shù)研究院，陜西 西安 710024)

脈沖堆又稱為TRIGA堆，可用于培訓、研究和同位素生產(chǎn)，且能為抗輻射加固技術(shù)和核參數(shù)測量技術(shù)提供輻射模擬平臺。西安脈沖堆(XAPR)是一座自然循環(huán)冷卻池式研究堆，堆芯采用鈾氫鋯燃料，由于鈾氫鋯具有較大的瞬發(fā)溫度負反饋系數(shù)，反應(yīng)堆既可穩(wěn)態(tài)運行，也可脈沖運行，其穩(wěn)態(tài)運行額定功率為2 MW，脈沖運行額定功率為4 300 MW。穩(wěn)態(tài)運行時堆內(nèi)熱量由冷卻劑自然循環(huán)帶到堆水池，堆池水溫度則依靠一回路給水保持水溫恒定，最終熱阱為大氣。盡管采用鈾氫鋯燃料的脈沖堆具有很好的固有安全特性，但若發(fā)生事故將對周圍環(huán)境和人員安全產(chǎn)生較大影響，因此有必要研究脈沖堆在典型事故工況下的瞬態(tài)安全特性。到目前為止對脈沖堆的安全研究主要針對的是失水事故[1-2]，而對反應(yīng)性引入事故的研究較少[3]。XAPR已運行多年，由于設(shè)備老化，發(fā)生控制棒失控彈出等反應(yīng)性引入事故的概率增加，且可能發(fā)生安全棒卡棒導致反應(yīng)堆無法實現(xiàn)緊急停堆。

為此，本文針對XAPR特有的堆芯結(jié)構(gòu)和運行工況，建立適用于XAPR瞬態(tài)安全特性分析的數(shù)學物理模型，并基于吉爾算法自主開發(fā)瞬態(tài)安全分析程序TSAC-XAPR，研究XAPR在反應(yīng)性引入事故下的安全特性，并為事故緩解提供參考依據(jù)。

1 物理模型與程序開發(fā)

XAPR采用自然循環(huán)冷卻，即堆外池水通過冷卻劑密度差產(chǎn)生的壓力驅(qū)動流經(jīng)堆芯，將堆芯產(chǎn)生的熱量帶出，堆池水溫則依靠一回路強迫循環(huán)給水系統(tǒng)保持恒定。為提高計算效率，基于單通道模型，將堆芯簡化為平均通道和熱通道并聯(lián)的模型，簡化后的堆芯如圖1所示?；诤喕蟮亩研窘Y(jié)構(gòu)及堆芯參數(shù)建立了詳細的物理模型，包括功率求解模型、換熱模型、自然循環(huán)流量模型及空泡份額模型。

1.1 功率求解模型

堆芯功率由兩部分組成：裂變功率和衰變功率。裂變功率采用具有6組緩發(fā)中子的點堆動力學模型求解，衰變功率采用保守的Glasstone關(guān)系式計算。反應(yīng)性反饋包括燃料溫度多普勒反饋、冷卻劑溫度反饋及冷卻劑空泡反饋。

圖1 XAPR簡化堆芯示意圖Fig.1 Simplified core schematic diagram of XAPR

Glasstone關(guān)系式[4]為：

0.87((t+2×107)-0.2-(t+t0+2×107)-0.2))

(1)

式中：Nβ，γ為停堆后的堆芯衰變功率，W；N0為停堆前的運行功率，W；t為停堆后的時間，s；t0為反應(yīng)堆在功率水平N0的運行時間，s。

1.2 換熱模型

冷卻劑通道內(nèi)的流型和換熱模型將隨反應(yīng)性引入事故的進程而變化，建立相應(yīng)于事故各階段的換熱模型，從而更準確地模擬事故動態(tài)行為。

1) 單相液體換熱

根據(jù)堆內(nèi)冷卻劑流量可將單相換熱模式分為池內(nèi)單相自然對流換熱和管內(nèi)強迫對流換熱以及二者之間的過渡換熱。

(1) 池內(nèi)單相自然對流換熱

當堆內(nèi)冷卻劑流量接近于0時，堆芯換熱模式為池內(nèi)單相對流換熱。

Nu=c(Gr·Pr)nGr/Re2≥10

(2)

式中：Nu為努塞爾數(shù)；Gr為格拉曉夫數(shù)；Pr為普朗特數(shù)；Re為雷諾數(shù)；c和n為常數(shù)，由實驗確定，文獻[5]給出其推薦值。

(2) 管內(nèi)強迫對流換熱

當冷卻劑流量足夠大時，堆芯內(nèi)就形成了明顯的宏觀流動。根據(jù)Re流動又可進一步分為層流和紊流。本文采用的換熱關(guān)系式分別為Sieder和Tate關(guān)系式與Gnielinski關(guān)系式[1]。

Re<2 300

(3)

Nu=0.012(Re0.87-280)·

(4)

式中：D為管道直徑，m；L為管道長度，m；μ為動力黏度，Pa·s；下標l表示液相，w表示壁面。

(3) 過渡階段混合對流換熱

當冷卻劑流量介于池內(nèi)自然對流和管內(nèi)強迫對流之間時，堆芯內(nèi)同時存在宏觀流動和微觀流動。文獻[5]建議混合對流換熱關(guān)系式采用下式計算：

(5)

式中：NuM為混合對流努塞爾數(shù)；NuF和NuN分別為采用流動和池內(nèi)自然對流關(guān)系式計算的努塞爾數(shù)，兩種流動方向相同時取正號，相反時取負號，指數(shù)n通常取3。

2) 過冷沸騰

當包殼表面溫度超過飽和沸騰溫度而冷卻劑還未達到飽和狀態(tài)時可能會發(fā)生過冷沸騰現(xiàn)象。為計算過冷沸騰階段的換熱系數(shù)，需采用合適的關(guān)系式確定充分發(fā)展過冷沸騰起始點。本文采用Saha-Zuber[6]模型計算過冷沸騰起始點，采用Jens-Lottes關(guān)系式[7]計算過冷沸騰工況下的包殼壁面溫度。

Saha-Zuber模型采用努塞爾數(shù)Nu和斯坦頓數(shù)St作為相似準則判斷充分發(fā)展泡核沸騰起始點FDB。

其中：

(6)

(7)

(8)

當Pe<70 000時：

(9)

當Pe>70 000時：

(10)

Jens-Lottes關(guān)系式為：

ΔTsat=22.65q0.5e-p/8.7

(11)

式中：Pe為貝克萊數(shù)；q為熱流密度，W/m2；kl為熱導率，W·(m·K)-1；ΔTsub為冷卻劑過冷度，℃；cp,l為比定壓熱容，J·(kg·K)-1；G為質(zhì)量流量，kg·(m2·s)-1；p為堆芯壓力；Pa；ΔTsat為壁面過熱度，℃。

3) 飽和流動沸騰

采用Chen關(guān)系式[8]計算飽和流動沸騰工況下的換熱系數(shù)。Chen關(guān)系式認為飽和沸騰和兩相強制對流蒸發(fā)區(qū)內(nèi)均存在兩種基本傳熱模式：泡核沸騰傳熱和強制對流傳熱，可表示為：

q=hmac(Tw-Tl)+hmic(Tw-Tsat)

(12)

(13)

(14)

式中：hmac、hmic分別為強制對流換熱系數(shù)和泡核沸騰換熱系數(shù)，W/(m2·K)；Tw、Tl和Tsat分別為壁面溫度、流體溫度和流體飽和溫度，℃；F為強化因子；kl為熱導率，W·(m·K)-1；ρl、ρg分別為液體和氣體密度，kg·m-3；hlg為汽化潛熱，J·kg-1；Δpsat為壁面過熱度ΔTsat所對應(yīng)的壓差，Pa；x為含氣率；S為滑速比；σ為表面張力，N·m-1。

4) 臨界熱流密度模型

采用Bernath關(guān)系式[9]計算研究堆低壓、低流量工況下的臨界熱流密度(CHF)：

qCHF=hCHF(Tw-Tl)

(15)

(16)

(17)

式中：qCHF為臨界熱流密度，W/m2；v為流體速度，m/s；De和Dh分別為水力學直徑和流體熱當量直徑，m。

5) CHF后換熱模型

當燃料元件表面熱流密度超過CHF時，換熱模式轉(zhuǎn)化為過渡沸騰或膜態(tài)沸騰工況。

過渡沸騰采用Cheng換熱關(guān)系式[10]：

exp(0.005G+0.018 8ΔTsub)

(18)

膜態(tài)沸騰采用Bromley換熱關(guān)系式[11]：

(19)

式中：kg為氣體熱導率，W·(m·K)-1；g為重力加速度，m·s-2；Tg為氣體溫度，℃；μg為氣體黏度，Pa·s。

6) 單相氣體換熱

采用Sieder-Tate關(guān)系式[12]計算單相氣體換熱：

(20)

1.3 自然循環(huán)流量模型

XAPR依靠堆芯自然循環(huán)進行冷卻，即堆芯冷卻劑流量取決于堆芯阻力壓頭與重力驅(qū)動壓頭之差。通過對動量守恒方程沿軸向進行積分可獲得自然循環(huán)流量計算公式。

(21)

式中：A為冷卻劑通道橫截面積，m2；f為摩擦阻力系數(shù)；K為形阻系數(shù)；l為冷卻劑通道長度，m；W為堆芯流量, kg·s-1；ρ為冷卻劑密度, kg·m-3；n為控制體數(shù)量；下標i表示控制體標識。式(21)左邊第1項為流量隨時間的變化率，第2項為摩擦阻力壓降和局部阻力壓降，第3項為重力壓降。

1.4 空泡份額模型

反應(yīng)性引入事故過程中可能發(fā)生過冷沸騰和飽和沸騰，需建立過冷沸騰和膜態(tài)沸騰工況下的空泡份額模型。本文采用Saha-Zuber模型[6]計算過冷沸騰起始點和過冷沸騰空泡份額，采用漂移流模型計算飽和沸騰下的空泡份額。

1.5 程序開發(fā)

基于以上數(shù)學物理模型開發(fā)了適用于XAPR的瞬態(tài)安全特性分析程序TSAC-XAPR。為更加精確計算堆外池水溫度，將堆外池水劃分為3個控制體，分別為堆池底部控制體、堆池下降段控制體及堆池上部控制體。最終通過數(shù)學建模獲得一系列用于模擬XAPR動態(tài)特性的偏微分方程，繼而采用適用于剛性問題求解的吉爾算法[13]進行求解。所開發(fā)的TSAC-XAPR程序流程圖如圖2所示。

圖2 TSAC-XAPR程序流程圖Fig.2 Flow chart of TSAC-XAPR code

2 結(jié)果討論

2.1 穩(wěn)態(tài)計算結(jié)果

采用TSAC-XAPR程序計算了XAPR的穩(wěn)態(tài)工況，并將穩(wěn)態(tài)計算結(jié)果與反應(yīng)堆穩(wěn)態(tài)運行值及PRTHA程序穩(wěn)態(tài)計算值進行比較，以驗證所開發(fā)程序的準確性和可靠性。PRTHA是一個子通道程序，由陳立新等[14]基于CPBRA-Ⅳ而開發(fā)，結(jié)合了XAPR的自身特點，適用于低溫、低壓研究堆的熱工水力分析。此外，穩(wěn)態(tài)計算結(jié)果還將作為瞬態(tài)計算的初始輸入?yún)?shù)。圖3示出TSAC-XAPR程序計算得到的穩(wěn)態(tài)關(guān)鍵熱工水力參數(shù)。表1列出2 MW穩(wěn)態(tài)工況下TSAC-XAPR程序計算值與XAPR穩(wěn)態(tài)運行值的比較，二者相對誤差在10%以內(nèi)?？紤]到XAPR每次開堆即使在相同功率下穩(wěn)態(tài)運行，由于實驗樣品、室溫及探測系統(tǒng)誤差，最大燃料溫度、堆池水溫等參數(shù)也存在一定的偏差，所以認為TSAC-XAPR程序計算值與運行值的誤差是可接受的。表2列出2 MW穩(wěn)態(tài)工況下TSAC-XAPR與PRTHA程序計算值的對比，二者非常接近，相對誤差在5%以內(nèi)，進一步證明本文模型和程序可用于XAPR瞬態(tài)安全分析。

圖3 穩(wěn)態(tài)關(guān)鍵熱工水力參數(shù)隨時間的變化Fig.3 Key thermal-hydraulic parameter variation with time under steady state

表1 穩(wěn)態(tài)工況下TSAC-XAPR程序計算值與運行值的比較Table 1 Comparison of TSAC-XAPR code calculation value and operating value under steady state

2.2 反應(yīng)性引入事故

反應(yīng)堆在高功率運行時堆芯燃料溫度較高，最小偏離泡核沸騰點較小，該工況下發(fā)生控制棒失控抽出事故會給反應(yīng)堆內(nèi)引入較大的正反應(yīng)性，功率及燃料溫度上升較快，容易超出安全限值，從而威脅到反應(yīng)堆安全，因此本文重點研究高功率運行工況下反應(yīng)堆控制棒失控抽出后所引入的正反應(yīng)性對反應(yīng)堆安全特性的影響。由于XAPR的控制棒速率可以在0～800 mm/min之間調(diào)節(jié)，控制棒可能以不同的速率失控連續(xù)抽出。此外，如果發(fā)生控制棒失控彈出事故，將會引入一階躍反應(yīng)性。考慮到上述因素，本文進一步分析不同反應(yīng)性引入方式對反應(yīng)堆安全特性的影響。

表2 穩(wěn)態(tài)工況下TSAC-XAPR與PRTHA程序計算值的比較Table 2 Comparison of TSAC-XAPR and PRTHA code calculation value under steady state

1) 不同運行功率下的反應(yīng)性引入事故

本文分別模擬了運行功率為2、2.2、2.65及2.7 MW 4種不同工況下的反應(yīng)性引入事故瞬態(tài)。計算時設(shè)置當堆芯在不同功率水平下穩(wěn)定運行1 000 s后以二階多項式曲線形式在14.25 s內(nèi)逐漸引入價值為1 889.2 pcm的正反應(yīng)性(XAPR穩(wěn)態(tài)運行所允許的提棒引入最大反應(yīng)性)。圖4為2 MW額定功率運行工況下，外部反應(yīng)性引入曲線及各種反應(yīng)性反饋曲線隨時間的變化。由圖4可看出，外部引入的反應(yīng)性很快被XAPR自身所具有的負反饋機制所抵消，其中起主要作用的是燃料多普勒反饋效應(yīng)。圖5示出不同運行功率水平下堆芯關(guān)

圖4 反應(yīng)性隨時間的變化Fig.4 Reactivity variation with time

鍵熱工水力參數(shù)對反應(yīng)性引入事故瞬態(tài)的響應(yīng)，包括堆芯相對功率、流量、出口空泡份額、最大包殼溫度、最大燃料溫度及最小偏離泡核沸騰點(MDNBR)，其中功率為瞬態(tài)功率與初始穩(wěn)態(tài)額定功率的比值，即相對功率P/P0。從圖5a可看出，當堆芯處于額定功率水平(2 MW和2.2 MW)，堆內(nèi)引入較大的反應(yīng)性時，反應(yīng)堆功率先急劇上升隨后有所下降，最終達到新的穩(wěn)定功率。這是因為引入的正反應(yīng)性使反應(yīng)堆功率迅速上升，且堆芯燃料溫度也隨之上升，而XAPR堆芯裝載的UHZr1.65燃料具有較大的瞬發(fā)負溫度系數(shù)，因而功率的上升引發(fā)瞬發(fā)負溫度效應(yīng)，產(chǎn)生負反應(yīng)性反饋，從而使反應(yīng)堆功率有所降低，最終穩(wěn)定在相對更低的水平。而當堆芯功率增加到2.65 MW和2.7 MW時，較大的功率使得堆芯冷卻劑溫度升高，并出現(xiàn)飽和沸騰，引起冷卻劑出口空泡份額振蕩。此時由于功率2.65 MW對應(yīng)的出口空泡份額較小，所以兩相沸騰換熱有利于自然循環(huán)流量的增加，且流量振蕩較小，所以并未引起最大燃料溫度、最大包殼溫度及MDNBR的振蕩。但當功率達到2.7 MW時，堆芯冷卻劑出口空泡份額進一步增加，堆芯出現(xiàn)膜態(tài)沸騰，導致自然循環(huán)流量下降，并劇烈振蕩，流量的振蕩進一步引起堆芯功率、最大包殼溫度、最大燃料溫度及MDNBR振蕩，且由于堆芯功率過高而堆芯流量降低，導致MDNBR低于安全限值1.3。

圖5 不同運行功率水平下堆芯關(guān)鍵熱工水力參數(shù)對反應(yīng)性引入事故瞬態(tài)的響應(yīng)Fig.5 Core key thermal-hydraulic parameter transient variation under reactivity insertion accident with different operating powers

2) 不同反應(yīng)性引入方式下的事故

圖6a為本文所模擬的5種不同方式的反應(yīng)性引入，分別對應(yīng)二階多項式曲線反應(yīng)性引入、階躍反應(yīng)性引入及3種不同速率反應(yīng)性引入。計算時同樣設(shè)置當堆芯在2 MW額定功率下穩(wěn)定運行1 000 s后以不同方式引入反應(yīng)性，而最終引入的反應(yīng)性均為1 889.2 pcm(提棒引入的最大反應(yīng)性)。圖6b為不同反應(yīng)性引入方式下的反應(yīng)堆總反應(yīng)性變化。由圖6b可看出，雖最終引入的反應(yīng)性相同，但由于引入方式不同導致反應(yīng)性在達到穩(wěn)定之前的變化過程不同，其中階躍引入方式由于引入反應(yīng)性速率過快，達到瞬發(fā)超臨界，從而引起反應(yīng)性急劇升高，但由于XAPR的UHZr1.65燃料具有較大的瞬發(fā)負溫度系數(shù)，能立即產(chǎn)生較大的負

反饋效應(yīng)，所以使反應(yīng)性引入過程呈現(xiàn)幅值較大的變化。而速率最慢的線性反應(yīng)性引入方式(117 s)對應(yīng)的反應(yīng)性變化最為平緩，但最終不同反應(yīng)性引入方式下的反應(yīng)性均趨于穩(wěn)定值0。

圖7為不同反應(yīng)性引入方式下堆芯關(guān)鍵熱工水力參數(shù)對反應(yīng)性引入事故瞬態(tài)的響應(yīng)。由圖7可見，在反應(yīng)性引入過程中，由于階躍引入方式下堆芯反應(yīng)性變化最為劇烈，從而導致堆芯相對功率、流量及最大燃料溫度等參數(shù)隨之劇烈變化，尤其相對功率變化的峰值最大，而反應(yīng)性引入速率最小的線性反應(yīng)性引入方式引起的參數(shù)變化最為平緩，但由于最終引入的反應(yīng)性一致，同時XAPR具有較大的燃料負反饋效應(yīng)，所以即使在不同方式反應(yīng)性引入工況下，最終堆芯熱工參數(shù)均趨于一致，且均在安全限值內(nèi)，證明XAPR具有良好的固有安全特性，且對不同方式反應(yīng)性引入事故有較好的緩解適應(yīng)機制。但考慮到XAPR的安全裕度和穩(wěn)定性，在啟堆或提升功率的過程中應(yīng)盡量避免快速階躍引入反應(yīng)性，防止瞬發(fā)超臨界現(xiàn)象的發(fā)生，應(yīng)當選擇合理的控制棒棒速，以避免過快的反應(yīng)性引入導致功率迅速上升，堆芯冷卻劑出現(xiàn)沸騰，從而引起堆芯流量、最大燃料溫度及MDNBR等參數(shù)的振蕩。

3 結(jié)論

本文采用自主開發(fā)的瞬態(tài)安全分析程序TSAC-XAPR計算了XAPR在反應(yīng)性引入事故下的關(guān)鍵熱工水力參數(shù)的動態(tài)響應(yīng)，并分析了反應(yīng)堆在不同功率運行水平和不同反應(yīng)性引入方式下的安全特性。程序的穩(wěn)態(tài)計算結(jié)果與參考值吻合較好，證明了所開發(fā)程序的準確性和可靠性。反應(yīng)性引入事故瞬態(tài)計算結(jié)果表明：在額定功率水平下發(fā)生反應(yīng)性引入事故時，反應(yīng)堆能依靠自身固有負反饋機制促使堆芯恢復穩(wěn)態(tài)并處于安全限值以內(nèi)。而當堆芯運行功率水平提高到2.7 MW時，由于功率過高，堆芯內(nèi)出現(xiàn)兩相飽和沸騰，引起空泡份額和堆芯流量等關(guān)鍵熱工水力參數(shù)的振蕩，導致堆芯MDNBR低于安全限值，從而威脅到反應(yīng)堆安全。此外，不同反應(yīng)性引入方式的模擬計算結(jié)果顯示：階躍引入反應(yīng)性會使堆芯參數(shù)的變化最為劇烈，而線性引入反應(yīng)性的速率越慢，堆芯的參數(shù)動態(tài)響應(yīng)越平緩，所以在啟堆或提升功率時應(yīng)盡量選擇合理的控制棒提升速率，避免快速階躍引入反應(yīng)性導致瞬發(fā)超臨界現(xiàn)象的發(fā)生。

圖6 引入反應(yīng)性(a)和總反應(yīng)性(b)隨時間的變化Fig.6 Reactivity insertion (a) and total reactivity (b) variations with time

圖7 不同反應(yīng)性引入方式下堆芯關(guān)鍵熱工水力參數(shù)對反應(yīng)性引入事故瞬態(tài)的響應(yīng)Fig.7 Core key thermal-hydraulic parameter transient variation under reactivity insertion accident by different insertion modes