0.   引言

在岩土工程建设与运营中，地基中土单元往往发生主应力轴循环旋转，比如离岸岩土结构受波浪荷载^[1]，路基受交通荷载作用等^[2]。大量试验表明，主应力轴循环旋转会使产生极其复杂的塑性效应，一方面，导致土体产生累积塑性变形，排水条件下产生累积塑性压缩体变^[3]，不排水条件下产生累积孔压^[4-5]；另一方面，主应力轴循环旋转将诱发显著的变形刚度弱化与非共轴性^[5-8]，即材料当前的主应力方向与塑性主应变率方向不一致。主应力轴循环旋转下的复杂塑性效应严重威胁着离岸岩土结构、路基等的安全建设与运营，因此，合理地考虑主应力轴循环旋转塑性效应具有重要的理论与实际意义。

传统的弹塑性本构理论，绝大多数模型都隐含了共轴性(应变率与应力的主方向一致性)的假设，无法模拟主应力循环旋转诱发的非共轴性。为模拟非共轴性，一些模型引入了非共轴流动法则^[9-10]。但由于纯主应力轴循环旋转加载中主应力大小保持不变，在应力不变量表示的空间中不会产生屈服，因此，这类模型无法模拟纯主应力轴循环旋转诱发的塑性变形和非共轴性。为模拟主应力轴循环旋转塑性效应，需引入额外的力学机制和本构框架。Wang等^[11]在边界面亚塑性框架基础上构建了砂土本构模型，模型中塑性应变率方向与应力率方向相关，模拟了砂土在纯主应力轴循环旋转下的不排水有效应力路径。基于下加载面模型，Tsutsumi等^[12]引入了能产生非共轴性的各向异性和切向加载效应。Li等^[13]建议了一种各向异性临界状态为基础的本构理论框架，模型引入了与组构各向异性及非共轴流动相关的非比例加载机制。基于这种本构框架，Gao等^[14]通过引入组构及其演化规律，建立了三维本构模型模拟砂土主应力轴循环旋转下的非共轴性，童朝霞等^[15]则将主应力幅值变化以及应力主轴旋转产生的塑性变形单独加以考虑，建立了可考虑应力主轴循环旋转效应的砂土本构模型。而Lashkari等 ^[16]在Gutierrez等^[17]总结的应力空间内建立了相应的边界面，并引入基于主应力旋转试验规律的非共轴流动方向，通过考虑土体组构和非共轴流动方向，模拟了砂土的非共轴变形特征。Yang等^[18]在一般应力空间建立了一种包含切向加载的考虑主应力轴旋转的运动硬化模型。随后，Tian等^[19]从材料各向异性的角度将UH模型推广至模拟土的主应力轴循环旋转效应，通过引入组构及其演化规律考虑了固有和应力诱发各向异性，模拟了土的纯主应力轴循环旋转效应。另外，陈洲泉等^[20]在Qian等^[21]的基础上，重新定义了非共轴流动方法，并考虑了边界面塑性映射法则和组构各向异性，模拟了砂土主应力轴旋转的非共轴性。应该说，上述研究采用不同的力学机制实现了土的共轴变形特性的本构模拟，但研究对象仍局限于砂性土。相比较，最新试验观察发现^[5-8]主应力偏转路径下，天然黏性土的非共轴特性更为复杂，如纯主应力轴循环旋转下孔压循环波动的累积效应、剪切刚度的循环弱化效应等，并且这些非比例加载变形特性与非共轴性有一定关系^[6-7]，因此，为更准确地实现对这些复杂的力学响应的模拟，有必要对既有非共轴本构模型加以改进与完善。

基于以上认识，本文建立了能够描述主应力轴循环旋转塑性效应的边界面模型。模型通过引入可移动映射法则，描述主应力轴循环旋转卸载情况下的塑性变形。通过考虑固有各向异性弹性，描述天然软黏土主应力轴循环旋转下循环波动的塑性累积行为。同时，修正了非共轴流动法则考虑循环过程中非共轴性变化。最后，通过对比温州软黏土的模型预测和试验结果，验证了模型的有效性。

1.   边界面模型的建立

根据弹塑性力学的基本假定，总应变率${\dot \varepsilon _{ij}}$可分解为弹性部分和塑性部分，即

