繞管式換熱器殼側汽化模擬時的液滴粒徑研究

吳志勇 張安昊 國慧敏 建偉偉 蔡偉華

1.遼寧石油化工大學石油天然氣工程學院 2.哈爾濱工業(yè)大學能源科學與工程學院

繞管式換熱器(spiral wound heat exchanger，SWHE)作為一種結構復雜的換熱器廣泛用于天然氣液化領域，具有抗高壓、耐低溫、多種介質可同時換熱的特點。繞管式換熱器內部構造見圖1[1-2]，直徑較小的換熱管以螺旋方式分層纏繞在中央芯筒上，相鄰兩層換熱管的纏繞方向不同，層間用金屬隔條控制間距[3]。

繞管式換熱器用于天然氣液化工藝時，天然氣以多股流方式在換熱管內自下向上地流動；而在換熱器殼側，烷烴制冷劑則自上而下地流過。管、殼兩側流動都存在相變換熱，其中烷烴制冷劑以液相進入殼側吸熱汽化，先后經歷降膜流、剪切流，并以過熱氣形式流出換熱器。

在繞管式換熱器殼側換熱性能研究方面，可借鑒的相關資料較少。Aunan[4]搭建了LNG繞管式換熱器殼側實驗裝置，對烷烴制冷劑在殼側的流動與換熱性能展開了研究，指出了干度、熱流密度、壓力、質量流率對換熱的影響程度。Wu等[5]對繞管式換熱器殼側降膜流沸騰工況做出了數值模擬，基于Ansys Fluent軟件下的VOF模型開展了傳質時間松弛系數的研究，給出了乙烷、丙烷冷劑下的傳質時間松弛系數的合理取值。李劍銳等[6]利用Fluent模擬出繞管式換熱器殼側降膜流下的沸騰換熱特性，指出冷劑干度對換熱系數影響大。季鵬等[7]對多個繞管式換熱器殼側換熱關聯(lián)式進行了對比分析，指出修正后的Abadzic關聯(lián)式在天然氣預冷段最為準確。Wu等[8]以數值模擬方法研究了繞管式換熱器殼側降膜沸騰時的換熱特性，分析了利用VOF模型模擬剪切沸騰的不適用性。Ren等[9]對滾動工況下的繞管式換熱器殼側降膜流沸騰換熱進行了數值模擬研究，認為滾動振幅與周期對換熱的影響大于質量流率和壓力的作用。

鑒于Wu等[8]指出了使用VOF模型模擬繞管式換熱器殼側剪切流沸騰換熱的不適用性，所以本研究建議使用Ansys CFX軟件下的歐拉模型對繞管式換熱器殼側剪切汽化流動進行數值模擬研究，并對歐拉模型下影響流動與換熱強度的液相平均粒子直徑這一關鍵參數的取值問題做深入研究，以便能夠準確模擬換熱器殼側剪切流沸騰換熱過程，進而掌握其流動與換熱特性。

1 繞管式換熱器殼側模型的建立

1.1 殼側幾何模型

實際應用的繞管式換熱器體積巨大，而內部結構卻細致緊湊，在對其進行數值模擬研究時必將導致幾何模型網格單元數量巨大，計算耗時長且難以收斂。因此，本研究按照文獻[4]中實驗裝置建立殼側幾何模型，并通過文獻[4]中的實驗數據對模擬結果予以驗證。

文獻[4]中的實驗裝置是對繞管式換熱器殼側的簡化，換熱管分三層螺旋纏繞，中間層采用斷面完整的換熱管，而內、外兩層使用斷面為半圓的換熱管，其結構見圖2。中間纏繞層下部并行的4根換熱管壁面為加熱壁面，由管內電加熱絲控制熱流密度，其他壁面均是絕熱壁面。烷烴冷劑由換熱器頂部30個分流孔進入殼側，經過換熱器內部多排換熱管緩沖后進入到底部4根換熱管測量段。測量段的加熱壁面溫度和流道斷面溫度由熱電偶多點實測獲得。該實驗裝置的幾何參數見表1。

表1 殼側模型幾何參數幾何參數參數值換熱管外徑/mm12.00換熱管軸向間距/mm13.94換熱管徑向間距/mm15.91換熱管纏繞角/(°)7.938換熱管纏繞層數/層3 纏繞層并管數目/根3,4,5纏繞直徑/mm96,127.82,159.63 模型高度/mm160 分流孔直徑/mm10 分流孔數目/個30

