浮式氢能储运过程中FLH2通道管外降膜流动的海上适应性强化机理

孙崇正，李玉星，许洁，韩辉，宋光春，卢晓

（1 山东科技大学储能技术学院，山东 青岛 266590；2 中国石油大学(华东)储运与建筑工程学院，山东省油气储运安全重点实验室，山东 青岛 266580；3 国家管网集团北京管道有限公司，北京 100101）

氢能具有清洁、高效的优点。在“碳达峰、碳中和”背景下，全球主要能源消费国相继把氢能提升到国家能源战略层面[1]。随着海上制氢和碳捕获、利用与封存（carbon capture, utilization and storage，CCUS）技术的发展，利用海上风电资源或天然气制备氢气，并通过储运技术送到氢能源市场，为解决海上风电并网和消纳的难题、促进深海天然气资源的低碳发展提供了可行的思路[2]。国际海洋工程界提出一种参考浮式液化天然气生产装置（FLNG）的浮式氢气液化装置（FLH2），浮式氢气液化装置通过将氢气直接在海上液化并装船运输，简化了海上氢能的储运过程，极大地提高了设备使用的灵活性和安全性[3-6]。但受海上风浪的影响，应用于浮式氢气液化装置的制冷剂在降膜流动换热过程中，会出现冷剂分配不均、薄液膜失稳等问题，进而影响液化工艺的性能指标，因此为强化浮式氢气液化装置的海上适应性，亟需开展相关领域的研究工作。

绕管式换热器具有传热系数高、冷剂充灌量少、结构紧凑等优点，被广泛应用于天然气和氢气液化领域。绕管式换热器主要由缠绕管路和壳体两部分组成，换热通过壳侧制冷剂在绕管管壁和管间的降膜沸腾来实现，为绕管内的工质提供冷量。Jian 等[7]搭建了整体绕管式换热器的实验装置，通过实验研究尺寸参数对换热器换热性能的影响。Zhang 等[8]通过低温实验和数值模拟相结合的方法研究了超低温氦在绕管式换热器的换热特性。Chien 等[9]通过实验研究了光管、翅片管和沸腾强化管的管外降膜沸腾换热特性。Hu 等[10]通过实验方法研究了乙烷/丙烷混合制冷剂在绕管式换热器壳侧流动沸腾的传热机理。Geir等[11]对应用于氢气液化工艺中板翅式换热器和绕管式换热器的换热性能进行理论研究，结果表明绕管式换热器中正-仲氢转化的 损失较小。GeniĆ等[12-13]通过实验研究优化了绕管式换热器壳侧换热计算方程。Lu 等[14-15]通过数值模拟方法研究了管间距、管外径、管层数以及中心筒直径对换热器性能的影响。Wang等和Jian等[16-18]通过多目标遗传算法分析了几何参数对绕管式换热器性能的影响。Yu 等[19-21]构建了数值模拟模型，研究了螺旋管内烷烃类工质的冷凝流动与换热特性。Li 等[22-24]分别研究了常规圆管、波纹管、方波管和螺旋强化管内烷烃类介质的冷凝传热特性，研究结果表明，相较于常规圆管，波纹管、方波管和螺旋强化管的管型结构均可提高管内冷凝换热特性。Ding 等[25-27]研究了烷烃类介质在换热器壳侧纯气相流动和降膜流动过程的换热与压降性能。

与陆上的应用条件不同，绕管式换热器在浮式海上平台应用时，受多自由度倾斜晃荡海况的影响，换热器内部会出现液体偏流的现象，导致管外制冷剂分布不均，管壁局部液膜较薄甚至出现干涸现象，引起换热恶化。为提高绕管式换热器的海上适应性，Zheng 等[28]对其分配器性能进行研究，提出一种新型气液分配器结构。Wang 等[29]建立了FLNG 绕管式换热器的三维数学模型，研究海上晃荡条件下换热器传热传质过程。Duan 等[30]建立了分相多股流动态模型，预测海况下FLNG绕管式换热器不同相态区域之间的相变换热规律。Li等[31]研究了晃荡条件下绕管式换热器盘管内烷烃类介质冷凝过程中的流动与换热特性。Ren 等[32-33]通过数值模拟研究发现，相较于纯气相流动换热，晃荡对降膜流动换热过程影响较大。

