呂玉鳳，趙民富，杜開文

(中國原子能科學(xué)研究院 反應(yīng)堆工程技術(shù)研究部，北京 102413)

冷卻劑喪失事故(LOCA)作為反應(yīng)堆的設(shè)計基準(zhǔn)事故之一，是反應(yīng)堆安全分析的重要方面。反應(yīng)堆運行在高溫、高壓條件下，發(fā)生LOCA時，冷卻水從破口噴出，回路迅速降壓，回路的高溫水急劇汽化，形成氣液兩相流，破口處將會發(fā)生臨界流動[1]。破口處的質(zhì)量流量決定著反應(yīng)堆冷卻水系統(tǒng)的裝量和能量喪失速率，進而強烈影響著系統(tǒng)壓力、堆芯冷卻劑溫度和燃料包殼溫度的變化規(guī)律，從而對事故后果起重要作用。在先進壓水堆中，自動卸壓系統(tǒng)(ADS)是一項重要的非能動安全技術(shù)[2]。ADS的流體排放也涉及臨界流動現(xiàn)象，此時的流動特性控制著事故的進程和安全系統(tǒng)的動作，與事故的緩解密切相關(guān)[3]。因此，清楚認(rèn)識臨界流動現(xiàn)象并建立合適的臨界流量計算方法，是核反應(yīng)堆安全分析程序發(fā)展中的重要課題。

迄今為止，有關(guān)兩相臨界流動的研究已進行了大量工作，并提出了很多分析模型和計算關(guān)系式，Abdollahian等[4]、Saha[5]對這方面工作做了非常詳盡的評述。文獻[6]中對現(xiàn)有臨界流模型進行評價，結(jié)果表明熱力學(xué)平衡模型(包括均勻平衡模型、滑移平衡模型)在預(yù)測滯止?fàn)顟B(tài)為飽和或過冷液體流過短管及噴嘴時，產(chǎn)生極大的誤差，有些工況誤差高達100%。對于短管及噴嘴的臨界流動，兩相膨脹非?？?，以致相間不能達到平衡狀態(tài)，汽化核心在過熱液體及通道的壁面上形成并長大。因此，需發(fā)展既反映相間速度差異(即力學(xué)非均勻)，又反映溫度差異(即熱力學(xué)非平衡)的數(shù)值模型來預(yù)測臨界流量。

熱力學(xué)非平衡模型可分為非平衡經(jīng)驗?zāi)Ｐ?如Henry-Fauske模型[7])和反映熱力學(xué)非平衡物理機制的兩流體模型。Henry-Fauske模型引入經(jīng)驗系數(shù)N來修正熱力學(xué)平衡含汽率xe而得到其真實含汽率，這相比熱力學(xué)平衡模型有了較大改善。但該系數(shù)僅是對實驗數(shù)據(jù)的關(guān)聯(lián)而未反映非平衡機理。Abdollahian等[4]采用Marviken數(shù)據(jù)對Henry-Fauske模型進行評價時，為使模型計算結(jié)果匹配實驗數(shù)據(jù)，對N進行了修改，N的取值范圍為7.1xe～100xe。兩流體模型試圖從臨界流動的機理出發(fā)，分別對氣液兩相的連續(xù)、動量、能量方程進行描述，并補充一些基本結(jié)構(gòu)關(guān)系式從而構(gòu)成封閉方程組，如兩相間的相互作用力、相間傳熱特性、相界面動量傳遞等。兩流體模型可反映相間的熱力學(xué)非平衡性，非平衡的一個關(guān)鍵問題是汽化(閃蒸)現(xiàn)象，包括閃蒸起始點、汽泡初始直徑、汽核密度及汽泡的增長規(guī)律。由于兩相流微觀機理的復(fù)雜性，目前對這些物理量的研究尚不充分。在Richter[8]、Dobran[9]、Schwellnus等[10]提出的考慮熱力學(xué)非平衡效應(yīng)的兩流體模型中，均部分地假設(shè)了未知參數(shù)數(shù)值來“符合”實驗數(shù)據(jù)。

