李玉奇^1,2, 罗忠^1,2, 栗江^1,2, 侯小捷^1,2
1. 东北大学 机械工程与自动化学院, 辽宁 沈阳 110819;
2. 东北大学 航空动力装备振动及控制教育部重点实验室, 辽宁 沈阳 110819
收稿日期：2018-04-20
基金项目：国家自然科学基金资助项目(11572082)；中央高校基本科研业务费专项资金资助项目(N170308028，N160312001)；辽宁省高等学校优秀人才支持计划项目(LJQ2015038)。
作者简介：李玉奇(1992-)，男，黑龙江嫩江人，东北大学博士研究生;
罗忠(1978-)，男，内蒙古商都人，东北大学教授，博士生导师。

摘要：转子系统中螺栓联接结构轴向联接刚度的不确定性对转子系统动力学特性具有重要影响, 为此, 在构建螺栓联接结构有限单元基础上, 建立了转子系统整体有限元模型, 采用非嵌入多项式混沌展开法分析轴向联接刚度不确定性对转子系统动力学特性的影响.结果表明:轴向联接刚度在一定范围内变化会导致临界转速及临界转速对应的稳态响应幅值偏离预期, 随着轴向刚度不确定性的标准差增大, 盘竖直方向稳态响应均值降低.研究结果可为螺栓联接转子的设计提供理论参考.
关键词：螺栓联接转子系统刚度不确定多项式展开动力学特性
Analysis of Bolted Joint Rotor System with Uncertain Axial Stiffness
LI Yu-qi^1,2, LUO Zhong^1,2, LI Jiang^1,2, HOU Xiao-jie^1,2
1. School of Mechanical Engineering & Automation, Northeastern University, Shenyang 110819, China;
2. Key Laboratory of Vibration and Control of Aero-Propulsion System, Ministry of Education, Northeastern University, Shenyang 110819, China
Corresponding author: LUO Zhong, E-mail: zhluo@mail.neu.edu.cn
Abstract: The stiffness uncertainty of bolted joint structure plays an important role in dynamic characteristics of the rotor system. The influences of the uncertainty of the axial stiffness on the vibration behaviors of the bolted joint rotor system are studied. Firstly, the motion equations of the bolted joint rotor system are established based on the finite element method, and then the stochastic dynamic model with uncertainty of axial stiffness of bolted joint structure is obtained by the nonintrusive polynomial expansion method. Finally, the dynamic characteristics is illustrated. The results show that the axial stiffness uncertainty will affect the critical speed and the steady response amplitude of the rotor system. The response amplitude of the disk in the vertical direction decreases as the increasing of axial stiffness uncertainty. The results can provide a theoretical reference for the design of the bolted joint rotor system.
Key words: bolted jointrotor systemstiffness uncertaintypolynomial expansiondynamic characteristics
螺栓联接因结构简单等特点而广泛用于航空发动机、燃气轮机等转子系统各零部件的联接, 由于螺栓预紧力在装配等过程中不可避免地存在一定损失, 造成联接刚度值有不确定性^[1].导致转子系统动力学特性理论分析与真实结果往往存在偏差.
目前, 在含螺栓联接结构转子动力学特性理论研究方面, Luan等^[2]在含螺栓法兰联接结构的航天器动力学分析中, 提出以分段线性刚度弹簧代替螺栓法兰联接结构的简化模型, 通过实验和理论分析, 验证了该模型在动力学分析中的可行性.Hu等^[3]基于Lagrange动力学方程建立了考虑非线性油膜力、不平衡质量及不平衡拉杆预紧力的拉杆紧固-轴承转子系统的非线性动力学模型, 采用四阶隆格库塔法求解方程, 分析该转子系统的动力学特性.Hei等^[4]建立了考虑陀螺效应的拉杆紧固-有限长轴承转子系统的动力学模型, 研究转子系统周期、倍周期、准周期等丰富的非线性特性.Qin等^[5-6]推导出盘毂联接结构时变弯曲刚度表达式, 并以简单盘毂螺栓联接转子为例, 分析了螺栓联接结构时变刚度对转子动力学特性的影响.以上研究均表明螺栓联接结构对转子系统动力学特性有重要影响, 但联接刚度参数不确定性对转子系统动力学特性的影响目前尚无明确量化结论.工程概率不确定性分析方法的推广, 使量化螺栓联接结构刚度不确定性对转子系统动力学特性影响的研究成为可能^[7].在工程概率不确定性分析领域, 王瑞利等^[8-9]采用多项式混沌展开法分析与辨识复杂工程模型中的误差与不确定度的来源, 并量化不确定度的影响.Li等^[10]基于非嵌入多项式混沌展开理论, 分析和量化了双转子系统中偏斜角、阻尼、非线性支承刚度等随机参数对转子系统动力学特性的影响.
不确定程度在工程问题分析中直接影响理论分析结果的准确程度, 量化随机参数的影响对含螺栓联接结构转子系统的设计有重要意义.本文基于有限元法构建螺栓法兰联接结构单元, 分析轴向联接刚度对转子动力学特性的影响, 并采用非嵌入混沌多项式展开法分析螺栓联接结构轴向刚度不确定性对转子系统动力学特性的影响.
1 含螺栓联接结构转子动力学模型为分析含螺栓联接结构转子系统中螺栓法兰轴向联接刚度不确定性对转子系统动力学特性的影响, 需先根据转子系统几何模型及相关参数构建转子系统整体有限元模型, 对其进行动力学特性分析, 进而完成轴向联接刚度参数不确定性的量化.
1.1 几何模型及参数说明为研究转子系统中螺栓法兰联接结构轴向刚度不确定性对转子系统振动特性的影响, 采用中间盘替换为螺栓法兰联接结构的Jeffcott转子模型, 两端转轴仅与相邻盘联接, 其三维模型如图 1所示.在转子左端点建立坐标系, z轴与转轴和盘轴线重合, y轴为垂直转子轴线竖直向上方向.转子左端轴长度为L1, 右端轴长度为L2, 两个完全相同的盘总厚度为d, 轴直径为d_r, 盘直径为D_r.
图 1(Fig. 1)

