储罐中心区填充多孔材料对LNG流动特征影响的数值模拟

邢志祥，张淑淑，汪李金，张 莹

(常州大学环境与安全工程学院，江苏 常州 213164)

多孔介质是一种多相物质共存、由固相和气相或液相组成的物质，固相作为固体骨架，气相或液相充满了孔隙空间，孔隙之间必须相互连通[1]。Nield等[2]最先开始了多孔介质的研究，认为多孔材料具有比表面积大、质量轻、体积小、导热性能好的特点。近年来，针对多孔介质的试验和模拟研究很多，研究主要集中在填充多孔材料后液体流动的情况以及多孔材料内部火焰传播的情况，而研究对象通常是管道、密闭容器等，但关于储罐内填充多孔材料的研究则相对较少。邢志祥等[3]模拟了圆柱形储罐内填充聚氨酯多孔材料后可燃气体火焰传播的情况，并与13 L圆柱形密闭储罐现场试验的结果进行了对比，结果表明：模拟结果与试验结果相吻合，填充物平均孔径越小，阻火抑爆性能越好；张健中等[4]分析了在加油站埋地油罐中是否有必要使用网状铝合金阻隔防爆材料等问题;田宏等[5-6]对填充在液化石油气储罐中的多孔金属材料做了系统的介绍，但并未通过试验来举证。

液化天然气(Liquified Natural Gas,LNG)储罐是储存液化天然气的重要设备，其运行处于超低温状态时容易引发各类安全事故，而LNG翻滚是其中一种事故类型。储罐内由于LNG密度差产生分层，上层液体与下层液体之间存在液液分界面，而漏热导致分界面被破坏进一步引发翻滚现象。分界面被破坏的原因有两种[7]：一种是储罐内壁的边界层穿透，导致分界面被破坏；另一种是上层液体在中心处向下的射流和下层液体在中心处向下的卷携，导致中心处的分界面被破坏，进一步造成了整个分界面被破坏。考虑到多孔材料的优点，本文将其应用到储罐中研究其是否能起到抑制LNG翻滚的作用，并选取储罐中心区填充多孔材料，对储罐内LNG的流动特征进行了数值模拟。

1 填充多孔材料的储罐内液体流动计算模型的建立

本文选取储罐中心区填充多孔材料，并分别设置多孔材料的填充厚度为1 m、1.5 m、2 m，通过建立储罐中心区填充多孔材料的物理模型和数学模型，对储罐内LNG的流动情况进行了数值模拟计算。

1. 1 物理模型

本文选取6 000 m3的圆柱形储罐，内罐直径为24 m，高度为17.4 m，设计液位高15.76 m，并取储罐的轴截面建立多孔介质中湍流流动的二维模型。坐标原点位于左侧壁面和底部壁面的交汇处，X轴向右为正，Y轴向上为正；模拟过程不考虑气相空间，只针对储罐内的液相部分来模拟LNG的流动情况，液体充装高度为14 m；多孔材料填充在储罐的中心部位；对储罐结构进行了简化，不考虑壁厚，并忽略储罐内部的管线、循环泵等部件，由此建立的储罐中心区填充多孔材料简化后的物理模型如图1所示。

图1 储罐中心区填充多孔材料简化后的物理模型Fig.1 Simplified physical model of the filling porous material of the storage tank central area

对储罐中心区填充多孔材料时储罐内LNG流动情况的物理模型(中间标明H的区域为多孔材料区域，左右两侧区域为纯流体区域)进行了如下的简化：①储罐内部流体分成两个高度相等的分层，分层高度为7 m，上分层流体密度为424 kg/m3，下分层流体密度为425 kg/m3，初始密度差为1 kg/m3，两个分层内的密度均匀一致；②储罐内流体密度符合Boussinesq假设；③选取金属铝多孔材料，且视多孔材料为均匀各向同性；④多孔材料的孔隙为球形空洞，孔隙之间相互连通，即为开孔结构；⑤多孔材料的物性参数为常数；⑥多孔材料固体骨架与流体处于局部热平衡，且无任何化学反应；⑦储罐内流体流动为湍流流动；⑧忽略流体流动过程中的黏性耗散；⑨不考虑辐射换热。

1. 2 数学模型

采用非稳态模型、VOF模型和标准k-ε湍流模型，并基于Boussinesq假设，计算过程中多孔区域动量方程采用Darcy-Brinkman-Forchheimer模型，能量方程采用局部热平衡模型，得到的控制方程如下：

