相对密度和棱径尺寸可控的纳米多孔介质建模方法及系统与流程

本发明属于多孔材料结构设计技术领域，特别涉及相对密度和棱径尺寸可控的纳米多孔介质建模方法及系统。

背景技术：

随着当代材料科学的不断进步，纳米多孔金属由于其相对较低的相对密度和极高的比表面积，引发了科学家们极大的关注。为了追求不断提高的比表面积带来的优异催化性能和降低相对密度的同时保证较高的力学稳定性，人们尝试通过模板法和去合金法制备具备越来越小棱径尺寸和孔隙的纳米多孔金属。

然而，对于实验可以制备的几纳米棱径的纳米多孔金属，很难通过实验的方法获得其局部的微观力学性能及原位观察局部材料变形行为。然而，有着大量的试验和理论分析作为基础，以数值模拟为方法的研究也许可以以更生动的方式对多孔材料的性能进行表征。

采用分子动力学模拟的方法，无须试验样本的准备，就可以通过计算机采用数值的方法得到尺度在纳米量级上面的金属材料体系的大量模拟实验数据，分子动力学模拟方法是一门结合了数学、物理和计算机等多科学的复合型研究方法。分子动力学以由第一性原理计算出的势能关系和分子力学为基石，应用微分方程，通过对体系原子热力学状态的不断迭代，从而可以求得体系变化过程的数值解。故对于微观材料，尤其是纳米材料的力学等性能预测有着指导性的作用。

在当前的纳米多孔金属分子动力学模拟研究中，不同棱径尺寸(均小于5nm)和相对密度的单晶纳米多孔金模型被单向拉伸、压缩等以得到其力学性能，并且与Gibson等的经典公式做了拟合。可见，分子动力学能够提供原子尺度的相对大规模的计算效率。

运用分子动力学进行模拟的前提是要建立良好的模型，将数学建模得到的模型通过分子动力学仿真软件进行仿真模拟计算，得到模型的位置、能量、应力等热力学参量。

技术实现要素：

本发明的目的是提供相对密度和棱径尺寸可控的纳米多孔介质建模方法及系统。

本发明相对密度和棱径尺寸可控的纳米多孔介质建模方法，包括步骤：

S100对立方体模型在三维方向上划分网格，获得三维网格模型；

S200以孔洞为一相、介质为另一相，构建三维网格模型的Cahn-Hilliard方程，并采用有限差分法求解Cahn-Hilliard方程，获得具三维双连续多孔结构的三维网格模型；

S300设置截断阈值，利用阶跃函数计算三维网络模型当前时刻各格点的阶跃函数值，并判断格点为孔格点或棱格点，计算三维网络模型当前时刻的相对密度；

S400在三维网格模型各棱格点的内部生成随机割线，使用双正态分布函数拟合随机割线的长度概率密度分布，长度概率密度分布最大的极值点对应的长度即三维网络模型的平均棱径；

S500通过调整截断阈值，重复迭代执行步骤S100～步骤S400，直至三维网格模型的相对密度和平均棱径满足预设要求，即获得纳米多孔介质模型。

进一步的，步骤S300中，计算三维网络模型当前时刻的相对密度，进一步包括：

S310计算各格点当前时刻对应的阶跃函数值，具体为：逐一比较各格点的质量uijk,s和截断阈值uc,s的大小，uijk,s≤uc,s的格点的阶跃函数值为0，该格点记为孔格点；uijk,s＞uc,s的格点的阶跃函数值为1，该格点记为棱格点；其中，uijk,s表示三维网格模型中(i,j,k)位置格点当前时刻s的质量；uc,s表示当前时刻s的截断阈值；

S320将三维网络模型中所有格点对应的阶跃函数值相加并除以格点数，得到相对密度。

进一步的，步骤S400进一步包括：

S410在三维网格模型各棱格点附近生产随机位置的大量起点，分别以各起点为中心，生产大量的随机方向；

