0 引言

飞机机体结构的绝大部分零件通过铆钉/螺栓等紧固件连接组合在一起。目前用于飞机结构精细化仿真的有限元模型多采用NASTRAN的CWELD/CFAST单元，或者Abaqus的FASTNER单元模拟紧固件，单元属性通常定义为线弹性^[1]，刚度计算采用理论分析、工程经验公式、试验数据拟合和有限元等方法^[2-4]确定。在进行结构的极限承载能力分析研究时，单钉钉载通常会超过其弹性极限，从而导致在钉群中载荷出现重新分配的情况。若此时再按照紧固件线弹性刚度计算出来的钉载进行评估，得到的强度结论是非常保守的，因此需要在失效分析中引入紧固件的非线性刚度，才能对紧固件承载、塑性、破坏的过程进行较为准确的数值仿真。

国内外学者针对连接区域钉载分布进行了大量的研究。张震等采用了解析方法来计算钉载分布^[5]，Johan采用三维有限元模型研究了多钉单剪钉载分布^[6]，赵美英采用二维有限元模型研究了多钉载荷分配及影响因素^[7]，吴存利等采用三维有限元模型研究了铆钉p-δ曲线有限元非线性数值分析方法^[8]，但关于紧固件的非线性刚度以及渐进损伤对结构承载能力的影响的研究并不多。本文采用Abaqus软件建立了三钉单剪的三维有限元模型，考虑了连接件之间的接触及材料的弹塑性，通过非线性计算得到全载荷历程的紧固件刚度。通过将计算得到的弹性阶段的柔度值与理论分析、工程经验公式计算结果进行对比，已验证了该有限元方法的可行性^[9]。但采用Abaqus手工建模方法计算紧固件的非线性刚度，在结构复杂或有限元模型比较大的时候，还需要分别对各种不同连接组合进行分析，出现重复建模分析的情况，效率低下，不利于后续结构强度的评估。

本文基于Abaqus软件形成的参数化建模方法，将三维模型分析的紧固件刚度曲线数据应用到二维壳单元模型中，采用三钉单剪模型和实际结构试验数据进行验证。最后，采用数值方法对复杂结构破坏试验中发生的钉载弹塑性重新分布以及钉断裂所引起的渐进损伤进行了非线性数值仿真，有限元计算结果与理论预期以及试验数据吻合良好。

1 参数化计算分析紧固件刚度

1.1 参数化建模方法

飞机结构中存在多种结构/紧固件排列组合情况，工程分析对计算效率要求较高。在标准材料库和紧固件库的基础上，本文对有限元建模、计算以及刚度提取方法^[9]，运用Abaqus Python脚本语言编程，固化了紧固件刚度计算的方法和流程，解决了工程分析中大量紧固件非线性刚度计算的需求^[10]。参数化建模及计算分析流程如图1所示。

如图2所示，针对单剪搭接拉伸试样模型，建立了一套参数化建模及计算分析程序模板，采用1/2对称模型进行分析。分析过程中考虑材料的弹塑性及接触等非线性因素，忽略摩擦影响。其中紧固件的几何尺寸信息按照牌号可从标准库读取。对于铆钉，假设上、下墩头对称；对于螺栓，八角形钉头和螺母简化为圆形。

1.2 材料参数及边界条件

上下板、紧固件的材料弹性参数如表1所示，基本材料参数来源于MMPDS-06^[11]：

在下板的右端（X方向）施加拉力。对称面部位采用对称边界条件设置，平板两侧约束Z方向自由度，上板左侧约束所有平动自由度。非线性刚度曲线数据的提取原理与文献^[1]中一致。

2 精细模型中非线性钉刚度的引入

2.1 钉载分配

参照实际工程计算中使用的精细模型^[12]的建模方法，建立如图3所示的壳单元模型。铆接的两块板采用S4R单元^[13]，材料、载荷与边界条件与三维有限元模型相同。铆钉采用三个Fastener单元。

定义铆钉的Connector Section时，采用其中的bushing来模拟。bushing具有六个自由度，两个剪切刚度采用三维模型的计算结果，其他4个刚度按梁理论估算。

以三维模型为基准，对比简化后的壳元模型的钉载分配系数，如表2所示。可以看出简化后壳元模型与三维基准模型的结果基本一致。

2.2 紧固件弹塑性及渐进损伤模拟

将紧固件的弹塑性刚度曲线引入到非线性Fastener元中^[14]，同时，将该曲线的载荷上限设定为钉的许用剪力，采用刚度折减的方法模拟紧固件断裂后的载荷释放。三颗钉在高载（17 136 N）下的计算结果如图4所示。图中Rivet1为第一颗铆钉，Rivet2为第二颗铆钉，Rivet3为第三颗铆钉。