连续性方程:

∂u∂x+∂v∂y=0

(1)

动量方程:

ρφ2∂u∂t+u∂u∂x+v∂u∂y=-∂p∂x+

μφ∂2u∂x2+∂2u∂y2-μKu+ρCFKu2+v2u

(2)

ρφ2∂v∂t+u∂v∂x+v∂v∂y=-∂p∂y+

μφ∂2v∂x2+∂2v∂y2-μKv+ρCFKu2+v2v-ρg

(3)

能量方程:

ρCp∂T∂t+u∂T∂x+v∂T∂y=λeff∂2T∂x2+∂2T∂y2

(4)

k方程：

∂∂t(ρk)+∂∂xi(ρkui)=∂∂xjμ+μtσk∂k∂xj〗+

Gk+Gb-ρε-YM+Sk

(5)

ε方程:

∂∂t(ρε)+∂∂xi(ρεui)=∂∂xjμ+μtσε∂ε∂xj〗+

C1εεk(Gk+C3εGb)-C2ερε2k+Sε

(6)

上式中：u为X轴方向的速度分量(m/s)；v为Y轴方向的速度分量(m/s)；ρ为液体密度(kg/m3)；t为时间(s)；p为压力(Pa)；μ为动力黏度(Pa·s)；g为重力加速度(9.81 m2/s)；φ为孔隙率；K为渗透率(m2)；CF为Forchheimer系数；Cp为定压比热容[J/(kg·K)]；T为温度(K)；λeff为多孔区域的有效导热系数[W/(m·K)],由流体的导热系数λf和固体的导热系数λs的体积平均值计算得到，即λeff=(1-φ)λs+φλf；k为湍动能(J)；ε为湍流扩散率；σk和σε分别为k方程和ε方程的普朗特数，σk=1.0，σε=1.3；μt为湍流黏性系数，μt=Cμρk2/ε，其中Cμ=0.09；Gk为平均速度梯度引起的湍动能k的产生项；Gb为浮力引起的湍动能k的产生项；YM为可压缩湍流流动中脉动扩张的贡献；C1ε、C2ε、C3ε为经验常数，C1ε=1.44，C2ε=1.92，C3ε=1；Sk和Sε为用户自定义源项。

当φ=1时，k→∞，λeff=λf，动量方程的源项趋近于0，能量方程也变为标准能量方程，则上述控制方程变为纯流体区域的控制方程。

1. 3 边界条件和初始条件

储罐侧壁面和底壁面取为无滑移壁面边界条件，并采用定热流的方式加热，热流密度恒定为30 W/m2；气液交界面取为压力出口；多孔区域采用金属铝多孔材料，孔隙率设定为0.95，并需设定黏性阻力系数和惯性阻力系数；多孔区域与纯流体区域的交界面设为内部边界。

初始时刻：初始表压为15 000 Pa，X和Y轴方向的速度分量都为0，分层区初始温度为111 K，主流区初始温度为111.5 K。

在计算过程中，多孔区域和纯流体区域通过φ的取值进行区分，并统一进行求解。

1. 4 相关参数介绍

(1) 孔隙率φ和孔径dp[8]：孔隙率φ是指孔隙体积占多孔材料总体积的比值；孔密度是指多孔材料每英寸上孔的数目(1 inch=0.025 4 m)，孔径dp=0.025 4/孔密度(m)。

(2) 渗透率K和Forchheimer系数CF[9]：渗透率K代表了多孔介质中孔隙的表面积和弯曲程度，通常由Vafai总结的经验公式确定，即

当多孔材料中的流体流动为湍流流动时，达西定律便不能描述流体流动的情况，此时可采用其修正模型进行描述，即在动量方程中加入源项，源项由两项组成，一项为黏性阻力项，一项为惯性阻力项，惯性阻力项的系数用CF表示，称为Forchheimer系数，可表示为

CF=1.75150φ3/2

(3) 压差Δp：压差指储罐内上下分层之间的压力差(Pa)，Δp=p2-p1，其中p2为底壁面平均压力(Pa)，p1为压力出口平均压力(Pa)。

2 储罐内LNG流动特征的数值模拟与分析

