基于有限元的不同粗糙度裂隙等效滲透性預(yù)測

王夢輝，涂福彬*，童俊，焦玉勇，鄭翔

(1.中國地質(zhì)大學(xué)(武漢)工程學(xué)院，武漢 430074；2.中建一局集團(tuán)建設(shè)發(fā)展有限公司，北京 100102)

滲透性是裂隙巖體的一個重要特性，是分析地下油氣運(yùn)移與污染物遷移等滲透問題的基礎(chǔ)。在巖體中由于裂隙的滲透性遠(yuǎn)遠(yuǎn)大于巖塊的滲透性，因此裂隙是巖體中流體運(yùn)移的主要通道，裂隙滲透性的研究足以確定裂隙巖體的滲透性。在實際巖體中，裂隙往往相互交錯形成網(wǎng)狀結(jié)構(gòu)，裂隙巖體的滲透性研究存在一定的困難，對巖體中單一裂隙滲透規(guī)律的研究是了解復(fù)雜裂隙巖體滲透性的基礎(chǔ)。

在巖體滲透性預(yù)測方面，Snow[1]通過光滑平行板裂隙水流實驗提出了著名的立方定律，即認(rèn)為通過裂隙的滲流量與裂隙寬度的三次方成正比。而自然界實際巖體的裂隙面往往是彎曲粗糙的，光滑平行板模型在表征裂隙表面時存在誤差[2-3]。任志善等[4]通過對數(shù)值分析模型研究發(fā)現(xiàn)裂隙傾角、間距、隙寬和跡長對巖體的滲透特性均有明顯影響。當(dāng)流體流經(jīng)粗糙裂隙時，由于節(jié)理表面粗糙度、接觸面積的影響，流體滲流路徑是曲折的[5]。朱紅光等[6]研究表明，立方定律只適用于壁面粗糙度足夠小的情況。為此，許多學(xué)者對立方定律進(jìn)行了修正。Rose等[7]提出了迂曲度的概念，隨后基于迂曲度修正的立方定律計算公式被用于預(yù)測含天然裂隙巖體的滲透性。夏才初等[8]采用節(jié)理表面分形維數(shù)和空腔組合形貌等參數(shù)來計算通過粗糙裂隙的滲流流量。肖維民等[9-12]通過歸納實驗結(jié)果，提出了考慮曲折效應(yīng)的粗糙節(jié)理滲流流量計算新公式。速寶玉等[13]對含天然裂隙的巖體滲流進(jìn)行了研究。覃源等[14]采用實驗手段研究了不同節(jié)理粗糙度對單裂隙滲流的影響。段慕白等[15]通過實驗研究了裂隙粗糙度及隙寬對滲流的影響。李崴等[16]通過數(shù)值模擬方法研究了輻射流情況下裂隙巖體的滲透性。王瑋等[17]采用基于正、負(fù)相關(guān)指標(biāo)的復(fù)合因子模型評估裂隙巖體的滲透系數(shù)。

采用數(shù)值方法求解不可壓縮黏性流體穩(wěn)態(tài)流動的Navier-Stokes方程是近年來精細(xì)確定裂隙巖體滲透性的主要手段。Nguyen等[18]采用快速傅里葉變換方法求解細(xì)觀表征體元內(nèi)不可壓縮黏性流體穩(wěn)態(tài)流動的Navier-Stokes方程，并通過平均流速和壓強(qiáng)梯度之比確定多孔介質(zhì)的滲透系數(shù)。王志良等[19]、張戈等[20]采用格子Boltzmann方法模擬巖體裂隙內(nèi)的流體流動。上述方法只模擬了一定長度裂隙內(nèi)的流體流動，忽略了裂隙延伸長度對滲流的影響。鐘振等[21]采用數(shù)值模擬手段發(fā)現(xiàn)粗糙單裂隙的滲透性存在明顯的尺寸效應(yīng)。Zhang等[22]通過微裂隙滲流實驗也發(fā)現(xiàn)裂隙的滲透性存在尺寸效應(yīng)。

