王 鵬，黃世軍，趙鳳蘭，趙大鵬

（1.中國(guó)石油大學(xué)（北京）石油工程教育部重點(diǎn)實(shí)驗(yàn)室，北京 102249；2.中國(guó)石油大學(xué)（北京）石油工程學(xué)院，北京 102249）

近幾十年來，瀝青質(zhì)沉積問題給石油生產(chǎn)、運(yùn)輸和煉油等領(lǐng)域帶來了許多工業(yè)挑戰(zhàn)［1］。瀝青質(zhì)沉積不僅會(huì)造成儲(chǔ)層孔喉堵塞，還會(huì)導(dǎo)致井筒和生產(chǎn)設(shè)備腐蝕結(jié)垢，嚴(yán)重影響油田的正常生產(chǎn)［2-3］。因此，油田開發(fā)過程中瀝青質(zhì)沉積的防治問題亟需解決。

瀝青質(zhì)結(jié)垢的清除成本高昂，如何規(guī)避或抑制瀝青質(zhì)沉積成為研究者們致力去尋找的解決方案。建立瀝青質(zhì)沉積的熱力學(xué)預(yù)測(cè)模型是研究和預(yù)防瀝青質(zhì)沉積的重要手段之一，而對(duì)瀝青質(zhì)宏觀尺度熱力學(xué)行為準(zhǔn)確建模的前提是需要了解其微觀尺度的聚集機(jī)理，這需要明確瀝青質(zhì)分子聚集涉及的分子間作用力［4-5］。然而，由于瀝青質(zhì)的分子結(jié)構(gòu)和官能團(tuán)尚未準(zhǔn)確測(cè)定，瀝青質(zhì)沒有精確的定義和統(tǒng)一的分子構(gòu)型。瀝青質(zhì)是石油中極性最強(qiáng)的重質(zhì)組分，目前應(yīng)用最廣泛的是根據(jù)溶解度將其定義為可溶于芳族溶劑（如甲苯）而又不溶于正構(gòu)烷烴（如正庚烷）的化合物［6］。溶解度定義中提及的瀝青質(zhì)一般特征不足以反映瀝青質(zhì)分子聚集涉及的分子間相互作用，關(guān)鍵還是要獲取具有代表性的瀝青質(zhì)分子構(gòu)型。

隨著計(jì)算科學(xué)理論和計(jì)算機(jī)技術(shù)的進(jìn)步，分子動(dòng)力學(xué)方法得到廣泛的應(yīng)用。分子動(dòng)力學(xué)方法是在分子和原子尺度上描述體系結(jié)構(gòu)與行為的計(jì)算機(jī)模擬方法，不僅可以用來預(yù)測(cè)物質(zhì)的性質(zhì)、材料微觀結(jié)構(gòu)以及微觀作用機(jī)理，還可以解決時(shí)間和成本的問題，具有物理模擬實(shí)驗(yàn)和理論計(jì)算不可比擬的優(yōu)勢(shì)［7］。丁勇杰依據(jù)三組分分離法，選擇基于核磁共振分析和紅外分析的兩種瀝青質(zhì)分子構(gòu)型，采用分子動(dòng)力學(xué)模擬方法，探究瀝青的化學(xué)結(jié)構(gòu)特征［8］；HEADEN 等開展經(jīng)典原子分子動(dòng)力學(xué)模擬研究，對(duì)比4種瀝青質(zhì)分子構(gòu)型、膠質(zhì)分子構(gòu)型及其混合物在甲苯或正庚烷中的聚集特征［9］。瀝青質(zhì)分子結(jié)構(gòu)對(duì)其聚集特征的影響較為顯著，且瀝青質(zhì)分子結(jié)構(gòu)多樣，每一個(gè)研究案例都需要單獨(dú)進(jìn)行建模。因此，有必要建立一個(gè)具有代表性的瀝青質(zhì)平均分子構(gòu)型。為此，采用分子動(dòng)力學(xué)方法，結(jié)合瀝青質(zhì)分子單體代表性結(jié)構(gòu)式、平均相對(duì)分子質(zhì)量和分子基團(tuán)特征，建立瀝青質(zhì)平均分子構(gòu)型，并將其應(yīng)用到瀝青質(zhì)模擬體系初始構(gòu)型的構(gòu)建，通過幾何構(gòu)型優(yōu)化和退火處理后，模擬瀝青質(zhì)在地層環(huán)境中的運(yùn)移過程，探究瀝青質(zhì)分子性質(zhì)和微觀聚集特征，進(jìn)一步揭示瀝青質(zhì)沉積的微觀作用機(jī)理。

