陳俊國，劉衛(wèi)群，梁浩楠

(1.中國礦業(yè)大學力學與建筑工程學院,江蘇徐州 221008；2.山東科技大學礦業(yè)與安全工程學院,山東青島 266590)

地下巖體裂隙蠕變滲流耦合分析

陳俊國1，2，劉衛(wèi)群1，梁浩楠1

(1.中國礦業(yè)大學力學與建筑工程學院,江蘇徐州 221008；2.山東科技大學礦業(yè)與安全工程學院,山東青島 266590)

震后斷層區(qū)巖體裂隙的愈合對于地震水力響應研究具有重要的作用。為研究裂隙巖體愈合對深部斷層區(qū)滲透率時空演化規(guī)律的影響,在離散裂隙網絡藕合模型基礎上加入裂隙蠕變效應,建立裂隙巖體流固藕合時空演化模型,并利用COMSOL Multiphysics對建立的藕合方程進行求解。結果表明:封堵之前,常規(guī)的藕合滲流達到穩(wěn)定狀態(tài),由于具有完整的通道,任意時刻的流固藕合并不能改變流體的壓力；隨著封堵發(fā)生,在蠕變效應下,裂隙開度減小,單元體滲透率降低,流體壓力增大。該研究成果為震后破裂帶巖體的愈合機理及滲透率演化分析提供了理論依據。

裂隙巖體；蠕變效應；滲透率；裂隙開度；流固藕合

2015,32(11):45-51

1 研究背景

由裂隙主導的地下巖體損傷區(qū)水力特性對于地震的產生和傳播具有重要的影響[1-3]。地下流體不僅參與地震孕育發(fā)生的全過程,而且充當?shù)卣鹦畔鞑サ拿浇?。美國勞倫斯伯克利國家實驗室的Montgomery和Manga[4]詳細探討了溪流量和井水水位對地震的響應,認為由地震引起的地殼變形以及地表振蕩通過固結地表沉積物、巖石裂隙擴展、含水層變形以及裂隙中充填物的阻塞與清除來影響地下裂隙巖體的滲透率,并且通過觀察認為溪流量和井水水位的變化幅度與地震震級有密切的關系；Wang (2003)[5]根據1999年臺灣井水對集集地震的響應情況提出了井水水位的4種響應形式,認為影響地震水力響應的因素是繁多復雜的；Wang等(2008)[6]通過觀察發(fā)現(xiàn),震后井水水位隨時間表現(xiàn)出較大范圍的波動；Liu等[7-8]在實驗室中研究了裂隙砂巖在動態(tài)載荷下滲透率變化的規(guī)律,并且在裂隙中添加充填物研究振動對裂隙含沙滲流的的影響；Geballe等[9]利用汶川地震過后臺灣劉家井的水位變化數(shù)據,提出了新的模型解釋滲透率隨時間的指數(shù)型遞減規(guī)律,與Brodsky等[10]提出的模型進行對比分析,顯示了較好的結果；Elkhoury等[11]應用井水水位對地球固體潮的響應測量地層滲透率,通過對20 a期間地震過后滲透率的測量,認為地震波能夠增加地層的滲透率,并認為此發(fā)現(xiàn)為提高煤層氣及石油的開采效率提供了一種可能的方法；Xue等[3]在汶川震后地層破裂區(qū)應用深部鉆孔實時監(jiān)測地下水位潮汐響應的方法測量滲透率隨時間的變化,在大部分情況下,當斷層區(qū)巖體愈合時,滲透率將會迅速下降,并且揭示了大地震過后,在遠震影響下斷層區(qū)巖體愈合和破裂交替進行的過程。

以往研究[9-10]均是根據實地測量結果和經驗給定滲透率模型,此類模型主要是從宏觀角度對現(xiàn)象進行解釋,并不能從微觀機理及各個影響因素層面進行定量分析,而地層滲透率的變化是由固體巖石的力學環(huán)境和流體流動相互作用的復雜過程。本文基于前人的現(xiàn)場監(jiān)測和實驗室測試結論,綜合考慮應力環(huán)境和流體流動特性,將裂隙巖體愈合能力歸因于裂隙的蠕變效應,建立考慮蠕變效應的裂隙網絡巖體藕合滲流模型,探討震后破裂區(qū)裂隙巖體滲透率的時空演化規(guī)律。

