陳 兵，康慶華，肖紅亮

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

隨著社會科學(xué)技術(shù)的進步和人們生活水平的提高，環(huán)保問題受到廣泛的關(guān)注，這不僅僅是我國面臨的一個嚴峻問題，也是全球亟待解決的問題。在這樣的國際大環(huán)境下，碳捕獲、利用與封存(Carbon Capture,Utilization and Storage,CCUS)這一項新興的、具有大規(guī)模CO2減排潛力和環(huán)境效益的技術(shù)[1-7]應(yīng)運而生。CO2的管道輸送是CCUS系統(tǒng)的中間環(huán)節(jié)，承擔著將CO2運輸?shù)焦ぷ鼽c的任務(wù)，是實現(xiàn)CCUS大規(guī)模工程化應(yīng)用的關(guān)鍵技術(shù)環(huán)節(jié)[8]。根據(jù)國家職業(yè)安全衛(wèi)生研究所(NIOSH)提供的數(shù)據(jù)顯示：CO2濃度達40 000 ppm(1 ppm=10-6)時，將直接危及生命[8]。盡管CO2管道輸送泄漏時無燃燒風險，但仍需將泄漏點周圍的超高壓風險當作泄漏風險來考慮，尤其是在發(fā)生預(yù)料之外的管道破裂時[9-10]。

管道破裂時會造成重大壓力損失和管道流體的損失[11]，在壓力釋放期間，減壓會使CO2氣體較之原先的體積以很大的膨脹系數(shù)擴散至大氣中，而且減壓會導(dǎo)致其發(fā)生相變[12-14]，隨著部分高蒸汽壓CO2汽化為氣體，剩下的CO2則變?yōu)橐后w和固體顆粒。通常在減壓結(jié)束后這種膨脹力會驅(qū)使這些CO2固體顆粒以很高的速度彈出，這會給人員和設(shè)備造成嚴重的危害[15-16]。此外，減壓致使CO2及其雜質(zhì)[17](如H2S)從管道中釋放出去，在釋放位置處，如果不與空氣迅速混合，那么在管道破裂位置相對較近范圍內(nèi)的環(huán)境中則具有毒性，當管道內(nèi)所有毒性物質(zhì)在失壓后全部脫離出來時，若不采取適當?shù)拇胧?，這些物質(zhì)可能會污染并沉淀在所在的區(qū)域，給暴露于該環(huán)境中的人造成危害，且在合適的溫度和壓力下，也會形成水合物[18-19]。CO2泄漏后直接排放在基質(zhì)上會出現(xiàn)干冰“噴砂效應(yīng)”，在失壓后，沙子、土壤、和其他固體碎片與生成的干冰相結(jié)合形成研磨能力更大的顆粒，增強了對泄漏位置和管道工程的腐蝕效應(yīng)[20]。CO2輸送管道破裂時，CO2流進入大氣后會擴散，在氣流消散前有不同的泄漏階段和主導(dǎo)作用，所以研究泄漏口附近的壓力隨時間的變化規(guī)律、壓力分布以及泄漏階段、持續(xù)時間和濃度等級是非常重要的，可以為管道放空CO2流的相態(tài)變化、干冰或者水合物的形成、CO2云擴散和管道完整性研究提供一定的參考價值。為此，本文選取某油田CO2管道輸送工況作為研究對象進行了動態(tài)數(shù)值模擬，得出的結(jié)論可為相關(guān)研究提供一定的技術(shù)支撐。

1 某油田管道輸送氣源指標及成分

經(jīng)調(diào)研，某油田含雜質(zhì)CO2管道輸送氣源指標見表1，CO2氣流成分見表2。

表1 某油田含雜質(zhì)CO2管道輸送氣源指標

表2 某油田管道輸送CO2氣流成分

2 超臨界CO2輸送管道泄漏點計算流體動力學(xué)(CFD)模擬研究

在設(shè)計和運行CO2管道輸送時，CO2泄漏模擬是定量風險分析的一部分，而分析隨后的CO2擴散規(guī)律將有助于對這類風險的影響因素進行評估。為了更深一步地了解CO2的擴散機理，本文重點研究CO2泄漏和擴散過程的特性，包括近場的壓力場和濃度場等，即選用Fluent軟件中的組元輸運模型，將網(wǎng)格“外疏里密”，建立結(jié)構(gòu)化網(wǎng)格，并利用某油田CO2管道輸送工況作為研究管道泄漏時CO2氣體擴散特性的入口邊界條件，選用湍流模型中的k-ε模型和真實氣體狀態(tài)P-R方程進行求解和計算，得到相應(yīng)的CO2濃度擴散云圖，最后經(jīng)數(shù)據(jù)處理后得到模擬結(jié)果。

