馬奧杰 陳頌佳 李玉秀 陳穎

(廣東工業(yè)大學材料與能源學院, 廣州 510006)

介觀尺度下顆粒布朗運動的摩擦系數符合黏性流體力學邊界條件, 當顆粒尺度減小至納米級別時, 邊界條件向滑移過渡; 另一方面, 隨著顆粒尺度的減小, 顆粒表面溶劑分子的吸附效應對顆粒水動力學半徑的影響不可忽略. 分子動力學模擬可以捕獲納米流體中顆粒與溶劑分子相互作用的微觀細節(jié)且計算精度高. 以剛性TIP4P/2005水分子模型為溶劑, 建立不同大小的Cu納米顆粒在水中擴散的全原子模型. 采用單顆粒追蹤方法對顆粒的平動和轉動擴散系數進行擬合, 將摩擦因子與黏性邊界條件和滑移邊界條件下的結果進行比較, 并研究了溶劑在顆粒表面的吸附特性. 研究發(fā)現納米顆粒的平轉動摩擦因子均在兩種邊界條件所預測的理論值之間, 顆粒尺寸越小, 溶劑分子的吸附越明顯, 顆粒表面的水分子層會增大顆粒的水動力學半徑使得摩擦因子的計算結果偏向黏性邊界, 納米顆粒尺寸約小于5倍溶劑分子尺寸時, 需要考慮顆粒表面的溶劑層對顆粒水動力半徑的影響.

1 引 言

布朗運動是指懸浮在流體中的納米至微米級顆粒所經歷的無休止的隨機運動, 在自然界中廣泛存在. 關于布朗運動的起源在很長的一段時間都備受爭論[1], 直到1905年Einstein[2]根據分子運動論原理提出了布朗動力學理論, 并建立了顆粒的擴散系數和摩擦系數的關系; Paul Langevin將隨機力和黏性阻力引入牛頓第二定律[3], 由此得到的Langevin方程以一種簡潔的方式再現了Einstein的結果. 一百多年來, 布朗運動的后續(xù)研究及其推廣對物理學、數學、生物學和材料科學等領域[4,5]的發(fā)展產生了深遠的影響.

Einstein的布朗動力學理論指出顆粒的隨機運動是周圍流體分子熱運動的結果. 顆粒在周圍流體分子的碰撞下產生隨機運動的同時, 又受到流體分子的阻力, 由碰撞所產生的隨機力和流體分子對顆粒的摩擦阻力通過漲落耗散定理約束. 顆粒的熱漲落和摩擦力的關系為D=kBT/f,D為顆粒的擴散系數,kB,T,f分別為玻爾茲曼常數、體系絕對溫度、顆粒的摩擦系數. 在黏性水動力學邊界條件下, 顆粒的平動摩擦系數ft=6πηa, 轉動摩擦系數fr=8πηa3,η為溶劑的動力學黏度,a為顆粒半徑.然而, 當研究顆粒的尺度由微米級降低至納米級別時, 相關學者[6]發(fā)現采用黏性邊界條件計算得出的旋轉弛豫時間比測量值大5—10倍. 顯然當研究的問題涉及微觀領域, 黏性流體力學邊界條件理論具有明顯的局限性. 如圖1(a)所示, 黏性流體力學邊界條件下, 在流體中運動的球形顆粒攜帶著與其相接觸的流體層而不發(fā)生滑移, 顆粒摩擦力來自兩個部分: 顆粒正前方流體分子阻礙顆粒運動而產生的阻力和由于流體的黏性而產生的平行于顆粒運動方向的粘滯力; 滑移邊界條件下(圖1(b)), 摩擦力僅來自于顆粒正前方流體分子所產生的的阻力[7].顆粒的旋轉摩擦系數取決于繞其對稱軸旋轉時對周圍流體產生的擾動, 黏性邊界條件下, 球形顆粒會攜帶周圍的流體旋轉, 摩擦系數為 8 πηa3; 而在滑移邊界下, 顆粒表面與流體分子的相互作用被忽略, 沒有流體分子阻礙顆粒的旋轉, 因此旋轉摩擦系數為0. 無論是黏性邊界還是滑移邊界條件,都對流體分子與顆粒的摩擦系數作了理想假設.Richardson[8]的研究表明, 在一個足夠粗糙的表面由滑移所產生的速度場和等效光滑表面由黏性產生的速度場相近, 這表明顆粒表面粗糙程度會影響流體與顆粒的相互作用, 而且隨著納米顆粒尺寸的降低, 顆粒與溶劑相互作用對顆粒有效水動力學半徑的影響不可忽略[9]. 顆粒表面特性和電解質溶液的濃度都會影響溶劑分子在顆粒表面的吸附從而改變顆粒的有效半徑[10,11], Velasco-Velez等[12]用X射線吸收方法和分子模擬方法研究了水分子在金顆粒表面的結構, 并指出顆粒表面懸空的氫鍵是水分子產生吸附的原因. Vasanthi等[13,14]以及何昱辰和劉向軍[15]采用粗粒化模型模擬了納米顆粒的平轉動擴散特性, 計算所得的擴散系數與滑移邊界條件下的結果相近, 然而粗?；Ｐ筒⒉荒芊从愁w粒表面特性和溶劑分子的真實結構. 在介觀尺度的模擬方面, 耗散動力學模型和布朗動力學模型將隨機力和黏性力引入顆粒的控制方程來約束顆粒的運動, 解決了由于顆粒與溶劑分子在時間和空間尺度上的差異造成計算成本過高的問題[16], 但是由于采用粗粒化模型或者隱式溶劑模型, 忽略了溶劑分子與顆粒相互作用的微觀細節(jié).

