楊 帆,陳保東,姜文全,賈永昌

(1.中國石油大學(xué)儲(chǔ)運(yùn)與建筑工程學(xué)院,山東青島266580；2.遼寧石油化工大學(xué)石油天然氣工程學(xué)院,遼寧撫順113001；3.中國石油青海油田采油一廠,青海海西州816400)

為了提高氣體的運(yùn)輸、儲(chǔ)存效率及降低運(yùn)行成本，天然氣“液態(tài)化”已成為天然氣行業(yè)發(fā)展的重要方向之一[1]。液化天然氣(LNG)的儲(chǔ)存和運(yùn)輸通常通過真空絕熱儲(chǔ)罐實(shí)現(xiàn)，一旦發(fā)生事故造成真空絕熱低溫容器外殼體破裂，絕熱夾層漏熱量的迅速增加使存儲(chǔ)的低溫LNG快速蒸發(fā)，進(jìn)而使低溫容器的壓力急劇增加，會(huì)使罐材料的微觀缺陷演變成裂紋并擴(kuò)展，給低溫容器的安全使用帶來極大威脅。有學(xué)者對(duì)絕熱夾層空間真空喪失后低溫容器的傳熱與升壓規(guī)律的問題進(jìn)行了試驗(yàn)與理論研究[2-3]，而對(duì)絕熱夾層空間完全真空喪失后的低溫容器罐體裂紋萌生及擴(kuò)展鮮有研究。為此，筆者應(yīng)用分子動(dòng)力學(xué)方法對(duì)絕熱夾層真空喪失過程中的LNG儲(chǔ)罐內(nèi)殼體裂紋萌生到擴(kuò)展生成過程進(jìn)行模擬研究，獲取在這種條件下的LNG儲(chǔ)罐裂紋擴(kuò)展的微觀機(jī)制，研究LNG低溫容器真空喪失過程中的壓力和應(yīng)力變化對(duì)罐體微裂紋擴(kuò)展的影響。

1 基本理論

采用分子動(dòng)力學(xué)法對(duì)帶有表面缺陷的LNG低溫儲(chǔ)罐常用材料0Cr18Ni9進(jìn)行研究，該方法通過允許原子和分子在一段時(shí)間內(nèi)交互作用，是一種根據(jù)已知的物理近似求解系統(tǒng)中所有粒子的運(yùn)動(dòng)方程，獲知原子運(yùn)動(dòng)過程的圖像的計(jì)算機(jī)模擬方法[4-7]。原子的運(yùn)動(dòng)軌跡和速度可由勢(shì)能和坐標(biāo)來確定，在給定初始條件下，選擇合理的原子間勢(shì)函數(shù)尤為重要。

1.1 原子間勢(shì)函數(shù)

采用鑲嵌原子勢(shì)(EAM)[8-10]描述鐵鎳合金(0Cr18Ni9)原子之間的相互作用，其理論表達(dá)式為

式中，E為總原子勢(shì)能；V為中心對(duì)勢(shì)；ρi為系統(tǒng)中所有其他原子在i原子處產(chǎn)生的局域背景電子密度；H(x)為Heaviside函數(shù)；參數(shù)ak、rk、Ak、Rk由擬合晶體的彈性常數(shù)、點(diǎn)陣常數(shù)、空位形成能、結(jié)合能以及壓強(qiáng)與體積關(guān)系確定。

1.2 計(jì)算方法

積分算法采用Swope等[11]提出的速度Verlet算法，其算法從當(dāng)前時(shí)刻t的位置矢量r(t)和速度矢量Vi(t)出發(fā)可得

在計(jì)算了ri(t+Δt)后，可由當(dāng)前時(shí)刻的力Fi(t)以及下一個(gè)時(shí)刻的力Fi(t+Δt)求得相應(yīng)的平均加速度，然后確定下一時(shí)刻速度為

可同時(shí)得出位置、速度和加速度等參數(shù)，應(yīng)用LAMMPS軟件編寫開放性程序計(jì)算得出模型的壓力、應(yīng)力和應(yīng)變。

2 分子動(dòng)力學(xué)模擬

2.1 模型建立

