LNG氣化器管內氣泡行為及兩相流型轉換

徐少杰, 高文學, 嚴榮松, 王 艷, 楊 林, 張 歡

(1.天津大學 環境科學與工程學院，天津300072;2.中國市政工程華北設計研究總院有限公司城市燃氣熱力研究院，天津300384)

1 概述

在中小型LNG供氣廠站中，空溫式氣化器是比較常見的一種氣化設備，其利用空氣作為熱源，對翅片管內的LNG進行加熱。目前關于LNG空溫式氣化器的研究，大多數集中于單根星型翅片管的傳熱問題，內容從管內的流動沸騰到管外的自然對流與濕空氣結霜[1-2]。雖然不同流型的轉換機理對于流動沸騰傳熱特性的研究具有重要意義，但是目前對于翅片管內天然氣氣泡行為以及氣液兩相流型轉換的研究遠遠不足。豎直管內的流動沸騰沿管程方向依次表現為泡狀流動、彈狀流動、環狀流動以及彌散狀流動[3]。李祥東等[4]通過分析研究兩相流型常用理論，認為常用理論很難準確地判斷低溫流體的流型變化。劉珊珊等[5]使用經典Lee模型[6]研究了兩相流型變化與局部傳熱系數之間的聯系，表明局部換熱系數隨氣相分數呈現先增后減的趨勢。文玨[7]的研究結果發現，與Lee模型相比，S.Hardt模型[8]對氣液的蒸發沸騰過程描述更為準確。

雖然近年來采用計算流體力學方法對空溫式氣化器內外場耦合傳熱特性的研究已有諸多報道[9]，但是更多是從宏觀層面探索空溫式氣化器的傳熱規律，然而，氣泡在沸騰表面的形成、生長以及脫離對于流動沸騰傳熱有非常顯著的影響。為了研究管內LNG沸騰過程氣泡的行為和兩相流型的轉化，本文借助于Fluent 19.0軟件，采用流體體積法(Volume of Fluid，簡稱VOF)模型并耦合S.Hardt模型，計算空溫式氣化器翅片管內的氣液兩相流動沸騰過程，分析不同LNG入口流速、LNG入口溫度以及管外表面傳熱系數等條件對管內流動沸騰換熱的影響。

2 問題描述及假定說明

2.1 問題描述

空溫式氣化器由呈陣列形式布置的鋁合金翅片管簇組成，LNG由翅片管底部進入，自下而上流動，在此過程中不斷從管壁吸收來自環境空氣的熱量。當LNG從入口溫度上升到飽和溫度后，管內流體進入蒸發沸騰階段，形成氣液兩相流動。隨著管內流體溫度升高，天然氣的體積分數不斷增大，管內氣液兩相流型將出現顯著變化。本文的研究對象為翅片管內流動沸騰過程中氣泡行為和兩相流型，分析不同LNG入口流速、LNG入口溫度以及管外表面傳熱系數等條件下，兩相流型的轉換特性。

2.2 假定說明

① 管道導熱為各向同性。

② 管道出口天然氣體積分數不隨時間改變時，認為管內流動沸騰傳熱過程發展到穩態。

③ 管外空氣表面傳熱系數沿管長方向為定值。

④ 管外空氣溫度為定值。

3 幾何模型描述

本文研究的重點為翅片管內流動沸騰過程中氣泡行為和兩相流型的變化，管外翅片對于管內流動沸騰過程的影響僅僅在于傳熱問題。因此，翅片所帶來的傳熱影響將通過等效的方式作用于管道外壁面，即將翅片管外的傳熱系數折算為豎直光管面積所對應的傳熱系數，從而將翅片管外壁處理為無翅片的圓光管，圓光管的直徑與翅片管外直徑相同。為了實現沸騰與固體基質的傳熱耦合，保留了管壁的厚度。

由于豎直光管具有對稱性，因此本文建立的幾何模型為二維模型，見圖1。LNG從管道的底部進入，受熱氣化后從頂部流出。幾何坐標軸的原點O在管道底部入口的中心處，x軸沿管道的徑向方向，y軸與管道中心軸線重合，方向與LNG的流動方向一致。管道的內直徑為0.02 m，外直徑為0.026 m，管道的長度為5 m。

