不同应力路径下孔隙比对于饱和砂土塑性应变影响的理论研究

段秋宇，戚承志，王皓楠，封焱杰，罗伊，马啸宇

（1.北京建筑大学北京未来城市设计高精尖创新中心，北京100044；2.北京建筑大学北京城市交通基础设施建设国际合作基地，北京100044；3.北京建筑大学土木与交通工程学院，北京100044）

1 引言

砂土作为自然界中常见的颗粒介质， 具有复杂的弹塑性变形特性， 且在其变形过程中， 表现出显著的应力路径相关性。 天然地基下的饱和砂土通常会受到复杂的加载卸载而发生变形，因此，研究砂土在不同应力路径下的变形特性，并应用于工程实际， 对于土力学的发展具有重要的理论意义和应用价值。

1976 年，LAMBE[1]首次提出应力路径方法（stress path method），定义应力路径概念为“以单元划分土体，土体单元最大剪切应力点从一种应力状态变换到另一状态的轨迹”。 随后，国内外学者围绕应力路径相关问题进行了大量试验研究，得出了一系列重要结论。 BRETH[2]研究了泊松比和初始弹性剪切模量受应力路径的影响规律。 LADE 和DUNCAN[3]的研究通过观察卸载再加载的过程，表明砂土应力-应变曲线受应力路径的影响。刘祖德等[4]通过对多种应力路径条件下砂土和黏土的三轴试验， 揭示了填土的应力应变特性受应力路径影响。沈珠江等[5]通过对4 种应力路径下标准砂的三轴试验结果分析， 验证了应变总量和应变增量方向的唯一性假设的适用性。

在砂土本构模型的研究方面，姚仰平等[6]在UH 模型（Unified Hardening，统一硬化模型）的基础上，引入压硬性参量、剪胀性参数、临界状态参量，建立了砂土的UH 模型，并通过数值模拟验证了模型的合理性， 该模型不同于依靠试验数据直接建立的本构模型，具有牢固的力学基础，易于被进一步扩展，但是循环荷载、部分排水等复杂应力路径对砂土应力应变关系的影响没有考虑。同时，罗汀等[7]提出双硬化概念，将影响砂土塑性变形的硬化参量分为耦合与非耦合两部分， 提出了饱和砂土的双硬化本构模型。 该模型可以分别考虑应力路径和应变路径对砂土应力应变特性的影响，除排水条件外，还能预测不排水、部分排水和部分吸水条件下的应力应变关系，但该模型对于应力路径从不同初始状态点出发对砂土应变影响规律的研究还未考虑。

不同应力增量方向是应力路径的一种特殊形式。 砂土从同一状态点出发，沿不同应力增量方向进行加载时，砂土的塑性应变不同， 这就是不同应力增量方向对砂土变形特性的影响。 对于岩土材料而言，应力路径影响塑性变形，弹性变形仅受应力状态的影响。 王硕等[8]提出的双屈服面本构模型，将传统双屈服模型中的应力洛德角方向的剪切屈服面和体积屈服面分别用作双屈服面， 该模型能够反映应力增量方向与塑性应变增量方向具有一定的相关性， 但对于二者之间具体的对应关系还需进一步研究。施维成等[9]针对传统弹塑性模型的塑性势面应用于岩土材料中的缺陷， 提出了连续渐变塑性势面模型和塑性势面区域的概念， 进一步反映了应力增量方向和塑性应变增量方向的相关性。 郭云鹏等[10]对不同应力增量方向这一特殊应力路径进行了福建标准砂的三轴试验， 得到了福建标准砂在不同应力增量方向时的变形规律。 以上的理论和试验研究表明， 应力增量方向对于砂土变形特性具有一定影响，但是相应的理论模型尚未建立。 因此，本文针对不同应力增量方向这一特殊应力路径进行研究， 考虑孔隙比对于饱和砂土塑性应变影响，以文献[10]中开展的关于福建标准砂的三轴试验结果为基础， 建立可以表示不同应力增量方向下饱和砂土塑性应变的理论模型，并使用MATLAB 软件进行模型的适用性分析。

2 饱和砂土变形特性的应力路径相关性研究

2.1 砂土变形的细观机制

由于砂土材料在变形过程中受各种复杂因素的影响，想要通过一种方式全面考虑各种因素对砂土变形特性的影响是不现实的，只能通过研究主要影响因素得出一定的结论。 姚仰平等[11]对日本Toyoura 砂三轴试验结果进行分析，得出结论：应力比η（η=q/p。p 为压力；q 为剪力。）是造成应力路径对塑性应变影响的主要因素。 文献[12]从砂土变形的细观机制出发，给出砂土变形细观机制简图（见图1）。简图取4 个单位长度的圆形钢性棒作为一个单元体进行研究， 分析了其变形的细观机制，定性地研究了砂土变形的宏观现象。

