LNG低温储罐夹层真空丧失过程内壳裂纹扩展研究

杨 帆,陈保东,姜文全,贾永昌

(1.中国石油大学储运与建筑工程学院,山东青岛266580；2.辽宁石油化工大学石油天然气工程学院,辽宁抚顺113001；3.中国石油青海油田采油一厂,青海海西州816400)

为了提高气体的运输、储存效率及降低运行成本，天然气“液态化”已成为天然气行业发展的重要方向之一[1]。液化天然气(LNG)的储存和运输通常通过真空绝热储罐实现，一旦发生事故造成真空绝热低温容器外壳体破裂，绝热夹层漏热量的迅速增加使存储的低温LNG快速蒸发，进而使低温容器的压力急剧增加，会使罐材料的微观缺陷演变成裂纹并扩展，给低温容器的安全使用带来极大威胁。有学者对绝热夹层空间真空丧失后低温容器的传热与升压规律的问题进行了试验与理论研究[2-3]，而对绝热夹层空间完全真空丧失后的低温容器罐体裂纹萌生及扩展鲜有研究。为此，笔者应用分子动力学方法对绝热夹层真空丧失过程中的LNG储罐内壳体裂纹萌生到扩展生成过程进行模拟研究，获取在这种条件下的LNG储罐裂纹扩展的微观机制，研究LNG低温容器真空丧失过程中的压力和应力变化对罐体微裂纹扩展的影响。

1 基本理论

采用分子动力学法对带有表面缺陷的LNG低温储罐常用材料0Cr18Ni9进行研究，该方法通过允许原子和分子在一段时间内交互作用，是一种根据已知的物理近似求解系统中所有粒子的运动方程，获知原子运动过程的图像的计算机模拟方法[4-7]。原子的运动轨迹和速度可由势能和坐标来确定，在给定初始条件下，选择合理的原子间势函数尤为重要。

1.1 原子间势函数

采用镶嵌原子势(EAM)[8-10]描述铁镍合金(0Cr18Ni9)原子之间的相互作用，其理论表达式为

式中，E为总原子势能；V为中心对势；ρi为系统中所有其他原子在i原子处产生的局域背景电子密度；H(x)为Heaviside函数；参数ak、rk、Ak、Rk由拟合晶体的弹性常数、点阵常数、空位形成能、结合能以及压强与体积关系确定。

1.2 计算方法

积分算法采用Swope等[11]提出的速度Verlet算法，其算法从当前时刻t的位置矢量r(t)和速度矢量Vi(t)出发可得

在计算了ri(t+Δt)后，可由当前时刻的力Fi(t)以及下一个时刻的力Fi(t+Δt)求得相应的平均加速度，然后确定下一时刻速度为

可同时得出位置、速度和加速度等参数，应用LAMMPS软件编写开放性程序计算得出模型的压力、应力和应变。

2 分子动力学模拟

2.1 模型建立

对LNG真空绝热储罐进行模拟分析，储罐内壳带有I型表面微裂纹，罐体结构及微裂纹方向如图1所示。储罐内壳体材料为0Cr18Ni9不锈钢，弹性模量为206 GPa，其他力学性能如表1所示。不锈钢基体相α-Fe是体心立方结构(bcc)，晶格常数为2.855 Å。对微裂纹尖端进行原子尺度研究，建立裂尖尺寸如图2所示。模型宽度为2b，裂纹长度为2a(b=114.2 Å，a=14.275 Å)。为了减小分子动力学模型的计算规模，利用模型沿x轴的对称性，仅建立分子动力学的1/2模型；裂纹面垂直于z轴，与xoy面平行；模型共包含7334个原子(图2(b))。

图1 LNG储罐结构Fig.1 LNG tank structure

表1 不锈钢0Cr18Ni9力学性能Table 1 Mechanical properties of stainless steel 0Cr18Ni9

图2 分子动力学模型Fig.2 Molecular dynamics model

2.2 计算过程

