張智奇 錢勝 王瑞金 朱澤飛
(杭州電子科技大學機械工程學院, 杭州 310018)
納米流體中懸浮的納米顆粒可以增強其導熱性能已經得到廣泛認可, 然而納米流體顆粒增強傳熱的機理目前尚不清楚. 研究表明, 納米顆粒的聚集是納米流體導熱系數增大的重要機制, 而且納米顆粒聚集的形態對納米流體的導熱系數有重要影響, 但是目前的導熱系數模型大多是建立在Maxwell有效介質理論的“靜態”和“均勻分散”假設基礎上. 本文用平衡分子動力學模擬Cu-Ar納米流體, 采用Green-Kubo公式計算導熱系數, 采用Schmidt-Ott關系式計算不同聚集形態下的分形維數. 對比導熱系數與分形維數可以發現: 在相同體積分數下, 較低的分形維數會有更高的導熱系數, 分析了分形維數與導熱系數的定量關系. 此外, 通過徑向分布函數可以看出納米顆粒緊密聚集與松散聚集的差異, 基液分子在納米顆粒附近的納米薄層中處于動態平衡狀態. 研究結果有助于理解納米顆粒聚集形態對導熱系數的影響機理.
近20年來, 納米流體作為傳熱介質能極大地提高傳熱性能, 引起了眾多學者的興趣. 1995年,Choi[1,2]發現在液體中加入納米顆粒大大提高了液體的導熱性能, 尤其是在低納米顆粒體積分數的情況下(一般小于0.01%)尤為明顯. 所謂納米流體就是納米顆粒穩定懸浮在液體中形成的膠體. 通常使用有效導熱系數來評估納米流體的導熱性能, 有效導熱系數模型最初是根據Maxwell有效介質理論建立的[3], 模型假設納米顆粒在液體中是靜態和均勻分散的. 此后, Hamilton 和 Crosser[4], Xuan等[5], 王補宣等[6]以及Prasher等[7]分別考慮了非球形納米顆粒、納米顆粒的布朗運動、納米顆粒外的納米薄層以及納米顆粒的聚集, 建立了新的有效導熱系數模型. 但是, 以上這些模型都是在Maxwell模型的框架基礎上得到的, 并不能擺脫納米顆粒在液體中是靜態和均勻分散的假設. 目前的實驗結果、理論結果和數值模擬結果之間尙有不少矛盾之處, 原因是迄今為止納米流體強化傳熱的微觀機理尚不清楚. Lin等[8]研究了球形顆粒和柱狀顆粒納米流體的導熱, 發現柱狀顆粒的納米流體導熱強于球形顆粒的納米流體. 李屹同等[9]針對ZnO-水納米流體, 研究了納米顆粒周圍液體分子的密度分布, 解釋了納米層提高納米流體導熱系數的機理;但是Xue等[10]的結論是納米顆粒周圍的納米層對熱傳遞沒有影響, 如果考慮界面熱阻, 納米層甚至會降低熱導率. Keblinsk等[11]對一系列實驗結果的綜合分析表明, 布朗運動、納米層和近場輻射都不是提高納米流體熱導率的關鍵原因, 最根本的是納米顆粒的聚集.
在實驗中很容易觀察到納米顆粒的聚集[12,13],顆粒聚集對于導熱系數的影響方面已經有不少成果[14-25]. Sedighi和 Mohebbi[14]通過分子動力學研究得出, 在一定的納米顆粒濃度下, 納米顆粒的聚集可以提高納米流體的導熱系數. Hong和Kim[15]的研究也表明納米顆粒的凝膠狀聚合可以大大增加導熱系數. 最令人驚訝的是含有Fe3O4納米顆粒的磁流體在外磁場的作用下形成直鏈狀聚集, 該方向上的導熱系數是基液的三倍[16], 比其他聚集狀態下的導熱系數大得多. 盡管如此, 一些研究工作表明, 由于聚集形態各異, 納米粒子的聚集并不一定能提高納米流體的導熱系數[17]. Thaseem和Christopher[18]也表示, 在不同聚集形態下, 5% 體積分數下的納米流體的導熱系數可以增加51%或減少32%. 為了探索顆粒聚集對有效導熱系數的影響規律, 研究者已經建立了一些導熱系數模型, 如Xuan等[19]較早在Maxwell有效介質理論[3]基礎上建立了考慮顆粒聚集的有效導熱系數模型;Feng等[20]建立了一個考慮納米顆粒聚集的新模型, 研究表明納米粒子越小, 聚集越多, 有效熱導系數越大; Wang等[21]在 Hamilton[4]模型的基礎上, 用聚集體等效半徑和分形維數建立了有效導熱系 數 模 型; Evans等[22]在 Bruggeman模 型 和Maxwell-Garnett理論的基礎上, 建立了一個考慮納米顆粒聚集的新模型, 并用蒙特卡羅直接模擬得到了驗證; Prasher等[23]將聚集的顆粒分成在骨架和在死端兩類, 分別考慮對導熱的影響后進行疊加, 得到了與聚集體慣性半徑、分形維數和化學維數相關的導熱系數模型, 最后與模擬結果比較時也只考慮了分形維數 df= 1.8 一種情況; Xiao 研究組[24,25]以分形理論為基礎建立了有效導熱系數的分形模型, 模型包含靜態部分和動態部分, 靜態部分根據Maxwell-Garnett模型得到, 動態部分考慮了旋轉橢球的Brownian運動引起的對流傳熱, 其中旋轉引起的對流傳熱和尺寸分布都與分形維數相關. 總之, 目前已經建立的考慮顆粒聚集的導熱系數模型, 依然以Maxwell有效介質理論為基礎,再考慮聚集對尺寸分布、納米層和布朗運動等的影響, 沒有涉及納米流體中聚集體的分形維數對導熱系數的影響的定量分析, 也沒有驗證靜態部分和動態部分簡單疊加的正確性.
聚集類型有兩種[26], 分別為克服第一勢能壁壘的緊密聚集 (reaction limited aggregation, RLA)和克服第二勢能壁壘的松散聚集(diffusion limited aggregation, DLA). Gaganpreet[27]從介觀角度分析了納米流體的傳熱, 結果表明: 松散聚集情況下納米顆粒聚集體分形維數較低、規模較大、導熱系數較高, 且松散聚集體生長迅速, 但容易被剪切稀化. Kang 等[28]和 Duan 等[29]的研究也表明, 聚集形態對導熱系數和黏度的影響很大, 更大規模的聚集會導致更大的導熱系數和更大的黏度, 產生非牛流特征, 而超聲處理后, 聚集結構被打破, 導熱系數減小, 黏度也會減小, 非牛頓流現象就會消失. Lin等[30]研究了Al2O3納米顆粒在湍流作用下聚集和破碎對導熱的影響, 得到了努賽爾數與聚集體平均尺寸的關系表達式.
分子動力學模擬是研究納米流體傳熱強化的有效方法, 但由于計算量龐大, 只適合于模擬小系統. 最近幾年有使用分子動力學模擬來探究納米粒子聚集行為對納米流體導熱系數影響的報道[14], 在分子動力學模擬中, 計算納米流體導熱系數的方法有: 平衡態分子動力學模擬和非平衡態分子動力學模擬. 平衡態分子動力學模擬基于Green-Kubo公式, 非平衡態分子動力學模擬基于導熱系數的定義. 此外, 不同分子之間的勢函數難以確定和評估,Cu-Ar納米流體的勢函數有經驗值, Kang等[31]利用Lennard-Jones (L-J)勢描述了氬原子之間的相互作用, 并采用嵌入原子勢 (embedded atom method,EAM)描述了銅原子之間的相互作用. Lee等[32]對納米流體中顆粒聚集和非聚集狀態下的熱導系數進行了分子動力學模擬, 結果表明納米顆粒聚集時的熱導系數比非聚集的情況下提高了35%.李凌等[33]研究了Cu-Ar納米流體時也得到了類似結論.
綜上可知, 直鏈狀聚集時導熱系數比其他聚集狀態大得多, 說明聚集形態會對導熱系數產生重要影響. 然而, 迄今為止還沒有納米顆粒聚集形態對導熱系數的影響的定量研究, 也未見導熱系數模型采用靜態部分和動態部分簡單疊加的正確性驗證.本文針對Cu-Ar納米流體, 用分子動力學模擬的方法, 定量分析分形維數對納米流體導熱系數的影響, 揭示添加納米顆粒增強納米流體傳熱系數的內在機制.
與文獻 [14, 34]一樣, 本文選用 Cu-Ar納米流體來研究導熱系數, 這是采用分子動力學手段研究不同聚集形態下納米流體導熱系數的理想選擇, 因為在Cu-Ar納米流體中僅需考慮兩體勢, 而不用像以水為基液時還需要考慮點電荷等因素, 大大減小了計算量. 本文模擬中采用L-J勢函數, 其與液氬的實驗數據非常吻合. 分子動力學是基于牛頓第二定律來計算模型中每個原子的運動狀態的一種數值模擬方法, 即

