利用分子尺度动态建模方法表征复杂聚合物材料沥青的粘弹性和疲劳损伤特性

《Advanced Industrial and Engineering Polymer Research》：Characterizing Viscoelasticity and Fatigue Damage of Complex Polymeric Material Asphalt Using Molecular-scale Dynamic Modelling

　　田伟伟|李辉|张玉清
东南大学交通学院，中国江苏省南京市211189

**摘要**
采用分子动力学模拟方法研究了基于聚合物的复合材料——沥青粘合剂的纳米尺度特性。然而，纳米尺度下的非线性粘弹性及疲劳损伤机制仍不甚明了，且分子模拟参数与宏观力学性能之间缺乏可靠的定量关联。由于缺乏能够捕捉关键参数（如动态模量）的分子尺度方法，沥青力学性能的跨尺度预测进展受到严重限制。本研究开发了一种动态加载建模算法，用于表征沥青材料在纳米尺度上的粘弹性和疲劳损伤行为。结果表明，该动态建模能够获得具有特征性粘弹性滞后行为的稳定应力-应变数据。模拟得到的动态剪切模量曲线能够反映不同分子组分之间的差异，并再现实验观察到的相对大小。尽管沥青粘合剂模型与实际材料在时间和空间尺度上存在九个数量级的差异，但高频模拟结果与SHRP报告中的宏观数据偏差小于0.5个数量级。在恒定温度和加载频率下，沥青分子模型表现出典型的力学响应阶段：线性粘弹性、非线性粘弹性和疲劳损伤。在AAM-1沥青粘合剂模型中，0.4 GPa和0.7 GPa的输入应力幅度分别对应这些阶段的阈值。在粘弹性阶段，能量耗散基本保持恒定；而在损伤阶段，能量耗散增加，同时伴随着微观结构的变化，如孔隙生长和自由体积分数的增加，这反映了结构的逐渐退化及变形抵抗力的降低。这些发现加深了对沥青粘合剂分子尺度动力学特性的理解，为特定阶段的粘弹性和损伤行为提供了定量洞察，并为宏观力学性能的跨尺度本构建模和预测奠定了基础。

**1. 引言**
沥青混凝土是道路工程中最广泛使用的材料之一[1]。沥青粘合剂在决定沥青混合料的结构稳定性、抗裂性和耐久性方面起着关键作用。研究人员通过流变分析、化学表征和形态分析等方法对沥青粘合剂的性能和内在机制进行了大量系统研究[2][3]。随着计算技术和材料科学的进步，分子动力学（MD）模拟已成为研究沥青材料的强大工具[4]。该方法在分析沥青材料的微观结构、物理化学性质和热力学特性方面展现了巨大潜力[5][6][7][8]。

分子动力学模拟能够计算分子层面的关键参数。在现有研究中，这些参数通常作为宏观性能的间接指标使用（见表1）。例如，扩散系数和松弛时间常用于近似流变行为，而结合能则用于估算沥青混合料中的界面粘结强度。然而，这些参数作为替代指标缺乏与宏观工程性能的直接物理关联。因此，建立微观模拟参数与宏观性能之间的定量关系仍然具有挑战性，导致当前分析主要停留在定性比较阶段，而无法进行严格的定量建模。此外，作为工程材料，沥青粘合剂的力学响应是其实际应用中性能和耐久性的关键因素。从物理角度来看，扩散系数和结合能等模拟结果并不适合表征沥青粘合剂及其混合料的力学响应。因此，缺乏直接计算沥青模型力学参数的分子动力学模拟方法，使得难以建立分子尺度材料属性与宏观力学性能之间的直接联系。

**表1. 现有研究中从MD模拟参数推断出的宏观性能**
| MD模拟参数 | 现有研究 | 宏观材料性能 |
|-----------------|-----------------|--------------------------------------------|
| 扩散系数 | You等人[9] | 变形抵抗力 |
| 粘度 | Ding等人[10] | |
| 自修复性能 | Cao等人[11] | |
| 耐老化性能 | Yaphary等人[8] | |
| 低温抗裂性 | Zadshir等人[12] | |
| 力学强度 | Sun等人[13] | |
| 热稳定性 | Du等人[14] | |
| 界面强度 | Gao等人[15] | |
| 储存稳定性 | Cui等人[16] | |
| 结合能密度 | Huang等人[17] | |
| 兼容性 | Xu等人[18] | |
| 黏附性 | Ding等人[19] | |
| 分子极性 | Meghwal等人[20] | |
| 自修复性能 | He等人[22] | |
| 松弛时间 | Sonibare等人[7] | |
| 粘度 | | |

**2. 方法论**
Li和Greenfield基于战略公路研究计划（SHRP）数据库[29]开发了三种典型沥青类型（AAA-1、AAK-1和AAM-1）的分子模型，并对其有效性进行了彻底验证。这些模型是目前分子模拟研究中最常用的沥青粘合剂模型。作为国际公认的核心参考数据库，SHRP包含了来自不同气候区和性能水平的样本。AAA-1、AAK-1和AAM-1模型能够反映沥青粘合剂的差异，因此非常适合用于验证本研究中模拟方法的一般适用性。选择这些模型有助于将模拟结果与宏观实验数据进行比较，提升方法的工程相关性。此外，这三种模型的广泛应用为将本研究的方法与以往工作联系起来提供了有效基础，为未来的优化和扩展提供了通用数据支持。