1 模擬方法

1.1 瀝青質(zhì)模擬體系的分子動(dòng)力學(xué)實(shí)驗(yàn)原理

分子動(dòng)力學(xué)方法是基于經(jīng)典力學(xué)原理，通過分子體系總勢(shì)能求解系統(tǒng)中所有原子核的牛頓運(yùn)動(dòng)方程數(shù)值解，分析原子核位置和速度隨時(shí)間的演化，來計(jì)算系統(tǒng)的動(dòng)力學(xué)和熱力學(xué)性質(zhì)［10］。量子力學(xué)方法從原子的電子結(jié)構(gòu)和分子的運(yùn)動(dòng)出發(fā)，通過求解不同亞原子粒子的波函數(shù)來描述粒子行為，是最精確的分子勢(shì)能計(jì)算方法［7］。但即使是很小的系統(tǒng)，量子力學(xué)計(jì)算也非常昂貴，且超出了當(dāng)前大多數(shù)計(jì)算機(jī)的計(jì)算能力。因此，研究者們建立了原子力場(chǎng)模型來近似計(jì)算分子勢(shì)能［10］。力場(chǎng)方法是在經(jīng)典力學(xué)的理論基礎(chǔ)上，借助經(jīng)驗(yàn)和半經(jīng)驗(yàn)參數(shù)計(jì)算分子結(jié)構(gòu)和能量的方法，它通過提供一組基于粒子相對(duì)位置的勢(shì)函數(shù)和力學(xué)參數(shù)來計(jì)算系統(tǒng)的勢(shì)能。

現(xiàn)有較為成熟的力場(chǎng)模型主要依據(jù)元素種類或原子所在官能團(tuán)種類對(duì)分子中的各原子賦予不同的力場(chǎng)參數(shù)，計(jì)算分子內(nèi)和分子間的相互作用，進(jìn)而研究分子性質(zhì)［11］。在力場(chǎng)模型中，化學(xué)鍵或所在官能團(tuán)不同的同種元素會(huì)看作是不同種類的原子，這里的種類并非指真實(shí)的原子種類，而是力場(chǎng)模型中表征不同原子特征的分類，如烷烴上的碳、苯環(huán)中的碳和卟啉中的碳在力場(chǎng)模型中隸屬于不同原子，會(huì)賦予不同的力場(chǎng)參數(shù)。本文研究的瀝青質(zhì)模擬體系，既包括有機(jī)大分子瀝青質(zhì)單體，又包括水和二氧化碳等無機(jī)分子，對(duì)于這樣一個(gè)復(fù)雜的有機(jī)-無機(jī)復(fù)合系統(tǒng)，明確分子中原子的成鍵類型、官能團(tuán)類型等結(jié)構(gòu)特征，選擇合適的力場(chǎng)模型賦予相應(yīng)的力場(chǎng)參數(shù)，是保證瀝青質(zhì)模擬體系分子動(dòng)力學(xué)模擬準(zhǔn)確性的基礎(chǔ)。

1.2 模擬體系初始構(gòu)型的建立

采用Accelrys公司研發(fā)的Materials Studio 軟件，建立瀝青質(zhì)平均分子構(gòu)型，并通過Amorphous Cell模塊組裝瀝青質(zhì)模擬體系（圖1）。Amorphous Cell是一種分子動(dòng)力學(xué)無定形組裝模塊，能建立復(fù)雜無定型系統(tǒng)中的代表性模型，以特定密度組裝多個(gè)分子，且各分子在模擬系統(tǒng)中的堆疊合理、仿真程度較高。以密度為1.2 g/cm3組裝4 個(gè)瀝青質(zhì)分子、48個(gè)水分子和10 個(gè)CO2分子，模擬瀝青質(zhì)在地層中的運(yùn)移環(huán)境，以二氧化硅分子層代表地層孔隙骨架，模擬盒尺寸為26.60 ?×9.83 ?×12.99 ?，在各方向上均采用周期性邊界條件。

