CN109885895A - 一种基于分子动力学的材料表面结冰形核过程的监测方法 - Google Patents
一种基于分子动力学的材料表面结冰形核过程的监测方法 Download PDFInfo
- Publication number
- CN109885895A CN109885895A CN201910074076.1A CN201910074076A CN109885895A CN 109885895 A CN109885895 A CN 109885895A CN 201910074076 A CN201910074076 A CN 201910074076A CN 109885895 A CN109885895 A CN 109885895A
- Authority
- CN
- China
- Prior art keywords
- icing
- molecular
- hydrone
- molecule
- forming core
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于分子动力学的材料表面结冰形核过程的监测方法，该方法借助分子动力学理论和仿真计算软件，对水在不同物理化学状态的材料表面上异质形核结冰过程进行模拟仿真，克服时间与空间尺度的约束，对水初始形核位置、形核结冰过程中的晶型变化进行监测，基于模拟结果分析不同材料表面结冰形核的过程，探索表面状态与异质形核结冰过程之间的关系，从而有利于从分子层面解释结冰机理，可为各种期望利用表面物理化学改性达到防结冰效果的实验探索、工艺优化与产品开发提供理论支持。
Description
技术领域
本发明涉及一种材料表面结冰形核过程的监测方法，具体涉及一种基于分子动力学的材料表面结冰形核过程的监测方法，属于形核动力学模拟技术领域。
背景技术
随着科学技术的发展，人类的活动范围已经不限于气候适宜的大陆，开始逐渐向高空、深海等气候恶劣的区域发展，其中，低温是恶劣环境中最具代表性的条件之一。材料表面结冰是低温环境中最容易发生的自然现象，结冰不仅会导致机器重量增加、外形改变、传动性能降低，还会干扰各种传感器的信号传输，所以低温环境造成的结冰现象往往是需要被避免的，例如在平流层中飞行的飞机往往需要克服-40℃～-50℃的低温，在这样的低温环境中，飞机的一些重要部件容易结冰，从而容易导致事故的发生。在低温条件下如何避免材料表面结冰现象的发生是近年来重要的研究课题。
若要想避免材料表面结冰现象的发生，就必须先了解结冰现象是如何发生的。研究者们发现，结冰过程分为形核过程和长大过程，其中，形核过程在结冰过程中起到了决定性作用，形核一旦完成，冰晶长大的速度是非常迅速的。所以，如何控制低温条件下水的形核过程，是解决结冰问题的关键。
水的结冰形核过程分为均质形核和异质形核，其中，均质形核发生的条件十分苛刻，需要很高的过冷度和纯净的水分子环境，现实条件下极少发生。所以，异质形核是现实生活中结冰形核的主要方式。异质形核易于发生主要是由于存在外界条件为形核过程提供能量，导致形核能垒降低，而各种材料表面是外界条件中最主要的部分。材料表面状态对异质形核过程的影响非常明显，例如，不同表面能和微结构分布会导致形核位置改变，同时会影响冰核的晶型分布。而控制结冰形核位置和冰的晶型对于防结冰材料的研究与开发是十分有利的。所以，研究材料表面状态对形核的影响机制对于新型防结冰材料的开发是十分必要的。但是，由于形核过程中初始晶核只有数百水分子，且在瞬间就有可能发生形核变化，这样的时间与空间尺度限制导致实际实验无法捕捉到水初始形核过程，因此，需要寻找新的理论研究方法。
分子动力学模拟技术的快速发展为微观形核机制的研究提供了重要的技术手段，从分子层面观察分析问题更加接近事物本质，不仅易于捕捉事物的微小变化，而且不需要耗费实际的资源，逐渐成为一种重要的研究方法。
发明内容
本发明的目的在于提供一种基于分子动力学的材料表面结冰形核过程的监测方法，该监测方法可以较好的捕捉到形核过程中的微小变化。
为了实现上述目标，本发明采用如下的技术方案：
一种基于分子动力学的材料表面结冰形核过程的监测方法，其特征在于，包括以下步骤：
Step1：用分子动力学仿真软件LAMMPS建立不同形状微结构与水的接触模型，或者，不设置微结构，直接建立平面基体与水的接触模型；
Step2：用分子动力学仿真软件LAMMPS对Step1中的接触模型进行初始化准备，得到平衡模型；
Step3：用分子动力学仿真软件LAMMPS对Step2得到的平衡模型进行可以使水分子发生结冰形核的动力学弛豫，在一个大气压下设定模拟温度为冰点以下的任意温度，同时设定合理的动力学弛豫时间，保证系统中的水分子完成形核结冰过程，设定合理的输出步长，利用dump命令输出各个状态点水分子与基体的坐标；
Step4：对于Step3得到的坐标信息，利用q6空间矢量判定程序判定每个分子的结冰状态，输出每一个状态点的结冰分子数量，得到数据后利用origin软件作图，得到结冰分子数量与时间的关系图，记录结冰分子数量结束震荡迅速上升的时间作为形核过程发生的时间，记录结冰分子数量停止上升趋于平缓的时间作为结冰过程结束的时间；
Step5：截取形核结冰过程发生时间段内的数据，输出每个状态点各个分子的结冰状态与坐标，得到分子状态与坐标信息Dat文件，将Dat文件输入TECPLOT软件中，绘制三维空间粒子分布云图，自定义不同的颜色分别表示三种不同的粒子，分析水分子形核的位置；
Step6：从Step5得到的Dat文件中提取已结冰的分子坐标，利用q3空间矢量判定程序判定每个结冰分子的晶型并分类输出晶型状态和坐标信息，统计立方晶型分子与六方晶型分子数量，得到分子状态与坐标信息Dat文件，将Dat文件输入TECPLOT软件中，绘制三维空间粒子分布云图，自定义不同的颜色分别表示立方晶型和六方晶型分子，统计两种晶型分子的数量并利用Origin软件绘图得到六方晶型分子与立方晶型分子数量随时间的变化趋势，分析结冰形核过程中晶型分布情况。
前述的基于分子动力学的材料表面结冰形核过程的监测方法，其特征在于，在Step1中，建立不同形状微结构与水的接触模型的方法具体如下：
采用三维周期盒子，首先在XY方向上建立无限长且厚度一定的平面基体，在平面基体上建立一定数量的微结构，固定微结构与平面基体的原子坐标形成材料基体，随后在微结构上建立一定数量的水分子，保证水分子与平面基体在截断半径内，最后进行能量最小化，使系统能量保持在合理范围。
前述的基于分子动力学的材料表面结冰形核过程的监测方法，其特征在于，微结构的形状和分布状态根据材料表面的微结构的形状和分布状态进行设定。
前述的基于分子动力学的材料表面结冰形核过程的监测方法，其特征在于，在Step1中，不设置微结构，直接建立平面基体与水的接触模型的方法具体如下：
采用三维周期盒子，首先在XY方向上建立无限长且厚度一定的平面基体，固定平面基体的原子坐标，随后在平面基体上建立一定数量的水分子，保证水分子与平面基体在截断半径内，最后进行能量最小化，使系统能量保持在合理范围。
前述的基于分子动力学的材料表面结冰形核过程的监测方法，其特征在于，在Step2中，对Step1中的接触模型进行初始化准备的方法具体如下：
赋予基体原子和水分子质量，粒子初始速度为常温下的粒子速度，采用SW粗粒化势描述基体与水分子之间的相互作用，根据模拟体系设定SW势函数参数，将系统处于300k、一个大气压下进行一定时间的可以使水分子发生结冰形核的动力学弛豫，系综选择为NPT，积分步长视情况而定，使系统达到平衡、密度合理。
前述的基于分子动力学的材料表面结冰形核过程的监测方法，其特征在于，SW势函数具体表现形式为：
式中，i、j、k代表三个相互作用的粒子；
ε为能量参数，表征粒子之间的相互作用强度；
σ为距离参数，表征截断半径；
θijk为i、j、k三个粒子所成的角度，θ0ijk为水分子的三个原子所成的角度；
r表示两个粒子之间的距离；
φ2描述二体势能量，φ3描述三体势能量，总能量等于二体势能量与三体势能量的和；
其余参数为无量纲的常数。
本发明的有益之处在于：
(1)克服了空间与时间尺度限制，可以直观敏锐捕捉到形核变化的发生，有利于从结冰源头理解结冰现象；
(2)通过监测形核位置和晶型分布，可以分析各种材料表面物理化学状态对形核结冰的影响，从而有助于实现从形核源头上控制结冰进程；
(3)给出了较为合理的分析方法，既可以保证获得正确的模拟结果，又可以减少计算量，提高了计算效率。
附图说明
图1是本发明提供的基于分子动力学的材料表面结冰形核过程的监测方法的流程图；
图2是基体与水的异质形核模型示意图；
图3是结冰原子数随时间变化示意图；
图4是Dat文件格式示意图；
图5是形核位置判断云图；
图6是初始晶核晶型判断云图；
图7是结冰过程中各晶型分子数量随时间变化示意图。
具体实施方式
为了对本发明的技术特征、目的和有益效果有更加清楚的理解和认识，现结合附图和具体实施例对本发明的技术方案进行以下详细说明，但不能理解为对本发明的保护范围的限定。
实施例1：监测粗糙碳表面水形核过程
参照图1，本发明提供的基于分子动力学的材料表面结冰形核过程的监测方法包括以下步骤：
Step1：建立平面基体与水的接触模型
用分子动力学仿真软件LAMMPS建立原子级粗糙碳基体与水的接触模型，具体的：
采用三维周期盒子，首先在XY方向上建立无限长且厚度为十层碳原子组成的平面基体，固定碳基体的原子坐标，随后在碳基体上建立3000个水分子，保证水分子与碳基体在截断半径内，最后进行能量最小化，使系统能量保持在合理范围。
Step2：对接触模型进行初始化准备
用分子动力学仿真软件LAMMPS对Step1中的接触模型进行初始化准备，得到平衡模型，具体的：
赋予碳原子(基体原子)和水分子质量，碳原子为12，水分子为18，粒子初始速度为常温(300k)下的粒子速度，采用SW粗粒化势描述碳基体与水分子之间的相互作用，SW势函数具体表现形式为：
式中，ε为能量参数，表征粒子之间的相互作用强度；σ为距离参数，表征截断半径；θijk为i、j、k三个粒子所成的角度，θ0ijk为水分子的三个原子所成的角度；r表示个粒子之间的距离；φ2描述二体势能量，φ3描述三体势能量，总能量等于二体势能量与三体势能量的和；其余参数为无量纲的常数。
根据模拟体系设定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。
将系统处于300k、一个大气压下进行30ns的平衡动力学弛豫，系综选择为NPT，积分步长为3fs，使系统达到平衡、密度合理。
碳基体与水的异质形核模型示意图见图2。
Step3：对平衡模型进行平衡动力学弛豫
用分子动力学仿真软件LAMMPS对Step2得到的平衡模型进行平衡动力学弛豫(结冰过程不限于使用平衡动力学方法，也可以使用降温动力学或其他可以使水分子发生结冰形核的动力学方法，具体方法可视使用者需求而定)，在一个大气压下设定模拟温度为冰点以下的任意温度(本实施例将模拟温度设定为230k)，同时设定动力学弛豫时间为150ns，保证系统中的水分子完成形核结冰过程，再设定每1000步输出一个状态点(步长越小，状态点越密集)，利用dump命令输出各个状态点水分子与碳基体的坐标。
Step4：判定结冰发生时间
对于Step3得到的坐标信息，利用q6空间矢量判定程序判定每个分子的结冰状态，输出每一个状态点的结冰分子数量，得到数据后利用origin软件作图，得到结冰分子数量与时间的关系图(见图3)，记录结冰分子数量结束震荡迅速上升的时间作为形核过程发生的时间，记录结冰分子数量停止上升趋于平缓的时间作为结冰过程结束的时间，结果为：系统在96.2ns时开始形核结冰，99.3ns时结冰完成。
Step5：分析形核位置
截取形核结冰过程发生时间段95ns-100ns内的数据，输出每个状态点各个分子的结冰状态与坐标，记录结冰分子为状态1、未结冰分子为状态2、基体原子为状态3，得到分子状态与坐标信息Dat文件(Dat文件中包含每个粒子的坐标信息和粒子的结冰状态)，如图4所示，将Dat文件输入TECPLOT软件中，绘制三维空间粒子分布云图，自定义结冰的水分子为黑色、未结冰的水分子为白色、基体为灰色，绘制得到的三维空间粒子分布云图见图5，通过分析发现：水分子初始形核的位置在碳基体表面。
截取形核结冰过程发生时间段内的数据可以减小计算量，截取数据时，在形核开始之前应多截取数纳秒数据，防止因读图不准错过关键状态点。
Step6：分析晶型分布和不同晶型分子数量
从Step5得到的Dat文件中提取已结冰的分子坐标，利用q3空间矢量判定程序判定每个结冰分子的晶型并分类输出晶型状态和坐标信息，统计立方晶型分子与六方晶型分子数量，水形核结冰时主要形成立方与六方两种晶型，其它晶型极少存在，所以只需要统计立方晶型分子与六方晶型分子数量，记立方晶型的分子为状态1、六方晶型的分子为状态2、基体原子为状态3，得到分子状态与坐标信息Dat文件，将Dat文件输入TECPLOT软件中，绘制三维空间粒子分布云图，自定义立方冰为黑色、六方冰为白色、基体为灰色，绘制得到的三维空间粒子分布云图见图6，统计两种晶型分子的数量并利用Origin软件绘图得到六方晶型分子与立方晶型分子数量随时间的变化趋势，绘制得到的变化趋势图见图7，通过分析发现：初始形核晶型为六方晶型且长大过程中晶型无变化。
在判定计算时需从Step5得到的Dat文件中截取已结冰粒子的坐标信息，再通过q3空间矢量判据进行晶型判断，否则将会产生大量的无效计算量。
实施例2：监测石墨烯表面水形核过程
参照图1，本发明提供的基于分子动力学的材料表面结冰形核过程的监测方法包括以下步骤：
Step1：建立平面基体与水的接触模型
用分子动力学仿真软件LAMMPS建立石墨烯基体与水的接触模型，具体的：
采用三维周期盒子，首先在XY方向上建立无限长且由四层石墨烯原子组成的石墨烯基体，固定石墨烯基体的原子坐标，随后在石墨烯基体上建立3000个水分子，保证水分子与石墨烯基体在截断半径内，最后进行能量最小化，使系统能量保持在合理范围。
Step2：对接触模型进行初始化准备
用分子动力学仿真软件LAMMPS对Step1中的接触模型进行初始化准备，得到平衡模型，具体的：
赋予碳原子(基体原子)和水分子质量，碳原子为12，水分子为18，粒子初始速度为常温(300k)下的粒子速度，采用SW粗粒化势描述石墨烯基体与水分子之间的相互作用，SW势函数具体表现形式为：
式中，ε为能量参数，表征粒子之间的相互作用强度；σ为距离参数，表征截断半径；θijk为i、j、k三个粒子所成的角度，θ0ijk为水分子的三个原子所成的角度；r表示两个粒子之间的距离；φ2描述二体势能量，φ3描述三体势能量，总能量等于二体势能量与三体势能量的和；其余参数为无量纲的常数。
根据模拟体系设定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。
将系统处于300k、一个大气压下进行30ns的平衡动力学弛豫，系综选择为NPT，积分步长为3fs，使系统达到平衡、密度合理。
Step3：对平衡模型进行平衡动力学弛豫
用分子动力学仿真软件LAMMPS对Step2得到的平衡模型进行平衡动力学弛豫，在一个大气压下设定模拟温度为冰点以下的任意温度(本实施例将模拟温度设定为230k)，同时设定动力学弛豫时间为150ns，保证系统中的水分子完成形核结冰过程，再设定每1000步输出一个状态点，利用dump命令输出各个状态点水分子与石墨烯基体的坐标。
Step4：判定结冰发生时间
对于Step3得到的坐标信息，利用q6空间矢量判定程序判定每个分子的结冰状态，输出每一个状态点的结冰分子数量，得到数据后利用origin软件作图，得到结冰分子数量与时间的关系图，记录结冰分子数量结束震荡迅速上升的时间作为形核过程发生的时间，记录结冰分子数量停止上升趋于平缓的时间作为结冰过程结束的时间，结果为：系统在97.3ns时开始形核结冰，100.1ns时结冰完成。
Step5：分析形核位置
截取形核结冰过程发生时间段96ns-101ns内的数据，输出每个状态点各个分子的结冰状态与坐标，记录结冰分子为状态1、未结冰分子为状态2、基体原子为状态3，得到分子状态与坐标信息Dat文件，将Dat文件输入TECPLOT软件中，绘制三维空间粒子分布云图，自定义结冰的水分子为黑色、未结冰的水分子为白色、基体为灰色，通过分析发现：水分子初始形核的位置在石墨烯基体表面。
Step6：分析晶型分布和不同晶型分子数量
从Step5得到的Dat文件中提取已结冰的分子坐标，利用q3空间矢量判定程序判定每个结冰分子的晶型并分类输出晶型状态和坐标信息，统计立方晶型分子与六方晶型分子数量，记立方晶型的分子为状态1、六方晶型的分子为状态2、基体原子为状态3，得到分子状态与坐标信息Dat文件，将Dat文件输入TECPLOT软件中，绘制三维空间粒子分布云图，自定义立方冰为黑色、六方冰为白色、基体为灰色，统计两种晶型分子的数量并利用Origin软件绘图得到六方晶型分子与立方晶型分子数量随时间的变化趋势，通过分析发现：初始形核晶型为六方晶型且长大过程中晶型无变化。
实施例3：监测石墨烯表面水形核过程
参照图1，本发明提供的基于分子动力学的材料表面结冰形核过程的监测方法包括以下步骤：
Step1：建立平面基体与水的接触模型
用分子动力学仿真软件LAMMPS建立石墨烯基体与水的接触模型，具体的：
采用三维周期盒子，首先在XY方向上建立无限长且由四层石墨烯原子组成的石墨烯基体，固定石墨烯基体的原子坐标，随后在石墨烯基体上建立3000个水分子，保证水分子与石墨烯基体在截断半径内，最后进行能量最小化，使系统能量保持在合理范围。
Step2：对接触模型进行初始化准备
用分子动力学仿真软件LAMMPS对Step1中的接触模型进行初始化准备，得到平衡模型，具体的：
赋予碳原子(基体原子)和水分子质量，碳原子为12，水分子为18，粒子初始速度为常温(300k)下的粒子速度，采用SW粗粒化势描述石墨烯基体与水分子之间的相互作用，SW势函数具体表现形式为：
式中，ε为能量参数，表征粒子之间的相互作用强度；σ为距离参数，表征截断半径；θijk为i、j、k三个粒子所成的角度，θ0ijk为水分子的三个原子所成的角度；r表示两个粒子之间的距离；φ2描述二体势能量，φ3描述三体势能量，总能量等于二体势能量与三体势能量的和；其余参数为无量纲的常数。
根据模拟体系设定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。
将系统处于300k、一个大气压下进行30ns的平衡动力学弛豫，系综选择为NPT，积分步长为3fs，使系统达到平衡、密度合理。
Step3：对平衡模型进行平衡动力学弛豫
用分子动力学仿真软件LAMMPS对Step2得到的平衡模型进行平衡动力学弛豫，在一个大气压下设定模拟温度为冰点以下的任意温度(本实施例将模拟温度设定为230k)，同时设定动力学弛豫时间为200ns，保证系统中的水分子完成形核结冰过程，再设定每1000步输出一个状态点，利用dump命令输出各个状态点水分子与石墨烯基体的坐标。
Step4：判定结冰发生时间
对于Step3得到的坐标信息，利用q6空间矢量判定程序判定每个分子的结冰状态，输出每一个状态点的结冰分子数量，得到数据后利用origin软件作图，得到结冰分子数量与时间的关系图，记录结冰分子数量结束震荡迅速上升的时间作为形核过程发生的时间，记录结冰分子数量停止上升趋于平缓的时间作为结冰过程结束的时间，结果为：系统在103.6ns时开始形核结冰，106.4ns时结冰完成。
Step5：分析形核位置
截取形核结冰过程发生时间段102ns-107ns内的数据，输出每个状态点各个分子的结冰状态与坐标，记录结冰分子为状态1、未结冰分子为状态2、基体原子为状态3，得到分子状态与坐标信息Dat文件，将Dat文件输入TECPLOT软件中，绘制三维空间粒子分布云图，自定义结冰的水分子为黑色、未结冰的水分子为白色、基体为灰色，通过分析发现：水分子初始形核的位置在石墨烯基体表面。
Step6：分析晶型分布和不同晶型分子数量
从Step5得到的Dat文件中提取已结冰的分子坐标，利用q3空间矢量判定程序判定每个结冰分子的晶型并分类输出晶型状态和坐标信息，统计立方晶型分子与六方晶型分子数量，记立方晶型的分子为状态1、六方晶型的分子为状态2、基体原子为状态3，得到分子状态与坐标信息Dat文件，将Dat文件输入TECPLOT软件中，绘制三维空间粒子分布云图，自定义立方冰为黑色、六方冰为白色、基体为灰色，统计两种晶型分子的数量并利用Origin软件绘图得到六方晶型分子与立方晶型分子数量随时间的变化趋势，通过分析发现：冰核晶型为六方晶型且长大过程中晶型无变化。
实施例4：监测表面具有规则阵列圆锥状微结构的石墨烯表面水形核过程
参照图1，本发明提供的基于分子动力学的材料表面结冰形核过程的监测方法包括以下步骤：
Step1：建立碳基体上圆锥状微结构与水的接触模型
用分子动力学仿真软件LAMMPS建立石墨烯基体上圆锥状微结构与水的接触模型，具体的：
采用三维周期盒子，首先在XY方向上建立无限长且由六层石墨烯原子组成的石墨烯基体，然后在石墨烯基体上构建规则阵列圆锥状微结构，该圆锥状微结构下表面半径为1.6nm、高度为1.2nm、间距为1.5nm，构建完圆锥状微结构后，固定石墨烯基体与圆锥状微结构的原子坐标形成材料基体，随后在圆锥状微结构上生成3000个水分子，保证水分子与石墨烯基体在截断半径内，最后进行能量最小化，使系统能量保持在合理范围。
Step2：对接触模型进行初始化准备
用分子动力学仿真软件LAMMPS对Step1中的接触模型进行初始化准备，得到平衡模型，具体的：
赋予碳原子(基体原子)和水分子质量，碳原子为12，水分子为18，粒子初始速度为常温(300k)下的粒子速度，采用SW粗粒化势描述石墨烯基体与水分子之间的相互作用，SW势函数具体表现形式为：
式中，ε为能量参数，表征粒子之间的相互作用强度；σ为距离参数，表征截断半径；θijk为i、j、k三个粒子所成的角度，θ0ijk为水分子的三个原子所成的角度；r表示两个粒子之间的距离；φ2描述二体势能量，φ3描述三体势能量，总能量等于二体势能量与三体势能量的和；其余参数为无量纲的常数。
根据模拟体系设定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。
将系统处于300k、一个大气压下进行30ns的平衡动力学弛豫，系综选择为NPT，积分步长为3fs，使系统达到平衡、密度合理。
对于材料基体与水接触的异质形核问题，水分子之间的相互作用由三体势描述，水分子与基体原子之间的相互作用由二体势描述，由于基体原子固定，不用设定基体原子间的作用势，所以可以减少计算量。
SW势函数中的参数可以视模拟体系的变化而变化，但总体而言，ε值越大，基体表面能越低，疏水效果越好，可通过改变ε值来改变基体的表面能状态。
Step3：对平衡模型进行平衡动力学弛豫
用分子动力学仿真软件LAMMPS对Step2得到的平衡模型进行平衡动力学弛豫，在一个大气压下设定模拟温度为冰点以下的任意温度(本实施例将模拟温度设定为230k)，同时设定动力学弛豫时间为200ns，保证系统中的水分子完成形核结冰过程，再设定每1000步输出一个状态点，利用dump命令输出各个状态点水分子与石墨烯基体的坐标。
Step4：判定结冰发生时间
对于Step3得到的坐标信息，利用q6空间矢量判定程序判定每个分子的结冰状态，输出每一个状态点的结冰分子数量，得到数据后利用origin软件作图，得到结冰分子数量与时间的关系图，记录结冰分子数量结束震荡迅速上升的时间作为形核过程发生的时间，记录结冰分子数量停止上升趋于平缓的时间作为结冰过程结束的时间，结果为：系统在137.5ns时开始形核结冰，141.4ns时结冰完成。
由于微结构有起伏，未结冰的水分子也可能被判定为结冰状态，但是总量很小，且不断震荡，所以结冰形核发生在结冰分子数量结束震荡开始迅速上升的时间段。
Step5：分析形核位置
截取形核结冰过程发生时间段136ns-142ns内的数据，输出每个状态点各个分子的结冰状态与坐标，记录结冰分子为状态1、未结冰分子为状态2、基体原子为状态3，得到分子状态与坐标信息Dat文件，将Dat文件输入TECPLOT软件中，绘制三维空间粒子分布云图，自定义结冰的水分子为黑色、未结冰的水分子为白色、基体为灰色，通过分析发现：水分子初始形核的位置远离石墨烯基体。
Step6：分析晶型分布和不同晶型分子数量
从Step5得到的Dat文件中提取已结冰的分子坐标，利用q3空间矢量判定程序判定每个结冰分子的晶型并分类输出晶型状态和坐标信息，统计立方晶型分子与六方晶型分子数量，记立方晶型的分子为状态1、六方晶型的分子为状态2、基体原子为状态3，得到分子状态与坐标信息Dat文件，将Dat文件输入TECPLOT软件中，绘制三维空间粒子分布云图，自定义立方冰为黑色、六方冰为白色、基体为灰色，统计两种晶型分子的数量并利用Origin软件绘图得到六方晶型分子与立方晶型分子数量随时间的变化趋势，通过分析发现：冰核晶型为立方晶型且长大过程中形成了部分六方晶型结构，立方晶型分子与六方晶型分子比例接近1：1。
实施例5：监测表面具有规则阵列圆台状微结构的石墨烯表面水形核过程
参照图1，本发明提供的基于分子动力学的材料表面结冰形核过程的监测方法包括以下步骤：
Step1：建立碳基体上圆锥状微结构与水的接触模型
用分子动力学仿真软件LAMMPS建立石墨烯基体上圆台状微结构与水的接触模型，具体的：
采用三维周期盒子，首先在XY方向上建立无限长且由六层石墨烯原子组成的石墨烯基体，然后在石墨烯基体上构建规则阵列圆台状微结构，该圆台状微结构上表面半径为0.4nm、下表面半径为1.6nm、高度为1.2nm、间距为1.5nm，构建完圆台状微结构后，固定石墨烯基体与圆台状微结构的原子坐标形成材料基体，随后在圆台状微结构上生成3000个水分子，保证水分子与石墨烯基体在截断半径内，最后进行能量最小化，使系统能量保持在合理范围。
Step2：对接触模型进行初始化准备
用分子动力学仿真软件LAMMPS对Step1中的接触模型进行初始化准备，得到平衡模型，具体的：
赋予碳原子(基体原子)和水分子质量，碳原子为12，水分子为18，粒子初始速度为常温(300k)下的粒子速度，采用SW粗粒化势描述石墨烯基体与水分子之间的相互作用，SW势函数具体表现形式为：
式中，ε为能量参数，表征粒子之间的相互作用强度；σ为距离参数，表征截断半径；θijk为i、j、k三个粒子所成角度，θ0ijk为水分子的三个原子所成的角度；r表示两个粒子之间的距离；φ2描述二体势能量，φ3描述三体势能量，总能量等于二体势能量与三体势能量的和；其余参数为无量纲的常数。
根据模拟体系设定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。
将系统处于300k、一个大气压下进行30ns的平衡动力学弛豫，系综选择为NPT，积分步长为3fs，使系统达到平衡、密度合理。
Step3：对平衡模型进行平衡动力学弛豫
用分子动力学仿真软件LAMMPS对Step2得到的平衡模型进行平衡动力学弛豫，在一个大气压下设定模拟温度为冰点以下的任意温度(本实施例将模拟温度设定为230k)，同时设定动力学弛豫时间为200ns，保证系统中的水分子完成形核结冰过程，再设定每1000步输出一个状态点，利用dump命令输出各个状态点水分子与石墨烯基体的坐标。
Step4：判定结冰发生时间
对于Step3得到的坐标信息，利用q6空间矢量判定程序判定每个分子的结冰状态，输出每一个状态点的结冰分子数量，得到数据后利用origin软件作图，得到结冰分子数量与时间的关系图，记录结冰分子数量结束震荡迅速上升的时间作为形核过程发生的时间，记录结冰分子数量停止上升趋于平缓的时间作为结冰过程结束的时间，结果为：系统在89.7ns时开始形核结冰，92.1ns时结冰完成。
Step5：分析形核位置
截取形核结冰过程发生时间段89ns-93ns内的数据，输出每个状态点各个分子的结冰状态与坐标，记录结冰分子为状态1、未结冰分子为状态2、基体原子为状态3，得到分子状态与坐标信息Dat文件，将Dat文件输入TECPLOT软件中，绘制三维空间粒子分布云图，自定义结冰的水分子为黑色、未结冰的水分子为白色、基体为灰色，通过分析发现：水分子初始形核位置在圆台微结构与石墨烯基体的交接处。
Step6：分析晶型分布和不同晶型分子数量
从Step5得到的Dat文件中提取已结冰的分子坐标，利用q3空间矢量判定程序判定每个结冰分子的晶型并分类输出晶型状态和坐标信息，统计立方晶型分子与六方晶型分子数量，记立方晶型的分子为状态1、六方晶型的分子为状态2、基体原子为状态3，得到分子状态与坐标信息Dat文件，将Dat文件输入TECPLOT软件中，绘制三维空间粒子分布云图，自定义立方冰为黑色、六方冰为白色、基体为灰色，统计两种晶型分子的数量并利用Origin软件绘图得到六方晶型分子与立方晶型分子数量随时间的变化趋势，通过分析发现：冰核晶型为立方晶型且长大过程中形成了少量六方晶型结构，六方晶型结构分子约占20％。
由此可见，本发明提供的基于分子动力学的材料表面结冰形核过程的监测方法，可以较好的捕捉到形核过程中的微小变化，通过晶型与形核位置的判定，可以分析材料表面状态对异质形核过程的影响，为材料表面形核结冰过程的研究提供了理论支撑，并为材料表面物理化学状态的设计提供了理论参考，具有较高的参考价值和理论指导意义。
综上所述，本发明提供的基于分子动力学的材料表面结冰形核过程的监测方法，借助分子动力学理论和仿真计算软件，对水在不同物理化学状态的材料表面上异质形核结冰过程进行模拟仿真，克服时间与空间尺度的约束，对水初始形核位置、形核结冰过程中的晶型变化进行监测，基于模拟结果分析不同材料表面结冰形核的过程，探索表面状态与异质形核结冰过程之间的关系，从而有利于从分子层面解释结冰机理，可为各种期望利用表面物理化学改性达到防结冰效果的实验探索、工艺优化与产品开发提供理论支持。
需要说明的是，上述实施例不以任何形式限制本发明，凡采用等同替换或等效变换的方式所获得的技术方案，均落在本发明的保护范围内。
Claims (6)
1.一种基于分子动力学的材料表面结冰形核过程的监测方法，其特征在于，包括以下步骤：
Step1：用分子动力学仿真软件LAMMPS建立不同形状微结构与水的接触模型，或者，不设置微结构，直接建立平面基体与水的接触模型；
Step2：用分子动力学仿真软件LAMMPS对Step1中的接触模型进行初始化准备，得到平衡模型；
Step3：用分子动力学仿真软件LAMMPS对Step2得到的平衡模型进行可以使水分子发生结冰形核的动力学弛豫，在一个大气压下设定模拟温度为冰点以下的任意温度，同时设定合理的动力学弛豫时间，保证系统中的水分子完成形核结冰过程，设定合理的输出步长，利用dump命令输出各个状态点水分子与基体的坐标；
Step4：对于Step3得到的坐标信息，利用q6空间矢量判定程序判定每个分子的结冰状态，输出每一个状态点的结冰分子数量，得到数据后利用origin软件作图，得到结冰分子数量与时间的关系图，记录结冰分子数量结束震荡迅速上升的时间作为形核过程发生的时间，记录结冰分子数量停止上升趋于平缓的时间作为结冰过程结束的时间；
Step5：截取形核结冰过程发生时间段内的数据，输出每个状态点各个分子的结冰状态与坐标，得到分子状态与坐标信息Dat文件，将Dat文件输入TECPLOT软件中，绘制三维空间粒子分布云图，自定义不同的颜色分别表示三种不同的粒子，分析水分子形核的位置；
Step6：从Step5得到的Dat文件中提取已结冰的分子坐标，利用q3空间矢量判定程序判定每个结冰分子的晶型并分类输出晶型状态和坐标信息，统计立方晶型分子与六方晶型分子数量，得到分子状态与坐标信息Dat文件，将Dat文件输入TECPLOT软件中，绘制三维空间粒子分布云图，自定义不同的颜色分别表示立方晶型和六方晶型分子，统计两种晶型分子的数量并利用Origin软件绘图得到六方晶型分子与立方晶型分子数量随时间的变化趋势，分析结冰形核过程中晶型分布情况。
2.根据权利要求1所述的基于分子动力学的材料表面结冰形核过程的监测方法，其特征在于，在Step1中，建立不同形状微结构与水的接触模型的方法具体如下：
采用三维周期盒子，首先在XY方向上建立无限长且厚度一定的平面基体，在平面基体上建立一定数量的微结构，固定微结构与平面基体的原子坐标形成材料基体，随后在微结构上建立一定数量的水分子，保证水分子与平面基体在截断半径内，最后进行能量最小化，使系统能量保持在合理范围。
3.根据权利要求2所述的基于分子动力学的材料表面结冰形核过程的监测方法，其特征在于，微结构的形状和分布状态根据材料表面的微结构的形状和分布状态进行设定。
4.根据权利要求1所述的基于分子动力学的材料表面结冰形核过程的监测方法，其特征在于，在Step1中，不设置微结构，直接建立平面基体与水的接触模型的方法具体如下：
采用三维周期盒子，首先在XY方向上建立无限长且厚度一定的平面基体，固定平面基体的原子坐标，随后在平面基体上建立一定数量的水分子，保证水分子与平面基体在截断半径内，最后进行能量最小化，使系统能量保持在合理范围。
5.根据权利要求1所述的基于分子动力学的材料表面结冰形核过程的监测方法，其特征在于，在Step2中，对Step1中的接触模型进行初始化准备的方法具体如下：
赋予基体原子和水分子质量，粒子初始速度为常温下的粒子速度，采用SW粗粒化势描述基体与水分子之间的相互作用，根据模拟体系设定SW势函数参数，将系统处于300k、一个大气压下进行一定时间的可以使水分子发生结冰形核的动力学弛豫，系综选择为NPT，积分步长视情况而定，使系统达到平衡、密度合理。
6.根据权利要求5所述的基于分子动力学的材料表面结冰形核过程的监测方法，其特征在于，SW势函数具体表现形式为：
式中，i、j、k代表三个相互作用的粒子；
ε为能量参数，表征粒子之间的相互作用强度；
σ为距离参数，表征截断半径；
θijk为i、j、k三个粒子所成的角度，θ0ijk为水分子的三个原子所成的角度；
r表示两个粒子之间的距离；
φ2描述二体势能量，φ3描述三体势能量，总能量等于二体势能量与三体势能量的和；
其余参数为无量纲的常数。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910074076.1A CN109885895B (zh) | 2019-01-25 | 2019-01-25 | 一种基于分子动力学的材料表面结冰形核过程的监测方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910074076.1A CN109885895B (zh) | 2019-01-25 | 2019-01-25 | 一种基于分子动力学的材料表面结冰形核过程的监测方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109885895A true CN109885895A (zh) | 2019-06-14 |
| CN109885895B CN109885895B (zh) | 2022-02-01 |
Family
ID=66926953
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910074076.1A Active CN109885895B (zh) | 2019-01-25 | 2019-01-25 | 一种基于分子动力学的材料表面结冰形核过程的监测方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109885895B (zh) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110782951A (zh) * | 2019-10-11 | 2020-02-11 | 安徽信息工程学院 | 基于分子动力学模拟的纳米液滴接触角获取方法 |
| CN111241680A (zh) * | 2020-01-10 | 2020-06-05 | 秒针信息技术有限公司 | 基于冷库温度变化的冷冻物品品质分析方法及装置 |
| CN113345530A (zh) * | 2021-07-02 | 2021-09-03 | 青岛科技大学 | 一种基于分子动力学的二元体系相互扩散系数模拟方法 |
| CN113643763A (zh) * | 2021-08-16 | 2021-11-12 | 南京航空航天大学 | 基于分子动力学的具有防结冰性能的图案化表面构建方法 |
| CN113836858A (zh) * | 2021-09-13 | 2021-12-24 | 深圳市紫光同创电子有限公司 | 芯片布局方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20120232685A1 (en) * | 2011-03-09 | 2012-09-13 | GM Global Technology Operations LLC | Systems and methods for computationally developing manufacturable and durable cast components |
| CN102750425A (zh) * | 2012-07-17 | 2012-10-24 | 哈尔滨工业大学 | 一种焊接过程热影响区组织演变的模拟方法 |
| US20170147723A1 (en) * | 2015-11-20 | 2017-05-25 | Metal Industries Research & Development Centre | Method of simulatively predicting a metal solidification microstructure for a continuous casting process |
| CN108664758A (zh) * | 2018-03-07 | 2018-10-16 | 苏州工业职业技术学院 | 一种低合金高强度钢碳化钛析出模拟仿真方法 |
-
2019
- 2019-01-25 CN CN201910074076.1A patent/CN109885895B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20120232685A1 (en) * | 2011-03-09 | 2012-09-13 | GM Global Technology Operations LLC | Systems and methods for computationally developing manufacturable and durable cast components |
| CN102750425A (zh) * | 2012-07-17 | 2012-10-24 | 哈尔滨工业大学 | 一种焊接过程热影响区组织演变的模拟方法 |
| US20170147723A1 (en) * | 2015-11-20 | 2017-05-25 | Metal Industries Research & Development Centre | Method of simulatively predicting a metal solidification microstructure for a continuous casting process |
| CN108664758A (zh) * | 2018-03-07 | 2018-10-16 | 苏州工业职业技术学院 | 一种低合金高强度钢碳化钛析出模拟仿真方法 |
Non-Patent Citations (3)
| Title |
|---|
| FITZNER M 等: "The Many Faces of Heterogeneous Ice Nucleation: Interplay Between Surface Morphology and Hydrophobicity", 《JOURNAL OF THE AMERICAN CHEMICAL SOCIETY》 * |
| SHEN Y 等: "Anti-icing potential of superhydrophobic Ti6Al4V surfaces: ice nucleation and growth", 《LANGMUIR》 * |
| 沈一洲: "Ti6Al4V超疏水表面的构建及其防/除冰机理研究", 《中国博士学位论文全文数据库工程科技Ⅱ辑》 * |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110782951A (zh) * | 2019-10-11 | 2020-02-11 | 安徽信息工程学院 | 基于分子动力学模拟的纳米液滴接触角获取方法 |
| CN110782951B (zh) * | 2019-10-11 | 2022-03-25 | 安徽信息工程学院 | 基于分子动力学模拟的纳米液滴接触角获取方法 |
| CN111241680A (zh) * | 2020-01-10 | 2020-06-05 | 秒针信息技术有限公司 | 基于冷库温度变化的冷冻物品品质分析方法及装置 |
| CN111241680B (zh) * | 2020-01-10 | 2023-12-08 | 秒针信息技术有限公司 | 基于冷库温度变化的冷冻物品品质分析方法及装置 |
| CN113345530A (zh) * | 2021-07-02 | 2021-09-03 | 青岛科技大学 | 一种基于分子动力学的二元体系相互扩散系数模拟方法 |
| CN113643763A (zh) * | 2021-08-16 | 2021-11-12 | 南京航空航天大学 | 基于分子动力学的具有防结冰性能的图案化表面构建方法 |
| CN113836858A (zh) * | 2021-09-13 | 2021-12-24 | 深圳市紫光同创电子有限公司 | 芯片布局方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109885895B (zh) | 2022-02-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109885895A (zh) | 一种基于分子动力学的材料表面结冰形核过程的监测方法 | |
| Ecke et al. | Turbulent rotating rayleigh–bénard convection | |
| Sulia et al. | Ice aspect ratio influences on mixed‐phase clouds: Impacts on phase partitioning in parcel models | |
| Fraedrich et al. | The Planet Simulator: Towards a user friendly model | |
| Petty et al. | Microwave backscatter and extinction by soft ice spheres and complex snow aggregates | |
| Harrington et al. | Cloud resolving simulations of Arctic stratus: Part II: Transition-season clouds | |
| Kritsuk et al. | On the density distribution in star-forming interstellar clouds | |
| Kuang | A moisture-stratiform instability for convectively coupled waves | |
| Andreas | Spray stress revisited | |
| Andrejczuk et al. | Numerical simulation of cloud–clear air interfacial mixing | |
| Kahnert et al. | Model particles in atmospheric optics | |
| Harrington et al. | A method for adaptive habit prediction in bulk microphysical models. Part II: Parcel model corroboration | |
| Gallagher et al. | On the influence of the Earth's rotation on geophysical flows | |
| Raymond | A model for predicting the movement of continuously propagating convective storms | |
| CN103294850A (zh) | 一种三维动态流体仿真算法智能匹配方法 | |
| Ovtchinnikov et al. | An investigation of ice production mechanisms in small cumuliform clouds using a 3D model with explicit microphysics. Part I: Model description | |
| CN109711094A (zh) | 一种基于分子动力学模拟微观结构表面防结冰性能的评价方法 | |
| Ničković et al. | Geostrophic adjustment on hexagonal grids | |
| Skamarock et al. | A fully compressible nonhydrostatic deep-atmosphere equations solver for MPAS | |
| Gonzalez et al. | Bubble nucleation in simple and molecular liquids via the largest spherical cavity method | |
| Rostami et al. | Evolution, propagation and interactions with topography of hurricane-like vortices in a moist-convective rotating shallow-water model | |
| Zhang et al. | The effects of surface kinetics on crystal growth and homogeneous freezing in parcel simulations of cirrus | |
| Martinez et al. | Thermodynamic consistency in the symmetric Poisson–Boltzmann equation for primitive model electrolytes | |
| Jones et al. | A very high resolution general circulation model simulation of the global circulation in austral winter | |
| Veron et al. | An Eulerian model for sea spray transport and evaporation |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |