饱和软黏土的ESF弹塑性模型及其验证

ESF Elastoplastic Model for Saturated Soft Clay and Model Validation

摘要

为描述饱和软黏土的强度变形特性并解决传统弹塑性模型中屈服面角点处塑性应变增量方向不够明确的问题，建立了一种饱和软黏土的蛋形屈服面弹塑性本构模型.首先，基于蛋形函数（ESF），采用相关联流动法则，以塑性体应变作为内变量，建立了饱和软黏土的ESF弹塑性模型.模型通过引入封闭光滑的蛋形屈服面，克服了传统屈服面的角点缺陷.通过调整蛋形形状控制参数，将屈服面灵活转化为包括子弹形、椭圆形在内的多种形式，体现了模型的可退化性及通用性.然后，开展了杭州地区饱和软黏土的应力路径三轴不排水剪切试验，分析了软黏土的力学特征，并据此研究了模型参数在不同围压下的稳定性；最后，编制隐式积分算法有限元程序，利用杭州地区饱和软黏土和相关文献中饱和高岭土在不同围压下的等压固结不排水三轴试验验证模型的有效性.对比结果表明：模型计算值与实测值吻合良好，说明该模型能够合理预测不同类型饱和软黏土的不排水剪切特性.

Abstract

An egg-shaped yield surface elastoplastic constitutive model suitable for saturated soft clay was proposed to describe the strength and deformation characteristics of saturated soft clay and overcome the defect that the direction of the plastic strain increments at the corner of the yield surface is insufficiently defined in the traditional elastoplastic model. Firstly， the ESF elastoplastic model was formulated based on the egg-shaped function （ESF）， adopting the associated flow law and employing the plastic volumetric strain as the internal variable. By introducing a closed and smooth egg-shaped yield surface， the corner problem in the traditional yield surface form was avoided. Additionally， by adjusting the value of egg-shaped parameters， the yield surface was flexibly transformed into various forms， including bullet-shaped and elliptical forms， which reflected the dependability and universality of this model. Subsequently， the undrained stress path triaxial tests were conducted on Hangzhou saturated soft clay to investigate its mechanical characteristics. The acquired data were used to analyze the stability of the model parameters under different confining pressures. Furthermore， a comprehensive analysis of the mechanical characteristics of the soft clay was performed， and the method for determining the model parameters was clarified. Finally， the finite element program using an implicit integration algorithm was developed， and its effectiveness was demonstrated through triaxial tests conducted on Hangzhou saturated soft clay and saturated kaolin clay from relevant literature under undrained conditions. The results revealed a satisfactory agreement between the simulated and experimental data， indicating that the proposed model can reasonably describe the behavior of various types of saturated soft clay under undrained conditions.

Graphical abstract

关键词

饱和软黏土 / 蛋形屈服面 / 弹塑性本构模型 / 不排水三轴试验 / 有限元程序

近年来，随着我国经济的快速发展和城市化进程的不断推进，饱和软黏土地区工程建设项目日益增多.饱和软黏土具有高含水量、高压缩性和低承载力等不良工程特性，因此开展其力学特性和本构理论研究对于解决实际工程问题具有重要价值.国内外众多学者提出了各类本构模型以研究饱和软黏土的弹塑性力学变形行为.而本构模型中屈服函数的选取对于模型能否准确预测土体的力学性质有着紧密关系.

早期弹塑性本构模型的屈服函数在子午面上多表现为线性开口形式，如扩展Mises准则、Mohr-Couloumb准则、Druck-Prager准则等^［1-4］.这类模型假定屈服强度会随静水压力增加而无限增大^［5］，这显然与岩土体材料的实际特性不符.为此，Drucker等^［6］考虑土体材料体积屈服特点，提出了加工硬化弹塑性理论，建立了盖帽模型.随后，Roscoe等^［7-8］基于临界状态理论建立发展了剑桥模型及其修正形式，这在土的弹塑性本构模型的发展过程中具有重要意义.然而针对不同土体表现出的不同屈服特性，构建屈服面形式更加广泛灵活的本构模型来预测土体的力学行为往往具有更加实际的工程价值.因此，数年间多数学者或基于修正剑桥模型的椭圆形式屈服函数引入形状控制参数，或基于岩土体材料试验数据，发展新的屈服函数形式.黄文熙^［9］根据试验资料，确定并提出了含有两个形状控制参数的椭圆形式屈服面的清华模型.Lagioia等^［10］建立了能够考虑砂土剪胀行为的弹塑性本构模型，在剑桥模型的基础上引入形状控制参数，并将其推广至三维应力空间.Kouretzis等^［11］针对松砂的不排水特性提出了泪滴形屈服面. Kan等^［12］、Moghadam等^［13］相继将该泪滴形屈服函数发展到边界面模型中.万征等^［14］基于超固结土UH模型通过引入状态参量构建水滴形屈服面，并将屈服面在p-q坐标系中整体左移来更好地描述黏土的结构性.Chen等^［15］为准确模拟超固结黏土的剪切响应，采用两个模型参数控制屈服面形状，提出了新的屈服面形式.相比之下，这类以修正剑桥为代表的椭圆形及其发展形式的屈服面能更为准确地描述土体屈服特性.然而由于这类屈服函数在偏应力空间中与临界状态破坏线的交点处存在奇异性^［16］，其塑性应变增量方向不明确，不利于数值计算.