圖1 地層環(huán)境下瀝青質(zhì)模擬體系的初始構(gòu)型Fig.1 Initial configuration of asphaltene simulation system in formation environment

1.3 力場(chǎng)參數(shù)的確定

在進(jìn)行分子動(dòng)力學(xué)模擬實(shí)驗(yàn)前，首先需要選擇合適的力場(chǎng)模型，并賦予相應(yīng)的力場(chǎng)參數(shù)來計(jì)算粒子之間的相互作用，而力場(chǎng)參數(shù)的選擇對(duì)分子動(dòng)力學(xué)模擬結(jié)果的準(zhǔn)確性至關(guān)重要。不同的力場(chǎng)模型適用的模擬體系有所不同，如AMBER（Assisted Model Building with Energy Refinement）力場(chǎng)主要適用于模擬生物大分子，而CHARMM（Chemistry at Harvard Macromolecular Mechanics）力場(chǎng)可以較好地?cái)M合小分子體系和溶劑化的大分子體系。因此，針對(duì)地層環(huán)境下瀝青質(zhì)模擬體系的特點(diǎn)，從以下3 個(gè)方面考慮來確定力場(chǎng)參數(shù)。

瀝青質(zhì)分子的自聚集行為 瀝青質(zhì)作為原油中極性最強(qiáng)的組分，具有很強(qiáng)的自聚集傾向，不同濃度的瀝青質(zhì)分子的聚集形態(tài)有所不同［4］。當(dāng)質(zhì)量分?jǐn)?shù)較低（小于10-4）時(shí)，瀝青質(zhì)以分散態(tài)的單分子形式存在；隨著體系中瀝青質(zhì)質(zhì)量分?jǐn)?shù)的不斷增加，瀝青質(zhì)分子開始發(fā)生聚集行為，聚集顆粒逐漸變大，依次形成納米聚集體（通常含有8～10個(gè)瀝青質(zhì)分子）、納米簇和黏彈性網(wǎng)絡(luò)［12］。在選擇與瀝青質(zhì)模擬體系相適應(yīng)的力場(chǎng)模型時(shí)，首先要考慮的是能夠較好地?cái)M合瀝青質(zhì)這種有機(jī)大分子自聚集行為。

CVFF（Consistent Valence Force Field）力場(chǎng)可適用于計(jì)算各種多肽、蛋白質(zhì)和大量的有機(jī)分子，在計(jì)算系統(tǒng)結(jié)構(gòu)與結(jié)合能方面最為準(zhǔn)確，可以較好地?cái)M合瀝青質(zhì)分子自聚集行為和瀝青質(zhì)-膠質(zhì)分子締合行為［13］。其勢(shì)能函數(shù)表達(dá)式為：

有機(jī)-無機(jī)復(fù)合體系 瀝青質(zhì)模擬體系中除了有機(jī)大分子瀝青質(zhì)，還包括水和二氧化碳等無機(jī)小分子，因此，需要選擇一種能兼顧有機(jī)-無機(jī)復(fù)合體系的力場(chǎng)模型。COMPASS（Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies）力場(chǎng)采用從頭計(jì)算與經(jīng)驗(yàn)方法相結(jié)合的方式計(jì)算力場(chǎng)參數(shù)，是第一個(gè)把有機(jī)分子體系的力場(chǎng)與無機(jī)分子體系的力場(chǎng)統(tǒng)一的力場(chǎng)模型［14］。其勢(shì)能函數(shù)公式為：

（2）式中右側(cè)前5 項(xiàng)為成鍵作用勢(shì)能，后2 項(xiàng)為非成鍵作用勢(shì)能，分別為靜電相互作用和范德華力，其計(jì)算公式分別為：

分子力場(chǎng)的遷移性問題 力場(chǎng)參數(shù)是通過大量擬合基分子的實(shí)驗(yàn)數(shù)據(jù)得到的，若擬合的實(shí)驗(yàn)數(shù)據(jù)是在常溫常壓條件下測(cè)定的，而模擬的瀝青質(zhì)在地層中的運(yùn)移環(huán)境往往是高溫高壓條件，此時(shí)應(yīng)用這些力場(chǎng)參數(shù)時(shí)分子力學(xué)的準(zhǔn)確性會(huì)降低［15］。Dreiding 力場(chǎng)是一種比較特殊的力場(chǎng)計(jì)算方法，可以計(jì)算一些缺乏實(shí)驗(yàn)數(shù)據(jù)的新結(jié)構(gòu)化合物，是典型的經(jīng)驗(yàn)力場(chǎng)模型。通過借助Dreiding 力場(chǎng)中的經(jīng)驗(yàn)參數(shù)進(jìn)行輔助參考，有助于提高分子動(dòng)力學(xué)模擬實(shí)驗(yàn)結(jié)果的準(zhǔn)確性。