圖1 球形顆粒在不同邊界條件下的平動摩擦系數, 顆粒向右運動 (a)黏性邊界條件; (b)光滑邊界條件Fig. 1. Translational friction coefficients of spherical particles under different boundary conditions and the particle move to the right: (a) Translational friction coefficient for stick boundary; (b) translational friction coefficient for slip boundary.

本文采用分子動力學模擬方法建立了單個Cu納米顆粒在水溶液中擴散的全原子模型, 在原子水平上研究流體分子與納米顆粒之間的相互作用, 在計算能力能夠承受的范圍內分析比較了4種不同半徑的納米顆粒在水溶液中的擴散特性. 采用單顆粒追蹤算法[17-19]并結合最小二乘法對顆粒的平轉動擴散系數進行擬合, 解決了擴散系數需要對大量顆粒系綜平均的問題. 將摩擦系數的計算結果與黏性邊界和滑移邊界條件下的結果進行比較, 討論了顆粒表面的溶劑吸附對顆粒水動力學半徑的影響.膠體系統的粗?；M過程中, 摩擦系數的選擇直接決定布朗力和黏滯力的大小, 本文的結果可以為摩擦因子的選擇提供支持, 而不必在黏性邊界與滑移邊界條件之間做出選擇, 為布朗動力學和耗散動力學模型的完善提供理論依據.

2 模型建立與驗證

2.1 全原子模型的構建

本文建立了不同大小的單個Cu納米顆粒在純水溶液中擴散的全原子模型(圖2(a)), Cu顆粒半徑分別為0.5, 1.0, 1.5和2.0 nm, 顆粒中原子按照FCC的晶格結構排列. 初始時刻顆粒位于模擬域的中心, 模擬域的大小和總原子數目需要根據顆粒的大小進行調整, 如表1所列. 邊界條件均采用周期性邊界條件, 從而保證顆粒的運動為自由布朗運動. 模擬開始時顆粒的初始速度符合麥克斯韋分布, 模擬時間步長為1 fs, 系統溫度均為298 K. 由于溶劑動力學黏度對顆粒摩擦系數的計算有直接的影響, 而不同水分子模型的動力學黏度相差較大, 本文中水分子的構建選用與實驗符合較好的剛性TIP4P/2005模型[20]. 除了兩個氫原子和一個氧原子, 該模型在H—O—H的鍵角平分線上增加了一個無質量的作用位點M, 水分子結構如圖2(b)所示,其中鍵長b= 0.9572 ?, 鍵角θ= 104.52°, O和H的電荷分別為—1.1128和+0.5564 eV, 水分子的半徑d= 2 ?, OM = 0.1546 ?. 水分子的鍵長和鍵角采用Shake算法固定, 因此鍵合勢為零, 長程靜電力采用PPPM算法處理[21]. 原子間的非鍵合相互作用包括范德瓦耳斯相互作用和靜電相互作用, 計算公式如下:

表1 各工況下模擬域大小與原子數目Table 1. The size of the simulation domain and number of atoms of different working conditions.

圖2 (a) Cu-H2O納米流體的物理模型; (b) TIP4P/2005水分子結構示意圖Fig. 2. (a) Physical model of Cu-H2O nanofluid; (b) schematic diagram of TIP4P/2005 water.

其中Uvdw代表原子間的短程引力勢,Uelectronstatic代表原子間的長程靜電勢,ε和σ分別代表原子的勢能參數和尺寸參數,qi和qj為水分子中不同原子的帶電量,rij代表原子間的距離.

顆粒與水分子的相互作用只考慮Cu原子與氧原子的范德瓦耳斯作用, Cu-H2O納米流體相互作用的勢能參數如表2所列, 計算過程不考慮氫原子和氧原子及銅原子的范德瓦耳斯作用. 顆粒中原子間的相互作用采用適用于金屬和合金的多體勢(EAM)[22], 這種EAM勢在描述金屬原子間的相互作用時比對勢更精確, 晶體的總勢能表示為

表2 Cu-H2O納米流體勢能參數Table 2. Potential energy parameter of Cu-H2O nanofluid.

其中,Fi為嵌入能, 是除了第ρi個原子之外的所有其他原子的核外電子在第i個原子處所產生的電子云密度之和,φij為對勢作用項,rij為原子i和原子j間的距離.

系統能量最小化后, 在NPT系綜下弛豫平衡,采用Nose-Hoover方法控溫控壓. 目標溫度和壓力分別為298 K和 1 bar (1 bar = 105Pa), 時間步長為1 fs, 經過20 ps的弛豫后, 壓力和溫度趨于穩(wěn)定. 體系穩(wěn)定后, 在NVE系綜下使顆粒做自由布朗運動, 總步數為6000 萬步, 對應總時長為60 ns.上述納米流體模型的建立以及模型驗證工作均采用LAMMPS軟件包完成.

2.2 水分子模型的驗證

為了得到穩(wěn)定的水分子體系, 需要確保系統在模擬過程中有足夠的時間達到平衡. 選取隨機分布的4096個水分子, 大小為49.6 ? × 49.6 ? × 49.6 ?的模擬域, 對弛豫過程中體系的溫度和總能進行監(jiān)測, 溫度和總能隨平衡時間的變化如圖3(a)和圖3(b)所示. 可以看出, 體系溫度和總能量在300 fs左右時已經達到平衡, 總體積和密度亦趨于穩(wěn)定.由于水分子的動力學黏度對摩擦系數的計算至關重要, 本文采用green-kubo公式(5)式計算水的黏度, 該方法將流體的剪切黏度和壓力張量中非對角元素的自相關函數聯系起來, 通過計算壓力張量的非主對角元素pxy,pxz,pyz的平均值得到TIP4P/2005剛性水分子模型的剪切黏度[23],

圖3 水分子模型的驗證 (a) NPT系綜下體系的溫度隨時間變化; (b) NPT系綜下體系的總能量隨時間變化Fig. 3. Validation of water molecular model: (a) Curve of temperature over time; (b) total energy of the system over time.

其中η為水分子的動力學黏度,V為計算域的體積,kB為玻爾茲曼常數,Pij為壓力張量的非主對角元素.

NPT系綜下平衡之后, 在NVT系綜下計算體系的黏度, 經過2 ns的模擬, 得到水分子的黏度為0.825 × 10—3Pa·s, 該黏度與實驗結果的誤差約為7%[24], 其他水分子模型(TIP3P, TIP4P, SPE/C)誤差與實驗結果在18%—64%之間[23], 因此TIP4P/2005水分子模型的動力黏度已明顯優(yōu)于其他水分子模型.

