姚寶龍，唐建峰，王國瑞，王劍琨，于 笑，許 濤

（1. 中國石油大學(xué)（華東） 儲運與建筑工程學(xué)院，山東 青島 266580；2. 中國石化天津液化天然氣有限責(zé)任公司，天津 300457；3. 中國石化青島液化天然氣有限責(zé)任公司，山東 青島 266580）

液化天然氣接收站作為LNG產(chǎn)業(yè)鏈中的重要環(huán)節(jié)，近年來在我國得到大力發(fā)展。截止2020年，LNG接收站數(shù)量已達(dá)22座[1]。接收站儲罐內(nèi)的LNG使用時經(jīng)過氣化器氣化后外輸。中間介質(zhì)式氣化器（Intermediate Fluid Vaporizer，IFV）由于其沒有結(jié)冰問題，對海水質(zhì)量要求低、占地面積小等優(yōu)點，在沿海LNG接收站及海上浮式氣化場景中有廣闊的應(yīng)用前景。

IFV結(jié)構(gòu)簡圖如圖1所示，根據(jù)換熱過程可將IFV分為蒸發(fā)區(qū)（E1）、凝結(jié)區(qū)（E2）和調(diào)溫區(qū)（E3）。IFV熱源介質(zhì)為海水，中間介質(zhì)為丙烷。LNG在E2段吸收管外丙烷的熱量后氣化成為NG，最后經(jīng)E3升溫后外輸，丙烷在E2段殼側(cè)釋放熱量凝結(jié)成液，在重力作用下進(jìn)入E1段與管內(nèi)海水換熱蒸發(fā)后返回至E1段，形成丙烷的“蒸發(fā)-冷凝”內(nèi)部循環(huán)過程。近些年隨著氣化器國產(chǎn)化進(jìn)程的加速，國內(nèi)已有部分廠家掌握了IFV的設(shè)計及制造，但現(xiàn)場應(yīng)用經(jīng)驗不足。國內(nèi)學(xué)者對IFV的研究更多集中在材料、結(jié)構(gòu)及運行操作等方面，對其傳熱機(jī)理研究較少。劉豐等[2]開展了IFV在海上浮動工況下的適用性研究，從結(jié)構(gòu)、材料等方面進(jìn)行了論證。陳雙雙等[3]對比了不同IFV的布置方式，采用MATLAB軟件分析了整體式IFV與分體式IFV的換熱性能。王博杰等[4]利用CFD數(shù)值模擬對LNG在IFV內(nèi)的流動傳熱特性進(jìn)行研究，得到了IFV內(nèi)部流動傳熱與溫度分布特性。IFV涉及多流體的不同換熱過程，其內(nèi)部復(fù)雜的換熱機(jī)理仍未明晰。其整個傳熱過程有兩個難點，一是中間介質(zhì)工作狀態(tài)的確定，二是凝結(jié)區(qū)管內(nèi)LNG的跨臨界流動傳熱過程。IFV在使用過程中工作壓力一般在6~10 MPa范圍內(nèi)，LNG在與管外丙烷換熱后由亞臨界態(tài)轉(zhuǎn)變至超臨界態(tài)。目前，針對流體跨臨界流動傳熱的研究主要集中在H2O和CO2等常規(guī)工業(yè)流體中[5-8]，針對特定組分LNG的研究較少，同時文獻(xiàn)中研究大多針對水平圓管和豎直管內(nèi)流體的跨臨界流動傳熱過程[9-11]，無法對IFV內(nèi)U型管工況提供參考。

圖1 IFV結(jié)構(gòu)Fig. 1 Structure of IFV

本文對LNG在IFV中的跨臨界流動換熱過程進(jìn)行數(shù)值模擬，對超臨界狀態(tài)下LNG熱物性進(jìn)行計算，探究管內(nèi)LNG在跨臨界升溫過程中的傳熱特性與流動特性，分析操作參數(shù)對IFV運行情況的影響。

1 超臨界態(tài)下LNG熱物性計算

