通過分子動力學方法研究瀝青質聚集體應變下勢能變化

葉正揚， 樊小哲， 龍 軍， 代振宇

(中國石化 石油化工科學研究院，北京 100083)

石油瀝青是一種重要的道路交通材料，其性能直接決定了瀝青路面的工程質量和結構形式。和水泥混凝土剛性路面不同，瀝青路面屬于柔性路面，相比具有行車舒適平穩、燃油消耗和機械磨損低的優點[1]。因此，改進瀝青材料性能，研究瀝青結構與性能之間的關系兼有經濟價值和戰略意義。目前，改進瀝青材料性能的主要方法是添加其他物質。最常見的添加物是SBS膠粉等高分子物質[2]，或添加石油樹脂等其他石油組分進行調和實現改性[3]，另外也有通過添加碳納米管等新型碳基材料進行瀝青改性的研究[4]。盡管這些改性方法都不同程度地實現了瀝青改性需求，但沒有研究瀝青材料自身組成結構與性能關系。

瀝青各相展現的形變狀態、剛度和設計允許形變范圍內的強度均與分子形變和分子間相互作用有關：對瀝青材料施加外力時，瀝青中分子間距和分子構型發生改變使材料整體應變，從而通過應力對抗外部作用。研究不同結構的分子在不同環境應變下的勢能變化情況，可以體現材料整體形變下分子結構和相互作用的變化，由此揭示不同結構分子對材料內應力的貢獻，從而為提升瀝青材料使用性能提供理論依據和方向。

瀝青的分子組成包含數萬種分子，復雜的分子組成導致其內部結構的復雜性。通常對瀝青化學成分的分析僅停留在擬組分分析，例如Corbett分析法[5]，這些方法難以關聯瀝青組成參數和材料性能之間的關系。通過原子力顯微鏡觀察發現，瀝青材料并非均相體系，而是在微米尺度由力學性質完全不同的多個相組成，這些由不同分子組成的相之間相互交錯滲透[6]。以瀝青質為例，在若干納米的尺度上具有自發相互聚集的特點，無論在減壓渣油[7]或小分子溶劑[8]中均可形成具有分形結構的納米聚集體，這些聚集體相互連接呈100 nm左右的小微粒，最終形成介觀尺度的分相狀態。

瀝青在外力作用下的形變情況也相對復雜，在介觀尺度展現出黏彈性。原子力相探測顯微鏡實驗顯示，瀝青中各介觀相在輸入的正弦外力信號作用下會輸出具有不同相位差的應變信號，這說明各相形變中彈性形變和塑性形變貢獻有所區別。常溫下瀝青主要分為3種介觀相：高低錯落相(Catana-phase，一些文獻中稱為蜂型結構，本文根據其載荷形變響應特征稱為高貯能模量相)、周相(Peri-phase)和鄰相(Para-phase)。高貯能模量相結構相對有序，主要由瀝青質組成，載荷下主要發生彈性形變；周相主要由軟瀝青中的芳香化合物構成，圍繞高貯能模量相分布，強度相對較大；鄰相主要由飽和分構成，柔軟且塑性形變特征最明顯[9]。各介觀相受力形變的響應特征不同，增加了宏觀尺度研究瀝青組成與材料性能關系的難度。

基于瀝青質自發相互聚集的特性，瀝青質富集在瀝青材料以彈性形變為主的高貯能模量相中，在材料整體中起彈性骨架作用[10]。通過分子動力學模擬的方法，構建無定形動力學盒子并將不同的瀝青質模型分子分別放入盒中，在模擬真實環境的溫度、壓強下構建瀝青質聚集體模型，模擬實際條件下瀝青質的存在環境。對這些盒子施加應變后統計分子構型和相互作用的變化，實現從分子尺度分析瀝青在形變過程中各種因素的貢獻，為開發更高性能的瀝青產品提供理論基礎，也為后續的研究提供可行的方案。該方法也可以擴展到瀝青其他介觀相分子的形變分析中。

