颗粒材料破碎行为数值分析方法研究综述1)

梁绍敏 冯云田 赵婷婷 王志华

* (太原理工大学机械与运载工程学院,太原 030024)

† (斯旺西大学辛克维奇计算工程中心,英国斯旺西 SA1 8EP)

引言

颗粒材料在自然界和工程领域普遍存在,如谷物、砂石和土壤等.在水利、港口和交通等工程领域中,堆石、砂石和砂土等颗粒材料得到了广泛应用,在堤坝、隧道和地基等工程中,大载荷作用下颗粒材料容易发生破碎现象,且颗粒材料的破碎行为会给工程和建筑建设带来极大的影响,如三峡水利枢纽工程中的填筑材料,其破碎率可能会达20%.因此,颗粒材料的破碎行为已经引起人们的关注,并且逐渐成为颗粒材料力学性质研究的新课题[1].此外,颗粒的破碎行为会引起颗粒系统级配的改变,从而影响其物理力学性质,大颗粒破碎后产生的小颗粒填充到颗粒系统孔隙中也会导致系统孔隙减小从而引起变形[2].同时,分析细观尺度下颗粒的破碎现象与宏观尺度下材料的非线性行为之间的关系,对颗粒破碎的宏、细观机理揭示具有重要意义[3].由此看来分析颗粒破碎过程既有实际的工程意义,又有理论的研究价值.

颗粒破碎是指颗粒材料在外载荷作用下发生开裂、崩解或磨耗,从而导致级配改变的现象,且受加载方式、外载荷大小以及颗粒形状、尺寸、自身强度等影响[4].颗粒材料的破碎行为会导致颗粒形状、尺寸、强度、配位数、级配以及密实度等发生变化,甚至引起颗粒系统运动状态的变化.颗粒的破碎过程存在以下3 个特征状态: 起始状态、临界状态和终止状态.把颗粒开始涉及破碎的初始状态定义为起始状态,一般破碎的起始状态可人为地合理指定,所以具有相对性;颗粒材料的应力和应变不再变化,而偏应变继续发生变化的状态称为临界状态,临界状态是客观的;当颗粒材料内部不再有破碎现象发生的极端状态称为终止状态,所以终止状态具有唯一性.颗粒破碎具有不可逆性、累加性、应力路径相关性和卸载不产生颗粒破碎的特点[4].

目前,对颗粒破碎的研究主要采用试验的方法,研究内容主要有破碎现象的描述、评价指标的建立、破碎影响因素分析以及破碎对颗粒力学特性的影响等内容.国内外已经开展了大量针对颗粒破碎的试验研究[5-6],并取得了一定的成果,如建立了颗粒破碎的评价指标,讨论了影响颗粒破碎的因素,并在试验结果的基础上建立了颗粒弹塑性本构模型[7].典型的工作有通过改变级配、相对密度、围压等条件对颗粒材料进行大型三轴剪切试验,分析其强度、变形及剪胀特性[8].

颗粒破碎过程是一个长期的、渐进的过程,而实验室试验仅仅可以模拟导致颗粒破碎的自然条件或漫长加载历史中的某一阶段[9],无法完全还原自然界中真实颗粒的破碎过程.此外,由于试验过程中无法满足试验材料和作用载荷的同一性,因此难以保证多次试验初始条件的一致性,导致结果的重复性较低,这样获得的研究结果对研究真实破碎过程意义不大.即使可以保证起始应力状态的一致性,不同加载历史或路径也将导致颗粒破碎后的大小以及破碎过程的能量变化不同.此外,试验研究还有成本高、耗时费力、难以揭示破碎的细观机理并存在一定的危险性等缺点.数值方法的快速发展和计算机技术的飞速提高,使得通过数值方法研究颗粒破碎成为可能.与实验室试验相比,数值模拟更容易获得物理参数或现象,如裂纹的萌生、扩展等,因此,数值模拟方法得到发展和推广.当然,实验室试验对验证理论和数值方法的有效性和准确性方面发挥了很重要的作用.

本文首先整理了目前研究颗粒破碎使用较多的数值分析方法,包括黏结−破碎法、碎片替代法、比例边界有限元法、组合有限元−离散元法和近场动力学法,包括方法的提出、发展、实现过程以及方法的优势和存在的问题.然后从目前的主要研究成果着手,详细讲述了破碎过程的细观机理研究和颗粒材料破碎在工程中的应用.最后对目前有关颗粒破碎的研究进行总结,并对未来的发展趋势进行展望.

1 基于离散元法的颗粒破碎数值算法

采用数值方法研究颗粒破碎过程时,考虑颗粒的离散特性,基于离散元方法理论则是最为合理和有效的方法.这一节将主要介绍基于离散元方法原理的黏结−破碎法和碎片替代法.

1.1 黏结−破碎法

1.1.1 黏结−破碎法的提出与发展

2000年Robertson[10]首次提出颗粒的黏结模型(bonded particle model,BPM),详细介绍了如何将单颗粒黏结成组合颗粒.2004 年Potyondy 等[11]对黏结颗粒通过二维和三维离散元软件PFC 对岩石的破碎过程进行了仿真,随后该方法被应用于脆性颗粒材料破碎过程的研究中.黏结−破碎模型将一定数量的基础颗粒通过黏结键形成颗粒簇,当基础颗粒间的黏结作用失效时,颗粒发生破碎行为.

黏结−破碎法在一定程度上推进了颗粒破碎的研究,并在岩土工程中解决了一些问题.如基于黏结−破碎法,构建土壤的颗粒簇,并通过压缩试验研究了颗粒材料破碎对土壤临界状态的影响[12];构造砂土的颗粒簇,对砂土颗粒在高应力下的破碎效应进行模拟,通过改变颗粒的细观参数、外载荷的加载速率、颗粒簇的数量以及黏结颗粒的平均粒径等分析了破碎现象的影响因素[13];构造堆石料的颗粒簇,通过双轴剪切数值试验,再次肯定了该方法模拟真实颗粒材料破碎特性的有效性[14];构造不规则形状的堆石颗粒簇,用以模拟堆石颗粒在破碎过程中的尺寸效应[15];采用球形单元黏结成不规则的碎石料颗粒簇,改变应力加载路径分析颗粒材料破碎对碎石料应力–应变关系及体应变特性的影响,再次证明了离散元方法能够较为准确地获得碎石颗粒的力学特征[16].Lu 等[17-18]采用不可破碎的“超级颗粒”构造真实形态颗粒,将一定数量的颗粒与“超级颗粒”通过平行黏结键连接,通过判断平行黏结失效与否来确定小颗粒与“超级颗粒”的合离,从而实现颗粒破碎过程模拟,并将其用于对道砟颗粒棱角断裂过程的分析,破碎过程如图1 所示.以上研究再现了颗粒材料破碎现象,对颗粒破碎过程的认识有一定的推动作用.

图1 黏结颗粒构造示意图[17-18]Fig.1 Schematic diagram of bonded particles[17-18]

