董益華　林俊光　張曦　羅海華

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

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

中圖分類號：TQ 015????? 文獻(xiàn)標(biāo)識碼： A? ???文章編號： 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è)計、過程模擬、科學(xué)研究及物質(zhì)檢測等方面都是不可或缺的^[1]。目前，物性數(shù)據(jù)的研究已經(jīng)到了模擬仿真的階段。范圍包括基本物性常數(shù)、熱力學(xué)性質(zhì)、傳遞性質(zhì)、微觀參數(shù)等。

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

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

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

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

1.1? 氣體分子的速率分布

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

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

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

（1）

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

（2）

即可得到分子的平均速率

（3）

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

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

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

（4）

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

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

（5）

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

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

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

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

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

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

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

（6）

根據(jù)分子內(nèi)的平均能量和分子數(shù)密度，利用分子平均自由程和碰撞頻率，可得在dt時間內(nèi)通過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ù)測結(jié)果

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

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

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

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

若粒子在 處x方向的動量是 ，通過泰勒級數(shù)展開，可以用y=y_r處的動量來表示它。

（10）

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

（11）

同理可得離開下側(cè)的動流量為，

（12）

表面從下部到上部的凈動流量為，

（13）

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

（14）

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

（15）

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

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

4? 結(jié)果分析

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

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

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

5 ?結(jié)論

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

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

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

參考文獻(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]姚凱，鄭會保，劉運(yùn)傳，等. 導(dǎo)熱系數(shù)測試方法概述[J].理化檢測-物理分冊，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]郭開文，代少軍，岳建鋒. 一類變導(dǎo)熱系數(shù)下三維溫度場解析模型[J]. 工程熱物理學(xué)報，2017，38（8）：1724-1730.

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

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

基金項目：國家自然科學(xué)基金（項目編號：51576066）。

收稿日期：2019-10-25

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