基于CALPUFF 模型的毒害性氣體泄漏事故三維動態(tài)模擬方法
——以開縣井噴事故為例

張溶倩 李梅 楊冬偶 劉暉

北京大學遙感與地理信息系統(tǒng)研究所, 北京 100871; ? 通信作者, E-mail: mli@pku.edu.cn

大氣擴散模型是一類研究污染物在大氣中擴散、運輸、轉化和沉降等過程的模型, 對解決工業(yè)發(fā)展帶來的污染物排放及危險品泄漏等問題具有重要意義。大氣擴散模型通常是獨立于 GIS 而在環(huán)境科學領域發(fā)展起來的, 在環(huán)境評價領域的應用較為廣泛, 對應用于應急管理中危險氣體擴散模擬及輔助決策模型的研究則相對較少, 因此對地形影響、三維擴散及計算速度等應急響應的特殊需求尚缺乏深入的研究。

近年來, 隨著我國經(jīng)濟建設的快速發(fā)展, 對油氣資源的需求越來越大, 伴隨著資源開采過程而來的危險事故和環(huán)境污染時有發(fā)生。我國含硫天然氣(體積濃度>1%)的儲量占全國天然氣儲量的 1/4,高含硫氣田的開發(fā)技術難度大, 風險高, 劇毒重氣體含量高, 一旦發(fā)生泄露事故, 就有可能帶來巨大的人員傷亡和財產(chǎn)損失。因此, 將大氣擴散模型應用于毒害性氣體泄漏事故的應急管理之中十分重要, 能夠為事故的應急預案及災后處理提供理論依據(jù)和決策支持。

目前已有學者對基于大氣擴散模型的應急管理系統(tǒng)進行研究和開發(fā), 不同研究中根據(jù)具體問題的需求選用了不同的大氣擴散模型。按基本原理擴散模型可以分為三大類: 數(shù)值模擬模型、高斯模型和拉格朗日模型。

數(shù)值模擬模型基于流體力學中納維爾–斯托克斯方程的大氣擴散偏微分方程, 常見的模型包括CMAQ, GEOS-Chem 和 Fluent 等 。 Feng 等[1]使 用CFD 模型建立地下空間有毒氣體泄漏風險評估模型, 并應用于廣州某地鐵中轉站。陳建國等[2]對開縣井噴事故進行數(shù)值模擬, 成功地模擬了有害氣體輸運的積聚和導向等現(xiàn)象; 練章華等[3]用 Fluent 軟件建立能反映實際硫化氫擴散過程的三維重氣流動和擴散的模型, 用 CFD 軟件模擬開縣事故源周邊的污染濃度變化。但是, 數(shù)值模擬方法計算精細、耗費大、計算時間較長, 不適用于應急情況下的實時計算。

高斯模型的基本原理是將大氣擴散過程遵循的基本偏微分方程——質量守恒方程在合理假設的基礎上推導轉化成一個簡單的數(shù)學公式, 稱為高斯煙羽分布公式, 常見的模型包括 ADMS[4], AERMOD[5]和 ALOHA 等。郭遵強[6]提出高斯模型與 GIS 組件 MapX 集成的應急系統(tǒng), 并對突發(fā)性大氣污染事故進行預測, 輸出結果主要用二維等值線和條形統(tǒng)計圖的方式表現(xiàn)。

拉格朗日模型的基本原理是追蹤并計算污染物粒子在確定性的風場及浮力和不確定性的湍流等共同作用下的軌跡常微分方程, 如 HYSPLIT[7]和 NAME等。學者們還提出一種利用拉格朗日軌跡來改進高斯模型的思路, 稱為拉格朗日–高斯煙團模型, 即將污染物按一定體積分割為若干個煙團, 使用拉格朗日方法計算煙團的軌跡, 煙團內部污染物的分布則使用高斯方法計算, 最終各個煙團進行疊加得到總濃度場, 常見的煙團模型包括 CALPUFF[8]和 RIMPUFF[9]等。S?rensen 等[10]提出基于拉格朗日方法的大氣應急模型 DERMA, 并與 ARGOS (accident reporting and guidance operational system)核事故決策支持系統(tǒng)進行集成。S?rensen 等[11]基于 RIMPUFF, 提出一種模擬病毒在大氣中擴散的應急模型。Li 等[12]開發(fā)了針對高含硫氣體應急響應的實時 GIS 平臺,實現(xiàn)在復雜地形下基于 CALPUFF 模型, 使用實時數(shù)據(jù)模擬氣體擴散的過程。謝穎斯[13]提出大氣擴散模型 HYSPLIT 和 CALPUFF 與 GIS 組件 DotSpatial集成的應急系統(tǒng)。本研究選用的 CALPUFF 模型屬于拉格朗日–高斯煙團模型, 往往用于較大尺度的實驗區(qū)域(大于 50 km), 由于該模型綜合考慮了復雜地形和氣象參數(shù)等多個方面, 擁有更高的計算精度[14], 被逐漸應用于應急響應 GIS 系統(tǒng)中。

