分相的多股流LNG繞管式換熱器動態模型

段鐘弟，任滔，丁國良，陳杰，浦暉，密曉光

（1上海交通大學制冷與低溫工程研究所，上海200240；2中海石油氣電集團技術研發中心，北京100028）

引 言

天然氣作為一種高效清潔的能源，近年來得到了廣泛的應用［1］。我國南海擁有豐富的待開發天然氣資源，已探明的儲量占我國天然氣總儲量的1／3［2］。開采海上天然氣的基礎裝置是浮式天然氣生產儲卸平臺 （LNG－FPSO）。LNG－FPSO的主低溫換熱器為繞管式換熱器，具有占地面積小、力學性能好和抗泄漏能力強等優點［3］。

繞管式換熱器由數以千計的管子按層排布螺旋纏繞在中心筒上，單個換熱器投資巨大，不能基于實驗進行設計，需要先進的設計方法［4］?；诜抡娴脑O計方法具有高效、低成本的優勢，能夠應用于繞管式換熱器設計［5］。LNG－FPSO是海上浮動平臺，會隨風浪在海上晃動，從而造成繞管式換熱器晃動，使得其內部流體的溫度、壓力隨時間發生波動［6］，因此需要能夠預測繞管式換熱器動態特性的模型。

因天然氣液化獨特的工藝和需求，繞管式換熱器的結構和工作具有如下特征，包括：①工質是由烴類和氮氣組成的混合冷劑；②工況范圍覆蓋包含亞臨界和超臨界的所有相區；③多股流體并行換熱；④單個換熱器內冷熱流體并發相變，同時存在蒸發和冷凝?；诩倕档亩喙闪鲹Q熱器動態模型［7－10］不能反映繞管式換熱器中的相變特征；基于分布參數的多股流動態模型［11－12］涉及大量微元的計算，計算耗時長［13］。因此，本研究采用具有良好精度和計算速度的移動邊界模型［14］來預測繞管式換熱器的動態特性?，F有的移動邊界模型可以對繞管式換熱器前兩項特征中的多組分工質、多相區模式進行計算，具體見表1；而對于后兩項特征中多股流并行換熱和冷熱流體并發相變的計算，還需要進一步的研究。

表1 已有研究中代表性的換熱器動態模型Table 1 Existing representative dynamic models of heat exchangers

多股流并行換熱中，一股流體同時與其他多股流體發生熱交換，使得各股冷熱流體之間的換熱關系復雜化［21］。同時，在移動邊界模型下流體的控制單元按照相區進行劃分，由于各個相區的長度不一，在空間分布上互相交錯，使得分屬不同流路的相區之間換熱關系復雜化。以上兩方面的原因導致移動邊界模型中的多股流并行換熱難以計算。針對特定多股流形式的建模方法［8－9］不能滿足繞管式換熱器各種不同流程布置的要求；采用基于集總參數的互連矩陣來描述流路換熱關系的方法［7］不適用于移動邊界模型下相區之間的換熱計算。因此，需要一種能夠描述繞管式換熱器多股流多相區的換熱關系，對任意的流路和相區分布形式進行計算的方法。

當換熱器中同時存在蒸發和冷凝相變過程時，管壁兩側的冷熱流體均存在相邊界。在動態過程中，兩側流體的相邊界同時發生移動，一方面會出現相邊界的 “消失”和 “重現”，另一方面還存在冷熱流體相邊界的 “交錯”，使得邊界的移動行為更加復雜。相邊界的 “消失”和 “重現”要求移動邊界模型能夠處理不同的相區模式并能夠光滑轉換，這一部分已有大量的文獻進行了研究［15－20，22－29］；相邊界的 “交錯”則會使得相邊界的排列形式多種多樣。在不同的排列形式下，管壁兩側的冷熱流體相區分布會發生變化，使得傳熱計算的邊界條件不同。因此，需要一種能夠計算熱交換兩側流體均存在相邊界時的傳熱，并且能夠處理相邊界 “交錯”的情形的方法。

