重質油穩定性的耗散粒子動力學模擬

關冬，張霖宙，趙鎖奇，徐春明

（中國石油大學（北京）重質油國家重點實驗室，北京 102249）

引 言

重質油在開采、儲運及加工過程中常因體系的穩定性下降而帶來一系列問題。重質油穩定性判定及預測對重質油加工過程具有重要指導意義[1]。傳統的重質油穩定性判定法依賴于實際經驗，如Stankiewicz 等[2]根據現場經驗與實驗觀察，將瀝青質與膠質的比值（asphaltene∕resin）及飽和分與芳香分的比值（saturate∕aromatic）作為判斷重質油穩定性的依據。膠體不穩定指數（C.I.I.）則將瀝青質與飽和分的加和（asphaltene+ saturate）同芳香分與膠質的加和（aromatic+ resin）的比值關系作為判斷重質油穩定性的依據[3-5]。上述經驗方法形式簡單，但存在一定的局限性，如重質油四組分（saturate aromatic resin and asphaltene，SARA）間比例關系不唯一。因此，進一步從理論層面分析，建立基于分子層次信息的重質油穩定性計算方法對于認識其化學本質和工業應用十分重要。

導致重質油穩定性下降的原因是重質油膠體分散體系被破壞，更深一層次的原因是瀝青質的聚沉。因此，研究瀝青質的微觀聚集態，得到重油分子微觀聚集模型，有助于從分子層面揭示重質油穩定性下降的原因。早在20 世紀60 年代，研究者對瀝青質的微觀聚集態就有了一定的認識。其中，Yen[6-8]提出了經典的瀝青質聚集體模型。Mullins[9-11]在Yen 模型的基礎上進一步提出了Yen-Mullins 模型，成為目前最受認可的瀝青質聚集模型。Yen-Mullins 模型認為，瀝青質分子形成最終聚集體的過程中，將經歷納米聚集體（nanoaggregate）和團簇（cluster）兩個聚集層次。當瀝青質濃度很低時，為單體分散狀態；隨著瀝青質濃度的增大，聚集數量增大，最終形成絮凝沉淀[9,12]。除此之外，很多學者對瀝青質聚集結構進行了定量表征[13]。Andreatta 等[14]用超聲波將瀝青質分散于甲苯溶液中，發現1 個納米聚集體約由5 個分子所構成，這與VPO 測量的相對摩爾質量一致。Indo 等[15]采用離心分離法，測量了瀝青質分子納米聚集體的直徑約為2 nm。同時，他們還發現納米聚集體通常由3～8 個瀝青質分子組成。Durand 等[16]采用1H 核磁共振擴散序譜（1H nuclear magnetic resonance diffusionordered spectroscopy，1H NMR DOSY）分析瀝青質納米聚集體時發現，納米聚集體的平均摩爾質量為6900 g∕mol，且主要分布在1500～8500 g∕mol 之間，其平均半徑約為1.56 nm。以上研究對瀝青質微觀聚集態給出了定性聚集構型與定量聚集體結構的描述。然而，重質油體系十分復雜，難以從實驗角度直接觀測到重質油的微觀聚集體結構。因此，研究人員試圖采用分子模擬技術得到重油體系的微觀聚 集 結 構[17]。 Headen 等[18]利 用 分 子 動 力 學（molecular dynamics，MD）模擬得到了平衡狀態下瀝青質在甲苯中的有序聚集結構，即納米聚集體結構。Frigerio[19-20]基于MD 模擬再現了瀝青質分子在不同溶劑中的微觀聚集態。然而，受限于MD 方法的模擬尺度，不少學者嘗試基于介觀尺度分子模擬方法得到重質油體系微觀狀態。Wang 等[21]基于粗粒化分子動力學方法再現了重質油中瀝青質在甲苯和正庚烷溶劑中的自組裝行為。Aguilera-Mercado 等[22]基于Monte Carlo（MC）法在介觀水平上建立了原油中瀝青質和膠質聚集的分子模型。相比于以上方法，耗散粒子動力學（dissipative particle dynamics，DPD）方法模擬尺度更大，且更適用于復雜膠體體系建模，將其應用于重質油體系尤為適宜。張勝飛等[23-30]經過多年研究，在DPD 模擬過程中引入剛體轉動算法，系統探究了重質油的微觀膠體結構與聚集行為。基于DPD 模擬結果，張勝飛等[29]認為瀝青質與膠質的比例存在一個極限，當大于這個極限時，重質油將出現聚沉。Ruiz-Morales等[31]采用不同于張勝飛等的建模方法，通過加大DPD 珠子間的彈性常數并限制其鍵長的方式，得到了瀝青質分子的剛性結構，研究了瀝青質分子在油水界面的擇優取向問題。此外，Rezaei 等[32]和Alvarez 等[33]對瀝青質分子按基團功能劃分，研究了瀝青質對油水乳狀液體系的影響。課題組前期工作中[34]建立了SU-DPD 框架，規范化了石油體系粗粒化建模方法并用于重質油體系，模擬結果符合理論與實驗規律，證明了SU-DPD框架的準確性。

