董益華　林俊光　張曦　羅海華

摘????? 要：物性數(shù)據(jù)是科學(xué)研究的基礎(chǔ)。根據(jù)分子的統(tǒng)計(jì)運(yùn)動(dòng)規(guī)律，以氮?dú)?、氫氣和氦氣為例，采用蒙特卡羅方法預(yù)測(cè)這三種氣體的導(dǎo)熱系數(shù)和黏性系數(shù)。建立了預(yù)測(cè)模型，并將預(yù)測(cè)結(jié)果進(jìn)行比較。結(jié)果表明：在常溫范圍，蒙特卡羅方法的預(yù)測(cè)值和NIST參考值、理論計(jì)算值非常接近。說(shuō)明蒙特卡羅方法在物性參數(shù)計(jì)算上是可行的。但是在低溫和高溫范圍，蒙特卡羅方法的預(yù)測(cè)值偏差較大。

關(guān)? 鍵? 詞：蒙特卡羅方法;導(dǎo)熱系數(shù);黏性系數(shù);NIST參考值;理論計(jì)算值

中圖分類號(hào)：TQ 015????? 文獻(xiàn)標(biāo)識(shí)碼： A? ???文章編號(hào)： 1671-0460（2020）07-1269-04

Prediction of Gas Thermophysical Properties Based on Monte Carlo Method

DONG Yi-hua1，2， LIN Jun-guang1，2， ZHANG Xi1，2，LUO Hai-hua1

（1. Zhejiang Energy Group Co.， Ltd.， Hangzhou Zhejiang 310007， China;

2. Zhejiang Energy Technology Research Institute Co.， Ltd.， Hangzhou Zhejiang 311121， China）

Abstract： Physical data are the basis of scientific research. According to the statistical movement of molecules， the thermal conductivity and viscosity of nitrogen， hydrogen and helium were predicted by Monte Carlo method. The prediction model was established， and the results were compared with other sources. The results showed that the predicted value of Monte Carlo method was very close to NIST reference value and theoretical calculation value in the range of normal temperature. So the Monte Carlo method is feasible in the calculation of physical parameters. But， at relative low temperature or high temperature， the predicted value of Monte Carlo method has large deviation.

Key words： Monte Carlo method; Thermal conductivity; Viscosity coefficient; NIST reference value; Theoretical calculation value

物性數(shù)據(jù)是科學(xué)研究的基礎(chǔ)，在工程設(shè)計(jì)、過(guò)程模擬、科學(xué)研究及物質(zhì)檢測(cè)等方面都是不可或缺的^[1]。目前，物性數(shù)據(jù)的研究已經(jīng)到了模擬仿真的階段。范圍包括基本物性常數(shù)、熱力學(xué)性質(zhì)、傳遞性質(zhì)、微觀參數(shù)等。

對(duì)于導(dǎo)熱系數(shù)、黏性系數(shù)等基礎(chǔ)熱物性，有Russell 模型^[2]、Maxwell-Eucken 模型^[3]、Y. Agari 模型等以及最小熱阻法、逾滲理論法、分形理論法等計(jì)算方法^[4]。有的還結(jié)合試驗(yàn)數(shù)據(jù)，研究并編制了多組分氣體的熱物性參數(shù)的計(jì)算模型軟件。這些方法的準(zhǔn)確性各不相同。

從微觀上看，物質(zhì)是由大量分子組成的。分子的大量無(wú)規(guī)則運(yùn)動(dòng)決定了物質(zhì)的微觀世界充滿了各種偶然和隨機(jī)現(xiàn)象。所以，從本質(zhì)上講，物理問(wèn)題都可以用隨機(jī)問(wèn)題來(lái)處理。蒙特卡羅方法的隨機(jī)抽樣理論恰好應(yīng)用于物質(zhì)的微觀計(jì)算，有著其他方法不可替代的計(jì)算手段^[5]。蒙特卡羅方法也已經(jīng)應(yīng)用在原子彈工程、氣體動(dòng)力學(xué)、物理工程、計(jì)算生物學(xué)等許多復(fù)雜的物理計(jì)算系統(tǒng)中^[6-8]。

因此本文從微觀角度，計(jì)算氣體的導(dǎo)熱系數(shù)和黏性系數(shù)，研究蒙特卡羅方法在氣體熱物性預(yù)測(cè)中的應(yīng)用。