本文选取储罐的轴截面建立二维模型，应用Fluent 6.3软件对储罐中心区无多孔材料填充和有多孔材料填充时储罐内LNG的流动进行了数值模拟，模拟断面均取自储罐的同一轴截面，并选取速度云图和流线图的模拟结果进行了分析。速度云图是描述流体速度大小分布情况的图；流线图是描述某一时刻流体运动趋势的图，可以通过流线图中的流线方向和流线疏密，判断出某一时刻流体的运动方向和区分流动强弱区域，流线图中箭头表示流体的运动方向，流线的疏密反映流速大小，流线越密，流体流速越大，反之亦然。

2.1 储罐中心区无多孔材料填充时储罐内LNG的流动情况

图2为储罐中心区无多孔材料填充时,t=100 s和t=1 000 s时刻储罐内LNG的流动情况模拟结果。

图2 储罐中心区无多孔材料填充时t=100 s和t=1 000 s 时刻的储罐内流体的速度云图和流线图Fig.2 Velocity cloud and stream traces of fluid in the storage tank with no filling porous materials at t=100 s and t=1 000 s

由图2可以看出：

(1) 在t=100 s时，由于外界环境的漏热，储罐壁面处的流体率先吸收热量导致密度减小，产生浮力驱动流，沿着壁面向上移动，形成流动边界层；由于上下两层流体存在温度差和密度差，因此位于分界面处的流体出现自然对流现象，且流体流速较储罐其他位置的流速大，而分界面处边界层的流体流速最大[见图2(a)]。

(2) 在t=1 000 s时，分界面处流体的自然对流发展到整个储罐内流体的自然对流，且上下两层流体的自然对流各自独立，即由液液分界面隔开，被限制在自身区域内。对于上层流体，由于储罐侧壁面的漏热，导致边界层处的流体受热膨胀密度减小，沿储罐壁面向上移动，当其到达气液交界面时，其中一部分流体蒸发变成气体释放到气相空间，其密度增大导致向下移动，由此形成了上层流体的自然能对流现象；对于下层流体，边界层处的流体受热膨胀密度减小，沿着边界层向上移动，当其到达液液分界面时，因浮力太小而无法穿过分界面进入上层，但会通过分界面向上层流体传递热量，因而温度降低，密度增大，导致向下移动，由此形成了下层流体的自然对流现象[见图2(b)]。

(3) 在纯自然对流运动中，流体的瑞利数Ra是判定由浮力产生的对流强度大小的标准[10]。当Ra小于临界值时，流体之间是热传导状态，不发生对流运动；当Ra大于临界值时，才会发生对流运动，从而在液体中出现宏观对流花纹，称其为Benard花纹[11]。储罐内无多孔材料填充时，在t=1 000 s时，上层和下层的自然对流结构明显不同，但流动都是由类似圆形的滚动圈组成，上层的滚动圈数为1个，为扁长形，而下层的滚动圈数为4个，分别为2个胖圆形和2个扁长形的滚动圈；相比于t=100 s时分界面处的流体流速较高，t=1 000 s时分界面处的流体流速反而较小，除此之外，在滚动圈中心处的流速也较小，称之为滞留区[见图2(c)]。

在模拟计算过程中，对二维截面上流体的最大流速进行了监测，得到流体的最大流速为0.195 m/s。

2.2 储罐中心区多孔材料填充厚度为1 m时储罐内LNG的流动情况

图3为储罐中心区多孔材料填充厚度为1 m时，t=100 s和t=1 000 s时刻储罐内LNG的流动情况模拟结果。

图3 储罐中心区多孔材料填充厚度为1 m时t=100 s和t=1 000 s时刻的储罐内流体的速度云图和流线图Fig.3 Velocity cloud and stream traces of fluid in the storage tank with the filling thickness of porous materials being 1 m at t=100 s and t=1 000 s

由图3可以看出：

(1) 在t=100 s时，储罐内流体流动发展情况相对于无多孔材料填充时比较滞后，分界面处的流体流速较无多孔材料填充时小，除边界层外，多孔材料填充区域两侧自然对流强度较高[见图3(a)]。

(2) 在t=1 000 s时，多孔材料填充区域相比相邻区域的流体流速小[见图3(b)]；相对于无多孔材料填充时，上层和下层的自然对流结构有明显不同，上层的自然对流结构由1个扁长形滚动圈变为几个眼泪形状的小滚动圈，且由于多孔材料增加了流体流动阻力，在多孔材料填充区域左侧流体流线竖直向下，穿过多孔材料后，便向气液交界面移动,而下层的自然对流结构由2个胖圆形和2个扁长形的滚动圈变为3个胖圆形的滚动圈，且中间的滚动圈比左右两侧的滚动圈大，滚动圈方向从左到右依次为逆时针、顺时针、逆时针，同时发现在壁面热边界条件和多孔材料作用下，罐体两边的流动强度较强，中间的流动强度较弱，中间滚动圈的滞留区面积比两侧滚动圈大[见图3(c)]。

