駱雙雙 包福兵 林建忠

(中國計(jì)量學(xué)院計(jì)量測試工程學(xué)院 杭州 310018)

1 引言

20世紀(jì)90年代以來，隨著納米技術(shù)的進(jìn)步，各種微/納機(jī)電系統(tǒng)獲得了長足發(fā)展。隨著器件尺寸的減小，壁面作用變得重要，流體會出現(xiàn)明顯不同于常規(guī)尺度下的流動現(xiàn)象。例如在7 nm薄膜孔徑的多壁碳納米管中，Majumder等［1］發(fā)現(xiàn)壓力驅(qū)動流的流速比用宏觀方法統(tǒng)計(jì)出來的流速高4到5個數(shù)量級;在內(nèi)徑為1.6 nm的雙重壁碳納米管中，Holt等［2］發(fā)現(xiàn)水流量比連續(xù)流體力學(xué)模型統(tǒng)計(jì)出的水流量高3個數(shù)量級;在43 nm的碳納米管道中，Whitby等［3］發(fā)現(xiàn)水、乙醇和癸烷的輸運(yùn)能力比理論值高出了45倍。Hummer等［3］利用分子動力學(xué)方法模擬了納米通道中氣體和液體的流動情況，發(fā)現(xiàn)兩者的流量均比理論預(yù)測的大;Majumder等［1］認(rèn)為這個現(xiàn)象可能是由于納米通道中流體的粘度較小引起的。因此，對納米通道中流動特性的研究，特別是對流體粘度的研究具有重要的意義。

流體的粘度是流體分子微觀作用的宏觀表現(xiàn)，是表征流體對形變的抵抗隨形變速率變化的一種性質(zhì)。流體的流動是一個連續(xù)變形過程，內(nèi)部各部分速度不同，流層間存在速度差，快速流層力圖加快慢速流層，而慢速流層則力圖減慢快速流層，這種相互作用隨層間的速度差的增加而加劇，這種力圖減小速度差的作用力稱為粘性。

流體的粘度一般可通過各種粘度計(jì)測得。但是，在納米尺度通道中，目前的各種粘度測量裝置難以使用。因此，許多學(xué)者采用基于分子模型的方法來計(jì)算流體的粘度。陳煜等人［5］采用分子動力學(xué)方法計(jì)算了納米通道中液態(tài)氬和氦的平均粘度值。何慶東［6］采用將宏觀尺度下的流體運(yùn)動規(guī)律和分子動力學(xué)方法相結(jié)合的方法計(jì)算出粘度，模擬得出不同半徑、不同溫度下圓截面納米通道中的粘度。張暉等［7］利用分子動力學(xué)方法模擬了微孔介質(zhì)中氬在不同溫度、密度和孔徑下的粘度，并提出了微孔介質(zhì)中流體粘度的計(jì)算模型。夏冬昱等［8］利用Green-Kubo公式對二維和三維系統(tǒng)下液態(tài)氬的粘度進(jìn)行比較，并總結(jié)出二維模型代替三維模型的粘度計(jì)算方案。Liu等［9］模擬計(jì)算出圓形通道半徑為16 nm和32 nm時，粘度在壁面附近的值，并比較了不同固壁和流體分子作用強(qiáng)度下粘度的變化。Strak等［10］采用分子動力學(xué)方法模擬了氬在高壓強(qiáng)、低溫度下的粘度，結(jié)果和實(shí)驗(yàn)值的誤差不超過4%。

但是，到目前為止，已有的研究基本停留在通道平均粘度的模擬上，對通道內(nèi)粘度的分布和各種因素對粘度造成的影響分析較少，特別是關(guān)于氣體粘度的分布研究較少。因此，本文利用分子動力學(xué)方法，研究納米通道中壁面對氣體粘度分布的影響，以及不同通道高度和溫度對粘度的影響。

2 模型及方法

Couette流計(jì)算模型如圖1所示，下壁面固定，上壁面以速度vx往右運(yùn)動，從而帶動中間的氣體分子運(yùn)動。上下壁面均由鉑原子組成，鉑原子按照FCC(100)晶格排列，晶格常數(shù)為0.392 nm。氣體分子采用氬原子，模擬中x，z方向采用周期性邊界條件。

圖1 計(jì)算模型示意圖Fig.1 Schematic diagram of simulation box

計(jì)算中，氬分子之間，以及氬分子和鉑原子之間的作用勢能采用截?cái)嗥揭芁ennard-Jones(L-J)12-6勢能函數(shù)模型，其形式如下:

