邱延峻， 苏 婷， 郑鹏飞， 丁海波,3

(1.西南交通大学 土木工程学院， 四川 成都 610031； 2.西南交通大学 道路工程四川省重点实验室， 四川 成都 610031； 3.加拿大女王大学 化学学院， 安大略 金斯顿 K7L)

沥青老化分为化学老化和物理老化.沥青化学老化是指沥青在贮存、运输、施工及使用过程中由于长时间暴露在空气中，在环境因素的作用下发生挥发、氧化、聚合，从而导致沥青内部分子结构发生变化，同时发生化学性质变化，导致其路用性能劣化的过程[1].玻璃态聚合物在低于玻璃化转变温度(Tg)的存放或使用过程中，其力学性能会随时间发生变化，但在此过程中，材料的化学结构未发生改变，且这种力学性能的改变可以通过对材料进行热处理(升温至玻璃化转变温度以上)来消除，这一改变称为物理老化[2].沥青在低于其Tg的温度下也会出现与玻璃态聚合物一致的物理老化现象.近年来，随着各种沥青技术的发展，沥青胶结料成分发生了明显的变化.首先，为了节约能源及降低拌和温度，蜡类温拌剂的使用越来越频繁；其次，为了提高废弃物资源化利用率，回收沥青路面材料(RAP)、废机油底渣及煤液化残渣等废弃物在沥青中的添加比例逐渐提高；再次，为了降低沥青成本，多聚磷酸(PPA)等低成本酸类添加剂得到更多商家的青睐；最后，全球对石油副产品需求量的增加使得石油炼制厂改进精炼工艺，从原油中蒸馏出更多的轻组分.前期研究表明[3-4]，这些因素均可能增加沥青物理老化的程度，从而使低温环境下沥青路面开裂的风险增加.因此，系统地开展沥青物理老化机理的研究，可用于指导沥青生产工艺的改进、控制各种废弃物的添加比例，以及选择优质的改性技术，对提高沥青混凝土结构的耐久性和改善其在低温下的路用性能具有重要意义.

物理老化过程普遍存在于各种玻璃材料和非平衡系统中，系统的特性在物理老化(或结构恢复)的过程中随着时间的推移逐渐向均衡方向发展[5].物理老化是一种可逆老化，人们通常把物理老化归因于自由体积或分子链构象[1].王秀娟[6]利用分子模拟分析了二烯类橡胶的老化在分子层面的表象.燕楠[7]研究了天然橡胶在热和压力作用下的物理老化行为.王兵等[8]利用分子模拟预测了静压力条件下黏弹性阻力材料的自由体积参数.而目前尚无专门针对沥青物理老化行为的模拟与分析.源自物理学的分子动力学(MD)在特定的时间段内可以追踪分子运动并预测材料的性质，被广泛应用于不同的研究领域来模拟原子尺度上的材料及其行为[9-10].目前，分子模拟在沥青材料中的研究主要集中于研究组分间相互作用、各组分聚集状态等方面.王鹏等[11]运用分子模型探究沥青蜂状结构成因.朱建勇[12]模拟了沥青的自愈合行为.王岚等[13]做了沥青老化前后与胶粉相容性的研究.目前尚无关于沥青模型物理老化过程的分子动力学模拟.沥青的玻璃化转变温度是沥青是否发生物理老化的分界点，曹丽萍等[14]提出可以用玻璃化转变温度来评价SBS改性沥青的低温性能，刘佳音[15]基于黏弹性模型提出了沥青玻璃化转变温度的确定办法，但目前尚无文献从微观层面分析和确定沥青的玻璃化转变温度.

微观结构变化到一定程度后会反映出宏观性能的变化，自由体积与聚合物的黏度、黏弹性、玻璃化转变、结构弛豫以及塑性屈服等物理化学性质密切相关.均方位移随外界条件的变化反映了聚合物分子活动随外界条件的变化，因此对均方位移的研究能够从微观上分析聚合物宏观性能的变化.对微观结构在自由体积、回转半径等方面的研究，能够帮助人们了解宏观物体的老化机理.本文通过改变温度使系统处于“非平衡”状态，对沥青的物理老化进行模拟，得到沥青的玻璃化转变温度，结合试验，观察特定温度下沥青的密度、自由体积随时间的变化情况；然后对沥青各组分进行相同条件下的分子动力学模拟，分析沥青出现物理老化的内在原因.

