自由分子區內納米顆粒的熱泳力計算*

崔杰 蘇俊杰 王軍? 夏國棟 李志剛

1) (北京工業大學能源與動力工程學院, 傳熱強化與過程節能教育部重點實驗室, 傳熱與能源利用北京市重點實驗室, 北京 100124)

2) (香港科技大學機械及航空航天工程學系, 香港)

基于非平衡態分子動力學模擬方法, 研究了自由分子區內納米顆粒的熱泳特性.理論研究表明, 納米顆粒與周圍氣體分子之間的非剛體碰撞效應會明顯地改變其熱泳特性, 經典的Waldmann 熱泳理論并不適用,但尚未有定量的直接驗證.模擬計算結果表明: 對于納米顆粒而言, 當氣?固相互作用勢能較弱或氣體溫度較高時, 氣體分子與納米顆粒之間的非剛體碰撞效應可以忽略, Waldmann 熱泳理論與分子動力學模擬結果吻合較好; 當氣?固相互作用勢能較強或氣體溫度較低時, 非剛體碰撞效應較為明顯, Waldmann 熱泳理論與模擬結果存在較大誤差.基于分子動力學模擬結果, 對納米顆粒的等效粒徑進行了修正, 并考慮了氣體分子與納米顆粒之間的非剛體碰撞效應, 理論計算結果與分子動力學模擬結果吻合較好.

1 引 言

微細顆粒在溫度不均勻的氣體中運動時, 會受到周圍氣體分子的差異性碰撞, 從而受到一個沿溫度梯度相反方向上力的作用, 即為熱泳力.顆粒在熱泳力作用下的運動, 稱為熱泳[1?4].通常情況下,來自高溫區域的氣體分子與顆粒之間的碰撞更為劇烈, 這便導致顆粒受到的熱泳力方向與溫度梯度方向相反, 指向冷端.在某些情況下, 也可能出現熱泳力與溫度梯度方向相同的情況, 即反向熱泳[2].對熱泳現象的研究在薄膜制備、微納米制造、環境科學、氣溶膠等領域中都有重要的意義[5?10].近年來, 針對顆粒輸運現象的理論與實驗研究越來越多, 但是, 相關現象的內在機理還有待深入研究和分析[11?18].

顆粒的輸運特性與其Knudsen 數Kn (Kn =λ/ L )密切相關.其中, λ 為氣體分子的平均自由程, L 表示顆粒的特征長度.對于球形顆粒來說,特征長度為顆粒半徑R, 并使用符號KnR表示.根據Knudsen 數的大小, 可以將顆粒的動力學行為粗略地劃分為連續介質區(continuum regime)、過渡區(transition regime)和自由分子區(free mole?cule regime).在連續介質區(KnR? 1), 流體分子的自由程很小, 分子間相互作用的頻率很高, 此時, 可以把流體看作連續介質處理, 使用經典的Navier?Stokes 方程來對顆粒的輸運特性進行分析計算.Epstein[19]基于此方法研究了連續介質區的熱泳現象, 但求解過程中采用的邊界條件并不合適.Brock[20]采用完備的滑移邊界條件得到了連續介質區熱泳力的一階近似解.其他學者們也嘗試開展了高階近似解的計算[21].然而, 此方法也受到了一定的質疑[22?24].在自由分子區(KnR? 1), 由于氣體分子的平均自由程遠大于顆粒的特征長度, 被顆粒反射后的分子需飛行很長的距離才會與其他分子相碰, 因此可以忽略入射氣體分子與被顆粒反射氣體分子間的相互作用.此時, 可假設氣體分子的速度分布函數不會受到懸浮顆粒的影響.1956 年,Waldmann[25]基于氣體分子與顆粒之間的剛體碰撞模型, 得到了自由分子區中顆粒在單原子氣體中所受熱泳力的計算公式:

