基于表面效应三维纳米多孔金属力学性能的有限元分析1)

2023-12-16张永超糜长稳苟晓凡

力学学报 2023年11期

张永超糜长稳苟晓凡 ,3)

* (河海大学力学与材料学院,南京 211100)

† (东南大学土木工程学院,南京 211102)

引言

近10 年来,伴随着先进增材制造技术的快速发展,纳米多孔金属代表着多孔金属材料极端形态,兼具传统多孔结构和纳米金属材料的双重特性,受到广泛地关注[1-4].纳米多孔金属不仅具有高孔隙率和高比表面积的结构特点,在表面效应的作用下,还具备异于传统多孔金属的力学性能[5-7].因此,针对纳米多孔金属材料,结合微结构表面效应,探究模型参数与力学性能的影响规律,具有重要的工程实际意义.

纳米多孔金属材料力学性能的研究与纳米表面理论发展密切相关,最早可以追溯到Gibbs[8]的开创性工作.为在力学框架内衡量表面效应,Gurtin 和Murdoch 等[9-10]将纳米表面视作一层能够承受面内拉伸的零厚度薄膜,并建立了考虑表面拉伸刚度的表面弹性模型.此后,针对不同形式的纳米多孔金属材料,微结构表面效应的影响研究逐渐发展起来.按结构的维度次序,Feng 等[11]较早地将表面效应引入一维纳米梁,结合欧拉-伯努利梁理论,提出可预测开孔纳米多孔金属材料有效杨氏模量的微观力学模型.Dai 等[12]和He 等[13]分别讨论了远场载荷作用下二维弹性平面和三维弹性介质中纳米孔洞附近的应力集中效应.Wang 等[14]则研究了两个具有表面张力的纳米椭圆孔间的相互作用.

分子动力学方法可通过追踪原子相互作用来捕捉纳米材料微结构表面效应[15-20].对特征尺寸处于几十纳米以内的多孔金属材料,分子动力学方法被证明具有一定的指导作用,但当结构尺寸上升至几百纳米甚至更高时,分子动力学方法的计算成本急剧增加,计算效率也显著降低[21-22].值得庆幸的是,在合理网格密度前提下,基于连续介质力学基本原理的有限元法,其计算效率并不明显受到制造工艺、结构建模和模型尺寸等因素的约束[23-26].然而,传统有限元方法因尚不具备模拟纳米尺度下材料和结构力学行为的公式体系,因此难以精准捕捉纳米材料微结构表面效应[27-28].

尽管如此,依然有部分学者通过为孔洞表面定制有限单元的方式,尝试捕捉纳米孔洞表面力学响应.二维纳米孔洞表面可离散为一维单元,因此,Gao等[29]在传统有限元程序中引入一维表面单元,并以此研究了二维纳米多孔金属材料力学性能.Wang 等[30]采用与Gao 等[29]类似的方法,讨论了二维纳米圆孔附近应力集中现象,不同的是,他们通过用户单元二次开发,进一步将所开发的一维纳米单元植入商用有限元软件.通过以上分析可以发现,以往采用有限元方法对纳米多孔金属力学性能的研究,将计算模型简化为较为简单的二维结构,因此无法真实刻画纳米多孔金属材料的力学特性.截至目前,通过直接定制表面单元的方式来捕捉三维纳米多孔金属表面效应的研究还未见报道,依然缺乏针对一般三维纳米多孔金属结构的计算模型和分析方法.

据此,本文基于能量最小原理,将Gurtin-Murdoch表面模型引入单元刚度矩阵和单元平衡方程,通过为纳米表面构建有限元单元的方式,结合微观模型非均匀性,构建针对一般三维纳米多孔金属力学性能的有限元模型.采用所开发的计算模型,揭示孔隙率、孔洞数量和表面参数等因素对杨氏模量、压缩屈服强度和吸能性等力学性能的影响规律,为纳米多孔金属的材料制备和力学性能预测提供科学依据.

1 纳米多孔金属有限元模型

1.1 系统能量方程

考虑包含任意形状纳米孔洞的弹性体V,其中,Ω代表纳米孔洞区域,Γ=∂Ω 代表纳米孔洞边界.假设弹性体V受体力b的作用,在弹性体外边界St上受面力t的作用,如图1 所示,则基体材料应变能UB、表面弹性能US、外力做功W和弹性系统的总能量П满足

图1 包含任意形状纳米孔洞的弹性体力学模型Fig.1 Mechanical model of elastomer containing arbitrarily shaped nanoholes

