霍飛鵬，閆大強(qiáng)，李京浩，王 捷

（1.清華大學(xué) 核能與新能源技術(shù)研究院，先進(jìn)核能技術(shù)協(xié)同創(chuàng)新中心，北京 100084；2.國家核電技術(shù)有限公司 北京軟件技術(shù)中心，北京 102209）

AP1000作為從國外引進(jìn)的第三代核電技術(shù)，是我國在建和未來的主要堆型之一。福島核事故以后，許多國家政府以及科研機(jī)構(gòu)對(duì)嚴(yán)重事故下核電站的應(yīng)急響應(yīng)措施予以高度重視。與我國已有的核電站堆型不同，AP1000在嚴(yán)重事故工況下采用堆內(nèi)熔融物滯留（IVR）策略，冷卻水在壓力容器外壁與保溫層之間吸熱汽化形成自然循環(huán)將熔融物的熱量傳出，即壓力容器外部冷卻（ERVC），從而保持壓力容器的完整性，使堆芯熔融物滯留其中?？梢奅RVC是IVR 能否成功實(shí)施的關(guān)鍵。

ERVC本質(zhì)上為壓力容器外壁上的過冷沸騰是否達(dá)到臨界熱流密度（CHF）的問題，目前國內(nèi)外對(duì)ERVC的研究多采用實(shí)驗(yàn)研究，用三維數(shù)值進(jìn)行模擬的較少。1999 年12 月，來自美國和歐洲的9個(gè)機(jī)構(gòu)組成的反應(yīng)堆壓力容器完整性評(píng)估（ARVI）研究小組開始研究關(guān)于如何保持壓力容器完整性的項(xiàng)目［1］，進(jìn)行了EC-FOREVER實(shí)驗(yàn)，并基于實(shí)驗(yàn)數(shù)據(jù)構(gòu)建分析程序模型進(jìn)而進(jìn)行結(jié)構(gòu)分析，自然對(duì)流和下封頭斷裂失效是其研究的重點(diǎn)。美國圣地亞國家實(shí)驗(yàn)室同樣以小比例實(shí)驗(yàn)裝置開展了圓筒沸騰（CYBL）實(shí)驗(yàn)［2］；韓國原子能研究院通過外部自然循環(huán)堆芯冷卻性能實(shí)驗(yàn)（HERMES）［3］，研究兩相自然循環(huán)流動(dòng)和傳熱機(jī)理。美國加州大學(xué)圣巴巴拉分校的全尺寸外部流道（ULPU-V）實(shí)驗(yàn)［4］針對(duì)AP1000進(jìn)行了CHF的實(shí)驗(yàn)，本文選取該實(shí)驗(yàn)進(jìn)行模擬分析。

三維數(shù)值模擬與實(shí)驗(yàn)相比，可獲取整體流場數(shù)據(jù)，能預(yù)測不同參數(shù)工況下流場的變化。本文擬對(duì)低壓過冷沸騰工況構(gòu)建三維流體力學(xué)模型，通過對(duì)Lee等的實(shí)驗(yàn)進(jìn)行模型驗(yàn)證，然后對(duì)AP1000ERVC 進(jìn)行數(shù)值模擬研究，并結(jié)合CHF模型分析ERVC的適用性。

1 數(shù)學(xué)模型

ERVC實(shí)質(zhì)上是兩個(gè)問題的疊加：一是加熱壁面上發(fā)生過冷沸騰的兩相流問題；二是CHF問題。目前用計(jì)算流體力學(xué)（CFD）方法研究兩相流在國內(nèi)外日漸增多，尤其對(duì)簡單流型，如泡狀流，CFD 研究較為成熟。而對(duì)于CHF問題普遍以實(shí)驗(yàn)研究為主，用CFD 研究的較少，且大多關(guān)注棒束等高壓工況，而非AP1000 ERVC 對(duì)應(yīng)的低壓工況。Bergles等［5］提出，CHF 發(fā)生前后流場并無明顯變化，即是否達(dá)到CHF 對(duì)流場影響很小。因此，本文先通過CFD 方法模擬AP1000ERVC 的兩相流流場，再結(jié)合成熟合理的CHF 模型計(jì)算CHF值。

