Z 字形斷層系統剪切應變陡降時空分布的數值模擬*

王學濱 馬 冰 呂家慶 顧 路

1)遼寧工程技術大學力學與工程學院，阜新 123000

2)中國地震局地質研究所(地震動力學國家重點實驗室)，北京100029

1 引言

巖石類材料的微破裂與地震活動具有很大的類似性。目前，在實驗室中主要通過聲發射技術獲取微破裂的時空分布規律，并由此進行斷層失穩前兆、地震機理等相關重要問題的探索［1－3］。

數值模擬作為一種研究手段，在破裂失穩研究方面發揮著重要作用。目前的一些數值模擬通過引入非均勻性，模擬破壞單元的數目的演變規律或能量釋放的時空分布規律［4－8］，以比擬物理實驗中的聲發射或天然的微震活動。眾所周知，在這些數值模型中，只有破壞的單元才有聲發射和釋放能量。也就是說，統計發生破壞單元的信息，是一種“一對一”的關系。所謂的“一對一”是指，當一個單元發生破壞時，即輸入是1，或者聲發射累計數增加1，或者計及這個單元釋放的能量，即輸出也是1。如果通過研究，能建立一種“一對多”的關系，即輸入是1時，輸出的是多個單元的信息，那么失穩的前兆可能會更加明顯，更易于識別。

本文運用二次開發的FLAC-3D 軟件，對含4 條斷層的巖石標本進行數值模擬，研究標本內部剪切應變陡降的時空分布規律和兩種能量釋放的演變規律，以及斷層失穩造成的標本整體失穩的前兆及失穩過程。

2 計算模型及參數取值

研究模型高、寬均為0.3 m，被剖分成尺寸相同的9 萬個立方體單元，其中斷層單元3 312 個。模型中有4 條斷層，其中3 條構成Z 字形斷層(圖1)。斷層A、B、C 及D 與水平方向的夾角分別為30°、30°、30°及60°。

首先，在靜水壓力σ1=σ3=－5 MPa 條件下，對模型進行計算(圖1(a))，阻尼力由FLAC-3D自行施加，迭代2 萬步后，節點的最大失衡力已經足夠小，說明模型已經達到靜力平衡狀態。然后，在模型的左、右邊界施加足夠小的、相反的速度v=2.5 ×1010m/時步(圖1(b))，即在水平方向上進行準靜態拉伸位移控制加載(σ1=－5 MPa)，計算出加載端的平均應力。在FLAC-3D 中，應力為負值代表壓縮，反之代表拉伸;σ1≤σ3。一個時步是指一次計算循環，從運動方程經由幾何方程，到本構方程，再到運動方程。循環一次，則時步數目增加1。

計算中，斷層單元和巖石單元都采用帶拉伸截斷的應變軟化莫爾-庫侖模型［6，7］。彈性模量、粘聚力和抗拉強度都服從Weibull 分布;非均質性參數m=9。巖石單元和斷層單元的泊松比為0.25 和0.2，彈性模量的均值為55 GPa 和5.1 GPa，粘聚力的均值為37.5 MPa 和5 MPa，抗拉強度的均值為24 MPa 和1.2 MPa，初始內摩擦角為50°和10°，擴容角均取為0°［7］。

圖1 Z 字形斷層構造的計算模型Fig.1 A computational model for a Z-shaped fault

3 結果分析

3.1 標本破壞的5 個階段

圖2 平均應力-時步曲線Fig.2 Average stress-timestep curve(a)and average stresstimestep curve in the vicinity of the peak stress(b)

圖2 為速度加載端的平均應力隨時步的演變規律。由圖2可以發現，剛開始時，左、右加載端的平均應力為－5 MPa，稍后，平均應力上升比較快，致使平均應力-時步曲線呈上凸現象，隨后，該曲線的斜率趨于定值。

