張小斌，厲勁風，邱利民（浙江大學制冷與低溫研究所，浙江省制冷與低溫技術(shù)重點實驗室，浙江 杭州 310027）

?

液化天然氣排放形成的羽流過程數(shù)值研究

張小斌，厲勁風，邱利民
（浙江大學制冷與低溫研究所，浙江省制冷與低溫技術(shù)重點實驗室，浙江 杭州 310027）

摘要：大量LNG溢出蒸發(fā)后將形成低溫云團，對風向下游地面人員可能造成凍傷、燃燒以及缺氧窒息等危險。基于計算流體力學方法，構(gòu)建了深低溫兩相多組分流動的Navier-Stokes方程以及湍流封閉方程，考慮了空氣中水蒸氣由于溫度降低到飽和溫度以下而相變傳質(zhì)過程。由于氣相中氧氣等非液化氣體的存在，修正了計算水蒸氣相變率的Hertz-Knudsen方程。詳細給出了計算過程邊界條件設置，評估了地球自轉(zhuǎn)引起的科里奧利力影響特性。數(shù)值模擬了美國LLNL完成的Burro系列LNG排放羽流實驗，發(fā)現(xiàn)考慮水蒸氣相變的計算結(jié)果，相比于未考慮水蒸氣相變的計算，更接近實驗結(jié)果。研究結(jié)果對LNG接收終端等安全環(huán)境評估及設計有指導意義。

關(guān)鍵詞：羽流；天然氣；低溫；凝結(jié)；數(shù)值分析

2015-05-29收到初稿，2015-10-16收到修改稿。

聯(lián)系人：邱利民。第一作者：張小斌(1976—)，男，教授。

Received date: 2015-05-29.

Foundation item: supported by the Outstanding Youth Foundation of Zhejiang Province (R15E060001) and the National Natural Science Foundation of China (51576169).

引　言

液化天然氣（LNG）在儲存和運輸中由于泄漏事故引起的安全問題受到廣泛關(guān)注[1-2]。LNG接收終端儲罐的儲量規(guī)格普遍為10萬多立方米，一旦破裂大量LNG（主要為甲烷）將頃刻溢出，于地面或海面快速相變蒸發(fā)后形成低溫云。低溫云團大概以相等或稍小于空氣速度向下游傳播，云團周邊氣體逐漸與空氣混合稀釋并升溫，同時將空氣中的水蒸氣液化成飽和液滴，形成羽流。云團核心區(qū)在較長時間內(nèi)保持低溫以及大于空氣的密度，因此云團向下游運動時高度可能下降到地面。這將對地面人員造成凍傷、燃燒等危險，也可能導致缺氧災難。據(jù)報道，當氧氣含量小于17.7%時，人將出現(xiàn)呼吸困難、頭昏乏力等現(xiàn)象[3]。云團到達地面時離氣源的距離、組分含量及覆蓋面積等取決于環(huán)境風速、濕度，氣源高度、溫度及流量，以及地形結(jié)構(gòu)等。因此，獲知不同條件下LNG儲罐溢流引起的羽流可能對地面人員造成危險的安全距離等參數(shù)是建設LNG儲罐必須的安全條件。

Roberts等[4-5]首先得到基于簡化假設的大氣擴散問題的理論解，該理論用以分析相關(guān)參數(shù)及源項的敏感性分析。高斯(Gaussian)羽流[6-7]被認為是最簡單的羽流情況，指點源連續(xù)無方向的向無限空間環(huán)境排放擴散組分。高斯羽流理論解已被廣泛應用于大工業(yè)氣體排放，工業(yè)及農(nóng)業(yè)養(yǎng)殖場氣味傳輸，火山爆發(fā)以及核、生化工業(yè)等的污染物排放預測，并最終形成污染物排放安全工業(yè)標準。