1.1 兩流體模型

CFD 方法模擬過冷沸騰兩相流一般采用歐拉兩流體模型，即對(duì)氣相和液相分別求解質(zhì)量、動(dòng)量和能量方程，對(duì)相界面質(zhì)量、動(dòng)量和能量交換構(gòu)筑模型來封閉整個(gè)控制方程組。

1）控制方程

兩流體模型的控制方程是基于三維流體力學(xué)方程分別對(duì)每項(xiàng)建立方程組并進(jìn)行系綜平均，或時(shí)間、空間平均。經(jīng)過平均后，舍棄相界面上階躍變化的物理信息。兩相守恒方程組和相界面的3個(gè)守恒方程構(gòu)成了兩流體模型的控制方程組。文獻(xiàn)［6］對(duì)這兩個(gè)方程組進(jìn)行了簡化，并認(rèn)為相界面是由液相和氣相構(gòu)成的厚度很薄的控制體。

體積分?jǐn)?shù)守恒方程：

質(zhì)量守恒方程：

動(dòng)量守恒方程：

能量守恒方程：

其中：α為體積分?jǐn)?shù)，下標(biāo)f和g分別代表液相和氣相，k代表兩相中的任一相；Γk為k 相因相變產(chǎn)生的質(zhì)量源；τk表示k 相切應(yīng)力張量；為k 相湍流應(yīng)力；vki表示相界面上k 相的速度；τki為相界面上k 相切應(yīng)力張量；pki為相界面上k 相的壓力；Mk為k 相界面上的動(dòng)量源；Mki為相間作用力；ek為內(nèi)能；qk為熱流率；為湍流產(chǎn)生的熱流率；Tk為壓力張量和切應(yīng)力張量構(gòu)成的應(yīng)力張量，Tk＝-pkI＋τk，I為單位張量；Ek為交界面上k 相的能量源；iki為相界面上k 相的比焓；qki為相界面作為控制體與相鄰的相傳熱的熱流密度；為相界面上的湍流作用做的功。

對(duì)于過冷沸騰兩相流問題，流型多以泡狀流為主，可假設(shè)相界面上的壓力、切向應(yīng)力與液相一致，相界面的比焓與同側(cè)相的一致，相界面的速度與氣泡速度一致。WTki近似為0。

簡化后，式（1）～（6）中除Γk、Mki、qki外的其他項(xiàng)均可由NS 方程的未知量，即速度、壓力、溫度等流場參數(shù)表示。因此需引入輔助模型將Γk、Mki、qki表示成流場參數(shù)的函數(shù)，從而封閉整個(gè)方程組。

國內(nèi)外對(duì)兩相流中相界面的質(zhì)量傳輸、相間作用力、相間傳熱、氣泡直徑等研究較多，本課題組前期曾對(duì)相關(guān)模型進(jìn)行了總結(jié)與對(duì)比，本文僅列出最優(yōu)化的模型。

2）相界面熱流密度（qki）模型

在過冷沸騰中，近似認(rèn)為氣相不存在過熱，即氣相溫度等于飽和溫度。則相界面熱流密度由下式表示：

其中：Tsat為飽和溫度；hi為相間換熱系數(shù)，通常用無量綱參數(shù)Nu 表示。Ranz等［7］給出了Nu和氣泡雷諾數(shù)ReB及液相普朗特?cái)?shù)Pr 的關(guān)系式：

3）相界面質(zhì)量傳輸（Γk）模型

相界面質(zhì)量傳輸是伴隨著蒸發(fā)或冷凝的能量傳輸過程進(jìn)行的，根據(jù)式（12）中的hi，Γk用下式表示：其中：hfg為蒸發(fā)潛熱；ai為相界面濃度。

相界面濃度ai是關(guān)于空泡份額、氣泡直徑、流道等效直徑的函數(shù)關(guān)系。Ishii等［6］給出了如下模型：

其中：αg為氣相體積分?jǐn)?shù)；αgs為小氣泡區(qū)域的體積分?jǐn)?shù)；Deq為等效直徑；dB為氣泡直徑。

