單汽泡沸騰過程數(shù)值模擬的研究

劉俊杰，王國清，張利軍，姚詩訓

（1. 中國石化 北京化工研究院，北京 100013；2. 卡內(nèi)基梅隆大學 機械系，美國 賓夕法尼亞州 匹茲堡，15213-3890）

單汽泡沸騰過程數(shù)值模擬的研究

劉俊杰1，2，王國清1，張利軍1，姚詩訓2

（1. 中國石化 北京化工研究院，北京 100013；2. 卡內(nèi)基梅隆大學 機械系，美國 賓夕法尼亞州 匹茲堡，15213-3890）

利用Comsol Multiphysics軟件中的Level Set方法對單汽泡沸騰過程進行了模擬，建立了微細結構網(wǎng)格，給出了邊界條件，求解了質(zhì)量、動量、能量和Level Set方程，得到了正確的模擬結果。分析了單汽泡沸騰過程中相含率、壓力場、速度場和溫度場隨時間的變化規(guī)律；考察了接觸角和壁面效應對汽泡脫離直徑和汽泡生長周期的影響。模擬結果表明，當接觸角小于12°時，汽泡脫離直徑維持不變；當接觸角大于12°時，接觸角越大，汽泡脫離直徑越大。汽泡生長周期隨接觸角的增大而延長。當加熱裝置的直徑小于汽泡脫離直徑時，汽泡受壁面效應的影響不易脫離；當加熱裝置的直徑大于汽泡脫離直徑時，汽泡脫離直徑不隨加熱裝置直徑的增大而變化；汽泡生長周期隨加熱裝置直徑的增大而縮短。

單汽泡；沸騰傳熱；數(shù)值模擬；接觸角；壁面效應；Comsol Multiphysics軟件；Level Set方法

沸騰傳熱在石油化工領域應用十分廣泛，單汽泡沸騰傳熱研究是其研究的基礎［1-7］。美國加州大學洛杉磯分校的Dhir教授和德國達姆施塔特工業(yè)大學的Stephan教授在單汽泡實驗研究和理論研究方面都取得了重要成果［8］。在實驗研究方面，他們設計了單汽泡實驗平臺，利用高速攝像機、熱電偶和紅外攝像機對單汽泡的動力學和熱力學特性進行了研究［9-12］。Dhir教授研究了微重力條件下單汽泡的脫離直徑和生長周期［13-14］。Stephan教授考察了汽液固三相接觸區(qū)的溫度場和熱通量隨時間的變化規(guī)律。在理論研究方面，他們對單汽泡的沸騰傳熱進行了過程模擬，耦合了汽液固三相接觸區(qū)與宏觀汽泡生長過程中的質(zhì)量、動量和熱量，深入揭示了單汽泡沸騰傳熱的機理。單汽泡沸騰傳熱的汽液相界面模擬主要有兩種數(shù)學方法：基于動網(wǎng)格的Volume of Fluid（簡稱VOF）方法和基于固定網(wǎng)格的Level Set方法。Dhir教授利用Fortran語言編寫了求解 Level-Set方程的程序，對單汽泡沸騰過程進行了模擬，但前處理軟件較簡單，模擬精度不高［15-16］。Stephan教授利用Open Foam 軟件中的VOF模型對單汽泡沸騰過程進行了模擬，根據(jù)實驗結果修改了汽液界面中的部分參數(shù)，提高了計算精度［17-18］。國內(nèi)微重力實驗室的科研人員也進行了相關研究，利用Level Set 方法研究了單汽泡沸騰過程中的汽泡動力學，但僅給出了初步模擬結果［19-23］。在上述模擬過程中，研究者均采用自己開發(fā)的程序，而沒有使用商用軟件，通用性較差，給后續(xù)研究帶來了障礙。佐治亞理工大學的Smith教授利用Cosmol Multiphysics軟件中的Arbitrary Lagrangian-Eulerian和Level Set方法對高溫噴射氣體在飽和沸水中的氣泡生長過程進行了模擬，Cosmol Multiphysics軟件自帶的例子也是采用相場模型對上述過程進行了模擬［24-27］。高溫氣體在沸水中的長大過程與沸水中汽泡長大過程的物理機制不同，盡管具有一定的借鑒意義，但兩者在模擬過程中存在較大差異。

