醇類有機物熱傳導的分子動力學模擬及微觀機理研究

劉萬強，楊帆，袁華，張遠達，易平貴，周虎

（湖南科技大學化學化工學院，理論有機化學與功能分子教育部重點實驗室，湖南省高校QSAR/QSPR重點實驗室，功能膜材料湖南省工程研究中心，湖南湘潭411201）

引 言

傳熱是化工生產領域最基本的問題之一，而熱導率(thermal conductivity, λ)是表征物質導熱性能的一個重要物性參數，λ 越大，導熱性能越好。根據熱導率的數值，可以推算出導熱介質所傳導或者交換的熱量。因此，熱導率是化工、石油和能源等工程領域需要的基礎數據。因此有機物熱導率的準確測量和計算具有非常重要的理論和實用價值[1-3]。醇類有機物是最常見的一類有機物，對醇類化合物的傳熱過程和熱導率的研究不僅具有實際意義，而且可以為其他有機物的傳熱研究提供參考。

然而目前熱導率的實驗測量耗時較長并且對設備要求高[4]。因此研究者們提出了一些估算有機物熱導率的方法。一類是從經驗發展起來的估算公式，但該方法適用的有機物結構和種類有限且計算誤差較大，并且在計算時常常需要該物質的其他性質的實驗數據[5-7]。另一類方法基于定量構效關系(quantitative structure-property relationship，QSPR)模型，這種方法計算誤差較小，但必須先建立分子參數和熱導率之間的定量模型才能使用，難以滿足實際需要[8-9]。

近年來，分子動力學(molecular dynamics，MD)模擬方法應用于液體的熱傳導以及微觀過程研究已取得一些進展[2,10-17]。表1 列出了近年來采用MD 方法對醇類有機物導熱過程的模擬和熱導率的計算結果。可以發現：(1)研究的化合物的種類較少；(2)通過傅里葉定律可以得知，溫度會對醇類有機物的熱傳導產生影響，但目前的研究沒有考察溫度對有機物導熱的影響；(3)對于液體導熱的微觀機理的研究雖取得一定進展，但仍然需要進一步深入研究[11-12]。

使用MD 方法對有機物熱傳導過程進行模擬，優點在于能夠獲得實驗中難以探測的微觀信息，例如原子水平的能量轉移和途徑等[10]，便于獲得原子間的相互作用[22-23]、化學結構[24]和振動模式[25]等對熱傳導的貢獻。因此本文采用非平衡分子動力學(non-equilibrium molecular dynamics，NEMD)方法初步探索了分子結構和溫度對醇類有機物傳熱的影響，模擬了8 種醇類有機物在4 種不同溫度下的熱傳導過程，通過對模擬過程數據的分析，初步探究不同微觀熱傳導模式之間的關系和醇類有機物熱傳導的微觀機理。

表1 部分醇類有機物的分子動力學模擬研究結果Table 1 Results of some MD simulation for calculating thermal conductivity of alcohols

1 理論基礎和模擬方法

1.1 分子模擬的構建

對一系列醇類有機物(從乙醇到正壬醇共8 種)的熱傳導過程進行了分子動力學模擬研究。模擬體系的初始構象采用Materials Studio 軟件中的Amorphous Cell 模塊構建。由于PCFF 力場(polymer consistent force field)適用于聚合物和有機物的分子動力學模擬[26-28]，因此采用PCFF 力場建立醇類有機物的分子模型。PCFF 力場相互作用的勢能主要分為鍵合能(bonded energy)和非鍵合能(non-bonded energy)兩類。其中，鍵合能包含鍵伸縮(bond stretching)、鍵 角彎 曲(angle bending)、二面 角扭 曲(torsion angle)、離面角作用力(inversion angle)和上述四項的交叉耦合項(combination cross)。非鍵合能包含分子間的庫侖(Coulomb) 作用項和9-6 Lennard-Jones勢能項(其截斷距離為12.5 ?，1 ?=0.1 nm)。這樣共有6類作用項。

圖1 非平衡分子動力學模擬模型示意圖Fig.1 Schematic diagram of NEMD simulation model

