VOF模型界面傳質與體積傳質的轉換方法

譚思超，趙富龍，李少丹，高璞珍

（1.哈爾濱工程大學核安全與仿真技術國防重點學科實驗室，黑龍江哈爾濱150001；2.清華大學核能與新能源技術研究院，北京100084）

VOF模型界面傳質與體積傳質的轉換方法

譚思超1，趙富龍2，李少丹1，高璞珍1

（1.哈爾濱工程大學核安全與仿真技術國防重點學科實驗室，黑龍江哈爾濱150001；2.清華大學核能與新能源技術研究院，北京100084）

針對流體體積模型（volume of fluid，VOF）在模擬相間傳質過程時需要將界面質量流密度轉換為單位體積傳質速率的問題，對VOF模型中的界面傳質與體積傳質轉換方法進行了改進。提出一種既解決了網格無關性問題又可以反映局部界面傳質特性的轉換方法，并給出了理論推導證明。推導出相應的網格無關性條件，通過劃分3組不同尺寸的網格模擬了汽泡冷凝問題，并在汽泡生長方面對該轉換方法加以簡單應用，模擬結果與理論分析及實驗結果符合良好，合理可行，可廣泛應用到多相流中界面傳熱傳質模擬的過程。

VOF模型；界面；局部傳質；理論推導；網格無關性；多相流

過冷沸騰由于其換熱效果強烈，能較好地滿足核工業體積小、能量密度高的要求，受到越來越多研究者的關注；窄通道由于結構緊湊、換熱性能好等特點，得到越來越多的重視。為了較好地研究強化換熱特性，大量研究者對窄通道內的汽泡行為進行了研究，大部分是基于實驗研究得出相關結論。汽泡行為的研究中，數值模擬是一種重要的手段，其中VOF模型憑借較好的界面追蹤特性得到越來越多研究者的青睞，Jeon等［1?6］分別采用VOF模型對過冷沸騰下汽泡的生長及冷凝行為進行數值模擬。

用VOF模型進行相變過程計算時，采用體積傳質速率進行計算，但一般已知多為傳質速率或者質量流密度，因此必然涉及到界面傳質速率與體積傳質速率間的轉換。目前，在使用VOF模型進行界面傳熱傳質特性計算時，主要有蒸發冷凝模型和平均處理方法。蒸發冷凝模型，是Lee［7］提出的體積傳質處理方法，主要思想是，溫度高于液相的飽和溫度時，液相蒸發，對應到相應的網格即為液相向汽相傳質；反之，汽相冷凝。該模型形象直觀，易于理解，但結果對網格尺寸有很強依賴性，假定界面層厚度不變，由于網格體積與尺寸之間為三次方關系，若劃分的網格越細，則界面上的網格數量呈三次方變化。該模型單個網格的體傳質速率保持不變，會使得相間界面總傳質速率隨網格數量呈現指數增加，在網格無關性方面存在一定的缺陷。平均處理方法是Jeon［1］在用VOF模型對過冷流動沸騰條件下汽泡冷凝行為模擬時提出的，進行汽泡冷凝傳質的計算時，先計算汽液界面總傳質量，隨后在每個網格進行平均。盡管考慮了網格尺寸的作用，解決了網格無關性問題，但由于涉及平均處理過程，不能很好地反映汽液界面的局部傳質過程。

綜上所述，現有VOF模型傳熱傳質計算方法均存在一定的不足。為此，在理論分析的基礎上提出一種VOF模型的界面傳質與體積傳質轉換的處理方法，用于計算涉及多相相間界面傳熱傳質問題，以便更好的用于VOF模型相變過程的計算。

1 汽泡描述方法及傳熱傳質模型

1.1 汽泡描述方法

VOF［8］通過各相的體積份額來描述不同的流體相，對于汽液兩相流，單元中汽相體積份額αv為0代表液相，αv為1代表汽相，介于兩者之間為汽液界面；在固定的歐拉網格下，求解各相體積份額的連續性方程來描述2種或者多種互不混合流體界面；各相共享單一動量、能量方程，各單元的密度、粘度等參數采用體積份額進行加權平均求得。

文中采用幾何重構方法描述汽液界面，即采用分段線性方法描述汽液界面，假定兩相間界面在每個單元內有一個線性斜面，使用此線性形狀計算穿過單元面的對流及位置信息，示意圖如圖1，白色弧線代表汽液界面，弧線左上方黑色區域代表液相，弧線右下方區域代表汽相。