式中，弹性应变$\dot \varepsilon _{ij}^{\text{e}}$可由弹性体积模量K和弹性剪切模量G来计算。根据剑桥模型的假定，K和G由当前的球应力p确定：

式中，$\nu$为泊松比，${\kappa ^*} = \kappa /(1 + {e_0})$，${e_0}$为初始孔隙比，$\kappa$为e–lnp空间中回弹曲线的斜率。增量型弹性应力应变关系为

式中，$\boldsymbol{D}_{i j k l}^{\mathrm{e}}$为弹性刚度矩阵，可由K和G计算。

1.1   各向异性边界面方程

Ling等^[22]建议的单面边界面模型能够合理描述软黏土各向异性应力应变特性^[23-24]，如图 1所示，本文仍沿用这一方法来定义各向异性边界面方程。

图  1  各向异性边界面

Figure  1.  Anisotropic bounding surface

下载: 全尺寸图片 幻灯片

考虑各向异性的边界面方程如下：

式中，$\bar p$和${\bar q_\alpha }$分别为平均应力和折减广义剪应力，${\bar q_\alpha } = \sqrt {3\bar s_{ij}^\alpha \bar s_{ij}^\alpha /2}$，其中$\bar s_{ij}^\alpha = {\bar s_{ij}} - {\bar \sigma _{kk}}{\alpha _{ij}}/3$为折减应力偏量；${\alpha _{ij}}$为各向异性张量，表示土体各向异性的大小，$\alpha = \sqrt {3{\alpha _{ij}}{\alpha _{ij}}/2}$为p–q空间中屈服面的倾角（如图 1所示），M为临界状态应力比，R为形状参数，参数的定义及其物理意义参见黄茂松等^[23]和Huang等^[24]。

1.2   可移动映射法则

本文关注的纯主应力轴循环旋转加载，又称纯环剪，是非比例加载的一种，加载中主应力方向连续循环旋转，但广义剪应力q和主应力大小保持不变^[25]。由于主应力保持不变，因此，在π平面中表示为边界面上一个位置不动的点。为更好地表示这种纯主应力轴旋转的情况，本文将其表征在偏平面${\tau _{z\theta }}$–${\text{(}}{\sigma _z} -$ ${\sigma _\theta }{\text{)/2}}$上，其中${\sigma _z}$，${\sigma _\theta }$和${\tau _{z\theta }}$分别为轴向应力、环向应力和剪切应力。在偏平面上，屈服面表示为绕坐标原点的圆，而纯主应力轴旋转应力路径如图 2所示。加载过程中，应力增量方向与屈服面外法线方向夹角始终大于90°，即${\text{(}}\partial f/\partial {\sigma _{ij}}{\text{)d}}{\sigma _{ij}}$ < 0，无塑性变形产生，为卸载情况。值得一提的是，这种加载因子小于0是普遍存在的，如超固结土在循环荷载条件下的变形^[26-27]，为描述这种卸载情况下塑性变形的产生，本文借鉴Wang等^[11]和Li^[28]提出的可移动映射法则，其在π平面上的表示见Li^[28]，为更好地表示这种映射法则对纯主应力轴旋转的情况，本文将其表征在偏平面上。

图  2  偏平面τ_zθ–(σ_z-σ_θ)/2上纯主应力旋转应力路径

Figure  2.  Stress path of pure principal stress rotation in τ_zθ–(σ_z-σ_θ)/2 space

下载: 全尺寸图片 幻灯片

根据边界面塑性理论^{[24, 29]}可知，塑性加载因子定义为

式中${{\partial F} \mathord{\left/ {\vphantom {{\partial F} {\partial {{\bar \sigma }_{ij}}}}} \right. } {\partial {{\bar \sigma }_{ij}}}}$表示塑性加载方向，为边界面虚应力处外法线方向；${H_{\text{p}}}$与${\bar H_{\text{p}}}$分别为当前应力点与虚应力点塑性模量，当前应力点与虚应力点重合时(即${\bar \sigma _{ij}} = {\sigma _{ij}}$)，${H_{\text{p}}} = {\bar H_{\text{p}}}$。

