复合材料加筋板轴压后屈曲及承载有限元分析

王维阳，李 伟，曹奇凯

(沈阳飞机设计研究所，辽宁 沈阳 110035)

1 引言

随着轻量化技术的发展，复合材料结构在航空、航天、汽车、船舶等领域中得到了广泛的应用。长桁复合材料壁板由于较好的稳定性、抗冲击性能、较高的纵向荷载传递效率在以上研究领域中得到了国内外广泛关注[1-3]。较高的工程应用价值及这类结构缺陷的敏感性，使得这类结构在应用过程中的损伤失效问题显得尤为重要。加筋复合材料壁板常见的失效形式为屈曲，但其屈曲后依然表现出较强的持续承载能力。为了兼顾结构的安全性能及使用效率，加筋复合材料的屈曲失效及后屈曲持续能力问题成为这类结构设计及分析的热点和难点[4]。

由于缺陷及干扰的敏感性，这类结构后屈曲行为的复杂多变，使得不同缺陷及干扰条件下的屈曲及后屈曲失效行为较难预测，这也加大了复合材料整体加筋板屈曲及后屈曲试验的难度和成本[5，6]。由于结构的复杂性，使得问题理论分析难以得到有效结果。目前，在这类问题数值模拟研究方面，学者们通常围绕不同类型缺陷的复合材料加筋板的屈曲及后屈曲行为试验，采用CAE软件对结构的屈曲及后屈曲持续能力进行计算，比较分析不同失效准则的计算精度，同时构建加筋复合材料诸多失效评估方法，比如渐进损伤分析方法、粘聚区模型。在描述单层复合材料失效方面，常见的失效准则有Hoffman准则、Tsai-Wu准则、Hashin准则、Puch准则、Chang-Chang准则等[5-7]。数值模拟往往与试验结果较难吻合，即使吻合较好，其破坏机理难以确定，另一方面，自主软件的缺失，使得深入完备的失效机理比较分析面临困境。

本文在以往研究基础上[1]，拟针对帽形长桁复合材料壁板屈曲、后屈曲失效承载过程进行非线性理论分析，基于Hoffman失效准则，采用层合板理论构建出帽形长桁复合材料壁板屈曲及后屈曲计算方法和自主程序。对比现有方法验证分析本文方法的有效性，进一步针对复杂帽形长桁复合材料壁板屈曲及后屈曲路径进行追踪，并进行失稳后持续承载特性、承载能力分析。

2 复合材料层合板基本控制方程及数值模拟方法

2.1 空间折板本构关系及Hoffman失效准则

本文采用Mindlin板理论对面板和长桁进行统一建模。对于一般平板，Mindlin板几何关系如下

u=u0-z(φx-γzx)

v=v0-z(φy-γzy)

w=w0

(1)

φx=w，x，φy=w，y

(1a)

其中，u，v，w是任一点x，y，z方向的位移，u0，v0，w0是面板中面位移，γzx，γzy是横向剪切角，φx，φy是中面转角。

中面应变

(2)

曲率变形

κx=-w，xx+γxz，x

κy=-w，yy+γyz，y

2κxy=-2w，xy+γxz，y+γyz，x

(3)

图1 帽形长桁结构

对于长桁结构，将如图1所示长桁结构视为空间折板，将其处理为5个相交的平板。针对每个板的中面分别建立局部坐标系x(i)，y(i)，z(i)(i=1～5)。在各自局部坐标系中，板的几何关系同式(1-3)，但各变量应加上标(i)以示区别，例如

(4)

(5)

(6)

由位移转换关系可以得到折板单元刚度阵从局部坐标系向总体坐标系的转换关系

(7)

其中

(8)

K′e是局部坐标系中的单刚，若i是转换节点Ri=R(α，e由长桁局部坐标系与总体坐标系之间确定)，若i不是转换节点Ri=I7(7阶单位矩阵)。

面板和长桁满足各自的单层板本构关系。各层的应力-应变关系为

(9)

其中k是层的编号。

假定面板破坏由单层破坏引起。复合材料单层破坏可采用Hoffman准则[8，9]

(10)