圖1 幾何重構方法描述的汽液界面形狀Fig.1 Liquid?vapor interface shape respresented by the geometric reconstruction scheme

按照幾何重構的方法，對汽泡的描述也是通過體積份額來刻畫，汽液界面內部即為典型的汽泡，如圖2所示。

圖2 VOF模型刻畫的汽泡剖面圖Fig.2 Bubble section map by VOF model

需要指出的是，典型汽泡的刻畫是需要多個網格來共同組合描述的，Rabha等［9］認為網格數／汽泡直徑（cells per buddle diameter，CPD）應至少為16，即當汽泡直徑與網格尺寸比值大于16時才能較好的刻畫汽泡，下文提出的體積轉換方法的網格要求也在這一范圍內，說明了改進方法的網格無關性條件符合VOF模型描述汽泡的基本要求。

1.2 汽液界面傳熱傳質模型及UDF

過冷流動沸騰汽泡行為方面，目前公認的汽泡生長動力為微液層蒸發與汽液界面蒸發冷凝綜合作用的結果［10］，本文采用經典的體現界面傳質特性的Hertz?Knudsen公式來進行汽液界面傳質過程計算［11］，其表達式如下：

式中：θ為界面蒸發／冷凝系數；R為氣體常數，J／mol-1·K-1；M為摩爾質量，kg／mol；Tl為液相溫度，K；ps為Tl對應的飽和壓力，MPa；pv為汽泡內壓力，MPa；Tv為汽相溫度，K。

根據上述數學模型編寫UDF程序，加入相應的質量源項和能量源項，能量源項是在質量源項的基礎上考慮汽化潛熱的作用得到的。

2 體積傳質轉換方法

2.1 基本假設

在理論推導之前作以下幾點假設：1）汽泡為已從加熱壁面脫離位于管道中央的孤立汽泡；2）所劃分的網格為正方體結構性網格，各個網格的邊長、表面積、體積均相等；3）傳質過程為均勻傳質過程；4）汽液界面厚度均勻；5）在計算過程中假設汽泡內為飽和狀態。

汽泡內壓力由Young?Laplace公式來確定：

式中：pv和pl分別為汽泡內外壓力，σ為液體表面張力，汽泡的半徑R是通過在計算過程中測定汽泡在x、y、z 3個方向的高度：

2.2 理論推導

VOF模型進行多相流計算時，汽液界面是通過多個網格間界面重構方法構建的。每個網格內都充滿流體，若流場中只含有汽液兩相，則

式中：αvi、αli分別為第i個網格中汽相和液相的體積份額。

從式（4）可知，用VOF模型對兩相流進行計算，蒸發與冷凝相當于2個逆向過程，即在VOF模型中適合蒸發過程的算法也適用于冷凝過程。若過程為蒸發過程，則相間傳質界面面積可近似為αviSi，其中Si為單個網格截面積。采用這種近似方法進行處理，是因為VOF模型各相靠體積份額描述的，汽液界面處的網格單元中既有汽相又有液相，相間有一定的接觸面，傳質過程就是在接觸面上進行的。已知網格表面積，若已知傳質界面面積與網格表面積的比例關系，便可得到傳質界面面積，但在VOF模型中傳質界面的面積不能通過上述方法得到，能得到的物理量只有各相的體積份額，為此采用體積份額對相間傳質界面面積近似處理，即認為相間傳質界面的面積約為αviSi。

按上述近似處理方法，則單個網格的傳質速率近似為αviGiSi，汽相所占的體積為αviVi，其中Vi為網格的體積，根據VOF模型中的體傳質速率的定義可得到單個網格的體傳質速率mi為

式中：Gi為質量流密度，Li為網格尺寸。式（5）即為本文中提出的將質量流密度轉換為體傳質速率的處理方法。下面在理論推導基礎上對提出的轉換方法進行網格無關性的分析。

汽液界面上總的傳質速率Mt表達式為

結合式（6）、（7）算汽液界面總傳質速率Mt：

從式（8）可以看出，當汽泡尺寸一定時，要以一定的精確度ε確保汽液界面上總的傳質速率Mt與網格尺寸無關，則需滿足下面的條件：

即

式（10）即為改進方法的網格無關性條件。可以看出要保證計算精度為10%，則要求汽泡半徑約為網格尺寸的10倍，下文用此精度進行方法驗證。

3 方法驗證

3.1 幾何模型及方法驗證

