前海软土蠕变特性试验及模型参数研究

梁荣康; 张庆军; 张成; 李栋; 苏栋; 李强

doi:10.11779/CJGE2021S2032

前海软土蠕变特性试验及模型参数研究 English Version

梁荣康^1,,
张庆军²,
张成³,
李栋⁴,
苏栋^{1, 5, ,},
李强¹

1.
深圳大学土木与交通工程学院,广东　深圳　518060
2.
中铁十五局集团有限公司,上海　200070
3.
珠海大横琴城市新中心发展有限公司,广东　珠海　519030
4.
珠海市规划设计研究院,广东　珠海　519000
5.
滨海城市韧性基础设施教育部重点实验室（深圳大学）,广东　深圳　518060

基金项目:

国家自然科学基金项目 51938008

国家自然科学基金项目 52090081

中铁十五局科研项目 CR15CG-ZH-HQMZSD-2021-00001（KJ）

详细信息

作者简介:
梁荣康(1997— ),男,硕士研究生,主要从事软土的力学特性方面学习研究工作。E-mail：2060471029@email.szu.edu.cn

通讯作者:
苏栋, E-mail: sudong@szu.edu.cn

中图分类号: TU447
计量
- 文章访问数: 202
- HTML全文浏览量: 49
- PDF下载量: 99
出版历程
- 收稿日期: 2021-08-16
- 网络出版日期: 2022-12-05
- 刊出日期: 2021-10-31

Creep characteristics and model parameters of Qianhai soft soil

1.
College of Civil and Transportation Engineering, Shenzhen University, Shenzhen 518060, China
2.
China Railway 15 Bureau Group Co., Ltd., Shanghai 200070, China
3.
Zhuhai Da Hengqin City New Center Development Co., Ltd., Zhuhai 519030, China
4.
Zhuhai Institute of Urban Planning & Design, Zhuhai 519000, China
5.
Key Laboratory for Resilient Infrastructures of Coastal Cities (MOE), Shenzhen University, Shenzhen 518060, China

摘要

摘要: 为研究沿海地区软弱土层的蠕变特性并建立相应的蠕变模型,对深圳市城市轨道交通12号线海上田园东站处淤泥质黏土进行室内直剪蠕变试验。结果表明,该软土试样的蠕变特性与所受应力状态有关,在400 kPa固结压力下,当试样所受剪应力较小时,蠕变特性不显著,随着剪应力的提高,蠕变特性逐渐显著。根据蠕变速率将软土的蠕变过程区分为减速和稳定两个阶段,随着剪应力的提高,该软土在各蠕变阶段的时间均逐渐延长,相应的应变增量也逐渐增大。采用修正的Singh-Mitchell模型对该软土的蠕变特性进行建模,该模型参数少,本构方程简单,适用性较强。通过采用不同参考时间的试验数据对模型参数进行校核,发现当参考时间t₁=1 h时模型计算结果与试验数据之间的误差最小,预测精度较高。
- 软土 /
- 蠕变 /
- 直剪试验 /
- 修正Singh-Mitchell模型
Abstract: In order to explore the creep characteristics of the soft soil layer in coastal areas and establish the corresponding creep model, a direct shear creep test is conducted on the silty clay at the Haishangtianyuan East Station of Shenzhen Metro Line 12. The results show that the creep characteristics of the soft soil samples are related to the stress state. Under the consolidation pressure of 400 kPa, when the shear stress of the samples is small, the creep characteristics are not significant. With the increase of shear stress, the creep characteristics gradually become more significant. Using the creep rate, the creep process of the soft soil is divided into two stages: deceleration and stability stages. With the increase of shear stress, the elapsed time for each creep stage grows, and the corresponding strain increment increases gradually. The modified Singh-Mitchell model is used to simulate the creep characteristics of the soft soil. This model is widely used because it has fewer parameters and its constitutive equation is simple. The model parameters are calibrated using the test data at different reference time. It is found that the error between the calculated results and the test data is the smallest when the reference time t₁=1 h, and the prediction accuracy is higher.
- soft soil /
- rheology /
- direct shear test /
- modified Singh-Mitchell model

HTML全文

0. 引言

由于特殊的沉积环境和形成方式，钙质砂具有形状不规则和孔隙率高等特点，使其在较低应力水平发生颗粒破碎，进而影响构筑物等的强度和变形，在工程设计、施工中，钙质砂颗粒破碎将直接威胁建筑物或构筑物的稳定和安全。因此，颗粒破碎对钙质砂力学特性的影响及考虑钙质砂颗粒破碎的本构模型研究一直是学者和工程建设人员关注的科学问题。

张家铭等^[1]研究了钙质砂在侧限压缩条件下的破碎特性，指出导致钙质砂在常压力下压缩特性与黏土类似的原因是发生了颗粒破碎。蔡正银等^[2]和吴杨等^[3]研究了钙质砂在三轴剪切条件下的变形、强度及破碎规律随密实度和围压的演化规律，并建立了颗粒破碎和围压及初始密实度的定量关系。研究表明，钙质砂变形与颗粒重排引起的剪胀效应和颗粒破碎引起的体缩密切相关，其剪切特性被认为是初始密实度、颗粒形状、尺寸、级配、应力水平和颗粒破碎等共同影响的结果，并表现出明显的状态相关特征^[4-6]。因此，诸多因素的影响导致对钙质砂力学特性的统一描述及建立合理的本构模型带来挑战。

Coop^[7]在高压三轴条件下发现，尽管钙质砂发生显著地颗粒破碎，其力学行为仍符合临界状态土力学的基本特征，但颗粒破碎对钙质砂的临界状态产生了明显影响。Bandini等^[8]指出钙质砂颗粒破碎不仅使压缩平面临界状态线(CSL)垂直移动，还产生旋转；此外，Yu^[9]还发现在 $p'$ -q平面CSL保持不变但与排水条件相关。王刚等^[10]利用离散元方法研究了钙质砂颗粒破碎对临界状态线的影响，指出临界状态线位置由破碎的最终颗粒级配决定。因此，颗粒破碎对临界状态线的影响规律是建立本构关系的关键。

