萬庭輝，張可倪，李占釗，王靜麗，于彥江

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

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

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

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

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

1 軟件介紹及研究方法

1.1 軟件簡介

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

1.2 研究方法

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

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

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

1.2.1 控制方程

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

ρgVg)+qg

(1)

(2)

(3)

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

(4)

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

Sw+Ss+Sg=1

(5)

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

1.2.2 方程離散

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

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

(6)

(7)

(8)

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

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

(9)

(10)

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

(11)

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

表1 主要變量的選擇

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

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

(12)

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

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

(13)

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

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

(14)

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

本程序采用這樣一個數(shù)值算法去建立式(13)的雅可比矩陣：對于全隱式項來說，用數(shù)值微分法去求解雅可比矩陣。對于隱式壓力顯式飽和法來說，用更簡單的半解析方法去求解雅可比矩陣。形成的線性方程可通過迭代法去求解。

1.2.3 賓漢流體流變模型

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

(15)

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

(16)

(17)

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

(18)

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

1.2.4 流體性質(zhì)計算

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

2 軟件應(yīng)用

2.1 模型建立

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

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

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

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

2.2 模擬結(jié)果

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

圖3 壓力場圖

圖4 砂相飽和度場圖

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

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

2.3 效果分析

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

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

3 結(jié)論

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