瀝青質分子間的范德華作用和靜電相互作用由分子動力學結果直接得出。但由于模擬體系中的瀝青質分子數量較多且分子結構復雜，如果將所有分子的每一個鍵長、鍵角、扭矩參數全部統計并直接用于描述每個分子的構型，巨大的數據量顯然難以統計，也無法直觀統計應變下眾多分子的形變情況。因此筆者采用統計每個盒子中瀝青質分子的形變相關勢能，通過相關各項勢能值的相對變化對分子的伸縮、鍵角變化和扭曲等形變類別進行統計，從而描述瀝青質分子在實際聚集狀態下隨形變發生時的構型變化。

1 研究對象與分子動力學模型構建

1.1 瀝青質模型分子

瀝青質是瀝青中不溶于正構烷烴(正戊烷或正庚烷)而溶于甲苯的擬組分。瀝青質作為一類相對分子質量在500～1500的稠環芳香化合物，分子種類繁多，因此無法對每種分子進行分離表征，但結構具有一定的共性，可以用Mullins的“島式結構”模型加以概括[11]。瀝青質分子結構可以概括為幾大特征要素，即相對分子質量、稠環連接方式、雜原子種類與含量、側鏈數量(即大小)等，這些結構因素可以概括不同結構類型瀝青質分子的差異。

瀝青質稠環核心的芳環縮合方式分為渺位縮合(cata-Condensed)和迫位縮合(peri-Condensed)2種，縮合方式除了決定稠環核心的構型外，也能影響分子的結構穩定性和稠環的共軛電子分布[12-13]。因此筆者以這一結構特征作為選擇瀝青質模型分子的標準。

以3種瀝青質分子作為模型化合物，分別記作As1、As2和As3，分子結構如圖1所示。3種瀝青質均來自石油瀝青，分子精細結構由Schuler等[14]通過原子力顯微鏡實驗掃描獲得。As1、As2分子符合大氣壓光電離傅里葉變換離子回旋共振質譜表征出的高含量瀝青質的結構特征，As3則相對含量較少，瀝青熒光發射光譜也同樣支持這一結論[15-16]。3種分子均為C50左右。As1的芳環全部為渺位縮合，帶有2個較短的側鏈；As2帶有1個由6個芳環迫位縮合而成的單元，其他芳環均為渺位縮合，并有若干個甲基取代基；As3則為高度迫位縮合的瀝青質，碳/氫原子比相比前2種分子明顯更高，且不含取代側鏈。

圖1 3種瀝青質分子結構模型Fig.1 Model molecules of three kinds of asphaltenes(a) As1—Totally cata-condensed asphaltene; (b) As2—Partly cata-condensed asphaltene; (c) As3—Totally peri-condensed asphaltene

1.2 瀝青質聚集體及形變特點

瀝青分子的不均勻分布有別于傳統的膠束模型，瀝青中各相遵循溶膠-凝膠(Sol-Gel)結構。對于瀝青質而言，4～8個瀝青質分子先形成盤狀的納米聚集體[8,17]，這些聚集體又進一步聚集形成100 nm以下、維數1.7左右的分形簇。介觀尺度上，部分分型簇“懸浮”在周相內部，形成類似膠束的結構；另一些則相互聚集形成連續結構并和其他相相互滲透，形成相邊界復雜的分形結構，HausDorff維數大約在1.7～1.9[10]。分子動力學的模擬尺度可以滿足模擬分形簇的要求，但遠小于分相出現的介觀尺度。

瀝青質所在的相以彈性形變為主，但依然具有黏彈性形變特征，不能簡單將該相形變特征等同于晶體等周期型結構，對比其他經典非晶體體系也具有一定的特殊性。以聚合物聚對苯二甲酰對苯二胺[18]和鋁酸三鈣水泥[19]的分子動力學實驗模擬結果為例進行對比，聚合物應力-應變曲線反映出典型的彈性形變-失效過程，水泥的分子動力學模擬圖像隨著應變發生出現明顯的斷裂層。由于彈性形變和塑性形變同時發生在瀝青質聚集體中，且隨著形變的發生2種形變的貢獻也會發生變化，瀝青質在實驗范圍內隨應變上升不會出現明顯的失效閾值。同時由于應變范圍相對較小，難以從分子模擬圖像上直接觀察到材料整體受到破壞的現象。

1.3 瀝青質形變勢能變化

