李 博，汪佳飛，劉日成，伍法權(quán)

(1.紹興文理學(xué)院 土木工程學(xué)院，浙江 紹興 312000；2.中國礦業(yè)大學(xué) 深部巖土力學(xué)與地下工程國家重點(diǎn)實(shí)驗(yàn)室，江蘇 徐州 221116)

探索巖體中流體的流動機(jī)理對解決隧道開挖和地下工程施工過程中的突涌水問題具有重要意義[1]。工程擾動引起的應(yīng)力、水壓變化會改變裂隙內(nèi)空腔的大小和幾何形態(tài)，進(jìn)而導(dǎo)致滲流特性的變化。這一變化在裂隙發(fā)生剪切錯位時(shí)顯得尤為突出[2]。在裂隙發(fā)育的富水巖體中，工程擾動容易使流體進(jìn)入非線性流動狀態(tài)，給滲流特性的精確評價(jià)帶來困難[3]。因此，通過基礎(chǔ)試驗(yàn)和數(shù)值模擬深入研究裂隙在彈性和彈塑性變形條件下的非線性滲流特性具有重要意義。

裂隙巖體一般包含有大量的單裂隙，研究單裂隙的變形與滲透特性是理解裂隙巖體水力耦合特性的先決條件。早期的研究將巖石裂隙簡化為2個(gè)光滑的平行板，利用Navier-Stocks方程可推導(dǎo)出表示流量與開度之間關(guān)系的立方定律。實(shí)際的天然裂隙由粗糙表面組成，開度e在空間上不是一個(gè)恒定值。因此，大量學(xué)者利用力學(xué)開度em、等效水力開度eh及表征裂隙面粗糙特征的參數(shù)對立方定律進(jìn)行修正[4]，在一定程度上拓展了立方定律的適用范圍。

對于單裂隙而言，施加于裂隙上的法向應(yīng)力所引起的變形是影響裂隙過流能力的主要原因。眾多學(xué)者對法向應(yīng)力作用下粗糙巖石裂隙的變形和滲流特性進(jìn)行了廣泛的試驗(yàn)和理論研究。Wang等[5]研究表明，巖石在受正應(yīng)力作用時(shí)，裂隙的閉合曲線為受裂隙表面特征控制的函數(shù)。Iwai[6]通過試驗(yàn)發(fā)現(xiàn)，裂隙面粗糙度對滲流特性的影響主要是通過改變其接觸面積體現(xiàn)。Kranzz等[7]通過室內(nèi)試驗(yàn)指出裂隙滲透系數(shù)與法向應(yīng)力呈冪函數(shù)遞減關(guān)系。張文泉等[8]根據(jù)采動過程中裂隙巖體的應(yīng)力變化，探討在恒定法向荷載和恒定法向剛度條件下裂隙巖體粗糙度和滲透水壓對試樣滲透系數(shù)的影響規(guī)律，發(fā)現(xiàn)粗糙度與巖體的剛度成反比，與水力開度最終穩(wěn)定值成正相關(guān)。

影響巖石裂隙水力耦合特性的另一個(gè)重要因素是剪切過程。一方面，剪切錯位造成裂隙內(nèi)部接觸狀態(tài)變化，直接影響裂隙在法向應(yīng)力作用下的變形特性；另一方面，剪切造成內(nèi)部空腔結(jié)構(gòu)變化，導(dǎo)致裂隙滲透系數(shù)及非線性滲流特性變化[9]。學(xué)者們通過剪切-滲流耦合試驗(yàn)，已初步掌握滲透系數(shù)在剪切過程中隨剪脹的三段式演化規(guī)律[10]。通過求解Navier-Stocks（NS）方程，發(fā)現(xiàn)了表征非線性滲流的臨界水力梯度和臨界雷諾數(shù)在剪切過程中的演化特征[11]。結(jié)果表明，隨著剪切位移的增大，法向位移不斷增大，臨界水力梯度和臨界雷諾數(shù)均逐漸增大，且增大速率均為先增大后減少。

