基于Newton-Raphson算法的COSINE熱工水力程序安注箱模塊研發(fā)

陳 俊，沙會娥，吳照國 ，孔曉寧，尹鐵男，林 萌,楊燕華,

(1.國家電投集團科學技術(shù)研究院有限公司 國家能源核電軟件重點實驗室， 北京 102209；2.上海交通大學 核能科學與工程學院，上海 200240)

核電設(shè)計軟件用于核電廠設(shè)計、安全分析及運行管理，是支撐我國先進核電技術(shù)持續(xù)創(chuàng)新發(fā)展的關(guān)鍵重大共性技術(shù)。COSINE軟件包是國家核電技術(shù)有限公司北京軟件技術(shù)中心自主研發(fā)的堆芯設(shè)計和安全分析一體化軟件包，其中熱工水力安全分析程序包括系統(tǒng)程序、子通道程序及安全殼程序。

非能動安全注入系統(tǒng)在非LOCA事故情況下，可對反應(yīng)堆冷卻劑系統(tǒng)(RCS)進行補水和硼化，在LOCA事故下可對RCS進行安全注入[1]。安注箱屬中壓安注，是核電廠重要的專設(shè)安全設(shè)施，目前國內(nèi)外針對他的實驗分析和模擬分析尚不成熟[2-4]。

安注箱模型主要與系統(tǒng)程序聯(lián)立計算。與其他程序?qū)Π沧⑾淠Ｐ吞厥庠O(shè)置不同[2，3,5-7]，COSINE系統(tǒng)程序安注箱模塊有獨立的守恒方程和求解算法，既可以與COSINE系統(tǒng)程序聯(lián)立求解，也可以安注箱模塊獨立程序求解。

本文首先介紹COSINE系統(tǒng)程序安注箱模塊的物理模型和數(shù)學算法；其次分析了COSINE安注箱模型與參考程序RELAP5結(jié)果的對比計算。

1 安注箱物理模型及算法設(shè)計

現(xiàn)有程序假設(shè)認為安注箱氣體空間(如氮氣)內(nèi)沒有水蒸氣；實際上氣空間與液空間處于動態(tài)平衡，氣相空間里存在部分通過相變轉(zhuǎn)換的蒸汽組分。本文假設(shè)氣體空間中存在蒸汽，并以飽和狀態(tài)均勻分散在氣相空間中；如果蒸汽含量足夠小，可以通過設(shè)置或者計算體積份額來判斷或者設(shè)置是否考慮蒸汽相。水蒸氣和不凝氣體共同占據(jù)氣相的體積，并具備相同的體積份額。由于液體比熱容大，密度變化可以忽略不計，整個安注過程中狀態(tài)基本不變；本文假設(shè)安注過程中液相溫度維持不變。

1.1 守恒方程

1.1.1 質(zhì)量守恒方程

質(zhì)量守恒方程主要考慮不凝氣體連續(xù)方程和液相連續(xù)方程。在一般電廠中，不凝氣體都采用氮氣，此處設(shè)計不凝氣體為后續(xù)擴展用；用戶可以選擇氮氣，也可以選其他不凝氣體作為安注氣體。相應(yīng)的質(zhì)量守恒方程表述如下。

不凝氣體的質(zhì)量守恒方程：

(1)

液相質(zhì)量守恒方程：

(2)

式中：αn——不凝氣體體積份額；

ρn——不凝氣體密度；

ρf——液相密度；

V——安注箱體積；

AL——安注箱流通面積；

Γf——液相質(zhì)量變化率；

u——安注速度。

1.1.2 能量守恒方程

由于液相在安注過程中狀態(tài)基本保持不變，認為是恒溫過程，只需建立氣相能量守恒方程：

(3)

式中：hn——不凝氣體焓；

hv——蒸汽焓；

pTK——安注箱壓力；

H——氣液相換熱系數(shù)；

QW——壁面?zhèn)鳠崃俊?/p>

1.1.3 動量守恒方程

安注過程可能安注氣相，也可能安注液相，或者氣液混合相。假設(shè)安注過程中氣液兩相具備相同流速，由于蒸汽和不凝氣體密度接近，為簡化計算，小量蒸汽部分用不凝氣體替代，相應(yīng)的動量守恒方程可以表示為：

(4)

式中：Lf,TK——安注箱內(nèi)液相高度；

