溶蝕白云巖含水層注漿漿水兩相運移特征LBM模擬研究

周禹良 袁東鋒 楊 雪 李 斌

(1.華北科技學院安全工程學院,河北 廊坊 065201;2.北京中煤礦山工程有限公司,北京 100013;3.礦山深井建設(shè)技術(shù)國家工程研究中心,北京 100013)

礦井水害具有強烈的致災(zāi)性,若發(fā)生涌水量過大或突水淹井事故,將造成巨大的經(jīng)濟損失和惡劣的社會影響[1-2]。注漿技術(shù)是一種主動防治礦井水害的有效手段[3-4],是通過一定壓力將漿液注入地下含水層,驅(qū)排地層中的地下水,漿液凝固后充填封堵地層導水通道形成注漿帷幕[5-6],提高巖體的黏聚力和內(nèi)摩擦角[7-8],以達到堵水和加固地層的目的。巖體注漿是典型的多相流問題[5-9],涉及的流體主要是漿液和地下水。通常,漿液黏度和密度均大于地下水,且不同類型漿液和地下水之間存在不同程度的混溶現(xiàn)象,造成漿液和地下水之間的相互作用復雜。注漿過程中,漿液在注漿壓力驅(qū)動下排出導水通道中的地下水,除了漿液本身運移產(chǎn)生的摩擦外,地下水流動也將產(chǎn)生一定壓力損失。如圖1所示,采用常規(guī)單相流模型將導致漿液擴散范圍分析存在偏差。同時,傳統(tǒng)單相流注漿分析方法忽略漿-水相互作用,難以考慮漿液和地下水的物性差異。

圖1 傳統(tǒng)單相流模型漿液擴散范圍偏差Fig.1 Deviation of grout diffusion range predicted by traditional single-phase flow model

為了解決上述問題,Zou等[9]建立了單裂隙中漿-水驅(qū)替模型,得到了相界面壓力與擴散距離的計算公式。Ozdemir等[10]采用向量體積法研究了孔隙介質(zhì)中化學漿注漿過程的兩相擴散規(guī)律。劉人太等[11-12]通過室內(nèi)試驗和數(shù)值模擬研究了不同條件下漿-水兩相在裂隙中的擴散規(guī)律。楊志全等[13]基于Comsol Multiphysics 達西滲流模塊研究了冪律型漿液在礫石土多孔介質(zhì)中的兩相擴散規(guī)律。研究表明巖體中漿-水兩相具有明顯的非均勻擴散特征,漿-水之間的黏度差和復雜的導水通道網(wǎng)絡(luò)易誘發(fā)漿液擴散產(chǎn)生指進現(xiàn)象,降低注漿驅(qū)水效率。近年來,隨著CO2地質(zhì)封存,非常規(guī)油氣資源采收和地熱開采等項目的實施,復雜地層中多相流研究得到極大的關(guān)注。JU等[14-15]采用格子玻爾茲曼方法(LBM)模擬了巖體介質(zhì)非均勻孔隙系統(tǒng)內(nèi)的非混相水油驅(qū)替過程。Zakirov等[16]研究了不同潤濕條件和毛細管壓力下孔隙的各向異性對兩相驅(qū)替行為的影響。由于LBM方法具有適合處理巖土多孔介質(zhì)復雜邊界和易于并行的優(yōu)點,基于數(shù)字巖芯的LBM模擬技術(shù)在巖體多相流及流固耦合方向具有良好的應(yīng)用前景。