1 ?蒙特卡羅方法預(yù)測(cè)的基本原理

1.1? 氣體分子的速率分布

蒙特卡羅方法的理論依據(jù)來(lái)自于大數(shù)定理和中心極限定理，兩者都具是漸進(jìn)性質(zhì)，需要進(jìn)行多次抽樣才能顯示出比較好的結(jié)果。

蒙特卡羅方法在物理學(xué)中的主要分為兩大類，分別是隨機(jī)性問(wèn)題（如粒子的運(yùn)輸問(wèn)題）和確定性問(wèn)題（如多重積分的計(jì)算問(wèn)題）。用蒙特卡羅方法求解問(wèn)題時(shí)，無(wú)論是哪一類問(wèn)題，均需要將其看作一個(gè)隨機(jī)事件，即在計(jì)算機(jī)上利用隨機(jī)變量進(jìn)行假想試驗(yàn)，當(dāng)試驗(yàn)次數(shù)足夠多時(shí)，得出事件的概率或算術(shù)平均值就可以看作是求解問(wèn)題的近似解。

氣體分子以不同的速率從不同的方向作無(wú)規(guī)則運(yùn)動(dòng)，并且頻繁發(fā)生碰撞從而改變分子的速率和方向。所以對(duì)于每個(gè)分子來(lái)說(shuō)，它的速率均是隨機(jī)的。平衡狀態(tài)下，大量分子的速率分布是有一定統(tǒng)計(jì)規(guī)律的。根據(jù)麥克斯韋速度分布率

（1）

引入球坐標(biāo)，可得速率分布函數(shù)

（2）

即可得到分子的平均速率

（3）

用蒙特卡羅方法計(jì)算積分 的數(shù)值解就是分子的平均速率，共可分為三個(gè)步驟：

1）將積分式改寫(xiě)成 的形式，依據(jù)概率分布 不斷的生成隨機(jī)數(shù)x，并依次計(jì)算f（x）的值;

2）累加這些計(jì)算值，并在[a， b]區(qū)間內(nèi)求平均值的近似值

（4）

3）設(shè)定生成N個(gè)隨機(jī)數(shù)x，到達(dá)停止條件后退出運(yùn)算。

當(dāng)N越大時(shí)，上述計(jì)算的近似值就越接近真實(shí)值。當(dāng)N足夠大時(shí)，用上述方法計(jì)算分子的平均速率非常逼近用公式（5）所得的計(jì)算值。

（5）

1.2 ?確定合適的粒子數(shù)N

由圖1可知，當(dāng)粒子數(shù)N越大時(shí)，隨機(jī)試驗(yàn)求得的值在真實(shí)值上下的浮動(dòng)值就越小。但當(dāng)粒子數(shù)N增大到10⁸時(shí)，雖然在真實(shí)值附近的浮動(dòng)非常小，但計(jì)算緩慢，費(fèi)時(shí)。所以粒子數(shù)確定為10⁶或10⁷比較合適，與真實(shí)值也非常接近，且計(jì)算簡(jiǎn)單、迅速。

綜合考慮計(jì)算精度和運(yùn)行時(shí)間，本文計(jì)算的粒子數(shù)最后確定為10⁷。

2 ?氣體導(dǎo)熱系數(shù)的預(yù)測(cè)

2.1 ?導(dǎo)熱系數(shù)的預(yù)測(cè)模型

導(dǎo)熱系數(shù)又稱熱導(dǎo)率，是導(dǎo)熱能力的標(biāo)志，其定義為：在穩(wěn)定傳熱條件下，1 m厚的材料，兩側(cè)表面的溫差為1 ℃，在1 h內(nèi)，通過(guò)1 m²面積傳遞的熱量。在數(shù)值上，它等于在單位溫度梯度的作用下物體內(nèi)的熱流密度^[9-10]。

根據(jù)傅里葉導(dǎo)熱定律，在dt時(shí)間內(nèi)通過(guò)dA面沿y方向傳輸?shù)臒崃繛椋?/p>

（6）

根據(jù)分子內(nèi)的平均能量和分子數(shù)密度，利用分子平均自由程和碰撞頻率，可得在dt時(shí)間內(nèi)通過(guò)dA面沿x方向傳輸?shù)臒崃縟Q為

（7）