为了减少模型尺度对模拟结果的影响[30][31]，模型中的分子数量已超出Li和Greenfield最初提出的模型[29]。为了评估分子尺度动态建模算法对组成变化的敏感性，构建了SARA组分的代表性分子模型，包括沥青质（As）、树脂（Re）、芳香烃（Ar）和饱和烃（Sa）。此外，还构建了不同尺度下的AAM-1沥青模型，以分析模型尺度和输入应力幅度对动态模量模拟结果稳定性的影响。所有模型均采用立方体几何结构。具体模型组成见表2，其中×n表示与Li和Greenfield最初提出的模型相比的分子数量扩展因子。表3总结了上述分析中使用的所有与粘弹性相关的模拟变量。

**表2. 沥青粘合剂及四组分模型的分子组成**
| 分子名称 | 化学配方 | 分子组成 |
|-----------------|-----------------|----------------------------------|
| AAA-1 (×10) | C42H54O30 | |
| AAK-1 (×10) | C66H81N20 | |
| AAM-1 (×n, n=1, 2, 4, 6, 8, 10) | C36H57N40 | |
| | C29H50O50 | |
| | C40H60S40 | |
| | C35H44 | |
| | C30H6240 | |
| | C40H59N40 | |
| | C18H10S2150 | |
| | C30H46 | |
| | C30H6240 | |
| | C35H6240 | |
| | C40H59N40 | |
| | C18H10S2150 | |
| | C30H46 | |
| | C30H6240 | |
| | C40H50O50 | |
| | C40H60S40 | |
| | C40H59N40 | |
| | C18H10S2150 | |
| | C30H44 | |
| | C30H6240 | |
| | C40H50O50 | |
| | C40H60S40 | |
| | C40H59N40 | |
| | C18H10S2150 | |
| | C30H46 | |
| | C30H6240 | |
| | C40H50O50 | |
| | C40H60S40 | |
| | C40H59N40 | |
| | C18H10S2150 | |
| | C30H44 | |
| | C30H6240 | |
| | C40H50O50 | |
| | C40H60S40 | |
| | C40H59N40 | |
| | C18H10S2150 | |
| | C30H46 | |
| | C30H6240 | |
| | C40H50O50 | |
| | C40H60S40 | |
| | C40H59N40 | |
| | C18H10S2150 | |
| | C30H44 | |
| | C30H6240 | |
| | C40H50O50 | |
| | C40H60S40 | |
| | C40H59N40 | |
| | C18H10S2150 | |
| | C30H46 | |
| | C30H6240 | |
| | C40H50O50 | |
| | C40H60S40 | |
| | C40H59N40 | |
| | C18H10S2150 | |
| | C30H44 | |
| | C30H6240 | |
| | C40H50O50 | |
| | C40H60S40 | |
| | C40H59N40 | |
| | C18H10S2150 | |
| | C30H46 | |
| | C30H6240 | |
| | C40H50O50 | |
| | C40H60S40 | |

**表3. 粘弹性模拟设置、加载参数和响应指标**
| 模拟集 | 变量类别 | 参数水平 | 响应指标 |
|-----------------|-----------------|----------------------------------|
| Set Ⅰ | 温度，°C | 25-185（步长40） | 动态模量 |
| | | 417, 208, 139, 104, 83, 60, 42, 30, 19 | |
| | | 2.4, 4.8, 7.2, 9.6, 12, 16.8, 24, 33.6, 52.8 | |
| | | 0.2 | |
| | | 10 | |
| | | 2400 | |
| Set Ⅱ | 温度，°C | 25 | 动态模量&耗散能量 |
| | | 208 | |
| | | 0.1–1.0（步长0.1） | |
| | | 500 | |
| | | 2400 | |
| Set Ⅲ | 模型尺度 | AAM-1×2–AAM-1×10（步长×2） | 动态模量 |
| | | 0.1–1.0（步长0.1） | |
| | | 25 | |
| | | 208 | |
| | | 4.8 | |
| | | 50 | |
| | | 240 | |
| | | 240 | |

总体而言，扩大模型尺度可以提高模拟结果的统计稳定性，但会显著增加计算成本和运行时间。根据可用计算资源和研究目标，设计了三种不同的模拟集来验证所提出的动态加载模拟方法的有效性。在表3的Set Ⅰ中，主要使用了大尺度的AAA-1×10、AAK-1×10、AAM-1×10及四组分沥青模型（As、Re、Ar、Sa），在低应力水平、不同的温度-频率组合和有限的加载循环次数（10次循环）下进行动态剪切加载模拟，以确保计算效率。动态剪切模量是在10次加载循环后平均得到的，并将不同温度下的频率依赖性模量曲线沿频率轴移动并拟合成主曲线，以评估时间-温度叠加原理在分子尺度上的适用性，并比较不同模型的响应。

