基于分形理論的巖石裂隙非線性滲流各向異性研究

王慧 蘇永軍 王鳳瑞

摘要：針對(duì)裂隙巖體非達(dá)西滲流問題，開展了不同粗糙裂隙非線性滲流特性的研究。將完整的方形巖塊采用劈裂法制成裂隙試件，用三維光學(xué)掃描系統(tǒng)量測(cè)裂隙表面形貌，并采用Kulatilake提出的自仿射分形維度方法計(jì)算表面的分形維度來表征粗糙表面的各向異性特征。根據(jù)裂隙流的曲折效應(yīng)，在Forchheimer定律基礎(chǔ)上，提出了水文彎曲度和表面分形冪律關(guān)系來表征非線性特點(diǎn)，并由此提出了新的粗糙裂隙非線性分形模型。對(duì)各組裂隙試件在2個(gè)方向（0°，90°）上進(jìn)行飽和滲流試驗(yàn)，發(fā)現(xiàn)流動(dòng)符合Forchheimer定律且存在各向異性特點(diǎn)，最后由分形模型，提出了區(qū)分達(dá)西流與Forchheimer流的新判據(jù)。

關(guān)鍵詞：巖石裂隙; 非線性滲流; 分形維度; Forchheimer定律; 各向異性; 臨界雷諾數(shù)

中圖法分類號(hào)：P641.2文獻(xiàn)標(biāo)志碼： ADOI：10.16232/j.cnki.1001-4179.2019.02.031

1研究背景

漫長(zhǎng)的地質(zhì)作用和人類活動(dòng)造成巖體被大量的斷層、裂隙切割，這些結(jié)構(gòu)面及其構(gòu)成的網(wǎng)絡(luò)成為地下水流動(dòng)的主要通道，并由此控制著巖體的滲透特性。在巖體水力學(xué)研究中，通常將結(jié)構(gòu)面概化為2塊光滑的平行板，通過理論和試驗(yàn)得到著名的立方定律。由于實(shí)際結(jié)構(gòu)面粗糙起伏、零星接觸或含有充填物等，許多學(xué)者據(jù)此提出各種修正模型。

在某些工程中，如河谷深厚覆蓋層建壩、低滲透油氣井開采、煤礦瓦斯突出等[1]都會(huì)出現(xiàn)高水力梯度現(xiàn)象，這時(shí)流體流動(dòng)的機(jī)制和規(guī)律將發(fā)生重大改變，用立方定律或其修正模型會(huì)造成較大的偏差。一般引用多孔介質(zhì)的非線性滲流模型Forchheimer定律[2]來描述這種滲流行為

▽P=AQ+BQ?2（1）

式中，▽P為單位滲流長(zhǎng)度的壓力差，Q為通過裂隙的流量，A和B分別為黏性系數(shù)和慣性系數(shù)。

Zimmerman等[3]通過試驗(yàn)和數(shù)值的方法，觀察到雷諾數(shù)Re>20時(shí)，粗糙裂隙的Forchheimer流現(xiàn)象;Zhang等[4]探討了不同圍壓下，粗糙裂隙線性和非線性流動(dòng)特性;Zhou等[5]利用不同圍壓的壓水試驗(yàn)，解釋了Forchheimer流系數(shù)A和B的物理意義，及內(nèi)部過渡機(jī)制，但對(duì)裂隙面粗糙性對(duì)非線性流動(dòng)的影響沒有詳細(xì)闡明。Chen等評(píng)價(jià)了Forchheimer判據(jù)方程的系數(shù)[6]。金毅等[7]從細(xì)觀層面上指出粗糙幾何對(duì)裂隙流的影響表現(xiàn)在三個(gè)方面：① 流體內(nèi)部的摩擦效應(yīng);② 裂隙面的曲折效應(yīng);③ 局部粗糙度效應(yīng)。Tsang等[8]認(rèn)為裂隙面的粗糙性會(huì)引發(fā)流動(dòng)的曲折性;肖維民等[9]引入曲折因子描述這種流動(dòng)的曲折現(xiàn)象。

