嚴重事故條件下壓力容器外部臨界熱流密度計算

王凡

國核自儀系統工程有限公司，上海 200241

近年來，非能動設計理念已經開始應用到大型先進壓水堆中，華龍一號(HPR1000)以及國和一號(CAP1400)均采用了非能動設計。發生嚴重事故時，該設計主要是通過實施反應堆壓力容器外部冷卻(external reactor vessel cooling，ERVC)達到熔融物堆內滯留(IVR)，保證放射性物質不釋放到環境中。IVR-ERVC 實施過程中，能否順利通過自然循環帶走多余的衰變熱是保障反應堆安全的核心要務。而自然循環措施成功實施的關鍵在于壓力容器外部的熱通量小于當地的臨界熱流密度(CHF)。

目前，國內外的學者已經對CHF 進行了大量的實驗和理論研究，取得了大量的實驗數據并且建立了相應的經驗關系式。劉瑞蘭等[1]對間隙1.0 mm 的環形狹縫中低質量流量下飽和沸騰強迫循環的換熱特性進行了實驗研究。楊震等[2]在微液層模型的基礎上，利用氣液守恒關系推導出兩相邊界層的初始厚度，選取最大夾帶系數時的空泡份額對下封頭飽和池式沸騰時的CHF 進行理論預測。周磊等[3]通過優化DRYOUT 模型中重要本構關系式組合，適當選取沉積率和夾帶率公式得到適用于矩形窄縫通道的DRTOUT 模型。張亞培等[4]分析了傾角變化的矩形通道內CHF的特點，并建立了適用于矩形窄縫通道的CHF 理論模型，該模型考慮了傾角變化對CHF 的影響。封坤等[5]以氣液兩相逆向對流限制機理為基礎，建立了球形窄縫通道的CHF 理論模型，并進一步分析了系統半徑、熔融物半徑和間隙尺寸等參數對CHF 的影響。陳薇等[6]建立了機理性試驗平臺對加熱表面朝下傾斜矩形通道內的流量波動對CHF 的影響進行試驗，總結出流量波動周期、振幅對CHF 的影響。郭銳[7]通過觀察不同傾角下氣泡的運動特征并結合試驗現象，在氣泡壅塞模型基礎上，建立了IVR-ERVC 條件下的CHF機理模型并對瞬態進程中的CHF 分布進行了計算，以評估嚴重事故瞬態進程中IVR 措施的有效性。Cheung 等[8]利用SBLB 試驗裝置進行了池式沸騰試驗并推導出適用于下封頭外表面的池式沸騰CHF 關系式。

本文在Cheung 等[8]的水動力不穩定性及宏液膜理論模型基礎上，建立了朝下表面矩形通道在過冷低質量流速條件下的CHF 機理模型，該模型將加熱壁面附近的流場分為4 個區域，同時通過能量及質量守恒求解出各層的流動參數，對IVR-ERVC 條件下的臨界熱流密度進行預測分析，重點研究了發生CHF 時，質量流速和流道深度對CHF 的影響，分析了不同質量流速及傾角變化對當地氣泡長度及氣泡速度的影響。

1 數學模型

1.1 基礎模型

壓力容器下封頭發生CHF 時，由于水位高度的原因會造成入口水溫低于當地的飽和溫度，隨著傾角的增加，兩相邊界層的厚度也會逐漸增加。當地發生CHF 時，貼近壁面區域產生的氣泡會有一部分被過冷主流區冷凝。為與實際情況相符合，本文建立的CHF 模型基于以下假設：

1)由于受力和聚合的作用，緊貼加熱壁面處存在較大的細長的變形氣泡；

2)加熱壁面與貼壁氣泡之間存在宏液膜，并且其厚度與過冷度存在關系；

3)兩相邊界層阻礙了主流液體與加熱壁面之間的換熱；

4)當上游流入和主流區擴散至加熱壁面的液體質量小于蒸發消耗量時，發生CHF；

5)宏液膜的厚度隨著當地熱流密度的增加而增加，當發生CHF 時，液膜厚度達到最大。

