复合材料加筋板后屈曲分析方法及实验验证

林国伟,李新祥

（中国飞机强度研究所 全尺寸飞机结构静力/疲劳航空科技重点实验室，西安 710065）

复合材料加筋板具有比强度和比模量大以及可设计性等优点，被广泛应用于飞机机翼和机身结构。飞机结构上的加筋板大多数为中长柱范围内的加筋板，在压缩载荷作用下会出现局部屈曲，屈曲后结构仍有很高的承载能力，即后屈曲承载能力。复合材料后屈曲承载能力的计算方法主要是有限元分析方法和工程计算方法。由于复合材料加筋板后屈曲的破坏过程较为复杂，现有的研究大多基于有限元方法进行数值模拟，主要的研究方法有虚裂纹闭合技术[1]，界面单元技术[2-7]，渐进损伤分析技术[8-12]等。虽然上述大多数有限元分析研究计算得到的屈曲载荷和破坏载荷吻合良好，但是没有给出应变分布和实验结果的对比，缺乏对有限元模型有效性进一步的验证。本工作应用渐进损伤分析方法（PDA）和黏聚区模型（CZM）分别模拟复合材料蒙皮、筋条的损伤以及它们之间的界面失效，同时通过应变分布和实验结果的对比来验证模型的有效性。

工程计算方法是根据实验数据和经验修正系数总结归纳的简便计算公式，该方法优点是计算简便，缺点是存在计算精度较低的问题；所以，工程上迫切需要较高精度的计算方法。目前估算加筋壁板轴压承载能力的工程计算方法主要有两种，分别是有效宽度法和分段处理法。有效宽度法的关键是蒙皮有效宽度的计算，陈金睿[13]总结了国内外常见的有限宽度计算公式，这些公式主要借鉴金属加筋板的有效宽度计算公式，适用范围有限，所以该方法工程上应用较少。目前工程实践中大多采用分段处理法。该方法的核心是要事先获取短加筋板的平均压损应力，获取方法目前有三种：第一种是通过短加筋板的压损实验得到平均压损破坏应力，这种方法成本较高；第二种是近似取蒙皮层合板的纯压缩破坏应力0.75 倍的值[14]，其中蒙皮层合板不考虑冲击的影响，该方法不考虑筋条形式及厚度，只用蒙皮破坏应力近似代替，故计算精度较低；第三种是通过有限元计算确定短加筋板的压损应力[15]，该方法需要进行有限元建模，工作量较大。由于短加筋板的后屈曲历程较短[16]，往往在筋条出现局部屈曲后结构很快破坏，故本工作提出取筋条各个板元的最小局部屈曲应力，近似作为短加筋板的压损应力。当板元局部屈曲应力超过其压缩极限应力时，则取压缩极限应力作为短加筋板的压损应力。本工作对上述基于短加筋板压损应力改进取法的工程计算方法进行了实验验证。

1 实验

1.1 实验设计

试件构型为复合材料T 型长桁加筋板。加筋板外形尺寸如图1所示。长桁截面尺寸如图2所示。试件数量为4 件。蒙皮铺层：[（± 45）/0/−45/90/45/0/−45/0/45/0/−45/0/45/90/−45/0/45/]S，（± 45）为 织物CF3011/BA9916，其他为单向带，材料为CCF300/BA9916。单层板和织物的材料属性分别如表1和表2所示，材料性能数据由材料供应商提供。单向带固化后单层厚度为0.125 mm，织物固化后单层厚度为0.2 mm。试件长度方向为纤维0°方向。腹板铺层：[（± 45）/0/−45/0/0/45/90/−45/0/0/0/45/0/0/45/0/0/0/−45/90/45/0/0/−45/0]S。缘条铺层：[（± 45）/0/−45/0/0/45/90/−45/0/0/0/45/0/0/45/0/0/0/−45/90/45/0/0/−45/0]。长桁和蒙皮之间采用J116B 胶膜共固化，材料属性如表3所示。

图1 试件外形和尺寸Fig.1 Shape and sizes of specimen

图2 试件长桁截面尺寸Fig.2 Shape and sizes of specimen

表1 CCF300/BA9916 材料属性Table 1 Mechanical properties of CCF300/BA9916

表2 织物CF3011/BA9916 材料属性Table 2 Mechanical properties of CF3011/BA9916

表3 J116B 胶膜材料属性Table 3 Mechanical properties of J116B

