均勻剪切流中粒子無(wú)規(guī)行走的分子動(dòng)力學(xué)模擬研究

， ,

(北京航空航天大學(xué)a.航空科學(xué)與工程學(xué)院;b.航空科學(xué)與工程學(xué)院；c.國(guó)家計(jì)算流體力學(xué)重點(diǎn)實(shí)驗(yàn)室， 北京 100191)

0 引 言

對(duì)于模擬真實(shí)的稀薄氣體分子與真實(shí)物面的相互作用[1]，要選擇一個(gè)計(jì)算區(qū)域。為保證流動(dòng)充分發(fā)展，這個(gè)計(jì)算區(qū)域至少要大于分子平均自由程[2]。分子平均自由程與分子擴(kuò)散系數(shù)成正比，而擴(kuò)散系數(shù)又與分子無(wú)規(guī)行走距離相關(guān)，但是在均勻剪切率情況下，平均自由程與剪切率是什么關(guān)系還不是很清處，因此要研究平均自由程和剪切率的關(guān)系，通過(guò)分子動(dòng)力學(xué)模擬的方法[3]來(lái)研究不同均勻剪切率Couette流與分子無(wú)規(guī)行走曲線的斜率的關(guān)系，這樣可以為該計(jì)算區(qū)域大小的選擇提供一個(gè)參考，節(jié)省計(jì)算量。

1 模擬細(xì)節(jié)

1.1 基本參數(shù)設(shè)定

選用一般力場(chǎng)中最常見(jiàn)的非鍵結(jié)形式為L(zhǎng)ennard-Jones(LJ)勢(shì)能[4]。此勢(shì)能又稱為12-6勢(shì)，其數(shù)學(xué)表達(dá)式為：

(1)

σ與ε為勢(shì)參數(shù)，令σ為無(wú)量綱特征長(zhǎng)度，并對(duì)其他長(zhǎng)度變量進(jìn)行無(wú)量綱化，r為兩粒子之間的距離。

模擬的物理模型如下如圖(1)所示：

圖1 分子動(dòng)力模擬Couette流的模型

A盒子中粒子是所要計(jì)算的一系列粒子，而B(niǎo)盒子與C盒子中的粒子是A盒子中的鏡像粒子，但是這并不會(huì)維持一個(gè)穩(wěn)定的Couette流動(dòng)，為了獲得沿著z軸方向的速度梯度，選擇讓B盒子沿著x軸正方向開(kāi)始以速度Vd運(yùn)動(dòng)，而C盒子沿著z軸負(fù)方向具有一個(gè)Vd的運(yùn)動(dòng)。在B盒子與C盒子的拖動(dòng)之下，形成了具有一定剪切率的Couette流動(dòng)。剪切率γ=Vd/Lz,其中Lz為A盒子的z軸高度。選擇Lz為40。

分子的初始位置按照面心立方(FCC)晶體結(jié)構(gòu)布置，粒子總數(shù)為2048個(gè)恒定。整個(gè)模擬的A盒子的尺寸大小為40*40*40。每個(gè)粒子初始速度給定。整個(gè)系統(tǒng)的溫度是不變的。采用蛙跳算法來(lái)計(jì)算粒子的加速度，速度和位置時(shí)。選用Andersen的定溫計(jì)算法[5]，運(yùn)用常用的分子熱運(yùn)動(dòng)公式Ct=3×k×T/2，其中εt為分子熱運(yùn)動(dòng)平動(dòng)動(dòng)能，k為波爾茨曼常數(shù)，調(diào)節(jié)每一步的系統(tǒng)溫度，使其保持定值。初始系統(tǒng)溫度T的無(wú)量綱值為0.75。

為了保持粒子數(shù)密度守恒，應(yīng)用改進(jìn)的周期性邊界條件[6]。而與該文章不同的是，沒(méi)有忽略粒子x方向的受力。

2 模擬結(jié)果與討論

2.1 無(wú)規(guī)行走曲線及斜率

y與z方向無(wú)規(guī)行走，可以通過(guò)這兩個(gè)方向的熱運(yùn)動(dòng)速度定義。當(dāng)粒子初始位置布置好之后，要給粒子一定的熱運(yùn)動(dòng)速度，再給予分子一有梯度的宏觀速度。讓其運(yùn)動(dòng)達(dá)到穩(wěn)定后，再來(lái)計(jì)算分子的無(wú)規(guī)行走曲線。研究選用粒子宏觀的運(yùn)動(dòng)速度為x軸正方向的。