在加熱壁面附近的細長的大氣泡下面有一個宏液膜。宏液膜由一系列液體薄膜組成，液膜內含有蒸汽柱。在穩定的飽和沸騰條件下，氣柱的質量流速等于泡核沸騰量。該模型將流場分為加熱壁面附近的宏液膜區、大氣泡所在的氣泡層、氣泡中部的兩相區和遠離壁的單相區4 個部分，如圖1 所示。

圖1 IVR-ERVC 條件下通道內流動沸騰示意

宏液膜、氣泡層和兩相區共同構成了兩相邊界層。

宏液膜中的質量流量可表示為

式中：ms為來自上游的液體補充速率，ρl為液體密度，ulm為宏液膜內的液相有效速度，Am為加熱表面和大氣泡之間的橫截面積。

在大氣泡下，宏觀液膜中的蒸發質量率可以表示為

式中：md為質量蒸發速率，qθ為核態沸騰熱通量切向分量，Aw為大氣泡下的加熱壁面積，hfg為汽化潛熱。

結合式(1)和式(2)，可以得到：

式中：

式中：δm為宏液膜的厚度，Lb為細長的氣泡長度，wb為壁面寬度。

將式(4)和式(5)代入式(3)可得

由式(6)可以看出，影響傳熱的獨立變量是宏液膜內的液相流速、宏液膜的厚度和大氣泡長度。

矩形通道中液相和氣相質量的守恒方程為

式中：ulm、ulb和uls分別為宏液膜內液相速度、氣泡中心液相速度、兩相區和單相區的液相速度，αm、αb和 αtw分別為宏液膜內空泡份額、氣泡層內空泡份額和兩相區內空泡份額，δ為兩相邊界層的厚度，H為當地矩形通道高度，G為當地截面平均質量流量。

當大氣泡的浮力超過了周圍流體的阻力和表面張力時，大氣泡開始脫離加熱壁面，這個循環過程在高熱通量下不斷重復，從氣泡開始產生到離開之間的時間稱為懸停時間。為了獲得切向分量，必須首先確定兩相邊界層的厚度。兩相邊界層中的氣相部分由兩相區和單相區擴散進入到宏液膜內的液相受熱蒸發而來，因此，單位時間內擴散入宏液膜內的液相質量有如下關系：

認為 τ時間內，擴散到兩相區熱量qr與兩相邊界層厚度 δ之間存在如下關系：

式中：τ為大氣泡懸停時間；qr為兩相邊界層中總相變熱量的E0倍，0≤E0≤1。

圖2 是宏液膜蒸干模型示意圖，宏液膜內部的波狀氣柱可以近似看作圓柱形，宏液膜內的空泡份額近似等于蒸汽柱面積Ag與氣泡下加熱壁面積的比值，于是：

圖2 宏液膜蒸干模型示意

Cheung 等[9]認為在池沸騰飽和工況發生CHF 時的空泡份額為0.915，所以在本文中認為發生CHF 時，大氣泡區域的空泡份額為0.915，即δm≤y＜δm+Db，單相區內的氣泡很少，空泡份額近似為0。各層平均空泡份額分布如下：

式中y為離壁面的距離。當y=δm+Db時，α=0.915；當y=δ時，α=0。

從宏液膜邊界 δm到變形氣泡中心處，空泡份額單調遞增。Bloch 等[10]通過實驗得出，發生CHF 時，空泡份額與壁面距離線性相關。這里認為，δm+Db到兩相邊界層邊界 δ處，空泡份額的值由0.915 單調遞減為0。

1.2 蒸汽柱面積比例計算

對于理想的無窮大壁面，氣泡從泰勒波的波峰處產生，臨界泰勒波長為最危險波長，臨界泰勒波長可表示為

假設每一個波長平方內的加熱壁面產生一個氣泡，則單位時間內提供給氣泡的蒸汽體積為

氣泡懸停時間和蒸汽體積有著如下關系[11]：

式中 系數11/16 是氣泡移動時所攜帶液體的體積比例。

