王海軍 ，任旭華，陶冉冉，張繼勛

(1.河海大學 水利水電學院，江蘇 南京，210098；

2.西澳大學 土木與資源工程學院，珀斯 6009，澳大利亞)

巖石記憶性是指巖石儲存外界對其影響信息并能在一定條件下通過某種物理量顯示這些信息的特性。巖石記憶性有多種，如變形記憶性[1]、聲發(fā)射記憶性(又稱Kaiser效應)、離子發(fā)射記憶性和導電率記憶性等[2]。變形記憶性是指可通過巖石變形數(shù)據(jù)獲取其記憶信息的特性，與Kaiser效應通過聲發(fā)射獲取記憶信息具有相似性。以單軸循環(huán)壓縮試驗為例，當后續(xù)加載應力超過先期加載應力最大值時，表現(xiàn)為應變曲線發(fā)生變化，由變化點得到先期加載最大應力信息，此為變形記憶性；若初始加載位于裂紋初始應力以上的應力區(qū)，應力達到先期加載最大應力時，聲發(fā)射現(xiàn)象明顯增加，此為Kaiser效應。地應力大小及方向的準確測量對于地震預報及各種地下工程的開挖設計具有重要意義[3]，初始地應力信息同樣儲存在巖芯中。變形率變化法(Deformation rate analysis，DRA)和聲發(fā)射法(Acoustic emission，AE)是分別基于巖石變形記憶性和Kaiser效應獲取地應力信息的方法。與ISRM推薦且應用最為廣泛的水力致裂法和套芯應力解除法不同，DRA法及AE法通過在室內(nèi)試驗室對巖芯試樣進行循環(huán)單軸壓縮試驗實現(xiàn)地應力測量。巖石試樣可來自前期探洞開挖積累下的巖塊或者巖芯。具有簡易、經(jīng)濟、高效率的優(yōu)勢；且可進行大量的實驗，用于統(tǒng)計分析以提高結果的精確度；同時，不受巖石參數(shù)各向異性的影響。巖石變形記憶性及DRA法被提出后，在日本[4-9]、俄羅斯[2]、澳大利亞[10-12]和中國臺灣[13-14]有大量研究與應用；韓國、中國大陸[15]有少量研究。其中絕大部分研究重點為其物理現(xiàn)象特征或地應力測量的應用，涉及其形成機理的研究較少。Yamamoto等[1]在1990年提出DRA法時推斷巖石變形記憶性的力學機理與Kaiser效應一樣，為巖石內(nèi)部新微裂紋的產(chǎn)生和原有微裂紋擴展。他們認為微裂紋擴展所造成的非線性應變包含可逆應變和不可逆應變2部分。在不改變微裂紋尺寸的情況下所產(chǎn)生的非線性應變都屬于可逆應變。在連續(xù)加載中，應力超過初始加載最大應力時，新裂紋產(chǎn)生，已有裂紋擴展，所造成的應變?yōu)椴豢赡鎽?，進而形成變形記憶性。部分研究者[16]沿用了此描述。基于Kuwahara等[17]的剪切裂紋理論模型，Tamaki等[18]指出當加載應力低于歷史應力最大值時，非彈性應變與應力呈線性關系。當加載應力超過該方向歷史應力最大值時，應力-應變曲線的斜率將產(chǎn)生變化，造成變形記憶性。Tamaki等[19]和Yamamoto等[5]采用了與Tamaki等[18]類似的解釋。Yamamoto[6]于2009年進一步指出當加載應力大于該方向初始加載應力最大值時，應變差對應力的微分值變?yōu)樨撝?，進而產(chǎn)生巖石變形記憶性，并認為其原因為：(1)當超過歷史最大應力值時，微裂紋將產(chǎn)生；(2)新微裂紋的產(chǎn)生會造成非彈性應變率增加。Hunt等[10,20]認為DRA折點是Kaiser效應的一個表現(xiàn)。他們基于PFC2D軟件，對Kaiser效應和巖石變形記憶性進行了模擬。基于模擬結果，確認Kaiser效應與微裂紋的產(chǎn)生和擴展有關，指出當初始加載應力小于微裂紋初始應力時，巖石變形記憶性并不存在。由以上研究可知：研究者基本遵從一種模式來解釋巖石變形記憶性的形成機理，伴隨加載應力增高，微裂紋的產(chǎn)生及擴展將造成非彈性應變率增加并表現(xiàn)在應力-應變曲線的變化上，從而形成巖石變形記憶性。但是，大量研究[1,2,10,13]表明：巖石變形記憶性在低于裂紋初始應力的低應力區(qū)仍然存在。在此低應力區(qū)域，微裂紋產(chǎn)生及擴展模型并不起作用。對物理現(xiàn)象內(nèi)在機理的正確把握是人們理解及利用這一外在現(xiàn)象的基礎。變形記憶性機理研究的缺失嚴重阻礙著人們對巖石變形記憶性的正確認識及對DRA法應用技術的改進與推廣。本文作者認為微裂紋的產(chǎn)生及擴展前的裂紋及顆粒接觸面間的摩擦滑動為低應力區(qū)變形記憶性的形成機理。為驗證此機理，首先綜述了裂紋初始應力值的相關研究及已有研究中低應力區(qū)存在巖石記憶性的證據(jù)，同時，通過物理實驗驗證了低應力區(qū)巖石記憶性的存在性?；诮佑|面摩擦滑動理論，采用Hook體和圣韋南體(St.V體)組合建立了單接觸面上摩擦滑動的基本單元；最終建立含有多裂紋及顆粒接觸面的理論模型，并分析其力學特性及巖石變形記憶性的表現(xiàn)，進而驗證摩擦滑動可以產(chǎn)生低應力區(qū)巖石記憶性。

