颗粒形状和中主应力影响钙质砂临界状态特性的离散元分析

Discrete Element Analysis About the Critical State Properties of Calcareous Sand with Considering the Influence of Particle Shape and Medium Principal Stress

摘要

钙质砂临界状态力学特性对其工程性能有着显著影响。作为南海岛礁吹填造陆工程的重要材料，钙质砂土体在天然状态下通常处于三维应力状态，存在明显的中主应力效应。因钙质砂颗粒易破碎，形态各异，对其宏细观临界状态力学特性有显著影响。为深入研究三维应力状态下颗粒形状对钙质砂临界状态特性的影响，在软件PFC^3D中采用碎片替换法构建颗粒破碎模型，开展了一系列真三轴排水剪切模拟，重点分析了颗粒形状和中主应力对数值试样临界状态特性的影响。研究发现：在宏观上，随着颗粒棱角度的增加，试样峰值强度随之增加，达到临界状态时的体积膨胀量更加显著，q‒p空间内临界状态线的斜率和e‒p空间内临界状态线的截距也随之增加（q为广义剪应力，p为平均主应力，e为孔隙比）。随着中主应力系数b的增加，试样峰值强度和临界状态体积膨胀量均逐渐降低，e‒p空间内临界状态线的截距随着b值的增加，呈现先减少后增加的趋势。在细观上，当试样达到临界状态时，随着颗粒棱角度的增加，颗粒体系的平均配位数逐渐减小而接触法向力各向异性逐渐增强，冗余比、孔隙比、平均主应力三者之间呈现良好的线性相关性，在三维空间内逐渐上移。随着b值的增加，试样法向接触力各向异性系数逐渐减小，而试样平均配位数逐渐增大。

Abstract

Objective This study addresses the mechanical challenges associated with calcareous sand in reef construction in the South China Sea. Using the discrete element method, the effects of particle shape and intermediate principal stress on the critical state behavior of calcareous sand are systemically investigated. The objective is to elucidate the underlying deformation and failure mechanisms, providing scientific guidance for foundation design and significant engineering value. Methods The PFC^3D (Particle Flow Code in 3-Dimension) discrete element method was employed to generate crushable particle models with various shapes representative of calcareous sand. Drained true triaxial compression simulations were conducted to evaluate the effects of particle shape and intermediate principal stress on critical state characteristics. Macroscopic properties, such as shear strength, volumetric strain, and critical state line distribution were evaluated. At the microscale, key descriptors such as coordination number, redundancy ratio, and fabric anisotropy were analyzed to characterize the evolving contact network and internal structure. The combined macro-micro analysis elucidates the mechanisms by which particle morphology and intermediate principal stress govern the critical state response of calcareous sand. Results and Discussion The effects of particle shape and angularity on granular material properties were investigated, focusing on stress‒strain relationships, volumetric deformation, critical state lines, and micromechanical responses. Four particle types with increasing angularity—particle a(α = 1.192), b(α = 1.255), c(α = 1.319), and d(α = 1.360)—were analyzed under varying confining pressures (100, 200, 300, 400, and 600 kPa) and intermediate principal stress ratios (b-values: 0, 0.2, 0.4, 0.8). Analysis of the peak stress ratio ((q/p)_max) revealed a significant dependence on particle angularity. At a confining pressure of 100 kPa, (q/p)_max increased from 1.58 for particle a to 1.76 for particle d. Each increment of 0.068 in angularity led to an average increase of 0.18 in (q/p)_max. Variations in the relationship between volumetric strain and major principal strain were also observed. Particle a demonstrated a peak dilatancy rate of 0.39 and a maximum volumetric expansion of 5.81%, whereas particle d had a peak dilatancy rate of 0.70 and a volumetric expansion of 12.71%. Each increase of 0.068 in angularity corresponded to an average increase of approximately 6.90% in volumetric expansion. Critical state line parameters, including the slope (M_cs) and intercept (e_Г​), were also influenced by particle shape and b-value. For b = 0.8, M_cs increased from 0.923 for particle a to 1.069 for particle d, corresponding to an approximate increase of 0.03 for every 0.068 increase in M_cs Similarly, e_Г​​ increased from 1.19848 for particle a to 1.374 67 for particle d, with an average increase of 0.02 per 0.068 increase in angularity. Micromechanical responses demonstrated that the average coordination number (Z_cs) was notably affected by particle shape at both 100 kPa and 400 kPa confining pressures. At 100 kPa, Z_cs decreased from 3.25 for particle a to 2.95 for particle d. At 400 kPa, these values decreased from 4.10 for particle a to 3.75 for particle d. Each increase of 0.068 in angularity led to a decrease of approximately 0.075 in Z_cs at 100 kPa and 0.085 at 400 kPa. The redundancy ratio (R_f) was also influenced by particle shape, with particle a exhibiting a higher redundancy ratio than that of particle d. An increase of 0.068 in angularity resulted in an average decrease of approximately 0.02 in R_f. Additionally, as the void ratio decreased from 0.90 at e_cs= 0.2 to 0.75 at e_cs =0.4, the redundancy ratio decreased. The coefficients of normal contact anisotropy (a_c) and normal contact force anisotropy (a_n ) varied with particle shape. For b = 0, particle a had an a_c of 0.25 and an a_n of 0.35, while particle d had an a_c of 0.22 and an a_n of 0.42. Each increase of 0.068 in angularity resulted in a decrease of approximately 0.007 in a_c and an increase of about 0.01 in a_n. These findings underscore the significant impact of particle shape and angularity on both macroscopic and micromechanical responses, providing new insights into granular material behavior. Conclusions The results indicate that increasing particle angularity enhances peak shear strength and volumetric expansion, while an increase in the intermediate principal stress coefficient (b-value) reduces these responses. With increasing particle angularity, the slope of the critical state line in q-p space and the intercept in the e-p space both increase, with varying trends depending on the b-value. Micromechanical analysis shows that higher particle angularity leads to increased normal contact force anisotropy and a decrease in the average coordination number. Conversely, increasing the b-value reduces anisotropy and promotes higher coordination numbers. A strong linear relationship exists between void ratio and redundancy ratio, with the slope of this relationship decreasing as particle angularity increases. Furthermore, the correlation among void ratio, mean principal stress, and redundancy ratio strengthens with increasing angularity, as evidenced by upward shift of the fitted relationship surface in three-dimensional space.

