地球物理学报  2015, Vol. 58 Issue (2): 474-480   PDF    
基于速度-状态摩擦本构定律的三维DDA方法
张龙1, 江在森1, 武艳强1,2, 邹镇宇1,3, 刘晓霞1,3, 魏文薪1    
1. 中国地震局地震预测研究所地震预测重点实验室, 北京 100036;
2. 中国地震局第一监测中心, 天津 300180;
3. 中国地震局地质研究所地震动力学国家重点实验室, 北京 100029
摘要:原有三维非连续变形分析(DDA)方法采用常摩擦系数的Mohr-Coulomb定律作为切向破坏准则,然而当描述更大尺度构造块体的运动与变形时,常摩擦系数不再适用.速度-状态摩擦本构定律能够定量描述地震周期各阶段断层面剪应力变化,解释发震断层行为.本文将速度-状态摩擦定律与三维DDA方法相结合,首先推导了计算摩擦系数的实用公式,随后通过滑动-保持-滑动实验与速度步进实验算例对改进的三维DDA方法进行了验证.结果表明,应用速度-状态摩擦本构定律的三维DDA方法能够比较准确地模拟静摩擦的时间依赖性与动摩擦的速度依赖性,解决了将三维DDA方法在地学中应用的基本问题.
关键词三维DDA     速度-状态摩擦本构定律     数值模拟    
Development of the 3D DDA method with rate- and state-friction laws
ZHANG Long1, JIANG Zai-Sen1, WU Yan-Qiang1,2, ZOU Zhen-Yu1,3, LIU Xiao-Xia1,3, WEI Wen-Xin1    
1. Key Laboratory of Earthquake Prediction, Institute of Earthquake Science, China Earthquake Administration, Beijing 100036, China;
2. First Crust Monitoring and Application Center, China Earthquake Administration, Tianjin 300180, China;
3. State Key Laboratory of Earthquake Dynamics, Institute of Geology, China Earthquake Administration, Beijing 100029, China
Abstract: The elastic medium in lithosphere should be considered as discontinuities since the complex tectonic background has produced many active faults as separation. The Mohr-Coulomb joint failure criterion with a constant friction coefficient, adopted in the original 3D discontinuous deformation analysis (DDA) method, cannot meet the requirement of highly accurate calculations for motion and deformation of tectonic block systems.
The rate- and state-friction laws, which are capable of reproducing virtually the entire range of observed fault behaviors, are combined into the 3D DDA method. Firstly, the formula of computing coefficient of friction on the interface with rate- and state friction laws is derived. In order to calculate the value of state variable θ, slip velocity V and friction coefficient μ in every time step, a first-order differential equation about V is deduced. The increment of V is calculated by the second-order Taylor series expansion in our scheme. Secondly, the evolution law is determined by the Runge-Kutta scheme with adaptive step-size control. The friction submatrix, which consists of discrete forms of V and θ, is rewritten and then combined into the 3D DDA method. Finally, on the basis of reasonable geometry and mechanical properties of the numerical model, slide-hold-slide tests and velocity stepping tests are designed to examine the accuracy of the modified method.
Suits of numerical slide-hold-slide tests are performed using hold time from 1 to 10000 seconds and then we take the numerical test of the 10-second hold time as an example to make a brief illustration. A rigid block moves along the base in uniform linear motion at the first five seconds, because of the equality between friction and point loading. At the fifth second, the point loading is set to zero, and the friction strength drops immediately. After ten seconds, the block is reset as the previous loading. Strength increases, reaches a peak value and returns to its previous steady-state value subsequently. The numerical results are consistent with the laboratory data measured by Beeler et al. (1994) with a 0.999-goodness of fits. Even though both data sets are mixed together, there is also a 0.997-goodness of fits. In addition, the friction healing rate β can be evaluated by the friction parameter b when the slip rate approximates to 0. The slope of these numerical data comes close to the b value of assumption with a 5% relative error. Besides, velocity stepping tests are also modeled using the improved 3D DDA. A slider keeps a steady sliding under a loading rate of 10 μm/s with a distance of 25 μm, and then decreases to 1 μm/s during a characteristic distance. The results show a strong velocity dependence of friction which is consistent with laboratory data.
Comparison between numerical results and laboratory data shows that the 3D DDA method in combination with rata- and state-dependent friction laws is capable of simulating velocity dependence of sliding friction and time dependence of static friction, which resolves a basic problem when using 3D DDA in geodynamics research. The improved 3D DDA method still has the limitation such as unbalanced embedment due to the uncertainty of stiffness on the discontinuous interface, which leads to inexact results of contact force determination. In the future, this modification can be used in quantitative simulation of regional crustal deformation in combination with observations, such as GPS data.
Key words: 3D DDA     Rate- and state-friction laws     Numerical simulation    

