高性能纤维增强材料具有轻质高强、可设计性好、耐腐蚀性能优异等特点，由其作为面层与轻质芯材复合而成的夹层结构，可充分发挥出纤维面层的力学性能优势[1―3]。刘伟庆等[4―6]将玻璃纤维复合材料(Glass Fiber Reinforced Plastics，下文简称为GFRP)、木材、竹材以及泡沫组合，采用真空导入工艺制作了GFRP增强的夹层结构，研究其力学性能，并应用于道路、桥梁、建筑等领域。
夹层结构的材料各向异性、分层性导致夹层结构的分析变得较为复杂。大多数理论基于经典弹性力学方法，通过假设不同的位移模式，引入层间应力条件、边界条件。常用理论有单层一阶理论、单层高阶理论、多层理论以及Zigzag理论[7―10]。但直接使用上述理论求解通常较繁锁，且一般仅能针对部分特定板和荷载形式获得精确解。实际分析和设计过程中常借助于：1) 等效截面法：引入了“折线”形式的轴向位移假设，体现了面层硬、芯材软的特点，将夹层结构横截面整体等效抗弯刚度等于各部分抗弯刚度之和，简单灵活[11]；2) 有限元法：其中三维有限元方法无疑能较好描述夹层结构内部应力状况，但建模复杂，计算工作量十分巨大，且大多数商业程序并不能很好支持各向异性材料的三维单元。如将各向异性材料等效成各向同性材料，则大大降低了该方法的准确度。因此基于各类板理论的有限元分析一直为关注的重点，计算成本低，精度相对较高。
经过约30年左右的发展，Zigzag理论已相对比较完善。该方法既避免了单层板理论中人为剪力修正系数的使用，高阶剪切变形理论导致的层间应变不连续问题，也极大降低了多层理论中由于层数的增加带来的计算工作量的增长过多的问题。同时，罚刚度乘子或者横向位移关于位置坐标导数的引入使得满足C0连续的有限单元就能满足要求，C1连续性的问题得到了较好地克服[12―19]。分析表明 Zigzag模型可以比较准确地预测从薄型到厚型复合夹层结构的位移与应力[20―22]。但基于软质芯材非线性位移模式的高阶剪切变形 Zigzag理论的复合材料夹层结构试验与理论分析之间的比较研究相对偏少。本文分别选取近年来文章作者完成的木质芯材—GFRP面板增强层、木芯材—GFRP竹混合增强层、以及泡沫芯材—GFRP面板增强层三种由真空导入工艺制作的夹层板抗弯实验结果作为参照，基于芯材层面外横向二次非线性位移而其它层为常数位移模式假设，以及面内各层沿整个夹层板厚度三次的高阶剪切变形Zigzag理论，运用有限元法对夹层结构挠度位移、特征点极限承载荷载、夹层板内部应力状态等方面进行了分析，并与等效截面法的结果进行了对比，评价了等效截面分析法的适用性。考虑到剪应力对于弱芯夹层结构具有重要的意义，但由于有限元法中直接由材料应力-应变关系与结点位移获得的夹层板内部横向剪应力精度相对较低，本文根据应力平衡方程，运用最小二乘法[23]通过有限元法获得的面板正应力来求解截面剪应力。
1 夹层板试验
图1 四点弯曲试验
Fig.1 Experiment schematic for four point bending test
图2 三点弯曲试验
Fig.2 Experiment schematic for four point bending test
木质芯材—GFRP面板、木质芯材—竹材和GFRP混合面板试件采用单向四点弯曲试验，泡沫芯材—GFRP面板采用单向三点弯试验，试件简支于刚性支敦上，具体加载方案如图 1、图 2(图中LVDT-位移计，SG-应变计)。四点弯曲试验以反力架提供支撑，采用200 kN量程的千斤顶逐级加载，每级荷载为1 kN，分配梁加载点之间间隔250 mm。在三点弯曲试验采用量程为200 kN的SANS试验机加载，加载速率为2 mm/min。夹层板试件的基本材料性能参数由材性试验获得，除在试件跨中和支座处设置位移计外，在试件上下表面层及GFRP层与竹层之间放置应变片。试件几何尺寸、极限状态下承载力与破坏形式见表1和图3，材性参数列于表2。
试验结果表明，横向木质芯材试样破坏形式表现为芯材树脂缝剪切破坏和顶层面板被局部压坏；纵向木质芯材试样破坏形式主要是底层面层纤维材料的强度破坏。横向木质芯材较纵向木质芯材有更高的刚度和承载能力，块状芯材之间，以及芯材与面层之间的树脂层是结构薄弱部位。泡沫芯材试件主要表现为顶层面板屈曲破坏；过于软弱的芯材则由于芯材的局部承压不足导致局部压陷。无GFRP面板夹层板加载达到一定荷载后表现出一定的非线性特征，最终由于下部面板的材料强度不足而破坏。总体而言，GFRP面板的存在有效提高了夹层板抗弯刚度，抑制了板侧向变形，推迟了初始损伤的发生，极大提高了承载能力，表现出较好的线性度，使得与弹性理论的分析结果具有较好的可比性。
表1 试件极限承载力及破坏模式
Table 1 Ultimate bearing capacity and failure modes
注：PS代表芯材为顺纹泡桐木竖直方向放置，通过树脂胶接，沿水平轴排布；P代表芯材为顺纹沿水平轴向放置泡桐木；J代表芯材为顺纹沿水平轴向放置花旗松胶合木；Z代表芯材顺纹沿水平轴向放置竹材；F代表芯材为泡沫；G代表GFRP面层，字母后面的数字代表相应材料层的厚度/mm。
图3 夹层板破坏模式
Fig.3 Failure modes of unidirectional sandwich panel
表2 基本材料性能
Table 2 Material performance
注：表中木材或竹材，下标1表示顺纹纵向，2表示横纹方向；对于GFRP则1表示沿板长度方向，2为横向；树脂缝1-竹与木树脂缝；树脂缝 2-其它树脂缝；泡沫 1-(密度 48 kg/m3)；泡沫 2-(密度199 kg/m3)；泡沫3-(密度413 kg/m3)；除泊松比μ12外，所有材料其它参数单位均为MPa；-表示无此项数据。
2 理论分析模型
2.1 位移模式
假设夹层板平面内位移为不同斜率函数的线性组合，沿夹层板厚度呈三次变化；假设夹层板芯材平面外位移沿厚度为二次抛物线变化，其它层为常数[12]。夹层板面内位移模式如图4(a)，其数学模型如式(1)；面外位移模式如图4(b)，其数学模型如式(2)。
图4 位移模式
Fig.4 Laminated-sandwich displacement mode
其中：u0代表夹层结构的中面；θx是夹层结构中面法线关于z轴的转角；nu、nl分别为中面以上夹层层数和中面以下的夹层层数；βx、ηx为未知高阶常数；
分别为中面上部第i层和中面下部第j层夹层的斜率，
为单位跃阶函数。
其中：wu、wl和w0分别为夹层结构上下面层和芯层中面的平面外位移；l1、l2和l3取相应的拉格朗日函数即，
2.2 本构模型
考虑到单向夹层板受力时，材料弹性主方向与坐标轴平行，所以仅取x-z平面内单层夹层本构关系即：
其中：
k表示夹层结构第k层。
不考虑夹层结构顶面和底面受其它剪应力影响以及各层间的缺陷等，则夹层板顶面层与底面层边界剪应力为零，各层层间剪应力相等。将该应力条件代入式(3)中，并整理为：
其中：
A阵是关于材料本构关系的矩阵。BA向量中最后两项顶层位移及底层位移的引入较好回避了板有限元中C1连续的要求，以C0连续的单元也可获得较好的计算精度。
利用式(4)将式(1)重新整理得：
式中：I=1、2、3或4。当I=1时，a=1；当I=2时，a=z；I=3、4时，a=0。取顶、底及中面层位移为参考广义位移，即：
线弹性范围内位移与应变的关系，即：
将式(2)和式(5)代入式(7)，并整理得：
其中，
H是由bI、lI或bI和lI关于z偏导数组成的矩阵，I取整数1、2、3或4。
3 有限元原理
由势能原理可得：
其中：Us为应变势能；Wl为横向外荷载做的功。基于式(3)、式(7)，应变势能可表达为：
其中，
横向外荷载作用下功可表达为：
其中：w为结构任意一点广义位移；p为对应点处横向外荷载。
结点自由度如式(6)所示，取
将式(9)改写成单元形式，即：
其中：
N为对应于各结点横向位移的非零形函数。ε用参考广义位移ε表达为：
其中，iδ为对应于各结点的位移向量，对应于单元节点i的Bi如下：
式中，Ni,x表示Ni对x的导数。
取
则可得：
将各单元刚度和荷载组装后，运用式(14)可获得对应的结点位移，并根据式(3)、式(7)、式(8)进一步计算各点应力。
4 最小二乘法剪应力分析
有限元方法中直接基于材料本构关系获得横向剪应力精度较低，本文采用最小二乘法[18]，利用单元内各点应力平衡的条件分析横向剪应力，可有效提高分析的精确度，这对于夹层结构芯材的准确分析有着重要意义。
不考虑体力和面力作用下，单向夹层板平面应力微分形式的方程如下式：
其中，σ和τ分别表示正应力和剪应力微分形式。
夹层板横向剪应力沿整个板高呈非线性分布，假设任意一点剪应力按照二次分布，在任意一夹层中可表示为多项式形式，如下式：
其中：p表示第p层；
是未知系数。认为层层间剪应力相等，并设夹层结构顶层与底层边界剪应力为零可得下式：
其中：
分别表示夹层板底部剪应力、沿高度方向剪力以及顶部剪应力；zpb和zpt分别为第p层的底部和顶部的坐标；
将式(17)写成矩阵形式，求逆后代入式(16)整理后得：
其中：
A为关于zpb、
的矩阵；
是关于夹层板高度坐标的形函数，即
其中，
是A中第i行、第j列元素。
设式(15)获得的夹层板第p层应力满足下式：
其中，
是根据有限元分析结果，直接应用应力-应变本构关系计算而得。按照最小二乘原则，式(20)成立，即：
利用式(18)、式(19)将式(20)展开得：
其中：
表示iφ对z的导数。将式(8)代入式(3)中，并对x求导，则获得
根据夹层结构层间剪应力连续的条件，利用式(21)单层剪力
组装成所需求解的包含所有层的拟求取单元的剪应力方程，即：
代入夹层板顶部和底部剪应力为零的条件，即可解得所有层剪力和剪应力。
5 等效截面法
等效截面法考虑夹层结构处于弹性，设横截面沿轴向变形满足平截面假定，将芯材以外其它层均以芯材弹性模量换算，则截面总惯性矩如下式：
其中：Ii为各材料层对自身层重心的惯性矩；iA为各材料层横截面面积；yi为各材料层截面重心点坐标。
其各点正应力如下式：
其中：M为截面弯矩；y为拟计算点相对中和轴坐标；Eα为芯材弹性模量相对于其它层弹性模量的比值。
将式(24)代入式(15)，等号两边沿夹层板厚度方向积分，整理后可得下式：
其中：V为拟计算点截面处剪力；Sc为拟计算点到相应边缘处的面积矩；tw为夹层板宽度。
6 分析结果对比
6.1 有限元分析模型
为提高分析精度，有限单元法采用基于拉格朗日插值的三结点二次等参杆单元。试算分析表明，挠度与正应力随着单元划分的增加迅速收敛，取10个单元划分时相对误差已达 1%以下；剪应力收敛相对较慢，取30个单元划分时相对误差达5%以下。本文分析过程中沿夹层板长度均匀划分40个单元，且有限元模型仅包含试验中两支座之间的板。边界采用简支条件，即一支座处节点的自由度u0、w0、wu、wl位移约束为0，θx、uu、ul自由；另一支座节点w0、wu、wl位移约束为0，u0、θx、uu、ul自由。
6.2 分析结果对比
6.2.1 跨中挠度
对应于木质芯材在20 kN、泡沫芯材在3 kN，以及各试件的试验极限荷载下，表3分别给出了试验实测、有限元和等效截面方法的夹层板底部跨中竖向挠度。图 5以 P48-Z6-G0、PS48-Z6-G0、P48-Z3-G3、P54-Z0-G3以及 F70-G3.6-2夹层板为例，给出了试验值、有限元和等效截面法计算获得的荷载-跨中挠度对比曲线。其中 F70-G3.6-2给出了试验和有限元法计算下上面层与底面层位移，以及等效截面法计算的等效位移。
表3 跨中底部点位移
Table 3 Bottom displacement of middle span
注：表3中等效截面法计算过程中泡沫芯材弹性模量等参数近似取拉伸和压缩平均值；“20(3)/kN”表示对于木质心材取20 kN计算，而泡沫芯材则取3 kN计算；泡沫夹层板挠度相对误差取试验的板底挠度作为比较对象。
图5 夹层板跨中挠度
Fig.5 Middle-span deflection of sandwich panel
表3和图5表明，相对于树脂胶结竖向放置的顺纹泡桐木芯材的夹层板，沿板轴向顺纹连续整体泡桐木芯材的夹层板线性度较好，基于软质芯材的高阶 Zigzag理论的有限元方法对于夹层板的跨中挠度的预测具有较好的准确度。设有GFRP面层的夹层板的线性度较好；未设有GFRP面层的树脂胶结竖向放置的顺纹泡桐木芯材的夹层板试件加载超过一定荷载后明显呈现非线性特征。基于线弹性理论分析方法的计算结果准确度相对较低。对于较硬的木质芯材夹层板，基于软质芯材的高阶Zigzag理论的有限元方法和等效截面法获得的位移的计算精度相当；但对于除GFRP面层外具有多层特征的夹层结构以及软质泡沫芯材的位移计算而言，基于高阶 Zigzag理论的有限元方法明显较等效截面的方法有较高准确度。
6.2.2 特征荷载
鉴于夹层板破坏形式(见表1)并结合加载的对称性特点，选取1、2、3、4、5五点作为分析参考点，具体如下：木芯、竹材与GFRP的夹层板，1、2、3、4、5位置见图1，分别对应于竹—木交界面树脂缝、竹材外侧受拉边缘、竹—GFRP层交界面树脂缝、载荷加载处底面GFRP层受拉外侧边缘、跨中底面GFRP层受拉外侧边缘；泡沫夹层板，1、2、3、4、5位置见图2，分别对应于跨中顶面GFRP外侧边缘、顶面GFRP层—泡沫芯材交界面泡沫或树脂缝、泡沫芯材中心、底面GFRP层—泡沫芯材交界面树脂缝、底面GFRP层外侧边缘。采用Mises屈服准则，以材料强度破坏作为特征控制分析了上述五点分别达到极限强度时的对应的荷载，如表4所示。
考虑到本文分析的试件由真空导入工艺制得，树脂实际作为单独的层存在于不同的芯材之间或芯材与面板之间。树脂层非常薄，且树脂硬化后的力学性能特征较难全面获得，因此直接采用 Mises屈服准则对其进行分析判断存在一定困难。但树脂层的应力状态实质上由树脂层两侧不同材性的材料的应力状态共同决定的，所以本文假设该层应力状态为两侧材料应力状态差值确定，即认为该应力状态差值达到屈服或极限状态后会引起树脂层某些点的破坏，如下式：
表4 各点极限荷载
Table 4 Ultimate load
注：表4中( )内表示有限元计算的极限荷载与试验值的相对误差百分比；Pu为试验极限荷载；—表示无此项数据。
其中：σs表示树脂缝处应力状态；σs_1和σs_2分别表示树脂缝两侧同一点位应力状态。
结合表1和图3试件破坏情况，以线弹性高阶剪切变形 Zigzag理论的有限元方法计算的木质芯材夹层板极限荷载与试验结果的相对误差大体在10%~20%，泡沫芯材极限荷载在33%以内。总体上有限元计算的荷载小于试验的极限荷载，这主要是由于荷载判断准则仅以个别边缘点是否达到极限状态作为准则，偏于保守。需要指出的是 F70-G3.6-1在3 kN集中荷载作用下，荷载作用点下方的泡沫受到的局部压应力已经达到3.82 MPa，超出了泡沫材料的抗压强度，加载点处的泡沫必然发生压陷。
表4表明夹层板在加载初期较低的荷载作用下，树脂缝处部分点位已经达到其极限强度，树脂缝处部分点位的破坏仅仅使得试件刚度少量下降，并不主要影响不同材料层间应力的传递，GFRP纤维面层对于芯材和树脂缝部分点位达到极限状态后的破坏具有有效的约束作用，提高了夹层板的刚度和承载能力。可以认为构件在该阶段仍然按照弹性状态工作，芯材和树脂缝层部分点位破坏不同程度的表现出弱非线性特征。所以对于具有一定厚度的GFRP增强层的夹层板按照线弹性来分析具有较好的准确度，基于高阶剪切变形Zigzag理论的有限元方法可以用于夹层板极限荷载的预估。
6.2.3 截面应力
表5给出了木质芯材夹层板在20 kN荷载作用下以及泡沫芯材在3 kN荷载作用下，分别由试验、有限元和等效截面方法计算的构件受拉区最大正应力和芯材中和轴处最大剪应力值的比较。其中“加载点附近”指由有限元法计算的中和轴处具有最大剪应力的截面。图6以P48-Z3-G3试件为例给出了截面沿轴向正应力和剪应力分布图，其中图6(a)为有限元跨中、分配梁加载点处及等效截面法求得的跨中加载点处的受拉应力分布；图6(b)为有限元分配梁加载点及加载点附近截面，以及等效截面求得的中和轴最大剪应力。由于数值方法本身精度误差问题，因此有限元截面剪应力计算过程中采用前文所述的最小二乘法。
表5 应力比较
Table 5 Stress compare
注：G为GFRP外边缘；B为竹材受拉外边缘；W为木芯外边缘。
图6 P48-Z3-G3截面正应力和剪应力
Fig.6 Normal stress and shear stress of P48-Z3-G3
表5和图6表明，对于木质芯材夹层板，除少数构件外基于高阶 Zigzag理论的有限元方法与等效截面法在截面正应力计算中是较为接近，且与试验结果有较好的一致度。但对于泡沫芯材夹层板截面正应力以及木质和泡沫芯材夹层板剪应力计算有着明显差异。等效截面法获得的泡沫芯材夹层板正应力明显偏小，而有限元法符合度较好。图7给出了泡沫芯材夹层板泡沫芯材沿板长度方向顶、底面层的应变试验值，在加载点附近明显存在突变。相对而言，本文有限元方法一定程度上较好的体现出了这部分特征。两种方法在中和轴处获得的最大剪应力数值上大体一致，但由等效截面法获得的剪应力在截面上的分布存在畸变。有限元法计算表明由于夹层板各层材料性质上的差异，最大剪应力不一定在中和轴处，这对于如泡沫等弱芯材有着重要意义。在基于高阶Zigzag理论的有限元方法计算基础上采用最小二乘法求得剪应力总体上更加合理。
图7 应变分布图
Fig.7 Strain distribution
7 结论
(1) 适当 GFRP增强层的存在可以抑制夹层板内部芯材破坏的发展，使得夹层板在正常工作状态下，总体呈现出较好的线弹性。采用基于线弹性理论的分析方法是合理的。
(2) 考虑芯材软质芯材非线性位移模式的高阶剪切变形的Zigzag理论有限元方法与试验结果有良好的符合度，适用于GFRP增强夹层板的精细化分析、设计。相对于等效截面法，本文有限元方法在软质芯材、多层正交异性夹层板，以及板内剪应力分析等方面更具有优势，为复杂夹层结构在复杂荷载工况下的分析和设计提供了非常有效的、实用的和通用的方法。
(3) 等效截面法对于硬质芯材的夹层板的分析具有一定准确性，适用于GFRP增强夹层板的初步设计和分析。其芯材等效过程简化了分析方法，提高了分析效率，但同时也降低了部分分析精度。
参考文献：
[1]刘伟庆, 方海.纤维增强复合材料及其在结构工程中的应用研究[C].2010年海峡两岸功能性复合材料论坛论文集.北京: 中国建筑工业出版社, 2010: 25―35.Li Weiqing, Fang Hai.Research on the application and development of fiber reinforced composites in structures[C]// Proceedings of the 2010 Cross Strait Forum on Functional Composite Materials.Beijing: China Architecture & Building Press, 2010: 25―35.(in Chinese)
[2]陈跃良, 雷园玲, 张勇, 等.T300/BPM316复合材料层合板的载荷放大系数确定[J].工程力学, 2016 (1):195―200.Chen Yueliang, Lei Yuanling, Zhang Yong, et al.Determination of load enhancement factor of T300/BPM316 composite laminates[J].Engineering Mechanics, 2016 (1): 195―200.(in Chinese)
[3]谭萍, 聂国隽.曲线纤维增强复合材料圆环板的非轴对称弯曲问题[J].工程力学, 2016 , 33 (3): 239―247.Tan Ping, Nie Guojun.Asymmetric bending of curvilinear-fiber reinforced composite annular plates [J].Engineering Mechanics, 2016, 33(3): 239―247.(in Chinese).
[4]Fang Hai, Sun Huiming, Liu Weiqing, et al.Mechanical performance of innovative GFRP-bamboo-wood sandwich beams: Experimental and modelling investigation [J].Composites Part B: Engineering, 2015,79: 182―196.
[5]Fang Hai, Mao Yifeng, Liu Weiqing, et al.Manufacturing and evaluation of Large-scale Composite Bumper System for bridge pier protection against ship collision [J].Composite Structures, 2016, 158:187―198.
[6]洪俊青, 刘伟庆, 方海.格构腹板增强型泡沫复合材料夹层结构腹板平面内屈曲分析[J].工程力学, 2015,32(08): 66―73.Hong Junqing, Liu Weiqing, Fang Hai.Web buckling behaviors under in-plane compression for web lattice-enhanced foam-sandwich structure [J].Engineering Mechenics, 2015, 32(08): 66―73.(in Chinese)
[7]Goyal V K, Kapania R K.A shear-deformable beam element for the analysis of laminated composites [J].Finite Elements in Analysis & Design, 2007, 43(6):463―477.
[8]Reddy J N.A simple higher-order theory for laminated composite plates[J].Journal of Applied Mechanics, 1984,51(4): 745―752.
[9]Reddy J N, Robbins D H.Theories and computational models for composite laminates[J].Applied Mechanics Reviews, 1994, 47(6): 155―199.
[10]Murakami H.Laminated composite plate theory with improved in-plane responses [J].Journal of Applied Mechanics, 1986, 53(3): 661―666.
[11]Allen H G.Analysis and design of structural sandwich panels [M].Oxford: Pergamon Press, 1969.
[12]任晓辉.复合材料层合板 C0型高阶锯齿理论[D].大连: 大连理工大学, 2014.Ren Xiaohui.C0-type higher-order zig-zag theory for laminated composite plates [D].Dalian: Dalian University of Technology, 2014.(in Chinese)
[13]任辉, 陈万吉.基于 C0型锯齿理论的软核夹层梁静力和稳定分析[J].工程力学, 2013, 30(4): 47―51.Ren Xiaohui，Chen Wanji.Static and stability analysis of soft-core sandwich beams in terms of C0-type zig-zag theory [J].Engineering Mechanics, 2013, 30(4): 47―51.(in Chinese)
[14]Sarangan S, Singh B N.Improved zigzag theories for laminated composite and sandwich plates with interlaminar shear stress continuity [J].Aerospace Science & Technology, 2016, 52: 243―255.
[15]Han J W, Kim J S, Cho M.Generalization of the C0-type zigzag theory for accurate thermomechanical analysis of laminated composites [J].Composites Part B, 2017, 122:173―191.
[16]Pandit M K, Sheikh A H, Singh B N.An improved higher order zigzag theory for the static analysis of laminated sandwich plate with soft-core [J].Finite Elements in Analysis & Design, 2008, 44(9/10): 602―610.
[17]Pandit M K, Sheikh A H, Singh B N.Buckling of laminated sandwich plates with soft core based on an improved higher order zigzag theory [J].Thin-Walled Structures, 2008, 46(11): 1183―1191.
[18]Chakrabarti A, Chalak H D, Iqbal M A, et al.A new FE model based on higher order zigzag theory for the analysis of laminated sandwich beam with soft core [J].Composite Structures, 2011, 93 (2): 271―279.
[19]Kapuria S, Kulkarni S D.An improved discrete Kirchhoff element based on third order zigzag theory for static analysis of composite and sandwich plates[J].International Journal for Numerical Methods in Engineering, 2007, 69(9): 1948―1981.
[20]Icardi U.Higher-order zig-zag model for analysis of thick composite beams with inclusion of transverse normal stress and sublaminates approximations [J].Composites Part B Engineering, 2001, 32(4): 343―354.
[21]Zhen W, Wanji C.An assessment of several displacement-based theories for the vibration and stability analysis of laminated composites and sandwich beams [J].Composite Structures, 2008, 84(4): 337―349.
[22]Icardi U, Sola F.Development of an efficient zig-zag model with variable representation of displacements across the thickness [J].Journal of Engineering Mechanics, 2014, 140(3): 531―541.
[23]Bhar A, Satsangi S K.Accurate transverse stress evaluation in composite/sandwich thick laminates using C0HSDT and a novel post-processing technique [J].European Journal of Mechanics A/Solids, 2011, 30 (1):46―53.