綜上分析,現(xiàn)有注漿模型雖然考慮了漿液黏度時變性、裂隙開度、注漿壓力等因素對漿液擴散范圍的影響,但對漿液驅(qū)替地下水效應(yīng)考慮不充分,難以處理復雜導水通道網(wǎng)絡(luò)中的漿-水相互作用,造成漿液有效擴散范圍分析結(jié)果偏大。漿-水兩相擴散范圍決定注漿帷幕交圈厚度和注漿帷幕堵水有效性,關(guān)系注漿堵水的成敗。目前,深部開采礦井水害防治問題越來越突出[17-18],漿-水兩相驅(qū)替運移規(guī)律是注漿堵水技術(shù)極為重要的基礎(chǔ)問題。因此,開展巖體注漿漿水兩相驅(qū)替機制和運移規(guī)律研究具有重大理論意義和工程價值。本文采用圖像處理方法識別溶蝕白云巖巖芯導水通道,利用三維重構(gòu)技術(shù)獲得導水通道三維數(shù)字模型。然后,基于Shan-Chen模型描述漿水兩相相互作用,采用LBM模擬巖體注漿過程中漿水兩相驅(qū)替行為,分析了漿相飽和度和殘余地下水的演化特征,揭示漿水驅(qū)替機制,為注漿防治礦井水害提供技術(shù)參考。

1 LBM基本原理及漿-水相互作用模型

1.1 LBM模型及原理

格子玻爾茲曼方法采用離散格式求解連續(xù)Boltzmann方程,是一種計算流體力學方法。該方法能夠適應(yīng)復雜邊界條件,同時具有程序?qū)嵤┫鄬唵魏鸵子诓⑿械葍?yōu)點[19]。如圖2所示,以二維格子模型D2Q9為例來介紹LBM的原理?？臻g流域離散為格子,時間離散為固定的時間步長。格子通過離散的節(jié)點聯(lián)系,節(jié)點包含粒子,粒子在節(jié)點內(nèi)運移并和相鄰節(jié)點產(chǎn)生接觸,實現(xiàn)動量和能量的交換。LBM并不追蹤單個分子,而是一群粒子,也即通過概率密度函數(shù)來代表粒子的運動。在節(jié)點上,流體的宏觀物理力學參數(shù)如密度和速度取決于分布函數(shù)和粒子遷移方向,計算公式為

圖2 D2Q9格子節(jié)點和遷移方向示意Fig.2 Sketch map of D2Q9 lattice and streaming direction of the particles

式中,ρ為格子節(jié)點上流體密度;fi為第i遷移方向上的分布函數(shù);u為格子節(jié)點上流體速度;ei為第i遷移方向的向量。

格子中粒子分布函數(shù)的遷移和碰撞過程如圖3所示,粒子分布函數(shù)演化方程表示為

圖3 D2Q9格子和分布函數(shù)碰撞遷移過程示意Fig.3 Sketch map of D2Q9 lattice and collide-streaming of the particles

式中,x為節(jié)點位置矢量;t為時間;Δt為時間步長;τ為松弛因子;feqi為平衡態(tài)分布函數(shù)。

LBM碰撞過程可以看作是分布函數(shù)向平衡狀態(tài)的松弛。D2Q9模型中,平衡狀態(tài)的分布函數(shù)feiq定義為

式中,c為基本格子速度,取1 lu·t/s;wi為第i遷移方向的權(quán)重,當i=0時,w0取4/9;當i=1,2,3,4時,wi取1/9;當i=5,6,7,8時,wi取1/36。

根據(jù)絕熱狀態(tài)方程,流體壓力取決于密度,計算公式為

式中,cs為格子聲速,cs=

采用D2Q9格子模型,流體的運動粘度為

當τ＞0.5時,代表實際物理粘度為正值。當τ→0.5時,會造成數(shù)值計算不穩(wěn)定。

1.2 漿水相互作用模型

Shan-Chen多組分模型是當前最常用的多相LBM模型之一,該模型通過引入偽勢來表征多相粒子間的相互作用。格子中,漿水兩相分別通過自身的粒子分布函數(shù)表示,當空間格子離散格式為DnQb時,組分k的分布函數(shù)演化方程為

式中,下標k代表流體組分,k=g、w,分別表示漿液和地下水;i為離散速度個數(shù),i=0,1,2,…,b-1;fki(x,t)為t時刻、位置x上組分k的粒子分布函數(shù);為組分k的平衡態(tài)分布函數(shù);τk為組分k的松弛時間,取決于流體的運動黏度,由υk=c2ks(τk-0.5)確定。

流體組分k的密度和速度采用對應(yīng)的分布函數(shù)計算,公式為

