繞管式換熱器殼側(cè)沸騰模擬時的時間松弛參數(shù)研究

吳志勇， 劉 洋， 高 陽， 蔡偉華， 姜益強

(1.遼寧石油化工大學 石油天然氣工程學院，遼寧 撫順 113001；2.哈爾濱工業(yè)大學 能源科學與工程學院，黑龍江 哈爾濱 150001；3.哈爾濱工業(yè)大學 市政與環(huán)境工程學院，黑龍江 哈爾濱 150090)

繞管式換熱器殼側(cè)沸騰模擬時的時間松弛參數(shù)研究

吳志勇1， 劉 洋1， 高 陽1， 蔡偉華2， 姜益強3

(1.遼寧石油化工大學 石油天然氣工程學院，遼寧 撫順 113001；2.哈爾濱工業(yè)大學 能源科學與工程學院，黑龍江 哈爾濱 150001；3.哈爾濱工業(yè)大學 市政與環(huán)境工程學院，黑龍江 哈爾濱 150090)

利用ANSYS Fluent模擬軟件下的流體體積(VOF)模型，對液化天然氣用繞管式換熱器殼側(cè)沸騰過程進行數(shù)值模擬，針對控制殼側(cè)沸騰強度的時間松弛參數(shù)取值問題展開研究，提出采用汽化當量法確定時間松弛參數(shù)。結(jié)果表明，殼側(cè)制冷劑為乙烷、丙烷時的時間松弛參數(shù)為3 s-1，并由沸騰換熱模擬偏差得到驗證。汽化當量法可用于各種沸騰仿真計算，對準確模擬傳質(zhì)過程具有使用價值。

天然氣液化； VOF模型； 沸騰傳質(zhì)； 汽化當量法； 汽相體積分數(shù)

1 繞管式換熱器介紹

繞管式換熱器作為一種特殊類型的換熱設(shè)備廣泛應用于天然氣液化工藝，其特點是結(jié)構(gòu)緊湊、耐低溫、抗高壓以及多種介質(zhì)可以同時換熱。這種換熱器是由許多小管徑換熱管以螺旋形狀分層纏繞在換熱器芯筒上而制成的，相鄰兩層換熱管的旋向相反，層間采用隔條保持間距[1-2]。繞管式換熱器內(nèi)部結(jié)構(gòu)見圖1[3]。

對于天然氣液化工藝，殼側(cè)烷烴制冷劑首先進入管側(cè)進行過冷，天然氣與需要過冷的制冷劑在不同的換熱管內(nèi)形成管側(cè)多股流，并一同向上流動。過冷后的管側(cè)制冷劑在換熱器外部經(jīng)節(jié)流后導入到殼側(cè)并向下流動，對管側(cè)天然氣和需要過冷的制冷劑實現(xiàn)制冷。烷烴制冷劑在殼側(cè)自上向下流動時，吸收管側(cè)熱量而汽化，先后經(jīng)歷兩相流動和過熱流動兩種狀態(tài)。

在繞管式換熱器流動與換熱研究領(lǐng)域，已公開發(fā)表的文獻比較少。文獻[1]介紹了繞管式換熱器結(jié)構(gòu)參數(shù)間的耦合關(guān)系，給出了殼側(cè)流通面積的計算方法。B.O.Neeraas[4]建立了繞管式換熱器管側(cè)實驗裝置，以實驗方式獲得了烷烴混合工質(zhì)在管側(cè)冷凝時的壓降、換熱數(shù)據(jù)，并對管側(cè)流動、換熱影響因素進行了分析。B.Aunan[5]建立了液化天然氣用繞管式換熱器殼側(cè)實驗臺，對乙烷、丙烷在殼側(cè)的流動與換熱特性進行了實驗研究。M.Moawed[6]以實驗方式研究了定熱流密度下纏繞管曲率比、扭曲比與殼側(cè)換熱系數(shù)之間的關(guān)系，指出當纏繞管曲率比、扭曲比變大時殼側(cè)單相換熱系數(shù)明顯增大。

圖1 繞管式換熱器內(nèi)部結(jié)構(gòu)