目前使用大氣擴散模型進行案例模擬和應急預案的相關研究中, 對氣體擴散的結果輸出和可視化大多局限于傳統(tǒng)的二維方法, 缺少三維方法的研究。二維表達可用于模擬地表層的污染擴散情形,卻無法提供垂直方向的信息, 難以模擬真實大氣環(huán)境污染的空間層級結構或展現(xiàn)氣體在空間中擴散的過程和趨勢[15]。雖然也有研究采用將氣體濃度計算結果按高度分層疊加顯示的方式, 在局部小場景中達到較好的三維可視化效果[16], 但將各層割裂開來, 難以直觀地體現(xiàn)氣體濃度在整個空間的連續(xù)分布情況, 并且對地表高程以及影像的分辨率要求較高, 限制了在大尺度或平原場景中的使用。

針對當前毒害性氣體泄漏應急管理中存在的問題和需求, 本研究基于 CALPUFF 模型設計并實現(xiàn)毒害性泄漏氣體的三維濃度計算方法。基于面繪制方法實現(xiàn)大氣擴散模型的三維動態(tài)模擬?！?2·23”開縣特大井噴事故是我國油氣開發(fā)歷史上傷亡最重的安全事故之一, 造成深遠的社會影響。在該案例中, 川東北地區(qū)丘陵溝壑眾多, 需考慮復雜地形對泄漏氣體擴散的影響; 應急救援指揮中需要快速的實時計算, 及時提供參考。本文對開縣事故進行模擬和分析, 將計算結果與遙感影像、地形數(shù)據(jù)等進行集成, 實現(xiàn)擴散的時空過程模擬。

1 研究區(qū)概況

羅家寨、滾子坪氣田位于四川省達州市宣漢縣及重慶市開縣境內。該地區(qū)高含硫氣田的氣藏儲量豐富, 但地形復雜, 人口稠密, 氣藏壓力高, 硫化氫含量高, 氣藏埋藏深, 勘探開發(fā)難度極大, 安全風險和環(huán)保風險高。

2003年12月23日, 羅家寨氣田開縣(現(xiàn)重慶市開州區(qū))境內羅家 16H 井因操作人員違章造成井噴失控事故, 確認 243 人死亡, 4000 多人受傷, 8.3 萬多人受災。羅家 16H 為高含硫(硫化氫)氣井, 屬于中國石油天然氣集團公司西南油氣田分公司川東北氣礦。該氣井為中國石油天然氣集團公司作為水平井開采新工藝而設計的鉆井, 設計井斜深 4322 m,垂深 3410 m, 水平移位 1586.5 m, 水平段長 700 m。2003年5月23日開鉆, 設計日產(chǎn) 100 萬 m3[17]。如圖1 所示, 羅家 16H 井所在的高橋鎮(zhèn)曉陽村位于重慶市開縣西北方向, 距離縣城約 80 km, 氣井位于黃泥埡口附近被山丘包圍著的一個平臺上。該平臺上還有羅家 14 號、15 號和 16 號井, 礦井周圍 30 m范圍內有六七戶人家, 周圍 1000 m 米范圍內村舍林立, 約有數(shù)百戶[18]。

圖1 研究區(qū)位置Fig.1 Location of studied area

2 研究方法

2.1 基于 CALPUFF 模型的技術路線

本研究提出的面向應急響應的氣體擴散三維濃度場計算和可視化框架分為 3個主要模塊: 底層數(shù)據(jù)庫模塊、三維計算模塊和三維可視化模塊, 總體架構如圖2 所示。

圖2 氣體擴散三維計算及可視化流程Fig.2 Flow chart of three dimensional calculation and visualization of gas diffusion

底層數(shù)據(jù)庫以文件數(shù)據(jù)庫的形式存儲輸入的地形高程數(shù)據(jù)(DEM)、土地利用數(shù)據(jù)(LULC)、污染源參數(shù)和氣象數(shù)據(jù)(高空數(shù)據(jù) UP.DAT 和地面數(shù)據(jù)SURF.DAT)。