1 引言

非连续变形分析(DDA)方法是岩土工程领域模拟节理岩体运动与变形的一种数值解析方法(Shi,1992).该方法将不连续面分割的岩体视为块体集合,以最小势能法建立总体平衡方程,利用单纯形积分与开闭迭代方法保证结果的精确性与计算的高效性,最终求得各块体的位移与变形(Shi,2001Wu et al., 2012).DDA方法因收敛速度高以及模拟结果精确等优势,广泛应用于边坡破坏(Sitar et al., 2005Wu,2010Wu and Chen, 2011)、洞室开挖(邬爱清等,2006Law and Lam, 2003)以及爆破(刘刚等,2003)等动力学模拟研究中.此外,DDA方法在模拟区域地壳变形与应变积累领域也得到了初步应用(陈兵等,2000陈祖安等,2008秦小军等,2001王辉等,2003Liu et al., 1996).

DDA方法中块体单元之间的接触界面是位移间断面,界面滑动与否由Mohr-Coulomb准则控制,当剪应力τ满足公式(1)时,接触状态由锁定变为摩擦滑动.

式中σn为接触边界上的正应力,φ和C分别是摩擦角和黏着力(Shi, 19922001),摩擦系数μ(=tanφ)在模拟过程中保持不变.由于岩土工程主要针对水坝、边坡等百米尺度、岩性相对均一的研究目标,采用常摩擦系数进行数值模拟已足够精确.然而,针对更大尺度的地学问题,实验室测量与野外调查结果均表明,地震成核、震前滑动、同震破裂、余震和震后余滑以及震间期断层的愈合过程中,断层面的摩擦系数并非常数(Marone,1998).因此,为合理解释地震周期各阶段发震断层行为及剪应力变化过程,需要在数值模拟时采用更符合实际的摩擦本构定律(张龙等,2013).

基于实验室摩擦实验结果得到的速度-状态摩擦本构定律能够模拟断层面上的成核过程以及发生在俯冲带的周期性慢滑移,结合数值模拟方法还可 直观给出地震周期的轮廓(Lapusta and Rice, 2003Liu and Rice, 2007). 该定律中,断层面的摩擦强度由滑动速度V和状态变量θ共同描述,其基本形式如下(Dieterich, 19791981Ruina,1983):

式中V0是参考速度,μ0是速度为V0时稳定状态下的摩擦系数(常数),a和b均为与实验有关的经验常数,Dc是临界滑动距离,即以某一速度稳定滑动到以另一速度稳定滑动过程中所需的滑动距离.

为了实现复杂介质的地震孕育、发生及震后调整全过程的数值模拟研究,首先需要解决速度-状态摩擦本构定律与动力学数值模拟方法的融合问题.本文选择三维DDA方法作为研究对象,首先推导了基于速度-状态摩擦本构定律的三维DDA模拟实用公式,随后设计了滑动-保持-滑动实验和速度步进实验两算例以验证模拟效果,并分析了数值模拟结果与实验室摩擦实验结果存在细微差别可能的原因. 2 计算方法

由公式(2)和公式(3)知,块体接触面上每一时间步的摩擦系数μ由当前时步的速度V和状态变量θ共同决定,其他参数0,V0,Dca和b)均为常数.下面详细推导了改进的三维DDA算法中计算摩擦系数的实用公式. 2.1 速度

应用三维DDA方法求解块体速度时,选择动力学计算模式,即每一时步结束时的速度为下一时步开始时的速度. D =(u v w rx ry rz εx εy εz γxy γyz γzx)T表示块体变形量,包含块体位移、旋转和应变.令Δt为时步长,Di和Di+1为t和t+Δt时刻的块体变形量,采用向前差分格式(Shi, 19922001),则有

在小步长情况下可假定每一时步内的加速度为常数,即

令当前时步初始变形量为0,即Di=0,公式(4)略去高阶无穷小并代入公式(5)得

匀加速直线运动速度计算表达式为

