一种气液相变的介观模拟方法

摘要

本发明涉及一种气液相变的介观模拟方法。该方法采用稠密气体状态方程刻画分子间短程排斥效应，采用成对相互作用力模仿分子间长程吸引效应。该方法基于双分布函数，其中密度分布函数用于描述并求解质量‑动量守恒定律，总动能分布函数用于描述并求解能量守恒定律。密度分布函数格子Boltzmann方程可恢复稠密气体状态方程和成对相互作用力，总动能分布函数格子Boltzmann方程可恢复稠密气体压力功、成对相互作用力做功、表面张力做功和粘性热耗散。该方法拥有明确的微观粒子图像和介观动理学理论背景，兼具概念及计算简洁性，并自然地满足热力学一致性，可实现气液相变过程的直接数值模拟，适用性广、可靠性高。

说明书

技术领域

本发明涉及流体力学、热力学和动理学领域，更具体地涉及一种气液相变的介观模拟方法。

背景技术

气液相变是一种基本的热物理现象，广泛存在于自然界和工程应用中。气液相变过程气液两相之间存在相界面，其位置未知且随时间动态演化。在气液相界面内，相态转变伴随着大量潜热的释放/吸收。气液两相通常存在极大的密度差异，其动力粘度、热导率等物性参数亦显著不同。气液相变在宏观上表现极为复杂，如成核、过热、过饱和、蒸发、沸腾、冷凝等。上述特征使得气液相变的数值模拟极具挑战，现有的宏观数值方法必须依赖于复杂的相界面捕捉与追踪算法，且在处理气液相变过程时唯象地采用了大量假设、近似及简化，无法反映气液相变过程底层的基本物理规律，数值方法的普适性和模拟结果的真实性亦难以保证。

气液相变虽然在宏观上表现极为复杂，但其微观机理则十分简单。从物理上讲，气液相变及与之相关的复杂现象均是微观分子间相互作用的自然体现。微观分子间相互作用通常由短程排斥效应和长程吸引效应构成，其中短程排斥效应为分子有限体积大小的直接体现，可采用稠密气体Enskog理论进行描述，而长程吸引效应则可利用平均场理论近似为局部点力。因此，从微观分子间相互作用的物理角度，对气液相变进行数值建模，具有概念清晰、计算简洁的独特优势，可反映相变过程的物理本质。格子Boltzmann方法是一种计算流体力学的介观方法，它起源于格子气自动机，亦可视为Boltzmann方程的一种特殊离散格式。格子Boltzmann方法兼具介观粒子属性和动理学理论背景，可考虑微观粒子间的相互作用，极为适用于气液相变的物理建模和数值模拟。

发明内容

基于格子Boltzmann方法的理论框架，本发明提供了一种气液相变的介观模拟方法。从物理上讲，该介观模拟方法采用稠密气体状态方程刻画分子间短程排斥效应，采用成对相互作用力模仿分子间长程吸引效应。从数值上讲，该介观模拟方法基于密度分布函数描述并求解质量-动量守恒定律，引入总动能分布函数描述并求解能量守恒定律。

本发明提供的气液相变的介观模拟方法包括如下步骤：

S1：选定实际气体状态方程及参数，确定初始温度，确定气相和液相的饱和密度，设定气液两相间的表面张力和相界面宽度；

S2：设定网格空间步长和模拟区域网格数，计算相互作用强度、格子声速、时间步长、定容比热容等模拟参数；

S3：初始化网格点处的密度、速度、总动能、温度、压力等物理量，根据密度场计算成对相互作用力，初始化密度分布函数和总动能分布函数；

S4：执行密度分布函数格子Boltzmann方程的局部碰撞过程，得到碰撞后的密度分布函数；执行总动能分布函数格子Boltzmann方程的局部碰撞过程，得到碰撞后的总动能分布函数；

S5：执行密度分布函数格子Boltzmann方程的线性迁移过程，得到下一时刻的密度分布函数；执行总动能分布函数格子Boltzmann方程的线性迁移过程，得到下一时刻的总动能分布函数；

S6：计算下一时刻的密度，根据密度场计算成对相互作用力，计算下一时刻的速度、总动能、温度、压力等物理量；

S7：根据实际边界条件确定边界网格点处的密度、速度、总动能、温度、压力等物理量，采用格子Boltzmann方法中的边界条件处理格式构造边界网格点处的密度分布函数和总动能分布函数；

S8：重复步骤S4～S7，直至气液相变结束或到达指定时刻。

本发明所述的总动能的热力学定义为：总动能ρe

本发明所述的总动能在介观层面的物理意义为：密度ρ、速度u、内动能ρ∈

本发明所述的内动能满足：内动能ρ∈

本发明所述的内位能的物理本质为：内位能ρ∈

本发明所述的内位能的处理方式为：内位能ρ∈

本发明所述的密度、速度、总动能、温度、压力的计算方式为：密度ρ和速度u根据密度分布函数计算，总动能ρe

本发明所述的演化方程满足：密度分布函数格子Boltzmann方程应可恢复稠密气体状态方程和成对相互作用力。

本发明所述的演化方程满足：总动能分布函数格子Boltzmann方程应可恢复稠密气体压力做功、成对相互作用力做功、表面张力做功和粘性热耗散。