對(duì)LNG真空絕熱儲(chǔ)罐進(jìn)行模擬分析，儲(chǔ)罐內(nèi)殼帶有I型表面微裂紋，罐體結(jié)構(gòu)及微裂紋方向如圖1所示。儲(chǔ)罐內(nèi)殼體材料為0Cr18Ni9不銹鋼，彈性模量為206 GPa，其他力學(xué)性能如表1所示。不銹鋼基體相α-Fe是體心立方結(jié)構(gòu)(bcc)，晶格常數(shù)為2.855 ?。對(duì)微裂紋尖端進(jìn)行原子尺度研究，建立裂尖尺寸如圖2所示。模型寬度為2b，裂紋長度為2a(b=114.2 ?，a=14.275 ?)。為了減小分子動(dòng)力學(xué)模型的計(jì)算規(guī)模，利用模型沿x軸的對(duì)稱性，僅建立分子動(dòng)力學(xué)的1/2模型；裂紋面垂直于z軸，與xoy面平行；模型共包含7334個(gè)原子(圖2(b))。

圖1 LNG儲(chǔ)罐結(jié)構(gòu)Fig.1 LNG tank structure

表1 不銹鋼0Cr18Ni9力學(xué)性能Table 1 Mechanical properties of stainless steel 0Cr18Ni9

圖2 分子動(dòng)力學(xué)模型Fig.2 Molecular dynamics model

2.2 計(jì)算過程

模擬采用Fe-Ni合金的EAM勢(shì)函數(shù)，邊界條件為pss(即載荷方向?yàn)橹芷谛赃吔?，其他方向?yàn)樽杂蛇吔?，通過強(qiáng)行約束對(duì)稱上邊界原子的位移實(shí)現(xiàn)對(duì)稱約束條件；系統(tǒng)采用NVE系綜，在x方向模型兩端原子上施加拉應(yīng)力載荷，時(shí)間步長為0.002 ps，運(yùn)行10 000步，每500步輸出一次結(jié)果。同時(shí)，通過陳篪[13]推導(dǎo)出的臨界應(yīng)力強(qiáng)度因子公式(其中，k≈0.5，w為壁厚)求得該材料臨界應(yīng)力強(qiáng)度因子。計(jì)算結(jié)果為KC=0.375 MPa·m1/2。

3 結(jié)果及其討論

在x向逐漸增加的拉應(yīng)力作用下，模擬了LNG低溫儲(chǔ)罐夾層真空喪失后殼體的升壓過程，分別求得壓力、應(yīng)力及應(yīng)變，如表2所示。

表2 計(jì)算結(jié)果Table 2 Calculation results

3.1 壓力變化

在模擬過程中，通過LAMMPS自編程序求得每500步對(duì)應(yīng)的模型所受壓力，其變化規(guī)律如圖3所示。在載荷加載初期(0～2 999步)時(shí)，儲(chǔ)罐內(nèi)壓(即模型內(nèi)表面壓力)迅速上升；當(dāng)達(dá)到3000步時(shí)，壓力為最大值5403.109 MPa；在3001～10000步之間，壓力開始回落，最后在4200 MPa上下波動(dòng)并趨于穩(wěn)定。此壓力變化規(guī)律與謝高峰、汪榮順、魏蔚等[12]的試驗(yàn)研究結(jié)果吻合較好。其原因分析可知：LNG低溫儲(chǔ)罐夾層真空喪失后，自然對(duì)流傳熱使得罐內(nèi)溫度升高，部分LNG氣化造成罐內(nèi)壓力急劇升高；待溫度穩(wěn)定后，壓力也趨于平穩(wěn)。在升壓過程中，殼體應(yīng)力也隨之升高；壓力達(dá)到極值時(shí)，應(yīng)力也是最大值，裂紋開始擴(kuò)展。

圖3 模擬過程中內(nèi)壓變化Fig.3 Change of inter pressure during simulation process

3.2 裂紋尖端原子變化

圖4 裂紋附近原子圖Fig.4 Atomic diagrams near crack

圖4給出了4個(gè)不同時(shí)間步的裂紋尖端原子變化。(a)為模型受載荷作用初期，裂紋尖端兩側(cè)出現(xiàn)位錯(cuò)；此時(shí)應(yīng)為最大值12.638 44 GPa。(b)所示當(dāng)載荷加載到4500步時(shí)，裂紋進(jìn)一步向兩側(cè)張開，第一個(gè)位錯(cuò)向前運(yùn)動(dòng)，同時(shí)有第二個(gè)位錯(cuò)產(chǎn)生與第一個(gè)位錯(cuò)疊加，模型應(yīng)力降低為11.098 65 GPa。(c)所示當(dāng)載荷加載到8500步時(shí)，裂紋尖端右側(cè)有多個(gè)位錯(cuò)發(fā)射，同時(shí)先前產(chǎn)生的位錯(cuò)都不同程度地向前運(yùn)動(dòng)，應(yīng)力值降低為9.233 6 GPa，并逐漸趨于穩(wěn)定。(d)所示當(dāng)載荷加載到10 000步時(shí)，由于位錯(cuò)的累積和運(yùn)動(dòng)，裂紋尖端右側(cè)原子鍵斷裂，裂紋沿位錯(cuò)運(yùn)動(dòng)方向擴(kuò)展?？梢姀奈⒂^角度分析裂紋擴(kuò)展主要是由于裂紋尖端原子有位錯(cuò)產(chǎn)生并累積和運(yùn)動(dòng)的結(jié)果。

