一种基于分子动力学模拟的Stribeck曲线计算方法与流程

一种基于分子动力学模拟的stribeck曲线计算方法
技术领域
1.本发明涉及摩擦学计算技术领域，尤其涉及一种基于分子动力学模拟的stribeck曲线计算方法。


背景技术：

2.早在1902年，德国学者stribeck r.在对轴承进行摩擦学实验研究中，探讨了润滑剂的粘度、相对速度以及法向载荷等参数与摩擦系数的关系，并得到了以“(润滑剂的粘度
×
相对速度)/法向载荷”物理量为横坐标，摩擦系数为纵坐标的函数图像，这就是摩擦领域中著名的stribeck曲线。
3.长期以来，stribeck曲线一直是研究人员用于定性分析摩擦系统性能的工具，并利用它将摩擦系统的润滑状态划分为动压润滑、混合润滑、薄膜润滑以及边界润滑等。同时，该曲线表明，摩擦系统随着法向载荷增大、相对速度和润滑剂粘度的降低，摩擦系数将会增大。完整的stribeck曲线横坐标采用对数坐标形式，其量级跨度较大，一般为6 ~ 14个数量级，因此在理论上难以通过单一数学模型完整且精准地表征出从零油膜到流体动压润滑等多种润滑状态；此外已有的实验设备难以满足如此大的测量工况范围。目前仅有四川大学的dong zhu团队利用准系统弹流润滑模型计算出较为完整的点接触stribeck曲线，并在宽转速摩擦试验台上得到良好的验证。
4.分子动力学是一种在微观尺度上采用牛顿力学理论开展系统性质研究的模拟方法，适用于金属、陶瓷、复合材料、有机物等材料体系。目前，该方法已在微观摩擦领域具有较为成熟的应用，可预测法向载荷、摩擦力、摩擦系数、温度及应力分布等摩擦学特征参量。因此，采用分子动力学方法开展摩擦副的stribeck曲线模拟计算是一种简单有效的方法。


技术实现要素：

