基于不同勢函數(shù)的正構(gòu)烷烴物理性質(zhì)分子動力學(xué)模擬

陳飛龍，張延志，何志霞

基于不同勢函數(shù)的正構(gòu)烷烴物理性質(zhì)分子動力學(xué)模擬

陳飛龍1，張延志1，何志霞2

(1. 大連理工大學(xué)能源與動力學(xué)院，大連 116024；2. 江蘇大學(xué)能源研究院，鎮(zhèn)江 212013)

采用分子動力學(xué)模擬詳細(xì)比較了4種常用勢函數(shù)對正構(gòu)烷烴物理性質(zhì)預(yù)測精度的影響，提出了一種新的修正力場并對其進(jìn)行廣泛的驗證，驗證主要包括飽和兩相密度、飽和蒸氣壓和自擴(kuò)散系數(shù)等熱物理和輸運性質(zhì)．之后，將改進(jìn)力場應(yīng)用于二元組分燃料熱物性研究中，探究了不同混合比例條件下通過每種組分加權(quán)平均的方法估算混合物飽和熱物性的準(zhǔn)確度．研究發(fā)現(xiàn)，改進(jìn)力場對正庚烷和正十二烷物理性質(zhì)的預(yù)測與實驗值較為吻合；在溫度較低、正庚烷摩爾分?jǐn)?shù)較高或者溫度較高、正庚烷摩爾分?jǐn)?shù)較低的工況下，通過加權(quán)平均得到的飽和物性值較為準(zhǔn)確．

分子動力學(xué)模擬；物理性質(zhì)；勢函數(shù)；二元組分

步入全球現(xiàn)代化的21世紀(jì)以來，化石燃料具有潛在的能源短缺的危機(jī)，燃燒化石燃料(主要是碳?xì)淙剂?而產(chǎn)生的二氧化碳是加快全球變暖的溫室氣體的主要來源之一．面對能源短缺和環(huán)境惡化的雙重危機(jī)，節(jié)約能源和清潔低碳化發(fā)展已成為全球共識．除了發(fā)展可再生能源和探索替代能源之外，提高碳?xì)淙剂系娜紵眯手陵P(guān)重要．實現(xiàn)這一目標(biāo)的基礎(chǔ)和前提是對碳?xì)浠衔锔黜椢锢硇再|(zhì)的定量估計，而對于燃料在一些極端工況如高溫高壓的環(huán)境條件下或者是不同組分混合物熱力學(xué)以及輸運性質(zhì)的實驗數(shù)據(jù)比較匱乏，且實驗成本較高，因此用理論預(yù)測的手段來補(bǔ)充很有必要．

狀態(tài)方程是預(yù)測流體熱物性的有效手段，但是先決條件是需要進(jìn)行各種假設(shè)、需要提供經(jīng)驗參數(shù)以及混合規(guī)則等，存在對于不同環(huán)境不同組分的適用性問題，同時計算效率和精度等問題也同樣需要考慮．為了避免傳統(tǒng)數(shù)值模擬和實驗所遇到的問題，分子動力學(xué)(molecular dynamics，MD)模擬技術(shù)近年來被越來越地多應(yīng)用于流體性質(zhì)的研究．MD模擬的優(yōu)勢是不需要引入對所研究物理過程的任何假設(shè)，只要有能正確描述物質(zhì)粒子間作用力的勢函數(shù)，適用于對各項熱力學(xué)以及輸運性質(zhì)的預(yù)測．