1.4 瀝青質(zhì)分子動(dòng)力學(xué)模擬的實(shí)驗(yàn)方法

瀝青質(zhì)分子動(dòng)力學(xué)模擬的實(shí)驗(yàn)方法主要包括以下4 個(gè)步驟：①基于本文建立的瀝青質(zhì)模擬體系初始構(gòu)型，選擇合適的力場(chǎng)參數(shù)來計(jì)算系統(tǒng)中粒子間的相互作用。②采用最陡下降法調(diào)整偏離平衡位置較大的結(jié)構(gòu)，再基于共軛梯度法對(duì)接近平衡狀態(tài)的結(jié)構(gòu)進(jìn)行優(yōu)化，使系統(tǒng)內(nèi)分子達(dá)到能量最小值和原子位置最佳分布，實(shí)現(xiàn)瀝青質(zhì)模擬體系的幾何構(gòu)型優(yōu)化。③在高溫條件（303.15～503.15 K）下進(jìn)行8個(gè)循環(huán)的退火處理，消除殘余應(yīng)力，退火過程持續(xù)8 ps。④采用玻爾茲曼分布函數(shù)，隨機(jī)賦予系統(tǒng)中各分子運(yùn)動(dòng)的初始速度（圖2），使用Forcite 模塊在NVT 系綜下模擬瀝青質(zhì)在地層巖石孔隙中的運(yùn)移過程，模擬地層溫度設(shè)為328.15 K，控溫方法選擇Nosé-Hoover 算法，時(shí)間步長(zhǎng)為1 fs，模擬時(shí)長(zhǎng)為15 000 ps。

圖2 系統(tǒng)中各分子運(yùn)動(dòng)初始速度的玻爾茲曼分布Fig.2 Boltzmann distribution of initial velocity of each molecule in system

2 結(jié)果與討論

2.1 瀝青質(zhì)平均分子構(gòu)型的構(gòu)建

瀝青質(zhì)分子的主體由含有若干個(gè)芳環(huán)組成的片狀單元結(jié)構(gòu)所構(gòu)成，這些片狀單元結(jié)構(gòu)的個(gè)體間結(jié)構(gòu)差異較大，根據(jù)芳核特征可將瀝青質(zhì)分子結(jié)構(gòu)分為大陸型和群島型［16］。由不同地區(qū)原油中瀝青質(zhì)的分子構(gòu)型（圖3）可以看出，不同地區(qū)的瀝青質(zhì)分子結(jié)構(gòu)存在差異，如遼河和華北原油中的瀝青質(zhì)分子的脂肪族支鏈較長(zhǎng)，塔里木原油中的瀝青質(zhì)分子則表現(xiàn)為標(biāo)準(zhǔn)的大陸型結(jié)構(gòu)，而伊朗和俄羅斯原油中瀝青質(zhì)分子的芳環(huán)數(shù)量較多［17］。雖然不同地區(qū)的瀝青質(zhì)分子結(jié)構(gòu)存在差異，但總體上看其結(jié)構(gòu)特征也具有一般性，這為建立代表性的瀝青質(zhì)平均分子構(gòu)型提供了參考。

圖3 不同地區(qū)原油中的瀝青質(zhì)分子構(gòu)型Fig.3 Molecular configurations of asphaltenes in crude oil samples from different regions