在沸騰流動(dòng)中，氣泡的尺寸分布范圍很大。Kurul［8］認(rèn)為氣泡直徑dB與當(dāng)?shù)剡^冷度ΔTsub（ΔTsub＝Tsat-Tf）呈正比。Koncar等［9］在Kurul模型的基礎(chǔ)上給出了低壓情況下的模型：

4）相間作用力（Mki）模型

在兩相流理論研究中，常將相間作用力拆分成5個(gè)部分，分別建立模型：

其中：MDg為拖曳力；MVg為虛擬質(zhì)量力；MLg為升力；MWg為壁面潤滑力；MTg為湍流擴(kuò)散力。

拖曳力MDg是定常流動(dòng)中兩相之間的主要相間作用力，意義在于當(dāng)兩相之間流動(dòng)速度不同時(shí)，液相會(huì)推動(dòng)氣相或拖拽氣泡。

其中：vr＝vg-vf；Ag為氣泡在流動(dòng)方向投影的橫截面積；VB為氣泡的體積；CD為拖曳力系數(shù)，采用Ishii等［10］給出的模型計(jì)算。

升力MLg是氣泡在垂直于相對(duì)運(yùn)動(dòng)速度方向的速度梯度產(chǎn)生的力。升力的作用是將氣泡推向流速較低的區(qū)域。Drew 等［11］提出的升力表達(dá)式如下：

其中，CL采用Tomiyama［12］的升力系數(shù)關(guān)系式計(jì)算。

其中：yW為與壁面的垂直距離；nW為壁面單位法向量；CW1和CW2為壁面潤滑力系數(shù)，在處理低壓工況時(shí)，Koncar［9］采用CW1＝-0.04，CW2＝0.08。

其中：μtf為液相湍動(dòng)黏度；σtf為液相湍流施密特?cái)?shù)，通常取1［14］。

1.2 CHF模型

CHF按照流場情況和觸發(fā)機(jī)理可分為兩類：干涸（Dry-out）和偏離泡核沸騰（DNB）。在過冷流動(dòng)和低空泡份額區(qū)域，發(fā)生的CHF 屬于DNB。本文關(guān)注的AP1000 ERVC 中的CHF屬于DNB類型。目前DNB 的物理機(jī)理尚未得到清楚的理解，得到較多認(rèn)可的有兩種理論模型：氣泡擁塞模型和液膜蒸干模型。

1）氣泡擁塞模型

Weisman等［15］提出在泡核沸騰起始點(diǎn)（ONB）后，加熱壁面上附著一層氣泡層。氣泡從上游進(jìn)入從下游流出。當(dāng)氣泡區(qū)中的氣泡密集到一定程度從而阻礙氣泡區(qū)與主流區(qū)之間的焓傳遞時(shí)，即認(rèn)為發(fā)生了CHF。CHF模型如下：

其中：G 為 質(zhì) 量 流 速；x1、x2分 別 為 主 流 區(qū) 和 氣泡區(qū)的質(zhì)量含汽率；hf為飽和液相焓；hl為主流區(qū)的液相焓；hld為氣泡附著點(diǎn)處的液相焓；ib為氣泡層和主流層接觸面上的湍流強(qiáng)度；ψ 為壁面附近的速度波動(dòng)系數(shù)。

2）液膜蒸干模型

Fiori等［16］認(rèn)為，在低質(zhì)量流率、低過冷度下，在壁面附近會(huì)出現(xiàn)大氣彈，與壁面間有一層液膜。當(dāng)液膜蒸干時(shí)發(fā)生CHF。Lee等［17］總結(jié)出的液膜蒸干CHF模型如下：

其中：Gm為進(jìn)入液體薄膜的液體流率；δm為液體薄膜的厚度；hm為液體薄膜的焓；Lm為亥姆霍茲臨界波長；a1為經(jīng)驗(yàn)系數(shù)。

2 過冷沸騰三維數(shù)值計(jì)算模型驗(yàn)證

2.1 過冷沸騰實(shí)驗(yàn)

