纳米多孔银力学性能表征分子动力学模拟∗

李杰杰 鲁斌斌 线跃辉 胡国明 夏热2)

1)(武汉大学动力与机械学院,水力机械过渡过程教育部重点实验室,武汉 430072)

2)(武汉大学动力与机械学院,水射流理论与新技术湖北省重点实验室,武汉 430072)

(2017年10月10日收到;2017年12月13日收到修改稿)

1 引 言

纳米多孔金属(nanoporous metals,NPMs)是一类具有高比表面积、纳米级孔棱直径和三维双连续微结构的新型功能性材料,兼具了传统多孔材料的结构特征以及纳米材料的优异性能.独特的微结构特征以及优良的物理、化学性能,赋予了纳米多孔金属在催化剂[1,2]、传感器[3]、执行器[4]、储能[5]、电化学驱动等[6,7]领域极大的应用潜力.而对其力学性能的研究和掌控,是推动纳米多孔金属功能化应用和发展的关键因素之一.

纳米多孔材料的相对密度、孔棱尺寸、拓扑几何结构等参数是决定其力学性能的重要因素,成为当前力学性能实验、理论及仿真研究的主要关注点之一.对于宏观孔隙结构的材料,Gibson和Ashby[8]给出了多孔固体力学模型和变形机理,涵括多孔材料的弹塑性形变、脆性坍塌、压缩变形等,定义了多孔材料力学性能随其几何参数变化的标度律.Liu等[9,10]探究了多孔材料在压缩载荷作用下的屈曲失效和剪切破坏模式.Diwu等[11]分析了多孔金属的应力应变关系及其位错发展规律.在纳米多孔材料研究中,Jin等[12]和Liu等[13]研究得出材料微结构是其力学性能的重要影响因素之一;发现实验制备的纳米多孔金样品的强度与硬度低于Gibson-Ashby标度律的预测值,并将这一现象归结于纳米多孔金样品中较低的孔棱连接性,同时引入有效相对密度概念,优化了Gibson-Ashby提出的有关强度与硬度的标度律,为合成高强度、高硬度的纳米多孔材料提供了思路.Zabihzadeh等[14]研究了纳米多孔多晶银的强度与孔隙率之间的相关性.Volkert等[15]通过纳米多孔金的单轴压缩实验分析纳米多孔金的力学性能,发现多孔材料的密度、结构以及孔棱和胞壁的绝对尺寸是影响材料力学性能的重要因素.Mangipudi等[16]将亚稳态分解过程形成的随机双连续结构、螺旋结构与真实纳米多孔金的整体力学性能进行对比,探究了拓扑结构和形态结构对标度律的影响.

然而,在现行纳米多孔金属力学性能分析中,理论及仿真所采用的模型各异,包括立方胞元模型[17]、球棍模型[18]、螺旋体模型[19]和随机双连续模型等[20].明晰不同微结构特征对力学响应的状态差异,以及模型选用对变形结果分析的影响趋势,将有助于建立更为统一、合理的力学分析体系.同时,通过调控几何结构使材料的力学性能达到最优是力学超材料的主要设计理念,而具有周期性结构的多孔银可视为力学超材料的一种,因此,基于Li和Gao[21]提出的“越小越强”的设计理念,探讨纳米尺度下微观拓扑结构对力学性能的影响,对于设计力学性能最优的纳米多孔金属具有重要指导意义,也将进一步推进其功能化应用.

基于此,本文采用分子动力学方法,分析了纳米多孔银(nanoporous sliver,NPS)在单轴拉伸载荷下的力学响应,着重考察了模型的拓扑结构差异引起的力学特性变化,并阐述了纳米多孔银力学性能对几何拓扑结构的依赖关系,在此基础上初步探讨了相对密度对其力学性能的影响.

2 模拟方案

2.1 模型建立

为分析相对密度对力学性能的影响,三种结构均取ρ=0.30—0.50之间五组不同相对密度值.相对密度ρ=ρ∗/ρb定义为孔棱处原子数与全体积原子总数的比值,其中ρ∗为纳米多孔材料的密度,ρb为块体材料的密度.三种结构的基体尺寸均取80a0×80a0×80a0,其中a0为银的晶格常数(a0=4.0898 Å,1 Å =0.1 nm).

图1 纳米多孔银3×3×3胞元结构(a)立方体结构;(b)金刚石结构;(c)螺旋体结构Fig.1.Assemblies of 3×3×3 unit cells of NPS:(a)Cube structure;(b)diamond structure;(c)gyroid structure.