3.3 應(yīng)力-應(yīng)變關(guān)系

應(yīng)用LAMMPS程序求得每間隔500步的應(yīng)力、應(yīng)變值，應(yīng)力-應(yīng)變關(guān)系如圖5所示。在0～2 999步之間，應(yīng)力-應(yīng)變近似呈線性上升，裂紋尖端無明顯變化；達(dá)3 000步時(shí)，應(yīng)力達(dá)到最大值12.638 44 GPa、應(yīng)變?yōu)?.061 2(表1)，裂紋尖端原子發(fā)生位錯(cuò)；在3001～10000步之間，應(yīng)變繼續(xù)增加，應(yīng)力呈下降趨勢(shì)，此過程中裂紋尖端原子位錯(cuò)不斷產(chǎn)生并向前運(yùn)動(dòng)，第10000步時(shí)應(yīng)力值為9.024731GPa，可算得應(yīng)力強(qiáng)度因子為，微裂紋沿位錯(cuò)運(yùn)動(dòng)方向發(fā)生擴(kuò)展。

圖5 應(yīng)力-應(yīng)變關(guān)系曲線Fig.5 Stress-strain relation curve

4 結(jié) 論

(1)LNG低溫儲(chǔ)罐夾層真空喪失后，自然對(duì)流傳熱使得罐內(nèi)溫度升高，部分LNG氣化造成罐內(nèi)壓力急劇升高，殼體應(yīng)力也隨之升高；待溫度穩(wěn)定后，壓力和應(yīng)力都趨于平穩(wěn)。

(2)載荷初期(0～2999步)，應(yīng)力增加，裂紋尖端原子無明顯變化；3 000步時(shí)，裂紋尖端原子出現(xiàn)位錯(cuò)；4500步時(shí)，第二個(gè)位錯(cuò)出現(xiàn)并累積；8 500步時(shí)，多個(gè)位錯(cuò)產(chǎn)生，并向前運(yùn)動(dòng)；10000步時(shí)，裂紋尖端原子鍵斷裂，裂紋宏觀擴(kuò)展。裂紋擴(kuò)展是在應(yīng)力、應(yīng)變作用下，裂紋尖端原子位錯(cuò)產(chǎn)生、累積和運(yùn)動(dòng)的結(jié)果。

(3)開始階段，應(yīng)力-應(yīng)變近似呈線性上升；應(yīng)力達(dá)到極值后，裂紋尖端原子產(chǎn)生位錯(cuò)，應(yīng)變繼續(xù)增加，應(yīng)力降低；同時(shí)，應(yīng)力強(qiáng)度因子KI＞KC，微裂紋已擴(kuò)展。

[1] 程進(jìn)杰，朱建炳，李正清.低溫容器高真空多層絕熱性能分析[J].低溫與超導(dǎo)，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].真空科學(xué)與技術(shù)學(xué)報(bào)，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液艙內(nèi)儲(chǔ)液晃動(dòng)特性的數(shù)值模擬[J].中國石油大學(xué)學(xué)報(bào)：自然科學(xué)版，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] 張晶晶，許輝，祝新發(fā)，等.反應(yīng)濺射HfC/Si3N4納米多層膜的微結(jié)構(gòu)與力學(xué)性能[J].真空科學(xué)與技術(shù)學(xué)報(bào)，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ī)理的原子級(jí)模擬研究[J].物理學(xué)報(bào)，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] 徐榮鵬，曾祥國，陳華燕，等.裂紋尖端位錯(cuò)發(fā)射與運(yùn)動(dòng)的分子動(dòng)力學(xué)模擬[J].四川大學(xué)學(xué)報(bào)：自然科學(xué)版，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，盧磊.薄膜材料-應(yīng)力、缺陷的形成和表面演化[M].北京：科學(xué)出版社，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].北京：科學(xué)出版社，2008：40-131.

[12] 謝高峰汪榮順魏蔚.高真空多層絕熱低溫容器完全真空喪失后升壓規(guī)律的試驗(yàn)研究[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] 陳篪.論裂紋擴(kuò)展的判據(jù)[J].金屬學(xué)報(bào)，1977，13(1/2)：57-72.CHEN Chi.On the criterion for crack extension[J].Acta Metallurgica Sinica，1977，13(1/2)：57-72.