因此，要建立繞管式換熱器移動邊界模型，需要：①多股流多相區的換熱關系的數學描述方法；②建立相邊界交錯時的換熱計算模型。

1 多股流多相區換熱關系的

LNG繞管式換熱器通常為級聯形式，級數目根據液化流程有所區別，換熱器結構如圖1所示。

圖1 繞管式換熱器級聯示意圖［5］Fig.1 Schematic diagram of SWHE bundles［5］

在單級管束中，管側流體自管束底部向上流動，殼側流體自管束頂部向下流動。管束內部各個相區長度不一，單個相區同時與多個相區發生換熱，如圖2所示。

圖2 單級管束相區分布示意圖Fig.2 Distribution of phase regions in single bundle

為了描述一個具有n股管側流路、1股殼側流體且每股流路中有多個相區的換熱器，可以通過二維矩陣的形式表示換熱計算中的相關參數。二維矩陣的行和列數目相同，一共由3（n＋1）×3（n＋1）個元素組成，每一行 （列）代表一股流路的單個相區。矩陣中的元素上角標表示當前行對應的流路編號和相區編號，下角標表示當前列對應的流路編號和相區編號。對于繞管式換熱器，殼程流路對應的編號為＃0，管程流路數目為n時對應的流路編號依次為＃1～＃n。每股流路中統一將過熱區編號為＃1，兩相區編號為＃2，過冷區編號為＃3。天然氣流路的超臨界區視為過熱區，編號為＃1。二維矩陣的表達形式如圖3所示。

圖3 多股流多相區二維矩陣Fig.3 Matrix of multiple streams and multiple zones

該二維矩陣可以視作一個分塊矩陣，由 （n＋1）×（n＋1）個子矩陣構成，每一個子矩陣由3×3個元素構成，子矩陣的表達形式如式 （1）所示。

矩陣中的元素可以用換熱計算中具體的物理參數代替。這些參數可以分為兩類：一類是與相區分布特征無關的參數，包括流體溫度Tf、單位長度換熱面積Funit、傳熱系數；另一類是與相區分布特征相關的參數，包括管壁溫度Tw、換熱長度L以及換熱量Q。

當矩陣表示與相區分布特征無關的參數時，矩陣為對角矩陣，只有在矩陣對角線上的元素具有物理意義，代表當前流路當前相區的物理參數。例如，當矩陣表示流體溫度時，矩陣表達形式如式（2）所示。

其中，子矩陣的表達形式如式 （3）所示。

當矩陣表示與相區分布特征相關的參數時，矩陣為互連矩陣，只有當流路i和流路j存在換熱關系時對應的子矩陣中的元素才有值。對角線上的子矩陣以及其他行列對應的流路間不存在換熱關系的子矩陣，里面的元素將沒有值。例如，當矩陣表示換熱長度時，矩陣表達形式如式 （4）所示。

其中，子矩陣Lij的表達形式如式 （5）所示。

這里使用的長度l為歸一化長度，表示當前段長度占總長度的比值。在長度矩陣中，任一元素）具有的物理意義為：流路i的第m個相區與流路j的第n個相區之間的換熱長度。通過單個長度子矩陣可以描述兩股流路中的相區分布情況，如圖4所示。

圖4 管殼側兩股流路的相區分布示意圖Fig.4 Distribution of phase regions in tube and shell side

對應圖4中的相區分布情況，長度子矩陣的表達形式如式 （6）所示。

當矩陣中的元素表示管壁溫度時，其物理意義為：流路i的第m個相區與流路j的第n個相區之間的管壁區間的平均溫度。對于圖3的相區分布情況，溫度子矩陣Tiw0的表達形式如式 （7）所示。

2 相邊界交錯下的傳熱計算模型

相變界交錯下的傳熱計算模型包括流體的控制方程、流體與管壁的傳熱方程和管壁的控制方程。其中，流體的控制方程根據相區的質量和能量守恒建立，采用 的 移 動 邊 界 模 型 的 假 設［15，17－18，20］包 括：①流體的流動視為一維流動；②單股流體內壓力均勻分布，忽略流動過程中的壓力損失；③忽略流體的軸向導熱；④忽略管壁內部的導熱。