在纳米多孔材料的研究进程中,基于立方胞元结构的研究为力学性能理论推导做出了巨大贡献[8,17].立方体多孔模型是在正六面体基础上,去除xyz三个方向中心对称的矩形柱得到最小重复单元,再利用最小重复单元3×3×3列阵得到.通过改变去除的矩形柱尺寸来控制立方体结构的相对密度,图1(a)为相对密度ρ=0.39的立方体结构纳米多孔银模型.

三重周期最小曲面(triply periodic minimal surface,TMPS)广泛应用于机械、化学、物理等造型研究,也为多孔结构模型的建立提供了基础[22].金刚石结构和螺旋体结构基于TMPS数学模型建立[23],由如下三角函数方程控制:

式中t为常数,用来控制相对密度大小;x,y,z为三个坐标方向;l为立方体晶胞边长.定义ϕ＞0部分为实体部分,ϕ＜0部分为孔隙部分,ϕ=0为实体部分与孔隙部分之间的分界面.图1(b)为相对密度ρ=0.40的金刚石结构纳米多孔银模型;图1(c)为相对密度ρ=0.40的螺旋体结构纳米多孔银模型.

2.2 计算方法

选用EAM(embedded atom method)势能函数[24],利用经典的分子动力学软件LAMMPS[25]开展纳米多孔银模型的单轴拉伸力学性能模拟.模拟过程中,运用Velocity-Verlet算法进行牛顿运动方程的数值积分,x,y,z三个方向均采用周期性边界条件,时间步长设为1 fs.加载前,所有模型采用共轭梯度法进行能量最小化,然后在等温等压(NPT)系综下使用Nosé-Hoover热浴使系统在300 K恒温条件下达到热平衡.平衡后,对纳米多孔银平衡体系进行单轴拉伸模拟.

单轴拉伸过程中,常用的系综为微正则(NVE)系综[26,27]和NPT系综[28,29],部分采用正则(NVT)系综[30].NVE系综下,体系能量保持不变,施加的单轴拉力对系统做功,相应原子的动能减少以抵消功的增加,从而保持系统能量恒定.拉伸过程中,系统纵向截面不发生变化,且温度有所下降,如果温度下降过快,材料将会发生脆断,应力应变过程较为接近实际情况.NPT系综下,在施加一个单轴拉力的基础上,再使用barostat来控制另两个非拉伸方向的压强为零,调整沿两个非拉伸方向的应力为零,且纵向截面随着拉伸方向均匀变化,不出现颈缩现象,整个模型在发生极大应变时仍未断裂,与实验中观察的现象存在一定的差异,不利于分析变形破坏机理.NVT系综下,在变形的同时利用thermostat保持温度恒定,在单轴拉伸模拟中使用较少.

基于此,在每个加载步中,x方向以109s−1的恒定应变率进行单轴拉伸,施加0.1%的工程应变增量.考虑到NVE系综下每个加载步后体系温度有所下降,每次拉伸完成后用Nosé-Hoover热浴在300 K的恒定温度下使系统松弛1 ps,松弛结束后进行下一步拉伸,循环反复,直至模型完全断裂.为减小热扰动波动等系统变量的影响,将每个加载步中最后100 fs的结果进行平均来分析处理数据.在模拟分析中,使用开放性可视工具OVITO(open visualization tool)[31],基于公共近邻分析法(common neighbor analysis,CNA)[32]分析纳米多孔银微结构的演化过程.

3 结果与分析

3.1 应力应变

图2 不同相对密度下的应力应变曲线 (a)立方体结构;(b)金刚石结构;(c)螺旋体结构Fig.2.Stress-strain curves of NPS with different relative densities:(a)Cube structure;(b)diamond structure;(c)gyroid structure.

图2(a)—(c)分别为不同相对密度立方体结构、金刚石结构和螺旋体结构纳米多孔银在单轴拉伸载荷下的应力应变曲线.曲线由三个主要部分组成,第一阶段为弹性形变区域,占整个应变过程较小一部分,应力随着应变的增加而急剧升高直至极限应力,应力应变呈现线性关系,斜率为模型的杨氏模量Es;第二阶段曲线进入塑性变形阶段,应力随着应变的增加而迅速降低,呈近似线性下降;最后阶段,应力随应变的增大而降低,下降速率减缓,此时,孔棱发生颈缩现象直至断裂.

