基于颗粒流程序的非定常西原体模型

工程科学学报，第37卷，第7期：831838,2015年7月 Chinese Journal of Engineering,Vol.37,No.7:831-838,July 2015 DOI:10.13374/j.issn2095-9389.2015.07.002:http://journals.ustb.edu.cn 基于颗粒流程序的非定常西原体模型 杨振伟，金爱兵四，高永涛，王凯，孙浩 北京科技大学土木与环境工程学院，北京100083 ☒通信作者，Email:jinaibing@vip.sina.com 摘要为实现岩石试样蠕变全过程的准确模拟，并从细观角度探究蠕变过程中微裂隙的发生和发展规律，在二维颗粒流程 序(P℉℃)中开发出具有黏弹塑性特征的西原体流变接触本构模型，进一步提出包含两种非定常元件的非定常西原体模型， 推导了模型本构关系和蠕变方程.在P℉C2”中调用自定义西原体流变模型，通过参数调试，获得与真实试样具有相同强度特 性的数值试样.以室内单轴压缩矯变试验数据为基础，在Mb中对模型非定常参数进行拟合反演分析.在此基础上，进行 单轴压缩蠕变试验的模拟，计算过程中分别采用定常和非定常两种模型，并对微裂隙进行监测.对比分析结果表明：定常模 型仅适用于衰减和稳定蠕变阶段：非定常模型也可用于描述加速蠕变阶段，从而准确模拟蠕变全过程：加速蠕变阶段主要是 由微裂隙的加速发展而产生，加速蠕变将导致试样剪切破坏 关键词岩石力学：蠕变：本构模型：颗粒流程序：裂隙 分类号TU452 Non-stationary Nishihara model in the particle flow code YANG Zhen-wei,JIN Ai-bing,GAO Yong-tao,WANG Kai,SUN Hao School of Civil and Environmental Engineering,University of Science and Technology Beijing,Beijing 100083,China Corresponding author,E-mail:jinaibing@vip.sina.com ABSTRACT In order to accurately simulate the whole process of rock creep and explore the generation and development of mirco- fractures from the mesoevel,a visco-elastoplastic Nishihara rheological constitutive model is developed in the two-dimensional particle flow code (PFC2).Then a non-stationary Nishihara model including two non-stationary elements is put forward in this article.A con- stitutive equation and a creep equation are derived based on the non-stationary Nishihara model.The user-defined Nishihara constitu- tive model is called in PFC2,and a numerical sample whose strength properties are the same as those of a real rock specimen is ac- quired by parameter testing.Using the data of uniaxial creep tests in laboratory,the non-stationary parameters are back analyzed in Matlab.At last,uniaxial creep tests are simulated using the stationary and non-stationary models,and micro-fractures are monitored. A comparison of these results show that the stationary model can only be used to describe the decay and steady stages,while the non- stationary model is also applicable to the accelerated stage,and thus can simulate the whole process of rock creep.Accelerated creep results from the accelerated development of mirco-fractures and leads to shear failure. KEY WORDS rock mechanics:creep:constitutive models;particle flow code;fractures 目前，对岩石流变特性的研究主要集中在室内试模型，所采用的流变参数大多为常量，无法描述岩石蠕 验·习、流变模型等方面.流变模型通常用三种理 变全过程.因此，许多学者开展了非线性和非定常流 想元件（弹性元件、黏性元件和塑性元件）组合而成， 变本构模型的研究，主要集中在两个方面：一是构造非 由于上述元件均为线性元件，组合出的模型均为线性 线性元件：二是将线性模型的流变参数视为非定常值. 收稿日期：2014-10-28 基金项目：国家自然科学基金资助项目(51074014)：中央直属高校基本科研业务费资助项目(FRF-SD-12O02A)

工程科学学报，第 37 卷，第 7 期: 831--838，2015 年 7 月 Chinese Journal of Engineering，Vol． 37，No． 7: 831--838，July 2015 DOI: 10． 13374 /j． issn2095--9389． 2015． 07． 002; http: / /journals． ustb． edu． cn 基于颗粒流程序的非定常西原体模型 杨振伟，金爱兵，高永涛，王 凯，孙 浩 北京科技大学土木与环境工程学院，北京 100083  通信作者，E-mail: jinaibing@ vip． sina． com 摘 要 为实现岩石试样蠕变全过程的准确模拟，并从细观角度探究蠕变过程中微裂隙的发生和发展规律，在二维颗粒流程 序( PFC2D ) 中开发出具有黏弹塑性特征的西原体流变接触本构模型，进一步提出包含两种非定常元件的非定常西原体模型， 推导了模型本构关系和蠕变方程． 在 PFC2D中调用自定义西原体流变模型，通过参数调试，获得与真实试样具有相同强度特 性的数值试样． 以室内单轴压缩蠕变试验数据为基础，在 Matlab 中对模型非定常参数进行拟合反演分析． 在此基础上，进行 单轴压缩蠕变试验的模拟，计算过程中分别采用定常和非定常两种模型，并对微裂隙进行监测． 对比分析结果表明: 定常模 型仅适用于衰减和稳定蠕变阶段; 非定常模型也可用于描述加速蠕变阶段，从而准确模拟蠕变全过程; 加速蠕变阶段主要是 由微裂隙的加速发展而产生，加速蠕变将导致试样剪切破坏． 关键词 岩石力学; 蠕变; 本构模型; 颗粒流程序; 裂隙 分类号 TU452 Non-stationary Nishihara model in the particle flow code YANG Zhen-wei，JIN Ai-bing ，GAO Yong-tao，WANG Kai，SUN Hao School of Civil and Environmental Engineering，University of Science and Technology Beijing，Beijing 100083，China  Corresponding author，E-mail: jinaibing@ vip． sina． com ABSTＲACT In order to accurately simulate the whole process of rock creep and explore the generation and development of mirco￾fractures from the meso-level，a visco-elastoplastic Nishihara rheological constitutive model is developed in the two-dimensional particle flow code ( PFC2D ) ． Then a non-stationary Nishihara model including two non-stationary elements is put forward in this article． A con￾stitutive equation and a creep equation are derived based on the non-stationary Nishihara model． The user-defined Nishihara constitu￾tive model is called in PFC2D，and a numerical sample whose strength properties are the same as those of a real rock specimen is ac￾quired by parameter testing． Using the data of uniaxial creep tests in laboratory，the non-stationary parameters are back analyzed in Matlab． At last，uniaxial creep tests are simulated using the stationary and non-stationary models，and micro-fractures are monitored． A comparison of these results show that the stationary model can only be used to describe the decay and steady stages，while the non￾stationary model is also applicable to the accelerated stage，and thus can simulate the whole process of rock creep． Accelerated creep results from the accelerated development of mirco-fractures and leads to shear failure． KEY WOＲDS rock mechanics; creep; constitutive models; particle flow code; fractures 收稿日期: 2014--10--28 基金项目: 国家自然科学基金资助项目( 51074014) ; 中央直属高校基本科研业务费资助项目( FＲF--SD--12--002A) 目前，对岩石流变特性的研究主要集中在室内试 验［1 － 3］、流变模型［4 － 6］等方面． 流变模型通常用三种理 想元件( 弹性元件、黏性元件和塑性元件) 组合而成， 由于上述元件均为线性元件，组合出的模型均为线性 模型，所采用的流变参数大多为常量，无法描述岩石蠕 变全过程． 因此，许多学者开展了非线性和非定常流 变本构模型的研究，主要集中在两个方面: 一是构造非 线性元件; 二是将线性模型的流变参数视为非定常值．

