基于速度－狀態(tài)摩擦本構(gòu)定律的三維DDA方法

張龍，江在森，武艷強(qiáng)，2＊，鄒鎮(zhèn)宇，3，劉曉霞，3，魏文薪

1中國地震局地震預(yù)測研究所地震預(yù)測重點實驗室，北京 100036

2中國地震局第一監(jiān)測中心，天津 300180

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

1 引言

非連續(xù)變形分析（DDA）方法是巖土工程領(lǐng)域模擬節(jié)理巖體運動與變形的一種數(shù)值解析方法（Shi，1992）.該方法將不連續(xù)面分割的巖體視為塊體集合，以最小勢能法建立總體平衡方程，利用單純形積分與開閉迭代方法保證結(jié)果的精確性與計算的高效性，最終求得各塊體的位移與變形（Shi，2001；Wu et al.，2012）.DDA方法因收斂速度高以及模擬結(jié)果精確等優(yōu)勢，廣泛應(yīng)用于邊坡破壞（Sitar et al.，2005；Wu，2010；Wu and Chen，2011）、洞室開挖（鄔愛清等，2006；Law and Lam，2003）以及爆破（劉剛等，2003）等動力學(xué)模擬研究中.此外，DDA方法在模擬區(qū)域地殼變形與應(yīng)變積累領(lǐng)域也得到了初步應(yīng)用（陳兵等，2000；陳祖安等，2008；秦小軍等，2001；王輝等，2003；Liu et al.，1996）.

DDA方法中塊體單元之間的接觸界面是位移間斷面，界面滑動與否由Mohr－Coulomb準(zhǔn)則控制，當(dāng)剪應(yīng)力τ滿足公式（1）時，接觸狀態(tài)由鎖定變?yōu)槟Σ粱瑒?

式中σn為接觸邊界上的正應(yīng)力，φ和C分別是摩擦角和黏著力（Shi，1992，2001），摩擦系數(shù)μ（＝tanφ）在模擬過程中保持不變.由于巖土工程主要針對水壩、邊坡等百米尺度、巖性相對均一的研究目標(biāo)，采用常摩擦系數(shù)進(jìn)行數(shù)值模擬已足夠精確.然而，針對更大尺度的地學(xué)問題，實驗室測量與野外調(diào)查結(jié)果均表明，地震成核、震前滑動、同震破裂、余震和震后余滑以及震間期斷層的愈合過程中，斷層面的摩擦系數(shù)并非常數(shù)（Marone，1998）.因此，為合理解釋地震周期各階段發(fā)震斷層行為及剪應(yīng)力變化過程，需要在數(shù)值模擬時采用更符合實際的摩擦本構(gòu)定律（張龍等，2013）.

基于實驗室摩擦實驗結(jié)果得到的速度－狀態(tài)摩擦本構(gòu)定律能夠模擬斷層面上的成核過程以及發(fā)生在俯沖帶的周期性慢滑移，結(jié)合數(shù)值模擬方法還可直觀給出地震周期的輪廓 （Lapusta and Rice，2003；Liu and Rice，2007）.該定律中，斷層面的摩擦強(qiáng)度由滑動速度V和狀態(tài)變量θ共同描述，其基本形式如下（Dieterich，1979，1981；Ruina，1983）：

為了實現(xiàn)復(fù)雜介質(zhì)的地震孕育、發(fā)生及震后調(diào)整全過程的數(shù)值模擬研究，首先需要解決速度－狀態(tài)摩擦本構(gòu)定律與動力學(xué)數(shù)值模擬方法的融合問題.本文選擇三維DDA方法作為研究對象，首先推導(dǎo)了基于速度－狀態(tài)摩擦本構(gòu)定律的三維DDA模擬實用公式，隨后設(shè)計了滑動－保持－滑動實驗和速度步進(jìn)實驗兩算例以驗證模擬效果，并分析了數(shù)值模擬結(jié)果與實驗室摩擦實驗結(jié)果存在細(xì)微差別可能的原因.

