0 引言

在飞机结构设计中，结构静强度是飞机结构设计成功的基础，在满足性能指标的基础上实现对结构更优的轻量化是飞机结构工程师不变的追求。飞机结构多采用加筋薄壳的设计，大量采用紧固件作为零件之间的机械连接方式，故连接件的强度分析是机体结构设计的重要组成部分。据统计，飞行器结构中有70%以上的破坏是由紧固件连接失效导致^[1-2]。钉载的研究始于NASA对金属双剪连接弹性段的研究，Rosenfeld^[3]据此建立了螺栓载荷的递推公式。在载荷作用下，大多数紧固件呈现出复合受载的情况，目前对紧固件刚度的主要研究方法有理论分析、工程经验公式和有限元等^[5]。

针对有限元方法，前人进行了一些研究。Gray^[6]等提出来用壳单元和三维单元混合的螺栓连接建模方法，以此获得螺栓内部的载荷分布情况以及真实的孔边应力状态。宋恩鹏^[7]等人通过试验和有限元仿真结果进行对比，提出了紧固件刚度会影响到多钉双剪连接的钉载分配，同时钉载的分配比例受到连接板材料刚度的影响。Ekh J^[8]等应用ABAQUS中的不同单元组合来模拟单排多钉的紧固件连接模式，在充分考虑接触非线性的影响下，通过与试验和三维细化模型的仿真结果进行对比，给出了单排多钉连接的单元选择建议。赵群^[9]等研究结构刚度对钉载分布的影响，得出适用于某型号飞机设计的钉载分布规律。

本文以某型飞机翼身组合体机翼破坏试验为研究对象，试验如图1所示。团队前期开展了一些翼身组合体破坏试验的仿真，为进一步提高仿真结果的精度，结合江蕊伶等^[10]的研究内容，将紧固件非线性对仿真结果的影响加以考虑。通过对紧固件刚度曲线的优化，将紧固件的刚度曲线引入翼组盒段失效仿真过程，得到了更为准确的金属机翼的承载能力评估结果。

1 机翼盒段有限元仿真

参照文献^[11]中的方法建立某型飞机机翼盒段的精细模型，精细有限元模型是在细化模型的基础上进行进一步的细分网格，以“全机打钉”的策略取代厚度融合的方法，建模思路更加直观。模型的网格划分更加精细，将结构更多的几何特征在模型中表现出来，实现对连接件的精确表达，避免了等效刚度简化对模型紧固件和结构件的省略。机翼盒段模型如图2所示，用该模型进行2.5 g载荷下的金属机翼承载能力计算。该模型主要由23万个壳单元和3 770个紧固件单元组成，节点数31万，总自由度数164万。在精细模型中，紧固件通过对连接单元添加刚度曲线实现对紧固件的模拟。

在ABAQUS中，紧固件有两种常见的表达方式。两种方式都采用连接单元对钉进行表达，具体差异在于紧固件头与板壳之间连接方式的表达。第一种方式是FASTENER，其特点在于通过COUPLING将紧固件端点与壳单元的节点进行连接，其连接范围通过设置影响半径来调整，如图3所示。第二种设置方式直接采用COUPLING进行设置，该方法属于离散定义，在连接单元节点与对应的壳面之间显示的定义COUPLING，当连接的壳单元范围控制在一个单元时，与NASTRAN中的CWELD单元表现基本一致，如图4所示。

两种数值表达中，都使用连接单元来表达紧固件钉杆，常见的连接单元有CARTESIAN-CARDAN和BUSHING。

在弯矩表达方面，两种连接单元的分配策略也有不同。BUSHING将单元弯矩平均分配到两个结点上，而CARTESIAN-CARDAN将单元弯矩全部算在第一结点上。对于紧固件连接壳单元的情况来说，连接单元的两个结点均通过将转角位移分配到受约束的从结点，进而模拟弯矩作用，因此连接单元内弯矩的处理对整体计算没有影响。

但BUSHING在计算时存在收敛性问题，主要源于其对转角自由度的定义方式。它的单元坐标轴由两个结点各自的坐标轴通过数学运算来建立，结点处的微小转动可能引起单元坐标轴的转角在正负之间突变，由此造成的弯矩突变将导致收敛困难。由于BUSHING单元本身数学模型的原因，收敛性比较差，故本文采用CARTESIAN-CARDAN作为表达紧固件的连接单元。

试验中和试验后的观测都显示，翼盒结构的最终破坏涉及到几何非线性（受压屈曲）、材料弹塑性以及接触等多个非线性现象，有限元仿真是一项极具挑战性的工作。为精确模拟紧固件的非线性传力特点，包括进入塑性段后的刚度变化和超过极限后的断裂（卸载）现象，本文根据江蕊伶等^[10]的研究内容，利用ABAQUS建立如图5所示三维单剪连接模型来计算钉的刚度曲线。下板采用固支约束，上板施加拉力，使用通用接触算法加入接触，忽略摩擦，对孔周边网格进行细化以提高计算精度。通过仿真可以得到紧固件的刚度曲线，如图6所示。