在繞管式換熱器數(shù)值模擬研究方面，文獻[2]研究了殼側(cè)氣相流動時的壓降與換熱特點，指出了換熱器優(yōu)化設(shè)計的方向。賈金才[7]建立了繞管式換熱器三維計算模型，研究了換熱管管徑、纏繞角、徑向比、軸向比對殼側(cè)換熱系數(shù)的影響，但所有模擬工況均為氣相流動，與殼側(cè)兩相流動為主的實際工況存在差別。X.Lu等[8]針對繞管式換熱器殼側(cè)氣態(tài)換熱過程展開了實驗與計算機仿真研究工作，建立了具有三層繞管的實驗裝置以及CFD仿真計算模型，實驗證實了所用CFD模型具有較高的計算精度。

縱觀以上文獻發(fā)現(xiàn)，繞管式換熱器流動與換熱性能研究可以通過實驗和數(shù)值模擬兩種方法來實現(xiàn)，但數(shù)值模擬均是在單相流動條件下進行的，其原因是在數(shù)值模擬條件下相變時的傳質(zhì)強度難以準確把握，導致相變換熱模擬難度加大。本文對控制繞管式換熱器殼側(cè)沸騰強度的時間松弛參數(shù)取值問題做深入研究，期望對相變模擬問題具有指導作用。

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

2.1 殼側(cè)幾何模型

液化天然氣用繞管式換熱器整體尺寸巨大，其外部直徑為3～5 m，高度為10～50 m，加之內(nèi)部結(jié)構(gòu)復雜，極不利于開展數(shù)值模擬研究。本文在研究過程中以文獻[5]中的實驗數(shù)據(jù)作為模擬計算檢驗標準，因此所建殼側(cè)幾何模型與文獻[5]中的實驗設(shè)備一致。

文獻[5]中的實驗設(shè)備為LNG繞管式換熱器的簡化模型，由3層換熱管以交替的纏繞方向纏繞而成，每層的換熱管并管纏繞數(shù)目由內(nèi)向外依次為3、4、5。由完整的換熱管纏繞而成的中間層及兩側(cè)流道為實驗研究的核心部位；內(nèi)、外兩側(cè)纏繞層均由半剖換熱管纏繞而成，并作為內(nèi)、外邊壁使用。換熱器簡化模型如圖2所示。文獻[5]對管-殼兩側(cè)換熱實行解耦，以給定熱流方式研究殼側(cè)換熱特性，其中內(nèi)、外纏繞層及中間纏繞層上部為絕熱壁面，中間纏繞層下部4排換熱管為受熱壁面，通過管內(nèi)電加熱絲實現(xiàn)定熱流加熱。烷烴制冷劑由頂部30個直徑為10 mm的分流孔進入殼側(cè)，經(jīng)過上方多排纏繞管束緩沖后進入到壓降與換熱測試區(qū)。

圖2 換熱器簡化模型

按照文獻[5]中的換熱器幾何參數(shù)，本文利用UG NX 6.0建模軟件建立了殼側(cè)模型，其幾何參數(shù)見表1。

表1 殼側(cè)幾何模型參數(shù)

由于模擬相變換熱需要精細的邊界層網(wǎng)格，因此在保證計算精度的前提下，殼側(cè)模型網(wǎng)格數(shù)量達到1 000萬以上，導致模擬計算耗時巨增?？紤]到該模型具有軸對稱性，本文將模型沿軸向切割36°作為模擬研究對象，使網(wǎng)格數(shù)量降至100萬左右，此時網(wǎng)格數(shù)量與計算精度能夠得以兼顧。切割下來的殼側(cè)36°幾何模型如圖3所示。

(a) 正視圖 (b) 側(cè)視圖

2.2 控制方程

烷烴制冷劑流入殼側(cè)時，換熱管管束被液膜覆蓋，殼側(cè)形成降膜流動，適宜采用兩相流VOF模型進行計算。ANSYS Fluent下的VOF模型控制方程見式(1)—(4)[9-10]。

連續(xù)方程：

(1)

(2)

動量方程：

(3)

式中，p為壓力，Pa；g為重力加速度，m/s2；μ為動力黏度，Pa·s。

