天然多孔材料如骨骼、贝壳、竹子、珊瑚、蜂窝等^[1]包含梯度、层级、螺旋、多孔以及自组装等结构，为新型结构设计及高性能材料研制提供了新的方法和途径。多孔蜂窝具有高强度密度比、独特的响应特性和较长的平台应力区，在受到外力作用时具有较好的吸能特性，因而成为理想的吸能材料，被广泛应用于建筑、汽车、机械、航空航天、化学工程、纳米制造以及生物医学等领域^[2]。传统蜂窝材料的静态与动态力学性能研究已经取得了丰硕的成果^[3-6]。

随着研究的不断深入，研究人员提出了更多的设计方式来进一步提升蜂窝结构的力学性能，如蜂窝填充泡沫^[7-8]、多层级蜂窝^[9]、梯度蜂窝^[10]、更换基体材料^[11-12]以及引入负泊松比^[13]等。其中蜂窝结构的梯度设计可以在不额外增加其他材料与蜂窝质量的前提下，仅仅依靠蜂窝质量分布的优化设计来提高其抗冲击性能，因此梯度蜂窝的抗冲击性能逐渐成为研究热点。研究蜂窝的能量吸收能力，需要首先掌握蜂窝的平台应力^[14]。研究表明，在动态载荷作用下，多孔结构的抗压强度随着加载速率的增大而显著提高^[15]。Reid等^[16]在研究木材类多孔材料时，基于冲击波模型建立了木材的动态压缩强度与准静态压缩强度之间的关系式。Ruan等^[17]通过模拟六边形蜂窝铝的面内冲击响应，研究了蜂窝铝动态压缩强度与冲击速度之间的关系。Qiao等^[18]对功能梯度蜂窝的动态冲击性能进行了研究，采用与Ruan等相似的方法，通过拟合参数得到了一种计算动态压缩强度的方法，该方法具有计算简单、便于推广等优点。

本研究受自然界王莲叶片脉络启发，提出一种兼顾强度与刚度的圆形辐射状仿王莲脉络梯度蜂窝芯层结构，利用ABAQUS有限元软件对其进行准静态与动态压缩数值模拟，分析其准静态压缩强度与相对密度的关系，并基于Ruan等^[17]和Qiao等^[18]提出的动态压缩理论，对面外梯度王莲仿生蜂窝结构在不同冲击速度下的动态压缩行为进行对比研究，探讨仿生蜂窝的动态压缩强度与相对密度、冲击速度之间的关系。

1.   王莲结构仿生蜂窝设计

王莲叶片极强的承载能力源于其独特的蜂窝状叶脉结构。叶脉相交形成蜂窝状，叶脉的粗壮程度和蜂窝状结构的尺寸大小从中心到边缘逐渐变化，构成既美观又能承重的整体，王莲叶脉分布如图1(a)所示。对王莲叶脉的测量和分析表明，王莲叶脉分布符合Rudwig植物形态学规律，各部分的长度比例满足黄金分割率^[19]。图1(a)与图1(b)显示了仿生王莲叶脉的拓扑演化过程，在图1(b)所示的简化模型中，王莲叶脉结构中的径向叶脉和环形叶脉分别用放射状的径向网和环形网表示。梯度蜂窝压缩模型的尺寸关系如图1(c)所示。

图  1  仿生蜂窝设计策略： (a)王莲叶脉， (b)王莲叶脉简化模型，(c)梯度蜂窝压缩模型

Figure  1.  Bionic honeycomb design: (a) natural Royal Water Lily vein structure, (b) a simplified model of natural Royal Water Lily vein, (c) compression model of gradient honeycomb

下载: 全尺寸图片 幻灯片

2.   梯度仿生蜂窝面外压缩模型的建立

2.1   有限元模型的建立

王莲仿生蜂窝的准静态压缩与动态压缩使用ABAQUS/Explicit软件模拟。根据图1(b)和图1(c)建立仿生芯层结构，冲击端面板与底部固定面板设置为刚体，壳单元的厚度为1 mm。在两个圆板的圆心处设置参考点，仿生蜂窝芯层放置在两板中间。仿生蜂窝半径R = 100 mm，各部分尺寸如图1(c)所示，总高度为21 mm，分为3层，即C1、C2、C3层，每层高7 mm。