Lee［18］對(duì)常壓環(huán)形管道強(qiáng)制對(duì)流過冷沸騰的體積份額、兩相速度徑向分布等進(jìn)行了實(shí)驗(yàn)研究。內(nèi)管外徑19 mm，外管內(nèi)徑37.5 mm，全長為2.376m，內(nèi)管加熱，其中加熱段長度為1.67m，如圖1所示。圖1中，RO為外管內(nèi)側(cè)半徑，Ri為內(nèi)管外側(cè)半徑。加熱段的前后管道長度分別為0.28m 和0.61 m，局部參數(shù)測量位置位于距離加熱段出口0.06m 處。

圖1 Lee實(shí)驗(yàn)管道示意圖Fig.1 Schematic of Lee's experimental pipe

根據(jù)AP1000ERVC 的 物 理 過 程［4］，壓 力容器外部的流道壓力約在0.13～0.16MPa之間，過冷度為13K 左右，因此本文選取表1中的3個(gè)實(shí)驗(yàn)工況進(jìn)行模擬計(jì)算。

表1 Lee［18］選取的實(shí)驗(yàn)工況Table 1 Experiment conditions of Lee［18］

2.2 計(jì)算結(jié)果及分析

1）網(wǎng)格無關(guān)性

由于實(shí)驗(yàn)管道是環(huán)形管道，中心對(duì)稱，所以計(jì)算域選取1／6的扇形棱柱。網(wǎng)格敏感性分析結(jié)果顯示，軸向網(wǎng)格和周向網(wǎng)格的精度均對(duì)計(jì)算結(jié)果影響不大，徑向網(wǎng)格的精度對(duì)計(jì)算結(jié)果影響較大。最終采用14（徑向）×12（周向）×100（軸向）的均分網(wǎng)格。

2）結(jié)果分析

計(jì)算采用ANSYS CFX 14.5，三維數(shù)值計(jì)算模型采用第1節(jié)中列舉的模型，計(jì)算結(jié)果與實(shí)驗(yàn)值的對(duì)比如圖2～4所示。其中，橫坐標(biāo)徑向位置是指距環(huán)形管軸線的距離，內(nèi)管外壁面徑向距離為0.009 5 m，外管內(nèi)壁面徑向距離為0.018 75m。

圖2 工況1計(jì)算結(jié)果與實(shí)驗(yàn)值對(duì)比Fig.2 Comparison between results of simulation and experiment for condition 1

由圖2～4可見，空泡份額分布趨勢及峰值的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果基本符合。尤其是在靠近加熱壁面附近區(qū)域，計(jì)算值與實(shí)驗(yàn)值符合最好。而ERVC中CHF與壁面附近的流動(dòng)情況關(guān)系最為密切，壁面附近區(qū)域計(jì)算結(jié)果是否準(zhǔn)確決定了CHF 值計(jì)算的準(zhǔn)確性。Lee實(shí)驗(yàn)計(jì)算結(jié)果表明，使用三維數(shù)值計(jì)算方法可得到較為準(zhǔn)確的近壁面流場。

圖3 工況2計(jì)算結(jié)果與實(shí)驗(yàn)值對(duì)比Fig.3 Comparison between results of simulation and experiment for condition 2

圖4 工況3計(jì)算結(jié)果與實(shí)驗(yàn)值對(duì)比Fig.4 Comparison between results of simulation and experiment for condition 3

3 AP1000ERVC三維數(shù)值計(jì)算模型驗(yàn)證

3.1 AP1000ERVC與Lee實(shí)驗(yàn)工況相似性比較

AP1000ERVC的工況參數(shù)和Lee實(shí)驗(yàn)工況的相同點(diǎn)是參數(shù)較為接近，皆為常壓，加熱壁面上發(fā)生過冷沸騰，且入口過冷度較相似。不同點(diǎn)是AP1000ERVC 為自然循環(huán)，Lee的實(shí)驗(yàn)為強(qiáng)迫對(duì)流。由于AP1000ERVC 重點(diǎn)在CHF發(fā)生時(shí)的工況，且ULPU-V 實(shí)驗(yàn)采取的是緩慢升溫的手段，可看作準(zhǔn)穩(wěn)態(tài)過程，將CHF發(fā)生時(shí)的流量作為穩(wěn)態(tài)流量計(jì)算。因此將前文總結(jié)的三維數(shù)值計(jì)算模型用在ERVC工況計(jì)算中。

3.2 ULPU-V實(shí)驗(yàn)

AP1000ERVC 的實(shí)驗(yàn)數(shù)據(jù)較少，本文采用ULPU-V 進(jìn)行研究。ULPU-V 是針對(duì)AP1000采用高度比為1比1的二維切片型實(shí)驗(yàn)裝置［19］，如圖5所示。用電加熱塊模擬熔池?zé)嵩?，加熱功率模擬AP1000嚴(yán)重事故熔池?zé)崃髅芏确植?。冷凝器壓力?01 325Pa，溫度為373K，管道均為絕熱，加熱表面為90°圓弧形切片，半徑為2.006m，流道寬度為0.15m。從冷凝器到進(jìn)口擋板的高度為6.136m。

圖5 ULPU-V 實(shí)驗(yàn)裝置示意圖Fig.5 Schematic of ULPU-V facility

將進(jìn)口處記為0°，加熱流道出口處記為90°。由于實(shí)驗(yàn)考慮了不同流道寬度、不同工質(zhì)對(duì)結(jié)果的影響，真正與AP1000相符的實(shí)驗(yàn)數(shù)據(jù)只有71°和83°兩個(gè)角度的CHF 值，且公開文獻(xiàn)中無局部流場參數(shù)。

3.3 三維數(shù)值計(jì)算結(jié)合CHF模型

計(jì)算模型采用第2節(jié)驗(yàn)證過的三維數(shù)值計(jì)算模型，傳熱邊界條件采用實(shí)驗(yàn)中發(fā)生CHF時(shí)加熱塊的熱流密度。

CHF發(fā)生時(shí)，流場并未發(fā)生明顯轉(zhuǎn)變，無法通過流型的過渡來判斷CHF 是否達(dá)到。而現(xiàn)有的關(guān)于DNB的CHF模型基本都表達(dá)為流場變量的函數(shù)。因此設(shè)想用三維數(shù)值模擬計(jì)算得到的流場數(shù)據(jù)結(jié)合CHF 模型推導(dǎo)出在此時(shí)流場變量下對(duì)應(yīng)的壁面位置的CHF 值，與三維數(shù)值計(jì)算的邊界條件中的熱流密度比較，如果反推出的CHF 值小于事先給定的熱流密度值，則認(rèn)為此邊界條件下會(huì)發(fā)生CHF。

由于AP1000ERVC兩相流流型以泡狀流為主，氣泡擁塞模型較液膜蒸干模型更為適合，因此本文選用Weisman等［15］提出的氣泡擁塞CHF模型，即式（23）。計(jì)算結(jié)果與實(shí)驗(yàn)值的對(duì)比列于表2。

表2 ULPU-V實(shí)驗(yàn)值與計(jì)算值比較Table 2 Comparison between results of experiment and simulation for ULPU-V

計(jì)算結(jié)果顯示71°和83°兩個(gè)角度計(jì)算值與實(shí)驗(yàn)值的相對(duì)偏差均在15%之內(nèi)，說明三維數(shù)值計(jì)算結(jié)合CHF模型的可行性。

71°和83°兩個(gè)角度的計(jì)算值均較實(shí)驗(yàn)值偏小，說明流場結(jié)果結(jié)合CHF 模型的計(jì)算值小于給定的熱流密度邊界條件即實(shí)驗(yàn)值，驗(yàn)證了前文的設(shè)想。

另一方面，71°和83°兩個(gè)角度的計(jì)算值均較實(shí)驗(yàn)值偏小，為AP1000的技術(shù)改進(jìn)提供了保守性建議。

應(yīng)當(dāng)注意的是，已有的CHF 模型都有一定的誤差，且前述設(shè)想的充分必要性沒有嚴(yán)格證實(shí)，需進(jìn)一步完善。

4 結(jié)束語

ERVC能否避免CHF 發(fā)生是AP1000嚴(yán)重事故響應(yīng)策略IVR 的關(guān)鍵。本文提出一種全新的方法，使用驗(yàn)證過的三維數(shù)值模擬方法得出流場結(jié)果，再結(jié)合已有的CHF 模型獲得流場結(jié)果對(duì)應(yīng)的CHF 值，與所給的熱流密度邊界條件比較得出此邊界條件下是否會(huì)發(fā)生CHF。計(jì)算結(jié)果表明這種方法有一定可行性，但需較為合適且準(zhǔn)確的CHF 模型，未來仍需詳細(xì)論證此設(shè)想的充分必要性。

［1］ SEHGAL B R，THEERTHAN A，GIRI A.Assessment of reactor vessel integrity（ARVI）［J］.Nuclear Engineering and Design，2003，221：23-53.

［2］ CHU T Y，BENTZ J H，SLEZASEDAG S E.Ex-vessel boiling experiments：Laboratory-and reactor-scale testing of the flooded cavity concept for in-vessel core retention，Part Ⅱ：Reactorscale boiling experiments of the flooded cavity concept for in-vessel core retention［J］.Nuclear Engineering and Design，1997，169：89-99.

［3］ PARK R J，HA K S，KIM S B，et al.Twophase natural circulation flow of air and water in a reactor cavity model under an external vessel cooling during a severe accident［J］.Nuclear Engineering and Design，2006，236：2 424-2 430.

［4］ THEOFANOUS T G，LIU C，ADDITION S，et al.In-vessel coolability and retention of a core melt［J］.Nuclear Engineering and Design，1997，169：1-48.

［5］ BERGLES A E，LOPINA R F，F(xiàn)IORI M P.Critical-heat-flux and flow-pattern observation for low-pressure water flowing in tubes［J］.J Heat Trans，1967，89：69-74.

［6］ ISHII M，HIBIKI T.Thermo-fluid dynamics of two-phase flow ［M］.New York：Springer，2011.

［7］ RANZ W E，MARSHALL W R.Evaporation from drops［J］.Chem Eng Prog，1952，48（3）：141-146.

［8］ KURUL N.Multidimensional effects in twophase flow including phase change［D］.US：Rensselaer Polytechnic Institute，1990.

［9］ KONCAR B，KLJENAK I，MAVKO B.Modeling of local two-phase flow parameters in upward subcooled flow boiling at low pressure［J］.Int J Heat Mass Transfer，2004，47：1 499-1 513.

［10］ISHII M，ZUBER N.Drag coefficient and rela-tive velocity in bubbly，droplet or particulate flows［J］.AIChE J，1979，25：843-855.

［11］DREW D A，LAHEY R T，Jr.The virtual mass and lift force on a sphere in rotating and straining flow［J］.Int J Multiphase Flow，1987，13：113-121.

［12］TOMIYAMA A.Struggle with computational bubble dynamics［C］∥Proceedings of Third Int Conf Multiphase Flow （CD-ROM）.Lyon：［s.n.］，1998.

［13］ANTAL S P，LAHEY R T，Jr，F(xiàn)LAHERTY J E.Ana1ysis of phase distribution in fully developed laminar bubbly two-phase flow［J］.Int J Multiphase Flow，1991，17：635-652.

［14］BURNS A D B，F(xiàn)RANK T，HAMILL I，et al.The favre averaged drag model for turbulent dispersion in eulerian multi-phase flows［C］∥Proceedings of 5th International Conference on Multiphase Flow.Yokohama：［s.n.］，2004.

［15］WEISMAN J，PEI B S.Prediction of critical heat flux in flow boiling at low qualities［J］.Int J Heat Mass Transfer，1983，26：1 463-1 477.

［16］FIORI M P，BERGLES A F.Model of critical heat flux in subcooled flow boiling［C］∥Heat Transfer 1970.Amsterdam：Elsevier Publishing Company，1970.

［17］LEE C H，MUDAWAR I.A mechanistic critical heat flux model for subcooled flow boiling based on local bulk flow conditions［J］.Int J Multiphase Flow，1988，14：711-728.

［18］LEE T H.Local flow characteristics of subcooled boiling flow of water in a vertical concentric annulus［J］.International Journal of Multiphase Flow，2002，28：1 351-1 368.

［19］DINH T，TU J P，SALMASSI T，et al.Limits of coolability in the AP1000-related ULPU-2400 configuration Ⅴfacility［R］.US：［s.n.］，2003.