一种管内超临界流体气化‑管外膜状冷凝传热计算方法

技术领域

本发明属于低温流体超临界传热技术领域，具体涉及一种管内超临界流体气化-管外膜状冷凝传热计算方法。

背景技术

将工业气体液化以液态形式储运、并在用户终端使其再气化的技术被广泛应用于液化天然气(LNG)、液氧、液氮、液氢、液态二氧化碳等领域，可大大节省储运空间、提高储运效率。低温流体气化器是实现该过程的关键装备之一。

由于全球范围内的能源危机和环境恶化，LNG工业快速发展。为满足LNG接收站气化低温LNG的需要，专利EP0550845B1提出一种开架式气化器(ORV)结构，通过将海水自上而下喷淋在板型换热管束外表面上，把能量传递给LNG使其气化；专利US3818893提出一种浸没燃烧式气化器(SCV)型式，将天然气燃烧后产生的高温烟气排入水浴池中加热并剧烈搅拌海水，浸泡在海水中的换热盘管内LNG被池水加热气化；US6164247提出一种中间流体式气化器(IFV)结构，将蒸发器、凝结器与调温器三组管壳式换热器组合起来，在蒸发器中用海水加热壳程的低沸点中间介质使其气化，在凝结器中间介质管外冷凝提供热量气化管内的LNG，在调温器中LNG被进一步加热以达到要求的温度。目前，ORV、SCV与IFV已成为世界范围内广泛应用的三种基本负荷型LNG气化器。

IFV是一种具有广泛应用潜力的低温流体气化设备。IFV相比于ORV，它避免了用海水直接气化LNG时可能带来的结冰问题，且适用于海水含沙量较高的水域，因而适应性更强。IFV相比于SCV，它不需要额外燃烧天然气，因而能效更好，且没有CO₂、NO_x等环境污染问题。由于IFV的组合式气化结构设计，它可用于海水水质较差的陆上LNG气化终端、空间和重量受限的海上浮式储存与气化平台、以及冷能发电系统等不同场合。

为进一步提升IFV在工程应用中的技术优势，专利US6164247提出在蒸发器/冷凝器共用的壳体内设置分隔板和支撑板，以适应海上平台晃荡工况，保持稳定的气化效率；专利US6367429B2提出当采用工业回收热水作为热源时，可将蒸发器设计成两管程或多管程，以使结构更紧凑；专利US6367265B1提出在凝结器管箱中设置LNG喷雾冷却装置，以降低IFV启动或低负荷运行时的温差应力；专利US20160146403A1提出凝结器采用管外带翅片或沟槽的强化传热管，以强化中间介质的管外冷凝传热过程。

由于IFV采用组合式气化结构、设备紧凑、LNG温度及物性变化范围大，无法采用现有商业软件直接设计。文献中主要采用分布式参数模型(DPM)进行设计型和校核型计算。该方法是将IFV各传热段沿管长方向划分成多个微元，对每个微元进行能量衡算和总传热系数核算，其假设整个换热过程是稳态的、不存在向环境的漏热。文献中一般将各传热段均按单管程、单壳程处理，并忽略压降。例如，文献[大型LNG中间介质气化器换热面积计算方法,天然气与石油,31,31-35,2013]提出了超临界工况下IFV各传热段换热面积计算方法；文献[Thermal performance analysis of intermediate fluid vaporizer for liquefiednatural gas,Applied Thermal Engineering,65,564-574,2014]对超临界工况IFV进行了校核型计算，分析了海水进口温度与流量、LNG进口压力与流量等工艺参数对传热特性的影响。

但上述现有技术忽略了管内超临界流体流经换热管的压力降、不适用于采用强化传热管时IFV的传热计算、且采用的管内超临界流体传热系数关联式预测精度不高，因而传热计算精度低、适用范围受到限制。

发明内容

为了解决上述技术问题，本发明提供一种管内超临界流体气化-管外膜状冷凝传热计算方法。

为了实现本发明的目的，本发明采用了以下技术方案：

