莊學(xué)強，孫 迪，高孝洪，李格升

(1.集美大學(xué) 輪機工程學(xué)院，福建 廈門 361021；2.福建省船舶與海洋工程重點實驗室，福建 廈門 361021；3.武漢理工大學(xué) 能源與動力工程學(xué)院，武漢 430070)

LNG船泄漏事故液池擴展計算及不確定性分析

莊學(xué)強1，2，孫 迪1，2，高孝洪3，李格升3

(1.集美大學(xué) 輪機工程學(xué)院，福建 廈門 361021；2.福建省船舶與海洋工程重點實驗室，福建 廈門 361021；3.武漢理工大學(xué) 能源與動力工程學(xué)院，武漢 430070)

根據(jù)小孔射流、慣性-質(zhì)量平衡擴展等理論搭建包括泄漏源強和液池半徑擴展在內(nèi)的LNG船泄漏液池擴展模型，分析影響液池擴展變化的重要參數(shù)及其不確定性，借助MATLAB編程計算這些不確定對預(yù)測結(jié)果的影響。研究結(jié)果表明：導(dǎo)致LNG船泄漏的事件不同，泄漏口的大小可能會有所不同(0～16 m2不等)，當(dāng)泄漏口面積<2 m2時，液池擴展對泄漏口大小的變化相當(dāng)敏感；受泄漏口形狀和流動狀態(tài)等因素影響，流量系數(shù)取值范圍較廣(0.4～1.0)，預(yù)測偏差可達(dá)2.5倍；因LNG-水?dāng)嚭铣潭燃胺序v形式不同，液池蒸發(fā)速率為0.042～0.25 kg/(m2·s)，預(yù)測的最大液池半徑差別較大。

水路運輸；液池；LNG船；泄漏

據(jù)統(tǒng)計,自2010年以來，我國液化天然氣(Liquefied Natural Gas,LNG)進(jìn)口量的年均增長率為30.5%。相關(guān)數(shù)據(jù)表明，到2020年末，我國很可能成為僅次于日本的全球第二大LNG進(jìn)口國。[1]LNG船作為天然氣供應(yīng)鏈中重要的海運裝備，承擔(dān)著遠(yuǎn)距離、大批量輸送天然氣的重任，直接關(guān)系到我國能源規(guī)劃的順利實施。然而,LNG船載有大量易燃易爆的液體，一旦因遭受碰撞、擱淺和恐怖襲擊等事故而泄漏，將導(dǎo)致災(zāi)難性后果。

LNG船舶泄漏事故的危害計算主要包括泄漏源強計算、液池擴展計算、可燃蒸氣云擴散計算、火災(zāi)爆炸危害計算等。[2]在這些計算中,前一步計算的結(jié)果是后一步計算的輸入,后一步計算的準(zhǔn)確度在很大程度上依賴于前一步的計算精度。目前研究最多的是LNG船泄漏蒸氣云擴散漂移，得到許多計算模型。但是,對LNG船泄漏液池擴展的研究較少，主要原因在于LNG船泄漏液池擴展是一個復(fù)雜的傳質(zhì)與傳熱耦合過程，影響因素有很多，存在諸多不確定性。

1 LNG船泄漏液池擴展計算模型建立

LNG船發(fā)生泄漏事故時，推動液池在水面擴展的主動力是泄漏入水的LNG，泄漏源強是液池擴展計算的必要輸入?yún)?shù)，因此一個完整的LNG船泄漏液池擴展計算模型應(yīng)包括泄漏源強子模型和液池擴展子模型。下面根據(jù)小孔射流、質(zhì)量平衡及慣性擴展理論建立微分方程。

1.1LNG船泄漏源強計算模型

定義Δt為時間步長，i為時間步數(shù)。根據(jù)伯努利方程可得LNG船泄漏速率qi為

(1)

(2)

ΔQ=qiΔt

(3)

Qi=Qi-1+ΔQ

(4)

