孔隙率對(duì)埋地超臨界CO2管道泄漏擴(kuò)散的影響

陳兵， 趙瓊， 郭煥煥

(西安石油大學(xué)機(jī)械工程學(xué)院， 西安 710065)

全球變暖對(duì)人們生活和生產(chǎn)的影響與日俱增，將生活及工業(yè)產(chǎn)生的廢氣捕集凈化后，加工提取其中的CO2并輸送至目的地利用和封存，這一高效減排CO2技術(shù)稱為CCUS(carbon capture、use and storage)。在CCUS過(guò)程中，碳運(yùn)輸是一個(gè)重要的中間環(huán)節(jié)，管道運(yùn)輸因其輸送量大、可長(zhǎng)期持續(xù)輸送等眾多優(yōu)點(diǎn)，成為最經(jīng)濟(jì)可行的方式[1-2]。為獲得更高的傳輸效率，管道輸送的CO2一般處于超臨界相態(tài)，超臨界態(tài)的CO2密度與液相密度接近，黏度系數(shù)非常小，接近于氣相且具有非常高的擴(kuò)散性[3-4]。管道的施工問(wèn)題或水合物形成而產(chǎn)生的局部腐蝕等原因易導(dǎo)致管道發(fā)生泄漏，泄漏的CO2會(huì)下沉到地面并聚集在低凹處，過(guò)高的CO2濃度會(huì)導(dǎo)致人和動(dòng)物窒息[5]。

國(guó)外在CO2管道方面的模擬研究與實(shí)際應(yīng)用已經(jīng)十分成熟。Cortis等[6]采用在穩(wěn)定狀態(tài)下的Navier-Stokes方程，建立二維系統(tǒng)研究CO2封存點(diǎn)泄露擴(kuò)散的密度效應(yīng)對(duì)其在大氣中擴(kuò)散的影響。Wen等[7]使用FOAM(field operation and manipulation)求解器，采用均勻的弛豫模式，基于湍流動(dòng)能k和耗散率ω的k-ω湍流模型對(duì)埋地低雷諾數(shù)流體管道穿孔時(shí)釋放致密相CO2并在大氣中的擴(kuò)散進(jìn)行模擬研究，并將模擬結(jié)果應(yīng)用于實(shí)際工程中。針對(duì)埋地CO2管道泄漏這一課題的探討，中國(guó)仍然處于發(fā)展階段，實(shí)驗(yàn)及模擬的相關(guān)成果也逐漸呈現(xiàn)。李康[8]通過(guò)自主搭建的地面超臨界CO2管道泄漏擴(kuò)散實(shí)驗(yàn)平臺(tái)，采用理論模擬與實(shí)驗(yàn)比對(duì)相結(jié)合的方法，研究不同影響因素對(duì)泄漏口附近參數(shù)的影響。劉正剛[9]通過(guò)自主搭建的三套CO2埋地管道輸送裝置，對(duì)發(fā)生穿孔和斷裂的CO2管道在氣相和液相這兩種輸送相態(tài)下，通過(guò)實(shí)驗(yàn)驗(yàn)證泄漏口徑和泄漏口壓力等對(duì)泄漏和泄放的影響。陳國(guó)龍等[10]認(rèn)為現(xiàn)有的泄漏事故均為小孔泄漏(泄漏孔徑小于20 mm)，搭建了泄漏孔徑為1 mm和3 mm的實(shí)驗(yàn)臺(tái)。郭曉璐等[11]對(duì)國(guó)內(nèi)外泄漏擴(kuò)散進(jìn)行了整理，認(rèn)為初期泄漏為小孔泄漏，持續(xù)發(fā)生泄漏則會(huì)逐漸變?yōu)榇罂仔孤?泄漏孔徑在20 mm至管徑之間)。劉少榮[12]根據(jù)工程多發(fā)事故對(duì)管道進(jìn)行了孔泄漏實(shí)驗(yàn)，得到了不同泄漏孔徑下的泄漏規(guī)律。孔泄漏和橫向斷裂是長(zhǎng)輸管網(wǎng)的兩種泄漏方式，但在實(shí)際中，管道一旦發(fā)生泄漏導(dǎo)致工況不穩(wěn)就會(huì)被監(jiān)測(cè)并搶修，所以工程上CO2管道泄漏事故多為小孔泄漏。