·832· 工程科学学报，第37卷，第7期 新型非线性元件的构造有多种方式：Zhang等切 1 定常西原体模型 通过将非线性弹性元件引入Burgers模型，得到了用于 描述水电站坝底破碎区域岩石流变特性的改进B- 1.1连续介质理论西原体模型 ges模型；：宋勇军等圆提出含分数阶微积分的软体元 西原体模型是一个五元件流变模型，由Bingham 件，并组合出四元件流变模型，可用于描述低应力水平 体和Kelvin体串联而成，其结构如图1所示.o,为塑 下初始蠕变和稳定蠕变，以及高应力水平下加速蠕变； 性强度，决定模型是否进入塑性阶段 徐卫亚等网在非线性黏塑性体基础上构建出可反映岩 石加速蠕变特性的河海模型，并通过蠕变试验进行了 验证：蒋昱州等@和赵延林等四也对非线性流变模 型进行了研究 此外，在线性模型中引入非定常参数的研究也有 Bingham Kelvin 很多：Yang和Cheng在黏弹性剪切蠕变模型的基础 图1西原体流变模型 上，提出了新的非定常流变模型，能够与试验数据更好 Fig.1 Nishihara theology model 地吻合：Zhao等图通过损伤变量和硬化函数的引入建 立了适用于软岩的非线性流变模型：潘晓明等提出 根据各元件流变特性，可推导出西原体模型本构 方程、蠕变方程、卸载方程和松弛方程四如下 非定常西原体流变模型，并在ABAQUS中开发出计算 (1)当o<g时， 程序，用于岩石边坡流变数值分析：阎岩等的在西原 2 体模型的基础上，研究模型参数随时间与应力的非线 E,+E0= E2E1+ 0+ E,+Ee+E,+E⊙， (1) 性变化规律，得到了变参数蠕变方程，并用于蠕变试验 研究：吕爱钟等a、余成学等切、康永刚等图均在线 s(0=+(1-e), (2) 性模型的基础上作了非定常参数分析. s0爱e-e (3) 以上对非线性流变模型的研究均将岩石视为连续 介质，而岩石是典型的非连续、非均质材料，使用连续 EE2 (4) 介质理论进行的模拟计算，无法监测试样内部微裂隙 0=2g8,1-。9)+Ee 的发生和发展过程，无法从微观结构研究试样在蠕变 式中，o为应力，e为应变，E,和E2为弹性系数，n和2 过程中的破坏.由Cundall和Strack网提出的颗粒流 为黏性系数，1'为卸载时间，σ。为恒定应力，8。为恒定 法(particle flow code,PFC)及PFC程序是求解非连续 应变 介质力学问题的一个重要数值分析方法，该程序以互 (2)当σ≥o,时， 相黏结的颗粒集合模拟试样，通过颗粒间黏结作用的 断裂和发展模拟加速蠕变阶段试样中微裂隙的发生和 -+(层+”+22=m+” 发展.目前，在这方面进行岩石试样流变试验的研究 (5) 却为数不多.王涛等圆开发了适用于P℉C2”的广义 0管++g-e， 71 E (6) Kelvin接触模型，研究了圆形巷道周围裂隙的发展规 律并以实际算例验证了该模型在流变模拟中的适 s)=-0r+(e-1e, (7) 71 用性. 本文在二维颗粒流程序中，开发出西原体流变 o() 本构模型，提出非定常弹性体和非定常黏性体，并构 建出非定常西原体模型.非定常弹性体的取值随时 后） (8) 间呈指数型增长，可反映衰减蠕变阶段瞬时变形模 式中， 量增加的“硬化”过程：非定常黏性体的取值随时间 P+H P-H 与应力呈幂律型减小，可描述加速蠕变阶段蠕变速 a=20B=20 率增加的“软化”过程.经过PFC2n反复调试和Mal- lab反演分析，确定了计算参数，通过FISH语言编写 P-2++公02=m-40 计算程序，使得西原体模型参数可随时间变化，实现 1.2P℉C中西原体本构模型开发 本构模型的非定常计算，从而有效地模拟岩石蠕变 由于西原体模型塑性元件的存在，在每步计算时 全过程. 需要判断两实体之间的接触是否进入塑性阶段，颗粒

工程科学学报，第 37 卷，第 7 期 新型非线性元件的构造有多种方式: Zhang 等［7］ 通过将非线性弹性元件引入 Burgers 模型，得到了用于 描述水电站坝底破碎区域岩石流变特性的改进 Bur￾gers 模型; 宋勇军等［8］提出含分数阶微积分的软体元 件，并组合出四元件流变模型，可用于描述低应力水平 下初始蠕变和稳定蠕变，以及高应力水平下加速蠕变; 徐卫亚等［9］在非线性黏塑性体基础上构建出可反映岩 石加速蠕变特性的河海模型，并通过蠕变试验进行了 验证; 蒋昱州等［10］和赵延林等［11］ 也对非线性流变模 型进行了研究． 此外，在线性模型中引入非定常参数的研究也有 很多: Yang 和 Cheng［12］在黏弹性剪切蠕变模型的基础 上，提出了新的非定常流变模型，能够与试验数据更好 地吻合; Zhao 等［13］通过损伤变量和硬化函数的引入建 立了适用于软岩的非线性流变模型; 潘晓明等［14］提出 非定常西原体流变模型，并在 ABAQUS 中开发出计算 程序，用于岩石边坡流变数值分析; 阎岩等［15］在西原 体模型的基础上，研究模型参数随时间与应力的非线 性变化规律，得到了变参数蠕变方程，并用于蠕变试验 研究; 吕爱钟等［16］、佘成学等［17］、康永刚等［18］均在线 性模型的基础上作了非定常参数分析． 以上对非线性流变模型的研究均将岩石视为连续 介质，而岩石是典型的非连续、非均质材料，使用连续 介质理论进行的模拟计算，无法监测试样内部微裂隙 的发生和发展过程，无法从微观结构研究试样在蠕变 过程中的破坏． 由 Cundall 和 Strack［19］提出的颗粒流 法( particle flow code，PFC) 及 PFC 程序是求解非连续 介质力学问题的一个重要数值分析方法，该程序以互 相黏结的颗粒集合模拟试样，通过颗粒间黏结作用的 断裂和发展模拟加速蠕变阶段试样中微裂隙的发生和 发展． 目前，在这方面进行岩石试样流变试验的研究 却为数不多． 王涛等［20］ 开发了适用于 PFC2D 的广义 Kelvin 接触模型，研究了圆形巷道周围裂隙的发展规 律并以实际算例验证了该模型在流变模拟中的适 用性． 本文在二维颗粒流程序中，开发出 西 原 体 流 变 本构模型，提出非定常弹性体和非定常黏性体，并构 建出非定常西原体模型． 非定常弹性体的取值随时 间呈指数型增长，可反映衰减蠕变阶段瞬时变形模 量增加的“硬化”过程; 非定常黏性体的取值随时间 与应力呈幂律型减小，可描述加速蠕变阶段蠕变速 率增加的“软化”过程． 经过 PFC2D反复调试和 Mat￾lab 反演分析，确定了计算参数，通过 FISH 语言编写 计算程序，使得西原体模型参数可随时间变化，实现 本构模型的非定常计算，从而有效地模拟岩石蠕变 全过程． 1 定常西原体模型 1. 1 连续介质理论西原体模型 西原体模型是一个五元件流变模型，由 Bingham 体和 Kelvin 体串联而成，其结构如图 1 所示． σs为塑 性强度，决定模型是否进入塑性阶段． 图 1 西原体流变模型 Fig． 1 Nishihara rheology model 根据各元件流变特性，可推导出西原体模型本构 方程、蠕变方程、卸载方程和松弛方程［21 － 22］如下． ( 1) 当 σ ＜ σs时， σ + η2 E2 + E1 σ · = E2E1 E2 + E1 ε + E1η2 E2 + E1 ε ·， ( 1) ε( t) = σ0 E1 + σ0 E2 ( 1 － e － E2 η2 t ) ， ( 2) ε( t) = σ0 E2 ( e E2 η2 t' － 1) e － E2 η2 t ， ( 3) σ( t) = E1E2 E2 + E1 ε0 ( 1 － e － E2 + E1 η2 t ) + E1ε0 e － E2 + E1 η2 t ． ( 4) 式中，σ 为应力，ε 为应变，E1和 E2为弹性系数，η1和 η2 为黏性系数，t'为卸载时间，σ0 为恒定应力，ε0 为恒定 应变． ( 2) 当 σ≥σs时， σ － σs ( + η1 E1 + η2 + η1 E ) 2 σ · + η1η2 E1E2 σ ·· = η1ε · + η1η2 E2 ε ·· ， ( 5) ε( t) = σ0 E1 + σ0 － σs η1 t + σ0 E2 ( 1 － e － E2 η2 t ) ， ( 6) ε( t) = σ0 － σs η1 t' + σ0 E2 ( e E2 η2 t' － 1) e － E2 η2 t ， ( 7) σ( t) = E1 α － [ ( β － E2 η1 + α ) e － αt ( + E2 η1 － β ) e － β ]t ε0 + σs ( H e － αt α － e － βt ) β + σs． ( 8) 式中， α = P + H 2Q ，β = P － H 2Q ， P = η1 E1 + η2 E2 + η1 E2 ，Q = η1η2 E1E2 ，H = P2 槡 － 4Q． 1. 2 PFC 中西原体本构模型开发 由于西原体模型塑性元件的存在，在每步计算时 需要判断两实体之间的接触是否进入塑性阶段，颗粒 · 238 ·

