基于分子動(dòng)力學(xué)模擬的聚合物類復(fù)合蓄冷劑相變機(jī)理研究

陳有良 盛 偉,2 王蕊蕊 李雪麗 郝曉茹 鄭海坤

（1 河南理工大學(xué)機(jī)械與動(dòng)力工程學(xué)院 焦作 454000；2 哈密豫新能源產(chǎn)業(yè)研究院 哈密 839000）

在碳達(dá)峰、碳中和的大背景下，相變材料在相變時(shí)通過(guò)潛熱的吸收、儲(chǔ)存和釋放，可以實(shí)現(xiàn)對(duì)現(xiàn)有能源的合理分配，并優(yōu)化和提高能源利用率[1-2]。蓄冷技術(shù)[3-6]是將含有蓄冷劑的裝置在低溫環(huán)境下進(jìn)行預(yù)冷，通過(guò)相變吸收并儲(chǔ)存大量冷量，再根據(jù)人們的需求在溫度較高的環(huán)境中能持續(xù)性的放出冷量，使周圍環(huán)境在一定時(shí)間內(nèi)保持恒定的低溫狀態(tài)，用于滿足人們對(duì)食品保鮮、冷藏或其他需求的技術(shù)。近年來(lái)，蓄冷技術(shù)在冷凍加工、低溫貯藏及長(zhǎng)途運(yùn)輸?shù)确矫鎽?yīng)用廣泛[7-8]，為人們的生活提供了極大便利。

相變蓄冷材料按照組成成分可以分為有機(jī)、無(wú)機(jī)以及復(fù)合相變蓄冷材料[9-10]。有機(jī)蓄冷材料具有較高的化學(xué)穩(wěn)定性,可重復(fù)循環(huán)使用，無(wú)腐蝕性及無(wú)過(guò)冷[11-12]等優(yōu)點(diǎn)。周孫希等[13]采用低共熔法研制了一種由十四烷和正辛酸按一定比例混合而成的相變溫度為0 ℃的二元有機(jī)相變蓄冷材料，但其主要缺點(diǎn)是導(dǎo)熱系數(shù)低、相變潛熱小，這極大限制了有機(jī)蓄冷材料的應(yīng)用。無(wú)機(jī)相變蓄冷材料相變潛熱大，導(dǎo)熱系數(shù)高，主要為化合物、水合鹽及金屬合金[14]。B. Han等[15]利用DSC測(cè)量了不同濃度的NaCl 水溶液的相變潛熱，發(fā)現(xiàn)在共晶濃度下的溶液效果最好，并測(cè)得共晶溫度為-22 ℃。但無(wú)機(jī)相變蓄冷材料也有過(guò)冷度大、相分離、腐蝕性較大等致命缺點(diǎn)，限制了在生活中的應(yīng)用。近年來(lái)，復(fù)合相變材料可以有效克服有機(jī)或無(wú)機(jī)蓄冷材料的缺點(diǎn)，而且復(fù)合相變材料[16]儲(chǔ)能密度大，具有極強(qiáng)的可研究性和良好的發(fā)展前景，得到廣大研究者們的關(guān)注。李靖等[17]以丙三醇和氯化鈉作為蓄冷材料研制出一種相變溫度為-31.5 ℃、相變潛熱為175.3 J/g的無(wú)毒低溫相變材料。Yang Ying等[18]研制了一種乙二醇、氯化銨和水以3∶2∶15的比例混合的相變溫度為-23 ℃、相變潛熱為175.9 kJ/kg的復(fù)合型相變蓄冷材料。隨著復(fù)合相變蓄冷劑的發(fā)展，高分子聚合物材料在物理性質(zhì)、化學(xué)性質(zhì)、穩(wěn)定性方面與無(wú)機(jī)、有機(jī)類相比，具有無(wú)毒、無(wú)腐蝕、比熱和融解熱較大、相變冰點(diǎn)可調(diào)等優(yōu)點(diǎn)，可以廣泛應(yīng)用于更多的復(fù)雜環(huán)境。但由于相變材料在相變前后受溫度、壓力等因素影響較大，目前對(duì)于含聚合物類復(fù)合相變蓄冷劑尤其是微觀機(jī)理方面的研究鮮有報(bào)導(dǎo)。