一种管内超临界流体气化-管外膜状冷凝传热计算方法，包括以下步骤：

S1、输入传热工艺参数与水平换热管参数，包括：

管内流体的组分及其含量、入口温度T_b⁽¹⁾、入口压力P_b⁽¹⁾、质量通量G_b、拟临界温度T_pc，其中管内流体的入口压力大于管内流体的临界压力；

管外流体的组分及其含量、饱和温度T_sat；

换热管的内直径d_i、外直径d_o、长度L、导热系数k_w、单列换热管束中沿竖直方向的换热管排数N_eff；

S2、根据换热管长度划分N个相等长度的传热微元，每个传热微元长度为dL＝L/N；所述传热微元数目N满足：管内流体进出换热管温度差/N≤20K；

假设所述换热管外壁面过冷度ΔT_sub的初值记为ΔT_sub,0，所述ΔT_sub,0取0.5～20K；

S3、计算传热微元j(j＝1,2,…,N)外壁面温度、管外冷凝传热系数及热负荷初值

根据换热管外壁面过冷度初值ΔT_sub,0，按公式(1)～(3)计算传热微元j外壁面温度T_wo、管外冷凝传热系数h_o及热负荷dQ的初值，并分别记为T_wo,0^(j)、h_o,0^(j)与dQ₀^(j)：

T_wo＝T_sat-ΔT_sub(1)

dQ＝h_oΔT_subπd_odL(3)

式中：g为重力加速度；

ρ_L、ρ_V、k_L、r、μ_L分别为管外流体在饱和温度T_sat下的液相密度、汽相密度、液相导热系数、蒸发焓、液相粘度；

m为管束效应修正系数，按如下方式进行取值：

在换热器样机中选取具有X排换热管的单列换热管束，按照GB/T27698.2-2011汽-液冷凝器用热交换器的测试程序，测试每根换热管内流体的进出口温度、压力和流量，管外流体的进出口温度、压力和流量，并计算沿竖直方向自上而下布置的第1排换热管外的冷凝传热系数h_o⁽¹⁾与第M排换热管外的冷凝传热系数h_o^(M)，其中X>M，按公式(4)回归得到m值：

h_o^(M)＝h_o⁽¹⁾·[M^1-m-(M-1)^1-m](4)

S4、计算传热微元j内壁面温度、管内超临界流体传热系数

根据传热微元j外壁面温度初值T_wo,0^(j)和热负荷初值dQ₀^(j)，按公式(5)计算内壁面温度T_wi^(j)；根据传热微元j的入口温度T_b^(j)、入口压力P_b^(j)，按公式(6)计算管内超临界流体传热系数h_i^(j)：

式中：

T_b以管内流体在传热微元j入口处的温度T_b^(j)代入；

ρ_b、C_pb、k_b、μ_b、H_b分别以管内流体在传热微元j入口温度T_b^(j)、入口压力P_b^(j)下的密度、定压比热、导热系数、粘度和焓值代入；

ρ_wi、H_wi分别以管内流体在传热微元j内壁面温度T_wi^(j)、入口压力P_b^(j)下的密度和焓值代入；

S5、核算传热微元j热负荷

根据传热微元j内壁面温度T_wi^(j)、管内流体传热系数h_i^(j)，按公式(7)核算传热微元j的热负荷dQ^(j)：

dQ＝h_i(T_wi-T_b)πd_idL(7)

式中，T_b以管内流体在传热微元j入口处的温度T_b^(j)代入；

S6、更新传热微元j外壁面温度、管外冷凝传热系数

以核算后的传热微元j热负荷dQ^(j)代入公式(1)～(3)，求解更新后的传热微元j外壁面温度T_wo^(j)、管外冷凝传热系数h_o^(j)；

S7、判断传热微元j传热计算是否收敛

根据传热微元j热负荷的初值dQ₀^(j)与核算值dQ^(j)，按公式(8)计算相对误差ε：

若ε≤0.001，则传热微元j传热计算收敛，进入S8；