鑒于該實驗裝置具有軸對稱性，因此本研究在幾何模型建立之后，將模型沿軸向切割36°作為仿真計算研究對象，其目的是降低網格數量，減少計算耗時。切割下來的殼側36°幾何模型見圖3。

1.2 控制方程

在繞管式換熱器殼側，當冷劑入口干度x>0.2時，由于氣、液相密度比很小，所以汽相占據殼側絕大部分空間，液相以離散相液滴形式存在于汽相之中以及出現(xiàn)在換熱管壁面上。在汽相吹掃作用下，換熱管壁面會出現(xiàn)暫時的干涸現(xiàn)象。此時，汽、液間作用力明顯，相間滑速比增大，在這種狀態(tài)下適宜采用歐拉模型進行模擬計算。

歐拉模型控制方程組如下[10]：

連續(xù)方程：

(1)

(2)

動量方程：

(3)

(4)

能量方程：

(5)

(6)

(1) 相間傳熱模型

單位控制單元內液相向汽相的傳熱量可通過粒子模型進行求解，將液相定為離散相粒子公式如下[10]：

Qlg=αlgAlg(Tl-Tg)

(7)

(8)

(9)

(10)

(11)

(12)

式中：αlg是汽、液相間換熱系數，W/(m2·K)；Alg是汽、液相界面面積濃度，m2/m3；dl是液相粒子直徑，m。

(2) 相間傳質模型

(13)

Qsg=αgAlg(Ts-Tg)

(14)

Qsl=αlAlg(Ts-Tl)

(15)

(16)

(17)

式中：Qsg、Qsl分別為單位控制單元內相界面向汽相和液相的傳熱量，W/m3；αg、αl分別為汽相和液相的換熱系數，W/(m2·K)；Ts為相界面溫度，K。當液相粒子采用零熱阻換熱模型時，有Ts=Tl。

(3) 相間作用力模型

(18)

CD=max{24(1+0.15Re0.687)/Re,0.44}

(19)

此外，歐拉模型在湍流計算時，汽相作為連續(xù)相使用標準k-ε湍流模型，液相以離散粒子形式使用零方程模型進行簡化處理，此時對應的湍流方程分別為[10]：

汽相湍動能方程：

(20)

汽相耗散率方程：

(21)

液相湍流黏度：

(22)

式中:σk=1.0，σε=1.3，c1ε=1.44，c2ε=1.92，σ=0.9。

1.3 邊界條件

對于繞管式換熱器殼側36°幾何模型，冷劑從上部3個分流孔進入，在底部端面流出，兩個軸向端面設為對稱邊界，中間纏繞層下部4根換熱管的外表面是加熱壁面，其余壁面屬于絕熱壁面，邊界條件設置見表2和圖3。

表2 邊界條件模型部位邊界條件類型頂部3個分流孔質量流量入口底部端面敞開流動軸向端面對稱邊界1、2中間層底部4根管壁加熱壁面給定熱流其余壁面絕熱壁面

2 離散液相平均粒子直徑推導

繞管式換熱器殼側剪切汽化流動模擬時，汽、液相分別設置成連續(xù)相和離散粒子，以便描述汽、液相間動量、能量交換，其中式(8)和式(18)的計算需要引入液相平均粒子直徑，但該參數尚無計算方法，給剪切流汽化模擬帶來困難，因此，需對液滴平均粒子直徑作理論研究。

殼側剪切汽化流動時，汽、液相溫度并不相同，汽相溫度高于液相溫度，流場內彌散狀液滴通過相界面從汽相獲取熱量實現(xiàn)汽化。液滴平均粒子直徑可以通過微元控制體能量平衡推導出來。流場內部微元控制體見圖4，圖中藍色部分表示液相，其余部分表示汽相。

假設沿x、y、z方向流入微元控制體的單位體積能量為hV,x、hV,y、hV,z，相應的流速為u、v、w。則x方向流入能量為：

Hx=hV,xudydz

(23)

同理，y、z方向流入能量為：

Hy=hV,yvdxdz

(24)

Hz=hV,zwdxdy

(25)

x方向流出能量為：

(26)

同理，y、z方向流出能量為：

(27)

(28)

