耗散粒子动力学

耗散粒子动力学

（DPD）

内容提要DPD发展历史DPD基本原理DPD边界条件DPD实例计算一二三四DPD发展历史耗散粒子动力学(dissipativeparticledynamics,简称DPD)Hoogerbrugge和Koelman于1992年结合了分子动力学和气体格子法的优点提出了DPD。1995年Espanol和Warren提出DPD中耗散力和随机力中的权函数必须满足涨落—耗散定理。Groot和Warren通过自由能将高分子系统的Flory-Huggins理论与DPD方法相联系。Espanol、Warren和Marsh等人的工作奠定了耗散粒子动力学在统计力学方面的基础。DPD可以用来对复杂流体的动态和静态行为进行模拟。DPD基本原理DPD理论模型DPD方法体系DPD参数取值DPD方法改进一二三四DPD无量纲化五2.1DPD理论模型在DPD系统中,基本的单元是一些离散的被称为“粒子”的动量载体(如图1所示),这些粒子在连续的空间和离散的时间上运动,每个粒子在一定范围内和周围粒子发生相互作用(如图2所示),粒子之间的相互作用力包括保守力、耗散力和随机力三种。2023/9/7图2粒子之间的相互作用范围图1DPD中的粒子

DPD模型中粒子的概念DPD系统中所谓的“粒子”具有粗粒化(coarse-graining)的概念,每颗粒子的运动代表的是大量分子（即所谓的流体粒子)的集体行为。DPD系统中的粒子在连续的空间而非离散的格子上运动,因此遵循Galilean不变性。另一方面,DPD又可以理解为宏观的微分流动控制方程在小尺度上的随机描述,从这个意义上说,DPD方法是连接微观分子动力学方法和宏观流体力学方法的一座桥梁,是一种真正的介观尺度的模拟技术。2023/9/7DPD模型中粒子大小的选取在耗散粒子动力学模型中，流体系统由一系列性质相同的粒子组成。这些粒子并非单个分子，而是由若干个分子组成。组成粒子的分子数目的多少与粒子大小，实际计算区域的几何尺寸，以及计算时间等等密切相关。如果组成粒子的分子数目很少，模型只能模拟较小的区域。极端的情形是粒子由单个分子构成，这时模型实际上就是带软作用力的分子动力学模型。而如果组成粒子的分子数目很多，模型能够充分发挥其优势，模拟较大的区域。因此耗散粒子动力学方法可以被视为一种粗粒化的分子动力学方法。DPD的优势DPD方法其优越性在于能够十分容易地实现对复杂流体系统的模拟。复杂流体中常常包含几种不同的组分,DPD已经成功地广泛用于高分子溶液、胶体以及多相流的模拟。与其它的介观尺度的模拟技术相比,DPD模拟的一个显著的优点就是在模拟一个具有不同组分的复杂流体系统时,并不比模拟单组分流体系统更复杂2023/9/72.2DPD方法体系

2.2.1控制方程作用力和运动方程2023/9/7运动方程：作用力：式中drj和dvj分别是位置和速度矢量保守力Fcij耗散力FDij，随机力Frij保守力、耗散力、随机力耗散力为：式中,ωD与ωR为权重因子,为团间距离r的函数。当r>rC＝1时，ωD与ωR均为零。θij(t)为高斯分布的随机函数，即随机力为：式中，分子团间的保守力为：耗散力与随机力之间的关系

耗散力方向与粒子间相对矢量的方向相反，因此减弱粒子间相互作用。其直接结果是减少系统的动能，降低系统的温度。而随机力引起粒子间的随机振动，增加系统的动能，提高系统的温度。耗散力和随机力的相互作用，在满足一定条件下，能使整个系统温度维持在基本恒定的水平上。

耗散力与随机力表达式Espanol与Warren建议权重因子的形式为：依此关系2023/9/7kBT表示DPD模型中的能量2.2DPD方法体系2.2.2DPD数值积分方法目前常用的修正的Velocity-Verlet算法λ=0.65温度稳定效果较好说明：用粒子当前的位置、速度和力来计算下一个时刻的位置和速度,然后再用新的位置和速度计算新的力,进而修正速度,每运行一步,力就更新一次。2.2DPD方法体系2.2.3粒子间作用及后处理:应力张量S的统计方法采用Irving-Kirkwood公式计算

这里mi是DPD粒子的质量，DPD方法中一般令其为1，N为统计的粒子数，V为体积,uia和uiβ是粒子的本动速度分量,其大小为uia=via-vα(x)。这里vα(x)指的是x处的流场速度，尖括号表示整体平均。