式(1)～式(4)中：Chole為泄漏孔的流量系數(shù)；Ahole為泄漏孔的流通面積，m2；PT為LNG貨艙液面上方的壓力，Pa；Patm為大氣壓力，101 325 Pa；ρL為LNG液體密度，kg/m3；hi為第i時間步時LNG液面距離泄漏孔中心線的長度，m；ΔQ為第i時間步內(nèi)的泄漏量,kg；Qi為總泄漏量,kg；AT，i-1為第i-1時間步對應(yīng)的液面橫截面積，m2。

1.2LNG船泄漏液池擴展模型

LNG船泄漏液池半徑擴展變化過程的計算式可表示為

Mi=Mi-1+qiΔt-Mevap

(5)

(6)

Vi=Mi/ρl

(7)

(8)

1) 若δi>δmin,則液池半徑擴展速度[3]為

(9)

(10)

2) 若δi≤δmin,則液池半徑為

(11)

式(5)～式(11)中：Ri為迭代步為i時LNG液池的半徑，m；Ks為擴展常數(shù)，取1.414；δi為迭代步為i時LNG液池的厚度，m；δmin為LNG液池最小穩(wěn)定厚度，m；g為重力加速度,9.8 m/s2；Mi為迭代步為i時LNG液池的質(zhì)量，kg；Vi為迭代步為i時LNG液池的體積，m3；Mevap為Δt時間步長內(nèi)的蒸發(fā)量，kg；mi為Δt時間步長內(nèi)LNG泄漏的速率,kg/s;mevap為單位面積內(nèi)的蒸發(fā)速率，kg/(m2·s)。

2 LNG船泄漏液池擴展不確定性分析

以上主要根據(jù)小孔射流、質(zhì)量平衡及慣性擴展理論建立LNG船泄漏液池擴展模型，國外許多著名機構(gòu)(如FAY,QUEST和ABS等)也根據(jù)類似理論建立相應(yīng)的模型。然而,LNG船泄漏及液池在水面擴展是一個非常復(fù)雜的過程，許多影響因素(如泄漏口大小、流量系數(shù)、液池蒸發(fā)速率及儲罐內(nèi)的壓力等)在突發(fā)事故中難以確定，因此采用目前的模型對LNG船泄漏液池進(jìn)行擴展模擬尚存在一定的不確定性。FAY,QUEST和ABS曾針對同一假定事故分別進(jìn)行模擬計算，得出的結(jié)果相差4倍多。[4]

2.1泄漏口大小的不確定性

對于LNG船而言,無論是碰撞、擱淺等操作性事故導(dǎo)致的泄漏，還是受到恐怖襲擊導(dǎo)致的泄漏，泄漏口大小均難以確定。美國Sandia國家實驗室[5]在建立一系列LNG船的三維有限元結(jié)構(gòu)模型分析各種船舶事故和恐怖襲擊可能造成的損害之后認(rèn)為：由碰撞、擱淺等事故造成的LNG液艙破損面積≤2 m2；由恐怖襲擊造成的LNG液艙破損面積為2～16 m2。由各行業(yè)組成的LNG危害評估小組[6]在研究之后認(rèn)為：由碰撞、擱淺等事故導(dǎo)致LNG船舶泄漏口的寬度為0～1.5 m，建議在進(jìn)行風(fēng)險評估時取平均值0.75 m作為最大可信泄漏孔徑。挪威船級社(Det Norske Veritas,DNV)認(rèn)為LNG船遭受恐怖襲擊的最大可信泄漏孔徑為1.5 m。由此可見,當(dāng)前在泄漏口大小的確定方面還存在一定的困難。

為定量分析泄漏口大小的不確定性對LNG船泄漏液池擴展的影響，以某138 000 m3LNG實體船為例，設(shè)定不同的泄漏口大小，結(jié)合已建立的液池擴展模型加以分析比較。該LNG實體船的貨艙為薄膜型貨艙，總艙容為36 650 m3，其剖面結(jié)構(gòu)尺寸見圖1。假設(shè)該船因故發(fā)生泄漏事故，泄漏口面積為0.5～16.0 m2不等，泄漏孔中心與液面的初始距離為12.35 m，LNG的密度為442.5 kg/m3，液池蒸發(fā)速率為0.2 kg/(m2·s)，最小液池厚度為5.2 mm。應(yīng)用MATLAB軟件對液池擴展模型進(jìn)行編程計算，計算結(jié)果見圖2。