目前，已知的大多數(shù)瀝青質(zhì)分子單體結(jié)構(gòu)包括1 個(gè)由4～10 個(gè)芳環(huán)組成的芳核和周圍長(zhǎng)度為3～7個(gè)碳原子的脂肪族鏈，這為構(gòu)建瀝青質(zhì)平均分子構(gòu)型的芳環(huán)數(shù)和支鏈長(zhǎng)提供了參考依據(jù)［4］。ROGEL提出的一種基于沉淀法測(cè)定的委內(nèi)瑞拉油樣瀝青質(zhì)餾分的平均結(jié)構(gòu)模型［18］被后續(xù)研究者沿用較多，其結(jié)構(gòu)式具有一定的代表性，是本文瀝青質(zhì)平均分子構(gòu)型構(gòu)建的基礎(chǔ)。而平均相對(duì)分子質(zhì)量會(huì)影響瀝青質(zhì)的密度和溶解度等物理性質(zhì)，可以作為構(gòu)建瀝青質(zhì)平均分子構(gòu)型的約束條件。隨著測(cè)定儀器精度的提高和測(cè)量方法的改進(jìn)，瀝青質(zhì)平均相對(duì)分子質(zhì)量的測(cè)定值發(fā)生了幾個(gè)數(shù)量級(jí)的變化，從最初測(cè)定的103～106到現(xiàn)在普遍認(rèn)可的750［19］。綜合多種瀝青質(zhì)分子結(jié)構(gòu)式、平均相對(duì)分子質(zhì)量（750）和瀝青質(zhì)分子基團(tuán)結(jié)構(gòu)的一般特征，構(gòu)建出如圖4 所示的瀝青質(zhì)平均分子構(gòu)型，該分子構(gòu)型的H/C 值為1.2，而瀝青質(zhì)分子的H/C 值一般為1.0～1.5，因此其結(jié)構(gòu)具有一定的代表性。

圖4 瀝青質(zhì)平均分子構(gòu)型Fig.4 Average molecular configuration of asphaltenes

2.2 瀝青質(zhì)平均分子構(gòu)型優(yōu)化與退火處理

利用商業(yè)分子設(shè)計(jì)軟件構(gòu)建分子構(gòu)型的原理往往不是第一性原理，而是依靠分子內(nèi)固有的化學(xué)鍵的鍵長(zhǎng)和鍵角，但這些數(shù)據(jù)往往是通過大量實(shí)驗(yàn)測(cè)得的平均數(shù)據(jù)，因此構(gòu)建的初始分子構(gòu)型往往不是目標(biāo)化合物的穩(wěn)定結(jié)構(gòu)。為了得到更加真實(shí)的分子幾何構(gòu)象，在分子動(dòng)力學(xué)中通常采用能量極小化方法，即在分子間作用力處于最小的狀態(tài)下，調(diào)整分子幾何構(gòu)象使分子鍵長(zhǎng)和鍵角接近“自然”狀態(tài)，進(jìn)而得到原子位置的最佳排布［10］。

從能量?jī)?yōu)化結(jié)果（圖5）可以看出，系統(tǒng)能量總體上隨著優(yōu)化時(shí)間步數(shù)的增加而降低，但變化過程中存在明顯的轉(zhuǎn)折點(diǎn)。其原因是：由于勢(shì)能面具有波動(dòng)性，數(shù)學(xué)上只能保證求得局部極小值，需要采用多種構(gòu)象搜索方法結(jié)合使用，來盡可能地達(dá)到全局最小值［20］。本文首先采用最陡下降法，其原理是根據(jù)能量梯度負(fù)值進(jìn)行一維搜索，是能量極小化初期最簡(jiǎn)單、快捷的搜索算法，適用于分子結(jié)構(gòu)偏離平衡位置較遠(yuǎn)［10］。但隨著不斷逼近極小值，其搜索方向總是垂直于上一步的特性使得搜索到的勢(shì)能值總是在極小值附近震蕩而無法收斂，因而可以看到優(yōu)化時(shí)間步數(shù)為0～380時(shí)能量衰減趨勢(shì)緩慢且波動(dòng)頻繁，此時(shí)需要采用更加精細(xì)的共軛梯度算法，通過不斷迭代，快速逼近系統(tǒng)能量的全局最小值，該過程對(duì)應(yīng)于優(yōu)化時(shí)間步數(shù)為380～500時(shí)系統(tǒng)能量衰減速度顯著增加，直到達(dá)到最低勢(shì)能構(gòu)象，即最穩(wěn)定的構(gòu)象。

圖5 瀝青質(zhì)模擬體系初始構(gòu)型能量?jī)?yōu)化Fig.5 Energy optimization of initial configuration of asphaltene simulation system

