液化天然氣泄漏和水面擴散過程模擬

何思念，常華偉，文科，舒水明

（華中科技大學(xué)能源與動力工程學(xué)院，湖北 武漢430074）

引 言

LNG（液化天然氣）在天然氣產(chǎn)業(yè)和能源市場中的地位日益顯著，僅2013年我國進口1700萬噸LNG，而2020年則預(yù)期進口量將會達到6000萬噸。LNG運輸量的增加導(dǎo)致水上泄漏的可能性也越來越大，為了避免LNG泄漏導(dǎo)致的危害，對LNG泄漏過程特性的研究也越來越受到關(guān)注。當?shù)蜏匾后w罐破裂且低溫液體流出時，巨大的溫差會引發(fā)快速相變 （RPT）爆炸，從而對周邊設(shè)施造成損害。即使沒有發(fā)生爆炸，擴散到水面或地面的低溫液體會從環(huán)境中吸收大量的熱，導(dǎo)致附近人員的凍傷，并且產(chǎn)生大量白霧。如果所泄漏的為可燃性物質(zhì)，如LNG或液化氫 （LH2），散發(fā)到空氣中的氣體可能被引燃乃至發(fā)生爆炸。LNG的泄漏過程危害評估主要包括泄漏、蒸發(fā)、擴展、氣體擴散和燃燒等過程［1－2］。其中，泄漏、蒸發(fā)和擴散的模擬結(jié)果是模擬可燃云擴散和建立池火模型的重要的前提條件。為了進一步提升LNG泄漏的預(yù)測模型的準確性，已經(jīng)有大量的理論工作用于研究LNG在水面泄漏和擴散的可能過程和現(xiàn)象［3－5］。

儲存在運輸船中的LNG處于一個大氣壓的飽和狀態(tài)［6］，其沸點約為－162℃，這使得LNG的儲運需要采用專門的容器。儲罐的結(jié)構(gòu)參數(shù)、裝載量及穿破位置及其孔徑都會影響到泄漏總體積、持續(xù)時間及流量變化過程。漏洞孔徑極大時的泄漏持續(xù)時間極短，這一過程可以視為瞬態(tài)泄漏；漏洞孔徑極小時的泄漏流量幾乎不隨時間而變化，這一過程可以視為恒速泄漏；其他情形下的則為變流速泄漏過程。由于LNG泄漏情形的變化，LNG在水面的蒸發(fā)和擴散現(xiàn)象也會有很大不同。自罐中泄漏出的LNG會在水面上形成LNG液池，該池會從周圍環(huán)境吸收熱量并汽化成天然氣蒸氣。這些低溫氣體會形成低臥云并隨風(fēng)而飄移。當天然氣在空氣中的濃度達到可燃濃度的上下限之間并與火源接觸，則可能發(fā)生劇烈燃燒或者爆炸。對于正在燃燒的LNG低溫液池，低溫池吸收的很大一部分熱量來自于火焰的輻射傳熱；如果沒有被引燃，則最主要的熱量來自水中。同時空氣被LNG池及低溫氣體冷卻，其溫度降低相對濕度增加，彌漫在空氣中的水蒸氣則會冷凝后產(chǎn)生大量白霧。

當只有少量LNG泄漏到水上，LNG將迅速汽化消失；反之，泄漏量大時，會在水面上形成低溫池，池徑與池深及總蒸發(fā)速率不斷變化。起初LNG池迅速擴張，同時由于水的阻力而導(dǎo)致其擴張速度降低。對于恒速泄漏，當液池的大小與形狀保持相對穩(wěn)定且泄漏率等于總蒸發(fā)率時，液池便達到了平衡狀態(tài)；對于瞬態(tài)泄漏，液池將會擴張到最大尺寸，之后便開始收縮或者形成環(huán)形液池；對于變流速泄漏，液池的大小始終會隨泄漏率的變化而變化。液池形狀還取決于堤岸與船體：如果LNG從運輸船或者岸邊泄漏到開闊的水域時，由于有船體和堤岸的阻擋，液面在水面的形狀就近似為半圓；如果沒有阻擋來限制液池擴展，則其形狀將接近圓形。由于重力及蒸發(fā)作用，池深在從源點向半徑擴展的方向上總體呈現(xiàn)出減小的趨勢，但不會低于一個最小值，這一最小值由表面張力決定。

關(guān)于低溫液體的泄漏的研究，目前多為國外學(xué)者的工作，國內(nèi)科研人員在相關(guān)方面研究為數(shù)不多［7－9］。

