基于MPCD方法计算流体的导热系数

2021-04-07王瑞金朱泽飞

杭州电子科技大学学报(自然科学版) 2021年2期

李龙，王瑞金，朱泽飞

(杭州电子科技大学机械工程学院，浙江杭州 310018)

0 引言

多粒子碰撞动力学是一种粗粒化分子模型的方法[1]，结合了分子动力学和格子-玻尔兹曼方法(Lattice-Boltzmann Method，LBM)的优点，能够解决流体的热涨落问题，计算量小，可以模拟复杂边界条件的体系[2]。很多学者使用MPCD方法进行研究，如陈文多等[3]采用MPCD方法研究了环形高分子的流体动力学并计算其回转半径，Akhter等[4]采用MPCD方法对微压缩低雷诺数条件下的颗粒流体进行数值模拟，Wysocki等[5]采用MPCD方法研究胶体系统在微米尺度狭缝中的水动力学不稳定性，Eggersdorfer等[6]采用MPCD方法进行流体中大颗粒聚集结构的研究。本文采用MPCD方法分析了影响流体导热系数计算的因素，并计算AR体系和水分子体系的导热系数。

1 MPCD方法

流体粗粒子采用MPCD进行数值模拟，一般分为流动步和碰撞步。在流动步中，粒子的运动速度不变且可得到粗粒子的新位移。在碰撞步中，同一个盒子中的粗粒子速度矢量绕着随机轴旋转一定的角度，得到粗粒子的新速度。流动过程中粒子属于平动，满足：

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

(1)

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

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

(2)

式中，ri(t)是粒子的位移矢量，vi是粒子速度矢量，Δt是时间步长，α是粒子速度矢量绕任意随机轴旋转的角度，Ω(α)是随机旋转矩阵。vcm(t)是每个盒子的质心速度，质心速度的计算公式如下：

(3)

式中，N是在一个盒子中的粒子数，M是粗粒子的质量。

2 导热系数的计算

(4)

(5)

(6)

3 模拟体系的大小及盒子划分

本文选择分子动力学中经常使用的LJ约化单位，故本文中的参数基本为无量纲参数。

选择模拟体系的尺寸为33.72×33.72×33.72，长度D=d/σ，d是标准单位下的长度，σ是LJ势函数中的单位距离参数，整个体系划分为20×20×20个盒子如图1所示，流体粗粒子体系如图2所示。

图1 体系盒子划分

图2 模拟体系

4 影响导热系数计算结果的因素

4.1 流体粗粒子质量

选取体系温度T*=0.71，划分盒子的尺寸为1.700，晶格常数f=1.55，整个模拟体系含有62 500个粗粒化粒子，每个盒子内含有的平均粒子数为9.112，粒子质量M分别为0.40，0.60，0.80，1.00，1.20。为了得出不同SRD时间步长下的流体导热系数，tSRD分别取0.25，0.30，0.35，0.40。导热系数随着粗粒子质量的变化情况如图3所示。从图3可以看出，导热系数随粗粒子质量的增大而减小。温度一定的情况下，质量越小，粒子的速度越大，所以粒子间碰撞几率增大。从图3还可以看出，导热系数随SRD时间步长的增大而增大。SRD时间步长越大，粒子移动的距离越远，不同盒子内粒子的动量和能量交换就更加频繁，间接增强了流体的传热。

图3 导热系数随着粗粒子质量的变化情况

4.2 盒子尺寸

选取参数粒子质量M=0.40，tSRD=0.35，温度T*=0.71，晶格常数f=1.55。计算体系在不同盒子尺寸和晶格常数下的导热系数。不同尺寸盒子内的平均粒子数不同，影响导热系数的计算结果。本文选取如表1所示的体系参数进行模拟，粒子密度ρ是指每个盒子内的平均粒子数，f采用的是LJ约化单位，盒子尺寸和f均是无量纲参数。标准单位下的晶格常数L与f成反比，即f越大，每个盒子内的平均粒子数就越多。在f不变的情况下，ρ随着盒子尺寸的增大而增大。

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

