微細(xì)通道內(nèi)R1234ze(E)和R134a冷凝特性的數(shù)值模擬

顧 昕，文 鍵，張 星，王斯民，厲彥忠

(1.西安交通大學(xué) 能源與動(dòng)力工程學(xué)院，陜西 西安 710049;2.西安交通大學(xué) 化學(xué)工程與技術(shù)學(xué)院，陜西 西安 710049)

1 前 言

由于科學(xué)技術(shù)的飛速發(fā)展，現(xiàn)代工業(yè)對(duì)換熱器緊湊性和高效性要求越來越高。采用微通道換熱器不僅滿足了緊湊性要求，還能利用管內(nèi)兩相流動(dòng)帶來的相變潛熱以獲得強(qiáng)大的換熱能力。此外，它還能減少系統(tǒng)的充灌量，滿足節(jié)能要求。探索微通道內(nèi)兩相流動(dòng)的規(guī)律對(duì)于換熱器的設(shè)計(jì)和優(yōu)化具有重要意義。目前已有大量的學(xué)者對(duì)微細(xì)通道內(nèi)冷凝流動(dòng)特性進(jìn)行了研究。CAVALLINI等[1]和AWAD等[2]對(duì)微細(xì)通道內(nèi)流型、換熱系數(shù)及壓降的實(shí)驗(yàn)研究進(jìn)行了文獻(xiàn)綜述。大部分研究基本都以氫氟烴作為研究工質(zhì)，如R134a和R410A。雖然這類工質(zhì)進(jìn)入大氣后不會(huì)對(duì)臭氧層造成破壞，但是它們的GWP (global warming potential)高，在大氣停留時(shí)間長(zhǎng)，長(zhǎng)期使用會(huì)加劇溫室效應(yīng)?？梢灶A(yù)見的是，類似于R134a和R410A這類高GWP的含氫氟烴類制冷劑終將會(huì)被淘汰。四氟丙烯(R1234ze(E))不僅ODP(ozone depletion potential)為0，GWP值也很低，其物理性能與R134a相似，被認(rèn)為是未來可替代R134a的新一代制冷劑。MOTA-BABILONI等[3]綜述了近年來關(guān)于 R1234ze(E)的研究進(jìn)展，發(fā)現(xiàn)目前的研究主要圍繞其物性、可燃性、與油的相容性以及蒸汽壓縮系統(tǒng)的性能展開。關(guān)于R1234ze(E)在水平管內(nèi)的冷凝兩相流動(dòng)的研究很少，且都是采用實(shí)驗(yàn)法進(jìn)行研究。WANG等[4]研究了R1234ze(E)在水力直徑為301.6 μm的微通道序列中的流型分布，觀察到了環(huán)狀流、注射流、波狀流和泡狀流。LIU等[5]研究了丙烷、R1234ze(E)和 R22在水平圓管和方管內(nèi)的冷凝換熱系數(shù)和壓降，發(fā)現(xiàn)丙烷和R1234ze(E)換熱系數(shù)和壓降都要大于R22。JIGE等[6]實(shí)驗(yàn)研究了R32、R134a、R410A 和 R1234ze(E)在水平多孔管內(nèi)的冷凝換熱系數(shù)及壓降，發(fā)現(xiàn) R1234ze(E)的壓降梯度要大于其他工質(zhì)。DEL COL等[7]對(duì)比了R1234ze(E)、R32、R134a 和R1234yf在水平圓管(Dh=1 mm)內(nèi)的冷凝換熱系數(shù)及絕熱壓降梯度。研究發(fā)現(xiàn)R1234ze(E)的壓降梯度要比R134a和R1234yf大。相同功耗下，R1234ze(E)和R1234yf所需要的換熱面積要比R134a分別大25%和15%。HOSSAIN等[8]對(duì)比了R1234ze(E)、R32和 R410A在水平圓管內(nèi)的冷凝換熱性能，發(fā)現(xiàn) R1234ze(E)的換熱系數(shù)要比 R32低20%～45%，但要比R410A高10%～30%。