Graphical abstract

关键词

颗粒形状 / 中主应力 / 颗粒破碎 / 临界状态 / 离散元分析

Key words

particle shape / medium principal stress / particle fragmentation / critical state / DEM analysis

本刊网刊

近年来，随着南海海洋资源开发及“一带一路”战略规划的实施，在南海海域的海上平台及岛礁建设工程逐渐增多。在“绿色工程，生态岛礁”的理念指导下，就地取材选取南海海域周边的珊瑚礁钙质砂作为吹填造陆工程的主要材料可以有效地缩短工期，节约陆地资源。但是由于特殊的生物质成因，钙质砂具有内孔隙丰富、孔隙比大^[1‒2]，棱角鲜明、形状不规则等性质^[3‒4]，因此，在实际工程中，钙质砂颗粒易碎且压缩性高，受波浪荷载和工程结构物的共同作用时易发生破碎与重排，从而显著降低剪切强度。这些特殊物理力学性质导致钙质砂的临界状态等特性与陆源砂相比，存在明显差异^[5]。因此，研究钙质砂的力学特性对实际工程具有重要的指导意义。

临界状态土力学框架的建立是现代土力学的一项重要成就^[6‒7]，研究者基于临界状态土力学开发了很多本构模型^[8‒9]，并将其成功应用于土体力学特性分析^[10‒11]，对实际工程具有重要指导意义。目前，关于颗粒材料临界状态的分析主要针对双轴或者常规三轴试样，这种试样里中主应力σ₂等于大主应力σ₁或小主应力σ₃，分别对应了常规三轴拉伸试验和三轴压缩试验。然而，实际土体通常处于三维应力状态，存在明显的中主应力效应，即在三轴应力状态下，中主应力σ₂不一定等于σ₁或σ₃，从而影响材料的力学行为，如屈服、破坏模式和强度。已有学者采用室内试验或离散元法（DEM）分析了中主应力对可破碎颗粒材料临界状态特性的影响，如：Xiao等^[12]对粗粒土进行的真三轴试验结果表明，e‒lg p（e为孔隙比，p为平均主应力）空间内的临界状态线随着中主应力系数b值的增大而向下移动；Liu等^[13]对由582个可破碎颗粒组成的离散元数值试样进行了真三轴模拟，结果表明，随着b值的增加颗粒破碎呈指数增长；Shi等^[14‒15]针对钙质砂进行了真三轴剪切模拟，讨论了b值和颗粒破碎对临界状态线的耦合影响，结果表明，e‒lg p空间内临界状态线的斜率和截距均随b值的增大而减小，q‒p空间内临界状态线的斜率随b值增大而减小，b为0.75时达到最小值（q为广义剪应力）。但上述研究均没有讨论颗粒形状对可破碎颗粒材料临界状态特性的影响，已有研究表明^[16‒17]，由于钙质砂颗粒多为块状、棒状、枝状等不规则形状^[18‒19]，颗粒形状将显著影响试样的强度和变形特性。

本文基于颗粒离散元法^[20]，通过PFC^3D数值平台构建不同颗粒形状的可破碎颗粒体来近似模拟钙质砂颗粒，对由不同颗粒形状颗粒体所构成的数值试样，开展排水条件下真三轴剪切数值模拟^[21]，重点分析颗粒形状和中主应力对数值试样临界状态特性的影响。通过分析剪切强度、体积变化和临界状态线分布等宏观特性，以及配位数、冗余比和组构各向异性系数等细观特性，探索宏细观之间的内在关联，揭示颗粒形状和中主应力影响钙质砂临界状态的宏微观力学特性关联机制，为本构模型中引入细观参数提供了物理机制，为南海钙质砂地基的设计和建设提供指导和借鉴。

1 离散元建模

1.1 颗粒形状的构建

为了近似体现实际钙质砂颗粒的不同形状，在PFC^3D中，采用团聚颗粒（Clump）方法，构建4种不同颗粒形状的颗粒，如图1所示。每种颗粒形状都由10个直径相同的基本球构成，且每个基本球两两相切。同时，保证4种颗粒形状不同的颗粒，其内孔隙大小保持一致。

