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

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

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

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

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

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

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

1.1 基本控制方程

1.1.1 孔隙滲流方程

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

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

線性儲(chǔ)水模型為：

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

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

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

1.1.2 裂隙滲流方程

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

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

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

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

1.2 模擬方案

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

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

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

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

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

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

2 模擬結(jié)果分析

2.1 裂隙流動(dòng)分析

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

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

圖2 不同時(shí)刻裂隙路徑速度分布

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

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

2.3 不同路徑裂隙流動(dòng)對(duì)比分析

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

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

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

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

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

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

3 結(jié) 論

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

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

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