膨胀土与水作用后会产生一系列的物理化学反应，引起膨胀土的膨胀效应和力学性能的改变^[1]，在实际工程中，常见到如基坑边坡膨胀土吸水膨胀而失稳、建筑物地基不均匀胀缩变形造成开裂等。因此，从理论和工程实践上解决膨胀土的工程灾害研究具有重大的工程意义。

非饱和膨胀土本构模型的研究起步于20世纪90年代，孙德安等^[2]和李舰^[3]在文献中回顾了非饱和膨胀土本构模型研究的发展历程：1990年Alonso等^[4-5]提出非饱和土的经典弹塑性模型 (BBM)，给出了非饱和土本构关系的基本理论框架。随后Gens等^[6]和Alonso等^[7]基于BBM提出了适用于膨胀性非饱和土的双尺度本构模型，并且将其称为BExM。卢再华等^[8]、曹雪山等^[9]和Sanchez等^[10]对本构模型进行了分析和改进。同时，Hoffmann等^[11]、Alonso等^[12]和Gens等^[13]将双尺度本构模型应用于解决膨胀性非饱和土的变形和渗流耦合的问题。李舰等^[14]结合非饱和膨胀土的双尺度本构模型BBxM和湿陷性非饱和土的毛细弹塑性变形耦合模型，建立了一个可预测非饱和膨胀土的毛细滞回和力学行为耦合的双尺度本构模型，2013年建立了适用于双孔隙结构非饱和膨胀土本构模型的理论框架^[15]，2014年从宏观角度建立一个能描述非饱和膨胀土的基坑力学性质以及吸力循环作用下土的行为的弹塑性本构模型^[16]。李吴刚等^[17]综合考虑G-A模型与SFG模型的优缺点，通过引入NL屈服面对传统SFG模型进行改进，提出更为简单易用的膨胀土本构模型。另外，Sheng^[18]、Chen^[19]从宏观角度建立了可以描述膨胀性非饱和土行为的力学和水力耦合的本构模型。

以上膨胀性非饱和土的弹塑性模型的研究并用没有直接考虑含水率或饱和度对膨胀土力学特性的影响，而是用吸力来表示非饱和的状态，工程应用难度较大。受热弹性力学理论启发，缪协兴等^[20]提出膨胀岩体中的弹性湿度应力场理论；随后朱珍德等^[21]在其基础上进行塑形修正，运用参变量变分原理对膨胀土塑形状态的本构模型进行了研究，提出了基于湿度应力场理论的膨胀岩弹塑性本构模型，并给出了相应的有限元形式，但是该模型并未对膨胀系数进行针对性的研究，应用于工程也较为复杂。陈茜等^[22]认为非饱和土的变形与含水率及其变化有关，并在非饱和土计算模型中引入含水率，建立了相关经验公式。本文旨在提出一个以湿度应力场理论为基础，具有工程实用价值的膨胀土弹塑性本构模型，简化目前膨胀土弹塑性本构研究中的塑性准则^[21-23]，在摩尔库仑准则的基础上，结合室内试验得到的含水量变化与变形、强度和膨胀参数变化之间的关系，提出基于摩尔库仑准则的膨胀土弹塑性本构模型，并通过FLAC3D软件所提供的二次开发程序接口实现自定义本构计算；并以成都东郊某膨胀土基坑边坡为实例，通过室内试验、渗流计算得到含水量变化与变形、强度和膨胀参数变化之间的关系以及湿度场分布，采用该本构模型进行数值计算，计算结果与现场监测结果相吻合，验证了该本构模型的正确性。

1 基于摩尔库仑准则的膨胀土弹塑性本构模型的建立
跳转到：本文顶部跳转到：本文顶部　1 基于摩尔库仑准则的膨胀土弹塑性本构模型的建立　    1.1 变形参数　    1.2 强度参数　    1.3 膨胀参数　2 膨胀土本构模型的二次开发　    2.1 FLAC3D的二次开发及程序流程图　    2.2 二次开发的实现　3 算例验证　    3.1 边坡含水量分布　    3.2 数值计算结果及分析　4 结论　参考文献

