LNG儲罐內(nèi)罐小孔泄漏工況下溫度場數(shù)值模擬

張 娜 李兆慈 郭志超

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

0 前言

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

1 LNG儲罐泄漏研究方法

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

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

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

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

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

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

2 內(nèi)罐泄漏溫度場計算模擬

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

在LNG儲罐內(nèi)罐發(fā)生泄漏時,LNG會沿有限孔緩慢流出,再沿著儲罐內(nèi)壁外的彈性玻璃纖維氈流入環(huán)形空間的絕熱層中。絕熱層由彈性氈、膨脹珍珠巖和泡沫玻璃磚組合而成[14]。以 20 000 m3雙金屬全容式LNG儲罐為例,內(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

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

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

表2 儲罐罐壁絕熱層數(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儲罐內(nèi)罐泄漏模擬可采用多孔介質(zhì)模型[15]。

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

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

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

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

(1)

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

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

(2)

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

(3)

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

再計算普朗特數(shù):

(4)

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

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

(5)

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

(6)

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

對于二維平面:

(7)

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

φ=φa+Cf(p-pa)

(8)

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

對于彈性液體:

ρ=ρaeCL(p-pa)=ρa[1+CL(p-pa)]

(9)

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

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

(10)

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

(11)

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

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

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

圖3 泄漏點下方壓力分布圖Fig.3 Pressure distribution below the leak point

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

(12)

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

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

3 Fluent軟件模擬驗證

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

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

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

a)100時步a)100 time steps

b)500時步b)500 time steps

c) 1 000 時步c)1 000 time steps

d) 1 000 時步放大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

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

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

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

4 結(jié)論

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

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