李 龍，王瑞金，朱澤飛

(杭州電子科技大學(xué)機(jī)械工程學(xué)院，浙江 杭州 310018)

0 引 言

多粒子碰撞動力學(xué)是一種粗粒化分子模型的方法[1]，結(jié)合了分子動力學(xué)和格子-玻爾茲曼方法(Lattice-Boltzmann Method，LBM)的優(yōu)點，能夠解決流體的熱漲落問題，計算量小，可以模擬復(fù)雜邊界條件的體系[2]。很多學(xué)者使用MPCD方法進(jìn)行研究，如陳文多等[3]采用MPCD方法研究了環(huán)形高分子的流體動力學(xué)并計算其回轉(zhuǎn)半徑，Akhter等[4]采用MPCD方法對微壓縮低雷諾數(shù)條件下的顆粒流體進(jìn)行數(shù)值模擬，Wysocki等[5]采用MPCD方法研究膠體系統(tǒng)在微米尺度狹縫中的水動力學(xué)不穩(wěn)定性，Eggersdorfer等[6]采用MPCD方法進(jìn)行流體中大顆粒聚集結(jié)構(gòu)的研究。本文采用MPCD方法分析了影響流體導(dǎo)熱系數(shù)計算的因素，并計算AR體系和水分子體系的導(dǎo)熱系數(shù)。

1 MPCD方法

流體粗粒子采用MPCD進(jìn)行數(shù)值模擬，一般分為流動步和碰撞步。在流動步中，粒子的運(yùn)動速度不變且可得到粗粒子的新位移。在碰撞步中，同一個盒子中的粗粒子速度矢量繞著隨機(jī)軸旋轉(zhuǎn)一定的角度，得到粗粒子的新速度。流動過程中粒子屬于平動，滿足：

ri(t+Δt)=ri(t)+vi(t)Δt

(1)

粗粒子速度計算公式如下：

vi(t+Δt)=vcm(t)+Ω(α)[vi(t)-vcm(t)]

(2)

式中，ri(t)是粒子的位移矢量，vi是粒子速度矢量，Δt是時間步長，α是粒子速度矢量繞任意隨機(jī)軸旋轉(zhuǎn)的角度，Ω(α)是隨機(jī)旋轉(zhuǎn)矩陣。vcm(t)是每個盒子的質(zhì)心速度，質(zhì)心速度的計算公式如下：

(3)

式中，N是在一個盒子中的粒子數(shù)，M是粗粒子的質(zhì)量。

2 導(dǎo)熱系數(shù)的計算

(4)

(5)

(6)

3 模擬體系的大小及盒子劃分

本文選擇分子動力學(xué)中經(jīng)常使用的LJ約化單位，故本文中的參數(shù)基本為無量綱參數(shù)。

選擇模擬體系的尺寸為33.72×33.72×33.72，長度D=d/σ，d是標(biāo)準(zhǔn)單位下的長度，σ是LJ勢函數(shù)中的單位距離參數(shù)，整個體系劃分為20×20×20個盒子如圖1所示，流體粗粒子體系如圖2所示。

圖1 體系盒子劃分

圖2 模擬體系

4 影響導(dǎo)熱系數(shù)計算結(jié)果的因素

4.1 流體粗粒子質(zhì)量

選取體系溫度T*=0.71，劃分盒子的尺寸為1.700，晶格常數(shù)f=1.55，整個模擬體系含有62 500個粗?；Ｗ樱總€盒子內(nèi)含有的平均粒子數(shù)為9.112，粒子質(zhì)量M分別為0.40，0.60，0.80，1.00，1.20。為了得出不同SRD時間步長下的流體導(dǎo)熱系數(shù)，tSRD分別取0.25，0.30，0.35，0.40。導(dǎo)熱系數(shù)隨著粗粒子質(zhì)量的變化情況如圖3所示。從圖3可以看出，導(dǎo)熱系數(shù)隨粗粒子質(zhì)量的增大而減小。溫度一定的情況下，質(zhì)量越小，粒子的速度越大，所以粒子間碰撞幾率增大。從圖3還可以看出，導(dǎo)熱系數(shù)隨SRD時間步長的增大而增大。SRD時間步長越大，粒子移動的距離越遠(yuǎn)，不同盒子內(nèi)粒子的動量和能量交換就更加頻繁，間接增強(qiáng)了流體的傳熱。