在y軸和z軸方向，每個(gè)粒子在某一個(gè)時(shí)刻的無(wú)規(guī)行走位移用以下公式計(jì)算：

(2)

其中，X1為在某一時(shí)刻某一粒子的位置，X2為該粒子初始時(shí)刻的位置。N為粒子總數(shù)。無(wú)規(guī)行走是所有粒子的運(yùn)動(dòng)位移平方相加取平均值。

畫(huà)出三個(gè)典型的速度梯度無(wú)規(guī)行走曲線與時(shí)間關(guān)系圖。由于y、z方向無(wú)規(guī)行走曲線類(lèi)似，都為直線，現(xiàn)只作出z方向無(wú)規(guī)行走曲線。如下圖所示為運(yùn)行106步的三種不同速度梯度Couette流動(dòng)的z方向無(wú)規(guī)行走曲線。

圖2 無(wú)規(guī)行走位移的方均值與時(shí)間的關(guān)系

利用最小二乘法擬合不同速度梯度下的無(wú)規(guī)行走曲線。

圖3 無(wú)規(guī)行走曲線斜率與速度梯度的關(guān)系

從圖(3)中可以看出，y、z方向的無(wú)規(guī)行走曲線斜率與速度梯度近似呈直線關(guān)系。隨著速度梯度的增加，y、z方向無(wú)規(guī)行走曲線斜率線性減小。無(wú)規(guī)行走曲線斜率亦即上文提到的2D，D為擴(kuò)散系數(shù)。這是因?yàn)榱Ｗ釉谶\(yùn)動(dòng)時(shí)，因?yàn)樗俣忍荻鹊脑龃?，粒子從速度較低層向速度較高層運(yùn)動(dòng)時(shí)，速度較高層粒子運(yùn)動(dòng)速度快，沿著x方向單位時(shí)間掃掠過(guò)的空間增大，更有可能與擴(kuò)散到此層的低速層粒子發(fā)生碰撞；反之，高速層粒子擴(kuò)散到低速度層，沿著x方向單位時(shí)間掃掠的空間更大，更有可能與低速層的粒子發(fā)生碰撞。這樣就會(huì)導(dǎo)致粒子在y、z兩個(gè)方向碰撞增多，所以平均位移減小。而同一個(gè)剪切率的y軸無(wú)規(guī)行走曲線與z軸無(wú)規(guī)行走曲線斜率的差別不大。因?yàn)檎{(diào)溫的作用，y方向與z方向的速度保持原來(lái)的量級(jí)不變，還是原來(lái)的熱運(yùn)動(dòng)速度，并且二者都是隨機(jī)的粒子向y軸運(yùn)動(dòng)就會(huì)以同樣的概率向著z軸運(yùn)動(dòng)，二者是協(xié)調(diào)相關(guān)的。由于空間的對(duì)稱性，可以證明沒(méi)有速度梯度的時(shí)候，稀薄氣體具有各項(xiàng)同性。而有速度梯度時(shí)這種各向同性趨勢(shì)被打破，稀薄氣體具有各項(xiàng)異性。從圖中即可看出隨著剪切率的增加，y方向無(wú)規(guī)行走曲線斜率比z方向的下降的更快。至于為什么會(huì)出現(xiàn)y方向無(wú)規(guī)行走曲線斜率比z方向的下降的更快的結(jié)果，有待進(jìn)一步去研究探索。

這樣來(lái)說(shuō)，粒子y軸和z軸粒子仍然是以熱運(yùn)動(dòng)速度在運(yùn)動(dòng)，而x軸的無(wú)規(guī)行走貌似不是很容易定義?？梢詮臏y(cè)量粒子無(wú)規(guī)行走位移的角度來(lái)定義無(wú)規(guī)行走速度，故需要在剪切流動(dòng)中建立速度運(yùn)動(dòng)坐標(biāo)系。本文定義的x方向無(wú)規(guī)行走速度是相當(dāng)于每個(gè)粒子每個(gè)時(shí)刻x軸方向的絕對(duì)速度在初始宏觀速度坐標(biāo)系下的相對(duì)速度。即公式(3)：

定義無(wú)規(guī)行走速度Vb(式3)為初始速度V0下的相對(duì)速度

Vb=Va-V0

(3)

Vb表示粒子在此時(shí)的無(wú)規(guī)行走速度，Va為粒子的絕度速度,V0為粒子在初始時(shí)刻的x方向宏觀運(yùn)動(dòng)運(yùn)動(dòng)速度。實(shí)際上粒子的無(wú)規(guī)行走速度相當(dāng)絕對(duì)速度在初始宏觀速度坐標(biāo)系下的相對(duì)速度。

從速度梯度0/40,到速度梯度5/40，每增加0.5/40的速度梯度計(jì)算無(wú)規(guī)行走，如下圖所示為運(yùn)行106步的速度梯度γ=1/40的Couette流動(dòng)的x方向無(wú)規(guī)行走曲線。

圖(4)B即是圖(4)A橫縱坐標(biāo)軸同時(shí)取對(duì)數(shù)的結(jié)果。從圖(4)A中可以看出，x方向無(wú)規(guī)行走曲線斜率隨著速度梯度的增大而增大，因?yàn)槎x的x方向無(wú)規(guī)行走速度是相當(dāng)絕對(duì)速度在初始宏觀速度坐標(biāo)系下的相對(duì)速度，當(dāng)粒子運(yùn)動(dòng)到宏觀速度更高的流層時(shí)，由于碰撞，它將獲得平均速度增量Δv=γz。這樣會(huì)導(dǎo)致x方向無(wú)規(guī)行走速度變大，因此無(wú)規(guī)行走位移增加。圖(4)B為圖(4)A在橫縱坐標(biāo)軸取對(duì)數(shù)之后的結(jié)果?？梢钥闯鰣D(4)B的為兩段直線。第一段直線的原因是由于流動(dòng)沒(méi)有充分發(fā)展，碰撞次數(shù)不夠多。位移與時(shí)間成線性關(guān)系，因此三種不同速度梯度的曲線嚴(yán)格重合,而且第一段直線的運(yùn)算時(shí)間只是在104的量級(jí),第二段直線是文中主要研究的,對(duì)第二段直線進(jìn)行最小二乘擬合，擬合出氣體分子無(wú)規(guī)行走位移的方均值與時(shí)間t的關(guān)系是：