2 巖體裂隙蠕變模型

裂隙巖體的蠕變特性一直是巖石力學研究的重要內容,目前對完整巖石和裂隙面蠕變模型的研究成果[12-14]較多。楊松林等[13]在Huang[15]給出的單一裂隙面的彈性本構關系基礎上將變形增量、柔度矩陣和應力增量表示為時間的函數(shù),給出了單一裂隙面的蠕變關系,即

或者

為了簡化裂隙蠕變模型,楊松林等[13]忽略裂隙面法向應力作用下的蠕變變形,將法向變形視為固定的彈性變形。實際上,當裂隙面較軟弱、法向力較大時,裂隙面的法向蠕變變形會比較明顯。因此,將法向剛度和切向剛度表示為時間的函數(shù),并代入式(2)得:

若裂隙面的法向和切向蠕變變形規(guī)律均符合廣義Kelvin模型,則有

式中:G1n,G2n和ηn分別為裂隙面的法向彈性模量、黏彈性模量和黏滯系數(shù)；G1s,G2s和ηs分別為裂隙面的切向彈性模量、黏彈性模量和黏滯系數(shù)。

本文建立的裂隙蠕變-滲流模型中,裂隙體積模量只與法向剛度有關,忽略剪脹效應,表征單元體內裂隙的變化只與裂隙和巖塊的體積模量有關。因此,本文采用裂隙面法向蠕變的廣義Kelvin模型作為裂隙蠕變的基本模型,探討裂隙巖體地層滲透率的時空演化規(guī)律。設裂隙的法向剛度為Kn,則有

3 裂隙蠕變-滲流耦合模型

本文基于雙重孔隙介質模型[16]建立裂隙蠕變-滲流藕合模型,模型中,REV表征單元體由巖塊和裂隙組成,斷裂區(qū)由裂隙主導,忽略巖塊的滲透率,水流在裂隙中進行。模型中綜合考慮了外在應力和流體壓力共同作用下巖塊的變形、裂隙的變形、裂隙的蠕變效應及裂隙-流體的藕合過程。模型推導基于以下假設:

(1)震后斷裂區(qū)裂隙密布,裂隙是流體流動的主要通道,忽略巖塊的滲透率。

(2)流體流動符合達西滲流法則。

(3)裂隙區(qū)充滿流體,即飽和滲流。

(4)蠕變只發(fā)生在裂隙面之間,忽略巖塊的蠕變。

3.1 裂隙網絡巖體變形的控制方程

根據假設,如圖1所示,巖體被概化為具有相同尺寸的巖塊和裂隙的結構[16]。巖塊為長度為a的立方體塊,巖塊之間為裂隙,寬度為b,σ為作用于巖體的外力。

圖1 斷層區(qū)的物理模型Fig.1 Physical model for fault zone

巖體變形的幾何方程為

式中εij為應變張量分量；ui為位移分量。

在忽略慣性力作用,巖體的平衡方程可表示為

式中:σij為應力張量的分量；fi為體力的分量。

根據多孔介質彈性理論[17],表征單元體的本構方程可表示為

裂隙的存在將會對巖體的力學性質產生很大的影響,本文根據雙重孔隙介質模型[1]中煤的彈性參數(shù)的定義對式(10)的參數(shù)進行如下定義:

式中:Em為巖塊的彈性模量；G為剪切模量；p為裂隙中流體壓力；K為巖體的體積模量；Km為巖塊的體積模量；Kn為單個裂隙的法向剛度；α和β分別為巖塊和裂隙的Biot系數(shù)；δij為Kronecker符號；ν為泊松比。

聯(lián)立式(8)至式(10),整理得到Navier型的公式為

式(18)即是裂隙巖體變形的控制方程。

3.2 流體流動的控制方程

地下水在裂隙巖體中流動的質量平衡方程為

式中:m為單位體積裂隙巖體中所含流體的質量；ρ為流體的密度；為Darcy定律的速度矢量；Qs為流體的補給源；t為時間。

