力学学报  2019 , 51 (1): 16-25 https://doi.org/10.6052/0459-1879-18-340

颗粒材料计算力学专题

颗粒材料破碎演化路径细观热力学机制1)

沈超敏, 刘斯宏2)

河海大学水利水电学院,南京 210098

EVOLUTION PATH FOR THE PARTICLE BREAKAGE OF GRANULAR MATERIALS: A MICROMECHANICAL AND THERMODYNAMIC INSIGHT1)

Shen Chaomin, Liu Sihong2)

Department of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing 210098, China

中图分类号:  TV16

文献标识码:  A

版权声明:  2019 力学学报期刊社 力学学报期刊社 所有

基金资助:  1) 国家自然科学基金(U1765205),国家重点研发计划(2017YFC0404800)以及中央高校基本科研业务费专项资金(2018B40914)资助项目.

作者简介:

作者简介: 2) 刘斯宏,教授,博导,主要研究方向:土石坝工程、粒状体力学与土工袋技术. E-mail: Sihongliu@hhu.edu.cn

展开

摘要

颗粒材料在高应力环境下会发生颗粒破碎现象,颗粒破碎不仅影响颗粒材料的力学特性,同时与大量工程问题密切相关.目前的相关研究主要集中在唯象地描述颗粒破碎的演化以及破碎对力学特性的影响层面,对颗粒破碎演化路径的物理机制研究较少.本文基于热力学框架,采用细观力学中细观-宏观的均匀化方法推导了颗粒体系弹性能和破碎能量耗散,并在最大能量耗散的假设下,在热力学框架内,建立了理想化的无摩擦球体颗粒等向压缩过程的弹性-破碎模型,阐述了颗粒材料破碎演化路径细观热力学机制.由于模型的推导不依赖任何唯象的经验公式,因此模型中包含的参数均有明确的物理意义.模型预测与前人试验结果对比表明,材料的初始级配对弹性压缩模量和破碎应力的影响并不相同:不同分形维数级配对应的弹性体变模量存在极大值,而破碎应力却随着分形维数的增大单调递增;颗粒破碎的演化符合最大能量耗散原理,且颗粒材料的压缩曲线可以分为弹性-破碎-拟弹性3个机制不同的阶段.

关键词: 颗粒材料 ; 颗粒破碎 ; 级配 ; 热力学 ; 细观力学

Abstract

Particle breakage of granular materials is ubiquitous in nature and engineering practices and often takes place under high stress levels. The phenomenon of particle breakage may not only influence the mechanical response of granular materials, resulting the contraction of the volume of the material and reduction of the shearing strength, but is also closely associated to a variety of engineering problems. The existing research is mainly focused on depicting the evolution of the particle breakage and uses a quantifiable parameter to relate the particle breakage to the subsequent mechanical response. However, little attention has been paid to exploring the underlying physics of the driving force that initiates and attenuates the particle breakage. In this study, we present the formulation of an elastic-breakage model for the isotropic compression of frictionless spheres in the framework of thermodynamics. In the model, both the elastic strain energy and the dissipation due to particle breakage are formulated using the micro-macro averaging procedure, which is often used in micromechanics of granular materials. The evolution path of the particle breakage is determined using the maximum energy dissipation hypothesis. As the modelling does not involve any empirical results, all the model parameters have concrete physical meanings. Comparison of the model prediction with the experimental data in the literature showed that the initial gradation has different effects on the elastic bulk modulus and the breakage stress: the bulk modulus increase initially and then decrease with the fractal dimension of the gradation, which implies that there is a peak bulk modulus for a certain value of the fractal dimension; while the breakage stress increases monotonically with the increase of the fractal dimension. In addition, both the bulk modulus and the breakage stress increase monotonically with the increase of the polydispersity of the particle sizes. The evolution path of the gradation due to particle breakage is found to indeed satisfy the maximum dissipation hypothesis. Both experimental results and model prediction show that the compression curve of granular materials can be divided into three stages: the elastic compression stage under low compressive stress, particle breakage stage and the pseudoelastic compression after sufficiently large amount of particle breakage.

Keywords: granular material ; particle breakage ; gradation ; thermodynamics ; micromechanics

0

PDF (2479KB) 元数据 多维度评价 相关文章 收藏文章

本文引用格式 导出 EndNote Ris Bibtex

沈超敏, 刘斯宏. 颗粒材料破碎演化路径细观热力学机制1)[J]. 力学学报, 2019, 51(1): 16-25 https://doi.org/10.6052/0459-1879-18-340

Shen Chaomin, Liu Sihong. EVOLUTION PATH FOR THE PARTICLE BREAKAGE OF GRANULAR MATERIALS: A MICROMECHANICAL AND THERMODYNAMIC INSIGHT1)[J]. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(1): 16-25 https://doi.org/10.6052/0459-1879-18-340

颗粒材料在一定外载荷条件下会发生颗粒破碎[1].颗粒破碎与大量工程问题有关,如高土石坝的沉降变形[2]和煤块的气动输送[3]等.同时,颗粒破碎也是影响颗粒材料本构关系的重要因素,例如,砂土或粗粒料的颗粒破碎会导致孔隙比和平均应力空间内的正常固结线以及临界状态曲线发生"下移",进而影响材料的剪胀变形[4- 5].

