劉碧帆,黃廣源,呂 欣,尹俊連,*,王德忠,趙劍剛,鄢夢琪

(1.上海交通大學(xué) 機械與動力工程學(xué)院,上海 200240;2.中廣核研究院有限公司,廣東 深圳 518026)

氮氣加壓具有非能動調(diào)節(jié)的優(yōu)點,廣泛應(yīng)用于反應(yīng)堆中。非能動氮氣安注箱是VVER[1]、AP1000[2]等第3代核電堆型的重要專設(shè)安全設(shè)施。失水事故(LOCA)中,當一回路壓力降至安注箱預(yù)設(shè)壓力時,氮氣安注箱可以非能動地向反應(yīng)堆內(nèi)注入冷卻水,確保反應(yīng)堆安全。氮氣穩(wěn)壓器在先進一體化小型壓水堆系統(tǒng)中具有極大應(yīng)用前景[3-4],其利用充入氮氣的可壓縮性非能動地調(diào)節(jié)一回路的壓力,具有結(jié)構(gòu)簡單可靠、運行能耗低、可快速啟動等優(yōu)點。但氮氣加壓帶來了新的問題與挑戰(zhàn)。氮氣加壓的冷卻劑中的溶解氮為飽和狀態(tài)。LOCA中,隨著壓力的持續(xù)下降,溶解氮進入過飽和狀態(tài),以微氣泡的形式大量析出積聚,導(dǎo)致傳熱惡化,阻滯安注流量,影響反應(yīng)堆安全[5-6]。為了保證反應(yīng)堆安全,防止溶解氮的析出,需要對溶解氮的遷移和析出過程進行建模。

本文針對氮氣析出的具體過程,建立氮氣壁面析出模型,解釋析出源項的物理意義。對氮氣的溶解平衡、對流擴散、壁面析出、氣泡上浮進行數(shù)學(xué)建模,基于MATLAB編寫一維集總參數(shù)氮氣遷移析出求解程序,對Sarrette實驗中SBLOCA泄壓工況下溶解氮遷移析出特性進行研究,獲得溶解氮濃度和含氣率變化。本研究建立的預(yù)測方法可為其他含氮系統(tǒng)內(nèi)氮氣析出的預(yù)測提供參考。

1 數(shù)學(xué)物理建模

本研究對氮氣遷移析出過程中的各物理過程進行了數(shù)學(xué)物理建模。LOCA泄壓工況中,氮氣遷移析出主要經(jīng)歷了4個物理過程:溶解平衡、對流擴散、壁面析出和氣泡上浮(圖1)。由于溶解平衡的影響,液面處的溶解氮濃度始終為對應(yīng)壓力的平衡濃度。隨著壓力逐漸降低,液面處的溶解氮濃度也逐漸降低。由于濃度梯度和液相脈動,溶解氮在液相中擴散。液相的運動對溶解氮起到對流輸運作用。在壁面附近,溶解氮向氣核內(nèi)析出氮氣。氮氣泡逐漸生長,最終脫離壁面上浮。

圖1 泄壓實驗中的物理過程Fig.1 Physical process in pressure relief experiment

1.1 氮氣質(zhì)量守恒方程

考慮如圖2所示的控制體V,在t時刻長度為δx、橫截面積為A的控制體V中氮氣的質(zhì)量為ρC(1-α)Aδx。圖2中:ρ為控制體中水的密度;C為溶解氮的濃度;α為空泡份額;Aδx為控制體的體積。在一維情況下,溶解氮只能在左右邊界上通過對流作用和擴散效應(yīng)進出控制體,且控制體內(nèi)氮氣滿足如下質(zhì)量守恒:

圖2 控制體中的氮氣質(zhì)量守恒Fig.2 Nitrogen mass conservation control volume

Aδx((1-αt+δt)ρt+δtCt+δt-(1-αt)ρtCt)=

(1)

式中:上標表示時間,下標L、R分別表示左、右邊界;u為液相速度;D為分子擴散系數(shù)Dm與渦流擴散系數(shù)Deddy之和。等號左側(cè)表示控制體內(nèi)氮氣質(zhì)量隨時間的變化率,等號右側(cè)第1項、第2項表示單位時間內(nèi)通過對流和擴散作用分別在左邊界進入及在右邊界離開控制體的氮氣質(zhì)量,第3項表示由于析出源項減少的氮氣質(zhì)量。

