過渡區(qū)內(nèi)納米顆粒的曳力特性模擬研究*

劉旺旺 張克學 王軍 夏國棟

(北京工業(yè)大學,傳熱強化與過程節(jié)能教育部重點實驗室,傳熱與能源利用北京重點實驗室,北京 100124)

在過渡區(qū)內(nèi),關(guān)于納米顆粒曳力計算及輸運特性的研究較為困難,通常會采用一些近似方法,將自由分子區(qū)或者連續(xù)介質(zhì)區(qū)的理論計算式進行修正,以適用于過渡區(qū),但是其準確性值得商榷.本文基于分子動力學模擬方法,研究了過渡區(qū)內(nèi)納米顆粒的曳力特性,并與相關(guān)理論進行對比.結(jié)果表明,氣-固分子間相互作用對納米顆粒的曳力具有顯著影響.當氣固結(jié)合強度較弱時,理論計算結(jié)果與分子動力學模擬值吻合較好;當氣固結(jié)合強度較強時,分子動力學模擬結(jié)果明顯大于理論值,這是由于氣體分子在納米顆粒表面的吸附所導致.基于氣體分子在顆粒表面的吸附特性,提出引入有效顆粒半徑修正,其過渡區(qū)內(nèi)曳力的理論計算結(jié)果與分子動力學模擬結(jié)果吻合較好.

1 引言

當顆粒以一定的相對速度在氣體中運動時,氣體會阻礙顆粒的運動,這個阻礙顆粒運動的力被稱為曳力[1-3].關(guān)于顆粒曳力特性的研究在很多領(lǐng)域都有重要意義,例如氣溶膠科學、環(huán)境科學、燃燒學、納米材料合成以及核安全等[4-11].曳力通常被認為是懸浮顆粒在流體中受到的最主要的作用力,根據(jù)Stokes-Einstein 輸運理論,顆粒的布朗擴散系數(shù)與其曳力系數(shù)成反比.近年來,關(guān)于納米顆粒在氣體中曳力計算及輸運特性的研究引起了廣泛的關(guān)注[12-18].

懸浮顆粒在氣體中的受力及其輸運特性與Knudsen 數(shù)密切相關(guān).Knudsen 數(shù)定義為氣體分子平均自由程λ 與顆粒半徑R的比值,Kn=λ/R.根據(jù)Knudsen 數(shù)的大小,可將顆粒的動力學行為劃分為自由分子區(qū)(Kn?1)、連續(xù)介質(zhì)區(qū)(Kn?1)和過渡區(qū)(Kn～1).在自由分子區(qū),氣體分子平均自由程遠大于顆粒半徑,因此可忽略被顆粒反射的氣體分子與其他入射分子之間的碰撞,即不考慮顆粒的存在對氣體分子速度分布函數(shù)的影響.基于氣體動理論方法,Epstein[19]推導了自由分子區(qū)內(nèi)球形顆粒的曳力.在連續(xù)介質(zhì)區(qū),氣體分子之間的碰撞頻率遠大于氣體分子與顆粒之間的碰撞頻率,因此,氣體介質(zhì)可視為連續(xù)流體.Stokes[20]提出了球形顆粒在連續(xù)介質(zhì)區(qū)內(nèi)曳力的表達式.在過渡區(qū),氣體分子間的碰撞頻率和顆粒與氣體分子之間的碰撞頻率相當,連續(xù)介質(zhì)假設(shè)不再適用,自由分子區(qū)的處理方法同樣也不再適用,因此對曳力的求解比較困難.通常會采用一些近似方法建立理論模型,如矩方法和Bhatnagar-Gross-Krook 方法等,但其準確性值得商榷[21].通過引入修正系數(shù),可以將連續(xù)介質(zhì)區(qū)或者自由分子區(qū)的理論結(jié)果拓展,以適用于過渡區(qū),例如Stokes-Cunningham 曳力修正公式[22]:

式中,μ為氣體的黏度,V為顆粒相對于氣體的速度,A,B,E為經(jīng)驗參數(shù).實際中常采用A=1.155,B=0.471,E=0.596[23].