流體組分k的平衡態(tài)粒子分布函數(shù)為

式中,為組分k的平衡態(tài)速度。為了模擬漿水的相互作用,假設(shè)漿液和地下水兩相粒子之間存在相互作用勢,相應(yīng)的勢函數(shù)[19]為

式中,x'=x+ciΔt為與當前節(jié)點緊臨的格子節(jié)點位置;(x,x')滿足(x,x')=(x',x),決定了流體組分k和之間的相互作用強度;ψk是流體組分k局部密度的函數(shù),表示組分k的有效密度。

一般地,函數(shù)(x,x')只考慮最臨近格點的相互作用。

式中,δ為格子間距;參數(shù)絕對值的大小表征組分k和的粒子間的相互作用的強弱,的符號決定粒子間相互作用是相互吸引還是排斥,負值是相互吸引而正值則為相互排斥。構(gòu)成一個對稱的相互作用矩陣。當對角元素都為負值,而非對角元素都為正值時,同組份相互吸引而異組分相互排斥。

根據(jù)相互作用勢原理,組分k的粒子在x處受到的長程作用力與x處和x'處的有效密度的乘積成正比,計算公式為

粒子間相互作用力Fk對分布函數(shù)的影響是通過對組分k的平衡態(tài)速度ueqk的影響來體現(xiàn)的。

式中,u(x)是混合流體宏觀速度,其計算公式為

為了保證流體系統(tǒng)的動量守恒,平衡態(tài)速度進行了重新定義為

兩相流體的分離通過選取足夠強G的和函數(shù)ψ(ρ)實現(xiàn)。Shan和Chen等建議采用指數(shù)形式的有效密度函數(shù)

式中,ρk0是組分k的參考密度。

2 溶蝕白云巖導水通道識別及三維數(shù)字模型

2.1 導水通道三維數(shù)字模型重構(gòu)流程

為了獲得溶蝕白云巖復雜導水通道的三維結(jié)構(gòu),采用X射線微納米CT掃描、圖像處理技術(shù)和三維重構(gòu)技術(shù)構(gòu)建巖樣導水通道三維數(shù)字模型,為漿水兩相運移規(guī)律數(shù)值模擬提供流場的幾何文件。如圖4所示,首先,采用X射線微納米CT掃描技術(shù)獲得巖樣高分辨率CT掃描圖片;然后,通過圖像處理技術(shù)并運用閾值分割識別巖樣CT掃描圖片中的導水孔隙和裂隙,將圖片像素分為固體基質(zhì)、接觸邊界和流體區(qū)域;最后,依據(jù)處理好的CT掃描圖像,采用Matlab代碼逐層重構(gòu)巖石導水通道三維數(shù)字模型,為巖樣中漿水兩相運移特征LBM模擬提供幾何文件。

圖4 巖石導水通道三維重構(gòu)示意Fig.4 Sketch map of 3D reconstruction of water conducting channels in rock

2.2 溶蝕白云巖試樣微納米CT掃描源文件

溶蝕白云巖作為特殊的多孔介質(zhì),其賦水特征和滲透特性與常規(guī)裂隙巖體不同。對于溶蝕孔隙型巖體,小孔、微孔是地下水賦存的主要存儲空間,而中等及以上溶孔構(gòu)成的孔隙網(wǎng)絡(luò)是地下流體運移的主要通道。本文采用的溶蝕白云巖巖樣的微納米CT掃描源文件包含288張CT掃描切片,分辨率為250 μm。x、y、z方向上像素分別為791、791、100。溶蝕白云巖巖樣孔隙切片及孔隙率分布情況如圖5所示。

圖5 溶蝕白云巖切片及孔隙率分布Fig.5 Slice of the dissolved dolomite sample and its porosity distribution

2.3 白云巖試樣導水通道三維數(shù)字模型

(1)圖像處理與導水通道識別。為了識別試樣導水通道,首先將圖片進行濾波去噪處理,然后利用圖片處理軟件ImageJ通過閾值分割法來識別巖樣中的孔隙結(jié)構(gòu)。如圖6(a)所示,通過選擇合適的閾值,準確識別切片中的孔隙,然后將圖片二值化,得到圖6(b)所示的二值化切片圖,除邊界輪廓外,圖片中白色區(qū)域部分為孔隙空間,灰度值為255,邊界輪廓內(nèi)部的黑色區(qū)域為巖石基質(zhì),灰度值為0。將二值化的CT掃描切片圖像重疊,就可得到溶蝕白云巖巖樣的三維孔隙網(wǎng)絡(luò)結(jié)構(gòu),如圖6(c)所示。

圖6 孔灰度圖片閾值分割與二值化Fig.6 Threshold segmentation and binarization of the grayscale image

(2)導水通道數(shù)字模型三維重構(gòu)。由于配備的工作站處理能力有限,選擇巖芯局部區(qū)域進行孔隙巖樣三維數(shù)字模型三維重構(gòu)。如圖7(a)所示,選擇的研究區(qū)域?qū)?50像素,高100像素,切片張數(shù)為100張。由于后續(xù)模擬軟件中的流動方向為x方向,也即垂直切片的縱向方向。因此,在三維重構(gòu)中,設(shè)置圖片寬度方向為y方向,圖片高度方向為z方向?？紫稁r樣數(shù)字模型三維重構(gòu)程序采用Matlab代碼,如果某一體素為流體區(qū)域,則標記為0;固體基質(zhì)體素則標記為2,而固體和流體的接觸層體素標記為1,獲得三維孔隙網(wǎng)絡(luò)數(shù)字模型如圖7(b)所示。