为解决角点奇异性问题，Yu^［17］和Wathugala等^［18］提出了屈服面角点的非平滑处理方式，但此类方法仍存在塑性应变不连续的问题，且涉及复杂数学表达式，难以在本构理论中得到广泛应用与发展.克服屈服面角点缺陷的另一种方法是采用封闭光滑形式的屈服面.徐日庆等^［19］针对岩土材料提出了蛋形屈服函数及蛋形边界面函数.该函数是一种处处光滑连续的封闭曲线，克服了组合式屈服面的角点问题，有利于数值计算.同时该函数在一定条件下可以灵活转化为包括子弹形、椭圆形在内的多种形式，可用于不同类型土体特性的预测.然而该模型参数的物理意义及确定方式并未得到准确研究.为解决这一问题，Xu等^［20］通过水泥土的剪切波速测试试验建立了蛋形形状控制参数与剪切波速的联系，提出了模型参数确定方法.蒋佳琪等^［21］和徐日庆等^［22］在此基础上，结合等效塑性功硬化原理构建了适用于超固结淤泥质软土的蛋形弹塑性本构模型，并开展了软土的排水三轴剪切试验，验证了模型的有效性. 徐日庆等^［23］通过三轴排水剪切试验验证了高围压条件下蛋形模型参数的稳定性.然而这些研究主要基于排水三轴试验开展，并未验证不排水条件下蛋形模型的有效性及模型参数稳定性，难以满足实际工程需要.此外，较多的模型参数在一定程度上也影响了模型的广泛应用.因此，有必要对蛋形模型进行进一步研究.

本文基于简化形式的蛋形函数（ESF），以塑性体应变作为等向硬化过程中的内变量，采用相关联流动准则，建立了能有效反映不排水条件下饱和软黏土力学特性的ESF弹塑性本构模型.随后开展了不同围压下饱和软黏土的应力路径三轴不排水剪切试验，分析了软黏土的力学特征，提出了模型参数确定方式并以三轴试验结果为依据验证了模型参数在不同围压下的稳定性.最后利用Fortran语言编制隐式积分算法有限元程序，获得了模型的数值解，通过与本文及相关文献中的试验数据对比，验证了该模型对于预测饱和软黏土力学变形特性的有效性.

1 ESF弹塑性本构模型

ESF弹塑性本构模型是基于弹塑性理论框架建立的，主要由蛋形屈服函数、蛋形塑性势函数、硬化准则、正交法则等组成.

1.1 屈服函数

蛋形屈服函数在p-q平面中的形式如图1所示，其表达式为：

F={\left(\frac{p-p_{\mathrm{0}}d}{a{p}_{\mathrm{0}}}\right)}^{\mathrm{2}}+{\left[\frac{a{p}_{\mathrm{0}}\left(\mathrm{1}-{\beta }^{\mathrm{2}}\right)}{a{p}_{\mathrm{0}}(\mathrm{1}+\beta )+\beta (p-p_{\mathrm{0}})}\right]}^{\mathrm{2}}{\left(\frac{q}{b{p}_{\mathrm{0}}}\right)}^{\mathrm{2}}-\mathrm{1}=\mathrm{0}（1）

式中：p、q分别为平均主应力和偏应力；p₀为土体先期固结应力；β为屈服面形状控制参数，与土体内摩擦角相关.

其中，ap₀/（1-β²）和bp₀/（1-β²） 分别表示蛋形屈服面的长短半轴；dp₀为函数几何中心的横坐标，其中d=1-a.

当β取不同值时，屈服面可以灵活演变为不同形态，如图2所示.当0<β<1时，屈服曲线形状呈蛋形；当-1<β<0时，曲线为子弹形；当β=0时，曲线还原为椭圆形，且当a=b时曲线退化为圆形；当β接近1时，曲线表现为角缘光滑的三角形.

1.2 流动法则

研究表明，相关联流动法则在饱和软黏土中具有较好的适应性^［24-25］，本文采用关联流动规则，即屈服面F和塑性势面Q采取相同形式：

F=Q（2）

塑性应变增量\mathrm{d}{\epsilon }_{ij}^{\mathrm{p}}表达式如下：