1 DRA法及低應力區(qū)變形記憶性

1.1 DRA法

巖石所儲存信息一般是通過實驗室內(nèi)單軸壓縮實驗獲得。最直接方法是檢測單軸壓縮實驗的應力-應變曲線的斜率變化，但是識別難度很高，常不可靠[1]。為解決此問題，Yamamoto等[1]于1990年提出了變形率變化法。首先定義單軸循環(huán)壓縮實驗中的應變差為：

式中：εi(σ)和 εj(σ)分別代表第i次和第j次加載中的軸應變；σ為相應的應力。壓應力及壓應變?yōu)檎?。?1)消除了連續(xù)2次壓縮應變曲線中的可逆應變部分，得到非可逆應變差。應變差及變形率變化法描述如圖1所示。從圖1可見：應變差-應力曲線(DRA曲線)出現(xiàn)1個折點(DRA折點)，對應的應力σDRA即為試樣在加載方向上所記憶的歷史最高應力。

圖1 應變差及變形率變化法Fig.1 Illustration of deformation rate analysis (DRA)

采用DRA法進行地應力測量時，一個基本假設是σDRA對應于試樣軸向方向的正應力分量。因此，只需測出同一區(qū)域巖芯6個不同方向的正應力分量，即可計算出地應力狀態(tài)。

地應力真實值一般未知，因此，在研究巖石記憶性時，研究者一般采用一次預加載(多為單軸壓縮)以形成可知的記憶信息(或?qū)Φ貞M行模擬)。此時，DRA法所得σDRA對應軸方向預加載的最大應力。

1.2 低應力區(qū)變形記憶性

巖石在單軸壓縮下的變形及破壞過程在近幾十年得到大量研究。一般認為，單軸壓縮下巖石應力-應變曲線可分為5個階段[21-22]：(1)裂紋閉合階段；(2)線彈性變形階段；(3)裂紋產(chǎn)生及穩(wěn)定擴展階段；(4)裂紋非穩(wěn)定擴展階段；(5)破壞或峰后軟化階段。(2)和(3)階段的分界點稱為裂紋初始應力。研究表明：裂紋初始應力一般為試樣單軸抗壓強度(Unconfined compression strength,σUCS)的30%～60%[21-23]。在低于裂紋初始應力的應力區(qū)，并不發(fā)生微裂紋的產(chǎn)生及擴展，此區(qū)域為低應力區(qū)。

