0.   引言

天然土体颗粒间通常具有胶结物质,比如碳酸钙、微生物、黏粒、甲烷水合物等,称为结构性土。结构性土可看成是一种特殊的胶结颗粒材料,颗粒间胶结物质的存在使结构性土的宏观力学性质不同于重塑土^[1-2]。为了分析结构性土工程问题,需要建立合理的结构性土本构模型。在临界状态土力学基础上,相关学者引入状态变量,分别建立了边界面模型^[3]、下负荷面模型^[4]、统一硬化模型^[5]等重塑土本构模型。在加荷过程中,结构性土颗粒间的胶结逐渐破坏并退出承担应力,土体向重塑土转化。因此,在重塑土本构模型基础上,通过引入表征结构性的参量及其演化规律建立结构性土本构模型是一条可行的研究途径。不少学者提出了结构性土本构模型,在模型中假定了不同的结构性表征参量。比如：①通过比较结构性土屈服面大小${\overline{p}}_{\text{c}}$和重塑土屈服面大小$p_{\text{c}}$的不同,定义$\Delta p_{\text{c}}={\overline{p}}_{\text{c}}-p_{\text{c}}$^[6]或$R^{*}=p_{\text{c}}/$${\overline{p}}_{\text{c}}$^[7]为结构性表征参量；②通过引入胶结组分（相对完整状态）和摩擦组分（完全调整状态）应力分担参量（扰动参量、结构破损参量）D^[8]或B₀^[9-10]表达结构性衰退程度。由于胶结破损演化难以通过试验测得,结构性表征参量的演化通常需要假设,缺乏微观试验结果验证。

离散单元法^[11]模拟结构性土弥补了室内微观测试在获取胶结破损演化方面的不足,已广泛应用于结构性土体宏微观力学性质研究。Jiang等^[10]通过开展结构性土的二维离散元模拟,对3种结构性土本构模型结构性参量演化假设的合理性进行了验证,基于离散元结果提出了新的胶结破损演化式,并应用到本构模型研究^[12],不足之处在于胶结破损演化式中的参数受加载路径影响。通过PFC软件植入完整胶结接触模型（考虑了颗粒间抗转动和抗扭转能力及胶结尺寸对刚度和强度的影响）建立了结构性土三维离散元试样,开展了侧限压缩、等向压缩、等应力比压缩以及常规三轴和真三轴试验的离散元数值模拟,分析了结构破损参量B₀随等效塑性应变的演化规律,并提出了应力路径相关性低的结构破损参量B_σ。

1.   胶结接触模型

离散元模拟采用软硬复合胶结接触模型框架^[13],考虑了颗粒间抗转动和抗扭转能力^[14]及胶结尺寸对刚度和强度的影响^[15]。对于通过胶结间接传力/力矩的接触（“厚胶结”）直接采用胶结接触公式计算接触力/力矩；对于颗粒直接接触并有胶结传力/力矩的接触（“薄胶结”）,按并联模式通过颗粒接触力/力矩和胶结接触力/力矩之和计算接触力/力矩。

1.1   颗粒接触

假设颗粒间通过一圆形截面接触,半径为

式中,$\beta$为抗转动系数,$R\text{=}2R_1{R}_2/(R_1+R_2)$,$R_1$,$R_2$为两颗粒半径。

接触间力/力矩的传递通过下列公式计算：

式中 $F_{\text{a}}$为粒间引力,比如毛细力、范德华力、电子引力等；$F_{\text{n}}^{\text{l}}$,$F_{\text{s}}^{\text{l}}$,$M_{\text{r}}^{\text{l}}$,$M_{\text{t}}^{\text{l}}$分别为颗粒法向力、切向力、弯矩和扭矩；$k_{\text{n}}$,$k_{\text{s}}$,k_r,k_t分别为颗粒法向、切向、抗弯、抗扭刚度,$k_{\text{r}}=0.25k_{\text{n}}{R}_{\text{c}}{}^2$,$k_{\text{t}}=0.5k_{\text{s}}{R}_{\text{c}}{}^2$；${\delta }_{\text{n}}$为颗粒重叠量,$\Delta {\delta }_{\text{s}}$,$\Delta {\theta }_{\text{r}}$和$\Delta {\theta }_{\text{t}}$分别为剪切位移、转动角和扭转角增量。

颗粒接触的法切向刚度可按PFC软件手册^[16]按下式计算：

式中,$A=\text{π}{R}^2$,$E^{*}$是颗粒接触的等效模量,${\kappa }^{*}$是颗粒接触的法切向刚度比,$L$为两颗粒半径之和。

