0 引言

钢结构节点的连接可采用多种方式，其中高强度螺栓连接施工简单^[1-2]、整体性刚度较好、耐疲劳，是钢桥安装的主要连接构件。高强度螺栓连接可分为摩擦型连接和承压型连接两种，其中，摩擦型高强度螺栓连接应用更为广泛^[3]。摩擦型连接是依靠连接件之间的摩阻力传递内力，并以外荷载引起的剪力不超过摩阻力作为设计准则。对于摩擦型高强度螺栓连接节点的承载能力研究可采用试验法和有限元仿真分析法。国内外对高强螺栓受力状态已进行了一些试验研究^[4-6]，但采用的试验方法时间花费较长，经济耗费较大。随着有限元方法的发展和计算机技术的提高，有限元仿真方法成为一种有效的研究手段。

目前，各国有关钢结构的设计规范中均给出了摩擦型高强螺栓承载力设计值计算公式，而未考虑其滑移变形，这是由于结构所承受的荷载未超过其设计承载力时，螺栓的滑移变形并不显著。但是在结构的极限承载力分析中，结构往往承受远大于其设计承载力的荷载，大量工程事故和结构试验表明，此时螺栓的滑移变形现象在钢结构的破坏中是非常普遍的^[7]。钢桁架桥梁的节点一般采用摩擦型高强螺栓实现抗剪连接，节点螺栓的滑移会对节点的转动刚度产生影响^[8]，进而影响整个桥梁的正常使用。因此准确地分析钢结构全过程的受力和变形性能，对摩擦型高强螺栓抗剪连接的滑移量及滑移过程进行深入研究是十分必要的，主桁下弦节点高强螺栓连接如图 1所示。

1 摩擦型高强螺栓抗剪滑移有限元模拟

摩擦型高强度螺栓抗剪连接的加载过程具有高度非线性特征，在仿真分析中需要多次反复迭代计算，计算时间较长且计算结果往往不收敛。采用大型通用有限元软件ANSYS模拟摩擦型高强螺栓一般有2种思路，一种是采用实体单元模拟钢板^[9]，接触面采用接触单元和目标单元形成接触对，高强螺栓预拉力通过预紧单元施加。由于螺栓本身尺寸较小，想要计算出理想的结果，就必须把单元划分得足够细小，这就对计算机硬件配置提出了很高的要求。虽然现阶段计算机性能得到了很大发展，但是对于比较复杂的螺栓连接节点，需要多次迭代计算，想要得到理想的计算结果仍比较困难。另一种方法是采用壳单元模拟钢板，接触面间的摩擦行为利用弹簧单元解决。为了模拟摩擦型高强螺栓连接的力学性能，已有许多学者建立了简化的数值模型。李海岗^[10]采用2个弹簧单元模拟1个螺栓的摩擦行为，可以精确捕捉螺栓开始滑移时的荷载，但是既不能模拟螺栓的实际滑移量，也不能很好反映芯板和盖板的应力规律。赵中伟等^[11]提出了一种简化的螺栓滑移数值模型，并且给出了计算摩擦单元实常数的计算公式，通过对简单螺栓连接实例的滞回分析，表明该方法可以精确模拟高强螺栓的滑移，结构的收敛性得到本质改善。

1.1 简化的螺栓滑移数值计算方法

运用ANSYS有限元软件模拟摩擦型高强螺栓拼接接头，钢板均采用SHELL181单元模拟，螺栓采用BEAM188单元模拟，采用COMBIN39弹簧单元模拟螺栓与盖板、盖板和芯板之间的摩擦行为，同时，COMBIN39单元也被用来模拟螺杆与螺栓孔壁之间的接触挤压行为。螺栓和COMBIN39弹簧单元之间采用MPC184刚性梁连接。

