石墨烯改性瀝青界面力學性能的分子動力學模擬*

周 昆，黃 君，鄧雅丹，黃立新，2

(1.廣西大學 土木建筑工程學院，南寧 530004；2.廣西大學 工程防災與結構安全教育部重點實驗室，南寧 530004)

0 引 言

瀝青作為工程領域常用的材料之一，廣泛應用于道路鋪裝、橋面鋪裝等。因其工作環境長期處于室外暴露環境下，極易受到溫度、日照、水損等多方面的影響。因此常添加其他材料對瀝青進行改性，以改善瀝青的路用性能。改性瀝青是向瀝青基體加入改性劑或其他外摻劑的方法，常用的瀝青改性劑包括橡膠、樹脂、碳納米材料、各類聚合物等。改性劑的加入可以提高瀝青彈性、軟化點等多種特性。石墨烯作為近幾年來熱門的碳納米材料，其優異的力學、熱學、電學性能獲得了研究者的廣泛關注[1]。但是由于石墨烯為二維材料，并且其彎曲剛度極小，這使得它極易產生褶皺和缺陷等問題。因此單獨將其作為結構材料具有較大困難。而石墨烯比表面積大的特點使得其能與其他材料充分結合，故近年來學者們開始研究石墨烯在復合材料中的運用[2]。研究發現，石墨烯對瀝青的改性能夠有效提高瀝青在高溫、低溫、抗老化方面的表現[3]，并且能夠提高瀝青路面的耐變形性能、減少車輛動荷載下的車轍效應[4]。但石墨烯作為單純瀝青改性劑，在微觀層面的穩定取決于石墨烯與瀝青之間的界面關系。界面的缺陷會嚴重影響瀝青的力學性能和破壞強度[5]，因此在改性瀝青的眾多性能中，界面力學性能的研究顯得尤為關鍵。

在瀝青復合材料界面性能實驗研究方面已經取得了許多成果。文獻[6-7]利用光學與電子顯微鏡對多孔瀝青混凝土微觀界面形狀進行了直接觀察。Liu[8]使用掃描電鏡、紅外光譜等方法研究了瀝青與鋼渣骨料的黏附、界面微觀結構等界面相關性能，從物理吸附和化學反應兩方面分析了瀝青與鋼渣骨料的黏附機理。Yin[9]使用動態剪切流變儀測試了瀝青與花崗巖的剪切流變性能，并利用X射線光電子能譜分析了二者之間的黏附性能。此外，文獻[10-11]采用通用試驗機和氣動黏附拉力試驗機進行改性瀝青的拉拔實驗，研究瀝青與集料在水分存在下的界面張力破壞。

雖然實驗能給出宏觀狀態下瀝青與瀝青混合料的相關粘結強度，但實驗費時費力，并且所測量到的界面行為受到了增強材料性質的影響，如孔隙率、水分、純度等。并且瀝青與骨料之間的界面黏附行為依賴于不同分子基團之間的分子相互作用[12]。因此近年來，文獻[13-14]采用分子動力學模擬(molecular dynamic，簡稱MD)計算結構的性質。

因此，本文采用MD模擬方法研究石墨烯改性瀝青的分離行為，旨在探討不同因素對石墨烯改性瀝青的界面力學性質影響。首先通過MD得到瀝青的物理性質，并與其他實驗進行對比確認模型的正確性。之后通過對模型施加速度來達到拉拔實驗的目的，進而得到界面應力分離曲線。并與宏觀實驗下得到的內聚區模型進行比較。最后分析了加載速率對石墨烯改性瀝青界面破壞的影響。

1 分子動力學模擬方法

1.1 模 型