定义加载方向从加载转变为卸载（塑性加载因子$\dot L$从正转变为负）时的上一步的应力点作为可移动的映射中心${\omega _{ij}}$，以连接映射中心${\omega _{ij}}$和当前应力点${\sigma _{ij}}$的直线与边界面的交点作为虚应力点${\bar \sigma _{ij}}$，如图 3所示，则构成了可移动映射法则。映射中心初始位置在偏平面的原点O。图 3中所示为一般加载情况，塑性加载因子$\dot L$＞0，则映射中心${\omega _{ij}}$保持不变。而对于一般卸载情况，如图 4（a）所示，塑性加载因子$\dot L$＜0，此时，映射中心${\omega _{ij}}$移动至上一步的应力点位置，如图 4（a）所示，映射中心${\omega _{ij}}$移动后，以连接映射中心${\omega _{ij}}$和应力增量${\text{d}}{\sigma _{ij}}$的直线与边界面的交点作为虚应力点${\bar \sigma _{ij}}$，如图 4（b）中分析，塑性加载因子重新变为$\dot L > 0$。即卸载时加载因子为负会导致映射中心${\omega _{ij}}$移动，而通过映射中心${\omega _{ij}}$移动会重新导致加载因子为正，这种可移动映射法则保证加卸载时有塑性变形产生。

图  3  偏平面τ_zθ–(σ_z-σ_θ)/2上加载映射法则

Figure  3.  Mapping rule in τ_zθ–(σ_z-σ_θ)/2 space during general loading situation

下载: 全尺寸图片 幻灯片

图  4  映射中心移动机制

Figure  4.  Relocatable mechanism of projection center

下载: 全尺寸图片 幻灯片

图 3，4中$\rho$为映射中心${\omega _{ij}}$到当前应力点${\sigma _{ij}}$的距离，$\bar \rho$为映射中心${\omega _{ij}}$到虚应力点${\bar \sigma _{ij}}$的距离，映射法则表示为

1.3   各向异性弹性的引入

由于天然软黏土的沉积使其呈现出横观各向同性弹性，又称各向异性弹性^[30-31]。主应力轴循环旋转下也会产生这种各向异性弹性效应，纯主应力轴循环旋转试验表明，不排水条件下的孔隙水压和排水条件下的体应变均呈波动形式累积^{[3-6, 8]}，如图 5所示。根据体应变可分为不可逆和可逆两种的定义^[32]，不排水试验中孔隙水压波动的下限对应于排水试验中永久累积塑性不可逆体应变，而上限和下限的差即为弹性波动可逆的孔隙水压，如图 5所示。早期的研究表明，孔隙水压波动这种塑性累积行为可能归因于天然软黏土弹性行为的各向异性^[30-31]。

图  5  纯主应力轴循环旋转孔隙水压累积中呈现的各向异性弹性

Figure  5.  Anisotropic elasticity in accumulation of excess pore pressure with cycles of pure principal stress rotation

下载: 全尺寸图片 幻灯片

为考虑这种波动塑性累积行为，模型中引入各向异性弹性，可表示为

式中，${E_{\text{h}}}$和${E_{\text{v}}}$分别为水平和竖向弹性模量，${\nu _{{\text{vh}}}}$和${\nu _{{\text{vv}}}}$分别为水平和竖向泊松比，${G_{{\text{vh}}}}$为剪切模量，满足${E_{\text{h}}} = {n^2}{E_{\text{v}}}$，${\nu _{{\text{vh}}}} = n{\nu _{{\text{vv}}}}$，${G_{{\text{vh}}}} = {{n{E_{\text{v}}}} \mathord{\left/ {\vphantom {{n{E_{\text{v}}}} {\left[ {2(1 + n{\nu _{{\text{vv}}}})} \right]}}} \right. } {\left[ {2(1 + n{\nu _{{\text{vv}}}})} \right]}}$，其中，$n$为各向异性因子。

对于应力控制的各向同性压缩，根据上述方程可知：

考虑到体积模量K可根据式(2)计算得到，竖向弹性模量可表示为

不同于各向同性弹性模型回弹指数$\kappa$、泊松比$\nu$，在各向异性弹性中还需要一个各向异性因子$n$。根据以上分析可知，${n^2}$为水平与竖向刚度之比，当$n$=1时，材料为各向同性。

1.4   非共轴流动法则

Qian等^[21]将非共轴应力率在当前应力状态方向上进行正交投影，给出了广义应力空间内的非共轴应力率：