CALPUFF 模型是由美國西格瑪公司研發(fā)的新一代非穩(wěn)態(tài)氣相和空氣質量建模系統(tǒng), 主要由 4 部分組成: 地理和氣象預處理器、氣象模塊 CALMET、擴散模塊 CALPUFF 和后處理器 CALPOST。三維計算模塊基于 CALPUFF 模型, 經(jīng)過數(shù)據(jù)交換、數(shù)據(jù)預處理、氣象處理、濃度計算和后處理,得到三維氣體擴散濃度場。

數(shù)據(jù)交換接口實現(xiàn)模型與 GIS 之間的數(shù)據(jù)傳輸, 進行基本的投影轉換、數(shù)據(jù)檢查和單位轉換等規(guī)范化處理。數(shù)據(jù)預處理步驟中, 特定格式的地形高程文件、土地利用文件、地表氣象數(shù)據(jù)文件和高

空氣象數(shù)據(jù)文件經(jīng)過 TERREL 地形高程數(shù)據(jù)預處理器、CTGPROC 土地利用數(shù)據(jù)預處理器、MAKEGEO 地理數(shù)據(jù)最終預處理器、SMERGE 地面站氣象數(shù)據(jù)預處理器和 READ62 探空氣象數(shù)據(jù)預處理器等處理后, 可輸出 CALMET 模塊所需的格式化文件 GEO.DAT, SURF.DAT 和 UP.DAT 等。氣象處理步驟中向 CALMET 子模塊輸入以上文件及包含定義模型運行信息的控制文件 CALMET.INP, 輸出報告文件 CALMET.LST 及包含模型產(chǎn)生的小時網(wǎng)格氣象數(shù)據(jù)的無格式的磁盤文件 CALMET.DAT 或PACOUT.DAT, 建立模擬區(qū)域的三維格點風溫場。在濃度計算步驟中, 由 CALPUFF 子模塊讀入 CALMET 生成的氣象場及用戶輸入的污染源參數(shù), 循環(huán)修改控制文件 CALPUFF.INP 中的受體(需要計算濃度值的點位)層高度, 計算該層受體點的濃度及干濕沉降通量, 輸出 CONC.DAT 文件, 循環(huán)結束后得到網(wǎng)格化的三維濃度場。后處理步驟中, 由 CALPOST 子模塊對這些輸出文件進行處理, 整理為用戶指定的格式并進行能見度模擬。本研究將上一步生成的 CONC.DAT 轉化為n個格式化的 grd 文件,n為系統(tǒng)模擬的擴散時間(分鐘級別), 即第i個 grd 文件描述開始擴散后第i分鐘的污染氣體濃度分布。

三維可視化模塊將受體點數(shù)據(jù)讀入, 利用 MarchingCubes 面繪制算法追蹤, 生成指定層數(shù)和濃度值的等值面時序三維模型。將模型讀入三維軟件中, 同時疊加高程和地表影像, 實現(xiàn)擴散場景的三維動態(tài)可視化。

2.2 改進的氣體擴散模型多層計算

CALPUFF 目前只能一次性地輸出二維濃度場,因此許多研究只關注近地表的氣體濃度分布。實際上, CALPUFF 模型中受體點的高度信息在 CALPUFF 模塊控制文件 CALPUFF.INP 中的“非網(wǎng)格(離散)受體信息”組中指定, 通過修改這一參數(shù), 就可以實現(xiàn)對不同高度層的濃度場計算及多層結果的輸出。

本研究完善了氣體擴散模型的計算流程, 將二維計算變?yōu)槿S計算。改進后的流程可以根據(jù)用戶指定的層高及層數(shù), 在一系列空中層設置受體點,并結合氣體類型、泄漏速率和泄漏點距地面高度等參數(shù), 在這些層上執(zhí)行濃度計算, 得到多層濃度值,最終整合生成三維濃度場, 為后續(xù)可視化提供基礎數(shù)據(jù)。進行多層計算改進后 CALPUFF 模型的計算流程如圖3 所示。

圖3 基于CALPUFF的三維濃度場計算模型Fig.3 Three dimensional concentration field calculation model based on CALPUFF

