張冉 謝文佳 常青 李樺

(國(guó)防科技大學(xué)空天科學(xué)學(xué)院,長(zhǎng)沙 410073)

(2017年12月21日收到;2018年1月16日收到修改稿)

1 引 言

近年來(lái),隨著微納機(jī)電系統(tǒng)(MEMS/NEMS)和微流動(dòng)系統(tǒng)的蓬勃發(fā)展,納米尺度氣體流動(dòng)現(xiàn)象受到廣泛關(guān)注[1?4].根據(jù)克努森數(shù)(Kn=λ/L,其中λ為氣體平均分子自由程,L為流動(dòng)特征長(zhǎng)度)的大小,錢(qián)學(xué)森[5]將流動(dòng)區(qū)域劃分為4個(gè)區(qū)域:連續(xù)流區(qū)域(Kn<0.001)、滑移流區(qū)域(0.00110).對(duì)于納米尺度氣體流動(dòng),一方面,由于其流動(dòng)特征長(zhǎng)度與分子平均自由程處于同一量級(jí),流動(dòng)處于過(guò)渡流區(qū)域,稀薄氣體效應(yīng)顯著,連續(xù)性假設(shè)早已不再適用;另一方面,由于流動(dòng)的面積與體積比明顯增大,壁面作用力對(duì)流動(dòng)的影響更加突出.

對(duì)于滑移流區(qū)域的微尺度流動(dòng),Navier-Stokes方程仍然是有效的,但是需要采用確定流動(dòng)邊界滑移的滑移邊界條件.在經(jīng)典的Maxwell滑移模型基礎(chǔ)上,一些能反映流體與固體界面流動(dòng)特性的滑移模型[6,7]的出現(xiàn)在很大程度上解決了滑移區(qū)微尺度流動(dòng)模擬的邊界條件問(wèn)題.然而,對(duì)處于過(guò)渡流區(qū)域的納米尺度氣體流動(dòng),Navier-Stokes方程附加滑移邊界條件的方法已經(jīng)不再適用.從氣體動(dòng)力學(xué)理論出發(fā)的求解玻爾茲曼方程方法[8?10]及直接模擬Monte Carlo(DSMC)方法[11,12]被認(rèn)為是研究過(guò)渡流區(qū)域氣體流動(dòng)的有效數(shù)值模擬方法,但這兩種方法在處理流動(dòng)邊界時(shí)都采取了一定的簡(jiǎn)化,認(rèn)為氣體分子與壁面碰撞后會(huì)處于一種特定的狀態(tài),這類(lèi)處理方法無(wú)一例外都忽視了氣體分子與壁面的相互作用力對(duì)流動(dòng)的影響.

分子動(dòng)力學(xué)模擬(MD)方法[13]作為一種確定論的模擬方法,采用分子間勢(shì)能函數(shù)準(zhǔn)確地描述氣體與氣體之間及氣體與壁面之間的相互作用.隨著計(jì)算機(jī)運(yùn)算能力的提升,MD方法已經(jīng)成為微納尺度流動(dòng)領(lǐng)域最重要的研究手段之一.MD方法在處理邊界條件時(shí)不附加任何假設(shè)條件,僅采用氣體-壁面勢(shì)能相互作用模型描述氣體與壁面相互作用的特點(diǎn)顯示出了其獨(dú)有的優(yōu)勢(shì).學(xué)者們采用MD方法研究了微納尺度通道中的液體流動(dòng)問(wèn)題[14,15],發(fā)現(xiàn)了流體在壁面附近發(fā)生的密度分層、速度滑移及溫度跳躍等現(xiàn)象,極大地拓展了人們對(duì)微納尺度流動(dòng)的認(rèn)識(shí).

針對(duì)納米尺度氣體流動(dòng),國(guó)內(nèi)外學(xué)者采用MD方法做了許多相關(guān)的研究.切向動(dòng)量適應(yīng)系數(shù)是表征氣體分子和表面動(dòng)量與能量的交換程度的重要參數(shù),它與表面性質(zhì)的關(guān)系已經(jīng)被廣泛地研究[16?21].氣體分子的切向動(dòng)量適應(yīng)系數(shù)對(duì)表面的條件比較敏感,壁面材料種類(lèi)、溫度、表面粗糙度都對(duì)適應(yīng)系數(shù)有著明顯的影響.曹炳陽(yáng)等[22,23]采用非平衡分子動(dòng)力學(xué)方法對(duì)納米通道內(nèi)的Couette流動(dòng)進(jìn)行了研究,獲得了不同壁面條件下通道內(nèi)流體的速度及密度分布,探討了速度滑移與勢(shì)能作用強(qiáng)度間的關(guān)系.Xie和Liu[24]模擬了混合氣體在納米通道中的Poiseuille流動(dòng),發(fā)現(xiàn)由于壁面作用力的影響,通道內(nèi)不同成分的氣體分布不再均勻的特殊現(xiàn)象.Barisik和Beskok[25,26]采用小壁分子動(dòng)力學(xué)(SWMD)方法研究了通道內(nèi)氣體的流動(dòng)規(guī)律,發(fā)現(xiàn)近壁區(qū)域內(nèi)氣體的速度、密度及壓力變化趨勢(shì)僅由壁面力場(chǎng)決定,不受通道高度和密度的影響.但SWMD方法中假定壁面原子靜止,雖然這種簡(jiǎn)化方法節(jié)省了計(jì)算量,卻忽視了壁面原子熱運(yùn)動(dòng)對(duì)氣體分子的影響.Bao等[27,28]針對(duì)有限長(zhǎng)度的納米通道內(nèi)壓力驅(qū)動(dòng)流動(dòng)進(jìn)行了分子動(dòng)力學(xué)模擬研究,其模擬結(jié)果也發(fā)現(xiàn)了壁面力場(chǎng)主導(dǎo)近壁面區(qū)域流動(dòng)特性的現(xiàn)象.從前人的研究來(lái)看,壁面力場(chǎng)是主導(dǎo)納米通道中氣體流動(dòng)規(guī)律的一個(gè)主導(dǎo)因素,但其影響近壁區(qū)域的流動(dòng)特性的物理機(jī)理還需要進(jìn)一步研究.在壁面作用力的影響區(qū)域,同時(shí)也是氣體分子和壁面發(fā)生相互作用的區(qū)域,氣體與壁面的碰撞和氣體分子在表面的吸附對(duì)納米元器件的力學(xué)性能有著顯著的影響[29?31].因此,了解該區(qū)域內(nèi)氣體的流動(dòng)特性對(duì)納米元器件的設(shè)計(jì)和優(yōu)化具有非常重要的意義.

