用于分子动力学模拟的生成粗糙表面并覆膜的建模方法

1.本发明属于分子动力学技术领域，具体涉及一种生成粗糙表面并覆膜的建模方法。


背景技术：

2.分子动力学模拟是研究如摩擦学等领域时常用的数值模拟方法，其中软件lammps是进行分子动力学模拟时常用的工具，该软件可以很好地模拟模型中各原子间相对运动时的力学特性。在进行摩擦学研究时，常见的分子动力学模型主要由压头和基体组成，其中基体为长方体结构，表面主要是理想的光滑表面，少数模型表面为有规则的粗糙表面或二维随机粗糙表面，但实际生产中，不可能生产出上述理想表面，而是有一定粗糙度的三维随机粗糙表面，因此有必要研究基体表面为三维随机粗糙表面时的各种力学性能。在零件表面附着润滑膜是减小相对运动机构摩擦磨损的常用方法，实际生产中，润滑膜主要通过溅射、喷涂等方式附着到基体表面，而现有分子动力学模型中很少将润滑膜附着于三维随机粗糙表面，因此有必要研究润滑膜附着于三维随机粗糙表面时的摩擦磨损问题，即需要建立基于上述实际问题的分子动力学模型，以便后续有关研究。
3.现有基于分子动力学的摩擦学研究主要将基体表面简化为光滑表面，少数模型会将基体表面简化为三维规则粗糙表面或二维随机粗糙表面，如2018年发表于期刊《journal of physics d:applied physics》第51卷第43期的文章《lubricity of graphene on rough au surfaces》，该文章所建立的模型由石墨烯、基体和压头组成，其中基体表面为二维随机粗糙表面，石墨烯附着于基体表面，研究了石墨烯附着于二维随机粗糙表面时的润滑性能，但实际工作条件下，基体表面应为三维随机粗糙表面。为克服上述不足之处，本发明建立的模型为长方体结构，模型由基体和润滑膜组成，基体表面为三维随机粗糙表面，润滑膜附着于基体表面，为实现从纳米尺度研究有关粗糙表面摩擦磨损机理奠定了基础。


技术实现要素：