式中, mg為氣體分子的質量, kB為Boltzmann 常數, T 為溫度, κ 為氣體熱導率, R 為納米顆粒半徑.? T 為氣體環境的溫度梯度.越來越多的理論與實驗研究表明, 隨著KnR增大, 顆粒所受的熱泳力值逐漸接近Waldmann 公式的計算結果[2], 因此, (1)式也被稱為熱泳力的自由分子極限.在連續介質區(KnR～ 1), 氣體分子的離散特性剛開始體現, 氣體分子的速度分布受顆粒運動的影響較大, 但又不夠達到可以看作連續介質的程度, 此時,求解Boltzmann 輸運方程是十分困難的.有學者曾嘗試將自由分子區的理論擴展到過渡區, 但結果并不理想[26].Talbot 等[27]通過改變幾個Brock 公式中的參數, 使Brock 的連續介質區理論公式在KnR? 1 時收斂于(1)式.

通常情況下, 氣體分子的平均自由程為100 nm量級.對于粒徑較小的納米顆粒而言, 其特征尺寸遠小于氣體分子的平均自由程, 因此納米顆粒在氣體中的動力學行為一般處于自由分子區.Waldmann公式是基于氣體分子與顆粒之間的剛體碰撞模型得到的, 但是, 隨著顆粒尺寸減小至納米尺度, 顆粒與氣體分子之間的非剛體碰撞效應將成為影響顆粒與氣體分子之間動量傳遞的重要影響因素之一[16].文獻[17]中, Li 和Wang 基于非剛體碰撞模型, 得到了自由分子區內納米顆粒所受熱泳力的理論計算公式:

式中, mp為顆粒的質量; mr= mgmp/(mg+ mp),為氣體分子與顆粒的約化質量; Ω(1,1)*與Ω(1,2)*為無量綱碰撞積分[11].對于剛體碰撞來說, Ω(1,1)*=Ω(1,2)*= 1, (2)式退化為(1)式.換言之, (1)式為(2)式在剛體碰撞假設下的特殊情形.但是, (2)式的準確性尚未在實驗或者模擬中得到定量的驗證.一方面受限于納米尺度下測量技術的不足, 另一方面顆粒的熱泳運動一般較弱, 很容易被其布朗運動所掩蓋, 因此關于納米顆粒熱泳現象的實驗研究往往比較困難[28].分子動力學模擬方法可以直觀地研究并分析顆粒的輸運及受力特性[28?33], 在模擬中, 可以在納米顆粒上人為地施加簡諧勢來對納米顆粒進行約束, 從而明顯降低了顆粒布朗運動對其熱泳力計算的影響.

對于實際納米顆粒而言, 其形狀多為非球形(例如碳納米管等).非球形顆粒在流場中的取向會對其受力產生重要的影響, 顆粒的瞬時熱泳力的大小和方向都與其形狀和取向有關, 較為復雜[34].在自由分子區, 氣體分子的平均自由程遠大于顆粒粒徑, 因此, 氣體分子與顆粒的隨機碰撞會使得顆粒高速旋轉.顆粒的宏觀受力體現為氣體分子與顆粒之間微觀動量傳遞量在長時間內的統計平均結果.如果顆粒沒有受到較強外勢場的作用, 顆粒取向分布較為均勻, 可以采用球形顆粒模型近似研究非球形顆粒的熱泳特性.

本文基于非平衡態分子動力學模擬方法, 研究了自由分子區內球形納米顆粒的熱泳現象, 計算得到了納米顆粒所受熱泳力的大小和方向.計算結果表明: 當氣?固相互作用勢能較強或氣體溫度較低時, 納米顆粒熱泳力的計算需要考慮氣體分子與納米顆粒之間的非剛體碰撞效應; 進行顆粒粒徑修正后, 基于Li?Wang 公式(2)的計算結果與分子動力模擬結果吻合較好.

2 分子動力學建模

2.1 分子動力學模擬系統