1 沥青四组分模型

1.1 力场与系综

力场是一种为分子动力学模拟中的每种原子提供参数，以计算系统能量的能量函数.Compass力场已被证实可以精确预测孤立体系或凝聚态体系中各种分子的结构、构象、振动以及热物理性质[16].系综是指通过算法对模拟环境设置不同的宏观条件以达到反映体系内分子运动的物理过程[6].Accelrys公司开发的Materials Studio(MS)软件中有正则系综(NVT)、微正则系综(NVE)、恒温恒压系综(NPT)和等焓等压系综(NPH)4种热力学系综可供选择，以确定不同的宏观特性.本文模拟使用以下2种热力学系综：NVT，用以保持系统体积(V)、温度(T)和系统中原子数(N)恒定；NPT，用以保持系统压力(P)、温度(T)和系统中原子数(N)恒定.

1.2 建立模型

运用单一的沥青分子平均模型难以描述沥青的特性，在此选用沥青的四组分(芳香分、饱和分、胶质、沥青质)模型[17]相对精细化地对沥青分子结构进行描述.分子模型的构建主要参考朱建勇[18]的研究工作，用四组分平均分子模型代表沥青质分子构型[19]，渺位缩合的分子模型代表胶质分子构型[20]，1,7-二甲基萘代表芳香分分子构型[21]，正二十二烷代表饱和分分子构型[21]，从而建立基质沥青的分子模型.

本文所进行的分子动力学模拟基于MS软件，首先利用Visualizer模块画出各个组分的结构单元，优化后的各组分分子构型如图1所示.

图1 优化后的各组分分子构型Fig.1 Optimized molecular structure of components

然后，利用Amorphous Cell模块构建无定形的沥青分子模型，设定相对密度1.0g/cm3，将其放入周期性单元格中，周期性边界条件用于控制最小化边缘效应和模拟宏观特性.沥青分子模型的基本数据如表1所示，根据某基质沥青中的四组分含量确定.

表1 沥青分子模型的基本数据

沥青分子模型较为复杂，在对其进行分子动力学模拟之前，需要利用能量最小化来消除系统中原子间的不良接触.对沥青分子模型进行结构优化、退火处理，得到如图2所示优化后的沥青分子模型.

图2 优化后的沥青分子模型Fig.2 Optimized asphalt molecular model

1.3 模型验证

为保证模拟结构具有代表性，需要对所建沥青分子模型的准确性进行验证.选取沥青在298K下的密度与实际沥青材料的性能指标进行对比.此时沥青模型的温度波动在控制范围内已达到相对平衡的状态.298K下实际沥青材料的密度在1.00～1.04g/cm3之间，所建沥青分子模型密度的算术平均值为1.018g/cm3，其差异控制在5%以内.

利用Forcite模块的Analysis对沥青分子模型分子内和分子间的原子径向分布函数g(r)进行分析(r为距中心分子的距离)，结果见图3.由图3可以看出：沥青分子模型分子内的原子径向分布函数峰后曲线平缓且趋近于0；分子间的原子径向分布函数曲线逐渐上升，达到最大值后，曲线平缓且趋近于1.沥青是典型的非晶体材料，原子在大分子内呈有序状态排列，近程有序，远程无序，且分子间作用力以范德华力为主[21].根据图3可以认为所建沥青分子模型能代表真实沥青分子结构.

图3 沥青分子模型分子内和分子间原子径向分布函数Fig.3 Radial distribution function(RDF) for distances of intra-molecular and inter-molecular atoms in asphalt molecular models

2 沥青模型的分子动力学模拟

2.1 玻璃化转变温度

选取195、215、235、255、273、303、333、363、393K这9个温度点，计算不同温度下沥青的密度和体积.在NPT系综下进行分子动力学模拟，压力为0.1MPa，当体系温度波动控制在10%以内时，认为分子动力学模拟到达平衡，提取沥青分子模型达到相对平衡状态后温度最接近设定温度298K的1帧作为分子动力学模拟的起点，来计算密度和体积的算术平均值[22].不同模拟温度下沥青分子模型的密度如表2所示.由表2可知，随着模拟温度的升高，沥青模型的密度不断减小.

表2 不同模拟温度下沥青分子模型的密度

