劉俊杰，王國(guó)清，張利軍，姚詩(shī)訓(xùn)

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

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

劉俊杰1，2，王國(guó)清1，張利軍1，姚詩(shī)訓(xùn)2

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

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

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

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

本工作嘗試?yán)没谟邢拊夹g(shù)［28-30］的Comsol Multiphysics軟件中的Level Set方法對(duì)單汽泡沸騰過(guò)程進(jìn)行初步模擬。

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

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

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

1.2 數(shù)值計(jì)算模型

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

能量方程：

質(zhì)量汽化速率：

汽液界面兩側(cè)切向應(yīng)力方程：

汽液界面兩側(cè)壓力方程：

物性參數(shù)方程：

1.2.2 求解方法

質(zhì)量、動(dòng)量、能量以及Level Set方程采用有限元方法離散。底部區(qū)域?yàn)楹銣責(zé)嵩矗吔鐪囟葹?83 K。與壁面接觸的熱邊界層厚度為0.625 mm，溫度分布為線性，此區(qū)域以上直至4 mm處為飽和液體，在4 mm之上為飽和汽體，出口邊界條件為0.1 MPa。模擬所采用的計(jì)算機(jī)為惠普Z(yǔ)420，含有8核Intel xeon（R） CPU，48 G內(nèi)存。一個(gè)算例所需的計(jì)算時(shí)間為24 h。

2 結(jié)果與討論

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

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

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

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

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

2.2 單汽泡的壓力場(chǎng)隨時(shí)間變化的規(guī)律

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

2.3 單汽泡的速度場(chǎng)隨時(shí)間變化的規(guī)律

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

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

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

2.4 單汽泡的溫度場(chǎng)隨時(shí)間變化的規(guī)律

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

2.5 接觸角對(duì)汽泡脫離直徑和生長(zhǎng)周期的影響

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

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

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

圖7 接觸角對(duì)汽泡生長(zhǎng)周期的影響Fig.7 Effect of contact angle on the bubble growing period.

2.6 壁面效應(yīng)對(duì)汽泡脫離直徑和生長(zhǎng)周期的影響

加熱裝置的直徑對(duì)汽泡脫離直徑和生長(zhǎng)周期的影響見(jiàn)圖8和圖9。

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

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

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

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

3 結(jié)論

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

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

3）當(dāng)加熱裝置的直徑小于汽泡脫離直徑時(shí)，汽泡受壁面效應(yīng)的影響不易脫離；當(dāng)加熱裝置的直徑大于汽泡脫離直徑時(shí)，汽泡脫離直徑保持不變。汽泡生長(zhǎng)周期隨加熱裝置直徑的增大而縮短。

符號(hào)說(shuō)明

C 汽液相變公式中的常量

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

Fst體積力

hlv汽液相變焓值，J/kg

I 單位向量

k 流體的熱導(dǎo)率，W/（m?K）

n 汽液界面的單位法向量

p 流體壓力，Pa

T 液體的溫度，K

Tsat飽和液體的溫度，K

t 單位時(shí)間，s

u 流體的速度，m/s

Vl液相的體積，L

γ 流場(chǎng)中速度的最大值，m/s

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

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

μ 流體黏度，Pa?s

ρ 流體密度，kg/m3

σ 表面張力，N/m

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

下角標(biāo)

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］李震東，趙建福，李晶，等.單氣泡池沸騰過(guò)程中的氣泡動(dòng)力學(xué)數(shù)值模擬研究［J］.力學(xué)與工程應(yīng)用，2010，13（增刊）：53-59.

［20］李震東，趙建福，秦文韜.氣泡/液滴運(yùn)動(dòng)的Level Set方法模擬研究［C］//中國(guó)力學(xué)學(xué)會(huì)學(xué)術(shù)大會(huì).鄭州：中國(guó)力學(xué)學(xué)會(huì)，2009：410-418.

［21］陳凡紅，王成，郝莉，等.用Level Set方法追蹤運(yùn)動(dòng)界面［J］.力學(xué)與實(shí)踐，2006，28（4）：22-26.

［22］管鵬.流動(dòng)沸騰中汽泡行為的理論與實(shí)驗(yàn)研究［D］.北京：北京交通大學(xué)，2013.

［23］薛龍昌.基于氣體動(dòng)力論界面質(zhì)量輸運(yùn)和微液層蒸發(fā)模型的過(guò)冷沸騰汽泡生長(zhǎng)特性研究［D］.重慶：重慶大學(xué)，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］曾攀.有限元基礎(chǔ)教程［M］.北京：高等教育出版社，2009：7-20.

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

［30］卡坦.MATLAB有限元分析與應(yīng)用［M］.北京：清華大學(xué)出版社，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—），男，山東省寧津縣人，博士，高級(jí)工程師，電話 010-59202725，電郵 liujj.bjhy@sinopec.com。