王學(xué)濱馬 冰呂家慶

1）中國地震局地質(zhì)研究所，地震動力學(xué)國家重點實驗室，北京 100029

2）遼寧工程技術(shù)大學(xué)力學(xué)與工程學(xué)院，阜新 123000

3）遼寧工程技術(shù)大學(xué)計算力學(xué)研究所，阜新 123000

實驗室尺度典型斷層系統(tǒng)破壞、前兆及粘滑過程數(shù)值模擬

王學(xué)濱1，2，3）馬 冰2）呂家慶2）

1）中國地震局地質(zhì)研究所，地震動力學(xué)國家重點實驗室，北京 100029

2）遼寧工程技術(shù)大學(xué)力學(xué)與工程學(xué)院，阜新 123000

3）遼寧工程技術(shù)大學(xué)計算力學(xué)研究所，阜新 123000

盡管實驗室尺度典型斷層模型的基本物理、力學(xué)參數(shù)已知，載荷、邊界條件明確，但其變形、破壞、前兆及穩(wěn)定性的綜合數(shù)值模擬并不容易。實驗室尺度模型的數(shù)值模擬可為大尺度模型的模擬積累經(jīng)驗，并提供某些必要的技術(shù)基礎(chǔ)。文中主要從數(shù)值模型的角度，評述實驗室尺度模型中引入非均質(zhì)性和粘滑過程模擬的各種方法，著重介紹了2種模型：斷層-巖石系統(tǒng)的非均質(zhì)應(yīng)變?nèi)趸Ｐ秃湍Σ翉娀?摩擦弱化模型。在第1種模型中，提出利用事件的頻次-能量釋放關(guān)系的斜率絕對值的不同演變規(guī)律評價斷層的不同應(yīng)力狀態(tài)。在擠壓雁列區(qū)貫通過程中，剪切應(yīng)變降對巖石破裂具有更為靈敏的指示作用，可能比聲發(fā)射更為有效。在第2種模型中，通過控制內(nèi)摩擦角在粘著和滑動兩階段不同的演變規(guī)律實現(xiàn)斷層粘滑過程的模擬。盡管在提出的摩擦強化-摩擦弱化模型中未引入速度，但由于對節(jié)點的運動方程進行求解，所以節(jié)點速度在粘著和滑動兩階段將發(fā)生變化。指出以下需要進一步研究的問題：1）亟待發(fā)展適于地質(zhì)體材料變形、破壞及運動整個過程模擬的計算力學(xué)模型；2）亟待發(fā)展復(fù)雜斷層系統(tǒng)中斷層穩(wěn)定性的可靠判別方法；3）亟待發(fā)展動力學(xué)方程的高效求解方法。

實驗室尺度模型 破壞 粘滑 前兆 摩擦 剪切應(yīng)變 雁列構(gòu)造

0 引言

近年來，隨著計算機軟、硬件技術(shù)水平的不斷發(fā)展，以有限元方法為代表的數(shù)值計算方法在地震研究中的應(yīng)用愈發(fā)得到重視。在過去很長一段時期內(nèi)，盡管數(shù)值方法在地殼應(yīng)力、應(yīng)變場研究方面和地震前兆、遷移、危險區(qū)預(yù)測研究方面發(fā)揮了積極的作用，但遠未像野外觀測、室內(nèi)實驗?zāi)菢拥玫矫芮嘘P(guān)注，數(shù)值計算一般聲譽較低?？陀^地講，數(shù)值計算具有其他研究手段不具備的先天優(yōu)勢，例如：節(jié)約人力、物力；正演、反演均可；可以復(fù)現(xiàn)過去，預(yù)測未來；各種參數(shù)提取方便，可以方便地開展參數(shù)的敏感性研究；可以模擬在當(dāng)前技術(shù)條件下難以開展的物理實驗；可以檢驗有關(guān)的假說，闡明復(fù)雜現(xiàn)象發(fā)生的機理及過程，等等。因此，人們探索及應(yīng)用數(shù)值計算方法的腳步從未停止。至今，尚未能找到量化預(yù)測強烈地震發(fā)生的運動學(xué)和動力學(xué)模型，也尚未發(fā)現(xiàn)控制其孕育、發(fā)生和發(fā)展的標(biāo)志性定量要素以及邊界條件和約束（滕吉文，2010）。因此，應(yīng)該倡導(dǎo)在多要素約束下提取初始模型進行數(shù)值模擬，為此在綜合分析的基礎(chǔ)上可試圖提出震級較強地震、強烈地震和大地震孕育和發(fā)生的強度（滕吉文，2010）。實現(xiàn)上述目的，不能單靠數(shù)值模擬，先進的數(shù)值計算方法應(yīng)該與詳實的觀測數(shù)據(jù)、地質(zhì)資料、巖石力學(xué)實驗、相關(guān)的理論模型、計算機技術(shù)等有機地結(jié)合。最近，國外學(xué)者對主震及同震、震間、震后現(xiàn)象、短周期慢速滑動及其發(fā)生間隔、長期和短期現(xiàn)象等進行了成功的數(shù)值模擬（Barbot et al．，2012；Shibazaki et al．，2012；Noda et al．，2013），給地震預(yù)測帶來了新的希望。

在地震的數(shù)值模擬研究中，常針對具體的大尺度問題進行建模和分析（陶瑋等，2000；鄭洪偉等，2006；李玉江等，2009；鄧志輝等，2011a，b；郭婷婷等，2011）。相比之下，少有研究人員針對實驗室尺度斷層標(biāo)本開展建模及計算分析工作。相反，從事力學(xué)、巖土力學(xué)的研究人員更熱衷于此。實際上，盡管實驗室尺度模型的基本物理、力學(xué)參數(shù)已知，邊界條件及加載條件明確，但其變形、破壞和穩(wěn)定性的綜合模擬并不容易。物理實驗中將發(fā)生一系列復(fù)雜的物理、力學(xué)現(xiàn)象，例如：聲發(fā)射，破裂擴展，斷層粘滑、穩(wěn)滑，變形局部化，紅外輻射及回彈等。一次伴隨著較大聲響的粘滑被認(rèn)為是一次實驗室地震。采用計算力學(xué)方法，部分地模擬上述復(fù)雜現(xiàn)象，是一項有意義的工作。實驗室尺度模型的模擬有助于充分了解計算模型的特點及優(yōu)劣；有助于深刻認(rèn)識實驗現(xiàn)象，支持或推翻一些假說和推測；有助于為大尺度模型的成功模擬積累經(jīng)驗，并提供某些必要的技術(shù)基礎(chǔ)；有助于發(fā)現(xiàn)新的力學(xué)量，從而探索地震預(yù)測的新途徑。這些力學(xué)量即使在現(xiàn)有技術(shù)條件下難以從實驗或野外觀測到，若能在不同的計算模型中得到反復(fù)的檢驗和確定，也將有廣闊的應(yīng)用前景。先在實驗室尺度模型的數(shù)值模擬方面取得若干進展，再將該模型進行適當(dāng)?shù)耐茝V和完善，開展大尺度問題的研究，是有益的嘗試。

本文主要從數(shù)值模型的角度，評述實驗室尺度模型中引入非均質(zhì)性和粘滑過程模擬的各種方法，著重介紹了2種模型：斷層-巖石系統(tǒng)的非均質(zhì)應(yīng)變?nèi)趸Ｐ停╓ang et al．，2010，2012，2013a）和摩擦強化-摩擦弱化模型（Wang et al．，2013b）。第1種模型適于模擬破壞及前兆，提出了事件的頻次-能量釋放關(guān)系的斜率絕對值（b0值）的計算方法，以2類雁列構(gòu)造為例進行了計算；在擠壓雁列區(qū)貫通過程中，發(fā)現(xiàn)了比聲發(fā)射和能量釋放的表現(xiàn)更敏感的力學(xué)指標(biāo)——剪切應(yīng)變陡降。第2種模型適于模擬粘滑過程，并得到了一個雙軸壓縮粘滑實驗的驗證。文章還指出了亟待研究的3個問題。

1 典型斷層系統(tǒng)破裂及前兆的數(shù)值模擬

1.1 斷層-巖石系統(tǒng)的非均質(zhì)應(yīng)變?nèi)趸Ｐ?/p>

地震活動具有明顯的區(qū)域性特征，在不同的觀察尺度上，大至地震帶，小到某個小震群的空間分布，無不顯示出地震活動在空間上的不均勻性（時振梁等，1995）。馬勝利等（2003）總結(jié)了中國實驗巖石力學(xué)與構(gòu)造物理學(xué)研究的若干新進展，認(rèn)為介質(zhì)和結(jié)構(gòu)的非均勻性以及時間過程的復(fù)雜性仍將是地震學(xué)和地球內(nèi)部物理學(xué)研究中的核心問題之一。劉力強等（1999）認(rèn)為，巖石構(gòu)造的非均勻程度可能是造成地震活動非均勻的重要因素之一。馬勝利等（2002）認(rèn)為，地震成核相的產(chǎn)生與斷層的非均勻性有關(guān)，發(fā)震斷層的非均勻性造成了地震破裂過程的非均勻性反應(yīng)。馬勝利等（2004）及Lei等（2004）發(fā)現(xiàn)，巖石的非均勻結(jié)構(gòu)是支配脆性巖石內(nèi)斷層形成過程的最重要因素；具有一定非均勻性的巖石在其破壞之前可以出現(xiàn)前兆性的異常微破裂活動特征。

在實驗室尺度模型中，考慮非均質(zhì)性的常見手段包括：1）引入細小單元力學(xué)參數(shù)的非均勻分布，例如：常見的Weibull和Gauss分布等（Cundall et al．，1979；Hobbs et al．，1989；唐春安等，1997；Fang et al．，2002；Hu et al．，2004；Liu et al．，2004）；2）考慮巖石的細觀結(jié)構(gòu)（Chen et al．，2004）；3）引入隨機材料缺陷（王學(xué)濱，2007）。在斷層力學(xué)和地震研究中，不同尺度的非均勻性的重要性怎么強調(diào)也不過分。研究表明，通過在連續(xù)介質(zhì)模型中引入某種程度上的空間非均勻性分布，即可呈現(xiàn)時空復(fù)雜的地震行為（Miyatake，1992；Bizzarri et al．，2001；Hillers et al．，2006）。甚至有人認(rèn)為，通過引入一個強烈的空間非均勻分布和強度自愈機制，速率和狀態(tài)依賴摩擦定律都不必考慮（Beroza et al．，1996）。在大尺度條件下，考慮一個塊體內(nèi)部細小單元力學(xué)參數(shù)非均勻性的意義不大，而考慮斷層的幾何非均勻性和摩擦系數(shù)等參數(shù)的非均勻性的意義更大。但是，對于實驗室尺度模型，可能并非如此。由于實驗室尺度模型需要模擬單元的破壞、破壞區(qū)或裂紋的擴展、聲發(fā)射、能量釋放及破壞前兆等，所以和考慮巖石的真實細觀結(jié)構(gòu)相比，考慮細小單元的非均質(zhì)性是一個有益的選擇。

FLAC-3D（3維連續(xù)介質(zhì)快速拉格朗日分析）是美國Itasca公司開發(fā)的一個商業(yè)軟件，在巖土力學(xué)、構(gòu)造地質(zhì)等領(lǐng)域具有世界聲譽和廣泛應(yīng)用（Strayer et al.，1997；Mckinnon et al.，1998；Erickson et al.，2001）。我們在FLAC-3D中引入了細小單元力學(xué)參數(shù)的Weibull分布，認(rèn)為斷層和巖石單元都是非均質(zhì)的（Wang et al．，2010，2012，2013a）。2種單元的本構(gòu)模型相同，但有些參數(shù)的取值不同。在彈性階段，2種單元均服從各向同性的線彈性模型；在破壞階段，2種單元均服從帶拉伸截斷的應(yīng)變?nèi)趸獱?庫侖模型。僅假定一部分材料參數(shù)（彈性模量、初始粘聚力和抗拉強度）是非均勻的，且它們彼此之間不相關(guān)，非均質(zhì)程度適中。其余參數(shù)被認(rèn)為是均勻的，例如泊松比、初始內(nèi)摩擦角等。

為了模擬應(yīng)變?nèi)趸?，除了需要?guī)定強度參數(shù)的初始值，還需規(guī)定弱化階段結(jié)束時的殘余強度參數(shù)和所對應(yīng)的塑性剪切或拉伸應(yīng)變的值。我們假定，隨著單元塑性變形的增加，單元的強度參數(shù)傾向于相同。也就是說，具有不同初始強度參數(shù)的單元的殘余強度參數(shù)及所對應(yīng)的塑性應(yīng)變相同。這也意味著，強度參數(shù)越高的單元峰后脆性越強。

研究發(fā)現(xiàn)，對于擠壓雁列構(gòu)造，破壞區(qū)除了會貫通雁列區(qū)，還會向外擴展（圖1a），而拉張雁列區(qū)由于斷層相互作用引起的破壞區(qū)僅位于雁列區(qū)內(nèi)部（圖1b），遠離拉張雁列區(qū)之外的破壞區(qū)和斷層的相互作用無關(guān)。上述結(jié)果與有關(guān)的野外觀測（Sibson，1985）和實驗（蔣海昆等，2002；馬瑾等，2007；馬勝利等，2008）結(jié)果吻合。

圖1 雁列構(gòu)造中的破壞區(qū)Fig.1 Failure zones in echelon faults.a擠壓雁列構(gòu)造；b拉張雁列構(gòu)造

1.2 事件的頻次-能量釋放關(guān)系的斜率絕對值（b0值）的計算方法

一個單元儲存的彈性應(yīng)變能依賴于3個主應(yīng)力、彈性模量、泊松比和單元的體積。一旦單元發(fā)生剪切破壞，則計算出的應(yīng)變能為儲存的剪切應(yīng)變能，反之，一旦單元發(fā)生拉伸破壞，則為儲存的拉伸應(yīng)變能。只有發(fā)現(xiàn)儲存的彈性應(yīng)變能在一定的時步間隔（時段）內(nèi)是下降的，才將其記錄下來，這一部分能量是釋放的彈性應(yīng)變能。對于釋放的彈性應(yīng)變能，可以利用多種方式進行統(tǒng)計，例如，統(tǒng)計其時空分布、總量的演變規(guī)律及頻次-能量釋放關(guān)系等。

下面，簡介事件的頻次-能量釋放關(guān)系的斜率絕對值（b0值）的計算方法（Wang et al．，2010，2012）：首先，將各單元釋放的彈性應(yīng)變能（一個單元釋放的彈性應(yīng)變能被看作是一個事件引起的）從小到大排列，最大值為Emax。由于在能量數(shù)據(jù)中，往往高值少而低值多。如將全部數(shù)據(jù)都保留，有時并不能獲得理想的b0值統(tǒng)計結(jié)果。為此，往往需要舍棄一部分?jǐn)?shù)據(jù)。設(shè)P為截斷指標(biāo)，將各單元的能量和PEmax相比較，舍棄那些超過PEmax的能量，僅對＜PEmax的能量進行統(tǒng)計。將剩余的數(shù)據(jù)又分為Q級，Q稱為分級指標(biāo)，確定每級的單元數(shù)目及其對數(shù)，這樣即可獲得頻次-能量釋放關(guān)系。最后，通過線性擬合即可獲得b0值。在任一時段之內(nèi)，均進行一次上述操作，即可獲得b0值的演變規(guī)律。

在下列2種條件下，無法計算出b0值：

（1）事件的能量分布規(guī)律反常，例如，小事件少，而大事件多；

（2）數(shù)據(jù)量不夠，例如，事件僅集中在1個能量級別之內(nèi)。

應(yīng)當(dāng)指出，在計算b0值時，未對能量取對數(shù)。這是由于，數(shù)值結(jié)果顯示，頻次-能量釋放關(guān)系基本上已呈線性，而非冪律形式，其原因可能如下：

（1）單元的尺寸偏大，使事件的數(shù)目總體上偏少，Q值偏小；