計算流體力學(CFD)方法也已被廣泛應用于LNG排放形成羽流擴散過程濃度的時域及地域預測。羽流擴散CFD計算一般基于如Fluent，CFX及Star-CD等商業(yè)軟件[8]。另外也有專門用于羽流排放過程而開發(fā)的數(shù)值計算軟件如SLAB、DEGADIS 和FEM3等，但這些程序目前只能計算穩(wěn)定邊界條件，不能計算組分濃度非穩(wěn)態(tài)波動[9]。其中DEGADIS假設組分濃度在垂直方向滿足指數(shù)定律分布，水平方向滿足修正的高斯分布，同時假設垂直方向風速滿足指數(shù)分布，則計算簡化為一維模型，只能用于分析重氣體排放到平坦、無障礙地形的情況。SLAB也是基于簡化的一維模型，DEGADIS也假設濃度的橫向和高度方向的分布函數(shù)，不能模擬建筑物等對組分傳播的影響[10-11]。在羽流數(shù)值分析中廣泛使用的商業(yè)軟件為FEM3，這是個基于有限元方法的模型包，由美國Lawrence Livermore National Laboratory (LLNL)開發(fā)，最初湍流模型采用混合長理論（一方程模型），現(xiàn)在已經(jīng)包含了κ-ε二方程模型。另外日本三菱重工也為羽流擴散的數(shù)值分析開發(fā)了CFD有限元軟件STD。

20世紀80年代LLNL進行了多次不同LNG排放條件的羽流實驗，包括Burro系列及Coyote系列[12]，這些實驗得到了不同位置非穩(wěn)態(tài)變化的濃度和溫度，成為驗證排放過程CFD計算結(jié)果準確性的基礎數(shù)據(jù)。如Sklavounos等[10-11]基于CFX軟件計算，通過實驗對比評估了基于標準κ-ε模型，標準κ-ω和SST湍流模型的計算結(jié)果，發(fā)現(xiàn)這些湍流模型都能產(chǎn)生相同精度的結(jié)果，但κ-ω低估了最大濃度，而其他模型高估。Gavelli等[13]基于Fluent計算了有建筑物時的Falcon系列實驗過程，對比了濃度和溫度分布，指出DEGADIS計算結(jié)果不能滿足工程設計要求。而Tauseef等[14]利用Fluent模擬了建筑物對組分離散的影響，對比了不同湍流模型結(jié)果，發(fā)現(xiàn)Real κ-ε模型產(chǎn)生了最合適和精確的濃度分布。

上述對于大流量LNG排放到大氣環(huán)境形成復雜兩相多組分擴散羽流過程的CFD研究，都沒有考慮空氣和低溫氣體強烈湍流混合后當?shù)販囟燃眲〗档?，導致水蒸氣飽和形成液滴的相變物理過程。然而LLNL的實驗明顯觀察到形成水霧特性，因此忽略水蒸氣相變過程的羽流CFD計算的準確性，需要驗證。

本文構(gòu)建了考慮水蒸氣相變的基于混合物模型(mixture model)的兩相多組分數(shù)值模型，并用Fluent14.5求解控制方程。數(shù)值模擬了Burro系列實驗中LNG排放過程。

1　兩相多組分流動過程的數(shù)值模型

1.1基本控制方程

對于混合物模型，將氣液兩相看成是混合物單相，其連續(xù)性方程為

假設形成的液滴與氣相運動速度相同，即兩者沒有動量交換，則動量方程為

式中，mm=alml+agmg為混合物黏度；P為壓力；g為重力加速度矢量。

能量方程為

有效熱導率keff=alkl+agkg，hk為k相焓值。

氣相體積含量ag方程為

考慮水蒸氣相變傳質(zhì)，組分i在氣相混合物中的質(zhì)量含量Ygi，由式(5)計算

上角標或下角標g、l分別為氣相和液相。對于湍流，式(5)中的擴散流量Ji由式(6)計算

Sct是湍流有效Schmidt數(shù)，μt/Sct是湍流引起的有效質(zhì)量擴散系數(shù)。因此，在湍流條件下的有效擴散包含作為已知條件輸入的層流擴散系數(shù)Di,m，加上增加的湍流擴散效應。湍流Sct表示動量和質(zhì)量的相對擴散，數(shù)量級大概為1，且對分子流體屬性不敏感，因此采用缺省的常數(shù)Sct=0.7是合理的。另外，湍流擴散效應遠遠大于層流擴散效應，采用常層流擴散系數(shù)，Di,m=2.8×10-5m·s-1也是可以接受的。

為封閉式(1)～式(5)，必須計算相變引起的水汽化和液化量

1.2相變傳質(zhì)模型

假設氣液界面溫度平衡，界面為幾個自由分子厚度形成的Knudsen層，水蒸氣為理想氣體分子，基于分子運動論可得到計算氣液相變率的Hertz-Knudsen關(guān)系[15]