圖7 溶蝕白云巖巖樣三維數(shù)字模型Fig.7 Three-dimensional digital model of dissolved dolomite rock sample

3 溶蝕白云巖注漿漿-水兩相運移特征

3.1 模擬軟件簡介

溶蝕白云巖注漿過程中漿水兩相運移特征采用Palabos軟件進行模擬分析,模擬運算采用Linux計算工作站,配置Xeon(R)Gold 6230 雙CPU,40核心,128G內(nèi)存。Palabos是以LBM為內(nèi)核的通用CFD開源分析軟件,該軟件編程接口界面采用C++編寫,內(nèi)置了單相流、多相流以及多場耦合模型,同時包含多種碰撞格式的類可供用戶調(diào)用,具有處理復雜邊界條件和并行計算效率高等優(yōu)點。由于該軟件沒有圖形界面,采用第三方軟件ParaView進行模擬結(jié)果的三維可視化分析。

3.2 初始和邊界條件及漿-水兩相基本參數(shù)

模擬分析的幾何文件采用前述溶蝕白云巖巖樣三維數(shù)字模型,每一個體素為一個格子單元,格子遷移速度模型采用D3Q19模型,如圖8(a)所示。漿-水兩相初始條件及邊界如圖8(b)所示。模擬流動方向為x方向,左側(cè)為流體入口,右側(cè)為流體出口,模型四壁和內(nèi)部巖石基質(zhì)格子為反彈邊界條件,也即無滑移邊界條件。模擬中,采用Shan-Chen模型描述漿-水兩相間相互作用,漿水兩相基本物理力學參數(shù)列于表1。

表1 巖體注漿模擬中漿-水物理力學參數(shù)Table 1 Pysical and mechanic perameters of the grout amd water

圖8 漿-水兩相運移模擬的格子模型以及初始和邊界條件Fig.8 Lattice model and initial & boundary conditions for simulation of grout water two-phase transport

3.3 模擬結(jié)果分析

