張 娜 李兆慈 郭志超

中國石油大學(xué)(北京)油氣管道輸送安全國家工程實(shí)驗(yàn)室/城市油氣輸配技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室, 北京 102200

0 前言

在LNG產(chǎn)業(yè)鏈中,儲(chǔ)存是重要的一環(huán)。天然氣低溫液化后經(jīng)遠(yuǎn)洋運(yùn)輸?shù)竭_(dá)世界各地LNG接收站,接收站既是海上LNG運(yùn)輸?shù)慕K點(diǎn)也是陸上天然氣供應(yīng)的起點(diǎn)[1]。儲(chǔ)罐是接收站內(nèi)最重要的儲(chǔ)存設(shè)備,它的安全至關(guān)重要。接收站的儲(chǔ)罐體積龐大,結(jié)構(gòu)復(fù)雜,內(nèi)罐工作在-162 ℃左右的低溫環(huán)境,因此儲(chǔ)罐在安全運(yùn)行方面存在諸多需要注意的問題[2]。隨著LNG儲(chǔ)存設(shè)施大型化趨勢(shì)和高安全性的要求,全容罐因?yàn)榫哂懈叩陌踩缘玫搅藦V泛應(yīng)用[3]。LNG儲(chǔ)罐作為接收站最核心的建筑物,是重大危險(xiǎn)源,一旦發(fā)生事故將會(huì)造成巨大損失。當(dāng)儲(chǔ)罐因?yàn)槿藶椴僮鞑划?dāng)、設(shè)備損壞、管道破裂等問題發(fā)生LNG泄漏時(shí),由于LNG是-162 ℃的深冷液體,泄漏出來的低溫液體會(huì)對(duì)設(shè)備造成低溫沖擊,危害儲(chǔ)罐結(jié)構(gòu)的安全[4]。如LNG泄漏于地面時(shí),會(huì)迅速從地面、周圍環(huán)境吸收熱量,蒸發(fā)為氣體。LNG全容罐分內(nèi)罐和外罐,在設(shè)計(jì)和建造上能容納所儲(chǔ)存的低溫液體以及液體蒸發(fā)后產(chǎn)生的低溫氣體。由于LNG接收站一般位于沿海地帶,可能發(fā)生地震等自然災(zāi)害,增加了LNG儲(chǔ)罐內(nèi)罐泄漏的風(fēng)險(xiǎn)。所以采用合理的方法對(duì)LNG儲(chǔ)罐內(nèi)罐泄漏進(jìn)行模擬研究,對(duì)于外罐壁厚的設(shè)置、制定發(fā)生泄漏時(shí)采取的應(yīng)急措施都具有重要意義。

1 LNG儲(chǔ)罐泄漏研究方法

因低溫導(dǎo)致設(shè)備材料斷裂,從而引起的儲(chǔ)罐泄漏,一般可根據(jù)泄漏口面積大小和泄漏處高度等因素分成兩類[5]:一類是大面積泄漏,在較短時(shí)間內(nèi)儲(chǔ)液大量漏出,瞬間產(chǎn)生大量BOG氣體,可能在儲(chǔ)罐內(nèi)形成瞬間高壓發(fā)生爆炸[6];另一類是有限孔泄漏,儲(chǔ)存的液體沿小泄漏口緩慢流入絕熱層,破壞絕熱層結(jié)構(gòu),對(duì)儲(chǔ)罐的正常運(yùn)行造成影響。內(nèi)罐小孔泄漏是目前最易發(fā)生的泄漏工況,本文主要針對(duì)LNG內(nèi)罐小孔泄漏工況進(jìn)行計(jì)算模擬。