\mathrm{d}{\epsilon }_{ij}^{\mathrm{p}}=\mathrm{d}\lambda \frac{\partial F}{\partial {\sigma }_{ij}}（3）

式中：\mathrm{d}\lambda是非负的塑性标量因子.

1.3 硬化规律

本模型选取塑性体应变{\epsilon }_v^{\mathrm{p}}作为各向同性硬化过程中的内变量，将土体的先期固结压力p₀作为硬化参数.推导出模型硬化规律如下：

p_{\mathrm{0},n+\mathrm{1}}=p_{\mathrm{0},n}\mathrm{e}\mathrm{x}\mathrm{p}\left\{\frac{\mathrm{1}+e_{\mathrm{0}}}{\lambda -\kappa }\left[({\epsilon }_v^{\mathrm{p}}{)}_{n+\mathrm{1}}-({\epsilon }_v^{\mathrm{p}}{)}_n\right]\right\}（4）

式中：e₀代表初始孔隙比；\lambda和\kappa分别表示e-lnp坐标下等向固结试验压缩曲线和回弹曲线的斜率；n代表增量步数.

1.4 增量弹塑性关系

根据增量弹塑性理论，土体总应变增量\mathrm{d}{\epsilon }_{ij}^{}为弹性应变增量\mathrm{d}{\epsilon }_{ij}^{\mathrm{e}}和塑性应变增量\mathrm{d}{\epsilon }_{ij}^{\mathrm{p}}之和，表示为：

\mathrm{d}{\epsilon }_{ij}=\mathrm{d}{\epsilon }_{ij}^{\mathrm{e}}+\mathrm{d}{\epsilon }_{ij}^{\mathrm{p}}（5）

弹性应变增量满足广义Hook定律，通过下式计算确定：

\left\{\begin{array}{c}\mathrm{d}p\\ \mathrm{d}q\end{array}\right\}=\left[\begin{array}{cc}K& \mathrm{0}\\ \mathrm{0}& \mathrm{3}G\end{array}\right]\left\{\begin{array}{c}\mathrm{d}{\epsilon }_v^{\mathrm{e}}\\ \mathrm{d}{\epsilon }_q^{\mathrm{e}}\end{array}\right\}（6）

式中：弹性体积模量K通常由等向压缩试验确定；弹性剪切模量G可以通过弹性体变模量K和泊松比v得到，具体见下式：

K=\frac{\mathrm{1}+e_o}{\kappa }p（7）

G=\frac{\mathrm{3}(\mathrm{1}-\mathrm{2}v)}{\mathrm{2}(\mathrm{1}+v)}K（8）

塑性应变增量分为塑性体应变增量\mathrm{d}{\epsilon }_v^{\mathrm{p}}和塑性剪应变增量\mathrm{d}{\epsilon }_q^{\mathrm{p}}：

\mathrm{d}{\epsilon }_v^{\mathrm{p}}=\mathrm{d}\lambda \frac{\partial F}{\partial p}（9）

\mathrm{d}{\epsilon }_q^{\mathrm{p}}=\mathrm{d}\lambda \frac{\partial F}{\partial q}（10）

该蛋形屈服函数的一致性方程如下所示：

\frac{\partial F}{\partial p}\mathrm{d}p+\frac{\partial F}{\partial q}\mathrm{d}q+\frac{\partial F}{\partial {\epsilon }_v^p}\mathrm{d}{\epsilon }_v^{\mathrm{p}}=\mathrm{0}（11）

其中，

\frac{\partial F}{\partial p}=\frac{\mathrm{2}(p-d{p}_{\mathrm{0}})}{(a{p}_{\mathrm{0}}{)}^{\mathrm{2}}}-

\frac{\mathrm{2}\beta (\mathrm{1}-{\beta }^{\mathrm{2}}{)}^{\mathrm{2}}}{[a{p}_{\mathrm{0}}(\mathrm{1}+\beta )+\beta (p-p_{\mathrm{0}}{)]}^{\mathrm{3}}}\frac{a^{\mathrm{2}}}{b^{\mathrm{2}}}{q}^{\mathrm{2}}（12）

\frac{\partial F}{\partial q}=\frac{\mathrm{2}(\mathrm{1}-{\beta }^{\mathrm{2}}{)}^{\mathrm{2}}}{[a{p}_{\mathrm{0}}(\mathrm{1}+\beta )+\beta (p-p_{\mathrm{0}}{)]}^{\mathrm{2}}}\frac{a^{\mathrm{2}}}{b^{\mathrm{2}}}q（13）

\frac{\partial F}{\partial {\epsilon }_v^{\mathrm{p}}}=-\frac{\mathrm{1}+e_{\mathrm{0}}}{\lambda -\kappa }{p}_{\mathrm{0}}\frac{\mathrm{2}(\mathrm{1}-{\beta }^{\mathrm{2}}{)}^{\mathrm{2}}(a-\beta +a\beta )}{[a{p}_{\mathrm{0}}(\mathrm{1}+\beta )+\beta (p-p_{\mathrm{0}}{)]}^{\mathrm{3}}}\frac{a^{\mathrm{2}}}{b^{\mathrm{2}}}\cdot