缪协兴受温度应力场理论的启发，提出了一种弹性湿度应力场理论。膨胀岩土在无约束条件下吸水会产生自由膨胀，给定含水量的变化ω(x, t)，x为位置坐标，t为时间，在弹性范围内的总应变为

${\varepsilon '_{ij}} = \alpha {\delta _{ij}}\omega$

(1)

式中：α为湿度膨胀系数；δ_ij为Kornecker记号；ε′_ij为含水量变化ω时产生的膨胀应变，可通过试验获得其关系曲线。

在有约束条件下膨胀岩土吸水时，ε′_ij不能自由发生，会产生膨胀应力，膨胀应力也要产生附加应变。因此，总应变变化为

${\varepsilon _{ij}} = {{\varepsilon ''}_{ij}} + {{\varepsilon '}_{ij}}$

(2)

式中：ε″_ij为附加应变，与膨胀应力之间服从Hooke定律，则总应变可表示为

${\varepsilon _{ij}} = \frac{{1 + \upsilon }}{E}{\sigma _{ij}}-\frac{\upsilon }{E}{\delta _{ij}}\sigma + \alpha {\delta _{ij}}\omega$

(3)

同时可写成总应力形式

${\sigma _{ij}} = \frac{{\upsilon E}}{{(1 + \upsilon )(1-2\upsilon )}}{\delta _{ij}}{\varepsilon _{kk}} + \frac{E}{{1 + \upsilon }}{\varepsilon _{ij}}-\alpha \frac{E}{{1-2\upsilon }}{\delta _{ij}}\omega$

(4)

式中：E、υ分别为弹性模量、泊松比，均为含水量的函数，可通过试验获得相关关系曲线。

式 (4) 即为弹性状态下的总应力表达式，等式右边前两项即为用E-v型模型表达的弹性模型，最后一项$\alpha \frac{E}{{1-2\upsilon }}{\delta _{ij}}\omega$即为膨胀效应产生的膨胀应力附加项，其中，$\frac{E}{{1-2\upsilon }}$可认为是膨胀效应时膨胀模量，是弹性模量E、泊松比υ的函数。

结合在Mohr-Coulomb模型中，当应力超过剪切、拉伸屈服准则，则进行塑形修正，其中剪切、拉伸屈服准则与膨胀岩土的强度参数c、φ有关，均为含水量的函数，可通过试验获得相关关系曲线，式 (5) 中f^s表示剪切准则，式 (7) 中f^t表示拉伸准则，σ^t表示拉伸强度。

${{f}^{\text{s}}}={{\sigma }_{1}}-{{\sigma }_{3}}{{N}_{\varphi }}+2c\sqrt{{{N}_{\varphi }}}$

(5)

${N_\varphi } = \frac{{1 + \sin \left( \varphi \right)}}{{1-\sin \left( \varphi \right)}}$

(6)

${{f}_{\text{t}}}={{\sigma }_{3}}-{{\sigma }^{\text{t}}}$

(7)

$\sigma _{{\text{max}}}^{\text{t}} = \frac{c}{{\tan \varphi }}$

(8)

上述方程再加上相关的边界条件、协调方程以及几何方程等就构成了非饱和膨胀土弹塑性本构方程，从而实现基于摩尔库仑模型的膨胀土弹塑性本构模型，该本构模型中变形、强度和膨胀参数与含水量的变化关系均可通过室内试验得到。

1.1 变形参数
跳转到：本文顶部跳转到：本文顶部　1 基于摩尔库仑准则的膨胀土弹塑性本构模型的建立　    1.1 变形参数　    1.2 强度参数　    1.3 膨胀参数　2 膨胀土本构模型的二次开发　    2.1 FLAC3D的二次开发及程序流程图　    2.2 二次开发的实现　3 算例验证　    3.1 边坡含水量分布　    3.2 数值计算结果及分析　4 结论　参考文献

当土中微元所受的应力较小时，由于塑形变形较小，可将土体视为弹性材料。此时的应力应变关系可通过Hooke定理进行计算，即σ=[D_e]{ε}。谢定义等^[24]对不同弹性参数表示的刚度矩阵进行了介绍，并讨论了常用的E-υ型模型和K-G型模型。韦秉旭^[23]通过GDS三轴试验对宁明膨胀土弹性模量E和泊松比υ进行了相关试验研究。本文在其基础上，采用GDS三轴仪对成都膨胀土进行k₀固结试验，试验土样取自成都东郊膨胀土分布区域，干密度1.7 g/cm³，天然含水量20.6%，缩限含水量13.2%，自由膨胀率50%。研究弹性模型以及泊松比随含水量ω和体积应力P的变化关系，试验结果图 1、图 2所示。

图 1

图 1 不同含水量下泊松比与体积应力的关系 Fig. 1 The relationship between poisson and volume stress under different water content

图 2

图 2 不同含水量下弹性模量与体积应力的关系 Fig. 2 The relationship between elasticity modulus and volume stress under different water content

根据不同含水量下泊松比、弹性模量与体积应力的变化曲线，进行多元线性回归，回归方程为

$\upsilon = 0.0015 \times P + 0.5 + \frac{{0.15-0.5}}{{1 + {{\text{e}}^{\frac{{\omega-0.205}}{{0.435}}}}}}$

(9)

$E = ({\text{0}}{\text{.63}}\omega {\text{ + 0}}{\text{.05}}) \times P-{\text{587}}{\text{.21}}\omega {\text{ + 166}}{\text{.26}}$

(10)

1.2 强度参数
跳转到：本文顶部跳转到：本文顶部　1 基于摩尔库仑准则的膨胀土弹塑性本构模型的建立　    1.1 变形参数　    1.2 强度参数　    1.3 膨胀参数　2 膨胀土本构模型的二次开发　    2.1 FLAC3D的二次开发及程序流程图　    2.2 二次开发的实现　3 算例验证　    3.1 边坡含水量分布　    3.2 数值计算结果及分析　4 结论　参考文献

弹塑性本构模型将应变分为弹性应变与塑性应变两部分，计算过程中先进行屈服判定，当应力超过剪切、拉伸屈服准则，产生的塑性应变按塑性理论计算。在摩尔库仑模型中屈服函数分为剪切屈服、拉伸屈服。由摩尔库仑模型屈服函数 (5)~(8) 可知，屈服判定与土中应力以及强度参数c、φ有关，通过不同含水量条件下成都膨胀土直剪试验，研究成都膨胀土强度参数随含水量的变化关系，试验结果见图 3~4。

图 3

图 3 不同初始含水率下c值的变化规律 Fig. 3 The relationship between cohesive force and initial water content

图 4

图 4 不同初始含水率下φ值的变化规律 Fig. 4 The relationship between frictional angle and initial water content

由图 4可知，强度参数c、φ随含水量的增加而降低，回归方程为：

$c =-255.9\omega + 99.9$

(11)

$\varphi =-301.3\omega + 98.7$

(12)

1.3 膨胀参数
跳转到：本文顶部跳转到：本文顶部　1 基于摩尔库仑准则的膨胀土弹塑性本构模型的建立　    1.1 变形参数　    1.2 强度参数　    1.3 膨胀参数　2 膨胀土本构模型的二次开发　    2.1 FLAC3D的二次开发及程序流程图　    2.2 二次开发的实现　3 算例验证　    3.1 边坡含水量分布　    3.2 数值计算结果及分析　4 结论　参考文献

膨胀土与普通粘土的不同在于其遇水膨胀、失水收缩的膨胀特性，而膨胀参数是膨胀土膨胀特性的力学指标。为了研究膨胀土膨胀特性的影响因素，解决工程应用问题，不少学者进行了大量的膨胀试验研究，研究表明，膨胀土的膨胀特性与土样的干密度、含水量有关^[25]。而在工程应用中，实际膨胀土工程干密度一定，导致工程出现变形破坏往往是因为降雨等涉水因素，降雨等导致土体含水量变化，但土体并未达到饱和状态，常规膨胀特性土工试验并不能满足工程应用。本文按照丁振洲^[26]提出的等同样试验方法对成都膨胀土进行膨胀率随含水量变化的试验研究，测得试验曲线如图 5所示。

图 5

图 5 膨胀率随过程含水量变化曲线 Fig. 5 The relationship between the expansive ratio and the process of water content

由图 5知，不同初始含水量条件下，土样自然膨胀力的增长趋势相近，对曲线形态进行近似拟合，见式 (13)，即为膨胀土弹塑性本构模型中膨胀应变的表达式。

${{\text{ }\!\!\varepsilon\!\!\text{ }}_{p}}=0.004{{\omega }_{0}}^{-1.863\ 6}\ln \frac{\text{e}(\omega -{{\omega }_{0}})+\left( 0.32-\omega \right)}{0.32-{{\omega }_{0}}}$

(13)

2 膨胀土本构模型的二次开发
跳转到：本文顶部跳转到：本文顶部　1 基于摩尔库仑准则的膨胀土弹塑性本构模型的建立　    1.1 变形参数　    1.2 强度参数　    1.3 膨胀参数　2 膨胀土本构模型的二次开发　    2.1 FLAC3D的二次开发及程序流程图　    2.2 二次开发的实现　3 算例验证　    3.1 边坡含水量分布　    3.2 数值计算结果及分析　4 结论　参考文献

2.1 FLAC3D的二次开发及程序流程图
跳转到：本文顶部跳转到：本文顶部　1 基于摩尔库仑准则的膨胀土弹塑性本构模型的建立　    1.1 变形参数　    1.2 强度参数　    1.3 膨胀参数　2 膨胀土本构模型的二次开发　    2.1 FLAC3D的二次开发及程序流程图　    2.2 二次开发的实现　3 算例验证　    3.1 边坡含水量分布　    3.2 数值计算结果及分析　4 结论　参考文献

目前，FLAC3D的自定义本构模型可采用Visual Studio 2005编程来创建。用户通过Visual Studio生成命令创建一个动态链接库文件 (后缀名为.dll)，这个动态链接库文件即为用来作为自定义本构模型的文件。在计算过程中FLAC3D程序会自动调用用户指定的动态链接库DLL文件，实现自定义本构模型的计算。

根据FLAC3D中摩尔库仑本构模型的编写过程，考虑基于摩尔库仑准则的膨胀土弹塑性本构模型程序流程图见图 6。

图 6

图 6 膨胀土本构模型程序流程图 Fig. 6 The program flow chart of the expansive soil constitutive model

2.2 二次开发的实现
跳转到：本文顶部跳转到：本文顶部　1 基于摩尔库仑准则的膨胀土弹塑性本构模型的建立　    1.1 变形参数　    1.2 强度参数　    1.3 膨胀参数　2 膨胀土本构模型的二次开发　    2.1 FLAC3D的二次开发及程序流程图　    2.2 二次开发的实现　3 算例验证　    3.1 边坡含水量分布　    3.2 数值计算结果及分析　4 结论　参考文献

膨胀土弹塑性本构模型的二次开发在Visual Studio 2005的环境中进行，主要开发工作包括头文件 (后缀为.h) 和C++源文件 (后缀为.cpp) 的修改。有3个头文件可不用修改，分别是stensor.h、axes.h和conmodel.h^[27-28]。其中stensor.h文件为张量头文件，用户根据文件中的定义可以得到当前单元应力应变关系得到当前应力张量及其应力增量张量以及应变增量张量等；Ases.h文件是坐标系头文件，主要用来定义坐标系统。Conmodel.h文件是本构模型结构体头文件，包含一个纯虚本构模型类以及两个结构体，主要用来描述子单元状态的变量。

头文件的修改主要包括模型编号以及私有变量的重新定义，包括模型的参数及迭代所需要的中间变量。源文件的修改是二次开发的关键所在。源文件中有几个关键函数，分别为Properties ()、SetProperty ()、Copy ()、Initialize ()、Run () 以及SaveRestore ()。最关键的两个函数是Initialize () 函数和Run () 函数。第1个函数是对模型计算中的变量进行初始化。在FLAC3D执行运行或执行大应变校正时，该函数执行一次。对于本文的膨胀土本构模型，需要初始化的参数有E、υ、K、G、c、φ。第2个函数是整个模型开发中最重要的函数，主要包括塑性状态判断、根据弹性状态下湿度应力应变关系 (式4、9、10) 计算三向主应力以及偏应力以及塑形判断与修正 (式5-8、11-12)。在FLAC3D在求解时会在每一个计算时间步内对每一个单元的子单元调用此函数。主程序通过重载第2个函数即为本构模型的实现过程。

采用Visual Studio 2005编程软件实现上述文件修改后，即可创建一个动态链接库文件。在数值计算过程中，通过主程序调用此动态链接库文件，即可实现自定义本构模型的计算。

3 算例验证
跳转到：本文顶部跳转到：本文顶部　1 基于摩尔库仑准则的膨胀土弹塑性本构模型的建立　    1.1 变形参数　    1.2 强度参数　    1.3 膨胀参数　2 膨胀土本构模型的二次开发　    2.1 FLAC3D的二次开发及程序流程图　    2.2 二次开发的实现　3 算例验证　    3.1 边坡含水量分布　    3.2 数值计算结果及分析　4 结论　参考文献

所选算例为成都东郊某膨胀土基坑边坡，所在区域为著名的成都粘土 (膨胀土) 的分布区域。基坑边坡土体主要以弱中等膨胀性的粘土层为主，粘土层天然含水率18%。所选边坡支护工程为单排桩，基坑长50 m，开挖深度6 m，锚固深度5 m；悬臂桩桩长11 m，桩径1 m；现场桩身变形测试点共3个，分别为1#、2#和3#。通过现场量测桩间距，建立FLAC3D基坑边坡数值计算模型如图 7所示。

图 7

图 7 基坑数值计算模型 Fig. 7 Numerical calculation model of foundation pit

3.1 边坡含水量分布
跳转到：本文顶部跳转到：本文顶部　1 基于摩尔库仑准则的膨胀土弹塑性本构模型的建立　    1.1 变形参数　    1.2 强度参数　    1.3 膨胀参数　2 膨胀土本构模型的二次开发　    2.1 FLAC3D的二次开发及程序流程图　    2.2 二次开发的实现　3 算例验证　    3.1 边坡含水量分布　    3.2 数值计算结果及分析　4 结论　参考文献

降雨和水位变动是导致膨胀土边坡失稳的主要外因。在本文算例中，由于降水井持续抽水，地下水位对基坑边坡的影响较小，因此，降雨是造成边坡变形的主要外因。现场监测结果表明，24 h持续大雨后，基坑边坡变形显著增加，以3#变形测试点为例，测试结果如图 8所示。

图 8

图 8 降雨前后边坡变形监测曲线 Fig. 8 Slope deformation monitoring curves before and after rainfall

袁俊平等^[29]采用有限元数值模拟方法分析了边坡地形、裂隙位置、裂隙开展深度及裂隙渗透特性等对边坡降雨入渗的影响，结果表明，坡上位置的裂隙对边坡入渗影响较大。结合在成都膨胀土地区裂隙统计调查结果，确定影响降雨入渗的裂隙深度为1 m左右。

结合以上研究成果以及边界条件，采用GeoStudio 2007渗流模块建立基坑边坡降雨入渗简化模型，计算结果如图 9所示，其中降雨强度10 mm/h，降雨24 h；初始体积含水量33%，饱和体积含水量53%，饱和渗透系数通过南55型渗透仪进行测定，测试结果饱和渗透系数为5.5×10^-8m/s，考虑裂隙分布饱和渗透系数为2×10^-7m/s，非饱和渗透系数通过土水特征曲线进行取值。

图 9

图 9 基坑边坡渗流计算结果云图 Fig. 9 The seepage calculation results nephogram of foundation pit

3.2 数值计算结果及分析
跳转到：本文顶部跳转到：本文顶部　1 基于摩尔库仑准则的膨胀土弹塑性本构模型的建立　    1.1 变形参数　    1.2 强度参数　    1.3 膨胀参数　2 膨胀土本构模型的二次开发　    2.1 FLAC3D的二次开发及程序流程图　    2.2 二次开发的实现　3 算例验证　    3.1 边坡含水量分布　    3.2 数值计算结果及分析　4 结论　参考文献

在进行降雨影响下数值计算之前，进行天然工况下的模拟计算，验证计算模型的正确性。提取数值模型变形测试点计算结果，与现场测试结果对比如图 10所示，由图可知，数值计算结果与现场变形相近，变形相差不超过±1 mm。

图 10

图 10 降雨前变形测试点变形对比曲线 Fig. 10 The contrast of the deformation monitoring curve before rainfall

根据3.1节中基坑边坡渗流计算结果，提取含水量的分布，赋值至基坑边坡数值计算模型，即为基坑边坡在大雨24 h后的湿度场，结合本文基于摩尔库伦模型的膨胀土本构模型便可进行降雨影响下边坡的数值分析，计算结果如图 11所示，由变形对比曲线可知，计算模型和现场边坡在降雨24 h后均产生了较大的变形，两者变形曲线相近，变形相差不超过±3 mm，验证了基于摩尔库仑模型的膨胀土本构模型的正确性。

图 11

图 11 降雨后变形监测点变形对比曲线 Fig. 11 The contrast of the deformation monitoring curve after rainfall

4 结论
跳转到：本文顶部跳转到：本文顶部　1 基于摩尔库仑准则的膨胀土弹塑性本构模型的建立　    1.1 变形参数　    1.2 强度参数　    1.3 膨胀参数　2 膨胀土本构模型的二次开发　    2.1 FLAC3D的二次开发及程序流程图　    2.2 二次开发的实现　3 算例验证　    3.1 边坡含水量分布　    3.2 数值计算结果及分析　4 结论　参考文献

1) 在湿度应力场理论基础上，考虑含水量变化与变形、强度和膨胀参数变化之间的关系，提出了基于摩尔库仑准则的膨胀土弹塑性本构模型。

2) 依据FLAC3D所提供的二次开发程序，结合基于摩尔库仑准则的膨胀土弹塑性本构模型，研究了数值软件二次开发程序运行的基本原理，给出了自定义本构模型的程序框图和代码编写中的几个关键技术。

3) 通过算例验证了二次开发的基于摩尔库仑模型的膨胀土本构模型程序的正确性与合理性。

