天然氣水合物非成巖出砂儲(chǔ)層砂流數(shù)值模擬

萬(wàn)庭輝，張可倪，李占釗，王靜麗，于彥江

(1.中國(guó)地質(zhì)調(diào)查局廣州海洋地質(zhì)調(diào)查局，廣州 510760；2.東北石油大學(xué)石油工程學(xué)院，大慶 163318；3.暨南大學(xué)地下水與地球科學(xué)研究院，廣州 510632)

天然氣水合物分布廣、埋藏淺、清潔無(wú)污染、儲(chǔ)量巨大，被視為油氣領(lǐng)域最具潛力的清潔能源。蘇聯(lián)、加拿大、美國(guó)、日本、中國(guó)等國(guó)家已先后在凍土層或海域成功進(jìn)行了天然氣水合物試驗(yàn)性開(kāi)采。2017年5月18日，由中國(guó)地質(zhì)調(diào)查局組織實(shí)施的南海北部神狐海域水合物試采成功，實(shí)現(xiàn)了持續(xù)產(chǎn)氣60 d和產(chǎn)氣總量超過(guò)30萬(wàn)m3的世界紀(jì)錄。2020年3月，由中國(guó)地質(zhì)調(diào)查局組織實(shí)施的中國(guó)海域天然氣水合物第二輪試采取得成功，并超額完成目標(biāo)任務(wù)。本輪試采1個(gè)月產(chǎn)氣總量86.14萬(wàn)m3、日均產(chǎn)氣量2.87萬(wàn)m3，是第一輪60 d產(chǎn)氣總量的2.8倍。攻克了深海淺軟地層水平井鉆采核心技術(shù)，實(shí)現(xiàn)了從探索性試采向試驗(yàn)性試采的重大跨越[1-12]。

天然水合物開(kāi)采過(guò)程中，水合物相變引起的儲(chǔ)層應(yīng)力行為的改變，可能會(huì)導(dǎo)致一系列地質(zhì)風(fēng)險(xiǎn)問(wèn)題，如井壁失穩(wěn)、井眼出砂、地層垮塌以及海床滑坡等。蘇聯(lián)麥索亞哈陸域水合物田(1967年)[13-14]、加拿大麥肯齊陸域水合物試采(2007年和2008年)[15-17]、阿拉斯加北坡陸域水合物試采(2008年和2012年)[17-20]和日本南海海槽海域水合物試采(2013年和2017年)[18-20]，這些試采工程均以降壓開(kāi)采為主，大部分地區(qū)的開(kāi)采過(guò)程中均有不同程度的出砂情況[20-22]；尤其是2013年和2017年日本南海海槽進(jìn)行的兩次海域天然氣水合物試采均因井筒大量出砂而被迫中止，這表明出砂問(wèn)題嚴(yán)重影響天然氣水合物開(kāi)采過(guò)程中的安全和效率，全面了解“砂”的行為，以期實(shí)現(xiàn)有效控砂，是水合物安全高效、長(zhǎng)期可控開(kāi)采的關(guān)鍵。

數(shù)值模擬是了解天然氣水合物開(kāi)采儲(chǔ)層出砂的有效方法。在模型研究方面，Ranjith等[23]綜述了2014年以前開(kāi)展的與出砂問(wèn)題相關(guān)的模型研究。認(rèn)為出砂的基本機(jī)理包括剪切和拉伸破壞、臨界壓力梯度、壓力下降和侵蝕。目前為止，連續(xù)體和離散元法都已應(yīng)用在出砂模型中。其中連續(xù)體法包括基于塑性理論、混合理論及微觀力學(xué)建立起來(lái)的方法。大部分方法都考慮了流體力學(xué)或剪切破壞機(jī)制，有些將它們結(jié)合起來(lái)。Uchida等[24-25]提出了一個(gè)水合物儲(chǔ)層出砂模型，綜合描述了水合物的出砂特征，包括顆粒的分離、運(yùn)移、沉積變形和水合物的離解。這些過(guò)程是傳熱-流動(dòng)-力學(xué)的耦合過(guò)程。顆粒脫離引起應(yīng)力降低，沉積物剪切變形引起顆粒脫離，顆粒流動(dòng)改變多相流體壓力和溫度分布。研究發(fā)現(xiàn)，在不同的操作方法中，減小降壓率是降低出砂量的最有效方法。Yan等[26]利用多場(chǎng)耦合理論系統(tǒng)分析了影響井周?chē)纸夥植?、塑性變形區(qū)及出砂規(guī)律的影響因素，結(jié)合TOUGH+HYDRATE和ABAQUS來(lái)分析出砂的多場(chǎng)耦合。Ning等[27]構(gòu)建了基于連續(xù)-離散介質(zhì)理論的水合物開(kāi)采儲(chǔ)層響應(yīng)與出砂數(shù)值模型，將連續(xù)介質(zhì)模型獲取的應(yīng)力、流速和水合物飽和度等數(shù)據(jù)作為邊界和初始條件傳遞給離散元模型，進(jìn)而分析出砂機(jī)理、規(guī)律并預(yù)測(cè)出砂量。