基于1.2節所述，對瀝青質聚集體受力形變過程的描述應當從分子的形變和分子間相互作用的變化入手。分子動力學模擬基于分子構形或相互作用相關的標準參考值和實際值的差值與能量的擬合公式，對模擬系統的勢能進行計算和統計，換言之，也可以通過勢能反過來對系統內分子總體構形變化進行描述。勢能項又可以分為分子內勢能(由分子構型決定)和分子間勢能(分子間相互作用)。在最簡單的力場中，分子內勢能可以分為鍵伸縮(Bond stretching)、鍵角彎曲(Bending)和二面角扭矩(Torsion)3類，分子間勢能是靜電相互作用和范德華相互作用之和。力場計算還考慮多個分子構型參數之間的相互作用，設立交叉項，這些項是對分子內勢能的修正，在數據處理中可以將這些項整合進3種基本分子內勢能中。同時考慮到sp2雜化的平面構型維持因素，統計出偏離平面構型附加勢能[20]。

在真實瀝青中，上千種瀝青質分子堆積成瀝青質聚集體進而成相。筆者采用勢能來研究具有不同稠環鏈接結構的瀝青質分子在應變過程中分子結構與分子間相互作用的變化，為了突出強調分子結構的影響，在每一個動力學模型盒子中只放置一種模型分子。考慮到瀝青的非晶態相態變化特點，在工作溫度范圍(250～353 K)內可能會出現玻璃態、高彈態甚至到黏流態的相態變化過程，故每種模型分子模擬250、293和353 K 3個溫度下的相態。為了模擬瀝青質在正常工作條件和極大形變情況下的勢能變化，應變分別選擇±0.01、±0.03、±0.05和±0.10共計8個應變值。

含有不同瀝青質的盒子在相同溫度下總勢能不同，含有同樣分子的盒子在不同溫度下總勢能也不相同，因此比較不同分子或溫度下盒子在形變過程中勢能的變化值并沒有太多意義。定義相對勢能變化值(Relative potential energy change)，即應變發生后系統勢能的變化值與應變前系統勢能的比值，來描述勢能隨應變發生的變化，從而可以描述瀝青分子構型和分子間相互作用隨應變發生的變化規律，也能直觀比較相同應變下不同模型分子或溫度系統的分子構型和分子間相互作用的變化情況。

2 實驗條件

采用Material Studio 2017R2進行分子動力學模擬。軟件支持無定形結構材料性質及相關過程的研究，筆者采用其中的Forcite模塊對瀝青質分子形成的介觀微結構進行模擬。利用預設力場參數，通過有限差分法，采用數值積分求解所構建體系中的牛頓運動方程，從而獲取每個分子隨時間變化的坐標、構象、瞬時勢能和動能，依據統計力學計算體系的整體性質。

實驗力場采用COMPASS Ⅱ[21]，相互作用截斷半徑為1.55 nm。動力學模擬分為模型構建和應變施加2個階段。

通過模型構建階段完成設定條件下瀝青質聚集體的模擬。在無定形盒子中分別放入100個相同的瀝青質模型分子作為瀝青質聚集體模型，3種瀝青質模型分子分別對應生成的3個盒子。模型構建階段對每個盒子進行如下流程的動力學模擬：①對盒子進行能量優化；②對盒子進行250～500 K、5次循環的退火過程；③在1 GPa、NPT系綜內對盒子進行體積收縮，每種瀝青質模型盒子分別固定在250、293和353 K溫度進行模擬，整個過程花費500 ps，步長1 fs，選擇Nosé-Hoover恒溫器和Andersen恒壓器；之后分別在相同溫度、正則系綜下進行200 ps、步長0.5 fs的動力學模擬以維持系統穩定，采用Nosé-Hoover恒溫器；此步完成時總計生成9個模型盒子；④各模型盒子的溫度設定與第3步中相同，在1×10-4GPa、NPT系綜內生成瀝青質聚集體模型，總模擬時間5 ns，步長1 fs，選擇Nosé-Hoover恒溫器和Andersen恒壓器；然后在相同溫度、正則系綜下進行200 ps、步長0.5 fs的動力學模擬穩定系統，采用Nosé-Hoover恒溫器。由此獲得瀝青質聚集體的模型。