为了定量描述颗粒形状，引入颗粒棱角度α。为了方便计算α，取图1中4种颗粒形状的外轮廓，如图2所示。

将棱角度定义为所有接触面夹角的外角的平均值，图3为其计算示意图。以图2（a）为例，给其顶点编上字母，如图3（a）所示，求其面ABD与其接触面ACD的夹角的外角β₁（图3（b））；依此类推，计算所有接触面夹角的外角，并取其平均值，即可得到棱角度α。经计算，图1中，从颗粒a到颗粒d，α的数值分别为1.192、1.255、1.319、1.360，逐次增加。

1.2 颗粒破碎建模

颗粒破碎建模采用颗粒破碎准则和碎片替代法，颗粒破碎准则采用八面体剪应力准则，将单粒压碎试验获取的特征拉应力设定为破碎准则的应力阈值。特征拉应力{\sigma }_{\mathrm{f}}可表示为^[22]：

{\sigma }_{\mathrm{f}}=\frac{F_{\mathrm{f}}}{d^{\mathrm{2}}}（1）

式中，F_f为单颗粒压缩试验中颗粒破碎时所受到的径向力，d为颗粒粒径。

作用在颗粒上的二阶应力张量{\sigma }_{ij}可表示为：

{\sigma }_{ij}=\frac{\mathrm{1}}{V}{\displaystyle \sum _{c=\mathrm{1}}^{N_c}}{f}_j^c{l}_i^c（2）

式中，{\sigma }_{ij}为应力张量，N_c 为颗粒的接触总数目，f_j^c为接触点c处j方向接触力分量，l_i^c为颗粒中心位置到接触点c向量在i方向的分量，V为颗粒体积。

作用在颗粒上的八面体剪应力q_oct可表示为：

q_{\mathrm{o}\mathrm{c}\mathrm{t}}=\frac{\mathrm{1}}{\mathrm{3}}\sqrt[]{({\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{2}}{)}^{\mathrm{2}}+({\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}}{)}^{\mathrm{2}}+({\sigma }_{\mathrm{2}}-{\sigma }_{\mathrm{3}}{)}^{\mathrm{2}}}（3）

一般来说，单个颗粒的特征强度遵循Weibull分布^[23]，可表达为：

{\sigma }_{\mathrm{n}}={\sigma }_{\mathrm{0}}{\left(\frac{d}{d_{\mathrm{0}}}\right)}^{-\frac{\mathrm{3}}{m}}（4）

式中：{\sigma }_{\mathrm{n}}为颗粒破碎强度，下标n表示破碎发生时的临界状态；d_{\mathrm{0}}为特征颗粒的直径；{\sigma }_{\mathrm{0}}为颗粒特征破碎强度；m为Weibull系数。

八面体剪应力准则中，颗粒的八面体剪应力强度q_{\mathrm{n}}与单颗粒破碎强度σ_n存在如下关系：

q_{\mathrm{n}}=\mathrm{0.9}{\sigma }_{\mathrm{n}}=\mathrm{0.9}{\sigma }_{\mathrm{0}}{\left(\frac{d}{d_{\mathrm{0}}}\right)}^{-\frac{\mathrm{3}}{m}}（5）

式中，当q_{\mathrm{o}\mathrm{c}\mathrm{t}}\ge q_{\mathrm{n}}时，颗粒发生破碎。

单颗粒破碎模式采用碎片替代法，如图4所示。首先，由10个相同粒径的球形颗粒组合形成一个团聚颗粒（Clump），代表初始非球形砂土颗粒（图4（a））；当初始非球形颗粒所承受的八面体剪应力满足前文所述的破碎准则时，颗粒发生破碎，破碎成10个球形子颗粒（图4（b））；在加载过程中，球形子颗粒还可以继续破碎，破碎形式为14个粒径不等且满足阿波罗填充模式的球形颗粒^[18]（图4（c））。需要说明的是，在颗粒发生2级及以下破碎时，会发生一定的体积损失，但本次模拟中发生2级及以下破碎的颗粒数量不超过颗粒总数的8.5%，由此带来的体积损失不超过颗粒总体积的2%，因此，这个体积损失可以忽略不计。

1.3 数值模型的验证

以蔡正银等^[24]的实际钙质砂室内三轴排水剪切试验结果为参照，对本文DEM模型进行验证。在PFC^3D中，由六面刚性墙围成一个初始边长为5.5 mm的立方体试样，试样由4种形状的颗粒（图1）等比例随机混合，粒径范围为0.25～2.00 mm，与文献[24]所用的实际钙质砂粒径基本一致，如图5所示。为了避免由颗粒内孔隙消散导致的孔隙比波动，本文DEM模型孔隙比通过固体颗粒体积与总体积计算得到，但是数值试样颗粒存在内孔隙，导致数值试样的孔隙比略大于相同级配曲线的室内试验结果。控制数值试样的初始孔隙比为0.931，与文献[24]实际钙质砂试样一致,生成的数值试样如图6所示。

数值模拟中颗粒密度设置为2 780 kg/m³，颗粒间的接触模型采用线性刚度模型，颗粒间法向和切向接触刚度均设为50 000 N/m。为确保试样中颗粒中的准静态平衡，局部阻尼系数和旋转阻抗系数均设为0.5，颗粒间摩擦系数设为0.42。颗粒破碎参数取自Ma等^[25]钙质砂单颗粒压碎试验，颗粒破碎强度{\sigma }_{\mathrm{0}}取为5.67 MPa，颗粒特征粒径取为4.17 mm。