3 結果與討論

3.1 單顆粒布朗運動平動擴散

在統計力學中, 均方位移是指顆粒隨時間運動的距離相對于參考位置的度量, Einstein在布朗運動理論中指出顆粒移動距離的平方的均值與時間呈線性關系:

其中 [ Δx(t)]2為顆粒當前位置沿x方向相對于參考位置距離的平方;〈〉 表示體系中所有顆粒的均值;D為顆粒的擴散系數, 反映了顆粒在溶劑中擴散能力的大小.

由于布朗運動的方向具有隨機性, 它的軌跡是一條由在空間隨機分布的點依次連接所構成的曲線, 該曲線處處連續(xù)且處處不可導[25]. 隨著顆粒尺寸的增加, 體系內的總原子數因模擬域的增大而成倍增長, 受到計算資源的限制, 所有體系中只放置一個顆粒. 顆粒在總模擬時長內軌跡如圖4所示,藍色代表軌跡的起點, 紅色表示軌跡的終點, 模擬溫度T= 298 K, 顆粒半徑a= 1 nm. 從圖4可以看出, 顆粒的運動方向并沒有發(fā)生明顯的變化, 但是軌跡同樣滿足處處連續(xù)且處處不可導的特征, 這是因為顆粒在運動過程中并未與其他顆粒發(fā)生碰撞. 對于單顆粒布朗運動軌跡, 采用傳統的計算方法得到的擴散系數隨計算時長的增加而增加, 為了準確地計算單個顆粒的均方位移, 本文采用單顆粒追蹤方法并結合最小二乘法對顆粒的平動擴散系數進行擬合. 下面以半徑a= 1 nm的顆粒為例說明擴散系數的擬合過程.

圖4 單顆粒布朗運動軌跡Fig. 4. Brownian motion trajectory of a single particle.

1)確定最佳關聯時間. 顆粒的軌跡為含有N=60000個數據點的連續(xù)曲線, 曲線的每一段的時間間隔為1 ps. 將該軌跡按不同的關聯時間(時間間隔)Nt分成Ns=N/Nt份, 對每一段軌跡, 根據下式計算顆粒的均方位移:

其中Nt為每一段軌跡的關聯時間,ri為顆粒的位置坐標.

根據所求均方位移進而確定擴散系數, 然后對每一個關聯時間下的Ns個擴散系數進行直方圖分析. 如圖5所示, 當關聯時間較小時, 擴散系數呈正偏態(tài)分布, 關聯時間很大時, 擴散系數呈負偏態(tài)分布, 一個呈正態(tài)分布的關聯時間即是最佳關聯時間[26].

圖5 半徑a = 1 nm的Cu顆粒在不同關聯時間下擴散系數分布Fig. 5. Diffusion coefficient distribution of Cu particles with radius a =1 nm at different correlation time.

2)確定最佳擬合點. 從(7)式可以看出, MSD曲線的第一個點代表每一段位移的平均, 而最后一個數據點沒有被平均. 計算擴散系數時統計點過多或者過少都會對計算結果帶來不利的影響, 因此本文計算了每一段軌跡采用不同擬合點時擴散系數D的變異系數(二階矩/一階矩), 變異系數最低的擬合點數目即為最佳擬合點[27]. 圖6為關聯時間t= 1363 ps時得到的變異系數隨擬合點變化的曲線, 可以看出, 擬合點數目在20左右時變異系數最低.

圖6 半徑a = 1 nm的Cu納米顆粒在采用不同擬合點擬合時擴散系數的變異系數, 關聯時間Nt = 1363 psFig. 6. Variation coefficient of the diffusion coefficient of Cu nanoparticles with radius a =1 nm when different fitting points were used for fitting, correlation time Nt = 1363 ps.

