摘?要：數(shù)值模擬是探究沸騰換熱機(jī)理的重要方法。本文對使用VOF方法模擬氣泡生長的數(shù)值方法進(jìn)行了調(diào)研，著重總結(jié)了模擬過程中所使用的相變模型，并對相變模型進(jìn)行了簡要的評析，使用其中兩個模型模擬了一維Stefan蒸發(fā)和冷凝問題，驗(yàn)證了模型準(zhǔn)確度。

關(guān)鍵詞：VOF方法；相變模型；數(shù)值模擬

1?概述

沸騰是指液體內(nèi)部生成氣泡或氣相的一種劇烈的汽化過程。沸騰換熱則指該過程中的熱量傳遞。沸騰換熱由于其優(yōu)秀的換熱能力，廣泛地應(yīng)用于制冷、發(fā)電、化工等領(lǐng)域。沸騰過程中伴隨著氣泡成核、生長、聚并等行為，這些氣泡行為體現(xiàn)了兩相間的質(zhì)量、能量、動量傳遞。對這些氣泡行為進(jìn)行研究，將有助于完善沸騰換熱機(jī)理，推廣沸騰換熱的工程應(yīng)用。

對于兩相流模擬，模型主要分為兩大類：高相分?jǐn)?shù)模型和界面類模型。前者適用于氣泡特別多的流動，不關(guān)注氣泡界面和氣泡形狀，著重關(guān)注含氣率；后者適用于需要捕捉相界面的情況。常用的界面捕捉方法包括Level?Set方法和VOF（Volume?of?Fluid）方法，其中VOF方法由于固有的質(zhì)量守恒特性，在許多CFD軟件中得到應(yīng)用。由于氣泡生長過程中存在相變，因此需要選擇合適的傳熱傳質(zhì)模型結(jié)合界面捕捉方法，才能夠較為準(zhǔn)確地模擬氣泡的生長。

從實(shí)際需求出發(fā)，本文對適用于VOF方法的相變模型進(jìn)行調(diào)研總結(jié)，并對不同模型進(jìn)行了分析比較，以便后續(xù)模擬研究中選擇合適的模型。

2?控制方程

在相變研究中，研究人員提出了許多相變模型來對氣泡生長進(jìn)行模擬。這些相變封閉模型每個都適合不同的應(yīng)用，但大多數(shù)都采用通用的應(yīng)用形式，即將相變源項(xiàng)應(yīng)用于質(zhì)量、能量和相分?jǐn)?shù)的控制方程［1］。

黏性不可壓縮流體的動量方程、質(zhì)量方程和能量方程如下：

ρu→t+u→·SymbolQC@

u→=-SymbolQC@

p+SymbolQC@

·μSymbolQC@

u→+SymbolQC@

u→T+F→（1）

SymbolQC@

·u→=1ρv-1ρlm·（2）

ρcpTt+u→·SymbolQC@

T=SymbolQC@

·（λSymbolQC@

T）-m·hlv（3）

上式中，ρ表示密度，u→表示速度，p表示壓力，μ表示動力黏度，F(xiàn)→包括重力和表面張力，m·表示相變速率，cp表示比熱容，T表示溫度，λ表示熱導(dǎo)率，hlv為液體的潛熱。下標(biāo)v和l分別表示氣相和液相。表面張力通過連續(xù)表面力（Continuum?Surface?Force，CSF）模型［2］，作用在相界面處：

Fσ=-σSymbolQC@

·SymbolQC@

αSymbolQC@

αSymbolQC@

α（4）

如前文中所述，需要使用VOF方法對界面進(jìn)行捕捉。VOF方法使用體積分?jǐn)?shù)α來表征網(wǎng)格單元中的相分布。α=0，表示該網(wǎng)格為純氣相；0<α<1，表示該網(wǎng)格為兩相混合，即為相界面所處的網(wǎng)格；α=1，表示該網(wǎng)格為純液相。求解體積分?jǐn)?shù)α的控制方程如下：

αt+SymbolQC@

·αu→=-m·1ρl（5）

3?相變模型

如上所述，針對相變研究需要選擇合適的相變模型來計算相變源項(xiàng)。本節(jié)將對不同學(xué)者使用VOF方法模擬氣泡生長時常用的相變模型進(jìn)行總結(jié)。

3.1?Lee相變模型

根據(jù)Lee［3］提出的相變模型，蒸發(fā)和冷凝過程傳質(zhì)速率分別表示為：

