基于Hyperworks和LS-DYNA的类蜂窝结构吸能效果分析和优化

Analysis and Optimization of Energy Absorption Effect of Honeycomb-Like Structures Based on Hyperworks and LS-DYNA

摘要: 蜂窝结构的材料具有相对密度低、韧性高、重量轻、隔热隔音效果好、吸能性能好等特点。因为蜂窝结构有这些优点，所以它被大量应用于航天、机械、包装、汽车等领域。以蜂窝结构为切入点，对其吸能性能进行研究时，从面内冲击、面外冲击两个方面进行了深入探讨。就前者而言，可将其优势充分反映出来，且在工程应用方面有着更大的价值。因此本文将研究蜂窝结构壁厚和折角对吸能性能的影响和优化分析。首先，基于有限元分析方法分析了蜂窝结构在面内的冲击载荷作用下，在固定区域内蜂窝结构胞体的壁厚、折角对蜂窝结构总吸能、总质量、比吸能的影响。最后，以结构壁厚和折角作为设计变量，以比吸能和质量作为优化目标，建立近似代理模型，运用第二代非支配排序遗传算法(NSGA-II)对多胞结构进行了多目标优化设计，得到了结构的pareto解集，经多目标优化后的结构达到了固定区域内最大总吸能、最小总质量。

Abstract: The materials of honeycomb structure have the characteristics of low relative density, high tough-ness, lightweight, good insulation and sound insulation effect, and good energy absorption perfor-mance. Due to these advantages, honeycomb structures are widely used in fields such as aerospace, machinery, packaging, and automobiles. When studying the energy absorption performance of honeycomb structures as a starting point, in-depth discussions were conducted from two aspects: In-plane impact and out-of-plane impact. As for the former, its advantages can be fully reflected and have greater value in engineering applications. Therefore, this article will study the impact and op-timization analysis of honeycomb structure wall thickness and bending angle on energy absorption performance. Firstly, based on the finite element analysis method, the influence of the wall thick-ness and bending angle of the honeycomb structure cell on the total energy absorption, total mass, and specific energy absorption of the honeycomb structure in a fixed area under the impact load in the plane of the honeycomb structure was analyzed. Finally, with the structural wall thickness and bending angle as design variables and the specific energy absorption and mass as optimization ob-jectives, an approximate surrogate model was established. The second generation non-dominated sorting genetic algorithm (NSGA-II) was used for the multi-objective optimization design of mul-ti-cell structures, and the pareto solution set of the structure was obtained. After multi-objective optimization, the structure achieved the maximum total energy absorption and minimum total mass in a fixed area.

1. 引言

为了减少交通事故造成的财产损失，现今人们将关注点放在汽车安全上。汽车在碰撞的过程中，大多通过汽车吸能部件塑变的方式，对汽车动能进行吸收，汽车的吸能部件，有保险杠、保险杠横梁、前纵梁、发动机罩、翼子板等，它们在汽车发生碰撞时会发生压溃变形。汽车出现碰撞时产生的能量，主要通过吸能部件压溃变形的方式进行吸收。就吸能部件而言，因受其变形空间大小的限制，所以通过压溃变形增大的方式进行碰撞能量吸收，这是不切合实际的。因此，为对汽车安全提供保障，可采取新型材料对既有传统材料进行替换 [1] [2] 。从汽车行业看，目前多孔材料逐步得到大范围推广使用，如蜂窝夹层材料、光洁材料，等等。这意味着此类材料有着相对突出的优势。