在 CALPUFF 模塊程序中, 通過循環(huán)修改 CALPUFF.INP 文件的“非網(wǎng)格(離散)受體信息”組, 指定相應的高度參數(shù), 開展不同層計算, 得到逐層擴散濃度結果, 該結果不斷寫入網(wǎng)格化的濃度場文件CONC.DAT。在 CALPOST 后處理模塊, 將 CONC.DAT 整理成格式化的輸出文件, 在輸出路徑中, 創(chuàng)建文件夾“Layer0”至“Layerx”, 分別保存各層的氣體濃度時序結果。每層的文件夾內, 則從氣體泄漏開始時刻起, 以 5 分鐘為間隔, 按時刻依次創(chuàng)建文本文件, 達到時空數(shù)據(jù)結構化存儲的目的。算法流程如下。

算法1基于 CALPUFF 的三維濃度場計算算法

輸入網(wǎng)格化逐時風場文件 CALMET.DAT; CALPUFF 控制文件 CALPUFF.INP; CALPOST 控制文件CALPOST.INP; 受體層數(shù)N; 受體層高 deltah

輸出三維氣體濃度場文件夾Layer0至Layer(N–1)

1.for eachn∈[0,N) do

2.計算此受體層高度h=n* deltah

3.修改 CALPUFF.INP 文件, 更新各受體高度為h

4.向 CALPUFF 模塊輸入 CALPUFF.INP 和 CALMET.DAT, 執(zhí)行計算

5.修改 conc.dat 文件

6.end for

7.向 CALPOST 模塊輸入 CALPOST.INP 和 conc.dat,執(zhí)行后處理

8.輸出格式化后的文件列表Layer0至Layer (N–1)

2.3 氣體擴散時序數(shù)據(jù)三維重建方法

三維可視化的方法分為兩大類: 面繪制和體繪制。面繪制是對三維空間目標的表面進行繪制顯示, 其數(shù)據(jù)處理基礎是序列二維圖像數(shù)據(jù), 再結合其空間結構關系, 還原三維空間結構。面繪制主要分為體素級重建方法和切片級重建方法, 其中, 體素級重建方法的本質是在一個三維的數(shù)據(jù)場中生成等值面的過程[19], 其實現(xiàn)多基于經(jīng)典的 Marching Cubes (移動立方體)算法。本研究選用 Marching Cubes 算法對氣體擴散過程進行三維可視化。

Marching Cubes 算法由 Lorenson[20]等提出, 以體素為單位, 制造出想要重構的部分和復雜背景的分界[21], 然后從體素中抽取三角面片來擬合此邊界, 大量的三角面片構成一個等值面。等值面的構造經(jīng)過兩個主要計算過程: 1) 體素中三角面片逼近等值面的計算; 2) 三角面片各頂點法向量的計算[22]。Marching Cubes 算法的基本假設是沿著體素的棱邊,數(shù)據(jù)場呈連續(xù)線性變化。構建某個氣體濃度值的擴散等值面時, 當搜索到相鄰兩受體點處的濃度值分別大于和小于等值面值時, 在該條邊上必有且僅有一點是這條邊與等值面的交點。根據(jù)這一原理, 就可以判斷所求等值面與哪些體素相交。每個體素包含 8個頂點, 根據(jù)頂點處數(shù)據(jù)值大于或小于等值面的值, 將頂點標記為 1 或 0 兩種狀態(tài), 故每個體素有 28= 256種組合狀態(tài), 根據(jù)互補對稱性和旋轉對稱性簡化為 15種。得到面片后, Marching Cubes 方法采用中心差分計算出體素各頂點處的梯度, 再次通過線性插值求出三角面片各頂點的梯度, 即各頂點處的法向量, 然后選擇適當?shù)墓庹漳Ｐ瓦M行計算,生成具有真實感的圖形[20]。

但是, Marching Cubes 算法存在連接方式的二義性問題, 相鄰體素連接不當時, 可能出現(xiàn)拓撲不一致, 從而使等值面出現(xiàn)孔隙。解決二義性問題的方法主要有二義面平均值判定法、子構型查找表法、梯度一致性準則法和漸近線判別法。本研究使用子構型查找表法, 即對于 15種基本構型中具有二義性的構型, 各自建立一個子查找表, 每個子查找表包含兩種三角剖分方式。此外, 還需存儲一個表來記錄這些三角剖分方式的相容性[23]。

本文模型首先讀取 CALPOST 輸出的逐時和逐層的氣體濃度文本文件, 以點數(shù)據(jù)結構記錄受體網(wǎng)格的三維坐標, 以浮點數(shù)組記錄對應位置的濃度值,每個受體點即為立方體體素的一個頂點, 進而從三維濃度場追蹤提取等值面。算法 2 為使用 Marching Cubes 算法追蹤等值面的流程。

