开级配大粒径沥青碎石OLSM-25细观结构的离散元分析初探

郭红兵，赵亚兰，陈拴发

（1.陕西交通职业技术学院公路工程系，陕西西安 710018；2.长安大学公路学院，陕西西安 710064；3.长安大学材料学院，陕西西安 710064）

开级配大粒径沥青碎石OLSM-25细观结构的离散元分析初探

郭红兵1,2，赵亚兰1,2，陈拴发3

（1.陕西交通职业技术学院公路工程系，陕西西安 710018；2.长安大学公路学院，陕西西安 710064；3.长安大学材料学院，陕西西安 710064）

沥青混合料是由沥青和粗细集料等粒状单元共同组成的多相复合材料，目前现有的沥青路面结构设计理论均将路面结构当作均匀连续、各向同性的弹性体进行计算分析，这与沥青混合料自身属于颗粒体非连续微细观结构的实际情况不相符．本文采用离散元方法，运用二维颗粒流程序（PFC2D），对比分析了两种颗粒生成方法的特点及其适用条件，以开级配大粒径沥青碎石混合料（Open-graded Large Stone Asphalt Mixes，OLSM）OLSM-25推荐级配的上限为1#级配，建立了1#级配OLSM-25混合料的二维离散元模型，验证了二维离散元方法用于沥青混合料颗粒体材料数值模拟的可行性，弥补了有限元方法的不足之处，为OLSM-25的细观结构分析开辟了新的解决思路和途径．

道路工程；开级配大粒径沥青碎石；细观结构；离散元法；数值分析

0 引言

沥青混合料是由沥青和粗细集料等粒状单元共同组成的多相复合材料，其工程性质相当复杂．在荷载作用下，组成沥青混合料的各颗粒体材料呈现出相对离散的运动状态．目前现有的沥青路面结构设计理论均将路面结构当作均匀连续、各向同性的弹性体进行计算分析，这与沥青混合料自身属于颗粒体非连续微细观结构的实际情况不相符[1]，离散元方法则为此类问题提供了解决思路与途径．本文采用离散元方法，运用二维颗粒流程序（PFC2D），对比分析了两种颗粒生成方法的特点及其适用条件，验证了二维离散元方法用于沥青混合料颗粒体材料数值模拟的可行性；根据OLSM-25混合料的大粒径多空隙结构特点，在OLSM-25推荐级配范围的基础上，建立了OLSM-25混合料的二维离散元模型．

图1 二维颗粒流程序(PFC2D)的计算过程Fig.1Calculation process of particle flow code in 2 dimensions(PFC2D)

图2 二维颗粒流程序（PFC2D）的两种接触形式Fig.2Two contact form of particle flow code in 2 dimensions(PFC2D)

1 离散单元法的基本理论

1.1 基本方程

离散元法是由CundalPA[2-4]（1971年）首先提出并应用于岩土体稳定性分析的一种数值分析方法．它是一种动态的数值分析方法，可以用来模拟边坡岩体的非均质、不连续和大变形等特点，因而，也就成为目前较为流行的一种岩土体稳定性分析数值方法．PFC2D（Particle Flow Code in 2 Dimensions）是基于二维圆盘单元的离散元颗粒流程序，它主要用于模拟大量颗粒元的非线性相互作用下的总体流动和材料的混合，含破损累计导致的破裂、动态破坏和地震响应等问题[5-7]．利用二维颗粒流程序（PFC2D）解决沥青混合料颗粒体材料的离散介质问题时，颗粒体材料之间的相互作用应满足平衡及物理方程[8-10]．二维颗粒流程序（PFC2D）计算过程见图1．

1.1.1 力-位移定律

力-位移定律通过“颗粒-颗粒”与“颗粒-墙体”2种接触形式，把相互接触的力与位移联系起来，2种接触形式[8,11-12]见图2．

对于颗粒-颗粒接触，其单位法向ni为

式中：xAi、xBi分别为颗粒A、颗粒B中心坐标矢量；d为两颗粒心之间的距离．

对于颗粒-颗粒接触，两颗粒心之间的距离d为

对于颗粒-墙体接触，其单位法向ni的确定方法见图3．

法向的相对接触位移定义为重叠量Un

式中：RA、RB分别为颗粒A、颗粒B的半径．

