复杂流动与流体行为的拟颗粒模拟

复杂流动与流体行为的拟颗粒模拟1

王利民1,2，葛蔚1，陈飞国1,2，侯超峰1,2，卢健新1,2，张家元1 1．中国科学院过程工程研究所多相复杂系统国家重点实验室，北京（100080）

2. 中国科学院研究生院，北京（100049）

E-mail：wge@

摘 要：复杂流动和流体行为中的多尺度结构对化学工程的理论和应用研究都是巨大的挑战。多尺度的粒子模拟是处理复杂界面、物性以及连续介质模型难以适用的纳微和极端条件下流动的有力手段之一。本文将重点介绍本实验室近年来以拟颗粒模拟为基础的微观粒子模拟研究及其在流态化和纳微流动中的应用，并在此基础上探索离散模拟通用化的可行性。关键词：化学工程；多尺度结构；拟颗粒模拟；多相系统；纳微流动

中图分类号：TQ 028 文献标识码：A

0. 引言

化工系统中普遍存在流动、传递与反应等多种过程的耦合、气液固等多相流动的耦合以及高温高压条件等，使系统具有显著的非平衡性，从而表现出典型的动态多尺度结构、多态性和突变现象、放大效应等[1]。要严格描述这些多尺度结构通常有两种途径：宏观结构的形成可用“自顶向下”(Top Down, TD)方法追溯到连续介质模型能够成立的最小尺度，即进行直接数值模拟(DNS)。传统上基于网格方法的DNS在处理运动颗粒或气泡等复杂边界时存在不少困难，比如规则网格中界面的描述问题和不规则网格中界面的追踪和网格映射问题。而对介观结构的形成，宏观连续介质理论已无法适用，目前还缺乏完善的理论，需要采取“由底而上”(Bottom Up, BU)方法, 通过微观模拟来探索进而检验。在这两方面，离散模拟都能发挥重要的作用。

离散模拟也称粒子方法或无网格方法，它的定义仍在发展中，但总体特征是系统由大量粒子组成、粒子间作用决定系统的行为、粒子间作用基本满足广义力的形式，具有可叠加性及近程性(渐远渐弱)。物质本来就是由大量相互作用的原子和分子组成的，因此离散模拟是微观模拟中最自然的方法，即分子动力学方法

MD(Molecular Dynamics)[2,3]。而宏观上，

离散元方法DEM(Discrete Element Method) [4]也是模拟颗粒物料的直观方法。但更多

的粒子方法可以理解为一种计算模型，使

用的是某种虚拟粒子，如微观的直接模拟

蒙特卡罗方法DSMC(Direct Simulation Monte Carlo)[5,6]，介观的耗散粒子动力学

方法DPD(Dissipative Particle Dynamics)[7]，格子气自动机LGA(Lattice Gas Automata)[8]，格子波尔兹曼方法LBM(Lattice Boltzmann Method)[9]，以及宏

观的光滑粒子动力学方法SPH(Smooth Particle Hydrodynamics)[10~12]和流体粒子

方法FPM(Fluid Particle Model)[13]。本文重

点介绍的拟颗粒模拟(Pseudo-particle modeling, PPM)[14,15]中也构造了一种变径

的微观虚拟粒子，即拟颗粒(PP)。

1. 拟颗粒模拟(Pseudo-Particle Modeling，PPM)

拟颗粒模拟最初是针对气固多相系统

的微观模拟提出的。拟颗粒采取硬粒子间

的相互作用方式与软粒子的自由运动方

式，从而结合了直接模拟蒙特卡洛方法的

效率、并行性和能量守恒优势以及分子动

力学模型的可靠性和微观表现力[16]。如图

1所示，每个拟颗粒具有4个属性：质量

1.本课题得到国家自然科学基金(20336040，20221603和20490201)的资助。

(m )、半径(r )、位置(P )和速度(v )。前两个属性在模拟中是一般恒定不变的。模拟中拟颗粒按相同的时间步长(　t )同步地演化。在各步间，拟颗粒先各自独立运动（可能受到外力或约束），然后在各步结束时，如果两个拟颗粒(如1和2)的距离小于或等于它们的半径的和并且互相靠近，即满

足|P 1- P 2|<(r 1+r 2)和(P 1- P 2)⋅ (V 1-V 2)<0，则它们会象两个光滑刚性球(或二维时的刚性碟)那样发生碰撞，碰撞的情形如图2所示。其中n 12为从拟颗粒1到2的质心连线的单位矢量；t 12为切向单位矢量，n 12⋅ t 12=0。

图1 拟颗粒模拟的示意图

Fig.1 Illustration of Pseudo-Particle Modeling

图2 两个拟颗粒碰撞的分析

Fig.2 Explanation of the collision between

two pseudo particles

拟颗粒模拟最初是针对气固多相系统的微观模拟提出的。拟颗粒采取硬粒子间的相互作用方式与软粒子的自由运动方式，从而结合了直接模拟蒙特卡洛方法的效率、并行性和能量

由动量守恒定律，可得出碰撞前后的速度关系为：

1122

1122m m m m +=+′′v v v v (1)

其中(v 1，v 2)和(v ′1，v ′2)分别是质量为m 1和m 2的两个拟颗粒在碰撞前和碰撞后的速度。

211221||−=

−P P n P P

(2)

1212()

=−v v v (3) ,12121212

()n =⋅n n v v

(4)

将式(2)和式(3)代人式(4)可得

1221,12212

21()()

()

||n −⋅−=

−−P P P P P v P v v (5)

由于拟颗粒之间的碰撞是完全弹性，故有

,12,12n n =−'v v (6) ,12,12

t t ='v v

(7)

由式(6)与式(7)相加可得