在结构装配中，螺栓是最常见的连接方式之一，在大型结构的冲击数值分析中，由于螺栓和螺栓孔在尺度上通常与被连接结构相差较大，全面精细建模复杂，计算工作量大，所以实际工程计算中需要对连接结构进行简化^[1-3]。目前，常采用弹簧或者梁单元简化螺栓连接，这种简化一方面不能准确地描述螺栓的复杂连接行为，另一方面这些弹簧或者梁单元的刚度也无法准确给出，这种简化的处理方法还不能计及螺栓预紧力与摩擦力的影响，而研究资料显示预紧力和摩擦力与连接刚度密切相关^[4]。因此，为了在大型结构的数值分析中既经济又准确描述螺栓行为，需要建立连接结构的连接本构模型。

虽然螺栓连接无处不在，但以往的文献对螺栓连接关系，即载荷位移曲线变化的物理机理研究不足^[5]。特别是螺栓在剪切载荷作用下的接触和连接是复杂的多尺度非线性问题，不同的外力作用下螺栓的载荷传递形式有很大差异。小载荷下振动问题的载荷传递通常是在螺栓预紧力作用下，被连接的两个部件之间的微观界面摩擦传递，在这个过程中界面产生微滑移，数值模拟中常用弹簧阻尼模型^[6-7]描述连接关系。桑迪国家实验室基于微滑移/微摩擦理论得到简洁的四参数的Iwan模型^[8]，主要关注部件相互间界面的能量耗散；H.Tian等^[9]、黄开放等^[10]通过假设微球体接触表面，计算得到了界面的刚度，但模型参数的给出较为困难。随着载荷的增加，部件相互间的界面发生了宏观滑移，但两个部件分别与螺母之间的界面保持黏结，载荷主要通过螺母与部件之间的界面和螺杆的弯曲变形传递，这方面研究较匮乏。随着载荷的继续增加，载荷主要通过螺杆的剪切变形传递直至材料失效破坏，这个过程研究人员多关注大载荷下螺栓连接的失效模式以及极限载荷，这方面研究资料很多^{[5, 11]}，而对连接刚度较少提及。综合发现螺栓连接在部件相互间的界面发生宏观滑移后到螺杆与螺孔壁完全接触前的载荷传递与变形状态基本没有文献进行讨论。

对工程实际中火箭导弹等结构，在受到冲击载荷的作用下，需要分析冲击产生的应力波对其中精密部件的影响^[12]。例如，火箭飞行达到一定高度后的整流罩解锁分离，分离动载荷经过含连接的结构传递，有时会使得精密部件如剪断锁紧弹射筒和两半罩的锁紧销提前断裂，导致弹射筒无法正常起爆^[13]。还有研究发现应力波的传播特性与螺栓的连接关系及其连接刚度密切相关，是否考虑螺栓的连接刚度，会使应力波传播问题中计算和实验在结构相同监测点的最大加速度相差25%以上^{[1, 14]}。由于螺栓连接的非线性，不同载荷阶段对应不同的载荷传递形式，相应的连接刚度也会发生变化，所以在数值分析中需要给出螺栓连接从低载荷到较高载荷过程中的本构关系。

对于上述问题，本文中采用实验和数值模拟相结合的分析手段，研究单搭接螺栓连接结构在受剪过程中的变形机理，在此基础上得到本构模型的基本形式，对连接刚度进行详细的理论推导，给出了本构模型各阶段的参数确定方法，最后针对不同结构与材料参数对本文建立的单搭接螺栓剪切本构模型及其参数进行了数值验证。

1.   螺栓连接剪切变形机理

首先对螺栓连接结构进行了实验研究。实验件设计如图 1所示，部件为45钢，在连接区的厚度为6 mm，表面粗糙度为1.6 μm，用8.8级M4的螺栓连接，对螺栓施加3.0 N·m的预紧力矩。实验用电子万能试验机夹持住实验件两端夹持区并拉伸，由试验机的力传感器得到载荷，由引伸计测量得到两个部件的相对位移。同时建立与实验相同边界与载荷条件的精细有限元模型，模型采用六面体实体单元，上下两个部件之间的界面设置0.2的摩擦因数，上下部件分别与螺母之间的界面设置0.1的摩擦因数，对螺栓施加通过《机械设计手册》查得3.0 N·m预紧力矩下约为4 kN预紧力，在模型一端夹持区约束，另外一端夹持区施加位移载荷，通过有限元结果分析界面相互作用情况。