應變施加階段中，對上一階段優化好的9個模型分別施加±0.01、±0.03、±0.05和±0.1的應變，施加應變后合計72個模型盒子。應變的施加在與各自上一階段相同的溫度、正則系綜下完成，整個過程花費500 ps，步長1 fs，選擇Nosé-Hoover恒溫器。

在勢能項的讀取中，將交叉項整合進鍵伸縮、鍵角彎曲和二面角扭矩3類勢能項中。將伸縮-伸縮交差勢整合到鍵伸縮勢能項，將伸縮-彎曲-伸縮交叉勢、彎曲-彎曲交叉勢和扭矩-彎曲-彎曲交叉勢整合到鍵角彎曲勢能項，將伸縮-扭矩交叉勢、伸縮-扭矩-伸縮交叉勢和彎曲-扭矩-彎曲交叉勢整合到二面角扭矩勢能項。

3 實驗結果

3.1 總勢能變化

圖2是3種模型分子As1、As2和As3在不同溫度下體系相對勢能變化與應變的關系。由圖2可見，3種模型分子聚集體的正應變和負應變的勢能變化圖像并不對稱，可見瀝青質聚集體的應變不能用簡單的彈性形變來描述。當應變超過±0.05時，體系的總勢能常能觀察到異常的變化，顯示為劇烈的上升或異常下降。為了避免應變中材料性能變化或內部結構破壞，應該保證體系內勢能在一定范圍內平穩變化，瀝青的設計應變不超過±0.03的經驗和模擬計算的結果相吻合。因此對相對勢能變化曲線應重點關注±0.03應變區間。

圖2中，在±0.03應變范圍內，250 K下3種分子聚集體模型的相對勢能變化范圍和趨勢相近，顯示出低溫下勢能相對變化和應變在一定范圍內的關系與分子構型無關。在293 K下，As2和As3模型的勢能變化曲線在±0.03應變范圍與各自在250 K下的類似，但As1模型則在293 K下相對勢能變化非常顯著，且隨應變絕對值上升而迅速上升。在353 K下，As1和As2的相對勢能變化曲線較為接近，與各自在293 K下的勢能變化曲線差別明顯；而As3在353 K下相對勢能變化曲線與250 K及293 K下的差別依舊不大。

圖2 瀝青質分子聚集體模型相對勢能變化與應變的關系Fig.2 Relationships between strains and relative potential energy changes for asphaltene aggregation models(a) As1—Totally cata-condensed asphaltene; (b) As2—Partly cata-condensed asphaltene; (c) As3—Totally peri-condensed asphaltene

綜上所述，瀝青質的稠環縮合方式不同，在不同的溫度下其聚集體隨應變勢能上升的趨勢也有所不同。在±0.03應變范圍內，迫位縮合的瀝青質聚集體隨應變勢能上升最不明顯，且溫度對曲線的影響不顯著；渺位縮合的瀝青質聚集體隨應變大幅上升，曲線形狀和趨勢受溫度影響很大；而局部渺位縮合的分子在較低溫度下與迫位縮合分子的曲線相近，高溫(353 K)下則和渺位縮合分子的曲線接近。

3.2 勢能變化類型分析

將瀝青質聚集體模型在應變中相對勢能變化分解為相對分子內勢能變化和相對分子間作用變化，前者也可以稱為相對形變勢能變化(Relative deformation potential energy change)。對各模型中2種相對勢能變化值和應變的關系進行分析。

圖3是3種瀝青質組成的模型在不同溫度下相對形變勢能變化值和應變的關系。由圖3可知，3種分子組成的模型相對形變勢能隨應變變化的大小和規律各不相同，渺位縮合的As1組成的模型相對形變勢能變化最大，As2次之，迫位縮合的As3相對形變勢能變化最小。相對形變勢能變化曲線的變化規律可以概括為2類：第一種是相對形變勢能隨應變絕對值增大而增加，無論是拉伸還是壓縮模型形變勢能均高于原始狀態，曲線呈現U字形趨勢；第二種是相對形變勢能隨應變增大而減小，即分子形變勢能隨壓縮形變增大而升高，隨拉伸形變增大而下降。也有一些形變勢能變化曲線是2種規律疊加的結果。

圖3中，在±0.01應變范圍即瀝青正常工作區內，293 K及以上溫度的As1模型、所有溫度下的As2模型和250 K時的As3模型均表現為U形趨勢，其余則表現為一致下降的趨勢。而在±0.03應變即正常工作應變邊界時，As1在3個溫度下的變化趨勢均與±0.01時的趨勢維持一致，As2在250 K 和293 K溫度下0.01～0.03區間內相對勢能曲線開始下降，而As3模型在250 K和293 K下在0.01～0.03應變區間也開始下降，353 K下卻開始回升?？偨Y來看，3種分子的模型的相對形變勢能變化曲線和溫度有如下關系：當溫度較低(250 K)時，整體趨勢傾向于一致下降的趨勢；而在較高溫度下(353 K)，相對形變勢能變化曲線均往U形趨勢變化。

圖3 瀝青質分子聚集體模型相對形變勢能變化值與應變的關系Fig.3 Relationships between strains and relative deformation potential energy changes for asphaltene aggregation models(a) As1—Totally cata-condensed asphaltene; (b) As2—Partly cata-condensed asphaltene; (c) As3—Totally peri-condensed asphaltene

圖4是3種模型分子組成的聚集體模型的相對分子間相互作用勢能變化值與應變的關系。如公式(1)所示，無應變和應變狀態下分子間相互作用勢能均為負值，因此若發生應變后分子間相互作用勢能上升即分子間相互作用勢能變化值為正，此時相對分子間相互作用勢能變化為負值；且應變后勢能上升越多，相對分子間相互作用勢能變化值越小(該值的絕對值此時會上升)。同理，若發生應變后分子間相互作用勢能下降，此時分子間相互作用勢能變化值為負，相對分子間相互作用勢能變化為正值；且應變后勢能下降越多，則相對分子間相互作用勢能變化值越大。

圖4 瀝青質分子聚集體模型相對分子間相互作用勢能變化值與應變的關系Fig.4 Relationships between strains and relative intermolecular potential energy changes for asphaltene aggregation models(a) As1—Totally cata-condensed asphaltene; (b) As2—Partly cata-condensed asphaltene; (c) As3—Totally peri-condensed asphaltene

(1)

分子間相互作用勢能變化=應變下分子間
相互作用勢-無應變分子間相互作用勢能

(2)

在分子動力學中，兩分子間的相互作用勢能遵循L-J勢能曲線。從圖4可以觀察到，一些曲線(例見圖4(c)中3條相對勢能曲線)呈現倒U形，表明分子間相互作用勢能隨應變的變化和L-J勢能曲線規律相似，倒U形的頂點對應的是分子間相互作用勢能的最低點。另外有一些相對勢能變化曲線單調遞減(見圖4(a)中293 K與353 K的相對勢能曲線)，這些案例意味著整個體系的相對分子間相互作用勢能最低點在應變范圍之外，隨著應變由負到正逐漸增大，體系分子間相互作用勢能逐漸增大。而在±0.03應變區間，所有的相對分子間相互作用勢能變化值均隨應變增大而下降，意味著在這一范圍內分子間相互作用勢能隨應變增大而單調遞增。

綜合來看，瀝青質聚集體模型隨應變絕對值增大勢能上升，是分子間相互作用勢能和分子構形勢能變化的共同結果，且在正應變(拉伸)和負應變(壓縮)下影響不同。隨著模型壓縮程度上升，負應變絕對值逐漸增大，分子相對形變勢能上升(見圖3)，相對分子間相互作用勢能可能略微下降，但不會抵消總勢能的上升趨勢。而隨著模型拉伸，正應變逐漸增大，分子相對形變勢能在不同的模型中上升或先升后降，而相對分子間相互作用勢能大幅上升，使總勢能依舊上升。

3.3 形變勢能類型分析

圖5是3種分子在250、293和353 K下形成的9個聚集體模型中的伸縮、彎曲、扭矩和偏離平面構型附加4項勢能變化。由圖5可知，各項形變勢能均處于規則不明顯的波動狀態，難以用簡單的規律進行描述。盡管如圖3所示，壓縮過程必然增加分子形變勢能，但依舊能從圖5中觀察到一些勢能變化值小于0的點。由此可見模型整體的壓縮并不會同步增大所有的分子形變勢能項，分子形變勢能的變化是多種形變形式共同協調作用的結果。

