管單元用于縮減巖體滲透張量計(jì)算規(guī)模

王俊奇，劉博偉，岳 瀟

(1.華北電力大學(xué) 水利與水電工程學(xué)院，北京 102206；2.山東電力工程咨詢院有限公司，濟(jì)南 250014)

裂隙巖體滲流是導(dǎo)致水利、土木工程事故的主要原因之一[1-2]，還關(guān)系到采油、核廢料污染的防治、煤田瓦斯危害的防治和填埋垃圾的污水下滲等方面，工程中需要進(jìn)行細(xì)致的考察和分析。

常用的裂隙巖體滲流分析數(shù)學(xué)模型有離散裂隙網(wǎng)絡(luò)模型、等效連續(xù)介質(zhì)模型、雙重介質(zhì)模型以及混合滲流模型[3]。因滲流主要在裂隙網(wǎng)絡(luò)內(nèi)發(fā)生，用離散裂隙網(wǎng)絡(luò)模型來模擬巖體滲流規(guī)律時精度較高，并且由于離散型裂隙網(wǎng)絡(luò)模型忽略了巖塊本身的孔隙系統(tǒng)及裂隙巖體之間的水交替過程，分析時更加簡便。二維離散裂隙網(wǎng)絡(luò)模型由于構(gòu)造簡單發(fā)展較為成熟，但計(jì)算結(jié)果和實(shí)際情況相比誤差較大。三維離散型裂隙網(wǎng)絡(luò)模型方面，李新強(qiáng)等[4]使用邊界元法計(jì)算裂隙網(wǎng)絡(luò)的滲透張量，并與有限元法流量分析結(jié)果比較驗(yàn)證方法的有效性。于青春等[5-7]基于管單元模型通過對裂隙網(wǎng)絡(luò)的滲透張量進(jìn)行求解來分析巖體的滲流情況，Ren等[8]提出用統(tǒng)一管網(wǎng)法對裂隙巖體內(nèi)的滲流情況進(jìn)行研究，Sun等[9]用3D統(tǒng)一管網(wǎng)法模擬含裂隙的多孔巖石中的裂隙情況。Wang等[6-10]從滲透張量的角度對裂隙網(wǎng)絡(luò)的REV尺寸進(jìn)行研究。王振偉等[11]通過巖體的滲透張量來確定滲流的REV尺寸，并研究了裂隙的密度和跡長對裂隙巖體的滲透張量和REV大小的影響。也有研究發(fā)現(xiàn)巖體裂隙存在優(yōu)勢水力路徑，Zhao等[12-13]發(fā)現(xiàn)在裂隙巖體中只有10%～20%的裂隙可以發(fā)生滲流，更多的裂隙僅起到存水或持水的作用。Zhang等[14-15]利用數(shù)值分析方法以及現(xiàn)場試驗(yàn)測試進(jìn)行滲流分析，對特定巖體中優(yōu)勢水流的形成機(jī)理及其路徑進(jìn)行了相關(guān)研究，林宏奕[16]利用質(zhì)點(diǎn)追蹤技術(shù)來對巖體裂隙中存在的優(yōu)勢水流路徑進(jìn)行分析，劉華梅等[17]通過對裂隙網(wǎng)絡(luò)滲流路徑進(jìn)行優(yōu)化，發(fā)現(xiàn)剔除對滲流無貢獻(xiàn)的裂隙，可以減少計(jì)算量和縮短計(jì)算時長，倪紹虎等[18-19]通過對裂隙巖體的滲透特性進(jìn)行研究發(fā)現(xiàn)在裂隙巖體中的滲流的確存在優(yōu)勢水力路徑。另外，也有一些學(xué)者采用不同的方法對砂巖油藏中優(yōu)勢滲流通道進(jìn)行識別和分析，如雷霆等[20]提出改進(jìn)的CRM模型結(jié)合生產(chǎn)動態(tài)數(shù)據(jù)來對優(yōu)勢滲透通道進(jìn)行識別描述。現(xiàn)有研究工作表明，確定滲透張量方面，一維管單元模型可以顯著縮減計(jì)算工作量，而且由于巖體中存在較大裂隙，溝通巖體范圍較大，一定程度也起到優(yōu)勢水力路徑作用，因此，在確定滲透張量時，一定程度上可以忽略掉連通范圍小的裂隙，但縮減的幅度沒有相關(guān)研究。

