基于分子动力学模拟的纳米液滴接触角获取方法与流程

本发明属于测量技术领域，提供了一种基于分子动力学模拟的纳米液滴接触角获取方法。

背景技术：

液滴润湿涵盖多个尺度(宏观、介观以及微观)，是基础研究和工程应用的重点研究方向，例如农药沉积、飞机机身防冰、内燃机内燃油黏壁以及微纳米机器的黏附失效等问题。在农药喷雾中，农药液滴越亲水，则农药与植物叶面的接触面积越大，除虫效果越明显。而对于飞机机身防冰、内燃机内燃油黏壁以及微纳米机器的黏附失效，液滴若能在壁面静止时呈现疏水或者超疏水形貌，则液滴在撞击壁面过程中容易发生反弹并离开壁面，这将有利于工程需求。近年来，采用分子动力学模拟方法,利用large-scaleatomic/molecularmassivelyparallelsimulator(lammps，一个分子动力学的模拟软件)软件研究纳米液滴在固体壁面上的润湿特性、自组织行为以及固体表面流体的凝结换热问题受到大批学者的关注。

液滴的静态润湿性采用其在固体表面所形成的接触角(该角是指气-固-液三相交界处，固液界面到液气表面所经过液滴内部的夹角)来区分。对于宏观液滴常采用接触角测量仪截取液滴平衡时的图片，确定三相接触线的液气界面，然后利用接触角测量仪内置软件选取切线，读出相应的角度即可。另一方面也有部分学者采用相场动力学数值方法捕捉相界面，然后利用界面拟合法或者界面切线法得到接触角。然而，随着液滴尺度的缩小，一方面受制于人眼生理的观察限制，另一方面得益于利用分子动力学研究润湿性的日益发展，对于纳米液滴数值仿真时纳米液滴静态接触角的获取还没很好的处理方法。大多数研究者常采用截取纳米液滴静态图片，利用软件量取角度的方法获取。此外，也有部分学者考虑液气界面存在一定的厚度，并且在接近固壁处液滴分子会形成“固溶区”，这两个因素将极大影响纳米液滴接触角的测量。综上所述，对纳米液滴静态接触角的测量既有人为因素影响又有其本身固有的难度，这导致现有方法测量误差较大，数据可信度较小。

技术实现要素：

本发明实施例提供一种基于分子动力学模拟的纳米液滴接触角获取方法，旨在解决人为测量静态接触角存在的测量误差大，数据可信度低的问题。

本发明是这样实现的，一种基于分子动力学模拟的纳米液滴接触角获取方法，所述方法具体包括如下步骤：

s1、利用lammps软件输出液滴分子平衡时的位置和数密度，并导入tecplot软件；

s2、tecplot软件将导入的数密度分成不同的数密度区域，将数密度中值所在的数密度等值线作为液气表面和固气表面，提取位于所述数密度等值线上的数据点，数据点由湿润长度及液滴高度组成，将提取的数据点导入matlab；

s3、利用matlab多项式拟合法对导入的数据点进行液气拟合曲线及固气拟合曲线拟合，并确定两拟合曲线的交点；

s4、确定上述交点在液气拟合曲线及固气拟合曲线上的倾斜角，两倾斜角差值的绝对值即为液滴的静态接触角。

进一步的，液滴分子在平衡时的位置和数密度获取方法具体如下：

s11、基于lappms软件模拟液滴在基底上润湿时，利用郎之万热浴法使液滴驰豫50ps；

s12、之后每隔1ps提取一组液滴分子的位置和数密度数据，连续提取l组；

s13、针对每组位置和数密度数据分别执行步骤s1至步骤s4，获取l个静态接触角，将l个静态接触角的平均值作为静态接触角进行输出。

进一步的，液气拟合曲线及固气拟合曲线的获取方法具体如下：

s31、确定多项式阶数s，4≤s≤6，固气拟合曲线及液气拟合曲线分别为p(x)和q(x)具体表达式如下：

s32、利用最小二乘法拟合多项式系数向量a＝[a0,a1,...,am]和b＝[b0,b1,...,bn]，其中，n及m取值为s。

