1 引言
蜂窝结构依赖于自身的拓扑形式可呈现出不同的力学特性和泊松比效应，由于其具备出色的比强度、比刚度，且具有可设计性强的特点，在航空、航天等领域得到广泛运用，并吸引了国内外众多工程技术人员研究[1-2]。
零泊松比蜂窝结构的零泊松比特性主要表现为单方向的变形时，正交方向不产生形变[3]。Olympio等[4]首先提出了零泊松比六边形蜂窝结构，通过理论分析推导出等效杨氏模量和剪切模量的理论表达式，并分析蜂窝参数对其力学性能的影响。刘卫东[5-6]、孙秦[7-8]等人对零泊松比手风琴蜂窝结构的面内外等效模量进行了理论建模与参数分析。Chen[9]基于剪纸工艺研制了弯曲型零泊松比SILICOMB蜂窝样件，并对其刚度性能进行了研究。Gong[10]提出了一种可实现2个正交方向变形的四角星形零泊松比蜂窝结构，且在变形时能保证曲率光滑。Huang[11-12]提出了一种将传统六边形蜂窝结构与薄板连接组装在一起的新型零泊松比构型，并建立了结构的面内拉伸模量及等效弯曲模量的理论表达式。Broccolo等[13]将内凹六边形和正六边形的蜂窝结构交替布置，得到混合零泊松比蜂窝阵列结构。艾森[14]等以铝合金和钢材为基体材料，针对零泊松比蜂窝结构的非线性变形行为进行了研究。
波纹结构同样作为一种典型的具备极端各向异性的结构，表现为横向柔性和纵向刚性。Gong[15]将正弦型剖面波纹板视为薄壳结构，考虑波纹板的轴向拉伸变形和弯曲变形，采用卡式定理对横向拉伸、弯曲刚度进行了理论建模。Winkler等[16]采用广义平面应变有限元法，研究了不同几何形状对波纹结构的等效刚度的影响。Bartolozzi[17]采用卡式定理，提出了表征波纹结构性能的一般解析公式，可适用于不同的波纹形状。
综上所述，目前对于各种零泊松比蜂窝结构以及波纹结构的力学性能已进行广泛研究，但对具有曲面壁板的蜂窝结构弯曲性能研究较少。本文针对具有正弦波纹壁板的零泊松比类蜂窝结构，基于能量法对其等效弯曲刚度进行了理论推导与有限元仿真，分析了结构参数对其性能的影响。本研究获得的相关结果可以为该结构的优化设计、材料选型等提供重要支撑，亦可服务于柔性变形结构设计。
2 新型结构力学建模与仿真
2.1 结构介绍
新型零泊松比类蜂窝结构的几何结构如图 1所示。该结构典型特征为以正弦波纹结构作为蜂窝结构的倾斜壁板，其基本方程为：
(1)
式(1)中：H表示幅值；L表示波长(即壁长)。此外采用直壁连接各个独立胞元，以使得该结构在纵向伸缩时，可呈现出典型的零泊松比特性。值得注意的是当波纹幅值H为0时，该结构退化为手风琴式蜂窝结构，可用于变形蒙皮；而当波纹幅值不为0时，由于波纹壁板较出色的拉伸变形能力，将使得该新结构比手风琴式蜂窝结构易形成更大的变形，这有利于如伸缩翼变形飞行器等需要大变形的场景。
取其单胞元结构进行分析，采用与蜂窝结构类似的参数定义方法，胞元倾斜壁直线长度为L，直壁长度为a，蜂窝角为β，直壁壁厚为th，正弦波纹壁厚tl，定义壁高度为b。
图1 新型零泊松比几何结构示意图
Fig.1 Geometry of the novel zero Poisson’s ratio structure
2.2 理论模型
本文基于能量等效法建立新型零泊松比类蜂窝结构的等效弯曲刚度理论模型。结构等效示意图如图2。
图2 结构等效示意图
Fig.2 Schematic diagram of the structural equivalent
将图2所示的胞元结构等效成为正交各向异性薄板，根据薄板的变形能公式，可以得到其变形能为[18]：
(2)
式(2)中：D11、D22、D12为弯曲刚度；D66为扭转刚度。
本研究主要是对弯曲刚度进行分析，而在弯矩作用下，薄板内任意一点的挠度可以表示为：
w=w0+w1x2+w2y2
(3)
其中：
(4)
取板中心为约束点，即w0=0，联立可得：
(5)
对于单胞元结构的变形能，可假设单胞元受到均布弯矩载荷
和
的作用，根据载荷平衡条件，则可以得到各个胞壁受到的弯矩和扭矩为：
(6)
针对各个胞壁建立如图3所示的局部坐标系oξηζ，根据坐标转换关系可求解获得在局部坐标系下的受力情况为：
(7)
式(7)中：下标l表示为倾斜壁板；h为直壁板。
图3 正弦曲线波纹壁板示意图
Fig.3 Skech of Local coordinate system of sinusoidal corrugated wall
根据参考文献[19]，可获得弯矩引起的转角和扭矩引起的扭转角表达式为：
(8)
式(8)中：Es为胞元材料弹性模量；t为胞元壁厚度；l为胞元壁长度；w为胞元壁的高度；q为薄板扭转引入的修正系数。
波纹壁板受力图如图4。值得注意的是由于正弦曲线薄板在端面受到力矩作用时，一方面在各个截面会产生弯矩和扭矩分量，因此在各微元段会出现弯扭复合的现象，另一方面各微元段的弯矩与扭矩分量也会随着截面法向量的变化而变化。为修正正弦波纹壁板的扭转变形，本研究采用商业有限元软件ABAQUS对正弦波纹壁板施加扭矩时的变形能进行提取，对扭转修正系数进行反解。
图4 波纹壁板受力图
Fig.4 Moments acting on corrugated wall
在此基础上，对于胞元单个壁板的变形能，可表示为：
(9)
因此，整个单胞元的变形能为：
(10)
接下来分别对波纹壁板及直壁板的变形能进行分析。将波纹壁板离散为多个微元板，如图4所示可得各处的弯矩和扭矩分别为：
(11)
其中根据波纹几何方程(1)可得：
(12)
结合式(1)、式(9)、式(11)、式(12)，可得波纹板壁板变形能为：
(13)
其中：
(14)
结合式(7)、式(9)、式(10)、式(13)，可得整个单胞元的变形能为：
(15)
对比式(5)和式(14)，根据能量等效，可得：
(16)
其中D12=0，结合薄板理论中关系式：D12=μ2D11=μ1D22，可以推导出该结构泊松比μ1、μ2也为0。一方面说明了将结构等效成为正交各向异性薄板，采用能量等效原理求解等效刚度的方法不会影响结构的零泊松比特性；另一方面也说明本研究提出的新结构在受到纯弯矩Mx和My作用时，两正交方向的弯曲变形不存在耦合效应，这从弯曲变形的角度也验证了该结构的零泊松比特点。此外，可以看出D11与D22并不相等，且与外载荷无关，主要取决于结构与材料参数；其中D11主要取决于直壁板部分，与直壁板厚度比(th/L)以及壁高、蜂窝角直接相关，而D22则主要取决于波纹壁板，与波纹度(H/L)、波纹壁厚、壁高、蜂窝角直接相关。
3 有限元模型
本文采用Python实现参数化建模，使用商业有限元软件ABAQUS对结构的弯曲刚度进行模拟，有限元模型如图5。在不同规模下进行弯曲刚度的仿真模拟，边界条件如表1所示，经对比可发现，不同规模下的弯曲刚度具有较好的一致性。因此后续仿真均为针对不同结构参数下的单胞元结构，并将反解结果与理论模型进行对比。
图5 有限元模型示意图
Fig.5 FEM model of honeycomb structure
表1 有限元模型边界条件
Table 1 FEM model boundary conditions
边界D11(单胞)D22(单胞)D11(整体)D22(整体)A自由弯矩自由弯矩B自由固定自由弯矩C弯矩自由弯矩自由D固定自由弯矩自由O--固定固定
本研究中尚不考虑材料的非线性，且采用各向同性材料进行数值模拟，材料参数选用弹性模量E=200 GPa和泊松比υ=0.3。此外，考虑到壁厚相对较小，应用S8R壳单元进行模拟，仿真模型如图5所示。根据薄板理论，采用弯矩和挠度的比值表征等效弯曲刚度，具体计算式为：
(17)
式(17)中，wx与wy分别为B、C端面的挠度值。
图6展示了使用壳单元仿真的结果。从变形图6(a)中可以看出，结构沿Y方向拉伸时，X方向不会产生变形；此外，对比弯曲变形图6(b)与图6(c)、图6(d)，其中图6(b)为典型正泊松比蜂窝结构，可以看出其在面外变形时呈现“马鞍效应”，符合正泊松比结构特点，而本研究中新型结构在单方向弯曲时并不会引起正交方向的耦合变形，弯曲面较为平顺。这些特点都进一步验证了新结构的零泊松比特性，也说明了壳单元不会对结构的泊松比特性产生影响。
图6 有限元仿真结果示意图
Fig.6 Schematic diagram of FEM results
4 结果讨论与分析
4.1 与现有理论模型的对比
本文结构中当波纹度(H/L)为0时，结构可退化为手风琴式蜂窝结构，可采用文献[8]中的理论模型进行求解。选取w=10 mm，tl=0.1 mm，th=0.2 mm，a=10 mm，b=20 mm时，本文理论解、有限元解与文献结果如图7所示。从图7中可以看出，本文理论解、有限元解以及文献解三者在2个正交方向的弯曲刚度上均吻合较好，验证了本文模型的正确性以及对于手风琴式蜂窝结构的适用性。此外可看出D11和D22随蜂窝角的增大呈现相反的变化规律，另一方面从理论推导过程可以看出D11主要取决于直壁，而D22取决于倾斜壁板，这导致D11在量级上大于D22，呈现出典型的正交异性特征。
4.2 结构参数对D11的影响
根据理论推导，D11主要与蜂窝角(β)、直壁厚度比(th/L)、以及壁高(w)相关。因此该部分选取L=10 mm，主要针对这3个结构参数对D11的影响进行分析。图8—图10分别展示了它们对D11的影响，可以看出理论值与仿真值吻合度较高，最大误差不超过0.7%，进一步验证了理论解的正确性与适用性。
图7 蜂窝角β对弯曲刚度的影响曲线(H=0)
Fig.7 Variations of the bending stiffness with angle β(H=0)
图8展示了在H=2 mm、tl=0.1 mm、w=2 mm时，th分别为0.2 mm、0.4 mm、0.6 mm下，蜂窝角β由20°变化至60°时D11的变化规律。从图8可以看出，不同的th下，D11均随着蜂窝角的增大呈现下降趋势，且随着th的增大而增大。图9进一步描述了β=60°时，不同壁高下th对D11的影响规律，由图9可以看出，D11与直壁厚呈现线性正相关的趋势，与理论推导相符，且随着w的增大，增长速度也逐渐增大。
图8 不同直壁厚th下蜂窝角β对D11的影响曲线
Fig.8 The bending stiffness D11versus the cell angles β for various thickness th
图9 不同壁高w下直壁厚th对D11的影响曲线
Fig.9 The bending stiffness D11versus the thickness th for various height w
图10对比了在th=0.2 mm、tl=0.1 mm、 β=60°时，不同波纹壁幅值下壁高对D11的影响情况。对比发现D11对w的变化较为灵敏，与其三次方呈正相关，即在w增加2倍时，D11将增加8倍，这与理论推导结论一致。此外，可以看出D11与H无关，这是由于D11主要取决于直壁的贡献，而波纹斜壁的引入并不会降低D11，这有利于同时满足变形与承载功能。
图10 不同波纹幅值H下壁高w对D11的影响曲线
Fig.10 The bending stiffness D11versus the height w for variousamplitude H
4.3 结构参数对D22的影响
根据理论推导，D22主要与波纹度(H/L)、蜂窝角(β)、波纹壁厚(tl)、以及壁高(w)相关。该部分同样在L=10 mm条件下，主要针对这4个结构参数对D22的影响进行分析。
图11描述了在w=2 mm、tl=0.1 mm、th=0.2 mm时，弯曲刚度D22在不同的β下，随H变化呈现出的不同变化规律。对比有限元结果与理论解，可以发现均吻合较好，误差均小于5%。此外，在β较小时，可以看出弯曲刚度对H的变化不敏感；而当β较大为60°时，D22随着波纹幅值的增大呈现明显下降趋势。H=2 mm时，对比不同蜂窝角对D22的影响如图12所示。
图11 不同蜂窝角β下波纹幅值H对D22的影响曲线
Fig.11 The bending stiffness D22versus the amplitude H for various cell angles β
图12 不同波纹壁厚tl下蜂窝角β对D22的影响曲线
Fig.12 The bending stiffness D22versus the cell angles β for various thickness tl
从图12中可以看出，D22整体呈现出随β的增大而增大的趋势，这是由于β的增大使得波纹壁变形逐渐以弯曲变形为主，扭转变形则逐渐减小，而胞壁承弯能力较强，因此单胞的弯曲刚度增大；此外在tl较小时，D22的增长较缓慢。不同壁厚tl，不同壁高w下的D22结果如图13所示。
图13 不同壁高w下波纹壁厚tl对D22的影响曲线
Fig.13 The bending stiffness D22versus the thickness tl for various height w
从图13中可以看出，D22随tl的增大而增大，这是由于tl的增大使得倾斜壁板的截面惯性矩增大，进而使得单胞的弯曲刚度增大。此外可以发现在w、tl均较大时，有限元与理论解，会出现一定误差，这是由于在这种情况下，基于薄板理论的推导过程中未得到充分考虑的剪切变形所致。等效弯曲刚度D22针对w的参数化分析如图14所示，设置为tl为0.1 mm不变，从图14中可以发现，D22在不同的波纹幅值情况下均随w的增大而增大，且相比于手风琴蜂窝结构(H=0)，可以看出波纹壁的引入，使得D22有一定程度的下降，尤其是在w较大时，D22对波纹度较为敏感。
图14 不同波纹幅值H下壁高w对D22的影响曲线
Fig.14 The bending stiffness D22versus the height w for various amplitude H
5 结论
针对具有正弦波纹壁板的新型零泊松比类蜂窝结构的弯曲性能，基于能量法建立了等效刚度模型，进行了仿真验证与参数分析，理论分析结果与有限元结果吻合度较高，证明了理论解推导正确，为后续结构参数的优化设计奠定了基础。得出以下结论：
1) 该结构D11远大于D22，呈现典型的正交各向异性，且D11随蜂窝角的增大而减小，D22变化趋势与之相反。
2) D11与波纹度无关，且在波纹壁厚tl与蜂窝角β较小时，D22对波纹度的变化不敏感；而在β为60°时，D22随波纹度的增大显著下降。
3) D11与直壁厚th呈线性正相关，D22随波纹壁厚tl的增大逐渐增大，两者均随着壁高w的增大而显著提升。
[1] Qi C,Jiang F,Yang S.Advanced honeycomb designs for improving mechanical properties:A review[J].Composites Part B:Engineering,2021,227:109393.
[2] Wang Z G.Recent advances in novel metallic honeycomb structure[J].Composites Part B:Engineering,2019,166:791-741.
[3] 杨德庆,钟山.零泊松比超材料设计的多评价点功能基元拓扑优化方法[J].复合材料学报,2020,37(12):3229-3241.
Yang D Q,Zhong S.Functional element topology optimization method based on multiple evaluation points for metamaterial design with zero Poisson’s ratio[J].Acta Materiae Composite Sinica,2020,37(12):3229-3241.
[4] OlympioK R,Gandhi F.Zero Poisson’s ratio cellular honeycombs for flex skins undergoing one-dimensional morphing[J].Journal of Intelligent Material Systems and Structures,2010,21(17):1737-1753.
[5] 刘卫东,李红林.零泊松比手风琴蜂窝等效模量[J].固体力学学报.2018,39(01):100-112.
Liu W D,Li H L.Equivalent modulus of accordion honeycomb with zero Poisson’s ratio[J].Acta Mechanica Solida Sinica,2018,39(01):100-112.
[6] Liu W D,Zhu H,Zhou S,et al.In-plane corrugated cosine honeycomb for 1D morphing skin and its application on variable camber wing[J].Chinese Journal of Aeronautics.2013,26(04):935-942.
[7] 董文俊,孙秦.手风琴式蜂窝材料的等效弹性模量分析[J].机械科学与技术,2011,30(07):1103-1106.
Dong W J,Sun Q.Equivalent elastic modulus analysis of accordion cellular material[J].Mechanical Science and Technology for Aerospace Engineering,2011,30(07):1103-1106.
[8] 高珂,孙秦,董文俊.手风琴式蜂窝材料等效弯曲和扭转刚度分析[J].机械科学与技术,2014,33(10):1579-1584.
Gao K,Sun Q,Dong W J.Analysis of the equivalent flexural and torsional stiffness for accordion cellular honeycomb[J].MechanicalScience and Technology for Aerospace Engineering,2014,33(10):1579-1584.
[9] Chen Y,Scarpa F,Remillat C,et al.Curved Kirigami Silicomb cellular structures with zero Poisson’s ratio for large deformations and morphing[J].Journal of Intelligent Material Systems and Structures.2013,25(06):731-743.
[10] GongX B,Huang J,Scarpa F,et al.Zero Poisson’s ratio cellular structure for two-dimensional morphing applications[J].Composite Structures,2015,134:384-392.
[11] Huang J,Gong X,Zhang Q,et al.In-plane mechanics of a novel zero Poisson’s ratio honeycomb core[J].Composites Part B-Engineering,2016,89:67-76.
[12] Huang J,Zhang Q H,Scarpa F,et al.Bending and benchmark of zero Poisson’s ratio cellular structures[J].Composite Structures,2016,152:729-736.
[13] Broccolo SD,Laurenzi S,Scarpa F.AUXHEX-A Kirigami inspired zero poisson’sratio cellular structure[J].Composite Structures,2017,176:433-441.
[14] 艾森,郭瑜超,聂小华,等.零泊松比蜂窝结构一维变形行为[J].南京航空航天大学学报,2021,53(04):629-636.
Ai S,Guo Y C,Nie X H,et al.One-dimensional deformation behavior of a honeycomb structure with zero Poisson’s ratio[J].Journal of Nanjing University of Aeronautics & Astronautics,2021,53(04):629-636.
[15] Gong X B,Liu L W,Scarpa F,et al.Variable stiffness corrugated composite structure with shape memory polymer for morphing skin applications[J].Smart Materials and Structures,2016,26(03):035052.
[16] Winkler M,Kress G.Influence of corrugation geometry on the substitute stiffness matrix of corrugated laminates[J].CompositeStructures,2012,94(09):2827-2833.
[17] Bartolozzi G,Baldanzini N,Pierini M.Equivalent properties for corrugated cores of sandwich structures:A general analytical method[J].Composite Structures,2014,108:734-746.
[18] 郭瑜超,王立凯,段世慧.复杂多边形蜂窝结构等效弯曲刚度研究[J].机械强度,2017,39(01):122-129.
Guo Y C,Wang L K,Duan S H.A study on equivalent bending stiffness of complex polygonal honeycomb structure[J].Journal of Mechanical Strength,2017,39(01):122-129.
[19] ChenD H.Equivalent flexural and torsional rigidity ofhexagonal honeycomb[J].Composite Structures,2011,93(07):1910-1917.