1.2 NEMD模擬模型、原理和過程

本文采用NEMD[29-30]方法對醇的熱傳導過程進行模擬，其模擬模型如圖1所示。

當系統平衡后，熱源(heat source)中粒子的熱能(thermal energy)將有規律地增加，即熱源粒子的速度將依照式(1)變化：

式中，vs是熱源粒子s 的速度；vc是熱源粒子的質心速度；Δε為常數，等于熱源中能量的增量；Erk是粒子相對于其在熱源中的質心的動能。

同樣地，冷源(heat sink)中的粒子s 的速度將按照式(2)有規律地減小：

每隔一定的時間間隔Δt，熱源和冷源的粒子將進行一次速度交換。根據冷源和熱源粒子之間速度能量轉換，可計算得到熱通量J(heat flux)，如式(3)所示：

式中，A 是模擬系統垂直于熱能傳遞方向的截斷面積。經歷一段時間的冷熱源粒子的速度交換，即可建立穩定的溫度梯度，此時可通過傅里葉定律(Fourier’s law)計算得到熱導率。如式(4)所示：

式中，dT/dz代表z方向的溫度梯度；負號表示熱能傳遞方向與溫度梯度方向相反。

基于Algaer 等[2]的研究發現，對于一個在傳熱方向上有幾納米長，并且含有上千個分子的分子液體模擬體系，模擬盒子尺寸大小對模擬結果的影響并不明顯。因此，本文不探究不同尺寸和模擬結果的相關性，并且將模擬盒子在x、y 和z 方向的初始長度分別設置為27、27 和216 ?。為了得到模擬系統的溫度梯度分布，模擬盒子沿z 方向被分為40 層。同時為了避免界面效應(surface effect)，在x、y 和z 方向都添加了周期性邊界條件(periodic boundary condition)。具體模擬過程及優化后的參數如下。

(1)使用Brown 等[31-37]針對熱通量計算修正后的LAMMPS(large-scale atomic/molecular massively parallel simulator)軟件進行分子動力學模擬。在模擬過程中系統的勢能和粒子的受力通過PCFF 力場計算得到。

(2)在弛豫前，對模擬體系進行結構預優化，避免不合理的原子間距離和受力，使體系的能量最小化。

(3)為了使體系在給定溫度下平衡穩定，首先分別在NPT、NVT 和NVE 系綜下進行了140 ps(200000步)的動力學弛豫來對體系進行預平衡。然后再分別在NVT 和NVE 系綜下分別進行700 ps(1000000步)的動力學弛豫。

(4)當體系弛豫平衡后，繼續在NVE 系綜下，根據非對稱eHEX (enhanced heat exchange)算法[38]，進行總共2 ns 的熱傳導模擬過程，并在模擬體系的z方向添加恒定熱流，其大小為Δε=0.003 kcal?mol-1?fs-1（1 cal=4.18 J）。

(5) 在弛豫過程中，根據Velocity-Verlet 方法進行數值迭代。

(6)由于模擬體系中H 原子的最小振動周期為10 fs，因此積分步長(time step)設置為0.7 fs。

(7)模擬過程在一個大氣壓下進行，并且溫度分別設置為273、288、298 和323 K。模擬體系的溫度和壓力分別由Nose-Hoover[39]熱浴法(thermostat)和壓浴法(barostat)進行調控。

(8) 體系的長程靜電相互作用(Long-ranged electrostatic interaction)根據Ewald法進行計算得到。

(9)模擬過程中當熱導率和溫度梯度收斂或穩定后，輸出模擬結果并分析數據。

1.3 熱流分解

根據不同的相互作用類型，由冷源和熱源之間動能交換得到的總熱通量Jtot，可分解為相應部分熱通量，并主要分為傳輸項(transport term，Jtrans)和相互作用項(interaction term，Jinteract)兩大類[18,23-24]。

傳輸項Jtrans主要包含由模擬體系中粒子遷移傳遞的能量，可分為動能項(kinetic energy term，Jkin)和勢能項(potential energy term，Jpot)：