上述理論均基于氣體分子與顆粒之間的剛體碰撞假設(shè)建立模型,此假設(shè)對于微米量級的大顆粒來說是適用的,但是當顆粒的尺寸逐漸降低到亞微米或納米量級時,剛體碰撞模型可能不再適用[24-26].此時需要考慮氣體分子與顆粒之間的非剛體碰撞效應(yīng),例如范德瓦耳斯力等[27-29].納米顆粒與氣體分子之間的非剛體相互作用會使氣體分子的運動軌跡發(fā)生偏轉(zhuǎn),從而影響碰撞過程中的動量傳遞,最終影響顆粒的受力特性.Li 和Wang[30]基于氣體動理論和非剛體碰撞模型,得到自由分子區(qū)內(nèi)納米顆粒所受的曳力表達式,并采用半經(jīng)驗的處理提出了過渡區(qū)內(nèi)納米顆粒的曳力表達式[31]:

式中,w是由實驗數(shù)據(jù)得到的擬合值,取為1.143[31],mr是約化質(zhì)量,kB是玻爾茲曼常數(shù),T是系統(tǒng)溫度,n是氣體數(shù)密度,是無量綱碰撞積分項.但是,(2)式的準確性尚未得到驗證.

基于分子動力學(molecular dynamics,MD)模擬方法,研究了過渡區(qū)內(nèi)納米顆粒的曳力特性.研究發(fā)現(xiàn),對于納米顆粒而言,模擬結(jié)果與(2)式并不相符.氣固結(jié)合強度對納米顆粒的曳力有顯著影響,當結(jié)合強度較強時,需要考慮氣體分子在顆粒表面的吸附效應(yīng).使用修正的有效顆粒半徑后計算(2)式的結(jié)果與分子動力學模擬結(jié)果吻合良好.

2 分子動力學模擬計算

2.1 分子動力學模型

本文采用分子動力學模擬方法計算納米顆粒在過渡區(qū)內(nèi)所受的曳力.分子動力學模型如圖1 所示,納米顆粒懸浮于一定量的氣體分子中,其中納米顆粒由NS個銅原子組成,氣體分子由NF個氬原子組成.模擬系統(tǒng)在x,y,z方向上的長度分別為Lx=270 ?,Ly=80 ?,Lz=80 ? (1 ?=10-10m),3 個方向均采用周期性邊界條件.在本文的模擬系統(tǒng)中,對銅原子之間的相互作用采用嵌入原子法(embedded atomic method,EAM)勢能[32]:

圖1 分子動力學模擬系統(tǒng)模型圖(紅色: 納米顆粒;藍色: 氣體分子)Fig.1.Configuration of MD simulation system (Red: nanoparticle;blue: gas molecules).

其中Ei,ρi,re分別代表嵌入能量、原子電子密度、最近鄰之間的平衡間距,P,ξ,Q,ψ 是4 個可調(diào)參數(shù),ω 和τ 是兩個額外的截止參數(shù).對氣-氣和氣-固原子之間的相互作用采用經(jīng)典的Lennard-Jones(LJ 12-6)勢能函數(shù):

其中,ε 表示LJ 勢阱深度,σ 是勢能為零時的距離,rij為原子i和j之間的距離.氣體與固體原子之間勢能參數(shù)的選取,遵循Lorentz-Berthelot 混合原則:

表1 列出了LJ 勢能函數(shù)的參數(shù)[33].本文通過改變氣-固分子之間勢函數(shù)的能量參數(shù)來改變氣固界面的結(jié)合強度:

表1 LJ 勢函數(shù)參數(shù)值[33]Table 1.LJ potential parameters[33].

本文使用LAMMPS (large-scale atomic/molecular massively parallel simulator)開源軟件進行分子動力學模擬計算[34].在模擬過程中,時間步長取為2 fs,截斷半徑rc=0.9 nm.在弛豫階段,系統(tǒng)首先在NVT 系綜下運行3 ns,使用Nosè-Hoover 熱浴使得系統(tǒng)的溫度保持在180 K.然后將模擬系統(tǒng)在x方向分為兩個區(qū)域,左端長度為2.7 nm 的區(qū)域設(shè)定為熱浴區(qū),右端長度為24.3 nm 的區(qū)域設(shè)定為外力區(qū),如圖1 所示.在熱浴區(qū)施加NVT 系綜使其維持在180 K,在外力區(qū)域給每個氣體分子添加一個外力,5 ns 后系統(tǒng)內(nèi)可以形成恒定的流速.最后在NVE 系綜下運行120 ns并計算相關(guān)物理量.基于遍歷原理,本文所給出的任何物理量都是在采樣時間范圍進行的長時間平均得到的.為了提高計算結(jié)果的準確性,本文采用不同初始微觀條件(即原子初始速度不同,但需滿足同一個Maxwell 速度分布律)下至少8 次獨立運行結(jié)果的算數(shù)平均值作為最終計算結(jié)果.

2.2 計算方法

