1 引言
    众所周知，纤维增强复合材料板壳结构对横向冲击作用十分敏感，较小的横向冲击作用可能导致面内压缩强度较大的损失^[1~3]，因此，这类结构的抗冲击能力一直是人们关注的焦点。求解层合板横向冲击下的瞬态响应和冲击破坏，一般都采用有限元法，并选用合适的本构理论^[4]。
    迄今为止，研究者们采用了不同的复合材料有限元模型，从不同角度去研究层合板的低速冲击性能。Hou^[5]利用LS-DYNA3D有限元程序，选用实体单元，单元弹性矩阵由单向板弹性常数确定，选用修正的Chang-Chang^[6]破坏准则判断纤维和基体的初始破坏以及Brewer-Lagace^[7]建议的破坏准则辨别分层,满足破坏判据的单元采用相应的应力更新方案,以此研究冲击过程中层合板的损伤问题。程小全等^[8]选用三维20节点等参元，单元弹性矩阵由单向板弹性常数确定，破坏判据选用Chang和Springer分层判据^[9]、Hashin基体开裂破坏判据^[10]，对层合板冲击后压缩全过程进行了分析。文献_^[11]利用3DIMPACT有限元程序，选用实体单元，单元弹性矩阵由单向板弹性常数确定、Choi基体开裂破坏判据^[12] 、Choi和Chang分层破坏判据^[13]，分析低速冲击荷载作用层合板的响应。文献[14]利用Abaqus有限元程序，选用壳单元，因壳单元不计平面外应力，在层合板模型中添加纯基体层，以该层破坏等效分层破坏，单元弹性矩阵由单向板弹性常数确定，破坏判据采用能量破坏判据、Hashin破坏判据。如同上述文献，其它大多数有限元模型的单元弹性刚度矩阵由单向板弹性常数确定，虽然选用的破坏判据和单元应力更新方案存在一些区别，但这都需要事先实验确定单向板动态弹性参数和强度参数。随着层合板铺设构造、纤维体积含量、冲击速度等因素的变化，实验将会变得相当烦琐。细观力学有限元模型直接通过两种组分的物理性能参数确定层合板冲击承载能力，仅需实验确定纤维和基体的动态弹性参数和强度参数，能够减少实验的工作量。然而，目前的有限元模型较少采用细观力学有限元模型，主要原因可能是细观力学有限元模型^[15~17]用于解决层合板冲击问题时存在计算量过大的难题。
    本文将复合材料细观力学桥联模型^[18]与有限元软件ABAQUS结合，分析层合板受低速冲击作用的极限冲击响应。即将细观力学本构模型、针对组分材料的破坏判据以及刚度衰减模式编制成用户子程序VUMAT,为ABAQUS求解提供一种自定义材料模型。该模型仅需输入纤维和基体的材料参数、纤维体积含量等有限数据并且无须单层板的实验数据，通过确定组分材料是否破坏来判断单层板是否破坏，对破坏后的单层实施一种常系数刚度衰减，就能使得复合材料层合板结构的冲击承载能力分析得以顺利实施。
