基于OpenFOAM的含不凝氣蒸汽氣泡冷凝模擬

陳家?guī)?曾義凱 都 宇

（1.西南交通大學(xué)機械工程學(xué)院 成都 610031；2.中國核動力研究設(shè)計院中核核反應(yīng)堆熱工水力技術(shù)重點實驗室 成都 610213）

0 引言

直接接觸冷凝具有較高的傳熱效率，在很多工業(yè)過程中都能看到他的身影，比如天然氣鍋爐煙氣余熱回收以及沸水堆核電站抑壓水池。近年來，人們對于直接接觸冷凝進行大量的研究[1-3]。在天然氣余熱回收領(lǐng)域，由于氣泡的冷凝行為直接影響冷凝熱水器和冷凝鍋爐等設(shè)備對煙氣中余熱的回收效率，因此在此類煙氣換熱器的研發(fā)設(shè)計中，直接接觸冷凝一直被視作核心問題而被研究[4-6]。含不凝結(jié)氣體的蒸汽氣泡表面發(fā)生傳熱和傳質(zhì)現(xiàn)象的同時，在氣泡內(nèi)部還存在蒸汽-不凝結(jié)氣體組成的擴散層，因此其運動行為和傳熱特性十分復(fù)雜，對直接接觸冷凝的研究是一個難度很大的挑戰(zhàn)。為了能夠更深入的了解含有不凝結(jié)氣體的蒸汽氣泡直接接觸冷凝現(xiàn)象，有必要對其冷凝行為進行深入的研究。

至今，已經(jīng)有許多學(xué)者對于氣泡冷凝行為進行了大量的實驗研究分析。在屈曉航[7]的研究中，通過對實驗數(shù)據(jù)的擬合得到了冷凝Nu 數(shù)的關(guān)聯(lián)式。在Kim 和Park[8]的研究中，通過實驗分析得到低壓條件下過冷沸騰中純蒸汽氣泡表面的換熱系數(shù)關(guān)聯(lián)式。隨著計算機技術(shù)的發(fā)展和數(shù)值模擬方法的改進，越來越多的研究者采用CFD 來研究流動和換熱類的問題[9,10]，兩相流的模擬得到了進一步的發(fā)展。很多研究者采用VOF 方法對蒸汽氣泡的冷凝過程進行數(shù)值模擬，并從氣泡大小、氣泡形狀、上升軌跡和氣泡速度等方面分析氣泡的行為[11,12]。Welch 和Wilson[13]、屈曉航[14]和Samkhaniani[15,16]等人的研究中檢驗了采用VOF 模型分析過冷液中氣泡冷凝的可行性。

本文主要研究了如何利用開源軟件包OpenFOAM 中的程序模擬空氣-水蒸氣混合氣泡的直接接觸冷凝過程。本文對OpenFOAM[17]中的多相流多組分求解器icoReactingMultiphaseInterfoam進行了改進以對含有不凝結(jié)氣體的蒸汽氣泡進行數(shù)值模擬。將數(shù)值模擬結(jié)果與實驗結(jié)果進行了比對，驗證了該求解器的適用性。通過該求解器，研究了不同初始條件下空氣-水蒸氣混合氣泡的冷凝特性，并與純蒸汽氣泡的冷凝特征對比，分析不凝結(jié)氣體對冷凝氣泡行為的影響。

1 數(shù)值模型與驗證

1.1 控制方程

在本文中，由Hirt 和Nichols 提出Volume of Fluid（VOF）界面追蹤方法[18]將被用來對含有干空氣的蒸汽氣泡在過冷水中冷凝這個兩相三組分系統(tǒng)中的氣-液相界面進行追蹤。VOF 方法的控制方程如下：

式中，αL= 1表示的是液相區(qū)域，αL= 0表示的是氣相區(qū)域，0 ＜ αL＜ 1表示氣相和液相交界的區(qū)域。氣泡的邊界位于氣相和液相交界的區(qū)域。在實際工程應(yīng)用中，通常將αL= 0.5的等值面作為氣相和液相的交界面。

隨著流體的流動，氣液交界面的位置隨之發(fā)生改變，既流場中氣相體積分數(shù)分布應(yīng)該同步發(fā)生改變，需要對相體積分數(shù)輸運方程進行求解。對于存在相變的多相流系統(tǒng)，其對應(yīng)的體積分數(shù)方程為：