當汽、液相間換熱采用兩熱阻模型且液相側使用零熱阻時，有αlg=αg，Ts=Tl。穩(wěn)態(tài)時，微元控制體內汽相向液相傳遞的熱量即是外界應向微元控制體輸入的汽化熱量，也即是微元體的凈導入熱量，該熱量為：

Qlg=αgAlg(Tg-Tl)dxdydz

(29)

穩(wěn)定流動下根據微元控制體能量平衡可得到：

Hx+Hy+Hz+Qlg=Hx+dx+Hy+dy+Hz+dz

(30)

因dx、dy、dz均趨于0，所以有：

hV,x=hV,y=hV,z=HV

(31)

將式(23)～式(29)及式(31)代入至式(30)，簡化后得到：

(32)

將式(8)及式(12)代入至式(32)中，得到：

(33)

在流場范圍Ω內對式(33)進行三重積分，有：

(34)

式(34)左側積分利用高斯定理可轉化為：

(35)

式中:n為流場邊界曲面數量，?x、?y、?z為邊界曲面外法線與坐標軸夾角。將殼側模型對稱面視為壁面并且以無滑移處理壁面后，式(35)右側可簡化為：

(36)

式中:Aheat為換熱管加熱壁面的面積，m2；q為換熱管壁面熱流密度，W/m2。

式(34)右側積分時利用中值定理可變?yōu)椋?/p>

(37)

(38)

將式(36)和式(38)代入至式(34)中，變形后得到：

(39)

3 熱相變比的應用

(40)

式中:Δx為殼側沸騰時的干度增加量。從定義式可以看出，熱相變比ζ表征了向殼側輸入熱量中，用于液滴汽化的熱量占總熱量的比例。

將式(40)變形后代入至式(39)，得到：

(41)

表3 乙烷液相平均粒子直徑及其工況參數壓力/MPa入口干度質量流率/(kg·(m2·s)-1)熱流密度/(W·m-2)熱相變比平均粒子直徑/m換熱偏差/%0.20.2144.947 839.130.967 60.000 082-15.540.3152.587 839.130.971 60.000 04010.760.4243.817 839.130.974 40.000 033-19.980.5361.067 839.130.978 10.000 022-16.470.6251.736 315.290.983 10.000 014-18.780.7251.164 822.870.991 60.000 007-15.180.8251.737 839.130.994 70.000 00414.390.9142.977 839.130.997 30.000 002-14.10.30.2161.347 839.130.947 30.000 12117.250.3260.774 744.320.958 40.000 06219.70.4260.773 157.650.959 50.000 044-2.980.5243.536 315.290.969 40.000 0335.030.6244.107 839.130.980 30.000 0202.270.7343.537 839.130.987 70.000 011-16.190.8742.687 839.130.990 30.000 006-17.890.40.2154.274 822.870.941 50.000 13713.050.3261.342 364.310.949 80.000 07614.750.4260.492 364.310.956 20.000 0565.630.5343.534 744.320.961 10.000 0408.310.6259.647 839.130.967 10.000 0245.480.7260.217 839.130.972 70.000 016-9.460.8157.957 839.130.984 20.000 009-13.970.9159.367 839.130.992 10.000 0041.66

本研究依據文獻[4]的實驗數據，研究了乙烷在變工況條件下做剪切汽化流動時熱相變比的取值問題。通過設定不同的ζ值，將模擬工況下的換熱結果與實驗數據做比對來確定合理的熱相變比。模擬工況共選取23個，涵蓋了變壓力、變干度、變熱流、變流量的情況，經過百余次計算，確定出來的ζ值見表3。

通過表3可以看出，殼側冷劑入口干度(x)對熱相變比的影響最為敏感，兩者間呈同向變化，即干度增大時要求熱相變比相應變大，從而減小液相平均粒子直徑，以滿足提升換熱系數的要求。而當質量流率(G)和熱流密度(g)變化時，兩因素對熱相變比的影響不明顯，即無法改變冷劑干度對熱相變比的影響走勢，因此，熱相變比的取值可以忽略質量流率和熱流密度的影響。

圖5是乙烷23個模擬工況按定壓條件繪制出的熱相變比變化圖。從圖5可以看出，熱相變比受壓力的影響較為明顯。在殼側入口干度相同的條件下，壓力越小則要求熱相變比越大。其原因在于，高壓向低壓轉變時,相對于其他物性參數而言，汽相密度變化最為明顯，汽相密度變小將導致汽相流速增大，隨之而來的相間動量、能量傳輸增強，必然要求液相平均粒子直徑進一步變小，即熱相變比需要相應增大。

