煤礦采空區碳封存CO2泄漏地表擴散規律研究

丁 洋,陳文彬,林海飛,朱 冰,譚軍紅

(1.西安科技大學 安全科學與工程學院,陜西 西安 710054;2.西安科技大學 西部礦井開采及災害防治教育部重點實驗室,陜西 西安 710054)

0 引 言

近年來,全球氣候變化加劇與生態環境不斷惡化引起各國對碳減排的高度重視[1-3]。中國作為碳排放大國,2020年9月22日,習近平總書記在第七十五屆聯合國大會一般性辯論上鄭重宣布中國CO2排放力爭2030年前達到峰值,努力爭取2060年前實現碳中和[4]。實現碳中和的根本路徑是減排增匯,碳捕獲、利用和封存(Carbon Capture,Utilization and Storage,CCUS)被認為是脫碳的關鍵戰略,而CO2封存是CCUS鏈的最后一環,地質封存是目前發展最為成熟的一種封存方式[5-9]。

煤炭長期以來是中國的主體能源,在煤炭資源開采過程中,形成了大量的地下采空區[10-11],由于各種因素的限制,在煤礦采空區依然殘留大量煤炭,煤基質中的CO2吸附優勢要強于CH4,CO2注入采空區會與殘留煤炭中的CH4產生競爭吸附[12-15],從而實現CO2以吸附態封存于采空區。但采空區原始應力遭到破壞,裂隙演化和滲流規律變得復雜[16-18],這導致了采空區封存CO2相較于其他封存方式可能具有更高的泄漏風險,而一旦泄漏,則可能對植物、動物以及人類造成嚴重威脅。因此,明確采空區封存CO2泄漏至地表后的擴散規律和濃度分布具有重要的安全和環保意義。

對于CO2泄漏至大氣后的擴散規律,前人利用數值模擬進行了大量的研究。MAZZOLDI等將Kit Fox試驗和高壓管道泄漏CO2的擴散模擬進行比較,研究表明CFD模型的準確性比高斯模型更高[19]。XING等利用CFD軟件模擬了CO2井噴事故的擴散情況,并將試驗數據進行了對比驗證,得到標準k-ε和SSTk-ω模型的計算結果與CHANG標準下的試驗結果吻合較好[20]。WEN等在開源CFD代碼OpenFOAM上開發了一個專門用于CO2在大氣中擴散的求解器,模擬了在無地形影響下CO2的泄漏情況[21-22]。WANG等研究了地形和釋放源溫度對復雜地形條件下埋地管道CO2泄漏擴散的影響[23]。眾多學者在CO2地表泄漏擴散方面取得了豐碩的研究成果,對于CO2的擴散模擬,目前大部分研究是高壓儲罐和運輸管道的泄漏,對于常壓泄漏擴散機制的研究較少,并且在研究地形對CO2擴散的影響方面,多是用單一的地形起伏來進行表征,對于真實地形條件下的復雜影響研究較少。

文中根據試驗礦井采空區地表真實地形提取地形特征,建立CFD模型,利用Fluent模擬研究不同影響因素下常溫常壓CO2在封存區上方地表濃度分布規律和危險的擴散范圍,為采空區封存CO2項目安全評價提供參考。

1 地形特征提取和數值建模

1.1 試驗礦井位置

陜西省神府礦區歷經多年的開采,形成了大量采空區,以淺埋煤層為主,并且有低透氣性的紅土層保證了低壓儲存CO2的可行性。同時,榆林能源化工基地每年會產生大量的純CO2,基地與礦區處于同一區域,便于項目的開展,降低成本。因此,項目選擇陜西省侏羅紀煤田神府礦區某煤礦作為典型淺埋煤層采空區吸儲CO2實施地點。項目的中試位置擬選定在試驗礦井首采25201工作面,試驗平面范圍初步確定為300 m×500 m,其封存位置如圖1所示。

圖1 CO2封存區域Fig.1 CO2 storage area

1.2 地形特征提取和數值建模

為了觀測CO2泄漏后封存區域及周邊的CO2分布,選擇封存地點周邊1 500 m×1 500 m的區域作為模擬區域。地形提取流程如圖2所示,將所選區域的瓦片數據下載拼接后導入ArcGIS軟件中,結合封存地點的Dem數據提取地形特征,再根據中央經線進行投影,將投影后的數據導入Global Mapper軟件得到了研究區域的高程圖(圖3),該區域最高海拔1 225 m,最低海拔1 091 m,海拔高差134 m。最后在Sketch UP軟件中完成建模。

圖2 建模流程Fig.2 Modeling process

圖3 研究區域高程圖Fig.3 Elevation in study area