图 1 含螺栓法兰联接结构的转子模型Fig.1 Schamatic of the rotor system with bolted joint structure

图 1转子系统中螺栓法兰联接结构由于其结构特性, 在受弯矩过程中角刚度呈分段线性^{[2, 4]}.小变形情况下拉、压刚度值相同, 随着载荷增大, 角位移曲线呈非线性变化.螺栓法兰联接结构简化受力模型如图 2所示^[13].
图 2(Fig. 2)

图 2 螺栓联接结构受载示意图Fig.2 Schamatic of the bolted joint structure under radial bending moment

由图 2可得螺栓法兰联接结构角刚度计算式:

(1)

式中:M_x为轴端承受的弯矩; D为相对安装的两螺栓间距; k_C, k_T分别为螺栓法兰联接结构的拉、压刚度.
当法兰间相对角位较小时, 可近似认为螺栓法兰联接结构拉、压刚度相等, 此时角刚度表达式:

(2)

1.2 构建螺栓法兰联接单元为分析螺栓法兰联接结构刚度参数对转子响应的影响, 基于有限单元法理论, 将图 1螺栓法兰联接结构构建为一个新的单元, 包括两个节点沿x与y方向的平动以及绕x与y轴的转动共8个自由度, 单元广义坐标为

(3)

式中:x, y, θ_x, θ_y分别表示x, y方向位移、绕x轴转角及绕y轴转角; 下标1, 2分别对应联接单元两个节点, 本节内容中下标含义相同.
根据给定的单元自由度, 可得该单元动能表达式:

(4)

式中, Ωm, J_d和J_p分别表示盘的极转角、转速质量、直径转动惯量和极转动惯量.
该单元势能表达式:

(5)