如前所述，大部分文獻將巖石變形記憶性的形成歸因于新裂紋的產(chǎn)生及已有裂紋擴展。若此模型成立，與Kaiser效應一樣，巖石變形記憶性應該存在于階段(3)和(4)中，即至少大于30%σUCS的應力區(qū)間。但是，實驗表明：在低應力區(qū)仍然存在巖石變形記憶性。

Yamamoto等[1]證實DRA法可以從花崗閃長巖的巖芯試樣中測出1～6 MPa的地應力；Seto等[16]的實驗結果表明：當歷史最大應力小于單軸抗壓強度的15%時，巖石仍然具有變形記憶性。Hunt等[10,13]的實驗表明：巖石記憶性存在于低于單軸抗壓強度的20%的應力區(qū)域。Yamshchikov等[2]指出巖石變形記憶性在彈性變形區(qū)仍然存在。

1.3 實驗驗證

試樣材料為砂巖，取自于澳大利亞北部一礦場76 mm標準巖芯。試樣為圓柱體，上下表面平整度為0.01 mm，如圖2所示。試樣從巖芯鉆取后，在常溫條件下放置時間為7 d。試樣密度為2.85 t/m3, 長度為39.9 mm, 直徑為18.3 mm，彈性模量為44 GPa，σUCS大于80 MPa。采用4對應變片采集應變數(shù)據(jù)。4對應變片在圓柱試樣側表面按照0°，90°，180°和270°的順序黏貼在試樣中間。加載過程由伺服加載控制系統(tǒng)控制，采用位移加載控制模式，為0.14 mm/min，加卸載速度相同。加載及數(shù)據(jù)采集系統(tǒng)如圖2所示。

圖2 砂巖試樣及加載儀器Fig.2 Sandstone sample and loading machine

圖3所示為實驗的加載方式。從圖3可以看出：實驗采用三次循環(huán)單軸壓縮。預加載最高應力σp為8.3 MPa(＜10.4%σUCS)，用來形成試樣記憶信息。第2和3次稱為測量加載，為DRA法中的Δε2,3(σ)提供應力-應變數(shù)據(jù)，用來獲取試樣對σp的記憶信息。測量加載一般為連續(xù)加載，最高應力σm相同，本實驗中為75 MPa。

應力-應變曲線如圖4所示，從應力-應變曲線上很難識別記憶信息。DRA曲線如圖5所示。從圖5可見：在8.3 MPa附近有清晰折點(箭頭標注)，由此證明在歷史應力小于10.4%σUCS的應力區(qū)域內(nèi)，巖石仍然有變形記憶性。

圖3 加載方式Fig.3 Loading cycles

圖4 砂巖試樣的應力-應變曲線Fig.4 Strain-stress curve of sandstone sample

圖5 砂巖試樣的DRA曲線Fig.5 DRA curve of sandstone sample

2 基于摩擦滑動的理論模型

2.1 接觸面的摩擦滑動

巖石是一種含有大量隨機分布微裂紋接觸面及顆粒接觸面的材料。考慮一個單位體積的此類材料的平均應力-應變狀態(tài)[24]為：

式中：σik為單位體積所受應力張量；ε0ik和σ0ik分別為彈性基質(zhì)的應變與應力張量，并不受微裂紋及顆粒接觸面的影響；ni為微裂紋及顆粒接觸面外法線單位矢量；Vi是沿接觸面的應變積分；α代表微裂紋或顆粒接觸面序號。

由式(2)可知：含有微裂紋及顆粒接觸面的彈性材料應變包含2部分：一部分為彈性基質(zhì)產(chǎn)生，另一部分為接觸面產(chǎn)生。影響第2部分的力學行為的因素很多，如微裂紋的方向和尺寸等。在低應力區(qū)，在加載作用下，裂紋及顆粒接觸面產(chǎn)生摩擦滑動，并不造成微裂紋的產(chǎn)生與擴展[25]。本研究將低應力區(qū)巖石變形記憶性的產(chǎn)生機理歸為摩擦滑動?；谝陨侠碚摚捎肏ook體及St.V體的組合構建模型，實現(xiàn)對裂紋及顆粒接觸面摩擦滑動的模擬。