在模擬中選取Ar作為基液, Cu納米顆粒為面心立方結構, 其晶格常數為 3.165 ?[34]. 為了簡便起見, 假設Cu納米顆粒中的所有Cu原子具有相同的速度和加速度[32]. Ar-Ar之間采用L-J勢函數:

對于Cu原子之間的相互作用, 可以采用更加精準的 EAM勢函數[35]. 經過試算比較, 發現Cu原子之間用EAM勢和用L-J勢時結果相差無幾, 但用L-J勢時計算量明顯減小, 故采用用LJ勢. Ar原子與Cu原子之間的作用可以用Berthlot混合法則求得[36]:

此外, 三種原子之間的作用關系如圖1所示,所有顆粒的截斷半徑統一設置為10 ?, 具體結果列于表1.

圖1 不同原子間的相互作用勢Fig. 1. Potential functions between various particles interaction.

表1 勢函數參數Table 1. Parameters of potential functions.
平衡分子動力學與非平衡分子動力學都可以計算納米流體的導熱系數, 本文選擇平衡分子動力學, 原因有: 1)在平衡分子動力學中, 可以用單次模擬來計算納米流體的導熱系數, 而在非平衡分子動力學中, 需要模擬多個體系求得多個溫度梯度;2) 當系統處于非平衡狀態時, 非平衡分子動力學的計算穩定性不如平衡分子動力學, 而且當熱通量過大時可能會引起模擬體系的崩潰; 3)與非平衡分子動力學方法相比, 平衡分子動力學方法計算納米流體所需要的模擬體系的體積更小, 可大大減少計算量[37].

