超固结非饱和黏性土三剪次加载面硬化模型

Triple-shear Subloading Surface Hardening Model for Over-consolidated Unsaturated Clays

摘要

将Bishop应力变量法和Fredlund应力变量法分别应用于三剪强度准则，推导出两种应力变量法下的非饱和土三剪破坏应力比及非饱和土三剪屈服面方程，并引入统一硬化参数来反映土的超固结性，得到了超固结非饱和土的屈服面方程。在此基础上，采用次加载面理论反映土在动静荷载作用下的变形特性，并以此建立超固结非饱和黏性土三剪次加载面硬化模型，该模型能够反映超固结非饱和黏性土的应变软化和剪胀特征及循环荷载作用下塑性应变的累积和滞洄特性。以江西南昌地区重塑红黏土为研究对象进行室内土工试验、压缩回弹试验，得到本构模型计算所需基本土性参数。利用非饱和三轴试验仪进行超固结非饱和土静三轴固结排水试验，利用动态土工三轴仪进行超固结饱和土循环三轴固结排水试验，得到在静、动荷载下超固结土的应力‒应变关系曲线。将本构模型计算结果分别与超固结土静、动三轴固结排水试验结果进行对比验证，结果表明：本文本构模型计算模拟结果与试验结果变化规律一致，可以较好地反映超固结土应变软化、剪胀性，以及动荷载下超固结黏性土棘轮效应和曼辛效应等力学特性；验证了所建立本构模型在模拟超固结土应力应变等基本力学特性方面的适用性。真三轴数值模拟结果表明，中间主应力影响系数、基质吸力和超固结比越大，常规静态三轴表现为应力峰值增大，动循环三轴表现为土体达到变形结束时的累计塑性应变减小。本文本构模型能够较好地反映超固结土体在真三轴固结排水状态下的力学特性。

Abstract

The Bishop stress variable method and the Fredlund stress variable method are applied to the triple-shear strength criterion to derive two triple-shear failure stress ratios and the related triple-shear yield surface equations for unsaturated clays under the two stress variable methods. The unified hardening parameter is then introduced into the original yield surface equations to obtain the over-consolidated yield surface equations for unsaturated soil. The subloading surface theory is applied to establish two triple-shear subloading surface hardening models (TSSH) for over-consolidated unsaturated clays to reflect the soil deformation characteristics under dynamic and static loads. These two proposed models reflect strain softening, dilatancy, as well as plastic strain accumulation and hysteresis under cyclic loading. Taking the remolded red clay in the Nanchang area of Jiangxi Province as the research object, the basic soil parameters required for the calculation of the TSSH are obtained by laboratory geotechnical tests. The unsaturated static triaxial CD (consolidated drained) tests and the saturated cyclic triaxial CD tests for the over-consolidated unsaturated clays are conducted. The stress‒strain relationship curves for the over-consolidated clays under static and cyclic loads are obtained, and the constitutive model results are compared to the static and cyclic triaxial CD tests. The results showed that the constitutive model outcomes are consistent with the experimental data, which better reflect the mechanical properties such as strain softening, dilatancy, and the ratchet and Masing effects of the over-consolidated clays under cyclic loads. The applicability of the proposed constitutive models in reflecting the basic stress‒strain properties of the over-consolidated clays is verified. The numerical simulations of the true triaxial show that when the influence coefficient of the intermediate principal stress, the matrix suction, and the over-consolidation ratio are larger, the peak stress of conventional static triaxial becomes larger, and the cumulative plastic strain of the cyclic triaxial decreases as the clays reach the end of deformation. The constitutive model effectively reflects the mechanical properties of over-consolidated soil under true triaxial consolidation and drainage.

Graphical abstract

关键词

超固结非饱和黏性土 / 次加载面硬化模型 / 三剪强度准则 / 统一硬化参数 / 江西红黏土

Key words

over-consolidated unsaturated clays / subloading surface hardening model / triple-shear strength criterion / unified hardening parameters / Jiangxi red clay

本刊网刊

21 网络出版日期:2024 ‒ 03 ‒ 25

国内外学者在常规弹塑性本构模型的基础上提出了边界面模型、下加载面模型、UH模型、两面模型等可以描述超固结土力学特性的本构模型。Dafalias等^[1‒4]提出的边界面模型能够描述应力在屈服面内变化而引起的塑性变形，不用记忆多个屈服面的大小、位置及塑性模量，但是其模型参数缺乏明确物理含义且确定方法十分复杂。Hashiguchi等^[5‒6]提出了下加载面模型。该模型可以反映土体在单调、循环荷载作用下的真实力学性能。Nakai等^[7]建立了适用于超固结土的t_ij 模型，该模型可以考虑中间主应力对应力路径的依赖性和对土应力应变的影响，可以较好地反映超固结土的3维力学特性，但是t_ij 模型也存在参数物理含义不明确且不能有效地描述超固结土剪胀特性的问题。Houlsby等^[8]在土体弹塑性本构模型基础上引入Hvorslev线，建立了适用于描述超固结土特性的弹塑性模型，但是该模型描述土体在3维状态下的力学特性并不是很理想。Mita等^[9]在其基础上进行修正，建立了可以描述超固结土在3维状态下的应力应变特性的Hvorslev面修正剑桥模型，但是该模型以Hvorslev面为屈服面，却以修正剑桥模型为塑性势面，所以其合理性需要进一步验证。李潇旋等^[10‒11]构建了超固结非饱和土弹塑性双面模型，考虑了土体超固结程度和基质吸力的影响。姚仰平等^[12]提出统一硬化参数H，构建了能够描述土体剪缩、剪胀的统一硬化模型，该模型中的材料参数物理意义明确且确定方法简单，可以通过对统一硬化参数形式的适当调整，实现对松砂、密砂、黏土等具有各向异性土体的应力应变特性的描述。

上述模型多针对超固结非饱和土在静荷载作用下的力学特性进行研究，但对动荷载作用下超固结非饱和土的力学特性还少有涉及。Hashiguchi^[13‒14]提出的次加载面模型可以较好地模拟和预测土体在循环荷载作用下所表现出的曼辛效应和棘轮效应等力学特性，同时克服了边界面模型中采用不严谨的试凑法来确定模型参数的缺陷，其模拟预测的土体在静、动荷载作用下的应力‒应变曲线连续、光滑。

目前，大多数的超固结非饱和土本构模型都是在临界状态土力学的理论框架下建立起来的，然而，屈服函数中的破坏应力比需要强度准则来确定。土体的强度准则主要有Mohr‒Coulomb准则、Druker‒Prager准则、Lade准则、SMP准则和双剪统一强度准则。Mohr‒Coulomb准则没有考虑中间主应力影响，Druker‒Prager准则没有反映出拉压差效应，Lade准则不能考虑黏聚力对土体破坏的影响，SMP准则因其模型特征单一而不能较好地反映土体力学特性的多样性。基于双剪统一强度准则^[15]，同时考虑十二面体单元上3个主剪应力共同作用所提出的三剪统一强度准则，在反映岩土材料的中间主应力效应特性、拉压不等效应特性及应力区间效应特性等强度特性方面具有自身优点^[16]。三剪统一准则已被应用于结构性土^[17]、含泥饱和砂土^[18]和正常固结非饱和黏性土的力学特性分析中^[19]。

综上所述，超固结非饱和黏性土力学特性较为复杂，常规临界状态模型无法合理地描述超固结非饱和土体应变软化、剪胀等特性；并且，已有的对于超固结土体在动态荷载作用下本构模型的研究较少，且大多数超固结土本构模型未考虑中间主应力效应影响。因此，本文基于非饱和土修正剑桥模型，引入统一硬化参数、次加载面理论和非饱和土三剪破坏应力比，建立超固结非饱和黏性土的三剪次加载面模型。该模型能够反映超固结土的硬化、软化、剪胀、剪缩等特性，可以体现超固结土在循环荷载作用下的棘轮效应（塑性应变的累积特性）和曼辛效应（滞回特性），还可以考虑中间主应力效应，能够反映超固结非饱和黏性土在更为复杂的真三轴状态下的应力应变特性。

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

1.1 三剪统一强度准则

三剪统一强度准则表达式为^[16‒18]：