圖3 導(dǎo)熱系數(shù)隨著粗粒子質(zhì)量的變化情況

4.2 盒子尺寸

選取參數(shù)粒子質(zhì)量M=0.40，tSRD=0.35，溫度T*=0.71，晶格常數(shù)f=1.55。計算體系在不同盒子尺寸和晶格常數(shù)下的導(dǎo)熱系數(shù)。不同尺寸盒子內(nèi)的平均粒子數(shù)不同，影響導(dǎo)熱系數(shù)的計算結(jié)果。本文選取如表1所示的體系參數(shù)進(jìn)行模擬，粒子密度ρ是指每個盒子內(nèi)的平均粒子數(shù)，f采用的是LJ約化單位，盒子尺寸和f均是無量綱參數(shù)。標(biāo)準(zhǔn)單位下的晶格常數(shù)L與f成反比，即f越大，每個盒子內(nèi)的平均粒子數(shù)就越多。在f不變的情況下，ρ隨著盒子尺寸的增大而增大。

表1 不同f和盒子尺寸下的粒子密度ρ

導(dǎo)熱系數(shù)隨著盒子尺寸的變化情況如圖4所示。從圖4可以看出，在f相同的情況下，導(dǎo)熱系數(shù)隨著盒子尺寸的增大而減小，雖然盒子尺寸增大使得每個盒子的粒子數(shù)增加，但粒子移出盒子的概率減少了，導(dǎo)致粒子間的自相關(guān)性增大，粒子不能和更遠(yuǎn)的盒子中的粒子發(fā)生碰撞。

4.3 體系溫度

選取參數(shù)M=1.00，tSRD=0.35，體系的大小為33.72×33.72×33.72，盒子尺寸為1.780，體系溫度T*分別為0.50，0.70，1.00，1.20，1.50，另外，晶格常數(shù)f分別為1.55，1.75，1.95。導(dǎo)熱系數(shù)隨著體系溫度的變化情況如圖5所示。從圖5可以看出，導(dǎo)熱系數(shù)隨著溫度的增大而增大，在體系粒子密度相同的情況下，溫度越高，粗粒子的動能越大，粒子的碰撞幾率就高，導(dǎo)熱系數(shù)越大。但是，一般情況下，大部分液體的導(dǎo)熱系數(shù)隨著溫度的增大而減小，因為本文是在溫度不同且粒子密度相同的情況下進(jìn)行模擬實驗的，并沒有考慮溫度對分子間距離的影響，即沒有考慮溫度對體系粒子密度的影響。大部分液體的導(dǎo)熱系數(shù)是隨著體系密度的減小而減小。溫度越大，體系的粒子密度是越小的，晶格常數(shù)與粒子密度成正比。從圖5還可以看出，T*=0.70，f=1.95時的導(dǎo)熱系數(shù)大于T*=1.00，f=1.55時的導(dǎo)熱系數(shù)，雖然溫度增加了，但是粒子密度卻降低了，導(dǎo)熱系數(shù)也會下降。

4.4 粒子旋轉(zhuǎn)角度

選取參數(shù)T*=0.71，M=0.40，盒子尺寸為1.780，f=1.95，選擇旋轉(zhuǎn)固定角度分別為為90°，270°，360°。為了準(zhǔn)確計算體系的導(dǎo)熱系數(shù)，使每個盒子內(nèi)的粗粒子在同一時間步內(nèi)旋轉(zhuǎn)不同的角度。