大氣泡阻礙主流區液體與巨液膜的接觸。因此，液膜的厚度由于蒸發而逐漸減薄。沸騰危機發生時，液膜完全干化的時間剛好等于氣泡懸停的時間。也就是說，氣泡離開時，液膜剛好完全干化。從能量平衡角度，在氣泡懸停的時間內，加熱壁面擴散進來的熱量等于液膜完全蒸發的潛熱。

將 τd和 δm代入式(7)可得：

Zuber 的水動力不穩定性臨界熱流密度模型適用于本實驗，且＜＜1，所以：

1.3 宏液膜厚度計算

高熱通量時的大氣泡與宏液膜已經被很多研究者證實發現[12-15]，宏液膜的厚度隨著熱通量的增加而變薄，大氣泡由于吸收宏液膜不斷蒸發的熱量，體積逐漸增大。當臨界熱通量發生時，在氣泡懸停的時間內，加熱壁面擴散進來的熱量必須等于宏液膜完全蒸發時的潛熱。

在宏液膜上的能量和質量平衡(假設來自加熱器表面的能量完全被用于宏液膜蒸發，導致蒸汽形成)存在以下關系：

式中：ρg為蒸汽密度；ρl為液體密度；vg為蒸汽噴口速度；vl為液體垂直流入加熱面的速度；hfg為汽化潛熱；qnb為當地核態沸騰熱通量，CHF 發生時qnb=qθ。

在τ時間內，宏液膜完全蒸干，則臨界熱流密度條件下的條件為t=τ，δm=0。所以：

為了消除式(10)中蒸汽速度vg，可以將亥姆霍茲不穩定性與宏觀層的初始厚度聯系起來。相對于氣相速度，可以忽略液體速度，并假設宏液膜厚度是四分之一亥姆霍茲不穩定性波長，即

又根據亥姆霍茲不穩定性波長為

將式(10)—式(12)聯立，可以得出：

從式(9)中可以看出，蒸汽柱面積和加熱表面的比值是氣相和液相密度比的函數。

將式(9)代入式(13)，可以得到：

當單相區的部分過冷液體擴散進兩相區，接觸大氣泡的上部時，一部分氣體在大氣泡內凝結，這會增加大氣泡在壁面的懸停時間[16-17]。

式中：Ja為雅各布數，它是液體相變(包括膜狀換熱和沸騰傳熱)時液體過熱的熱量與潛熱的一種度量，；k為經驗系數，這里值取為1。

在式(14)中我們可以看到氣泡懸停時間隨著過冷水平的增加而增加。

1.4 變形氣泡速度及長度計算

宏觀液膜存在于反應堆壓力容器壁外表面上的大氣泡下方。在高熱通量下，大氣泡吸收一些從上游滑移來的氣泡聚結成更大的細長團狀變形氣泡。考慮切向力平衡，液相速度和氣相速度將由阻力和浮力的平衡決定。

式中：Db是蒸氣層的厚度，m；Cd是阻力系數；θ是傾斜角。阻力系數Cd是Db的函數[18]，定義為

聯立式(15)—式(18)可得：

在單相紊流中，液相速度計算公式為[19]

式中：Re為雷諾數，G為質量流速，δm為微液膜厚度，Db為氣泡直徑，μl為液相動力黏度。

實際過程是兩相流動，使用經驗常數0.758 并不準確，因此氣泡中心液相速度計算公式為

式中C是與半徑有關的經驗系數，C=0.15R2。

自然兩相流動過程中，浮力是主要驅動力之一，沸騰產生的氣泡拖曳周圍的液相流動。因此，氣泡層中心處的液相速度最高，氣泡層內空泡份額較高，在氣泡層內的液相速度可以認為等效于氣泡中心處液相速度，然后向兩邊逐漸降低。由于大氣泡緊挨宏液膜，宏液膜內的液相速度ulm與氣泡中心處液相速度ubl呈正相關。

在流動沸騰的過程中，大氣泡在加熱的表面上形成，從壁上分離，然后再次形成，氣泡的形狀和大小會不斷變化。為了便于計算，本文使用氣泡的平均尺寸作為氣泡的脫離直徑，并且假設靠近壁的氣泡的平均直徑與氣泡的分離直徑相同，這樣就可以獲得流體力與氣泡上的表面力之間的平衡關系。計算氣泡直徑采用Cole 等[20]推導出來的關系式：