实验在YY200A 压力试验机上进行，在试件端部沿长桁轴向施加压缩载荷，并在翼肋（距试件中面各350 mm）的位置上由刀口对试件提供支持，限制该处沿板厚度方向的位移，为防止试件受压两端压劈和保证压力分布均匀，在试件加工过程中已采用端部灌胶处理。

实验分两步进行。首先进行小载荷调试，测量应变，根据测量结果调整试件压心；然后进行正式实验，以70 kN 的加载级差加载至280 kN，接着以20 kN 加载级差加载至破坏，测量每一级载荷下的应变，观测并记录变形与破坏情况。

试件加载到1180 kN 后发出连续响声直至1300 kN 破坏（以第一件为例），试件中部靠上位置处整体折断，蒙皮撕裂。屈曲载荷和破坏载荷如表4所示。

2 有限元分析方法

2.1 渐进损伤分析方法

渐进损伤分析方法能够预测层合板的失效过程和破坏强度，该方法主要包括应力分析和破坏分析两个方面。首先建立有限元模型，给定外载和边界条件，在每个载荷增量步下通过非线性有限元分析确定各个单元的应力状态，然后判断单元是否满足失效准则，若满足，对单元进行刚度折减，继续增加载荷重复上述过程，直至结构出现大面积的单元失效，载荷无法施加。

2.2 失效准则

渐进损伤分析方法的核心是选取合适的复合材料失效准则。目前复合材料后屈曲的有限元分析多采用Hashin 准则或Puck 准则等基于失效模式的准则，但是目前没有研究表明，对于后屈曲分析，这些准则比其他准则能提供更高的计算精度。林国伟等[17]采用常见的不同失效准则对L 型加筋条后屈曲强度进行计算，并与实验对比，发现Hoffman准则[18]精度最高。本工作进行有限元分析时采用Hoffman 准则，其具体形式如下：

式中：XT为轴向拉伸强度；XC为轴向压缩强度；YT为横向拉伸强度；YC为横向压缩强度；S12为面内剪切强度；σ1为纤维方向应力；σ2为基体方向应力；τ12为面内剪切应力。

2.3 黏聚区模型

黏接元是一种基于黏聚区模型的零厚度界面元，由于不需要预制裂纹，特别适合用来模拟界面分层失效。因此本工作采用黏接元来模拟蒙皮和长桁之间界面的分层扩展。分别采用二次应力强度准则和二次能量释放率准则判断界面的损伤萌生与分层扩展，具体形式如下：

应力强度准则：

能量释放率准则：

式中：σI、σII和σIII为黏接元法向和两个切向的应力；σu，I、σu，II和σu，III为黏接元法向和两个切向的强度；GIC、GIIC和GIIIC分别为Ⅰ型、Ⅱ型和Ⅲ型断裂模式下的断裂韧度；GI、GII和GIII分别为Ⅰ型、Ⅱ型和Ⅲ型断裂模式下的能量释放率。

2.4 材料性能退化方案

当某个单元在分析过程中，应力状态满足失效准则，则根据材料性能退化方案对该单元的刚度进行折减。目前，准确给出失效以后单元的材料性能是很困难的。国内外研究者提出了很多材料性能退化方法，主要可分为两类：（1）一次性的刚度折减方案，即按照预先定义的方案将材料的工程弹性常数进行一次性的退化，一般是根据经验将其更新为原来的0 到0.01 之间；（2）连续性的刚度折减方案，通过定义损伤状态变量，并引入到单元的刚度矩阵来实现。损伤状态变量一般是应力或应变的连续函数，定义损伤状态变量时需要纤维方向或者基体方向的拉伸/压缩断裂韧性等材料属性，但是这些材料属性一般很难全部获得。

本工作采用第一种材料性能退化方案，满足失效准则时将工程弹性常数更新为原来的0.01。

3 有限元分析结果

3.1 有限元模型

采用ABAQUS 有限元分析软件进行有限元分析。长桁和蒙皮采用壳元（S4R）。长桁和蒙皮损伤演化过程通过Hoffman 失效准则判断结构中是否有单元出现损伤，若满足失效准则就将相应的单元材料属性折减为原来的0.01。长桁、蒙皮之间加入黏接元模拟脱粘分层失效。单元尺寸为10 mm。模型共计14580 个单元、18748 个节点。因为端部是钢框灌胶，支持性较强，为减少模型复杂度，忽略试件两端的端头，用边界位移条件代替，将试件底部固支，顶部施加轴向压缩位移载荷，并约束其他自由度。刀口支持部位约束面外位移。

