基于CFD的独立C型低温液舱蒸发特性分析

，

(大连理工大学 船舶工程学院，辽宁 大连 116024)

对于LNG运输船上的低温液舱而言，外部热量会通过绝热维护系统进入液舱内部，会导致低温液体蒸发，造成舱内压力升高及液货损失。目前，国内大多数的研究是从热力学角度通过解析公式对低温液舱蒸发率进行了计算[1-3]，但没有直观揭示低温液舱温度、压力及流场的分布规律；或者是利用公式计算了壁面热流密度，并利用二维模型研究了储罐内部流场、温度场以及压力场的分布规律[4]。国外有学者通过4种方法从二维模型到三维模型数值计算了液舱绝热维护系统的温度场[5]，但是没有研究液舱内部流体相变及其流场的分布规律。为此，考虑采用三维模型，将绝热维护系统及液舱内部低温流体蒸发耦合进行计算，得到整个系统的流场、温度场、压力场等分布规律。

1 数值模型和数值方法

1.1 VOF方法

采用VOF方法来追踪相间界面，求解一套守恒方程并通过求解控制单元内单相或多相体积分数的连续性方程来确定相间界面，对于第q相来说，方程如下。

▽·(αqρqVq)=

(1)

方程(1)中主项的体积分数求解式为

(2)

1.2 蒸发模型

在液舱内部饱和状态下，克拉佩龙方程描述了压强和温度的相互关系。

(3)

式中：ρg、ρl分别表示气、液相的密度。一般情况下，饱和状态在交界面附近，而气相部分的温度高于饱和温度，此时压强的变化率与气相温度之间的关系如下[6]。

(4)

式中：R为气体常数;V为气相体积;Pi、Tgi分别为初始压强和初始温度。

壁面漏热量会引起液舱内部流体能量的变化，同时引起传热传质过程，在气液交界面处的传热传质过程可以描述为[7]

(5)

选择Lee模型来确定蒸发过程，模型中气相和液相的质量转移由下面的输运方程控制。

(6)

在引入饱和温度时，质量转移用下式描述。

Tl>Tsat时，

(7)

Tg

(8)

1.3 控制方程

1)质量守恒方程。

▽·(ρV)=Sm

(9)

式中：V是速度矢量；Sm为源相。

2)动量守恒方程。

▽·(ρVV)=

-▽p+▽·(τ)+ρg+F

(10)

式中：g为用在微元体上的重力；F为其他外部体积力，同时还包括其他模型的相关源相;τ为因分子粘性作用而产生的作用在微元体表面上的粘性应力张量。

1.4 能量守恒方程

式中：keff为有效导热系数；hj为j组分的焓；Jj为j组分的扩散通量；τeff为粘性应力张量；Sh为能量源相。

2 低温物理模型实验

1.5 m3卧式独立C型低温液舱见图1，液舱尺寸见表1，绝热维护系统所用材料及其属性见表2。

图1中①、②、③、④处分别安装温度传感器，在管路处安装压力传感器。

实验装置包括实验液舱及支撑装置、管路系统及控制阀件、温度传感器、压力传感器、压力表，数据采集系统、地磅、空调、温度湿度计、电脑、梯子、千斤顶等辅助工具。实验介质为液氮(标准大气压下，沸点：-195.8 ℃)。实验工况为充装率49.13%。

表1 液舱尺寸 mm

表2 材料属性

实验流程。在实验装置准备完毕后，进行仪表标定调试，并设置空调参数保证室内温度，然后开始对液舱进行预冷，同时检查阀门、仪表及采集系统是否正常工作。预冷工作就绪后，进行加注，然后关闭所有阀门憋压静置。当整个系统温度场达到相对稳定，打开气相管子出口阀门泄压，保持阀门一直开放，进行0表压(即基本与大气压强相等)蒸发过程的监测实验。实验过程中用温度、压力传感器进行温度、压力监测，并实时采集重量数据。

处理实验数据，提取充装率49.13%工况下，①、②、③、④处温度传感器，以及管路处压力传感器对应的测量值见表3。

表3 充装率49.13%实验结果

3 数值计算

选用Fluent软件进行数值模拟，计算分为两个阶段。

1)固体域传热计算。只考虑固体热传导，对绝热层、鞍座、内壳及管路系统进行计算，得到相对稳定的温度场，然后输出各处所对应的温度值。

2)整个系统耦合计算。整理第1阶段得到的数据，并编写初始化UDF，以初始条件的方式将第1阶段的计算结果加载到计算模型固体域上；液舱内部流体，分别对气相、气液界面附近、液相3个区域给定不同的初始温度，将绝热维护系统以及低温液舱内的流体耦合起来进行计算。

3.1 建模及网格划分