Gandhi与Atli以蜂窝结构为切入点，通过有限元仿真法的运用，对其单个胞体性能进行了探讨，以面内冲击为条件下，就其比吸能与夹角(直边和斜边之间的夹角)的关联性进行了分析 [3] 。Reid与Zou等人则以不一样的时段动态为切入点，选择了不同速度与变形模式两个维度，就其与蜂窝材料力学变化之间的关系进行了研究 [4] 。胡玲玲等运用显示动力学有限元分析方法研究了不同胞元结构的六边形蜂窝结构在面内冲击载荷下的力学性能 [5] 。在Griese与Schultz对蜂窝材料进行研究的过程中，在确保其相对密度不发生变化的基础上，通过响应面方法的运用，以高速面内冲击载荷影响为前提条件，就夹角与材料能量吸收二者的关系进行了优化分析 [6] 。Wu和Zheng等 [7] [8] 提出一种新型的功能梯度结构，研究其在动态载荷冲击工况下的吸能特性，同时运用基于径向基神经网络的近似模型，选用第二代非支配排序遗传算法对多胞管进行多目标优化设计。Wang和Peng [9] 等运用简化超折叠单元理论推导高速列车五胞吸能结构的平均压溃力理论模型，并基于优化拉丁超立方试验设计法，采用响应面模型寻找了结构的最优解，相比于原始结构，吸能量表现得到提升。

本文选择蜂窝材料进行面内冲击研究，选取总吸能、比吸能作为评价指标，对蜂窝的壁厚和折角的进行多目标优化。

2. 蜂窝结构面内吸能特性理论分析研究

2.1. 蜂窝材料变形机制理论分析

Figure 1. Hexagonal honeycomb material

图1. 六边形蜂窝材料

六边形蜂窝材料(如图1所示)的面内(X-Y平面内)的刚度和强度相比于面外的刚度和强度是低很多，由于面内冲击的原因而出现的力，造成孔壁出现弯曲弹性变形。本文以蜂窝结构为研究对象，以面内(X-Y内)冲击载荷影响为前提，基于此对其吸能性能进行研究，因此未就方向(面外)冲击进行分析。

从蜂窝材料看，若存在压缩或面内冲击的情况，先是会存在线弹性区域，且孔壁将表现为弯曲变形形态，基于此出现线弹性形变；若比其蜂窝临界应变值大，则孔穴将坍塌，且孔壁由于脆性断裂或者塑性屈服、弹性弯曲较大而出现坍塌。就这些变化而言，主要是由蜂窝材料的材质所决定的 [10] 。若当孔壁材料为线弹性材料时，孔壁会发生弹性弯曲，这是可以恢复；当孔壁材料为弹塑性材料时，孔壁会发生塑性屈服，不可恢复；当孔壁材料为脆性材料时，孔壁会发生脆性断裂，不可恢复。最终，当处在高应变时，孔穴坍塌至相邻的孔壁开始相互发生接触或是孔壁之间发生断裂堆积聚集在一起。孔穴因受力被挤压靠拢在一起，这导致了后段应力–应变的曲线急剧上升，如图2所示。

Figure 2. Compressive stress-strain curve of honeycomb materials

图2. 蜂窝材料的压缩应力–应变曲线

2.2. 蜂窝材料吸能评价参数

蜂窝材料在受到外部冲击的时候，冲击在蜂窝材料外的力就会做功，然后蜂窝材料自身产生变形吸收大量的冲击能量，来减缓碰撞带来的冲击。就此冲击能而言，最后会进行转化，得到动能与内能两种。此外，还有部分能量将由于撕裂、热效应与摩擦等原因而被消耗掉 [11] 。就蜂窝结构胞体而言，往往有着相对光滑的内部孔壁，因此如果冲击速度不大，可不考虑热效应及摩擦导致出现的消耗。从有限元仿真实验了解到，对于弹性或塑性材料而言，如果有着良好的韧性，出现断裂的几率并不大，因此也可不考虑由于断裂原因而消耗的能量。

这里针对蜂窝吸能结构的评价指标主要有吸能量(Energy Absorption, EA)，比吸能(Specific Energy Absorption, SEA)。

由力–位移曲线积分可得吸能量EA的表达式为：

$\text{EA}={\displaystyle {\int }_0^{d}F\left(x\right)\text{d}x}$ (1)

式中，d表示压溃位移；F(x)表示结构的压溃力随压溃位移变化的大小。吸能量越大，结构的吸能效果越好。

对于蜂窝材料，从有限元法计算方面看，其内能I、动能K求解方程为(2)，(3)。