其中,u表示外力作用引起的位移张量,γB和γS分别表示基体应变能密度和表面弹性能密度.弹性基体应变能密度和Gurtin-Murdoch 表面模型[9]定义的表面弹性能密度分别为

其中,γ0为表面自由能密度,εB,εS和τS分别为基体应变、表面应变和表面残余应力张量,DB和DS分别为基体和表面弹性刚度张量.

对式(1)进行变分后根据能量最小原理可得

将式(3)带入式(2)后进行变分,可得基体应变能、表面弹性能和外力做功的变分式

1.2 单元平衡方程

为刻画表面效应,对纳米多孔金属表面定制有限元单元.其中,基体单元的位移uB、应变 εB和表面单元的位移uS和应变 εS可分别表示为

其中,NB和NS分别为基体单元和表面单元的形函数,BB和BS分别为基体单元和表面单元的应变矩阵,分别为基体单元和表面单元节点位移.在孔洞表面,基体单元与表面单元的节点位移之间满足

其中,T为坐标转换矩阵,可将单元节点位移由全局坐标转换为局部坐标.将式(5)～式(8)代入式(4),可以得到系统能量方程的有限元格式

因此,离散系统的单元平衡方程可表示为

对于基体单元满足:c1=1,c2=0;对于表面单元满足:c1=0,c2=1.式(10)表示的单元平衡方程由基体单元刚度矩阵、表面单元刚度矩阵和单元节点力f3 个部分组成.其中,表面单元刚度矩阵部分包含表面弹性刚度DS的贡献,而表面残余应力(Gurtin-Murdoch 表面本构模型中应变无关项的表面残余应力)则与体力和面力一同被包含进单元节点力部分.因此,只要分别对基体和表面制定合适的单元,按式(11)计算出基体单元刚度矩阵、表面单元刚度矩阵和单元节点力f,再带入系统单元平衡方程进行迭代,即可计算出纳米表面力学响应.

1.3 表面单元定制

对于基体材料,ABAQUS 内置了丰富的体单元类型,因此不需要进行特别定制.三维纳米多孔金属的孔洞表面经离散后可视作二维单元,由式(11)可知,Gurtin-Murdoch 表面模型将纳米表面视为一层能够承受拉伸的零厚度薄膜,因此变形模式与二维面单元类似.经常使用的二维面单元有3 节点单元和四节点单元.其中,4 节点单元计算精度较高,3 节点单元则对复杂边界具有较强的适应能力,同时也更有利于有限元计算的快速收敛,可以达到节约计算成本的目的,因此,本文以3 节点三角形面单元为基础,来对Gurtin-Murdoch 型纳米表面进行离散,基体材料则采用四面体单元进行离散.如图2(a)所示,3 节点三角形面单元的i,j和k节点分别与4 节点四面体单元的I,J和K节点互为共有节点.

图2 表面单元和基体单元Fig.2 Surface element and bulk element

3 节点三角形面单元是二维单元,在利用FORTRAN语言对ABAQUS 进行二次开发时,可直接获取有限元系统全局坐标系下表面单元节点坐标,但无法直接获取局部坐标系下表面单元节点坐标,因此可分三步对表面单元进行定制:

(1) 将全局坐标系下表面单元节点坐标转换到局部坐标系中;

(2) 在局部坐标系下推导表面单元应变矩阵和单元刚度矩阵;

(3) 将局部坐标下应变矩阵和单元刚度矩阵转换到全局坐标系中,建立表面参数和单元参数关联机制,实现表面单元与三维基体单元的数据交换.

将全局坐标系下的单元节点坐标代入式(12),可以得到投影ars,随后根据式(13)(坐标转换方程)可得3 节点三角形面单元在局部坐标系下的节点坐标值.

其中,a1,a2和a3分别为局部坐标系原点O′在全局坐标系下的坐标.

在局部坐标系和全局坐标系下,3 节点三角形面单元的单元节点位移可分别表示为

二者之间满足坐标转换关系

其中,坐标转换矩阵满足

需要指出的是,在孔洞表面处,全局坐标系下表面单元节点位移和基体单元节点位移满足:对3 节点三角形面单元的应变矩阵在空间中进行扩充可得

其中,单元面积Atri为

常数bi和ci为

式中的 (i,j,m) 表示下标轮换操作,即i→j,j→m,m→i.因此,通过坐标转换可获得全局坐标系下应变矩阵

