朱光昱，郭　勇，元一單，劉宇生，李 煒

氣泡微細化沸騰觸發(fā)溫度數(shù)值模擬研究

朱光昱1,2，郭勇1，元一單1，劉宇生2，李 煒1

（1. 中國核電工程有限公司核電安全研究中心，北京 100840；2. 生態(tài)環(huán)境部核與輻射安全中心，北京 100082）

氣泡微細化沸騰（MEB）現(xiàn)象具有極高的換熱能力，成功工程化應用后將極大提升核電廠中高熱負荷設備的安全裕量。本文參照以往研究獲得的可視化研究結(jié)果，采用Fluent建立相關模型，綜合考慮氣膜附近Marangoni對流、蒸發(fā)冷凝作用以及溫度對物性參數(shù)的影響，結(jié)合數(shù)值模擬手段和沸騰不穩(wěn)定性分析對MEB現(xiàn)象的發(fā)生機理進行了研究。結(jié)果表明，在不同過冷度下，汽液界面處的蒸汽平均流速隨著壁溫升高而增大。蒸汽平均流速達到該過冷度下Helmholtz失穩(wěn)極限速度時對應的壁溫與在實驗獲得的MEB觸發(fā)壁溫十分接近，說明Helmholtz失穩(wěn)可能是導致MEB現(xiàn)象中氣膜發(fā)生破裂的原因。

氣泡微細化沸騰；Marangoni對流；數(shù)值模擬；Helmholtz失穩(wěn)

為了提高核電機組的經(jīng)濟性指標，大功率反應堆成為了現(xiàn)階段三代堆型的主要發(fā)展方向，這使得嚴重事故后堆芯熔融物增多，提高了熔融物施加在壓力容器內(nèi)壁的熱流密度，極大的降低了熔融物堆內(nèi)滯留技術（In-vessel Retention，IVR）的安全裕量。近幾十年發(fā)現(xiàn)的氣泡微細化沸騰現(xiàn)象（Microbubble Emission Boiling，MEB），由于具有高于一般池沸騰臨界熱流密度（critical heat flux，CHF）的換熱能力，被認為是提高堆內(nèi)熔融物滯留安全裕量的理想手段。

自MEB現(xiàn)象被發(fā)現(xiàn)后，很多學者對氣泡微細化沸騰的換熱特性和沸騰現(xiàn)象進行了研究。Shoji等[1]以鉑絲作為加熱源成功觀察到MEB現(xiàn)象。Tang等以及Suzuki等[2-4]分別以直徑10 mm的銅作為加熱源完成了一系列池式MEB實驗。在他們的研究中發(fā)現(xiàn)，沸騰達到臨界后，銅加熱面上會形成一層氣膜導致壁溫飛升，當壁溫達到一定溫度時氣膜迅速破裂從而使沸騰狀態(tài)進度到MEB階段。根據(jù)最近的實驗研究[2]結(jié)果，隨著液體過冷度的升高，觸發(fā)MEB發(fā)生的壁溫會逐漸降低。由于MEB現(xiàn)象十分劇烈，無法在實驗中直接獲取氣膜或周圍流場的信息，因此本文采用Fluent軟件建立相應幾何模型，結(jié)合數(shù)值模擬方法和沸騰不穩(wěn)定性分析，對MEB現(xiàn)象的觸發(fā)溫度與Helmholtz失穩(wěn)的關系進行了研究，探討了MEB現(xiàn)象的發(fā)生機理。

1　幾何模型

文獻[4]中MEB實驗臺架的加熱面為直徑10 mm的圓形銅面。實驗過程中，通過傾斜角度為18.5°的攝影儀拍攝沸騰現(xiàn)象，由此獲得的剛進入MEB階段時的沸騰現(xiàn)象如圖1所示，加熱面上的不穩(wěn)定氣膜在約1 ms的時間內(nèi)處于與膜態(tài)沸騰狀態(tài)類似的準穩(wěn)態(tài)，隨后迅速發(fā)生破裂。因此可以建立一個氣膜模型，通過數(shù)值模擬手段分析準穩(wěn)態(tài)階段氣膜內(nèi)的流場，從而分析MEB的產(chǎn)生機理。