目前已有大量研究通過(guò)差示掃描量熱法（differential scanning calorimetry,DSC）等熱分析測(cè)試手段來(lái)測(cè)定相變溫度、相變潛熱等相變蓄冷材料的宏觀熱物性[19-20]。根據(jù)相變熱力學(xué)理論，當(dāng)物質(zhì)由一個(gè)相態(tài)轉(zhuǎn)變到另一相態(tài)時(shí)，其關(guān)聯(lián)函數(shù)（密度、擴(kuò)散系數(shù)等）均會(huì)在相變點(diǎn)出現(xiàn)轉(zhuǎn)折或突躍，該溫度點(diǎn)即為相變溫度。依據(jù)此原理，分子動(dòng)力學(xué)模擬（molecular dynamics simulation）已經(jīng)廣泛應(yīng)用于研究材料的玻璃化轉(zhuǎn)變及相變溫度預(yù)測(cè)。目前采用分子動(dòng)力學(xué)方法研究和預(yù)測(cè)相變蓄冷材料的報(bào)道較少。為進(jìn)一步探究聚合物復(fù)合相變蓄冷劑的微觀機(jī)理，本文將采用分子動(dòng)力學(xué)方法，針對(duì)低溫冷鏈物流的進(jìn)一步應(yīng)用需要，保證較低的相變溫度和較高的潛熱值，以高潛熱值的冰水相變?yōu)橹?，采用丙三醇、氯化銨作為成核劑和降溫劑，并在此基礎(chǔ)上添加質(zhì)量分?jǐn)?shù)為0%、3%、5%的聚合物PEG，對(duì)其相變過(guò)程進(jìn)行模擬，從成鍵與構(gòu)象、鍵長(zhǎng)與鍵角、徑向分布函數(shù)、擴(kuò)散系數(shù)和密度等參數(shù)進(jìn)行分析并預(yù)測(cè)其相變溫度，為開發(fā)應(yīng)用于低溫冷鏈物流領(lǐng)域的中低溫復(fù)合相變蓄冷劑材料提供參考。

1 模型建立與參數(shù)設(shè)定

1.1 力場(chǎng)選擇與模擬參數(shù)的設(shè)定

在所有的力場(chǎng)中，COMPASSII力場(chǎng)（凝聚相，能準(zhǔn)確模擬凝聚態(tài)體系的結(jié)構(gòu)和性質(zhì)，用于粒子模擬研究的最優(yōu)化分子勢(shì)）可將無(wú)機(jī)分子體系和有機(jī)分子體系結(jié)合在一起，能夠模擬出體系中各個(gè)獨(dú)立分子的結(jié)構(gòu)、熱力學(xué)性質(zhì)和振動(dòng)頻率等。COMPASSII力場(chǎng)用于液相分子相互作用的函數(shù)形式可定義為[21]：

（1）

模擬主要采用 Materials Studio （MS） 軟件包中Amorphous、Forcite以及Discover模塊中的 Dynamics來(lái)完成[22-23]。體系溫度的控制采用Andersen方法，壓力的控制采用Berendsen方法。設(shè)置模型體系分別在223.15～293.15 K之間進(jìn)行階段性降溫模擬，設(shè)定整個(gè)模擬過(guò)程處于一個(gè)標(biāo)準(zhǔn)大氣壓和 NPT 系綜[24]。采用Velocity-Verlet算法計(jì)算分析求解體系的運(yùn)動(dòng)方程，時(shí)間步長(zhǎng)設(shè)置為1 fs，取樣時(shí)間間隔為1 ps。為使模擬體系達(dá)到充分弛豫，在正式模擬之前，首先對(duì)體系進(jìn)行600 ps的預(yù)平衡模擬計(jì)算。待體系充分弛豫后，進(jìn)行模擬時(shí)長(zhǎng)為2 ns的模擬，前1 ns用于體系的平衡計(jì)算，后1 ns的模擬結(jié)果則用于統(tǒng)計(jì)分析。

1.2 平衡狀態(tài)的判定

在動(dòng)力學(xué)模擬過(guò)程中，隨著模擬時(shí)間步長(zhǎng)的增加，體系的動(dòng)能、勢(shì)能、總能量、非鍵能等所呈現(xiàn)的漲落現(xiàn)象來(lái)判定模擬體系是否達(dá)到平衡狀態(tài)。圖1（a）所示為T=263.15 K時(shí)3%PEG-10%丙三醇-20%氯化銨-水溶液體系（本文均指質(zhì)量分?jǐn)?shù)）的總能量、勢(shì)能、動(dòng)能、非鍵能等隨模擬時(shí)長(zhǎng)的變化，可知該體系的各參數(shù)在一定的平衡位置呈上下均勻波動(dòng)，表明此時(shí)模型體系已經(jīng)達(dá)到平衡狀態(tài)。圖1（b）所示為該模型在T=263.15 K時(shí)該體系在時(shí)間步長(zhǎng)為0、200、400、600 ps體系的模擬結(jié)果，可知，在進(jìn)行到600 ps之后溶液已經(jīng)達(dá)到充分弛豫的狀態(tài)。