否则，返回S4，依据更新后的外壁面温度和热负荷，重新计算内壁面温度、管内超临界流体传热系数；

S8、计算传热微元j管内流体压力降

根据传热微元j的入口温度T_b^(j)、入口压力P_b^(j)，按公式(9)～(10)计算传热微元j管内流体压力降ΔP_b，并记为ΔP_b^(j)：

式中：f为直管摩擦阻力系数；

R为U形换热管的弯曲半径；当采用直换热管时，取R＝0；

S9、计算传热微元j管内流体出口压力与温度

根据传热微元j入口处管内流体的压力P_b^(j)和流经传热微元j的压力降ΔP_b^(j)，按公式(11)计算传热微元j出口处管内流体的压力P_b^(j+1)；根据传热微元j入口处管内流体的焓值H_b^(j)，按公式(12)计算传热微元j出口处管内流体的焓值H_b^(j+1)：

根据传热微元j出口处管内流体的焓值H_b^(j+1)和压力P_b^(j+1)，通过查询物质理化性能手册或采用热力学方程计算方法，得到传热微元j出口处管内流体的温度T_b^(j+1)；

S10、判断是否计算至最后一个传热微元，即当前j值是否等于N

若是，则进入步骤S12；若否，则进入步骤S11；

S11、进行下一传热微元的传热计算

进入下一传热微元j+1，并执行步骤S3；

S12、完成传热计算。

进一步，所述步骤S3中m值根据换热管的类型按照如下方式进行取值：

对于普通光滑管，m＝0.1667～0.25；

对于二维外翅片管，m＝0.01～0.1；

对于三维非连续翅片管，m＝0.23～0.36。

进一步，所述换热器由多列、多排换热管组成时，所述步骤S3中的管外冷凝传热系数h_o按公式(13)计算平均值：

式中，y表示第y列换热管束；Y表示换热器共有Y列换热管束；h_0,y表示第y列换热管束的管外冷凝传热系数；N_eff,y表示第y列换热管束中沿竖直方向的换热管排数。

本发明的有益效果在于：

本发明提出了一种管内超临界流体气化-管外膜状冷凝传热计算方法，尤其是一种卧式换热器中管内超临界流体气化-管外膜状冷凝传热过程的校核型计算方法，考虑了管内超临界流体的压力降、管外膜状冷凝传热的管束效应影响，采用了高精度的传热系数关联式，适用于流体进出口温差大、采用普通光滑换热管或不同类型强化传热管工况，计算精度高。

附图说明

图1为本发明中管内超临界流体气化-管外膜状冷凝传热计算流程示意图。

图2为本发明实施例1计算的某LNG气化过程中沿管长方向的温度分布与管内超临界流体传热系数分布结果。

图3为本发明实施例2计算的某LNG气化过程中沿管长方向的温度分布与管内超临界流体传热系数分布结果。

图4为本发明实施例3计算的某LNG气化过程中沿管长方向的温度分布与管内超临界流体传热系数分布结果。

具体实施方式

下面结合结合附图和实施例对本发明技术方案做出更为具体的说明：

如附图1所示，一种管内超临界流体气化-管外膜状冷凝传热计算方法，适用于低温流体在超临界压力下气化时卧式换热器的传热校核计算，包括以下步骤：

S1、输入传热工艺参数与水平换热管参数，包括：

管内流体的组分及其含量、入口温度T_b⁽¹⁾、入口压力P_b⁽¹⁾、质量通量G_b、拟临界温度T_pc，其中管内流体的入口压力大于管内流体的临界压力；

管外流体的组分及其含量、饱和温度T_sat；

换热管的内直径d_i、外直径d_o、长度L、导热系数k_w、单列换热管束中沿竖直方向的换热管排数N_eff；

所述管内流体的组分及其含量、管外流体的组分及其含量主要用于计算相应流体的物理性质参数，比如密度、粘度、导热系数、焓等。

所述含量可以是质量含量或体积含量或摩尔含量。

