楊新飛，郭延柱，任麗，王海蕊，于朋玲

(1.山東省環(huán)境保護科學(xué)研究設(shè)計院，山東 濟南 250013；2.山東省化工研究院，山東 濟南 250014)

?

降膜蒸發(fā)器傳熱傳質(zhì)與流動過程數(shù)值模擬

楊新飛1，郭延柱1，任麗1，王海蕊1，于朋玲2

(1.山東省環(huán)境保護科學(xué)研究設(shè)計院，山東 濟南 250013；2.山東省化工研究院，山東 濟南 250014)

隨著計算機技術(shù)的發(fā)展，計算流體力學(xué)成為深入研究流動傳熱傳質(zhì)過程的重要方法。文章基于流體動力學(xué)基本原理和VOF多相流模型，建立了水平管外降膜流動與相變傳熱傳質(zhì)過程數(shù)學(xué)模型，針對降膜蒸發(fā)復(fù)雜的相變過程，通過用戶自定義函數(shù)UDF，將編制的計算程序嵌入FLUENT軟件相應(yīng)模塊，對降膜蒸發(fā)過程進行了模擬研究，闡明了管外液膜、速度、溫度分布及其局部傳熱傳質(zhì)特性沿圓周方向的變化規(guī)律。結(jié)果表明：在X=0.5～0.8時，液膜厚度較薄，有利于過程傳熱，圓管底部與頂部局部傳熱系數(shù)較大，并在管頂部出現(xiàn)傳熱系數(shù)最大值；相界面過余溫度沿周向距離的增大而增加，圓管下方具有最大的相界面換熱溫差和最高的蒸發(fā)速率；此數(shù)值模擬方法和傳熱傳質(zhì)模型可以較為合理地進行水平管外降膜蒸發(fā)過程的理論分析。

降膜; 蒸發(fā); 數(shù)值模擬; 傳熱; 傳質(zhì)

0　引言

在海水淡化過程中，流動介質(zhì)的工作溫度通常在70 ℃以下，開發(fā)小溫差下的高效蒸發(fā)器逐漸成為該領(lǐng)域研究的主流方向[1]。近年來，由于低溫多效蒸發(fā)系統(tǒng)中的水平管降膜蒸發(fā)器在小溫差低流量下具有高效的傳熱傳質(zhì)性能而備受關(guān)注[2]，許多學(xué)者對降膜蒸發(fā)器內(nèi)傳熱和流動過程進行了理論分析和實驗研究[3-8]。

數(shù)值模擬作為傳熱、傳質(zhì)研究的重要方法，已全面應(yīng)用到相關(guān)領(lǐng)域，作為理論分析、實驗研究或工程設(shè)計應(yīng)用的的重要手段，解決了一些實驗室或?qū)嶒炑b備由于條件受限而難以做到的情景模擬問題。

Jafar等通過CFD軟件進行了水平管二維數(shù)值模擬并研究了管排數(shù)對介質(zhì)流動形態(tài)的影響[9]。隨后，針對管間距、管束布置方式，管子直徑的模擬研究逐步展開[10-12]，噴淋密度、入口介質(zhì)溫度等對傳熱的影響也得到了進一步研究[13-14]。為了獲得更好的管外成膜特性，諸如橢圓管、蛋形管等異形降膜蒸發(fā)管也得到眾多學(xué)者的關(guān)注[15-16]，并且，計算也不再僅僅考慮流動，包括傳熱過程對膜厚和傳熱性能的影響也得到了深入研究分析[2.17-18]。然而，液膜蒸發(fā)這一相變過程由于涉及傳質(zhì)和傳熱，增加了數(shù)值模擬的復(fù)雜程度，因此，傳統(tǒng)的降膜蒸發(fā)管流動仿真通常不考慮相變蒸發(fā)傳熱過程，從而造成一定程度的計算誤差。

文章基于前人的研究，采用ANSYS Fluent6.3.26軟件，通過其二次開發(fā)接口編制了用戶自定義程序udf，將蒸發(fā)相變及傳熱過程考慮進來，分析了降膜蒸發(fā)器成膜過程、速度場、溫度場及不同噴淋密度下液膜厚度和傳熱系數(shù)的分布規(guī)律。模擬結(jié)果同已有文獻模擬及實驗結(jié)果進行了分析對比。

