梅 苑，帥 健，任 飛

(中國石油大學(北京)安全與海洋工程學院，北京 102249)

隨著世界能源結構的調整，城市天然氣管網迅速發(fā)展。天然氣管網的建成給人們的生活及生產帶來了巨大的便捷，但也帶來了相應的安全隱患[1-2]。如：2010年9月9日，加利福尼亞州圣布魯諾天然氣管道發(fā)生破裂事故，事故后果造成維護費用高達1 350萬美元；2018年6月10日，中緬天然氣輸氣管道發(fā)生泄漏燃爆事故，造成1人死亡、23人受傷，直接經濟損失達2 145萬元。因此，從安全角度來看，對于天然氣管道泄漏擴散規(guī)律的研究具有重要意義。

現場試驗、風洞試驗和數值模擬方法是研究氣體泄漏擴散的重要手段，多年來學者們基于以上方法對可燃氣體泄漏擴散特性展開了大量的研究工作。如：Krogstad等[3]通過風洞試驗研究了重氣連續(xù)泄漏后矩形建筑物對重氣云羽擴散的影響，結果發(fā)現重氣云羽擴散時出現分叉現象，且建筑物底部產生的漩渦將會減小其周圍的重氣濃度，對重氣云羽的擴散造成了極大的影響；Konig-Langlo等[4]通過風洞試驗對比了瞬時泄漏源與連續(xù)泄漏源在穩(wěn)定大氣流動和湍流條件下對重氣擴散的影響，證實了重氣密度比對其擴散起次要作用，并得到不利大氣條件下重氣可燃距離的計算方法；Cowan[5]模擬了簡單建筑物側面發(fā)生氣體泄漏的情況，并將數值模擬結果與風洞試驗結果進行了對比分析；張甫仁等[6]利用計算流體力學(CFD)軟件對建筑物群外空間內城市燃氣連續(xù)泄漏擴散過程進行了數值模擬，對比分析了環(huán)境溫度、濕度對燃氣泄漏擴散的影響和濃度場的變化規(guī)律；高煒等[7]借助高斯擴散模型研究了不同大氣環(huán)境下天然氣管道泄漏擴散的危險范圍，結果發(fā)現提高風速會顯著降低對天然氣管道泄漏事故后果的影響；周寧等[8]采用fluent平臺針對不同風速下丁烷泄漏擴散過程進行了數值模擬研究，結果發(fā)現風速增大有利于減小丁烷高濃度區(qū)域的面積，從而減小事故危害。

目前，學者們針對氣體泄漏擴散的試驗研究大多采用縮比性試驗平臺，較難實現大空間場地內天然氣泄漏擴散過程的研究，而在使用模擬方法時忽略了CFD軟件在流場預測方面的潛力，僅從氣云的表象變化展開分析。鑒于此，本文基于CFD軟件對不同風速下天然氣管道泄漏擴散過程進行了數值模擬研究，重點研究了風速作用下建筑區(qū)內流場變化特性對于天然氣氣云擴散的影響。

1 數值建模和工況設置

1.1 控制方程

Fluent軟件在氣體泄漏擴散模擬仿真過程中設置了更多的監(jiān)測點，可使測量數據更為全面、有效[9-11]，部分學者[12-14]通過試驗驗證了其在天然氣管道泄漏擴散預測方面的可靠性，因此本文采用Fluent軟件作為研究手段。利用Fluent軟件進行模擬計算時需要將擴散氣體的質量、能量等守恒方程作為控制方程，并利用Realizablek-ε湍流模型和組分運輸模型來模擬氣體泄漏擴散過程。在流體力學中流動運動的連續(xù)性方程、動量方程和能量守恒方程可通過統(tǒng)一的方程表示如下[15]：

(1)

式中:φ代表某一變量；Γ表示流體的擴散系數；S表示源項；Γ和S均對應特定的變量φ。

公式(1)從左到右的4項分別為時間項、對流項、擴散項和源項，取不同的φ、Γ和S，則可得到相應的連續(xù)性方程、動量方程和能量方程。

1.1.1 湍流動能方程

Realizablek-ε湍流模型比標準k-ε湍流模型在濃度分布上有更好的精度。對氣體泄漏擴散湍流問題采用Realizablek-ε湍流模型，模型中k和ε方程如下：

k方程為

(2)

ε方程為

(3)