MD計算中，勢函數(shù)作為描述粒子間作用力近似數(shù)學(xué)表達(dá)方法的選取極為關(guān)鍵．碳?xì)浠衔锏某Ｓ梅肿幽Ｐ椭饕譃槿幽Ｐ秃吐?lián)合原子模型，前者視每個原子為一個位點，后者將每個甲基CH3和亞甲基CH2基團(tuán)分別視為一個聯(lián)合原子來簡化計算．最早出現(xiàn)的有Jorgensen等[1-4]提出的OPLS-AA和OPLS-UA力場，廣泛應(yīng)用于蛋白質(zhì)、烴類和醇類等復(fù)雜的大分子物質(zhì)．之后Siepmann等[5-6]提出的SKS力場，更適用于模擬鏈長較短的烷烴．除此之外，又相繼出現(xiàn)了Martin等[7]通過擬合烷烴的臨界溫度和飽和液相密度的實驗數(shù)據(jù)提出的TraPPE-UA力場以及Nath等[8]提出的對短長鏈烷烴都適用的NERD力場．鄧?yán)赱9]比較了OPLS-AA力場以及一種聯(lián)合原子力場對正庚烷液滴的蒸發(fā)特性的影響，結(jié)果表明，在采用截斷半徑法的計算下所選用的聯(lián)合原子力場比OPLS-AA力場對正庚烷粒子間作用力的描述更優(yōu)，而且全原子模型計算成本較高，尤其對于長鏈正構(gòu)烷烴．然而，到目前為止，詳細(xì)比較不同常用聯(lián)合原子勢函數(shù)對正構(gòu)烷烴物理性質(zhì)預(yù)測精度的研究比較缺乏．

本文選取了內(nèi)燃機(jī)燃料中典型組分正庚烷和正十二烷為研究對象，比較了不同勢函數(shù)對正構(gòu)烷烴物理性質(zhì)(包括飽和及輸運物性)的影響，并提出了一種新的修正力場并對其進(jìn)行了烷烴飽和兩相密度、飽和蒸氣壓和自擴(kuò)散系數(shù)等熱物理性質(zhì)的驗證．另外，將改進(jìn)力場應(yīng)用于雙組分燃料(正庚烷和正十二烷)熱物性的研究中，主要探究了不同混合比例條件下通過將每種組分飽和物性按摩爾比例加權(quán)平均的方法來估算混合物飽和物性的準(zhǔn)確度．

1?模擬方法

1.1?勢函數(shù)模型

分子動力學(xué)模擬中，勢函數(shù)或者力場模型是作為粒子間作用力的近似描述，力場的選取對于整個模擬的過程可謂重中之重，因此在進(jìn)行模擬前需要選擇準(zhǔn)確合適的力場．每種物質(zhì)有各自不同適合于本身的力場參數(shù)，同一種物質(zhì)也有適合于表述各種性質(zhì)而對應(yīng)的不同的力場參數(shù)．為了節(jié)約計算成本，本文采用了常用的4種適用于碳?xì)浠衔锏穆?lián)合原子力場模型OPLS-UA、SKS、TraPPE-UA以及NERD．此外，基于SKS非鍵勢能、鍵伸縮勢能以及鍵角彎曲勢能和另一種五參數(shù)約束的二面角扭曲勢能[10]，本文提出一種新的修正力場(命名為SKS-New)．

聯(lián)合原子力場的非鍵勢能可以用截斷的Lennard-jones12-6勢來描述：

在最初的SKS和TraPPE-UA的力場模型中，將所有鍵視為固定長度的剛性鍵，然而LAMMPS平臺不提供正構(gòu)烷烴等大分子的鍵長約束，因此本研究中統(tǒng)一使用諧波勢來代替剛性鍵拉伸描述鍵伸縮勢能：

同樣用一種諧波勢來描述鍵角彎曲勢能：

用OPLS二面角扭曲勢[4]描述前4種力場的二面角扭曲勢能：

改進(jìn)力場的二面角扭曲勢能則由余弦多項式表示的多諧波勢[10]來描述：

各聯(lián)合原子模型力場的4項勢能參數(shù)列于表1．

1.2?熱物性統(tǒng)計方法及計算公式

表1?不同力場的勢參數(shù)

Tab.1?Potential parameters for different force fields

1.3?模擬設(shè)置

本文比較了5種力場對從低溫到接近臨界溫度(正庚烷臨界溫度為540K、正十二烷臨界溫度為658K)下氣液平衡模擬純組分正構(gòu)烷烴飽和物性和純液相烷烴的自擴(kuò)散特性的影響以及研究了不同摩爾比(正庚烷：正十二烷)的雙組分烷烴在正庚烷臨界溫度以下的飽和物性．系統(tǒng)的初始構(gòu)型如圖1所示，初始構(gòu)型氣液氣構(gòu)型用于飽和物性的計算；純液相構(gòu)型用于自擴(kuò)散系數(shù)的計算．因考慮實際計算效率和對結(jié)果的影響，不同的模擬過程選擇不同的適合的分子數(shù)以及模擬盒尺寸，詳細(xì)模擬細(xì)節(jié)在表2中總結(jié).