2 計算方法

由公式（2）和公式（3）知，塊體接觸面上每一時間步的摩擦系數(shù)μ由當(dāng)前時步的速度V和狀態(tài)變量θ共同決定，其他參數(shù)（μ0，V0，Dc，a和b）均為常數(shù).下面詳細(xì)推導(dǎo)了改進(jìn)的三維DDA算法中計算摩擦系數(shù)的實用公式.

2.1 速度

應(yīng)用三維DDA方法求解塊體速度時，選擇動力學(xué)計算模式，即每一時步結(jié)束時的速度為下一時步開始時的速度.D＝ （uvwrxryrzεxεyεzγxyγyz γzx）T表示塊體變形量，包含塊體位移、旋轉(zhuǎn)和應(yīng)變.令Δt為時步長，Di和Di＋1為t和t＋Δt時刻的塊體變形量，采用向前差分格式（Shi，1992，2001），則有

在小步長情況下可假定每一時步內(nèi)的加速度為常數(shù)，即

令當(dāng)前時步初始變形量為0，即Di＝0，公式（4）略去高階無窮小并代入公式（5）得

勻加速直線運動速度計算表達(dá)式為

部分公民環(huán)保意識薄弱，碳信息披露可以促進(jìn)企業(yè)可持續(xù)發(fā)展，但企業(yè)對此認(rèn)識不夠深刻，這些都使碳披露難以發(fā)揮其對經(jīng)營管理的指導(dǎo)作用。這就需要社會廣泛宣傳低碳經(jīng)濟(jì)理念，提高公民的環(huán)保意識，進(jìn)而逐步上升到企業(yè)管理層面，使管理者提高環(huán)保責(zé)任意識，加大對碳信息的重視。通過學(xué)習(xí)，掌握必要的碳披露知識，推動經(jīng)濟(jì)社會可持續(xù)發(fā)展。

將公式（6）代入公式（7），整理得速度迭代表達(dá)式

令Vi和Vi＋1分別表示，則（8）式可表示為

2.2 狀態(tài)變量

根據(jù)公式（2）或（3）的演化方程計算得到狀態(tài)變量θ.以公式（2）為例，將演化方程用4階Runge－Kutta方法離散，即

式中θi和θi＋1分別為t和t＋Δt時刻的狀態(tài)變量，Kj（j＝1，2，3，4）為中間變量.初始時刻的θ值由公式（2）（或（3））的主方程求得.

根據(jù)公式（9）和（10）分別計算得到下一時步的速度V和狀態(tài)變量θ，代入公式（2）（或（3））可計算得到下一時步的摩擦系數(shù)μ，完成該時間步計算.將上述算法嵌入三維DDA程序中，建立了應(yīng)用速度－狀態(tài)摩擦本構(gòu)定律的三維DDA方法.

3 數(shù)值模擬

為驗證結(jié)合速度－狀態(tài)摩擦本構(gòu)定律的三維DDA方法在計算地震周期中斷層行為問題時的精確性，本文設(shè)計了滑動－保持－滑動實驗和速度步進(jìn)實驗兩個算例.

3.1 滑動－保持－滑動實驗

研究表明，相比于其他形式的本構(gòu)定律，Dieterich law能夠準(zhǔn)確描述應(yīng)力降與保持時間的關(guān)系，即靜摩擦的時間依賴性（Dieterich，1978）.最初，實驗室做滑動－保持－滑動實驗使用“三明治”直剪實驗機(jī)（圖1a），Scholz（2002）將其簡化為圖1b形式.

圖1 （a）“三明治”直接剪切實驗機(jī)；（b）Scholz簡化的直接剪切實驗機(jī)Fig.1 （a）Schematic diagram of“sandwich”direct shear apparatus；（b）Direct shear apparatus modifid by Scholz