计算模型具有对称性，因此采用四分之一模型进行数值模拟，见图2，计算区域包括内部流体域、内壳及管路系统、绝热层、鞍座。计算模型坐标原点位于横向半圆切面的圆心处，x轴沿模型纵向方向，y轴沿模型竖直方向，z轴沿模型横向方向。采用ICEM软件划分网格，所有区域均采用结构网格，整体网格见图3；对内壳及管路、鞍座进行网格的细化，以满足计算要求，见图4。

3.2 数值模拟设置

第1阶段。由文献可知，独立C型低温液舱内部气相的温度高于气液界面处的温度，气液界面处的温度高于液相的温度。因此，在给定初始条件时，分别在内壳及管路系统与气相、气液界面处、液相接触处设置不同的温度，具体见表4。绝热层与空气接触表面设置为293.15 K的固定温度壁面条件。

表4 初始温度设置(充装率49.13%) K

第2阶段。将第1阶段所得到的数据以初始化UDF的形式加载到整体模型算例中，液舱内气相、气液界面处、液相初始温度设置与表4一致。初始压强给定为0.2 MPa，舱内部流体采用液氮及氮气作为模拟工质，液氮在标准大气压下的物理属性见表5。

表5 标准大气压下液氮物理属性

采用压力基瞬态计算，采用VOF多相流模型，并包含Implicit Body Force，采用标准k-ε双方程湍流模型，以及标准壁面函数，同时计算Viscous Heating。选择Lee模型作为蒸发冷凝模型，并设置饱和温度为77.35 K，表面张力选择Continuum Surface Stress模型，并设置表面张力为0.008 59 N/m。绝热层与空气接触壁面设置为293.15 K的固定温度壁面条件，其他区域接触面设置为coupled热边界条件，内壳、管路系统与流体接触面设置为无滑移速度条件。管子出口设置为pressure outlet。压力速度耦合采用coupled算法。

4 结果分析

图1中①～④处传感器的位置(相对于数值模型坐标系)见表6。

表6 传感器位置坐标 m

提取位置①～④处数值计算结果以及管路内的压强与实验值比较，见表7。

表7 实验与数值计算结果对比

由表7可以看出，数值计算结果与实验结果相对误差均小于3%，进而也说明了数值计算方法的可靠性。

图5a)为设置完初始温度并加载完UDF后的温度场，可以看出初始化UDF已将温度参数加载到整体模型上。图5b)可以看出储罐内部流体存在着温度梯度，整体表现是自下而上，温度逐渐升高，对比图5a)可以说明，气相和液相通过气液界面存在传热现象，且气相内部也存在着热量交换，这也造成了初始 “过热”气体温度逐渐降低，且越靠近气液界面处气相温度降低越明显。

图6显示了管路处温度分布，由图6a)可以看出管路的存在使得周围的温度明显低于同一高度其他各处的温度。这是由于蒸发气体不断地通过管路排出，而气体相对于绝热层及管路来说温度更低，且管路的导热系数相对较大，管路和气相接触，导致管路的温度迅速降低，进而影响了管路附近区域温度分布。由图6b)可以看出，在管路左侧0.02 m处(图6a)中黑线)，由于低温储罐及管路的影响，沿竖直方向在0.51～0.83 m范围内，温度由104 K变化到111 K,温度变化较小，而从0.83～0.91 m范围内，温度由111 K变化到225 K,温度变化明显。

另外由于鞍座材料的导热系数相对较大，所以鞍座温度也会很快降低，进而影响其周围温度分布，见图7。

由图7可见，鞍座的存在使得其附近区域温度低于较远区域的温度，且鞍座处温度最低，大约为210 K，在距鞍座0.25 m范围内最大温差为30 K。

图8显示了不同切面处温度分布，从图8a)可以直观地看出液舱内部温度整体趋势为上高下低，但从图8c)、图8d)可以看出在底部近壁面处，温度稍高于其上部温度。从图8b)、图8c)、图8d)可以看出液体在壁面处的温度都高于同一水平内部流体温度，这是因为外部传热使壁面处流体温度升高。

由图9可以看出，在近壁面处，液相温度由77.35 K逐步增加到近壁面处的81.5 K。由图8d)可以看出，在出口处气体温度有所升高，这是由于管路在此处直接与外界接触，且气体也直接与外界接触，传热过程加剧所致。

从图10可以具体地看出，在液相整体温度相对均匀，温度梯度非常小，而对于气液界面附近及气相区域温度梯度相对较大。

图11显示出了内壳及管路系统温度场。

由图11a)可以看出内壳与管路接触处温度明显低于其周边温度，这是由于低温气体在此处排出，会带走部分热量，致使其温度降低；图11b)显示了内壳与气液界面接触处，由液相温度过渡到气相温度的温度梯度；由图11c)可以看出，内壳与鞍座接触处温度高于周边温度，但温度差异不明显。

图12显示了气体的体积分数。