q^{\mathrm{2}}-\frac{\mathrm{1}+e_{\mathrm{0}}}{\lambda -\kappa }{p}_{\mathrm{0}}\frac{\mathrm{2}p(p-p_{\mathrm{0}}+a{p}_{\mathrm{0}})}{a^{\mathrm{2}}{p_{\mathrm{0}}}^{\mathrm{3}}}（14）

结合一致性方程式（11）与塑性应变增量表达 式（3），可得该模型的塑性模量K_p以及应力应变增量关系：

K_{\mathrm{p}}=-\frac{\partial F}{\partial {\epsilon }_v^{\mathrm{p}}}\frac{\partial F}{\partial p}（15）

\left\{\begin{array}{c}\mathrm{d}{\epsilon }_v\\ \mathrm{d}{\epsilon }_q\end{array}\right\}=\left[\begin{array}{cc}\frac{\mathrm{1}}{K}+\frac{\mathrm{1}}{K_{\mathrm{p}}}\frac{\partial F}{\partial p}\frac{\partial F}{\partial p}& \frac{\mathrm{1}}{K_{\mathrm{p}}}\frac{\partial F}{\partial p}\frac{\partial F}{\partial q}\\ \frac{\mathrm{1}}{K_{\mathrm{p}}}\frac{\partial F}{\partial p}\frac{\partial F}{\partial q}& \frac{\mathrm{1}}{\mathrm{3}G}+\frac{\mathrm{1}}{K_{\mathrm{p}}}\frac{\partial F}{\partial q}\frac{\partial F}{\partial q}\end{array}\right]\left\{\begin{array}{c}\mathrm{d}p\\ \mathrm{d}q\end{array}\right\}（16）

根据有效应力原理，土体的孔隙水压力增量du为平均主应力增量与平均有效主应力增量之差，表示为：

\mathrm{d}u=\mathrm{d}{\sigma }_{{}_{\mathrm{3}}}^{\mathrm{t}}+\frac{\mathrm{1}}{\mathrm{3}}\mathrm{d}q-\mathrm{d}p（17）

式中：\mathrm{d}{\sigma }_{\mathrm{3}}^{\mathrm{t}}表示围压增量，t表示总应力形式.

不排水状态下的土体应变增量为零：

\mathrm{d}{\epsilon }_v=\mathrm{d}{\epsilon }_v^{\mathrm{p}}+\mathrm{d}{\epsilon }_v^{\mathrm{e}}=\mathrm{0}（18）

将式（18）与弹性应变增量表达式（6）代入式（17），得到：

\mathrm{d}u=\mathrm{d}{\sigma }_{{}_{\mathrm{3}}}^{\mathrm{t}}+\frac{\mathrm{1}}{\mathrm{3}}\mathrm{d}q-K\mathrm{d}{\epsilon }_v^{\mathrm{e}}（19）

将一致性方程表达式（11）和塑性模量K_p表达式（15）代入式（19），得到孔隙水压力增量表达式为：

\mathrm{d}u=\mathrm{d}{\sigma }_{{}_{\mathrm{3}}}^{\mathrm{t}}+\frac{\mathrm{1}}{\mathrm{3}}\mathrm{d}q+K\frac{\partial F}{\partial p}\left[\frac{\mathrm{1}}{K_{\mathrm{p}}}\left(\frac{\partial F}{\partial p}\mathrm{d}p+\frac{\partial F}{\partial q}\mathrm{d}q\right)\right]（20）

2 固结不排水剪应力路径试验

2.1 试验材料与试样制备

试验所用土样取自杭州某工程软黏土，其物理力学指标如表1所示.

考虑原状黏土样易扰动，难以获取，本试验采用重塑土样.根据《土工试验方法标准》（GB/T 50123—2019）将土样烘干、碾碎、过筛、制备土样；标准固结试样高度20 mm，直径61.8 mm；三轴试样高度 80 mm，直径39.1 mm.过程中控制试样的含水率、孔隙比与原状土一致.将制备好的土样放置真空饱和容器中进行抽气饱和.

2.2 试验方案

2.2.1 标准固结试验

为获取土体压缩回弹特性，得到模型中的\lambda和\kappa等材料参数，本研究利用一维固结仪开展标准固结试验，设置试验加载压力共6级，分别为25 kPa、 50 kPa、100 kPa、200 kPa、400 kPa、800 kPa，逐级加压后逐级卸压，每级荷载维持24 h.

2.2.2 固结不排水剪应力路径试验