式中:K_x为单元x方向位移刚度; K_y为单元y方向位移刚度; K_θx为单元绕x轴转动角刚度; K_θy为单元绕y轴转动角刚度.
由单元动能和势能表达式, 可通过拉格朗日能量法得出单元质量矩阵、刚度矩阵、陀螺矩阵:

(6)

(7)

(8)

式中:k为转子系统两盘间的径向刚度; k_θ为螺栓法兰联接结构受弯矩作用下的角刚度，k_θ的值可通过式(2)得出.本文假设螺栓法兰联接结构之间的轴向刚度为线性刚度, 即k_T取定值.
1.3 动力学方程及求解根据1.1节中所述几何模型及1.2节构建的螺栓法兰联接单元, 基于Rayleigh梁理论, 建立图 3所示转子有限元模型.其中, 盘左侧转轴划分为9个梁单元, 右侧转轴划分为10个梁单元, 将联接单元置于两段轴之间, 最终建立的含螺栓法兰联接结构的转子简化模型.
图 3(Fig. 3)

图 3 含螺栓联接转子简化模型Fig.3 Simplified model of the rotor with bolted joint

根据图 3建立转子系统整体运动微分方程:

(9)

式中:M为转子系统整体质量矩阵; C为转子系统整体阻尼矩阵; G为陀螺力矩阵; K为转子系统整体刚度矩阵; F为系统外力矩阵, 此算例中为盘偏心力.
转子系统中轴和盘的材料密度为7 800 kg/m³、弹性模量为210 GPa、泊松比为0.3, 取L1=450 mm, L2=500 mm, d=20 mm, d_r=20 mm, D_r=100 mm, 并给定转子系统两端x, y方向支承刚度为5×10⁷ N/m, 阻尼为500 (N·m^-1)/s, 瑞利梁阻尼参数α=5, β=3×10^-8.在图 3所示联接单元处添加2.5 kg·mm的不平衡质量, 采用Newmark法求解有限元方程, 得到联接单元轴向刚度分别为1×10⁶, 1×10⁷, 1×10⁸和1×10⁹ N/m时联接单元y方向位移随转速变化曲线如图 4所示, 对应的临界转速及盘y方向稳态响应最大值如表 1所示.
图 4(Fig. 4)

图 4 不同轴向刚度盘竖直方向响应Fig.4 Vertical response of the disk with different axial stiffness

表 1(Table 1)

表 1 不同刚度下转子临界转速及盘竖直方向最大响应值Table 1 Critical speeds and maximum response values at critical speeds with different stiffness 轴向刚度/(N·m-1) 1×106 1×107 1×108 1×109 临界转速/(r·min-1) 1 300 1 800 1 900 1 900 最大幅值/mm 43.43 11.04 11.58 9.042

表 1 不同刚度下转子临界转速及盘竖直方向最大响应值 Table 1 Critical speeds and maximum response values at critical speeds with different stiffness

由图 4和表 1可知, 螺栓法兰联接结构的轴向刚度对转子系统临界转速及临界转速状态下盘y方向稳态响应最大幅值有明显影响.图 5为不同轴向刚度对应的盘在临界转速工作状态下竖直方向稳态响应幅值曲线.
图 5(Fig. 5)

图 5 不同轴向刚度下盘竖直方向最大幅值曲线Fig.5 Response curves of the disk in vertical direction

图 5可进一步说明螺栓法兰联接结构轴向联接刚度在一定范围内刚度越小, 转子系统临界转速状态下盘稳态响应幅值越大, 但当轴向联接刚度达到一定值后, 该影响变小.当螺栓因装配不当导致预紧力分布不均、应力松弛过大等问题时, 螺栓法兰联接结构轴向联接刚度会因此降低, 这些因素往往无法准确预测和计算, 因此, 联接结构轴向刚度可视为不确定性因素, 针对该刚度值不确定性对转子系统振动特性影响的研究具有重要意义.
2 不确定参数动力学分析不确定度量化是根据统计学方法, 量化不确定参数的影响.目前, 工程中多采用Monte Carlo法和多项式混沌展开法分析量化参数不确定的影响^[8-9].本文采用多项式混沌展开法分析航空发动机转子螺栓联接结构轴向刚度不确定性对系统响应的影响^[10].
将不确定参数转换为随机高斯型参数, 可表示为