2.1 Fluent軟件簡介

Fluent是世界領(lǐng)先的計算流體動力學(xué)(Computational Fluid Dynamics,CFD)軟件，在流體建模中被廣泛應(yīng)用。該軟件針對各種復(fù)雜流動的物理現(xiàn)象，采用不同的離散格式和數(shù)值方法，以期在特定的流域內(nèi)使計算速度、穩(wěn)定性和精度等方面達到最佳組合，從而高效地解決各個領(lǐng)域的復(fù)雜流體流動計算問題[21]。Fluent軟件主要基于N-S方程和有限體積法，適用于復(fù)雜的流體流動計算，其組元輸運模型可用于計算氣體的擴散過程，且可操作性較好，因此本文選用該軟件對超臨界CO2輸送管道泄漏過程進行數(shù)值模擬。

2.2 數(shù)值建模

超臨界CO2輸送管道泄漏過程主要分為射流階段和擴散階段兩個階段，建立的數(shù)值模型見圖1。

2.3 基本控制方程和狀態(tài)方程

CO2在泄漏或放空過程中，其流體流動仍遵循動量、能量、質(zhì)量三大守恒定律，并由此推導(dǎo)出在直角坐標系下CO2流體適用的最基本的連續(xù)方程、運動方程和能量方程如下。

(1) 連續(xù)方程：

(1)

式中：ρ為氣體密度；u、v、ω分別為x、y、z軸的分量，該方程本身守恒。

(2) 運動方程：

(2)

式中：pij為定點壓強的張量表示法，xi、xj為坐標矢量的張量表示形式，i、j表示自由指標，分別取1、2。

一般對于氣體流動而言，質(zhì)量力fi可以忽略，即認為fi=0，再將牛頓黏性應(yīng)力公式代入，可得到N-S方程(納維斯-托克斯方程)：

(3)

(4)

式中：ui、uj為速度矢量的張量表示形式。

(3) 能量方程：

(5)

該式表明，單位流體微團能量的增量等于單位時間流體面力所做功加上它所獲取的能量。

(4) 狀態(tài)方程：

本文在計算CO2流體物性參數(shù)和CO2流體擴散過程時均選用P-R方程：

(6)

式中：ω為偏心因子,ω=-1;μ為動力黏性系數(shù),μ=2；a、b為物質(zhì)特有常數(shù)，與臨界參數(shù)有關(guān)，其關(guān)系式如下：

式中：Tc、Pc為臨界溫度和臨界壓力，R為通用氣體常數(shù)。

2.4 網(wǎng)格劃分

建立好滿足要求的幾何模型后，需要進行網(wǎng)格劃分。網(wǎng)格劃分是CFD 流體力學(xué)計算中最為關(guān)鍵的一個環(huán)節(jié)，網(wǎng)格質(zhì)量的優(yōu)劣對最終求解結(jié)果的準確度有直接的影響。若網(wǎng)格劃分過疏，會使計算結(jié)果誤差過大或者所得結(jié)果根本無法收斂；若網(wǎng)格劃分過密，盡管會提高計算結(jié)果的準確度，但計算精度增幅不大，反而會極大地延長計算時間，造成計算成本的增加。

某油田CO2管道輸送工況中CO2濃度為5 000 ppm的擴散區(qū)域很大，考慮到氣體在大氣中的擴散呈三維對稱，可將三維模型簡化成在主風向上垂直于地面的一個剖面，并使用GAMBIT軟件進行網(wǎng)格劃分(見圖2)。在進行網(wǎng)格劃分時，為了保證計算速度和計算精度，將建立結(jié)構(gòu)化網(wǎng)格，即選用四邊形網(wǎng)格劃分區(qū)域，由于泄漏口射流速度較大，因此整個區(qū)域網(wǎng)格由外層向內(nèi)層逐步細化，由外到內(nèi)共分為5層。左邊設(shè)置為風速進口(VELOCITY_INLET)，上邊和右邊設(shè)置為壓力出口(PRESSURE_OUTLET)，底邊除進口外設(shè)置為墻(WALL)。

圖2 網(wǎng)格劃分示意圖Fig.2 Grid generation diagram

2.5 基本設(shè)定