颗粒-颗粒、颗粒-墙体接触点的坐标可表示为

接触力Fi可分解为法向分量和切向分量

式中：Fni、Fsi分别为接触力的法向分量和切向分量．

1.1.2 运动定律

在二维颗粒流程序（PFC2D）中，颗粒体的运动方程[8,13-14]可表示为合力（线性运动）与合力矩（旋转运动）的组合形式

式中：Fi为作用于颗粒体上的外部合力；m为颗粒体质量；xi为颗

粒体加速度；gi为重力加速度；Mi为合力矩；Hi为转动角动量．

取颗粒体的主惯性轴作为局部坐标系，则该局部坐标系下的欧拉运动方程为式中：xiP为墙体上P点的坐标位置，xiP为墙体上点P速度，其他符号含义同上．

图3 颗粒-墙体接触中单位法向ni的确定Fig.3Determination of the unit normal vector niin ball-wall contact

在二维颗粒流程序（PFC2D）中，通过初始化命令设定颗粒体的初始速度，常量外力（矩）施加于颗粒体重心位置，并随着程序运行不断累加至当前值之上．

1.1.4 时步的确定

在二维颗粒流程序（PFC2D）中，通过有限差分法进行积分，且程序运行时步不能超过关键时步，程序在每一次计算循环开始之前自行确定关键时步数．

1.2 接触模型

二维颗粒流程序（PFC2D）提供了3种接触模型[1,8]，即接触刚度模型：接触力与相对位移呈弹性关系；滑动模型：允许颗粒体在限定范围内发生相对移动；粘接模型：允许颗粒体在粘接强度范围内相互接触，包括接触粘接、平行粘接两种．对于沥青混合料的离散元数值模拟问题，文献[1]表明平行粘接模型适合于模拟沥青胶浆裹覆在粗集料之间的粘弹性本构行为（见图4）．在沥青混合料材料组成中，沥青与矿粉形成沥青胶浆，填充于粗、细集料之间的空隙结构中，浸润在粗集料表面上形成稳定结构．采用平行粘接模型可以实现对包括粗集料、细集料、矿粉、沥青在内的沥青混合料的离散元数值模拟，弥补了有限元方法的不足之处．

1.3颗粒流程序的数值分析步骤

运用二维颗粒流程序（PFC2D）进行数值模拟分析的步骤[8,18-20]如下：1）依据颗粒流的运动定律和力-位移定律，定义模拟对象；2）建立力学模型的基本概念，根据颗粒间的接触本构模型，选取模拟的力学接触模型；3）构造并运行简化的力学模型；4）补充数值模型分析所需的数据资料；5）数值模拟分析之前的准备，选取及确定模拟参数，包括合理确定分析步、计算时长、循环次数等；6）运行计算模型；7）进行后处理，分析解释计算结果．

图4 二维颗粒流程序（PFC2D）的平行粘结模型Fig.4Parallel-bond model of particle flow code in 2 dimensions(PFC2D)

2 OLSM-25离散元试件的颗粒生成

在二维颗粒流程序（PFC2D）中，为了生成不规则排列的颗粒体试件，通常采用的方法是先限定颗粒生成区域，然后在限定区域内填充一定数量的颗粒，使之符合预期的密实度和目标空隙率要求，并且在限定区域内生成的颗粒体试件必须达到一个初始应力平衡状态．在以下几种颗粒生成方法中，颗粒生成区域的墙体边界始终固定不变，通过程序不断地“微调”，最终生成符合预期密实度和目标空隙率要求，并达到初始应力平衡状态的颗粒体试件．

2.1 扩大颗粒半径法

在二维颗粒流程序（PFC2D）中，GENERATE命令用于在限定区域内随机生成一定数量的颗粒，但是，在程序随机生成颗粒的过程中，当新生颗粒与已有颗粒或墙体发生位置重叠时，GENERATE命令将终止，不再生成新的颗粒．扩大颗粒半径法（GenerationbyRadiusExpansion）则克服了这一难题，该法先在限定区域内生成大量“小”半径颗粒，此时颗粒体试件尚未达到目标空隙率要求；然后通过扩大颗粒半径的方法，直到颗粒体试件达到目标空隙率要求时程序终止．这一过程通过PFC2D程序内置的FISH子程序实现，扩大颗粒半径法涉及以下程序参数．

