复合材料加筋壁板的抗冰雹冲击动力响应及损伤预测*

张晓晴 丁铁 龙舒畅 姚小虎

(华南理工大学 土木与交通学院, 广东 广州 510640)

复合材料加筋壁板的抗冰雹冲击动力响应及损伤预测*

张晓晴 丁铁 龙舒畅 姚小虎

(华南理工大学 土木与交通学院, 广东 广州 510640)

为了解复合材料加筋结构受到冰雹冲击后的损伤特性,基于软件Abaqus建立了复合材料加筋结构的有限元模型,采用光滑粒子法模拟了冰雹撞击的过程,分析了冲击位置、冲击能量及冲击入射角对复合材料加筋壁板冲击动力响应的影响,探讨了面板与筋条的分层和面板的损伤情况，通过对比数值模拟与试验结果验证了文中所建有限元模型的准确性.在各种损伤形式中,基体的拉伸损伤和面板与筋条的分层损伤最为显著;冲击点为筋条的边缘时,面板与筋条的分层损伤最为严重.

复合材料加筋壁板;冰雹;冲击;动力响应;损伤；有限元模型

复合材料因其优良的材料性能而被越来越多地应用到飞机设计中,在波音公司最新研制的梦幻客机787中,复合材料所占比例已超过50%[1].飞机在服役期间往往会遇到冰雹等恶劣天气,其表面结构面临着冰雹冲击的威胁.在世界范围内,每年因冰雹和鸟体撞击造成的损失高达30亿美元[2].复合材料加筋壁板在飞机结构中逐渐作为主承力构件使用,在整个结构体系中更为重要.如果主承力构件失去承载能力,将会对飞机的使用造成严重的安全隐患,因此十分有必要对复合材料加筋壁板的抗冰雹冲击能力进行研究,为复合材料加筋壁板的设计提供一定的参考依据.

Found等[3- 4]通过试验研究了碳纤维增强复合材料加筋壁板受冲击荷载后的动力响应,但选取的冲击落点仅仅为两筋条中间的区域,并没有考虑筋条边缘和筋条上受冲击后的损伤.Faggiani等[5]建立了复合材料加筋壁板的低速冲击模型,其计算结果与试验结果吻合较好,但计算效率太低.邵青等[6]进行了复合材料加筋壁板的冲击试验,但没有涉及数值模拟方面的研究.当前的研究主要集中在刚性体低速冲击复合材料结构,而将可变形体如冰雹作为冲击源的研究很少.冰雹因其独特的晶体结构,在高速冲击下会发生破碎,同时由于应变率效应而呈现出很强的非线性特征.冰雹的这种材料特性大大提高了建立准确、有效冰雹模型的难度.虽然目前已有冰雹冲击复合材料结构方面的研究,但这些研究主要局限于层合板结构.如刘洋等[7]通过模拟不同速度下的冰雹冲击复合材料层合板,研究了复合材料层合板的临界破坏速度和破坏形式;胡宗文等[8]研究了在冰雹高速冲击下复合材料板壳的结构参数对冲击结果的影响规律.

为此,文中采用光滑粒子法建立冰雹模型,在复合材料结构中引入损伤演化和层间损伤,模拟了冰雹冲击复合材料加筋壁板的过程;选取了3个典型的冲击位置,研究不同冲击速度和冲击角度下结构的动力响应.

1 复合材料损伤模型

1.1 层内损伤模式

复合材料的损伤可分为层内损伤和层间损伤两类.层内损伤又可分为纤维拉伸断裂、纤维挤压破坏、基体拉伸开裂和基体压缩破坏4种形式.Hashin失效准则能比较准确地判断层内损伤[9],在业界得到了广泛的应用.该失效准则认为,层内损伤的4种形式相互独立,当某一单元的应力满足其中一项表达式时,该单元就发生相应形式的损伤.其具体的判断条件如下:

(1)

(2)

(3)

(4)

图1 层内损伤模型

1.2 层间损伤模式

层间损伤模型是用来模拟层与层之间界面的损伤与损伤演化.分层的起始和扩展通常是在混合模式载荷作用下发生的,采用混合模式下的界面损伤准则比单一模式下的损伤准则更为精确.这里采用二次名义应力准则来判断损伤的起始[10],当其值大于1时,认为界面发生损伤:

(5)

(6)

采用BK能量准则[11]来控制层间损伤的演化:

(7)

2 模型验证