圖1 幾何模型

4 網格劃分

本文將計算域分為LNG流體域和鋁合金管道固體域，計算域網格劃分見圖2。首先在ICEM 19.0軟件中建立幾何模型，并創建塊文件，通過兩次切分，在幾何模型與塊之間形成映射，最終生成結構化四邊形網格。對鋁合金管內外壁面附近法線方向上的LNG流體區域進行邊界層加密。經過網格無關性檢驗后確定，在管內側采用BiGeometric分布規律，節點數為131，第1層的高度為1.0×10-4m，增比Ratio=1.02；鋁合金管壁內的網格采用Uniform節點平均分布，節點數為5，每1層的高度為7.5×10-4m。沿鋁合金管高度方向，網格高度相同，每層高度為2.5×10-3m。最終所得的網格總數為1 379 655。

圖2 計算域網格劃分

5 各項設置

采用ANSYS Fluent 19.0軟件進行計算，求解器設定為：Double precision、Parallel processing(12 processes)、Pressure-based、Absolute、Unsteady、2D。

模型設置：求解質量、動量以及能量控制方程，多相流模型選擇VOF模型。湍流模型為標準k-ε湍流模型、近壁處為標準壁面函數。沸騰相變模型選擇S. Hardt源項模型，通過用戶自定義(UDF)的方式將S. Hardt源項模型加入控制方程中。選用分段線性構造界面(Piecewise Linear Interface Construction)的幾何重構法捕捉氣液兩相界面。

材料設置：LNG和天然氣的物性數據采用文獻[10]中的方法，分別對氣液兩相的密度、比熱容、動力黏度和熱導率進行計算，同樣進行多項式擬合，在Fluent中創建LNG和天然氣兩種物質，其中LNG中各烷烴組分的摩爾分數分別為，甲烷88%，乙烷8%，丙烷4%。固體域鋁合金的物性參數采用默認設置，其密度為2 623.3 kg/m3, 比熱容為983 J/(kg·K),熱導率為227.95 W/(m·K)。

操作條件設置：打開重力場，沿管道軸線方向，即y軸方向，設置重力加速度為-9.81 m/s2。

邊界條件設置：對于鋁合金管內，入口處設置為速度入口(velocity inlet)，LNG的速度在0.05～0.20 m/s范圍內，質量分數為1，入口溫度為111～126 K。翅片管出口采用壓力出口邊界(pressure outlet)，出口處壓力設置為0.5 MPa，液相的回流參數設置為0。鋁合金管外壁設置為wall邊界條件，其表面傳熱系數變化范圍為100～300 W·m-2·K-1，來流環境溫度設定為300 K。管內壁面會被Fluent軟件自動識別流固耦合傳熱，在讀入網格時通常為流固界面生成一個對應的shadow面，并將這兩個面在傳熱問題上確定為耦合(coupled)關系。

計算工況的設計見表1。

表1 計算工況

續表1

6 求解

各方程的離散格式：Pressure-Velocity Coupling Scheme采用PISO，Gradient選擇Least Squares Cell Based，Pressure選擇PRESTO，其余方程均采用Second Order Upwind格式進行離散。

欠松弛因子設定：為保證計算過程的穩定性，在計算前期先調低所有欠松弛因子，待計算穩定后再將其設定default，即Pressure、Density、Body Forces、Momentum、Vaporization Mass、Slip Velocity、Volume Fraction、Turbulent Kinetic Energy、Turbulent Dissipation Rate、Turbulent Viscosity、Energy的取值分別為0.3、1、1、0.7、1、0.1、0.5、0.8、0.8、1、0.91。

殘差設定：Energy殘差為1×10-6，其余為1×10-3。

初始化：先選用Standard Initialization，Relative to Cell Zone方式，Gauge Pressure、X velocity、Y velocity、Z velocity均為0，Turbulent Kinetic Energy為8.437 5×10-5，Turbulent Dissipation Rate為1.657 8×10-3。由于本文采用瞬態模擬計算，需要確保初始化的結果盡可能貼近實際物理現象，因此，通過Patch功能分別設置管內計算域的溫度為300 K，壓力為5×105Pa，天然氣的體積分數為1；固體鋁合金管的溫度為300 K。