另一方面，在表3的Set Ⅱ中，使用小尺度的AAM-1×1模型在500次加载循环下计算动态力学响应。动态模量和能量耗散的演变被分析出来，以表征在不同输入应力幅度下的粘弹性行为和疲劳损伤进展。在表3的第三组中，不同的AAM-1沥青模型在0.1到1 GPa的输入应力幅度下经历了50个周期的动态剪切加载。这项研究的目的是调查模型尺度和应力幅度对动态模量模拟结果稳定性的耦合效应。在本研究中，模型的粘弹性响应自然地来源于所采用的计算力场控制的原子间相互作用。关键的粘弹性特征，如应力-应变滞后，本质上源自动态机械加载下的分子动力学。没有明确施加或规定任何粘弹性本构模型来定义机械响应。影响粘弹性响应的关键变量包括温度、加载频率（或周期）、施加的应力幅度和总加载时间。在加载过程中，使用Nose-Hoover-Langevin（NHL）恒温器将温度保持在目标值。通过改变每个正弦周期内的应力增量次数（每个应力增量模拟0.2 ps）来调整加载频率。正弦应力是使用Souza-Martins应力控制方法施加的，控制原理将在后续部分详细说明。实际应力响应是从体积平均原子应力张量（维里应力）计算得出的，然后用于计算动态模量。模型的空间几何结构，主要是立方尺寸和三维示意图，以及动态加载模拟的总体框架，在图1中进行了说明。

图1. 动态加载模拟的工作流程。

2.2. 初始模型构建和平衡
本研究采用了COMPASS力场，该力场已经针对小分子、聚合物、沥青和某些无机系统进行了彻底验证[32]，确保所提出的动态加载模拟方法可以可靠地扩展到其他材料系统。COMPASS最初由Sun等人开发[33]，并经过持续改进，形成了一个成熟的参数集，并不断更新以提高计算精度。本研究中施加的应力幅度（0.1 GPa-1 GPa）完全在COMPASS的传统范围内，远低于之前报道的极端应力水平。例如，Chowdhury等人在纤维模拟中施加了35 GPa[34]，而Zaminpayma等人在复合材料中使用了15 GPa[35]。在整个研究中，没有观察到力场失效的迹象，如能量发散或异常原子位移。COMPASS力场通过授权的Materials Studio提供，其核心参数被加密，只能在软件内部访问。为了确保透明度和可重复性，开放和加密的参数文件都已上传到公共仓库（链接如下）。COMPASS适用于常规应力水平和短时间尺度，但无法捕捉超高应力下的键断裂（例如，冲击加载）或有效处理长时间尺度模拟。因此，在进行超出常规条件的模拟时需要仔细评估。

采用周期性边界条件以消除边界效应。采用具有平移对称性的P1空间群，以避免引入额外的假设对称性。分别使用Lennard-Jones势和Coulomb势来描述粒子之间的范德华力和静电相互作用。根据表2中的分子组成，构建了初始立方晶胞模型，初始密度为0.7 g/cm3。模型使用Smart算法进行了100,000次几何优化步骤，然后在300K到650K的温度范围内进行了2.5 ns的退火循环，以减少不合理的局部结构。在特定温度下进行了初步的1 ns动态松弛，随后进行了5 ns的分子动力学平衡。上述所有过程都在等温-等压系综（NPT系综）中进行，使用Nosé-Hoover恒温器和Berendsen压强计分别控制温度和压力。平衡阶段的最终帧被作为后续研究的结构，代表一个充分平衡的模型。由于空间限制，这里没有包括额外的验证分析。先前的研究表明，该模型在相似的松弛和平衡时间内可以达到充分平衡的状态[31]、[36]、[37]、[38]。

2.3. 由目标模拟温度驱动的初始原子速度配置
在动态加载之前，初始化原子速度以确保分子模型在整个模拟过程中达到并保持目标温度。每个原子的初始速度基于Maxwell-Boltzmann分布[39]分配。首先，在[-1, 1]范围内随机生成两个实数，并均匀分布。这些值被检查以满足方程（1）的要求。(1)
根据生成的值和，计算一个转换因子，将均匀分布的随机数转换为标准正态分布值。修正后的玻尔兹曼常数显示在方程（4）中。(2)(3)(4)
其中是玻尔兹曼常数1.380658×10-23 J/K，是对应于1 amu的质量，1.6605402×10-27 kg，是将?/ps转换为m/s的转换因子，其值为102。
接下来，计算温度转换因子，并确定每个原子在X、Y和Z方向上的速度分量。(5)(6)
其中是目标温度，是原子的质量，和，代表遵循标准正态分布的随机数。
随机分配的原子速度可能会引入非零的总动量，导致不必要的平移漂移。为了消除这种效应，计算质心速度并从每个原子中减去，确保系统的总动量为零。(7)(8)(9)
其中代表所有原子在方向上的总动量，而表示系统中的原子总数。变量对应于原子的质量，代表原子在方向上的速度分量。项表示在方向上的质心速度，表示在方向上经过动量调整后的原子速度分量。
由于质心运动不贡献于系统的内部能量，因此可以使用方程（10）计算包含个原子的系统的有效自由度，减去与平移运动相关的三个自由度。系统的总动能是所有原子动能的总和，由质量和速度决定。然后根据气体动力学理论中的总动能和温度关系计算系统的绝对温度。(10)(11)(12)(13)
其中代表系统的总动能，是原子总数，表示原子的质量。原子在三个空间方向上的速度分量分别由给出。系统的有效自由度表示为，而表示修改后的玻尔兹曼常数，是系统的绝对温度。
接下来，通过缩放速度将系统温度调整到目标温度。使用方程（14）计算温度缩放因子，并使用方程（15）相应地更新系统中所有原子的速度分量。(14)(15)
其中表示温度缩放因子，表示不同方向上原子的更新后的速度分量。

