周禹良 袁東鋒 楊 雪 李 斌

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

礦井水害具有強烈的致災性,若發(fā)生涌水量過大或突水淹井事故,將造成巨大的經(jīng)濟損失和惡劣的社會影響[1-2]。注漿技術是一種主動防治礦井水害的有效手段[3-4],是通過一定壓力將漿液注入地下含水層,驅排地層中的地下水,漿液凝固后充填封堵地層導水通道形成注漿帷幕[5-6],提高巖體的黏聚力和內(nèi)摩擦角[7-8],以達到堵水和加固地層的目的。巖體注漿是典型的多相流問題[5-9],涉及的流體主要是漿液和地下水。通常,漿液黏度和密度均大于地下水,且不同類型漿液和地下水之間存在不同程度的混溶現(xiàn)象,造成漿液和地下水之間的相互作用復雜。注漿過程中,漿液在注漿壓力驅動下排出導水通道中的地下水,除了漿液本身運移產(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]建立了單裂隙中漿-水驅替模型,得到了相界面壓力與擴散距離的計算公式。Ozdemir等[10]采用向量體積法研究了孔隙介質中化學漿注漿過程的兩相擴散規(guī)律。劉人太等[11-12]通過室內(nèi)試驗和數(shù)值模擬研究了不同條件下漿-水兩相在裂隙中的擴散規(guī)律。楊志全等[13]基于Comsol Multiphysics 達西滲流模塊研究了冪律型漿液在礫石土多孔介質中的兩相擴散規(guī)律。研究表明巖體中漿-水兩相具有明顯的非均勻擴散特征,漿-水之間的黏度差和復雜的導水通道網(wǎng)絡易誘發(fā)漿液擴散產(chǎn)生指進現(xiàn)象,降低注漿驅水效率。近年來,隨著CO2地質封存,非常規(guī)油氣資源采收和地熱開采等項目的實施,復雜地層中多相流研究得到極大的關注。JU等[14-15]采用格子玻爾茲曼方法(LBM)模擬了巖體介質非均勻孔隙系統(tǒng)內(nèi)的非混相水油驅替過程。Zakirov等[16]研究了不同潤濕條件和毛細管壓力下孔隙的各向異性對兩相驅替行為的影響。由于LBM方法具有適合處理巖土多孔介質復雜邊界和易于并行的優(yōu)點,基于數(shù)字巖芯的LBM模擬技術在巖體多相流及流固耦合方向具有良好的應用前景。

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

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

1.1 LBM模型及原理

格子玻爾茲曼方法采用離散格式求解連續(xù)Boltzmann方程,是一種計算流體力學方法。該方法能夠適應復雜邊界條件,同時具有程序實施相對簡單和易于并行等優(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遷移方向的權重,當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的密度和速度采用對應的分布函數(shù)計算,公式為

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

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

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

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

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

根據(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ù)字模型重構流程

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

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

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

溶蝕白云巖作為特殊的多孔介質,其賦水特征和滲透特性與常規(guī)裂隙巖體不同。對于溶蝕孔隙型巖體,小孔、微孔是地下水賦存的主要存儲空間,而中等及以上溶孔構成的孔隙網(wǎng)絡是地下流體運移的主要通道。本文采用的溶蝕白云巖巖樣的微納米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通過閾值分割法來識別巖樣中的孔隙結構。如圖6(a)所示,通過選擇合適的閾值,準確識別切片中的孔隙,然后將圖片二值化,得到圖6(b)所示的二值化切片圖,除邊界輪廓外,圖片中白色區(qū)域部分為孔隙空間,灰度值為255,邊界輪廓內(nèi)部的黑色區(qū)域為巖石基質,灰度值為0。將二值化的CT掃描切片圖像重疊,就可得到溶蝕白云巖巖樣的三維孔隙網(wǎng)絡結構,如圖6(c)所示。

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

(2)導水通道數(shù)字模型三維重構。由于配備的工作站處理能力有限,選擇巖芯局部區(qū)域進行孔隙巖樣三維數(shù)字模型三維重構。如圖7(a)所示,選擇的研究區(qū)域寬150像素,高100像素,切片張數(shù)為100張。由于后續(xù)模擬軟件中的流動方向為x方向,也即垂直切片的縱向方向。因此,在三維重構中,設置圖片寬度方向為y方向,圖片高度方向為z方向。孔隙巖樣數(shù)字模型三維重構程序采用Matlab代碼,如果某一體素為流體區(qū)域,則標記為0;固體基質體素則標記為2,而固體和流體的接觸層體素標記為1,獲得三維孔隙網(wǎng)絡數(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)置了單相流、多相流以及多場耦合模型,同時包含多種碰撞格式的類可供用戶調用,具有處理復雜邊界條件和并行計算效率高等優(yōu)點。由于該軟件沒有圖形界面,采用第三方軟件ParaView進行模擬結果的三維可視化分析。

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

模擬分析的幾何文件采用前述溶蝕白云巖巖樣三維數(shù)字模型,每一個體素為一個格子單元,格子遷移速度模型采用D3Q19模型,如圖8(a)所示。漿-水兩相初始條件及邊界如圖8(b)所示。模擬流動方向為x方向,左側為流體入口,右側為流體出口,模型四壁和內(nèi)部巖石基質格子為反彈邊界條件,也即無滑移邊界條件。模擬中,采用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 模擬結果分析

(1)孔隙巖樣滲透率分析。采用溶蝕白云巖巖樣的三維數(shù)字模型可進行巖樣絕對滲透率分析。滲透率是表征巖土介質傳輸流體能力的關鍵參數(shù)。巖樣絕對滲透率與孔隙度、孔隙幾何形狀、迂曲情況、巖土顆粒大小及排列方式等因素有關。圖9為白云巖巖樣孔隙連通情況和滲透率模擬結果。從圖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 白云巖巖樣滲透率模擬結果Fig.9 Simulated permeability of Guelph dolomite rock sample

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

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

4 結 論

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

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

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