當時間微元δt與控制體尺寸δx足夠小,且AL=AR時,得到一維情況下描述溶解氮質(zhì)量守恒方程:

(2)

式中:等號左側(cè)表示單位長度的控制體內(nèi)溶解氮在單位時間內(nèi)的質(zhì)量變化;等號右側(cè)第1項為對流引起的質(zhì)量變化,其物理意義為溶解氮由于液相運動而跟隨運動的溶解氮質(zhì)量,第2項為擴散引起的質(zhì)量變化,第3項為析出引起的質(zhì)量變化。

同理可得描述氮氣對流和析出的質(zhì)量守恒方程為:

(3)

式中:ρg為氣相密度;ug為氣相速度;等號左側(cè)表示單位長度的控制體內(nèi)氮氣在單位時間內(nèi)的質(zhì)量變化;等號右側(cè)第1項為對流引起的質(zhì)量變化,第2項為析出引起的質(zhì)量變化。

1.2 溶解模型

氮氣溶解的動態(tài)平衡主要存在于氣液界面。Jirka[8]通過實驗測量發(fā)現(xiàn),在十分靠近界面的區(qū)域中溶解氣體的濃度滿足亨利定律:

pB=HBCB

(4)

式中:pB為平衡壓力;HB為亨利系數(shù);CB為飽和濃度。

溶解度數(shù)值由Baranenko等[9]的數(shù)據(jù)插值得到。

1.3 擴散模型

D為分子擴散系數(shù)Dm與渦流擴散系數(shù)Deddy之和:

D=Dm+Deddy

(5)

分子擴散系數(shù)表達式為:

(6)

式中:kB為玻爾茲曼常數(shù);T為熱力學(xué)溫度;nSE為Stokes-Einstein常數(shù);η為水的黏度;a為氮氣分子的動力學(xué)半徑。

渦流擴散系數(shù)的表達式為:

(7)

式中:νt為湍流黏度;Sct為湍流施密特數(shù),Sct取值1.36[10]。

1.4 壁面析出模型

1) 析出的主要機理

在泄壓過程中壓力變化較大,溶解氮氣達到過飽和狀態(tài),因此需要額外考慮氮氣的析出及其對氮氣遷移的影響。初步估算表明,泄壓過程中最大過飽和度不會超過10,遠不足以克服在連續(xù)的液相介質(zhì)中形成穩(wěn)定的氣液界面所需的能量壁壘。因此本文僅考慮第Ⅳ類非均相核化過程[11],即認為系統(tǒng)固體壁面上已經(jīng)穩(wěn)定存在著大于臨界尺寸的微氣核,當系統(tǒng)處于過飽和狀態(tài)后,液相中溶解的氮氣通過擴散作用進入氣核并使其長大直至脫離壁面。下面針對這一過程進行物理建模。

2) 氣泡成核條件

固體壁面上穩(wěn)定存在的微氣核半徑往往在微納米量級,其形貌可以近似簡化為球帽。下面推導(dǎo)微氣核生長所需的過飽和度。

根據(jù)楊-拉普拉斯公式,氣泡內(nèi)外壓差為:

(8)

式中:σ為對應(yīng)溫度壓力下的表面張力系數(shù);R為氣泡半徑;pg為氣泡內(nèi)氣體壓力;p0為液相壓力(系統(tǒng)壓力)。根據(jù)亨利定律,此時氣液界面上的平衡濃度Cg為:

(9)

式中:Hg為對應(yīng)溫度壓力下的亨利系數(shù);C0=p0/Hg,為系統(tǒng)當?shù)販囟葔毫ο聦?yīng)的平衡氮氣濃度,即該處的氮氣溶解度。

此時過飽和度ζ為:

(10)

初始氣核半徑一般在微米量級[12-14]。將SBLOCA中壓安注工況(p0=5.5 MPa,σ=0.070 78 N/m)代入式(10)可得:在初始氣核半徑為100 μm時,過飽和度僅為2.57×10-4,當初始氣核半徑為1 μm時,臨界過飽和度為2.57×10-2?？梢?在壁面微氣泡尺度下,過飽和析出所需的過飽和度非常小,不同初始氣核半徑對臨界過飽和度的絕對值影響可忽略,可近似認為當系統(tǒng)一旦達到過飽和狀態(tài)后,壁面上就會發(fā)生核化析出。

