王笑雨，張可霓，郭朝斌，陳開軍，李 毅

(1.北京師范大學(xué)水科學(xué)研究院，北京 100875;2.北京師范大學(xué)地下水污染控制與修復(fù)教育部工程研究中心，北京 100875;3.同濟(jì)大學(xué)機(jī)械與能源工程學(xué)院，上海 200092;4.杰瑞能源服務(wù)有限公司，山東煙臺 264003)

巖屑回注技術(shù)將鉆屑研磨成細(xì)小顆粒與水或其他液體以及稠化劑制漿并注入到地層中［1～2］，從而實現(xiàn)鉆井零排放，并有效規(guī)避運(yùn)輸途中的泄露風(fēng)險。相比于以往的鉆井廢水無害化處理［3～4］方式，該技術(shù)經(jīng)濟(jì)環(huán)保，世界許多地區(qū)均采用該技術(shù)處理鉆井廢棄物，如墨西哥灣、阿拉斯加等［5］。由于地質(zhì)信息的局限性，實際回注工程仍存在不確定性及風(fēng)險，例如:漿液竄流對地下環(huán)境造成污染;攜帶能力小的漿液會導(dǎo)致井筒或者孔隙堵塞［6］;注入率過低時產(chǎn)生沉淀床堵塞地層，過高時導(dǎo)致地層壓力急劇提升;鉆井回注能力難以確定等。通過數(shù)值模擬技術(shù)評估上述不確定性，制定合理工程計劃，可有效規(guī)避巖屑回注實際工程中的風(fēng)險。

本文針對水基鉆井液巖屑回注工程，建立賓漢流體型漿液的巖屑回注模型。目前，賓漢流體的數(shù)值模擬方式主要有:多相態(tài)賓漢流體的數(shù)值模擬［7］;分形樹狀網(wǎng)格孔隙介質(zhì)中的賓漢流體模擬［8］;離散網(wǎng)格玻爾茲曼方法模擬［9］;賓漢流體混合表述模型［10］;基于IPSAR算法的直管內(nèi)含顆粒賓漢流體兩相湍流模型［11］;圓管中牛頓流體與賓漢流體分層流動模型［12］等。由于賓漢流體的多樣性以及地下水流動系統(tǒng)的復(fù)雜性，上述模型難以直接應(yīng)用到巖屑回注技術(shù)上。

本文在TOUGH2［13］基礎(chǔ)上建立三維全隱式積分有限差分模型，并通過Fortran編程實現(xiàn)三維非結(jié)構(gòu)網(wǎng)格、非均質(zhì)孔隙/裂隙介質(zhì)中的漿液運(yùn)移情況模擬。針對賓漢流體型漿液，考慮其非牛頓流動過程，及其與地下水流動系統(tǒng)的相互作用，包括漿液與水混合流動及沉淀溶解過程，沉淀溶解對地層巖性的影響。通過理想模型展示該軟件的部分模擬功能，探討漿液在地層中遷移規(guī)律以及間歇回注對達(dá)到地層破裂壓力時間的影響。

1 軟件介紹及研究方法

1.1 軟件功能

本文開發(fā)的軟件在TOUGH2基礎(chǔ)上進(jìn)行改編，建立孔隙/裂隙介質(zhì)中巖屑回注的三維全隱式積分有限差分模型，并通過流體狀態(tài)方程實現(xiàn)賓漢流體型漿液的運(yùn)移模擬。

通過積分有限差分法進(jìn)行空間離散，軟件可以更精確地描述復(fù)雜幾何外形的地層，從而處理三維非結(jié)構(gòu)網(wǎng)格，并可以實現(xiàn)多重孔隙度模型，從而模擬裂隙介質(zhì)中的滲流過程。在時間上采取全隱式一階向后差分，保證求解方程的收斂性和穩(wěn)定性。在多相流達(dá)西定律的基礎(chǔ)上，軟件結(jié)合漿液的非牛頓流動過程、漿液與水的混合流動過程、漿液與水混合流體的沉淀與稀釋過程，描述巖屑回注漿液在地下的運(yùn)移過程。此外，軟件還充分考慮了沉淀對地層孔隙度、滲透率的影響，以及滲透率改變對漿液流動的影響。從而實現(xiàn)對三維非結(jié)構(gòu)網(wǎng)格、非均質(zhì)孔隙/裂隙介質(zhì)中的漿液運(yùn)移情況進(jìn)行模擬。