本工作嘗試利用基于有限元技術［28-30］的Comsol Multiphysics軟件中的Level Set方法對單汽泡沸騰過程進行初步模擬。

1 單汽泡沸騰過程數(shù)值模型的建立

1.1 單汽泡計算區(qū)域?qū)嶓w建模

二維軸對稱實體建模的直徑為5 mm，高5 mm；人工汽穴的直徑為10 μm，高15 μm。網(wǎng)格尺寸從4 μm到20 μm，總自由度是40萬。

1.2 數(shù)值計算模型

1.2.1 質(zhì)量、動量和能量方程的數(shù)學模型

能量方程：

質(zhì)量汽化速率：

汽液界面兩側切向應力方程：

汽液界面兩側壓力方程：

物性參數(shù)方程：

1.2.2 求解方法

質(zhì)量、動量、能量以及Level Set方程采用有限元方法離散。底部區(qū)域為恒溫熱源，邊界溫度為383 K。與壁面接觸的熱邊界層厚度為0.625 mm，溫度分布為線性，此區(qū)域以上直至4 mm處為飽和液體，在4 mm之上為飽和汽體，出口邊界條件為0.1 MPa。模擬所采用的計算機為惠普Z420，含有8核Intel xeon（R） CPU，48 G內(nèi)存。一個算例所需的計算時間為24 h。

2 結果與討論

2.1 單汽泡的相含率隨時間的變化規(guī)律

單汽泡的三維模擬形態(tài)見圖1。

圖1 單汽泡的三維模擬形態(tài)Fig.1 Three dimension simulation of a single bubble.

由于接觸角較小，僅為25°，汽泡形狀為球形。為了更好地顯示汽泡隨時間的變化規(guī)律，圖2給出了汽泡逐漸長大、脫離和再生成的過程。汽泡由最初的0.1 mm逐漸長到1.2 mm，汽泡根部直徑則由小到大，再逐漸由大到小，直至變?yōu)?。在前10 ms，汽泡生長速度較快，直徑從0.01 mm迅速長到0.7 mm，隨后汽泡生長速度減慢，在37 ms時汽泡脫離了加熱底部，直徑達到1.2 mm。此結果與文獻［14］報道的實驗結果和模擬結果一致。汽泡逐漸長大是由于汽泡周圍過熱液體變成蒸汽，汽液之間存在傳質(zhì)。在汽泡成長的初期階段，汽泡的直徑較小且完全處于熱邊界層中，過飽和液體的溫度相對較高，液體變成汽體的質(zhì)量較多，因而汽泡生長速度較快；當汽泡直徑長到一定程度時，部分汽泡位于飽和液體中，汽液兩側溫度相等，汽液之間無質(zhì)量傳遞發(fā)生，僅靠汽泡下部過熱液體提供汽泡成長動力，汽泡生長速度降低。

圖2 單汽泡的相含率隨時間變化的規(guī)律Fig.2 Single bubble holdup vs.time.

2.2 單汽泡的壓力場隨時間變化的規(guī)律

單汽泡的壓力場隨時間變化的規(guī)律見圖3。由圖3可見，汽液相變導致汽泡內(nèi)的壓力大于外部液相壓力，汽泡逐漸增大。液相壓力沿高度逐漸降低，且隨時間變化較小。在汽泡成長的初始階段，汽泡內(nèi)壓力增加較快，汽相壓力與外部液相壓力相差較大，隨加熱時間的延長，汽泡內(nèi)外的壓差逐漸縮小，當汽泡脫離加熱壁面時，汽泡內(nèi)外壓力趨于一致。

2.3 單汽泡的速度場隨時間變化的規(guī)律

單汽泡的速度場隨時間變化的規(guī)律見圖4，等勢線為汽液界面。單汽泡沸騰傳熱過程中汽液兩相的雷諾數(shù)均相對較小，均處于層流狀態(tài)。從圖4可看出，汽泡內(nèi)的汽相速度遠大于汽泡外的液相速度。汽泡內(nèi)部流場顯示：一部分液體從汽泡根部汽化，在汽泡長大過程中汽泡內(nèi)部的汽體逐漸形成環(huán)流。隨著汽泡脫離壁面，汽泡周圍的液體迅速填補汽泡離開時的空白，在汽泡周圍逐漸形成尾渦。