自分形幾何被B.B.Mandelbrot提出以后，謝和平等[10]首先將其引入到裂隙粗糙度的描述，后來又用來描述巖體斷裂面滲流特征[11];Murata等[12]研究了分形參數(shù)對(duì)曲折效應(yīng)的影響;王剛等[13]提出了考慮分形特征的節(jié)理面滲流公式;Ju等[14]對(duì)不同分形維度粗糙單裂隙物理模型進(jìn)行水滲流試驗(yàn)，闡明了粗糙結(jié)構(gòu)對(duì)滲流的影響;Develi等[15]對(duì)7種人工張拉型裂隙面進(jìn)行飽和滲流試驗(yàn)，并用分形維度描述粗糙度，研究了粗糙度、各向異性和法向應(yīng)力對(duì)滲流特征的影響。

本文從曲折效應(yīng)對(duì)裂隙流非線性作用強(qiáng)度和方式入手，并結(jié)合裂隙自仿射分形特點(diǎn)，推導(dǎo)了粗糙裂隙非線性分形模型。并結(jié)合滲流試驗(yàn)，分析了Forchheimer流的規(guī)律和各向異性，并驗(yàn)證了新模型。

2粗糙裂隙非線性分形模型

Forchheimer流動(dòng)定律包含兩個(gè)部分：線性部分AQ和非線性部分BQ?2，體現(xiàn)出達(dá)西滲流和非達(dá)西滲流。在低流速階段，一般可用立方定律描述流量和壓力的線性關(guān)系，則黏性系數(shù)A可表示為

A=12μwe?3h（2）

式中，μ為流體的黏性系數(shù)，w為裂隙的寬度，eh為水力開度，是在初始階段根據(jù)立方定律反求而得，eh=（12μQ/w/▽P）1/3。

慣性系數(shù)B體現(xiàn)了滲流曲線偏離線性段的程度。Schrauf等[2]通過量綱分析，提出

B=bDρe?3hw?2（3）

式中，ρ為流體密度，bD是與裂隙面結(jié)構(gòu)特征有關(guān)的參數(shù)。Chen等[6]用絕對(duì)粗糙度描述裂隙面幾何特征，提出的雙參數(shù)模型為

bD=a（ξ2eh）?b（4）

式中，a和b為擬合參數(shù)。但僅僅用絕對(duì)粗糙度ξ表征裂隙面粗糙程度的影響，不能反映流動(dòng)的曲折性和各向異性。郭建春等[16]通過數(shù)值方法研究了裂隙流的非達(dá)西效應(yīng)，結(jié)果表明表面越粗糙，隙寬分布越復(fù)雜，流動(dòng)時(shí)流態(tài)將越不穩(wěn)定，曲折效應(yīng)將增加，在高流速下會(huì)出現(xiàn)渦流和回流現(xiàn)象，增大慣性阻力。為了表征曲折性的影響，提出如下冪律關(guān)系

bD=cτ?aτ?bs（5）

式中，a，b和c為擬合參數(shù);τ為水文彎曲度，定義為流體實(shí)際滲流路徑長(zhǎng)度Lt與沿裂隙水頭壓力梯度方向的水平長(zhǎng)度Lc之比。τs為裂隙面曲折率，定義為裂隙面表面積As與投影面積Ac之比。對(duì)于具有分形特征的裂隙，測(cè)度F（δ）與測(cè)量尺度δ存在如下關(guān)系

F（δ）=F0δ（6）

式中，與分形維數(shù)D有關(guān)，即=f（D），D∈（1，3）;F0為直觀測(cè)度。邊長(zhǎng)為δ1的正方形網(wǎng)格覆蓋的裂隙表面總的面積AS與δ1存在如下關(guān)系

AS=F1δ2-DS1（7）

式中，DS為裂隙面面積分形維度，對(duì)于正方形裂隙，當(dāng)δ1=Lc時(shí)，將（7）式代入AS（Lc）=Ac得：

F1=ADS2c（8）

由τs的定義，有

τs=AS（δ1）Ac=（δ?21Ac）2-DS（9）

同理，對(duì)于滲流路徑長(zhǎng)度Lt同觀測(cè)尺度δ2存在如下分形關(guān)系：

Lt（δ2）=F2δ1-DT2（10）

式中，DT為滲流路徑長(zhǎng)度分形維度。

當(dāng)δ2=Lc時(shí)，將（10）式代入Lt（Lc）=Lc得

F2=LDTc（11）

由τ定義得

τ=Lt（δ2）Lc=（δ2Lc）1-DT（12）

分形創(chuàng)始人Mandelbrot建議表面分形維數(shù)可用表面某一剖線的維數(shù)加1的方法計(jì)算，故對(duì)于剖面線長(zhǎng)度分形維數(shù)DL與面積分形維度DS存在如下關(guān)系

