遼寧科技大學材料與冶金學院  劉  鄧
中國科學院大學工程熱物理所  武廣龍

       【摘  要】利用有限容積法，針對純制冷劑R134a和非共沸混合制冷劑R407c，在水平管和U型管組合的物理模型內沸騰換熱過程進行數值模擬，得到了溫度、速度、干度及換熱系數的分布，并分析了管內工質的流速和熱流密度的變化對沸騰換熱系數的影響。此外，還對R134a在U型回轉式彎頭的水平光滑管內的沸騰換熱進行非穩態的數值模擬，可以清晰地觀察到氣泡的產生、躍離、生長和破碎的氣液兩相流型變化，并且模擬結果與Baker流型圖也較吻合。

       【關鍵詞】氣液兩相流 沸騰傳熱 數值模擬

Abstract:Adopting the method of Finite Volume Method to simulate the saturated boiling heat transfer in a smooth horizontal tube of R134a and R407c, and to get the distribution of temperature, velocity, dryness and heat transfer coefficient of R134a and R407c in the process of boiling heat transfer in tubes under the condition of different velocity and different heat flux density. In addition, the boiling heat transfer in tubes of R134a under the condition of unsteady three-dimensional numerical simulation is conducted separately, the changing process of the gas-liquid two phase flow pattern in tubes is observed, and the simulations give good agreement with the flow regimes expected from the Baker chart.

Key words:Gas liquid two-phase flow, Boiling heat transfer, Numerical simulation

1 引言

       在80年代，科學界確認氟氯烴（CFC）是引起臭氧破壞和溫室效應的危害物質之一。淘汰了R11、R12、R113、R114等普遍使用的CFC類制冷劑，即使R22等HCFC類對環境的污染破壞影響相對較小一些，最終也將被禁止使用。所以，研究并開發環保的新型制冷劑在制冷領域內是一個重要的研究課題。R134a是第一個被提出來的非臭氧破壞物質（ODP=0），具有較好的熱物理性質，其溫室效應潛能值（GWP）很低，不易燃，無毒，在冰箱、冷柜、汽車空調等制冷設備中已經開始使用。R407c具有清潔性、低毒性和不燃性等特點，尤其不破壞臭氧層（ODP=0），是目前使用較多的中高溫環保制冷劑。所以研究制冷劑R134a和R407c的換熱特性對今后的制冷領域有著重要作用。

       隨著在社會生產力提高的同時，也開始引起人們對研究氣液兩相流的注意。郭雷^[1]等著重研究在豎直矩形的細通道內水的沸騰換熱情況；歐陽新萍^[2]等在水平管外R404a和R407c的沸騰換熱進行了實驗研究；李靜^[3]模擬研究的是在豎直圓管內的沸騰換熱；吳曉敏^[4]等對CO₂/丙烷混合工質在水平管內的沸騰換熱進行分析。而本文選用的是將水平管和U型管兩者組合起來的物理模型，針對純制冷劑R134a和非共沸混合制冷劑R407c在帶有翅片的光滑管內沸騰換熱過程進行數值模擬，并進行對比分析。此外，還對在U型回轉式彎頭的水平光滑管內的沸騰換熱過程單獨做非穩態的數值模擬，以此觀察氣液兩相流的流型變化，并與Baker流型圖做了對比。

2 數值模擬的模型

       2.1 物理模型及參數

       圖2.1為模擬用的水平彎管的示意圖。其中管直徑D為25 mm，壁厚為0.5 mm，管長為12.89 m。選取R134a和R407c兩種工質。由于本文重點研究的是制冷劑在管內的流動換熱，所以創建幾何模型時，忽略管壁，只是針對管內區域劃分網格，然后進行了相關網格的獨立性檢驗，最后得到了符合條件的計算網格。本文采用固定熱流密度（30 kW/m²）的方法，對管壁使用壁面邊界條件，并且為0.5 m/s的入口流速，自由出流的出口流速。

圖2.1 水平彎管示意圖

       非穩態模擬采用的邊界條件與穩態模擬相同。此外，利用軟件Refprop8.0計算分別得到了R134a和R407c的熱物性參數，分別如表2.1和表2.2所示。

表2.1 R134a的熱物性參數