圖3 單汽泡的壓力場隨時間變化的規(guī)律Fig.3 Single bubble pressure field vs.time.

圖4 單汽泡的速度場隨時間變化的規(guī)律Fig.4 Single bubble velocity field vs.time.

2.4 單汽泡的溫度場隨時間變化的規(guī)律

單汽泡的溫度場隨時間變化的規(guī)律見圖5。由圖5可見，底部液體區(qū)域存在著一個熱邊界層，在此熱邊界層區(qū)域，液體處于過熱狀態(tài)。由于底部加熱，液相溫度場從底部向上呈線性衰減。熱邊界層區(qū)域為汽泡長大提供熱源。在加熱過程中熱邊界層逐漸發(fā)展，但由于加熱時間較短，熱邊界層變化較小，液體溫度基本保持不變。汽泡與熱邊界層接觸區(qū)域隨汽泡形狀的變化而變化。

2.5 接觸角對汽泡脫離直徑和生長周期的影響

接觸角對汽泡脫離直徑和生長周期的影響見圖6和圖 7。所謂接觸角是指在汽液固三相交點處所作的汽液界面的切線穿過液體與固液交界線之間的夾角。當接觸角小于12°時，汽泡脫離直徑均為1 mm左右，接觸角對汽泡脫離直徑的影響可忽略不計。當接觸角大于12°時，汽泡脫離直徑隨接觸角的增大而顯著增大。汽泡生長周期隨接觸角的增大而延長，從16 ms逐漸延長至36 ms。這是由于隨接觸角的增大，汽泡與壁面的作用力逐漸增強，即液體表面張力逐漸增強，汽泡不易脫離壁面。

圖5 單汽泡的溫度場隨時間變化的規(guī)律Fig.5 Single bubble temperature field vs.time.

圖6 接觸角對汽泡脫離直徑的影響Fig.6 Effect of contact angle on the bubble departure diameter.

圖7 接觸角對汽泡生長周期的影響Fig.7 Effect of contact angle on the bubble growing period.

2.6 壁面效應對汽泡脫離直徑和生長周期的影響

加熱裝置的直徑對汽泡脫離直徑和生長周期的影響見圖8和圖9。

當加熱裝置的直徑小于在此溫度下的汽泡脫離直徑時，由于壁面的阻擋效應，汽泡不易長大，數(shù)值模擬過程發(fā)生周期性的震蕩。當加熱裝置的直徑大于此溫度下的汽泡脫離直徑時，隨加熱裝置直徑的增大，汽泡脫離直徑未發(fā)生任何變化。這是由于無壁面阻擋作用且初始溫度場保持不變時，液體熱邊界層保持穩(wěn)定，汽泡外部液體的汽化動力未變。汽泡生長周期則隨加熱裝置直徑的減小而逐漸延長，從最初的38 ms延長到61 ms。這是由于汽泡上升過程中受到壁面效應的影響，汽泡所受曵力增大，而其升力保持不變。

圖8 不同加熱裝置直徑下單汽泡脫離壁面時的相含率Fig.8 Bubble holdup in the bubble separation with the wall at different heating foil diameter.

圖9 加熱裝置的直徑對汽泡生長周期的影響Fig.9 Effect of the heating foil diameter on the bubble growing period.

在單汽泡沸騰模擬過程中，本工作未考慮汽液固三相接觸線區(qū)域的液膜蒸發(fā)過程，將在后續(xù)工作中展開。

3 結論

1）利用Comsol Multiphysics軟件對單汽泡沸騰過程進行了模擬，考察了單汽泡逐漸長大的過程，給出了單汽泡的相含率、壓力場、速度場和溫度場隨時間變化的規(guī)律，并進行了合理的解釋。

2）當接觸角小于12°時，汽泡脫離直徑維持不變；當接觸角大于12°時，隨接觸角的增大，汽泡脫離直徑逐漸增大。汽泡生長周期隨接觸角的增大而延長。

3）當加熱裝置的直徑小于汽泡脫離直徑時，汽泡受壁面效應的影響不易脫離；當加熱裝置的直徑大于汽泡脫離直徑時，汽泡脫離直徑保持不變。汽泡生長周期隨加熱裝置直徑的增大而縮短。

符號說明

C 汽液相變公式中的常量