本研究采用三維裂隙網(wǎng)絡(luò)模型對巖體裂隙滲流進(jìn)行分析，根據(jù)巖體露頭裂隙統(tǒng)計(jì)參數(shù)，用蒙特卡羅法生成巖體的圓盤裂隙網(wǎng)絡(luò)樣本，將滲流通道簡化為管單元模型。一定程度上，巖體主干裂隙控制滲透性能，通過對裂隙規(guī)模的研究，探索不同縮減規(guī)模對滲透張量和巖體REV大小的影響，以及計(jì)算規(guī)模可縮減的幅度，并通過實(shí)際工程案例對此方法進(jìn)行了驗(yàn)證。

1 三維裂隙網(wǎng)絡(luò)生成與滲透張量求解

1.1 模型的建立與參數(shù)生成

采用Baecher模型[21]，建立離散網(wǎng)絡(luò)模型時采用三維整體坐標(biāo)系，將生成的所有裂隙通過該坐標(biāo)系確定。將坐標(biāo)系的主軸方向與地理位置一一對應(yīng)，即X軸正半軸指向東，Y軸正半軸指向北，Z軸正半軸的方向指向上(與高程一致)，按照右手法則來進(jìn)行方向的角度旋轉(zhuǎn)。生成裂隙所需要確定的幾何參數(shù)包括裂隙的方向(傾向α,傾角β)、大小(半徑r)、位置(圓盤中心點(diǎn)坐標(biāo))、隙寬b，即(α,β,r,x,y,z,b)，每個參數(shù)變量的生成都遵循其自身的概率密度函數(shù)，利用積分法產(chǎn)生隨機(jī)變量。

裂隙的位置用裂隙面上有代表性的點(diǎn)來描述，在統(tǒng)計(jì)過程中用裂隙圓盤圓心來確定裂隙在空間中的位置，在三維空間中即圓心點(diǎn)(x,y,z)。假設(shè)圓心點(diǎn)在研究域各處出現(xiàn)的概率是相等的，即裂隙圓盤中心點(diǎn)的生成服從均勻分布。在生成裂隙網(wǎng)絡(luò)時，應(yīng)先確定生成域的大小，記錄圓心在該生成域內(nèi)的所有裂隙，然后在生成域內(nèi)再創(chuàng)建一定范圍內(nèi)的研究域，記錄圓心在該研究域內(nèi)的所有裂隙和與研究域邊界相交的裂隙。

生成裂隙三維裂隙網(wǎng)絡(luò)后，求出各圓盤的交線，傳統(tǒng)方法是在此基礎(chǔ)上對裂隙圓盤剖分成平面單元進(jìn)行裂隙網(wǎng)絡(luò)水力學(xué)分析，但是一個研究域的裂隙多達(dá)上萬條，用平面單元分析會產(chǎn)生十幾萬的自由度，使求解效率低下。為簡化分析計(jì)算過程，Tsang[22]和Cacas等[23]提出水流在裂隙面內(nèi)以溝槽流形態(tài)運(yùn)動，因此，可以對裂隙滲流的分析計(jì)算進(jìn)行簡化，如圖1所示。張有天[1]認(rèn)為可以用一維管單元進(jìn)行三維裂隙網(wǎng)絡(luò)水力學(xué)分析。具體簡化模型為將相交裂隙節(jié)理圓盤的圓心和交線中心相互連接起來，然后生成管單元，只在管單元內(nèi)產(chǎn)生水流，圓心與交線的中點(diǎn)即為節(jié)點(diǎn)。對于滲流區(qū)域內(nèi)所有節(jié)點(diǎn)建立有限元方程組，用矩陣表示為

圖1 三維裂隙網(wǎng)絡(luò)簡化分析概念簡圖Fig.1 Conceptual simplified analysis diagram of 3D fracture network

KH=Q

(1)

式中：K為總滲透矩陣，在三維空間中有9個分量；H為節(jié)點(diǎn)水頭矩陣；Q為節(jié)點(diǎn)流量矩陣。由于是恒定流，方程可解。

1.2 滲透張量的數(shù)值求法

