分子動力學模擬汽油組分在不同脫硫膜材料中的擴散行為

侯影飛， 蔣 馳， 劉 敏， 劉熠斌， 李 鵬， 牛青山

(1.中國石油大學 重質油國家重點實驗室 化學工程學院， 山東 青島 266580;2.新疆應用職業技術學院 石油與化學工程系， 新疆 奎屯 833200)

?

分子動力學模擬汽油組分在不同脫硫膜材料中的擴散行為

侯影飛1,2， 蔣 馳1， 劉 敏1， 劉熠斌1， 李 鵬1， 牛青山1

(1.中國石油大學 重質油國家重點實驗室 化學工程學院， 山東 青島 266580;2.新疆應用職業技術學院 石油與化學工程系， 新疆 奎屯 833200)

為探究滲透氣化脫硫膜材料脫硫性能的決定因素，采用分子動力學(MD)模擬了典型汽油組分(噻吩、正庚烷、環己烷、環己烯和甲苯)在乙基纖維素(EC)、聚乙二醇(PEG)和聚二甲基硅氧烷(PDMS)中擴散過程，繪制小分子在聚合物中擴散軌跡并分析了其擴散方式。計算求得不同分子在不同膜材料中的擴散速率。通過分析擴散分子與聚合物體系的相互作用能、擴散分子三維尺寸、聚合物自由體積和聚合物鏈段運動對分子擴散速率的影響，考察了分子在聚合物中的微觀擴散機理。結果表明，在3種膜材料中噻吩的擴散速率均大于其余4組分，與實驗結果相一致；擴散分子與聚合物相互作用是噻吩類硫化物在脫硫膜材料中優先擴散的決定因素，且其余3種因素相互協同同樣影響著其在聚合物中的擴散速率。

分子動力學； 分子擴散； 汽油脫硫； 膜分離

隨著世界各國對環境保護的日益重視，以及環保法規的日益嚴格，生產低硫及超低硫汽油已成迫切需求[1]。FCC汽油是汽油的主要來源，其所含的硫化物以硫醇、硫醚、二硫化物和噻吩類這4種有機硫化物為主，其中噻吩類硫含量占總硫含量的60%以上，因此，FCC汽油脫硫主要以脫除噻吩硫為主[2]。與常規加氫脫硫技術相比，滲透氣化汽油脫硫技術具有投資和操作費用低、可深度脫硫、辛烷值損失小等優點，在實驗室和工業應用層面都取得了長足的進展[3]。目前，有關滲透氣化汽油脫硫技術的研究主要集中在對膜材料的選擇、評價和工藝開發，對汽油組分在膜材料中的微觀擴散行為和機理以及硫化物在脫硫膜材料中優先擴散的決定因素的研究較少。

隨著分子模擬技術的發展，其已成為研究高分子膜材料的分子結構以及滲透組分在膜中溶解擴散過程的重要工具[4]。分子動力學(MD)可以在分子水平上精確計算小分子在高分子膜中的擴散性。Li等[5]采用分子動力學模擬研究了苯和水在聚二甲基硅氧烷(PDMS)和杯[4]芳烴(CA)共混膜中的擴散行為，分析了不同共混比例膜體系中PDMS和CA相互作用能、聚合物鏈段運動性、膜材料自由體積，并計算了苯和水在共混體系中的擴散系數。結果表明，分子擴散速率不僅與膜自由體積有關，而且受到CA與擴散分子的相互作用力的影響，模擬結果與實驗相符。Ling等[6]采用分子動力學模擬了甲基萘和二苯并噻吩在不同膜材料中的擴散速率，結果與實驗相一致；考察了影響分子擴散速率的因素，結果表明，擴散分子與膜材料相互作用能是影響分子擴散速率的最主要因素。Amani等[7]采用分子動力學模擬研究了不同氨基甲酸酯基/脲基比例的聚氨酯-聚脲(PUUs)膜的納米結構，并考察了 O2、N2、CO2、CH4、H2S氣體分子在不同PUUs膜中的分離性能，結果與實驗相一致。然而到目前為止，對多汽油組分模型在聚合物中擴散行為的分子模擬研究鮮見報道。