目前關(guān)于R1234ze(E)在微細(xì)通道內(nèi)冷凝流動(dòng)的研究非常有限，而從制冷劑替代的發(fā)展趨勢(shì)來看，這方面的研究很有必要?，F(xiàn)有的研究主要以實(shí)驗(yàn)為主，缺乏相應(yīng)的理論分析。本文通過數(shù)值模擬手段，研究了R1234ze(E)在水平管內(nèi)的冷凝流動(dòng)，并與R134a進(jìn)行了對(duì)比，討論了液膜分布規(guī)律以及剪切力、重力和表面張力的相對(duì)作用，并就現(xiàn)有關(guān)聯(lián)式對(duì)R1234ze(E)換熱系數(shù)和壓降的適用性進(jìn)行了評(píng)估和討論。

2 數(shù)學(xué)模型

2.1 計(jì)算域

如圖1所示，本文采用的幾何模型是管徑為1 mm的水平微細(xì)圓管。管道入口為飽和蒸汽(入口質(zhì)量流速的變化范圍為400～800 kg?m-2?s-1)，隨著冷凝過程進(jìn)行，氣相減少而液相增多，管內(nèi)依次會(huì)出現(xiàn)霧狀流、環(huán)狀流以及間歇流。本文通過NEMA等[9]提出的流型轉(zhuǎn)化準(zhǔn)則來判斷流型變化，保證管道的出口干度處于環(huán)狀流范圍，從而確定管長(zhǎng)為130 mm。由于本文采用三維穩(wěn)態(tài)模擬，入口段的霧狀流假設(shè)為環(huán)狀流。計(jì)算域進(jìn)口為充分發(fā)展的速度入口邊界條件，出口為壓力出口。管壁設(shè)定為光滑無滑移且壁面溫度為常數(shù)(303 K)。工質(zhì)的物性參數(shù)來源于REFPROP 9.0。計(jì)算域全部采用六面體網(wǎng)格填充，網(wǎng)格從壁面到管道中心逐漸加密。近壁面網(wǎng)格第一層高度為0.8 μm且Y+小于1。對(duì)網(wǎng)格進(jìn)行無關(guān)性驗(yàn)證，當(dāng)網(wǎng)格數(shù)為798萬和470萬時(shí)，換熱系數(shù)的誤差在0.5%以內(nèi)，因而最終網(wǎng)格數(shù)確定在470萬左右。

圖1 幾何示意圖及邊界條件Fig.1 Geometric schematic diagram with boundary conditions

2.2 VOF模型

通過求解每一相的體積分?jǐn)?shù)φ，VOF模型可有效追蹤互不相混不可壓縮流體之間的相界面，現(xiàn)已被成功用于模擬管內(nèi)的冷凝流動(dòng)。對(duì)于VOF模型，每一個(gè)控制體內(nèi)各項(xiàng)的體積分?jǐn)?shù)之和必須為1。

VOF模型將各相混合物看成一種流體，而流動(dòng)變量和物性參數(shù)都取體積平均值

由于各相看成一種流體，因此不管是動(dòng)量方程還是能量方程都是單一求解。對(duì)于環(huán)狀流動(dòng)，可假設(shè)冷凝過程為穩(wěn)態(tài)。求解方程如下：

體積分?jǐn)?shù)方程：

連續(xù)性方程：

動(dòng)量方程：

為考慮表面張力的影響，本文采用BRACKBILL等[10]提出的CSF模型(continuum surface force)，使表面張力可以表達(dá)成體積力的形式從而以源項(xiàng)的方式添加到動(dòng)量方程，其表達(dá)式為：

能量方程：

式中

2.3 湍流模型

在本文的研究工況范圍內(nèi)，ReG的范圍為17 127～252 591，而ReL的范圍是210～5 293。從各相雷諾數(shù)來看，氣相可以確定為完全湍流，而液相在入口段由于干度很大，可能處于層流，而層流向湍流轉(zhuǎn)換的位置不能確定。本文假定氣液兩相一直處于湍流狀態(tài)。根據(jù) DA RIVA等[11]的研究，低雷諾數(shù)的 SST(shear stress transport)k-ω模型一方面可有效處理近壁區(qū)層流到湍流的過渡，另一方面，其引入了可避免渦黏性被過高預(yù)測(cè)的模型。

2.4 相變模型