本文采用MD模擬方法,針對(duì)三維納米通道內(nèi)的Couette氣體流動(dòng)展開(kāi)模擬,改變氣體-壁面相互作用勢(shì)能強(qiáng)度,獲得流場(chǎng)內(nèi)速度、密度、正應(yīng)力及剪切應(yīng)力分布,并著重探討壁面力場(chǎng)對(duì)近壁區(qū)域氣體流動(dòng)的影響規(guī)律.進(jìn)一步根據(jù)剪切應(yīng)力與克努森數(shù)的關(guān)系,獲得了不同氣體與表面相互勢(shì)能強(qiáng)度下的切向動(dòng)量適應(yīng)系數(shù).

2 三維分子動(dòng)力學(xué)模擬

2.1 物理模型

Couette流動(dòng)的幾何結(jié)構(gòu)如圖1所示,單原子氣體氬均勻分布在距離為H的兩無(wú)限平板之間,在流動(dòng)系統(tǒng)的流向(X)和展向(Y)施加周期性邊界條件.

上下壁面以相同的速度UW沿流向做相對(duì)運(yùn)動(dòng),為了避免流動(dòng)出現(xiàn)壓縮性現(xiàn)象,UW選取0.2Ma.對(duì)于單原子氣體氬,其馬赫數(shù)Ma=其中氬原子氣體比熱比γ=5/3,玻爾茲曼常數(shù)kB=1.3806×10?23J·K?1,氣體溫度T=298 K,氬原子質(zhì)量mAr=6.63×10?26kg.

模擬中,采用截?cái)郘ennard-Jones(L-J)6—12勢(shì)能函數(shù)[31]模擬氣體-氣體及氣體-壁面間的勢(shì)能相互作用,公式如下:

為了方便計(jì)算及分析,本文中壁面原子的質(zhì)量及勢(shì)能參數(shù)與氣體分子一致,即σAr-Ar=σAr-Surface,εAr-Ar=εAr-Surface,氬原子直徑σAr=3.405× 10?10m,勢(shì)能參數(shù)εAr=1.67× 10?21J,rij為分子間的中心距離.C為氣體與壁面之間的勢(shì)能相互作用系數(shù),下文簡(jiǎn)稱(chēng)為壁面勢(shì)能系數(shù),通過(guò)調(diào)節(jié)其大小可以改變氣體分子與壁面原子之間的相互作用強(qiáng)度,考察不同種類(lèi)氣體和壁面材料之間勢(shì)能作用強(qiáng)弱對(duì)流動(dòng)的影響.公式中勢(shì)能截?cái)喟霃絩c為3.0σAr,超過(guò)該距離的勢(shì)能作用相對(duì)于氣體分子的動(dòng)能可以忽略不計(jì).模擬過(guò)程中采用元胞關(guān)聯(lián)法(link cell method)處理分子間的勢(shì)能作用力的計(jì)算.

對(duì)于壁面原子,其排布按FCC(100)晶格排布,晶格常數(shù)為1.6σAr.考慮壁面溫度對(duì)氣體分子運(yùn)動(dòng)的影響,壁面模型采用愛(ài)因斯坦固體模型[33],該模型認(rèn)為每個(gè)壁面原子都以相同的頻率在各自的平衡點(diǎn)附近做簡(jiǎn)諧振動(dòng),本文中壁面原子簡(jiǎn)諧振動(dòng)的彈性系數(shù)為36.75 N/m.

2.2 模擬控制細(xì)節(jié)與統(tǒng)計(jì)方法

將流場(chǎng)區(qū)域沿Z軸劃分為若干相同高度(0.34 nm)的子區(qū)域,采用統(tǒng)計(jì)平均的方法獲得各子區(qū)域的平均速度,將其標(biāo)記為對(duì)應(yīng)氣體分子的宏觀速度.在計(jì)算溫度的過(guò)程中,去除每個(gè)氣體分子的宏觀速度,從而獲得氣體的真實(shí)溫度.模擬系綜選取正則系綜(NVT),采用Nose-Hoover熱浴處理方法[34]將系統(tǒng)溫度控制在298 K.模擬開(kāi)始時(shí),氣體分子均勻分布在流場(chǎng)中,氣體分子的初始速度根據(jù)給定溫度下的平衡態(tài)分布抽樣獲得.速度積分采用Velocity-Verlet算法,時(shí)間步長(zhǎng)?t=0.002τ(≈4 fs),系統(tǒng)運(yùn)行2×106步后達(dá)到穩(wěn)定,穩(wěn)定后采用較長(zhǎng)的時(shí)間(3×107?t)對(duì)宏觀物理量(數(shù)密度,速度及壓力張量)進(jìn)行時(shí)間平均.將模擬盒子沿Z方向且平行于XY平面劃分為等間距(?H)的若干層,用于統(tǒng)計(jì)流場(chǎng)中各宏觀物理量的分布,為了獲取近壁面區(qū)域的流動(dòng)細(xì)節(jié),層間距取較小的高度(0.2σ=0.068 nm).

