闭孔泡沫金属几种不同建模方法的对比性研究

2020-08-24郭亚周刘小川白春玉郑志军王计真

航空材料学报 2020年4期

郭亚周, 刘小川, 白春玉, 郑志军, 王计真

（1.中国飞机强度研究所结构冲击动力学航空科技重点实验室，西安 710065；2.中国科学技术大学中国科学院材料力学行为与设计重点实验室，合肥 230026）

泡沫金属作为一种多孔、轻质的结构功能型材料，在受压缩过程中能够在较大范围内保持相对恒定的压缩应力，因此具有较好的能量吸收特性。目前泡沫金属作为一种理想的轻质吸能防护工程材料已经广泛应用于航空航天、汽车、建筑、舰船等领域中[1-3]。

有限元数值模拟作为一种解决工程实际问题的有利计算工具，能够在合适的建模方法基础上较好地还原物理实验现象，并且随着建模和分析方法的成熟能够在一定程度上逐步预测实验结果，节约研究成本，缩短研究周期[4]。为研究闭孔泡沫金属的压缩力学性能，摆脱实验工况的限制，探讨泡沫金属的细观胞元结构对整体力学性能的影响，不少学者逐步开展了闭孔泡沫金属的细观建模方法研究。张健等[5]基于微CT扫描影像信息建立泡沫金属二维细观有限元模型，研究了在高速加载下泡沫金属的压缩变形机理。袁本立等[6]从不同类型的泡沫金属单一胞元模型出发，研究了结构较为简单的单心立方模型、面心立方模型和体心立方模型表征下泡沫金属的力学性能，发现面心立方模型和体心立方模型能够相对较好地预测泡沫金属性质。Gibson等[7]提出构建了Gibson-Ashby泡沫金属模型，该模型的最大特点在于其为包含12条边的立体框架结构，通过模型分析后发现影响泡沫金属性质的因素主要是孔隙率、基体和泡孔结构。刘培生[8]提出了一种泡沫金属八面体模型，主要是为了克服Gibon-Ashby模型的不足，该模型结构的六个顶点分别是长方体两个平行平面和其他四条垂直棱边的中点，目前该模型在扭矩作用下泡沫金属的工程设计强度准则等方面具有较为广泛的应用[9]。随着泡沫金属细观建模方法的发展，不少学者发现Kelvin模型能够较为合理有效地模拟出泡沫金属的结构和性质，Kelvin模型是由十四条等长的棱长构成的正十四面体，该模型由6个四边形和8个六边形组成，能够在空间中实现无空隙堆叠[10]。宋延泽等[11]基于Kelvin模型研究了泡沫金属的动态压缩性能，研究结果表明，随着相对密度的变化，冲击速率和致密应变能够对变形模式产生较大影响。为了弥补Kelvin模型是由胞元规则堆叠的不足，学者逐渐提出了Voronoi模型，Voronoi模型包括二维和三维两种类型。Sun等[12]开发出了二维Voronoi模型用来描述泡沫金属的微观胞元结构。Zhang等[13]基于Abaqus分析软件用Voronoi模型研究了软件参数和质量缩放量、单元类型、单元大小对模型的模拟结果。

虽然目前关于泡沫金属的细观建模方法较多，但是每种建模方法对于不同场合的适用性却各不相同，且不同建模方法的合理性和表征准确性也各有差异。目前的研究中，更多的只是聚焦单一模型定性或定量的研究泡沫金属力学性能，而针对不同建模方法适用性和准确性的差异对比研究则相对较少。因此，本工作挑选出目前较为常用的三种模型来模拟泡沫金属的压缩过程，对比分析三种模型的变形模式、应力响应和压缩力学性能表征能力的差异性，并针对不同模型特点分别给出其研究适用场合建议。

1 模型构建

1.1 二维随机胞元模型

由于泡沫金属本身呈现近似各向同性的结构特征，且在通常的压缩实验中并不能够观察到试件内部的变形情况，因此为了较为直观地研究泡沫金属压缩内部结构响应过程，不少学者通过构建泡沫金属二维模型来研究其在不同工况下的变形模式和应力响应规律。