將最佳擬合點下的擴散系數取平均可以得到半徑a= 1 nm的Cu顆粒的擴散系數, 采用同樣的方法可求得a= 0.5, 1.5, 2.0 nm時顆粒的擴散系數. 圖7(a)表示4種顆粒的平動擴散系數, 紅色曲線和藍色曲線分別表示黏性邊界條件和滑移邊界條件下的預測結果, 從圖7(a)可以看出, 平動擴散系數均在預測值之間. 平動擴散的摩擦因子可由α=kBT/Dtπηa求得, 如圖7(b)所示, 當半徑a=1.0, 1.5, 2.0 nm時, 摩擦因子更接近滑移邊界條件; 半徑a= 0.5 nm時, 摩擦因子更接近黏性邊界條件. 隨著顆粒尺寸的減小, 摩擦因子應逐漸趨近于4, 但對于顆粒半徑為0.5和1.0 nm的顆粒, 摩擦因子有逐漸增大的趨勢, 這是因為吸附在顆粒上的水分子增大了顆粒的水動力學半徑, 導致計算的摩擦因子偏向黏性邊界.

圖7 納米顆粒的平動擴散特性 (a)不同尺寸Cu納米顆粒的平動擴散系數; (b)不同尺寸Cu納米顆粒的平動摩擦因子Fig. 7. Translational diffusion characteristics of nanoparticles: (a) Translational diffusion coefficients of Cu nanoparticles with different sizes; (b) translational friction factors of Cu nanoparticles of different sizes.

3.2 單顆粒布朗運動的旋轉擴散

與平動擴散均方位移的定義類似, 顆粒的均方旋轉角度可按下式定義:

其中, MSR表示顆粒的均方旋轉角度, Δθ(t) 代表顆粒從初始時刻到t時刻繞某一旋轉軸轉過的角度,Dθ即為顆粒的旋轉擴散系數. 如圖8所示,顆粒在一個時間間隔內轉過的角度為O-O軸與O'-O'軸的夾角, 由于球形顆粒的旋轉具有各向同性, 不計顆粒繞O-O軸的旋轉對計算結果沒有影響. 四元數可以很方便地表征空間向量在三維空間中的旋轉, 本工作通過顆粒在每一個位置的四元數來計算顆粒旋轉的角度. 假設顆粒在初始時刻的單位四元數q= (a1,a2,a3,a4), 將其轉化為旋轉矩陣R, 如(9)式所示. 在球面上任取一單位向量u,將該單位向量右乘旋轉矩陣R即可得到初始時刻的單位向量u1, 同樣的方法可求得向量u在其他時刻的方向向量.

圖8 顆粒旋轉擴散系數計算的物理模型Fig. 8. Physical model for particle rotation diffusion coefficient calculation.

不同于平動擴散系數, 旋轉擴散系數的擬合對關聯時間的選擇不敏感, 因此可以直接對顆粒的MSR曲線進行擬合. 顆粒的均方旋轉角度如圖9所示, MSR曲線隨時間單調遞增, 對該曲線進行線性擬合即得到顆粒的旋轉擴散系數. 圖10(a)表示顆粒的旋轉擴散系數隨顆粒尺寸的變化, 紅色曲線代表由Stokes-Einstein關系得到的顆粒在黏性邊界條件下的擴散系數, 從圖10(a)可以看出, 隨顆粒尺寸的增加, 擴散系數逐漸趨近于黏性邊界調下的結果. 圖10(b)給出了旋轉摩擦因子的大小, 旋轉摩擦因子β可由β=kBT/πηa3求得. 從圖10(b)可以看出, 旋轉摩擦因子β均在黏性邊界與滑移邊界所規(guī)定的范圍內, 介于4—8之間, 因此更接近黏性邊界條件. 顆粒尺寸越大, 邊界條件越接近黏性邊界, 但半徑a= 0.5 nm和a= 1.0 nm的顆粒, 摩擦因子相比a= 1.5 nm和a= 2.0 nm的顆粒反而有所升高, 這是顆粒表面溶劑分子吸附的結果.

圖9 不同半徑的Cu納米顆粒的均方旋轉角度Fig. 9. Mean-square rotation angles of Cu nanoparticles with different sizes.

圖10 顆粒旋轉擴散特性 (a)顆粒的轉動擴散系數; (b) 顆粒轉動擴散摩擦因子Fig. 10. Rotational diffusion characteristics of particles: (a) Rotational diffusion coefficients of particle; (b) rotational friction factors of particle.