式中:rij是兩原子之間的距離，ε是勢阱深度，σ是分子直徑，rc為截?cái)喟霃?。在模擬中，采用以下廣泛采用的 L-J參數(shù):εAr=1.65 ×10-21J，σAr=0.34 ×10-9m，εPt-Ar=8.95 ×10-22J，σPt-Ar=0.31 ×10-9m，截?cái)喟霃綖?σAr。

在分子動力學(xué)方法中，采用鄰居列表法來計(jì)算原子間的受力，采用Velocity-Verlet方法來積分原子速度和位置，采用Nose-Hoover熱浴法來控制氣體的溫度，使系統(tǒng)保持為正則系綜(NVT)，時間步長為10×10-15s。要獲得通道中的速度梯度，要求速度剖面比較光滑，因此需要較大的采樣。在計(jì)算中，前500萬步使系統(tǒng)達(dá)到平衡，后10 000萬步統(tǒng)計(jì)相關(guān)的物理量。統(tǒng)計(jì)宏觀量時，沿y方向按dy為0.05 nm將通道分層，計(jì)算出每層的密度、速度、剪切應(yīng)力等值，就可得出這些量沿y方向的分布。

在牛頓流體中，剪切應(yīng)力和剪切速率之間存在線性的本構(gòu)關(guān)系，粘度與剪切速率無關(guān)。但是，在納米尺度通道的流動中，由于壁面原子與氣體分子的相互作用，靠近壁面的流體分子偏離熱力學(xué)平衡態(tài)。此時，牛頓定律中的線性關(guān)系可能不再成立。為了與常規(guī)粘度定義一致，可以如下定義等效粘度:

在分子動力學(xué)模擬中，可以獲得通道中沿y方向的Pxy和vx分布，因此根據(jù)式(2)，就可以獲得等效粘度沿著y方向的分布。

3 模型結(jié)果及討論

3.1 納米通道中粘度的分布

模擬中，上壁面速度設(shè)為221 m/s，下壁面靜止。通道高度為70 nm，溫度保持在300 K，氣體密度為7.22 kg/m3。根據(jù)計(jì)算可得氣體壓強(qiáng)為450 kPa，可得氣體的平均分子自由程為15.36 nm，根據(jù)通道高度，可得流動的Kn數(shù)為0.22，流動屬于滑移區(qū)。圖2給出了沿著y方向的速度剖面。從圖中可以看出，通道中間區(qū)域的氣體速度隨著高度線性變化，基本符合宏觀Couette流的特點(diǎn);而在上下壁面附近，速度呈現(xiàn)非線性變化。為了研究壁面附近的特性，圖2也給出了壁面附近氣體速度的分布。在非?？拷诿鎱^(qū)域(前3層)，由于壁面原子的排斥作用，這些區(qū)域不存在氣體原子，因此速度為0。而在第4層(距離壁面0.2 nm處)，存在一個速度躍變，存在速度滑移現(xiàn)象。在距離壁面0.2 nm到0.5 nm之間，氣體速度存在一定的波動，這與氣體分子在壁面附近的吸附有關(guān)。在0.5 nm之后，氣體速度逐漸增大。上壁面附近的原子速度也存在類似的現(xiàn)象。

圖3給出了沿著y方向剪切應(yīng)力的分布。在分子動力學(xué)方法中，剪切應(yīng)力由兩部分組成，一部分與分子的運(yùn)動速度相關(guān)，另一部分與分子之間的相互作用勢能相關(guān)。從圖中可以看出，剪切應(yīng)力在通道中呈對稱分布。在通道中間區(qū)域，剪切應(yīng)力保持在33.6 kPa。在這一區(qū)域，由于分子間距較大，氣體分子間的相互作用可以忽略，因此，分子的動能部分占據(jù)主導(dǎo)。而在壁面附近，剪切應(yīng)力變化劇烈，存在一個剪切應(yīng)力最大值，最大值為82.4 kPa，出現(xiàn)在0.4 nm 附近，這主要是由于壁面與氣體分子之間的相互作用引起的。在距離壁面0.5 nm以外，剪切應(yīng)力保持為33.6 kPa。

圖2 通道中的速度剖面Fig.2 Velocity profile across channel

圖3 通道中不同位置的剪切應(yīng)力Fig.3 Distribution of shear stress across channel