目前中外水合物開(kāi)采出砂數(shù)值模擬研究大多是建立溫度-水力-應(yīng)力-化學(xué)耦合出砂模擬模型，使用TOUGH+HYDRATE水合物模擬器獲取開(kāi)采時(shí)流場(chǎng)如壓力、飽和度、產(chǎn)氣產(chǎn)水速率等參數(shù)變化以及溫度場(chǎng)如溫度和熱量變化規(guī)律，然后導(dǎo)入應(yīng)力場(chǎng)模擬軟件FLAC3D，獲取地層中的應(yīng)力和位移演化情況，結(jié)合顆粒流軟件PFC3D、工程模擬有限元軟件ABAQUS等[28-39]，對(duì)水合物儲(chǔ)層出砂機(jī)制和出砂規(guī)律進(jìn)行研究。此類(lèi)研究方法考慮影響因素全面，雖然有通過(guò)力學(xué)耦合描述“砂”的行為，但都沒(méi)有涉及“砂”的流動(dòng)過(guò)程。針對(duì)泥質(zhì)粉砂型天然氣水合物儲(chǔ)層，水合物分解后儲(chǔ)層多孔介質(zhì)中“砂”流動(dòng)過(guò)程的數(shù)值模擬，現(xiàn)通過(guò)通用儲(chǔ)層滲流模擬器GPSFLOW，實(shí)現(xiàn)水、氣(甲烷)、砂三相三組分模擬，把砂當(dāng)作一種非牛頓流體，建立賓漢流體型砂流模型，把“砂”當(dāng)成賓漢流體進(jìn)行模擬，并用理想模型進(jìn)行驗(yàn)證。

1 軟件介紹及研究方法

1.1 軟件簡(jiǎn)介

GPSFLOW(general-purpose subsurface flow simulator)軟件是一套考慮等溫復(fù)雜介質(zhì)系統(tǒng)中油、氣、水三相的牛頓和非牛頓流體多相流動(dòng)的商業(yè)數(shù)值模擬程序，允許用戶設(shè)定液相組分為牛頓流體或非牛頓流體。GPSFLOW假設(shè)非牛頓流體流動(dòng)符合冪律、賓漢流體等流變模型，黏度隨著剪切速率和孔隙流速的變化而變化。程序可應(yīng)用于一維、二維或三維的、各向異性的、非均質(zhì)孔隙裂縫系統(tǒng)。

1.2 研究方法

非牛頓流體的表觀黏度通常取決于孔隙流速或壓力梯度。不同非牛頓流體在孔隙介質(zhì)中流動(dòng)的流變函數(shù)不同。圖1給出了在孔隙介質(zhì)流動(dòng)中經(jīng)常遇到的非牛頓流體的剪切應(yīng)力和剪切速率間的典型關(guān)系。

τ為剪切應(yīng)力；γ為剪切速率；μ為流體動(dòng)力黏性系數(shù)(即黏度)；μb為賓漢流體的塑性黏度系數(shù)；μa為假塑性流體的塑性黏度系數(shù)；A、B、C、D分別為流體的擬合曲線

在天然氣水合物儲(chǔ)層中，水合物作為膠結(jié)物能有效增加儲(chǔ)層強(qiáng)度，水合物分解后，儲(chǔ)層膠結(jié)程度降低，地層應(yīng)力重新分布，當(dāng)有效壓力梯度達(dá)到臨界值時(shí)，孔徑較小的砂體顆粒被水合物分解產(chǎn)生的水氣攜帶并開(kāi)始流動(dòng)，表現(xiàn)為儲(chǔ)層出砂；當(dāng)有效壓力梯度不滿足條件時(shí)，微小砂粒度又會(huì)靜止沉淀，反過(guò)來(lái)影響儲(chǔ)層孔隙度、滲透率等，這一現(xiàn)象和非牛頓流體中賓漢流體的流動(dòng)特性非常相似。通過(guò)模擬器GPSFLOW，在軟件中將“非水”組分變?yōu)椤吧啊?，并且把“砂”?dāng)成賓漢流體模擬。在多相流達(dá)西定律的基礎(chǔ)上，軟件結(jié)合砂流的非牛頓流動(dòng)過(guò)程，砂、水與氣的混合流動(dòng)過(guò)程，描述水合物分解后儲(chǔ)層中“砂”的運(yùn)移過(guò)程。