以上研究結果表明，基于DPD 方法可以準確計算重質油的微觀聚集結構，為重質油穩定性的判斷提供基礎。本文基于SU-DPD 框架，建立四組分含量不同重質油體系，計算重質油體系的微觀聚集態，統計瀝青質的聚集率以判斷重質油的穩定性。將瀝青質聚集率與C.I.I.指數及穩定性判定圖對比，以證明模擬體系的準確性，并討論C.I.I.指數及穩定性判定圖的局限性。基于DPD 模擬結果提出改進的重質油穩定性判定圖，為實際重油加工過程提供指導。

1 模擬方法

1.1 DPD基本理論

耗散粒子動力學方法最先由Hoogerbrugge等[35-36]提出，最初應用于復雜流體力學領域，是一種全新的基于粗粒化粒子之間軟相互作用的介觀尺度模擬方法。隨后，Kong 等[37]在耗散粒子動力學模擬體系中引入珠子-彈簧模型，使耗散粒子動力學方法可以模擬大分子體系。Espa?ol 等[38-39]關聯了DPD 中隨機力與耗散力的權函數，證明了兩者之間必須滿足漲落-耗散定理，夯實了DPD 的理論基礎。Groot 等[40]證明了保守力參數與Flory-Huggins 參數之間存在對應關系，為保守力參數的計算提供理論依據。至此，耗散粒子動力學方法已經具備了一定的理論基礎，可以用于模擬如重質油這樣的復雜體系。

傳統耗散粒子動力學方法為了簡化計算采用無量綱單位制，模擬的基本結構單位稱作DPD 珠子，DPD 珠子與分子結構一一對應[35-36]。通過規定珠子的單位質量、單位長度和單位能量可以計算體系中的所有物理量。隨著計算機計算能力的加強，如今可以采用真實單位制以達到更高的模擬精度。耗散粒子動力學方法力場形式相對簡單，每個珠子所受合力由5個分力組成，包括保守力、耗散力和隨機力三種非成鍵力以及彈性力和角作用力兩種成鍵力[35-36,40]：

其中，aij為保守力參數；γij為耗散力參數；σij為隨機力參數；eij為單位向量；ξij=ξji為平均數為0、方差為1 的隨機數。隨機力與耗散力可以通過隨機-耗散定理關聯，如式（5）和式（6）所示，兩者動態平衡以達到對體系控溫的目的[40]。

其中，kB為Boltzmann 常數。一旦溫度T和γij確定了之后，σ隨即確定。研究結果表明，在無量綱單位下，γij= 4.5 往往可以獲得精確的模擬結果，其值的改變對模擬結果影響不大[38-40]。

DPD 模擬過程中最重要的參數為保守力參數，保守力參數可由Flory-Huggins 參數關聯得到。當體系密度接近1 g∕cm3時可以采用式（7）計算保守力參數：

其中，aAB代表A 珠子與B 珠子之間的保守力參數。Flory-Huggins 參數χ可由溶解度系數求得，如式（8）所示。

其中，δA和δB分別為A 珠子和B 珠子的溶解度參數；v為珠子的平均體積。可以看出DPD 模擬體系中除保守力參數外，其余參數均可直接確定，保守力參數及體系的粗粒化方法的確定是DPD 方法準確性的關鍵所在。

1.2 SU-DPD框架