2.4. 应力控制的动态加载和应力-应变响应数据收集方法
每个周期的输入应力由方程（16）定义。每次应力调整后进行200个动态模拟步骤，前100步用于系统平衡，剩余100步用于数据收集。时间步长为1 fs。加载周期由每个周期内的应力变化次数控制。(16)
其中表示输入应力值，表示应力输入的幅度，表示当前应力值在一个周期内的索引，表示应力值根据正弦曲线变化的次数。和是时间依赖函数。由于连续周期共享一个零应力点，因此可以从方程（17）确定周期周期。例如，当一个周期内的应力变化次数为61时，所施加应力的周期周期为12 ps。(17)(18)
其中表示周期周期，表示每个应力变化的系统松弛时间，表示每个应力变化的松弛步骤数，表示每个步骤的松弛时间（1 fs）。
施加外部应力的主要目的是为模型设置一个目标应力。实际上，由于分子动力学的非瞬时性和异质性，分子的实际应力响应可能会偏离这个目标。因此，需要计算模型对输入应力的实际响应。单个原子上的力可以通过应力张量来描述，其表达式和计算方法在方程（19）、（20）[40]中提供。通过沿给定方向对系统中每个原子的应力张量分量求和并除以系统总体积，得到实际响应应力，如方程（21）所示。(19)(20)(21)
其中表示原子应力张量，，等表示其在不同方向上的分量；表示原子在方向和上的应力张量分量，表示原子的质量，和是原子在方向和上的速度分量，表示原子在方向上施加的力的分量；表示模型在方向和上的宏观应力，表示系统的总体积。
基于连续介质力学理论，使用变形梯度来描述模型中的应变。根据方程（22）、（23）构建晶胞矩阵，其中变形梯度是从参考晶胞（模拟的初始模型输入）和变形晶胞（受应力加载的模型）之间的基矢量变化得出的，如方程（24）所示。使用方程（25）[41]制定应变张量。最后，使用方程（26）、（27）计算剪切或拉伸-压缩方向上的应变。(22)(23)(24)
其中表示施加应力载荷下的晶胞矩阵，而表示初始晶胞矩阵。输出模型晶胞向量在方向上的分量由给出。类似地，，，表示初始晶胞向量的分量。变形梯度表示为。(25)(26)(27)
其中表示应变张量，表示单位矩阵。平面或轴向方向上的应变分量由给出；表示平面工程应变，表示轴向工程应变。
如上所述，模拟包括一个100 fs的结构松弛阶段，然后是每个目标应力点的100 fs的生产阶段。在生产阶段，每5 fs记录一次结构快照，共得到20帧。这些帧被分成四个块，每五个帧一组。为了评估模拟结果的统计可靠性，对生产阶段收集的轨迹帧进行了块平均处理。该过程在方程（28）、（29）、（30）中进行了形式化。表示每个应力持续时间内的结构帧总数，结构轨迹被分成块。然后给出每个轨迹块内的结构帧数：(28)
其中表示每个轨迹块内的结构帧数。该值是通过应用函数获得的，以避免索引错误并确保数据均匀划分。
计算在给定应力下整个生产时间内收集的数据的平均值。方程（30）的结果是对指定输入应力值的实际响应。(29)(30)
其中表示轨迹块的数据平均值，表示块的起始帧索引，表示下一个块的起始帧索引，表示轨迹块内结构帧的收集数据。表示生产时间内收集的数据的平均值。
作为一个示例，图2显示了AAA-1×10、AAK-1×10和AAM-1×10粘合剂模型在25 °C、12 ps加载周期和0.2 GPa应力幅度下的动态剪切加载下的剪切应力响应。该图表展示了目标应力峰值（0.2 GPa）处轨迹帧的原始应力值、相应的块平均数据以及各自的变异系数（CV）。对于块平均数据，首先计算了每个块的平均响应应力，将每个值视为独立样本来计算总体平均值和变异系数。值得注意的是，尽管所有三种沥青模型的原始应力响应数据都表现出相当大的波动，变异系数范围从0.27到0.31，但块平均数据的变异性显著降低，变异系数低于0.05。结果表明，在给定的加载方法下，平均机械响应在统计上是明确且可重复的，从概率意义上反映了稳定性。

下载：下载高分辨率图像（406KB）
下载：下载全尺寸图像

图2. 使用块平均方法对剪切应力响应的统计表征。

动态模量是通过傅里叶变换[42]在基频下计算应力与应变幅度的比值得到的。

2.5. 其他模拟细节
在动态加载模拟过程中，使用Souza-Martin方法来驱动系统的瞬态响应[43]。由于周期性边界条件（PBCs）通过在空间中无限复制模拟单元来消除边界效应，因此不存在物理边界原子，也无法直接在固定边界上施加应力。Souza-Martins方法迭代调整模拟单元的变形，使得整个系统的应力接近预定的目标值。它将单元变形的虚构动力学（通过度量张量表达）与哈密顿框架内的实际原子运动结合起来，从而确保严格的能量守恒并便于精确的应力调节。这种方法还避免了传统晶格矢量控制方法中固有的关键问题，包括虚假旋转和单元变换下的不变性丧失。图3展示了Souza-Martins方法的闭环工作流程，包括目标应力输入、单元变形调整和应力-应变响应输出。由于空间限制，详细的公式可以在Souza和Martins的原始工作中找到，这里不再赘述[43]。

下载：下载高分辨率图像（908KB）
下载：下载全尺寸图像

图3. Souza-Martins方法的应力控制示意图。

在进行动态加载模拟之前，首先进行5000 fs的零应力平衡，以最小化残余应力对数据收集结果的影响。零应力平衡后，使用模型的晶胞作为参考单元。