根據上述假設，微分形式的流體控制方程如式（8）和式 （9）所示。

式中，ρ表示密度，A表示流體的流通面積，m表示質量流量，h表示比焓，p表示壓力，Q表示流體的換熱量，τ表示時間，z表示流體的流動方向。

將上述偏微分方程沿流體的相區長度積分，可以得到不同展開形式的微分代數方程［30］。對展開后的微分代數方程進行求解，可以得到流體相區的壓力、溫度等狀態參數以及相區長度隨時間的變化。對于相變界交錯下的多股流繞管式換熱器，求解上式的關鍵在于建立單個流體相區與管壁之間的傳熱模型和管壁溫度的計算模型。

2.1 流體與管壁的傳熱模型

流體與管壁的傳熱模型，可以根據第1節中建立的流體溫度、管壁溫度、換熱長度的二維矩陣，使用統一的表達形式進行計算。

2.1.1 兩個相區之間的換熱計算 假定相區A為流路i的第m個相區，相區B為流路j的第n個相區。由相區A傳遞給相區B的換熱量可以通過式（10）進行計算

反之，由相區B傳遞給相區A的熱量為

2.1.2 單個相區的總換熱量計算 假定相區A為流路i的第m個相區，則相區A的總換熱量為相區A與其他所有相區的換熱量之和，可以通過式（12）進行計算

2.2 管壁溫度的計算模型

熱交換兩側的冷熱流體均發生相變時，管壁邊界的排列有多種形式，圖5是其中兩種代表性的邊界排布示意圖。對于管壁邊界，最多可能出現的邊界數為6個，如圖5（a）所示。對于只有兩個相區或者單個相區的情況，保留所有的邊界編號，對于其中實際不存在的相區，將其邊界視為重疊在一起，如圖5（b）所示。

圖5 管壁邊界排布示意圖Fig.5 Boundary arrangement of tube wall

對于其中的單個管壁區間 （區間k），其能量傳遞過程如圖6所示。根據圖6，并應用萊布尼茲積分法則，得到微分形式的控制方程，如式 （13）所示。

圖6 管壁熱量傳遞示意圖Fig.6 Heat transfer in tube wall

式 （19）中，lw，m＋n、Tw，m＋n根據式 （20）計算得到

為了使模型能夠在相邊界發生 “交錯”時光滑切換，虛擬區間的管壁溫度應保持與對應實際區間的管壁溫度一致，通過式 （21）進行計算。

式中，Kw為松弛因子，Tw，track為對應的實際區間的管壁溫度值。

3 計算實例和驗證

3.1 計算實例

選取一個典型流路布置的繞管式換熱器進行實例計算，其流路方案如圖7所示。

設置具體的工藝方案：流路1和流路2為熱流體，流路3為冷流體；流路1的天然氣運行在超臨界區，流路2的混合冷劑運行在亞臨界區并由氣態轉變為兩相，流路3的殼側冷劑運行在亞臨界區并由兩相轉變為氣態。初始穩態時的具體工藝方案見表2。

圖7 計算實例的流路方案Fig.7 Flow arrangement in case study

表2 計算實例的初始工藝方案Table 2 Initial parameters of case

設置擾動條件為：殼側流體的傳熱系數隨時間周期性波動。實例計算中流體與管壁的傳熱計算、管壁溫度的計算采用本研究提出的模型，流體相區的控制方程采用文獻中的模型［19］，傳熱系數采用換熱關聯式進行計算 （表3）。

表3 計算實例中的換熱關聯式Table 3 List of applied heat transfer correlations

3股流體的壓力、出口溫度和相區長度隨時間的變化分別如圖8～圖10所示。

圖8 各股流路的壓力隨時間的變化曲線Fig.8 Graph of pressure in each stream with time