利用改進(jìn)的三維DDA方法數(shù)值模擬該實驗，需要設(shè)計形狀合適的模型并合理選擇材料的物理參數(shù).考慮到摩擦強(qiáng)度變化與接觸面積大小無關(guān)，為簡化模擬，對模型幾何形狀略作修改，如圖2所示.模型的介質(zhì)參數(shù)、數(shù)值控制參數(shù)及摩擦本構(gòu)參數(shù)見表1.

表1 滑動－保持－滑動實驗數(shù)值參數(shù)表Table 1 Numerical parameters of slide－h(huán)old－slide tests for 3DDDA forward modeling simulation

為與實驗室摩擦實驗數(shù)據(jù)（Beeler et al.，1994）進(jìn)行對比，本文分別做了保持時間為1s、10s、100s、1000s和10000s的滑動－保持－滑動數(shù)值實驗.限于篇幅，僅以圖3為例給出保持時間為10s的滑動－保持－滑動數(shù)值實驗?zāi)M結(jié)果.

圖2 滑動－保持－滑動實驗的模型幾何Fig.2 Numerical model for slide－h(huán)old－slide tests

圖3 滑動－保持－滑動數(shù)值實驗?zāi)M結(jié)果（保持時間為10s）Fig.3 Numerical results of slide－h(huán)old－slide tests（hold period is 10s）

如圖3所示，前5s滑塊在荷載力和摩擦力的共同作用下做勻速直線運動，即“滑動”過程；5～15s為“保持”階段，即令荷載速度為零，滑塊在摩擦力的作用下做減速運動，此時摩擦系數(shù)也發(fā)生相應(yīng)變化；15～20s仍為“滑動”階段，15s時重載初始時刻荷載速度，可明顯觀察到靜摩擦系數(shù)瞬間達(dá)到峰值，隨著時間的增加不斷減小，最終恢復(fù)到初始水平.不同保持時長的滑動－保持－滑動數(shù)值實驗?zāi)Σ料禂?shù)演化過程與該結(jié)果類似.

圖4給出了靜摩擦變化量Δμs（峰值摩擦系數(shù)與穩(wěn)態(tài)值的差值）與保持時間的關(guān)系.實驗室摩擦實驗結(jié)果與地震學(xué)對斷層愈合速率的估計均表明，準(zhǔn)靜態(tài)接觸時摩擦愈合量Δμs與保持時間的對數(shù)呈線性關(guān)系（Beeler et al.，1994），即所謂靜摩擦?xí)r間依賴性.圖4表明，模擬結(jié)果與實驗室測量結(jié)果（Beeler et al.，1994）基本吻合.當(dāng)保持時間為101～104s時，模擬數(shù)據(jù)與實驗數(shù)據(jù)對數(shù)擬合后呈線性變化且斜率相同，擬合優(yōu)度均為0.999，將兩部分?jǐn)?shù)據(jù)混合再進(jìn)行對數(shù)擬合的擬合優(yōu)度依然高達(dá)0.997.由此得出結(jié)論，應(yīng)用速度－狀態(tài)摩擦本構(gòu)定律的三維DDA方法能夠比較準(zhǔn)確地描述靜摩擦的時間依賴性.另外，摩擦愈合速率β（＝Δμs每十年的變化量）.在速度接近零時可由摩擦本構(gòu)參數(shù)b衡量（Beeler et al.，1994）.換言之，用該方法可求出Dieterich law中本構(gòu)參數(shù)b.圖4中保持時間為101～104s的數(shù)值模擬結(jié)果斜率為0.00856，與摩擦本構(gòu)參數(shù)b（＝0.009）十分接近，相對誤差約為5%，表明應(yīng)用速度－狀態(tài)摩擦本構(gòu)定律的三維DDA方法進(jìn)行數(shù)值模擬可以基本滿足精度要求.