对图6中的数值试样施加100 kPa的均等围压，随后在保持侧向100 kPa围压不变的条件下，对顶面和底面墙体同时施加相对运动的速度来实现三轴加载。图7对比了数值试样与室内钙质砂试样的偏应力‒轴应变、体应变‒轴应变关系。图7中，体应变正值表示体积膨胀，负值表示体积收缩。由图7可知，三轴数值试样表现出的力学响应与实际钙质砂试样基本吻合，DEM试样的峰值强度比室内试样的强度稍低，原因在于数值模拟中使用的颗粒形状相对规则且只有4种，而实际钙质砂颗粒的形状更加复杂多样。此外，采用碎片替代法时，子颗粒破碎引起的体积损失也会对强度产生影响。同时，数值试样与实际钙质砂的级配存在一定差异，这进一步导致了强度的差异。

为了量化颗粒破碎程度，Hardin^[26]提出了相对破碎率B_r的概念：

B_{\mathrm{r}}=\frac{B_{\mathrm{t}}}{B_{\mathrm{p}}}（6）

式中，B_t为总破碎势，B_p为初始破碎势。相对破碎势B_r的下限值为0，表示没有颗粒破碎发生；理论上限值为1，表示所有颗粒都将破碎为0.074 mm以下的细颗粒。

对比本文数值试样与蔡正银等^[24]实际钙质砂三轴试验的破碎率，当轴向应变均达到25%时，实际钙质砂的B_r为0.104 4，数值试样的B_r为0.110 2，两者较为接近，也进一步验证了本文颗粒破碎模型的有效性。

1.4 真三轴排水剪切建模

中主应力系数b可表示为：

b=\frac{({\sigma }_{\mathrm{2}}-{\sigma }_{\mathrm{3}})}{({\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}})}（7）

如图6所示，在剪切过程中上部墙体和底部墙体沿着z方向以0.05 m/s的恒定速度向内运动，大主应力{\sigma }_{\mathrm{1}}由上部墙体和底部墙体的平均应力计算得到。在z方向的加载过程中，x方向和y方向始终采用伺服控制。围压{\sigma }_{\mathrm{3}\mathrm{c}}作为x方向的目标伺服应力以确保整个加载过程中小主应力不变，即{\sigma }_{\mathrm{3}}={\sigma }_{\mathrm{3}\mathrm{c}}，中主应力{\sigma }_{\mathrm{2}}可由式（7）计算得到，将计算所得的{\sigma }_{\mathrm{2}}作为y方向的目标伺服应力，从而保证在剪切过程中中主应力系数b值保持恒定。

本次数值模拟中，首先构建4种不同颗粒形状的试样；对于每种颗粒形状的试样，施加5种不同的围压，即100、200、300、400、600 kPa；对于每级围压，均开展了b值为0、0.2、0.4、0.8共4种不同应力条件的排水剪切模拟。

剪切过程中，试样的广义剪应变ε_d、广义剪应力q、平均主应力p可表示为：

{\epsilon }_{\mathrm{d}}=\frac{\sqrt[]{\mathrm{2}}}{\mathrm{3}}\sqrt[]{({\epsilon }_{\mathrm{1}}-{\epsilon }_{\mathrm{2}}{)}^{\mathrm{2}}+({\epsilon }_{\mathrm{2}}-{\epsilon }_{\mathrm{3}}{)}^{\mathrm{2}}+({\epsilon }_{\mathrm{3}}-{\epsilon }_{\mathrm{1}}{)}^{\mathrm{2}}}（8）

q=\frac{\mathrm{1}}{\sqrt[]{\mathrm{2}}}\sqrt[]{({\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{2}}{)}^{\mathrm{2}}+({\sigma }_{\mathrm{2}}-{\sigma }_{\mathrm{3}}{)}^{\mathrm{2}}+({\sigma }_{\mathrm{3}}-{\sigma }_{\mathrm{1}}{)}^{\mathrm{2}}}（9）

p=\frac{{\sigma }_{\mathrm{1}}+{\sigma }_{\mathrm{2}}+{\sigma }_{\mathrm{3}}}{\mathrm{3}}（10）

式（8）中，ε₁、ε₂、ε₃分别为试样的大、中、小主应变。

2 结果分析

2.1 宏观力学响应

2.1.1 应力‒应变关系

图8为不同b值下，4种不同颗粒形状数值试样的广义剪应力比与大主应变的变化关系，限于篇幅，仅给出围压100 kPa时的结果。由图8可知，颗粒形状会影响试样的峰值剪切强度，无论何种b值，广义剪应力比的峰值(q/p)_max均随着颗粒棱角度的增加而增大。以b=0情况为例，当棱角度α从1.192增加至1.360时，(q/p)_max从1.58增大至1.76，增幅达11.39%。Azema等^[27]认为颗粒外表面的面数越多，颗粒棱角度就会变大，剪切强度也会随之增加。本文中，4种颗粒形状面数分别为4、9、12、14，因此棱角度也依次增加。随着b值的增加，试样(q/p)_max逐渐降低；以试样α=1.192为例，当b从0增加至0.8时，(q/p)_max从1.58降至1.12，降幅约30%。