由于流體只在裂隙中流動,單元體內流體的質量可以表示為

式中φf為裂隙巖體的裂隙度,表示單位體積裂隙巖體中,裂隙體積所占的比值。

由前文假設可知,地下水在裂隙巖體中的流動為穩(wěn)態(tài)達西滲流,根據達西定律,流體的速度矢量為

式中:k為裂隙巖體滲透率；μ為流體的動力黏滯系數(shù)。

下面探討裂隙巖體中單元體變形特性及滲透率的表達形式。由立方定律可知,裂隙的滲透率取決于裂隙開度,巖塊和裂隙的變形將會引起裂隙寬度的變化,為研究這些變化對滲透率的影響,本文從從離散裂隙網絡模型中取出一個表征單元體(REV)進行研究。如圖2所示,實線表示REV變形前的形態(tài),虛線表示變形后的形態(tài)。

如圖2所示REV的結構,任意方向上的變形由巖塊的變形和裂隙變形組成,則

圖2 表征單元體(REV)變形示意圖Fig.2 Schematic diagram of deformation of representative elementary volume

其中,

由于b?a,則d=a+b≈a,式(23)的應變公式可變形為

式中:d為巖塊和裂隙的寬度之和；a為巖塊的寬度；b為裂隙開度；ε為REV在任意一方向上的應變之和；εm,εf分別為巖塊和裂隙的體積應變,其表達式分別為

式中:Kf為裂隙的等效體積模量,Kf=a0Kn；Δσm為裂隙對巖塊的作用力；Δσf為巖塊對裂隙的作用力。Δσm=Δσf且與有效應力的變化Δσe是相等的。

假設在外力和流體壓力的作用下,REV的體積應變是εv,則

聯(lián)立以上各式得

式中:a0為巖塊的初始寬度；b0為裂隙的初始張開度。

將式(29)代入式(26)得

在REV變形過程中,裂隙開度為

聯(lián)立式(30)和式(31)得

同理:

對于如圖1所示的裂隙系統(tǒng),由立方定理表示的滲透率表達式[16]為

則將式(32)和式(33)代入式(34)得

表征單元體的裂隙孔隙度可以由裂隙的間距a (即巖石塊的寬度)和裂隙開度b計算得到,即

則

3.3 地下巖體裂隙蠕變-滲流耦合分析的控制方程

裂隙巖體變形的控制方程為

式中p為裂隙中流體壓力,由流體流動的達西定律的初始條件和邊界條件確定。

對于流體流動的控制方程,本文根據斷層區(qū)裂隙巖體的在不同階段具有不同的水力響應特性,將流體流動的控制方程根據時間劃分成2個階段,分別為:

(1)0≤t≤t0,地震剛過時,斷層區(qū)裂隙巖體滲透率較高且還沒有發(fā)生堵塞,可認為裂隙具有蠕變效應的穩(wěn)態(tài)滲流,即在初期,求解域(斷層區(qū)某一區(qū)域)具有完整的流通性,根據壓差可自由流出邊界,因此流體壓力不會隨時間發(fā)生變化,則有

至此,建立了斷層區(qū)在發(fā)生堵塞之前的蠕變-滲流藕合的控制方程。

(2)t0＜t≤t′0,裂隙經過一段時間的擠壓和愈合,區(qū)域的邊界被堵塞,不能發(fā)生流體的流動,此時流體具有壓縮性,設ρ=c0p,則由質量守恒方程:

得流體流動的控制方程為

結合斷層區(qū)表征單元體的變形控制方程(18)得裂隙愈合堵塞后水力響應的藕合模型。

至此,建立了地下巖體裂隙蠕變-滲流藕合的控制方程。下面探討模型的邊界條件和初始條件。

(1)對于變形控制方程,位移邊界條件為

(2)對于流體流動方程,采用Dirichlet和Neumann邊界條件,即:

求解域位移初始條件為

式中u0是求解域的初始位移。

求解域流體壓力初始條件為

式中p0是流體初始壓力。

4 模型求解

本文采用COMSOL Multiphysics軟件來進行數(shù)值模擬。COMSOL Multiphysics是一款大型高級數(shù)值仿真軟件,因其高效的計算性能和杰出的多場直接藕合分析能力實現(xiàn)了任意多物理場的高度精確的數(shù)值仿真,在工程中得到了廣泛的應用。本文的數(shù)值模型均是在該軟件上建立并運算。

4.1 數(shù)值模型與模擬參數(shù)

如圖3所示,模型尺寸為10 m×10 m,在模型的中點處設置一個圓形裂口,半徑為0.01 m,其目的就是通過裂口邊界的開合,決定流體的自由流通和流體堵塞在區(qū)域裂隙中。左、下、右邊界為簡支約束,上邊界受6 MPa均布力。模型的滲流邊界設置為四周均為無滲流邊界,根據時間控制內部井口邊界的開合,求解域內初始水壓為p0=2 MPa。模型求解所需要的參數(shù)見表1。

圖3 數(shù)值模型Fig.3 Numerical simulation model

表1 數(shù)值模擬參數(shù)Table 1 Parameters of numerical simulation model

4.2 數(shù)值模擬結果及分析

由于求解域根據受力情況及邊界條件是一個對稱結構,本文只對區(qū)域內點進行分析。設置時間0≤t≤50 d時,為有蠕變效應的穩(wěn)態(tài)流固藕合；50 d＜t≤1 000 d時,流體流動的邊界被封堵,此時無流體運動,裂隙中流體壓力隨蠕變進行發(fā)生變化。

由滲透率的推導過程可知,滲透率與體應變有直接的關系,在t＞0的時間里,蠕變一直進行著,因此區(qū)域內體應變是隨時間的一個變化量,區(qū)域內點A(5 m,4 m)的變化如圖4所示。

圖4 求解域中A點體應變隨時間的變化Fig.4 Volume strain vs.time of point A in solution domain

在t=0時刻,即裂隙蠕變效應尚未發(fā)生時,常規(guī)裂隙巖體流固藕合所產生的體應變,可認為是震后裂隙巖體的流固藕合產生的體應變。在t＞0時,蠕變發(fā)生后,由控制方程可知,體應變主要受蠕變效應的影響,蠕變時裂隙剛度是隨時間的變化量,因而體應變隨時間的變化關系如圖4所示,可以看出在t=1 000 d時,體應變絕對值較蠕變之前增加了1倍,影響較大。由式(32)和式(35)可知,裂隙開度和滲透率均是體應變的函數(shù),圖5為觀察點A(5 m,4 m)的裂隙開度和滲透率的變化曲線。

圖5 求解域中A點裂隙開度和滲透率隨時間的變化關系Fig.5 Fracture aperture vs.time and permeability vs. time of point A in solution domain

由圖5(a)可知,當t=0時,在法向應力和孔隙壓力的藕合作用下,裂隙的開度已經降低；當t＞0時,由于裂隙蠕變,裂隙開度隨剛度的降低和外力、孔隙壓力共同作用下,裂隙開度隨時間降低,即:在裂隙愈合過程中,由于應力擠壓的蠕變作用下,斷層區(qū)中裂隙開度會慢慢降低；繼而導致與裂隙開度密切相關的裂隙巖體滲透率的降低。由圖5(b)可知,滲透率隨時間降低,其曲線理應是圖5(a)的三次方關系,只是圖中并不能明顯地表達出來。由于裂隙開度與體應變之間的線性關系,對比圖4和5(a)可以發(fā)現(xiàn)二者隨時間的變化曲線是相似的。