早在1972年，當石油被廣泛運輸時，F(xiàn)annelop和Waldman即提出淺水理論來研究原油的泄漏與擴散過程［10］。該模型將液池簡化為圓柱體，其半徑的增加由下式給出

式中，Δ＝ρLNG／ρw。

由于這一式子中只考慮了重力，因而這一模型被稱為重力模型。Fay成功地將這一模型應(yīng)用于LNG擴散過程的模擬中［11］，此后這一模型得到了廣泛的發(fā)展和使用［12－13］。為了將擴散過程中的阻力考慮在內(nèi)，Webber［14］發(fā)展出了基于淺水理論的平均深度模型。Webber的模型將池深描述為積分形式，因而在本文中將其簡稱為積分模型。積分模型的數(shù)值求解過程相對簡單，具有良好的理論基礎(chǔ)，它得到了廣泛的使用與推廣［15－17］。雖然積分模型在LNG泄漏過程的仿真建模中使用廣泛，但其缺點在于，使用平均高度來簡化液池的高度曲線后，池深的不均勻性被忽略了。為了更準確地描述池深及速度變化，改進后的淺水模型采用微分形式來表示質(zhì)量守恒與動量守恒定律［18］。為了便于描述，在本文中改進后的淺水模型被稱為微分模型。

從已有的文獻資料可以看出，雖然已經(jīng)有了很多學(xué)者關(guān)注于LNG泄漏相關(guān)研究，但泄漏時及其隨后的現(xiàn)象的很多細節(jié)仍需要大力研究。首先，幾乎所有文獻中的模擬方案都采用恒速泄漏或瞬態(tài)泄漏，對泄漏細節(jié)描述得很少。其次，僅有的考慮了變化的泄漏速率的數(shù)據(jù)都是基于立方形儲罐的。本文將模擬LNG儲罐水面泄漏的情形，并給出更多相關(guān)細節(jié)，為此建立了一種基于微分模型并整合了各種儲罐的計算方法。首先采用使用該計算方法來計算不同孔徑下的泄漏與擴散過程，然后將MOSS儲罐與立方形儲罐二者泄漏后的擴散過程相比較，最后還分析討論了LNG液池的形狀。

1 數(shù)學(xué)模型和計算方法

1.1 LNG儲罐結(jié)構(gòu)

目前的LNG運輸船主要有兩種類型：MOSS和GTT儲罐。MOSS型LNG儲罐是標準球形，由鋁合金制成。如圖1（a）所示，它采用了立式提升技術(shù) （VST）來擴展容量。忽略泵塔與其他一些內(nèi)部結(jié)構(gòu)的體積，MOSS儲罐的橫截面積可以表示為

GTT膜儲罐的橫截面通常為八角形，如圖1（b）所示。忽略泵塔與其他一些內(nèi)部結(jié)構(gòu)的體積，其橫截面積可以表示為

圖1 LNG儲罐的結(jié)構(gòu)Fig.1 Construction of LNG tanks

如果簡化為立方形儲罐，其橫截面積At則為常數(shù)。

1.2 LNG泄漏過程

當儲罐被穿破后，LNG在重力作用下會泄漏出來，其流率為

其中Cd的推薦值在0.6～0.7之間［19］。

LNG液面下降速率等于源流率除以橫截面積：

1.3 LNG水面擴散過程

LNG在水面擴散時，其圓形的擴散行為可以用與質(zhì)量守恒和動量守恒相關(guān)的微分方程來描述。

連續(xù)性方程：

動量守恒方程：

式中，r為徑向坐標；h和u分別為當?shù)爻厣詈推骄疃葟较蛩俣?；v為液體蒸發(fā)率，v＝q／ρLiL；w為LNG向水中泄漏的流率，它與總泄漏流量Q的關(guān)系為

1.4 模擬策略

圖2展示了對LNG泄漏過程進行模擬分析的策略。首先是選擇低溫液體的種類，包括LNG、液氫和液氮等低溫液體。然后是選擇泄漏模型，這里共有3種模型：瞬態(tài)泄漏、恒速泄漏、儲罐泄漏。根據(jù)選擇會有3種不同的計算模型分別運行。其中瞬態(tài)泄漏只需設(shè)置泄漏總量。在恒速模型中需要將泄漏速率設(shè)置為常量。儲罐泄漏模型中的泄漏速率根據(jù)儲罐類型和裝載量以及泄漏孔徑來計算得出。在這3種模型中，儲罐模型是最復(fù)雜也是最準確的。接下來設(shè)置泄漏模型的技術(shù)參數(shù)、初始條件和液池形狀數(shù)。計算完離散方程后，輸出不同時間步長的h和u還有高度剖面。還可以輸出液池半徑隨時間變化關(guān)系曲線。

