劉紅巖,李俊峰,裴小龍
(1.中國地質大學(北京)工程技術學院,北京 100083;2.武警黃金部隊第四支隊,遼寧 遼陽 111000)
目前巖體工程如隧道、邊坡及地下廠房等在爆炸、地震等動載下的力學響應及破壞機理已引起了極大的關注。而巖體的強度、破壞模式及變形等力學特性主要受控于節理等天然缺陷[1-3]。由于節理發育規模一般遠小于斷層,延伸長度大多為數米,但是其數目眾多,因而無法一一考慮其對巖體力學性質的弱化效應。為此,20世紀80年代Kyoya等[4]最早提出了巖體損傷理論的概念,視節理為巖體的初始損傷,從而建立了節理巖體損傷模型。目前巖體損傷張量定義方法主要有2種。一種是Kawamoto等[5]和Swoboda等[6]提出的幾何損傷張量定義方法,他們采用節理密度、表面積、法向矢量等參數來定義損傷張量。但是該方法中僅考慮了節理的長度和方位等幾何參數對巖體力學特性的影響,而未考慮節理摩擦角及黏聚力等強度參數的影響,認為應力無法在節理位置傳遞,因而該方法更適用于巖體拉伸破壞的模擬,而無法用于巖體壓縮力學特性的分析。這是因為在受壓條件下,節理面將產生壓縮閉合或剪切滑移,從而傳遞部分壓、剪應力[7-8]。針對這一問題,Li等[9]和Liu等[10]基于損傷及斷裂力學中的能量原理,提出了一個新的能夠同時考慮節理的幾何參數(如節理長度、傾角、組數等)及強度參數(如節理摩擦角)的巖體損傷張量計算方法。相對于前述的幾何損傷理論而言,第2種方法減小了由于節理傳壓及傳剪系數的人為選取所帶來的誤差,因而具有較強的合理性[10-11]。
然而,盡管相對于幾何損傷理論而言,第2種方法已經有了較大的發展,但是仍不夠完善。這主要是由于目前學術界為了更準確地刻畫節理的物理力學性質,提出了描述節理物理力學性質的3類參數,即幾何參數、強度參數和變形參數(如法向和切向剛度)。而Prudencio等[12]通過對節理巖體的壓縮實驗,發現巖體的力學特性同樣受節理變形參數的影響。而節理變形參數在目前的節理巖體損傷張量計算中均未得到很好地考慮,這顯然會導致一定的計算誤差。因此,本文中擬基于Li等[9]和Liu等[10]的研究方法,在現有的、能夠考慮節理幾何及強度的損傷張量計算方法的基礎上,考慮節理變形參數的影響,最終在經典巖石細觀動態損傷本構模型—Taylor-Chen-Kuszmaul (TCK)模型的基礎上,提出能夠考慮宏觀缺陷(如節理)與細觀缺陷(如微裂紋)耦合的斷續節理巖體單軸壓縮動態損傷本構模型,并采用相關算例研究節理幾何參數(如節理長度)、強度參數(如節理摩擦角)和變形參數(如節理法向及切向剛度)對巖體單軸壓縮動態應力應變曲線的影響。
1986年,Taylor等[13]提出了著名的巖石動態損傷本構模型,即TCK模型,其相應的單軸動態壓縮條件下軸向應力σ與應變ε的關系[14]為:
(1)
Li等[9]推導了當巖體含有一條如圖1所示的傾角為α的斷續節理時,平面應力條件下,由其引起的加載方向上的巖體損傷變量:
(2)
式中:KⅠ、KⅡ分別為節理尖端的Ⅰ、Ⅱ型應力強度因子;V為試件體積;A為節理表面積;單個節理時,A=Ba(單邊節理)或2Ba(中心節理);多個節理時,A=NBa(單邊節理)或2NBa(中心節理);N為節理個數,B為節理深度,a為節理半長。
節理在壓縮荷載下往往會在已有節理尖端產生翼裂紋并近似沿外力方向擴展[15-17],如圖1所示。
對于完整巖石,在單軸壓縮下,傾角為α的斜面上的正應力和切應力分別為:

(3)
(4)
而對于節理巖體,即當在傾角為α的斜面上存在一條長度為2a的斷續節理時,該面上的正應力σα和切應力τα將明顯受到節理力學性質的影響,可表示為:
σα=(1-Cn)σcos2α
(5)
τα=(1-Ct)σsinαcosα
(6)
由于節理面黏聚力對試件強度影響較小,因此可忽略不計,同時若設節理摩擦角為φ,節理面摩擦系數為tanφ,那么在單軸壓縮下,由式(5)、(6)可得節理處的下滑力為:
(7)
根據Lee等[19]的研究,則可得:
(8)
式中:引入l*=0.27a,以保證l=0時,KI和KII不存在奇異;其余參數[15-17]如圖1所示。
對于節理未擴展的初始狀態,即翼裂紋長度l=0時,則有:
(9)
把式(9)代入式(2)即可得到由初始節理所引起的巖體損傷變量,進而通過張量化處理即可得到相應的損傷張量。可知由該方法求出的損傷張量同時考慮了節理的幾何參數(如節理長度、傾角等)、強度參數(如節理摩擦角)和變形參數(如節理法向及切向剛度),因此將更加符合實際情況。
當巖體含有如圖2所示的單排及多排斷續平行節理時(其中,b、d分別為相鄰兩節理中心距及排間距),節理尖端應力強度因子的計算方法已在相關文獻中有較多闡述[20-21],這里不再重述。
當巖體含單組節理時,把式(5)~(6)、(9)及相應的節理尖端應力強度因子計算公式代入式(2),進而通過求積分可得相應的損傷變量。然而由于被積函數比較復雜,式(2)沒有相應的解析解,因此這里采用Matlab軟件求其近似解。
由2.4節求得的損傷變量僅是加載方向上的損傷值,為反映節理對巖體損傷的奇異性,必須對其進行張量化處理。這里采用Chen等[22]提出的方法,即引入損傷張量:
(10)
式中:D0即為由式(2)求得的損傷值。
這里引用Liu等[10]提出的宏細觀缺陷耦合損傷張量計算公式:
(11)
基于連續介質力學損傷理論,把式(1)所示的完整巖石損傷變量ω(t),用式(11)所示的宏細觀缺陷耦合損傷變量Ω12替換,即可得到相應的單軸壓縮斷續節理巖體動態損傷本構方程:
(12)
式中:E為完整巖石的彈性矩陣,其余參數同前。