2 有限元模拟
   纤维增强复合材料层合板受低速横向冲击作用，该冲击系统包含冲击头、层合板和两者之间的相互作用。层合板的单元选择与低速冲击的研究目的密切相关。如研究目的为冲击后剩余强度,分层则是主要研究对象之一，选择三维实体单元以及三维破坏判据是相当必要的。本文的研究目的为层合板低速冲击作用下的极限承载能力,根据基体和纤维的平均应力判断基体开裂、纤维断裂，对破坏后的单  元组分材料的物理参数进行折减，而不在于识别分层破坏。一方面分层往往受面内基体开裂的诱发^[19],如果已经对发生基体开裂的层进行刚度衰减，就不需要因分层而再次折减该层刚度。另一方面分层往往出现在冲击背面的受拉区^[11]，分层对复合材料面内受拉强度影响并不明显,这种做法的实质是把分层吸收的能量包含在基体开裂中。基于上述理由并减少计算量，本文选择壳单元。
    ABAQUS提供了4节点有限应变壳单元S4R^[20]，每个节点3个平移自由度、3个转角自由度。为了获得单层的应力和应变信息，每个单层采用Simpson 3点积分法，如图1所示，在层合板中面布置节点并划分网格，壳截面广义力沿厚度积分得到。冲击锤假定为刚性体，选择R3D4刚性单元模拟。

              
    图1中，节点1、2、3、4在壳参考平面上，层合板(5层)与S4R复合材料截面定义平面内1点积分；“x”表示厚度上Simpson积分点位置。[-page-] 
    纤维复合材料层合板虽属于薄壳结构，仍需要适当考虑横向剪切，而且对于一阶剪切变形的壳单元，横向剪切作用对动量守恒或虚功方程以及离散结点的动力方程有一定影响。不正确设置剪切刚度将导致收敛困难和计算时间漫长。对称布置的层合板可视为正交各向异性材料，其横向剪切刚度通过下面的简式计算：

         
    其中，下标1,2,3为壳单元坐标方向；D₁₁，D₂₂，D₁₂为复合材料壳截面的抗弯刚度。
    复合材料层合板与冲击头之间的相互作用通过接触对算法实现。接触对中，冲击头为主面，受冲击的层合板为从面，并指明与冲击头表面发生接触的方向。另外从面网格细化须能满足解的稳定性。由于不考虑两者间的切向摩擦，选择主从面动态强化约束接触算法^[20]。
3 本构模型
3.1 应力更新方案
    桥联模型细观力学方法给出了准静态荷载下基体、纤维与单向板的应力-应变之间的增量关系式，但在瞬态动力分析中需要材料率形式的本构方程。本文假定纤维、基体及复合材料为率无关材料。将准静态性能简单推广到动态性能，纤维和基体以及单向复合材料在小变形下率形式的本构方程分别为：

       
    其中，“・”表示对时间求导；V_f和V_m分别为纤维和基体的体积含量；[S]、[S^f]、[S^m]依次为单向板、纤维和基体的柔度矩阵。类似文献[18]的推导，单向复合材料应力-应变率形式的本构方程为：

         
    式中，[A_ij]为桥联矩阵；[I_ij]为单位矩阵。由于在VUMAT中{△εi}已知，材料点在当前积分步的应力增量和总应力为：

              
3.2 破坏判据与刚度衰减
    采用大正应力破坏准则，当纤维或基体满足下列不等式：

             
    则认为与之对应的单层板发生了破坏。其中,上标f,m分别为纤维和基体;下标1、2分别为分量坐标方向;u、c分别为极限拉伸和压缩强度。
    单层板破坏后须进行刚度衰减，后续调用该单元时使用衰减的刚度矩阵。刚度衰减格式及应力更新策略采用经验方法^[5,21]，针对组分材料进行。一旦纤维发生破坏，则认为该层完全丧失承担后续荷载的能力，有：

           
    若基体发生破坏，认为基体开裂，但由于纤维和上下层的约束，该层还能够继续承受荷载作用，泊松比保持不变，只是纤维和基体的当前模量变为：

            