除了流體本身的性質(zhì)之外，裂隙中的接觸面和粗糙度，即裂隙內(nèi)的空腔結(jié)構(gòu)特征是導(dǎo)致非線性滲流的重要原因。流速較小時(shí)，流體一般為達(dá)西流，可以用達(dá)西定律計(jì)算裂隙的滲透系數(shù)；但隨著流速的增大，裂隙兩端的水頭差與流量之間逐漸演變?yōu)榉蔷€性關(guān)系，這一轉(zhuǎn)換受到裂隙內(nèi)部空腔結(jié)構(gòu)特征的控制。Bear[12]提出用Forchheimer定律表示該非線性關(guān)系，已得到大量試驗(yàn)和數(shù)值模擬結(jié)果的驗(yàn)證；該方程中的慣性系數(shù)B受到諸多因素的影響，難以定量確定，在實(shí)際工程應(yīng)用中還存在較大困難。為此，Chen等[13]研究了可變形粗糙裂隙中圍壓對Forchheimer定律中系數(shù)A和B的影響；Zhou等[14]整理了大量室內(nèi)和現(xiàn)場滲流試驗(yàn)結(jié)果和孔隙尺度滲流模型計(jì)算結(jié)果，定量分析了Forchheimer定律中的黏性滲透系數(shù)kv和慣性滲透系數(shù)ki，發(fā)現(xiàn)兩種滲透系數(shù)存在ki=β·kvα/2的關(guān)系；Liu等[15]通過數(shù)值方法生成大量粗糙裂隙表面，并模擬了裂隙面剪切滲流過程，認(rèn)為隨著剪切的進(jìn)行，系數(shù)A和B具有相同的變化規(guī)律，但系數(shù)B的數(shù)值比A的高約6個(gè)數(shù)量級；Zhang等[16]發(fā)現(xiàn)系數(shù)B隨著微裂縫開度的增大而減??；Li[17]和Zhang[18]等都發(fā)現(xiàn)了當(dāng)開度場固定時(shí)，增加剪切位移和粗糙度可以提高Forchheimer定律中二次項(xiàng)的比例。同一裂隙在不同剪切位移和受力環(huán)境下的幾何形態(tài)不同，從而改變了裂隙交叉口的幾何特征，Liu等[19]發(fā)現(xiàn)在復(fù)雜裂隙網(wǎng)絡(luò)中的裂縫通道長度和寬度是影響比流體分流比的主要因素。

以往研究一般分別考慮法向應(yīng)力、剪切過程和非線性滲流，尚未綜合建立關(guān)鍵控制參數(shù)，例如應(yīng)力、位移、開度、慣性系數(shù)等之間的關(guān)系，且相關(guān)模型也較少得到試驗(yàn)驗(yàn)證。為定量研究法向應(yīng)力和剪切錯位對裂隙非線性滲流特性的影響規(guī)律，針對含單裂隙花崗巖試樣開展了應(yīng)力-滲流耦合試驗(yàn)，基于自主開發(fā)的裂隙受壓變形模擬程序計(jì)算裂隙變形，獲取裂隙內(nèi)部空腔幾何數(shù)據(jù)；然后，利用COMSOL軟件模擬裂隙內(nèi)的滲流，分析同一裂隙在不同剪切錯動條件下的變形破壞與非線性滲流特性，找出各控制參數(shù)之間的數(shù)學(xué)關(guān)系。

1 試驗(yàn)方法

試驗(yàn)中主要采用了夾持器、TELEDYNE ISCO 柱塞泵、VR-3100 3D 輪廓測量儀、高精度差壓計(jì)等儀器。柱塞泵的流量量程為0.001至200 mL/min，分辨率為0.001 mL/min；差壓計(jì)的量程為10 Pa至100 kPa，分辨率為10 Pa。試驗(yàn)選取3個(gè)采自同一地點(diǎn)的花崗巖標(biāo)準(zhǔn)試樣（G-1、G-2、G-3），直徑均為50 mm，高度均為100 mm，將試樣放置于具有楔形壓頭的萬能試驗(yàn)機(jī)上，通過加壓沿軸線劈裂制成裂隙，見圖1。

圖1 花崗巖試樣制備和裂隙表面形貌Fig. 1 Preparation of granite samples and their surface morphology