以上研究主要是采用球体作为基础单元通过黏结键构造颗粒簇,为了更接近真实颗粒形态,以正四面体为核通过晶胞繁衍的方法生成堆石颗粒,对三轴剪切作用下堆石料的破碎过程进行研究[19].针对复杂颗粒形状,湘潭大学龙志林团队[20]采用凸多面体块作为基本颗粒,通过可逆函数法对颗粒破碎强度变异性进行建模,如图2 所示.此外,采用传统的修正“巴西”准则作为破碎准则,再根据颗粒的接触点和质心可以确定颗粒内部最终断裂,一旦目标颗粒满足断裂准则,它就会被一系列与最终断裂相一致的虚拟切割面切割成多个碎片.这样可以避免局部应力的产生及质量和体积的不守恒.这种方法不必预定义碎片的模式[20-23].该研究将黏结基础单元由球体发展到多面体,这使黏结颗粒的几何形状与真实颗粒更加接近.

图2 破碎过程的实现[20]Fig.2 Realization of breakage process[20]

目前国内外已经发展了多种离散元模拟软件,其中采用黏结−破碎法实现颗粒的破碎分析的主要有PFC、大连理工大学季顺迎课题组研发的SDEM、南京大学刘春课题组的MatDEM 等.在处理多颗粒破碎问题时,EDEM 也是基于黏结−破碎模型实现的.

1.1.2 黏结−破碎法对颗粒破碎细观力学机理的研究

Griffith 强度理论认为固体材料中包含众多细小裂纹,在外力作用下这些裂纹的尖端会出现应力集中现象,且这个应力始终为拉应力,称之为诱导应力,当这个应力大于材料的抗拉强度时,尖端的裂纹将发生扩展并导致材料最终的破坏[3],这就是材料破碎机理研究的起源.借鉴材料的破坏机理研究方法,考虑颗粒材料的破碎机理,需要讨论破碎产生的原因,包括引起破碎的外因、内因,破碎发生的初始状态、发展过程以及最终的状态.同时,要对破碎进行描述和衡量,还需要提出合理的破碎指标[24].颗粒的破碎机理直接决定破碎准则的建立,对破碎过程的模拟具有决定性作用[25].

通过关注颗粒内部接触力的变化情况以及破碎发生的位置,从而获得颗粒破碎的演化规律,进而揭示颗粒破碎的细观力学机理.Cil 等[26]采用数值方法记录了黏结颗粒内裂纹的产生和扩展过程,并与X-ray 扫描结果进行对比,发现模拟结果与扫描结果基本一致,如图3 所示,这在一定程度上将颗粒材料破碎的现象研究发展到破碎过程的机理研究.Wang等[27-28]研究了颗粒材料破碎与颗粒材料剪切破坏行为之间的关系,如图4 所示,通过微观力学分析和机理演示,深入研究了平面应变剪切条件下颗粒破碎对致密可压试样的应力比、体应变、塑性变形和剪切破坏行为的影响.首先,基于单个土壤颗粒真实断裂行为对颗粒试样进行建模,通过对比压板压缩模拟结果与物理实验室试验结果证明了模型的有效性,如图5 所示;然后研究了影响颗粒材料破碎的主要因素,结果表明,体积膨胀率和峰值应力比、塑性变形机制和剪切破坏模式随土壤可压性的变化是影响颗粒材料破碎的主要因素;最后,对破碎的机理进行了讨论,从宏观和微观力学现象出发,提出剪切带和大规模体积收缩是致密试样的两种端部破坏模式,且这两种模式分别由颗粒重排引起的膨胀和颗粒材料破碎引起的压缩控制,在土壤可压性和围压的中等范围内,两种破坏模式是组合和竞争关系,如图6 所示[27].

图3 单颗粒、三颗粒柱的SMT 和DEM 图像以及不同载荷水平下实验室规模的压缩试验[26]Fig.3 SMT and DEM images of single-particle and three-particle columns as well as laboratory-scale compression tests under different load levels[26]

图4 具有膜型柔性侧边界的DEM 试件[27]Fig.4 DEM specimens with membrane-type flexible side boundaries[27]

图5 单颗粒压裂典型模拟结果[27]Fig.5 Typical simulation results of single particle fracturing[27]

图6 黏结单元及其在20%轴向应变下的颗粒位移和旋转场[27]Fig.6 Bonding element and particle displacement and rotation field under 20% axial strain[27]

Wang 等[29]建立了水泥的椭球颗粒模型,采用黏结−破碎法分析了单颗粒的破碎机理,如图7 所示.研究结果表明,高强度粒子碎裂成几片,而低强度粒子碎裂成两大块.这说明颗粒的断裂形态随脆性的不同而变化,导致高强度颗粒出现脆性损伤,低强度颗粒出现延性损伤.对于高抗压强度材料,由于颗粒脆性,损伤后无法承受载荷.对于低抗压强度材料,由于颗粒的延展性,颗粒破碎不会立即完成,这些颗粒在损伤后仍能承受载荷.颗粒在剪应力作用下膨胀,而孔隙率同时增加.这种剪应力下的膨胀本质是一种不稳定的高势能态.颗粒的剪胀变形增大导致其结构不稳定,从而可能引起局部骨架结构的坍塌.剪胀变形的积累是造成水泥椭球颗粒局部骨架结构坍塌的原因之一.因此,在研究颗粒的破碎机理时,可通过单个颗粒的细观和宏观力学特性进行研究.

图7 椭球体颗粒整体破碎过程演化(俯视图,颗粒间线代表接触黏结)[29]Fig.7 Evolution of ellipsoidal particle disintegration process as a whole(top view,interparticle lines represent contact bonding)[29]

1.1.3 黏结−破碎法在工程中的应用

由于黏结颗粒法可以构造复杂的颗粒形态,在实现颗粒破碎的模拟中可以对颗粒强度进行表征,因此,在分析此类工程问题时可以采用该方法.

铁路道砟一般采用离散元方法建立数值模型,在列车的作用下会发生破碎.Lim 等[30]采用黏结颗粒法构造道砟的三维非规则几何形态颗粒,并进行了单颗粒压碎试验和侧限压缩试验,研究结果与物理实验进行对比从而证明模拟结果的准确性,如图8 所示.

图8 加载前黏结道砟颗粒的分布情况[30]Fig.8 Distribution of bonded ballast particles before loading[30]

对海冰与海洋结构相互作用的研究,一般采用离散元方法建立海冰的数值模型,在海洋结构的作用下海冰会发生破碎,因此,海洋工程中海冰的破碎是颗粒材料破碎的重要应用领域.Lubbad 等[31]针对冰船相互作用过程,考虑了浮冰破碎后形成的小的浮冰大多会被推到一边、旋转或淹没的复杂情况,建立了一个实时模拟船冰相互作用过程的数值模型,并将新的解析闭合解用于表示破冰过程.利用PhysX 求解计算域内所有浮冰6 自由度刚体运动方程.图9 为锥体结构的破冰过程.其中浮冰采用多面体黏结单元模拟,破碎后的黏结单元不发生二次破碎,但也不会被删除,将参与到离散元计算中并对结构产生作用.针对海洋工程结构物与海冰作用下产生的海冰破碎过程,季顺迎团队发展了球体[32]和多面体[33]离散元方法,建立了海冰的黏结破碎模型,并将其应用于海洋工程的破冰过程中,如图10所示.