為了能準確模擬重質油體系，課題組前期在傳統DPD 模擬框架下，提出了SU-DPD 框架[34]，用于石油體系粗粒化分子建模及模擬參數計算，本文僅做簡要介紹。SU-DPD 方法中對分子的粗粒化方法提出了完整構建規范，其構建思想源自結構導向集總（structure-oriented lumping，SOL）法[41-42]。在此基礎上，本課題組對SOL 法進行了修改和補充，提出了結構單元耦合鍵電矩陣（structural unit and bondelectron matrix，SU-BEM）法，實現對石油中的復雜分子的準確描述，SU-DPD 法繼承了SU-BEM 方法中分子的劃分方法[43]。SU-DPD 方法中共定義了3類珠子，即烴類珠子、雜原子珠子和溶劑珠子，如圖1（a）所示，不同的字母代表結構不同的珠子。N 類和R 類珠子中碳原子以碳碳單鍵相連，N2～N6 珠子用于組成環烷環結構，R2 和R3 珠子用于組成烷基側鏈結構，數字的大小代表碳原子的個數。雜原子片段珠子S、NI5、NI6 和O 珠子分別對應噻吩、吡咯、吡啶和苯酚環狀分子片段。溶劑珠子中T為甲苯珠子，W 為水珠子，其中一個甲苯分子劃分為一個T珠子，三個水分子劃分為一個W 珠子。隨著模擬體系的不同，SU-DPD 方法中珠子的數目可以進一步擴充。

圖1 SU-DPD珠子與分子表達示意圖Fig.1 Schematic diagram of SU-DPD beads and molecular expression

基于SU-DPD 珠子的分子粗粒化過程示意圖如圖1（b）所示。粗粒化過程從稠合芳香環結構開始，鎖定A6 珠子后，在其周圍間隔地劃分A4 與A2 珠子，以此類推。隨后加入環烷環系珠子，連接方法與稠合芳香環珠子一致。對于鏈狀分子或者烷基側鏈，每3個相鄰碳原子劃分入1個R3珠子，劃分完所有R3 珠子之后若剩余2 個碳原子則劃分為R2 珠子。對于含雜原子分子片段，需要在劃分稠合芳香環及環烷環時一并劃分。

SU-DPD 框架中珠子間保守力參數采用關聯Flory-Huggins 參數獲得。基于DPD 珠子所對應的分子結構可以采用分子動力學模擬計算其平衡構型，進而統計其溶解度參數，得到每個珠子的溶解度參數后可以計算兩兩珠子間的Flory-Huggins 參數χ，進而計算兩兩珠子間的保守力參數，如式（7）、式（8）所示[34]。珠子之間的保守力參數如表1 所示。區別于傳統無量綱DPD 方法，本文采用真實單位DPD 方法，即不同的珠子具有不同的質量與體積，珠子的質量和體積與所對應的結構片段保持一致。模擬過程中的截斷半徑以水珠子為基準，取6.46 ?（1 ?=0.1 nm），與文獻報道的參數設置一致，且水珠子半徑介于所有SU-DPD 珠子之間，較為適宜用作截斷半徑基準[21-28]。耗散力參數取0.08854 amu∕fs。模擬體系采用NVT 系綜，盒子的邊長一般大于10 nm，分子填入量一般大于10000 個。為保持稠合芳香環系的剛性特征，粗粒化珠子之間彈性鏈的彈簧常數取300 kcal∕(mol·?2)（1 cal=4.18 J），平衡距離為3 ?[44]。具體模擬步驟為：（1）體系能量最小化；（2）NVT 系綜模擬；（3）在步驟（2）的最終結構基礎上進行NVT系綜模擬。

表1 DPD珠子間的保守力參數Table 1 Conservative force parameters between DPD beads

1.3 模擬體系

為了得到重質油的微觀聚集結構與分子行為的關系，由于不能窮盡重質油的所有分子結構，過去的基礎性研究工作普遍基于四組分（飽和分、芳香分、膠質和瀝青質）模型展開。采用四組分平均分子結構模型表達重質油體系不僅可以降低體系計算量，還可以較好地與文獻結果對比。因此，本文選用重質油四組分平均分子結構構建重質油體系，所采用的重質油四組分平均分子結構與前期工作一致，符合重質油分析數據與理論模型，如圖2所示[34]。

圖2 重質油四組分粗粒化平均分子結構Fig.2 Coarse-grained model for SARA component

對于重質油的穩定性與其四組分含量的關系，Stankiewicz 等[2]根據現場經驗與實驗觀察，總結了230 種重質油的穩定性規律，最終形成重質油穩定性判定圖，如圖3（a）所示。位于圖中曲線下方的實驗點為穩定重質油，曲線上方的實驗點為不穩定重質油。然而，圖中實驗點僅能確定重質油四組分中兩種組分的比值關系，四組分之間的比值關系不唯一。與此同時，膠體不穩定指數（C.I.I.）也同樣存在四組分比值不唯一的問題，C.I.I.值計算公式如下[3]：

式中,w為質量分數，%。分子部分為瀝青質與飽和分質量分數之和，為膠體不穩定組分；分母部分為芳香分和膠質質量分數之和，為膠體穩定組分。Asomaning 等[3]指出，C.I.I.值小于0.5 時，大部分情況下重質油十分穩定，C.I.I.值大于1.0 時，大部分情況下重質油不穩定，當C.I.I.值在0.5～1.0 時，重質油處于亞穩狀態。