为了验证模型的有效性,将数值模拟的结果与文献[12- 13]的一组试验数据进行对比.该组试验是将一个直径为42.7mm的冰雹以73.2m/s的初速度冲击一块尺寸为304.7mm×304.7mm的复合材料层板,层板四周简支,厚度为2.44mm,铺层顺序为[0/45/90/-45]2S.

冰雹在高速冲击下会表现出一定的流体特性,传统的有限元法难以真实地模拟冰雹的冲击过程.目前模拟冰雹的方法有拉格朗日法、任意拉格朗日-欧拉法和光滑粒子法,Anghileri等[14]比较了这3种方法,发现光滑粒子法的效果最好.光滑粒子法是一种无网格化的拉格朗日计算方法,该方法将连续结构体离散成等质量的质点,通过求解质点组的动力学方程及每个质点的运动轨迹来求出整体结构的力学行为[15].文中将冰雹模型离散成67 048个等质量的质点,冰雹的材料属性为:ρ=846kg/m3,E=9 203.6MPa,σyield=10.3MPa,ν=0.33,σcutoff=4MPa.定义冰雹的破坏需要在Abaqus的input文件中添加“Tensilefailure”语句,语句的内容为“Elementdeletion=no,Press=ductile,Shear=brittle”.如果不定义这个语句,冰雹在模拟过程中就不会发生破碎.冰雹在高速撞击条件下会发生硬化现象,高应变率的硬化作用改变了冰雹的力学性质,屈服强度发生动态上升.由于屈服强度并不唯一,故文中对塑性阶段屈服强度值的设定,采用了基于应变率的强度设置,屈服强度比与应变率的关系如表1所示[16].

表1 冰雹的屈服强度比与应变率的关系[16]

Table1Relationshipbetweenyieldstrengthratioandstrainrateofhail

屈服强度比应变率/s-11．0000．01．0100．11．4960．51．7091．02．2055．02．41810．02．91450．03．127100．0屈服强度比应变率/s-13．623500．03．8361000．04．3325000．04．54510000．05．04150000．05．254100000．05．751500000．05．963100000．0

图2为复合材料层合板冰雹冲击试验数据与数值模拟结果的比较,二者吻合较好,可见文中所建立的冰雹冲击复合材料结构模型是合理的.图3为数值模拟给出的各种能量分配,动能、内能和损伤耗散能满足能量守恒定律,伪应变能较小,为总能量的7%,从能量角度来说,计算结果是合理的.

图2 位移时程曲线数值模拟与试验结果比较

Fig.2 Comparison of displacement versus time history curves between numerical and experimental results

图3 几种能量的时程曲线

3 复合材料加筋壁板模型

复合材料加筋壁板模型采用SPH-FEM耦合算法,其基本思想是对问题域中的大变形区域采用光滑粒子法计算,而对小变形区域采用有限元法计算.这样就使算法既具有良好的大变形模拟能力,又具有较高的计算效率.文中采用固定耦合算法,即在计算的初始时刻就确定采用有限元法和光滑粒子法计算的区域,在后续过程中固定不变.

图4 工字型复合材料加筋壁板示意图(单位：mm)

Fig.4 Schematic diagram of I-shaped composite stiffened panel(unit：mm)

表2 界面的材料属性

为了提高计算效率,复合材料加筋壁板模型采用局部细化方式,将面板分为两个区域:冲击区域和其他区域.在一个80 mm×80 mm的正方形区域内使用 6个部件来模拟面板的冲击区域,区域外的面板使用1个部件来模拟,两个区域间使用绑定连接,采取两端简支的边界条件,在部件与部件间使用零厚度的单元来模拟层间分层,有限元模型及其细节如图5所示.文中采用光滑粒子法建立冰雹模型,冰雹被离散为67 048个粒子.复合材料加筋结构采用连续壳单元进行离散,在整个模型中共使用了27 038个连续壳单元.面板冲击区域的网格局部细化,边长为 2.5 mm,共有6 144个单元,其他区域网格边长为5 mm,面板其他区域与筋条的单元数目分别为6 854和14 040个.经验证,此时网格已经收敛.

图5 冲击位置A的有限元模型

冰雹的冲击速度(v)设置在75～165 m/s之间,最低速度为75 m/s,速度间隔为15 m/s,共7个速度工况.此速度范围基本上可以覆盖飞机在起飞、降落以及正常飞行时所遇到的冰雹速度.

4 结果和讨论

4.1 冲击速度的影响