式中，h為液化或蒸發(fā)系數(shù)。式(7)由純組分條件推導得到，當考慮到氣相中存在非液化的氧氣和氮氣時，Pv指的是水蒸氣的分壓力。其他非液化氣體對水蒸氣相變率的影響，主要通過修正通用氣體常數(shù)R來體現(xiàn)[16]

R=8314.34 J·kmol-1·K-1，Mi為組分i的分子量。多組分氣相中，水蒸氣分壓PH2O由式(9)計算

式中，yH2O為水蒸氣摩爾含量。在實際計算中，通常以溫度作為判斷是否發(fā)生相變的標準。利用Clausius-Clapeyron方程，式(8)可重寫為

式中，Tsat為飽和溫度，為水蒸氣分壓PH2O的函數(shù)；rg為飽和水蒸氣密度；L為潛熱。需注意的是傳質(zhì)率m單位為kg·m-2·s-1。在數(shù)值計算中，需要轉(zhuǎn)變?yōu)閱挝惑w積的傳質(zhì)率。假設生成的液滴直徑相同，則界面積密度為

Vcell為網(wǎng)格體積，d為液滴直徑。由式(10)和式(11)得到

T*為當?shù)貙嶋H溫度。引入系數(shù)C（L·s-1）

則式(12)最后簡化為

需要說明的是，液化與蒸發(fā)過程系數(shù)C絕對值不相等，具體值通過和實驗對比得到。由于Tsat為水蒸氣分壓力的函數(shù)，因此式(14)需通過自定義函數(shù)(UDF)的方法，耦合到Fluent的程序中進行迭代計算。

為簡化計算，并沒有考慮可能發(fā)生的水液固相變?？紤]到液態(tài)水在混合物中體積含量很小，因此雖然液固相變時密度從1000 kg·m-3減小到冰的900 kg·m-3，對流場的影響認為可忽略不計。同樣的原因，雖然比熱容從水的4200 J·kg-1·K-1減小到冰的2100 J·kg-1·K-1，對溫度場的影響也認為可忽略不計。

1.3湍流模型

κ-ε二方程湍流模型已經(jīng)被成功地用于大氣羽流的數(shù)值模擬[17-18]。文獻[8]指出，缺省的κ-ε二方程模型參數(shù)對計算結(jié)果的濃度分布誤差很小。Sklavounos等[10]利用軟件CFX建模了Coyote實驗，也發(fā)現(xiàn)使用缺省的湍流參數(shù)與實驗結(jié)果吻合很好。因此，本文計算采用缺省的Realizable κ-ε二方程湍流模型。

1.4科里奧利力影響

由于地球自轉(zhuǎn)，大氣流動坐落于有加速的參考坐標中，運動方程有必要考慮科里奧利力的影響。量綱1參數(shù)Rossby數(shù)(Ro)比較了慣性力與科里奧利力的大小[17]

式中，V、D和f分別為特征速度、特征長度和科里奧利參數(shù)f=2?sin(φ)，φ為緯度角。由于慣性力在空氣流動過程一般占主導作用，因此，Ro通常用來判斷是否需要考慮科里奧利力。如果Ro較小，則科里奧利力應該考慮。以溢出10萬立方米 LNG到水面為例，如果風速為2 m·s-1，f約10-4s-1，D取天然氣云最大值約1 km，那么Ro約20，說明科里奧利力相比慣性力不重要。本文忽略科里奧利力影響，目前也未見考慮科里奧利力影響的羽流CFD數(shù)值計算報道。

2　CFD數(shù)值模型邊界條件

為驗證上述數(shù)理方程組及相變模型的準確性，選擇Burro系列[19]實驗數(shù)據(jù)進行對比。實驗過程獲得了不同位置的溫度、風速和CH4組分非穩(wěn)態(tài)變化，表1給出了實驗中相關(guān)參數(shù)。

表1　Burro 9系列LNG排放實驗參數(shù)[19]Table 1　Test conditions of LNG dispersion in Burro 9[19]

為將上述物理過程轉(zhuǎn)變?yōu)镃FD數(shù)值模型，進行了與文獻[10]同樣的簡化，包括：