m·e=reαρlT-TsatTsat蒸發(fā)（6）

m·c=rc（1-α）ρvT-TsatTsat冷凝（7）

上式中，re和rc分別表示為蒸發(fā)和冷凝的傳質(zhì)強(qiáng)度因子，Tsat為飽和溫度。re和rc一般按照經(jīng)驗(yàn)選取，過大會導(dǎo)致數(shù)值收斂問題，取得過小則會導(dǎo)致界面溫度與飽和溫度的顯著偏差。

3.2?Sato相變模型

Sato等人［4］提出了一個簡單而直接的相變模型，從溫度場計算出的相變速率直接作為體積守恒的源項(xiàng)。相變速率m·定義為：

m·=M·Ai/V（8）

上式中，M·表示界面相變速率（kg/m2s），V表示網(wǎng)格單元的體積。因此對于不包含相界面的網(wǎng)格，其相變速率m·為0。界面相變速率M·定義為：

M·=（λlSymbolQC@

Tl·n→+-λvSymbolQC@

Tv·n→）/hlv（9）

上式中，n→為界面的法向量，由氣相指向液相。假設(shè)氣液相界面溫度為飽和溫度，用于計算溫度梯度。

3.3?Onishi相變模型

Onishi［5］使用溫度恢復(fù)法，將相變的速率與模擬的時間步長進(jìn)行關(guān)聯(lián)，得到如下的相變模型：

m·=ρcpTcell-Tsathlv·Δt（10）

上式中Tcell為相界面所在網(wǎng)格單元的平均溫度，ρ為相界面所在網(wǎng)格單元的平均密度。Rattner［6］認(rèn)為應(yīng)當(dāng)對式（10）的相變速率進(jìn)行限制。在一個時間步內(nèi)，網(wǎng)格單元內(nèi)蒸發(fā)的質(zhì)量不應(yīng)超過液相質(zhì)量，冷凝的質(zhì)量不應(yīng)超過氣相的質(zhì)量，對應(yīng)的蒸發(fā)以及冷凝傳質(zhì)速率如下：

m·e，lim=αρlΔt蒸發(fā)（11）

m·c，lim=-（1-α）ρvΔt冷凝（12）

類似于庫朗數(shù)的限制，Rattner認(rèn)為［6］，為了保證模擬過程的數(shù)值穩(wěn)定性，一個時間步內(nèi)蒸發(fā)產(chǎn)生的氣體體積或冷凝產(chǎn)生的液體體積不應(yīng)超過該網(wǎng)格的體積，對應(yīng)的蒸發(fā)和冷凝速率如下：

m·e，CFL=1Δt1ρv-1ρl-1蒸發(fā)（13）

m·c，CFL=-1Δt1ρv-1ρl-1冷凝（14）

因而，模擬相變時，相變速率模型的選取有如下限制：

m·e=minm·，m·e，lim，m·e，CFL蒸發(fā)（15）

m·c=maxm·，m·c，lim，m·c，CFL冷凝（16）

對于模擬過程中的時間步長，除了滿足CFL條件外，還應(yīng)滿足熱擴(kuò)散穩(wěn)定條件：

ΔtSymbolcB@

lminΔ2/λρcpeffi（17）

上式中，l為一個自定義的約束，對于二維情況取1/4，對于三維情況取1/6；Δ為網(wǎng)格的最小長度。

4?模型分析驗(yàn)證

上文中所述三個相變模型是目前VOF方法模擬相變時，使用較為廣泛的三個模型。Lee模型由于簡單而被廣泛使用，但它是一種經(jīng)驗(yàn)?zāi)Ｐ?，模型中的傳質(zhì)強(qiáng)度因子對于不同的工況需要選擇不同的值。Onishi模型以及改良后的Rattner模型，采用溫度恢復(fù)法將相變速率與時間步長關(guān)聯(lián)，但需要對傳熱速率以及模擬的時間步長進(jìn)行限制，不然結(jié)果將顯著失真。Sato模型根據(jù)界面兩側(cè)的溫度梯度差確定相變速率，結(jié)果更為準(zhǔn)確；但是該方法需要對相界面進(jìn)行幾何重構(gòu)，得到界面面積、方向、位置等詳細(xì)的信息，在OpenFOAM中使用該方法較為困難。因此出于后續(xù)研究的需求，不考慮使用Sato模型。

目前Lee模型和Rattner模型已經(jīng)被一些研究人員將其植入開源程序OpenFOAM中。本節(jié)將分別使用Lee模型和Rattner模型模擬一維Stefan蒸發(fā)問題和冷凝問題，考察兩模型的適用性與準(zhǔn)確性。