参考文献跳转到：本文顶部跳转到：本文顶部　1 基于摩尔库仑准则的膨胀土弹塑性本构模型的建立　    1.1 变形参数　    1.2 强度参数　    1.3 膨胀参数　2 膨胀土本构模型的二次开发　    2.1 FLAC3D的二次开发及程序流程图　    2.2 二次开发的实现　3 算例验证　    3.1 边坡含水量分布　    3.2 数值计算结果及分析　4 结论　参考文献跳转到：本文顶部跳转到：本文顶部　1 基于摩尔库仑准则的膨胀土弹塑性本构模型的建立　    1.1 变形参数　    1.2 强度参数　    1.3 膨胀参数　2 膨胀土本构模型的二次开发　    2.1 FLAC3D的二次开发及程序流程图　    2.2 二次开发的实现　3 算例验证　    3.1 边坡含水量分布　    3.2 数值计算结果及分析　4 结论　参考文献

[1]	FREDLUNDD G, RAHARDJOH. 非饱和土土力学[M]. 北京: 中国建筑工业出版杜, 1997. FREDLUND D G, RAHARDJO H. Unsaturated soil mechanics[M]. Beijing: China Architecture & Building Press, 1997. (in Chinese)
[2]	孙德安. 非饱和土的水力和力学特性及其弹塑性描述[J]. 岩土力学, 2009, 30(11): 3217–3231. SUN D A. Hydro-mechanical behaviours of unsaturated soils and their elastoplastic modelling[J]. Rock and Soil Mechanics, 2009, 30(11): 3217–3231. (in Chinese)
[3]	李舰. 膨胀性非饱和土的本构模型的研究[D]. 北京: 北京交通大学, 2014. LI J. Research on the constitutive model for unsaturated expansive soils[D]. Beijing:Beijing Jiaotong University, 2014. (in Chinese)
[4]	ALONSO F E, GENS A, JOSA A. A constitutive model for partially saturated soils[J]. Géotechnique, 1990, 40(3): 405–430. DOI:10.1680/geot.1990.40.3.405
[5]	ALONSO E E, GENS A, JOSA A. A constitutive model for partially saturated soils[J]. Géotechnique, 1990, 40(3): 405–430. DOI:10.1680/geot.1990.40.3.405
[6]	GENS A, ALONSO E E. A framework for the behaviour of unsaturated expansive clays[J]. Canadian Geotechnical Journal, 1992, 29(6): 1013–1032. DOI:10.1139/t92-120
[7]	ALONSO E E, VAUNAT J, GENS A. Modelling the mechanical behaviour of expansive clays[J]. Engineering Geology, 1999, 54(1/2): 173–183.
[8]	卢再华, 王权民, 陈正汉. 非饱和膨胀土本构模型的试验研究及分析[J]. 地下空间, 2001, 21(5): 379–385. LU Z H, WANG Q M, CHEN Z H. Research on the constitutive model of unsaturated expansive soil[J]. Underground Space, 2001, 21(5): 379–385. (in Chinese)
[9]	曹雪山. 非饱和膨胀土的弹塑性本构模型研究[J]. 岩土工程学报, 2007, 27(7): 832–837. CAO X S. Elastoplastic constitutive model of unsaturated expansive soils[J]. Chinese Journal of Geotechnical Engineering, 2007, 27(7): 832–837. (in Chinese)
[10]	SÁNCHEZ M, GENS A, NASCIMENTO GUIMARAES L, et al. A double structure generalized plasticity model for expansive materials[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2005, 29(8): 751–787. DOI:10.1002/(ISSN)1096-9853
[11]	HOFFMANN C. Caracterización hidromecánica de mezclas de pellets de bentonita. Estudio experimental y constitutive [D]. Barcelona: Universitat Politècnica de Catalunya, 2005.
[12]	ALONSO E E, Romero E, HOFFMANN C. Hydromechanical behaviour of compacted granular expansive mixtures: experimental and constitutive study[J]. Géotechnique, 2011, 61(4): 329–344. DOI:10.1680/geot.2011.61.4.329
[13]	GENS SOLÉ A, VÁLLEJAN B, SÁNCHEZ M, et al. Hydromechanical behaviour of a heterogeneous compacted soil: experimental observations and modelling[J]. Geotechnique, 2011, 61(5): 367–386. DOI:10.1680/geot.SIP11.P.015
[14]	李舰, 赵成刚, 黄启迪. 膨胀性非饱和土的双尺度毛细-弹塑性变形耦合模型[J]. 岩土工程学报, 2012, 34(11): 2127–2133. LI J, ZHAO C G, HUANG Q D. Constitutive modeling with double-scale pore structure for coupling of capillary hysteresis and stress-strain behaviors in unsaturated expansive soils[J]. Chinese Journal of Geotechnical Engineering, 2012, 34(11): 2127–2133. (in Chinese)
[15]	李舰, 赵成刚, 蔡国庆, 等. 双孔隙结构膨胀性非饱和土功的表达式及其本构模型的理论框架[J]. 科学通报, 2013, 58(13): 1255–1263. LI J, ZHAO C Q, CAI G Q, et al. The input work expression and the thermodynamics-based modeling framework for unsaturated expansive soils with double porosity[J]. China Science Bull, 2013, 58(13): 1255–1263. (in Chinese)
[16]	李舰, 赵成刚, ASREAZADS. 适用于吸力循环作用的膨胀性非饱和土本构模型[J]. 岩土工程学报, 36(1): 132–139. LI J, ZHAO C G, ASREAZAD S. Constitutive model for unsaturated expansive soils under suction cycles[J]. Chinese Journal of Geotechnical Engineering, 36(1): 132–139. (in Chinese)
[17]	李吴刚, 杨庆, 刘文化, 等. 基于SFG模型的非饱和膨胀土本构模型研究[J]. 岩土工程学报, 2015, 37(8): 1449–1453. LI W G, YANG Q, LIU W H, et al. Constitutive model for unsaturated expansive clays based on SFG model[J]. Chinese Journal of Geotechnical Engineering, 2015, 37(8): 1449–1453. (in Chinese)
[18]	SHENG D, SLOAN S W, GENS A. A constitutive model for unsaturated soils: thermomechanical and computational aspects[J]. Computational Mechanics, 2004, 33(6): 453–465. DOI:10.1007/s00466-003-0545-x
[19]	CHEN R. Experimental study and constitutive modelling of stress-dependent coupled hydraulic hysteresis and mechanical behaviour of an unsaturated soil[D]. HongKong:HongKong University of Science and Technology, 2007.
[20]	缪协兴, 杨成永, 陈至达. 膨胀岩体中的湿度应力场理论[J]. 岩土力学, 1993, 14(4): 49–55. MIAO X X, YANG C D, CHEN Z D. The humidity in the expansion of rock mass stress field theory[J]. Rock and Soil Mechanics, 1993, 14(4): 49–55. (in Chinese)
[21]	朱珍德, 张爱军, 张勇, 等. 基于湿度应力场理论的膨胀岩弹塑性本构关系[J]. 岩土力学, 2004, 25(5): 700–702. ZHU Z D, ZHANG A J, ZHANG Y, et al. Elastoplastic constitutive law of swelling rock based on humidity stress fieldtheory[J]. Rock and Soil Mechanics, 2004, 25(5): 700–702. (in Chinese)
[22]	陈茜, 郭鸿. 非饱和土本构模型中考虑含水率的应力-应变曲线的构建方法[J]. 土木建筑与环境工程, 2015, 37(4): 97–104. CHEN X, GUO H. Method of stress-strain relationship curve considering water contents in unsaturated soil constitutive model[J]. Journal of Civil, Architectural & Environmental Engineering, 2015, 37(4): 97–104. (in Chinese)
[23]	韦秉旭. 宁明灰黑色膨胀土实用型本构模型研究[J]. 公路, 2008, 3: 151–156. WEI B X. The research on practical constitutive model for ningming expansive soils[J]. Highway, 2008, 3: 151–156. (in Chinese)
[24]	谢定义, 姚仰平, 党发宁. 高等土力学[M]. 北京: 高等教育出版社, 2008. XIE D Y, YAO Y P, DANG F N. Advanced Soil Mechanics[M]. Beijing: Higher Education Press, 2008. (in Chinese)
[25]	谭罗荣, 孔令伟. 膨胀土膨胀特性的变化规律研究[J]. 岩土力学, 2004, 25(10): 1555–1559. TAN L R, KONG L W. Study on variation regularity of swelling behavior of expansive soil[J]. Rock and Soil Mechanics, 2004, 25(10): 1555–1559. (in Chinese)
[26]	丁振洲, 郑颖人, 李利晟. 膨胀力变化规律试验研究[J]. 岩土力学, 2007, 28(7): 1328–1332. DING Z Z, ZHENG Y R, LI L S. Trial study on variation regularity of swelling force[J]. Rock and Soil Mechanics, 2007, 28(7): 1328–1332. (in Chinese)
[27]	褚卫江, 徐卫亚, 杨圣奇, 等. 基于FLAC3D岩体黏弹塑性流变模型的二次开发研究[J]. 岩土力学, 2006, 27(11): 2005–2010. CHU W J, XX W Y, YANG S Q, et al. Secondary development of a viscoelasto-plastic rheological constitutive model of rock based on FLAC3D[J]. Rock and Soil Mechanics, 2006, 27(11): 2005–2010. (in Chinese)
[28]	Itasca Consulting Group, Inc. Fast agrangian analysis of continua in three dimensions (version 2.1), user's manual[R]. Minnesota: Itasca Consulting Group, Inc., 2003.
[29]	袁俊平, 殷宗泽. 考虑裂隙非饱和膨胀土边坡入渗模型与数值模拟[J]. 岩土力学, 2004, 25(10): 1581–1586. YUAN J P, YIN Z Z. Numerical model and simulation of expansive soils slope infiltration considered fissures[J]. Rock and Soil Mechanics, 2004, 25(10): 1581–1586. (in Chinese)

X