邱慶剛，呂 多

（大連理工大學(xué) 能源與動力學(xué)院，遼寧 大連 116024）

0 引 言

水平管外液體的降膜流動具有傳熱溫差小、傳熱效率高、低溫傳熱性能優(yōu)良等優(yōu)點，因此被廣泛應(yīng)用于海水淡化、制冷系統(tǒng)、石油化工等行業(yè)中，在節(jié)能、環(huán)保方面有重大的意義［1－2］.國內(nèi)外許多學(xué)者對降膜流動過程進(jìn)行了卓有成效的研究.Ribatski等［3］的研究表明在水平管外液體降膜流動充分發(fā)展的層流區(qū)域，液膜與管壁的換熱是以導(dǎo)熱過程為主的對流換熱，液膜厚度與表面?zhèn)鳠嵯禂?shù)成反比關(guān)系，因此，研究液膜厚度對于了解液膜流動機(jī)理以及提高換熱系數(shù)有重大意義.Nusselt［4］通過理論分析給出了液膜厚度的經(jīng)典表達(dá)式ρl（ρl－ρg）gsinβ，Rogers［5］應(yīng)用動量方程積分方法得出了液膜厚度與噴淋密度、黏度、管周向角度關(guān)系的表達(dá)式；Gstoehl等［6］應(yīng)用光學(xué)測量法測得了沿管壁周向80°范圍內(nèi)的液膜厚度分布，并發(fā)現(xiàn)水平管外周向下半?yún)^(qū)域同理論預(yù)測值不符，液膜厚度偏?。辉S莉等［7］測量了無換熱情況下管外液體成膜的平均厚度及膜的幾率分布與流率、管徑的關(guān)系；何茂剛等［8］通過實驗測量得出了流量及布液高度對液膜厚度的影響規(guī)律；Harikrishnan等［9］通過二維數(shù)值計算方法計算了R134a水平管外降膜吸收過程的傳熱傳質(zhì)特性；Li等［10］通過實驗方法測試了雷諾數(shù)在21.6～108.1的水平管束降膜傳熱特性；Hou等［11］利用實驗測量了水平管壁周向不同角度的液膜厚度；Kostoglou等［12］利用統(tǒng)計學(xué)方法研究了湍流狀態(tài)下液體降膜流動的表面形態(tài)及液膜重組問題.

上述的研究工作雖取得了各種成果，但仍有許多工作要做，比如水平管外降膜流動過程中結(jié)構(gòu)及流動參數(shù)對液膜厚度與分布影響規(guī)律的相關(guān)文獻(xiàn)就較少.本文在前人研究基礎(chǔ)上，首先建立二維水平管外降膜流動模型，根據(jù)Rogers對水平管外降膜流動的分析，把液膜分為發(fā)展區(qū)與充分發(fā)展區(qū)，然后利用VOF 方法模擬不同雷諾數(shù)狀態(tài)下的水平管外降膜流動過程，通過精確讀取管壁周向180°的液膜厚度，分析管間距、雷諾數(shù)等結(jié)構(gòu)及流動參數(shù)對液膜厚度與分布的影響規(guī)律.

1 數(shù)值模擬方法

1.1 模型與基本參數(shù)

本文所研究的水平管外液體降膜流動模型參見圖1.圖2為水平管降膜蒸發(fā)器內(nèi)部管束排列示意圖，管束采用等邊三角形排列，管束上方設(shè)有布液器.本文主要研究雷諾數(shù)及管間距對液膜的影響，故選取垂直于管排方向一橫截面作為計算區(qū)域，為消除布液方式對液膜的影響，取第2根管讀取周向液膜厚度；為消除管排方式對液膜的影響，管間距s為變化值，布液器孔徑為3mm，布液高度為3mm，管徑R為19mm，數(shù)值計算區(qū)域為圖中陰影部分.

圖1 水平管外液體降膜流動示意圖Fig.1 Diagram of liquid falling film flow out a horizontal tube

圖2 管排方式及計算區(qū)域Fig.2 Arrangement of tube bundles and computational domain

1.2 流體物性與基本假設(shè)