1　模型參數(shù)與計算方法

1.1物理模型網(wǎng)格及邊界條件

文章選取工業(yè)常用管徑為19 mm的水平圓管作為研究對象。物理模型及其坐標表示如圖1所示。選取流動介質(zhì)為60 ℃的飽和水，外部空間充滿此溫度下的飽和水蒸氣，介質(zhì)物性參數(shù)見表1。

水平管外降膜流動Reynolds數(shù)由式(1)表示為

Re=4Γ/μ

(1)

式中：Re為Reynolds數(shù)；Γ為單位管長上單側(cè)質(zhì)量流量，kg/(m·s)；μ為流體黏性系數(shù)，pa·s。

文章設(shè)定噴淋密度為0.29 kg/(m·s)，此時Re為2469，研究表明，該數(shù)下流動處于層流狀態(tài)[13]。模擬用二維模型，網(wǎng)格及具體邊界條件設(shè)置如圖1所示，網(wǎng)格采用四面體劃分方式，壁面進行邊界層處理。網(wǎng)格總數(shù)量為18973，網(wǎng)格無關(guān)性驗證結(jié)果表明[15]，此時已獲得網(wǎng)格無關(guān)性解。邊界設(shè)定速度入口，入口寬度為1 mm，飽和水布液高度為9 mm，管壁為無滑移壁面，計算采用標準壁面函數(shù)進行修正，恒壁溫Tw為333 K；設(shè)定管壁和飽和水接觸角為20 °。

圖1　圓管幾何模型、網(wǎng)格及邊界條件圖

物質(zhì)密度ρ/(kg·m-3)黏度μ/(kg·m-3·s-1)表面張力σ/(N·m-1)飽和水983.14.699×10-40.06622飽和蒸汽0.13021.0424×10-5-

1.2數(shù)學(xué)模擬和計算方法

模擬過程為常物性二維不可壓縮非穩(wěn)態(tài)層流兩相流動，過程具有明顯氣液相界面，而Volume of Fluid模型(VOF)相對于其他多相流模型更能準確的捕捉兩相界面。其引入了流體體積分數(shù)α的概念，各項流體體積分數(shù)之和等于1，由式(2)表示為

αl+αv=1

(2)

式中:αl、αv分別為飽和水與蒸汽的體積分數(shù)。

兩相體積分數(shù)控制方程可由式(3)表示為

(3)

由Lee提出的相變傳質(zhì)模型[19]可知,蒸發(fā)過程由式(4)表示為

(4)

動量方程由式(5)表示為

(5)

其中，ρ、μ分別由式(6)、(7)表示為

ρ=ρlαl+ρvαv

(6)

μ=μlαl+μvαv

(7)

式中:ρv為蒸汽的密度，kg/m3;μl、μv分別為蒸汽和飽和水的黏度，kg/(m3·s)。

(8)

式中：σ為表面張力,N/m;界面曲率κl由式(9)表示為

(9)

(10)

式中：h為單位質(zhì)量混合相焓值， J/kg；λ為介質(zhì)導(dǎo)熱系數(shù)， W/(m·K)；Sh為單位控制體積能量源項，J/(m3·s)。

混合相焓值h由式(11)表示為

(11)

其中，飽和水和蒸汽的焓值hl和hv分別由式(12)、(13)表示為

(12)

式中：Cp,l、Cp,v分別為飽和水與蒸汽的定壓熱容，J/(kg·k)。

介質(zhì)導(dǎo)數(shù)系數(shù)λ由式(13)表示為

λ=λlαl+λvαv

(13)

單位控制體積能量源Sh由式(14)表示為

(14)

式中：hfg為汽化潛熱， J/kg。

為方便結(jié)果的對比分析，相對坐標X及無量綱過余溫度Θ分別由式(15)、(16)表示為

(15)

(16)

式中：θ為柱坐標弧度值；T、Tw及T0分布為介質(zhì)溫度、管壁溫度及入流介質(zhì)溫度， K。

為衡量管外液膜傳熱特性，采用局部傳熱系數(shù)hw由式(17)計算

(17)

式中：qw為局部傳熱量，J/m2。