在三种工况下分别开展杭州地区饱和软黏土的应力路径三轴试验：初始固结压力p₀分别为100 kPa、200 kPa和300 kPa.各工况下分别设计3组应力路径试验，包括减压三轴压缩RTC路径、等p三轴压缩CMS路径以及常规三轴压缩CTC路径，如图3所示.具体试验方案见表2，其中应力增量比\mathrm{\Delta }\eta =\mathrm{\Delta }q/\mathrm{\Delta }p.考虑剪切过程中不同剪切控制方式会对孔压变化规律造成影响，应变控制相比应力控制更易产生孔压滞后现象^［26］，为避免孔压滞后，试验中均采用应力控制剪切速率，为0.2 kPa/min.

2.3 试验结果

2.3.1 标准固结试验结果

图4为饱和软黏土正常固结一维压缩试验的e-lnp曲线.试验所得数据点连线呈两条斜率不同直线，分别为压缩曲线和回弹曲线，曲线斜率分别代表正常固结黏土的压缩指数C_c和回弹指数C_e.

模型参数\lambda和\kappa分别表示重塑土e-lnp坐标下等向压缩段和回弹段的曲线斜率，并不等于上述试验得到的压缩指数C_c和回弹指数C_e.而相关研究表明^［27］，参数\lambda和\kappa可以通过一维压缩试验确定的压缩指数C_c以及回弹指数C_e进行转化，其转化关系可以表示为：

\begin{array}{l}\lambda =\mathrm{0.434}{C}_{\mathrm{c}}\\ \kappa =\mathrm{0.434}{C}_{\mathrm{e}}\end{array}（21）

通过计算，得到本文杭州地区软黏土的模型参数\lambda和\kappa分别为0.116 49和0.012 98.

2.3.2 不排水剪应力路径试验结果

图5~图7分别为土样在不同围压、不同应力路径下的偏应力-轴向应变关系曲线、孔隙水压力-轴向应变关系曲线和有效应力路径曲线.

图5为土样剪切过程中偏应力-轴向应变关系曲线，试验结果表明：试样在等压固结后不同剪切路径下得到的应力-应变曲线形态相似，总体上呈现明显的应变硬化特征.在同一初始固结压力下，减p、等p、增p三条剪切路径得到的土体峰值强度依次递增；初始固结压力增加时，土样的偏应力峰值也随之增加，同时递增幅值也有所增大；初始固结压力越大，剪切初期偏应力增加越显著，体现在应力-应变曲线上即剪切初期的曲线斜率越大，这主要是由于初始固结压力不同，引起试样抗剪强度在剪切初期存在差异.

在三轴不排水剪切试验过程中，试样通过内部压力的调整即超静孔隙水压力的变化来调整有效围压，以保持体积不变.由图6可知，在相同应力路径下，孔隙水压力（绝对值）随初始固结压力的增加而增加；在增p和等p应力路径试验中，孔隙水压力均为正值，试验初期的孔压迅速增大，土样表现出明显剪缩特性，试样屈服后孔隙水压力保持稳定或有所降低，表现出一定的剪胀性；在减p应力路径试验中，孔隙水压力为负值，在试样屈服前孔隙水压力变化较小，在屈服后负值明显扩大，主要考虑压应力减小引起的弹性变形和剪切引起的塑性变形的综合影响.

有效应力路径反映了试样在剪切过程中应力状态的变化.本试验中试样的有效应力路径如图7所示.不同初始固结压力下试样的有效应力路径曲线形态基本相同，在增p和等p应力路径试验中，有效应力路径曲线均处于总应力路径左侧，在试样屈服后有效应力路径曲线向靠近总应力路径方向发生偏转；减p应力路径试验中，在试样屈服前有效应力路径与总应力路径曲线基本一致，在试样屈服后向远离总应力路径方向偏转.同一初始固结压力下，三种不同应力路径三轴试验的有效应力路径形态基本一致，验证了正常固结土p^{\text{'}}-q-e的唯一性关系，即土体在等压固结至同一密度（e相同）后，在不排水条件下进行三轴剪切试验，由于不排水剪切过程中土体体积不变，孔隙比保持恒定，因此在p^{\text{'}}-q平面上对应唯一的有效应力路径.

3 模型参数确定及稳定性分析

ESF弹塑性模型共有7个模型参数，包括反映土体基本性质的4个参数e_{\mathrm{0}}\mathrm{、}\lambda \mathrm{、}\kappa \mathrm{、}\nu和反映屈服面形状的3个参数a\mathrm{、}b\mathrm{、}\beta.具体参数确定方法如下.

e_{\mathrm{0}}\mathrm{、}\lambda \mathrm{、}\kappa \mathrm{、}v为反映土体基本物理特性的材料参数.其中，e_{\mathrm{0}}为土体的初始孔隙比，可由物理性质指标计算；\lambda \mathrm{、}\kappa可以通过等向压缩回弹试验得到；土体泊松比v可根据实际情况选取，一般为0.3~0.5^［28-29］.