表2.2 R407c的熱物性參數

       2.2  數學模型

       控制方程可以表示成如下的通用形式^[5]：

     

       式中ф為通用變量；г_ф為廣義擴散系數；S_ф為廣義源項。

       制冷劑液體在管內蒸發時，流體被加熱達到所在壓力下的飽和溫度以后，會從液相變為氣相，所以兩相間必然會有質量傳遞。質量傳遞過程采用De Schepper等^[6]提出來的方程，不同相之間的質量源項如下：

       當T_liq≥T_sat時，制冷劑液體蒸發

       液相：  

       氣相：

       氣液兩相間的能量傳遞是伴隨著質量傳遞過程進行的，蒸發過程中能量源項的計算公式為：

       

       其中，T_sat為蒸發溫度；T_liq為制冷劑液相的溫度；α_liq為液相容積率；α_vap為氣相容積率；ρ_liq為制冷劑液相密度；為單位時間單位體積內向液相傳遞的質量，為單位時間單位體積內向氣相傳遞的質量；ΔH為汽化潛熱。以上模型公式針對的是純制冷劑工質，而對于混合制冷劑工質，擬采用Lee模型^[7]定義不同相之間的質量傳遞：

       液相：

       氣相：

       制冷劑在蒸發過程中，氣液兩相間的能量傳遞與純制冷劑工質的能量源項的計算公式類似：

       

       其中，T_b為泡點溫度；T_d為露點溫度；τ_l為液相的時間松弛系數，τ_v為氣相的時間松弛系數，建議取值為0.1^[7]。

3  數值模擬結果及分析

       3.1  R134a非穩態模擬結果分析

       管內沸騰換熱的兩相體積分數分布如圖3.1所示，其中液相的體積分數用紅色表示，氣相的體積分數用藍色表示。氣彈狀流、細泡狀流、波狀流和分層流的兩相流特定流型能在圖中很明顯的展現出來。分析：液體在進入剛管內時被加熱，所以換熱管開始部分并沒有氣泡，只能觀察到液相。首先被加熱到飽和溫度的液體開始蒸發，在這些部分首先開始形成氣泡，受重力影響會逐漸上升到管子頂部，隨著細泡狀流的形成，一系列其他兩相流流型也開始逐漸發展。液相的速度小于氣相的速度，當氣相的數量逐漸增加時，管內的兩相流動也變得愈發混亂。

圖3.1  管內沸騰換熱的兩相體積分數分布

       當液體繼續向管內流入時，氣相的體積分數和速度都會增加，也會觀察到氣彈狀流型。氣彈狀流區域的混亂程度比細泡狀流和氣塞狀流的混亂程度大，但其兩相交界面曲率與后兩者相比較會變緩和。第二管道后部，第三和四管道分別出現了出現氣彈狀、波狀流和分層流流型。許多采用水平管件的換熱器的設計必須用180°彎頭（回轉型彎頭）。由于回轉型彎頭的出現，氣相在管道彎頭處會重新分布，如圖3.2（a）所示。在一些絕熱彎頭處，還可能出現氣相的再冷凝現象。兩相的溫度場如圖3.2（b）所示，液體在管內加熱，當達到沸點283.15 K時，液體開始蒸發，氣相逐漸形成。就如圖3.2（b）所示，累積的氣相由于過熱其溫度有所上升，然而由于液體的替換流動和氣相在管頂部，所以只是壁溫有所上升。氣液兩相的速度分布具有明顯的差別如圖3.2（c）所示。而且，氣相速度往往大于液相速度，速度的差別導致了氣彈狀流和波狀流的形成和發展。

圖3.2  氣液兩相在管內的沸騰換熱過程

       在上述對非穩態模擬結果分析的基礎之上，本文利用相關的相應的流型圖與模擬得到的結果進行比對。圖3.3為在管道中心界面上選取的流型不同的6個點，與目前使用較多的Baker流型圖較為吻合，如圖3.4所示。

圖3.3  選取的6個不同流型點

圖3.4  Baker流型圖^[8]

       3.2  R134a穩態模擬結果分析

       由于換熱管較長，故選取具有特征的截面。如圖3.5所示可以看出，當換熱管管壁為恒定的熱流密度時，R134a從換熱管入口流動到出口進行沸騰換熱，管壁溫度高于工質的溫度，工質被加熱，其溫度逐漸升高，近壁面處的工質溫度高于工質中心部分的溫度，故溫度與溫度梯度沿管熱管管壁面的法線方向逐漸逐漸減小。

圖3.5  R134a在截面上的的溫度場分布

       如圖3.6所示，R134a從換熱管入口流動到出口進行沸騰換熱時，在入口速度一定的情況下，由于受管壁加熱的作用，工質由粘性系數相對較高的液相制冷劑逐漸沸騰轉化為粘性系數相對較低的氣相制冷劑，工質的沿程阻力也相應地逐漸減小，而且在同一壓力梯度下，蒸汽的密度較小，速度較大。在沿壁面法線方向，由于壁面處是首先發生沸騰換熱的地方，傳質阻力影響比較明顯，阻力沿法線方向起初較大，但逐漸減小。在邊界層內出現很明顯的速度梯度變化現象，即中心部分流速大，壁面附近流速相對較小，故出現沿管長方向R134a速度逐漸增大的現象。

