- 1 -
中国科技论文在线
三维超固结土弹粘塑性本构模型#
姚仰平，孔令明*
基金项目：高等学校博士学科点专项科研基金（编号 20091102110030
作者简介：姚仰平（1960-），男，教授，主要研究方向：土的本构关系
（北京航空航天大学，交通科学与工程学院）
5 摘要：分析蠕变对土超固结程度的影响，基于蠕变规律，建立老化时间与超固结参数的关
系；采用 UH 模型的再加载公式计算瞬时变形，提出超固结土的一维弹粘塑性应力应变关
系；分析一维弹粘塑性应力应变关系的演化规律，定义特征速率并给出其与超固结参数的
关系；在 UH 模型屈服面方程中引入折算时间，建立三维弹粘塑性本构模型。该模型在 UH
模型的基础上考虑时间效应，不仅能够描述土的剪切蠕变、应变率效应等粘性规律，而且10
能够反映土的超固结性；与修正剑桥模型相比，仅增加一个参数来反映土的蠕变规律；该
模型在超固结土瞬时加载时退化为原始 UH 模型。
关键词：蠕变；粘塑性；应力应变；超固结；三维
中图分类号：TU43
15
A three-dimensional elastic-viscoplastic constitutive model
for overconsolidated clays
YAO Yangping, KONG Lingming
(School of Transportation Science and Engineering , Beihang University)
Abstract: The influences of creep on clays are discussed first. Based on creep laws, the relation 20
between aging time and overconsolidation parameter is established. With the reloading function of
UH model (Unified Hardening model) which is used to calculate the instant compression, a
one-dimensional elastic viscoplastic stress-strain relationship is proposed. Furthermore, evolution
of the relationship is analyzed, and critical rate which is a function of overconsolidation parameter
is defined. Then a three-dimensional elastic viscoplastic model is proposed by incorporating 25
equivalent time into the yield function of UH model. The new model considering time effects on
the basis of UH model, can describe not only shear creep, rate effect and other viscous phenomena,
but also characteristics of overconsolidated clays. Compared with the modified Cam-clay model,
the proposed model requires only one additional parameter to reflect creep rules. The model would
be changed into the original UH model under instant loading. 30
Key words: Creep; Viscoplasticity; Stress-strain; Overconsolidation; Three-dimension
0 引言
土应力应变关系的时间相关性被称作粘性。与土粘性相关的现在主要包括蠕变和应变率
效应。它们能够改变土的工程性质从而造成地基沉降和滑坡等工程问题。因此，人们对粘性35
进行了大量研究并建立了许多本构模型来描述土的粘性。
目前建立考虑粘性的本构模型的方法主要有三种：经验方法、粘弹性方法、粘塑性方法。
以经验方法建立的模型在使用时受到边界和加载条件的限制，而以粘弹性方法建立的线性模
型不适用于土发生大应变的情况。因此，越来越多的研究集中在以粘塑性方法建立本构模型
上[1-3]。然而，迄今为止以粘塑性方法建立的考虑粘性的本构模型大多只能描述正常固结土40
或轻微超固结土，并不适用于描述超固结程度较大的土。尽管部分模型做出改进以反映超固
结土的粘性规律，但这些模型公式复杂、参数较多[1]。所以，有必要建立一个新的、简单、
实用的粘塑性模型来描述超固结程度较大的土在时间影响下的应力应变关系。
- 2 -
中国科技论文在线
UH 模型[4]是采用统一硬化参数建立的三维超固结土弹塑性本构模型，能在没有比修正
剑桥模型增加任何参数的基础上能够描述土体的剪缩、剪胀、硬化、软化和应力路径依赖性。45
本文拟在 UH 模型基础上，引入时间效应参量，建立合理、简单、实用的弹粘塑性模型以同
时描述土的超固结特性与粘性规律。本文模型的粘性部分基于蠕变规律建立，非粘性部分基
于 UH 模型建立。蠕变规律的描述只需一个描述蠕变难易程度的参数，蠕变规律与应力应变
耦合后可以反映土的加载率效应。
1 一维弹粘塑性应力应变关系 50
固结过程
固结过程在传统上被划分为主固结与次固结。然而，对渗透性较小的软粘土而言，其排
水速率过慢导致其主固结过程较长，主固结中也伴随有不可忽视的蠕变产生。因此，Bjerrum[5]
提出应将固结过程划分为瞬时压缩、延时压缩，如图 1 所示。
瞬时正常压缩线
正常固结线
e
A
B
F
C
C'
B'
1 年
10 年
100 年
lnp O
55
图 1 固结过程的划分
Fig. 1 Division of consolidation progress
由于土骨架的压缩过程中不可避免的存在孔隙水的排出，因此，针对饱和粘土并不存在
瞬时压缩线，一般固结试验中所得的正常固结线是土骨架经历了孔隙水压力消散时间的压缩
曲线。不过，Bjerrum 通过试验证明了瞬时压缩线在理论上存在的正确性。如图 1 所示，从60
正常固结线上的不同点进行等时间蠕变，其最终状态的连线为平行于正常固结线的直线，即
BC=B'C'。随着蠕变时间的变化，会形成一簇平行于正常固结线的直线，如“1 年”、“10
年”、“100 年”所示。依次类推，则必然存在一条蠕变时间为“0”的直线，该线即为瞬
时压缩线。这条瞬时压缩线与正常固结线平行，表明土表现出正常固结土的特征，本文称之
为瞬时正常压缩线。 65
根据瞬时、延时压缩理论，土的压缩变形可表示为：
s t se sp tp e p
v v v v v v v v
(1)
式中， sv 为由应力作用产生的应变，即瞬时压缩，包括应力产生的弹性应变
se
v 和塑性应变
sp
v ；
tp
v 为时间作用产生的应变，即延时压缩，时间作用只产生塑性应变，故
t tp
v v ；
e
v 为
弹性应变，只由应力产生； pv 为塑性应变，包括应力塑性应变
sp
v 和粘塑性应变
tp
v 。 70
蠕变规律
根据瞬时、延时压缩理论，只需知道瞬时压缩 AB 及延时压缩 BC 的计算公式即可得土
的总压缩变形。
- 3 -
中国科技论文在线
瞬时压缩
正常固结土的瞬时压缩线与正常固结线平行，在 e-lnp 空间内可用公式(2)表示： 75
ln
p
e e
p
(2)
公式(2)中各变量含义如图 2 所示，λ 为瞬时正常压缩线的斜率，p 为平均正应力，pλ0为初始
平均正应力，eλ 为瞬时正常压缩线上当前应力 p 所对应的孔隙比，eλ0 为瞬时正常压缩线上
初始应力 pλ0所对应的初始孔隙比。图 2 中的其余变量，e 为孔隙比，κ 为弹性回弹线斜率，
pc为前期固结压力。 80
瞬时压缩
图 2 中 FC 是有效应力恒定的蠕变段，通过固结试验可得 FC 的变形公式，即为蠕变公
式。图 3 中点线为固结试验数据[6]。分析可知，在 e-ln(t+1)空间内固结曲线的 FC 近似直线，
因此蠕变公式可采用对数形式，如下:
β ln( 1)e e t (3) 85
公式(3)中，t 为蠕变时间，β 为 e-ln(t+1)坐标系中蠕变曲线的斜率，是一个粘性参数，反映
了土蠕变的难易程度。
瞬时正常压缩线
瞬时压缩
延时压缩
A
B
D
F
C
λ
a
1
1
κ
lnp pc p pλ0 O
e
eλ
eλ0
e
图 2 正常固结土压缩曲线
Fig. 2 Compression line of normally consolidated clays 90
1 10 100 1000 10000 100000
延时压缩
瞬时压缩
孔隙水消散 有效应力恒定
A
B
F
C β 1
e
ln(t+1): min
图 3 固结试验曲线
Fig. 3 Consolidation curve obtained from test
- 4 -
中国科技论文在线
老化效应 95
图 2 中，BC 描述了土的蠕变过程。土经过蠕变后由正常固结状态进入超固结状态。可
知，时间对土的作用相当于使土老化变硬。天然土沉积后，在自重的作用下经历长时间的蠕
变，其强度会有所增加，表现出一定的超固结性，体现了时间的老化效应。这种时间老化效
应所对应的时间可称为老化时间。因此，公式(3)中的蠕变时间可以用老化时间 ta 代替。土
的超固结度代表了土老化效应的强弱，通过图 2，可以建立老化时间与超固结参数的对应关100
系：
α
a
1t R (4)
公式(4)中，ta 为老化时间，α = (λ - κ) / β。R 是 UH 模型中描述土超固结程度的参数，指的
是等应力比条件下当前应力与参考应力的比值，将在 节中详述。在一维条件下，R 与
OCR（超固结度）呈倒数关系，即 R=1/OCR。 105
土的超固结程度表明了土所处的密实状态。而蠕变率仅取决于土所处的状态，与应力历
史与应力路径无关。因此，蠕变率取决于土的超固结参数 R。土的老化程度除受时间的影响
外，还受应力影响。例如，图 2 所示正常压缩再卸载路径 BDC，虽然没有经历时间，却使
土产生了超固结，对土的作用等同于 BC。故老化效应是一种广义效应，凡是使土具有超固
结性的效应均为老化效应。 110
一维应力应变时间的关系
时间老化效应使土产生超固结。因此，研究土的应力应变关系时间相关性不可避免的要
涉及土的超固结。UH 模型是目前能够合理描述土超固结性的、最简单的本构模型。本文模
型通过将蠕变效应引入 UH 模型建立，可以同时描述土的应力应变关系时间相关性和超固结
特性。 115
瞬时应变增量
正常固结土的瞬时变形采用瞬时正常压缩线计算。对于超固结土，本文采用 UH 模型的
再加载线计算，如图 4 所示。
O lnp p0
再加
载线
κ 线
瞬时正常压缩线 e
e
λ
e
0
e
B
D
E
C
0p
图 4 瞬时加载线 120
Fig. 4 Instant compression line
经典的弹塑性本构模型，如修正剑桥模型，认为处于超固结状态的土的加载曲线为弹性
压缩线，即 κ 线。试验则表明超固结土的加载仍然为弹塑性变形，加载线在 e-lnp 坐标系内
为如图 4 所示的曲线。该加载曲线的斜率介于 λ 与 κ 之间，当超固结度减小时，曲线斜率趋
近于 λ，反之则趋近于 κ。UH 模型基于以上思想建立了再加载线公式，其斜率如式(5)，其125
- 5 -
中国科技论文在线
实质是对加载线斜率在 λ 和 κ 之间进行插值，使其符合试验规律。
4
4
f
d M
= κ + (λ κ)
dln M
e
p
(5)
其中 M 为临界状态应力比，是正常固结土的强度参数，Mf是根据超固结程度对 M 进行修正
所得的超固结土潜在强度，采用公式(6)计算[7]。
fM 6 1
R R R
(6a) 130
2M
χ =
12(3 M)
(6b)
其中 R 的计算可见公式(16)。由公式(5)，超固结土瞬时应变增量的计算公式如下：
4
s
v 4
0 0 f
κ λ κ M d
d = +
1+ 1+ M
p
ε
e e p
(7)
公式(7)中，e0 为土的初始孔隙比。
延时应变增量 135
通过公式(4)可知，初始超固结土相对于初始正常固结土，经历的老化时间较长，两者
蠕变计算的区别在于初始老化效应的不同。由公式(3)和公式(4)得：
tp α
v
0 a 0
β d β
d = = d
1+ +1 1+
t
ε R t
e t e
(8)
公式(8)中 R 反映老化效应对蠕变速率的影响。R 愈小，老化效应愈强，粘塑性应变愈小，
表明时间对土的影响随着老化效应的增强而逐渐减弱的性质。 140
应力-应变-时间的关系
结合公式(1)、(7)、(8)得超固结土一维弹粘塑性应力-应变-时间关系：
4
α
v 4
0 0 f 0
κ d λ κ M d β
d = + + d
1+ 1+ M 1+
p p
ε R t
e p e p e
(9)
当瞬时加载，即没有时间效应时，公式(9)退化为 UH 模型的一维应力-应变关系。
应变率效应 145
图 5 中以不同应变率加载所得的压缩曲线最终均平行于正常压缩线。当应变率较大时，
如图5中的1×10-3/min速率的压缩线，土表现出超固结性。当应变率较小时，如图5中的2×10-5
速率的压缩线，处于超固结状态的土不仅没有出现超固结土的压缩曲线，反而表现出类似欠
固结的性质。由此可见土的超固结与否不仅受应力历史影响，更与应变率息息相关。所以，
必然存在有一个特定的应变率，当以此应变率加载时，土恰巧表现出正常固结土的性质，本150
文称这一特定的应变率为特征速率。特征速率与土的状态一一对应，土的每个状态点均有其
所对应的特征速率，具有相同特征速率的状态点的连线被称为特征速率线。如果当前对土所
施加应变率不等于其特征速率，则土的加载线会不断变化，直至其最后趋近于和当前应变率
等值的特征速率线。
- 6 -
中国科技论文在线
10 100 1000
∞
fast
slow
veryslow
瞬时
1×10-3
2×10-5
1×10-6
松弛
瞬时正常压缩线
应变率：1/min
蠕
变
变速率
特征速率
e
logp:kPa
155
图 5 不同应变率下的压缩曲线
Fig. 5 Compression curves at various strain rates
特征速率线上的土的超固结度是恒定的，即土的超固结度与特征速率是一一对应的。根
据公式(9)，可以推得超固结参数 R 与特征速率 v,crε 的关系式：
-1
4
α
v,cr 4
0 f
λ β M
= 1
1+ λ κ M
ε R
e
(10) 160
式中 Mf的计算请见公式(6)。
2 三维弹粘塑性模型
在 UH 模型基础上，结合土的蠕变规律，将时间效应参量引入屈服函数，采用相关联流
动法则，建立三维弹粘塑性本构模型，使其可以描述土在剪切作用下的超固结特性和粘性规
律。 165
当前屈服面
A
B
C
F
参考屈服面
xtpxtp
( , )p q
( , )p q
M
1
p
G
t
D
E
xpxpO
q
图 6 当前屈服面和参考屈服面
Fig. 6 Current yield surface and reference yield surface
屈服条件
当前屈服面 170
如图 6 所示，土在时间效应作用下的当前屈服面为当前应力点 ( , )p q 所在的屈服面。该
屈服面是通过在 UH 模型当前屈服面的基础上引入时间效应参量获得，采用统一硬化参数
H 作为硬化参量（H 将在 节介绍）。当前屈服面函数可表示为：
- 7 -
中国科技论文在线
2
2
xt0 p
1
ln ln 1 d 0
M c
p
f t H
p
(12a)
4
αf
4
Mβ
= d
λ κ M
t R t
(12b) 175
式中 pxt0 为当前屈服面与 p 轴的初始交点， cp=(λ−κ)/(1+e0)，η=q/p 为应力比， t 的表达式
如公式(12b)。 t 不是真实的时间，而是根据土的超固结状态修正所得的时间参量，代表了
时间对土的影响程度，本文定义 t 为折算时间。
参考屈服面
如图 6 所示，参考屈服面是参考应力点 ( , )p q 所在的屈服面，( , )p q 为 ( , )p q 在等应力180
比条件下在参考屈服面上的对应点。在 UH 模型中，当前屈服面相当于图 4 中的 C 点，代
表土的超固结状态，参考屈服面相当于 D 点，代表土的正常固结状态。引入参考屈服面后，
可借用当前屈服面与参考屈服面的关系来描述土的超固结程度。UH 模型采用修正剑桥模型
（正常固结土模型）的屈服面作为参考屈服面，本文模型也采用该屈服面，其函数如公式
(13)： 185
2
p
v2
xt0 p
1
ln ln(1 ) 0
M c
p
f
p
(13)
式中
xt 0p 为参考屈服面与 p 轴的初始交点。
图 6 显示了 p - q - t 空间的三维屈服面。空间曲面 ABC-OGF 内的空间为应力加载弹性
区，其外为应力加载弹塑性区。
超固结参数 190
当前屈服面与参考屈服面的关系反映土的超固结程度，这一关系即为超固结参数 R：
p
R
p
(14)
由(13)式，有：
p2
v
x0 2 2
p
M
exp
M c
p p
(15)
联立（14）、（15）式，可得： 195
p2
v
2
p
1 exp
M c
p
R
p
x0
(16)
p
v 是加载过程中产生的塑性体积应变，包括应力塑性体积应变和粘塑性体积应变。
硬化定律
姚仰平等[4]建议采用与潜在强度和临界状态应力比相关的硬化参数 H：
4 4
p pf
v v4 4
M 1
d d d
M
H H
(17) 200
式中， 4 4 4 4f(M ) / (M ) 。Mf在 节已介绍，在三维应力状态，Mf具有更明确的
物理意义，是指在当前密度以及时间效应、应力条件下，超固结土的潜在抵抗破坏能力。
- 8 -
中国科技论文在线
弹粘塑性本构关系
根据屈服函数，将时间增量作为等效荷载，在 p-q 空间中的应力应变关系可以表示为：
1 21
2 32
3
3 3
v
d
K A KG A ddp B dt
KG A G A ddq B dt
(18) 205
式中，
3(1 2 )
E
K
为弹性体积模量，
2(1 )
E
G
为弹性剪切模量，A1、A2、A3 为不同的塑
性影响因素，B1、B2 为不同的时间影响因素，可以推导得出 A1、A2、A3、B1、B2 的表达式：
4 4 2
f p
1 2
4 4 2 2 2
f p p
M 12 c
M 12 c c M
p G
A
p G K
(19a)
2 2
p
2 2
4 4 2 2 2
f p p
2c M
M 12 c c M
A
p G K
(19b)
2
4 4 2 2
f p
3 2
4 4 2 2 2
f p p
M c M
M 12 c c M
p K
A
p G K
(19c) 210
4 4
p
1 2
4 4 2 2 2
f p p
c M
M 12 c c M
p K
B
p G K
(19d)
2 2
p
2 2
4 4 2 2 2
f p p
6 c M
M 12 c c M
q G
B
p G K
(19e)
其中，
4
αf
4
Mβ
λ κ M
R
。
对超固结土而言，该模型在没有时间效应时退化为原始 UH 模型；当土既不受时间影响
又没有超固结度时，该模型退化为修正剑桥模型。若令 t 1d dp B t 为时间效应引起的等效平215
均主应力增量， t 2d dq B t 为时间效应引起的等效广义剪应力增量，则可将时间效应等效为
应力来考虑时间对土的影响。
三维弹粘塑性模型的演化规律
剪切蠕变
本文模型虽然基于一维蠕变规律建立，但在当前屈服面中与剪力耦合，因此可以描述土220
的三维蠕变。图 7、8 所示为土的不排水剪切持荷应力路径计算结果。其中图 7 为在低应力
比（q/p<M）下的剪切持荷，图 8 显示的为在高应力比下的剪切持荷（q/p>M）。
图 7 表明土在剪切持荷的作用下，最终会缓慢趋于破坏。其机理可以解释如下：时间效
应产生的粘塑性变形不断增大，而不排水条件下总体积不变，因此需要弹性变形减小来抵消
粘塑性体积变形，产生应力弹性卸载。卸载过程中 p 不断变小，抗剪强度随之减小，直到发225
生破坏。图 8 与图 7 的持荷过程类似，但发生在土的剪胀过程中，因此出现的路径是应力弹
性加载。
恒应变率剪切
对初始超固结度 OCR=8 的土以不同的应变率进行剪切，图 9、图 10 分别显示了应力比
和应力路径的变化。预测结果表明模型可以描述应变率效应，即应变率越高，不排水抗剪强230
- 9 -
中国科技论文在线
度越高，峰值应力比越大。并且，图中“2×10-6/min”、“2×10-5/min”两种速率大于特征
速率，以该两种速率剪切，土表现出超固结性，峰值应力比高于临界状态应力比，应力路径
越过临界状态线。而以小于特征速率的“2×10-4/min”速率剪切时，土则表现出欠固结性。
图 7 低应力比剪切持荷应力路径 235
Fig. 7 Stress paths at lower stress ratios
图 8 高应力比剪切持荷应力路径
Fig. 8 Stress paths at higher stress ratios 240
-04
-05
-06
剪应变
2×10-6
2×10-5
2×10-4
OCR=8
剪应变率：1/min η
图 9 初始超固结土加载应力路径
Fig. 9 Stress-strain relationships of overconsolidated clays
q
:k
P
a
0
20
40
60
80
100
120
140
160
180
0 50 100 150 200
OCR=2
OCR=4
OCR=8
临界状态线
p:kPa
0
20
40
60
80
100
120
140
0 20 40 60 80 100 120
OCR=1
OCR=2
OCR=4
p:kPa
q
:k
P
a
临界状态线
- 10 -
中国科技论文在线
0
50
100
150
200
250
0 50 100 150 200
-04
-05
-06
q:kPa
2×10-6
2×10-5
2×10-4
剪应变率：1/min
OCR=8
临界状态线
q
:k
P
a
245
图 10 初始超固结土加载应力路径
Fig. 10 Stress paths of overconsolidated clays
3 模型验证
本文将模型预测结果与香港海滩土的试验数据进行对比。香港海滩土的材料参数[6]如表
1 所示： 250
表 1 香港海滩土的土性参数
Tab. 1 Parameters of Hong Kong Marine Deposits
eλ0 pλ0/kPa λ κ M v β
50
与修正剑桥模型相同 粘性参数
注：表 1 中的 eλ0 、 pλ0两项参数用来表示瞬时正常压缩线的位置。本文模型总共需要 7 个参数，与修
正剑桥模型相比，只增加一个粘性参数 β。
255
图 11、12 是对初始状态为正常固结的土进行的常规三轴压缩试验预测模拟。模拟采用
三种剪切应变率加载，分别为 %/h、%/h、15%/h。土体所受的初始围压为 400kPa。
预测结果对试验数据能够合理反映。
0
50
100
150
200
250
300
350
瞬时
15%/h
%/h
%/h
15%/h
%/h
%/h
预测
剪应变
q
:k
P
a
试验
图 11 应力-应变关系对比 260
Fig. 11 Comparison of stress-strain relationships
- 11 -
中国科技论文在线
0
100
200
300
400
500
0 100 200 300 400 500 600
瞬时
15%/h
%/h
%/h
15%/h
%/h
%/h
q
:k
P
a
剪应变
预测
试验
图 12 应力路径对比
Fig. 12 Comparison of stress paths
4 结论 265
在一维蠕变规律的基础上，通过采用 UH 模型的再加载线，建立了一维弹粘塑性应力应
变关系。之后将折算时间引入 UH 模型的当前屈服面函数，建立了考虑剪切的三维弹粘塑性
本构模型。该模型能够描述土的剪胀、软化等超固结性以及蠕变、应变率效应等粘性规律，
同时与修正剑桥模型相比，仅增加一个参数。
[参考文献] (References)270
[1] KUTTER B L, SATHIALINGAM N. Elastic viscoplastic modelling of the rate-dependent behaviour of
clays[J]. Geotechnique, 1992, 42(3): 427-441.
[2] YIN J H, ZHU J G, GRAHAM J. A new elastic viscoplastic model for time-dependent behaviour of normally
and overconsolidated clays: theory and verification[J]. Can Geotech J, 2002, 39(1): 157-173.
[3] YIN Z Y, CHANG C S, KARSTUNEN M, et al. An anisotropic elastic-viscoplastic model for soft clays[J]. Int 275
J Solids Struct, 2010, 47(5): 665-677.
[4] YAO Y P, HOU W, ZHOU A N. UH model: three-dimensional unified hardening model for overconsolidated
clays[J]. Geotechnique, 2009, 59(5): 451-469.
[5] BJERRUM L. Engineering geology of Norwegian normally-consolidated marine clays as related to settlements
of buildings[J]. Geotechnique, 1967, 17(2): 83-118. 280
[6] ZHU J G. Experimental study and elastic visco-plastic modelling of the time-dependent stress-strain behaviour
of Hong Kong marine deposits[D]. Hong Kong: The Hong Kong Polytechnic niversity, 2000.
[7] YAO Y P, GAO Z W, ZHAO J D, et al. Modified UH model: constitutive modeling of overconsolidated clays
based on a parabolic Hvorslev envelope[J]. J Geotech Geoenviron Eng, 2012, 138(7): 860-868.
285