通過設(shè)置不同巖性、漿液性質(zhì)、注入條件等，可以分析漿液運(yùn)移規(guī)律、地層壓力分布、蓋層封閉性，預(yù)測各回注層位注入量以及地層壽命，確定最優(yōu)回注方案，為實際巖屑回注工程提供技術(shù)支持，預(yù)測回注能力以及評估潛在風(fēng)險。

1.2 研究方法

巖屑回注過程中的地下水流動系統(tǒng)由固液兩相組成。液相由漿液和水兩組分組成，其鹽度由漿液的質(zhì)量分?jǐn)?shù)描述，固相為漿液釋水沉淀物。該方案假設(shè)漿液和水可以混溶。回注早期，漿液溶解于地下水流動系統(tǒng)的液相中，達(dá)到一定濃度后視液相為賓漢流體。當(dāng)賓漢流體的壓力梯度小于臨界值時，流體停止流動，并在數(shù)小時后逐漸沉淀;當(dāng)壓力梯度足夠大時，沉淀物被攪動起來再次以漿液形式融入地下流體流動。沉淀物的形成會改變地層孔隙結(jié)構(gòu)，減小孔隙度以及滲透率，并進(jìn)一步影響流體流動。

1.2.1 控制方程

地下多相流體流動系統(tǒng)由固液兩相、水和漿液兩組分組成。漿液組分在固液兩相中通過釋水沉淀和吸水溶解轉(zhuǎn)換。在等溫系統(tǒng)中，該系統(tǒng)可由兩個質(zhì)量守恒方程描述:

式中:Vn——模型中一個任意形狀的子區(qū)域n的體積/m3;

Гn——子區(qū)域n的邊界;

M——質(zhì)量累積量/kg;

κ——組分標(biāo)識，水或漿液;

F——質(zhì)量通量/(kg·s-1·m-2);

q——源匯項/(kg·s-1);

n——曲面微元dГn的內(nèi)法向量。

組分κ的質(zhì)量累積量Mκ包括各相態(tài)內(nèi)組分κ的質(zhì)量:

式中:Φ——孔隙度;

β——相態(tài)標(biāo)識，溶液相或固相;

Sβ——β 相態(tài)的飽和度;

ρβ——β 相密度/(kg·m-3);

Xκβ——組分κ在β相內(nèi)的質(zhì)量分?jǐn)?shù)。

組分κ的通量Fκ包括各相態(tài)內(nèi)的通量:

根據(jù)多相流達(dá)西定律:

式中:k——絕對滲透率/Darcy;

krβ——β 相相對滲透率;

μβ——黏度/(Pa·s);

Pβ——β 相壓力/Pa;

g——重力加速度/(m·s-2)。

1.2.2 空間及時間離散方案

為求解控制方程，軟件采用積分有限差分的空間離散方案［14］(圖1)，從而實現(xiàn)對非結(jié)構(gòu)網(wǎng)格的數(shù)值模擬。任意網(wǎng)格內(nèi)流體質(zhì)量可以用其平均質(zhì)量Mn表示:

任意網(wǎng)格表面通量積分可以近似描述成該網(wǎng)格各個表面上的通量和:

式中:Anm——網(wǎng)格n和m接觸面面積/m2;

Fnm——通過n和m接觸面的通量F的平均量/(kg·s-1·m-2)。

用網(wǎng)格Vn和Vm上的平均參數(shù)表示離散化的流量，則達(dá)西定律為:

式中:nm——通過插值計算的網(wǎng)格 n和 m的平均值量;

Dnm——網(wǎng)格n和m的結(jié)點距離/m。

將式(5)和(6)帶入式(1)，得到一階常微分方程:

對其進(jìn)行采用隱式一階向后差分的時間離散方案以保證方程的收斂性和穩(wěn)定性，則對于任意網(wǎng)格Vn，在tk+1=tk+Δt時刻的控制方程為:

對于整個模型來說，需求解的非線性方程個數(shù)為組分個數(shù)×網(wǎng)格個數(shù)，可通過Newton/Raphson迭代法求解上述非線性方程組。

