冷岳峰 王嵩博 劉運(yùn)宇 李 強(qiáng)

(遼寧工程技術(shù)大學(xué)機(jī)械工程學(xué)院 遼寧阜新 123000)

分子動力學(xué)(MD)模擬方法是研究熱輸運(yùn)特性最有效的工具之一，特別是對于基于晶格動力學(xué)方法的復(fù)雜結(jié)構(gòu)，其計算能力非常強(qiáng)大。LAMMPS是一種可對處于三態(tài)的粒子團(tuán)進(jìn)行建模的分子動力學(xué)代碼，它可以使用各種原子間勢和邊界條件對原子、聚合物、金屬、陶瓷、氧化物、顆?；蚝暧^系統(tǒng)進(jìn)行建模[1]。

原子結(jié)構(gòu)是影響金剛石、石墨烯、碳納米管等重要的因素，利用原子間原子勢的原子模擬方法是研究碳原子結(jié)構(gòu)最有效的方法之一。但必須要有可靠的原子間勢才能研究發(fā)現(xiàn)各種基本物理性質(zhì)，如結(jié)構(gòu)性質(zhì)、表面特性、熱特性等?，F(xiàn)有的文獻(xiàn)表明[2-3]，納米顆粒的含量、分散納米顆粒的大小和形狀、納米流體的溫度以及分散納米顆粒的分布等因素都會影響納米流體的熱性能。這些因素中的每一個都是納米流體傳熱現(xiàn)象廣泛研究的來源，從而促進(jìn)了該領(lǐng)域知識的發(fā)展。

納米流體的概念最早是由DAI[4]在20世紀(jì)90年代中期提出的，作為一種替代微米級固體的方法，這些固體被用于傳統(tǒng)的冷卻劑中以提高其導(dǎo)熱性。研究人員致力于進(jìn)一步了解支配納米流體的物理學(xué)，包括使用不同的方法計算各種類型納米流體的熱導(dǎo)率，以及研究熱的機(jī)制。與傳統(tǒng)流體相比，納米流體具有更高的導(dǎo)熱系數(shù)、更好的對流換熱以及更低的壓降，這使得納米流體成為一種具有良好傳熱潛力的新技術(shù)。

SUKKAR等[5]將球形氧化銅納米粒子和氧化鈦粒子添加到潤滑油中，研究了不同質(zhì)量分?jǐn)?shù)的納米粒子在潤滑油中的導(dǎo)熱性，結(jié)果表明，納米氧化銅顆粒對潤滑油導(dǎo)熱性能的提升優(yōu)于氧化鈦；當(dāng)納米粒子質(zhì)量分?jǐn)?shù)為0.1%時，含納米氧化銅和氧化鈦潤滑油的導(dǎo)熱系數(shù)分別增加了7.27%和4.54%。XIE等[6]研究了碳納米管體積分?jǐn)?shù)對納米流體熱導(dǎo)率的影響，發(fā)現(xiàn)碳納米管使納米流體導(dǎo)熱系數(shù)提高了7%～20%。ASSAEL等[7]研究了在2種不同分散劑作用下碳納米管水懸浮液的熱導(dǎo)率，發(fā)現(xiàn)十六烷基三甲基溴化銨(CTAB)表面活性劑使得懸浮液導(dǎo)熱系數(shù)增長了8%。HWANG等[8]研究了質(zhì)量分?jǐn)?shù)0.5%碳納米管在礦物油基底液中的熱導(dǎo)率，結(jié)果表明潤滑油的熱導(dǎo)率增加了9%。GARG等[9]采用超聲波研究了碳納米管對水基納米流體傳熱性能的影響，發(fā)現(xiàn)質(zhì)量分?jǐn)?shù)1%碳納米管可使水基納米流體的導(dǎo)熱系數(shù)提高20%。SINGH等[10]研究了碳納米管對水-乙二醇納米流體熱導(dǎo)率的影響，發(fā)現(xiàn)質(zhì)量分?jǐn)?shù)0.4%的碳納米管使得納米流體熱導(dǎo)率提高了72%。

李春風(fēng)等[11]對富勒烯、石墨烯、金剛石、碳納米管和納米石墨的摩擦學(xué)性能進(jìn)行了歸納總結(jié)，發(fā)現(xiàn)1%的富勒烯便可有效降低石蠟油摩擦因數(shù)33%；含5%富勒烯的石蠟油潤滑下鋼球磨斑直徑可降低33%，摩擦因數(shù)降低了20%。ZHANG等[12]研究了棕櫚油、大豆油、菜籽油分別作為MQL基礎(chǔ)油的潤滑性能，證明了3種植物油的潤滑效果差異很小，摩擦因數(shù)只有0.04的差異，與液體石蠟相比都具有較低的摩擦因數(shù)。但是大豆油相比其他2種植物油，具有較低的黏度，更適合作為NMQL的基礎(chǔ)油。