所有模拟过程都使用自开发的Perl脚本算法执行。算法执行平台是Materials Studio 2020。

2.6. 实验验证方法
为了系统地验证本研究中提出的动态加载模拟方法的可靠性，将使用两组实验数据进行验证分析。首先，将权威SHRP-A-369报告中公开的两种沥青类型（AAK-1和AAM-1）在25°C下的动态剪切模量主曲线数据与分子模拟结果进行比较，以验证其可靠性。

此外，还将进行动态剪切流变仪（DSR）测试（型号：Superpave 102e），包括线性幅度扫描（LAS）和时间扫描（TS）测试，以表征实际沥青粘合剂在不同应力幅度下的线性和非线性粘弹性行为以及损伤演变。这些测试将用于验证宏观尺度实验和分子模拟在粘弹性和损伤阶段之间的一致性。LAS和TS测试是在江苏新材料有限公司生产的沥青粘合剂上进行的。根据中国技术标准JTG 3410-2025，沥青粘合剂的基本技术指标见表4。LAS测试按照AASHTO TP 101标准进行。所有测试采用一致的参数：10 Hz的正弦剪切振荡加载、25°C的测试温度以及8 mm直径的平行板和2 mm的间隙。详细分析将在后续部分展示。

表4. 70#基础沥青的技术指标。

测试项目 单位 技术规格 测试结果
penetration (25°C, 100g, 5s) 0.1mm 60-80 71
softening point (°C) ≥46 48
ductility (5cm/min, 10°C) cm ≥15 55
rotational viscosity (60°C) Pa·s ≥180 266
RTFO % ≤0.8 0.03
residual ductility (5cm/min, 10°C) cm ≥67 66

3. 结果与讨论
3.1. 基于第2节中平衡的沥青粘合剂模型的循环剪切加载下的动态应力-应变响应
使用AAA-1×10模型作为代表，分析其在动态剪切加载下的应力-应变响应特性。动态剪切加载模拟在25°C下进行，输入应力幅度为0.2 GPa，加载周期为12 ps。图4(a)连续收集了10个加载周期的输入应力、响应应力和应变的波形数据，图4(b)和4(c)显示了第6-10个周期的应力-应变滞后环。

如图4(a)所示，应变曲线相对于输入应力表现出典型的相位滞后，反映了沥青粘合剂的粘弹性。响应应力幅度低于输入值，主要是由于系统的响应时间有限。每个应力点下的短动态模拟持续时间阻止了系统与施加的应力达到完全的瞬时平衡。图4(b)和4(c)中观察到明显的椭圆形滞后环，这与实际沥青粘合剂的动态行为一致。总体而言，所提出的动态建模算法成功再现了关键的纳米尺度粘弹性特征，具有稳定且连续的应力-应变响应，可以支持进一步的机械分析。

3.2. 沥青粘合剂动态剪切模量的分子尺度上的时间-温度叠加
动态剪切模量是表征沥青粘合剂粘弹性行为的关键流变参数。然而，分子动力学模拟的固有时间尺度限制使得可访问的加载频率远高于宏观动态剪切流变仪实验中使用的频率，直接进行定量比较具有挑战性。根据时间-温度叠加原理，可以在高温和短时间或低温和长时间下观察到等效的机械松弛行为。因此，通过提高模拟温度来加速分子松弛，可以在分子尺度上获得的高频加载下的动态剪切模量通过构建主曲线转换为等效的低频行为。这种方法旨在扩展纳米尺度动态剪切模量分析的时间尺度，缩小模拟和实验加载频率之间的差距。

按照第2节中描述的模拟方法，在XY平面上对AAA-1×10、AAK-1×10、AAM-1×10和四组分沥青粘合剂模型施加了应力控制的动态剪切加载，加载条件总结在表3（Set I）中。对于每个加载条件，动态剪切模量在10个加载周期内取平均值以获得最终值。在不同温度下获得的动态剪切模量沿频率轴水平移动，并使用方程(31)中的S形函数模型进行拟合，通过最小化平方差之和来优化拟合。所得参数用于构建参考温度下沥青粘合剂分子模型的主曲线。温度移动因子使用温度依赖的移动因子方程(32)、(33)确定。这些模型在25°C参考温度下的动态剪切模量和相应的主曲线显示在图5(a)–(g)中。表5总结了温度移动因子、上述方程的拟合参数以及相应的皮尔逊相关系数。主曲线和温度移动因子方程的皮尔逊相关系数均高于0.99，确认时间-温度叠加原理在分子尺度上有效适用。

(31)
(32)
(33)

其中 表示降低后的频率， 表示相应的动态剪切模量； 是加载频率（Hz）， 表示模拟温度，而 是温度 下的移动因子。参数 、 、 和 是动态剪切模量主曲线方程的拟合参数，而 、 和 是移动因子方程的拟合参数。

下载：下载高分辨率图像（2MB）
下载：下载全尺寸图像

图5. 分子模型的动态剪切模量、主曲线和统计分析。

表5. 不同温度下的时间-温度移动因子、方程拟合参数和皮尔逊相关系数。

模型 不同温度下的时间-温度移动因子 移动因子方程参数 相关系数
AAA-1×10 10.76 0.32 0.60 47 0.45 0.33 67 0.40 50 33 67 0.45 0.34 31 > 0.99
AAK-1×10 10.76 0.76 0.60 47 0.59 21 40.41 92 30 34 51 84 24 0.34 31 > 0.99
AAM-1×10 10.73 23 50.58 86 0.43 77 60 0.31 75 24 0.34 31 > 0.99
As1 10.86 0.86 27 10.55 76 80 0.38 48 83 98 65 0.05 41 > 0.99
Re 10.83 82 10.65 32 20 51 79 0.40 68 74 0.38 48 83 98 65 0.05 41 > 0.99