2.2 基本單元模型

圖6所示為含有單接觸面的單位體積巖石的基本單元模型。其包括串聯(lián)的2部分：一是為上部Hook體，用來模擬彈性基質(zhì)的變形；二是Hook體與St.V體并聯(lián)，用于模擬接觸面的摩擦滑動對總應變的貢獻。為了簡便，本研究用“||”表示并聯(lián)，第2部分稱之為“Spr||St.V”體。在“Spr||St.V”體中，St.V體用來模擬摩擦滑動，當其承受應力超過黏聚力時，將開始滑動并保持應力值不變。k2為頂部串聯(lián)Hook體剛度；k1為并聯(lián)Hook體剛度；c1為St.V體的黏聚力。

圖6 基本單元模型Fig.6 Basic element

Hook體 2 與“Spr||St.V”體串聯(lián)，因此，兩者應力相同，模型應變?yōu)閮烧邞冎停?/p>

對于St.V體，有靜止和滑動2種狀態(tài)：滑動狀態(tài)：

靜止狀態(tài)：

式中：σ為加載應力；ε為模型總應變；σ2和ε3分別為Hook體2的應力和應變；σ1和ε1分別為“Spr||St.V”體的應力和應變。σ1_fric是St.V體的應力，當|σ1_fric|＜c1時，“Spr||St.V”體被St.V體鎖定，保持初始應變ε1(0)不變。

數(shù)值試驗的加載方式如圖3所示，加載應力為時間的函數(shù)。預加載最高應力值σp，用以形成記憶信息；后2次為測量加載，最大值為σm(σm＞σp)。2次測量加載的應變差Δε2,3(σ)曲線為DRA曲線。

數(shù)值實驗結果表明：基本單元中Hook體并不影響其變形記憶性，隨著加載區(qū)域的不同，基本單元的DRA曲線最多可出現(xiàn)2個折點，其對應的應力值按其出現(xiàn)順序記為σDRA1和σDRA2，如圖7和8所示。從圖7和圖8可見：

(1)當σp＜σm＜c1時，“Spr||St.V”體未滑動，無DRA折點，模型無記憶性；當2c1＜σp＜σm時，2次測量加載應力-應變曲線為完全重合的循環(huán)加卸載滯回環(huán)，此時也沒有DRA折點。

(2)當c1＜σp＜σm＜2c1時，如圖7 所示，DRA曲線有1個折點，此折點對應應力值σDRA1=σp。在此加載區(qū)域里，DRA法精確測出了預加載應力值。

(3)當c1＜σp＜2c1＜σm時，DRA曲線如圖8所示，DRA曲線存在2個折點：第1個折點為σDRA1=σp，是對預加載應力的記憶；第2個折點為σDRA2=2c1，是對2倍黏聚力的記憶，此折點與變形記憶性無關。

圖7 σp/c1=1.5和σm/c1=1.8時的DRA曲線Fig.7 DRA curve under condition of σp/c1=1.5 and σm/c1=1.8

圖8 σp/c1=1.5和σm/c1=3時的DRA曲線Fig.8 DRA curve under condition of σp/c1=1.5 and σm/c1=3

2.3 多接觸面的摩擦滑動模型

含有n個裂紋及顆粒接觸面的單位體積巖石材料的理論模型如圖9所示。不考慮接觸面之間的相互影響，由式(2)可知：對基本單元采用串聯(lián)方式組合進行含有多裂紋及顆粒接觸面巖石材料的模擬：

圖9 含多接觸面的理論模型Fig.9 Model for rock with many interfaces

式中σ為施加應力；α為接觸面序列數(shù)。每個接觸面的力學性質(zhì)符合式(3a)～(3f)。

在數(shù)值分析中，選擇n=500。由基本單元可知：黏聚力對DRA折點有控制作用，本研究選用均勻分布與正態(tài)分布2種分布的黏聚力。其均值和標準差設置相同，用以比較2種分布對結果的影響。多接觸面模型的其余參數(shù)都相同。加載方式如圖3所示。黏聚力均勻分布的范圍為0～5 MPa；2種分布均值為2.505 MPa，標準差為1.444 8 MPa；正態(tài)分布如圖10所示。