① 由于Burro 9中主要組分皆為甲烷氣體(CH4)，且所占比重較大（Burro 9為87.4%，Burro 9 為83.1%），忽略少量乙烷及丙烷不會對流體物性產(chǎn)生較大影響。因此不考慮LNG組分組成，蒸發(fā)的天然氣近似為100% 甲烷氣體；

② 本文模擬LNG在水池中泄漏并快速相變之后以氣態(tài)NG在環(huán)境中擴散的過程，因此忽略LNG在水池中泄漏并快速相變的過程，整個水池液面(d=58 m)為氣體CH4進口邊界，溫度T=111.5 K，速度根據(jù)表1給定的LNG排放量計算，水池面積則CH4流速為

③ 忽略現(xiàn)實中存在的風向隨機變化，羽流關(guān)于水池中心垂直平面對稱。

基于上述假設，構(gòu)建了如圖1所示的平面對稱CFD數(shù)值模型，該模型和文獻[10]幾乎一致，除了后者沒有利用平面對稱條件。

圖1　Burro系列實驗的CFD建模邊界條件Fig.1　CFD geometrical model and boundary conditions of Burro serial LNG spill tests

上述計算域有7個邊界，邊界條件設置分別如下：

① 右邊出口為定壓邊界條件，設置相對壓力P=0，平均壓力為大氣壓；

② 計算域頂面高度足夠大，使得底部CH4的進入對頂部流場影響可忽略不計。因此，理論上，可設置為定壓邊界，但頂部定壓邊界將改變左邊空氣入口的速度分布。由此，在忽略風向隨機變化的簡化條件下，頂面設置為對稱邊界條件，滑移速度為0，邊界速度與流場速度相等，也沒有流體穿過對稱面；

③ 由假設③，圖中計算域正面為對稱平面；

④ 計算域的背面，保證計算域足夠?qū)?250 m)，則CH4進入對該面的影響可忽略不計，因此相同于計算域頂面，背面也設置為對稱面；

⑤ 計算域底面（模擬地面），由于實際LNG排放時間較短（約100 s），與地面換熱只影響到很薄的邊界層，對CH4宏觀擴散影響可忽略不計，因此設置為絕熱面(Q=0)；

⑥ CH4進口為速度進口邊界條件，其值為假設②計算。入口湍流強度及黏度比都為4%。入口I值實際是未知值，因此有必要進行I對計算結(jié)果影響的敏感性分析。

⑦ 計算域左邊空氣進口條件最為復雜，因為空氣速度、溫度和濕度都是高度的函數(shù)。M-O相似律[12]給出了精確計算方法，但該方法較復雜。另外一種更簡單的計算高度方向速度分布的方法是指數(shù)定律

式中，zr為參考高度，Ur為在高度zr的水平風速，z為高度坐標，指數(shù)p為按照Pasquill- Gifford分類方法由大氣穩(wěn)定性條件決定的系數(shù)[20]。該方法雖然沒有M-O相似律精確，但在低高度范圍內(nèi)兩者吻合很好，因此，已廣泛應用于該類問題的CFD0.1計5算[13,21]。由表1可知，對在計算中，利用UDF方法將上式公式指定為邊界上的速度分布。

入口空氣溫度分布，假設為表1中所給的常數(shù)，文獻[8,10]也作同樣處理。

入口空氣濕度分布，也假設為常數(shù)，根據(jù)表1中相對濕度r計算得到體積含量

式(1)～式(5)以及式(14)加上Realizable κ-ε湍流模型，構(gòu)成了封閉的基于混合物模型的兩相多組分數(shù)值框架，利用CFD商業(yè)軟件Fluent14.5進行迭代求解。其中連續(xù)性方程和動量方程通過二階迎風方案離散，利用Couple方案計算速度與壓力修正。計算收斂標準為能量方程殘差為10-6，其他方程殘差為10-3。計算過程，純組分氧氣、氮氣、CH4及水蒸氣的熱導率、比熱容及黏度都為溫度的函數(shù)，液態(tài)水的物性為常數(shù)，物性參數(shù)來自Refprop 8.0。整個計算域超過17.73萬個網(wǎng)格，將CH4進口附近局部加密后網(wǎng)格超過19.12萬個，參照Burro 9實驗，以相對排放源中心為原點，坐標x=397 m、y=48.75 m處的T3傳感器為比較對象，網(wǎng)格獨立性檢驗計算表明基于17.73萬個網(wǎng)格數(shù)的CH4濃度值幾乎和后者一致，見圖2。