图9 锥体结构的破冰过程[31]Fig.9 Ice breakage process of cone structure[31]

图10 平整冰与海洋工程结构物互相作用并发生破碎Fig.10 Interaction and fragmentation between flat ice and marine engineering structures

此外,徐永福[34]采用离散元软件PFC-2D 通过黏结单元法构造了粗粒土模型,并讨论了直剪试验中基础颗粒间的黏结力、黏结颗粒的孔隙率和粗粒土试样的孔隙率等因素对颗粒材料体应变和剪切强度的影响.

1.1.4 黏结−破碎法的优势及存在的问题

黏结颗粒法能够构造复杂的颗粒形状[15],并实现对颗粒破碎行为的力学模拟,再现了可破碎颗粒复杂的力学响应[16],可以对颗粒强度的尺寸效应和Weibull 分布特性进行一定程度的分析[35].然而,由于每个颗粒簇是由多个基础颗粒黏结而成,模拟颗粒破碎过程将涉及大量的基础颗粒从而影响计算效率.因此,基于黏结−破碎法的颗粒破碎分析方法在单颗粒或小规模颗粒集合体的破碎研究中具有一定的优势.黏结−破碎法中破碎后的颗粒无法发生二次破碎,这与真实的破碎过程不符.此外,在离散元模拟中,黏结颗粒的物理参数需要通过大量的试验进行校正[14],在模拟颗粒破碎时需要预先设定破碎路径,为了与真实破碎现象接近,需要大量的基础颗粒数[36].因此,黏结−破碎法面临的主要问题是计算效率问题,且该问题是无法通过方法优化解决的,目前主要解决办法是发展GPU 并行算法,通过提高计算机硬件来提高计算效率.

1.2 碎片替代法

1.2.1 碎片替代法的提出与发展

碎片替代法(fragment replacement method,FRM)是strm 等[37]在1998 年提出的.该方法主要实现过程如下: 当颗粒的受力情况满足预设的破碎准则时,开始碎片替换,在将要发生破碎的颗粒所占空间中生成碎片并删除原颗粒从而完成替换,如图11 所示.替换后的子颗粒在非常短的时间步内以较低的速度膨胀,这可能导致局部应力增大.因此,需要控制子颗粒的膨胀量,这主要通过局部颗粒重叠量和颗粒速度来判定[38-39].

图11 破碎替代法的破碎模式[39]Fig.11 Breakage mode of breakage substitution method[39]

碎片替代法研究颗粒破碎时,常用的破坏准则是八面体剪应力破坏准则.该准则用八面体剪应力来表示

式中,σ1,σ2,σ3分别为第1、第2、和第3 主应力.离散元方法中一个颗粒的应力张量 σij可定义为

式中,V为颗粒的体积,nc为该颗粒接触的总数,为接触力,为接触中心的支向量,σ1,σ2和 σ3则分别对应于 σ11,σ22,σ33.在模拟过程中,当离散元计算得到的八面体剪应力超过设定的容许八面体剪应力q时,则认为颗粒破碎.De Bono 等[25]指出q=0.9σf=0.9Ff/d2,其中,σf为 颗粒的破碎强度,Ff为单颗粒径向压缩试验中颗粒的峰值破碎力,d为颗粒的粒径.根据Weibull 分布理论,每个颗粒的破碎强度 σf可表示为

式中,Ps(d),m,σ0,d0分别为颗粒存活概率、颗粒破碎强度的Weibull 模量、颗粒特征破碎强度和颗粒特征破碎强度对应的颗粒粒径.其中颗粒破碎强度的尺寸效应通过 (d/d0)−3/m体现.

Åström 等[37]提出了3 个构建碎片替代模式的标准: 尽量减少碎片数量,这是为了避免多个破碎发生后,碎片数目的迅速增加导致计算成本增大;合理布置碎片颗粒,以达到破碎时局部应力突减的效果;替代破碎颗粒的碎片的粒径分布应尽可能与实际情况相符.Åström 指出这3 个准则很难同时满足.

Tsoungui 等[40]建立的碎片替换模式包含4 种尺寸12 个子颗粒,并采用二维离散元数值方法对侧限压缩试验进行了模拟,结果显示,该碎片替代模式能够很好地模拟颗粒破碎的过程.Lobo-Guerrero等[41]采用碎片替代法,通过用许多不同大小颗粒的组合替换一个在拉力中失效的颗粒,模拟颗粒材料在不同双轴应力组合作用下的破碎演化,并将颗粒材料在这些条件下的破碎过程可视化.

此外,Ben-Nun 等[42]采用两种新的替代模式对颗粒集合体在侧限压碎时的拓扑结构对自组织特性的影响进行了分析.杨贵等[43]基于Ben-Nun 提出的替代模式发展了两种新的碎片替代模式用于粗粒料破碎行为的研究.以上研究证明基于碎片替代法的颗粒破碎模拟可有效分析可破碎颗粒的力学响应.然而,以上均为二维替代模式,真实颗粒均处于三维环境中,且三维情况下颗粒材料的力学行为更加复杂,开展三维颗粒破碎研究具有更加重要意义.

近些年,三维情况下基于碎片替代法的颗粒破碎分析方法也得到了发展.例如,McDowell 等[44]通过在初始颗粒边界内放入有一定重叠量的轴对称的等粒径小颗粒的方式,构建3 种碎片替代模式,用于侧限压缩试验中颗粒细观力学特性分析,如图12 所示.Marketos 等[45]采用构建的替代模式研究了砂岩颗粒的开裂过程,重点关注了压缩带的萌生和扩展过程.Tavares 等[46-47]根据应力的施加位置人为设置母颗粒内部的颗粒碎片,其中较大的碎片在垂直于引起破碎的应力方向上,并且碎片颗粒在该方向上具有重叠量,较小的碎片排列在剩余的空间,通常较小的碎片与较大的碎片之间也有重叠量,如图13 所示.初始状态的重叠量是为了保证体积守恒,但这会引起母颗粒产生爆炸力,为了消除这个爆炸力,其中一个方法是使用松弛因子来限制由粒子重叠量计算得来的法向力.另一种简单的方法是人为消除碎片生成时产生的重叠量,只允许母颗粒之间真正重叠量的存在.

图12 碎片替代法研究压缩作用下颗粒破碎的演化过程[44]Fig.12 Fragment substitution method for studying the evolution process of particle fragmentation under compression[44]

图13 碎片的替代方案[47]Fig.13 Alternative scheme for fragments[47]