中國(guó)對(duì)CO2管道泄漏擴(kuò)散的研究大多停留在氣相和液相短程架空管道的實(shí)驗(yàn)及模擬。而輸送量大、輸送效率高且溫度壓力更適合CCUS工程要求的超臨界CO2長(zhǎng)輸管道相關(guān)研究和應(yīng)用還沒(méi)有深入展開(kāi)。由于超臨界CO2管道基本為埋地管道，只有少數(shù)管段會(huì)因地理因素架空，所以針對(duì)埋地CO2管道泄漏過(guò)程的擴(kuò)散規(guī)律展開(kāi)研究在工程上具有更高的實(shí)際價(jià)值。根據(jù)國(guó)外學(xué)者的研究方法，基于標(biāo)準(zhǔn)k-ε高雷諾數(shù)模型，建立埋地管道三維模型，與中國(guó)相關(guān)的實(shí)驗(yàn)數(shù)據(jù)對(duì)比驗(yàn)證模型的正確性，對(duì)超臨界CO2埋地輸送管道因腐蝕等因素導(dǎo)致管道發(fā)生初期泄漏(小孔泄漏)后，模擬土壤孔隙率對(duì)埋地管道不同截面處CO2濃度擴(kuò)散規(guī)律的影響。

1 理論基礎(chǔ)與控制方程

以流體力學(xué)[13]為基本理論，根據(jù)多孔介質(zhì)理論建立埋地管道發(fā)生小孔泄漏空間模型[14],在此基礎(chǔ)上利用GAMBIT軟件采用疏密結(jié)合的方式建立網(wǎng)格；結(jié)合流體力學(xué)基本方程以及埋地管道的Fick擴(kuò)散模型，使用專業(yè)模擬軟件FLUENT選用湍流模型中的Realizablek-ε模型[15-18]以及合適的控制方程并設(shè)置合理的邊界條件進(jìn)行數(shù)值計(jì)算，經(jīng)迭代收斂得到模擬結(jié)果。

1.1 多孔介質(zhì)

1.1.1 孔隙率

由于土壤這一多孔介質(zhì)的內(nèi)孔隙參差錯(cuò)落，故采用孔隙率來(lái)度量其空隙性質(zhì)，即

(1)

式(1)中：Vp、Vt為各孔隙的體積和整體單元體積，并默認(rèn)在垂直方向孔隙率是均勻的，以便于數(shù)值模擬計(jì)算。

1.1.2 達(dá)西速度

可用達(dá)西定律來(lái)描述一維流體運(yùn)動(dòng)，即

(2)

式(2)中：vx為x軸方向的達(dá)西速度；K為流體滲透率；μ為流體的黏度；p為流體壓力。

1.1.3 多孔介質(zhì)的狀態(tài)方程

只考慮彈性變形范圍之內(nèi)的情況，當(dāng)壓力差不大時(shí)，孔隙率變化的狀態(tài)方程為

φ=φ0(1+CφΔp)

(3)

式(3)中：φ0為參考孔隙率；Cφ為孔隙壓縮系數(shù)；Δp為相鄰孔隙之間的壓力差。

1.2 流體力學(xué)的基本理論

埋地CO2管道泄漏擴(kuò)散的影響因素有流速、壓力、密度以及溫度，所以在直角坐標(biāo)系下適用于CO2流體的三大守恒定律。

1.2.1 連續(xù)方程

可壓縮流體在多孔介質(zhì)中的連續(xù)性方程為

(4)

式(4)中：t為時(shí)間；ρ為密度；v為速度矢量。

1.2.2 運(yùn)動(dòng)方程

