氣泡微細化沸騰觸發溫度數值模擬研究

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

氣泡微細化沸騰觸發溫度數值模擬研究

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

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

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

氣泡微細化沸騰；Marangoni對流；數值模擬；Helmholtz失穩

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

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

1　幾何模型

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

圖1　氣泡微細化沸騰現象

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

圖2　計算域和網格劃分

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

表1　網格敏感性分析

2　數學物理模型

2.1　計算模型和物性設置

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

2.2　邊界條件設置

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

fg——汽化潛熱；

v——蒸汽密度；

v——蒸汽溫度；

v——蒸汽壓力；

參考以往MEB現象仿真經驗，水作為工質時的換熱系數可設置為230 000 W/m2·K[9]。

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

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

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

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

根據池式MEB現象的可視化研究結果，氣相的活躍程度遠大于液相，可以認為準穩態階段，表面張力梯度產生的Marangoni Stress完全作用于氣相之上，則公式（2）可簡化為：

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

3　模擬結果分析

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

圖4　MEB沸騰特性曲線

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

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

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

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

圖6　氣膜幾何特征參數

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

圖7　模擬結果與實驗值對比

4　結論

本文根據MEB可視化研究結果建立了準穩態下的計算模型，通過引入氣泡等效半徑對MEB現象觸發溫度進行了研究，得到如下結論：

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

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

［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］ 邵天，杜亞威，劉燕，等．蒸汽噴射器的三維數值模擬研究［J］．真空科學與技術學報，2014，34（3）：305-311．

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

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

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

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

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

［11］唐繼國，閻昌琪，孫立成，等．氣泡微細化沸騰過程中氣泡冷凝破裂現象［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

國家重點研發計劃資助（2018YFB1900100）

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