由于碎片替代法中替代颗粒的尺寸直接影响模拟效果,因此,众多的专家学者对颗粒破碎过程的尺寸效应进行了研究.Ciantia 等[48]采用碎片替代法进行了侧限压缩数值试验,研究了颗粒破碎的尺寸效应.Zhou 等[49]借鉴等粒径三碎片替代方法研究了主应力对可破碎颗粒材料宏细观力学响应的影响.为了使碎片粒径分布尽可能地符合实际情况,Zhou 等[50]又构建了新的三维碎片替代模式,对可破碎颗粒破碎过程的尺寸效应进行了研究.Cil 等[51]采用所提出的替代模式,进一步研究了颗粒破碎能量的尺寸效应.通过以上研究,对颗粒破碎过程中的尺寸效应有了更深入的理解和探索.Tavares 等[46]认为颗粒的破碎程度可以通过一个参数t10来表征,t10表示小于母颗粒尺寸 1/10 的碎片的比例,并通过以下表达式来计算颗粒的破碎程度,其中,A和b′为由实验数据拟合的模型参数,E50b为发生破碎的颗粒的断裂能的中位数,为破碎过程中颗粒可用的能量.t10值越高,子代颗粒尺寸分布越细.

EDEM 在考虑单颗粒破碎时是基于碎片替代法实现的.EDEM 在研究单颗粒破碎过程时可直接获取粒径分布、破碎率等数据,直观展示破碎效果,模拟速度非常快,但模拟结果目前无法体现破碎后的物料料型.

1.2.2 碎片替代法的破碎准则研究

目前,采用碎片替代法研究颗粒破碎过程时,破碎准则并没有统一的标准.现存在的破碎准则均基于一定的假设并做了相应的简化.Åström 等[37]认为基于碎片替代法的破碎准则主要可以归结为两类:一是基于应力的判定准则,即通过对颗粒的应力设定阈值,当应力大于所设定的阈值时颗粒发生破碎;二是基于力的判定准则,即通过对颗粒的接触力设定阈值,当最大接触力大于所设定的阈值时颗粒发生破碎.其他的破碎准则均是在这两种准则的框架下发展起来的.

假设最大接触力引起颗粒局部应力集中从而导致颗粒破碎[52-53],那么最大接触力可以作为破碎准则.Lobo-Guerrero 等[41,54-55]在二维模拟中对颗粒受力情况进行了简化,认为当接触数目不大于3 时,颗粒的受力情况可以类比于巴西圆盘试验.基于以上假定,采用颗粒的特征拉应力 σt=2P1/(πLD) 作为颗粒破碎的判定条件,即 σt＞σtmax,也 即P1＞σtmaxπLD/2,其中P1为作用力,L为单位长度,D为盘直径.

由于配位数对颗粒破碎有一定的影响[40,56],在力判定准则中引入配位数的影响,以上假定不再适用,即配位数大于3 时颗粒仍可能发生破碎[57].根据物理试验结果建立球形颗粒在径向加载时的最大流动剪切强度公式[58],当配位数在6～12 时该公式同样适用[59].通过以上公式Ciantia 等[48]求得颗粒最大接触力阈值,进而建立三维模拟中颗粒破碎的判定准则,即kmob≤k,其中kmob和k分别为颗粒的流动强度和固有强度,当流动强度大于固有强度时颗粒发生破碎.选取3 种不同破碎难易程度的材料,张科芬等[60]对其进行破碎模拟,证明以上破碎准则适用于离散元中球形颗粒的破碎模拟,并且能够反映颗粒的真实破碎过程.Cil 等[51]基于颗粒材料特征强度的3 种尺寸效应公式发展了3 种不同的力判定标准,并通过侧限固结试验研究了采用上述3 种力判定准则时颗粒宏细观力学行为的不同.

根据均匀化思想发展了应力判定准则假设外载荷作用下颗粒内部的应力均匀分布,是颗粒受力状态的整体考量.半径为R的颗粒受到任意力作用时,其平均应力张量为是颗粒的接触数目,分别为接触点(c)处的单位法向量沿i方向的分量和沿j方向的接触力,Vp为颗粒的体积.其中,以等效力张量为破碎准则的本质是通过颗粒的特征强度判定颗粒是否发生破碎[57].颗粒的受力情况可以简化为静水压力和偏应力,采用等效二维Drucker-Prager 准则判定颗粒的破碎情况[40].然而,由于颗粒间的接触十分复杂,当存在多个接触力时颗粒可能处于高静水压力、低偏应力的状态,在这种情况下颗粒的破碎现象不容易发生,此时平均应力判定准则不在适用.针对以上情况,八面体剪应力破碎准则被提出并通过单颗粒压碎试验获得该准则中特征拉应力的阈值[44,61];带截断的Mohr-Coulomb 准则也可有效解决平均应力判定准则不适用的情况[49-50,62];杨贵等[43]通过对颗粒的拉裂破坏和剪切破坏情况提出了不同的破碎准则.

以上两种判定准则在一定程度上可以反映颗粒的破碎机制,然而,由于其均在一定假设的前提下提出的,且进行了一定的简化,所以无法完全描述颗粒复杂的破碎行为.此外,对破碎准则的选取并没有特定的标准,选取不同的破碎准则获得的颗粒系统孔隙率和级配的演化过程的正确性也并未得到验证.此外,二维和三维问题存在较大差异.因此.针对三维问题提出一套全面的、准确的选取破碎准则的方法是十分必要的.

1.2.3 碎片替代法在工程中的应用

由于碎片替代法可以实现二次破碎,且对于不考虑体积守恒和对破碎颗粒形态要求不严格的工程问题,可以采用该方法实现颗粒的破碎研究,同时,碎片替代法在大规模的破碎过程问题研究中也有较大的优势,因此在以下工程问题中得到广泛应用.

碎片替代法在岩土工程中具有广泛的应用,如通过模拟岩体破碎的过程来预测岩体的破碎区域、岩块大小分布、岩体的稳定性等[63];通过模拟地下开挖过程中岩土体的破碎行为,预测岩土体的破碎区域、位移和变形情况,这对于开挖工程的设计和施工具有重要意义[64];该方法还应用于模拟岩土体在振动载荷下的破碎行为,如在振动台试验中模拟土体的破碎过程,研究振动对土体强度和稳定性的影响;在模拟地下爆破引起的岩土体破碎过程中可获得爆破对周围岩土体的影响,如爆破震动对周围地基的影响、爆破引起的岩层破碎区域等;此外,碎片替代法可以用于预测岩土材料的力学性质,如弹性模量、断裂强度等,这对于岩土工程设计和材料选取具有重要意义[65].

Lobo-Guerrero 等[55]发展的碎片替代模式包含3 种不同大小共8 个子颗粒,并通过二维离散元数值方法研究了由于颗粒材料破碎引起的铁路道砟级配的变化对铁沉降变形的影响,如图14 所示.

图14 碎片替代法在铁路道砟颗粒材料破碎中的应用[55]Fig.14 Application of fragment substitution method in railway ballast particle breakage[55]

张科芬等[60]基于局部应力集中的点载荷破碎准则,采用碎片替代法开展了3 种不同破碎难易程度材料的数值模拟.其中,为了保证破碎前、后颗粒之间的平衡,利用阿波罗填充和膨胀法,并引入尺寸因子来表征不同粒径的颗粒强度.