為此，有必要施加周期性邊界條件(periodic boundary conditions,PBC)來消除裂隙延伸長度所引起的誤差。首先采用不可壓縮單元離散不可壓縮流體的穩(wěn)態(tài)納維-斯托克斯方程，并施加周期性邊界條件構(gòu)造裂隙巖體滲透性的有限元分析方法。接著采用該方法分析了不同節(jié)理面粗糙度系數(shù)(joint roughness coefficient,JRC)裂隙的滲透性。

1 裂隙滲流的有限元分析方法

1.1 控制方程

通過求解Navier-Stokes方程，可以獲得流體在裂隙內(nèi)的流動速度[18]。對于滲流問題，流體流動往往處于恒定層流狀態(tài)。假定流體不可壓縮，忽略體積力，取長度為L的表征體元(representative volume element,RVE)，如圖1所示。

圖1 無限長平行板周期性結(jié)構(gòu)示意圖Fig.1 Periodic structure diagram of infinite length parallel plate

不可壓縮黏性流體穩(wěn)態(tài)流動的Navier-Stokes方程為

(1)

?·u=0

(2)

對表征體元施加周期性邊界條件，則有

u(x+L,y)=u(x,y)

(3)

(4)

式中：x為平行于無限長平行板的方向；y為垂直于無限長平行板的方向。

1.2 有限元離散

采用四節(jié)點等參數(shù)單元對速度和周期性分布的壓強(qiáng)進(jìn)行插值，并得離散方程為

(5)

式(5)中：K為剛度矩陣；G為耦合矩陣；F為外力等效節(jié)點力向量。

(6)

(7)

(8)

式(5)的數(shù)值穩(wěn)定形式為

(9)

式(9)中：π為保證數(shù)值求解穩(wěn)定而引入的偽變量;L為周期性壓強(qiáng)對應(yīng)的剛度矩陣;Q為周期性壓強(qiáng)和穩(wěn)定變量的耦合剛度矩陣;M為穩(wěn)定變量對應(yīng)的剛度矩陣。

式(9)可以簡化為

(10)

(11)

(12)

(13)

基于式(10)～式(13)，編寫了ABAQUS子單元程序?qū)崿F(xiàn)裂隙滲流的有限元求解，并通過Python腳本施加周期性邊界條件。

2 不同JRC值裂隙滲流模型構(gòu)建

采用不同JRC值來描述裂隙表面粗糙度，探究裂隙面粗糙度系數(shù)、裂隙寬度以及裂隙兩表面相互錯動對滲流的影響。國際巖石力學(xué)學(xué)會推薦給出的10級JRC值對應(yīng)的粗糙度剖面曲線如圖2[23]所示，JRC值為0表示裂隙面為光滑平面。

0～10 cm為裂隙的長度圖2 不同JRC值范圍的粗糙度剖面[23]Fig.2 Roughness profiles with different JRC values[23]

地質(zhì)構(gòu)造的作用極易引起裂隙兩表面的相互錯動，而錯動將影響裂隙的寬度，進(jìn)而改變裂隙內(nèi)部流體的流速，影響滲透性。JCR值為18～20時，兩表面相互錯動的裂隙模型如圖3所示。

圖3 兩表面相互錯動的裂隙模型Fig.3 Fracture model with dislocation of two surfaces

Rose等[7]提出了迂曲度的概念來表征滲流路徑的曲折特性，迂曲度τ為沿著裂隙走向由一端走至末尾的實際長度Le與裂隙兩端直線距離L之比，可表示為

(14)

將立方定律公式采用迂曲度修正之后得

(15)

3 結(jié)果與討論

3.1 有限元模型驗證

現(xiàn)從無限長的光滑平行板中選取表征體元，如圖1所示。通過施加周期性邊界條件來計算無限長的光滑裂隙間的滲流。并將有限元計算的結(jié)果與泊肅葉流動的解析解進(jìn)行對比。泊肅葉流動的流速分布為

(16)

選取參數(shù)如表1所示。將有限元方法計算得到的結(jié)果與泊肅葉流動的解析解進(jìn)行對比，如圖4所示，有限元計算結(jié)果與解析解相吻合。可見，所提出的裂隙滲流有限元計算方法是可行的。

表1 模型參數(shù)Table 1 Parameters of the model