(4)

(5)

(6)

由于z方向具有速度梯度，當(dāng)粒子從x方向宏觀速度較低層向x方向較高層運(yùn)動(dòng)時(shí)，由于粒子之間的相互碰撞，假定粒子獲得了Δv=γz的速度，由此式帶入得到(6)：

(7)

(8)

將(8)式兩邊平方得到

(9)

此式中的γ2* 2 *Dz即為上文中的系數(shù)a，可見(jiàn)它并不單單與擴(kuò)散系數(shù)有關(guān)。

在不同速度梯度下計(jì)算a的值，用tecplot作圖并作出最佳擬合曲線：

圖5 無(wú)規(guī)行走曲線斜率與速度梯度的關(guān)系

從圖(5)中可以看出，隨著速度梯度的增大，比例系數(shù)a值呈現(xiàn)非線性增加。至于為什么會(huì)出現(xiàn)這樣的結(jié)果有待進(jìn)一步去探索研究。

3 結(jié) 論

通過(guò)分子動(dòng)力學(xué)模擬的方法來(lái)研究不同均勻剪切流與分子無(wú)規(guī)行走曲線的斜率的關(guān)系。應(yīng)用分子動(dòng)力學(xué)的模擬方法，模擬出了均勻剪切流，利用分子無(wú)規(guī)行走曲線計(jì)算公式，計(jì)算出不同均勻剪切流的無(wú)規(guī)行走曲線，計(jì)算無(wú)規(guī)行走擬合曲線斜率。計(jì)算結(jié)果表明，在相同時(shí)間條件下，速度梯度越大，x軸方向無(wú)規(guī)行走曲線斜率增大，y軸z軸兩個(gè)方向的無(wú)規(guī)行走曲線斜率越小，速度梯度越大對(duì)無(wú)規(guī)行走位移壓縮越嚴(yán)重。而同一剪切率下的y軸無(wú)規(guī)行走曲線與z軸無(wú)規(guī)行走曲線斜率的差別不大。在x軸方向分子無(wú)規(guī)行走位移的方均值與時(shí)間t的三次方成正比，比例系數(shù)a與分子擴(kuò)算系數(shù)Dz以及速度梯度γ有關(guān)，并隨著速度梯度增大呈現(xiàn)非線性趨勢(shì)；在y軸、z軸方向分子無(wú)規(guī)行走的方均值與時(shí)間t成正比，其比例系數(shù)即擴(kuò)散系數(shù)的二倍即2*D隨著速度梯度增加而線性減小。