圖3.6  R134a在截面上的速度場分布

       圖3.7為R134a在截面上的含氣率分布。由于管道的幾何的對稱性，氣體的體積分數的分布因此也具有很好的對稱性。R134a從換熱管入口流動到出口時，在加熱作用下，壁溫升高迅速，當獲得能夠沸騰的過熱度時，管內開始產生氣泡，同時氣泡的數量越來越多，故氣相體積分數也逐漸增大。在剛剛開始沸騰時，同一截面上，由于浮力和拽曳力等外力，壁面上產生的氣泡從壁面脫躍離，逐漸進入到主流中，故貼近管壁面處含氣率比管中心部分較高。隨著沸騰換熱過程的停止，大部分液相都轉化為氣相。此外，從圖中可以看出，在彎頭部位，其截面的體積含氣率往往比其前后的截面含氣率較低，出現氣相再冷凝的現象。

圖3.7  R134a在截面上的含氣率分布

       本文截取典型區域的圖像來展示管壁面上的換熱系數分布情況如圖3.8所示。由于工質在進口時已經處于飽和狀態，所以直接就開始飽和沸騰，沒有了過冷沸騰階段。飽和沸騰換熱時，由于管道上壁面的液相體積率遠小于氣相，隨著換熱系數的降低，換熱過程逐漸減弱。

圖3.8  R134a的局部換熱系數分布

       3.3  R407c模擬結果分析

       如圖3.9所示，當換熱管管壁為恒定的熱流密度時，R407c從換熱管入口流動到出口沸騰換熱時，管壁溫度高于工質的溫度，工質被加熱，其溫度逐漸升高，近壁面處的工質溫度高于工質中心部分的溫度，故換熱管中心區的溫度分布較為均勻，出現溫度與溫度梯度沿管熱管管壁面的法線方向逐漸減小的現象。

圖3.9  R407c在截面上的溫度分布

       如圖3.10所示，R407c從換熱管入口流動到出口進行沸騰換熱時，在入口速度一定的情況下，由于管壁的加熱作用，工質受熱開始發生沸騰，由液相向氣相開始轉變，氣相體積分數逐漸增大；在管長方向上，速度也隨之逐漸增大。在沿管徑方向，管壁受加熱作用溫度迅速上升，壁面處的液體工質首先出現氣泡，自然其處的傳質阻力的影響比較明顯，但是沿管徑方向逐漸減小，導致了中心部分流速大，壁面附近流速相對較小，故出現沿管長方向R407c速度逐漸增大的現象。

圖3.10  R407c在截面上的速度分布

       圖3.11為R407c的截面含氣率分布。由于管道的幾何的對稱性，除彎頭部分外，氣體的體積分數的分布也具有很好的對稱性。R407c從換熱管入口流動到出口時，在加熱作用下，壁溫升高迅速，當獲得能夠沸騰的過熱度時，管內開始產生氣泡，同時氣泡的數量越來越多，故氣相體積分數也逐漸增大。在沸騰開始階段，同一截面上，氣泡在壁面上產生，在各種力的作用下，從壁面躍離，變大，逐漸進入到主流中，故貼近管壁面處含氣率比管中心部分較高。隨著沸騰換熱過程的停止，大部分液相都轉化為氣相。此外，與圖3.7相比較，可以看出，由于流速增加，液體的替代作用增強，相似橫截面含氣率分布的截面出現的位置向后延遲出現。在彎頭部分，兩相會重新分布并且會有一定的再冷凝現象。

圖3.11  R407c在截面上的含氣率分布

       圖3.12是模擬得到的R407c的局部換熱系數沿管長方向上的分布和變化情況。由于換熱管細長，所以本文只是截取典型區域的圖像來展示管壁面上的換熱系數分布情況。由于工質在進口時已經處于飽和狀態，所以直接就開始飽和沸騰，沒有了過冷沸騰階段。飽和沸騰換熱時，由于管道上壁面的液相體積率遠小于氣相，隨著換熱系數的降低，換熱過程逐漸減弱。