大約7.5 萬時步時，平均應力已由負值變為正值，表明標本左、右加載端的平均應力已成為拉應力。繼續加載，平均應力-時步曲線開始向下彎曲，即又呈現上凸現象，隨后，平均應力達到最高點。然后，平均應力開始下降。當時步達到14 萬時，平均應力下降開始變慢，即平均應力-時步曲線越來越平緩，表明標本已接近進入殘余變形階段，殘余應力為壓應力。這一點似乎是難于理解的。實際上，當水平方向的拉應力降至很低時(接近于零)，標本將在σ1的驅動下在水平方向上運動，標本內部的斷層相應地錯動，因而標本速度加載端的平均應力可以成為壓應力。

總體上，標本的平均應力-時步曲線包括5 個典型階段:

1)上凸階段。與位移控制加載有關;

2)直線階段，即線性穩態階段［9］。在這一階段發生位移強化及應力累積;

3)上凸階段，即偏離線性穩態階段［9］。預示局部應力開始釋放，總體上仍表現為應力累積階段;

4)軟化階段，為亞失穩階段和失穩滑動階段［9］;

5)殘余階段，即調整階段［9］。

圖2(b)為應力峰及附近放大結果。由圖2(b)可見，在應力峰之前，平均應力呈非線性上升特點，但有許多應力降，應力降的幅度比應力峰之后的小。

圖2(b)中，帶箭頭的虛線為平均應力突然下降的時刻，該時刻標志著標本整體失穩的開始，標本整體失穩是斷層失穩錯動造成的。

在應力峰至標本失穩開始時刻之間，總體上平均應力是下降的，但存在局部的上升和下降，應力降盡管比應力峰之前的大，但遠小于標本整體失穩之后的應力降。這一階段應對應于地震的短臨階段或亞失穩階段［9］。觀察短臨階段平均應力的變化規律可以發現，在大約前2/3 階段(亞失穩Ⅰ階段)，應力降的幅度不太大，盡管應力釋放，但并不占優;在后1/3 階段(亞失穩Ⅱ階段)，存在較大的應力降，表明應力釋放開始占優勢。

3.2 剪切應變陡降的時空分布規律

在8 ～22 萬時步共劃分70 個時段，每個時段持續2 000 個時步。圖3 為剪切應變陡降量在不同時段的分布規律，本文僅統計超過5 ×10－6的剪切應變降(或稱之為剪切應變陡降)。為了與上述結果進行對比，圖4 給出了拉伸應變能及剪切應變能釋放［6，7］的時空分布規律。限于篇幅，僅給出了有代表性的部分時段的結果。

第1 ～13 時段(圖3(a，b))，剪切應變陡降量較小，零星地分布在除斷層A 上靠近標本左端面的局部區段外的各條斷層上。在第1 時段時，標本所受的平均應力尚處于線性階段，在第13 時段時，平均應力早已進入非線性階段，即應力峰之前的上凸階段。

第15 ～19 時段(圖3(c ～e))，剪切應變陡降量明顯增大，主要集中在斷層A、C 上，特別是在斷層A與D 相交匯的位置。顯然，斷層A 與D 形成了一種拐折結構。拐折部位的破壞為其他部位的破壞創造了條件。由圖4(a)可以發現，在斷層A 上有一個單元釋放了較高的拉伸應變能;此外，在拐折部位可以發現一些單元釋放能量。在第19 時段時，標本所受平均應力已接近應力峰。

第21 時段(圖3(f))，標本的平均應力已處于應力峰附近，許多發生剪切應變陡降的單元云集于斷層A 與D 的拐折部位，這一位置即將完全破壞。

第23 時段(圖3(g))，標本的平均應力已處于短臨階段。此階段觀察到的各種現象應與標本整體失穩的前兆有關。發生剪切應變陡降的單元已發生在斷層B 上，這是此前少見的現象。中間斷層B 的破壞或錯動，導致了標本的整體失穩。由圖4(b)可以發現，在斷層C 上有一個單元釋放了非常高的剪切應變能;此外，在拐折部位和中間斷層B 上，都可以發現大量的單元釋放能量。

