软黏土不排水循环应力应变关系的数值模拟.pdf

第 3期 2 0 1 5年 6月 水利水运 工程学报 AND ENGD旺 RI NG No 3 J u n 2 0 1 5 D O I ：1 0 1 6 1 9 8 j c n k i 1 0 0 9 6 4 0 X 2 0 1 5 0 3 0 1 2 王哲学 ， 王建化，程星磊 软黏土不排水循环应力应变关系的数值模拟 J 水利水运工程学报， 2 0 1 5 ( 3 ) ：8 1 8 7 ( WA N G Z h e x u e ， WA N G J i a n h u a ，C H E N G X i n g l e i N u me r i c a l s i mu l a t i o n o f u n d r a i n e d c y c l i c s t r e s s s t r a i n r e s p o n s e o f s o f t c l a y J H y d r o S c i e n c e a n d E n g i n e e r i n g ，2 0 1 5 ( 3 ) ： 8 1 8 7 ) 软黏 土不排水循环应 力应 变关 系 的数值模拟 王哲学 一，王建华 一，程星磊 ， ( 1 天津大学 水利工程仿真与安全国家重点实验室， 天津3 0 0 0 7 2 ； 2 天津大学 岩土工程研究所，天津 3 0 0 0 7 2 ) 摘 要 ：对程星磊等提出的总应力形式增量弹塑性本构模型进行二次开发， 以模拟复杂应力状态下软黏土的响 应。通过 N e w t o n R a p h s o n算法， 对材料非线性问题进行迭代求解； 针对本构模型中应力反向等一系列关键性问 题， 应用欧拉切线算法编写了有限元程序， 并结合子增量方法提高了计算精度。预测了软黏土在轴对称应力状 态下 的响应 ， 得到了应力应变关系 曲线 ， 将其与单元预测结果进 行 比较 ， 二者 趋于一致 ， 从而 验证 了有 限元程序 编写的合理性。利用该程序模拟三轴不固结不排水试验， 模拟结果与试验吻合良好， 表明该本构模型有限元程 序可以反映轴对称应力状态下软黏土的不排水应力应变特性 ， 可应用于更加复杂边值问题的模拟计算。 关 键 词 ： 循环荷载； 饱和软黏土； 应力应变； 本构模型； 欧拉切线算法； 三轴试验 中图分类号： T U 4 3 文献标志码 ： A 文章编号 ： 1 0 0 9 6 4 0 X( 2 0 1 5 ) 0 3 - 0 0 8 1 - 0 7 嵌固在海洋软黏土中的锚 固基础在风 、 浪等荷载作用下 ， 为维持上部结构工作性能， 必须保持 自身稳定 ， 因此需要合理分析软黏土中锚 固基础在循环荷载作用下的变形失稳过程。锚固基础的变形不仅与基础本身 有关 ， 海洋中软土地基对其影响更为重要。因此研究循环荷载下软黏土的不 固结不排水应力应变关系具有 重要意义。J H P r 6 v o s t 1 - 2 J 提 出了基于运动硬化塑性理论 的多面模型 ， 但模型的变形分析很复杂 ， 需要记忆 所有假设的多个次加载面形状 ， 并分析它们的移动 ， 导致在数值应用 中计算过程复杂 。Y F D a f a l i a s 等 _ 3 j 建立了包括 1 个正常屈服面和 1 个次加载面的两面模型 ， 大大简化 了多面模型 ； 后来他们又将两面模型中的 次加载面简化为一点 ， 提出了单面模型 。两面和单面模型都 叫边界面模型_ 6 J 。由于改进初期 的边界面 模型映射中心 固定在原点 ， 无法合理描述土单元 的应力应变滞 回曲线 J ， 且硬化模 量场 的演化过程表述较 为复杂， 不便 于应用到数值分析中。程星磊等 。 。 采用总应力形式在偏应力空问中构建硬化模量场 ， 硬化模 量插值函数数学表达简洁 ， 引入模型参数较少 ， 硬化模量场 的演化规律表述简单 ， 跟踪循环应力路径所需记 忆 的参数较少 ， 便于计算应用 。通过上述分析 ， 针对程星磊等建立 的软黏土循环 弹塑性本构模型进行二次开 发 ， 并对循环三轴试验进行数值模拟 ， 通过计算结果与试验结果的比较 ， 验证该本构模型数值分 析方法的可 行 性 1 本构模型简述 采用 M i s e s 屈服准则构建边界面约束方程 ： F = 1 5 i ； 一A =0 收稿 日期 ： 2 0 1 4 - 0 9 3 0 基金项 目：国家 自然科学基金资助项 目( 5 1 1 7 9 1 2 0 ) ； 博 士点基金项 目( 2 0 1 3 0 0 3 2 1 1 o o 4 5 ) 作者简介：王哲学( 1 9 8 9 一) ，男， 河北承德人 ， 硕士研究生 ， 主要从事土动力学与海洋岩土工程研究。 E - m a i l ： w a n g z h e x u e 0 1 1 7 1 2 6 e o m通信作者 ： 王建华 ( t d w j h t j u e d u c n ) 学兔兔 w w w .x u e t u t u .c o m 8 2 水 利 水 运 工 程 学 报 2 0 1 5年 6月 式 中： ； 为边界面上偏应力张量 ； A 。为边界面半径。 依据经典弹塑性理论 ， 有总偏应变增量计算式 ： d d s + 3 S ij- z ( 2 ) 式中： d e 为总偏应变增量 ； G为弹性剪切模量 ； H 为塑性硬化模量 ； d s 为当前应力点的偏应力增量 。 引入弹塑性硬化模量 ， 令 ： 1 H =1 ( 2 G)+1 H ( 3 ) 在初始加载或应力反 向时， 塑性变形模量为无穷大 ： H =H =2 G ( 4 ) 图 1给出了应力空 间中加载和卸载时确定弹塑性 模量的映射规则 图解 。这里边界面针对九维应力空 间 中的超 曲面。图中 s 为映射 中心点 ， s 为 当前应力点 ， S 为通过映射 中心点和当前应力点 的射线 与边界 面的 交点 ， 此点称为像应力点 ， s 为应力反向点。当发生应 力反向时 ， 如图 1 ( b ) ， s 则变为映射中心 ， 其他各点定 义不变 ， 此时会出现新的应力反 向点 s 。定义 和 分别为当前应 力点 和初始加 载点 到像 应力点 的距离 。 弹塑性模量插值规则为 ： ( a ) 加载时 ( h ) 卸载时 图 l 弹塑性模量 的映射规则 图解 F i g 1 El a s t o - p l a s t i c ma pp i n g r u l e f o r l o a d i ng a n d u nl o a di n g H=( t3 6 o ) H ( 5 ) 式中： 日 为弹塑性模量最大值 ， 取 2 G； 为反映剪应变累积速率和大小的模型参数 。 在不固结不排水条件下 ， 边界面半径 A 。 和弹性模量 G可通过不 固结不排水循环三轴试验得到。在三轴 应力状态下 ， 式( 1 ) 简化为下式 ： F=( ： 一 ) 一A =0 ( 6 ) 式中 ： 为试样轴向应力 ； 为试样径向应力。 由于初始应力应变响应均遵循同一初始加荷 曲线 ， 且可用式( 7 ) 拟合 ， 则拟合曲线 的渐近线所对应 的强 度就是边界面半径 A 。 ， 其值等于 1 b； 曲线 的初始切线模量记为 E， 其值等于 1 a。另外 ， 弹性 剪切模量 G 可 由式 ( 8 ) 求得 ： = ( 口+6 ) ( 7 ) G=E 2 ( 1+ ) ( 8 ) 式 中： 口和 b 为拟合 函数 中的待定系数 ； 为泊松比， 不排水时取 0 5 。 参数 可以用下式表示 ： = x 0 d y 8 。 ( 9 ) 式 中： 为反映累积偏应变大小的待定系数 ； d y 。 为八面体累积剪应变增量 ， 以此描述偏应变累积变化 的速 率和趋势。对试验数据采用式( 1 0 ) 拟合 ： 8 。=c ( 1 0 ) 式中： c ， d为待定系数 ； 为循环次数 。则有 d 7 。=c d N d ( 1 1 ) 下面结合某一具体三轴试验 ， 给出待定系数 c， d， I x 的确定方法。图 2给出了该三轴试验得到的静应 力 比1 8 = 和 5时八面体累积剪应变随振次的变化关系。用式 ( 1 0 ) 拟合试验数据 ， 如 图中虚线所 ,aT 0 3 0 示 ， 得到不同静应力 比T S ,a T 不同循环应力 比 8 c 下 的待定系数 c 和 d值 。通过式 ( 1 1 ) 可得应力循 环一周八面体累积剪应变增量 。 。 学兔兔 w w w .x u e t u t u .c o m 第 3 期 王哲学 ， 等 ： 软黏土不排水循环应力应变关系 的数值模拟 8 3 1 00 2 0 0 3 0 0 4 O 0 5 0 0 6 0 0 7 0 0 8 0 0 9 00l 0 0 0 _7v ( a ) 5=0 3 1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 6 0 0 7 0 0 8 0 0 9 0 0 1 0 0 0 _7v ( h ) 5 = 0 5 图 2 三轴压缩时八面体累积剪应变随循环次数的变化关系曲线 F i g 2 R e l a t i o n c u r v e s o f o c t a h e d r a l c u mu l a t i v e s h e a fi n g s t r a i n a n d c y c l e t i me s i n t r i a x i a l c o mp r e s s i o n 进一步 ， 可通过试算方法确定 。 ， 即利用式 ( 9 ) 将 8 D 代人式 ( 5 ) ， 给定一个 。 ， 通过式( 2 ) 计算得到破 坏振次下的八面体累积剪应变。这里规定按照静偏应力作用下试样的八面体静剪应变与静偏应力与循环应 力共 同作用下八面体循环累积剪应变之和达到 1 4 的标准确定循环破坏次数 。若模型计算所得破坏振 次 下的八面体累积剪应变与试验结果一致 ， 那么此时的 值 即为该试验应力条件下待定系数 的取值。 这样 ， 不 同静应力 比 T S 8 ， ， 不同循环应力比 ， 下 的参数 c ， 。 都 已给出( 见表) 。则对任意 循环应力 比下的参数值 ， 可通, a 过 T 线性插值方法确定 。 8 ,cyT d 1 表 1 不同静应力比、 循环应力比下拟合函数中的待定系数 C ， d及p。 取值 T a b 1 Un d e t e r mi n e d c o e f f i c i e n t s C a n d d o f fi t t i n g f u n c t i o n a n d 0 w i t h d i f f e r e n t s t a t i c s t r e s s r a t i o s a n d c y c l i c s t r e s s r a t i o s f8 ，d T s J T 8 ，c y s d C d I z o T 8 r s f T 8 8 j C d 0 0 5 3 0 0 2 4 7 4 3 1 2 O1 5 3 2 7 2 0 0 3 9 2 0 5 9 4 9 4 0 7 02 0 3 4 8 2 0 0 5 2 3 0 1 5 8 8 4 1 1 5 2 O 6 3 7 4 0 0 3 6 3 1 0 2 3 3 6 0 48 7 3 3 6 8 3 0 0 3 0 5 0 4 2 3 0 O 0 4 0 4 1 5 9 4 0 2 7 8 6 0 0 3 3 4 0 3 5 9 1 0 0 6 03 3 4 9 8 4 0 0 3 0 0 0 7 03 0 1 0 3 61 3 8 2 5 8 4 0 0 2 4 2 1 1 8 8 9 5 0 3 3 3 21 l 9 0 4 0 至此 ， 软黏土不 固结不排水条件下循环本构关系及模型参数的确定方法 已完全给出。 2 本构模型二次开发 2 1 平衡迭代过程 增量形式 的有 限元格式 中需要根据本步或迭 代计 算得到的位移增量 ， 决定 s， A o - ， 进而确定本步或迭代 结束时的弹塑性状态。 非线性平衡 方程采用 N e w t o n R a p h s o n迭代算 法 求解 。在迭代过程中， 把整个加载过程分成若干分段 ， 荷 载以增量形式施加在结构上， 在增量 内部进行若干迭代 ， 直至达到所需精度为止。整个流程融合 了增量与迭代 的 过程 。在增量内的迭代过程( 见图 3 ) 作如下说明 ： 根据当前的应力状态 ， 确定弹塑性矩阵D ， 进而确 定刚度矩阵 。 由 A u=K A R ， 得 到位移增量 。其 R 司 l R 图 3 增量 内迭代过程 Fi g 3 I t e r a t i v e p r o c e s s i n i nc r eme n t 中， A R为外加荷载增量 ， Au为位移增量。由 A e=B Au， 得到应变增量。采用应力更新算法 ， 通过当前 2 1 O 9 8 7 6 5 4 3 2 l O 2 1 O 9 8 7 6 5 4 3 2 l 0 学兔兔 w w w .x u e t u t u .c o m 8 4 水 利 水 运 工 程 学 报 2 0 1 5 年 6月 应力状态 以及上一步得到的应变增量计算应力增量 A o， 进而得到新的应力状态 ” 。判断 内力 、 外力是 否平衡 ， 位移是否满足收敛标准 ： 若不满足， 返 回 ， 继续进行迭代 ； 若满足， 进入下一增量计算。 2 2 应力更新算法 对第步的应力更新算法做出详细阐述， 由应力应变关系式( 2 ) ， 可得 ： d s =2 G d e 一1 5 ( 2 GH)s S k l d e 由上述关 系式可得到偏应力增量。静水压力增量可由弹性虎克定律得到 ： d o - =3 K d 6 由此得到应力增量大小 ： d o - =d o- +d 将式 ( 8 ) 和( 9 ) 代入式( 1 0 ) ， 有 ： ( 1 2 ) ( 1 3 ) ( 1 4 ) d o - =3 K 8 6 +2 G d e 一1 5 ( 2 G一日) s A s d e w ( 1 5 ) 得到增量形式的弹塑性应力应变关系 ， 整理成矩阵形式 ： d o =D印 d e ( 1 6 ) 式 中： d 为应力张量 ； d e为应变张量 ； D 为弹塑性矩 阵， 由于在应力应变关系 曲线上 D叩 为切线方 向， 因此 也 叫切线矩阵。 因此， 本构关系的积分式为： A o=I D 印 d e ( 1 7 ) J 在计算过程中， 采用式( 1 6 ) 将式 ( 1 7 ) 线性化 。但是这种计算方法必然带来偏差 。可采取减小增量大小 的方法减小这种偏差。 王勖成等 提出子增量方法将总的应变增量分成若干个子增量。每一个子增量结束时的弹塑性状态 作为下一个子增量的初始状态 ， 所有子增量计算结束 时的弹塑性状态作 为本步增量 的最终状态 。虽然子增 量法提高了计算精度 ， 加快了后继迭代 的收敛速度 ， 但是子增量数 目的增加使得计算量也相应增加。因此确 定适当子增量数 目很重要。本程序根据增量步开始 时的应力偏量 s 和增量步结束时的弹性预测应力偏量 s 之间的夹角 确定子增量数 ： N=1+ ， 其 中： k可由计算精度确定 ， 这里 k=0 O 1 。 当应力状态趋于边界面附近时 ， 由式 ( 5 ) 得硬化模量趋于零 ， 由式 ( 1 7 ) 得偏应力增量趋于零 ， 因此应力 状态点将始终保持在边界面以内。但实际计算 中， 由于采用式( 1 6 ) 所示的欧拉切线方法对式 ( 1 7 ) 进行线性 化 ， 而欧拉方法只有在偏应变增量趋于无限小的时候才能足够准确 ， 因此在边界 面附近时， 应力增量 d o - 有 可能大于零 ， 从而使应力状态超 出边界面。但边界面是破坏面 ， 土体 的应力状态不允许超过边界面。这将导 致计算出错。采用子增量方法虽然大大减小了偏应变增量 的大小 ， 但是也不足以保证应力状态不超过边 界 面。为此 ， 将超出边界面的点沿边界面的半径方向径向拉回至边界面 ， 并认为这一点的应力状态即为当前的 应力状态 。具体做法是 ， 计算应力空间原点与当前应力点之间距离 ， 按照边界面半径与该距离 比例对当前应 力点各个应力分量进行线性缩放 ， 则缩放后的应力点必然落在边界面上。进一步利用修正后 的应力状态建 立 与之 相适 应 的弹塑 性矩 阵 D 。这 样 就保 证 了任 意 时刻 的应 力 状 态 只能 在 边 界 面 以 内 或 者 在边 界 面 之 上 ； 而到达边界面的点在持续加载的条件下 ， 应力状态始终保持在边界面上 ， 类似 于理想弹塑性模型 中的塑 性加载状态 。与之相适应的弹塑性矩阵， 则保证了应力与应变状态 的一致性 ， 有利于后续分析的收敛 。 与单调加载相比， 循环荷载作用下土体的应力应变关系 ， 应力状态反向判断标准 的确定很重要。经过应 力状态的反 向， 土体单元的应力状态由加载变为卸载 ， 或者由卸载变为加载。本模型中应力状态反向对硬化 模量的大小演化有重要影响 ： 当应力反 向时， 前点变为映射 中心 ， 硬化模量变为最大值。有 了正确的应力状 态反向的确定标准 ， 才能保证后续计算 的准确性 。轴对称三轴应力状态下 ， 轴向应力方 向的改变 即是应力状 态 的反向。对于一般应力状态 ， 采用应力增量 s 与当前边界面外法线方 向 s 在偏应力空间中夹角 0的余 弦值来判断应力状态的改变 。 ( 1 8 ) 如果这一余弦值为正或零 ， 说明应力增量张量与当前边界面外法线张量呈锐角或直角， 则此时没有应力 学兔兔 w w w .x u e t u t u .c o m 第 3期 王哲学，等： 软黏土不排水循环应力应变关系的数值模拟 8 5 状态改变； 如果这一余弦值为负， 说明应力增量张量与当 前边界面外法线张量呈钝角 ， 则认为此时应力状态发生 改变。三轴应力状态 只有夹角为 0 。 ( 应力状态不变) 和 1 8 0 。 ( 应力状态反向) 两种应力状态的改变形式 。 3 有 限元程序验证 把有限元程序和单元预测循环三轴不 固结不排水试 验得到的应力应变 曲线进行对 比( 图 4 ) 。实线 为有限元 计算结果 ， 虚线为单元预测结果 。围压为 1 0 0 k P a ， 静强 度大小为 丁 8 ， =1 6 k P a ， 静剪应力 比 T 8 8 ， = ， 循环 ,a z 0 5 剪应力 比T 8 7 - 8 ， = ， 其他参数为 ： A o = ，G = ,a 0 3 9 2 31 k Pa 7 45 7 6 k Pa，C=0 8 61，d=0 6 0 6， o=8。 0 1 2 3 4 5 6 y 预测结果 测结果 图 4 有 限元预测结果 与单元预测结果 比较 F i g 4 C o mp a r i s o n b e t we e n F E M r e s u l t s a n d u n i t p r e d i c t e d r es u l t s 图 4表明数值算法结果与基于此本构模型进行 的单元预测结果吻合 良好 ， 说 明数值算法程序能够正确 反映本构模型的应力应变关系， 此程序编写是合理的。 4 三轴试验模拟 针对前文所述三轴试验 ， 静应力比 T 8 = 的三轴 压缩试验和三轴拉伸试验分别进行数值计, a 算 z 。根 0 据 5 表 1 ， 通 过插值的方法确定模型所采用 的参数( 表 2 ) 。 将循环八面体累积剪应变随循环次数的变化关系进行 比 较 ( 图 5 ) 。对于图5 ( a ) 中 f z ， =0 3 3 4的试验模拟 ， 当循环 应变大于 1 0 时 ， 计算得 到的曲线与试验 曲线 略有差 异 ， 但 表 2 三轴试验模型参数 Ta b 2 Mo de l pa r a me t e r s o f t r i a xi a l t e s t 循环累积应变小于 1 0 时 ， 二者吻合 良好。如上所述 ， 当静应变与循环累积应变之和达到 1 4 时 ， 即视为土 体破坏。这里静应变约为 4 ， 因此循环累积应变达到 1 0 时， 试样 已经破坏 ， 即在土样破坏之前可 以较好 地模拟试验结果。图 5 ( b ) 中 r =0 1 8 0的试验模拟 ， 循环次数达 9 0 0次 ， 模拟结果与试验有较好的一致 性 ， 且 曲线趋于平稳 ， 因此该模型可以描述循环次数较大时土体的应力应变响应 。综上 ， 相同振次下的累积 变形预测结果与试验结果较为一致 ， 此本构模型的有限元程序可以较好地反应循环变形的累积特性。 1 o o 2 0 HD 3 0 0 4 0 0 5 o o 6 o o 7 0 0 8 0 0 9 0 o ， 次 ( b ) 三轴拉伸 图 5 循环 累积剪应变 随振次 的变 化 F i g 5 C y c l i c c u mu l a t i v e s h e a r s t r a i n V S t h e n u mb e r o f c y c l e s 将静剪应力比 T s , 7 s I， =0 5 、 循环剪应力比 r s x y T s ， =0 3 9 2的应力应变 曲线进行 比较 ( 图6 ) 。结果表 明， 模型可较好地描述循环加载时应力应变曲线的滞 回特性及应变累积特性。 6 4 2 O 8 6 4 2 21O98 765432l 0 学兔兔 w w w .x u e t u t u .c o m 学兔兔 w w w .x u e t u t u .c o m 第 3期 王哲学 ， 等 ： 软黏土不排水循环应力应变 关系的数值模拟 8 7 A n a l y s i s a n d c o m p a ri s o n o f s e v e r a l c y c l i c p l a s t i c m o d e l s J N i n g x i a E n g i n e e r i n g T e c h n o l o g y ，2 0 0 3 ，2( 1 ) ：1 9 2 3 ( i n C h i n e s e ) ) 程 星磊 ， 王建华 ， 李书兆 软黏土 不排 水循 环应 力应 变 响应 的弹塑 性模 拟 J 岩 土工 程学 报 ， 2 0 1 4 ，3 6 ( 5 ) ： 9 3 3 9 4 1 ( C H E N G X i n g l e i ， WA N G J i a n h u a ， L I S h u z h a o E l a s t o p l a s t i c s i mu l a t i o n f o r u n d r a i n e d c y c l i c s t r e s s s t r a i n r e s p o n s e o f s o f t c l a y J C h i n e s e J o u rna l o f G e o t e c h n i c a l E n g i n e e ri n g ， 2 0 1 4 ， 3 6 ( 5 ) ： 9 3 3 9 4 1 ( i n C h i n e s e ) ) 1 1 张学言， 闫澍旺岩土塑性力学基础 M 天津： 天津大学 出版社， 2 0 0 3 ：8 8 1 1 5 ( Z H A N G X u e y a n ， Y A N S h u w a n g F u n d a me n t a l s o f g e o t e c h n i c s p l a s t i c i t y M T i a n j i n ： T i a n j i n U n i v e r s i t y P r e s s ， 2 0 0 3 ： 8 8 1 1 5 ( i n C h i n e s e ) ) 1 2 王军祥， 姜谙男 完全隐式返回映射算法对岩土地基问题的求解 J 工程力学 ， 2 0 1 3 ， 3 6 ( 8 ) ： 8 3 8 9 ( WA N G J u n x i a n g ， J I A N G A n n a n F u l l y i m p l i c i t r e t u r n m a p p i n g a l g r i t h m fo r s o l v i n g t h e p r o b l e m s o f g e o t e c h n i c a l fo u n d a t i o n J E n g i n e e r i n g Ma c h a n i c s ， 2 0 1 3 ，3 6 ( 8 ) ： 8 3 8 9 ( i n C h i n e s e ) ) 1 3 朱伯芳 有限单元法原理与应用 M 北京 ： 中国水利水 电出版社，1 9 9 8 ： 2 9 2 3 0 2 ( Z H U B o f a n g T h e fi n i t e e l e m e n t m e t h o d t h e o r y a n d a p p l i c a t i o n s M B e i j i n g ： C h i n a Wa t e r P o w e r P r e s s ， 1 9 9 8 ： 2 9 2 3 0 2 ( i n C h i n e s e ) ) 1 4 王勖成 ， 邵敏 有 限单元法基本原理和数值方 法 M 北京 ： 清华 大学 出版社 ， 1 9 9 7 ： 5 1 0 5 1 1 ( WA N G X u c h e n g ，S H A O Mi n T h e b a s i c p ri n c i p l e s o f t h e f i n i t e e l e m e n t m e t h o d a n d n u me ri c a l me t h o d s M B e i j i n g ： T s i n g h u a U n i v e r s i t y P r e s s ，1 9 9 7 ： 5 1 0 5 1 1 ( i n C h i n e s e ) ) HU Cu n，LI U Ha i x i a oI mpl i c i t a n d e x pl i c i t i n t e g r a t i o n s c h e me s i n t he a n i s o t r o pi c bo u n di n g s ur f a c e p l a s t i c i t y mo d e l f o r c y c l i c b e h a v i o u r s o f s a t u r a t e d c l a y J C o m p u t e r s a n d G e o t e c h n i c s ， 2 0 1 4 ，5 5 ： 2 7 4 1 Nume r i c a l s i mul a t i o n o f u nd r a i ne d c y c l i c s t r e s s - s t r a i n r e s po n s e o f s o f t c l a y W ANG Z h e x u e 一，WANG J i a n h u a ，C HENG Xi n g 1 e i ， ( 1 S t a t e K e y L a b o r a t o r y o f H y d r a u l i c E n g i n e e r i n g S i m u l a t i o n a n d S a f e t y ，T i a n j i n U n i v e r s i t y ，T i a n j i n 3 0 0 0 7 2 ， C h i n a ： 2 G e o t e c h n i c a l E n g i nee r i n g I n s t i t u t e ， T i a n j i n U n i v e r s i t y ，T i a n j i n 3 0 0 0 7 2 ，C h i n a ) Abs t r a c t：An i n c r e me n t a l e l a s t o - p l a s t i c c o n s t i t u t i v e mo de l i n t he f o r m o f t o t a l s t r e s s de v e l o p e d b y CHENG Xi n g - l e i e t a l c a n e f f e c t i v e l y r e fle c t t h e u nd r a i ne d c y c l i c s t r e s s s t r a i n r e s p o n s e o f t h e s a t u r a t e d s o ft c l a y I n o r d e r t o s i mu l a t e t h e r e s p o ns e o f t he s o ft c l a y u n de r c o mp l e x s t r e s s s t a t e s， a s e c o n d a r y d e v e l o p me n t for t h e c o n s t i t u t i v e mo d e l i s c o n d u c t e d b y us e o f t h e fin i t e e l e me n t me t h o dTh e Ne wt o n Ra p hs o n me t h o d i s a p pl i e d t o s o l v e i t e r a t i v e l y t he ma t e ria l n o n l i n e a r pr o b l e m ； c o n s i d e rin g a s e rie s o f k e y p r o b l e ms s u c h a s s t r e s s r e v e r s i n g， i n p r o g r a mmi n g t he mo de l ，Eu l e r t a n g e nt a l g o r i t h m i s a d o p t e d， wi t h a s u b i n c r e me n t a l me t ho d i mp r o v i n g c a l c u l a t i o n a c c u r a c y Th e r e s p o ns e o f t h e s o ft c l a y u n d e r a x i s y mme t r i c s t r e s s s t a t e s i s p r e d i c t e d，a n d t h e s t r e s s s t r a i n r e l a t i o n s h i p C HI V E S a r e r e c e i v e d，whi c h i s c o mp a r e d wi t h t h e r e s u hs g i v e n b y t h e u n i t p r e d i c t i ng The a n a l y s i s c o mpa ris o n r e s u l t s s h o w t h a t i t i s i n g o o d a g r e e me n t ，v e rif y i n g t h e r a t i o n a l i t y o f t h e fi ni t e e l e me n t p r o g r a mTh e fini t e e l e me n t p r o g r a m i s a p p l i e d t o s i mu l a t e u n c o n s o l i d a t e d u n d r a i n e d t r i a x i a l t e s t s，a nd t h e r e s u l t s a r e i n g o o d a g r e e me n t wi t h t h e e x p e r i me n t s ， wh i c h i nd i c a t e s t h a t t he fi ni t e e l e me n t p r o g r a m o f t h e c o n s t i t u t i v e mo d e l c a n d e s c r i b e t h e un d r a i n e d s t r e s s s t r a i n b e h a v i o r o f t h e s o f t c l a y u n d e r t r i a x i a l s t r e s s s t a t e s Fu r t h e r mo r e，i t c a n b e a p p l i e d t o s o l v e mo r e c o mpl e x b o u n d a r y v a l u e pr o b l e ms Ke y wo r d s：c y c l i c l o a d；s a t u r a t e d s o ft c l a y；s t r e s s s t r a i n c u r v e； a c o n s t i t u t i v e mo d e l ；Eu l e r t a n g e n t a l g o rit h m ； t r i a x i a l t e s t s 学兔兔 w w w .x u e t u t u .c o m