劉曉光

(盤錦河海土木工程咨詢有限公司，遼寧 盤錦 124010)

天然巖體結(jié)構(gòu)面的存在使得巖體的水力特性復(fù)雜多變，其中節(jié)理面的粗糙起伏對(duì)結(jié)構(gòu)面的力學(xué)、滲透特性具有重要影響[1]。近年來，隨著大型水利工程建設(shè)數(shù)量的不斷增加，結(jié)構(gòu)面強(qiáng)度和滲流問題成為工程設(shè)計(jì)和建設(shè)必須要面對(duì)和解決的課題。例如，法國在20世紀(jì)中期修建Malpasset拱壩過程中，由于沒有充分考慮結(jié)構(gòu)面滲透的影響，造成大壩設(shè)計(jì)標(biāo)準(zhǔn)過低，進(jìn)而誘發(fā)了重大潰決事故[2]。由此可見，結(jié)構(gòu)面力學(xué)和滲透特性研究對(duì)保證相關(guān)工程的建設(shè)和運(yùn)行安全具有十分重要的意義和價(jià)值。

水利工程中的巖體滲流問題是一個(gè)十分復(fù)雜的問題，其中結(jié)構(gòu)面形成的孔隙較天然巖石的其他孔隙小一個(gè)數(shù)量級(jí)，但是其造成的滲流卻要大幾個(gè)數(shù)量級(jí)[3- 5]。因此，結(jié)構(gòu)面滲透特性研究是天然巖體滲透特性研究的關(guān)鍵與核心。在傳統(tǒng)的巖土工程研究領(lǐng)域，均以線性達(dá)西定律對(duì)結(jié)構(gòu)面滲透性進(jìn)行評(píng)估。但是，近年來的諸多研究成果顯示，結(jié)構(gòu)面滲流過程中的流速與水力梯度之間并不是簡(jiǎn)單服從線性達(dá)西定律，而是呈現(xiàn)比較復(fù)雜的非線性關(guān)系[6]，基于傳統(tǒng)的達(dá)西定律進(jìn)行結(jié)構(gòu)面滲流分析容易導(dǎo)致工程安全隱患[7]。天然巖體中的節(jié)理等結(jié)構(gòu)面的剪切位移與錯(cuò)動(dòng)變形十分常見，特別是受到施工擾動(dòng)時(shí)，這種變形幾乎是不可避免的。這種變形不僅會(huì)對(duì)工程的穩(wěn)定性造成影響，同時(shí)與滲透性強(qiáng)烈相關(guān)[8]。在節(jié)理的剪切錯(cuò)動(dòng)變形中，其開度、曲折率以及壁面接觸均會(huì)發(fā)生變化，但是這些變化難以通過物理試驗(yàn)觀測(cè)和界定，因此剪切錯(cuò)動(dòng)變形如何影響節(jié)理滲流方面一直缺乏深入研究。在節(jié)理滲流的早期研究中，一般采取立方定律，但是受巖石中天然裂隙邊界曲折的影響，該定律的研究結(jié)果往往存在較大誤差。陳益峰等基于結(jié)構(gòu)面的可壓縮性，同時(shí)考慮節(jié)理裂隙面的粗糙度，提出了復(fù)雜應(yīng)力條件下的廣義立方定律，成為目前節(jié)理變形條件下滲流特征研究的主要理論[9]。

根據(jù)流體力學(xué)理論，節(jié)理流體運(yùn)動(dòng)過程可以通過N-S方程描述[10- 11]，但由于N-S方程是一組非線性偏微分方程，采用傳統(tǒng)的流體計(jì)算動(dòng)力學(xué)方法對(duì)其求解存在諸多問題。因此，本文通過COMSOL軟件求解N-S方程，對(duì)節(jié)理剪切錯(cuò)動(dòng)過程中的開度、曲折度以及接觸面變化對(duì)非線性滲流的影響進(jìn)行研究，以期進(jìn)一步探究結(jié)構(gòu)面剪切錯(cuò)動(dòng)節(jié)理非線性滲流機(jī)理。

1 計(jì)算模型與邊界條件

1.1 幾何模型

研究采用的物理模型為耦合的上下兩個(gè)節(jié)理面構(gòu)成，其示意圖如圖1所示。其中，上節(jié)理面的長(zhǎng)度為50mm，下節(jié)理面的長(zhǎng)度為100mm，寬度均為25mm。其中，上節(jié)理面受到剪切荷載的作用，沿著下節(jié)理面做剪切滑動(dòng)，節(jié)理面的粗糙系數(shù)在4～6之間。研究中選取節(jié)理面的上下面所對(duì)應(yīng)的部分構(gòu)建裂隙滲流幾何模型。采取自由三角形網(wǎng)格對(duì)節(jié)理面進(jìn)行網(wǎng)格剖分，最終獲得1417002個(gè)計(jì)算單元。其中，最小單元尺寸為0.0256mm，最大單元尺寸為0.320mm。