設(shè)定首相為飽和水，次相為飽和水蒸氣，氣液界面采用Geo-Reconstruct重構(gòu)方案。操作壓力為0.02 MPa，操作密度為0.1302 kg/m3?？刂品匠滩捎糜邢摅w積差分法進行離散，動量及能量方程的離散為二階差分格式，壓力離散為Body Force Weighted格式，速度壓力耦合方式為PISO。

2　結(jié)果及分析

2.1液膜分布

在層流降膜流動過程中，液膜厚度對傳熱具有顯著影響。圖2為液膜形成過程及鋪展情況，布液高度為9mm時，液膜從觸及圓管頂部到覆蓋整個圓管表面耗時約為0.06s。圖3為不同噴淋密度對液膜厚度的影響，由圖知，沿流動方向液膜厚度有逐漸減小的趨勢，圓管上側(cè)液膜厚度較下側(cè)大；當(dāng)X在0.5至0.8之間時，液膜厚度較薄；當(dāng)X大于0.8時，即液膜在圓管末端由于重力與表面張力的共同作用，液膜厚度增加。隨著噴淋密度的增大，管表面液膜厚度增加，噴淋密度過大，液膜厚度太大，不利于傳熱的進行；噴淋密度同樣不易選擇太小，較小的噴淋密度易引起下排管束的局部燒干[11]。

圖2　液膜形成過程及鋪展情況圖

圖3　不同噴淋密度對液膜厚度的影響

2.2速度分布

膜內(nèi)流動介質(zhì)速度分布對液膜導(dǎo)熱具有重要影響，較高的膜內(nèi)速度有利于對流換熱的進行。圖4為圓管外側(cè)速度分布，取值圖右側(cè)0.5m/s矢標。壁面處速度值為0，沿法向迅速增大，具有較大的速度梯度，液膜表面速度值最大，由于重力作用，液膜表面速度沿軸向逐漸增加，最大值出現(xiàn)在管子底部，約為0.75m/s，該位置處液體湍流擾動強烈，傳熱較強。受表面切應(yīng)力的影響，界面處飽和蒸汽速度值也較高，然后，沿法向減小至0。

2.3溫度分布

液膜表面溫度能夠體現(xiàn)液膜內(nèi)溫度的滲透性，越薄的液膜，溫度升高越快，越有利于界面處飽和水的蒸發(fā)。圖5為不同噴淋密度下周向相界面過余溫度Θi的分布。由圖知，噴淋密度越小，同位置處液膜表面溫度越高，由式(4)，蒸發(fā)速率越大。上側(cè)管外液膜蓄積能量，使得溫度在圓管底部得到顯著提升。液膜內(nèi)溫度的滲透點(即界面液體受壁面加熱溫度開始升高的點)隨噴淋密度的增加逐漸向后推移，噴淋密度為0.49kg/(m·s)時，發(fā)生在X=0.8位置處;噴淋密度為0.39kg/(m·s)時，發(fā)生在X=0.6位置處;噴淋密度為0.29kg/(m·s)時，發(fā)生在X=0.25左右。

圖4　圓管外側(cè)速度分布圖

圖5　不同噴淋密度下周向相界面過于溫度分布圖

2.4蒸發(fā)速率

管外降膜蒸發(fā)速率直接關(guān)系到蒸發(fā)量和降膜器效率，是降膜蒸發(fā)過程研究最為關(guān)心的問題之一。圖6為管外降膜蒸發(fā)速率。當(dāng)X<0.8時，管外的降膜蒸發(fā)速率較小;當(dāng)X>0.8時，蒸發(fā)速率迅速升高，最大蒸發(fā)速率約為0.3kg/(m3·s)。從圖5知，Γ=0.29kg/(m·s)，當(dāng)X>0.8時，氣液截面過余溫度值迅速升高，由式(4)知，較大的傳熱溫差可以提高傳質(zhì)速率，因此，降膜管底端蒸發(fā)速率始終要明顯高于上側(cè)氣液截面。

圖6　管外降膜蒸發(fā)速率圖

2.5傳熱特性