第25 ～29 時段(圖3(h ～j))，標本的平均應力已處于應力峰之后，剪切應變陡降比較密集地分布在各條斷層附近，這意味著各條斷層上的高強度單元被相繼錯斷，造成了標本整體的失穩。由圖4(c)可以發現，能量釋放遍布在各條斷層上。

第31 ～63 時段(圖3(k，l))，發生剪切應變陡降的單元變少，回歸平靜，零星地散布在各條斷層上，標本的平均應力處于緩慢的下降階段，隨著變形的繼續，逐漸進入殘余階段。能量釋放繼續發生(圖4(d))。

3.3 3 種量的統計結果分析

圖5是對釋放的剪切應變陡降、拉伸應變能及剪切應變能統計結果，分別統計這些量的累計(在任一時段之內，將所有單元的某種量求和)、最大值及發生剪切應變陡降或能量釋放的單元數目。在圖5(b，e，h)中，可以觀察到上述3 種量的最大值隨時段的演變規律。

圖3 不同時段內剪切應變陡降的時空分布(圓圈的半徑代表剪切應變陡降的大小，圓圈的中心標明了發生剪切應變陡降的單元位置)Fig.3 Spatiotemporal distribution of the shear strain with steep drop in different timestep intervals

圖4 不同時段彈性應變能釋放的時空分布(圓圈的半徑代表兩種能量釋放的大小，圓圈的中心標明為能量釋放單元的位置;黑色和紅色圓圈的半徑分別代表發生剪切和拉伸破壞單元能量釋放的大小)Fig.4 Spatiotemporal distribution of the liberated elastic strain energy in different timestep intervals

第1 ～13 時段，由圖5(a ～c)可以發現，與剪切陡降有關的3 種量隨著時段的增加，并無大的變化。但是，在與拉伸應變能及剪切應變能釋放有關的一些量上，卻能發現增加的現象(圖5(d ～g，i))。在圖5(h)中，發生拉破壞單元能量釋放的最大值處于較高的水平上，這應該與斷層單元的強度比巖石單元的低，先發生拉張破壞有關。

第15 ～19 時段，隨著時段的增加，很多量都在增加(圖5(a ～g，i))。圖5(h)的結果是個例外，在第18 ～19 時段內，發生拉張破壞單元能量釋放的最大值下降到較低的水平，表現為異常的平靜，盡管此時發生拉張破壞單元的數目及它們能量釋放的累計并不小。這說明，盡管事件比較多，但沒有大事件。

第21 時段，圖5(h)中的結果繼續表現為低值。

第23 時段，圖5(a，c ～f)可以觀察到與剪切應變陡降及能量釋放有關的量的加速現象;圖5(h)中，結果仍然處于低值;圖5(i)中，結果已經由上凸拐至水平。

圖5 9 種力學指標在70 個時段中的表現(虛線標明了標本整體失穩的開始)Fig.5 Evolution of nine mechanical parameters during 70 timestep intervals(The dashed line denotes the onset of the overall instability of the sample)

第25 ～29 時段，圖5(a ～c)中的3 個結果都表現為高值，特別是在圖5(a，c)中，兩種結果表現為全局的最高值。在圖5(d，e，g)中，各種結果也都表現為高值;在圖5(d，f，g)中，隨著時段的增加，各種結果上升的趨勢比較明顯，這與圖5(i)中結果的表現正好相反;在圖5(h)中，結果從低值有所恢復。從圖5(d，f，g，i)可以發現，第25 ～29 時段的結果在此前和此后的結果之間起到過渡的作用，即經過25 ～29 時段，各種結果從一種規律轉變成另一種規律，這反映了標本變形模式的轉變。

第31 ～63 時段，圖5(a，c)中的兩種結果降到比第1 時段還低的程度;圖5(d)中，結果在高水平上波動;圖5(f)中，結果不斷增加;圖5(g)中，結果與圖5(d)中的結果具有一定的類似性;圖5(e ～h)中，每種結果存在一定的波動，但小于全局的最大值;在圖5(i)中，結果維持不變。