式中：μt為湍動黏度(kg·m·s)；E為時均應變率(%)；ρ為流體密度(kg/m3)；ui為時均速度(m/s)；xi為控制單元體長度(m)；k為湍流動能(J)；ε為湍流動能耗散率(%)；C1、C1ε、C3ε為3個常量；σk=1.0，σε=1.2，C2=1.9；Gk為因速度梯度產生的湍流動能源項；Gb為因浮力產生的湍流動能源項；YM為在可壓縮湍流中波動擴張引起的耗散項；θ為運動黏度(m2/s)。

1.1.2 組分運輸方程

組分運輸方程如下：

(4)

式中：u、v、w為流體速度分別在x、y、z軸上的速度分量(m/s)；mi為不同組分所占的質量比例；Γi為湍流擴散系數。

1.2 模擬工況設置

本文模擬了5種風速條件下天然氣管道泄漏擴散過程，模擬工況設置的詳細信息見表1。

表1 模擬工況設置信息

1.3 數值模型建立與參數設置

從對氣體泄漏擴散影響的角度來看，建筑物的布局方式可分為封閉式、斑塊式和街道峽谷式，而封閉式的建筑物排列方式不僅在城區(qū)中最為常見，也是對氣體泄漏擴散影響最大的建筑物布局方式[16]。因此，本文在計算域中加入封閉式建筑物排列。為了保證計算域內流場充分發(fā)展，減小邊界層對流場發(fā)展的影響，計算域的尺寸需盡可能大，但這也相應地增加了計算成本。根據阻塞率原則(即建筑群的阻塞率應≤3%)，本文確定計算域的空間尺寸為530 m×380 m×150 m(x，y，z)，建筑物的尺寸為15 m×30 m×30 m(x，y，z)，此時阻塞率為2.3%，滿足阻塞率原則，建立的封閉式建筑群內天然氣管道泄漏擴散的物理幾何模型， 見圖1。天然氣與環(huán)境溫度均設置為296.3 K，環(huán)境壓力設置為標準大氣壓。天然氣管道泄漏口設置于建筑群空腔區(qū)內，為直徑為105 mm的圓形泄漏孔，其邊界條件設置為質量入口邊界。泄漏的天然氣管道管徑為200 mm，運行壓力為0.35 MPa，屬于中壓輸送水平。利用大孔泄漏計算模型[17]得出泄漏孔徑為105 mm條件下天然氣管道泄漏速率約為5.0 kg/s。風入口面及出口面的邊界條件分別設置為速度入口、壓力出口；地面及建筑的邊界條件采用壁面邊界條件；其余邊界條件皆采用對稱邊界條件。采用PISO算法對壓力場和速度場進行耦合，對流項離散采用二階逆風格式，擴散項采用二階中心差分格式。

圖1 封閉式建筑群內天然氣管道泄漏擴散的物理幾何模型

1.4 網格劃分與模型驗證

本文借助ICEM軟件對計算域進行結構化網格劃分，從而提升整體網格質量。針對建筑區(qū)及泄漏源附近的網格進行局部加密，但網格數量過大，會導致計算成本增加。因此，本文共劃分4組數量梯度的網格來分析驗證模型的獨立性，網格數量依次為145萬、198萬、241萬和296萬，并將各數量梯度網格模擬計算結果與試驗數據[16]進行了對比(見圖2)，結果發(fā)現：隨著天然氣泄漏時間的增加，因網格數量產生的計算誤差會逐步拉大，且當網格數量達到241萬和296萬時，網格數量對天然氣氣云擴散范圍和濃度分布的影響較小，可以達到較高的計算精度。因此，考慮到計算成本，最終確定計算網格數量為241萬。

圖2 不同風格數量模擬結果與試驗結果的對比

2 數值模擬結果與分析

2.1 封閉式建筑群風場模擬

風場的發(fā)展會直接影響天然氣管道泄漏事故后果。為了獲得可靠的預測結果，本文采用兩步模擬的方法，即首先進行風場模擬，預算300 s，以便獲得較好的風場初始值；然后開啟泄漏源，進行天然氣管道泄漏事故后果模擬。大氣流動是一種復雜的湍流運動，大氣邊界層的流動特性通常采用平均風廓線和湍流參數來描述。根據Davenport的實測數據[18]，采用冪律方程來描述入流邊界處的風速分布，即：

(5)

式中:z0為參考高度(m)，取10 m；z為距離地面的高度(m);α為地面粗糙度指數；v0為距離地面10 m高度處的平均風速(m/s)，文中提及的風速皆指距離地面10 m高度處的平均風速；vz為距離地面高度z處的平均風速(m/s)。

如圖3所示，本文選取10 m高度處風速為6 m/s時計算域內風場發(fā)展情況進行分析。

圖3 建筑物影響下的風場中風速分布云圖(風速v為6 m/s)

