硬化土本构模型在FLAC3D中的开发及应用
第33卷 第1期
岩石力学与工程学报 V ol.33 No.1
2014年1月 Chinese Journal of Rock Mechanics and Engineering Jan. ，2014
硬化土本构模型在FLAC 中的开发及应用
王春波12，丁文其12，乔亚飞1
，
，
，2
3D
(1. 同济大学 地下建筑与工程系，上海 200092；2. 同济大学 岩土及地下工程教育部重点实验室，上海 200092)
摘要：基于塑性增量理论，推导硬化土本构模型的有限差分格式。在VC++编程环境下，利用FLAC 3D 提供的二次开发平台，编制硬化土本构模型的有限差分程序，实现硬化土本构模型在FLAC 3D 中的二次开发，并给出二次开发的步骤、方法和编程要点。通过不同应力路径室内试验结果与采用不同本构模型的FLAC 3D 数值计算结果进行对比分析，验证硬化土本构模型在FLAC 3D 中二次开发的正确性和必要性。对比分析表明：采用硬化土本构模型的FLAC 3D 数值计算结果与试验结果吻合较好，而且硬化土本构模型的模型参数可从常规三轴试验获得，模型参数简单，硬化土本构模型不仅能够反映土体的非线性特性，还能够反映深基坑工程复杂的应力路径，与其他土体本构模型相比，硬化土本构模型更适合用于深基坑工程的计算分析。硬化土本构模型在FLAC 3D 中二次开发的 实现，扩大了FLAC 3D 的适用范围，在一定程度上弥补FLAC 3D 在分析岩土工程尤其是深基坑工程方面上的不足。关键词：基坑工程；复杂应力路径；硬化土本构模型；FLAC 3D ；二次开发；应力路径试验 中图分类号：TU 47 文献标识码：A 文章编号：1000–6915(2014)01–0199–10
DEVELOPMENT AND APPLICATION OF HARDENING SOIL
CONSTITUTIVE MODEL IN FLAC3D
WANG Chunbo12，DING Wenqi12，QIAO Yafei12
，
，
，
(1. Department of Geotechnical Engineering，Tongji University，Shanghai 200092，China ；
2. Key Laboratory of Geotechnical and Underground Engineering of Ministry of Education，Tongji University，
Shanghai 200092，China )
Abstract ：The finite difference scheme of hardening soil constitutive model is derived based on the increment theory of plastic. The finite difference program of hardening soil constitutive model is established so as that the secondary development of hardening soil constitutive model is realized with the platform for secondary development offered by FLAC3D software in VC++ environment. And then，the steps，methods and program essentials of constitutive model′s secondary development based on FLAC3D software are given. The comparison analysis between test results under different stress paths and FLAC3D numerical simulation results with different constitutive model is made so as to verify that the development of hardening soil constitutive model in FLAC3D is correct and necessity. Comparative analysis shows that the FLAC3D numerical simulation results with hardening soil constitutive model accord with test results. The parameters of hardening Soil constitutive model is simple and can be easily determined by conventional triaxial test. Hardening soil constitutive model can not only reflect the nonlinear characteristics of soil，but also can reflect complex stress path of deep foundation pit engineering. So，
收稿日期：2013–01–30；修回日期：2013–05–15
基金项目：国家自然科学基金资助项目(50878149)；国家青年科学基金项目(51008156)；无锡市科技计划资助项目(092167)
作者简介：王春波(1982–) ，男，2006年毕业于河北工程大学土木工程专业，现为博士研究生，主要从事深基坑工程流固耦合理论及岩土工程数值计算方面的研究工作。E-mail ：[email protected]
• 200 • 岩石力学与工程学报 2014年
hardening soil constitutive model is more suitable for computing analysis of deep foundation pit engineering compared with other constitutive model of soil. Implementation of secondary development that hardening soil constitutive model is developed in FLAC3D software expands the scope of FLAC3D application，to some extent，it makes up for disadvantages of FLAC3D in analysis of geotechnical engineering，especially deep foundation pit engineering.
Key words：foundation pit engineering；complex stress paths；hardening soil constitutive model；FLAC 3D ；secondary development；stress path tests
1 引 言
基坑工程尤其是深基坑工程应力路径复杂，土体本构模型的选取是进行数值计算分析的关键，目前用于基坑工程分析的土体本构模型主要有：邓肯–张(DC)模型、莫尔–库仑(MC)模型、修正剑桥(MCC)模型、硬化土(hardening soil)模型(简称HS 模型) 模型以及小应变硬化土(hardening soil model with small-strain)模型(简称HSS 模型) 等[1]。DC 模型为非线性弹性模型，可以反映土体应力、应变的非线性特性，但却不能反映土体的塑性应变，也不能反映基坑工程不同的应力路径；MC 模型会过高估计基坑坑底的隆起量；MCC 模型和HSS 模型模型参数较多且较难确定[2]。HS 模型不仅可以反映土体应力、应变的非线性特性和深基坑工程复杂的应力路径，而且模型参数可以从常规三轴试验获得，模型参数简单，因此成为基坑工程分析最广泛，也是最准确的土体本构模型之一[3-5]。目前，只在Plaxis 和Zsoil 软件中嵌有该模型。由ITASCA 咨询公司开发的三维有限差分程序FLAC 3D 采用差分有限元，不形成刚度矩阵，因此在大型岩土工程以及动力计算方面具有独特的优势。FLAC 3D 软件本身带有较多的本构模型，但在岩土工程领域广泛应用的土体本构模型(如DC 模型、HS 模型、HSS 模型等) 并没有包含在FLAC 3D 中[6-7]，这就使得软件在分析岩土工程时受到限制。基于此，众多科研人员对FLAC
3D
进行了二次开发工作，陈育民和刘汉龙
[8]
将DC 本构模型开发到FLAC 3D 软件中，并对其正确性进行了验证分析；蓝 航等[9]基于FLAC 3D 提供的二次开发平台，实现了节理岩体采动损伤本构模型开发；李英杰等[10]考虑变形模量在屈服后劣化的的应变软化本构模型，详细推导劣化损伤本构模型在FLAC 3D 软件中的差分格式，编制相应的开发流
程，实现了二次开发。
本文基于塑性增量理论，推导HS 模型的有限差分格式，并在VC++环境下编制HS 模型的有限差分程序，生成动态链接库文件，成功导入FLAC 3D 软件中。通过室内不同应力路径试验与采用不同本构模型的FLAC 3D 数值计算结果进行对比分析，以验证HS 模型在FLAC 3D 中二次开发的正确性和必要性。HS 模型在FLAC 3D 中二次开发的实现，扩大了FLAC 3D 软件的适用范围，在一定程度上弥补了FLAC 3D 在岩土工程尤其是深基坑工程计算分析方面的不足。
2 HS模型
R. L. Kondner[11]根据大量土的三轴试验的应力–应变关系曲线，提出可用双曲线拟合土在三轴试验加载下竖向应变ε1和偏应力q 之间的关系，后来这种关系被用在著名的双曲线模型即DC 模型[12]中。T. Schanz等[13]将双曲线形式的非线性弹性DC 模型扩展为弹塑性模型，并对其进行了修改完善。
当q ＜q f (其中，
q f 为土体强度，q f =(σ1−σ3) f ，由莫尔–库仑强度理论有q f =(2c cos ϕ+2σ3sin ϕ) / (1−sin ϕ) ) 时，土体处于弹性阶段，竖向应变ε1和偏应力q 之间满足双曲线关系：
εq a σ1−σ31=
2E =1q
(1a)
50q a −(σ1−σ3) 2E 501−q /q a
其中，
m
E 50=E
ref ⎛⎜σ3+c cot ϕ⎞50
⎝σref +c cot ϕ⎟⎠
(1b)
式中：E ref 50为加载模量；σ为相关应力，一般取
100 kPa；E ref
ref 50为相关应力σ时的加载模量；σ3为小主应力；c 为土体黏聚力；ϕ为土体内摩擦角；m
为幂指数；q a 为极限偏应力，q a =(σ1−σ3) ult 。设R f 为破坏比，且R f =(σ1−σ3) f /(σ1−σ3) ult =q f /q a (下标“utl ”表示极限，“f ”表示破坏) ，通过设定破
第33卷 第1期 王春波等：硬化土本构模型在FLAC 3D 中的开发及应用 • 201 •
坏比R f 取值(默认为0.9) 可求得极限偏应力q a 。
当q ≥q f 时，土体进入塑性阶段，产生塑性变形，随着硬化参数的变化，HS 模型屈服面也在不断的变化。
(1) 屈服函数
HS 模型为双屈服函数，包括剪切屈服和帽盖屈服函数，其屈服面[13]如图1所示。
o
(a) 主应力空间
c c
(b) p -q 平面
图1 HS模型屈服面
Fig.1 Yield surfaces of hardening soil constitutive model
① 剪切屈服函数表达式为
F q a q s =
E −2q
−γp (2)
50q a −q E ur
式中：E ⎛σ+c cot ϕ⎞m
ref 3ur 为回弹模量，E ur =E ur ⎜⎝σref +c cot ϕ⎟，⎠
E ref p ur 为相关应力σref 时的回弹模量；γ为塑性剪应变。
剪切屈服面以塑性剪应变γp 为硬化参数：
γp =εp p
p 1−ε2−ε3 (3)
式中：εp p
1，ε2
，εp 3分别为塑性第一、第二、第三
主应变。
② 帽盖屈服函数，也称体积屈服函数：
F =q
%2c M
2+p 2−p 2c (4)
式中：M 为摩擦常数，M =
6sin ϕ
3−sin ϕ
；p c 为前期固
结应力；p 为平均主应力，p =(σ1+σ2+σ3) /3；q %=σ1+(δ−1) σ2−δσ3，δ=(3+sin ϕ) /(3−sin ϕ) 。
帽盖屈服面以状态变量p c 为硬化参数，且有
m
d p ⎛⎜p c +c cot ϕ⎞p
c =H ⎝σref +c cot ϕ⎟⎠
d εv (5)
式中：H 为计算参数，H =
K s K c
K ，其中，K s 为
s −K c
回弹体积模量，且K E ref ur
s =3(1−2ν)
，ν为土体泊松比，
K 为割线体积模量，K =(1+2k 0) ref ref
c c 3
E oed ，E oed
为相关应力σref 时的割线模量，k 0为静止侧压力系数；εp v 为塑性体应变。
(2) 塑性势函数
① 剪切屈服面塑性势函数，采用非相关流动法
则，有
Q 1−σ3
σ1+σ3
s =
σ2
−
2
sin ψm (6) 式中：ψsin ϕm −sin ϕcv
m 为机动剪胀角，sin ψm =
1−sin ϕ，
m sin ϕcv
其中，ϕm 为机动摩擦角，sin ϕ1−σ3
m =
σσ，
1+σ3−2c cot ϕϕsin ϕsin ϕ−sin ψ
cv 为临界摩擦角，cv =1−sin ϕsin ψ
，ψ为土
体固有剪胀角。 ② 帽盖屈服面塑性势函数，采用相关联流动
法则，有
2Q q
%c =M
2+p 2 (7)
3 HS模型在FLAC 3D 中的二次开发
3.1 FLAC3D 计算步骤及二次开发环境
FLAC 3D 采用显示差分求解算法，即给定t 时
刻应力状态和Δt 时间步总应变增量，从而确定
t+Δt 时刻的应力状态[14]，首先根据Gauss 定律，由
• 202 • 岩石力学与工程学报 2014年
节点的速度求解单元应变率，利用应力–应变关系(即本构方程) ，计算单元应力，并对其进行积分，
可得到作用在节点上的应力矢量，然后再根据平衡方程，由节点力求解节点速率和位移，其计算流程如图
2所示。
图2 FLAC3D
计算流程图
Fig.2 Calculation procedure diagram of FLAC3D
FLAC 3D 采用面向对象的开发环境，其所有的
本构模型都以动态链接库文件(.DLL)形式提供给用户。自定义本构模型的主要功能是对给出的应变增量计算，得到新的应力[15]，模型文件的编写主要包括5部分内容：(1) 基类(class constitutive model)的描述；(2) 成员函数的描述；(3) 模型的注册；(4) 模型与FLAC 3D 之间的信息交换；(5) 模型状态指示器的描述。将编写好的头文件(H)和源文件(CPP)，通过编译、调试，生成.DLL 动态链接文件，将其复制到FLAC 3D
软件中即可通过命令加载和调用。 3.2 HS模型有限差分格式及开发流程
HS 模型在FLAC 3D
中的实现用到了主应力σ1，
σ2，σ3。主应力和主方向由应力张量确定，并按
顺序排列(压应力为负) ：
σ1＜σ2＜σ3 (8)
相应的主应变增量Δεi 可分解如下：
ΔεΔεe p
i =i +Δε i
(i = 1， 2 ，3) (9) 式中：上标“e ”和“p ”分别为弹性和塑性，且塑性分量仅在塑型流动阶段不为0，塑性应变由剪切屈服和体积屈服产生：
Δεp =Δεps pv i i +Δεi
(i = 1， 2 ，3) (10) 式中：上标“ps ”和“pv ”分别为塑性剪切和塑性体；Δε
pv
为塑性体应变，且Δε
pv
=(Δε
pv pv 1
+Δε
2
+
Δεpv 3) 为塑性体应变，以负号表示。
弹性部分应力增量可由广义Hooke 定律求得
Δσ(Δε e
i =S i i ) (i =
1， 2，3) (11) 即
Δσe +αe e
1=α1Δε12(Δε2+Δε3) ⎫
⎪Δσe
+αe ⎪⎪2=α1Δε22(Δε1+Δεe 3) ⎬ (12)
⎪
ΔσΔεe Δεe e ⎪3=α13+α2(1+Δε2) ⎪⎭
式中：α1，α2均为土体参数，且α1=K +4G /3，
α2=K −2G /3，其中，K 和G 分别为体积模量和
剪切模量，K =E t /3(1−2v ) ，G =E t /2(1+v ) ，加载时E t 取切线弹性模量E ，卸载时E t 取回弹模量
E ur 。
HS 模型给出的是加载模量E 50，增量形式的广
义Hooke 定律应用的为切线弹性模量E ，因此需要建立切线弹性模量E 和加载模量的E 50之间的关系式，在ε1-q 坐标平面内：
−2
E =d q 1ε1⎞d ε=⎛⎜1
+⎟ (13)
12E 50⎝2E 50q a ⎠
将式(1)代入式(13)可得
E =2E 50(1−q /q a ) 2 (14)
塑性剪切和体积屈服的流动法则分别为
Δεps ∂Q s ⎫i =λs ∂σ⎪i ⎪
⎬(i = 1， 2 ，3) (15) Δεpv
i =λv ∂Q c ⎪
∂σi ⎪⎭
式中：λs 为塑性剪应变系数，λv 为塑性体应变系数。
对式(6)求偏导可得各主方向的塑性剪应变：
Δεps ⎛11⎞⎫
1=λs ⎜⎝2−2sin ψm ⎟⎠⎪
⎪
Δεps
⎪2=0⎬ (16)
Δεps 11⎪
⎛⎞⎪
3=λs ⎜⎝−2−2sin ψm ⎟⎪
⎠⎭
塑性剪应变增量为
d γp =λs ⎛⎜∂Q s ∂σ−∂Q s −∂Q s ⎞
⎟=λs (17)
⎝1∂σ2∂σ3⎠
同理，对式(7)求偏导可得各主方向上的塑性体应变：
第33卷 第1期 王春波等：硬化土本构模型在FLAC 3D 中的开发及应用 • 203 •
Δεpv
=λv ⎛⎜2q
%2p ⎞⎫1⎝M
2+3⎟
⎠⎪⎪⎪
Δεpv
2
=λv ⎡%⎢2q
2p ⎤⎪⎪⎣M 2(δ−1) +3⎥⎦⎬ (18) ⎪
⎪
Δε
pv =λv ⎡⎢2q %2p ⎤⎪3
⎣M 2(−δ) +3⎥⎦⎪
⎪⎭
总的塑性体应变增量为
d εp
=λv ∂Q c
v
∂p
=λv (2p ) (19) 将式(19)代入式(5)可得硬化参数p c 的增量表达式：
m
d p ⎛⎜p c +c cot ϕ⎞
c =λv 2H ⎝σref +c cot ϕ⎟⎠
p (20)
如果发生剪切屈服，则剪切屈服函数F s (σI 1，
σI 3) ＜0，
将式(16)代入式(9)得
Δεe εs ⎛11⎞⎫
1=Δ1−λ⎜⎝2−2sin ψm ⎟⎠⎪
⎪
Δεe
Δε⎪2=2⎬ (21)
⎪
Δεe
εs ⎛11⎞⎪3=Δ3−λ⎜⎝−2−2sin ψm ⎟⎪
⎠⎭
基于弹性假设，由上标“I ”表示原来的应力σo
i 加上由总应变增量计算所得应力增量，则有
σI 1=σo 1+α1Δε1+α2(Δε2+Δε3) ⎪
⎫
σI o
+α⎪⎪2=σ21Δε2+α2(Δε1+Δε3) ⎬ (22)
⎪σI =σo ⎪
33+α1Δε3+α2(Δε1+Δε2) ⎪⎭
由硬化定律有F s (σo i +Δσi ) =0，若σo i 在屈服面上，即F o s (σi ) =0，则有
F s (Δσi ) = 0 (i = 1， 2 ，3) (23)
由总应变表示的应力增量为
Δσ⎛∂Q ⎞
i =S i (Δεi ) −λs S i ⎜s 1， 2 ，3) (24) ⎝∂σ⎟(i =
i ⎠
将式(24)代入式(23)可求得塑性剪应变系数：
λs = f s (σI i )
q −(i = 1， 2 ，3) (25)
a 122(1−2E q −
)
50a −(α1−α2) E ur
则剪切塑性修正后的应力为
σN I ⎡
⎛1
1⎞⎫
1=σ1−λs ⎢⎣
α1⎜⎝2
−2
sin ψm ⎟⎠
+⎪⎪
⎪ α⎛11⎞⎤⎪2⎜⎝−2−2sin ψm ⎟⎠⎥⎪⎦
⎪⎪σN I
2=σ2−λs αψ⎪2(−sin m ) ⎬ (26)
⎪σN I −λs ⎡⎛11⎞⎪
3=σ3⎢⎪⎣α1⎜⎝−2−2
sin ψm ⎟⎠+
⎪⎪ α⎛11⎞
⎤⎪
⎪2⎜⎝2−2sin ψm ⎟⎠⎥⎦⎪
⎭
式中：σN i (i = 1，2，3) 为当前时步主应力；
σI i =σo i +Δσi (i = 1，2，3) 为初始主应力。
如果发生体积屈服，即F I
I v (σI 1，σ2，σ3) ＜0，将
式(18)代入式(9)得
Δεe =Δε⎛2q %2p ⎞⎫
11−λv ⎜⎝M
2+3⎟
⎠⎪⎪⎪
Δεe
⎡⎢2q %2p ⎤⎪⎪2=Δε2−λv ⎣M 2(δ−1) +3⎥⎦⎬ (27) ⎪⎪
Δεe
=Δε⎡2q %2p ⎤⎪33−λv ⎢⎣M 2(−δ) +3⎥⎦⎪⎪⎭
塑性体积应变修正后的应力为
σN ⎛2q
⎫
1=σI 1−λv ⎣⎢⎡
α1⎜
%⎝M 2
+2p ⎞3⎟⎠
+⎪
⎪
α⎛2q %4p ⎞⎤⎪
2⎜⎝−M 2+3⎟⎠⎥⎦K ⎪
s ⎪⎪
σN =σI v ⎡⎛2q %2p ⎞⎪
22−λ⎢⎣
α1⎜⎪⎝M 2(δ−1) +3⎟⎠+⎪⎪⎬ (28) α⎡%⎪⎢4p ⎣3+2q ⎤
2M 2(1−δ) ⎥⎦K s ⎪
⎪
σN =σI ⎧⎡2q %2p ⎤
⎪
33−λv ⎨⎩α1⎢⎣M
2(−δ) +3⎦+⎪
⎪⎪ α⎛4p 2q %⎞⎫⎪
2⎜⎝3+M 2δ⎟⎪⎠⎬⎭K s
⎪⎭
同理可求塑性体应变参数λv 为
• 204 • 岩石力学与工程学报 2014年
λv =
f v (σI
n )
16q %2
(α 1−α2) 2(1−δ+δ2) 2M
+49p 2
(α1+2α2) 2
(29)
根据上述推导的HS 模型差分格式，图3给出了基于FLAC 3D 的HS 模型开发的流程图，可以用Visual Studio 2008编译器进行程序的编制，限于篇
幅，本文不给出程序代码。
图3 HS模型程序开发流程图
Fig.3 Program flow chart of hardening soil constitutive model
4 不同应力路径试验验证分析
为了验证HS 模型在FLAC 3D 中开发的正确性，本文以不同应力路径试验结果与数值计算结果进行
对比分析，并选取不同本构模型进行研究，以说明
HS 模型在深基坑工程分析中的优越性，深基坑工
程中不同区域土体应力路径如图4所示。
图4 深基坑工程不同区域土体应力路径示意图
Fig.4 Sketch of stress path of soil in different zones of
foundation pit engineering
4.1 试验仪器
试验仪器采用英国GDS 公司生产的STDTTS+ UNSAT(7 kN/1 700 kPa)型号高级应力路径三轴测
试系统，该系统主要由B&W
压力室、GDS 压力/体积数字控制器、霍尔效应传感器三部分组成。与
常规三轴仪相比，其主要优点如下：(1) 配备先进的压力与体积控制器；(2) 配备高精度的霍尔效应传感器；(3) 荷载传感器内置，可以实现复杂的应 力路径试验；(4) 采用独特的试样帽加压方式，可以进行土体拉伸试验，试验采用GDSLAB 软件进行试验控制和数据采集，因此系统在控制方式和数据采集方面全部均可自动控制。 4.2 试验土样及试验步骤
(1) 试验土样
土样取自无锡地铁1号线胜利门广场③1黏土层，土样埋深10 m。③1层黏土为第四系全新统上 段河湖相沉积物，可塑～硬塑，含铁锰结核，干强度高，韧性高，土样密度1.9×103 kg/m3，含水量
24.47%，天然孔隙比0.71。
(2) 试验步骤
① 测定K 0值。开启GDS 系统，选用Saturation
Ramp 试验模块，设置围压105 kPa，反压100 kPa，
当孔压系数B = Δu /Δσ = 0.98时，饱和完成。再选用K-ZERO 模块分别选取有效固结围压100，150 kPa ，测得相应的侧压力系数K 1 = 0.58，K 2 = 0.46，
其平均值K 0= 0.52。而现场扁铲原位测试所得K 0 = 0.48，因此本试验侧压力系数K 0 取0.5。
第33卷 第1期 王春波等：硬化土本构模型在FLAC 3D 中的开发及应用 • 205 •
② 土样固结。选用Consolidation 试验模块，loading 试验模块，保持围压不变，以0.2 kPa/min
保持反压100 kPa，设置有效围压100 kPa，进行等压固结。选用Advance loading试验模块，保持反压110 kPa，根据试验要求设定固结轴压，并由应力比
减小轴压至试样破坏，各围压下的应力路径试验及如表1所示。
④ 卸下土样：将所有压力回0，待排除压力室中的水后，拆去压力室外罩，取下橡皮套，记录试样破坏形状，测定试验后土样含水量。 4.3 数值分析模型
K 0=σ3/σ1设定固结围压，以0.3 kPa/min的速率加
载预设固结压力值，进行K 0固结。待压力设置稳定后打开反压阀门排水，直至试样固结度达到98%，固结完毕。
③ 剪切试验。土样固结完成后关闭排水阀门，
由试验确定的各模型计算参数如表2所示。计算模型与试验土样模型大小一致，直径38 mm、高76 mm，共划分单元1 900个，模型底部为竖向位
进行固结不排水下的3种应力路径试验：常规三轴压缩应力路径试验CUCTC ，选择标准三轴试验模块，保持围压不变，以0.05 mm/min的轴向位移速率增加轴向变形至试样破坏；减压三轴压缩应力路径试验CURTC ，选择Advance loading试验模块，保持轴压不变，以0.2 kPa/min减小围压至试样破坏；常规三轴轴向减压试验CURTE ，选择Advance
移约束，周边采用应力边界条件，数值计算模型如图5所示，力学模型分别采用本文开发的HS 模型、FLAC 3D 自带的修正剑桥模型(MCC模型) 和莫尔–库仑模型(MC模型) ，计算收敛不平衡力比例取105。
－
4.4 试验结果与数值计算结果对比分析
表1 试验应力路径 Table 1 Test stress paths
固结方式
固结压力/kPa
试验应力路径
模拟基坑工程位置 基坑超载处 基坑侧壁 基坑坑底 基坑超载处 基坑侧壁 基坑坑底 –
CUCTC80：围压不变，轴向加载
围压80 轴压160
CURTC80：围压减小，轴压不变 CURTE80：围压不变，轴向卸载
K 0固结
CUCTC150：围压不变，轴向加载
围压150 轴压300
CURTC150：围压减小，轴压不变 CURTE150：围压不变，轴向卸载
等压固结 围压100 轴向加、卸载
表2 各本构模型计算参数
Table 2 Calculation parameters of various constitutive models
本构模型名称 硬化土模型
ref
E 50/MPa
ref
E oed /MPa
ref
E ur /MPa
模型参数
m c /kPaϕ/(°) ψ/(°) σref /kPa OCR
R f
v
(HS模型)
修正剑桥模型
10.8 5.4 32.4 0.8 41.4 9.4 0.0 100 1.0 0.9 0.34
M
C c
λ
C s
κ
P c /kPa
(MCC模型)
莫尔–库仑模型
1.16 0.086 0.037 2 0.019 0.008 25 65
E /MPa
v
c /kPa
ϕ/(°) ψ/(°)
σt /kPa
(MC模型)
32.4 0.34 41.4 9.4 0.0 0.0
注：OCR 为超固结比，M 为摩擦常数，C c 为压缩指数，λ为正常固结坡度，C s 为回弹指数，κ为正常膨胀坡度，P c 为前期固结压力，σt 为抗拉强度。
• 206 • 岩石力学与工程学报 2014年
不同应力路径试验结果与FLAC 3D 选用不同本构模型的计算结果对比如图6所示。
图5 FLAC3D 计算模型 Fig.5 Calculation model of FLAC3D
图6(a)，(b)分别为K 0固结后轴向压缩试验值与采用FLAC 3D 不同本构模型计算值对比分析曲线，由图可知，HS 模型在轴向加载作用下，轴向应变与偏应力满足双曲线关系，而且随着轴向应变的增加，偏应力增大，表现出硬化特性，适用于硬土的计算分析，由于
HS 模型也采用莫尔–库仑强
a P k /力应偏CUCTC80试验值 HS 模型(FLAC3D ) MCC 模型(FLAC3D ) MC 模型(FLAC3D )
轴向应变/%
(a)
K 0固结后轴向压缩，σ3=80 kPa，σ1=160 kPa
a P k /力应偏CUCTC150试验值 HS 模型(FLAC3D ) MCC 模型(FLAC3D ) MC 模型(FLAC3D )
轴向应变/%
(b) K 0固结后轴向压缩，σ3=150 kPa，σ1=300 kPa
a P k /力应偏CURTC80试验值
HS 模型(FLAC3D ) MCC 模型(FLAC3D ) MC 模型(FLAC3D )
轴向应变/%
(c) K 0固结后侧向卸载，σ3=80 kPa，σ1=160 kPa
a P k /力应偏CURTC150试验值
HS 模型(FLAC3D ) MCC 模型(FLAC3D ) MC 模型(FLAC3D )
轴向应变/%
(d) K 0固结后侧向卸载，σ3=150 kPa，σ1=300 kPa
CURTE80试验值 HS 模型(FLAC3D ) MCC 模型(FLAC3D ) MC 模型(FLAC3D )
a P k /力应偏
轴向应变/%
(e) K 0固结后轴向卸载σ3=80 kPa，σ1=160 kPa
CURTE150试验值
HS 模型(FLAC3D ) MCC 模型(FLAC3D 3D ) MC 模型(FLAC)
a P k /力应偏
轴向应变/%
(f) K 0固结后轴向卸载，σ3=150 kPa，σ1=300 kPa
第33卷 第1期 王春波等：硬化土本构模型在FLAC 3D 中的开发及应用 • 207 •
a P k /力应偏试验值
HS 模型(FLAC3D ) MCC 模型(FLAC3D ) MC 模型(FLAC3D )
轴向应变/% (g) σ1=σ3=100 kPa
图6 试验值与计算值对比
Fig.6 Contrast between test value and calculation value
度准则，因此其最大偏应力与MC 模型的偏应力接近；MC 模型为线弹塑性模型，应力–应变曲线为线性，MCC 模型基本上能够反映土体的应力–应变非线性关系，但硬化特性不明显，最大偏应力较试验值偏小。
图6(c)，(d)分别为K 0固结后侧向卸载的试验结果与FLAC 3D 计算结果对比分析曲线，
在侧向卸载、轴压不变的应力路径下，偏应力与轴向应变之间仍然满足双曲线关系，HS 模型采用剪应变、固结压力为硬化参数，应力–应变关系完全表现为应变硬化特性，其硬化特性明显，故本文开发的HS 模型计算结果与试验结果吻合较好。MC 模型为理想线弹塑性模型，不能反映应力–应变的非线性特性，当土体进入塑性阶段以后，其应力大小保持恒定。MCC 模型采用塑性体应变和剪应变为硬化参数，当
塑性体应变增量为正(即体积压缩) 时，屈服面向外扩张，表现为应变硬化特性；当塑性体应变为负(即体积膨胀) 时，屈服面向内收缩，则表现为应变软化特性，故MCC 模型为应变硬化软化模型，即应力随应变的增加而增加，而后应力逐渐降低，其硬化特性不如HS 模型明显。
MCC 模型更适合软黏土的计算分析，由此模型计算得到的偏应力较试验值偏低。
图6(e)，(f)分别为K 0固结后轴向卸载试验结果与采用FLAC 3D 不同本构模型计算结果对比分析曲线。由轴向卸载的试验曲线可知，随着轴压的减少，会产生负的轴向应变(即土样被拉伸) ，卸载初期，
轴向应变变化速率较小，当轴向应力小于围压后，偏应力为负值，随着轴压的进一步减小，轴向应变突然增加，最后在偏应力不变时，轴向应变还会继
续增大，说明轴向卸载时轴向应变会有一个延迟过程。HS 模型的计算结果与试验值接近；MC 模型偏
应力随轴向应变变化缓慢；
MCC 模型为硬化软化模型，即轴向应力随轴向应变的增加有先增大后减小，所以当土样破坏时，由MCC 模型计算的应力较MC
模型、HS 模型计算的应力值要小，轴向卸载时
σ1＜σ3，偏应力σ1−σ3为负值，故偏应力的绝对值
较试验值偏大，这与轴向压缩、侧向卸载试验中由MCC 模型计算得到的偏应力(此时σ1＞σ3，σ1−σ3
为正值) 较试验值偏低的结论是一致的。
图6(g)为等压固结后轴向加、卸载试验结果与数值计算结果对比曲线，由图可知，采用本文开发的HS 模型的计算结果与试验结果吻合较好，加载时竖向应变与偏应力成双曲线关系，卸载、再加载时具有明显的“滞回圈”，可见，HS 模型比较适合分析具有复杂应力路径的岩土工程问题。MC 模型计算得到的竖向应变与偏应力基本成线性关系，偏应力达到一定值后便保持不变，这与MC 模型采用的理论基础是一致的。
MCC 模型弹性变形为非线性弹性变形，塑性阶段以塑性体积变形表示，初始阶段MCC 模型与试验值拟合较好，但应力–应变曲线硬化特性不明显，逐渐表现出硬化软化特性。
5 结 论
(1) 基于塑性增量理论，推导了HS 模型的有限
差分格式。
(2) 在VC++编程环境下，利用FLAC 3D 提供的
二次开发平台，编制了HS 模型的有限差分程序。
(3) 通过室内不同应力路径试验结果，验证了HS 模型在FLAC 3D 中二次开发的正确性和必要性，
不同本构模型的计算结果对比表明，HS 模型较其他本构模型更适合于深基坑工程的计算分析。
(4) HS模型在FLAC 3D 中二次开发的实现，扩
大了FLAC 3D 的适用范围，在一定程度上弥补了FLAC 3D 在分析岩土工程尤其是深基坑工程方面的
不足。
• 208 • 岩石力学与工程学报 2014年
参考文献(References)：
[1] 李广信. 高等土力学[M]. 北京：清华大学出版社，2004：32–
113.(LIU Guangxin. Advanced soil mechanics[M]. Beijing：Tsinghua University Press，2004：32–113.(in Chinese))
[2] 管 飞. 基于HSS 本构模型的软土超大型深基坑3D 数值分析[J].
岩土工程学报，2010，32(增1) ：177–180.(MGUAN Fei. 3D
numerical analysis for a ultra-large deep excavation in soft clay based on HSS constitutive model[J]. Chinese Journal of Geotechnical Engineering ，2010，32(Supp.1)：177–180.(in Chinese))
[3] 刘国彬，王卫东. 基坑工程手册[M]. 北京：中国建筑工业出版社，
2009：213–226.(LIU Guobin，WANG Weidong. The manual of foundation pit engineering[M]. Beijing：China Architecture and Building Press，2009：213–226.(in Chinese))
[4] 赵锡宏，李 蓓，杨国祥，等. 大型超深基坑工程实践与理论[M].
北京：人民交通出版社，2005：116–141.(ZHAO Xihong，LI Bei，
YANG Guoxiang，et al. Practice and theory for specially big and deep excavation engineering[M]. Beijing：China Communications Press，2005：116–141.(in Chinese))
[5] 马晓文，艾英钵. 基坑开挖土体卸荷特性研究进展[J]. 岩土工程学
报，2011，33(增1) ：182–186.(MA Xiaowen，AI Yingbo. Review of
soil behaviors under excavation unloading of foundation pits[J]. Chinese Journal of Geotechnical Engineering，2011，33(Supp.1)：182–186.(in Chinese))
[6] 陈育民，刘汉龙. FLAC/FLAC3D 基础与工程实例[M]. 北京：中国
水利水电出版社，2009：51–60.(CHEN Yumin，LIU Hanlong.
Fundamentals of FLAC/FLAC3D with engineering example[M]. Beijing ：China Water Power Press，2009：51–60.(in Chinese)) [7] 张书伟，林 杭，任连伟. FLAC3D 在岩土工程中的应用[M]. 北京：
中国水利水电出版社，2011：83–107.(ZHANG Shuwei，LIN Hang，
REN Lianwei. The Application of FLAC3D to geotechnical engineering[M]. Beijing：China Water Power Press，2011：83–107.(in Chinese))
[8] 陈育民，刘汉龙. 邓肯–张本构模型在FLAC 3D 中的开发与实现[J].
岩土力学，2007，28(10)：2 123–2 127.(CHEN Yumin，LIU Hanlong.
Development and implementation of Duncan-Chang constitutive model in FLAC3D [J]. Rock and Soil Mechanics，2007，28(10)：2 123– 2 127.(in Chinese)) [9]
蓝 航，姚建国，张华兴，等. 基于FLAC 3D 的节理岩体采动损伤本构模型的开发及应用[J]. 岩石力学与工程学报，2008，27(3)：
572–579.(LAN Hang，YAO Jianguo，ZHANG Huaxing，et at. Development and application of constitutive model of jointed rock mass damage due to mining based on FLAC3D [J]. Chinese Journal of Rock Mechanics and Engineering，2008，27(3)：572–579.(in Chinese))
[10] 李英杰，张顶立，刘保国. 考虑变形模量劣化的应变软化模型在
FLAC 3D 中的开发与验证[J]. 岩土力学，2011，32(增2) ：647–653.(LI Yingjie，ZHANG Dingli，LIU Baoguo. Development and verification of strain-softening model considering deformation modulus degradation in FLAC3D [J]. Rock and Soil Mechanics，2011，32(Supp.2)：647–653.(in Chinese))
[11] KONDNER R L. Hyperbolic stress-strain response：cohesive soils[J].
Journal of the Soil Mechanics and Foundations Division，1963，89(1)：115–143.
[12] DUNCAN J M，CHANG C Y. Nonlinear analysis of stress and strain
in soils[J]. Journal of the Soil Mechanics and Foundations Division，ASCE ，1970，96(5)：1629–1653.
[13] SCHANZ T，VERMEER P A，BONNIER P G. The hardening soil
model ：formulation and verification[C]// Beyond 2000 in Computation Geotechnics-10 Years of PLAXIS. Balkema：ISBN90 5809 040X，1999：281–296
[14] 刘 波，韩彦辉. FLAC原理、实例与工程应用[M]. 北京：人民交
通出版社，2005：7–14.(LIU Bo，HAN Yanhui. The theory，examples
and engineering application for FLAC[M]. Beijing：China Communications Press ，2005：7–14.(in Chinese))
[15] 张传庆，周 辉，冯夏庭. 统一弹塑性本构模型在FLAC 3D 中的计
算格式[J]. 岩土力学，2008，29(3)：596–602.(ZHANG Chuanqing，
ZHOU Hui，FENG Xiating. Numerical format of elastoplastic constitutive model based on the unified strength theory in FLAC3D [J].
Rock and Soil Mechanics，2008，29(3)：596–602.(in Chinese))