模拟采用Fe-Ni合金的EAM势函数，边界条件为pss(即载荷方向为周期性边界，其他方向为自由边界)，通过强行约束对称上边界原子的位移实现对称约束条件；系统采用NVE系综，在x方向模型两端原子上施加拉应力载荷，时间步长为0.002 ps，运行10 000步，每500步输出一次结果。同时，通过陈篪[13]推导出的临界应力强度因子公式(其中，k≈0.5，w为壁厚)求得该材料临界应力强度因子。计算结果为KC=0.375 MPa·m1/2。

3 结果及其讨论

在x向逐渐增加的拉应力作用下，模拟了LNG低温储罐夹层真空丧失后壳体的升压过程，分别求得压力、应力及应变，如表2所示。

表2 计算结果Table 2 Calculation results

3.1 压力变化

在模拟过程中，通过LAMMPS自编程序求得每500步对应的模型所受压力，其变化规律如图3所示。在载荷加载初期(0～2 999步)时，储罐内压(即模型内表面压力)迅速上升；当达到3000步时，压力为最大值5403.109 MPa；在3001～10000步之间，压力开始回落，最后在4200 MPa上下波动并趋于稳定。此压力变化规律与谢高峰、汪荣顺、魏蔚等[12]的试验研究结果吻合较好。其原因分析可知：LNG低温储罐夹层真空丧失后，自然对流传热使得罐内温度升高，部分LNG气化造成罐内压力急剧升高；待温度稳定后，压力也趋于平稳。在升压过程中，壳体应力也随之升高；压力达到极值时，应力也是最大值，裂纹开始扩展。

图3 模拟过程中内压变化Fig.3 Change of inter pressure during simulation process

3.2 裂纹尖端原子变化

图4 裂纹附近原子图Fig.4 Atomic diagrams near crack

图4给出了4个不同时间步的裂纹尖端原子变化。(a)为模型受载荷作用初期，裂纹尖端两侧出现位错；此时应为最大值12.638 44 GPa。(b)所示当载荷加载到4500步时，裂纹进一步向两侧张开，第一个位错向前运动，同时有第二个位错产生与第一个位错叠加，模型应力降低为11.098 65 GPa。(c)所示当载荷加载到8500步时，裂纹尖端右侧有多个位错发射，同时先前产生的位错都不同程度地向前运动，应力值降低为9.233 6 GPa，并逐渐趋于稳定。(d)所示当载荷加载到10 000步时，由于位错的累积和运动，裂纹尖端右侧原子键断裂，裂纹沿位错运动方向扩展。可见从微观角度分析裂纹扩展主要是由于裂纹尖端原子有位错产生并累积和运动的结果。

3.3 应力-应变关系

应用LAMMPS程序求得每间隔500步的应力、应变值，应力-应变关系如图5所示。在0～2 999步之间，应力-应变近似呈线性上升，裂纹尖端无明显变化；达3 000步时，应力达到最大值12.638 44 GPa、应变为0.061 2(表1)，裂纹尖端原子发生位错；在3001～10000步之间，应变继续增加，应力呈下降趋势，此过程中裂纹尖端原子位错不断产生并向前运动，第10000步时应力值为9.024731GPa，可算得应力强度因子为，微裂纹沿位错运动方向发生扩展。

图5 应力-应变关系曲线Fig.5 Stress-strain relation curve

4 结 论

(1)LNG低温储罐夹层真空丧失后，自然对流传热使得罐内温度升高，部分LNG气化造成罐内压力急剧升高，壳体应力也随之升高；待温度稳定后，压力和应力都趋于平稳。

(2)载荷初期(0～2999步)，应力增加，裂纹尖端原子无明显变化；3 000步时，裂纹尖端原子出现位错；4500步时，第二个位错出现并累积；8 500步时，多个位错产生，并向前运动；10000步时，裂纹尖端原子键断裂，裂纹宏观扩展。裂纹扩展是在应力、应变作用下，裂纹尖端原子位错产生、累积和运动的结果。

(3)开始阶段，应力-应变近似呈线性上升；应力达到极值后，裂纹尖端原子产生位错，应变继续增加，应力降低；同时，应力强度因子KI＞KC，微裂纹已扩展。