目前的文獻(xiàn)都是針對一些水為基礎(chǔ)液進(jìn)行的分子動力學(xué)仿真，對基礎(chǔ)油液加入金剛石納米顆粒及其他碳結(jié)構(gòu)的納米顆粒研究較少。本文作者以蓖麻油酸為基礎(chǔ)潤滑劑，以金剛石納米顆粒為添加劑，采用Tersoff勢描述金剛石原子之間的勢函數(shù)，進(jìn)行了分子動力學(xué)仿真，研究了加入不同體積分?jǐn)?shù)和不同粒徑的納米顆粒對潤滑性能以及熱導(dǎo)率的影響。

1 研究方法和數(shù)學(xué)模型

1.1 研究方法及理論基礎(chǔ)

熱傳導(dǎo)是溫度梯度引起能量傳遞的現(xiàn)象，屬于非平衡過程。另一種計算熱導(dǎo)率的方法就是基于非平衡分子動力學(xué)的模擬[13]。在該方法中，首先通過某種方式建立一個溫度梯度，并等待足夠長的時間，使得系統(tǒng)達(dá)到穩(wěn)態(tài)。在達(dá)到穩(wěn)態(tài)后，對系統(tǒng)的溫度梯度和非平衡穩(wěn)態(tài)熱流進(jìn)行測量，然后就可以根據(jù)傅里葉定律來計算熱導(dǎo)率。

在此非平衡態(tài)研究熱導(dǎo)率方法中，系統(tǒng)的總能量近似為

(1)

在參數(shù)優(yōu)化的過程中，需要有正確的生成能和晶格常數(shù)。然而，僅用目前的方法很難得到具有正確晶格參數(shù)值的碳結(jié)構(gòu)。此外，考慮到金剛石鍵合特性各不相同，試圖用一種常見的勢函數(shù)來描述內(nèi)部原子的相互作用是不可取的。因此，有必要把金剛石和外部油液的相互作用分開，用不同的形式來描述它們。因此，最終決定用Tersoff形式描述金剛石原子之間的相互作用；用Lennard-Jones勢來描述層間相互作用。根據(jù)Lennard-Jones作用勢，原子間的相互作用[14]如下:

(2)

式中：εij表示力的強(qiáng)度的參數(shù)；σij表示原子大小的參數(shù)；r表示原子之間的距離。

基于文獻(xiàn)[3,15]的研究結(jié)果，文中熱導(dǎo)率的研究方法運(yùn)用了速度重標(biāo)法。當(dāng)系統(tǒng)達(dá)到穩(wěn)態(tài)之后，開始逐漸地增加熱源中粒子的速率，并且對熱源中的每一個粒子i，作如下速度轉(zhuǎn)換

(3)

式中：vc是熱源粒子的質(zhì)心速度，

(4)

(5)

Δε是一個常數(shù)，該常數(shù)等于熱源中能量的增量，

(6)

如果每 Δt作一次速度變換，則產(chǎn)生的熱流密度大小為

(7)

其中，A是模擬系統(tǒng)在垂直于傳導(dǎo)方向的橫截面積。當(dāng)使用周期性邊界條件時，應(yīng)當(dāng)將該熱流除以 2，因為當(dāng)在某一位置注入熱量時，熱量將從熱源區(qū)域兩邊流出，每一邊流出的能量是總量的1/2。在系統(tǒng)達(dá)到平衡穩(wěn)態(tài)之后，開始測量溫度梯度，然后根據(jù)傅里葉定律計算非平衡態(tài)熱導(dǎo)率。

1.2 物理模型和參數(shù)配置

圖1給出了以蓖麻油酸為潤滑劑的分子動力學(xué)模型，模型分為上下壁面，上壁面為金剛石，下壁面為金屬鎳，中間為蓖麻油酸分子以及內(nèi)部的納米金剛石顆粒?，F(xiàn)以下壁面Ni原子層固定不動，上壁面的金剛石層以10 m/s的速度向右側(cè)移動，同時對上壁面施加一個穩(wěn)定的壓力，以10 m/s的速度壓縮，對摩擦過程進(jìn)行模擬，用來模擬壓縮和滑移時的過程。