Cp流體的比熱容，J/（kg?K）

Fst體積力

hlv汽液相變焓值，J/kg

I 單位向量

k 流體的熱導率，W/（m?K）

n 汽液界面的單位法向量

p 流體壓力，Pa

T 液體的溫度，K

Tsat飽和液體的溫度，K

t 單位時間，s

u 流體的速度，m/s

Vl液相的體積，L

γ 流場中速度的最大值，m/s

δ 汽液界面的厚度，1/m

ε 模型中網(wǎng)格尺寸的最大值，m

μ 流體黏度，Pa?s

ρ 流體密度，kg/m3

σ 表面張力，N/m

φ Level Set的相界面函數(shù)

下角標

l 液相

v 汽相

［1］Goldste in R J，Eckert E R G，Ibele W E，et al.Heat Transfer：A Review of 1999 Literature［J］.Int J Heat Mass Transfer，2001，44（19）：3579-3699.

［2］Plesset M S，Zwick S A.The Growth of Vapor Bubbles in Superheated Liquids［J］.J Appl Phys，1954，25（4）：493-500.

［3］Jiang Yuyan，Osada Hiroshi，Inagaki Masahide，et al.Dynamic Modeling on Bubble Growth，Detachment and Heat Transfer for Hybrid-Scheme Computations of Nucleate Boiling［J］.Int J Heat Mass Transfer，2013，56（1/2）：640-652.

［4］Yi H，Tipton J，Kihm K D，et al.Effect of Disjoining Pressure（Ⅱ） on Multi-Scale Modeling for Evaporative Liquid Metal（Na） Capillary［J］.Int J Heat Mass Transfer，2014，78（6）：137-149.

［5］Son G.A Numerical Method for Bubble Motion with Phase Change［J］.Numer Heat Transfer，Part B，2001，39（5）：509-523.

［6］Dhavaleswarapu H K，Garimella S V，Murthy J Y.Microscale Temperature Measurements Near the Triple Line of an Evaporating Thin Liquid Film［J］.J Heat Transfer，2009，131（6）：061501-1-061501-7.

［7］Nam Youngsuk，Wu Jinfeng，Warrier G，et al.Experimental and Numerical Study of Single Bubble Dynamics on a Hydrophobic Surface［J］.J Heat Transfer，2009，131（12）：121004-1-121004-7.

［8］Dhir V K，Warrier G R，Aktinol E.Numerical Simulation of Pool Boiling：A Review［J］.J Heat Transfer，2013，135（6）：061502-1-061502-17.

［9］Wagner E，Stephan P.High-Resolution Measurements at Nucleate Boiling of Pure FC-84 and FC-3284 and Its Binary Mixtures［J］.J Heat Transfer，2009，131（12）：121008-1-121008-12.

［10］Migliaccio C P，Dhavaleswarapu H K，Garimella S V.Temperature Measurements Near the Contact Line of an Evaporating Meniscus V-Groove［J］.Int J Heat Mass Transfer，2011，54（7/8）：1520-1526.

［11］Nama Youngsuk，Aktinola E，Dhir V K，et al.Single Bubble Dynamics on a Superhydrophilic Surface with Artificial Nucleation Sites［J］.Int J Heat Mass Transfer，2011，54（7/8）：1572-1577.

［12］Son G，Dhir V K，Ramanujapu N.Dynamics and Heat Transfer Associated with a Single Bubble During Nucleate Boiling on a Horizontal Surface［J］.J Heat Transfer，1999，121（3）：623-631.

［13］Aktinol E，Dhir V K.Numerical Simulation of Nucleate Boiling Phenomenon Coupled with Thermal Response of the Solid［J］.Microgravity Sci Technol，2012，24（4）：255-265.

［14］Mukherjee A.Numerical and Experimental Study of Lateral Merge of Vapor Bubbles Formed on a Horizontal Surface During Nucleate Pool Boiling［D］.Los Angeles：University of California，2004.

［15］Kunkelmanna C.Numerical Modeling and Investigation of Boiling Phenomena［D］.Darmstadt：Technische Universit?t Darmstadt，2011.

［16］Son G，Dhir V K.A Level Set Method for Analysis of Film Boiling on an Immersed Solid Surface［J］.Numer Heat Transfer，Part B，2007，52（2）：153-177.