對于懸浮在氣體中的納米顆粒,其布朗運動會給曳力的計算帶來較大的困難.因此,本文在納米顆粒上施加了一個額外的簡諧回復勢能Uhar[35-37]:

式中,khar表示簡諧勢能的彈性系數(shù),xc和xo分別表示納米顆粒的質(zhì)心位置和初始位置.當系統(tǒng)達到穩(wěn)定后,納米顆粒所受的曳力與簡諧力平衡,納米顆粒相對于它的初始位置會具有一個相對位移Δx,此時,納米顆粒所受的曳力FD可以通過下式計算:

為研究外加簡諧勢彈性系數(shù)對納米顆粒所受曳力的影響,計算了曳力隨彈性系數(shù)khar的變化情況.圖2 計算結(jié)果表明,彈性系數(shù)的大小對納米顆粒的曳力計算結(jié)果幾乎沒有影響,這是由于本文采用的曳力計算方法是基于顆粒所受曳力與簡諧回復力之間的平衡.當彈性系數(shù)khar增大時,相對距離Δx會相應(yīng)減小,但是二者的乘積基本保持不變.因此,曳力不會受到彈性系數(shù)khar的影響,而是依賴于顆粒與流體的物性及它們之間的相互作用.文獻[38]的研究結(jié)果表明,只要khar＞0.1εF/,彈性系數(shù)對統(tǒng)計結(jié)果不會產(chǎn)生影響,這與本文的結(jié)論一致.本文選取khar=30εF/.

圖2 納米顆粒的曳力隨彈性系數(shù)khar 的變化Fig.2.Drag force on nanoparticles versus elastic coefficient khar.

本文采用Green-Kubo 公式對氣體黏度進行計算[39]:

式中,V為模擬盒子的體積,〈···〉為系綜平均值,Jγν為γν應(yīng)力張量的分量 (γν=xy,xz,yz),N,m,v分別為氣體分子的數(shù)量、質(zhì)量、速度.

球形納米顆粒的半徑可使用下式進行計算[40]:

式中,xj為固體原子j的坐標.

3 納米顆粒的曳力計算

3.1 納米顆粒的有效粒徑

在過渡區(qū)內(nèi),納米顆粒的曳力受Knudsen 數(shù)的影響較大.通過選取適當?shù)臍怏w分子數(shù)來改變系統(tǒng)的Knudsen 數(shù),使得模擬區(qū)域處于過渡區(qū)內(nèi).圖3 顯示了3 種氣固結(jié)合強度下曳力隨著Kn的變化.可以看到,納米顆粒的曳力隨著Kn而逐漸降低.這是由于當Kn增大時,氣體分子數(shù)逐漸減少,氣體分子與顆粒之間的碰撞頻率逐漸降低,氣固之間的動量傳遞隨之減弱,因而納米顆粒受到的曳力逐漸減小.此外,當氣固結(jié)合強度較大時,氣體分子與顆粒之間的相互作用較強,氣固之間的動量傳遞較強,因而顆粒會受到更大的曳力.

圖3 MD 模擬值與理論公式在3 種氣固結(jié)合強度下的對比結(jié)果 (a) kSF=0.01;(b) kSF=0.1;(c) kSF=0.5Fig.3.Comparison of MD simulation results of the drag force with theoretical values at various kSF: (a) kSF=0.01;(b) kSF=0.1;(c) kSF=0.5.