2 數值模擬

2.1 模擬方案

文中利用Fluent數值模擬軟件,分析了不同泄漏速度、不同風速下CO2在目標區域的擴散規律,并且在平坦無起伏的地形模型上進行了對照模擬,以便觀察分析地形影響。泄漏點設在圖3紅色圓圈區域內,該區域位于計算域的中部位置,更有利于觀測CO2泄漏擴散規律,并且注入井也位于封存區域的中部位置,泄漏點設在此處能模擬出井筒泄漏時的CO2擴散情況,同時該區域四周海拔都相對較高,便于研究地形的影響。設定泄漏口為直徑1 m圓孔,泄漏時間3 600 s,泄漏速度用CO2泄漏的質量流量來表征。模擬方案見表1。

表1 模擬方案Table 1 Simulation scheme

采空區封存CO2泄漏特點與管道或儲罐泄漏特點不同。管道或儲罐泄漏是由于高壓突然釋放,因此泄漏速度快、但并不持續;而采空區泄漏屬于常壓泄漏,泄漏速度緩慢、但持續泄漏。因此為了充分模擬這種泄漏特點,將泄漏速度根據梯度由小到大設置,并且對于每種情況,都模擬了從泄漏開始3 600 s后的CO2泄漏擴散結果。為了充分了解CO2泄漏的危害范圍,選擇較高的泄漏速度即10 kg/s作為初始條件模擬其他因素變化對CO2擴散的影響。根據調查,封存區域所在地風速常年小于5級(5級風速為8～10.7 m/s),因此選擇0～10 m/s的風速進行數值模擬。由圖3可知,研究區域總體海拔是沿著y軸正方向逐漸降低,因此為了研究地形和風速的相互影響,研究方案中將風速變化時的風向設為沿y軸負方向。

2.2 基本假設與流動方程

2.2.1 基本假設

對于采空區封存CO2泄漏釋放到地表的氣體,在抵達地表后,其壓力溫度已經和環境值基本一致,所以假設其泄漏壓力溫度與大氣壓和環境溫度相同,即常溫常壓泄漏,因此用于擴散模擬的源流體為氣態CO2。模擬研究基于以下假設：①氣體為連續介質;②氣體流動過程為不可壓縮;③氣體流動不考慮傳熱與化學反應;④地表組分傳輸模型只考慮空氣和CO2。

2.2.2 流動方程[24]

1)質量守恒方程。

(1)

2)動量守恒方程。慣性(非加速)參考系中的動量守恒由下式描述

(2)

(3)

式中μ為分子粘度,Pa·s;I為單位張量,右側的第二項是體積膨脹的影響。

3)能量守恒方程。

(4)

上述3個方程流體運動的基本控制方程,主要是將流體運動時必須遵循的物理定律用方程的形式表達出來,從而解出流體流動的未知量。

4)組分輸運模型。組分輸運模型是模擬混合物各組分之間或與其他相之間的相互作用,對于CO2泄漏擴散,主要是模擬CO2和空氣的相互作用。選擇求解組分的守恒方程時,Icepak通過解第i個組分的對流擴散方程來預測每個組分的局部質量分數Yi。這個守恒方程的一般形式如下。

(5)

式中Yi為物質i的質量分數,即質量擴散矢量,而Si為物質i的凈生產速率。

湍流中的質量擴散

(6)

該方程主要表征CO2在地表擴散的行為。其中Di,m是組分i的質量擴散系數,Sct是湍流施密特數,默認設置為0.7。

5)標準k-ε模型的傳輸方程。湍流是自然界普遍的流動形式,湍流運動的特征是在流動過程中流體各質點具有的互相混合的特性,速度和壓力等物理量在空間、時間上具有隨機性的脈動值,而湍流模型就是表征湍流運動特點的方程。湍流動能k及其耗散率ε由以下輸運方程獲得

(7)

(8)

在這些方程中,Gk為由于平均速度梯度產生的湍流動能;Gb為由于浮力產生的湍流動能,根據k-ε模型中浮力對湍流的影響計算。YM為可壓縮湍流中的脈動膨脹對總耗散率的貢獻。C1ε,C2ε,C3ε為常數;σk,σε分別為k和ε的湍流普朗特數;Sk,Sε為用戶定義的源項。

2.3 計算域和邊界條件

根據圖3高程數據建立采空區地表模型,通過紅圈所在位置確定泄漏口中心點坐標為(715,725,47.22)。為了適應地表134 m的高程差異,以模型在z軸上的點即(0,0,114.48)作為基準點,建立高度為300 m的立方體模型。模型計算域如圖4(a)所示。