Γ=0.29kg/(m·s)時，管外局部傳熱系數(shù)沿周向變化曲線如圖7所示，由圖7可知，本次模擬得到平均傳熱系數(shù)為6.58kW/(m2·s)；而圓管底部與頂部局部傳熱系數(shù)較大，傳熱性能較好，并在管頂部出現(xiàn)傳熱系數(shù)最大值為11.51kW/(m2·s)，這是由于，在管子頂端X介于0～0.3，傳熱處于非正規(guī)狀況階段，局部傳熱量大，由式(17)知，局部傳熱系數(shù)較大；傳熱系數(shù)沿周向逐漸降低，當(dāng)X>0.9時，液膜流至圓管低端并開始脫落，由于兩側(cè)液膜對沖增大了膜內(nèi)液體擾動，強化了管壁側(cè)的傳熱特性，傳熱量增加，局部傳熱系數(shù)增大。

圖7　局部傳熱系數(shù)沿周向變化曲線圖

為了驗證本研究所用流動與傳熱傳質(zhì)模型，將現(xiàn)有的相近工況下研究結(jié)果同繪于圖7中進行比較發(fā)現(xiàn)，壁面?zhèn)鳠嵯禂?shù)沿周向的變化趨勢同已有研究結(jié)果基本一致，說明本模型及數(shù)值方法的合理性和可行性，文章模擬結(jié)果在X<0.5時同Ouldhadda等[13]模擬結(jié)果和Liu等[22]實驗結(jié)果相符，但是，當(dāng)X>0.5時，模擬結(jié)果低于Liu等[22]的實驗值。這可能與控制相變強度因子r的取值有關(guān)，它可以通過相應(yīng)的實驗結(jié)果不斷修正，進而獲取符合降膜流動過程的控制相變強度因子r。該模型可以為后續(xù)針對降膜蒸發(fā)過程及蒸發(fā)量的理論研究提供了一定程度的參考。

3　結(jié)論

通過上述研究可知：

(1) 管外液膜厚度沿流動方向逐漸減小，當(dāng)X大于0.8時，液膜厚度受重力與表面張力的共同作用而增加；

(2) 周向距離X大于0.8時，氣液界面過余溫度顯著升高，相應(yīng)的蒸發(fā)速率明顯高于上側(cè)氣液截面；

(3) 圓管局部傳熱系數(shù)沿周向逐漸降低，當(dāng)X>0.9時，由于擾動現(xiàn)象，強化了管壁側(cè)的傳熱特性，局部傳熱系數(shù)增大，尾部擾動有利于降膜蒸發(fā)過程的進行；

(4) 文章采用的數(shù)值模擬方法及傳熱傳質(zhì)模型可行并且較為合理，可以為后續(xù)針對降膜蒸發(fā)過程及蒸發(fā)量的理論研究提供了一定程度的參考，但如何選取合適的控制相變強度因子r尚需進一步的試驗及模擬修正。

[1]KronenbergG.,LokiecF..Low-temperaturedistillationprocessesinsingle-anddual-purposeplants[J].Desalination,2001,136(1):189-197.

[2]LuoL.C.,ZhangG.M.,PanJ.H., et al.Flowandheattransfercharacteristicsoffallingwaterfilmonhorizontalcircularandnon-circularcylinders[J].JournalofHydrodynamics,2013(3):404-414.

[3]范延品. 水平管降膜蒸發(fā)實驗研究[D].大連:大連理工大學(xué),2006.

[4]高國紅. 水平管降膜蒸發(fā)器傳熱性能研究[D].天津:天津大學(xué),2013.

[5]RifertV.G.,PodbereznyV.I.,PutilinJ.V., et al.Heattransferinthinfilm-typeevaporatorwithprofiletubes[J].Desalination,1989,74:363-372.

[6]HuX.,JacobiA.M..Theintertubefallingfilm:Part1—Flowcharacteristics,modetransitions,andhysteresis[J].JournalofHeatTransfer,1996,118(3):616-625.

[7]HuX.,JacobiA.M..Theintertubefallingfilm:Part2—Modeeffectsonsensibleheattransfertoafallingliquidfilm[J].JournalofHeatTransfer,1996,118(3):626-633.

[8]RibatskiG.,JacobiA.M..Falling-filmevaporationonhorizontaltubes—acriticalreview[J].InternationalJournalofRefrigeration,2005,28(5):635-653.