［17］Kunkelmanna C，Stephan P.Numerical Simulation of the Transient Heat Transfer During Nucleate Boiling of Refrigerant HFE-7100［J］.Int J Refrigerat，2010，33（7）：1221-1228.

［18］Kunkelmanna C，Stephan P.CFD Simulation of Boiling Flows Using the Volume-of-Fluid Method Within Open Foam［J］.Numer Heat Transfer，Part A，2009，56（8）：631-646.

［19］李震東，趙建福，李晶，等.單氣泡池沸騰過程中的氣泡動力學數(shù)值模擬研究［J］.力學與工程應用，2010，13（增刊）：53-59.

［20］李震東，趙建福，秦文韜.氣泡/液滴運動的Level Set方法模擬研究［C］//中國力學學會學術大會.鄭州：中國力學學會，2009：410-418.

［21］陳凡紅，王成，郝莉，等.用Level Set方法追蹤運動界面［J］.力學與實踐，2006，28（4）：22-26.

［22］管鵬.流動沸騰中汽泡行為的理論與實驗研究［D］.北京：北京交通大學，2013.

［23］薛龍昌.基于氣體動力論界面質(zhì)量輸運和微液層蒸發(fā)模型的過冷沸騰汽泡生長特性研究［D］.重慶：重慶大學，2014.

［24］Forster C J，Smith M K.The Transient Modeling of Bubble Pinch-Off Using an ALE Moving Mesh［EB/OL］.［2015-09-15］.https：//br.comsol.com/paper/download/62616/smith_paper.pdf.

［25］Forster C J，Smith M K.The Transient Modeling of Single-Bubble Nucleate Boiling in a Sub Cooled Liquid Using an ALE Moving Mesh［EB/OL］.［2015-09-15］.www.comsol.com/paper/download/83905/forster_paper.pdf.

［26］Zhou Hongwei，Gué A M.Simulation Bubble Nucleation and Bubble Growth of a Thermal Bubble Microejector［EB/ OL］.［2015-09-15］.https：//ch.comsol.com/paper/download/45165/Hongwei.pdf.

［27］Cinzia Iacovelli.Boiling Water，a Heat Transfer and CFD Problem［EB/OL］.［2015-09-15］.http：//www.comsol.com/ blogs/boiling-water-heat-transfer-and-cfd-problem/.

［28］曾攀.有限元基礎教程［M］.北京：高等教育出版社，2009：7-20.

［29］畢超.計算流體力學有限元方法及其編程詳解［M］.北京：機械工業(yè)出版社，2013：1-10.

［30］卡坦.MATLAB有限元分析與應用［M］.北京：清華大學出版社，2004：1-9.

（編輯 王 萍）

Preliminary Study on Numerical Simulation of Single Bubble Boiling

Liu Junjie1，2，Wang Guoqing1，Zhang Lijun1，Yao Shi-chune2

（1.SINOPEC Beijing Research Institute of Chemical Industry，Beijing 100013，China；2.Department of Mechanical Engineering，Carnegie Mellon University，Pittsburgh PA 15213-3890，USA）

The single bubble boiling process was simulated by means of the Level Set equation in the Comsol Multiphysics software，the microstructure network was established，boundary conditions were provide，and the mass，momentum，energy and Level Set equation were solved.The satisfied simulation results were acquired.The changes of the holdup，pressure field，velocity field and temperature field with time were studied.The influences of contact angle and wall effect on the bubble departure diameter and bubble growing period were investigated.It was showed that，when the contact angle was bigger than 12°，the bigger the contact angle，the bigger the departure diameter and the longer the bubble growing period.But the bubble departure diameter was constant when the contact angle was less than 12°.As the heating foil diameter was less than the bubble departure diameter，the bubble was not easy to grow up and leave the wall.When the heating foil diameter was more than the bubble departure diameter，the bubble departure diameter was constant and the bubble growing period decreased with the heating foil diameter increasing.

single bubble；boiling heat transfer；numerical simulation；contact angle；wall effect；Comsol Multiphysics software；Level Set equation

1000-8144（2015）11-1295-07

TQ 018

A

2015-09-22；［修改稿日期］2015-10-08。

劉俊杰（1976—），男，山東省寧津縣人，博士，高級工程師，電話 010-59202725，電郵 liujj.bjhy@sinopec.com。