在上下兩層壓縮面之間建立蓖麻油酸分子團(tuán)尺寸為20 nm×10 nm×10 nm的基礎(chǔ)流體，上下壁面為2 nm。金剛石納米粒子的直徑為4 nm，建立的模型總的原子數(shù)為22 584個。

模擬計算采用聯(lián)合力場勢函數(shù)，在x、y方向上采用周期性邊界條件，在z方向上采用非周期性邊界條件。在分子動力學(xué)模擬中，在有限溫度下進(jìn)行的實驗中經(jīng)常發(fā)現(xiàn)，金剛石結(jié)構(gòu)在高溫破壞或轉(zhuǎn)變?yōu)槠渌Y(jié)構(gòu)之前應(yīng)該是穩(wěn)定的。

圖1 低體積分?jǐn)?shù)金剛石納米顆粒潤滑油的分子動力學(xué)模型

在具有周期性邊界條件的直徑4 nm的球形金剛石結(jié)構(gòu)樣品上進(jìn)行了模擬。在模擬中，所有的性能計算都是在253 K左右溫度下進(jìn)行的。在253 K的條件下，經(jīng)過300 ps時間步后，進(jìn)行徑向分布函數(shù)(RDF)的分析。在以前的研究中已經(jīng)發(fā)現(xiàn)，如果一個結(jié)構(gòu)不是穩(wěn)定的，在幾千個MD步驟內(nèi)，原子結(jié)構(gòu)很快會發(fā)生崩潰。因為布朗運(yùn)動的作用，納米粒子周圍表現(xiàn)出較高的不穩(wěn)定性，納米粒子受到的布朗運(yùn)動更加劇烈，對減摩性能有較大影響。

模擬過程中需要進(jìn)行3個步驟：弛豫、壓縮和剪切。首先進(jìn)行弛豫，控制在253 K的溫度下，采用聯(lián)合力場，金剛石納米粒子采用Tersoff勢函數(shù)，鎳單質(zhì)采用EAM勢函數(shù)，截斷半徑選取7.0，充分弛豫。施加壓力，上壁面以10 m/s的速度壓向下壁面。采用周期性的邊界條件，模擬采用的計算方法為Verlet法，時間步長選定為0.1 fs。采用Nose-Hoover控溫的方法，在NVT系統(tǒng)下對所有的原子進(jìn)行弛豫，通過改變壓力與溫度以及納米顆粒的粒徑大小，文中重點考察不同條件下對薄膜潤滑的原子分布、滑移特性和分子團(tuán)聚特性的影響。

2 模擬驗證及其討論

2.1 數(shù)密度分析

當(dāng)系統(tǒng)中存在固液氣三態(tài)相交時，數(shù)密度可以很好地給出分界附近的結(jié)構(gòu)性特點。粒子數(shù)密度的不均勻性導(dǎo)致粒子的輸運(yùn)也就是擴(kuò)散現(xiàn)象。圖2所示為中心壁面區(qū)域油膜厚度的原子數(shù)密度分布，可以發(fā)現(xiàn)近壁面區(qū)域的原子呈層狀散布，這說明近壁面附近存在油膜吸附層。

圖2 中心壁面區(qū)域沿y軸的原子數(shù)密度分布

從圖3中可以看出，金剛石顆粒周圍存在吸附層，納米粒子周圍的密度較大，潤滑劑附著在金剛石納米顆粒上，當(dāng)納米顆粒被壓扁時，這層薄膜也未被破壞。由此，可以認(rèn)定金剛石顆粒周圍的吸附層會對潤滑劑分子產(chǎn)生影響，有可能會改變局部的密度等。另外，如果油膜相對比較厚，這種情況下的影響基本可以忽略，但對于比較薄的潤滑膜，這種影響還是比較大的，不能忽視。

圖3 壓縮過程中納米粒子的壓扁狀態(tài)和邊界油膜吸附層

從圖3中還可以看出,當(dāng)位移在2 nm左右時，分子運(yùn)動加快，大量的金剛石納米顆粒周圍分子被擠入前面吸附層。此時，金剛石納米顆粒與壁面間存在一層潤滑油分子?？梢酝茰y如果上下壁面進(jìn)一步靠近，納米顆粒將會與壁面直接接觸，此時，薄膜就會發(fā)生破裂，從而就會發(fā)生固體潤滑劑的作用。

2.2 徑向分布函數(shù)分析

徑向分布函數(shù)(RDF)是描述系統(tǒng)結(jié)構(gòu)特點非常有效的方法，可以表征固體或液體微觀結(jié)構(gòu)的無序化程度[1]。圖4所示為基礎(chǔ)油蓖麻油酸的徑向分布函數(shù)。