2.3 壓力張量的計(jì)算

本文采用了Irving-Kirwood(I-K)公式[35,36]計(jì)算系統(tǒng)中氣體的應(yīng)力張量,對(duì)于一個(gè)包含N個(gè)氣體分子,體積為Vol的粒子系統(tǒng),其表達(dá)式為

式中等式右邊第一項(xiàng)為應(yīng)力張量的動(dòng)力學(xué)項(xiàng),wkl為分子間作用力產(chǎn)生的位力項(xiàng).動(dòng)力學(xué)項(xiàng)中,mi為粒子i質(zhì)量,k和l分別代表笛卡爾坐標(biāo)系中的不同方向,分別為粒子i在k和l方向上對(duì)應(yīng)的速度分量,分別為粒子i所在位置所對(duì)應(yīng)的k和l方向上氣體的平均速度.對(duì)于位力項(xiàng),為粒子i與j之間的距離矢量在k方向上的分量,則是粒子j在l方向上對(duì)i粒子的分子間作用力分量.應(yīng)力張量為三階反對(duì)稱(chēng)張量,可表示為

其中sxx,syy,szz為正應(yīng)力,其余各項(xiàng)為剪切應(yīng)力,對(duì)處于靜止的平衡態(tài)理想氣體(無(wú)宏觀速度,忽略分子間作用力),剪切應(yīng)力為零,根據(jù)理想氣體定律,氣體的壓力P和正應(yīng)力的關(guān)系為

不難看出,I-K公式是對(duì)理想氣體狀態(tài)公式的修正,它考慮了氣體分子間相互作用力對(duì)氣體應(yīng)力張量的貢獻(xiàn).對(duì)于稠密氣體和液體而言,由于流體分子間間距較小,分子間的相互作用變得不可忽視,此外,壁面附近的流體分子不可避免地會(huì)受到表面力場(chǎng)影響,因此在計(jì)算壓力張量時(shí)必須要考慮分子間位力項(xiàng)的貢獻(xiàn).Couette流動(dòng)為剪切應(yīng)力驅(qū)動(dòng)的流動(dòng),剪切應(yīng)力τ=sxz=?szx?=0,其余剪切應(yīng)力均為零.

2.4 模擬區(qū)域設(shè)置

分子平均自由程(λ)是衡量氣體稀薄程度的重要參數(shù),假設(shè)氣體分子為硬球分子,其氣體動(dòng)力學(xué)定義[11]為

式中σHS為硬球分子直徑,n為分子數(shù)密度,Murat Barisik[28]采用MD方法計(jì)算了納米通道內(nèi)氣體氬在298 K時(shí)的硬球分子直徑為1.06σAr.

流動(dòng)在某個(gè)方向上采用周期性邊界條件時(shí),為了保證氣體分子間的碰撞充分,該方向上的模擬尺度至少應(yīng)為一個(gè)分子平均自由程.本文中的模擬盒子在X,Y方向的長(zhǎng)度均為λ,即模擬盒子的體積為λ·λ·H.Couette流動(dòng)的特征長(zhǎng)度為通道的高度H,流動(dòng)的克努森數(shù)Kn=λ/H,針對(duì)過(guò)渡流領(lǐng)域的氣體流動(dòng),線性玻爾茲曼方程的求解[11,12]中廣泛采用了修正克努森數(shù)的概念,為了方便分子動(dòng)力學(xué)模擬結(jié)果與之對(duì)比,本文中也采用了這種表示方法.

3 結(jié)果與討論

3.1 k=5.0時(shí)Couette流動(dòng)的模擬結(jié)果與分析

Couette流動(dòng)是以相同的速度沿相反方向分別移動(dòng)上下壁面驅(qū)動(dòng)氣體流動(dòng)而獲得的.固定溫度下,通道中不同的密度會(huì)導(dǎo)致氣體分子平均自由程的變化,通過(guò)改變通道的高度使流動(dòng)的克努森數(shù)保持一致.為了考察壁面力場(chǎng)對(duì)不同通道高度流動(dòng)的影響設(shè)置第一組算例,模擬條件細(xì)節(jié)見(jiàn)表1.