圖9 各股流路的出口溫度隨時間的變化曲線Fig.9 Graph of outlet temperature in each stream with time

圖10 各股流路的相區長度隨時間的變化曲線Fig.10 Graph of zone length in each stream with time

3股流體中，流路2和流路3均發生相變，其相變界面處的邊界隨時間移動的曲線如圖11所示。圖中，兩條曲線發生交叉時，表示流路2和流路3的相邊界發生 “交錯”的情況。

圖11 相界面隨時間的移動曲線Fig.11 Graph of interface move with time

3.2 精度驗證

目前公開文獻中尚未發表過LNG或FLNG繞管式換熱器的相關實驗數據和實際工廠的動態運行數據，因此沒有實驗數據可以用來與模型比對?？紤]到本研究做的創新工作主要集中在如何描述和計算多股流的多相區耦合的問題，對于其中的任意一股流體，本研究用的模型與已有文獻上的分相移動邊界模型一致，因此本研究開發的模型的正確性可以通過已有的單股流的分相移動邊界模型來驗證。

驗證時流體建模采用與文獻中模型相同的P－ˉγ法，流體與管壁的換熱計算和管壁溫度的計算采用本研究提出的模型。驗證實例采用如下工藝方案：管側為單股流體形式，工質采用C2和C3組成的混合冷劑，并設置管側為過熱氣體、殼側為兩相到過熱的相變流體。驗證時傳熱系數取為常數，這樣避免了關聯式誤差對驗證精度的影響，但在實際應用中可以采用表3列出的換熱關聯式或其他關聯式進行計算?？张菹禂的Ｐ筒捎镁嗄Ｐ?，同樣可以由其他模型代替。流體的物性計算采用GERG2008［35］方 法，通 過 REFPROP 6.01 生 成。具體的工藝參數和結構參數見表4。

設置擾動條件為：在初始時換熱器處于穩態，隨后保持管側和殼側的進口焓值不變，增大殼側流量，此處增量取為5%。

仿真結果顯示，兩種模型下殼側及管側的壓力和出口焓值的計算結果吻合良好，最大偏差小于4%，如圖12和圖13所示。

表4 繞管式換熱器精度驗證參數Table 4 SWHE parameters of validation case

圖12 殼側壓力和出口焓值隨時間的變化曲線Fig.12 Graph of shell－side pressure and outlet enthalpy with time

圖13 管側壓力和出口焓值隨時間的變化曲線Fig.13 Graph of tube－side pressure and outlet enthalpy with time

4 結 論

（1）建立的二維矩陣數學描述方法能夠對繞管式換熱器多股流路中各個相區之間的換熱關系進行描述，實現了移動邊界下多股流并行換熱的計算。

（2）建立的傳熱計算模型能夠對冷熱流體并發相變的情況進行計算，并可以處理相邊界的 “交錯”，實現了移動邊界下蒸發和冷凝并存的計算。

（3）實例表明本研究提出的建模方法可以應用于多股流繞管式換熱器的動態仿真，與文獻中模型的仿真結果吻合良好，最大偏差小于4%，滿足精度要求。

符 號 說 明

Aw——管壁橫截面積，m2

Cw——管壁比熱容，J·kg－1·K－1

Funit——單位長度換熱面積，m2·m－1

Kw——時間常數的倒數，s－1

L——換熱長度矩陣

L——總換熱長度，m

l——歸一化的換熱長度

lw——歸一化的管壁區間長度

Q——換熱量，W

Tf——流體溫度矩陣

Tf——流體溫度，K

Tint——管壁區間交界面溫度，K

Tw——管壁溫度矩陣

Tw——管壁溫度，K

α——傳熱系數，W·m－2·K－1

ρw——管壁材料密度，kg·m－3

τ——時間，s

χl——管壁區間的下邊界位置

χu——管壁區間的上邊界位置

χw——管壁區間的邊界位置

下角標

f——流體

in——向內傳遞的

int——管壁區間交界面

i，j——流路編號

k——管壁區間編號

l——管壁區間的下邊界