算法 2基于 Marching Cubes 的三維濃度場等值面追蹤算法

輸入受體點位置坐標 points; 濃度場數(shù)組 conc;受體網(wǎng)格維度 dims[3]; 等值面數(shù) contNum; 等值面值數(shù)組 contValue [contNum]; 體素中面片連接方式表triangleCases[256]

輸出等值面集合 contours

1.for eachk∈[0, dims[2]) do

2.for eachj∈[0, dims[1]) do

3.for eachi∈[0, dims[0]) do

4.由(i,j,k), (i+1,j,k), (i+1,j+1,k), (i,j+1,k), (i,j,k+1), (i+1,j,k+1), (i+1,j+1,k+1), (i,j+1,k+1)組成當前受體點體素的 8個頂點, 從 conc 中取出頂點處濃度值, 記入數(shù)組s[8]

5.if 該受體點處梯度未計算 then

6.調用中心差分法函數(shù)計算 8 頂點梯度, 記入數(shù)組gradient[8]

7.for eachn∈[0, contNum) do

8.當前追蹤的等值面值 value =contValues[n]

9.根據(jù) 8個頂點值 s[8]與 value的關系確定當前體素的索引下標值index

10.得到該體素中面片連接方式triangleCases[index]

11.構造該體素中的三角面片 triCase

12.for each edge∈triCase -> edges do

13.由中心差分法計算三角面片各頂點的坐標、濃度值、梯度值

14.將當前三角面片插入等值面集合contours中

15.end for

16.end for

17.end for

18.end for

19.end for

圖4(a)為某個時刻原始濃度場的一部分, 為 20行×14 列×11 層的矩陣。在該模擬設置中, 受體網(wǎng)格的水平間隔為 500 m, 垂直間隔為 20 m, 矩陣中心紅色點位處表示氣體濃度較高, 即污染源位置,藍色點位處則表示氣體濃度較低。圖4(b)為由該濃度場重建的等值面模型, 可以看出氣體由濃度最高的污染源處向東、向外層擴散, 濃度逐漸下降, 與原始濃度場表現(xiàn)的變化趨勢一致。另外, 追蹤氣體等值面時, 為保證三維表面的封閉性, Marching Cubes 算法會將原始數(shù)據(jù)范圍向外擴展, 最終氣體模型頂部達到地表以上約 600 m, 超出受體網(wǎng)格的垂直范圍, 因此該擴展部分的數(shù)值主要由地表起伏趨勢和臨近受體點的數(shù)值決定。

圖4 由原始濃度場追蹤得到等值面Fig.4 Isosurfaces from the original concentration field

為評估等值面追蹤的時間效率, 本研究針對同一氣體擴散案例, 改變受體網(wǎng)格分辨率, 對同一區(qū)域進行擴散過程模擬。選取一個特定濃度值的等值面進行構建, 對相關效率指標進行統(tǒng)計分析。表1選取擴散初期至后期的 9個時刻, 列出當時的非零值受體點數(shù)、構成等值面的頂點數(shù)、三角面片數(shù)以及等值面追蹤所需的時長。其中, 網(wǎng)格分辨率為 200, 500 和 1000 m 時的數(shù)據(jù)體維度分別為 100×100×11、40×40×11 和 20×20×11。程序運行的電腦配置為 Windows10 操作系統(tǒng), 16 G 內存, CPU采用 Intel(R) Core(TM) i7-1065G7 1.30 GHz, 基于Pycharm Community Edition 2020 平臺, 使用 Python語言編寫測試程序。

由表1 可知, 對于某一固定分辨率, 每個時刻的等值面構建時長主要受非零值受體點數(shù)影響, 隨著擴散時間增加, 模擬區(qū)域內濃度值超過零值的受體點數(shù)增加, 所需的構建時間也增加。當分辨率為200 m 時, 在 20~180 分鐘, 非零值受體點數(shù)由 110增至3355, 所需時長由 152.36 ms 遞增至998.68 ms。三角面片數(shù)與等值面具體形態(tài)相關, 對構建時長的影響不顯著。

表1 不同分辨率下等值面構建效率指標對比Table 1 Comparison of efficiency indexes of isosurface construction under different resolutions