為避免重質油四組分比例不唯一的問題，本文結合兩種方法，提出模擬體系，如圖3（b）所示。該圖是在重質油穩定性判定圖的基礎上加入了等C.I.I.值的一系列曲線，在曲線上取一系列試驗點以覆蓋重質油穩定性判定圖。圖中各系列直線（虛線）的通式為：

圖3 重油穩定性判定圖與模擬體系選取示意圖Fig.3 Diagram of heavy oil stability determination and selection of simulation system

其中，x為wasph∕wres；y為wsat∕warom，與重質油穩定性判定圖中的橫縱坐標一致。

該圖具有以下特點：

（1）當C.I.I.為定值時，所有斜率不同但C.I.I.值相同的曲線過定點（C.I.I.，C.I.I.），稱過定點（C.I.I.，C.I.I.）與原點的直線為固定點直線；

（2）當C.I.I.值大于1.00 時，大多數試驗點位于重質油穩定性判定圖中的不穩定區；

（3）當C.I.I.值小于0.50 時，大多數試驗點位于重質油穩定性判定圖中的穩定區；

（4）當C.I.I.值在0.50～1.00 時，大多數試驗點位于重質油穩定性判定圖中的亞穩區。

以上特點說明，C.I.I.與重質油穩定性判定圖在一定程度上判定效果相當。根據此圖，確定本文所采用模擬體系如表2 所示，共99 組模擬體系。表中的第1～16 組模擬體系對應著固定點直線上的4 個數據點，C.I.I.值從0.25 到1.50，分別取M=0.5，1.0，2.0 和3.0，用于將模擬結果與C.I.I.值進行對比。其余模擬組，為C.I.I.值為0.25～2.00 時，圖中虛線部分間隔一段距離選取的模擬點，用于將模擬結果與重質油穩定性判定圖作比較。

表2 模擬體系Table 2 Simulation systems

2 結果與討論

2.1 DPD計算結果與C.I.I.值的比較

圖4 展示了第1～16 組模擬體系重質油微觀聚集結構，為了便于觀察瀝青質的聚集情況，圖中僅顯示了瀝青質分子。可以看出，隨著C.I.I.值的增大瀝青質的聚集程度增加，這與C.I.I.值的判定效果是一致的。然而，C.I.I.的定義式中沒有考慮到膠質與芳香分的比值對瀝青質聚集率的影響，從重質油的微觀聚集結構中可以觀察到，隨著膠質與芳香分比值的增加，瀝青質的聚集程度增加。

圖4 第1～16組模擬體系重質油微觀聚集結構Fig.4 Heavy oil aggregation structure of simulation systems 1—16

前期工作表明，在當前模擬體系下瀝青質發生聚集時的分子間距離為4.4 ?[34]。因此，本文將分子間距小于4.4 ?作為統計瀝青質聚集率的依據，統計結果如圖5 所示。同種瀝青質聚集率隨C.I.I.值變化曲線表明，瀝青質的聚集率與C.I.I.值成正比，當C.I.I.值大于0.50 時聚集率均大于5%，當C.I.I.值大于1.00時聚集率均大于15%，這與C.I.I.的預測效果也是相符的。可以認為，在當前模擬體系下，瀝青質的聚集率大于15%時，重質油處于不穩定狀態，瀝青質聚集率小于5%時，重質油處于穩定狀態，瀝青質聚集率介于5%與15%之間時，重質油處于亞穩狀態。膠質與芳香分的比值對瀝青質聚集率影響存在極限，在膠質與芳香分比值小于2.0 時，瀝青質的聚集率與膠質與芳香分比值成正比。當體系膠質與芳香分的比值大于2.0 時，隨著膠質與芳香分比值的加大瀝青質的聚集率并未明顯增大。

圖5 瀝青質聚集率隨C.I.I.的變化Fig.5 Change of asphaltene aggregation rate with C.I.I

2.2 DPD計算結果與穩定性判定圖的比較

為了將DPD 計算結果與穩定性判定圖進行比較，本文統計了17～99組重質油微觀聚集結構。圖6展示了具有代表性的重質油微觀聚集結構，可以看出，當瀝青質含量達到一定值時，隨著瀝青質含量的增加，瀝青質聚集率增大，當瀝青質含量不高時，瀝青質的聚集率與四組分的含量有關。統計四組分含量對瀝青質聚集率的影響，如圖7 所示。結果表明，當瀝青質含量達到一定值時，瀝青質的含量直接影響其聚集率，符合Yen-Mullins模型。當瀝青質含量較低時，不可以忽略除瀝青質外的其他組分對瀝青質聚集率的影響。