能量方程：

(4)

式中，μg、μl為汽、液相動力黏度，kg/(m·s)；E為比內(nèi)能，J/kg；Γe為能量傳輸源項，W/m3。

控制方程組中的連續(xù)方程和能量方程源項通過下列模型進行計算：

(1)相變傳質(zhì)模型。根據(jù)W.H.Lee[11]的傳質(zhì)模型，當流體溫度高于汽化溫度即T≥Tsat時，單位控制單元內(nèi)由液相轉(zhuǎn)變?yōu)槠嗟膫髻|(zhì)量(汽化量)即為Γm，并有：

(5)

反之，當T

(6)

式中，βm為相變傳質(zhì)過程中的時間松弛參數(shù)，s-1，其為控制汽化強度的決定性參數(shù)，其取值問題是本文研究的重點。

(2)能量傳輸模型。當汽、液相間傳質(zhì)量Γm確定后，相應的能量傳輸源項Γe可用式(7)來表達：

Γe=-hLHΓm

(7)

式中，hLH為流體汽化潛熱，J/kg，式中的負號表示需要外界向控制單元輸入熱量。

另外，VOF模型用于湍流時，本文使用RNGk-ε湍流應力模型作為封閉方程，其形式見式(8)—(9)[12]。

湍動能方程：

(8)

耗散率方程：

(9)

式中，k為湍動能，m2/s2；ε為耗散率，m2/s3；Rk為平均速度梯度引起的湍動能產(chǎn)生項，kg/(m·s3)；Eij為湍流時均應變率，s-1；cμ、ck、cε、c1ε、φ0、φ0均為常數(shù)，且cμ=0.084 5，ck=1.39，cε=1.39，c1ε=1.42，c2ε=1.68，φ0=4.377。

2.3 邊界條件

對于殼側(cè)36°幾何模型，流體從頂部3個分流孔流入，在底部端面流出，兩側(cè)剖面均近似處理為對稱邊界，中間纏繞層底部4排換熱管為加熱壁面，其余壁面均為絕熱壁面，邊界條件見表2。

表2 邊界條件

2.4 網(wǎng)格劃分及無關(guān)性驗證

使用ANSYS ICEM CFD網(wǎng)格劃分軟件包對殼側(cè) 36°幾何模型進行網(wǎng)格劃分，采用非結(jié)構(gòu)化四面體及混合網(wǎng)格，先生成面網(wǎng)格，進而自適應生成體網(wǎng)格。在模型中，纏繞管軸向縫隙 1.94 mm為最小幾何尺寸，故以此確定面網(wǎng)格尺寸大小。在確保軸向管縫處存在 3 行網(wǎng)格且網(wǎng)格邊長比近似于 2 的前提下，做出3套網(wǎng)格，網(wǎng)格參數(shù)見表 3。此外，為了準確計算換熱系數(shù)，在換熱管受熱壁面設(shè)置精細邊界層，邊界層初始厚度 0.005 mm，增長因子1.2，共計10 層邊界層。

表3 殼側(cè)模型網(wǎng)格方案的網(wǎng)格參數(shù)

為了考察網(wǎng)格疏密程度對數(shù)值求解的影響并確定網(wǎng)格劃分方案，本文以文獻[5]中的丙烷降膜沸騰流動實驗數(shù)據(jù)作為依據(jù)，進行模擬結(jié)果對比。實驗工況參數(shù)為：p=0.2 MPa；質(zhì)量流率G=48.6 kg/(m2·s)；入口干度x=0；q=3 951 W/m2。數(shù)值模擬的換熱系數(shù)及其偏差即網(wǎng)格無關(guān)性驗證數(shù)據(jù)見表 4。

表4 網(wǎng)格無關(guān)性驗證數(shù)據(jù)

由表 4 可知，隨著網(wǎng)格數(shù)量的增多，換熱系數(shù)模擬偏差略有下降，但變化趨勢趨緩。網(wǎng)格方案 1 的網(wǎng)格數(shù)量明顯高于網(wǎng)格方案 2，導致網(wǎng)格方案 1 的計算時間遠大于網(wǎng)格方案 2 所需時間。所以，在綜合考量計算精度與耗時的前提下，本文將網(wǎng)格方案 2 作為最終的網(wǎng)格劃分方案。