S2、根据换热管长度划分N个相等长度的传热微元，每个传热微元长度为dL＝L/N；所述传热微元数目N满足：管内流体进出换热管温度差/N≤20K；

假设所述换热管外壁面过冷度ΔT_sub的初值记为ΔT_sub,0，所述ΔT_sub,0取0.5～20K；

S3、计算传热微元j(j＝1,2,…,N)外壁面温度、管外冷凝传热系数及热负荷初值

根据换热管外壁面过冷度初值ΔT_sub,0，按公式(1)～(3)计算传热微元j外壁面温度T_wo、管外冷凝传热系数h_o及热负荷dQ的初值，并分别记为T_wo,0^(j)、h_o,0^(j)与dQ₀^(j)：

T_wo＝T_sat-ΔT_sub(1)

dQ＝h_oΔT_subπd_odL(3)

式中：g为重力加速度；

ρ_L、ρ_V、k_L、r、μ_L分别为管外流体在饱和温度T_sat下的液相密度、汽相密度、液相导热系数、蒸发焓、液相粘度；

m为管束效应修正系数，按如下方式进行取值：

在换热器样机中选取具有X排换热管的单列换热管束，按照GB/T27698.2-2011汽-液冷凝器用热交换器的测试程序，测试每根换热管内流体的进出口温度、压力和流量，管外流体的进出口温度、压力和流量，并计算沿竖直方向自上而下布置的第1排换热管外的冷凝传热系数h_o⁽¹⁾与第M排换热管外的冷凝传热系数h_o^(M)，按公式(4)回归得到m值：

h_o^(M)＝h_o⁽¹⁾·[M^1-m-(M-1)^1-m](4)

换热器管束外冷凝传热过程存在边界效应，即末排换热管外冷凝传热系数不降反升。为避免边界效应对冷凝传热试验过程的影响，所述换热器样机中单列换热管束的换热管总排数X需大于待测换热管所在排数M。

S4、计算传热微元j内壁面温度、管内超临界流体传热系数

根据传热微元j外壁面温度初值T_wo,0^(j)和热负荷初值dQ₀^(j)，按公式(5)计算内壁面温度T_wi^(j)；根据传热微元j的入口温度T_b^(j)、入口压力P_b^(j)，按公式(6)计算管内超临界流体传热系数h_i^(j)：

式中：

T_b以管内流体在传热微元j入口处的温度T_b^(j)代入；

ρ_b、C_pb、k_b、μ_b、H_b分别以管内流体在传热微元j入口温度T_b^(j)、入口压力P_b^(j)下的密度、定压比热、导热系数、粘度和焓值代入；

ρ_wi、H_wi分别以管内流体在传热微元j内壁面温度T_wi^(j)、入口压力P_b^(j)下的密度和焓值代入；

所述传热微元j的入口温度T_b^(j)、入口压力P_b^(j)，是由步骤S9计算得到的传热微元(j-1)的出口温度和出口压力；对于传热微元1即是指整根换热管的入口温度T_b⁽¹⁾、入口压力P_b⁽¹⁾；

S5、核算传热微元j热负荷

根据传热微元j内壁面温度T_wi^(j)、管内流体传热系数h_i^(j)，按公式(7)核算传热微元j的热负荷dQ^(j)：

dQ＝h_i(T_wi-T_b)πd_idL(7)

式中，T_b以管内流体在传热微元j入口处的温度T_b^(j)代入；

S6、更新传热微元j外壁面温度、管外冷凝传热系数

以核算后的传热微元j热负荷dQ^(j)代入公式(1)～(3)，求解更新后的传热微元j外壁面温度T_wo^(j)、管外冷凝传热系数h_o^(j)；

S7、判断传热微元j传热计算是否收敛

根据传热微元j热负荷的初值dQ₀^(j)与核算值dQ^(j)，按公式(8)计算相对误差ε：

若ε≤0.001，则传热微元j传热计算收敛，进入S8；