Lf,L——安注管內(nèi)液柱長度；

Lg,L——安注管內(nèi)氣體長度；

λ——摩擦因子；

KL——局部損失系數(shù)。

1.2 數(shù)值算法

守恒方程采用半隱式離散方法。首先根據(jù)動量守恒方程得到速度與壓力之間的關(guān)系式，并利用此關(guān)系式得到臨時速度，用于質(zhì)量守恒方程和能量守恒方程的更新計算。速度的壓力表達式可以表述為：

(5)

式中：C——速度系數(shù)；

A——方程顯式項；

B——壓力系數(shù)。

其次聯(lián)立質(zhì)量守恒方程和能量守恒方程，建立質(zhì)量守恒方程和能量守恒方程的線性方程組；利用Newton-Raphson算法計算下一時刻氣相體積份額，不凝氣體焓，以及不凝氣體壓力。該線性方程組包括不凝氣體質(zhì)量殘差方程，液相質(zhì)量殘差方程以及氣相能量殘差方程。質(zhì)能方程求解時選取氣體體積份額，不凝氣體焓以及不凝氣體壓力為主變量。相應(yīng)的控制體線性方程組可以表示為：

(6)

式中：εc,n——不凝氣體質(zhì)量殘差表達式；

εc,f——液相質(zhì)量殘差表達式；

εe,g——氣相能量殘差表達式；

δαn——氣相體積份額變化率；

δhn——不凝氣體焓變化率；

δpn——不凝氣體壓力變化率。

2 結(jié)果與討論

根據(jù)安注箱物理模型和數(shù)學算法，將COSINE計算結(jié)果與RELAP5相對比。為排除系統(tǒng)等其他因素對安注箱計算程序影響，程序系統(tǒng)一回路處理為邊界參數(shù)輸入安注箱。安注箱主要有球形和柱形，本文選取柱形安注箱為對象。

選取定背壓工況和變背壓工況作為安注箱邊界條件進行結(jié)果驗證，其柱形安注箱幾何形狀如圖1所示，幾何參數(shù)如表1所示。

表1 柱形安注箱幾何參數(shù)Table 1 Geometry of cylindrical accumulator

圖1 柱形安注箱示意圖Fig.1 Tank of accumulator

2.1 定背壓工況

對定背壓邊界工況，設(shè)置安注箱初始壓力為6 MPa，邊界出口背壓分別設(shè)置為1 MPa和3 MPa兩種情況。安注箱的表面粗糙度設(shè)置為0.000 6，局部阻力損失系數(shù)設(shè)置為20.0。

圖2、圖3給出不同定背壓條件下安注箱流量隨時間變化趨勢圖。虛線為COSINE計算結(jié)果，實線為RELAP5計算結(jié)果(下同)。當背壓為1 MPa時(見圖2)，RELAP5計算的安注箱安注入量在110 s是趨近于零，而COSINE在120 s之后；當背壓為3 MPa時(見圖3)，RELAP5安注箱注入量在25 s左右趨于0，而COSINE在30 s左右趨于0。另對比圖2、圖3中的剛開始幾秒，RELAP會發(fā)生劇烈震蕩，COSINE計算結(jié)果穩(wěn)定性要稍好于RELAP5。

圖2 安注箱流量變化趨勢圖(背壓為1 MPa)Fig.2 Flow rate of accumulator(Pback=1 MPa)

圖3 安注箱流量變化趨勢圖(背壓為3 MPa)Fig.3 Flow Rate of accumulator(Pback=3 MPa)

圖4、圖5給出不同定背壓下安注箱內(nèi)壓力隨時間變化趨勢圖。當背壓為1 MPa時(見圖4)，RELAP5和COSINE安注箱壓力在120 s左右趨于零；當背壓為3 MPa時(見圖5)，RELAP5和COSINE安注箱壓力在30 s左右趨于0。

圖4 安注箱壓力變化趨勢圖(背壓為1 MPa)Fig.4 Tank pressure of accumulator(Pback=1 MPa)

圖5 安注箱壓力變化趨勢圖(背壓為3 MPa)Fig.5 Tank pressure of accumulator(Pback=3 MPa)