窄通道一般是指間隙尺寸相對很小，進出口長度與常規通道相當的通道。Kandlikar［12］認為當量直徑在200 μm～3 mm的通道稱為窄通道或小通道。而根據工程實際，通常將當量直徑在1～3 mm的通道稱為窄縫通道。考慮到進行三維數值計算網格量及計算精度要求，本文選取窄矩形通道幾何模型的尺寸為2 mm× 4 mm×2 mm，流體豎直向上流動，重力方向豎直向下。劃分3組不同疏密程度的六面體結構性網格：1）50× 100×50，2）64×128×64，3）80×160×80，相應的六面體網格的邊長分別為0.04、0.031 25、0.025 mm，并對3組不同網格尺寸的汽泡的冷凝速率及冷凝過程中的汽泡形狀圖像進行對比。通過對汽泡冷凝過程的三維數值模擬結果分析，對提出的轉換方法進行驗證。

邊界條件的選擇，進口邊界條件用速度入口，為穩定工況，出口邊界條采用壓力出口。離散方法的選擇，考慮到求解速度、穩定性、精度等因素，采用雙精度進行計算，壓力速度耦合采用PISO算法，使用PRESTO方法離散壓力，動量和能量方程使用二階迎風格式求解，瞬態方程的離散采用一階迎風格式求解，容積比率方程用幾何重構方法求解，非穩態時間步長為5 μs。

3.2 模擬結果及工況

為了證明模型的準確性，參考潘良明等［3］的實驗工況進行了模擬，所選擇的計算條件如下表1所示。選擇該文獻進行對比主要是因為其采用VOF模型對過冷流動沸騰條件下汽泡冷凝過程進行數值模擬，其模擬結果與相應條件下的實驗結果吻合較好，誤差在±20%之內，表明能夠反映真實參數變化過程。

表1 方法驗證模擬工況參數Table1 Simulation parameters for method verification

計算過程中，所有的液相參數根據過冷水的溫度和相應的壓力查得，汽相參數根據相應的飽和溫度查得，隨著冷凝過程的進行汽泡的直徑會逐漸減少，為保證計算精度，同時驗證網格無關性條件，需要對網格進行實時加密，保證網格與汽泡半徑的比值與初始時刻的比值相同。所得到的冷凝過程中的3組不同疏密程度網格對應的汽泡冷凝曲線如下圖3所示，其中，1、2、3分別代表上述3組網格對應的模擬結果，上標“'”表示文獻［3］的實驗結果。b工況冷凝過程中汽泡的形狀變化圖像如圖4所示，在圖4中（a）、（b）、（c）網格結構分別為：50×100×50、64×128×64、80×160×80。

圖3 汽泡冷凝曲線Fig.3 Bubble condensation curves

圖4 模擬得到的汽泡生長過程形狀Fig.4 The simulated shapes during bubble growth

從圖3中，可以看出本文的計算結果與PAN等［5］的實驗結果趨勢基本一致，冷凝時間差別不大，氣泡直徑最大誤差為14.8%，表明本文所提出的界面傳質與體積傳質的轉換方法準確。從圖3可以看出，汽泡冷凝過程中，不同疏密網格計算的汽泡直徑變化吻合良好，最大誤差在±6.3%以內。對于b工況，通過對圖3中“b?1，b?2，b?3”汽泡半徑隨著時間的變化曲線的對比，可以看出不同疏密的網格計算結果吻合良好，冷凝過程一致，說明了上述的VOF模型中汽液兩相相間傳質的轉換方法中網格無關性條件正確。

圖3中的模擬結果表現出線性特性，主要是因為數值模擬中的相間傳熱傳質過程的實現是按照經驗公式人為編寫程序加入的，而現有的經驗公式考慮的因素有限，最終會導致結果變現為一定的線性關系；另外，為了保證文中的網格無關性條件，在方法驗證過程中，汽泡并沒有完全冷凝，如圖3中的曲線所示，要想進一步計算則需要劃分尺寸更小的網格。

方法驗證雖是在均勻傳質假設下進行的，但同樣適用于非均勻傳質過程，如過冷流動沸騰，雖涉及微液層蒸發和經典界面傳質公式2種模型作用，但可將整個汽泡劃分為2個部分球缺形，分別進行單獨計算，同樣適用于上述的論證。對于非結構性網格，可根據劃分非結構性網格的算法，通過計算給出相應的權重系數來進行分配。

4 方法應用

在轉換方法驗證的基礎上，加以簡單應用，將轉換方法應用到豎直窄矩形通道中過冷沸騰條件下汽泡在熱力生長特性的三維模擬仿真。