圖4 滑動－保持－滑動實驗的數(shù)值結(jié)果與實驗室結(jié)果圖中三角形表示三維DDA的數(shù)值實驗結(jié)果，圓點和菱形表示Beeler等（1994）實驗室摩擦實驗結(jié)果.實線代表用Dieterich law對滑動－保持－滑動實驗數(shù)值模擬的結(jié)果.Fig.4 Comparison between numerical results and laboratory measurements for slide－h(huán)old－slide tests In this figure，numerical results performed by 3DDDA are showed as triangles.Circles and diamonds represent the laboratory data measured by Beeler et al.（1994）.The curve represents numerical results performed by traditional integral method with Dieterich law.

3.2 速度步進(jìn)實驗

速度步進(jìn)實驗可用來研究擴(kuò)容作用及斷層泥對動摩擦速度依賴性的影響 （Dieterich and Kilgore，1994）.研究表明，Ruina law能夠準(zhǔn)確描述摩擦系數(shù)對速度變化的響應(yīng)，即動摩擦的速度依賴性（何昌榮，1999；Dieterich and Kilgore，1994），因而在該算例中選擇公式 （3）計算摩擦系數(shù).本實驗選擇與滑動－保持－滑動實驗相同的幾何模型（圖2），介質(zhì)參數(shù)、數(shù)值控制參數(shù)及摩擦本構(gòu)參數(shù)見表2.

模擬結(jié)果如圖5b所示，荷載點速度為10μm／s，滑塊穩(wěn)定滑動；運動25μm后將荷載點速度突然降至1μm／s，運動一段距離后滑塊恢復(fù)穩(wěn)定滑動狀態(tài)，重復(fù)上述過程一次.

圖5表明，應(yīng)用速度－狀態(tài)摩擦本構(gòu)定律的三維DDA方法能夠有效模擬動摩擦的速度依賴性.選擇與摩擦實驗相同的加載方式和加載條件，三維DDA方法能夠比較準(zhǔn)確地描述摩擦系數(shù)對速度由高值到低值再到高值的響應(yīng)，曲線形態(tài)與實驗結(jié)果基本吻合.值得注意的是，數(shù)值模擬曲線中摩擦系數(shù)在速度突變時的響應(yīng)較實驗曲線銳利，通過設(shè)置不同的數(shù)值彈簧剛度及介質(zhì)材料參數(shù)，并未使得模擬曲線峰尖變緩.分析表明，物理實驗的加載（速度變化）做不到數(shù)值模擬那樣的瞬時響應(yīng)，因而存在這種差別是合理的.

表2 速度步進(jìn)實驗數(shù)值參數(shù)表Table 2 Numerical parameters of velocity stepping tests for 3DDDA forward modeling simulation

圖5 速度步進(jìn)摩擦實驗結(jié)果與數(shù)值模擬結(jié)果Fig.5 Comparison between numerical results and laboratory measurements for velocity stepping tests

4 結(jié)論

本文通過公式推導(dǎo)、模型驗證等過程，拓展了基于速度－狀態(tài)摩擦本構(gòu)定律的三維DDA方法，取得如下認(rèn)識和結(jié)論：

（1）從速度－狀態(tài)摩擦本構(gòu)定律的基本形式出發(fā)，通過推導(dǎo)給出了融合該定律的三維DDA方法計算摩擦系數(shù)的實用公式，拓展了三維DDA方法.

（2）滑動－保持－滑動實驗和速度步進(jìn)實驗?zāi)M結(jié)果表明，應(yīng)用速度－狀態(tài)摩擦本構(gòu)定律的三維DDA方法可以比較準(zhǔn)確地模擬實驗中觀察到的靜摩擦?xí)r間依賴性和動摩擦速度依賴性，進(jìn)而表明該方法模擬接觸面剪應(yīng)力動態(tài)變化是有效的.

