ERVC數(shù)值模擬研究

霍飛鵬，閆大強，李京浩，王 捷

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

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

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

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

1 數(shù)學模型

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

1.1 兩流體模型

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

1）控制方程

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

體積分數(shù)守恒方程：

質量守恒方程：

動量守恒方程：

能量守恒方程：

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

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

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

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

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

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

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

3）相界面質量傳輸（Γk）模型

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

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

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

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

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

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

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

拖曳力MDg是定常流動中兩相之間的主要相間作用力，意義在于當兩相之間流動速度不同時，液相會推動氣相或拖拽氣泡。

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

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

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

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

其中：μtf為液相湍動黏度；σtf為液相湍流施密特數(shù)，通常取1［14］。

1.2 CHF模型

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

1）氣泡擁塞模型

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

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

2）液膜蒸干模型

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

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

2 過冷沸騰三維數(shù)值計算模型驗證

2.1 過冷沸騰實驗

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

圖1 Lee實驗管道示意圖Fig.1 Schematic of Lee's experimental pipe

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

表1 Lee［18］選取的實驗工況Table 1 Experiment conditions of Lee［18］

2.2 計算結果及分析

1）網格無關性

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

2）結果分析

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

圖2 工況1計算結果與實驗值對比Fig.2 Comparison between results of simulation and experiment for condition 1

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

圖3 工況2計算結果與實驗值對比Fig.3 Comparison between results of simulation and experiment for condition 2

圖4 工況3計算結果與實驗值對比Fig.4 Comparison between results of simulation and experiment for condition 3

3 AP1000ERVC三維數(shù)值計算模型驗證

3.1 AP1000ERVC與Lee實驗工況相似性比較

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

3.2 ULPU-V實驗

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

圖5 ULPU-V 實驗裝置示意圖Fig.5 Schematic of ULPU-V facility

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

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

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

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

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

表2 ULPU-V實驗值與計算值比較Table 2 Comparison between results of experiment and simulation for ULPU-V

計算結果顯示71°和83°兩個角度計算值與實驗值的相對偏差均在15%之內，說明三維數(shù)值計算結合CHF模型的可行性。

71°和83°兩個角度的計算值均較實驗值偏小，說明流場結果結合CHF 模型的計算值小于給定的熱流密度邊界條件即實驗值，驗證了前文的設想。

另一方面，71°和83°兩個角度的計算值均較實驗值偏小，為AP1000的技術改進提供了保守性建議。

應當注意的是，已有的CHF 模型都有一定的誤差，且前述設想的充分必要性沒有嚴格證實，需進一步完善。

4 結束語

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

［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.