COMBIN39称为非线性弹簧单元，是一种具有非线性的广义力-变形曲线(简称F-D曲线)的单向单元，适用于各种分析^[12]。该单元具有轴向或扭转功能，轴向弹簧为单轴拉压行为，每个节点有3个自由度，即坐标系x，y，z方向的平动位移。COMBIN39单元参数KEYPOT(1)为卸载路径的选择，设置为0时表示卸载路径与加载路径相同，设置为1时表示卸载路径与加载路径的原点段平行。模拟接触挤压行为时，将参数KEYPOT(1)设置为0，模拟摩擦行为时，将参数KEYPOT(1)设置为1，有摩擦行为的系统为非保守系统，即系统有能量的损失。将COMBIN39单元的轴向设置为沿螺孔的径向，由螺孔圆心向外(见图 2)。参数KEYPOT(3)的设置与建立COMBIN39单元时所使用的总体坐标系相关，为了便于编制APDL宏文件建立单个螺栓，建议将KEYPOT(3)设置为0。APDL宏文件是集成一系列的命令后缀为MAC的可调用文件，用来完成大批量的重复问题^[13]。将建立单个螺栓的命令流编制成APDL宏文件BOLTZ. MAC，采用循环指定要建立螺栓的位置，反复调用BOLTZ.MAC可以方便实现全部螺栓的建立，大幅缩短建模时间。

COMBIN39单元模拟摩擦行为时，F-D曲线如图 3(a)所示。COMBIN39单元模拟接触挤压行为时，F-D曲线如图 3(b)所示。为了保证COMBIN39单元在计算时的收敛性，图 3中θ不能等于90°，为了保证计算精度，尽量使θ接近90°。将接触单元在受压时的承载力设置为10⁶ N已足够大，可以防止螺栓杆与螺栓孔钢板之间的侵入。

非线性迭代计算时，随着施加外力的逐渐增大，某个螺栓位置摩擦单元首先达到最大静摩擦力F_max。外力继续增大直到所有螺栓位置摩擦单元均达到最大静摩擦力值，螺栓群整体滑移，滑移最大值为(d₁-d₂)/2，d₁和d₂分别为螺栓孔直径和螺杆直径。外力继续增大，螺杆和螺孔接触，按照所给定的本构模型发生非线性变化。摩擦单元所能承受最大内力F_max计算公式为^[11]：

(1)

式中，μ为滑动摩擦系数；m为摩擦面数量；P为螺栓预紧力；n为沿螺栓孔壁的节点个数。

COMBIN39单元模拟接触行为时，实常数定义为：R，4，- (d₁-d₂)/2，-10⁶，- (d₁-d₂)/2+0.001，-10²，0，0。COMBIN39单元模拟摩擦行为时，实常数定义为：R，5，0，0，(d₁-d₂)/2/100，F_max/10，(d₁-d₂)/2，F_max，RMORE，(d₁-d₂)/2×10，F_max。实常数定义中，模拟接触行为时单元实常数号为4，模拟摩擦行为时单元实常数号为5。

1.2 单行螺栓试验验证

为验证本研究数值计算方法的准确性，首先建立单行螺栓有限元模型，将有限元计算结果与文献[14]中的单行两列螺栓实例试验结果进行对比。文献中单行螺栓实例各项参数见表 1，实例尺寸如图 4所示。

简化的单行螺栓滑移有限元模型如图 5所示。分析时，一端固定，另一端施加拉力。

单行螺栓有限元结果和文献[14]试验结果对比如图 6所示。单行螺栓实例所能承受的最大摩擦力F_sum=μ×m×P×n=155 kN。单行螺栓有限元模型发生显著滑移时的拉力为154 kN。由图 6可见，试验荷载-位移曲线与有限元荷载-位移曲线吻合较好，弹塑性阶段试验对应的荷载低于有限元结果，这是由于钢材本构模型与试验存在误差。

1.3 多行螺栓数值模型的建立

本研究多行螺栓连接采用文献[15]中试件Ⅰ建立简化螺栓抗剪滑移模型，模型中拼接板采用Q235，摩擦系数取0.45，采用M16螺栓，栓孔直径为17 mm，预紧力取文献中数值为80 kN。试件尺寸如图 7所示。