圖1?初始構(gòu)型

表2?不同模擬過程的初始設(shè)置

Tab.2?Initialsettings for different simulation processes

氣液氣構(gòu)型中模擬盒子為長方體，初始時刻液相分子位于盒子中心，兩端為氣相區(qū)．為消除模擬體系的規(guī)模限制帶來的邊界效應(yīng)，在模擬盒3個方向上均采用周期性邊界條件，整個模擬過程均使用正則(NVT)系綜對系統(tǒng)溫度進(jìn)行恒定．采用velocity-verlet算法對原子運動方程進(jìn)行積分，原子初始速度按高斯分布隨機(jī)分配．取時間步長為2×10-6ns，模擬總時長為3ns，其中前2ns為弛豫時間，后1ns為統(tǒng)計時間．當(dāng)整個系統(tǒng)的溫度、壓力、總能和勢能4個指標(biāo)僅在固定值附近波動時，則認(rèn)定達(dá)到平衡狀態(tài)．純液相模型為立方模擬盒，選擇合適的模擬盒尺寸，使初始構(gòu)型的液體密度在不同溫度下接近真實值，本次模擬中初始密度取自NIST數(shù)據(jù)庫[13]．先用NPT系綜運行1ns使整個系統(tǒng)壓力穩(wěn)定在常壓附近，然后再使用NVT系綜運行2ns，同樣取最后1ns作為統(tǒng)計時間．

2?結(jié)果與討論

2.1?純正庚烷的氣液平衡MD模擬

氣相與液相的相交處為氣液界面，體系平衡狀態(tài)下即可獲得穩(wěn)定界面．MD模擬很好地再現(xiàn)了氣液兩相平衡特征，各力場下氣液相密度分布曲線如圖2所示，在整個體系達(dá)到平衡之后，根據(jù)這些分布，可以識別出密度較高的液相區(qū)域、密度較低的氣相區(qū)域和氣液界面．5種聯(lián)合原子力場在低于473K以下時都有較好的模擬表現(xiàn)，能夠完整地模擬正庚烷從293～473K的飽和過程，有明顯的密度梯度存在，溫度的升高導(dǎo)致液相的密度下降，而導(dǎo)致氣相的密度增加；而在高于473K時，因接近于臨界溫度(540K)，密度梯度縮小而界面區(qū)慢慢消失，5種力場中只有OPLS-UA和改進(jìn)力場SKS-New能完整模擬所有給定溫度下的氣液平衡過程，在高于473K時也能有明顯的氣液密度梯度，SKS力場在533K時密度梯度消失，而TraPPE-UA和NERD力場在513K和533K時模擬效果較差，已全部變?yōu)闅庀啵?/p>

2.2?對純組分正庚烷和正十二烷的熱物性驗證

2.2.1?飽和物性

由于界面輪廓隨時間波動，氣液界面的位置在不同的模擬中有所不同．盡管有些波動，氣液界面的統(tǒng)計位置依然穩(wěn)定．平均統(tǒng)計了5種力場、兩種純組分在各溫度下兩邊氣相區(qū)密度和中間液相區(qū)密度，來與NIST數(shù)據(jù)庫[13]的兩相密度進(jìn)行對比．所得氣液共存曲線如圖3所示．正庚烷氣液共存曲線5種力場的模擬結(jié)果總體上大都與NIST值相符，其中TraPPE-UA和NERD力場在413K以下的低溫區(qū)氣液兩相區(qū)吻合較好，在高于413K時誤差增大，氣相區(qū)密度與NIST值相比偏大而液相區(qū)密度偏小，由于二者在513K以上之后就全部變?yōu)闅庀?，因此無513K以上的氣液共存密度，SKS力場在533K時亦是如此．OPLS-UA力場對所有溫度下的液相密度模擬結(jié)果都略微偏大，而SKS及提出的改進(jìn)力場的模擬結(jié)果與NIST值最為吻合，但由于新力場在533K時還能保持良好的氣液兩相平衡，說明改進(jìn)力場對臨界溫度的預(yù)測效果最好．因為選取精確的勢函數(shù)是表征模擬準(zhǔn)確性的前提，力場描述的粒子間吸引力越小，在給定溫度下從中心液核弛豫到氣相區(qū)的分子越多，使得液相區(qū)密度偏小而氣相區(qū)密度偏大，從模擬結(jié)果來看，OPLS-UA力場表述的正構(gòu)烷烴分子吸引力偏大，TraPPE-UA和NERD力場則偏小，而SKS和改進(jìn)力場相對而言最為準(zhǔn)確．