应用于浮式氢气液化装置的新型绕管式换热器通道管内由于填充了正-仲氢转化催化剂，导致通道管内流体区域狭小，而海况对有自由表面的管外流体降膜流动过程影响较大，进而影响其传热传质过程，降低FLH2的海上适应性。因此，为强化FLH2装置的海上适应性，本文基于高速摄像系统和晃荡平台系统以及3D 建模打印技术，同时基于耦合的VOF 和Level Set 两相流模型和动网格技术，搭建了降膜流动浮式可视化实验装置和数值模拟模型。采用浮式微观可视化实验和数值模拟相结合的方法，研究FLH2中降膜流动过程的海上适应性强化机理，提出降膜流动强化措施。

1 降膜流动可视化实验装置

1.1 实验装置建立

氢气液化工艺中制冷循环由预冷循环和深冷循环两部分组成，深冷循环以氦等工质的膨胀制冷为主，预冷循环中冷剂种类较多，包括氮气，甲烷、乙烷、丙烷、正丁烷等烷烃类工质，将深冷循环中的氦气和原料氢气从常温25℃冷却至-193℃左右，预冷循环的稳定运行对整个氢气液化系统非常重要。预冷循环的非共沸制冷剂在降膜换热过程中，存在薄液膜界面，海况对含有自由液膜的流动过程影响较大，通过影响降膜流动状态与液膜分布形态，进而影响预冷循环中绕管式换热器的换热效率，最终降低氢气液化系统的液化率与比功耗，因此本文主要研究海况对烷烃类工质降膜流动与液膜的影响情况。但是常规低温液相制冷剂温度较低且组分配比复杂，采用液相混合冷剂作为实验工质时需要对换热器进行较为复杂的保冷和吹扫试压等操作，不利于观测和研究不同管型结构对降膜流动特性的影响，因此本文选用常温下为液体的工质进行替代。常温下水的密度为1007kg/m3、黏度为0.8904×10-3Pa·s、表面张力为72.10mN/m，乙醇密度为787.6kg/m3、黏度为0.9626×10-3Pa·s、表面张力为22.10mN/m，戊烷密度为620.7kg/m3、黏度为0.2204×10-3Pa·s、表面张力为15.39mN/m。与常规低温液相制冷剂相比，戊烷的物性相差较小。同时戊烷的量纲为1 的伽利略数Ga1/4（559.3）与常规低温液相制冷剂的Ga1/4（558.9）数值相近，而与乙醇（178.3）和水（497.5）的伽利略数Ga1/4则相差较大。因此优选戊烷作为降膜流动可视化实验工质，解决了液氢用预冷循环制冷剂中组分较轻不利于可视化实验，而常规的降膜流动可视化实验工质（水与乙醇）的物性与制冷剂相差较大的问题。

降膜流动浮式微观可视化实验流程与设备实物如图1 所示，实验装置包括：测试流体循环系统、高速显微摄像与图像采集处理系统、海上倾斜晃荡平台系统。其中测试流体循环系统主要由动力装置、储罐、调节阀组和降膜流动测试装置组成。测试流体由动力系统增压，经流量计、调节阀组后，流入降膜流动测试装置，在测试装置内部降膜流动后进入储罐，测试流体从储罐底部流回动力装置的入口，实现流体循环。

图1 降膜流动浮式微观可视化实验流程(a)与设备实物图(b)