圖2顯示了k=5.0時(shí)不同通道高度下的Couette流動(dòng)系統(tǒng)快照?qǐng)D.為了方便研究近壁面區(qū)域氣體的流動(dòng)特性,將Z軸坐標(biāo)原點(diǎn)移到壁面第一層原子的圓心位置.由于流場(chǎng)密度相對(duì)于通道中心為對(duì)稱(chēng)分布,且在主流區(qū)域趨于定值,在圖中只截取了距下壁面2 nm范圍內(nèi)流場(chǎng)的密度分布.圖3(a)為k=5.0時(shí),三個(gè)高度分別為5.45,10.90,16.35 nm的納米通道中流場(chǎng)的密度分布.圖中Z=0 nm處的圓圈代表壁面第一層原子,Z=1.02 nm處的虛線為近壁區(qū)域與主流區(qū)的分界線.本文后續(xù)的其他圖中也采用了相同的表示方法.由圖3(a)可知,隨著氣體分子進(jìn)入近壁區(qū)域后,在離壁面約0.8 nm處,氣體分子密度開(kāi)始增大,在0.31 nm處達(dá)到峰值后急劇下降,直至變?yōu)榱?壁面力場(chǎng)是導(dǎo)致氣體分子密度出現(xiàn)上述分布規(guī)律的主要原因.由于壁面力場(chǎng)的作用,氣體分子在壁面附近停留的時(shí)間變長(zhǎng),導(dǎo)致氣體分子在壁面附近出現(xiàn)集聚現(xiàn)象.隨著氣體分子距離壁面越來(lái)越近,氣體分子受到的壁面吸引力逐漸增大,氣體密度出現(xiàn)快速增大的趨勢(shì).當(dāng)氣體分子距離壁面很近時(shí),壁面原子對(duì)氣體的吸引力轉(zhuǎn)變?yōu)榕懦饬?氣體密度隨之減小.圖3(b)為采用各通道中心的密度(3.73,1.86,1.25 kg/m3)分別對(duì)各個(gè)流場(chǎng)的密度進(jìn)行歸一化處理后得到的歸一化密度分布,發(fā)現(xiàn)盡管通道內(nèi)的密度不同,但是歸一化密度的分布卻是一致的,這說(shuō)明近壁面區(qū)域的密度分布變化規(guī)律是壁面力場(chǎng)作用的結(jié)果,密度的變化趨勢(shì)只取決于氣體與壁面的相互作用勢(shì)能模型,而與通道高度和氣體密度無(wú)關(guān).

表1 固定克努森數(shù)下Couette流動(dòng)模擬條件Table 1.Simulation details of Couette flows under determined k.

圖2 k=5.0時(shí)不同通道高度下的Couette系統(tǒng)快照?qǐng)DFig.2.Snapshots of various height channel simulation domains at k=5.0.

圖3 (a)不同通道高度的流場(chǎng)密度分布;(b)分別采用通道中心密度對(duì)各個(gè)流場(chǎng)的密度進(jìn)行歸一化處理得到的歸一化密度分布Fig.3.(a)Density and(b)normalized density distributions for k=5.0 flows inside 5.35,10.9,and 16.35 nm height channels,respectively.Normalizations are made using the channel center densities of 3.73,1.86,1.25 kg/m3.

圖4 通道高度為10.90 nm,密度為1.86 kg/m3情況下流場(chǎng)內(nèi)的(a)正應(yīng)力的動(dòng)力學(xué)分量分布,(b)正應(yīng)力的位力分量分布,(c)正應(yīng)力分布Fig.4.Kinetic(a)and virial(b)components of the normal stress distribution(c)for k=5.0 flows inside 10.9 nm height channels with the density of 1.86 kg/m3.

為了了解Couette流動(dòng)中正應(yīng)力的分布規(guī)律,對(duì)通道高度為10.9 nm,密度為1.86 kg/m3條件下通道內(nèi)的正應(yīng)力分布進(jìn)行詳細(xì)的分析.圖4(a)—圖4(c)分別為動(dòng)力學(xué)項(xiàng)、位力項(xiàng)及二者相疊加后得到的正應(yīng)力各分量分布曲線.由圖可知,在通道的主流區(qū),由于氣體密度相對(duì)較小,導(dǎo)致氣體分子間間距較大,同時(shí)又不受壁面作用力的影響,正應(yīng)力的位力項(xiàng)趨近于零,正應(yīng)力主要由動(dòng)力學(xué)項(xiàng)貢獻(xiàn).在近壁面區(qū)域,由于壁面作用力的影響,位力項(xiàng)增大到和動(dòng)力學(xué)項(xiàng)同一量級(jí).此時(shí),位力項(xiàng)和動(dòng)力學(xué)項(xiàng)共同決定了近壁面區(qū)域的正應(yīng)力分布.正應(yīng)力分布在近壁區(qū)域呈現(xiàn)出各向異性的特征,其sxx,syy分量分布一致,但與szz分量的分布不同.表面力場(chǎng)的各向異性分布通過(guò)影響位力項(xiàng)的分布進(jìn)而導(dǎo)致正應(yīng)力各向異性.

圖5 (a)不同通道高度的流場(chǎng)內(nèi)正應(yīng)力sxx分量分布;(b)分別采用通道中心壓力對(duì)各個(gè)流場(chǎng)的正應(yīng)力進(jìn)行歸一化處理得到歸一化正應(yīng)力sxx分量分布Fig.5.(a)sxxand(b)normalized sxxdistributions for k=5.0 flows inside 5.35,10.9,and 16.35 nm height channels,respectively.Normalizations are made using the channel center normal stress of 222.87,109.2,76.23 kPa.

圖5(a)與圖6(a)分別為不同通道高度(5.45,10.90,16.35 nm)下正應(yīng)力的sxx及szz分量的分布曲線,分別采用各通道中心區(qū)域的壓力(222.87,109.2,76.23 kPa)對(duì)正應(yīng)力各分量進(jìn)行歸一化處理得到圖5(b)與圖6(b)所示的歸一化正應(yīng)力分布.結(jié)果顯示近壁區(qū)域內(nèi)不同通道高度下的正應(yīng)力分布趨勢(shì)一致,這表明相同克努森數(shù)條件下,正應(yīng)力分布規(guī)律是由表面力場(chǎng)決定的,與通道高度和密度無(wú)關(guān).