液體降膜流動雷諾數(shù)定義為Re＝4Γ／μ，其中Γ為液體的噴淋密度，為管徑方向單位長度上流體的質(zhì)量流量，μ為液體的動力黏度.本文所計算的Re范圍均小于1 000，因此可認(rèn)為本文所涉及的液膜流動均在穩(wěn)定層流的范圍內(nèi).

假設(shè)蒸發(fā)器內(nèi)氣相全部為空氣，液相為飽和水，基本物性參數(shù)見表1.

表1 流體熱物性參數(shù)Tab.1 Thermo physical properties parameters of fluid

1.3 網(wǎng)格及邊界條件

本文采用精度較高的四邊形網(wǎng)格對計算區(qū)域進(jìn)行網(wǎng)格劃分，同時在氣液交界面及液相區(qū)進(jìn)行了網(wǎng)格加密.計算域的網(wǎng)格劃分如圖3所示，計算中比較了網(wǎng)格數(shù)分別為2 127、5 083、12 149的計算模型，結(jié)果差異較小，最終選取網(wǎng)格數(shù)為2 127的計算模型，在達(dá)到精度要求的前提下節(jié)省了計算時間；同時比較了時間步長分別為0.05、0.10、0.20ms的計算結(jié)果，經(jīng)比較0.10ms的時間步長能較為準(zhǔn)確且高效率地得出計算結(jié)果，殘差收斂于10－5.

圖3 網(wǎng)格劃分及局部網(wǎng)格放大圖Fig.3 Sketch of grid and local magnified grid

圖中速度入口為不同噴淋密度下的流體入口流速，壓力條件為1.013×105Pa，其余為對稱邊界.由于模擬條件無換熱過程，邊界條件中無溫度及熱流密度邊界.

1.4 控制方程及數(shù)值計算方法

微元控制體質(zhì)量守恒方程為

式中：ρ為流體密度；t為時間；x、y、z分別表示坐標(biāo)方向；u、v、w分別為速度矢量u在x、y、z方向的分量.由于ρ為常數(shù)，模型簡化為二維，故該式可轉(zhuǎn)化為＝0.

本文中流體為牛頓流體，N－S方程可簡化為

式中：p為壓強(qiáng)；μ為流體的動力黏度；ρgx、ρgy為微元體上的體積力，即液膜重力在各方向的分力.

本文采用VOF方法追蹤兩相流自由界面的流動情況，VOF方法的基本思想是引入某一相的體積分?jǐn)?shù)αi，對于某一計算單元，αi＝0說明第i相在該單元內(nèi)是空的；αi＝1說明第i相充滿該單元；0＜αi＜1說明該單元包含了多相流界面.αi的控制方程為

本文中設(shè)定空氣為第一相，水為第二相.計算單元中密度及黏度值為體積分?jǐn)?shù)加權(quán)平均值，由下式?jīng)Q定：

2 結(jié)果及分析

通過數(shù)值計算得出不同參數(shù)下的氣液兩相分布圖與速度云圖，并讀取管壁不同位置的液膜厚度值.圖4為管間距s＝6.4mm，Re＝574時的氣液兩相分布圖與速度云圖.從氣液兩相分布圖中可以看出此時液膜流動狀態(tài)較穩(wěn)定，液膜鋪展均勻連續(xù)；同時從速度云圖中也可看出管壁下方135°區(qū)域液膜流速相對較大.

圖4 s＝6.4mm，Re＝574時的氣液兩相分布圖與速度云圖Fig.4 Vapor－liquid phase and velocity distribution chart while s＝6.4mm，Re ＝574

圖5 液膜厚度隨角度變化圖Fig.5 Film thickness versus angle