图9为不同b值下，4种不同颗粒形状数值试样的体应变随大主应变的变化关系，本文仅给出了围压100 kPa时的结果。由图9可知，无论何种b值，试样的峰值剪胀率(dε_v/dε₁)_max和达到残余状态（ε₁=40%）时的体积膨胀量均随着颗粒棱角度α的增加而增加。以b=0为例，当α从1.192增加至1.360时，试样的峰值剪胀率从0.39增至0.70，达到残余状态时的体积膨胀量从5.81%增至12.71%。吴杨等^[28]曾以石英砂与玻璃圆珠或玻璃碎珠采用不同比例相混合，并对试样施加真三轴剪切，试验结果表明，与玻璃圆珠相比，石英砂与玻璃碎珠相混合的试样，其剪胀性更为显著，其机理为在剪切过程中，颗粒会翻越与之相邻的颗粒，不规则颗粒在这种翻越过程中，错动幅度大，从而诱发更强的剪胀。本文采用4种颗粒形状，颗粒外轮廓面数逐渐增加，面数越多形状也越不规则，从而导致试样的剪胀越显著。同时，随着b值的增加，试样的临界状态体积膨胀量逐渐降低。以α=1.192试样为例，当b值从0增加至0.8时，试样的体积剪胀量从13.12%降低至8.81%，降幅约33%。

2.1.2 临界状态线

参考Zhou等^[29]研究结果，可近似认为当轴向应变达到40%时，试样达到临界状态。图10给出了不同b值下，q‒p空间内的临界状态线分布。由图10可知，临界状态时的偏应力q_cs与平均主应力p_cs之间存在良好线性关系，可表示为：

q_{\mathrm{c}\mathrm{s}}=M_{\mathrm{c}\mathrm{s}}{p}_{\mathrm{c}\mathrm{s}}（11）

式中，M_cs为临界状态线斜率，也称为临界状态应力比。

由图10可知，颗粒棱角度对q‒p空间内临界状态线的分布有一定影响。以b=0.8为例，α为1.192、1.255、1.319、1.360时，4种试样的临界状态线斜率M_cs分别为0.923、0.937、0.951、1.069，依次递增。但中主应力系数b对上述规律有一定的耦合影响，当b值越小时，颗粒棱角度α对M_cs的影响越小。

按照Li等^[30]的建议，在e‒p空间内，临界状态线可表示为：

e_{\mathrm{c}\mathrm{s}}=e_{\mathrm{\Gamma }}-\lambda {\left(\frac{p_{\mathrm{c}\mathrm{s}}}{p_{\mathrm{a}}}\right)}^{\xi }（12）

式中：e_cs为试样临界状态时的孔隙比；e_Г为当p_cs为0时，试样的临界状态孔隙比，即e‒p平面临界状态线的截距；λ为e‒p平面的临界状态线的斜率；e_Г和λ两个参数决定着临界状态线所处的位置；p_a为大气压强，本文取101 kPa；ξ为材料常数，本文取0.7。

图11为不同b值下，4种试样达到临界状态时的孔隙比e_cs与平均主应力p_cs之间的关系，由式（12）计算获得。由图11可知，式（12）对试样临界状态点拟合效果较好，不同颗粒形状试样之间的临界状态线近似接近平行且临界状态线的截距与棱角度有一定关联。表1为不同b值下的e‒p空间内临界状态线的截距e_Г数值。由表1可知，4种颗粒形状试样的临界状态线的截距随着颗粒棱角度的增加而增大，且这一规律不受b值大小的影响。

图12为e‒p空间内临界状态线的截距及斜率与b值之间的变化关系。由图12（a）可知，随着b值的增加，4种颗粒形状试样的截距基本都呈现出先降低后增加的趋势。在b值为0.2至0.4时达到最小值，这一趋势表明颗粒形状对临界状态线截距有显著影响。由图12（b）可知，随着b值的增加，4种颗粒形状试样的临界状态线的斜率呈现逐渐降低的趋势。不同颗粒形状试样在相同b值下的临界状态线斜率存在差异，例如在b值为0.4时，斜率分别为0.049、0.050、0.046、0.048，表明了颗粒形状和中主应力对试样临界状态线分布的耦合影响。

2.2 细观响应

2.2.1 平均配位数

配位数是指试样中某一颗粒与其相邻颗粒的接触数目，是定量描述颗粒材料内部结构和稳定性的重要指标。颗粒体系中平均每个颗粒的接触数目可以用平均配位数Z表示：

Z=\frac{\mathrm{2}{N}_c}{N_{\mathrm{P}}}（13）

式中，N_p是总颗粒数目。