储罐内沿中心区填充多孔材料厚度为1 m时，二维截面上流体的最大流速为0.190 m/s。

2.3 储罐中心区多孔材料填充厚度为1.5 m时储罐内LNG的流动情况

图4为储罐中心区多孔材料填充厚度为1.5 m时，t=100 s和t=1 000 s时刻储罐内LNG的流动情况模拟结果。

图4 储罐中心区多孔材料填充厚度为1.5 m时t=100 s和 t=1 000 s时刻的储罐内流体的速度云图和流线图Fig.4 Velocity cloud and stream traces of fluid in the storage tank with the filling thickness of porous materials being 1.5 m at t=100 s and t=1 000 s

由图4可以看出：

(1) 在t=100 s时，仍是边界层处流体和分界面处流体的流速较储罐其他位置的流速大，分界面处出现自然对流现象，但与同时期无多孔材料填充时相比，分界面上层流体的流速较高，在多孔材料填充区域两侧各出现了3个“火焰”形的高流速滚动圈，其速度分布同“火焰”的温度分布一样，中心速度最大，向四周减弱[见图4(a)]。

(2) 在t=1 000 s时，多孔材料填充区域的流体流速最小，且下层的自然对流结构明显沿储罐中心线呈对称分布[见图4(b)]；上层和下层的自然对流结构与无多孔材料填充时有明显的不同，上层的自然对流结构与多孔材料填充厚度为1 m时一样，仍为几个大小不一、眼泪形状的小滚动圈，而且在多孔材料填充区域内部有一个滚动圈，而下层的自然对流结构沿储罐中心线呈对称分布，由2个近似相等的胖圆形和2个扁长形的滚动圈组成，且中间的滚动圈比左右两侧的滚动圈大，中心线两侧的滚动圈方向依次为顺时针、逆时针、顺时针、逆时针，顺时针与逆时针交替出现，同时发现与多孔材料填充厚度为1 m时相似，中间滚动圈的滞留区面积比两侧滚动圈大，两侧滚动圈的滞留区面积几乎为0[见图4(c)]。

储罐内沿中心区填充多孔材料厚度为1.5 m时，二维截面上流体的最大流速为0.371 m/s。

2.4 储罐中心区多孔材料填充厚度为2 m时储罐内LNG的流动情况

图5为储罐中心区多孔材料填充厚度为2 m时，t=100 s和t=1 000 s时刻储罐内LNG的流动情况模拟结果。

图5 储罐中心区多孔材料填充厚度为2 m时t=100 s和t=1 000 s时刻的储罐内流体的速度云图和流线图Fig.5 Velocity cloud and stream traces of fluid in the storage tank with the filling thickness of porous materials being 2 m at t=100 s and t=1 000 s

由图5可以看出：

(1) 在t=100 s时，储罐内流体流动发展情况与同时期多孔材料填充厚度为1 m时较相似，分界面处的流体流速较无多孔材料填充时小，除边界层外，多孔材料填充区域两侧自然对流强度较高[见图5(a)]。

(2) 在t=1 000 s时，相对于无多孔材料填充时，低速区域范围变大，高速区域范围变小，且速度场与填充厚度为1.5 m时的速度场较相似[见图5(b)]；上层的自然对流结构为2个小的滚动圈，下层的自然对流结构主要为2个大的滚动圈，滚动圈方向依次为顺时针、逆时针，此时的流场结构不再沿储罐中心线呈对称分布，而是以中心线为分界线，左右两侧流场具有相似性，说明多孔材料填充厚度为2 m时，极大地增加了流动阻力，好像一堵墙，起到了分割流场的作用[见图5(c)]。

储罐内沿中心区填充多孔材料厚度为2 m时，二维截面上流体的最大流速为0.189 m/s。

2.5 储罐中心区多孔材料填充厚度对储罐内LNG流动情况的影响分析

图6为多孔材料不同填充厚度(H=0、1 m、1.5 m、2 m)时储罐内二维截面上流体平均流速的变化曲线。