其中，σL，σT，τLT是单层中沿纤维方向应力，垂直纤维方向应力和剪应力，Xt，Xc是纤维方向拉伸和压缩强度，Yt，Yc是垂直纤维方向拉伸和压缩强度，S是面内剪切强度。当时上式满足时，该层破坏、刚度降为0；若大量单层破坏，导致层板刚度大幅度下降，即可认为层板失效。

2.2 协调层板单元基本假定及控制方程

基于广义变分原理和拟协调方法[10]，采用拟协调三角形层板单元对帽形长桁结构进行离散，节点自由度为u0，v0，w，w，x，w，y，γzx，γzy。

(11)

其中ε*αβ，κ′*αβ，φ*α是试函数，可用节点位移表示线性应变、弯曲曲率和转角的离散变量，对位移自由度建立单刚，从而得到拟协调位移模式几何非线性三角形层板单元。

为了模拟实际边界在后屈曲变形中保持直线的特点，在面板边界引入节点自由度为6的空间梁元，因此边界节点自由度变为8

u0，v0，w，w，x，w，y，γzx，γzy，φ

φ是板的法线自转角，即边梁在板平面内的转角，通过罚因子(梁在面板中面内的弯曲刚度)使该转角趋于0。

2.3 屈曲及后屈曲数值模拟方法

系统的平衡条件是

δP=0

(12)

P=U-W是位能(W是已知外力的功)。上式给出在小有限转动条件下的后屈曲状态有限元平衡方程。为了求解非线性方程(12)，需要采用增量迭代法，上式的有限元增量迭代形式为

(13)

(14)

具体的基本计算步骤如下：

形成线性刚度阵K0和参考载荷向量F，求解方程K0η=F，给出线性解η0；

形成几何刚度矩阵K1(η0)，用反幂法求解特征值问题[K0+λK1(η0)]Δη=0，得到临界载荷因子λcr和屈曲形态η1；

用渐近后屈曲有限元分析，形成分支路径转换初值(λb，aη1)、切线刚度矩阵KT和不平衡力R，求解增量方程(13)，通过迭代计算后屈曲路径上的第1个解点；

在求出第i步位移和应力后，进行失效分析，确定本步各单元单层破坏情况和i+1步初值预测，修改切线刚度阵KT，迭代求解方程(13)，继续跟踪后屈曲路径，直到结构失效。

需要注意的是为了监视二次屈曲、发现和越过后屈曲路径上的奇异点，在每步收敛之后需要判断解路径上本载荷步的平衡稳定性。根据稳定性能量准则，可利用该点切线刚度阵(KT)LU分解结果，如果L阵对角元素皆为正，即

Lii>0

(15)

则该点是稳定的；若至少有一个元素为负值，则不稳定；如果相邻两个计算点，发生从稳定平衡到不稳定平衡的转变，则两点间必存在奇异点。

3 算例分析

【算例1】为验证所构建方法的正确性，将如图2所示工字截面长桁壁板(360×480)轴压后屈曲有限元分析，计算结果与平板/梁组合模型有限元分析结果符合，见表1。组合板模型中筋的下翼板对于面板局部屈曲的约束比梁更强，所以临界载荷较高。

图2 算例1 模型

表1 工字长桁壁板屈曲有限元分析结果

【算例2】 采用层板单元离散面板和长桁，整个面板作为一个区域。荷载及结构截面图如图3所示。面板铺层：16层(45/-45/0/-45/0/45/0/90)s。长桁铺层：11层 (45/-45/0/-45/0/90)s，单层厚度均为δ=0.125mm。面板和长桁材料为T300/BA9913。材料力学性能参数为E11=127GPa，E22=9.435GPa，ν12=0.305，G12=G13=4.680GPa，G23=3.300GPa，XT=1569，XC=1140，YT=54.7，YC=156，S=118(MPa) 。

如图4所示，分别采用稀疏网格和密集网格对结构进行离散，其中稀疏网格模型包括5472个板元(其中长桁2736个)、148个梁元(边梁)、2262个节点，15982个自由度。密集网格模型包括8600个板元(其中长桁4200个)、188个梁元(边梁)、3366个节点，2259个坐标转换节点，23750个自由度。

为模拟实际边界，假设后屈曲时边界保持直线，在面板四边引入边梁，加载边边梁单元在板面内的弯曲刚度(EJ′)很大，作为罚因子，取EJ′=1010，其余刚度为0。边界条件为总体简支和固支两种，令板中点A和y=0边中点B点的位移uA=vA=uB=0。