将公式(6)代入公式(7),整理得速度迭代表达式

令Vi和Vi+1分别表示ii+1,则(8)式可表示为

2.2 状态变量

根据公式(2)或(3)的演化方程计算得到状态变量θ.以公式(2)为例,将演化方程用4阶Runge-Kutta方法离散,即

式中θi和θi+1分别为t和t+Δt时刻的状态变量,Kj(j=1,2,3,4)为中间变量.初始时刻的θ值由公式(2)(或(3))的主方程求得.

根据公式(9)和(10)分别计算得到下一时步的速度V和状态变量θ,代入公式(2)(或(3))可计算得到下一时步的摩擦系数μ,完成该时间步计算.将上述算法嵌入三维DDA程序中,建立了应用速度-状态摩擦本构定律的三维DDA方法. 3 数值模拟

为验证结合速度-状态摩擦本构定律的三维DDA方法在计算地震周期中断层行为问题时的精确性,本文设计了滑动-保持-滑动实验和速度步进实验两个算例. 3.1 滑动-保持-滑动实验

研究表明,相比于其他形式的本构定律,Dieterich law 能够准确描述应力降与保持时间的关系,即静摩擦的时间依赖性(Dieterich,1978).最初,实验室做滑动-保持-滑动实验使用“三明治”直剪实验机(图 1a),Scholz(2002)将其简化为图 1b形式.

图 1(a)“三明治”直接剪切实验机;(b)Scholz简化的直接剪切实验机 Fig. 1(a)Schematic diagram of “s and wich” direct shear apparatus;(b)Direct shear apparatus modifid by Scholz

利用改进的三维DDA方法数值模拟该实验,需要设计形状合适的模型并合理选择材料的物理参数.考虑到摩擦强度变化与接触面积大小无关,为简化模拟,对模型几何形状略作修改,如图 2所示.模 型的介质参数、数值控制参数及摩擦本构参数见表 1.

图 2 滑动-保持-滑动实验的模型几何 Fig. 2 Numerical model for slide-hold-slide tests

表 1 滑动-保持-滑动实验数值参数表 Table 1 Numerical parameters of slide-hold-slide tests for 3D DDA forward modeling simulation

为与实验室摩擦实验数据(Beeler et al., 1994)进行对比,本文分别做了保持时间为1 s、10 s、100 s、1000 s和10000 s的滑动-保持-滑动数值实验.限于篇幅,仅以图 3为例给出保持时间为10 s的滑动-保持-滑动数值实验模拟结果.

图 3 滑动-保持-滑动数值实验模拟结果(保持时间为10 s) Fig. 3 Numerical results of slide-hold-slide tests(hold period is 10 s)

图 3所示,前5 s滑块在荷载力和摩擦力的共同作用下做匀速直线运动,即“滑动”过程;5~15 s为“保持”阶段,即令荷载速度为零,滑块在摩擦力的作用下做减速运动,此时摩擦系数也发生相应变化;15~20 s仍为“滑动”阶段,15 s时重载初始时刻荷载速度,可明显观察到静摩擦系数瞬间达到峰值,随着时间的增加不断减小,最终恢复到初始水平.不同保持时长的滑动-保持-滑动数值实验摩擦系数演化过程与该结果类似.

图 4给出了静摩擦变化量Δμs(峰值摩擦系数与稳态值的差值)与保持时间的关系.实验室摩擦实验结果与地震学对断层愈合速率的估计均表明,准静态接触时摩擦愈合量Δμs与保持时间的对数呈线性关系(Beeler et al., 1994),即所谓静摩擦时间依赖性.图 4表明,模拟结果与实验室测量结果(Beeler et al., 1994)基本吻合.当保持时间为101~104 s时,模拟数据与实验数据对数拟合后呈线性变化且斜率相同,拟合优度均为0.999,将两部分数据混合再进行对数拟合的拟合优度依然高达0.997.由此得出结论,应用速度-状态摩擦本构定律的三维DDA方法能够比较准确地描述静摩擦的时间依赖性.另外,摩擦愈合速率β(=Δμs每十年的变化量).在速度接近零时可由摩擦本构参数b衡量(Beeler et al., 1994).换言之,用该方法可求出 Dieterich law中本构参数b.图 4中保持时间为101~104 s 的数值模拟结果斜率为0.00856,与摩擦本构参数b(=0.009)十分接近,相对误差约为5%,表明应用速度-状态摩擦本构定律的三维DDA方法进行数值模拟可以基本满足精度要求.