目前，已有不少学者在试验基础上从不同理论出发开展了关于钙质砂本构模型的研究。孙吉主等^[11]在状态相关理论下，建立了考虑状态相关和剪胀性的双屈服面模型。蔡正银等^[12]基于临界状态理论，建立了考虑钙质砂颗粒破碎临界状态的本构模型。曾凯锋等^[13]基于颗粒破碎对钙质砂的影响修正了邓肯-张模型。王兆南等^[14]在边界面理论框架下，建立了考虑钙质砂颗粒破碎的循环边界面本构。砂土弹塑性本构模型对于塑性流动的描述多采用非相关联流动法则，传统做法是假定一个不同于屈服面的塑性势面，对该塑性势面微分得到相应的塑性流动方向，显然增加了模型参数。而Sumelka^[15]尝试引入分数阶微分到弹塑性理论中考虑非相关联流动法则并取得较好结果。孙逸飞等^[16]利用分数阶塑性流动法则建立了砂土的本构模型。Liang等^[17]和梁靖宇等^[18]基于特征应力空间建立了分数阶临界状态模型，减少了以往模型的物理参数。由此可见，分数阶塑性流动法则具有很好的优势和应用。但基于分数阶塑性流动法则建立考虑钙质砂颗粒破碎的弹塑性本构模型尚未有尝试。

本文将基于临界状态理论，考虑颗粒破碎和密实度对临界状态的影响，建立颗粒破碎临界状态线；结合分数阶塑性流动法则及边界面塑性理论，建立考虑钙质砂状态相关及颗粒破碎特性的边界面本构模型。

1. 颗粒破碎临界状态线

土体在剪切作用下达到的稳定状态，即体变增量，应力增量或孔压增量为0，剪应变增量不为0的状态，通常被称为土体的临界状态，并被广泛应用于土的本构模型中。对于钙质砂，试样在剪切过程中，破碎程度不同的土体其级配状态、内部排列方式、接触、孔隙都发生变化，其CSL还将受颗粒破碎的影响，导致CSL的位置和形状发生改变。Russell等^[19]采用分段函数描述砂土从低应力到高应力水平的CSL并建立本构模型，但形式复杂、参数较多。Daouadji等^[20]提出用一组平行直线描述颗粒破碎对CSL的影响，并被Wood等^[21]用级配状态指数表征破碎导致的级配变化和对CSL的影响。基于此，Xiao等^[22-24]提出颗粒破碎临界状态面概念，建立了孔隙比-应力-相对破碎指数三者的关系，很好地描述了粗粒土颗粒破碎效应。Tong等^[25]采用双对数方法和参考状态线，描述颗粒破碎对CSL的影响。因此，本文基于以上方法，引入破碎应力的概念到CSL，反映加载过程颗粒破碎对试样CSL的影响，建立颗粒破碎和破碎应力的关系，得到颗粒破碎临界状态线(BCSL)：

$e_{\mathrm{cs}}=\mathit{\Pi}-\lambda \ln \left[\left(p^{\prime}+p_{\mathrm{b}}\right) / p_{\mathrm{r}}\right]\ \ 。$

(1)

式中： $e_{\mathrm{cs}}$ 为临界状态孔隙比； $\mathit{\Pi }$ ， $\lambda$ 为压缩平面CSL的截距和斜率； $p^{\prime}$ 为平均有效应力； $p_{\mathrm{b}}$ 为破碎应力，与破碎程度相关， $p_{\mathrm{b}} = p_{\mathrm{b}}(B_{\mathrm{r}})$ ， $B_{\mathrm{r}}$ 为量化土体破碎程度因子； $p_{\mathrm{r}}$ 为参考应力，取 $p_{\mathrm{r}}$ =1 kPa。

如图 1，图中虚线为BCSL，只有当颗粒破碎发生时，式（1）所表示的BCSL随颗粒破碎的增加不断向下移动，与已有试验现象较为吻合^[8-9]。

图 1 颗粒破碎临界状态线示意图

Figure 1. Schematic diagram of critical state line of breakage

下载: 全尺寸图片幻灯片

2. 状态相关非正交边界面模型

2.1 弹塑性本构关系

根据弹塑性本构关系的基本思想，将材料的变形看作体应变和剪应变之和，其中体应变和剪应变又可分别用弹性部分和塑性部分之和表示，其增量形式可表达为

$d\varepsilon = d{\varepsilon _{\mathrm{v}}} + d{\varepsilon _{\mathrm{s}}} \text{，}$

(2)

$\begin{array}{l}d{\epsilon }_{\mathrm{v}}=d{\epsilon }_{\mathrm{v}}^{\mathrm{e}}+d{\epsilon }_{\mathrm{v}}^{\mathrm{p}}\text{ }\text{，}\\ d{\epsilon }_{\mathrm{s}}=d{\epsilon }_{\mathrm{s}}^{\mathrm{e}}+d{\epsilon }_{\mathrm{s}}^{\mathrm{p}}\text{ }\ 。\end{array}\}$

(3)

式中： $d\varepsilon$ 为总应变增量； $d{\varepsilon _{\mathrm{v}}}$ 为体应变增量； $d{\varepsilon _{\mathrm{s}}}$ 为剪应变增量； $d\varepsilon _{\mathrm{v}}^{\mathrm{e}}$ ， $d\varepsilon _{\mathrm{s}}^{\mathrm{e}}$ 分别为弹性体应变增量和弹性剪应变增量部分； $d\varepsilon _{\mathrm{v}}^{\mathrm{p}}$ ， $d\varepsilon _{\mathrm{s}}^{\mathrm{p}}$ 分别为塑性体应变增量和塑性剪应变增量部分。其中弹性部分可根据胡克定律求得

$\left.\begin{array}{l} \mathrm{d} \varepsilon_{\mathrm{v}}^{\mathrm{e}}=\mathrm{d} p^{\prime} / K \quad, \\ \mathrm{~d} \varepsilon_{\mathrm{s}}^{\mathrm{e}}=\mathrm{d} q /(3 G) \quad 。 \end{array}\right\}$