DS=DL+1（13）

金毅等[7]研究表明，DL≈DT，故將其與式（9）和式（12）代入式（5）中，得

bD=c（δ2Lc）a（1-Dt）（δ?21Ac）b（1-Dt）（14）

取δ1和δ2的測(cè)度為eh，同時(shí)有Ac=L?2c，代入式（14）得

bD=c（ehLc）（a+b）（1-DL）（15）

重新整理式（15），得到B新的表達(dá)式為

B=a（ehLc）b（1-DL）ρe?3hw?2（16）

式中，a和b為待定常數(shù)，可由裂隙滲透試驗(yàn)確定。計(jì)算時(shí)，首先獲取裂隙表面數(shù)據(jù)，計(jì)算分形維度，然后通過裂隙滲透試驗(yàn)得到流量與壓力梯度的曲線，擬合得到A和B參數(shù)，并進(jìn)一步獲得a和b。

3裂隙面測(cè)量與分形計(jì)算

3.1裂隙面制備

為了研究裂隙面粗糙度對(duì)流體流動(dòng)的影響，需要對(duì)不同粗糙程度的巖體裂隙面進(jìn)行飽和滲流試驗(yàn)。而天然裂隙面較難獲取，因此，本次采用人工張拉裂隙試件研究裂隙中流體流動(dòng)特征。將選自湖北省境內(nèi)某采石場(chǎng)的天然花崗巖在室內(nèi)加工成尺寸為150 mm×150 mm×150 mm方形試塊，然后用巴西劈裂試驗(yàn)方法將試塊劈開，制成人工張拉性節(jié)理試件。最終制備成5組不同粗糙程度的裂隙面（F1，F(xiàn)2，F(xiàn)3，F(xiàn)4和F5）。

3.2裂隙表面形貌測(cè)量

試驗(yàn)采用便攜式3D CaMega光學(xué)三維掃描系統(tǒng)對(duì)裂隙表面三維形貌進(jìn)行量測(cè)。3D CaMega光學(xué)三維掃描系統(tǒng)是一臺(tái)非接觸高精度光學(xué)掃描儀（見圖1）。該系統(tǒng)主要由掃描控制單元、掃描鏡頭、轉(zhuǎn)臺(tái)和三腳架組成;掃描鏡頭安放在三角架上，三腳架可自由伸縮旋轉(zhuǎn)，方便調(diào)整位置，其他部分通過USB方式連接。系統(tǒng)通過獲取投射到物體表面的可見光光柵圖像，按照條紋曲率變化形狀利用相位法和三角法精確確定每一點(diǎn)空間坐標(biāo)（X，Y，Z），形成三維點(diǎn)云數(shù)據(jù)，具有快速、精度高（測(cè)量精度25μm）、非接觸測(cè)量等的特點(diǎn)。

圖13D光學(xué)掃描系統(tǒng)示意Fig.1Stereotopometric scanning system

在實(shí)際測(cè)量過程中，測(cè)量環(huán)境（光線、粉塵等）和測(cè)量方法等對(duì)三維形貌數(shù)據(jù)的精度都會(huì)產(chǎn)生影響，故在獲取裂隙面三維數(shù)據(jù)后，運(yùn)用自帶軟件CloudForm對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行降噪、去除無關(guān)點(diǎn)、濾除數(shù)據(jù)波紋以及補(bǔ)面等處理。同時(shí)，原每個(gè)裂隙面均由幾十萬個(gè)離散點(diǎn)構(gòu)成，平均點(diǎn)距0.025 mm，數(shù)據(jù)量龐大且排列雜亂無序。為了便于分析數(shù)據(jù)，用matlab編寫程序，對(duì)測(cè)量點(diǎn)重新刪減和排序，新得到的裂隙面共計(jì)22 801個(gè)點(diǎn)，平均點(diǎn)距1 mm，如圖2所示，測(cè)得的裂隙表面形貌參數(shù)見表1。

3.3裂隙分形描述與計(jì)算