4 算例
    考虑一个周边夹支的[0_n/90_m/0_n]层合板,2n＋m＝10层，密度为1.678kg・m^-3，有效直径为160mm，单层厚为0.18mm,纤维体积含量为V_f＝0.43，中心承受横向低速冲击作用。文献[22]给出的实验数据：各向同性纤维，E_f＝74GPa、G_f＝30GPa、μ_f＝0.25；各向同性基体，E_m＝3GPa、G_m＝1.1GPa、μ_m＝0.4；冲击头质量2.3Kg，头部直径25mm的半圆球型。
    有限元模型见图2，中心区域网格加密见图3，图中数字为单元编号。模型相对几何中心对称，材料性能则相对x、y轴对称，取半圆进行分析，径向设置为对称边界条件。

          
    给定冲击能量1J，即冲击速度为0.933m/s。从图4可以看出，实验值和预报值吻合较好，说明本文采用的细观力学有限元模型计算有效。图5给出了三种不同铺层的层合板在冲击载荷下的中点挠度，彼此差别不大，反映了目前工况下层合板承受冲击作用与层合板的铺设角影响不大，该结论与实验数据一致。[-page-] 
    为了清晰反映基体应力沿径向的变化趋势，选取有代表性的单元3,1,80,144,192,具体位置见图3。上述单元的纤维和基体平均应力变化曲线如图6、7所示。

    图6(a)、(b)表明,在冲击背面0°铺层内基体受双向拉应力，冲击作用点附近的应力峰值远高于其它部分；随着远离冲击中心，基体的应力峰值迅速变小；根据单元3,1,80,144,192的基体应力逐步变小的变化趋势，该层距离中心比单元192单元更远的区域，基体应力峰值较小。图6(e)、(f)表明，冲击正面内的基体受双向压应力，且应力峰值小于冲击背面0°铺层中的应力的峰值；随着远离冲击中心，基体的应力峰值迅速变小；该层距离中心比192单元更远的区域，基体应力峰值小。另外，根据基体拉压性能的差异，容易判断冲击背面将先出现基体拉伸破坏。图6(c)、(d)表明90°铺层内的应力值峰值小且随径向变化平缓、应力梯度小；其它区域基体应力峰值小。
    纤维平均应力随时间变化曲线呈现的规律与基体类似。图7(a)、(b)表明，在冲击背面0°铺层内纤维受双向拉应力，冲击作用点附近的应力峰值远高于其它部分；随着远离冲击中心，纤维的应力峰值迅速变小；根据单元3,1,80,144,192的基体应力逐步变小的变化趋势，该层距离中心比单元192单元更远的区域，纤维应力峰值较小。图7(e)、(f)表明，单元3,1,80冲击正面内的纤维受双向压应力，且应力峰值小于冲击背面0°铺层中的应力的峰值；144、192单元应力随时间变化的过程体现不同特点，先是压应力而后是拉应力。如将冲击作用等效为横力弯曲，一般来说反弯点靠近支座，该点应力应该为压应力，这从侧面反映了冲击作用与横力弯曲的区别。图7(c)、(d)表明90°铺层内纤维的应力值峰值小且随径向变化平缓、应力梯度小；该层距离中心比单元192单元更远的区域，纤维应力峰值很小。
    计算结果(图4、5)表明，三种铺层[0₄/90₂/0₄],[0₃/90₄/0₃],[0₂/90₆/0₂]横向反作用合力曲线差别很小，与实验结论一致。将冲击能量提高到27J，组分材料的强度取值，基体拉伸强度为80MPa、压缩强度120MPa，纤维拉伸强度为2150MPa、压缩强度1450MPa，其余参数同前。达到破坏后，纤维和基体性能按(12)式进行折减，结果见图8，可知计算值与实验值比较吻合。

  
    随着冲击能量不断提高，冲击过程持续的时间越来越短，大冲击合力越来越大，如图9所示。当冲击速度为10m/s时，冲击作用合力-时间曲线发生剧变(图10)，合力的峰值很小。究其原因，是因为大量单元的基体发生破坏，按(9)、(10)式衰减刚度，导致整个结构刚度急剧变小，计算结果并不合理，至少冲击作用合力-时间曲线组没有体现出连续性，这需要对刚度衰减和应力更新方案做进一步研究。尽管如此，从中还是可以判断层合板承受9~10m/s冲击力作用将发生较大程度的破坏。
