董益华　林俊光　张曦　罗海华

摘      要：物性数据是科学研究的基础。根据分子的统计运动规律，以氮气、氢气和氦气为例，采用蒙特卡罗方法预测这三种气体的导热系数和黏性系数。建立了预测模型，并将预测结果进行比较。结果表明：在常温范围，蒙特卡罗方法的预测值和NIST参考值、理论计算值非常接近。说明蒙特卡罗方法在物性参数计算上是可行的。但是在低温和高温范围，蒙特卡罗方法的预测值偏差较大。

关  键  词：蒙特卡罗方法;导热系数;黏性系数;NIST参考值;理论计算值

中图分类号：TQ 015      文献标识码： 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

物性数据是科学研究的基础，在工程设计、过程模拟、科学研究及物质检测等方面都是不可或缺的^[1]。目前，物性数据的研究已经到了模拟仿真的阶段。范围包括基本物性常数、热力学性质、传递性质、微观参数等。

对于导热系数、黏性系数等基础热物性，有Russell 模型^[2]、Maxwell-Eucken 模型^[3]、Y. Agari 模型等以及最小热阻法、逾渗理论法、分形理论法等计算方法^[4]。有的还结合试验数据，研究并编制了多组分气体的热物性参数的计算模型软件。这些方法的准确性各不相同。

从微观上看，物质是由大量分子组成的。分子的大量无规则运动决定了物质的微观世界充满了各种偶然和随机现象。所以，从本质上讲，物理问题都可以用随机问题来处理。蒙特卡罗方法的随机抽样理论恰好应用于物质的微观计算，有着其他方法不可替代的计算手段^[5]。蒙特卡罗方法也已经应用在原子弹工程、气体动力学、物理工程、计算生物学等许多复杂的物理计算系统中^[6-8]。

因此本文从微观角度，计算气体的导热系数和黏性系数，研究蒙特卡罗方法在气体热物性预测中的应用。

1  蒙特卡罗方法预测的基本原理

1.1  气体分子的速率分布

蒙特卡罗方法的理论依据来自于大数定理和中心极限定理，两者都具是渐进性质，需要进行多次抽样才能显示出比较好的结果。

蒙特卡罗方法在物理学中的主要分为两大类，分别是随机性问题（如粒子的运输问题）和确定性问题（如多重积分的计算问题）。用蒙特卡罗方法求解问题时，无论是哪一类问题，均需要将其看作一个随机事件，即在计算机上利用随机变量进行假想试验，当试验次数足够多时，得出事件的概率或算术平均值就可以看作是求解问题的近似解。

气体分子以不同的速率从不同的方向作无规则运动，并且频繁发生碰撞从而改变分子的速率和方向。所以对于每个分子来说，它的速率均是随机的。平衡状态下，大量分子的速率分布是有一定统计规律的。根据麦克斯韦速度分布率

（1）

引入球坐标，可得速率分布函数

（2）

即可得到分子的平均速率

（3）

用蒙特卡罗方法计算积分 的数值解就是分子的平均速率，共可分为三个步骤：

1）将积分式改写成 的形式，依据概率分布 不断的生成随机数x，并依次计算f（x）的值;

2）累加这些计算值，并在[a， b]区间内求平均值的近似值

（4）

3）设定生成N个随机数x，到达停止条件后退出运算。

当N越大时，上述计算的近似值就越接近真实值。当N足够大时，用上述方法计算分子的平均速率非常逼近用公式（5）所得的计算值。

（5）

1.2  确定合适的粒子数N

由图1可知，当粒子数N越大时，随机试验求得的值在真实值上下的浮动值就越小。但当粒子数N增大到10⁸时，虽然在真实值附近的浮动非常小，但计算缓慢，费时。所以粒子数确定为10⁶或10⁷比较合适，与真实值也非常接近，且计算简单、迅速。

综合考虑计算精度和运行时间，本文计算的粒子数最后确定为10⁷。

2  气体导热系数的预测

2.1  导热系数的预测模型

导热系数又称热导率，是导热能力的标志，其定义为：在稳定传热条件下，1 m厚的材料，两侧表面的温差为1 ℃，在1 h内，通过1 m²面积传递的热量。在数值上，它等于在单位温度梯度的作用下物体内的热流密度^[9-10]。

根据傅里叶导热定律，在dt时间内通过dA面沿y方向传输的热量为：

（6）

