紧固件连接是飞机结构中广泛存在的一种连接方式，关系着整个结构的安全性和可靠性，虽然构造简单，但受力机理复杂。
在紧固件三维实体有限元建模方面，Kim[1]对单个连接紧固件的结构特征进行了详细建模。三维实体建模方法的好处是能将接紧固件的所有细节进行建模，但缺点是建模过程繁琐且需要大量计算时间。因此，建立具有高效性、准确性的非三维实体有限元模型就显得尤为重要。
在非三维实体有限元建模方面，目前绝大多数研究是将螺栓连接简化为梁单元、弹簧单元或者是有限元连接器单元[2-6]，这样简化处理在保证计算精度的前提下，能快速建立有限元模型，节约计算成本。另外，连接件之间的接触也是一个关键因素，徐梓雯等[7]利用 GAP 单元模拟了螺栓连接件之间的弹性接触，建立了更加合理的力学简化模型，但需要借助有限元外部程序进行干预。
紧固件刚度对非三维实体有限元模型的计算精度有极大影响，紧固件刚度直接影响结构的安全性[8]。在紧固件载荷分析以及结构疲劳计算时，赋予紧固件实际的剪切刚度值是最佳的选择，但紧固件的刚度分析较为复杂。目前，人们往往采用半经验公式计算紧固件刚度值。例如Tate & Rosenfeld、Huth-Schwarman、Swift和Grumman等半经验公式。为验证半经验公式的可靠性和准确性，学者们对半经验公式进行了大量评估工作[9-10]，研究结果表明Huth-Schwarman半经验公式在飞机结构上具有较好的效果[11]。但半经验公式毕竟是工程经验得出的一种表达形式，对新材料和新研产品是否适用还有待实践检验，在复杂结构中所得到的刚度值还有待商榷。另外，林海彬等[12]将螺栓的工程算法与有限元仿真的结果做出了对比，发现有限元的计算结果相对工程算法结果更符合螺栓的实际受力状态。
在飞机结构设计过程中，不但要校核紧固件自身强度是否满足设计要求，同时还要给出紧固件的载荷分布。载荷分布的准确性是疲劳寿命预测精度的关键之一，若紧固件载荷分布不合理，将导致局部应力过大，影响结构安全和结构疲劳寿命。采用非三维实体有限元建模方法广泛应用于紧固件载荷分布的计算[13-14]。例如，Naruse等[15]通过对螺栓的力学分析，研究了螺栓的静态失效和疲劳失效。姚卫星[16]曾用杆-板-Spring单元计算了飞机壁板螺栓连接的载荷分布，并利用SSF(Stress Severity Factor Stress Serious Factor)法预测了结构的疲劳寿命。
针对非三维实体有限元建模方法，本文依托ABAQUS软件，采用2种非三维实体紧固件有限元模型，并对紧固件载荷分布进行了研究分析，并通过与试验数据和实际工程案例对比，验证了紧固件有限元模型的有效性和高效性。
1 紧固件有限元模型
紧固件有限元模型有2种，分别为ABAQUS中的Springs/Dashpots连接器以及Bushing连接器，如图1所示。
图1 紧固件有限元模型示意图
Springs/Dashpots连接器需要输入载荷方向的剪切刚度。当考虑预紧力时需要与Cartesian连接器一起使用，此时需要输入紧固件轴向刚度。紧固件轴向刚度为
(1)
式中： S为紧固件横截面积；L为紧固件长度；Ef为紧固件弹性模量。
Bushing连接器需要输入紧固件3个方向刚度，可以直接考虑预紧力。3个方向的刚度值分别为K1、K2和K3，在有限元模型中刚度方向如图2所示。Bushing连接器可以很好的模拟出Springs/Dashpots连接器与Cartesian连接器一起使用的效果。
图2 紧固件刚度示意图
对于连接器与连接件之间采用较为柔性的Continuum distributing约束来取代常用的MPC刚性单元，主要原因有2个：一是由于紧固件刚度较大，再加上连接件之间的接触，容易在节点处出现局部刚度过大的问题；二是因为在紧固件与连接件之间一般会有垫片，Continuum distributing弹性约束影响的区域更接近真实的连接效果。
2 半经验公式
紧固件刚度在载荷分布、强度校核以及疲劳分析时，都是极其重要的影响因素。Siddabathuni[10]的研究结果表明：在荷载传递和变形方面，紧固件的刚度(柔性)决定了荷载从一个部件传递到另一个部件的方式。假设有无限刚度的紧固件，那么两边的紧固件将承受所有载荷，中间紧固件不承受载荷；假设是无限柔度紧固件，那么载荷将均匀分配在每个紧固件上面。Siddabathuni[10]给出的紧固件载荷分布如图3所示，由此可见紧固件刚度对载荷分布的影响是极其重要的。
图3 刚度影响下的紧固件载荷分布示意图
紧固件刚度K2和K3在有些情况下是没有实测值的。因此，刚度值可借助半经验公式计算。
Tate & Rosenfeld半经验公式[17]为
C=
(2)
式中：Ef为紧固件弹性模量； t1、t2分别为板厚度； E1、E2分别为板弹性模量；C、K分别为紧固件柔度及刚度；D为紧固件直径。
Swift半经验公式[18]为
(3)
Grumman半经验公式[19]为
(4)
Huth-Schwarman半经验公式[11]为
(5)
式中：n、a和b为常数，取值如表1所示。
表1 参数取值
符号描述n连接为单剪接头时n=1,连接为双剪接头时n=2a紧固件为螺栓时a=2/3,紧固件是铆钉时a=2/5b紧固件为金属螺栓时b=3,紧固件为复合螺栓时b=4.2,紧固件为铆钉时b=2.2
从半经验公式可以看出：所有的半经验公式都考虑了板材料弹性模量、板尺寸、紧固件直径和紧固件材料等因素的影响。不同的是Huth-Schwarman半经验公式还考虑了单剪、双剪以及紧固件类型的影响，在飞机结构中应用较为广泛。
3 实例验证
3.1 初始参数
该实例来源于文献[11，13]，边界条件为沿金属板横向施加单位载荷，预紧力为3 600 N，试验和计算模型如图4所示。板材料为2024-T3铝合金，紧固件材料为Ti 6Al4V钛合金，具体参数见表2。
图4 连接件模型示意图
表2 紧固件以及板参数
符号说明数值D/mm紧固件直径5t/mm板厚度5.1E/MPa板弹性模量72 000Ef/MPa紧固件弹性模量110 300K1/(N·mm-1)紧固件轴向刚度212 327K2、K3/(N·mm-1)紧固件剪切刚度45 535
3.2 载荷分布数值计算
Gunbring[13]利用三维实体建模，得到的紧固件载荷分布如图5所示。2种连接器紧固件模型的计算结果见图6和图7。可以看出：3种模型求得的载荷分布大体相同，两边紧固件承受的载荷略大于中间紧固件承受的载荷，与Siddabathuni[10]给出的紧固件载荷分布(图3)具有相同的趋势。
图5 采用三维实体模型计算的载荷分布示意图
图6 采用Springs/Dashpots连接器计算的载荷分布示意图
图7 采用Bushing连接器计算的载荷分布示意图
3.3 对比分析
将Huth给出的实验数据[11]、Gunbring利用三维实体建模得到的计算结果[13]与本文方法的计算结果汇总于表3。可以发现，三维实体单元、Springs/Dashpots连接器以及Bushing连接器与试验数据均具有较好的吻合度，说明连接器紧固件模型具有满意的计算精度。然而在计算效率方面，连接器紧固件模型耗时大约1 min，而三维实体模型大约耗时半个小时，连接器紧固件模型的计算效率明显高于三维实体模型。
表3 计算结果
孔号实验[11]三维实体模型[13]本文(Springs/Dashpots)本文(Bushing)135%34.4%35.39%34.85%230%31.2%29.31%30.14%335%34.4%35.30%35.01%时间/h-0.50.020.01
4 工程应用
应力严重系数法是针对飞机结构连接件而发展起来的一种疲劳寿命估算方法。应力严重系数可通过钉传载荷和旁路载荷求得，如果应力严重系数的求解精度不足，将直接影响疲劳危险部位的判断和疲劳寿命的预测精度。
应力严重系数法属于名义应力法的范畴，大体步骤如下：通过飞机结构的有限元细节分析，得到钉孔处的钉传载荷和旁路载荷；求出孔边的SSF值，以SSF值和名义应力 Snom的乘积最大为依据确定疲劳危险孔；将SSF看作理论应力集中系数KT； 利用缺口件的S-N曲线求出当前应力集中系数和应力水平下的S-N曲线；应用疲劳损伤累积理论(如Miner线性累积损失理论)计算疲劳寿命。
某飞机机翼下表面典型螺栓连接试验件如图8所示，连接件材料为LY12-CZ铝合金，螺栓剪切刚度为9 921 kg/mm，其他参数和载荷谱详见文献[16]。
图8 螺栓连接接头试验件
有限元细节分析如图9所示。由于对称关系，在有限元分析时取试验件的1/2，即取1～6号孔钉传载荷以及旁路载荷即可。由于试验件只受面内拉压载荷且不考虑预紧力，因此采用Springs/Dashpots连接器。本文方法计算的钉传载荷与旁路载荷结果见图10和图11，与文献[16]的对比见表4。
图9 有限元细节分析示意图
图10 带板有限元分析结果示意图
图11 接头有限元分析结果示意图
表4 有限元细节分析结果
钉孔号带板123接头123文献[16]旁路载荷0.651 60.350 3000.348 40.649 7钉传载荷0.348 40.301 30.350 30.348 40.301 30.350 3SSF3.741 84.283 96.840 55.443 73.638 03.255 1文献[16]旁路载荷0.643 40.3634000.356 60.636 6钉传载荷0.356 60.280 00.363 40.356 60.280 00.363 4SSF3.700 64.15436.840 25.443 33.562 03.305 3
由表4可以看出：本文方法所得旁路载荷与钉传载荷与文献[16]结果略有不同。对于旁路载荷，1号钉、6号钉降低，2号钉、5号钉提高，而3号钉、4号钉不变，均为0。对于钉传载荷，1号钉、3号钉、4号钉、6号钉增加，2号钉、5号钉降低。
通过将SSF值和名义应力Snom乘积可知，1号孔的 Snom×SSF=2.525 3×3.700 6=9.345 1最大，即疲劳危险部位为1号孔，与文献[16]给出的疲劳危险位置相同。然后将SSF=3.700 6看作理论应力集中系数KT，查询相应的S-N曲线，利用Miner线性累积损伤理论计算的疲劳寿命见表5。由表5可以看出：本文方法预测的疲劳寿命更接近试验值，相比文献[16]采用杆-板-Spring单元的计算精度提高了8.08%。
为了分析半经验公式给出的紧固件刚度值对疲劳寿命的影响，本文在紧固件有限元模型中，螺栓剪切刚度分别取文献[16]给出的真实剪切刚度(9 921 kg/mm)和由Huth-Schwarman半经验公式给出的剪切刚度(12 739 kg/mm)进行计算，计算结果见表5。由表5可以看出：采用Huth-Schwarman半经验公式剪切刚度的计算结果略次于真实剪切刚度的计算结果，说明紧固件剪切刚度对紧固件载荷分布和疲劳寿命具有一定影响。因此，在有条件下需要开展紧固件剪切刚度的试验。
表5 疲劳寿命计算结果
预测飞行次数实测飞行次数[16]误差/%文献[16]6 9138 42117.91本文方法(真实刚度)7 5938 4219.83本文方法(经验刚度)7 3838 42112.33
5 结论
文中给出的紧固件有限元模型适用于计算紧固件结构载荷分布，其计算结果与实验值吻合较好，相比三维实体模型，建模简单，计算效率高。
基于本文给出的紧固件有限元模型，钉传载荷和旁路载荷的计算精度有所提高，提高了利用SSF法估算连接件结构疲劳寿命的预测精度。紧固件剪切刚度对载荷分布具有一定的影响，影响疲劳寿命的预测精度。应尽量获得紧固件真实剪切刚度，如无法获取真实剪切刚度，可借助半经验公式获取紧固件剪切刚度进行初步估算。
[1] KIM J，YOON J C，KANG B S.Finite element analysis and modeling of structure with bolted joints[J].Applied Mathematical Modelling，2007，31(05):895-911.
[2] LIU L，CHEN K.Modeling method for deformation analysis in thick laminate mechanical joints[J].Journal of Aircraft,2012，49(06):1974-1981.
[3] GRAY P J，MCCARTHY C T.A global bolted joint model for finite element analysis of load distributions in multi-bolt composite joints[J].Composites Part B，2010，41(04):317-325.
[4] BELARDIA V G，FANELLIB P，VIVIOA F.FE analysis of single-bolt composite bolted joint by means of a simplified modeling technique[J].Procedia Structural Integrity，2019，24(11):888-897.
[5] HONG H，LEE H，JEONG N，et al.A study on an equivalent model of the threaded fasteners in complex structures through tightening and loosening analysis[J].Journal of Mechanical Science & Technology，2020，34(06):1195-1205.
[6] 张晓頔，郝志明，李一堃，等.含螺栓连接结构的非线性静力学拓扑优化[J].兵器装备工程学报，2018，39(08):74-78.
[7] 徐梓雯，那景新，张志远，等.螺栓连接有限元模型的弹性接触研究[J].中国机械工程，2012，23(15):1830-1832.
[8] HAIDAR N，OBEED S，JAWAD M.Mathematical representation of bolted-joint stiffness:A new suggested model[J].Journal of Mechanical Science and Technology，2011，25(11):2827-2834.
[9] SÖDERBERG J.A finite element method for calculating load distributions in bolted joint assemblies[D].Sweden:Link pings University，2012.
[10] SIDDABATHUNI V G U，MEKAM S，SASTRY S S S.Parametric studies on the effect of four types of fastener modeling in channel type tension fitting[J].International Journal of Recent Advances in Mechanical Engineering，2015，41(01):15-27.
[11] HUTH H.Influence on fastener flexibility on the prediction of load transfer and fatigue life for multi-row joints[M].Fatigue in mechanically fastened composite and metallic joints ASTM STP927 Philadelphia:ASTM，1986:221-250.
[12] 林海彬，陈雷，孙英超，等.吊挂横梁螺栓拉力计算结果差异性分析[J].兵器装备工程学报，2019，40(10):216-220.
[13] GUNBRING F.Prediction and modelling of fastener flexibility using FE[D].Sweden:Linköpings University，2008.
[14] KONKONG N，PHUVORAVAN K.An analytical method for determining the load distribution of single-column multibolt connection[J].Advances in Civil Engineering，2017(01):1-19.
[15] NAEUSE T，KAWASAKI T，HATTORI T.Simple modelling and strength evaluation methods for bolt joints using shell elements and beam elements[J].Journal of Computational Science and Technology，2009，3(01):34-45.
[16] 姚卫星.结构疲劳寿命分析[M].北京:科学出版社，2019.
[17] TATE M B，ROSENFELD S J.Preliminary investigation of the loads carried by individual bolts in bolted joints[R].NACA-TN-1051，1946.
[18] SWIFT T.Fracture analysis of stiffened structure.In:Damage tolerance of metallic structures:Analysis methods and application.ASTM STP842[M].Philadelphia:ASTM，1984.
[19] JARFALL L.Shear Loaded Fastener Installations[J].International Journal of Vehicle Design，1986，7(3/4):337-380.