杨振伟等：基于颗粒流程序的非定常西原体模型 ·833· 间的塑性阈值记为∫， (20) (l)当∫<f时数值积分方法.对于Kelvin体部 分，可以推得： 式中， -E2u2±f E2△ 42=- (9) A=1+2m 8=1器 式中，，为Kelvin体相对位移，f为颗粒间接触力，± 对于Bingham体，可以推得 分别表示法向和切向. ∫- 利用有限差分法，取2和∫的平均值可以得到 =±官小 (21) *2 利用有限差分法，取∫的平均值，可以得到 .(10) △21 2 (22) 式中，△：为时间步长，上标t为当前时步计算值，上标 △ 201 t+1为下一时步计算值 对上式进行整理，可以得到山，的计算公式： 对上式进行整理，可以得到下一时步山2的计算 =±±+24+.23) 公式： E 2m1 71 对于整个西原体模型： ={±六"+小 (11) u=u2+1, (24) 式中， h1-n=u1-u;+- (25) E,△l A=1+2 E2△1 ,B=1-2m2 颗粒之间的接触力∫为 对于Bingham体部分，可以推得 "=±{-i+(-月)±] 71 =±官 (26) (12) 式中， 式中，山，为Bingham体相对位移. G总官+六0，扇店+瑞 △11，△ 利用有限差分法，取f的平均值，可以得到 (13) 因此，获得颗粒之间下一时步的接触力∫： △1 f1= 对上式进行整理，可以得到下一时步山，的计算公式： :±+听 (14) ±2[-t+(1-)0f] f<f), E 对于整个西原体模型： +±[-d+1-只）0±n] (f≥f). u=L2+u1, (15) u1-u=-+*1- (16) (27) 颗粒之间的接触力∫为 将上式按照P℉C中本构模型开发流程进行编程， 在Visual Studio2005中进行编译，将所开发的模型存 f=±[-+(-)吗年D小 储到动态链接库文件(DLL)中，供P℉C程序调用. (17) 2非定常西原体模型 式中， G点+官0，瑞营 △1 典型的岩石蠕变曲线分为三个阶段，如图2所示， AB段为衰减蠕变阶段，BC为稳定蠕变阶段，CD为加 (2)当≥f时数值积分方法.对于Kelvin体，可 速蠕变阶段.三个阶段应变率随时间变化曲线如图3 以推得 所示. 根据范庆忠网的研究，岩石蠕变过程同时存在 (18) “硬化”和“软化”两种现象.在衰减蠕变阶段，应变率 利用有限差分法，取山2和∫的平均值可以得到 迅速减小，岩石瞬时变形模量呈增大趋势，即为“硬 "-4川5心2] 化”；进入稳定蠕变阶段后，应变率和瞬时变形模量均 △t2 2 保持不变；直至加速蠕变阶段，应变率逐渐增加，瞬时 (19) 变形模量相应减小，即为“软化”.由简单元件组成的 对上式进行整理，可以得到2的计算公式： 线性模型无法模拟蠕变曲线的三个阶段

杨振伟等: 基于颗粒流程序的非定常西原体模型 间的塑性阈值记为 fs． ( 1) 当 f ＜ fs时数值积分方法． 对于 Kelvin 体部 分，可以推得: u · 2 = － E2 u2 ± f η2 ． ( 9) 式中，u2 为 Kelvin 体相对位移，f 为颗粒间接触力，± 分别表示法向和切向． 利用有限差分法，取 u2 和 f 的平均值可以得到 ut + 1 2 － ut 2 Δt = 1 η [ 2 － E2 ( ut + 1 2 + ut 2 ) 2 ± f t + 1 + f t ] 2 ． ( 10) 式中，Δt 为时间步长，上标 t 为当前时步计算值，上标 t + 1 为下一时步计算值． 对上式进行整理，可以得到下一时步 u2 的计算 公式: ut + 1 2 = 1 [ A But 2 ± Δt 2η2 ( f t + 1 + f t ] ) ． ( 11) 式中， A = 1 + E2Δt 2η2 ，B = 1 － E2Δt 2η2 ． 对于 Bingham 体部分，可以推得 u · 1 = ± f · E1 ． ( 12) 式中，u1 为 Bingham 体相对位移． 利用有限差分法，取 f 的平均值，可以得到 ut + 1 1 － ut 1 Δt = ± f t + 1 － f t E1Δt ． ( 13) 对上式进行整理，可以得到下一时步 u1 的计算公式: ut + 1 1 = ± f t + 1 － f t E1 + ut 1 ． ( 14) 对于整个西原体模型: u = u2 + u1， ( 15) ut + 1 － ut = ut + 1 2 － ut 2 + ut + 1 1 － ut 1 ． ( 16) 颗粒之间的接触力 f t + 1为 f t + 1 = ± 1 C [ 1 ut + 1 － ut ( + 1 － B ) A ut 2D1 f ]t ． ( 17) 式中， C1 = Δt 2η2A + 1 E1 ，D1 = Δt 2η2A － 1 E1 ． ( 2) 当 f≥ fs时数值积分方法． 对于 Kelvin 体，可 以推得 u · 2 = － E2 u2 ± f η2 ． ( 18) 利用有限差分法，取 u2 和 f 的平均值可以得到 ut + 1 2 － ut 2 Δt = 1 η [ 2 － E2 ( ut + 1 2 + ut 2 ) 2 ± f t + 1 + f t ] 2 ． ( 19) 对上式进行整理，可以得到 u2 的计算公式: ut + 1 2 = 1 [ A But 2 ± Δt 2η2 ( f t + 1 + f t ] ) ． ( 20) 式中， A = 1 + E2Δt 2η2 ，B = 1 － E2Δt 2η2 ． 对于 Bingham 体，可以推得 u · 1 = ± f · E1 ± f － fs η1 ． ( 21) 利用有限差分法，取 f 的平均值，可以得到 ut + 1 1 － ut 1 Δt = ± f t + 1 － f t E1Δt ± f t + 1 + f t 2η1  fs η1 ． ( 22) 对上式进行整理，可以得到 u1 的计算公式: ut + 1 1 = ± f t + 1 － f t E1 ± Δt( f t + 1 + f t ) 2η1  Δt·fs η1 + ut 1 ． ( 23) 对于整个西原体模型: u = u2 + u1， ( 24) ut + 1 － ut = ut + 1 2 － ut 2 + ut + 1 1 － ut 1 ． ( 25) 颗粒之间的接触力 f t + 1为 f t + 1 = ± 1 C [ 2 ut + 1 － ut ( + 1 － B ) A ut 2D2 f t ± Δt·fs η ] 1 ( 26) 式中， C2 = Δt 2η2A + 1 E1 + Δt 2η1 ，D2 = Δt 2η2A － 1 E1 + Δt 2η1 ． 因此，获得颗粒之间下一时步的接触力 f t + 1 : f t + 1 = ± 1 C [ 1 ut + 1 － ut + 1 － ( B ) A ut 2D1 f ]t ( f ＜ fs) ， ± 1 C [ 2 ut + 1 － ut + 1 － ( B ) A ut 2D2 f t ± σsΔt η ] 1 ( f≥fs { ) ． ( 27) 将上式按照 PFC 中本构模型开发流程进行编程， 在 Visual Studio 2005 中进行编译，将所开发的模型存 储到动态链接库文件( DLL) 中，供 PFC 程序调用． 2 非定常西原体模型 典型的岩石蠕变曲线分为三个阶段，如图 2 所示， AB 段为衰减蠕变阶段，BC 为稳定蠕变阶段，CD 为加 速蠕变阶段． 三个阶段应变率随时间变化曲线如图 3 所示． 根据范庆忠［23］ 的研究，岩石蠕变过程同时存在 “硬化”和“软化”两种现象． 在衰减蠕变阶段，应变率 迅速减小，岩石瞬时变形模量呈增大趋势，即为“硬 化”; 进入稳定蠕变阶段后，应变率和瞬时变形模量均 保持不变; 直至加速蠕变阶段，应变率逐渐增加，瞬时 变形模量相应减小，即为“软化”． 由简单元件组成的 线性模型无法模拟蠕变曲线的三个阶段． · 338 ·

·834· 工程科学学报，第37卷，第7期 本文提出参数值与时间呈指数关系的非定常弹性 D 体，其参数取值与时间的关系如下： E()=E。[a+b×e-w]. (29) 式中：E。为初始弹性系数：a、b和c均为量纲一的常 数，可由室内试验结果拟合得到；t为蠕变时间：。为时 间参数(t和。均取为量纲一的常数) 2.2非定常黏性体 对于黏性系数，孙钧网研究指出，在非线性黏弹 性问题中，应变率符合一种幂律型的关系： E=Aσr"t". (30) 式中，为蠕变应变率，o为等效应力，A、n和m为拟 图2岩石蠕变全过程曲线 合蠕变参数 Fig.2 Whole process of rock creep 熊良宵等网研究认为非线性黏滞体的黏滞系数 可以表达为 o0=(侣)广() (31) 式中，。为初始黏滞系数，σ1为参考应力，1为参考 安 时间. 本文提出参数值与时间、应力存在指数关系的非 定常黏性体，其参数取值与时间、应力的关系如下： B n(σ，l)=- (32) *侣 式中：o为初始黏性系数pqn和m均为量纲一的常 图3岩石蠕变过程应变率变化 数，可由室内试验结果拟合得到：σ1为应力参数，为 Fig.3 Strain rate of rock creep 时间参数（σ、g,、1和t,取为量纲一的常数）. 然而，在P℉℃中根据定常西原体本构方程所开发 将本文提出的非定常弹性体和非定常黏性体替换 的定常线性模型，参数不随时间变化，只能体现前两个 部分定常元件，组合成非定常西原体模型见图4. 蠕变阶段，对于加速蠕变阶段无法体现.大量研究证 E 明岩石的流变参数随时间和应力变化非常明 显48242.因此，岩石蠕变过程中流变参数应该为 非定常值. E0) 70) 2.1非定常弹性体 图4非定常西原体流变模型 鉴于在衰减蠕变阶段瞬时变形模量有随时间增加 Fig.4 Non-stationary Nishihara rheology model 的趋势，丁志坤和吕爱钟网认为弹性系数与时间存在 指数关系，并将H一K体弹性系数的变化趋势总结为 根据熊良宵等四对非定常Burgers模型的蠕变方 Ex(t）=p1+P2e. (28) 程的推导方法，推导出非定常西原体模型的本构方程 式中，P1P2和p3为参数，由试验拟合得到. 如下： 72 E,E(） 【o+E,+EG)°=E,+EOe+E,+Ee (o<), o-a+[20,+2a0]b+2a=o08+o6 (33) E(t)E2 (o≥o）. E2 E2 非定常西原体模型的蠕变方程如下： 0+管u-e (g<o), e(t)= 。(o-0)+(1-e) (34) E0+n(c,)+2 E (σ≥o)

工程科学学报，第 37 卷，第 7 期 图 2 岩石蠕变全过程曲线 Fig． 2 Whole process of rock creep 图 3 岩石蠕变过程应变率变化 Fig． 3 Strain rate of rock creep 然而，在 PFC 中根据定常西原体本构方程所开发 的定常线性模型，参数不随时间变化，只能体现前两个 蠕变阶段，对于加速蠕变阶段无法体现． 大量研究证 明岩 石 的 流 变 参 数 随 时 间 和 应 力 变 化 非 常 明 显［14 － 18，24 － 28］． 因此，岩石蠕变过程中流变参数应该为 非定常值． 2. 1 非定常弹性体 鉴于在衰减蠕变阶段瞬时变形模量有随时间增加 的趋势，丁志坤和吕爱钟［27］认为弹性系数与时间存在 指数关系，并将 H--K 体弹性系数的变化趋势总结为 EK ( t) = p1 + p2 ep3t ． ( 28) 式中，p1、p2和 p3为参数，由试验拟合得到． 本文提出参数值与时间呈指数关系的非定常弹性 体，其参数取值与时间的关系如下: E( t) = E0［a + b × ec( t － t0) ］． ( 29) 式中: E0 为初始弹性系数; a、b 和 c 均为量纲一的常 数，可由室内试验结果拟合得到; t 为蠕变时间; t0为时 间参数( t 和 t0均取为量纲一的常数) ． 2. 2 非定常黏性体 对于黏性系数，孙钧［28］研究指出，在非线性黏弹 性问题中，应变率符合一种幂律型的关系: ε · = Aσn t m ． ( 30) 式中，ε · 为蠕变应变率，σ 为等效应力，A、n 和 m 为拟 合蠕变参数． 熊良宵等［29］研究认为非线性黏滞体的黏滞系数 可以表达为 η( σ，t) = η0 ( σ σ ) 1 ( n t t ) 1 m ． ( 31) 式中，η0 为初 始 黏 滞 系 数，σ1 为参 考 应 力，t1 为参 考 时间． 本文提出参数值与时间、应力存在指数关系的非 定常黏性体，其参数取值与时间、应力的关系如下: η( σ，t) = η0 p + ( q σ σ ) 1 ( n t t ) 1 m ． ( 32) 式中: η0为初始黏性系数; p、q、n 和 m 均为量纲一的常 数，可由室内试验结果拟合得到; σ1为应力参数，t1 为 时间参数( σ、σ1、t 和 t1取为量纲一的常数) ． 将本文提出的非定常弹性体和非定常黏性体替换 部分定常元件，组合成非定常西原体模型见图 4． 图 4 非定常西原体流变模型 Fig． 4 Non-stationary Nishihara rheology model 根据熊良宵等［29］对非定常 Burgers 模型的蠕变方 程的推导方法，推导出非定常西原体模型的本构方程 如下: σ + η2 E2 + E( t) σ · = E2E( t) E2 + E( t) ε + E( t) η2 E2 + E( t) ε · ( σ ＜ σs) ， σ － σs [ + η( σ，t) E( t) + η2 + η( σ，t) E ] 2 σ · + η( σ，t) η2 E( t) E2 σ ·· = η( σ，t) ε · + η( σ，t) η2 E2 ε ·· ( σ≥σs { ) ． ( 33) 非定常西原体模型的蠕变方程如下: ε( t) = σ0 E( t) + σ0 E2 ( 1 － e － E2 η2 t ) ( σ ＜ σs) ， σ0 E( t) + ( σ0 － σs) η( σ，t) t + σ0 E2 ( 1 － e － E2 η2 t ) ( σ≥σs { ) ． ( 34) · 438 ·

杨振伟等：基于颗粒流程序的非定常西原体模型 835· 即 o E。(a+be-w)+ (1-e) E (o<o）, e(t)= ,-,)[p+9()(片)门 (35) LEo (a+be-) 3实例验证 杨彩红等四对天然页岩进行了单轴压缩蠕变试 验，本文选取其中轴压为55MPa的试验数据用于非定 常西原体模型研究，其研究过程如下 (1)模型建立.为从微细观角度研究试样在蠕变 过程中微裂隙的发展过程，以及加速蠕变过程的发展 机理，本文在PFC2D中通过Fishtank生成50mm×lO0 mm的数值试样，见图5. (2)定常参数调试.调用西原体本构模型DLL文 图5颗粒流数值试样 件，并将所有接触替换为西原体模型，通过单轴和三轴 Fig.5 Numerical sample in the particle flow code 压缩试验对参数进行反复调试，最终获得与天然页岩 具有相同强度的数值试样，并用于蠕变试验研究.所 采用的细观力学参数见表1. 表1数值试样细观力学参数 Table 1 Meso-mechanical parameters of the numerical sample 最小颗粒半径， 颗粒体密度， 粒间摩擦系数4 颗粒弹性模量： 颗粒法向/切向刚 最大/最小粒径比 Ri/mm el(kg'm-3) E/GPa 度比，kk 0.4 1.66 2500 0.5 18.0 3.0 平行黏结法向强度 平行黏结切向强度 平行黏结半径 平行黏结弹性 平行黏结法向/切 系数，A 模量，E。/GPa 向刚度比，k。k 平均值， 标准差， 平均值， 标准差， Comem /MPa a-e /MPa T/MPa Te /MPa 1.0 18.0 3.0 80.0 ±25.0 80.0 ±25.0 Eo/GPa no /GPa E2/(GPas) 2/(GPa-s) o./MPa 摩擦系数 11 800 10 1000 50 0.3 (3)定常模型计算.采用FISH语言编写蠕变试 1.0r 验命令流，删去侧墙约束，利用伺服机制控制沿y轴加 0.9 ,定常拟合 载应力55MPa.轴向应变将随时间推移而增加，从而 08 ★实验数据 ★ 实现单轴蠕变试验的模拟. 0.7 利用定常模型进行蠕变试验模拟时，E。oE~ 0.6 ★ ,0,摩擦系数等参数均为定常参数，见表1.模拟计 0.5 55 MPa 算结果与室内试验的拟合效果见图6. 0.4 由图6可知，定常西原体模型对衰减蠕变阶段和 0.3 稳定蠕变阶段拟合效果良好，但是不能模拟加速蠕变 0.2* 阶段 0.1 50100150200250300350 (4)非定常参数反演.根据式(35)中蠕变方程的 时间h 图6定常西原体模型拟合曲线 具体形式，在Matlab的拟合工具Cftool中，对55MPa Fig.6 Fitting curves by the stationary Nishihara model 的试验数据进行拟合分析，反演得到的各参数如表2 所示，相关性系数为0.953，表明该模型可用于岩石蠕 中，非定常参数E()=E。(a+be-W),n(σ，t）=no/ 变全过程的模拟. p+g(σ/o,)“(t,)]均与时间有关，在蠕变试验命 (5)非定常模型计算.非定常西原体模型参数有 令流中需设置每个计算步时更新E(t)和η（σ，t)的取 六个，分别是E(t）、n(σ，)、E2、2、o,和摩擦系数.其 值：E。no、E2n2σ，、摩擦系数等参数由PFC数值试验

杨振伟等: 基于颗粒流程序的非定常西原体模型 即 ε( t) = σ0 E0 ( a + bec( t － t0) ) + σ0 E2 ( 1 － e － E2 η2 t ) ( σ ＜ σs) ， σ0 E0 ( a + bec( t － t0) ) + t( σ0 － σs [ ) p + ( q σ0 σ ) 1 ( n t t ) 1 ] m η0 + σ0 E2 ( 1 － e － E2 η2 t ) ( σ≥σs)        ． ( 35) 3 实例验证 杨彩红等［30］对天然页岩进行了单轴压缩蠕变试 验，本文选取其中轴压为 55 MPa 的试验数据用于非定 常西原体模型研究，其研究过程如下． ( 1) 模型建立． 为从微细观角度研究试样在蠕变 过程中微裂隙的发展过程，以及加速蠕变过程的发展 机理，本文在 PFC2D中通过 Fishtank 生成 50 mm × 100 mm 的数值试样，见图 5． ( 2) 定常参数调试． 调用西原体本构模型 DLL 文 件，并将所有接触替换为西原体模型，通过单轴和三轴 压缩试验对参数进行反复调试，最终获得与天然页岩 具有相同强度的数值试样，并用于蠕变试验研究． 所 图 5 颗粒流数值试样 Fig． 5 Numerical sample in the particle flow code 采用的细观力学参数见表 1． 表 1 数值试样细观力学参数 Table 1 Meso-mechanical parameters of the numerical sample 最小颗粒半径， Ｒmin /mm 最大/最小粒径比 颗粒体密度， ρ /( kg·m － 3 ) 粒间摩擦系数，μ 颗粒弹性模量， Ec /GPa 颗粒法向/切向刚 度比，kn /ks 0. 4 1. 66 2500 0. 5 18. 0 3. 0 平行黏结半径 系数，λ 平行黏结弹性 模量，Ec /GPa 平行黏结法向/切 向刚度比，kn /ks 平行黏结法向强度 平行黏结切向强度 平均值， σn-mean /MPa 标准差， σn-dev /MPa 平均值， τs-mean /MPa 标准差， τs-dev /MPa 1. 0 18. 0 3. 0 80. 0 ± 25. 0 80. 0 ± 25. 0 E0 /GPa η0 /GPa E2 /( GPa·s) η2 /( GPa·s) σs /MPa 摩擦系数 11 800 10 1000 50 0. 3 ( 3) 定常模型计算． 采用 FISH 语言编写蠕变试 验命令流，删去侧墙约束，利用伺服机制控制沿 y 轴加 载应力 55 MPa． 轴向应变将随时间推移而增加，从而 实现单轴蠕变试验的模拟． 利用定常模型进行蠕变试验模拟时，E0、η0、E1、 η1、σs、摩擦系数等参数均为定常参数，见表 1． 模拟计 算结果与室内试验的拟合效果见图 6． 由图 6 可知，定常西原体模型对衰减蠕变阶段和 稳定蠕变阶段拟合效果良好，但是不能模拟加速蠕变 阶段． ( 4) 非定常参数反演． 根据式( 35) 中蠕变方程的 具体形式，在 Matlab 的拟合工具 Cftool 中，对 55 MPa 的试验数据进行拟合分析，反演得到的各参数如表 2 所示，相关性系数为 0. 953，表明该模型可用于岩石蠕 变全过程的模拟． ( 5) 非定常模型计算． 非定常西原体模型参数有 六个，分别是 E( t) 、η( σ，t) 、E2、η2、σs和摩擦系数． 其 图 6 定常西原体模型拟合曲线 Fig． 6 Fitting curves by the stationary Nishihara model 中，非定常参数 E( t) = E0 ( a + bec( t － t0) ) ，η( σ，t) = η0 / ［p + q( σ/σ1 ) n ( t /t1 ) m ］均与时间有关，在蠕变试验命 令流中需设置每个计算步时更新 E( t) 和 η( σ，t) 的取 值; E0、η0、E2、η2、σs、摩擦系数等参数由 PFC 数值试验 · 538 ·

·836· 工程科学学报，第37卷，第7期 调试得出，见表1：其余非定常参数由室内试验数据反 演得到，见表2. 表2非定常西原体参数取值 Table 2 Parameters of the non-stationary Nishihara model a b C to 9 m 0.017 8.7×10-10 -0.099 98.3 103.7 3.507 1.36 53.75 12.35 213.4 采用非定常西原体模型与试验结果进行拟合，效 为集中；到280h,基本出现与水平方向呈60°左右的 果见图7.由图7可知，非定常西原体模型不仅能有效 “X型”剪切破坏 地模拟岩石蠕变的衰减阶段和稳定阶段，而且对加速 蠕变阶段也能有很好的体现，模拟结果与试验数据误 差较小. 1.0r 非定常拟合 0.9F ★实验数据 0.8 07 0.6 100h 150h 0.5 55 MPa 0.4 0.3 02 0.1 50 100150200250300350 时间h 图7非定常西原体模型拟合曲线 200h 250h 280h Fig.7 Fitting curves by the non-stationary Nishihara model 图9镯变试验中微裂隙发展过程 Fig.9 Development of micro-fractures in creep test (6)不同应力水平验证.在55MPa试验研究基础 上，进行45MPa和65MPa应力水平下蠕变试验验证， 4 结论 所得结果与室内试验结果基本一致，见图8 1.0 定常流变模型只能描述岩石蠕变曲线的前两个阶 0.9 段，对于加速蠕变阶段则无法体现.本文提出非定常 0.8 65 MPa 西原体模型，并在离散单元法计算程序P℉℃中实现模 07 型开发和非定常计算，可用于岩石蠕变全过程的模拟. 0.6 主要研究结论如下， 0.5 55MP: (1)基于颗粒流程序的线性本构模型，通过FISH 0.4 兴4★★★★★★★4平￥￥米 02444*+45MP 语言编程，可实现流变参数的实时更新，从而实现定常 0.3 一非定常拟合 模型的非定常计算. ★实验数据 (2)非定常弹性体的取值随时间呈指数型增长， 01 0 50100150200250300350 可用于反映衰减蠕变阶段瞬时变形模量增加的“硬 时间h 化”过程.非定常黏性体的取值随时间和应力呈幂律 图8不同应力水平下非定常模型拟合曲线 Fig.8 Fitting curves by the non-stationary model under different 型减小，可用于描述加速蠕变阶段蠕变速率增加的 stress levels “软化”过程. (3)分别采用定常与非定常西原体模型对室内试 (7)微裂隙监测与分析.在蠕变试验过程中，通 验数据进行拟合对比分析.结果显示：定常模型仅适 过调用Crk init函数实现对试样微裂隙数量和位置的 用于衰减和稳定蠕变阶段：而非定常模型还可用于描 监测，不同时刻微裂隙的发展见图9.由图9可知：在 述加速蠕变阶段，实现对蠕变全过程的准确模拟. 蠕变试验前期(100h以内)，颗粒间的黏结破坏数量较 (4)通过对蠕变过程中微裂隙的监测，可知蠕变 少，裂纹发展速度缓慢，分布范围广泛，无明显规律：但 过程中加速蠕变阶段主要是由裂隙的加速发展产生， 是随着加载时间的增加，微裂隙加速发展，局部裂隙较 加速蠕变将导致试样被剪切破坏

工程科学学报，第 37 卷，第 7 期 调试得出，见表 1; 其余非定常参数由室内试验数据反 演得到，见表 2． 表 2 非定常西原体参数取值 Table 2 Parameters of the non-stationary Nishihara model a b c t0 p q n σ1 m t1 0. 017 8. 7 × 10 － 10 － 0. 099 98. 3 103. 7 3. 507 1. 36 53. 75 12. 35 213. 4 采用非定常西原体模型与试验结果进行拟合，效 果见图 7． 由图 7 可知，非定常西原体模型不仅能有效 地模拟岩石蠕变的衰减阶段和稳定阶段，而且对加速 蠕变阶段也能有很好的体现，模拟结果与试验数据误 差较小． 图 7 非定常西原体模型拟合曲线 Fig． 7 Fitting curves by the non-stationary Nishihara model ( 6) 不同应力水平验证． 在 55 MPa 试验研究基础 上，进行 45 MPa 和 65 MPa 应力水平下蠕变试验验证， 所得结果与室内试验结果基本一致，见图 8． 图 8 不同应力水平下非定常模型拟合曲线 Fig． 8 Fitting curves by the non-stationary model under different stress levels ( 7) 微裂隙监测与分析． 在蠕变试验过程中，通 过调用 Crk_init 函数实现对试样微裂隙数量和位置的 监测，不同时刻微裂隙的发展见图 9． 由图 9 可知: 在 蠕变试验前期( 100 h 以内) ，颗粒间的黏结破坏数量较 少，裂纹发展速度缓慢，分布范围广泛，无明显规律; 但 是随着加载时间的增加，微裂隙加速发展，局部裂隙较 为集中; 到 280 h，基本出现与水平方向呈 60°左右的 “X 型”剪切破坏． 图 9 蠕变试验中微裂隙发展过程 Fig． 9 Development of micro-fractures in creep test 4 结论 定常流变模型只能描述岩石蠕变曲线的前两个阶 段，对于加速蠕变阶段则无法体现． 本文提出非定常 西原体模型，并在离散单元法计算程序 PFC 中实现模 型开发和非定常计算，可用于岩石蠕变全过程的模拟． 主要研究结论如下． ( 1) 基于颗粒流程序的线性本构模型，通过 FISH 语言编程，可实现流变参数的实时更新，从而实现定常 模型的非定常计算． ( 2) 非定常弹性体的取值随时间呈指数型增长， 可用于反映衰减蠕变阶段瞬时变形模量增加的“硬 化”过程． 非定常黏性体的取值随时间和应力呈幂律 型减小，可用于描述加速蠕变阶段蠕变速率增加的 “软化”过程． ( 3) 分别采用定常与非定常西原体模型对室内试 验数据进行拟合对比分析． 结果显示: 定常模型仅适 用于衰减和稳定蠕变阶段; 而非定常模型还可用于描 述加速蠕变阶段，实现对蠕变全过程的准确模拟． ( 4) 通过对蠕变过程中微裂隙的监测，可知蠕变 过程中加速蠕变阶段主要是由裂隙的加速发展产生， 加速蠕变将导致试样被剪切破坏． · 638 ·

杨振伟等：基于颗粒流程序的非定常西原体模型 837 参考文献 (潘晓明，杨钊，许建聪.非定常西原黏弹塑性流变模型的 应用研究.岩石力学与工程学报，2011,30（增刊1）：2640) [1]Zhang Z L,Xu W Y,Wang W,et al.Triaxial creep tests of rock 15] Yan Y,Wang S J,Wang E Z.Creep equation of variable param- eters based on Nishihara model.Rock Soil Mech,2010,31 from the compressive zone of dam foundation in Xiangjiaba Hydro- (10):3025 power Station.Int J Rock Mech Min Sci,2012,50(1):133 2]Maranini E,Brignoli M.Creep behaviour of a weak rock:experi- (阀岩，王思敬，王恩志.基于西原模型的变参数蠕变方程 岩土力学，2010,31(10)：3025) mental characterization.Int J Rock Mech Min Sci,1999,36(1): 127 [16]Li AZ.DingZ K,Jiao C M,et al.Identification of non-station- B]Li Y S,Xia CC.Time-dependent tests on intact rocks in uniaxial ary creep constitutive models of rock.Chin I Rock Mech Eng, 2008,27(1):16 compression.Int J Rock Mech Min Sci,2000,37(3):467 [4]Zhang H B,Wang Z Y,Zheng Y L,et al.Study on tri-axial (吕爱钟，丁志坤，焦春茂，等.岩石非定常蠕变模型辨识 岩石力学与工程学报，2008,27(1)：16) creep experiment and constitutive relation of different rock salt. Saf Sci,2012,50(4):801 [7]She C X.Research on nonlinear viscoelasto-plastic creep model [5]Liu L,Wang G M,Chen J H,et al.Creep experiment and rheo- of rock.Chin J Rock Mech Eng,2009,28(10):2006 logical model of deep saturated rock.Trans Nonferrous Met Soc (佘成学.岩石非线性黏弹塑性孀变模型研究.岩石力学与 China,2013,23(2):478 工程学报，2009,28(10)：2006) [6]Shao J F,Zhu QZ.Su K.Modeling of creep in rock materials in [18]Kang Y G,Zhang X E.A non-stationary model for rock creep terms of material degradation.Comput Geotech,2003,30(7): based on Burgers model.Rock Soil Mech,2011,32 (Suppl 1): 549 424 7]Zhang Y,Xu W Y,Gu J J,et al.Triaxial creep tests of weak (康永刚，张秀娥.基于Burgers模型的岩石非定常蠕变模 sandstone from fracture zone of high dam foundation.Cent South 型.岩土力学，2011,32（增刊1）：424) Uniw,2013,20(9）:2528 19] Cundall P A,Strack O D L.A discrete numerical model for [8]Song Y J,Lei S Y,Han T L.A new nonlinear viscoelastoplastic granular assemblies.Geotechnique,1979,29(1):47 rheological model for rocks.Rock Soil Mech,2012,33(7):2076 [20]Wang T,Lii O,Li Y,et al.Development of contact model in (宋勇军，雷胜友，韩铁林.一种新的岩石非线性黏弹塑性流 particle discrete element method.Chin J Rock Mech Eng,2009, 变模型.岩土力学，2012,33(7)：2076) 28(Suppl2):4040 Xu WY,Yang SQ,Chu W J.Nonlinear viscoelasto-plastic rheo- (王涛，吕庆，李杨，等.颗粒离散元方法中接触模型的开 logical model (Hohai model)of rock and its engineering applica- 发.岩石力学与工程学报，2009,28（增刊2）：4040) tion.Chin J Rock Mech Eng,2006,25(3):433 21]Tao B,Wu F Q,Guo G M,et al.Flexibility of visco-elastoplast- (徐卫亚，杨圣奇，赭卫江.岩石非线性黏弹塑性流变模型 ic model to rheological characteristics of rock and solution of (河海模型)及其应用.岩石力学与工程学报，2006,25(3)： rheological parameter.Chin Rock Mech Eng,2005,24(17): 433) 3165 [10]Jiang YZ,Zhang MM.Li L Q.Study on nonlinear viscoelasto- (陶波，伍法权，郭改梅，等.西原模型对岩石流变特性的适 plastic creep model of rock and its parameter identification.Chin 应性及其参数确定.岩石力学与工程学报，2005,24(17)： JRock Mech Eng,2008,27(4):832 3165) (蒋昱州，张明鸣，李良权.岩石非线性黏弹塑性蠕变模型 [22] Tian H M,Chen WZ,Zhao W S,et al.Analysis of triaxial 研究及其参数识别.岩石力学与工程学报，2008,27(4)： stress relaxation properties of red silty mudstone of Yichang-Ba- 832) dong highway.Rock Soil Mech,2013,34(4):981 [11]Zhao Y L,Cao P,Wen Y D,et al.Elastovisco-plastic rheologi- (田洪铭，陈卫忠，赵武胜，等.宜一巴高速公路泥质红砂岩 cal experiment and nonlinear rheological model of rocks.Chin 三轴应力松弛特性研究.岩士土力学，2013,34(4)：981) Rock Mech Eng,2008,27(3):477 D3]Fan Q Z.Experimental Study on Creep and its Disturbed Effect of (赵延林，曹平，文有道，等.岩石弹黏塑性流变试验和非线 Rocks [Dissertation].Qingdao:Shandong University of Science 性流变模型研究.岩石力学与工程学报，2008,27(3)： and Technology,2006:30 477) (范庆忠.岩石蠕变及其扰动效应试验研究[学位论文].青 [2]Yang SQ,Cheng L Non-stationary and nonlinear visco-lastic 岛：山东科技大学，2006：30) shear creep model for shale.Int J Rock Mech Min Sci,2011,48 [24 Xu H F.Time depend behaviors of strength and elasticity modu- (6):1011 lus of weak rock.Chin J Rock Mech Eng,1997,16(3):246 [13]Zhao Y L,Cao P,Wang W J,et al.Viscoelastoplastic rheologi- (许宏发.软岩强度和弹模的时间效应研究.岩石力学与工 cal experiment under circular increment step load and unload and 程学报，1997,16(3)：246) nonlinear creep model of soft rocks.J Cent South Unie Technol, 25]Song DZ,Sun J.Nonlinear rheological property of rock materials 2009,16(3):484 and mechanic model.J Tongi Unir,1991,19(4):395 [14]Pan X M,Yang Z,Xu J C.Application study of nonstationary (宋德彰，孙钧。岩质材料非线性流变属性及其力学模型 Nishihara viscoelasto-plastic theological model.Chin Rock 同济大学学报，1991,19(4)：395) Mech Eng,2011,30(Suppl 1):2640 6]Li X Q.The creep calculation method of variable parameters on