否则，返回S4，依据更新后的外壁面温度和热负荷，重新计算内壁面温度、管内超临界流体传热系数；

S8、计算传热微元j管内流体压力降

根据传热微元j的入口温度T_b^(j)、入口压力P_b^(j)，按公式(9)～(10)计算传热微元j管内流体压力降ΔP_b，并记为ΔP_b^(j)：

式中：f为直管摩擦阻力系数；

R为U形换热管的弯曲半径；当采用直换热管时，取R＝0；

S9、计算传热微元j管内流体出口压力与温度

根据传热微元j入口处管内流体的压力P_b^(j)和流经传热微元j的压力降ΔP_b^(j)，按公式(11)计算传热微元j出口处管内流体的压力P_b^(j+1)；根据传热微元j入口处管内流体的焓值H_b^(j)，按公式(12)计算传热微元j出口处管内流体的焓值H_b^(j+1)：

根据传热微元j出口处管内流体的焓值H_b^(j+1)和压力P_b^(j+1)，通过查询物质理化性能手册或采用热力学方程计算方法，得到传热微元j出口处管内流体的温度T_b^(j+1)；

S10、判断是否计算至最后一个传热微元，即当前j值是否等于N

若是，则进入步骤S12；若否，则进入步骤S11；

S11、进行下一传热微元的传热计算

进入下一传热微元j+1，并执行步骤S3；

S12、完成传热计算。

所述步骤S3中m值根据换热管的类型按照如下方式进行取值：

对于普通光滑管，m＝0.1667～0.25；

对于二维外翅片管，m＝0.01～0.1；

对于三维非连续翅片管，m＝0.23～0.36。

换热管通常分为普通光滑管和强化传热管两种类型，所述强化传热管按表面结构的不同，又分为二维外翅片管、异形强化管、三维非连续翅片管。其中所述普通光滑管为不带翅片且壁面光滑的换热管；所述二维外翅片管是将铝带、铜带或钢带构成的翅片缠绕在换热管表面，所述二维外翅片管的翅片表面规则；所述三维非连续翅片管的翅片比二维外翅片管的翅片具有更不规则的扩展表面。

所述换热器由多列、多排换热管组成时，所述步骤S3中的管外冷凝传热系数h_o按公式(13)计算平均值：

式中，y表示第y列换热管束；Y表示换热器共有Y列换热管束；h_0,y表示第y列换热管束的管外冷凝传热系数；N_eff,y表示第y列换热管束中沿竖直方向的换热管排数。

实施例1

以某LNG气化过程为例，依据本发明计算方法，完成该气化过程中换热管的校核型传热计算，具体步骤如下：

S1、输入传热工艺参数与换热管参数

LNG气化工况1中，LNG组成见表1，传热工艺参数见表2，换热管参数见表3。

表1 LNG气化工况1的LNG组分及其含量

表2 LNG气化工况1的传热工艺参数

表3 LNG气化工况1的换热管参数

S2、根据换热管长度划分N个相等长度的传热微元

传热微元数目N＝200；每个传热微元长度为dL＝9/200＝0.045m；

假设换热管外壁面过冷度ΔT_sub的初值，记为ΔT_sub,0，取10K；

S3、计算传热微元j(j＝1,2,…,N)外壁面温度、管外冷凝传热系数及热负荷初值

以j＝1为例：

按公式(1)计算，T_wo,0⁽¹⁾＝260.15K；

对于光滑换热管，取m＝0.1667，按公式(2)计算，h_o,0⁽¹⁾＝2126.95W·m^-2·K^-1；

按公式(3)计算，dQ₀⁽¹⁾＝57.13W；

S4、计算传热微元j内壁面温度、管内超临界流体传热系数

以j＝1为例：

按公式(5)计算，T_wi⁽¹⁾＝257.20K；

传热微元1入口温度为T_b⁽¹⁾＝111.15，入口压力P_b⁽¹⁾＝12000kPa，按公式(6)计算，此时T_b⁽¹⁾<T_pc<T_wi⁽¹⁾，因此n＝0.4+0.2(T_wi/T_pc-1)＝0.4093，h_i⁽¹⁾＝957.14W·m^-2·K^-1；