由于屈曲对几何缺陷很敏感，而实际结构的几何缺陷又很难测定，所以在进行后屈曲的非线性分析时，可以使用一阶屈曲模态代替初始缺陷。分析过程分为两步。第一步，先进行特征值线性屈曲分析，将一阶模态的节点变形结果乘以扰动系数作为结构初始几何缺陷引入有限元模型中，扰动系数通常取厚度的1%；第二步，对该有限元模型进行非线性有限元分析，分析过程中考虑复合材料的失效。

利用配置为Intel E3-1231 8 核3.40 GHz CPU、32 G 内存的PC 机进行计算分析，耗时8 h 左右。

3.2 结果分析

图3是有限元分析得到屈曲模态和典型部位的应变片编号（背对背贴片，背面应变片编号为正面应变片编号加200，贴片方向均为沿试件长桁方向），长桁之间的蒙皮沿试件长度方向有6 个半波。图4和图5是蒙皮屈曲后波峰位置的应变片、以及波峰和波谷之间的反节点线位置的应变片的载荷-应变曲线。图6是试件长桁腹板中间典型位置的应变片的载荷-应变曲线。计算曲线和实验曲线在线性阶段均吻合良好，表明有限元模型能较好地模拟试件加载方向的刚度。载荷加到860 kN 左右时，图4的各曲线出现拐折，图5和图6各曲线也发生明显分叉，表明这些应变测量点处发生较大的面外位移，蒙皮和长桁腹板几乎同时出现了局部屈曲。计算曲线的分叉和拐折出现时间略晚于实验曲线。当结构屈曲后进入非线性状态，分析和实验曲线变化趋势基本一致，但存在一定偏差，这是由分析时几何缺陷的假设与实际并不完全一致，以及非线性分析求解过程中数值误差的累积等因素造成的。

图3 屈曲模态和典型部位应变片位置Fig.3 Finite element model and typical location of gauges

图4 20 号和220 号应变片位置的实验值和计算值对比Fig.4 Strain-load curves of gauge No.20 and No.220

图5 23 号和223 号应变片位置的实验值和计算值对比Fig.5 Strain-load curves of gauge No.23 and No.223

图6 70 号和270 号应变片位置的实验值和计算值对比Fig.6 Strain-load curves of gauge No.70 and No.270

表5和表6是加载到−800 kN 时蒙皮、长桁的应变实验值和计算值的比较。蒙皮的计算值和实验值吻合良好，而长桁数据误差略大一点，这是由于试件T 型长桁填充区（又称捻子条）填充的是单向带，而有限元模型不考虑填充区细节，直接用腹板铺层简化代替填充区的0°单向带，包含0°、± 45°和90°的腹板铺层在加载方向的刚度上小于铺贴方向均为0°的单向带，导致有限元模型长桁的刚度以及整体刚度较实际结构都偏小，在同样的载荷作用下反映到应变的计算值上就会比实验值偏大一些。表7是屈曲载荷实验值和计算值对比，有限元分析结果与实验值误差为5.9%，因为试件端部连接的是灌胶铁盒，其对试件的支持介于简支和固支之间，但计算分析时用固支边界条件来简化代替，其对试件的支持要强于实际情况，故计算得到的屈曲载荷要偏大一些。表8是破坏载荷实验值和计算值对比，有限元分析结果与实验值误差为1.2%。

表5 加载到−800 kN 时蒙皮应变的实验值和计算值对比Table 5 Comparison of skin strain between FEM results and experimental data（−800 kN）

表6 加载到−800 kN 时长桁应变的计算值和实验值对比Table 6 Comparison of stringer strain between FEM results and experimental data（−800 kN）

表7 屈曲载荷有限元计算值和实验值对比Table 7 Comparison of buckling load between FEM results and experimental data

表8 破坏载荷有限元计算值和实验值对比Table 8 Comparison of damage load between FEM results and experimental data

图7是有限元模型破坏情况，和实验结果一致，破坏位置在试件长度方向的中间靠上位置处。图8是试件破坏时长桁和蒙皮界面脱粘情况，可见界面脱粘轻微，不是引起试件破坏的主要因素。试件主要是长桁和蒙皮局部屈曲后变形过大引起的材料破坏。

图7 有限元模型破坏位置Fig.7 Failure location of finite element model