筆者選用乙基纖維素(EC)、聚乙二醇(PEG)和聚二甲基硅氧烷(PDMS)3種常用于汽油膜分離脫硫研究的材料[8]，采用分子動力學模擬技術研究噻吩、正庚烷、環己烷、環己烯和甲苯5種汽油組分在不同膜材料中的擴散行為，觀察分子的擴散方式，并通過Einstein方程計算汽油組分在不同聚合物中的擴散速率；通過考察擴散分子與不同聚合物體系的相互作用、分子三維尺寸、聚合物自由體積和聚合物分子鏈段運動性，探究影響分子在聚合物中擴散性的因素。

1 汽油組分在膜材料中擴散行為的模擬方法與模型構建

分子動力學方法是以經典力學為基礎的模擬計算方法。建立合適體系后，對體系中每個原子求解牛頓運動方程，計算各原子受力及加速度來確定每個時刻各原子的位置及其速度，進而通過體系運動軌跡計算宏觀物理性質[6]。筆者使用Materials Studio的Amorphous Cell和Forcite模塊進行分子動力學模擬。模擬過程中使用COMPASS力場，體系氫鍵計算使用DREIDING力場，能量最小化使用Smart方法，靜電相互作用力計算選擇Ewald方法，范德華作用力選擇Atom-based方法[6]。

分別以EC、PEG和PDMS與5種典型汽油組分混合構建聚合物體系。不同汽油組分按照實驗常用模擬汽油比例[9]填充，且在3種聚合物體系中汽油組分占總質量的15%。按照初始密度0.6 g/cm3構建10個晶胞結構并進行優化，選取能量最低結構，在298 K、101.325 kPa下進行200 ps的NPT系綜動力學模擬，使體系接近真實密度；然后對體系進行200 ps的NVT系綜動力學模擬使結構弛豫，再對其在101.325 kPa下進行500 ps的NPT系綜退火循環加速體系平衡，溫度從298 K到498 K，溫度梯度為25 K。退火后的結構幾何優化后進行500 ps的NPT系綜動力學模擬，使體系達到298 K下真實密度，晶胞結構如圖1所示，分子模型采用球棍結構；繼續升溫到330.15 K，依次進行500 ps的NPT動力學模擬和500 ps的NVT系綜動力學模擬，使體系在目標溫度下充分平衡，然后進行2500 ps的NVT動力學計算，模擬汽油組分擴散過程。模擬數據列于表1。

圖1 3種聚合物體系的分子模型

表1 3種聚合物體系分子模擬數據

2 結果與討論

2.1 小分子在聚合物中的擴散方式

筆者以噻吩作為擴散分子，通過Materials Studio腳本分析動力學軌跡文件，獲得噻吩分子每一幀的重心坐標，繪制出噻吩分子在一個聚合物晶胞體積內的運動軌跡，如圖2所示。模擬結果揭示了小分子在聚合物中的2種運動形式，一是小分子在聚合物空穴范圍內震蕩，二是小分子在聚合物空穴間的跳躍式遷移。從圖2可以看出，噻吩分子先在空穴內震蕩，當噻吩分子有足夠的能量，或者當聚合物鏈運動產生了合適的路徑后，發生空穴間遷移。通過這種方式，噻吩分子在聚合物中進行傳質擴散。分子擴散軌跡所展現的擴散方式與通常認為的分子在聚合物中跳躍擴散的形式相一致[7]。從圖2還可以看出，噻吩分子在EC體系中跳躍次數最少，大部分在固定的空穴中震蕩，可以認為在EC中分子更難以擴散。

圖2 噻吩分子在3種聚合物中的擴散軌跡

2.2 小分子在聚合物中的擴散速率

在MD模擬中，分子擴散速率常采用Einstein方法計算，Einstein法將擴散速率(D)與分子的均方位移(MSD)相關聯，如式(1)[5]所示。

(1)