圖6 (a)不同通道高度的流場(chǎng)內(nèi)正應(yīng)力的szz分量分布;(b)分別采用通道中心壓力對(duì)各個(gè)流場(chǎng)的正應(yīng)力進(jìn)行歸一化處理得到的歸一化正應(yīng)力的szz分量分布Fig.6.(a)szzand(b)normalized szzdistributions for k=5.0 flows inside 5.35,10.9 and 16.35 nm height channels,respectively.Normalizations are made using the channel center normal stress of 222.87,109.2,76.23 kPa.

采用壁面速度UW(63.48 m/s)及通道高度H(5.45,10.90,16.35 nm)分別對(duì)流動(dòng)速度分布和高度進(jìn)行歸一化處理,得到k=5.0時(shí)各流場(chǎng)的無(wú)量綱速度分布如圖7所示.由于流場(chǎng)速度分布為反對(duì)稱(chēng)分布,在圖中只展示了通道下半部分的無(wú)量綱速度分布.此外,給出了TMAC為0.76和1.00時(shí)的線性玻爾茲曼方程(LBE)的理論解[11]與模擬結(jié)果相對(duì)比.結(jié)果顯示不同通道高度下的無(wú)量綱速度分布在主流區(qū)域呈線性分布且趨于重合,與TMAC為0.76時(shí)的理論解符合得比較好.壁面作用力的影響范圍為1.02 nm,由于通道的高度不同,導(dǎo)致近壁面區(qū)域所占通道的比例也不同.因此,通道越窄,壁面作用力對(duì)流動(dòng)的影響也就越大.和理論預(yù)測(cè)的結(jié)果不同,近壁面區(qū)域內(nèi)氣體的速度迅速增大并在表面達(dá)到最大值(約0.75UW),呈現(xiàn)出與解析解截然不同的分布規(guī)律.為了更加詳細(xì)地闡明壁面作用力的影響,截取了距下壁面2 nm范圍內(nèi)流場(chǎng)的速度分布如圖8所示.從圖中可以看到,盡管通道的高度和剪切率不同,但氣體速度在近壁面區(qū)域內(nèi)的分布卻是一致的.但隨著氣體遠(yuǎn)離壁面,壁面不再對(duì)氣體分子施加作用力,剪切應(yīng)力開(kāi)始主導(dǎo)氣體的流動(dòng),由于流動(dòng)的剪切率不同,速度分布逐漸出現(xiàn)分離.

圖7 k=5.0時(shí),MD計(jì)算的流場(chǎng)無(wú)量綱速度分布及TMAC=0.76與1.00下的線化玻爾茲曼方程解析解Fig.7.k=5.0 flow normalized velocity profiles of MD results and linearized Boltzmann solutions at TAMC=0.76 and 1.00.

圖8 k=5.0時(shí)距下壁面2 nm范圍內(nèi)流場(chǎng)的速度分布Fig.8.Velocity profiles of k=5.0 flows as a function of channel height within 2 nm distance from the bottom walls.

圖9(a)為不同通道高度及密度情況下流場(chǎng)內(nèi)的剪切應(yīng)力分布,由于流場(chǎng)剪切應(yīng)力的分布相對(duì)于通道中心為對(duì)稱(chēng)分布且在主流區(qū)域趨于定值,在圖中只截取了距下壁面2 nm范圍內(nèi)流場(chǎng)的剪切應(yīng)力分布.從圖9(a)中可以看到,大部分區(qū)域的剪切應(yīng)力趨于定值,只是在近壁面約一個(gè)原子直徑(σAr=0.3405 nm)處出現(xiàn)突變.分別采用通道中心的切應(yīng)力(?25.74,?12.74,?8.68 kPa)對(duì)各個(gè)流場(chǎng)的剪切應(yīng)力進(jìn)行歸一化處理,得到圖9(b)所示的歸一化剪切應(yīng)力分布,結(jié)果顯示對(duì)于不同通道高度,歸一化剪切應(yīng)力分布一致.這表明相同克努森數(shù)條件下,剪切應(yīng)力的分布規(guī)律由表面力場(chǎng)決定,與通道高度和密度無(wú)關(guān).

圖9 (a)不同通道高度的流場(chǎng)內(nèi)剪切應(yīng)力分布;(b)分別采用通道中心的剪切應(yīng)力對(duì)各個(gè)流場(chǎng)的剪切應(yīng)力進(jìn)行歸一化處理得到的歸一化剪切應(yīng)力分布Fig.9. (a)Dimensional and(b)normalized shear stress distributions of for k=5.0 flows inside 5.35,10.9,and 16.35 nm height channels,respectively.

為了解釋表面力場(chǎng)對(duì)剪切應(yīng)力的影響規(guī)律,圖10給出了通道高度為10.9 nm、密度為1.86 kg/m3情況下流場(chǎng)內(nèi)剪切應(yīng)力的動(dòng)力學(xué)項(xiàng)、位力項(xiàng)及二者相疊加后得到的剪切應(yīng)力分布曲線.從圖中可以看到,剪切應(yīng)力的位力項(xiàng)在大部分區(qū)域?yàn)榱?剪切應(yīng)力的大小由動(dòng)力學(xué)項(xiàng)決定.但是在近壁面約一個(gè)原子直徑(σAr=0.3405 nm)的范圍內(nèi),剪切應(yīng)力的位力項(xiàng)突然增大,同時(shí)動(dòng)力學(xué)項(xiàng)減小,值得注意的是,在這個(gè)區(qū)域內(nèi)壁面排斥力是占主導(dǎo)地位的.