(4)

式中： $K$ ， $G$ 分别为弹性体积模量和弹性剪切模量，

$\left.\begin{array}{l} K=\left(1+e_0\right) p^{\prime} / \kappa\ \ , \\ G=3 / 2(1-2 v) /(1+v) K \quad 。 \end{array}\right\}$

(5)

式中： ${e_0}$ 为材料初始孔隙比； $\kappa$ 为材料在压缩平面回弹线斜率； $\nu$ 为材料泊松比。

塑性部分按照塑性理论可由下式计算:

$\left.\begin{array}{l} \mathrm{d} \varepsilon_{\mathrm{v}}^{\mathrm{p}}=\left\langle\left(n_{\mathrm{p}} \mathrm{d} p^{\prime}+n_{\mathrm{s}} \mathrm{d} q\right) / H\right\rangle m_{\mathrm{p}},\\ \mathrm{d} \varepsilon_{\mathrm{s}}^{\mathrm{p}}=\left\langle\left(n_{\mathrm{p}} \mathrm{d} p^{\prime}+n_{\mathrm{s}} \mathrm{d} q\right) / H\right\rangle m_{\mathrm{s}^{\circ}} \end{array}\right\}$

(6)

式中：( ${n_{\mathrm{p}}},{n_{\mathrm{s}}}$ )为加载方向单位向量；( ${m_{\mathrm{p}}},{m_{\mathrm{s}}}$ )为塑性流动方向单位向量； $H$ 为归一化塑性模量； $\left\langle L \right\rangle$ 为Macaulay括号，当 $L > 0$ ， $\left\langle L \right\rangle = L$ ；当 $L \leqslant 0$ ， $\left\langle L \right\rangle = 0$ 。对于塑性流动方向单位向量、加载方向单位向量和统一塑性模量在后面小节详细介绍。

2.2 考虑破碎的边界面硬化规律

本文边界面方程 $\bar f$ 选取修正椭圆方程^[26-28]，其表达式为

$\begin{gathered} \bar{f}=M_{\mathrm{c}}^2 \beta^2\left(\bar{p}^{\prime}-\beta \bar{p}_0^{\prime}\right)^2+(1-\beta)^2 \bar{q}^2- \\ M_{\mathrm{c}}^2 \beta^2(1-\beta)^2 \bar{p}_0^{\prime 2}=0\ \ \ 。 \end{gathered}$

(7)

式中： ${M_{\mathrm{c}}}$ 为应力平面临界状态应力比； ${\bar p'_0}$ 为边界面与平均有效应力轴的交点，控制边界面的大小，其演化由塑性硬化规律决定； $\beta$ 为边界面形状的长宽比，控制边界面形状，即 $\beta = {\bar p'_A}/{\bar p'_0}$ ， ${\bar p'_A}$ 为椭圆上顶点，交于应力平面CSL；见，( $\bar p',\bar q$ )为边界面上的像应力点，对应屈服面上当前真实应力点( $p',q$ )，并定义：

$\chi = p'/\bar p' = q/\bar q \ \ 。$

(8)

图 2 屈服面、边界面与硬化规律示意图

Figure 2. Schematic diagram of yield surface, bounding surface and hardening law

下载: 全尺寸图片幻灯片

表征当前应力点和像应力点关系。结合图 2另有

$\eta = \bar q/\bar p' = q/p' \ \ 。$

(9)

将式（8），（9）代入式（7）得

$\begin{aligned} \chi= & \frac{\beta^3 M_{\mathrm{c}}^2}{(\beta-1)^2 \eta^2+\beta^2 M_{\mathrm{c}}^2}+ \\ & \frac{\sqrt{(\beta-1)^2 \beta^2 M_{\mathrm{c}}^2\left((1-2 \beta) \eta^2+\beta^2 M_{\mathrm{c}}^2\right)}}{(\beta-1)^2 \eta^2+\beta^2 M_{\mathrm{c}}^2} \end{aligned}。$

(10)

对于硬化规律，本文仍选取塑性体应变 $\varepsilon _{\text{v}}^{\text{p}}$ 作为硬化参量，与修正剑桥模型类似，并参考Yin等^[29]考虑颗粒破碎对压缩特性的影响进行修正：

${\bar p'_0} = {\bar p'_{0i}}\exp \left( {\frac{{1 + {e_0}}}{{\lambda - \kappa }}\varepsilon _{\mathrm{v}}^{\mathrm{p}}} \right)\exp ( - \rho {B_{\mathrm{r}}})\ \ 。$

(11)