（2）未對Weibull分布作任何截斷；

（3）不同強度單元的脆性不同；

（4）僅統(tǒng)計每個能量級別之內(nèi)的事件數(shù)目，沒有統(tǒng)計比其大的事件。

通過2種雁列構(gòu)造的計算和分析，觀察到了b0值的3種表現(xiàn)：下降、先降后升及不變。不是每個引起能量釋放的事件均能引起b0值的下降，b0值的下降只發(fā)生在一定的變形階段。斷層系統(tǒng)的結(jié)構(gòu)變化（或出現(xiàn)新的破壞區(qū)）不一定能引起b0值的改變。能引起雁列構(gòu)造b0值下降的主要階段是斷層活化的階段以及擠壓雁列區(qū)的貫通階段（圖2a，c，d），而在拉張雁列區(qū)內(nèi)部的貫通階段和擠壓雁列區(qū)外側(cè)的破壞擴展階段中（圖2a，c），b0值未見明顯的變化。在拉張雁列構(gòu)造的應(yīng)力峰附近，b0值出現(xiàn)回升（圖2b）。

在擠壓雁列區(qū)的貫通階段，b0值發(fā)生了明顯的下降?？偟膩碚f，擠壓雁列區(qū)貫通階段破壞區(qū)的擴展、連接的增強，發(fā)生大量釋放高能量的事件，使b0值下降。在擠壓雁列區(qū)貫通階段b0值的下降應(yīng)該出現(xiàn)在擴展斷層附近區(qū)域，這與斷層活化階段的表現(xiàn)范圍有很大差別，分清2個階段的差別對于分析區(qū)域變形狀態(tài)有重要意義。擠壓雁列構(gòu)造的計算結(jié)果與破裂型地震或混合型地震之前的異常較多的事實類似，而拉張雁列構(gòu)造的計算結(jié)果與走滑型地震之前的異常較少的事實類似。

1.3 剪切應(yīng)變降——一個更為靈敏的前兆指標(biāo)

圖2 兩類雁列構(gòu)造在不同變形階段b0值的演變規(guī)律（Wang et al．，2012）Fig.2 Evolution of b0-value during different deformational stages for two kinds of echelon faults（Wang et al．，2012）.a拉張雁列構(gòu)造的斷層活化和巖橋破壞階段，虛線代表巖橋內(nèi)部開始出現(xiàn)破壞單元；b拉張雁列構(gòu)造的應(yīng)力峰之前及之后階段，虛線代表應(yīng)力峰；c擠壓雁列構(gòu)造的斷層活化和破壞區(qū)向外擴展階段，點線代表破壞區(qū)開始向外擴展；d擠壓雁列構(gòu)造的巖橋破壞階段，點線代表巖橋發(fā)生貫通，實線代表局部的應(yīng)力峰值

在數(shù)值模型中，可以方便地模擬出破壞單元數(shù)目的演變規(guī)律或能量釋放的時空分布規(guī)律，以比擬物理實驗中的聲發(fā)射或天然的地震活動。然而，只有破壞的單元才有聲發(fā)射。也就是說，統(tǒng)計聲發(fā)射時只統(tǒng)計那些破壞單元的信息，這是一種“一對一”的關(guān)系（當(dāng)一個單元發(fā)生破壞時，即輸入是1，或者聲發(fā)射累計數(shù)增加1，或者計及這個單元釋放的能量，因此，輸出也是1），而破壞單元的數(shù)量通常并不多，這對于前兆的研究不利，因為畢竟發(fā)生破壞單元的數(shù)目和總量相比，通常并不占優(yōu)勢。若能建立一種“一對多”的關(guān)系，即輸入是1時，輸出的是多個單元的信息，那么，前兆可能會更加明顯，易于識別。

FLAC-3D提供了計算單元剪切應(yīng)變增量的功能，剪切應(yīng)變增量和偏應(yīng)變張量的增量有關(guān)。通過對FLAC-3D的二次開發(fā)，在標(biāo)本全場內(nèi)可以對大于某一臨界值（如5×10-6）的剪切應(yīng)變降進行捕捉，并對其進行各種統(tǒng)計（Wang et al．，2013a），例如，計算剪切應(yīng)變降的累計（在任一時段之內(nèi)，將所有單元的剪切應(yīng)變降量求和）、最大值及發(fā)生剪切應(yīng)變降的單元數(shù)目等。數(shù)值模擬可以獲取剪切應(yīng)變降全場的信息，也方便進行各種統(tǒng)計，而實驗和野外應(yīng)變的點式觀測獲取的信息則相對少得多，更難以統(tǒng)計。

圖3 不同斷層間距a時擠壓雁列構(gòu)造破壞過程中2種量的演變規(guī)律（Wang et al．，2013a）Fig.3 Evolution of two quantities during the jog failure for compressive echelon faults with different fault spacing a（after Wang et al．，2013a）.

在擠壓雁列區(qū)貫通過程中，與剪切應(yīng)變降有關(guān)的統(tǒng)計量表現(xiàn)比較突出（圖3a），而與聲發(fā)射或能量釋放有關(guān)的量的表現(xiàn)則不然（圖3b），變化相對較小。從理論上看，這是由于一個發(fā)生破壞（應(yīng)變?nèi)趸┑膯卧?，只能記為一個聲發(fā)射，但是它卻能引起其周圍大量未破壞單元的反應(yīng)（卸載），形成強烈的前兆信號。對于不同的斷層系統(tǒng)，一些與剪切應(yīng)變降有關(guān)的量的演變規(guī)律是否存在共性和差異，值得進一步探討。

2 斷層粘滑的數(shù)值模擬

2.1 數(shù)值模型評述

目前，多種模型和方法已被用于粘滑過程模擬，例如：彈簧-滑塊模型、邊界元方法、有限元法、有限差分方法、離散元方法及格構(gòu)方法等。粘滑過程的模擬一般需要引入各種摩擦定律，尤其是速率和狀態(tài)依賴摩擦定律，但不盡然。在連續(xù)介質(zhì)模型中，引入空間或幾何上的非均勻性，或者考慮一些量（例如，孔隙壓力等）隨深度的變化規(guī)律，或者采用一些非連續(xù)介質(zhì)模型（例如，離散元和格構(gòu)方法等），均可產(chǎn)生時空復(fù)雜的地震行為。

彈簧-滑塊模型經(jīng)常被引用來檢驗新定律和研究新現(xiàn)象（Miyatake，1992；Scholz，1998；Ma et al．，2001；何昌榮，2003；Kato et al．，2003）?；瑝K受到重力、外力、彈簧力、摩擦力和平面或斜面的約束反力等力的作用。Miyatake（1992）采用1個二維的彈簧-滑塊（其周圍有4個彈簧）模型模擬了粘滑過程，并且考慮了摩擦系數(shù)的非均勻性。Mitsui等（2004）采用5個彈簧和滑塊模擬了一條狹長斷層的不同區(qū)段，并且將每段參數(shù)的取值盡量與實際相符。彈簧-滑塊模型雖形式上簡單，但所含物理意義豐富，特別適于地震機理、過程及各種因素的影響研究；然而，其考慮因素畢竟有限。

邊界元方法廣泛應(yīng)用于各種斷層的粘滑過程模擬研究（Lapusta et al．，2000；Bizzarri et al．，2001）。在計算中，可以采用多種摩擦定律，也可以通過規(guī)定不同的臨界滑動距離以引入斷層的非均質(zhì)性（Hillers et al．，2006）。邊界元方法只需要在邊界上建立方程，區(qū)域內(nèi)部采用積分方法，計算量小，適用于彈性動力學(xué)問題研究，但不允許介質(zhì)發(fā)生破壞。