本文采用LEE[12]模型模擬氣液兩相之間的傳質(zhì)過程。該模型假定氣液相界面為飽和溫度，當(dāng)控制體的溫度小于飽和溫度且氣相分?jǐn)?shù)大于0時(shí)，氣相向液相傳遞質(zhì)量。當(dāng)氣相體積分?jǐn)?shù)為0時(shí)，傳質(zhì)過程不發(fā)生。

本文通過 UDF(user defined function)的形式將上述模型以源項(xiàng)的方式添加到各相的體積分?jǐn)?shù)方程及能量方程中。關(guān)于系數(shù)r的確定，DA RIVA等[13]指出r太大會(huì)導(dǎo)致計(jì)算不收斂，太小會(huì)使相界面度嚴(yán)重偏離飽和溫度。本文通過試錯(cuò)法來調(diào)整r的取值，并通過將模擬的換熱系數(shù)與實(shí)驗(yàn)值對(duì)比來驗(yàn)證r的準(zhǔn)確性。根據(jù)工況和幾何結(jié)構(gòu)的不同，本文r取值為 3×105～1.5×106s-1，每種情況下相界面的溫度與飽和溫度的差值都小于1 K。

圖2 換熱系數(shù)模擬值與實(shí)驗(yàn)數(shù)據(jù)對(duì)比Fig.2 Comparison of simulated HTC with experimental data of DEL COL et al.[7]

本文采用壓力基穩(wěn)態(tài)求解，采用SIMPLE算法實(shí)施壓力與速度的耦合，VOF方程采用隱式格式離散，操作密度設(shè)為氣相密度。動(dòng)量、能量、湍動(dòng)能和湍動(dòng)黏度方程采用三階MUSCL格式離散。采用bounded gradient maximization (BGM)算法捕捉相界面。

3 結(jié)果與討論

3.1 模型驗(yàn)證

為驗(yàn)證數(shù)學(xué)模型的準(zhǔn)確性，本文將數(shù)值計(jì)算得到的結(jié)果與DEL COL等[7]的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行比較。圖2分別展示了相同工況下?lián)Q熱系數(shù)以及壓降的模擬值與實(shí)驗(yàn)測(cè)量值隨干度的變化關(guān)系。由圖可知，本文模擬得到的換熱系數(shù)和壓降梯度不管是趨勢(shì)上還是數(shù)值上都能和實(shí)驗(yàn)數(shù)據(jù)較好地符合，最大偏差分別控制在16.8% 和22.8% 以內(nèi)，說明了本文數(shù)學(xué)模型的準(zhǔn)確性。

3.2 換熱系數(shù)

圖3展現(xiàn)了不同質(zhì)量流速下R1234ze(E)和R134a的換熱系數(shù)隨干度的變化趨勢(shì)。對(duì)于R1234ze(E)和R134a，換熱系數(shù)都隨干度和質(zhì)量流速的增大而增大。當(dāng)干度小于0.75時(shí)，換熱系數(shù)與干度幾乎成線性關(guān)系。相同質(zhì)量流速和干度下，R134a的換熱系數(shù)要稍大于R1234ze(E)，特別是在高干度區(qū)和高質(zhì)量流速的時(shí)候。當(dāng)G= 800 kg?m-2?s-1時(shí)，R1234ze(E)的換熱系數(shù)平均要比R134a低8.87%。

圖3 換熱系數(shù)隨干度的變化Fig.3 Profiles of condensation heat transfer coefficients as a function of vapor quality

圖4 液膜厚度隨干度的變化Fig.4 Profiles of liquid film thickness as a function of vapor quality

3.3 液膜厚度

