分相的多股流LNG绕管式换热器动态模型

段钟弟，任滔，丁国良，陈杰，浦晖，密晓光

（1上海交通大学制冷与低温工程研究所，上海200240；2中海石油气电集团技术研发中心，北京100028）

引 言

天然气作为一种高效清洁的能源，近年来得到了广泛的应用［1］。我国南海拥有丰富的待开发天然气资源，已探明的储量占我国天然气总储量的1／3［2］。开采海上天然气的基础装置是浮式天然气生产储卸平台 （LNG－FPSO）。LNG－FPSO的主低温换热器为绕管式换热器，具有占地面积小、力学性能好和抗泄漏能力强等优点［3］。

绕管式换热器由数以千计的管子按层排布螺旋缠绕在中心筒上，单个换热器投资巨大，不能基于实验进行设计，需要先进的设计方法［4］。基于仿真的设计方法具有高效、低成本的优势，能够应用于绕管式换热器设计［5］。LNG－FPSO是海上浮动平台，会随风浪在海上晃动，从而造成绕管式换热器晃动，使得其内部流体的温度、压力随时间发生波动［6］，因此需要能够预测绕管式换热器动态特性的模型。

因天然气液化独特的工艺和需求，绕管式换热器的结构和工作具有如下特征，包括：①工质是由烃类和氮气组成的混合冷剂；②工况范围覆盖包含亚临界和超临界的所有相区；③多股流体并行换热；④单个换热器内冷热流体并发相变，同时存在蒸发和冷凝。基于集总参数的多股流换热器动态模型［7－10］不能反映绕管式换热器中的相变特征；基于分布参数的多股流动态模型［11－12］涉及大量微元的计算，计算耗时长［13］。因此，本研究采用具有良好精度和计算速度的移动边界模型［14］来预测绕管式换热器的动态特性。现有的移动边界模型可以对绕管式换热器前两项特征中的多组分工质、多相区模式进行计算，具体见表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.