接触抗剪强度、抗弯强度和抗扭强度分别为^[14-15]

式中,$\mu$是摩擦系数,${\zeta }_{\text{c}}$取2.1。颗粒接触力和力矩达到以上强度后颗粒将进入滑移、转动和扭转状态。

1.2   胶结计算

假定胶结物质呈圆形截面连接两颗粒,则胶结法向力$F_{\text{n}}^{\text{b}}$,切向力$F_{\text{s}}^{\text{b}}$,弯矩$M_{\text{r}}^{\text{b}}$和扭矩$M_{\text{t}}^{\text{b}}$分别为^[16]

式中 $k_{\text{n}}^{\text{b}}$,$k_{\text{s}}^{\text{b}}$,$k_{\text{r}}^{\text{b}}$,$k_{\text{t}}^{\text{b}}$分别为胶结法向、切向、抗弯和抗扭刚度,$k_{\text{t}}^{\text{b}}=0.5k_{\text{s}}^{\text{b}}{R}_{\text{b}}{}^2$,$k_{\text{r}}^{\text{b}}=0.25k_{\text{n}}^{\text{b}}{R}_{\text{b}}{}^2$；$\Delta {\delta }_{\text{n}}^{}$是法向位移增量；$R_{\text{b}}={\lambda }_{\text{b}}R$为胶结半径,${\lambda }_{\text{b}}$为胶结半径系数。

为了考虑胶结厚度和宽度对胶结刚度的影响,用下式计算胶结法向刚度（图1）^[15]：

图  1  胶结刚度计算示意图

Figure  1.  Schematic diagram for calculation of bond stiffness

下载: 全尺寸图片 幻灯片

式中,${\delta }_{A1}$,${\delta }_{B1}$,${\delta }_{C1}$分别为颗粒1的A1_,B1和C1三段的线柔度,

式中,$E_{\text{1}}$,$E_{\text{2}}$分别为骨架材料和胶结材料的等效变形模量,$R_1$为颗粒1的半径,$h_0^{\text{b}}$为胶结厚度,同理可计算${\delta }_{A2}$,${\delta }_{B2}$。在模型使用中,$E_{\text{1}}={\overline{E}}^{*}$为胶结等效模量,$E_{\text{2}}={\eta }_{\text{E}}{\overline{E}}^{*}{}_{}$,${\eta }_{\text{E}}$为模量折减系数。${\xi }_{\text{E}}$,${\Lambda }_{\text{E}}$计算如下：

胶结切向刚度根据法切向刚度比${\kappa }_{\text{b}}$计算：

采用考虑胶结厚度和宽度的胶结强度计算方法^[15,17]。胶结抗压和抗拉强度计算如下：

式中,${\sigma }_{\text{c}}^{\text{b}}$,${\sigma }_{\text{t}}^{\text{b}}$分别为胶结材料抗压和抗拉强度；$c_1=1.45$,$c_2=-0.7$和$c_3=1/6$是反映胶结半径和厚度对抗压强度影响的参数（定值）。

在一定的法向力下（未达到胶结抗压和抗拉强度）,胶结可能承受剪切力、弯矩和扭矩,复合荷载作用下胶结强度符合椭球包面：

式中,$R_{\text{s}}^{\text{b}}$,$R_{\text{r}}^{\text{b}}$,$R_{\text{t}}^{\text{b}}$分别为胶结抗剪、抗弯和抗扭强度,

式中,$S_{\text{s}}^{\text{b}}$,$S_{\text{r}}^{\text{b}}$,$S_{\text{t}}^{\text{b}}$为胶结抗剪、抗弯和抗扭包面的形状,

式中,下标i是s,r和t,$m_{\text{s}}=0.5$,$m_{\text{r}}=0.3$和$m_{\text{t}}=0.36$为参数（定值）,$f_{\text{n}}^{\text{b}}=(F_{\text{n}}^{\text{b}}+R_{\text{nt}}^{\text{b}})/(R_{\text{nc}}^{\text{b}}+R_{\text{nt}}^{\text{b}})$为正则化胶结法向力。

胶结的存在对结构性土体的宏微观力学性质有很大影响。胶结达到强度包面后会发生胶结破坏,胶结破坏后胶结部分不再承担力和力矩,接触力和力矩退化到残余值（颗粒接触部分）。

2.   离散元模拟步骤

2.1   离散元试样制备