圖4展現(xiàn)了G= 400 kg?m-2?s-1時(shí)，R1234ze(E)和R134a局部液膜厚度隨干度的變化關(guān)系。由圖可知，液膜厚度隨干度的增大而減小。在高干度區(qū)，管內(nèi)冷凝的初始階段，氣相速度很大而液相速度很小。由于氣液兩相之間較大的相對(duì)速度差，產(chǎn)生了很大的相間剪切力，減薄了液膜厚度。相對(duì)重力和表面張力，剪切力在高干度區(qū)起主導(dǎo)作用。隨著干度的降低，由于氣液兩相之間的動(dòng)量交換，氣相速度逐漸降低而液相速度逐漸增大，相間剪切力減小，液膜變厚。此外，盡管R134a的換熱系數(shù)較高，相同工況下其液膜厚度也更厚。 R1234ze(E)的液膜厚度平均要比 R134a薄 15.7%。由于 R1234ze(E)的氣液密度比相對(duì)R134a要大，因此相同質(zhì)量流速和干度下，R1234ze(E)的氣液相速度比更大，相間剪切力增強(qiáng)從而減薄了液膜厚度。本文的模擬研究假設(shè)壁溫為常數(shù)，因而局部換熱系數(shù)的大小和局部的熱流量成正比。DA RIVA等[11]討論了湍流效應(yīng)對(duì)R134a在水平圓管內(nèi)換熱特性的影響，發(fā)現(xiàn)在高質(zhì)量流速下，由于液相內(nèi)部的湍流效應(yīng)，局部熱流量的大小不再只取決于液膜厚度，還取決于液相有效熱導(dǎo)率的大小。

3.4 有效熱導(dǎo)率

圖5展現(xiàn)了不同干度下有效熱導(dǎo)率的分布。圖中圓形符號(hào)和方形符號(hào)分別代表 R134a和R1234ze(E)氣液相界面的位置。對(duì)任意橫截面(XY平面)，Y= 0 m代表圓管的中心位置，Y= -0.005 m代表圓管的底部。由圖可知，相同情況下R134a的有效熱導(dǎo)率要大于 R1234ze(E)。當(dāng)G= 400 kg?m-2?s-1，x= 0.4和θ= 180°時(shí)，R134a液相的平均有效熱導(dǎo)率比 R1234ze(E)大 28%。這或許可以解釋即使R134a的液膜厚度更厚但是其換熱系數(shù)仍要大于R1234ze(E)。本文的結(jié)果說明，在衡量不同工質(zhì)在水平圓管內(nèi)的冷凝換熱特性時(shí)，液膜厚度的大小并不是決定性因素，有效熱導(dǎo)率的影響不可忽視。

圖5 有效熱導(dǎo)率沿Y軸的變化Fig.5 Profiles of effective thermal conductivity along Y-axis

3.5 液膜分布特征

圖6展現(xiàn)了冷凝過程中不同干度下R134a和R1234ze(E)相界面的分布。整體來看，R1234ze(E)的相界面分布特性與R134a相似。在高干度區(qū)時(shí)，剪切力占主導(dǎo)，重力的影響可以忽視，液膜沿圓周均勻分布。隨著干度的降低，在圓管頂部的冷凝液由于重力的作用積聚到底部，逐漸開始分層。

圖6 R1234ze(E)和R134a氣液相界面形狀Fig.6 Liquid–vapor interfaces of R1234ze(E) and R134a

圖7 局部液膜厚度沿周向位置的變化Fig.7 Profiles of local film thickness as a function of angular coordinate

圖7展現(xiàn)了R1234ze(E)和R134a局部液膜厚度隨周向位置的變化特性。由圖可知，當(dāng)干度比較高時(shí)(x=0.8)，R1234ze(E)和R134a不同周向位置處的液膜厚度幾乎一樣。且液膜厚度幾乎不隨周向位置發(fā)生變化。隨著干度的減小(x=0.5)，R1234ze(E)的液膜厚度明顯比 R134a要薄，主要體現(xiàn)在圓管的頂部和底部位置。此外，當(dāng)干度較小時(shí)，重力的影響增強(qiáng)，兩種工質(zhì)的液膜厚度都是隨著角度的增大而逐漸增大。

3.6 模擬結(jié)果與換熱關(guān)聯(lián)式對(duì)比

由于目前沒有可預(yù)測(cè)R1234ze(E)換熱性能的關(guān)聯(lián)式提出，本文嘗試將現(xiàn)有的關(guān)聯(lián)式去預(yù)測(cè)R1234ze(E)的換熱系數(shù)。本文通過定義平均絕對(duì)誤差(MAD)來衡量關(guān)聯(lián)式的預(yù)測(cè)效果，其計(jì)算方式如下：