1.2.1 控制方程

對(duì)于一個(gè)包含氣、水、砂三種物質(zhì)成分的三相流體等溫系統(tǒng)，任意流動(dòng)區(qū)域的三相物質(zhì)平衡方程如下。

ρgVg)+qg

(1)

(2)

(3)

β相(β=g為氣，β=w為水，β=s為砂)的流動(dòng)服從多相流達(dá)西定律，即

(4)

三種流體相在任意時(shí)刻的孔隙中均滿足如下約束條件：

Sw+Ss+Sg=1

(5)

砂作為特殊流體，假設(shè)其相對(duì)滲透率總是1.0，與水相和氣相之間不存在虹吸現(xiàn)象，即毛細(xì)壓力為0。水氣之間按常規(guī)的兩相流處理其毛細(xì)壓和相對(duì)滲透率的計(jì)算。

1.2.2 方程離散

GPSFLOW程序利用有限體積法對(duì)連續(xù)性方程進(jìn)行空間離散，用一系列積分有限差分方程表達(dá)油氣水的物質(zhì)平衡方程；利用全隱式的穩(wěn)定性和時(shí)間步長(zhǎng)較長(zhǎng)的特點(diǎn)，或利用自適應(yīng)隱式方法(AIM)可以加快模擬和降低空間存儲(chǔ)要求的特點(diǎn)去求解這些離散非線性方程，對(duì)于簡(jiǎn)單的問(wèn)題，也可以采用隱式壓力顯式飽和度(IMPES)方法，以提高計(jì)算速度；一般通過(guò)確定有限子域或網(wǎng)格塊的屬性來(lái)描述流體和巖石的熱力學(xué)性質(zhì)；通過(guò)有限差分逼近法去計(jì)算通過(guò)連通網(wǎng)格塊表面部分的質(zhì)量流量；應(yīng)用牛頓-拉甫森迭代程序求解這些離散的非線性的有限差分物質(zhì)平衡方程。

對(duì)上述質(zhì)量守恒方程用積分有限差分方法把連續(xù)性方程空間離散化，用向后差分方法實(shí)現(xiàn)該方程的時(shí)間離散化。寫(xiě)成殘差形式為

(6)

(7)

(8)

式中：i=1,2,…,N；m可以是時(shí)間點(diǎn)n或n+1，如果m=n，即單元i是應(yīng)用隱式壓力顯式飽和度法(IMPES)，如果m=n+1，即單元i是應(yīng)用全隱式方法；N為網(wǎng)格的總數(shù)量；n為前時(shí)間點(diǎn)；n+1為有待解決的現(xiàn)有時(shí)間點(diǎn)；Vi為成分i的體積;R為組分在當(dāng)前時(shí)間點(diǎn)質(zhì)量守恒計(jì)算的殘余量；λ為β相流度，定義為流體的相對(duì)滲透率與黏度之比；γ為導(dǎo)流系數(shù)；ψ為網(wǎng)格的總相壓力；Q為源匯項(xiàng)；Δt為時(shí)間步長(zhǎng);ηi包含臨近塊(j)組或單元i的結(jié)點(diǎn)，這樣直接包含了單元i；下腳標(biāo)ij+1/2為單元i和j分界面的一個(gè)恰當(dāng)平均值，與β相流度(相對(duì)滲透率與黏度之比)有關(guān)。

流動(dòng)邊界的傳導(dǎo)系數(shù)定義為

(9)

(10)

式中：Aij為連通單元i和j的共同界面面積；di為從單元i的中心到單元i和j交界面的距離；dj為從單元j的中心到單元i和j交界面的距離；kij+1/2為沿著單元i和j連通處的平均絕對(duì)滲透率；Di為單元i中心的深度；單元i的匯點(diǎn)/源點(diǎn)項(xiàng)定義為

(11)