将式(16)～式(19)代入式(20)可得

其中,D和ttri分别代表3 节点三角形面单元的弹性刚度矩阵和单元厚度.

在使用FORTRAN 对ABAQUS 用户子程序进行二次开发时,单元应变矩阵和单元刚度矩阵是参与计算迭代和数据更新中最为重要的两个指标,可分别通过式(21)和式(22)计算得到.为准确构建纳米表面单元,还必须建立表面参数与3 节点三角形面单元材料参数之间的关联.

3 节点三角形面单元的物理方程可表示为

在平面应力问题中,弹性刚度矩阵可表示为

其中,Etri和vtri分别为单元的弹性模量和泊松比.

根据Gurtin-Murdoch 表面模型定义,纳米表面物理方程可表示为

其中,σS和 εS为表面应力张量和表面应变张量,IS为表面单位张量.将式(23)与式(25)进行对比,不难得出以下关系

利用式(21)、式(22)和式(26)构造表面单元刚度矩阵,随后连同表面残余应力,代入单元平衡方程即可完成对Gurtin-Murdoch 型纳米表面单元的构建.

1.4 有限元模型和材料参数

为纳米多孔金属选取合适的微观结构是多尺度力学框架的关键一环,分别选取含1/8 球孔、单球孔和随机多球孔的3 种模型进行研究,其中,1/8 球孔模型用来验证所开发有限元程序的有效性,单球孔模型用来分析表面效应对材料吸能性的影响,以及模型参数对材料杨氏模量和压缩屈服强度的影响,多球孔模型则包含一组具有相同孔径并随机分布的圆形孔洞,用来探究三维多孔金属材料的单轴压缩应力应变曲线、杨氏模量和压缩屈服强度.需要注意的是,本文所定义压缩屈服强度为压缩应力应变曲线屈服平台应力波动段中应力最低点对应的应力.图3 展示了上述3 种有限元模型的示意图.

图3 有限元模型Fig.3 Finite element model

为与参考文献[13]无限大空间中孔洞附近应力场作对比,将1/8 球孔模型的孔洞半径R0设为10 nm,边长设为20 倍孔径,以消除模型边界影响.如图3(a)所示,分别约束面N0N2N4N3沿x方向的位移、面N0N1N5N3沿y方向的位移和面N0N1N6N2沿z方向的位移.采用与参考文献[13]相同的材料参数,即基体材料为线弹性铁材,杨氏模量E和泊松比v分别取为177.33 GPa 和0.27,表面残余应力τ0为1.70 N/m,表面拉梅常数λ0和μ0分别取为-8.00 N/m 和2.50 N/m.

对单球孔和随机多球孔模型则采用周期性边界条件,按下式分别对立方体代表性体元边界面、棱和顶点处的节点位移进行约束[31]

2 结果与讨论

2.1 有限元程序校验

为验证所开发有限元程序的有效性,选取1/8球孔模型进行有限元仿真,沿z轴施加100 MPa 的单轴拉伸应力,如图4(a)所示.边界条件和材料参数按1.4 节设置,表面和基体单元分别按1.3 小节设置,单元总数约5 万个.

图4 1/8 球孔模型附近应力分布Fig.4 Stress distribution near 1/8 spherical hole model

图4(b)～图4(d)所示为分别考虑经典和Gurtin-Murdoch 表面模型的球孔附近(φ=0 和 φ=π) 径向、纬向和经向应力分布曲线,其中散点代表利用本文所开发计算模型得到的数值解,曲线代表He 等[13]所发展的解析解.通过比较可以发现,图中散点能够较为精准地落在对应的曲线上,这表明所开发的计算模型能够有效地捕捉球孔附近应力响应.此外,考虑表面效应的模型在孔洞表面出现了不同程度的应力集中,随着观测位置远离孔洞,考虑表面效应的应力分布曲线逐渐与经典工况重合,这说明表面效应对弹性体应力场的影响只在数倍于孔径的范围内较为明显.

图5 为球孔附近von Mises 应力云图.在经典工况下,应力集中出现在xy平面与孔洞表面的交汇处,而xz平面、yz平面和孔洞表面的交汇处出现了应力凹陷.相对于经典工况,考虑表面效应的孔洞表面整体应力水平有所提升,但应力集中发生在xz平面、yz平面和孔洞表面的交汇处(经典工况的应力凹陷处).由此可见,表面效应不仅加剧纳米多孔金属表面应力集中,还改变应力集中出现的位置.

