蜂窝是一种由平面上单个胞元进行某种特定排布拉伸形成的多孔结构。具有重量轻、比刚度高、比强度高、吸能能力强、孔隙率高等优点[1]，在吸能结构[2–4]、热防护结构[5]、医学[6]、化学工程[7]及飞机结构[8–9]等领域得到了广泛的应用。
变体飞行器能够根据飞行工况改变自己的形状，从而改变气动性能，在整个飞行包线中获得更加优异的操纵能力或者更加经济的油耗[10]。蜂窝结构可以通过蜂窝壁的弯曲变形来实现平面内的大变形[11]，此结构受到众多变体领域研究者的关注。随着超材料技术的发展[12–13]，蜂窝的性能得到了极大的开发，零泊松比蜂窝、负泊松比蜂窝等蜂窝结构相继被提出，并在变体飞行器领域得到运用。Budarapu 等[14]提出了一种设计飞机机翼结构的框架，并分析了具有手性结构的变形翼型，该翼型的特点是能够在有限的应变下产生较大的位移变形，且泊松比接近–1。Bornengo等[15]提出了一种负泊松比为–1 的六边形手性蜂窝结构，并研究了其面内拉伸性能。Airoldi 等[16–17]制造了由复合材料薄层板制成的六韧带手性蜂窝，并对其拉伸性能进行了有限元分析及试验对比。Gong 等[18]提出一种用于二维机翼变形的零泊松比蜂窝并进行了有限元仿真，分析了其面外横向刚度。Huang 等[19]提出一种用于折叠翼梢变形的新型零泊松比蜂窝，并研究了其等效弯曲刚度及弯曲强度。对于飞机进气道、雷达整流罩等采用的鼓包构型，负泊松比蜂窝的面外变形由于其负泊松比效应，呈现出双曲面 （与鼓包形状相似），可以很好地适应鼓包变形，在这些构型的变体结构中有不错的应用前景。但为了使蜂窝变形后与鼓包贴合，通常需要对蜂窝弯曲刚度分布进行设计优化，其中蜂窝结构面外变形需要进行大量求解计算，仅靠有限元方法求解往往因时间成本高而不具有可行性，因此需要一种能够快速准确对蜂窝面外变形进行求解计算的方法。
本文针对四韧带反手性 （反四手性）负泊松比蜂窝面外变形求解的问题，首先利用有限元分析软件对均匀蜂窝的面外变形进行求解，使用最小二乘法对求解结果拟合，得到蜂窝对称轴上的拟合曲线；然后对求解结果以及蜂窝尺寸参数进行二次拟合，得到蜂窝尺寸和变形曲线的关系，利用该关系对其他尺寸下蜂窝的变形曲线进行预测，并与仿真值进行对比；最后对非均匀的蜂窝，利用之前的拟合结果构建一个新的拟合函数，求得拟合结果，并把拟合结果与仿真结果进行了对比。
1 均匀反四手性蜂窝面外变形数值方法分析
如图1 所示，反四手性蜂窝的胞元单位由1 个圆及4 条肋构成，L 为水平方向的肋长；L*为竖直方向肋长。对于L = L*的胞元，4 条肋关于圆心中心对称；而对于L≠L*的胞元，1、3 方向的肋关于圆心中心对称，2、4 方向的肋关于圆心中心对称。整个蜂窝结构有两条对称轴，把1、3 方向的对称轴设为x 轴，2、4 方向的对称轴设为y 轴。
图1 反四手性蜂窝示意图
Fig.1 Diagram of anti-tetrachiral honeycomb
1.1 面外变形挠度仿真计算
运用ABAQUS 有限元软件计算反四手性蜂窝在固定载荷下的变形。引入变形系数k = L*/L 表征蜂窝的形状特征，以k = 0.6 为例介绍有限元仿真过程。采用mm–ton–N 量纲系统，假设蜂板四周固定，受到指向z 轴的均布压力。胞元半径r = 10 mm，壁厚t = 1.0 mm，y 方向固定肋长L = 50 mm，x 方向肋长L*为变化值（本文中变化范围为[30 mm，70 mm]），蜂窝结构的宏观尺寸为1000 mm×500 mm×20 mm。由于y 方向胞元肋长固定，因此y 方向胞元数固定，在蜂窝边界处，为了不出现胞元被裁剪的情况，在边界处添加了一定长度的肋长使整个蜂窝结构尺寸完整。对于x 方向，由于该方向蜂窝肋长不固定，因此该方向上的胞元数及边界肋长是一个随k 变化的数值，具体取值见表1。蜂窝的物理参数设为 （200 MPa，0.35）。在蜂窝上端添加一个相同尺寸的薄膜用于加载压强载荷，其弹性模量设为蜂窝的千分之一，薄膜厚度为1.0 mm。蜂窝四周固定，薄膜四周铰接，蜂窝与薄膜之间采用绑定 （TIE）约束。薄膜受到一个指向蜂窝的合力，大小为1000 N 的压强。
表1 不同k 值下蜂窝的尺寸参数
Table 1 Dimensional parameters of honeycomb with different k values
kx 方向胞元数/个x 方向边界长/mmy 方向胞元数/个y 方向边界长/mm 0.633101015 0.72817.51015 0.825101015 0.92217.51015 1.020151015 1.11822.51015 1.217101015 1.3162.51015 1.414351015
使用静态隐式分析步进行计算，开启几何非线性选项，初始增量0.01，最小增量1E–7，最大增量0.1，最大增量步2000。在蜂窝两个方向的对称轴上分别取点作为位移挠度输出，如图2 所示，输出量为COOR和U。使用S4R 壳单元对结构进行网络划分，蜂窝网格种子大小为10，薄膜种子大小为50 mm。蜂窝结构包含10540 个S4R 单元、 13863 个节点；薄膜结构包含200 个S4R 壳单元、 231 个节点。其余参数采用默认设置。
图2 选取的位移输出点
Fig.2 Selected displacement output points
由于载荷的对称性，蜂窝的挠度分别在x 轴、y 轴两条对称轴上取得最大值。因此以这两条对称轴的挠曲线作为典型特征研究蜂窝的变形。提取出x 轴和y 轴在z 轴方向上的挠度。图3 为不同k 值下蜂窝沿x方向和y 方向对称轴的变形挠曲线，从图3（a）可以看出，挠曲线呈现出一个开口朝下的单峰曲线，曲线关于x = 0 对称且与x 轴相交于固定约束点 （–500 mm，0）、 （500 mm，0）。随着k 值的增大，蜂窝的变形也增大。当k 较大时，沿着x 方向，蜂窝的挠度先增大，接近最大值后增长趋势变缓，形成一个类平台形状，维持一段距离后降低到0；当k >1.2 时，挠曲线的光滑性降低，这是由于当k 值增大时，蜂窝单胞尺寸也随之增大，胞元数量降低，采样节点减少，连续性降低。图3（b）为不同k 值下蜂窝沿y 方向对称轴的变形挠曲线，可以看出y 方向与x 方向的挠曲线特征基本相同，曲线关于y = 0 对称且与x轴相交于固定约束点 （–250 mm，0）、（250 mm，0）。主要区别在于y 方向在挠曲线最大值处没有形成类似台阶的形状，这是因为蜂窝结构在两个方向上的尺寸不同，x 方向尺寸大于y 方向尺寸，因此x 轴方向的挠曲线相比y 轴方向有一个更平缓的峰值。
图3 不同k 值下沿2 个方向蜂窝对称轴的挠曲线
Fig.3 Deflection curves along two-direction honeycomb symmetry axis at different k values
1.2 面外变形挠度拟合
在得到蜂窝受力变形的仿真结果后，提取出x 轴和y 轴在z 轴方向的挠度，以k = [0.6，0.8，1.0，1.2，1.4]时的结果为拟合数据，拟合出一条曲线，得到挠度与坐标的关系，再次拟合得到系数与蜂窝尺寸的关系。
建立坐标系，使坐标系原点与蜂窝中点重合，x 轴、y 轴、z 轴分别与蜂窝x 轴、y 轴、z 轴重合。从图3（a）中可以看出，不同k 值下的挠曲线w有边界条件 （w）x = ± 500 = 0，且曲线关于x = 0 （y 轴）对称，设x 方向的拟合方程形式为w （x）= a1（x2–5002）3 + b1（x2–5002）2 + c1（x2–5002），其中包括a1、b1、c1 3 个待定系数。拟合结果如表2 所示。在保证拟合精度 （R2>0.99）的情况下，对系数的值进行了微调。
表2 不同k 值下x 方向挠曲线拟合结果
Table 2 Fitting results of the deflection in x-direction with different k values
k 方程形式a1/E–15b1/E–9c1/E–5R2 4.51.5–2.640.9988 0.85.01.7–3.050.9991 1.05.51.9–3.440.9967 1.26.02.1–3.720.9969 1.46.52.3–3.870.9963 0.6 w（x）= a1（x2–5002）3+ b1（x2–5002）2+ c1（x2–5002）
从表2 中可以看出，随着k 值的增大，a1 和b1 的值呈现出单调递增的态势，而c1 呈现出单调递减的趋势。其中，a1 和b1 的线性程度非常明显，c1 的线性程度稍弱。此外，虽然a1 的值的数量级达到了E–15，但是这并不意味a1 值对应的项a1×（x2–5002）3 可以忽略。为了兼顾方程的物理意义，拟合时并没有对挠度数据进行中心化和缩放，事实上，在k 值为0.8 情况下，当x = 0 时，a1×（x2–5002）3 项对应的值为–78.125，而此时整个方程w（0）=35.75，显然不能忽略。b1、c1 值同理。
图4 为不同k 值下沿x 方向蜂窝挠曲线拟合情况，可以看出，当k值较小时，拟合曲线与仿真结果吻合度较高，而随着k 值的增大，仿真结果波动也越来越明显，这是由于k 值变大以后，蜂窝胞元的尺寸也增大，在宏观上看蜂窝结构的均匀度降低，其物理参数的波动也更加剧烈。挠曲线的系数a1、b1、c1 以k 为自变量进行拟合，结果如表3 所示。
表3 x 方向挠曲线参数拟合结果
Table 3 Parameter fitting results of deflection curve in x-direction
参数拟合方程R2 a1a1（k）=（2.5k + 3）E–151 b1b1（k）=（k + 0.9）E–91 c1c1（k）=（1.1k2 – 3.765k – 0.767）E–50.9987
图4 不同k 值下沿x 方向蜂窝对称轴的挠曲线拟合情况
Fig.4 Fitting curves of the deflection in x-direction at different k values
同理，确定y 方向的拟合方程形式为w（y）= a2（y2–2502）3 + b2（y2–2502）2+ c2 （y2–2502），a2、b2、c2 为待定系数。对仿真结果进行最小二乘法拟合，得到拟合结果如表4 及图5所示。从表4 可以看出，对于y 方向，拟合方程的三次项a2×（y2–2502）3 为0，随着k 值的增大，b2 的值也呈现出单调递增的态势，而c2 的值呈现出单调递减的趋势。
表4 不同k 值下y 方向挠曲线拟合结果
Table 4 Fitting results of the deflection in y-direction at different k values
k 方程形式a2/E–15b2/E–9c2/E–5R2 0.6 w（y）=a2（y2–2502）3+b2（y2–2502）2+c2（y2–2502）0 7.2–3.010.9994 0.808.2–6.840.9991 1.009.2–10.420.9991 1.2010.2–12.820.9994 1.4011.2–14.710.9993
表5 y 方向挠曲线参数拟合结果
Table 5 Parameter fitting results of deflection curve in y-direction
参数拟合方程R2 a2a2（k）= 01 b2b2（k）=（5k + 4.2）E–91 c2c2（k）=（10k2 – 34.69k + 14.33）E–50.9993
图5 不同k 值下沿y 方向蜂窝对称轴的挠曲线拟合情况
Fig.5 Fitting curves of the deflection in y-direction at different k
values
比较x 方向和y 方向的拟合结果 （表3 和5），最明显的差别是x 方向的挠曲线有 （x2–5002）3 项，而y 方向的挠曲线没有对应的项，这是因为x 方向蜂窝的尺寸比y 方向大，在挠曲线中表现为在最大值处有一个台阶，而y 方向的挠曲线顶部较为尖锐，因此可以认为， （x2–5002）3项与蜂窝的宏观尺寸相关。而对于（x2–5002）2 项和 （y2–2502）2 项，两个方向的b 值拟合曲线均为关于k 的一次多项式，这互相印证了拟合结果的有效性。系数c 的情况与系数b 类似，只是拟合结果是关于k 的二次多项式。此外，y 方向的挠曲线拟合效果优于x 方向，这是因为当k值较大时，x 方向仿真结果的波动更加明显。
1.3 拟合结果预测情况
利用k = [0.7，0.9，1.1，1.3] 时的仿真结果验证，根据1.1 节中的拟合结果，当k 值分别取0.7、0.9、1.1、1.3时，拟合曲线仿真如表6 所示。
表6 不同k 值下蜂窝挠度的预测曲线
Table 6 Prediction curves of honeycomb deflection at different k values
k x 轴挠度y 轴挠度0.7w（x）= 4.75E–15×（x2–5002）3+1.60E–9×（x2–5002）2–2.80E–5×（x2–5002）w（y）=7.5E– 9×（y2–2502）2– 5.05E– 5×（y2–2502）0.9w（x）=5.25E–15×（x2–5002）3+1.80E–9×（x2–5002）2–3.18E–5×（x2–5002）w（y）=8.5E– 9×（y2–2502）2– 8.79E– 5×（y2–2502）1.1w（x）=5.75E–15×（x2–5002）3+2.0E–9×（x2–5002）2–3.49E–5×（x2–5002）w（y）=9.5E– 9×（y2–2502）2– 11.73E– 5×（y2–2502）1.3w（x）=6.25E–15×（x2–5002）3+2.2E–9×（x2–5002）2–3.69E–5×（x2–5002）w（y）=10.5E– 9×（y2–2502）2– 13.87E– 5×（y2–2502）
图6 为两个方向下拟合预测曲线与仿真结果的对比。可以看出，预测值与仿真值有较好的一致性，这说明拟合结果是有效的。当k 值较小时，二者的一致性更好，这是因为k值较大时，仿真结果出现了一定程度的偏差，该偏差在x 方向更为明显。
图6 x 方向和y 方向拟合预测曲线与仿真值对比
Fig.6 Comparison of x-direction and y-direction prediction curves and simulation values
2 非均匀反四手性蜂窝面外变形性能分析
本节对非均匀的反四手性蜂窝结构整体面外变形进行分析，根据第1 节的理论方法，采用数值分析方法，在不改变y 轴肋长变量的同时，将x 轴方向肋长分为两段,如图7 所示。每一段的肋长取不同值，分析其变形曲线。在–500 mm≤x≤0范围内，x 方向肋长L1 是固定值50 mm，在0≤x≤500 mm 部分肋长L2是变化值，变化范围为[30 mm，70mm]，同样引入变形系数k* = L2/L1 表征蜂窝的形状特征。由于沿y 方向挠度在本节不是主要关注对象，因此忽略不做分析。
图7 分段蜂窝示意图（mm）
Fig.7 Diagram of segmented honeycomb (mm)
2.1 仿真结果分析
与第1 节类似，对长1000 mm、宽500 mm、厚度20 mm 的蜂窝薄板进行力学仿真分析，仿真过程与第1节相似。对于不同的k*值，蜂窝结构的尺寸参数如表7 所示。
表7 不同k*值下蜂窝的尺寸参数
Table 7 Dimensional parameters of honeycomb with different k* values
k*x 方向胞元数/个x 方向边界长/mmy 方向胞元数/个y 方向边界长/mm 0.626151015 0.72412.51015 0.822201015 0.921151015 1.020151015 1.119201015 1.218301015 1.31812.51015 1.417301015
仿真结果如图8 所示。可以看出，挠曲线呈现出双峰开口向下的曲线特征。不同k*值下，对蜂窝固定肋长部分的挠曲线影响较小，主要还是影响变肋长部分的曲线。同样，k*值较大时，仿真数据抖动比较明显，可以发现k*=1.3 时挠度的最大值甚至超过了k*=1.4 时的最大挠度。
图8 非均匀情况下不同k*值沿x 方向对称轴的蜂窝挠度
Fig.8 Honeycomb deflection along the x-direction symmetry axis at different k* values in the non-uniform case
2.2 面外变形挠度拟合
用uk（x）表示第1 节中求得的沿x 方向挠曲线方程，合理猜测，本节2 段非均匀蜂窝挠曲线方程由2段各自k 值对应的挠曲线方程加权求和得到，即
式中，g1（x）和g2（x）为uk（x）的权重函数，根据对称性其需满足
引入Logistic 函数，其基本形式为
该函数图像如图9 所示，显然，该函数满足式 （2）中要求。
图9 Logistic 函数图像
Fig.9 Graphics of Logistic function
设
式中，m 为拟合变量，则
此方程为非均匀蜂窝挠曲线拟合方程，该方程仅有m 一个未知数。
运用最小二乘法，对式 （5）进行拟合，得到m 的拟合值为– 0.016（R2 = 0.9949），其拟合曲线与仿真结果如图10 所示。可以看出，拟合结果与仿真曲线相符程度较好，证明对w（x）的猜测是正确的。
图10 k*= 0.6 时蜂窝沿x 方向对称轴挠度的仿真结果与拟合曲线
Fig.10 Simulation results and fitted curves of honeycomb deflection along the x-direction symmetry axis for k*=0.6
同理，对k*= 0.6、0.8、1.0、1.2、1.4对应的w（x）分别进行拟合，得到结果如表8 所示。可以看出，其拟合结果的R2 均大于0.99，说明相符程度较好，但是随着k*值的增加，R2 有降低的趋势，这依然是因为当k*较大时，仿真结果的抖动较为剧烈。m 的值关于k*=1.0 对称。当k*=1.0 时，式（5）变为w（x）=uk=1.0（x），m 被消去，因此没有拟合结果，但此时仍有拟合曲线，其与本小节中的拟合结果相同。
表8 非均匀蜂窝拟合结果
Table 8 Non-uniform honeycomb deflection curve fitting results
注：*表示此时m 被消去。
k*方程形式mR2 0.6 w（x） = g1（x） uk = 1.0（x） + g2（x）uk=k*（x）g1（x） =1 1+emx g2（x） =1 1 1 1+ -emx– 0.0160.9949 0.8– 0.0130.9974 1.0—*0.9967 1.2– 0.0130.9968 1.4– 0.0160.9923
对表8 中的m 值进行拟合，得到非均匀蜂窝挠曲线方程系数拟合结果m（k）= – 0.025k2 + 0.05k – 0.037（R2=1），m 关于k 是一个二次函数，且有对称轴k = 1.0。至此，本文所涉及非均匀蜂窝挠曲线已解出。
运用k*= 0.7、0.9、1.1、1.3 时的仿真结果对拟合结果进行验证，结果如图11 所示，可以发现，当k*值较小时，预测曲线与仿真结果较为吻合。当k*= 1.3 时，预测蜂窝的挠度最大值与仿真结果偏差较大，这与之前的拟合分析相呼应。
图11 对非均匀蜂窝挠曲线的验证结果
Fig.11 Comparison of prediction curves and simulation values for non-uniform honeycomb
3 结论
本文以反四手性蜂窝为研究对象，通过数值分析与仿真分析相结合的方法，分析其面外的变形性能，得出以下结论。
（1）对于蜂窝结构的仿真计算，胞元的大小及数量会影响其变形曲线的连续性，当胞元过大或者胞元数量过少时，蜂窝变形连续性较差。此时对蜂窝变形的分析应该从蜂窝结构上入手，利用等效模量的材料分析方法结果与实际差别较大。
（2）反四手性蜂窝在均布压力下的变形曲线与胞元的肋长有关，在一个方向胞元肋长不变的情况下，另一个方向胞元肋长越大，相同载荷下的变形越大。运用曲线拟合的方式可以较为准确地求出其对称轴上的挠曲线。
（3）对于两种不同尺寸的反四手性蜂窝拼接非均匀蜂窝，其对称轴上的挠曲线可以由两种蜂窝各自的挠曲线加权求和得到，加权函数为Logistic 函数。该拟合函数只有一个变量，但是却能取得R2>0.99 的拟合结果。
（4）在与有限元仿真结果误差不大的情况下，拟合曲线可以快速对蜂窝面外变形进行求解。
（5）该数值拟合方法可以用于四手性蜂窝变形预测与设计，使其在自适应变形结构中发挥更大的作用。
[1] 马成, 刘方军. 蜂窝材料加工工艺研究进展[J]. 航空制造技术, 2016, 59(3): 48–54.
MA Cheng, LIU Fangjun. Research progress in processing technology of honeycomb materials[J]. Aeronautical Manufacturing Technology, 2016, 59(3): 48–54.
[2] LIU J F, WANG Z G, HUI D. Blast resistance and parametric study of sandwich structure consisting of honeycomb core filled with circular metallic tubes[J]. Composites Part B: Engineering, 2018, 145: 261–269.
[3] ZHAI J Y, LIU Y F, GENG X Y, et al. Energy absorption of pre-folded honeycomb under in-plane dynamic loading[J]. Thin-Walled Structure, 2019, 149: 106356.
[4] ZHANG D H, FEI Q G, ZHANG P W.In-plane dynamic crushing behavior and energy absorption of honeycombs with a novel type of multi-cells[J]. Thin-Walled Structures, 2017,117: 199–210.
[5] 薛华飞, 姚秀荣, 程海明, 等. 热防护用轻质烧蚀材料现状与发展[J]. 哈尔滨理工大学学报, 2017, 22(1): 123–128.
XUE Huafei, YAO Xiurong, CHENG Haiming, et al. Current situation development of lightweight ablation materials for thermal protection[J]. Journal of Harbin University of Science and Technology, 2017, 22(1): 123–128.
[6] TEJAVIBULYA N, YOUSSEF J, BAO B, et al. Directed self-assembly of large scaffold-free multi-cellular honeycomb structures[J]. Biofabrication, 2011, 3(3):034110.
[7] MEI H, LI C Y, LIU H, et al.Simulation of catalytic combustion of methane in a monolith honeycomb Reactor1[J]. Chinese Journal of Chemical Engineering, 2006, 14(1):56–64.
[8] OLYMPIO K, GANDHI F. Zero-v cellular honeycomb flexible skins for onedimensional wing morphing[C]//Proceedings of the 48th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference. Honolulu: AIAA, 2007.
[9] 郝新超, 胡杰. Nomex 蜂窝夹层结构侧向变形机理及蜂窝稳定化[J]. 航空制造技术, 2020, 63(13): 69–74, 82.
HAO Xinchao, HU Jie. Lateral crush mechanism and honeycomb stabilization of Nomex honeycomb sandwich structure[J].Aeronautical Manufacturing Technology, 2020,63(13): 69–74, 82.
[10] ZHANG C, LIU R H, TAO C C,et al. Design, construction, and modeling of aircraft door sealing plate based on SMAs[J].International Journal of Smart and Nano Materials, 2022, 13(3): 481–503.
[11] 尹维龙, 石庆华. 变体飞行器蒙皮材料与结构研究综述[J]. 航空制造技术, 2017,60(17): 24–29.
YIN Weilong, SHI Qinghua. Review of material and structure for morphing aircraft skin[J]. Aeronautical Manufacturing Technology, 2017, 60(17): 24–29.
[12] ZHANG Q Q, DONG J Q, ZHAO Y T, et al. Three-dimensional meta-architecture with programmable mechanical properties[J].International Journal of Smart and Nano Materials, 2022, 13(1): 152–165.
[13] YIN Y Q, YU Y, LI Bo. Notch flexure as kirigami cut for tunable mechanical stretchability towards metamaterial application[J]. International Journal of Smart and Nano Materials, 2022, 13(2): 203–217.
[14] BUDARAPU P R, SUDHIR SASTRY Y B, NATARAJAN R. Design concepts of an aircraft wing: Composite and morphing airfoil with auxetic structures[J].Frontiers of Structural and Civil Engineering,2016, 10(4): 394–408.
[15] BORNENGO D, SCARPA F,REMILLAT C. Evaluation of hexagonal chiral structure for morphing airfoil concept[J].Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, 2005, 219(3): 185–192.
[16] AIROLDI A, BETTINI P,PANICHELLI P, et al. Chiral topologies for composite morphing structures—Part I:Development of a chiral rib for deformable airfoils[J]. Physica Status Solidi (b), 2015,252(7): 1435–1445.
[17] AIROLDI A, BETTINI P,PANICHELLI P, et al. Chiral topologies for composite morphing structures—Part II: Novel configurations and technological processes[J].Physica Status Solidi (b), 2015, 252(7): 1446–1454.
[18] GONG X 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.
[19] 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.