粗糙裂隙表面分形維度計(jì)算方法多種。Clarke等[17]首次提出了以空間表面面積為測(cè)度，以正方形網(wǎng)格為測(cè)量尺度的三角形棱柱表面積法;謝和平等[18]提出投影覆蓋法。這兩種方法屬于碼尺法（driver method）范疇，計(jì)算得到的維數(shù)是相似分形維數(shù)，而不是真正意義上的幾何分形維數(shù)。后來周宏偉等[19]提出計(jì)算三維表面分形維度的計(jì)盒維數(shù)法（box-counting method）：立方體覆蓋法和改進(jìn)的立方體覆蓋法。以上計(jì)算理論都是基于統(tǒng)計(jì)意義上的自相似性。事實(shí)上，Brown、Kulatilake等[20]學(xué)者認(rèn)為巖石粗糙裂隙面符合自仿射模型特征。

對(duì)于自仿射模型，Kulatilake[21]提出以二維剖面線Z（x）變異函數(shù)值2γ（x，h）為測(cè)度，間距h為測(cè)量尺度的變異函數(shù)法（Variogram method）。

第一步：生成不同方向二維輪廓線數(shù)據(jù)。首先，將獲得的裂隙表面數(shù)據(jù)劃分成1 mm×1 mm網(wǎng)格，選取0°方向，其他方向逆時(shí)針按照θ=15°×k（k=1，2，…，11）選擇，見圖3。然后，計(jì)算一條方向線上高度數(shù)據(jù)Z（x，y），如對(duì)15°方向線P上坐標(biāo)Q（x，y）表面高度數(shù)據(jù)ZQ對(duì)該坐標(biāo)半徑1 mm范圍內(nèi)的高度按照（6）式計(jì)算，循環(huán)獲得該條輪廓線數(shù)據(jù)。再以相同的方法獲得相隔10 mm的下一組輪廓線數(shù)據(jù)。最后反復(fù)循環(huán)計(jì)算獲得各個(gè)方向上輪廓線數(shù)據(jù)。

Z（x，y）=ni=1Zirini=11ri（17）

式中，Zi為與Q點(diǎn)相距1 mm半徑范圍內(nèi)第i個(gè)點(diǎn)的高度，ri為i點(diǎn)到Q點(diǎn)距離。

第二步：計(jì)算不同方向二維分形維度。對(duì)每個(gè)方向，考慮邊界影響，剔除兩邊的二維輪廓線;然后計(jì)算各輪廓線的變異函數(shù)值2γ（x，h）

2γ（x，h）=1NNi=1[Z（xi）-Z（xi+h）]?2（18）

式中，N為輪廓線上所有間距h的數(shù)據(jù)點(diǎn)數(shù)。對(duì)于自仿射模型，存在如下關(guān)系

2γ（x，h）=KVh2（2-D）（19）

式中，KV為比例系數(shù)，D為二維輪廓線分形維度。對(duì)每個(gè)輪廓線計(jì)算至少7個(gè)不同間距h下的變異函數(shù)，并由式（19）擬合得到分形維度D，對(duì)一個(gè)方向下所有輪廓線的分形維度求平均，得到該方向下的分形維度D。

為了使裂隙面粗糙度各向異性特征更加直觀，圖4給出了5個(gè)不同裂隙面（F1～F5）分形維度的玫瑰花圖。各圖呈現(xiàn)出近似“橢圓”對(duì)稱狀，但分形維度隨方向隨機(jī)分布，裂隙面具有不規(guī)則的粗糙度各向異性結(jié)構(gòu)特征。通過比較，F(xiàn)4裂隙面分形維度較大，其中90°方向最大，分形維度為1.60，由此可知，F(xiàn)4具有較大的粗糙度。

4粗糙裂隙非線性流動(dòng)特性

4.1粗糙裂隙飽和滲流試驗(yàn)

試驗(yàn)中，裂隙采用上節(jié)提到的5組巖石裂隙面作為原型，用硅膠拷貝其表面的形貌，然后基于硅膠模型，用環(huán)氧樹脂AB膠復(fù)制其表面形態(tài)。環(huán)氧樹脂AB膠是一種液型、雙組份硬性膠，固化以后形成高透明、高強(qiáng)度的硬性體，在試驗(yàn)過程中可實(shí)時(shí)觀察水流流動(dòng)情況。

采用自行設(shè)計(jì)的裝置（見圖5）對(duì)5組裂隙進(jìn)行飽和滲流試驗(yàn)，每組裂隙進(jìn)行2種流向（0°，90°）測(cè)試。測(cè)試過程記錄不同流量Q下的壓力差△P，試驗(yàn)的流量范圍在0～100 mL/s之間。在上下裂隙表面匹配后，利用Kulatilake[21]提出的非接觸方法確定每個(gè)裂隙開度值，F(xiàn)1～F5開度值分別為1.58，2.20，1.50，2.45 mm和2.10 mm。