图6给出了冲击点A、B在不同冲击速度下的冲击力时程曲线.冲击点A、B的冲击力峰值随着冲击速度的增大而增大,而冲击的持续时间变化不大.图 6(a)中的冲击力达到峰值后,下降迅速;图6(b)中冲击速度为165 m/s的曲线在冲击力达到峰值之后有一明显的下降过程,之后为一平台段.图6中冲击力时程曲线振荡十分剧烈,这是由于采用光滑粒子法建立冰雹模型,通过网格划分将模型离散为多个相互作用具有质量的粒子,这种方法可以很好地模拟冰雹在高速撞击过程中的流体特性,但会造成冲击力时程曲线振荡.当模型离散为等质量的粒子后,在冲击的过程中,由于粒子与结构接触有先后,使得结构受到的力不连续,从而出现振荡,可以通过网格的疏密来控制振荡的程度.使用不同冲击速度的冰雹撞击点A、B、C,发现冲击力峰值与冲击能量呈线性关系,如图7所示,这与文献[12]得出的规律一致.

图6 不同冲击速度下的冲击力时程曲线

Fig.6 Impact force versus time history curves at different impact velocities

图8为冲击点A、B在不同冲击速度下的位移时程曲线.随着冲击速度的增大,冲击点位置的最大位移也在变大,而冲击持续时间变化不大.图8(a)中冲击速度为135 m/s的位移时程曲线存在两个波峰,其原因有:①由于高速冲击下冰雹表现为脆性,与面板接触后不断破碎,接触面积不断变大,腹板两侧面板会受到冰雹碎粒的作用;②冲击点A区域的刚度不连续,当点A处的位移达到第一个峰值后,腹板开始回弹,而腹板两侧面板上的点由于冰雹碎粒的作用继续向下移动,当腹板回弹的力小于两侧冰雹碎粒的作用力时,腹板随着面板再次向下移动,点A的位移达到第二个峰值.当冲击速度较小时,腹板回弹的力大于两侧冰雹碎粒的作用力,此时位移时程曲线不存在这个现象.当冲击速度较大时,面板会延滞筋条的回弹,第一个波峰掩盖了第二个波峰.图8(b)的位移时程曲线在低速时需要较长的回弹时间,这是由于冰雹在回弹时没有完全碎裂,未破坏部分阻碍了结构的回弹.

图7 冲击力峰值与冲击能量的关系

图8 不同冲击速度下的位移时程曲线

Fig.8 Displacement versus time history curves at different impact velocities

冲击位置为筋条边缘时,不同冲击速度下面板的基体拉伸损伤情况如图9所示.由冰雹冲击所产生的损伤面积随着速度的增加而变大,冲击点两侧的筋条限制了基体拉伸损伤横向扩展,使得损伤向纵向扩展.图10为冲击点B在不同冲击速度下面板与筋条间的分层损伤情况.当冲击速度小于105 m/s时,面板与筋条间没有发生分层损伤.点B处的分层形状呈对称的半椭圆形,椭圆的长轴为工字型筋条上翼缘边缘,随着冲击速度的增加,损伤沿长轴扩展.面板受到冰雹的挤压后,冲击点两侧面板向上鼓起,面板与筋条的变形不一致,从而导致分层,如图11所示.

4.2 冲击位置的影响

当冲击速度为135 m/s时,冲击点A、B、C的冲击力时程曲线如图12所示,3个位置的冲击力时程曲线在最初阶段的斜率吻合.虽然3个冲击位置的刚度不同,但此时的冲击速度高,使得冲击力在很短的时间内达到峰值,这也是冲击力时程曲线在最初阶段斜率非常陡峭的原因.冲击点为A时的冲击力在达到峰值之后,因为A点的正下方为筋条的腹板,冲击点的刚度提高,位移减小,冰雹在冲击过程中更早的碎裂,从而使冲击力迅速下降.

图9 不同冲击速度下基体的拉伸损伤

图10 不同冲击速度下面板与筋条间的分层损伤

图11 点B处面板与筋条翼缘的变形图

Fig.11 Deformation of the panel and stiffener flange at pointB

图12 不同冲击位置的冲击力时程曲线

Fig.12 Impact force versus time history curves at different impact locations

冲击速度为135 m/s的冰雹在冲击点A、B、C的位移时程曲线如图13所示.对比发现,3个位置的位移差别很大.冲击点的刚度越小,位移越大,冲击的持续时间越长.