對(duì)瀝青質(zhì)模擬體系初始構(gòu)型進(jìn)行幾何構(gòu)型優(yōu)化后，完成系統(tǒng)溫度從303.15 K 到503.15 K 的8 個(gè)循環(huán)退火處理，退火過程中系統(tǒng)的各種能量變化如圖6 所示。由于退火過程中系統(tǒng)溫度波動(dòng)，總動(dòng)能和總勢(shì)能的變化也存在一定波動(dòng)，但相比于初始構(gòu)型，隨著系統(tǒng)溫度升高，分子熱運(yùn)動(dòng)加劇，系統(tǒng)總動(dòng)能和總勢(shì)能增加。隨著分子間距增加，分子間的范德華力相互作用減弱，靜電相互作用將占據(jù)系統(tǒng)非成鍵作用中的主導(dǎo)地位。在高溫條件下，瀝青質(zhì)分子內(nèi)各原子熱運(yùn)動(dòng)加劇，能夠充分暴露其內(nèi)部結(jié)構(gòu)，隨后通過退火處理，系統(tǒng)溫度緩緩降低，在每個(gè)溫度都達(dá)到結(jié)構(gòu)的平衡態(tài)，如此循環(huán)往復(fù)退火過程，消除殘余應(yīng)力，使瀝青質(zhì)分子結(jié)構(gòu)達(dá)到能量基態(tài)。

圖6 退火過程中系統(tǒng)能量變化Fig.6 Energy change during annealing

2.3 瀝青質(zhì)分子構(gòu)型的徑向分布特性

徑向分布函數(shù)是表征粒子微觀結(jié)構(gòu)特征的物理量，其物理意義為任意指定的中心粒子周圍距離為r處出現(xiàn)其他粒子的概率，如圖7所示。它反映了凝聚態(tài)物質(zhì)粒子聚集的徑向有序性，其特點(diǎn)為“以宏表微”，即只需測(cè)定瀝青質(zhì)的宏觀數(shù)值就可以得到其微觀細(xì)節(jié)［21-23］。

圖7 徑向分布函數(shù)概率分布示意Fig.7 Probability distribution of radial distribution function

徑向分布函數(shù)可以理解為系統(tǒng)的局部密度與平均體密度的比值，其計(jì)算式［23］為：

分子動(dòng)力學(xué)模擬獲得的系統(tǒng)各種統(tǒng)計(jì)學(xué)平均性質(zhì)，主要通過統(tǒng)計(jì)分析函數(shù)來計(jì)算軌跡文件中關(guān)于各種性質(zhì)的數(shù)據(jù)并繪制分布圖，如總能量、總勢(shì)能、總動(dòng)能、速度和密度等隨時(shí)間的變化。g(r)同樣也可以從軌跡數(shù)據(jù)文件中獲得，其實(shí)質(zhì)就是原子間矢量長(zhǎng)度的球型平均分布［24］。對(duì)非晶體而言，當(dāng)r值很大時(shí)，g(r)的函數(shù)值趨近于1，表明距離參考粒子為r的位置的粒子分布無規(guī)律性；而當(dāng)r值相對(duì)較小時(shí)，g(r)的函數(shù)值會(huì)在小范圍內(nèi)出現(xiàn)幾個(gè)波峰（極大值），即這些位置處的局部密度與平均體密度的比值較大，說明波峰位置處出現(xiàn)粒子的概率要明顯高于其他位置?？傮w上，非晶體的徑向分布特征表現(xiàn)為近程有序、遠(yuǎn)程無序［25］。

從瀝青質(zhì)平均分子構(gòu)型的徑向分布函數(shù)（圖8）可以看出，g(r)總體上呈現(xiàn)近程波動(dòng)劇烈且分散、遠(yuǎn)程波動(dòng)平緩且密集的特征。當(dāng)r＜0.25 nm 時(shí)，g(r)表現(xiàn)出較強(qiáng)的規(guī)律性，分別在0.11，0.14，0.17 和0.22 nm等處出現(xiàn)波峰，這些波峰表明從參考粒子到該位置處出現(xiàn)粒子的概率要遠(yuǎn)高于其他位置，主峰位置一般標(biāo)志著聚集結(jié)構(gòu)的出現(xiàn)，這可能與瀝青質(zhì)分子的稠環(huán)結(jié)構(gòu)有關(guān)；而當(dāng)r＞0.5 nm時(shí)，g(r)的函數(shù)值趨近于常數(shù)1，說明粒子的分布沒有規(guī)律，呈無序狀態(tài)。瀝青質(zhì)是典型的非晶體，在成鍵作用下分子內(nèi)各原子呈有序排列狀態(tài)，表現(xiàn)為近程有序；隨著r值增大，瀝青質(zhì)的分子聚集狀態(tài)方位會(huì)對(duì)瀝青質(zhì)的出現(xiàn)頻率產(chǎn)生影響，在常見的方位配置中，具有非面對(duì)面配置方式的瀝青質(zhì)分子徑向分布規(guī)律性較差，其出現(xiàn)頻率無規(guī)則，分子堆積更松散，因此遠(yuǎn)程表現(xiàn)為無序狀態(tài)。