粒子以不同概率旋轉(zhuǎn)不同角度下的體系導(dǎo)熱系數(shù)計算結(jié)果如表2所示。旋轉(zhuǎn)360°意味著在一個SRD時間步內(nèi)速度矢量是不變的，這和增大SRD時間步長所起到的效果相同。從表2可以看出，旋轉(zhuǎn)360°的概率占比越大，導(dǎo)熱系數(shù)越大，但概率占比不能為1，否則體系會細(xì)?；?，計算結(jié)果毫無意義。

圖4 導(dǎo)熱系數(shù)隨著盒子尺寸的變化情況

圖5 導(dǎo)熱系數(shù)隨著體系溫度的變化情況

表2 粒子以不同概率旋轉(zhuǎn)不同角度下的體系導(dǎo)熱系數(shù)計算結(jié)果

5 Ar原子和水分子體系的導(dǎo)熱系數(shù)計算

在計算體系的導(dǎo)熱系數(shù)時，溫度通常是已定的。本文通過計算Ar原子和水分子體系的導(dǎo)熱系數(shù)來驗證MPCD方法是否可以用于流體導(dǎo)熱系數(shù)的計算。

選取參數(shù)M=1.00，T*=0.71(轉(zhuǎn)化為標(biāo)準(zhǔn)單位為86 K)，f=1.75。整個體系的大小為33.72×33.72×33.72，密度為1 406 kg/m3，盒子尺寸取1.700，體系含有32 000個粒子，整個體系劃分成了20×20×20個盒子。使Ar的粗?；Ｗ釉诿恳粋€時間步內(nèi)以不同的概率旋轉(zhuǎn)不同的角度(90°，270°，360°)，導(dǎo)熱系數(shù)隨著SRD時間步長的變化情況如圖6所示。當(dāng)粗粒子每一個時間步旋轉(zhuǎn)90°，270°，360°的概率分別是1/6，1/6，4/6，且時間步長tSRD=1.00，T*=0.71時，模擬計算的導(dǎo)熱系數(shù)為0.136 5 W·(m·K)-1，與文獻(xiàn)[8]得到的0.132 0 W·(m·K)-1較接近，誤差為3.4%。Ar原子體系導(dǎo)熱系數(shù)隨著時間的變化情況如圖7所示。

圖6 導(dǎo)熱系數(shù)隨著時間步長的變化情況

圖7 Ar原子體系導(dǎo)熱系數(shù)隨時間的變化情況

另外，本文還計算了T*=1.00時的導(dǎo)熱系數(shù)，Ar體系中粒子數(shù)變?yōu)?3 328，因為溫度變大，粒子間的實際晶格常數(shù)變大，體系中所含的粒子數(shù)就減少了。選取盒子尺寸為1.983，通過增大盒子尺寸，可以保證在體系溫度不同時的盒子所含的平均粒子數(shù)一樣。本次模擬的計算結(jié)果為0.101 6 W·(m·K)-1。

同理，在模擬水分子體系時，一個水分子粗粒化成一個粗粒子，模擬體系尺寸為33.72×33.72×33.72，模擬參數(shù)選取M=0.45(分子量為18)，T*=2.44(T=298 K)。導(dǎo)熱系數(shù)計算結(jié)果為0.599 0 W·(m·K)-1，與實驗值0.608 4 W·(m·K)-1相比[9]，誤差是1.5%。

本文進(jìn)行模擬計算時，使用的是8核的Intel-i7服務(wù)器，僅用3 min左右即可完成Ar原子或水分子體系的運(yùn)動模擬和導(dǎo)熱系數(shù)的計算。相同條件下，采用MD方法，大約需要約6 h[10]，計算量大大減少。

6 結(jié)束語

本文使用MPCD方法模擬流體的運(yùn)動分析了影響流體導(dǎo)熱系數(shù)計算的因素，并計算Ar體系和水分子體系的導(dǎo)熱系數(shù)。相比MD方法，計算效率得到較大提高。但是，導(dǎo)熱系數(shù)的計算易受各種參數(shù)的影響，故需要合理選取體系參數(shù)。若要預(yù)測未知流體的導(dǎo)熱系數(shù)，需要先研究各種影響因素與導(dǎo)熱系數(shù)的關(guān)系，這是本文后期研究的重點。