大部分的试验和离散单元法数值模拟以颗粒材料的压缩作为探究颗粒破碎机制最简单有效的途径.针对单一颗粒的压缩试验和数值模拟结果表明,颗粒破碎应力与其尺寸负相关[6],且随着有效接触点数目的增加而提高[7].\,\,而颗粒材料作为颗粒集合体,其破碎特性更为复杂:颗粒材料的破碎与应力大小、应力状态、颗粒岩性、颗粒形状、颗粒尺寸和级配等多种因素有关,并进而导致颗粒材料的压缩曲线呈现不同的形态[7].然而,这些形态各异的压缩曲线可以分为三个共性的阶段:低围压下的弹性压缩阶段、颗粒破碎阶段和颗粒充分破碎后的拟弹性阶\段[8].离散单元法的迅速发展为颗粒破碎过程的细观机制的揭示提供了新契机.Cheng等[9]在Robertson[10]的开创性工作的基础上,采用绑定基本颗粒的团簇法,模拟了颗粒破碎过程中的能量演化,并指出了颗粒材料发生破碎伴随着弹性能、摩擦能量耗散和破碎能量耗散的变化.而后,团簇法被广泛运用于颗粒形状、级配等变量对破碎的影响机制研究[11-12].另一种有效模拟颗粒破碎的方法是定义颗粒层面的破碎准则,并用"子颗粒"代替大颗粒并继续参与计算的分裂法[13- 15].Ben-Nun等[16]采用了该方法研究了颗粒破碎的演化规律,认为颗粒破碎在接触力分布空间存在"吸引子",并最终导致颗粒破碎演化为唯一的分形分布.近年来,FDEM[17]和随机虚拟开裂DEM[18]等方法被广泛提出,用于模拟更复杂条件下的颗粒破碎.

虽然颗粒破碎的宏观响应和细观机制已经有了深入的研究,颗粒破碎过程宏细观尺度间的定量联系尚未明确.McDowell和Bolton[19-20]通过修正颗粒破碎的概率分布函数,重新推导了考虑颗粒破碎的能量守恒方程,并提出了颗粒破碎引起的材料破碎硬化模型.Einav[21-22]在热力学框架内提出了破碎力学理论,其中颗粒破碎伴随的内能与能量耗散的变化与当前级配相对于初始与最终级配的相对位置建立了定量关系.Zhang和Buscarnera[23]在破碎力学理论基础上,进一步考虑了环境变量对颗粒材料破碎演化的影响.尽管破碎力学理论能够较好地预测颗粒材料的级配演化及其对宏观力学响应的影响,理论中破碎演化路径的确定仍是基于试验结果,而没有较好的理论支撑.

考虑到颗粒破碎影响因素复杂,本文从无摩擦球形颗粒的破碎入手,在热力学与细观力学理论基础上,借鉴了破碎力学的框架,推导了宽粒径无摩擦脆性球体颗粒集合体等向压缩过程中的破碎演化,并将预测结果与试验资料对比验证,实现了破碎理论从描述颗粒破碎到预测颗粒破碎的突破.尽管严格的无摩擦颗粒在自然与工程界并不存在,但是考虑到颗粒摩擦引起的能量耗散远低于破碎断裂引起的能量耗散[21],忽略摩擦耗散是有效且合理的途径.另外,由于颗粒破碎时,球形颗粒与其他形状颗粒仅在颗粒配位数与比表面积上存在定量差异,因此本文的理论对理解更为复杂的颗粒材料的破碎行为也具有参考价值.

1 破碎表征

颗粒破碎会伴随着材料级配曲线的变化.为了表征颗粒破碎的程度,学者们提出了一些破碎指标用以表征破碎的程度,其中Hardin[24]的破碎指标$B_{r}$以及Einav[21]在此基础上的改进形式应用最普遍.这两种破碎指标可以用以下广义的形式表达

$$\varGamma = L\left( {\varGamma _0 ,\varGamma _{u} ,B_{r} }\right)(1)$$

式中,$\varGamma $是当前级配,$\varGamma _0 $和$\varGamma _{u}$分别是最初和最终级配,$B_{r}$是无量纲的相对破碎指标,$L$代表$B_{r}$的线性映射.可以看出,采用相对破碎指标表征颗粒破碎的优点显著:在已知材料的最初和最终级配的前提下,式(1)既定义了破碎率指标,同时通过线性映射$L$指定了级配演化的路径,因此仅通过一个参数$B_{r}$即可唯一确定当前级配.然而,人为强制给定级配的演化路径不仅无法解释颗粒破碎的物理内涵,也使得相应的破碎力学理论中模型参数依赖初始级配,不能预测初始级配对破碎演化的影响.

为了使级配曲线有选择演化路径的自由度,本文采用以下函数拟合任意时刻颗粒材料的级配曲\线[21]

$$\frac{M_{\left( {L < r} \right)} }{M_{T} } = \frac{\varLambda^{3-D}{-}\lambda ^{3-D}}{\varLambda ^{3-D}{ -}1}(2)$$

式中,$M_{\left( {L < r} \right)}

$是粒径小于$r$的颗粒总质量;$M_{T} $是试样总质量;$\lambda{ = }r /{r_{\max } }$为粒径与最大粒径的比值;$\varLambda$为$\lambda $的最小值,表征粒径的均一度,当$\varLambda =1$时,则代表该试样为均一粒径,当$\varLambda =0$时,则试样趋于分形分布,此时参数$D$为级配的分形维数[25-26].值得说明的是,首先式(2)并非描述任意的自然或工程界存在的颗粒材料级配的精确方法,只是一种近似的拟合;其次,式(2)是土体级配曲线中应用较为广泛的拟合函数之一,若改用其他拟合函数,本文接下来的理论推导依然可以进行.根据式(2)可知,级配曲线的演化可以通过参数$\varLambda$,$r_{\max }$和$D$的变化表征,记为以下向量形式

$$\varGamma _i = \left( {r_{\max } ,\varLambda ,D}\right), i= 1,2,3(3)$$

对于均质球形颗粒,可通过式(2)进一步推导得到颗粒数目的概率分布函数

$$p\left( r \right) = \frac{{-D}}{{{r_{\max }}\left( {1-{\varLambda ^{-D}}} \right)}}{\lambda ^{-1 -D}}(4)$$

2 颗粒破碎的细观热力学过程

2.1 热力学框架

任意一个率无关的热力学过程的能量守恒可以写作

$$\dot {W}=\dot {u}+\dot {\varPhi ,\dot {\varPhi}\ge {0}}(5)$$