圖1 某薄膜型LNG船貨艙橫剖面尺寸

圖2 不同泄漏口面積的液池擴展情況比較

由圖2可知，泄漏口的面積不同，LNG液池的擴展情況(主要指液池最大半徑和存在時間)是不同的。液池擴展情況對不同泄漏口面積呈現(xiàn)出的敏感度也是不同的:當(dāng)泄漏口面積<2 m2時，液池擴展對泄漏口大小的變化相當(dāng)敏感;而當(dāng)泄漏口面積>7 m2時，液池擴展隨著泄漏口大小的變化相對不明顯。例如:當(dāng)泄漏口面積為16 m2時，液池最大半徑和存在時間分別為292.6 m及436 s;當(dāng)泄漏口面積為7 m2時，液池最大半徑和存在時間分別為284.4 m及463.5 s;當(dāng)泄漏口面積為0.5 m2時，液池最大半徑和存在時間分別為97.6 m及4 676 s;當(dāng)泄漏口面積為2 m2時，液池最大半徑和存在時間分別為188.6 m及1 263 s。導(dǎo)致船舶泄漏事故發(fā)生的主要原因仍是碰撞和擱淺，這就意味著LNG船發(fā)生泄漏事故時，其泄漏口面積通常仍<2 m2，即液池擴展變化對泄漏口面積的變化很敏感。因此,在無法確定泄漏口大小的情況下，LNG船泄漏液池擴展計算有很大的不確定性。

2.2流量系數(shù)的不確定性

流量系數(shù)用來表征流體流過孔洞時受到的摩擦損失的大小。在對LNG船泄漏事故風(fēng)險進(jìn)行評價時，不同機構(gòu)采用的流量系數(shù)不盡相同(例如：DNV和Sandia取0.6;FAY和ABS取1.0)[7]，因此計算得到的危害后果差別很大。實際上，LNG船泄漏事故中流量系數(shù)很難確定,這是因為：

1) 泄漏口的形狀往往難以確定。泄漏口的形狀不同，其流量系數(shù)的建議取值是不同的。表1為不同泄漏口形狀液體泄漏流量系數(shù)建議值。[8]

表1 不同泄漏口形狀液體泄漏流量系數(shù)建議值

2) LNG船的液貨在泄漏到外部環(huán)境之前會經(jīng)過內(nèi)薄膜、絕熱層、外薄膜和船舶的雙層殼等多層阻礙，流經(jīng)距離至少2 m，正常達(dá)到3～4 m，整個流動過程十分復(fù)雜，流量系數(shù)難以確定。

3) 根據(jù)SOTERIOU等[9]的研究結(jié)論，LNG泄漏射流的流量系數(shù)主要由泄漏孔內(nèi)部的流動狀態(tài)(單相流動、空穴流動和返流流動)來決定,而泄漏孔內(nèi)部的流動狀態(tài)主要受泄漏孔直徑、長度、上游壓力、下游壓力、液體黏度、密度及飽和蒸汽壓等因素的影響。

圖3為射流流動狀態(tài)影響因素。在泄漏過程中，隨著LNG儲罐自由液面高度不斷下降，泄漏孔上游壓力不斷下降，泄漏孔內(nèi)部的流動狀態(tài)將發(fā)生變化。此外,流量系數(shù)在整個泄漏過程也是不斷變化的。以上述LNG船為例，假設(shè)泄漏面積為2 m2，流量系數(shù)分別取0.4,0.6和1.0，其他條件不變，計算得到液池擴展情況比較見圖4。

圖3 射流流動狀態(tài)影響因素

圖4 不同流量系數(shù)的液池擴展情況比較

由圖4可知：流量系數(shù)的取值對液池擴展變化的影響較大。在其他條件相同的情況下，當(dāng)流量系數(shù)取0.4時，液池的最大半徑和存在時間分別為123.3 m及2 965 s;而當(dāng)流量系數(shù)取1.0時，液池最大半徑和存在時間分別為188.7 m及1 263 s。由此可見，在LNG船泄漏流量系數(shù)難以確定的情況下，液池擴展變化計算也存在著不確定性。