圖3.12  R407c的局部換熱系數分布

       3.4  R134a和R407c沸騰換熱特性的對比分析

       3.4.1  R134a和R407c在管內沸騰換熱的溫度分布

       圖3.13中為制冷劑R134a和R407c溫度隨干度的變化情況。從（a）可以看出，在發生沸騰后流體的溫度并不是保持不變，而是有一定的升高，大約有2 K，這是由于此時傳熱的的熱量大部分用于補充制冷劑發生相變時的所需要的潛熱。雖然理論上講，流體發生相變時溫度保持不變，但是此結果與文獻[8]的結果有較好的吻合。從圖（b）可以看出，隨著干度的增加，流體溫度也逐漸增加，由于混合制冷劑R407c的非共沸的特性使其與純制冷劑相比具有非等溫特性，所以其溫度增加的幅度要遠大于同條件下的R134a的溫升。

圖3.13  制冷劑的溫度隨干度的變化情況

       3.4.2  R134a和R407c的沸騰換熱系數對比

       如圖3.14所示，R134a和R407c隨著干度的增加沸騰換熱系數先增加后減小。

圖3.14  制冷劑的沸騰換熱系數隨干度的變化情況

       對比分析圖3.15（a）、（b）得，R407c的平均沸騰換熱系數大于R134a的平均沸騰換熱系數。

圖3.15  制冷劑在不同流速下的平均沸騰換熱系數

       3.4.3  不同的流速與沸騰換熱系數的關系

       隨著制冷劑流速的增加，由于氣泡沿壁面滑移速度逐漸變快，使得氣泡從壁面躍離，故平均沸騰換熱系數隨著制冷劑流速的增加而增加，如圖3.16所示。

圖3.16  制冷劑在不同流速下的沸騰換熱系數

       3.4.4  不同的熱流密度與沸騰換熱系數的關系

       在泡狀沸騰時，熱流密度對換熱的影響要大一些，使得沸騰換熱系數增大，故出現沸騰換熱系數會隨著熱流密度的增加而增加的現象，如圖3.17所示。

圖3.17 制冷劑在不同熱流密度下的沸騰換熱系數

4  結論

       本文對R134a和R407c在水平光滑管內的沸騰換熱過程進行數值模擬。結果表明：

       （1）由于水平光滑管內除彎頭部分的橫截面的對稱性，使得R134a的速度、干度及溫度分布也具有對稱性。沿流動方向，工質的溫度逐漸升高，從非穩態模擬結果中可以看出氣泡數量逐漸增加，并且不斷躍離壁面、生長、破碎進入主流區。壁面邊界層內，由于粘性，速度由0沿管徑方向逐漸增加到主流速度，而速度梯度則與之相反。在管彎頭部分，由非穩態模擬結果可以看出，氣相和液相會重新分布，而且還可能會伴有氣相的再冷凝現象。

       （2）另外，混合制冷劑R407c的非共沸的特性使其與純制冷劑相比具有非等溫特性。

       （3）當制冷劑流速或熱流密度增大時，換熱管內的沸騰換熱系數也會增加。

       （4）將R134a在非穩態條件下的模擬結果與Baker流型圖進行了驗證，兩者結果比較吻合，模擬計算滿足了一定的合理性，為進一步模型優化提供了理論依據。

參考文獻：

       [1] 郭雷,張樹生,程林. R134a豎直矩形細通道內誰沸騰換熱的數值模擬[J].熱能動力工程, 2011,(1):31-35.
       [2] 歐陽新萍,劉超,林夢. R404A和R407C在水平強化管外的凝結換熱實驗研究[J].制冷學報, 2015,(4):72-77.
       [3] 李靜.三維豎直圓管內流動沸騰的數值模擬[J].廣東化工, 2011, 38(3): 75-76.
       [4] 吳曉敏,魏兆福,莫少嘉,等.CO2/丙烷混合工質管內流動沸騰換熱的數值模擬[J].工程熱物理學報, 2011, (6):1019-1022.
       [5] 陶文栓.數值傳熱學[M].第2版.西安:西安交通大學出版社,2001:5.
       [6] De Schepper S C K, Heynderickx G J, Marin G B. Modeling the evaporation of a hydrocarbon feedstock in the convection section of a steam cracker[J]. Computers & Chemical Engineering, 2009, 33(1): 122~132.
       [7] Lee W H. Pressure iteration scheme for two-phase flow modeling[J]. Multiphase Transport: Fundamentals, Reactor Safety, Applications, 1980: 407-432.
       [8] De Schepper S C K, Heynderickx G J, Marin G B. CFD modeling of all gas–liquid and vapor–liquid flow regimes predicted by the Baker chart[J]. Chemical Engineering Journal, 2008, 138(1): 349-357.

       備注：本文收錄于第21屆暖通空調制冷學術年會論文集。版權歸論文作者所有，任何形式轉載請聯系作者。