本文研究的內容是從LNG進入鋁合金豎直翅片管開始氣化到管內的氣化狀態達到穩定的過程，因此采用瞬態模擬計算，時間步長采用自適應的方式，將步長控制在1×10-5～1×10-4s。每個時間步長內最大迭代次數為200。Fluent實時監測管道出口天然氣的體積分數變化，當體積分數變化達到穩定時，即可認為此時管內流動沸騰達到穩定狀態。

7 計算結果及分析

7.1 管內流場特性

本文定義工況0：LNG的入口流速為0.15 m/s，入口溫度為126 K，管外表面傳熱系數為300 W/(m2·K)。圖3～6展示了工況0條件下，豎直翅片管管內LNG氣化過程達到穩定狀態后，不同高度區域內氣相體積分數云圖和速度矢量分布(軟件截圖)，其中，色標右側的標值為氣相體積分數，速度矢量用黑色箭頭線表示。沿管道豎直向上(y軸正方向)，管內不同位置所表現出的流場特性有明顯差異。根據天然氣氣泡的大小、位置、形狀以及密集程度，可以劃分不同的兩相流型，例如泡狀流、彈狀流以及攪拌狀流等。從圖3～6可以看出，當LNG剛進入管道后，管內主流區域的速度場分布較為均勻，但是由于受到熱壁面的加熱，管壁附近有零星的氣泡產生，氣泡處的流動出現輕微紊亂。隨著LNG向上流動，受到管壁持續加熱，管內氣泡數量增加，氣泡直徑增大。在氣泡存在的區域，氣相的速度比液相的速度大，氣液兩相之間有明顯的速度差，液相受到氣相的擾動更強。在圖5和圖6中，一些兩相流動區域出現了局部回流的現象。

圖3 工況0穩定狀態時管內y=[0.015 m, 0.047 m] 區域內氣相體積分數云圖和速度矢量分布(軟件截圖)

圖4 工況0穩定狀態時管內y=[0.145 m, 0.177 m] 區域內氣相體積分數云圖和速度矢量分布(軟件截圖)

圖5 工況0穩定狀態時管內y=[0.385 m, 0.417 m] 區域內氣相體積分數云圖和速度矢量分布(軟件截圖)

圖6 工況0穩定狀態時管內y=[0.725 m, 0.757 m] 區域內氣相體積分數云圖和速度矢量分布(軟件截圖)

7.2 氣泡行為

氣泡的生長過程受到多個力的共同作用，主要有體積力、表面張力、來自氣液界面的黏性應力、氣泡內的蒸氣壓力以及熱壁面對氣泡的反作用力。起初，氣化核心的半徑較小，盡管熱壁面的液膜溫度較高，但是氣泡內的壓力大于液相壓力，其對應的飽和溫度比較高。此時液膜與氣泡之間的溫差較小，將其視為等溫生長過程，在這個階段氣泡生長的推動力主要是氣泡內部的蒸氣壓力與液相作用在氣液界面上的應力之和減去氣液之間表面張力。當氣泡直徑逐漸增大時，其受到的體積力迅速增大，氣泡內的蒸氣壓力逐漸減小，溫度隨之降低，液膜與氣泡之間的溫差逐漸變大，此后氣泡的生長過程便受液膜與氣泡之間的傳熱過程主導。

圖7～9為工況0條件下，豎直翅片管管內高度為y=[0.075 m, 0.210 m]區域內熱壁面上氣泡間的聚合現象。其中，色標右側的標值為氣相體積分數。為了便于觀察，使用白色方框框選圖中選定的氣泡，使用粉紅色箭頭指示氣泡行為的變化歷程。

圖7 工況0條件下熱壁面上氣泡間的聚合行為 (軟件截圖)