沥青模型体积-温度(V-T)曲线如图4所示.由图4可知，沥青模型的体积在195～255K和303～363K范围内基本呈线性增长,但高温阶段的斜率明显大于低温阶段.在255～273K范围内，高温阶段斜率发生明显变化，体积出现波动.沥青体系在玻璃态和黏弹态下具有不同的属性，在玻璃化转变温度附近，密度、体积等属性会出现较大的变化[15].因此，可以推断沥青分子模型的玻璃化转变温度在255～273K范围内.

图4 沥青模型体积-温度曲线Fig.4 Volume-temperature relationship diagram of asphalt model

由此，可以通过在曲线拐点两侧分别拟合线性回归曲线的方法来预测体系的玻璃化转变温度[23]，得到分段线性拟合沥青模型体积-温度曲线.根据线性拟合可以得出2条拟合线的交点对应温度为263K，即该沥青模型的玻璃化转变温度为263K.

从分子动力学模拟结果中导出不同温度(195、215、235、255、273、298、303、333、363K)下沥青模型的均方位移(MSD)曲线，如图5所示.均方位移是作为时间函数来计算并表征平均原子位移的，即运动粒子位移平方的平均值[14]，可以采用均方位移来分析沥青模型在不同温度下分子动力学模拟的运动轨迹的区别.

图5 不同温度下沥青模型的MSD曲线Fig.5 MSD curves of asphalt model at different temperatures

从图5可以看出：温度的升高使分子链的活动性增强，均方位移随温度升高而不断增加；在195～255K区间内，均方位移缓慢上升；在255～363K区间内，均方位移的变化较快.可以推测在255～273K区间存在1个过渡点，使得均方位移改变的速率出现较大变化.这一情况对应了玻璃化转变温度所在区间.

为了验证利用分子动力学模拟计算的基质沥青玻璃化转变温度的可靠性，对四组分比例与沥青模型相近的某基质沥青采用美国TA仪器公司的DSC 2500差示扫描量热分析仪进行了玻璃化转变温度的测试.消除热历史后，快速降温至-88℃，调制后温度变化速率为0.5℃/min，再从-88℃升高温度至140℃，得到基质沥青的可逆热容与其温度的1阶导数曲线，如图6所示.图6中曲线最高点对应的温度即试验所得基质沥青的玻璃化转变温度，其值为259K，与模拟所得玻璃化转变温度(263K)非常接近.但模拟所得玻璃化转变温度高于试验值，这是由于分子动力学模拟过程中的降温速率较快所致.

图6 基质沥青的可逆热容与其温度的1阶导数Fig.6 Reversing heat capacity and its temperature derivative of asphalt binder

2.2 沥青的物理老化模拟

在玻璃化转变温度附近，会发生低温条件下沥青模量变化的物理老化现象[6].根据试验和玻璃化转变温度预测，本文将沥青物理老化模拟的温度设定为255K.先在目标温度下对起点模型进行100ps 的NVT系综下的分子动力学模拟，沥青分子模型完全弛豫.然后再对该过程的最后1帧模型进行800ps的NPT系综下的分子动力学模拟，得到体系在该温度下性质随时间变化的情况.

根据密度随时间的变化关系，选定密度变化的特别点，利用软件中Atom Volumes & Surfaces选项对体系的自由体积进行测算，得到如图7所示沥青模型的密度与自由体积随时间变化曲线.由图7可以看出：沥青的密度随时间的增加而逐渐增加，但密度的改变速率不断降低，当增加到一定程度时沥青密度保持相对稳定；自由体积随时间的增加而降低，然后在一定区间内波动变化，波动维持在1个相对稳定的区间；密度的升高对应着自由体积的降低，密度与自由体积的变化呈现反比关系.

图7 沥青模型的密度与自由体积随时间变化曲线Fig.7 Density and free volume variation curves with time of asphalt model

自由体积的存在是非平衡非晶态的直接结果，对聚合物的密度和力学性能有着重要的影响.当系统处于非平衡状态时，系统的特性会在物理老化或结构恢复过程中随着时间的推移逐渐向均衡方向发展[14].由自由体积理论可知，自由体积的下降是通过聚合物空穴的重新分布产生的，其下降速率取决于分子链段的运动能力.体积松弛过程是非线性的，自由体积的减小使分子链段的运动能力下降.而分子链段的运动能力和自由体积的松弛速率都随着老化时间的增加而降低，自由体积的松弛速率只能无限接近但不等于零.这与密度的变化具有一致性.