压力由应力张量的迹得到。

2.3DPD系统中的参数

在DPD模拟中模型粒子是一种软粒子,相互之间可以重叠,它不再对应于真实原子或几个原子组成的基团,它可能对应于流体中的微小区域,或高分子的数个链段,也可以是系统在较长时间上的平均效果,即它们是粗粒化的粒子。采用不同的粗粒化方法它有着不同的解释。DPD模拟中选择的势能是一个较软的势能模型,它正是考虑到在介观层次上这些微区之间可以相互重叠这一事实。然而这种势能毕竟是一种理想化了的模型,如何得到它与具体的分子系统间的映射关系就成了DPD模拟中的核心问题。

DPD系统中的参数

DPD模拟方法现在已经成为介观模拟的最重要工具之一。它既消除了分子动力学模拟中对系统描述的过多的模拟细节,使得模拟可以在较大的时间和空间尺度上得以进行,同时也引入了局部的流体力学作用,这种流体力学作用据证实对于特定结构的形成是不可缺少。在DPD模拟中模型粒子是一种软粒子,相互之间可以重叠,它不再对应于真实原子或几个原子组成的基团,它可能对应于流体中的微小区域,或高分子的数个链段,也可以是系统在较长时间上的平均效果,即它们是粗粒化的粒子。

2.3DPD系统中的参数

采用不同的粗粒化方法它有着不同的解释。DPD模拟中选择的势能是一个较软的势能模型,它正是考虑到在介观层次上这些微区之间可以相互重叠这一事实。然而这种势能毕竟是一种理想化了的模型,如何得到它与具体的分子系统间的映射关系就成了DPD模拟中的核心问题。DPD流体Flory-Huggins理论自由能模拟体系的相互作用与Flory-Huggins理论的映射相联系

：

1、最直接的方法是用Flory-Huggins模型关联高分子系统的相平衡或其他热力学性质的实验结果而得到参数；

2、利用分子力学方法得到粗粒化粒子系统的混合能而获得参数。保守力参数

FH理论适用于研究液体-液体、液体-固体的高分子混合系统。依据FH理论,两成分的混合系统,其平均自由能可表示为：χ为A与B间的作用参数。即：

将FH理论与DPD方法比较,可推得参数α与χ间的关系为：

随机力与耗散力参数

耗散力和随机力的值通过涨落-耗散定理相关联,两者中只能任取一个,我们把这个参数取为耗散力。Groot和间发现高斯噪声与均匀噪声对体系的运行几乎没有影响,因此选用简单的均匀噪声,当噪声幅度大于σ=8时,积分运行结果不稳定,当σ=3时,在温度区间kBT=1-10内,弛豫过程快而合理。因此,在选择耗散力参数时,所选的温度和参数不能使噪声比σ=3大,否则模拟结果不可靠。由式（6）可得当kBT=1时，γ=4.5。2.4DPD方法改进2.4.1耗散力的改进为更好地改进

DPD模拟系统的动力学行为，我们将耗散力权函数修正为更通用形式：耗散黏度可以表示为关于s的函数如（下）式所示2.4.1耗散力的改进

耗散力权函数上图显示了权函数在不同s值时的变化趋势。当s=2时权函数是与传统DPD模拟中的权函数表达式是一致的.权函数wD(r)以及它的斜率在r/rc=1时都是连续的，然而当s<1时权函数在r/rc=1c时仍连续，但是它的斜率不是连续的.并且s越小跳动越大.但权函数斜率的不连续对系统的影响不明显.故我们只要求权函数自身的连续，并允许在较小范围内权函数斜率的不连续性.2.4.1耗散力的改进S对粘度的影响2.4.1耗散力的改进三种DPD的动力特性从表中可以看出把耗散力权函数的指数从

2.0变到

0.5就能够明显增加系统

Schmidt数目和耗散黏度。表中也表明了耗散系数γ和截断半径rc对系统动力学性质的决定作用,Sc与γ2和rc8成比例。然而增加Sc对系统温度会造成大的波动。假使系统温度得到有效控制，可以通过增加Sc来改进系统性质，但是实施这种方法有很多困难。很显然最有效的增加

Sc的方法是增加rc，又考虑到对耗散力和随机力的计算要求是它们与rc3成比例。模拟复杂系统时计算耗时又是一个考虑的重要因素。在目前研究中，我们结合修正的权函数和适当地增大截断半径来达到合理的计算耗时。

2.4.2保守力权函数的改进