由圖3可見，建筑群的存在會極大地影響風場中風速的分布情況。在圖3(a)中，風速隨高度形成速度梯度，由于建筑物及大氣湍流特性的影響，整個風場的風速梯度呈現一種不均勻分布；另外結合圖3(b)可以發(fā)現，建筑物的迎風面會形成速度滯留區(qū)(如紅色實線矩形框所示)，風速大幅下降，而在建筑物的背面則會形成低速回流區(qū)(如紅色虛線矩形框所示)，此時氣流形成反向運動。對應的流場變化為風速分布情況的分析提供了依據。圖4為計算域內建筑物影響下的流場變化情況。

圖4 建筑物影響下的流場變化圖(風速v為6 m/s)

由圖4可見，建筑的阻塞作用使得風從建筑物側方及頂部進行繞流，形成反向運動的渦流，導致低速回流區(qū)的形成。本文選取圖4中典型的特征區(qū)域流場進行了局部放大分析。在圖4(a)中，當風與建筑物相遇時，在建筑物底部形成一個小回流區(qū)，由此往上氣流沿建筑物“爬上”頂部與主流相遇匯合，造成速度梯度上移，形成速度滯留區(qū)[見圖4(a)-Ⅰ]；當風離開建筑物頂端時，部分氣流脫體朝斜下流動折回地面，在建筑物空腔內產生與風速相反的反向封閉渦旋，促成低速回流區(qū)I形成[見圖4(a)-Ⅱ]；建筑物頂端氣流并未折回地面形成反向渦流[圖4(a)-Ⅲ]，由此可知頂部氣流運動并非是低速回流區(qū)Ⅱ形成的主要原因。在圖4(b)中，風側向繞過建筑物在背風面形成的反向渦團及渦對是低速回流區(qū)形成的重要因素；而風在離開建筑群后，流線恢復正常，風的運動也趨于穩(wěn)定。

2.2 天然氣管道泄漏初期天然氣泄漏擴散特性模擬

10 m高度處風速為6 m/s條件下天然氣管道泄漏初期天然氣氣云擴散云圖，見圖5。

圖5 天然氣管道泄漏初期天然氣泄漏擴散云圖(甲烷濃度為1%，風速v為6 m/s)

建筑物作為天然的障礙物，使得空間中的風場需要繞流而行，從而導致建筑物背面及建筑群空腔區(qū)域內形成湍流渦，影響天然氣氣云的擴散傳播。由圖5可見：當天然氣管道發(fā)生泄漏時，由于Coanda效應，天然氣氣云首先緊貼建筑物A向上擴散，且在天然氣氣云上升過程中，側向繞流形成的渦會影響天然氣氣云發(fā)展的結構，天然氣氣云頂部逐漸形成分叉結構(如紅色矩形框所示)；在泄漏時間t為90 s，當天然氣氣云抵達建筑物區(qū)域A頂端時，由于建筑物頂端風場繞流的作用，天然氣氣云開始平行向下風口傳播，另外建筑群空腔區(qū)域內的反向渦旋迫使天然氣氣云在下風口傳播過程中折回地面，并在建筑群空腔區(qū)域內積聚；隨著天然氣在建筑群空腔區(qū)域內積聚，天然氣氣云體積不斷增大，最終在風場的作用下，建筑群空腔區(qū)域內的天然氣氣云團將從建筑物B的側面及頂端繞行通過。

2.3 天然氣管道泄漏后期不同風速下建筑群內天然氣氣云積聚效應模擬

天然氣管道泄漏后期(泄漏時間t分別為500 s、550 s和600 s時)在順風向擴散距離內不同風速下天然氣濃度(體積分數)分布及氣云形狀變化，見圖6。其中，天然氣泄漏擴散云圖中的天然氣氣云濃度邊界為1%，天然氣氣云顏色表示x軸向速度。