本发明提供了一种的可信的液滴静态接触角的测量方法，同时通过捕捉合适的液气和固气界面使得计算液滴接触角更加可信，解决了人为测量静态接触角存在的测量误差大，数据可信度低的问题。

附图说明

图1为本发明实施例提供的分子动力学模拟的纳米液滴接触角获取方法流程图；

图2为本发明实施例提供的tecplot处理后是数密度云图；

图3为本发明实施例提供的matlab处理后的固气和液气拟合曲线图相交图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

图1为本发明实施例提供的分子动力学模拟的纳米液滴接触角获取方法流程图，该方法具体包括如下步骤：

s1、利用lammps软件中输出液滴分子在平衡时的位置和数密度数据，存入txt文本中，然后添加相应的文件头，导入进tecplot中；

s2、tecplot将数密度数据划分不同的数密度区域，将数密度中值所在的数密度等值线作为为液气表面和固气界面，提取位于所述数密度等值线上的数据点(xi,yi)，xi为湿润长度，yi为液滴高度组成，如图2所示，将提取的数据点存为txt文本后导入到matlab中；

s3、基于matlab多项式拟合法对导入的数据点进行液气拟合曲线和固气拟合曲线的拟合，确定两拟合曲线的交点；

s4、定上述交点在液气拟合曲线及固气拟合曲线上的倾斜角，两倾斜角差值的绝对值即为液滴的静态接触角。

在本发明实施例中，步骤1中输出相应的液滴分子在平衡时的位置和数密度数据，具体操作包含以下操作：

s11、使用lappms软件模拟液滴在基底上润湿时，利用郎之万热浴法使液滴驰豫50ps。然后，在in文件后添加以下命令：

computecc1watchunk/atombin/cylinderzlower2000.070.035discardyesunitsbox；

fix6watave/chunk100011000cc1density/numberaveonefilevel.profile；

computecc1表示对计算盒子进行分割取片，计算数密度；fix6表示输出坐标信息和数密度信息；

s12、之后每隔1ps提取一组液滴分子的位置和数密度数据，连续提取l组；l组位置和数密度数据作为可信数据，导入tecplot中进行数据的处理，有利于降低随机性，取平均后增加数据可信度；

s13、针对每组位置和数密度数据分别执行步骤s1至步骤s4，获取l个静态接触角，将l个静态接触角的平均值作为静态接触角进行输出。

在本发明实施例中，步骤2中txt文本中添加相应的文件头具体操作包含以下操作：

s22、tecplot具有官方的数据文件格式，导入外部数据时数据源文件必须按照相应的格式填写，否则tecplot提示出错。具体文件头表达如下：

title＝"datafornumberdensity",表示题目，一般用作提示作用；

variables＝"x","y","p"表示变量名称，x为横坐标(液滴湿润长度)，y为纵坐标(液滴高度)，p为数密度(用来进行云图显示)；

zonei＝65,j＝40,f＝point为控制行，定义了数据的排列格式,例如i＝65，j＝40表示数据是按ij顺序组织的，其中指标i从1到65变化，指标j从1到40变化,i和j的值应当根据实际情况选择；

s23、导入tecplot后选择相应的云图显示即可。

在本发明实施例中，步骤3中的具体操作包含以下操作：

s31、定义液气拟合曲线和固气拟合曲线的拟合的数据名称，令固气拟合曲线横坐标和纵坐标分别为x1＝[x10,x11,...,x1n]、y1＝[y10,y11,...,y1n]；令液气拟合曲线横坐标和纵坐标分别为x2＝[x20,x21,...,x2m]、y2＝[y20,y21,...,y2m]，n和m不一定相等；

s32、在matlab中新建m(脚本)文件，在进行曲线拟合之前需要对原数据进行绘图，在该m文件中具体指令如下：

plot(x1,y1,'r*')；％以红色散点图绘出固气等值线；

xlabel('x轴')；ylabel('y轴')；title('散点图')；holdon；％添加坐标轴信息；

plot(x2,y2,'b*')；holdon；％以蓝色散点图绘出液气等值线；