由于夾持器內(nèi)部橡膠圈直徑有限，表面起伏大的裂隙在錯動后無法放置，所以試驗(yàn)選取裂隙表面起伏較小的試樣（G-1）作為后期滲流試驗(yàn)的對象和數(shù)值對象，G-2和G-3掃描后作為數(shù)值模擬對象。為更精確地控制試樣的剪切位移量，根據(jù)試樣的幾何特征，利用3D打印機(jī)打印了2、4、6、8 mm 4種厚度的高強(qiáng)度樹脂墊片各2片。這些墊片底部設(shè)置了鋸齒狀結(jié)構(gòu)，不影響水的流動；將墊片用膠水粘貼于試樣端面后，為橡膠套提供支撐。

圖2為水力耦合試驗(yàn)裝置的示意圖。如圖2所示，將試樣置于夾持器中，夾持器的入口端連接柱塞泵C和差壓計(jì)入口端，出口端連接電子秤上的廢液盒和差壓計(jì)出口端；圍壓孔連接柱塞泵A，軸壓孔連接柱塞泵B。另外，準(zhǔn)備3個(gè)標(biāo)準(zhǔn)試樣在MTS815試驗(yàn)機(jī)上開展單軸壓縮試驗(yàn)，測得該種花崗巖的平均單軸抗壓強(qiáng)度UCS=137.80 MPa，彈性模量E=41.20 GPa，泊松比ν=0.29。由于研究主要關(guān)注不同剪切錯位下法向應(yīng)力對滲透特性的影響效果，所以未采用直接剪切的形式，而是在一個(gè)剪切位移測試結(jié)束后取出試樣，掃描表面；然后，人工錯位到下一個(gè)剪切位移，再裝入試樣進(jìn)行法向加載和滲流測試，以保證所有的表面損傷皆源自法向加載。

圖2 水力耦合試驗(yàn)裝置示意圖Fig. 2 Schematic diagram of the coupled stress-flow testing apparatus

首先，在開始滲流試驗(yàn)前先施加0.50 MPa的圍壓和軸壓，從而保證試樣的密閉性；然后，緩慢注水以排出內(nèi)部氣體使各個(gè)導(dǎo)管與試樣內(nèi)部完全充滿水；隨后，關(guān)閉柱塞泵C，待差壓計(jì)穩(wěn)定顯示為0后，對試樣依次施加不同的法向應(yīng)力，在每種應(yīng)力下進(jìn)行不同流速的滲流測試，并用計(jì)算機(jī)實(shí)時(shí)記錄差壓計(jì)的讀數(shù)。體積流速范圍為1.67×10-7至3.33×10-6m3/s。首次進(jìn)行滲流測試的法向應(yīng)力為2 MPa，并按每次2 MPa不斷增大。剪切位移變大后，試樣懸空的部分易發(fā)生斷裂破壞，為保證試樣的完整性，在2、4、6、8 mm 4種剪切位移時(shí)施加的最大法向應(yīng)力分別為10、8、8、6 MPa。在每個(gè)剪切位移測試結(jié)束后，將上下裂隙表面進(jìn)行沖水清理，去除附著于裂隙面上的碎屑。使用VR-3100高精度3維輪廓儀（點(diǎn)間距分辨率為0.10 μm）掃描得到壓縮試驗(yàn)后的3維裂隙表面，將測試前后的裂隙表面幾何模型相減即可得到表面的破壞區(qū)域。

2 數(shù)值模擬方法

2.1 裂隙面的受壓變形破壞數(shù)值模擬

通過高精度輪廓儀掃描0（初始狀態(tài)）、2、4、6、8 mm 5種剪切位移條件下的試樣表面，并建立3維數(shù)值模型。考慮到整個(gè)試驗(yàn)過程使用的是同一個(gè)試樣，每個(gè)剪切位移的法向加載測試后，裂隙表面都會發(fā)生一定的破壞并不斷累積，所以下一個(gè)剪切位移的表面數(shù)值模型使用前一個(gè)位移試驗(yàn)后掃描所得的裂隙面作為原型。