降膜流动测试装置安装在海上倾斜晃荡平台上，在平台操作系统的作用下，随着平台的运动模拟海上的倾斜晃荡工况。测试装置外壳为透明亚克力材料，由于其内部流体流动速度较快，普通摄像机不能对其流动过程进行精确捕捉，因此采用高速摄像显微技术，将高速摄像机安装在海上倾斜晃荡平台上，对海况条件下测试装置内部降膜流动过程进行动态捕捉。不同的管型结构尺寸参数对管外降膜流动过程影响较大，为了更加准确地研究不同管型降膜流动流型演化机理，基于3D建模打印技术，制作了降膜流动可视化实验测试管，包括光圆管、螺纹波管、方波管、波纹管（图2）。

图2 降膜流动可视化实验测试3D打印管

1.2 可视化实验研究结果

通过可视化实验发现流体的雷诺数（Re）影响管外降膜流动流型，随着雷诺数的增加，降膜流动流型依次为滴状流流型、柱状流流型、扇状流流型及其过渡流型。相邻流型之间的Re与Ga有关，表达式如式(1)所示。滴状流流型下冷剂流量不足，管壁表面会出现干涸区域，换热效率较低，扇状流流型下制冷剂流量较大，冷剂大量浪费，导致液化系统比功耗较低，而柱状流流型下制冷剂的经济性较好。基于高速显微摄像与图像处理系统，通过微观可视化实验得到管外径8mm、管间距4mm 光圆管、螺纹波管、波纹管和方波管的降膜流动流型演化规律，如图3所示，Re与工质流量直接相关，本文所研究的Re范围为0～900。可以看出常规圆管和波纹管比螺纹波管和方波管的柱状流范围较大，流体更易形成柱状流流型。流型划分图的提出为工程中合理选择制冷剂流量范围提供理论依据。

图3 不同结构的管外降膜流动流型划分

式中，a是降膜流型划分因子；Ca是毛细数，Γ是换热管单侧单位长度降膜流动液体质量流量；μ是液体的动力黏度；σ是液体的表面张力；ρ是液体的密度；g是重力加速度。

水平静止条件下光圆管外降膜流动状态如图4所示，可以看出，在水平静止条件下降膜流动流型较为稳定，液膜均匀覆盖在管壁表面，管间液柱在液体入口中心，液柱形态稳定，液膜均匀铺展。在海况3°倾斜条件下光圆管外降膜流动状态如图5所示，实验结果表明，海况3°条件对降膜流动流型影响较小，管间液柱和液膜铺展都较为稳定，因此在海况3°条件下，常规光圆管的降膜流动流型有较好的海上适应性。

图4 水平静止条件下光圆管外降膜流动状态

图5 海况3°条件下光圆管外降膜流动状态

基于高速摄像系统和晃荡平台系统，研究了海况6°条件下光圆管外降膜流动的动态过程。如图6(a)所示，倾斜一侧流型以扇状流流型为主，流体流量过剩，而另一侧没有液膜覆盖；随着时间的推移，如图6(b)～(e)，另一侧流体流量不足，降膜流动流型以滴柱状流为主，降膜流动和液膜铺展的稳定性均较差；最后如图6(f)所示，管间降膜流动流型又恢复到图6(a)的形态。海况9°条件下光圆管外降膜流动状态如图7 所示，可以看出，管间降膜流动状态与海况6°条件下流动状态相似，管间液体偏移更为明显，液滴快速倾斜降落，在倾斜侧聚并形成液扇。这是因为海况主要是对重力起主导作用的降膜流动过程影响较大，倾斜海况下，重力的轴向分力会产生轴向加速度，从而影响液膜沿轴向的铺展过程，导致管间降膜流型多以局部扇状流和滴状流为主，因此在海况6°和9°条件下常规光圆管的海上适应性较差。

图6 海况6°条件下光圆管外降膜流动状态

图7 海况9°条件下光圆管外降膜流动状态