圖4 計算域和局部網格特征Fig.4 Computational domain and local mesh features

圖4(b)顯示了擴散模型的面網格特征,由于幾何形狀比較復雜,計算區域主要離散為六面體單元,其中分別對地表和泄漏源附近的網格進行了不同程度的加密。在最終仿真之前,進行了網格獨立性研究,以確保網格對結果的影響可以忽略不計,最終確定的網格大約由400萬個單元組成。

擴散模型邊界條件設置如下：①CO2源：質量流量入口,環境壓力和溫度,質量流率根據模擬方案所述設定;②進風口：速度入口,環境壓力和溫度,風速根據模擬方案所述設定;③出口：具有環境壓力和溫度的壓力出口,模擬無風條件時,進風口也設為壓力出口;④地面和頂部：溫度等于環境溫度的防滑等溫墻。

2.4 模型驗證

利用上述k-ε模型和數值方法模擬了平地上泄漏速度10 kg/s,風速6 m/s的CO2擴散。WEN等利用SSTk-ω模型模擬了CO2和空氣的混合物從泄漏源以50 kg/s的泄漏速度、2 m/s的風速泄漏至大氣的擴散情況[21]。圖5比較了距離泄漏源下風向200 m處CO2濃度的試驗和兩者的模擬結果,k-ε和SSTk-ω模型的模擬結果比較表明,兩者的變化趨勢具有高度的一致性,k-ε模擬的較大風速造成了CO2到達200 m處的時間更早,SSTk-ω模擬的較大泄漏速度造成了同一時刻更高的CO2濃度,符合實際情況,因此本研究所用的k-ε模型具有較高的可靠性。

圖5 距離泄漏源下風向200 m處CO2濃度的試驗和模擬結果的比較Fig.5 Comparison of experimental and simulated results of CO2 concentration at 200 m downwind from the leakage source

3 結果與討論

研究的主要關注點是CO2濃度在10 000 ppm(即1%)水平的范圍。根據參考文獻[20]的研究結果,當CO2濃度超過1%時,可能引發人員出現氣喘、頭痛、眩暈等不適癥狀。為了分析CO2在濃度范圍為0%～1%時的擴散規律,繪制了CO2濃度的云圖,并將1%濃度水平的等值面作為對人體可能產生危害的范圍進行描述。這一范圍的描述包括在縱向(y軸方向)和橫向(x軸方向)的影響距離以及等值面的面積。需要明確的是,上述危害范圍僅表示對人體有害的范圍,并非CO2擴散的范圍(同樣適用于后續描述)。在豎直方向(即z軸方向)上,由于重力的作用,CO2的擴散距離較小,主要分布在地表空間內,并且相對穩定。因此,本研究不考慮豎直方向上的擴散,而將地表危害區域作為立體空間內的危害區域加以考慮。

3.1 平坦無特征地形的泄漏

為了評估真實地形對危害范圍的影響,提前模擬了平坦無特征地形上的CO2擴散作為對照組。圖6顯示了無風條件下泄漏速度為10 kg/s時,泄漏不同時間的CO2濃度分布云圖。無風時,CO2泄漏后往四周均勻擴散,擴散范圍逐漸增大,但高濃度區域相對于低濃度區擴散較為緩慢。

有風的條件下,整個模擬分為2個步驟進行：①穩態模擬：根據模擬方案提供的風速在地形上建立風場,為下一步的模擬提供初始條件;②瞬態模擬：在風場下模擬CO2在計算域內的泄漏擴散。圖7顯示了在平坦地形上,CO2在泄漏速度10 kg/s、風速6 m/s、風向沿y軸負方向時,泄漏不同時間的濃度分布云圖。此時CO2主要沿著風向擴散,在190 s左右,其計算域內的濃度分布云圖就趨于穩定,高濃度區域擴散相對于低濃度區域擴散也較為緩慢。

圖7 無風時不同泄漏速度下CO2泄漏不同時間的濃度分布(0%～1%)云圖Fig.7 Concentration distribution(0%～1%)of CO2 leakage at different time under different leakage velocity without wind

表2顯示了以上2種情況泄漏3 600 s后的1%濃度等值面數據,藍色區域代表可能對人體產生危害的區域。由表可見,無風和有風條件下的影響面積存在顯著差異,這是因為無風時CO2均勻擴散,有風時CO2沿風向擴散,導致縱向影響距離增加,橫向影響距離縮小,并且風力加速了CO2的擴散,使CO2累積變得困難,所以相比無風條件高濃度區域面積急劇縮小。

3.2 真實地形條件下泄漏速度的影響