杨振伟等: 基于颗粒流程序的非定常西原体模型 参 考 文 献 ［1］ Zhang Z L，Xu W Y，Wang W，et al． Triaxial creep tests of rock from the compressive zone of dam foundation in Xiangjiaba Hydro￾power Station． Int J Ｒock Mech Min Sci，2012，50( 1) : 133 ［2］ Maranini E，Brignoli M． Creep behaviour of a weak rock: experi￾mental characterization． Int J Ｒock Mech Min Sci，1999，36( 1) : 127 ［3］ Li Y S，Xia C C． Time-dependent tests on intact rocks in uniaxial compression． Int J Ｒock Mech Min Sci，2000，37( 3) : 467 ［4］ Zhang H B，Wang Z Y，Zheng Y L，et al． Study on tri-axial creep experiment and constitutive relation of different rock salt． Saf Sci，2012，50( 4) : 801 ［5］ Liu L，Wang G M，Chen J H，et al． Creep experiment and rheo￾logical model of deep saturated rock． Trans Nonferrous Met Soc China，2013，23( 2) : 478 ［6］ Shao J F，Zhu Q Z，Su K． Modeling of creep in rock materials in terms of material degradation． Comput Geotech，2003，30 ( 7 ) : 549 ［7］ Zhang Y，Xu W Y，Gu J J，et al． Triaxial creep tests of weak sandstone from fracture zone of high dam foundation． J Cent South Univ，2013，20( 9) : 2528 ［8］ Song Y J，Lei S Y，Han T L． A new nonlinear viscoelasto-plastic rheological model for rocks． Ｒock Soil Mech，2012，33( 7) : 2076 ( 宋勇军，雷胜友，韩铁林． 一种新的岩石非线性黏弹塑性流 变模型． 岩土力学，2012，33( 7) : 2076) ［9］ Xu W Y，Yang S Q，Chu W J． Nonlinear viscoelasto-plastic rheo￾logical model ( Hohai model) of rock and its engineering applica￾tion． Chin J Ｒock Mech Eng，2006，25( 3) : 433 ( 徐卫亚，杨圣奇，褚卫江． 岩石非线性黏弹塑性流变模型 ( 河海模型) 及其应用． 岩石力学与工程学报，2006，25( 3) : 433) ［10］ Jiang Y Z，Zhang M M，Li L Q． Study on nonlinear viscoelasto￾plastic creep model of rock and its parameter identification． Chin J Ｒock Mech Eng，2008，27( 4) : 832 ( 蒋昱州，张明鸣，李良权． 岩石非线性黏弹塑性蠕变模型 研究及其参数识别． 岩石力学与工程学报，2008，27 ( 4) : 832) ［11］ Zhao Y L，Cao P，Wen Y D，et al． Elastovisco-plastic rheologi￾cal experiment and nonlinear rheological model of rocks． Chin J Ｒock Mech Eng，2008，27( 3) : 477 ( 赵延林，曹平，文有道，等． 岩石弹黏塑性流变试验和非线 性流变模 型 研 究． 岩 石 力 学 与 工 程 学 报，2008，27 ( 3 ) : 477) ［12］ Yang S Q，Cheng L． Non-stationary and nonlinear visco-elastic shear creep model for shale． Int J Ｒock Mech Min Sci，2011，48 ( 6) : 1011 ［13］ Zhao Y L，Cao P，Wang W J，et al． Viscoelasto-plastic rheologi￾cal experiment under circular increment step load and unload and nonlinear creep model of soft rocks． J Cent South Univ Technol， 2009，16( 3) : 484 ［14］ Pan X M，Yang Z，Xu J C． Application study of nonstationary Nishihara viscoelasto-plastic rheological model． Chin J Ｒock Mech Eng，2011，30( Suppl 1) : 2640 ( 潘晓明，杨钊，许建聪． 非定常西原黏弹塑性流变模型的 应用研究． 岩石力学与工程学报，2011，30( 增刊 1) : 2640) ［15］ Yan Y，Wang S J，Wang E Z． Creep equation of variable param￾eters based on Nishihara model． Ｒock Soil Mech，2010，31 ( 10) : 3025 ( 阎岩，王思敬，王恩志． 基于西原模型的变参数蠕变方程． 岩土力学，2010，31( 10) : 3025) ［16］ Lü A Z，Ding Z K，Jiao C M，et al． Identification of non-station￾ary creep constitutive models of rock． Chin J Ｒock Mech Eng， 2008，27( 1) : 16 ( 吕爱钟，丁志坤，焦春茂，等． 岩石非定常蠕变模型辨识． 岩石力学与工程学报，2008，27( 1) : 16) ［17］ She C X． Ｒesearch on nonlinear viscoelasto-plastic creep model of rock． Chin J Ｒock Mech Eng，2009，28( 10) : 2006 ( 佘成学． 岩石非线性黏弹塑性蠕变模型研究． 岩石力学与 工程学报，2009，28( 10) : 2006) ［18］ Kang Y G，Zhang X E． A non-stationary model for rock creep based on Burgers model． Ｒock Soil Mech，2011，32( Suppl 1) : 424 ( 康永刚，张秀娥． 基于 Burgers 模型的岩石非定常蠕变模 型． 岩土力学，2011，32( 增刊 1) : 424) ［19］ Cundall P A，Strack O D L． A discrete numerical model for granular assemblies． Geotechnique，1979，29( 1) : 47 ［20］ Wang T，Lü Q，Li Y，et al． Development of contact model in particle discrete element method． Chin J Ｒock Mech Eng，2009， 28( Suppl 2) : 4040 ( 王涛，吕庆，李杨，等． 颗粒离散元方法中接触模型的开 发． 岩石力学与工程学报，2009，28( 增刊 2) : 4040) ［21］ Tao B，Wu F Q，Guo G M，et al． Flexibility of visco-elastoplast￾ic model to rheological characteristics of rock and solution of rheological parameter． Chin J Ｒock Mech Eng，2005，24( 17) : 3165 ( 陶波，伍法权，郭改梅，等． 西原模型对岩石流变特性的适 应性及其参数确定． 岩石力学与工程学报，2005，24( 17) : 3165) ［22］ Tian H M，Chen W Z，Zhao W S，et al． Analysis of triaxial stress relaxation properties of red silty mudstone of Yichang-Ba￾dong highway． Ｒock Soil Mech，2013，34( 4) : 981 ( 田洪铭，陈卫忠，赵武胜，等． 宜--巴高速公路泥质红砂岩 三轴应力松弛特性研究． 岩土力学，2013，34( 4) : 981) ［23］ Fan Q Z． Experimental Study on Creep and its Disturbed Effect of Ｒocks ［Dissertation］． Qingdao: Shandong University of Science and Technology，2006: 30 ( 范庆忠． 岩石蠕变及其扰动效应试验研究［学位论文］． 青 岛: 山东科技大学，2006: 30) ［24］ Xu H F． Time depend behaviors of strength and elasticity modu￾lus of weak rock． Chin J Ｒock Mech Eng，1997，16( 3) : 246 ( 许宏发． 软岩强度和弹模的时间效应研究． 岩石力学与工 程学报，1997，16( 3) : 246) ［25］ Song D Z，Sun J． Nonlinear rheological property of rock materials and mechanic model． J Tongji Univ，1991，19( 4) : 395 ( 宋德彰，孙钧． 岩质材料非线性流变属性及其力学模型． 同济大学学报，1991，19( 4) : 395) ［26］ Lü X Q． The creep calculation method of variable parameters on · 738 ·