圖1 物理模型示意圖

1.2 流體控制方程

根據(jù)李睿劬等學(xué)者的研究成果，黏性不可壓流體在復(fù)雜裂隙中的穩(wěn)定運(yùn)動(dòng)是流體質(zhì)點(diǎn)運(yùn)動(dòng)的平衡過程。這種力的平衡關(guān)系可以通過數(shù)學(xué)形式反映在控制方程N(yùn)-S方程[12]，其數(shù)學(xué)表達(dá)式如下：

(1)

式中，ρ—流體的密度；μ—流體的動(dòng)力黏滯系數(shù)；U—流體質(zhì)點(diǎn)的流速矢量；P—外加壓力梯度；F—體力矢量。

不可壓縮流體連續(xù)方程的表達(dá)式為：

U=0

(2)

由于本次研究對(duì)象為水體滲流，因此按照25℃條件下水的密度和動(dòng)力黏滯系數(shù)取值，其中ρ為0.997g/cm3；μ為0.89×l0-3Pa·s。

2 節(jié)理錯(cuò)動(dòng)方案與開度變化特征

本次模擬研究主要考慮上節(jié)理面剪切位移過程中產(chǎn)生的剪脹行為，根據(jù)幾何模型的尺寸，設(shè)計(jì)了1、2、3、5mm四種不同的剪切步長(zhǎng)，同時(shí)保證各個(gè)剪切步長(zhǎng)條件下上下節(jié)理面的接觸面積不大于1%。

節(jié)理的張開度是節(jié)理滲流的重要影響因素[13]，研究中對(duì)不同剪切狀態(tài)下的上下節(jié)理面的張開度進(jìn)行統(tǒng)計(jì)分析，結(jié)果顯示節(jié)理的最大開度和最小開度分別為0.53、0mm，平均開度為0.213mm，標(biāo)準(zhǔn)差為0.056。進(jìn)一步對(duì)節(jié)理開度進(jìn)行頻率統(tǒng)計(jì)分析，獲得如圖2所示的不同開度頻率分布直方圖。由圖2可知，節(jié)理絕大部分位置的開度介于0.15～0.30mm之間，0.10mm以下的小開度和0.40mm以上的大開度所占的比例極小。通過分析開度頻率分布特征，發(fā)現(xiàn)Gaussian曲線可以對(duì)上述特征進(jìn)行較好的擬合，其表達(dá)式為：

(3)

式中，f(x)—節(jié)理面開度統(tǒng)計(jì)頻率；u—平均開度；σ—標(biāo)準(zhǔn)差。

圖2 節(jié)理開度分布直方圖

按照上述方法對(duì)四種不同剪切步長(zhǎng)下的開度進(jìn)行統(tǒng)計(jì)計(jì)算，獲得平均開度與標(biāo)準(zhǔn)差計(jì)算結(jié)果，見表1。從表1中的計(jì)算結(jié)果可知，節(jié)理的開度會(huì)隨著剪切位移的增加而增大，但是接觸面積的變化并不明顯，呈現(xiàn)出比較穩(wěn)定的特征。從開度的標(biāo)準(zhǔn)差計(jì)算結(jié)果來看，隨著節(jié)理面剪切位移的不斷增大，標(biāo)準(zhǔn)差由0.056逐漸增加到0.476，說明開度分布離散程度不斷增加。

表1 不同剪切步長(zhǎng)節(jié)理開度特征

3 模型計(jì)算結(jié)果與分析

3.1 主要滲流特征量分析

為了對(duì)節(jié)理滲流的非線性特征進(jìn)行模擬，研究中選擇雷諾數(shù)變化范圍為5～30。利用數(shù)值模擬的方法對(duì)不同剪切步長(zhǎng)和雷諾系數(shù)條件下的節(jié)理順流向流速進(jìn)行計(jì)算。結(jié)果顯示，當(dāng)剪切步長(zhǎng)為0mm時(shí)，節(jié)理面與接觸面上的順流方向流速均為0，但是在開度較大的局部流速分量較大，而局部開度變化比較明顯的部位，順流向流速分量較小，究其原因主要是這些部位的滲流通道較為曲折。此外，隨著模型計(jì)算中入口流量的增大，裂隙內(nèi)順流方向的流速分量由正變負(fù)，并不斷減小，說明滲流通道內(nèi)存在局部回流。究其原因，主要是節(jié)理裂隙內(nèi)存在局部漩渦，并且隨著雷諾數(shù)的增加，裂隙內(nèi)的漩渦也逐漸增大，因而造成逐漸增加的慣性力損失。由此可見，節(jié)理裂隙內(nèi)的滲流并不服從線性達(dá)西定律。

