一种基于分子动力学模拟微观结构表面防结冰性能的评价方法与流程

本发明属于材料表面性能评价技术领域，涉及一种基于分子动力学模拟微观结构表面防结冰性能的评价方法。本发明所提供的分子动力学模拟评定方法意在对高效防结冰表面的微观结构设计提供理论指导。

背景技术：

结冰是一种常见的自然现象，但是在某些特殊工业条件下，结冰会带来严重的安全危害，例如发动机叶片，机械轴承，传输电线上结冰会导致十分严重的事故，如何避免特殊工况下的结冰现象是当今颇为重要的研究课题。

材料表面微观结构对于冷凝结冰过程有着非常重要的影响，主要体现在两方面：一方面，水在基体表面结冰需要先润湿基体表面，表面特殊微观结构能够将空气“封锁”在结构内部，这些空气能够有效地隔离开水和基体，使得固液接触面积减小，接触角增大，影响水在基体上的润湿过程，从而抑制结冰形核的发生。另一方面，微观结构自身的角度与曲率等几何因素也会影响水在微观结构上的形核结冰过程。较大的曲率，特殊的角度等因素都会抑制水分子的形核结冰，判别并设计有效抑制形核结冰的微观结构对于防结冰技术领域的探索至关重要。

现有实验方法可以通过试错法来构建大量不同的微观结构，逐一检测它们的防结冰性能，从而得到一些规律。但试错法需要大量的时间成本和材料成本，并且无法穷尽各种微观结构形貌。而且实验无法捕捉到结冰形核发生时的相关情况，实验观察到结冰现象时，形核结冰其实早开始。因此，需要寻找新的理论研究方法。

分子动力学方法由于其是从分子层面观察分析问题，更加接近事物本质，易于捕捉事物的微小变化，而且不需要耗费实际的资源，逐渐成为一种重要的研究方法。由于形核是一个随机过程，虽然分子动力学技术可以通过统计学方法计算形核率，但需要耗费大量的时间与计算资源。

技术实现要素：

本发明的目的在于提供一种基于分子动力学模拟微观结构表面防结冰性能的评价方法，以期对防结冰表面微观结构的设计提供理论指导。

技术方案：

一种基于分子动力学模拟微观结构表面防结冰性能的评价方法，包括以下步骤：

1)基于分子动力学仿真软件lammps建立基体表面不同形状微结构与水的接触模型；采用三维周期盒子，首先在xy轴上建立无限长、厚度一定的平面基体，然后在平面基体上构建规则阵列微结构，控制微结构上下面的形状与尺寸得到特定形状的微结构，同时定义微结构高度与间距；微结构创建后，固定基体与微结构的原子坐标；在微结构上生成水分子，保证水分子与基体距离保持在截断半径内，并进行能量最小化设置；

2)基于分子动力学仿真软件lammps对步骤1)中的模型进行初始化准备；具体步骤为：首先赋予原子质量与粒子初始速度，采用sw粗粒化势函数表征水与基体原子之间的相互作用，将系统处于室温，一个大气压下进行数纳秒平衡弛豫，系综选择为npt系综，积分步长可自行选定，使系统达到平衡，密度合理；

3)基于分子动力学仿真软件lammps对步骤2)得到的平衡态模型进行降温动力学弛豫；在一个大气压下，设定温度从298k均匀变化到200k或以下温度，降温速度可因精确度而定，每隔固定步长设定一个状态点，利用dump命令输出各个状态点的水分子坐标；

4)基于分子动力学仿真软件lammps对步骤2)得到的平衡态模型进行低温平衡态动力学弛豫；在一个大气压下，设定保温温度，进行平衡态动力学弛豫；每隔固定步长设定一个状态点，利用dump命令输出各个状态点的水分子坐标；

5)利用脚本程序读取降温动力学过程得出的数据结果，判定各个分子的结冰情况，并输出每个状态点的结冰分子数量；得到数据后利用origin软件作图，得到结冰粒子数与温度变化之间的关系图，取结冰分子数量开始明显上升的温度作为结冰温度；