在進(jìn)行數(shù)值計算時需要準(zhǔn)確獲得實際LNG組分下各物性隨溫度的變化規(guī)律。LNG通常含有含量（物質(zhì)的量分?jǐn)?shù)，下同）為90%以上的甲烷，此外還有少量的C2H4、C3H6、C5+及N2等。表1為國內(nèi)沿海某接收站貧LNG組分。根據(jù)表1中LNG組分，利用工質(zhì)物性計算軟件Refprop 9.1對所需工況下LNG的熱物性進(jìn)行計算，以9 MPa為例計算結(jié)果如圖2所示。由圖2(a)可知，LNG在跨臨界過程中各項熱物性隨溫度連續(xù)變化，氣液分界面消失，不存在明顯的相變現(xiàn)象。由定壓比熱容曲線可以看出，LNG在跨臨界升溫過程中定壓比熱容存在明顯的極值點，該點對應(yīng)的溫度稱為擬臨界點。由圖2(b)可知，壓力越高，定壓比熱容曲線越平緩且峰值越低。本文采用變物性的湍流模型對此過程進(jìn)行描述分析[12]，利用上述計算得到的LNG熱物性隨溫度變化的規(guī)律，通過FLUENT軟件自帶的線性差值（Piecewise-liner）方式輸入到計算過程中。

圖2 9 MPa時LNG熱物性隨溫度的變化(a)和壓力對LNG定壓比熱容的影響(b)Fig. 2 Variation of LNG thermophysical properties with temperature at 9 MPa (a) and effect of pressure on specific heat capacityof LNG (b)

表1 貧LNG組分Table 1 Lean LNG components

2 數(shù)值計算模型

2.1 物理模型

國內(nèi)沿海某LNG接收站使用的國產(chǎn)化IFV，具體尺寸參數(shù)及性能指標(biāo)見表2。其凝結(jié)區(qū)內(nèi)有多根LNG傳熱管，各傳熱管結(jié)構(gòu)相同，選取單根傳熱管作為研究的計算域。IFV冷凝區(qū)內(nèi)傳熱管規(guī)格如表3所示。利用ICEM對計算域進(jìn)行建模及結(jié)構(gòu)化網(wǎng)格劃分。對彎管部分進(jìn)行網(wǎng)格加密，以更好地捕捉彎管區(qū)域流體流動特性。網(wǎng)格劃分如圖3所示，其中Z軸起始點為管段入口處。

表2 IFV設(shè)計參數(shù)Table 2 IFV design parameters

表3 凝結(jié)區(qū)傳熱管規(guī)格Table 3 Specification of heat transfer tube in condensation area

圖3 彎管及其橫截面網(wǎng)格劃分Fig. 3 Bend and its cross section meshing

LNG在傳熱管內(nèi)的流動受到加熱壁面的約束，在近壁區(qū)會形成邊界層，F(xiàn)LUENT為湍流模型提供了不同的壁面函數(shù)來處理邊界層內(nèi)的流動與傳熱。本文采用標(biāo)準(zhǔn)壁面函數(shù)處理近壁面區(qū)的流動，標(biāo)準(zhǔn)壁面函數(shù)需要將第一層網(wǎng)格節(jié)點置于對數(shù)律層，引入無量綱化的壁面距離y+定義第一層網(wǎng)格高度，對數(shù)律層要求y+的范圍在30~300之間[13]，y+的定義式如式(1)。通過式(1)可估算第一層網(wǎng)格高度y。

式中，y為第一層網(wǎng)格高度，m；τw為壁面的切應(yīng)力，MPa；ν為流體的平均動力黏度，Pa·s；ρ為流體密度，kg/m3。

在符合上述網(wǎng)格劃分要求下，劃分了4組不同網(wǎng)格數(shù)量的換熱管模型并進(jìn)行網(wǎng)格獨立性檢驗，結(jié)果如圖4所示。由圖4可知，隨著網(wǎng)格數(shù)量的提升，計算結(jié)果逐漸趨于穩(wěn)定，考慮到計算時間成本，選擇網(wǎng)格數(shù)量為1264784的網(wǎng)格進(jìn)行后續(xù)計算。

圖4 網(wǎng)格獨立性驗證Fig. 4 Grid independence verification

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