3) 氣泡生長模型

氣泡生長主要受擴散作用控制,用無限大過飽和流體中氣泡的自由長大過程近似地描述[12]:

(11)

(12)

根據(jù)理想氣體狀態(tài)方程,氣泡內(nèi)的氣相密度ρg可表示為:

(13)

式中,ρg0為p0壓力下對應(yīng)的氣相密度。

將式(8)、(9)、(13)代入式(12)得:

(14)

該式具有R′=f(t,R)的形式,可通過顯式的龍格-庫塔方法對R進行數(shù)值求解。

為了分析半徑R(t)隨時間變化的趨勢,可考慮將(C∞-Cg)/ρg視為常數(shù),獲得形式簡單的解析解。如圖3所示,數(shù)值計算結(jié)果表明,氣泡生長周期尺度(數(shù)十毫秒)遠小于析出過程的時間尺度(數(shù)十秒)。過飽和度在1個氣泡周期內(nèi)幾乎沒有變化。氣泡尺寸R在10-6量級,2σ/(p0R)?1,ρg=ρg0(1+2σ/p0R)≈ρg0。因此可將(C∞-Cg)/ρg視為常數(shù),從而獲得式(12)的解析解:

圖3 SBLOCA工況中不同深度的氣泡生長曲線Fig.3 Bubble growth curves at different depths in SBLOCA

(15)

4) 氣泡脫落模型

考慮Jones等[11]和Klausner等[14]提出的力平衡模型。由于事故工況下回路的剪切效應(yīng)較弱,認為析出氣泡的脫離主要是浮升力引起的。

如圖4所示,在豎直固體壁面上,氣泡脫離時浮升力與三相接觸線上的表面張力在垂直方向上達到平衡:

圖4 氣泡在豎直壁面上的受力平衡Fig.4 Force balance of a bubble on vertical wall

Fb=Fσ

(16)

(17)

Fb=(ρl-ρg)gV

(18)

式中:θ1、θ2分別為氣液固三相接觸線的前進接觸角和后退接觸角;dw為三相接觸線的直徑。

當核化點上的氣泡析出長大達到臨界直徑脫落后,將其重置為初始尺寸,重新開始新一輪的析出長大以及脫落。

5) 析出質(zhì)量源項

析出氮氣質(zhì)量源項S的物理意義為單位時間內(nèi)單位體積的控制體中由于氣泡在壁面析出而變化的氮氣質(zhì)量,用單位體積內(nèi)的氣泡數(shù)目乘以氣泡體積的變化率表示為:

(19)

若忽略表面張力系數(shù)σ隨時間的變化,有:

(20)

將式(14)代入可得:

(21)

式中,R為該計算時刻的析出氣泡半徑。

1.5 氣泡上浮模型

氣泡析出之后在水中上浮,受到浮力、曳力、虛擬質(zhì)量力、壁面潤滑力等力的作用。由于阻力的影響,氣泡進入主流后作加速度逐漸減小的變加速運動,最后以某一特定速度勻速上升。Kamp和Kraume等[15]根據(jù)力平衡推導(dǎo)出氣泡運動方程:

(22)

式中:u為氣泡運動速度;ρc為連續(xù)相密度;ρd為離散相(氣泡)密度;d為氣泡直徑;CVM為虛擬質(zhì)量系數(shù),取值為0.785[16];CD為曳力系數(shù)[17],可表示為:

(23)

2 數(shù)值求解與程序編制

2.1 數(shù)值方法

質(zhì)量守恒方程的對流項采用一階迎風(fēng)格式離散。質(zhì)量守恒方程的擴散項采用兩次中心差分格式進行離散。質(zhì)量守恒方程的析出質(zhì)量源項由式(19)、(21)和氣泡半徑R代數(shù)求解得到,R通過式(14)由顯式的龍格-庫塔方法數(shù)值求解得到。式(14)是非線性的常微分方程,具有R′=f(t,R)的形式,通過顯式的龍格-庫塔方法對R進行數(shù)值求解,累積誤差為Δt4階,在時間步長較小時計算準確度很高。

2.2 程序編制