顆粒聚集體通常是由一定數量的基礎顆粒聚集成隨機的微結構. 由于聚集過程復雜, 納米流體中的大多數聚集體具有較廣的尺寸分布[38]. 盡管它們的大小、形狀、回轉半徑和顆粒密度各不相同,但這些具有復雜幾何形狀的聚集體可以用分形維數去衡量[39].顆直徑為的顆粒組成的聚集體,其回轉半徑為, 則Schmidt-Ott公式可表示為

在模擬中先使用NVT進行弛豫, 再在微正則系綜(NVE)下收集數據, 并在三個方向上均使用周期邊界條件[41]. 每個時間步各原子的位置與速度采用Verlet算法, 這一算法在Ar和Cu體系中具有良好的穩定性[42]. 在分子動力學模擬中, 時間步長是影響計算量、準確性和穩定性的關鍵, 太小的時間步長會導致計算量巨大, 而太大的時間步長會導致系統崩潰, 對Cu-Ar納米流體在85 K的恒溫條件下進行經過一系列試驗計算后, 確定時間步長設置為 10 fs.
在計算Cu-Ar納米流體的導熱系數之前, 先計算液氬的導熱系數, 以驗證計算方法和程序.Sarkar和 Selvam[43]研究表明, 當顆粒數目大于500時, 系統的導熱系數與系統中的顆粒數無關.本文中, 4000顆Ar原子以面心立方結構(FCC)分布在 5.72 nm × 5.72 nm × 5.72 nm 的模擬盒子中. 在弛豫階段, 采用Nose-Hoover熱浴使得系統的溫度保持在85 K. 在收集數據計算導熱系數階段, 撤去熱浴, 以免影響熱流的計算[39]. 熱流計算公式如下:

由于模擬是針對的離散時間的, 故(3)式可以離散為

如圖2所示, 歸一化的熱流自關聯函數很快收斂到0, 這樣就可以進行數據采樣了. 圖3所示為每次計算得出的導熱系數, 圖中每根細線代表一次計算, 粗線代表最終的平均結果. 由于氬的熔點和沸點非常接近, 分別為84和87.5 K, 為此本文模擬85 K 下的液氬的導熱系數, 數值為 0.131 W/(m·K),與文獻[30]中液氬在85 K下實驗得到的0.129 W/(m·K)相比, 誤差僅為 1.55%.