m，n——相區編號

out——向外傳遞的

unit——單位長度

w——管壁

［1］ Momirlan M，Veziroglu T N.The properties of hydrogen as fuel tomorrow in sustainable energy system for a cleaner planet ［J］.InternationalJournalofHydrogenEnergy，2005，30 （7）：795－802.

［2］ Zhang Fengjiu （張鳳久）.Outlook for natural gas development in the South China Sea ［J］.Natur.GasInd.（天然氣工業），2009，29 （1）：17－20.

［3］ Bukowski J，Liu Y N，Boccella S，etal.Innovations in natural gas liquefaction technology for future LNG plants and floating LNG facilities／／Proceedings of the International Gas Union Research Conference ［C］.2011.

［4］ Pacio J C，Dorao C A.A review on heat exchanger thermal hydraulic models for cryogenic applications ［J］.Cryogenics，2011，51 （7）：366－379.

［5］ Wang T，Ding G，Duan Z，etal.A distributed－parameter model for LNG spiral wound heat exchanger based on graph theory ［J］.AppliedThermalEngineering，2015，81（25）：102－113.

［6］ Lex T，Ohlig K，fredheim A O，etal.Ready for floating LNG－qualification of spiral wound heat exchangers／／Proceedings 15th International Conferences on Liquefied Natural Gas ［C］.Barcelona，Spain，2007：PO－27.1－PO－27.14.

［7］ Luo X，Guan X，Li M，etal.Dynamic behaviour of onedimensional flow multistream heat exchangers and their networks ［J］.InternationalJournalofHeatandMass Transfer，2003，46 （4）：705－715.

［8］ Mishra M，Das P K，Sarangi S.Dynamic behavior of threefluid crossflow heat exchangers ［J］.JournalofHeat Transfer，2008，130 （1）：118－125.

［9］ Das S K，Murugesan K.Transient response of multipass plate heat exchangers with axial thermal dispersion in fluid［J］.InternationalJournalofHeatandMassTransfer，2000，43 （23）：4327－4345.

［10］ Zhang Qin （張勤），Li Zhixin （李志信），Liang Xingang（梁新剛）.Field synergy analysis of the dynamic process of a multi－stream heat exchanger ［J］.JournalofEngineering forThermalEnergy＆Power（熱能動力工程），2009，24（6）：782－786.

［11］ Hammer M.Dynamic simulation of a natural gas liquefaction plant ［D］.Trondheim：Norwegian University of Science and Technology，2004：55－108.

［12］ Xu Yifeng （徐益峰），Cai Zuhui（蔡祖恢）.Research on dynamic model of multi－fluid plate－fin heat exchangers ［J］.JournalofChemicalIndustryandEngineering（China）（化工學報），1998，49 （6）：721－728.

［13］ Lei Z，Zaheeruddin M.Dynamic simulation and analysis of a water chiller refrigeration system ［J］.AppliedThermal Engineering，2005，25 （14）：2258－2271.

［14］ Ding G L.Recent developments in simulation techniques for vapour－compression refrigeration systems ［J］.International JournalofRefrigeration，2007，30 （7）：1119－1133.

［15］ Munch Jensen J，Tummescheit H.Moving boundary models for dynamic simulations of two－phase flows／／Modelica’2002 Proceedings ［C］.2002：235－244.

［16］ Rasmussen B P，Alleyne A G.Control－oriented modeling of transcritical vapor compression systems ［J］.Journalof DynamicSystems，Measurement，andControl，2004，126（1）：54－64.

［17］ Zhang W J，Zhang C L.A generalized moving－boundary model for transient simulation of dry－expansion evaporators under larger disturbances ［J］.InternationalJournalof Refrigeration，2006，29 （7）：1119－1127.

［18］ McKinley T L，Alleyne A G.An advanced nonlinear switched heat exchanger model for vapor compression cycles using the moving－boundary method ［J］.InternationalJournalof Refrigeration，2008，31 （7）：1253－1264.