瀝青是一種高度復雜的材料，包含105～106種不同的分子[15]。因此選取能夠體現瀝青物理化學性質的分子模型是MD模擬的基礎。Li[16]選擇了美國戰略公路研究計劃(strategic highway research program)的3種瀝青體系：AAA-1、AAK-1、AAM-1。它們分別代表來自Canadian/Lloydminster(AAA-1), Venezuelan/Boscan(AAK-1)和 USA/West Texas(AAM-1)的原油來源。在對這3種瀝青系統進行元素分析、分子類別選取、以及各組分質量分數適配后，確定了瀝青的4組分模型。4個組分分別為：瀝青質(asphaltene)、飽和烴(saturates)、環烷芳烴(Naphthene aromatics)、極性芳烴(Polar aromatics)。這4個組分總共包含12個分子，單個分子結構如圖1所示。因為該模型在研究瀝青結合料內聚性能和粘合性能上具有明顯優勢[17]，所以本文選用該模型進行石墨烯/瀝青界面力學性能的研究。瀝青模型各個分子的信息列于表1。

圖1 瀝青四組分模型

表1 瀝青四組分模型參數

石墨烯單胞尺寸為a=0.246 nm，b=0.426 nm，c=3 nm，α=β=γ=90°，如圖2(a)所示。為建立石墨烯/瀝青復合材料界面模型，石墨烯片尺寸須同瀝青尺寸相近。因此，導入石墨烯單胞模型后，以A=16，B=10，C=1的超晶胞范圍(supercell range)建立超晶胞。并以(0 0 1)為切割面法向量進行切割，從而得到所需要的單片石墨烯表面，如圖2(b)所示。二維石墨烯片尺寸為u=3.935 nm，v=4.256 nm，θ=90°。然后添加1.0 nm真空層將石墨烯片從二維變為三維。最后聯合石墨烯與瀝青建立界面模型。

圖2 模擬采用的石墨烯片參數

1.2 力場

石墨烯與瀝青之間的作用力通過力場來定義。本研究采用OPLS-AA力場(Optimized Potentials for Liquid Simulations all-atom)。該力場對有機液體、生物大分子等物質有更為精確的擬合效果[18-20]。OPLS-AA力場包含了價態項和非鍵作用。價態項主要包括鍵伸縮、角彎曲、二面角扭轉。非鍵作用由靜電相互作用和van der Waals相互作用組成。OPLS-AA力場各種相互作用的勢函數形式為：

Etotal=Ebonds+Eangles+Edihedral+Enon-bonds

(1)

(2)

(3)

(4)

(5)

式中，Kr、Kθ、Vn和φ分別表示與鍵長、鍵角、二面角有關的經驗參數。共價鍵伸縮勢Ebonds和鍵角彎曲勢Eangles采用諧振子勢函數，二面角扭曲勢Edihedral取Fourier展開式前4項。van der Waals相互作用采用Lennard-Jones 12-6勢函數，采用Lorentz-Berthelot混合規則計算交叉項參數。靜電相互作用采用庫倫勢函數，當i，j為1，4時，fij=0.5(表示分子間和分子內1，4相互作用)，否則取1.0。本文中van der Waals相互作用與靜電相互作用的截斷半徑均取1.0 nm，超出10 ?的庫侖力使用PPPM(particle-particle particle-mesh)方法進行計算。

1.3 模擬設置

宏觀狀態下，石墨烯改性瀝青處于非均勻狀態，石墨烯隨機分布在瀝青基體當中，如圖3(a)所示。石墨烯/瀝青復合材料界面的復雜性主要是由于瀝青無定形性質的結果，并且沒有特定的模型來描述微觀狀態下瀝青與石墨烯平衡時的界面狀態。對這樣的體系進行分析時，不可能將所有夾雜結構進行分析。為降低模擬復雜程度，因此須取一個如圖3(b)所示的代表體積單元(representative volume element,簡稱RVE)。該微元界面處的分離行為可以看作復合材料下有效的界面本構關系，實際上是以均勻的微元來替代宏觀狀態下非均勻狀態。在本項工作中，因僅研究石墨烯改性瀝青界面分離響應，無關其他物理化學性質，因此實驗假定石墨烯是平坦無缺陷的。并且石墨烯與瀝青基體之間不存在化學鍵，最終所取RVE如圖4所示。

圖3 代表體積單元的選取

圖4 石墨烯改性瀝青代表體積單元