圖2 歸一化的熱流自關聯函數Fig. 2. Normalized heat current autocorrelation function.
由于Cu-Ar納米流體沒有實驗得到的導熱系數, 故通過與已發表的數值模擬結果對比來判斷模型的有效性. 如圖4所示, 在模擬盒子的中心位置,用一顆由327顆銅原子組成的粒徑為1.926 nm的銅納米顆粒代替原本位置上共計79顆的Ar原子,則納米顆粒的體積分數為2%. 銅原子堆積方式為面心立方堆積. 圖5為Cu-Ar納米流體的模型. 采用與計算液氬導熱系數相同的方式計算前文所描述的2%體積分數的Cu-Ar納米流體, 其結果為0.165 W/(m·K), 與文獻 [30]中得到的結果相吻合, 驗證了計算模型和程序的有效性.

圖3 液氬在 85 K 時的導熱系數Fig. 3. Thermal conductivity of Ar at 85 K.

圖4 銅納米顆粒模型Fig. 4. Model of Cu nanoparticle.

圖5 Cu-Ar納米流體模型Fig. 5. Cu-Ar nanofluid model.
納米顆粒周圍的納米薄層被認為是納米流體導熱系數的重要影響因素之一, 但仍有不少研究人員提出相反的觀點[9,10,14,34]. 為了揭示納米顆粒周圍納米薄層增強導熱的內在機理, 采用分子動力學模擬來獲取在能量最小化后運行10個時間步的Cu納米顆粒周圍Ar原子的分布, 如圖6所示. 由于Cu原子與Ar原子之間更強的相互作用, 部分Ar原子被吸引聚集在銅納米顆粒周圍. 納米顆粒周圍0.5 nm的Ar原子分布如圖7(a)所示, 在該時間點2000時間步后, 這些原子的位置如圖7(b)所示. 由圖可見部分Ar原子會離開原來的納米薄層, 但是會有大約相同數量的Ar原子被吸附, 納米薄層達到一種動態平衡. 這可能是納米薄層增強傳熱的關鍵機理, 因為這種動態平衡下納米顆粒與Ar原子間可以不斷產生能量交換, 也因此用界面熱阻的觀點解釋顯得不合理[10].

圖6 能量最小化后的納米薄層結構Fig. 6. Nanolayer at minimizational system energy.

圖7 銅納米顆粒周圍的 Ar原子初始分布 (a)和 2000 步時的分布(b)Fig. 7. Initial distribution (a) and at 2000th time-step (b) of Ar-atoms around Cu-nanoparticle.
當納米顆粒間的斥力較小而引力較大時, 顆粒會產生聚集. 如引言中所述, 聚集方式有緊密聚集RLA和松散聚集DLA兩種[26]. 當電子濃度足夠高時, 靜電引力克服顆粒之間的斥力時會產生緊密聚集, 而當電子濃度中等情況下, 不能完全克服顆粒之間的斥力則會產生松散聚集. 在本文的模擬體系中, 采用兩顆相同大小的納米銅顆粒進行模擬.圖8(a)和圖8(b)分別顯示了緊密聚集和松散聚集兩種聚集, 從圖8可以清楚地看出: 松散結構中的兩個銅納米顆粒之間存在一層Ar原子, 而在緊密聚集中的兩個納米顆粒之間有直接接觸.
徑向分布函數可以解釋為局部密度(local density)與總體密度 (bulk density)之比, 寫作:


圖8 納米顆粒聚集 (a)緊密聚集; (b)松散聚集Fig. 8. Aggregation of nanoparticles: (a) Compact aggregation; (b) loose aggregation.
當模擬體系中存在兩顆納米顆粒時, 情況發生了變化. 圖10(a)所示為兩顆納米顆粒未發生聚集時的對分布函數, 圖10(b)所示為兩顆納米顆粒發生松散聚集時的對分布函數, 觀察可見兩者幾乎相同, 而且與圖10(b)中存在單顆納米顆粒的情形十分類似. 在圖10(b)中, 兩顆納米顆粒中心的距離經計算為22.358 ?, 而圖中第一個峰出現的位置距納米顆粒的中心為11.142 ?. 前者略大于后者的2倍. 這一現象表明, 松散聚集可以理解為納米顆粒及其周圍的納米薄層作為一個整體產生聚集,其中納米薄層占6.358 ?, 納米顆粒的粒徑為16 ?.而在如圖10(c)所示的緊密聚集中, 兩顆納米顆粒質心的距離為15.237 ?, 略小于納米顆粒的粒徑,這與圖8(a)和圖8(b)所示的情形一致.
一些研究表明, 在分散狀態下納米顆粒的形狀會對于強化傳熱產生影響[7,44-46]. 然而, 在聚集狀態下納米顆粒形態對傳熱的影響尚無定量的研究.通過分子動力學模擬計算了Cu-Ar納米流體的導熱系數, 并根據模擬得到的納米顆粒的質心位置計算分形維數. 為了減少計算量, 本文選用包含3—6顆納米顆粒的小系統. 在模擬中, 納米顆粒的大小根據其體積分數計算確定, 以便在只改變體積分數的情況下, 模擬盒中存在相同數量的納米粒子. 計算結果列于表2和表3中.