5.本发明所要解决的技术问题是提供一种简单、快速的基于分子动力学模拟的stribeck曲线计算方法。
6.为解决上述问题，本发明所述的一种基于分子动力学模拟的stribeck曲线计算方法，包括以下步骤：
⑴
在分子动力学模拟软件中创建包括固定件、润滑剂以及运动件的摩擦系统微观结构的几何模型；
⑵
利用数值计算方法完成几何模型的能量最小化；
⑶
设定模拟系统的计算参数；
⑷
对运动件和润滑剂分别施加法向载荷p1和法向载荷p2，然后在分子动力学模拟软件中给定计算步数，完成加载过程的模拟；
⑸
给定运动件沿运动方向的运动速度v，在分子动力学模拟软件中完成摩擦过程的模拟，并计算运动件的摩擦力fy；
⑹
计算stribeck曲线中涉及的物理量；
⑺
绘制stribeck曲线图。
7.所述步骤
⑴
中固定件为填充有待被摩擦材料原子的长方体盒子ⅰ，并从长方体盒子ⅰ底面起自下而上依次定义两个l型区域，即固定层与恒温层，其余区域定义为牛顿层。
8.所述步骤
⑴
中润滑剂为构成润滑剂的基础分子通过给定润滑剂密度ρ和基础分子数量比例构建出填充润滑油分子的长方体盒子ⅱ，该长方体盒子ⅱ与所述固定件的长方体盒子ⅰ的长、宽大小一致，且位于所述固定件的上表面。
9.所述步骤
⑴
中运动件为填充有摩擦材料原子，半径为r的球体，该球体位于所述润滑剂的长方体盒子ⅱ上方，且二者间距为10 ~ 15
ꢀå
。
10.所述步骤
⑵
中数值计算方法为共轭梯度法、最速下降法、牛顿-拉森法中的一种。
11.所述步骤
⑵
中几何模型的能量最小化按下述方法完成：在分子动力学仿真软件中给定几何模型在迭代过程中的最大迭代步数与迭代收敛精度，然后采用数值计算方法得到稳定的摩擦副微观几何模型。
12.所述最大迭代步数是指模型受力的最大迭代步数与迭代收敛精度和模型能量的最大迭代步数，且两者最大迭代步数为1000 ~ 5000步。
13.所述迭代收敛精度为1
×
10-10 ~ 1
×
10-15
。
14.所述步骤
⑶
中模拟系统的计算参数包括模拟系统的长方体盒子ⅲ三个正交方向的边界条件、模拟时间步长、系统初始温度t0、系综和原子之间的势函数。
15.所述边界条件为周期性边界条件、自由边界条件或固定边界条件中的一种。
16.所述模拟时间步长为被研究材料的分子运动最小振动周期的1/30 ~ 1/10。
17.所述系统初始温度t0为给定系统内全部原子在0时刻的温度，其中恒温层采用系统控温方法进行控温。
18.所述系综为nvt系综、nve系统、npt系综或nvh系综中的一种或多种。
19.所述原子之间的势函数是指设定系统为混合力场，并给出系统内任意两种原子的力场；所述力场为l-j力场、l-j库伦力场、l-j长程库伦力场、l-j短程库伦力场、mores力场、eam力场、meam力场、opls力场、pcff力场、cvff力场、tersoff力场、reax/c力场中的一种或多种。
20.所述步骤
⑷
中法向载荷p1为固定件与运动件的相互作用力，且p1为0.1 ~ 10 kcal/(mol
·
å
)。
21.所述步骤
⑷
中法向载荷p2为润滑剂所处环境的边界压力，且边界压力为标准大气压。
22.所述步骤
⑷
中计算步数为能够使系统弛豫后达到动态平衡所需要的计算步数，且计算步数为1
×
10
4 ~ 1
×
107步。
23.所述步骤
⑸
中运动速度v为运动件与固定件间的相对运动速度，且运动速度为0.0001 ~ 1000 m/s。
24.所述步骤
⑹
中涉及的物理量分别为润滑剂温度t、润滑剂粘度η、系统的摩擦系数μ和stribeck曲线横坐标参量u。
25.本发明与现有技术相比具有以下优点：1、本发明通过建立固定件、润滑剂以及运动件的摩擦系统几何模型，开展在给定温度、载荷和相对运动速度工况下，系统在摩擦过程中的分子动力学模拟，得到摩擦系数，
进而得到摩擦系统的stribeck曲线。
26.2、本发明在理论上可以得到外界工况与系统的摩擦系数的定量关系，解决了难以通过同一种计算模型完整地描述出从零油膜到流体动压润滑多种润滑状态的stribeck曲线计算问题，可为材料摩擦、润滑等领域的润滑状态定量划分研究提供依据。
27.3、本发明方法简单、快速、易于实现，适用于点、线、面等各种接触形式的摩擦系统。
附图说明
28.下面结合附图对本发明的具体实施方式作进一步详细的说明。
29.图1为本发明聚α烯烃基础油的等规分子模型。
30.图2为本发明聚α烯烃基础油的间规分子模型。
31.图3为本发明聚α烯烃基础油构成的润滑剂模型。
32.图4为本发明实施例建立初始状态的微观几何模型。
33.图5为本发明实施例加载后达到平衡的微观几何模型。
34.图6为本发明实施例通过分子动力学计算方法得到的stribeck曲线。
具体实施方式
35.一种基于分子动力学模拟的stribeck曲线计算方法，包括以下步骤：
⑴
在分子动力学模拟软件中创建包括固定件、润滑剂以及运动件的摩擦系统微观结构的几何模型。
36.其中：固定件为填充有待被摩擦材料原子的长方体盒子ⅰ，并从长方体盒子ⅰ底面起自下而上依次定义两个l型区域，即固定层与恒温层，其余区域定义为牛顿层。
37.润滑剂为构成润滑剂的基础分子通过给定润滑剂密度ρ和基础分子数量比例构建出填充润滑油分子的长方体盒子ⅱ，该长方体盒子ⅱ与固定件的长方体盒子ⅰ的长、宽大小一致，且位于所述固定件的上表面。
38.运动件为填充有摩擦材料原子，半径为r的球体，该球体位于润滑剂的长方体盒子ⅱ上方，且二者间距为10 ~ 15
ꢀå
。
39.⑵
利用数值计算方法完成几何模型的能量最小化。具体过程如下：在分子动力学仿真软件中给定几何模型在迭代过程中的最大迭代步数与迭代收敛精度，最大迭代步数是指模型受力的最大迭代步数与迭代收敛精度和模型能量的最大迭代步数，且两者最大迭代步数为1000 ~ 5000步；迭代收敛精度为1
×
10-10 ~ 1
×
10-15
。然后采用数值计算方法消除因建模不合理或局部应力过大导致的结构不稳定，得到稳定的摩擦副微观几何模型。
40.数值计算方法为共轭梯度法、最速下降法、牛顿-拉森法中的一种。
41.⑶
设定模拟系统的计算参数，计算参数包括模拟系统的长方体盒子ⅲ三个正交方向的边界条件、模拟时间步长、系统初始温度t0、系综和原子之间的势函数。
42.其中：边界条件为周期性边界条件、自由边界条件或固定边界条件中的一种；优选地，边界条件为周期性边界条件。
43.模拟时间步长为被研究材料的分子运动最小振动周期的1/30 ~ 1/10；优选地，时
间步长为0.1 ~ 1000 fs。
44.系统初始温度t0为给定系统内全部原子在0时刻的温度，其中恒温层采用系统控温方法进行控温；优选地，温度为100 ~800 k。
45.系综为nvt系综、nve系统、npt系综或nvh系综中的一种或多种。优选地，固定件与润滑剂的系综为nvt系综，运动件的系综为nve系综。
46.原子之间的势函数是指设定系统为混合力场，并给出系统内任意两种原子的力场；所述力场为l-j力场、l-j库伦力场、l-j长程库伦力场、l-j短程库伦力场、mores力场、eam力场、meam力场、opls力场、pcff力场、cvff力场、tersoff力场、reax/c力场中的一种或多种。优选地，金属材料内部的原子力场为eam力场或meam力场，有机物内部的原子选择pcff力场或cvff力场，金属材料之间的原子力场或金属材料与有机物之间的原子力场为l-j力场。
47.⑷
对运动件和润滑剂分别施加法向载荷p1和法向载荷p2，然后在分子动力学模拟软件中给定计算步数，完成加载过程的模拟。
48.其中：法向载荷p1为固定件与运动件的相互作用力，且p1为0.1 ~ 10 kcal/( mol
·
å
)。
49.法向载荷p2为润滑剂所处环境的边界压力，且边界压力为标准大气压。
50.计算步数为能够使系统弛豫后达到动态平衡所需要的计算步数，且计算步数为1
×
10
4 ~ 1
×
107步。
51.⑸
给定运动件沿运动方向的运动速度v，运动速度v为运动件与固定件间的相对运动速度，且运动速度为0.0001 ~ 1000 m/s。在分子动力学模拟软件中完成摩擦过程的模拟，并计算运动件的摩擦力fy。
52.⑹
计算stribeck曲线中涉及的物理量，所涉及的物理量分别为润滑剂温度t、润滑剂粘度η、系统的摩擦系数μ和stribeck曲线横坐标参量u。
53.系统的摩擦系数μ通过公式（1）计算求取。
54.ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（1）式中，fy为运动件的摩擦力，单位为kcal/(mol
·
å
)；fz为运动件的法向受力，单位为kcal/(mol
·
å
)。
55.润滑剂粘度η由公式（2）与公式（3）计算求取。
56.ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（2）
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（3）式中，a和b为润滑剂的材料常数，由润滑剂生厂商给出或通过粘度计测量，其中a的单位为cst，b的单位为cst/k；ν为润滑剂的运动粘度，单位为cst。
57.润滑剂的温度t可由公式（4）计算求取。
58.ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（4）式中，ek为润滑剂的分子动能，单位为kcal/mol。kb为玻尔兹曼常数，单位为kcal/
( mol
·
k)。
59.stribeck曲线横坐标参量u可由公式（5）计算得到。
60.ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（5）
⑺
绘制stribeck曲线图。
61.通过改变步骤
⑶
中初始温度t或步骤
⑷
中施加在运动件上的法向载荷p1或步骤
⑸
中给定运动件的运动速度v，得到不同的u值与摩擦系数μ。以参量u为横坐标，以系统的摩擦系数μ为纵坐标绘制图像，即得stribeck曲线。
62.实施例 以聚α烯烃基础油为例，一种基于分子动力学模拟的stribeck曲线计算方法，包括以下步骤：
⑴
建立摩擦系统的分子动力学几何模型：首先，建立润滑剂的微观几何模型。在materials studio软件中绘制聚α烯烃基础油的等规分子模型和间规分子模型，分别如图1和图2所示。采用amorphous cell模块，设定聚α烯烃基础油的密度为0.85 g/cm3，等规分子与间规分子的数量比为1:1，填充在边长分别为90
ꢀå
、140
ꢀå
、35
ꢀå
的长方体盒子ⅱ。使用msi2lmp命令，将润滑剂几何模型转化成data文件。在lammps中使用read_data命令读取data文件，并利用change_box命令扩展润滑剂长方体盒子ⅱ高度方向的边长自35
ꢀå
扩展到195
ꢀå
。将润滑剂置于x方向0~90
ꢀå
、y方向0~140
ꢀå
、z方向100~135
ꢀå
的范围内，并建立如图3所示的坐标系。
63.其次，建立固定件的微观几何模型。在扩展后的长方体盒子ⅰ内，在x方向0~90
ꢀå
、y方向0~140
ꢀå
、z方向0~100
ꢀå
范围内的区域填充fe原子，其中定义在x方向0~90
ꢀå
、y方向0~20
ꢀå
、z方向0~100
ꢀå
的区域内以及在x方向0~90
ꢀå
、y方向20~140
ꢀå
、z方向0~20
ꢀå
的区域内为固定层；定义在x方向0~90
ꢀå
、y方向20~40
ꢀå
、z方向20~100
ꢀå
的区域内以及在x方向0~90
ꢀå
、y方向40~140
ꢀå
、z方向20~40
ꢀå
的区域内为恒温层；定义在x方向0~90
ꢀå
、y方向0~20
ꢀå
、z方向0~100
ꢀå
的区域内为牛顿层。
64.最后，建立摩擦系统的分子动力学几何模型，如图4所示。在扩展后的长方体盒子ⅲ内，选择坐标点（45
ꢀå
，95
ꢀå
，165
ꢀå
）为圆心，15
ꢀå
为半径建立一个球体，并在该区域内填充cr原子。
65.⑵
几何模型的能量最小化：在lammps中设定系统能量与受力的最大迭代步数与迭代收敛精度分别为5000和1
×
10-15
，并利用min_stlye命令选择共轭梯度法对摩擦系统的微观几何模型开展能量最小化处理，得到优化的几何模型。
66.⑶
设定模拟系统的计算参数：在lammps中利用boundary命令设定系统的x、y、z三个方向均为周期性边界。
67.利用timestep命令设定模拟步长为1 fs。
68.利用velocity命令设定系统的初始温度t0为300 k。
69.利用fix命令设定系统的系综为nvt，并利用fix_modify命令控制恒温层的温度为300 k。
70.利用fix rigid/nve命令设定运动件为刚体。
71.利用pair_style命令将模拟系统的势函数设定为hybrid类型。其中，润滑剂采用lj/cut/coul/cut势函数、固定件采用eam/fs势函数，运动件采用meam/c势函数，润滑剂与固定件之间、润滑剂与运动件之间、固定件与运动件之间采用lj/cut势函数。
72.⑷
加载过程的模拟：首先，施加载荷。在lammps中利用addforce命令对固定件系统施加大小为0.337 kcal/(mol
·
å
)，方向沿z轴负方向的载荷p1；对润滑剂系统施加大小为0.137 kcal/(mol
·
å
)，方向沿z轴负方向的载荷p2。
73.其次，系统弛豫。设定计算步数为50000步，采用nvt系综进行弛豫。弛豫后沿长方体盒子ⅲ的x方向中面的剖视图如图5所示。
74.⑸
摩擦过程的模拟：首先，赋予速度。在lammps中利用fix命令给定运动件沿y轴负方向的运动速度v为0.001 m/s。
75.其次，系统弛豫。设定计算步数为100000步，与步骤
⑶
的系综类型一致进行弛豫。
76.⑹
计算stribeck曲线中涉及的物理量：在lammps中，利用compute命令计算运动件的受力大小，并利用variable命令定义运动件y方向与z方向的受力为fy、fz。
77.根据公式（1）利用variable命令定义系统的摩擦系数μ并计算在摩擦过程中的数值。
78.根据公式（4）利用variable命令定义润滑剂的温度t并计算在摩擦过程中的数值。
79.⑺
绘制stribeck曲线：通过改变步骤
⑶
中初始温度t0或步骤
⑷
中施加在运动件上的法向载荷p1或步骤
⑸
中给定运动件的运动速度v，并根据公式（2）、（3）、（5）得到不同的stribeck曲线横坐标参量u以及纵坐标参量摩擦系数μ。以参量u为横坐标，摩擦系数μ为纵坐标，绘制曲线。图6为改变步骤
⑸
中的运动速度v的stribeck曲线。
80.应该理解，这里讨论的实施方式只是为了说明，对熟悉该领域的人可以提出各种改进和变化，这些改进和变化将包括在本技术的精神实质和范围以及所附的权利要求范围内。