圖10 通道高度為10.90 nm,密度為1.86 kg/m3情況下流場(chǎng)內(nèi)的剪切應(yīng)力的動(dòng)力學(xué)分量分布、剪切應(yīng)力的位力分量分布及二者相加后的剪切應(yīng)力分布Fig.10.Kinetic and virial components of the shear stress distribution for k=5.0 flows inside 10.9 nm height channels with the density of 1.86 kg/m3.

總結(jié)上述密度、速度及應(yīng)力的分布規(guī)律,發(fā)現(xiàn)在相同克努森數(shù)條件下,近壁面區(qū)域內(nèi)的密度、速度及應(yīng)力分布規(guī)律與通道的高度及密度無(wú)關(guān),僅與氣體與表面的相互勢(shì)能作用相關(guān).由于氣體密度相對(duì)較小,當(dāng)單個(gè)氣體分子與壁面相互作用時(shí),周?chē)鷼怏w分子對(duì)該分子的分子間作用力相對(duì)于壁面對(duì)其的作用力幾乎可以忽略,所以在近壁區(qū)域內(nèi),氣體的宏觀參數(shù)如密度、速度及應(yīng)力僅表現(xiàn)出對(duì)壁面力場(chǎng)的依賴(lài)性關(guān)系.基于氣體動(dòng)力學(xué)理論,處于自由分子流領(lǐng)域(Kn→∞)的Couette流動(dòng)的剪切應(yīng)力[37]為

式中ρ為來(lái)流密度,R為理想氣體常數(shù),TW為壁面溫度(298 K);α為切向動(dòng)量適應(yīng)系數(shù)(TMAC),它代表氣體分子在壁面發(fā)生漫反射的部分,而剩下的(1?α)部分則發(fā)生鏡面反射.根據(jù)之前近壁區(qū)域內(nèi)的流動(dòng)規(guī)律,我們推斷切向動(dòng)量適應(yīng)系數(shù)也應(yīng)與通道的高度與密度無(wú)關(guān),其大小僅取決于氣體與壁面的相互勢(shì)能用.對(duì)于過(guò)渡流領(lǐng)域的流動(dòng),由于需要考慮有限克努森數(shù)對(duì)剪切應(yīng)力的影響,不能直接用(6)式計(jì)算切向動(dòng)量適應(yīng)系數(shù)的大小,因此采用(7)式對(duì)處于過(guò)渡流領(lǐng)域的Couette流動(dòng)的剪切應(yīng)力進(jìn)行修正[38],該修正被證明對(duì)過(guò)渡流和滑移流領(lǐng)域的流動(dòng)都是有效的.

將(6)式與(7)式聯(lián)立,結(jié)合分子動(dòng)力學(xué)模擬計(jì)算的剪切應(yīng)力值,得到氣體分子對(duì)壁面的TMAC值.對(duì)于本文中的三個(gè)算例,通道中心的剪切應(yīng)力分布分別為?47.43,?23.71,?15.89 kPa.通過(guò)計(jì)算,對(duì)應(yīng)的TMAC值分別為0.76,0.75,0.76,結(jié)果表明TMAC的大小趨于定值,且與前文中線性玻爾茲曼方程所預(yù)測(cè)的TMAC值0.76相符合,直接證明了切向動(dòng)量適應(yīng)系數(shù)大小與通道高度及密度無(wú)關(guān)的推論的正確性.

3.2 勢(shì)能系數(shù)對(duì)流動(dòng)的影響

為了研究氣體與壁面勢(shì)能作用強(qiáng)度對(duì)納米通道內(nèi)氣體流動(dòng)特性及切向動(dòng)量適應(yīng)系數(shù)的影響規(guī)律,針對(duì)通道高度為10.9 nm的Couette氣體流動(dòng)展開(kāi)模擬,改變壁面勢(shì)能系數(shù)的大小,同時(shí)保持通道中心區(qū)域的氣體密度一致,統(tǒng)計(jì)獲得流場(chǎng)內(nèi)密度、速度及剪切應(yīng)力分布.模擬條件設(shè)置見(jiàn)表2.

圖11為不同壁面勢(shì)能系數(shù)下Couette氣體流動(dòng)系統(tǒng)快照?qǐng)D.圖12為不同壁面勢(shì)能系數(shù)下的納米通道中距下壁面2 nm范圍內(nèi)流場(chǎng)的密度分布.結(jié)果顯示不同壁面勢(shì)能系數(shù)下,主流區(qū)的密度分布一致.在近壁面區(qū)域內(nèi),隨著壁面勢(shì)能系數(shù)的增大,壁面原子對(duì)氣體分子的吸引力增強(qiáng),氣體分子的密度呈指數(shù)增長(zhǎng),特別是在C=4.0時(shí),在距壁面約一個(gè)原子直徑處,密度達(dá)到峰值453.29 kg/m3,表明壁面上出現(xiàn)了明顯的氣體分子吸附層,在流動(dòng)系統(tǒng)快照?qǐng)D中也可以明顯地觀察到吸附層.

