納米通道內氣體剪切流動的分子動力學模擬?

張冉 謝文佳 常青 李樺

(國防科技大學空天科學學院,長沙 410073)

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

1 引 言

近年來,隨著微納機電系統(MEMS/NEMS)和微流動系統的蓬勃發展,納米尺度氣體流動現象受到廣泛關注[1?4].根據克努森數(Kn=λ/L,其中λ為氣體平均分子自由程,L為流動特征長度)的大小,錢學森[5]將流動區域劃分為4個區域:連續流區域(Kn<0.001)、滑移流區域(0.00110).對于納米尺度氣體流動,一方面,由于其流動特征長度與分子平均自由程處于同一量級,流動處于過渡流區域,稀薄氣體效應顯著,連續性假設早已不再適用;另一方面,由于流動的面積與體積比明顯增大,壁面作用力對流動的影響更加突出.

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

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

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

本文采用MD模擬方法,針對三維納米通道內的Couette氣體流動展開模擬,改變氣體-壁面相互作用勢能強度,獲得流場內速度、密度、正應力及剪切應力分布,并著重探討壁面力場對近壁區域氣體流動的影響規律.進一步根據剪切應力與克努森數的關系,獲得了不同氣體與表面相互勢能強度下的切向動量適應系數.

2 三維分子動力學模擬

2.1 物理模型

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

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

模擬中,采用截斷Lennard-Jones(L-J)6—12勢能函數[31]模擬氣體-氣體及氣體-壁面間的勢能相互作用,公式如下:

為了方便計算及分析,本文中壁面原子的質量及勢能參數與氣體分子一致,即σAr-Ar=σAr-Surface,εAr-Ar=εAr-Surface,氬原子直徑σAr=3.405× 10?10m,勢能參數εAr=1.67× 10?21J,rij為分子間的中心距離.C為氣體與壁面之間的勢能相互作用系數,下文簡稱為壁面勢能系數,通過調節其大小可以改變氣體分子與壁面原子之間的相互作用強度,考察不同種類氣體和壁面材料之間勢能作用強弱對流動的影響.公式中勢能截斷半徑rc為3.0σAr,超過該距離的勢能作用相對于氣體分子的動能可以忽略不計.模擬過程中采用元胞關聯法(link cell method)處理分子間的勢能作用力的計算.

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

2.2 模擬控制細節與統計方法

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

2.3 壓力張量的計算

本文采用了Irving-Kirwood(I-K)公式[35,36]計算系統中氣體的應力張量,對于一個包含N個氣體分子,體積為Vol的粒子系統,其表達式為

式中等式右邊第一項為應力張量的動力學項,wkl為分子間作用力產生的位力項.動力學項中,mi為粒子i質量,k和l分別代表笛卡爾坐標系中的不同方向,分別為粒子i在k和l方向上對應的速度分量,分別為粒子i所在位置所對應的k和l方向上氣體的平均速度.對于位力項,為粒子i與j之間的距離矢量在k方向上的分量,則是粒子j在l方向上對i粒子的分子間作用力分量.應力張量為三階反對稱張量,可表示為

其中sxx,syy,szz為正應力,其余各項為剪切應力,對處于靜止的平衡態理想氣體(無宏觀速度,忽略分子間作用力),剪切應力為零,根據理想氣體定律,氣體的壓力P和正應力的關系為

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

2.4 模擬區域設置

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

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

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

3 結果與討論

3.1 k=5.0時Couette流動的模擬結果與分析

Couette流動是以相同的速度沿相反方向分別移動上下壁面驅動氣體流動而獲得的.固定溫度下,通道中不同的密度會導致氣體分子平均自由程的變化,通過改變通道的高度使流動的克努森數保持一致.為了考察壁面力場對不同通道高度流動的影響設置第一組算例,模擬條件細節見表1.

圖2顯示了k=5.0時不同通道高度下的Couette流動系統快照圖.為了方便研究近壁面區域氣體的流動特性,將Z軸坐標原點移到壁面第一層原子的圓心位置.由于流場密度相對于通道中心為對稱分布,且在主流區域趨于定值,在圖中只截取了距下壁面2 nm范圍內流場的密度分布.圖3(a)為k=5.0時,三個高度分別為5.45,10.90,16.35 nm的納米通道中流場的密度分布.圖中Z=0 nm處的圓圈代表壁面第一層原子,Z=1.02 nm處的虛線為近壁區域與主流區的分界線.本文后續的其他圖中也采用了相同的表示方法.由圖3(a)可知,隨著氣體分子進入近壁區域后,在離壁面約0.8 nm處,氣體分子密度開始增大,在0.31 nm處達到峰值后急劇下降,直至變為零.壁面力場是導致氣體分子密度出現上述分布規律的主要原因.由于壁面力場的作用,氣體分子在壁面附近停留的時間變長,導致氣體分子在壁面附近出現集聚現象.隨著氣體分子距離壁面越來越近,氣體分子受到的壁面吸引力逐漸增大,氣體密度出現快速增大的趨勢.當氣體分子距離壁面很近時,壁面原子對氣體的吸引力轉變為排斥力,氣體密度隨之減小.圖3(b)為采用各通道中心的密度(3.73,1.86,1.25 kg/m3)分別對各個流場的密度進行歸一化處理后得到的歸一化密度分布,發現盡管通道內的密度不同,但是歸一化密度的分布卻是一致的,這說明近壁面區域的密度分布變化規律是壁面力場作用的結果,密度的變化趨勢只取決于氣體與壁面的相互作用勢能模型,而與通道高度和氣體密度無關.