應(yīng)用牛頓-拉甫森迭代法去求解一個(gè)流動(dòng)系統(tǒng)的方程，每個(gè)單元都有油氣水三成分的物質(zhì)平衡方程，表現(xiàn)為3×N個(gè)耦合非線性方程。每個(gè)單元選擇3個(gè)主要變量(x1、x2、x3)，分別是水壓力、含砂飽和度和含氣飽和度，如表1所示。

表1 主要變量的選擇

根據(jù)這3個(gè)主要變量，利用牛頓-拉甫森迭代方法得

[xk,p+1-xk,p]=0

(12)

式(12)中：指數(shù)k=1，2，3分別代表主要變量x1、x2、x3的下角標(biāo)1，2，3；p為迭代點(diǎn)。式(12)可表示為

k=1,2,3，i=1,2,…,N

(13)

迭代中產(chǎn)生的主要變量的增量表達(dá)式為

δxk,p+1=xk,p+1-xk,p

(14)

式(13)代表3×N個(gè)未知δxk,p+1的一系列3×N個(gè)線性方程。

本程序采用這樣一個(gè)數(shù)值算法去建立式(13)的雅可比矩陣：對(duì)于全隱式項(xiàng)來(lái)說(shuō)，用數(shù)值微分法去求解雅可比矩陣。對(duì)于隱式壓力顯式飽和法來(lái)說(shuō)，用更簡(jiǎn)單的半解析方法去求解雅可比矩陣。形成的線性方程可通過(guò)迭代法去求解。

1.2.3 賓漢流體流變模型

賓漢流體是一種黏塑性材料，在低應(yīng)力下表現(xiàn)為剛性，在高應(yīng)力下像黏性流體一樣流動(dòng)。砂流將按非牛頓流體類(lèi)型中的賓漢流體進(jìn)行處理。由于其切應(yīng)力與剪切速率呈線性關(guān)系(圖1)，因此在高應(yīng)力下賓漢流體具有牛頓流體性質(zhì)[40-41]。

(15)

式(15)中：τ為剪應(yīng)力，Pa；dv/dy為剪切速率，s-1；η為運(yùn)動(dòng)黏性系數(shù)，Pa·s；τ0為屈服應(yīng)力或屈服值，Pa。在數(shù)值模擬過(guò)程中，對(duì)賓漢流體的有效壓力梯度進(jìn)行修正比描述其表觀黏度有更高的運(yùn)算效率[42]，賓漢流體服從的達(dá)西定律為

(16)

(17)

式中：μb為賓漢流體的塑性黏度系數(shù)，與砂質(zhì)量分?jǐn)?shù)相關(guān)，Pa·s；?Φe為有效壓力梯度，Pa/m；k為絕對(duì)滲透率，D；kr為相對(duì)滲透率；v為達(dá)西流速，m/s；G為使得賓漢流體流動(dòng)的最小壓力梯度，Pa/m[43]，滿足：

(18)

式(18)中：α為經(jīng)驗(yàn)系數(shù)，或者擬合參數(shù)。

1.2.4 流體性質(zhì)計(jì)算

模型中儲(chǔ)層流動(dòng)系統(tǒng)由水、氣和砂3個(gè)組分組成。流體性質(zhì)計(jì)算主要包括流體的密度、黏度、相對(duì)滲透率和毛細(xì)壓等，這些參數(shù)一般是溫度與壓力的函數(shù)，氣體(甲烷)參數(shù)可通過(guò)PVT表輸入進(jìn)行插值，模型暫時(shí)只考慮等溫情況。砂的黏度和密度被認(rèn)為是常數(shù)，不受溫壓變化影響，相對(duì)滲透率總是1，毛細(xì)壓為0。

2 軟件應(yīng)用

2.1 模型建立

建立均值理想模型，模擬開(kāi)采區(qū)域長(zhǎng)寬均為1 000 m，厚10 m；模型四周為恒壓邊界，開(kāi)采井井筒半徑為0.15 m，位于模型中心處。均值儲(chǔ)層原始滲透率k=2 mD，孔隙度Φ=0.41，儲(chǔ)層初始?jí)毫?5 MPa，儲(chǔ)層初始溫度14.1 ℃；假定水合物分解后儲(chǔ)層多孔介質(zhì)中砂飽和度為Ss=0.1，氣飽和度為Sg=0.2，其余為水。砂流啟動(dòng)壓力梯度為1 MPa/m，其中殘余砂飽和度為0.067，理想模型儲(chǔ)層特征參數(shù)如表2所示。

表2 理想模型儲(chǔ)層特征參數(shù)