S5、核算传热微元j热负荷

以j＝1为例：

按公式(7)核算，dQ⁽¹⁾＝296.44W；

S6、更新传热微元j外壁面温度、管外冷凝传热系数

以j＝1为例：

将dQ⁽¹⁾＝296.44W带入公式(1)～(3)，求解得更新后的T_wo⁽¹⁾＝180.32K，h_o⁽¹⁾＝1228.58W·m^-2·K^-1；

S7、判断传热微元j传热计算是否收敛

以j＝1为例：

按公式(8)计算，相对误差ε＝0.807>0.001，传热计算不收敛，返回S4，按T_wo,0⁽¹⁾＝180.32K、dQ₀⁽¹⁾＝296.44W，重新计算内壁面温度、管内超临界流体传热系数；

重复上述迭代过程，最终传热计算收敛时结果为：

T_wo⁽¹⁾＝216.14K，h_o⁽¹⁾＝1477.59W·m^-2·K^-1，dQ⁽¹⁾＝214.34W，T_wi⁽¹⁾＝205.04K，h_i⁽¹⁾＝1076.60W·m^-2·K^-1；

S8、计算传热微元j管内流体压力降

以j＝1为例：

按公式(9)～(10)计算，ΔP_b⁽¹⁾＝0.0034kPa；

S9、计算传热微元j管内流体出口压力与温度

以j＝1为例：

按公式(11)计算，P_b⁽²⁾＝11999.997kPa；

按公式(12)计算，H_b⁽²⁾＝-18310.97J/kg；

根据传热微元1出口处管内流体的焓值H_b⁽²⁾和压力P_b⁽²⁾，通过调用美国国家标准与技术研究院开发的工质物性计算软件(NIST>b⁽²⁾＝113.15K；

S10、判断是否计算至最后一个传热微元

当前j值为1，不等于N，即未计算至最后一个传热微元，进入S11；

S11、进入下一传热微元

令j＝2，进入下一传热微元，并执行S3；

如此依序计算每个传热微元，直到完成所有传热微元的传热计算。

S12、完成传热计算。

本实施例中，利用本发明方法计算的管内LNG沿管长方向的温度分布与管内超临界流体传热系数分布结果如图2所示。LNG最终的出口温度为255.12K，LNG进、出换热管总的压力降为0.79kPa，管内超临界流体传热系数为1076.60～1558.19W·m^-2·K^-1。

实施例2

LNG气化工况2，LNG组成与传热工艺参数同表1与表2，换热管参数见表4。

表4 LNG气化工况2的换热管参数

按本发明的传热计算方法，划分传热微元数目N＝200，取外壁面过冷度初值ΔT_sub,0＝0.5K，取管束效应修正系数m＝0.1667。

本实施例中，利用本发明方法计算的管内LNG沿管长方向的温度分布与管内超临界流体传热系数分布结果如图3所示。LNG最终的出口温度为248.84K，LNG进、出换热管总的压力降为7.08kPa，管内超临界流体传热系数为1094.65～1577.07W·m^-2·K^-1。

实施例3

LNG气化工况3，LNG组成如表5所示，传热工艺参数见表6，换热管参数见表7。

表5 LNG气化工况3的LNG组分及其含量

表6 LNG气化工况3的传热工艺参数

表7 LNG气化工况3的换热管参数

按本发明的传热计算方法，划分传热微元数目N＝200，取外壁面过冷度初值ΔT_sub,0＝1K，取管束效应修正系数m＝0.04。

本实施例中，利用本发明方法计算的管内LNG沿管长方向的温度分布与管内超临界流体传热系数分布结果如图4所示。LNG最终的出口温度为252.66K，LNG进、出换热管总的压力降为0.89kPa，管内超临界流体传热系数为1176.25～1567.55W·m^-2·K^-1。