6)读取低温平衡态动力学弛豫得出的数据结果，利用q6判据判定各个分子的结冰情况，并输出每个状态点的结冰分子数量；得到数据后利用origin软件作图，得到结冰分子数与时间变化之间的关系图，取结冰分子数量开始明显上升的温度作为结冰时间；

7)利用步骤5)、6)得到的结冰温度与结冰时间数据，评价微结构形状抑制结冰性能的优劣。

所述的评价方法，在步骤1)中，建立不同微结构与水接触的模型：微结构为周期性阵列结构，控制微观结构形状与尺寸，微观结构高度和单个微观结构之间的距离根据需求定义，高度定为1nm-4nm，间距定为1nm-4nm；微结构的形状是圆台，立方体，棱台之一。

所述的评价方法，在步骤2)中，参数初始化设置，基体元素种类可以改变，sw势函数适用于描述粗粒化水分子与基体的相互作用，具体表现形式为：

式中ε为能量参数，表征原子之间的相互作用强度；σ为距离参数，表征截断半径；θ为水中三个原子所成的角度；r表示两原子之间的距离；其余参数为无量纲的常数。其中φ2描述二体势能量，φ3描述三体势能量，系统总能量e等于二体势能量与三体势能量的加和，对于微结构与水接触的异质形核问题，水分子之间的相互作用由三体势描述，水与基体原子间的相互作用由二体势描述，由于基体原子固定，则不用设定基体原子间的作用势，以减少计算量；由于基体原子不同，水分子与基体原子的势函数值也会有所不同，具体可通过赋予不同的ε值改变基体的元素种类。

所述的评价方法，在步骤3)中，降温速度，降温范围与输出时间步长均可自行设定。降温速度一般在0.1k/ns-1k/ns之间，降温范围从室温到200k附近任意值。

所述的评价方法，在步骤4)中，设定的平衡温度可以设定为冰点以下的任意温度，优选选择230k，保温时间需要随保温温度改变，温度较高时，保温时间需要延长，温度较低时保温时间缩短。

所述的评价方法，在步骤5)中，结冰温度为结冰分子数目结束震荡，开始明显上升的温度。

所述的评价方法，在步骤6)中，结冰时间为结冰分子数目结束震荡，开始明显上升的时间。

所述的评价方法，在步骤7)中，结冰温度越低，结冰时间越长，微结构对于结冰的抑制能力越优。

采用本发明的材料微观结构表面防结冰性能的评价方法具有以下优势：

排除了实际实验时外界不稳定因素的影响，可直接判断微观结构几何参数这一自变量对于结冰难易程度的影响，避免了试错法的资源浪费，大大降低了实验成本。

规避了计算量需求巨大的形核率计算，选择了同样能表征结冰难易程度的结冰温度，结冰时间参数来表征微观结构表面防结冰性能，大幅度减小了计算量。

给出了较为合理的模拟参数范围，既可以保证得到结冰结果，又可以尽量减小计算量，提高了运算效率。

附图说明

图1为基体微观结构模型的示意图；

图2为基体微观结构与水接触的模型示意图；

图3为结冰粒子数与温度的变化关系图；

图4为结冰粒子数与时间的变化关系图；

图5为本发明方法流程图；

具体实施方式

以下结合具体实施例，对本发明进行详细说明。

实施例1：

本发明的一种基于分子动力学模拟微观结构表面防结冰性能的评价方法如下：

1)基于分子动力学仿真软件lammps建立碳基体上圆柱状微结构与水的接触模型。采用三维周期盒子，首先建立在xy轴上无限长，厚度为6个原子层的平面基体，后在基体上构建规则阵列微结构。固定微结构上下表面为圆形，上表面半径为1.6nm，下表面半径为1.6nm，微结构高度为1.2nm，使微结构呈圆柱状，间距为1.5nm。微结构创建后，固定基体与微结构的原子坐标。在微结构上生成1000个水分子，