图  1  单搭接剪切试样

Figure  1.  Single lap shear specimen

下载: 全尺寸图片 幻灯片

实验载荷位移曲线如图 2所示，随着载荷增加，曲线斜率发生几次明显的变化。发生变化的原因难以直接通过实验分析，为此我们利用有限元仿真计算研究曲线斜率变化的原因。通过螺栓连接安装和受力过程的有限元仿真计算得到了载荷位移曲线，如图 2所示，与实验曲线吻合较好。根据有限元计算的载荷位移曲线斜率可将曲线分为4个阶段，每个阶段的界面接触状态分别如图 3(a)~(d)所示。每阶段的分图中，左边表示部件间的界面，右边表示螺母与上部件间的界面，螺母与下部件的界面与其对称，红色表示界面没有发生滑移，绿色表示界面发生了滑移。图 3(a)所示为第1阶段：界面黏结阶段，该阶段位移几乎不变，连接刚度较大，部件相互间的界面和螺母与部件之间的界面都没有发生滑移，载荷主要通过两个部件相互间的界面传递。图 3(b)所示为第2阶段：部分滑移阶段，该阶段内刚度逐渐降低，部件之间的界面全部发生滑移，螺母与部件之间的界面仍有很大部分黏结在一起，这个阶段的载荷主要通过螺母与部件之间界面摩擦和螺杆的弯曲刚度传递。图 3(c)所示为第3阶段：整体滑移阶段，该阶段内载荷保持不变，刚度几乎为零，部件相互间的界面和螺母与部件之间的界面全部发生了滑移。图 3(d)所示为第4阶段：螺杆剪切阶段，该阶段刚度突然增大，螺杆和部件的孔壁完全接触，载荷主要通过螺杆的剪切刚度传递。对于火箭导弹等结构爆炸分离时产生的冲击波幅值由于衰减较快通常在被考察部件附近不会导致螺栓变形进入图 2所示的第4阶段^[13]，所以本文中针对剪切螺栓连接载荷位移曲线的前三阶段开展研究工作。

图  2  载荷位移曲线

Figure  2.  Load-displacement curves

下载: 全尺寸图片 幻灯片

图  3  接触面状态

Figure  3.  Status of interface

下载: 全尺寸图片 幻灯片

2.   本构模型及其参数确定

基于上节对螺栓剪切载荷位移曲线每个阶段的受力分析，初步确定单搭接螺栓剪切本构模型的表达式见式(1)，基本形式见图 4，其中δ为部件之间在螺孔中心处的相对位移，F为螺栓传递的钉载，μ₁为两个部件之间的摩擦因数，μ₂为部件与螺母之间的摩擦因数，P_pre为螺栓的预紧力，k₀为界面黏结阶段的刚度，k₁为部分滑移阶段的刚度，k₂为整体滑移阶段刚度。

图  4  本构模型

Figure  4.  Constitutive model

下载: 全尺寸图片 幻灯片

可以看出，上述本构模型计及摩擦因数和预紧力，它们影响在多大载荷下螺栓进入部分滑移阶段和整体滑移阶段。另外，确定该本构模型还需各阶段的刚度，其中整体滑移阶段的刚度k₂为0，k₀和k₁采用下列方法计算。

2.1   k₀的确定

k₀为界面黏结阶段的连接刚度，主要为部件之间的界面切向刚度。根据分形模型与Hertz接触理论得出在小载荷下界面的切向刚度为^[15]：

式中：${\bar G}$为界面等效剪切模量，与两接触材料的剪切模量和泊松比相关：

D为分形维数，它与界面加工方法以及粗糙度R_a相关^[16]，对于磨削表面：