式中，m&表示單位體積冷凝率，kg/(m3·s)。數(shù)值擴散的存在導(dǎo)致求解上述方程得到的相體積分數(shù)場中兩相的交接區(qū)域過大，相界面不夠尖銳。前人已提出多種可用于VOF 模型中以保證界面尖銳的方法。在OpenFOAM 中采用Weller[19]提出的在體積分數(shù)輸運方程中添加人工對流項的方法，來對兩相界面附近的相分數(shù)進行壓縮，削弱數(shù)值耗散的影響，保證求得相界面足夠銳利。這個額外的對流項為 ?· （α1（1 - α1）Uc ），其只在兩相的交接區(qū)域起作用，在純液相區(qū)域和純氣相區(qū)域該項的值為零。加入人工對流項后體積分數(shù)輸運方程變?yōu)椋?/p>

Cα表示壓縮系數(shù)，當壓縮系數(shù)的取值在1 到4 之間時，對于大多數(shù)工況都能得出較好的結(jié)果，但是對于某些工況壓縮系數(shù)的值需要大于4，本文中的壓縮系數(shù)Cα取值為1。

為了簡化邊界條件的定義，在OpenFOAM 中動量方程的總壓（P）轉(zhuǎn)換為了動壓（Prgh）。求解器中求解的動量方程如下：

等式左邊的最后一項表示氣相和液相之間的表面張力。表面張力通過連續(xù)表面應(yīng)力模型[20]計算。σ表示表面張力系數(shù)，k 表示氣液交界面的曲率。

由于相變伴隨能量傳遞，需要求解能量方程。在OpenFOAM 的多相流求解器中，為了削弱相界面兩側(cè)氣相和液相物性參數(shù)的不平衡性，將能量方程表示為方程（7）所示的形式，以增強數(shù)值算法的魯棒性。

式中，HLG表示單位質(zhì)量蒸汽冷凝所釋放的相變潛熱，J/kg。

氣泡中水蒸氣在干空氣中的擴散，由如下方程描述：

等式右邊第二項是組分擴散方程的源項，既蒸汽的冷凝量。由于混合氣體中只有水蒸氣發(fā)生冷凝，水蒸氣的冷凝量即是氣相質(zhì)量的變化量。除了不凝氣-蒸汽混合氣體密度Gρ 是由理想氣體公式計算，其余的物性參數(shù)根據(jù)蒸汽和干空氣的質(zhì)量分數(shù)加權(quán)求得。

為了閉合控制方程，需要采用合適的相變模型計算相界面出的質(zhì)量通量。本文中采用的相變模型是Lee 模型，其對應(yīng)的數(shù)學(xué)方程為：

式中，rc是一個經(jīng)驗系數(shù)稱為傳質(zhì)強度因子，通常是通過和實驗數(shù)據(jù)對比得到。本文通過和已有的實驗數(shù)據(jù)對比將其取值為500。

1.2 模型驗證

為了驗證本文所建立的數(shù)值模型的合理性，本文對屈曉航等人的實驗工況進行數(shù)值模擬，并將數(shù)值模擬結(jié)果與實驗結(jié)果比對。對含有不凝結(jié)氣體的蒸汽氣泡冷凝的模擬流程如圖1 所示，具體步驟如下：

圖1 直接接觸氣泡冷凝模擬流程Fig.1 Simulation flow of direct contact bubble condensation

第一步：在三維笛卡爾坐標系下，采用blockMesh 生成幾何模型和網(wǎng)格。如圖1 中所示，創(chuàng)建的計算域的尺寸為30×30×60mm3氣泡具有向上的初始速度。計算域的上邊界設(shè)置為壓力出口邊界，下邊界和側(cè)面則設(shè)置為壁面邊界。

第二步：采用OpenFOAM 中自帶的小工具setFields 對計算域中各個物理量進行初始化。對于氣相區(qū)域的蒸汽質(zhì)量分數(shù)、溫度和速度等物理量的初始條件根據(jù)實驗中測得的數(shù)據(jù)進行設(shè)置。所要模擬的兩個工況的初始條件如表1 中所示，氣泡1 對應(yīng)的過冷度為14K，氣泡2 對應(yīng)的過冷度為10K。

表1 含不凝氣蒸汽氣泡的初始條件Table 1 Initial conditions of bubbles air-vapor mixture bubble