其中,$W$是外界对单位体积试样做的功;$u$和$\varPhi$分别是弹性能和能量耗散密度;顶标“$\cdot$”代表该变量关于时间的偏导$\partial/{\partial t}$.对于处于等向压缩应力条件下的无摩擦球形颗粒体系,弹性能$u$和能量耗散率$\dot{\varPhi }$可以写成以下形式

$$u{ = }u\left( {\sigma _{m},\varGamma } \right),\dot {\varPhi } = \dot {\varPhi }\left({\varGamma ,\dot {\varGamma }} \right)(6)$$

式中,$\sigma _{m} $是施加的球应力,$\varGamma$为式(3)定义的级配演化指标.

颗粒破碎往往伴随着颗粒重排引起的摩擦耗散,然而由于本文考虑了无摩擦的颗粒,因此不考虑颗粒在摩擦作用下趋于稳定过程对应的变形量.同时Einav[21]也指出,相较于摩擦耗散,颗粒体系的主要能量耗散由颗粒破碎造成.在忽略摩擦耗散的条件下,总应变可以认为等于弹性应变,即$\dot{\varepsilon }_{v} { = }\dot {\varepsilon }_{v}^{e} $.因此外界输入功率$\dot {W}$可以写成

$$\dot {W} = \sigma _{m} \dot {\varepsilon }_{v} = \sigma _{m} {\dot {\sigma}_{m} } /{K_{B}^{e} }(7)$$

式中,$\varepsilon _{v} $为体应变,$K_{B}^{e}$为弹性体积模量.将式(6)和式(7)代入式(5)可得

$$\left( {{\sigma_{m} } /{K_{B}^{e} }{-}\frac{\partial u}{\partial \sigma _{m} }} \right)\dot {\sigma }_{m} =\frac{\partial u}{\partial \varGamma }\dot {\varGamma } + \dot{\varPhi }(8)$$

式(8)为颗粒材料压缩过程的一般热平衡方程,假设较低应力条件下颗粒未发生破碎,即$\dot{\varGamma }{ = }\dot {\varPhi }{ =}0$,则该过程为可逆过程,式(8)可以简化为

$${\sigma _{m} }/{K_{B}^{e} }{-}\frac{\partial u}{\partial \sigma_{m} }{ = }0(9)$$

类似地,假设某颗粒材料在某一恒定应力条件下发生破碎,即$\dot{\sigma }_{m} { = }0$,则可得到颗粒破碎引起的耗散部分

$$\frac{\partial u}{\partial \varGamma }\dot {\varGamma } + \dot {\varPhi}{ = }0(10)$$

尽管式(9)与式(10)似乎仅在附加条件下成立,但考虑到应力变化率$\dot{\sigma }_{m} $和破碎率$\dot {\varGamma}$取值独立于状态量,若要使式(8)恒成立,则需要式(9)与式(10)同时成立.根据热力学第二定律可知,能量耗散非负,即$\dot{\varPhi } \ge0$.因此,结合式(10)可以进一步得到$\dfrac{\partial u}{\partial \varGamma }\dot {\varGamma } \le0$,也就是说颗粒破碎会沿着降低系统弹性能的方向演化.

2.2 宏细观推导

本节主要基于细观颗粒尺度的推导,并通过均匀化过程建立式(8)$\sim$式(10)中的弹性能与破碎耗散的表达式.本文采用的均匀化方法是采用颗粒体积加权平均法,即对于任意颗粒尺度变量$X^{p}$,其试样尺度的平均值按照如下公式定义

$$X = \frac{\overline{X^{p}V_{p} } }{\left( {1 + e} \right)\overline {V_{p}} }(11)$$

式中,上标p代表该变量为颗粒尺度变量;$X$是试样尺度的平均值;$e$是试样的孔隙比,顶标“-”代表试样内变量的平均值,即$\dl\int_{{r}_{\min }}^{r_{\max } } {p\left( r \right)} {d}r$.

(1) 弹性能密度

准静态颗粒单一颗粒的细观应力$\sigma _{ij}^{p} $可以写作

$$\sigma _{ij}^{p} = \frac{1}{V_{p} }\sum\limits_{c \in V_{p} } {f_i^c } l_j^c(12)$$

式中,上标p代表颗粒尺度,上标$c$代表接触点,$f_i^c $和$l_j^c$分别为接触力矢量和接触支向量,$V_{p} $是颗粒体积,$c \inV_{p} $代表遍历了颗粒$V_{p}$的所有接触点.对于半径为$r$的无摩擦的球形颗粒,有$V_{p} =\dfrac4 3{\pi}r^3$,$\left| {{{l}}^c} \right| =r$,且接触力方向与支向量重合,因此式(12)可以简化为

$$\sigma _{m}^{p} = \frac{f\left( r \right)C\left( r \right)}{4{\pi}r^2}(13)$$

式中,$f$为该颗粒上的平均接触力,$C$是颗粒的配位数.

结合式(11)定义的均匀化方法和式(13),试样尺度的球应力$\sigma _{

m} $为

$$\sigma _{m} { = }\frac{\overline {C\left( r\right)f\left( r \right)r} }{4{\pi }\left( {1 + e}\right)\overline {r^3} }(14)$$

Hertz接触条件下颗粒的接触力与位移的关系\ 为[27]

$$\varDelta = \varOmega r^{-1 / 3}f^{2 / 3}(15)$$

式中,$\varDelta $为颗粒接触点的变形位移;$\varOmega$是颗粒本身的材料参数,与杨氏模量$E$和泊松比$\nu $的关系为

$$\varOmega ^3 = \frac{9}{4}\left( {\frac{1-\nu ^2}{{\pi }E}}\right)(16)$$

对于任意Hertz接触,通过对式(15)的积分可以得到该接触点中储存的弹性能

$$E_c = \int_0^\varDelta {f\left( \varDelta \right)} {d}\varDelta {= }\frac{2}{5}\varOmega r^{-1 / 3}f^{5 / 3}(17)$$