圖1 平衡狀態(tài)時(shí)溶液體系的能量、溫度隨時(shí)間的變化Fig.1 Variation of energy and temperature of the solution system in the equilibrium state with time

圖2 強(qiáng)氫鍵作用O—H…O和C—O…H（紅色）與弱氫鍵作用C—H…O和N—H…O（藍(lán)色）Fig.2 Strong hydrogen bonding O—H…O and C—O…H （red line） and weak hydrogen bonding C—H…O and N—H…O （blue line）

2 結(jié)果與討論

2.1 成鍵與構(gòu)象分析

在10%丙三醇-20%氯化銨-水溶液中，主要存在強(qiáng)氫鍵和弱氫鍵兩種作用力，如圖2所示。研究發(fā)現(xiàn)[25]，在含羥基醇類水溶液中強(qiáng)氫鍵作用O—H…O和C—O…H對(duì)溶液的性質(zhì)起著決定性作用，而C—H…O和N—H…O氫鍵作用相對(duì)較弱，對(duì)相變蓄冷材料的性質(zhì)和狀態(tài)影響較小。

圖3所示為在PEG-丙三醇-氯化銨-水溶液體系中，添加劑分子與水分子形成的氫鍵示意圖。在相變過(guò)程中，伴隨著大量的潛熱被釋放出來(lái)，水分子與水分子之間通過(guò)氫鍵締結(jié)成冰晶微核，進(jìn)而會(huì)有更多的水分子黏附在上面形成更大的冰晶核，隨著溫度的進(jìn)一步降低最終將導(dǎo)致產(chǎn)生相變。而溶液中的添加劑分子通過(guò)與水分子之間的氫鍵作用，這些氫鍵在體系中呈立體網(wǎng)狀結(jié)構(gòu)，對(duì)水分子運(yùn)動(dòng)有一定的約束作用，消耗溶液體系中的自由水分子，使其轉(zhuǎn)變?yōu)榻Y(jié)合水即“束縛水”，溶液中的自由水含量減少，水分子之間的碰撞幾率變小，從而達(dá)到抑制冰核的形成，降低相變溫度的目的。

圖3 PEG-丙三醇-氯化銨-水溶液體系中添加劑與水分子成鍵示意圖Fig.3 Schematic diagram of the bonding between additives and water molecules in the PEG-glycerol-ammonium chloride-water system

圖4 不同溫度下3%PEG-10%丙三醇-20%氯化銨-水溶液構(gòu)象圖Fig.4 3% PEG-10% glycerol-20% ammonium chloride-water conformation diagram at different temperatures

圖4所示為在不同溫度下3%PEG-10%丙三醇-20%氯化銨-水溶液體系模型。從整個(gè)動(dòng)力學(xué)過(guò)程來(lái)看，溫度越低，分子動(dòng)能越小，分子運(yùn)動(dòng)越平緩。在243.15～273.15 K溫度區(qū)間內(nèi)，溫度對(duì)于丙三醇和水分子的影響較顯著。丙三醇分子與水分子之間通過(guò)氫鍵結(jié)合傾向于形成小團(tuán)簇，分子整體運(yùn)動(dòng)較為平緩。溫度對(duì)于水分子之間的結(jié)合力影響較大。當(dāng)T=273.15 K時(shí)，分子運(yùn)動(dòng)劇烈，分子的擴(kuò)散速度較大。隨著溫度的降低，水分子的運(yùn)動(dòng)速度減小。當(dāng)T=243.15 K時(shí)，水分子的運(yùn)動(dòng)速度相較于273.15 K時(shí)大幅降低。

2.2 鍵長(zhǎng)與鍵角隨溫度的變化