對(duì)于LNG儲(chǔ)罐泄漏問題,早期通過現(xiàn)場(chǎng)模擬實(shí)驗(yàn),在野外搭建與真實(shí)事故完全一致的場(chǎng)景,重現(xiàn)當(dāng)時(shí)的后果影響,得到最真實(shí)的數(shù)據(jù)資料。但這種現(xiàn)場(chǎng)試驗(yàn)隨機(jī)性很大,且需花費(fèi)大量人力物力,得到的結(jié)果僅對(duì)某一特定的事故有幫助,但對(duì)類似事故參考意義不大,因此需要更簡(jiǎn)單的具有普適性的方法[7]。通過物理模型構(gòu)建類似實(shí)際情況,可在實(shí)驗(yàn)室完成的小型泄漏模擬實(shí)驗(yàn)來研究實(shí)際LNG泄漏擴(kuò)散規(guī)律。研究LNG泄漏擴(kuò)散過程的物理模擬一般采用風(fēng)洞試驗(yàn)[8]。20世紀(jì)70、80年代,MeroneyR N[9]已經(jīng)模擬出LNG高空和地面釋放情況下的LNG擴(kuò)散。

數(shù)值模擬是應(yīng)用比較普遍的一種借助計(jì)算機(jī)完成的LNG儲(chǔ)罐泄漏擴(kuò)散模擬方法。這種方法無需設(shè)置實(shí)驗(yàn),可借助不同的數(shù)學(xué)模型完成,經(jīng)濟(jì)省時(shí),可以得到詳盡的結(jié)果。隨著計(jì)算機(jī)計(jì)算能力的不斷提高,加上有限容積法、有限元分析法等數(shù)值計(jì)算方法的發(fā)展,基于數(shù)值計(jì)算的計(jì)算流體力學(xué)(CFD)方法得到了蓬勃發(fā)展[10]。

LNG泄漏后的擴(kuò)散過程可通過速度、濃度和溫度進(jìn)行描述,其控制方程包括連續(xù)性方程、動(dòng)量方程、能量方程、擴(kuò)散方程、湍流動(dòng)能方程、湍流動(dòng)能耗散方程[11]等。

現(xiàn)場(chǎng)實(shí)驗(yàn)會(huì)耗費(fèi)較大的人力、物力和財(cái)力,同時(shí)外界的環(huán)境變化復(fù)雜,在某些情況下難以獲取理想的數(shù)據(jù)。風(fēng)洞實(shí)驗(yàn)在低風(fēng)速和低雷諾數(shù)條件下進(jìn)行比較困難[11]。幾種方法互有長短,實(shí)際研究中可綜合運(yùn)用,但鑒于數(shù)值模擬的耗時(shí)短、成本低以及高度可重復(fù)性,近幾年成為此類研究中最廣泛使用的方法[12]。

對(duì)比LNG泄漏擴(kuò)散的各種研究方法可知,CFD模擬具有較好的適用性,但要準(zhǔn)確模擬LNG泄漏擴(kuò)散過程,對(duì)初始條件、邊界條件設(shè)置要求較高[13]。同時(shí)CFD模擬還存在一些問題需要解決,如過分簡(jiǎn)化處理源項(xiàng)和地面?zhèn)鳠岱绞?對(duì)于關(guān)鍵影響因素的敏感性分析不足等。

2 內(nèi)罐泄漏溫度場(chǎng)計(jì)算模擬

隨著內(nèi)罐小孔泄漏的持續(xù),內(nèi)罐受影響區(qū)域逐漸向保溫層外側(cè)和底部發(fā)展,危及LNG外罐的安全。

在LNG儲(chǔ)罐內(nèi)罐發(fā)生泄漏時(shí),LNG會(huì)沿有限孔緩慢流出,再沿著儲(chǔ)罐內(nèi)壁外的彈性玻璃纖維氈流入環(huán)形空間的絕熱層中。絕熱層由彈性氈、膨脹珍珠巖和泡沫玻璃磚組合而成[14]。以 20 000 m3雙金屬全容式LNG儲(chǔ)罐為例,內(nèi)罐高度30 m,內(nèi)罐罐徑(直徑)30 m,外罐罐壁高度31 m,外罐罐徑(直徑)32.378 m，外罐拱頂高度4.66 m。內(nèi)罐由底板、頂板及10層壁板組成,外罐由底板、頂板、10層壁板及梯子欄桿組成。

由于內(nèi)罐整體高度為30 m,所以每層壁板高度為3 m。內(nèi)、外罐壁每層壁板厚度見表1。

表1 內(nèi)、外罐各層壁板厚度表

Tab.1 Thickness of each layer of the tank innerwall and outerwall