针对二维泡沫金属模型，一部分研究人员采用了二维Voronoi模型来表征泡沫金属的细观结构，虽然二维Voronoi模型能够实现较高的孔隙率，但是二维情况下的Voronoi模型与实际的孔型和相关分布具有一定差异性。基于此，本工作根据一个闭孔泡沫铝样件截面实际胞元的形状和分布情况（图1），构建了泡沫金属的二维随机胞元模型（图2）。模型的尺寸与泡沫铝试件的截面尺寸保持一致，平面尺寸为30 mm × 30 mm，厚度设为0.02 mm，模型的平均孔径为4 mm。其中若取泡沫金属基体材料密度为ρs，泡沫金属的密度为ρ，则泡沫金属的相对密度为：ρ* =ρ/ρs，此时二维随机胞元模型的相对密度为0.3。

图 1 泡沫金属模型（a）泡沫铝试样；（b）二维几何模型Fig. 1 2D metal foam model （a）aluminum foam sample；（b）2D geometric model

图 2 泡沫金属二维有限元网格模型Fig. 2 2D metal foam finite element mesh model

1.2 Voronoi模型

图3为三维Voronoi模型。三维Voronoi技术能够生成具有均匀胞壁厚度的闭孔泡沫金属模型，在给定的体积为V的区域内投放N个核点，其中任意两个核点之间的距离被约束为大于给定的最小距离tmin，最小距离定义为：

式中：tmin为任意两个核点之间的最小距离；k为胞元不规则度；t0为正十四面体任意两个核点之间的最小距离；N为给定的体积V的区域内投放的核点个数。通过改变k和t0的值即可生成不同体积和不同规则度的Voronoi模型。

图 3 Voronoi模型Fig. 3 Voronoi model

Voronoi模型的生成方式和泡沫金属发泡过程类似，因此三维的Voronoi模型相比于二维更能够还原真实泡沫金属胞元间孔棱和孔壁关系，已有的研究表明，当闭孔泡沫金属的相对密度低于0.3时，真实泡沫金属的拓扑结构与Voronoi模型基本相同[14]，因此Voronoi模型能够较好地表征泡沫金属。

本工作模型尺寸为30 mm × 30 mm × 30 mm，胞元直径为4 mm，胞元不规则度k为0.5。由于Voronoi模型全部由壳体构成，通过改变壳体厚度来改变模型的相对密度，如式（3）所示：

式中：ρs为泡沫金属基体材料密度；ρ为泡沫金属密度；ρ*为泡沫金属相对密度；V为模型体积；为模型所有胞壁面积总和；t为模型胞壁厚度。

1.3 Kelvin模型

如图4所示，Kelvin模型是由正十四面体堆叠阵列而得，该模型构建方法比Voronoi模型简单，目前仍然有不少学者采用这种模型来表征泡沫金属的结构。Kelvin模型本质上是不规则度k为0的Voronoi模型，因此Kelvin模型和Voronoi最大的区别在于Voronoi在规则胞元的基础上结合实际的泡沫金属发泡生成情况后对模型添加了扰动。虽然Kelvin模型并不具备胞元的随机性，但是就目前相关学者的研究结果上来看，Kelvin在某些工况下仍然和实验具备较好的一致性。

图 4 Kelvin模型Fig. 4 Kelvin Model

2 模型压缩分析方法

采用LS-DYNA软件对模型进行压缩模拟，其中二维随机胞元模型主要采用四边形网格，单元尺寸为0.5 mm，单元总数量是12353，模型在平面上施加对称约束。Voronoi模型全部由壳体构成，主要采用四边形网格，单元网格尺寸为0.5 mm，单元网格总数量是186743。Kelvin模型与Voronoi模型网格划分规则相同，单元网格总数量分别为229927和202528。在模型的加载方向两端放置两个刚性板，下刚性板固定，上刚性板加载恒定速率，设定压缩量为24 mm，当泡沫金属被压缩至应变为0.8时，计算终止。上下两个刚性板变形可以忽略，采用*MAT_RIGID模型，泡沫金属的基体材料选用*MAT_PLASTIC_ KINEMATIC来描述其力学行为，该模型为双线性应变硬化模型，另泡沫金属的基体材料为铝，铝基体材料参数如表1所示，泡沫金属模型采用*CONTACT_AUTOMATIC_SINGLE_SURFACRE接触，泡沫铝模型与刚性板之间采用*CONTACT_ AUTOMATIC_NODES_TO_SURFACE接触，摩擦系数设置为0.1，Song等[15]验证了该参数对模型的模拟结果几乎没有影响。