因此,对于一个有若干个接触点的颗粒,其接触点中存储的弹性能密度为

$$u^{p}\left( r \right) = \frac{{1}}{V_{p}}\sum\limits_{c \in V_{p} } {E_c } { = }\frac{2}{5V_{p} }\varOmega r^{-1 / 3}\left[ {f\left( r \right)} \right]^{5 /3}C\left( r \right)(18)$$

由式(11)和式(18)可知某一颗粒试样均匀化后的弹性能密度为

$$u{{ =}}\frac{{3\varOmega }}{{10{{\pi}}\left( {1 + e}\right)}}\frac{{\overline {{r^{-1/3}}{{\left[ {f\left( r\right)} \right]}^{5/3}}C\left( r \right)} }}{{\overline {{r^3}}}}(19)$$

值得说明的是,式(14)和式(19)建立的宏观球应力和弹性能的细观表达式不仅包含细观常数,也包含细观变量(配位数$C$和接触力$f)$.由于细观变量在宏观层面无法测量,因此还需要进一步将这些细观变量与宏观变量建立联系.

为了建立颗粒配位数$C$与当前级配的关系,本文借鉴了Shaebani等[28]针对多分散颗粒体系配位数的预测理论,并进行了适当改进.图1为某一颗粒在多分散体系中的接触示意图,其中$r$为待求配位数的颗粒半径,$R$为周围接触颗粒的粒径. $S\left( {r,R}

\right)$为颗粒$R$在颗粒$r$外表面投影的面积,不难计算得到

$$S\left({r,R} \right) = 2{\pi }r^2\left( {1-\sqrt {\frac{2Rr +r^2}{\left( {{r} + R} \right)^2}} } \right)(20)$$

考虑到周围接触颗粒的粒径并非常数,而是应当符合与颗粒体系级配相同的分布,因此式(20)对应的数学期望值$S\left(r \right)$应为

$$S(r) = 2{\pi }r^2\int_{r_{\min } }^{r_{\max }} {\left( {1-\sqrt {\frac{2Rr + r^2}{\left( {{r} + R}\right)^2}} } \right)p(R){d}R}(21)$$

颗粒配位数应与颗粒$r$的表面积除以式(21)给出的单颗粒占据的投影面积成正比,表示为

$$C(r) = \frac{4{\pi }r^2}{c_{s} S(r)} = {2}\Big/\left[{c_{s}\int_{r_{\min } }^{r_{\max } } {\left( {1-\sqrt {\frac{2Rr +r^2}{\left( {r + R} \right)^2}} } \right)p(R){d}R}}\right](22)$$

式中,参数$c_{s}$又称为颗粒体系的线性容积率[28],由于实际情况下,颗粒$r$的表面不一定被与其接触的颗粒投影完全覆盖,1/$c_{s}$表征其表面积的覆盖率.文献[28]的模拟结果表明,$c_{s}$在一定条件下可以认为是常数.

图1   多分散体系颗粒接触示意图...

Fig. 1   Contact information for a typical particle in a polydisperse granular media

式(22)建立了多分散颗粒体系内的配位数与级配的关系,因此式(19)中的未确定的细观变量只有颗粒的接触力$f$.尽管颗粒体系的接触力分布已经有了广泛的讨论并有相应的统计分布表达式[29-30],然而统计描述不可避免地会引入物理意义不明确的拟合参数,因此本文提出了另一种处理方法.本文假定某一应力条件下,颗粒体系内的$C\left(r \right)f\left( r \right)$的取值与颗粒的表面积${4\pi}r^2$成正比,且该比值与尺寸无关.这一假设的本质是认为颗粒尺度的细观应力$\sigma_{m}^{p} $为尺度无关量,即

$$\sigma _{m}^{p} =\frac{C\left( r \right)f\left( r \right)}{4{\pi }r^2} = {const}.(23)$$

由于$\sigma _{m}^{p}$为尺度无关量,因此根据式(11)的定义,宏观球应力$\sigma _{m}$可以与$\sigma _{m}^{p} $建立更简单的联系

$$\sigma _{m} = \frac{\sigma _{m}^{p} }{1 + e}(24)$$

将式(24)代入式(23)可得

$$C\left( r \right)f\left( r \right){ =4\pi }\left( {{1 + e}} \right)\sigma _{m} r^2(25)$$

通过式(25),即可将颗粒的接触力$f$与宏观应力$\sigma _{m}$和配位数$C$建立联系.

结合式(19)和式(25),则试样的弹性能密度为

$$u{ = }\left[ {{4\pi }\left( {1{ + }e} \right)} \right]^{5 /3}\frac{3\varOmega }{10{\pi }\left( {1 + e} \right)}\sigma_{m} ^{5 / 3}\frac{\overline {{r^3} /{\left[ {C\left( r\right)} \right]}^{2 / 3}} }{\overline {r^3} }(26)$$

进一步地,若将式(4)的颗粒分布代入式(26),则弹性能密度最终可以表示为

$$u = \varTheta \sigma _{m}^{{5 / 3}} \eta \left( {\varLambda ,D}\right)(27)$$

式中$\varTheta $为颗粒材料和密实程度有关的参量

$$\varTheta { =}\frac{6}{5}\sqrt[3]{4}\varOmega \left( {{\pi}c_{s} }\right)^{2 / 3}\left( {1 + e} \right)^{2 / 3}(28)$$

$\eta$ 为当前级配的函数

$$\eta \left( {\varLambda ,D} \right) =\left( {\dfrac{-D}{1-\varLambda ^{-D}}} \right)^{2 /3}\int_\varLambda ^1 \lambda ^{2-D}\Bigg[\int_{\varLambda}^1\Bigg( 1-\\ \sqrt{\frac{2{\lambda }'\lambda + \lambda ^2}{\left({\lambda + {\lambda }'} \right)}^2}\Bigg){\lambda '}^{-1 -D}{d}{\lambda }'\Bigg] ^{2 / 3}{d}\lambda\Bigg/{\int_{\varLambda} ^1 {\lambda ^{2-D}{d}\lambda }}(29)$$

