基于分子動力學模擬的聚合物類復合蓄冷劑相變機理研究

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

（1 河南理工大學機械與動力工程學院 焦作 454000；2 哈密豫新能源產業研究院 哈密 839000）

在碳達峰、碳中和的大背景下，相變材料在相變時通過潛熱的吸收、儲存和釋放，可以實現對現有能源的合理分配，并優化和提高能源利用率[1-2]。蓄冷技術[3-6]是將含有蓄冷劑的裝置在低溫環境下進行預冷，通過相變吸收并儲存大量冷量，再根據人們的需求在溫度較高的環境中能持續性的放出冷量，使周圍環境在一定時間內保持恒定的低溫狀態，用于滿足人們對食品保鮮、冷藏或其他需求的技術。近年來，蓄冷技術在冷凍加工、低溫貯藏及長途運輸等方面應用廣泛[7-8]，為人們的生活提供了極大便利。

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

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

1 模型建立與參數設定

1.1 力場選擇與模擬參數的設定

在所有的力場中，COMPASSII力場（凝聚相，能準確模擬凝聚態體系的結構和性質，用于粒子模擬研究的最優化分子勢）可將無機分子體系和有機分子體系結合在一起，能夠模擬出體系中各個獨立分子的結構、熱力學性質和振動頻率等。COMPASSII力場用于液相分子相互作用的函數形式可定義為[21]：

（1）

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

1.2 平衡狀態的判定

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

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

圖2 強氫鍵作用O—H…O和C—O…H（紅色）與弱氫鍵作用C—H…O和N—H…O（藍色）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 結果與討論

2.1 成鍵與構象分析

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

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

圖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%氯化銨-水溶液構象圖Fig.4 3% PEG-10% glycerol-20% ammonium chloride-water conformation diagram at different temperatures

圖4所示為在不同溫度下3%PEG-10%丙三醇-20%氯化銨-水溶液體系模型。從整個動力學過程來看，溫度越低，分子動能越小，分子運動越平緩。在243.15～273.15 K溫度區間內，溫度對于丙三醇和水分子的影響較顯著。丙三醇分子與水分子之間通過氫鍵結合傾向于形成小團簇，分子整體運動較為平緩。溫度對于水分子之間的結合力影響較大。當T=273.15 K時，分子運動劇烈，分子的擴散速度較大。隨著溫度的降低，水分子的運動速度減小。當T=243.15 K時，水分子的運動速度相較于273.15 K時大幅降低。

2.2 鍵長與鍵角隨溫度的變化

圖5所示為丙三醇分子與水分子構象示意圖。隨著溫度的降低，物質分子的微觀形態也會發生變化，如鍵長與鍵角也會隨著相變過程溫度的變化而發生變化，通過分析相變過程中分子的鍵長與鍵角的變化來進一步探究相變過程機理。

圖5 體系中丙三醇分子和水分子鍵長鍵角構象圖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鍵長隨溫度的變化Fig.6 Variation of H—O—H bond angle and O—H bond length with temperature in the system

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

2.3 徑向分布函數（RDF）

徑向分布函數（radial distribution function, RDF）是整個模型體系的區域密度和平均密度的比值,是描述以某一中心原子為中心周圍出現另一種原子的概率，可用來表征物質的有序性。gαβ（r）數學表達式為：

（2）

當r較小時，gαβ（r）可用來表征原子的堆積情況和各個原子對的成鍵距離，峰值高且尖，則有序性強，原子之間聯系較緊密。在溶液體系中，存在于分子間的非鍵合力通常可分為：氫鍵作用力和范德華力，前者的作用范圍為 2.6～3.1 ?，后者的作用范圍 3.1～5.0 ?。

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

2.4 體系能量隨時間步長的變化

根據分子力學原理，在COMPASSII力場中的總能量項可以表示為動能、勢能、非鍵能之和：Etotal=Epotential energy+Enon-bond energy+Ekinetic energy，由范德華力和庫侖力組成非成鍵相互作用能。為進一步分析在相變過程中體系能量變化對相變溫度的影響, 將探討力場能量項對聚合物類復合蓄冷劑的影響機理。

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

圖8（a）所示為不同質量分數的PEG體系的非鍵能隨時間的變化，體系非鍵能受不同質量分數的PEG影響較大，體系能量大小關系為：0%PEG>5%PEG>3%PEG。圖8（b）所示為3%PEG模型體系在不同溫度下的非鍵能示意圖，可知，非鍵能隨溫度的降低而降低，即體系的氫鍵等作用力隨溫度的降低而減弱，這是因為隨著溫度的降低，分子所獲得的動能減少，流動性減弱，產生新氫鍵的幾率幅減小，當蓄冷劑發生相變時，由液態變為固態，分子流動性幾乎為零。