传统的DPD方法采用一个简单线性形式的保守力权函数ωC(r)=1-r,该式所描述的是一个纯排斥性质的软作用力。采用此种纯排斥性的保守力权函数的DPD模型理论上只能模拟类似气体的特性，因为纯排斥力导致DPD粒子相互排斥离散开来，占据整个计算区域（或容器体积）。虽然传统的DPD模型被广泛应用于多组分多相系统，却不能够模拟单组分多相（气-液-固）系统，因为粒子间只有排斥力，而没有吸引力将粒子结合在一起。因此传统的DPD模型不能模拟带有自由表面的流体流动，中空带气泡的液体，液滴动力学特性等等含气-液甚至气-液-固共存的多相系统。2.4.2保守力权函数的改进因此传统的DPD模型不能模拟带有自由表面的流体流动，中空带气泡的液体，液滴动力学特性等等含气-液甚至气-液-固共存的多相系统。考虑到气-液-固共存及相互转化的物理本质,在保守力势函数中引入吸引力是必需的。引入在DPD粒子非常趋近时表现为排斥形式，而当DPD粒子间距离趋远时表现为吸引形式的保守力势函数能够描述气-液共存的流体系统2.4.2保守力权函数的改进引入在DPD粒子非常趋近时表现为排斥形式，而当DPD粒子间距离趋远时表现为吸引形式的保守力势函数能够描述气-液共存的流体系统。构造这种近距排斥远程吸引的保守力权函数的方式较多。一种可能的方式是把不同作用强度和截距的光滑粒子动力学的光滑函数叠加，从而组合出满足特定要求的保守力权函数。光滑粒子动力学中最常用的是三次样条光滑函数2.4.2保守力权函数的改进采用光滑粒子动力学中最常用的三次样条光滑函数该式对应的系数在一维、二维、三维系统中分别为3/2h、10/7πh2、1/πh32023/9/72.4.2保守力权函数的改进三次样条光滑函数所定义的是非负的单调递减函数，在原点及截距光滑过渡。采用不同的作用强度系数及截距，能够组合构造出如下部分区域为正，部分区域为负的函数

U（r）以作为DPD模型中粒子间相互作用的保守力势函数。式中W1（r）和W2（r）是非正则化的三次样条光滑函数，rc1和rc2，以及A和B分别是对应于W1（r）和W2（r）的截距及作用强度系数。2.4.2保守力权函数的改进不同的rc1和rc2以及A和B组合可以得到不同形式的保守力势函数及相应的权函数。值得注意的是，为了得到近距排斥远程吸引的保守力势函数和权函数，排斥作用的强度应该大于吸引作用的强度，即A＞B,而排斥作用的作用范围应该小于吸引作用的作用范围，即rc1＜rc2。合理选择rc1和rc2以及A和B,就能得到描述特定流体的含近距排斥远程吸引的保守力势函数和权函数。2.4.2保守力权函数的改进所构造的不同的保守力势函数UC(r)=18.75(2W1(r,0.8)-W2(r,1.0))及传统的势函数UC(r)=0.5-(r-0.5r2)2.5DPD单位无量纲化无量纲化：用若干个基本量的乘方之积来表示一个物理量时导出的表达式，叫做此物理量的量纲式，简称为量纲。是用在选定了单位制之后的基本物理量表达出的式子，在模型编制中有量纲的物理量都可以进行无量纲化处理。2023/9/72.5无量纲化为了简化在

DPD模拟时的计算量，我们用对比单位来表示各个参数量纲，用截断半径rc来表示长度单位，用KBT来表示能量单位。DPD体系中与分子动力学不同之处在于,粒子之间的相互作用是通过一些唯象的力来实现的。其中选择粒子间相互作用截断半径与粒子的质量作为一个单位，无量纲化由所有情况通过除以粒子半径与质量来实现，实现无量纲化后利用无量纲的单位来定量系统中的物理量。2.5DPD无量纲化以

为单位来度量时间，同时利用能量均分定理来规一化所有粒子的速度，因此，这样处理后的所有量都是无量纲的。在耗散粒子动力学(DPD)模拟中，物理单位上的长度、速度和时间：r、v和

t，无量纲化定义为下式，这样处理之后，一个

DPD计算模拟的结果可以跟许多种情形的物理体系相互联系起来。周期性边界条件LE边界条件固体壁面边界条件一二三3边界条件3.1周期性边界条件