表5. 主曲线模型拟合参数和相关系数

考虑MD模拟中的固有热力学波动，对模拟的动态模量进行了统计分析，以确保结果的可靠性。图5(h)展示了在本研究定义的边界条件下，AAA-1、AAK-1和AAM-1×10沥青模型经过十个加载周期后获得的模拟动态剪切模量。结果涵盖了最低和最高的模拟温度（298.15 K和458.15 K）以及最小和最大的加载频率，重点关注动态响应最敏感且容易波动的边界条件。此外，图5(k)比较了两次独立重复模拟得到的循环平均动态剪切模量。

在图5(h)中，所有变异系数均低于5%，表明模拟保持了较强的循环间稳定性。图5(k)进一步显示，两次重复模拟之间的平均动态剪切模量的决定系数（R2）超过0.99，确认了独立运行之间的一致性。总体而言，模拟框架可以提供高统计可靠性和可重复的输出，为后续的定量分析提供了坚实的基础。在图5中，所有分子模型都表现出温度和频率依赖的粘弹性。动态剪切模量随温度升高而降低，随加载频率增加而增加，这与宏观行为一致。

在图5(m)中，沥青质显示出最高的动态剪切模量，其次是AAK-1×10、AAM-1×10和AAA-1×10。树脂和芳香族的曲线几乎相同，而饱和烃的模量最低。沥青粘合剂模型的相对排名与SHRP-A-645报告中的宏观实验结果一致[44]。例如，在60°C和10 rad/s的加载频率下，实验动态剪切模量从AAK-1（5.4 kPa）降至AAM-1（3.8 kPa）和AAA-1（2.8 kPa）。沥青质的高动态剪切模量反映了它们较大的分子量和刚性的芳香族结构，而饱和烃主要由灵活的长链分子组成，表现出最弱的剪切变形抵抗力。此外，树脂和芳香族的几乎重叠的主曲线表明它们的流变响应相似，支持将这些组分等同处理的建模方法[45]。图5(n)展示了从两次模拟和SHRP-A-369报告[46]获得的AAA-1和AAM-1沥青在25°C下的动态剪切模量主曲线。由于缺乏AAA-1的实验主曲线数据，比较仅限于这两种粘合剂。它关注模拟和实验数据重叠的高频区域。尽管时间和空间尺度相差九个数量级，模拟的动态剪切模量再现了实验观察到的趋势，模拟和实验结果之间的差异小于半个数量级。

由于沥青粘合剂的组成复杂性、时间和空间尺度的限制、力场的代表性以及边界条件的理想化，分子模拟和宏观实验之间的绝对数值差异是不可避免的。尽管如此，本研究的主要目的不是实现精确的一一对应，而是建立一个分子尺度框架，捕捉关键变量对动态机械响应的物理合理影响，并验证模拟方法的内部一致性。结果表明，尽管存在数值差异，模拟的主曲线成功再现了实验观察到的温度-频率依赖性、相对大小和粘合剂模量的相对排序，模拟值与实验数据处于同一数量级。因此，当前研究为未来将粗粒化和离散元模拟技术结合起来实现动态模量主曲线的跨尺度集成提供了方法论基础。

3.3. 沥青粘合剂的分子尺度粘弹性和损伤行为
粘弹性材料在连续加载下表现出不同的机械响应机制，其主导模式由施加的应力水平决定。对于沥青粘合剂，宏观行为通常分为三种机制：线性粘弹性、非线性粘弹性和疲劳损伤，基于动态模量的应力依赖性和时间演变。在线性粘弹性机制下，动态模量与应力幅度无关，并在加载过程中保持恒定。在非线性粘弹性范围内，模量随着应力幅度的增加而减小，但随时间仍然没有衰减，表明材料处于未受损状态。当动态模量随着应力幅度的增加和加载时间的延长而下降时，材料已经进入了疲劳损伤阶段，其特征是结构发生不可逆的退化。图6展示了实际沥青粘合剂的线性幅度扫描和时间扫描测试结果。如图6(a)所示，在10 kPa到90 kPa的应力幅度范围内，动态模量基本保持恒定。根据图6(b)中的时间扫描测试结果，在30 kPa时，动态剪切模量随着应力幅度的增加而减小，但随着时间的推移基本保持不变，这证实了非线性粘弹性的开始。一旦应力幅度达到50 kPa，模量开始随着加载时间的延长而逐渐减小，表明材料已经过渡到疲劳损伤阶段。基于这些实验确定的粘弹性和损伤特性，后续部分将使用开发的动态加载模拟方法来验证其在分子尺度上捕捉沥青粘合剂动态力学响应的应力依赖性和时间演变的能力。

为了实现长时间循环加载模拟，采用了AAM-1×1模型以确保在现有计算资源下的计算效率。对AAM-1×1模型的Z轴施加正弦循环拉伸-压缩应力，周期为4.8 ps，模拟温度为25°C，应力幅度范围从0.1到1 GPa。对于每个应力幅度，进行了500次连续的加载循环以得出动态模量，并进行了两次独立的模拟以确保统计稳健性。如图7(a)–7(j)所示，两次模拟的原始动态模量数据在加载循环中紧密重叠，LOWESS平滑后的演变曲线显示出几乎相同的趋势。在图7(m)中，两次独立模拟的500次加载循环的平均动态模量的决定系数（R2）达到0.9984。结合第3.2节的统计分析，所提出的模拟框架在分子组成、加载频率、温度和应力幅度的变化下保持了出色的统计稳定性和可重复性。