2.3蒸發(fā)速率的不確定性

LNG液池蒸發(fā)速率的大小主要取決于水與LNG液池間的熱通量。根據(jù)液體與下墊面間的溫度差，液體的沸騰形式有核沸騰(Nucleate Boiling)、過渡沸騰(Transition Boiling)和膜沸騰(Film Boiling)，沸騰的狀態(tài)不同，熱通量是不同的(見圖5)。研究[3]表明，泄漏入水的LNG在液池內(nèi)大部分時間都處于膜沸騰狀態(tài)，只有當(dāng)液池蒸發(fā)進(jìn)入后期 (蒸發(fā)量達(dá)69%～84%)時才轉(zhuǎn)為過渡沸騰狀態(tài)(見圖5中的虛箭頭)。因此，許多蒸發(fā)模型經(jīng)常將膜沸騰熱通量28 kW/m2作為LNG液池與水的熱通量，與之對應(yīng)的蒸發(fā)速率為0.042 kg/(m2·s)。該蒸發(fā)速率是假定泄漏入水的LNG對水面沒有擾動的情況下得到的理論值，而實際上泄漏入水的LNG或多或少都會對水面造成擾動，LNG與水的攪合將使熱交換急劇增加，蒸發(fā)速率也隨之大大提高。表2[10]為不同泄放(攪合)方式下LNG蒸發(fā)速率試驗數(shù)據(jù)統(tǒng)計,可看出泄放的方式不同，液池蒸發(fā)速率有較大差別。然而,表2中的數(shù)據(jù)僅僅是通過小尺度試驗得到的，對于真正的LNG船泄漏事故，泄漏到水面的LNG速度可高達(dá)15 m/s,射流可至水面以下13 m處[11]，LNG與水?dāng)嚭陷^為激烈，其蒸發(fā)速率更加難以確定。正因為如此，不同機構(gòu)推薦的LNG水面蒸發(fā)率是不同的。例如:DNV使用的PHAST模型假設(shè)蒸發(fā)率為0.182 kg/(m2·s)；ABS假設(shè)蒸發(fā)率為0.072 kg/(m2·s)；QUEST軟件假設(shè)蒸發(fā)率為0.2 kg/(m2·s)[12]。圖6為泄漏面積為2 m2，流量系數(shù)為1.0，蒸發(fā)速率分別為0.042 kg/(m2·s)(理論值),0.072 kg/(m2·s),0.182 kg/(m2·s)和0.200 kg/(m2·s)時，LNG船泄漏液池擴展半徑計算結(jié)果比較。由圖6可知:按理論蒸發(fā)速率計算會過高地估算最大液池半徑和存在時間；采用不同機構(gòu)不同的蒸發(fā)率取值，計算結(jié)果將有較大差別(例如,按ABS蒸發(fā)率計算出液池最大半徑和存在時間為293.7 m及1 532 s，而按 QUEST蒸發(fā)率計算出液池最大半徑和存在時間為188.7 m及1 263 s)。蒸發(fā)率不確定將直接導(dǎo)致液池擴展計算不確定。

圖5 典型的沸騰-熱通量曲線

表2 不同泄放方式下LNG蒸發(fā)速率試驗數(shù)據(jù)統(tǒng)計

圖6 不同蒸發(fā)速率的液池擴展情況比較

2.4貨艙壓力不確定分析