(2) 破碎能量耗散密度

对于无摩擦颗粒,颗粒破碎的能量耗散主要是颗粒的断裂能损耗,根据断裂力学的基本假定,颗粒发生破碎后会新增表面积,断裂能与该表面积的增量成正比,表示为$${\dot \varPhi ^{{p}}} = {G_{f}}{\dot S^{{p}}}(30)$$式中,$\varPhi ^{p}$为单颗粒的破碎能量耗散;$S^{p}$为颗粒的单位体积具有的表面积,对于球形颗粒有$S^{p} = 3 /r$; $G_{f}$为颗粒的应变能释放率.

类似地,结合式(11)和式(30)可以得到颗粒试样尺度均匀化的断裂能量耗散密度

$$\dot {\varPhi }{ = }G_{f} \dot {S}{ = }\frac{{3}}{1{ +e}}G_{f} \frac{\partial \left( {{\overline {r^2} } /{\overline{r^3} }} \right)}{\partial t}(31)$$

进一步将式(4)的粒径分布代入式(31),则断裂能量耗散密度可以最终表达为

$$\dot {\varPhi } = \frac{G_{f} }{1 + e}\dot {S}(32)$$

其中$S$为级配的函数

$$S = 3\frac{3-D}{2-D}\frac{1-\varLambda^{2-D}}{1-\varLambda ^{3-D}}\frac{1}{r_{\max } }(33)$$

3 弹性-破碎演化

结合式(9)和式(27),则可以建立弹性阶段体变模量与宏观应力的关系

$$K_{B}^{e} = \frac{3}{5\varTheta \eta }\sigma _{m}^{{1 / 3}}(34)$$

根据定义,弹性体变模量可以进一步与孔隙比变化建立关系,式(34)可以进一步写成以下微分形式

$${-}\frac{{d}e}{1 + e} = \frac{5}{3}\varTheta \eta \sigma_{m}^{{-1 / 3}} {d}\sigma _{m}(35)$$

定义$\chi _i $为单位破碎速率$\dot {\varGamma}$引起的破碎能量耗散率;同样,定义$\bar {\chi }_i$为单位破碎指标$\varGamma $引起的内能损失,两者的定义式为

$$\chi_i { = }\frac{\partial \dot {\varPhi }}{\partial \dot{\varGamma }_i },\quad\bar {\chi }_i { = }\frac{\partial u}{\partial \varGamma _i }(36)$$

由于本文假定的颗粒破碎是率无关的过程,即认为颗粒破碎引起的耗散是破碎速率$\dot{\varGamma }$的一阶齐次式,因此根据齐次函数的欧拉定理可知

$$\frac{\partial \dot {\varPhi }}{\partial \dot {\varGamma }}\dot {\varGamma}{ = }\dot {\varPhi }(37)$$

将式(36)和式(37)代入式(10)可得

$$\left( {\chi _i + \overline{\chi}_i } \right)\dot {\varGamma }_i = 0(38)$$

由于$\chi _i $, $\overline{\chi }_i $和$\dot {\varGamma}$都是向量,因此式(e8)仅能表明向量$\left( {\chi _i +\overline{\chi }_i } \right)$和$\dot {\varGamma }_i$正交.然而,Ziegler[31]认为热力学中的最大耗散功原理允许比式(38)更严格的条件,写作

$$ {\chi _i + \overline{\chi }_i } = 0(39)$$

最大耗散能假设可以通过颗粒材料级配演化的规律验证:由于式(26)中的弹性能是尺寸无关的,即${\partial u} /{\partial r_{\max } } = 0$,因此结合式(39)可得${\partial \dot{\varPhi }} / {\partial \dot {r}_{\max } } =0$.最大耗散功假设的推论表明级配演化过程中$r_{\max }$不发生变化,这个推论与一系列针对颗粒材料压缩的试验结果吻合[8,32].事实上,这一结果并不是显而易见的,在Einav之前,主流观点认为在足够大的压力作用下,颗粒体系会破碎成粒径小于0.074mm的粉\ 末[31]. 由于颗粒破碎过程中$r_{\max }$保持不变,因此破碎指标$\varGamma _i $降维成

$$\varGamma _i =\left( {\varLambda ,D} \right)(40)$$

借鉴Collin和Houlsby[34]基于热动力学的本构建模理论,颗粒破碎速率$\dot{\varGamma }_i $和单位破碎速率引起的破碎能量耗散率$\chi _i$可以通过Legendre变换互换,而一阶齐次函数的Legendre变换对应着$\chi_i$恒为零的函数即为材料的屈服.因此,将能量耗散率表达式(32)和式(33)代入式(39),可得

$$\lambda _{B} y_{B} \equiv F_i \chi _i-G_{f} = 0(41)$$

式中,为方便起见,定义${F_i} = {{\left( {1 + e} \right)} /{\left({\partial S/\partial {\varGamma _i}} \right)}}$为非负的乘子.

式(41)的Legendre变换可以得到流动法则

$$\dot {\varGamma }_i =\lambda _{B} \frac{\partial y_{B} }{\partial \chi _i } =\lambda _{B} F_i(42)$$

对式(41)微分可得

$$F_i \dot {\chi }_i + \frac{\partial F_i}{\partial \varGamma _j }\dot {\varGamma }_j \chi _i = 0(43)$$

同时,结合式(36)和式(39)可得

$$\chi _i { =-}\bar {\chi }_i{ =-}\frac{\partial u}{\partial \varGamma _i}(44)$$

由于弹性能为应力和级配的函数,因此式(44)的进一步全微分可得

$$\dot {\chi }_i =-\frac{\partial ^2u}{\partial\varGamma _i \partial \varGamma _j }\dot {\varGamma }_j -\frac{\partial ^2u}{\partial \varGamma _i \partial \sigma _{m}}\dot {\sigma }_{m}(45)$$