數(shù)值求法的基本原理為：在模型的某一相對面施加等水頭邊界，兩個面之間存在一定的水頭差，另外兩對相對面施加變水頭邊界[24](如圖2所示)，滲流邊界的節(jié)點(diǎn)水頭分布情況如圖3所示。利用該模型共模擬了3組實(shí)驗(yàn)，第一組實(shí)驗(yàn)中令其中一對相對面施加等水頭，在另外兩對相對面施加變水頭，進(jìn)行計(jì)算并保存相應(yīng)的計(jì)算結(jié)果。另外兩組實(shí)驗(yàn)只改變等水頭邊界施加的方向，其他不變。由于對稱性，每組實(shí)驗(yàn)的滲流量只計(jì)算6個即可，3組實(shí)驗(yàn)共產(chǎn)生的流量數(shù)據(jù)為18個。根據(jù)滲透張量的定義得[21]

圖2 模型的邊界條件Fig.2 Boundary condition of the model

圖3 裂隙網(wǎng)絡(luò)滲流水頭分布示意Fig.3 Schematic of seepage head distribution in fracture network

(2)

式中：kij為滲透張量在j方向滲透坡降下產(chǎn)生的i方向分量，qij為滲流量在j方向滲透坡降下產(chǎn)生的i方向分量，ΔPi為在i方向的水力坡降。

對于均質(zhì)體，在各個方向上的入滲量和出滲量都相等。然而，對于節(jié)理巖體，由于巖體中裂隙的分布具有很強(qiáng)的隨機(jī)性，巖體的入滲流量和出滲流量并不總是相等。

2 計(jì)算規(guī)?？s減探索

于青春等[25]研究發(fā)現(xiàn)，滲流主要發(fā)生在較大的裂隙中，而較小的裂隙對滲流的影響十分微小，因此，在進(jìn)行裂隙滲流計(jì)算時可以在精度允許的情況下將較小的裂隙刪除掉，從而減少計(jì)算量，縮短計(jì)算時間。本文主要研究了在進(jìn)行滲流計(jì)算時可縮減的規(guī)模及其相應(yīng)的誤差，并用實(shí)際工程進(jìn)行了驗(yàn)證，當(dāng)然在實(shí)際工程設(shè)計(jì)中，應(yīng)根據(jù)工程的精度要求選擇合適的縮減量。在此約定：縮減規(guī)模為0.1即先將滲流裂隙網(wǎng)絡(luò)中的裂隙直徑按照由小到大的順序進(jìn)行排列，然后去掉排在前1/10直徑的裂隙，只用剩下的裂隙來進(jìn)行計(jì)算。

規(guī)模的縮減一般為全排列縮減，即對裂隙網(wǎng)絡(luò)空間內(nèi)所有組數(shù)的裂隙從小到大進(jìn)行整體的規(guī)?？s減。例如,若縮減規(guī)模為0.1，則將3組裂隙從小到大排列，刪除排在前10%的裂隙，保留剩下的90%進(jìn)行計(jì)算。全排列縮減方式的弊端在于容易較多地刪除直徑最小的那組裂隙，可能會降低滲透張量的準(zhǔn)確性?？紤]到力學(xué)成因，此組裂隙可能開度較大，滲透能力并不弱，因此，提出另一種縮減方法，即對各組裂隙按指標(biāo)分別從小到大規(guī)?？s減，例如,若縮減規(guī)模為0.1，分別對每組裂隙從小到大排列，再分別刪除每組前10%的裂隙，然后將各組的剩余裂隙重新統(tǒng)計(jì)在一起，進(jìn)行滲透張量計(jì)算。該方法在計(jì)算滲透張量前保證了各組裂隙縮減的平衡性，使各組裂隙數(shù)量均勻縮減。

三維裂隙網(wǎng)絡(luò)的參數(shù)信息采用李新強(qiáng)等[4]研究的三維裂隙網(wǎng)絡(luò)滲流模型，如表1所示，在該區(qū)域共產(chǎn)生的裂隙組數(shù)為3。

表1 三維裂隙網(wǎng)絡(luò)模型參數(shù)Tab.1 Model parameters of 3D fracture network