式中，Us、ms和vs分別代表粒子s的勢能、質量和速度；VCV代表冷源和熱源之間的過渡區域(control volume)，求和表示過渡區域中粒子動能或勢能的總和。

相互作用項Jinteract主要由粒子間通過相互作用項進行能量傳輸。相互作用項包含非鍵合作用部分的分子間相互作用項(intermolecular interaction，Jinter)和鍵合部分的分子內相互作用項(intramolecular interaction，Jintra)。分子間相互作用項Jinter可分為范德華相互作用項JvdW和庫侖相互作用(Coulomb interaction)項JCl。分子內相互作用項Jintra包含了四個相互作用項：Jstr、Jang、Jtor和Jinv，分別對應于鍵伸縮(bond stretching)、角 彎 曲(angle bending)、二 面 角(torsion angle)和離面角(inversion angle)。在液體中，相互作用項一般占主要部分[18,24]。

式中，fX代表作用于粒子s 的某一作用類型X 的相互作用。

根據式(5)～式(7)的分析，總的熱通量Jtot可以分解為：動能項Jkin、勢能項Jpot、范德華力項JvdW、庫侖力JCl、鍵伸縮Jstr、角彎曲Jang、二面角Jtor和離面角Jinv共8個部分。根據不同類型的熱通量JX，熱導率也可被分 為 8 類 部 分 熱 導 率 (partial thermal conductivity，λX)。

除超過一定距離的粒子間相互作用對外，每一種相互作用力都可視為一條粒子間傳熱路徑(atomistic heat path)[24]，而熱能沿這個路徑進行傳輸。每一種相互作用力X 的對應的部分熱導率λX等于傳熱路徑密度(heat path density，ρX)與傳熱效率(efficiency，ΛX)的乘積。

其中，傳熱效率指單位時間內通過單位面積熱能傳遞速率。對于分子內相互作用項(X=鍵伸縮、角彎曲、二面角和離面角)，每一個多體相互作用可視為一條傳熱路徑，其密度為：

2 結果與討論

2.1 熱導率的計算

本文采用NEMD 模擬了8 種液態醇類有機物在273、288、298、323 K 下的熱傳導。表2 列出了根據式(3)、式(4)計算得到的32個熱導率計算值與實驗值的比較。從表2 可以看出，熱導率的計算值與實驗值的平均相對偏差為3.77%，最大相對偏差為4.94%，計算值與實驗值基本吻合，表明采用的NEMD 模擬方法及優化的參數用于模擬計算醇類有機物的熱導率較為理想。

從表2 可以看出，雖然熱導率的實驗值與計算值之間存在偏差，但本文的熱導率計算值與實驗值的變化趨勢一致。另一方面，熱導率的計算值和實驗值的偏差，說明模擬時存在部分偏差如質量密度和力場等仍需改進。根據質量密度的模擬值與實驗值(附錄表A1)的比較可知，最終得到的模擬值低于實驗值，但是兩者隨鏈長和溫度的變化趨勢相同。本文使用的PCFF 力場屬于全原子(all-atom)力場，分子中某些自由度(包含氫原子的鍵伸縮項、角彎曲項等)在實際中屬于基態(ground states)振動模型，而在模擬中做柔性鍵伸縮項或角彎曲項處理。因此在模擬中這類鍵伸縮項或角彎曲項可用來傳遞能量，并使熱導率的計算值與實驗值產生偏差。雖然上述誤差使熱導率的計算存在不確定性，但本文對分子傳熱機理與分子鏈長和溫度的關系的研究仍然具有意義。

表2 NEMD模擬計算得到的熱導率與實驗值［1］比較Table 2 Comparison of thermal conductivity of experimental values and simulated values of alcohols obtained by NEMD simulation

2.2 熱流分解

根據式(5)～式(7)對熱通量進行了分解，得到了不同相互作用類型的熱通量數據(數據請見附錄表A2)。

圖2 273 K(a)、288 K(b)、298 K(c)和323 K(d)8種醇的不同作用類型相互作用熱通量相對大小Fig.2 The proportion of heat flux of alcohols at 273 K(a),288 K(b),298 K(c)and 323 K(d)respectively