式(1)中，(ri(t)-ri(0))2為i組分的均方位移，ri(t)為t時間分子i的坐標，ri(0)為分子i的初始坐標。該方程假設分子在聚合物中自由運動，MSD與D呈線性關系，相應的擴散速率是均方位移對時間的變化率。通常，在分子動力學模擬的初始階段會發生非正常擴散，分子擴散受限；在足夠長的擴散時間后，擴散分子從非正常擴散轉變為正常擴散。是否達到正常擴散可以通過(ri(t)-ri(0))2∝tn檢驗，n為時間t的指數；n<1表明在聚合物內小分子為非正常擴散，n=1表明在聚合物內小分子為正常擴散[10]。圖3為噻吩在3種膜材料中MSD與t雙對數曲線。可以看出，隨著模擬時間增長，斜率n趨近于1，表明噻吩在3種聚合物膜材料中均由非正常擴散轉變為正常擴散，因此2500 ps的模擬時間對筆者考察的膜材料體系來說足夠長。通過式(1)求取汽油組分在3種聚合物膜材料中的擴散速率，結果示于圖4。這3種聚合物滲透氣化膜的脫硫性能列于表2。

圖3 噻吩在不同聚合物體系中MSD與模擬時間(t)的對數曲線

從圖4可以看出，在PEG和PDMS體系中，噻吩、正庚烷、環己烯、環己烷、甲苯分子的擴散速率依次降低，而在EC體系中擴散速率降低的順序變為噻吩、甲苯、正庚烷、環己烯、環己烷。在EC、PEG和PDMS體系中，噻吩分子的擴散速率均大于其余4組分，表現出優先擴散性，且其在EC中擴散選擇性最大，與表2所示EC、PEG和PDMS滲透氣化膜的脫硫性能相符合[8]。同一擴散分子在EC體系中擴散速率最小，在PDMS體系中擴散速率最大，與圖2觀察分子擴散軌跡所得結論相一致。

圖4 模擬汽油組分在不同聚合物體系中的擴散速率(D)

表2 實驗所得3種聚合物膜材料的脫硫性能

1) Composed ofn-heptane, cyclohexene, cyclohexane, toluene and thiophene; 2) The sulfur enrichment factor, which represents the selectivity of sulfur removal.

2.3 影響小分子在聚合物體系中擴散的因素

2.3.1 擴散分子與聚合物相互作用能

小分子在聚合物體系中的擴散過程中，其與體系之間的非鍵相互作用具有很大的影響，可以通過擴散分子與體系的相互作用能ΔE來表征，如式(2)所示。

ΔE=(Emol+Esys)-Etol

(2)

式(2)中，Emol為擴散分子能量，Esys為去除擴散分子后混合體系的能量，Etol為總混合體系能量，單位均為kJ/mol。ΔE越大表明相互作用能越強，分子擴散離開接觸面需要克服更高的能量勢壘[6]。

本研究中，擴散分子與聚合物僅是物理混合，并無化學鍵作用，因此使用COMPASS力場通過公式(2)計算得到擴散分子與體系非鍵相互作用能ΔE，結果列于表3。從表3可見，在3種聚合物膜材料中，噻吩的ΔE均低于其余4種擴散分子的，噻吩擴散離開接觸面需要的能量最低，最易于擴散，這與計算的擴散速率結果相符，表明擴散分子與聚合物相互作用是噻吩類硫化物在脫硫膜材料中優先擴散的決定因素，與Ling等[6]所得結論相一致。

表3 擴散組分與聚合物的相互作用能(ΔE)

在3種聚合物膜材料中，噻吩、環己烷、環己烯和甲苯均體現出ΔE越大分子擴散速率越低的規律，而正庚烷在EC和PDMS中的ΔE最大，但其擴散速率卻不是最低；就同一擴散分子在3種聚合物體系中的ΔE而言，與PDMS的ΔE最大，而其在PDMS體系中擴散速率卻大于在EC和PEG體系中的(見圖4)，說明ΔE是影響擴散速率的重要因素，但不是單一決定因素。

2.3.2 擴散分子尺寸