當可用的非零網(wǎng)格點過少時, 得到等值面的形態(tài)可能失真, 故越高的網(wǎng)格分辨率越利于三維場景的重建。綜合考慮時間效率, 使用 200 m 分辨率時,每個時刻的構建時長基本上控制在 1 秒以內, 可以滿足應急場景的需求。故本研究使用 200 m 濃度場網(wǎng)格。

3 氣體泄漏事故模擬案例研究

3.1 實驗參數(shù)設置

事故發(fā)生于 2003年12月23日晚 21 時 55 分, 24日 15 時發(fā)生井噴的主管道成功封堵, 16 時放噴管線實施點火, 共持續(xù) 18 小時。當日該地氣溫為 4.6~8.0℃, 相對濕度為 94%~99%, 平均風速為 0.13 m/s,最大風速為 0.7 m/s, 風向西北偏西。根據(jù)現(xiàn)場調研,羅家 16H 氣井地層壓力高達 40 MPa 以上, 井口噴出的氣體富含 H2S, 預計無阻流量為 400×104~1000×104m3/d, 天然氣中 H2S 含量約為 151 mg/m3[24]。事故造成井場周圍居民和井隊職工共 243 人死亡, 上萬人因硫化氫中毒接受治療, 九萬余人疏散轉移, 經(jīng)濟損失上億元[25]。

本研究關注毒害性氣體泄漏后, 對污染源周圍人居環(huán)境和自然環(huán)境的影響。研究區(qū)地形起伏多變, 以污染源為中心的 50 km×50 km 范圍內海拔為50~1738 m, 平均 996 m, 因而氣體擴散情況復雜,垂直方向的變化明顯。我們在z軸方向設置從地表向上, 以 20 m 為間距的 11 層受體點。另外, 在突發(fā)氣體泄漏事故的應急場景中, 由于初始噴發(fā)速度高、噴發(fā)氣體量大, 有毒氣體濃度在時間和空間上的變化迅速, 常規(guī)大氣監(jiān)測中以小時為單位的模擬時間步長無法滿足應急場景的需求。一項關于 CALMET 時空分辨率對 CALPUFF 模擬濃度場影響的研究結果[26]顯示, 10 min 的時間步長即可達到較好的模擬結果, 能夠滿足實際應用的需求, 本研究選擇更精細的 5 min 時間步長。在 CALPUFF 中輸入的污染源、氣象條件及模擬設置參數(shù)如表2 所示。

根據(jù) H2S 氣體的性質, 本研究將追蹤濃度分為4 級, 由擴散外層向內依次為 0.01, 6.9, 15 和 100 mg/m3。表3 為氣體濃度超標閾限和對人體產(chǎn)生不同程度傷害的閾值。

表3 硫化氫對人體的影響及閾限Table 3 Effects and threshold of hydrogen sulfide on human body

3.2 模擬結果分析

3.2.1 三維時序擴散過程

將表2 中參數(shù)輸入計算系統(tǒng), 經(jīng)過模型計算和三維重建, 得到氣體泄漏后快速擴散階段的動態(tài)模擬過程(圖5), 分別表示事故發(fā)生后 15, 30, 45, 60,75, 90, 120, 150 和 180 分鐘的情形。

表2 模擬參數(shù)Table 2 Simulation parameters

在本次事故中, 氣井所在的開縣高橋鎮(zhèn)及相鄰的麻柳鄉(xiāng)、正壩鎮(zhèn)、天和鄉(xiāng) 4個鄉(xiāng)鎮(zhèn) 9.3 萬余人受災。其中, 受災最重的是高橋鎮(zhèn)的曉陽和高旺兩個村, 共 2419 人受到影響。90%的受影響群眾居住在距離氣井 2 km 外, 但由于暴露在高硫化氫環(huán)境中,90%的遇難者均位于氣井周圍 1.5 km 以內。由圖5 的模擬結果可以觀察到, 由于井噴位置處于低洼地帶, 且當天氣象條件較穩(wěn)定, 風速小且空氣濕度大, 擴散范圍總體上呈現(xiàn)偏心半球形, 有毒氣體從污染源水平向外和垂直向上擴散, 氣體半球的中心在偏西風作用下位于污染源下風的東南方向。從整體上看, 高濃度硫化氫聚集在污染源附近一定范圍內, 這些特點與實際受災情況基本上相符。

圖5 氣體擴散過程三維動態(tài)可視化Fig.5 Three dimensional dynamic visualization of gas diffusion process

3.2.2 典型時刻二維和三維對比