圖1　氣泡微細化沸騰現(xiàn)象

圖2所示為根據(jù)沸騰現(xiàn)象建立的計算模型，在準穩(wěn)態(tài)階段的氣膜可以簡化為一個扁平的橢球體，橢球體長軸為11.2 mm，短軸為3.52 mm，底部與圓形加熱面相接。與加熱面結(jié)合后，橢球體剩余高度為2.55 mm?？紤]到加熱面和氣膜的對稱性，在數(shù)值模擬過程中只對幾何模型的六分之一進行建模以減少計算量。

圖2　計算域和網(wǎng)格劃分

在預計算階段，分別對比了水溫343.15 K壁溫438.15 K時，網(wǎng)格數(shù)為206 668、287 036、391 101的計算模型得到的蒸汽計算平均速度和最大速度，其結(jié)果如表1所示。其中，網(wǎng)格數(shù)為287 036、391 101模型得到的計算結(jié)果變化很小，因此最終采用的網(wǎng)格數(shù)為287 036。整體網(wǎng)格質(zhì)量的Equi-size skew不超過0.51。

表1　網(wǎng)格敏感性分析

2　數(shù)學物理模型

2.1　計算模型和物性設置

參考以往蒸汽數(shù)值模擬研究[5-6]，湍流模型采用標準k-epsilon模型。由于最終采用網(wǎng)格寬度在0.04 mm左右，導致邊界層處Y Plus值小于15，因此壁面函數(shù)采用Enhanced Wall Treatment。本文采用水作為工質(zhì)進行計算，在計算過程中通過將工質(zhì)的密度設置為不可壓縮理想氣體來模擬密度差導致的自然對流[7]。參考以往仿真經(jīng)驗[8]，模擬過程中考慮了溫度對流體物性的影響，蒸汽的比熱容、導熱系數(shù)和粘度均設置為Fluent材料庫自帶的物性溫度函數(shù)。

2.2　邊界條件設置

計算域內(nèi)，加熱面和汽液界面均設置為壁面邊界條件。加熱面的為恒溫壁面，汽液界面采用對流邊界作為熱力學邊界條件，其中，傳質(zhì)換熱系數(shù)采用公式（1）計算[9]：

fg——汽化潛熱；

v——蒸汽密度；

v——蒸汽溫度；

v——蒸汽壓力；

參考以往MEB現(xiàn)象仿真經(jīng)驗，水作為工質(zhì)時的換熱系數(shù)可設置為230 000 W/m2·K[9]。

由于MEB現(xiàn)象僅在液體過冷度和壁面過熱度足夠高的工況下發(fā)生，目前普遍認為劇烈和蒸發(fā)冷凝作用，以及汽液界面附近存在由表面張力梯度引起的Marangoni對流是導致MEB發(fā)生時氣膜破碎的原因[9]。圖3為汽液界面附近的Marangoni對流的示意圖，在溫度差導致的表面張力梯度作用下，氣相和液相分別受剪切力v和l作用，此時受力平衡可表示為：

式中：d/d——表面張力隨溫度的導數(shù)；

d/d——汽液界面上的溫度梯度。

圖3　汽液界面附近Marangoni對流示意圖

根據(jù)池式MEB現(xiàn)象的可視化研究結(jié)果，氣相的活躍程度遠大于液相，可以認為準穩(wěn)態(tài)階段，表面張力梯度產(chǎn)生的Marangoni Stress完全作用于氣相之上，則公式（2）可簡化為：