4.2非線性流動(dòng)規(guī)律

圖6給出了0°方向各組裂隙流量和壓力梯度的關(guān)系。當(dāng)流量較小（Q<10 mL/s）時(shí)，黏滯作用占主導(dǎo)地位，壓力梯度隨著流量呈線性增加;當(dāng)流量較大（Q>10 mL/s）時(shí)，慣性作用增強(qiáng)，流動(dòng)從線性階段過渡至非線性階段，壓力梯度將隨著流量增加呈非線性遞增。為了表達(dá)這種關(guān)系，結(jié)合最小二乘法，用Forchheimer公式擬合試驗(yàn)數(shù)據(jù)。同時(shí)擬合90°方向滲流數(shù)據(jù)，擬合結(jié)果見表2。各組裂隙的擬合精度都較高，說明Forchheimer公式能夠定量描述非線性流動(dòng)規(guī)律，這與Zimmerman等[3]研究結(jié)論是一致的。

為了分析流動(dòng)的各向異性特征，提出了一個(gè)新的描述各向異性的參數(shù)：各向異性度▽，定義為90°向與0°向壓力梯度之差與90°向壓力梯度的比值。

▽=▽P90°-▽P0°▽P90°（20）

式中，▽P90°為90°向壓力梯度，▽P0°為0°向壓力梯度。圖7給出了各向異性度隨流量的變化曲線，從圖可以看出：總體上，隨著流量的增加，各組裂隙各向異性度在緩慢減小，但變化不明顯;各組裂隙各向異性度大多在0附近“波動(dòng)”，且各值有一定的差別，說明裂隙流存在各向異性特征，且與裂隙形貌和開度分布有關(guān)。

4.3非線性分形模型的驗(yàn)證

4.3.1系數(shù)a和b的確定

為了求解待定常數(shù)a和b，首先獲取F1、F2和F3裂隙形貌數(shù)據(jù)，并根據(jù)上節(jié)方法計(jì)算分形維度;然后利用滲流試驗(yàn)結(jié)果（見表2），根據(jù)非線性擬合Levenberg-Marquardt法對(duì)式（16）進(jìn)行擬合，得到a和b。擬合結(jié)果見圖8，得到a=0.246，b=-0.964，代入式（1），（2）和（16），得

▽P=12μwe?3hQ+0.246（ehL）-0.964（1-D）ρe?3hw?2Q?2（21）

圖8B與eh/L的擬合關(guān)系Fig.8Fitting relationship between B and eh/L

4.3.2非線性分形模型與其他模型的對(duì)比

為了驗(yàn)證模型的正確性，將非線性分形模型與滲流試驗(yàn)數(shù)據(jù)及Chen的雙參數(shù)模型[6]對(duì)比。Chen的公式如下

▽P=12μwe?3hQ+0.0138 7ρξ0.666w?2e3.666hQ?2（22）

針對(duì)裂隙F4和F5，分別按式（21）和（22）計(jì)算壓力梯度，并將結(jié)果與實(shí)驗(yàn)值進(jìn)行比較，結(jié)果見圖9。由圖可知：采用非線性分形模型計(jì)算結(jié)果與實(shí)測(cè)值較為接近，相對(duì)誤差大多數(shù)在20%以內(nèi)，表明非線性分形模型能較好地描述裂隙介質(zhì)中非線性滲流情況。而Chen的公式計(jì)算誤差較大，原因在于僅僅用絕對(duì)粗糙度難以量化表面形貌對(duì)裂隙流的影響。

5分析與討論

Forchheimer定律在裂隙介質(zhì)非線性滲流中有著廣泛的應(yīng)用價(jià)值，但對(duì)于線性流過渡到非線性流的機(jī)制還需進(jìn)一步討論。

Zimmerman等[3]首次提出裂隙內(nèi)流體流動(dòng)存在線性流到非線性流的轉(zhuǎn)變過程，并用非線性達(dá)西效應(yīng)因子區(qū)分兩種流態(tài)，認(rèn)為=0.1。

=BQ?2AQ+BQ?2（23）