$I=\frac{1}{2}{\alpha }^{\text{T}}{K}^{\ast }\alpha$ (2)

$K=\frac{1}{2}{\alpha }^{\text{T}}M\alpha$ (3)

上述公式中，K*、M分别表示单元刚度、质量矩阵， $\alpha$ 表示单元节点位移向量。

从蜂窝材料看，对吸能性能要求比较高，一方面要做到轻量化，另一方面要可以充分吸能。为对其吸能性能进行直观分析，对对蜂窝材料进行研究时，以单位质量为前提，基于此对其吸能情况进行分析，也就是比吸能(SEA)。其能够将在受到冲击的作用下蜂窝材料吸收能量过程所使用材料的比率反映出来。就比吸能而言，其计算公式如下：

$\text{SEA}=\frac{\text{EA}}{M}$ (4)

M为蜂窝材料的总质量，EA为在面内对蜂窝材料发生冲击时，发生位移S变形量吸收的总能量。其中，总质量M有公式为：

$M=4ABh{t}_h+4\left(A-1\right)\left(B-1\right)h{t}_h+16\left(A-1\right)Bl{t}_l$ (5)

$A=L_1/\left(2l\sin \left(\theta \pi /180\right)\right)+1$ (6)

$B=\left(L_2-h+2l\cos \left(\theta \pi /180\right)\right)/\left(2h-2l\cos \left(\theta \pi /180\right)\right)+1$ (7)

上式(5)，(6)，(7)中：t_h为竖直边壁厚，t₁为斜边壁厚，h为竖直边的长度，l为斜边长度，斜边与竖直边的夹角为θ，L₁为固定填充区域的长，L₂为固定填充区域的宽。如下图3，图4所示。

3. 蜂窝结构建模和有限元分析

3.1. 蜂窝模型的建立

通过CATIA进行建模然后将其导入Hypermesh中，如下图3所示。

Figure 3. Geometric model of honeycomb structure

图3. 蜂窝结构的几何模型

Figure 4. Geometric dimensions of honeycomb material cell body

图4. 蜂窝材料胞体几何尺寸

蜂窝材料单个胞体的尺寸都是相同的，单个胞体的尺寸边长h、倾斜方向的边长尺寸l，两边(尺寸边长h与倾斜方向的边长尺寸l)构成的折角为θ，胞体沿着纵向Z轴的长度为d，胞体的壁厚为t (斜边壁厚t₁，横直方向的壁厚t_h)，如图4所示。模型的长(X轴方向)为74个胞体，宽(Y轴上的方向)是21个胞体，在Y方向上L₂的尺寸为138 mm，沿着X方向L₁的尺寸为418 mm。对于单个蜂窝胞体，其相关参数具体如下：h = 4 mm，l = 4 mm，θ = 120˚，d = 4 mm，为确保其相对密度ρ* = 0.1，将模型壁厚设置成0.346 mm。

3.2. 蜂窝模型的前处理

在Hypermesh前处理软件对蜂窝结构进行前处理。对于该几何模型，选择了能够分析大变形的Belytschko-Tsay壳单元来划分网格。就该单元而言，是属于缺省的壳单元公式，通过面内单点进行积分，能够快速计算，耗费时间短，一般应用在大变形问题中，是最稳定有效的公式。对于六边形胞体而言，划分其每个边得到2单元的网格密度，基于此便能够确保模型能够更加精准地进行仿真处理。

Belytschko-Tsay壳单元相关参数设置如下：边长等于2 mm × 2 mm。模型由Belytschko-Tsay壳单元19,404个，节点24,447个组成，为确保收敛及精确地进行求解，单元厚度方向上采用5个积分节点，面内采用1个积分点，面内选择缩减积分；为避免模型之间发生穿透现象，刚性墙和蜂窝结构采用automatic_single_surface自动单面接触，刚性墙与蜂窝结构间采用automatic_surface_to_surface自动面面接触。静摩擦和动摩擦系数设置为0.1，以防止刚性块在压缩蜂窝结构时滑开。刚性墙均采用LS-DYNA中的MAT20 MAT_RIGID刚性材料，材料属性为默认值。