仿生蜂窝和面板使用S4R壳单元进行网格划分，网格大小为1 mm，共219 240个单元，单元类型为有限薄膜应变，缩减积分，伪应变能控制设置为增强模式；接触设置为通用接触，用于避免前后面板与蜂窝发生穿透；摩擦系数设置为0.2；下面板参考点处设置为全约束，上面板参考点处设置一个固定压缩速度，如图2所示。

图  2  面外压缩模型

Figure  2.  Schematic diagram of out-of-plane compression model

下载: 全尺寸图片 幻灯片

准静态压缩的压缩速度为1 m/s，上面板以1 m/s的速度均匀压缩仿生蜂窝，下面板保持固定。动态压缩采用相似的设置，压缩速度v设置为10～250 m/s。

2.2   材料参数与相对密度

仿生蜂窝基体材料选择6060T4铝合金^[20]，具体的材料参数如表1所示，其中：$\,\rho$为蜂窝材料的密度，E为杨氏模量，$\;\mu$为泊松比，$\,\sigma$_cY为屈服强度，E_tan为切线模量。由于此类铝合金材料的应变率敏感性较低，因此有限元模拟过程中未考虑应变率效应。

表  1  6060T4铝合金材料参数

Table  1.  Material parameters of 6060T4 aluminum alloy

$\,\rho$/(kg·m−3)	E/GPa	$\,\mu$	$\,\sigma$cY/GPa	Etan/GPa
2 700	70	0.3	0.08	0.07

下载: 导出CSV 

| 显示表格

相对密度是蜂窝材料最主要的特征参数之一，用表观密度与基体密度的比值表示。仿生蜂窝按分叉的位置划分为5个部分，每个部分相对密度的计算公式为

式中：N_i为径向蜂窝壁的数量；R_{i, max}为第i部分的半径，当i = 0时，R_{0, max} = 0；t为壁厚，上标r表示径向，c表示环向。

表2为3种类型的王莲仿生蜂窝结构参数，分别为面外负梯度仿生蜂窝Case 1、面外正梯度仿生蜂窝Case 2以及均匀仿生蜂窝UG，3种仿生蜂窝的面内梯度均为均匀，面外平均相对密度均为5%。其中：Case 1仿生蜂窝的相对密度梯度从上到下分别为2%、5%和8%；Case 2仿生蜂窝的相对密度梯度从上到下分别为8%、5%和2%；而对于均匀仿生蜂窝UG，不同的相对密度可以通过改变UG-5%的壁厚获得。

表  2  模型类型与相关参数

Table  2.  Model and related parameters

Model	Wall thicknesses of C2 layer/mm						Relative density of out-of-plane/%
	Part 1	Part 2	Part 3	Part 4	Part 5		C1	C2	C3
Case 1	0.190	0.180	0.175	0.150	0.120		2	5	8
Case 2	0.190	0.180	0.175	0.150	0.120		8	5	2
UG-5%	0.190	0.180	0.175	0.150	0.120		5	5	5

下载: 导出CSV 

| 显示表格

2.3   有限元方法验证

评估模拟是否产生了正确的模态响应，最重要的验证方式是研究模型中的各种能量。采用显式动态算法来求解壳单元的准静态与动态压缩问题，需要进行速度敏感性验证以及能量平衡验证，同时需要检查以下能量：

（1）系统动能$E_{\rm{k}}$与内能$E_{\rm{i}}$的比值在5% ～10%之间；

（2）能量平衡，总能量为$E_{\rm{t}}$，即表示为$E_{\rm{t}}-(E_{\rm{i}} + E_{\rm{k}})\approx 0$；

（3）伪应变能$E_{\rm{a}}$与系统内能$E_{\rm{i}}$的比值对于一般的工程问题可以放宽到5%^[21]。

相对密度为2%的均匀仿生蜂窝的能量验证如图3所示。由图3(a)可以看出，准静态压缩时，系统总能量为系统内能与动能之和，动能约为内能的0.3%，表明整个压缩过程中模型均满足准静态模拟的要求。图3(b)给出了相对密度为2%的均匀蜂窝在冲击速度为250 m/s时的动态压缩过程中的能量平衡情况。动态压缩时伪应变能$E_{\rm{a}}$与内能$E_{\rm{i}}$的比值低于5%，可以认为模型是有效的。

图  3  相对密度为2%的均匀仿生蜂窝的能量验证： (a) 准静态压缩， (b) 250 m/s动态压缩