如图1 所示， 当颗粒单元体由状态1 变到状态2 再变到状态3 时，孔隙比先增大后减小，对应宏观砂土土体先剪胀后剪缩。设与主应力σ1相对应的外力为F1、F2，与主应力σ3相对应的外力为F3、F4，如图2 所示，其中F1=F2，F3=F4。

图2 颗粒单元体受力分析简图

2.1.1 主应力比σ1/σ3与θ 之间的关系及孔隙比e 的确定

文献[12]根据颗粒的受力平衡条件、摩擦定律和牛顿定律，推导出主应力比σ1/σ3与θ 之间的关系，其中θ=∠O1O2O3，如图2 所示。该角度大小可以描述单元体的结构构成，且反映单元体的应力比和孔隙比的大小，其与主应力比之间的关系如下：

式中，tanφm为颗粒表面摩擦系数， 对于给定砂土，tanφm为常量。

文献[12]根据颗粒单元体计算孔隙比e。假设单位长的4 根圆棒轴心线围成的柱体为单元体，半径为r，单元体总体积v为：

此时圆棒体积为v1：

因此，根据式（2）和式（3）可得孔隙比为：

式中，v2为孔隙体积。

可得出：

2.1.2 应力比η 与孔隙比e 的关系

根据三轴对称情况可得：

由式（6）可以看出，应力比η 与主应力比σ1/σ3有着一一对应的关系。 由式（1）可以看出，主应力比σ1/σ3与θ 相关；由式（4）和式（5）可以看出，孔隙比e 与θ 相关。 由此得出结论：孔隙比e 与主应力比σ1/σ3一一相关， 又根据主应力比σ1/σ3与应力比η 的相关关系， 可得出孔隙比e 与应力比η 的对应关系。 因此，宏观砂土发生塑性变形的主要原因可理解为：颗粒单元体中θ 发生改变，导致孔隙比e 发生变化。 由此可得，砂土变形特性的应力路径相关性主要原因在于孔隙比发生改变。

密砂和中密砂的变形过程是先发生剪缩后发生剪胀，先硬化后软化。 根据式（1）、式（5）和式（6）推算出砂土发生剪缩和剪胀时，应力比η 与θ 的关系式函数f，函数f 体现了应力比η 与θ 的对应关系：

根据式（7）可求出某一应力比η 对应的θ，再由式（4）可知此时对应孔隙比es，应力比与孔隙比关系见式（8）：

2.2 不同应力路径下饱和砂土应力应变特性理论模型的研究

2.2.1 饱和砂土的双硬化弹塑性本构模型

罗汀等认为， 在运用弹性理论计算岩土材料的弹性应变时，球应力仅引起弹性体应变增量，剪切应力仅引起弹性剪应变增量，而岩土材料的塑性应变增量不具有简单的对应关系，要考虑球应力与剪切应力分别对剪切和压缩的耦合作用，即球应力不仅会引起塑性体应变增量， 也会引起塑性剪应变增量，剪切应力不仅会引起塑性剪应变增量，也会引起塑性体应变增量。 剑桥模型通过建立剪胀方程表示这种耦合作用。NAKAI[13]的试验表明，复杂应力路径下剪胀方程与应力路径密切相关。 因此，为考虑应力应变特性的应力路径相关性，假定塑性体应变分为与塑性剪应变有耦合作用和无耦合作用两部分。 在UH 模型的基础上将硬化参量H 分为耦合硬化参量Hc和非耦合硬化参量Hunc，如式（9）所示[14]：

其中：

式中，M 为特征状态点应力比；Mf为峰值应力比；为耦合塑性体应变增量；为非耦合塑性体应变增量；p0为初始应力；e0为初始孔隙比；μ 为耦合系数，为土的特性指标，可根据不排水条件下的特殊应力路径求出。

由此， 塑性应变增量可以表示为由非耦合硬化产生的非耦合塑性应变和耦合硬化产生的耦合塑性应变二者之和：

又通过对式（11）进行微分得非耦合塑性体应变增量为：

分析式（14），影响非耦合塑性体应变增量的因素只有平均正应力，与剪应力不相关。而非耦合塑性剪应变增量为0，即：

故塑性应变为：

此时剪胀方程为：