表 1 泡沫铝基体材料参数Table 1 Parameters of foam aluminum matrix material

3 模型压缩结果分析

3.1 二维随机胞元模型

由于二维细观模型孔壁孔棱与实际三维孔之间堆叠后的孔壁孔棱之间具有几何差异性，同时为了保证二维随机胞元模型的几何和网格特性，其孔棱不可能非常小，小的孔棱在几何上较难实现，同时也会导致网格尺寸和计算精度受到影响，因此二维随机胞元模型的相对密度 ρ*数值会偏大。

如图5所示，通过仿真计算后可得相对密度为0.3的泡沫金属在准静态下的压缩变形过程。由压缩变形图中可以看出，二维随机胞元模型能够便捷清晰地看出泡沫金属在压缩过程中内部胞孔和胞壁的变形情况。

在泡沫金属受到压缩的初级阶段，试件的整体变形较为均匀，试件内部胞壁的部分薄弱结构逐渐开始出现应力集中，进而在这些部位开始集中出现了相应的塑性屈曲，众多屈曲部位的塑性铰在宏观上表现为一条或者多条不规则的变形带。

随着压缩量的增大，变形带逐渐向周边延伸，变形带由条状逐渐向区域块状转变，此时孔壁之间会相互接触摩擦，下部孔壁为上部分坍塌孔壁提供支撑，使得整体试件的支撑能力能够在较大的应变内维持在一个较为稳定的水平，此压缩过程即为试样的平台压缩阶段。

当压缩应变增大到一定值之后，试样内部的胞元孔隙都基本上被压实，整个试件上胞孔的上下孔壁之间起到互相支撑的作用，这个阶段试件的压缩应力开始急剧上升并开始进入致密化阶段。

如图6所示，提取出二维随机胞元模型的应力应变曲线与文献[5]相同相对密度下的泡沫金属试样压缩实验结果对比，利用Miltz等[16]提出的吸能效率法来计算出响应的平台应力和致密应变，结果如表2所示。二维随机胞元模型与实验结果具有良好的一致性，但又由于内部胞元结构的差异性，导致其在压缩初期的弹性变形阶段具有一定的误差，从而使得平台段产生了差异性。

图 5 二维随机胞元模型准静态压缩响应过程Fig. 5 Quasi-static compression processes of 2D stochastic cell model （a）ε = 0；（b）ε = 0.15；（c）ε = 0.4；（d）ε = 0.6；（e）ε = 0.75

3.2 Kelvin和Voronoi模型

如图7所示，通过仿真计算后可得相对密度为0.1的Kelvin模型泡沫金属在准静态下的压缩变形过程，如图8所示为相对密度0.1的Voronoi模型泡沫金属在准静态下的压缩变形过程。由压缩变形图中可以看出，三维模型能够便捷清晰地看出泡沫金属在压缩过程中试样整体的变形和宏观局部失效情况。

在压缩初始阶段，由图7（b）所示，Kelvin模型的变形多集中在接触端，而中间位置变形较少。由图8（b）所示，Voronoi模型在相同应变下虽然在接触端仍然由一定的变形，但是在模型的中间部位则已经出现了较为明显的变形带，这与泡沫金属压缩实验现象更相符。

表 2 二维随机胞元模型与实验压缩平台应力和致密应变Table 2 2D random cell model and experimental compression platform stress and compact strain

图 6 二维随机胞元模型（ρ* = 0.3）和实验压缩结果对比Fig. 6 Comparison of 2D model（ρ* = 0.3）and experimental compression results

图 7 Kelvin模型准静态压缩响应过程Fig. 7 Quasi-static compression response processes of Kelvin model （a）ε = 0；（b）ε = 0.15；（c）ε = 0.4；（d）ε = 0.6；（e）ε = 0.75