圖3 雷諾數(shù)為10時(shí)不同位移步長(zhǎng)節(jié)理流線分布圖

在不同剪切步長(zhǎng)條件下，對(duì)節(jié)理裂隙中的滲流路徑進(jìn)行模擬研究，獲得如圖3所示的雷諾數(shù)為10時(shí)不同位移步長(zhǎng)節(jié)理流線分布圖。由圖3可知，節(jié)理裂隙中的流速不均勻分布特征十分明顯，流線也呈現(xiàn)出強(qiáng)烈的曲折性特點(diǎn)，究其原因主要是開度分布不均。另一方面，隨著剪切步長(zhǎng)的增加，節(jié)理面的接觸率隨之增大，并直接導(dǎo)致了溝槽優(yōu)勢(shì)流的形成，溝槽流效應(yīng)對(duì)節(jié)理面滲流有増強(qiáng)作用，該模擬結(jié)果與相關(guān)學(xué)者的研究結(jié)果一致[14- 15]。

圖4 雷諾數(shù)為10時(shí)不同位移步長(zhǎng)節(jié)理滲流壓力分布圖

對(duì)不同位移步長(zhǎng)條件下的節(jié)理滲流壓力分布進(jìn)行模擬，如圖4所示。結(jié)果顯示，在位移步長(zhǎng)為0mm條件下，壓力分布的主要特征是模型進(jìn)口到出口逐漸減小的特征，存在比較明顯的壓力損失。此外，在壓力的空間分布上看，模型進(jìn)口部位損失較快而出口部位損失較慢，呈現(xiàn)出明顯的空間分布不均勻性。其他位移步長(zhǎng)條件下的壓力分布具有相似的特點(diǎn)，但是壓力分布的不均勻性隨著位移步長(zhǎng)的增加逐漸增大。究其原因，主要是節(jié)理面的接觸使?jié)B流路徑受到阻擋，進(jìn)而產(chǎn)生愈加明顯的繞流效應(yīng)，這也成為接觸面附近壓力降低的重要原因。

表2 擬合線性和非線性參數(shù)以及臨界雷諾數(shù)的計(jì)算結(jié)果

3.2 非線性特征分析

為了進(jìn)一步研究節(jié)理剪切錯(cuò)動(dòng)條件下的非線性滲流特點(diǎn)，研究中利用模型模擬的方法對(duì)5種剪切步長(zhǎng)下的滲流模型進(jìn)行25次數(shù)值模擬，根據(jù)模擬結(jié)果對(duì)模型的出口壓力梯度與流量之間的關(guān)系進(jìn)行統(tǒng)計(jì)。試驗(yàn)結(jié)果顯示Forchheimer方程能夠有效地描述這一滲流過程，擬合結(jié)果如圖5和圖6所示。由圖5和圖6可以看出結(jié)構(gòu)面剪切錯(cuò)動(dòng)節(jié)理滲流具有顯著的非線性特征。

圖5 位移步長(zhǎng)0mm條件下流量水力梯度關(guān)系曲線

圖6 其余位移步長(zhǎng)條件下流量水力梯度關(guān)系曲線

對(duì)每一個(gè)剪切步長(zhǎng)，采用Forchheimer方程擬合實(shí)驗(yàn)數(shù)據(jù)，得到Forchheimer方程線性參數(shù)a與非線性參數(shù)b，結(jié)果見表2。從表2中的結(jié)果可知，隨著位移步長(zhǎng)的增加，線性參數(shù)a與非線性參數(shù)b均呈現(xiàn)出減小的趨勢(shì)，并且減小幅度均在兩個(gè)數(shù)量級(jí)左右，同時(shí)兩個(gè)系數(shù)的減小呈現(xiàn)出先快后慢的特點(diǎn)。這說明隨著位移步長(zhǎng)的增加，節(jié)理張開度的變化對(duì)滲透性的影響逐漸減小。根據(jù)Zeng和Grigg的研究成果[13]，在每一個(gè)剪切步長(zhǎng)計(jì)算節(jié)理中線性與非線性滲流的臨界雷諾數(shù)Rec的計(jì)算方法見公式(4)，計(jì)算結(jié)果見表2。

(4)

由表2中的結(jié)果可知，線性參數(shù)a的變化范圍為5.10×108～5.37×1010kg·s-1·m-5；非線性參數(shù)b的變化范圍為1.80×1014～1.23×1016kg·m-8；臨界雷諾數(shù)的變化范圍為13.5～19.2。由此可見，在將來的工作中，建立一個(gè)能夠考慮節(jié)理形貌的臨界雷諾數(shù)模型對(duì)于界定節(jié)理非線性流有著重要的工程意義。

