杨振伟等：基于颗粒流程序的非定常西原体模型 ·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 ·