3.3 納米顆粒表面溶劑分子的吸附特性

在滑移邊界條件假設下, 納米顆粒運動過程中表面不會攜帶流體分子. 然而考慮到納米顆粒的粗糙度和固-液相互作用, 顆粒周圍往往會吸附溶劑分子, 這部分溶劑分子與納米顆粒一起做隨機運動, 因此吸附的水分子層會對顆粒的有效水動力學半徑產生影響. 本文統計了納米顆粒周圍水分子的分布[28], 如圖11(a)所示, 通過統計距離納米顆粒表面為r, 厚度為dr的球殼內的原子數目, 得到了水分子的徑向分布函數(RDF), RDF可用下式來計算:

圖11 納米顆粒周圍水分子的RDF (a) RDF計算的物理模型; (b)不同尺寸的顆粒周圍水分子的RDF曲線Fig. 11. Radial distribution function (RDF) of water molecules around nanoparticles: (a) Physical model for RDF calculation; (b) RDF curves of water molecules around particles of different sizes.

其中n(r) 表示落在厚度為 dr的球殼內的原子數目,Vr為球殼的體積,ρ0表示無窮遠處水分子體系的原子數密度, 隨著r的增加, 球殼內原子的數密度等于溶劑分子的平均數密度,g(r)趨近于1.

圖11 (b)為納米顆粒周圍原子數目的徑向分布函數. 可以看出, 不同尺寸的納米顆粒周圍原子數分布均在1.7 ? (略小于水分子半徑)左右時達到峰值, 峰值越大, 代表溶劑分子的吸附程度越緊密, 曲線中第一個峰的高度隨著顆粒尺寸的增加依次降低; 第二個峰值出現在距離表面4.3 ?處(略大于一個水分子的直徑), 且峰的高度沒有明顯差別, 之后水分子數密度逐漸趨于標準水分子密度.因此, 在距離顆粒表面一個水分子直徑的區(qū)域內水分子排列呈現結構化, 導致顆粒周圍水分子的密度大于水分子的平均密度, 固體和水分子之間懸空的氫鍵是導致水分子在顆粒表面產生吸附的原因[12],氫鍵存在于水分子中, 但水分子不能和金屬原子之間形成氫鍵, 因此水分子趨于向顆粒表面靠近. 水分子在顆粒表面的吸附增大了顆粒的有效半徑, 而且顆粒尺寸越小, 表面水分子層相對于顆粒的尺寸越大, 對顆粒水動力學半徑的影響越大, 這是導致在轉動摩擦因子的計算中, 尺寸a= 0.5和1.0 nm的顆粒摩擦因子向非滑移邊界移動的原因.

4 結 論

本文采用分子動力學模擬方法建立了不同尺寸的單個Cu納米顆粒在水中擴散的全原子模型,通過單顆粒軌跡分析, 采用最小二乘法擬合了顆粒的平動和轉動擴散系數, 并將摩擦因子與黏性邊界與滑移邊界條件下的結果作對比. 顆粒的平動和轉動擴散系數隨著顆粒尺寸的減小而增加, 且摩擦因子均在黏性邊界條件與滑移邊界條件所規(guī)定的范圍內. 由于水分子之間存在氫鍵網絡而水分子與金屬原子之間無法形成氫鍵, 所以水分子會吸附在金屬表面, 吸附在顆粒上面的水分子會增大顆粒的水動力學半徑. 無論顆粒大小, 水分子的吸附現象始終存在, 對于大的顆粒, 由水分子吸附所增加的半徑遠小于顆粒的真實半徑, 所以對顆粒有效半徑的影響不夠明顯; 對于尺寸較小的顆粒, 由水分子吸附所增加的半徑與顆粒的真實半徑接近, 所以對小顆粒有效半徑的影響較為顯著. 因此對半徑a=0.5, 1.0 nm的顆粒, 其平動和轉動摩擦因子比預測值偏大, 因為顆粒的表面水分子的吸附能力較強, 當顆粒的尺寸小于水分子半徑的5倍時, 水分子吸附對顆粒有效水動力學半徑的影響不可忽略.