4.1?Stefan蒸發(fā)問題

圖1為一維Stefan蒸發(fā)問題的示意圖。初始階段整個計算域?yàn)轱柡鸵后w，左壁保持高于飽和溫度的恒定溫度，受左壁影響，液體開始蒸發(fā)，蒸汽層厚度逐漸增加。對于一維Stefan蒸發(fā)問題，其解析解為：

δt=2η?λvtρvcp，v（18）

上式中η由下式獲得：

ηexpη2erfη=cp，vTwall-Tsat?πhlv（19）

該解析解作為精確解用于與模擬結(jié)果對比，作為相變模型誤差分析的參考。

在本研究中，模擬的工質(zhì)為水，飽和溫度為373.15K，熱壁溫度為373.15K。計算域長度為0.8mm，總模擬時長為1s。模擬中所使用工質(zhì)的相關(guān)物性參數(shù)如下表所示。

兩種相變模型時，相界面位置隨時間的變化以及與精確解的對比。曲線顯示數(shù)值結(jié)果與精確解吻合良好。

4.2?Stefan冷凝問題

圖3為一維Stefan冷凝問題的示意圖。初始階段整個計算域?yàn)轱柡驼羝蟊诒３值陀陲柡蜏囟鹊暮愣囟?，受左壁影響，蒸汽開始冷凝，液體層厚度逐漸增加。對于一維Stefan冷凝問題，其解析解為：

δt=?2tλlρlcp，l12+hlvcp，lTsat-Twall-1（20）

模擬中左側(cè)為冷壁，溫度為363.15K。其余計算參數(shù)及物性參數(shù)與4.1中相同。

圖4分別顯示了使用兩種相變模型時，相界面位置隨時間的變化以及與精確解的對比。結(jié)果顯示，Rattner模型與精確解吻合較好，Lee模型在初始階段冷凝速率偏大，導(dǎo)致液層厚度偏大，之后逐漸與精確解吻合。

結(jié)語

本文對采用VOF方法模擬氣泡生長時所使用的相變模型進(jìn)行了概述，總結(jié)其中十分重要的相變模型，并對模型進(jìn)行了評析與對比。結(jié)果發(fā)現(xiàn)，Sato模型較為精確，但其尚未應(yīng)用到后續(xù)研究所需使用的軟件OpenFOAM中，Lee模型和改進(jìn)后的Rattner模型在選擇合適的參數(shù)時均可得到較為準(zhǔn)確的結(jié)果。但從實(shí)際出發(fā)，Lee模型的蒸發(fā)冷凝傳質(zhì)強(qiáng)度因子應(yīng)當(dāng)對于不同工況甚至不同網(wǎng)格取不同的值，取合適的值較為困難。因此更推薦在后續(xù)的研究中選擇Rattner模型。

參考文獻(xiàn)：

［1］Nabil?M，Rattner?A?S.InterThermalPhaseChangeFoamA?framework?for?twophase?flow?simulations?with?thermally?driven?phase?change［J］.Softwarex，2016，5（C）：216226.

［2］Brackbill?J?U，Kothe?D?B，Zemach?C.A?continuum?method?for?modeling?surface?tension［J］.Journal?of?Computational?Physics，1992，100（2）：335354.

［3］Lee?W?H.A?PRESSURE?ITERATION?SCHEME?FOR?TWOPHASE?FLOW?MODELING［M］.Washington，USA：Hemisphere?Publishing，1980.

［4］Sato?Y，Nieno?B.A?New?Conservative?Phase?Change?Model?for?Nucleate?Boiling［C］//.20th?International?Conference?Nuclear?Energy?for?New?Europe?2011，2011.

［5］Onishi?H，Kawamura?M，Tada?Y，et?al.Numerical?Analysis?on?Heat?Transfer?Characteristics?of?Looped?Minichannel?Using?PhaseChange?VOF?Method［C］//.Asme?International?Conference?on?Nanochannels，2013.

［6］Rattner?A?S，Garimella?S.Simple?Mechanistically?Consistent?Formulation?for?VolumeofFluid?Based?Computations?of?Condensing?Flows［C］//.Asme?International?Mechanical?Engineering?Congress?&?Exposition，2014.

作者簡介：王子威（1998—?），男，漢族，江蘇宿遷人，碩士研究生，研究方向?yàn)榉磻?yīng)堆熱工水力。