1.2.4 碎片替代法的优势及存在的主要问题

基于碎片替代法的颗粒破碎数值方法在一定程度上再现了颗粒的破碎过程,且该方法能够实现颗粒的二次破碎.碎片替代法为较大规模颗粒破碎过程分析和细观尺度的机理研究提供了一定的思路.然而,由于碎片替代模式的主观性导致模拟结果与真实的颗粒形态及破碎过程存在较大的差异;三维模拟中颗粒破碎的判定准则并不完善.对以上问题的有效解决将促进碎片替代法的发展和应用.

2 基于离散元−有限元耦合方法的颗粒破碎数值算法

当参考连续介质力学对材料破碎的研究方法并结合颗粒材料的离散特性,可采用离散元−有限元耦合的方法对颗粒破碎进行数值研究.以下将介绍比例边界有限元方法、组合有限元−离散元方法以及典型的内聚力模型.

2.1 比例边界有限元法

2.1.1 比例边界有限元法的提出与发展

比例边界有限元法是由Song 等[66]提出的一种新的半解析方法,该方法是基于有限单元和边界元方法提出的,同时拥有有限元法和边界元法的优点,已经在一些工程领域中的动力学和静力学问题研究中得到了应用.比例边界有限元法相对于有限元法和边界元法计算过程更简单有效,并可得到工程满意的结果.Ooi 等[67]发展了适用于任意边数的凸多边形比例边界有限元法.凸多边形比例边界有限元仅用一个多边形单元就可以模拟具有任意边数的凸多边形颗粒,并得到其变形和应力状态.

比例边界有限元法实现颗粒的破碎过程时需要对每个颗粒的应力−应变进行完整分析,然后根绝准则来判定颗粒内部的破坏点,当破坏点的数量达到一定比例后,认为该颗粒发生破碎.比例边界有限元法认为颗粒破碎的路径是直线,所以采用加权最小二乘法对破坏点进行拟合,将破碎点连接起来获得破碎路径,从而使颗粒破碎为两部分.该方法不需要对破坏路径的位置和方向进行预定义.破碎后产生的颗粒在下一时间步中参与计算.罗滔等[36]采用比例边界有限元方法与离散元方法结合,对堆石料的双轴压缩实验进行了模拟,结果如图15 所示.

图15 加载结束后不同围压下的颗粒破碎分布[36]Fig.15 Particle fragmentation distribution under different confining pressures after loading[36]

2.1.2 比例边界有限元法的破碎准则

为了研究颗粒的破碎行为,需要建立相应的破碎准则,由于比例边界有限元方法已经计算获得了每个颗粒内部的应力场,因此,采用基于应力的破坏准则判断颗粒是否破坏.判断破坏的流程为: 首先在颗粒内部选取具有代表性的样本点,如图16 所示,基于Hoek-Brown 破碎准则判断样本点是否破坏;当破坏样本点的数量达到一定比例时颗粒破坏;比例边界有限元方法假设颗粒的破坏路径为直线,并通过加权最小二乘法对破坏点进行拟合得到破坏路径;沿着破坏路径将颗粒分为两个颗粒[68].

图16 样本点的分布Fig.16 Distribution of sample points

为了判断样本点是否破坏,应用Hoek-Brown 破坏准则[68]

式中,σt为抗拉强度;σ1f是由式(4)得出的第1 主应力,如果安全系数Fs＜1 则表示此点破坏,由Hoek-Brown 准则得到的破坏点如图17 所示.

图17 由Hock-Brown 准则得到的破坏点Fig.17 Failure points obtained from the Hock-Brown criterion

2.1.3 比例边界有限元法的优势及存在的主要问题

在一定程度上,比例边界有限元方法从细观尺度揭示了颗粒破碎机制.比例边界有限元方法计算颗粒内部的应力分布,从而实现颗粒的破碎过程,颗粒破碎产生的新的颗粒将直接参与到离散元和比例边界有限元的计算中,比例边界有限元方法对颗粒破碎的研究无需预定义任何子颗粒或者网格重划分[36].然而,该方法需要人为设计样本点的分布,且其破碎路径均为直线,因此,比例边界有限元法的主观性很强,很难客观反映颗粒破碎真实过程,发展从细观尺度反映颗粒破碎真实过程的数值方法迫在眉睫.

2.2 组合有限元−离散元方法

2.2.1 组合有限元−离散元方法的提出与发展

组合有限元−离散元方法(FDEM)是Munjiza等[69]在 20 世纪90 年代为解决混凝土加载中裂纹的萌生和扩展问题发展的数值方法.组合有限元−离散元方法首先对颗粒进行单元划分,并通过接触面单元对其进行连接,离散元法计算获得颗粒之间的作用力作为每个颗粒应力计算的边界条件,从而实现单个颗粒的应力−应变分析,获得颗粒的应力状态和变形情况.组合有限元−离散元法中颗粒沿有限单元的边界破坏.

FDEM 采用FEM 理论作为颗粒的变形、压裂破坏标准,通过DEM 概念来检测新的接触点,解决离散单元的平移、旋转和相互作用问题[70].FDEM利用分布接触力方法和罚函数方法解决接触力问题,有效地处理了复杂形状的变形问题[71].FDEM 在非线性弹性断裂力学(NLEFM)框架下模拟裂缝的起裂和扩展过程,此外,在FDEM 中可以实现重网格技术,但在建模可破碎颗粒材料时,计算成本非常高.由于在模拟过程中没有执行重网格,网格拓扑也没有更新,因此在加载过程中粒子不能继续变小.

基于FDEM 发展的颗粒破碎仿真软件主要有中科院力学所李世海团队的GDEM,该软件采用有限元和块体离散元耦合的方法对岩土工程中的颗粒材料破碎问题实现了数值分析.基于多物理场断裂分析软件MultiFracS 的基本理论也是有限元−离散元耦合算法,该软件在岩土材料断裂及破碎数值建模方面有独特优势,在隧道、边坡、采矿、水利水电、爆破、非常规油气及地热开采等领域有重要应用价值.

2.2.2 组合有限元−离散元方法在工程中的应用

由于FDEM 方法在模拟颗粒破碎的过程时,对每个颗粒进行应力分析,因此破碎过程与实际更加相符.在对单个颗粒的破碎分析中可有效揭示颗粒的破碎机理.因此针对考虑破碎机理的工程问题,以及需要模拟的破碎过程更加接近真实情况时,可采用该方法实现.

组合有限元−离散元方法在岩土工程中得到广泛的应用,如严成增等[72]将其与数值图像技术结合,系统分析了岩体的破裂机制;Ma 等[73]采用组合有限元−离散元方法建立了岩石颗粒的数值模型,如图18 所示,并对岩石颗粒在轴压作用下的断裂过程进行了模拟,系统研究了破碎前颗粒的形貌对破碎后碎片的形态和尺寸分布的影响,如图19 所示.此外,在组合有限元−离散元方法中引入基于应力的破坏准则,并对堆石料的破碎过程进行模拟,模拟过程中无需预先定义或者假设颗粒的破坏路径[74-75].Mahabadi等[76]使用三维FDEM 方法研究Opalinus clay 黏土页岩地层中隧道周围的开挖损伤区的发展.通过巴西圆盘试验和单轴压缩试验,对三维FDEM 模型进行校准,结果如图20 所示.