擴散分子的尺寸是影響分子擴散速率的重要因素。隨著擴散分子尺寸的增大，其發生空穴間遷移的幾率減小，傳質的阻力增大[11]。利用Connoly Surface計算了擴散分子范德華體積和表面積，分子體積大小依次是正庚烷、環己烷、甲苯、環己烯、噻吩。然而分子體積僅能體現分子總體積的大小，不能描繪分子的空間形態和結構。為了更精確、全面地表達分子的形狀和尺寸，使用Corey-Pauling-Koltun(CPK scale：1.0)模型模擬分子最低能量構象，將其作三維長方體模型處理。圖5為正庚烷最低能量構象三視圖。用分子在3個垂直方向上的長度a×b×c表示分子的三維尺寸，且a

由表4可以看出，雖然正庚烷的體積最大，但其最小截面尺寸明顯小于其余4種分子，所以正庚烷沿著最小截面垂直方向擴散時阻力更小，需要的自由體積更小，更易于擴散[12]。雖然正庚烷擴散時需要克服的能量勢壘要大于環己烷和環己烯，但其最小截面尺寸遠小于環己烷和環己烯，擴散時有更多的空間和路徑，因此表現出的擴散速率大于環己烷和環己烯。可以看出，擴散分子尺寸和分子與體系相互作用能2種因素相互協同，共同影響分子在聚合物中的擴散性。

圖5 正庚烷最低能量構象三視圖

表4 擴散分子尺寸

1) Molecular three-dimensional dimension; 2) Minimum cross-sectional dimension

2.3.3 聚合物自由體積

根據自由體積理論[13]，聚合物總體積由占有體積和自由體積構成。聚合物自由體積孔穴的分布和大小對擴散分子在體系中的擴散行為具有重要影響[6]。普遍采用正電子湮沒技術來測定材料的自由體積大小，然而實驗方法經常存在較大誤差，且無法精確描繪出自由體積的微觀形貌。分子動力學方法采用特定原子作為探針來探測致密膜的自由體積，結果較為準確[5]。體系內被探測出的自由體積占總體積的比例稱為自由體積分數(FFV)，可由式(3)計算。

(3)

式(3)中，Vf和Vo分別為聚合物的自由體積和聚合物鏈的占據體積，nm3。

膜材料中體積較小的自由體積孔穴對膜的傳遞性質影響很小，能夠容納擴散分子的自由體積空穴對分子擴散有較大影響[14]。筆者利用噻吩為探測分子，用動力學直徑為0.53 nm的硬球表示，考察不同膜材料容納噻吩分子的自由體積分布及其所占比例。根據式(3)求得EC、PEG和PDMS自由體積分數分別為0.739%、0.217%和2.574%。聚合物自由體積空穴的模擬形貌示于圖6，其中灰色和藍色部分為自由體積空穴。可以看出，EC和PEG體系中自由體積空穴分散，而PDMS中自由體積空穴較為連續、緊密。聚合物的FFV越大，能為噻吩提供的擴散空間就越大， FFV從大到小的聚合物體系依次為PDMS、EC、PEG，而噻吩在這3種聚合物中的擴散速率大小卻依次為PDMS、PEG、EC。雖然EC具有比PEG更多擴散空穴，但還要考慮擴散分子在聚合物中擴散方式，在EC中擴散分子主要在空穴內震蕩，很難發生空穴間遷移。可見，不同膜材料自由體積大小并不是同一分子擴散速率大小的決定因素。

圖6 3種聚合物自由體積空穴的模擬形貌

2.3.4 聚合物分子鏈運動性

聚合物分子鏈運動一定程度上決定了自由體積空穴間通道的形成。聚合物分子鏈運動越強，小分子越容易發生空穴間躍遷，表現出越高的擴散速率。聚合物分子鏈的運動性強弱也可以用分子鏈的自擴散速率(Ds)表示[15](見表1)。PDMS、PEG、EC 3種聚合物的Ds依次降低，與擴散組分在這3種材料中擴散速率大小順序一致。雖然EC的自由體積分數比PEG大，但是EC分子鏈運動性遠小于PEG，EC中自由體積之間很難形成有效通道，擴散分子在固定空穴中震蕩，導致分子在EC膜中擴散速率小于在PEG中的。自由體積和聚合物分子鏈運動的重要性也從側面驗證了分子在空穴間跳躍擴散的擴散方式。