圖3?不同力場下的氣液共存曲線

不同力場下正十二烷從450～650K的模擬結(jié)果與正庚烷大致趨勢相同，同樣是SKS和改進(jìn)力場的模擬結(jié)果與NIST值較為符合，同樣SKS、TraPPE-UA和NERD這3種力場在溫度升高時氣相密度會偏大、液相密度會偏小，而OPLS-UA力場則相反.從氣液共存密度的模擬結(jié)果來看，改進(jìn)力場SKS-New對正十二烷的描述優(yōu)于其他4種力場，改進(jìn)力場在每個溫度下的模擬結(jié)果能相對比較準(zhǔn)確地模擬出氣液相密度，且在650K即接近臨界溫度(658K)時也能有明顯氣液界面，在5種力場中甚至比SKS力場對高溫度下的描述更為準(zhǔn)確，預(yù)測臨界點的效果更好．這可能是由于改進(jìn)力場勢能參數(shù)所用的二面角扭轉(zhuǎn)勢能參數(shù)比原SKS力場更為準(zhǔn)確，同時在鏈更長的烷烴里所包含的二面角更多，因此對二面角扭曲勢能更為準(zhǔn)確的描述會使得模擬結(jié)果更為精確．

燃料的飽和蒸氣壓對于發(fā)動機(jī)內(nèi)混合氣形成以及后續(xù)燃燒過程有著重要影響，因此為驗證此改進(jìn)力場參數(shù)是否能準(zhǔn)確表述烷烴其他物性，用氣液氣模型計算了烷烴的不同溫度下的飽和蒸氣壓．將氣液平衡狀態(tài)下兩端飽和氣相區(qū)的壓力統(tǒng)計平均值即可得到飽和蒸氣壓，MD模擬正庚烷和正十二烷的計算結(jié)果與實驗值[14-15]和NIST數(shù)據(jù)庫[13]的對比如圖4所示．兩種純組分不同力場下的飽和蒸氣壓都隨溫度升高而增大，模擬結(jié)果與前面兩相密度的效果大致相同，5種力場中還是SKS和改進(jìn)力場SKS-New的模擬結(jié)果相對較為準(zhǔn)確，其他3種力場與實驗值或數(shù)據(jù)庫數(shù)據(jù)偏差較大，尤其在對正十二烷的模擬中更加明顯．SKS和改進(jìn)力場模擬的飽和蒸氣壓都隨溫度升高而誤差增大，而對于正十二烷改進(jìn)力場的飽和蒸氣壓在高溫區(qū)也能較好地吻合且比SKS力場表現(xiàn)更好，這也印證了上面對于改進(jìn)力場勢能參數(shù)更適用于表述中長鏈烷烴的推斷．

圖4?不同力場下的飽和蒸氣壓

2.2.2?純液相的自擴(kuò)散系數(shù)

除了對氣液平衡體系下正構(gòu)烷烴進(jìn)行不同力場的MD模擬驗證外，還通過建立純組分純液相的立方初始構(gòu)型來模擬計算液相的自擴(kuò)散系數(shù)，用不同力場模擬的正庚烷和正十二烷結(jié)果與實驗值[16-17]對比分別如圖5(a)和(b)所示．自擴(kuò)散系數(shù)的總體趨勢隨著溫度升高而增大，與前面氣液模型計算的兩相密度和飽和蒸氣壓效果不同，在５種力場中最為相符的是OPLS-UA力場，其次為SKS-New改進(jìn)力場和SKS力場，而另外兩種力場還是與實驗值相比偏大．結(jié)合上一小節(jié)的研究來看，這可能是因為OPLS-UA力場描述的粒子間吸引力在5種力場中最大，因此模擬過程中分子在自由運動中的均方位移最小，通過均方位移計算的自擴(kuò)散系數(shù)也就最小，而其他力場對應(yīng)的描述的粒子間吸引力越小，則計算的自擴(kuò)散系數(shù)越大．

