《沥青混合料离散颗粒流数值模拟与方法》系统地讲述了离散颗粒流软件PFC2D/3D的基础理论方法以及实例分析介绍，主要内容包括：颗粒流程序PFC2D/3D介绍，集料成型与级配骨架理论实例分析，沥青混合料疲劳失效实例分析，沥青混合料永久变形实例分析，胶粉改性沥青混合料断裂实例分析以及双层排水沥青混合料空隙堵塞与衰变实例分析。《沥青混合料离散颗粒流数值模拟与方法》*次系统且全面介绍了离散颗粒流在道路工程领域中的应用原理与技术要点。不仅阐明了各类道路材料的细观力学机制，也针对沥青路面各类病害提出针对性的跨尺度分析方法。通过《沥青混合料离散颗粒流数值模拟与方法》的学习，读者可以从宏观尺度更加深刻地理解道路材料力学强度形成规律，从而掌握相应的细观力学分析方法并将其应用于专业科研。

目录
前言
第1章　绪论 1
1.1　离散颗粒流的发展与应用 1
1.1.1　数值计算方法与多尺度划分 1
1.1.2　离散颗粒流基本思想 4
1.1.3　离散颗粒流的数值计算软件 9
1.2　离散颗粒流实际工程应用可行性 11
1.2.1　离散颗粒流在岩土工程中的应用 11
1.2.2　离散颗粒流在水利工程中的应用 12
1.2.3　离散颗粒流在农业工程中的应用 14
1.2.4　离散颗粒流在道路工程中的应用 15
参考文献 16
第2章　颗粒流程序PFC2D/3D 17
2.1　离散元方法与颗粒流程序PFC 17
2.1.1　颗粒流程序PFC的基本元素 17
2.1.2　颗粒流程序的特点与优势 20
2.1.3　颗粒流程序运行的计算机硬件要求 21
2.2　颗粒流程序PFC的组成与工具 22
2.2.1　PFC的工具组成 22
2.2.2　PFC的常用菜单 24
2.2.3　PFC的常用命令 26
2.2.4　PFC的文件类型 30
2.3　颗粒流程序PFC的Fish语言 32
2.3.1　Fish语言规则、变量和函数 32
2.3.2　Fish语句 35
2.3.3　Fish和PFC的关联 37
2.3.4　PFC的帮助文档 39
参考文献 39
第3章　基于离散颗粒流的集料成型与级配骨架理论研究 41
3.1　集料形貌扫描及数值成型方法 41
3.1.1　AIMS骨料扫描技术 41
3.1.2　基于半径扩大法和内部填充法的集料简化二维形态成型方法 46
3.1.3　基于棱角度指标的集料真实二维形态成型方法 47
3.1.4　基于随机算法的集料不规则三维形态成型方法 49
3.2　沥青混合料二维集料骨架力学响应数值模拟分析 52
3.2.1　单一、组合粒径集料单轴贯入试验数值模拟 52
3.2.2　二维骨架成型方法细观接触力状态分析 54
3.2.3　二维骨架成型方法计算效率对比分析 57
3.3　沥青混合料三维集料骨架理论数值模拟分析 58
3.3.1　沥青混合料局部骨架理论 58
3.3.2　沥青混合料全局骨架理论 62
3.3.3　沥青混合料级配骨架理论验证 64
3.3.4　三维骨架级配细观接触力学特性分析 67
3.3.5　三维骨架级配配位数与力链形成机制 71
参考文献 78
第4章　基于离散颗粒流的沥青混合料疲劳失效研究 79
4.1　沥青混合料数值成型方法 79
4.1.1　沥青混合料二维数值试件级配的确定 79
4.1.2　沥青混合料二维数值试件成型 83
4.1.3　沥青混合料三维数值试件成型 89
4.2　沥青混合料虚拟疲劳试验数值模拟 93
4.2.1　疲劳损伤离散元表征 93
4.2.2　四点小梁弯*疲劳试验数值模拟 97
4.2.3　虚拟疲劳数值试验验证 102
4.3　沥青混合料细观结构对疲劳性能的影响 106
4.3.1　集料宏观特征影响性分析 106
4.3.2　集料细观特征影响性分析 112
4.3.3　空隙结构特征影响性分析 118
4.4　移动荷载下沥青路面材料疲劳特征多尺度数值模拟 121
4.4.1　移动荷载下路面结构响应分析 121
4.4.2　路面结构层模量和厚度影响性分析 131
4.4.3　路面结构层间黏结影响性分析 134
4.4.4　移动荷载大小和速度影响性分析 136
4.4.5　移动荷载间歇时间影响性分析 138
参考文献 139
第5章　基于离散颗粒流的沥青混合料永久变形研究 140
5.1　力学-经验车辙预测模型 140
5.1.1　基于应力-应变的两阶段车辙预估模型 140
5.1.2　基于黏弹性回弹应变的三阶段车辙预估模型 142
5.2　车辙试验数值模拟 144
5.2.1　数值车辙试件成型与分析 144
5.2.2　数值车辙试件内部接触力分布 148
5.2.3　虚拟车辙试件内部位移分布 151
5.2.4　粗集料空间运动轨迹变化 155
5.2.5　粗集料转角变化 158
5.2.6　粗集料接触状态变化 160
5.3　沥青混合料结构对永久变形性能的影响 161
5.3.1　集料宏观特征影响性分析 161
5.3.2　集料细观特征影响性分析 166
5.3.3　空隙结构特征影响性分析 172
参考文献 178
第6章　基于离散颗粒流的胶粉改性沥青混合料断裂研究 179
6.1　室内断裂试验设计 179
6.1.1　胶粉改性沥青断裂试验 179
6.1.2　胶粉改性沥青砂浆断裂试验 184
6.1.3　胶粉改性沥青混合料断裂试验 188
6.2　胶粉改性沥青及混合料数值成型方法 191
6.2.1　边缘识别分水岭算法 191
6.2.2　胶粉颗粒随机分布算法 194
6.2.3　本构添加与非均匀赋参 195
6.3　胶粉改性沥青断裂数值模拟分析 200
6.3.1　宏观断裂响应行为特性分析 200
6.3.2　虚拟断裂过渡区演化分析 203
6.3.3　微裂纹扩展与胶粉增韧机制 204
6.3.4　宏观断裂能与微观耗散能分析 209
6.4　胶粉改性沥青砂浆断裂数值模拟分析 211
6.4.1　宏观断裂响应行为特性分析 211
6.4.2　微裂纹扩展与胶粉增韧机制 213
6.4.3　宏观断裂能与微观耗散能分析 216
6.4.4　断口形貌与失效机制分析 218
6.5　胶粉改性沥青混合料断裂数值模拟分析 226
6.5.1　DC(T)数值模拟 226
6.5.2　I-FIT数值模拟 228
6.5.3　IDEAL-CT数值模拟 232
6.5.4　界面黏结强度影响性分析 235
6.5.5　空隙率影响性分析 238
参考文献 241
第7章　基于离散颗粒流的双层排水沥青混合料空隙堵塞与衰变研究 242
7.1　双层排水沥青混合料空隙结构仿真 242
7.1.1　空隙结构提取 242
7.1.2　空隙结构分类 243
7.1.3　空隙结构指标 245
7.2　双层排水沥青混合料空隙结构形成数值模拟分析 248
7.2.1　空隙率与连通空隙率 248
7.2.2　公称*大粒径对空隙结构形成的影响 248
7.2.3　干扰粒径集料对空隙结构形成的影响 252
7.2.4　空隙结构形成机制验证 254
7.3　双层排水沥青混合料空隙结构衰变规律数值模拟分析 255
7.3.1　排水沥青混合料车辙试验数值模拟 255
7.3.2　室内车辙试验对比验证 257
7.3.3　空隙结构衰变规律 258
7.4　双层排水沥青路面空隙堵塞数值模拟分析 260
7.4.1　双层排水沥青路面结构成型 260
7.4.2　堵塞物模型构建 261
7.4.3　重力和流水作用下的空隙堵塞试验数值模拟 263
参考文献 268