圖2　不同網(wǎng)格數(shù)量下計算的Burro 9 T3傳感器的CH4濃度非穩(wěn)態(tài)變化Fig.2　Simulated CH4concentration of Burro 9 T3 sensor with different mesh number

3　數(shù)值結(jié)果分析

首先對式（14）中系數(shù)C對相變率的定量影響進行評估。實驗結(jié)果顯示，液化系數(shù)Cc一般要大于蒸發(fā)系數(shù)Ce[15]，Rubel等[22]的測量更進一步指出Ce≈1.2Cc。采用兩組系數(shù)對，Ce=1，Cc=1.2及Ce=100，Cc=120，兩者相差100倍，對Burro 9實驗進行了數(shù)值計算，結(jié)果如圖3所示。由圖可見，兩者計算結(jié)果吻合較一致。因此，計算中采用液化和蒸發(fā)系數(shù)對：Ce=100和Cc=120。

圖3　不同系數(shù)C時計算的Burro 9在不同位置的CH4濃度非穩(wěn)態(tài)變化Fig.3　Simulated CH4concentration of Burro 9 with different evaporation and condensation coefficient pairs

圖4和圖5分別給出了考慮相變及未考慮相變時，計算的Burro 9溫度和CH4濃度分布與實驗測量的對比[19，23]。根據(jù)實驗布置，圖中傳感器位置：以圖1中CH4進口中心為坐標原點: G6: x=137.9 m，y=24.3 m；G15: x=385.25 m，y=107.55 m；T3: x=397 m，y=48.75 m；T4: x=132.32 m，y=46.47 m。這些坐標值根據(jù)實際傳感器位置以及風向平均偏轉(zhuǎn)角修正計算得到。由圖可見，無論溫度還是濃度，實驗結(jié)果呈無規(guī)律波動，這是由實驗過程上游風速和風向隨機變化導致。而數(shù)值計算固定風向和風速，因此結(jié)果曲線光滑?？傮w上，從實驗數(shù)據(jù)依然能捕捉到濃度和溫度的宏觀變化，因此可用于檢驗數(shù)值結(jié)果。在圖4中，CFD計算的不同位置的溫度，無論是考慮相變或不考慮相變，和實驗值趨勢一致，但最低溫度都明顯低于實驗值，這和數(shù)值模型與實驗邊界條件不完全一致有關(guān)。首先，實際風向及風速隨機高頻率變化導致實驗時云團中心平面與計算模型的云團中心平面不一致，因此實驗中傳感器位置與對應的計算模型中的坐標，相對各自云團中心平面位置卻不一致（因為實驗時云團中心平面位置隨風向在變化）。顯然云團中心平面上溫度最低，導致測量溫度比計算值高。另外，傳感器離地面1米，云團和地面的換熱以及地面的輻射對會導致云團溫度上升?？傮w上，CFD計算結(jié)果能夠得到可接受的溫度非穩(wěn)態(tài)變化。另外從圖中可發(fā)現(xiàn)，考慮相變的CFD模型，得到的溫度總是大于未考慮相變的CFD模

圖4　計算的Burro 9 溫度分布與實驗測量對比(h=3 m)Fig.4　Comparison of simulated temporal temperature curves with experimental records for Burro 9 (h=3 m)

圖5　計算的Burro 9 CH4濃度分布與實驗對比(h=1 m)Fig.5　Comparison of simulated temporal CH4concentration curves with experimental records for Burro 9 (h=1 m)

對CH4濃度變化，由圖5可見，未考慮相變的CFD模型，幾乎在受云團影響的時間內(nèi)，在所有位置比考慮相變模型的CFD計算結(jié)果高。這是由于后者溫度普遍比前者高（圖4），使得氣體黏度降低，Re增加，云團中湍流強度更加強烈，由此導致氣體分子運動更活躍，CH4組分擴散及流動傳質(zhì)加劇，因此CH4更容易被空氣稀釋，濃度相對降低?？傮w上，在所有位置計算的CH4濃度值和實驗測量的平均值趨勢和大小基本一致，但G15位置計算值始終要小于實驗值。由于G15距離CH4源最遠，風向偏移造成的影響最大。