通道內(nèi)的無(wú)量綱速度分布如圖13(a)所示,在整個(gè)流動(dòng)區(qū)域氣體速度隨著壁面勢(shì)能系數(shù)的增大而增大;在流動(dòng)的主流區(qū)域,速度為線性分布,線性分布的斜率隨著壁面勢(shì)能系數(shù)的增大而增大.近壁面區(qū)域的速度分布如圖13(b)所示,由于壁面力場(chǎng)的作用,導(dǎo)致氣體速度急劇增大,直到在距壁面約一個(gè)分子直徑處達(dá)到峰值.氣體速度的峰值隨著勢(shì)能系數(shù)的增大而增大,當(dāng)C=4.0時(shí),由于氣體分子被壁面大量吸附而導(dǎo)致速度峰值和壁面速度達(dá)到了一致,壁面速度滑移消失.

速度線性分布的斜率即Couette流動(dòng)的剪切率,剪切率和剪切應(yīng)力的大小成正比關(guān)系.不同壁面勢(shì)能系數(shù)下的剪切應(yīng)力分布如圖14所示,剪切應(yīng)力在主流區(qū)域?yàn)槎ㄖ?其大小隨著勢(shì)能系數(shù)的增大而增大,而且越來(lái)越逼近切向動(dòng)量適應(yīng)系數(shù)為1.0時(shí)的理論值(21.34 kPa).在靠近壁面約一個(gè)原子直徑(0.3405 nm)的區(qū)域內(nèi),剪切應(yīng)力突然增大到峰值,該峰值也隨著勢(shì)能系數(shù)的增大而增大.值的注意的是,當(dāng)C=4.0時(shí),由于壁面的氣體分子吸附層的存在,導(dǎo)致了剪切應(yīng)力在近壁面區(qū)域的分布出現(xiàn)了一定的扭曲現(xiàn)象.

表2 不同壁面勢(shì)能系數(shù)下Couette流動(dòng)模擬條件Table 2.Simulation details of Couette flows under dif f erent potential coefficients.

圖11 C=1.0—4.0時(shí)的Couette流動(dòng)系統(tǒng)快照Fig.11.Snapshots of Couette flows with various C(1.0–4.0)values.

圖12 不同壁面勢(shì)能系數(shù)下納米通道中的密度分布Fig.12.Density distribution for k=5.0 flows with various C values inside 10.9 nm height channels with the density of 1.86 kg/m3in bulk flow.

圖13 (a)不同壁面勢(shì)能系數(shù)下流場(chǎng)歸一化速度分布;(b)不同壁面勢(shì)能系數(shù)下距下壁面2 nm范圍內(nèi)流場(chǎng)的速度分布Fig.13.(a)Normalized velocity profiles of flows with various gas-surface potential strength ratios as a function of normalized channel height(Z/H);(b)velocity profiles of flows with various gas-surface potential strength ratios as a function of channel height within 2 nm distance from the bottom walls.

根據(jù)前文中剪切應(yīng)力與克努森數(shù)的關(guān)系,得到不同壁面勢(shì)能系數(shù)下氣體分子對(duì)壁面的切向動(dòng)量適應(yīng)系數(shù)如圖15所示.由圖可知,切向動(dòng)量適應(yīng)系數(shù)隨壁面勢(shì)能系數(shù)的增大而增大.當(dāng)C=0.5時(shí)TMAC=0.63,按照切向動(dòng)量適應(yīng)系數(shù)的定義,有63%的氣體分子在壁面發(fā)生漫反射;當(dāng)C=4.0時(shí)TMAC=0.96,此時(shí)幾乎所有氣體分子都會(huì)在壁面發(fā)生漫反射.上述結(jié)果說(shuō)明勢(shì)能系數(shù)越大,壁面對(duì)氣體分子的勢(shì)能作用力越強(qiáng),則氣體分子就越容易在壁面發(fā)生漫反射.

圖14 不同壁面勢(shì)能系數(shù)下流場(chǎng)內(nèi)剪切應(yīng)力分布Fig.14.Shear stress distributions of flows with various gas-surface potential strength ratios as a function of channel height(Z)within 2 nm distance from the bottom walls.

圖15 不同壁面勢(shì)能系數(shù)下的切向動(dòng)量適應(yīng)系數(shù)Fig.15.TMAC values for k=5.0 flows under dif f erent gas-surface potential strength ratios.

4 結(jié) 論

采用分子動(dòng)力學(xué)模擬方法研究了處于過(guò)渡流領(lǐng)域的納米通道內(nèi)剪切力驅(qū)動(dòng)的氣體流動(dòng),著重對(duì)近壁面區(qū)域內(nèi)氣體流動(dòng)特性進(jìn)行了分析,并探討了壁面勢(shì)能系數(shù)對(duì)流動(dòng)的影響規(guī)律,得到以下結(jié)論.

1)由于壁面作用力的影響,在距壁面約1 nm的范圍內(nèi)形成了不同于主流區(qū)的近壁區(qū)域.在該區(qū)域內(nèi),氣體的密度、正應(yīng)力由于壁面力場(chǎng)的影響出現(xiàn)劇烈的變化.壁面力場(chǎng)的影響主導(dǎo)了正應(yīng)力在近壁區(qū)域的各向異性分布特性.在相同的克努森數(shù)條件下,對(duì)于不同通道高度及密度的情況,將密度及正應(yīng)力歸一化處理后呈現(xiàn)出分布的一致性;氣體在壁面上的速度滑移明顯小于理論預(yù)測(cè)值,且近壁區(qū)域的速度分布不受通道高度及氣體密度的影響.