图 4 滑动-保持-滑动实验的数值结果与实验室结果
图中三角形表示三维DDA的数值实验结果,圆点和菱形表示Beeler等(1994)实验室摩擦实验结果.实线代表用Dieterich law 对滑动-保持-滑动实验数值模拟的结果.
Fig. 4 Comparison between numerical results and laboratory measurements for slide-hold-slide tests
In this figure,numerical results performed by 3D DDA 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 速度步进实验

速度步进实验可用来研究扩容作用及断层泥对动摩擦速度依赖性的影响(Dieterich and Kilgore, 1994).研究表明,Ruina law 能够准确描述摩擦系数对速度变化的响应,即动摩擦的速度依赖性(何昌荣,1999Dieterich and Kilgore, 1994),因而在该算例中选择公式(3)计算摩擦系数.本实验选择与滑动-保持-滑动实验相同的几何模型(图 2),介质参数、数值控制参数及摩擦本构参数见表 2.

表 2 速度步进实验数值参数表 Table 2 Numerical parameters of velocity stepping tests for 3D DDA forward modeling simulation

模拟结果如图 5b所示,荷载点速度为10 μm/s,滑块稳定滑动;运动25 μm后将荷载点速度突然降至1 μm/s,运动一段距离后滑块恢复稳定滑动状态,重复上述过程一次.

图 5 速度步进摩擦实验结果与数值模拟结果 Fig. 5 Comparison between numerical results and laboratory measurements for velocity stepping tests

图 5表明,应用速度-状态摩擦本构定律的三维 DDA方法能够有效模拟动摩擦的速度依赖性.选择与摩擦实验相同的加载方式和加载条件,三维DDA方法能够比较准确地描述摩擦系数对速度由高值到低值再到高值的响应,曲线形态与实验结果基本吻 合.值得注意的是,数值模拟曲线中摩擦系数在速度突变时的响应较实验曲线锐利,通过设置不同的数值弹簧刚度及介质材料参数,并未使得模拟曲线峰尖变缓.分析表明,物理实验的加载(速度变化)做不到数值模拟那样的瞬时响应,因而存在这种差别是合理的. 4 结论

本文通过公式推导、模型验证等过程,拓展了基于速度-状态摩擦本构定律的三维DDA方法,取得如下认识和结论:

(1)从速度-状态摩擦本构定律的基本形式出发,通过推导给出了融合该定律的三维DDA方法 计算摩擦系数的实用公式,拓展了三维DDA方法.

(2)滑动-保持-滑动实验和速度步进实验模拟结果表明,应用速度-状态摩擦本构定律的三维DDA方法可以比较准确地模拟实验中观察到的静摩擦时间依赖性和动摩擦速度依赖性,进而表明该方法模拟接触面剪应力动态变化是有效的.

(3)同原有方法相比,改进后的三维DDA方法解决了应用于大尺度构造应力演化和区域应变积累的基本问题,即Coulomb定律中的常摩擦系数不能反映地震周期中断层面剪应力变化.然而,改进的 三维DDA方法在模拟中依然存在由于不连续面刚度的不确定性,发生嵌入深浅不一的问题,往往会影响求解接触力或应力的精度.相信随着三维DDA凹块体接触算法不断完善,本文这种改进将有助于今后与观测资料(GPS等)相结合,定量描述区域地壳变形特征.

致谢 感谢石根华先生、九州大学陈光齐教授提供最新版本的三维DDA代码,感谢中国地震局地质研究所何昌荣研究员关于数值模拟技术的有益指导,感谢两位匿名审稿专家为完善本文提出的宝贵意见.