·838· 工程科学学报，第37卷，第7期 rock mass engineering.J Wuhan Unie Technol,1987(2):201 (孙钧.岩石流变力学及其工程应用研究的若干进展.岩石 (侣学清.岩体工程变参数流变计算方法.武汉工业大学学 力学与工程学报，2007,26(6)：1081) 报，1987(2)：201) 29]Xiong L X.Yang L D,Zhang Y.Non-stationary Burgers model 27]Ding Z K.Li A Z.Parameter identification of non-stationary for rock.J Cent South Unig Sci Technol,2010,41 (2):679 creep equation of rocks.Rock Soil Mech,2004,25 (Suppl 1): (熊良宵，杨林德，张尧.岩石的非定常Bugs模型.中南 37 大学学报：自然科学版，2010,41(2)：679) (丁志坤，吕爱钟.岩石粘弹性非定常蠕变方程的参数辨 B0]Yang C H,Mao J,Li J G.Modified creep model and its stabili- 识.岩土力学，2004,25（增刊1）：37) ty.J Jilin Unir Earth Sci Ed,2008,38(1):92 28]Sun J.Rock rheological mechanics and its advance in engineer- (杨彩红，毛君，李剑光改进的蠕变模型及其稳定性.吉林 ing applications.Chin J Rock Mech Eng,2007,26(6):1081 大学学报：地球科学版，2008,38(1)：92)

