薄板的后屈曲损伤分析与疲劳寿命预估

肖浩，胡伟平*，张淼，孟庆春

(1.北京航空航天大学 航空科学与工程学院，北京100191;2.中国空间技术研究院，北京100094)

在航空航天、管道运输等工程领域，结构的稳定性问题一直受到重点关注.结构的屈曲问题早在18世纪由Euler和Lagrange提出，在这之后，Duncan［1］对不均匀柱进行了研究，为屈曲临界载荷推导了封闭形式的表达式.此后，对于封闭解的研究较少.1999 年，Elishakoff和 Rollot［2］才发现了新的封闭解.但工程结构的实际屈曲载荷明显要低于由理论解给出的屈曲临界载荷，当时的理论无法解释这一现象.

1941 年，von Kármán 和 Tsien［3］根据非线性大挠度理论，提出了板壳后屈曲分析的一般方法以及非线性屈曲理论.Donnell和 Wan［4］于1950年将非线性大挠度分析推广到非完善圆柱薄壳.Koiter［5］用摄动法研究了弹性结构的初始后屈曲性态，导出了临界压力与缺陷参数之间的关系，并提出了初始缺陷敏感度的概念.1968年，Stein［6］提出了非线性前屈曲一致理论，解释了经典线性理论与实验之间差异的原因.非线性屈曲分析考虑了结构屈曲前的变形，计算结果更符合实际情况.

以von Kármán大挠度理论为依据，对于单一的屈曲和后屈曲问题已有大量的研究成果，但在实际工程应用中，当薄壁构件受到循环载荷作用时，可能会出现屈曲问题与疲劳问题的耦合［7-8］，即结构可能在反复进入后屈曲状态下发生疲劳破坏，或是在产生一定的疲劳损伤累积后发生屈曲或后屈曲破坏.

目前，国内外很多学者主要研究构件在循环载荷下屈曲后屈曲过程中模态的改变［9-10］，以及临界载荷和极限载荷随着循环次数的增加而产生的变化［11-12］，很少有研究将屈曲分析与疲劳损伤分析进行耦合，并进一步预估构件的疲劳寿命.

本文基于连续损伤力学理论与方法，研究结构在后屈曲情况下的损伤以及在疲劳载荷作用下屈曲与疲劳损伤的耦合特性.首先采用有限元方法，通过线性屈曲分析得到屈曲临界载荷和屈曲模态，进而采用大挠度理论，将线性屈曲的一阶屈曲模态作为初始扰动，进行薄板的非线性屈曲分析，得到屈曲临界载荷.其次，根据损伤力学理论与方法建立薄板材料在单次加载过程中的损伤演化方程以及参数识别方法，并根据非线性屈曲分析结果进行后屈曲损伤分析.最后，考虑疲劳载荷的作用，基于损伤力学理论，采用有限元数值方法进行求解，考虑每次加载引起的疲劳损伤与后屈曲应力应变场的耦合作用，通过反复迭代计算，给出结构疲劳寿命.本研究为结构在后屈曲情况下的疲劳特性分析提供了一种新的分析方法和实现途径，可以应用于飞行器中某些特殊受力状态下薄壁构件的疲劳寿命分析.

1 薄板的线性屈曲分析

板的线性屈曲分析是一种最常规的分析方法，可以给出结构的各阶屈曲临界载荷以及相应的屈曲模态.

1.1 薄板的有限元模型

考虑一四边简支且在某两对边受均匀压力的方形薄板，板边长a=40mm，板厚h=0.6mm，材料为LC9CgS1［13］，其几何形状及载荷如图1所示.

图1 结构几何模型Fig.1 Geometry model of the structure

载荷大小为Px=40 N/mm，四边简支板的边界条件为

式中w为结构变形挠度.

在ANSYS中建立薄板的有限元模型.选择Shell181单元，该单元适用于薄到中等厚度的壳结构，单元有4个节点，每个节点有6个自由度，适用于结构的几何非线性分析.

1.2 特征值屈曲分析