第1章　绪论
　　现今，我国的道路建设和技术研究正处于一个蓬勃发展的时期，但沥青路面的各类病害，仍然对车辆行驶的安全性和舒适性造成了很大影响。当前国内外学者对于沥青混合料的研究主要聚焦于沥青混合料的宏观和物理特性等室内试验方面，而从微细观的角度对沥青混合料的内部组成结构及力学特性的探究相对薄弱。因此本章主要对道路工程沥青及其混合料的研究现状以及离散颗粒流的发展与实际工程应用进行简要的介绍，以便在未来的研究中应用先进的计算手段、力学分析工具，定量地、多尺度地分析沥青混合料的力学特征，这对于深刻理解沥青混合料的力学机制、改进沥青混合料设计方法，具有十分重要的意义。
　　？1.1　离散颗粒流的发展与应用
　　？1.1.1　数值计算方法与多尺度划分
　　当对道路工程中的常见材料进行弹性简化后，路面结构可以被视为半无限空间体上的层状弹性体系，如图1-1所示，在此基础上，运用弹性力学理论可以获取在双圆荷载作用下各层路面结构的受力及变形状态，E为各层模量，h为各层厚度， 为各层泊松比，*后根据与材料、结构性能相关的经验公式验算结构的服役性能。以上步骤是当前路面结构力学-经验设计的主要流程。然而，路面材料的力学行为、路面结构的服役状态与层状弹性体系理论偏离较大，受诸多因素影响，主要表现在以下几个方面。
　　（1）材料本构的差异：道路材料多种多样，包括沥青混合料、水泥混凝土、无机结合料稳定类材料、土、岩石、级配碎石等，简单的弹性力学本构无法充分表征材料的力学行为，例如，路面车辙与沥青混合料的黏塑性变形、粒料类材料的塑性变形相关。
　　（2）材料组成的差异：宏观力学往往建立在均匀性的假设条件之上，但实际工程材料并非为单一介质材料，而是由多种介质组成的复合材料，变异性较高。
　　（3）荷载差异：从几何角度分析，车轮-路面之间的接触形状与接触力分布并非规则的圆形荷载，从时间角度分析，车轮-路面之间接触的应力值、面积和车辆速度关联紧密；此外，路面结构还将承受温湿度变化下的冻融、干湿循环等作用，同时部分病害处将产生动水冲刷作用等。
　　图1-1　弹性层状体系理论计算示意图
　　（4）边界差异：路面结构往往由于既有病害，如路基不均匀沉降、裂缝、层间黏结失效等，使得结构的支承条件与边界条件发生复杂变化。
　　实际路面结构受力分析影响因素如图1-2所示。
　　图1-2　实际路面结构受力分析影响因素
　　基于以上特点，道路工程中的力学问题无法直接使用简单的物理-力学模型进行描述，更不可能得到数学上的解析解，因此，需要使用数值计算的方法对复杂工程问题进行仿真模拟。在数值计算中，沥青路面及材料可以分为以下四个尺度，分别为沥青胶浆（asphalt mastic）尺度、沥青砂浆（asphalt mortar）尺度、沥青混合料（asphalt mixture）尺度以及沥青路面（asphalt pavement）尺度。其中，沥青胶浆尺度是在沥青中加入直径数量级在10？6～5×10？5m的石灰粉、矿粉等材料；沥青砂浆尺度是在沥青胶浆中加入直径数量级在7.5×10？5～10？3m的细集料等材料；沥青混合料尺度是在沥青砂浆中加入直径数量级在2×10？3～5×10？2m的粗集料等材料；沥青路面尺度则是将不同种类的沥青混合料加入路面结构中进行设计。掌握并理解不同尺度下的沥青路面及其材料的宏观与微（细）观变形特性，有助于沥青路面材料性能及结构的优化设计，目前该类研究的数值分析手段主要包括有限单元法、有限差分法、边界元法等。
　　？1. 有限单元法
　　有限单元法是数值分析方法的一种，其基本思想是将连续的求解区域离散为一组有限个且按一定方式相互联结在一起的单元组合体。单元能按不同的联结方式进行组合，且单元本身又可以有不同的形状，因此可以模拟几何形态复杂的求解域。有限单元法作为数值分析方法的另一个重要特点，是利用每一个单元内假设的近似函数来分片地表述全求解域上待求的未知场函数。单元内的近似函数，通常由未知场函数以及各个节点的数值来表达。如此针对一个问题的有限元分析，未知场函数及其导数在各个节点上的数值就成为新的未知量（即自由度），从而使一个连续的无限自由度问题转化为离散的有限自由度求解。当单元尺寸满足收敛要求时，数值计算的近似解将收敛于精确解。
　　？2. 有限差分法
　　有限差分法是计算机数值模拟*早采用的方法之一，其以Taylor级数展开等方式，将控制方程中的导数用网格节点上函数值的一阶微分代替进行离散，从而建立以网格节点值为未知数的代数方程组。该方法概念直观，表达简单，其与有限单元法的差异主要体现在以下两个方面：①在网格划分上，有限元网格划分不需要规则，各种单元可以混合使用，而有限差分法划分的网格必须是规则的，从而对方程进行离散化；②在计算效率上，虽然有限单元法与有限差分法都是求解一组方程组，但有限单元法通常采用隐式、矩阵解算方法，而有限差分法采用显式、时间递步法解算代数方程，一般而言，有限差分法的计算效率要高于有限单元法，计算量小于有限单元法。
　　？3. 边界元法
　　边界元法是在有限单元法后发展起来的一种较精确、有效的工程数值分析方法，它以定义在边界上的边界积分方程为控制方程，通过对边界元进行插值离散，化为代数方程组求解。其与基于偏微分方程的区域解法相比，由于降低了问题的维数，从而显著地降低了自由度数，边界的离散也比区域的离散方便得多。特别是对于边界变量变化梯度较大的问题，如应力集中问题或边界变量出现奇异性的裂纹问题，边界元法公认比有限单元法更精确高效。边界元法所利用的微分算子基本解能够自动满足无限远处的条件，因此边界元法特别便于处理无限域和半无限域问题。
　　？1.1.2　离散颗粒流基本思想
　　弹性力学的方法要求研究对象满足连续性、均匀性、各向同性、线弹性、小变形五大假设，其通过平衡微分方程、物理方程、几何方程及边界条件所求解的结构受力、变形问题的结果是*为精确的。其对于复杂的结构却很难给出微分方程组的解析解，虽然运用变分法可以巧妙地将微分求解转换为积分求解，但结果的准确性仍与假定的变形假设有关。有限单元法在此基础上，通过将求解对象的网格离散为单元，基于本构模型构建刚度矩阵，可以得到力和位移的方程组，再借助边界条件开展求解。然而在利用有限元求解一系列非线性问题时，单元之间必须连续且保持变形协调。 与有限单元法不同，近年来，新兴的离散元方法（DEM）允许各单元在运动过程中相互分离，因此在处理非均质颗粒体系以及结构大变形上具有一定的优势。离散元方法的思想源于较早的分子动力学，是研究不连续体力学行为的一种新数值方法。其基本思想是把散粒群体简化成具有一定形状和质量颗粒的集合，赋予接触颗粒间及颗粒与接触边界（机械部件）间某种接触力学模型和模型中的参数，以考虑颗粒之间及颗粒与边界间的接触作用和散粒体与边界的不同物理机械性质。离散元方法采用牛顿第二定律、动态松弛法和时步迭代求解每个颗粒的运动速度和位移，因而特别适用于求解非线性问题。其特点是在分析高度复杂的系统时，无论颗粒还是边界均不需要做大的简化；当赋予接触颗粒间不同的接触模型时，还可以分析颗粒结块、颗粒群聚合体的破坏过程、多相流动甚至可以包括化学反应和传热等问题。离散元方法的一般求解过程为将求解空间离散为离散元单元阵，并根据实际问题用合理的连接元件将相邻两单元连接起来；单元间相对位移是基本变量，由力与相对位移的关系可得到两单元间法向和切向的作用力；对单元各方向所受到的作用力与其他单元间的作用力以及其他物理场对单元作用所引起的外力求合力和合力矩，根据牛顿第二定律可以求得单元的加速度；对其进行时间积分，进而得到单元的速度和位移，从而得到所有单元在任意时刻的速度、加速度、角速度、线位移和转角等物理量。
　　1. 离散元方法的基本方程
　　连续介质问题中需要满足的方程为平衡方程、变形协调方程（几何方程）、物理方程（本构方程）以及边界条件。离散元方法一开始就将对象离散化，并不赋予连续性假设，则变形协调方程可以不满足。为了进一步简化计算，离散元方法引入如下假设：①单元是刚性体；②单元间的接触范围十分小；③单元间的接触为柔性接触，因此允许单元接触时出现一定的“重叠”，但重叠量相对于单元的尺寸很小。由以上三点基本假设可知，在离散元方法内单元自身是不允许变形的，单元间的接触力与“重叠”的量有关且为柔性接触。这意味着本构方程及边界条件只需在边界接触上成立，且必须满足单元的平衡方程，柔性接触的变形量以及单元间接触力的大小则由细观接触模型完全确定。综上可知，离散元方法的基本方程应包括接触本构方程、平衡方程、运动方程这三类。
　　1）接触本构方程
　　接触本构方程即代表接触力与位移之间的关系，通过接触模型的选择，可以模拟单元间弹性、塑性、黏性等力学行为，离散元内的接触本构方程形式多样，包括接触刚度、接触滑动、伯格斯模型等。
　　2）平衡方程
　　在离散元中，由于单元是运动的，其平衡方程为动力平衡方程。平衡方程包括力的平衡方程以及力矩的平衡方程，具体如下。采用式（1-1）～式（1-3）计算各个单元的不平衡力及不平衡力矩，并将这些结果代入运动方程中求解。
　　（1-1）
　　（1-2）
　　（1-3）
　　式中， 为体力在x方向上的分量； 为体力在y方向上的分量； 为第i个接触力在x方向上的分量； 为第i个接触力在y方向上的分量； 为不平衡力在x方向上的分量； 为不平衡力在y方向上的分量； 为接触数量； 为第i个接触力的大小； 为第i个接触力到质心的距离； 为不平衡力矩。
　　3）运动方程
　　运动状态包括平动状态和转动状态。以二维状态下的x方向为例，设在时间t0时，单元的质量为m，转动惯量为I，运动方程如下：
　　（1-4）
　　（1-5）
　　已知 状态下的运动状态，采用的是中心有限差分的方法向前推进。因此，在 时，颗粒在x方向的平动速度和转动速度如式（1-6）和式（1-7）所示：
　　（1-6）
　　（1-7）
　　在 时，颗粒在x方向的位移为
　　（1-8）
　　计算到这一步时，接触力已进行更新，相应的不平衡力以及不平衡力矩也已分别更新为 、 ，故有
　　（1-9）
　　（1-10）
　　2. 离散元方法的求解方法
　　在离散元基本原理方程的基础上，以PFC5.0为例，介绍在一个时间步循环内各方程的求解过程，总体求解步骤如图1-3所示。
　　1）确定时间步长
　　为了使上述牛顿时间步的显式解稳定，则要求所取的时间步小于等于全局时间步的*小值。一般而言，*小时间步可以通过遍历模型所有单元，考虑运动速度、加速度的*不利形式组合，并通过约束的运动位移反算求得。以平动为例，利用牛顿第二定律，计算平动加速度的一个分量：
　　（1-11）