本文從氣液兩相的連續(xù)方程、動量方程及能量方程出發(fā)，建立兩流體六方程熱力學(xué)非平衡模型，對現(xiàn)有計算閃蒸起始點位置和過熱度的公式進行比較分析，最終選取較合理的計算公式。在泡狀流區(qū)域，將汽泡增長方程與基本方程聯(lián)立求解，可較準(zhǔn)確地反映汽泡增長規(guī)律。

1 兩流體模型

本模型由6個方程構(gòu)成：氣相和液相的質(zhì)量守恒方程、氣相和液相的動量守恒方程、兩相混合物的總能量守恒方程、氣相的能量守恒方程。另外在泡狀流區(qū)補充1個汽泡增長方程。主要假設(shè)有：1) 一維絕熱流動；2) 在兩相流動中液相是飽和或過熱的，氣相是飽和的，溫度等于當(dāng)?shù)貕毫ο碌娘柡蜏囟龋?) 忽略兩相間的黏性應(yīng)力及黏性耗散。

1.1 基本控制方程

Dobran[9]給出了一維、穩(wěn)態(tài)的兩相流動守恒方程。

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

(1)

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

(2)

液相動量守恒方程：

Fwl+Flg-Fli-(1-α)ρlgsinθ

(3)

氣相動量守恒方程：

-Fwg-Flg-Fgi-αρggsinθ

(4)

總能量守恒方程：

(5)

氣相能量守恒方程：

(6)

式中：A為通道截面積；α為空泡份額；cpl為液相比定壓熱容；Fwl、Fwg分別為通道壁面與液相、氣相摩擦引起的壓力損失；Flg為相界面相互作用力；Fli、Fgi為相界面上液相、氣相由于質(zhì)量交換引起的動量傳遞；g為重力加速度；G為質(zhì)量流速；hl、hg分別為液相和氣相的比焓；p為壓力；Tl為液相溫度；ul、ug分別為液相和氣相的速度；qi為單位體積內(nèi)界面?zhèn)鬟f的熱量；θ為通道軸向與水平方向夾角；ρl、ρg分別為液相、氣相的密度；x為含汽率；z為通道軸向距離。

上述方程涉及到6個變量：α、x、ug、ul、p、Tl。在泡狀流區(qū)域(0<α≤0.3)，引入第7個變量，汽泡直徑db。假設(shè)平均汽核密度為N，則有：

(7)

將式(7)與式(1)聯(lián)立，可得：

(8)

當(dāng)處于泡狀區(qū)時，將式(1)～(6)與式(8)聯(lián)立，構(gòu)成七方程組，db作為第7個變量增加到方程組中，db隨x和α的變化而變化；當(dāng)不在泡狀區(qū)時，方程組直接由式(1)～(6)求解，變?yōu)?個守恒方程。

在Richter[8]、Dobran[9]提出的兩流體模型中，假設(shè)汽核密度N為一常數(shù)。而在本模型中將式(1)～(6)與式(8)聯(lián)立求解，N與db、x及α相互關(guān)聯(lián)，是不斷變化的，因而能更好、更真實地描述泡化過程。

為求解上述變量，需給出Fwl、Fwg、Flg、Fli、Fgi、qi結(jié)構(gòu)關(guān)系式。狀態(tài)參數(shù)ρg、ρl、hg、hl、cpl、?ρg/?p、?ρl/?p、?hg/?p可由水物性計算軟件得到。

1.2 結(jié)構(gòu)關(guān)系式

1)Fwl、Fwg

Fwl通過兩相摩擦壓降倍率與液相單獨流過通道時的壓降梯度聯(lián)系，即：

(9)

(10)

式中，fl0為單相摩擦系數(shù)，由下式確定：

(11)

式中：Re為雷諾數(shù)；ε/D為壁面粗糙度，由于氣相和壁面不是直接接觸，故可忽略氣相和壁面的摩擦力，即Fwg=0。

2)Flg

Flg主要由兩方面組成：相界面上的拖拽力FD和兩相間相對加速引起的虛擬質(zhì)量力FVM[12]。

Flg=FVM+FD

(12)

FVM=CVMρlαVM

