破碎产生碎片的现象广泛存在于自然界和人类生产生活中，材料在撞击、爆破、风化、研磨和压缩等作用下产生的碎片大小不一、形状各异. 了解破碎过程及其产生碎片的大小和形状，对安全防御、爆破开采、级配设计和制粉效率等具有重要意义^[1]. 碎片形状在太空探索^[2]、地貌演变、采矿挖掘、堆石坝填筑^[3]、冲击破碎^[4]和药粉研磨等方面备受关注. 了解碎片形状，有利于评估空间碎片对航天器和卫星的影响，揭示岩石碎片风化剥落、磨损坍塌过程，控制粉体研磨能量需求和碎片数量. 对于岩石工程，可以提高岩石破碎效率，减少能量浪费，生产出粒度和形态可控的碎片，节省二次加工成本.

在堆石坝工程中，爆破开挖的棱角料和开采的卵石料须经人工破碎至设计粒径，满足级配要求才能上坝填筑，颗粒往往因形状和尺寸差异表现出不同的力学特性. 在振动碾压、上部自重和水压力等荷载作用下，堆石颗粒发生磨损、断裂和破碎，导致颗粒级配和形状发生变化，引起颗粒间局部重新排列，进而影响堆石体强度和变形特性^[5-6]. 由于堆石颗粒粒径过大，须按照一定的方法缩尺后进行室内三轴试验，缩尺前后颗粒形状存在差异^[7]，如何评价颗粒形状变化对堆石料变形特性的影响是十分重要的^[8]. 已有研究^[9-10]表明缩尺后的堆石料强度偏高，可能与缩尺方法和颗粒形状有关.

在颗粒破碎过程中会产生大量碎片，受试验方式和测量技术限制，扫描和识别每个碎片须耗费大量时间，鲜有研究关注颗粒破碎后碎片的形状，碎片综合特征描述仍然是难点和挑战^{[8, 11-12]}. X射线CT三维成像技术已用于多种材料的变形和断裂研究，获取原位加载过程中的三维数字化图像数据，进而对颗粒和破碎形成的碎片进行定量研究^[13-16].

研究不同尺寸和材质的岩石颗粒破碎行为，通常是将砂子^[12]和岩石骨料^[7-8]放置在2个刚性板之间进行压缩. 在进行室内试验的同时，使用扩展有限元法（extended finite element method，XFEM）^[17]、离散元法（discrete element method，DEM）^[3]、分子动力学（molecular dynamics，MD）^[18]和连续-离散耦合方法（finite-discrete element method，FDEM）^[1,4,19-20]等进行数值模拟，以扩展对断裂机制的认识.

本研究通过使用X射线CT扫描获取真实岩石颗粒的初始形状，利用FDEM模拟颗粒的断裂和破碎，应用数字图像处理技术表征碎片形状，探索单个岩石颗粒破碎后的碎片形状特征，有助于分析岩石颗粒初始形状对破碎的影响，研究碎片尺寸分布和形状特征，从颗粒和碎片形状角度为破碎过程和碎片现象提供一定参考依据.

为了模拟实际颗粒破碎行为，在数值模型中须具备完整颗粒的初始形状^[15]. 不同成像技术均可用于颗粒形状获取和模型组装重建，例如X射线CT扫描^{[12, 17]}、同步加速断层扫描^[21]和三维激光扫描^[22]等. 本研究选取9个从采石场获得的棱角状花岗岩颗粒，研究其在平板压缩下产生的碎片形状特征，部分扫描颗粒如图1所示.

重点关注碎片的整体形状，忽略断裂引起的碎片表面局部起伏和粗糙变化^{[11, 16]}，采用主成分分析确定颗粒和碎片的主轴和长度^{[4, 11, 14, 16]}. 在获得颗粒和碎片表面积、体积、主轴长度等几何信息后，采用几个简单实用的形状指标量化颗粒和碎片的形状特征，并分析这些指标的统计分布和演变规律. 三维物体的形状描述一般可以采用伸长率 ${\rm{EI}} = {I / L}$和扁平率 ${\rm{FI }}= {S / L}$这2个参数. 同时使用Domokos因子 ${S_{\rm{f}}}$、圆度 ${\psi _{{\rm{3D}}}}$和凸度 ${C_{\rm{x}}}$3个形状指标进行定量分析，定义分别为

(1) ${S_{\rm{f}}} = {{({1 / L} + {1 / I}{{ + 1} / S})\sqrt {{L^2} + {I^2} + {S^2}} } }\left/{ {\sqrt 3 }}\right.,$