為了便于將模擬結(jié)果與理論公式進行比較,將(1)式和(2)式的計算結(jié)果也繪制在圖3 中.當kSF=0.01 時,氣體分子與納米顆粒之間的碰撞過程近似于剛體碰撞,此時MD 模擬結(jié)果與(1)式和(2)式的結(jié)果均吻合良好.當kSF=0.1 和kSF=0.5 時,氣體分子與納米顆粒之間的非剛體碰撞效應(yīng)不可忽略,相比于(1)式,MD 模擬值更接近于(2)式的結(jié)果,這表明考慮氣-固間非剛體相互作用的(2)式對納米顆粒的曳力具有更好的預測.此外,當kSF=0.1 時,MD 模擬結(jié)果略大于(2)式的預測值,而當kSF=0.5 時,MD 模擬結(jié)果明顯大于(2)式的結(jié)果(見空心圓點).由于氣體分子與納米顆粒之間的相互作用,氣體分子通常會被吸附在顆粒表面,并在顆粒表面隨機行走,當氣體分子獲得足夠多的能量后,就會從顆粒表面離開,在顆粒表面達到一個吸附與脫附的動態(tài)平衡[41,42],如圖4 所示.此外,氣體分子在顆粒表面的吸附特性可由一個無量綱參數(shù)來表征,即/(kBT),其中=kSFεSF.該比率體現(xiàn)了氣固結(jié)合能與氣體分子動能之間的相對大小[42-44].當/(kBT) 較小時,氣固結(jié)合能低于氣體分子的動能,氣體分子容易克服界面勢壘而脫附.當/(kBT)逐漸增大時,氣固之間的結(jié)合能逐漸增大或者氣體分子的動能逐漸降低,因而大量氣體分子被吸附于顆粒表面.圖4 展示了4 種kSF下氣體分子在納米顆粒表面吸附的模擬快照.當kSF=0.01 時,/(kBT)?1,幾乎沒有氣體分子能被吸附在顆粒表面,如圖4(a)所示,此時氣體分子與顆粒之間的相互作用接近于剛體碰撞,因此MD 模擬結(jié)果與(1)式及(2)式的預測吻合良好;當kSF=0.1 時,/(kBT)略小于1,部分氣體分子可吸附在顆粒表面,見圖4(b);當kSF=0.5 時,/(kBT)＞1,大量氣體分子被吸附,并在顆粒表面形成一個吸附層,見圖4(c);當kSF=2.0 時,/(kBT)?1,顆粒表面吸附達到飽和.這個吸附層會增大顆粒的尺寸,并增大了氣體分子與顆粒之間的碰撞面積,更有利于二者之間的動量傳遞,從而增大了納米顆粒所受的曳力.因此,當kSF=0.1 和kSF=0.5 時,MD 模擬值會大于(2)式的預測結(jié)果.

圖4 氣體分子在納米顆粒表面吸附狀況的模擬快照Fig.4.Snapshots of the adsorption of gas molecules (blue)on the nanoparticle (red) surface.

考慮到氣體分子在顆粒表面的吸附效應(yīng),本文提出需要引入有效顆粒半徑對(2)式進行修正.為了計算有效顆粒半徑,首先需要評估吸附層厚度.吸附層厚度與界面處氣體分子的分布情況緊密相關(guān).當kSF較小時,顆粒表面吸附的氣體分子較少,其等效的吸附層厚度較小;kSF逐漸增大時,顆粒表面吸附的氣體分子增多,吸附層厚度增大,顆粒的有效半徑增大.本文引入一個吸附參數(shù)α 來表征納米顆粒表面氣體分子的吸附強弱,α(kSF)=Na(kSF)/Nm,其中Na(kSF) 為各種kSF下吸附的氣體分子數(shù),Nm為kSF=2 時的吸附氣體分子數(shù).圖5 展示了吸附率α 隨kSF的變化情況.可以看到,吸附率α 隨著kSF而逐漸增大,并在kSF較大時逐漸趨于定值.當氣固結(jié)合強度較弱時,只有很少的氣體分子會被吸附在顆粒表面,因而隨著氣固結(jié)合強度的增強,吸附參數(shù)增大較快.對于較強的氣固結(jié)合強度,由于納米顆粒表面已經(jīng)吸附了較多的氣體分子并覆蓋在顆粒表面,入射氣體分子難以與納米顆粒發(fā)生相互作用,而更可能與已被吸附的氣體分子發(fā)生碰撞,此時入射氣體分子不容易被顆粒吸附.

圖5 吸附率α 隨氣固結(jié)合強度kSF 的變化Fig.5.Adsorption ratio α versus gas-particle coupling parameters kSF.

吸附率α 可表示各種kSF下的吸附層厚度與顆粒表面被完全覆蓋時吸附層厚度的比值,吸附層厚度與原始顆粒半徑之和即為有效顆粒半徑.有效顆粒半徑可以近似基于下式計算:

式中,δ 是界面吸附層的有效厚度,δm是當kSF=2 時的吸附層厚度.kSF=2 時的吸附層厚度可以在圖4 中得到.考慮到氣體分子在納米顆粒表面的吸附對Knudsen 數(shù)的計算也會有影響,因此本文對Kn也進行了修正,Kne=λ/Re.把修正后的有效半徑Re和Knudsen 數(shù)Kne代入(2)式中,將計算的曳力值重新繪制在圖3 中(見實心圓點).可以看到,修正后的理論公式結(jié)果與MD 模擬值非常吻合.

3.2 顆粒尺寸對納米顆粒曳力的影響