左端为带有速度V的运动刚性墙，右端为固定的刚性墙，如下图5所示，模型的上下两侧自由。蜂窝所用材料为铝，密度ρ为2. 7 × 10³ kg/m³，弹性模量69 GPa，泊松比 $\lambda =0.3$ ，屈服强度 ${\sigma }_y$ 为130 MPa，选用MAT24 MAT_PIECEWISE_LINEAR_PLASTICITY建立模型，因为应变率基本上不会影响铝，故模型构建过程中，无须就应变率与材料参数间的关系考虑在内。

Figure 5. Finite element model of honeycomb materials

图5. 蜂窝材料有限元模型

3.3. 有限元仿真结果分析

3.3.1. 有限元模型验证

通常选择高斯单元积分，就此积分而言，尽管能够让计算更加高效，省时，且能够更加精准、可靠地分析大变形，然而，选择该单点积分计算法，也存在不足之处，极有可能转化为零能模式，也就是沙漏模式 [12] 。主要原因是，选择单点积分计算方法时，单元中间部位出现积分点，就单元节点而言，位移为零的几率不大，然而在插值计算之后，便可获得等于零的应变值，导致内能出现为零的结果。这意味着，单元未出现变形。沙漏模式(零能模式)的出现导致了模型的能量不守恒，最终会影响有限元的计算精度甚至会出现计算错误，因此要将沙漏控制在一定的合理范围。判断标准是：因沙漏模式产生沙漏能不能超过系统总内能的10%。当超过这一判断标准时，我们就需要对沙漏能进行控制。所以下面将对建模完成的蜂窝结构进行沙漏能的验证。

提交有限元求解软件LS-DYNA进行计算，输出能量变化曲线如图6所示，仿真过程中，内能随仿真时间的增加而不断增加，沙漏能占总能量的比例为8.9%，能量变化合理，且曲线变化光滑，没有发生突变，因此该模型符合一定的精度要求。

Figure 6. Cellular energy and hourglass energy

图6. 蜂窝内能和沙漏能

3.3.2. 蜂窝结构的变形

下面将是对蜂窝结构模型施加一个速度为10 m/s面内冲击，如图7所示。

Figure 7. Deformation process of honeycomb structure under 10 m/s impact

图7. 蜂窝结构在10 m/s冲击下的变形过程

当速度在10 m/s时，在加载的初始0.015 s时刻，模型最左侧胞元的胞壁达到屈服应力开始失稳重叠，前面出现模糊“V”型，随着时间推移模糊的“V”型消失不见，取而代之的是“一”字型，胞元被压溃，堆积成类似于“一”字型，继续随着时间推移“一”字型压溃区域不断扩大直至将整个模型完全压溃。

Figure 8. Incomplete collapse

图8. 未完全压溃

Figure 9. Collapse

图9. 压溃

在蜂窝压缩的过程中，其倾斜的胞壁首先被前方横直的胞壁压至竖直，因为胞壁弯矩达到塑性极限弯矩时，就会发生塑性坍塌(即形成塑性铰)，如图8，随后更多的胞元进入塑性变形区，最后压溃整个胞体结构，直到压缩密实。如图9。

3.3.3. 壁厚对吸能效果的影响

首先按照表1中的数据完成蜂窝材料在刚性墙不同冲击速度下的仿真计算。

Table 1. Simulation parameter data table for different wall thicknesses

表1. 不同壁厚的仿真参数数据表

研究蜂窝胞体几何参数中壁厚对蜂窝材料总吸能EA和比吸能SEA的影响。将壁厚的大小作为变量，其余几何参数都不发生变化。蜂窝材料的总质量如表2所示。

Table 2. Total mass of honeycomb material M

表2. 蜂窝材料的总质量M

将总吸能和壁厚关系数据汇总如下图10。

Figure 10. The relationship between total energy EA and wall thickness t under different in-plane impact velocities

图10. 不同面内冲击速度下总能量EA与壁厚的t的关系

下面是关于在不同面内冲击速度下壁厚对蜂窝材料比吸能的影响，如图11所示。