EC分子鏈的運動性遠差于PEG和PDMS的，為分析造成這種差異的原因，筆者模擬計算了這3種材料的氫鍵作用力和玻璃化轉變溫度(見表1)。通常情況下，氫鍵是最強的分子間作用力，雖然其小于化學鍵，但其對材料的微觀形貌起重要作用[7]；玻璃化轉變溫度(Tg)體現聚合物分子鏈的軟硬程度，是對比聚合物分子鏈運動性的重要指標[7]，玻璃化轉變溫度越低，聚合物分子鏈運動性越強。計算結果表明，EC膜體系中氫鍵作用力遠大于PEG和PDMS的，EC分子間大量氫鍵制約其分子鏈運動；對比聚合物玻璃化溫度與其分子鏈運動性發現，在體系模擬溫度(330.15 K)未達到EC的玻璃化轉變溫度時，EC分子鏈幾乎不運動，導致很難形成有效的分子擴散路徑。

3 結 論

(1) 采用分子動力學方法和COMPASS力場模擬了噻吩、正庚烷、環己烷、環己烯和甲苯分子在EC、PEG和PDMS膜脫硫材料中的擴散方式，并計算了分子的擴散速率。結果表明，噻吩在EC、PEG和PDMS膜材料中均表現出優先擴散性，與這3種材料表現出的脫硫性能相符。

(2) 分析了擴散分子與聚合物的相互作用能、擴散分子尺寸、聚合物自由體積分數和聚合物分子鏈運動性對分子擴散的影響，考察了擴散分子在聚合物中微觀擴散機理。結果表明，擴散分子與聚合物相互作用是噻吩類硫化物在脫硫膜材料中優先擴散的決定因素，且擴散分子與體系的相互作用能越低、擴散分子尺寸越小、聚合物自由體積分數越大、聚合物分子鏈運動性越強，越有利于擴散分子在聚合物中的擴散。各因素相互協同共同影響擴散分子在聚合物中的擴散速率。

[1] 孔瑛， 盧福偉， 呂宏凌， 等. 滲透汽化膜法汽油脫硫技術研究進展[J].膜科學與技術， 2011， 31(3)： 162-171.(KONG Ying， LU Fuwei， LU Hongling， et al. The research progress of gasoline desulfurization by pervaporation[J].Membrane Science and Technology， 2011， 31(3)： 162-171.)

[2] 林立剛， 張玉忠. 膜技術在汽油脫硫中應用的新進展[J].石油學報(石油加工)， 2010， 26(3)： 476-485.(LIN Ligang， ZHANG Yuzhong. Recent advances in sulfur removal from gasoline by membrane process[J].Acta Petrolei Sinica (Petroleum Processing Section)， 2010， 26(3)： 476-485.)

[3] 盧福偉， 孔瑛， 呂宏凌， 等. 分子動力學模擬噻吩/正庚烷在PEG膜中溶解-擴散行為[J].高分子材料科學與工程， 2012， 28(3)： 114-117.(LU Fuwei， KONG Ying， Lü Hongling， et al. Molecular dynamics simulation studies on sorption and diffusion properties of thiophene/heptane in peg membranes[J].Polymer Materials Science and Engineering， 2012， 28(3)： 114-117.)

[4] 王俊， 朱宇， 陸小華. 分子模擬在聚合物膜研究中的應用[J].現代化工， 2003， 23(10)： 59-62.(WANG Jun， ZHU Yu， LU Xiaohua. Application of molecular simulation in polymeric membrane science[J].Modern Chemical Industry， 2003， 23(10)： 59-62.)

[5] LI Ben， PAN Fusheng， FANG Zhiping， et al. Molecular dynamics simulation of diffusion behavior of benzene/water in PDMS-Calix[4]arene hybrid pervaporation Membranes[J].Ind Eng Chem Res， 2008， 47(13)： 4440-4447.

[6] LING Changjian， LIANG Xiaoyi， FAN Feichao， et al. Diffusion behavior of the model diesel components in different polymer membranes by molecular dynamic simulation[J].Chemical Engineering Science， 2012， 84： 292-302.