以Burro 9的G6和T3傳感器所得到的CH4濃度分布實驗數(shù)據(jù)為參考，圖6給出了Fluent模擬結(jié)果與FEM3軟件模擬結(jié)果[8]的對比。由圖6可見，基于Fluent的模擬結(jié)果要比基于FEM3的模擬結(jié)果更接近實驗數(shù)據(jù)。

圖6　計算的Burro 9 CH4濃度分布與實驗及FEM3軟件結(jié)果對比 (h=1 m)Fig.6　Comparison of Fluent simulated temporal CH4concentration curves with experimental records and simulation results of FEM3 for Burro 9 (h=1 m)

4　結(jié)　論

本文構(gòu)建了兩相多組分流動的Navier-stokes方程以及湍流封閉方程，給出了基于氣體動力學的羽流過程水蒸氣相變計算的修正模型。計算了Burro系列實驗中LNG排放羽流實驗，得到的溫度及濃度非穩(wěn)態(tài)變化與實驗進行了對比，得出如下結(jié)論。

（1）大流量低溫流體如LNG等排放形成羽流過程的CFD數(shù)值模型，空氣中水蒸氣飽和形成雨滴的物理過程不可忽略。考慮水蒸氣相變過程的計算結(jié)果，相比于未考慮水蒸氣相變過程的計算，與實驗結(jié)果相比更一致。

（2）基于Fluent的數(shù)值計算，相對于FEM3等傳統(tǒng)的擴散模擬軟件的模擬計算，與實驗結(jié)果擬合度更高。

（3）構(gòu)建的基于mixture model兩相多組分數(shù)值模型及基于氣體動力學的水蒸氣相變計算模型，能準確模擬羽流過程的水蒸氣相變過程。若能在模型中增加來流空氣隨機性波動的考慮，將進一步提高數(shù)值計算結(jié)果的準確性。

References

[1] 錢新明，劉牧，劉振翼. 隧道內(nèi)液化天然氣管道泄漏火災溫度場的數(shù)值模擬 [J]. 化工學報，2009，60(12): 3184-3188. QIAN X M，LIU M，LIU Z Y. Numerical simulation of LNG leakage fire temperature distribution from pipeline in tunnel [J]. CIESC Journal，2009，60(12): 3184-3188.

[2] 唐建峰，蔡娜，郭清，等. LNG垂直噴射源連續(xù)泄漏擴散的模擬[J]. 化工學報，2013，64(3): 1124-1131. DOI: 10.3969/j.issn. 0438-1157.2013.03.048. TANG K F，CAI N，GUO Q，et al. Simulation of LNG diffusion: a continuous vertical jet release [J]. CIESC Journal，2013，64(3): 1124-1131. DOI: 10.3969/j.issn.0438-1157.2013.03.048.

[3] FLYNN T M. Cryogenic Engineering[M]. New York: Marcel Dekker Inc.，2005.

[4] ROBERTS O F T. The theoretical scattering of smoke in a turbulent atmosphere[C]//Proceedings of The Royal Society of London Series A-Containing Papers of A Mathematical and Physical Character. London，England: Royal Soc London，1923，104: 640-654.

[5] SUTTON O G. A theory of eddy diffusion in the atmosphere [C]//Proceedings of The Royal Society of London Series A-Containing Papers of A Mathematical And Physical Character. London，England: Royal Soc London，1932，135: 143-165.

[6] STOCKIE J M. The mathematics of atmospheric dispersion modeling [J]. SIAM Review，2011，53(2): 349-372.

[7] SEINFELD J H. Atmospheric Chemistry and Physics: from Air Pollution to Climate Change[M]// PANDIS S N. New York: JohnWiley & Sons Inc.，1998.

[8] LUKETA-HANLIN A，KOOPMAN R，ERMAK D L. On the application of computational fluid dynamics codes for liquefied natural gas dispersion [J]. Journal of Hazardous Materials，2007，140: 504-517.

[9] OHBA R，KOUCHI A，HARA T，et al. Validation of heavy and light gas dispersion models for the safety analysis of LNG tank [J]. Journal of Loss Prevention in the Process Industries，2004，17: 325-337.

[10] SKLAVOUNOS S，RIGAS F. Simulation of Coyote series trials(Ⅰ): CFD estimation of non-isothermal LNG releases and comparison with box-model predictions [J]. Chem. Eng. Sci.，2006，61: 1434- 1443.