圖1 積分有限差分的空間離散方案Fig.1 Space discretization in the integral finite difference method

1.2.3 賓漢流體流動過程

賓漢流體是一種黏塑性材料，在低應(yīng)力下表現(xiàn)為剛性，在高應(yīng)力下像黏性流體一樣流動。在鉆井工程中和漿液處理中常被當(dāng)作泥漿流動模型。由于其切應(yīng)力與剪切速率呈線性關(guān)系，因此在高應(yīng)力下賓漢流體具有牛頓流體性質(zhì)。

式中:τ——剪應(yīng)力/Pa;

dv/dy——剪切速率/s-1;

η——運(yùn)動黏性系數(shù)/(Pa·s);

τ0——屈服應(yīng)力或屈服值/Pa。

在數(shù)值模擬過程中，對賓漢流體的有效壓力梯度進(jìn)行修正比描述其表觀黏度有更高的運(yùn)算效率［5］。賓漢流體服從的達(dá)西定律為:

式中:μb——賓漢流體的塑性黏度系數(shù)，與漿液質(zhì)量分?jǐn)?shù)相關(guān)/(Pa·s);

▽Φe——有效壓力梯度/(Pa·m-1);

k——絕對滲透率/Darcy;

kr——相對滲透率;

v——達(dá)西流速/(m·s-1);

G——使得賓漢流體流動的最小壓力梯度/(Pa·m-1)，滿足［15］:

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

流體屈服值主要受添加劑和介質(zhì)密度影響［16］。根據(jù)實驗記錄繪制屈服值隨黃原膠每桶密度變化圖(圖2)。屈服值YP變化可以采用線性模型:YP=b×ρXG+a，其中 ρXG為黃原膠密度。

圖2 屈服值隨黃原膠濃度(kg·m－3)變化圖Fig.2 Change in the yield value with the concentration of xanthan gum(kg·m －3)

1.2.4 混合流體性質(zhì)

模型中地下水流動系統(tǒng)的流體為液相，由漿液和水兩個組分組成。漿液和水混合流體中，認(rèn)為漿液體積和水的體積是加性的，并假設(shè)漿液與水的擴(kuò)張系數(shù)和壓縮系數(shù)相同，從而混合流體密度為:

式中:ρw——水的密度/(kg·m-3);

ρsly——漿液密度/(kg·m-3);

ρmix——混合流體密度/(kg·m-3);

Xsly——漿液質(zhì)量分?jǐn)?shù)。

給定漿液在P0、T0條件下的密度參考值，對于任意溫壓條件，通過將式(15)求解該條件下漿液密度，從而帶入式(14)即可計算混合流體密度。

混合流體黏度與漿液的質(zhì)量分?jǐn)?shù)相關(guān)［17］:

式中:μmix——混合流體黏度/(Pa·s);

μw——水的黏度/(Pa·s);

v1，v2，v3——黏度模型系數(shù)。

1.2.5 沉淀模型

漿液進(jìn)入地層中，與地下流體混合流動。當(dāng)流體壓力梯度小于臨界值時，流體停止流動并在數(shù)小時后逐漸沉淀;當(dāng)壓力梯度足夠大時，沉淀物被攪動起來再次以漿液形式融入地下流體流動。假設(shè)漿液內(nèi)懸浮物與水的體積是可加的，則漿液質(zhì)量分?jǐn)?shù)與沉淀物固相飽和度滿足如下關(guān)系:

若Δt時刻內(nèi)有ΔXsly的漿液沉淀，則tk+1=tk+Δt時刻時固相飽和度:

漿液質(zhì)量分?jǐn)?shù):

若Δt時刻內(nèi)有Δss沉淀物吸水溶解，則tk+1=tk+Δt時刻時固相飽和度:

漿液質(zhì)量分?jǐn)?shù):

式中:ss——固相飽和度;

ρs——漿液內(nèi)懸浮物密度/(kg·m-3);

ssly(liq.)，Xsly——液相漿液體積分?jǐn)?shù)和質(zhì)量分?jǐn)?shù);

ss(sly.)，Xs——漿液內(nèi)懸浮顆粒體積分?jǐn)?shù)和質(zhì)量分?jǐn)?shù)。

1.2.6 滲透率變化模型