罐壁板層內(nèi)罐壁板厚度/mm外罐壁板厚度/mm第1層1612第2層1212第3層810第4層810第5層68第6層68第7層68第8層68第9層68第10層68

儲(chǔ)罐絕熱層結(jié)構(gòu)見圖1，內(nèi)壁絕熱材料的厚度及導(dǎo)熱系數(shù)等參數(shù)見表2。

圖1 罐壁絕熱層結(jié)構(gòu)示意圖Fig.1 The wall insulation of the tank

表2 儲(chǔ)罐罐壁絕熱層數(shù)據(jù)表

Tab.2 Wall insulation of the tank

材料厚度/mm導(dǎo)熱系數(shù)/(W·m-1 ·K-1)熱容/(J·kg-1 ·K-1)玻璃纖維氈3000.043792.00珠光砂7000.030753.74泡沫玻璃磚1500.048837.49

絕熱材料是由氣相與固相構(gòu)成的多孔性組織,固體材料內(nèi)部存在互相連通的氣相空間,這給LNG泄漏提供了一條滲流通道,因此LNG儲(chǔ)罐內(nèi)罐泄漏模擬可采用多孔介質(zhì)模型[15]。

在內(nèi)罐泄漏條件下,對(duì)儲(chǔ)罐的溫度邊界條件作以下假設(shè):外罐壁初始溫度取正常操作條件下的環(huán)境溫度20°C,外罐壁與空氣的對(duì)流換熱系數(shù)取25 W/(m2·K)；儲(chǔ)罐內(nèi)LNG溫度為-162°C,內(nèi)罐溫度與儲(chǔ)液溫度相同。不考慮LNG泄漏時(shí)的相變，外罐壁和罐底部絕熱層的下部均為絕熱壁面,外罐壁頂部也設(shè)為絕熱壁面。由此得到穩(wěn)態(tài)時(shí)罐壁絕熱層的溫度分布,見圖2。

圖2 穩(wěn)態(tài)溫度分布圖Fig.2 Temperature distribution in steady state

儲(chǔ)罐壁高30 m,故i的取值范圍在0～3 000,在距離罐底10 m(i=1 000)處設(shè)置泄漏工況,泄漏孔徑為10 mm,泄漏點(diǎn)流速為0.02 m/s。

考慮泄漏點(diǎn)處對(duì)流換熱,設(shè)置泄漏點(diǎn)(i=1 000)處溫度為LNG的溫度,與上下邊界進(jìn)行對(duì)流換熱,由傅立葉定律和牛頓冷卻公式可知其熱流量相等,即：

(1)

式中:h為對(duì)流換熱系數(shù),W/(m2·K);Tw為壁面溫度,℃;Tf為流體溫度,℃;T為絕熱層溫度，℃;λ為介質(zhì)導(dǎo)熱系數(shù),W/(m·K)；Y為內(nèi)壁高度方向。

其中,對(duì)流換熱系數(shù)h是一個(gè)與流體密度、比熱容、流速、黏度、導(dǎo)熱系數(shù)有關(guān)的復(fù)雜函數(shù),h=f(ρ,cp,λ,u,l,μ),可通過計(jì)算雷諾數(shù)Ref、努賽爾數(shù)Nuf和普朗特?cái)?shù)Prf得到[16],具體計(jì)算按式(2)。

(2)

式中:u為流體流速,m/s;l為當(dāng)量長度,m;μ為流體黏度,m2/s。對(duì)于非圓形截面的槽道,l取當(dāng)量直徑作為特征長度de。

(3)

式中:Ac為槽道的流動(dòng)截面積,m2;P為潤濕周長,即槽道壁與流體接觸面的長度,m。

再計(jì)算普朗特?cái)?shù):

(4)

式中:Cp為流體比熱容,J/(kg·K);a為流體導(dǎo)熱系數(shù),W/(m·K)。

最后再由努賽爾數(shù)求得流體的對(duì)流換熱系數(shù)h:

(5)

將內(nèi)罐泄漏后的罐壁傳熱簡(jiǎn)化成二維非穩(wěn)態(tài)對(duì)流換熱問題,先求得泄漏流體在罐壁內(nèi)的流場(chǎng)分布。對(duì)彈性多孔介質(zhì)單相不可壓縮流體不穩(wěn)定滲流問題,其運(yùn)動(dòng)方程為:

(6)

式中:ν為流體泄漏速度,m/s;K為介質(zhì)滲透率,μm2。

對(duì)于二維平面:

(7)

對(duì)于彈性孔隙介質(zhì)，其狀態(tài)方程(多孔介質(zhì)和液體都是可壓縮的):

φ=φa+Cf(p-pa)

(8)

式中:φ為孔隙介質(zhì)度；Cf為流體的壓縮系數(shù)；φa為標(biāo)準(zhǔn)大氣壓下孔隙介質(zhì)度；pa為標(biāo)準(zhǔn)大氣壓，Pa;p為流體壓力，Pa。

對(duì)于彈性液體:

ρ=ρa(bǔ)eCL(p-pa)=ρa(bǔ)[1+CL(p-pa)]

(9)

式中:ρ為液體密度,kg/m3;ρa(bǔ)為標(biāo)準(zhǔn)大氣壓下液體密度,kg/m3;CL為液體彈性壓縮系數(shù)。

單相流體的連續(xù)性方程為:

(10)

將運(yùn)動(dòng)方程和狀態(tài)方程代入連續(xù)性方程可得:

(11)

通過對(duì)設(shè)置的儲(chǔ)罐內(nèi)壁絕熱層進(jìn)行網(wǎng)格劃分,確定出邊界點(diǎn)、角點(diǎn)和內(nèi)點(diǎn),再對(duì)二維非穩(wěn)態(tài)滲流壓力方程進(jìn)行離散,最后采用TDMA算法加快迭代計(jì)算速度,由此計(jì)算模擬可得到LNG泄漏至內(nèi)罐壁的流場(chǎng)分布,其橫坐標(biāo)為內(nèi)壁厚度,縱坐標(biāo)為內(nèi)壁高度,罐底至泄漏點(diǎn)即i=1 000下方的壓力分布見圖3。

由圖3可知,在計(jì)算步長內(nèi),其壓力在泄漏點(diǎn)處最大,逐漸向下削弱。通過對(duì)x和y方向上的壓力求偏導(dǎo),可計(jì)算得到二維速度場(chǎng)的分布,圖4為罐壁泄漏點(diǎn)(i=1 000)下方二維速度分布矢量圖。

由模擬結(jié)果可以看出,由于罐壁絕熱層材料物理特性不同,其滲透擴(kuò)散速率也不同。將模擬計(jì)算得到的流場(chǎng)分布與溫度結(jié)合,將模擬計(jì)算得到的速度分布作為對(duì)流擴(kuò)散項(xiàng)添加到二維非穩(wěn)態(tài)對(duì)流換熱方程中,即完成流場(chǎng)溫度場(chǎng)的耦合。

圖3 泄漏點(diǎn)下方壓力分布圖Fig.3 Pressure distribution below the leak point

圖4 速度分布矢量圖Fig.4 Speed distribution vetor

(12)

對(duì)式(12)進(jìn)行離散,采用中心差分形式,對(duì)角點(diǎn)、邊界點(diǎn)和內(nèi)點(diǎn)分開離散,再用Gauss-Seidel迭代法計(jì)算內(nèi)壁絕熱層區(qū)域內(nèi)的內(nèi)點(diǎn),得到LNG液體泄漏后在罐壁內(nèi)的溫度分布(泄漏點(diǎn)i=1 000)。

由圖5的溫度分布云圖可知,隨著計(jì)算時(shí)間的推進(jìn),泄漏液體向下擴(kuò)散程度大,且越靠近泄漏點(diǎn)處泄漏速度越大,溫度變化梯度也隨之增大。

3 Fluent軟件模擬驗(yàn)證

Fluent軟件是采用有限容積法對(duì)控制方程進(jìn)行離散并求解的,離散方程有清晰的物理解釋,其內(nèi)置各類模型及方法可用于求解不同流體力學(xué)問題,針對(duì)滲流問題,亦采用多孔介質(zhì)模型進(jìn)行模擬。