对于同一结构,随着相对密度ρ的增加,极限应力σu不断增大.如图2(a)所示,立方体结构相对密度ρ从0.28增大到0.50时,极限应力σu从0.47 GPa增大到1.01 GPa,提升了114.9%.对于同一结构,在弹性应变阶段,随着相对密度ρ的增大,应力应变曲线不断上移,表明杨氏模量不断增加.如图2(b)所示,金刚石结构相对密度ρ从0.30增大到0.50时,杨氏模量Es从7.86 GPa增长到19.17 GPa,增长了143.9%.分析可知,纳米多孔银与宏观孔棱结构的多孔材料类似,基本力学性能在很大程度上依赖于相对密度.

3.2 变形行为

图3—图5分别为三种结构拉伸过程中拉伸方向的切片图,可以看出三种结构的拉伸变化过程大体相似.为探讨纳米多孔银的变形特征,以相对密度ρ=0.40时金刚石模型的切片为例,研究单个模型的变形行为,如图4所示,其中原子按照共同近邻原子法染色,灰色原子为表面原子,绿色原子为面心立方(face-centered cubic,FCC)原子,红色原子为密排六方(hexagonal closepacked,HCP)原子.图4(a)为模型初始状态的形貌,经过平衡处理,模型中存在极少位错.随着应变的增加,开始产生少量位错,HCP原子比例增加.在弹性阶段,孔棱的变形是完全可逆的.图4(b)中应变ε=0.062,此时存在大量HCP原子,纳米孔洞的半径开始增加,部分孔棱出现颈缩和层错.但大部分孔棱依旧能够承受载荷增长,此时材料表现为应力硬化行为,直至应力到达极限强度.应变ε=0.116时,HCP原子继续增加,此时应力已经越过峰值,出现大量层错及颈缩,模型即将从中间孔棱区域部分断裂,如图4(c)所示.随着应变的增加,最弱的竖直孔棱首先被拉断,导致模型净截面积减小,孔棱的断裂使应力重新分布,从而导致多孔材料连续变形和破坏.从图4(d)可以看到,当应变达到0.221时,所有竖直的孔棱均有破坏.拉伸变形过程中,塑性变形主要由孔棱的屈服主导,破坏发生在孔棱区域.

图3 立方体结构纳米多孔银单轴拉伸时原子的运动过程 (a)ε=0;(b)ε=0.062;(c)ε=0.209;(d)ε=0.391Fig.3.Atoms’movement in cube NPS under uniaxial tension:(a)ε=0;(b)ε=0.062;(c)ε=0.209;(d)ε=0.391.

图4 金刚石结构纳米多孔银拉伸时的原子运动过程 (a)ε=0;(b)ε=0.062;(c)ε=0.116;(d)ε=0.221Fig.4.Atoms’movement in diamond NPS under uniaxial tension:(a)ε=0;(b)ε=0.062;(c)ε=0.116;(d)ε=0.221.

图5 螺旋体结构纳米多孔银拉伸时原子的运动过程 (a)ε=0;(b)ε=0.062;(c)ε=0.209;(d)ε=0.310Fig.5.Atoms’movement in gyroid NPS under uniaxial tension:(a)ε=0;(b)ε=0.062;(c)ε=0.209;(d)ε=0.310.

图6为螺旋体结构纳米多孔银拉伸过程的切片示意图.图6(a)为初始状态,两螺旋状孔棱间夹角为111.5°;图6(b)为应变ε=0.138时的切片图,孔棱间夹角为120.1°,在拉伸过程中夹角变大,表明纳米多孔银在单轴拉伸变形过程中螺旋状孔棱逐渐被拉直.对比分析可以看出,螺旋体结构的高应力平台持续时间更久(图2(c)),塑性更好.此现象与螺旋体结构中螺旋状的孔棱有关,在单轴拉伸的过程中,螺旋形式的孔棱在受力拉直过程中抵抗变形,表现出较好的塑性.

图6 不同应变下螺旋体结构纳米多孔银螺旋状孔棱示意图 (a)ε=0;(b)ε=0.138Fig.6. Schematics of the typical spiral ligament in gyroid NPS under different strains:(a) ε =0;(b)ε=0.138.

3.3 力学性能的标度律

Gibson和Ashby[8]系统研究了多孔材料的力学性能,给出多孔材料的杨氏模量与相对密度之间的标度律:

式中Es为多孔结构的杨氏模量;Eb为块体模量,选用银的块体模量Eb=83 GPa[33];ρ为相对密度;CE和n为常数,大小取决于材料微结构,可以由实验得到.由(3)式可知杨氏模量Es与相对密度ρ的关系为Es=aρn.

图7 双对数坐标下不同结构的杨氏模量与相对密度之间的拟合关系Fig.7.Fitting curves of Young’s modulus as a function of the relative density for different morphological architectures in double logarithmic plots.

图7所示为三种结构在不同相对密度下杨氏模量Es的变化曲线图,在双对数坐标下,相对密度ρ与杨氏模量Es的拟合线为直线.根据拟合直线可得,立方体结构杨氏模量与相对密度之间满足Es1=51.5ρ1.71;金刚石结构满足Es2=50.6ρ1.58;螺旋体结构满足Es3=55.5ρ1.68.Gibson-Ashby给出关于杨氏模量的标度律,只考虑弯曲变形不考虑拉伸的贡献时,Es/Eb∝ρ2;考虑弯曲变形时,可以得到Es/Eb=(d/l)2,其中l为孔棱的长度,d为孔棱的直径.同时,相对密度与孔棱的尺寸有几何关系ρ∝(d/l)2,因此对于拉伸变形的贡献有Es/Eb∝ρ.对于开孔多孔材料而言,当相对密度较低时,以弯曲变形为主导;相对密度较高时,主要是拉伸变形在控制.根据幂指数关系拟合了多孔银模量与相对密度之间的关系,立方体结构的幂指数为1.71,金刚石结构的幂指数为1.58,螺旋体结构的幂指数为1.68,三者的值均介于1—2之间,表明在变形中既有弯曲变形也有拉伸变形.立方体结构和螺旋体结构的幂指数值较为接近,表明两种结构变形过程中弯曲变形和拉伸变形的影响程度相近;金刚石结构的幂指数值较接近1,变形过程中拉伸变形的影响程度相对其他两种结构稍大.由图7可知,金刚石结构和螺旋体结构的拟合直线几乎重合,两者的杨氏模量大小和变化趋势较为相似,原因在于两者均由三重周期最小曲面构成,两结构具有较高的相似性.同时,螺旋体结构和金刚石结构的杨氏模量大于立方体结构,主要是由于立方体结构简单,孔棱形式单一,抵抗拉伸变形的能力较弱.由此可得出,同一相对密度下,拓扑结构是影响纳米多孔银杨氏模量的重要因素之一,孔棱分布复杂的结构比形式单一的结构抵抗变形的能力更强,模量更高.

图8为不同模型的极限强度σu随相对密度的变化曲线图,强度与相对密度满足线性关系.图中直线为拟合直线,满足公式σu=Bρ+C,其中B和C是待定常数.根据图中拟合直线,立方体结构拟合参数B1=2.7,C1=−0.34;金刚石结构拟合参数B2=3.0,C2=−0.36;螺旋体结构拟合参数B3=2.4,C3=−0.31.对于单轴拉伸作用下的多孔材料,Gibson和Ashby指出其塑性变形主要来自两种典型的变形机理,一是连接点处的塑性坍塌,二是孔棱的屈服.当连接点处的塑性坍塌主导时,强度与相对密度呈二次幂指数关系;当孔棱的屈服为主导时,强度与相对密度呈线性关系.三种模型的强度与相对密度呈线性关系,表明纳米多孔银的屈服行为主要为孔棱的屈服.从3.2节可以看出三种结构的纳米多孔银的塑性变形主要依赖于孔棱,破坏发生在孔棱处,这与上述拟合得到的结论一致.由图8可以看出,在同一相对密度条件下,金刚石结构的极限强度最大,立方体结构次之,螺旋体结构最小.金刚石结构中存在大量相互交错的孔棱,孔棱与孔棱之间形成类似三角骨架的结构,如图9所示,该结构具有较好的稳定性,展现出相对高的极限强度.

图8 不同结构的极限强度与相对密度之间的拟合关系Fig.8.Fitting curves of ultimate strength as a function of the relative density for different morphological architectures.

图9 金刚石结构纳米多孔银三角骨架结构示意图Fig.9.Schematic of triangular skeleton structure of diamond NPS.

4 结 论

纳米多孔银的极限强度和杨氏模量均随相对密度的增大而增大.极限强度与相对密度近似为线性关系,塑性变形和破坏主要依赖于孔棱;杨氏模量与相对密度近似为幂指数关系,孔棱变形过程中不仅包含弯曲变形,还存在拉伸变形.

