ALONSO等^[1]在修正剑桥模型的基础上提出了非饱和土的巴塞罗纳模型(也称为BBM模型)。BBM模型采用净应力和基质吸力作为模型中的应力状态变量，其特点在于提出的LC屈服线可以很好地描述非饱和土的屈服应力随基质吸力变化的性质。SHENG等^[2]基于双应力变量法提出了一个新的本构模型—SFG模型。该模型认为净应力和基质吸力对体应变的影响需要分开考虑，同时该模型采用积分的方法来计算LC屈服面方程。SFG模型可以很好地表现基质吸力对土体压缩性的影响，但该模型无法从非饱和状态退化到饱和状态，且无法考虑饱和度对非饱和土性质的影响。刘艳^[3]基于广义有效应力原理，采用连续介质力学中热力学的方法建立了固液气三相耦合的非饱和土本构模型，其特点在于可综合考虑固液气三相对非饱和土力学性质的影响，还考虑到了饱和度的影响。目前，大多数的非饱和土本构模型都是在临界状态土力学理论框架下建立起来的^[4-6]，而破坏应力比是一个非常重要的参数，由土的强度准则得出，强度准则大多采用基于某个有效应力原理表示的摩尔-库仑强度准则^[7] ，但该准则没有考虑中间主应力的影响。为了克服上述弊端，一些研究者相继提出了非饱和土的SMP准则^[8]和双剪统一强度准则^[9]等。其中，SMP准则也因其模型特征单一不能反映土体的多样性。而双剪统一强度准则虽然能较全面地表征非饱和土的强度特性，但存在双重破坏角的现象^[10]。而三剪统一强度准则^[10-11]则能够较全面地反映土的拉压不等效应、中间主应力效应和应力区间效应等特性，同时还不会出现双重破坏角的问题，而且该强度准则是一个强度准则系列，可以对其他准则作非线性逼近，由此以反映土体的多样性。在非饱和土的有效应力原理方面，采用较多的是Bishop应力变量法^[12]和Fredlund应力变量法^[13]。但是，Bishop应力变量法^[12]存在有效应力系数难以确定并且和饱和度对应关系不唯一的问题^[14-15]，Fredlund应力变量法^[13]最明显的缺点就是无法考虑饱和度的影响^[16]且无法描述土体从非饱和状态到饱和状态间的平稳过渡^[17-18]，原因是双应力变量有效应力原理在非饱和土基质吸力为零时不能与饱和土的有效应力原理相衔接。而赵成刚等^[19-21]提出的非饱和土广义有效应力原理可以克服上述缺点。

综上所述，本文的研究内容主要有：

1)将广义有效应力变量法与三剪强度准则^[22-23]相结合，得到非饱和土广义有效应力三剪强度准则，并进行验证。为了更好地在修正剑桥模型中反映非饱和黏性土黏聚力的影响，采用文献[24-25]提出的等量代换法和坐标平移法推导出相应的破坏应力比。

2)将所得破坏应力比与非饱和土修正剑桥模型相结合得到屈服函数。

3)在弹塑性理论的框架下，建立弹塑性本构模型。

4)用江西正常固结非饱和重塑红黏土对模型进行验证。验证内容为不同压实度、不同净围压、不同基质吸力条件下土在常规三轴固结排水试验中的偏应力和轴向应变、体应变和轴向应变的关系。

5)用所建立的本构模型对该红黏土进行真三轴固结排水数值模拟，并探究中间主应力的影响。

1   非饱和土三剪强度准则与破坏应力比

将非饱和土广义有效应力^[19]代入到三剪强度准则^[10-11]中，得到基于广义有效应力的非饱和土三剪强度准则为：