目前許多模型都假設(shè)在LNG船泄漏過程中作用于貨艙LNG液面上的壓力保持不變。實際上，LNG貨艙是相對封閉的空間，艙內(nèi)壓力的變化主要取決于LNG泄漏體積率、真空破壞閥補氣率和LNG體積蒸發(fā)率三者間的相互作用。當(dāng)LNG船正常卸貨或泄漏孔洞較小時，卸貨/泄漏體積率等于補氣閥補氣率與LNG體積蒸發(fā)率之和，艙內(nèi)壓力基本上保持不變;但當(dāng)泄漏孔洞較大時，尤其是在泄漏初期，LNG泄漏體積率明顯大于貨艙的真空破壞閥補氣率與貨艙內(nèi)LNG的蒸發(fā)率之和，貨艙內(nèi)作用于LNG液面的壓力顯著下降，出現(xiàn)較高的真空度，艙內(nèi)壓力的下降反過來又會影響LNG泄漏速率和艙內(nèi)LNG的蒸發(fā)率，泄漏率和蒸發(fā)率改變又將使艙內(nèi)壓力發(fā)生新的變化，這種相互作用使得貨艙內(nèi)的壓力變得難以確定。為說明LNG船泄漏過程中貨艙內(nèi)壓力的變化情況，根據(jù)理想氣態(tài)方程和Hashemi-Wesson定理[13]建立LNG船泄漏時LNG泄漏體積率、真空破壞閥補氣率及貨艙內(nèi)LNG體積蒸發(fā)率相互作用微分方程組(限于篇幅此處不再給出)，以前述薄膜型貨艙為例，設(shè)定泄漏孔徑為4 m，求解貨艙內(nèi)壓力隨泄漏時間的變化過程，計算結(jié)果見圖7。

由圖7可知:在LNG泄漏過程中，艙內(nèi)壓力并不是穩(wěn)定不變的，而是會發(fā)生復(fù)雜變化，在泄漏剛開始時急劇下降，隨后逐漸回升;貨艙在整個泄漏過程中基本上都處于真空狀態(tài)。艙內(nèi)壓力的這種變化將使得平均泄漏速率下降，液池的最大半徑變小，液池的持續(xù)時間變長。

圖7 艙內(nèi)壓力隨泄漏時間變化過程

3 結(jié)束語

為對LNG船泄漏液池擴展變化進(jìn)行計算，根據(jù)小孔射流伯努利方程、質(zhì)量平衡及液體慣性擴展理論建立包括泄漏源強和液池半徑擴展在內(nèi)的LNG船泄漏液池擴展模型。由于缺乏事故統(tǒng)計數(shù)據(jù)及事故尺寸的現(xiàn)場試驗數(shù)據(jù)，模型中有些重要參數(shù)尚存在一定的不確定性，影響LNG船泄漏液池擴展計算。主要的不確定參數(shù)有泄漏口大小、流量系數(shù)、液池蒸發(fā)速率和儲罐內(nèi)壓力，其不確定性表現(xiàn)為：導(dǎo)致LNG船泄漏的事件不同(恐怖襲擊、碰撞或擱淺)，泄漏口大小可能在0～16 m2不等；流量系數(shù)受泄漏口形狀、LNG泄漏流動狀態(tài)和評估機構(gòu)主觀因素等影響，取值范圍較廣(0.4～1.0)，計算結(jié)果可相差2.5倍；LNG在水面大部分時間內(nèi)都處于膜沸騰狀態(tài)，當(dāng)液池蒸發(fā)快結(jié)束時轉(zhuǎn)為過渡沸騰狀態(tài)，蒸發(fā)速率在液池存在時間內(nèi)是變化的，LNG泄放入水的方式不同，與水的攪合激烈程度不同，蒸發(fā)速率會有很大不同(常取值0.042～0.25 kg/(m2·s))；LNG船貨艙是相對封閉的空間，當(dāng)發(fā)生事故性泄漏時(尤其是泄漏口徑較大時)，儲罐內(nèi)會產(chǎn)生一定的真空度，真空度的大小及其對泄漏速率的影響是很難確定的。此外，水面波浪、風(fēng)向及風(fēng)速等其他許多因素對LNG液池擴展的影響目前也很難確定。對LNG液池擴展變化的研究有待進(jìn)一步深入進(jìn)行。

[1] 趙敏,黃東風(fēng),佘孝云.從BP世界能源統(tǒng)計年鑒看中國能源發(fā)展[J].能源與環(huán)境,2014(6):16-17.

[2] 莊學(xué)強,鄭超瑜.液化天然氣船水上泄漏事故源強計算方法[J].大連海事大學(xué)學(xué)報,2010,36(4):125-128.

[3] HISSONG D W.Keys to Modeling LNG Spills on Water [J].Journal of Hazardous Materials,2007,140(3):465-477.

