纳米流体液滴的耗散粒子动力学方法模拟Word格式.docx
《纳米流体液滴的耗散粒子动力学方法模拟Word格式.docx》由会员分享，可在线阅读，更多相关《纳米流体液滴的耗散粒子动力学方法模拟Word格式.docx（26页珍藏版）》请在冰豆网上搜索。
通常，悬浮在流体中的纳米颗粒会受到诸如流动阻力、布朗运动、粒子间扩散及重力等内外因素的影响，其运动规律极其复杂。
1
前人对于纳米流体性质的各个方面做了广泛的研究。
1993年日本东北大学的Masuda等人2在水中添加平均粒子直径为13nm的
和27nm的
粒子，制备了不同体积浓度的纳米颗粒胶体并测量了胶体的导热系数;
1995年，美国Argonne国家实验室的Choi等人以一定的方式和比例在液体中添加纳米级金属或金属氧化物粒子，并称之为纳米流体。
3其中的氧化物粒子包括
、
等，另外还有一些金属粒子和碳化物等4。
在1995年Choi3之后，国内外的学者纷纷对纳米流体展开了深入的研究。
在国内方面，范庆梅等人5对纳米流体的热导率和粘度进行了计算。
李云翔等人6为纳米流体的研究进展做了一个总结，包括：
纳米流体稳定性的研究、纳米流体物性的研究、纳米流体传热特性的研究，其中既包括实验方面的研究进展也对纳米流体物性以及传热特性的理论研究进行了系统的总结。
国外方面，B.Davidovitch等人7研究了考虑了热扰动的粘性液滴在基板上的过程。
杰出的研究人员如Kim8-11、Vassallo12、Truong13等人也在纳米流体的性质研究方面做出了杰出的贡献，这里不再一一赘述。
分子动力学模拟（MolecularDynamics,MD）是一种在微观尺度下模拟原子和简单分子运动的方法。
但由于超出这个尺度的时候，MD方法不适用，我们需要采用另外的方法。
对于介于宏观和微观之间的尺度上的流体动力学行为的研究是目前学术界的热点问题。
这个中间的尺度通常被称为介观尺度，通常指的是10-1000nm和1ns-10ms的尺度。
本文所采用的介观尺度下模拟流体的动力学行为的方法是耗散粒子动力学方法（DissipativeParticleDynamics，简称DPD）14。
此方法通过对模拟区域内的粒子进行粗粒化，以减少计算代价，在更短的计算时间内计算更大的区域。
DPD方法中的相互作用系数（DPDinteractionparameter）是由内在的原子之间的相互作用决定的，在DPD方法中表现为DPD粒子之间的相互作用。
AmiteshMaiti等人15提出，DPD的作用系数和DPD粒子的尺寸有关。
接触角是衡量材料表面性能的一个重要参数，通过测量接触角可获得不同物质在界面之间相互作用的很多信息。
国内外的学者们对接触角的测量以及性质也已经有了一些成熟的研究。
DeckerE等人16对宏观和微观的接触角提出了合适的定义，并且给出了部分测量接触角的方法。
MarmurA等人17总结了测量接触角方法的一些缺陷，并指出了如何规避这些缺陷。
SaeidVafaei等人18更是研究出了与体积无关的更加可以准确反映物性的渐近接触角。
宁乔等人19利用照相机捕捉到的液滴灰度图像，利用边界提取等方法得到液滴的轮廓，采用轮廓拟合的方法进而求得液滴的接触角。
徐志钮等人20认为，如何准确测量接触角对相关学科如材料、医药、半导体、生命科学、油墨工业、电力系统等领域都具有重要意义，并指出，当接触角较小时需要使用基于圆的拟合算法；
否则使用基于椭圆的算法。
对接触角的研究文献很多，在这里不再一一赘述。
在本文中，我们主要用到的测量接触角的方法是测量表观接触角，即通过提取液滴表面粒子进行球拟合，得到液滴在固体表面上的接触角。
具体做法是：
对于一个由粒子组成的液滴，首先根据粒子所处位置的粒子密度来判断粒子是否处于界面上，然后将所有处于界面上的粒子提取出来，找到液滴界面所在的球面，通过拟合得到准确的接触角。
除了测量表观接触角之外，以后还有可能会考虑采用测量局部接触角的方法。
局部接触角一般用来描述液滴的局部状态，表征液滴的动态接触行为。
获得局部接触角常用的方法是切线法：
首先根据液滴局部区域粒子密度来判断粒子是否处于界面上，然后提取该处界面上的粒子坐标，一般可采用三次多项式来拟合界面轮廓，求导后即可获得该时刻液滴的局部接触角。
2.耗散粒子动力学方法
2.1基本理论
在1998年，PBWarren21说明了DPD粒子的对象，可以指代分子、原子、高分子蛋白质聚团、带电离子团等等，指出了DPD方法的广泛的适用性。
DPD算法中粒子采用和MD算法一样的Lagrange坐标描述，体系中的粒子运动满足牛顿运动方程：
（1）
其运动方程中
为i粒子的位移和速度，
为i粒子受到的3种合力，包括颗粒间的有势作用力
，还包括类似于布朗动力学中的随机力
和耗散力
。
这些有势力都属于软势力：
（2）
其中，
即i、j粒子间的距离，
为和
方向相同的单位矢量，
为保守力系数，反映粒子i和j之间的最大排斥力，
和
分别为耗散力和随机力幅值，
，
为高斯白噪声项，满足以下统计学特性：
（3）
为各力的权重函数，由于热平衡（恒温系统）的需要，其形式一般采用如下形式：
（4）
其中
为截断半径，这3种力都只在两粒子相距小于
时存在，都是短程力。
一般选取
的无量纲长度为1。
Espanol和Warren22将耗散涨落定理引入耗散粒子动力学方法中，为耗散力和随机力提供了统计热力学上的平衡条件。
DPD模拟算法的首次提出是由Hoogerbrugge和Koelman14提出的，他们采用了欧拉算法，属于一阶时间精度算法。
为了改进这种算法，Pagonabarraga等人23在蛙跳算法的基础上提出了一种自洽算法。
还有另外一种积分算法是由Groot等人24提出的修正的Velocity-Verlet算法：
（5）
为可调参数，当
时和标准的Velocity-Verlet算法25一致。
目前绝大多数DPD模拟都采用这种修正的Velocity-Verlet算法，本文也将采用这种算法。
在DPD体系中我们将单个DPD粒子的质量
、截断半径
和能量
选取为系统的参考单位。
如果以水为例，将
个水分子粗粒化为一个DPD粒子，得到一个DPD粒子的质量为
，继而可以得到无量纲单位1对应的真实物理长度：
（6）
为DPD粒子的数密度。
通过长度、质量、能量三个无量纲基本量，我们可以分别导出速度、时间和力的无量纲量表达式：
（7）
是玻尔兹曼常数，一般地来说，
；
，即单个DPD粒子的质量。
如无特别说明，在本文接下来的内容中的模拟中采用的参数均以DPD无量纲单位给出，各种参数与物理单位的转化在具体算例中不再逐一说明。
2.2三组分系统下的固固作用系数
在20世纪40年代，PaulFlory26和MauriceHuggins27分别独立地研究出了基于格子模型的可以应用于聚合物溶液的理论，叫做Flory-Huggins模型。
简单地来说，就是推出了一个化学上的混合物聚合物的熵增的公式：
其中，R是理想气体常数，
是溶剂的摩尔分数，
分别对应溶剂和溶质所占据的格子数目。
我们也可以把
改写成
，其中
分别代表物质1和物质2的体积分数。
对高分子聚合物来说，同理，该式可以改成
对于焓的变化，有
（8）
又因为自由能的表达式
，将上述两个表达式代入可得，对于双组分系统，有如下的公式：
（9）
1997年，RobertD.Groot和PatrickB.Warren24建立了一个连接DPD和宏观的Gibbs自由能的桥梁，即
ToshiyukiKataoka、YoheiNagao28等人在2007年指出了三组分情况下自由能密度的表达式：
（10）
其中，T是绝对温度，k是玻尔兹曼常数，
是总格子数，
是第i中组分的体积分数，
是Flory-Huggins作用系数。
经过一系列推导之后，我们可以得到在三组分系统当中自由能密度的两种表达形式：
根据上面的两个式子，只需令上式=下式，我们也可以推出在三组分系统中的固固作用系数（
）的表达式，即：
其中，星号*代表第一个式子的等号右半边。
以上是关于作用系数的理论研究情况。
本文主要着眼于数值模拟，通过调整不同参数来揭示不同物质的作用系数对系统的一些影响。
3.具体算例
3.1建模
基本的模型图如下：
图1
从上图可以看出，系统由三种物质组成，分别是水（液滴）、固壁和悬浮在液滴中的二氧化硅纳米颗粒。
为了消除壁面附近的密度波动，本文使用了随机固壁模型。
具体建模过程为：
先在长宽高分别为xyz的体积为x×
y×
z的长方体计算域中用面心立方排布充满DPD固体粒子，经过第一次热平衡之后通过规定粒子的位置确定粒子的种类和系统内各个组分的形状。
该过程我们称之为“熔铸”；
熔铸完成之后得到的固壁即原来热平衡之后的固体DPD粒子形成。
随机固壁可以有效地减弱密度波动。
我们的壁面形状较为简单，壁面的反弹条件设置同FCC（面心立方，FaceCenteredCubic）固壁，是反弹边界条件（BounceBack）。
图2
3.2计算方法及理论
在传统DPD模型中最重要的参数便是各种粒子间的保守力参数，该参数的求得可通过与Flory-Huggins理论的对应得到，通过求解描述相互溶解程度的Flory-Huggins参数，进而解得保守力参数。
标准的DPD方法在介观模拟方面已经取得了很大的成功，而多体DPD（Multi-bodyDissipativeParticleDynamics,MDPD）方法则可以解决标准的DPD方法所无法解决的一些问题，即由于在状态方程中缺失密度的高次项使得该软排斥力的方法不能模拟气液共存的问题24,29。
在本文中，为了模拟气液共存的状态，必须使状态方程中出现密度的高次项，所以我们采用多体耗散粒子动力学（Multi-bodyDissipativeParticleDynamics,MDPD）来对液滴进行圆拟合。
MDPD方法是在保留了耗散力、随机力的基础上，通过引入吸引相互作用，对保守力进行了修正，修正后的保守力的表达形式为
（11）
上式中，
（12）
（13）
式中，
分别是吸引和排斥截断半径。
P.B.Warren29改进了状态方程的形式，得出了
，并做了一系列的数值实验，选取了两组参数做了进一步的研究，并最终确定了系数矩阵A和B的具体的数值。
在确定一些其他项的时候更多的是进行了人为地、经验地修正。
3.3参数设置
计算域设为30×
30×
20的立方体空间，x、y方向使用周期性边界条件，在z方向上下边界为熔铸后切割出来的固壁，厚度为1.8，计算域是30×
（20-1.8）的长方体空间,熔铸出液滴和纳米颗粒，液滴下表面使用反弹边界条件将进入固壁的水弹回。
液滴初始位置设于z=0处，半径Rball为7，构筑方法同为熔铸。
下表是详细的关键参数列表：
表1
参数
符号
大小
初始数密度
6.0
吸引力截断半径
1.0
排斥力截断半径
0.75
时间步长
0.01
系统温度
随机力幅值
3.0
耗散力幅值
4.5
调节因子
0.65
排斥力系数
（i,j=1,2,3）
25
吸引力系数
（i=1,2,3）
-40
模拟中所涉及的两种物质间的保
守力的排斥参数为
，具体的吸引力参数
为下表：
物质
水
固壁
二氧化硅纳米颗粒
-28.88
-34.8
表2
其中，当i=j时，也就是对角线上的作用系数都是-40；
水和二氧化硅纳米颗粒的吸引力系数采用的是MaterialsStudio软件中的blend模块进行计算的结果，具体我们选用COMPASS力场，温度设定为333K。
水和固壁的作用系数是对应初始接触角为70°
的
（对于本课题组自主开发的圆拟合测量接触角的程序，相对而言测量比较大的角度更加精确，并且该结果与上海大学的狄勤丰课题组的实验结果相对应）；
，即
，是固壁和二氧化硅纳米颗粒的吸引力作用系数，或称为固固作用系数，是本文需要研究的对象之一。
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的理论规定，
，其中也包括固固作用系数，但我们需要得到更加精确的
的取值范围为后来的计算模拟等工作做准备。
从之前的研究中可以得出，
的选取只会影响初始的接触角。
如果从理论上进行分析，这个结论是可以接受的，因为
的含义就是水和固壁的相互作用，与纳米颗粒无关。
虽然在整个系统内部，三个组分会两两相互作用，但是就
而言，我们尚未研究其与纳米颗粒的位置、速度等物理量之间的关系。
接下来的两幅图从可视化的角度说明了改动参数之后会出现的“悬浮在液滴”、“沉降在固壁”这两种现象。
图5
为了更加清楚直观地观察二氧化硅纳米颗粒的分布，我们取消了DPD水粒子的显示，只显示固壁和纳米颗粒。
可以从以上两幅图中清楚地看到，左边的图中，纳米颗粒是悬浮在液滴中的；
而在右边的图中，绝大多数纳米颗粒都沉降在固壁。
正是固固作用系数的不同取值，才出现了上两幅图的现象。
通过大量的算例可以得到，对于这种固壁——水——二氧化硅纳米颗粒的三组分系统来说，固固作用系数取值在合理区间内的话，呈现的就是如左图的样子；
如果固固作用系数超过了一定范围的话，便会得到右图的情况。
如果纳米颗粒的物质种类变化了，即
的值产生了变化，那么
的取值范围也会发生变化。
举个例子，如果
的值从-34.8（水和二氧化硅纳米颗粒的作用系数）变成-50的话，那么相应的
的取值范围就会变成
如果
的话，那么右图的情况也会出现。
一般地来说，我们有
对于本次模拟来说，我们有
3.6标准MDPD下参数取值对接触角的影响
确定了参数的取值范围之后，接下来我们需要研究的就是参数取值对接触角的影响。
所谓接触角，是指液体在水平固体表面铺展达到平衡时，在气、液、固三相交点处所作的气-液界面的切线穿过液体与固-液交界线之间的夹角，如下图所示。
从物性上来说，接触角的大小通常由固-液、液-气、固-气之间表面张力决定。
当液滴在理想表面上达到三相平衡时，可由Young方程（下式）给出各相关表面张力与接触角的函数关系。
图6接触角示意图
为固体-气体间的表面张力，
为固体-液体间的表面张力，
为液体-气体间的表面张力，
为液、固、气三相平衡时的接触角。
通常，
时液体完全润湿固体表面，固体表面具有超亲水性；
时液体部分润湿固体表面，固体表面具有亲水性；
时液体不能润湿固体表面，固体表面具有疏水性，
时，固体表面具有超疏水性。
在计算机模拟中，我们采用“圆拟合”的方法求液滴静态表观接触角。
首先，根据粒子所处位置的密度来确定界面，并根据界面上所有的粒子拟合球面（对于二维问题为圆），然后根据接触角的定义求得接触角。
下图为液滴表面水粒子的示意图。
图7
通过提取出的液滴表面粒子，进行整体球拟合（2D状态下为圆拟合），然后再求出接触角。
30
3.6.1固固作用系数与接触角的关系
在科学研究当中，为了研究某个点，我们需要考虑的就是暂时固定其他可以变动的因素并尽量消除一些外在因素的影响，集中考虑需要研究的问题。
在本文中，要研究固固作用系数，我们就需要避免牵扯除了固固作用系数之外的因素。
在标准MDPD方法中，可以变动的保守力系数不多，可以变动的只有
，分别对应的是水和固壁的吸引作用系数、固壁和纳米颗粒的吸引作用系数以及水和纳米颗粒的吸引作用系数。
和接触角的关系已经由黄汝佳等人给出31；
水和二氧化硅的物质已经确定，
可以由MaterialStudio计算得出。
所以我们采用固定
的方法，仅仅考察
和接触角的关系。
基本的参数参考之前规定的系数；
下表列出了
-5
-10
-15
-20
接触角
65.08
66.43
65.41
66.99
表3
从上表可以看出，
的变动对接触角的影响不大。
所以我们得到的结论之一就是，固固作用系数
对纳米流体液滴的接触角影响非常小。
3.6.2浓度与接触角的关系
从上面的分析可以得知，接触角对固固作用系数的变化不是很敏感，由作用系数引起的接触角变化很小。
除了前人得出的固液作用系数可以显著地影响接触角大小的结论外，对于纳米流体液滴，我们还需要考虑一些新的可以影响接触角的大小的因素。
SVafaei等人32经过研究发现，在纳米流体的浓度为0.005‰到5‰左右的范围内，纳米流体的接触角与纳米颗粒的浓度、尺寸以及基板的性质有关。
具体的结果图示如下：
图8
其中，第一行的两幅图是在相同的基板下不同尺寸的纳米颗粒对应的纳米流体的接触角关于浓度的点图，第二行的是另外一种性质的基板下不同尺寸的纳米颗粒对应的纳米流体的接触角关于浓度的点图。
由于DPD方法空间尺度和时间尺度的限制，暂时无法模拟浓度变化如此大的情况。
3.7MDPD参数
3.7.1关于MDPD系数取值的进一步探讨
P.B.Warren29在2003年提出了MDPD算法用来解决模拟气液共存遇到的困难，并给出了基本的假设。
对于气液共存的情况，我们需要令
以使得在物态方程EOS（EquationofState）里面存在范德华循环，使得粒子成对吸引，并加入作用半径比吸引力更小的排斥力项。
EOS在气液共存的时候可以被简化为如下的形式：
（14）
当
的时候，也就是界面压力为0的时候可以得到气液界面。
根据下图，可以估计出
，并且
其中纵坐标为等号的左边p，横坐标为密度。
图9
在这之后作者选出了两组数据进行分析对比，如下表。
图10
对于温度，根据之前挑选出的两组数据，作者给出了一组图：
图11
纵坐标为温度，横坐标为密度，分别对应的是
的点。
低温可以制造一个泾渭分明的界面，但是如果温度太低，界面靠近液体的那一侧会出现振荡（Oscillation）。
所以
根据以上的分析，可以得到一些标准的MDPD数据：
如果要选择数密度为6的话，那么
3.7.2修改系数之后的算例
在计算过程中，作者发现了一个奇特的现象。
在二氧化硅纳米颗粒数大于某个值的时候，即
（
为一个特定值，大约等于50）的时候，整个系统的平衡会遭到破坏，出现如下图的现象。
图12
在这四幅小图里，左边两幅图是破碎的体系的完全展示；
右边两幅图是不显示水的时候的样子。
可以看到，整个固壁——大液滴——液滴中悬浮的二氧化硅纳米颗粒的体系变得支离破碎，不过有趣的是二氧化硅纳米颗粒仍然都均匀地分散在各自的小系统里面。
因此，当纳米颗粒浓度达到一定程度时，悬浮在液滴中的二氧化硅纳米颗粒对液滴与液滴之间的作用影响已不可忽视，或者说是水的DPD粒子之间吸引力不够大or排斥力过于大才导致了如此现象的出现，之前所提及的标准参数的取值需要改进。
所以，我们调整了之前文献中规定的参数（set2，
）并进行了新的计算，得到了更加合理的结果。
事实上，无论是从MDPD的理论上来看还是从实际情况来看，对于任意一对粒子（假设为i粒子和j粒子）来说，周围的环境变化都会对这对粒子的相互作用产生一定的影响。
我们所说的i粒子和j粒子之间的相互作用
也是不考虑多体作用的成对相互作用假设下的做法。
我们寻找
中比较敏感的系数。
如果单纯地调节吸引力系数
，得到的结果如下：
图13a图13b
升级会员 