在前文建立模型時,設定了區(qū)域封堵的時間,限制模型求解域中流體出口,使流動不能發(fā)生,持續(xù)的蠕變會造成裂隙開度降低,其中的流體被擠壓,流體壓力將會發(fā)生變化。如圖6所示,無論在何時進行封堵,流體壓力均會增大,且從與圖5(a)的對比中可以看出,流體壓力的增大與裂隙開度降低成反對應關系,這些結論與文獻資料中對于愈合過程中流體壓力增大是相吻合的。t＜t0時,流體可根據達西定律流動,且有完整的通道,任意時刻的流固藕合并不能改變流體的壓力,因此在這期間,求解域中流體壓力穩(wěn)定在2 MPa；t≥t0時,流體壓力增大并呈相同的增長趨勢,這與蠕變造成的裂隙開度降低形式不隨外界有效應力變化相吻合。另外,從圖6可以看出,封堵發(fā)生的時間越早,流體壓力達到某一閾值的時間就會越早。在實際中,由于地震發(fā)生后,一般會有大量的瑣碎巖石顆?；蛘吣嗌愁w粒流入斷層區(qū)造成封堵,這是愈合前期的封堵,會在地震過后數(shù)天內完成；另外一個封堵機制是在長期的愈合過程中,礦物質在各種物理、化學作用下結晶析出,附著在裂隙表面造成封堵,因此實際的封堵機制比較復雜難測,本節(jié)提供了一種封堵的效果模擬。

圖6 不同初始封堵時間對流體壓力的影響Fig.6 Influence of initial plugging time on fluid pressure

5 結 語

(1)推導了裂隙巖體流固藕合模型的一般表達式。假定本文表征單元體是由巖塊和圍繞著巖塊的裂隙組成,其變形為基質和裂隙的變形之和,并由此建立了基于立方定理的單元體滲透率表達式。

(2)將巖體裂隙愈合現(xiàn)象歸因于裂隙蠕變效應,并將廣義Kelvin模型作為裂隙蠕變的基本模型,取得較好的效果。但是廣義Kelvin模型并不能認為是深部裂隙巖體愈合的真實模型,因為在實驗室條件下很難模擬深部巖體所處的應力環(huán)境及地震作用的破裂形式,而廣義Kelvin模型是在實驗條件下得到的,與真實的深部裂隙巖體蠕變形式必然不同。并且由裂隙蠕變表達式可知,裂隙剛度并不能無限減小,這就使得裂隙不能達到完全愈合,與真實震后巖體能夠完全愈合不符,因此,探索真實的深部裂隙巖體愈合模型是未來的研究方向。

(3)本文從影響地層滲透率的機理分析,給出了綜合各種地質因素的數(shù)學模型,為地震水力響應提供了合理的解釋,具有重要的參考價值。

[1]RICE J R.Heating and Weakening of Faults During Earthquake Slip[J].Journal of Geophysical Research,2006, 111(B05):311.

[2]REMPEL A W,RICE J R.Thermal Pressurization and Onset of Melting in Fault Zones[J].Journal of Geophysical Research,2006,111(B09):535-540.

[3]XUE L,LI H B,BRODSKY E E,et al.Continuous Permeability Measurements Record Healing Inside the Wenchuan Earthquake Fault Zone[J].Science,2013,340: 1555-1559.

[4]MONTGOMERY D R,MANGA M.StreamFlow and Water Well Responses to Earthquakes[J].Science,2003,300: 2047-2049.

[5]WANG C Y,DREGER D S,WANG C H,et al.Field Relations among Coseismic Ground Motion,Water Level Change and Liquefaction for the 1999 Chi-Chi(Mw=7.5) Earthquake,Taiwan[J].Geophysical Research Letters, 2003,30(17):1890-1893.

[6]WANG C Y,CHIA Y.Mechanism of Water Level Changes During Earthquakes:Near Field versus Intermediate Field[J].Geophysical Research Letters,2008,35(12) :86-109.

[7]LIU W Q,MANGA M.Changes in Permeability Caused by Dynamic Stresses in Fractured Sandstone[J].Geophysical Research Letters,2009,36(20):146-158.

[8]朱 立,劉衛(wèi)群,王甘林.振動對充填裂隙滲透率影響的實驗研究[J].實驗力學,2012,27(2):201-206. (ZHU Li,LIU Wei-qun,WANG GAN-lin.Experimental Study the Influence of Vibration on the Permeability of Fractured Sandstone with Sediment Particles[J].Journal of Experimental Mechanics,2012,27(2):201-206.(in Chinese))

[9]GEBALLE Z M,WANG C Y,MANGA M.A Permeabilitychange ModelforWater-levelChangesTriggeredby Teleseismic Waves[J].Geofluids,2011,11(3):302-308.