采用ANSYS经典界面命令流方式建模，命令流分为主程序部分和可重复调用APDL宏文件BOLTZ.MAC两部分。主程序部分定义参数，建立拼接板有限元模型，其中钢板弹性模量取E_s=2.06×10⁵ MPa，假设钢材为各向同性材料，泊松比为0.3，材料本构选用双线性等向强化模型，采用双直线描述材料的应力-应变关系。APDL宏文件BOLTZ.MAC中建立栓杆、刚性梁，根据主程序中孔壁划分的节点个数建立摩擦单元和接触单元，旋转摩擦单元和接触单元的节点坐标系方向，从螺栓中心沿径向向外，施加每个螺栓边界条件。分析中将被连接钢板的一端固定，另一端施加拉力。根据对称性，本研究只研究受拉力一侧螺栓传力规律，受拉力一侧螺栓编号见图 7。图 8为本研究建立的有限元模型。

2 螺栓连接的传力性能

多行螺栓拼接接头在螺栓滑移之前处于摩擦阶段，第一个螺栓松动到所有螺栓达到最大滑移量为滑移阶段，当所有螺栓达到最大滑移量即所有栓杆和孔壁接触后，模型进入弹塑性变化阶段。由于螺栓屈服强度较大，当抗剪螺栓连接达到极限承载力时，会出现螺栓孔局部挤压破坏，与受力较大的螺栓接触挤压的螺栓孔区域首先进入塑性状态。按照文献[16]，一个高强螺栓的抗剪承载力设计值为：

(2)

式中，n_f为传力摩擦面数目；P_d为一个高强螺栓的预拉力；μ为摩擦面的抗滑移系数；N_vd^b=0.9×2×0.45×80=64.8 kN。模型中有9个螺栓，全部螺栓抗剪承载能力设计值为N_G=583.2 kN，当外荷载超过N_G时，所有栓杆和孔壁接触模型进入弹塑性变化阶段。文献[15]中试件Ⅰ为单侧3行3列共9个螺栓，对多于3行的螺栓连接，具有相同的传力规律，本研究限于篇幅，只给出3行螺栓连接研究结果。

2.1 传力比分析

在螺栓拼接接头中每列螺栓传递的外力与全部螺栓所传递的总外力之比称为传力比。各列螺栓传力比-荷载曲线如图 9所示。

(1) 摩擦阶段，即螺栓松动之前，接头第1，2，3列螺栓传力比分别为0.35，0.31，0.33，表现为首列、末列螺栓传力比较大，中间列螺栓传力比较小的趋势。

(2) 滑移阶段，当外荷载达到375 kN时传力比最大的第1列螺栓首先开始松动，之后出现微小滑移，记为N_F=375 kN。随着荷载增大，出现滑移的螺栓所传递的剪力值逐渐减小，第2列螺栓所传递的剪力值逐渐增大，各列螺栓的传力比向同一数值逼近。

(3) 弹塑性阶段，当外荷载达到N_G附近时所有栓杆和孔壁接触，各列螺栓的传力比达到同一数值。随着外荷载继续增大，螺杆与孔壁接触挤压，进入弹塑性阶段，各列螺栓的传力比由同一数值逐渐出现和摩擦阶段相同的分布规律，即第1列和第3列螺栓传力比较大，第2列传力比较小，并且随着荷载的增大，各列螺栓的传力比逐渐均匀。第2列螺栓传力比逐渐增大，第1列螺栓传力比逐渐减小，当外荷载达到905 kN时，第1列螺栓传力比到达最小值，此时第1列螺栓位置芯板区域大面积进入塑性，可以认为芯板已经达到极限承载能力状态，之后随着第1列螺栓位置芯板区域塑性进一步发展孔壁发生挤压坡坏。

2.2 芯板应力分析

根据各列螺栓传力比变化规律，图 10给出芯板在摩擦滑移阶段、第1列螺栓开始松动、全部螺栓达到最大滑移量和芯板达到极限承载能力状态下的受拉方向(轴向)应力云图。芯板两块钢板应力以拼接中心呈对称分布趋势，证明本研究取受拉侧研究是可列的。受拉侧芯板轴向应力分布具有以下特征：