a\mathrm{、}b\mathrm{、}\beta为屈服面形状参数，可以结合不排水三轴试验结果和应力路径法^［30］进行确定.参数具体确定方式如下：根据式（22）将不排水三轴试验中有效应力路径上的点通过一系列等应力比线映射得到屈服轨迹，如图8所示；再利用蛋形屈服函数进行拟合，得到屈服面参数.

P_B=P_C{\left(\frac{P_C}{P_A}\right)}^{\frac{\kappa }{\lambda -\kappa }}（22）

根据杭州地区软黏土不排水三轴试验结果，有效应力路径与孔隙比具有唯一性关系，即同一围压下不同应力路径三轴试验的有效应力路径几乎相同.因此本文采用应力路径法得到的屈服轨迹也具有唯一性.为验证蛋形屈服面中模型参数在不同围压下的稳定性，本文采用应力路径法计算得到不同围压下的屈服轨迹，分别利用蛋形屈服函数拟合，得到不同围压下的屈服面参数.拟合结果如图9所示，图中数据点均为根据应力路径法计算得到的屈服轨迹上的点.拟合得到杭州地区软黏土的蛋形模型参数见表3.

从表3中可以发现，不同围压条件下得到的ESF模型参数变化较小，验证了该模型参数的稳定性.将3种不同围压下所得到的模型参数取平均值作为本文杭州地区饱和软黏土的ESF模型参数，用于后续模型验证.

4 ESF模型验证

利用Fortran语言编制了隐式积分算法有限元程序，获得模型的数值解.通过与本文试验及相关文献中的试验数据对比，验证ESF弹塑性模型的适用性和有效性.

4.1 饱和软黏土不排水三轴试验模拟

利用ESF弹塑性模型对杭州地区饱和软黏土三轴不排水剪切试验进行了模拟预测，得到不同围压下的模拟值.模拟结果和试验结果对比如图10、 图11所示.所采用的模型参数见表4.

图10为试验和预测结果的应力-应变曲线对比.可以看出，不同围压下的模拟曲线与实测曲线的趋势大致相同，均表现出应变硬化特性，在轴向应变较小时，模拟值与试验值基本吻合，随着轴向应变的增加，偏应力也逐渐增大.

图11为有效应力路径预测和试验结果.可以看出，模型能较好地预测试样加载过程中的应力状态变化，当试样加载逐渐达到屈服状态，偏应力随平均有效应力的减小而逐渐增大，有效应力路径预测结果与试验结果基本吻合.达到屈服状态后，试样逐渐呈现剪胀特性，因此有效应力路径向右偏转，而由于模型的硬化参量中尚未考虑剪胀的因素，导致有效应力路径的模拟与实测数据存在一定差异.整体而言，不同围压下的模型预测结果与试验结果符合良好，模拟曲线和试验点的变化规律基本一致，可以认为该模型能有效预测饱和软黏土的不排水剪切特性.

4.2 kaolin（高岭）黏土的等压固结三轴试验

Stipho等^［31］针对饱和高岭土分别开展了围压值为207 kPa和414 kPa的等压固结三轴不排水试验，现用其试验结果来进一步验证ESF弹塑性本构模型的合理性.模型各参数取值如表4所示，模拟值与实测值的对比如图12~图14所示.

图12~图14为利用本文模型对Stipho等^［31］开展的饱和高岭土等压固结常规三轴试验的模拟值及试验值对比.可以看出在不同围压下模型模拟结果与试验结果基本吻合，随着轴向应变增加，偏应力和孔隙水压力逐渐增加并趋向稳定，最终稳定时模拟值与试验值基本一致.同时模型能较好预测土体在剪切过程中应力状态的变化，其有效应力路径与试验值基本吻合.表明该ESF弹塑性模型程序具有较好的计算精度和数值稳定性，可以有效预测不排水条件下饱和高岭土的力学特性.

5 结 论

为准确反映饱和软黏土的力学特性，本文简化了蛋形屈服函数的形式，建立了ESF弹塑性本构模型，并以三轴试验为依据验证了模型参数在不同围压下的稳定性，采用Fortran语言编制隐式积分算法有限元子程序，获得了模型的数值解，通过与杭州地区饱和软黏土及饱和高岭土的不排水三轴试验结果进行对比，得到以下结论：

1）蛋形屈服函数是一种处处光滑连续的封闭曲线，克服了组合式屈服面的角点奇异性缺陷，有利于数值计算.同时该函数可以灵活转化为包括子弹形、椭圆形在内的多种形式，具有较强的普适性，能反映不同类型软黏土的强度和变形特性.