式中：σ為表面張力，Cpl為液相比熱容，Tsat為當地的飽和溫度。

2 計算結果與分析

2.1 模型計算流程

新建立CHF 理論模型計算流程如圖3 所示。對于給定的管道尺寸、質量流量和入口水溫，首先，計算微液膜的厚度；然后，通過質量守恒方程計算矩形通道內的氣相速度和液相速度；最后，獲得宏液膜蒸發帶走的熱量和單相及兩相區過冷流體擴散入壁面帶走的熱量，兩者相加可得當地的臨界熱流密度。

圖3 CHF 模型計算流程

2.2 模型對比驗證

為驗證本模型對CHF 的預測效果，將模型預測結果和試驗數據進行對比，本文采用Theofanous等[21]有關ULPU-2 400 臺架的試驗數據。圖4 為不同角度計算值與試驗值的比較。由圖4 可以看出，本文建立的CHF 模型隨角度的變化規律與實際試驗相近，這說明建立的CHF 模型對嚴重事故條件下壓力容器外部的沸騰臨界預測是有效的。

圖4 不同角度實驗值與計算值的比較

此外，由圖4 還可以看出，CHF 值隨著角度的增加先是增加較快，到達一定角度(65°～70°)后，變化開始趨于平緩，甚至略有下降。這一方面是由于氣泡上游的累積效應造成壁面傳遞到兩相邊界層和單相區的熱量減少；另一方面，隨著傾角的增加，氣相速度有所增加，同一時間段內帶走的熱量增多。CHF 最終的變化趨勢是由上述2 種因素綜合作用的結果。

由圖5 可以看出，本模型與試驗數據的計算精度在15%以內，模型計算結果與試驗數據的一致性較高，并能較準確地預測CHF 變化趨勢。

圖5 CHF 試驗值與計算值對比

2.3 影響因素分析

圖6 為流道深度0.15 m、壓力容器半徑2 m時，不同質量流速對CHF 的影響。由圖6 可以看到，質量流速的增加對可以顯著提升CHF。

圖6 質量流速對CHF 的影響

圖7 為壓力容器半徑2 m、質量流速600 時，不同流道深度對CHF 的影響。由圖7 可以看出：1)在質量流速一定的情況下，適當增加流道深度(0.01 m)可以提升CHF，但CHF 提升幅度相比質量流速的影響要小。2)流道深度的變化對小角度CHF 影響稍小，對大角度CHF 的影響稍大。這是因為氣泡受浮升力與阻力的共同作用，小角度氣相流速相對大角度較小，大角度通過相變帶走的熱量更多。

圖7 流道深度對CHF 的影響

圖8 和圖9 選取的是壓力容器半徑為2 m、流道深度為0.15 m 條件下的數據。由圖8 和圖9 可以看出：質量流速低時，氣泡的尺寸較大，氣相的速度也較低；相同質量流速時，氣泡的尺寸隨著傾角的增加而變短；氣泡的尺寸隨著質量流速的增加而減小；氣相的速度則隨著質量流速的增加或傾角的增加而增大。

圖8 氣泡長度隨傾角的變化

圖9 氣相速度隨傾角的變化

3 結論

本文在分析現有試驗數據的基礎上，基于水動力不穩定性模型及宏液膜理論，建立了適用于嚴重事故條件下壓力容器外部CHF 理論模型。

1)對于不同區域，必須考慮空泡份額分布的影響，不同區域平均空泡份額在0.013～0.915 變化，進而求解出兩相邊界層厚度。

2)單相紊流中液相速度的求解大體上可應用于兩相流動，但直接應用會造成較大的誤差，相關的經驗系數需要進行修正。

3)對嚴重事故條件下壓力容器外部的CHF進行了計算。發現其他條件相同時，質量流速、傾角和流道深度的增加均對CHF 產生有益的影響。

將計算結果與Dinh 等的試驗數據進行了對比分析。結果表明，模型計算結果與實際試驗數據誤差在15%以內，符合較好。該模型對改善IVR-ERVC 流道結構設計及CHF 理論研究都有一定借鑒意義。