如圖8所示，現(xiàn)有關(guān)聯(lián)式的預(yù)測(cè)誤差大概在±30%以內(nèi)。其中，SHAH 等[14]提出的關(guān)聯(lián)式預(yù)測(cè)效果最好，平均絕對(duì)誤差為10.14%。THOME等[15]和CAVALLINI等[16]提出的關(guān)聯(lián)式低估了本文的模擬值，MAD分別為17.27%和14.48%。而DOBSON和CHATO關(guān)聯(lián)式[17]一定程度上高估了模擬值，平均絕對(duì)誤差為11.48%。

3.7 壓降特性

對(duì)于管內(nèi)的兩相流動(dòng)，兩相流動(dòng)壓降主要由以下幾部分組成：摩擦壓降、加速壓降和重力壓降。其中，加速壓降是由于進(jìn)出口氣液含量變化而導(dǎo)致的壓降損失。由于本文研究的是水平管內(nèi)的冷凝流動(dòng)，因此重力壓降可以忽略。摩擦壓降等于總的壓降減去加速壓降，其表達(dá)式如下：

圖8 換熱系數(shù)模擬值與關(guān)聯(lián)式預(yù)測(cè)值對(duì)比Fig.8 Comparison of simulated heat transfer coefficients with predicted data

圖9展現(xiàn)了不同流量下R134a和R1234ze(E)摩擦壓降梯度隨干度的變化趨勢(shì)。由圖9可知，R134a和 R1234ze(E)摩擦壓降梯度隨干度和質(zhì)量流速的增大而增大。相同情況下，R1234ze(E)的摩擦壓降要大于 R134a。由于 R1234ze(E)的氣相密度要小于 R134a，相同質(zhì)量流速下，其氣相速度要大于 R134a。與此同時(shí)，R134a和R1234ze(E)的液相黏度幾乎一致。因此相同工況下，R1234ze(E)的流動(dòng)阻力大于R134a。

圖9 壓降梯度隨干度的變化Fig.9 Pressure gradients versus vapor quality

圖10 壓降梯度模擬值與關(guān)聯(lián)式預(yù)測(cè)值對(duì)比Fig.10 Comparison of simulated pressure gradients with predicted data

3.8 模擬結(jié)果與壓降關(guān)聯(lián)式對(duì)比

圖10展現(xiàn)了現(xiàn)有壓降關(guān)聯(lián)式對(duì)模擬得到的R1234ze(E)壓降梯度值得預(yù)測(cè)效果。很明顯，現(xiàn)有的關(guān)聯(lián)式對(duì) R1234ze(E)的壓降存在嚴(yán)重的低估，平均絕對(duì)誤差幾乎都在 30%左右。其中最好的是 FRIEDEL[18]提出的關(guān)聯(lián)式，MAD為27.06%。

4 結(jié) 論

采用VOF模型對(duì)R1234ze(E)在水平圓管內(nèi)的冷凝換熱及阻力特性進(jìn)行了數(shù)值模擬研究，討論了質(zhì)量流量、干度和物性對(duì)換熱系數(shù)和壓降梯度的影響。對(duì)于R1234ze(E)和R134a，換熱系數(shù)和壓降梯度都隨干度和質(zhì)量流速的增大而增大。相同情況下，R1234ze(E)換熱系數(shù)小于R134a，但壓降大于R134a。當(dāng)G= 800 kg?m-2?s-1時(shí)，R1234ze(E)的換熱系數(shù)平均要比R134a低8.87%。R1234ze(E)的液膜厚度平均要比R134a薄 15.7%。當(dāng)氣液兩相都為湍流，在衡量不同工質(zhì)在水平圓管內(nèi)的冷凝換熱特性時(shí)，有效熱導(dǎo)率的影響要大于液膜厚度。R1234ze(E)在管內(nèi)的液膜分布特性整體上和R134a相似。SHAH等提出的關(guān)聯(lián)式能準(zhǔn)確地預(yù)測(cè)R1234ze(E)的冷凝換熱系數(shù)，平均絕對(duì)誤差為10.14%。FRIEDEL，KIM and MUDAWAR，SUN and MISHIMA以及ZHANG and WEBB提出的壓降關(guān)聯(lián)式不能準(zhǔn)確地預(yù)測(cè)R1234ze(E)的壓降梯度，預(yù)測(cè)值偏小，誤差在30%左右。

符號(hào)說明：