圖4 有限元方法與解析解對比Fig.4 Comparison between result of finite element method and analytical solution

3.2 周期性邊界條件對滲流的影響

對JRC值10～12隙寬1 mm的一條裂隙在水力坡度為0.01的情況下進(jìn)行有限元分析，比較有無周期性邊界條件計算結(jié)果的差異，如圖5所示。不施加周期性邊界條件裂隙左右兩側(cè)細(xì)觀上的流速不同，這與沿走向延伸較長的裂隙中滲流不符合，施加周期性邊界條件對延伸較長的裂隙非常重要。

圖5 有無周期性邊界條件結(jié)果對比Fig.5 Comparison of results with and without periodic boundary conditions

3.3 JRC值對滲流的影響

通過模擬隙寬3 mm的不同JRC值裂隙內(nèi)流體流動，探究相同水力坡度下節(jié)理面粗糙程度對通過裂隙的流體平均流速和最大流速的影響。不同JRC值裂隙內(nèi)滲流的最大流速和平均流速如圖6所示。隨著裂隙表面粗糙程度的增大，滲流的平均流速整體上呈逐漸減小的趨勢。細(xì)觀上滲流的最大流速并不是隨著節(jié)理面粗糙程度的增大而逐漸增大，而是與節(jié)理的局部起伏角相關(guān)，這與覃源等[14]的研究結(jié)果相似。

圖6 不同JRC值裂隙內(nèi)流體流動情況(隙寬3 mm)Fig.6 Fluid flow in fracture with different JRC values (crack width 3 mm)

對不同水力坡度下隙寬3 mm的不同JRC值的裂隙進(jìn)行有限元分析，計算不同水力坡度下通過不同JRC值裂隙的單寬流量值，結(jié)果如圖7所示。通過不同JRC值裂隙的單寬滲流量與水力坡度成正比，符合達(dá)西定律。

3.4 迂曲度對滲流的影響

計算隙寬3 mm的不同JRC值裂隙在水力坡度為0.01時的滲流量，將有限元計算結(jié)果與引入迂曲度修正后的立方定律公式計算所得滲流量進(jìn)行比較，如圖8所示。在節(jié)理面粗糙程度較低時，單寬流量的迂曲度修正值與有限元計算結(jié)果較為接近，但是當(dāng)節(jié)理面粗糙程度較大時采用迂曲度修正的立方定律公式高估了單寬流量的值，在JRC值為18～20時約高估了19.4%。通過裂隙的單寬流量整體上隨著節(jié)理面粗糙程度的增加而降低。

圖8 數(shù)值計算結(jié)果與迂曲度修正值對比(隙寬3 mm)Fig.8 Comparison between the numerical results and the modified value of tortuosity (crack width 3 mm)

3.5 JRC值及裂隙寬度對滲流的影響

取水力坡度為0.01，對不同JRC值的裂隙取不同的隙寬分別進(jìn)行有限元計算，得到通過裂隙的流體平均流速如圖9所示。裂隙寬度為1 mm時，和光滑平行板裂隙內(nèi)流體的平均流速相比，JRC值為18～20的裂隙內(nèi)流體的平均流速下降了27.66%，當(dāng)裂隙寬度為3 mm和5 mm時，這一值分別下降了21.38%和19.71%。由此可見，隨著裂隙隙寬的增大，裂隙表面的粗糙程度對滲流的影響減小。當(dāng)裂隙寬度由1 mm分別變化到3 mm和5 mm時，光滑平行板裂隙內(nèi)的流體平均流速分別增加了9倍和25倍，而JRC值為18～20的裂隙內(nèi)流體的平均流速分別增加了9.64倍和27.74倍。由此可見，JRC值越大，裂隙隙寬變化對滲流的影響越顯著。這是由于粗糙度即裂隙表面的起伏作用范圍有限，隨著裂隙隙寬的逐漸增大，粗糙度對滲流的影響逐漸減小。

圖9 不同寬度裂隙內(nèi)流體平均流速對比Fig.9 Comparison of average velocity of fluid in different width fractures