圖6 第17～99組模擬體系代表性重質油微觀聚集結構Fig.6 Typical heavy oil aggregation structure of simulation systems 17—99

圖7 四組分含量對瀝青質聚集率的影響Fig.7 Effect of SARA content on asphaltene aggregation ratio

為了將模擬結果與重質油穩定性判定圖進行比較，統計了17～99 組模擬體系的瀝青質聚集率如圖8 所示。在低Asphaltene∕Resin 值時，瀝青質的聚集率均比較小，其中完全不聚集的模擬點多位于重質油穩定性判定圖中分界線附近。但是當Asphaltene∕Resin 值增大時，模擬結果與重質油穩定性判定圖存在不一致的情況。模擬結果表明，在高Asphaltene∕Resin 值且位于重質油穩定性判定圖曲線下方的點出現了瀝青質聚集率較高的情況。然而，原始的重質油穩定性判定圖在這一部分沒有試驗數據，說明DPD 模擬結果證明了重質油穩定性判定圖的局限性。從理論角度分析，當Asphaltene∕Resin值較大時，體系的瀝青質濃度和瀝青質聚集率也會較大，說明DPD 模擬結果是合理的。總地來說，DPD 模擬的結果與重質油穩定性判定圖的試驗數據基本吻合，同時彌補了重質油穩定性判定圖的不足。

圖8 第17～99組模擬體系瀝青質聚集率統計結果Fig.8 Statistical results of asphaltene aggregation rate of simulation systems 17—99

2.3 重質油穩定性經驗判據的改進

基于上述討論，發現DPD 模擬結果可以用于判定重質油的穩定性，將瀝青質的聚集率5%作為穩定與亞穩狀態分界點，瀝青質的聚集率15%作為亞穩與不穩定狀態分界點是合理的。DPD 模擬結果可有效補充重質油穩定性判定圖與C.I.I.的不足。因此，基于本文模擬結果提出改進的重質油穩定性經驗判定圖，如圖9所示。圖中兩條曲線，藍線為原重質油穩定性判定圖的曲線，綠線為根據模擬結果統計，平行于藍線畫出的分界線，兩條線將圖中分成3個區域。對于重質油穩定性的判斷可參照下述判定法：

圖9 改進的重質油穩定性判定圖Fig.9 Improved stability judgment diagram of heavy oil

（1）當C.I.I.值小于0.50 且位于圖中Ⅰ區時，重質油體系絕對穩定；

（2）當C.I.I.值大于1.00 且位于圖中Ⅲ區時，重質油體系絕對不穩定；

（3）當C.I.I.值小于1.00且大于0.50或者位于圖中Ⅱ區時，重質油的穩定性不確定，可以采用DPD模擬得到重質油的穩定性情況。

3 結 論

本文基于耗散粒子動力學法模擬了不同四組分含量下重質油體系的微觀聚集結構，并由微觀聚集結構統計瀝青質分子聚集率，用于重質油穩定性的判定。結果表明，基于DPD 模擬得到的瀝青質聚集率與C.I.I.值成正比，證明了本文選用的模擬體系的準確性。將瀝青質聚集率與C.I.I.值對比得到重質油穩定性判定標準，當體系瀝青質聚集率大于15%時，重質油不穩定，瀝青質聚集率小于5%時，重質油穩定，瀝青質聚集率介于5%與15%之間時，重質油處于亞穩狀態。選取83 組模擬結果與重質油穩定性判定圖對比，結果表明，本文模擬體系得到的瀝青質聚集率與重質油穩定性判定圖結果較為吻合，DPD 模擬可有效補充穩定性判定圖的不足。基于DPD 模擬結果，本文提出了改進的重質油穩定性判據，用于重質油穩定性的快速判定，本文提出的重質油穩定性判定方法有望實現工業應用。

符 號 說 明

kB——Boltzmann常數，J∕K

M——膠質與芳香分的質量比

rij——i與j珠子間的距離，?

T——熱力學溫度，K

Δt——時間差，fs

V——珠子體積，?3

vij——i與j珠子間的速度差，?∕fs

γij——i與j珠子間的耗散力參數，amu∕fs

δ——溶解度參數，MPa1∕2

ξij,ξji——隨機變量

σij——i與j珠子間的隨機力參數，amu∕fs

χ——Flory-Huggins參數

ωC,ωD,ωR——權函數