图13 不同冲击位置的位移时程曲线

Fig.13 Displacement versus time history curves at different impact locations

冲击速度为165 m/s时,冲击点A、B、C的基体拉伸损伤和面板与筋条间的分层损伤情况如图14、15所示.相比之下,冲击点A处的各种形式的损伤面积最小.在冲击位置B、C处,面板产生的损伤基本上被限制在冲击点两侧的筋条区域.对于面板与筋条间的分层损伤,冲击点A的分层面积为585 mm2,冲击点B的分层面积为4 863 mm2,冲击点A、B的损伤面积差别很大.当冲击点为点C时,面板与筋条间没有产生损伤,这说明冲击位置不同,结构的吸能模式也不同.

图14 不同冲击位置处的基体拉伸损伤对比

图15 不同冲击位置处面板与筋条的分层损伤对比

Fig.15 Comparison of delamination damage between panel and stiffener at different impact locations

4.3 冲击角度的影响

图16为不同入射角度(冰雹速度方向与复合材料表面之间的夹角)冰雹冲击下点A的冲击力时程曲线.模拟中选取的入射角为30°、45°、60°、75°、90°,冰雹的冲击速度为135 m/s.随着入射角的增加,冲击力的峰值在不断增大,当入射角度大于45°后,除了冲击力时程曲线的峰值段有区别外,上升段和下降段的走势趋于一致.入射角为30°时冰雹的变形如图17所示,在冰雹与面板接触0.4 ms时,冰雹与面板未接触部分的形状仍基本保持完好;在冰雹与面板接触0.6 ms时冰雹并没有完全破碎,此时冰雹已经不在最初的冲击位置,冰雹的动能没有完全耗散完.失效后的冰雹形成碎末状粒子,呈现出流体特性.

图16 不同入射角度下的冲击力时程曲线

Fig.16 Impact force versus time history curves at different incident angles

图18为冲击点A处在不同冲击角度下的位移时程曲线,此时冰雹的冲击速度为135 m/s.随着入射角的增加,冲击点的位移不断增大,但冲击持续时间基本不变.冲击角度为75°和90°时的位移时程曲线趋势一致,同时曲线出现两个波峰,这说明冲击速度满足一定条件后位移时程曲线才会出现两个波峰,同时较大的斜撞角度对冲击点位移的影响与正撞差别不大.

图17 入射角为30°时冰雹的变形图

图18 不同冲击角度下的位移时程曲线

Fig.18 Displacement versus time history curves at different impact angles

5 结论

文中采用光滑粒子法建立了冰雹模型,在复合材料加筋壁板模型中考虑了损伤演化,研究了在冰雹冲击下冲击速度、冲击位置、入射角度对复合材料加筋壁板的影响,得出以下主要结论:

(1)冲击力峰值与冲击能量呈线性关系.

(2)当冰雹冲击筋条正上方时,冲击速度大于某一临界速度后,位移时程曲线出现两个波峰,随着冲击速度的增加,两个波峰间的界限逐渐模糊.

(3)复合材料加筋壁板在冰雹冲击下的主要损伤模式为基体的拉伸损伤和面板与筋条的分层.当冰雹冲击筋条的边缘时,面板与筋条的分层损伤最为严重,分层呈对称的半椭球形,随着冲击速度的增加,分层沿纵向扩展.

(4)对于冰雹冲击,正撞比斜撞对结构损伤的影响更为严重.

[1] 范玉青,张丽华.超大型复合材料机体部件应用技术的新进展:飞机制造技术的新跨越 [J].航空学报,2009,30(3):534- 543. FAN Yu-qing,ZHANG Li-hua.New development of extra large composite aircraft components application technology:advance of aircraft manufacture technology[J].Acta Aeronautica et Astronuatica Sinica,2009,30(3):534- 543.

[3] FOUND M S,HOWARD I C,PARAN A P.Impact perforation of thin stiffened CFRP panels [J].Composite Structures,2000,48(1/2/3):95- 98.

[4] FOUND M S,HOWARD I C,PARAN A P.Impact beha-viour of stiffened CFRP sections [J].Composite Structures,1997,39(3/4):229- 235.

[5] FAGGIANI A,FALZON B G.Predicting low-velocity impact damage on a stiffened composite panel [J].Compo-sites Part A:Applied Science and Manufacturing,2010,41(6):737- 749.