式中： $\rho$ 为破碎相关参数； ${\bar p'_{0i}}$ 为初始边界面尺寸，其确定方法可由压缩平面BCSL和卸载回弹线(SL，取与BCSL相同形式)来确定^[28]。

如所示，A( ${\bar p'_A},{e_A}$ )既在BCSL上，又在SL上，因此，分别代入BCSL和SL方程可得

$\begin{array}{l}{e}_{A}={e}_{N0}-\kappa \mathrm{ln}\left[({\overline{{p}^{\prime }}}_{A}+{p}_{\mathrm{b}})/{p}_{\mathrm{r}}\right]\text{，}\\ {e}_{A}=\Pi -\lambda \mathrm{ln}\left[({\overline{{p}^{\prime }}}_{A}+{p}_{\mathrm{b}})/{p}_{\mathrm{r}}\right]。\end{array}\}$

(12)

式中： ${e_{N0}}$ 为回弹线在压缩平面截距；同时，点( ${p_{i \mathrm{c}}^{\prime}},{e_0}$ )也在SL上，可表示为

$e_0=e_{N 0}-\kappa \ln \left[\left(p_{i \mathrm{c}}^{\prime}+p_{\mathrm{b}}\right) / p_{\mathrm{r}}\right] \ \ 。$

(13)

联立式（12），（13）以及边界面上顶点A( ${\bar p'_A},\bar q$ )与右端点( ${\bar p'_{0i}},0$ )的关系 ${\bar p'_A} = \beta {\bar p'_{0i}}$ ，可得

$\left. {\begin{array}{*{20}{l}} {\bar p_{0i}^\prime = 1/\beta .}\\ {{p_{\rm{r}}}\exp \left[ {\frac{{\mathit{\Pi } - {e_0} - \kappa \ln \left( {\frac{{p_{i{\rm{c}}}^\prime + {p_{\rm{b}}}}}{{{p_{\rm{r}}}}}} \right)}}{{(\lambda - \kappa )}}} \right] - {p_{{\rm{b}}^\circ }}} \end{array}} \right\}$

(14)

由式（11），（14）可知，在初始边界面的确定和边界面演化过程均考虑了颗粒破碎的影响。

2.3 状态相关非正交流动法则

根据传统的边界面塑性理论^[30]，加载方向可由边界面方程得到

$\begin{array}{l}\partial \overline{f}/\partial \overline{{p}^{\prime }}=2{\overline{{p}^{\prime }}}_{0}{M}_{\text{c}}^{\text{2}}{\beta }^{2}(\chi -\beta )\text{，}\\ \partial \overline{f}/\partial \overline{q}=2{\overline{{p}^{\prime }}}_{0}{(}^{1}\chi \eta 。\end{array}\}$

(15)

令

${L}_{f}^{1}=\sqrt{{(\partial \overline{f}/\partial \overline{{p}^{\prime }})}^{2}+{(\partial \overline{f}/\partial \overline{q})}^{2}}。$

(16)

为加载方向的模，则加载单位向量可表示为

$\begin{array}{l}{n}_\mathrm{p}=2{\overline{{p}^{\prime }}}_{0}{M}_\mathrm{c}^{2}{\beta }^{2}(\chi -\beta )/{L}_\mathrm{f}^{1}\text{，}\\ {n}_\mathrm{s}=2{\overline{{p}^{\prime }}}_{0}{(}^{1}\chi \eta /{L}_\mathrm{f}^{1}。\end{array}\}$

(17)

对于塑性流动方向，传统的边界面塑性理论中或采用相关联流动法则，使塑性流动方向和加载方向相同，但随着人们对土体认识和研究深入，发现砂土等粗粒料其塑性流动方向并不与塑性加载方向重合。如前所述，近年来，分数阶微分方法不断引起学者们的关注，采用分数阶微分方法，可以直接通过原边界面方程得到一个与加载方向不重合的分数阶塑性流动方向^[16-18]，即非正交流动法则。

由文献[16~18]可知，常用的两种分数阶微分定义中，Caputo型微分定义简单，具有弱奇异性，应用广泛，但因其对常数微分为0，对于椭圆方程，当微分阶次大于1时，具有不稳定性^[31]，因此，本文采用左Riemann-Liouville型分数阶微分，

${ }_a^{\mathrm{RL}} D_t^\alpha u(t)=\frac{1}{\mathit{\Gamma}(n-\alpha)} \frac{\mathrm{d}^n}{\mathrm{~d} t^n}\left(\int_a^t(t-\xi)^{n-\alpha-1} u(\xi) \mathrm{d} \xi\right) \ \ 。$

(18)

式中： $a < t$ ，为积分上下限( $a,t$ )； $n$ 为大于微分阶数 $\alpha$ 的最小整数； $\mathit{\Gamma} ( \bullet )$ 为Gamma函数，

$\mathit{\Gamma} (z) = {\int_0^\infty {\exp ( - \xi )\xi } ^{z - 1}}d\xi \ \ 。$

(19)

因此，根据Riemann-Liouville型分数阶微分定义及边界面方程，可以得到分数阶塑性流动方向为

$\begin{aligned} & { }_0^{\mathrm{RL}} D_{\bar{p}^{\prime}}^\alpha f\left(\bar{p}^{\prime}\right)=\frac{\partial^\alpha \bar{f}}{\partial \bar{p}^{\prime \alpha}}= \\ & \frac{2 M_{\mathrm{c}}^2 \beta^2 \chi^{1-\alpha} \bar{p}_0^{\prime 2-\alpha}[\chi-(2-\alpha) \beta]}{\mathit{\Gamma}(3-\alpha)}+ \\ & \frac{[2(1-\alpha)(2-\alpha)(1-\beta)] M^2 \beta^3 \chi^{-\alpha} \bar{p}_0^{\prime 2-\alpha}}{\mathit{\Gamma}(3-\alpha)}+ \\ & \frac{(1-\alpha)(2-\alpha)(1-\beta)^2 \eta^2 \chi^{2-\alpha} \bar{p}_0^{\prime 2-\alpha}}{\mathit{\Gamma}(3-\alpha)}\ \ , \end{aligned}$

(20)

$\begin{gathered} { }_0^{R L} D_{\bar{q}}^\alpha f(\bar{q})=\frac{\partial^\alpha \bar{f}}{\partial \bar{q}^\alpha}=M^2 \beta^2\left(\chi^2-2 \chi \beta+2 \beta-1\right) \cdot \\ \frac{(1-\alpha)(2-\alpha) \eta^{-\alpha} \chi^{-\alpha} \bar{p}_0^{\prime 2-\alpha}}{\mathit{\Gamma}(3-\alpha)}+ \\ \frac{2(1-\beta)^2 \eta^{2-\alpha} \chi^{2-\alpha} \bar{p}_0^{\prime 2-\alpha}}{\mathit{\Gamma}(3-\alpha)} \ \ 。 \end{gathered}$

(21)

令

$L_f^\alpha = \sqrt {{{({\partial ^\alpha }\bar f/\partial {{\bar p'}^\alpha })}^2} + {{({\partial ^\alpha }\bar f/\partial {{\bar q}^\alpha })}^2}} \ \ 。$

(22)

为塑性流动方向的模，则流动单位向量可表示为

$\begin{array}{l}{m}_{p}=({\partial }^{\alpha }\overline{f}/\partial {\overline{{p}^{\prime }}}^{\alpha })/{L}_{f}^{\alpha }\text{ }\text{，}\\ {m}_{s}=({\partial }^{\alpha }\overline{f}/\partial {\overline{q}}^{\alpha })/{L}_{f}^{\alpha }\text{ }。\end{array}\}$

(23)

根据土力学中剪胀比的定义可得

$\begin{aligned} & d=\frac{\mathrm{d} \varepsilon_{\mathrm{v}}^{\mathrm{p}}}{\mathrm{d} \varepsilon_{\mathrm{s}}^{\mathrm{p}}}=\frac{m_{\mathrm{p}}}{m_{\mathrm{s}}}=\left\{2 M_{\mathrm{c}}^2 \beta^2 \chi[\chi-(2-\alpha) \beta]+\right. \\ & \left.\left(2-3 \alpha+\alpha^2\right)\left[(2-\beta) M_{\mathrm{c}}^2 \beta^3+(1-\beta)^2 \eta^2 \chi^2\right]\right\} / \\ & \left\{M_{\mathrm{c}}^2 \beta^2\left(\chi^2-2 \chi \beta+2 \beta-1\right)\left(2-3 \alpha+\alpha^2\right) \eta^{-\alpha}+\right. \\ & \left.\quad 2(1-\beta)^2 \eta^{2-\alpha} \chi^2\right\} \ \ 。 \end{aligned}$

(24)

从式（24）可知，如果令分数阶阶数 $\alpha {\text{ = 1}}$ ，则

${d}_{\alpha =1}=\frac{{M}_\mathrm{c}{}^{2}{\beta }^{2}\text{(}\chi -\beta \text{)}}{{(1-\beta )}^{2}\chi \eta }=\frac{{n}_\mathrm{p}}{{n}_\mathrm{s}}\ \ 。$

(25)

分数阶剪胀方程和整数一阶微分求得的加载方向重合。表明，通过分数阶微分方法可以统一相关联和非相关联流动法则。

由于钙质砂等粗粒土对初始密实度和围压的相关性，本文采用Been等^[32]定义的状态参数 $\psi$ 并修正为 ${\psi _b}$ 来考虑钙质砂考虑破碎的状态相关特性：

${\psi }_\mathrm{b}=e-{e}_\mathrm{cs}\ \ 。$

(26)

式中： ${e_\mathrm{cs}}$ 为BCSL上的孔隙比。将状态参数引入塑性流动法则可以很好地考虑土体的状态相关行为^[33]：

$\alpha = \exp ({\psi _\mathrm{b}}) \ \ 。$

(27)

当 ${\psi _\mathrm{b}} < 0$ 时， $\alpha < 1$ ，土体主要表现为应变软化和体积膨胀行为；当 ${\psi _\mathrm{b}} > 0$ 时， $\alpha > 1$ ，土体主要表现为应变硬化和体积收缩行为。因此，状态相关流动法则能很好地描述土体剪切过程土体所处状态对塑性流动方向的影响，当达到临界状态时， $\eta ={M}_\mathrm{c}，$ ${\psi }_\mathrm{b}=0，\text{ }$ $\alpha = {\text{1}}$ ，此时流动方向与加载方向相同，且剪胀公式（24）自然满足临界状态 $d{\text{ = 0}}$ 的条件。

2.4 塑性模量

由边界面上的一致性条件：

$df = \frac{{\partial \bar f}}{{\partial \bar p'}}d\bar p' + \frac{{\partial \bar f}}{{\partial \bar q}}d\bar q + \frac{{\partial \bar f}}{{\partial {{\bar p'}_0}}}\frac{{\partial {{\bar p'}_0}}}{{\partial \varepsilon _\mathrm{v}^\mathrm{p}}}d\varepsilon _\mathrm{v}^\mathrm{p} = 0 。$

(28)

和塑性流动法则：

${\rm{d}}\mathit{\boldsymbol{\varepsilon }}_{ij}^{\rm{p}} = \mathit{\Lambda }\frac{{{\partial ^\alpha }\bar f}}{{\partial \mathit{\boldsymbol{\sigma }}_{ij}^\alpha }}\ \ 。$

(29)

式中： $\partial^\alpha \bar{f} / \partial \boldsymbol{\sigma}_{i j}^\alpha$ 为塑性流动方向； $\mathit{\Lambda }$ 为塑性乘子，其与边界面像应力塑性模量 ${\bar K_\mathrm{p}}$ 的关系为

$\mathit{\Lambda}=\frac{1}{\bar{K}_{\mathrm{p}}} \frac{\partial \bar{f}}{\partial \mathit{\boldsymbol{\sigma}}_{i j}} \mathrm{~d} \mathit{\boldsymbol{\sigma}}_{i j} \ \ 。$

(30)

将式（29），（30）代入式（28）得

${\bar K_\mathrm{p}} = - \frac{{\partial \bar f}}{{\partial {{\bar p'}_0}}}\frac{{\partial {{\bar p'}_0}}}{{\partial \varepsilon _\mathrm{v}^\mathrm{p}}}\frac{{{\partial ^\alpha }\bar f}}{{\partial {{\bar p'}^\alpha }}} \ \ 。$

(31)

由此，边界面上像应力归一化的塑性模量可表示为

${\bar H_\mathrm{b}} = {\bar K_\mathrm{p}}/(L_f^1L_f^\alpha ) \ \ \ 。$

(32)

根据边界面塑性理论，当前应力点落在边界面上时，真实应力塑性模量 $H = {\bar H_\mathrm{b}}$ ；达到临界状态时， $H = 0$ ；因此通过插值函数法，以应力平面坐标原点为映射中心，引入表征当前面与边界面距离的模量^{[28, 30]}：

${H_\mathrm{c}} = {h_0}p'\frac{{1 + {e_0}}}{{\lambda - \kappa }}\frac{\delta }{{{\delta _\mathrm{m}} - \delta }} \ \ \ 。$

(33)

式中： ${h_0}$ 为硬化模量参数； $\delta$ ， ${\delta _\mathrm{m}}$ 为像应力点到真实应力点和映射中心的距离。所以当前应力状态塑性模量可以表示为

$H = {\bar H_\mathrm{b}} + {H_\mathrm{c}} \ \ \ 。$

(34)

3. 模型验证

3.1 模型参数和确定方法

所建立的本构模型共有9个参数，其中主要包括颗粒破碎临界状态参数 $\Pi$ ， ${p_\mathrm{b}}$ ， $\lambda$ ， ${M_c}$ ，边界面形状参数 $\beta$ ，破碎硬化参数 $\rho$ ，塑性硬化参数 ${h_0}$ 和弹性参数 $\kappa$ ， $\nu$ 。其中临界状态参数和弹性参数同修正剑桥模型，将三轴试验结果绘制于压缩平面，通过拟合得到；边界面参数 $\beta$ 和破碎硬化参数 $\rho$ 可通过一组三轴数据分别绘制( $d - \eta$ )曲线^{[26, 28]}和( ${\varepsilon _\mathrm{v}} - {\varepsilon _\mathrm{a}}$ )曲线，利用最佳拟合方法得到。塑性模量参数 ${h_0}$ 可由( ${\varepsilon _\mathrm{s}} - \eta$ )曲线，通过最小二乘法拟合得到^[26]，或通过特征状态点（峰值点），塑性模量应为0的条件确定，即当加载到峰值状态时，令 $H = 0$ ，解出

${h_0} = - \frac{{\bar K_\mathrm{p}^\mathrm{peak}(\lambda - \kappa )(\delta _\mathrm{m}^\mathrm{peak} - {\delta ^\mathrm{peak}})}}{{{\delta ^\mathrm{peak}}(1 + {e_0})({{p'}^\mathrm{peak}}L_\mathrm{f}^1L_\mathrm{f}^\alpha )}} \ \ 。$

(35)

式中：上标“peak”表示峰值状态；该方法得到的 ${h_0}$ 对不同围压会有差异，需取最佳拟合结果。破碎应力 ${p_\mathrm{b}}$ 需通过试验结果的破碎指标值和压缩平面的BCSL拟合得到，建立相关关系。

3.2 模型预测结果和试验对比

为了验证模型对钙质砂颗粒破碎特性和状态相关特性模拟的适用性，选取了南海某吹填岛礁钙质砂^{[3, 34]}的三轴排水试验结果。Wu等^[34]针对3种不同初始密实度( ${D_\mathrm{r}}$ =30%， ${D_\mathrm{r}}$ =60%， ${D_\mathrm{r}}$ =90%)钙质砂在围压100，200，400，800，1600 kPa下分别进行三轴试验，本文将选取 ${D_\mathrm{r}}$ =30%的试验结果作为模型参数标定数据，以 ${D_\mathrm{r}}$ =60%和 ${D_\mathrm{r}}$ =90%试验数据验证建立的模型的适用性。

吴杨等^[3]在试验前后测得的颗分曲线，获得每组试验结束时的破碎量，并给出了颗粒破碎和塑性功的关系：

${B_\mathrm{r}} = \frac{{{W_\mathrm{p}}}}{{5518.85 + {W_\mathrm{p}}}} \ \ 。$

(36)

其中塑性功 ${W_\mathrm{p}}$ 可为

${W_\mathrm{p}} = \int {p'd\varepsilon _\mathrm{v}^\mathrm{p}} + \int {qd\varepsilon _\mathrm{s}^\mathrm{p}} \ \ 。$

(37)

如图 3，考虑颗粒破碎的影响，给出了相对密实度为30%^[34]时颗粒破碎临界状态线的演化：可以看出，当初始状态点在BCSL下方时，试样表现出软化型特征，当试样初始状态在BCSL上方时，表现出硬化型特征，与状态相关理论吻合。

图 3 压缩平面颗粒破碎临界状态线的演化

Figure 3. Evolution of BCSL in compression plane

下载: 全尺寸图片幻灯片

如，结合中BCSL的演化，按指数函数拟合了破碎应力 ${p_\mathrm{b}}$ 和颗粒破碎指数 ${B_\mathrm{r}}$ 的关系：

图 4 破碎应力与颗粒破碎指数拟合关系

Figure 4. Relationship between breakage stress and index

下载: 全尺寸图片幻灯片

${p_\mathrm{b}} = \mathit{\mathbb{Z}}\left[ {1 - \exp ( - \vartheta {B_\mathrm{r}})} \right] \ \ 。$

(38)

式中： $\mathit{\mathbb{Z} }$ ， $\vartheta$ 为拟合系数，当颗粒破碎指数为0时，破碎应力为0，BCSL退化为CSL，且由于破碎指数的极限值为1，破碎应力存在极限值，与粗粒土颗粒破碎不会全部碎到最小粒径，而是达到一种极限分形状态相一致。

在偏应力-平均有效应力平面，对于易破碎粗粒土，有学者建议采用非线性CSL^[35-36]， ${M_\mathrm{c}}$ 的值与围压、密实度或破碎相关；也有学者研究发现，砂颗粒破碎不会导致 $q$ - $p'$ 平面CSL变化^[9]。Wu等^[34]在研究中给出了 $q$ - $p'$ 平面的线性CSL，且与密实度相关，如图 5。

图 5

$q$ -

$p'$ 平面临界状态线^[34]

Figure 5. CSLs in

$q$ -

$p'$ plane^[34]

下载: 全尺寸图片幻灯片

根据相对密实度为30%的试验数据，结合3.1节给出的模型参数确定方法，对所建立模型的参数进行标定，得到模型参数如所示。需要指出的是，其中临界状态参数 $\Pi$ 与相对密实度相关^{[12, 26, 34, 37]}，临界状态应力比 ${M_\mathrm{c}}$ 按照原参考文献也取与密实度相关。

表 1 模型参数

Table 1. Model parameters

弹性参数	颗粒破碎临界状态参数	边界面形状参数	硬化参数
$\kappa = 0.005$ $\nu = 0.3$	$\mathit{\Pi} = 2.4$ ^，2.332^，2.285^** $\lambda = 0.19$ ${M_\mathrm{c}} = 1.54$ ^, 1.592^，1.68^** $\mathbb{Z}$ =671 $\vartheta$ =5.88	$\beta = 0.58$	$\rho = 0.1$ ${h_0} = 1.2$
注：上标“^，^，^**”分别代表相对密实度30%，60%，90%的取值。

下载: 导出CSV

| 显示表格

为探究所建立考虑钙质砂颗粒破碎和状态相关的分数阶本构模型在模拟应力应变关系中颗粒破碎的影响，选取相对密实度为30%时的模型参数，见，选取了低围压200 kPa和高围压1600 kPa两种情况下，分别研究模型考虑颗粒破碎（即令破碎应力 ${p_\mathrm{b}}$ 与破碎指数满足式（38）关系，在计算中破碎指数不断更新演化）与不考虑颗粒破碎（即在计算中直接令破碎指数为0）时的应力应变曲线的区别。如图 6可以看出，模型可以模拟材料在高围压下的应变硬化和体积收缩行为，以及低围压下材料的应变软化和体积膨胀行为。同时，无论对于硬化型还是软化型材料，当考虑颗粒破碎影响时，所预测的偏应力均小于不考虑颗粒破碎时的结果，而体应变均大于不考虑颗粒破碎时的结果，与颗粒破碎的发生会导致材料强度降低及其剪胀性降低且增加压缩性的试验现象相吻合^{[3, 5-6]}，且在低围压时，两种预测结果相差较小，其主要原因是低围压下破碎量较小。另一方面，该对比结果表明，工程应用中考虑颗粒破碎影响的重要性：如对于地基承载力和变形的预测中，忽略颗粒破碎的影响所计算得到的强度高，变形小，设计或施工中将高估其强度，而低估了其变形，这意味着在工程设计和施工中应关注破碎对建筑物稳定性和安全性的影响。

图 6 模型参数影响分析

Figure 6. Analysis of model parameters on simulations

下载: 全尺寸图片幻灯片

同时，在考虑颗粒破碎的基础上，探究了破碎硬化参数 $\rho$ 对材料破碎硬化过程的影响。如中，对比破碎硬化参数 $\rho = 1$ 和 $\rho = 10$ 结果可知，破碎硬化参数增大，会导致强度降低，体变（压缩性）增大，尤其是在高围压下，该现象更加明显，主要是由于高围压下颗粒破碎量更大，对材料的变形和强度影响更大，说明对颗粒破碎影响材料硬化规律的修正能够合理地反映颗粒破碎对材料的应力应变响应的影响。

图 7~9给出了钙质砂分别在初始相对密实度为30%，60%，90%时，试样在不同围压下偏应力-轴向应变、体积应变-轴向应变及孔隙比-平均有效应力的关系和模型模拟及预测结果的对比。从中可以看出，

图 7 模型模拟结果和试验结果^[33] (

${D_r}$ =30%)对比

Figure 7. Comparison between simulations and test data^[33]with

${D_r}$ =30%

下载: 全尺寸图片幻灯片

图 8 模型模拟结果和试验结果^[33] (

${D_r}$ =60%)对比

Figure 8. Comparison between simulations and test data^[33]with

${D_r}$ =60%

下载: 全尺寸图片幻灯片

图 9 模型模拟结果和试验结果^[33] (

${D_r}$ =90%)对比

Figure 9. Comparison between simulations and test data^[33]with

${D_r}$ =90%

下载: 全尺寸图片幻灯片

对于3种密实度的钙质砂在较低围压下如100 kPa，均表现出明显的应变软化和体积膨胀行为，随着围压的增加，土体也从庆变软化行为逐渐转为应变硬化行为。图 7为密实度30%时试验结果与模拟结果的对比，该组曲线主要用来标定表 1中部分模型参数。从模拟结果可以看出，建立的模型能够反映钙质砂在不同围压条件下的应力应变响应，能合理表征钙质砂在破碎影响下随围压增加逐渐由应变软化（体积膨胀）转为应变硬化（体积收缩）行为。利用表 1中标定的参数，对相对密实度为60%和90%的钙质砂试样进行了模拟，如，所示，模拟结果在 $q - {\varepsilon _\mathrm{a}}$ ， ${\varepsilon _\mathrm{v}} - {\varepsilon _\mathrm{a}}$ 和 $e$ - $\ln p'$ 平面均与试验结果吻合较好，合理地反映了不同围压和密实度条件，钙质砂在颗粒破碎影响下的力学响应。总体上，所有模拟结果误差均在可接受范围内，验证了所建立的本构模型能够考虑钙质砂的颗粒破碎行为和状态相关行为。

4. 结论

本文基于临界状态理论，引入破碎应力建立一系列破碎临界状态线；基于分数阶微分定义和状态相关理论，建立了状态相关分数阶塑性流动法则，并基于边界面塑性理论框架，推导了考虑颗粒破碎的硬化规律，建立了考虑钙质砂颗粒破碎和状态相关的非正交流动边界面本构模型，并通过一组钙质砂三轴试验数据标定参数，并对比另两组试验数据和模拟结果，得出4点结论。

（1）钙质砂颗粒破碎影响临界状态线的位置，因此，引入破碎应力建立了颗粒破碎临界状态线，能够很好地描述颗粒破碎导致压缩平面临界状态线位置发生向下移动的现象。

（2）基于颗粒破碎临界状态线和边界面塑性硬化，推导了考虑颗粒破碎影响的边界面尺寸演化规律和硬化规律，其能够反映钙质砂颗粒破碎导致的强度和剪胀性降低，压缩性增加现象。

（3）结合状态相关理论和分数阶微分定义建立了状态相关的塑性流动法则，其能够统一地描述试样在剪切过程中状态变化对塑性流动方向的影响，并能够统一相关联和非相关联流动法则。

（4）基于边界面塑性理论、状态相关理论和临界状态理论建立的分数阶边界面模型，能够很好地模拟钙质砂考虑颗粒破碎影响和不同初始密实度和围压下的应变软化和应变硬化特性。

图 1 ZLB-1型三联流变直剪仪

Figure 1. ZLB-1 triple rheological direct shear apparatus

下载: 全尺寸图片幻灯片

图 2 不同剪切荷载下的应变-时间关系的蠕变曲线

Figure 2. Creep curves of strain-time relationship under different shear loads

下载: 全尺寸图片幻灯片

图 3 各级荷载作用下的蠕变应变增量

Figure 3. Creep strain increments under various shear loads

下载: 全尺寸图片幻灯片

图 4 剪应力-剪应变的等时曲线

Figure 4. Isochronous curves of shear stress-shear strain

下载: 全尺寸图片幻灯片

图 5 剪应变-时间的双对数曲线

Figure 5. Double-logarithmic graph of shear strain-time

下载: 全尺寸图片幻灯片

图 6 不同参考时间所对应的误差百分比

Figure 6. Error percentage corresponding to different time

下载: 全尺寸图片幻灯片

图 7 修正的Singh-Mitchell模型预测效果

Figure 7. Fitting effect of modified Singh-Mitchell model

下载: 全尺寸图片幻灯片

表 1 土样的基本物理性质

Table 1 Basic physical properties of soil

含水率w/% 重度/(kN·m^-3) 液限/% 塑限/% 孔隙比e

48 17.2 48.6 31.6 1.28

下载: 导出CSV

表 2 土样施加水平剪应力方案

Table 2 Scheme of applying horizontal shear stress to soil samples

固结压力/kPa 水平剪应力大小与级数/kPa
1 2 3 4 5

400 32 71 103 134 155

下载: 导出CSV

表 3 不同剪应力的 $\ln ε$ -lnt的拟合直线斜率n

Table 3 Fitting straight line slopes n of $\ln ε$ -lnt for different shear stresses

剪应力/kPa 32 71 103 134 155

n值 0.0098 0.0135 0.0134 0.0112 0.0139
n平均值 0.0124

下载: 导出CSV

表 4 不同的参考时间下模型参数b和C值

Table 4 Values of model parameters b and C at corresponding time

参考时间t₁/h b（直线斜率） C（直线截距）

0.25 0.514 2.040
0.5 0.515 2.054
1 0.519 2.063
2 0.524 2.067
3 0.526 2.072
6 0.531 2.078
9 0.532 2.092
12 0.538 2.107
24 0.531 2.119
36 0.528 2.137
48 0.526 2.147

下载: 导出CSV

参考文献(7)

[1]	蒙韵, 周树. 深圳地区软土的工程特性及软基处理措施探讨[J]. 资源环境与工程, 2016, 30(3): 450-452. doi: 10.16536/j.cnki.issn.1671-1211.2016.03.047 MENG Yun, ZHOU Shu. Engineering characteristics of soft soil in Shenzhen area and discussion on soft foundation treatment measures[J]. Resources Environment and Engineering, 2016, 30(3): 450-452. (in Chinese) doi: 10.16536/j.cnki.issn.1671-1211.2016.03.047
[2]	李妤妤. 南沙软土宏微观蠕变试验及本构模型研究[D]. 广州: 广州大学, 2020. LI Yu-yu. Macro and Micro Creep Test and Constitutive Model Research of Nansha Soft Soil[D]. Guangzhou: Guangzhou University, 2020. (in Chinese)
[3]	朱铨雯. 珠海软土宏微观蠕变试验及非线性模型研究[D]. 广州: 广州大学, 2020. ZHU Quan-wen. Research on Macro and Micro Creep Test and Nonlinear Model of Zhuhai Soft Soil[D]. Guangzhou: Guangzhou University, 2020. (in Chinese)
[4]	肖谊. 广州黄埔地区软土蠕变试验及蠕变模型选择[J]. 现代城市轨道交通, 2017(9): 28-31. https://www.cnki.com.cn/Article/CJFDTOTAL-XDGD201709008.htm XIAO Yi. Creep test of soft soil in Huangpu area of Guangzhou and selection of creep model[J]. Modern Urban Rail Transit, 2017(9): 28-31. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-XDGD201709008.htm
[5]	刘柱辉. 软土直剪蠕变试验分析与本构模型研究[D]. 广州: 广州大学, 2018. LIU Zhu-hui. Direct Shear Creep Test Analysis and Constitutive Model Study of Soft Soil[D]. Guangzhou: Guangzhou University, 2018. (in Chinese)
[6]	SINGH A, MITCHELL J K. General stress-strain-time function for soils[J]. Journal of the Clay Mechanics and Foundation Division, 1968, 94(SM1): 21-46.
[7]	王常明, 王清, 张淑华. 滨海软土蠕变特性及蠕变模型[J]. 岩石力学与工程学报, 2004(2): 227-230. doi: 10.3321/j.issn:1000-6915.2004.02.010 WANG Chang-ming, WANG Qing, ZHANG Shu-hua. Creep characteristics and creep model of coastal soft soils[J]. Chinese Journal of Rock Mechanics and Engineering, 2004(2): 227-230. (in Chinese) doi: 10.3321/j.issn:1000-6915.2004.02.010

施引文献(6)

期刊类型引用(6)

1.	蒋承辉. 基于BOTDA的光纤位移传感器在建筑基坑沉降监测中的实验研究. 企业科技与发展. 2025(01): 98-101 . 百度学术
2.	吴哲辉，何宁，姜彦彬，孔洋，槐伟. 堆石坝内部二维变形一体化分布式监测技术试验研究. 水利信息化. 2025(02): 56-63 . 百度学术
3.	管翰林，蒋陵，王驭扬，张灿. 基于DOFS技术的电力光缆多维度故障预警研究. 电力电子技术. 2025(05): 78-81+105 . 百度学术
4.	赵建勋，高冰，高丽娟. 分布式光纤传感技术下架空输变电线路监测. 电子设计工程. 2025(11): 72-76 . 百度学术
5.	黄惇汉. 基于分布式光纤传感技术的车辆撞击道路护栏检测研究. 运输经理世界. 2024(22): 122-124 . 百度学术
6.	马刚，艾志涛，郭承乾，李少林，陈华，周伟. 高土石坝变形监测研究进展. 水利学报. 2024(10): 1174-1186 . 百度学术