本文采用分子動力學模擬方法計算納米顆粒在自由分子區內所受熱泳力.分子動力學模擬系統如圖1 所示, 直徑和質量分別為DP和mP的納米顆粒被浸入到氣體分子數為N 的模擬區域中, 選擇固體分子與氣體分子質量相同, 均為m.在系統y 和z 方向采用周期性邊界條件.在x 方向的邊界上采用鏡面反射邊界條件, 并在鄰近邊界的兩端分別建立溫度控制區域(高低溫熱浴), 其對應的溫度分別為Th和Tl.溫度控制區域在x 方向上的尺寸為lT= 0.1lx, 此處, lx為x 方向的系統尺寸.

圖1 分子動力學模型圖Fig.1.Snapshot for the MD simulation system.

采用經典的Lennard?Jones(L?J)勢函數來模擬氣?氣、氣?固、固?固分子間的相互作用:

式中, 參數ε 和σ 分別為L?J 勢參數, rij為原子i 和原子j 之間的距離.系統中的原子(固體原子、氣體原子)參數均采用Ar 原子的參數[35]: σ =3.405 ?, ε/kB= 114 K.截斷半徑rc= 2.5σ.研究氣?固結合強度對顆粒所受熱泳的影響發現, 氣?固結合強度εGP可以在一定范圍內變化,

其中kij為可在一定范圍內變化的參數.此外, 在固體顆粒內部相鄰原子之間額外引入finite exten?sible nonlinear elastic(FENE)勢函數以保持納米顆粒為準球形, 其勢能表達式為[30]

式中, 勢能參數設定為R0= 1.5σ; kFENE= 30ε/σ2,為FENE 勢能的彈性系數.kFENE的大小決定了納米顆粒的導熱特性, 文獻[30]對熱泳力隨kFENE變化情況進行了研究分析, 結果表明納米顆粒的熱導率對其所受熱泳力影響不大.球形納米顆粒的直徑使用下式進行計算[30]:

式中, xc為納米顆粒的質心, NS為組成納米顆粒的固體原子數量, 研究過程中選擇納米顆粒的固體原子數量至少為40, 納米顆粒的質量mP= NSm.Galliero 和Volz[30]的研究表明, 當mP/m ＞ 10 時,顆粒的熱泳特性將不再受到顆粒質量的影響.相比于氣體分子, 納米顆粒的質量與尺寸已足夠大, 因此可忽略顆粒質量對熱泳的影響[36,37].

由于納米顆粒的布朗運動, 顆粒會在系統中隨機行走, 并受到一個瞬時的曳力(一般情況下曳力的大小比熱泳力高2—3 個數量級), 為熱泳力的計算帶來較大的困難.因此, 本文在顆粒質心與系統中心引入額外簡諧回復勢能Uhar:

式中, rc和ro分別為顆粒質心與系統中心坐標,khar為簡諧回復勢能的彈性系數.作用在納米顆粒上的簡諧回復力將均勻施加在納米顆粒的每一個原子上.簡諧回復勢在系統三個方向上給顆粒所施加的回復力分別為F = (Fx, Fy, Fz).當系統達到穩定后, 可得顆粒所受到的熱泳力為

在熱泳力的計算中引入簡諧勢, 可以一定程度限制顆粒的隨機運動, 消除顆粒運動過程中所受曳力對其熱泳力的影響.文獻[32]研究了熱泳力隨彈性系數khar的變化情況, 結果表明, 只要彈性系數khar大于0.1ε/σ2, 其大小對統計結果影響就可以忽略不計.本文取khar= 30ε/σ2.

本文中, 除特別說明, 所有參數將采用無量綱形式(用上標“*”表示), 采用下列各物理量進行無量綱化: 長度x0= σ, 溫度T0= ε/kB, 時間t0=(mσ2/ε)0.5, 力F0= ε/σ.例如, 溫度的無量綱形式為T*= T/T0.

根據氣體分子運動論, 氣體分子的平均自由程為

式中, n 為氣體的分子數密度, d 為氣體分子的有效硬球直徑.根據參考文獻[38], 溫度T*= 2.63時氬原子氣體分子的有效硬球直徑約為0.94σ.為了驗證模擬系統中氣體的狀態, 從模擬系統中選擇一組氣體的狀態參數與氣體的狀態方程進行對比,