4.1 數學模型

過冷流動沸騰下汽泡生長動力一方面來自汽液界面的傳熱傳質，另一方面汽泡生長過程中微液層不斷蒸發也會成為汽泡生長動力。為此，在計算中加入界面傳熱傳質模型及微液層蒸發模型，界面傳熱傳質模型采用式（1）來計算。微液層蒸發模型，采用Addlesee等［13］對滑移氣泡的流場分析得到的滑移汽泡微液層厚度表達式：

式中：vl、H、u分別為液相運動粘度、汽泡高度和滑移速度。則轉化的質量流密度為

式中：λ為液相導熱系數，hγ為汽化潛熱。

傳質量以及傳熱量通過提出的轉化方法轉化為UDF程序語言，加入到計算過程中。

4.2 結果分析

由于不能實現自主產生汽泡，需要初始給定汽泡，初始直徑為0.1 mm。考慮計算速度及精度要求選擇管道尺寸為1 mm×1 mm×1 mm，不同區域劃分不同粗細網格，既保證滿足網格無關性條件，又能適當減少計算量，經過測試網格數量為173 781。

工況參數為工作壓力0.101 325 MPa、液體過冷度5 K、壁面過熱度10 K、流速0.1 m／s。邊界條件及離散方法等與“方法驗證”一節中相同。計算得到的汽泡生長過程如圖4所示，重力方向豎直向下。

對汽泡生長曲線進行擬合，得到的擬合曲線見圖5，擬合曲線的相關性系數為0.994，表明模擬的汽泡生長遵循指數規律，與大量實驗得出的汽泡后期生長階段規律一致。Thorncroft等［14］對FC?87工質進行的大量實驗發現過冷沸騰下汽泡生長的指數介于1／3～1／2之間，本文數值模擬結果中指數為0.416 39，吻合較好，該指數n小于0.5，是因為初始給定汽泡尺寸，忽略了汽泡的慣性生長過程。

從圖4、5可以看出，本文汽泡傳質模型能較好的模擬汽泡熱力生長過程；但前期的慣性控制階段差別較大，是因為在計算時初始給定汽泡尺寸，不是經過自然核化點核化來產生汽泡，無法準確模擬汽泡慣性生長階段的行為。另外，由于文中更加關注兩相間的傳質過程，而對汽泡的受力考慮不是很充分，不能直觀看出汽泡的滑移運動。

上述簡單的應用，表明提出的VOF模型界面傳質與體積傳質轉換方法合理可行，可較好應用到沸騰條件下汽泡演化特性的模擬仿真過程中。

圖5 汽泡生長曲線Fig.5 The bubble growth curve

5 結束語

本文在理論分析基礎上，提出了一種適用于VOF模型的界面傳質與體積傳質的轉換方法，數值模擬結果與理論推導得到的結果吻合較好，誤差在±6.3%以內。該轉換方法可以通過理論推導出相應的網格無關性條件，解決了網格無關性問題又，并可以描述局部界面傳質特性，是較為適用于VOF模型的界面傳質與體積傳質轉換的處理方法。

通過對豎直窄矩形通道中過冷沸騰條件下汽泡在熱力生長階段行為的三維數值模擬，對方法進行了簡單的應用，與現有的方法相比精度較高，為今后多相相間界面相變過程提供了一種準確的數值仿真方法。

［1］JEON S，KIM S，PARK G.Numerical study of condensing bubble in subcooled boiling flow using volume of fluid model［J］.Chemical Engineering Science，2011，66（23）：5899?5909.

［2］魏敬華，潘良明，袁德文，等.過冷流動沸騰相變過程汽泡特性的VOF方法模擬［J］.核動力工程，2012，33（6）：65?71.WEI Jinghua，PAN Liangming，YUAN Dewen，et al.VOF simulation of bubble characteristics of subcooled flow boiling［J］.Nuclear Power Engineering，2012，33（6）：65?71.

［3］潘良明，譚智威，閆曉，等.窄流道內過冷流動沸騰汽泡凝結過程的數值模擬研究［J］.核動力工程，2011，32（6）：86?90. PAN Liangming，TAN Zhiwei，YAN Xiao，et al.Numerical investigation of bubble condensation of subcooled boiling in narrow rectangular channels［J］.Nuclear Power Engineer?ing，2011，32（6）：86?90.