以往相關(guān)研究中，通常將巖石假設(shè)為純彈性體，裂隙的變形模擬一般通過剛性地移動裂隙2個(gè)表面計(jì)算接觸體上的彈性變形來實(shí)現(xiàn)。事實(shí)上，裂隙表面凸體在受壓過程中不僅會產(chǎn)生彈性變形，還可能發(fā)生破壞，所以在數(shù)值模擬計(jì)算過程中，應(yīng)充分考慮到裂隙表面接觸體的彈性變形和彈塑性變形破壞過程。裂隙中的流體在2個(gè)粗糙表面形成的復(fù)雜空腔結(jié)構(gòu)中流動，變形計(jì)算的精度直接影響著模型內(nèi)部空腔結(jié)構(gòu)的建模精度。由于加壓試驗(yàn)過程中試樣無法取出，所以為同時(shí)獲得裂隙面破壞情況和重現(xiàn)裂隙在夾持器內(nèi)受壓時(shí)的狀態(tài)，進(jìn)而實(shí)現(xiàn)不同受壓狀態(tài)下的裂隙滲流模擬，采用自行開發(fā)的基于接觸力學(xué)變分原理框架[20]的受壓變形計(jì)算程序解決這一問題。該方法以求解接觸系統(tǒng)的最小勢能為基礎(chǔ)，將塑性變形引起的能量耗散添加到Tian和Bhushan提出的總互補(bǔ)勢能公式中，得到以下的計(jì)算接觸應(yīng)力與裂隙變形之間的關(guān)系式[21]：

式中： Ω為接觸域； σn為法向應(yīng)力，MPa；(x,y)為接觸區(qū)域內(nèi)的復(fù)合表面法向位移，mm；(x,y)為接觸區(qū)域內(nèi)2個(gè)接觸微凸體基于幾何干涉的總位移，mm；（mm）為法向應(yīng)力達(dá)到巖石基質(zhì)硬度（即 σn>H，巖石基體硬度H取0.1E，E為楊氏模量，MPa）時(shí)塑性變形區(qū)域的復(fù)合增量表面位移。本方法充分考慮了接觸體的變形破壞行為，并通過與實(shí)際裂隙壓縮試驗(yàn)的對比驗(yàn)證了其有效性[21]。

計(jì)算時(shí)，根據(jù)剪切位移量錯動2個(gè)表面，設(shè)定巖石的力學(xué)參數(shù)E和υ，按照Δσn=0.10 MPa逐步加載，并記錄法向位移和法向應(yīng)力及每步加載中產(chǎn)生的塑形破壞量。根據(jù)裂隙表面的彈塑性計(jì)算結(jié)果可得到裂隙內(nèi)部空腔的高精度幾何數(shù)據(jù)。計(jì)算模型的點(diǎn)間距為0.20 mm，滿足高精度受壓變形計(jì)算的需要[21]。

2.2 裂隙滲流數(shù)值模擬

COMSOL Multiphysics是一款以有限元為基礎(chǔ)的多物理場耦合模擬仿真軟件，廣泛應(yīng)用于巖石裂隙滲流計(jì)算，其有效性得到了充分的檢驗(yàn)[22]。將受壓變形計(jì)算得到的上、下裂隙面處理后以dxf.格式導(dǎo)入COMSOL中，隨后在軟件內(nèi)建立裂隙空腔模型，如圖3所示。

圖3 裂隙空腔數(shù)值模型Fig. 3 Numerical model of void spaces in a fracture

模型長100 mm，寬50 mm，高度用顏色表示，在0～8 mm之間。通過求解COMSOL中的層流模塊自帶的NS方程獲取裂隙內(nèi)的3維流場分布。對于不可壓縮牛頓流體，穩(wěn)態(tài)流動由NS方程控制，可寫成[22]：

式中：u為流速矢量，m/s；t為時(shí)間，s ；P為壓力矢量，Pa；ρ為流體的密度，kg/m3；T為剪應(yīng)力矢量，Pa；f為體力矢量，m/s2。根據(jù)試驗(yàn)條件，設(shè)裂隙中的流體為等溫、穩(wěn)定的不可壓縮牛頓流體，密度為999.70 kg/m3，動力黏度為1×10-3N·s/m2。