式中, P*為NVT 系綜的壓力, Bi(T*)是維里系數.取ΔT*=0.92.P*隨n 的變化情況如圖2 所示.維里系數的取值分別為= 0.3734,= 0.0421.

圖2 壓力P*隨數密度n 變化圖Fig.2.The pressure P* versus the gas molecular number density n.

分子動力學模擬中, 使用Velocity?Verlet 算法對顆粒的運動方程進行積分, 為了保證系統能量的穩定, 取時間步長Δt*= 10—3(Δt = 2 fs).初始時刻, 氣體分子隨機分布于模擬系統中, 初始速度為溫度T*下的麥克斯韋分布.系統首先在NVT 系綜下進行弛豫, 在高低溫熱源溫度的算術平均溫度下運行直到t*== 103, 達到穩定.然后, 采用速度校正法[39]將高低溫熱浴的溫度分別控制為和, 直到t*== 1.5 × 104, 從而得到一個恒定的溫度差ΔT*, 系統的采樣時間將在t*=到時間段內進行(= 8 × 105).本文所給出的任何物理量都是在此時間范圍進行的長時間統計平均得到的.為了提高結果的準確性, 在相同的熱力學宏觀狀態下, 取不同初始速度分布, 至少8 次獨立模擬結果的平均值作為最終結果.

3 有限空間效應

Waldmann 熱泳力公式的一個前提為與顆粒相互作用的氣體分子是不受空間尺寸束縛的(空間體積無窮大).在這種情況下, 氣體熱物性不會受到空間邊界的影響.然而, 在模擬研究中, 氣體分子的活動范圍總是被具有一定空間尺寸的邊界所束縛, 在自由分子區中, 氣體分子的平均自由程可能會與空間尺寸的量級相當, 此時必須考慮有限的體積效應對氣體熱物性的影響[40].取氣體環境所處的空間尺寸L 作為特征長度, 此時努森數可用KnL表示.在圖1 所示的分子動力學模型中, L 為x 方向高低溫熱浴的間距, 即L = 0.8lx.1972 年,Phillips[41]借助Chapman?Enskog 分布函數對Bo?ltzmann 輸運方程進行了求解, 得到了有限空間內顆粒的熱泳力表達式:

式中, αmt和αmn分別為顆粒表面切向和法向的動量協調系數, α1和α2分別為高低溫壁面的溫度協調系數.在本文的分子動力學模擬系統中, 可以認為氣體分子是完全協調的, 即= α2≈ 1 以及αmt= αmn≈ 1.

分 別 采 用 Phillips (11)式 和 Waldmann(1)式進行熱泳力理論計算, 并與分子動力學模擬計算結果進行對比分析.Phillips 公式的計算結果表示為FT,Phillips, 分子動力學模擬結果為FT,MD.Waldmann 公式將分別使用氣體的宏觀熱導率κ 和有效熱導率κ'計算, 計算結果分別用FT,B和FT,E表示.圖3(a)與圖3(b)分別表示在kij= 5.26 (強氣?固相互作用下)和kij= 1.0 (弱氣?固相互作用下)情況下FT/FT,B隨l 的變化圖.由圖3 可以看出, 模擬系統中“壁面”間距是影響顆粒所受熱泳力的重要因素之一, 這說明有限空間效應確實存在.FT,E與FT,Phillips吻合較好, 表明Phillips 對Wald?mann 熱泳力計算式的修正是可以通過對氣體熱導率的修正來實現的.

表1 分子動力學模擬系統的幾何特征參數Table 1.Geometric and characteristic parameters of the simulation systems.

圖3 FT, MD/FT, B, FT, E/FT, B 和FT, Phillips/FT, B 隨模擬空間尺寸的變化 (a) kij = 5.26; (b) kij =1Fig.3.Influence of the size of simulation box on FT, MD/FT, B,FT, E/FT, B and FT, Phillips/FT, B: (a) kij =5.26, (b) kij =1.