（3）同原有方法相比，改進(jìn)后的三維DDA方法解決了應(yīng)用于大尺度構(gòu)造應(yīng)力演化和區(qū)域應(yīng)變積累的基本問題，即Coulomb定律中的常摩擦系數(shù)不能反映地震周期中斷層面剪應(yīng)力變化.然而，改進(jìn)的三維DDA方法在模擬中依然存在由于不連續(xù)面剛度的不確定性，發(fā)生嵌入深淺不一的問題，往往會影響求解接觸力或應(yīng)力的精度.相信隨著三維DDA凹塊體接觸算法不斷完善，本文這種改進(jìn)將有助于今后與觀測資料（GPS等）相結(jié)合，定量描述區(qū)域地殼變形特征.

致謝感謝石根華先生、九州大學(xué)陳光齊教授提供最新版本的三維DDA代碼，感謝中國地震局地質(zhì)研究所何昌榮研究員關(guān)于數(shù)值模擬技術(shù)的有益指導(dǎo)，感謝兩位匿名審稿專家為完善本文提出的寶貴意見.

Beeler N M，Tullis T E，Weeks J D.1994.The roles of time and displacement in the evolution effect in rock friction.Geophys.Res.Lett.，21（18）：1987－1990.

Chen B，Jiang Z S，Wang S X，et al.2000.Preliminary study on block movement and its stress field by discontinuous deformation analysis.CrustalDeformationandEarthquake（in Chinese），20（1）：38－42.

Chen Z A，Lin B H，Bai W M.2008.3－D numerical simulation on influence of 1997Mani earthquake occurrence to stability of tectonic blocks system in Qingzang and Chuandian zone using DDA＋FEM method.ChineseJ.Geophys.（in Chinese），51（5）：1422－1430，doi：10.3321／j.issn：0001－5733.2008.05.015.Dieterich J H.1978.Time－dependence friction and the mechanics of stick－slip.PureAppl.Geophys.，116（4－5）：790－806.

Dieterich J H.1979.Modeling of rock friction：1.Experimental results and constitutive equations.J.Geophys.Res.，84（B5）：2161－2168.

Dieterich J H.1981.Constitutive properties of faults with simulated gouge.／／Carter N L，F(xiàn)riedman M，Logan J M eds.Mechanical Behavior of Crustal Rocks：The Handin Volume.Washington，DC：Am.Geophys.Union：24：103－120.

Dieterich J H，Kilgore B D.1994.Direct observation of frictional contacts：new insights for state－dependent properites.Pure Appl.Geophys.，143（1－3）：283－302.

He C R.1999.Comparing two types of rate and state dependent friction laws.SeismologyandGeology（in Chinese），21（2）：137－146.

Lapusta N，Rice J R.2003.Nucleation and early seismic propagation of small and large events in a crustal earthquake model.J.Geophys.Res.，108（B4），doi：10.1029／2001JB00793.Law H K，Lam I P.2003.Evaluation of seismic performance for tunnel retrofit project.J.Geotech.Geoenviron.Eng.，129（7）：575－589.

Liu G，Shu D Q，Jiang Q H.2003.Application of discontinuous deformation analysis method in analyzing shaft well stabilization.Blasting（in Chinese），20（4）：38－44.

Liu L B，Linde A T，Sacks I S，et al.1996.A seismic fault slip and block deformation in north China.PureAppl.Geophys.，146（3－4）：717－740.

Liu Y J，Rice J R.2007.Slow slip predictions based on granite and gabbro friction data compared to GPS measurements in northern Cascadia.J.Geophys.Res.，114（B9）：B09407，doi：10.1029／2008JB006142.

Marone C.1998. Laboratory－derived friction laws and their application to seismic faulting.Annu.Rev.EarthPlanet.Sci.，26：643－696.

Qin X J，Zhou S Y，Zhao L Q，et al.2001.Research on horizontal crustal movement from 1994to 1998in Jiashi and the northeast to Pamir ＆its relation to strong earthquake swarm on the basis of GPS data.CrustalDeformationandEarthquake（in Chinese），21（4）：1－7.