本文采用COMSOL中的層流模塊，用常規(guī)自由四面體對模型進(jìn)行劃分，由于裂隙面接觸點(diǎn)處的幾何形狀復(fù)雜，因此，在這些位置進(jìn)行適當(dāng)?shù)木W(wǎng)絡(luò)細(xì)化，以提高計(jì)算的精度。以剪切位移2 mm的模型為例，模型包含了1 841 775個(gè)單元，在i7內(nèi)核的計(jì)算機(jī)上單次計(jì)算時(shí)長為9 h。入口邊界條件定義為考慮流體慣性和定流量的“充分發(fā)展的流動”，體積流速范圍設(shè)為10-8～4×10-5m3/s，出口邊界條件為“0 Pa”，即開放的邊界。該范圍相對于試驗(yàn)測試的范圍更寬，便于充分考察線性和非線性滲流特征。對于不可壓縮流體，針對非線性流動的Forchheimer方程和線性流動的滲透率計(jì)算公式[23]可以分別寫為：

式中：Q為入口流量，m3/s；A為線性項(xiàng)系數(shù)，Pa·s·m-4；B為非線性項(xiàng)的系數(shù)，也即慣性系數(shù)，Pa·s2·m-7；kv為黏性滲透率，m2；ki為慣性滲透率，m；υ為體積通量，m/s；μ為動力黏度，Pa·s；K為裂隙滲透率，m2；S為試樣橫截面面積，m2。

對于裂隙流，雷諾數(shù)Re可以定義為[24]：

式中，w為裂隙沿壓力梯度方向的寬度，m。

3 試驗(yàn)和數(shù)值模擬結(jié)果分析

3.1 裂隙變形和破壞區(qū)域分析

圖4為裂隙力學(xué)開度em與法向應(yīng)力σn、剪切位移量d的關(guān)系。

圖4 裂隙力學(xué)開度與法向應(yīng)力、剪切位移量關(guān)系Fig. 4 Relationships among fracture mechanical aperture, normal stress and shear displacement

如圖4（a）所示，隨著裂隙面上所施加的法向應(yīng)力的增加，力學(xué)開度非線性減小，具有式（6）所示的冪函數(shù)關(guān)系：

式中，a和b為系數(shù)。

系數(shù)a和b的值都隨著剪切位移的增加而增加，這主要是由裂隙面上凸起部位的隨機(jī)分布引起的。隨著剪切位移d從2增加到8 mm，系數(shù)a由1.35增加到2.02，增加了90.45%；系數(shù)b由-0.11增加到-0.06，增加了45.45%。裂隙面錯動后，接觸面逐漸集中，在法向應(yīng)力作用下發(fā)生破壞，導(dǎo)致上下表面逐漸接近，從而產(chǎn)生新的接觸，增大裂隙的法向剛度。隨著法向壓力的增大，開度的變化速率減緩，裂隙的變形過程趨于穩(wěn)定。

如圖4（b）所示，隨著剪切位移的增加，力學(xué)開度呈非線性關(guān)系增大，但增加速率逐漸降低，具有式（7）所示的對數(shù)函數(shù)關(guān)系：

式中，c和n為系數(shù)。

隨著法向應(yīng)力從2增加到10 MPa，系數(shù)c變化較小，穩(wěn)定于0.51至0.53之間；系數(shù)n從0.92降低到0.69，降低了25%，其變化趨勢和Bandis等[25]的試驗(yàn)結(jié)果一致。

為進(jìn)一步探究裂隙開度和裂隙面接觸位置的破壞形態(tài)，通過改變應(yīng)力大小，對2、4、6、8 mm 4種剪切位移發(fā)生后的裂隙進(jìn)行相應(yīng)的數(shù)值模擬，并繪制裂隙在不同應(yīng)力下的塑性破壞高程圖，見圖5（a）。隨著剪切位移增大，法向位移也隨之增大，上下表面的相關(guān)度降低，導(dǎo)致接觸點(diǎn)個(gè)數(shù)和塑性破壞點(diǎn)數(shù)量減少。當(dāng)保持恒定的剪切位移時(shí)，隨著法向應(yīng)力的增加，塑性破壞點(diǎn)位置不會發(fā)生顯著的變化。2 MPa法向應(yīng)力下所產(chǎn)生的塑性破壞區(qū)域多為尖狀的微凸體，其面積隨著法向應(yīng)力的增加而增加，這進(jìn)一步證實(shí)了圖4（a）中應(yīng)力-開度曲線逐漸趨于平緩的原因。