圖7顯示了無風時不同泄漏速度下CO2濃度為0%～1%的擴散云圖。在無風條件下,CO2泄漏后呈不規則分布,泄漏速度越大,CO2的橫向和縱向擴散范圍越大。泄漏速度較小時,CO2主要往地勢較低的區域擴散,當CO2泄漏速度足夠大時,其擴散也能克服部分地形影響,往地勢高的區域擴散,但相對于地勢低的方向,擴散有限?？傮w而言,CO2擴散的趨勢是一致的,泄漏速度的變化只會影響擴散范圍的大小。無風條件下,CO2擴散主要受2個因素的影響,一是地形因素,由于CO2是重氣,在重力的作用下會下沉,因此會往地勢低的區域擴散;二是濃度因素,CO2會由濃度高的區域往濃度低的區域擴散。

在相同的泄漏速度和擴散時間內,平坦地形的CO2擴散范圍大于真實地形,因為CO2在平坦地形是從泄漏點往四周均勻擴散,而在真實地形CO2主要沿著y軸正方向運動,地勢的升高限制了其他方向的擴散。然而,在y軸負方向,真實地形的擴散距離大于平坦地形。

表3顯示了無風條件下CO2泄漏3 600 s后濃度為1%的等值面數據。由表可見,隨著泄漏速度的增加,CO2擴散的橫縱向影響距離以及影響面積都逐漸增加,因此泄漏速度是影響CO2擴散的關鍵因素,直接決定了CO2泄漏之后的危害范圍。當泄漏速度為10 kg/s時,與平地相比,真實地形下的影響面積更大,說明在無風條件下封存區域地表整體地形更有利于CO2的局部累積,從而造成CO2濃度1%以上的區域面積較大,這一點在預防真實地形發生CO2泄漏方面需要引起重視。

表3 無風情況下真實地形泄漏3 600 s后的1%CO2濃度等值面數據Table 3 Isosurface data of 1% CO2 concentration after 3 600 s of real terrain leakage under no wind condition

3.3 真實地形條件下風速的影響

圖8顯示了在泄漏速度為10 kg/s,風向沿y軸負方向時,不同風速下的CO2濃度分布云圖。在CO2云圖前段,結合圖3高程圖所示,CO2有明顯的沿山谷擴散的情況,這是由于擴散云圖前段風向與山谷走向幾乎一致,但在后段山谷走向發生了偏移,此時,風為CO2提供的力克服了地形影響,CO2不再沿山谷擴散。在可見區域(計算域)內,CO2在風力作用濃度分布趨于穩定,即隨時間變化幅度很小或不變,并且風速越大,其達到穩定越快。

圖8 泄漏速度10 kg/s時不同風速下CO2泄漏不同時間的濃度分布(0%～1%)云圖Fig.8 Concentration distribution(0%～1%)of CO2 leakage at different time under different wind speeds at the leakage rate of 10kg/s

風速2 m/s時,在近泄漏源處,高濃度區域呈現沿x軸負方向偏移的傾向,這是由于地形的影響,從圖3中可以看出,在該區域右側是陡峭的山坡,上下落差超過50 m,所以阻擋了大部分CO2沿x軸正方向的橫向擴散。當風速逐漸增大后,以上情形不再出現。風速相同時,CO2在真實地形的擴散達到穩定狀態所需的時間比平坦地形更長,因為當風沿y軸負方向吹時真實地形下的CO2擴散需要克服地形對它的阻礙作用,而平坦地形則不需要。

表4顯示了不同風速下CO2泄漏3 600 s后濃度為1%的等值面數據。由表可知,風速增大,CO2擴散的危害范圍并沒有隨之擴大,其橫向影響距離隨風速的增大而減小,但其縱向影響距離是先增大后減小,這可能是由于風速到達一定值時,會加速CO2與空氣的混合,從而造成了泄漏CO2的快速稀釋,使高濃度區域的覆蓋范圍減小。

表4 有風情況下真實地形泄漏3 600 s后的1%CO2濃度等值面數據Table 4 Isosurface data of 1% CO2 concentration after 3 600 s of real terrain leakage under windy conditions

4 結 論

1)試驗礦井采空區封存CO2的地表泄漏速度越大,擴散范圍越大,危害范圍也越大;風速越大,擴散速度越快,橫向擴散越小。

2)在無風和風速較小(小于2 m/s)時,CO2會主要沿著地勢低的方向擴散。

3)當風速大于等于2 m/s時,風和地形對CO2擴散都具有重要作用,并且隨著風速逐漸增加,CO2的擴散會慢慢受風的主導,因此一旦發生泄漏需要注意泄漏源下風向的CO2濃度變化。