在图7(n)中，随着输入应力幅度从0.1 GPa增加到0.6 GPa，动态模量随着加载水平的增加而单调减小，而在加载循环中基本保持不变。这反映了沥青粘合剂的典型非线性粘弹性响应。当输入应力幅度达到0.7 GPa时，动态模量开始随着加载循环的增加而逐渐减小，表明损伤演变的开始。增加的应力幅度加剧了微观结构的改变，从而促进了微裂纹的起始和传播。因此，0.7 GPa可以被视为AAM-1×1模型中粘弹性阶段与疲劳损伤阶段的分界阈值。

如图7(k)和7(l)所示，在低输入应力幅度（0.1 GPa–0.4 GPa）下，动态模量随着加载循环的出现而明显波动。这些波动主要源于两个因素。首先，由于分子模型的尺寸相对较小，固有的各向异性和热力学波动更加明显，这自然增加了应力-应变响应的变异性。其次，当施加的应力幅度较小时，施加的应力水平更接近模型的固有热力学应力波动，使得测量响应对这些内部扰动更加敏感。统计和计算误差可能掩盖了模型的线性粘弹性行为。基于此，将在以下小节中系统地研究模型尺寸尺度和应力幅度对动态力学响应稳定性的影响。

本节将系统地分析模型尺寸尺度和应力幅度对动态模量结果稳定性的耦合效应。根据先前的研究[31]，增加模型尺寸尺度可以显著减少固有的各向异性。因此，在本节中，使用了表2中的五个放大尺度模型（AAM-1×2至AAM-1×10）在0.1 GPa到1 GPa的输入应力幅度范围内进行动态加载模拟。模拟遵循上一节中使用的确切条件，在可用计算资源下进行了50次加载循环。如图8(a)–8(e)所示，随着模型尺寸尺度和输入应力幅度的增加，沥青模型的应变曲线逐渐从不稳定和不规则的图案演变为定义明确、周期性且幅度稳定的响应。模型尺寸尺度的增大显著减少了模拟输出中的统计噪声。增加施加的应力幅度放大了模型的真实变形信号，增加了有效应力/应变在总响应中的比例，进一步减少了噪声干扰。在图8(f)和8(g)中，1 GPa输入应力幅度下的动态模量波动比0.1 GPa时小十倍。在相同的0.1 GPa输入应力幅度下，AAM-1×10模型的动态模量波动大约是AAM-1的五分之一。

图8(m)总结了在0.1 GPa到1 GPa应力幅度范围内，不同尺寸尺度模型的50次加载循环的平均动态模量。随着应力幅度的增加超过0.4 GPa，模型之间的平均动态模量差异减小。在0.6 GPa以上，模量几乎不受模型尺寸尺度的影响。在0.1–0.3 GPa的低应力范围内，模型尺度对模拟的动态模量有明显的影响。随着模型尺度的增加，平均动态模量逐渐趋于一致。对于AAM-1×10模型，在0.1 GPa到0.3 GPa范围内，平均动态模量不再随应力幅度显著变化。图8(n)展示了在不同应力幅度下AAM-1×10模型50次加载循环的相应LOWESS平滑动态模量。在0.1 GPa到0.3 GPa范围内，动态模量基本保持不变。因此，0.4 GPa被确定为AAM-1×1模型从线性粘弹性行为转变为非线性粘弹性行为的应力幅度阈值。图8(m)的结果表明，当输入应力幅度达到0.6 GPa时，动态模量基本上不受模型尺度的影响。因此，0.7 GPa仍然是从非线性粘弹性到损伤转变的阈值，这与第3.3节的研究结果一致。

总体而言，观察到的动态模量的尺度依赖性主要源于由固有数值波动引起的统计偏差，通过增加模型尺度或使用更高的应力幅度可以有效减少这些偏差。未来的改进可以考虑优化模拟框架，例如使用粗粒度模型来简化原子间相互作用，或结合机器学习力场来扩展可访问的时间和空间尺度。这样的发展将允许在更大的时间和空间尺度上进行模拟，为沥青粘合剂的粘弹性和损伤演变的定量分析提供坚实的数值基础。

疲劳损伤、裂纹和其他失效行为本质上是能量耗散过程。由于使用大尺度AAM-1×1模型可靠地识别了线性粘弹性阶段，因此当前分析重点关注非线性粘弹性和损伤阶段，在这些阶段中，能量耗散随着应力水平和加载循环的进行而显著变化。选择了AAM-1×1在不同输入应力幅度（0.4 GPa、0.6 GPa、0.8 GPa和1 GPa）下的应力-应变数据来计算不同机械响应阶段的耗散能量。使用梯形法则通过方程(34)、(35)计算了代表能量耗散的滞后环面积，并显示在图9(a)中。计算出的能量耗散结果显示在图9(b)–9(e)中。图9(f)展示了不同输入应力幅度下的能量耗散的LOWESS平滑比较。