沉淀物對地層孔隙度的影響較為直觀與簡單，然而孔隙度變化對滲透率的影響卻非常復(fù)雜。實驗發(fā)現(xiàn)孔隙度的輕微減小可能引起滲透率大幅下降［18］。當(dāng)孔吼被沉淀堵住時，部分孔隙無法與孔隙通道連接形成無效孔隙。因此，滲透率變化的影響因素不止是孔隙度大小，還有孔隙空間的幾何形狀與沉淀物在孔隙內(nèi)的分布。

1.2.6.1 管道系列模型［19］

考慮天然孔隙通道由半徑大小交替的管道空間組成(圖3)，根據(jù)該模型可知滲透率的相對變化為:

式中:k0——原滲透率/Darcy;

k——沉淀后地層滲透率/Darcy;

Ss——固體飽和度;

Г——孔隙主體占孔隙管道的長度比例(圖3);

Φr——滲透率降低為0時對應(yīng)的孔隙度減小比例。

圖3 孔隙管道縮放模型Fig.3 Model for converging-diverging pore channels

1.2.6.2 平行板模型［19］

用平行板代替管道系列模型中的管道，則得到類似的滲透率相對變化模型:

1.2.6.3 直管模型［19］

對于管道半徑相同的直管模型，Г=0，Φr=0，從而滲透率的相對變化模型簡化為:

1.2.6.4 模型對比

對比三種滲透率變化模型，取Г=0.8，Φr=0.8，繪制滲透率相對變化隨固相飽和度變化圖(圖4)。在同樣參數(shù)下，直管模型是最寬松的滲透率變化模型，由于沒有孔吼限制，其滲透率下降遠(yuǎn)小于其他兩種模型。平板模型對地層滲透率下降做最保守的估計，隨沉淀產(chǎn)生地層堵塞嚴(yán)重，在固相飽和度為0.2時徹底堵塞，地層滲透率降低為0。

圖4 滲透率變化模型對比圖Fig.4 Comparison of different permeability change models

2 軟件應(yīng)用

2.1 模型設(shè)置

本文建立理想模型，包含三個砂巖含水層并由泥巖蓋層分隔，其中中部砂巖層(-1925～-1935 m)為回注目標(biāo)層，模型四周及頂部底部邊界為無流量邊界。模型由半徑10 km的徑向網(wǎng)格組成，厚60 m，回注目標(biāo)層厚10 m，各泥巖隔水層厚15 m，上下各砂巖層厚10 m(圖5)，只在回注目標(biāo)層射孔。模型參數(shù)見表1。

圖5 模型示意圖Fig.5 Diagram of the model

表1 模型參數(shù)設(shè)置Table 1 The parameters of the model

2.2 基礎(chǔ)模型

在注入層段以2.0 kg/s定速注入漿液，漿液密度1250 kg/m3，塑性黏度21 MPa·s，稠化劑(黃原膠)濃度為1.170 kg/m3。采用直管滲透率變化模型，模擬連續(xù)注入30 d，隨后停注30 d的漿液運(yùn)移過程。模擬結(jié)果如圖6～8所示。

由圖6可知，巖屑回注的最大壓力積聚出現(xiàn)在注入層。回注初期壓力迅速上升后變緩，最大達(dá)到靜水壓1.33倍，停住時壓力迅速傳遞從而急劇下降至稍高于靜水壓，致使注入期間漿液持續(xù)流動沉淀較少(圖7前30天模擬結(jié)果)，而停注后產(chǎn)生大量沉淀(圖7模擬結(jié)果31～60 d)。由于蓋層滲透率很低，沉淀先在頂部蓋層和底層出現(xiàn)(圖8a);地層壓差從井周向邊界逐漸減小，因此在注入層漿液暈出現(xiàn)沉淀并積聚，該范圍即為漿液的影響范圍(圖8b);停注后壓力迅速傳遞、地層壓差變小，漿液流動減緩直至停止并逐漸沉淀。由于井周漿液含量高于遠(yuǎn)處漿液含量，停注后井周沉淀含量高于遠(yuǎn)處沉淀含量;在整個模擬期間漿液暈處壓差低于漿液流動的最小壓差，因此持續(xù)產(chǎn)生沉淀，沉淀含量較多(圖8c、d)。