造成COSINE與RELAP5計算結(jié)果之間差異的主要可能原因為RELAP5采用了三個守恒方程：質(zhì)量守恒方程，動量守恒方程和能量守恒方程；而COSINE的安注箱模型采用四個守恒方程，兩個質(zhì)量守恒方程，動量守恒方程及能量守恒方程。具體分析如下：

(1)質(zhì)量守恒方程

RELAP5質(zhì)量守恒方程認為空氣中氣相不考慮任何相變，其質(zhì)量為常數(shù)，不考慮流動，忽略一切其他因素影響：

Mc,n=const=ρc,nV

(7)

本文安注箱程序采用如公式(1)、公式(2)所示的兩個質(zhì)量守恒方程。首先，COSINE程序安注箱模塊考慮液體的相變，當相變量很小時兩者差別并不大，一旦發(fā)生相變則RELAP5會產(chǎn)生較大誤差。

其次，COSINE考慮質(zhì)量守恒方程的對流作用；由于RELAP5沒有考慮對流項的影響，如果系統(tǒng)初始給的速度不合理，程序計算容易發(fā)散；COSINE則不會遇到此問題。因此考慮對流項可使程序不同狀態(tài)之間變換更具平滑性和穩(wěn)定性。

最后，COSINE還給出了安注箱液相質(zhì)量守恒方程，是程序計算更具適用性。

(2)能量守恒方程

RELAP5的能量守恒方程沒有考慮蒸汽部分的能量變化，而COSINE將不凝氣體(如氮氣等)和蒸汽同時考慮成氣相，并合并成統(tǒng)一能量守恒方程。此改動與質(zhì)量守恒方程考慮蒸汽的原因一致。

(3)求解算法

RELAP5將質(zhì)量守恒方程、能量守恒方程以及動量守恒方程，結(jié)合相應(yīng)的狀態(tài)方程將三個方程化簡，結(jié)合相應(yīng)的假設(shè)折合成一等價公式，并利用此公式進行計算。COSINE利用選取的四個守恒方程，在基于Newton-Raphson算法的基礎(chǔ)上結(jié)合四個守恒方程設(shè)計具備自主知識產(chǎn)權(quán)的求解方法；此方法利用現(xiàn)有的守恒方程，不進行任何假設(shè)，采用純數(shù)學的方法設(shè)計，具有通用性和擴展性。

2.2 變背壓工況

圖6給出背壓隨時間的變化趨勢，其背壓從初始4.38 MPa變化到150 s附近的0.9 MPa。

圖6 安注箱背壓變化趨勢圖Fig.6 Back pressure of accumulator

變背壓工況下安注箱注入量隨時間的變化趨勢如圖7所示。從結(jié)果上看兩者變化趨勢一致，RELAP5計算的波動性要大于COSINE計算結(jié)果。具體原因如前面所述，守恒方程選取、求解方法不同以及假設(shè)等因素所致。

圖7 安注箱注入量變化趨勢Fig.7 Flow rate of accumulator(variable back pressure)

圖8給出了安注箱壓力隨時間變化趨勢。COSINE與RELAP5計算趨勢一致，但過程中存在差異，相對誤差并不大。

圖8 安注箱壓力變化趨勢圖Fig.8 Pressure of accumulator(variable back pressure)

3 結(jié)論

本文建立安注箱的物理模型和數(shù)學算法，并在此基礎(chǔ)上編寫成程序。通過對該模型的分析，以及不同條件下與RELAP5的計算結(jié)果對比分析，其結(jié)論如下：

(1)該安注箱程序既可作為安注箱設(shè)備程序獨立計算，也可以嵌入COSINE程序作為安全箱模塊與COSINE程序整體耦合計算；

(2)通過設(shè)計不同工況下算例，將此程序與現(xiàn)有的系統(tǒng)程序RELAP5計算結(jié)果進行對比，整體趨勢吻合較好，在局部情況下較RELAP5更具穩(wěn)定性和收斂性。造成此差異的可能主要原因有以下兩點：

1)兩者守恒方程的選取和模型簡化；

2)兩者算法的設(shè)計原理不同。

本文的計算算法是COSINE熱工水力程序的外擴，為具備自主知識產(chǎn)權(quán)的COSINE熱工水力程序的開發(fā)和應(yīng)用奠定了一定的基礎(chǔ)。

致謝

感謝項目國家科技重大專項“核動力廠安全分析用計算機軟件評估基準題及共享平臺開發(fā)”資助(2019ZX06005001)。