圖8 瀝青質(zhì)平均分子構(gòu)型的徑向分布函數(shù)Fig.8 Radial distribution function of average molecular configuration of asphaltenes

2.4 瀝青質(zhì)在地層巖石中的運(yùn)移模擬

瀝青質(zhì)在地層環(huán)境運(yùn)移過程中，系統(tǒng)內(nèi)各分子運(yùn)動(dòng)的初始速度由玻爾茲曼分布函數(shù)隨機(jī)產(chǎn)生，模擬過程中系統(tǒng)溫度在設(shè)定值328.15 K上下波動(dòng)。從瀝青質(zhì)模擬體系在地層中的運(yùn)移終態(tài)（圖9）可以看出：大量水分子向模擬盒的兩端運(yùn)移，并最終傾向于吸附到二氧化硅分子層表面，這驗(yàn)證了二氧化硅的強(qiáng)親水性；而二氧化碳分子在整個(gè)分子動(dòng)力學(xué)模擬過程中位移量較小，這可能與其是非極性分子有關(guān)，且目前文獻(xiàn)中也尚未有記載關(guān)于從分子間的相互作用的角度去解釋二氧化碳對(duì)瀝青質(zhì)沉積的直接影響；瀝青質(zhì)分子中的芳核主體以及部分脂肪族鏈均不同程度的向二氧化硅分子層表面靠近，并傾向于逐漸被吸附到二氧化硅分子層表面，瀝青質(zhì)分子單體或聚集體在二氧化硅表面的吸附造成孔隙內(nèi)壁的非均質(zhì)性會(huì)改變接觸角，是導(dǎo)致孔隙潤(rùn)濕性改變最可能的原因之一，這一點(diǎn)在ANDREW 等［26］和LI 等［27］的實(shí)驗(yàn)研究中都得到了驗(yàn)證。巖石表面的潤(rùn)濕性由親水性向親油性轉(zhuǎn)變，是瀝青質(zhì)沉積導(dǎo)致原油采收率降低的重要機(jī)理之一。

此外，從圖9 中亦可觀察到瀝青質(zhì)分子間的聚集狀態(tài)，聚集的驅(qū)動(dòng)力來自于瀝青質(zhì)分子中芳核的π-π 堆積作用，而脂肪族支鏈引起的位阻則是瀝青質(zhì)分子聚集的阻力。π-π 堆積作用是由于芳香族環(huán)的負(fù)π 電子云和正σ框架之間的靜電吸引而產(chǎn)生的，常見的堆積方式包括面對(duì)面或平行、錯(cuò)位堆疊或平行錯(cuò)位、以及邊對(duì)面或T形［4］。其中，圖10中所示的堆積方式為面對(duì)面堆積，是最穩(wěn)定的苯二聚體構(gòu)型，它能夠最大程度地減少負(fù)π 電子云之間的相互排斥作用。瀝青質(zhì)分子中芳核發(fā)生的π-π 堆積作用引發(fā)了瀝青質(zhì)分子聚集，最終導(dǎo)致瀝青質(zhì)絮凝并從原油中沉積出來，堵塞巖石孔喉并導(dǎo)致原油采收率降低。

圖9 瀝青質(zhì)模擬體系運(yùn)移終態(tài)Fig.9 Final migration state of asphaltene simulation system

圖10 瀝青質(zhì)分子二聚體構(gòu)型示意Fig.10 Configuration of asphaltene molecular dimer

3 結(jié)論