3 沸騰模擬下時間松弛參數(shù)的確定

VOF模型中的連續(xù)方程含有質(zhì)量源項，由式(5)—(6)可知，傳質(zhì)量的多少主要取決于時間松弛參數(shù)。數(shù)值模擬時，若選取的時間松弛參數(shù)過小，則沸騰汽化傳質(zhì)量少；反之，若選取的時間松弛參數(shù)過大，則必將導致汽化量增大。因此，合理確定烷烴類的時間松弛參數(shù)是準確模擬繞管式換熱器殼側(cè)沸騰這一實際過程的關(guān)鍵。

縱觀國內(nèi)外文獻，目前尚無沸騰傳質(zhì)時間松弛參數(shù)相關(guān)研究內(nèi)容。鑒于此，本文提出采用汽化當量法確定時間松弛參數(shù)。汽化當量法是指相對于任一沸騰工況，存在一個汽化當量工況(Equivalent-Evaporation Simulation， EES)與之對應，該工況是以沸騰工況下的理論汽化量及兩相總流量為依據(jù)而構(gòu)造出來的，屬于絕熱工況，其汽、液相入口質(zhì)量流量分別為：

(10)

(11)

汽化當量法的涵義是：在沸騰工況下，工質(zhì)汽化后在殼側(cè)出口的汽化量應等于由熱平衡計算出來的理論汽化量；并且，在兩相總流量不變的前提下，沸騰工況下的殼側(cè)出口汽體量應與汽化當量工況下的汽體量相同。

在VOF模型下，殼側(cè)任一截面處關(guān)于分相質(zhì)量流量及干度的計算數(shù)據(jù)無法讀取，僅能得到汽、液相體積分數(shù)計算結(jié)果，所以汽化當量法只能通過比對出口端流道截面汽相體積分數(shù)的方式來使用。具體而言，在同一沸騰工況下，選取不同的時間松弛參數(shù)進行沸騰模擬計算，當汽化量不同時流道出口處汽相體積分數(shù)必定不同，而能反映真實汽化程度的時間松弛參數(shù)值可通過汽化當量工況下的流道出口汽相體積分數(shù)來決定。

本文參照文獻[5]的實驗工況，對乙烷、丙烷的時間松弛參數(shù)進行研究。對于同一沸騰工況，時間松弛參數(shù)分別取0.1、1.0、2.0、3.0、4.0、5.0 s-1，經(jīng)過動態(tài)模擬并對收斂結(jié)果求均值后在流道截面2處(見圖3)得到不同的汽相體積分數(shù)，然后與汽化當量工況下的對應值做比較從而確定時間松弛參數(shù)值。

模擬工況分為兩類，一類用于確定時間松弛參數(shù)值；另一類覆蓋多壓力、多熱流、多流率情況，用于檢驗時間松弛參數(shù)值的通用性。在兩類工況下，殼側(cè)入口干度均為0，其余參數(shù)見表5—6。

表5 時間松弛參數(shù)確定工況

表6 時間松弛參數(shù)檢驗工況

丙烷和乙烷在流道截面2處的汽相體積分數(shù)隨時間松弛參數(shù)(βm)的變化曲線如圖4所示。丙烷與乙烷相應的EES模擬結(jié)果也示于圖4中。由圖4可知，出口汽相體積分數(shù)隨著βm的增大而變大，即流道內(nèi)汽化強度增強；βm=3 s-1對丙烷比較適宜，而βm=4 s-1對乙烷較適宜。

圖4 出口汽相體積分數(shù)隨βm變化曲線

當時間松弛參數(shù)的取值不同時，加熱壁面汽相體積分數(shù)及壁溫也會不同，由此可能會導致模擬沸騰換熱系數(shù)發(fā)生相應的變化，因此可通過考核沸騰換熱系數(shù)的方法輔助確定時間松弛參數(shù)。βm分別為2、3、4 s-1時換熱系數(shù)偏差如圖5所示。換熱系數(shù)偏差是根據(jù)模擬結(jié)果和實驗數(shù)據(jù)而得到的，換熱系數(shù)計算公式見式(12)。