有限元和有限差分法通過將計算區(qū)域離散成若干單元的方式進行數(shù)值求解。一般認(rèn)為，前者的適用性更廣，而后者的網(wǎng)格要求規(guī)則。其實不然，對于非規(guī)則網(wǎng)格，有限差分法亦能方便處理。在FLAC-3D中，僅運動方程的求解用到了有限差分方法，網(wǎng)格可以是6面體，不要求必須是立方體。上述2種方法可以處理復(fù)雜的斷層系統(tǒng)，例如，模型中可以包含不同尺度、方位及幾何特性的各種斷層，可以允許斷層之外單元發(fā)生破壞，可以引入各種摩擦定律，可以方便地考慮不同區(qū)域巖性的差異，可以逼真地刻畫斷層的結(jié)構(gòu)（例如，斷層面的凹凸不平、光滑的斷層面和含斷層泥的斷層面等），可以考慮力-熱-流體的耦合作用，可以考慮斷層之間的相互影響和作用，因此，被認(rèn)為是最有前景的數(shù)值方法。到目前為止，已有一些研究人員采用上述方法模擬了斷層的粘滑行為（Xing et al．，2003，2004；朱守彪等，2008，2009；薛霆虓等，2009；Ozaki，2011）。斷層通常有3種處理方式：1）將斷層視為材料屬性區(qū)別于其他部分的實體單元；2）將斷層視為各種類型的界面接觸單元；3）將斷層視為接觸邊界。Xing等（2004）指出，將速率和狀態(tài)依賴摩擦定律引入到有限元法中的結(jié)果未見報道，大多數(shù)人在模擬中都使用了斷層面固定的摩擦系數(shù)，不因條件的變化而改變，因而無法模擬地震現(xiàn)象。他們將“節(jié)點到點”的接觸單元策略引入到靜態(tài)顯式算法中研究了5.6°拐折斷層對斷裂形核、終止及重新開始的影響，雖然引入了一個簡單的速率和狀態(tài)依賴摩擦定律，但是，未見模擬出應(yīng)力的周期性變化。Ozaki（2011）將速率和狀態(tài)依賴摩擦定律引入到動態(tài)有限元法中，研究了2個塊體之間的接觸問題及不同參數(shù)時系統(tǒng)的動力粘滑行為，并且指出了庫侖摩擦定律的不足。

顆粒流、離散元、細胞自動機、格構(gòu)方法等純粹的非連續(xù)介質(zhì)方法也被用于粘滑過程模擬。Lorig等（1990）采用離散元模擬了彈簧-滑塊系統(tǒng)，引入了速率和狀態(tài)依賴摩擦定律，但是模擬出的所謂的粘滑失穩(wěn)現(xiàn)象（圖4a）與實驗結(jié)果差距較大。顆粒流方法不必引入速率和狀態(tài)依賴摩擦定律，一般僅需要引入庫侖摩擦定律就能模擬出類似粘滑的現(xiàn)象，例如，Mora等（1994）采用顆粒模型模擬了粘滑行為（圖4b），沒有引入任何摩擦特性；Abe等（2006）將斷層看作離散的顆粒的組合，顆粒之間的摩擦系數(shù)都是相同的，從他們模擬出的所謂的粘滑現(xiàn)象（圖4c）可以發(fā)現(xiàn)，模擬結(jié)果與實驗結(jié)果差別較大。將巖石簡化為數(shù)量龐大的圓形顆粒的顆粒流方法，不適于模擬大尺度問題，即使是使用并行計算技術(shù)也很難勝任。Doz等（2000）采用格構(gòu)模型，即將立方體單元離散成沿棱桿和對角桿，模擬了位于一個平面上的塊體的滑動過程，并沒有引入任何摩擦定律，也模擬出了他們聲稱的粘滑現(xiàn)象（圖4d），但是，與實驗結(jié)果差距較大。

2.2 兩種摩擦定律

在滑動弱化定律中，摩擦力僅依賴于滑動距離（Andrews，1985；Harris et al．，1993；He，1995；Ide et al．，1997；Ampuero et al．，2008；Tinti et al．，2009）。由于這一模型本質(zhì)上忽略了時間和速率效應(yīng)，因此它不能復(fù)現(xiàn)一個以上地震循環(huán)，除非當(dāng)滑動停止時，重新設(shè)置摩擦力的水平。Ohnaka等（1989）和Ohnaka（1993）提出了考慮初始滑動強化的滑動弱化模型。Bizzarri等（2001）討論了初始滑動強化的意義。

速率和狀態(tài)依賴摩擦定律亦有不同版本（Dieterich，1979，1992，1994；Ruina，1983；Kato et al．，2003；Kato et al．，2007），例如，慢度定律、滑動定律和復(fù)合定律。這一定律被認(rèn)為是斷層行為的更完全的本構(gòu)描述。然而，目前還沒有哪一個定律被公認(rèn)能描述所有的實驗數(shù)據(jù)，似乎也不清楚哪一個更適合描述地震循環(huán)（Kato et al．，2003）。Okubo（1989）和Bizzarri等（2001）詳細比較了滑動弱化定律與速率和狀態(tài)依賴摩擦定律的結(jié)果差別。Bizzarri等（2001）強調(diào)，速率和狀態(tài)依賴摩擦定律包括了滑動弱化模型的最突出特點。Okubo（1989）指出，在速率和狀態(tài)依賴定律中進行高速度截斷，則斷層的行為就像滑動弱化定律的一樣。

圖4 純粹的非連續(xù)介質(zhì)方法的結(jié)果Fig.4 Results of pure discontinuousmedium methods.a Lorig et al.，1990；b Mora et al.，1994；c Abe et al．，2006；d Doz et al.，2000

2.3 一種摩擦強化-摩擦弱化模型

FLAC-3D在每個節(jié)點上建立運動方程，采用中心差分方法求解。節(jié)點在節(jié)點力的作用下，將發(fā)生平動。節(jié)點力包括單元應(yīng)力引起的彈性力、阻尼力、重力和外力等。在每個時步之內(nèi)，幾何方程、本構(gòu)方程和運動方程會被依次執(zhí)行1遍，直到獲得滿意結(jié)果為止。FLAC-3D中的應(yīng)力、應(yīng)變等量存儲在單元上，而速度、位移、節(jié)點質(zhì)量等存儲在節(jié)點上。單元和節(jié)點之間的信息通過幾何方程和虛功原理溝通。

FLAC-3D既適于動力分析，也適于擬靜力分析。后者尤其適于模擬實驗室中的準(zhǔn)靜態(tài)加載實驗和長期的構(gòu)造運動。在模擬靜力分析時，時步長度取為1。相應(yīng)地，節(jié)點質(zhì)量必須放大，以滿足數(shù)值上的穩(wěn)定性條件。對于動力分析，時步長度取值較小，這是為了獲得穩(wěn)定的數(shù)值解答。動力分析適于模擬瞬態(tài)行為，并不適于模擬實驗室中的準(zhǔn)靜態(tài)加載實驗，因為這將需要極長的計算時間。

采用摩擦強化-摩擦弱化模型（圖5），Wang等（2013b）模擬了典型斷層的粘滑現(xiàn)象，該模型在FLAC-3D中執(zhí)行，可以復(fù)現(xiàn)起源于斷層面上的應(yīng)力的周期性變化，后傳遞至標(biāo)本端部。作為一個關(guān)鍵，在一個粘滑循環(huán)之內(nèi)，我們引入了一個增量的塑性剪切應(yīng)變。它在弱化過程中從零開始增加，然而一旦達到最大值后被清零。上述變化是為了控制斷層單元內(nèi)摩擦角φ的變化。在粘著和滑動2個階段，僅要求φ的演變規(guī)律不同，運動方程求解方法及參數(shù)完全相同。

為了正確理解上述摩擦強化-摩擦弱化模型，需要區(qū)分2種塑性剪切應(yīng)變增量。一種是專門用于控制不同階段內(nèi)摩擦角的不同演變規(guī)律的，記為Δγp，其最大值記為Δγpm。Δγp在FLAC-3D軟件中原本并不存在，需要對其二次開發(fā)后引入。另一種是FLAC-3D中原有的，記為Δεp，它用于描述單元的塑性變形程度，并用于決定當(dāng)前的弱化參數(shù)。

FLAC-3D每次調(diào)用本構(gòu)方程，都會根據(jù)流動法則計算出一個Δεp。顯然，斷層單元當(dāng)前總的塑性剪切應(yīng)變εp應(yīng)為此前所有Δεp之和。

