劉萬強，楊帆，袁華，張遠(yuǎn)達(dá)，易平貴，周虎

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

引 言

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

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

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

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

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

1 理論基礎(chǔ)和模擬方法

1.1 分子模擬的構(gòu)建

對一系列醇類有機(jī)物(從乙醇到正壬醇共8 種)的熱傳導(dǎo)過程進(jìn)行了分子動力學(xué)模擬研究。模擬體系的初始構(gòu)象采用Materials Studio 軟件中的Amorphous Cell 模塊構(gòu)建。由于PCFF 力場(polymer consistent force field)適用于聚合物和有機(jī)物的分子動力學(xué)模擬[26-28]，因此采用PCFF 力場建立醇類有機(jī)物的分子模型。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 非平衡分子動力學(xué)模擬模型示意圖Fig.1 Schematic diagram of NEMD simulation model

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

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

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

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

同樣地，冷源(heat sink)中的粒子s 的速度將按照式(2)有規(guī)律地減?。?/p>

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

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

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

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

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

(2)在弛豫前，對模擬體系進(jìn)行結(jié)構(gòu)預(yù)優(yōu)化，避免不合理的原子間距離和受力，使體系的能量最小化。

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

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

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

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

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

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

(9)模擬過程中當(dāng)熱導(dǎo)率和溫度梯度收斂或穩(wěn)定后，輸出模擬結(jié)果并分析數(shù)據(jù)。

1.3 熱流分解

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

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

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

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

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

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

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

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

2 結(jié)果與討論

2.1 熱導(dǎo)率的計算

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

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

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

2.2 熱流分解

根據(jù)式(5)～式(7)對熱通量進(jìn)行了分解，得到了不同相互作用類型的熱通量數(shù)據(jù)(數(shù)據(jù)請見附錄表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在其中所占比例較大。這說明在醇類分子中主要通過分子的運動、分子內(nèi)成鍵原子間振動和分子間庫侖作用傳熱。②以圖2(a)為例，273 K 下從乙醇到壬醇，動能項Jkin所占的比例從18.2%逐漸增大到24.4%，但勢能項Jpot變化并不明顯。由于分子鏈的增長，使得分子間碰撞的概率增加，因此動能項Jkin所占比例增加較顯著。③分子內(nèi)相互作用鍵伸縮Jstr，鍵角彎曲Jang，二面角Jtor的熱通量所占比例隨分子鏈的增長而增大。其中二面角項Jtor所占的比例最大且增加最為顯著。273 K 下三者貢獻(xiàn)值的總和從26.6%升高到44.4%。④除乙醇外，隨碳原子數(shù)目的增加，分子間庫侖作用項JCl和范德華作用項JvdW所占比例逐漸減小，特別是庫侖作用項JCl變化尤為顯著。273 K 下兩者所占比例之和從37.5%降低到21.2%。⑤對于乙醇分子，傳輸項中的勢能項Jpot的熱通量比例遠(yuǎn)大于其他醇類分子。乙醇與其他分子相比較，其分子極性更大，且易形成分子間氫鍵，分子間的作用更強，分子遷移的勢能項更大，運動更劇烈，因而勢能項Jpot的熱通量比例更高。

根據(jù)上述分析發(fā)現(xiàn)，在醇類有機(jī)物中，當(dāng)分子鏈的碳原子數(shù)為5～6 時，熱通量的主要部分將由分子間相互作用項Jinter轉(zhuǎn)為分子內(nèi)相互作用項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)，可見由于溫度升高，分子的整體運動速度加快，因此動能項的熱通量顯著增加。但是模擬體系隨著溫度升高密度減小(具體數(shù)值見附錄表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)呈比例關(guān)系。因此，隨著溫度的升高，原子間距離增大，范德華作用項與庫倫相互作用均減小，但庫侖作用項的變化較小。③在分子內(nèi)相互作用項中，隨著溫度的升高，原子的振動加劇，導(dǎo)致鍵伸縮和鍵角彎曲項作用力增大，使得鍵伸縮項Jstr(111.14～128.15 MW·m-2)和鍵角彎曲項Jang(161.58～183.90 MW·m-2)都增加，二面角項Jtor減小(272.23～258.53 MW·m-2)。

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

2.3 原子傳熱路徑分析

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

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

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

圖4 醇類有機(jī)物在273 K下各分子內(nèi)相互作用項的部分熱導(dǎo)率λ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 下，不同醇分子內(nèi)鍵伸縮項、鍵角彎曲項、二面角和離面角作用項的傳熱效率(efficiency,Λintra)。由于模擬體系質(zhì)量密度隨分子鏈增長而增大，原子間碰撞概率增加，使得單位時間內(nèi)通過單位面積的熱能傳遞速率增大。因此，從圖中可以看出，除了乙醇分子外，鍵伸縮項Λstr、鍵角彎曲項Λang和二面角項Λtor的傳熱效率都隨分子鏈的增長而增大。此結(jié)果與Matsubara 等[24]對液態(tài)烷烴的研究結(jié)果相似。此外，當(dāng)分子鏈較短時，如乙醇，分子的柔性較小，二面角的作用更強，因此二面角項Λtor較大。但當(dāng)分子鏈增長時，二面角項Λtor逐漸減小，最后與鍵伸縮項的Λstr幾乎相等，鍵角彎曲項Λang次之。

3 結(jié) 論

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

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

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

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