圖2 模擬流程圖Fig.2 Framework chart of simulation

1.5 數(shù)值計算方法

對微分方程進行求解，采用一階有限差分顯式離散格式和交錯網(wǎng)格法。將時間步長設(shè)置為1ms，網(wǎng)格間距為0.05m，通過C＋＋語言編程來進行求解。低溫液體的熱力學(xué)參數(shù)由REFPROP提供。本文采用Verfondern等［20］的液氫在水池中擴散蒸發(fā)試驗驗證差分淺水方程對預(yù)測低溫液體在水面擴散問題的可行性。

圖3是5L·s－1的速率泄漏的LH2在水面擴散液池的半徑隨著時間變化情況。從圖中可以看出，LH2一旦泄漏后就迅速擴散。在達到峰值后，液池就處于動態(tài)平衡中。圖3中所顯示出的實驗與仿真結(jié)果的吻合說明了此方法的可行性。

2 模擬結(jié)果與討論

2.1 LNG水面泄漏的基本參數(shù)

假設(shè)有12500m3LNG從立方形儲罐中通過直徑1m的孔洞泄漏出來，并在水上形成一個半圓形LNG池。

圖3 液氫擴散實驗和模擬結(jié)果對比圖Fig.3 Validation of numerical results for spreading of liquid hydrogen with experimental results

LNG性質(zhì)：

成分：甲烷

密度：422.5kg·m－3

潛熱：510.82kJ·kg－1泄漏參數(shù)：

容積：25000m3

排放口系數(shù)：0.65

儲罐類型：立方體

橫截面積：961.5m2

總泄漏量：12500m3

孔徑：1m

孔以上液體初始高度：13m

液池形狀：半圓形

LNG水面換熱熱通量：85kW·m－2

圖4顯示了LNG池在泄漏伊始迅速擴散，之后液池開始萎縮至全部蒸發(fā)，在第341秒時池徑達到峰值108m。通過方程 （1）知，當LNG泄漏到水面上時，其在重力作用下迅速橫向擴散。圖4中的趨勢線清楚地描述了隨著池深逐步減小，LNG的擴散速率也在降低。當LNG池達到最大時，泄漏流量與蒸發(fā)流量二者近似平衡：

隨著儲罐內(nèi)的液面越來越低，LNG的泄漏流量也越來越小，液池半徑也就隨著時間逐漸減小，直至完全消失。

圖5中顯示了計算所得的最大LNG池的池深及平均高度曲線。開始時池深迅速減小，而后穩(wěn)定在0.1m左右。在R＝30m處存在一個波峰，之后池深隨著半徑的增大不斷減小到0。深度平均值為0.0873m，其標準偏差為0.1005m。當LNG在水上擴散時，由于蒸發(fā)與重力作用，其深度通常會變小但是當其速度大于當?shù)厣疃忍幍闹亓Σǖ乃俣葧r，會發(fā)生水躍，使得高度截面突然上升［21］，影響液池蒸發(fā)和擴散過程。

圖4 12500m3 LNG從孔徑為1m的儲罐中泄漏后水面擴散半徑計算結(jié)果Fig.4 Calculated result of pool radius for 12500m3 LNG releasing from a 1mdiameter breach hole

圖5 LNG液池在341s時的高度曲線Fig.5 Height profile of LNG pool at 341s

2.2 漏洞孔徑對LNG擴散過程的影響

改變基本情形下的泄漏口直徑，保持其他條件不變，重新計算這一系列的泄漏過程。從圖6中能夠看出，隨著漏洞孔徑增加，微分模型計算的最大池徑與積分模型的計算結(jié)果有著相同的變化趨勢，最大半徑跟漏洞孔徑呈現(xiàn)相似的Boltzmann非線性關(guān)系，即隨著孔徑增加，最大液池尺寸也急劇增加，當泄漏孔徑達到5m時，液池半徑穩(wěn)定至一個漸進值。但這兩個模型結(jié)果的漸近值不同，微分模型的漸近值約為265m，而積分模型的漸近值約為350m。Qiao等［19］的長期瞬態(tài)泄漏理論很好地解釋了這一現(xiàn)象。然而從圖中可以看出，采用積分模型和采用微分模型得到的漸近線之間存在差異?？梢酝茢嘁好娓叨惹€的不均勻性導(dǎo)致了微分模型和平均深度積分模型計算結(jié)果之間明顯的差異。

圖6 分別采用微分模型和積分模型模擬LNG從不同孔徑泄漏過程Fig.6 Calculation of LNG pool behavior of LNG spilling from various hole of different size with intergral model and differential model.