以饱和结构性黄土^[15]为基准制备离散元试样和选取接触参数,为提高研究的适用性开展了不同胶结材料强度离散元试样的数值分析,下文仅以一种胶结材料强度为例进行分析。采用分层欠压法^[18]分5层制备离散元试样,颗粒数目42180个,如图2（a）所示。结构性土离散元试样制备采用的颗粒级配曲线如图2（b）所示。试样生成后在12.5 kPa竖向压力下压缩稳定以反映原位K₀应力状态。试样稳定后若两个颗粒之间的间距小于$g_{\text{c}}{R}_{\text{b}}$则生成胶结,$g_{\text{c}}$为胶结临界厚度系数。

图  2  离散元试样及其粒径级配

Figure  2.  DEM sample and grain-size distribution curve

下载: 全尺寸图片 幻灯片

颗粒间引力作用能使试样以更松散状态稳定赋存^[19],试样制备过程施加弱吸引力以制备松散试样^[15]：

式中,d₅₀为中值粒径,${\sigma }_{\text{a}}$为弱吸引力系数,根据典型黄土的最大稳定孔隙比和饱和重塑土黏聚力,通过试算取4 kPa。颗粒与墙体间仅承担法向力作用,忽略摩擦。颗粒密度=2710 kg/m³,颗粒局部阻尼系数=0.7。颗粒之间的接触参数：颗粒等效模量=800 MPa,法切向刚度比=1.5,摩擦系数=0.5,抗转动系数=0.21,胶结等效模量=200 MPa,胶结模量折减系数=0.2,胶结法切向刚度比=2,胶结半径系数=0.35,胶结材料抗压强度=25MPa,胶结材料抗拉强度=2.5MPa,胶结临界厚度系数=0.1。

2.2   不同应力路径离散元模拟步骤

为了研究不同应力路径结构性土的胶结破损演化,对上述结构性土离散元试样进行了分级加载的侧限压缩试验、等向压缩试验和等应力比压缩试验,以及应变加载的常规三轴试验和等p真三轴试验。

对分级加载压缩试验,上下墙体从12.5 kPa开始分级施加荷载,每次施加的荷载约为上级荷载的$\sqrt{2}$倍；侧墙伺服应力目标通过应力比要求计算（侧限压缩试验使4个侧墙固定）。

对三轴试验,首先将离散元试样等向压缩到不同的围压,然后对侧墙进行伺服控制,使上下两墙相向运动进行剪切试验。在常规三轴试验中,侧墙的伺服目标为围压。在真三轴试验中通过平均应力和中主应力系数b不变的条件计算侧墙伺服目标应力。

3.   离散元模拟结果

为验证离散元数值分析结果的有效性,在研究结构性土的胶结破损演化规律之前,需确保离散元试样能够反映结构性土的主要压缩和剪切特性。考虑到本文针对一般结构性土体（不同胶结强度的胶结颗粒材料）进行破损规律数学建模,因此,下面主要进行定性对比,未对某一特定的结构性土进行标定。

图3给出了结构性土离散元试样在不同压缩路径下的压缩曲线。图3中$\eta =q/p$为应力比,$\text{CSR}={\sigma }_3/{\sigma }_1$为主应力比。随着平均应力的增加,离散元试样压缩线先缓慢下降而后逐渐转变为快速下降,存在结构屈服应力,离散元试样能够定性再现结构性土体的主要压缩特性^[20]。离散元试样体积变形特征与真实结构性土定量上差别较大,主要是因为颗粒形状、颗粒破碎机制等微观结构的差别。

图  3  结构性土压缩性质

Figure  3.  Compression properties of structured soils

下载: 全尺寸图片 幻灯片

图4给出了结构性土离散元试样常规三轴试验应力–应变关系和孔隙比–应变关系。离散元分析能够定性再现结构性土体三轴剪切性质^[2,21]。试样在剪切初始段模量很大,近似弹性段。随着围压的增加,试样抗剪强度增加,孔隙比减小,试样从倾向于软化和剪胀向倾向于硬化和剪缩发展。离散元反映结构性土剪切破坏特性的效果优于反映压缩变形特性。

图  4  结构性土常规三轴剪切性质

Figure  4.  Conventional triaxial test results of structured soils

下载: 全尺寸图片 幻灯片

图5给出了结构性土离散元试样真三轴试验应 力-应变关系和孔隙比–应变关系（以p=300 kPa为例）。可见离散元模拟结果能够反映室内试验中中主应力系数b（$b=({\sigma }_2-{\sigma }_3)/({\sigma }_1-{\sigma }_3)$）的影响^[22-23]。同样平均应力下随着b的增加,破坏偏应力下降。b=0时试样倾向于剪缩,随着b的增加试样剪缩程度减弱。天然土试样的变形模量相对离散元试样和人工制备结构性土试样较低,这可能与试样扰动有关。

图  5  结构性土真三轴剪切性质