(1) 芯板各列螺栓孔受拉力一侧的孔壁，从第1列至第3列拉应力逐渐减小，芯板各列螺栓孔固定侧孔壁，从第1列至第3列压应力逐渐增大。由此可见，芯板内力是逐渐传递到各列螺栓的。

(2) 芯板应力沿板横向呈波浪形分布，在沿轴向螺栓连线上应力较小而在螺栓列间较大。

(3) 第1列螺栓开始松动之前，芯板应力较小，第1列螺栓松动时孔壁附近最大应力为142.7 MPa。所有螺栓达到最大滑移量时，即栓杆与孔壁接触后。第1列螺栓孔壁附近最大应力245.5 MPa，可见当外荷载达到N_G时，第1列螺栓孔壁附近应力发生较大变化。

(4) 当外荷载达到905 kN时，图 10(d)所示芯板第1列螺栓截面大部分进入塑性，发生孔壁挤压破坏。

2.3 盖板应力分析

由于上下盖板受力情况完全相同，图 11给出上盖板在摩擦滑移阶段、第1列螺栓开始松动、全部螺栓达到最大滑移量和芯板达到极限承载能力状态下的受拉方向(轴向)应力云图。可见盖板应力以拼接中心呈对称分布趋势，受拉侧右半边盖板轴向应力分布具有以下特征：

(1) 与芯板应力相反，盖板各列螺栓孔受拉力一侧孔壁，从第1列至第3列压应力逐渐减小，盖板各列螺栓孔固定侧孔壁，从第1列至第3列拉应力逐渐增大。

(2) 盖板应力沿板横向呈波浪形分布，在沿轴向螺栓连线上应力较小而在螺栓列间较大。这与芯板应力沿板横向分布规律一致。

(3) 随着外荷载的增大，盖板应力逐渐增大，相对于芯板来说，各个阶段应力值均较小。当外荷载达到905 kN时，由图 11(d)所示，第3列8号螺栓孔壁局部最大应力值为253.7 MPa，而查看Von Mises应力为231.2 MPa，这说明盖板没有进入屈服状态。

3 多行螺栓滑移曲线的模拟 3.1 多行螺栓抗剪滑移曲线

图 12为简化螺栓抗剪滑移模型芯板受拉端荷载位移曲线，根据曲线特征，整个模型经历了3个阶段：摩擦阶段、滑移阶段和弹塑性阶段。在摩擦阶段，由各列螺栓传力比分析可知外荷载N小于N_F时，荷载与受拉端位移为线性关系；在滑移阶段，当外荷载N大于N_F而小于N_G时，根据各列螺栓传力比变化规律，传力比较大的个别螺栓首先达到抗剪承载能力设计值N_vd^b，螺栓松动，发生微小滑移，此时大部分螺栓并未出现滑移，荷载与受拉端位移基本保持线性关系，当外荷载N到达N_G附近时，所有螺栓发生滑移并达到最大滑移量；当外荷载N大于N_G时，进入弹塑性阶段。

在忽略连接钢板本身在外力作用下的伸缩变形以及栓杆的弯曲变形的前提下，根据文献[17]中单个高强螺栓抗剪滑移计算公式，结合本研究中简化模型计算结果，多行高强螺栓抗剪滑移(总滑移量)可以分成3段按下式进行模拟：

(3)

(4)

式中，Δ_s0为螺栓孔的间隙；d₁和d₂分别为螺栓孔直径和螺杆直径；S₁为芯板螺栓孔壁变形均值；S₂为盖板螺栓孔壁变形均值；N为施加于受拉侧外力，对于3行螺栓连接N_F=375 kN为第1个螺栓开始松动时的外力；N_G为螺栓群抗剪承载能力设计值，N_G=583.2 kN。对于3行及以上螺栓连接，式(3)同样适用，螺栓连接个数不同，芯板和盖板孔壁变形均值会有差异，各列螺栓的芯板和盖板孔壁变形规律基本一致。

3.2 芯板和盖板孔壁变形分析