S420对各起点分别执行：以预设的长度步长，将起点向随机方向逐步增长，一旦到达距离任一孔格点的位置，即停止增长；之后，将起点向随机方向的反方向逐步增长，一旦到达距离任一孔格点的位置，即停止增长，从而获得随机割线；

S430存储所有随机割线的长度；

S440统计所有随机割线的长度的概率密度分布，并使用双正态分布函数拟合随机割线的长度概率密度分布，即双正态分布的长度概率密度分布，双正态分布的长度概率密度分布最大处的极值点所对应的随机割线长度，即平均棱径。

本发明相对密度和棱径尺寸可控的纳米多孔介质建模系统，包括：

第一模块，用来对立方体模型在三维方向上划分网格，获得三维网格模型；

第二模块，用来以孔洞为一相、介质为另一相，构建三维网格模型的Cahn-Hilliard方程，并采用有限差分法求解Cahn-Hilliard方程，获得具三维双连续多孔结构的三维网格模型；

第三模块，用来设置截断阈值，利用阶跃函数计算三维网络模型当前时刻各格点的阶跃函数值，并判断格点为孔格点或棱格点，计算三维网络模型当前时刻的相对密度；

第四模块，用来在三维网格模型各棱格点的内部生成随机割线，使用双正态分布函数拟合随机割线的长度概率密度分布，拟合出长度概率分布最大的极值点对应的长度，即三维网络模型的平均棱径；

第五模块，用来通过调整截断阈值，使得第一模块、第二模块、第三模块、第四模块重复迭代工作，直至三维网格模型的相对密度和平均棱径满足预设要求，即获得纳米多孔介质模型。

本发明具有如下特点和有益效果：

采用本发明方法可以获得指定相对密度和棱径尺寸的纳米多孔介质模型，通过对所获得的模型进行单向拉伸分子动力学模拟，可以预测纳米多孔材料的微观力学变形行为、弹塑性变形行为、缺陷演变行为等数据，从而为纳米多孔材料的优化设计提供指导数据，对于缩短纳米多孔材料的研发周期、减少研发费用等具有借鉴性意义。

附图说明

图1为本发明方法的具体流程图；

图2为Cahn-Hilliard方程在1000步和25000步时的格点状态，其中，图(a)和图(b)分别为在1000步和25000步时的格点状态；

图3为随机割线的长度概率密度分布图；

图4为实施例中获得的固定棱径、变化相对密度的一组模型；

图5为实施例中获得的固定相对密度、变化棱径的一组模型。

具体实施方式

为了更清楚地说明本发明和/或现有技术中的技术方案，下面将对照附图说明本发明的具体实施方式。显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图，并获得其他的实施方式。

应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

见图1所示，本具体实施方式以纳米多孔金属介质为例进行说明，具体步骤如下：

S100对立方体模型在三维方向上划分网格，获得三维网格模型。后文将“纳米多孔金属介质”简记为“介质”。

本具体实施方式中，所述纳米多孔金属介质模型为一立方体模型，将该模型分别在三维方向上分成N等分，即获得N³个小立方体，一小立方体即一格点。本具体实施方式中，N取100。

S200对三维网格模型采用有限差分法求解Cahn-Hilliard方程(卡恩-希利亚德方程)，获得有足够随机性的理想三维双连续多孔结构。

Cahn-Hilliard方程如下：

式(1)中：

u'代表某一相的浓度，则另一相的浓度为1-u'，-1≤u'≤1，u'的极值±1；本发明的三维网络模型中包括两相，一相为孔洞，另一相为介质；