图8 胶层破坏位置Fig.8 Failure location of adhesive layer

4 工程计算方法

4.1 蒙皮局部屈曲载荷计算

四边简支正交各向异性层合板的轴压屈曲载荷Nxcr可以按经典理论公式[14]计算，计算时必须对板的纵向半波数m 求极小值。公式如下：

式中：a 是板元的长度；b 是板元的宽度；Dij是层合板的弯曲刚度系数（i，j=1，2，6）。

根据蒙皮板元宽度b 的三种可能取法（如图9，b1对应于筋间蒙皮净宽，b2对应于传统组合壁板筋间铆接线间距，b3对应于竖筋根部间距）及两种侧边支持情况（简支或固支）进行蒙皮局部屈曲计算。其中，b2=（b1+b3）/2。根据已有的研究[19-20]，板元宽度取b2，同时侧边支持取简支条件。

蒙皮屈曲应力σcr、蒙皮屈曲应变εcr和加筋板屈曲载荷Pcr计算公式如下：

如表9所示，按上述工程计算方法得到的屈曲载荷为801.4 kN，实验结果为865 kN，误差为–7.3%。

表9 屈曲载荷工程计算值和实验值对比Table 9 Comparison of buckling load between engineering calculation results and experimental data

4.2 加筋板承载能力计算

加筋板试件的有效长细比 L′/ρ为40.7，有效长细比的定义为，L 为加筋板长度，c 为加筋板受压端的端部支持系数，根据文献[14]推荐，取1.56。I 为加筋板的剖面对其形心轴的弯曲惯性矩。对于有效长细比20～60 之间的加筋板，可用分段处理法的中长柱计算公式进行计算，加筋板的平均破坏应力公式如下：

式中：σcr为蒙皮局部屈曲应力，由式（4）和（5）确定。

式（8）中σr为不考虑蒙皮或筋条局部屈曲后刚度降的影响，按欧拉公式计算的加筋板总体屈曲应力，由式（9）～（11）确定。

式中：G 为筋条中竖直腹板的等效剪切模量；λ 为形状系数，对于剖面形状不同的筋条，选用不同的λ 值；对于矩形剖面筋条，取λ=1.2；为筋条中竖直腹板的剖面积（腹板不垂直于蒙皮时，取垂直于蒙皮方向的投影值）；（EI）为加筋板相对于剖面中性轴的弯曲刚度；Pe为不考虑横向剪切效应的加筋板欧拉屈曲载荷；Pcr为考虑横向剪切效应的加筋板欧拉屈曲载荷。

通过经典层合板理论的计算极限强度的方法[21]计算各个板元的压缩极限应力。该方法的基础是计算每一层单层板的应力状态。根据各单层板的材料性能，计算层合板的刚度系数和柔度系数，然后求各单层材料主方向上应力和外载荷之间的关系，将各单层应力分别代入失效准则进行比较，确定哪一层先失效。如某一层满足失效准则，将该层的所有刚度系数设为零，仍保持其他单层板的几何位置，且失效层退化后整个层合板仍按经典层合板理论计算刚度系数和柔度系数。直至当层合板的各层全部失效，此时的载荷即为层合板极限强度。上述计算过程的层合板失效判断采用Hoffman 失效准则。通过该方法可以确定长桁各个板元的压缩极限应力，并标记为σf1，σf2，…，σfm。

按上述工程计算方法得到加筋板破坏应力，再乘以加筋板横截面积即可到破坏载荷为1246.6 kN，实验结果为1280.0 kN，误差为–2.6%。同时和文献[14]的短加筋板的压损破坏应力取法的计算结果进行比较，如表10所示，结果表明本工作提出方法误差更小。

表10 破坏载荷工程计算值和实验值对比Table 10 Comparison of damage load between engineering calculation results and experimental data

工程计算方法通过Pascal 语言编程进行计算，总用时0.4 h 左右，为有限元分析的1/20，计算效率提高明显。

5 结论

（1）有限元分析方法和工程计算方法误差均在10%以内，都能较准确地预测T 型加筋板的屈曲载荷和破坏载荷。

（2）工程计算方法比有限元分析方法的计算效率更高，计算用时仅需有限元分析方法的1/20，适合对T 型加筋板屈曲载荷和破坏载荷的快速预测。

（3）本工作提出的改进的工程计算方法比已有的方法精度更高，可供结构设计人员参考使用。