(1)孔隙巖樣滲透率分析。采用溶蝕白云巖巖樣的三維數(shù)字模型可進行巖樣絕對滲透率分析。滲透率是表征巖土介質(zhì)傳輸流體能力的關(guān)鍵參數(shù)。巖樣絕對滲透率與孔隙度、孔隙幾何形狀、迂曲情況、巖土顆粒大小及排列方式等因素有關(guān)。圖9為白云巖巖樣孔隙連通情況和滲透率模擬結(jié)果。從圖9(a)可以看出,白云巖試樣孔隙發(fā)育,孔隙大小和形狀分布復雜,通過軟件統(tǒng)計分析得到研究區(qū)域的平均孔隙率為0.42,但在x方向上孔隙連通程度相對較差,存在較多的孤立孔隙區(qū)域。圖9(b)為試樣中模擬的流動速度分布情況,可以看出,流體選擇沿連通性較好的大孔隙運移,同時,孔隙中部流速明顯大于孔隙壁面附近流速,孔隙壁面摩擦對流體流動產(chǎn)生明顯的抑制。試樣的真實物理滲透率是格子滲透率乘以空間分辨率的平方。模擬得到的格子單位滲透率為kl=0.095 8 lu2,空間分辨率為23 μm/像素,換算后可以得到研究區(qū)域的真實絕對滲透率為kp= 50.70 μm2,屬高滲透率巖樣。

圖9 Guelph 白云巖巖樣滲透率模擬結(jié)果Fig.9 Simulated permeability of Guelph dolomite rock sample

(2)孔隙巖樣中漿水兩相運移特征。設(shè)漿相到達出口的時間點為t0,如圖10所示,選擇4個不同的相對時刻來分析漿水兩相運移分布情況。在注漿壓力作用下,漿液開始驅(qū)替地下水。初始階段,入口處孔隙結(jié)構(gòu)和孔徑大小對漿相運移具有重要影響。漿液易沿大孔中部向前推進,形成平緩的相界面。由于孔徑網(wǎng)絡(luò)連通情況的不同,漿相前沿推進距離不同。小孔和連通性較差的孔隙中,漿相前沿明顯滯后。驅(qū)替初始階段,孔隙壁面附近殘存一定地下水,隨著注漿時間的增加,孔隙壁面附近地下水逐步減少。漿液驅(qū)替地下水行為具有明顯的非均質(zhì)性(各向異性),漿液擴散具有明顯的選擇性,連通性較好的孔隙為漿相主要運移通道,漿相最先達到出口。由于漿液不能完全排出孔隙壁面附近地下水,同時,后續(xù)漿液固結(jié)過程中可能會析出少量地下水,造成漿液固結(jié)體與孔隙壁面不能有效接觸,是注漿堵水的薄弱環(huán)節(jié)?？紫稁r樣中孤立孔隙中的地下水幾乎不能被漿液驅(qū)替。被束縛的地下水降低了整個孔隙巖樣的滲透性,后續(xù)漿液沿著已經(jīng)驅(qū)替形成的優(yōu)勢通道運移擴散,降低了注漿驅(qū)水效率。

圖10 Guelph 白云巖巖樣中漿-水兩相驅(qū)替行為模擬結(jié)果Fig.10 Displacement behavior of grout and water in Guelph dolomite rock sample

4 結(jié) 論

(1)溶蝕白云巖試樣軸向方向上孔隙連通程度相對較差,存在較多的孤立孔隙區(qū)域,試樣平均孔隙率為0.42。注漿過程中,漿液選擇沿連通性較好的大孔隙運移,孔隙中部流速明顯大于孔隙壁面附近流速,孔隙壁面摩擦對流體流動產(chǎn)生明顯的抑制。

(2)驅(qū)替初始階段,孔隙壁面附近殘存一定地下水,隨著注漿時間的增加,孔隙壁面附近地下水逐步減少。漿液驅(qū)替地下水行為具有明顯的非均質(zhì)性,漿液擴散具有明顯的選擇性,連通性較好的孔隙為漿相主要運移通道,漿相最先達到出口,孤立孔隙區(qū)域中的地下水基本不能被漿液排出置換。

(3)被束縛的地下水降低了整個孔隙巖樣的滲透性,后續(xù)漿液沿著已經(jīng)驅(qū)替形成的優(yōu)勢通道運移擴散,降低了注漿驅(qū)水效率。單次注漿漿液很難完全排出孔隙壁面附近地下水,而且漿液固結(jié)過程中也會析出少量地下水,造成漿液固結(jié)體與孔隙壁面之間很難有效接觸,是注漿堵水的薄弱部位。