為使得數(shù)值計算結(jié)果更精確，現(xiàn)對LNG在IFV換熱管中的實際氣化傳熱過程進(jìn)行如下假設(shè)：（1）將LNG平均分配到每根傳熱管，且LNG經(jīng)過管線到達(dá)氣化器入口時處于充分發(fā)展段；（2）當(dāng)氣化器殼側(cè)工況達(dá)到穩(wěn)定后，管外壁溫度與丙烷氣相主體溫度相差不大，忽略該部分溫差，設(shè)置模型中傳熱管的外壁溫度為恒定值；（3）管程流體的壓降（0.04 MPa）在數(shù)量級上遠(yuǎn)遠(yuǎn)小于工作壓力，可忽略不計?；谏鲜黾僭O(shè)，LNG在傳熱管內(nèi)的流動傳熱過程滿足連續(xù)性方程、動量守恒方程和能量守恒方程，如式(2)~式(5)。

連續(xù)性方程：

式中，為Hamilton算子；為速度向量，m/s。

對于不可壓縮牛頓流體，N-S方程如下：

式中，E為熱力學(xué)能，J；keff為流體的有效導(dǎo)熱系數(shù)，W/(m·K)；方程右邊前兩項分別是由導(dǎo)熱引起的能量傳遞和粘性耗散引起的能量傳遞；τeff為粘性應(yīng)力張量；Sh為體積熱源項，W/m3。

整個傳熱過程還應(yīng)包括換熱管壁中的熱傳導(dǎo)：

式中，xi代表i的方向，i= 1, 2, 3；λ為固體壁面的熱導(dǎo)率，W/(m·K)。

IFV內(nèi)部LNG傳熱管是U型傳熱管，而RNGk-ε模型在模擬強(qiáng)旋流、彎曲壁面或彎曲流線流動時表現(xiàn)出較好的吻合性。本研究選取RNGk-ε模型，壁面函數(shù)選取Standard Wall Functions，該模型的輸運方程如式(6)~式(7)。

湍流動能方程：

式中，k為流體的湍動能，m2/s2；αk為k方程的湍流普朗特數(shù)；ui為速度，m/s；μeff為有效動力黏度，Pa·s；Gk為由平均速度梯度引起的湍流動能，m2/s2；Gb為由浮力引起的湍流動能，m2/s2；ε為湍流耗散率，m2/s3；YM為可壓縮湍流流動脈動膨脹對總的耗散率的影響；Sk為k方程中的源項。

湍流耗散率方程：

式中，αε為ε方程的湍流普朗特數(shù)；G1ε和G2ε均為模型常量；G3ε為決定ε方程受浮力影響程度的量；Sε為ε方程中的源項。

本文采用商業(yè)CFD軟件FLUENT1 9.0對數(shù)值過程進(jìn)行求解。求解方法采用基于壓力-速度耦合的SIMPLEC算法，壓力項采用Body Force Weighted，動量方程和能量方程采用QUICK格式進(jìn)行離散，時間步長設(shè)置為0.001 s，設(shè)置迭代步數(shù)為2000步。在計算過程監(jiān)測出口截面的平均溫度，當(dāng)出口截面平均溫度不再變化后，認(rèn)為計算收斂。

2.3 模型準(zhǔn)確性驗證

以某LNG接收站夏季某日現(xiàn)場運行數(shù)據(jù)為基礎(chǔ)，對數(shù)值模型的準(zhǔn)確性進(jìn)行驗證，該日實際運行數(shù)據(jù)如表4所示，驗證結(jié)果如表5所示。通過對比模擬數(shù)值與實際監(jiān)測數(shù)值可知，模擬計算得到的NG出口溫度與實際運行數(shù)值誤差為0.1%，出口壓力與實際數(shù)值誤差為4.9%，考慮到模型中海水溫度、LNG入口溫度等與實際運行存在一定的偏差，因此誤差在合理范圍之內(nèi)，可利用此數(shù)值模擬方法進(jìn)行進(jìn)一步的研究。

表4 某LNG接收站氣化單元IFV某日實際運行數(shù)據(jù)Table 4 Actual operation data of IFV of gasification unit in a LNG terminal on a certain day