式中，${S_{ij}} = {s_{ik}}{s_{kj}} - {{2{J_2}{\delta _{ij}}} \mathord{\left/ {\vphantom {{2{J_2}{\delta _{ij}}} 3}} \right. } 3} - {{3{s_{ij}}{J_3}} \mathord{\left/ {\vphantom {{3{s_{ij}}{J_3}} {2{J_2}}}} \right. } {2{J_2}}}$。相应的非共轴应变率定义为

式中，${H_{\text{t}}}$为塑性模量。

大量黏土试验结果表明^{[8, 33-34]}，随着当前剪切应力增大，非共轴应变率与应力方向不一致逐渐减小，至临界状态时趋于共轴，这种非共轴特性不能被方程（11）定义的非共轴特性反映^{[16, 20]}。另外，方程（11）中的非共轴参数为切向塑性模量，在循环剪切过程中，切向塑性模量呈现显著的弱化效应，若忽略其弱化效应则无法准确地预测非共轴塑性变形。因此，为更合理地模拟非共轴流动性，本文修正了Qian等^[21]非共轴流动理论，具体表述为

式中，$n_{ij}^{nc}$为和$\dot s_{ij}^n$具有相同方向的非共轴应力率，表述为

其中，$n_{ij}^{\dot s} = {{{{\dot s}_{ij}}} \mathord{\left/ {\vphantom {{{{\dot s}_{ij}}} {\left\| {{{\dot s}_{ij}}} \right\|}}} \right. } {\left\| {{{\dot s}_{ij}}} \right\|}}$为单位偏应力增量，$\eta$和$M$分别为当前应力比和临界状态应力比，当$\eta$=$M$时，$\dot \varepsilon _{ij}^{np} = 0$，$\gamma$，${k_{\text{n}}}$均为材料参数，控制非共轴塑性应变随当前应力比变化程度及量级，$\dot L$为塑性加载因子。

Macaulay括号内表示随着应力比增加，非共轴程度减小，至临界状态时，应变增量方向与应力方向共轴，这与试验结论是符合的。值得一提的是，与Qian等^[21]非共轴不同的是，方程（13）中非共轴应变率与共轴应变率用了相同的塑性加载因子，这是由于非共轴与共轴应变率均为不可恢复塑性变形，由相同的塑性机制驱动，与Li等^[35]讨论的一致，仅仅是方向不同。另外，虽然用相同的塑性加载因子，但非共轴与共轴应变率强度却不同，$\gamma$和${k_{\text{n}}}$共同控制非共轴塑性应变相对于共轴塑性应变的量级。通过以上分析可知，本文修正使非共轴应力率无量纲化，然后非共轴流动可以和共轴流动通过塑性加载因子耦合起来，并考虑应力比影响。

共轴塑性应变采用相关联流动法则，即塑性流动方向与塑性加载方向重合，且等于边界面在虚应力点处的外法线方向。考虑共轴塑性应变和非共轴塑性应变，可得总的流动法则：

1.5   增量型应力应变关系

模型采用两个独立的硬化内变量${p_{\text{c}}}$以及${\alpha _{ij}}$, 其中${p_{\text{c}}}$反映应力历史的影响，${\alpha _{ij}}$定义屈服面倾角的大小，反映各向异性的程度。硬化规律的具体计算见Huang等^[24]。边界面虚应力点的塑性模量${\bar H_{\text{p}}}$可根据一致性条件求出：

当前应力点的塑性模量${H_{\text{p}}}$可由插值函数计算^[24]：

式中，${p_{\text{a}}}$为标准大气压（101.325 kPa），$\varepsilon _{\text{s}}^{\text{p}} = \int {\dot \varepsilon _{\text{s}}^{\text{p}}}$为塑性偏应变的累积值，$\zeta$，${\psi _0}$，$\xi$为模型参数。

增量型弹性应力应变关系见式（3），$\boldsymbol{D}_{ijkl}^{\text{e}}$为弹性刚度矩阵，可由式（7）计算，将式（1），（3），（5），（14）代入一致性条件得

最终可得如下弹塑性应力应变增量关系：

其中，$\boldsymbol{D}_{ijkl}^{{\text{ep}}}$为弹塑性刚度矩阵，表达为