由邊界元法計(jì)算確定初值，再用現(xiàn)場壓水試驗(yàn)校核得到滲透主系數(shù)(m/s)和滲透主方向?yàn)閇4]

(3)

應(yīng)用自編的C++有限元滲流程序，控制管徑與隙寬的比值為11.506 9[26]，應(yīng)用上述數(shù)據(jù)進(jìn)行滲流模擬，選取100 m×100 m×100 m的立方體區(qū)域作為生成域，在生成域的中心選擇20 m×20 m×20 m的立方體作為研究域。對生成的裂隙進(jìn)行如上所述兩種方式的規(guī)?？s減，表2為每組裂隙分別按相同比例的規(guī)?？s減，表3為3組裂隙全排列情況下的規(guī)模縮減。表2和表3中滲透主值的誤差為

(4)

表2和表3中滲透主方向的誤差不可以直接求出，將計(jì)算所得滲透主值方向與現(xiàn)場壓水試驗(yàn)主值方向的夾角作為主方向的絕對誤差。兩滲透主值方向的夾角余弦為

表2 每組裂隙分別按相同比例的規(guī)?？s減及相應(yīng)誤差Tab.2 Scale reduction of each fracture group according to the same proportion and corresponding errors

表3 3組裂隙全排列規(guī)模縮減及相應(yīng)誤差Tab.3 Scale reduction of three fracture groups with full arrangement and corresponding errors

(5)

將未縮減時所得到的滲透主系數(shù)和滲透主方向與式(3)進(jìn)行比較，相應(yīng)的誤差如下：

(6)

式中：εi為主值相對誤差，θi為主方向絕對誤差，i=1,2,3。

由式(6)模擬可見,模擬計(jì)算的結(jié)果與現(xiàn)場壓水試驗(yàn)校正結(jié)果的誤差并不太大，在工程實(shí)踐可接受的范圍內(nèi)。

為了便于比較兩種排列方式對滲透張量的影響，將表2和表3中數(shù)據(jù)分別繪制兩種方式的滲透主值和滲透主方向的誤差圖像，結(jié)果如圖4和圖5所示。

圖4 兩種縮減方式滲透主值相對誤差對比Fig.4 Comparison of relative errors of principal permeability for two scale reduction methods

圖5 兩種縮減方式滲透主方向絕對誤差對比Fig.5 Comparison of absolute errors of permeability directions for two scale reduction methods

由圖4可知，兩種縮減方式滲透主值誤差在縮減規(guī)模0.1～0.7基本一致。當(dāng)縮減規(guī)模達(dá)到0.8時，分組排列的誤差絕對值為11%，全排列的誤差絕對值為9.7%。當(dāng)縮減規(guī)模達(dá)到0.8以后，誤差絕對值開始驟然增大，分組排列的誤差絕對值為32.6%，全排列的誤差絕對值為41.3%，此時縮減效果已成病態(tài)，不適宜再縮減裂隙規(guī)模進(jìn)行計(jì)算。同樣根據(jù)圖5可知，兩種縮減方式滲透主方向均值在縮減規(guī)模0.1～0.7波動較小。當(dāng)縮減規(guī)模達(dá)到0.7以后，分組排列的角度誤差絕對值為2.03°，全排列的誤差絕對值為2.60°。當(dāng)縮減規(guī)模達(dá)到0.8以后，誤差絕對值增幅增大，此時分組排列的誤差絕對值為9°，全排列的誤差絕對值為4.9°，此時方向誤差已經(jīng)不夠準(zhǔn)確，不適宜再縮減裂隙規(guī)模進(jìn)行計(jì)算。

綜上，對于兩種縮減方式效果基本一致，實(shí)際工程中可以采用任意一種縮減方式。當(dāng)最大縮減規(guī)模為0.7時，最為適宜本例，當(dāng)縮減規(guī)模大于0.7時，裂隙滲流網(wǎng)絡(luò)中重要的裂隙將被刪除，導(dǎo)致滲流分析的結(jié)果與裂隙巖體的實(shí)際滲流情況差別較大。

3 巖體表征單元體(REV)的確定