式中：R0、R分别为颗粒的初始半径和扩大后半径，R=mR0，m为颗粒半径扩大系数；n0、n分别为颗粒体的初始空隙率和目标空隙率．

式中：R为颗粒平均半径；N为颗粒总数量．

由此可得GENERATE命令中颗粒半径的上限RHI、下限RLO及最大与最小颗粒的半径之比r为

在程序中设置以上参数，通过GENERATE命令生成颗粒体试件，并调用PFC2D程序内置的FISH子程序get_poros测定颗粒体试件的空隙率，然后循环运行程序直至颗粒体试件达到初始应力平衡状态之后，即可运行计算模型进行数值模拟分析，并分析解释计算结果．

2.2 激增颗粒数量排斥法

在二维颗粒流程序（PFC2D）中，BALL命令用于在限定区域内指定圆心位置生成指定半径的颗粒，激增颗粒数量排斥法（Generation byExplosiveRepulsion）正是应用BALL命令在限定区域内生成足够数量的、规定半径的颗粒，直至颗粒体试件达到目标空隙率为止．但是，在应用BALL命令生成颗粒的过程中，颗粒在限定区域内是随机分布的，这有可能导致部分颗粒相互重叠，并且重叠颗粒之间存在较大的相互作用力，这些颗粒间的作用力导致颗粒产生较大的初速度，颗粒的初速度大到足以使得部分颗粒从限定的墙体边界“溢出”．为了防止出现颗粒溢出现象，在程序运行的最初几个循环中反复多次将颗粒的动能减小为零，这样以来，程序才能逐步运行直至稳定平衡状态．当采用激增颗粒数量排斥法生成符合目标空隙率要求的颗粒体试件时，应用BALL命令生成颗粒，只需指定颗粒半径的上下限值，而不指定所需生成的颗粒总数．这样以来，应用BALL命令每次只生成一个颗粒，然后测定已生成颗粒体试件的空隙率，当实测空隙率远远大于目标空隙率时，程序继续循环生成颗粒，直至实测空隙率接近目标空隙率时循环终止．这一过程需要多次程序循环和能量消减才能达到最终的颗粒体平衡状态．

对比分析以上2种方法，采用“扩大颗粒半径法”虽然可以得到预期的空隙率，但是由于颗粒半径扩大，不能严格控制配合比（即各级筛孔颗粒的数量）；而采用“激增颗粒数量排斥法”既可得到预期的空隙率，也能严格控制配合比．故本文采用激增颗粒数量排斥法生成1#级配OLSM-25混合料的二维离散元模型．2.3OLSM-25混合料离散元模型

OLSM-25混合料由沥青、粗集料、细集料及矿粉组成，本文以国内外OLSM-25参考级配为基础，依据贝雷法采用0.22NMPS（取4.75 mm）作为粗、细集料的分界点，应用变I法原理设计OLSM-25级配，结合空隙率（15%～20%）要求，提出OLSM-25推荐级配范围，以OLSM-25推荐级配的上限、中值、下限为1#～3#级配，见表1．

表11 #～3#OLSM-25级配组成Tab.1Composition of 1#～3#OLSM-25 gradation

在二维颗粒流程序（PFC2D）中，采用激增颗粒数量排斥法生成1#级配OLSM-25混合料二维离散元模型．即按照1#级配OLSM-25混合料通过率，在限定矩形区域（50 cm×25 cm）内应用BALL命令随机分布生成颗粒，只需指定各级筛孔颗粒半径的上下限值，而不指定所需生成的颗粒总数；应用BALL命令每次只生成一个颗粒，然后测定已生成颗粒体试件的空隙率，当实测空隙率远远大于目标空隙率（15%）时，程序继续循环生成颗粒，直至实测空隙率接近目标空隙率时循环终止．

表2 采用平行粘结模型的OLSM-25计算参数Tab.2Calculation parameters of OLSM-25 with the parallel-bond model

在二维颗粒流程序（PFC2D）中，采用平行粘结模型模拟OLSM-25混合料中沥青裹覆并浸润集料的本构行为，采用平行粘结模型的OLSM-25计算参数见表2．