建立有限元模型后，在ANSYS软件中进行特征值屈曲分析.计算得到薄板的一阶屈曲临界载荷为34.6 N/mm.

通过特征值屈曲分析能够给出结构屈曲模态，能为结构设计提供参考，并为后续的非线性屈曲分析提供初始扰动.

2 板结构的非线性后屈曲分析

在工程实际中，由于线性屈曲分析忽略了屈曲前变形的影响，导致过高地估计了结构的临界载荷.因此，为了提高分析的精度，更多地采取非线性屈曲分析方法进行结构稳定性分析.

2.1 非线性屈曲分析方法

非线性屈曲分析的目的是求解结构从稳定平衡过渡到不稳定平衡的临界载荷，结构的载荷临界点可以分为两种类型:分叉临界点和极值临界点［14］.

分叉临界点的特征是:结构在基本载荷-位移平衡路径Ⅰ的附近还存在另一分叉平衡路径Ⅱ.当载荷达到临界值Pcr时，如果结构或载荷受到一微小扰动，载荷-位移曲线将沿分叉平衡路径发展，其载荷-位移曲线如图2(a)所示.

图2 屈曲载荷临界点Fig.2 Critical point of buckling load

极值临界点的特征是:当载荷达到临界值时，如果载荷位移有微小变化，将分别发生位移的跳跃和载荷的快速下降，其载荷-位移曲线如图2(b)所示.

利用有限元软件进行结构的非线性屈曲分析，需要在施加沿板面内压缩载荷的同时，对结构施加一横向小扰动，通过失稳点的载荷-位移曲线计算屈曲临界载荷.

2.2 非线性屈曲分析结果

以特征值屈曲分析得到的第一阶屈曲模态作为初始几何扰动进行非线性屈曲分析，得到失稳点的离面位移随外载荷的变化曲线，如图3所示.从该曲线可以判断出非线性屈曲的临界载荷为31.5 N/mm.当外载荷达到非线性屈曲临界载荷时，结构进入后屈曲，屈曲前和屈曲后对应结构的两个不同的平衡路径，平衡路径的转换也是因为结构上微小扰动的缘故.图4为结构失稳部位的von Mises应力随外载荷的变化曲线，可以看出，微小扰动在前屈曲阶段对结构的应力影响很小，当外载荷达到屈曲临界载荷后，结构发生屈曲并进入另一个平衡状态，此时结构的应力随外载荷增加而显著增大.

图3 失稳点离面位移随外载的变化曲线Fig.3 Variation curve of deflection of instability point with external load

图4 失稳点的von Mises应力变化曲线Fig.4 Variation curve of von Mises stress of instability point

3 板结构的后屈曲损伤分析

通过板的后屈曲分析可知，当板所受载荷超过屈曲临界载荷时，板进入后屈曲状态，此时结构失稳部位的应力迅速增加.当板承受最大载荷超过屈曲临界载荷的循环载荷时，板在每一次载荷循环中都将经历一次后屈曲，因此为了预估板的疲劳寿命，首先需要进行板在后屈曲状态下的损伤分析.下面拟采用连续损伤力学理论与方法，建立结构在后屈曲状态下的损伤分析模型.

3.1 损伤演化模型

根据连续损伤力学理论，在疲劳载荷作用下，材料的损伤将引起材料有效承载面积的下降.在各向同性损伤的情况下，引入损伤度D描述这一物理过程［15］，令

在单向受力状态下，建立材料单元的力平衡方程:

对于受损伤的材料单元，其真实应变为

在三维受力情况下，对于各向同性材料，根据应变等价原理，得到含损伤线弹性体本构关系为

式中，λ和μ为拉梅常数;σij和εij分别为应力张量和应变张量分量.

单轴受力情况下，以应变为基本变量，含损伤材料应变能密度为

根据热力学原理，引入损伤驱动力［16］:

单轴受力情况下，有

对处于复杂受力状态下的结构而言，引入损伤力学等效应变:

式中εm为平均应变.

采用如下形式的损伤演化方程［17］:

式中，Ymax为损伤力学驱动力峰值;Yth为损伤力学驱动门槛值;N为载荷循环次数;mk为材质参数.式(12)可以写成如下应变形式的表达式:

式中，β和εth为材质参数;εmax为最大损伤力学等效应变.式(13)需要通过材料标准件疲劳试验结果识别.

3.2 参数识别

3.1节已经建立了材料在疲劳载荷作用下的损伤演化方程，方程中包含有β，mk，εth等材质参数，这些参数需要根据材料标准件的疲劳试验结果来确定.

对于光滑试样，根据单轴受力的本构关系可得

式中，D0，1为应力集中系数KT=1试件的初始损伤;σ0，max为交变载荷引起的最大应力;σ1，th为KT=1时材料发生疲劳损伤的应力门槛值.

对损伤演化方程积分，可得

式中，Nf为裂纹萌生寿命;m为材质参数.式(15)可以写为

式(16)即为KT=1时材料的裂纹萌生寿命和最大应力之间的关系.根据这一应力-寿命表达式，结合材料的疲劳试验结果，就可以确定式(16)中的材质参数 β，mk，D0，1及 σ1，th.LC9CgS1材料的静力力学性能如表1所示.当KT=1，应力比R=0.1时，LC9CgS1材料标准件的疲劳寿命如表2所示.

表1 LC9CgS1材料的静力力学性能参数Table1 Static mechanical properties of LC9CgS1

表2 KT=1，R=0.1时 LC9CgS1 材料的疲劳试验结果Table2 Fatigue test results of LC9CgS1 when KT=1，R=0.1

根据以上疲劳试验结果并假设实验结果中离中值寿命曲线最远点的实验点对应的材料初始损伤度为0，可以得到损伤演化方程中的材质参数如表3所示.

表3 KT=1，R=0.1时LC9CgS1的材质参数Table3 Material parameters of LC9CgS1 when KT=1，R=0.1

计算结果与实验点的对比如图5所示.

图5 KT=1，R=0.1时计算结果与实验点对比曲线Fig.5 Comparison of experimental points and calculation results when KT=1，R=0.1

从图5可以看出KT=1，R=0.1时拟合的材质参数是合理的.

4 板结构疲劳寿命预估

当板受到的最大载荷超过屈曲临界载荷的疲劳载荷作用时，每一次载荷循环中，板都要进入后屈曲状态，造成局部材料的损伤，这一损伤又会降低下一次加载时板的屈曲临界载荷，这种屈曲与疲劳损伤的耦合特性难以用闭合形式的方程描述.因此，基于损伤力学理论，采用有限元数值解法，通过反复迭代的后屈曲分析和疲劳损伤分析来预估板的疲劳寿命.

4.1 损伤力学-有限元数值解法

利用ANSYS软件中的APDL语言控制求解过程，将ANSYS软件计算得到的应力场代入损伤演化方程，计算结构的损伤场，给出各单元损伤度;根据单元损伤度计算各单元的等效弹性模量，在ANSYS软件中给各单元重新赋值，再次计算下一个载荷作用下的应力场.如此迭代计算，在ANSYS平台中，基于损伤力学理论，采用有限元数值解法对构件进行疲劳寿命预估.

具体做法如下:

1)利用ANSYS软件计算结构在某一载荷循环下的应力应变场.

2)根据上一步计算得到的应力应变场，通过APDL语言对ANSYS软件进行二次开发，计算每个单元的损伤力学等效应力.

3)利用APDL语言，在ANSYS软件中计算每个单元在ΔN次载荷循环之后，损伤度的增量:

式中，βk为材质参数;εe，e，max为单元的最大损伤力学等效应变;εk，th为损伤力学应变门槛值.

那么，可以得到当前各个单元的损伤度:

4)由上述计算得到的各个单元的损伤度可知每个单元的等效弹性模量的变化量，在ANSYS软件中重新赋予各个单元弹性模量，在上次计算结果基础上，进行下一个加载循环的应力应变分析.

5)重复步骤2)～步骤4)的过程，直到某一单元损伤度达到或超过1，此时的ΔN的累积值即为疲劳裂纹萌生寿命.