DPD采用分子动力学中同样的周期性边界条件使用周期性条件的目的就在于让我们用有限的计算资源来模拟无限的空间。这样体系中的粒子数位就会被保持在一个水平常量上。而且在可以使“无限”长的效果使用周期性条件的方向上得以实现。周期性条件被采用在模拟体系之后，这样系统实际上是在周期性边界条件的方向上由模拟单元重复排列构成。从而在模拟计算中只计算基本单元即可得到整个模拟系统的相关数据域性质，剩下的所有其它单元都是基本单元的镜像。基本单元的镜像问题包括两个具体方面的内容：一方面与粒子位移轨迹相关，另一个是最小影像条件。3.1周期性边界条件如左图所，周期性条件被应用在其中的两个方向上，其余八个附近的正方形这样就被当做模拟单元的镜像立方体。在我们模拟体系中的基本模拟单元有一个粒子运动出时(如图2.2所示)，在基本单元的右侧粒子可能会移动出去。此时，会有一个粒子从左侧相邻镜像体系中进入模拟，以此保证粒子数不变。3.1周期性边界条件在右图中所示，粒子

3与其所在的基本单元中，粒子

1、粒子

2与粒子

3的距离都超过了它的截断半径，这时

DPD系统计算时，就会计算与粒子

3的距离在截断半径内的镜像中的粒子，2.3图中计算的是基本单元的下方的镜像中粒子

1与粒子

2对其总的力的作用。然而特别需要注意的是，通常在模拟计算时，模拟体系的边长应远大于截断半径。

3.2DPD中Lees－Edwards边界条件

粒子间作用力的计算

在DPD中如何恰当选取方法计算粒子间相互作用力将在很大程度上影响整个模拟所花费的时间，元胞分割法（cellubdivision）不失为一种高效的方法．如图１所示，计算区域被划分为一系列的元胞区域，每个元胞的边长均大于rC，计算区域内的各ＤＰＤ粒子处于不同的元胞区域内．显然，粒子之间的相互作用仅仅发生在处于同一个元胞内的粒子之间以及分别处于２个紧邻的元胞内的粒子之间；否则，粒子之间的相互距离必定大于rC，从而可以不考虑其相互作用。图1

LE边界条件可模拟周期性的简单剪切流动．图２为周期性系统中简单剪切流动示意图，图中Ｌ为周期性区域在Z方向的长度，为ux方向的流体剪切速度．假定所研究的区域内仅有２个粒子．在粒子运动的过程中，受到处于本区域内其他粒子及周期性镜像区域内粒子的作用．由于剪切率，镜像区域各自角点的x方向速度ｕ(r,t)与ｚ成正比，即ｕ(r,t)＝iDxz,其中i为x方向单位矢量．图2LE边界示意图3.3固体壁面边界条件目前DPD方法大多基于周期性边界条件，用于研究简单区域中复杂流体的流动现象。传统DPD方法把规则分布的DPD粒子冻结在固体边界上代表固体壁面。结合不同反弹形式实施对应的壁面边界条件。流动DPD粒子达到固体壁面后反弹回流动区域的形式包括:1)镜面反射2)弹跳反射3)麦克斯韦反射可以把上述三种壁面反弹形式结合起来形成新的壁面处理方法3.3.1固体壁面边界处理方法在DPD数值模拟中，固体边界处理有三方面的要求:1.精确实施固体壁面无滑移条件，2.尽量降低靠近固体壁面的流动区域中的温度扰动以及其他参数因此产生的波动，3.防止流动DPD粒子非物理性地穿透固体边界。对于复杂区域的流动问题，固体边界处理还应该满足对复杂边界的适应性和计算效率方面的要求。早期DPD方法一般在固体壁面上仅使用一层规则分布的DPD粒子。因为DPD方法使用软球形式的作用势，一层边界DPD粒子不足以阻止某些流动DPD粒子穿透固体壁面，3.3.3边界处理改进方法1具体措施：1)把整个计算区域用规则背景网格覆盖2)根据所需的粒子密度在计算区域中随机布置DPD粒子3)冻结固体障碍区域中DPD粒子4)根据每个边界DPD粒子所处网格及邻近网格种类计算粒子的法线方向及切线方向，法向及切向均指向流动区域内部。5)在临近固液边界的流动区域中设置流动反弹层6)当流动DPD粒子进入流动反弹层后切向速度大小不变，方向相反。和传统的DPD边界处理方法相比，这种新的边界处理方法具有三个优点:1)新的边界处理方法可以采用随机分布的粒子作为边界粒子，能够方便地表征复杂流动区域。2)新的边界处理方法结合无规则分布的边界粒子和流动反弹层，可以确保流动DPD粒子不穿透固体壁面。3)新的边界处理方法中边界粒子密度与流动区域中流动粒子密度一致，流体粒子反弹时流动DPD粒子法向速度根据系统温度给定重新设置，满足给定麦克斯韦