(10)

式中:f为螺栓预紧力均值; ζ为随机变量; δ为对应的标准差.
将系统响应展开为混沌多项式:

(11)

式中:X(τ), Y(τ)为展开多项式系数; ψ_k(ζ)为高斯艾尔米特展开多项式; p为高斯艾尔米特多项式混沌展开阶数.
多项式系数X(τ), Y(τ)可由系统真实响应X(τ, ζ_i)计算得出.最后, 由式(12)计算得出系统均值响应:

(12)

式中:ζ为高斯积分点, 与式(10)中的随机变量对应; ρ(ζ)为ζ的概率密度函数.
为分析螺栓预紧力不确定性对转子系统动力学特性的影响, 给定式(10)中螺栓轴向联接刚度为均值1×10⁶N/m, 标准差分别为5%, 10%, 15%及20%的高斯型随机数.将转子系统响应展开为一维9阶高斯艾尔米特多项式, 前6阶展开式基函数表达式及其内积如表 2所示.
表 2(Table 2)

表 2 一维展开式基函数及其内积Table 2 One-dimensional expansion base function and its inner product k ψk(ζ) E[ψk(ζ)2] 0 1 1 1 ζ 1 2 ζ2－1 2 3 ζ3－3ζ 6 4 ζ4－6ζ2+3 24 5 ζ(ζ4－10ζ2+15) 120

表 2 一维展开式基函数及其内积 Table 2 One-dimensional expansion base function and its inner product

由转子系统动力学模型及多项式混沌展开法理论, 可得转子均值临界转速曲线如图 6所示.
图 6(Fig. 6)

图 6 不同标准差对应的均值临界转速曲线Fig.6 Average critical speed curve with different standard deviations

由图 6可知, 螺栓法兰联接结构轴向刚度标准差较小时, 转子系统均值临界转速与预期临界转速基本吻合, 这是由于当标准差较小时, 刚度分布范围较小, 此范围内轴向联接刚度对转子系统临界转速及临界转速对应的幅值影响程度相近; 当标准差较大时, 刚度分布范围较大, 此时整个刚度值分布范围内的螺栓联接结构轴向刚度值对转子系统临界转速及临界转速对应的幅值影响程度不同, 且差别较大, 因此, 不同轴向联接刚度标准差对应的转子系统临界转速均值曲线差别明显.
由于不同联接刚度对应的振动频率不同, 因此在不确定刚度对系统响应均值影响的分析中, 难以用不同时间下对应的响应值做计算, 因此, 后续分析中, 以点数Δt作为响应曲线的横坐标.为具体分析、量化上述影响, 计算得出转速为1 300 r/min, 轴向联接刚度标准偏差分别为5%, 10%, 15%及20%时盘竖直方向的均值响应如图 7所示.
图 7(Fig. 7)

图 7 转速1 300 r/min时盘竖直方向均值响应Fig.7 Average response of the disk in the vertical direction under 1 300 r/min

由图 7可知, 预期临界转速工作状态下, 盘竖直方向均值响应随着标准差增大而降低, 这是由于标准差增大, 螺栓联接结构轴向联接刚度分布范围变大, 而较小的轴向联接刚度对系统响应的影响程度明显大于轴向联接刚度较大时的情况.为进一步分析螺栓法兰轴向联接刚度不确定性对转子系统低阶临界转速工作状态下稳态响应的影响, 根据图 6所得4种轴向联接刚度标准差下所得转子最大均值响应对应的临界转速, 分别计算盘竖直方向均值响应如图 8所示.
图 8(Fig. 8)

图 8 均值临界转速下盘竖直方向均值响应Fig.8 Average response of the disk in the vertical direction under the mean critical speeds