Figure  3.  Energy verification of uniform bionic honeycomb with relative density 2%: (a) quasi-static compression, (b) 250 m/s dynamic compression

下载: 全尺寸图片 幻灯片

3.   结果及分析

3.1   均匀仿生蜂窝准静态面外压缩

通过改变各部位的壁厚，得到不同相对密度的均匀仿生蜂窝。作为参考，首先对均匀仿生蜂窝的准静态面外压缩行为进行研究，主要研究准静态压缩下均匀仿生蜂窝的压缩强度。根据文献[22]，能量吸收效率$\eta \left( \varepsilon \right)$定义为

与面内压缩的仿生蜂窝和泡沫材料的密实化应变的定义不同，本研究中仿生蜂窝面外压缩密实化应变$\varepsilon$_D定义为将能量吸收效率-应变曲线的最大值所在的应变作为密实化应变的值

对于均匀仿生蜂窝而言，压缩强度$\,\sigma$_nY可以采用应力-应变积分与密实化部分的应变的比值求得，即

式中：$\varepsilon_0$为应力首次达到峰值时的应变。

图4中的黑色曲线显示了均匀仿生蜂窝在准静态加载下的典型压缩应力-应变曲线，应力采用名义应力计算：$\,\sigma$= F/S，其中F为冲击端的反力，S为仿生蜂窝的初始面积，S = $\text{π}$R²。应变$\varepsilon$ = $\delta$/h_c，其中$\delta$表示压缩距离，h_c表示仿生蜂窝的高度。由图4可知$\varepsilon_0$很小，由式(4)可以求出相对密度为2%的均匀仿生蜂窝的准静态压缩强度$\,\sigma$_nY为1.45 MPa。

图  4  均匀仿生蜂窝准静态压缩曲线与能量吸收效率

Figure  4.  Quasi-static compression curve and energy absorption efficiency of uniform bionic honeycomb

下载: 全尺寸图片 幻灯片

为了进一步研究均匀仿生蜂窝的准静态压缩强度，对不同相对密度仿生蜂窝的准静态压缩行为进行分析。王莲仿生蜂窝的胞元形状类似于四边形蜂窝，四边形蜂窝压缩强度与相对密度之间的关系可以表示为^[23]

式中：m = 1，$\overline \rho$为仿生蜂窝的相对密度，${\sigma _{\rm{Y}}}$为蜂窝基体材料的屈服强度。

面外压缩的六边形蜂窝在经历初始峰值之后具有软化的准静态响应，类似Ⅱ型结构的吸能特性^[24]，从图5可以看出，本研究建立的仿生蜂窝在经历了初始软化相之后，响应更类似Ⅰ型结构的吸能特征，即经过短暂的Ⅱ型结构响应之后，后续结构响应基本上与I型结构响应类似。同时，采用类似四边形蜂窝的压缩强度与相对密度的关系，根据图5可以求出不同相对密度对应的m值，可以发现m值随着相对密度的增大而增大，并且随着相对密度的增大，m值逐渐趋近于1.1，如图6所示。王莲仿生蜂窝的面外压缩强度优于四边形蜂窝、六边形蜂窝以及菱形蜂窝的面外压缩强度。考虑到计算误差，本研究计算采用的m值介于 0.9～1之间。

图  5  不同相对密度均匀仿生蜂窝准静态压缩曲线

Figure  5.  Quasi-static compression curves of bionic honeycomb with different relative densities

下载: 全尺寸图片 幻灯片

图  6  不同相对密度的m值与准静态压缩强度

Figure  6.  m and quasi-static compressive strengths of different relative densities

下载: 全尺寸图片 幻灯片

3.2   均匀仿生蜂窝面外动态压缩

动态压缩强度是体现多孔材料动态冲击性能的一个重要参数。Ruan等^[17]在研究六边形蜂窝的面内冲击时，基于Reid等^[16]建立的描述木材动态压缩强度的冲击波模型，建立了动态压缩强度与准静态压缩强度之间的关系式

式中：$\sigma _{\rm{nY}}^{\rm{D}}$为动态压缩强度，${\sigma _{\rm{nY}}}$为准静态压缩强度，$\,\rho$₀为仿生蜂窝的密度，${\varepsilon _{\rm{D}}}$为蜂窝的密实化应变，v为冲击速度。Qiao等^[18]将Ruan等建立的公式进一步改写为