在加载之前先对试件施加一个足够小的初始锁定应力，在程序运行的最初几个循环中反复多次将颗粒的动能减小为零，有效防止了颗粒溢出现象．经多次程序循环和能量消减最终达到平衡状态，1#级配OLSM-25混合料离散元模型平均不平衡力及平均接触力见图5．由图5可知，在加载之前对试件所施加的足够小的初始锁定应力及足够多的程序循环次数，使得1#级配OLSM-25混合料离散元模型达到最终的颗粒平衡状态，确保了离散元数值模拟的收敛性及稳定性．

1#级配OLSM-25混合料离散元模型及其骨架接触力分布见图6．由图6可知，1#级配OLSM-25混合料属于骨架空隙结构，空隙率15%，由于粗集料含量较多，细集料含量较少，表现为骨架接触力大小及其分布不均匀（骨架接触力分布图中线宽表示接触力大小，位置表示接触力分布），在大粒径粗集料颗粒接触点处接触力较大，而细集料接触点处接触力较小，从而由大粒径粗集料之间良好的石-石接触形成骨架结构，能够承受重载并具有良好的抵抗车辙性能．

图51 #级配OLSM-25离散元模型平均不平衡力及平均接触力Fig.5The average imbalance force and the average contact force of discrete element model of No.1 grade OLSM-25

图61 #级配OLSM-25离散元模型及其骨架接触力分布Fig.6Discrete element model and skeleton contact force distribution of No.1 grade OLSM-25

3 结束语

1）根据OLSM-25混合料的大粒径多空隙结构特点，以国内外OLSM-25参考级配为基础，依据贝雷法采用0.22 NMPS（取4.75 mm）作为粗、细集料的分界点，应用变I法原理设计OLSM-25级配，结合空隙率（15%～20%）要求，提出OLSM-25推荐级配范围．

2）采用离散元方法，运用二维颗粒流程序，对比分析了两种颗粒生成方法的特点及其适用条件，以OLSM-25推荐级配的上限为1#级配，建立了1#级配OLSM-25混合料的二维离散元模型，验证了二维离散元方法用于沥青混合料颗粒体材料数值模拟的可行性，弥补了有限元方法的不足之处，为开级配大粒径沥青碎石OLSM-25细观结构分析开辟了新的解决思路和途径．

[1]王端宜，赵熙．沥青混合料单轴压缩试验的离散元仿真[J]．华南理工大学学报（自然科学版），2009，37（7）：37-41．

[2]Cundal P A，Strack O D L．A distcrete numerical model for granular assemblies[J]．Geotechnique，1979，29（1）：47-65．

[3]Cundall P，BoardM，TinucciJ．Numericalmodelingrelatedtorockburstresearch[C]//AppliedRockburstResearch(Proceedings oftheInternational Workshop May 1994)．Cereceda J C，Van Sint Jan M．Santiago：Editec Ltda，1995：43-59．

[4]Cundall P A．A computer model for simulating progressive large scale movements in blocky rock systems[C]//Proceedings of the Symposium of the International Society for Rock Mechanics(Nancy France)．Vol 1．1971：II-8．

[5]Cundall P A．A discontinuous future for numerical modelling in geomechanics[J]．Geotech Eng，2001，149（1）：41-47．

[6]Abbas A R．Simulation of the micromechanical behavior of asphalt mixtures using the discrete element method[D]．Pullman W A：Department of Civil and Environmental Engineering，Washington State University，2004．

[7]Cundall P A，Strack O D L．Modeling of microscopic mechanisms in granular material[C]//Mechanics of Granular Materials:New Models and Constitutive Relations．Jenkins J T，Satake M．Amsterdam：Elsevier Scientific Publications B V，1983：137-149．

[8]Cundall P A．PFC2D user's manual:Version 3．1[Z]．Minnesota:Itasca Consulting Group Inc，1999．

[9]Cundall P A．Computer simulationsof densesphereassemblies[C]//Micromechanics of GranularMaterials．SatakeM，JenkinsJ T．Amsterdam：Elsevier Science Publishers B V，1988：113-123．

[10]Cundall P A．Distinct element models of rock and soil structure[C]//Analytical and Computational Methods in Engineering Rock Mechanics，Chapter 4．Brown E T．London：George Allen and Unwin，1987：129-163．