5 结语
    本文利用复合材料细观力学有限元模型对层合板低速冲击进行了数值模拟，并与实验数据进行对比。分析结果中，除了得到冲击荷载时间曲线、中点位移时间曲线、铺层对抗冲击能力的影响和变形特点等常规结论之外，还直接得到纤维、基体平均应力随时间变化及其在空间上的分布。综上所述，本文所述细观力学方法，能够反应层合板受横向低速冲击响应主要特点，计算过程简洁，具有一定的可行性。
         参考文献
[1] Abrate S. Impact on laminated composites：recent Advance[J]．Ap-pl Mech Rev，1991,44（4）：155-190.
[2] Abrate S. Impact on laminated composites：rencent Advance[J]．Appl Mech Rev，1994,47(11);517-544.
[3] 陈浩然,白瑞祥.低速冲击后含损复合材料夹层板剩余强度研究进展[J].力学进展,2002,32(32):402-413．
[4] 黄争鸣，张华山．纤维增强复合材料强度理论的研究现状与发展趋势-“破坏分析奥运会”评估综述[J]．力学进展,2007,37(1):80-98.
[5] Hou JP，Petrinic N，Ruiz R. Prediction of impact damage in com-posite plates[J]．Composite Science and Technology，2006,60:273-281．
[6] Chang F，Chang K. A progressive damage model for laminated com-posites containing stress concentrations[J]．Journal of Composite Material，1987，21：834-55．
[7] Brewer JC,Lagace PA. Quadratic stress criterion for initiation of delami-nation[J].Journal of Composite Materials，1988,22(12):1141-55.
[8] 程小全,郦正能.复合材料层合板低速冲击后压缩的损伤积累模型[J].应用数学和力学,2005,26(5):569-576.
[9] Chang FK，Springer GS. The strength of fiber reinforced composite bends[J]．Journal of Composite Materials，1986,20：30-45.
[10] Hashin Z. failure ctiteria for unidirectional fiber composites[J].Jourual of Applied Mechanics，1980,47:329-334.
[11] Aslan A，Karakuzu R，Okutan B．The response of laminated composite plates under low-velocity impact loading[J]．Compos-ite Structures, 2003,59：119-127.
[12] Choi HY，Wu HYT，Chang FK. A new approach toward understand-ing damage mechanisms and mechanics of laminated composites due to low-velocity impact：part II -analysis[J]．J Compos Mater，1991,25:1012-38．
[13] Choi HY,Chang FK. A mode for predicting damage in graphite/ep-oxy laminated composites resulting from low velocity point impact[J].Journal of Composite Materials,1992,26(14):2134-69．
[14] Johnson HE，Louca LA，Mouring SE. Current research into model-ling of shock damage to large scale composite panels[J]．Journal of Materials Science，2006，41：1573-1590.
[15] Mayes SJ，Hansen AC. Composite laminate failure analysis using multicontinuum theory[J]．Compos Sci Tech，2004，64;379-394.
[16] Gotsis PK，Chamis CC，Minnetyan L. Prediction of composite lami-nate fracture：micromechanics and progressive frature[J]．Compos Sci Technol，1998，58：1137-1150.
[17] 田峰,胡更开,杨改英.塑性细观力学模型的比较分析[J].北京理工大学学报,1996,16(5):477-482.
[18] 黄争鸣．复合材料细观力学引论[M]．北京：科学出版社，2004；49-60.
[19] 徐颖，温卫东，崔海坡．复合材料层合板低速冲击逐渐积累损伤预测方法[J]．材料科学与工程学报,2006,24(1):77-81．
[20] Abaqus Analysis User Manual，version 6.4.
[21] Hou JP，Petrinic N，RuiZ C. A delamination criterion for laminated composites under low-velocity impact[J]．Composites Science and Technology，2001，61：2069-2074.
[22] Collombet F，Lalbin X,Latailade JL. Impact of laminated compos-ites：physical basis for finite element analysis[J]．Composi Sci and Tech，1998,58:463-478．