圖2為273、288、298和323 K時8種醇的部分熱通量相對大小示意圖。從圖2 可以看出：①動能項Jkin、二面角項Jtor和庫侖作用項JCl在其中所占比例較大。這說明在醇類分子中主要通過分子的運動、分子內成鍵原子間振動和分子間庫侖作用傳熱。②以圖2(a)為例，273 K 下從乙醇到壬醇，動能項Jkin所占的比例從18.2%逐漸增大到24.4%，但勢能項Jpot變化并不明顯。由于分子鏈的增長，使得分子間碰撞的概率增加，因此動能項Jkin所占比例增加較顯著。③分子內相互作用鍵伸縮Jstr，鍵角彎曲Jang，二面角Jtor的熱通量所占比例隨分子鏈的增長而增大。其中二面角項Jtor所占的比例最大且增加最為顯著。273 K 下三者貢獻值的總和從26.6%升高到44.4%。④除乙醇外，隨碳原子數目的增加，分子間庫侖作用項JCl和范德華作用項JvdW所占比例逐漸減小，特別是庫侖作用項JCl變化尤為顯著。273 K 下兩者所占比例之和從37.5%降低到21.2%。⑤對于乙醇分子，傳輸項中的勢能項Jpot的熱通量比例遠大于其他醇類分子。乙醇與其他分子相比較，其分子極性更大，且易形成分子間氫鍵，分子間的作用更強，分子遷移的勢能項更大，運動更劇烈，因而勢能項Jpot的熱通量比例更高。

根據上述分析發現，在醇類有機物中，當分子鏈的碳原子數為5～6 時，熱通量的主要部分將由分子間相互作用項Jinter轉為分子內相互作用項Jintra。圖3為己醇在不同溫度時部分熱通量的變化。

圖3 不同溫度下己醇的各類相互作用熱通量比較Fig.3 Comparison of heat flux of different interactions of hexanol at different temperatures

從圖3 可以看出：①在傳輸項中，動能項Jkin熱通量最大且隨著溫度的升高顯著增大(314.87～363.52 MW·m-2)，可見由于溫度升高，分子的整體運動速度加快，因此動能項的熱通量顯著增加。但是模擬體系隨著溫度升高密度減小(具體數值見附錄表A1)，分子間距離增大，勢能項Jpot急劇減小(149.18～63.02 MW·m-2)。②隨著溫度升高，范德華作用項JvdW減小(85.98～72.88 MW·m-2)，庫侖相互作用 項JCl略 有 增 加(262.54～272.83 MW·m-2)。 在PCFF 力場中，庫侖作用勢能反比于原子之間距離r，而范德華作用為9-6 Lennard-Jones 勢，與(1/r12-1/r6)呈比例關系。因此，隨著溫度的升高，原子間距離增大，范德華作用項與庫倫相互作用均減小，但庫侖作用項的變化較小。③在分子內相互作用項中，隨著溫度的升高，原子的振動加劇，導致鍵伸縮和鍵角彎曲項作用力增大，使得鍵伸縮項Jstr(111.14～128.15 MW·m-2)和鍵角彎曲項Jang(161.58～183.90 MW·m-2)都增加，二面角項Jtor減小(272.23～258.53 MW·m-2)。

綜上所述，隨著分子鏈增長，對熱流的貢獻由分子間相互作用項Jinter轉為分子內相互作用項Jintra，這與Matsubara 等[18,24]和Ohara 等[40]的研究結果基本一致。其中，二面角項Jtor和動能項Jkin所占比例隨分子鏈增長而增加，并占主導地位。此外，隨著溫度的升高，分子的整體運動速度加快，模擬系統的質量密度降低，因此，在己醇中，傳輸項的熱通量Jtrans(包含Jpot、Jkin)隨溫度升高而升高，而相互作用項Jinteract(包含Jstr、Jang、Jtor、Jinv、JvdW、JCl)則降低(具體數值見附錄表A2)。因此，最終的凈效應是己醇的熱導率隨著溫度的升高而降低。