2)基于分子动力学仿真软件lammps对步骤1)中的模型进行初始化准备。具体步骤为：采用sw粗粒化势函数表征水与碳基体原子之间的相互作用，根据模拟体系，确定势函数参数。其中，水分子之间参数为ε＝6.189kcal/mol，水和碳基体之间ε＝0.13kcal/mol，其余参数如下：a＝1.80,λ＝23.15,γ＝1.20，cosθ＝-0.333333，a＝7.04955，b＝0.60222，p＝4.0,q＝0.0。将系统处于298k，一个大气压下进行平衡弛豫，系综选择为npt，积分步长为3fs，使系统达到平衡，密度合理。

3)基于分子动力学仿真软件lammps对步骤2)得到的平衡模型进行降温动力学弛豫。在一个大气压下，设定温度从298k变化到200k，降温速度为1k/ns，使系统温度从298k室温匀速降到200k，利用dump命令输出各个状态点的水分子坐标，输出步长为1000步。

4)基于分子动力学仿真软件lammps对步骤2)得到的平衡态模型进行低温平衡态动力学弛豫。在一个大气压下，设定温度为230k，进行100ns平衡态动力学弛豫。利用dump命令输出各个状态点的水分子坐标，输出步长为1000步。

5)对于降温动力学过程得出的数据结果，利用q6判据(q6空间矢量判定脚本程序)，判定各个分子的结冰情况，并输出每个状态点的结冰分子数量。得到数据后利用origin软件作图，得到结冰粒子数与温度变化之间的关系图，取结冰分子数量开始明显上升的温度作为结冰温度，结冰温度为221k。

6)对于低温平衡态动力学得出的数据结果，利用q6判据判据判定各个分子的结冰情况，并输出每个状态点的结冰分子数量。得到数据后利用origin软件作图，得到结冰粒子数与时间变化之间的关系图，取结冰分子数量开始明显上升的温度作为结冰时间，结冰时间为89ns。(q6判据参考文献：tianshuli,davidedonadio,giovannarusso,giuliagalli，homogeneousicenucleationfromsupercooledwater[j]phys.chem.chem.phys2011,13,19807–19813)

7)本例中圆柱状微结构上水结冰温度为221k，结冰时间为89ns，相比于在平面基体上的结冰过程，结冰温度较高，结冰时间缩短，说明本例中圆柱状微结构比平面基体更有利于结冰形核。

实施例2：

本发明的一种基于分子动力学模拟微观结构表面防结冰性能的评价方法如下：

1)基于分子动力学仿真软件lammps建立碳基体上圆台状微结构与水的接触模型。采用三维周期盒子，首先建立在xy轴上无限长，厚度为6个原子层的平面基体，后在基体上构建规则阵列微结构。固定微结构上下表面为圆形，上表面半径为0.4nm，下表面半径为1.6nm，微结构高度为1.2nm，使微结构呈圆台状，间距为1.5nm。微结构创建后，固定基体与微结构的原子坐标。在微结构上生成1000个水分子，

2)基于分子动力学仿真软件lammps对步骤1)中的模型进行初始化准备。具体步骤为：。采用sw粗粒化势函数表征水与碳基体原子之间的相互作用，根据模拟体系，确定势函数参数。其中，水分子之间参数为ε＝6.189kcal/mol，水和碳基体之间ε＝0.13kcal/mol，将系统处于298k，一个大气压下进行平衡弛豫，系综选择为npt，积分步长为3fs，使系统达到平衡，密度合理。

3)基于分子动力学仿真软件lammps对步骤2)得到的平衡模型进行降温动力学弛豫。在一个大气压下，设定温度从298k变化到200k，降温速度为1k/ns，使系统温度从298k室温匀速降到200k，利用dump命令输出各个状态点的水分子坐标，输出步长为1000步。

4)基于分子动力学仿真软件lammps对步骤2)得到的平衡态模型进行低温平衡态动力学弛豫。在一个大气压下，设定温度为230k，进行120ns平衡态动力学弛豫。利用dump命令输出各个状态点的水分子坐标，输出步长为1000步。

5)对于降温动力学过程得出的数据结果，利用q6判据，判定各个分子的结冰情况，并输出每个状态点的结冰分子数量。得到数据后利用origin软件作图，得到结冰粒子数与温度变化之间的关系图，取结冰分子数量开始明显上升的温度作为结冰温度，结冰温度为213k。

