巖體孔隙-裂隙雙滲流數(shù)值模擬研究

邵建立，周 斐，薛彥超，杜后謙

（山東科技大學 礦業(yè)與安全工程學院，山東 青島266590）

巖體滲流一直是礦山、水利水電、建筑等巖土工程的重要問題，流體在巖體裂隙中快速運移，也會相對緩慢的通過周圍基質(zhì)塊中微小孔隙遷移。采取合理的防滲措施是防控巖體滲流危害的有效手段，而準確地選取理論模型進行計算和模擬是預防和消除巖體滲流影響的關(guān)鍵[1-2]。基于巖體具有裂隙和基質(zhì)雙重滲流過程，研究孔隙-裂隙雙重介質(zhì)滲流場發(fā)展變化規(guī)律尤為重要。

國內(nèi)外學者已經(jīng)針對裂隙巖體滲流特征進行了許多相關(guān)的研究。朱斌[3]等結(jié)合開灤趙各莊礦14 水平開拓東大巷揭露的薄層煤巖體滲流演化過程進行數(shù)值模擬，通過調(diào)節(jié)滲透系數(shù)，獲得了薄層煤巖體裂隙-孔隙雙滲流在時間和空間上的孔隙水壓變化過程；速寶玉[4-6]等通過實驗研究裂隙巖體滲流應(yīng)力耦合情況，闡明了單裂隙面的各種經(jīng)驗公式、間接公式及其適用條件，分析了裂隙巖體滲流應(yīng)力耦合模型優(yōu)缺點及目前工程應(yīng)用情況。李琛亮[7]等研制的雙重介質(zhì)滲流水力特性試驗系統(tǒng)，研究了基于雙重介質(zhì)模型的水量交換以及滲流場的水壓分布規(guī)律以及雙重介質(zhì)的水力性態(tài)和滲流機制，得出孔隙-裂隙雙重介質(zhì)水交換影響因子對雙重介質(zhì)水交換的影響能力；國外Barenblatt 提出均質(zhì)、各向同性的孔隙-裂隙雙重介質(zhì)概念[8]，后續(xù)學者們開展了孔隙-裂隙雙重介質(zhì)模型及其解析和數(shù)值算法[9-10]。Samardzioska[11]比較了巖體等效介質(zhì)模型、裂隙網(wǎng)絡(luò)模型和裂隙－孔隙雙重介質(zhì)模型的滲流演化規(guī)律，獲得了不同介質(zhì)假設(shè)下巖體滲流演化對比研究成果。

由于巖體內(nèi)部不可視性和裂隙網(wǎng)絡(luò)錯綜復雜，學者們難以可視化地揭露內(nèi)部孔隙-裂隙滲流規(guī)律變化。基于此，采用多孔介質(zhì)滲流數(shù)值模擬的方法，建立斷裂的多孔介質(zhì)塊模型，通過不同形狀路徑的裂隙，對孔隙-裂隙雙重介質(zhì)滲流場分布進行模擬研究，以期為揭露巖體滲流規(guī)律、滲流危害預測與防治提供一定的理論支持。

1 孔隙-裂隙滲流數(shù)值模擬

1.1 基本控制方程

1.1.1 孔隙滲流方程

流體在均質(zhì)的基質(zhì)塊中滲流遵循達西定律[9]，隨時間變化的方程為：

式中：u 為速度矢量，m/s；p 為孔隙水壓力，Pa；εp為基質(zhì)塊的孔隙率；S 為基質(zhì)塊儲水系數(shù)，1/Pa；ρ為密度，kg/m3；t 為時間。

線性儲水模型為：

式中：Xf為流體壓縮率，1/Pa；Xp為基質(zhì)塊等效壓縮率，1/Pa。

在塊內(nèi)，內(nèi)建的速度變量u 給出達西速度，達西速度是多孔介質(zhì)單位面積的體積流量。

式中：k 為基質(zhì)塊的滲透率，m2；μ 為流體動力黏度，Pa·s。