典型應力-應變曲線如圖11(a)和12(a)所示，所對應的DRA曲線如圖11(b)和12(b)所示。結果如下：

(1)多裂紋及顆粒接觸面模型可形成下折的DRA曲線，即可產(chǎn)生變形記憶性。同時，黏聚力的分布范圍的增大拓寬了模型變形記憶性的可記憶范圍，DRA折點存在于最小黏聚力到2倍最大黏聚力之間。在此區(qū)間，DRA法可精確測出預加載應力。

(2)與基本單元模型相比，多接觸面模型中由于接觸面數(shù)目增加，DRA曲線在折點前為水平線，下折后的曲線部分呈現(xiàn)光滑下凹的特性，此特性和試驗中DRA曲線特性相吻合(圖5)。

(3)黏聚力的正態(tài)分布及均勻分布并不影響DRA曲線的形狀，如圖11(b)和12(b)所示。

圖11 黏聚力均勻分布下的結果Fig.11 Results under uniformly distributed cohesions

圖12 黏聚力為正態(tài)分布時的結果Fig.12 Results under normally distributed cohesions

3 結論

(1)巖石微裂紋產(chǎn)生及擴展并不能解釋低應力巖石變形記憶性，本研究給出低應力區(qū)巖石存在變形記憶性的實驗證據(jù)，并提出了巖石內(nèi)部原有微裂紋及顆粒接觸面上的摩擦滑動為低應力區(qū)形成變形記憶性的機理，揭示了巖石對于所承受加載應力包括地應力具有普遍記憶性的深層原因。

(2)基于摩擦滑動理論，構建了模擬單接觸面摩擦滑動的基本單元。結果表明：在一定條件下，基本單元可以形成變形記憶性，記憶范圍位于1～2倍黏聚力之間的應力區(qū)域內(nèi)。

(3)基于基本單元，構建多接觸面理論模型，結果證明模型可形成記憶性。DRA曲線在DRA折點前為水平線，折點后為光滑下凹曲線，與物理實驗中DRA曲線特性吻合。同時，黏聚力不同分布規(guī)律并不影響DRA曲線形狀。

(4)低應力區(qū)微裂紋及顆粒接觸面的摩擦滑動可以形成巖石變形記憶性。下一步工作將在模型中納入時間因子，以研究變形記憶性的時間效應如失憶性。

致謝：

本文作者感謝國家留學基金委提供國家建設高水平大學公派項目獎學金，感謝澳大利亞地質(zhì)力學中心(Australian Centre for Geomechanics)提供獎學金支持作者在西澳大學完成相關研究。驗證實驗由西澳大學Ariel Hsieh完成，數(shù)據(jù)由西澳大學Ariel Hsieh，Arcady Dyskin和Phil Dight提供，在此表示感謝。

[1]Yamamoto K, Kuwahara Y, Kato N, et al.Deformation rate analysis: A new method for in situ stress estimation from inelastic deformation of rock samples under uni-axial compressions[J].Tohoku Geophysical Journal, 1990, 33:127-147.

[2]Yamshchikov V S, Shkuratnik V L, Lavrov A V.Memory effects in rocks (review)[J].Journal of Mining Science, 1994, 30(5):463-473.

[3]Fairhurst C.Stress estimation in rock: A brief history and review[J].International Journal of Rock Mechanics and Mining Sciences, 2003, 40(7/8): 957-973.

[4]Yabe Y, Yamamoto K, Sato N, et al.Comparison of stress state around the Atera fault, central Japan, estimated using boring core samples and by improved hydraulic fracture tests[J].Earth,Planets, and Space, 2010, 62: 257-268.

[5]Yamamoto K, Yabe Y.Stresses at sites close to the Nojima Fault measured from core samples[J].Island Arc, 2001, 10(3/4):266-281.

[6]Yamamoto K.A theory of rock core-based methods for in-situ stress measurement[J].Earth Planets Space, 2009, 61(10):1143-1161.

[7]Yabe Y, Song S R, Wang C Y.In-situ stress at the northern portion of the Chelungpu fault,Taiwan, estimated on boring cores recovered from a 2-km-deep hole of TCDP[J].Earth Planets Space, 2008, 60(8): 809-819.