圖5 不同剪切位移下裂隙面塑性破壞演化規(guī)律Fig. 5 Evolutions of plastic failure zones under different shear displacements

圖5（b）為掃描試樣獲取的不同剪切位移下裂隙塑性破壞高程圖及相應(yīng)的數(shù)值模擬結(jié)果，試驗(yàn)和數(shù)值模擬計(jì)算得到的塑性破壞空間位置分布吻合較好，但試驗(yàn)所得的破壞面積和高度略大于數(shù)值模擬結(jié)果。這主要是因?yàn)榛◢弾r屬于脆性材料，由大量的礦物顆粒組成，部分顆粒的損壞會導(dǎo)致周邊的顆粒隨之脫落；數(shù)值計(jì)算過程僅將局部應(yīng)力大于臨界應(yīng)力的接觸面定義為破壞區(qū)域，忽略了對鄰近巖石基質(zhì)的影響。上述結(jié)果驗(yàn)證了裂隙彈塑性接觸模型能夠準(zhǔn)確捕捉粗糙裂隙損傷行為，重現(xiàn)了裂隙表面的塑性破壞區(qū)域及其演化規(guī)律。

3.2 滲流特性分析

圖6為不同剪切位移條件下水力梯度和滲透率的關(guān)系。

如圖6所示：數(shù)值模擬得到的滲透率K隨水力梯度J的變化規(guī)律和室內(nèi)滲流試驗(yàn)結(jié)果基本一致，均呈現(xiàn)逐漸減少的變化趨勢，且降低速率均先增大后減少。隨著剪切位移的增加，接觸點(diǎn)個(gè)數(shù)減少且開度增大，滲透率K也逐漸增大。試驗(yàn)得到的滲透率K略小于數(shù)值模擬得到的滲透率K，主要是因?yàn)樵趻呙璧倪^程中會損失一部分精度，導(dǎo)致實(shí)際裂隙表面比掃描得到的表面更粗糙，同時(shí)錯動過程也會帶來一定的誤差。

圖6 不同剪切位移條件下水力梯度和滲透率的關(guān)系Fig. 6 Relationships between hydraulic gradient and permeability under different shear displacement

雖然應(yīng)力和剪切位移會對滲透率和水力梯度曲線產(chǎn)生一定的影響，但總體變化趨勢相同。在流速較小時(shí)，流體處于達(dá)西流動狀態(tài)，裂隙巖體的滲透率K變化較??；當(dāng)水力梯度J增加到一個(gè)臨界值后，流體在慣性力作用下進(jìn)入了非線性滲流階段，滲透率K降低較快，且隨著流速的繼續(xù)增大，所有工況的K值都逐漸集中并趨于平緩。以d= 6 mm和σn= 10 MPa為例，裂隙在雷諾數(shù)Re為0.2、2.0、20.0和200.0時(shí)的流線分布如圖7所示，入口均設(shè)置了200條流線。

圖7 不同雷諾數(shù)下裂隙內(nèi)的流線分布Fig. 7 Streamline distributions of fluid flow under different Re values

由圖7可知，隨著Re的增加，出口流線數(shù)量減少，且流體所流經(jīng)的區(qū)域也減少，這與流體慣性力的增大有關(guān)。流體在高流速情況下會繞開大開度裂隙附近的部分小裂隙，而且更易被裂隙的接觸和突起等幾何結(jié)構(gòu)變化所擾動，導(dǎo)致數(shù)值模擬中流線的不連續(xù)和數(shù)量減少的現(xiàn)象。隨著剪切位移的增大，接觸點(diǎn)個(gè)數(shù)減少，但單個(gè)接觸點(diǎn)面積增大，流體集中在若干優(yōu)勢流動通道中，對裂隙滲流能力起主要控制作用。