隨著氣泡本身的不斷生長，當氣泡直徑增大到一定范圍時，氣泡之間就會出現聚合現象，將這種沿管長方向的聚合稱為縱向聚合。從圖7中可以看出，氣泡的聚合行為最先出現在熱壁面附近，上游氣化核心處產生的氣泡多數并未直接脫離，而是受到浮力和液相的拖曳力，沿著壁面向上滑移，同時由于持續受熱，其直徑不斷增大，相鄰兩個氣泡之間的距離減小，從而逐漸發生融合。

從圖8可以看到，在時間為0.65 s時，兩個氣泡間的液膜已經開始逐漸融合，最終在時間為0.67 s時，兩個氣泡在靠近管道中心軸的主流區實現了完全融合，成為體積更大的氣泡。在氣泡脫離熱壁面進入主流區以后，由于浮力和拖曳力的作用，氣泡上升速度不斷加快，但是對于處在上下相鄰位置的兩個氣泡而言，下方的氣泡會受到上方氣泡的尾流作用，其上升加速度大于上方氣泡的加速度，使得兩個氣泡之間的距離不斷減小。

圖8 工況0條件下主流區氣泡縱向聚合(軟件截圖)

在氣泡重點聚合區，其數量逐漸增多導致間距變小，在主流區出現了氣泡間的橫向聚合，見圖9。在時間為1.23 s時，圖9中標出的兩個氣泡基本處于同一水平高度，且兩者之間尚存在較厚的液膜，但是在時間為1.26 s時，兩個氣泡之間的距離縮短，逐漸發生聚合行為。值得注意的是，與熱壁面處氣泡的上升速度相比，圖9所示的主流區的氣泡上升速度較小，原因是熱壁面附近氣相體積分數較大，導致混合相流速高于主流區。

圖9 工況0條件下主流區氣泡橫向聚合(軟件截圖)

在主流區，氣泡失去熱壁面的橫向束縛，常以球狀或者橢球狀向上運動。上升氣泡的形態和大小與流體本身的物性有較大的關系，主要受到表面張力和上升浮力的影響。奧托斯數[11]正是表征氣泡上升過程中表面張力和浮力的相對影響，其計算公式見式(1)。

(1)

式中Eo——奧托斯數

de——與氣泡體積相等的圓球的直徑，m

g——重力加速度，m/s2

ρl——液相密度，kg/m3

ρg——氣相密度，kg/m3

σ——表面張力系數，N/m

當離散相為空氣，連續相為水時，Eo臨界值為40，Eo低于臨界值時氣泡主要呈現球狀或者橢球狀，高于臨界值時則由帽狀轉變為不規則的形狀[11]。由于Eo臨界值是混合相物性參數的函數，因此，本文根據空氣-水混合相的Eo臨界值，以及LNG-天然氣混合相的密度和表面張力系數，推算出LNG-天然氣混合相的Eo臨界值為7 443。天然氣氣泡的體積在上升過程中不斷增大，其等效直徑增大后，從式(1)中可以看出，等式右側分子部分所代表的浮力將逐漸增大，并開始占據主導地位。當Eo高于7 443時，氣泡形狀隨之發生相應變化，不再呈現規則的球狀或橢球狀。從圖9可以看出，氣泡的體積越大，其形狀越扭曲。此外，從圖4可以發現，氣泡在靠近熱壁面側和主流側的速度相差較大，氣泡緊貼熱壁面側的速度大于主流區的速度，因此，形成了拖曳的效果，導致圖4和圖9中熱壁面上的大體積氣泡多呈現“半個箭頭”形狀。

7.3 氣液兩相流型

本文模擬結果表明，LNG在管內的流動沸騰過程中出現的流型有泡狀流、彈狀流、攪拌流、霧狀流，并未出現環狀流。在對流型的劃分中，本文部分借鑒了文獻[5]的劃分方法，具體的劃分方法見表2。與文獻[5]不同的是，將氣相體積分數為0～0.18階段分為了單液相流和泡狀流兩種流型。對于單液相流，主流溫度盡管未達到飽和溫度，但由于存在過冷沸騰的原因，依然會有氣泡產生，極少量的小氣泡脫離壁面后，受到過冷液體的冷卻而快速湮滅，因此，將其作為兩相流區的特殊階段。