将式（14）～（16）代入剪胀方程式（17），得到该模型的总剪胀方程为：

由式（18）可以看出，总剪胀方程中由耦合硬化引起的剪胀方程部分与应力路径无关， 但是非耦合硬化引起的塑性体应变增量与塑性剪应变增量之比体现了应力路径相关性。

根据UH 模型[15]，将分解后的耦合硬化参量[式（10）]、非耦合硬化参量[式（11）]两部分代入其屈服函数中，得双硬化弹塑性本构模型的屈服函数如下：

式中，λ 为压缩系数；κ 为回弹系数。

2.2.2 不同应力路径下饱和砂土塑性应变理论模型的建立

文献[7]提出的饱和砂土双硬化弹塑性本构模型，只能在确定的初始应力状态点下推导出不同应力路径下砂土的应力应变关系，而对于任意初始应力状态点下饱和砂土的应力应变特性尚未给出具体处理方法。 根据本文2.1 节砂土变形的细观机制分析可知，应力比的变化引起砂土变形的细观原因在于不同应力比条件下孔隙比发生变化，因此，孔隙比是影响砂土变形的主要原因。 本文考虑孔隙比因素的影响建立能够描述砂土塑性应变的理论模型。

根据式（7）、式（8）可知，应力比与孔隙比具有一一对应关系，当得知某一状态点的应力比，就可推算出此应力比对应的孔隙比es。 为了考虑孔隙比es变化对于材料塑性变形的影响，本文将由式（7）、式（8）得出的与应力比对应的孔隙比es代替初始孔隙比e0引入双硬化弹塑性本构模型的屈服函数式（19）中，得出新屈服函数：

式中，ps为任意状态点对应的正应力；es为ps对应的孔隙比。

改进屈服函数后耦合硬化引起的塑性应变增量计算如下。

对式（20）进行全微分：

得出改进后双硬化本构模型塑性应变增量：

由式（7）、式（8）计算任意状态应力比对应的孔隙比的方法，相对于运用压缩回弹曲线计算孔隙比更加简便。 上述改进后的双硬化弹塑性本构模型， 可以表示从任意初始状态点出发不同应力路径下饱和砂土的应力应变关系， 为研究不同初始状态点或不同应力路径下饱和砂土的应力应变行为提供理论基础。

3 模型预测及与试验结果比较

文献[10]对饱和福建标准砂进行了固结完全排水三轴试验，加压属于低围压范围。进行本文模型预测时以工况4（对应初始应力p=450 kPa，q=500 kPa）为例，取应力方向角为、两种情况进行模拟比较。 文献[10]中饱和福建标准砂基本土性参数为：比重G=2.60 g/cm3，泊松比υ=0.32，等向压缩指数λ=0.0089，回弹指数κ=0.00667，峰值应力比Mf=1.45，特征状态点应力比M=1.16，取颗粒表面摩擦系数tanφm=0.4；工况4 中B 点对应应力比η=1.11。 运用MATLAB 软件编程，代入上述参数，对本文提出的新模型方程式（22）进行数值分析，得出结果如下：图3、图4 分别为应力增量方向角为120°、145°两种情况下， 塑性剪应变和塑性体应变随应力比变化的本文模型预测结果与试验结果比较图。

图3 关系结果对比

图4 关系结果对比

由图3 和图4 可知， 本文模型模拟预测结果与饱和福建标准砂的三轴试验结果基本吻合。 对比结果说明，本文提出的考虑不同应力比对应不同孔隙比的模型能够合理预测饱和砂土在不同应力路径下的塑性应变。

4 结论

本文通过分析砂土变形的细观机制， 在砂土的双硬化弹塑性本构模型基础上，考虑到应力比与孔隙比的对应关系，提出一种新的砂土弹塑性本构模型， 并与福建标准砂三轴试验结果进行了对比。 主要结论如下：

1）通过简化砂土细颗粒单元体，得到宏观应力比与孔隙比之间的对应关系。 相较于利用压缩回弹曲线计算孔隙比，该方法更简便、直接。

2）本文模型的数值模拟、预测结果与试验结果符合，本文模型能够合理预测饱和砂土在不同应力路径下的塑性变形。

现有的关于饱和砂土变形特性应力路径依赖性理论模型的研究，更多的是针对某几种特定的比较简单的应力路径，如常规三轴不排水和三轴排水应力路径， 对于其他复杂应力路径下砂土变形特性的理论模型研究还需进一步深入。 另外，考虑其他因素，如动静荷载、高低围压等对于不同应力路径下砂土变形特性的影响也值得更加深入地研究。