f(u')代表自由能函数；

t代表过程进行的时间，即时间长度；

θ代表两相过渡区的宽度；

表示哈密顿算子，又称nabla算子。

应用有限差分法求解Cahn-Hilliard方程，有限差分法将微分式转化为差分式，同时利用双井势f(u')＝(u'²-1)²/4代替Cahn-Hilliard方程中的f(u')，使用迎风式差分对方程Cahn-Hilliard进行向前差分：

式(2)中：

和表示(i,j,k)位置格点分别在时间步数m和m+1时的浓度，其中，(i,j,k)为笛卡尔坐标；

τ代表迭代的时间步长。

时间步数和时间步长的乘积即相对零时的时间长度t。

在三维上采用正方形网格，拉普拉斯算子在笛卡尔坐标系上对空间采用二阶中心差分进行离散的格式为：

式(3)中：

b代表空间步长；

分别表示(i+1,j,k)、(i-1,j,k)、(i,j+1,k)、(i,j-1,k)、(i,j,k+1)、(i,j,k-1)位置格点在时间步数m的浓度。

(i+1,j,k)、(i-1,j,k)、(i,j+1,k)、(i,j-1,k)、(i,j,k+1)、(i,j,k-1)位置格点即(i,j,k)位置格点的邻域格点。

对某一特定的离散点(ib,jb,kb)，代表在t＝mτ时刻下该点的浓度状态。为了实现该离散方程的收敛，θ、b、τ分别设置为0.1、0.1、0.0001。采用周期性边界条件，数值求解三维体系中的格点，经过时间维度上的150000次迭代，每1000次输出一次其解的状态。图2是第1000次和第25000次迭代后的格点状态。

S300在三维网格模型某个时刻的状态中，设置截断阈值，计算当前三维网格模型的相对密度。

纳米多孔金属介质的相对密度定义为棱处原子的质量除以具有相同宏观体积的密实块体的质量。对三维网格模型s时刻的状态中，设置s时刻的截断阈值uc,s，设置三维网格模型中(i,j,k)位置格点在s时刻的质量为uijk,s，格点的质量和浓度呈线性关系，定义阶跃函数H(uijk,s-uc,s)，令：

Gijk,s表示阶跃函数H(uijk,s-uc,s)的值。当uijk,s≤uc,s时，H(uijk,s-uc,s)＝0，表示格点(i,j,k)为孔洞，该格点记为孔格点；当uijk,s＞uc,s时，H(uijk,s-uc,s)＝1，表示格点(i,j,k)为介质，该格点记为棱格点。

将网格中所有的N³个格点对应的阶跃函数值相加并除以N³，即得到三维网格模型的相对密度ρ，通过调整截断阈值uc,s，即可调整模型样品的相对密度：

S400在三维网格模型棱格点的内部生成随机割线，使用双正态分布函数拟合随机割线的长度概率密度分布，拟合出长度概率密度分布最大的极值点对应的长度，即具有随机双连续孔结构的三维网络模型的平均棱径。

如果将纳米多孔金属介质棱的形状看成一长度为l、直径为D的实心圆柱体，那么该圆柱体内的随机割线长度分布有两个尖锐的概率极值，分别对应l和D的长度。第一高极值点对应圆柱体的直径，第二高极值点对应圆柱体的长度。如果棱的长度与棱的直径相近的话，那么长度所对应的棱的第二高极值点会淹没在直径所对应的第一高极值点内。在这里引入使用两个正态分布的加和作为一种折中的方案来拟合任何统计方法得到的割线长度概率分布，拟合曲线的最大极值处所对应的割线长度对应棱的名义直径，见图3。

具体而言，本步骤进一步包括以下子步骤：

S410在三维网格模型的随机棱格点附近生成随机位置的大量起点，以起点为中心，生成足够多的随机方向；

S420以起点为中心，按随机方向以每步0.01倍格点间距d的增加量将随机割线长度增长，一旦到达距离任一孔格点的位置，即停止增长；以起点为中心，按随机方向的反方向，以每步0.01倍格点间距d的增加量将随机割线长度增长，一旦到达距离任一孔格点的位置，即停止增长；

S430将随机割线的长度添加至存储链表或数组；

S440统计存储链表或数组中所有长度的概率密度分布，并使用双正态分布拟合该概率密度分布，获得双正态分布的长度概率密度分布，双正态分布的长度概率密度分布最大处的极值点对应的随机割线长度，即平均棱径。

S500通过调整截断阈值uc,s，重复迭代执行步骤S100～步骤S400，直到生成的三维网格模型的相对密度和平均棱径满足设定要求。

S600输出相对密度和平均棱径满足要求的三维网格模型，即多孔金属介质模型。

本发明迭代结束后，可以得到满足要求的多孔金属介质的三维模型，图4中给出了算法输出的定棱径、变相对密度的纳米多孔金属介质模型，其中，图(a)～图(c)所示的模型分别表示棱径为4.3nm、相对密度分别为32％、38％、44％的多孔金属介质模型。图5给出了定相对密度、变化棱径的多孔金属介质模型，其中，图(a)～图(c)所示的模型分别表示相对密度为38％、棱径分别为3.6nm、4.3nm、5.0nm的多孔金属介质模型。

利用图4和图5的模型来生成单晶纳米多孔铜，运用美国桑迪亚国家实验室开发的LAMMPS软件包对生成的模型来进行仿真模拟计算。在加载之前，每一个样品通过共轭梯度法实现能量的最小化，然后通过NPT系综在300K下平衡50ps来实现内应力的最小化，用于共轭梯度法的相邻两个采样点的最大能量与应力变化设置值为10-15。单晶纳米多孔铜的拉伸采用循环的方式进行，在每一个循环内，首先对样品以10的工程应变速率进行总量为0.001的拉伸，然后在NPT系综下平衡1ps来稳定应力。应力值的提取采用每一个循环内对平衡时整个过程的应力值得平均值。整个拉伸过程循环往复，直至总工程应变达到0.30。

表1中所示的为模拟计算的单晶纳米多孔铜结构参数(相对密度、迟豫后的体系的尺寸、棱径)和力学性能参数(杨氏模量、抗拉强度、屈服强度)列表。其中样本I-V指的是固定棱径、变化相对密度的一组数值模拟实验，样本VI-IX指的是固定相对密度和变化棱径的一组实验。数值模拟实验的原子量为1,292,260至1,763,268个原子。杨氏模量的提取是将应力应变曲线应变为0-0.035阶段的部分做线性拟合得到的，屈服强度的提取是采用工程上常用的斜率为0.9倍杨氏模量直线和应力应变去线的交点所对应的应力值，抗拉强度的提取是直接统计应力应变曲线应力的最大值。可以注意到的是，整个体系均发生了较小的体积变化，均在0.6％以内，这种体系体积减少的现象来源于纳米多孔铜极大的比表面积带来的表面张力作用(毛细作用)，引起了整个体系的体积收缩。所模拟的体系的杨氏模量在几GPa量级上，屈服强度和抗拉强度在几百MPa量级上。

表1单晶纳米多孔铜结构及力学性能参数

可以发现的是，对于固定棱径，变化相对密度的一组数值模拟实验，样本的杨氏模量、屈服强度、抗拉强度均随着相对密度的变化发生了较大的变化，随着相对密度的增加，杨氏模量、屈服强度、抗拉强度均呈现非线性的增加，如相对密度从32.51％增加到44.18％，杨氏模量从3.261GPa增加到6.970GPa。对于固定相对密度，变化棱径的一组数值模拟实验，样本的杨氏模量随棱径的变化而变化幅度很小，屈服强度和抗拉强度随着棱径的增加而减小。

需要指出的是，分子动力学数值模拟的结果的得到的杨氏模量、屈服强度和抗拉强度等力学能参数可能比实验得到值要大，这是因为分子动力学中所采用的应力应变速率普遍大于试验中所采用的应力应变速率。然而，由分子动力学数值模拟得到的单晶纳米多孔铜的应力应变曲线的趋势和变形的机制是能够为我们理解纳米多孔金属的变形行为及力学响应提供洞悉和指导。