最终,将式(42)和式(45)代入式(43)可得流动法则中非负乘子为

$$\lambda _{B} =-\dfrac{F_i \dfrac{\partial^2u}{\partial \varGamma _i\partial \sigma _{m} }\dot {\sigma }_{m} }{F_i\dfrac{\partial ^2u}{\partial \varGamma _i \partial \varGamma _j}F_j + \dfrac{\partial u}{\partial \varGamma _i }\dfrac{\partia lF_i }{\partial B_j }F_j }(46)$$

4 模型预测与讨论

本文提出的颗粒破碎模型包含4个物理概念明确的参数,包括3个与单颗粒的材料相关参数:杨氏模量$E$,泊松比$\nu$和应变能释放率$G_{f}$,以及一个颗粒体系组构相关参数:线性容积率$c_{s}$.由于该模型的推导没有基于任何唯象的试验或者数值模拟结果,因此颗粒破碎演化的模拟是完全基于细观力学和热力学的预测.

4.1 弹性压缩

式(34)预测认为颗粒体系的弹性体变模量与应力的1/3次方成正比,即$K_{B} \sim \sigma _{m}^{{1 / 3}}$.图2对比了模型预测结果与两种颗粒材料[35]的一维压缩模量随着竖向压力$\sigma'$的变化规律,其中压缩模量$K_{B}$和压力$\sigma'$均除以了大气压力$P_{a}$.可以看出,两种颗粒材料的压缩模量随着施加压力的增大,在双对数坐标内几乎线性增大,其斜率与本文预测的结果吻合较好,直到应力超过6MPa后出现显著拐点(颗粒开始大量破碎).

图2   弹性体变模量的应力相关性...

Fig.2   Pressure-dependent elastic bulk modulus

同时,结合式(29)和式(34),可以得出级配曲线(通过参数$\varLambda$和$D$表征)对弹性体变模量的影响.图3($a)$给出了模型预测的弹性体变模量随着参数$\varLambda$和$D$的演变趋势,其中为了消除材料参数和球应力的影响,图中压缩模量$K_{B}$乘以了系数$\varTheta / {\sigma _{m}^{{1 / 3}}}$.模型预测结果表明,随着粒径均一度的降低($\varLambda$减小),材料的弹性模量提高;随着级配分形维数$D$的提高,弹性模量先增大后减小.压缩模量随着级配分形维数先增后减的趋势与Minh和Cheng[36]的DEM模拟结果大致相同(见图3(b)),其中,DEM模拟采用土力学中常用的压缩指标$C_{c}$的倒数定性表征压缩模量.

图3   级配曲线对弹性体变模量的影响...

Fig.3   Influence of the gradation on the elasticcompressibilityof granular materials

值得说明的是,理论预测的最大压缩模量对应的分形维数$D\approx 1.65$小于DEM模拟结果$D \approx2.1$,原因一方面是DEM模拟中采用了与本文不同的级配定义函数,另一方面可能由于本文的理论模型中对配位数的预计方法的精确程度还有待进一步提高.

4.2 破碎应力

将式(27)、式(33)和式(36)代入屈服函数(41),可以得到破碎应力$\sigma_f $

$$\sigma _{f} = \left( {-\dfrac{G_{f} }{\varThetaH_i \dfrac{\partial \eta }{\partial \var Gamma _i }}} \right)^{3 /5}\left( {r_{\max } } \right)^{-3 / 5}(47)$$

式中,记$F\left( {r_{\max } ,\varLambda ,D} \right) = H\left({\varLambda ,D} \right)r_{\max }$.从式中可以看出,颗粒破碎应力与颗粒本身的应变能释放率$G_{f}$正相关,与材料参数$\varTheta$负相关.级配对破碎应力的影响主要通过$H_i $,$\eta$和$r_{\max}$表征:首先,式(47)表明破碎应力与试样的尺寸有关,且正比于$\sigma_{f} \propto r_{\max }^{-3 / 5}$.图4对比了大量的单颗粒的硅砂的破碎应力和模型的预测结果,可以看出,试验结果在双对数坐标系内可以用一条直线拟合,且斜率与理论预测的值基本吻合.

图4   破碎应力的尺寸效应...

Fig.4   Influence of the particle size on the breakage stress

进一步地,图5给出了根据式(47)预测的级配对颗粒破碎应力的影响,其中破碎应力$\sigma_{f} $乘以了系数$\left( {{r_{\max } \varTheta } /{G_{f} }}\right)^{3 /5}$排除材料参数和颗粒尺寸的影响.可以从图中看出,与压缩模量不同,破碎应力随着级配的分形维数$D$增大而增大,且随着粒径均一度$\varLambda$的降低而增大.且当材料趋于分形分布时($\varLambda \to0)$,破碎应力大幅提高,材料后续破碎的难度随之增大,这也解释了颗粒材料在充分破碎后级配会止步于分形分布的普遍现象.

图5   预测的级配对破碎应力的影响规律...

Fig.5   Prediction of the effect of the gradation on thebreakagestress

4.3 级配演化

