纳米流体液滴的耗散粒子动力学方法模拟

Copyright © 2014
版权所有 中国力学学会
地址: 北京市北四环西路15号 邮政编码:100190 Address: No.15 Beisihuanxi Road, Beijing 100190
第八届全国流体力学学术会议 2014年9月18～21日 甘肃兰州
文章编号： CSTAM2014-B01-0333
标题：纳米流体液滴的耗散粒子动力学方法模
拟
作者：沈世元，周哲玮
单位：上海市应用数学与力学研究所
上海大学
第八届全国流体力学学术会议 2014年9月18-21日 甘肃 兰州
CSTAM2014-A26-BS10029
纳米流体液滴的耗散粒子动力学方法模拟
沈世元1,2，周哲玮1,2
（1上海市应用数学与力学研究所，上海闸北区 200070）
（2上海大学，上海闸北区 200070）
摘要 纳米流体是指把直径范围从10nm —100nm 的金属或非金属纳米颗粒分散到水、醇、油等传统物质中形成的新型流体。

纳米颗粒的尺寸和浓度会对纳米流体的表面张力、润湿性和导热性等产生很大的影响，是近年来材料、物理、化学、传热学等众领域的研究热点。

DPD （耗散粒子动力学， Dissipative Particle Dynamics ）是研究介观尺度下粒子运动的有力工具，其算法中的参数与物理系统的关系是研究热点之一。

本文利用DPD 方法模拟介观尺度的纳米液滴，根据纳米液滴的接触角确定DPD 方法中的参数，研究了固壁、液体和纳米颗粒之间的相互作用系数。

由于各种参数可以根据实测的数据来确定，此方法适用于研究实际工程中的问题。

关键词 纳米流体；耗散粒子动力学；接触角；相互作用系数
1.引言
1995年，美国Argonne 国家实验室的Choi 3等人提出了一个崭新的概念—— 纳米流体。

随着纳米技术日益深入人心，相关研究逐渐成为一个热点，并在许多工业领域中得到拓展，比如含有表面活性剂的纳米流体可用来增加石油开采量，改良油污后的土壤；由于其易于浸入固体表面的特性，还常被用于对材料进行优化和改良。

此外，纳米流体的有效导热系数高于相应纯流体，这使其传热性明显增强，因而多用于芯片散热的液冷系统中。

通常，悬浮在流体中的纳米颗粒会受到诸如流动阻力、布朗运动、粒子间扩散及重力等内外因素的影响，其运动规律极其复杂。

1
前人对于纳米流体性质的各个方面做了广泛的研究。

1993年日本东北大学的Masuda 等人2在水中添加平均粒子直径为13nm 的32O Al 和27nm 的2O T i 粒子，制备了不同体积浓度的纳米颗粒胶体并测量了胶体的导热系数;1995年，美国Argonne 国家实验室的Choi 等人以一定的方式和比例在液体中添加纳米级金属或金属氧化物粒子，并称之为纳米流体。

3
其中的氧化物粒子包括O C u 、2O S i 、32O Al 等，另外还有一些金属粒子和碳化物等4。

在1995年Choi 3之后，国内外的学者纷纷对纳米流体展开了深入的研究。

在国
内方面，范庆梅等人5
对纳米流体的热导率和粘度进行了计算。

李云翔等人6
为纳米流体的研究进展做了一个总结，包括：纳米流体稳定性的研究、纳米流体物性的研究、纳米流体传热特性的研究，其中既包括实验方面的研究进展也对纳米流体物性以及传热特性的理论研究进行了系统的总结。

国外方面，B. Davidovitch 等人7研究了考虑了热扰动的粘性液滴在基板上的过程。

杰出的研究人员如Kim 8-11、Vassallo 12、Truong 13等人也在纳米流体的性质研究方面做出了杰出的贡献，这里不再一一赘述。

分子动力学模拟（Molecular Dynamics, MD ）是一种在微观尺度下模拟原子和简单分子运动的方法。

但由于超出这个尺度的时候，MD 方法不适用，我们需要采用另外的方法。

对于介于宏观和微观之间的尺度上的流体动力学行为的研究是目前学术界的热点问题。

