基于塑性损伤的黏土岩本构模型及其数值实现.pdf

第4 5 卷第2 期 2 0 2 0 年2 月 煤炭 学 报 J O U R N A L0 FC H I N AC O A LS O C I E T Y V 0 1 ．4 5N o ．2 F e b ．2 0 2 0 移动阅读 李翻翻，陈卫忠，于洪丹，等．基于塑性损伤的黏土岩本构模型及其数值实现[ J ] ．煤炭学报，2 0 2 0 ，4 5 2 6 3 3 6 4 2 ．d o i 1 0 ．1 3 2 2 5 ／j ．c nk i ．j c c s ．2 叭9 ．0 1 6 7 L IF a n f a n ，C H E NW e i z h o n g ，Y UH o n g d a n ，e ta 1 ．C o n s t i t u t i v em o d e lo fc l a y s t o n eb a s e do np l a s t i cd 帅a g ea n di t sn u m e r i c a li m p l e m e n t a t i o n [ J ] ．J o u m a lo fc h i n ac o a ls o c i e t y ，2 0 2 0 ，4 5 2 6 3 3 6 4 2 ．d o i 1 0 ．1 3 2 2 5 ／j ．c n k i ．j c c s ．2 0 1 9 ． 0 1 6 7 基于塑性损伤的黏土岩本构模型及其数值实现 李翻翻1 ’2 ，陈卫忠1 ，于洪丹1 ’2 ，马永尚1 ’2 ，雷 江1 ’2 1 ．中国科学院武汉岩土力学研究所岩土力学与工程国家重点实验室，湖北武汉4 3 0 0 7 l ；2 ．中国科学院大学，北京1 0 0 0 4 9 摘要黏土岩是深部地下工程中常见的岩土介质，在施工过程中为确保巷道安全开挖，急需开展 与其力学特性相关的研究。已有试验研究结果表明黏土岩在三轴压缩试验过程中表现出明显的塑 性硬化、软化以及流动特性，同时发现试样在加载过程中内部会产生明显的缺陷，这些缺陷的存在 会在一定程度上削弱其力学性能。为建立能描述黏土岩力学现象的模型，以微裂纹作为损伤基元， 构建了合理的塑性损伤变量并将其与塑性硬化变量一同引入至修正D m c k e r P r a g e r 帽盖模型中， 建立了新的本构模型。通过A B A Q u S 平台及其u M A T 子程序对该模型进行数值实现，子程序在向 后欧拉本构积分算法的基础上引入塑性硬化和损伤变量参与迭代计算，使得屈服面在迭代过程中 与应力一同更新 屈服面的大小会随着塑性损伤与硬化变量的更新而变化 ，直至应力回到屈服面 上。采用该模型对常温条件下的黏土岩三轴压缩试验进行数值模拟，并将新模型的计算结果与试 验结果进行对比，对比结果显示模拟结果与试验结果大致吻合，说明本文所建立的本构模型可以较 好的反映黏土岩的力学特性。最后建立黏土岩隧硐开挖平面模型，采用新模型对开挖过程进行模 拟，并对整个模拟过程中损伤变量、应力以及孔隙水压力的演化规律进行分析，模拟结果验证了新 模型在地下硐室开挖数值模拟上的适用性。 关键词黏土岩；D m c k e r P r a g e r 准则；本构模型；积分算法；u M A T 中图分类号T u 4 5 文献标志码A 文章编号0 2 5 3 9 9 9 3 2 0 2 0 0 2 0 6 3 3 1 0 C o n s t i t u t i V em o d e lo fc l a y s t o n eb a s e do np l 硼t i cd a m a g ea n di t s n u m e r i c a l i m p l e m e n t a t i o n L IF a n f a n l ”，C H E NW e i z h o n 9 1 ，Y UH o n g d a n l ’- ，M AY o n g s h a n 9 1 ”，L E IJ i a n 9 1 ’2 1 ．S 衄￡eK 叮k 6 0 m ￡o 叮∥ 洳，懈地n 妇。蒯＆o ￡e 以n 池f E n 彰聊e 一愕，』瑚砒吡e 谚月。幽Ⅱn ds o 以朋e c 施n 如，啪i 珊船e A ∞d e ，Ⅳ∥ 曲n 鲫，ⅣⅡ舰n4 3 叫7 l ， 吼i ∞；2 ．‰妇m 蚵0 ，矾i M eA c o d e 彤o ，＆如n c e ，B e 拗g1 0 0 0 4 9 ，醌i M A b s t r a c t C 1 a y s t o n ei sac o m m o ng e o t e c h n i c a lm e d i u mi nd e e pu n d e 唱r o u n de n g i n e e r i n gp r o j e c t ．I no r d e rt o e n s u r et l l e s a f e t yo ft u n n e ld u r i n gi t se x c a v a t i o n ，t h ew o r kw h i c hi sr e l a t e dt oi t sm e c h a n i c a lp r o p e r t i e si sr e q u i r e dt ob es t u d i e d ． T h ee x i s t i n gt e s tr e s u l t ss h o wt h a tt h ec l a y s t o n ee x h i b i t ss o m eo b v i o u ss t r a i n - h a r d e n i n g ，s 缸a i n s o f t e n i n ga n dp l a s t i c n o wc h a r a c t e r i s t i c sd u r i n gt h et r i a X i a lc o m p r e s s i o nt e s t ．A tt h es a m et i m e ，i ti sf o u n dt h a tt h es a m p l e sp r o d u c es o m e o b v i o u sd e f b c t sd u r i n gt h el o a d i n gp r o c e s s ，a n dt h ee x i s t e n c eo ft h e s ed e f b c t sw o u l dw e a k e nt h em e c h a n i c a lp r o p e r t i e s ． I no r d e rt oe s t a b l i s ham o d e lw h i c hc a nd e s c r i b et h em e c h a n i c a lp h e n o m e n ao fc l a y s t o n e ，t h em i c r o c r a c ki su s e da s t h ed a m a g ee l e m e n tt oc o n s t r u c tar e a s o n a b l ep l a s t i cd a m a g eV a r i a b l e ．T h ed a m a g eV a r i a b l ei si n t r o d u c e di n t ot h em o d 一 收稿日期2 0 1 9 一0 2 一1 2修回日期2 0 1 9 一0 6 2 3责任编辑郭晓炜 基金项目国家自然科学基金资助项目 5 1 4 7 9 1 9 0 ；中国科学院青年创新促进会基金资助项目 作者简介李翻翻 1 9 9 2 一 ，男，安徽六安人，博士。E m a j l 1 3 6 6 5 6 9 9 3 2 6 s i n a ．c n 通讯作者陈卫忠 1 9 6 8 一 ，男，江苏南通人，教授。E m a i l w z c h e n w h r s m ．a c ．c n 万方数据 煤炭学报 2 0 2 0 年第4 5 卷 i f i e dD m c k e r - P m g e rc a pm o d e lt o g e t h e rw i t ht h ep l a s t i ch a r d e n i n gv a r i a b l et oc o n s t I u c tan e wc o n s t i t u t i V em o d e l ．T h e m o d e li sn u m e r i c a l l yi m p l e m e n t e db yt h eA B A Q U Sp l a t f o 瑚a n di t sU M A Ts u b r o u t i n e ．I nt h es u b r o u t i n e ，t h ep l a s t i c h a r d e n i n ga n dd a m a g ev 耐a b l e sa r ei n t r o d u c e di n t ot h ei t e r a t i v ec a l c u l a t i o nb a s e do nt h eb a c k w a r dE u l e rc o n s t i t u t i v e i n t e 铲a la l g o r i t h m ，s ot h a tt h ey i e l ds u d A c ei su p d a t e dw i t ht h es t r e s si nt h ei t e r a t i v ep r o c e s su n t i lt h es t r e s sr e t u m st o t h ey i e l ds u I f a e e T h es i z eo ft h ey i e l ds u r f a c ec h a n g e sw i t ht h ep l a s t i cd a m a g ea n dh a r d e n i n gv a r i a b l e ．T h i sm o d e li s u s e dt os i m u l a t et h et r i a x i a lc o m p r e s s i o nt e s to fc l a y s t o n eu n d e rn o 瑚a lt e m p e r a t u r e ，a n dt h es i m u l a t i o nr e s u l t sa r e c o m p a r e dw i t ht h ee x p e r i m e n t a lr e s u l t s ．T h ec o m p a r i s o nr e s u l t ss h o wt h a tt h es i m u l a t i o nr e s u l t sa r ei ng o o da g r e e m e n t w i t ht h ee x p e r i m e n t a lr e s u l t s ，i n d i c a t i n gt h a tt h ec o n s t i t u t i v em o d e le s t a b l i s h e di nt h i sp 印e rc a nb e t t e rr e n e c tt h em e - c h a n i c a lp m p e n i e so fc l a y s t o n e ．F i n a U y ，t h ep l a n ee x c a v a t i o nm o d e lo fc l a y s t o n et u n n e li se s t a b l i s h e d ．T h en e wm o d e l i su s e dt os i m u l a t et h ee x c a v a t i o np r o c e s s ．T h ee v o l u t i o no fd a m a g ev a r i a b l e s ，s t r e s sa n dp o r ew a t e rp r e s s u r ed u r i n gt h e w h o l es i m u l a t i o np m c e s si sa n a l y z e d ．T h es i m u l a t i o nr e s u l t sv e r i f yt h ea p p l i c a b i l i t yo ft h en e wm o d e li nn u m e r i c a ls i m u l a t i o no fu n d e r g r o u n dt u n n e le x c a v “o n ． K e yw o r d s c l a y s t o n e ；D m c k e r P r a g e rc r i t e r i o n ；c o n s t i t u t i V em o d e l ；i n t e g r a t i o na l g o r i t h m ；U M A T 黏土岩是深部地下工程中常见的岩土介质，这一 介质具有非常复杂的物理力学性质⋯。为在施工过 程中选择合理的支护措施L 2 J ，需对其力学特性进行 研究。L I ∞o 通过c T 扫描技术对黏土岩试样进行观 测后发现，在各向等压渗透试验后黏土岩内部会产生 大量的微裂纹。因此构建合理的本构模型来反应其 力学特性对深部地下工程施工尤为重要。 近些年损伤力学的发展为这种加载过程中内部 产生缺陷的岩土材料的力学性能研究提供了一种方 法。细观损伤力学采用平均化的方法将材料的细观 损伤机制反映到宏观的力学行为中，这种方法既避免 了复杂的统计学计算，又可以提供很好的物理背景。 陈新等Mo 将岩土材料假设成含微孔洞的球形单 元集合体，将微孔洞作为损伤基元，通过微孔洞所占 比例大小来衡量损伤程度，提出了适用于花岗岩的塑 性损伤模型。韦立德等”1 综合考虑剪胀效应对岩体 变形和渗流的影响以及损伤对于岩体渗流的影响提 出了岩石损伤本构模型。谢和平等M 1 认为岩石的变 形与破坏是一个与外界能量交换的过程，基于能量耗 散原理建立了损伤演化的宏一细一微观多层次岩石力 学体系。s H A 0 等“ o 对岩石的损伤演化与微裂纹发 展进行了系统的研究，建立连续损伤模型，分析非饱 和条件下损伤对于岩石渗透系数的影响并建立流固 耦合有限元模型。w A N G 等旧。通过分析强度准则和 残余强度对岩土材料的软化影响建立了统计损伤模 型并考虑了软化效应。王军详等一1 在弹塑性理论和 损伤理论的基础上引入修正有效应力原理建立适用 于脆性岩石的弹塑性损伤本构模型。袁小平等刨认 为软化分为塑性软化和损伤软化，损伤软化是由于岩 石内部微裂纹的扩展引起。基于D r u c k e r P m g e r 准 则建立适用于脆性岩石的弹塑性损伤本构模型。陈 亮等u 基于北山花岗岩的试验结果，通过岩石的宏 观力学性能推导岩石在外力作用下微裂纹的发展规 律，在热力学框架的基础下提出了适用于花岗岩的唯 象弹塑性损伤本构模型。 笔者以黏土岩为研究对象，将黏土岩内部的微裂 纹作为损伤基元在修正D m c k e r P r a g e r 帽盖模型的 基础上引入了塑性损伤和硬化，构建了适用于黏土岩 的本构模型。基于有限元软件A B A Q u s ，在u M A T 子 程序中将塑性硬化和损伤引入迭代算法中编制了该 模型的本构程序。通过对常温条件下黏土岩的三轴 压缩试验进行数值模拟，并与试验结果进行比较，对 比结果表明新模型不仅可以很好地反映黏土岩的硬 化和软化现象，还很好的反映了黏土岩的塑性流动特 性。最后将新模型应用在黏土岩隧硐开挖算例上，验 证了新模型在地下工程施工数值模拟上的适用性。 1 基于塑性损伤的黏土岩本构模型 龚哲，于洪丹等2 。1 副等对黏土岩进行了一系列 的三轴压缩试验，试验结果表明黏土岩是一种具有明 显硬化和软化现象的介质，且整个试验过程的应力一 应变曲线可以划分为4 个阶段，分别为弹性阶段、塑 性硬化阶段、塑性软化阶段以及塑性流动阶段。 1 ．1屈服函数与塑性势函数 B A L D I 等4o 在总结大量试验结果的基础上得到 了黏土岩的屈服面，发现D 1 1 l c k e r P r a g e r 对于黏土岩 的力学行为描述更加准确。在修正D m c k e r P r a g e r 准则的基础上引入“帽盖模型”来描述黏土岩在静水 压力下屈服的特性5 | 。修正D r u c k e r P r a g e r 帽盖模 型的屈服函数由剪切屈服面和帽盖屈服面两段组成， 连接处光滑，在p L q 平面内的屈服面如图1 所示，图 1 中，6 叩。 七 。 万方数据 第2 期 李翻翻等基于塑性损伤的黏土岩本构模型及其数值实现 6 3 5 面 图1修正D r u c k e r - P r a g e r 帽盖模型 F i g ．1 M o d i f i e dD r u c k e r P r a g e rc a pm o d e l 其中剪切屈服面F 。的表达式为 F 。 q 一叩’一尼 0 1 式中，p ，_ t r 盯 ；g 争s 2 ，P a ；t r 为一种运 算符号，求解矩阵的迹，即对角线元素之和，s 盯一 筇，仃为应力张量，s 为偏应力张量，6 为K r o n e c k e r 函 数；o 和庇分别为剪切屈服面在p7 一g 平面上的斜率和 截距，其表达式为 6 s i n D 旷苹 2 忌．鱼 一 式中，妒，c 分别为内摩擦角 。 和黏聚力 P a 。 帽盖屈服面F 。为一段椭圆弧，其表达式为 耻√9 2 华 一 叩。“Ⅲ 3 其中，R 为控制帽盖几何形状的参数；p 。的表达为 旷篇 ㈩p a2 而 4 式中，p 。为先期固结应力，P a 。 传统塑性力学中假设塑性势函数与屈服函数采 用相同的数学表达式，即在计算过程中采用关联流动 法则。但是关联流动法则并不适用于岩土材料，已经 有证据证明岩土材料的屈服函数与塑性势函数并不 相同，因此塑性计算时应当采用非关联流动法则6 f 。 本文假设塑性势函数与屈服函数形状相似但数 学表达式不相同口7 | ，其两段表达式G 。和G 。分别为 G 。 g a ，p ’一尼7 0 5 G 。 √ p7 一p 。 2 尺q 2 6 式中，G 。，G 。分别为剪切屈服面F 。对应的塑性势函数 和帽盖屈服面F 。对应的塑性势函数；o ’和七’表达式 为 式中，9 为膨胀角， 。 。 1 ．2 硬化变量及损伤变量演化方程 本文考虑使用损伤理论来描述黏土岩内部的 缺陷8 ‘1 9 ] 。损伤理论从细观力学的角度解释了材 料力学性能的弱化，黏土岩在加载过程中会产生很 多微裂纹啪≈2 | ，微裂纹在加载过程中不断地扩展与 聚集，从而引起材料力学性能的下降心3 | 。已有试验 研究也证实了黏土岩在加载过程中存在微裂纹扩 展的现象3 ’2 4 ] 。笔者将以微裂纹作为损伤基元构 建损伤变量，并引入到本构模型中来描述黏土岩的 力学特性。 1 ．2 ．1 塑性损伤变量 在塑性阶段认为塑性损伤是存在于整个塑性过 程中的，选择合理的损伤演化方程对于描述黏土岩的 塑性软化和塑性流动十分重要。随着塑性应变的不 断累积，塑性损伤不断增大。然而，由于明显的塑性 流动阶段的存在，假设的塑性损伤不能随着塑性应变 一直增大，而应存在一个上限值，接近上限值时塑性 损伤几乎不变，材料进入塑性流动阶段。塑性损伤变 量D 。演化方程为 D ， D ，一D ，e x p 一卢，9 1 8 式中，D 。为塑性损伤变量D 。的上限值；d ，和戽为模 型参数；亭。为等效塑性应变。 1 ．2 ．2 塑性硬化变量 黏岩的应力一应变曲线的塑性阶段在宏观上 先表现为硬化。因此，为了描述其塑性硬化力学特性 引入合理的硬化法则十分重要。本文采用双曲线型 塑性硬化方程，塑性硬化变量九。表达式为 危， 危田 i { 立畿 9 式中，庇p 0 和6 。为模型参数。 S C H M E R T M A N n 和O S T E R B E R G [ 25 | ，I A M B E [ 2 6 ] 以及H A J I A B D O L M A J I D 等‘”1 对不同的岩土材料研 究发现在三轴压缩试验过程中，当材料进入塑性阶段 后内摩擦角随着应变的累积呈负指数的形式升高，此 后并没有明显的降低，而黏聚力在整个塑性过程中随 着应变的累积先增大后减小。 因此模型中黏聚力c 和内摩擦角‘；p 的演化方程 采用 f c 1 一D 。 p c o i 9 a 。c t a n ，t a n 妒。 ‘1 0 ’ 式中，c 。为黏聚力的初始值，P a ；‘；p 。为内摩擦角的初 始值， o 。 一9 一日 9一n n ．m一．5； ，“ 一．川 鱼一j 一 ～一～∞一～ 万方数据 6 3 6 煤炭学报 2 0 2 0 年第4 5 卷 先期固结应力p 。的演化受塑性体应变的影响， 其表达式为 p 。 e x p c p 8 ； p c o 1 1 式中，p 棚为先期固结应力的初始值，P a ；c 。为模型参 数；s 为塑性体应变。 2 本构积分算法 进入塑性阶段时，材料本身的屈服面会随着塑性 应变的累积 硬化变量和损伤变量的变化 不断地发 生变化 硬化过程屈服面扩张／软化过程屈服面收 缩 ，且在屈服面变化的整个过程应力始终位于屈服 面之上。第1 节给出了塑性阶段屈服面的变化规律， 本节将介绍在屈服面变化过程中如何通过数值方法 使得应力始终位于屈服面之上。 塑性阶段计算时，先假设只发生弹性变形，计算 弹性试探应力，当试探应力超出屈服面时，需要通过 本构积分算法心8 。3 0 1 将超出屈服面的应力拉回到屈服 面上。 本文采用的积分算法除了考虑塑性硬化变量之 外，还引入了塑性损伤变量，在算法迭代过程中塑性 硬化变量和塑性损伤变量将随着塑性应变的变化与 应力一起进行调整，并通过更新后的塑性硬化变量和 塑性损伤变量对屈服面进行更新，从而使得每一个增 量步更新后的应力始终位于更新后的屈服面上 图 2 。 后鼍K 徭算法 三鬻雾蓬劳磊鋈 然 增量步开始、稿滗、 时的屈服面l 憾，褥l a 硬化阶段 b 软化阶段 图2 本构积分算法示意 F 培．2 S c h e m a t i cd i a g r a mo fc o n s t i t u t i v ei m e g r a la l g o r i t h m 2 ．1 向后欧拉式本构积分算法 由弹塑性增量理论可知，应变增量可以写成 △s △s 。 △8 9 1 2 式中，△￡，△￡。，△矿分别为总应变增量、弹性应变增量 和塑性应变增量的二阶张量。 本文采用非关联流动法则，则塑性应变增量的形 式可以写成 △8 p d A i 业 1 3 ‘a 盯 其中，d | L 。为塑性因子；盯为应力二阶张量；G 为塑性 势函数。塑性势函数求导如下 丝一a ，型 亟 1 4 2 一a 。。十。 Ll 斗J 皇鱼羔一旦二里 鱼芷。 1 。一一T 汩、厄丁i y 了琢矿a 盯 塑一亟 1 5 、厄厂i 了了丽a 仃 。’ 式 1 5 中平均主应力p7 和偏应力q 对应力张量 的求导结果写成向量矩阵的形式 鼍 ⋯1 ，1 ’0 ，0 ，0 ] ’ 粤氧 力圳2 丁。丁，] t ‘1 6 ’ d 盯 2 盯 式中，盯 以，其中，厶为偏应力第二不变量；s 。 盯，一 盯。 盯， 盯 ／3 ；s ， 盯y 一 盯。 盯y 盯 ／3 ；5 盯一 盯。 盯， 盯； ／3 ，盯，，盯，，盯，丁w ，丁。，丁“ 为6 个应力 分量。 当增量步开始时，假设该增量步内的总应变增量 全部是弹性应变增量，则可根据下式计算第n 1 增 量步的弹性试探应力。暑。 以1 呒 D 。△￡川 1 7 式中，D 。为刚度矩阵的张量形式；叽为上一增量步 结束时应力张量；△s 川为第n l 增量步的应变增 量张量。 当试探应力超出屈服面时，需要通过向后欧拉算 法将其拉回到屈服面上。进入塑性阶段后实际应力 为 叽 1 D _ D 。 △￡。 】一△s 1 1 8 式中，△￡ 。为第乃 1 增量步的塑性应变增量张量。 将式 1 3 ， 1 7 代入式 1 8 可得 ％2 叱“k ，D e ．老L 1 9 式中，△． L 州为更新后的塑性因子。 万方数据 第2 期李翻翻等基于塑性损伤的黏土岩本构模型及其数值实现 6 3 7 F [ p 7 川，q 川， D ， 川， h ， 川] 0 2 0 将式 1 3 代入式 8 ， 9 中可得到塑性硬化变 式中， p ’ 川，q 川， D ， 川和 危。 川分别为更新后的 量和塑性损伤变量 。矗，，。。田兰_二二二_二_二』兰圣兰三兰三二竺竺兰』l兰[≠兰兰兰兰兰戛 。2 。， ‘”㈨ 6 p √扣心川瓢。 ． 订‰甄。’ 卜危p 0 6 r √亏 ￡ △入川嘉㈠ ￡ △入川荔∽ ‰- 瓦一瓦唧h 层∽Ⅷ川瓢。 ． “ ‰甄．1 2 2 将式 2 1 和式 2 2 代入式 2 0 中可以消除未知 量 D 。 川和 ， 川。则式 1 9 和 2 0 变成求解 吒 ，和△入川的隐式方程组，采用牛顿法求解得 △． L 川和更新后的应力叽 ．。将△入川代回到式 2 1 和 2 2 更新硬化变量和损伤变量再带回式 2 0 更新 屈服面。 2 ．2 一致性切线刚度矩阵 由前人的研究结果可知阳，在本构积分算法中 一致切线刚度矩阵可以在提高计算精度和稳定性的 同时加速收敛。 向后欧拉算法式 1 9 可以写成矩阵形式 盯 。 盯l 。 i 一△A 。 。D 。 生} 2 3 对式 2 3 进行微分可得 训川诅乩。叫川。。老L 一 扒 1 D 。耋I 扎 。 2 4 对式 2 4 进行化简可以得到 扎。 I 心川。。害 ～” ㈠川一儿蔫 2 5 式中，J 为单位矩阵。 根据一致性准则可得 叫烈 。】1 扎。 差卺 蓑鲁 天k l √手l 羞[ 。 2 6 将式 2 5 代入式 2 6 并化简可得 天 j 拉 2 7 J L } 竺L 一 2 7 z 里 M a 盯 I 式中小蚓州】7 J △k - 。。妾 ～既M 舞 装 周磊。＼a D 。a 孝，a p a 亭，／√3 a 盯Ⅱl 。 l 。 将式 2 7 代入式 2 5 中可以得到 扎． J 心州。。妾 ～ 互一』生 z 』里 M a 盯 d 盯 l 2 8 则一致切线刚度矩阵的表达式为 ．2 门．一l D 。。 f J △入。 l D 。i ≥气1 D 。 、 a o 。n l ， f ，1 一墨．』生、I I ‘z 旦 Ma 盯圳 2 9 L a 盯 1 ／ 式中，D 。为一致性刚度矩阵。 2 ．3I 阳性A T 子程序调用流程 A B A Q u S 提供了用F o r t r a n 语言编写的U M A T 子 程序接口，供用户二次开发之用。A B A Q u s 每个增量 步计算流程大致如下 图3 1 增量步开始时，A B A Q u S 会将应力、应变、应 变增量以及时间增量传输给u M A T 子程序，同时通 过P R O P s n 将固体参数导入子程序中，通过子程序 计算得到新的应力、应变以及一致切线刚度矩阵返回 主程序中。 2 A B A Q u s 主程序再调用流体参数并通 过s O I L S 分析步进行平衡迭代计算。 3 如若迭代不收敛或者迭代达到最大次数，则 减小增量步时长重新开始计算；如若迭代收敛则更新 应力、应变以及其他参数，进入下一增量步开始计算。 3 黏土岩三轴压缩试验数值模拟 编写新模型的u M A T 子程序，经编译调试后对 常温条件下黏土岩的三轴压缩试验进行有限元模拟， 并与室内试验结果进行比较验证该子程序的有效性。 3 ．1 试验结果 常温条件下，黏土岩三轴压缩试验在双联动软岩 万方数据 6 3 8 煤炭 学报 2 0 2 0 年第4 5 卷 图3A B A Q u S 调用u M A T 子程序计‘算流程 F i g ．3 C a l c u l a t i o nn o wc h a r to fA B A Q U Sc a U sU M A T s u b m u t i n e 渗流一应力耦合流变仪上进行3 I 。试验考虑不同初 始围压条件 4 ．7 ，3 ．7M P a ，初始孔压1 ．2M P a ，采用 不排水压缩 压缩速率2 0 n ∥m i n ，得到应力一应变 曲线和超孑L 压曲线，如图4 所示。从图4 可以看出， 黏土岩进入塑性阶段之后，除了有非常明显的塑性硬 化外，还有明显的塑性软化以及塑性流动现象；整个 加载过程中试样的孔隙水压力在逐渐增大，且进入塑 性阶段后，其增长速率逐渐减小。 3 ．2 模型参数反演及模拟结果对比 建立黏土岩试样有限元模型 图5 ，模型直径为 3 8m m 、高度为7 6m m ，模型采用H M 耦合4 节点轴对 称减缩积分单元 c A x 4 R P ，约束试样底部轴向自由 度和轴心径向自由度。模型四周为不排水边界。按 照试验流程施加相应的围压和孔压，再对模型顶部施 加轴向位移。 表1 中给出了黏土岩数值模拟的主要物理参数。 这些参数是通过大量的现场试验、室内试验以及现场 声波探测得到的结果、32 | 。模型中有关参数主要通过 反演获得。反演过程中的对比试验数据为三轴压缩 试验所得应力一应变曲线，反演采用贾善坡。33 。提出的 N e l d e r M e a d 法与有限元联合反演法，将有限元程序 作为一个单独模块嵌入N e l d e r M e a d 算法程序中，通 a 应变曲线 轴向应变 b 超孔压曲线 图4 黏土岩不同围压下三轴压缩试验曲线⋯l F i g ．4 T r i a x i a lc o m p r e s s i o nt e s tc u r v e so fc l a y s t o n eu n d e r d i f k r e n tc o n n n i n gp r e s s u r e s [ 1 3 一 轴． 、 对．j一 称1 固定约束 JU |『J 法向约束 围 压 图5 试验模拟有限兀模型 F i g ．5 F i n i t ee l e n l e l l tI n o d e lo ft e s ts i m u l a t i o n 过数值算例对反演结果进行求解。在整个反演过程 中，需要对比的数据是试验与计算所得的相同应变下 对应的应力值与超孔隙水压力值，为了同时利用这两 种数据，需要为两种数据设置一个权重关系，权重的 大小关系与试验过程中测得该类数据的可靠度有关。 其表达式为 1 一入 妒l 入 D 2 一m i n 3 0 式中，人为两种误差函数所占比重；‘P ，， p 分别为试 验值与计算值在偏差应力和孔压上的误差函数；当 入 1 时，则认为第1 种物理量所测数据不可信，当 A 0 时，则认为第2 种物理量所测数据不可信。 万方数据 第2 期李翻翻等基于耀性损伤的黏土岩本构模型及其数值实现 6 3 9 表l黏土岩基本力学参数‘3 2 T a b I elB a s i cm e c h a n i c a lp a r a m e t e r sO fc l a y s t O n e [ 3 2 物理参数数值 弹性模量／M P 8 剪切模矗} ／M P a 黏聚力／M I ，a 先期固结应力／M P 8 膨胀角／ 。 内摩擦角／ o 渗透系数／ m s “ 泊松比 干密度／ k g m 。 孔隙比／％ f L 隙水密度／ k g m _ 3 7 L 隙水体积模量／G P a E h o 14 0 0 ，E v 0 7 0 0 G 。1 ．o 2 8 0 0 ．4 5 ．O O 1 8 血 6 1 0 1 2 ，丘、 3 1 0 一1 2 ‰ O ．2 5 ，“ 、l ， 0 ．1 2 5 l6 4 0 7 0 9 9 9 2 ．0 由计算所得和试验所得的偏差应力数据构成的 误差函数表达式为 一 ， 妒， D 。，，1 8 p ， 。。，⋯⋯ ∑ 盯罗m 一盯≯1 2 3 1 】 式中，盯≯为不同的轴向应变所对应的偏应力试验 值；盯∥为相应的反演分析计算值；z 为参与反演分 析的试验值个数。 由计算所得和试验所得的超孔隙水压力数据构 成的误差函数表达式为 一 f ‘P D “ ％危⋯⋯。 ∑ 蒯。一搿 2 3 2 k l 式中，p 为不同的轴向应变所对应的试验孔压值； p 滞为相应的反演分析计算值；叫为孔隙水压力。 模型参数的反演结果见表2 [ 2 2 ] 。 表2 黏土岩本构模型参数m 3 T a b l e2 C I a y s t o n ec o n s t i t u t i v em o d e lp a r a m e t e r s [ 2 2 模型参数数值模型参数数值 肛M O0 3 5，0 ．0 5 D 。 O ．3 5 J 3 。 23 8 0 n ．2 ．8 9尺4 ．4 2 图6 给出了新模型计算得到的不同围压条件下 不排水三轴压缩试验与模拟应变曲线的对比结果。 从图6 可以看出塑性阶段模拟结果与试验结果相似 表现为先硬化再软化；塑性阶段后期试验结果与新模 型的模拟结果均表现出明显的塑性流动特性，这是因 为塑性损伤存在损伤上限，试样内的损伤不可能无限 制的增加，当试样内的损伤达到上限值时，宏观上表 现为塑性流动阶段。 图7 给出了试验和模拟的超孔压曲线对比图。 从图7 可以看出，模拟的超孑L 压结果可以较好的反应 超孔压的演化趋势，即进入塑性阶段后，超孔压增大 的速率逐渐减小。 径向应变 轴向应变 a 围压47 M P a 4 ／／一、、3 ／，一～、、 一／7 垩j ／、＼