(2) ${\psi _{{\rm{3D}}}} = {{{A_{\rm{S}}}} / A} = {{\sqrt[3]{{36{\rm{ {\text{π}}}}{V^2}}}} } \Big/{A}\Big.,$

(3) ${C_{\rm{x}}} = {V / {{V_{{\rm{ch}}}}}}.$

式中： $L$、 $I$、 $S$分别为颗粒或碎片长轴、中轴和短轴长度，对于类圆球状， $L \approx I \approx S$， ${S_{\rm{f}}} \approx 3$，对于扁平状， ${S_{\rm{f}}} > 3$^[11]； $A$、 $V$分别为颗粒或碎片的表面积和体积； ${A_{\rm{S}}}$为与颗粒或碎片体积相等的球体的表面积； ${V_{{\rm{ch}}}}$为包络颗粒或碎片凸包的体积.

9个颗粒的初始形状指标如表1所示.

Dugdale^[23]和Barenblatt^[24]分别基于钢铁和脆性材料的开裂提出了适用性更广的弹塑性断裂力学模型—内聚力模型（cohesive zone model，CZM）. 该理论认为裂纹尖端存在断裂过程区（fracture process zone，FPZ），具有非线性力学行为，在断裂面完全形成前能够传递力. 通过无厚度界面单元（cohesive interface element，CIE）的损伤模拟裂纹扩展，如图2所示（为形象表示，按一定厚度显示），可以反映任意断裂路径下裂缝从萌生、扩展直至完全失效的全过程，并消除裂纹尖端的应力奇异性^{[19, 25]}.

根据界面单元的局部受力情况，断裂可能发生在模式Ⅰ、模式Ⅱ和混合模式Ⅰ/Ⅱ条件下，如图3所示. 图中， $D$为损伤因子， ${k_{\rm{n}}}$、 ${k_{\rm{s}}}$分别为法向、切向刚度. 当法向应力 ${t_{\rm{n}}}$达到界面单元抗拉强度 ${f_{\rm{t}}}$时，模式Ⅰ受拉起裂，随着损伤演化， ${f_{\rm{t}}}$逐渐减小至零，并形成不受力界面. 根据滑动弱化模型，当切向应力 ${t_{\rm{s}}} = \sqrt {t_{{\rm{s1}}}^2 + t_{{\rm{s2}}}^2} $达到界面单元抗剪强度 ${f_{\rm{s}}}$时，模式Ⅱ受剪起裂. 假定 ${f_{\rm{t}}}$为常数， ${f_{\rm{s}}}$由带拉伸截断的Mohr-Coulomb准则确定：