本文提出的颗粒破碎模型可以根据破碎屈服和流动准则预测颗粒破碎过程中级配的演化.图6给出了Hagerty等[8]针对BlackBeauty矿渣料的一维压缩试验得到的不同竖向压力条件下级配曲线.在预测模型中,本文假定了压缩的侧限系数(水平向应力除以竖向应力比值)$K_0=0.4$.模型参数上,由于岩石材料泊松比变化幅度不大,且其取值对模型影响较小,本文酌情取$\nu{ = }0.09$;参数$c_{s}$通过宽级配料的DEM模拟结果可得$c_{s} \approx 3.0^{[28].根据式(34)的体变模量公式,结合Hagerty试验中材料的初始级配与其对应的体变模量$K_{B}=$237 MPa,可以计算得到$E \approx 100{GPa}$;根据式(47)破碎应力公式,结合Hagerty试验得到的材料破碎应力$\sigma_{f} { = }18{MPa}$以及该应力条件下的材料级配,可以计算得到应变能释放率$G_{f}$取值约为$300{J} /{m}^{2}$.

根据以上参数,图6对比了模型预测的级配曲线演变和一维压缩试验过程中级配曲线演变.从图中可以看出,在不预设级配演化路径的前提下,仅通过最大耗散能假设,预测得到的级配曲线演化与试验结果吻合度较高.图6中插图为级配分形维数$D$的演变规律,可以看出,试验和模型预测都表明,随着颗粒破碎的继续,级配的分形维数增大的趋势逐渐趋缓.

图7给出了模型预测的压缩曲线($e-\sigma _{m}$曲线),由于本文提出的模型未考虑塑性体应变,因此模型预测结果无法直接与试验结果对比.不过,总体上可以看出,颗粒体材料的压缩存在3个阶段,在低应力条件下的弹性压缩阶段、颗粒破碎阶段和高围压下颗粒破碎基本停止的拟弹性阶段.这一现象与大量试验结果[8,32]吻合,事实上,颗粒破碎会增大材料体应变,导致压缩曲线出现拐点,而当级配接近分形分布时,颗粒破碎应力提高,颗粒破碎量降低,破碎变形回归拟弹性阶段.

图6   颗粒破碎过程中级配的演化路径...

Fig.6   Evolution of the gradation due to particle breakage

图7   模型预测的压缩曲线...

Fig.7   Predicted compression curve

5 结论

本文基于热力学和细观力学理论,建立了理想化的无摩擦球体颗粒等向压缩过程的弹性-破碎模型,探索了颗粒材料破碎演化路径细观热力学机制.与破碎力学理论采用相对破碎指标表征破碎程度和破碎演化路径不同,本文采用级配曲线的特征参数的变化(最大粒径$r_{\max} $,分形维数$D$和粒径均一度$\varLambda)$表征颗粒破碎率,颗粒破碎的演化路径则基于最大耗散功假设推导得到.由于模型推导中没有采用任何试验或模拟资料拟合,本文提出的弹性-破碎模型包含的四个参数(颗粒的杨氏模量$E$,泊松比$\nu$和应变能释放率$G_{f}$以及颗粒系统的线性容积率$c_{s})$均有明确的细观物理意义.

采用该模型预测颗粒的弹性体变模量、屈服应力和破碎演化趋势,并与现有的试验和数值模拟结果对比,得到了以下与试验及模拟结果吻合的现象:在弹性压缩阶段,颗粒材料的体变模量与球应力的1/3次方成正比;弹性体变模量随着粒径均一度的减小而增大,随着级配分形维数的增大先增大后减小,存在某一分形维数对应最大的弹性压缩模量;颗粒材料的破碎应力与颗粒试样的尺寸有关,且与最大粒径的$-$3/5次方成正比;破碎应力随着分形维数的增大而增大,且随着粒径均一度的减小而增大;最终,理论预测的级配演化规律与试验结果吻合度较高.

另外,本文的理论模型也揭示了颗粒材料破碎的机理:当级配接近分形分布时,破碎应力的大幅增加是"颗粒材料在高应力条件下破碎趋于停止"的主要原因;在没有预先设定破碎演化路径的前提下,本文模型预测的颗粒破碎演化路径与试验结果吻合较好,表明了颗粒破碎过程确实沿着能量耗散最大化的方向发展;尽管因为没有考虑塑性体应变,从本模型预测的压缩曲线上也可以看出颗粒材料的压缩分为弹性-破碎-拟弹性3个不同阶段.

由于本文的关于颗粒材料的推导是基于无摩擦球体的假设,因此未能考虑塑性体应变和颗粒形状的影响;针对本文提出的破碎模型进一步的DEM和室内试验验证还有待进一步开展.另外,由于热力学框架的通用性,在解决应力-组构定量关系问题和应变局部化能量分布问题的基础上,本文提出的理论也可以拓展到考虑任意应力路径的颗粒破碎.

The authors have declared that no competing interests exist.


参考文献

[1] 王增会, 李锡夔.

基于介观力学信息的颗粒材料损伤--愈合与塑性宏观表征

. 力学学报, 2018, 50(2): 284-296

URL      [本文引用: 1]     

(Wang Zenghui, Li Xikui.

Meso-mechanically informed macroscopic characterization of damage-healing-plasticity for granular materials

. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(2): 284-296 (in Chinese))

URL      [本文引用: 1]     

[2] Fu Z, Chen S, Wang T.

Predicting the earthquake-induced permanent deformation of concrete face rockfill dams using the strain-potential concept in the finite-element method

. International Journal of Geomechanics, 2017, 17(11): 04017100

[本文引用: 1]     

[3] Zhou JW, Liu Y, Du CL, et al.

Effect of the particle shape and swirling intensity on the breakage of lump coal particle in pneumatic conveying

. Powder Technology, 2017, 317: 438-48

[本文引用: 2]     

[4] 姚仰平, 张民生, 万征.

基于临界状态的砂土本构模型研究

. 力学学报, 2018, 50(3): 589-598

URL      [本文引用: 1]     

(Yao Yangping, Zhang Minsheng, Wan Zheng, et al.

Constitutive model for sand based on the critical state

. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(3): 589-598 (in Chinese))

URL      [本文引用: 1]     

[5] 万征, 孟达.

复杂加载条件下的砂土本构模型

. 力学学报, 2018, 50(4): 929-948

URL      [本文引用: 1]     

(Wan Zheng, Meng Da.

A constitutive model for sand under complex loading conditions

. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(4): 929-948(in Chinese))

URL      [本文引用: 1]     

[6] Nakata A, Hyde M, Hyodo H, et al.

A probabilistic approach to sand particle crushing in the triaxial test

. Géotechnique, 1999, 49(5): 567-583

[本文引用: 1]     

[7] Zheng W, Tannant DD.

Grain breakage criteria for discrete element models of sand crushing under one-dimensional compression

. Computers and Geotechnics, 2018, 95: 231-239

[本文引用: 2]     

[8] Hagerty M, Hite D, Ullrich C, et al.

One-dimensional high-pressure compression of granular media

. Journal of Geotechnical Engineering, 1993, 119(1): 1-18

[本文引用: 4]     

[9] Cheng Y, Bolton M, Nakata Y.

Grain crushing and critical states observed in DEM simulations

. Powders & Grains, 2005, 2: 1393-1397

[本文引用: 1]     

[10] Robertson D.

Computer simulations of crushable aggregates. [PhD Thesis]

. Cambridge: University of Cambridge, 2000

[本文引用: 1]     

[11] Liu Y, Liu H, Mao H.

DEM investigation of the effect of intermediate principle stress on particle breakage of granular materials

. Computers and Geotechnics, 2017, 84: 58-67

[本文引用: 1]     

[12] Wang J, Yan H.

DEM analysis of energy dissipation in crushable soils

. Soils and Foundations, 2012, 52(2): 644-657

[本文引用: 1]     

[13] Ciantia M, Arroyo Alvarez De Toledo M, Calvetti F, et al.

An approach to enhance efficiency of DEM modelling of soils with crushable grains

. Géotechnique, 2015, 65(2): 91-110

[本文引用: 1]     

[14] Ciantia MO, Arroyo M, Calvetti F, et al.

A numerical investigation of the incremental behavior of crushable granular soils

. International Journal for Numerical and Analytical Methods in Geomechanics, 2016, 40(13): 1773-1798

[15] Einav I.

Fracture propagation in brittle granular matter

. Proceedings of the Royal Society A: Mathematical,Physical and Engineering Sciences, 2007, 463(2087): 3021-3035

[本文引用: 1]     

[16] Ben-Nun O, Einav I, Tordesillas A.

Force attractor in confined comminution of granular materials

. Physical Review Letters, 2010, 104(10): 108001

[本文引用: 1]     

[17] Ma G, Zhou W, Regueiro RA, et al.

Modeling the fragmentation of rock grains using computed tomography and combined FDEM

. Powder Technology, 2017, 308: 388-397

[本文引用: 1]     

[18] Zhou M, Song E.

A random virtual crack DEM model for creep behavior of rockfill based on the subcritical crack propagation theory

. Acta Geotechnica, 2016, 11(4): 827-847

[本文引用: 1]     

[19] Mcdowell G, Bolton M.

On the micromechanics of crushable aggregates

. Géotechnique, 1998, 48(5): 667-679

[本文引用: 1]     

[20] Mcdowell G, Bolton M, Robertson D.

The fractal crushing of granular materials

. Journal of the Mechanics and Physics of Solids, 1996, 44(12): 2079-2101

[本文引用: 1]     

[21] Einav I.

Breakage mechanics-Part I: Theory

. Journal of the Mechanics and Physics of Solids, 2007, 55(6): 1274-1297

[本文引用: 5]     

[22] Einav I.

Breakage mechanics-Part II: Modelling granular materials

. Journal of the Mechanics and Physics of Solids, 2007, 55(6): 1298-1320

[本文引用: 1]     

[23] Zhang Y, Buscarnera G.

Breakage mechanics for granular materials in surface-reactive environments

. Journal of the Mechanics and Physics of Solids, 2018, 112: 89-108

[本文引用: 1]     

[24] Hardin BO.

Crushing of soil particles

. Journal of geotechnical engineering, 1985, 111(10): 1177-1192

[本文引用: 1]     

[25] Zhang X, Hu W, Scaringi G, et al.

Particle shape factors and fractal dimension after large shear strains in carbonate sand

. Géotechnique Letters, 2018, 8(1): 73-79

[本文引用: 1]     

[26] Tyler SW, Wheatcraft SW.

Fractal scaling of soil particle-size distributions: analysis and limitations

. Soil Science Society of America Journal, 1992, 56(2): 362-369

[本文引用: 1]     

[27] Johnson K, Kendall K, Roberts A.

Surface energy and the contact of elastic solids

. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 1971, 324(1558): 301-313

[本文引用: 1]     

[28] Shaebani MR, Madadi M, Luding S, et al.

Influence of polydispersity on micromechanics of granular materials

. Physical Review E, 2012, 85(1): 011301

[本文引用: 4]     

[29] Sufian A, Russell A, Whittle A.

Anisotropy of contact networks in granular media and its influence on mobilised internal friction

. Géotechnique, 2017, 67(12): 1067-1080

[本文引用: 1]     

[30] Mueth DM, Jaeger HM, Nagel SR.

Force distribution in a granular medium

. Physical Review E, 1998, 57(3): 3164

[本文引用: 1]     

[31] Ziegler H.

An Introduction to Thermomechanics

. Elsevier, 2012

[本文引用: 2]     

[32] Mesri G, Vardhanabhuti B.

Compression of granular materials

. Canadian Geotechnical Journal, 2009, 46(4): 369-392

[本文引用: 2]     

[33] 尹振宇, 许强, 胡伟.

考虑颗粒破碎效应的粒状材料本构研究: 进展及发展

. 岩土工程学报, 2012, 34(12): 2170-2180

URL     

(Yin Zhenyu, Xu Qiang, Hu Wei.

Constitutive relations for granular materials considering particle crushing: Review and development

. Chinese Journal of Geotechnical Engineering, 2012, 34(12): 2170-2180(in Chinese))

URL     

[34] Collins IF, Houlsby GT.

Application of thermomechanical principles to the modelling of geotechnical materials

. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 1997, 453(1964): 1975-2001

[本文引用: 1]     

[35] Pestana JM, Whittle A.

Compression model for cohesionless soils

. Géotechnique, 1995, 45(4): 611-632

[本文引用: 1]     

[36] Minh N, Cheng Y.

A DEM investigation of the effect of particle-size distribution on one-dimensional compression

. Géotechnique, 2013, 63(1): 44-53

[本文引用: 1]     

/