[6] 邵青,何宇廷,张腾,等.复合材料加筋板低速冲击损伤及剩余压缩强度试验研究 [J].复合材料学报,2014,31(1):200- 206. SHAO Qing,HE Yu-ting,ZHANG Teng,et al.Experimental research on low-velocity impact and residual compre-ssive strength of composite stiffened panels [J].Acta Materiae Compositae Sinica,2014,31(1):200- 206.

[7] 刘洋,王富生,岳珠峰.冰雹冲击复合材料层合板失效模式分析与数值模拟 [J].振动与冲击,2011,30(9):150- 154. LIU Yang,WANG Fu-sheng,YUE Zhu-feng.The failure model analysis and simulation of hailstone impacting with laminated composite plate [J].Journal of Vibration and Shock,2011,30(9):150- 154.

[8] 胡宗文,袁国青.基于SPH模型的冰雹高速冲击复合材料结构响应规律研究 [J].力学季刊,2016,37(1):104- 113. HU Zong-wen,YUAN Guo-qing.Response regularity of composite structures under high velocity hail impact based on SPH model [J].Chinese Quarterly of Mechanics,2016,37(1):104- 113.

[9] HASHIN Z,ROTEM A.A fatigue failure criterion for fiber reinforced materials [J].Journal of Composite Materials,1973,7(4):448- 464.

[10] HE W,GUAN Z,LI X,et al.Prediction of permanent indentation due to impact on laminated composites based on an elasto-plastic model incorporating fiber failure [J].Composite Structures,2013,96(4):232- 242.

[11] BENZEGGAGH M L,KENANE M.Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus [J].Composites Science & Technology,1996,56(4):439- 449.

[12] KIM H,KEDWARD K T.Modeling hail ice impacts and predicting impact damage initiation in composite structures [J].AIAA Journal,2000,38(7):1278- 1288.

[13] KIM H.The damage resistance of composite structures to high velocity ice impacts and their tolerance to impact damage [D].San Diego:University of California,1998.

[14] ANGHILERI M,CASTELLETTI L M L,INVERNIZZI F,et al.A survey of numerical models for hail impact ana-lysis using explicit finite element codes [J].Internatio-nal Journal of Impact Engineering,2005,31(8):929- 944.

[15] 姚小虎,何艳斌,张晓晴.高速磁悬浮车辆复合板抗冰雹冲击动力响应研究 [J].北京理工大学学报,2013,33(2):116- 120. YAO Xiao-hu,HE Yan-bin,ZHANG Xiao-qing.Dynamic response of sandwich plates of maglev train under Hail impact [J].Transaction of Beijing Institute of Technology,2013,33(2):116- 120.

[16] FUNAI S.Hail ice impact on composite structures at glancing angles [D].San Diego:University of California,2012.

Dynamic Response and Damage Prediction of Composite Stiffened Panel Under Hail Impact

ZHANGXiao-qingDINGTieLONGShu-changYAOXiao-hu

(School of Civil Engineering and Transportation, South China University of Technology, Guangzhou 510640, Guangdong, China)

In order to discover the damage characteristics of composite stiffened panel under hail impact, firstly, a finite element model of composite stiffened structure is established in light of Abaqus. Secondly, the process of hail impact is simulated by means of smooth particle hydrodynamic approach. Then, the influences of impact location, impact energy and impact angle on the dynamic response of composite stiffened panels are analyzed. Moreover, the panel damage and the delamination between panel and stiffener are also discussed. Finally, a comparison between the numerical results and the experimental ones is made to verify the accuracy of the proposed finite element model. It is found that, among various damage forms, tensile damage of the matrix and delamination damage between the panel and the stiffener are most significant; and that, when the impact point is at the edge of stiffener, the delamination damage between the panel and the stiffener becomes the most serious.

composite stiffened panel; hail; impact; dynamic response; damage; finite element model

2016- 08- 17

国家自然科学基金资助项目(11472110,11372113);华南理工大学中央高校基本科研业务费专项资金资助项目(2015ZP042);广东省普通高校特色创新项目(2014ktscx015);爆炸科学与技术国家重点试验室开放基金资助项目(KFJJ15- 20M) Foundation items: Supported by the National Natural Science Foundation of China(11372113,11472100)

张晓晴(1974-),女,教授,博士生导师,主要从事复合材料力学和冲击动力学研究.E-mail:tcqzhang@scut.edu.cn

1000- 565X(2017)05- 0120- 09

TB 33

10.3969/j.issn.1000-565X.2017.05.017