水平管外降膜流動液膜厚度數值模擬

邱慶剛，呂 多

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

0 引 言

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

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

1 數值模擬方法

1.1 模型與基本參數

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

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

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

1.2 流體物性與基本假設

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

假設蒸發器內氣相全部為空氣，液相為飽和水，基本物性參數見表1.

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

1.3 網格及邊界條件

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

圖3 網格劃分及局部網格放大圖Fig.3 Sketch of grid and local magnified grid

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

1.4 控制方程及數值計算方法

微元控制體質量守恒方程為

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

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

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

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

本文中設定空氣為第一相，水為第二相.計算單元中密度及黏度值為體積分數加權平均值，由下式決定：

2 結果及分析

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

圖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°范圍內的液膜厚度值（d），并同文獻［4］的理論值及文獻［6］的實驗值進行對比，如圖5所示.從圖中可以看出，當Re為574時，模擬值變化趨勢同實驗值比較吻合，液膜平均厚度隨管間距s的增加而減小，其原因可能是隨著布液高度的增大，液體流速變大，液膜形成初期滴狀流變得更加不穩定，產生的液滴噴濺對膜流量存在較大影響，這也驗證了液膜厚度與液體流速成反比關系；同一條件下沿管壁周向液膜分布趨勢同實驗值較吻合，同文獻［4］的理論值在90°以上范圍內存在一定偏差，液膜厚度總體上在沿管壁周向180°范圍內先減小后增大，趨勢基本符合文獻［4］的拋物線理論值，但是在管壁下端，液膜厚度增長的趨勢較小，其原因可能是管壁上方90°區域屬于液膜層流充分發展區，液膜成膜較好，但波動較大；而在管壁下方90°區域內，體積力占主導作用，液膜在黏滯力、表面張力與體積力共同作用下流動趨于穩定，其中黏滯力使流體能均勻地附著在管壁上，表面張力維持液膜的形態，體積力保持液膜的流動，這可能是由文獻［4］分析中忽略液體動量的變化這一條件所造成的.同時，20°～45°液膜厚度波動相對較大，在45°～90°液膜厚度變化較為平穩，這一點同文獻［5］把管外液體降膜流動分為發展區與充分發展區相符合，可近似地認為在該種情況下45°以上區域為液膜流動的發展區，液膜厚度波動較大；之后區域為充分發展區，液膜厚度平穩減小.

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

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

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

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

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

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

3 結 論

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

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

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

（4）一定雷諾數與管間距條件下，管壁下部會出現“干區”，“干區”的范圍隨著雷諾數及管間距的增大而增大，工程實際中應在適當范圍內盡量減小液膜雷諾數及管間距來避免“干區”的出現.

［1］ 王小飛，何茂剛，張 穎.水平管降膜蒸發器管外液體流動數值模擬［J］.工程熱物理學報，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］ 楊 麗，王 文，白云飛，等.水平管降膜蒸發器傳熱優化研究［J］.工程熱物理學報，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］ 許 莉 王世昌 王宇新 等.水平管外壁液膜流動狀態及其對傳熱的影響 ［J］.化工學報，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］ 何茂剛，王小飛，張 穎.制冷用水平管降膜蒸發器的研究進展及新技術［J］.化工學報，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.