水平管降膜蒸發器殼程流場和溫度場數值研究

侯昊，畢勤成，張曉蘭

(西安交通大學 動力工程多相流國家重點實驗室，陜西 西安，710049)

淡水資源匱乏已成為當今世界威脅人類生存的嚴重問題。在地球的水資源總量中海水占97.5 %，從某種意義上說，海水是取之不盡、用之不竭的主要水源。因而，發展海水淡化技術，向海洋要淡水是解決淡水資源供需矛盾的有效途徑。海水淡化也就成為開發淡水資源的重要途徑之一。目前，已商業化的海水淡化方法有多種[1-2]，其中，低溫多效蒸餾海水淡化技術不僅熱能利用率高、產水率高，而且節能，近年來發展迅速，裝置規模日益擴大。另外，由于采用廉價材料，成本也日益降低，使其成為較為理想的海水淡化方式。因此，這種海水淡化技術備受各國學者的關注，成為當前海水淡化領域的研究熱點之一[3-5]。作為低溫多效蒸餾海水淡化系統的重要裝置，水平管降膜蒸發器受管內蒸汽和管外海水的熱工參數影響，殼程流場處于復雜的三維流動狀態，每一局部區域都經歷變質量的流動、熱力和傳熱過程，并且流動與換熱相互耦合，因此整個熱力過程錯綜復雜。另外，幾何結構復雜，又封閉在殼體內部，無法觀察和測量蒸發器內部所有位置的流場數據，一方面使研究者很難全面了解其內部工作過程，準確掌握整體以及局部區域真實流動和傳熱狀況的細觀信息；另一方面，又制約了試驗數據對于蒸發器結構改進和運行參數優化的指導作用。隨著多物理場耦合技術的發展，分布參數法因其能夠描述實際物理過程，獲得更加真實準確的參數分布，受到越來越多的重視。國內外一些學者已經將該方法應用于某些換熱器的熱力性能分析和優化設計[6-10]，而將其用于水平管降膜蒸發器的研究卻鮮有報道。為此，本文作者構建了水平管降膜蒸發器實體的三維分布參數模型，對蒸發器殼程區域的流場和溫度場特性進行了數值研究，分析和總結了該區域內流體流動和傳熱的細觀規律，直觀地刻畫了蒸發器內部的熱力過程，為深入了解復雜流動與傳熱現象提供了理論指導。

1 物理模型

圖1所示為水平管降膜蒸發器示意圖，其中管束分2個管程布置，1管程布置在下部，2管程在上部。換熱過程分為管內蒸汽冷凝和管外海水蒸發2部分：一方面，加熱蒸汽首先進入1管程并大部分冷凝，剩余蒸汽經管箱折返后進入2管程直至完全冷凝；另一方面，海水通過噴頭分布到管束上，自上而下降膜流動并逐漸吸收熱量，由過冷變為飽和狀態后開始蒸發。產生的二次蒸汽由管束周邊逸出，然后穿過絲網除沫器；在去除蒸汽中夾帶的海水液滴后，進入殼體兩側的蒸汽通道，并沿管束軸向流出。分管程布管在熱力學上是合理的，可以減小冷熱流體傳熱溫差，減少熵增引起的熱損失，另外，管內蒸汽壓力高于管外二次蒸汽壓力，即使傳熱管因腐蝕穿透也能保證產品水水質。

圖1 水平管降膜蒸發器示意圖Fig.1 Schematic of horizontal-tube falling-film evaporator

2 數學模型

2.1 基本假設

為了使模型推導方便，模型建立基于以下假設：(1) 忽略管排間液體的飛濺，即上排管底部落下的液量等于下排管頂部的初始液量；(2) 汽液界面處于熱力學平衡狀態，不考慮蒸汽的剪切力作用；(3) 不考慮傳熱管的污垢熱阻。

2.2 計算單元劃分

對于蒸發器內部不同區域的管束而言，管內蒸汽干度沿管子軸向逐漸減小，管外噴淋海水流量沿管排方向變化，從而造成蒸發器的換熱性能沿管子軸向和管排方向均有變化。管束示意圖如圖2所示。將蒸發器內管束從i，j，k3個方向分別劃分為Nx，Ny，Nz等份，其中i，j，k分別表示水平、垂直和軸向方向，這樣整個管束劃分為Nx×Ny×Nz個計算單元。這種方法的優點在于數學模型的尺度基于物理模型的實際拓撲結構，具有明確的物理意義，避免陷入CFD方法過于繁瑣的計算，極大地提高計算速度和加快程序收斂。針對每一單元建立傳熱與流動控制方程，使該單元出口參數成為下一單元的入口參數，采用迭代法求解出所有單元進出口的溫度、壓力等熱力參數。每一計算單元的熱負荷為：