圖6 注入層壓力變化圖Fig.6 Pressure of the injection layer

圖7 系統(tǒng)內(nèi)漿液及沉淀變化圖Fig.7 Mass of slurry and precipitate in the subsurface system

圖8 固相飽和度分布圖Fig.8 Distribution of solid saturation

2.3 間歇注入對地層壽命的影響

在基礎(chǔ)模型的方案上，考慮間歇注入方案(表2)對地層壓力的影響。

模擬結(jié)果見圖9～11。圖9中，停注時地層壓力迅速傳遞從而急劇下降，產(chǎn)生大量沉淀。地層孔隙被沉淀堵塞，滲透性下降，因此在隨后的注入過程中壓力提升加劇，最終到達(dá)地層破裂壓力。在相同注入速率條件下(2 kg/s)，連續(xù)注入383 d才會使最大壓力積聚的地層到達(dá)地層破裂壓力，而在間歇注入方案中，最長需要270 d(120 d周期方案)，最短只需要82 d(10 d周期方案)。而在相同注入量條件下(平均注入速率相同，為1 kg/s)，連續(xù)注入方案在10 a內(nèi)都不會到達(dá)地層破碎壓力。

圖9 注入點壓力變化圖(虛線為到達(dá)地層破裂壓力時間)Fig.9 Pressure of injection point(the time when the formation breaks down is marked by dotted line)

圖10 沉淀質(zhì)量變化圖Fig.10 Mass of precipitate in the subsurface system

沉淀主要產(chǎn)生于停注期，因此系統(tǒng)中沉淀量隨間歇注入周期階梯型增長(圖10)，漿液質(zhì)量隨間歇注入周期性變化(圖11)。由于長周期間歇注入對地層損害較小，滲透性降低較少，在相同注入及停注時間時(120 d)，長周期間歇注入沉淀量明顯少于短周期間歇注入沉淀量(圖10)，漿液含量則高于短周期間歇注入的漿液含量(圖11)。

圖11 漿液質(zhì)量變化圖Fig.11 Mass of slurry in the subsurface system

3 結(jié)論

本文詳細(xì)探討了漿液在地下的流動過程，開發(fā)模擬器實現(xiàn)賓漢型漿液的巖屑回注數(shù)值模擬技術(shù)。模擬器充分考慮了漿液的非牛頓流動過程、漿液與水的混合流動過程、漿液與水混合流體的沉淀過程及沉淀對地層孔隙度滲透率的影響。該模擬器可實現(xiàn)三維非結(jié)構(gòu)網(wǎng)格、非均質(zhì)孔隙/裂隙介質(zhì)中的漿液運(yùn)移情況的模擬，通過隱式積分有限差分求解模型，具有良好的穩(wěn)定性和收斂性。

通過理想模型可知，漿液在回注期產(chǎn)生少量沉淀，而在停注期迅速產(chǎn)生大量沉淀。因此在隨后的注入過程中地層壓力提升加劇，存在到達(dá)破裂壓力的風(fēng)險。短周期間歇注入的壓力提升與沉淀量都遠(yuǎn)高于長周期間歇注入方案與連續(xù)注入方案。無論在壓裂回注還是低于破裂壓力的巖屑回注工程中，均應(yīng)特別注意短周期間歇注入對地層的損害。

［1］ 胡小剛，康濤，柴占文，等.國外鉆井巖屑處理技術(shù)與國內(nèi)應(yīng)用研制分析［J］.石油機(jī)械，2009，37(9):159-161.［HU X G，KANG T，CHAI Z W，et al.The abroad technology of drilling cuttings disposal and the analysis of domestic application and development［J］.China Petroleum Machinery，2009，27(9):159-161.(in Chinese)］

［2］ Piper W，Harvey T.Waste management Drilling fluids processing hand book［M］.Burlington，MA:ASME Shale Shaker Committee，2005:367-406.

［3］ 朱敬濤，韓志勇，馮欣，等.黃土—微生物系統(tǒng)處理含油廢水的靜態(tài)模擬試驗研究［J］.水文地質(zhì)工程地質(zhì)，2011，38(4):121-124.［ZHU J T，HAN Z Y，F(xiàn)ENG X，et al.Static simulation experimental research on the treatment of oily wastewater by natural loess-microorganism system［J］. Hydrogeology ＆Engineering Geology，2011，38(4):121-124.(in Chinese)］