6)对于低温平衡态动力学得出的数据结果，利用q6判据，判定各个分子的结冰情况，并输出每个状态点的结冰分子数量。得到数据后利用origin软件作图，得到结冰粒子数与时间变化之间的关系图，取结冰分子数量开始明显上升的温度作为结冰时间，结冰时间为96ns。

7)本例中圆台状微结构上水结冰温度为213k，结冰时间为96ns，相比于在平面基体上的结冰过程，结冰温度和结冰时间均无明显区别，说明本例中圆台状微结构没有明显的抑制结冰能力。

实施例3：

本发明的一种基于分子动力学模拟微观结构表面防结冰性能的评价方法如下：

1)基于分子动力学仿真软件lammps建立碳基体上圆柱状微结构形状与水的接触模型。采用三维周期盒子，首先建立在xy轴上无限长，厚度为6个原子层的平面基体，后在基体上构建规则阵列微结构。固定微结构呈圆锥状，，下表面半径为1.6nm，微结构高度为1.2nm，微结构间距为1.5nm。微结构创建后，固定基体与微结构的原子坐标。在微结构上生成1000个水分子，

2)基于分子动力学仿真软件lammps对步骤1)中的模型进行初始化准备。具体步骤为：。采用sw粗粒化势函数表征水与碳基体原子之间的相互作用，根据模拟体系，确定势函数参数。其中，水分子之间参数为ε＝6.189kcal/mol，水和碳基体之间ε＝0.13kcal/mol，将系统处于298k，一个大气压下进行平衡弛豫，系综选择为npt，积分步长为3fs，使系统达到平衡，密度合理。

3)基于分子动力学仿真软件lammps对步骤2)得到的平衡模型进行降温动力学弛豫。在一个大气压下，设定温度从298k变化到200k，降温速度为1k/ns，使系统温度从298k室温匀速降到200k，利用dump命令输出各个状态点的水分子坐标，输出步长为1000步。

4)基于分子动力学仿真软件lammps对步骤2)得到的平衡态模型进行低温平衡态动力学弛豫。在一个大气压下，设定温度为230k，进行150ns平衡态动力学弛豫。利用dump命令输出各个状态点的水分子坐标，输出步长为1000步。

5)对于降温动力学过程得出的数据结果，利用q6判据，判定各个分子的结冰情况，并输出每个状态点的结冰分子数量。得到数据后利用origin软件作图，得到结冰粒子数与温度变化之间的关系图，取结冰分子数量开始明显上升的温度作为结冰温度，结冰温度为204k。

6)对于低温平衡态动力学得出的数据结果，利用q6判据，判定各个分子的结冰情况，并输出每个状态点的结冰分子数量。得到数据后利用origin软件作图，得到结冰粒子数与时间变化之间的关系图，取结冰分子数量开始明显上升的温度作为结冰时间，结冰时间为137ns。

7)本例中圆锥状微结构上结冰温度为204k，结冰时间为137ns，相比于在平面基体上的结冰过程，结冰温度降低，结冰时间延长，说明本例中圆锥状微结构具有明显的抑制结冰能力。

实施例4：

本发明的一种基于分子动力学模拟微观结构表面防结冰性能的评价方法如下：

1)基于分子动力学仿真软件lammps建立碳基体上圆柱状微结构形状与水的接触模型。采用三维周期盒子，首先建立在xy轴上无限长，厚度为6个原子层的平面基体，后在基体上构建规则阵列微结构。固定微结构呈圆锥状，，下表面半径为1.6nm，微结构高度为0.6nm，微结构间距为1.5nm。微结构创建后，固定基体与微结构的原子坐标。在微结构上生成1000个水分子，

2)基于分子动力学仿真软件lammps对步骤1)中的模型进行初始化准备。具体步骤为：。采用sw粗粒化势函数表征水与碳基体原子之间的相互作用，根据模拟体系，确定势函数参数。其中，水分子之间参数为ε＝6.189kcal/mol，水和碳基体之间ε＝0.13kcal/mol，将系统处于298k，一个大气压下进行平衡弛豫，系综选择为npt，积分步长为3fs，使系统达到平衡，密度合理。