復合材料界面粘結強度主要由拉拔實驗確定，通過石墨烯與瀝青分離的過程可以確定材料破壞機制。因此，在基于RVE的基礎上，使用MD模擬方法模擬拉拔實驗，以研究原子尺度上石墨烯與瀝青分離的機制。確定RVE以后，為了考慮除RVE外的其他地方的影響，模擬時在xyz3個方向均為周期性條件。運用周期性邊界條件的缺點是沒有足夠距離實現界面脫離，因此需要在拉伸方向添加足夠的真空層，本文中在拉伸方向添加了2.0 nm真空層。

建立好模型后，使用CG(conjugate gradient)方法進行能量最小化處理。然后對模型進行足夠的弛豫，弛豫時間取500 ps。為了得到模擬溫度下的構型，在NVT系綜下使用Nose-Hoover恒溫器進行500 ps的弛豫。拉拔模擬時忽略石墨烯變形，因此將其設置為剛體。為避免石墨烯同瀝青一同運動，需要將底部一定范圍內的瀝青進行固定。給石墨烯設置一個恒定速度向上移動，從而實現分離界面分離。整個分離過程記錄石墨烯所受到的牽引力。通過將該牽引力除以石墨烯面積得到所需要的拉拔應力。

MD模擬中，過小的時間步長會導致計算資源的浪費，過大則會導致積分計算不穩定。文獻[21]指出，最合適的時間步長是通過NVE的能量波動來選定的。因此在拉拔模擬過程中，以1 fs的時間步長對運動方程進行積分，并選擇使用阻尼系數為10 fs的langevin恒溫器控制體系溫度。阻尼系數以時間為單位，決定著體系松弛速度。本研究選取10 fs作為阻尼系數，這代表著體系溫度每10 fs松弛一次，這樣能很好地控制溫度。

瀝青與石墨烯建模均在Material Studio中進行，并結合LigParGen[22]在線網站確定瀝青在OPLS-AA力場下的參數。另外，石墨烯力場參數取自文獻[23]。之后采用自行開發的MATLAB程序進行文件改寫，并輸入至LAMMPS(large-scale atomic/molecular massively parallel simulator)進行運算。模型與整個界面脫離過程使用可視化軟件OVITO進行觀察。

2 瀝青的物理性質

界面拉拔測試的首要前提是確定模型的正確性，因此首先確定模型密度是否與真實瀝青接近。首先以0.1 g/cm3的初始密度構建瀝青無定形模型。之后使用NPT系綜(p=1 atm，T=298 K)對該模型進行弛豫。弛豫過程中，隨著模擬時間的增長，瀝青晶胞體積縮小，分子間間隙逐漸減小，密度也因此上升。密度變化如圖5所示。SHRP給出的AAA-1瀝青的密度在常溫常壓下為1.03～1.04 g/cm3。本研究所得密度為1.006 g/cm3。瀝青模型與真實模型密度上的差距，是由于MD模擬中僅考慮了12種分子類型，與真實瀝青相比少了許多其他組分。密度僅相差0.024～0.034 g/cm3，屬于合理范圍內。

圖5 NPT系綜(T=298 K，p=1atm)下瀝青密度

采用比體積曲線方法確定瀝青玻璃化轉換溫度Tg。首先對體系進行能量最小化處理，而后使用NPT系綜分別得到不同溫度下的穩定構型。取208～478 K的溫度區間，以10 K為間隔記錄了10個溫度下的瀝青比體積，并根據不同溫度下繪制比體積散點圖，如圖6所示。利用最小二乘法確定兩條斜率不同的回歸方程，兩條回歸方程交點即為所求Tg。計算得到的Tg為299 K，該溫度將溫度-比體積曲線劃分為兩個區間，左側區域內比體積變化較緩，右側區域內則變大。本文計算得到的Tg與文獻的結果比較情況列于表2。從表中可以發現文獻[24]使用彎曲梁流變儀測得的Tg為223～303 K，說明Tg屬于一個溫度范圍，而不僅僅是一個特定的點，這主要取決于原油來源和測量方法。只要所測得的結果在該范圍內，都是正確的。綜合模擬所得密度與Tg，并與實驗值對比后，驗證了瀝青模型的正確性以及OPLS-AA力場的適用性。