圖5所示為丙三醇分子與水分子構(gòu)象示意圖。隨著溫度的降低，物質(zhì)分子的微觀形態(tài)也會(huì)發(fā)生變化，如鍵長(zhǎng)與鍵角也會(huì)隨著相變過(guò)程溫度的變化而發(fā)生變化，通過(guò)分析相變過(guò)程中分子的鍵長(zhǎng)與鍵角的變化來(lái)進(jìn)一步探究相變過(guò)程機(jī)理。

圖5 體系中丙三醇分子和水分子鍵長(zhǎng)鍵角構(gòu)象圖Fig.5 Bond angle conformation diagram of the bond length of glycerol molecule and water molecule in the system

圖6 體系中H—O—H鍵角和O—H鍵長(zhǎng)隨溫度的變化Fig.6 Variation of H—O—H bond angle and O—H bond length with temperature in the system

圖6 分別統(tǒng)計(jì)了∠H—O—H類型的鍵角和O—H型羥基鍵長(zhǎng)隨溫度的變化。由圖6（a）可知，溫度對(duì)于∠H—O—H型鍵角的影響較小，其概率峰值會(huì)隨著溫度的降低而呈增加趨勢(shì)，但增加量較小。而對(duì)于O—H型鍵長(zhǎng)，由圖6（b）可知，隨著溫度的降低，鍵長(zhǎng)變化較為顯著，其概率峰值隨著溫度的降低而增大，鍵長(zhǎng)向變小的方向變化，且鍵長(zhǎng)變短的概率隨溫度的降低而增加。上述結(jié)果表明，在蓄冷材料相變過(guò)程中，溫度對(duì)鍵角的影響較小，對(duì)鍵長(zhǎng)的影響較大，羥基鍵長(zhǎng)變短的概率增加。

2.3 徑向分布函數(shù)（RDF）

徑向分布函數(shù)（radial distribution function, RDF）是整個(gè)模型體系的區(qū)域密度和平均密度的比值,是描述以某一中心原子為中心周圍出現(xiàn)另一種原子的概率，可用來(lái)表征物質(zhì)的有序性。gαβ（r）數(shù)學(xué)表達(dá)式為：

（2）

當(dāng)r較小時(shí)，gαβ（r）可用來(lái)表征原子的堆積情況和各個(gè)原子對(duì)的成鍵距離，峰值高且尖，則有序性強(qiáng)，原子之間聯(lián)系較緊密。在溶液體系中，存在于分子間的非鍵合力通?？煞譃椋簹滏I作用力和范德華力，前者的作用范圍為 2.6～3.1 ?，后者的作用范圍 3.1～5.0 ?。

圖7 含3%、5% PEG模型體系徑向分布函數(shù)示意圖Fig.7 Schematic diagram of the radial distribution function of the 3% and 5% PEG model system

2.4 體系能量隨時(shí)間步長(zhǎng)的變化

根據(jù)分子力學(xué)原理，在COMPASSII力場(chǎng)中的總能量項(xiàng)可以表示為動(dòng)能、勢(shì)能、非鍵能之和：Etotal=Epotential energy+Enon-bond energy+Ekinetic energy，由范德華力和庫(kù)侖力組成非成鍵相互作用能。為進(jìn)一步分析在相變過(guò)程中體系能量變化對(duì)相變溫度的影響, 將探討力場(chǎng)能量項(xiàng)對(duì)聚合物類復(fù)合蓄冷劑的影響機(jī)理。

圖8 模擬體系中非鍵能隨模擬時(shí)長(zhǎng)的變化Fig.8 Variation of non-bonding energies with simulation time in the simulated system

圖8（a）所示為不同質(zhì)量分?jǐn)?shù)的PEG體系的非鍵能隨時(shí)間的變化，體系非鍵能受不同質(zhì)量分?jǐn)?shù)的PEG影響較大，體系能量大小關(guān)系為：0%PEG>5%PEG>3%PEG。圖8（b）所示為3%PEG模型體系在不同溫度下的非鍵能示意圖，可知，非鍵能隨溫度的降低而降低，即體系的氫鍵等作用力隨溫度的降低而減弱，這是因?yàn)殡S著溫度的降低，分子所獲得的動(dòng)能減少，流動(dòng)性減弱，產(chǎn)生新氫鍵的幾率幅減小，當(dāng)蓄冷劑發(fā)生相變時(shí)，由液態(tài)變?yōu)楣虘B(tài)，分子流動(dòng)性幾乎為零。

2.5 均方根位移（MSD）與自擴(kuò)散系數(shù)