4 非線性滲流機(jī)理與討論

根據(jù)上文的模擬計(jì)算結(jié)果，在節(jié)理剪切錯(cuò)動(dòng)穩(wěn)定滲流狀態(tài)下，節(jié)理開度較小時(shí)體力可以忽略不計(jì)，流速不會(huì)隨著時(shí)間的變化而變化。因此，N-S方程可以進(jìn)一步簡(jiǎn)化，簡(jiǎn)化后的表達(dá)式如下：

P=μ2U+ρU·U

(5)

由此可以看出，水力梯度損失主要包括線性項(xiàng)和非線性項(xiàng)兩部分。其中，線性損失主要為黏性力損失；非線性損失主要為慣性力損失。根據(jù)流體力學(xué)理論，滲流過程中的慣性力損失大小主要取決于流體質(zhì)點(diǎn)的運(yùn)動(dòng)速度。因此，在流速較小的情況下，此時(shí)節(jié)理面中流體的流動(dòng)變?yōu)槿鋭?dòng)流，相應(yīng)地，慣性力損失影響變得十分微小，滲流的線性特征就比較顯著。另一方面，非線性慣性力的微分表達(dá)式可以進(jìn)一步表達(dá)為如下形式：

U·U=U·(ux，uy，uz)

(6)

顯然，慣性力的大小主要受流速U及其空間位置變化率U兩個(gè)因素的影響。研究中對(duì)雷諾數(shù)為30，不同剪切位移步長(zhǎng)下的節(jié)理內(nèi)流速進(jìn)行模擬計(jì)算，并對(duì)計(jì)算結(jié)果進(jìn)行平均。結(jié)果顯示，當(dāng)位移步長(zhǎng)為0、5mm時(shí)，平均流速分別為0.119m/s和0.022m/s。由此可見，當(dāng)位移步長(zhǎng)增加時(shí)平均流速會(huì)大幅減小，說明U不變時(shí)，模型節(jié)理面剪切位移量的增加有助于減小滲流流體的慣性力損失。

另一方面，針對(duì)雷諾數(shù)為30，不同剪切位移步長(zhǎng)下的節(jié)理內(nèi)流線分布進(jìn)行模擬計(jì)算分析，如圖7所示。結(jié)果顯示，隨著位移步長(zhǎng)的增加，節(jié)理內(nèi)流線變得更為曲折，說明慣性力在流速不變的情況下，位移量的增加可以導(dǎo)致慣性力空間位置變化率的增大，對(duì)增加慣性力的損失較為有利。

綜合上述，天然巖石中的節(jié)理開度均勻，那么節(jié)理滲流的非線性特征主要是由流速增大慣性力損失造成的；而節(jié)理內(nèi)部的流速較低時(shí)，節(jié)理滲流的非線性特征主要是由滲流路徑曲折造成的。

圖7 雷諾數(shù)為30時(shí)不同位移步長(zhǎng)的節(jié)理滲流流線分布圖

對(duì)上下節(jié)理面接觸部位的縱截面沿x方向的流速分量進(jìn)行計(jì)算。結(jié)果顯示：該方向上的流速有正有負(fù)，并分別代表正向滲流和反向回流。從具體計(jì)算結(jié)果來看，位移步長(zhǎng)為0、5mm情況下接觸部位的最小流速分別為-9.14×10-3、-3.34×10-4m/s。說明節(jié)理內(nèi)部的渦旋與流速之間存在比較密切的關(guān)系，流速越大，越容易引發(fā)回流渦旋，從而造成慣性力損失，進(jìn)而導(dǎo)致節(jié)理滲流的非線性特征出現(xiàn)。

5 結(jié)論

(1)節(jié)理裂隙中的流速不均勻分布特征十分明顯，流線也呈現(xiàn)出強(qiáng)烈的曲折性特點(diǎn)，導(dǎo)致節(jié)理滲流存在著比較明顯的溝槽流和優(yōu)勢(shì)流效應(yīng)。

(2)隨著雷諾數(shù)的增加，節(jié)理滲流逐漸表現(xiàn)出顯著的非線性特征，流量和水力梯度之間的關(guān)系可以用Forchheimer方程擬合。在剪切位移量從0mm增加到5mm的過程中，線性參數(shù)a與非線性參數(shù)b均減小了兩個(gè)數(shù)量級(jí)左右。

(3)節(jié)理滲流的非線性特征主要是由流速增大慣性力損失造成的；當(dāng)節(jié)理內(nèi)部的流速較低時(shí)，節(jié)理滲流的非線性特征主要是由滲流路徑曲折造成的。