(13)

(14)

(15)

λ=2(1-α)

(16)

式中：CVM為虛擬質(zhì)量力系數(shù)；λ為系數(shù)。

在泡狀流區(qū)域，F(xiàn)D為：

(17)

(CD)1-α=CD(1-α)-4.7

(18)

CD=0.44Reb>1 000

(19)

式中，CD為汽泡與汽泡之間的相互作用力系數(shù)。

對于環(huán)狀流區(qū)域，F(xiàn)D為：

(20)

式中，界面拖拽力系數(shù)Cfi為：

Cfi=0.005[1+75(1-α)]

(21)

CD和Cfi之間的關(guān)系為：

(22)

對于彈狀流區(qū)，兩個轉(zhuǎn)換點的數(shù)值變化很大，采用指數(shù)插值形式。

3)Fli、Fgi

Fli、Fgi如下：

(23)

(24)

式中，系數(shù)η取0.5。

4) 相界面間的傳熱qi

單位體積內(nèi)的相界面?zhèn)鳠岜磉_式為：

qi=aihi(Tl-Tg)

(25)

式中：ai為單位體積的氣液兩相相間面積；Tl為液相溫度；hi為相間傳熱系數(shù)。

對流傳熱系數(shù)在不同的流型中不同，在泡狀區(qū)由下式?jīng)Q定：

(26)

式中：Nu為努塞爾數(shù)；kl為液相導(dǎo)熱系數(shù)；Pr為普朗特數(shù)。

對于環(huán)狀流區(qū)域，相界面?zhèn)鳠嵯禂?shù)表達式為：

(27)

對于彈狀流的傳熱系數(shù)，可通過線性插值的方法得到。

1.3 邊界條件

Fraser等[13]的研究表明，閃蒸點的位置和過熱度對臨界流量有顯著的影響。對于面積不變的噴管，發(fā)生閃蒸起始點的位置并不是固定的，而是與上游滯止參數(shù)、通道結(jié)構(gòu)、壁面條件、凝結(jié)氣體等因素有關(guān)。

在兩流體計算模型中[8,10]采用的計算閃蒸起始點過熱度的公式為：

(28)

式中：Tsat為飽和溫度；σ為該溫度下的表面張力；vlg為氣、液兩相比容之差；slg為氣、液兩相比熵之差；文獻中一般設(shè)置閃蒸起始點初始直徑db為25 μm。

但按照該公式計算的閃蒸起始點過熱度一般很低。文獻[6]根據(jù)試驗數(shù)據(jù)對現(xiàn)有閃蒸點過熱度的公式進行評價表明，Alamgir & Lienhard公式[14](式(29))計算的過熱度與短噴管的試驗數(shù)據(jù)符合較好。并且，該式可體現(xiàn)上游滯止參數(shù)、卸壓速度等對過熱度的影響，更能反映真實情況。因此，本文首先采用Alamgir & Lienhard公式計算閃蒸起始點的壓力：

pFi=psat(T0)-

(29)

式中：pFi為閃蒸起始點壓力；T0為初始溫度；Tr=T/Tc，Tc為臨界溫度，Tc=647.3 K；k為波爾茲曼常數(shù)，k=1.380 54×10-23；Σ′為卸壓速度。

求得pFi之后，可按T0-Tsat(pFi)計算閃蒸起始點處的過熱度。然后由伯努利方程及質(zhì)量守恒方程，求得閃蒸起始點的位置zFi：

(30)

1.4 求解過程

根據(jù)上述的控制方程、狀態(tài)方程和基本結(jié)構(gòu)關(guān)系式，經(jīng)整理可得如下矩陣方程：

(31)

2 模型驗證

2.1 與中國原子能科學(xué)研究院(CIAE)試驗比較