表2 LNG豎直翅片管管內氣液兩相流型劃分

從LNG進入管道一直到管內完成氣化，即氣相天然氣的體積分數從0達到1，這個區域統稱為兩相流區域。各流型在管內分布長度的統計依據是流型的特征和對應的體積分數。各流型在管內分布長度與兩相流區域總長度的比值，稱為兩相流中各流型占比。在上述眾多流型中，攪拌流所對應的流動沸騰傳熱系數較大[5]，故而對不同LNG入口流速、LNG入口溫度以及管外表面傳熱系數等條件下的管內流動沸騰進行模擬研究，目的在于擴大攪拌流區域，增大流動沸騰傳熱系數，提高空溫式氣化器的利用效率。

7.3.1 LNG入口溫度對兩相流型分布的影響

LNG一般經泵加壓后進入空溫式氣化器，進入空溫式氣化器的LNG并不會即刻蒸發沸騰。若進口LNG的溫度太低，將使得氣化器管長增加，導致設備體積過大。圖10為豎直翅片管管內流動沸騰達到穩定狀態后LNG入口溫度對兩相流型分布的影響，其計算條件是，LNG的入口流速為0.15 m/s，管外表面傳熱系數為300 W/(m2·K)，入口溫度為111～126 K。

圖10 LNG入口溫度對兩相流型分布的影響

從圖10可以看出，在不同LNG入口溫度的條件下，管內兩相流中各流型占比不同。當空溫式氣化器的入口處溫度為111 K時，管內單液相流區域較長，占氣液兩相流區域的31.4%；其次是泡狀流區域，其占比為18.1%。單液相流和泡狀流的區域均隨著LNG入口溫度的升高而減小，其中單液相流區域的變化幅度最大，當溫度升高到126 K時，單液相流區域的占比下降到8.2%。這是因為當入口LNG溫度較低時，低溫的LNG需要在管內流動更長的距離、吸收更多的熱量，才能進入泡狀流區域。對于泡狀流區域而言，在管壁和管道中軸線之間存在較大的溫度梯度，靠近熱壁面處氣泡密度較大，中軸線附近的氣泡密度小，由此使得泡狀流區域受到一定程度的拉長。隨著LNG入口溫度的增高，攪拌流區域的占比獲得了較大的提高。當入口溫度為126 K時，攪拌流和霧狀流區域占比的總和達到了70.6%，表明LNG的入口溫度越高，管內的流動沸騰換熱系數越大，換熱效果越好。

7.3.2 LNG入口流速對兩相流型分布的影響

隨著空溫式氣化器負荷的變化，入口處的LNG流速會隨之發生變化。當LNG的入口流速為0.05 m/s，對應的天然氣出口流量(折合成標準狀態)為34 m3/h。圖11為豎直翅片管管內流動沸騰達到穩定狀態后LNG入口流速對兩相流型分布的影響，其計算條件是，LNG的入口流速為0.05～0.20 m/s，入口溫度為126 K，管外表面傳熱系數為300 W/(m2·K)。

圖11 LNG入口流速對兩相流型分布的影響

由圖11可見，隨著LNG入口流速的增大，單液相流、攪拌流、霧狀流在兩相流中的占比逐漸增大；泡狀流與彈狀流區域則隨著入口流速的增大而減小，且彈狀流區域減小的幅度最大。當入口流速為0.05 m/s時，單液相流段較小，在兩相流中的占比為3.8%。這是由于較小的入口流速，意味著管內的質量流量較小，主流區的LNG經過較短的管程距離便達到了泡狀流階段。當LNG的入口流速從0.05 m/s增大到0.20 m/s時，泡狀流區域和彈狀流區域的占比(絕對值)分別減少了9.5%和22.9%。這是因為隨著入口流速的增加，氣液兩相之間的相對速度不斷增大，氣液界面之間的黏性應力增加，使得氣泡尺寸被拉長，且管內氣液兩相的流速增大后，兩相流的湍流強度增大，氣泡之間更容易擠壓碰撞，由此導致管內的泡狀流和彈狀流區域不斷縮小。與霧狀流相比，攪拌流的占比僅從28.8%增加到33.7%，增加幅度較小。原因主要為隨著氣相體積分數的增大，氣相的流速劇烈升高，較高的氣流速度使得霧狀流區域中的液相體積分數升高，增大了霧狀流區域在管內的分布長度。