Δεp和Δγp關(guān)系密切。在弱化階段，當(dāng)φ從最大值φm突然降至最小值φ0的過程中，將導(dǎo)致Δγp的不斷增加。這一階段將經(jīng)歷若干次計算循環(huán)。在每1次計算循環(huán)中，由于存在Δφ，都將產(chǎn)生一個Δεp。在本次粘滑循環(huán)內(nèi)，將從弱化階段開始的各個計算循環(huán)內(nèi)的Δεp求和即可得到Δγp。在隨后的硬化階段，由于位移控制加載的繼續(xù)進行，Δγp將被不斷提升，從而使φ從φ0逐漸增至φm。我們假定，φ與Δγp呈線性關(guān)系。

總之，在弱化階段，φ的降低導(dǎo)致了Δγp的增加；在硬化階段，Δγp的繼續(xù)增加，導(dǎo)致了φ的增加。上述計算循環(huán)不斷地進行下去即可模擬粘滑現(xiàn)象。φ的下降將導(dǎo)致節(jié)點速度的增加，而φ的上升將導(dǎo)致節(jié)點速度的降低。所以，該模型能反映出節(jié)點速度的變化，這是由于對節(jié)點運動方程進行求解的緣故。也就是說，速度被作為一種響應(yīng)處理，而非作為輸入量。

該模型經(jīng)由繆阿麗等（2010）開展的一個雙軸壓縮粘滑實驗進行了檢驗。含60°單一斷層標(biāo)本的差應(yīng)力-時間曲線見圖6a?？梢园l(fā)現(xiàn)，加載速度v=0.1μm／s時的粘滑周期大約是v= 0.05μm／s時的一半；v=0.1μm／s時的6次粘滑事件共經(jīng)歷了290μm的縱向位移增量。在實驗中，為了避免壓頭之間的干涉，標(biāo)本的4個角被切去。在數(shù)值模擬中，完全真實地考慮標(biāo)本的形狀并無必要。所以，僅對標(biāo)本中央的一部分進行建模，去掉兩側(cè)，選擇四邊形單元。圍壓、標(biāo)本高度、斷層寬度等參數(shù)取值都盡量與實驗條件一致。加載方式選擇準(zhǔn)靜態(tài)位移控制加載，在平面應(yīng)變條件下進行計算。

由圖6a，b可以發(fā)現(xiàn)，2種速度條件下的數(shù)值和實驗結(jié)果吻合較好；v0=1×10-9m／時步時的粘滑周期是v0=2×10-9m／時步時的2倍。

為了保證數(shù)值和實驗結(jié)果的一致，必須要保證數(shù)值和實驗中的一次粘滑事件所對應(yīng)的標(biāo)本縱向位移增量相同。在準(zhǔn)靜態(tài)加載時，加載速度應(yīng)足夠小，為了模擬v=0.1μm／s時的結(jié)果，加載速度取為v0=2×10-9m／時步。這樣，上述6次粘滑事件共需要7.3×104個時步，對應(yīng)2,900s。因此，1個時步對應(yīng)0.04s，或者1s對應(yīng)25個時步。

特別需要指出，在數(shù)值計算中，在標(biāo)本上、下兩端，均進行位移控制加載。和僅在一端進行位移控制加載的情形相比，這將節(jié)省一半的計算時間。為了模擬上述6次粘滑事件，在普通配置的微機上僅需計算8min。但是，若采用動力加載方式，即加載速度取為實際值，則至少需要1,800d，這是一個不可接受的時間。所以，在計算中，在不影響計算精度和效率的前提下，采用準(zhǔn)靜態(tài)加載方式模擬某些物理實驗是一個經(jīng)濟、可行的選擇。

圖5 摩擦強化-摩擦弱化模型Fig.5 A frictional-hardening and frictional-softeningmodel.

圖6 雙軸壓縮條件下含單一斷層標(biāo)本的差應(yīng)力的演變規(guī)律Fig.6 Evolution of the principal stress difference for samples with a single straight fault under biaxial compression.a繆阿麗等，2010；b Wang et al．，2013b

圖7a給出了雙軸壓縮條件下4個45°單一斷層標(biāo)本的應(yīng)力-時步數(shù)目曲線的數(shù)值模擬結(jié)果，其中一個標(biāo)本被劃分成若干正方形，而另外3個標(biāo)本均被劃分成若干四邊形。由此可以發(fā)現(xiàn)：

（1）和平坦的斷層面相比，不平坦的斷層面使滑動開始和結(jié)束時的應(yīng)力水平有所提高，但并不改變粘滑的周期；

（2）斷層厚度w越小，粘滑周期越小，應(yīng)力降越小。這是由于，在不變的前提下，斷層厚度越小，則意味著臨界滑動距離越小，所以粘滑周期越小。

圖7b給出了一個針對45°單一斷層標(biāo)本的動力計算結(jié)果，其中，還給出了位于斷層上的1個節(jié)點的速度演變規(guī)律。由此可以發(fā)現(xiàn)，在斷層快速滑動過程中，速度達到最大值；在摩擦硬化過程中，速度有所下降，在此過程中，可以發(fā)現(xiàn)若干小事件。

圖7 應(yīng)力及節(jié)點速度的演變規(guī)律Fig.7 Evolution of stress and nodal velocity with the number of timesteps.

Wang等（2013b）共開展了40個算例的計算，阻尼采用局部自適應(yīng)阻尼，由于它不需明確頻率，所以比瑞利阻尼簡單。局部自適應(yīng)阻尼力的大小和節(jié)點的不平衡力有關(guān)，而方向和速度方向有關(guān)，這與粘性阻尼力不同。標(biāo)本中包含單一、拐折、交叉和雁列等典型斷層構(gòu)造。單元形狀包括正方形和四邊形。對于前者，斷層面不可避免地呈現(xiàn)鋸齒形，這會增加斷層錯動的阻力；而對于后者，斷層面保持平整。加載方式包括準(zhǔn)靜態(tài)和動力加載。變形模式包括小變形和大變形。重點考察了加載條件（圍壓、加載速度）、斷層寬度、斷層幾何對粘滑行為的影響。

3 展望

斷層系統(tǒng)破裂和粘滑過程的綜合模擬是極富挑戰(zhàn)性的力學(xué)、地學(xué)難題，欲成功解決尚有大量工作需要進行：

（1）亟待發(fā)展適于地質(zhì)體材料變形、破壞及運動整個過程模擬的計算力學(xué)模型。目前的計算模型對于破壞和粘滑問題往往單獨處理。對于破壞問題，力學(xué)界較為關(guān)注；對于粘滑問題，地震界尤為重視。針對每一問題都發(fā)展出了不少較為有效的計算模型。一般地，地震粘滑過程本質(zhì)上是非均勻的（李普春等，2013）。在提出的摩擦強化-摩擦弱化模型中引入單元力學(xué)性質(zhì)或臨界滑動距離的非均勻性并不困難，斷層之外巖石的破壞亦方便考慮。這樣，一系列復(fù)雜的地震或礦震現(xiàn)象（例如，前兆、斷裂擴展、斷層運動等）都可望能在一定程度上進行模擬，因而具有較為廣闊的應(yīng)用前景。

（2）亟待發(fā)展復(fù)雜斷層系統(tǒng)中斷層穩(wěn)定性的可靠判別方法。斷層的滑動有穩(wěn)定滑動和失穩(wěn)滑動之區(qū)別。根據(jù)斷層位置應(yīng)力或應(yīng)變足夠高，或者滑動速率足夠大等經(jīng)驗方法難以準(zhǔn)確地判別斷層的穩(wěn)定性。斷層不一定平直，斷層之間存在相互影響和作用。一般可根據(jù)其應(yīng)力狀態(tài)或幾何特性等劃分成若干區(qū)段。不同區(qū)段不一定同時失穩(wěn)。斷層的一次大失穩(wěn)事件，必然是若干區(qū)段同時失穩(wěn)的結(jié)果。最近的實驗觀測（云龍等，2011；郭彥雙等，2012；馬瑾，2012）發(fā)現(xiàn)，斷層失穩(wěn)即將啟動之前，斷層呈現(xiàn)協(xié)同變形特征。有必要通過數(shù)學(xué)、力學(xué)手段建立在數(shù)值模型中方便使用的可靠的斷層穩(wěn)定性判別方法，其正確性需經(jīng)物理實驗的檢驗。