基于海上倾斜晃荡平台、高速摄像和3D 打印技术，分别研究了螺纹波管、方波管和波纹管的降膜流动海上适应性。水平静止条件下方波管外降膜流动状态如图8 所示，可以看出，液柱多从方波管的大管径处降落，液柱的流动过程较为稳定，液膜铺展性也较好。但当处于倾斜9°海况时，如图9(a)～(f)所示，由于大管径和小管径处突变的影响，液体无法在轴向均匀铺展，出现了大量的干涸区域，因此方波管外降膜流动的海上适应性较差。水平静止0°条件下螺纹波管外降膜流动状态如图10 所示。从图10(a)、(b)、(f)可以看出，流体多从螺纹波管的螺纹缝隙降落，会出现少量的干涸区域。当螺纹波管处于倾斜6°海况时，如图11 所示，与常规光圆管相似，在倾斜一侧以扇状流流型为主，而另一侧以滴状流流型为主，受螺纹缝隙阻隔的影响，滴状流流型下管壁表面出现了大量的干涸区域，海上适应性较差[图11(c)～(f)]。因此通过可视化实验发现，在海况条件下，由于管壁表面突变导致了较强的阻隔效应，螺纹波管和方波管外降膜流动和液膜铺展的均匀性、稳定性均较差。

图8 水平静止条件下方波管外降膜流动状态

图9 海况9°条件下方波管外降膜流动状态

图10 水平静止条件下螺纹波管外降膜流动状态

图11 海况6°条件下螺纹波管外降膜流动状态

图12 给出了水平静止0°条件下波纹管外降膜流动过程，可以看出，液柱多从波纹管的突出部分降落，液柱较为稳定，液膜铺展性也较好。当处于倾斜6°海况时，如图13 所示，管壁表面仍未出现干涸区域，管间流型多以柱状流流型为主，虽然液柱倾斜，但管间流型依然稳定，液膜的均匀稳定性较好。当海上倾斜角度增加到9°时，如图14所示，与6°海况类似，波纹管管间液柱虽有少量倾斜，但整体流型较为稳定，且液膜铺展的均匀性和稳定性均较好。因此在海况条件下，从管间流型稳定性、液膜分布均匀稳定性、管壁干涸情况、液膜纵向偏移和柱状流流型区间等多方面考虑，相较于常规光圆管、螺纹波管和方波管，推荐采用波纹管作为FLH2中降膜流动换热段的管型结构，以强化浮式氢气液化装置的海上适应性。

图12 水平静止条件下波纹管外降膜流动状态

图13 海况6°条件下波纹管外降膜流动状态

图14 海况9°条件下波纹管外降膜流动状态

2 波纹管降膜流动数值模拟模型建立

2.1 控制方程与计算方法

为了进一步定量研究波纹管的管外降膜流动特性，得到薄液膜厚度分布等微观尺度参数，本文建立了波纹管降膜流动数值模型。模型的连续性方程如式(2)。

动量方程如式(3)～式(6)所示。

F为动量源项，如式(7)所示，包括重力源项（FG）、表面张力源项（Fσ）、气液拖曳力源项（FLG），如式(8)、式(9)。

本文采用耦合的VOF[volume of fluid，如式(10)]和Level Set两相流模型进行降膜流动数值模拟。数值模拟模型采用压力速度耦合方法选择PISO，动量方程空间离散采用精度较高的QUICK，压力离散选用PRESTO!，体积分数选用Geo-Reconstruct，Level Set采用QUICK。

2.2 网格划分及模型验证

波纹管管外降膜流动数值模型的网格划分如图15 所示，圆周角定义为α。在管壁周围区域进行网格加密，贴近管壁处网格尺寸为0.015mm，以捕捉管壁附近的薄液膜等微观尺寸参数。模型为质量较高的六面体结构化网格，时间步长设置为1×10-5s。入口边界条件设置为质量流量，左右两侧边界为对称边界条件（symmetry）。为了验证波纹管管外降膜流动数值模拟结果的可靠性，模拟得到的戊烷工质（密度620.7kg/m3、黏度0.2204×10-3Pa·s、表面张力15.39mN/m）降膜流动流型与可视化实验结果进行了对比，如图16所示，通过对比发现，构建的降膜流动数值模拟模型与微观可视化实验结果相似性较好，可以准确模拟管外降膜流动流型。

图15 网格划分示意图

图16 降膜流动流型模拟结果与实验结果对比

3 数值模拟结果