為解決流體流動的非線性，定義非線性因子Nf[26]為：

該方程量化了非線性分量在總量中的比例，且被廣泛接受。當(dāng)Nf=0.1時(shí)，流體流動從線性狀態(tài)過渡到非線性狀態(tài)，臨界點(diǎn)相應(yīng)的水力梯度定義為臨界水力梯度Jc。

Forchheimer方程作為描述裂隙中流體非線性流動特性的方程，慣性系數(shù)B與非線性滲流特性密切相關(guān)，將法向應(yīng)力σn和剪切位移d對臨界水力梯度Jc和慣性系數(shù)B的影響均做了分析，如圖8所示。

圖8 臨界水力梯度和慣性系數(shù)隨剪切位移的變化規(guī)律Fig. 8 Evolutions of the critical hydraulic gradients and the coefficients under different shear displacement

由圖8可知：隨著剪切位移的增大，臨界水力梯度Jc和慣性系數(shù)B均呈冪函數(shù)關(guān)系逐漸減少，減少速率逐漸降低。隨著法向應(yīng)力的增大，Jc和B均逐漸增大。當(dāng)法向應(yīng)力在2到10 MPa范圍內(nèi)變化時(shí)，隨著剪切位移的增大，Jc和B的波動范圍均急劇減少。例如：當(dāng)剪切位移從2增加到8 mm時(shí)，Jc的波動范圍從2.07×10-2-1.46×10-2=6.10×10-3減少到4.90×10-3-3.70×10-3=1.20×10-3，降低了80.32%；B的波動范圍從4.63×1014-1.66×1014=2.97×1014Pa·s2·m-7減少到7.09×1013-4.66×1013=2.43×1013Pa·s2·m-7，降低了91.28%。

當(dāng)剪切位移保持恒定時(shí)，增加法向應(yīng)力會導(dǎo)致裂隙面的接觸區(qū)域增加，進(jìn)而使得裂隙內(nèi)部的幾何更加復(fù)雜，流體的路徑也隨之復(fù)雜化，增強(qiáng)了裂隙內(nèi)滲流的非線性能力。如圖9所示，隨著法向應(yīng)力的增加，Jc和B均呈近似線性的上升趨勢，且隨著剪切位移的增大，Jc和B的增加速率逐漸降低。

圖9 臨界水力梯度和慣性系數(shù)隨法向應(yīng)力的變化規(guī)律Fig. 9 Evolutions of the critical hydraulic gradients and the coefficients under different normal stresses

例如：當(dāng)剪切位移從2增加到4 mm時(shí)，臨界水力梯度Jc的增加速率從7.58×10-4降低到4.58×10-4MPa-1，減少了39.58%；慣性系數(shù)B的增加速率從3.71×1013降低到8.84×1012Pa·s2·m-7·MPa-1，減少了76.17%。當(dāng)剪切位移從4增加到8 mm時(shí)，臨界水力梯度Jc的增加速率從4.58×10-4降低到0.16×10-4MPa-1，減少了65.94%；慣性系數(shù)B的增加速率從8.84×1012降低到3.04×1012Pa·s2·m-7·MPa-1，減少了65.60%。結(jié)果表明，剪切位移越小，相同應(yīng)力作用下裂隙開度越小，Jc和B對法向應(yīng)力越敏感。

相對標(biāo)準(zhǔn)偏差（RSD）是刻畫裂隙內(nèi)開度分布特征的常用指標(biāo)，定義為局部開度的標(biāo)準(zhǔn)偏差除以平均開度。通過數(shù)據(jù)擬合，建立裂隙內(nèi)部幾何特征與非線性滲流控制參數(shù)（Jc和B）的關(guān)系，結(jié)果如圖10所示，其中，U為剪切位移與裂隙面總長度之比（即剪切應(yīng)變）。

圖10 臨界水力梯度Jc和慣性系數(shù)B和裂隙幾何參數(shù)RSD、U的數(shù)據(jù)擬合結(jié)果Fig. 10 Data fitting results for critical hydraulic gradient Jc, inertia coefficient B and fracture geometric parameters RSD and U

由圖10可知：Jc和B與RSD之間都存在類似的冪函數(shù)關(guān)系。這一結(jié)果表明，應(yīng)力和位移是影響裂隙滲流特性的外因，內(nèi)部空腔的幾何特征是控制非線性滲流的首要因素，其中，剪切位移可造成空腔形態(tài)的顯著變化，是一個(gè)不可忽略的因素。圖10中包含了G-1、G-2、G-3這3個(gè)試樣的數(shù)值計(jì)算結(jié)果，得到的規(guī)律具有普適性。在實(shí)際應(yīng)用中，可通過統(tǒng)計(jì)裂隙開度分布特征，基于圖10中的擬合公式計(jì)算裂隙的非線性滲流控制參數(shù)。

4 結(jié) 論

通過試驗(yàn)和數(shù)值模擬的相互驗(yàn)證，建立了一套評價(jià)法向應(yīng)力和剪切錯位條件下裂隙非線性滲流特性的有效方法。力學(xué)方面，通過對應(yīng)力、位移及裂隙表面形貌的高精度測試及數(shù)值計(jì)算，精確獲取了裂隙的宏觀變形行為及細(xì)觀接觸破壞規(guī)律，從而建立了裂隙內(nèi)部空腔結(jié)構(gòu)的3維高精度模型；滲流方面，通過對水壓、流速的精確測控加上精細(xì)化3維數(shù)值模擬，定量描述了裂隙力學(xué)、幾何參數(shù)與非線性滲流控制參數(shù)之間的關(guān)系，揭示了壓剪作用下裂隙內(nèi)部滲流場的演化規(guī)律。主要結(jié)論如下：

1）裂隙力學(xué)開度em和法向應(yīng)力σn呈冪函數(shù)關(guān)系em=aσnb，隨著剪切位移d從2增加到8 mm，系數(shù)a由1.35增加到2.02，系數(shù)b由-0.11增加到-0.06。裂隙力學(xué)開度em和剪切位移d呈對數(shù)函數(shù)關(guān)系em=clnd+n，隨著法向應(yīng)力σn從2增加到10 MPa，系數(shù)c變化較小，穩(wěn)定于0.51至0.53之間，系數(shù)n從0.92降低到0.69。

2）在線性流動區(qū)域內(nèi)，滲透率不受水力梯度J的影響，保持恒定的常數(shù)；但在非線性流動區(qū)域內(nèi)，隨著水力梯度J的增加，滲透率依次經(jīng)歷了緩慢降低、急劇降低和緩慢降低的過程。滲流試驗(yàn)結(jié)果和數(shù)值模擬結(jié)果吻合良好，證明了該數(shù)值模擬技術(shù)的可靠性。

3）隨著剪切位移的增大，區(qū)分線性流動和非線性流動的臨界水力梯度Jc和慣性系數(shù)B的波動范圍均急劇減小。例如，當(dāng)剪切位移從2增加到8 mm時(shí)，Jc降低了80.32%，B降低了91.28%。

4）隨著法向應(yīng)力增加，Jc和B均呈近似線性的上升趨勢，且隨著剪切位移增大，Jc和B的增加速率逐漸降低。如：當(dāng)剪切位移從2增加到4 mm時(shí)，Jc的增加速率減少了39.58%，B的增加速率減少了76.17%。當(dāng)剪切位移從4增加到8 mm時(shí)，Jc的增加速率減少了65.94%，B的增加速率減少了65.60%。結(jié)果表明，剪切位移越小，裂隙開度越小，Jc和B對法向應(yīng)力越敏感。

5）臨界水力梯度Jc和慣性項(xiàng)系數(shù)B與裂隙空腔幾何特征參數(shù)RSD之間存在冪函數(shù)關(guān)系，且受到剪切應(yīng)變U的影響，所建立的方程可為裂隙內(nèi)非線性滲流特性評價(jià)提供定量的參考。

目前的研究尚未考慮實(shí)際巖體中裂隙長度、空間分布、巖石種類等的影響。這些因素對裂隙滲流具有重要作用，需要在今后的工作中逐漸完善。