参考文献
[1] 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.
[2] Chen B, Jiang Z S, Wang S X, et al. 2000. Preliminary study on block movement and its stress field by discontinuous deformation analysis. Crustal Deformation and Earthquake (in Chinese), 20(1): 38-42.
[3] Chen Z A, Lin B H, Bai W M. 2008. 3-D numerical simulation on influence of 1997 Mani earthquake occurrence to stability of tectonic blocks system in Qingzang and Chuandian zone using DDA+FEM method. Chinese J. Geophys. (in Chinese), 51(5): 1422-1430, doi: 10.3321/j.issn:0001-5733.2008.05.015.
[4] Dieterich J H. 1978. Time-dependence friction and the mechanics of stick-slip. Pure Appl. Geophys., 116(4-5): 790-806.
[5] Dieterich J H. 1979. Modeling of rock friction: 1. Experimental results and constitutive equations. J. Geophys. Res., 84(B5): 2161-2168.
[6] Dieterich J H. 1981. Constitutive properties of faults with simulated gouge. //Carter N L, Friedman M, Logan J M eds. Mechanical Behavior of Crustal Rocks: The Handin Volume. Washington, DC: Am. Geophys. Union: 24: 103-120.
[7] 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.
[8] He C R. 1999. Comparing two types of rate and state dependent friction laws. Seismology and Geology (in Chinese), 21(2): 137-146.
[9] 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.
[10] Law H K, Lam I P. 2003. Evaluation of seismic performance for tunnel retrofit project. J. Geotech. Geoenviron. Eng., 129(7): 575-589.
[11] 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.
[12] Liu L B, Linde A T, Sacks I S, et al. 1996. A seismic fault slip and block deformation in north China. Pure Appl. Geophys., 146(3-4): 717-740.
[13] 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.
[14] Marone C. 1998. Laboratory-derived friction laws and their application to seismic faulting. Annu. Rev. Earth Planet. Sci., 26: 643-696.
[15] Qin X J, Zhou S Y, Zhao L Q, et al. 2001. Research on horizontal crustal movement from 1994 to 1998 in Jiashi and the northeast to Pamir & its relation to strong earthquake swarm on the basis of GPS data. Crustal Deformation and Earthquake (in Chinese), 21(4): 1-7.
[16] Ruina A L. 1983. Slip instability and state variable friction law. J. Geophys. Res., 88(B12): 10359-10370.
[17] Scholz C H. 2002. The Mechanics of Earthquakes and Faulting. 2nd ed. New York: Cambridge University Press.
[18] 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.
[19] 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.
[20] 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.
[21] 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. Earthquake Research in China (in Chinese), 19(3): 243-254.
[22] 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. Chinese J. Rock Mech. Engng. (in Chinese), 25(1): 1-8.
[23] Wu J H. 2010. Seismic landslide simulations in discontinuous deformation analysis. Comp. Geotech., 37(5): 594-601.
[24] Wu J H, Chen C H. 2011. Application of DDA to simulate characteristics of the Tsaoling landslide. Comp. Geotech., 38(5): 741-750.
[25] 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.
[26] Zhang L, Jiang Z S, Wu Y Q. 2013. Review of rate-and state-dependent friction laws and their applications to seismic faulting. Progress in Geophysics (in Chinese), 28(5): 2352-2362, doi: 10.6038/pg20130517.
[27] 陈兵, 江在森, 王双绪等. 2000. 利用非连续变形数值方法研究块体运动及其应力场初探. 地壳形变与地震, 20(1): 38-42.
[28] 陈祖安, 林邦慧, 白武明. 2008. 1997年玛尼地震对青藏川滇地区构造块体系统稳定性影响的三维DDA+FEM方法数值模拟. 地球物理学报, 51(5): 1422-1430, doi: 10.3321/j.issn:0001-5733.2008.05.015.
[29] 何昌荣. 1999. 两种摩擦本构关系的对比研究. 地震地质, 21(2): 137-146.
[30] 刘刚, 舒大强, 姜清辉. 2003. 应用DDA方法分析竖井爆振稳定问题. 爆破, 20(4): 38-44.
[31] 秦小军, 周硕愚, 赵齐乐等. 2001. 依据GPS数据研究伽师及帕米尔东北侧1994—1998年地壳水平运动及其与强震群的关系. 地壳形变与地震, 21(4): 1-7.
[32] 王辉, 张国民, 吴云等. 2003. 中国大陆活动地块变形与地震活动的关系. 中国地震, 19(3): 243-254.
[33] 邬爱清, 丁秀丽, 陈胜宏等. 2006. DDA方法在复杂地质条件下地下厂房围岩变形与破坏特征分析中的应用研究. 岩石力学与工程学报, 25(1): 1-8.
[34] 张龙, 江在森, 武艳强. 2013. 速度-状态摩擦本构定律及其在地震断层中的应用研究进展. 地球物理学进展, 28(5): 2352-2362, doi: 10.6038/pg20130517.