表1 固定克努森數下Couette流動模擬條件Table 1.Simulation details of Couette flows under determined k.

圖2 k=5.0時不同通道高度下的Couette系統快照圖Fig.2.Snapshots of various height channel simulation domains at k=5.0.

圖3 (a)不同通道高度的流場密度分布;(b)分別采用通道中心密度對各個流場的密度進行歸一化處理得到的歸一化密度分布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情況下流場內的(a)正應力的動力學分量分布,(b)正應力的位力分量分布,(c)正應力分布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流動中正應力的分布規律,對通道高度為10.9 nm,密度為1.86 kg/m3條件下通道內的正應力分布進行詳細的分析.圖4(a)—圖4(c)分別為動力學項、位力項及二者相疊加后得到的正應力各分量分布曲線.由圖可知,在通道的主流區,由于氣體密度相對較小,導致氣體分子間間距較大,同時又不受壁面作用力的影響,正應力的位力項趨近于零,正應力主要由動力學項貢獻.在近壁面區域,由于壁面作用力的影響,位力項增大到和動力學項同一量級.此時,位力項和動力學項共同決定了近壁面區域的正應力分布.正應力分布在近壁區域呈現出各向異性的特征,其sxx,syy分量分布一致,但與szz分量的分布不同.表面力場的各向異性分布通過影響位力項的分布進而導致正應力各向異性.

圖5 (a)不同通道高度的流場內正應力sxx分量分布;(b)分別采用通道中心壓力對各個流場的正應力進行歸一化處理得到歸一化正應力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)下正應力的sxx及szz分量的分布曲線,分別采用各通道中心區域的壓力(222.87,109.2,76.23 kPa)對正應力各分量進行歸一化處理得到圖5(b)與圖6(b)所示的歸一化正應力分布.結果顯示近壁區域內不同通道高度下的正應力分布趨勢一致,這表明相同克努森數條件下,正應力分布規律是由表面力場決定的,與通道高度和密度無關.

圖6 (a)不同通道高度的流場內正應力的szz分量分布;(b)分別采用通道中心壓力對各個流場的正應力進行歸一化處理得到的歸一化正應力的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)分別對流動速度分布和高度進行歸一化處理,得到k=5.0時各流場的無量綱速度分布如圖7所示.由于流場速度分布為反對稱分布,在圖中只展示了通道下半部分的無量綱速度分布.此外,給出了TMAC為0.76和1.00時的線性玻爾茲曼方程(LBE)的理論解[11]與模擬結果相對比.結果顯示不同通道高度下的無量綱速度分布在主流區域呈線性分布且趨于重合,與TMAC為0.76時的理論解符合得比較好.壁面作用力的影響范圍為1.02 nm,由于通道的高度不同,導致近壁面區域所占通道的比例也不同.因此,通道越窄,壁面作用力對流動的影響也就越大.和理論預測的結果不同,近壁面區域內氣體的速度迅速增大并在表面達到最大值(約0.75UW),呈現出與解析解截然不同的分布規律.為了更加詳細地闡明壁面作用力的影響,截取了距下壁面2 nm范圍內流場的速度分布如圖8所示.從圖中可以看到,盡管通道的高度和剪切率不同,但氣體速度在近壁面區域內的分布卻是一致的.但隨著氣體遠離壁面,壁面不再對氣體分子施加作用力,剪切應力開始主導氣體的流動,由于流動的剪切率不同,速度分布逐漸出現分離.