本文在Re為574、744 時計算了s分別為6.4、9.5、19.4 mm 條件下管外液體降膜流動情況，讀取了管壁周向180°范圍內(nèi)的液膜厚度值（d），并同文獻(xiàn)［4］的理論值及文獻(xiàn)［6］的實驗值進(jìn)行對比，如圖5所示.從圖中可以看出，當(dāng)Re為574時，模擬值變化趨勢同實驗值比較吻合，液膜平均厚度隨管間距s的增加而減小，其原因可能是隨著布液高度的增大，液體流速變大，液膜形成初期滴狀流變得更加不穩(wěn)定，產(chǎn)生的液滴噴濺對膜流量存在較大影響，這也驗證了液膜厚度與液體流速成反比關(guān)系；同一條件下沿管壁周向液膜分布趨勢同實驗值較吻合，同文獻(xiàn)［4］的理論值在90°以上范圍內(nèi)存在一定偏差，液膜厚度總體上在沿管壁周向180°范圍內(nèi)先減小后增大，趨勢基本符合文獻(xiàn)［4］的拋物線理論值，但是在管壁下端，液膜厚度增長的趨勢較小，其原因可能是管壁上方90°區(qū)域?qū)儆谝耗恿鞒浞职l(fā)展區(qū)，液膜成膜較好，但波動較大；而在管壁下方90°區(qū)域內(nèi)，體積力占主導(dǎo)作用，液膜在黏滯力、表面張力與體積力共同作用下流動趨于穩(wěn)定，其中黏滯力使流體能均勻地附著在管壁上，表面張力維持液膜的形態(tài)，體積力保持液膜的流動，這可能是由文獻(xiàn)［4］分析中忽略液體動量的變化這一條件所造成的.同時，20°～45°液膜厚度波動相對較大，在45°～90°液膜厚度變化較為平穩(wěn)，這一點同文獻(xiàn)［5］把管外液體降膜流動分為發(fā)展區(qū)與充分發(fā)展區(qū)相符合，可近似地認(rèn)為在該種情況下45°以上區(qū)域為液膜流動的發(fā)展區(qū)，液膜厚度波動較大；之后區(qū)域為充分發(fā)展區(qū)，液膜厚度平穩(wěn)減小.

當(dāng)Re為744時，液膜厚度總體上增大，同時波動也增大，隨著雷諾數(shù)及管間距的增大，氣液兩相區(qū)邊界液膜厚度波動較大.在此雷諾數(shù)下，當(dāng)s＝19.4mm 時，液膜厚度沿管壁方向呈現(xiàn)較大的波動趨勢，模擬值普遍比實驗值以及理論推導(dǎo)值大，其原因可能主要有兩點：一是本文中采用的是二維簡化模型，計算時忽略了沿管長方向液膜間的相互影響；二是計算中邊界條件較為理想，并未估計擾動值，導(dǎo)致液體成膜情況較好，與實驗值存在一定差異.

Re為100～1 000，管壁周向角45°、90°、135°處液膜厚度隨Re變化如圖6所示.

圖6 液膜厚度隨雷諾數(shù)變化Fig.6 Film thickness versus Re

從圖中可以看出，管外壁周向某一角度處液膜厚度隨雷諾數(shù)的增大而增大，增長的速率近似為線性，管壁周向角為45°時，液膜厚度隨管間距的增大而減小，液膜增長曲線較平滑，其原因可能是液膜流動處于充分發(fā)展區(qū)，流動穩(wěn)定；周向角為90°時，液膜厚度增長曲線開始出現(xiàn)波動，管間距較大時波動較為明顯，其原因可能是此時液膜流速增大，流動不穩(wěn)定性增大；周向角為135°時，液膜厚度波動加劇，局部出現(xiàn)了“干區(qū)”，可見液膜厚度波動最為劇烈，此時，從速度云圖中可以看出液膜流速較大.重力與表面張力共同作用，導(dǎo)致液膜厚度波動較大，因此，水平管降膜流動液膜最薄處應(yīng)該出現(xiàn)在管壁周向90°～135°.

當(dāng)管間距大于12mm，液膜雷諾數(shù)大于200時，管壁下部會出現(xiàn)“干區(qū)”，如圖7所示，并且“干區(qū)”的范圍隨雷諾數(shù)及管間距的增大而增大.在水平管降膜蒸發(fā)器中，“干區(qū)”的出現(xiàn)會造成換熱系數(shù)降低，局部過熱引起設(shè)備的損壞，在工程中應(yīng)當(dāng)避免.可以通過在適當(dāng)范圍內(nèi)減小液膜雷諾數(shù)及管間距來避免“干區(qū)”出現(xiàn).

圖7 水平管壁“干區(qū)”示意圖Fig.7 The dryout region of horizontal tube

3 結(jié) 論