將所研究的CO2流體視作可壓縮流體，在忽略質(zhì)量力的情況下，流動(dòng)氣體的N-S(Navier-Stokes)方程為

(5)

式(5)中：ν為速度矢量；μ′與μ分別為體積與動(dòng)力黏性系數(shù)；D為變形張量。

1.2.3 能量方程

要使流體在多孔介質(zhì)中擴(kuò)散的模擬過(guò)程遵循物理規(guī)律，應(yīng)使用能量方程加以限制，即

(6)

式(6)中：e為內(nèi)能；f為質(zhì)量力矢量；P為壓力張量；λ為熱傳導(dǎo)系數(shù)；T為絕對(duì)溫度；q為輻射熱量。

1.2.4 狀態(tài)方程

將CO2流體的壓縮性與多孔介質(zhì)中的空腔結(jié)構(gòu)相結(jié)合，得到合適的狀態(tài)方程。

壓縮系數(shù)為

(7)

當(dāng)壓力差不大時(shí)，式(9)可表示為

ρ=ρ0(1+cfΔp)

(8)

式中：V為流體體積；cf為流體壓縮系數(shù)；ρ0為參考密度。

1.3 埋地管道的Fick擴(kuò)散模型

根據(jù)Fick擴(kuò)散模型，氣體在土壤中的非穩(wěn)態(tài)泄漏擴(kuò)散濃度為

(9)

式(9)中：C為擴(kuò)散物質(zhì)的濃度，mol/L；M為擴(kuò)散物質(zhì)的質(zhì)量，kg；D為擴(kuò)散系數(shù)，m2/s；x為擴(kuò)散層的深度，m；t為物質(zhì)擴(kuò)散的時(shí)間，s。

1.4 k-ε湍流模型

由于CO2黏滯系數(shù)較低可以忽略，可以假設(shè)流體在多孔介質(zhì)中為完全湍流。故選用廣泛應(yīng)用于完全湍流過(guò)程的k-ε標(biāo)準(zhǔn)模型，其湍動(dòng)能k和耗散率ε方程形式為

(10)

(11)

式中：速度梯度與浮力對(duì)k的影響系數(shù)分別為Gk和Gb；YM為由脈動(dòng)膨脹引起的ε的變化；μt為湍流黏性系數(shù)；1/σk、1/σε為k、ε的湍流普朗特?cái)?shù)，C1ε、C2ε、C3ε為常參數(shù)。

1.5 建立理想化模型并設(shè)置邊界條件

FLUENT軟件在流體力學(xué)基礎(chǔ)理論的支持下集成了多種解決湍流、傳熱與相變、多相流等理論模型，可以利用有限體積法等算法計(jì)算模擬復(fù)雜流體的動(dòng)態(tài)過(guò)程，通過(guò)其組分輸運(yùn)模型可以得到流體各組分的擴(kuò)散特性，所以軟件非常適用于土壤中泄漏以及擴(kuò)散的研究。合理建立理想化模型有利于軟件的定量計(jì)算。

土壤簡(jiǎn)化為各向同性的多孔介質(zhì)，孔內(nèi)為標(biāo)壓空氣； 由于小孔泄漏過(guò)程所以忽略泄漏過(guò)程對(duì)土壤結(jié)構(gòu)的改變；由于多孔介質(zhì)導(dǎo)熱系數(shù)很低，小孔泄漏規(guī)模較小，故忽略氣體與土壤之間的熱交換；管道輸送CO2濃度為99%。

在理想模型的基礎(chǔ)上采用疏密結(jié)合的方式進(jìn)行網(wǎng)格劃分，其入口邊界條件、出口邊界條件以及管壁邊界條件的具體設(shè)置如表1所示。

表1 模型邊界條件

將網(wǎng)格參數(shù)導(dǎo)入FLUENT，考慮重力影響因而啟動(dòng)非定常求解器，由于處理過(guò)程為完全湍流過(guò)程故采用k-ε湍流模型，開(kāi)啟組元輸運(yùn)即可追蹤C(jī)O2的泄漏擴(kuò)散過(guò)程。以土壤顆粒直徑Dp=1.6×10-3mm，孔隙率φ=0.45為例，根據(jù)Ergun公式計(jì)算可得黏性阻力和慣性阻力系數(shù)分別為C1=101 035 302.91 m-2，C2=13 203.02 m-1。