(4) ${f_{\rm{s}}} = \left\{ \begin{array}{l} c - {t_{\rm{n}}}\tan\; {\varphi _i},\;{t_{\rm{n}}} < {f_{\rm{t}}} ;\\ c - {f_{\rm{t}}}{\rm{tan}}\;{\varphi _i},\;{t_{\rm{n}}} \geqslant {f_{\rm{t}}}. \end{array} \right.$

式中： $c$为凝聚力； ${\varphi _i}$为内摩擦角； ${t_{\rm{n}}}$为垂直于剪切方向的法向力，此处拉应力为正，压应力为负.

当 ${t_{\rm{n}}} > {f_{\rm{t}}}$时，拉伸截断将自动产生，假设 ${t_{\rm{s}}}$逐渐减小到纯摩擦阻力 ${f_{\rm{r}}} = - {t_{\rm{n}}}{\rm{tan}}\;{\varphi _{\rm{f}}}$， ${\varphi _{\rm{f}}}$为插入界面单元后的内摩擦角. 然而界面单元往往会在拉、剪应力的复合作用下发生模式Ⅰ/Ⅱ断裂，采用二次应力准则作为拉剪混合破坏的起裂准则：

(5) ${\left\{ {{{\left\langle {{t_{\rm{n}}}} \right\rangle }}/{{{f_{\rm{t}}}}}} \right\}^2} + {\left\{ {{{{t_{\rm{s}}}}}/{{{f_{\rm{s}}}}}} \right\}^2} \geqslant 1.$

式中： $\left\langle {} \right\rangle $为Macaulay括号，当 ${t_{\rm{n}}} > 0$时， $\left\langle {{t_{\rm{n}}}} \right\rangle {\rm{ = }}{t_{\rm{n}}}$，当 ${t_{\rm{n}}} < 0$时， $\left\langle {{t_{\rm{n}}}} \right\rangle {\rm{ = 0}}$.

从能量角度来看，裂纹扩展须克服应力，在断裂过程中会消耗能量，通常使用幂律准则预测混合模式损伤演变：

(6) $ {\left({{G}_{{\rm{I}}}}/{{G}_{{\rm{I}}{\rm{C}}}}\right)}^{2}+{\left({{G}_{{\rm{II}}}}/{{G}_{{\rm{II}}{\rm{C}}}}\right)}^{2}=1.$

式中： $ {G}_{{\rm{IC}}}$、 $ {G}_{{\rm{II}}{\rm{C}}}$分别为Ⅰ型和Ⅱ型临界断裂能.

本研究通过界面单元强度的统计分布和有限元网格的非结构化来反映岩石颗粒微观结构和微观缺陷的非均质性，抗拉强度 ${f_{\rm{t}}}$服从对数正态分布， ${\varphi _i}$和 ${{{f_{\rm{c}}}} / {{f_{\rm{t}}}}}$分别设置为40°和15. FDEM数值模拟涉及的参数如表2所示. 表中， $\rho $为密度， $E$为弹性模量， $\nu $为泊松比， $ {G}_{{\rm{I}}}$、 $ {G}_{{\rm{II}}}$分别为Ⅰ型和Ⅱ型断裂能， $\;\mu $为滑动摩擦系数. 这些参数选自相关研究^{[1, 19, 25]}中的典型值，用于代表典型的岩石颗粒.

刚度系数是内聚力模型中的一个重要参数，为了防止颗粒整体弹性模量降低，界面刚度应充分大以保证实体单元不相互侵入，但是刚度过大将产生如应力伪震荡之类的计算问题. 根据Turon等^[26]提出的计算公式可以近似估算刚度系数：

(7) ${k_{\rm{n}}} = {{\alpha E} / {{l_{\rm{e}}}}},$

(8) ${k_{\rm{s}}} = {{\alpha E} / {\left[ {2\left( {1 + \nu } \right){l_{\rm{e}}}} \right]}}.$

式中： ${l_{\rm{e}}}$为实体单元平均边长； $\alpha $为比例系数，Turon等^[26-27]建议取值 $\alpha \geqslant 50$.

本试验加载速率设置为0.01 m/s（应变速率为1.00~1.67 s⁻¹）以模拟物理试验的准静态加载过程^[27]，时间步长为1.0×10⁻⁷ s.

在模拟颗粒断裂时，由于裂纹尖端前应力梯度较高，只能获得近似的应力场和应变场^[28]. 为了准确反映裂纹尖端的应力和变形状态，有限元网格尺寸必须比断裂过程区的实际尺寸小得多. 由Muskhelishvili解和Westergaard解分别计算出断裂过程区长度的上限值和下限值：

(9) $ {l}_{\rm{fpz}}^{\rm{max}}=3E{G}_{{\rm{I}}}/\left(4{f}_{\rm{t}}^{2}\right),$

(10) $ {l}_{\rm{fpz}}^{\rm{min}}=3{{ {\text{π}}}} E{G}_{{\rm{I}}}/\left(32{f}_{\rm{t}}^{2}\right).$

根据表2中相关参数，可以估算断裂过程区的长度限值为9.74~24.79 mm.

对其中一个岩石颗粒进行有限元网格密度的敏感性分析，4种网格密度对应的四面体单元平均边长分别为5.18、4.48、3.10、2.21 mm，如图4所示. 每个实体单元边缘一般都有2个界面单元，可以根据 ${{2l_{{\rm{fpz}}}^{\min }} / {{l_{\rm{e}}}}}$大致估算断裂过程区中界面单元最小数量 $n_{{\rm{CIE}}}^{\min }$，如表3所示. 表中， ${n_{\rm{n}}}$为节点数， ${n_{\rm{e}}}$为实体单元数. 并对每种网格密度进行10次独立的FDEM单颗粒压缩试验.

由于颗粒形状的不规则和材料的非均质性，荷载位移曲线（ $F$- ${\delta _{\rm{n}}}$）出现了锯齿形波动，如图5所示. 当球形颗粒破碎时，荷载位移曲线通常只有一个明显的峰值，当不规则形状颗粒破碎时通常有几个局部峰值，在每次峰值后，荷载会急剧下降，代表一次局部破碎事件. 分析将颗粒压缩破碎所需要的外力功 $W$（即荷载位移曲线阴影面积），以表征颗粒破碎的难易程度.

在4种网格密度下，同一颗粒独立模拟10次的外力功箱形图如图6所示. 在同一网格密度下，10次独立的FDEM单颗粒压缩试验的外力功具有明显的离散性，这与界面单元强度参数的统计分布特性有关，合理反映出材料的非均质性和微观缺陷的随机分布. 当网格尺寸逐渐减小时，破碎所需外力功逐渐趋于稳定. 从统计学上讲，颗粒破碎所需外力功随单元尺寸的减小而减小，当网格尺寸足够小时，外力功趋于稳定.

Guo等^[29]认为沿着断裂过程区至少有5、6个界面单元就可以消除网格的敏感性，否则会由于断裂过程区内部没有足够的单元而无法准确捕获裂纹尖端前方的应力梯度，导致计算结果不准确. 结合上述网格密度敏感性分析，选择单元平均边长为5.18 mm，断裂过程区中界面单元数量至少为5、6个的网格密度进行后续单颗粒破碎数值模拟.

本研究关注颗粒破碎过程中碎片形状的统计分布和演化规律，以及碎片形状和尺寸之间的相关性. 在FDEM数值试验中，由于采用了较密的网格，颗粒破碎充分，碎片的数量从一百到数百不等，为碎片形状研究提供了足够多的统计数据.

以其中一个颗粒为例介绍单颗粒破碎过程. 颗粒G1荷载位移曲线如图7所示，荷载F与失效界面单元数量 $n_{{\rm{CIE}}}^{\rm{f}}$ 2条曲线均呈锯齿状，两者有较强的正相关性. 小的峰值表示较小的破裂和碎片剥离，大的峰值表示颗粒破裂面急剧扩展贯穿和产生了大块碎片. 荷载位移曲线和失效界面单元数量的变化趋势基本一致，但是失效界面单元数量的演化有一定的滞后现象，说明颗粒破碎不是瞬间完成的，裂纹扩展和碎片重新排列导致延迟破碎现象.

在加载过程中，界面单元不断失效产生新的裂纹面，并产生越来越多的碎片. 失效界面单元累积百分数 $F\left( {n_{{\rm{CIE}}}^{\rm{f}}} \right)$与碎片数量 ${N_{\rm{f}}}$之间呈现出较强的相关性，如图8所示.

加载过程中体积最大的10个碎片如图9所示，通过不同颜色区分，以更加直观地观察断裂模式. 由于碎片间的相互作用，颗粒在压缩过程中，内部发生拉伸、剪切和弯折等不同的断裂机制，产生交叉和分支的裂纹，并延伸发展至颗粒表面形成裂纹面，从而导致断裂破碎，产生少量的大体积碎片和大量的碎屑.

采用碎片的等体积球直径 $d$表征碎片的尺寸大小，分析破碎过程中碎片的尺寸分布，采用Rosin-Rammler（RR）模型拟合碎片尺寸累积分布曲线. RR模型也称为两参数Weibull分布，适用于描述研磨、粉碎、爆破和冲击等破碎过程产生的碎片尺寸分布^[4,30-31]. 正如Åström^[30]所指出的一样，假设颗粒内部和表面的裂缝与现有碎片边缘无关，可以采用RR模型描述碎片尺寸累积分布，该假设与界面单元强度的统计分布较吻合. RR模型为

(11) $F\left( d \right) = 1 - \exp \;\left[ { - {{\left( {d/{d_{\rm{c}}}} \right)}^n}} \right].$

式中： $F\left( d \right)$为尺寸小于 $d$的碎片体积百分数； $n$为分布宽度，较大的 $n$对应较窄范围的碎片尺寸分布； ${d_{\rm{c}}}$为体积分数1−1/e对应的碎片尺寸，即当d=d_c时，F_d=1−1/e，可以用来衡量碎片中碎屑的体积分数，如果碎屑体积分数较高，则 ${d_{\rm{c}}}$期望值较低.

颗粒G1在加载过程中的碎片尺寸累积分布如图10所示，并使用RR模型拟合数据. RR模型拟合效果较好，通过RR模型拟合参数可以进一步分析碎片尺寸分布特征. 随着加载和破碎进行，大尺寸碎片逐渐断裂破碎，小尺寸碎片和碎屑逐渐增多， $n$和 ${d_{\rm{c}}}$逐渐降低，碎片尺寸分布曲线两端逐渐平缓. 表明随着破碎程度加深，碎片尺寸分布更加均匀连续，大尺寸碎片数量逐渐减少.

由于小尺寸碎片和碎屑中含有较少的实体单元，碎片形状主要取决于实体单元的连接方式，在统计分析碎片形状时仅考虑尺寸大于2 mm的碎片. 在加载结束时碎片Domokos因子 ${S_{\rm{f}}}$的概率密度分布如图11所示，其分布特性与物理试验中观察到的相似，概率密度分布 $p\left( {{S_{\rm{f}}}} \right)$在所有试验中都遵循对数正态分布.

在加载结束时碎片的扁平率 ${S / L}$和圆度 ${\psi _{{\rm{3D}}}}$的频率分布 $p(S/L) $、 $p({\psi _{{\rm{3D}}}} )$如图12所示，两者都服从正态分布，均值分别为0.48、0.63，标准差分别为0.19、0.08，可知圆度 ${\psi _{{\rm{3D}}}}$的分布更为集中. 分析碎片凸度 ${C_{\rm{x}}}$的统计分布和演变规律，发现只有尺寸足够大的碎片才会显示局部凹陷， ${C_{\rm{x}}}$的平均值为0.80，小于表1中颗粒初始值.

分析加载过程中的碎片形状的演化规律，当加载位移 ${\delta _{\rm{n}}} = 1\;{\rm{mm}}$时，颗粒G1已经破碎产生48个大于2 mm的碎片. 如图13所示为加载过程中碎片扁平率 ${S / L}$和圆度 ${\psi _{{\rm{3D}}}}$的累积分布 $F(S/L) $、 $F({\psi _{{\rm{3D}}}}) $. 碎片的2个形状指标都比完整颗粒的相应值小，表明破碎产生新的碎片比完整颗粒更具棱角状和更细长. Altuhafi等^[13]和Zhao等^[16]的颗粒破碎研究中也得出了类似结论. 值得关注的是，尽管颗粒在后续加载过程中产生了几次大块破碎，但是碎片形状指标的分布规律没有明显变化，表明破碎过程中断裂机制变化较小，呈现某种自相似性，Zhao等^[16]将此现象归因于岩石高度异质的微观结构. 在颗粒破碎后，大多数碎片的扁平率 ${S / L} \in \left[ {0.37,0.56} \right]$，圆度 ${\psi _{{\rm{3D}}}} \in \left[ {0.58,0.66} \right]$.

颗粒G2、G3和G5在加载过程中荷载位移曲线和碎片数量如图14所示，两者存在明显的相关性，每一个荷载峰值处都对应了一次明显的破碎，碎片数量随破碎程度的加深而增多.

根据荷载位移曲线，将单颗粒破碎分为3种典型类型^[32]：1）第1类（如G2）. 在加载初期，颗粒调整位置达到新的稳定的承载状态，颗粒一旦破碎，便从上到下贯穿，然后彻底丧失承载能力，一般产生少量的大尺寸碎片，属于不完全破碎. 2）第2类（如G3）. 颗粒在加载过程中存在多个荷载峰值，甚至后续峰值会超过第1峰值，在加载时碎片不断重排达到新的承载状态，破碎后产生许多细小碎片，属于完全破碎. 3）第3类（如G5）. 是前两者的过渡形态. 一般认为圆度较好的颗粒加载后即会产生贯穿性的子午面裂纹，继续加载不会再次大规模破碎，因此碎片数量较少，碎片尺寸较大. 本研究发现圆度与碎片数量、碎片尺寸的关系并不完全对应，颗粒G2、G3和G5的圆度相近，但G3完全破碎，G2和G5不完全破碎，说明颗粒初始形状确实会影响碎片形状，但并不起决定性作用.

9个颗粒加载结束后碎片形状指标的累积分布如图15所示，发现碎片的圆度 ${\psi _{{\rm{3D}}}}$、扁平率 ${S / L}$、Domokos因子 ${S_{\rm{f}}}$和凸度 ${C_{\rm{x}}}$的累积分布曲线由紧密变稀疏，数值波动幅度逐渐增大，而初始颗粒形态各异，说明碎片的 ${\psi _{{\rm{3D}}}}$、 ${S / L}$、 ${S_{\rm{f}}}$和 ${C_{\rm{x}}}$对初始颗粒形状的敏感性逐渐增强，当颗粒形状变化时，碎片的凸度 ${C_{\rm{x}}}$变化最明显.

为了量化颗粒形状与碎片形状之间的关系，采用2种函数模型拟合图15中形状指标的累积分布. ${S_{\rm{f}}}$的累积分布采用指数函数拟合：

(12) $F\left( {{S_{\rm{f}}}} \right) = 1 - {\rm{exp}}\;\left[ { - q\left( {{S_{\rm{f}}} - 3} \right)} \right].$

式中： $q$为拟合参数，反映拟合曲线的形状.

其他3个形状指标的累积分布采用Weibull累积分布函数（CDF）拟合：

(13) $F\left( x \right) = 1 - {\rm{exp}}\;\left[ { - {{\left( {{x / a}} \right)}^b}} \right].$

式中： $a$为比例参数， $b$为分布的形状参数. 当 $b > 3$时，式（13）类似于高斯分布；当 $b = 2$时，式（13）为瑞利分布；当 $b = 1$时，式（13）为指数分布.

进行最小二乘法拟合，拟合参数结果如表4所示，颗粒G1加载结束时形状指标的拟合情况如图16所示. 分析颗粒形状指标与表4拟合参数之间的相关性，大多数相关系数不超过0.35，表明颗粒形状与碎片形状之间的相关性较小. 碎片形状除了受颗粒形状影响外，还与破碎类型和颗粒微观结构有关.

进一步探索碎片形状与碎片大小的关系，将9个颗粒破碎产生的共计1 884个碎片混合并按粒径进行分组. 不同粒径组碎片圆度 ${\psi _{{\rm{3D}}}}$的累积分布如图17（a）所示，可以观察到粒径组之间的形状差异明显. 对不同粒径组的数据进行Kruskal - Wallis检验，如图17（b）所示. 图中，绘制了一条水平线，代表每个碎片粒径组中位数99%的置信区间；R为秩，表示所有碎片完整数据中每个形状指标的相对位置；Y轴为碎片粒径组的尺寸范围. 如果置信区间没有重叠，则有99%的置信度表明2组数据的中位数不同. 粒径组 $2 < d < 5$和 $ d>20$的99%置信限用2条垂线表示，置信区间没有重叠，因此Kruskal - Wallis检验以1%的显著水平拒绝2组样本源自相同分布的零假设. 粒径组 $10 < d < 20$和 $ d>20$之间没有显著差异，表明碎片大小与 ${\psi _{{\rm{3D}}}}$没有明确关系. 通过扁平率 ${S / L}$、Domokos因子 ${S_{\rm{f}}}$的累积分布和Kruskal - Wallis检验也能得到相同的结论. 同时可以观察到，大尺寸碎片 ${\psi _{{\rm{3D}}}}$相对较小， ${S / L}$相对较大，总体上表现出较大的碎片往往更具有棱角状，这与Le等^[33]的试验结论一致.

碎片不同粒径组的Zingg形状分类图如图18所示. 该图以 $I/L$和 $S/I$的比值绘制，以0.67为界限将碎片形状分为扁平状、球状、叶片状和柱状4类. 碎片叶片状较少，多数呈球状，整体分布集中于右上侧. 同一粒径组碎片形状分布散乱，无明确趋势表明大小不同的碎片会呈现出某种特定形状.

结合三维光学扫描、数字图像处理技术和FDEM数值模拟，研究不同形状的岩石颗粒在加载过程中的破碎，对颗粒破碎后碎片的大小和形状进行定量分析，主要结论如下.

（1）颗粒有限元网格密度的敏感性表明，断裂过程区至少需要5、6个界面单元以准确描述裂纹尖端的应力梯度和损伤演化，减弱单元尺寸对破碎结果的影响，为采用FDEM方法模拟颗粒破碎提供参考依据.

（2）颗粒加载过程中的3种典型荷载位移曲线和碎片数量表明，颗粒圆度与碎片数量、碎片尺寸的关系并不是完全对应，并使用RR模型拟合碎片尺寸累积分布. 对碎片形状指标进行统计分析，发现初始形状各异的颗粒破碎产生碎片的形状指标累积分布服从同一概率分布函数.

（3）碎片的圆度 ${\psi _{{\rm{3D}}}}$、扁平率 ${S / L}$、Domokos因子 ${S_{\rm{f}}}$和凸度 ${C_{\rm{x}}}$对颗粒初始形状的敏感性逐渐增强. 不同粒径组碎片形状指标的累积分布表明，大小不同的碎片形状差异明显，较大的碎片往往更具有棱角状，但无明确趋势会呈现出某种特定形状.