图 7、图 8对比可知, 螺栓法兰轴向联接刚度标准差为5%时, 对盘竖直方向稳态响应影响较小, 当螺栓法兰轴向联接刚度标准差为10%, 15%和20%时, 由于刚度分布范围较大, 该影响显著.
3 结论1) 在一定变化范围内, 螺栓法兰轴向联接刚度减小, 转子系统临界转速降低且对应的盘稳态响应幅值增大, 当轴向联接刚度增加到一定值后, 联接结构近似等效于整体结构, 此时该影响变小.
2) 螺栓法兰联接结构轴向刚度标准差较小时, 由于刚度分布范围较小, 转子系统均值临界转速与预期临界转速基本吻合.当标准差较大时, 刚度分布范围较大, 转子系统均值临界转速与预期临界转速差别明显.
3) 预期临界转速工作状态下, 标准差较大时, 由于螺栓联接结构轴向联接刚度分布范围变大, 盘竖直方向均值响应随着标准差增大而降低, 标准差较小时, 情况相反.
参考文献

[1]	Abid M, Khan A, Nash D H, et al. Optimized bolt tightening strategies for gasketed flanged pipe joints of different sizes[J].International Journal of Pressure Vessels & Piping, 2016, 139/140: 22–27.
[2]	Luan Y, Guan Z, Cheng G, et al. A simplified nonlinear dynamic model for the analysis of pipe structures with bolted flange joints[J].Journal of Sound and Vibration, 2012, 331(2): 325–344.DOI:10.1016/j.jsv.2011.09.002
[3]	Hu L, Liu Y, Zhao L, et al. Nonlinear dynamic behaviors of circumferential rod fastening rotor under unbalanced pre-tightening force[J].Archive of Applied Mechanics, 2016, 86(9): 1–11.
[4]	Hei D, Lu Y, Zhang Y, et al. Nonlinear dynamic behaviors of rod fastening rotor-hydrodynamic journal bearing system[J].Archive of Applied Mechanics, 2015, 85(7): 1–21.
[5]	Qin Z Y, Han Q K, Chu F L. Bolt loosening at rotating joint interface and its influence on rotor dynamics[J].Engineering Failure Analysis, 2016, 59: 456–466.
[6]	Qin Z Y, Han Q K, Chu F L. Analytical model of bolted disk-drum joints and its application to dynamic analysis of jointed rotor[J].Proceedings of the Institution of Mechanical Engineers, Part C:Journal of Mechanical Engineering Science, 2014, 228(4): 646–663.DOI:10.1177/0954406213489084
[7]	Hosder S, Grossman B, Haftka R, et al. Quantitative relative comparison of CFD simulation uncertainties for a transonic diffuser problem[J].Computers & Fluids, 2006, 35(10): 1444–1458.
[8]	王瑞利, 梁霄. 多项式混沌方法在偶然不确定度量化中的应用[J].数学建模及其应用, 2016, 5(3): 17–24, 81. ( Wang Rui-li, Liang Xiao. Application of polynomial chaos method in quantifying the aleatory uncertainty[J].Mathematical Modeling and Its Applications, 2016, 5(3): 17–24, 81.DOI:10.3969/j.issn.2095-3070.2016.03.007)
[9]	刘全, 王瑞利, 林忠. 非嵌入式多项式混沌方法在拉氏计算中的应用[J].固体力学学报, 2013, 33.
[10]	Liu Quan, Wang Rui-li, Lin Zhong. Uncertainty quantification for Lagrangian computation using non-intrusive polynomial chaos[J].Chinese Journal of Solid Mechanics, 2013, 33(sup1): 224–233.
[11]	刘卓乾, 曹树谦, 郭虎伦, 等. 含螺栓联接转子系统非线性振动特性研究[J].振动与冲击, 2016, 35(22): 10–16, 37. ( Liu Zhuo-qian, Cao Shu-qian, Guo Hu-lun, et al. Vibration characteristics of rotor systems with bolt joints[J].Journal of Vibration and Shock, 2016, 35(22): 10–16, 37.)