螺纹连接是工程结构中最常用的紧固方式之一，作为连接部件，螺栓螺母的紧固能力关系着包括被连接件在内整个结构的安全性与可靠性^[1−3]。其松动脱落等力学问题一直是国内外相关学者研究重点之一^[4−6]。目前对螺纹连接力学性能及松动机理的研究方法主要有实验检测、理论分析和数值计算。

在对螺纹结构进行实验检测时，由于螺纹加工精度和表面粗糙度不同，导致实验结果较为分散，实验重复性差^[7−8]。

随着计算机技术的发展，国内外学者运用有限元法对螺纹连接结构的性能进行了广泛研究，并取得了一定成果。在用有限元法分析螺纹连接的力学行为时，由于螺纹升角结构复杂，难以划分优质的有限元网格。为了简化计算，许多研究人员通常假定有限元模型的螺纹部分具有轴对称几何特性或直接忽略模型的螺纹结构。如：Verwaerde等^[9]采用3维实体单元模型对螺栓连接结构进行了有限元分析，虽然在一定程度上表征了螺纹连接结构上的应力分布，但模型忽略了连接处的螺纹，无法反映螺纹连接处的真实应力分布。Noda等^[10]利用轴对称模型进行螺栓连接结构的应力分析，得到了相对满意的结果，但轴对称模型虽然考虑了螺纹的作用，却忽略了螺纹升角对结果的影响。

使用带升角结构的精确螺纹有限元模型开展数值计算，是准确预测螺纹处3维应力分布的有效方法，但由于工业生产中实际的需求不同，使用到的螺栓螺母紧固件种类繁多。因此，在进行有限元数值仿真计算时，如果需要针对不同型号的螺纹结构进行有限元数值研究，建立模型和网格划分将耗费大量精力，同时将出现网格质量难以保证、数值计算不易收敛、仿真成本高等问题^[11−13]。故寻找一种计算结果准确、建模便捷的螺纹紧固件有限元建模方法具有十分重要的意义。

本文采用ISO给定的标准螺旋线公式建立垂直于螺纹中心轴线的螺纹螺旋线几何形状，以此螺纹螺旋线几何形状为基本几何，将其分别沿螺纹径向和轴向，顺次构建等比例缩放和固定角度旋转的平面螺纹螺旋线相似几何；然后，将单个节距螺纹上所得到的全部螺旋线等间距离散成若干节点，并将这些离散节点按位置关系组合成精度较高的3维8节点（C3D8）单元，即获得单个节距螺纹3维有限元模型；再将此单个节距的模型，沿螺纹中心轴线堆叠适当个数，得到多节距的螺纹连接结构3维有限元模型（该建模具体实施过程可查阅文献[14]）。最后，通过数值算例对建立的模型进行了验证，并在Qt框架下开发螺栓螺母有限元参数化建模平台，为实际工程提供了更为科学的支撑和更高效的设计工具。

1 有限元模型

本文螺纹轮廓线的规格由文献[15]给出，为了防止螺纹齿根处发生过度的应力集中现象，要求齿根半径不得小于0.125 $P$ ， $ P $ 为节距。图1为沿螺栓中心轴线的螺纹截面，其中， $ r $ 为径向坐标， $ \theta $ 为周向坐标， $ {d_1} $ 为螺纹小径， $ d $ 为螺纹公称直径， $ H $ 为螺纹工作高度，A、B、C、D和A′、B′、C′、D′分别为轮廓线上关于 $ \theta = 0 $ 相互对称的关键点。由A点开始，从 $ \theta = 0 $ 到 $ \theta = \text{π}$ 的螺栓外螺纹轮廓线可以被分为3个部分，即0～ $ {\theta _1} $ 对应AB （螺纹齿根）， $ {\theta _1} $ ～ ${\theta _2} $ 对应BC（螺纹齿侧）， $ {\theta _2} $ ～ $\text{π} $ 对应CD（螺纹顶部）。

将这3部分的轮廓线延展成绕螺栓轴线的平面，投影到垂直于螺栓轴的平面上，可得外螺纹轮廓线，如图2所示。

此轮廓线轨迹 ${r_1} $ 可由式（1）给出：