图6曲线上的点为模型每个收敛步结果计算所得的紧固件剪力和变形。补充紧固件的许用值信息后，图6所示的非线性刚度曲线就可以带入精细模型中模拟紧固件在高载下的塑性和断裂卸载特性，从而完成结构承载力的仿真计算。但是实际的计算表明，非线性迭代在远小于试验破坏载荷时，大约四分之三，即无法收敛，非线性仿真失败。

采用非线性有限元定义结构的弹塑性刚度曲线，多数软件要求曲线原点后的第一个数值点为弹性段的最后一个点。上述紧固件刚度曲线的弹性段定义了过多的中间数值点，导致数值迭代过程发散，删除“中间点”即可解决该数值问题。处理方法虽然简单，但是由于结构模型中包含大量的紧固件，弹塑性曲线的“清理”工作人工完成效率太低，需要建立一个可靠的算法由程序自动完成。

2 中点误差的计算方法

考虑曲线上任意连续的三个数据点$\left(x_1，y_1\right)，\left(x_2，y_2\right)，\left(x_3，y_3\right)$，上述清理中点的数学表述即是以P₁和P₃连线上的线性插值点P′₂近似点P₂。中点误差定义为e:

对于已有曲线上的数值点，除首点和尾点外，计算每个点的插值误差，当误差小于门槛值，即认为该三点线性度良好，中间点可删除。上式中的e有清晰的数学意义，门槛值的确定可采用试算的方法，既保证一定的线性化精度，又尽可能多地去除曲线中的线性部分。经过大量刚度曲线的比对和试算，设定门槛值e=5%可取得良好数值效果。

3 典型紧固件刚度曲线

以机翼连接后梁和蒙皮的一颗螺栓为例，螺栓直径为6.35 mm，后梁厚度为6.57 mm，蒙皮厚度为6 mm。材料参数如表1所示。

在ABAQUS中建立图6所示的实体单元模型，计算该紧固件的刚度曲线，如图8所示。

有限元计算所得的优化前刚度曲线如图8中蓝色图例所示，包含50个数据点，橙色图例为经过上述优化计算后的刚度曲线，包括6个点。两条曲线的吻合度非常好，通过优化可以减少刚度曲线数据点的数量，去除原有的“线性段”，但同时不会降低刚度曲线的数值精度。

4 仿真流程

模型的仿真计算流程，拆解为两步。第一步采用通用分析步（Static，General），大量试算表明采用静力分析步，在不打开稳定控制下，计算结果只能计算模型有单元进入塑性左右，计算终点表现为长桁失稳，此后进入第二步计算。第二步采用准静态分析步（Quasi-Static），稳定计算到结构破坏。

对于模型中紧固件的处理采用如图9所示的仿真流程。通过不断对模型添加弹塑性刚度曲线的方式，将模型中首先出现高钉载的钉加入弹塑性刚度曲线，当紧固件进入塑性段，紧固件的承载增量降低，为了满足结构的承载要求，载荷会重新分配到高载荷紧固件周围载荷较低的紧固件上，从而降低紧固件高应力，同时这种仿真更加符合实际情况，可以提高仿真精度。

5 仿真结果

采用上述方法，对翼盒精细模型中相关的232个紧固件的刚度曲线进行简化，在不更改其他非线性设置的基础上，重复了非线性有限元分析，结构开始发生不稳定破坏的状态如图10所示。承载能力计算结果如表2所示，其中试验测得的结构承载能力以F表示。最终数值结果与试验的误差为7%左右，与金属材料性能的保守取值有关。分析达到工程计算的精度要求。

该非线性计算过程中需要反复调用紧固件刚度曲线进行迭代运算。优化前由于刚度曲线定义不合理，且涉及的紧固件数量较多，导致结果不容易收敛; 当刚度曲线优化后，数值迭代的速度显著提高，原有的收敛性问题得到解决。

仿真结果显示在0.90F的加载条件下，后梁与蒙皮之间的紧固件最先出现高钉载，此时的长桁已经出现了损伤，长桁的承载增量降低，前后梁将承载大部分的载荷增量，蒙皮在载荷作用下向上拉，紧固件作为蒙皮和后梁之间的传力零件，当载荷加载到0.91F时，原来高钉载的钉出现了断裂。当首个紧固件发生断裂后，结构的载荷会传递到周围的紧固件上，载荷会分配到周围的紧固件上，接着周围的钉会出现拉链式的破坏，如图11所示，其破坏形式与图12所示试验件的破坏形式基本一致。

6 结论

在实际复杂飞机结构非线性仿真中，非线性刚度曲线是影响载荷传递准确性的重要因素之一，而复杂模型通常含有大量紧固件牌号和不同的连接板厚度组合。本文在三维有限元模型计算紧固件刚度的基础上，通过控制中点误差值，能够实现对紧固件刚度曲线的自动优化，成功解决了非线性刚度曲线所引起的数值收敛性问题，为实现复杂模型下大量考虑非线性紧固件提供了基础。

参考文献