在60%载荷以前，三颗铆钉钉载线性增加；之后3号铆钉进入塑性，钉载增加明显滞后；当载荷增加到90%后，3号铆钉断裂并卸载，1、2号铆钉钉载迅速增加。钉载变化与预期一致。

3 复杂结构紧固件渐进损伤模拟

3.1 试验介绍

为研究某型号飞机龙骨梁部段结构的破坏损伤机制，曾进行过图5所示的结构静力破坏试验。试验件由龙骨梁两个盒段、后延伸段及部分机身蒙皮组成，盒段前端设计有加强盒，在后延伸段上的两个框上端分别设计有槽型铝型材加强。试验件结构及支持条件示意图如图5所示。图中Mx、My和Mz分别代表绕X、Y、Z坐标轴的扭转力矩。

现场试验件上端固定，下端通过作动筒进行加载，约束及加载系统如图6所示。

破坏后，试验件后延伸段上缘条和腹板发现失稳破坏，两侧上缘条和腹板对称向外侧鼓出，如图7所示。拆卸后发现多个铆钉及高锁螺栓断裂，最终试验件在103%极限载荷时破坏。

3.2 仿真结果及对比

参照文献^[12]中的建模方案对试验件进行精细有限元建模：结构件采用S4R壳单元建模，单元尺寸10 mm，盒段中间的圆形开口处采用B31梁单元模拟孔边加强，紧固件采用Fastener进行模拟。此外，在高应力区补充材料的弹塑性定义，试验件的精细有限元模型共包括55 902个单元，其中330个紧固件单元，含32个非线性紧固件单元，剪切刚度采用三维体元计算。

首先进行线弹性分析以验证模型。载荷达到40%极限载荷时应力如图8所示，可以看出结构的左右两侧的应力分布基本一致，重点关注部位龙骨梁后延伸段上、下缘条以及腹板的位置应力较大，与实测应变片的对比如表3所示，最大误差5%。放大变形效果后，可以发现后延伸段有向外侧翻的趋势，与最终屈曲失稳形式一致。

施加120%极限载荷，采用Standard求解器进行非线性计算。计算到52%极限载荷时，后延伸段上缘条迅速侧翻，结构失稳，局部应力迅速增加并进入塑性，实测失稳点为50%左右；计算到116%极限载荷时，由于过大的塑性变形导致求解发散，该理论承载能力与试验值103%相差13%。需要说明的是，模型中目前尚未引入壳单元的失效机制，试验中结构本身的撕裂未能模拟，这是计算承载能力高于试验值的主要原因。试验件在100%极限载荷下的应力结果和变形状态如图9所示。

铆钉的破坏过程是本次计算的重点，区后延伸段下缘条内侧与机身蒙皮连接的部分铆钉及高锁螺栓断裂情况见图10。分析中第一颗铆钉断裂发生在47%，位置为图10中的蓝色铆钉（钉载为零，显示断裂）右边第一颗，试验中第一次断裂声响发生在65%，第一颗钉断裂前附近应变片数值最高。分析在100%载荷下共有8颗铆钉断裂，与实测一致。这8颗钉的钉载变化过程如图11所示。图中Rivet1-L、Rivet2-L······Rivet8-L分别为龙骨梁后延伸段的8颗断裂铆钉。

分析中所采用的多项输入数据，如材料的屈服极限和紧固件的许用载荷，都是在大量试验数据基础上的统计值，无论是A基准还是B基准都偏保守，这造成计算所得铆钉断裂载荷早于实际结构。但是，本模型的分析，还是很好地模拟了结构屈曲失稳、钉群载荷重新分布、紧固件渐次断裂等复杂非线性现象，其价值不仅在于预测结构失效载荷的下限和失效位置，更在于揭示初始破坏发生后，内载重新分布，再平衡，再破坏，直至最终完全丧失承载能力的过程，为优化结构的细节设计提供了理论依据。

4 结论

本文结合三维实体和二维壳元有限元模型，采用非线性有限元方法，对飞机结构中大量使用的紧固件的建模与分析方法进行了研究，结论如下：

1）将完整的紧固件非线性剪切刚度曲线，应用于壳元模型的失效仿真当中；

2）利用Abaqus参数化建模方法，编制了一套计算紧固件非线性剪切刚度的有限元模型程序，满足了实际工程分析中大批量紧固件计算的需求，同时与简化后的壳元模型的钉载分布进行了对比验证；

3）综合上述成果，应用于实际复杂结构的失效分析当中，对结构屈曲失稳及紧固件渐进损伤等高度非线性过程进行了仿真分析，计算结果与试验所得吻合良好，为结构渐次损伤过程的深入研究提供了一套有效的数值方法。

参考文献