图6 多孔材料不同填充厚度时储罐内二维截面上流体平均流速的变化曲线Fig.6 Graph of average flow velocity of the fluid at the two-dimensional cross-section with different thicknesses of porous materials

由图6可见，600 s是一个分界点，在600 s之前，填充多孔材料厚度为1.5 m的流体平均流速最高，这也与速度云图得到的结果相一致，在600 s之后，流体的自然对流发展到整个储罐，有多孔材料填充的储罐内流体的平均流速比无多孔材料填充的要低。

在气液交界面有蒸发气体进入到气相空间时，气相空间压力会增大，当其超过安全阀设定压力时，安全阀开启，释放蒸发气体到大气空间。当储罐内液体发生翻滚时，会瞬间产生大量蒸发气体，不仅会带来经济损失还会发生安全事故。图7为多孔材料不同填充厚度时储罐内压力出口的质量流量变化曲线。

图7 多孔材料不同填充厚度时储罐内压力出口的质量流量变化曲线Fig.7 Graph of mass flow rate at the pressure outlet of the storage tank with different thicknesses of porous materials

由图7可见，不管是哪种情况，储罐内压力出口的质量流量都为负值，但填充多孔材料的储罐均能不同程度地减少压力出口的质量流量，这样不仅可减少经济损失，还可减弱事故发生时的剧烈程度，起到一定的抑制作用。

图8为多孔材料不同填充厚度时储罐内上下分层之间的压差变化曲线。

由图8可见，在储罐中心区内填充多孔材料可以减小储罐内上下分层之间的压差，即减小储罐内翻滚事故发生时释放的能量，从而在一定程度上抑制了翻滚事故的发生。

图8 多孔材料不同填充厚度时储罐内上下分层之间的压差变化曲线Fig.8 Graph of pressure differentials of layers in the storage tank with different thicknesses of porous materials

3 结 论

通过对储罐中心区无多孔材料填充和有多孔材料填充时储罐内液化天然气的流动情况进行数值模拟与分析，得出以下结论：

(1) 在储罐中心区内填充多孔材料，多孔材料填充部分的流体流速明显比相邻区域小，同时在储罐壁面热边界条件和多孔材料作用下，罐体两边的流动强度较强，中间的流动强度较弱，中间滚动圈的滞留区面积比两侧滚动圈大，且滚动圈的方向是顺时针与逆时针交替出现。

(2) 由于储罐内液体翻滚的发生需要孕育时间和受计算机的限制，在此前对无多孔材料填充时储罐内液体翻滚进行模拟时需要很长的时间，因此本文在对储罐中心区填充多孔材料后储罐内液体的流动情况进行模拟时并未对整个过程进行模拟，只进行了前1 000 s时的模拟计算及分析，发现储罐中心区填充多孔材料后，可减小储罐内流体的平均流速，减少压力出口的质量流量，并减小储罐内上下分层之间的压差。

[1] 林瑞泰.多孔介质传热传质引论[M].北京:科学出版社,1995:39-47.

[2] Nield D A,Bejan A.ConvectioninPorousMedia[M].New York:Springer Velag,1999:23-65.

[3] 邢志祥,杜贞,张成燕,等.密闭储罐内填充非金属多孔材料后预混可燃气体火焰传播的数值模拟[J].安全与环境学报,2014,14(6):91-95.

[4] 张健中,许光,周金广,等.网状铝合金阻隔防爆材料功效及应用探讨[J].中国安全科学学报,2014,24(3):42-46.

[5] 田宏,王旭,高永庭.用于石油液化气体储罐填充的多孔金属材料的防火防爆机理及应用[J].消防技术与产品信息,2000(1):29-30.

[6] 田宏,王旭,高永庭.多孔填充材料的防火防爆机理及应用[J].工业安全与防尘,2000(4):43-46.

[7] 孙福涛,蒲亮,齐迪.大型LNG储罐分层破坏及翻滚过程的数值研究[J].低温工程,2017(1):47-53.

[8] 付全荣.泡沫金属填充套管换热器内流体流动和传热研究[D].太原：太原理工大学,2010.

[9] 于明跃.多孔固体构架与气流对流换热特性数值研究[D].哈尔滨：哈尔滨工业大学,2010.

[10] 李泊然.晃荡条件下LNG液货舱分层与翻滚现象的数值模拟[D].哈尔滨：哈尔滨工业大学,2016.

[11] 秦允毫.热学[M].2版.北京:高等教育出版社,2004:154-156.