［19］ Li B，Alleyne A G.A dynamic model of a vapor compression cycle with shut－down and start－up operations ［J］.InternationalJournalofRefrigeration，2010，33 （3）：538－552.

［20］ Cecchinato L，Mancini F.An intrinsically mass conservative switched evaporator model adopting the moving－boundary method ［J］.InternationalJournalofRefrigeration，2012，35 （2）：349－364.

［21］ Picon－Nunez M，Polley G T，Medina－Flores M.Thermal design of multi－stream heat exchangers ［J］.Applied ThermalEngineering，2002，22 （14）：1643－1660.

［22］ Willatzen M，Pettit N，Ploug－S?rensen L.A general dynamic simulation model for evaporators and condensers in refrigeration （Ⅰ ）： Moving－boundary formulation of twophase flows with heat exchange ［J］.InternationalJournal ofRefrigeration，1998，21 （5）：398－403.

［23］ Ding G，Zhang C，Lu Z.Dynamic simulation of natural convection bypass two－circuit cycle refrigerator－freezer and its application（Ⅰ）：Component models ［J］.AppliedThermalEngineering，2004，24 （10）：1513－1524.

［24］ Fu L，Ding G，Zhang C.Dynamic simulation of air－to－water dual－mode heat pump with screw compressor ［J］.Applied ThermalEngineering，2003，23 （13）：1629－1645.

［25］ Schurt L C，Hermes C J L，Neto A T.A model－driven multivariable controller for vapor compression refrigeration systems ［J］.InternationalJournalofRefrigeration，2009，32 （7）：1672－1682.

［26］ Pettit N，Willatzen M，Ploug－S?rensen L.A general dynamic simulation model for evaporators and condensers in refrigeration （Ⅱ）：Simulation and control of an eva－porator［J］.InternationalJournalofRefrigeration，1998，21 （5）：404－414.

［27］ Beghi A，Cecchinato L.A simulation environment for dryexpansion evaporators with application to the design of autotuning control algorithms for electronic expansion valves［J］.InternationalJournalofRefrigeration，2009，32 （7）：1765－1775.

［28］ Liang N，Shao S，Tian C，etal.Dynamic simulation of variable capacity refrigeration systems under abnormal conditions ［J］.AppliedThermalEngineering，2010，30（10）：1205－1214.

［29］ Bendapudi S，Braun J E，Groll E A.A comparison of moving－boundary and finite－volume formulations for transients in centrifugal chillers ［J］.InternationalJournal ofRefrigeration，2008，31 （8）：1437－1452.

［30］ Zhao Dan （趙 丹），Ding Guoliang （丁 國 良），Hu Haitao（胡 海 濤 ）.An intrinsically mass and energy conservative model for transient simulation of evaporator ［J］.Journalof Refrigeration（制冷學報），2015，36 （1）：76－83.

［31］ Shah M.A general correlation for heat transfer during film condensation inside pipes ［J］.InternationalJournalofHeat＆MassTransfer，1979，22 （4）：547－556.

［32］ Boyko L D，Kruzhilin G N.Heat transfer and hydraulic resistance during condensation of steam in a horizontal tube and in a bundle of tubes ［J］.InternationalJournalofHeat＆MassTransfer，1967，10 （67）：361－373.

［33］ Neeraas B O，Fredheim A O，Aunan B.Experimental shellside heat transfer and pressure drop in gas flow for spiralwound LNG heat exchanger ［J］.InternationalJournalof Heat＆MassTransfer，2004，47 （2）：353－361.

［34］ Neeraas B O，Fredheim A O，Aunan B.Experimental data and model for heat transfer，in liquid falling film flow on shell－side，for spiral－wound LNG heat exchanger ［J］.InternationalJournalofHeat＆MassTransfer，2004，47（14）：3565－3572.

［35］ Kunz O，Wagner W.The GERG－2008wide－range equation of state for natural gases and other mixtures：an expansion of GERG－2004 ［J］.JournalofChemical＆Engineering Data，2012，57 （11）：3032－3091.