这个中间的尺度通常被称为介观尺度，通常指的是10-1000nm 和1ns-10ms 的尺度。

本文所采用的介观尺度下模拟流体的动力学行为的方法是耗散粒子动力学方法（Dissipative Particle Dynamics ，简称DPD ）14。

此方法通过对模拟区域内的粒子进行粗粒化，以减少计算代价，在更短的计算时间内计算更大的
区域。

DPD 方法中的相互作用系数（DPD interaction parameter ）是由内在的原子之间的相互作用决定的，在DPD 方法中表现为DPD 粒子之间的相互作用。

Amitesh Maiti 等人15提出，DPD 的作用系数和DPD 粒子的尺寸有关。

接触角是衡量材料表面性能的一个重要参数，通过测量接触角可获得不同物质在界面之间相互作用的很多信息。

国内外的学者们对接触角的测量以及性质也已经有了一些成熟的研究。

Decker E 等人16对宏观和微观的接触角提出了合适的定义，并且给出了部分测量接触角的方法。

Marmur A 等人17总结了测量接触角方法的一些缺陷，并指出了如何规避这些缺陷。

Saeid Vafaei 等人18更是研究出了与体积无关的更加可以准确反映物性的渐近接触角。

宁乔等人19利用照相机捕捉到的液滴灰度图像，利用边界提取等方法得到液滴的轮廓，采用轮廓拟合的方法进而求得液滴的接触角。

徐志钮等人20认为，如何准确测量接触角对相关学科如材料、医药、半导体、生命科学、油墨工业、电力系统等领域都具有重要意义，并指出，当接触角较小时需要使用基于圆的拟合算法；否则使用基于椭圆的算法。

对接触角的研究文献很多，在这里不再一一赘述。

在本文中，我们主要用到的测量接触角的方法是测量表观接触角，即通过提取
液滴表面粒子进行球拟合，得到液滴在固体表面上的接触角。

具体做法是：对于一个由粒子组成的液滴，首先根据粒子所处位置的粒子密度来判断粒子是否处于界面上，然后将所有处于界面上的粒子提取出来，找到液滴界面所在的球面，通过拟合得到准确的接触角。

除了测量表观接触角之外，以后还有可能会考虑采用测量局部接触角的方法。

局部接触角一般用来描述液滴的局部状态，表征液滴的动态接触行为。

获得局部接触角常用的方法是切线法：首先根据液滴局部区域粒子密度来判断粒子是否处于界面上，然后提取该处界面上的粒子坐标，一般可采用三次多项式来拟合界面轮廓，求导后即可获得该时刻液滴的局部接触角。

2.耗散粒子动力学方法 2.1 基本理论
在1998年，PB Warren 21说明了DPD 粒子的对象，可以指代分子、原子、高分子蛋白质聚团、带电离子团等等，指出了DPD 方法的广泛的适用性。

DPD 算法中粒子采用和MD 算法一样的Lagrange 坐标描述，体系中的粒子运动满足牛顿运动方程：
i i i i f dt v d v dt r d ==, ∑≠++=i
j R
ij D ij C ij i F F F f )( (1)
其运动方程中i r 、i v
为i 粒子的位移和
速度，i f
为i 粒子受到的3种合力，包括
颗粒间的有势作用力c
ij F ，还包括类似于布
朗动力学中的随机力R F ij 和耗散力D
F ij 。

这些有势力都属于软势力：
ij ij C ij C ij e r a F )(ω=
ij ij ij ij D D ij e v e r F ))((⋅-=γω
ij ij ij R R ij e r F
θσω)(= (2)
其中，ij r 即i 、j 粒子间的距离，ij e
为和ij r
方向相同的单位矢量，ij a 为保守力系数，反映粒子i 和j 之间的最大排斥力，γ和σ分别为耗散力和随机力幅值，
j i v v v -=ij
，ij θ为高斯白噪声项，满足以
下统计学特性：
)(=t ij θ
)
'()()'()(t t t t jk il jl ik kl ij -+=δδδδδθθ (3)
C ω、
D ω、R ω为各力的权重函数，
由于热平衡(恒温系统)的需要，其形式一
般采用如下形式：
⎪⎩
⎪⎨⎧><-===)(,0)(,)/1()()(2
2
2c ij c ij c ij R D C R r R r R r ωωω (4)
其中c R 为截断半径，这3种力都只在两粒子相距小于c R 时存在，都是短程力。

一般选取c R 的无量纲长度为1。

Espanol 和Warren 22将耗散涨落定理引入耗散粒子动力学方法中，为耗散力和随机力提供了统计热力学上的平衡条件。

DPD 模拟算法的首次提出是由
Hoogerbrugge 和Koelman 14提出的，他们采用了欧拉算法，属于一阶时间精度算法。

为了改进这种算法，Pagonabarraga 等人
23
在蛙跳算法的基础上提出了一种自洽算法。

还有另外一种积分算法是由Groot 等人24提出的修正的Velocity-Verlet 算法：
))
()((2/1)()())(~),(()()()()(~)
()(2/1)()()(2t t f t f t t v t t v t t v t t r f t t f t tf t v t t v
t f t t tv t r t t r i i i i i
i i i
i i i i i i ∆++∆+=∆+∆+∆+=∆+∆+=∆+∆+∆+=∆+λ （5）
其中，λ为可调参数，当5.0=λ时和标准的Velocity-Verlet 算法
25
一致。

目前
绝大多数DPD 模拟都采用这种修正的Velocity-Verlet 算法，本文也将采用这种算法。

在DPD 体系中我们将单个DPD 粒子
的质量i m 、截断半径c R 和能量T k B 选取为系统的参考单位。

如果以水为例，将m N 个水分子粗粒化为一个DPD 粒子，得到一个DPD 粒子的质量为O H m i m N m 2⋅=，继而可以得到无量纲单位1对应的真实物理长
度
：
3
*2
2
O
H O
H m DPD m N L ρρ⨯⨯= （6）
其中，DPD ρ为DPD 粒子的数密度。

通过长度、质量、能量三个无量纲基
本量，我们可以分别导出速度、时间和力的无量纲量表达式：
第八届全国流体力学学术会议 2014年9月18-21日 甘肃 兰州
CSTAM2014-A26-BS10029
1）基金资助项目 *
*****
*////L T k F M T k L t M T k V B B B === （7） 其中，B k 是玻尔兹曼常数，一般地来说，0.1=T k B ；m N M m **
=，即单个DPD 粒子的质量。

如无特别说明，在本文接下来的内容中的模拟中采用的参数均以DPD 无量纲单位给出，各种参数与物理单位的转化在具体算例中不再逐一说明。

2.2 三组分系统下的固固作用系数
在20世纪40年代，Paul Flory 26和Maurice Huggins 27分别独立地研究出了基于格子模型的可以应用于聚合物溶液的理论，叫做Flory-Huggins 模型。

简单地来说，就是推出了一个化学上的混合物聚合物的熵增的公式：
)ln ln (2211x x x x R S m +-=∆
其中，R 是理想气体常数，1x 是溶剂的摩尔分数，
2
11
1n n n x +=
其中1n 和2n 分别对应溶剂和溶质所占据的格子数目。

我们也可以把
)ln ln (2211x x x x R S m +-=∆改写成)ln ln (2211φφx n k S m +-=∆，其中1φ和
2φ分别代表物质1和物质2的体积分数。

对高分子聚合物来说，同理，该式可以改成)ln ln (2
11∑=+
-=∆N
i i
i
m n n k S φφ。

对于焓的变化，有
1221χφkTN H m =∆ （8）
又因为自由能的表达式
m m m S T H G ∆-∆=∆，将上述两个表达式
代入可得，对于双组分系统，有如下的公式：
)ln ln (21122211φχφφn n n kT G m ++=∆ （9）
1997年，Robert D. Groot 和Patrick B. Warren 24建立了一个连接DPD 和宏观的Gibbs 自由能的桥梁，即
)(B A V f F ρρ+=。

Toshiyuki Kataoka 、Yohei Nagao 28等人在2007年指出了三组分情况下自由能密度的表达式：
S
B BS S A AS B A AB S S B B
B
A A
A
m m m N N f kTf F φφχφφχφφχφφφφφφ+++++
=
Ω=ln ln ln （10）
其中，T 是绝对温度，k 是玻尔兹曼常数，Ω是总格子数，i φ是第i 中组分的体积分数，ij χ是Flory-Huggins 作用系数。

经过一系列推导之后，我们可以得到在三组分系统当中自由能密度的两种表达形式：
T
k a a a a a a N N N N f N N f B B S BS S A AS B A AB S SS B BB A AA S B B
A A S S
B B B A A A m S B BS S A AS B A AB S S B B
B A A A m )222(ln ln ln (*)ln ln ln 222ρρρρρρρρραρρρρρρρρρφφχφφχφφχφφφφ
φφ++++++
---++=+++++=
根据上面的两个式子，只需令上式=下
式，我们也可以推出在三组分系统中的固
固作用系数（SS a ）的表达式，即：
]
)222(ln ln ln [*2
2
2T
k a a a a a N N N N T k a B B S BS S A AS B A AB B BB A AA S B B A A S S B B B A A A S B SS ρρρρρρρραρρ
ρρρρρρραρ+++++---++-=
其中，星号*代表第一个式子的等号右半
边。

以上是关于作用系数的理论研究情
况。

本文主要着眼于数值模拟，通过调整不同参数来揭示不同物质的作用系数对系统的一些影响。

3.具体算例 3.1 建模
基本的模型图如下：
图1
从上图可以看出，系统由三种物质组成，分别是水（液滴）、固壁和悬浮在液滴中的二氧化硅纳米颗粒。

为了消除壁面附近的密度波动，本文使用了随机固壁模型。

具体建模过程为：先在长宽高分别为xyz 的体积为x×y×z 的长方体计算域中用面心立方排布充满DPD 固体粒子，经过第一
次热平衡之后通过规定粒子的位置确定粒子的种类和系统内各个组分的形状。

该过程我们称之为“熔铸”；熔铸完成之后得到的固壁即原来热平衡之后的固体DPD 粒子形成。

随机固壁可以有效地减弱密度波动。

我们的壁面形状较为简单，壁面的反弹条件设置同FCC （面心立方，Face
Centered Cubic ）固壁，是反弹边界条件
（Bounce Back ）。

图2
3.2 计算方法及理论
在传统DPD 模型中最重要的参数便是各种粒子间的保守力参数，该参数的求得可通过与Flory-Huggins 理论的对应得到，通过求解描述相互溶解程度的Flory-Huggins 参数，进而解得保守力参数。

标准的DPD 方法在介观模拟方面已经取得了很大的成功，而多体DPD （Multi-body Dissipative Particle Dynamics, MDPD ）方法则可以解决标准的DPD 方法所无法解决的一些问题，即由于在状态方程中缺失密度
的高次项使得该软排斥力的方法不能模拟气液共存的问题24,29。

在本文中，为了模拟气液共存的状态，必须使状态方程中出现密度的高次项，所以我们采用多体耗散粒子动力学（Multi-body Dissipative Particle Dynamics, MDPD ）来对液滴进行圆拟合。

MDPD 方法是在保留了耗散力、随机力的基础上，通过引入吸引相互作用，对保守力进行了修正，修正后的保守力的表达形式为
ij ij d j i ij ij ij C ij ij C e r w e r w F
)()()(ρρβα++= （11）
其中，0<ij α，0>ij β。

上式中，
⎩⎨
⎧>≤-=)(,0)(,/1)(c ij c ij c ij ij C r r r r r r r ω，⎩
⎨⎧>≤-=)(,0)
(,/1)(d ij d ij d ij ij d r r r r r r r ω （12）
2
3
)/1(215)(,)(d ij d ij i
j ij i r r r r r -==∑≠πωωρρρ （13）
式中，c r 和d r 分别是吸引和排斥截断半径。

P.B.Warren 29改进了状态方程的形式，
得出了
)
(22342d c Br A T k p d B +-++=ρραραρ，并做了一系列的数值实验，选取了两组
参数做了进一步的研究，并最终确定了系数矩阵A 和B 的具体的数值。

在确定一些其他项的时候更多的是进行了人为地、经验地修正。

3.3 参数设置
计算域设为30×30×20的立方体空间，x 、y 方向使用周期性边界条件，在z
方向上下边界为熔铸后切割出来的固壁，厚度为1.8，计算域是30×30×（20-1.8）的长方体空间,熔铸出液滴和纳米颗粒，液滴下表面使用反弹边界条件将进入固壁的水弹回。

液滴初始位置设于z=0处，半径Rball 为7，构筑方法同为熔铸。

下表是详细的关键参数列表：
表1
模拟中所涉及的两种物质间的保
守力的排
斥参数为25ij =B ，具体的吸引力参数ij A 为下表：
表2
其中，当i=j 时，也就是对角线上的作用系数都是-40；水和二氧化硅纳米颗粒的吸引力系数采用的是Materials Studio 软件中的blend 模块进行计算的结果，具体我们选用COMPASS 力场，温度设定为333K 。

水和固壁的作用系数是对应初始接触角为70°的
ij A
（对于本课题组自主开发的圆拟合测量接触角的程序，相对而言测量比较大的角度更加精确，并且该结果与上海大学的狄勤丰课题组的实验结果相对应）；
ss A ，即solid solid A ，是固壁和二氧化硅纳
米颗粒的吸引力作用系数，或称为固固作用系数，是本文需要研究的对象之一。

3.4 热平衡步数的选择
在正式计算之前，我们需要把通过人为设定的初始条件给出的系统进行热平衡，以最大程度地消除人为因素带来的影
响。

首先，我们需要确定热平衡的步数。

如果热平衡步数太少，那么理论上来说可能不会基本消除人为因素带来的影响；如果热平衡步数太多，则浪费计算时间。

在建模过程中，我们采取了二次热平衡的步骤。

第一次热平衡：为了实现随机固壁，我们采用熔铸的方法。

具体地来说，就是先在计算域中充满DPD 固壁粒子，然后在第一次平衡步的时间内进行随机热平衡，达到稳定之后截取随机固壁。

图3是经过1000步热平衡步之后得到的dat 文件通过可视化软件Tecplot 查看得到。

一般来说，不同大小的计算域需要的热平衡步也不一样，一次热平衡的平衡步在8000~10000步之间。

保险起见，我们选择10000步的平衡步。

图3 在长方体计算域中，充斥着固壁粒子
第二次热平衡：第一次热平衡结束、熔铸好了固壁和液滴并添加了二氧化硅纳米颗粒之后，需要进行二次热平衡，使得二氧化硅纳米颗粒均匀地分布在液滴中。

作者在二次热平衡的时候也设置了10000步的平衡步，下面两幅图分别是step=400和step=6000时打出的图像。

图4蓝色的粒子为二氧化硅纳米颗粒，黑色的粒子为固壁，绿色的粒子为液滴
上面两幅图分别是step=400和
step=6000时的图。

可以看出，在二次热平
衡时间不长不够大的时候，二氧化硅纳米
颗粒都飞出了计算域；在二次热平衡的步
数足够多的时候，二氧化硅纳米颗粒都悬
浮到了液滴中（二氧化硅纳米颗粒被液滴
挡住导致看不到）。

为了确保足够均匀地
分布在了液滴中，我们还是采用前人采用
的10000步作为二次热平衡的步数。

综上所述，一次热平衡和二次热平衡
都采用10000步的时间步。

3.5 参数取值范围
在数值计算中，每个参数都有其合理
的取值范围。

同理，固固作用系数也有其
适用范围。

MDPD的理论规定，0
ij
A，
其中也包括固固作用系数，但我们需要得
到更加精确的
ij
A的取值范围为后来的计算
模拟等工作做准备。

从之前的研究中可以得出，
12
A的选取
只会影响初始的接触角。

如果从理论上进
行分析，这个结论是可以接受的，因为
12
A
的含义就是水和固壁的相互作用，与纳米
颗粒无关。

虽然在整个系统内部，三个组
分会两两相互作用，但是就
12
A而言，我们
尚未研究其与纳米颗粒的位置、速度等物
理量之间的关系。

接下来的两幅图从可视化的角度说明
了改动参数之后会出现的“悬浮在液
滴”、“沉降在固壁”这两种现象。

图5
为了更加清楚直观地观察二氧化硅纳
米颗粒的分布，我们取消了DPD水粒子的
显示，只显示固壁和纳米颗粒。

可以从以
上两幅图中清楚地看到，左边的图中，纳
米颗粒是悬浮在液滴中的；而在右边的图
中，绝大多数纳米颗粒都沉降在固壁。

正
是固固作用系数的不同取值，才出现了上
两幅图的现象。

通过大量的算例可以得到，对于这种固壁——水——二氧化硅纳米颗粒的三组分系统来说，固固作用系数取值在合理区间内的话，呈现的就是如左图的样子；如果固固作用系数超过了一定范围的话，便会得到右图的情况。

如果纳米颗粒的物质种类变化了，即13A 的值产生了变化，那么
ss A 的取值范围也会发生变化。

举个例子，
如果13A 的值从-34.8（水和二氧化硅纳米颗粒的作用系数）变成-50的话，那么相应的
ss A 的取值范围就会变成
530-<<-ss A 。

如果30-<ss A 的话，
那么右图的情况也会出现。

一般地来说，
我们有5152313-≤≤+A A ；对于本次模拟来说，我们有51823-≤≤-A 。

3.6 标准MDPD 下参数取值对接触角
的影响
确定了参数的取值范围之后，接下来我们需要研究的就是参数取值对接触角的影响。

所谓接触角，是指液体在水平固体表面铺展达到平衡时，在气、液、固三相交点处所作的气-液界面的切线穿过液体与固-液交界线之间的夹角，如下图所示。

从物性上来说，接触角的大小通常由固-液、液-气、固-气之间表面张力决定。

当液滴在理想表面上达到三相平衡时，可由Young 方程（下式）给出各相关表面张力与接触角
的
函
数
关
系。

lv sl sv γγγθ/)(cos -=
图6 接触角示意图
其中，sv γ为固体-气体间的表面张力，sl γ为固体-液体间的表面张力，lv γ为液体-气体间的表面张力，θ为液、固、气三相平衡时的接触角。

通常， 0=θ时液体完全润湿固体表面，固体表面具有超亲水性；
900<<θ时液体部分润湿固体表面，固体表面具有亲水性；
180
90<<θ时液体不能润湿固体表面，固体表面具有
疏水性，
150>θ时，固体表面具有超疏
水性。

在计算机模拟中，我们采用“圆拟合”的方法求液滴静态表观接触角。

首先，根据粒子所处位置的密度来确定界面，并根据界面上所有的粒子拟合球面（对于二维问题为圆），然后根据接触角的定义求得接触角。

下图为液滴表面水粒子
的
示
意
图。

图7
通过提取出的液滴表面粒子，进行整体球拟合（2D 状态下为圆拟合），然后再求出接触角。

30
3.6.1 固固作用系数与接触角的关系
在科学研究当中，为了研究某个点，我们需要考虑的就是暂时固定其他可以变动的因素并尽量消除一些外在因素的影响，集中考虑需要研究的问题。

在本文中，要研究固固作用系数，我们就需要避免牵扯除了固固作用系数之外的因素。

在
数不多，可以变动的只有12A 、23A 和13A ，分别对应的是水和固壁的吸引作用系数、固壁和纳米颗粒的吸引作用系数以及水和纳米颗粒的吸引作用系数。

其中，12A 和接触角的关系已经由黄汝佳等人给出31
；水和二氧化硅的物质已经确定，13A 可以由Material Studio 计算得出。

所以我们采用固定12A 和13A 的方法，仅仅考察23A 和接触角的关系。

基本的参数参考之前规定的系数；下表3
从上表可以看出，23A 的变动对接触角的影响不大。

所以我们得到的结论之一就是，固固作用系数23A 对纳米流体液滴的接触角影响非常小。

3.6.2 浓度与接触角的关系
从上面的分析可以得知，接触角对固固作用系数的变化不是很敏感，由作用系数引起的接触角变化很小。

除了前人得出的固液作用系数可以显著地影响接触角大
小的结论外，对于纳米流体液滴，我们还需要考虑一些新的可以影响接触角的大小的因素。

S Vafaei 等人32
经过研究发现，在纳米流体的浓度为0.005‰到5‰左右的范围内，纳米流体的接触角与纳米颗粒的浓度、尺寸以及基板的性质有关。

具体的结果图示如
下
：
第八届全国流体力学学术会议 2014年9月18-21日 甘肃 兰州
CSTAM2014-A26-BS10029
1）基金资助项目
图8
其中，第一行的两幅图是在相同的基板下不同尺寸的纳米颗粒对应的纳米流体的接触角关于浓度的点图，第二行的是另外一种性质的基板下不同尺寸的纳米颗粒对应的纳米流体的接触角关于浓度的点图。

由于DPD 方法空间尺度和时间尺度的限制，暂时无法模拟浓度变化如此大的情况。

3.7 MDPD 参数 3.7.1 关于MDPD 系数取值的进一步探讨
P.B.Warren 29在2003年提出了MDPD 算法用来解决模拟气液共存遇到的困难，并给出了基本的假设。

对于气液共存的情况，我们需要令0<A 、0>B 以使得在物态方程EOS （Equation of State ）里面存在范德华循环，使得粒子成对吸引，并加入作用半径比吸引力更小的排斥力项。

EOS 在气液共存的时候可以被简化为如下的形式：
)(22342d c Br A T k p d B +-++=ρραραρ （14）
其中，)1(101.0=α，)2(16.4=c ，)1(18=d 。

当0)(=L p ρ的时候，也就是界面压力为0的时候可以得到气液界面。

根据下
图，可以估计出5~L ρ，并且30~~B A -。

其中纵坐标为等号的左边
p ，横坐标为密度。

图9
在这之后作者选出了两组数据进行分析对比，如下表。

图10
对于温度，根据之前挑选出的两组数据，作者给出了一组图：
图11
纵坐标为温度，横坐标为密度，分别对应的是0)(=L p ρ的点。

低温可以制造一个泾渭分明的界面，但是如果温度太低，界面靠近液体的那一侧会出现振荡（Oscillation ）。

所以1=T k B 。

根据以上的分析，可以得到一些标准的MDPD 数据：1=T k B ；75.0=d r ；40-=A 。

如果要选择数密度为6的话，
那么25=B 。

3.7.2 修改系数之后的算例
在计算过程中，作者发现了一个奇特的现象。

在二氧化硅纳米颗粒数大于某个值的时候，即cri SiO c N >2（cri c 为一个特定值，大约等于50）的时候，整个系统的平衡会遭到破坏，出现如下图的现象。

图12
在这四幅小图里，左边两幅图是破碎的体系的完全展示；右边两幅图是不显示水的时候的样子。

可以看到，整个固壁——大液滴——液滴中悬浮的二氧化硅纳米颗粒的体系变得支离破碎，不过有趣的是二氧化硅纳米颗粒仍然都均匀地分散在各自的小系统里面。

因此，当纳米颗粒浓度达到一定程度时，悬浮在液滴中的二氧化硅纳米颗粒对液滴与液滴之间的作用影响已不可忽视，或者说是水的DPD 粒子之间吸引力不够大or 排斥力过于大才导致了如此现象的出现，之前所提及的标准参数的取值需要改进。

所以，我们调整了之前文献中规定的参数（set2，25,40=-=ij ii B A ）并进行了新的计算，得到了更加合理的结果。

事实上，无论是从MDPD 的理论上来看还是从实际情况来看，对于任意一对粒子（假设为i 粒子和j 粒子）来说，周围的环境变化都会对这对粒子的相互作用产生一定的影响。

我们所说的i 粒子和j 粒子之间的相互作用ij A 也是不考虑多体作用的成对相互作用假设下的做法。

我们寻找ij ii B A ,中比较敏感的系数。

如果单纯地调节吸引力系数11A ，得到的结果
如
下
：
图13a
图13b
图13c
上图是在别的参数都不变的情况下，把11A 调小得到的两个结果，其中图13a 是
4511-=A 时的结果，图13b 是50
11-=A 的结果，图13c 是6011-=A 时的结果。

如果单纯地调节排斥力系数11B 的话，得
到
的
结
果
如
下
：
图14a 图14b 图14c
在上图中，图14a 是1511=B 时的情
况；图14b 是2011=B 时的情况；图14c。