导热系数随着盒子尺寸的变化情况如图4所示。从图4可以看出，在f相同的情况下，导热系数随着盒子尺寸的增大而减小，虽然盒子尺寸增大使得每个盒子的粒子数增加，但粒子移出盒子的概率减少了，导致粒子间的自相关性增大，粒子不能和更远的盒子中的粒子发生碰撞。

4.3 体系温度

选取参数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，另外，晶格常数f分别为1.55，1.75，1.95。导热系数随着体系温度的变化情况如图5所示。从图5可以看出，导热系数随着温度的增大而增大，在体系粒子密度相同的情况下，温度越高，粗粒子的动能越大，粒子的碰撞几率就高，导热系数越大。但是，一般情况下，大部分液体的导热系数随着温度的增大而减小，因为本文是在温度不同且粒子密度相同的情况下进行模拟实验的，并没有考虑温度对分子间距离的影响，即没有考虑温度对体系粒子密度的影响。大部分液体的导热系数是随着体系密度的减小而减小。温度越大，体系的粒子密度是越小的，晶格常数与粒子密度成正比。从图5还可以看出，T*=0.70，f=1.95时的导热系数大于T*=1.00，f=1.55时的导热系数，虽然温度增加了，但是粒子密度却降低了，导热系数也会下降。

4.4 粒子旋转角度

选取参数T*=0.71，M=0.40，盒子尺寸为1.780，f=1.95，选择旋转固定角度分别为为90°，270°，360°。为了准确计算体系的导热系数，使每个盒子内的粗粒子在同一时间步内旋转不同的角度。

粒子以不同概率旋转不同角度下的体系导热系数计算结果如表2所示。旋转360°意味着在一个SRD时间步内速度矢量是不变的，这和增大SRD时间步长所起到的效果相同。从表2可以看出，旋转360°的概率占比越大，导热系数越大，但概率占比不能为1，否则体系会细粒化，计算结果毫无意义。

图4 导热系数随着盒子尺寸的变化情况

图5 导热系数随着体系温度的变化情况

表2 粒子以不同概率旋转不同角度下的体系导热系数计算结果

5 Ar原子和水分子体系的导热系数计算

在计算体系的导热系数时，温度通常是已定的。本文通过计算Ar原子和水分子体系的导热系数来验证MPCD方法是否可以用于流体导热系数的计算。

选取参数M=1.00，T*=0.71(转化为标准单位为86 K)，f=1.75。整个体系的大小为33.72×33.72×33.72，密度为1 406 kg/m3，盒子尺寸取1.700，体系含有32 000个粒子，整个体系划分成了20×20×20个盒子。使Ar的粗粒化粒子在每一个时间步内以不同的概率旋转不同的角度(90°，270°，360°)，导热系数随着SRD时间步长的变化情况如图6所示。当粗粒子每一个时间步旋转90°，270°，360°的概率分别是1/6，1/6，4/6，且时间步长tSRD=1.00，T*=0.71时，模拟计算的导热系数为0.136 5 W·(m·K)-1，与文献[8]得到的0.132 0 W·(m·K)-1较接近，误差为3.4%。Ar原子体系导热系数随着时间的变化情况如图7所示。

图6 导热系数随着时间步长的变化情况

图7 Ar原子体系导热系数随时间的变化情况

另外，本文还计算了T*=1.00时的导热系数，Ar体系中粒子数变为23 328，因为温度变大，粒子间的实际晶格常数变大，体系中所含的粒子数就减少了。选取盒子尺寸为1.983，通过增大盒子尺寸，可以保证在体系温度不同时的盒子所含的平均粒子数一样。本次模拟的计算结果为0.101 6 W·(m·K)-1。

同理，在模拟水分子体系时，一个水分子粗粒化成一个粗粒子，模拟体系尺寸为33.72×33.72×33.72，模拟参数选取M=0.45(分子量为18)，T*=2.44(T=298 K)。导热系数计算结果为0.599 0 W·(m·K)-1，与实验值0.608 4 W·(m·K)-1相比[9]，误差是1.5%。

本文进行模拟计算时，使用的是8核的Intel-i7服务器，仅用3 min左右即可完成Ar原子或水分子体系的运动模拟和导热系数的计算。相同条件下，采用MD方法，大约需要约6 h[10]，计算量大大减少。