s33、确定多项式阶数s，(一般选择6阶多项式就有相当高的精度要求，即s＝6)，则固气拟合曲线及液气拟合曲线分别为p(x)和q(x)，其具体表达式如下：

s34、利用最小二乘法拟合多项式系数向量a＝[a0,a1,...,a6]和b＝[b0,b1,...,b6]，以求多项式系数向量a为例，在该m文件中具体指令如下：

n＝t；s＝6；％t为n的具体数值；

a＝zeros(s+1)；％初始化系数矩阵；

forj＝1:s+1

fori＝1:s+1

fork＝1:n+1

a(j,i)＝a(j,i)+x1(k)^(j+i-2)

end

end

end；

b＝zeros(1,s+1)；％初始化非齐次向量矩阵；

forj＝1:s+1

fori＝1:n+1

b(j)＝b(j)+y1(i)*x1(i)^(j-1)

end

end；

s35、写出正规方程，求出多项式系数向量a＝[a0,a1,...,a6]，在该m文件中具体指令如下：

b＝b'；

a＝inv(a)*b；

s36、采用步骤s34和s35求出多项式系数向量b；

步骤307：画出固气拟合曲线和液气拟合曲线，在该m文件中具体指令如下：

z1＝a(0)+a(1)*x1+a(2)*x1.^2+a(3)*x1.^3+a(4)*x1.^4+a(5)*x1.^5+a(6)*x1.^6；

plot(x1,z1,'r')；holdon；％以深灰曲线图绘出固气拟合曲线

z2＝b(0)+b(1)*x2+b(2)*x2.^2+b(3)*x2.^3+b(4)*x2.^4+b(5)*x2.^5+b(6)*x2.^6；

plot(x2,z2,'b')；holdon；％以浅灰曲线图绘出液气拟合曲线

步骤310：取三次平均，减少随机性，将平均值作为液滴的最终静态接触角。

所述利用最小二乘法拟合多项式系数向量的方法具体如下：

步骤401：对于给定的n+1个数据点(xi,yi)，假设取m(通常取6)阶拟合多项式p(x)使得平方和i(a0,a1,...,am)最小，其中p(x)和i(a0,a1,...,am)的具体表达如下：

其中a0,a1,...,am为待求的系数未知数。

步骤402：为了平方和i(a0,a1,...,am)最小，根据多元函数求极值问题，则该问题化为以下表达式：

进而有：

用矩阵表示如下：

简写矩阵为aa＝b，其中a为m×n阶矩阵，a和b均为n维列向量；只要给出数据点(xi,yi)和拟合多项式的最高阶数m，编写m文件程序就可以得到相应的拟合多项式系数向量。

步骤308：确定两条拟合曲线的交点，在该m文件中具体指令如下：

symsxy

s＝solve(y＝＝a(0)+a(1)*x+a(2)*x^2+a(3)*x^3+a(4)*x^4+a(5)*x^5+a(6)*x^6,y＝＝b(0)+b(1)*x+b(2)*x^2+b(3)*x^3+b(4)*x^4+b(5)*x^5+b(6)*x^6,x,y)；％利用slove求解交点

x＝double(s.x)；y＝double(s.y)；

text(x,y)；％在图中显示该点

步骤309：利用倾斜角之差，得到液滴的静态接触角，在该m文件中具体指令如下：

p’＝ployder(a)；％求固气拟合曲线的导数系数向量

q’＝ployder(b)；％求液气拟合曲线的导数系数向量

sita1＝atan(ployval(p’,x))；％求固气拟合曲线的在交点处的倾斜角

sita2＝atan(ployval(q’,x))；％求液气拟合曲线的在交点处的倾斜角

sita＝sita1-sita2；

本发明提供了一种的可信的液滴静态接触角的测量方法，同时通过捕捉合适的液气和固气界面使得计算液滴接触角更加可信，解决了人为测量静态接触角存在的测量误差大，数据可信度低的问题。此外，图3为在matlab中进行处理后的固气和液气拟合曲线图以及相应的交点。经过数据处理后的接触角测量值为104°，而实际接触角为102±2°，因此本发明对测量液滴在壁面上的接触角具有很好精度值，数据可信度高，人为因素少，随机性小。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。