表5 IFV模擬數(shù)值與實際數(shù)值對比Table 5 Comparison between IFV simulation value and actual value

3 結(jié)果與討論

3.1 管內(nèi)溫度場及速度場分析

以表4中的工況為例，對LNG在管內(nèi)跨臨界流動換熱的瞬態(tài)數(shù)值模擬結(jié)果進(jìn)行溫度場分析和速度場分析，結(jié)果分別如圖5和圖6。由圖5可知，管程的前半部分LNG處于液態(tài)，與管壁之間溫差較大，此時流體升溫較快，在溫度跨過擬臨界點后，流體的熱導(dǎo)率迅速降低，且流體溫度與壁面溫度溫差減小，升溫變緩。由圖6可知，在入口速度為0.48 m/s的工況下，沿管程流體速度不斷增大，出口處截面平均速度達(dá)到10.35 m/s，出現(xiàn)了明顯的流動加速效應(yīng)。管內(nèi)LNG沿程溫度及密度變化見圖7。由圖7可知，LNG在管內(nèi)吸熱后密度迅速降低，由424 kg/m3下降至77.93 kg/m3，過程中LNG不斷膨脹產(chǎn)生軸向推動力使得流體不斷加速。

圖5 管內(nèi)LNG沿程溫度場分布Fig. 5 Distribution of LNG temperature field along pipe

圖6 管內(nèi)LNG沿程速度場分布Fig. 6 Distribution of LNG velocity field along pipe

圖7 管內(nèi)LNG溫度及密度沿程變化Fig. 7 Variation of LNG temperature and density along pipe

3.2 擬臨界點傳熱特性分析

由圖2(b)可知，LNG在不同壓力下擬臨界溫度不同，查詢制冷劑熱物性計算軟件可知，該組分LNG在入口壓力為9.26 MPa下的擬臨界溫度為218.400 K。圖8給出了管內(nèi)LNG沿軸向2~8 m各處截面熱物性參數(shù)的分布云圖。從管內(nèi)LNG沿程的溫度分布與密度分布可以看出，流體在進(jìn)入換熱管后物性變化較為劇烈，在此區(qū)域內(nèi)LNG經(jīng)歷了跨臨界的狀態(tài)轉(zhuǎn)變過程。在Z= 2 m處，由圖8(a)可知，圓管頂部小部分的LNG較早達(dá)到擬臨界溫度，由于超臨界壓力下LNG的熱物性對溫度有較強(qiáng)的依賴性，由圖8(b)可知，在Z= 2 m處，靠近圓管上頂部LNG密度下降較快，截面上密度梯度較大，由于受到重力的作用，圓管各截面的密度并不是呈現(xiàn)“同心圓”式的分布，而是出現(xiàn)了明顯的分層現(xiàn)象。隨著LNG向前流動，當(dāng)管程超過6 m后，LNG主體溫度均超過擬臨界溫度，截面上密度變化開始趨緩。圖8(c)為不同位置截面上定壓比熱容變化云圖，由圖8(c)可知，定壓比熱容的峰值首先出現(xiàn)在橫截面上部，隨著LNG主體溫度達(dá)到擬臨界溫度，截面上定壓比熱容的峰值不斷向截面中心及下部移動。LNG進(jìn)入傳熱管后靠近壁面的部分被加熱而密度降低，受到重力的作用，LNG在截面上開始重新分布。

圖8 沿軸截面管內(nèi)LNG熱物性分布云圖Fig. 8 LNG thermophysical property distribution cloud chart on pipe section

圖9為距離LNG入口2 m處截面流線輪廓。從圖9可知，靠近管內(nèi)壁，較輕的流體受到下部較重流體的擠壓后沿管壁被推至上部，在圓管截面上形成徑向?qū)ΨQ的“二次環(huán)流”，使得管內(nèi)流體沿軸向加速的過程中，截面上物性也在重新分布。

圖9 距入口2 m處截面流線輪廓Fig. 9 Section streamline at 2 m from entrance

3.3 彎管處流動換熱分析

LNG在彎管處受到壁面的阻礙，流動方向發(fā)生突變，在浮升力和離心力的共同作用下，彎管處的流動換熱情況與直管段有著較大的區(qū)別，建立如圖10所示的3個觀測面（1-1'、2-2'和3-3'），分析LNG在彎管處的流動情況。

圖10 彎管處觀測面Fig. 10 Observation surface at bend

彎管處3個觀測面LNG溫度分布如圖11所示。由圖11(a)可知，相比于彎管外側(cè)區(qū)域，靠近彎管內(nèi)側(cè)區(qū)域流體溫度更高。由圖11(b)可知，進(jìn)入彎管前，由于物性的重新分布，較重的LNG在管內(nèi)下部流動，進(jìn)入彎管后，受到管壁的阻礙以及離心力的作用，較重的LNG沿彎管外側(cè)繼續(xù)流動。在截面內(nèi)，Y方向的分速度逐漸增大，外側(cè)較重的LNG受到管壁和上升部分LNG的擠壓沿管壁向內(nèi)側(cè)流動，在截面上形成了兩個對稱的“迪恩渦”[14]。當(dāng)LNG運動至2-2'截面時，Z方向的分速度逐漸趨于0，X方向和Y方向的分速度達(dá)到最大值，此時在橫截面上出現(xiàn)了4個對稱的“迪恩渦”，截面流場被完全重構(gòu)，流體邊界層厚度被進(jìn)一步削弱，顯著加強(qiáng)了彎管區(qū)域的局部對流換熱系數(shù)。

圖11 彎管觀測面處LNG溫度分布及流線輪廓Fig. 11 LNG temperature distribution and streamline at observation surface of bend

3.4 入口流量對傳熱的影響

實際生產(chǎn)運行中，由于受到下游用戶季節(jié)性需求不同的影響，單臺IFV流量波動范圍較大。選定單臺IFV，以外輸量50~200 t/h作為研究區(qū)間，根據(jù)IFV凝結(jié)區(qū)LNG換熱管尺寸及數(shù)量，換算為單管質(zhì)量流量，如表5所示。針對LNG不同入口質(zhì)量流量進(jìn)行模擬計算，分析不同入口流量對LNG流動傳熱的影響。

表5 單臺IFV氣化量對應(yīng)的入口質(zhì)量流量Table 5 Inlet mass flow corresponding to gasification volume of single IFV

LNG入口壓力為9.26 MPa，丙烷飽和溫度為288.15 K時，入口質(zhì)量流量對管內(nèi)LNG對流換熱的影響如圖12所示。由圖12(a)可知，同一截面LNG溫度隨著入口質(zhì)量流量的增加而減小，這是因為隨著質(zhì)量流量的增加，相同質(zhì)量的LNG經(jīng)過截面的時間減少，使得LNG與管壁的熱量交換減少。由圖12(b)可知，同一質(zhì)量流量下，局部對流傳熱系數(shù)呈現(xiàn)先增大后減小的趨勢，在彎管處達(dá)到最大值。這是因為在彎管前，LNG主體升至擬臨界點溫度，定壓比熱容達(dá)到最大值，使得局部傳熱系數(shù)出現(xiàn)峰值。質(zhì)量流量越大，局部傳熱系數(shù)的峰值越靠后，這是因為隨著質(zhì)量流量的增加，LNG流速增大，相同質(zhì)量的LNG吸熱達(dá)到擬臨界溫度點的換熱時間增加，導(dǎo)致峰值位置出現(xiàn)延后。從圖12(b)還可以看出，LNG入口質(zhì)量流量越大，局部對流換熱系數(shù)也越大，這是因為質(zhì)量流量變大，增加了管內(nèi)流體的湍流強(qiáng)度，有效抑制了管壁處流體邊界層增厚，降低了對流熱阻，從而增強(qiáng)了流體的對流傳熱能力。

圖12 質(zhì)量流量對管內(nèi)LNG沿程傳熱的影響Fig. 12 Effect of mass flow rate on heat transfer of LNG along pipe

3.5 丙烷飽和溫度對傳熱的影響

在凝結(jié)區(qū)中，氣態(tài)丙烷在管外凝結(jié)成液，其飽和溫度的大小直接關(guān)系到LNG在傳熱管內(nèi)的流動換熱情況，進(jìn)而影響到IFV的氣化能力。在LNG入口溫度為125 K、入口質(zhì)量流量為0.042 kg/s下，選 擇268 K、278 K、288 K和298 K 4種丙烷飽和溫度工況進(jìn)行模擬計算，結(jié)果如圖13所示。

圖13 丙烷飽和溫度對管內(nèi)LNG沿程傳熱的影響Fig. 13 Effect of propane saturation temperature on heat transfer of LNG along pipe

由圖13(a)可知，丙烷飽和溫度越高，相同截面處流體平均溫度越高，這在管程的后半段表現(xiàn)更為明顯。這是由于LNG在擬臨界溫度前主要是液態(tài)，物性對溫度的變化并不敏感，流體升溫較快，到達(dá)擬臨界溫度點之后LNG定壓比熱容隨著溫度的升高迅速降低，同時較高的丙烷氣相溫度會增大流體與管壁間的溫差，同一截面處的LNG平均溫度也隨之增加。由圖13(b)可知，丙烷的飽和狀態(tài)對管內(nèi)LNG的局部傳熱系數(shù)有較大的影響。丙烷的飽和溫度一定時，局部對流傳熱系數(shù)先升高，達(dá)到峰值后迅速降低，原因可能是當(dāng)LNG跨過擬臨界溫度后，發(fā)生了傳熱惡化的現(xiàn)象。結(jié)合圖 8(a)可知，頂部母線內(nèi)壁附近流體會比底部母線附近流體更早達(dá)到擬臨界溫度，由于LNG熱物性對溫度有強(qiáng)烈的依賴性，頂部母線附近的LNG密度迅速下降，體積比底部母線附近的流體膨脹的更快。此時由于密度差及傳質(zhì)的影響，LNG向底部區(qū)域流動，頂部沿軸向的質(zhì)量通量開始減小，同時由圖8(c)可以看出，隨著溫度的升高，頂部母線附近LNG熱導(dǎo)率迅速下降，導(dǎo)致傳熱開始惡化，且丙烷飽和溫度越高，這種現(xiàn)象越明顯。

4 結(jié)論

采用數(shù)值模擬研究了LNG在IFV U型管中跨臨界流動傳熱的過程，探究了不同入口質(zhì)量流量、不同丙烷飽和溫度對LNG流動傳熱的影響，分析了LNG在擬臨界點附近的傳熱機(jī)理以及在彎管處的流動傳熱特性。主要結(jié)論如下：

（1）升溫過程中，LNG在超臨界壓力下沿軸向物性變化劇烈，出現(xiàn)了明顯的“流動加速”效應(yīng)。在截面上，由于溫度分布不均，密度會出現(xiàn)明顯的分層現(xiàn)象，從而形成對稱的徑向“二次環(huán)流”，加強(qiáng)了流體的換熱能力。

（2）在彎管處，受到離心力的作用及管壁阻礙，在截面上形成了對稱的“迪恩渦”，不同觀測面上流線分布差異較大，當(dāng)LNG運動至彎管的中心截面時，徑向流動進(jìn)一步加強(qiáng)，在截面上形成了4個對稱的“迪恩渦”，內(nèi)部流場被完全重構(gòu)，流體邊界層被進(jìn)一步削弱，顯著加強(qiáng)了彎管處的局部換熱系數(shù)。

（3）局部對流傳熱系數(shù)沿管程呈現(xiàn)先增大后減小的趨勢，在擬臨界點附近達(dá)到最大值。隨著入口質(zhì)量流量的增加，局部對流傳熱系數(shù)的峰值明顯增大，峰值出現(xiàn)的位置向后移動。同一截面處LNG溫度呈現(xiàn)減小的趨勢，但對出口NG溫度影響不大。

（4）隨著丙烷飽和溫度的升高，相同截面處流體平均溫度也越高，但對局部傳熱系數(shù)的影響相反，局部傳熱系數(shù)的峰值隨著丙烷飽和溫度的升高而降低。隨著丙烷飽和溫度的升高，U型管內(nèi)外介質(zhì)溫差增大，容易發(fā)生傳熱惡化的現(xiàn)象。