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

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

0 引言

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

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

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

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

1.1 控制方程

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

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

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

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

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

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

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

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

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

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

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

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

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

1.2 模型驗(yàn)證

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

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

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

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

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

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

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

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

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

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

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

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

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

2 幾何模型

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

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

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

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

3 結(jié)果討論

3.1 單個(gè)含不凝氣蒸汽氣泡運(yùn)動(dòng)特征分析

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

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

3.2 單個(gè)含不凝氣蒸汽氣泡冷凝行為分析

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

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

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

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

3.3 多氣泡分析

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

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

4 結(jié)論

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

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

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

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