圖4給出了等效粘度沿著y方向的分布。從圖中可以看出，等效粘度在通道中呈對稱分布，中間部分等效粘度保持在一個定值，而在壁面附近，等效粘度逐漸變小，越靠近壁面，等效粘度值越小。根據(jù)等效粘度的分布特性，可以把流場分成兩個區(qū)域，中心區(qū)域和壁面區(qū)域。根據(jù)計(jì)算結(jié)果可以發(fā)現(xiàn)，壁面的影響范圍為20 nm左右。在中心區(qū)域，等效粘度的平均值為2.169 ×10-5Pa·s左右，和實(shí)驗(yàn)值2.276 ×10-5Pa·s［11］的相對誤差為 4.70% 。越靠近壁面，氣體的等效粘度越小。當(dāng)通道尺寸小于40 nm時，整個通道的等效粘度都將小于實(shí)驗(yàn)值。通道尺寸越小，通道的整體粘度越小。流體等效粘度的降低，導(dǎo)致了流體在通道中流動阻力的減小。因此，Majumder等［1］，Holt等［2］以及Whitby等人［3］在碳納米管的實(shí)驗(yàn)中觀察到了遠(yuǎn)大于宏觀方程預(yù)測到的流量。

圖4 通道中不同位置的等效粘度Fig.4 Distribution of effective viscosities across channel

3.2 通道高度的影響

為了對比不同通道高度下等效粘度的變化規(guī)律，比較了H=50、70、100和200 nm等4種通道高度的情況。計(jì)算中，氣體溫度為300 K，壓強(qiáng)為450 kPa。

對比不同通道高度下的等效粘度分布，發(fā)現(xiàn)存在著類似的現(xiàn)象，流場都可以分為兩個區(qū)域，中心區(qū)域和壁面區(qū)域，而且不同通道高度下壁面區(qū)域的厚度均為20nm左右。圖5給出了不同通道高度時中心區(qū)域的等效粘度和實(shí)驗(yàn)結(jié)果的對比，從圖中可以看出，該區(qū)域等效粘度隨著通道高度的增加增幅很小，最大值和最小值之間的相對誤差為6.43%，這可能是由于壁面厚度要略大于20 nm所致。最大值、最小值與實(shí)驗(yàn)值之間的相對誤差分別為0.55%，5.92%。在本文的計(jì)算中，氣體密度保持一致，通道高度越大，計(jì)算中用到的原子數(shù)越多，計(jì)算量也越大。為了獲得比較穩(wěn)定的結(jié)果，需要進(jìn)行大量的采樣(當(dāng)前模擬中均采用了10000萬步的結(jié)果取平均)，因此，沒有計(jì)算更大的通道高度。

圖6給出了不同通道高度時下壁面附近等效粘度的分布，從圖中可以看出，隨著與壁面距離的增大，等效粘度迅速增大，然后逐漸趨向于一個定值(實(shí)驗(yàn)粘度值)。不同通道高度時，等效粘度分布幾乎一致，說明等效粘度在壁面附近的分布與通道高度無關(guān)，主要是氣體分子與壁面原子的相互作用所致。

3.3 溫度的影響

圖5 中心區(qū)等效粘度隨通道高度的變化Fig.5 Effective viscosities of bulk region at different channel heights

圖6 不同通道高度時壁面附近等效粘度的分布Fig.6 Effective viscosities near the wall at different channel heights

圖7給出了通道高度是70 nm時，不同溫度(T=200、300、400和500 K)下通道中心區(qū)域等效粘度的變化，并與實(shí)驗(yàn)結(jié)果進(jìn)行了對比。從圖中可以看出，隨著溫度的增大，等效粘度也逐漸增大。圖8給出了不同溫度時壁面附近的等效粘度的分布。從圖中可以看出，在非?？拷诿娴膮^(qū)域，由于壁面原子和氣體分子的相互作用占主導(dǎo)地位，因此不同溫度下等效粘度的分布幾乎一致。而隨著氣體分子遠(yuǎn)離壁面，壁面對氣體分子的影響逐漸減弱，而氣體分子之間的間距較大，其相互作用可以忽略，此時，分子的動能部分占主導(dǎo)地位;溫度越高，氣體熱運(yùn)動越劇烈，氣體分子速度越快，越容易與其它分子碰撞，進(jìn)行動量的交換，從而使流層間的內(nèi)摩擦加劇，等效粘度越大。

4 總結(jié)與展望

圖7 等效粘度隨氣體溫度的變化Fig.7 Variation of bulk viscosities with gas temperature