對于強氣?固相互作用(圖3(a)), 分子動力學結果FT,MD與理論結果FT,E和FT,Phillips存在較大誤差, 這種偏差是由于氣體分子與納米顆粒之間的非剛體碰撞所引起的.對于弱氣?固相互作用(圖3(b)), 氣體分子與納米顆粒之間的非剛體碰撞效應較弱, 此時, 基于剛體碰撞模型假設的Wald?mann 理論仍然適用于納米顆粒, FT,E和FT,Phillips都與分子動力學結果FT,MD吻合較好.圖4 為熱泳力隨氣體有效熱導率κ' 的變化曲線.與前文結果類似, 對于弱氣固耦合的情況(kij= 1.0), FT,E和與分子動力學結果FT,MD吻合較好; 對于強氣固耦合(kij= 5.26), 分子動力學模擬得到的熱泳力則明顯高于基于有效熱導率的理論計算結果.

圖4 熱泳力FT,MD 和FT,E 隨氣體有效熱導率κ' 的變化Fig.4.Influence of the effective thermal conductivity of the media gas on the thermophoretic forces.

4 結果與討論

4.1 系統溫度(氣-固結合強度)影響

由前文所述可知, 將氣體的有效熱導率引入到Waldmann 公式中, 可得到有限空間中氣體的Waldmann 自由分子極限.但是, 對于納米顆粒來說, 特別是當氣?固相互作用較強時, Waldmann 熱泳力計算式的誤差較大.基于分子動力學模擬結果,熱泳力隨系統平均溫度氣?固結合強度kij和溫度梯度 ? T 的變化情況如圖5 所示(圖5(a)中的插圖為在相同系統參數下本文分子動力學結果與文獻[28]中結果的對比).由圖5 可以看出, 對于納米顆粒來說, Waldmann 熱泳力計算式的適用性受系統平均溫度與氣?固結合強度kij的影響, 當系統溫度較低或氣?固結合強度較大時, Waldmann公式給出的熱泳力計算結果誤差較大; 隨著系統溫度的升高或氣?固結合強度的減小, 納米顆粒所受熱泳力的分子動力學計算結果逐漸收斂于Wald?mann 公式的預測值.這是由氣體分子與納米顆粒間非剛體碰撞效應的強弱受溫度或氣?固結合強度影響所導致的.當結合強度較大或溫度較低時, 勢能在氣體分子與納米顆粒進行動量交換的運動軌跡中占主導地位; 而當結合強度較小或溫度較高時, 氣?固相互作用勢能的影響十分微弱, 剛體碰撞模型近似成立.

圖5 不同參數下熱泳力的分子動力學結果FT,MD 與Waldmann 公式結果FT,E 比較圖 (a) 環境溫度T*; (b) 氣?固結合強度kij; (c) 溫度梯度?T?Fig.5.The variation of thermophoretic force FT between present MD simulation result FT, MD and Waldmann equa?tion result FT, E under different parameters: (a) The environ?ment temperature ; (b) the intensity of gas?particle in?teraction kij; (c) temperature gradient ? T?.

4.2 顆粒尺寸影響

顆粒尺寸的大小不僅影響著顆粒在高低溫兩側的動量交換, 還影響著顆粒與氣體分子之間的碰撞模型, 是影響顆粒所受熱泳力的重要因素之一.為了便于比較, 定義無量綱熱泳力:

式中, FT為熱泳力的理論計算結果, FT= FT,E為基于有效熱導率的Waldmann 理論計算結果, FT=FT,LW為考慮非剛體碰撞效應之后(2)式的計算結果, 對應的無量綱熱泳力分別表示為和圖6 展示無量綱熱泳力和隨顆粒直徑的變化圖.模擬過程中系統的參數為:= 55,ΔT*= 0.92, n*=0.0054.顆粒粒徑變化范圍為 DP?≈ 2.5 到≈9.98, 且模擬系統處于自由分子區.由圖6(a)可以看出, 隨著顆粒尺寸逐漸接近分子尺寸(= 1),無量綱熱泳力單調增加, 而對于顆粒尺寸較大的納米顆粒, 顆粒所受熱泳力的分子動力學結果則逐漸收斂于Waldmann 公式的預測值(由于計算能力的限制, 并沒有給出更大尺寸顆粒的計算結果).對于納米顆粒而言, 由于納米顆粒與氣體分子之間的非剛體碰撞效應, 熱泳力結果并不能收斂于經典的Waldmann 理論值.(2)式對Waldmann理論進行了修正, 修正后的無量綱熱泳力隨顆粒尺寸變化結果如圖6(b)所示.相對于Wald?mann 理論的計算結果, 無量綱熱泳力計算誤差明星降低, 這驗證了(2)式的準確性.但是, 不同氣?固結合強度下的計算結果還存在一定的誤差,這是由于當氣?固結合強度較大時, 氣體分子的動能不足以克服氣?固結合勢能的束縛, 氣體分子會被吸附于納米顆粒表面, 從而引起納米顆粒尺寸的增加, 造成理論與模擬結果的偏差.

圖6 無量綱熱泳力 隨顆粒直徑 變化圖 (a)剛體碰撞模型; (b) 非剛體碰撞模型FT = FT, LW1; (c) 非剛體碰撞模型FT = FT, LW2(考慮氣體分子吸附引起的顆粒尺寸變化)Fig.6.The dimensionless thermophoretic force as a function of for different gas?particle interaction intensity:(a) The rigid body collision model FT = FT, E; (b) the non?rigid body collision model FT = FT, LW1; (c) FT = FT, LW2(wherein the vary of particle size is considered).

圖7 對顆粒表面氣?固相互作用勢能以及氣體分子動能分布進行了統計, 同時計算了顆粒表面氣體分子密度分布函數g.圖7 中, 黑色與藍色虛線分別為顆粒表面第一層與第二層原子的位置, r'=.由圖7 可以看出, 當氣?固結合強度較小時(圖7(a)), 勢井深度較淺, 氣體分子所攜帶的動能足以克服勢能的束縛.此時, 納米顆粒表面的氣體分子數密度與環境中的氣體分子數密度相當, 并沒有發生吸附現象.當氣?固結合強度較大時(圖7(b)), 顆粒表面第一層內的勢井深度較深,而氣體分子的動能相對太低, 使得該層氣體分子無法擺脫勢能的束縛, 從而被吸附在納米顆粒表面,這就造成了納米顆粒表面的第一層內的氣體分子數密度遠遠高于環境中的氣體分子數密度.通過對納米顆粒的粒徑進行修正, 可以得到修正后的等效粒徑, 將等效粒徑代入 (2)式, 無量綱后的結果表示為如圖6(c)為無量綱熱泳力隨顆粒尺寸的變化結果.結果表明, 計算結果與分子動力學模擬結果吻合較好.

圖7 納米顆粒表面氣體原子的勢能Ep, 動能Ek 及密度函數g 分布圖Fig.7.The potential energy Ep, kinetic energy Ek and dens?ity function g on the surface of nanoparticles.

5 結 論

本文研究了自由分子區內納米顆粒的熱泳特性, 基于非平衡態分子動力學模擬方法計算了納米顆粒所受熱泳力.考慮有限空間效應對氣體熱物性的影響, 引入氣體的有效熱導率計算得到了較為準確的Waldmann 自由分子極限, 并與分子動力學模擬結果進行了對比.結果表明, 對于納米顆粒來說, 當氣?固相互作用勢能較弱或氣體溫度較高時,Waldmann 熱泳力公式與分子動力學模擬結果吻合較好; 當氣?固相互作用勢能較強或氣體溫度較低時, 由于氣體分子與納米顆粒之間的非剛體碰撞效應較為明顯, 基于剛體碰撞模型假設的Wald?mann 計算式與模擬結果存在較大誤差.基于Li?Wang 計算式, 并考慮由于納米顆粒對氣體分子的吸附效應引起顆粒等效尺寸的變化, 計算結果與分子動力學模擬結果較為吻合.