[1] 程进杰，朱建炳，李正清.低温容器高真空多层绝热性能分析[J].低温与超导，2013，41(2)：11-14.CHENG Jin-jie，ZHU Jian-bing，LI Zheng-qing.Analysis on the performance of high vacuum multilayer insulation for cryogenic storage vessel[J].Cryogenics and Superconductivity，2013，41(2)：11-14.

[2] 陈树军，谭粤，杨树斌，等.低温绝热气瓶漏放气性能的研究[J].真空科学与技术学报，2012，32(5)：447-451.CHEN Shu-jun，TAN Yue，YANG Shu-bin，et al.Study of leakage and out-gassing of cryogenic insulated cylinder[J].Journal of Vacuum Science and Technology，2012，32(5)：447-451.

[3] 陈海阳，李玉星，孙法峰，等.LNG FPSO液舱内储液晃动特性的数值模拟[J].中国石油大学学报：自然科学版，2011，35(4)：134-139.CHEN Hai-yang，LI Yu-xing，SUN Fa-feng，et al.Numerical simulation of liquid sloshing characteristics in LNG FPSO containers[J].Journal of China University of Petroleum(Edition of Natural Science)，2011，35(4)：134-139.

[4] 张晶晶，许辉，祝新发，等.反应溅射HfC/Si3N4纳米多层膜的微结构与力学性能[J].真空科学与技术学报，2011，31(5)：574-578.ZHANG Jing-jing，XU Hui，ZHU Xin-fa，et al.Growth and characterization of super-hard HfC/Si3N4 multilayers[J].Journal of Vacuum Science and Technology，2011，31(5)：574-578.

[5] 张杨，张建华，文玉华，等.含圆孔纳米薄膜在拉伸加载下变形机理的原子级模拟研究[J].物理学报，2008，57(11)：7094-7099.ZHANG Yang，ZHANG Jian-hua，WEN YU-Hua，et al.The deformation mechanism of nanofilm with void under tensile loading：an atomistic simulation study[J].Acta Physica Sinica，2008，57(11)：7094-7099.

[6] BENNAN H Honmid，TAKADOUM J.Finite element model of elastic stresses in thin coatings submitted to applied forces[J].Surface and Coatings Technology，1999，111：80-85.

[7] 徐荣鹏，曾祥国，陈华燕，等.裂纹尖端位错发射与运动的分子动力学模拟[J].四川大学学报：自然科学版，2011，48(3)：612-616.XU Rong-peng，ZENG Xiang-guo，CHEN Hua-yan，et al.The molecular dynamics simulation of the dislocation emission and movement around a crack tip[J].Journal of Sichuan University(Natural Science Edition)，2011，48(3)：612-616.

[8] FREUND L B，SURESH S，卢磊.薄膜材料-应力、缺陷的形成和表面演化[M].北京：科学出版社，2006.

[9] SHEN Lu-ming.Combined grain size，strain rate and loading condition effects on mechanical behavior of Nan crystalline Cu under high strain rates[J].Chinese Journal of Theoretical and Applied Mechanics，2012，28(4)：1125-1132.

[10] SWOPE W C，ANDERSEN H C，BERENS P H，et al.A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules：application to small water clusters[J].Chem Phys，1982，76：637.

[11] 范镜泓.材料变形与破坏的多尺度分析[M].北京：科学出版社，2008：40-131.

[12] 谢高峰汪荣顺魏蔚.高真空多层绝热低温容器完全真空丧失后升压规律的试验研究[J].压力容器，2011，28(4)：12-17.XIE Gao-feng，WANG Rong-shun，WEI Wei.Experimental investigation of the pressure rise progress of HVMLI cryogenic tanks after catastrophic loss of insulating vacuum[J].Pressure Vessel Technology，2011，28(4)：12-17.

[13] 陈篪.论裂纹扩展的判据[J].金属学报，1977，13(1/2)：57-72.CHEN Chi.On the criterion for crack extension[J].Acta Metallurgica Sinica，1977，13(1/2)：57-72.