據(jù)此，眾多學(xué)者[22]提出用臨界雷諾數(shù)Re0描述這種過渡機(jī)制。對(duì)于裂隙介質(zhì)，雷諾數(shù)[5]可用下式表達(dá)

Re=ρQμw（24）

將式（2）、（16）和式（23）代入式（24），并令=0.1，得到：

Rec=5.42（ehLc）-0.964（DL-1）（25）

式（25）表明，臨界雷諾數(shù)與水力開度、裂隙面的分形維度及流動(dòng)的方向性密切相關(guān)，水力開度越小，裂隙面越粗糙，臨界雷諾數(shù)也越小，裂隙流越容易發(fā)生非線性流動(dòng)。計(jì)算得到的臨界雷諾數(shù)見表3，大小在30～60之間，遠(yuǎn)小于王媛等[23]認(rèn)為的2300。進(jìn)一步表明，裂隙隙寬分布不一，表面的粗糙所造成流動(dòng)的渦流及回流效應(yīng)，使得流動(dòng)變得曲折，慣性阻力增加，從而導(dǎo)致在低的雷諾數(shù)時(shí)轉(zhuǎn)變成非線性流。

6結(jié) 論

（1） 根據(jù)裂隙流具有曲折效應(yīng)的特點(diǎn)，并用水文彎曲度和表面曲折率來描述，結(jié)合裂隙分形特征，提出了粗糙裂隙非線性滲流新公式。

（2） 運(yùn)用了3D CaMega光學(xué)三維掃描系統(tǒng)獲取了裂隙面點(diǎn)云數(shù)據(jù)，并采用了Kulatilake提出的自仿射模型分形維度計(jì)算方法，分析了裂隙面粗糙性的各向異性特征。

（3） 對(duì)5種不同粗糙裂隙分別從0°和90°方向進(jìn)行滲流試驗(yàn)。試驗(yàn)結(jié)果表明，在高流速下，裂隙流符合Forchheimer定律，且具有各向異性特點(diǎn)。并對(duì)比了新模型、Chen的雙參數(shù)模型和試驗(yàn)結(jié)果。結(jié)果表明，新模型與實(shí)驗(yàn)結(jié)果更接近。

（4） 根據(jù)新模型，提出了區(qū)分達(dá)西流與Forchheimer流的判據(jù)，即臨界雷諾數(shù)的表達(dá)式。分析表明，臨界雷諾數(shù)與水力開度、裂隙面的粗糙度及流動(dòng)的方向性密切相關(guān)。

（5） 本文提出的模型沒有考慮相反方向滲流各向異性特征，而這有賴于新的粗糙度參數(shù)，這將是今后研究的重點(diǎn)。

參考文獻(xiàn)：

[1]許凱，雷學(xué)文，孟慶山，等.非達(dá)西滲流系數(shù)研究[J].巖石力學(xué)與工程學(xué)報(bào)，2012，31（1）：164-170.

[2]Schrauf T W，Evans D D.Laboratory studies of gas flow though a single natural fracture[J].Water Resources Research，1986，22（7）：1038-1050.

[3]Zimmerman R W，Al-Yaarubi A，Pain C C，et al.Nonlinear regimes of fluid flow in rock fractures[J].International Journal of Rock Mechanics and Ming Sciences，2004，41（3）：163-169.

[4]Zhang Z， Nemcik J.Fluid flow regimes and nonlinear flow characteristics in deformable rock fractures[J].Journal of Hydrology，2013（477）：139-151.

[5]Zhou J Q，Hu S H，F(xiàn)ang S，et al.Nonlinear flow behavior at low Reynolds numbers though rough-walled fractures subjected to normal compressive loading[J].International Journal of Rock Mechanics and Ming Sciences，2015（80）：202-218.

[6]Chen Y F，Zhou J Q，Hu S H，et al.Evaluation of Forchheimer equation coefficients for non-Darcy flow in deformable rough-walled fractures[J].Journal of Hydrology，2015（529）：993-1006.

[7]金毅，鄭軍領(lǐng)，董佳斌，等.自仿射粗糙割理中流體滲流的分形定律[J].科學(xué)通報(bào)，2015，60（21）：2036-2047.

[8]Tsang Y W.The effect of tortuosity on fluid flow through a single fracture[J].Lubric Technology，1984，20（9）：1209-1215.