式中：hoveral，hi和ho分別為總傳熱系數、管內和管外傳熱系數；A為傳熱面積；λ為導熱系數；Ti和To分別為管內、外工質熱力學溫度；di和do分別為管內、外直徑。

圖2 管束示意圖Fig.2 Schematic of tube bundle

為了計算每一單元的熱負荷和確定出口干度，需考慮蒸汽沿管子軸向和海水沿管排方向各熱工參數的變化，將物性、傳熱、流動等條件實現局部化，從而統一管內外側相變與非相變情況下的計算，分別獲得管內、外傳熱系數分布。另外，管內蒸汽的壓力、溫度變化與兩相壓降有關，因而還需計算管內汽液兩相壓降。

2.3 管內、外傳熱系數

2.3.1 管內傳熱系數

蒸汽進入管束時為飽和狀態，隨著熱交換的進行，經歷汽液兩相流和單相液體狀態。當蒸汽處于汽液兩相時，Chato[11]在Nusselt冷凝換熱理論基礎上進行解析，得到了水平圓管周向平均傳熱系數，Jaster等[12]又對其進行了修正，由此得到管內傳熱系數公式如下：

式中：ρL和ρG分別為工質液相和氣相的密度；r為汽化潛熱；μL為液相動力黏度；g為重力加速度；Tw為管子內壁面熱力學溫度；x為汽液兩相流干度。

當管內蒸汽完全冷凝后，管內傳熱系數采用齊德-泰勒公式[13]進行計算：

式中：l為管子長度；μw的定性溫度為管子內壁面熱力學溫度。

2.3.2 管外傳熱系數

海水進入蒸發器為過冷狀態，沿管排降膜并逐漸吸熱后達到飽和狀態。對于管外液膜處于過冷狀態，采用如下公式[14]進行計算：

對于管外液膜處于飽和狀態，采用如下公式[15]進行計算：

2.4 兩相壓降

水平管內汽液兩相流動壓降由重位壓降 Δpstatic、加速壓降Δpmom和摩擦壓降Δpfrict組成，即：

由于是水平管，重位壓降Δpstatic=0，加速壓降為：

式中：Gtotal為液汽兩相流總質量流速；ε為空泡份額；下標in和out表示計算單元的進口和出口。

摩擦壓降采用Friedel公式[16]進行計算：

式中：ΔpL為液相摩擦壓降；φ2為兩相摩擦倍率。

3 程序計算方法

根據上述計算模型，再加上對熱力過程其他一些現象的數學描述及求解，開發了水平管降膜蒸發器三維穩態過程的數值計算程序，能夠模擬不同幾何結構和流體物性參數條件下蒸發器內部流動與傳熱過程，具體計算步驟如下。

步驟 1 輸入已知條件：包括蒸發器幾何結構參數(傳熱管規格以及管排布置結構)、加熱蒸汽和供入海水各熱工參數和相關物性參數；

步驟2 對蒸發器管束劃分計算單元；

步驟3 預估管外飽和壓力、2管程蒸汽進口壓力和流量；

步驟4 依次求解2管程和1管程各排管的熱工參數，針對任一排管，先按照軸向順序分別計算各個單元，然后沿管排方向轉至下一排管進行計算，直到蒸發器的底排管；

步驟5 根據算得的管外飽和壓力、2管程蒸汽進口壓力和流量代替舊值，重新迭代計算直至收斂。

4 計算結果及分析

4.1 模型驗證

為了驗證模型及程序的正確性與可靠性，本文對某低溫多效海水淡化系統的水平管降膜蒸發器進行了模擬，計算中一共劃分單元26 970個，滿足網格獨立性要求，并將計算結果與該蒸發器的實際運行工況進行比較。圖3所示為計算值與實際值的對比。由圖3可見：在加熱蒸汽和海水參數相同的條件下，本文計算所得二次蒸汽參數與實際工況吻合良好，兩者之間的相對誤差在-4.2%～1.6 %范圍內，證明了本文所用模型具有很高的準確度。

4.2 海水流場及噴頭布置的影響

圖3 計算值與實際值的對比Fig.3 Comparison between calculated and actual results

圖4 海水流場分布Fig.4 Distribution of seawater flow field

圖5 改進噴頭布置后海水流場分布Fig.5 Distribution of seawater flow field for changing placement of spray nozzles