2）基于蛋形屈服函数，建立了能有效反映不排水饱和软黏土力学特性的ESF弹塑性本构模型.阐述了模型参数的确定方法，结合三轴试验结果和应力路径法得到了不同围压下的模型参数，结果表明模型参数在不同围压下具有较好的稳定性.

3）开展了饱和软黏土在不同围压下的应力路径三轴不排水剪切试验.在不同初始围压下，杭州地区软黏土均呈现应变硬化特征，初始围压越大，剪切初期偏应力增加越明显；在不同应力路径下，软黏土的有效应力路径形态基本一致，验证了正常固结土 p'-q-e的唯一性关系及模型参数确定方式的合理性.

4）利用软黏土和饱和高岭土的相关试验数据对本研究提出的模型进行了有效性验证，结果表明本文所建立的ESF弹塑性本构模型能够合理预测不同类型饱和软黏土的不排水剪切特性.

参考文献

[1] DRUCKER D C， PRAGER W ．Soil mechanics and plastic analysis or limit design［J］．Quarterly of Applied Mathematics，1952，10（2）：157-165．

[2] DRUCKER D C ．Limit analysis of two and three dimensional soil mechanics problems［J］．Journal of the Mechanics and Physics of Solids，1953，1（4）：217-226．

[3] 罗汀，姚仰平，侯伟 ．土的本构关系［M］．北京：人民交通出版社，2010．

[4] LUO T， YAO Y P， HOU W ．Soil constitutive models［M］．Beijing：China Communications Press，2010．（in Chinese）

[5] 程星磊，王建华，李书兆 ．软黏土不排水循环应力应变响应的弹塑性模拟［J］．岩土工程学报，2014，36（5）：933-941．

[6] CHENG X L， WANG J H， LI S Z ．Elastoplastic simulation for undrained cyclic stress-strain responses of soft clays［J］．Chinese Journal of Geotechnical Engineering，2014，36（5）：933-941．（in Chinese）

[7] MILLER T W， CHEATHAM J B ．A new yield condition and hardening rule for rocks［J］．International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts，1972，9（4）：453-474．

[8] DRUCKER D C， GIBSON R E， HENKEL D J ．Soil mechanics and work-hardening theories of plasticity［J］．Transactions of the American Society of Civil Engineers，1957，122（1）：338-346．

[9] ROSCOE K H， SCHOFIELD A N， THURAIRAJAH A ．Yielding of clays in states wetter than critical［J］．Géotechnique，1963， 13（3）： 211-240．

[10] ROSCOE K H， BURLAND J B. On the generalized stressstrain behavior of wet clay［J］. Engineering plasticity， 1968， 535-609.

[11] 黄文熙 ．土的弹塑性应力-应变模型理论［J］. 岩土力学， 1979（1）：1-20．

[12] HUANG W X. Theory of elastoplastic stress strain model for soil［J］．Rock and Soil Mechanics，1979（1）：1-20．（in Chinese）

[13] LAGIOIA R， PUZRIN A M， POTTS D M ．A new versatile expression for yield and plastic potential surfaces［J］．Computers and Geotechnics，1996，19（3）：171-191．

[14] KOURETZIS G P， SHENG D C， WANG D ．Numerical simulation of cone penetration testing using a new critical state constitutive model for sand［J］．Computers and Geotechnics，2014，56：50-60．

[15] KAN M E， TAIEBAT H A， KHALILI N ．Simplified mapping rule for bounding surface simulation of complex loading paths in granular materials［J］．International Journal of Geomechanics，2014，14（2）：239-253．

[16] MOGHADAM S I， TAHERI E， AHMADI M，et al ．Unified bounding surface model for monotonic and cyclic behaviour of clay and sand［J］．Acta Geotechnica，2022，17（10）：4359-4375．

[17] 万征，秋仁东，宋琛琛 ．基于统一硬化参量的原状饱和黏土的结构性本构模型［J］．岩石力学与工程学报，2019，38（9）：1905-1918．

[18] WAN Z， QIU R D， SONG C C ．A structural constitutive model of undisturbed saturated clays based on the uniform hardening parameter［J］．Chinese Journal of Rock Mechanics and Engineering，2019，38（9）：1905-1918．（in Chinese）

[19] CHEN Y N， YANG Z X ．A family of improved yield surfaces and their application in modeling of isotropically over-consolidated clays［J］．Computers and Geotechnics，2017，90：133-143．

[20] JU L Y， YAN Z H， WU M M，et al ．Implicit integration of an anisotropic egg-shaped elastoplastic model forsaturatedsoft clay［J］．Engineering Computations，2023，40（3）：517-537．

[21] YU H S ．A closed-form solution of stiffness matrix for Tresca and Mohr-Coulomb plasticity models［J］．Computers & Structures，1994，53（3）：755-757．