需要说明的是，载荷循环次数步长ΔN需通过不断减小ΔN的值来进行收敛性验证.

损伤力学有限元法计算流程如图6所示.

图6 损伤力学有限元法计算流程图Fig.6 Calculation flow chart of damage mechanics with finite element method

4.2 计算结果

给定板受到的疲劳载荷为:最大载荷为40 N/mm，载荷比 R=0.1.由于最大载荷大于板的屈曲临界载荷31.5 N/mm，因此每一次加载循环中板都将进入后屈曲阶段.

薄板在后屈曲状态下，失稳部位局部发生损伤，这一损伤会导致板的总体刚度矩阵发生变化，从而影响下一次加载时的屈曲临界载荷以及结构位移.图7显示了不同损伤阶段失稳部位的载荷位移曲线的变化.从图7可以看出，板的屈曲临界载荷随着损伤的增加而减小，失稳点的离面位移随着损伤的增加而增大.

图7 不同损伤度下结构失稳点位移变化曲线对比Fig.7 Comparisons of deflection changing curves of instability points in different damage degrees

通过计算，板在上述疲劳载荷作用下的疲劳寿命为6×105次，失效的单元位于失稳点，失稳点周围的区域是损伤比较大的区域，因为随着损伤的增大，结构的屈曲临界载荷会逐渐变小，失稳区域会扩大，这一部分区域损伤也会因此而变大.

4.3 初始几何缺陷对结构裂纹萌生寿命的影响

在前面考虑后屈曲的疲劳寿命分析中，为了模拟屈曲这一物理过程，在每一次屈曲分析时都对结构施加了扰动.

采用的施加扰动方法是引入初始几何缺陷，即将第一阶线性屈曲位移模态作为初始缺陷施加在板上.然而这里存在一个问题，即作为初始缺陷的一阶线性屈曲位移模态的幅值如何选取，不同幅值对板的疲劳寿命的影响程度如何.

因此计算了不同缺陷幅值比(幅值比=几何缺陷幅值/板厚，无量纲)情况下板的疲劳寿命，如图8所示.

图8 裂纹萌生寿命随几何缺陷幅值比的变化曲线Fig.8 Changing curve of crack initiation life with different geometric imperfection amplitude ratios

从图8中可以看出，初始几何缺陷的大小对板的疲劳寿命影响很大，因此在预估实际结构疲劳寿命时，应首先测量结构的几何外形，根据测量结果尽量准确估计结构的初始几何缺陷，以便得到更加准确的疲劳寿命分析结果.

5 结论

基于损伤力学理论构建了薄板后屈曲状态下的裂纹萌生寿命预估方法，为工程结构考虑后屈曲损伤的疲劳寿命分析提供了一种新的研究方法和手段.

1)首先对板进行了线性屈曲分析，得到了屈曲临界载荷以及屈曲位移模态，为后续的非线性屈曲分析提供了位移扰动形式.

2)以线性屈曲分析得到的一阶屈曲位移模态作为板的初始几何缺陷，对板进行非线性屈曲分析，得到了板的非线性屈曲载荷，其值要低于特征值屈曲分析得到的屈曲临界载荷.

3)建立了板后屈曲损伤分析模型，并根据板材料的标准件疲劳试验结果确定损伤演化方程中的各材质参数.

4)建立了考虑后屈曲过程的板的疲劳寿命分析方法，通过损伤力学-有限元数值解法，可以考虑加载循环过程中板的后屈曲与疲劳损伤的耦合作用.

5)分析初始几何缺陷的大小对板裂纹萌生寿命的影响，结果表明，初始几何缺陷对结构的疲劳寿命有很明显影响，因此，一方面在制造过程中应尽量提高加工精度，减小几何缺陷;另一方面，在预估结构疲劳寿命时应根据实际测量结果正确估计几何缺陷的大小，以便给出更符合实际情况的寿命预估结果.

References)

［1］ Duncan W J.Galerkin’s method in mechanics and differential equations，ADA951862［R］.London:Aeronautical Research Council，1937.