Figure  5.  True triaxial test results of structured soils

下载: 全尺寸图片 幻灯片

4.   胶结破损演化规律

4.1   胶结颗粒材料微观力学

土体可以看成一种颗粒集合体,重塑土（无黏结）可以看成是一种散粒体材料而结构性土可以看成是一种胶结颗粒材料。对于理想的散粒体材料,土中的应力通过接触力进行传递,而土中的应变由颗粒间的相对变形产生。基于此,学者通过研究土颗粒之间的接触关系和微观变形机理,建立了理想散粒体材料宏观应力和应变与微观接触力和变形之间的基本方程^[24-25]。Jiang等^[10]在此基础上建立了胶结颗粒材料的宏微观理论框架。

对于颗粒材料,荷载通过颗粒之间的接触力进行传递,因此颗粒材料的代表性单元的平均应力张量可根据下式求得

式中,V为代表性单元的体积,N为代表性单元内颗粒间接触点的总数,$F_i^k$为第k个接触点处接触力的分量,$l_j^k$为第k个接触点处方向向量的分量。

对于结构性土体,根据颗粒间接触点是否存在胶结,土体代表性单元平均应力张量${\overline{\sigma }}_{ij}^{}$可分解为胶结接触分担应力${\overline{\sigma }}_{ij}^{\text{in}}$与无胶结接触分担应力${\overline{\sigma }}_{ij}^{\text{d}}$之和,

式中,$N^{\text{in}}$,$N^{\text{d}}$分别为代表性单元内胶结接触点和无胶结接触点的数目,且$N=N^{\text{in}}+N^{\text{d}}$。$V^{\text{in}}$,$V^{\text{d}}$分别为代表性单元内胶结颗粒所占体积与无胶结颗粒所占体积（含相应的孔隙体积）,且$V=V^{\text{in}}+V^{\text{d}}$；$\lambda =V^{\text{d}}/V$为体积破损率,即代表性单元内无胶结组分所占体积与代表性单元总体积的比值；${\overline{\sigma }}_{ij}^{\text{in}}$和${\overline{\sigma }}_{ij}^{\text{d}}$可由下式计算：

为了研究结构性土体的结构性演化,需要在本构模型中引入（或蕴含）胶结破损参量。体积破损率就是一个简单直观的胶结破损参量,考虑到胶结组分和无胶结组分的体积难以划分,可用无胶结接触数量与总接触数量比值近似表示,

无胶结组分分担的应力与总应力的比值也可以作为胶结破损参量,

胶结破损参量表征胶结破损程度或结构性衰退程度。当胶结破损参量为0时,土体结构性完好,土体应力全部由胶结组分承担；随着胶结破损参量的增加,土体应力由胶结组分和无胶结组分共同分担（胶结组分承担应力逐渐变小）；当胶结破损参量为1时,土体结构性完全破坏,土体应力全部由无胶结组分承担。

胶结破损参量$\lambda$,$B_0$难以通过室内试验获得,离散单元法可以用来研究$\lambda$,$B_0$的演化。尽管$\lambda$,$B_0$参量的大小有一定不同,但二者演化规律是相似的,因此本文以破损参量$B_0$为例进行研究。

4.2   胶结破损演化规律

随着塑性应变的发展,胶结逐渐破损。由于弹性应变较小,塑性应变采用如下简单方法计算。离散元试样的弹性能$E_{\text{k}}$为颗粒弹性能$E_{\text{k}}^{\text{p}}$和胶结弹性能$E_{\text{k}}^{\text{b}}$之和：

式中,$V_{\text{DEM}}$为离散元试样的体积,${\epsilon }_{\text{v}}^{\text{e}}$为弹性体应变,${\epsilon }_{\text{s}}^{\text{e}}$为弹性偏应变,假设体积模量K和剪切模量G相等,联立以上各式可求得${\epsilon }_{\text{v}}^{\text{e}}$,${\epsilon }_{\text{s}}^{\text{e}}$,从而求出塑性体应变和塑性偏应变。

图6给出了破损参量随塑性体应变和塑性偏应变的变化规律。对于三轴剪切试验,塑性体应变较小且变化复杂,而对于等应力比压缩试验（应力比较小时）,塑性偏应变较小,二者均不适合单独用于描述试样破损参量演化。因此,胶结的破损既与塑性体应变有关,又与塑性偏应变有关。

图  6  不同应力路径下破损参量B₀随塑性体应变和塑性偏应变的演化规律

Figure  6.  Evolution of degradation variable B₀ versus plastic volumetric and deviator strain for DEM sample under different stress paths

下载: 全尺寸图片 幻灯片