當泄漏時間為500 s時，由圖6(a)可以明顯看出隨著下風口距離的增加，天然氣濃度逐漸降低，但在建筑群空腔區(qū)域內及建筑物B背風面附近，大風速條件下天然氣濃度明顯高于小風速條件，如在距離泄漏源為56.5 m處(即x=200 m時)，風速為2 m/s條件下該距離處天然氣濃度為10.23%，而風速為10 m/s條件下該距離處天然氣濃度卻能保持在18.81%，高于前者83.87%。結合對應的天然氣泄漏擴散云圖對上述現象進行分析，發(fā)現：當風速較小時(風速為2 m/s)，天然氣氣云從建筑物頂端通過建筑區(qū)沿下風口飄散，天然氣氣云整體顏色分布均勻，表示天然氣氣云在x軸向速度大體保持一致，穩(wěn)定在2～4 m/s范圍內；而隨著風速的提高(風速為6 m/s)，在建筑群空腔區(qū)域內，由于建筑物A頂端形成的反向渦旋導致天然氣氣云下沉，在建筑群空腔區(qū)域內積聚，之后天然氣氣云開始繞過建筑物B向下風口傳播。而由上述風場分析可知，建筑物B背風面在風場作用下形成了規(guī)模較大的渦對，對天然氣氣云產生較強的卷吸作用，導致天然氣氣云在建筑物B背風面發(fā)生二次下沉、積聚。通過分析天然氣氣云顏色分布可知，天然氣氣云整體產生較大的速度差值，下沉部分天然氣氣云x軸向速度較低，維持在3 m/s，而天然氣氣云上端仍保持較高的擴散速度，達到12 m/s左右。而當風速進一步增大時，即風速為8 m/s、10 m/s時，天然氣氣云將會牢牢地積聚在建筑群空腔區(qū)域內，停止向下風口擴散。與風速為6 m/s時相比，天然氣氣云上端x軸向速度及分布范圍都相應減小，保持在9 m/s。這是因為：大風速條件下建筑群空腔區(qū)域內及建筑物B背風面的渦團發(fā)展規(guī)模將會進一步增大，渦能更高，對于天然氣氣云的積聚效應也會更好，天然氣氣云整體受到更強的減速效果。

圖6 天然氣管道泄漏后期不同風速下建筑群內天然氣氣云積聚效應(甲烷濃度為1%)

不同風速下建筑區(qū)背風面渦團發(fā)展規(guī)模，見圖7。

圖7 不同風速下建筑區(qū)背風面渦團發(fā)展規(guī)模

當泄漏時間達到550 s時，由圖6(b)可知雖然天然氣氣云的擴散形狀并未發(fā)生明顯改變，但不同風速下建筑群空腔區(qū)域附近的天然氣濃度差距變得更加顯著，遠距離處天然氣濃度略有上升(如紅色實線框所示)，說明隨著泄漏時間的延長，小風速條件下(風速為2 m/s、4 m/s)天然氣氣云仍會持續(xù)向下風口傳播；當泄漏時間為600 s時，由圖6(c)可知除天然氣濃度差異明顯外，天然氣氣云積聚作用的影響區(qū)域也在擴大，當泄漏時間由500 s延長至600 s時，積聚區(qū)Ⅱ的x軸向最遠范圍由x=300 m擴張至x=350 m。

綜上分析可知，風速的提升會加劇建筑群區(qū)域內天然氣氣云的累積。

3 結 論

(1) 本文使用UDF(用戶定義函數)在Fluent軟件中對不同高度的平均風速與湍流參數進行了修改，并采用兩步模擬的方法，對風場進行了預先計算，獲得了穩(wěn)定的風場初始值。建筑物作為天然的障礙物，使得空間中的風場需要繞流而行，導致建筑物背風面及建筑群空腔區(qū)域內形成湍流渦，影響天然氣的泄漏擴散。

(2) 在天然氣管道泄漏初期，建筑群空腔區(qū)域內，氣流頂端繞行形成的反向封閉渦旋導致天然氣氣云在建筑群空腔區(qū)域內積聚。隨著天然氣氣云在建筑群空腔區(qū)域內積聚，天然氣氣云體積不斷增大，最終在風場的作用下，建筑群空腔區(qū)域內的天然氣氣云將從建筑物B的側面及頂端繞行通過。

(3) 在天然氣管道泄漏后期，建筑群與風場作用下流場運動對于天然氣會產生積聚作用，隨著風速的增大，繞流形成的渦團規(guī)模逐漸增加，其卷吸效果也逐漸增強，導致天然氣氣云整體在建筑區(qū)附近發(fā)生下沉。因此，天然氣管道泄漏事故發(fā)生時，根據風速條件和泄漏源附近的居民建筑物分布情況，可以判斷事故的主要影響區(qū)域及其嚴重程度，從而制定合理的應急預案。

(4) 隨著泄漏時間的延長，小風速條件下(風速為2 m/s、4 m/s)天然氣氣云仍會持續(xù)向下風口傳播；而大風速條件下(風速為6 m/s、8 m/s、10 m/s)天然氣氣云仍大量積聚于建筑群區(qū)域內，建筑物B背風面的積聚區(qū)Ⅱ范圍略微擴大，當泄漏時間由500 s延長至600 s時，該積聚區(qū)Ⅱ的x軸向最遠范圍由x=300 m擴張至x=350 m。