[11]Cundall P A，Drescher A，StrackO D L．Numerical experiments ongranular assemblies，measurements and observations[C]//Deformation and Failure of Granular Materials．Vermeer P A，Luger H J．Rotterdam：A A Balkema，1982：355-370．

[12]Cundall P A，Hart R D．Numerical modeling of discontinua[J]．Engineering Computations，1992，9（2）：101-113．

[13]Cundall PA，JenkinsJT，IshibashiI．Evolutionofelasticmoduliinadeforminggranularassembly[C]//PowdersandGrains．BiarezJ，Gourvés R．Rotterdam:A A Balkema，1989：319-322．

[14]CundallP A，KonietzkyH，PotyondyDO．PFCeinneueswerkzeugfürnumerischemodellierungen[J]．Bautechnik，1996，73（8）：492-498．

[15]CundallP A．Rockstabilityandlong-termdamageassessment inthemine-bytunnelatURL，Canada[C]//Proceedingsof theInternationalWorkshop on Reversibility:Scientific and Technical Bases for the Reversibility of Geological Disposal(Paris，November 1998)，2nd Session，Chatenay-Malabry，France：ANDRA，1998：1．

[16]Cundall P A，Potyondy D O，Lee C A．Micromechanics-based models for fracture and breakout around the mine-by tunnel[C]//Proceedings of theExcavationDisturbedZoneWorkshoponDesigningtheExcavationDisturbedZoneforaNuclearRepositoryinHardRock(Winnipeg，September 1996)．Martino J B，Martin C D．Toronto:1996 Int Conf on Deep Geological Disposal of Radioactive Waste，1996：113-122．

[17]CundallPA，StrackODL．Thedevelopmentofconstitutivelawsforsoil usingthedistinct elementmethod[C]//NumericalMethodsinGeomechanics，Vol 1．Rotterdam：A A Balkema，1979：289-317．

[18]魏星．云南高海拔地区沥青混合料性能仿真研究[D]．重庆：重庆交通大学，2009．

[19]田莉．基于离散元方法的沥青混合料劲度模量虚拟试验研究[D]．西安：长安大学，2008．

[20]常明丰．多孔沥青路面微观力学特性与空隙衰变行为研究[D]．西安：长安大学，2009．

[责任编辑 杨屹]

Preliminary discrete element analysis on mesostructure of open-graded large stone asphalt mixes(OLSM-25)

GUO Hongbing1,2，ZHAO Yalan1,2，CHEN Shuanfa3

(1.Department of Highway Engineering,Shaanxi College of Communication Technology,Shaanxi Xi'an 710018,China;2.School ofHighway,Chang'anUniversity,Shaanxi Xi'an710064,China;3.School ofMaterialScienceandEngineering,Chang'anUniversity, Shaanxi Xi'an 710064,China)

Asphalt mixture is a kind of multiphase composite material composed of granular cells such as asphalt,coarse aggregate and fine aggregate,in the existing structural design theory of asphalt pavement at present,asphalt pavement structure has been taken as an uniform continuous,isotropic and elastic body to be analysed and calculated,which is inconsistent with the fact that asphalt mixture belongs to the discontinuous particles microstructure.In this paper,by using of the discrete element method and the two-dimensional particle flow code(PFC2D),comparative analysis on the characteristics and its applicable conditions have been done for two kinds of particles generation methods,the upper limit of the recommended gradation range of Open-graded Large Stone Asphalt Mixes(OLSM-25)has been taken as 1#gradations, the two-dimensional discrete element model of 1#gradation of OLSM-25 mixes has been established,therfore,it verified that the two-dimensional discrete element method is feasible for the numerical simulation of asphalt mixes composed of particle materials,and made up the deficiencies of finite element method,all of these provide a new idea and solution to do analysis on mesostructure of OLSM-25.

road engineering;open-gradedlarge stone asphaltmixes;mesostructure;discreteelementmethod;numerical analysis

U414.1

A

1007-2373(2015)05-0090-07

10.14081/j.cnki.hgdxb.2015.05.019

2014-12-17

陕西省教育厅专项科研计划（14JK1067）；中央高校基本科研业务费专项资金（CHD 2010JC011）

郭红兵（1976-），男（汉族），副教授．通信作者：赵亚兰（1979-），女（汉族），副教授，329639783@qq.com.