式中：$\alpha = \overline \rho /{\varepsilon _{\rm{D}}}$，为蜂窝材料相对密度和几何参数的函数，$\;\rho$_s为基体材料的密度。研究结果表明：随着冲击速度增大，$\alpha$值的改变主要依赖于其相对密度，而对几何参数不敏感。

为了研究王莲仿生蜂窝的动态压缩强度，考虑到${\sigma _{\rm{nY}}}$与相对密度$\overline \rho$之间的关系，即$\displaystyle\frac{{{\sigma _{\rm{nY}}}}}{{{\sigma _{\rm{Y}}}}} = m\overline \rho$，将式(7)改写为

令$Y = \displaystyle\frac{{\sigma _{\rm{{nY}}}^{\rm{D}}}}{{\overline \rho }} - m{\sigma _{\rm{Y}}}$，$A = \displaystyle\frac{{{\rho _{\rm{s}}}}}{{{\varepsilon _{\rm{D}}}}}$。采用m-s-Pa单位制拟合不同冲击速度下仿生蜂窝的动态压缩强度与冲击速度的关系，本研究采用m = 1时A的值进行计算。图7显示当m = 1时，A的拟合推荐值为2 977，单位为kg/m³。根据A的定义可知：${\varepsilon _{\rm{D}}} = {\rho _{\rm{s}}}/A \approx 0.9$，此时的密实化应变与中高速冲击时的密实化应变接近，${\varepsilon _{\rm{D}}} \approx 0.85$，而低速冲击时的密实化应变要小一些。

图  7  动态参数A的拟合曲线

Figure  7.  Fitting curves of dynamic parameter A

下载: 全尺寸图片 幻灯片

综上所述，通过结合动态压缩强度的两种表示形式，可以得到相对密度、冲击速度共同影响下均匀仿生蜂窝的面外压缩强度公式，方便后续的理论计算。式(8)的关键参数为A，对A进行分析可知，图7的拟合参数A实际上隐含了Qiao等的研究结论，即随着冲击速度的增大，A对几何参数不敏感，表明式(8)在中高速冲击下是精确的，但是对于低速冲击，利用式(8)得到的压缩强度会产生一定的误差。对相对密度为4%与8%的均匀仿生蜂窝在不同冲击速度下的模拟结果与理论结果进行对比，结果如图8所示。由图8可以看出，均匀仿生蜂窝动态压缩理论值与模拟值对比具有较好的一致性，模拟解与理论解的相对偏差基本都在10%左右。理论解在较低冲击速度下与模拟解的差异较大，随着冲击速度的增大，两者之间的差异逐渐减小。另外，随着仿生蜂窝相对密度的增大，理论解与模拟解的相对偏差逐渐缩小。

图  8  均匀仿生蜂窝动态压缩强度理论值与模拟值对比：(a)相对密度4%，(b)相对密度8%

Figure  8.  Comparison between theoretical calculations and simulation values of dynamic compressive strengths ofuniform bionic honeycomb: (a) relative density 4%, (b) relative density 8%

下载: 全尺寸图片 幻灯片

3.3   面外梯度仿生蜂窝动态压缩

为了研究面外梯度仿生蜂窝的动态压缩性能，对由2%、5%、8%相对密度组成的两类面外梯度仿生蜂窝进行了动态压缩，仿生蜂窝分别表示为面外负梯度（C1、C2、C3层的密度梯度分别为2%、5%和8%）仿生蜂窝Case 1与面外正梯度（C1、C2、C3层的密度梯度分别为8%、5%和2%）仿生蜂窝Case 2。作为参考，对相对密度为5%的均匀仿生蜂窝UG-5%的动态压缩行为也进行了研究。同时，利用式(8)，对不同冲击速度下不同相对密度的仿生蜂窝的动态压缩强度进行了计算，并与模拟结果进行了对比。面外梯度仿生蜂窝在不同冲击速度下的压缩强度曲线与能量吸收效率曲线如图9所示。

图  9  面外梯度仿生蜂窝在不同冲击速度下的压缩强度曲线与能量吸收效率曲线

Figure  9.  Compressive strength curves and energy absorption efficiency curves of out-of-plane gradient bionic honeycombs under different impact velocities

下载: 全尺寸图片 幻灯片