圖8 不同溫度時壁面附近等效粘度的分布Fig.8 Effective viscosities near wall at different temperature

納米尺度通道中，流體會出現(xiàn)明顯不同于常規(guī)尺度的流動規(guī)律，這主要是壁面原子與流體分子的相互作用所致。本論文提出了等效粘度的概念，利用分子動力學(xué)方法，獲得了納米尺度Couette流中氣體的速度和剪切應(yīng)力分布，然后利用牛頓定律來求得通道中不同位置處氣體的等效粘度。模擬發(fā)現(xiàn)，在納米尺度通道中，氣體的等效粘度分布可以分為兩個區(qū)域，中心區(qū)域和壁面區(qū)域。在中心區(qū)域，由于距離壁面較遠(yuǎn)，此處的氣體粘度不受壁面的影響，計(jì)算得到的等效粘度與實(shí)驗(yàn)測得的粘度值符合得很好。而在壁面區(qū)域，氣體分子由于受到壁面分子的作用而偏離熱力學(xué)平衡狀態(tài)，等效粘度不是一個定值，越靠近壁面，等效粘度越小，壁面區(qū)域的厚度在20 nm左右。同時發(fā)現(xiàn)，通道高度對氣體的等效粘度幾乎沒有影響，而溫度越高，氣體的等效粘度越大。

利用等效粘度的概念，可以利用常規(guī)的Navier-Stokes方程結(jié)合獲得的等效粘度分布來模擬納米尺度通道中的流動特性，這樣可以擴(kuò)大Navier-Stokes方程的使用范圍，極大的提高納米尺度通道流動模擬的效率。

1 Majumder M，Chopra N，Andrews R，et al.Nanoscale hydrodynamics:enhanced flow in carbon nanotubes［J］.Nature，2005，438(7064):44-44.

2 Holt J K，Park H G，Wang Y，et al.Fast mass transport through sub-2-nanometer carbon nanotubes［J］.Science，2006，312(5776):1034-1037.

3 Whitby M，Cagnon L，Thanou M，et al.Enhanced fluid flow through nanoscale carbon pipes［J］.Nano letters，2008，8(9):2632-2637.

4 Hummer G，Rasaiah J C，Noworyta J P.Water conduction through the hydrophobic channel of a carbon nanotube［J］.Nature，2001，414(6860):188-190.

5 陳煜，陳碩，巨永林，等.超臨界L-J流體粘度的分子動力學(xué)模擬［J］.低溫工程，2008(4):32-37.Chen Yu，Chen Shuo，Ju Yonglin，et al.Molecular dynamics simulation of viscosities for supercritical L-J fluids［J］.Cryogenics，2008(4):32-37.

6 何慶東.納米尺度流體粘度計(jì)算的新方法研究［D］.秦皇島:燕山大學(xué)，2012.He Qingdong.Research on a new method to calculate flow viscosity in nanoflow problems［D］.Qinhuangdao:Yanshan University，2012.

7 張暉，張秉堅(jiān)，梁世強(qiáng)，等.微孔中簡單流體粘度的分子動力學(xué)模擬及關(guān)聯(lián)模型［J］.物理化學(xué)學(xué)報(bào)，2003，19(4):352-355.Zhang Hui，Zhang Bingjian，Liang Shiqiang，et al.Shear viscosity of simple fluids in porous media:molecular dynamics simulations and correlation models［J］.Acta Physico-Chimica Sinica，2003，19(4):352-355.

8 夏冬昱，王美芬，胡建軍.Green-Kubo公式在不同維數(shù)模擬系統(tǒng)中的比較分析［J］.武漢理工大學(xué)學(xué)報(bào)，2012，34(10):35-39.Xia Dongyu，Wang Meifen，Hu Jianjun.The comparison between 2-D and 3-D fluid models in molecular dynamics simulations［J］.Journal of Wuhan University of Technology，2012，34(10):35-39.

9 Liu Qixin.The study of Argon＇s viscosity near the wall of nanopore by molecular dynamics simulations［C］.Proceedings of the ASME 2013 4th International Conference on Micro Nanoscale Heat and Mass Transfer，Hong Kong:2013.

10 Pawel S，Stanislaw K.Molecular dynamic simulations of viscosity of argon at high pressures［J］.High Pressure Research，2008，28(4):469-476.

11 陳國邦，黃永華，包銳.低溫流體熱物理性質(zhì)［M］.北京:國防工業(yè)出版社，2006.