随着压缩应变的增大，由图7（c）和图7（d）所示，Kelvin模型的变形带逐渐延伸至试样的中部，且由于Kelvin模型是由正十四面体堆叠而成，其试样内部薄弱结构较为集中且规律，因此压缩过程中呈现了较为规律的屈曲变形带，变形区域沿着变形带逐渐变大。由图8（c）和图8（d）所示，Voronoi模型与Kelvin模型略有不同，Voronoi模型由于内部胞元更具备随机性，其生成过程与泡沫金属的实际发泡过程类似，胞元形状和实际试样的孔型胞元结构特征相符，使得试样内部薄弱结构较为分散且呈现随机性，因此随着压缩应变的增大，Voronoi模型的变形带相比于Kelvin模型具有更明显的随机性。

在试样基本上被压实阶段，由图7（e）所示，Kelvin模型压缩过程每一层的胞元被较为规律地压实，最终呈现明显的分层现象。而如图8（e）所示，Voronoi模型则更加符合变形带的随机性原则，试样的各个部位基本上都发生了较大的变形。

如图9所示，分别提取出Kelvin模型和Voronoi模型的应力应变曲线与相同相对密度下的泡沫金属试样压缩实验结果对比。由图9就可以明显看出，Voronoi模型比Kelvin模型更接近实验结果，Kelvin模型屈服段与致密段与实验结果相比于Voronoi模型具有较大的误差。通过计算得出三维模型的平台应力和致密应变如表3所示，Voronoi模型与实验之间的误差几乎可以忽略不计，因此更能表征泡沫金属的压缩性能。

图 8 Voronoi模型准静态压缩响应过程Fig. 8 Quasi-static compression response processes of Voronoi model （a）ε = 0；（b）ε = 0.15；（c）ε = 0.4；（d）ε = 0.6；（e）ε = 0.75

图 9 三维模型（ρ* = 0.1）和实验压缩结果对比Fig. 9 Comparison of 3D model（ρ* = 0.1）and experimental compression results

表 3 三维模型与实验压缩平台应力和致密应变Table 3 3D model exprimental compression platform stress and compact strain

3.3 差异性分析

三维模型和二维模型具有一定的差异性，其中最主要的差异性主要体现在以下三个方面：（1）相对密度；（2）结构特征；（3）压缩变形响应。

从相对密度角度来讲，三维模型由于全部都是由壳体构成，因此通过改变壳体厚度即可在相对密度小于0.3的范围内任意改变泡沫金属模型的相对密度，而二维随机胞元模型由于受到孔壁孔棱的几何限制导致其相对密度往往不会小于0.25。在实际闭孔泡沫金属的使用场景中，为了达到质轻和高效吸能的目的，采用的闭孔泡沫金属相对密度大多较小。因此就相对密度方面来说，二维随机模型更偏向于泡沫金属的压缩规律性研究，而三维模型则更偏向于泡沫金属的应用优化研究。

从结构特征角度来讲，二维随机胞元模型以泡沫金属试件横截面为基准构建，由于泡沫金属具有一定的各向同性属性，因此在不考虑相对密度范围的情况下，二维随机胞元模型在可接受的误差范围内能够比三维模型更容易构建且能够更便捷地表征泡沫金属的压缩性能。三维模型更适合表征真实的试样状况，能够更好地表征泡沫金属的宏观压缩性能，但是其构建方法较为复杂，且Kelvin模型由于是由规则正十四面体堆叠而成，导致其压缩后呈现规则的分层坍塌，使得其虽然平台应力误差不大，但是压缩致密应变误差却较大。因此就结构特征方面来说，二维随机胞元模型和Voronoi模型都符合真实泡沫金属胞元的随机性，致密应变误差较小；而Voronoi模型和Kelvin模型整体空间结构与真实泡沫金属更相似，平台应力误差较小。

从压缩变形响应角度来讲，三维模型相比于二维随机胞元模型能够更全面地观察泡沫金属压缩的变形带产生与演变以及泡沫金属整体的变形模式变化，同时可以在多工况下考察模型不同应变下的应力场分布，其中Voronoi模型由于生成过程与闭孔泡沫金属发泡过程类似，内部结构更趋真实合理，其变形相比于Kelvin模型的分层压实来说更能体现变形带的随机性和真实性。二维随机胞元模型则较为直观地观察出泡沫金属内部胞元和胞壁的变形响应情况以及变形带的整体演变状况，但是就宏观变形来讲，较三维模型仍有较大不足。