[8]Yabe Y, Omura K.In-situ stress at a site close proximity to the Gofukuji Fault, central Japan, measured using drilling cores[J].Island Arc, 2011, 20(2): 160-173.

[9]Seto M, Villaescusa E, Utagawa M, et al.In situ stress evaluation from rock cores using AE method and DRA[J].Shigen-to-Sozai,1998, 114(12): 845-855.

[10]Hunt S P, Meyers A G, Louchnikov V.Modelling the Kaiser effect and deformation rate analysis in sandstone using the discrete element method[J].Computers and Geotechnics, 2003,30(7): 611-621.

[11]Dight P M, Dyskin A V.On the determination of rock anisotropy for stress measurements[C]//SHIRMS 2008.Perth, Australia,2008: 575-585.

[12]Dight P.Stress states in open pits.Keynote lecture-Slope Stability in Mining and Civil Engineering(Report)[R].2011.Perth.

[13]Chan S C.Investigating the laboratory experiments to estimate pre-stress on Changchikeng sandstone[D].Tainan: National Cheng Kung University, 2008: 57-85.

[14]Wu J H, Jan S C.Experimental validation of core-based pre-stress evaluations in rock: A case study of Changchikeng sandstone in the Tseng-wen reservoir transbasin water tunnel[J].Bulletin of Engineering Geology and the Environment, 2010,69(4): 549-559.

[15]謝強, 邱鵬, 余賢斌, 等.利用聲發(fā)射法和變形率變化法聯(lián)合測定地應力[J].煤炭學報, 2010, 35(4): 559-564.XIE Qiang, QIU Peng, YU Xian-bin, et al.Initial-stress measurements with AE and DRA combined technique[J].Journal of China Coal Society, 2010, 35(4): 559-564.

[16]Seto M, Nag D K, Vutukuri V S.In-situ rock stress measurement from rock cores using the acoustic emission method and deformation rate analysis[J].Geotechnical and Geological Engineering, 1999, 17(3): 241-266.

[17]Kuwahara Y, Yamamoto K, Hirasawa T.An experimental andtheoretical study of inelastic deformation of brittle rocks under cyclic uniaxial loading[J].Tohoku Geophys J, 1990, 33(1): 1-21.

[18]Tamaki K, Yamamoto K, Furuta T, et al.An experiment of in situ stress estimation on basaltic rock core samples from Hole 758A,Ninetyeast Ridge, Indian Ocean[C]//Proc ODP, 121: College Station.TX(Ocean Drilling Program), 1991: 697-717.

[19]Tamaki K, Yamamoto K.Estimating in situ stress field from basaltic rock core samples of Hole 794C, Yamato Basin, Japan Sea[C]//Proceedings of the Ocean Drilling Program.San Francisco, 1992: 1047-1059.

[20]REN Xu-hua, WANG Ha-jun, ZHANG Ji-xun.Numerical study of the AE and DRA method in sandstone and granite in orthogonal loading directions[J].Water Science and Engineering,2012, 5(1): 93-104.

[21]Martin C D.The strength of massive Lac du Bonnet granite around underground openings[D].Manitoba: University of Manitoba, 1993.

[22]Eberhardt E, Stead D, Stimpson B, et al.Identifying crack initiation and propagation thresholds in brittle rock[J].Canadian Geotechnical Journal, 1998, 35(2): 222-233.

[23]Alkan H, Cinar Y, Pusch G.Rock salt dilatancy boundary from combined acoustic emission and triaxial compression tests[J].International Journal of Rock Mechanics and Mining Sciences,2007, 44(1): 108-119.

[24]Salganik R L.Overall effects due to cracks and crack-like defects[C]//Proceedings of First International Symposium on Defects and Fracture.Tuczno, Poland: 1982: 199-208.

[25]David E C, Brantut N, Schubne A, et al.Sliding crack model for nonlinearity and hysteresis in the uniaxial stress–strain curve of rock[J].International Journal of Rock Mechanics & Mining Sciences, 2012, 52: 9-17.