由于在天然巖體中存在的裂隙是隨機(jī)分布的，國內(nèi)外學(xué)者在進(jìn)行巖體力學(xué)參數(shù)選取時經(jīng)常會對巖體的表征單元體(REV)進(jìn)行分析。巖體的各項(xiàng)等效參數(shù)具有十分明顯的尺寸效應(yīng)，當(dāng)巖體的試樣體積增加到一定值V0時，各項(xiàng)滲透參數(shù)趨于穩(wěn)定，這一V0稱為巖體滲透性代表單元體積(REV)，確定巖體REV的尺寸是進(jìn)行裂隙滲流研究的一個重要方面。目前，國內(nèi)外學(xué)者對于REV尺寸有從等效力學(xué)參數(shù)進(jìn)行分析確定的，也有從水力學(xué)參數(shù)方面對REV進(jìn)行研究的，還有從裂隙塊體幾何參數(shù)對巖體的REV進(jìn)行研究的?；诓煌瑤r體參數(shù)的尺寸效應(yīng)而得到的表征單元體的尺寸都不太相同，周創(chuàng)兵等[27]提出了3種確定巖體REV的方法，即數(shù)值模擬法、地質(zhì)統(tǒng)計(jì)法和能量疊加法。

本研究運(yùn)用數(shù)值模擬的方法來計(jì)算巖體的滲透張量，根據(jù)滲透張量3個滲透主系數(shù)的變化情況來確定巖體的REV尺寸，通過對不同縮減規(guī)模情況下裂隙網(wǎng)絡(luò)REV的研究，進(jìn)一步探索規(guī)?？s減對裂隙巖體REV的影響。在縮減規(guī)模確定的情況下，對研究域從9 m×9 m×9 m開始逐漸增加到23 m×23 m×23 m，每次研究域的邊長增加1 m。以計(jì)算規(guī)模未進(jìn)行縮減時的情況為例，管徑與隙寬的比值為11.506 9，控制生成域的大小為100 m×100 m×100 m，通過不斷改變研究域的大小，得到不同研究域的滲透張量，將所得的裂隙網(wǎng)絡(luò)滲透主系數(shù)繪制成折線圖，如圖6所示?？梢钥闯?，當(dāng)所生成的裂隙網(wǎng)絡(luò)研究域的邊長增大到15 m后，隨著研究域邊長的增大，巖體滲透張量的滲透主系數(shù)變化情況基本趨于穩(wěn)定，故當(dāng)縮減規(guī)模為0時的REV尺寸可取15 m×15 m×15 m。同理，對其他縮減規(guī)模下的REV尺寸進(jìn)行分析，并將得到的REV的邊長繪制在折線圖上，如圖7所示?？梢钥闯觯S著裂隙網(wǎng)絡(luò)縮減規(guī)模數(shù)值的增大，巖體的REV尺寸逐漸變大，說明裂隙密度對巖體的表征單元體大小有一定的影響，REV尺寸隨裂隙密度減小而增大，該認(rèn)識與李錦輝等[28]對裂隙數(shù)量與REV尺寸關(guān)系[11]的研究結(jié)論相一致。

圖6 計(jì)算規(guī)模不縮減時的滲透主系數(shù)Fig.6 Principal permeability curves without calculation scale reduction

圖7 不同縮減規(guī)模下的REV邊長Fig.7 REV side lengths for different scale reductions

4 工程驗(yàn)證

中國水利水電科學(xué)研究院的科研人員在新疆對某水利樞紐的左岸邊坡巖體進(jìn)行現(xiàn)場調(diào)研后，根據(jù)巖體中結(jié)構(gòu)面的分布和組合情況得到如表4所示的統(tǒng)計(jì)結(jié)果[29]。

表4 巖體結(jié)構(gòu)面幾何參數(shù)統(tǒng)計(jì)結(jié)果Tab.4 Statistical results of geometric parameters of rock mass structural surface

對此工程實(shí)例，用自編的程序進(jìn)行數(shù)值模擬，將生成域的大小設(shè)為40 m×40 m×40 m，研究域的大小為15 m×15 m×15 m。不縮減時在生成的研究域內(nèi)裂隙共有8 129條，在裂隙網(wǎng)絡(luò)中節(jié)點(diǎn)數(shù)為53 650，單元數(shù)為99 659，經(jīng)過數(shù)值仿真模擬得到相應(yīng)的滲透主系數(shù)及滲透主方向?yàn)?/p>