與傅里葉導(dǎo)熱定律類比，導(dǎo)熱系數(shù)λ為

（8）

其中

（9）

式（9）表明，氣體的導(dǎo)熱系數(shù)λ與分子的平均速率 、平均自由程l和氣體的密度ρ、定體比熱 有關(guān)，因此λ取決于氣體的性質(zhì)和狀態(tài)。

2.2 ?導(dǎo)熱系數(shù)的預(yù)測(cè)結(jié)果

用蒙特卡羅方法分別計(jì)算氮?dú)?、氫氣和氦氣的?dǎo)熱系數(shù)。這三種氣體的導(dǎo)熱系數(shù)結(jié)構(gòu)參數(shù)如表1所示。計(jì)算結(jié)果見(jiàn)圖2—圖4。

3? 氣體黏性系數(shù)的預(yù)測(cè)

3.1? 黏性系數(shù)的預(yù)測(cè)模型

流體的黏性是工質(zhì)熱物性研究的重要內(nèi)容之一，它是化工、制冷、能源、材料等應(yīng)用領(lǐng)域中不可缺少的基礎(chǔ)數(shù)據(jù)。黏性是流體（液體或氣體）的一個(gè)重要性質(zhì)^[11]，是流體抵抗流動(dòng)的度量。實(shí)際流體都具有黏性，都產(chǎn)生摩擦力，而氣體的黏度是表征氣體內(nèi)摩擦力的參數(shù)。

若粒子在 處x方向的動(dòng)量是 ，通過(guò)泰勒級(jí)數(shù)展開(kāi)，可以用y=y_r處的動(dòng)量來(lái)表示它。

（10）

而由上到下的分子流是 ，將它乘以 得到離開(kāi)上側(cè)的動(dòng)量流，

（11）

同理可得離開(kāi)下側(cè)的動(dòng)流量為，

（12）

表面從下部到上部的凈動(dòng)流量為，

（13）

根據(jù)流體力學(xué)可知，此一維流動(dòng)有

（14）

而 就是 ，進(jìn)而可得到黏性系數(shù) 是

（15）

3.2 ?黏性系數(shù)的預(yù)測(cè)結(jié)果

氮?dú)?、氫氣和氦氣的黏性系?shù)結(jié)構(gòu)參數(shù)如表2所示。計(jì)算結(jié)果見(jiàn)圖5—圖7。

4? 結(jié)果分析

圖2-7的分析中，NIST參考值來(lái)自美國(guó)國(guó)家標(biāo)準(zhǔn)與技術(shù)研究院（National Institute of Standards and Technology，NIST），該研究院直屬美國(guó)商務(wù)部，從事物理、生物和工程方面的基礎(chǔ)和應(yīng)用研究，以及測(cè)量技術(shù)和測(cè)試方法方面的研究，提供標(biāo)準(zhǔn)、標(biāo)準(zhǔn)參考數(shù)據(jù)及有關(guān)服務(wù)，在國(guó)際上享有很高的聲譽(yù)。

理論計(jì)算值均來(lái)自文中公式。從圖2—圖7可知，無(wú)論是導(dǎo)熱系數(shù)還是黏性系數(shù)，在常溫條件下，NIST參考值、理論計(jì)算值和蒙特卡羅方法的預(yù)測(cè)值均非常接近。誤差很小。

但是在在低溫（10¹K量級(jí)）及高溫（大于10³K量級(jí)）情況下，計(jì)算結(jié)果與參考值相差較大。這是因?yàn)殡S著溫度的升高，分子自由度被逐步激發(fā)。以氫氣為例，在低溫下，只有3個(gè)平動(dòng)自由度;在常溫下，有3個(gè)平動(dòng)自由度和2個(gè)轉(zhuǎn)動(dòng)自由度;而在高溫下，除了平動(dòng)和轉(zhuǎn)動(dòng)自由度外，還有一個(gè)振動(dòng)自由度。

5 ?結(jié)論

蒙特卡羅方法是一種高效、經(jīng)濟(jì)、方便、精確度高、容易實(shí)現(xiàn)的隨機(jī)模擬方法。與其他的數(shù)值計(jì)算方法相比，蒙特卡羅方法有其自己的優(yōu)點(diǎn)，如計(jì)算出結(jié)果所用的時(shí)間與待解問(wèn)題的維數(shù)無(wú)關(guān);受到問(wèn)題條件限制的影響小;程序簡(jiǎn)單，結(jié)構(gòu)清晰易懂，在計(jì)算機(jī)上容易實(shí)現(xiàn)蒙特卡羅方法，便于編制和檢驗(yàn);尤其是對(duì)于中子輸運(yùn)等物理問(wèn)題有著不可替代的作用。