不同JRC值及隙寬下裂隙內(nèi)流體最大流速對比如圖10所示。相同隙寬的裂隙，流體最大流速并沒有隨著裂隙面JRC值的變化呈現(xiàn)出規(guī)律性的變化趨勢。流體最大流速與裂隙面的局部粗糙度以及粗糙作用范圍有關(guān)，在裂隙寬度相同的情況下，裂隙面局部曲折程度越大流體最大流速越大。最大流速反映了裂隙內(nèi)流體流動速度的不均勻性，裂隙寬度越小節(jié)理面粗糙程度對最大流速的影響越大。當(dāng)裂隙寬度為5 mm時，各JRC值的裂隙中流體最大流速差別已較小，這說明節(jié)理面JRC值對裂隙內(nèi)滲流的影響與節(jié)理面最大起伏和裂隙寬度的比值有關(guān)，當(dāng)這一比值較小時，節(jié)理面JRC值對裂隙內(nèi)滲流的影響也較小。

圖10 不同寬度裂隙內(nèi)流體最大流速對比Fig.10 Comparison of relative maximum velocity of fluid in different width fractures

3.6 裂隙兩表面相互錯動對滲流的影響

選取JRC值為18～20、初始隙寬5 mm的一條裂隙，水力坡度值取為0.01，計算裂隙兩表面相對錯動不同距離時流體的滲流速度。裂隙兩表面相互錯動5 mm時的細(xì)觀流速分布如圖11所示。在隙寬最小處，流體流速最大，這是由于通過該裂隙的流體是穩(wěn)定流，各時刻通過各橫截面的流量相等，在隙寬最窄處流體的滲流斷面面積最小，因而流速最大。在裂隙中靠近裂隙面的位置流體流速最小，裂隙中央流體流速最大。

圖11 裂隙兩表面相對錯動5 mm時的流速分布Fig.11 Distribution of fluid velocity with relative displacement of 5 mm between two surfaces of fracture

裂隙兩表面不同錯動距離下滲流速度如圖12所示。裂隙內(nèi)滲流的平均流速隨著裂隙兩表面的相對錯動距離的增加而逐漸降低，錯動距離增加，裂隙最小寬度變窄，滲流量降低。通過裂隙的最大滲流速度并未隨著裂隙兩表面相對錯動距離的增大呈現(xiàn)規(guī)律性的變化。

圖12 裂隙兩表面相對錯動狀態(tài)下的滲流速度(JRC為18～20)Fig.12 Seepage velocity at relative dislocation of fracture surfaces (JRC is 18～20)

選取JRC值為18～20、初始隙寬5 mm的一條裂隙，計算裂隙兩表面相互錯動不同距離時不同水力坡度下的單寬流量值，計算結(jié)果如圖13所示。發(fā)現(xiàn)在裂隙兩表面相互錯動后，通過裂隙的單寬滲流量仍然與水力坡度呈線性關(guān)系，即符合達(dá)西定律。裂隙兩表面相互錯動距離越大，裂隙的滲透系數(shù)越小。裂隙兩表面相互錯動后，裂隙形狀變得復(fù)雜，流體的滲流路徑也變得更加曲折，因此裂隙的滲透系數(shù)減小。

圖13 裂隙兩表面不同錯位情況下單寬流量隨水力坡度的變化Fig.13 Variation of unit width seepage flow with hydraulic gradient under different dislocation of fracture surfaces

4 結(jié)論

使用ABAQUS自定義不可壓縮單元來求解不可壓縮黏性流體穩(wěn)態(tài)流動的Navier-Stokes方程，引入周期性邊界條件求得壓強(qiáng)梯度作用下表征體元的滲流速度，數(shù)值預(yù)測JRC值對滲流的影響，得出以下結(jié)論。

(1)僅選用一小段裂隙作為表征體元而忽略裂隙延伸長度計算得到的裂隙滲透性有一定的誤差，施加周期性邊界條件能有效地消除這種誤差。

(2)節(jié)理面粗糙度系數(shù)JRC值對裂隙的滲透性有影響，隨著裂隙寬度的增大，影響減小。當(dāng)裂隙寬度和水力坡度一定時，裂隙的滲透性整體趨勢上隨著JRC值的增大而減小。

(3)隨裂隙錯動距離的增大，裂隙的滲透性減小。