（3）巖石動態(tài)斷裂和摩擦失穩(wěn)等動力學(xué)研究應(yīng)是解決地震預(yù)測難題最根本的出路（劉啟元，2005）。亟待發(fā)展動力學(xué)方程的高效求解方法。在本文所述的準(zhǔn)靜態(tài)計算中，沒有時間的概念，時步長度取為1，這給構(gòu)造運動及地震孕育等長期過程的模擬帶來方便，但不利于斷層滑動失穩(wěn)等瞬間發(fā)生現(xiàn)象的模擬。在動態(tài)計算中，為了確保數(shù)值穩(wěn)定性，時步長度往往保持恒定，且極小。上文提過，采用動態(tài)計算模型模擬一個常見的雙軸壓縮實驗是不可行的。Lapusta等（2000）發(fā)展了依賴于最大滑動速度的變時步中心差分方法，值得借鑒和探索。自始至終都采用顯式計算是否有必要？顯式計算和隱式計算能否實現(xiàn)高效率、高精度的耦合是值得深入思考和探索的問題之一。

本文著重對實驗室尺度典型斷層系統(tǒng)破壞、前兆及粘滑過程的數(shù)值模擬進行了簡要的回顧和討論，并介紹了第1作者近年來在博士后合作導(dǎo)師指導(dǎo)下的積極探索。相信，這種探索和介紹有助于引起學(xué)者們對斷層系統(tǒng)破裂和粘滑過程的綜合模擬問題的關(guān)注，并有助于上述問題的最終解決。

致謝 衷心感謝本文第1作者的博士后合作導(dǎo)師馬瑾院士的指導(dǎo)，感謝審稿專家提出的寶貴修改意見與建議。

鄧志輝，宋鍵，孫君秀，等.2011a.數(shù)值模擬方法在地震預(yù)測研究中應(yīng)用的初步探討（Ⅰ）［J］.地震地質(zhì)，33（3）：660—669.doi：10.3969／j.issn.0253-4967.2011.03.015.

DENG Zhi-hui，SONG Jian，SUN Jun-xiu，et al．2011a.Preliminary study on application of numerical simulation methods to earthquake prediction research（Ⅰ）［J］.Seismology and Geology，33（3）：660—669（in Chinese）.

鄧志輝，胡勐乾，周斌，等.2011b.數(shù)值模擬方法在地震預(yù)測研究中應(yīng)用的初步探討（Ⅱ）［J］.地震地質(zhì)，33（3）：670—683.doi：10.3969／j.issn.0253-4967.2011.03.016.

DENG Zhi-hui，HU Meng-qian，ZHOU Bin，et al．2011b.Preliminary study on application of numerical simulation methods to earthquake prediction research（Ⅱ）［J］.Seismology and Geology，33（3）：670—683（in Chinese）.

郭婷婷，徐錫偉.2011.有限元法在構(gòu)造應(yīng)力場與地震預(yù)測研究中的應(yīng)用與發(fā)展［J］.地震研究，34（2）：246—253.

GUO Ting-ting，XU Xi-wei.2011.App lication and progress of the finite-elementmethod in structural stress field and earthquake prediction［J］.Journal of Seismological Research，34（2）：246—253（in Chinese）.

郭彥雙，馬瑾，卓燕群.2012.平直斷層失穩(wěn)前應(yīng)變場演化特征的實驗研究［J］.國際地震動態(tài)，（6）：244.

GUO Yan-shuang，MA Jin，ZHUO Yan-qun.2012.Experimental study of evolution of strain field prior to instability for a straight fault［J］.Recent Developments in World Seismology，（6）：244（in Chinese）.

何昌榮.2003.速度與狀態(tài)依賴性摩擦本構(gòu)關(guān)系下雙滑塊系統(tǒng)的相互作用［J］.中國科學(xué)（D輯），33（增）：53—59.

HE Chang-rong.2003.Interaction between two sliders in a system with rate-and state-dependent friction［J］.Science in China（Ser D），33（Suppl）：53—59（in Chinese）.

蔣海昆，馬勝利，張流，等.2002.雁列式斷層組合變形過程中的聲發(fā)射活動特征［J］.地震學(xué)報，24（4）：385—396.

JIANG Hai-kun，MA Sheng-li，ZHANG Liu，et al．2002.Spatio-temporal characteristics of acoustic emission during the deformation of rock samples with compressional and extensional en-echelon fault pattern［J］.Acta Seismologica Sinica，24（4）：385—396（in Chinese）.

李玉江，陳連旺，葉際陽.2009.數(shù)值模擬方法在應(yīng)力場演化及地震科學(xué)中的研究進展［J］.地球物理學(xué)進展，24（2）：418—431.

LIYu-jiang，CHEN Lian-wang，YE Ji-yang.2009.Application and development of numerical simulation in stress field development and seismology［J］.Progress in Geophys，24（2）：418—431（in Chinese）.

劉力強，馬勝利，馬瑾，等.1999.巖石構(gòu)造對聲發(fā)射統(tǒng)計特征的影響［J］.地震地質(zhì)，21（4）：377—386.

LIU Li-qiang，MA Sheng-li，MA Jin，et al．1999.Effect of rock structure on statistic characteristics of acoustic emission［J］.Seismology and Geology，21（4）：377—386（in Chinese）.

李普春，劉力強，郭玲莉，等.2013.粘滑過程中的多點錯動［J］.地震地質(zhì)，35（1）：125—137.doi：10.3969／j. issn.0253-4967.2013.01.011.

LIPu-chun，LIU Li-qiang，GUO Ling-li，et al．2013.Multi-point dislocation in stick-slip process［J］.Seismology and Geology，35（1）：125—137（in Chinese）.

劉啟元.2005.地震預(yù)報研究的主攻方向：動力數(shù)值預(yù)測［J］.國際地震動態(tài)，5：63—68.

LIU Qi-yuan.2005.Main task of earthquake prediction research：Numerical dynam ic forecast［J］.Recent Developments in World Seismology，5：63—68（in Chinese）.

馬瑾.2012.地震新主體地區(qū)和失穩(wěn)危險階段的探尋［J］.國際地震動態(tài)，（6）：231.

MA Jin.2012.Study of the new seism ic major region and hazardously unstable stages［J］.Recent Development in World Seismology，（6）：231（in Chinese）.

馬瑾，劉力強，劉培洵，等.2007.斷層失穩(wěn)錯動熱場前兆模式：雁列斷層的實驗研究［J］.地球物理學(xué)報，50（4）：1141—1149.

MA Jin，LIU Li-qiang，LIU Pei-xun，et al．2007.Thermal precursory pattern of fault unstable sliding：An experimental study of en echelon faults［J］.Chinese JGeophys，50（4）：1141—1149（in Chinese）.

馬勝利，陳順云，劉培洵，等.2008.斷層階區(qū)對滑動行為影響的實驗研究［J］.中國科學(xué)（D輯），38（7）：842—851.

MA Sheng-li，CHEN Shun-yun，LIU Pei-xun，et al．2008.Effect of fault jogs on frictional behavior：An experimental study［J］.Science in China（Ser D），38（7）：842—851（in Chinese）.

馬勝利，雷興林，劉力強.2004.標(biāo)本非均勻性對巖石變形聲發(fā)射時空分布的影響及其地震學(xué)意義［J］.地球物理學(xué)報，47（1）：127—131.

MA Sheng-li，LEI Xing-lin，LIU Li-qiang.2004.Effects of heterogeneity in rock samples on spatial and temporal distribution of acoustic emission and their significance in seismology［J］.Chinese JGeophys，47（1）：127—131（in Chinese）.

馬勝利，馬瑾.2003.中國實驗巖石力學(xué)與構(gòu)造物理學(xué)研究的若干新進展［J］.地震學(xué)報，25（5）：453—464.