α=q/[Tw,heat-0.5(T1+T2)]

(12)

式中，Tw,heat為加熱管的壁溫，K；T1、T2分別為流道截面1、2的平均溫度，K。

圖5 沸騰換熱系數(shù)偏差隨βm變化

由圖5可知，當βm=3 s-1時，丙烷與乙烷均取得了最小換熱系數(shù)偏差，鑒于乙烷在βm=3 s-1時出口汽相體積分數(shù)與EES對應值偏差不大，故將丙烷和乙烷的時間松弛參數(shù)統(tǒng)一定為βm=3 s-1，其值將在檢驗工況下被進一步考察。

βm=3 s-1時，丙烷、乙烷在各檢驗工況及其相應EES的出口汽相體積分數(shù)如圖6所示；丙烷、乙烷在各檢驗工況下出口汽相體積分數(shù)與EES對應值之間的偏差如圖7所示。相對偏差以EES出口汽相體積分數(shù)為基準。由圖7可知，丙烷的相對偏差為-15%～20%，乙烷的相對偏差為-20%～25%，表明βm=3 s-1比較理想。

(a) 丙烷

(b) 乙烷

(a) 丙烷

(b) 乙烷

丙烷和乙烷在各檢驗工況下的沸騰換熱系數(shù)偏差如圖8所示。由圖8可知，兩種工質(zhì)的換熱系數(shù)偏差均控制在-20%～10%，說明從換熱角度來看βm=3 s-1也是合適的。

(a) 丙烷

(b) 乙烷

經(jīng)過汽化當量法以及換熱偏差兩方面檢驗后，本文認為βm=3 s-1對以乙烷和丙烷為代表的殼側(cè)冷劑是適宜的，并且在變壓力、變熱流密度以及變流率情況下具有通用性，在VOF模型下使用此值能夠比較真實地反映殼側(cè)沸騰情況。

4 沸騰模擬下的殼側(cè)可視化圖像

對各沸騰工況設(shè)置時間松弛參數(shù)并模擬求解后，即可得到殼側(cè)不同部位的可視化圖像。以P1和E1工況為例，在時間松弛參數(shù)取值不同時，殼側(cè)不同部位的汽相體積分數(shù)分布見表7?？梢暬课话〒Q熱管中間纏繞層、殼側(cè)軸向剖面以及流道周向剖面。

從表7可以看出，隨著時間松弛參數(shù)變大，殼側(cè)汽化量增加，即時間松弛參數(shù)控制殼側(cè)沸騰時的汽化強度；同時，通過對流道出口汽相體積分數(shù)以及沸騰換熱系數(shù)偏差的分析，可以確定βm=3 s-1的取值比較適宜。因此，可以認為βm=3 s-1時的可視化圖像能夠反映實際沸騰過程。

表7 βm變化時殼側(cè)不同部位汽相體積分數(shù)分布

5 結(jié) 論

(1)汽化當量法能夠用于確定時間松弛參數(shù)，通過比對模擬工況與汽化當量工況下的汽化量可以明確時間松弛參數(shù)的取值。

(2)數(shù)值模擬乙烷、丙烷低溫制冷劑在繞管式換熱器殼側(cè)沸騰時，經(jīng)過比較模擬工況與汽化當量工況的汽相體積分數(shù)，可以確定時間松弛參數(shù)βm=3 s-1比較適宜，并且該值適用于變壓力、變熱流、變流率的不同工況。

(3)時間松弛參數(shù)取值變化時，由于汽化強度不同，加熱壁面的壁溫也會不同，由此導致模擬的沸騰換熱系數(shù)發(fā)生相應變化。因此，可以通過考核沸騰換熱系數(shù)的方法來輔助確定時間松弛參數(shù)。在時間松弛參數(shù)βm=3 s-1的條件下，不同工況下的沸騰換熱系數(shù)模擬偏差均為-20%～10%，進一步證實了乙烷、丙烷時間松弛參數(shù)βm=3 s-1的合理性。