图3 截面形状、荷载及结构分布图

图4 结构网格离散

首先，针对不同约束进行结构屈曲分析，屈曲载荷在表2中给出。从计算结果可以看出计算简支屈曲载荷稀网格已经满足要求，而固支情况由于屈曲形态变化更加剧烈，需要采用密网格。

表2 屈曲载荷Pcr

图5和图6分别给出稀疏网格模型简支和密集网格固支屈曲形态。从图5结果图可以看出，简支边界时壁板为总体屈曲(上凸)，纵向为1个半波、横向为1个半波与4个微小半波之和，在加载端附近面板波形幅度增加较快，在中部面板和长桁单向上凸。由长桁屈曲形态可见，受面板转角(w，y)的影响，帽形长桁两个侧板外凸显著。由图6结果图可见，固支边面板的屈曲变形遍及大部分板，长桁间面板幅度较大，剧烈起伏，这也是稀网格临界载荷误差大的原因，长桁波形与面板协调，幅度略小。

图5 简支约束条件下稀疏网格屈曲形态图

图6 固支约束条件下密集网格屈曲形态图

从分支点(λcr=33.465)开始跟踪后屈曲路径。图7为载荷-位移路径图，其路径为0-1-2-3-4-5-6。0是分支点(λ=33.465)。图8和图9为典型的两个荷载步19步(λ=31.35)和75步(λ=58.73)时单层应力分布图。

图7 面板中点挠度随荷载变化图

图8 单层失效发生后长桁顶部和面层R分布

图9 第75荷载步长桁顶部和面层R分布

路径1-2(向长桁一侧凸出状态)是不稳定的，对于完善(无缺陷)结构实际上不能发生；路径2-3(向面板一侧凸出状态)是稳定的，并且未发现有其它邻近的稳定平衡状态，因此可以认为路径2-3能够发生。实际上，受缺陷影响，屈曲前的基本状态会偏离理想路径上升，并转向路径2-3，在理想情况下可达到载荷极值λs=63.97(Ps≈640kN)，壁板发生不稳定失效。在特定缺陷影响下，可能发生3根长桁同时向上弯曲的后屈曲波形，即第7步之前、与屈曲形态类似的变形，此时的承载能力应为λ′s≈33。

从图8也即是单层失效发生后的第2步R的分布情况可见失效单层只在长桁2#的顶部各层出现，共628个单元层。图9为结构接近失稳破坏的高载荷下应力水平参数及单层破坏分布，情况与图8类似，破坏单层只在2#长桁的顶部，仍无各层均破坏的单元，但面板应力水平显著提高(R=0.309)。

从后屈曲损伤计算结果可以看出，壁板失效前未发生强度破坏。3根帽形长桁加强壁板可以有效地提高局部屈曲载荷，降低弯曲应力水平、改善结构强度。但后屈曲行为十分复杂，伴随波形的突然变化，对缺陷敏感程度可能加大。壁板在较高载荷下发生总体屈曲和失稳破坏，预测的临界压力约为Pcr≈335kN，极限承载能力约为压力Ps≈640kN。受缺陷因素影响承载能力可能更低一些，在特定缺陷情况下甚至低得多、与临界载荷相当。

4 结论

本文基于Hofman失效准则和层合板理论构建了帽形长桁加筋结构板屈曲及后屈曲分析方法。通过工字长桁壁板算例验证了方法的有效性，同时针对典型的3根帽形加筋板壁屈曲临界荷载和后屈曲行为进行了求解，预测了临界压力和极限承载能力。计算结果表明帽形长桁加强壁板可以有效地提高局部屈曲载荷，降低弯曲应力水平、改善结构强度，另外，帽形长桁壁板后屈曲失效发生在中间长桁顶部，壁板失效前未发生强度破坏，具有较大的持续承载能力。本文方法的构建为复杂的长桁加筋结构板屈曲临界荷载预测及后屈曲承载能力预测提供了有力的理论及软件技术支持。进一步为类似复杂加筋结构设计和优化提供可靠的数值依据，可以有效降低结构设计的试验成本，并提高结构设计的效率。