［4］潘良明，譚智威.過冷流動沸騰汽泡凝結變形及流場特性的數值模擬［J］.重慶大學學報，2012，35（6）：53?57.PAN Liangming，TAN Zhiwei.Numerical investigation of bubble deformation and flow field characteristics of subcooled boiling during condensation［J］.Journal of Chongqing Uni?versity，2012，35（6）：53?57.

［5］PAN L，TAN Z，CHEN D，et al.Numerical investigation of vapor bubble condensation characteristics of subcooled flow boiling in vertical rectangular channel［J］.Nuclear Engi?neering and Design，2012，248：126?136.

［6］袁德文，潘良明，陳德奇.豎直窄流道內過冷流動沸騰的單汽泡生長模型［J］.化工學報，2009，60（11）：2723?2728.YUAN Dewen，PAN Liangming，CHEN Deqi.Model for sin?gle bubble growth of subcooled flow boiling in vertical narrow rectangular channel［J］.Journal of Chemical Industry and Engineering，2009，60（11）：2723?2728.

［7］LEE W H.Pressure iteration scheme for two?phase flow mod?eling［J］.Multiphase Transport：Fundamentals，Reactor Safety，Applications，1980，1：407?432.

［8］ANSYS Inc.FLUENT theory guide［Z］.Pittsburg：ANSYS Inc，2012.

［9］RABHA S S，BUWA V V.Volume?of?fluid（VOF）simula?tions of rise of single／multiple bubbles in sheared liquids［J］.Chemical Engineering Science，2010，65（1）：527?537.

［10］郭烈錦.兩相與多相流動力學［M］.西安：西安交通大學出版社，2002：323?410.

［11］王遵敬，陳民，過增元.蒸發與凝結現象的分子動力學研究［J］.西安交通大學學報，2001，35（11）：1126?1130.WANG Zunjing，CHEN Min，GUO Zengyuan.Molecular dynamics study on evaporation and condensation［J］.Jour?nal of Xi'an Jiaotong University，2001，35（11）：1126?1130.

［12］KANDLIKAR S G.Fundamental issues related to flow boil?ing in minichannels and microchannels［J］.Experimental Thermal and Fluid Science，2002，26（2?4）：389?407.

［13］ADDLESEE A J，CORNWELL K.Liquid film thickness a?bove a bubble rising under an inclined plate［J］.Chemical Engineering Research and Design，1997，75（7）：663?667.

［14］THORNCROFT G E，KLAUSNER J F，MEI R.An experi?mental investigation of bubble growth and detachment in vertical upflow and downflow boiling［J］.International Jour?nal of Heat and Mass Transfer，1998，41（23）：3857?3871.

The transformation method of mass flux and mass transfer rate per volume at the interface in VOF model

TAN Sichao1，ZHAO Fulong2，LI Shaodan1，GAO Puzhen1
（1.National Defense Key Subject Laboratory for Nuclear Safety and Simulation Technology，Harbin Engineering University，Harbin 150001，China；2.Institute of Nuclear and New Energy Technology，Tsinghua University，Beijing 100084，China）

The interface mass flux needs to be transferred to mass transfer rate per volume，when simulating mass transfer process between the phases by volume of fluid（VOF）model.In order to solve this problem，the method for transformation from interface mass flux to mass transfer rate per volume in the VOF model is improved.A new trans?formation method was proposed in this paper，which could solve the network independence.However，it also re?flects the mass transfer characteristics at a local interface.The transformation method was proved by theoretical in?ference.The corresponding condition that it is independent of mesh size was deduced and subsequently the simula?tion of bubble condensation was conducted by dividing three groups of meshes with different sizes.The transforma?tion method was simply applied in the aspect of bubble growth.The simulation results coincided with the theoretical analyses and experimental studies very well，proving that this method is feasible and can be widely applied in the simulation of the mass and heat transfer processes at the interface between the different phases.

VOF model；interface；local mass transfer；theoretical derivation；mesh independence；multiphase flow

10.3969／j.issn.1006?7043.201312053

http：／／www.cnki.net／kcms／detail／23.1390.U.20150109.1525.013.html

TL331

A

1006?7043（2015）03?0317?05

2013?12?17.網絡出版時間：2015?01?09.

核反應堆系統設計技術重點實驗室基金資助項目（KZA?KA?1101）；教育部留學歸國基金資助項目（2012?1707）；黑龍江省青年學術骨干支持計劃資助項目（1254G017）.

譚思超（1979?），男，教授，博士，博士生導師.

譚思超，E?mail：tansichao＠hrbeu.edu.cn.