納米顆粒所受曳力與其顆粒尺寸有關(guān),本文通過改變固體原子數(shù)來改變納米顆粒的大小.固體原子數(shù)分別選取NS=141,225,369,555,767 和1061,根據(jù)(13)式可計算出對應(yīng)顆粒半徑R分別為0.724,0.869,1.014,1.158,1.303 和1.448 nm.同時,氣體分子數(shù)密度維持在1.97 nm-3,以確保流動區(qū)域處于過渡區(qū).圖6 展示了當kSF=0.01,0.1 和0.5 時納米顆粒的曳力隨著顆粒半徑的變化圖,其中氣體流速維持在大約1.5 m/s.結(jié)果表明,納米顆粒的曳力隨顆粒半徑的增大而增大.這是由于當顆粒半徑增大時,其表面積也會隨之越大,氣體分子更容易與顆粒發(fā)生碰撞,從而促進了氣體分子與納米顆粒之間的動量傳遞,因此顆粒會受到更大的曳力.此外,為了便于與理論公式進行比較,(2)式的計算值也被繪制在圖6 中.可以看到,當kSF=0.01 時,MD 模擬結(jié)果與(2)式的預測值相當吻合,然而當kSF=0.1 和kSF=0.5 時,MD 模擬結(jié)果則大于(2)式的計算結(jié)果.如上所述,當氣固結(jié)合強度較大時,氣體分子在顆粒表面的吸附作用不容忽視,需要對顆粒半徑進行修正.根據(jù)(14)式可得到修正的有效顆粒半徑,代入(2)式重新計算了曳力,并繪制于圖6 中(見藍色點線).由圖6 可知,此時考慮顆粒粒徑修正之后的計算結(jié)果與MD 模擬值更加吻合.

圖6 MD 模擬值與(2)式在3 種氣固結(jié)合強度下的對比結(jié)果Fig.6.Comparison of MD simulation results of the drag force with Eq.(2) at various kSF.

3.3 氣體流速對納米顆粒曳力的影響

氣體流速也是影響納米顆粒所受曳力的重要因素之一.在本文的模擬中,通過調(diào)整外力區(qū)域內(nèi)施加在每個氣體分子上的外力從而形成不同的氣體流速.圖7 展示了在3 種氣固結(jié)合強度下曳力隨著氣體流速的變化圖.結(jié)果表明,納米顆粒的曳力隨著氣體流速而增大.當氣體流速增大時,氣體分子與顆粒的碰撞更加劇烈,加強了二者之間的動量傳遞,因而顆粒具有更大的曳力.此外,(2)式也被繪制在圖7 中并與MD 模擬值進行對比.當kSF=0.01 時,氣體分子難以吸附在顆粒表面,MD 模擬值能與(2)式的結(jié)果良好吻合.而當kSF=0.1 和kSF=0.5 時,氣體分子與納米顆粒之間的相互作用較強,氣體分子容易被吸附在顆粒表面,并增大了納米顆粒的半徑,使得MD 模擬值大于(2)式的預測值,這與上文結(jié)論一致.考慮到氣體分子的吸附,需要采用修正的有效顆粒半徑.在過渡區(qū),由于氣體分子比較稀薄,可以忽略流速對氣體分子吸附效應(yīng)的影響.根據(jù)(14)式可計算得到有效顆粒半徑,代入(2)式并將計算結(jié)果繪制在圖7 中(見藍色點線).不難看出,修正后的(2)式與MD 模擬結(jié)果吻合更好.

圖7 MD 模擬值與(2)式在三種氣固結(jié)合強度下的對比結(jié)果Fig.7.Comparison of MD simulation results of the drag force with Eq.(2) at various kSF.

4 結(jié)論

本文基于分子動力學模擬方法,計算了過渡區(qū)內(nèi)納米顆粒所受的曳力,并與Li 和Wang[31]提出的納米顆粒曳力計算式進行對比.結(jié)果表明,對于納米顆粒而言,當氣固結(jié)合強度較弱時,分子動力學模擬結(jié)果可以與理論計算結(jié)果吻合較好;當氣固結(jié)合強度較強時,分子動力學模擬結(jié)果明顯大于理論計算的結(jié)果.進一步分析表明,此差異是由氣體分子在納米顆粒表面的吸附行為所導致的.本文基于氣體分子在顆粒表面的吸附特性,提出引入有效顆粒半徑修正,其過渡區(qū)內(nèi)曳力的理論計算結(jié)果與分子動力學模擬結(jié)果可以更好地吻合.因此,對于納米顆粒而言,氣體分子在顆粒表面的吸附效應(yīng)不可忽略.