MA Sheng-li，MA Jin.2003.Recent progress in studies of experimental rock mechanics and tectonophysics in China［J］.Acta Seismologica Sinica，25（5）：453—464（in Chinese）.

馬勝利，馬瑾，劉力強.2002.地震成核相的實驗證據(jù)［J］.科學(xué)通報，47（5）：387—391.

MA Sheng-li，MA Jin，LIU Li-qiang.2002.Experimental evidence for seismic nucleation phase［J］.Chinese Science Bulletin，47（5）：387—391（in Chinese）.

繆阿麗，馬勝利，周永勝.2010.硬石膏斷層帶摩擦穩(wěn)定性轉(zhuǎn)換與微破裂特征的實驗研究［J］.地球物理學(xué)報，53（11）：2671—2680.

MIAO A-li，MA Sheng-li，ZHOU Yong-sheng.2010.Experimental study on frictional stability transition and microfracturing characteristics for anhydrite fault zones［J］.Chinese JGeophys，53（11）：2671—2680（in Chinese）.

時振梁，王健，張曉東.1995.中國地震活動性分區(qū)特征［J］.地震學(xué)報，17：20—24.

SHIZhen-liang，WANG Jian，ZHANG Xiao-dong.1995.The regional characteristics of seismic activity in China［J］. Acta Seismologica Sinica，17：20—24（in Chinese）.

唐春安，傅宇方，趙文.1997.震源孕育模式的數(shù)值模擬研究［J］.地震學(xué)報，19（4）：337—346.

TANG Chun-an，F(xiàn)U Yu-fang，ZHAO Wen.1997.A new approach to numerical simulation of source development of earthquake［J］.Acta Seismologica Sinica，19（4）：337—346（in Chinese）.

陶瑋，洪漢凈，劉培洵.2000.中國大陸及鄰區(qū)強震活動主體地區(qū)形成的數(shù)值模擬［J］.地震學(xué)報，22（3）：271—277.

TAOWei，HONG Hai-jing，LIU Pei-xun.2000.The finite elementmodeling of strong earthquake activities in major regions in Chinese mainland and its vicinity［J］.Acta Seismologica Sinica，22（3）：271—277（in Chinese）.

滕吉文.2010.強烈地震孕育與發(fā)生的地點、時間及強度預(yù)測的思考與探討［J］.地球物理學(xué)報，53（8）：1749—1766.

TENG Ji-wen.2010.Ponder and research on the genesis and occurrence of strong earthquakes and the prediction of their place，time and intensity［J］.Chinese JGeophys，53（8）：1749—1766（in Chinese）.

王學(xué)濱.2007.擴容角對初始隨機材料缺陷巖石試樣破壞過程及力學(xué)行為的影響［J］.中國有色金屬學(xué)報，17（7）：1065—1071.

WANG Xue-bin.2007.Effects of dilation angle on failure process and mechanical behavior for rock specimen withrandom material imperfections［J］.The Chinese Journal of Nonferrous Metals，17（7）：1065—1071（in Chinese）.

薛霆虓，傅容珊，林峰.2009.幾何彎曲斷層活動性的模擬［J］.地球物理學(xué)報，52（10）：2509—2518.

XUE Ting-xiao，F(xiàn)U Rong-san，LIN Feng.2009.The simulation of activities of bend fault［J］.Chinese JGeophys，52（10）：2509—2518（in Chinese）.

云龍，郭彥雙，馬瑾.2011.5°拐折斷層在粘滑過程中物理場演化與交替活動的實驗研究［J］.地震地質(zhì)，33（2）：356—368.doi：10.3969／j.issn.0253-4967.2011.02.009.

YUN Long，GUO Yan-shuang，MA Jin.2011.An experimental study of evolution of physical field and the alternative activities during stick-slip of 5°bend fault［J］.Acta Seismologica Sinica，33（2）：356—368（in Chinese）.

鄭洪偉，李廷棟，高銳，等.2006.數(shù)值模擬在地球動力學(xué)中的研究進展［J］.地球物理學(xué)進展，21（2）：360—369.

ZHENG Hong-wei，LITing-dong，GAO Rui，et al．2006.The advance of numerical simulation in geodynamics［J］. Progress in Geophys，21（2）：360—369（in Chinese）.

朱守彪，邢會林，謝富仁，等.2008.地震發(fā)生過程的有限元法模擬：以蘇門答臘俯沖帶上的大地震為例［J］.地球物理學(xué)報，51（2）：460—468.

ZHU Shou-biao，XING Hui-lin，XIE Fu-ren，et al．2008 Simulation of earthquake processes by finite elementmethod：The case ofmegathrust earthquakes on the Sumatra subduction zone［J］.Chinese JGeophys，51（2）：460—468（in Chinese）.

朱守彪，張培震.2009.2008年汶川MS8.0地震發(fā)生過程的動力學(xué)機制研究［J］.地球物理學(xué)報，52（2）：418—427.

ZHU Shou-biao，ZHANG Pei-zhen.2009.A study on the dynamicalmechanisms of theWenchuan MS8.0 earthquake，2008［J］.Chinese JGeophys，52（2）：418—427（in Chinese）.

Abe S，Latham S，Mora P.2006.Dynamic rupture in a 3-D particle-based simulation of a rough planar fault［J］.Pure and Applied Geophysics，163：1881—1892.

Ampuero J P，Ben-Zion Y.2008.Cracks，pulses and macroscopic asymmetry of dynamic rupture on a bimaterial interface with velocity-weakening friction［J］.Geophys J Int，173：674—692.

Andrews D J.1985.Dynamic plane-strain shear rupture with a slip-weakening friction law calculated by a boundary integralmethod［J］.Bull Seismol Soc Am，75：1—21.

Barbot S，Lapusta N，Avouac JP.2012.Under the hood of the earthquakemachine：Toward predictivemodeling of the seismic cycle［J］.Science，336（6082）：707—710.

Beroza G C，M ikumo T.1996.Short slip duration in dynamic rupture in the presence of heterogeneous fault properties［J］.JGeophys Res，101（B10）：22449—22460.

Bizzarri A，Cocco M，Andrews D J，et al．2001.Solving the dynamic rupture problem with different numerical approaches and constitutive laws［J］.Geophys J Int，144：656—678.

Chen S，Yue Z Q，Tham L G.2004.Digital image-based numericalmodeling method for prediction of inhomogeneous rock failure［J］.Int JRock Mech Min Sci，41（6）：939—957.

Cundall P A，Strack O D L.1979.A discrete numericalmodel for granular assemblies［J］.Geotechnique，29（1）：47—65.

Dieterich JH.1979.Modeling of rock friction，1：Experimental results and constitutive equations［J］.JGeophys Res，84：2161—2168.

Dieterich J H.1992.Earthquake nucleation on faults with rate-and state-dependent strength［J］.Tectonophysics，211—134.

Dieterich JH.1994.A constitutive law for rate of earthquake production and its application to earthquake clustering［J］.JGeophys Res，99（B2）：2601—1618.

Doz G N，Riera J D.2000.Towards the numerical simulation of seismic excitation［J］.Nuclear Engineering and Design，196：253—261.

Erickson SG，Strayer L M，Suppe J.2001.Initiation and reactivation of faults during movement over a thrust-fault ramp：Numericalmechanicalmodels［J］.JStruct Geol，23：11—23.

Fang Z，Harrison J P.2002.Development of a local degradation approach to the modeling of brittle fracture in heterogeneous rocks［J］.Int JRock Mech M in Sci，39（4）：443—457.

Harris R A，Day S M.1993.Dynamic of fault interaction：Parallel strike-slip faults［J］.J Geophys Res，98 B3：4461—4472.

He C.1995.Slip-weakening constitutive relation and the structure in the vicinity of a shear crack tip［J］.Pure Appl Geophys，145：147—157.

Hillers G，Ben-Zion Y，Mai P M.2006.Seism icity on a fault controlled by rate-and state-dependent friction with spatial variations of the critical slip distance［J］.JGeophys Res，111（B01403）：1—23.

Hobbs B E，Ord A.1989.Numerical simulation of shear band formation in a frictional-dilational material［J］. Ingenigeur-Archiv，59（2）：209—220.