2.2 单球孔模型分析

纳米多孔金属材料是一种优良的吸能材料,在压缩载荷作用下,经过较小的弹性变形阶段之后,会进入一个较长的应力平台阶段,在这个阶段,应力水平相差不大,但应变却快速增加,因此材料能在较为稳定的应力水平下吸收大量的能量,此后随应变进一步增加,应力迅速上升,材料进入密实化阶段.为衡量纳米多孔金属材料的吸能性,分别采用应变能密度vε和能量吸收率pe作为吸能指标

其中,εm代表任意应变,σm为对应的应力.由定义可知,应变能密度即单位体积材料所储存的应变能,反映纳米多孔金属的能量吸收能力,而能量吸收率为材料所吸收能量与对应应力之间的比值,当能量吸收率达到最大值时,表明在该应力水平下,材料的吸能效率最高.

首先对结构形式较为简单的单球孔模型力学性能进行研究,选取模型为10 nm×10 nm×10 nm 的立方体,孔隙率 ρ 分别设置为10%,20%,30% 和40%.沿z轴负方向施加最大7 nm 的位移载荷,按1.4 节设置周期性边界条件和对应的材料参数,表面和基体单元分别按1.3 小节设置,为保证计算精度,每个模型中单元数量不低于150 万个.对孔洞表面进行自接触设置,表面法向接触为硬接触,切向接触的摩擦系数为0.1.

图6(a)为单轴压缩下不同模型和孔隙率的纳米单球孔模型应变能密度随应变的改变.可以发现,曲线分为线性和非线性两区域,其中,当应变较低(小于0.4)时,由于对应的应力-应变曲线处于应力平台阶段,应力水平较为平稳,因此应变能密度随应变呈线性增长.当应变较高(大于0.4)时,由于模型进一步被压实,进入密实化阶段,压缩应力迅速增加,致使材料吸收更多的应变能,因此应变能密度随应变呈非线性增长.

图6 单轴压缩下不同表面模型和孔隙率的纳米单球孔模型应变能密度,能量吸收率随应变的改变以及孔隙率对表面效应的影响Fig.6 Variation of strain energy density energy absorption rate with seain,and effect of porosity on surface effects for nano single sphere hole model with different surface models and porosities under uniaxial compression

图6(b)为单轴压缩下不同模型和孔隙率的纳米单球孔模型能量吸收率随应变的改变.能量吸收率曲线可分为上升和下降两个区域,当应变量达到0.4 左右时,材料的能量吸收率达到最大,说明在应力平台即将结束时,材料的吸能效率最高.在密实化阶段,虽然材料吸能性得到进一步提升,但吸能效率却呈现出下降趋势,这同样也是由于压缩应力进一步增加导致的.

图6(c)展示了不同孔隙率下表面效应对模型应变能密度和能量吸收率的提升比例.相对于经典工况,表面效应降低了模型应变能密度,但提高了模型能量吸收率,产生这一现象的原因是考虑表面效应模型的应力-应变曲线中应力水平低于经典工况[33],因此吸收的总能量也低于经典工况,也正是因为较高的应力水平,导致模型吸能效率的降低(见式(28)).此外,随着孔隙率的增加,表面效应对应变能密度和能量吸收率的影响都逐渐增强,这是由于高孔隙率时模型比表面积较大,导致表面效应更加明显.

图7(a)为单轴压缩下纳米单球孔模型的杨氏模量和压缩屈服强度随孔隙率的变化.考虑表面效应的杨氏模量小于经典工况,而压缩屈服强度却与经典工况接近.在高孔隙率下表面效应对杨氏模量的降低效果更加明显,例如,当孔隙率为10%时,考虑表面效应的杨氏模量相对于经典工况降低了0.7%,而当孔隙率为40%时则降低了2.9%,这一规律与图6 中应变能密度随孔隙率变化类似,同样也是因为高孔隙率模型的比表面积更大导致的.

图7 纳米单孔模型杨氏模量或压缩屈服强度随孔隙率变化Fig.7 Variation of Young's modulus or compressive yield strength with porosity for nano single sphere hole model

如图7 中图例所示,为研究表面参数对纳米单孔模型力学性能的影响规律,以2.2 节中纳米铝表面参数为基准,选取不同量级的数值进行研究.采用控制变量法,在研究表面残余应力的影响时,令表面拉梅常数为零,在研究表面拉梅常数的影响时,令表面残余应力保持0.91 N/m 不变.