基于动网格技术，编写了UDF（user defined function）程序导入数值模拟模型中，研究不同角度的海况条件对降膜流动过程的影响规律。在倾斜6°海况条件下，常规光圆管外降膜流动随时间的变化情况如图17 所示。在48.581ms 时，第一根管壁表面开始覆盖液膜，液膜逐渐铺展，到65.581ms时，第一根管壁表面完全覆盖，液柱开始逐渐形成，随后在100.301ms时，液柱逐渐向下拉伸，与第二根管壁表面接触，但在倾斜海况和重力的作用下，倾斜一侧液柱先降落，重力的轴向加速度导致倾斜侧管壁表面液膜较厚，另一侧管壁表面液膜铺展性较差。在174.791ms 时，出现了大量干涸区域，干涸区域面积随着时间的增加虽有所减小，但仍然存在（t=210.110～252.681ms）。

图17 海况6°条件下光圆管外降膜流动随时间变化过程

水平静止0°条件下波纹管外降膜流动随时间变化过程如图18 所示。在25.086ms 时，液体开始降落，在35.086ms 时，降膜流动液体逐渐在轴向汇聚形成波峰，随后59.876ms 时，第一根管壁表面已经完全铺展液膜，并在突出部分汇聚液体。在103.016ms 时，第二根管壁表面液膜覆盖约一半，到123.016ms，液膜几乎完全覆盖波纹管，第二根管只有底部少量管壁表面出现干涸区域。随着时间的推移，到149.341ms时，液膜完全铺展覆盖在波纹管壁表面。

图18 水平静止条件下波纹管外降膜流动随时间变化过程

图19 给出了海况6°条件下波纹管外降膜流动随时间变化过程，在31.303ms时，液膜开始覆盖管壁，逐渐沿着轴向和周向铺展，到51.303ms时，液膜完全覆盖第一根波纹管壁表面。在71.303ms时，液体以液柱的形式开始降落，随后到91.303ms，液柱逐渐与第二根管壁表面相接触，在重力轴向加速度的作用下，倾斜一侧液体较多，覆盖的区域也较宽。在117.880ms时，倾斜一侧管壁几乎完全覆盖，而另一侧管壁还出现少量干涸区域，但在183.223ms时，波纹管壁表面被液膜完全覆盖。

图19 海况6°条件下波纹管外降膜流动随时间变化过程

海况9°条件下波纹管外降膜流动随时间变化过程如图20 所示。在20.701ms 时，液体降落到第一根管壁表面，在65.712ms时，液体逐渐向轴向铺展，完全覆盖第一根波纹管并在突出部位汇聚。随后108.652ms时，与海况6°时相似，管间液体在海况的影响下逐渐偏移。在138.652ms时，液膜几乎覆盖了第二根波纹管的表面，但在倾斜另一侧仍旧有少量的干涸区域，到152.852ms，干涸区域逐渐缩小。随着时间的推移，在251.412ms时，液膜完全覆盖第二根波纹管管壁表面。通过图19 和图20研究发现，在海况条件下，波纹管管间液柱虽有少量偏移，但降膜流动流型总体稳定，液膜也较为均匀地覆盖在管壁表面，因此波纹管有较好的海上适应性。

图20 海况9°条件下波纹管外降膜流动随时间变化过程

为了与水平静止0°时的降膜液膜厚度进行定量对比，采用齐次坐标技术来描述空间的各点坐标及其变换，将通过动网格技术倾斜的降膜流动数值模型恢复到水平静止0°状态，转换方法如式(11)～式(14)所示。