本文以氮?dú)?、氫氣和氦氣為例，?jì)算這三種氣體的導(dǎo)熱系數(shù)和黏性系數(shù)。結(jié)果表明在常溫范圍，蒙特卡羅方法的預(yù)測(cè)值和NIST參考值、理論計(jì)算值非常接近。說(shuō)明蒙特卡羅方法在物性參數(shù)計(jì)算上是可行的。

但是在低溫和高溫范圍，蒙特卡羅方法的預(yù)測(cè)值偏差較大。而在當(dāng)今，不少問(wèn)題已超出了現(xiàn)有的試驗(yàn)條件，此時(shí)用蒙特卡羅方法進(jìn)行數(shù)值計(jì)算和計(jì)算機(jī)模擬雖然有一定優(yōu)勢(shì)，但仍需修正，以改進(jìn)預(yù)測(cè)精度。

參考文獻(xiàn)：

[1]GINER-SANZ J J， ORTEGA E M， P?REZ-HERRANZ V.? Application of a Montecarlo based quantitative Kramers-Kronig test for linearity assessment of EIS measurements[J].ElectrochimicaActa， 2016， 209（2）： 254-268.

[2]姚凱，鄭會(huì)保，劉運(yùn)傳，等. 導(dǎo)熱系數(shù)測(cè)試方法概述[J].理化檢測(cè)-物理分冊(cè)，2018，54（10）：741-747.

[3]HE Y. Rapid thermal conductivity measurement with a hot disk sensor， Part 1. Theoretical considerations [J].Themochimica A cta， 2005， 436 （1）：122-129.

[4]CARSON J K， LOVATT S J. Experimental measurements of the effective thermal conductivity of a peudo-porous food analogue over a range of porosities and mean pore sizes [J].Journal of Food Engineering， 2004， 63 （1）：87-95.

[5]雷桂媛. 關(guān)于蒙特卡羅及擬蒙特卡羅方法的若干研究[D]. 浙江：浙江大學(xué)，2003：55-64.

[6]GROPE B O H，? ZACHERLE T， NAKAYAMA M， et al. Oxygen ion conductivity of doped ceria： A Kinetic Monte Carlo study[J]. Solid State Ionics， 2012， 225（10）： 476-483.

[7]HUANG B F， FAN F X， LI Y S， et al. Numerical prediction of ultrasonic attenuation in concentrated emulsions and suspensions using Monte Carlo method[J]. Ultrasonics， 2019， 94（4）： 218-226.

[8]Sahar Abdolbaghi， Ali Barati-Harooni， Adel Najafi-Marghmaleki. Improving the prediction ability of reference correlation for viscosity of carbon dioxide[J]. Journal of CO2 Utilization， 2019， 31（5）： 106-114.

[9]郭開(kāi)文，代少軍，岳建鋒. 一類變導(dǎo)熱系數(shù)下三維溫度場(chǎng)解析模型[J]. 工程熱物理學(xué)報(bào)，2017，38（8）：1724-1730.

[10]任玉鴻. 圓柱坐標(biāo)系下非穩(wěn)態(tài)導(dǎo)熱問(wèn)題改進(jìn)數(shù)值求解方法[J]. 當(dāng)代化工，2015，44（7）：1634-1637

[11]常勇強(qiáng)，曹子棟，趙振興，等. 多組分氣體熱物性參數(shù)的計(jì)算方法[J]. 動(dòng)力工程學(xué)報(bào)，2010， 30 （10）：772-776.

基金項(xiàng)目：國(guó)家自然科學(xué)基金（項(xiàng)目編號(hào)：51576066）。

收稿日期：2019-10-25

作者簡(jiǎn)介：董益華（1979-），男，浙江奉化人，高級(jí)工程師，碩士，2004年畢業(yè)于東南大學(xué)，研究方向：綜合能源系統(tǒng)性能測(cè)試與節(jié)能優(yōu)化研發(fā)。E-mail：dongyihua109@sina.com。