據(jù)此，汽液界面處采用Marangoni Stress作為動量邊界條件。水的Marangoni Stress邊界條件應設置為-0.000 18 N/（m·K）[9]。數(shù)值模擬過程中，壓力速度耦合采用SIMPLC算法，其中的畸變修正和壓力方程亞松弛因子均設置為1，其余松弛因子保持默認值[7]。離散格式選項中，梯度采用Least Squares Cell Based，動量、湍流等其余項均采用二階迎風格式。

3　模擬結(jié)果分析

圖4展示了Tang等[2]以水作工質(zhì)的實驗中獲得的MEB沸騰特性曲線。其中，CHF對應的壁面過熱度和MEB觸發(fā)溫度分別隨著過冷度升高呈現(xiàn)升高和降低的趨勢，導致沸騰臨界后的壁溫飛升的幅度逐漸減小。由于MEB觸發(fā)溫度主要集中在壁面過熱度45 K至60 K的區(qū)間內(nèi)，在數(shù)值模擬過程中，以上述區(qū)間作為加熱面溫度的設置范圍，分別計算不同過冷度和壁面過熱度下的氣膜速度場，從而分析MEB現(xiàn)象的產(chǎn)生機理。

圖4　MEB沸騰特性曲線

圖5所示為60 K過冷度下，壁溫為（a）428.15 K和（b）438.15 K時，汽液界面和對稱面上的速度場。在自然對流和Marangoni Stress的作用下，氣膜內(nèi)的蒸汽形成了一個循環(huán)流場，隨著壁溫提升汽液界面上的蒸汽流速也逐漸升高，但流場的形式?jīng)]有發(fā)生變化。通過對比局部流場可以發(fā)現(xiàn)，汽液界面上的蒸汽流速明顯高于其他區(qū)域，因此本文以汽液界面上的蒸汽平均速度作為循環(huán)流場的特征速度。

根據(jù)Helmholtz失穩(wěn)條件[10]，沸騰過程中，導致蒸汽擁塞的氣相極限速度與蒸汽脫離通道波長w的關系可采用式（4）計算。

一般情況下，w等于脫離加熱面氣泡的直徑，在池式MEB沸騰現(xiàn)象中Marangoni對流速度場是軸對稱的，因此采用等效半徑作為氣膜的特征尺寸。根據(jù)以往研究結(jié)果，較小加熱面上的不同階段的氣泡可以近似按軸對稱的回轉(zhuǎn)體考慮。

對于幾何形態(tài)類似圖6所示的準穩(wěn)態(tài)下的氣膜，其等效半徑可以按公式（5）來計算[11]：

圖6　氣膜幾何特征參數(shù)

圖7示出了采用水作為計算工質(zhì)時，數(shù)值模擬結(jié)果與圖4中MEB觸發(fā)溫度的對比，其中共示出了：不同過冷度下采用公式（4）和公式（5）計算的Helmholtz失穩(wěn)極限速度；不同過冷度和壁面過熱度下通過數(shù)值模擬得到的汽液界面上的蒸汽平均速度；在不同過冷度下Helmholtz極限速度曲線上標注了實驗獲得的MEB觸發(fā)溫度。如圖所示，汽液界面附近的蒸汽計算平均速度和Helmholtz失穩(wěn)極限速度均隨著液相過冷度或壁面過熱度的增加而增加。以水作為計算工質(zhì)時，在不同的液相過冷度下，氣相計算平均速度超過Helmholtz極限速度時對應的壁面過熱度與實驗過的MEB觸發(fā)溫度基本一致。說明Helmholtz失穩(wěn)可能是導致氣膜無法繼續(xù)穩(wěn)定存在的原因。

圖7　模擬結(jié)果與實驗值對比

4　結(jié)論

本文根據(jù)MEB可視化研究結(jié)果建立了準穩(wěn)態(tài)下的計算模型，通過引入氣泡等效半徑對MEB現(xiàn)象觸發(fā)溫度進行了研究，得到如下結(jié)論：

（1）在準穩(wěn)態(tài)階段，在自然對流和表面張力梯度驅(qū)動下氣膜內(nèi)形成了一個循環(huán)流場。汽液界面處的蒸汽平均速度隨著壁面過熱度或液相過冷度的升高而升高。