相对于单个螺栓来说^[18]，多行螺栓抗剪滑移分析比较复杂，假设在摩擦阶段和滑移阶段孔壁不发生变形，当栓杆与孔壁接触后开始变形。螺栓孔壁沿受拉方向的变形可定义为：

(5)

式中u_节点1和u_节点1为图 13中节点1和2沿受拉方向的位移。芯板各螺栓孔壁荷载变形曲线如图 14(a)所示，盖板各螺栓孔壁荷载变形曲线如图 14(b)所示。

由图 14可见，芯板各螺栓孔壁变形与盖板各螺栓孔壁变形规律基本一致，当荷载较小时，孔壁变形和荷载保持线性关系，当荷载超过某个确定值时，孔壁变形和荷载为非线性关系。图 14(a)中第1列1^#~3^#螺栓孔壁变形趋势基本一致，1^#和3^#螺栓同步变形，2^#螺栓孔稍有差别。第2~3列4^#，6^#，7^#，9^#螺栓孔壁变形基本一致，第2列5^#螺栓和第3列8^#螺栓孔壁变形与其他螺栓差别较大。2^#，5^#，8^#螺栓为同一列螺栓中处于中间位置的螺栓，这也说明，在同一列螺栓中处于中间位置的螺栓分配到的剪力较小，从第1列到第3列这种现象越来越明显。图 14(b)中，盖板各螺栓孔壁变形与芯板相比，变形量较小，可见螺栓拼接接头中芯板首先发生屈服。

图 14(a)中，当荷载小于704.5 kN时，荷载与变形保持线性关系，704.5 kN以后表现为非线性关系。为了简单实用，采用4次多项式进行拟合已足够精确。芯板孔壁变形S₁与外荷载N的关系式为：

(6)

图 14(b)中, 当荷载小于884.5 kN时，荷载与变形保持线性关系，884.5 kN以后表现为非线性关系，采用3次多项式进行拟合已足够精确。盖板孔壁变形S₂与外荷载N关系式为：

(7)

将式(6)对应的S₁-N曲线与芯板各螺栓孔壁变形平均值有限元结果进行对比，如图 15(a)所示，可见二者符合较好。将式(7)对应的S₂-N曲线与盖板各螺栓孔壁变形平均值有限元结果进行对比，如图 15(b)所示，可见二者符合较好。当外荷载大于905 kN时，拟合公式仅有参考意义。

多行螺栓连接中各行芯板和盖板孔壁变形规律具有一般性规律，但是对于不同行列的螺栓连接来说，每个螺栓的芯板和盖板孔壁变形具有差异，本研究中拟合公式仅适用3行3列9个螺栓连接的情况。

4 结论

通过简化的螺栓滑移数值计算方法可以精确模拟多行螺栓抗剪滑移过程，得到以下结论：

(1) 多行螺栓的抗剪滑移过程可以分为摩擦阶段、滑移阶段、弹塑性变形阶段。摩擦阶段中，各列螺栓所传递的力并不均匀，第1列和最末列传递力较大，中间列较小；滑移阶段中，各列螺栓所传递的力逐渐变化，在滑移阶段末尾达到均匀分布；弹塑性变形阶段中，栓杆与孔壁接触挤压，各列螺栓又表现为摩擦阶段的传力特征，随着外力的增大，各列螺栓所传递的力趋于均匀。

(2) 芯板应力，从第1列到第末列拉应力逐渐减小，压应力逐渐增大。盖板应力从第1列到第末列压应力逐渐减小，拉应力逐渐增大。当外荷载较大时，芯板第1列螺栓截面大部分进入塑性，发生孔壁挤压破坏。

(3) 根据三阶段滑移过程，得出螺栓群总滑移量计算公式，该公式适用于多行多列螺栓连接。根据孔壁变形特点，采用分段函数拟合出3行3列9个螺栓连接的S₁-N曲线和S₂-N曲线，并与有限元结果进行对比，二者符合较好。S₁-N曲线和S₂-N曲线仅适用于3行3列螺栓连接，对于多行多列螺栓连接仅供参考。