图18 岩石颗粒的FDEM 数值模型[73]Fig.18 Numerical model of rock particles[73]

图19 FDEM 方法研究岩石颗粒在轴压作用下的断裂结果[73]Fig.19 FDEM method for studying the fracture process of rock particles under axial compression[73]

图20 采用FDEM 方法进行巴西圆盘和单轴压缩试验[76]Fig.20 Brazilian disc and uniaxial compression tests using FDEM method[76]

Rougier 等[77]将离散裂纹与塑性变形相结合,并采用FDEM 方法模拟二维泰勒杆的断裂过程,如图21 所示.其中变形采用FDEM 方法描述,塑性变形在材料嵌入坐标系中表述,乘法分解和塑性流动在拉伸空间中求解.该工作将FDEM 中的裂纹和破碎标准有效结合.

图21 FDEM 方法模拟二维泰勒杆的断裂过程[77]Fig.21 FDEM method simulation of the fracture process of twodimensional Taylor bars[77]

针对脆性材料的断裂和破碎响应,Chen 等[78]采用FDEM 方法对夹胶玻璃在硬体冲击下的断裂和破碎过程进行了模拟,如图22 所示,并讨论了玻璃、层间和层间界面的破坏模型,证明了该模型在模拟夹层玻璃破裂时是可靠的.Lilja 等[79]采用FDEM 方法建立海冰的三维数值模型,如图23 所示,该模型由同转动、黏性阻尼Timoshenko 梁有限元组成的平面梁晶格与形成实际海冰的刚性离散单元的质量质心相连接.通过质心Voronoi 镶嵌程序生成网格,由于内部存在阻尼,基于晶格的结构的机械响应依赖于应变速率和尺寸.

图22 夹胶玻璃在硬体冲击下的断裂和破碎响应[78]Fig.22 Fracture and fragmentation response of laminated glass under hard impact[78]

图23 采用FDEM 方法建立海冰的三维数值模型[79]Fig.23 Establishing a three-dimensional numerical model of sea ice using FDEM method[79]

2.2.3 组合有限元−离散元方法的优势及存在的主要问题

FDEM 方法可有效模拟颗粒的破碎过程,由于该方法对每个颗粒进行应力分析,因此破碎过程与实际更加相符.在对单个颗粒的破碎分析中可有效揭示颗粒的破碎机理.然而,由于模拟过程中需要对每个颗粒进行有限单元划分,且颗粒的破碎只能沿着网格进行,为了获得更为真实的破碎过程,需要网格精细化,这均导致计算成本非常大[36].这将限制了其在工程中的应用,对实际工程问题的模拟无法起到推动作用.

2.3 内聚力模型

2.3.1 内聚力模型的提出与发展

FDEM 方法中,有限元离散化以后,将黏性界面单元插入到与每个粒子相关的有离散单元中,以考虑粒子的破碎过程.本节将重点介绍黏性界面的内聚力模型(cohesive zone model,CZM)[9,80-82].内聚力模型是1960 年Dugdale 等提出的[83-84],假设裂纹的尖端是由两个裂纹界面组成的很小的内聚区,该内聚区的本构关系就是界面上作用牵引力T与两裂纹面间相对位移U之间的关系.内聚力模型和裂纹尖端内聚区的分布[85]如图24 所示.

图24 内聚力模型和裂纹尖端内聚区[85]Fig.24 Cohesion model and crack tip cohesive zone[85]

内聚力模型从未完全承载的点A开始,T随着U的增加而增加,随之达到应力最大值Tmax的点C,此时材料的应力承载达到最大值,材料在该点处开始出现初始损伤.随着界面位移的继续增大,应力开始下降,该阶段为材料的损伤扩展阶段,点E为材料裂纹界面完全分离点,此时材料的承载力降为0[85].内聚区应力的变化受内聚力模型和裂纹界面位移的影响,因此,针对不同的材料,应该选择不同的内聚力模型,同时通过选取适当的参数,从而反映界面的强度、韧度等力学性能.

内聚力模型可分为基于位移的内聚力模型和基于势能的内聚力模型.在基于位移的内聚力模型中,裂纹上下表面之间的有效牵引力可由有效分离位移的函数表示,通常称之为牵引力分离法则.常用到的牵引力分离法则包括线性软化、双线性软化、指数和梯形法则等.不同法则的区别在于与之间函数关系的不同.将有效牵引力与内聚强度σmax归一化处理的结果如图24(b)所示.基于位移的内聚力模型出现的裂纹扩展问题,可以在基于势能的通用内聚力模型中得到解决.基于势能的内聚力模型采用三次多项式表示法向牵引力,采用线性关系式表示切向牵引力[85].

早期非线性断裂研究一般认为,当内聚区尺寸小于裂纹和试样尺寸时,内聚力模型理论与Griffith的能量平衡理论等效.Dugdale[83]认为内聚应力的分布在数值上与材料的屈服强度相等,但这有悖于实际物理事实,Barenblatt[84]考虑内聚力分子尺度的特征,把内聚应力看作内聚区裂纹面各点处裂纹张开位移的函数,但具体的解析式无法准确给出,在大多数情况下内聚力仍被假设为常数.基于以上工作,Hillerborg 等[86]考虑颗粒的拉伸强度,采用有限元方法建立脆性材料的内聚力模型并对其断裂过程进行模拟,实现了原有裂纹的增长、新裂纹的萌生和扩展研究,对断裂细节进行了完整的描述.Needleman[87]通过高次多项式函数对延性材料的断裂情况进行了模拟.Kolher 等[88]采用了分段函数描述内聚力模型的方法对镍铝合金的剪切断裂性能进行了模拟.

为模拟脆性岩体的断裂行为,清华大学徐文杰团队[89]提出了一种针对多面体离散元颗粒设计的内聚断裂模型(CFM),如图25 所示.岩体被离散成一系列刚性多面体块体,这些块体沿边界面的法向和剪切方向黏结在一起.黏性准则规定了CFM 的法向和剪切断裂强度.通过在Blaze-DEM 框架中对多个图形处理单元的并行处理,提高了CFM 在多面体模拟中的计算效率.

图25 断裂前、后的接触模型[89]Fig.25 Model after fracture[89]

2.3.2 内聚力模型在工程中的应用

内聚力模型在处理非线性和大变形问题时具有非常强的适应性.随着内聚力模型的改进和发展,在多种材料的断裂问题上具有独特的优势,因此,在处理此类工程问题时可采用该方法.

Gui 等[90]在通用离散元编码(UDEC)中,引入了结合拉伸、压缩和剪切材料特性的内聚断裂模型,用于模拟岩石动力试验中的压裂过程,通过缺口半圆弯曲试验和巴西盘试验两个数值算例,验证了内聚力模型对岩石动态破坏模拟的有效性.