第三步：為了保證結(jié)果的合理性，需要保證氣相區(qū)域有足夠數(shù)量的網(wǎng)格，并且為了避免數(shù)值擴散對相界面的過度涂抹往往需要保證界面區(qū)域的網(wǎng)格尺寸小于氣泡直徑的1/16。本文為以較小的數(shù)值消耗計算出較為尖銳的氣液兩相界面，在界面區(qū)域采用了自適應(yīng)網(wǎng)格技術(shù)。

第四步：在三維笛卡爾坐標系下，采用本文改進后的icoReactingMultiphaseInterfoam 求解器對單個含有不凝結(jié)氣體的蒸汽氣泡的冷凝過程進行模擬計算。

圖2 定性的展示了在兩組不同的初始工況下，數(shù)值模擬結(jié)果和實驗結(jié)果的對比，包括冷凝過程中氣泡形狀，大小以及穿透距離等。氣泡在上升的過程中隨著氣泡中蒸汽的冷凝，其體積減小、外形也隨之變化。如圖2 中所示，對于直徑稍小的氣泡，其外形由圓球形變?yōu)榘肭蛐谓又俎D(zhuǎn)變?yōu)闄E圓形，對于初始直徑較大的氣泡，其外形由初始時的球形，隨后下表面內(nèi)凹，接著轉(zhuǎn)變?yōu)榘肭蛐危詈筠D(zhuǎn)變?yōu)闄E球形。這是由于，隨著氣泡直徑的增大，表面張力對氣泡變形的影響力減弱。而最后都變?yōu)闄E球形則是由于受力的原因。由于數(shù)值模擬的初始條件和實驗的真實條件有區(qū)別，故而數(shù)值模擬的結(jié)果和實驗結(jié)果在氣泡外形、氣泡穿透距離等方面有細微的差異。但是鑒于實驗中初始條件的不確定性，可以認為本文的數(shù)值模型計算得到的結(jié)果和實驗結(jié)果是一致的。

圖2 數(shù)值模擬和實驗的結(jié)果對比Fig.2 Comparison of numerical simulation and experimental results

圖3 定量展示了氣泡的體積隨時間變化的規(guī)律。數(shù)值模型計算出的體積變化曲線和實驗數(shù)據(jù)基本吻合。對于含不凝結(jié)氣體的蒸汽氣泡，隨著蒸汽在氣泡表面處凝結(jié)，氣泡中不凝氣體含量增高，組分擴散層中的濃度梯度越發(fā)明顯，傳熱和傳質(zhì)阻力也越發(fā)的大，最終會達到穩(wěn)定的狀態(tài)，氣泡的體積基本維持穩(wěn)定。由圖2 和圖3 展示的結(jié)果可看出，本文建立的數(shù)值模型中的質(zhì)量源項，貼合實際相變過程。

圖3 氣泡體積變化對比Fig.3 Comparison of bubble volume changes

圖4 為圖2 中冷凝氣泡的橫截面處溫度分布和蒸汽質(zhì)量分數(shù)分布。在圖4 中紅實線外的區(qū)域，紅實線和紅虛線之間的區(qū)域，以及紅虛線以內(nèi)的區(qū)域分別表示過冷水區(qū)，蒸汽-干空氣的組分擴散層以及氣泡內(nèi)部蒸汽聚集的區(qū)域。如圖所示，含不凝結(jié)氣體的蒸汽氣泡中溫度連續(xù)分布，而純蒸汽氣泡在冷凝過程中的溫度分布在相界面處有明顯的間斷。這是由于在蒸汽-干空氣混合的氣泡中，從氣泡中心到氣泡邊界，水蒸氣的質(zhì)量分數(shù)不斷的下降，水蒸氣的分壓也隨之下降，進而導(dǎo)致對應(yīng)的飽和溫度也逐漸降低，且低于氣泡中心的飽和溫度。

圖4 含不凝結(jié)氣體蒸汽氣泡的溫度分布和質(zhì)量分數(shù)分布Fig.4 Temperature distribution and mass fraction distribution of air-vapor bubble

2 幾何模型