圖4 基礎(chǔ)油的徑向分布函數(shù)(d=4 nm)

圖5表示的是金剛石納米潤滑油分子團(tuán)的RDF對比，可以發(fā)現(xiàn)納米潤滑油與蓖麻油酸基礎(chǔ)油的RDF在2～4 nm之間有些許偏差，當(dāng)上壁面位移超出這部分范圍時，二者的曲線才逐漸趨于平穩(wěn)，這說明金剛石納米顆粒周圍稀疏的吸附層對油膜分子排列的緊密程度和邊界條件下潤滑膜強(qiáng)度具有一定的影響，但是影響不大。

當(dāng)油膜的厚度下降到某一程度后，油膜的各部位厚度就會發(fā)生轉(zhuǎn)變，轉(zhuǎn)變?yōu)椤吧舷卤诿嫖綄雍退闹芪綄右约敖饎偸{米顆粒”，三者之間的彼此作用提高了邊界潤滑膜的強(qiáng)度和承載力?？梢钥吹郊{米粒子并沒有明顯的變形，只是隨著油液一起向下運(yùn)動，此時，油液分層現(xiàn)象更加明顯，分層數(shù)也有一定的減少，單層油膜厚度也在增加，層數(shù)減小。但是，金剛石納米顆粒只能在一定載荷范圍內(nèi)提高邊界潤滑的承載能力，載荷超過這一范圍時邊界膜將會發(fā)生破裂，導(dǎo)致邊界潤滑膜的承載能力失效。

圖5 金剛石納米潤滑油分子團(tuán)與基礎(chǔ)油液徑向分布函數(shù)的對比(d=4 nm)

2.3 納米顆粒的角速度分析

如圖6 所示，改變周期性邊界條件中的顆粒數(shù)量和粒徑，使系統(tǒng)內(nèi)的顆粒體積分?jǐn)?shù)發(fā)生變化。當(dāng)顆粒體積分?jǐn)?shù)一定時，粒徑變小，潤滑膜承載力越高，金剛石顆粒的個數(shù)也會增加，納米顆粒對基礎(chǔ)潤滑油的布朗運(yùn)動效果也會加強(qiáng)，無規(guī)則運(yùn)動也會進(jìn)一步增強(qiáng)。

圖6 高體積分?jǐn)?shù)的納米顆粒潤滑油分子動力學(xué)模型

如圖7—9所示，可以發(fā)現(xiàn)，粒徑越小，納米粒子的轉(zhuǎn)動角速度越高，這個現(xiàn)象進(jìn)一步說明了納米粒子大小的擾動問題。此外，納米粒子的直徑擴(kuò)大，增大了粒子體積效應(yīng)的影響效果。因而，隨著粒徑的減小，納米粒子潤滑膜的承載能力增強(qiáng)，液固轉(zhuǎn)化的壓力變大，摩擦力降低。

圖7 不同粒徑納米顆粒繞x軸的轉(zhuǎn)動角速度分量

圖8 納米顆粒的線速度

圖9 納米顆粒繞x軸的轉(zhuǎn)動角速度

2.4 系統(tǒng)的熱導(dǎo)率

根據(jù)熱力學(xué)第二定律，一個系統(tǒng)在不受外力作用時，若其內(nèi)部有熱力學(xué)性質(zhì)的不均勻性，則它一定處于非平衡的狀態(tài)，并有向平衡態(tài)靠近的趨勢。這種由熱力學(xué)性質(zhì)的不均勻性導(dǎo)致的熱力學(xué)過程叫做輸運(yùn)過程。將一個系統(tǒng)置于2個溫度不同的熱源之間，最終會在系統(tǒng)內(nèi)建立一個穩(wěn)定的溫度分布?？烧J(rèn)為這樣的系統(tǒng)處于一個穩(wěn)態(tài)，但不處于一個平衡態(tài)。

如圖10所示為非平衡態(tài)熱傳導(dǎo)模擬，在非平衡態(tài)熱導(dǎo)率(NEMD)方法中，首先建立一個溫度梯度，并弛豫足夠長的時間，使得系統(tǒng)達(dá)到穩(wěn)態(tài)。在達(dá)到穩(wěn)態(tài)后，對系統(tǒng)的溫度梯度和非平衡穩(wěn)態(tài)熱流進(jìn)行測量，然后就可以根據(jù)傅里葉定律來計算熱導(dǎo)率。將模擬系統(tǒng)在熱傳導(dǎo)方向分割為若干相等塊，將其中的兩塊分別設(shè)置為熱源區(qū)域和熱匯區(qū)域。采用周期邊界條件，塊的個數(shù)一般取為偶數(shù)2N，并可以取第1塊為熱源，取第N+ 1塊為熱匯，或者反過來。