为了模拟岩石破裂过程,广东工业大学高伟团队提出一种新的内聚模型本构.在该本构中,采用Mohr-Coulomb 准则作为裂纹萌生准则,建立压缩和拉伸状态下的牵引−分离模型,描述了岩石强压力依赖性的力学特性.在粘结区完全破坏后,采用基于库仑公式的摩擦模型来确定裂纹自由面之间的接触相互作用.采用所提出的零厚度黏合单元本构,模拟了巴西盘、单轴压缩两个岩石破裂实例.比较仿真结果与实验结果发现以上模型是可靠的[9],如图26～图27 所示.

图26 0 厚度黏合单元本构[9]Fig.26 Constitutive law of zero thickness adhesive element[9]

图27 采用内聚力单元模拟巴西盘破裂和单轴压缩岩石破裂实例[9]Fig.27 Simulation of brazilian plate rupture and example of rock fracture using cohesive elements[9]

2.3.3 内聚力模型的优势及存在的主要问题

内聚力模型的理论基础是弹塑性断裂力学理论,该模型认为断裂问题是非线性边值问题,不必考虑断裂的起裂和扩展准则,因此,适应性强,可以解决非线性和大变形问题.随着计算机技术和有限元方法的发展,内聚力模型得到了改进和发展,已经可以解决多种材料的断裂问题.然而,针对颗粒型材料,仍然存在计算效率的问题,因此,该方法在研究颗粒型材料的破碎问题中仍没有得到推广.

3 近场动力学方法及其他方法

3.1 近场动力学方法

3.1.1 近场动力学方法的提出与发展

近场动力学(peridynamics,PD)是一种新兴的非局部性理论,针对模型中的非连续问题,该方法不需要任何处理,对模型中从连续到非连续的过程可以很好地处理[91-92],因此,该方法在研究离散介质的破碎问题时有一定的优势.近场动力学方法是Silling等[93]于2004 年提出的,该方法基于长程非局部作用思想建立力学模型,通过求解空间积分型运动方程来描述物质力学行为,实现了对连续介质和非连续介质力学行为的统一描述.此外,由于损伤概念的引入,近场动力学理论对材料的损伤和裂纹的形成过程可以有效描述.因此,在处理裂纹自发萌生和扩展等非连续问题时,近场动力学理论具有一定的优势[94].

近场动力学方法将固体材料离散为一系列具有体积和质量的物质点,每个物质点与其周围一定区域内所有其他物质点之间存在相互作用力.如图28所示.

图28 近场动力学理论的局部接触式作用和长程非局部作用[91]Fig.28 Local contact type action and long-range non local action in peridynamics theory[91]

近场动力学理论是一种非局部连续介质力学理论,采用积分的形式表示邻域半径 δ 内物质点x和x′之间的相互作用力,即键力.物质点x的动力学方程为

近场动力学理论通常被分为键型(bond-based)和态型(state-based)两种.其中键型近场动力学理论认为两物质点之间的相互作用仅包含两物质点间的键力.而态型近场动力学理论在分析两物质点间的相互作用时考虑了邻域内所有物质点对键力的影响.显然,键型近场动力学理论是态型近场动力学理论的一种特例[91].

计算程序编写促进了近场动力学理论的发展和应用.Silling 团队在1998 年采用Fortran 90 语言开发了近场动力学方法的代码EMU,并对其进行了测试,这是第一次将近场动力学理论程序化;Parks 等基于分子动力学软件LAMMPS 开发了软件PDLAMMPS,扩展了近场动力学理论的应用;Macek 等在经典动力学分析软件ABAQUS 中引入近场动力学算法,促进了近场动力学理论的发展;此外,LS-DYNA 中也引入了近场动力学方法扩展程序包,并命名为DYNAPeridynamics,有效地结合了现有商业软件和近场动力学方法在处理动力学问题的优势;近期,桑迪亚国家实验室基于离散单元法开发了开源软件Peridigm,该软件兼容了Cubit 网格生成工具和Paraview 可视化程序,将近场动力学方法平台化,具有人机交互性好、计算稳定、后处理方便、实用性强等优点[91];赵吉东团队[95]开发了PD-NSCD 混合法的代码,并应用于颗粒材料连续破碎过程的研究,自主开发代码实现近场动力学理论部分,使用开源代码Project Chrono 进行基于NSCD 的离散建模.同时,通过MATLAB,Python,C 语言,Fortran 等语言进行近场动力学数值计算与分析也是部分学者在处理实际问题时采用的手段[96-97].

赵吉东团队[95,98-99]开发的PD-NSCD 混合法,采用非光滑接触动力学模拟颗粒系统的刚体运动以及颗粒间相互作用,利用近场动力学理论分析单个颗粒的破碎.基于近动力学和物理引擎提出一种新的混合计算框架来模拟机械载荷下的可破碎颗粒材料.在此框架下,利用近场动力学方法模拟单个粒子的破碎,利用基于非光滑接触动力学方法的物理引擎模拟粒子的刚体运动和粒子间的相互作用.该混合框架可以对颗粒破碎过程进行模拟,并实现连续破碎过程中不规则颗粒形状的表征和模拟,克服了许多现有方法面临的缺陷和挑战.其工作为研究颗粒破碎的复杂过程提供了一种有效的方法,为今后研究可破碎颗粒材料的微观结构行为提供了参考.图29 显示了近场动力学分析初始化和构造粒子的完整过程[98],其中,图29(a)为受接触力作用的颗粒,图29(b)为对离散粒子进行的近动力学分析,图29(c)未颗粒破碎的近动力学分析结果,图29(d)为物理引擎中的模拟粒子,图29(e)为破碎粒子的拆分视图.图30～图31 为砂石颗粒连续破碎过程的数值模拟结果[95,98].

图29 颗粒破碎建模过程示意图[98]Fig.29 Schematic diagram of particle breakage modeling process[98]

图30 球形颗粒模拟砂石颗粒连续破碎过程[98]Fig.30 Simulation of continuous breakage process of sand and stone particles with spherical particles[98]

图31 多面体颗粒模拟砂石颗粒连续破碎过程[95]Fig.31 Simulation of continuous breakage process of sand and gravel particles using polyhedral particles[95]

3.1.2 近场动力学方法在工程中的应用

目前,近场动力学理论已经在多个领域得到了应用[100].如混凝土压碎[93,101]、金属冲击[102]、子弹穿透[103-104]、土体[105-107]和陶瓷破碎[108]、薄膜撕裂[109]以及玻璃冲击[110]等.目前,特别是在岩石力学研究中,近场动力学方法已经成功模拟了水压致裂破坏过程[111],对于预制裂纹的萌生、扩展和破坏过程模拟,近场动力学方法也是有效的.通过近场动力学理论,Ha 等[112]、朱其志等[113]和王振宇等[114]研究了单条预制裂纹在单轴压缩条件下裂纹的演化过程;Lee等[115-116]研究了含有预制裂纹岩石的翼型裂纹和次生裂纹的形成过程;Rabczuk 等[117]采用双邻域近场动力学模型对含多条随机裂纹的岩石的破坏过程进行了研究.此外,Zhou 等[118-121]引入摩尔库伦准则和极限张拉破坏准则,对预制裂纹的巴西圆盘[120-121]、三点弯曲[119,121]、单轴压缩[119]以及双轴张拉[118]进行了数值模拟.