2.3 LNG從MOSS儲罐中泄漏的擴散過程

一個標準LNG船的載運量為1×105～2×105m3，每個船內(nèi)含4～6個儲罐。這里研究的是內(nèi)置圓柱型體積為25000m3的MOSS儲罐，其主要的技術(shù)結(jié)構(gòu)參數(shù)列在表1中。除了儲罐的技術(shù)結(jié)構(gòu)參數(shù)外，其余所有計算條件都與基本情形保持一致。表2表示了MOSS儲罐與立方形儲罐的LNG泄漏后的擴散過程特性。從表2所顯示的結(jié)果中可以看出，MOSS儲罐的最大池半徑比立方形儲罐大1m，但池存在周期則少14s。

表1 容積為25000m3的MOSS LNG儲罐的主要技術(shù)參數(shù)Table 1 Main technical parameters of 25000m3 MOSS LNG carrier

表2 LNG從MOSS儲罐和立方體儲罐泄漏的模擬結(jié)果Table 2 Computated results for LNG release from MOSS spherical tank and cubic tank

基于式 （5）可知，當液壓頭相同時，儲罐橫截面積越小則液面降低越快。由于MOSS儲罐的橫截面積隨深度而變化，其LNG泄漏率的變化與立方形儲罐的泄漏率變化不一致，流量的不同變化會影響到LNG在水面的擴散過程。不過，從表2可以看出，MOSS儲罐和立方形儲罐的泄漏過程的差別并不是很大，其中影響最大的是達到最大擴散半徑的時間，其相差也不過2%，而最大半徑以及池壽命相差都在1%之內(nèi)。由此可以得出，在要求不高的情形下，用立方形儲罐替代其他不規(guī)則形狀的儲罐也是可行的。

2.4 液池形狀對LNG擴散的影響

如圖7所示，因船身的阻擋，LNG池的形狀在圓形與半圓形之間變化。一般情況下，液池形狀總體上近似為半圓形，少部分情形下則近似于圓形。表3表征了分別采用半圓形和圓形液池來模擬液池擴散過程的結(jié)果。結(jié)果表明，半圓形液池的最大半徑是圓形液池的1.385倍 （約等于），這說明兩種形狀下的液池最大面積幾乎相同。從表中還可以看出，液池的生命周期也幾乎相同，但由于船體的阻擋，液池擴散速度相對較慢，這就導(dǎo)致使用半圓形液池達到最大直徑的時間比圓形的晚了71s。

圖7 LNG在水面擴散俯視圖Fig.7 Top view of LNG spreading on water

表3 采用圓形和半圓形液池的計算結(jié)果Table 3 Computational results with semicircle and circular pool

3 結(jié) 論

LNG運輸船泄漏是一個比較危險而且較復(fù)雜的物理過程。在眾多的危害預(yù)測模型中，淺水微分模型相對更復(fù)雜和準確。本文選擇淺水微分模型對泄漏過程中的各種影響因素進行了評估，并對比了積分模型和微分模型的結(jié)果。結(jié)果表明，在LNG擴散的過程中，液池會形成階躍，由此導(dǎo)致液池高度曲線與平均值的離線度達到1.15。由于考慮了液池高度的不均勻，微分模型計算的液池半徑比積分模型的要小。但兩者的泄漏孔徑的影響趨勢是嚴格一致的Boltzmann非線性關(guān)系。通過計算MOSS儲罐發(fā)現(xiàn)，采用標準立方體儲罐來替代MOSS等不規(guī)則儲罐來對LNG的泄漏擴散過程進行模擬是合理的。最后，對比圓形和半圓形液池擴散過程發(fā)現(xiàn)，船體的阻擋會減緩擴散過程，但對液池的最大擴散面積的幾乎沒有影響。

符 號 說 明

Ap——液池面積，m2

At——LNG上表面面積，m2

a0～a2，b0，b1，l——GTT儲罐基本結(jié)構(gòu)參數(shù)，m

Cd——泄漏系數(shù)，1

d——距離，m

F——摩擦力，kg

g——重力加速度，m·s－2

Hh——泄漏孔高度，m

Hi——垂直加高高度，m

Ht——LNG在儲罐內(nèi)高度，m

h——液池高度，m

ˉh——平均高度，m

R——液池半徑，m

Rt——MOSS儲罐的半徑，m

r——半徑，m

Q——體積流率，m3·s

t——時間，s

u——速度，m·s－1

Vp——液池體積，m3

v——單位面積汽化速率，m·s－1

w——源項，m·s－1

Δ——相對密度，1