[1] 吳志勇，陳杰，浦暉，等.LNG繞管式換熱器結(jié)構(gòu)與流通參數(shù)計算[J].煤氣與熱力，2014，34(3)：42-47.

[2] 吳志勇，陳杰，浦暉，等.LNG繞管式換熱器殼側(cè)過熱態(tài)流動的數(shù)值模擬[J].煤氣與熱力，2014，34(8)：6-11.

[3] Steffen H.Near-optimaloperation of LNG liquefaction processes by means of regulation[D].Berlin：Berlin Institute of Technology，2011.

[4] Neeraas B O.Condensation of hydrocarbon mixtures in ciol-wound LNG heat exchangers tube-side heat transfer and pressure drop[D].Trodheim:The Norwegian Institute of Technology，1993.

[5] Aunan B.Shell-sideheat transfer and pressure drop in coil-wound LNG heat exchanger[D].Trodheim：Norwegian University of Science and Technology，2000.

[6] Moawed M.Experimental study of forced convection from helical coiled tubes with different parameters[J].Energy Conversion and Management，2011，52(2)：1150-1156.

[7] 賈金才.幾何參數(shù)對繞管式換熱器傳熱特性影響的數(shù)值研究[J].流體機械，2011，39(8)：33-37.

[8] Lu X，Du X P，Zeng M，et al. Experimental and numerical investigation on shell-side performance of multilayer spiral-wound heat exchangers[J].Chemical Engineering Transactions，2013，35：445-450.

[9] Yang Z，Peng X F，Ye P.Numerical and experimental investigation of two phase flow during boiling in a coiled tube[J].International Journal of Heat and Mass Transfer，2008，51(5)：1003-1016.

[10] 李霞，馬貴陽，杜明俊，等. LNG球罐分層翻滾汽化過程的數(shù)值分析[J].遼寧石油化工大學學報，2010，30(4)：42-44.

[11] Lee W H.A pressure iteration scheme for two-phase flow modeling[M]. Washington,DC：Hemisphere Publishing，1980：407-431.

[12] Orszag S A，Yakhot V，F(xiàn)lannery W S，et al.Renormalization group modeling and turbulence simulations[C]. Proceedings of the International Conference on Near-wall Turbulent Flows.Arizona：Elsevier Science Publishers B.V.，1993.

An Investigation on the Time Relaxation Parameter Used in Boiling Simulation of Helical Coil Heat Exchanger Shell-Side

Wu Zhiyong1, Liu Yang1， Gao Yang1， Cai Weihua2, Jiang Yiqiang3

(1.CollegeofPetroleumEngineering，LiaoningShihuaUniversity，F(xiàn)ushunLiaoning113001，China；2.SchoolofEnergyScienceandEngineering，HarbinInstituteofTechnology，HarbinHeilongjiang150001，China；3.SchoolofMunicipalandEnvironmentalEngineering，HarbinInstituteofTechnology，HarbinHeilongjiang150090，China)

The volume of fluid (VOF) model which was a multi-phase flow model in ANSYS Fluent software was used to simulate the shell-side boiling of helical coil heat exchanger (HCHE), the time relaxation parameter which dominated mass transfer was studied by equivalent-evaporating method, and the value of time relaxation parameter was recommended as 3 s-1for ethane and propane. The value of time relaxation parameter was also verified by heat transfer deviation. The equivalent-evaporating method which was proposed could be used in all kinds of boiling simulations, and had good application value for mass transfer calculation.

Nature gas liquefaction； VOF model； Boiling mass transfer； Equivalent-evaporating method； Vapor volume fraction

1672-6952(2017)05-0053-08

投稿網(wǎng)址：http://journal.lnpu.edu.cn

2016-10-27

2016-11-28

工信部高科技項目“大型LNG繞管式換熱器研制”(2013418)。

吳志勇(1977-)，男，博士，講師，從事天然氣液化技術(shù)與設(shè)備研究；E-mail：1123589558@qq.com。

TE821；TK172

A

10.3969/j.issn.1672-6952.2017.05.011

(編輯 宋錦玉)