[9]JafarF.,ThorpeG.,TuranO.F..Liquidfilmfallingonhorizontalcircularcylinders[C].16thAustralasianFluidMechanicsConference(AFMC),Sydney.Australia:TheUniversityofQueensland,2007.

[10]何茂剛,王小飛,張穎,等. 水平管降膜蒸發(fā)器管外液體流動研究及膜厚的模擬計算[J].熱科學(xué)與技術(shù),2007(4):319-325.

[11]王小飛,何茂剛,張穎. 水平管降膜蒸發(fā)器管外液體流動數(shù)值模擬[J].工程熱物理學(xué)報,2008(8):1347-1350.

[12]邱慶剛,陳金波. 管束排列方式及管間距對水平管外液體成膜情況的影響分析[J].熱科學(xué)與技術(shù),2011(2):117-122.

[13]OuldhaddaD.,IdrissiA.I.,AsbikM..Heattransferinnon-Newtonianfallingliquidfilmonahorizontalcircularcylinder[J].HeatandMassTransfer,2002,38(7):713-721.

[14]SunF.,XuS.,GaoY..Numericalsimulationofliquidfallingfilmonhorizontalcirculartubes[J].FrontiersofChemicalscienceandEngineering,2012,6(3):322-328.

[15]羅林聰,張冠敏,田茂誠,等. 異形管降膜傳熱與流動數(shù)值模擬研究[J].工程熱物理學(xué)報,2013(4):710-714.

[16]呂多. 水平橢圓管外液體降膜流動數(shù)值模擬[D].遼寧,大連:大連理工大學(xué),2013.

[17]羅林聰, 潘繼紅, 田茂誠, 等. 管形對水平管降膜圓周膜厚和Nusselt數(shù)的影響[J]. 化工學(xué)報, 2013(8): 2760-2768.

[18]羅林聰. 水平異形管降膜蒸發(fā)流動與傳熱強化機理及實驗研究[D].濟南:山東大學(xué),2014.

[19]LeeW.H. .Apressureiterationschemefortwo-phaseflowmodeling[M].Washington:HemispherePublishing,1980.

[20]SandraC.K,,SchepperD.,GeraldineJ.H., et al.Modelingtheevaporationofahydrocarbonfeedstockintheconvectionsectionofasteamcracker[J].Computers&ChemicalEngineering,2009,31(1):122-132.

[21]孫東亮,徐進良,王麗. 求解兩相蒸發(fā)和冷凝問題的氣液相變模型[J].西安交通大學(xué)學(xué)報,2012(7):7-11.

[22]LiuP.J..TheEvaporatingFallingFilmonHorizontalTubes[D].Masison:UniversityofWisconsin-Madison,1975.

(學(xué)科責(zé)編：吳芹)

Numerical study on flow, heat and mass transfer characteristics of horizontal-tube falling-film evaporators

Yang Xinfei, Guo yanzhu, Ren li, et al.

(Shandong Academy of Environmental Science, Jinan 250013, China)

A mathematical model with phase transition was set up based on the principles of fluid dynamics and VOF multiphase flow model. The falling film evaporation procedure code was embedded in FLUENT through the user define functions (UDF). The principles in the circumferential direction of wall film thickness, velocity, temperature, mass transfer rate and local heat transfer coefficient was analyzed. Results show that the thickness of falling film was rather thin when the X was about 0.5～0.8, and this is good for the process of the heat transfer. The bottom and the top of the tube have the large value of local heat transfer coefficient, and the biggest value occurs at the top of the tube. In the circumferential direction, the surplus temperature increases with the distance grows, causing the biggest surplus temperature and evaporation rate at the bottom of the tube. The numerical result is in good agreement with the experimental and predicted data in literature, which means that the mathematical model with heat and mass transfer is reasonable and more comprehensive, and can be used in studying the flow, heat and mass transfer characteristics of falling water film outside horizontal tubes.

falling film; evaporation; numerical simulation; heat transfer; mass transfer

2015-07-05

楊新飛(1980-)，男，工程師，碩士，主要從事大氣污染防治等方面的研究.E-mail：52789044@qq.com

1673-7644(2016)01-0014-05

X512

A