在分子動(dòng)力學(xué)計(jì)算中，自擴(kuò)散系數(shù)反映體系內(nèi)粒子流動(dòng)性強(qiáng)弱，可用于描述模擬體系在不同溫度下的擴(kuò)散特性。自擴(kuò)散系數(shù)可由速度自相關(guān)函數(shù)（velocity autocorrelation function,VACF）或均方根位移[26-27]（mean square displacement, MSD）來(lái)表示，本文采用均方根位移來(lái)分析自擴(kuò)散系數(shù)。

均方根位移：

（3）

擴(kuò)散系數(shù)：

（4）

氫鍵的存在會(huì)對(duì)蓄冷劑相變過(guò)程產(chǎn)生重要的影響，如圖9所示。隨著時(shí)間步長(zhǎng)的增加，PEG分子MSD曲線呈上升趨勢(shì)。隨著溫度的降低，MSD曲線逐漸降低，可知，在263.15 K時(shí)MSD曲線發(fā)生突變，隨著溫度的不斷降低，體系由液體狀態(tài)逐漸變?yōu)楣腆w狀態(tài)，當(dāng)溫度降至253.15 K時(shí)，曲線逐漸趨于平緩，當(dāng)溫度低于243.15 K時(shí)，體系MSD曲線斜率趨向于零，說(shuō)明此時(shí)溶液已完全由液態(tài)轉(zhuǎn)為固態(tài)。當(dāng)T=243.15 K時(shí)，MSD值隨時(shí)間增加，3%PEG-10%丙三醇-20%氯化銨水溶液最終穩(wěn)定在6 ?附近；在5%PEG-10%丙三醇-20%氯化銨水溶液中，最終穩(wěn)定在7.5 ?附近。在非相變溫度變化時(shí)，MSD隨溫度的降低有減小的趨勢(shì)。隨著溫度的進(jìn)一步降低，分子的能量太小，在較低的溫度下無(wú)法打破氫鍵。當(dāng)溶液達(dá)到結(jié)晶狀態(tài)時(shí)，MSD值達(dá)到上限。整個(gè)冷卻過(guò)程中體系的相變溫度在253.15～263.15 K之間。

圖9 3%、5% PEG體系中PEG分子MSD隨時(shí)間步長(zhǎng)的變化Fig.9 Relationship between the MSD of PEG molecules in 3% and 5% PEG systems with time step changes

圖10 0%、3%、5% PEG-10%丙三醇-20%氯化銨水溶液體系擴(kuò)散系數(shù)隨溫度的變化Fig.10 Diffusion coefficient of 0%, 3%, 5% PEG-10% glycerol-20% ammonium chloride aqueous solution system changes with temperature

在此基礎(chǔ)上加入3%PEG、5%PEG，由圖10可知，此時(shí)模型體系的相變溫度分別為254.3、253.8 K，加入聚合物后，溶液的相變溫度分別降低1.7、2.2 K，這是因?yàn)镻EG聚合物分子一方面由于具有很長(zhǎng)的碳鏈，在溶液中可以限制水分子以及其它分子流動(dòng)，增加體系的黏度；另一方面PEG聚合物分子中的羥基氧以及醚基氧能夠與水分子形成氫鍵，越來(lái)越多的水分子聚集在聚合物分子附近形成小團(tuán)簇，自由水轉(zhuǎn)變?yōu)榻Y(jié)合水，使溶液的相變溫度降低。

圖11 0%、3%、5% PEG-10%丙三醇-20%氯化銨-水溶液體系密度隨溫度的變化Fig.11 0%、3%、5% PEG-10% glycerol-20% ammonium chloride water system density changes with temperature

2.6 體系密度隨溫度的變化