圖4和圖5所示為不同噴頭布置方式下蒸發器內海水流場云圖，直觀地給出了海水噴淋密度沿管排方向分布的完整數據。由圖4和5可知：蒸發器內噴頭布液效果良好，降膜傳熱過程順暢進行。由圖4可知：蒸發器采用均勻噴頭布置方式時，海水噴淋密度由上部的63 g/(ms)減小到下部的 40 g/(ms)，呈現出“上高下低”的分布規律，其原因在于過冷的海水需要2管程蒸汽加熱后才能達到飽和狀態，而且由于2管程蒸汽量較少，其換熱量也較少，即使海水已經飽和，產汽量也非常小，因此海水噴淋密度在2管程時大體上保持不變。到達1管程后，海水開始蒸發，產生二次蒸汽。在1管程中部海水噴淋密度由60 g/(ms)逐漸降低到46 g/(ms)，而且整體變化較為平穩，說明此區域換熱效果良好，產汽量保持穩定；而在1管程進口的下半區海水噴淋密度已經很小，其最小值為40 g/(ms)，但仍能夠保證傳熱管束周向完全潤濕，不會引起液膜破裂。同時，在該區域海水噴淋密度沿管長方向由40 g/(ms)逐漸增大到50 g/(ms)，表現出分布不均勻的特點，不利于避免“干斑”現象發生。因此，可以采用“前密后疏”的噴頭布置方式來改善這一區域的海水流動情況。如圖5所示，2管程上部區域海水分布不均，噴淋密度對應地由 70 g/(ms)逐漸減小到 60 g/(ms)，而1管程下部區域海水分布則更加均勻，噴淋密度大體上保持不變，約為48 g/(ms)。

4.3 海水溫度場

圖6所示為蒸發器內海水溫度場云圖，過冷的海水由50 ℃逐漸達到飽和溫度57 ℃，說明2管程已經達到了對海水的預熱作用。在2管程出口的上半區海水溫度始終保持在50 ℃，這是因為管內蒸汽已經完全冷凝，兩相換熱變為單相換熱，過冷海水無法得到足夠熱量使其溫度明顯上升，還要繼續吸熱才能達到飽和狀態。根據能量守恒原理，加熱蒸汽帶入的熱量等于進出口海水變化的熱量與二次蒸汽帶走的熱量之和，因此，可以通過對海水進行預熱提高海水的初始溫度來有效減小其過冷度，進而增大二次蒸汽所占的熱量份額，提高加熱蒸汽的熱利用率，從而增加單效蒸發器的二次蒸汽產量。

圖6 海水溫度場分布Fig.6 Distribution of seawater temperature field

4.4 傳熱特征

蒸發器內部的傳熱特征如圖7～9所示。由圖7可知：蒸汽進入管束后沿管長方向逐漸冷凝，通過釋放潛熱將熱量傳給管外海水，干度也逐漸降低且變化較為平穩，說明蒸汽沿管長流動過程中并未出現急劇冷凝現象，整體傳熱狀況良好。在1管程出口處平均干度約為0.15，剩余蒸汽經管箱折返后進入2管程直至完全冷凝。另外，總傳熱系數也可以表征各個計算單元的換熱效果，與管內、外傳熱系數密切相關，特別受蒸汽質量流速、干度以及海水噴淋密度影響較大。由圖8和圖9可知：在1管程和2管程蒸汽進口處總傳熱系數最大，約為3.3 kW/(m2K)，而后沿管長方向逐漸減小，說明進口區換熱較強，特別是2管程表現得更加明顯，因為該區域蒸汽干度大，管內擾動強烈，使得管內傳熱系數較大；而后沿管長方向冷凝換熱逐漸減弱，總傳熱系數也發生較大變化，管內單相流動時總傳熱系數很小，僅相當于管內蒸汽冷凝情況下總傳熱系數的7%左右。在1管程進口的下半區總傳熱系數略有減小，原因在于上半區海水已經達到飽和狀態，開始蒸發汽化，噴淋密度沿管排方向逐漸減小，管外液膜波動減緩，削弱了管外換熱，進而影響了總體傳熱。此外，整個蒸發器總傳熱系數的平均值為2.8 kW/(m2K)，已經達到了設計參數要求。

圖7 蒸汽干度分布Fig.7 Distributions of steam quality

圖8 2管程總傳熱系數分布Fig.8 Distributions of overall heat transfer coefficient for second tube pass

圖9 1管程總傳熱系數分布Fig.9 Distributions of overall heat transfer coefficient for first tube pass

4.5 海水鹽度分布