（2）在不同過冷度下，汽液界面處蒸汽平均速度超過Helmholtz失穩(wěn)極限速度時對應的壁溫與實驗中獲得的MEB觸發(fā)溫度一致，說明Helmholtz失穩(wěn)可能是導致氣膜無法繼續(xù)穩(wěn)定存在的原因。

［1］ Shoji M，Yoshihara M．Burnout heat flux of water on a thin wire［C］．Proceeding of 28th National Heat Transfer Symposium of Japan，Japan，1991：121-123．

［2］ Tang J，Mo Z，Sun L，et al．An experimental study on Microbubble Emission Boiling in a subcooled pool：Heat transfer characteristics and visualized presentation［J］．Experimental Thermal and Fluid Science，2017，80：40-52．

［3］ Suzuki K，Inagaki F，Hong C．Subcooled boiling in the ultrasonic field：on the cause of microbubble emission boiling［J］．Heat Transfer Engineering，2011，32（7-8）：673-682．

［4］ Zhu G，Sun L，Tang J，et al．A visualized study of micro-bubble emission boiling［J］．International Communications in Heat and Mass Transfer，2014，59：148-157．

［5］ 邵天，杜亞威，劉燕，等．蒸汽噴射器的三維數(shù)值模擬研究［J］．真空科學與技術學報，2014，34（3）：305-311．

［6］ 史立群，楊建道，楊銳，等．耦合末級葉片的汽輪機排汽缸氣動數(shù)值模擬［J］．動力工程學報，2011，31（9）：655-658．

［7］ 李鵬飛，徐敏義，王飛飛．精通CFD工程仿真與案例實戰(zhàn)［M］．北京：人民郵電出版社，2011．

［8］ 彭偉頔，鄭樂樂，盧川，等．泡核沸騰兩相CFD模擬的參數(shù)敏感性分析與模型驗證［J］．核科學與工程，2018，38（02）：194-203．

［9］ 唐繼國．氣泡微細化沸騰現(xiàn)象極其形成機制研究［D］．哈爾濱工程大學，2016．

［10］魯鐘琪．兩相流與沸騰傳熱［M］．北京：清華大學出版社，2002．

［11］唐繼國，閻昌琪，孫立成，等．氣泡微細化沸騰過程中氣泡冷凝破裂現(xiàn)象［J］．化工學報，2014，65（8）：2902-2907．

Simulation Study on Triggering Temperature of Microbubble Emission Boiling

ZHU Guangyu1,2，GUO Yong1，YUAN Yidan1，LIU Yusheng2，LI Wei1

（1. Center for Nuclear Power Safety，China Nuclear Power Engineering Co.，Ltd.，Beijing 100840，China；2. Nuclear and Radiation Safety Center，MEE，Beijing 100082，China）

Due to its extremely high heat transfer capability，Microbubble Emission Boiling（MEB）is considered as the ideal means to promote heat exchanger in nuclear power plant．In this paper，the factor triggering MEB was studied by numerical work and boiling instability analysis．A numerical mode with the consideration of Marangoni stress，heat transfer on the vapor-liquid interface and the thermophysical property of test fluids was built based on the visualized study data．The results showed that，the average velocity near the vapor-liquid interface increases while the wall temperature increases．Furthermore，The wall temperature which the average vapor velocity equals to Helmholtz instability limiting speed is consistent with the temperature triggering MEB in the experiment．Thus，Helmholtz instability may be the reason for the collapsing of vapor film when MEB occurs．

Microbubble emission boiling；Marangoni convection；Numerical simulation；Helmholtz instability

TL334

A

0258-0918（2021）03-0639-05

2020-11-27

國家重點研發(fā)計劃資助（2018YFB1900100）

朱光昱（1989—），男，黑龍江人，工程師，碩士研究生，現(xiàn)主要從事核電廠嚴重事故及緩解措施方面研究