Figure 11. The relationship between SEA and wall thickness t under different in-plane impact velocities

图11. 不同面内冲击速度下SEA与壁厚的t的关系

从图10可以看出，在10 m/s面内冲击速度下，随着壁厚t的增加，蜂窝材料结构的总吸能呈现着一种不断增大的趋势。导致这种不断增大的趋势原因是因为蜂窝材料的总质量M也再不断的增加。从图10可以看出，在50 m/s的面内冲击速度下，总吸能数值是大于在10 m/s面内冲击速度下的总吸能数值。随着面内冲击速度不断增大下，总吸能也在不断增多的。这主要的原因在于刚性墙的面内冲击速度不断增大，冲撞产生的能量转换成蜂窝材料的内能部分越来越多，内能的产生主要发生在胞壁弯曲的过程中，同时动能也因为速度的增大而增加，最后累加在总吸能EA上越明显，所以50 m/s面内冲击速度总吸能的“斜率”是比10 m/s面内冲击速度下的总吸能的“斜率”要大。

蜂窝结构的比吸能(SEA)如图11，随着壁厚t的增大不断的增大，这是因为总质量M在增加的时，总吸能EA也再增加。比吸能的大小等于总吸能比上总质量，所以比吸能SEA与壁厚的关系呈一种正比例增大关系。

3.3.4. 折角对吸能效果的影响

按照表3中的数据完成蜂窝材料在刚性墙不同冲击速度下的仿真计算。

Table 3. Simulation parameter data table for different bending angles

表3. 不同折角的仿真参数数据表

蜂窝材料的总质量M，如图12。

Figure 12. The relationship between the total mass M and the bending angle of honeycomb materials

图12. 蜂窝材料总质量M和折角的关系

从图12可以看出，蜂窝材料的总质量M随着折角增大呈现出先减小，到达120˚开始有缓慢增大的趋势，这是因为随着折角的增大，蜂窝胞体的整体宽度和高度都在不断的增大，导致在固定填充面积下蜂窝模型中胞体的数量减少，最后蜂窝模型的总质量M也随之减少。在折角达到120˚时，总质量变化出现了一个拐点，蜂窝材料的总质量M开始不断的增大，但由于蜂窝胞元高度不断减少，同时宽度也在不断的增大，所以蜂窝材料的总质量M增大的趋势缓慢。

将总吸能和折角关系数据汇总如图13所示。

Figure 13. The relationship between total energy absorption EA and bending angle under different in-plane impact velocities

图13. 不同面内冲击速度下总吸能EA与折角的关系

比吸能和折角关系如图14所示。

Figure 14. The relationship between specific energy absorption SEA and bending angle under different in-plane impact velocities

图14. 不同面内冲击速度下比吸能SEA与折角的关系

分析图13可以看出，在面内冲击速度为10 m/s和50 m/s下，随着折角的从80˚到90˚增大时，蜂窝结构的总吸能EA是呈现一种增大的趋势；在面内冲击速度为10 m/s折角在90˚到140˚之间整体呈现一种减小的趋势；在在面内冲击速度为50 m/s折角达到120˚时，相对于面内冲击速度10 m/s总吸能EA变化出现了一个拐点，蜂窝结构的总吸能EA开始有缓慢的增加的趋势。主要原因是：1) 受到蜂窝结构的总质量M的影响。2) 在刚性墙的不同冲击速度下，随着速度的增大，蜂窝结构的总能量也在不断增加，将刚性墙带来的大量冲击能量转化成蜂窝结构的总能量。

从图14可以看出，不管在冲击速度10 m/s，还是在50 m/s的冲击速度下，蜂窝结构的比吸能SEA都是随着斜边和竖直边间的折角增大，都显示出一种先增大后减小的趋势。当折角是90˚的时，蜂窝材料的比吸能SEA最大。

4. 蜂窝结构多目标优化