根據(jù)氮氣的數(shù)學(xué)物理模型和相應(yīng)的數(shù)值方法,采用MATLAB編寫了一維集總參數(shù)氮氣遷移析出求解程序,計算流程如圖5所示。該程序考慮了氮氣的的溶解平衡、對流擴散、壁面析出和氣泡上浮過程,可獲得一維管道中每一時刻每一深度的溶解氮濃度和空泡份額。輸入初始條件和邊界條件后,程序開始計算。根據(jù)氮氣壁面析出模型,求解當前時間步氣泡半徑,若氣泡半徑小于氣泡脫離半徑,則基于當前氣泡半徑和質(zhì)量守恒方程求解溶解氮濃度和空泡份額。若氣泡半徑大于等于氣泡脫離半徑,則判定該氣泡脫離,將氣泡半徑重置為初始半徑之后再求解溶解氮濃度和空泡份額。到達設(shè)定的仿真時間后停止計算。

圖5 計算流程Fig.5 Calculation process

3 驗證工況

選擇Sarrett等[7]的SBLOCA泄壓實驗進行驗證。Sarrett的實驗在1個豎直放置的密閉圓管中進行,圓管的頂部與蓄壓箱、泄壓閥相連。壓力測點和熱電偶測點布置如圖6所示。圓管高8.8 m、內(nèi)徑139.7 mm、水深5 780 mm,圓管上部分充有初始壓力為6 MPa、溫度為27 ℃的氮氣。水的初始溫度為33 ℃。實驗在靜水工況下進行,可認為水的速度為0。

圖6 泄壓實驗裝置示意圖[7] Fig.6 Pressure relief experimental section geometry[7]

Loviisa核電站中氮氣蓄壓安注箱的預(yù)設(shè)壓力為5.5 MPa,注入上排管路時的水溫為40 ℃。為模擬該堆型SBLOCA下安注管道的泄壓過程,實驗過程中水的溫度維持在32 ℃,在700～787 s,氣體壓力從5.5 MPa快速降至2.0 MPa(圖7)。

圖7 泄壓實驗段壓力變化曲線[7]Fig.7 Pressure profile of experimental section under pressure relief condition[7]

4 驗證計算

4.1 參數(shù)選擇

如前文所述,對氮氣遷移析出計算影響較大且目前尚不確定的參數(shù)包括:氣泡脫離直徑、氣相運動速度、計算的時間步長Δt、初始氣核尺寸。

根據(jù)前述氣泡脫離的力平衡方程,初步計算得析出氣泡脫離直徑在10-1mm量級,與SBLOCA泄壓實驗測得氣泡脫離直徑(0.1 mm)基本吻合。取0.1 mm作為氣泡脫離直徑。

根據(jù)前述顆粒運動方程,初步計算得氣泡在主流運動的終速度在10-2m/s量級,與SBLOCA泄壓實驗測得氣相平均速度(0.01 m/s)基本吻合。對于10-1mm量級的析出氣泡,計算得到變加速運動階段在0.1 s以內(nèi),相比于氣泡從產(chǎn)生到到達液面所需的時間(數(shù)百秒)而言十分短暫。因此可忽略變加速階段,近似地認為氣相以終速度作勻速運動,取0.01 m/s作為氣相運動速度。

如圖8所示,在滿足顯式計算庫朗特數(shù)的前提下,隨著計算時間步長的減小,氮氣析出變化過程逐漸收斂,可認為,在計算的時間步長小于0.01 s后氮氣的析出過程計算是準確的,能夠獲得時間無關(guān)解。

圖8 計算時間步長對氮氣析出的影響Fig.8 Effect of time step on nitrogen release

在經(jīng)過一段時間t之后,某核化點析出的氮氣總質(zhì)量m可表示為:

(24)

式中:T為氣泡生長周期;ρg為這段時間內(nèi)的氣泡中的氣體平均密度。

(25)

假設(shè)氣泡半徑均在微米量級,滿足1 mm>Rd>R1>R2>1 μm,且Rd=k1R1,R1=k2R2,其中k1>1,k2>1。

將R1、R2代入上式得:

(26)

由k1>1,k2>1,且:

(27)

(28)

可得到不等式關(guān)系:

(29)

