600 岩 土 力 学 2011年

足的缺点，也拓展了ABAQUS的适应范围，缩短程序开发周期，极大降低了开发和维护的成本，是解决复杂土工数值分析的有效途径[7

－8]

ε(t)+ ε(t)。

④主程序进行迭代次数的检查，如超过限定的迭代次数，主程序自动减小时间增量并返回，重新迭代计算，然后主程序将σ(t+ t)代入系统平衡方程，进行平衡判断，若满足归一化系统平衡方程，结束本次迭代计算，进入下一步计算，否则，主程序放弃更新，重新回到增量步初始时刻的值，并进行本次增量步的下一步平衡迭代，直至满足系统平衡方程。

。徐远杰等[9]

完成了Duncan-Chang本构模型在ABAQUS中二次开发；岑威钧等[10]完成了土石料本构模型二次开发。由此可见，基于ABAQUS的二次开发是目前土工数值模拟的发展趋势，因而有必要把目前岩土工程界常用的土体本构模型嵌入到ABAQUS原程序之中。

本文利用ABAQUS提供的二次开发用户子程序，编制了模型库接口程序，完成了土体本构模型库的开发研究。实际工程应用表明，程序运行稳定，模拟结果更接近实际工程。

3 本构模型库

借助ABAQUS强大非线性计算能力，本文完成了目前国内常用土体本构模型库的开发，主要包括在国内岩土工程领域广泛采用的一些土体本构，如Duncan-Chang弹性非线性模型、南水弹塑性模型及基于状态相关的砂土本构模型。这些本构模型有的在实际工程中得到了广泛的应用，积累了丰富的实践经验，有的很好地反映相关土体的一些特征，如砂土的剪胀性等。 3.1 Duncan-Chang本构模型

Duncan-Chang等根据康纳（Kondner）的建议，将三轴剪切试验中当σ3等于常数时的(σ1 σ3)-

2 用户子程序及算法

ABAQUS主程序包括42个用户子程序接口和13个应用程序接口，可以定义边界条件、载荷条件、接触条件、材料特性及利用子程序与其他软件进行数据交换等。

用户材料子程序（user-defined material mech- anical behavior），简称umat，通过与主求解程序的接口实现与ABAQUS的数据交流，在输入文件中，使用关键词*umat，表示自定义用户材料。在umat中，可以实现定义材料的本构关系，也可以应用依赖于解的状态变量及结合应用程序来定义场变量。

ABAQUS计算程序与umat子程序之间是一个动态交互传递数据、协同工作的过程。umat子程序作为计算主程序的一个接口，在单元积分点上调用，增量步开始时，计算程序通过接口进入子程序，单元积分点必要变量的初始值随即传递给子程序相应变量，umat子程序计算单元积分点的雅克比矩阵，并更新应力张量和状态变量，最后将这些变量的更新值通过接口返回计算程序。

本文基于广义中心积分算法，通过平衡迭代保证应力状态始终位于当前屈服面，因而具有较高的计算精度和较快收敛速度[11

－13]

ε1关系和ε3-ε1关系近似地用双曲线来表示为

ε1

σ1 σ3= （1）

a+bε1

ε3=

hε1

（2） 1 dε1

在上述2个基本假定的基础上，Duncan-Chang

E-B模型的表达式为

R(1 sin )(σ1 σ3) σ3 Et= 1 f Kpa

2ccos +2σ3sin pa σ dp

=Kbpa 3 Bt=dεv pa (σ1 σ3)f=Rf=

m

2

n

。计算过程如下：

2ccos +2σ3sin

1 sin

①平衡时刻t，主程序提供给子程序总应变

ε(t)和总应力σ(t)，同时自动生成一个时间增量

t。

②abaqus主程序提供一个总应变增量 ε(t)，调用子程序，计算相应的应力增量 σ(t)、为非线性迭代求解过程，迭代收敛后更新应力（应变）张量σ(t+ t)=σ(t)+ σ(t)。

③子程序将更新后的应力（应变）值返回给计算主程序，同时自动更新总的应变ε(t+ t)=

(σ1 σ3)f(σ1 σ3)ult

（3）

卸荷和重复加载时的回弹模量采用下式表示：

σ

Eur=Kurpa 3 （4）

pa

式中：c为黏聚力； 为内摩擦角；(σ1 σ3)ult为剪应力极限值；(σ1 σ3)f为破坏剪应力；Kur为回弹模量指数，一般取Kur=(1.5～3)K，K及n为切线

n