通过前面章节有限元仿真结果分析，可知蜂窝结构胞体的壁厚和折角对蜂窝材料的吸能性能有不同的影响。本章节对蜂窝结构进行优化分析，选取蜂窝材料中胞体的几何尺寸(壁厚和折角)为设计变量，设置变量的变化区间。其余胞体几何尺寸l、h、d为不变量。设定蜂窝结构在固定填充区域L₁ × L₂ = 418 mm × 138 mm内尽可能填充，蜂窝结构的总质量不超过0.5 kg，刚性墙的冲击速度为10 m/s，研究蜂窝结构压缩到80%时最大总吸能EA和在规定总质量范围内最大的比吸能。表达如式8所示：

$\{\begin{array}{l}\min \left(-\text{EA},M\right)\\ \text{s}\text{.t}\text{.}0.2\text{mm}\le t\le 0.6\text{mm}\\ {70}^{\circ }\le \theta \le {140}^{\circ }\\ 0\text{kg}\le M\le 0.5\text{kg}\end{array}$ (8)

根据上面章节的采集到的样本数据进行建立近似代理模型，选取常见的响应面模型(RSM)、克里格模型(KRG)、正交多项式模型(OPM)和径向基神经网络模型(RBF) [13] ，建立四种代理模型，对比不同代理模型的精度，如表4所示。

代理模型的精度通常由决定系数R2和均方根值RMSE评估，其中R2 < 0.2，RMSE > 0.9，则可认为拟合的模型具有一定的可靠性。通过表4四种代理模型对比可知，响应面模型对于模型相对具有较高的精度。

Table 4. Error accuracy evaluation of proxy model

表4. 代理模型的误差精度评估

代理模型建立之后，需选取优化算法对目标函数寻找最优解。第二代非支配排序遗传算法(NSGA-II)具有运算速度快、解集收敛性好的优点，目前的应用较为广泛。选择该算法对结构进行多目标优化。设置种群数为20，进化代数为50，最大迭代数为1000，交叉概率为0.9。求解得到多目标NSGA-II算法优化的帕累托前沿解，下图中红点时为最优解，如图15所示。

Figure 15. Pareto solutions for honeycomb structures

图15. 蜂窝结构的pareto解

计算模型理论值与仿真值的误差，依据图中最优点对应得到的壁厚t和折角θ修改模型。如表5所示，可以看出，EA与SEA的误差值均控制在± 5%以内，因此通过建立近似代理模型，利用优化算法寻找目标最优解的方法能够很好的适用于结构优化设计。

Table 5. Finite element simulation verification of pareto solution

表5. pareto解有限元仿真验证

5. 结论

本文对蜂窝材料在面内冲击作用下的吸能性能进行了研究，首先分析了六边形蜂窝结构的变形机制和蜂窝结构的吸能原理。通过有限元软件完成了有限元分析，分析了蜂窝结构中胞体的壁厚和折角对蜂窝结构的总质量、总吸能、和比吸能的影响，分析了不同冲击速度下对蜂窝结构的总吸能和比吸能的影响。最后对蜂窝结构依据前面得到的结论进行优化分析，以壁厚和折角为变量，以总质量最小和总吸能为最大为目标，在10 m/s面内冲击速度进行优化。本文研究得到的结论如下：

1) 从速度上来看，面内冲击速度的增加，蜂窝材料的总吸能和比吸能也随之增多，吸能效果越好。这主要的原因在于面内冲击速度增大，冲撞产生的能量转换成蜂窝材料的内能部分越来越多，同时动能也因为速度的增大而增加。

2) 从壁厚上来看，蜂窝壁厚对蜂窝材料的总吸能和比吸能的影响呈一种正相关性。随着壁厚的增大，总吸能和比吸能表现出一种不断增加的趋势。

3) 基于响应面模型(RSM)近似代理模型，采用第二代非支配排序遗传算法(NSGA-II)对结构进行了多目标优化设计，分别得出了帕累托前沿解，并且对结构在冲击速度为10 m/s时进行模型优化，优化完成同时得到有限元仿真结果与代理模型的误差在± 5%以内，说明基于近似代理模型的多目标优化设计能够有效应用于结构性能优化。

参考文献