h为界面的切向载荷与法向压力比，μ为界面摩擦因数，ψ为微接触面积分布的域扩展因子。a_c为微接触临界接触面积，意为微球体发生屈服时的接触面积：

其中Г为粗糙表面的幅度系数，D=1.5时，其与表面粗糙度R_a的关系为^[17]R_a=0.811Г^1/2，κ为较软材料的硬度Η与屈服强度σ_y之间的系数，γ为较软材料的屈服强度σ_y与材料当量弹性模量E之间的比例系数。a_m为微接触最大接触面积，与界面的法向压力相关^[18]：

式(2)~(6)显示通过两接触材料的弹性模量、剪切模量、屈服强度、表面粗糙度以及界面的摩擦因数可求得k₀的值，在随后的第3节进行算例验证。

2.2   k₁的确定

k₁为部分滑移阶段的连接刚度，即k₁=ΔF/Δδ，这个阶段增加的载荷ΔF通过部件与螺母间的界面作用到螺栓上，受力与变形分析见图 5。将螺杆作为短粗梁，梁两端的力偶(ΔF，-ΔF)会让螺杆产生弯曲变形和剪切变形，从而分别使A、B两点在水平方向增加位移Δδ_b和Δδ_s，对应的刚度为k_b和k_s。另外，这对力偶会在螺杆根部产生附加力矩ΔFL/2，这个力矩会使与螺杆相接的螺母部件产生Δθ_ebn的变形。这个附加力矩通过螺母传递到部件上，又会使部件产生的Δθ_ebc变形，这两个变形导致螺杆发生Δθ的偏转，进而导致A、B两点在水平方向出现由螺母变形导致的位移Δδ_ebn和由部件变形导致的位移Δδ_ebc，因此，k₁由螺杆弯曲刚度k_b、螺杆剪切刚度k_s、部件对螺母的支撑刚度k_ebc以及螺母对螺杆的支撑刚度k_ebn组成，其表达式为：

图  5  螺栓变形示意图

Figure  5.  Diagram of bolt deformation

下载: 全尺寸图片 幻灯片

(1) k_b代表螺杆弯曲刚度，通过经典梁弯曲变形理论得到：

(2) k_s代表螺杆的剪切刚度，由于螺杆一般是短粗梁，所以横向剪切变形不能被忽略，剪切刚度为：

式(8)~(9)中:E和G分别为螺栓的杨氏模量和剪切模量，L为螺杆的长度，ζ为Timoshenko剪切修正系数，对于圆柱，其值为6(1+ν)/(7+6ν)，ν为材料的泊松比。

(3) k_ebc代表部件对螺母的支撑刚度。部件的支撑刚度一直是螺栓连接计算中的难题，螺栓受拉情况下的支撑刚度是一个轴对称问题，S.A.Nassar等^[19]、杨国庆等^[20]等通过部件的应力假设以及轴向位移积分得到部件轴向支撑刚度。而螺栓剪切是非轴对称问题，难以简化为二维问题直接求位移。在线弹性假设下，这里通过应变能法求部件的剪切支撑刚度。

假设螺栓为刚体，由能量平衡可知，外力做功等于部件应变能的增加量，k_ebc有如下表达形式：

所以，只要得到在部分滑移阶段由ΔF导致的应变能U，就可以得到部件对螺母的支撑刚度。应变能可由应力在整个变形体内积分求得：

为了得到合理的应力假设，对部件的应力分布进行有限元数值分析。在图 6中以螺栓轴线为中心线建立圆柱坐标，图中z向应力σ_z是支撑螺母的主要应力，分析中忽略其他应力分量的影响。

图  6  连接结构的应力假设

Figure  6.  Hypothesis of stress for connection structure

下载: 全尺寸图片 幻灯片

在预紧力作用下，假设σ_z分布为锥体形状，如图 6所示，其中虚线为锥体的母线，其半顶角为α，α一般取30°~45°^[20-21]，R_N为螺母与部件接触区半径，R_B为螺杆半径，R_H为螺孔半径，H为螺母厚度。

假设σ_z沿着环向φ的分布为：