［4］ 馮欣，韓志勇，羅維剛，等.黃土對含油廢水的吸附作用研究［J］.水文地質(zhì)工程地質(zhì)，2010，27(6):121-125.［FENG X，HAN Z Y，LUO W G，et al.Sorption of oily wastewater on loess in Lanzhou Area［J］.Hydrogeology＆ Engineering Geology，2010 ，27(6):121-125.(in Chinese)］

［5］ Guo Q，Geenhan T.An overview of drill cuttings reinjection–Lessons learned and recommendations［C］.11th International Petroleum Environmental Conference Albuquerque.New Mexico:2004:1-10.

［6］ MiSwaco. Cuttings re-injection: Site-specific processes to meet the challenges，issues or limitations of any waste-disposal injection site［M/OL］. ［2003］［2014-11-13］http://www.slb.com/～/media/Files/miswaco/brochures/cuttings _ reinjection_ms06182.pdf.

［7］ Wu Y S， PruessK.A numericalmethod for simulating non-Newtonian fluid flow and displacement in porous media［J］.Advances in Water Resources，1998，21:351-362.

［8］ Wang S F，Yu B M.A fractal model for the starting pressure gradient for Bingham fluids in porous media embedded with fractal-like tree networks［J］.International Journal of Heat and Mass Transfer，2011，54:4491-4494.

［9］ Tang G H，Wang S B，Ye P X，et al.Bingham fluid simulation with the incompressible lattice Boltzmann model［J］. Journal of Non-Newtonian Fluid Mechanics，2011，166:145-151.

［10］ Aposporidis A，Haber E，Olshanskii M A，et al.A mixed formulation of the Bingham fluid flow problem:Analysis and numericalsolution［J］. Computer Methods in Applied Mechanics and Engineering，2011，200:2434-2446.

［11］ 亢力強(qiáng)，曾卓雄，姜培正.賓漢流體與顆粒間的密相兩相湍流模型［J］.西安交通大學(xué)學(xué)報，2002，36(7):693-696.［KANG L Q，ZENG Z X，JIANG P Z.Study on dense two-phase turbulent flow of Bingham fluid and particles［J］.Journal of Xi’an Jiaotong University，2002，36(7):693-696.(in Chinese)］

［12］ 賀成才.牛頓流體—賓漢流體的圓管分層層流的數(shù)值模擬［J］.天然氣與石油，2006，24(2):15-18.［HE C C.Numerical simulation of stratified laminar flow of Newton fluid and Bingham fluid in pipe［J］.Natural Gas and Oil，2006，24(2):15-18.(in Chinese)］

［13］ Pruess K.TOUGH2-A general purpose numerical simulator for multiphase fluid and heat flow［R］.Berkeley，CA:Lawrence Berkeley Laboratory，1991.

［14］ Narasimhan T N，Witherspoon P A.An integrated finite difference method for analyzing fluid flow in porous media［J］.Water Resour Res ，1976，12(1):57-64.

［15］ Pascal H.A theoretical analysis on stability of a moving interface in a porous medium for bingham displacing fluids and its application in oil displacement mechanism［J］.Can J Chem Eng，1986，64:375–379.

［16］ Measurement understanding yield value home formulations［EB/OL］.［2002-01］［2014-11-13］.http:// www. lubrizol. com/Home-Care/Documents/Technical-Data-Sheets/TDS - 244 -Measurement-Understanding-Yield-Value-Personal-Care-Formulations.pdf.

［17］ Pruess K，Oldenburg C，Mcridis G.TOUGH2 User’s Guide，Version 2.0［R］.Berkeley，CA:Lawrence Berkeley National Laboratory，1999.

［18］ Verma A K，Pruess K，Tsang C F，et al.A study of two-phase concurrent flow of steam and water in an unconsolidated porousmedium ［C］. Proc23rd National Heat Transfer Conference， Am.Denver，CO:Society of Mechanical Engineers，1985:135-143.

［19］ Verma A，Pruess K.Thermohydrological conditions and silica redistribution near high-level nuclear wastes emplaced in saturated geological formations［J］.of Geophys Res，1988，93(B2):1159-1173.