（1）在一定雷諾數(shù)下，水平管壁周向180°范圍內(nèi)液膜厚度的變化趨勢為先減小后增大，同實驗值及理論值吻合較好，證明了本文模型與算法的正確性和可靠性.

（2）隨著雷諾數(shù)的增大，液膜厚度波動增大，管壁周向某一角度液膜厚度逐漸增大，增大速率近似線性.

（3）管壁周向某一角度的液膜厚度隨著管間距的增大而減小.

（4）一定雷諾數(shù)與管間距條件下，管壁下部會出現(xiàn)“干區(qū)”，“干區(qū)”的范圍隨著雷諾數(shù)及管間距的增大而增大，工程實際中應(yīng)在適當(dāng)范圍內(nèi)盡量減小液膜雷諾數(shù)及管間距來避免“干區(qū)”的出現(xiàn).

［1］ 王小飛，何茂剛，張 穎.水平管降膜蒸發(fā)器管外液體流動數(shù)值模擬［J］.工程熱物理學(xué)報，2008，29（8）：1347－1350.WANG Xiao－fei，HE Mao－gang，ZHANG Ying.Numerical simulation of the liquid flowing outside the tube of the horizontal tube falling film evaporator ［J］.Journal of Engineering Thermophysics，2008，29 （8）：1347－1350.（in Chinese）

［2］ 楊 麗，王 文，白云飛，等.水平管降膜蒸發(fā)器傳熱優(yōu)化研究［J］.工程熱物理學(xué)報，2009，30（11）：1913－1916.YANG Li，WANG Wen，BAI Yun－fei，etal.Heat transfer optimization of the horizontal tube bundles in falling film evaporators ［J］.Journal of Engineering Thermophysics，2009，30（11）：1913－1916.（in Chinese）

［3］ Ribatski G，Jacobi A M.Falling－film evaporation on horizontal tubes — a critical review ［J］.International Journal of Refrigeration，2005，28（5）：635－653.

［4］ Nusselt W.Die oberflchenkondensation des wasserdampfes zeitschr ［J］.Verein Deutscher Ingenieure，1916，60（2）：541－546.

［5］ Rogers J T.Laminar falling film flow and heat transfer characteristics on horizontal tubes［J］.The Canadian Journal of Chemical Engineering，1981，59（2）：213－222.

［6］ Gstoehl D，Roques J F，Crisinel P，etal.Measurement of falling film thickness around a horizontal tube using a laser measurement technique［J］.Heat Transfer Engineering，2004，25（8）：28－34.

［7］ 許 莉 王世昌 王宇新 等.水平管外壁液膜流動狀態(tài)及其對傳熱的影響 ［J］.化工學(xué)報，2002，53（6）：555－559.XU Li，WANG Shi－chang，WANG Yu－xin，etal.Flowing state of liquid films over horizontal tubes and its influences on heat－transfer characteristics［J］.Journal of Chemical Industry and Engineering，2002，53（6）：555－559.（in Chinese）

［8］ 何茂剛，王小飛，張 穎.制冷用水平管降膜蒸發(fā)器的研究進(jìn)展及新技術(shù)［J］.化工學(xué)報，2008，59（2）：23－28.HE Mao－gang，WANG Xiao－fei，ZHANG Ying.Review of prior research and new technology for horizontal－tube falling－film evaporator used in refrigeration［J］.Journal of Chemical Industry and Engineering，2008，59（2）：23－28.（in Chinese）

［9］ Harikrishnan L，Shaligram T，Maiya M P.Numerical study of heat and mass transfer characteristics on a falling film horizontal tubular absorber for R－134a－DMAC ［J］.International Journal of Thermal Sciences，2011，50（2）：149－159.

［10］ LI Wei，WU Xiao－yu，ZHONG Luo，etal.Heat transfer characteristics of falling film evaporation on horizontal tube arrays［J］.International Journal of Heat and Mass Transfer，2011，54（9－10）：1986－1993.

［11］ HOU Hao，BI Qin－cheng，MA Hong，etal.Distribution characteristics of falling film thickness around a horizontal tube［J］.Desalination，2012，285：393－398.

［12］ Kostoglou M，Samaras K，Karapantsios T D.Reconstruction of film thickness time traces for wavy turbulent free falling films［J］.International Journal of Multiphase Flow，2010，36（3）：184－192.