基于分子動(dòng)力學(xué)方法，結(jié)合多種瀝青質(zhì)分子結(jié)構(gòu)式、平均相對(duì)分子質(zhì)量和瀝青質(zhì)分子基團(tuán)結(jié)構(gòu)的一般特征，建立了瀝青質(zhì)平均分子構(gòu)型。該構(gòu)型具有一定的代表性，并用于構(gòu)建瀝青質(zhì)模擬體系初始構(gòu)型來模擬瀝青質(zhì)在地層巖石中的運(yùn)移。

瀝青質(zhì)作為典型的非晶體，其徑向分布特征總體表現(xiàn)為近程有序、遠(yuǎn)程無序。當(dāng)周圍概率粒子到參考粒子的距離較小時(shí)，在成鍵作用下分子內(nèi)各原子呈有序排列狀態(tài)；隨著周圍概率粒子到參考粒子的距離增大，具有非面對(duì)面配置方式的瀝青質(zhì)分子徑向分布規(guī)律性較差，分子堆積更松散。退火處理過程中系統(tǒng)的各種能量變化表明，退火優(yōu)化可使分子鏈充分松弛，有利于模擬該分子結(jié)構(gòu)的能量初態(tài)。隨著分子間距增加，分子間的范德華力相互作用減弱，靜電相互作用將占據(jù)系統(tǒng)非成鍵作用中的主導(dǎo)地位。

從瀝青質(zhì)模擬體系的運(yùn)移終態(tài)可以看出，瀝青質(zhì)分子的芳核主體及部分脂肪族鏈傾向于逐漸被吸附到二氧化硅分子層表面，這會(huì)使巖石表面的潤(rùn)濕性由親水性向親油性轉(zhuǎn)變；瀝青質(zhì)芳核的π-π 堆積作用是引發(fā)瀝青質(zhì)分子聚集的驅(qū)動(dòng)力，從而導(dǎo)致瀝青質(zhì)絮凝并從原油中沉積出來。潤(rùn)濕性反轉(zhuǎn)和巖石孔喉堵塞是瀝青質(zhì)沉積造成原油采收率降低的重要機(jī)理，明確瀝青質(zhì)分子的微觀聚集特征，能夠有效預(yù)防原油開采過程中潛在的儲(chǔ)層傷害風(fēng)險(xiǎn)，尤其是在注氣提高采收率技術(shù)應(yīng)用過程中。

符號(hào)解釋

a，Dr——鍵伸縮能計(jì)算參數(shù)；

Aij——范德華相互作用參數(shù)，代表原子間的排斥作用；

b，b0——分子的鍵長(zhǎng)和平衡鍵長(zhǎng)，nm；

Bij——范德華相互作用參數(shù)，代表原子間的吸引作用；

dN——球形殼內(nèi)的粒子數(shù)目，個(gè)；

dr——球形殼的厚度，nm；

g(r)——徑向分布函數(shù)；

i，j——體系中任意兩個(gè)不同原子；

kθ——鍵角彎曲能計(jì)算參數(shù)；

K?——鍵扭轉(zhuǎn)能計(jì)算參數(shù)；

kx——鍵角面外彎曲能計(jì)算參數(shù)；

n——鍵扭轉(zhuǎn)能、鍵角面外彎曲能參數(shù)；

qi——原子i的電荷量，C；

qj——原子j的電荷量，C；

r——周圍概率粒子到參考粒子的距離，nm；

rij——原子i與原子j之間的相對(duì)距離，nm；

——原子i與原子j間范德華勢(shì)能為0時(shí)的距離，nm；

U——總勢(shì)能，kJ/mol；

Ub——鍵伸縮能，kJ/mol；

Ucross——耦合能，kJ/mol；

Ux——鍵角面外彎曲能，kJ/mol；

Uθ——鍵角彎曲能，kJ/mol；

U?——鍵扭轉(zhuǎn)能，kJ/mol；

Ucoul——靜電勢(shì)能，kJ/mol；

Uvdw——范德華勢(shì)能，kJ/mol；

x，x0——分子的離平面振動(dòng)角度、平均離平面振動(dòng)角度；

θ，θ0——分子的鍵角、平衡鍵角；

?，?0——分子的二面角、平衡二面角；

ρ——密度，g/cm3；

ε0——介電常數(shù)；

εij——?jiǎng)葳迳疃取?/p>