·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 ·