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

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

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

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

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

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

1 繞管式換熱器殼側(cè)模型的建立

1.1 殼側(cè)幾何模型

實際應(yīng)用的繞管式換熱器體積巨大，而內(nèi)部結(jié)構(gòu)卻細(xì)致緊湊，在對其進行數(shù)值模擬研究時必將導(dǎo)致幾何模型網(wǎng)格單元數(shù)量巨大，計算耗時長且難以收斂。因此，本研究按照文獻(xiàn)[4]中實驗裝置建立殼側(cè)幾何模型，并通過文獻(xiàn)[4]中的實驗數(shù)據(jù)對模擬結(jié)果予以驗證。

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

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

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

1.2 控制方程

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

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

連續(xù)方程：

(1)

(2)

動量方程：

(3)

(4)

能量方程：

(5)

(6)

(1) 相間傳熱模型

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

Qlg=αlgAlg(Tl-Tg)

(7)

(8)

(9)

(10)

(11)

(12)

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

(2) 相間傳質(zhì)模型

(13)

Qsg=αgAlg(Ts-Tg)

(14)

Qsl=αlAlg(Ts-Tl)

(15)

(16)

(17)

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

(3) 相間作用力模型

(18)

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

(19)

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

汽相湍動能方程：

(20)

汽相耗散率方程：

(21)

液相湍流黏度：

(22)

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

1.3 邊界條件

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

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

2 離散液相平均粒子直徑推導(dǎo)

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

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

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

Hx=hV,xudydz

(23)

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

Hy=hV,yvdxdz

(24)

Hz=hV,zwdxdy

(25)

x方向流出能量為：

(26)

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

(27)

(28)

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

Qlg=αgAlg(Tg-Tl)dxdydz

(29)

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

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)

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

(34)

式(34)左側(cè)積分利用高斯定理可轉(zhuǎn)化為：

(35)

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

(36)

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

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

(37)

(38)

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

(39)

3 熱相變比的應(yīng)用

(40)

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

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

(41)

表3 乙烷液相平均粒子直徑及其工況參數(shù)壓力/MPa入口干度質(zhì)量流率/(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

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

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

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

鑒于變壓力時，汽相密度對熱相變比有較大的影響，因此本研究對熱相變比進行公式擬合時自變量包含了干度x和汽相密度ρg兩個物理量，根據(jù)模擬工況下的計算數(shù)據(jù)擬合出的熱相變比公式見式(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作為基準(zhǔn)進行熱相變比修正。上式擬合程度較好，擬合相關(guān)系數(shù)R=0.993 9。

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

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

4 結(jié)論

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

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

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

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

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

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

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