啟動非定常求解器，重力加速度為9.81 m/s2，方向為y軸方向，壓力單位為kPa，環(huán)境溫度為280.95℃，大氣壓為101.325 kPa，開啟能量方程，選用湍流模型中的k-ε模型，利用水力直徑和湍流強度，開啟組元輸運模型，組分選用某油田CO2管道輸送組分(但這里僅以CO2為研究對象，故只追蹤CO2)。

2.6 入口邊界條件

考慮到輸送管道中CO2流體為超臨界/密相狀態(tài)，泄漏后必然會發(fā)生相變，因此計算得到孔口處流量、溫度、壓力隨時間的變化規(guī)律，并用C語言編寫UDF導(dǎo)入Fluent軟件作為入口邊界條件。

本文使用ASPEN HYSYS V8.4創(chuàng)建Depressuring-Dynamics單元，以儲罐表示管道，閥門表示孔口，對某油田CO2管道輸送工況進行動態(tài)數(shù)值模擬。該工況為CO2流體超臨界/密相輸送，絕熱過程，輸送壓力為10 MPa，初始溫度為60℃，最優(yōu)管徑為φ273 mm×10.0 mm。模擬得到孔口處流量、壓力和溫度隨時間的變化曲線，詳見圖3至圖5。

圖3 孔口處流量隨時間的變化曲線Fig.3 Curve of flow at the pipe orifice with time

圖4 孔口處壓力隨時間的變化曲線Fig.4 Curve of pressure at the pipe orifice with time

圖5 孔口處溫度隨時間的變化曲線Fig.5 Curve of temperature at the pipe orifice with time

經(jīng)擬合，可得到孔口處流量(Q)、壓力(P)和溫度(T)與時間(t)的關(guān)系式如下：

Q=14.2+1 263.4e(-t/92.3);

P=156.4+9 895.0e(-t/119.7);

T=-93.3-144.6e(-t/307.3).

2.7 無關(guān)性檢測

由于模型所涉及的擴散區(qū)域過大，因此需要對模型無關(guān)性區(qū)域進行檢測，即對CO2濃度波動進行評價。初始場域定為200 m×1 800 m，場域過大會使網(wǎng)格劃分極為復(fù)雜，為了避免網(wǎng)格奇異、過大縱橫比，網(wǎng)格劃分基本遵循“外疏里密”的原則，泄漏口位置根據(jù)左側(cè)場域變化進行調(diào)整。經(jīng)模擬驗證，左側(cè)CO2擴散變化范圍為120 m左右，在擴散時間為1 000 s時，上方130 m處、右側(cè)750 m處CO2擴散濃度不再發(fā)生變化，因此模型擴散區(qū)域定為150 m×1 000 m，泄漏口離左側(cè)邊界為200 m。

2.8 模擬結(jié)果與分析

利用Fluent軟件的模擬計算中，由于未考慮干冰的形成，在離開泄漏口，溫度會很快接近大氣溫度，高壓CO2射流決定了其擴散形貌，所以本次模擬不考慮速度場、溫度場，主要分析濃度場和近場的壓力場[22]。

2.8.1 CO2濃度的變化

本文利用Fluent軟件模擬得到的某油田管道泄漏CO2濃度云圖，見圖6。

圖6 某油田管道泄漏CO2濃度云圖Fig.6 Cloud map of leaking CO2 concentration of an oil field

由圖6可以看出：

(1) 在泄漏最初階段，CO2濃度云圖并未發(fā)生傾斜，可見在此階段初始泄漏高壓起主導(dǎo)作用，對應(yīng)模型中的射流階段。

(2) 隨著時間的推移，CO2擴散高度的增加，動能不斷減小，CO2濃度云圖開始向右傾斜，即風速開始起主導(dǎo)作用，大氣湍流作用增強，云團變化變得復(fù)雜化，對應(yīng)模型中的擴散階段。

(3) CO2擴散以順風風向為主，在泄漏源停止泄漏之前(800 s之前)，沿順風方向越遠，CO2擴散濃度越低，泄漏源停止泄漏之后，CO2會繼續(xù)沿順風方向運移，其濃度基本保持在10-2數(shù)量級(比實際略大)。

2.8.2 壓力的變化

本文利用Tecplot軟件對Case和Data文件進行后處理，將泄漏口附近3個時刻的壓力云圖與Tecplot導(dǎo)出的壓力等值線圖、流線圖集于一圖，得到3個時刻泄漏口附近的壓力云圖，見圖7至圖9。

圖7 10 s時泄漏口附近的壓力云圖局部放大圖Fig.7 Partial enlargement of pressure cloud near the leak hole at 10 s