圖5?不同力場下的液相自擴(kuò)散系數(shù)

2.3?正庚烷/正十二烷二元組分的飽和物性

式中：為純組分的摩爾分?jǐn)?shù)，為純組分物性的模擬值．結(jié)果如圖6所示，混合物物性并不是單組分物性的簡單加權(quán)平均，由于正十二烷分子對正庚烷分子的吸附作用，本該在接近臨界溫度時全部變?yōu)闅庀嗟恼榉肿颖阮A(yù)計更少地從液核弛豫到氣相區(qū)，因此液相密度的加權(quán)值始終比實際模擬值偏小，相應(yīng)地，氣相密度和飽和蒸氣壓始終偏大，且總體趨勢似乎是隨著溫度增大而偏差越多．

圖7?相對誤差隨溫度、混合比的變化關(guān)系

3?結(jié)?論

采用MD模擬比較了不同勢函數(shù)對不同溫度下純組分正庚烷和正十二烷燃料物理性質(zhì)的影響，并提出了一種新的修正力場并對其進(jìn)行了烷烴飽和兩相密度、飽和蒸氣壓和自擴(kuò)散系數(shù)等熱物性的驗證．還將改進(jìn)力場應(yīng)用于二元組分烷烴的模擬中，主要研究了不同摩爾比(1∶2、1∶1、2∶1)的雙組分烷烴在正庚烷臨界溫度下的飽和過程．本文的研究結(jié)果如下.

(1)建立了單組分烷烴氣液平衡的MD模型．比較了不同勢函數(shù)對于預(yù)測兩種正構(gòu)烷烴兩相飽和物性的精度差異，分析了溫度對于飽和狀態(tài)下兩相密度以及飽和蒸氣壓的影響．隨溫度升高，液相密度減小，氣相密度與飽和蒸氣壓增大，5種力場皆能良好地模擬烷烴在近臨界溫度下的飽和過程．

(2)提出了一種改進(jìn)的適用于正構(gòu)烷烴的修正力場，比較了與其他4種常用力場對預(yù)測飽和兩相密度、飽和蒸氣壓以及自擴(kuò)散系數(shù)的精度．其中SKS力場對于正庚烷的飽和性質(zhì)模擬最為準(zhǔn)確，但對于臨界點的預(yù)測不如改進(jìn)力場SKS-New；改進(jìn)力場對于正十二烷的飽和性質(zhì)模擬相較于SKS力場要更優(yōu)，且對兩種正構(gòu)烷烴的臨界溫度的預(yù)測更準(zhǔn)確；OPLS-UA力場對于自擴(kuò)散系數(shù)的模擬最為準(zhǔn)確，其次為改進(jìn)力場．綜合來看，改進(jìn)力場對正構(gòu)烷烴的粒子間作用力描述最優(yōu)且更適用于中長鏈烷烴．

(3)研究了雙組分烷烴飽和過程．探究了不同混合比和溫度條件下通過將純組分飽和物性按不同摩爾比加權(quán)平均的方法估算混合物飽和物性的準(zhǔn)確度．在估算飽和物性時，選擇溫度較低、正庚烷比重較高或者是溫度較高、正庚烷比重較低的工況下加權(quán)平均更為準(zhǔn)確．

［1］ Jorgensen W L，Tirado-Rives J. The OPLS[optimized potentials for liquid simulations] potential functions for proteins，energy minimizations for crystals of cyclic peptides and crambin[J].，1988，110(6)：1657-1666.

［2］ Jorgensen W L，Maxwell D S，Tirado-Rives J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids[J].，1996，118(45)：11225-11236.

［3］ Jorgensen W L. Optimized intermolecular potential functions for liquid alcohols [J].，1986，90(7)：1276-1284.