[22] WATHUGALA G W， DESAI C S ．Constitutive model for cyclic behavior of clays.I：theory［J］. Journal of Geotechnical Engineering，1993，119（4）：714-729．

[23] 徐日庆，杨林德，龚晓南 ．土的边界面应力应变本构关系［J］．同济大学学报（自然科学版），1997，25（1）： 29-33．

[24] XU R Q， YANG L D， GONG X N ．Stress-strain constitutive relationship of boundary surface of soils［J］．Journal of Tongji University （Natural Science），1997，25（1）：29-33．（in Chinese）

[25] XU R， WANG X， RONG X N. Macro-and meso-scale parametric analysis of egg-shaped yield surface of structural soil［J］．International Journal of Offshore and Polar Engineering，2014，24：301-308．

[26] 蒋佳琪， 徐日庆， 裘志坚， 等 ．超固结土的蛋形弹塑性本构模型［J］．浙江大学学报（工学版），2021，55（8）：1444-1452．

[27] JIANG J Q， XU R Q， QIU Z J，et al ．Egg-shaped elasto-plastic constitutive modeling for over-consolidated clay［J］. Journal of Zhejiang University （Engineering Science），2021，55（8）：1444-1452．（in Chinese）

[28] 徐日庆，蒋佳琪， 冯苏阳，等 ．一种旋转塑性势面模型及非关联塑性流动法则［J］．岩土力学，2020，41（5）：1474-1482．

[29] XU R Q， JIANG J Q， FENG S Y，et al ．A rotational plastic potential model and non-associated plastic flow rule［J］．Rock and Soil Mechanics，2020，41（5）：1474-1482．（in Chinese）

[30] 徐日庆，朱黄鼎，孙华俊，等 ．基于高围压三轴试验的蛋形弹塑性模型参数分析［J］．地基处理，2023，5（1）：9-18．

[31] XU R Q， ZHU H D， SUN H J，et al ．Parameter analysis of the egg-shaped elastoplastic model based on triaxial tests of high confining pressure［J］．Journal of Ground Improvement， 2023， 5（1）： 9-18．（in Chinese）

[32] 沈恺伦，王立忠 ．天然软黏土屈服面及流动法则试验研究［J］．土木工程学报，2009，42（4）：119-127．

[33] SHEN K L， WANG L Z ．Experimental study on the yield surface and flow rule of natural clays［J］．China Civil Engineering Journal，2009，42（4）：119-127．（in Chinese）

[34] 徐日庆，鞠露莹，俞建霖，等 ．饱和软黏土的蛋形单面边界面模型［J］．岩土工程学报，2020，42（12）：2170-2179．

[35] XU R Q， JU L Y， YU J L，et al ．Egg-shaped bounding surface model for saturated soft clay［J］．Chinese Journal of Geotechnical Engineering，2020，42（12）：2170-2179．（in Chinese）

[36] 曾玲玲，陈晓平 ．软土在不同应力路径下的力学特性分析［J］．岩土力学，2009，30（5）：1264-1270．

[37] ZENG L L， CHEN X P ．Analysis of mechanical characteristics of soft soil under different stress paths［J］．Rock and Soil Mechanics，2009，30（5）：1264-1270．（in Chinese）

[38] YIN J H， GRAHAM J ．Elastic viscoplastic modelling of the time-dependent stress-strain behaviour of soils［J］．Canadian Geotechnical Journal，1999，36（4）：736-745．

[39] 陈树峰，孔令伟，黎澄生 ．低幅应变条件下粉质黏土泊松比的非线性特征［J］．岩土力学，2018，39（2）：580-588．

[40] CHEN S F， KONG L W， LI C S ．Nonlinear characteristics of Poisson’s ratio of silty clay under low amplitude strain［J］．Rock and Soil Mechanics，2018，39（2）：580-588．（in Chinese）

[41] 汤大明， 曾纪全，胡应德，等. 关于泊松比的试验和取值讨 论［J］．岩石力学与工程学报，2001，20（）：1772-1775．

[42] TANG D M， ZENG J Q， HU Y D，et al ．Discussion on the experiment and value of Poisson’s ratio［J］．Chinese Journal of Rock Mechanics and Engineering，2001，20（Sup.1）：1772-1775．（in Chinese）

[43] 施维成，朱俊高 ．建立土体屈服函数的方法研究［J］．岩土工程学报，2010，32（4）：624-629．

[44] SHI W C， ZHU J G ．Methods for establishing yield function of soils［J］．Chinese Journal of Geotechnical Engineering，2010， 32（4）： 624-629．（in Chinese）

[45] STIPHO S A S A ．Experimental and theoretical investigation of the behaviour of anisotropically consolidated Kaolin［D］．Cardiff，Wales，UK：Cardiff University，1978

基金资助

浙江省重点研发计划资助项目（2019C03103）