Hu 等[110]采用近场动力学方法分析了玻璃板的损伤模式,其中冲击速度从61～200 m/s 不断增加.将实验结果与简化模型的近场动力学模拟结果进行了比较,如图32 所示.实验观察到的主要断裂模式是由3 种冲击速度的每一种测试的近场动力学模型捕获的.该方法通过进一步改进边界条件以实现对多层体系的冲击造成的脆性损伤更准确的模拟.近场动力学计算模型揭示了玻璃层复杂脆性损伤演化的早期阶段,以及边界条件对动态破裂过程的影响.

图32 撞击速度为61 m/s 的破碎图,以上分别是玻璃正面、底面和实验结果[110]Fig.32 Fragmentation diagram with an impact velocity of 61 m/s,showing the front,bottom,and experimental results of the glass[110]

Rabczuk 等[117]提出了一种用于颗粒状和岩石状材料断裂的双层近场动力学(DH-PD)方法.与离散裂纹方法相比,DH-PD 方法不需要任何裂纹表面的表示形式,也不需要任何标准来处理复杂的裂纹形态,裂纹路径是模拟的自然结果.在该方法中作者提出了一种新的模拟压缩裂缝接触和约束侵彻条件的惩罚方法.该方法被应用于地质力学的一些基准问题,图33 为采用双层近场动力学方法研究方形板裂纹扩展的数值模拟结果.

图33 双层近场动力学研究方形板中裂纹的扩展[111]Fig.33 Double layer Peridynamics study on crack propagation process in square plates[111]

Shen 等[101]采用近场动力学理论分析混凝土结构的损伤和渐进破坏.建立了混凝土柱单轴受压非局部近动力模型,如图34 所示.模拟结果表明,混凝土近动力体的裂缝是自发形成的.此外,模拟结果也展现了基于近场动力学理论建模混凝土损伤累积和渐进破坏的优势.Shen 等的研究为分析复杂的不连续问题提供了一种新的思路.

图34 单轴压缩混凝土柱的近场动力学模拟[101]Fig.34 Peridynamics simulation of uniaxial compression concrete columns[101]

3.1.3 近场动力学方法的优势及存在的主要问题

近场动力学方法克服了基于局部作用思想建模和求解空间微分方程在处理非连续问题时遇到的奇异性困难,突破了经典分子动力学方法在计算尺度上的局限,在宏观和微观非连续力学问题分析中具有明显的优势[91].仅经过十几年的发展,近场动力学在模拟材料非连续破坏方面已展现出明显的优越性.因此,随着近场动力学理论与应用研究的不断深入将促进岩石力学的进一步发展.

目前,近场动力学模型中的“边界效应”问题是需要解决的一个重要问题[122].对近场动力学模型中自由边界附近的物质点进行积分时,会出现邻域不完整和受力不平衡的现象,这导致了“边界效应”问题的产生.无外力作用下发生自由膨胀时,由于积分域不完整导致模型自由边界附近物质点的键力分布不对称,导致物质点的受到非平衡力,降低了计算精度.键型近场动力学模型中的微模量c被视为常数也会引起“边界效应”问题.目前,主要采用以下办法处理近场动力学模型中“边界效应”问题: (1)对自由边界附近物质点进行体积修正或虚构边界层;(2)进行力密度修正;(3)能量密度修正[91].然而,上述方法仅仅削弱了“边界效应”的影响,而没有解决边界问题的根本.

近场动力学方法的计算效率相对传统数值计算方法较低,主要原因在于其在邻域半径范围内进行积分,且求解过程以显式计算为主.因此,计算时间步需要足够小才能获得稳定解,这将导致计算量增大.虽然近场动力学方法在解决非连续问题上具有优势,但其不完善性也限制了其在大规模问题上的应用.为了提高计算效率,目前研究者们将近场动力学方法与传统数值方法相结合.同时,随着计算机技术的快速发展和近场动力学理论的不断完善,计算效率问题正在逐步得到解决.

3.2 其他方法

除了以上论述的数值方法外,还有一些方法本文未详细介绍,如Harmon 等[123]提出的水平集切割法,武汉大学王桥团队[124-125]和周伟团队[126]发展的相场理论方法,刘彪等[68]采用的边界元方法,以及物质点法等.这些方法都在颗粒破碎的研究过程发挥了重要的推动作用,为今后发展更为可靠、高效的数值方法奠定了基础.

4 总结与展望

4.1 总结

颗粒破碎是颗粒材料研究中的关键内容.目前,颗粒破碎的研究已经取得了一定的成果,并在工程领域得到广泛的应用,解决了一定的实际问题.然而,现有的方法大多是从现象出发,缺乏对颗粒破碎机理的明确研究,因此,截止目前,对颗粒材料的破碎机理认识仍不清晰.此外,现有的方法一般预设了破坏路径,这与真实的破碎过程存在一定的差异;有的方法仅可以破碎一次,破碎后的颗粒无法实现二次破碎,这与实际也是不符的;若不提前预设破坏路径,且颗粒可以发生多次破碎,那么计算效率将会非常低.面对以上问题,亟待提出一种从机理出发不需要预设破坏路径且可以发生多次破碎的研究颗粒破碎的高效数值方法.

4.2 展望

针对颗粒破碎的数值研究,目前仍有一些问题存在,这将是今后研究的主要内容和方向.

首先,为了探究颗粒破碎的内在机理,则需要从破碎的本质出发,考虑颗粒真实破碎的内在原因.从以上内容不难看出,目前有关颗粒破碎的机理研究内容较少,为了揭示颗粒破碎的内在机理,必须将机理研究作为首要研究内容.

其次,当从颗粒破碎的现象描述出发,则需要考虑颗粒破碎程度的度量指标问题.对颗粒的破碎程度以及破碎后颗粒进行描述,这就需要有合理的描述指标来表示颗粒的破碎程度以及破碎后颗粒的相关性质.

然后,工程中颗粒材料破碎的计算效率问题: 从以上内容不难看出,目前的研究方法均存在计算效率问题,尤其是面向工程实际时,需要考虑的颗粒数量往往会达到百万甚至千万级,面对如此巨大数量的颗粒,如何保证计算效率是今后工作的重点,也是将颗粒材料破碎应用于工程的必要条件.

最后,颗粒破碎有效模拟的计算软件问题: 面对机理研究以及工程应用,如何更有效地研究颗粒破碎过程,研发有效的计算软件将是解决这一问题最为直接和有效的方法.因此,研发具有自主知识产权的模拟颗粒破碎的计算软件,也将是今后工作的重点,当然,这也是难点.