Hu N，Molinari JF.2004.Shear bands in densemetallic granularmaterials［J］.Journal of the Mechanics and Physics of Solids，52（3）：499—531.

Ide S，Takeo M.1997.Determination of constitutive relations of fault slip based on seismic wave analysis［J］.J Geophys Res，102：27379—27391.

Kato N，Tullis T E.2003.Numerical simulation of seismic cycles with a composite rate-and state-dependent friction law［J］.Bulletin of the Seismological Society of America，93（2）：841—853.

Kato N，Lei X，Wen X.2007.A synthetic seismicity model for the Xianshuhe Fault，Southwestern China：Simulation using a rate-and state-dependent friction law［J］.Geophys J Int，169：286—300.

Lapusta N，Rice J R，Ben-Zion Y，et al．2000.Elastodynamic analysis for slow tectonic loading with spontaneous rupture episodes on faults with rate-and state-dependent friction［J］.JGeophys Res，105 B10：23765—23789. Lei X L，Masuda K，Nishizawa O，et al．2004.Three typical stages of acoustic emission activity during the catastrophic fracture of heterogeneous faults in jointed rocks［J］.J Struc Geol，26：247—258.

Lorig L J，Hobbs B E.1990.Numericalmodeling of slip instability using the distinct elementmethod with state variable friction laws［J］.Int JRock Mech Min Sci，27：525—534.

Liu H Y，Kou SQ，Lindqvist P A，et al．2004.Numerical studies on the failure process and associated microseism icity in rock under triaxial compression［J］.Tectonophysics，384：149—174.

Ma S，He C.2001.Period doubling as a result of slip complexities in sliding surfaces with strength heterogeneity［J］. Tectonophysics，337：135—145.

Mckinnon SD，De La Barra IG.1998.Fracture initiation，growth and effect on stress field：A numerical investigation［J］.JStruct Geol，20（12）：1673—1689.

Mitsui N，Hirahara K.2004.Simple spring-mass model simulation of earthquake cycle along the Nankai trough in Southwest Japan［J］.Pure and Applied Geophysics，161：2433—2450.

Miyatake T.1992.Numerical simulation of the three-dimensional faulting processes with heterogeneous rate-and statedependent friction［J］.Tectonophys，211：223—232.

Mora P，Place D.1994.Simulation of the frictional stick-slip instability［J］.Pure Appl Geophys，143：61—87.

Noda H，Lapusta N.2013.Stable creeping fault segments can become destructive as a result of dynam ic weakening［J］.Nature，493：518—521.

Ohnaka M.1993.Critical size of the nucleation zone of earthquake rupture inferred from immediate foreshock activity［J］.J Phys Earth，41：45—56.

Ohnaka M，Yamashita T.1989.A cohesive zone model for dynamic shear faulting based on experimentally inferred constitutive relation and strong motion source parameters［J］.JGeophys Res，94：4089—4104.

Okubo P G.1989.Dynamic rupture modeling with laboratory-derived constitutive relations［J］.J Geophys Res，94：12321—12335.

Ozaki S.2011.Finite element analysis of rate-and state-dependent frictional behavior［J］.Key Eng Mat，462-463：547—552.

Ruina A L.1983.Slip instabilities and state variable friction laws［J］.JGeophys Res，88：10359—10370.

Scholz C H.1998.Earthquakes and friction laws［J］.Nature，391（1）：37—42.

Shibazaki B，Obara K，Matsuzawa T，et al．2012.Modeling of slow slip events along the deep subduction zone in the Kii Peninsula and Tokai regions，southwest Japan［J］.JGeophys Res，117（B6）.doi：10.1029／2011JB009083. Sibson R H.1985.Stopping of earthquake ruptures at dilatational jogs［J］.Nature，316：248—251.

Strayer L M，Hudleston P J.1997.Numericalmodeling of fold initiation at thrust ramps［J］.JStruct Geol，19：551—566.

Tinti E，Cocco M，F(xiàn)ukuyama E，et al．2009.Dependence of slip weakening distance（Dc）on final slip during dynam ic rupture of earthquakes［J］.Geophys J Int，177：1205—1220.

Wang X B，Ma J，Liu L Q.2010.Numerical simulation of failed zone propagation process and anomalies related to the released energy during a compressive jog intersection［J］.Journal of Mechanics of Materials and Structures，5（6）：1007—1022.

Wang X B，Ma J，Liu L Q.2012.A comparison ofmechanical behavior and frequency-energy relations for two kinds of echelon fault structures through numerical simulation［J］.Pure Appl Geophys，169：1927—1945.

Wang X B，Ma J，Liu L Q.2013a.Numerical simulation of large shear strain drops during jog failure for echelon faults based on a heterogeneous and strain-softening model［J］.Tectonophysics，608：667—684.

Wang X B，Ma J，Pan Y S.2013b.Numerical simulation of stick-slip behaviours of typical faults in biaxial compression based on a frictional-hardening and frictional-softening model［J］.Geophys J Int，194（2）：1023—1041.

Xing H L，Makinouchi A.2003.Finite element modeling of frictional instability between deformable rocks［J］. International Journal for Numerical and Analytical Methods in Geomechanics，27（12）：1005—1025.

Xing H L，Mora P，Makinouchi A.2004.Finite element analysis of fault bend influence on stick-slip instability along an intra-plate fault［J］.Pure and Applied Geophysics，161：2091—2102.

NUM ERICAL SIMULATION OF FAILURES，PRECURSORS AND STICK-SLIP PROCESSES FOR TYPICAL FAULT STRUCTURES AT A LABORATORY SCALE

WANG Xue-bin1，2，3）MA Bing2）LüJia-qing2）

1）State Key Laboratory of Earthquake Dynamics，Institute of Geology，China Earthquake Administration，Beijing 100029，China
2）College of Mechanics and Engineering，Liaoning Technical University，F(xiàn)uxin 123000，China
3）Institute of Computational Mechanics，Liaoning Technical University，F(xiàn)uxin 123000，China

For typical faults at a laboratory scale，with basically known physical and mechanical parameters and specific loading and boundary conditions，the integrative numerical simulation of their deformations，failures and stabilities is difficult.The numerical simulation at a laboratory scale can lay a foundation for that at a large scale，and can provide some necessary technical conditions.In this paper，from the viewpoint of numerical modeling，we reviewed the methods of introducing heterogeneities and reproducing stick-slip processes at a laboratory scale.Two models are specially emphasized：the heterogeneous strain-softening model with faults and rock blocks，and the frictional hardening and frictional softening model.In the firstmodel，to assess the state of stress for faults，we present the method to calculate the absolute value of the slope of the relation between frequency of events and released energy.During the jog failure for the compressive echelon faults，some quantities related to large shear strain drops are particularly sensitive indicators for rock fractures，and can be more effective than acoustic em issions.In the second model，the internal friction angle follows different laws at the stick and slip stages.Though no velocity is introduced in the model，the velocity changes at the stick and slip stages because of solving themotion equations.Three problems worthy of further investigation are pointed out.Firstly，a computationalmodel is needed to reproduce the entire process of geomaterials：deformation，failure and motion.Secondly，a reliable method to assess the stability of faults for a complex system is necessary.Finally，high efficient solving algorithm for the motion equations is especially needed.

laboratory-scale model，failure，stick-slip，precursor，friction，shear strain，echelon fault

P313

A

0253-4967（2014）03-0845-17

王學(xué)濱，男，1975年生，2006年在遼寧工程技術(shù)大學(xué)獲得工程力學(xué)專業(yè)博士學(xué)位，2009—2012年在中國地震局地質(zhì)研究所做博士后研究，教授，主要從事巖石力學(xué)、計算力學(xué)等方面的研究，電話：0418-3350924，E-mail：wxbbb＠263.net。

10.3969／j.issn.0253-4967.2014.03.023

2014-03-26收稿，2014-09-10改回。

地震動力學(xué)國家重點實驗室開放課題（LED2014B03）、高等學(xué)校博士學(xué)科點專項科研基金（20132121110006）和國家自然科學(xué)基金（51374122）共同資助。