圖11所示為0%、3%、5%PEG-10%丙三醇-20%氯化銨-水溶液體系密度隨溫度的變化?？芍?，在NPT系綜下，隨著體系溫度的升高，溶液密度變小，原因是溶液系統(tǒng)被加熱后會(huì)經(jīng)歷一個(gè)宏觀的膨脹過(guò)程，體系分子的運(yùn)動(dòng)強(qiáng)度越來(lái)越大，會(huì)使它們彼此分離。在液態(tài)體系中，分子的運(yùn)動(dòng)較快，遷移率也較大，當(dāng)溫度改變時(shí)，液態(tài)體系的密度會(huì)隨時(shí)改變，即密度對(duì)溫度的響應(yīng)較為敏感，曲線的斜率也較大。當(dāng)體系轉(zhuǎn)變?yōu)楣腆w時(shí)，體系中各個(gè)分子的相對(duì)位置基本固定。隨著溫度的進(jìn)一步降低，分子只能從平衡位置處稍稍偏移，遷移率很小，此時(shí)密度對(duì)溫度的響應(yīng)較為遲鈍，曲線斜率也較小。相變前后密度隨溫度變化曲線的拐點(diǎn)，即為相變點(diǎn)。由圖11可知，在丙三醇-氯化銨-水溶液體系中，相變溫度的拐點(diǎn)為255.2 K，這與文獻(xiàn)[28]中的實(shí)驗(yàn)結(jié)果255.5 K吻合良好。而在添加3%、5%PEG聚合物后，預(yù)測(cè)的溶液體系相變溫度分別為254.1、253.5 K，分別下降1.1、1.7 K。

3 結(jié)論

本文在NPT系綜下采用分子動(dòng)力學(xué)方法對(duì)10%丙三醇-20%氯化銨-水模型體系在223.15～293.15 K進(jìn)行了模擬和動(dòng)力學(xué)分析，添加0%、3%、5%PEG聚合物后并預(yù)測(cè)了該體系的相變溫度。得到結(jié)論如下：

1）通過(guò)對(duì)溶液體系的成鍵與構(gòu)象分析，溫度越低，分子動(dòng)能越小，分子運(yùn)動(dòng)越平緩，體系中的O—H…H型的強(qiáng)氫鍵作用是降低蓄冷劑相變溫度的主要驅(qū)動(dòng)力。溫度對(duì)于鍵角的影響較小，而對(duì)于鍵長(zhǎng)的影響較大，O—H型羥基鍵長(zhǎng)變短的概率增加，且鍵長(zhǎng)朝著變小的方向變化。

2）計(jì)算體系徑向分布函數(shù)，隨著溫度的升高體系的有序性逐漸降低，O—H…H型氫鍵鍵長(zhǎng)位于2.3～2.5 ?之間，OW—HW、OW—HA、OW—HPEG、OW—HG、NA—HW氫鍵在溶液中占據(jù)主要位置，而OG—HW、OPEG—HW兩種類型的成鍵能力較小。

3）通過(guò)對(duì)溶液體系的MSD、自擴(kuò)散系數(shù)及密度等分析，在10%丙三醇-20%氯化銨水溶液體系中，相變溫度的拐點(diǎn)為255.2～256.0 K,與實(shí)驗(yàn)結(jié)果255.5 K吻合良好。而在添加3%、5%PEG聚合物后，預(yù)測(cè)的溶液體系相變溫度分別為254.1～254.3 K、253.5～253.8 K，分別下降1.1～1.7 K、1.7～2.2 K。

符號(hào)說(shuō)明

b——鍵長(zhǎng)，?

b0——鍵長(zhǎng)平均值，?

K——前置因子

qi、qj——原子電荷，e

φ——二面角

θ——鍵角，（°）

θ0——鍵角平均值，（°）

n——核數(shù)

D——擴(kuò)散系數(shù), ?2/ps

gαβ（r） ——表示以α原子為中心，在r和r+Δr范圍內(nèi)單位體積出現(xiàn)的β粒子的平均概率

mi、vi——離子i的質(zhì)量和瞬時(shí)速度

N——溶液體系中的離子總數(shù)

Nα、Nβ——α、β兩個(gè)原子的粒子總數(shù)

Nm——原子的總數(shù)

niβ（r） ——代表以第i個(gè)α原子為中心，在r和r+Δr范圍內(nèi)所有可能的β粒子數(shù)目

rij——原子i與j之間的距離，nm

rj（t） ——原子j在時(shí)刻t的真實(shí)的位移矢量

T——熱力學(xué)溫度，K

Vs——模擬體系的體積，nm3

εij——?jiǎng)菽軈?shù)

σij——分子間電勢(shì)變?yōu)榱愕挠邢蘧嚯x，nm

本文受河南省自然科學(xué)基金項(xiàng)目（162300410112）、河南援疆資金、哈密市科技計(jì)劃支持項(xiàng)目、新疆維吾爾自治區(qū)人才專項(xiàng)課題支持。（The project was supported by the Provincial Natural Science Foundation of Henan （No.162300410112） and Henan Aid Xinjiang Funds， Hami City Science and Technology Plan Support Project, Xinjiang Uygur Autonomous Region Talent Special Project.）