2.5 均方根位移（MSD）與自擴散系數

在分子動力學計算中，自擴散系數反映體系內粒子流動性強弱，可用于描述模擬體系在不同溫度下的擴散特性。自擴散系數可由速度自相關函數（velocity autocorrelation function,VACF）或均方根位移[26-27]（mean square displacement, MSD）來表示，本文采用均方根位移來分析自擴散系數。

均方根位移：

（3）

擴散系數：

（4）

氫鍵的存在會對蓄冷劑相變過程產生重要的影響，如圖9所示。隨著時間步長的增加，PEG分子MSD曲線呈上升趨勢。隨著溫度的降低，MSD曲線逐漸降低，可知，在263.15 K時MSD曲線發生突變，隨著溫度的不斷降低，體系由液體狀態逐漸變為固體狀態，當溫度降至253.15 K時，曲線逐漸趨于平緩，當溫度低于243.15 K時，體系MSD曲線斜率趨向于零，說明此時溶液已完全由液態轉為固態。當T=243.15 K時，MSD值隨時間增加，3%PEG-10%丙三醇-20%氯化銨水溶液最終穩定在6 ?附近；在5%PEG-10%丙三醇-20%氯化銨水溶液中，最終穩定在7.5 ?附近。在非相變溫度變化時，MSD隨溫度的降低有減小的趨勢。隨著溫度的進一步降低，分子的能量太小，在較低的溫度下無法打破氫鍵。當溶液達到結晶狀態時，MSD值達到上限。整個冷卻過程中體系的相變溫度在253.15～263.15 K之間。

圖9 3%、5% PEG體系中PEG分子MSD隨時間步長的變化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%氯化銨水溶液體系擴散系數隨溫度的變化Fig.10 Diffusion coefficient of 0%, 3%, 5% PEG-10% glycerol-20% ammonium chloride aqueous solution system changes with temperature

在此基礎上加入3%PEG、5%PEG，由圖10可知，此時模型體系的相變溫度分別為254.3、253.8 K，加入聚合物后，溶液的相變溫度分別降低1.7、2.2 K，這是因為PEG聚合物分子一方面由于具有很長的碳鏈，在溶液中可以限制水分子以及其它分子流動，增加體系的黏度；另一方面PEG聚合物分子中的羥基氧以及醚基氧能夠與水分子形成氫鍵，越來越多的水分子聚集在聚合物分子附近形成小團簇，自由水轉變為結合水，使溶液的相變溫度降低。

圖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系綜下，隨著體系溫度的升高，溶液密度變小，原因是溶液系統被加熱后會經歷一個宏觀的膨脹過程，體系分子的運動強度越來越大，會使它們彼此分離。在液態體系中，分子的運動較快，遷移率也較大，當溫度改變時，液態體系的密度會隨時改變，即密度對溫度的響應較為敏感，曲線的斜率也較大。當體系轉變為固體時，體系中各個分子的相對位置基本固定。隨著溫度的進一步降低，分子只能從平衡位置處稍稍偏移，遷移率很小，此時密度對溫度的響應較為遲鈍，曲線斜率也較小。相變前后密度隨溫度變化曲線的拐點，即為相變點。由圖11可知，在丙三醇-氯化銨-水溶液體系中，相變溫度的拐點為255.2 K，這與文獻[28]中的實驗結果255.5 K吻合良好。而在添加3%、5%PEG聚合物后，預測的溶液體系相變溫度分別為254.1、253.5 K，分別下降1.1、1.7 K。

3 結論

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

1）通過對溶液體系的成鍵與構象分析，溫度越低，分子動能越小，分子運動越平緩，體系中的O—H…H型的強氫鍵作用是降低蓄冷劑相變溫度的主要驅動力。溫度對于鍵角的影響較小，而對于鍵長的影響較大，O—H型羥基鍵長變短的概率增加，且鍵長朝著變小的方向變化。

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

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

符號說明

b——鍵長，?

b0——鍵長平均值，?

K——前置因子

qi、qj——原子電荷，e

φ——二面角

θ——鍵角，（°）

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

n——核數

D——擴散系數, ?2/ps

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

mi、vi——離子i的質量和瞬時速度

N——溶液體系中的離子總數

Nα、Nβ——α、β兩個原子的粒子總數

Nm——原子的總數

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

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

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

T——熱力學溫度，K

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

εij——勢能參數

σij——分子間電勢變為零的有限距離，nm

本文受河南省自然科學基金項目（162300410112）、河南援疆資金、哈密市科技計劃支持項目、新疆維吾爾自治區人才專項課題支持。（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.）