圖6 GMA玻璃化轉換溫度的確定

表2 瀝青玻璃化轉換溫度的文獻對比

3 拉拔過程中的界面行為

3.1 界面分離過程

拉拔模擬的結果如圖7所示，橫坐標定義為石墨烯片向上移動的距離，縱坐標定義為石墨烯的分離應力。采用移動平均濾波法進行應力數據的收集，每100步(即100 fs)統計一次應力數據，每5個數據進行一次數據平均。該例子記錄了溫度為298 K時的分離情況，石墨烯片移動速率為10 m/s。記錄了20萬步的分離情況，這意味著石墨烯一共移動2.0 nm的距離。拉拔2.0 nm后，石墨烯與瀝青已完全分離破壞。a點為初始狀態，此時石墨烯片應力大于0，這是石墨烯與瀝青之間相互粘結產生的吸引力控制的。而后隨著基底向上移動，應力顯著增加。曲線到b點(x=0.09 nm)后，達到峰值應力，隨后在[0.09 nm，0.2 nm]區間內開始振蕩。可以觀察到此時瀝青還未發生破壞，瀝青間的作用力仍保持著瀝青的完整性。隨著位移增大，瀝青與石墨烯片之間相互作用力減弱(c點)，界面應力迅速減小，開始振蕩并最終減小為0(d點)。

圖7 石墨烯改性瀝青的界面分離響應

當溫度為常溫，加載速率為10 m/s時，可以發現當石墨烯開始脫離瀝青基體時，瀝青首先是被拉伸，之后內部出現微小空隙，但這些空隙并未引起瀝青基體內部結構失效。隨著石墨烯片遠離瀝青，而不對其產生作用力時，此時定義為復合材料的破壞。因此10 m/s加載下，瀝青與石墨烯界面破壞的機制是界面粘結破壞，而非內聚破壞。

3.2 內聚規律

模擬得到的應力-位移曲線與準脆性材料的內聚區模型(cohesive zone model，簡稱CZM)所給出的粘結規律相似。CZM是對瀝青破壞機制參數化的一種模型，對于描述瀝青的破壞機制給出了詳細解釋[29-30]。理想的CZM將材料斷裂限制在兩個零厚度的界面之間，其余區域內是沒有破壞的。作用在裂紋上的牽引力隨界面的分離距離而產生非線性變化，這與本研究所得結果一致。因此可以將MD所得規律同跨尺度模型結合在一起進行考慮。

CZM可以采用線性、雙線性、梯形、光滑梯形、指數等不同的方法來描述[30]。本研究選取指數模型進行非線性擬合，因為這一模型與實驗數據貼合較好，且能更好地反映峰值點過后軟化過程。以應力最大值點σc，以及發生最大應力時的分離距離xc作為參數進行設置。函數模型如下：

圖8給出了常溫常壓下(p=1 atm,T=298 K)使用MD模擬所得出的結果，以及使用指數內聚區模型擬合所得的曲線。擬合后得到應力在分離距離為0.1662 nm時達到最大值，即647.95 MPa。可以發現，指數內聚區模型與MD模擬所得數據吻合較好，內聚規律認為界面破壞發生在界面應力達粘結強度時。在達到峰值后，瀝青進入軟化狀態。隨著損傷的累積，界面應力逐漸降低，并減小至0。

圖8 擬合MD數據得到的內聚區模型

特別地，內聚規律其實是一種不可逆的規律。這一不可逆規律表現在加載-卸載-再加載這一過程中，再加載后所達到的界面應力值無法達到初次加載時的界面應力。這一過程用MD模擬所得結果如圖9所示。圖9(a)給出了在未達到界面峰值應力之前進行卸載的過程，可以看到再加載時峰值應力與初加載時的應力相近。圖9(b)展示在達到界面峰值應力之后進行卸載的過程，可以發現再加載后無法達到原有的峰值應力。這意味著在界面應力達到粘結強度后發生了損傷，瀝青內部結構的損傷導致其無法再保持原有強度。當然，這里需要指出，MD模擬因時間尺度和模型尺度上的限制，加載速率高于實驗所能達到的速率，所以得到的界面應力遠遠超出實驗值[31]。因此MD模擬僅能從微觀角度解釋瀝青界面破壞機制，之后如要進行跨尺度研究，需要將微觀機制與實驗測量結合起來才能深入了解瀝青的軟化過程。