通过式(11)～式(14)的计算分析，得到了水平静止0°、海况6°和9°倾斜条件下，波纹管外降膜流动液膜厚度沿纵向管长分布。如图21 所示，在圆周角为45°时，水平静止、海况6°和9°的主波峰厚度分别为0.482mm、0.519mm、0.546mm，海况6°和9°的主波峰纵向位置分别偏移了1.237mm 和1.425mm。圆周角增加到80°时，水平静止0°、海况6°和9°的主波峰厚度分别为0.532mm、0.557mm、0.535mm，主波峰的纵向位置分别偏移了1.328mm和1.895mm。进一步增加圆周角至110°时，倾斜海况6°的波峰位置偏移了1.565mm，倾斜海况9°的波峰位置偏移了1.853mm。圆周角为145°时，水平静止0°、海况6°和9°的主波峰厚度分别为0.688mm、0.701mm、0.766mm，倾斜6°和9°海况的主波峰位置分别偏移了1.140mm 和1.378mm。倾斜海况导致了波峰沿纵向管长的偏移，倾斜角度越大，偏移量越大；但随着圆周角增加到145°时，波峰纵向位置偏移量逐渐减小，可以看出，随着流体沿波纹管降膜铺展，波纹管的管型结构对倾斜流体的修正作用不断增强，偏移程度逐渐降低，进一步验证了波纹管较好的海上适应性。

基于数值模拟结果，通过修正Hou 等[34-35]提出的液膜厚度计算关联式(15)，得到波纹管外降膜流动平均液膜厚度的计算关联式。

式中，α为圆周角；ρG和ρL为气体和液体的密度；μL为液体的黏度；c和n通过研究结果进行拟合得到；S为管间距；管径D=( )Dmin+Dmax2。本文根据雷诺数Re和管径D范围划分计算关联式的适用范围：低（Re/D）<70mm-1，中（Re/D）为70mm-1≤（Re/D） <170mm-1， 高（Re/D） ≥170mm-1。 在 低（Re/D）的条件下，c和n值与Hou模型的数值相同。在中（Re/D）条件下，圆周角<90°时，c为1.15，n为-0.07；圆周角>90°时，c和n值与Hou模型的数值一致。在高（Re/D）条件下，圆周角<90°时，c和n值分别为1.96和-0.07；圆周角>90°时，c和n值分别为1.60和0.07。优化后波纹管平均液膜厚度关联式计算值与模拟值的误差分析如图22所示，可以看出误差都在±25%以内，大部分误差在±20%以内。图23给出了海况下平均液膜厚度计算关联式的精度，计算关联式也可以较好地预测海况下波纹管管外降膜流动的平均液膜厚度。平均液膜厚度计算关联式的提出为工程中合理选择波纹管型结构和流体雷诺数提供理论依据。

图22 优化后理论计算关联式的波纹管液膜厚度与模拟值对比

图23 海况下理论计算关联式的波纹管液膜厚度与模拟值对比

4 结论

基于高速摄像系统、晃荡平台系统和3D 建模打印技术，搭建了研究常规光圆管、螺纹波管、方波管和波纹管的降膜流动浮式可视化实验装置，进行了海况条件下降膜流动的微观可视化实验研究；基于耦合的VOF和Level Set两相流模型以及动网格技术，进行FLH2中降膜流动的数值模拟研究，揭示降膜流动海上适应性的强化机理，主要结论如下。

（1）构建的降膜流动数值模拟模型与微观可视化实验结果吻合性较好，可以准确模拟管外降膜流动流型；绘制了常规光圆管、螺纹波管、波纹管和方波管管外降膜流动的流型划分图。

（2）通过降膜流动浮式微观可视化实验和数值模拟研究结果可知，在海况条件下，从管间流型稳定性、液膜分布均匀稳定性、管壁干涸情况、液膜纵向偏移和柱状流流型区间等多方面考虑，相较于常规光圆管、螺纹波管和方波管，推荐采用波纹管作为FLH2中降膜流动换热段的管型结构，以强化浮式氢气液化装置的海上适应性。

（3）通过实验和数值模拟相结合的方法，得到了波纹管降膜流动平均液膜厚度计算关联式，根据雷诺数和管径比值划分计算模型的适用范围。基于Hou理论模型，拟合了不同雷诺数Re和管径D下的波纹管外降膜流动的平均液膜厚度计算模型，误差都在±25%以内，海况条件下液膜厚度计算关联式也有较高的计算精度。