基于MPCD方法計(jì)算流體的導(dǎo)熱系數(shù)

李 龍，王瑞金，朱澤飛

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

0 引 言

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

1 MPCD方法

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

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

(1)

粗粒子速度計(jì)算公式如下：

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

(2)

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

(3)

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

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

(4)

(5)

(6)

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

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

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

圖1 體系盒子劃分

圖2 模擬體系

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

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

選取體系溫度T*=0.71，劃分盒子的尺寸為1.700，晶格常數(shù)f=1.55，整個(gè)模擬體系含有62 500個(gè)粗粒化粒子，每個(gè)盒子內(nèi)含有的平均粒子數(shù)為9.112，粒子質(zhì)量M分別為0.40，0.60，0.80，1.00，1.20。為了得出不同SRD時(shí)間步長(zhǎng)下的流體導(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時(shí)間步長(zhǎng)的增大而增大。SRD時(shí)間步長(zhǎng)越大，粒子移動(dòng)的距離越遠(yuǎn)，不同盒子內(nèi)粒子的動(dòng)量和能量交換就更加頻繁，間接增強(qiáng)了流體的傳熱。

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

4.2 盒子尺寸

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

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

導(dǎo)熱系數(shù)隨著盒子尺寸的變化情況如圖4所示。從圖4可以看出，在f相同的情況下，導(dǎo)熱系數(shù)隨著盒子尺寸的增大而減小，雖然盒子尺寸增大使得每個(gè)盒子的粒子數(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òng)能越大，粒子的碰撞幾率就高，導(dǎo)熱系數(shù)越大。但是，一般情況下，大部分液體的導(dǎo)熱系數(shù)隨著溫度的增大而減小，因?yàn)楸疚氖窃跍囟炔煌伊Ｗ用芏认嗤那闆r下進(jìn)行模擬實(shí)驗(yàn)的，并沒(méi)有考慮溫度對(duì)分子間距離的影響，即沒(méi)有考慮溫度對(duì)體系粒子密度的影響。大部分液體的導(dǎo)熱系數(shù)是隨著體系密度的減小而減小。溫度越大，體系的粒子密度是越小的，晶格常數(shù)與粒子密度成正比。從圖5還可以看出，T*=0.70，f=1.95時(shí)的導(dǎo)熱系數(shù)大于T*=1.00，f=1.55時(shí)的導(dǎo)熱系數(shù)，雖然溫度增加了，但是粒子密度卻降低了，導(dǎo)熱系數(shù)也會(huì)下降。

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

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

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

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

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

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

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

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

選取參數(shù)M=1.00，T*=0.71(轉(zhuǎn)化為標(biāo)準(zhǔn)單位為86 K)，f=1.75。整個(gè)體系的大小為33.72×33.72×33.72，密度為1 406 kg/m3，盒子尺寸取1.700，體系含有32 000個(gè)粒子，整個(gè)體系劃分成了20×20×20個(gè)盒子。使Ar的粗?；Ｗ釉诿恳粋€(gè)時(shí)間步內(nèi)以不同的概率旋轉(zhuǎn)不同的角度(90°，270°，360°)，導(dǎo)熱系數(shù)隨著SRD時(shí)間步長(zhǎng)的變化情況如圖6所示。當(dāng)粗粒子每一個(gè)時(shí)間步旋轉(zhuǎn)90°，270°，360°的概率分別是1/6，1/6，4/6，且時(shí)間步長(zhǎng)tSRD=1.00，T*=0.71時(shí)，模擬計(jì)算的導(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ù)隨著時(shí)間的變化情況如圖7所示。

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

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

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

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

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

6 結(jié)束語(yǔ)

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