鑒于變壓力時，汽相密度對熱相變比有較大的影響，因此本研究對熱相變比進行公式擬合時自變量包含了干度x和汽相密度ρg兩個物理量，根據模擬工況下的計算數據擬合出的熱相變比公式見式(42)：

(42)

a1=-0.155，a2=0.088，a3=7.606，a4=-0.207，a5=0.034，a6=7.564；

b1=8.122，b2=-14.247，b3=7.879，b4=-13.859

式中：ρ/7.36為變壓力時汽相密度發(fā)生改變而做出的修正，以乙烷0.4 MPa下汽相密度7.36 kg/m3作為基準進行熱相變比修正。上式擬合程度較好，擬合相關系數R=0.993 9。

圖6和圖7分別為乙烷于換熱器殼側做剪切汽化流動時在軸向剖面、中間纏繞層壁面處的汽相體積分數模擬云圖，其模擬工況參數p=0.4 MPa、G=59.36 kg/(m2·s)、x=0.91、q=7 839.13 W/m2。

通過圖6和圖7可以看出，殼側冷劑高干度條件下的兩相流動趨于汽態(tài)流動，液相是以極細的液滴粒子形式彌散于汽相當中，因此，流道內每個網格單元的汽相體積分數接近于1。在中間纏繞層壁面上，由于液相潤濕效應而使壁面局部位置處的液相體積分數略有增加。在近似于汽態(tài)流動、換熱的條件下，只有將熱相變比增至0.992 1，液滴平均粒子直徑減小到0.000 004 m，才能準確地模擬出換熱系數高達12 153 W/(m2·K)的汽化工況，并使模擬偏差低至1.66%。

4 結論

(1) 兩相流歐拉模型適用于繞管式換熱器殼側剪

切汽化流動時的數值模擬研究，通過理論推導以及引入熱相變比參數，可以確定出離散液相的平均粒子直徑，進而實現(xiàn)對殼側剪切流汽化換熱的準確模擬，其換熱模擬偏差在±20%范圍內。

(2) 以乙烷在繞管式換熱器殼側剪切汽化流動為例，通過換熱模擬結果與實驗數據相對比的方式，確定出變工況條件下的熱相變比取值。經分析發(fā)現(xiàn)，殼側冷劑入口干度對熱相變比的影響最為強烈，冷劑工作壓力對熱相變比的影響次之，而殼側質量流率、熱流密度對熱相變比的影響可以忽略不計。

(3) 在分析熱相變比影響因素之后，選取殼側冷劑入口干度和汽相密度作為自變量對熱相變比進行公式擬合，給出了相應的計算式，其擬合程度較好，相關系數R=0.993 9。

(4) 數值模擬研究結果表明，通過熱相變比確定離散相平均粒子直徑的方法，可以使歐拉模型準確地模擬出兩相流相變換熱過程，解決了相變模擬失真這一突出矛盾。該方法對兩相流仿真計算具有實用價值。

符號表符號物理量名稱單位Alg汽、液相界面面積濃度m2/m3Aheat換熱管受熱面積m2Cp定壓比熱J/(kg·K)d離散相粒子直徑mF相間作用力N/m3g重力加速度m/s2h比焓J/kg

續(xù)符號表符號物理量名稱單位hV單位體積焓J/m3H單位時間總焓Wk湍動能m2/s2m殼側質量流量kg/m單位控制容積相間傳質量kg/(m3·s)Nu努塞爾準則p壓力PaPr普朗特準則q熱流密度W/ m2Q傳熱量W/ m3Re雷諾準則T溫度Ku、v、w直角坐標系下速度分量m/sv速度矢量m/sV體積m3x干度α換熱系數W/(m2·K)β體積分數?x、?y、?z曲面外法線與坐標軸夾角(°)ρ密度kg/m3λ導熱系數W/(m·k)ε湍流耗散率m2/s3τ時間sζ熱相變比μ動力黏度N·s/m2下角標符號含義符號含義g汽相gl汽相作用于液相l(xiāng)液相l(xiāng)g液相作用于汽相model殼側幾何模型s相界面sat 飽和t湍流