[10]BRODSKY EE,ROELOFFS E,WOODCOCK D,et al.A Mechanism for Sustained Groundwater Pressure Changes Induced by Distant Earthquakes[J].Journal of Geophysical Research,2003,108(B8):2390-2392.

[11]ELKHOURY J E,BRODSKY E E,AGNEW D C.Seismic Waves Increase Permeability[J].Nature,2006,441: 1135-1138.

[12]徐衛(wèi)亞,楊圣奇.節(jié)理巖石剪切流變特性試驗與模型研究[J].巖石力學與工程學報,2005,24(增2):5536-5542.(XU Wei-ya,YANG Sheng-qi.Experiment and Modeling Investigation on Shear Rheological Property of Joint Rock[J].Chinese Journal of Rock Mechanics and Engineering,2005,24(Sup.2):5536-5542.(in Chinese))

[13]楊松林,張建民,黃啟平.節(jié)理巖體蠕變特性研究[J].巖土力學,2004,25(8):1225-1228.(YANG Song-lin, ZHANG Jian-min,HUANG Qi-pin.Analysis of Creep Model of Jointed Rock[J].Rock and Soil Mechanics, 2004,25(8):1225-1228.(in Chinese))

[14]熊良霄,楊林德.考慮節(jié)理面法向蠕變的節(jié)理巖體蠕變模型[J].中南大學學報(自然科學版),2009,40(3): 814-821.(XIONG Liang-xiao,YANG Lin-de.Creep Model for Rock Mass Considering Normal Creep of Rock Joint Plane[J].Journal of Central South University(Science and Technology),2009,40(3):814-821.(in Chinese))

[15]HUANG T H.Elastic Moduli for Fractured Rock Mass[J]. Rock Mechanics&Rock Engineering,1995,28(3):135-144.

[16]WU Y,LIU J S,ELSWORTH D,et al.Dual Poroelastic Response of Coal Seam to CO2Injection[J].International Journal of Greenhouse Gas Control,2010,4(4):668-678.

[17]DETOURNAY E,CHENG A H D.Fundamentals of Poroelasticity[M]//FAIRHURST C.Comprehensive Rock Engineering:Principles,Practice and Projects:Vol.2.Oxford:Pergamon Press,1993:113-171.

(編輯:黃 玲)

Creep-seepage Coupled Analysis of Underground Fractured Rock Mass

CHEN Jun-guo1,2,LIU Wei-qun1,LIANG Hao-nan1
(1.School of Mechanics&Civil Engineering,China University of Mining and Technology, Xuzhou 221008,China；2.College of Mining&Safety Engineering,Shandong University of Science&Technology,Qingdao 266590,China)

Crack-healing of rock mass in fault zones after earthquake plays an important role for hydraulic response for the earthquake.In order to study the healing effect of fractured rock mass on the spatio-temporal evolution law of permeability in deep fault zones,we introduced creep effect of fracture to the network coupled model of discrete fractures.On this basis,a new spatio-temporal evolution model for fractured rock mass based on fluid-solid coupling was built,and the coupled equations were solved with the software of COMSOL Multiphysics.The results show that, before sealing,common coupled seepage achieves a steady state.Due to complete seepage channel,the fluid-solid interaction can't change fluid pressure at any given time.As sealing happens,under the influence of creep,fracture aperture and permeability of unit body decrease,but fluid pressure increases.The conclusions can provide a theoretical basis for healing mechanism of rock mass and evolution analysis of permeability in fractured zones after earthquake.

fractured rock mass；creep effect；permeability；fracture aperture；fluid-solid coupling

TU45

A

1001-5485(2015)11-0045-07

10.11988/ckyyb.20140499

2014-06-20；

2014-08-10

國家自然科學基金項目(41074040)；國家重點基礎研究發(fā)展計劃(973)項目(2009CB219605)

陳俊國(1981-),男,山東濰坊人,博士研究生,主要從事地下裂隙巖體滲流力學方面的研究,(電話)13406482700(電子信箱)chenjunguo2002@163.com。