將氣體等值面三維模型導出, 疊加地形和影像,并與地表二維擴散結果進行逐時對比。

在事故發(fā)生 15 分鐘時, 污染源東南方向 500 m范圍內, 近地表的硫化氫氣體濃度已超過 15 mg/m3(圖6(a)中橙色區(qū)域)。曉陽村位于污染源西南方向約 500 m 處, 此時硫化氫濃度約 3 mg/m3; 高旺村位于其東偏南方向約 800 m 處, 濃度達到約 6 mg/m3。據(jù)調查, 離井口最近的曉陽村村民聽到井噴剛發(fā)生時的轟然巨響后, 不到 10 分鐘, 便聞到一股臭雞蛋味, 即已明顯超過嗅覺閾值。在井噴開始后約10 分鐘, 曉陽村和高旺村的硫化氫濃度均迅速從 0升至40 mg/m3[2-3]?？紤]到 CALPUFF 模型多用于長時間尺度的環(huán)境評價領域, 常規(guī)應用中一般采用模式默認的模擬方案和參數(shù)設置, 如逐時風場、默認的湍流參數(shù)化方案和地表參數(shù)等[27], 在應急情景中的模擬結果容易偏小, 但從數(shù)值量級的角度看,本方法的短時模擬結果仍然具有一定的參考意義。

事故發(fā)生 30 分鐘時, 污染源周圍 1000 m 以內的大部分區(qū)域硫化氫濃度已達到 15 mg/m3(圖6(c)中橙色區(qū)域), 該濃度已經(jīng)能夠對人眼產(chǎn)生刺激性,并伴隨明顯的異味。位于泄漏源東南方向約 500 米的高橋鎮(zhèn)已完全處于藍色區(qū)域, 根據(jù)基于 Fluent 的數(shù)值模擬結果[24], 井噴后半小時的硫化氫濃度已超過人的嗅覺閾值, 與本研究結果一致。

事故 1 小時后, 兩座村莊已經(jīng)完全被濃度約100 mg/m3的硫化氫氣體覆蓋。該濃度的硫化氫氣體會刺激人的眼睛及呼吸道, 甚至危害生命安全。從二維結果中可以觀察到, 近地表氣體被劃分為以污染源為圓心的兩個扇形區(qū)域。對應圖6(f)的三維地形可知, 在污染源北、西南、東南方向均為狹長山谷, 氣體從靠近西側山丘的山腳處噴發(fā), 同時在偏西風作用下, 沿山谷東側蔓延。由于周圍山體的阻擋, 高濃度硫化氫氣體在污染源附近的洼地中不斷聚集。另外, 從三維地形觀察到, 高旺村位于一座低山的西坡, 相對于污染源的海拔高度更高, 但此處整個山坡都被高濃度氣體覆蓋, 可見在有風、離污染源較近的情況下, 地勢較高處也非安全地帶。辛保泉等[28]通過風洞實驗也發(fā)現(xiàn), 泄漏源下游為山峰條件時, 氣體可能與山體碰撞或在山體背后的山坡上形成高濃度的累積。

圖6 事故發(fā)生15, 30, 60, 90分鐘的氣體濃度分布Fig.6 Gas concentration distribution at 15, 30, 60 and 90 min after accident

事故后 1.5 小時內, 污染源周圍 5 km 均已受到濃度超標的硫化氫氣體影響(圖6(g)中藍色區(qū)域)。該范圍也是本次事故中報道受災最重的區(qū)域, 居民均需要被疏散轉移。

3.2.3 擴散模型剖面分析

由圖6 可以看出, 三維空間模型可以完整、直觀地反映事故發(fā)展過程和危害影響范圍。為了呈現(xiàn)隨空間位置變化的更細致的氣體擴散特征, 可以配合垂直方向二維剖面圖輔助來綜合分析。圖7 為事故發(fā)生 60 分鐘時的擴散二維圖像, 沿線 AA′, BB′,CC′和 DD′得到下風向 0, 500, 1000 和 1500 m 處橫向擴散的氣體模型剖面圖(圖8)。

圖7 事故發(fā)生60分鐘的二維濃度分布圖及剖面線位置Fig.7 Two-dimensional concentration distribution and locations of profile lines at 60 min after accident