本文所有计算都通过Bardet等^[36]提出的积分算法实现。

1.6   模型参数

本文模型共需要14个参数和初始状态参数${e_0}$，如表 1所示，模型相关参数根据其功能可分为6组。相对于黄茂松等^[23]和Huang等^[24]采用的边界面模型，模型中减少了结构性参数，引入了非共轴参数$\gamma$、${k_{\text{n}}}$和各向异性弹性参数$n$，除上述3个参数以外的其他参数，其标定方法可以参考黄茂松等^[23]和Huang等^[24]的论述，本文不再重复。非共轴参数$\gamma$控制随着应力比增加非共轴程度减小的速率，${k_{\text{n}}}$控制非共轴塑性应变相对于共轴塑性应变的大小。这两个非共轴特性参数$\gamma$和${k_{\text{n}}}$可由两组${\alpha _\sigma }$≠0°和${\alpha _\sigma }$≠90°的主应力轴固定单调剪切试验来标定，根据式（12）可得

表  1  模型参数

Table  1.  Model parameters

临界状态参数	形状参数	硬化参数	非共轴参数	各向异性弹性参数	边界面插值参数
$\lambda$=0.36 $\kappa$=0.04 ${M_{\text{e}}}$=0.81 ${M_{\text{c}}}$=1.11 $\nu$=0.2	R=3.2	$\mu$=60 $\beta$=0.7	$\gamma$=0.5 ${k_{\text{n}}}$=50	n=1.1	$\zeta$=40.0 ${\psi _0}$=12.0 $\xi$=6.0

下载: 导出CSV 

| 显示表格

式中，(${\text{d}}\varepsilon _{z\theta }^{{\text{p1}}}$, ${\eta _1}$)和(${\text{d}}\varepsilon _{z\theta }^{{\text{p2}}}$, ${\eta _2}$)为${\text{d}}\varepsilon _{z\theta }^{\text{p}}$–$\eta$曲线上两个不同的点，参数${k_{\text{n}}}$可在$\gamma$确定后通过拟合其中一条$\varepsilon _{z\theta }^{\text{p}}$–$\eta$曲线来标定。各向异性弹性参数$n$为各向异性因子，是天然黏土的各向异性弹性参数。${n^2}$为水平与竖向刚度之比，可从固结后水平与竖向切片三轴试验或弯曲元试验获得的弹性模量或泊松比之比来标定。如果试验数据不足，也可通过试错法模拟可用数据，然后通过最佳拟合来确定，特别是通过模拟不排水试验孔压波动曲线或排水试验中体应变波动曲线来标定。

2.   模型验证

模型通过温州软黏土在纯主应力轴循环旋转下的不排水空心圆柱扭剪试验进行验证。Wang等^[8]对温州天然软黏土开展了一系列纯主应力轴循环旋转试验，首先对试样进行各向同性固结，初始固结压力p′=150 kPa，然后增加剪应力q至固定值，保持各向异性固结，随后进行纯主应力轴循环旋转不排水试验，如2.2节中描述，保持平均主应力p、中主应力系数b和剪应力q不变，连续旋转主应力方向，试验具体条件和方案见Wang等^[8]，试验初始条件为${p_0}$=100 kPa和${e_0}$=1.570。土体计算参数如表 1所示，需要说明的是，该试验并未提供临界状态参数$\lambda$和$\kappa$，根据Nakase等^[37]的研究，临界状态参数$\lambda$和$\kappa$与土体物理性质指标液限w_L或塑性指数I_P存在线性关系，但其关系需根据土的特性来确定，参考王立忠等^[38]对温州黏土的研究，对比其土体物理性质指标液限w_L和塑性指数I_P，综合确定的温州软黏土临界状态参数见表 1。

纯主应力轴循环旋转下，温州饱和软黏土不排水试验的孔隙水压变化与模拟结果如图 6所示。由图 6可知，虽然每个循环广义剪应力q保持不变，但孔压随主应力轴旋转呈波动形式累积，本文本构模型能够合理模拟孔压波动这种塑性积累行为。由图 6还可看出，b=0的孔压累积小于b=0.5，1.0，在循环的初级阶段，孔压累积较快，随循环数增加，孔压累积变慢，而本文模型也能合理模拟中主应力系数和循环阶段对孔压累积的影响。图 6中虚线表示b=1.0条件下采用各向同性弹性n=1的模拟，即不考虑各向异性弹性，结果显示没有孔压波动。通过对比可知，孔压波动这种塑性累积行为的模拟要归结于本构模型中引入的各向异性弹性。

