陳飛龍，張延志，何志霞

基于不同勢函數的正構烷烴物理性質分子動力學模擬

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

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

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

分子動力學模擬；物理性質；勢函數；二元組分

步入全球現代化的21世紀以來，化石燃料具有潛在的能源短缺的危機，燃燒化石燃料(主要是碳氫燃料)而產生的二氧化碳是加快全球變暖的溫室氣體的主要來源之一．面對能源短缺和環(huán)境惡化的雙重危機，節(jié)約能源和清潔低碳化發(fā)展已成為全球共識．除了發(fā)展可再生能源和探索替代能源之外，提高碳氫燃料的燃燒利用效率至關重要．實現這一目標的基礎和前提是對碳氫化合物各項物理性質的定量估計，而對于燃料在一些極端工況如高溫高壓的環(huán)境條件下或者是不同組分混合物熱力學以及輸運性質的實驗數據比較匱乏，且實驗成本較高，因此用理論預測的手段來補充很有必要．

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

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

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

1?模擬方法

1.1?勢函數模型

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

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

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

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

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

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

各聯合原子模型力場的4項勢能參數列于表1．

1.2?熱物性統計方法及計算公式

表1?不同力場的勢參數

Tab.1?Potential parameters for different force fields

1.3?模擬設置

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

圖1?初始構型

表2?不同模擬過程的初始設置

Tab.2?Initialsettings for different simulation processes

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

2?結果與討論

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

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

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

2.2.1?飽和物性

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

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

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

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

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

2.2.2?純液相的自擴散系數

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

圖5?不同力場下的液相自擴散系數

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

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

圖7?相對誤差隨溫度、混合比的變化關系

3?結?論

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

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

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

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

［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)境下氣液界面特性的分子動力學模擬[D]. 大連：大連理工大學能源與動力學院，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．

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

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

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

(責任編輯：隋韶穎)