[11] SKLAVOUNOS S，RIGAS F. Validation of turbulence models in heavy gas dispersion over obstacles [J]. Journal of Hazardous Materials，2004，A108: 9-20.

[12] ERMAK D L，OHAPMAN R C V，GOLDWIRE H C，et al. Heavy gas dispersion test summary report[R]. 1989，ESL-TR-88-22.

[13] GAVELLI F，BULLISTER E，KYTOMAA H. Application of CFD (Fluent) to LNG spills into geometrically complex environments [J]. Journal of Hazardous Materials，2008，159: 158-168.

[14] TAUSEEF S M，RASHTCHIIAN D，ABBASI S A. CFD-based simulation of dense gas dispersion in presence of obstacles [J]. Journal of Loss Prevention in the Process Industries，2011，24: 371-376.

[15] MAREK R，STRAUB J. Analysis of the evaporation coefficient and the condensation coefficient of water [J]. International Journal of Heat and Mass Transfer，2001，44: 39-53.

[16] GRAYSON G，LOPEZ A，CHANDLER F. CFD modeling of helium pressurant effects on cryogenic tank pressure rise rates in normal gravity[C]// AIAA. Cincinnati，Ohio: 2007，MSFC-448.

[17] ALINOT C，MASSON C. k-ε model for the atmospheric boundary layer under various thermal stratifications [J]. Solar Eng.，2005，127: 438-443.

[18] HUSER A，NILSEN P J，SKATUN H. Applications of k-ε model to the stable ABL: pollution in complex terrain [J]. Wind Eng. Ind. Aero，1997，67/68: 425-436.

[19] KOOPMAN R P，BAKER J，CEDERWALL R T. Burro series data report LLNL/NWC 1980 LNG spill tests[R]. 1982，UCID-19075.

[20] CCPS. Consequence analysis of chemical releases[C]// AIChE. New York: Wiley-Blackwell，1999: 80-83

[21] MCBRIDE M A，REEVES A B，VANDERHEYDEN M D，et al. Use of advanced techniques to model the dispersion of chlorine in complex terrain [J]. Process Safety & Environmental Protection，2001，79(2): 89-102.

[22] RUBEL G O，GENTRY J W. Measurement of the kinetics of solution droplets in the presence of adsorbed mono-layers: determination of water accommodation coefficients [J]. Phys. Chem.，1984，88(14): 3142-3148.

[23] CHAN S T，RODEAN H C，ERMAK D L. Numerical simulations of atmospheric releases of heavy gases over variable terrain[C]//Air Pollution Modeling and Its Applications Ⅲ. New York: Plenum Press，1984，5: 295-328.

Numerical study on plume characteristics of liquefied natural gas spills

ZHANG Xiaobin，LI Jingfeng，QIU Limin
（Key Laboratory of Refrigeration and Cryogenic Technology of Zhejiang Province，Institute of Refrigeration and Cryogenics，Zhejiang University，Hangzhou 310027，Zhejiang，China）

Abstract:When a large number of LNG suddenly leaks and evaporates，it will form a low temperature cloud in the wind downstream，which may cause frostbite，burn and oxygen deficit hazard to the ground personnel. Based on the computational fluid dynamics (CFD) method，the two-phase multicomponent flow mathematical framework together with the turbulence closure are built to model the cryogenic flow. The phase-change mass transfer of water vapor in the air due to the temperature depression is considered. Because of the existence of non-liquefied gases such as oxygen，the Hertz-Knudsen equation for calculating the mass transfer rate of water vapor is modified. Detailed methods for setting the boundary conditions of the computational domain are presented and the influence of Coriolis force caused by the earth rotation is evaluated. The Burro experimental series of LNG released by LLNL are simulated and the results are used to evaluate the numerical models. It is found that the results using models with the phase change of water vapor are closer to the experimental results than that without the phase change. The studies are of directive significance for the safety environment assessment and design of LNG received terminal.

Key words:plume flow; NG; cryogenic; condensation; numerical analysis

DOI：10.11949/j.issn.0438-1157. 20150740

中圖分類號：TQ 021.4

文獻標志碼：A

文章編號：0438—1157（2016）04—1225—08

基金項目：浙江省杰出青年基金項目（R15E060001）；國家自然科學基金項目（51576169）。

Corresponding author:Prof. QIU Limin，liminqiu@zju.edu.cn