圖9 內聚規律

3.3 不同加載速率下的界面行為

對于宏觀狀態下的瀝青路面，工作過程中所受力的大小不一樣，因此本小節討論不同加載速率對破壞機制的影響。考察10、20、30 m/s等3種加載速率的影響，即每一步(t=1 fs)石墨烯片分別移動1×10-5、2×10-5、3×10-5nm的情況下，石墨烯改性瀝青的破壞過程。圖10給出了不同加載速率下的應力-位移曲線圖。結果表明，低速率下(10 m/s)，即使在對數據進行移動平均處理后，曲線振蕩仍然較為明顯。主要是因為低速率拉拔時，仍然有部分瀝青粘結在石墨烯上。并且應力達到峰值后，有較長一段的分離過程，屬于延性破壞。而在高速率(20和30 m/s)下，曲線基本沒有波動，到達一定位置界面應力就下降為0，這意味著界面的徹底破壞，屬于界面粘結破壞。其次，30 m/s的破壞位置相較于20 m/s的破壞位置更為提前，并且到達峰值點的速度更快，應力峰值更大。這一趨勢與宏觀尺度下的實驗得到的內聚規律一致[32]。以上3種速率加載時，均為界面粘結破壞，即破壞時是石墨烯與瀝青的分離。

圖10 不同加載速率下的拉拔過程

進一步考慮更低加載速率時石墨烯改性瀝青的破壞情況，取5 m/s的加載速率。圖11展示了石墨烯改性瀝青的破壞過程，從圖中可以看到低速率下所測得的最大拉伸應力會降低，且應力-位移曲線振蕩幅度更加顯著。這是由于石墨烯位移的增加使得瀝青內部空隙逐漸增多，其余未破壞部分的瀝青熱運動就更加顯著，從而大部分瀝青產生較大的蠕變變形，并且黏附在石墨烯片上的瀝青比高速率下黏附的要更多。不同于高速率下的界面粘結破壞，低速率下的拉拔破壞則為更嚴重的內聚破壞。

圖11 低加載速率下(5 m/s)的應力分離響應

4 結 論

采用分子動力學模擬方法對石墨烯改性瀝青的界面力學性能進行了研究。首先通過模擬得到常溫常壓(p=1 atm,T=298 K)下AAA-1瀝青密度為1.006 g/cm3，與實驗值相差不大。并通過計算不同溫度下的比體積，使用最小二乘法進行線性擬合得到瀝青玻璃化轉換溫度Tg=299 K。與實驗值對比，Tg=299 K在正確的范圍之內，驗證了瀝青模型的正確性。在此基礎上，開展了石墨烯改性瀝青分子動力學模擬的拉拔實驗，得到如下結論：

(1)通過MD模擬得到10 m/s速率下的應力-位移關系，發現瀝青與石墨烯的破壞屬于界面粘結破壞，并且這一關系同宏觀實驗下得到的內聚規律趨勢一致，應力先隨位移增大而迅速上升，當達到峰值應力后則開始緩慢下降。而這一緩慢下降過程(軟化過程)與指數型內聚區模型吻合最好。

(2)在達到界面粘結強度之前進行卸載時，石墨烯改性瀝青保持彈性狀態，此時的變形可逆。但在達到界面粘結強度之后卸載時，因瀝青內部出現損傷，導致瀝青變形不可逆，無法達到之前的粘結強度。

(3)通過對不同速率的拉拔實驗得知，高加載速率下石墨烯改性瀝青的界面破壞屬于界面粘結破壞。而當加載速率下降到一定程度時，瀝青在拉拔過程中出現大塊空隙，形成孔洞，屬于內聚破壞。

從上述方面可以確定，MD模擬所得結果與實驗值偏差屬于容許范圍內，并且可以與更大高尺度模型相結合，為全面分析石墨烯改性瀝青的界面力學性能提供更加全面的解釋。