[7] AMANI M， AMJAD-IRANAGH S， GOLZAR K， et al. Study of nanostructure characterizations and gas separation properties of poly(urethane-urea)s membranes by molecular dynamics simulation[J].Journal of membrane Science， 2014， 462： 28-41.

[8] RYCHLEWSKA K， KONIECZNY K， BODZEK M. Pervaporative desulfurization of gasoline separation of thiophene/n-heptane mixture[J].Archives of Environmental Protection， 2015， 41(2)： 3-11.

[9] 林立剛. 滲透汽化膜分離法脫除汽油中有機硫化物的研究[D].青島： 中國石油大學(華東), 2008.

[10] GOLZAR K， AMJAD-IRANAGH S， AMANI M， et al. Molecular simulation study of penetrant gas transport properties into the pure and nanosized silica particles filled polysulfone membranes[J].Journal of Membrane Science， 2014， 451： 117-134.

[11] 盧福偉. 共聚高分子膜材料微相結構與汽油脫硫性能關系研究 [D].青島： 中國石油大學(華東), 2011.

[12] 袁帥， 龍軍， 田輝平， 等. 烴分子尺寸及其與擴散能壘關系的初步研究[J].石油學報(石油加工)， 2011， 27(3)： 376-380.(YUAN Shuai， LONG Jun， TIAN Huiping, et al. Molecular dimensions of hydrocarbons and a primary study on its relationship to diffusion energy barriers in zeolites[J].Acta Petrolei Sinica (Petroleum Processing Section)， 2011， 27(3)： 376-380.)

[13] LIPATOV Y. The iso-free-volume state and glass transitions in amorphous polymers[J].Advances in Polymer Science， 1978， 26： 63-104.

[14] 潘福生. 超薄復合膜的理性設計、結構調控及脫濕性能研究[D].天津： 天津大學, 2009.

[15] PAN Fusheng， MA Jing， CUI Li， et al. Water vapor/propylene sorption and diffusion behavior in PVA-P(AA-AMPS) blend membranes by GCMC and MD simulation[J].Chemical Engineering Science， 2009， 64(24)： 5192-5197.

Diffusion Behavior of the Model Gasoline Components inDifferent Polymer Membranes Studied by Molecular Dynamic Simulation

HOU Yingfei1,2, JIANG Chi1, LIU Min1, LIU Yibin1, LI Peng1, NIU Qingshan1

(1.StateKeyLaboratoryofHeavyOilProcessing,ChinaUniversityofPetroleum,Qingdao266580,China；2.DepartmentofPetroleumandChemicalEngineering,XinjiangVocationalandTechnicalCollegeofApplication,Kuitun833200,China)

To explore the main factor in desulfurization performance of pervaporation membranes, molecular dynamic (MD) simulation was used to investigate the diffusion behavior of model gasoline components (thiophene,n-heptane, cyclohexane, cyclohexene, and toluene) in three different membranes of ethyl cellulose (EC), polyethylene glycol (PEG), and polydimethylsiloxane (PDMS). Trajectories of molecules diffusing were drawn by obtaining the barycentric coordinates of gasoline components in three polymer membranes. Diffusion coefficients of five gasoline components in different membranes were calculated, and the micro diffusion mechanism of diffusion molecule inside polymer membrane was investigated by analyzing the interaction energy between diffusion molecule and polymer membrane, the size of diffusion molecule, fractional free volume (FFV) of polymer membrane and the mobility of polymer chains. It is obviously discovered that the diffusion coefficients of thiophene were considerably greater than other components in three different membranes, which agreed well with the experimental results. Overall, the diffusion coefficients were attributed to the combined effect of the four factors.

molecular dynamic; molecular diffusion; gasoline desulfurization; membrane separation

2015-12-29

國家自然科學基金項目(51173203)、中央高校基本科研業務費專項資金項目(15CX05014A)、山東省自然科學基金項目(ZR2012BL15)資助

侯影飛，男，副教授，博士，從事石油與天然氣加工、膜分離、油田化學品等研究；Tel： 0532-86984702；E-mail： houyf@upc.edu.cn

1001-8719(2016)06-1121-07

TQ028.8

A

10.3969/j.issn.1001-8719.2016.06.006