圖9 對分布函數 (a)液氬體系; (b) Cu-Ar 納米流體Fig. 9. Pair correlation function of fluid Ar (a) and nanofluid Cu/Ar (b).

圖10 兩顆納米顆粒下納米流體的對分布函數 (a)分散; (b)松聚集; (c)緊聚集Fig. 10. Pair correlation function of nanofluid with two nanoparticles: (a) Dispersed nanoparticle; (b) loose aggregation; (c) compact aggregation.
從表2中可以看出, 分形維數在1—2之間,這與大多數實驗結果不符(大部分為1.4—2.7范圍), 其主要原因是本文模擬體系中的納米顆粒數少(最多6顆), 而實驗測量中是有成百上千的納米顆粒聚集在一起. 另外從表2還可以看出, 體積分數越大, 導熱系數越大, 但其增量隨體積分數的增大而減小, 這與已發表的研究成果相一致[47]. 另外,分形維數越小, 對應的導熱系數越大, 但由于納米顆粒大小的不同, 會存在一些例外, 如對比表2與表3中4%體積分數下的結果, 可以發現分形維數為1.60的6個顆粒的導熱系數為0.278 W/(m·K),而分形維數為1.59的4個顆粒的導熱系數為0.246 W/(m·K), 其原因是同體積分數下, 6 個納米顆粒的粒徑更小. 從具有相同體積分數的3個顆粒和4個顆粒的比較中也可以得出同樣的結論. 從圖11中可更清楚地看出, 在體積分數和納米顆粒尺寸一定的情況下, 隨著分形維數的增加, 熱導系數幾乎線性減小, 而且這種關系隨著體系中顆粒數的增大更加明顯.

表2 3—4 顆納米顆粒下的導熱系數與分形維數Table 2. Thermal conductivity for various aggregation morphologies with 3?4 particles.

表3 6顆粒4%體積分數時導熱系數與分形維數Table 3. Thermal conductivity for various aggregation morphologies with 6 particles.

圖11 分形維數與導熱系數對比 (3 顆為實線, 4 顆為虛線)Fig. 11. Thermal conductivity at various fractal dimensions(3-p solid line, 4-p dash line).
現有的研究成果表明, 納米顆粒的聚集對納米流體的導熱性有著很大的影響. 采用平衡態分子動力學模擬并通過Green-Kubo公式計算熱導系數;隨后, 從模擬中獲取納米顆粒的質心位置并通過Schmidt-Ott公式計算出聚集體的分形維數; 此外,還分析了單顆納米顆粒和兩顆納米顆粒情況下的對分布函數. 本文得出以下結論:
1)納米顆粒的聚集方式可分為緊湊聚集(RLA)和松散聚集(DLA)兩種, 松散聚集情況下, 兩顆納米顆粒之間存在著納米薄層, 而緊密聚集情況下兩顆納米顆粒直接接觸;
2)納米顆粒周圍納米薄層的厚度可以從徑向分布函數中讀出, 此外, 納米薄層中的基液分子達到動態平衡, 這應該是納米層增強傳熱的內在機理之一, 也說明界面熱阻來解釋邊界層傳熱是不合理的;
3)納米顆粒體積分數越大, 納米流體的導熱系數就越大, 但其增量會隨體積分數的增大而減小,如果納米顆粒的體積分數和粒徑均固定, 導熱系數隨聚集體分形維數的減小而幾乎線性增大, 且較小尺寸的納米顆粒聚集更有利于提高納米流體的導熱系數.