［4］ Jorgensen W L，Madura J D，Swenson C J. Optimized intermolecular potential functions for liquid hydrocarbons[J].，1984，106(22)：6638-6646.

［5］ Siepmann J I，Karaborni S，Smit B. Vapor-liquid equilibria of model alkanes[J].，1993，115(14)：6454-6455.

［6］ Siepmann J I，Karaborni S，Smit B. Simulating the critical behaviour of complex fluids[J].，1993，365：330-332.

［7］ Martin M G，Siepmann J I. Transferable potentials for phase equilibria(1)：United-atom description of n-alkanes[J].，1998，102(14)：2569-2577.

［8］ Nath S K，Escobedo F A，de Pablo J J. On the simulation of vapor-liquid equilibria for alkanes[J].，1998，108(23)：9905-9911.

［9］ 鄧?磊. 亞/超臨界環(huán)境下氣液界面特性的分子動力學(xué)模擬[D]. 大連：大連理工大學(xué)能源與動力學(xué)院，2015.

Deng Lei. Molecular Dynamics Simulation of Liquid-Gas Interface in Sub/Supercritical Surroundings[D]. Dalian：School of Energy and Power Engineering，Dalian University of Technology，2015(in Chinese).

［10］ Kalyanasundaram V，Spearot D E，Malshe A P. Molecular dynamics simulation of nanoconfinement induced organization of n-decane[J].，2009，25(13)：7553-7560.

［11］ Mo G，Qiao L. A molecular dynamics investigation of n-alkanes vaporizing into nitrogen：Transition from subcritical to supercritical[J].，2017，176：60-71.

［12］ Lee S-H. Pressure analyses at the planar surface of liquid-vapor argon by a test-area molecular dynamics simulation[J].，2012，33(9)：3039-3042.

［13］ NIST. NIST Chemistry WebBook[EB/OL]. https://?webbook.nist.gov/chemistry/，2022-03-16.

［14］ Weber L A. Vapor pressure of heptane from the triple point to the critical point[J].，2000，45(2)：173-176.

［15］ Morgan D L，Kobayashi R. Direct vapor pressure measurements of ten n-alkanes in the C10-C28 range[J].，1994，97：211-242.

［16］ Fishman E. Self-diffusion in liquid n-pentane and n-heptane[J].，1955，59(5)：469-472.

［17］ Harris K R，Ganbold B，Price W S. Viscous calibration liquids for self-diffusion measurements[J].，2015，60(12)：3506-3517.

Molecular Dynamics Simulation of Physical Properties of n-Alkane Based on Different Potential Models

Chen Feilong1，Zhang Yanzhi1，He Zhixia2

(1. School of Energy and Power，Dalian University of Technology，Dalian 116024，China；2. Institute for Energy Research，Jiangsu University，Zhenjiang 212013，China)

In this paper，molecular dynamics simulation was used to compare in detail the effect of four commonly used potential models on the prediction accuracy of the physical properties of n-alkane，and a new modified model was proposed and extensively validated，which mainly includes the thermophysical and transport properties such as saturated two-phase density，saturated vapor pressure and self-diffusion coefficient. After that，the modified potential model was applied to the study of the thermophysical properties of binary-component fuels，and the accuracy of estimating the saturated physical properties of the mixture by the weighted average method of each component at different mixing ratios was explored. It was found that the predictions from the modified model show good agreement with the measurements for n-heptane and n-dodecane. When the temperature is lower and the mole fraction of n-heptane is higher，or when the temperature is higher and the mole fraction ofnheptane is lower，the saturated physical property value obtained by the weighted average is more accurate.

molecular dynamics simulation；physical property；potential model；binary-component

10.11715/rskxjs.R202305023

TK401

A

1006-8740(2023)04-0451-09

2022-03-16．

國家自然科學(xué)基金資助項目(52106154)；江蘇省自然科學(xué)基金資助項目(BK20190855).

陳飛龍(1998—??)，男，碩士研究生，fl_chen98@163.com.

張延志，男，博士，副教授，zhangyanxzhi@dlut.edu.cn.

(責(zé)任編輯：隋韶穎)