3)基于分子动力学仿真软件lammps对步骤2)得到的平衡模型进行降温动力学弛豫。在一个大气压下，设定温度从298k变化到200k，降温速度为1k/ns，使系统温度从298k室温匀速降到200k，利用dump命令输出各个状态点的水分子坐标，输出步长为1000步。

4)基于分子动力学仿真软件lammps对步骤2)得到的平衡态模型进行低温平衡态动力学弛豫。在一个大气压下，设定温度为230k，进行150ns平衡态动力学弛豫。利用dump命令输出各个状态点的水分子坐标，输出步长为1000步。

5)对于降温动力学过程得出的数据结果，利用q6判据，判定各个分子的结冰情况，并输出每个状态点的结冰分子数量。得到数据后利用origin软件作图，得到结冰粒子数与温度变化之间的关系图，取结冰分子数量开始明显上升的温度作为结冰温度，结冰温度为210k。

6)对于低温平衡态动力学得出的数据结果，利用q6判据，判定各个分子的结冰情况，并输出每个状态点的结冰分子数量。得到数据后利用origin软件作图，得到结冰粒子数与时间变化之间的关系图，取结冰分子数量开始明显上升的温度作为结冰时间，结冰时间为101ns。

7)本例中圆锥状微结构上结冰温度为210k，结冰时间为101ns，相比于在平面基体上的结冰过程，结冰温度略有降低，结冰时间略微延长，说明本例中圆锥状微结构具有一定的抑制结冰能力。

实施例5：

本发明的基体微结构形状影响防结冰性能的研究方法如下：

1)基于分子动力学仿真软件lammps建立碳基体上圆柱状微结构形状与水的接触模型。采用三维周期盒子，首先建立在xy轴上无限长，厚度为6个原子层的平面基体，后在基体上构建规则阵列微结构。固定微结构呈圆锥状，，下表面半径为1.6nm，微结构高度为0.6nm，微结构间距为2.5nm。微结构创建后，固定基体与微结构的原子坐标。在微结构上生成1000个水分子，

2)基于分子动力学仿真软件lammps对步骤1)中的模型进行初始化准备。具体步骤为：。采用sw粗粒化势函数表征水与碳基体原子之间的相互作用，根据模拟体系，确定势函数参数。其中，水分子之间参数为ε＝6.189kcal/mol，水和碳基体之间ε＝0.13kcal/mol，将系统处于298k，一个大气压下进行平衡弛豫，系综选择为npt，积分步长为3fs，使系统达到平衡，密度合理。

3)基于分子动力学仿真软件lammps对步骤2)得到的平衡模型进行降温动力学弛豫。在一个大气压下，设定温度从298k变化到200k，降温速度为1k/ns，使系统温度从298k室温匀速降到200k，利用dump命令输出各个状态点的水分子坐标，输出步长为1000步。

4)基于分子动力学仿真软件lammps对步骤2)得到的平衡态模型进行低温平衡态动力学弛豫。在一个大气压下，设定温度为230k，进行150ns平衡态动力学弛豫。利用dump命令输出各个状态点的水分子坐标，输出步长为1000步。

5)对于降温动力学过程得出的数据结果，利用q6判据，判定各个分子的结冰情况，并输出每个状态点的结冰分子数量。得到数据后利用origin软件作图，得到结冰粒子数与温度变化之间的关系图，取结冰分子数量开始明显上升的温度作为结冰温度，结冰温度为213k。

6)对于低温平衡态动力学得出的数据结果，利用q6判据，判定各个分子的结冰情况，并输出每个状态点的结冰分子数量。得到数据后利用origin软件作图，得到结冰粒子数与时间变化之间的关系图，取结冰分子数量开始明显上升的温度作为结冰时间，结冰时间为95ns。

7)本例中圆锥状微结构上结冰温度为213k，结冰时间为95ns，相比于在平面基体上的结冰过程，结冰温度与结冰时间无明显差异，说明本例中圆锥状微结构不具有的抑制结冰能力。

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。