沿着径向r分布：由于最大接触应力在螺母与部件接触的边缘即r=R_N处，所以在r<R_N和r>R_N区域内有不同的分布表达式，具体为：

式中：R_z为在z处锥体的半径:

沿着轴向z的分布为:

假设分布(Hyp.)和有限元结果(FEM)对比如图 7所示，表明假设分布形式能够比较好地描述实际应力分布。

图  7  部件假设应力分布与有限元应力分布比较

Figure  7.  Stress comparison of hypothesis with FEM on component

下载: 全尺寸图片 幻灯片

因此，部件的三维应力分布可表示为：

式中：$a = R_{\rm{N}}^2 - R_{\rm{H}}^2,b = {R_{\rm{N}}}{\rm{tan}}\alpha {\rm{,}}{c}{\rm{ = }}\frac{1}{3}{\rm{ta}}{{\rm{n}}^2}\alpha$。

利用z=0处的力矩平衡关系：

同时将式(17)右端积分为代数表达式，求解可获得三维应力分布中F_z的表达式：

利用三维应力分布计算部件的弹性应变能，这个过程中对为r<R_N和r>R_N区域分别积分如下：

式中:${I_0} = \int_0^L {\frac{1}{{{{\left( {a + bz + c{z^2}} \right)}^2}}}{\rm{d}}z} , \;{I_1} = \int_0^L {\frac{z}{{{{\left( {a + bz + c{z^2}} \right)}^2}}}{\rm{d}}z} , \;{I_2} = \int_0^L {\frac{{{z^2}}}{{{{\left( {a + bz + c{z^2}} \right)}^2}}}{\rm{d}}z}$。

将(19)代入(10)并消去ΔF，最终得到部件对螺母的支撑刚度:

(4) k_ebn代表螺母对螺杆的支撑刚度。

与求k_ebc方法相同，首先根据有限元分析结果得到螺母的假设应力分布为：

式中:n取5。

在每个坐标系方向假设分布与有限元对比如图 8所示，在φ向和z向假设分布形式能够较好描述实际应力分布，r向的应力分布在r较小时有一定误差，但由于应变能是应力的平方在体积上的积分，假设应力分布与真实应力分布误差主要处于小半径、低应力区域，对于应变能的结果影响有限。

图  8  螺母假设应力分布与有限元应力分布比较

Figure  8.  Stress comparison of hypothesis with FEM on nut

下载: 全尺寸图片 幻灯片

利用三维应力分布计算得部件的弹性应变能，这个过程中对为r<R_B和r>R_B区域分别积分如下：

进而得到螺母对螺杆的支撑刚度：

至此，通过部件与螺栓的几何参数、弹性模量、硬度、屈服强度以及界面粗糙度采用式(2)、(7)分别得到界面黏结阶段刚度k₀和部分滑移阶段k₁，结合预紧力和摩擦因数以及已知为0的整体滑移阶段刚度k₂就能完整给出本文中提出的螺栓剪切模型。

3.   验证

图 2显示，精细有限元模型能够准确模拟连接结构的力学行为，所以采用精细有限元模型对本文中建立的单搭接螺栓连接剪切本构模型进行验证。

模型采用8.8级M6螺栓连接，螺栓的弹性模量为200 GPa；两个被连接部件为45钢，弹性模量为210 GPa，屈服强度为360 MPa，硬度为200 MPa，厚度均为12 mm。对螺栓施加10 kN的预紧载荷，该载荷为螺栓屈服载荷的55%，部件之间以及部件与螺母之间的摩擦因数为0.1。部件的表面粗糙度R_a为1.6，取扩展因子ψ=1.2，由式(2)计算得到在低切向载荷下，k₀=3.54×10⁷ N/mm。由式(7)计算得到k₁=7.86×10³ N/mm。

有限元计算得到的载荷位移曲线与预测载荷位移曲线如图 9所示，位移为上下螺母与部件接触区域中心点的位移差。

图  9  预测结果与有限元结果对比

Figure  9.  Comparison of prediction with FEM

下载: 全尺寸图片 幻灯片