\begin{array}{l}\left[({\overline{\sigma }}_{\mathrm{1}}^{\mathrm{\text{'}}}-{\overline{\sigma }}_{\mathrm{3}}^{\mathrm{\text{'}}}{)}^{\mathrm{2}}+b({\overline{\sigma }}_{\mathrm{1}}^{\mathrm{\text{'}}}-{\overline{\sigma }}_{\mathrm{2}}^{\mathrm{\text{'}}}{)}^{\mathrm{2}}+b({\overline{\sigma }}_{\mathrm{2}}^{\mathrm{\text{'}}}-{\overline{\sigma }}_{\mathrm{3}}^{\mathrm{\text{'}}}{)}^{\mathrm{2}}\right]-(\mathrm{1}+b)\cdot \\ \left[({\overline{\sigma }}_{\mathrm{1}}^{\mathrm{\text{'}}}{)}^{\mathrm{2}}-({\overline{\sigma }}_{\mathrm{3}}^{\mathrm{\text{'}}}{)}^{\mathrm{2}}\right]\mathrm{s}\mathrm{i}\mathrm{n}{\phi }^{\text{'}}=\mathrm{2}{c}^{\text{'}}(\mathrm{1}+b)({\overline{\sigma }}_{\mathrm{1}}^{\mathrm{\text{'}}}-{\overline{\sigma }}_{\mathrm{3}}^{\mathrm{\text{'}}})\mathrm{c}\mathrm{o}\mathrm{s}{\phi }^{\text{'}}\end{array}（1）

式中，{\overline{\sigma }}_{\mathrm{1}}^{\mathrm{\text{'}}}、{\overline{\sigma }}_{\mathrm{2}}^{\mathrm{\text{'}}}、{\overline{\sigma }}_{\mathrm{3}}^{\mathrm{\text{'}}}分别为最大、中间和最小有效主应力，b为中间主应力影响系数，c^{\text{'}}为土体的有效黏聚力，{\phi }^{\text{'}}为土体的有效内摩擦角。

文献[20]分别将Bishop应力变量法和Fredlund应力变量法与三剪统一强度准则相结合，提出了Bishop应力变量法和Fredlund应力变量法的非饱和土三剪强度准则。

Bishop应力变量法的非饱和土三剪统一强度准则表达式为^[20]：

q=A(p^{\text{'}}+\chi s)\mathrm{s}\mathrm{i}\mathrm{n}{\phi }^{\text{'}}+A{c}^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }^{\text{'}}（2）

式中：q为广义剪应力；p^{\text{'}}为平均净主应力；\chi为Bishop系数；s为基质吸力；A的表达式为

\begin{array}{l}A=\left[\mathrm{6}(\mathrm{1}+b)\mathrm{c}\mathrm{o}\mathrm{s}\left(\theta -\frac{\mathrm{\pi }}{\mathrm{6}}\right)\right]/\left\{\mathrm{2}\sqrt[]{\mathrm{3}}\left[\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\left(\theta -\frac{\mathrm{\pi }}{\mathrm{6}}\right)+\right.\right.\\ \left.\left.b\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\left(\theta +\frac{\mathrm{\pi }}{\mathrm{6}}\right)+b\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\theta \right]-(\mathrm{1}+b)\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{2}\theta +\frac{\mathrm{\pi }}{\mathrm{6}}\right)\mathrm{s}\mathrm{i}\mathrm{n}{\phi }^{\text{'}}\right\}\end{array}，

其中，\theta为应力角。

Fredlund应力变量法的非饱和土三剪统一强度准则表达式为^[20]：

q=A{p}^{\text{'}}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }^{\text{'}}+A(c^{\text{'}}+c_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{s}{\phi }^{\text{'}}（3）

式中，c_{\mathrm{s}}为非饱和土的吸附黏聚力，其表达式为c_{\mathrm{s}}=s\cdot \mathrm{t}\mathrm{a}\mathrm{n}{\phi }^{\mathrm{b}}，φ^b为吸力角。

1.2 破坏应力比

将Bishop应力变量法的非饱和土三剪统一强度准则作为破坏准则，采用坐标平移法推导出新的破坏应力比。由文献[20]可知，非饱和土Bishop应力变量法坐标平移后的破坏应力比{\widehat{M}}_{\mathrm{1}}为：

{\widehat{M}}_{\mathrm{1}}=\frac{\widehat{q}}{{\widehat{p}}^{\mathrm{\text{'}}}}=\frac{q}{p^{\text{'}}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s}=A\mathrm{s}\mathrm{i}\mathrm{n}{\phi }^{\text{'}}（4）

式中，{\widehat{p}}^{\mathrm{\text{'}}}、\widehat{q}分别为坐标平移后的平均净主应力和偏应力。

平移前后土的压缩指数\lambda、\widehat{\lambda }的关系为：

\widehat{\lambda }={\Lambda }_{\mathrm{1}}\lambda（5）

式中，{\Lambda }_{\mathrm{1}}=\frac{\mathrm{l}\mathrm{n}(c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s)}{\mathrm{l}\mathrm{n}\left[\frac{\mathrm{2}(c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s)}{\mathrm{1}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s}\right]}。

平移前后土的回弹指数\kappa \mathrm{、}\widehat{\kappa }间的关系为：

\widehat{\kappa }={\Lambda }_{\mathrm{1}}\kappa（6）

与Bishop应力变量法类似，Fredlund应力变量法中坐标平移后的非饱和土破坏应力比{\widehat{M}}_{\mathrm{2}}为：

{\widehat{M}}_{\mathrm{2}}=\frac{\widehat{q}}{{\widehat{p}}^{\mathrm{\text{'}}}}=\frac{q}{p^{\text{'}}+(c^{\text{'}}+c_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}}=A\mathrm{s}\mathrm{i}\mathrm{n}{\phi }^{\text{'}}（7）

采用坐标平移法时Fredlund应力变量法中的比例因子为：

{\Lambda }_{\mathrm{2}}=\frac{\mathrm{l}\mathrm{n}\left[(c^{\text{'}}+c_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}\right]}{\mathrm{l}\mathrm{n}\left[\frac{\mathrm{2}(c^{\text{'}}+c_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}}{\mathrm{1}+(c^{\text{'}}+c_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}}\right]}（8）

2 超固结非饱和黏性土三剪次加载面模型

2.1 三剪正常屈服面方程

饱和土的修正剑桥模型屈服面方程为：

f=q^{\mathrm{2}}+M^{\mathrm{2}}{p^{\text{'}}}^{\mathrm{2}}-M^{\mathrm{2}}{p^{\text{'}}}^{\mathrm{2}}{p}_{\mathrm{y}}^{\text{'}}(\mathrm{0})=\mathrm{0}（9）

式中，f 为屈服函数，p_{\mathrm{y}}^{\text{'}}(\mathrm{0})为有效屈服应力，M为破坏应力比。

采用孙德安^[21]提出的LC屈服线，并根据其修正后的LC屈服线得到非饱和状态时土的净屈服应力，其表达式为：

p_{\mathrm{y}}^{\text{'}}\left(s\right)=p_{\mathrm{n}}^{\text{'}}{\left(\frac{p_{\mathrm{y}}^{\text{'}}\left(\mathrm{0}\right)}{p_{\mathrm{n}}^{\text{'}}}\right)}^{\frac{\lambda (\mathrm{0})-\kappa (\mathrm{0})}{\lambda (s)-\kappa (s)}}（10）

式中，p_{\mathrm{y}}^{\text{'}}\left(s\right)为非饱和土的净屈服应力，p_{\mathrm{n}}^{\text{'}}为当非饱和土吸湿时不发生湿化变形的应力，\lambda (\mathrm{0})和\kappa (\mathrm{0})分别为饱和土的压缩指数和回弹指数，\lambda (s)和\kappa (s)分别为非饱和土的压缩指数和回弹指数。

由文献[20]可知，非饱和土压缩指数\lambda (s)、回弹指数\kappa (s)与基质吸力s之间的关系为：

\left\{\begin{array}{l}\lambda (s)=\lambda (\mathrm{0})-\frac{{\lambda }_{\mathrm{s}}s}{p_{\mathrm{a}}+s},\\ \kappa (s)=\kappa (\mathrm{0})+{\kappa }_{\mathrm{s}}s\end{array}\right.（11）

式中，{\lambda }_{\mathrm{s}}为反映非饱和土压缩指数\lambda (s)随基质吸力大小变化的参数，{\kappa }_{\mathrm{s}}为反映非饱和土回弹指数\kappa (s)随着基质吸力大小变化的参数，p_{\mathrm{a}}代表标准大气压。

用p_{\mathrm{y}}^{\text{'}}(s)对饱和土修正剑桥模型中的有效屈服应力p_{\mathrm{y}}^{\text{'}}(\mathrm{0})进行替换，即将式（10）代入式（9）中，可得非饱和土修正剑桥模型中的屈服函数为：

f=q^{\mathrm{2}}+M^{\mathrm{2}}{p^{\text{'}}}^{\mathrm{2}}-M^{\mathrm{2}}{p}^{\text{'}}{p}_{\mathrm{n}}^{\text{'}}{\left(\frac{p_{\mathrm{y}}^{\text{'}}(\mathrm{0})}{p_{\mathrm{n}}^{\text{'}}}\right)}^{\frac{\lambda (\mathrm{0})-\kappa (\mathrm{0})}{\lambda (s)-\kappa (s)}}=\mathrm{0}（12）

土在饱和状态时的硬化函数为：

\mathrm{l}\mathrm{n}\left[\frac{p_{\mathrm{y}}^{\text{'}}(\mathrm{0})}{p_{\mathrm{0}}^{\text{'}}}\right]=\frac{\mathrm{1}+e_{\mathrm{0}}}{\lambda (\mathrm{0})-\kappa (\mathrm{0})}{\epsilon }_{\mathrm{v}}^{\mathrm{p}}（13）

式中，p_{\mathrm{0}}^{\text{'}}为初始平均净主应力，e_{\mathrm{0}}为土的初始孔隙比，{\epsilon }_{\mathrm{v}}^{\mathrm{p}}为土的塑性体应变。

采用式（4）中的{\widehat{M}}_{\mathrm{1}}和式（7）中的{\widehat{M}}_{\mathrm{2}}分别对式（12）中的M进行替换，并结合式（13）中的硬化函数，可得在Bishop应力变量法和Fredlund应力变量法中非饱和土三剪正常屈服面方程分别为：

f=(p^{\text{'}}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s)+\frac{q^{\mathrm{2}}}{{\widehat{M}}_{\mathrm{1}}^{\mathrm{2}}(p^{\text{'}}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s)}-p_{x\mathrm{1}}（14）

f=(p^{\text{'}}+(c^{\text{'}}+c_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}})+\frac{q^{\mathrm{2}}}{{\widehat{M}}_{\mathrm{2}}^{\mathrm{2}}(p^{\text{'}}+(c^{\text{'}}+c_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}})}-p_{x\mathrm{2}}（15）

式（14）～（15）中，p_{x\mathrm{1}}和p_{x\mathrm{2}}分别为Bishop应力变量法和Fredlund应力变量法中的等效平均净应力，即非饱和土三剪正常屈服面与p′ 轴的交点。

正常屈服面定义为^[22]：

f({\sigma }_{\mathrm{y}})=F(\overline{H})（16）

式中，{\sigma }_{\mathrm{y}}为当前应力状态\sigma在正常屈服面上的对偶应力，F为加载屈服面，\overline{H}为各向同性硬化参量。

引入文献[23]中的修正统一硬化参数H，使本构模型能够反映非饱和土在超固结状态下的应力软化特性和剪胀特性：

H=\int \frac{M_{\mathrm{f}}^{\mathrm{4}}-{\eta }^{\mathrm{4}}}{M^{\mathrm{4}}-{\eta }^{\mathrm{4}}}\mathrm{d}{\epsilon }_{\mathrm{v}}^{\mathrm{p}}（17）

式中：M_f为潜在强度；\eta为应力比，\eta =q/p^{\text{'}}。

根据式（14）～（17），可得超固结非饱和黏性土在Bishop应力变量法和Fredlund应力变量法中的三剪正常屈服面方程分别为：

\left\{\begin{array}{l}f(p^{\text{'}},q,\theta ,s)=(p^{\text{'}}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s)+\frac{q^{\mathrm{2}}}{{\widehat{M}}_{\mathrm{1}}^{\mathrm{2}}(p^{\text{'}}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s)},\\ F(H)=p_{x\mathrm{1}}=(p_{\mathrm{n}}^{\text{'}}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s){\left[\frac{p_{\mathrm{0}}^{\text{'}}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s}{p_{\mathrm{n}}^{\text{'}}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s}\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{\mathrm{1}+e_{\mathrm{0}}}{{\Lambda }_{\mathrm{1}}\left[\lambda (\mathrm{0})-\kappa (\mathrm{0})\right]}H\right)\right]}^{\frac{\lambda (\mathrm{0})-\kappa (\mathrm{0})}{\lambda (s)-\kappa (s)}}\end{array}\right.（18）

\left\{\begin{array}{l}f(p^{\text{'}},q,\theta ,s)=\left[p^{\text{'}}+(c^{\text{'}}+c_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}\right]+\frac{q^{\mathrm{2}}}{{\widehat{M}}_{\mathrm{2}}^{\mathrm{2}}(p^{\text{'}}+(c^{\text{'}}+c_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}})},\\ F(H)=p_{x\mathrm{2}}=\left[p_{\mathrm{n}}^{\text{'}}+(c^{\text{'}}+c_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}\right]{\left[\frac{p_{\mathrm{0}}^{\text{'}}+(c^{\text{'}}+c_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}}{p_{\mathrm{n}}^{\text{'}}+(c^{\text{'}}+c_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}}\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{\mathrm{1}+e_{\mathrm{0}}}{{\Lambda }_{\mathrm{2}}\left[\lambda (\mathrm{0})-\kappa (\mathrm{0})\right]}H\right)\right]}^{\frac{\lambda (\mathrm{0})-\kappa (\mathrm{0})}{\lambda (s)-\kappa (s)}}\end{array}\right.（19）

2.2 超固结非饱和黏性土的三剪次加载面方程

次加载面理论假定应力空间中在正常屈服面（传统屈服面）的内部存在一个与其几何相似的次加载面，其定义为^[13]：

f(\overline{\sigma })=RF(H)（20）

式中，\overline{\sigma }为次加载面上的当前应力，R为次加载面和正常屈服面尺寸的比值（相似比）。将式（18）、（19）分别代入式（20）可得超固结非饱和黏性土在Bishop应力变量法和Fredlund应力变量法中的三剪次加载面方程分别为：

\left\{\begin{array}{l}f=(p^{\text{'}}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s-{\overline{\alpha }}_{\mathrm{p}})+\frac{(q-{\overline{\alpha }}_{\mathrm{q}}{)}^{\mathrm{2}}}{{\widehat{M}}_{\mathrm{1}}^{\mathrm{2}}(p^{\text{'}}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s-{\overline{\alpha }}_{\mathrm{p}})}-R{p}_{x\mathrm{1}},\\ p_{x\mathrm{1}}=(p_{\mathrm{n}}^{\text{'}}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s){\left[\frac{p_{\mathrm{0}}^{\text{'}}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s}{p_{\mathrm{n}}^{\text{'}}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s}\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{\mathrm{1}+e_{\mathrm{0}}}{{\Lambda }_{\mathrm{1}}\left[\lambda (\mathrm{0})-\kappa (\mathrm{0})\right]}H\right)\right]}^{\frac{\lambda (\mathrm{0})-\kappa (\mathrm{0})}{\lambda (s)-\kappa (s)}}\end{array}\right.（21）

\left\{\begin{array}{l}f=(p^{\text{'}}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s-{\overline{\alpha }}_{\mathrm{p}})+\frac{(q-{\overline{\alpha }}_{\mathrm{q}}{)}^{\mathrm{2}}}{{\widehat{M}}_{\mathrm{1}}^{\mathrm{2}}(p^{\text{'}}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s-{\overline{\alpha }}_{\mathrm{p}})}-R{p}_{x\mathrm{2}},\\ p_{x\mathrm{2}}=\left[p_{\mathrm{n}}^{\text{'}}+(c^{\text{'}}+c_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}\right]{\left[\frac{p_{\mathrm{0}}^{\text{'}}+(c^{\text{'}}+c_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}}{p_{\mathrm{n}}^{\text{'}}+(c^{\text{'}}+c_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}}\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{\mathrm{1}+e_{\mathrm{0}}}{{\Lambda }_{\mathrm{2}}\left[\lambda (\mathrm{0})-\kappa (\mathrm{0})\right]}H\right)\right]}^{\frac{\lambda (\mathrm{0})-\kappa (\mathrm{0})}{\lambda (s)-\kappa (s)}}\end{array}\right.（22）

式（21）～（22）中，({\overline{\alpha }}_{\mathrm{p}},{\overline{\alpha }}_{\mathrm{q}})为三剪次加载面的几何中心。

相似比R的微分dR采用Hashiguchi^[13]提出的方程：

\mathrm{d}R=U‖\mathrm{d}{\epsilon }^{\mathrm{p}}‖（23）

式中：{\epsilon }^{\mathrm{p}}为塑性应变；U为单调递减的函数，假定为：^[13]

U=-u\mathrm{l}\mathrm{n}R（24）

将式（23）和（24）代入到式（17）中，可得超固结下相似比R的演化规律为：

\mathrm{d}R=U‖\mathrm{d}{\epsilon }^{\mathrm{p}}‖=-u\mathrm{l}\mathrm{n}R‖\mathrm{d}{\epsilon }_{\mathrm{v}}^{\mathrm{p}}+\mathrm{d}{\epsilon }_{\mathrm{q}}^{\mathrm{p}}‖=-u\mathrm{l}\mathrm{n}R‖\mathrm{d}H+\mathrm{d}{\epsilon }_{\mathrm{q}}^{\mathrm{p}}‖（25）

式（24）～（25）中：u为材料常数，其作用是控制应力点向正常屈服面靠近的速率；\mathrm{d}{\epsilon }^{\mathrm{p}}为塑性应变增量；\mathrm{d}{\epsilon }_{\mathrm{v}}^{\mathrm{p}}为塑性体应变增量；\mathrm{d}{\epsilon }_{\mathrm{q}}^{\mathrm{p}}为塑性剪切应变增量。

2.3 相似中心演化规律

相似中心面方程为^[24] ：

f(\widehat{S})=R_{\mathrm{s}}F(H)（26）

式中，\widehat{S}为考虑正常屈服面背应力的相似中心，R_s为相似中心面和正常屈服面尺寸的比值。

将式（18）、（19）分别代入式（26），可得Bishop应力变量法和Fredlund应力变量法的三剪相似中心面方程分别为：

f=(p^{\text{'}}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s)+\frac{q^{\mathrm{2}}}{{\widehat{M}}_{\mathrm{1}}^{\mathrm{2}}(p^{\text{'}}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s)}-p_{x\mathrm{1}}（27）

f=\left[p^{\text{'}}+(c^{\text{'}}+c_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}\right]+\frac{q^{\mathrm{2}}}{{\widehat{M}}_{\mathrm{2}}^{\mathrm{2}}\left[p^{\text{'}}+(c^{\text{'}}+c_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}\right]}-p_{x\mathrm{2}}（28）

相似中心的封闭条件为^[13]：

\mathrm{d}S-\mathrm{d}\alpha -\frac{\mathrm{d}F}{F}\widehat{S}\le \mathrm{0}（29）

式中，S为相似中心，\alpha为正常屈服面的几何中心。

根据相似中心的封闭条件并结合修正统一硬化参数H，可得在超固结状态下相似中心的演化规律为：

\mathrm{d}S=C‖\mathrm{d}{\epsilon }^{\mathrm{p}}‖\left(\frac{\sigma -\overline{\alpha }}{R}-\frac{S-\alpha }{w}\right)+\frac{\mathrm{1}+e_{\mathrm{0}}}{\Lambda (\lambda -\kappa )}(S-\alpha )\mathrm{d}{\epsilon }_{\mathrm{v}}^{\mathrm{p}}（30）

式中：C为次加载面的材料常数，C\ge0；w为次加载面的材料常数，0 \lew＜1，且为相似中心比R_s的最大值；\Lambda取{\Lambda }_{\mathrm{1}}时为使用Bishop应力变量法，取{\Lambda }_{\mathrm{2}}时为使用Fredlund应力变量法。

2.4 模型加卸载准则及演化机理

三剪次加载面模型的加卸载准则如下^[13]。

加载（\mathrm{d}{\epsilon }^{\mathrm{p}}\ne \mathrm{0}）：

\frac{\partial f}{\partial {\overline{\mathit{\sigma }}}^{\text{'}}}\mathrm{d}{\mathit{\sigma }}^{\text{'}}>\mathrm{0}（31）

式中，\frac{\partial f}{\partial {\overline{\mathit{\sigma }}}^{\text{'}}}为次加载面对其当前应力向量{\overline{\mathit{\sigma }}}^{\text{'}}的偏导向量，\mathrm{d}{\mathit{\sigma }}^{\text{'}}为净应力增量向量。

加载时，次加载面尺寸不断接近正常屈服面，同时R应满足的如下条件：

（32）

卸载（\mathrm{d}{\epsilon }^{\mathrm{p}}=\mathrm{0}）：

\frac{\partial f}{\partial {\overline{\mathit{\sigma }}}^{\text{'}}}\mathrm{d}{\mathit{\sigma }}^{\text{'}}\le \mathrm{0}（33）

卸载时，相似比R减小，次加载面逐渐变小，该阶段只产生弹性变形。次加载面几何中心表达式为：

\overline{\alpha }=(\mathrm{1}-R)S+R\alpha（34）

\left\{\begin{array}{l}{\overline{\alpha }}_{\mathrm{p}}=(\mathrm{1}-R)S_{\mathrm{p}}+R{\alpha }_{\mathrm{p}},\\ {\overline{\alpha }}_{\mathrm{q}}=(\mathrm{1}-R)S_{\mathrm{q}}+R{\alpha }_{\mathrm{q}}\end{array}\right.（35）

式中，(S_p，S_q)为相似中心面的几何中心。

将式（20）代入式（35）中，可得模型在卸载过程中的相似比R为：

R=\frac{-K_{\mathrm{2}}\pm \sqrt[]{K_{\mathrm{2}}^{\mathrm{2}}-\mathrm{4}{K}_{\mathrm{1}}{K}_{\mathrm{3}}}}{\mathrm{2}{K}_{\mathrm{1}}}（36）

其中：对于Bishop应力变量法有

K_{\mathrm{1}}=S_{\mathrm{p}}^{\mathrm{2}}+\frac{S_{\mathrm{p}}^{\mathrm{2}}}{{\widehat{M}}_{\mathrm{1}}^{\mathrm{2}}}-p_{x\mathrm{1}}{S}_{\mathrm{p}}，

K_{\mathrm{2}}=\mathrm{2}(p^{\text{'}}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s-S_{\mathrm{p}})S_{\mathrm{p}}+\frac{\mathrm{2}(q-S_{\mathrm{q}})S_{\mathrm{q}}}{{\widehat{M}}_{\mathrm{1}}^{\mathrm{2}}}-

(p^{\text{'}}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s-S_{\mathrm{p}})p_{x\mathrm{1}}，

K_{\mathrm{3}}=(p^{\text{'}}+c^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}+\chi s-S_{\mathrm{q}}{)}^{\mathrm{2}}+\frac{(q-S_{\mathrm{q}}{)}^{\mathrm{2}}}{{\widehat{M}}_{\mathrm{1}}^{\mathrm{2}}}；

对于Fredlund应力变量法有

K_{\mathrm{1}}=S_{\mathrm{p}}^{\mathrm{2}}+\frac{S_{\mathrm{p}}^{\mathrm{2}}}{{\widehat{M}}_{\mathrm{2}}^{\mathrm{2}}}-p_{x\mathrm{2}}{S}_{\mathrm{p}}，

K_{\mathrm{2}}=\mathrm{2}(p^{\text{'}}+(c^{\text{'}}+c_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}-S_{\mathrm{p}})S_{\mathrm{p}}+\frac{\mathrm{2}(q-S_{\mathrm{q}})S_{\mathrm{q}}}{{\widehat{M}}_{\mathrm{2}}^{\mathrm{2}}}-

(p^{\text{'}}+(c^{\text{'}}+c_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}-S_{\mathrm{p}})p_{x\mathrm{2}}，

K_{\mathrm{3}}={\left[p^{\text{'}}+(c^{\text{'}}+c_{\mathrm{s}})\mathrm{c}\mathrm{o}\mathrm{t}{\phi }^{\text{'}}-S_{\mathrm{q}}\right]}^{\mathrm{2}}+\frac{(q-S_{\mathrm{q}}{)}^{\mathrm{2}}}{{\widehat{M}}_{\mathrm{2}}^{\mathrm{2}}}。

2.5 塑性模量

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

\left\{\begin{array}{l}\mathrm{d}{\epsilon }_{\mathrm{v}}^{\mathrm{p}}=〈L〉\frac{\partial f}{\partial p^{\text{'}}},\\ \mathrm{d}{\epsilon }_{\mathrm{q}}^{\mathrm{p}}=〈L〉\frac{\partial f}{\partial q}\end{array}\right.（37）

式中：L为加载指数；〈〉表示一种运算，当L>\mathrm{0}时，〈L〉=L，当L\le \mathrm{0}时，〈L〉=\mathrm{0}。

L和塑性模量K_p之间的关系为：

L=\frac{\mathrm{1}}{K_{\mathrm{p}}}·\frac{\partial f}{\partial {\mathit{\sigma }}^{\text{'}}}\mathrm{d}{\mathit{\sigma }}^{\text{'}}（38）

将式（38）代入到式（37）中可得：

\left\{\begin{array}{l}\mathrm{d}{\epsilon }_{\mathrm{v}}^{\mathrm{p}}=\frac{\mathrm{1}}{K_{\mathrm{P}}}\left(\frac{\partial f}{\partial p^{\text{'}}}\mathrm{d}{p}^{\text{'}}+\frac{\partial f}{\partial q}\mathrm{d}q\right)\frac{\partial f}{\partial p^{\text{'}}},\\ \mathrm{d}{\epsilon }_{\mathrm{q}}^{\mathrm{p}}=\frac{\mathrm{1}}{K_{\mathrm{P}}}\left(\frac{\partial f}{\partial p^{\text{'}}}\mathrm{d}{p}^{\text{'}}+\frac{\partial f}{\partial q}\mathrm{d}q\right)\frac{\partial f}{\partial q}\end{array}\right.（39）

假设正常屈服面几何中心为原点，由次加载面方程一致性条件可得：

\frac{\partial f}{\partial \overline{\mathit{\sigma }}}\mathrm{d}\overline{\mathit{\sigma }}-F\mathrm{d}R-R\mathrm{d}F=\mathrm{0}（40）

式中，\mathrm{d}\overline{\mathit{\sigma }}=\mathrm{d}\mathit{\sigma }-(\mathrm{1}-R)\mathrm{d}\mathit{S}+\mathit{S}\mathrm{d}R。

由式（40）进一步可得：

\frac{\partial f}{\partial p^{\text{'}}}\mathrm{d}{p}^{\text{'}}+\frac{\partial f}{\partial q}\mathrm{d}q+\frac{\partial f}{\partial {\overline{\alpha }}_{\mathrm{p}}}\mathrm{d}{\overline{\alpha }}_{\mathrm{p}}+\frac{\partial f}{\partial {\overline{\alpha }}_{\mathrm{q}}}\mathrm{d}{\overline{\alpha }}_{\mathrm{q}}+\frac{\partial f}{\partial R}\mathrm{d}R+\frac{\partial f}{\partial p_x}\mathrm{d}{p}_x=\mathrm{0}（41）

式中：\mathrm{d}{\overline{\alpha }}_{\mathrm{p}}=(\mathrm{1}-R)\mathrm{d}{S}_{\mathrm{p}}-S_{\mathrm{p}}\mathrm{d}R，\mathrm{d}{\overline{\alpha }}_{\mathrm{q}}=(\mathrm{1}-R)\mathrm{d}{S}_{\mathrm{q}}-S_{\mathrm{q}}\mathrm{d}R；\mathrm{d}R=-u\mathrm{l}\mathrm{n}R‖\mathrm{d}H+\mathrm{d}{\epsilon }_{\mathrm{q}}^{\mathrm{p}}‖；\mathrm{d}{p}_x=p_x·\frac{\mathrm{1}+e_{\mathrm{0}}}{\Lambda \left(\lambda (s)-\kappa (s)\right)}\mathrm{d}H，p_x 为等效平均净应力。

采用相关联流动法则，可得超固结非饱和黏性土的三剪次加载面塑性模量K_{\mathrm{p}}为：

\begin{array}{l}{K}_{\mathrm{p}}=-\frac{\partial f}{\partial p^{\text{'}}}\left[\frac{\partial f}{\partial {\overline{\alpha }}_{\mathrm{p}}}(\mathrm{1}-R)\frac{\partial S_{\mathrm{p}}}{\partial {\epsilon }_{\mathrm{v}}^{\mathrm{p}}}+\frac{\partial f}{\partial {\overline{\alpha }}_{\mathrm{q}}}(\mathrm{1}-R)\frac{\partial S_{\mathrm{q}}}{\partial {\epsilon }_{\mathrm{v}}^{\mathrm{p}}}-\left(\frac{\partial f}{\partial {\overline{\alpha }}_{\mathrm{p}}}{S}_p+\frac{\partial f}{\partial {\overline{\alpha }}_{\mathrm{q}}}{S}_q-\frac{\partial f}{\partial R}\right)\frac{\partial R}{\partial {\epsilon }_{\mathrm{v}}^{\mathrm{p}}}+\frac{\partial f}{\partial p_x}·\frac{\partial p_x}{\partial {\epsilon }_{\mathrm{v}}^{\mathrm{p}}}\right]-\\ \frac{\partial f}{\partial q}\left[\frac{\partial f}{\partial {\overline{\alpha }}_{\mathrm{q}}}(\mathrm{1}-R)\frac{\partial S_{\mathrm{q}}}{\partial {\epsilon }_{\mathrm{q}}^{\mathrm{p}}}-\left(\frac{\partial f}{\partial {\overline{\alpha }}_{\mathrm{p}}}{S}_p+\frac{\partial f}{\partial {\overline{\alpha }}_{\mathrm{q}}}{S}_q-\frac{\partial f}{\partial R}\right)\frac{\partial R}{\partial {\epsilon }_{\mathrm{q}}^{\mathrm{p}}}\right]\end{array}（42）

2.6 弹塑性本构关系

由弹塑性应变增量的定义可知：

\mathrm{d}\mathit{\epsilon }=\mathrm{d}{\mathit{\epsilon }}^{\mathrm{e}}+\mathrm{d}{\mathit{\epsilon }}^{\mathrm{p}}（43）

式中，\mathrm{d}\mathit{\epsilon }为总应变增量，\mathrm{d}{\mathit{\epsilon }}^{\mathrm{e}}为弹性应变增量，\mathrm{d}{\mathit{\epsilon }}^{\mathrm{p}}为塑性应变增量。

对式（43）两边同乘以弹性刚度矩阵{\mathit{D}}^{\mathrm{e}}有

{\mathit{D}}^{\mathrm{e}}\mathrm{d}\mathit{\epsilon }={\mathit{D}}^{\mathrm{e}}\mathrm{d}{\mathit{\epsilon }}^{\mathrm{e}}+{\mathit{D}}^{\mathrm{e}}\mathrm{d}{\mathit{\epsilon }}^{\mathrm{p}}（44）

由胡克定律可得：

\mathrm{d}{\mathit{\sigma }}^{\text{'}}={\mathit{D}}^{\mathrm{e}}\mathrm{d}{\mathit{\epsilon }}^{\mathrm{e}}={\mathit{D}}^{\mathrm{e}}(\mathrm{d}\mathit{\epsilon }-\mathrm{d}{\mathit{\epsilon }}^{\mathrm{p}})（45）

又因为

\mathrm{d}{\mathit{\epsilon }}^{\mathrm{p}}=〈L〉\frac{\partial f}{\partial {\mathit{\sigma }}^{\text{'}}}，L=\frac{\mathrm{1}}{K_{\mathrm{p}}}·\frac{\partial f}{\partial {\mathit{\sigma }}^{\text{'}}}\mathrm{d}{\mathit{\sigma }}^{\text{'}}（46）

将式（46）代入到式（45）可得加载指数〈L〉的表达式为：

〈L〉=\frac{{\left\{\frac{\partial f}{\partial {\mathit{\sigma }}^{\text{'}}}\right\}}^{\mathrm{T}}{\mathit{D}}^{\mathrm{e}}\mathrm{d}\mathit{\epsilon }}{K_{\mathrm{p}}+{\left\{\frac{\partial f}{\partial {\mathit{\sigma }}^{\text{'}}}\right\}}^{\mathrm{T}}{\mathit{D}}^{\mathrm{e}}\frac{\partial f}{\partial {\mathit{\sigma }}^{\text{'}}}}（47）

将式（47）和（46）代入到式（45）得：

\mathrm{d}{\mathit{\sigma }}^{\text{'}}=\left({\mathit{D}}^{\mathrm{e}}-\frac{{\mathit{D}}^{\mathrm{e}}\frac{\partial f}{\partial {\mathit{\sigma }}^{\text{'}}}{\left\{\frac{\partial f}{\partial {\mathit{\sigma }}^{\text{'}}}\right\}}^{\mathrm{T}}{\mathit{D}}^{\mathrm{e}}}{K_{\mathrm{p}}+{\left\{\frac{\partial f}{\partial {\mathit{\sigma }}^{\text{'}}}\right\}}^{\mathrm{T}}{\mathit{D}}^{\mathrm{e}}\frac{\partial f}{\partial {\mathit{\sigma }}^{\text{'}}}}\right)\mathrm{d}\mathit{\epsilon }={\mathit{D}}^{\mathrm{e}\mathrm{p}}\mathrm{d}\mathit{\epsilon }（48）

式中，{\mathit{D}}^{\mathrm{e}\mathrm{p}}为弹塑性刚度矩阵。

3 单调加载下模型三轴固结排水试验验证

3.1 三剪次加载面模型计算参数的确定

采用江西南昌地区非饱和重塑红黏土为试验土样，利用室内土工试验得到土样的相对重度、最大干密度、最优含水率、液限和塑限、泊松比和初始孔隙比等土工参数。利用GDS非饱和常规三轴试验仪进行控制基质吸力的压缩回弹试验，得到不同基质吸力下对应的非饱和土压缩指数和回弹指数，结果见表1。

根据不同基质吸力所对应的\lambda (s)和\kappa (s)值，结合式（21）做非线性拟合，得到回归参数\lambda (\mathrm{0})=0.057 39和{\lambda }_{\mathrm{s}}=0.022 41，\kappa (\mathrm{0})=0.009 88和{\kappa }_{\mathrm{s}}=-3.10×10^-5。

使用GDS非饱和常规三轴试验仪进行超固结非饱和土静三轴试验，采用SDT‒10温控动态土工三轴仪进行超固结饱和土循环三轴固结排水试验，得到在静、动荷载下超固结土的应力‒应变关系曲线。利用GDS非饱和三轴仪进行固结排水试验，由s为0时的摩尔圆和强度包线，可得饱和土的有效黏聚力和有效内摩擦角。上述试验所得模型计算参数见表2。

Bishop应力变量法有效应力参数采用文献[25]所提出的方法来确定，计算公式为：

\chi ={\left[\frac{\mathrm{1}}{\mathrm{1}+{\left(\tilde{\alpha }s\right)}^n}\right]}^m（49）

式中，拟合参数值\tilde{\alpha }=1.368，n=0.08，m=3.76。

因与参考文献采用的土样相同，故可直接将上述拟合参数代入到式（49）中求得不同基质吸力值所对应的Bishop有效应力参数值。

3.2 三剪次加载面模型的三轴固结排水试验验证

在非饱和超固结土本构模型的基础上，采用Fortran语言编制该模型的计算程序，将超固结非饱和土常规三轴固结排水试验结果和本构模型理论计算结果进行对比分析，以验证新建立的本构模型的正确性。在单调加载条件下，试样中的基质吸力分别为0、50、100、150 kPa，先期固结压力为800 kPa，超固结比R_{\mathrm{O}\mathrm{C}}分别为16、8、4。试样以较高的先期固结应力固结完成后，保持基质吸力不变，根据试验所需要的超固结程度将围压卸载到所设定的值并保持12 h，从而使土样达到超固结状态。限于篇幅且考虑到结果均类似，故只列出s=\mathrm{0}、R_{\mathrm{O}\mathrm{C}}=\mathrm{16}，s=\mathrm{150} kPa、R_{\mathrm{O}\mathrm{C}}=\mathrm{8}和s=\mathrm{150} kPa、R_{\mathrm{O}\mathrm{C}}=\mathrm{4}这3种情况下轴向应变{\epsilon }_{\mathrm{1}}与偏应力q、体应变{\epsilon }_{\mathrm{v}}的关系曲线，如图1～3所示。由图1～3可知，本文本构理论模型计算结果和室内三轴固结排水试验结果的变化规律一致。当超固结比较大时，土体试样在剪切过程中的偏应力‒轴向应变曲线出现明显的峰值强度且在峰值强度后出现软化现象。本文本构模型可以较好地反映非饱和土在超固结状态下的应变软化特性，验证了理论模型的可靠性。

Bishop应力变量法所得结果与试验结果相比偏小5%～8%；Fredlund应力变量法在当轴向应变在10%以内时与三轴试验结果基本吻合，轴向应变超过10%以后模拟结果略高于试验结果。采用Fredlund应力变量法所得结果较Bishop应力变量法偏大，同时Fredlund应力变量也更为接近实际试验结果。原因是Fredlund应力变量法采用净应力和吸力两个应力变量来描述非饱和土的应力特征量，故可以更加准确地反映土体的力学特性。文献[19]中也有相似的结论表明Fredlund应力变量法较Bishop应力变量法更为接近实际试验结果。

4 循环荷载作用下模型的三轴固结排水试验验证

4.1 超固结饱和江西红黏土

将所得的Frdlund应力变量法理论计算结果与饱和土动三轴固结排水试验数据进行对比分析。u、C、w为次加载面模型的3个特有参数，其中：u可以控制次加载面向正常屈服面靠近的速率参数，可由中等变形速率的应力‒应变曲线的斜率初步确定；参数C通过控制曼辛系数的大小来影响滞洄曲线的宽度；w为极限相似中心面和正常屈服面的比值，其值不超过1^[6]。选取初始压实度为90%、前期固结压力为600 kPa、净围压为400 kPa、R_{\mathrm{O}\mathrm{C}}=1.5、振幅q_d为50 kPa的循环动三轴固结排水试验数据作为参数标定对象，采用文献[13]所用的方法对模型参数进行初步确定，依据加载—卸载—再加载应力路径的三轴试验结果对其进行不断地调整，直至数值计算结果能较好地拟合超固结非饱和黏性土的应力‒应变曲线为止。最终确定次加载面动荷载参数为u=220，C=370，w=0.9。

试样轴向应变与偏应力、循环次数N间关系对比曲线如图4所示。由图4可知，本文所建立的本构理论模型计算结果与超固结饱和土循环动三轴固结排水试验结果的变化规律一致，取得了较好的模拟预测结果，可以较好地反映出在循环荷载作用下超固结饱和土随着循环加载次数不断增加，滞洄圈面积变小的棘轮效应（塑性应变的累积特性）和曼辛效应（滞洄特性）等循环塑性特征。

4.2 超固结非饱和土

选取工况基质吸力为300 kPa、净围压为400 kPa、R_{\mathrm{O}\mathrm{C}}=1.38、振幅q_{\mathrm{d}}为250 kPa，模型所需参数为e₀=0.705、v=0、p_{\mathrm{y}}^{\text{'}}(\mathrm{0})=50 kPa、p_{\mathrm{n}}^{\text{'}}=43 kPa、\lambda (\mathrm{0})=\mathrm{0.078}、\kappa (\mathrm{0})=\mathrm{0.01}、\lambda (s)=\mathrm{0.073}、\kappa (s)=\mathrm{0.015}，基质吸力为300 kPa时两种模型参数u、C和w分别均取50、270和0.7。对Fredlund应力变量法和Bishop应力变量法所预测的结果与文献[10]中的双面模型预测结果进行对比验证，结果如图5所示。

由图5可知，本文所建立的本构模型计算值与李潇旋等^[10]所建立的超固结非饱和土的弹塑性双面模型模拟结果吻合较好，可以较好地体现明显的滞洄圈及滞洄圈随着循环次数的增加而不断减小、所产生的塑性变形也逐渐减小等循环荷载作用下的超固结非饱和土力学特性。在相同工况条件下，Fredlund应力变量法达到稳定时的塑性累积应变相对Bishop应力变量法的结果更大，也更接近实际试验结果。

5 三剪次加载面模型的真三轴数值模拟

5.1 超固结非饱和土静真三轴数值模拟

基于超固结非饱和土三剪次加载面本构模型，利用Fortran语言编写可体现中间主应力影响系数b值影响的超固结非饱和土固结排水条件下真三轴静、动循环的数值计算程序。静真三轴数值模拟的对象是江西红黏土，以在Fredlund应力变量法下采用坐标平移法的次加载面模型为例，所选取的工况为：{\sigma }_{\mathrm{2}}^{\text{'}}=200 kPa、{\sigma }_{\mathrm{3}}^{\text{'}}=100 kPa，基质吸力s=50 kPa，中间主应力系数b=0、0.25、0.5、0.75、1.00，超固结比为8，先期固结压力为800 kPa。轴向应变{\epsilon }_{\mathrm{1}}与主应力差{\sigma }_{\mathrm{1}}^{\text{'}}-{\sigma }_{\mathrm{2}}^{\text{'}}、体应变{\epsilon }_{\mathrm{v}}间的数值模拟关系曲线如图6所示。由图6可知，中间主应力影响系数会对超固结非饱和土的强度峰值及体应变产生影响。在轴向应变较小时，b值对曲线影响较小，不同b值下的模拟曲线较为接近；随着轴向应变的增大，偏应力及体应变曲线差值逐渐增大。当保持基质吸力、围压和超固结比等条件不变时，随着中间主应力影响系数的增大，曲线的峰值强度随之增大。

5.2 超固结非饱和土循环动真三轴模拟分析

超固结非饱和土循环动真三轴试验存在仪器成本高、耗费时间长、操作不当会导致试验结果误差较大等问题。数值试验的方法则具有成本低、耗时短等优势，还能够确保生成试样的同一性、加载方式和模型边界的多样性，弥补因缺乏条件而无法进行室内试验的缺陷。本文采用PFC3D真三轴颗粒流模拟代替室内真三轴试验对所建本构模型进行动真三轴试验验证^[26‒27]。在PFC3D颗粒流建模中，采用Hill接触模型为颗粒间直接赋予基质吸力，采用柔性边界面模拟室内土工试验橡皮模对土样的包裹。

5.2.1 细观参数标定

由文献[28]可知，阻尼系数对试验影响很小，故可取常数0.7。相对于砂土颗粒，黏性土颗粒具有较大黏结强度，选取黏性土体颗粒间黏结强度为\mathrm{2.0}\times {\mathrm{10}}^{\mathrm{6}} N^[29]。数值模拟试验的超固结土土样压实度为90%，前期固结压力为800 kPa，基质吸力s为50 kPa，当前围压p^{\text{'}}为50 kPa，超固结比为16。在常规静三轴固结排水条件下，将模拟结果与土工试验结果进行对比并进行细观参数标定。

由以往经验可知，黏性土体的杨氏模量基本上在\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{6}}～\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{8}} Pa范围内，本文选取杨氏模量为\mathrm{5.0}\times {\mathrm{10}}^{\mathrm{7}}、\mathrm{3.0}\times {\mathrm{10}}^{\mathrm{7}}、\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{7}}、\mathrm{5.0}\times {\mathrm{10}}^{\mathrm{6}} Pa进行敏感性分析。在PFC3D中，黏性土泊松比取值一般为0～0.5，本文选取泊松比为0.30、0.35、0.40进行参数标定。摩擦系数值在0到1之间，根据经验选取摩擦系数为0.2、0.3、0.4进行标定。柔性膜颗粒之间的接触并不需要重点研究，只需将膜颗粒与土体颗粒之间的黏结强度设为一个较大值（\mathrm{2.0}\times {\mathrm{10}}^{\mathrm{6}} N），保证膜不会破裂即可。刚性墙体接触刚度选取\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{7}} N/m^[30]，膜颗粒密度选取\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{3}}\mathrm{k}\mathrm{g}\cdot {\mathrm{m}}^{-\mathrm{3}}^[30]，膜颗粒和刚性墙体的摩擦系数选取0.7^[30]，同时将膜颗粒的接触刚度设置为\mathrm{5.0}\times {\mathrm{10}}^{\mathrm{3}}、\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{4}}和\mathrm{2.0}\times {\mathrm{10}}^{\mathrm{4}} N/m进行敏感性分析。图7～10为不同杨氏模量、泊松比、摩擦系数和不同膜颗粒接触刚度条件下的PFC3D数值模拟结果与试验结果对比。

由图7～10可知，当颗粒杨氏模量为\mathrm{3.0}\times {\mathrm{10}}^{\mathrm{7}} Pa、泊松比为0.35、摩擦系数为0.3、膜颗粒的接触刚度为\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{4}} N/m时，模拟结果与试验结果最为匹配，由此标定出一组适用于超固结非饱和黏性土的PFC3D数值模拟细观参数，见表3和4。

5.2.2 建模及结果对比分析

PFC3D真三轴模型为80.0 mm×39.1 mm×39.1 mm的长方体，如图11所示。

荷载施加方法：1）首先对六面体单元的3个方向同时施加最小净主应力{\sigma }_{\mathrm{3}}^{\text{'}}，固结排水；2）固结完成后，一个方向应力{\sigma }_{\mathrm{3}}^{\text{'}}保持不变，其余两个方向应力增至中间净主应力{\sigma }_{\mathrm{2}}^{\text{'}}，固结排水；3）一个方向应力{\sigma }_{\mathrm{2}}^{\text{'}}继续保持不变，另一个方向应力在 {\sigma }_{\mathrm{2}}^{\text{'}}基础上施加循环荷载，排水。当轴向应变达到15%或者试样破坏时认定试验结束。

施加的循环荷载振幅为50 kPa，振动频率为0.5 Hz，先期固结压力为800 kPa。由于土质、固结排水条件、先期固结压力相同，循环动真三轴数值试验采用前述标定过的参数。将本构模型动真三轴计算结果与PFC3D模拟结果进行对比，如图12所示。由图12可知，本构理论计算结果与PFC3D数值模拟结果所表现出的趋势基本一致。当其他条件一定时，Fredlund应力变量法相对于Bishop应力变量法所得到的累计塑性应变较大，也更为接近PFC3D模拟结果。总体而言，所建立的本构模型可以体现出固结排水条件下超固结非饱和土体在循环动真三轴条件下的曼辛效应、棘轮效应等力学特性，可以反映中间主应力、最小主应力、振动幅值和基质吸力等对土体在复杂条件下力学性能的影响，验证了所建立模型的合理性。

6 结 论

本文在非饱和土修正剑桥模型基础上，引入统一硬化参数、次加载面理论和非饱和土三剪破坏应力比，建立了超固结非饱和黏性土三剪次加载面模型。将江西南昌地区重塑红黏土作为试验土样，采用GDS非饱和常规三轴试验仪进行超固结非饱和土固结排水静三轴试验，采用SDT‒10温控动态土工三轴仪进行超固结饱和土固结排水循环三轴试验，并将本构模型计算结果与试验结果进行了对比验证，结论如下：

1）该模型能够反映超固结土的硬化、软化、剪胀、剪缩等特性，还可以体现超固结土在循环荷载作用下的棘轮效应（塑性应变的累积特性）和曼辛效应（滞洄特性）。

2）本构模型计算结果与试验结果变化规律一致，验证了所建模型的正确性。在超固结土静三轴固结排水试验中，当轴向应变在10%以内时，Fredlund应力变量法与三轴试验结果基本吻合；当应变大于10%时，Fredlund应力变量法模拟结果略高于试验结果。Bi-shop应力变量法所得结果与试验结果相比偏小5%～8%。在循环动三轴试验中，在相同工况条件下，Fredlund应力变量法达到稳定时的塑性累积应变相对于Bishop应力变量法的计算结果更大，也更为接近实际试验结果。

3）真三轴数值模拟结果表明：中间主应力影响系数越大，剪应力峰值和体应变越大，土体抵抗破坏和变形的能力增强；常规静态三轴体现在应力峰值的增大，动循环三轴体现在土体达到稳定时的累计塑性应变减小。

参考文献

[1] Dafalias Y F, Popov E P.Plastic internal variables formalism of cyclic plasticity[J].Journal of Applied Mechanics,1976,43(4):645‒651. doi:10.1115/1.3423948

[2] Dafalias Y F.Plastic spin:Necessity or redundancy?[J].International Journal of Plasticity,1998,14(9):909‒931. doi:10.1016/s0749-6419(98)00036-9

[3] Dafalias Y F, Manzari M T, Akaishi M.A simple anisotropic clay plasticity model[J].Mechanics Research Communications,2002,29(4):241‒245. doi:10.1016/s0093-6413(02)00252-5

[4] Dafalias Y F.Bounding surface plasticity.I:Mathematical foundation and hypoplasticity[J].Journal of Engineering Mechanics,1986,112(9):966‒987. doi:10.1061/(asce)0733-9399(1986)112:9(966)

[5] Hashiguchi K, Chen Z.Elastoplastic constitutive equation of soils with the subloading surface and the rotational hardening[J].International Journal for Numerical and Analytical Methods in Geomechanics,1998,22(3):197‒227. doi:10.1002/(sici)1096-9853(199803)22:3<197::aid-nag914>3.0.co;2-t

[6] Hashiguchi K P, Ozaki S, Okayasu T.Unconventional friction theory based on the subloading surface concept[J].International Journal of Solids and Structures,2005,42(5/6):1705‒1727. doi:10.1016/j.ijsolstr.2004.08.006

[7] Nakai T, Hinokio M.A simple elastoplastic model for normally and over consolidated soils with unified material parameters[J].Soils and Foundations,2004,44(2):53‒70. doi:10.3208/sandf.44.2_53

[8] Houlsby G T, Wroth C P, Wood D M.Predictions of the results of laboratory tests on a clay using a critical state model[C]//Proceedings of the International Symposium Constitutive Laws of Soils.Grenoble,1982:99‒121.

[9] Mita K A, Dasari G R, Lo K W.Performance of a three-dimensional hvorslev-modified cam clay model for overconsolidated clay[J].International Journal of Geomechanics,2004,4(4):296‒309. doi:10.1061/(asce)1532-3641(2004)4:4(296)

[10] Li Xiaoxuan, Li Tao, Li Jian.Elastoplastic two-surface model for overconsolidated unsaturated soils[J].Journal of Hydraulic Engineering,2020,51(10):1278‒1288.

[11] 李潇旋,李涛,李舰.超固结非饱和土的弹塑性双面模型[J].水利学报,2020,51(10):1278‒1288.

[12] Li Xiaoxuan.Study on elastic-plastic double-sided model of unsaturated soil under static and cyclic loads[D].Beijing:Beijing Jiaotong University,2020.

[13] 李潇旋.静态与循环荷载下非饱和土的弹塑性双面模型研究[D].北京:北京交通大学,2020.

[14] Yao Yangping, Hou Wei, Luo Ting.Unified hardening model for soils[J].Chinese Journal of Rock Mechanics and Engineering,2009,28(10):2135‒2151.

[15] 姚仰平,侯伟,罗汀.土的统一硬化模型[J].岩石力学与工程学报,2009,28(10):2135‒2151.

[16] Hashiguchi K.Foundations of Elastoplasticity:Subloading surface model[M].Cham:Springer International Publising,2017. doi:10.1007/978-3-319-48821-9_22

[17] Hashiguchi K.Exact formulation of subloading surface model:Unified constitutive law for irreversible mechanical phenomena in solids[J].Archives of Computational Methods in Engineering,2016,23(3):417‒447. doi:10.1007/s11831-015-9148-x

[18] Yu Maohong, Yang Songyan, Fan Shouchang,et al.Twin shear unified elasto-plastic constitutive model and its applications[J].Chinese Journal of Geotechnical Engineering,1997,19(6):2‒10.

[19] 俞茂宏,杨松岩,范寿昌,等.双剪统一弹塑性本构模型及其工程应用[J].岩土工程学报,1997,19(6):2‒10.

[20] Hu Xiaorong, Yu Maohong.Research on triple-shear yield criterion for materials[J].Engineering Mechanics,2006,23(4):6‒11.

[21] 胡小荣,俞茂宏.材料三剪屈服准则研究[J].工程力学,2006,23(4):6‒11.

[22] Dong Xiaolong, Hu Xiaorong, Chen Xiaoyu.Triple-shear u-nified structured constitutive model for undisturbed saturated clay in normal consolidation[J].Journal of Civil and Environmental Engineering,2018,40(3):86‒93.

[23] 董肖龙,胡小荣,陈晓宇.正常固结原状饱和黏性土的三剪统一结构性本构模型[J].土木建筑与环境工程,2018,40(3):86‒93.

[24] Hu Xiaorong, Wang Ritang, Dong Xiaolong.Approaches to the triple-shear elasto-plasticity bounding surface model for saturated sands(1)—Model theory[J].Chinese Journal of Applied Mechanics,2020,37(2):580‒588.

[25] 胡小荣,汪日堂,董肖龙.饱和砂土三剪弹塑性边界面模型研究(一)——模型理论[J].应用力学学报,2020,37(2):580‒588.

[26] Hu Xiaorong, Cai Xiaofeng, Li Chunbo,et al.Triple-shear elasto-plastic constitutive models of unsaturated clays in normal consolidation[J].Engineering Mechanics,2022,39(1):164‒174.

[27] 胡小荣,蔡晓锋,李春博,等.正常固结非饱和黏性土的三剪弹塑性本构模型[J].工程力学,2022,39(1):164‒174.

[28] Hu Xiaorong, Cai Xiaofeng, Chen Hao.Triple-shear bounding surface models for unsaturated clays in normal consolidation[J].Engineering Mechanics,2021,38(6):175‒190.

[29] 胡小荣,蔡晓锋,陈昊.正常固结非饱和黏性土的三剪边界面模型研究[J].工程力学,2021,38(6):175‒190.

[30] Sun Dean.Hydro-mechanical behaviours of unsaturated soils and their elastoplastic modelling[J].Rock and Soil Mechanics,2009,30(11):3217‒3231.

[31] 孙德安.非饱和土的水力和力学特性及其弹塑性描述[J].岩土力学,2009,30(11):3217‒3231.

[32] 罗汀,姚仰平,侯伟.土的本构关系[M].北京:人民交通出版社,2010.

[33] Kong Liang, Zheng Yingren, Yao Yangping.Subloadoing su-rface cyclic plastic model for soil based on generalized plasticity (Ⅱ):Constitutive equation and identification[J].Rock and Soil Mechanics,2003,24(3):349‒354.

[34] 孔亮,郑颖人,姚仰平.基于广义塑性力学的土体次加载面循环塑性模型(Ⅱ):本构方程与验证[J].岩土力学,2003,24(3):349‒354.

[35] Hashiguchi K.Subloading surface model in unconventi-onal plasticity[J].International Journal of Solids and St-ructures,1989,25(8):917‒945. doi:10.1016/0020-7683(89)90038-3

[36] Liu Xiaowen, Ye Yunxue.Experimental study of the soil-water characteristic curve of unsaturated laterite under different affecting factors[J].Hydrogeology & Engineering Geology,2015,42(2):97‒104.

[37] 刘小文,叶云雪.不同影响因素下非饱和红土土‒水特征曲线的试验研究[J].水文地质工程地质,2015,42(2):97‒104.

[38] Zhu W C, Tang C A.Numerical simulation of Brazilian di-sk rock failure under static and dynamic loading[J].International Journal of Rock Mechanics and Mining Sciences,2006,43(2):236‒252. doi:10.1016/j.ijrmms.2005.06.008

[39] Lee H, Moon T, Haimson B C.Borehole breakouts induced in arkosic sandstones and a discrete element analysis[J].Rock Mechanics and Rock Engineering,2016,49(4):1369‒1388. doi:10.1007/s00603-015-0812-0

[40] Li Yao, Fu Yingpeng, Liao Hongjian,et al.Discrete element model of unsaturated soil considering suction and influence analysis of meso-parameters[J].Chinese Journal of Geotechnical Engineering,2021,43(Supp1):246‒250.

[41] 李瑶,伏映鹏,廖红建,等.考虑吸力作用的非饱和土离散元模型及细观参数影响分析[J].岩土工程学报,2021,43():246‒250.

[42] Zha Chonglun.Study on constitutive model of overconsolidated unsaturated soil and particle flow simulation[D].Nanchang:Nanchang University,2021.

[43] 查崇伦.超固结非饱和土本构模型及颗粒流模拟研究[D].南昌:南昌大学,2021.

[44] Hu Xiaorong, Cai Xiaofeng, Liu Cao.Approaches to the triple-shear elasto-plasticity bounding surface model for saturated sands(3)—PFC3D numerical test verifications[J].Chinese Journal of Applied Mechanics,2022,39(2):324‒335.

[45] 胡小荣,蔡晓锋,刘操.饱和砂土的三剪弹塑性边界面模型研究(三)——PFC3D数值试验验证[J].应用力学学报,2022,39(2):324‒335.

基金资助

国家自然科学基金项目(10862004)

国家自然科学基金项目(11362016)