Ruina A L.1983.Slip instability and state variable friction law.J.Geophys.Res.，88（B12）：10359－10370.

Scholz C H.2002.The Mechanics of Earthquakes and Faulting.2nd ed.New York：Cambridge University Press.

Shi G H.1992.Discontinuous deformation analysis：A new numerical model for the statistics and dynamics of deformable block structures.Engng.Comp.，9（2）：157－168.

Shi G H.2001. Three－dimensional discontinuous deformation analysis.／／Elsworth D ed.Proceedings of the 38th U.S.Rock Mechanics Symposium.Taylor ＆Francis，Washington，D.C.，USA，1421－1428.

Sitar N，MacLaughlin M M，Doolin D M.2005.Influence of kinematics on landslide mobility and failure mode.J.Geotech.Geoenviron.Eng.，131（6）：716－728.

Wang H，Zhang G M，Wu Y，et al.2003.The deformation of active crustal－blocks on the Chinese Mainland and its relation with seismic activity.EarthquakeResearchinChina（in Chinese），19（3）：243－254.

Wu A Q，Ding X L，Chen S H，et al.2006.Researches on deformation and failure characteristics of an underground powerhowse with complicated gelolgical conditions by DDA method.ChineseJ.RockMech.Engng.（in Chinese），25（1）：1－8.

Wu J H.2010.Seismic landslide simulations in discontinuous deformation analysis.Comp.Geotech.，37（5）：594－601.

Wu J H，Chen C H.2011.Application of DDA to simulate characteristics of the Tsaoling landslide.Comp.Geotech.，38（5）：741－750.

Wu Y Q，Chen G Q，Jiang Z S，et al.2012.The algorithm of simplex integration in three－dimension and its characteristic analysis.Int.J.Adv.Comp.Tech.，4（10）：246－256.

Zhang L，Jiang Z S，Wu Y Q.2013.Review of rate－and statedependent friction laws and their applications to seismic faulting.ProgressinGeophysics（in Chinese），28（5）：2352－2362，doi：10.6038／pg20130517.

附中文參考文獻(xiàn)

陳兵，江在森，王雙緒等.2000.利用非連續(xù)變形數(shù)值方法研究塊體運動及其應(yīng)力場初探.地殼形變與地震，20（1）：38－42.

陳祖安，林邦慧，白武明.2008.1997年瑪尼地震對青藏川滇地區(qū)構(gòu)造塊體系統(tǒng)穩(wěn)定性影響的三維DDA＋FEM方法數(shù)值模擬.地球物理學(xué)報，51（5）：1422－1430，doi：10.3321／j.issn：0001－5733.2008.05.015.

何昌榮.1999.兩種摩擦本構(gòu)關(guān)系的對比研究.地震地質(zhì)，21（2）：137－146.

劉剛，舒大強(qiáng)，姜清輝.2003.應(yīng)用DDA方法分析豎井爆振穩(wěn)定問題.爆破，20（4）：38－44.

秦小軍，周碩愚，趙齊樂等.2001.依據(jù)GPS數(shù)據(jù)研究伽師及帕米爾東北側(cè)1994—1998年地殼水平運動及其與強(qiáng)震群的關(guān)系.地殼形變與地震，21（4）：1－7.

王輝，張國民，吳云等.2003.中國大陸活動地塊變形與地震活動的關(guān)系.中國地震，19（3）：243－254.

鄔愛清，丁秀麗，陳勝宏等.2006.DDA方法在復(fù)雜地質(zhì)條件下地下廠房圍巖變形與破壞特征分析中的應(yīng)用研究.巖石力學(xué)與工程學(xué)報，25（1）：1－8.

張龍，江在森，武艷強(qiáng).2013.速度－狀態(tài)摩擦本構(gòu)定律及其在地震斷層中的應(yīng)用研究進(jìn)展.地球物理學(xué)進(jìn)展，28（5）：2352－2362，doi：10.6038／pg20130517.