1.1.2 裂隙滲流方程

使用COMSOL 裂縫流邊界條件，允許沿著內(nèi)部邊界或裂隙定義流動。在這種邊界條件下，裂隙速度方程遵循基質(zhì)塊內(nèi)速度方程（即達西定律）的修正形式。考慮到裂隙對流動阻力較小，裂隙厚度較小，使得裂縫與基體的尺寸一致性，對達西定律修正，得到以下方程：

式中：Sf為裂隙儲水系數(shù)，1/Pa；kf為裂隙的滲透率，m2；df為裂隙厚度，m；▽T為裂隙切向平面的梯度算子。

由于裂縫流動方程中含有厚度，內(nèi)建變量uf給出了裂隙單位長度的體積流量：

式中：uf為裂隙速度矢量，m/s。

1.2 模擬方案

使用COMSOL 數(shù)值模擬軟件對孔隙-裂隙雙滲流進行數(shù)值模擬。巖體孔隙-裂隙雙滲流數(shù)值模擬模型如圖1。物理模型為4 種斷裂的均質(zhì)的多孔介質(zhì)塊，塊體每邊長度為1 m。塊中為不同形狀路徑的裂隙，依次為90°夾角型、45°夾角型、135°夾角型、圓角型。與孔隙滲流相比，裂隙對流體的滲透性更強，同時裂隙厚度為0.1 mm，遠小于塊的尺寸。流體從右向左移動，通過塊進入裂隙下部邊界并從上部邊界離開。流體最初不在塊內(nèi)滲流。出口邊界處的壓力隨時間下降，而入口邊緣處的壓力在整個模擬過程中保持初始壓力。

除了在裂隙邊界之外，基質(zhì)塊的壁是不可滲透的。沿塊體的所有面應(yīng)用0 流量邊界條件：

式中：n 為向外指向邊界的法向量。

在裂隙入口和裂隙出口，采用壓力邊界條件，關(guān)系式如下（0≤t≤1 000 s）：

式中：p 為裂隙出口壓力，Pa；p0為裂隙入口壓力，Pa；t 為時間，s；a 為壓力變化率，Pa/s。

本次數(shù)值模擬的相關(guān)參數(shù)和賦值見表1。

2 模擬結(jié)果分析

2.1 裂隙流動分析

以90°角裂隙模擬為例分析裂隙流動變化。不同時刻裂隙路徑速度分布如圖2。隨著時間的變化，出口邊界的壓力呈線性減少，而入口邊界的壓力保持恒定p0，產(chǎn)生的壓力差驅(qū)使流體流動，裂隙上的速度分布逐漸改變。可以發(fā)現(xiàn)，計算剛開始速度分布均勻，隨著壓力差增大，速度場也發(fā)展增大。當計算達到1 000 s，入口邊界和出口邊界壓力差達到最大值，速度場也發(fā)展到最大值。其次，流體從裂隙入口到出口，其速度是線性連續(xù)的，其中在向上導升過程中部分動能轉(zhuǎn)化為勢能，速度變化微小，而在入口邊界和出口邊界的速度始終是流場中的最大值。

表1 數(shù)值模擬相關(guān)參數(shù)和賦值

圖2 不同時刻裂隙路徑速度分布

2.2 基質(zhì)孔隙流動分析

以90°角裂隙模擬為例分析內(nèi)部流動。不同時刻基質(zhì)塊內(nèi)部壓力等值面分布如圖3。隨著時間變化，出口壓力逐漸減小，整個基質(zhì)塊內(nèi)部孔隙壓力重新分布，壓力等值面分級增多，高壓等值面靠近入口，低壓等值面靠近出口，基質(zhì)塊內(nèi)部壓力梯度越來越顯著。同時等壓面穿過裂隙，表明壓力分布在斷裂的基質(zhì)塊中是連續(xù)的，因此裂隙上壓力分布也是連續(xù)。但在壓力等值面與裂隙相交處具有不同程度的彎折，聯(lián)合圖2 可發(fā)現(xiàn)，裂隙速度越大壓力等值面彎曲越明顯，表明流體在裂隙和孔隙中的流動狀態(tài)不同。