工程科学学报，第 37 卷，第 7 期 rock mass engineering． J Wuhan Univ Technol，1987( 2) : 201 ( 吕学清． 岩体工程变参数流变计算方法． 武汉工业大学学 报，1987( 2) : 201) ［27］ Ding Z K，Lü A Z． Parameter identification of non-stationary creep equation of rocks． Ｒock Soil Mech，2004，25( Suppl 1) : 37 ( 丁志坤，吕爱钟． 岩石粘弹性非定常蠕变方程的参数辨 识． 岩土力学，2004，25( 增刊 1) : 37) ［28］ Sun J． Ｒock rheological mechanics and its advance in engineer￾ing applications． Chin J Ｒock Mech Eng，2007，26( 6) : 1081 ( 孙钧． 岩石流变力学及其工程应用研究的若干进展． 岩石 力学与工程学报，2007，26( 6) : 1081) ［29］ Xiong L X，Yang L D，Zhang Y． Non-stationary Burgers model for rock． J Cent South Univ Sci Technol，2010，41( 2) : 679 ( 熊良宵，杨林德，张尧． 岩石的非定常 Burgers 模型． 中南 大学学报: 自然科学版，2010，41( 2) : 679) ［30］ Yang C H，Mao J，Li J G． Modified creep model and its stabili￾ty． J Jilin Univ Earth Sci Ed，2008，38( 1) : 92 ( 杨彩红，毛君，李剑光． 改进的蠕变模型及其稳定性． 吉林 大学学报: 地球科学版，2008，38( 1) : 92) · 838 ·