7.3.3 管外表面傳熱系數對兩相流型分布的影響

由于環境空氣變化以及管外結霜的原因，管外的表面傳熱系數會有明顯的變化。圖12為豎直翅片管管內流動沸騰達到穩定狀態后管外表面傳熱系數對兩相流型分布的影響，其計算條件是，LNG的入口流速為0.15 m/s，入口溫度為126 K，管外表面傳熱系數為100～300 W/(m2·K)。

圖12 管外表面傳熱系數對兩相流型分布的影響

由圖12可見，隨著管外表面傳熱系數的提高，單液相流和泡狀流在兩相流中的占比快速下降，攪拌流和霧狀流的占比則是迅速上升，而彈狀流占比的變化較小。當管外的表面傳熱系數為100 W/(m2·K)時，單液相流在氣液兩相流中的占比為30.2%，泡狀流的占比為33.4%，遠高于彈狀流和攪拌流的占比。這是因為，管外表面傳熱系數較低時，在管內壁處的熱流密度較小，使得液相中氣泡的直徑和密度均比較小，氣泡之間融合的概率就比較低，氣泡之間基本呈現彼此孤立的狀態。當管外表面傳熱系數從100 W/(m2·K)增大到300 W/(m2·K)時，熱壁面上氣泡受到較高的熱流密度加熱，無論是氣泡直徑還是脫離熱壁面的頻率都會增大，使得攪拌流的占比從8.6%增加到32.3%。

鑒于攪拌流所對應的流動沸騰傳熱系數較大，增大攪拌流區域，有利于提高管內的流動沸騰傳熱系數，進而提高空溫式氣化器的利用效率。對比圖10～12中的數據可知，當入口溫度從111 K增大到126 K時，攪拌流區域占比增大了14.3%；入口速度從0.05 m/s增大到0.20 m/s時，攪拌流區域占比僅增大了4.8%；而改善管外的表面傳熱系數使其從100 W/(m2·K)增大到300 W/(m2·K)時，則可使攪拌流區域占比增大23.7%。因此，與增大LNG入口溫度、LNG入口流速相比，通過改善管外的表面傳熱強度更能夠提高空溫式氣化器的換熱效率。

8 結論

為了揭示LNG空溫式氣化器管內流動沸騰特征，對沸騰過程中氣泡行為和兩相流型的分布和轉換進行了研究。研究對象為豎直鋁合金翅片管，并將其等效簡化為豎直光管，其高度為5 m，管內直徑為0.02 m，管外直徑為0.026 m。管內LNG的入口流速為0.05～0.2 m/s，入口溫度為111～126 K，管外表面傳熱系數變化范圍為100～300 W/(m2·K)，環境溫度為300 K，環境壓力為101.325 kPa。采用ICEM 19.0軟件建立豎直光管的二維模型并劃分四邊形結構化網格，利用Fluent 19.0軟件對管內液化天然氣或天然氣計算域和管壁固體計算域進行傳熱模擬，研究不同LNG入口流速、LNG入口溫度以及管外表面傳熱系數等條件對沸騰過程中氣泡行為和兩相流型分布和轉換的影響。采用流體體積法對氣液兩相界面追蹤捕捉，利用用戶自定義的方式將表征沸騰過程的S. Hardt模型引入計算流程中。研究結果表明：

① 氣泡的聚合及沿熱壁面的滑移是在表面張力、浮力以及氣液之間剪切應力作用下的重要運動形式。

② LNG入口溫度增加可以減小單液相流在兩相流中的比例，同時增大攪拌流和霧狀流區域；LNG入口流速增大時，兩相流中泡狀流和彈狀流區域快速減少，攪拌流區域小幅增加；管外表面傳熱系數增大，可顯著提高攪拌流的占比。

③ 與增大LNG入口溫度、LNG入口流速相比，通過改善管外的表面傳熱強度更能夠提高空溫式氣化器的換熱效率。