超固结土塑性次加载面模型的数值实施及应用

黄 杰，李持庆

(1.河海大学力学与材料学院，江苏南京 210098；2.山东省城乡建设勘察院，山东 济南250031)

0 引 言

近年来，随着大型商业有限元软件的发展和完善，美国HKS公司的软件产品ABAQUS[1-4]，它以其强大的分析能力和高度的非线性问题模拟能力而受到各界研究机构的广泛应用。

目前，ABAQUS有限元计算软件虽然包含了包括土木、机械、电子、材料等各个专业的多种本构模型，但是在我国大部分地区，尤其是东南沿海一带，都处在软土地基上。软土的受力特性非常复杂，不能简单认为某个土体单元加载或卸载，这就必然考虑土的固结等非线性问题。姚仰平等[5-7]在考虑土的复杂力学特性基础上，提出了可以反映土的剪胀、剪缩及硬化软化特性的塑性次加载面模型(当前屈服面和参考屈服面)。本文借助ABAQUS有限元软件提供的用户子程序UMAT接口，实现了此模型的二次开发并对工程算例进行了测试。计算结果反映出了较好的稳定性和较高的精度，符合实际情况。

1 基于改进伏斯列夫线的超固结土本构模型

塑性次加载面模型是在剑桥模型的基础上，通过引入统一硬化参数H而提出的。此模型弥补了剑桥模型不能反映土体剪胀性的问题[6-7]。

当前屈服面函数:

参考屈服面函数:

统一硬化参数:

超固结参数:

潜在强度:

从此本构模型中，我们可以看出超固结参数、潜在强度、硬化参数之间相互影响相互制约。也就是说，模型通过潜在强度，将硬化参数和超固结参数之间建立了一定的联系。当超固结土受一定荷载时，当前屈服面与参考屈服面共同决定了超固结参数的变化，进而影响了潜在强度及硬化参数。

2 模型计算流程

在程序的每一步开始时，ABAQUS都会自动调用UMAT子程序[8]，并且传入相应的时间增量、荷载增量和应变增量，再通过给定的状态变量计算相应的应力增量和jacobian矩阵，最终形成整体刚度矩阵[9]。其流程如图1所示。

3 算 例

考虑到地下土体开挖过程中的对称性，选用1/2土体进行分析计算，土体右方和下方采用无限区域，模型计算范围选取地表以下50 m，横向范围30 m，隧洞直径8 m，隧洞中心距离地表20 m，衬砌厚度0.15 m。

计算模型分为两个部分:周围土体部分和衬砌部分。周围土体部分，划分为177个单元，采用平面应变单元类型，即CPE4；衬砌部分，划分为12个单元，采用非协调模式平面应变单元，即CPE4I。模型如图2所示。

图1 模型计算流程

图2 隧洞开挖计算模型

模型计算参数见表1。

表1 饱和粘土材料参数

模型加载采用分步分级加载的方式，按照应力释放法逐步降低衬砌与土体单元相接触部分结点的反力，达到应力释放的效果。

自然状态下，土体材料会呈现一定的超固结状态。随着土体深度的增加，超固结特性逐渐减弱，最终达到OCR=1的正常固结状态。在图3和图4中，我们分别分析了土体在正常固结状态下，进行深隧洞开挖和支护过程的水平位移和竖直沉降结果。从实际角度考虑，模拟计算的沉降结果都偏大，且过于理想化，无法满足实际工程土体都存在一定超固结性的要求，接下来我们按照土体所具有的不同超固结比来进一步分析，并做相应比较。

图3 各结点水平位移云图

图4 各结点竖直位移云图

五种不同超固结状态的土体进行隧洞开挖结束后，表层土沉降均呈现平滑曲线形式，如图5所示，随着距离隧洞中心水平距离的增大，沉降量逐渐减少，达到一定距离时，土体出现隆起的现象。且随着超固结性越高(即OCR值越大)，表层土体最大沉降值越小，也越早出现土体的隆起上涌现象。对比曲线如图5所示。

图5 表层土各结点竖直位移对比

隧洞正上方土体竖直位移分布如图6所示。由图6可看出，随着深度的增加，沉降量逐步累积，在隧洞上方沉降量达到最大值。随着超固结性越高(即OCR值越大)，隧洞正上方土体累计沉降量越小。

图6 隧洞正上方土各结点竖直位移对比

隧洞上方表层土体水平位移分布如图7所示。由图7可看出，随着表层土体距隧洞中心水平距离的增大，水平位移呈现负向先增后减的状态。在距离11.2 m附近，负向水平位移达到最大值，当距离增大到30 m以上时，水平位移减到非常小，可以忽略不计。考虑土体所具有的不同超固结性的影响，随着超固结比的增大，负向水平位移值逐步减小[10]。不同超固结状态土体水平位移对比如图7。

距离隧洞中心18.52 m处土体水平位移如图8所示。由图8可看出，最大负向位移在地表处，随着深度的增大，负向位移逐步减小，并最终变为正向位移，在隧洞中心右侧时，正向位移达到最大值。随着深度的进一步增大，位移逐步减小并趋于0。

图7 不同超固结度表层土各结点水平位移对比

图8 不同超固结度隧洞两侧18.52m处土各结点水平位移对比

4 结 语

经典的剑桥模型，其屈服和破坏均是在基于准则的基础上，即假设屈服面在平面上的投影为圆形。这大大高估了土体的三轴拉伸强度，与实际数据存在较大差异。本文借助ABAQUS平台所开发的UMAT子程序，将超固结比引入模型中，解决了此类问题。

当然，一个本构模型的提出离定量地解决实际工程还有较大的差距。塑性次加载面模型虽然能够近似反映土体的强度和变形特性，但由于土体复杂的特性，还需要进一步改善。此子程序仅考虑土体基本特性，但未引入地下水渗流造成的影响，进一步的工作还需要继续进行。

[1]Hibbitt，Karlsson，Sorenson，INC.ABAQUS/Standard:User'sManual[M].America:Hibbitt，Karlasscn&Screnson，2002.

[2]Hibbitt，Karlsson，Sorenson，INC.ABAQUS/Analysis User's Manual[M].America:Pawtucket，RI，2000.

[3]石亦平，周玉蓉.ABAQUS有限元分析实例详解[M].北京:机械工业出版社，2006.

[4]王金昌，陈页开.ABAQUS在土木工程中的应用[M].杭州:浙江大学出版社，2006.

[5]谢定义，姚仰平，党发宁.高等土力学[M].北京:高等教育出版社，2008.

[6]姚仰平，李自强，侯 伟，等.基于改进伏斯列夫线的超固结土本构模型[J].水利学报，2008，39(11):1244-1250.

[7]姚仰平，侯 伟，周安楠.基于Hvorslev面的超固结土本构模型[J].中国科学(E 辑)，2007，37(11):1417-1429.

[8]Lee K，Ghosh S.Small deformation multi-scale analysis of heterogenous materialswith the Voronoi cell finite model and homogenization theory[J].Computational Materials Science，1996，7(1-2):131-146.

[9]白 云.FORTR AN90程序设计[M].上海:华东理工大学出版社，2003.

[10]Wang X，Wang L B，Xu L M.Formulaton of the return mapping algorithm for elastoplastic soilmodels[J].Computers and Geotechnics，2004 ，31(4):315-338.