[9]肖維民，夏才初，王偉，等.考慮曲折效應(yīng)的粗糙節(jié)理滲流計(jì)算公式研究[J].巖石力學(xué)與工程學(xué)報(bào)，2011，30（12）：2416-2425.

[10]謝和平.分形幾何及其在巖土力學(xué)中的應(yīng)用[J].巖土工程學(xué)報(bào)，1992，14（1）：14-24.

[11]謝和平，周宏偉.巖體斷裂面滲流特性的分形研究[J].煤炭學(xué)報(bào)，1998，23（6）：585-589.

[12]Murata S，Saito T.Estimation of tortuosity of fluid flow through a single fracture[J].Journal of Canadian Petroleum Technology，2003，42（12）：39-45.

[13]王剛，黃娜，蔣宇靜，等.考慮分形特征的節(jié)理面滲流計(jì)算模型[J].巖石力學(xué)與工程學(xué)報(bào)，2014，33（S2）：3397-3405.

[14]Ju Y，Zhang Q G，Yang Y M，et al.An experimental investigation on the mechanism of fluid flow through single rough fracture of rock[J].Science China Technological Sciences，2013（56）：2070-2080.

[15]Develi K，Babadagli T.Experimental and visual analysis of single-phase flow through rough fracture replicas[J].International Journal of Rock Mechanics and Ming Sciences，2015（73）：139-155.

[16]郭建春，莊園，劉超.考慮非達(dá)西效應(yīng)酸蝕裂縫流場(chǎng)數(shù)值模擬[J].巖土力學(xué)，2015，36（11）：3315-3321.

[17]Clarke K C.Computation of the fractal dimension of topographic surfaces using the triangular prism surface area method[J].Computers and Geosciences，1986，12（5）：713-722.

[18]Xie H P，Wang J A.Direct fractal measurement of fractal surfaces[J].International Journal Solids and Structures，1999（36）：3073-3084.

[19]周宏偉，謝和平，Kwansniewski M A.粗糙表面分維計(jì)算的立方體覆蓋法[J].摩擦學(xué)學(xué)報(bào)，2000，20（6）：455-459.

[20]Kulatilake P H S W，Balasingam P，Park J，Morgan R.Natural rock joint roughness quantification through fractural techniques[J].Geotechnical and Geological Engineering，2006（24）：1181-1202.

[21]Kulatilake P H S W，Park J，Balasingam P，et al.R.Quantification of aperture and relations between aperture， normal stress and fluid flow for natural single rock fractures[J].Geotechnical and Geological Engineering， 2008（26）：269-281.

[22]Javadi M，Sharifzadeh M，Shahriar K，et al.Crital Reynolds number for nonlinear flow through rough-walled fractures：The role of shear processes[J].Water Resources Research，2014，50（2）：1789-1804.

[23]王媛，顧智剛，倪小東.光滑裂隙高流速非達(dá)西滲流運(yùn)動(dòng)規(guī)律的試驗(yàn)研究[J].巖石力學(xué)與工程學(xué)報(bào)，2010，29（7）：1404-1407.

引用本文：王慧，蘇永軍，王鳳瑞.基于分形理論的巖石裂隙非線性滲流各向異性研究[J].人民長(zhǎng)江，2019，50（2）：174-180.

Study on non-linear flow anisotropy behavior in rough rock fractures based on fractal theory

WANG Hui?，SU Yongjun?，WANG Fengrui?2

（1.Hebei University of Water Resources and Electric Engineering， Cangzhou 061001，China;2. Hebei Cangzhou Hydrology and Water Resources Survey Bureau，Cangzhou 061000， China）

Abstract： Aiming at the issue of Non-Darcy seepage flow in fracture， the seepage characteristics in different rough fractures were studied . The test speciment of fracture rock mass were obtained by splitting the integrated cubic rock mass， and the fracture surface topographies of rock mass were measured by stereotopometric scanning system. Surface roughness anisotropy was quantified by self-affine fractal dimension D using measured data. According to tortuosity effect of flow in fracture， the non-linear flow was characterized by a power law relationship of hydraulic tortuosity and fractal dimension based on Forchheimer law. Fracture seepage tests were conducted from two incidence directions （0°and 90°）. The results showed that non-linear flow with anisotropy could be well described by Forchheimer law. Finally， a new formula was proposed， which can distinguish the linear flow and Forchheimer flow.

Key words：fracture; non-linear flow; fractal dimension; Forchheimer law; anisotropy; critical Reynolds number