图13（a）和（b）分别为围压σ_3c=100和400 kPa下，试样达到临界状态时，不同棱角度颗粒试样平均配位数Z_cs随b值的变化规律。分析图13可得几点规律。1）围压对Z_cs的影响为，随着围压的增加，Z_cs显著增大；当σ_3c=100 kPa时，Z_cs的范围为5.259～5.997；当σ_3c=400 kPa时，Z_cs的范围增至6.185～6.793，表明较高的围压条件下，试样的平均配位数会显著增加。围压对Z_cs影响较大，围压越大，试样达到临界状态时的平均配位数越大。2）b值对Z_cs的影响为，在不同围压下，Z_cs均随b值的增加而增大。以试样α=1.360为例，当σ_3c=100 kPa时，b值从0增加到0.8，Z_cs从5.259增加到5.832；在σ_3c=400 kPa时，Z_cs由6.185增至6.733。这一结果显示，随着b值增大，试样的平均配位数增加。3）颗粒棱角度α对Z_cs的影响：Z_cs随着颗粒棱角度α的增加而减小。在b值为0.4且σ_3c=400 kPa的条件下，α从1.192增加到1.360，Z_cs从6.793降至6.524。较大的颗粒棱角度导致更高的破碎率，从而降低了Z_cs。

2.2.2 冗余比

冗余比常被应用于颗粒材料微观结构分析，以描述颗粒之间的接触行为和配位数的演化。对于类固体颗粒材料，冗余比的值大于1，冗余比为1是类固体材料和类液体间的界限。冗余比R_f可表示为^[30]：

R_{\mathrm{f}}=\frac{Z(\mathrm{5}{R}_{\mathrm{n}\mathrm{s}\mathrm{r}}+\mathrm{3}{R}_{\mathrm{s}}+\mathrm{3}{R}_{\mathrm{r}}+R_{\mathrm{s}\mathrm{r}})}{\mathrm{12}}（14）

式中，R_nsr为既没有旋转又没有滑移的接触数目与接触总数的比值，R_s为滑移的接触数目与接触总数的比值，R_r为旋转的接触数目与接触总数的比值，R_sr为既滑移又旋转的接触数目与接触总数的比值。

图14显示了试样达到临界状态时，不同颗粒形状试样的冗余比R_fcs和孔隙比e_cs之间的关系。史旦达等^[15]指出，对于给定的DEM试样，达到临界状态时其冗余比和孔隙比的关系不受加载路径即b值的影响，因此，在这里只给出了b=0.8的结果。由图14可知，随着e_cs的增加，R_fcs逐渐减小，以试样α=1.192为例，e_cs从0.824增加到1.075，R_fcs从1.771降至1.385，这表明孔隙比的增加使得冗余比明显下降；不同颗粒形状试样的R_fcs和e_cs基本均可拟合成线性关系，且拟合直线的斜率随着颗粒棱角度α的增加而减小，如棱角度从1.192增加至1.360时，拟合直线的斜率从-1.18降低至-1.45，降幅约23%。

图15给出了临界状态时，孔隙比e_cs、平均主应力p和冗余比R_fcs之间的关系。

分析图15可知，当试样达到临界状态时，上述三者在三维空间内呈现良好的线性相关性；且随着颗粒棱角度的增加，拟合直线在空间内逐渐上移。

2.2.3 组构各向异性

除了冗余比外，由剪切荷载导致的组构各向异性演化是另一个需要研究的重要细观现象。通常情况下，接触各向异性可以用接触法线及分支矢量各向异性来描述，力各向异性可由切向力各向异性和法向力各向异性来描述。接触法向二阶张量{\varphi }_{ij}^{}为：

{\varphi }_{ij}={\int }_{\mathrm{\Theta }}E\left(\mathrm{\Theta }\right)n_i{n}_j\mathrm{d}\mathrm{\Theta }=\frac{\mathrm{1}}{N_c}\sum _{c\in N_c}{n}_i{n}_j（15）

式中，n_i 为单位接触法向在i方向上的分量，n_j 为单位接触法向在j方向上的分量，E(Θ)为概率密度函数，\mathrm{\Theta }为积分域。

使用二阶傅里叶展开：

E\left(\mathrm{\Theta }\right)=\frac{\mathrm{1}}{\mathrm{4}\mathrm{\pi }}(\mathrm{1}+a_{ij}^{\mathrm{c}}{n}_i{n}_j)（16）

式中，a_{ij}^{\mathrm{c}}为接触法向各向异性张量,可进一步表示为：