如图9所示，在非线性粘弹性阶段，耗散能量并没有随着加载循环次数的增加而显著增加，而是水平波动。在未受损的非线性粘弹性阶段，沥青粘合剂模型内的分子链在循环加载下经历了可逆的重排，允许变形发生而不引起内部损伤。能量耗散主要来自分子间摩擦和局部结构调整，某些区域表现出组分流动。随着输入应力幅度的增加，能量耗散逐渐增加。这种增加主要是由于在更高应力水平下分子链的运动更加剧烈。升高的应力幅度导致更强的内部摩擦和更频繁的分子链及局部结构的重排，即使没有发生内部损伤。这些分子相互作用增加了每个加载循环中损失的能量。

在0.8 GPa和1 GPa输入应力幅度的损伤阶段，耗散能量随着加载循环的增加而逐渐增加。此外，每个循环的耗散能量显著高于粘弹性阶段。这种增加可以归因于在高应力加载下沥青粘合剂模型内部发生的不可逆损伤。损伤和裂纹过程需要更多的能量来打破加载过程中的局部分子相互作用。此外，随着损伤的加深，模型内部裂纹腔的发展速度加快，导致每个循环的耗散能量进一步增加。

在500次加载循环结束时，获得了AAM-1×1模型在0.4 GPa、0.6 GPa、0.8 GPa和1 GPa输入应力幅度下的结构快照，以评估输入应力幅度对模型结构配置的影响。图10显示了AAM-1×1模型的三维配置和加载结束时的自由体积。图10中的粉红色区域代表Connolly表面的等值面（半径1 ?），这些区域表示具有1 ?半径的探针球体的可访问空间边界。这些区域代表了材料的固有分子空隙，由于力场控制的分子间相互作用，这些空隙是自然分散且不连续的。在0.4 GPa和0.6 GPa的循环加载下，分子重排导致更多的连通空隙路径，其大小相对稳定。因此，结构达到了动态平衡，机械性能不会随着进一步循环而退化。在高应力幅度（0.8–1.0 GPa）下，循环加载导致相互连接的空隙网络和更连续的等值面分布，反映了分子尺度损伤的起始和传播。同时，更高的应力幅度导致加载结束时自由体积的增加，这是由于更明显的分子重排和局部不稳定性。在宏观尺度上，这对应于分子接触网络的减弱和结构紧凑性的丧失，导致动态模量的降低。结果支持使用分数自由体积作为沥青模型中变形抵抗力的指标。

本研究开发了一种分子尺度的动态建模算法，用于研究沥青粘合剂的动态力学行为。使用典型的沥青模型（AAA-1、AAK-1和AAM-1），该算法分析了动态模量对温度、频率、应力幅度和加载时间的依赖性。具体结论如下：(1)所提出的动态建模算法成功再现了沥青粘合剂的典型粘弹性行为，包括应变相位滞后和滞后环特性。此外，所得到的应力和应变响应数据的连续性和稳定性得到了很好的保持。(2)时间-温度叠加原理在分子尺度上也是适用的。分子动力学（MD）模拟能够准确捕捉动态剪切模量的温度和频率依赖性趋势，并再现了不同模型之间的相对排序，这与实验观察结果一致。由于沥青粘合剂的成分复杂性以及时间和空间尺度上的九个数量级差异，与实验结果的绝对数值差异是不可避免的。尽管如此，模拟得到的动态剪切模量与实验结果的偏差小于0.5个数量级。(3)在恒定温度和加载频率下，沥青分子模型的动态模量在加载循环过程中表现出三个特定阶段的响应：线性粘弹性、非线性粘弹性和损伤演化。增加输入应力幅度和模型尺度可以提高动态模量模拟的稳定性。对于AAM-1沥青粘合剂模型，0.4 GPa和0.7 GPa的应力幅度是区分这三个阶段的阈值。(4)在粘弹性阶段，每个循环中耗散的能量几乎保持不变。在损伤阶段，随着加载循环的进行，耗散能量增加，同时伴随着孔隙的扩展和粘性流动的增强。能量耗散主要克服了分子间的作用力，推动了裂纹的扩展和孔隙的形成，这与实际沥青粘合剂中观察到的损伤演化是一致的。(5)较高的应力幅度会加速结构演化，增加自由体积的比例并降低动态模量。这反映了结构紧凑性的丧失和变形抵抗力的下降。在微观尺度上，内部孔隙连通性的增加表明在循环加载下逐渐发生疲劳损伤。4.2. 展望虽然本研究已经达到了既定的研究目标，但在分子到宏观尺度上提高动态力学性能的预测准确性仍有改进的空间。未来工作的潜在方向如下。(1)开发一种用于沥青的粗粒度模型。简化分子结构可以降低计算成本，并允许在更大的空间和时间尺度上进行模拟，同时保持与宏观材料的结构相关性。(2)开发基于机器学习的力场。将机器学习纳入力场参数化可以提高计算效率和准确性，从而实现大规模、长时间的模拟。(3)整合多尺度模拟技术。将分子尺度上的时间-温度等效性与介观方法（例如离散元模型）相结合，可以将动态模量的预测扩展到更低的频率，并建立一个涵盖分子到宏观尺度的框架。CRediT作者贡献声明：田伟伟：撰写——审阅与编辑、撰写——原始草稿、监督、软件开发、资金获取、正式分析、数据管理、概念化；李辉：撰写——审阅与编辑、监督、方法论、资金获取、正式分析；张玉清：撰写——审阅与编辑、验证、项目管理、资金获取、正式分析。数据可用性：数据将根据请求提供。访问COMPASS力场文件的链接为：https://doi.org/10.6084/m9.figshare.30860330。