纳米多孔银的力学性能与拓扑结构具有紧密的相关性.螺旋体结构和金刚石结构均由三重周期最小曲面构成,结构存在一定的相似性,抵抗变形的能力相近,表现为模量值接近;立方体结构的孔棱形式单一,分布形式简单,模量值较小.金刚石结构的孔棱交错排列,形成类似三角骨架结构,塑性变形过程中具有一定的稳定支撑作用,表现出相对较高的极限强度.

螺旋体结构因存在螺旋形式的孔棱,在受力拉直的过程中抵消部分变形,表现出更好的塑性.

[1]Zhai X,Ding Y 2017Acta Phys.-Chim.Sin.33 1366(in Chinese)[翟萧，丁轶2017物理化学学报 33 1366]

[2]Wittstock A,Zielasek V,Biener J,Friend C M,Bäumer M 2010Science327 319

[3]Zhang L,Chang H,Hirata A,Wu H,Xue Q K,Chen M 2013ACS Nano7 4595

[4]Detsi E,Onck P R,de Hosson J T M 2013Appl.Phys.Lett.103 193101

[5]Ding Y,Zhang Z 2016Nanoporous Metals for Advanced Energy Technologies(Berlin:Springer Cham)pp83–131

[6]Ye X L,Liu F,Jin H J 2014Acta.Metall.Sin.50 252(in Chinese)[叶兴龙,刘枫,金海军2014金属学报50 252]

[7]Jin H J,Wang X L,Parida S,Wang K,Seo M,Weissmüller J 2010Nano Lett.10 187

[8]Gibson L J,Ashby M F 1997Cellular Solids:Structure and Properties(2nd Ed.)(Cambridge:Cambridge University Press)

[9]Liu P S 2010Acta Phys.Sin.59 8801(in Chinese)[刘培生2010物理学报59 8801]

[10]Liu P S 2010Acta Phys.Sin.59 4849(in Chinese)[刘培生2010物理学报59 4849]

[11]Diwu M J,Hu X M 2015Acta Phys.Sin.64 170201(in Chinese)[第伍旻杰，胡晓棉 2015物理学报 64 170201]

[12]Jin H J,Weissmüller J 2011Science332 1179

[13]Liu L Z,Ye X L,Jin H J 2016Acta Mater.118 77

[14]Zabihzadeh S,van Petegem S,Holler M,Diaz A,Duarte L I,van Swygenhoven H 2017Acta Mater.131 467

[15]Volkert C A,Lilleodden E T,Kramer D,Weissmüller J 2006Appl.Phys.Lett.89 061920

[16]Mangipudi K R,Epler E,Volkert C A 2016Acta Mater.119 115

[17]Feng X Q,Xia R,Li X,Li,B 2009Appl.Phys.Lett.94 011916

[18]Huber N,Viswanath R N,Mameka N,Markmann,J,Weißmüller J 2014Acta Mater.67 252

[19]Pia G,Brun M,Aymerich F,Delogu F 2017J.Mater.Sci.52 1106

[20]Sun X Y,Xu G K,Li X,Feng X Q,Gao H 2013J.Appl.Phys.113 023505

[21]Li X,Gao H 2016Nat.Mater.15 373

[22]Abueidda D W,Al-Rub R K A,Dalaq A S,Lee D W,Khan K A,Jasiuk I 2016Mech.Mater.95 102

[23]Yoo D J 2011Int.J.Precis.Eng.Man.12 61

[24]Daw M S,Baskes M I 1984Phys.Rev.B29 6443

[25]Pavia F,Curtin W A 2015Model.Simul.Mater.Sc.23 055002

[26]Yuan F,Huang L 2012J.Non-Cryst.Solids358 3481

[27]Vu-Bac N,Lahmer T,Keitel L,Zhao J,Zhuang X,Rabczuk T 2014Mech.Mater.68 70

[28]Luo J,Shi Y 2015Acta Mater.82 483

[29]Shen X,Lin X,Jia J,Wang Z,Li Z,Kim J K 2014Carbon80 235

[30]Pedone A,Malavasi G,Menziani M C,Segre U,Cormack A N 2008Chem.Mater.20 4356

[31]Stukowski A 2010Model.Simul.Mater.Sc.18 015012

[32]Faken D,Jónsson H 1994Comp.Mater.Sci.2 279

[33]Smith D R,Fickett F R 1995J.Res.Natl.Inst.Stand.Technol.100 119