4.为了克服现有技术的不足，本发明提供了一种分子动力学模拟的生成粗糙表面并覆膜的建模方法，模型为长方体结构，由基体和润滑膜组成，基体表面为三维随机粗糙表面，润滑膜附着于基体表面。主要包括以下步骤：步骤一：输出表面为三维随机粗糙表面的模型文件；步骤二：建立含有随机粗糙表面和润滑膜的分子动力学模型；步骤三：压头下压，使润滑膜附着于基体表面；步骤四：输出基体和润滑膜的原子信息，建立目标模型。与传统的分子动力学模型相比，本发明更加符合实际，操作简单，为实现从纳米尺度研究有关粗糙表面摩擦磨损机理奠定了基础。
5.本发明解决其技术问题所采用的技术方案包括如下步骤：
6.步骤1：输出表面为三维随机粗糙表面的模型文件；
7.步骤1-1：利用高斯随机表面生成程序生成一组包含坐标信息的随机点，设随机点所组成的集合为a，其中集合a中点的z坐标为复数；
8.步骤1-2：取集合a中每个点z坐标的实部，得到纯实部集合b1，依据如下公式(1)计算得到集合b：
9.zb(xi,yj)＝h0+k(z
b1
(xi,yj)-z
b0
)
ꢀꢀꢀꢀꢀꢀꢀ
(1)
10.式中，h0为建立粗糙表面模型的初始基准；k为放大系数；z
b1
(xi,yj)为集合b1中点(xi,yj)对应的z坐标值；z
b0
为集合b1中所有点z坐标值的平均值；zb(xi,yj)为集合b中点(xi,yj)对应的z坐标值；xi和yj分别表示每个点的x和y坐标，i＝1,2,
……
,m，j＝1,2,
……
,n，m和n分别表示集合b1在x和y方向原子的个数；
11.步骤1-3：通过原子排布规律，计算得到粗糙表面模型中每个原子的三维坐标；设粗糙表面模型中任一原子的z坐标为z(xi,yj)；将每个原子的z坐标值与集合b中对应点的z坐标值进行对比，如果z(xi,yj)≤zb(xi,yj)，则将该原子的原子序号、原子类型和原子三维坐标信息输出到文件“initialsubstrate.txt”中，反之，如果z(xi,yj)》zb(xi,yj)，则将该原子的原子序号、原子类型和原子三维坐标信息输出到文件“initialindenter.txt”中；
12.步骤1-4：利用最小二乘法，计算文件“initialsubstrate.txt”所对应模型的基准面方程，进而计算得到所建模型表面粗糙度值；若粗糙度值不符合设定的建模要求，则返回步骤1-1重新开始建模；
13.步骤2：建立含有随机粗糙表面和润滑膜的分子动力学模型；
14.编写运行lammps软件所需的in文件，将文件“initialindenter.txt”包含的原子信息读取到in文件中，所对应的分子动力学模型设为压头，将该部分原子上移d的距离；将文件“initialsubstrate.txt”包含的原子信息读取到in文件中，所对应的分子动力学模型设为基体，在压头和基体之间建立润滑膜；所述基体和压头均由固定层、恒温层和牛顿层组成，所述润滑膜设为单独一个组；
15.步骤3：压头下压，使润滑膜附着于基体表面；
16.步骤3-1：设定模拟区域的边界条件、原子间相互作用势函数、原子的初始速度、系统温度和积分步长；设定完成后，模型首先在正则系综nvt下弛豫，之后转为在微正则系综nve下弛豫；
17.步骤3-2：压头下压，下压距离为d1：
18.d1＝d-h
f-c1，d1》0
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
19.式中，d表示压头与基体之间的距离；hf表示润滑膜的厚度；c1是为了确保压头下压过程中不会对基体造成影响而引入的常量；
20.步骤3-3：压头下压完成后，将模型在微正则系综nve下弛豫预设时间长度；
21.步骤4：输出基体和润滑膜的原子信息，建立目标模型；
22.步骤4-1：分别输出基体和润滑膜原子信息，原子信息包括原子序号、原子类型和原子三维坐标；设输出只含基体原子信息的文件为“finalsubstrate.txt”，只含润滑膜原子信息的文件为“finalfilm.txt”，对输出的信息文件中原子序号重新排序，确保原子序号连续，输出文件的结构与“initialsubstrate.txt”相同；
23.步骤4-2：编写运行lammps软件所需的最终in文件，将文件“finalsubstrate.txt”包含的原子信息读取到最终所需的in文件中，所对应的模型设为基体，将文件“finalfilm.txt”包含的原子信息读取到最终所需的in文件中，所对应的模型设为润滑膜，即完成用于分子动力学模拟的润滑膜附着于三维随机粗糙表面的模型建立。
24.优选地，步骤1-3中输出到的文件由记录模型信息的部分和记录各原子信息的部分组成：第一部分，即记录模型信息的部分，第一行为空行，或以“#”开头的内容，第二行记录了模型所含原子总数，第三行记录了模型所含原子类型总数，第四到六行分别记录了模型在x、y、z方向的尺寸；第二部分，即记录各原子信息的部分，第一行为“atoms”，空一行之后的每行分别记录一个原子的原子信息，每行的信息依次为原子序号、原子类型、原子x坐标、原子y坐标和原子z坐标，每个数据之间用空格隔开，第一部分与第二部分之间不需要空行。
25.优选地，所述步骤1-4中利用最小二乘法，计算文件“initialsubstrate.txt”所对应模型的基准面方程，进而计算得到所建模型表面粗糙度值具体为：
26.设集合c为文件“initialsubstrate.txt”中每个(xi,yj)所对应z坐标值最大的点所组成的集合，设基准面方程为：
27.f(x,y)＝a+bc+cy
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
28.利用公式推导公式计算基准面方程中的三个未知数a，b，c，推导公式为：
[0029][0030]
式中zc(xi,yj)表示集合c中点(xi,yj)对应的z坐标值，分别表示集合c中所有点在x，y，z方向坐标值的平均值；
[0031]
计算三维随机粗糙表面的粗糙度参数——表面算术平均偏差sa：
[0032][0033]
优选地，所述c1的取值范围为0.2-0.3nm。
[0034]
本发明的有益效果如下：
[0035]
与传统的分子动力学建模方法相比，本发明更加贴合实际工作条件，提高了模拟结果的准确性。实际工作条件下，相对运动零部件表面不可能是光滑表面、理想粗糙表面或二维随机粗糙表面，但现有模拟方法大多集中于将模型表面简化为上述形式，很少将模型表面建立为三维随机粗糙表面，因此该建模方法更加贴合实际；考虑到实际工作条件下，相对运动的零部件在工作过程中常会采取一定的润滑措施(如在表面覆上一定的固体润滑剂)，在本发明中，结合这些因素，建立的最终模型将润滑膜附着于表面为三维随机粗糙表面的基体上，与实际情况相吻合，为后续的研究奠定了基础。
附图说明
[0036]
图1为本发明流程图。
[0037]
图2为本发明步骤1的流程图。
[0038]
图3为本发明实施例中“initialsubstrate.txt”文件结构。
[0039]
图4为本发明实施例中建立的分子动力学模型。
[0040]
图5为本发明步骤3中模型运动流程图。
[0041]
图6为本发明实施例提取的基体模型。
[0042]
图7为本发明实施例提取的润滑膜模型。
[0043]
图8为本发明实施例最终建立的分子动力学模型。
具体实施方式
[0044]
下面结合附图和实施例对本发明进一步说明。
[0045]
现有基于分子动力学的摩擦学研究主要将基体表面简化为光滑表面，少数模型会将基体表面简化为三维规则粗糙表面或二维随机粗糙表面，但在实际生产中，基体表面应该是三维随机粗糙表面，同时考虑到在基体表面附着润滑膜是减小相对运动零部件摩擦磨损常见方式，而现有分子动力学模型中很少将润滑膜附着于三维随机粗糙表面，因此本发明提出了一种用于分子动力学模拟的生成粗糙表面并覆膜的建模方法，通过此方法可以建立基体表面为三维随机粗糙表面，润滑膜附着于基体表面的分子动力学模型，以便后续的计算研究。
[0046]
一种分子动力学模拟的生成粗糙表面并覆膜的建模方法，包括如下步骤：
[0047]
步骤1：输出表面为三维随机粗糙表面的模型文件；
[0048]
步骤1-1：利用高斯随机表面生成程序生成一组包含坐标信息的随机点，设随机点所组成的集合为a，其中集合a中点的z坐标为复数；
[0049]
步骤1-2：取集合a中每个点z坐标的实部，得到纯实部集合b1，依据如下公式(1)计算得到集合b：
[0050]
zb(xi,yj)＝h0+k(z
b1
(xi,yj)-z
b0
)
ꢀꢀꢀꢀꢀꢀꢀ
(1)
[0051]
式中，h0为建立粗糙表面模型的初始基准；k为放大系数；z
b1
(xi,yj)为集合b1中点(xi,yj)对应的z坐标值；z
b0
为集合b1中所有点z坐标值的平均值；zb(xi,yj)为集合b中点(xi,yj)对应的z坐标值；xi和yj分别表示每个点的x和y坐标，i＝1,2,
……
,m，j＝1,2,
……
,n，m和n分别表示集合b1在x和y方向原子的个数；
[0052]
步骤1-3：通过原子排布规律，计算得到粗糙表面模型中每个原子的三维坐标；设粗糙表面模型中任一原子的z坐标为z(xi,yj)；将每个原子的z坐标值与集合b中对应点的z坐标值进行对比，如果z(xi,yj)≤zb(xi,yj)，则将该原子的原子序号、原子类型和原子三维坐标信息输出到文件“initialsubstrate.txt”中，反之，如果z(xi,yj)》zb(xi,yj)，则将该原子的原子序号、原子类型和原子三维坐标信息输出到文件“initialindenter.txt”中；
[0053]
输出到的文件由记录模型信息的部分和记录各原子信息的部分组成：第一部分，即记录模型信息的部分，第一行为空行，或以“#”开头的内容，第二行记录了模型所含原子总数，第三行记录了模型所含原子类型总数，第四到六行分别记录了模型在x、y、z方向的尺寸；第二部分，即记录各原子信息的部分，第一行为“atoms”，空一行之后的每行分别记录一个原子的原子信息，每行的信息依次为原子序号、原子类型、原子x坐标、原子y坐标和原子z坐标，每个数据之间用空格隔开，第一部分与第二部分之间不需要空行。
[0054]
步骤1-4：利用最小二乘法，计算文件“initialsubstrate.txt”所对应模型的基准面方程，进而计算得到所建模型表面粗糙度值；
[0055]
设集合c为文件“initialsubstrate.txt”中每个(xi,yj)所对应z坐标值最大的点所组成的集合，设基准面方程为：
[0056]
f(x,y)＝a+bx+cy
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0057]
利用公式推导公式计算基准面方程中的三个未知数a，b，c，推导公式为：
[0058][0059]
式中zc(xi,yj)表示集合c中点(xi,yj)对应的z坐标值，分别表示集合c中所有点在x，y，z方向坐标值的平均值；
[0060]
计算三维随机粗糙表面的粗糙度参数——表面算术平均偏差sa：
[0061][0062]
若粗糙度值不符合设定的建模要求，则返回步骤1-1重新开始建模；
[0063]
步骤2：建立含有随机粗糙表面和润滑膜的分子动力学模型；
[0064]
编写运行lammps软件所需的in文件，将文件“initialindenter.txt”包含的原子信息读取到in文件中，所对应的分子动力学模型设为压头，将该部分原子上移d的距离；将文件“initialsubstrate.txt”包含的原子信息读取到in文件中，所对应的分子动力学模型设为基体，在压头和基体之间建立润滑膜；所述基体和压头均由固定层、恒温层和牛顿层组成，所述润滑膜设为单独一个组；
[0065]
步骤3：压头下压，使润滑膜附着于基体表面；
[0066]
步骤3-1：设定模拟区域x向和y向为周期性边界条件，用以减小尺寸效应，通过系统内不同元素原子之间相互作用势函数，获得原子之间的相互作用，根据maxwell-boltzmann能量分布函数设定系统内所有原子的初始速度，该速度对应于其各自平衡状态的温度；设定系统温度为300k，积分步长选择为2fs(1fs＝1e(-15)s)。在完成上述设定之后，在正则系综(nvt)下，依据所设定系统内各原子的初始位置和初始速度，根据选取的势函数计算系统中各原子的势能、位置及速度，对压头和基体原子进行平衡约束，以使系统初始模型达到平衡状态；之后，在微正则系综(nve)下对系统进行内部平衡态演化，使系统达到平衡状态。
[0067]
步骤3-2：压头下压，下压距离为d1：
[0068]
d1＝d-h
f-c1，d1》0
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0069]
式中，d表示压头与基体之间的距离；hf表示润滑膜的厚度；c1是为了确保压头下压过程中不会对基体造成影响而引入的常量；c1的取值范围为0.2-0.3nm；
[0070]
步骤3-3：压头下压完成后，将模型在微正则系综nve下弛豫一段时间；
[0071]
步骤4：输出基体和润滑膜的原子信息，建立目标模型；
[0072]
步骤4-1：分别输出基体和润滑膜原子信息，原子信息包括原子序号、原子类型和
原子三维坐标；设输出只含基体原子信息的文件为“finalsubstrate.txt”，只含润滑膜原子信息的文件为“finalfilm.txt”，对输出的信息文件中原子序号重新排序，确保原子序号连续，输出文件的结构与“initialsubstrate.txt”相同；
[0073]
步骤4-2：编写运行lammps软件所需的最终in文件，将文件“finalsubstrate.txt”包含的原子信息读取到最终所需的in文件中，所对应的模型设为基体，将文件“finalfilm.txt”包含的原子信息读取到最终所需的in文件中，所对应的模型设为润滑膜，即完成用于分子动力学模拟的润滑膜附着于三维随机粗糙表面的模型建立。
[0074]
具体实施例：
[0075]
参见图1所示，本实施例提供了用于分子动力学模拟的生成粗糙表面并覆膜的建模方法，具体步骤如下：
[0076]
步骤一：输出表面为三维随机粗糙表面的模型文件；
[0077]
参见图2所示，该步骤包括以下子步骤：
[0078]
步骤1：利用高斯随机表面生成程序生成一组包含坐标信息的随机点，设随机点所组成的集合为a，其中集合a中点的z坐标为复数。
[0079]
步骤2：取集合a中每个点z坐标的实部，设得到纯实部集合为b1，依据如下公式计算得到集合b。
[0080]
zb(xi,yj)＝h0+k(z
b1
(xi,yj)-z
b0
)
[0081]
式中，h0为建立粗糙表面模型的初始基准，取决于所建粗糙表面模型在z方向的尺寸大小；k为放大系数；z
b1
(xi,yj)为集合b1中点(xi,yj)对应的z坐标值；z
b0
为集合b1中所有点z坐标值的平均值；zb(xi,yj)为集合b中点(xi,yj)对应的z坐标值；xi和yj分别表示每个点的x和y坐标，i＝1,2,
……
,m，j＝1,2,
……
,n，m和n分别表示集合b1在x和y方向原子的个数，取决于所建粗糙表面模型在x和y方向尺寸的大小。
[0082]
步骤3：通过原子排布规律，计算得到完整模型中每个原子的三维坐标，设完整模型中任一原子的z坐标为z(xi,yj)，其中xi和yj分别表示该原子的x和y坐标，i和j的取值与集合b相应符号完全相同，将每个原子的z坐标值与集合b中对应点的z坐标值进行对比，如果z(xi,yj)≤zb(xi,yj)，则将该原子的原子序号、原子类型和原子三维坐标等信息输出到文件“initialsubstrate.txt”中，反之，则将该原子相应信息输出到文件“initialindenter.txt”中，输出文件由记录模型信息的部分和记录各原子信息的部分组成。
[0083]
步骤4：利用最小二乘法，计算文件“initialsubstrate.txt”所对应模型的基准面方程，进而计算得到所建模型表面粗糙度值。设集合c为文件“initialsubstrate.txt”中每个(xi,yj)所对应z坐标值最大的点所组成的集合，设基准面方程为：
[0084]
f(x,y)＝a+bx+cy
[0085]
利用公式推导公式计算基准面方程中的三个未知数a，b，c，推导公式为：
[0086][0087]
式中zc(xi,yj)表示集合c中点(xi,yj)对应的z坐标值，分别表示集合c中所有点在x，y，z方向坐标值的平均值，i、j、m和n的取值与集合b完全相同。
[0088]
将集合c中每个点的坐标信息代入推导公式中，计算得到基准面方程的具体表达式，经计算由文件“initialsubstrate.txt”生成的三维随机粗糙表面的基准面方程为：
[0089]
f(x,y)＝103.1111-0.0194x+0y
[0090]
计算三维随机粗糙表面的粗糙度参数——表面算术平均偏差sa。
[0091][0092]
在本实施例中通过计算得到的粗糙度值为若粗糙度值不符合建模要求，则返回步骤1重新建模。
[0093]
参见图3所示，以文件“initialsubstrate.txt”为例，所述步骤一中步骤3输出文件由记录模型信息的部分和记录各原子信息的部分组成，第一部分，即记录模型信息的部分，第一行为空行，或以“#”开头的内容，第二行记录了模型所含原子总数，第三行记录了模型所含原子类型总数，第四到六行分别记录了模型在x、y、z方向的尺寸；第二部分，即记录各原子信息的部分，第一行为“atoms”，空一行之后的每行分别记录一个原子的原子信息，每行的信息依次为原子序号、原子类型、原子x坐标、原子y坐标和原子z坐标，每个数据之间用空格隔开，第二部分与第一部分之间不需要空行。
[0094]
步骤二建立含有随机粗糙表面和润滑膜的分子动力学模型，该步骤方法如下：
[0095]
编写运行lammps软件所需的in文件，将文件“initialindenter.txt”包含的原子信息读取到in文件中，所对应的分子动力学模型设为压头，将该部分原子上移d的距离；将文件“initialsubstrate.txt”包含的原子信息读取到in文件中，所对应的分子动力学模型设为基体，在压头和基体之间建立润滑膜；如图4所示。对模型进行分组，所述基体和压头均由固定层、恒温层和牛顿层组成，将基体底部和压头顶部一个晶胞范围内的原子设为固定层，将基体固定层上方和压头固定层下方一个晶胞范围内的原子设置为恒温层原子，用以吸收运动过程中产生的热量，将基体和压头其余原子设为牛顿层，其中恒温层和牛顿层原子均遵循牛顿第二定律，将润滑膜为设为单独一个组。在本实施例中，所述基体含161618个原子，压头含97582个原子，润滑膜含9282个原子。
[0096]
参见图5，步骤三压头下压，使润滑膜附着于基体表面，该步骤包括以下子步骤：
[0097]
步骤1：设定模拟区域x向和y向为周期性边界条件，用以减小尺寸效应，通过系统内不同元素原子之间相互作用势函数，获得原子之间的相互作用，在本实施例中共包含三种原子，分别是金(au)、硫(s)和钼(mo)，所述基体和压头均由金原子组成，润滑膜为二硫化钼，由硫原子和钼原子组成，其中基体金原子之间和压头金原子之间均采用eam势，基体金
原子与压头金原子之间采用morse势，润滑膜的硫原子和钼原子之间采用rebo势，压头和基体的金原子与二硫化钼中硫原子和钼原子均采用l-j势。根据maxwell-boltzmann能量分布函数设定系统内所有原子的初始速度，该速度对应于其各自平衡状态的温度；设定系统温度为300k，积分步长选择为2fs。在正则系综(nvt)下，依据所设定系统内各原子的初始位置和初始速度，根据选取的势函数计算系统中各原子的势能、位置及速度，对压头和基体原子进行平衡约束，以使系统初始模型达到平衡状态；完成之后，在微正则系综(nve)下对系统进行内部平衡态演化，使系统达到平衡状态。
[0098]
步骤2：压头下压，下压距离为d1。
[0099]
d1＝d-h
f-c1，d1》0
[0100]
式中，d表示压头与基体之间的距离，hf表示润滑膜的厚度，c1是为了确保压头下压过程中不会对基体造成影响而引入的量，取0.2-0.3nm为宜。
[0101]
步骤3：压头下压完成后，将模型在(nve)系综下弛豫一段时间，在本实施例中，弛豫步长为10000步，即20000fs。
[0102]
步骤四输出基体和润滑膜的原子信息，建立目标模型，该步骤包括以下子步骤：
[0103]
步骤1：分别输出基体和润滑膜原子信息，原子信息包括原子序号，原子类型和原子三维坐标；设输出只含基体原子信息的文件为“finalsubstrate.txt”，只含润滑膜原子信息的文件为“finalfilm.txt”，对输出的信息文件中原子序号重新排序，确保原子序号连的，输出文件的结构与“initialsubstrate.txt”相同，如图6和图7所示分别为文件“finalsubstrate.txt”和文件“finalfilm.txt”所对应的模型。
[0104]
步骤2：编写运行lammps软件所需的最终in文件，将文件“finalsubstrate.txt”包含的原子信息读取到最终所需的in文件中，所对应的模型设为基体，将文件“finalfilm.txt”包含的原子信息读取到最终所需的in文件中，所对应的模型设为润滑膜，即完成用于分子动力学模拟的润滑膜附着于三维随机粗糙表面的模型建立，如图8所示。