為了深刻理解含不凝氣蒸汽氣泡冷凝過程中各種因素的影響，本文在直角坐標系下，對圖5 中所示的二維氣泡進行了數(shù)值模擬。在計算域底部釋放含不凝氣蒸汽氣泡，研究各種參數(shù)對其冷凝行為的影響。情形一，如圖5（a）所示，研究不凝結(jié)氣體含量、氣泡初始直徑、流場速度以及過冷度等因素對氣泡冷凝行為的影響。情形二，如圖5（b）和圖5（c）所示，將兩個氣泡在計算域底部并排或者豎排放置，研究相鄰氣泡對冷凝行為的影響。

圖5 干空氣-蒸汽混合氣泡冷凝模擬示意圖Fig.5 Schematic diagrams of air-vapor mixture bubble condensation simulations

為了求得不受網(wǎng)格密度影響的數(shù)值解，在五組不同的網(wǎng)格密度下，對圖5（a）中的工況進行了數(shù)值計算。圖6 展示了在不同的網(wǎng)格密度下，計算所得的氣泡等效直徑變化曲線以及氣泡在同一時刻的外形。通過對五組網(wǎng)格計算結(jié)果的比較，可以看出250×500 所對應(yīng)的網(wǎng)格密度計算出的結(jié)果是較為合理的，因此其被用于本文的后續(xù)計算。

圖6 網(wǎng)格無關(guān)性驗證Fig.6 Obtaining a grid independence solution

3 結(jié)果討論

3.1 單個含不凝氣蒸汽氣泡運動特征分析

為了研究流場速度以及不凝氣含量對于遷移軌跡的影響，在過冷度為15K 的條件下，對初始直徑為8mm 的具有不同的初始蒸汽質(zhì)量的氣泡在不同的流場速度下模擬。如圖7 中所示，圖中對五種工況進行了計算，（a）、（b）、（c）、（d）、（e）對應(yīng)的初始蒸汽質(zhì)量和流場速度分別為0.9、0.7、0.7、0.6、0.6 和0m/s、0m/s、0.1m/s、0.1m/s、0.2m/s。對每個工況每間隔0.01s 提取一次氣泡的邊界坐標。對于含不凝結(jié)氣體的蒸汽氣泡，在相同的過冷度下，不凝結(jié)氣體的含量主要影響氣泡體積基本穩(wěn)定后氣泡的變形，對于氣泡遷移軌跡影響甚微。流場速度對于氣泡的遷移軌跡影響較大，但是增大流場速度不能夠加快氣泡中的蒸汽凝結(jié)。

圖7 氣泡遷移軌跡Fig.7 Bubble migration trajectory

3.2 單個含不凝氣蒸汽氣泡冷凝行為分析

對于含不凝氣蒸汽氣泡直接接觸冷凝，需要注意其完成冷凝需要的時間，本文將氣泡體積基本不變視為完成冷凝，此時氣泡對應(yīng)的等效直徑稱為終端等效直徑。冷凝時間受到氣泡初始尺寸，氣泡中初始蒸汽質(zhì)量分數(shù)以及過冷度等因素的影響。由圖8（a）可知，隨著直徑的增大，氣泡的冷凝時間增大，但是氣泡的冷凝時間與氣泡初始直徑之間不呈現(xiàn)明顯的線性關(guān)系，這點與純蒸汽氣泡不同。隨著過冷度的增大，冷凝時間隨之增大，并且隨著氣泡直徑的增大，在一定的過冷度范圍內(nèi)，氣泡的冷凝時間分布范圍減小。如圖8（b）所示，終端等效直徑隨過冷度的增大而減小，隨氣泡初始直徑的增大而增大，并且與氣泡的初始直徑呈線性關(guān)系。此外，在一定的過冷度范圍內(nèi)，隨著氣泡直徑的增大，其對應(yīng)的終端等效直徑的范圍也隨之擴大。已有文獻指出，純蒸汽氣泡的冷凝時間與氣泡的初始直徑在一定程度上呈現(xiàn)線性關(guān)系，而對于含不凝氣蒸汽氣泡，不同初始直徑或者不同初始蒸汽質(zhì)量分數(shù)的氣泡，其在冷凝過程中形成的不凝氣-蒸汽組分擴散層也隨之出現(xiàn)差異，進而對蒸汽冷凝阻礙程度不同，導(dǎo)致氣泡的冷凝時間與氣泡尺寸和不凝結(jié)氣體的含量呈非線性關(guān)系。

圖8 過冷度對含不凝氣蒸汽氣泡冷凝的影響Fig.8 Influence of subcooling degree on air-vapor bubble condensation