Fluent軟件中涉及到傳熱問題時(shí),其基本守恒方程為質(zhì)量守恒、動(dòng)量守恒以及能量守恒方程[17]。多孔介質(zhì)模型本質(zhì)上其實(shí)僅在動(dòng)量方程上附加源項(xiàng);若涉及傳熱問題,可以只修改能量方程的傳導(dǎo)通量及瞬態(tài)項(xiàng)[17]。實(shí)際上,Fluent軟件并不存在獨(dú)立的多孔介質(zhì)區(qū)域,模型假設(shè)多孔介質(zhì)固體區(qū)域?yàn)榱黧w區(qū)域,對(duì)流入其中的流體設(shè)置阻力系數(shù)以達(dá)到滲流效果[18];如涉及熱傳導(dǎo),設(shè)置多孔區(qū)域孔隙率,通過孔隙率換算出有效熱導(dǎo)率,將其添加于能量方程中以達(dá)到熱量傳遞效果[19]。計(jì)算得出絕熱材料的阻力系數(shù)見表3。

使用ANSYS workbench工作平臺(tái)建立二維模型,合理劃分網(wǎng)格后,將參數(shù)條件導(dǎo)入Fluent軟件中進(jìn)行模擬,在上述相同的溫度邊界條件下,通過穩(wěn)態(tài)求解的方法得到罐壁的溫度場(chǎng)分布，見圖6。

a)100時(shí)步a)100 time steps

b)500時(shí)步b)500 time steps

c) 1 000 時(shí)步c)1 000 time steps

d) 1 000 時(shí)步放大d)1 000 time steps amplification

表3 多孔介質(zhì)模型各絕熱材料參數(shù)表

Tab.3 Insulation materials parameters in porous model

系數(shù)膨脹珍珠巖各相同性材料彈性氈x方向(徑向)y方向(垂直)孔隙率0.080.080.08黏性阻力系數(shù)/m-21.413×10112.5×10115×1010慣性阻力系數(shù)/m-14.746×1063.5×1067×105

圖6 罐體穩(wěn)態(tài)溫度分布圖Fig.6 Temperature distribution of the tank in steady state

再模擬泄漏工況,泄漏點(diǎn)為距儲(chǔ)罐底部10 m處,泄漏孔直徑10 mm,泄漏熱溫度條件與上述一致。LNG從泄漏口流出,采用多孔介質(zhì)模型,按表3設(shè)置模擬參數(shù),獲得不同時(shí)刻罐體溫度分布,結(jié)果見圖7。

圖7 LNG儲(chǔ)罐泄漏工況下溫度分布圖Fig.7 Temperature distribution diagram of LNG tank under leakage conditions

由圖7可以看出,內(nèi)罐泄漏時(shí),漏液先向下聚集,靠近內(nèi)罐壁的彈性氈快速降溫至LNG溫度,當(dāng)?shù)蜏匾后w到達(dá)罐底時(shí),低溫LNG液體富集底部,而后沿罐底逐漸滲透到膨脹珍珠巖中,最后擴(kuò)散到整個(gè)絕熱層。經(jīng)計(jì)算,泄漏后約26 h低溫液體滲透至全部絕熱層中。

4 結(jié)論

通過Fluent軟件模型和CFD計(jì)算得到了相似的結(jié)果,在計(jì)算時(shí)長和邊界條件的設(shè)置難度上,Fluent模型要優(yōu)于CFD直接計(jì)算,所以如果有模型符合應(yīng)用場(chǎng)景時(shí),可以選取模型進(jìn)行計(jì)算。

LNG內(nèi)罐發(fā)生泄漏時(shí),漏液先向下聚集,呈不對(duì)稱分布,溫度分布表現(xiàn)為靠近泄漏源處的彈性氈首先降溫到LNG溫度,且此處的溫度梯度較大;而后沿罐底逐漸滲透到膨脹珍珠巖中,隨著泄漏時(shí)間的不斷推進(jìn),漏液最后擴(kuò)散至整個(gè)內(nèi)壁絕熱層。但由于次容器與泡沫玻璃磚的存在,外罐壁的溫度與外界溫度相差不到1 ℃,起到了良好的絕熱保冷效果。