2.3 不同路徑裂隙流動對比分析

對于4 種不同斷裂形狀的多孔介質(zhì)塊，保持相同初始參數(shù)和邊界條件模擬，分別在4 種模型種選取x=0.5 m 處yz 截面，以1 000 s 計算為例，繪制速度等值線并填充（圖4）。

圖3 不同時刻基質(zhì)塊內(nèi)部壓力等值面分布

由圖4 可以看出，速度場基本在截面上呈中心對稱分布，且不同路徑的裂隙對孔隙滲流速度影響是不同的，在45°、90°、135°角形狀裂隙路徑中，夾角處均出現(xiàn)了較小的相對高速區(qū)域，圖4 圓角則沒有出現(xiàn)這種相對高速區(qū)域，原因是在夾角處流體的流動方向突然發(fā)生改變，產(chǎn)生了局部阻力損失，而緩和的路徑減少了這種能量損失，圓角型裂隙流體在出口處速度比其他路徑裂隙出口速度都高，達到2.59×10-4m/s。

圖4 不同斷裂形狀的基質(zhì)塊截面速度分布

沿不同形狀路徑裂隙出口邊界上，不同時刻流體邊界通量分布如圖5。流體邊界通量是指單位時間內(nèi)流經(jīng)邊界單位面積的物質(zhì)量，是表示輸送強度的物理量。本次模擬中，除了裂隙形狀不同，其余參數(shù)保持相同，因此出口邊界通量可以反映不同路徑裂隙情況下流體流動的強度或速度。從圖5 中看出，起始階段不同形狀裂隙的邊界通量相差微小甚至相交，因為初始階段入口和出口的壓力差不大，流體運動緩慢，隨著時間增加，速度場逐漸發(fā)展，不同路徑裂隙滲流的能量損失逐漸明顯，同時不同路徑裂隙的邊界通量也呈線性增長趨勢。當計算至1 000 s時，圓角型和135°夾角型路徑裂隙具有較高的邊界通量，表示滲流過程中流體能量損失較小，而圓角型路徑裂隙邊界通量最高，說明圓角型路徑裂隙流體能量損失最小；90°夾角型和45°夾角型路徑裂隙邊界通量較小，表示這2 種路徑滲流過程能量損失較大，其中45°夾角型邊界通量最小，說明滲流過程能量損失最多。

圖5 不同時刻4 種裂隙出口邊界通量變化

3 結(jié) 論

1）裂隙是滲流的主要途徑，在壓力充足的情況下，流體在裂隙流動發(fā)展最充分，滲流場在裂隙路徑上連續(xù)分布。在礦山深部和地下巖土工程等具有高水壓威脅的地方，裂隙滲流的影響不可忽視，應(yīng)從裂隙滲流角度防控災(zāi)害。

2）各向同性的多孔介質(zhì)內(nèi)部孔隙壓力梯度分布均勻，斷裂的多孔介質(zhì)塊壓力分布是連續(xù)的，裂隙和孔隙存在流體交換，這取決于流場內(nèi)壓力分布。由于裂隙和孔隙具有不同的滲透率等因素，流體在裂隙和孔隙的流動狀態(tài)明顯不同。

3）不同形狀路徑的裂隙滲流場也有差異。保持相同初始參數(shù)，4 種路徑裂隙出口邊界流體通量在起始階段相差微小，隨著時間變化速度場充分發(fā)展，這種差距愈加明顯，到1 000 s 時，4 種裂隙出口邊界通量的大小關(guān)系為：圓角型＞135°夾角＞90°夾角＞45°夾角，說明4 種路徑裂隙的能量損失大小關(guān)系為：45°夾角＞90°夾角＞135°夾角＞圓角型。