可以發(fā)現(xiàn), 在該時刻, 在下風向 1500 m 以內不同距離處, 氣體橫向擴散的范圍大致相同, 但不同濃度等級的圈層范圍逐漸變化。0 m 處, 即經(jīng)過泄漏源點的剖面中, 濃度 100mg/m3以上的占比很大(圖8(a)中紅色范圍), 表明此處災害嚴重, 而沿下風向則逐漸減輕, 橙色范圍占比增加。不同位置處剖面的輪廓形態(tài)也差異很大, 可以觀察到它們受地形起伏影響顯著, 高濃度氣體易在低洼處聚集。對比圖8(a)~(d)的地表剖面線, 每隔 500 m 的地形趨勢幾乎截然不同, 可見事故發(fā)生地地形復雜, 會給災害預測和應急救援帶來很大的不確定性和難度。

圖8 事故發(fā)生60分鐘的下風向不同位置處氣體剖面圖Fig.8 Gas profiles at different locations downwind at 60 min after accident

3.2.4 高橋案例總結分析

通過對開縣案例擴散過程的模擬分析, 可以得到如下結論。

1) 該案例中, 氣體的擴散主要受地形和氣象條件影響, 危害區(qū)域分布在順風的山坡前部或山谷洼地, 該區(qū)域不適宜居民居住, 一旦發(fā)生突發(fā)安全事故, 后果不堪設想, 應當在氣井選址和民房選址時予以考慮。高排放率的有毒氣體在小風速和低洼地形作用下迅速覆蓋污染源周邊, 并以偏心半球形態(tài)不斷擴大蔓延, 在高壓、高含硫且有風的情況下,污染源附近會受到高濃度毒害氣體影響, 向附近高地疏散并非安全的應對措施。

2) 污染源附近是發(fā)生事故后氣體濃度迅速上升且最大濃度持續(xù)出現(xiàn)的區(qū)域, 在川東北地區(qū), 常年平均風速較小, 且丘陵地貌不利于氣體的水平輸送和擴散, 人員撤離半徑可以酌情劃定為距離污染源 5 km。

4 結論與展望

本研究針對傳統(tǒng)的氣體擴散計算和可視化問題, 設計并實現(xiàn)基于 CALPUFF 模型的氣體濃度計算和三維可視化框架, 將模擬維度從近地表的二維擴展為動態(tài)的考慮垂直分布的空間三維加時間維度, 能夠更直觀、方便地進行應急模擬和事故發(fā)生機理的分析。將該方法應用到開縣井噴事故案例中, 驗證了其可行性和有效性, 并得到如下主要結論。

氣體泄漏后在不同高度層的擴散分布不同, 通過三維可視化的方式更易于從多個維度確定氣體擴散的范圍, 能夠對氣體的空間分布有連續(xù)的和直觀的認識; 疊加三維地形場景后, 便于分析各地的受災程度以及擴散形態(tài)的形成原因; 加入時間維度后,從動態(tài)擴散過程更能反映災情的發(fā)展變化趨勢, 對于應急事件的預演以及宏觀把握具有重要的指導意義。

本研究還存在一些不足之處, 未來應進行如下方面的研究。

1) 在模型方面, CALPUFF 模型是為長時間、大范圍的大氣質量模擬評估所設計, 應用于應急場景時存在一些難以處理的問題: 應急場景具有極強的突發(fā)性和危害性, 要求模型必須同時具有高計算速度和高計算精度, 而現(xiàn)有的 CALPUFF 等大氣擴散模型空間分辨率為十米級, 氣象輸入?yún)?shù)為秒級,時空分辨率不能自洽; 進行模擬參數(shù)設定時, 氣體排放率、溫度和風速等輸入?yún)?shù)為恒定值, 但在實際場景中, 這些條件往往在不斷的變化, 導致模擬結果與現(xiàn)實的差異。

2) 在三維計算方面, 多層受體、高網(wǎng)格分辨率的設定導致計算量增大, 加上 CALPUFF 模型考慮因素的復雜性, 計算時長通常在分鐘級。若要提高計算效率, 可以考慮應用 CPU 多核并行計算以及GPU 加速等技術。

3) 在體視化效果方面, Marching Cubes 算法存在不完善之處, 如數(shù)據(jù)量較小時等值面形態(tài)失真、追蹤范圍不夠精確等, 可以通過算法的改進得到修正。根據(jù)模擬目的的不同, 如關注事故對近地表居民區(qū)的影響, 關注毒害物質對山間動植物的危害、關注氣體擴散對高層大氣的污染, 等等, 可以選擇性地設定受體層的數(shù)量和層高, 聚焦于特定的三維空間范圍, 以期提高計算模擬效率。