由图9可以看出，不同面外梯度仿生蜂窝在低速（10 m/s）、中速（50 和100 m/s）和高速（200 m/s）冲击下的变形模式与压缩强度不同。在不同冲击速度下，Case 1与UG-5%蜂窝均为渐进压缩模式。在10 m/s冲击速度下，两种面外梯度仿生蜂窝的压缩强度曲线基本一致，表现为阶梯状模式，表明两者的变形模式类似，即准静态抗压强度较低的芯层依次被压缩。Case 1蜂窝表现为C1-C2-C3层渐进压缩，而Case 2蜂窝则是靠近固定端的C3层先被压缩，表现为C3-C2-C1层渐进压缩。当冲击速度增大到50 m/s时，两种面外梯度仿生蜂窝的压缩强度曲线开始不同，在初始冲击阶段，Case 2蜂窝会产生较高的初始应力峰值，而Case 1蜂窝的初始应力峰值较低。这是由于在惯性作用下，Case 2蜂窝中准静态压缩强度最高的C1层会优先被轻微压缩，产生大约6%的应变，而Case 1蜂窝的C1层准静态压缩强度最低，因此产生的初始峰值也低。随后Case 2与Case 1蜂窝的应力-应变曲线再次重合，表明压缩过程重新变为准静态抗压强度较低的芯层依次被压缩。在100 m/s冲击速度下，仿生蜂窝的变形模式与冲击速度为50 m/s时基本一致，但Case 2蜂窝的初始应力峰值更高，C1层持续压缩时间变长且压缩量增至18%。在200 m/s冲击速度下，两种面外梯度仿生蜂窝压缩强度曲线的差异明显变大，表现为Case 1蜂窝压缩强度曲线依旧保持为阶梯状，而Case 2蜂窝的压缩强度曲线呈现“U”形变化。这是由于在高速冲击下，惯性效应进一步增强，芯层初始压溃呈现出渐进压溃模式，当冲击波传播至远端时，各层的压溃和密实化取决于其静态压缩强度，密实化依次出现在强度较低的芯层。对于Case 2蜂窝，当C1层压缩量接近50%时，C1层不再压缩，C2与C3层发生耦合变形并开始被压缩，但准静态压缩强度较低的C3层优先被压缩，然后C2层被压缩，最终C1层被压缩，Case 2蜂窝整体达到密实化，这也使得Case 2的压缩强度由高到低再变高。另外，由图9还可以看出，采用式(8)计算得到的相对密度为5%的仿生蜂窝在不同冲击速度下的压缩强度（图9中蓝色直线）与模拟值（图9中绿色曲线）差异较小，同时对比两种梯度芯层在不同冲击速度下各层的压缩强度与理论值，表明各层的动态压缩强度可以由均匀仿生蜂窝动态压缩强度分别计算得出。

4.   结　论

基于王莲脉络模型，设计了一种新型仿生蜂窝结构，通过调节结构各部分的相对密度，对均匀蜂窝与面外梯度蜂窝的准静态压缩、动态压缩性能进行了研究，得出以下主要结论。

（1）利用ABAQUS有限元软件对新型王莲仿生蜂窝结构进行了准静态与动态压缩数值模拟，通过对均匀仿生蜂窝的静态压缩研究，得到了此仿生蜂窝准静态压缩强度与基体材料屈服应力之间的关系，基于现有的动态压缩理论，探讨了仿生蜂窝的动态压缩强度与相对密度、冲击速度之间的关系，结果表明，动态压缩强度的计算值与有限元模拟值之间的差异随着冲击速度的增大逐渐减小。

（2）面外梯度王莲仿生蜂窝结构在不同冲击速度下的动态压缩行为研究表明，低速冲击时（10 m/s），两种面外梯度仿生蜂窝的压缩强度曲线基本一致，表现为阶梯状，两者的变形模式类似，即准静态抗压强度较低的芯层依次被压缩。在高速冲击下（200 m/s），梯度芯层压溃模式与芯层梯度分布方式密切相关，Case 1蜂窝压缩强度曲线依旧保持为阶梯状渐进压缩模式，而Case 2蜂窝的压缩强度呈现“U”形变化，初始压溃呈现出渐进压溃模式，当冲击波传播至远端时，各层的压溃和密实化取决于其静态压缩强度，密实化依次出现在强度较低的芯层。梯度芯层每层的动态压缩强度可以由均匀仿生蜂窝动态压缩强度分别计算得出，此结果可为王莲仿生蜂窝结构的进一步优化设计提供参数依据。

感谢中国科学院西双版纳热带植物园园林园艺部吴福川副部长提供的王莲高清图片与有益的讨论。