[4] KYTOMAA H,GAVELLI F.Studies of LNG Spills Over Water Point up Need for Improvement [J].Oil & Gas Journal,2005,103(18):61-65.

[5] HIGHTOWER M M,LUKETA-HANLIN A J,ATTAWAY S W.Breach and Safety Analysis of Spills Over Water from Large Liquefied Natural Gas Carriers [C]∥Listmail.Energy.Ca.Gov,2008.

[6] PITBLADO R M,BAIK J,HUGHES G J,et al.Consequences of LNG Marine Incidents[C].CCPS Conference Orlando,2004:1-19.

[7] BAIK J,RAGHUNATHAN V,WITLOX H.Consequence Modeling of LNG Marine Incidents[C].7th Professional Development Conference & Exhibition,2006:1-7.

[8] 王琳,季潔,賀寶龍,等.LNG地面泄漏蒸發(fā)速率的計算[J].天然氣與石油,2014,32(5):29-32.

[9] SOTERIOU C,ANDREWS R,SMITH M.Direct Injection Diesel Sprays and the Effect of Cavitation and Hydraulic Flip on Atomization[J].Simulation & Modeling,1995,120(16):84-93.

[10] WOODWARD J L,PITBLADO R M.LNG Risk Based Safety-Model & Consequence Analysis [M].John Wiley & Sons,Inc.,Hoboken,New Jersey,2010:134-145.

[11] DAHLSVEEN J,KRISTOFFERSON R,SAETRAN L R.Jet Mixing of Cryogen and Water[C].Second Symposium on Turbulence and Shear Flow Phenomena,Stockholm,Sweden,2001:329-334.

[12] 于桂峰,張彬,吳宛青.LNG船舶泄漏水面液池擴展的數(shù)值分析研究[J].大連海事大學(xué)學(xué)報,2008,34(S2):170-172.

[13] KHELIFI-TOUHAMI M S,BENBRIK A,LEMONNIER D,et al.Laminar Natural Convection Flow in a Cylindrical Cavity Application to the Storage of LNG[J].Journal of Petroleum Science and Engineering,2010,71(3-4):126-132.

PoolSpreadCalculationofLNGCarrierSpillsandUncertaintyAnalysis

ZHUANGXueqiang1,2，SUNDi1,2，GAOXiaohong3，LIGesheng3

(1.Marine Engineering Institute,Jimei University,Xiamen 361021,China; 2.Fujian Provincial Key Laboratory of Naval Architecture and Ocean Engineering,Xiamen 361021,China；3.School of Energy and Power Engineering,Wuhan University of Technology,Wuhan 430070,China)

A model for pool spread of LNG carrier spills is proposed based on the theory of orifice flow and inertia-gravity balance.The uncertainties of important parameters are analyzed,and the degree of the impacts of them on predicting results is calculated by means of the software MATLAB.The research reveals the following facts:Different kind of accidents (collision，grounding or terrorist attack)may result different LNG carrier leakage and the possible breach size can vary from 0 to 16 m2; The pool spread is sensitive to the breach size when it is less than 2 m2; The discharge coefficient,which is dependent on the breach shape and the liquid property,varies significantly,and the prediction deviation can be great,possibly up to 250%; The average evaporation rates also show wide variations (0.042～0.25 kg/(m2·s)) depending on the effects of the LNG-water turbulence and the boiling region.

waterway transportation; liquid pool; LNG carrier; release

U698.6;U674.13+3.3

A

2017-01-15

福建省科技重大平臺資助項目(2014H2001)；福建省科技(重點)項目(2017Y0065)；廈門市科技計劃項目(3502Z20173030)；廈門南方海洋研究中心項目(14GZB66NF30)；集美大學(xué)博士啟動基金(ZQ2013006)

莊學(xué)強(1974—),男,福建泉州人，副教授，碩士生導(dǎo)師，博士，從事水上危險品運輸安全與防污染研究。E-mail:13015917291@163.com

孫 迪(1979—)，女，遼寧沈陽人，講師，博士，主要研究方向為柴油機關(guān)鍵摩擦副摩擦監(jiān)測技術(shù)。E-mail:sundi6329@sina.com

1000-4653(2017)02-0073-05