式中：$c'$为饱和有效黏聚力；$\varphi '$为饱和有效内摩擦角；$b$为中间主应力影响系数；${\widetilde \sigma '_1}$、${\widetilde \sigma '_2}$、${\widetilde \sigma '_3}$分别为最大、中间和最小广义有效主应力。

非饱和土广义有效应力具体表达式为^[3]：

式中：${\widetilde \sigma _{ij}}$为广义有效应力；${s_{\text{r}}}$为饱和度；$s$为基质吸力；${p_{\text{a}}}$为孔隙气压力；${\delta _{ij}}$为Kronecker符号。

其广义有效主应力为：

式中：${\sigma _1}$、${\sigma _2}$和${\sigma _3}$分别为最大、中间和最小总主应力；${\sigma '_1}$、${\sigma '_2}$和${\sigma '_3}$分别为最大、中间和最小净主应力。

将式(3)代入式(1)得到用净应力表示的基于广义有效应力的非饱和土三剪强度准则为：

净主应力及其相关的应力不变量关系式为^[26-27]：

式中：$p'$为平均净主应力；$\rho$为$\pi$平面上极限线上的点到静水压力轴的垂直距离；$q$为广义剪应力；$\theta$为应力角(${0 \leqslant \theta \leqslant {{60} }}$°)。

由式(4)和式(5)得到强度准则$p'$-$q$子午线表达式为：

其中：

由式(6)得到基于等量代换法的破坏应力比为：

修正剑桥模型中临界状态极限线是过坐标原点的线，黏聚力不等于零的黏性土并不符合其条件。采用坐标平移法使子午线经过原点，即将坐标轴向左平移使其经过原点。基于坐标平移法下新的非饱和土三剪破坏应力比为：

坐标平移前后其它应力状态量间的关系为：

式中：${\overline I'_1}$、${\overline J_2}$、${\overline J_3}$和$\overline \theta$分别为坐标平移后的非饱和土净应力的第一应力不变量、第二应力偏张量不变量、第三应力偏张量不变量和Lode角。

坐标平移前后的土体压缩指数关系为：

式中，$\lambda$和$\overline \lambda$分别为坐标平移前和坐标平移后的土体压缩指数。

坐标平移前后的土体回弹指数关系为：

式中，$\kappa$和$\overline \kappa$分别为坐标平移前和坐标平移后的土体回弹指数。其中：

文献[28]中的非饱和黄土在$s = 100\;{\text{kPa}}$时所对应的饱和度${s_{\text{r}}} = 34.5\text{%}$，有效黏聚力$c' = 5.3\;{\text{kPa}}$，有效内摩擦角$\varphi '{\text{ = }}27.92$°。图1为本文新建立的非饱和土广义有效应力三剪强度准则计算结果与文献[28]中的真三轴试验数据对比验证图。

由图1可知，取$b{\text{ = }}0.25$的非饱和土广义有效应力三剪强度准则计算预测结果与试验结果较为吻合，验证了本文所推导的新准则的正确性。新准则采用了广义有效应力原理，具备广义有效应力原理的优点，考虑了饱和度的影响，可以更全面地描述非饱和土的力学特性。

图  1  广义有效应力变量法计算结果和试验数据对比图

Figure  1.  Comparison diagram of calculation results of generalized effective stress variable method and test data

下载: 全尺寸图片 幻灯片

由图2可知，将广义有效应力变量法预测的$p'$-$q$子午线与试验点相比对，广义有效应力变量法预测的$p'$-$q$子午线基本与试验点吻合，新准则取得很好的预测效果。

图  2  子午线计算结果与试验数据对比图

Figure  2.  Comparison diagram of meridian calculation results and experimental data

下载: 全尺寸图片 幻灯片

文献[29-30]中的试验所用的非饱和黏土砂有效黏聚力$c' = 0\;{\rm{kPa}}$，有效内摩擦角$\varphi ' = 30$°，基质吸力$s = 200\;{\rm{kPa}}$，饱和度${s_{\text{r}}} = 37.5\text{%}$，平均净主应力$p' = 100\;{\rm{kPa}}$。图2为本文新建立的非饱和土广义有效应力三剪强度准则计算结果与文献[29-30]中的真三轴试验数据对比验证图。

由图3可知，将本文所推导的破坏准则预测结果与单应力变量法和双应力变量法的三剪强度准则预测结果相比较，可以得到广义有效应力变量法最为吻合试验结果，同时和双应力变量法比较接近。单应力变量法预测误差更大。新准则虽未对双应力变量法三剪强度准则产生明显的优势，但其采用了广义有效应力原理，能更全面地描述非饱和土的力学特性。

图  3  广义有效应力变量法计算结果和试验数据对比图

Figure  3.  Comparison diagram of calculation results of generalized effective stress variable method and test data

下载: 全尺寸图片 幻灯片

由图4可知，将广义有效应力变量法预测的$p'$-$q$子午线与试验点相比对，文献中的试验数据点比较少，试验点基本与$\theta = {0}$°的广义有效应力变量法下的$p'$-$q$子午线吻合，新准则取得很好的预测效果。

图  4  子午线计算结果与试验数据对比图

Figure  4.  Comparison diagram of meridian calculation results and experimental data

下载: 全尺寸图片 幻灯片

2   非饱和黏性土三剪统一弹塑性本构模型

在弹塑性力学理论框架下，将第1节所推导的破坏应力比引入到非饱和土修正剑桥模型中，建立基于广义有效应力变量法的非饱和土三剪统一弹塑性本构模型。采用相关联的流动法则，由一致性条件和硬化规律推导出本构模型中的塑性模量和弹塑性刚度矩阵各元素的具体表达式。

2.1   屈服函数

饱和土修正剑桥模型表达式为：

式中：$f$为屈服函数；${p'_x}$为饱和土的有效屈服应力；$p'$为饱和土的平均有效主应力；$q$为偏应力。

用非饱和土的净屈服应力对饱和土修正剑桥模型中的有效屈服应力进行替换得到非饱和土修正剑桥模型为：

式中：$p'$为非饱和土净平均主应力；$M$为破坏应力比；${p'_y}\left( s \right)$为非饱和土的屈服应力。

文献[31]曾提出过非饱和土的LC屈服曲线，但其中土的回弹指数为一定值，没有考虑基质吸力的影响。文献[32]通过试验研究发现回弹指数并非是一定值，而是随基质吸力变化的量，并将LC屈服曲线改写为：

式中：${p'_y}\left( s \right)$为非饱和土的净屈服应力；${p'_y}\left( 0 \right)$为饱和土的有效(净)屈服应力；${p'_{\rm{n}}}$为吸湿不发生湿化变形的应力；$\lambda \left( 0 \right)$为饱和土的压缩指数；$\kappa \left( 0 \right)$为饱和土的回弹指数；$\lambda \left( s \right)$为基质吸力为s的非饱和土的压缩指数；$\kappa \left( s \right)$为基质吸力为s的非饱和土的回弹指数。由文献[32]可知，非饱和土压缩指数和回弹指数与基质吸力间的关系为：

式中：${\lambda _{\rm{s}}}$为表达非饱和土压缩指数随基质吸力大小变化的参数；${\kappa _{\rm{s}}}$为表达非饱和土回弹指数随基质吸力大小变化的参数；${p_{\text{1}}}$为大气压。

将广义有效应力变量法得到的非饱和土三剪统一破坏比代入到非饱和土的修正剑桥模型中，得到$p'$-$q$子午面上的非饱和土屈服线函数为：

饱和土的应变硬化函数为：

式中：${p'_0}$为初始平均净主应力；$\varepsilon _{\rm{v}}^{\rm{p}}$为塑性体应变。

将式(17)代入式(16)，将式(16)写成椭圆方程形式，得到基于广义有效应力的非饱和土屈服函数表达式如下：

1) 采用等量代换法时，

2) 采用坐标平移法时，

2.2   流动法则

采用相关联的流动的流动法则，屈服函数与塑性势函数相同，且塑性流动方向与加载方向一致。

式中，${L_{ij}}$为加载方向。

2.3   塑性模量

塑性应力和应变之间的增量关系为：

式中：${\rm{d}}\varepsilon _{\rm{v}}^{\rm{p}}$为塑性体应变增量；${\rm{d}}\varepsilon _{\rm{q}}^{\rm{p}}$为塑性剪应变增量。

加载指数和塑性模量之间的关系为：

式中，${K_{\rm{p}}}$为塑性模量。

由一致性条件可得：

式中：${\rm{d}}{\rm{}}\varepsilon _{\rm{v}}^{\rm{p}}$为塑性体应变增量；${\rm{d}}p' = \dfrac{{\partial p'}}{{\partial {\sigma '_{ij}}}}{\rm{d}}{\sigma '_{ij}}$；${\rm{d}}q = \dfrac{{\partial q}}{{\partial {\sigma '_{ij}}}}{\rm{d}}{\sigma '_{ij}}$。

由式(21)和式(22)得：

将式(24)代入式(23)得：

根据复合函数求导法则得：

将式(26)代入式(25)得：

由式(27)求得塑性模量的表达式为：

其中：

1) 使用等量代换法时，

2) 使用坐标平移法时，

2.4   弹塑性本构关系

由弹塑性应变增量可知：

式中：${\rm{d}}{\varepsilon _{ij}}$为应变增量；${\rm{d}}\varepsilon _{ij}^{\rm{e}}$为弹性应变增量；${\rm{d}}\varepsilon _{ij}^{\rm{p}}$为塑性应变增量。

又因为：

式中：${\rm{d}}{\sigma '_{ij}}$为净应力增量；${{\boldsymbol{D}}^{\rm{e}}}$为弹性刚度矩阵。

将式(21)进行变换得式(32)：

将式(30)和式(31)代入式(32)得：

由此得到加载指数的表达式为：

将式(31)和式(34)代入式(30)得：

式中，${{\boldsymbol{D}}^{{\rm{ep}}}}$为弹塑性刚度矩阵。

从而得到弹塑性刚度矩阵为：

式(36)中的具体表达式见附录。

3   正常固结非饱和江西红黏土试验

采用江西南昌地区非饱和重塑红黏土，通过室内土工试验得到试验土样的相对重度、最大干密度、最优含水率、液限和塑限等土工参数。制备压实度分别为80%、85%和90%三种土样，在不同基质吸力条件下进行压缩回弹试验，得到本文试验土样的压缩指数和回弹指数。通过土水特征曲线试验得到本文试验土样的土水特征曲线方程，由土水特征曲线方程计算饱和度结果为：初始压实度为90%的条件下，基质吸力等于100 kPa和300 kPa时所对应的土体的饱和度分别为83.9%和69.9%；初始压实度为85%的条件下，基质吸力等于50 kPa和100 kPa时所对应的土体的饱和度分别为79.3%和72.3%；初始压实度为80%的条件下，基质吸力等于100 kPa和200 kPa时所对应的土体的饱和度分别为63.7%和51.2%。最后在压实度分别为80%、85%和90%；净围压分别为100 kPa、200 kPa和300 kPa；不同基质吸力条件下进行非饱和土固结排水试验。由试验结果得到非饱和土的有效内摩擦角、有效黏聚力和偏应力和轴向应变、体应变和轴向应变间的关系曲线。将上述试验所得模型计算参数汇总表如表1所示。

表  1  模型计算参数汇总表

Table  1.  Summary table of model calculation parameters

压实度/(%)	饱和土 压缩指数$\lambda \left( 0 \right)$	非饱和土 材料参数${\lambda _{\rm{s}}}$	饱和土 回弹指数$\kappa \left( 0 \right)$	非饱和土 材料参数${\kappa _{\rm{s}}}$	饱和土 有效黏聚力$c'$/kPa	饱和土 有效内摩擦角$\varphi '$/(o)	泊松比$\nu$	初始 孔隙比${e_0}$
90	0.0666	0.01930	0.00639	−2.640×10−6	26.90	31	0.35	0.56
85	0.1195	0.06602	0.01369	−2.767×10−5	26.76	29	0.35	0.71
80	0.1285	0.04024	0.01588	−2.965×10−5	25.10	26	0.35	0.81

下载: 导出CSV 

| 显示表格

4   本构模型的常规三轴固结排水试验验证

为验证新建立的弹塑性本构模型的正确性，采用Fortran语言编制本文模型的计算程序，对非饱和土固结排水试验进行数值模拟。将等量代换法和坐标平移法下的本构模型计算结果与非饱和土三轴固结排水试验结果进行对比，进一步验证新建立的本构模型的正确性。对比结果如图5和图6所示。

图  5  剪应力模拟值与试验值对比图

Figure  5.  Comparison between simulated shear stress and experimental stress

下载: 全尺寸图片 幻灯片

图  6  小变形体应变模拟值与试验值对比图

Figure  6.  Comparison between simulated and experimental strain values of small deformed body

下载: 全尺寸图片 幻灯片

由图5(a)可知，在广义有效应力变量法下，等量代换法与试验数据差值最大不超过8%，坐标平移法与试验数据差值最大不超过12%；由图5(b)可知，等量代换法与试验数据差值最大不超过10%，坐标平移法与试验数据差值最大不超过15%；由图5(c)可知，等量代换法与试验数据差值最大不超过9%，坐标平移法与试验数据差值最大不超过14%；由图5(d)可知，等量代换法与试验数据差值最大不超过14%，坐标平移法与试验数据差值最大不超过18%。总体来看，两种方法模拟数据与试验结果均符合较好。不同参数下所有的偏应力模拟值与试验值对比图有着类似的变化规律。本构模型模拟曲线和试验点的变化规律是一致的，说明该本构模型能够较好地反映土体在固结排水状态下的强度特性。同时本构模型计算值与试验值比较接近，对试验结果取得了较好的模拟预测，进一步验证了该模型的正确性。等量代换法模拟值和坐标平移法模拟值在轴向应变较小时比较接近，随着轴应变的增大，二者之间的差值逐渐增大，最终等量代换法计算结果大于坐标平移法计算结果。在轴向应变超过4%以后，坐标平移法的模拟效果相对较差。采用等量代换法计算结果比坐标平移法结果更加精确。土样在剪切的过程中，轴向应变较小时，偏应力随着轴向变形的增大而迅速增大，直至后期逐渐达到稳定值后缓慢增长。在初始压实度和基质吸力相同的条件下，净围压越大，土样的抗剪强度越大；在初始压实度和净围压相同的条件下，基质吸力增大，土样抗剪强度增大；在净围压和基质吸力相同的情况下，初始压实度越大，土样的抗剪强度也会越大。导致计算结果产生误差的原因可能与试验结果误差有关，也有可能与模型中的土样参数值存在误差有关。

由图6可知，本构模型模拟曲线和试验点的变化规律是一致的，说明该本构模型能够较好地反映体应变和轴向应变之间的变化规律。同时该模型计算的体应变值接近于试验值，对试验取得了较好的模拟预测结果，进一步验证了该模型的正确性。在轴向应变较小时，等量代换法计算值和坐标平移法计算值比较接近，在体应变逐渐达到平稳状态时，最终等量代换法计算值大于坐标平移法计算值，等量代换法计算值更符合试验值结果。在土体剪切过程中，体应变在轴向变形较小时增长迅速，直至后期逐渐接近平稳状态时增长缓慢。在初始压实度和基质吸力相同的条件下，净围压越大，土样达到稳定时的体应变值越大；在净围压和基质吸力相同的情况下，初始压实度越大，土样达到稳定时的体应变值越小。

5   固结排水条件下真三轴试验数值模拟

为了反映真实的土体受力情况，本节用Fortran语言根据本文所建立的本构模型，采用等量代换法和坐标平移法对长方体试样进行真三轴非饱和固结排水试验条件下的数值模拟。本节所采用的参数仍旧是江西红黏土的相关土性参数，试验应力加载路径为先对试样三个主应力方向上都施加最小主应力固结，等第一阶段固结完成以后逐步增大其中一个水平方向上的主应力增大至中间主应力大小，直至固结完成。固结阶段完成以后保持中间主应力大小和最小主应力大小不变，逐步增大轴向应力，使试样开始剪切，直至产生15%轴向应变发生剪切破坏。其中，中间主应力影响系数分别取0、0.25、0.50、0.75、1.00，最小主应力取200 kPa，最大主应力取300 kPa。真三轴试样尺寸为：$140\;{\rm{mm}} \times 70\;{\rm{mm}} \times 70\;{\rm{mm}}$。在真三轴固结排水试验条件下，本构模型模拟所得到的偏应力与轴向应变关系曲线、体应变与轴向应变关系曲线对比图如图7~图10所示。

图  7  压实度90%、基质吸力200 kPa偏应力和轴向应变关系图

Figure  7.  Relationship between the 200 kPa deviant stress and axial strain of the matric suction with 90% compactness

下载: 全尺寸图片 幻灯片

图  8  压实度85%、基质吸力200 kPa偏应力和轴向应变关系图

Figure  8.  Relationship between the 200 kPa deviant stress and axial strain of the matric suction with 85% compactness

下载: 全尺寸图片 幻灯片

图  9  压实度85%、基质吸力200 kPa体应变和轴向应变关系图

Figure  9.  Relationship between 200 kPa volume strain and axial strain of 85% compacted matric suction

下载: 全尺寸图片 幻灯片

图  10  压实度80%、基质吸力200 kPa体应变和轴向应变关系图

Figure  10.  Relationship between 200 kPa volume strain and axial strain of 80% compacted matric suction

下载: 全尺寸图片 幻灯片

由图7~图10可知，在真三轴固结排水条件下，中间主应力影响系数对偏应力和体应变都产生影响。在其他条件都相同的情况下，中主应力影响系数值越大，所对应的偏应力和体应变就越大，表明中主应力影响系数可以提高土体的抗剪强度。剪切初期，在轴向应变较小时，不同b值模拟的曲线比较接近，此时偏应力和体应变增长较快；随着轴向应变的增加，不同b值计算的偏应力和体应变值间的差值逐渐增大，曲线增长速度减缓逐渐达到平稳状态。

6   结论

本文将非饱和土广义有效应力原理与三剪强度准则相结合，提出了非饱和土广义有效应力三剪强度准则，并推导出了能够反映土体全应力状态的破坏应力比$M\left( {p',\theta ,s,{s_{\rm{r}}}} \right)$和$M\left( \theta \right)$，将新推导的破坏应力比引入到非饱和土修正剑桥模型中，建立了新的非饱和土三剪统一屈服函数。在弹塑性理论的框架下，建立了基于广义有效应力变量法的非饱和土三剪统一弹塑性本构模型。将本构模型计算结果与非饱和土三轴固结排水试验结果进行了对比验证。结论如下：

(1) 新的强度准则能够很好地反映土体的真实情况，有效地拓宽了非饱和土破坏强度理论，同时它也继承了三剪强度准则的优势，是一系列统一破坏准则，通过取不同的b值可对其他准则作非线性逼近。相比于原来破坏应力比为定值的修正剑桥模型中的屈服函数，新的屈服函数可以更好地反映土体全应力状态、中间主应力效应、拉压不等效应。

(2) 新的本构模型该模型适用于正常固结非饱和黏性土，能够体现中间主应力、黏聚力、基质吸力和饱和度对非饱和土强度和变形特性的影响。可以反映单调静荷载条件下非饱和土体的强度特性、变形特性、屈服特性、压硬性、剪缩特性。

(3) 将数值计算结果与固结排水试验结果进行了对比验证，结果表明：模拟曲线和试验点的变化规律是一致的，且计算值与试验值比较接近，说明新的本构模型能够较好地反映土体在固结排水状态下的强度和变形特性，验证了新的本构模型的正确性。在轴向变形较小时，等量代换法和坐标平移法模拟结果比较接近，随着轴向变形逐渐增大直至剪应力达到平稳状态的过程中，等量代换法计算结果要大于坐标平移法模型计算结果，且更接近于试验值，说明等量代换法能更好地反映土体强度和变形特性，而坐标平移法模拟预测效果相对较差。

附录：${\;\;}$

$\small \dfrac{{\partial f}}{{\partial q}}$和$\small \dfrac{{\partial f}}{{\partial \theta }}$的计算式为：

1) 使用等量代换法时，

2) 使用坐标平移法时，

$\small {B_1}$的计算式为：

1) 使用等量代换法时，

2) 使用坐标平移法时，