根据分子内的平均能量和分子数密度，利用分子平均自由程和碰撞频率，可得在dt时间内通过dA面沿x方向传输的热量dQ为

（7）

与傅里叶导热定律类比，导热系数λ为

（8）

其中

（9）

式（9）表明，气体的导热系数λ与分子的平均速率 、平均自由程l和气体的密度ρ、定体比热 有关，因此λ取决于气体的性质和状态。

2.2  导热系数的预测结果

用蒙特卡罗方法分别计算氮气、氢气和氦气的导热系数。这三种气体的导热系数结构参数如表1所示。计算结果见图2—图4。

3  气体黏性系数的预测

3.1  黏性系数的预测模型

流体的黏性是工质热物性研究的重要内容之一，它是化工、制冷、能源、材料等应用领域中不可缺少的基础数据。黏性是流体（液体或气体）的一个重要性质^[11]，是流体抵抗流动的度量。实际流体都具有黏性，都产生摩擦力，而气体的黏度是表征气体内摩擦力的参数。

若粒子在 处x方向的动量是 ，通过泰勒级数展开，可以用y=y_r处的动量来表示它。

（10）

而由上到下的分子流是 ，将它乘以 得到离开上侧的动量流，

（11）

同理可得离开下侧的动流量为，

（12）

表面从下部到上部的净动流量为，

（13）

根据流体力学可知，此一维流动有

（14）

而 就是 ，进而可得到黏性系数 是

（15）

3.2  黏性系数的预测结果

氮气、氢气和氦气的黏性系数结构参数如表2所示。计算结果见图5—图7。

4  结果分析

图2-7的分析中，NIST参考值来自美国国家标准与技术研究院（National Institute of Standards and Technology，NIST），该研究院直属美国商务部，从事物理、生物和工程方面的基础和应用研究，以及测量技术和测试方法方面的研究，提供標准、标准参考数据及有关服务，在国际上享有很高的声誉。

理论计算值均来自文中公式。从图2—图7可知，无论是导热系数还是黏性系数，在常温条件下，NIST参考值、理论计算值和蒙特卡罗方法的预测值均非常接近。误差很小。

但是在在低温（10¹K量级）及高温（大于10³K量级）情况下，计算结果与参考值相差较大。这是因为随着温度的升高，分子自由度被逐步激发。以氢气为例，在低温下，只有3个平动自由度;在常温下，有3个平动自由度和2个转动自由度;而在高温下，除了平动和转动自由度外，还有一个振动自由度。

5  结论

蒙特卡罗方法是一种高效、经济、方便、精确度高、容易实现的随机模拟方法。与其他的数值计算方法相比，蒙特卡罗方法有其自己的优点，如计算出结果所用的时间与待解问题的维数无关;受到问题条件限制的影响小;程序简单，结构清晰易懂，在计算机上容易实现蒙特卡罗方法，便于编制和检验;尤其是对于中子输运等物理问题有着不可替代的作用。

本文以氮气、氢气和氦气为例，计算这三种气体的导热系数和黏性系数。结果表明在常温范围，蒙特卡罗方法的预测值和NIST参考值、理论计算值非常接近。说明蒙特卡罗方法在物性参数计算上是可行的。

但是在低温和高温范围，蒙特卡罗方法的预测值偏差较大。而在当今，不少问题已超出了现有的试验条件，此时用蒙特卡罗方法进行数值计算和计算机模拟虽然有一定优势，但仍需修正，以改进预测精度。

参考文献：

[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]姚凯，郑会保，刘运传，等. 导热系数测试方法概述[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]雷桂媛. 关于蒙特卡罗及拟蒙特卡罗方法的若干研究[D]. 浙江：浙江大学，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]郭開文，代少军，岳建锋. 一类变导热系数下三维温度场解析模型[J]. 工程热物理学报，2017，38（8）：1724-1730.

[10]任玉鸿. 圆柱坐标系下非稳态导热问题改进数值求解方法[J]. 当代化工，2015，44（7）：1634-1637

[11]常勇强，曹子栋，赵振兴，等. 多组分气体热物性参数的计算方法[J]. 动力工程学报，2010， 30 （10）：772-776.

基金项目：国家自然科学基金（项目编号：51576066）。

收稿日期：2019-10-25

作者简介：董益华（1979-），男，浙江奉化人，高级工程师，硕士，2004年毕业于东南大学，研究方向：综合能源系统性能测试与节能优化研发。E-mail：dongyihua109@sina.com。