［2］ Elishakoff I，Rollot O.New closed-form solutions for buckling of a variable stiffness column by Mathematica［J］.Journal of Sound and Vibration，1999，224(1):172-182.

［3］ von Kármán T，Tsien H.The buckling of thin cylindrical shells under axial compression［J］.Journal of the Aeronautics Science，1941，8(6):303-312.

［4］ Donnell L H，Wan C C.Effect of imperfections on buckling of thin cylinders and columns under axial compression［J］.Journal of Applied Mechanics-transactions of the ASME，1950，17(1):73-83.

［5］ Koiter W T.The stability of elastic equilibrium，AD0704124［R］.Standford:Department of Aeronautics and Astronautics of Standford University，1970.

［6］ Stein M.Some recent advances in the investigation of shell buckling［J］.AIAA Journal，1968，6(12):2339-2345.

［7］ Kong C，Bang J，Sugiyama Y.Structural investigation of composite wind turbine blade considering various load cases and fatigue life［J］.Energy，2005，30(11):2101-2114.

［8］ Rhead A T，Butler R，Hunt G W.Post-buckled propagation model for compressive fatigue of impact damaged laminates［J］，International Journal of Solids and Structures，2008，45(16):4349-4361.

［9］ 李世荣，刘平.弹性地基上Euler-Bernoulli梁的热屈曲模态跃迁特性［J］.应用力学学报，2011，28(1):90-94.Li S R，Liu P.Thermal buckling mode transition characteristics of Euler-Bernoulli beam on the elastic foundation［J］.Journal of Applied Mechanics，2011，28(1):90-94(in Chinese).

［10］ 贺尔铭，张钊，胡亚琪，等.单轴压缩载荷下夹层梁结构屈曲及皱曲模拟研究［J］.西北工业大学学报，2012，30(5):668-674.He E M，Zhang Z，Hu Y Q，et al.Buckling and wrinkling study of sandwich beam under uniaxial compression load［J］.Journal of Northwestern Polytechnical University，2012，30(5):668-674(in Chinese).

［11］ 左志亮，蔡健，朱昌宏.带约束拉杆L形钢管混凝土短柱的偏压试验研究［J］.东南大学学报:自然科学版，2010，40(2):346-351.Zuo Z L，Cai J，Zhu C H.Bias experimental study on the binding L-shaped concrete bars filled with steel tubes［J］.Journal of Southeast University:Natural Sciences，2010，40(2):346-351(in Chinese).

［12］ 夏盛来，何景武，海尔瀚.大展弦比机翼屈曲及后屈曲分析［J］.机械强度，2012，33(6):907-912.Xia S L，He J W，Hai E H.Buckling and post-buckling analysis of the high aspect ration wing［J］.Mechanics Strength，2012，33(6):907-912(in Chinese).

［13］ 吴学仁.飞机结构金属材料力学性能手册:静强度疲劳/耐久性［M］.北京:航空工业出版社，1996:225-353.Wu X R.Manual of metal material mechanics performance of the plane structure:static fatigue strength/durability［M］.Beijing:Aviation Industry Press，1996:225-353(in Chinese).

［14］ 王勖成.有限单元法［M］.北京:清华大学出版社，2003:649-651.Wang X C.The finite element method［M］.Beijing:Tsinghua University Press，2003:649-651(in Chinese).

［15］ Lemaitre J，Desmorat R.Engineering damage mechanics:ductile，creep，fatigue and brittle failures［M］.Berlin:Springer，2005:3-6.

［16］ 张行，崔德渝，孟庆春.断裂与损伤力学［M］.北京:北京航空航天大学出版社，2006:329-331.Zhang X，Cui D Y，Meng Q C.Fracture and damage mechanics［M］.Beijign:Beijing University of Aeronautics and Astronautics Press，2006:329-331(in Chinese).

［17］ 张淼.构件疲劳寿命预估的多变量损伤力学方法［D］.北京:北京航空航天大学，2009.Zhang M.Multivariate fatigue damage mechanics methods of fatigue life prediction of the structure［D］.Beijing:Beijing University of Aeronautics and Astronautics，2009(in Chinese).