將工況參數(shù)(p0=5.5 MPa,σ=0.070 78 N/m,1 mm>Rd>R1>R2>1 μm)代入,當初始氣核半徑比脫離半徑小1倍,滿足k1>2、k2>1時,代入上式可知,1.026>m2/m1>0.731,析出氮氣質(zhì)量最多只相差26.4%。當初始氣核半徑比脫離半徑小1個量級,滿足k1>10、k2>1時,1.026>m2/m1>0.965,析出氮氣質(zhì)量最多只相差3.5%?？梢?由于初始氣核半徑一般都在微米量級[13,18-19],其具體數(shù)值的選取并不會顯著影響氮氣析出質(zhì)量。出于保守設(shè)計考慮,取初始氣核半徑R0=1 μm。

4.2 初始條件與邊界條件

一維計算液面邊界節(jié)點處給定Dirichlet邊界條件以描述氮氣的溶解平衡,即C=CB。

一維計算圓管底部邊界節(jié)點的氣相控制方程離散時,不考慮流入控制體的氮氣質(zhì)量,只考慮流出控制體的氮氣質(zhì)量和析出的氮氣質(zhì)量,即滿足:

(30)

根據(jù)SBLOCA泄壓實驗的實驗數(shù)據(jù),初始時刻的空泡份額為0.3%。初始時刻的溶解氮濃度為飽和濃度。

4.3 計算結(jié)果

程序計算結(jié)果如下。溶解氮平均濃度,平均空泡份額均與Sarrett實驗值符合良好(圖9)。如圖9所示,不考慮氮氣析出,僅考慮擴散作用的情況下,溶解氮的濃度幾乎沒有變化。

圖9 SBLOCA泄壓工況溶解氮濃度和含氣率計算結(jié)果Fig.9 Dissolved nitrogen concentration and void fraction under pressure relief condition during SBLOCA

如圖10所示,對于圓管中各位置,隨著壓力的下降,濃度逐漸下降,濃度的下降速率逐漸增加。平均濃度從815 ppm下降至約656 ppm,下降了19%。各時刻,圓管中濃度隨深度的增加而線性增加,但變化幅度較小,不超過1%。

圖10 SBLOCA泄壓工況濃度-時間、濃度-深度數(shù)值模擬Fig.10 Numerical simulation results of concentration-time and concentration-depth under pressure relief condition during SBLOCA

如圖11所示,除去圓管底部,圓管中各位置隨著壓力的下降,含氣率逐漸上升,含氣率的上升速率逐漸增加。含氣率從約0.003升至約0.013。管道中析出了1.152 L的氮氣。圓管底部的含氣率-時間曲線與其他位置表現(xiàn)出不同的變化趨勢。這是因為圓管底部流入的氮氣質(zhì)量更少。在降壓初期,由于流出的氮氣質(zhì)量大于析出的氮氣質(zhì)量,圓管底部的含氣率逐漸降低;在降壓后期,由于過飽和度逐漸增加,析出速率增加,流出的氮氣質(zhì)量小于析出的氮氣質(zhì)量,圓管底部的含氣率逐漸上升。從圖11還可看出,各時刻,圓管中含氣率隨深度的增加而減小。含氣率-深度曲線可分為線性下降段和加速下降段。由于圓管底部控制體節(jié)點的含氣率求解采用了特殊的邊界條件(不考慮對流流入項),該邊界條件的影響會隨著時間逐漸從上游向下游傳遞。線性下降段是未受影響的區(qū)域,加速下降段則是已經(jīng)受到影響的區(qū)域。

圖11 SBLOCA泄壓工況含氣率隨時間和深度的變化Fig.11 Void fraction vs. time and depth under pressure relief condition during SBLOCA

5 結(jié)論

本研究針對LOCA泄壓工況下氮氣在冷卻劑中的遷移析出過程進行了理論建模及數(shù)值計算分析,計算結(jié)果與驗證實驗數(shù)據(jù)吻合良好,可得到以下結(jié)論。

2) LOCA泄壓工況下,氮氣的壁面析出效應(yīng)顯著,溶解氮的擴散過程相比析出過程十分微弱。析出作用下,溶解氮的平均濃度下降了19.5%,含氣率變?yōu)樵瓉淼?.33倍。

本研究建立的氮氣遷移析出模型得到SBLOCA泄壓實驗數(shù)據(jù)的驗證,對于LOCA工況中的氮氣遷移及析出特性的研究具有重要意義。