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

工程科学学报，第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 ·