导出沥青分子模型在255K下0～800ps时间段内分子动力学模拟的连续动画，图8所示为其第1帧和最后1帧图像.由图8可见：从模拟开始到模拟结束，沥青分子模型的体积减小，同时，四组分的表现各不相同；沥青质和胶质为核心，饱和分和芳香分为外侧结构，饱和分和芳香分的分布相对呈游离态；在分子动力学模拟过程中，沥青质和胶质逐渐聚集，但活动缓慢，芳香分和饱和分运动相对活跃.

图8 沥青模型动力学模拟输出图像Fig.8 Molecular simulation output image of asphalt model

2.3 沥青各组分变化情况

由分子动力学模拟输出的图7、8可以看出，沥青的化学组成与其物理老化过程中性能变化有着密切的联系.利用Amorphous Cell模块建立沥青四组分各自的分子模型，研究微观尺度下各组分的物理老化情况，推断其在沥青物理老化过程中所扮演的角色.各模型中的相对分子质量与沥青模型中该组分的相对分子质量一致，建模相对密度均设定为1.0g/cm3，其余分子动力学模拟过程与沥青模型的分子动力学模拟过程一致.得到沥青模型与各组分模型分子动力学模拟结果，如表3所示.由表3可见：从298K到255K,各模型密度均增大.其中，饱和分模型与芳香分模型在298、255K时的密度均小于建模密度，这与饱和分和芳香分的结构有关.

表3 沥青模型与各组分模型分子动力学模拟结果

物理老化是通过链段的微布朗运动实现的[24]，沥青和沥青各组分模型在255K、0.1MPa、NPT系综下的均方位移(MSD)曲线如图9所示.由图9可以看出：沥青和沥青各组分模型的均方位移均随温度升高而增大；芳香分模型的均方位移曲线明显大于其他组分，其次是饱和分模型，这2种模型的均方位移均大于沥青模型；沥青质模型的均方位移曲线小于胶质模型，两者的均方位移均小于沥青模型.

图9 沥青与沥青各组分模型的MSD曲线Fig.9 MSD curves of asphalt model and component models

温度增加，沥青模型的体积和自由体积均发生膨胀，分子间距离增大，分子链段的活动性增强[8].芳香分的结构简单，相对分子质量小，同时其拥有最少的键组成，不具有长链结构，键和角很少，能量变化也最小，因此芳香分的活动能力更大.芳香分具有非极性的芳香环结构，为沥青质和胶质组分提供了可塑性.饱和分是非极性链烃分子，具有较高的氢碳比，也具有较强的活动能力，为整个沥青分子提供流动性.这也可以用来解释芳香分和饱和分的密度小于建模密度的原因.而沥青质分子具有较大的分子结构，与其他组分相比具有更高的键拉伸能，为沥青提供了黏性和胶体特性，故其分子活动能力最弱.沥青质分子被胶质分子包围，从而赋予沥青延展性.可以认为，沥青质和胶质抑制了沥青的活动性，从而直接导致物理老化，当沥青模型在255K下进行链段的微布朗运动时，主要的分子活动力来自芳香分和饱和分.

3 结论

(1)利用分子动力学模拟了沥青分子模型密度随温度变化的规律，并通过体积-温度曲线得到了沥青模型的玻璃化转变温度，通过试验验证了所得玻璃化转变温度的可靠性.

(2)沥青在255K下分子动力学模拟过程中，密度增大，自由体积减小，出现了物理老化的情况.

(3)低温条件下沥青分子中的沥青质与胶质组分逐渐聚集、分子活动性降低，从而导致了物理老化现象，该过程中微布朗运动的分子活动力主要由芳香分和饱和分承担.沥青分子的物理老化行为与各组分的性质有关，更与材料间形成的聚集态有关.本文仅从各组分本身出发，简要分析了各组分与沥青物理老化的关系.

(4)分子动力学模拟发生在纳米时间尺度上，而物理老化实际发生在宏观时间尺度上，表现为沥青低温性能改变，此处用分子动力学模型模拟物理老化过程无法完全得到由时间因素决定的沥青材料的物理性质；四组分沥青模型和有限数量的分子也不能完全代表复杂的沥青材料.但是分子动力学模拟为沥青的性质预测提供了新的解决方法，后期可以通过不断解决这些问题来进行更接近的预测，进一步改善模型.