2 模擬方案可行性驗(yàn)證

劉正剛[9]將發(fā)生小孔泄漏的氣體CO2管道埋入沙箱中，經(jīng)過(guò)實(shí)驗(yàn)得到了管道發(fā)生泄漏時(shí)，不同因素的影響下CO2的擴(kuò)散規(guī)律。結(jié)合上述理論模型開(kāi)展與沙箱實(shí)驗(yàn)相對(duì)應(yīng)的仿真模擬過(guò)程，模擬結(jié)果與實(shí)驗(yàn)結(jié)果相比對(duì)，從而驗(yàn)證所建模型的可靠性。

2.1 數(shù)值模擬結(jié)果

設(shè)置好管道泄漏擴(kuò)散模型的入口邊界條件、出口邊界條件以及管壁邊界條件，選定設(shè)置的基本類型，在表2所示工況下，監(jiān)測(cè)泄放口上方地表的固定監(jiān)測(cè)點(diǎn)C1(距泄放口水平距離和豎直距離分別為0.2 m和0.15 m)的CO2濃度(體積分?jǐn)?shù))隨時(shí)間的變化情況，并對(duì)比在不同泄漏口徑下模擬數(shù)據(jù)以及實(shí)驗(yàn)測(cè)量數(shù)據(jù)，對(duì)比結(jié)果如圖1所示。

表2 實(shí)驗(yàn)工況參數(shù)

圖1 不同泄漏口徑下C1點(diǎn)二氧化碳濃度隨時(shí)間的變化Fig.1 Changes of CO2 concentration at C1 point under different leakage caliber with time

2.2 誤差計(jì)算

將圖1所示的實(shí)驗(yàn)數(shù)據(jù)與模擬數(shù)據(jù)進(jìn)行比對(duì)，分析其中誤差，結(jié)果如表3所示。

表3 實(shí)驗(yàn)、模擬數(shù)據(jù)對(duì)比

2.3 結(jié)果分析

由圖1和表3可得實(shí)驗(yàn)和模擬在不同因素作用下固定點(diǎn)測(cè)得的CO2濃度隨時(shí)間變化曲線是相接近的。實(shí)驗(yàn)時(shí)泄漏產(chǎn)生的膨脹波會(huì)使CO2濃度發(fā)生波動(dòng)，膨脹波會(huì)隨著泄漏孔徑的增大而更加明顯、更加劇烈。但本文中模擬的是小孔泄漏的穩(wěn)態(tài)擴(kuò)散過(guò)程，CO2濃度曲線較為平穩(wěn)，故認(rèn)為模擬數(shù)據(jù)和實(shí)驗(yàn)結(jié)果相比偏差較小，因而反映了本文所建立模型的正確性。

3 孔隙率對(duì)泄漏CO2擴(kuò)散的影響

3.1 建立埋地管道模型

3.1.1 輸送CO2管道參數(shù)

以中國(guó)某油田超臨界CO2埋地管道為研究對(duì)象，將其輸送工況參數(shù)轉(zhuǎn)化為模擬實(shí)驗(yàn)?zāi)Ｐ偷膮?shù)，數(shù)據(jù)如表4所示。

3.1.2 管道-小孔泄漏模型

研究小孔泄漏初期擴(kuò)散分析[19]的具體模型如圖2所示。如圖3所示，結(jié)合圖2以及表4中的管輸參數(shù)建立了4 m×3 m×3 m的埋地輸氣管道小孔泄漏模型，將直角坐標(biāo)系的XOZ面設(shè)置為地表，X軸的正反向?yàn)镃O2輸送方向，則泄漏口的坐標(biāo)為(0，-1.8，0)m。

表4 CO2輸送管道參數(shù)匯總表