圖9（a）中展示了過冷度為15K，初始直徑為6mm 的不同蒸汽質(zhì)量分數(shù)的氣泡在冷凝過程中的體積變化規(guī)律。圖9（b）則展示了不凝氣體含量以及氣泡初始直徑對于氣泡終端等效直徑的影響。隨著氣泡中不凝氣體含量的增加，氣泡的體積衰減速率減小，終端等效直徑增大。并且，隨著氣泡初始直徑的增大，在一定的過冷度范圍內(nèi)，氣泡的終端等效直徑的分布范圍增大。不凝氣體含量大的蒸汽氣泡隨著初始直徑的增大，對應(yīng)的終端等效直徑增長的更快。

圖9 不凝氣含量對冷凝的影響Fig.9 Influence of non-condensable gas content on condensation

3.3 多氣泡分析

圖10 展現(xiàn)了兩個并列和豎排的含不凝氣蒸汽氣泡同時冷凝時，相互之間的影響。圖中列出了四中工況的模擬結(jié)果。情形一，模擬豎排的氣泡對的冷凝過程，氣泡中心之間的間距為1.1D0，如圖5（c）所示。情形二和情形三分別是情形一中兩個氣泡的冷凝過程的單獨模擬，以和情形作對照，探究相鄰氣泡對冷凝過程的影響。情形四是模擬并列的氣泡對的冷凝過程，氣泡中心的間距仍為1.1D0，如圖5（b）所示。情形二中氣泡在0.05s 對應(yīng)的等效直徑為5.223mm，情形三中氣泡在0.05s 對應(yīng)的等效直徑為5.175mm。這是由于情形三中氣泡所處的深度較情形二深，其氣泡內(nèi)部的壓力較大，進而氣泡內(nèi)蒸汽對應(yīng)的飽和溫度越高，冷凝加快。另外，綜合對比情形一、情形二和情形三，位于下側(cè)的氣泡的蒸汽冷凝率要低于其單氣泡時的冷凝率。當氣泡對豎直放置時，且氣泡間隔較小時，下部氣泡位于上部氣泡的熱邊界層中，氣泡中的蒸汽溫度與氣泡外側(cè)流體溫度的差值減小，蒸汽冷凝減緩。由情形四可知，當兩個含不凝氣蒸汽氣泡并排放置時，在冷凝過程中受到來自過冷液的力的作用將向兩側(cè)旋轉(zhuǎn)分離。

圖10 多氣泡冷凝過程Fig.10 Multi-bubble condensation process

4 結(jié)論

（1）含不凝氣蒸汽氣泡在冷凝過程中的變形過程為：球形，月牙形，半球形，橢圓形，當氣泡的體積不再變化后，氣泡的形狀會發(fā)生震蕩。對于初始直徑較小的含不凝氣蒸汽氣泡，其在冷凝過程中不會出現(xiàn)下表面內(nèi)凹。氣泡邊界處存在不凝氣-蒸汽組分擴散層，存在明顯的溫度梯度。

（2）氣泡中不凝氣含量決定了氣泡完成冷凝后的體積，進而直接影響氣泡的變形，對氣泡的移動軌跡影響相對較小。而流場速度主要影響氣泡的遷移軌跡。

（3）含不凝氣蒸汽氣泡表面存在蒸汽冷凝的時長隨過冷度的增大而增大，比如對于初始蒸汽質(zhì)量分數(shù)為0.7，初始直徑為6mm 的氣泡，當過冷度為5K 時，冷凝時長為50ms，而當過冷度為20K時，冷凝時長為60ms。此外，含不凝結(jié)氣體的蒸汽氣泡的冷凝時長與氣泡初始直徑以及過冷度均呈現(xiàn)非線性關(guān)系，而終端等效直徑則是與氣泡的初始直徑呈線性正相關(guān)，并隨過冷度的增大而減小。初始不凝氣含量較多的氣泡對應(yīng)的蒸汽冷凝率較小，其對應(yīng)的終端等效直徑隨氣泡初始直徑的增大而增長較快。

（4）當兩個含不凝氣蒸汽氣泡豎直排列，且間隔較小時，上面的氣泡加熱的液體位于下面的氣泡的上部，進而對下面氣泡的冷凝會起到一個抑制，故而對于利用直接接觸氣泡冷凝來進行換熱的換熱器，要保證生成的氣泡之間有足夠的間隔。