可以看出自气液界面向下，气体的体积分数在逐渐减小，且在与壁面接触处气体的体积分数并无显著变化，这说明蒸发现象主要发生在气液界面附近。这是由于壁面附近的边界层流对流循环会将热量携带到气液界面处，再加上气液界面处的壁面漏热，以及气相与液相之间的传热，使得气液界面处液体温度升高相对较大。因此，蒸发现象明显，而在壁面附近由于液体及固体壁面温差相对来说并不大，漏热量并不显著，同时壁面边界层的对流循环也会将壁面漏入的热量带走，因而蒸发并不显著。此外，分析可以看出，在正常的低温液体存储下，流体与固体壁面之间的温差相对较小，此状态时的漏热量并不会引起低温液体的沸腾，而传质过程以表面蒸发的方式发生。

压强分布见图13。

由图13可见，当处于自由蒸发状态时，液舱内压强上低下高，且管路内部压强基本处于大气压强，而气相部分压强稍高于大气压强。由于液体压强的存在导致液舱底部压强有明显的增大，相对气相来说压强增大4013.3 Pa，这在结构设计时应予以注意。

整体速度矢量见图14(管路内气体速度明显高于其他部分，为显示出速度分布差异，对显示范围作了相应的调整得到图14)，可以看出气相流体的速度高于液相流体速度，且壁面附近流体速度高于内部流体速度，气体的比热容比液体小，所以传入相同的热量，气相温度升高的更快，温差增大，因而气相流体流动更剧烈。

由图15和图16速度矢量图可以看出壁面处的流体沿壁面向上运动，这是由于壁面处的流体随着热量的不断流入，温度升高，形成热边界层，且在浮升力作用下沿着壁面向上流动。

流线见图17，局部放大见图18。

由图17、18可见，随着外界热量的进入，在壁面附近形成的边界层流沿壁面向上流动，在气液界面处向中心流动，在此过程中一部分由于蒸发进入到气相，一部分滞留在气液界面附近区域，在中心区域形成中心射流向下运动，继续参与液体内部的对流循环运动。液体的这一对流运动过程，将热量不断地带入自由液面区域，使得自由液面处的液体流速增大，对流强度增大。同理，在气相也会发生类似的对流循环过程，沿壁面运动的边界层流一部分从管路排出，一部分向下运动，形成回流区域。图19所示的整体流线图，也表明了在气液界面处及气相靠近壁面处，对流强度较大，换热明显。

5 结论

1)对于固体域温度分布而言，鞍座、管路与液舱内壳接触处，温度分布会有变化，这在设计建造时应予以重视；对于液舱内部的流体而言，整体温度呈下低上高分布，同时壁面附近处温度高于内部区域温度，且相对液相来说，气液界面处以及气相区域温度梯度较大。

2)相对壁面附近而言，蒸发现象在气液界面处更加明显。

3)自由蒸发时，液舱内部压强呈现上低下高分布，气相区域压强基本均匀，且稍高于大气压强，而液相由于液体压强的存在，随着深度的增加，压强逐渐增大。

4)气体的速度高于液体，且在气相近壁面处及管路内尤为明显，液相沿壁面向上流动的边界层流一部分蒸发进入气相，一部分沿界面向中心流动且在界面下方形成回流区域。

5)液舱内部气相近壁面处及气液界面附近自然对流现象相对于其他区域更为明显。

相对于目前的研究而言，本文所述将低温液舱绝热维护系统及舱内低温流体相变耦合进行计算的方法更符合实际情况，更便于分析整个系统的相应规律，但随着低温液舱的大型化，前期初始条件加载需要处理的温度数据会显著增多，在进一步应用时可以考虑将数据分组回归出函数表达式的形式进行加载。本研究所得到的蒸发特性规律，对低温液舱的设计建造具有重要的参考价值，同时在得到整个系统的温度场分布时，可以进一步提取温度数据，对结构的热应力进行深入分析，从而更好地优化设计。

[1] 刘文华,陆晟.中小型LNG船C型独立液货舱蒸发率计算[J].船舶设计通讯,2012(1):25-28.

[2] 时光志,盛苏建.独立C型液货舱的传热分析及蒸发率计算[J].船海工程,2013,42(1):65-69.

[3] 田竹刚,尹奇志,王冰.船用LNG燃料储罐蒸发量与维持时间计算方法[J].交通信息与安全,2016,34(3):57-63.

[4] 姚寿广,闫兴武,叶勇,等.某型船用LNG储罐热响应仿真研究[J].江苏科技大学学报(自然科学版),2016,30(2):136-141.

[5] MIANA M, LEGORBURO R, DEZ D, et al. Calculation of boil-off rate of liquefied natural gas in mark III tanks of ship carriers by numerical analysis[J]. Applied Thermal Engineering,2016,93:279-296.

[6] HOPFINGER E J, DAS S P. Mass transfer enhancement by capillary waves at a liquid vapour interface[J]. Experiments in Fluids,2009,46(4):597-605.

[7] GROTLE E J, et al. Non-isothermal sloshing in marine liquefied natural gas fuel tanks[C]. The Twenty-sixth International Ocean and Polar Engineering Conference,2016:879-885.