蒸發器內海水鹽度的變化情況如圖10所示。與圖4所表述的流場分布一一對應，表現出“上小下大”的趨勢。由于2管程海水的噴淋密度基本保持不變，其鹽度也維持在36 g/kg。當海水到達1管程后開始蒸發，產生二次蒸汽，噴淋密度逐漸減小，鹽度則隨之升高。在1管程中部海水鹽度整體變化較為平穩，從36 g/kg逐漸增加到48 g/kg。在1管程進口的下半區海水鹽度明顯升高，最大值達到56 g/kg，再次說明此區域換熱強烈，生成較多的二次蒸汽。但是海水鹽度過高，會促使其在蒸發過程中所含鹽分受熱分解并轉化成不同形態物質沉淀在傳熱管表面形成結垢，因此在工程設計中應當格外引起注意，控制海水濃縮比不超過2.5，最高溫度低于70 ℃可以避免發生上述現象。

圖10 海水鹽度分布Fig.10 Distributions of seawater salinity

5 結論

(1) 建立了水平管降膜蒸發器實體的三維分布參數模型，數值解與實際值的相對誤差在-4.2%～1.6%范圍內，驗證了該模型的可靠性，該模型適用于單元復疊循環計算且成本低廉。

(2) 通過模擬蒸發器內部的熱力過程，獲得了大量流場、溫度場和傳熱參數分布信息，真實地反映了蒸發器內部復雜的流動與傳熱特性，以細觀角度揭示出局部區域流體流動和傳熱的深層次問題。

(3) 三維分布參數模型是進行水平管降膜蒸發器數值研究的一種合理且實用的模型，對蒸發器的機理研究及整體的優化設計具有指導意義。

[1]周赤忠,李焱. 當前海水淡化主流技術的分析與比較[J]. 電站輔機,2008,29(4) : 1-5.ZHOU Chi-zhong,LI Yan. Analysis and comparison between current major technologies of seawater desalination[J]. Power Station Auxiliary Equipment,2008,29(4): 1-5.

[2]宿翠霞,倪華秋,李凌霄,等. 海水淡化技術及其應用現狀[J].中國資源綜合利用,2009,27(7) : 31-32.SU Cui-xia,NI Hua-qiu,LI Ling-xiao,et al. Seawater desalination technology and its application[J]. China Resources Comprehensive Utilization,2009,27(7): 31-32.

[3]Alasfour F N,Darwish M A,Bin Amer A O. Thermal analysis of ME-TVC+MEE desalination systems[J]. Desalination,2005,174(1): 39-61.

[4]Darwish M A,El-Dessouky H. The heat recovery thermal vapor-compression desalting system: A comparison with other thermal desalination processes[J]. Applied Thermal Engineering,1996,16(6): 523-537.

[5]Dey P K,Hawlader M N A,Chou S K,et al. Performance of a single-effect desalination system operating with different tube profiles and materials[J]. Desalination,2004,166(15): 69-78.

[6]Zhang L N,Yang C X,Zhou J H. A distributed parameter model and its application in optimizing the plate-fin heat exchanger based on the minimum entropy generation[J]. International Journal of Thermal Sciences,2010,49(8): 1427-1436.

[7]Rodriguez M,Diaz M S. Dynamic modeling and optimisation of cryogenic systems[J]. Applied Thermal Engineering,2007,27(7):1182-1190.

[8]Wang F Q,Maidment J F,Missenden J F,et al. A novel special distributed method for dynamic refrigeration system simulation[J]. International Journal of Refrigeration,2007,30(5):887-903.

[9]Sultana P,Wijeysundera N E,Ho J C,et al. Modeling of horizontal tube-bundle absorbers of absorption cooling system[J].International Journal of Refrigeration,2007,30(4): 709-723.

[10]Zavala-Rio A,Santiesteban-Cos R. Reliable compartmental models for double-pipe heat exchangers: An analytical study[J].Applied Mathematical Modelling,2007,31(9): 1739-1752.

[11]Chato J C. Laminar condensation inside horizontal and inclined tubes[J]. American Society of Heating Refrigerating and Airconditioning Engineers Journal,1962,4: 52-60.

[12]Jaster H,Kosky P G. Condensation heat transfer in a mixed flow regime[J]. International Journal of Heat and Mass Transfer,1976,19(1): 95-99.

[13]楊世銘,陶文銓. 傳熱學[M]. 北京: 高等教育出版社,1998:168-169.YANG Shi-ming,TAO Wen-quan. Heat transfer[M]. Beijing:Higher Education Press,1998: 168-169.

[14]Sernas V. Heat transfer correlation for subcooled water films on horizontal tubes[J]. Transactions of the ASME,1979,101(1):176-178.

[15]Parken W H,Fletcher L S,Sernas V,et al. Heat transfer through falling-film evaporation and boiling on horizontal tubes[J].Journal of Heat Transfer,1990,112(3): 744-750.

[16]Ould Didi M B,Kattan N,Thome J R. Prediction of two-phase pressure gradients of refrigerants in horizontal tubes[J].International Journal of Refrigeration,2002,25(7): 935-947.