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

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

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

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

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

安注箱模型主要與系統程序聯立計算。與其他程序對安注箱模型特殊設置不同[2，3,5-7]，COSINE系統程序安注箱模塊有獨立的守恒方程和求解算法，既可以與COSINE系統程序聯立求解，也可以安注箱模塊獨立程序求解。

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

1 安注箱物理模型及算法設計

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

1.1 守恒方程

1.1.1 質量守恒方程

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

不凝氣體的質量守恒方程：

(1)

液相質量守恒方程：

(2)

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

ρn——不凝氣體密度；

ρf——液相密度；

V——安注箱體積；

AL——安注箱流通面積；

Γf——液相質量變化率；

u——安注速度。

1.1.2 能量守恒方程

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

(3)

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

hv——蒸汽焓；

pTK——安注箱壓力；

H——氣液相換熱系數；

QW——壁面傳熱量。

1.1.3 動量守恒方程

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

(4)

式中：Lf,TK——安注箱內液相高度；

Lf,L——安注管內液柱長度；

Lg,L——安注管內氣體長度；

λ——摩擦因子；

KL——局部損失系數。

1.2 數值算法

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

(5)

式中：C——速度系數；

A——方程顯式項；

B——壓力系數。

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

(6)

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

εc,f——液相質量殘差表達式；

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

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

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

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

2 結果與討論

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

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

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

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

2.1 定背壓工況

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

圖2、圖3給出不同定背壓條件下安注箱流量隨時間變化趨勢圖。虛線為COSINE計算結果，實線為RELAP5計算結果(下同)。當背壓為1 MPa時(見圖2)，RELAP5計算的安注箱安注入量在110 s是趨近于零，而COSINE在120 s之后；當背壓為3 MPa時(見圖3)，RELAP5安注箱注入量在25 s左右趨于0，而COSINE在30 s左右趨于0。另對比圖2、圖3中的剛開始幾秒，RELAP會發生劇烈震蕩，COSINE計算結果穩定性要稍好于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給出不同定背壓下安注箱內壓力隨時間變化趨勢圖。當背壓為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計算結果之間差異的主要可能原因為RELAP5采用了三個守恒方程：質量守恒方程，動量守恒方程和能量守恒方程；而COSINE的安注箱模型采用四個守恒方程，兩個質量守恒方程，動量守恒方程及能量守恒方程。具體分析如下：

(1)質量守恒方程

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

Mc,n=const=ρc,nV

(7)

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

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

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

(2)能量守恒方程

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

(3)求解算法

RELAP5將質量守恒方程、能量守恒方程以及動量守恒方程，結合相應的狀態方程將三個方程化簡，結合相應的假設折合成一等價公式，并利用此公式進行計算。COSINE利用選取的四個守恒方程，在基于Newton-Raphson算法的基礎上結合四個守恒方程設計具備自主知識產權的求解方法；此方法利用現有的守恒方程，不進行任何假設，采用純數學的方法設計，具有通用性和擴展性。

2.2 變背壓工況

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

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

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

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

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

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

3 結論

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

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

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

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

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

本文的計算算法是COSINE熱工水力程序的外擴，為具備自主知識產權的COSINE熱工水力程序的開發和應用奠定了一定的基礎。

致謝

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