$ {r_1} = \left\{ \begin{array}{l} \dfrac{d}{2} - \dfrac{7}{8}H + 2\rho - \sqrt {{\rho ^2} - \dfrac{{{P^2}}}{{4{\text{π} ^2}}}{\theta ^2}} {\text{ }},0 \le \theta \le {\theta _1} ; \\ \dfrac{H}{\text{π} }\theta + \dfrac{d}{2} - \dfrac{7}{8}H{\text{ }},{\theta _1} \le \theta \le {\theta _2} ; \\ \dfrac{d}{2}{\text{ }},{\theta _2} \le \theta \le \text{π} \end{array} \right. $

(1)

式中： ${\theta _1} = \dfrac{{\sqrt 3 \text{π} }}{P}\rho$ ； ${\theta _2} = \dfrac{7}{8}\text{π}$ ； $H = \dfrac{{\sqrt 3 }}{2}P$ ； $\;\rho \le \dfrac{{\sqrt 3 }}{{12}}P$ ， $\;\;\rho$ 为外螺纹齿根半径的上限。

同样地，可获得螺母的内螺纹轮廓线轨迹 ${r_2} $ 的数学表达：

$ {r_2} = \left\{ \begin{array}{l} \dfrac{{{d_1}}}{2},0 \le \theta \le {\theta _1} ;\\ \dfrac{H}{\text{π} }\theta + \dfrac{d}{2} - \dfrac{7}{8}H,{\theta _1} \le \theta \le {\theta _2} ;\\ \dfrac{d}{2} + \dfrac{1}{8}H - 2{\rho _{\text{n}}}{\text{ + }} \hfill \\ \quad \sqrt {{\rho _{\text{n}}}^2 - \dfrac{{{P^2}}}{{4{\text{π} ^2}}}{{(\text{π} - \theta )}^2}} ,{\text{ }}{\theta _2} \le \theta \le \text{π} \hfill \\ \end{array} \right. $

(2)

式中： ${\theta _1} = \dfrac{\text{π} }{4}$ ； ${\theta _2} = \text{π} \left(1 - \dfrac{{\sqrt 3 {\rho _{\text{n}}}}}{P}\right)$ ； ${\;\rho _{\text{n}}} \le \dfrac{{\sqrt 3 }}{{24}}P$ ， $\;{\rho _{\text{n}}}$ 为内螺纹齿根半径的上限。

为了实现对不同型号螺纹连接结构有限元模型的便捷调用，在Qt平台上，利用OpenGL图形库开发了模型建立软件，软件界面如图3所示。

图4为单位节距长度 $P$ =2、螺纹轴向单元数 $ne$ =96、匹配螺杆处单元层数 $t$ =0、螺纹向内层数 $j$ =3、节距数量 $n$ =6条件下，螺栓（LS）螺母（LM）公称直径 $d$ 分别为10 mm、16 mm和20 mm的有限元模型对比图。

2 螺纹连接中的力学关系

螺纹连接结构拧紧和受载时，力学行为的变化是一个强烈的状态非线性过程，该过程通常伴有复杂的几何非线性（大位移）和材料非线性（塑性）影响；螺纹连接结构处于不同状态时，其力学行为有着明显的差别。为便于计算分析，采用下述假定：

1）螺纹结构体的材料是连续均匀且完全弹性的；

2）作用在螺纹连接结构上的位移和形变是微小的，且结构体内无初始应力；

3）摩擦系数在整个运动受力过程中不发生变化，系统能量守恒。

2.1 螺栓轴向力与拧紧转角关系

从螺纹的几何特征入手，当螺母绕螺栓转动时，螺母的轴向位移量 $ s $ 为：

$ s = P\vartheta /2\text{π} $

(3)

式中， $\vartheta$ 为螺母转角， $ P $ 为节距。如果对螺母的位移做刚性的约束，同时限制螺栓位移，则转动螺母时螺栓将伸长变形，其伸长量 $ {\delta _{\rm{b}}} $ 为：

$ {\delta _{\rm{b}}} = P\vartheta /(2\text{π}) $

(4)

然而，实际上，被连接件不可能刚性约束螺母的位移，所以在转动螺母，螺母发生位移的同时，螺栓也被拉长，螺栓和被连接件两者被拉长的占比，取决于二者刚度 $ {K_{\rm{b}}} $ 和 $ {K_{\rm{j}}} $ 。转动螺母一方面压缩被连接件，一方面拉伸螺栓，相当于串联两个弹簧。串联弹簧的系统刚度为 $ {K_{\rm{s}}} = {K_{\rm{b}}}{K_{\rm{j}}}/({K_{\rm{b}}} + {K_{\rm{j}}}) $ 。由胡克定律中形变和力的关系，可得螺栓上的轴向力 $ F $ 为：

$ F = {K_{\rm{s}}}P\vartheta /2\text{π} {\text{ }} = {K_{\rm{b}}}{K_{\rm{j}}}P\vartheta /[2\text{π} ({K_{\rm{b}}} + {K_{\rm{j}}})] $

(5)

实际拧紧一个螺母时，在螺母和被连接件未发生接触的旋转期间，螺栓上不会产生轴向力。一旦螺母接触被连接件后，即会产生轴向力，由于表面接触区域面积较小或被连接件与周围构件的摩擦力影响等原因，轴向力起初会很小，但后期增长很迅速，这个过程称为被连接件的贴紧过程。贴紧过程要转多少角度无法预知，即使使用同种型号的螺栓和被连接件，贴紧过程中螺母转动的角度也会不一样^[16]。被连接件贴紧后，从式（5）可知轴向力 $ F $ 与转角 $ \vartheta $ 呈线性关系。继续转动螺母，螺栓继续按比例伸长，当达到材料的屈服点后，同样的转角增量下，螺栓伸长量的增量增大，轴向力增量减小，轴向力 $ F $ 和转角 $ \vartheta $ 由直线关系变成曲线关系。图5给出了拧紧螺母全过程的螺栓轴向力 $ F $ 和 $ \vartheta $ 之间的关系曲线（OQ 为“空转”螺栓，QR为贴紧过程，RS为线性段，在S点开始屈服，ST为屈服后的曲线段）。

2.2 螺纹牙的轴向力分布推导

图6为旋和螺纹部分的轴向力F示意图，其中， $ L $ 为内外螺纹啮合长度。螺栓杆受拉伸力 $ {F_{\rm{b}}} $ 的作用，如以螺母顶面位置为原点，在 $ x $ 位置处将产生一个垂直截面的轴向力 $ F $ ，则在 $ x $ 位置将产生一个垂直于螺纹啮合面的单位力 $ p $ ，此时，在 $ x $ 位置处螺栓的伸长长度 $ {\varepsilon _{\rm{b}}} $ 和螺母的缩短长度 $ {\varepsilon _{\rm{n}}} $ 可由式（6）求出：

$ \left\{ {\begin{array}{l} {{\varepsilon _{\rm{b}}} = \displaystyle\int_0^x {\dfrac{F}{{{A_{\rm{b}}} \cdot {E_{\rm{b}}}}}} {\rm{d}}x}, \\ {{\varepsilon _{\rm{n}}} = \displaystyle\int_0^x {\dfrac{F}{{{A_{\rm{n}}} \cdot {E_{\rm{n}}}}}} {\rm{d}}x} \end{array}} \right. $

(6)

式中， $ {A_{\rm{b}}} $ 和 $ {A_{\rm{n}}} $ 分别为螺栓和螺母的垂直截面面积， $ {E_{\rm{b}}} $ 和 $ {E_{\rm{n}}} $ 分别为螺栓和螺母纵向的弹性模量。

设作用于 $ x $ 与 $ x + {\rm{d}}x $ 之间的某螺纹牙上的轴向力为 $ {\rm{d}}F $ ，则有：

$ p\cos \;\alpha \cdot {\rm{d}}x\cot\; \varphi = {\rm{d}}F $

(7)

式中， $ \varphi $ 为螺纹升角， $\alpha$ 为螺纹斜角，由式（7）可得：

$ p\cos \;\alpha = \frac{{{\rm{d}}F}}{{{\rm{d}}x}}\tan \;\varphi $

(8)

图7描述了由各种原因产生的螺纹牙变形。如图7（a）所示：对牙宽为 $a $ 的梯形截面横梁，在 $ y = c $ ，跨度为 $b $ 的位置上，每单位宽度作用一个大小为 $ p $ 的力。基于山本晃^[17]对无摩擦螺纹模型弹性变形的理论分析可推导出在 $ y = c $ 位置处横梁弯曲引起的变形长度 $ {\delta _1} $ 、剪力引起的变形长度 $ {\delta _2} $ 、牙根倾斜引起的变形长度 $ {\delta _3} $ 、牙根剪切引起的变形长度 $ {\delta _4} $ 和径向分力引起的变形长度 $ {\delta _5} $ 。

把上述所求得的变形分量相互累加，即获得外螺纹和内螺纹变形的总量 $ {\delta _{\rm{b}}} $ 和 $ {\delta _{\rm{n}}} $ ，分别为：

$ \left\{ {\begin{array}{l} {{\delta _{\rm{b}}} = {\delta _{1{\rm{b}}}} + {\delta _{2{\rm{b}}}} + {\delta _{3{\rm{b}}}} + {\delta _{4{\rm{b}}}} + {\delta _{5{\rm{b}}}} = \dfrac{{{k_{\rm{b}}}}}{{{E_{\rm{b}}}}}p\cos\;\alpha}, \\ {{\delta _{\rm{n}}} = {\delta _{1{\rm{n}}}} + {\delta _{2{\rm{n}}}} + {\delta _{3{\rm{n}}}} + {\delta _{4{\rm{n}}}} + {\delta _{5{\rm{n}}}} = \dfrac{{{k_{\rm{n}}}}}{{{E_{\rm{n}}}}}p\cos\;\alpha} \end{array}} \right. $

(9)

式中， ${k_{\rm{b}}} $ 和 ${k_{\rm{n}}} $ 分别为与材料弹性模量和泊松比相关的系数。

将式（8）代入式（9），得：

$ \left\{ {\begin{array}{l} {{\delta _{\rm{b}}} = \dfrac{{{k_{\rm{b}}}}}{{{E_{\rm{b}}}}} \cdot \tan\;\varphi \dfrac{{{\rm{d}}F}}{{{\rm{d}}x}}}, \\ {{\delta _{\rm{n}}} = \dfrac{{{k_{\rm{n}}}}}{{{E_{\rm{n}}}}} \cdot \tan\;\varphi \dfrac{{{\rm{d}}F}}{{{\rm{d}}x}}} \end{array}} \right. $

(10)

取图6中 $ L = x $ ，在 $ {\varepsilon _{\rm{b}}} $ 、 $ {\varepsilon _{\rm{n}}} $ 、 $ {\delta _{\rm{b}}} $ 、 $ {\delta _{\rm{n}}} $ 之间可得出如下关系：

$ {({\varepsilon _{\rm{b}}} + {\varepsilon _{\rm{n}}})_{x = x}} = {({\delta _{\rm{b}}} + {\delta _{\rm{n}}})_{x = x}} - {({\delta _{\rm{b}}} + {\delta _{\rm{n}}})_{x = 0}} $

(11)

把式（7）和（10）代入式（11），并对 $ x $ 求微分，可得：

$ \frac{{{{\rm{d}}^2}F}}{{{\rm{d}}{x^2}}} = {\lambda ^2}F $

(12)

式（12）中：

$ \lambda = \sqrt {\dfrac{{\dfrac{1}{{{A_{\rm{b}}}{E_{\rm{b}}}}} + \dfrac{1}{{{A_{\rm{n}}}{E_{\rm{n}}}}}}}{{\left(\dfrac{{{k_{\rm{b}}}}}{{{E_{\rm{b}}}}} + \dfrac{{{k_{\rm{n}}}}}{{{E_{\rm{n}}}}}\right)\tan\; \varphi }}} = \sqrt {\dfrac{{\dfrac{1}{{{A_{\rm{b}}}}} + \dfrac{{{E_{\rm{b}}}}}{{{E_{\rm{n}}}}}\dfrac{1}{{{A_{\rm{n}}}}}}}{{\left({k_{\rm{b}}} + \dfrac{{{E_{\rm{b}}}}}{{{E_{\rm{n}}}}}{k_{\rm{n}}}\right)\tan\; \varphi }}} $

(13)

式（12）的通解为：

$ F = {c_1}\cosh (\lambda x) + {c_2}\sinh (\lambda x) $

(14)

将边界条件 ${F_{x = 0}} = 0$ ， ${F_{x = L}} = {F_{\rm{b}}}$ 代入式（12）得：

$ \left\{ {\begin{array}{l} {{c_1}\cosh \;0 + {c_2}\sinh \;0 = 0}, \\ {{c_1}\cosh (\lambda L) + {c_2}\sinh (\lambda L) = {F_{\rm{b}}}} \end{array}} \right. $

(15)

式中， ${c_1} $ 和 ${c_2} $ 为方程待定系数。

双曲正弦和双曲余弦的函数式为：

$ \left\{ \begin{array}{l} \sinh \;x = \dfrac{{{{\rm{e}}^x} - {{\rm{e}}^{ - x}}}}{2}, \\ \cosh \;x = \dfrac{{{{\rm{e}}^x} + {{\rm{e}}^{ - x}}}}{2} \end{array} \right. $

(16)

联立式（15）与（16）可得：

$ \left\{ \begin{array}{l} {c_1} = 0,\\ {c_2} = \dfrac{{{F_{\rm{b}}}}}{{\sinh (\lambda L)}} \end{array} \right. $

(17)

将 $ {c_1} $ 、 $ {c_2} $ 的值代入式（14），则式（14）变为：

$ \frac{F}{{{F_{\rm{b}}}}} = \frac{{\sinh (\lambda x)}}{{\sinh (\lambda L)}} $

(18)

将式（18）代入式（8）可得：

$ p\cos \;\alpha = \lambda \frac{{\cosh (\lambda x)}}{{\sinh (\lambda L)}} \cdot {F_{\rm{b}}}\tan\; \varphi $

(19)

第 $ i $ 圈螺纹承受的轴向载荷 $F_i$ 为：

$ \begin{aligned}[b] F_i =& \int_{(i - 1)P}^{iP} {p\cos \;\alpha \cot \;\varphi {\rm{d}}x} =\\& \int_{(i - 1)P}^{iP} {\lambda \frac{{\cosh (\lambda x)}}{{\sinh (\lambda L)}}{F_{\rm{b}}}{\rm{d}}x},i = 1,2,3, \cdots, L/P \end{aligned}$

(20)

定义 $ {\eta _i} $ 为第 $ i $ 圈工作螺纹所承受的轴向载荷占总轴向力 $ {F_{\rm{b}}} $ 的比例，则有：

$ {\eta _i} = \frac{{F_i}}{{{F_{\rm{b}}}}} = \left. {\frac{{\sinh (\lambda x)}}{{\sinh (\lambda L)}}} \right|_{(i - 1)P}^{iP} $

(21)

式（21）表明，对于螺纹结构式连接，轴向力在旋合各圈螺纹牙间所占比例满足双曲正弦函数关系，其轴向载荷比例分布如图8所示。

以M12标准螺栓螺母为研究对象，对螺栓轴向力分布进行计算。螺纹的几何参数见表1，螺栓、螺母弹性材料参数见表2，塑性应变行为见表3。

根据表2：取螺栓材料35CrMn，其弹性模量 ${E_{\rm{b}}}$ 为213000 MPa，泊松比 ${\nu _{\rm{b}}}$ 为0.286；取螺母材料45钢，其弹性模量 ${E_{\rm{n}}}$ 为209000 MPa，泊松比 ${\nu _{\rm{n}}}$ 为0.269。M12螺纹中径 ${d_{\rm{p}}}$ 、螺母大径 ${D_{\rm{O}}}$ 取值见表1。可求得外螺纹螺栓的变形分量分别为：

$ \left\{ {\begin{array}{l} {{\delta _{1{\rm{b}}}} = 0.034\left(\dfrac{{1 - {\nu _{\rm{b}}}^2}}{{{E_{\rm{b}}}}}\right)p\cos \;\alpha = 0.031\;22\dfrac{p}{{{E_{\rm{b}}}}}\cos \;\alpha}, \\ {{\delta _{2{\rm{b}}}} = 1.08\left(\dfrac{{1 + {\nu _{\rm{b}}}}}{{{E_{\rm{b}}}}}\right)p\cos \;\alpha = 1.388\;88\dfrac{p}{{{E_{\rm{b}}}}}\cos \;\alpha}, \\ {{\delta _{3{\rm{b}}}} = 0.229\left(\dfrac{{1 - {\nu _{\rm{b}}}^2}}{{{E_{\rm{b}}}}}\right)p\cos \;\alpha = 0.210\;27\dfrac{p}{{{E_{\rm{b}}}}}\cos \;\alpha}, \\ {{\delta _{4{\rm{b}}}} = 1.18\left(\dfrac{{1 - {\nu _{\rm{b}}}^2}}{{{E_{\rm{b}}}}}\right)p\cos \;\alpha = 1.083\;48\dfrac{p}{{{E_{\rm{b}}}}}\cos \;\alpha}, \\ {\delta _{5{\rm{b}}}} = {u_{\rm{b}}}\tan \;\alpha = (1 - {v_{\rm{b}}})\dfrac{{{{\tan }^2}\alpha}}{2} \cdot \dfrac{{{d_{\rm{p}}}}}{P} \cdot \dfrac{{p\cos \;\alpha}}{{{E_{\rm{b}}}}} {\text{ }} = \\\qquad\; 0.738\;68\dfrac{p}{{{E_{\rm{b}}}}}\cos \;\alpha \end{array}} \right. $

(22)

同样地，可求得内螺纹螺母的变形分量分别为：

$ \left\{ {\begin{array}{l} {{\delta _{1{\rm{n}}}} = 0.073\left(\dfrac{{1 - {\nu _{\rm{n}}}^2}}{{{E_{\rm{n}}}}}\right)p\cos \;\alpha = 0.067\;72\dfrac{p}{{{E_{\rm{n}}}}}\cos \;\alpha}, \\ {{\delta _{2{\rm{n}}}} = 1.15\left(\dfrac{{1 + {\nu _{\rm{n}}}}}{{{E_{\rm{n}}}}}\right)p\cos \;\alpha = 1.459\;35\dfrac{p}{{{E_{\rm{n}}}}}\cos \;\alpha}, \\ {{\delta _{3{\rm{n}}}} = 0.294\left(\dfrac{{1 - {\nu _{\rm{n}}}^2}}{{{E_{\rm{n}}}}}\right)p\cos \;\alpha = 0.272\;73\dfrac{p}{{{E_{\rm{n}}}}}\cos \;\alpha} ,\\ {{\delta _{4{\rm{n}}}} = 1.14\left(\dfrac{{1 - {\nu _{\rm{n}}}^2}}{{{E_{\rm{n}}}}}\right)p\cos \;\alpha = 1.057\;51\dfrac{p}{{{E_{\rm{n}}}}}\cos \;\alpha}, \\ {\delta _{5{\rm{n}}}} = {u_{\rm{n}}}\tan\;\alpha = \left(\dfrac{{{D_{\rm{O}}}^2 + {d_{\rm{p}}}^2}}{{{D_{\rm{O}}}^2 - {d_{\rm{p}}}^2}} + {\nu _{\rm{n}}}\right)\dfrac{{{{\tan }^2}\alpha}}{2} \cdot \dfrac{{{d_{\rm{p}}}}}{P} \cdot \dfrac{{p\cos \;\alpha}}{{{E_{\rm{n}}}}} = \\ \qquad\; 2.306\;20\dfrac{p}{{{E_{\rm{n}}}}}\cos \;\alpha \end{array}} \right. $

(23)

将式（22）和（23）代入式（9），得：

$ \left\{ {\begin{array}{l} {\delta _{\rm{b}}} = {\delta _{1{\rm{b}}}} + {\delta _{2{\rm{b}}}} + {\delta _{3{\rm{b}}}} + {\delta _{4{\rm{b}}}} + {\delta _{5{\rm{b}}}}= \\ \qquad \dfrac{{{k_{\rm{b}}}}}{{{E_{\rm{b}}}}}p\cos \;\alpha = 3.452\;53\dfrac{p}{{{E_{\rm{b}}}}}\cos \;\alpha, \hfill \\ {\delta _{\rm{n}}} = {\delta _{1{\rm{n}}}} + {\delta _{2{\rm{n}}}} + {\delta _{3{\rm{n}}}} + {\delta _{4{\rm{n}}}} + {\delta _{5{\rm{n}}}} = \\ \qquad \dfrac{{{k_{\rm{n}}}}}{{{E_{\rm{n}}}}}p\cos \;\alpha = 5.163\;51\dfrac{p}{{{E_{\rm{n}}}}}\cos \;\alpha \end{array}} \right. $

(24)

由式（24）可得， ${k_{\rm{b}}} = 3.452\;53$ ， ${k_{\rm{n}}} = 5.163\;51$ 。将表1中螺栓、螺母的垂直截面面积 $ {A_{\rm{b}}} $ 、 $ {A_{\rm{n}}} $ ，螺纹的螺旋升角 $ \varphi $ 和求得的 $ {k_{\rm{b}}} $ 、 $ {k_{\rm{n}}} $ 值代入式（13），得：

$ \begin{gathered} \lambda = \sqrt {\dfrac{{\dfrac{1}{{{A_{\rm{b}}}{E_{\rm{b}}}}} + \dfrac{1}{{{A_{\rm{n}}}{E_{\rm{n}}}}}}}{{\left(\dfrac{{{k_{\rm{b}}}}}{{{E_{\rm{b}}}}} + \dfrac{{{k_{\rm{n}}}}}{{{E_{\rm{n}}}}}\right)\tan\; \varphi }}} = \sqrt {\dfrac{{\dfrac{1}{{{A_{\rm{b}}}}} + \dfrac{{{E_{\rm{b}}}}}{{{E_{\rm{n}}}}}\dfrac{1}{{{A_{\rm{n}}}}}}}{{\left({k_{\rm{b}}} + \dfrac{{{E_{\rm{b}}}}}{{{E_{\rm{n}}}}}{k_{\rm{n}}}\right)\tan \;\varphi }}} {\text{ }} = 0.196\;11 \end{gathered} $

(25)

取螺母节距个数为6，单位长度为1.75 mm，则螺母总高度L=10.5 mm。因此，将以上结果代入式（18）和（21），可以获得第 $ i $ 圈工作螺纹所承受的轴向载荷 $ {F_i} $ 与总轴向力 $ {F_{\rm{b}}} $ 的比值 $ {F_i}/{F_{\rm{b}}} $ ；通过求相临螺纹之间 $ {F_i}/{F_{\rm{b}}} $ 的差值，得各圈螺纹承受轴向力的占比 $ {\eta _i} $ ，如表4所示。

图9展示了沿螺母轴线向下，第 $ i $ 圈螺纹上轴向力与总轴向力之比 $ {\eta _i} $ 。

由图9看出，第1圈承受了约30.2%的轴向力，第2圈承受了约21.9%的轴向力，前3圈螺纹约承受全部轴向力的68.3%。

3 有效性验证

选用性能等级为5.8的M12粗牙标准螺纹紧固件有限元模型，在商业有限元计算软件ABAQUS中进行数值计算，将数值结果与第2节理论计算结果和以往工程实验结果进行对比分析，验证该螺纹连接结构有限元建模方法的有效性。

建立的整体模型结构如图10所示，共包含上层板、下层板、螺母和螺栓4个部件。这些部件全部采用C3D8R单元划分，模型共计节点87245个，单元77712个。

在部件接触界面之间定义以下4组接触副：螺栓头部下表面和上层板的上表面之间（Cont1）、上层板的下表面和下层板的上表面之间（Cont2）、下层板的下表面和螺母的上表面之间（Cont3）、螺栓外螺纹表面和螺母内螺纹表面之间（Cont4），将以上4组接触副全部定义为表面与表面约束（standard）。接触行为设为切向和法向，其中，螺栓与螺母表面法向行为设为Exponential软接触约束，其余接触表面的法向行为设为硬接触约束。切向行为均采用罚公式接触约束。通过在螺母外部加载扭转载荷，使螺母向上转动从而夹紧两块形状相同的板材。模型中螺栓、螺母为弹塑性材料，板为弹性材料。各部件弹性材料参数见表2，螺栓、螺母塑性应变行为见表3。

针对本模型，在ABAQUS中进行如下求解设置。建立静力通用分析步，加载时间为1 s。最大增量步数设为100，初始增量步为0.01，模型考虑受几何非线性的影响。场输出变量选择应力、应变等。对于螺栓轴向力加载，ABAQUS中提供了一种模拟螺栓轴向力的加载方式，即在螺栓截面上施加螺栓载荷（bolt load）。这种加载方式将导致在下一分析步做扭转激励加载时，螺栓轴向力始终维持螺栓载荷所施加的荷载值而不发生变化。显然这与实际情况并不相符。因此，本文提出在螺母外侧施加旋转角度，以此模拟螺栓的预紧过程。图11中分别为ABAQUS中螺栓载荷加载方式与本文提到的扭转螺母加载方式在扭转激励作用下轴向力随时间的变化曲线。对比两条曲线发现：1） 由于螺栓载荷忽略了螺栓在预紧过程所受到的摩擦力矩和部件间实际接触状态，预紧过程呈线性增长，与工程实际中螺栓预紧情况不符。2） 在1～2 s阶段的扭转激励作用下，螺栓载荷加载方式使得轴向力维持不变，无法真实反映轴向力随扭转载荷变化的机理。因此，本文选用更加贴近实际工程的扭转螺母加载方式来完成螺栓轴向力的加载。

3.1 拧紧转角和轴向力关系验证

在初始分析步（Step 0）中，对模型（图10）施加如下边界条件：在上层板上表面施加6个自由度的固定约束，下层板x方向侧面施加除z方向平移外的5个自由度固定约束，螺栓头部仅对绕z旋转的方向施加固定约束。

在维持初始分析步不变的前提下，在分析步（Step 1）中，对螺母施加0.3 rad转角，使结构处于拧紧状态。

在以上边界条件下，拧紧转角与轴向力的数值计算结果如图12所示，可以看到，拧紧转角和轴向力之间存在非线性关系。由于起初受被连接件与周围构件的摩擦力影响等原因，轴向力会很小，但轴向力的增长很迅速。当继续转动螺母，螺栓继续按比例伸长，此时，拧紧转角和轴向力呈线性关系。当达到材料的屈服点后，在同样的转角增量下，螺栓伸长量的增量增大，轴向力增量减小，轴向力 $ F $ 和拧紧转角 $ \vartheta $ 之间由直线关系变成曲线关系。数值计算的曲线变化规律与第2.1节中的理论分析相吻合。

3.2 轴向力分布验证

为重现第2.2节中推导的无摩擦模型^[17]在弹性理论下螺栓的轴向力分布结果，在不改变模型接触定义的情况下，将部件间的摩擦系数全设为0，并将图10中所有部件的塑性行为全部移除，其弹性材料参数见表2；对于性能等级为5.8的M12螺栓按照文献^[18]中的轴向力施加标准，不改变第3.1节Step 0中设置的初始边界条件，在Step 1中对螺母重新进行扭转预紧，使螺栓获得螺栓轴向力20 kN，得到M12螺栓轴向力分布Mises应力云图（判断共用节点是否平均取值的变量变化梯度限定值为75%，下文同），如图13所示。

在ABAQUS中创建自由体切片，获得图13中横截面序号1～7指向的各截面j的截面轴向力 ${F_j}$ ，并计算得到 ${F_j}$ 与横截面轴向合力的比值 ${F_j}/{F_{\rm{b}}}$ ，以及各圈螺纹承受轴向力的占比 ${\eta _i}$ ，见表5。图14对比了以上数值计算结果与第2.2节中理论推导出的计算结果。

从图14可知：相对轴向力的解析结果和数值模拟结果的最大误差发生在第5圈螺纹处，且误差小于5%；相对轴向力的解析结果和数值结果的平均误差在3%以内；相对轴向力的解析结果和数值结果变化的整体趋势具有一致性。

3.3 螺纹牙底应力集中系数验证

现有研究资料^[19-20]表明，螺纹连接结构最大应力发生在螺栓螺纹根部，位于一个螺距内的螺母受载表面。使用第1节中获得的带升角的螺纹有限元模型研究螺纹根部周围的应力集中情况。螺纹在受载时，连接件中的载荷分布非常不均匀，有些螺纹可能已经达到塑性阶段，而其他螺纹仅处于弹性状态。螺纹根部存在的严重应力集中往往会导致连接件疲劳失效。应力集中系数 $ {K_{\rm{t}}} $ 为缺口平面上局部峰值应力 ${\sigma _{\max }} $ 与名义应力 ${\sigma _{\rm{n}}}$ 之比，是衡量应力集中严重性的重要指标^[21]，表达式为：

$ {K_{\rm{t}}} = \frac{{{\sigma _{\max }}}}{{{\sigma _{\rm{n}}}}} $

(26)

由于螺纹紧固件的名义应力是螺栓在外力（轴向拉力）作用下缺口处产生的平均应力值，故可将式（26）改写为：

$ {K_{\rm{t}}} = \frac{{{\sigma _{\max }}}}{{{F_{{\rm{sum}}}}/A}} $

(27)

式中， $F_{{\rm{sum}}}$ 为横截面上的合力， $ A $ 为横截面面积。选取图15中序列号1′～6′指向的螺栓和螺母齿根位置处的应力峰值作为 $ {\sigma _{\max }} $ 的取值。

根据式（26）、（27），计算各螺纹圈上对应的螺纹齿根处应力集中系数 $ {K_{\rm{t}}} $ ，如图16所示。从图16可知，计算的应力集中系数 $ {K_{\rm{t}}} $ 沿螺纹轴线上背离连接件的方向上逐渐减小，该分布特点与文献[22]和[23]中的理论和数值解一致。结果表明：螺纹根部是最危险的应力集中区，峰值应力 $ {\sigma _{\max }} $ 首先发生在此区域；螺栓螺母啮合的第1圈工作螺纹应力集中系数最大，随着螺纹圈数的增加，应力集中系数逐渐减小，螺栓的应力集中值比螺母的应力集中值略大。

3.4 螺纹表面应力分布

在不改变模型其他设置的情况下，将螺栓/螺母的材料参数加入塑性应变，并重新进行数值计算，得到螺栓上Mises应力分布和最大主塑性应变分布，如图17、18所示。

由图17可知，前3圈螺纹承受载荷最大，且应力的较大值主要集中在螺纹的根部位置，这与第2.2节的理论推导结果及以往的实验结果^[24]相符合。由图18可知，塑性形变首先发生在螺栓的第1圈螺纹的根部，而实际工程中螺栓的疲劳断裂也常发生在第1圈螺纹的根部^[25-26]。

4 结　论

本文基于螺纹紧固件垂直于螺纹轴线的任意横截面形状相同的事实，选用ISO给定的标准螺旋线公式建立了完整的螺纹3维有限元模型。依据该建模方案，在Qt平台上借助OpenGL图形库，开发了一套螺栓螺母有限元参数化建模软件；通过此建模软件对螺纹结构进行参数控制，可以便捷地获取不同型号的螺纹连接件有限元模型。然后，通过数值实验的方法模拟了螺栓的拧紧过程，将数值实验结果与相应理论计算或以往研究进行对比分析，验证了本文建立的模型的有效性和保真性。