根據上述參數,計算結果如圖4所示。可以看出:(1)采用本文中提出的模型所計算出的斷續節理巖體動態應力應變曲線也包含了應力線性上升段、極值點及應力下降段等,很好地反映了巖體變形全過程;(2)斷續節理巖體的峰值強度及彈性模量分別為25.65 MPa和8.32 GPa,而相應完整巖石的動態峰值強度及彈性模量則分別為32.7MPa和10.8GPa,前者分別約為后者的78.44%和77.04%,均低于完整巖石的峰值強度,即斷續節理會顯著削弱巖體的動態強度及剛度;(3)另外,采用同樣的算例,劉紅巖等[11]計算出的相應巖體峰值強度及彈性模量分別為10.28 MPa和3.25 GPa,分別約為相應完整巖石的動態峰值強度及彈性模的31.44%和30.09%,如圖4所示,即遠小于本文計算結果。劉紅巖等[11]在計算由斷續節理對巖體造成的損傷時考慮了節理傾角的影響,認為式(12)中的E應為斷續節理巖體的彈性矩陣,進而視巖體為巖塊與節理組成的復合體系來計算E,在此過程中,同時又引入了節理傾角對彈性模量的影響,因此可以認為劉紅巖等[11]考慮了2次節理傾角對斷續巖體動態應力應變曲線的影響,即重復計算了一次節理傾角對斷續巖體動態應力應變曲線的影響,因此導致其計算出的巖體動態應力偏小,因而計算結果不甚合理。而相比而言,在本文所提出的模型中,對節理的3類參數,即幾何參數、強度參數及變形參數等均只考慮了一次,不存在重復計算的問題。
計算模型如圖3所示,設節理摩擦角分別為0°、15°和30°,其余計算參數同4.1節,計算結果如圖5所示。由圖5可知:(1)由巖體動態應力應變曲線斜率所反映的動態彈性模量是隨著節理摩擦角的增大而增大的,即3種不同節理摩擦角條件下的巖體動態彈性模量分別為7.24、8.32和9.61 GPa;(2)當節理摩擦角逐漸增大時,相應試件的動態峰值強度也由22.34 MPa分別增大到25.65和29.42 MPa。這是因為試件宏觀損傷程度隨節理摩擦角的增大而降低,進而導致試件強度提高,即節理摩擦角對試件的動態力學特性有較大影響。
計算模型仍為圖3所示,首先取節理法向剛度Kn分別為2、20和200 GPa/cm,其余計算參數同4.1節;其次取節理切向剛度Ks分別為0.8、8和80 GPa/cm,其余計算參數仍同4.1節,計算結果如圖6所示。可以看出:(1)當節理法向剛度由2 GPa/cm依次增大到20和200 GPa/cm時,即節理法向剛度依次增大到原來的10倍及100倍時,試件的峰值強度則由25.65 MPa依次降低到23.35和22.34 MPa,這說明隨著節理法向剛度的增大,試件動態峰值強度是降低的。這是因為在同樣的垂直應力下,隨著節理法向剛度的增大,節理面上的法向應力會相應減小,進而節理面上抵抗節理發生滑移的摩擦力就會減小,最終導致節理更易發生滑移破壞,進而峰值強度降低。(2)當節理切向剛度由0.8 GPa/cm依次增大到8和80 GPa/cm時,即節理切向剛度依次增大到原來的10倍及100倍時,試件的峰值強度則由25.65 MPa依次增大到31.82和32.32 MPa,這說明節理切向剛度對試件的動態峰值強度有較大影響,且隨著節理切向剛度逐漸增大,節理切向剛度對巖體動態峰值強度的影響則逐漸趨于定值。(3)如前所述,由圖6(a)、(b)分別看出,隨著節理法向及切向剛度的增大,試件的動態峰值強度分別減小及增大。然而Prudencio等[12]研究發現隨著節理法向及切向剛度的增大,試件動態峰值強度則是增大的。分析認為對圖6(a)、(b)兩種工況,節理的法向及切向剛度的選取是任意的。而節理法向及切向剛度則首先是由Goodman等[23-24]提出用以描述法(切)向應力與法(切)向位移之間的關系。Bandis等[25]發現節理的切向剛度遠小于其法向剛度,并研究了Kn/Ks的值隨法向應力的變化關系。這說明二者應該是按比例變化的。然而圖6(a)、(b)中所給出的結果卻并沒有考慮這一因素。為此,下面就研究二者按比例變化時巖體的動態應力應變曲線。取節理法向及切向剛度的初始值分別為2和0.8 GPa/cm,然后分別增大到10倍及100倍,結果如圖6(c)所示。可以看出,巖體的動態峰值強度則由26.65 MPa分別增大到30.36和32.64 MPa。這說明隨著節理法向及切向剛度的增大,巖體動態峰值強度是增大的,這與Prudencio等[12]的研究成果則是一致的。(4)從應力應變曲線的斜率還可以看出,節理法向及切向剛度對試件的動態彈性模量也均有一定的影響。總之,節理法向及切向剛度對試件動態峰值強度及彈性模量均有一定影響。
計算模型如圖3所示,取節理長度2a分別為10、20和30 mm,其余計算參數同4.1節,計算結果如圖7所示。由圖7可知,當2a由10 mm依次增大到20和30 mm時,試件峰值強度則由30.57 MPa依次降低到25.65和20.63 MPa,即分別減小為完整巖石峰值強度的93.49%、78.44%和63.09%,即隨節理長度的增大,試件動態峰值強度明顯下降。這是由于節理長度的增大引起巖體損傷變量的增大及節理傳壓和傳剪系數的變化,進而導致強度降低。同時巖體動態彈性模量也隨節理長度的增大而減小。
計算模型如圖3所示,取載荷應變率分別為100、200和300 s-1,其余計算參數同4.1節,計算結果如圖8所示。由圖8可知,當載荷應變率由100 s-1依次升高到200和300 s-1時,試件峰值強度則由25.65 MPa依次增大到31.65和35.71 MPa,即試件具有明顯的應變率效應。關于應變率對巖體動態力學特性的影響機理已有較多闡述[11-26]。
(1)針對目前節理巖體宏觀損傷張量計算中僅考慮了節理幾何及強度參數的不足,根據前人的研究成果,提出了能夠同時考慮節理幾何參數(如節理長度及傾角等)、強度參數(如節理摩擦角)及變形參數(如節理法向及切向剛度)的宏觀損傷張量計算方法,這將使得宏觀巖體損傷張量的計算更符合實際。
(2)基于完整巖石動態TCK模型及宏細觀損傷耦合的觀點,建立了考慮宏細觀缺陷耦合的非貫通節理巖體動態單軸壓縮損傷本構模型。該模型中同時考慮了節理幾何、強度及變形參數對非貫通節理巖體動態力學特性的影響。
(3)算例分析表明,節理的存在降低了巖體的動態彈性模量,削弱了巖體的峰值強度。同時利用該模型討論了節理摩擦角、節理法向及切向剛度、節理長度及應變率等對巖體動態力學特性的影響,計算結果與已有研究成果吻合較好,說明該模型是合理性的。
[1] AKIN M. Slope stability problems and back analysis in heavily jointed rock mass: A case study from Manisa, Turkey[J]. Rock Mechanics and Rock Engineering, 2013,46(2):359-371.
[2] 張力民,呂淑然,劉紅巖.綜合考慮宏細觀缺陷的巖體動態損傷本構模型[J].爆炸與沖擊,2015,35(3):428-436.
ZHANG Limin, Lü Shuran, LIU Hongyan. A dynamic damage constitutive model for rock mass by comprehensively considering macroscopic and mesoscopic flaws[J]. Explosion and Shock Waves, 2015,35(3):428-436.
[3] 劉紅巖,楊艷,李俊峰,等.基于TCK模型的斷續節理巖體動態損傷本構模型[J].爆炸與沖擊,2016,36(3):319-325.
LIU Hongyan, YANG Yan, LI Junfeng, et al. Dynamic damage constitutive model for rock mass with non-persistent joints based on the TCK model[J]. Explosion and Shock Waves, 2016,36(3):319-325.
[4] KYOYA T, ICHIKAWA Y, KAWAMOTO T. A damage mechanics theory for discontinuous rock mass[C]∥Proceedings of the 5th International Conference on Numerical Methods in Geomechanics. Nagoya, 1985:469-480.
[5] KAWAMOTO T, ICHIKAWA Y, KYOYA T. Deformation and fracturing behavior of discontinuous rock mass and damage mechanics theory[J]. International Journal for Numerical Analysis Method in Geomechanics, 1988,12(1):1-30.
[6] SWOBODA G, SHEN X P, ROSAS L. Damage model for jointed rock mass and its application to tunneling[J]. Computers and Geotechnics, 1998,22(3/4):183-203.
[7] YUAN X P, LIU H Y, WANG Z Q. An interacting joint-mechanics based model for elastoplastic damage model of rock-like materials under compression[J]. International Journal of Rock Mechanics & Mining Sciences, 2013,58(9):92-102.
[8] SWOBODA G, YANG Q. An energy-based damage model of geomaterials:Ⅰ: Formulation and numerical results[J]. International Journal of Solids and Structures, 1999,36(9):1719-1734.
[9] LI N, CHEN W, ZHANG P, et al. The mechanical properties and a fatigue-damage model for jointed rock mass subjected to dynamic cyclical loading[J]. International Journal of Rock Mechanics & Mining Sciences, 2001,38(7):1071-1079.
[10] LIU Hongyan, ZHANG Limin. A damage constitutive model for rock mass with non-persistently closed joints under uniaxial compression[J]. Arabian Journal for Science and Engineering, 2015,40(1):3107-3117.
[11] 劉紅巖,王新生,張力民,等.斷續節理巖體單軸壓縮動態損傷本構模型[J].巖土工程學報,2016,38(3):426-436.
LIU Hongyan, WANG Xinsheng, ZHANG Limin, et al. A dynamic damage constitutive model for rock mass with non-persistent joints under uniaxial compression[J]. Chinese Journal of Geotechnical Engineering, 2016,38(3):426-436.
[12] PRUDENCIO M, JAN M V S. Strength and failure modes of rock mass models with non-persistent joints[J]. International Journal of Rock mechanics & Mining Sciences, 2007,46(6):890-902.
[13] TAYLOR L M, CHEN E P, KUSZMAUL J S. Microjoint induced damage accumulation in brittle rock under dynamic loading[J]. Computer Method in Applied Mechanics & Engineering, 1986,55:301-320.
[14] GRADY D E, KIPP M E. Continuum modeling of explosive fracture in oil shale[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1987,17(3):147-157.
[15] HUANG C, SUBHASH G, VITTON S J. A dynamic damage growth model for uniaxial compressive response of rock aggregates[J]. Mechanics of Materials, 2002,34(5):267-277.
[16] HUANG C, SUBHASH G. Influence of lateral confinement on dynamic damage evolution during uniaxial compressive response of brittle solids[J]. Journal of the Mechanics and Physics of Solids, 2003,51(6):1089-1105.
[17] PALIWAL B, RAMESH K T. An interacting micro-joint damage model for failure of brittle materials under compression[J]. Journal of the Mechanics and Physics of Solids, 2008,56(3):896-923.
[18] LIU Taoying, CAO Ping, LIN Hang. Damage and fracture evolution of hydraulic fracturing in compression-shear rock cracks[J]. Theoretical and Applied Fracture Mechanics, 2014,74:55-63.
[19] LEE S, RAVICHANDRAN G. Joint initiation in brittle solids under multiaxial compression[J]. Engineering Fracture Mechanics, 2003,70(13):1645-1658.
[20] 李建林,哈秋瓴.節理巖體拉剪斷裂與強度研究[J].巖石力學與工程學報,1998,17(3):259-266.
LI Jianlin, HA Qiuling. A study of tensile-shear joint and strength related to jointed rock mass[J]. Chinese Journal of Rock Mechanics and Engineering, 1998,17(3):259-266.
[21] 范景偉,何江達.含定向閉合斷續節理巖體的強度特性[J].巖石力學與工程學報,1992,11(2):190-199.
FAN Jingwei, HE Jiangda. The strength behavior of rockmasses containing oriented and closed intermittent joints[J]. Chinese Journal of Rock Mechanics and Engineering, 1992,11(2):190-199.
[22] CHEN W, LA BORDERIE C, MAUREL O, et al. Simulation of damage-permeability coupling for mortar under dynamic loads[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2014,38(5):457-474.
[23] GOODMAN R E. The mechanical properties of joints[C]∥Proceeding of the 3rd Congress ISRM. Denver, 1974,I(A):127-140.
[24] GOODMAN R E, TAYLOR R L, BREKKE T. A model for the mechanics of jointed rock[J]. Journal of Soil Mechanics and Foundations Division, 1968,94:637-659.
[25] BANDIS S C, LUMSDEN A C, BARTON N R. Fundamentals of rock joint deformation[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1983,20:249-268.
[26] KUMAR A. The effect of stress rate and temperature on the strength of basalt and granite[J]. Geophysics, 1968,33(3):501-510.