a_{ij}^{\mathrm{c}}=\mathrm{7.5}{\varphi }_{ij}^{\text{'}}（17）

式中，{\varphi }_{ij}^{\text{'}}为{\varphi }_{ij}^{}偏张量的部分。

接触法向、法向接触力、切向接触力各向异性系数可统一表示为：

a_{\mathrm{*}}=\sqrt[]{\frac{\mathrm{3}}{\mathrm{2}}{a}_{ij}^{\mathrm{*}}{a}_{ij}^{\mathrm{*}}}（18）

式中，*可为c、n和t，即a_c为接触法向的各向异性系数，a_n为法向接触力的各向异性系数，a_t为切向接触力各向异性系数。

图16为围压σ_3c=100 kPa时，不同b值条件下4种颗粒形状试样的各向异性系数（a_c、a_n）演化规律。由于a_t的值很小，其最大值仅为0.03，对颗粒材料抗剪强度的影响可以忽略，因此，没有给出a_t的演化过程。分析图16可得几点规律：1）剪切过程中，a_n的值总是大于a_c值，说明在颗粒材料体系剪切强度发挥过程中，法向接触力各向异性占主导作用；2）颗粒形状对a_n的演化规律有一定影响，a_n基本随颗粒棱角度α的增加而增大，这与第2.1.1节中试样剪切强度随着α的增加而增大的结论相一致；3）与a_n值相比，试样的a_c的值受颗粒形状的影响不大，例如，不同棱角度下，a_c值差距稳定在0.2以内；4）随着中主应力系数b值的增加，试样的a_n值逐渐减小，从3.0降至2.0附近，而a_c值基本不受b值变化的影响，一直稳定在1.0至1.5之间。

3 结 论

针对实际钙质砂不规则形状，本文基于PFC^3D构建了非球形颗粒破碎模型，研究了颗粒形状和中主应力对钙质砂临界状态的影响及机理，主要结论有：

1）试样峰值抗剪强度随着颗粒棱角度的增加而增大，随着中主应力系数b值的增加而减小；试样的峰值剪胀率和达到临界状态时的体积膨胀量均随着颗粒棱角度的增加而增大。

2）随着颗粒棱角度的增加，q‒p空间内临界状态线的斜率和e‒p空间内临界状态线的截距均逐渐增加；随着b值的增加，e‒p空间内临界状态线的截距呈现先减少后增加的趋势，斜率逐渐降低。

3）达到临界状态时，随着颗粒棱角度的增加，试样法向接触力各向异性系数逐渐增大，而试样平均配位数逐渐减小；随着b值的增加，试样法向接触力各向异性系数逐渐减小，而试样平均配位数逐渐增大。

4）达到临界状态时，随着试样细观冗余比与宏观孔隙比之间存在良好的线性关系，且随着颗粒棱角度的增加，拟合直线的斜率逐渐减小；孔隙比、平均主应力和冗余比三者之间呈现良好的线性相关性，且随着颗粒棱角度的增加，拟合直线在三维空间内逐渐上移。

参考文献

[1] Xiong Hao, Qiu Yuanyi, Liu Jiaying,et al.Macro-microscopic mechanism of suffusion in calcareous sand under tidal fluctuations by coupled CFD-DEM[J].Computers and Geotechnics,2023,162:105676. doi:10.1016/j.compgeo.2023.105676

[2] Haibo Lyu, Gu Jianxiao, Zhou Jianrong,et al.Mechanical behavior and particle breakage of calcareous sand in triaxial test[J].Marine Geophysical Research,2023,44(3):18. doi:10.1007/s11001-023-09525-0

[3] Chen Shasha, Zhang Junhui, Long Zhilin,et al.Effects of particle size on the particle breakage of calcareous sands under impact loadings[J].Construction and Building Materials,2022,341:127809. doi:10.1016/j.conbuildmat.2022.127809

[4] Chen Huibin, Hou Heying, Jiang Pengming,et al.Study on modified dynamic constitutive model for calcareous sand considering particle breakage[J].Geofluids,2022,2022(1):2755102. doi:10.1155/2022/2755102

[5] Xu Liangjie, Wang Xinzhi, Wang Ren,et al.Physical and mechanical properties of calcareous soils:A review[J].Marine Georesources & Geotechnology,2022,40(6):751‒766. doi:10.1080/1064119x.2021.1927270

[6] Ayala J, Fourie A, Reid D.A unified approach for the analysis of CPT partial drainage effects within a critical state soil mechanics framework in mine tailings[J].Journal of Geotechnical and Geoenvironmental Engineering,2023,149(6):04023036. doi:10.1061/jggefk.gteng-10915

[7] Fedorowicz J, Fedorowicz L, Kadela M.Description of limit states in the subsurface layer of loosened subsoil in view of critical state soil mechanics[J].Materials,2021,14(23):7288. doi:10.3390/ma14237288

[8] Liu Haoyuan, Nagula S, Petter Jostad H,et al.Considerations for using critical state soil mechanics based constitutive models for capturing static liquefaction failure of tailings dams[J].Computers and Geotechnics,2024,167:106089. doi:10.1016/j.compgeo.2024.106089

[9] Rahman M M, Nguyen H B K, Fourie A B,et al.Critical state soil mechanics for cyclic liquefaction and postliquefaction behavior:DEM study[J].Journal of Geotechnical and Geoenvironmental Engineering,2021,147(2):04020166. doi:10.1061/(asce)gt.1943-5606.0002453

[10] Su Yonghua, Wang Dong.An analysis of direct shear test results of sand-gravel mixture based on the discrete element method[J].Hydrogeology & Engineering Geology,2021,48(6):97‒104. doi:10.16030/j.cnki.issn.1000-3665.202012029

[11] 苏永华,王栋.基于离散元法的砂石混合体直剪试验结果分析[J].水文地质工程地质,2021,48(6):97‒104. doi:10.16030/j.cnki.issn.1000-3665.202012029

[12] Cherif Taiba A, Mahmoudi Y, Azaiez H,et al.Impact of the overall regularity and related granulometric characteristics on the critical state soil mechanics of natural sands:A state-of-the-art review[J].Geomechanics and Geoengineering,2023,18(4):299‒308. doi:10.1080/17486025.2022.2044076

[13] Xiao Yang, Liu Hanlong, Desai C S,et al.Effect of intermediate principal-stress ratio on particle breakage of rockfill material[J].Journal of Geotechnical and Geoenvironmental Engineering,2016,142(4):06015017. doi:10.1061/(asce)gt.1943-5606.0001433

[14] Liu Yiming, Liu Huabei, Mao Haijun.DEM investigation of the effect of intermediate principle stress on particle breakage of granular materials[J].Computers and Geotechnics,2017,84:58‒67. doi:10.1016/j.compgeo.2016.11.020

[15] Shi D D, Cao D, Deng Y B,et al.DEM investigations of the effects of intermediate principal stress ratio and particle breakage on the critical state behaviors of granular soils[J].Powder Technology,2021,379:547‒559. doi:10.1016/j.powtec.2020.10.094

[16] Shi Danda, Cao Dong, Xue Jianfeng,et al.DEM studies on the effect of particle breakage on the critical state behaviours of granular soils under undrained shear conditions[J].Acta Geotechnica,2022,17(11):4865‒4885. doi:10.1007/s11440-022-01580-y

[17] Wang Xing, Wang Ying, Liao Cheng,et al.Particle breakage mechanism and particle shape evolution of calcareous sand under impact loading[J].Bulletin of Engineering Geology and the Environment,2022,81(9):372. doi:10.1007/s10064-022-02868-9

[18] Touiti L, Kim T, Jung Y H.Analysis of calcareous sand particle shape using Fourier descriptor analysis[J].International Journal of Geo-Engineering,2020,11(1):15. doi:10.1186/s40703-020-00122-3

[19] Li Xue, Liu Jiankun, Li Jinze.Correction to:Fractal dimension,particle shape,and particle breakage analysis for calcareous sand[J].Bulletin of Engineering Geology and the Environment,2022,81(5):197. doi:10.1007/s10064-022-02704-0

[20] Rui Shengjie, Guo Zhen, Si Tongling,et al.Effect of particle shape on the liquefaction resistance of calcareous sands[J].Soil Dynamics and Earthquake Engineering,2020,137:106302. doi:10.1016/j.soildyn.2020.106302

[21] Kodicherla S P K.Discrete element modelling of granular materials incorporating realistic particle shapes[J].International Journal of Geo-Engineering,2023,14(1):15. doi:10.1186/s40703-023-00193-y

[22] Zhang Chao, Fu Jianxin, Wang Yu.Numerical simulation of true triaxial unilateral unloading effect of fractured rock[J].Engineering Computations,2023,40(9/10):2110‒2128. doi:10.1108/ec-07-2022-0477

[23] McDowell G R, Bolton M D.On the micromechanics of crushable aggregates[J].Géotechnique,1998,48(5):667‒679. doi:10.1680/geot.1998.48.5.667

[24] Weibull W.A statistical distribution function of wide applicability[J].Journal of Applied Mechanics,1951,18(3):293‒297. doi:10.1115/1.4010337

[25] Cai Zhengyin, Hou Heying, Zhang Jinxun,et al.Critical state and constitutive model for coral sand considering particle breakage[J].Chinese Journal of Geotechnical Engineering,2019,41(6):989‒995. doi:10.11779/CJGE201906001

[26] [蔡正银,侯贺营,张晋勋,等.考虑颗粒破碎影响的珊瑚砂临界状态与本构模型研究[J].岩土工程学报,2019,41(6):989‒995. doi:10.11779/CJGE201906001

[27] Ma Linjian, Li Zeng, Wang Mingyang,et al.Effects of size and loading rate on the mechanical properties of single coral particles[J].Powder Technology,2019,342:961‒971. doi:10.1016/j.powtec.2018.10.037

[28] Hardin B O.Crushing of soil particles[J].Journal of Geotechnical Engineering,1985,111(10):1177‒1192. doi:10.1061/(asce)0733-9410(1985)111:10(1177)

[29] Azéma E, Radjai F, Dubois F.Packings of irregular polyhedral particles:Strength,structure,and effects of angularity[J].Physical Review E,Statistical,Nonlinear,and Soft Matter Physics,2013,87(6):062203. doi:10.1103/physreve.87.062203

[30] Wu Yang, Rong Haojun, Wang Jinlian,et al.A true triaxial experimental study on the coupled effect of particle shape and intermediate principal stress on the mechanical properties of sand[J].Chinese Journal of Rock Mechanics and Engineering,2023,42(2):497‒507. doi:10.13722/j.cnki.jrme.2022.0721

[31] [吴杨,容浩俊,王金莲,等.颗粒形状和中主应力对砂土力学特性耦合影响的真三轴试验研究[J].岩石力学与工程学报,2023,42(2):497‒507. doi:10.13722/j.cnki.jrme.2022.0721

[32] Zhou Wei, Liu Jiaying, Ma Gang,et al.Three-dimensional DEM investigation of critical state and dilatancy behaviors of granular materials[J].Acta Geotechnica,2017,12(3):527‒540. doi:10.1007/s11440-017-0530-8

[33] Li X S, Wang Y.Linear representation of steady-state line for sand[J].Journal of Geotechnical and Geoenvironmental Engineering,1998,124(12):1215‒1217. doi:10.1061/(asce)1090-0241(1998)124:12(1215)

基金资助

国家自然科学基金项目(52571297)

国家自然科学基金项目(52471290)

国家自然科学基金项目(52301327)

上海市自然科学基金项目(24ZR1427900)

上海市自然科学基金项目(23ZR1426200)