2.3 原子傳熱路徑分析

從2.2 節分析可以看出，隨著分子鏈長的增加，分子內作用項Jintra顯著增加并占主導地位。因而進一步考察了分子內各作用項的熱導率(λX)，并統計了模擬體系各項的密度(density, ρX)，計算了各項的傳熱效率(efficiency,ΛX)。

圖4(a)為273 K 下不同醇分子內鍵伸縮項、角彎曲項、二面角項和離面角項對應的部分熱導率(partial thermal conductivity, λintra)。可以看出，鍵伸縮項部分熱導率λstr、角彎曲項λang和二面角項λtor都隨分子鏈的增長而增加。三者對整體熱導率的貢獻依次為λtor＞λang＞λstr。由于在鏈狀分子中，幾乎不存在離面角，因而λinv貢獻幾乎為零。此結論與烷烴的研究結果相似[24]。

圖4(b)為273 K 下，模擬體系中鍵伸縮項、鍵角彎曲項、二面角和離面角的密度(density,ρintra)隨分子鏈長的變化。由于分子鏈增長，分子內共價鍵數目增加；此外，通過模擬體系質量密度的分析發現，分子鏈增長模擬體系質量密度也增加了(詳見附錄表A1)。因此如圖4(b)所示，鍵伸縮項ρstr、角彎曲項ρang、二面角項ρtor和離面角項密度ρinv均隨著分子鏈的增長而增大，其大小依次為ρtor＞ρang＞ρinv＞ρstr。

圖4 醇類有機物在273 K下各分子內相互作用項的部分熱導率λintra(a)，密度ρintra(b)和傳熱效率Λintra(c)Fig.4 The partial thermal conductivity λintra(a),the density ρintra(b),and the efficiency Λintra(c)of intramolecular interaction paths for alcohols at 273 K,respectively

圖4(c)為273 K 下，不同醇分子內鍵伸縮項、鍵角彎曲項、二面角和離面角作用項的傳熱效率(efficiency,Λintra)。由于模擬體系質量密度隨分子鏈增長而增大，原子間碰撞概率增加，使得單位時間內通過單位面積的熱能傳遞速率增大。因此，從圖中可以看出，除了乙醇分子外，鍵伸縮項Λstr、鍵角彎曲項Λang和二面角項Λtor的傳熱效率都隨分子鏈的增長而增大。此結果與Matsubara 等[24]對液態烷烴的研究結果相似。此外，當分子鏈較短時，如乙醇，分子的柔性較小，二面角的作用更強，因此二面角項Λtor較大。但當分子鏈增長時，二面角項Λtor逐漸減小，最后與鍵伸縮項的Λstr幾乎相等，鍵角彎曲項Λang次之。

3 結 論

本文采用NEMD 模擬方法，模擬了8 種常見醇類有機物在不同溫度下的熱傳導過程。得到的熱導率的計算值與實驗值的相對平均偏差為3.77%，最大相對偏差小于5.00%。

通過熱流分解、部分熱導率和傳熱效率的分析，表明醇類有機物的分子結構對其熱能傳輸有較大影響。當分子鏈較短時，對醇類熱傳導起主要作用的是分子的傳輸項(動能和勢能項)及分子內的二面角項。隨著分子鏈增長，分子內的共價鍵數目增多，熱流的主要傳導方式逐漸由分子間相互作用過渡到分子內相互作用項(鍵伸縮項、角彎曲項、二面角項)傳導。

此外，溫度對醇類有機物的熱傳導具有一定影響，隨著溫度的升高，與分子和原子運動和振動相關的作用項，如動能項、二面角項、鍵伸縮項等傳導的熱流增加，而與勢能有關的作用項，如勢能項、范德華力項傳導的熱流減小，最終使得醇類有機物的熱導率隨溫度升高而減小。

本文對醇類有機物的分子動力學模擬研究和熱導率的計算，為研究液體的熱傳導過程提供了微觀依據，并為熱導率的計算提供了新的途徑，同時也有助于理解液體熱傳導的途徑和方式。