圖7 k=5.0時,MD計算的流場無量綱速度分布及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時距下壁面2 nm范圍內流場的速度分布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)為不同通道高度及密度情況下流場內的剪切應力分布,由于流場剪切應力的分布相對于通道中心為對稱分布且在主流區域趨于定值,在圖中只截取了距下壁面2 nm范圍內流場的剪切應力分布.從圖9(a)中可以看到,大部分區域的剪切應力趨于定值,只是在近壁面約一個原子直徑(σAr=0.3405 nm)處出現突變.分別采用通道中心的切應力(?25.74,?12.74,?8.68 kPa)對各個流場的剪切應力進行歸一化處理,得到圖9(b)所示的歸一化剪切應力分布,結果顯示對于不同通道高度,歸一化剪切應力分布一致.這表明相同克努森數條件下,剪切應力的分布規律由表面力場決定,與通道高度和密度無關.

圖9 (a)不同通道高度的流場內剪切應力分布;(b)分別采用通道中心的剪切應力對各個流場的剪切應力進行歸一化處理得到的歸一化剪切應力分布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.

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

圖10 通道高度為10.90 nm,密度為1.86 kg/m3情況下流場內的剪切應力的動力學分量分布、剪切應力的位力分量分布及二者相加后的剪切應力分布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.

總結上述密度、速度及應力的分布規律,發現在相同克努森數條件下,近壁面區域內的密度、速度及應力分布規律與通道的高度及密度無關,僅與氣體與表面的相互勢能作用相關.由于氣體密度相對較小,當單個氣體分子與壁面相互作用時,周圍氣體分子對該分子的分子間作用力相對于壁面對其的作用力幾乎可以忽略,所以在近壁區域內,氣體的宏觀參數如密度、速度及應力僅表現出對壁面力場的依賴性關系.基于氣體動力學理論,處于自由分子流領域(Kn→∞)的Couette流動的剪切應力[37]為

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

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

3.2 勢能系數對流動的影響

為了研究氣體與壁面勢能作用強度對納米通道內氣體流動特性及切向動量適應系數的影響規律,針對通道高度為10.9 nm的Couette氣體流動展開模擬,改變壁面勢能系數的大小,同時保持通道中心區域的氣體密度一致,統計獲得流場內密度、速度及剪切應力分布.模擬條件設置見表2.

圖11為不同壁面勢能系數下Couette氣體流動系統快照圖.圖12為不同壁面勢能系數下的納米通道中距下壁面2 nm范圍內流場的密度分布.結果顯示不同壁面勢能系數下,主流區的密度分布一致.在近壁面區域內,隨著壁面勢能系數的增大,壁面原子對氣體分子的吸引力增強,氣體分子的密度呈指數增長,特別是在C=4.0時,在距壁面約一個原子直徑處,密度達到峰值453.29 kg/m3,表明壁面上出現了明顯的氣體分子吸附層,在流動系統快照圖中也可以明顯地觀察到吸附層.

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

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

表2 不同壁面勢能系數下Couette流動模擬條件Table 2.Simulation details of Couette flows under dif f erent potential coefficients.

圖11 C=1.0—4.0時的Couette流動系統快照Fig.11.Snapshots of Couette flows with various C(1.0–4.0)values.

圖12 不同壁面勢能系數下納米通道中的密度分布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)不同壁面勢能系數下流場歸一化速度分布;(b)不同壁面勢能系數下距下壁面2 nm范圍內流場的速度分布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.

根據前文中剪切應力與克努森數的關系,得到不同壁面勢能系數下氣體分子對壁面的切向動量適應系數如圖15所示.由圖可知,切向動量適應系數隨壁面勢能系數的增大而增大.當C=0.5時TMAC=0.63,按照切向動量適應系數的定義,有63%的氣體分子在壁面發生漫反射;當C=4.0時TMAC=0.96,此時幾乎所有氣體分子都會在壁面發生漫反射.上述結果說明勢能系數越大,壁面對氣體分子的勢能作用力越強,則氣體分子就越容易在壁面發生漫反射.

圖14 不同壁面勢能系數下流場內剪切應力分布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 不同壁面勢能系數下的切向動量適應系數Fig.15.TMAC values for k=5.0 flows under dif f erent gas-surface potential strength ratios.

4 結 論

采用分子動力學模擬方法研究了處于過渡流領域的納米通道內剪切力驅動的氣體流動,著重對近壁面區域內氣體流動特性進行了分析,并探討了壁面勢能系數對流動的影響規律,得到以下結論.

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

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

3)隨著壁面勢能系數的增大,近壁面區域內的密度隨之增大,當壁面系數增大到一定程度時,氣體分子吸附在表面形成了吸附層;壁面勢能系數的變化對流場內的速度分布影響是全局性的,近壁區域的速度隨著勢能系數的增大而增大,并最終和壁面速度一致,剪切應力隨著勢能系數的增大而增大,并導致主流區的剪切率隨之增大;切向動量適應系數隨著壁面勢能系數的增大而增大,氣體分子在勢能相互作用強的壁面更容易發生動量適應.

[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)[曹炳陽,陳民,過增元 2006物理學報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