理想模型xy平面網(wǎng)格數(shù)量為3 293個(gè)(圖2所示)，z方向上分為10個(gè)網(wǎng)格。整個(gè)模擬過(guò)程所要求解的方程數(shù)為32 930個(gè)，近井筒周?chē)W(wǎng)格進(jìn)行了加密處理。

圖2 網(wǎng)格劃分示意圖

2.2 模擬結(jié)果

上述理想模型以井底流壓Pwf=0.5 MPa開(kāi)采10 d后，隨即以井底流壓Pwf=2 MPa開(kāi)采50 d，該模型開(kāi)采10、30、50、60 d的壓力和砂飽和度場(chǎng)圖分別如圖3和圖4所示。

圖3 壓力場(chǎng)圖

圖4 砂相飽和度場(chǎng)圖

模擬開(kāi)采期間，砂的流動(dòng)呈現(xiàn)賓漢流體的流動(dòng)特征，砂的產(chǎn)出速率和累積產(chǎn)砂量如圖5所示。開(kāi)采初期(0～10 d)，壓力迅速傳遞，生產(chǎn)壓差小于啟動(dòng)壓力梯度，砂相不能流動(dòng)，產(chǎn)砂速率為0 m3/d，累積產(chǎn)砂0 m3；開(kāi)采后期(10～60 d)，生產(chǎn)壓差大于啟動(dòng)壓力梯度，砂相持續(xù)流動(dòng)，由于井周砂相飽和度不斷增加，井周滲透率不斷減小，產(chǎn)砂速率由最初的0.6 m3/d，逐漸減小至0.000 01 m3/d，累積產(chǎn)砂1.62 m3。

圖5 產(chǎn)砂速率及累積產(chǎn)砂

2.3 效果分析

為探究泥質(zhì)粉砂型天然氣水合物儲(chǔ)層，水合物分解后儲(chǔ)層多孔介質(zhì)的宏觀砂流運(yùn)移規(guī)律，結(jié)合通用儲(chǔ)層滲流模擬器GPSFLOW，創(chuàng)新性地提出把“砂”當(dāng)作非牛頓流體，通過(guò)建立賓漢流體型砂流模型，實(shí)現(xiàn)水、氣(甲烷)、砂三相三組分模擬，并用理想模型進(jìn)行了驗(yàn)證。

模擬結(jié)果表明該砂流模擬方法簡(jiǎn)單高效，并可提供試采場(chǎng)地規(guī)模、長(zhǎng)時(shí)間尺度的砂流運(yùn)移規(guī)律可視化顯示；后續(xù)通過(guò)耦合水合物開(kāi)采滲流模型，結(jié)合試采場(chǎng)地天然氣水合物儲(chǔ)層層間非均質(zhì)地質(zhì)模型、各小層砂流啟動(dòng)壓力梯度、降壓條件等實(shí)際生產(chǎn)數(shù)據(jù)，未來(lái)有望實(shí)現(xiàn)對(duì)水合物開(kāi)采過(guò)程中場(chǎng)地規(guī)模砂流運(yùn)移規(guī)律的動(dòng)態(tài)監(jiān)測(cè)，有力支撐天然氣水合物試采工程現(xiàn)場(chǎng)決策。

3 結(jié)論

探討了“砂”在非成巖儲(chǔ)層的流動(dòng)過(guò)程，使用通用儲(chǔ)層滲流模擬器GPSFLOW實(shí)現(xiàn)賓漢型砂流的數(shù)值模擬。模擬器充分考慮了“砂”的非牛頓流動(dòng)過(guò)程，可實(shí)現(xiàn)三維非結(jié)構(gòu)網(wǎng)格、非均質(zhì)孔隙/裂隙介質(zhì)中的“砂”運(yùn)移情況的模擬，通過(guò)隱式積分有限差分求解模型，具有良好的穩(wěn)定性和收斂性。提出的砂流模擬方法簡(jiǎn)單高效，后續(xù)通過(guò)耦合水合物開(kāi)采滲流模型，有望實(shí)現(xiàn)對(duì)水合物開(kāi)采過(guò)程中砂流的準(zhǔn)確模擬。通過(guò)設(shè)置不同巖性、砂的性質(zhì)、降壓條件等，可以分析水合物開(kāi)采潛力、砂的運(yùn)移規(guī)律、地層壓力分布，預(yù)測(cè)地層出砂情況和地層穩(wěn)定性，為實(shí)際試采工程提供技術(shù)支撐，預(yù)測(cè)出砂情況和評(píng)估潛在風(fēng)險(xiǎn)。