(7)

本工程實(shí)例采用縮減規(guī)模為0.7，經(jīng)縮減后共剩2 394條裂隙，在裂隙網(wǎng)絡(luò)中共有15 986個節(jié)點(diǎn)、29 716個單元，經(jīng)數(shù)值仿真模擬得到相對應(yīng)的滲透主系數(shù)及滲透主方向?yàn)?/p>

(8)

為形象化滲透特性,根據(jù)滲透橢球方程

(9)

將式(7)、(8)用MATLAB軟件繪制縮減前后的滲透橢球圖，如圖8所示,由于橢球具有對稱性，為更好地顯示，取橢球的一半進(jìn)行觀察。

圖8 縮減前后滲透橢球?qū)Ρ菷ig.8 Comparison of permeability ellipsoid before and after reduction

由式(7)和(8)對比及圖8可以看出，縮減前與縮減規(guī)模為0.7的滲透張量及滲透橢球相差不大，說明采用縮減規(guī)模為0.7時，計(jì)算出的滲透主系數(shù)和滲透主方向與未縮減前計(jì)算的結(jié)果比較接近，這一方法對于工程實(shí)際具有一定的參考價(jià)值。

對縮減前后裂隙巖體模型的REV尺寸進(jìn)行計(jì)算得：

1)未縮減時裂隙巖體的REV邊長為4 m，相應(yīng)的滲透主系數(shù)和滲透主方向?yàn)?/p>

(10)

2)縮減規(guī)模為0.7時的REV邊長為6 m，相應(yīng)的滲透主系數(shù)及滲透主方向?yàn)?/p>

(11)

將式(11)與(10)進(jìn)行比較，滲透主系數(shù)及滲透主方向的誤差分別為

(12)

式中:εi為主值相對誤差，θi為主方向絕對誤差，i=1,2,3。

由式(12)可知，未縮減時邊長為4 m的REV與縮減為0.7時邊長為6 m的REV滲透主系數(shù)相對誤差的最大值為0.257，滲透主方向絕對誤差的最大值為4.62°，縮減前后REV的滲透主系數(shù)和滲透主方向誤差均較小。

5 結(jié) 論

1)通過刪減次要裂隙、保留主干裂隙進(jìn)行滲透張量計(jì)算，結(jié)果表明，裂隙網(wǎng)絡(luò)確實(shí)存在優(yōu)勢路徑，保留主干裂隙、縮減計(jì)算規(guī)模0.7時，計(jì)算得到的滲透張量仍然具有較好的精度，可以滿足工程需要。

2)對全排列和分組排列兩種方法的規(guī)?？s減程度與結(jié)果精度進(jìn)行分析，裂隙的滲透主值與主方向與現(xiàn)場壓水試驗(yàn)校核后的結(jié)果誤差在可以接受的范圍內(nèi)，而且兩種方法的縮減效果基本一致。

3)裂隙巖體REV尺寸的大小會隨著縮減規(guī)模的增大而變大，通過對巖體表征單元體的討論，發(fā)現(xiàn)縮減規(guī)模從0變化到0.7時，裂隙巖體的REV體積由15 m×15 m×15 m逐漸增大到21 m×21 m×21 m，說明REV的尺寸受縮減規(guī)模的影響。工程實(shí)例REV的邊長隨著縮減規(guī)模的增加從4 m增大到了6 m。

4)通過實(shí)例驗(yàn)證縮減計(jì)算的實(shí)用性。巖體裂隙網(wǎng)絡(luò)滲流滿足逾滲條件，存在REV情況下，對于那些對精度要求不高、地質(zhì)條件不復(fù)雜的情況，滲透張量可以采用縮減規(guī)模數(shù)為0.7進(jìn)行計(jì)算。對于工程精度要求較高、地質(zhì)條件復(fù)雜的重要工程項(xiàng)目，應(yīng)結(jié)合現(xiàn)場的實(shí)測資料與試驗(yàn)，選擇合適的縮減規(guī)模數(shù)進(jìn)行分析。在誤差允許的條件下相對地減少最大縮減量，可以達(dá)到減少計(jì)算量的目的，提高三維裂隙離散網(wǎng)絡(luò)模型在工程中的實(shí)用性。