ρLNG——LNG密度，m3·kg－1

ρw——水的密度，m3·kg－1

［1］ Saraf S，Melhem G.Modeling LNG pool spreading and vaporization／／Proceedings of the AIChE Spring Meeting，ioMosaic Corporation ［C］.Atlanta，GA，2005.

［2］ Woodward J L，Pitbaldo R.LNG Risk Based Safety：Mode－ling and Consequence Analysis ［M］.New York：John Wiley ＆Sons，2010.

［3］ Fay J A.Model of spills and fires from LNG and oil tankers［J］.JournalofHazardousMaterials，2003，96 （2）：171－188.

［4］ Lehr W，Simecek－Beatty D.Comparison of hypothetical LNG and fuel oil fires on water ［J］.JournalofHazardous Materials，2004，107 （1）：3－9.

［5］ Pitblado R，Baik J，Hughes G，F(xiàn)erro C，Shaw S.Consequences of liquefied natural gas marine incidents ［J］.Process SafetyProgress，2005，24 （2）：108－114.

［6］ Spaulding M L，Swanson J C，Jayko K，Whittier N.An LNG release，transport，and fate model system for marine spills ［J］.JournalofHazardousMaterials，2007，140（3）：488－503.

［7］ Zhang Lijing （張禮敬），Ding Li（丁立），Tao Gang （陶剛）.Simulation of accident consequences of liquefied chlorine pipeline leakage ［J］.JournalofChemicalIndustryand Engineering（China）（化 工 學(xué) 報 ），2008，59 （5）：1333－1337.

［8］ Liu Zhongliang（劉中良），Luo Zhiyun （羅志云），Wang Jieteng（王皆騰），Zhang Jian （張建），Sun Xiaochun （孫曉春）.Analytical calculation of leakage rate for natural gas pipelines ［J］ .JournalofChemicalIndustryand Engineering（China）（化 工 學(xué) 報 ），2008，59 （8）：2121－2126.

［9］ Pan Xu （潘旭），Jiang Juncheng （蔣 軍 成）.Research on discharging dispersion of chemical dangerous gases and its influence factors ［J］.JournalofNanjingUniversityof ChemicalTechnology（南京化工大學(xué)學(xué)報），2001，23（1）：19－22.

［10］ Fannelop T K，Waldman G D.Dynamics of oil slicks ［J］.AmericanInstituteofAeronauticsandAstronautics Journal，1972，10 （4）：506－510.

［11］ Fay J A.Spread of large LNG pools on the sea ［J］.Journal ofHazardousMaterials，2007，140 （3）：541－551.

［12］ Ivings M，Lea C，Webber D，Jagger S，Coldrick S.A protocol for the evaluation of LNG vapour dispersion models［J］.JournalofLossPreventionintheProcessIndustries，2013，26 （1）：153－163.

［13］ Parihar A，Vergara C，Clutter J K.Methodology for consequence analysis of LNG releases at deepwater port facilities［J］.SafetyScience，2011，49 （5）：686－694.

［14］ Webber D.Source terms ［J］.JournalofLossPreventionin theProcessIndustries，1991，4 （1）：5－15.

［15］ Brambilla S，Manca D.Accidents involving liquids：a step ahead in modeling pool spreading，evaporation and burning［J］.JournalofHazardousMaterials，2009，161 （2／3）：1265－1280.

［16］ Subramani A，Jayanti S，Shet U S P，Selvaraj P.Dynamics of liquid sodium pool spreading under sodium fire conditions［J］.NuclearEngineeringandDesign，2009，239 （7）：1354－1361.

［17］ Khajehnajafi S，Pourdarvish R.Correlations for mass transfer from a liquid spill：comparisons and recommendations ［J］.ProcessSafetyProgress，2011，30 （2）：178－184.

［18］ Verfondern K，Dienhart B.Pool spreading and vaporization of liquid hydrogen ［J］.InternationalJournalofHydrogen Energy，2007，32 （13）：2106－2117.

［19］ Qiao Yuanhua，West H H，Mannan M S，Johnson D W，Cornwell J B.Assessment of the effects of release variables on the consequences of LNG spillage onto water using FERC models［J］.JournalofHazardousMaterials，2006，130（1／2）：155－162.

［20］ Verfondern K，Dienhart B.Experimental and theoretical investigation of liquid hydrogen pool spreading and vaporization［J］.InternationalJournalofHydrogenEnergy，1997，22（7）：649－660.

［21］ Raj P K.Evaporating liquid flow in a channel （an integral model based on shallow water flow approximation）［J］.JournalofLossPreventionintheProcessIndustries，2011，24（6）：886－899.