崔飞蝶，董焱章，2，李青峰
（1.湖北汽车工业学院 汽车工程学院，湖北 十堰 442002；
2.汽车动力传动与电子控制湖北省重点实验室，湖北 十堰 442002）
负泊松比材料具有较大的屈服强度、剪切模量、断裂韧性和较低的疲劳裂纹扩展，对冲击能量的吸收能力远大于传统材料，因此负泊松比结构广泛应用于冲击、碰撞和爆炸等领域［1-5］。Sigmund等［6］首次运用逆均匀化方法求解拓扑优化问题，得到了具有负泊松比特性的多孔结构；
Zhang 等［7］运用基于能量均匀化方法得到了具有极限负泊松比特性的微结构；
赵颖等［8］采用遗传算法设计出新型内凹三角形负泊松比微结构B柱，相较于无夹芯碳纤维B 柱与碳纤维B 柱，其侧向弯曲位移显著降低；
刘宇等［9-10］通过拓扑优化方法得到较为完整的负泊松比基元优化构型，相比传统凹角六边形结构，具有更好的吸能效果；
孙晓旺等［11］结合多目标优化方法将内凹六边形负泊松比蜂窝材料作为车辆底部防护组件的夹芯，有效降低了防护组件的挠度和动能；
杜义贤等［12］构建了以负泊松比和刚度为优化目标的多目标拓扑优化模型，比较了不同等效刚度的负泊松比结构的性能；
张会凯［13］使用拓扑优化方法设计了新的多材料手性拉涨（负泊松比）超材料，同时具备扭转特性。现有文献大多借助优化方法设计负泊松比材料单个微结构，然后进行一定规则的布置得到负泊松比吸能宏观结构。文中提出了多型负泊松比材料微结构设计理念，从微观材料和宏观结构层面对负泊松比吸能结构进行优化设计，分别建立负泊松比的材料吸能微结构和承载宏观结构拓扑优化模型，并实现材料-结构协同优化设计，最终对优化结果进行冲击性能分析。
负泊松比材料拥有优异的冲击吸能性能，将拓扑优化后的负泊松比材料微结构作为单胞，经过阵列得到单一型负泊松比微结构（single-microstructure negative poisson’s ratio，SNPR）蜂窝吸能材料，但负泊松比材料的刚度往往较低，为了提高负泊松比材料的刚度性能，同时保证冲击吸能特性，引入材料-结构协同优化设计替代原有的周期性布置方案，探究由多型负泊松比微结构（multi-microstructure negative poisson’s ratio，MNPR）组成的负泊松比吸能宏观结构的冲击性能。整体研究分为3个步骤，首先建立负泊松比材料微结构拓扑优化模型，以微结构泊松比最小为优化目标，微观材料设计区域各单元的相对密度为设计变量，优化得到在多种体分比约束下对应的最佳负泊松比微结构材料分布，通过特征提取获得微结构优化构型；
其次对吸能宏观结构进行刚度优化，以结构宏观刚度最大为优化目标，宏观设计区域各单元的相对密度为设计变量，优化得到在特定体分比约束下的宏观结构最佳材料密度分布；
最后以宏观结构优化结果的材料分布情况来指导各负泊松比材料微结构间的排列组合，建立包含多型微结构的负泊松比承载吸能材料模型。如图1所示。
1.1 负泊松比材料微结构拓扑优化模型
结构泊松比表达式为
式中：εx为材料在弹性范围内纵向加载时产生的横向应变；
εy为纵向应变。若直接以式（1）为目标函数，求解具有负泊松比特性的拓扑构型时存在高度非线性问题，故采用线性表达式近似拟合结构的泊松比［13］，考虑数值的稳定性，引入Q2222，等效数学函数表示为
式中：Q1122为纵向拉伸时横向产生的能量；
Q1111为纵向拉伸时纵向产生的能量；
Q2222为横向拉伸时横向产生的能量；
β取0.8；
l为优化迭代次数。
负泊松比微结构拓扑优化目标为结构泊松比最小，将微观设计域离散为100×100 个设计单元，以划分单元密度为设计变量，在指定体分比的约束条件下，优化得到具有最小泊松比的微结构，拓扑优化优化列式为
式中：c1为结构的等效泊松比；
f为全局体积约束；
U和F分别为宏观尺度下的整体位移矢量和全局载荷矢量；
K为整体刚度矩阵；
xe为微观划分单元的相对密度，取值0～1；
N为设计域中单元的总数；
V为单元体积；
V0为设计域体积。为避免迭代过程中矩阵奇异，取xmin为0.001。
1.2 负泊松比宏观承载结构拓扑优化模型
结构刚度越大，承载能力越强，因此在宏观层面以结构刚度最大或柔度最小为优化目标，设计区域各单元的相对密度为设计变量，设定体分比约束，优化得到具有最优承载能力的宏观材料分布。宏观结构拓扑优化列式为
式中：c2为结构的柔度；
ui为单元的位移矢量；
k0为单元刚度矩阵；
f取0.4；
xi为宏观划分单元相对密度；
N为设计域中单元的总数；
a和b为相对密度的下限和上限，分别取0.3和0.5。
1.3 负泊松比材料微结构的边界连接设置
进行负泊松比材料微结构拓扑优化时，通过设置不同的体分比约束，可得到不同体分比下的最优微结构。虽然优化过程中施加了周期性边界条件，但考虑到相邻微结构之间的连通性，在微观设计域内设置部分非设计区域便于后期连接，如图2 所示。非设计域分布于设计域的4个边界。
2.1 负泊松比材料微结构拓扑优化设计
依据上述研究方法和材料-结构拓扑优化模型，基于能量法和优化理论，设置材料相对弹性模量E0为1，泊松比μ为0.3，微观设计域划分为100×100 个4 节点离散单元，对初始构型进行中间挖孔来诱导优化加快形成负泊松比结构。在不同体分比约束条件下，优化得到的部分负泊松比微结构最优拓扑构型如图3 所示。当体分比在［0.30～0.37）区间时，最终优化结果均为同一种构型，即体分比为0.3 时所对应的优化拓扑构型，其他区间也表现出类似情况。考虑到对拓扑结果进行特征提取时微小差异可以被简化，选取同一拓扑优化结果下的代表性体分比来代替对应的体分比区间，最终选取3 个代表性体分比来对应在体分比约束区间为0.3～0.5 优化得到的3 种拓扑微结构。图3 中黑色部分表示有材料，白色部分表示无材料。
图3 不同体分比下负泊松比材料微结构最优拓扑构型
2.2 负泊松比宏观承载结构拓扑优化设计
宏观承载吸能结构简化为矩形防冲击结构，顶部施加均布载荷，底部固定。优化程序中，设置体分比约束为0.4，设计域离散为15×15 个4 节点单元，设计区域各单元材料相对密度上下限设为0.3和0.5，以结构柔度最小为优化目标，优化结果见图4a，白色区域相对密度为0.3，黑色区域相对密度为0.5，灰色区域相对密度为0.3～0.5，按照宏观划分单元密度值与负泊松比微结构体分比数值相同的原则，将相应的负泊松比材料微结构优化构型进行排列组合，最终得到由多型负泊松比微结构组成的冲击防撞结构最优拓扑构型，如图4b所示。
图4 矩形冲击防撞结构最优拓扑构型
3.1 有限元模型建立与抗冲击评估指标
采用静力学与动力学仿真来模拟试样拉胀过程，测试MNRP 材料模型的静态承载能力，并与同体分比下的SNRP 蜂窝材料进行对比，探究MNRP材料模型的抗冲击性能。
有限元模型使用Abaqus 建立，包含中间的负泊松比蜂窝材料和位于蜂窝上下端的2个刚性板。负泊松比模型采用金属材料，用理想弹塑性模型材料预测负泊松比材料的塑性变形，材料密度为2700 kg·m-3，杨氏模量为68 GPa，屈服应力为255 MPa，材料泊松比为0.3；
上下刚性板杨氏模量为210 GPa，材料密度为7800 kg·m-3。负泊松比材料和上下刚性板的有限元网格模型采用4节点线性减缩积分传统壳单元（S4R），综合模型计算时间和精度，设置单元网格尺寸为0.5 mm，其中MNRP模型的单元数量为629 680 个，SNRP 模型的单元数量为548 400 个，厚度方向设置5 个积分点以确保计算收敛和结果的准确性。
体分比为0.4的多型微结构负泊松比材料和单一型微结构负泊松比材料有限元模型如图5所示，蜂窝尺寸为150 mm×150 mm，面内厚度为10 mm。仿真过程中，仅赋予上端刚性板Y方向的移动自由度，并设置恒定的压溃速度（准静态压缩）；
下端刚性板被完全约束；
约束蜂窝Z方向的移动自动度和绕X轴、Y轴的转动自由度。采用一般接触算法定义模型的全局相互作用，包括摩擦系数为0.2的切向摩擦相互作用和用来避免相邻部件之间相互渗透的法向硬接触。
图5 MNPR和SNPR结构的有限元模型
负泊松比材料的抗冲击防撞性能评价指标主要包含模型总吸能量、比质量吸能、碰撞峰值力以及平台力等。
3.2 准静态压缩有限元仿真分析
设置准静态压缩速度为10 mm·min-1，对比分析MNRP 材料和SNRP 材料的等效弹性模量大小，只考虑材料在准静态压缩的弹性阶段。图6是2种结构的准静态压缩应力-应变曲线对比图，可以看出，2种结构在准静态压缩过程中均有应力增强现象；
在相同应变下，MNPR 的应力高于SNPR，曲线斜率近似看作结构的等效弹性模量，因此MNPR的等效弹性模量高于SNPR，意味着MNPR 有更好的静态以及小载荷承载能力。
图6 MNPR和SNPR准静态压缩应力-应变曲线
3.3 动态冲击有限元仿真分析
对MNRP 材料和SNRP 材料进行冲击仿真分析，冲击速度设为10 m·s-1、50 m·s-1以及100 m·s-1，对比分析2种结构的抗冲击吸能性能。
图7～9给出了MNPR和SNPR在不同冲击速度下的变形模式。从图7 可以看出：低速冲击作用时，结构的变形主要集中在冲击端和固定端；
随着结构变形的增加，两端被压溃的层数逐渐增多，同时MNPR 内部也逐渐表现出从两侧向内收缩的变形模式；
随着结构变形的进一步增加，MNPR 中间层结构开始主导整个结构的变形，呈现出类似X型变形模式，展现了结构的负泊松比效应，而SNPR依旧表现为从两端依次向中间溃缩的变形模式，中间层结构相比于结构两端的变形量较小。从图8可以看出：随着冲击速度的增大，在结构变形初始阶段，2 种结构的变形均主要集中在冲击端，同时伴随着固定端小范围变形模式；
随着变形的增加，2种结构均呈现“一”字型逐层溃缩变形模式，同时2种结构的固定端变形进一步小范围增加，且SNPR固定端变形量率多于MNPR；
当应变逐渐增加到0.6时，MNPR中间结构变形逐渐明显，左右两侧向中间靠拢，且呈现中间快两侧慢的变形趋势，SNPR仍呈现“一”字型逐层溃缩变形模式，同时2种结构底部固定端压溃量也逐渐增多。从图9可以看出：当冲击速度增加到100 m·s-1时，惯性效应进一步增强，2 种结构的变形带基本集中在冲击端，均呈现出从冲击端到固定端“一”字型逐层压溃变形模式；
在压缩过程中，MNPR较SNPR中间部位略有不同，左右两侧向内收缩的现象更加明显。综上，在冲击载荷作用下，2 种结构均表现出负泊松比现象，即受压结构两侧向内收缩；
此外冲击速度对蜂窝材料的变形模式有很大影响，冲击速度越高，结构变形模式越趋向于由冲击端到固定端的逐层溃缩；
并且MNPR 相对于SNPR 中间部位的收缩更加明显和复杂。
图7 冲击速度10 m·s-1时MNPR和SNPR变形模式
图8 冲击速度50 m·s-1时MNPR和SNPR变形模式
图9 冲击速度100 m·s-1 时MNPR和SNPR变形模式
不同冲击速度下2种结构的应力-应变曲线如图10 所示。在冲击开始阶段，与一般多胞材料类似，MNRP 材料和SNRP 材料的动态响应均近似于线性，短时间内应力水平达到较高峰值，随着材料被逐渐压缩，应力进入较平稳的平台区，直到结构被压缩至密实阶段，应力迅速增加，材料吸能结束。从图10 中可以看出：冲击速度对材料应力水平的影响较大，随着冲击速度的增加，峰值应力和平台区应力均增加；
相同冲击速度下2种结构的应力水平不同，冲击速度为10 m·s-1时MNPR 和SNPR 峰值应力近似相等，而平台阶段MNPR 应力大于SNPR（图10a）；
冲击速度为50 m·s-1时，SNPR峰值应力略高于MNPR，而平台阶段MNPR应力仍大于SNPR（图10b）；
当速度增加到100 m·s-1时，SNPR达到峰值应力晚于MNPR，且峰值应力高于MNPR（图10c），而平台阶段SNPR应力水平依旧低于MNPR。综上，3种冲击速度下MNPR峰值应力均低于或近似等于SNPR，而较低的峰值应力意味着被保护结构在承受冲击时受到更小的冲击载荷，有利于被保护结构的完好性，同时MNPR 展现出比SNPR更高的平台应力，有利于结构对冲击能量的吸收，因此多型微结构负泊松比材料的抗冲击性能优于单一型微结构负泊松比材料。
图10 不同冲击速度下MNPR和SNPR的应力-应变曲线
在不同冲击速度下2种结构的能量吸收-时间曲线如图11 所示，可以看出：3 种冲击速度下MNPR的能量吸收均高于SNPR；
由于2种结构的体分比都为0.4，所以同等质量下，MNPR的能量吸收高于SNPR，因此MNPR 的抗冲击吸能能力优于SNPR；
冲击速度越大，材料防护时效越短，冲击过程中材料能量吸收量越多。
图11 不同冲击速度下MNPR和SNPR的能量吸收-时间曲线对比
表1列出了2种负泊松比材料模型在3种冲击速度下产生相同应变时对应的吸能量大小。冲击速度较低时，MNPR 的能量吸收是SNPR 的2 倍以上；
冲击速度为50 m·s-1时，MNPR的能量吸收接近SNPR 的2 倍；
冲击速度为100 m·s-1时，小应变下MNPR 的能量吸收量是SNPR 的164%，ε为0.4 时MNPR的能量吸收量是SNPR的177%。
表1 相同应变下2种模型的吸能量 kJ
文中系统分析了基于材料-结构协同优化获得的多型微结构负泊松比承载吸能材料在准静态压缩以及面内冲击载荷作用下的动力学性能，并与同等体分比下的单一型微结构均匀负泊松比材料进行了对比分析。准静态压缩仿真结果表明，多型微结构负泊松比材料有更好的静态承载能力；
面内冲击仿真结果表明，冲击速度对材料应力水平的影响较大，随着冲击速度的增大，2 种材料的固定端应力-应变曲线中包括峰值应力和平台区应力均增大；
同时基于多型微结构的负泊松比材料的抗冲击能力优于同等体分比下的单一型微结构均匀负泊松比材料，表现为峰值碰撞力低、平台区应力高等；
此外，多型微结构负泊松比材料有更高的能量吸收能力，同等应变下多型微结构负泊松比材料的能量吸收是单一型微结构负泊松比材料的1.6倍以上。从2种材料的冲击变形模式来看，多型微结构负泊松比材料在纵向受压时横向收缩幅度更大，与单一型微结构均匀负泊松比材料的受冲击逐层压溃变形模式相比，更加复杂和集中，归因于结构的宏观刚度优化使材料分布更加合理，从而冲击吸能量和平台区应力等显著高于单一型微结构负泊松比材料，证实了基于多型微结构负泊松比材料-结构协同优化设计方法的可行性和有效性，为设计具有轻质高承载和良好抗冲击性能的吸能防护材料提供了参考。