圖2 埋地管道小孔泄漏擴(kuò)散模型Fig.2 Model of small hole leakage diffusion in buried pipeline

3.1.3 網(wǎng)格劃分

采用GAMBIT軟件網(wǎng)格劃分埋地管道三維模型，計(jì)算基于穩(wěn)態(tài)條件地表CO2濃度在不同單元數(shù)網(wǎng)格下的擴(kuò)散情況，驗(yàn)證網(wǎng)格的無(wú)關(guān)性，結(jié)果如圖4所示。

圖4 網(wǎng)格無(wú)關(guān)性驗(yàn)證Fig.4 Grid independence verification

由圖4可知，四種單元數(shù)網(wǎng)格中196萬(wàn)網(wǎng)格數(shù)量級(jí)就可滿足模擬要求，網(wǎng)格劃分示意圖如圖5所示。

圖5 網(wǎng)格劃分示意圖Fig.5 Map of grid

3.2 模擬分析孔隙率對(duì)CO2泄漏擴(kuò)散的影響

3.2.1 擴(kuò)散云圖

埋地管道一旦發(fā)生泄漏其輸送物質(zhì)的擴(kuò)散介質(zhì)就是土壤，中國(guó)土壤資源十分豐富且種類繁多，土壤因疏松程度的差異可分為砂土、黏土和壤土，因滲流特性的不同也可分為均勻土壤和非均勻土壤。在模擬中把土壤視作是各向同性的多孔介質(zhì)，其疏松程度用孔隙率來(lái)表征[20]，具體數(shù)據(jù)如表5所示。

表5 砂土類別

由中國(guó)土壤的特點(diǎn)選取孔隙率φ分別為0.35、0.45、0.55、0.65，導(dǎo)入表4中的數(shù)據(jù)。各個(gè)孔隙率下截面XOY的 CO2濃度云圖在t=100、200 s時(shí)如圖6、圖7所示，在時(shí)刻t=400、500 s如圖8、圖9所示。

圖7 t=200 s時(shí)XOY截面CO2濃度云圖Fig.7 The CO2 concentration cloud of XOY section when t=200 s

圖8 t=400 s時(shí)地表CO2濃度云圖Fig.8 The surface CO2 concentration cloud when t=400 s

(1)由圖6、圖7所示，在相同泄漏擴(kuò)散時(shí)間內(nèi)孔隙率的增大使得高濃度范圍隨之增大，CO2會(huì)以更短的時(shí)間從土壤擴(kuò)散到地表；隨著泄漏擴(kuò)散時(shí)間的增加，土壤中CO2濃度云圖以半圓形逐步向外擴(kuò)散，由于Y軸上CO2擴(kuò)散速率最大，故監(jiān)測(cè)軸線上CO2的濃度便可測(cè)得CO2濃度距地表距離。

(2)從圖8、圖9所示，擴(kuò)散到地面的CO2會(huì)以平行于地表且泄漏口為圓心按照中間高邊緣低的濃度梯度呈圓形分布，在相同泄漏擴(kuò)散時(shí)間內(nèi)，圓形區(qū)域的直徑隨著孔隙率的增大而增大；隨著泄漏擴(kuò)散時(shí)間的增加，圓形區(qū)域中的高濃度范圍也越來(lái)越大，確定濃度在Y=0軸線上的擴(kuò)散距離即可求出危險(xiǎn)區(qū)域的范圍。

3.2.2 模擬結(jié)果及分析

濃度高于5%的CO2足以在短時(shí)間內(nèi)致人窒息，故將5%確定為危險(xiǎn)臨界濃度[21]。處理云圖數(shù)據(jù)并繪制如圖10所示曲線。

繪制各孔隙率Y軸上臨界濃度隨時(shí)間變化曲線如圖10所示，可以看到，曲線可以明顯分為兩段，即地面下0～1.0 m段5%濃度節(jié)點(diǎn)位置與孔隙率的相關(guān)性很好且與預(yù)期相符，即土壤越松散擴(kuò)散速率越大，CO2濃度達(dá)到臨界值的時(shí)間越短；在地面以下1.0～1.8 m段5%濃度節(jié)點(diǎn)位置與孔隙率的相關(guān)性并不符合預(yù)測(cè)規(guī)律。

考慮不符合規(guī)律的1.0～1.8 m段，繪制該段不同時(shí)刻Y軸上壓力分布曲線如圖11所示。

分析圖11可得，地面以下1.0～1.8 m段在泄漏發(fā)生的初始時(shí)刻高壓使CO2在孔隙中以射流的形式擴(kuò)散，由于多孔介質(zhì)結(jié)構(gòu)的會(huì)降低射流的方向及動(dòng)量，故在泄漏口處隨著壓力分層而出現(xiàn)負(fù)壓區(qū)域；壓力在泄漏后一段時(shí)間內(nèi)變小并接近大氣壓，云圖呈均勻擴(kuò)散。

圖9 t=500 s時(shí)地表CO2濃度云圖Fig.9 The surface CO2 concentration cloud when t=500 s

圖10 各孔隙率Y軸上5%濃度節(jié)點(diǎn)隨時(shí)間變化曲線Fig.10 Changes of 5% concentration nodes in each porosity Y axis with time

圖11 不同時(shí)刻Y軸上的壓力分布圖Fig.11 Distribution of pressure on Y axis at different times

圖12 各孔隙率下地表5%等濃度線距Y軸線 距離與時(shí)間的關(guān)系Fig.12 Correlation between the distance of the surface 5% isoconcentration line to Y axis and diffusion time under each porosity

各孔隙率下地表5%等濃度線距Y軸線距離與時(shí)間的關(guān)系曲線如圖12所示。當(dāng)土壤孔隙率為0.35時(shí)，對(duì)應(yīng)擴(kuò)散到地表的圓形區(qū)域半徑最大為1.69 m，當(dāng)土壤孔隙率為0.65時(shí)，對(duì)應(yīng)擴(kuò)散到地表的圓形區(qū)域半徑最大為1.975 m，所以以5%為危險(xiǎn)濃度可以將已經(jīng)確定泄漏位置2 m地表范圍定為危險(xiǎn)區(qū)域。

4 結(jié)論與建議

基于中國(guó)某油田超臨界CO2埋地管道的輸送工況，利用FLUENT軟件模擬預(yù)測(cè)了超臨界CO2管道發(fā)生小孔泄漏后在不同孔隙率的土壤中的擴(kuò)散特性，得出以下結(jié)論。

(1)土壤松散程度越高即孔隙率越大，管道泄漏的CO2擴(kuò)散速率越快，故而越松散的土壤在相同時(shí)間內(nèi)形成的危險(xiǎn)區(qū)域的面積就越大。

(2)超臨界管道為高壓輸送管道，CO2會(huì)在泄漏初期數(shù)秒內(nèi)在土壤孔隙中以射流擴(kuò)散，多孔介質(zhì)結(jié)構(gòu)會(huì)改變射流的方向并降低射流動(dòng)量，在泄漏口處由于壓力分層而出現(xiàn)負(fù)壓區(qū)域，但壓力會(huì)在極短的時(shí)間內(nèi)接近大氣壓，擴(kuò)散形式變?yōu)榫鶆驍U(kuò)散。故可以將均勻擴(kuò)散視為為超臨界CO2埋地管道發(fā)生小孔泄漏后擴(kuò)散的主要形式。

(3)以計(jì)算模擬得到的擴(kuò)散規(guī)律為依據(jù)，將5%確定為CO2危險(xiǎn)臨界濃度，泄漏時(shí)間大于800 s時(shí)地表CO2濃度分布基本趨于穩(wěn)定，此時(shí)可確定發(fā)生泄漏后地表的危險(xiǎn)區(qū)域。以最大土壤孔隙率0.65為例，模擬可知其對(duì)應(yīng)的危險(xiǎn)區(qū)域半徑為2 m。