图 9显示，该模型预测结果与有限元结果主要在部分滑移阶段有偏差，偏差的大小可用这个阶段的刚度k₁的误差描述。下面分别分析螺栓和部件不同选材、不同螺栓直径、不同部件厚度对k₁的误差影响。

表 1中分别对螺栓和部件选用工程中常用的材料铝合金和钢进行验证，两个被连接部件等厚，均为15 mm，螺栓直径为10 mm，钢的弹性模量为210 GPa，铝的弹性模量为70 GPa。分析结果显示，材料的弹性模量对k₁的误差影响不大，都在9%左右。

表  1  不同结构选材k₁验证

Table  1.  Verification of k₁ for different material

螺栓	部件	k1/(N·mm-1)		误差/%
		预测值	有限元值
钢	铝	1.83×104	1.68×104	8.8
铝	钢	9.69×103	8.88×103	9.1
铝	铝	8.44×103	7.70×103	9.5

下载: 导出CSV 

| 显示表格

表 2对不同的螺栓直径进行验证，两个被连接部件的厚度均为15 mm，部件和螺栓的材料都为钢材，其模量为210 GPa，分析结果显示，螺栓直径影响k₁的误差，直径增大会导致误差变大。螺栓直径小于10 mm时，误差小于10%，直径为16 mm时，误差增大至28%。

表  2  不同螺栓直径k₁验证

Table  2.  Verification of k₁ for different bolt diameters

螺栓直径/mm	长细比	k1/(N·mm-1)		误差/%
		预测值	有限元值
6	5.0	4.29×103	4.13×103	3.8
10	3.0	2.53×104	2.30×104	9.9
16	1.9	1.17×105	9.14×104	28.3

下载: 导出CSV 

| 显示表格

表 3对不同部件厚度的结构进行了分析，表中螺栓直径为10 mm，部件和螺栓的材料都为钢材，其模量为210 GPa，分析结果显示部件越厚，k₁误差越小。当被连接两部件总厚度大于40 mm时，误差小于10%，小于20 mm时，误差接近15%。

表  3  不同部件厚度k₁验证

Table  3.  Verification of k₁ for different thicknesses

部件厚度/mm	长细比	k1/(N·mm-1)		误差/%
		预测值	有限元值
20	2.0	6.57×104	5.76×104	14.0
40	4.0	1.23×104	1.03×104	8.3
70	7.0	2.76×103	2.64×103	4.7

下载: 导出CSV 

| 显示表格

综合表 1~3的验证结果显示，k₁的预测值都大于有限元值，其中螺栓和部件的选材对误差影响较小，螺栓和部件的尺寸对误差有影响。表 2~3显示，误差与部件厚度/螺栓直径(长细比)显著相关，长细比大于2.0的情况下，给出的k₁计算方法能够得到较为准确的结果(误差小于14%)。由于推导k_ebn和k_cbn过程中，仅主要方向的应力σ_z，计算中的应变能小于实际应变能，导致k_ebn和k_cbn的预测结果比有限元支撑刚度偏大，又由于k_ebn和k_cbn是k₁的分量，并且长细比越小，k_ebn和k_cbn的影响越大，所以k₁随着长细比减小，与有限元刚度偏差明显。

4.   结论

(1) 通过实验和计算得到单搭接螺栓剪切结构的载荷位移曲线，根据曲线斜率对不同阶段进行了物理解释：连接载荷首先通过部件相互间的界面传递；当部件相互间的界面发生滑移后，载荷通过螺母与部件之间的界面以及螺栓杆的弯曲传递；当螺母和部件之间的界面发生滑移后，连接结构都处于整体滑移阶段，载荷不再继续增加。

(2) 提出了计及摩擦因数、预紧力、螺栓与部件尺寸、材料等参数的单搭接连接结构的剪切本构模型。给出了模型中参数的获取方法。该模型适用于大型结构在中低载冲击应力传播数值分析中大量螺栓连接简化。

(3) 本构模型针对部分滑移阶段刚度k₁的计算采用应变能法获得了部件对螺母的支撑刚度k_ebc以及螺母对螺杆的支撑刚度k_ebn的影响，这种方法能够解决三维非轴对称变形体对梁的支撑刚度问题。