2)在不受壁面作用力影響的主流區(qū)域,氣體的流動(dòng)規(guī)律符合氣體動(dòng)力學(xué)定律,密度、正應(yīng)力,剪切應(yīng)力均為恒定值.氣體速度分布符合應(yīng)力-應(yīng)變的線性響應(yīng)關(guān)系,通過(guò)剪切應(yīng)力與克努森數(shù)的關(guān)聯(lián),預(yù)測(cè)了氣體分子對(duì)特定壁面條件下的切向動(dòng)量適應(yīng)系數(shù).將得到的切向動(dòng)量適應(yīng)系數(shù)應(yīng)用到線化玻爾茲曼方程的求解中,得到的速度分布理論解在主流區(qū)域與分子動(dòng)力學(xué)模擬結(jié)果符合得很好,證明通過(guò)剪切應(yīng)力計(jì)算切向動(dòng)量適應(yīng)系數(shù)是一種有效的預(yù)測(cè)方法.

3)隨著壁面勢(shì)能系數(shù)的增大,近壁面區(qū)域內(nèi)的密度隨之增大,當(dāng)壁面系數(shù)增大到一定程度時(shí),氣體分子吸附在表面形成了吸附層;壁面勢(shì)能系數(shù)的變化對(duì)流場(chǎng)內(nèi)的速度分布影響是全局性的,近壁區(qū)域的速度隨著勢(shì)能系數(shù)的增大而增大,并最終和壁面速度一致,剪切應(yīng)力隨著勢(shì)能系數(shù)的增大而增大,并導(dǎo)致主流區(qū)的剪切率隨之增大;切向動(dòng)量適應(yīng)系數(shù)隨著壁面勢(shì)能系數(shù)的增大而增大,氣體分子在勢(shì)能相互作用強(qiáng)的壁面更容易發(fā)生動(dòng)量適應(yīng).

[1]Ekinci K L,Roukes M L 2005Rev.Sci.Instrum.76 061101

[2]Verbridge S S,Craighead H G,Parpia J M 2008Appl.Phys.Lett.92 013112

[3]Boettcher U,Li H,Callafon R A,Talke F E 2011IEEE T.Magn.47 1823

[4]Song H Q,Yu M X,Z H U W Y,Zhang Y,Jiang S X 2013Chin.Phys.Lett.30 014701

[5]Tsien H S 1964J.Aero.Sci.13 653

[6]Zhang Z Q,Zhang H W,Ye H F 2009Appl.Phys.Lett.95 154101

[7]Zhang H W,Zhang Z Q,Zheng Y G,Ye H F 2010Phys.Rev.E81 066303

[8]Sone Y,Takata S,Ohwada T 1990Euro.J.Mech.B:Fluids9 273

[9]Taheri P,Torrilhon M,Struchtrup H 2009Phys.Fluids21 017102

[10]Dehdashti E 2016Chin.Phys.B25 024702

[11]Bird G A 1994Molecular Gas Dynamics and the Direct Simulation of Gas Flows(Oxford:Oxford University Press)pp199–206

[12]Fan J,Shen C 1999J.Comput.Phys.167 393

[13]Rapaport D C 2004The Art of Molecular Dynamics Simulation(New York:Cambridge University Press)pp4–5

[14]Thompson P A,Troian S M 1997Nature389 360

[15]Zhu Y X,Granick S 2002Phys.Rev.Lett.88 106102

[16]Cao B Y,Sun J,Chen M,Guo Z Y 2009Int.J.Mol.Sci.10 4638

[17]Finger G W,Kapat J S,Bhattacharya A 2007J.Fluid Eng.-T.ASME129 31

[18]Sun J,Li Z X 2008Mol.Phys.106 2325

[19]Sun J,Li Z X 2010Comput.Fluids39 1645

[20]Markvoort A J,Hilbers P A J,Nedea S V 2005Phys.Rev.E71 066702

[21]Arya G,Chang H C,Maginn E J 2003Mol.Simul.29 697

[22]Cao B Y,Chen M,Guo Z Y 2005Appl.Phys.Lett.86 091905

[23]Cao B Y,Chen M,Guo Z Y 2006Acta Phys.Sin.55 5305(in Chinese)[曹炳陽(yáng),陳民,過(guò)增元 2006物理學(xué)報(bào)55 5305]

[24]Xie H,Liu C 2012AIP Adv.2 042126

[25]Barisik M,Beskok A 2011Microfluid Nanofluid11 269

[26]Barisik M,Beskok A 2012Microfluid Nanofluid13 789

[27]Bao F B,Huang Y L,Qiu L M,Lin J Z 2014Mol.Phys.113 561

[28]Bao F B,Huang Y L,Zhang Y H,Lin J Z 2015Microfluid Nanofluid18 1075

[29]Yang Y T,Callegari C,Feng X L,Roukes M L 2011Nano Lett.11 1753

[30]Zhang W M,Meng G,Zhou J B,Chen J Y 2009Sensors9 3854

[31]Wu L,Bogy D B 2002J.Tribol.-T.ASME124 562

[32]Allen M P,Tildesley D J 1991Computer Simulation of Liquids(Oxford:Oxford University Press)pp145–146

[33]Hook J R,Hall H E 1991Solid State Physics(Chichester:Wiley)pp96–106

[34]Evans D J,Hoover W G 1986Annu.Rev.Fluid Mech.18 243

[35]Irving J H,Kirkwood J G 1950J.Chem.Phys.18 817

[36]Todd B D,Evans D J,Davis P J 1995Phys.Rev.E52 1627

[37]Fukui S,Shimada H,Yamane K,Matsuoka H 2005Microsyst.Technol.11 805

[38]Bahukudumbi P,Park J H,Beskok A 2003Microscale Thermophy.Eng.7 291