As1—Totally cata-condensed asphaltene；As2—Partly cata-condensed asphaltene；As3—Totally peri-condensed asphaltene圖5 As1、As2和As3分子聚集體模型在250、293和353 K各類形變勢能變化值與應變的關系Fig.5 Relationships between strains and several deformation potential energy changes of As1, As2 andAs3 models at 250, 293 and 353 K(a) Stretching; (b) Bending; (c) Torsion; (d) Out-of-plane

從伸縮勢能項可以觀察到，壓縮形變時As3模型中該項的變化十分可觀，而As2維持在一個較低的變化水平；拉伸形變時則正好相反，As2拉伸形變勢能變化水平上升，As3的變化反而較低(見圖5(a))。彎曲勢能項變化值由高到低依次是As1、As2、As3，As3分子一直保持在很低的彎曲勢能變化水平上，可見迫位縮合程度越高的瀝青質分子鍵角在應變中變化越小(見圖5(b))。3種分子模型的扭矩勢能變化值顯示，As3模型在應變中的扭矩變化遠小于前兩者(見圖5(c))。綜上所述，渺位縮合分子在聚集體應變過程中的形變主要以分子內部不同平面之間的夾角扭曲來實現，而不僅僅是將化學鍵進行拉伸、壓縮或分子平面鍵角發生變化；迫位縮合的瀝青質分子鍵角彎曲勢能和二面角扭曲勢能隨應變的變化較小，分子整體呈現相對剛性，分子形變勢能反而主要由鍵伸縮貢獻。As1模型的3種形變勢能變化項均最為明顯，這可能是由渺位縮合的稠環鏈接方式和短飽和側鏈2個因素共同所致。

瀝青質的稠環芳核有維持平面構型的趨勢。當整個模型被壓縮時，分子構型會扭曲偏離平面構型，偏離平面附加勢能會有一定程度的上升如圖5(d)所示。所有的模型在不同溫度下均符合隨應變增大而下降的趨勢。在壓縮形變過程中，偏離平面構型附加勢能在形變勢能上升中起到不可忽視的作用。

4 結 論

(1)所有瀝青質聚集體模型隨應變的勢能變化基本呈現U字形趨勢。在瀝青工作的±0.03應變范圍內，全迫位縮合瀝青質聚集體模型的相對勢能變化值最小，全渺位縮合分子聚集體則最大。這意味著全迫位縮合的分子在近似彈性響應的狀態下材料整體更容易發生形變，在設計形變范圍內強度最低；而全渺位縮合分子正好相反。

(2)迫位縮合瀝青質聚集體的相對勢能變化值隨應變的變化趨勢受溫度影響最小，渺位縮合分子聚集體受影響則相對更大。由此可見，迫位縮合瀝青質分子盡管對材料整體強度貢獻較低，但能夠降低材料整體的感溫性，保證在工作溫度范圍材料強度性能變化幅度較小。

(3)瀝青質聚集體模型在應變下的勢能增加主要由形變勢能和分子間相互作用勢能貢獻，兩類勢能隨應變變化的趨勢各有特點。分子間相互作用勢能在±0.03應變范圍內隨應變增加單調上升；分子形變勢能則存在U字形或單調遞減2種可能的趨勢，同一種分子組成的模型在不同溫度下變化趨勢也可能不同。這意味著瀝青質分子在材料拉伸和壓縮狀態下的貢獻是不同的，在一些場合下，材料的形變也可能會使分子平均構型反而更接近能量最低的狀態。

(4)不同構型的瀝青質分子在聚集體應變中的分子形變類型也有所不同。渺位縮合分子主要依靠分子扭矩的變化提高分子形變勢能。迫位縮合分子在壓縮過程中則主要由于鍵伸縮帶來勢能變化，分子扭矩和鍵角彎曲對勢能的影響均不明顯。此外，由于瀝青質分子均有盡力維持平面構型的稠合芳環結構，因此所有分子模型的平面維持附加勢能均隨著應變上升而下降。由此可見，瀝青材料形變時瀝青質分子立體構型提供了可觀的內應力。