總之，從一些與剪切應變陡降和能量釋放有關的結果的演變規律上看，在第25 時段稍前，能看到有的加速，有的拐平，有的處于低谷等異常，這說明標本整體失穩的前兆是有的，但表現各異，明顯程度也不相同。似乎，一些與剪切應變陡降有關的量具備一定的優勢，因為在失穩稍前及失穩過程中(加載端平均應力快速下降)，它們都處于高值，遠高于此前及此后的背景值。這有利于準確識別標本所處的應力狀態和整體失穩的時刻。

4 結論

1)對于由4 條斷層構成的Z 字形斷層系統，斷層A 與D 形成了拐折斷層，拐折部位的破壞及錯動為其他斷層的破壞及錯動創造了條件，兩條平行斷層A 與C 之間的斷層B 的破壞及錯動造成了標本整體的失穩。

2)在標本整體失穩稍前，能觀察到一些與剪切應變陡降和能量釋放有關的量發生了異常的變化，例如，加速、拐平及處于低谷。與剪切應變陡降有關的兩種量在標本整體失穩稍前及應力急劇下降過程中，處于高值，顯著地區別于此前及此后的低背景值。

3)與剪切應變陡降有關的一些量在失穩過程中表現突出的原因是由于一個發生破壞的單元，能引起其周圍大量未破壞單元的反應，因而更有利于斷層失穩前兆的識別，相當于對前兆起到了放大的作用。

4)本文的數值模擬研究是初步的，尚有許多問題值得探討，例如，不同的斷層系統在失穩過程中，一些與剪切應變陡降有關的量是否存在共性和差異;各條斷層的失穩是否也能通過剪切應變陡降反映出來;亞失穩階段的仔細考察，等等。

致謝衷心感謝導師馬瑾院士的指導!

1 來興平，等.基于AE 的煤巖破裂與動態失穩特征實驗及綜合分析［J］.西安科技大學學報，2006，26(3):289－292，305.(Lai Xingping，et al.Comprehensive characteristics analysis of crack propagation and dynamical destabilization of coal-rock based on acoustic emission experiment［J］.Journal of Xi’an University of Science and Technology，2006，26(3):289－292，305)

2 李爭榮，等.巖體失穩過程中的聲發射監測研究［J］.金屬礦山，2012，(1):28－32.(Li Zhengrong，et al.Study on AE monitoring for rockmass failure processing［J］.Metal Mine，2012，(1):28－32)

3 馬勝利，等.斷層階區對滑動行為影響的實驗研究［J］.中國科學D 輯:地球科學，2008，38(7):842－851.(Ma Shengli，et al.Effect of fault jogs on frictional behavior:An experimental study［J］.Science in China(Ser D)，2008，38(7):842－851)

4 Liu H Y，et al.Numerical studies on the failure process and associated microseismicity in rock under triaxial compression［J］.Tectonophysics，2004，384:149－174.

5 Fang Z and Harrison J P.Development of a local degradation approach to the modeling of brittle fracture in heterogeneous rocks［J］.Int J Rock Mech.Min.Sci.，2002，39(4):443－457.

6 Wang X B，Ma J and Liu L Q.Numerical simulation of failed zone propagation process and anomalies related to the released energy during a compressive jog intersection［J］.J Mech Mat Struct.，2010，5:1 007－1 022.

7 Wang X B，Ma J and Liu L Q.A comparison of mechanical behavior and frequency-energy relations for two kinds of echelon fault structures through numerical simulation［J］.Pure Appl Geophys，2012，169:1 927－1 945.

8 王學濱.缺陷數目對巖樣聲發射及應變能降低的影響［J］.中國有色金屬學報，2008，18(8):1 541－1 549.(Wang Xuebin.Influence of imperfection number on acoustic emissions and elastic strain energy decrease of rock specimens with initially random imperfections［J］.The Chinese Journal of Nonferrous Metals，2008，18(8):1 541－1 549)

9 馬瑾.地震新主體地區和失穩危險階段的探尋［J］.世界地震動態，2012，(6):231.(Ma Jin.Study on the new seismic major region and hazardously unstable steps［J］.Recent Development in World Seismology，2012，(6):231)