圖10 非平衡態(tài)熱傳導(dǎo)模擬

在基礎(chǔ)液油中，導(dǎo)熱系數(shù)可以運(yùn)用導(dǎo)熱儀測得，如DTC-300，僅需比較小的樣品，就可以通過比較小的容器測得，薄膜也可以使用多層技術(shù)準(zhǔn)確得到測量值。但是在加入納米金剛石顆粒之后，混合油模擬過程中的瞬時熱導(dǎo)率只能通過計算機(jī)模擬得到。根據(jù)圖11所示的體系溫度分布，在體系的兩側(cè)添加熱源并且繪制初始狀態(tài)的溫度分布情況。圖12示出了模擬得到的熱導(dǎo)率隨模擬時間點的變化，熱導(dǎo)率會在其準(zhǔn)確值上下進(jìn)行浮動，最后進(jìn)行數(shù)據(jù)處理可以得到相對準(zhǔn)確的數(shù)據(jù)。

圖13所示是253 K溫度下，納米顆粒體積分?jǐn)?shù)分別為0、0.6%、1.2%、1.8%和2.4%時潤滑油的熱導(dǎo)率。結(jié)果表明，納米粒子體積分?jǐn)?shù)的提高會極大地影響潤滑油的熱傳導(dǎo)，使納米流體的熱導(dǎo)率呈近似線性增加。其主要原因是，納米粒子的傳熱性大于液體分子，并且每一個加入的納米粒子對液體都會有較大的擾動，從而增強(qiáng)了納米流體的熱導(dǎo)率。因此，在保證納米流體液體性質(zhì)的前提下，適當(dāng)增大納米粒子體積分?jǐn)?shù)可以有效提高納米流體的熱導(dǎo)率。

圖11 系統(tǒng)局部溫度熱流

圖12 熱導(dǎo)率隨模擬時間點的變化

圖13 不同體積分?jǐn)?shù)納米顆粒潤滑油的熱導(dǎo)率

3 結(jié)論

(1)納米粒子加入基礎(chǔ)油液中，會與基礎(chǔ)油之間發(fā)生吸附作用，在納米粒子表面形成吸附層，當(dāng)載荷在納米管潤滑膜的可承受范圍內(nèi)時是可以存在邊界潤滑的。納米粒子在布朗運(yùn)動的作用下也會產(chǎn)生能量交換，促進(jìn)納米流體液體間的能量傳遞。因為布朗運(yùn)動的作用，納米粒子周圍表現(xiàn)出較高的不穩(wěn)定性。納米流體與基礎(chǔ)流體相比有更好的減摩性能，納米粒子受到的布朗運(yùn)動更加劇烈，對減摩性能有較大影響。

(2)在形成邊界油膜后，金剛石納米顆?？梢蕴岣呋A(chǔ)油液的承載能力，但是由于壁面和納米顆粒周圍的吸附層在壓縮過程中不斷靠近，納米顆粒在邊界油膜被破壞之前就發(fā)生了形變，在此之前，邊界油膜能很好地保護(hù)納米顆粒。

(3)當(dāng)系統(tǒng)中的納米顆粒粒徑發(fā)生變化時，在壓縮和滑動過程中，納米顆粒的轉(zhuǎn)動角速度會發(fā)生變化，粒徑越小，轉(zhuǎn)動速度就越高；反之，轉(zhuǎn)動速度就越小，相應(yīng)地摩擦力就會發(fā)生變化。

(4)利用非平衡態(tài)熱導(dǎo)率的方法研究了納米金剛石潤滑油的熱導(dǎo)率，發(fā)現(xiàn)納米顆粒的加入會使得基礎(chǔ)油液熱導(dǎo)率提高，隨納米粒子體積分?jǐn)?shù)的增加，納米流體的熱導(dǎo)率呈近似線性增加。

(5)文中研究方法也可用于其他系統(tǒng)的傳熱研究，在獲得穩(wěn)定的溫度分布后，通過線性擬合的方法來確定溫度梯度。所有的模擬都在熱傳輸?shù)木€性響應(yīng)范圍內(nèi)，進(jìn)一步驗證了傅里葉定律求得NEMD熱導(dǎo)率方法的可行性。在模擬實驗中可以發(fā)現(xiàn)，在比較長的樣品模擬中，溫度梯度應(yīng)該很小，否則靠近熱源的溫度將偏離目標(biāo)溫度，導(dǎo)致溫度分布非線性。