本发明的有益效果：本发明涉及一种气液相变的介观模拟方法。该方法采用稠密气体状态方程刻画分子间短程排斥效应，采用成对相互作用力模仿分子间长程吸引效应。该方法基于双分布函数，其中密度分布函数用于描述并求解质量-动量守恒定律，总动能分布函数用于描述并求解能量守恒定律。密度分布函数格子Boltzmann方程可恢复稠密气体状态方程和成对相互作用力，总动能分布函数格子Boltzmann方程可恢复稠密气体压力功、成对相互作用力做功、表面张力做功和粘性热耗散。该方法拥有明确的微观粒子图像和介观动理学理论背景，兼具概念及计算简洁性，并自然地满足热力学一致性，可实现气液相变过程的直接数值模拟，适用性广、可靠性高。

附图说明

图1是本发明的气液相变的介观模拟方法的计算流程图。

图2是实施例二维空间液滴蒸发示意图。

图3是实施例二维空间液滴蒸发无量纲液滴直径的平方随无量纲时间的演化，以及四个时刻液滴附近的局部密度场和温度场。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明作进一步的详细说明。

更具体地，本发明的气液相变的介观模拟方法，包括如下步骤，如图1所示：

S1：选定实际气体状态方程及参数，确定初始温度，确定气相和液相的饱和密度，设定气液两相间的表面张力和相界面宽度；

S2：设定网格空间步长和模拟区域网格数，计算相互作用强度、格子声速、时间步长、定容比热容等模拟参数；

S3：初始化网格点处的密度、速度、总动能、温度、压力等物理量，根据密度场计算成对相互作用力，初始化密度分布函数和总动能分布函数；

S4：执行密度分布函数格子Boltzmann方程的局部碰撞过程，得到碰撞后的密度分布函数；执行总动能分布函数格子Boltzmann方程的局部碰撞过程，得到碰撞后的总动能分布函数；

S5：执行密度分布函数格子Boltzmann方程的线性迁移过程，得到下一时刻的密度分布函数；执行总动能分布函数格子Boltzmann方程的线性迁移过程，得到下一时刻的总动能分布函数；

S6：计算下一时刻的密度，根据密度场计算成对相互作用力，计算下一时刻的速度、总动能、温度、压力等物理量；

S7：根据实际边界条件确定边界网格点处的密度、速度、总动能、温度、压力等物理量，采用格子Boltzmann方法中的边界条件处理格式构造边界网格点处的密度分布函数和总动能分布函数；

S8：重复步骤S4～S7，直至气液相变结束或到达指定时刻。

根据上述步骤，本发明实施例中的总动能的热力学定义为：总动能ρe

本发明实施例中的总动能在介观层面的物理意义为：密度ρ、速度u、内动能ρ∈

本发明实施例中的内动能满足：内动能ρ∈

本发明实施例中的内位能的物理本质为：内位能ρ∈

本发明实施例中的内位能的处理方式为：内位能ρ∈

本发明实施例中的密度、速度、总动能、温度、压力的计算方式为：密度ρ和速度u根据密度分布函数计算，总动能ρe

本发明实施例中的演化方程满足：密度分布函数格子Boltzmann方程应可恢复稠密气体状态方程和成对相互作用力。

本发明实施例中的演化方程满足：总动能分布函数格子Boltzmann方程应可恢复稠密气体压力做功、成对相互作用力做功、表面张力做功和粘性热耗散。

具体应用实施例

1)本实施例以图2所示的二维空间液滴蒸发为例，模拟计算液滴蒸发过程液滴直径随时间的变化，以及密度场和温度场随时间的演化。

2)采用如下Carnahan-Starling状态方程

其中

3)对于二维情形，采用标准D2Q9离散速度集，9个离散速度为

其中c＝δ

格子声速设置为

其中比例因子K

4)网格点处的密度初始化为

其中x

其中

其中

m

F

n

q

此处，

5)在矩空间执行密度分布函数格子Boltzmann方程的局部碰撞过程

其中m＝M[f

Q

此处，

在矩空间执行总动能分布函数格子Boltzmann方程的局部碰撞过程

其中n＝M[g

Y用于恢复粘性热耗散

碰撞后的总动能分布函数为

6)在速度空间执行密度分布函数格子Boltzmann方程的线性迁移过程，得到下一时刻的密度分布函数

在速度空间执行总动能分布函数格子Boltzmann方程的线性迁移过程，得到下一时刻的总动能分布函数

7)计算下一时刻的密度

根据方程(6)计算成对相互作用力F

根据关系式

8)图2所示的模拟区域四边的边界条件为外流出边界条件、恒压边界条件、恒温边界条件，由此可确定边界网格点处的密度、速度、总动能、温度、压力等物理量，采用格子Boltzmann方法中的边界条件处理格式构造边界网格点处的密度分布函数和总动能分布函数。

9)重复步骤5～8，直至无量纲时间

10)图3给出了无量纲液滴直径的平方(D/D

以上所述是本发明在二维情形下的一种实施方式，应当指出，对于本领域的技术人员来说，在不脱离本发明精神和原则的前提下，还可以做出若干修改、替换和改进等，这些修改、替换和改进也视为本发明的保护范围。