CIAE開展了穩(wěn)態(tài)臨界流動試驗[15]，試驗段直徑為2 mm，長徑比(L/D)分別為1、3、6、12、20。試驗壓力范圍為3～16 MPa，入口溫度范圍為30～345 ℃，入口含汽率為0～0.5。圖1a為L/D=1的噴管在壓力為10～16 MPa、過冷度為0～40 ℃范圍內(nèi)質(zhì)量流速測量值與模型計算值的比較。對于大部分實驗工況，程序預(yù)測結(jié)果偏低，這是因為該區(qū)域熱力學(xué)非平衡效應(yīng)非常嚴(yán)重，最大相對偏差約為10%。圖1b為L/D=12噴管的模型計算值與試驗結(jié)果比較，同樣模型預(yù)測質(zhì)量流速與試驗結(jié)果符合較好。

圖2示出了兩相流基本參量沿管長的變化。試驗工況為：初始壓力為10 MPa，初始過冷度為20.1 ℃，L/D=12?？梢?，沿管長方向α、x、氣相和液相速度是不斷增加的，在接近通道出口處，α及兩相速度急劇增加。

圖1 模型計算結(jié)果與CIAE試驗值的比較Fig.1 Comparison of model prediction and CIAE experiment results

圖2 兩相流基本參量沿管長的變化Fig.2 Variation of two-phase flow parameter along tube

'

2.2 與Fauske實驗比較

1965年，F(xiàn)auske[16]采用銳角入口的噴管試驗段開展了臨界流動實驗，試驗段直徑為6.35 mm，滯止壓力范圍為0.5～15 MPa，工質(zhì)為飽和水或過冷水。圖3a為Fauske采用的銳角入口試驗段在不同L/D下實驗結(jié)果與模型計算結(jié)果的比較?？梢姵绦蝾A(yù)測的質(zhì)量流速與實驗值符合很好，最大相對偏差不超過10%。圖3b為不同滯止壓力條件下Fauske實驗數(shù)據(jù)與模型計算結(jié)果的比較，此處Fauske在實驗中采用圓角入口的試驗段，試驗段直徑為6.35 mm，長度為254 mm，即L/D=40，實驗工質(zhì)為飽和水。可見，在不同的滯止壓力下模型計算結(jié)果與實驗結(jié)果是比較符合的。

圖3 模型計算值與Fauske實驗值的比較Fig.3 Comparison of model prediction and Fauske experiment results

2.3 與Celata等實驗的比較

圖4為本模型計算的壓力沿通道分布與Celata等[17]的實驗數(shù)據(jù)比較。試驗采用的噴管直徑為0.012 5 m，L/D=97，圓角入口。圖4a的實驗工況為：滯止壓力2.0 MPa，過冷度ΔTsub=22.7 ℃，實驗測量得到的臨界質(zhì)量流速為23 000 kg/(m2·s)，程序計算值為23 654 kg/(m2·s)，二者相對偏差為2.8%。程序計算的壓力分布與實驗結(jié)果符合很好。圖4b的實驗工況為：滯止壓力2.23 MPa，滯止?fàn)顟B(tài)為飽和水。實驗測量得到的臨界質(zhì)量流速與程序計算值相對偏差為5.2%。程序計算的壓力較測量值略高，可能是因為管道入口摩擦壓降和兩維效應(yīng)對汽泡生長有較大影響，與程序中關(guān)于汽泡生長的模型描述不一致導(dǎo)致的。

圖4 模型計算壓力分布與Celata等實驗值的比較Fig.4 Comparison of pressure distributions between prediction and Celata et al. experiment

3 結(jié)論

本文從兩相臨界流動的熱力學(xué)非平衡機制出發(fā)，建立了兩流體臨界流模型。模型做出的主要改進如下。

1) 引入了經(jīng)過實驗驗證的計算閃蒸起始點過熱度和閃蒸點位置的計算公式，可體現(xiàn)上游滯止參數(shù)、卸壓速度等對過熱度的影響，更能反映真實情況。

2) 在泡狀流區(qū)域，將汽泡增長方程與控制方程聯(lián)立求解，得到的汽泡直徑與壓力、含汽率和空泡份額有關(guān)，能比較準(zhǔn)確地反映汽泡的增長規(guī)律，更接近實際情況。

該模型在很寬的試驗范圍內(nèi)成功地預(yù)測了臨界質(zhì)量流速和壓力分布，驗證了模型的通用性。