图  6  孔隙水压累积试验与模拟

Figure  6.  Tests and simulation of accumulation of pore water pressure

下载: 全尺寸图片 幻灯片

图 7为温州软黏土在纯主应力轴循环旋转下b=1.0系列的应变分量随主应力轴旋转角度的关系。由图 7可知，尽管广义剪应力保持不变，纯主应力轴旋转仍然产生了显著的塑性应变累积，模型合理的模拟这种塑性应变累积，这主要由于边界面模型中采用了可移动映射法则。

图  7  应变分量变化规律试验与模拟

Figure  7.  Tests and simulation of variation of strain components

下载: 全尺寸图片 幻灯片

纯主应力轴循环旋转b=1.0下，温州饱和软黏土不排水试验的剪切应力应变变化与模拟结果如图 8所示，图 8（a），8（b）分别为5个循环的试验和模拟结果。如图 8（a）所示，尽管广义剪应力保持不变，剪切应力应变曲线仍呈现滞回特性并伴随着累积塑性的发展，滞回圈刚度显著弱化且滞回圈形状呈开放型性发展，对比图 8（a），8（b）可知，尽管有一定差距，但模型合理地模拟了试验中这种剪切应力应变关系。

图  8  剪切应力应变关系试验与模拟

Figure  8.  Tests and simulation of shear stress strain relationship

下载: 全尺寸图片 幻灯片

图 9为温州软黏土在纯主应力轴循环旋转下b=0.5的第三个循环（N=3）应变增量矢量试验和模拟结果，相应的非共轴角随主应力轴旋转角变化的试验和模拟如图 10所示。由图 9可知，应变增量在循环加载过程中具有明显的分段特征，当${\alpha _\sigma }$在[0°，45°]和[90°，135°]变化时，非共轴角度减小，而当主应力轴旋转角${\alpha _\sigma }$在[45°，90°]和[135°，180°]变化时，非共轴角增大，由模拟结果可知，模型计算出的分段特征大致与试验数据吻合。由图 10的对比结果可知，模型计算的非共轴角随循环变化规律也与试验结果一致，说明模型能合理模拟黏土随主应力轴循环旋转表现的非共轴性。图 10还给出了不考虑非共轴流动法则${k_{\text{n}}}$=0，即共轴模型的预测结果，由于采用的可移动映射法则的映射中心并不在坐标中心，因此塑性流动方向并不等于应力方向角${\alpha _\sigma }$，所以非共轴角并不等于0，但是几乎不随主应力轴旋转变化。因此，对比可知本文引入的修正非共轴流动法则在非共轴性预测中起显著作用，合理模拟了非共轴角随循环变化的规律。

图  9  非共轴特性试验与模拟(N=3)

Figure  9.  Tests and simulation of noncoaxial behavior (N=3)

下载: 全尺寸图片 幻灯片

图  10  非共轴角试验与模拟

Figure  10.  Tests and simulation of noncoaxial angle

下载: 全尺寸图片 幻灯片

3.   结论

本文在考虑各向异性边界面模型的基础上，通过引入可移动映射法则，处理主应力轴循环旋转卸载下的塑性效应。通过考虑固有各向异性弹性，模拟主应力轴循环旋转下循环波动的塑性累积行为。同时，修正了非共轴流动法则反映非共轴性的循环变化，建立了一种考虑主应力轴循环旋转塑性效应的本构模型。最后，对温州天然饱和软黏土在纯主应力轴循环旋转下的不排水行为进行了模拟，得到以下3点结论。

（1）在纯主应力轴循环旋转试验中，因为模型引入了可移动映射法则，使模型可描述这种卸载情况下的塑性变形。

（2）由于模型中引入了固有各向异性弹性，从而模型可以描述纯主应力轴循环旋转中循环波动这种塑性累积行为。

（3）修正的非共轴理论中考虑了应力比影响并将非共轴流动和共轴流动耦合起来，从而可以描述纯主应力轴循环旋转中非共轴性随循环变化规律。