为了考虑两种塑性应变对胶结破损参量的影响,采用下式简单等效塑性应变研究胶结破损演化。

式中,${\epsilon }_{\text{v}}^{\text{p}}$为塑性体应变,${\epsilon }_{\text{s}}^{\text{p}}$为塑性偏应变。

图7给出结构性土离散元试样胶结破损参量$B_0$随等效塑性应变的变化规律。从图7中可见,对不同应力路径,破损参量$B_0$随等效塑性应变的增加先较快增加,而后缓慢增加。但不同应力路径的破损参量$B_0$演化在数值上有较大不同。对不同应力比压缩试验,随着应力比的增加,$B_0$表征的胶结破损速率加快；对不同中主应力系数b的真三轴试验,随着b的增加,$B_0$表征的胶结破损速率减小；对不同围压常规三轴试验,围压也对$B_0$表征的胶结破损速率影响显著。

图  7  不同应力路径下破损参量B₀随等效塑性应变的演化规律

Figure  7.  Evolution of degradation variable B₀ versus equivalent plastic strain for DEM sample under different stress paths

下载: 全尺寸图片 幻灯片

考虑到胶结破损跟胶结强度相关,在宏观上,胶结破损速率跟目前应力状态和强度包线之间的相对关系相关,结构参数$B_0$破损的速率随应力比的增加而增加。因此,将应力比和强度包线引入胶结破损参量$B_0$,通过应力比与峰值应力比的比值$\eta /M_{\text{f}}$调整胶结破损速率,构建一个新的胶结破损参量如下式,从而蕴含应力路径对胶结破损的影响,

式中,$M_{\text{f}}=q_{\text{f}}/p$为峰值应力比,本文通过常规三轴试验结果采用莫尔–库仑准则得到峰值偏应力$q_{\text{f}}=$$0.9275p+123.8\text{ kPa}$。

图8给出了结构性土离散元试样胶结破损参量$B_{\sigma }$随等效塑性应变的破损规律。对不同应力路径,胶结破损参量$B_{\sigma }$随等效塑性应变逐渐增加,但破损速率逐渐降低。除等向压缩试验（$\eta =0$）外,试验应力比、围压和中主应力系数对$B_{\sigma }$表征的胶结破损演化影响很小。由此可见,胶结破损参量$B_{\sigma }$应力路径相关性明显降低,可作为状态变量引入结构性土体二元介质类本构模型（比如扰动状态概念模型、损伤模型等）表征胶结破损情况,使模型的结构破损规律更为合理。通过数据分析,建议采用指数函数拟合破损参量$B_{\sigma }$与等效塑性应变的关系,

图  8  不同应力路径下破损参量B_σ随等效塑性应变的演化规律

Figure  8.  Evolution of degradation variable B_σversus equivalent plastic strain for DEM sample under different stress paths

下载: 全尺寸图片 幻灯片

式中,$c_{\text{a}}$,$c_{\text{b}}$为胶结破损相关的参数,可通过侧限压缩或常规三轴试验确定。

5.   结论

通过将相对完备的三维胶结接触模型植入粒间接触,反映结构性土粒间胶结的关键特征,建立了结构性土离散元模型；开展了侧限压缩、等向压缩、等应力比压缩以及常规三轴、真三轴等不同加载路径下的离散元数值分析；再现了结构性土体的主要宏观力学性质（尤其是结构性土的破坏特性）；分析了胶结破损参量$B_0$,$B_{\sigma }$随等效塑性应变的演化规律。

（1）胶结破损参量（$\lambda$,$B_0$,$B_{\sigma }$）表征胶结破损程度或结构性衰退程度。当胶结破损参量为0时,土体结构性完好,土体应力全部由胶结组分承担；当胶结破损参量为1时,土体结构性完全破坏,土体应力全部由无胶结组分承担。

（2）胶结破损参量$B_0$的演化受应力比、围压和中主应力系数等的影响,具有应力路径相关性。将应力比和强度包线引入胶结破损参量$B_0$,构建了一个新的胶结破损参量$B_{\sigma }$,应力路径相关性明显降低。在本构模型应用中,胶结破损参量$B_{\sigma }$随等效塑性应变的演化规律可用指数函数表示,并通过侧限压缩或常规三轴试验确定相关参数。

不同应力路径下$B_{\sigma }$的演化仍然具有一定的离散性,今后将考虑改进等效塑性应变的定义和强度包面公式。因为胶结破损规律$B_{\sigma }$基于离散元结果提出,适用于理想胶结颗粒材料（理想结构性土体）,对真实结构性土体的适用性需要进一步验证。