由图7(a)可知,表面效应对纳米单孔模型压缩屈服强度的影响较小,因此图7(b)～图7(d)仅展示采用不同表面残余应力和表面拉梅常数条件下纳米单孔模型杨氏模量随孔隙率的改变.可以发现,表面拉梅常数对杨氏模量的影响较小,相比之下孔洞表面残余应力对杨氏模量的影响较为显著.此外,从图7(b)可以发现,表面残余应力对杨氏模量的影响强烈依赖于外载荷方向,当加载方向为拉伸时,表面残余应力对材料杨氏模量有增强效果,而当加载方向为压缩时,表面残余应力则降低了材料杨氏模量.由于正向的表面残余应力有促使孔洞收缩的效果[13,34],因此在压缩载荷作用下,孔洞收缩促进了材料的压缩变形,这导致材料杨氏模量的减小.在拉伸载荷作用下,孔洞收缩抵抗了材料拉伸变形,这导致材料杨氏模量的增加.值得注意的是,宏观材料的杨氏模量属于材料固有属性,一旦确定结构形式和材料类型,其杨氏模量不受加载条件等外界因素影响,而通过以上分析可知,纳米结构的杨氏模量不仅与表面参数相关,还与加载方向密切相关,不同的加载方向甚至可以改变表面效应对杨氏模量影响趋势.

2.3 随机多球孔模型分析

为考虑模型非均匀性的影响,选取具有相同孔径的随机多球孔模型进行研究,模型仍取10 nm×10 nm×10 nm 的立方体,孔隙率 ρ 分别为10%,20%和30%,孔洞数量分别为5,10,15 和20 个,如图8(a)所示.沿z轴负方向施加最大值为6 nm 的单轴压缩位移载荷,边界条件、单元类型、单元数量和接触设置与2.2 节相同.

图8 单轴压缩下不同孔隙率和孔洞数的随机多球孔模型及应力-应变曲线Fig.8 Under uniaxial compression,random multispherical hole models with different porosity and number of holes,and stress-strain curves for models

图8(b)～图8(d)展示了单轴压缩下不同孔隙率的纳米随机多球孔模型单轴压缩应力-应变曲线.可以发现,单轴压缩应力-应变曲线出现了与传统多孔材料类似的弹性、应力平台和密实化3 个阶段.当孔隙率为10%时,考虑表面效应的应力-应变曲线和经典工况几乎重合.当孔隙率为20%和30%时,考虑表面效应的曲线应力水平明显低于经典工况,并且从放大图中可以观察到其与孔洞数量成负相关,而在经典工况下却无法观察到类似的规律.这说明相同孔隙率下,纳米多孔金属材料具有尺寸效应,减小孔洞尺寸会增强材料的表面效应,而宏观材料的应力-应变曲线是独立于孔洞尺寸的.

图9 为单轴压缩下随机多球孔模型杨氏模量和压缩屈服强度随孔隙率和孔洞数量的变化.可以观察到,表面效应不同程度上降低了材料的杨氏模量和压缩屈服强度,这与图8 中表面效应对曲线应力水平的影响规律类似.此外,在相同的孔隙率下,经典工况的模型杨氏模量和压缩屈服强度并不随孔洞数量改变,而考虑表面效应的模型杨氏模量和压缩屈服强度随孔洞数量的增加逐渐下降,这一现象与图7(a)保持一致,需要指出的是,图7(a)中表面效应对模型杨氏模量的影响并不像图9 那样明显,这是由于相同孔隙率下多孔模型具有更高的比表面积,因此表面效应更为显著.

图9 单轴压缩下随机多球孔模型的杨氏模量和压缩屈服强度随孔隙率和孔洞数量的变化Fig.9 Young's modulus and compressive yield strength of the random multi-sphere holes model in uniaxial compression vary with porosity and number of holes

由于缺乏可对比的实验数据,本文仅采用1/8球孔模型、单球孔模型和随机多球孔模型来对纳米多孔金属材料力学性能进行研究,所开发的有限元模型基于纳米表面理论和有限元原理,因此该方法并不受微观建模限制,仍可适用于一般纳米多孔金属结构.可以预见,对于具有复杂微观结构的纳米多孔金属材料,可借助专业建模程序(如SOLIDWORKS)和网格划分程序(如HYPERMESH)进行建模和离散,在此基础上结合本文所构建的有限元表面单元即可完成对复杂纳米结构力学性能的预测.