圖8 50 s時泄漏口附近的壓力云圖局部放大圖Fig.8 Partial enlargement of pressure cloud near the leak hole at 50 s

圖9 200 s時泄漏口附近的壓力云圖局部放大圖Fig.9 Partial enlargement of pressure cloud near the leak hole at 200 s

根據(jù)上述給出的3個時刻泄漏口附近的壓力云圖，現(xiàn)將泄漏口上方0.1 m處的壓力導(dǎo)出，得到局部壓力變化曲線，見圖10。

圖10 局部壓力變化曲線Fig.10 Curve of local pressure

由圖7至圖10可以看出：

(1) 高壓CO2流體釋放到大氣中會形成自由射流，在泄漏口處壓力很高，有明顯的壓力分層，并形成負壓區(qū)域，50 s時壓力達到最高，300 s后趨于穩(wěn)定，即泄漏口壓力趨于大氣壓。

(2) 泄漏口附近壓力場幾乎呈對稱分布，流線呈輻射狀分布，如同CO2濃度云圖在泄漏口附近并未發(fā)生傾斜，此階段初始泄漏高壓起主導(dǎo)作用，對應(yīng)模型中的射流階段。

(3) 隨著時間的推移，泄漏口上方0.1 m處的壓力分布變得復(fù)雜化，負壓區(qū)區(qū)域由大變小而后又變大，這是由于高壓CO2流體會在泄漏初期快速膨脹又迅速減小；當壓力快速上升進入減壓階段時，減壓會使CO2氣體較之原先的體積以很大的膨脹系數(shù)擴散至大氣中，具體膨脹程度由管道工況壓力與溫度協(xié)同作用決定，此階段對應(yīng)模型中的擴散階段。

(4) 減壓時使負壓區(qū)區(qū)域變大，時間持續(xù)較長，使得擴散出來的高壓CO2流體與空氣更易混合且復(fù)雜化。結(jié)合圖6可知，在300 s之內(nèi)，CO2濃度等級特別高(從釋放位置起，會持續(xù)一段距離)，減壓致使CO2及其毒性物質(zhì)(如H2S)從管道中全部脫離出來，使管道破裂位置一定范圍內(nèi)的環(huán)境具有較大的毒性，會對意外進入該范圍內(nèi)的人類或動物造成嚴重的危害。為了安全，應(yīng)在危險濃度區(qū)域邊界放置警告牌，即：在橙色區(qū)域附近應(yīng)放置職業(yè)濃度警告牌，當人或動物進入該范圍時可能會出現(xiàn)危險；在紅色區(qū)域附近應(yīng)放置危險濃度警告牌[23]，當人或動物進入該范圍時會立即出現(xiàn)危險。

3 結(jié) 論

(1) 本文利用Fluent軟件，選用組分輸運模型對含雜質(zhì)CO2管道輸送泄漏進行了數(shù)值模擬計算，得到了CO2濃度和壓力云圖。

(2) 對應(yīng)模型中的泄漏階段劃分，由CO2濃度云圖可以看出CO2管道輸送泄漏主要分為兩個階段：①射流階段,該階段輸送壓力起主導(dǎo)作用，范圍小，持續(xù)時間較短；②擴散階段,該階段風速、大氣湍流起主導(dǎo)作用，范圍大，持續(xù)時間較長。

(3) 在泄漏口附近會形成壓力膨脹層，壓力值迅速上升又下降，并由高速氣流造成負壓區(qū)，該區(qū)域的流線和壓力場呈對稱分布。

(4) 通過數(shù)值模擬得到的含雜質(zhì)CO2管道輸送泄漏擴散規(guī)律基本符合實際，但由于Fluent軟件的可操作性較大，模擬得到的具體數(shù)值與實際略有差別，仍需參照實驗值并設(shè)定具體函數(shù)進一步研究CO2管道輸送泄漏過程特性。

(5) 為了安全，根據(jù)CO2擴散特性和壓力響應(yīng)應(yīng)在泄漏點一定范圍內(nèi)(警示區(qū)域、危險區(qū)域)放置安全警示牌，即：當人或動物進入該范圍時可能會出現(xiàn)危險時放置職業(yè)濃度警告牌；當人或動物進入該范圍時會立即出現(xiàn)危險時放置危險濃度警告牌。在控制好危險范圍的同時，更要做好氣源地的除雜工作，嚴格限制氣源中H2S、CO、SO2等雜質(zhì)的含量。