连接局部迟滞非线性的时频域动力学降阶方法

王东 张周锁

摘要： 针对连接结构时域和频域的非线性振动问题，提出一种基于局部非线性转化的动力学降阶方法。采用Newmark法和多谐波平衡法分别将时域和频域的非线性动力学方程转化为非线性代数方程组，将整体结构的非线性动力学响应降阶到仅与连接非线性相关的自由度上进行求解，通过减小迭代过程中Jacobian矩阵的维数来提高计算效率。采用Iwan模型描述连接界面的局部迟滞非线性，利用三自由度质量弹簧振子和连接梁结构研究非线性动力学分析方法的计算精度和效率。结果表明，本文建立的降阶方法预测的非线性动力学响应与现有的未降价方法吻合较好，提出的非线性动力学降阶方法能够有效地减少迭代过程的计算耗费，提高计算效率，时域方法约提高30%，频域方法提高近70倍。

关键词： 迟滞非线性; 连接界面; 动力学降阶; 多谐波平衡法; Newmark法

中图分类号： O322    文献标志码： A    文章编号： 1004-4523（2021）03-0559-08

DOI：10.16385/j.cnki.issn.1004-4523.2021.03.013

引  言

连接是复杂机械装配结构中不可或缺的重要组成单元，广泛地应用于部组件之间载荷和能量传递[1]。连接界面的存在引起复杂的动力学特征，也是装配结构中的薄弱环节，严重地影响机械设备的可靠性、稳定性和实用性能。动力学响应的预测对机械结构的设计、优化、控制、健康监测有重要的作用[2?3]。

非线性动力学响应的预测方法主要分为时域和频域两大类方法[4]。时域方法又分为直接数值积分和半解析方法。直接时程积分方法包括中心差分法、Runge?Kutta法、Newmark法、Wilson?θ法等[5]。采用时程积分方法求解非线性动力学问题需要耗费巨大的计算资源[6]。Oldfield等采用Iwan、Bouc?Wen模型描述单螺栓连接结构的迟滞非线性行为，利用四阶Runge?Kutta法对非线性动力学方程进行求解[7]。Song等[8?9]和Gaul等[10]也采用时程积分法对螺栓连接结构的瞬态动力学响应进行求解。半解析法可以直接获得结构的动力学响应，将非线性系统截断为多个考虑不同初值的线性系统，但需要在每个时程截断点求解过渡方程，仅适用于特殊的非线性系统，如分段线性化模型[11]。在求解小阻尼结构稳态响应时，计算瞬态响应仍需要花费大量的计算资源，尤其是高频激励载荷。因此，研究者倾向于在频域求解非线性动力学方程[12]。

针对非线性动力学方程的稳态响应，研究者提出了多种半解析的频域近似方法，如摄动法、平均法、渐近法、多尺度法、谐波平衡法、增量谐波平衡法等[4，13?16]。其中，谐波平衡法求解过程简单，已经广泛地被用于求解非线性系统的周期稳态响应。谐波平衡法将动力学响应展开成一系列含未知系数的傅里叶级数，将非线性偏微分方程转化为非线性代数方程，通过匹配方程的谐波系数获得系统的稳态非线性响应[10，17?18]。利用谐波平衡法求解非线性动力学问题时，一般将研究对象简化为单自由度或低自由度的等效模型，截取的谐波阶数对计算精度有重要的影响，以上的研究在求解过程中普遍保留一阶或较少阶的谐波项。随着谐波阶次的增加，计算精度变高但计算量迅速增大，计算效率降低。

对非线性动力学方程的非线性项进行傅里叶变换时，可采用混合时频域方法对结果进行求解。Cameron?Griffin首先将该方法应用于求解含迟滞非线性的单自由度系统，基于时频交替变换和谐波平衡法求解稳态非线性动力学响应。这种方法融合了频域求解振动方程的高效性和时域判断非线性力的便捷性，通过快速离散傅里叶正?逆变换，反复迭代获得系统的稳态非线性动力学响应[6，18?22]。Petrov等开发的谐波平衡法求解器FORSE （Force Response Suite）已被广泛地应用于摩擦系统的非线性动力学响应求解，如叶盘、发动机系统[23?26]。Krack等研发的开源程序NLVIB （Nonlinear Vibration）也可以用来求解非线性方程组[27?28]。两种求解器均采用Newton's迭代法搜尋匹配的谐波系数[29?30]。由于Jacobian矩阵的病态（奇异矩阵），求解非线性响应的迭代过程往往是条件收敛的[31?33]。在迭代过程中，Jacobian矩阵的构造与求逆（或特征值求解）耗费了大量的计算资源。因此，发展有效的非线性动力学降阶方法，减小Jacobian矩阵的规模，提高迭代过程的计算效率是非常必要的。

本文提出一种基于局部非线性转化的动力学降阶方法。针对连接结构时域和频域非线性振动问题，分别采用Newmark法和多谐波平衡法将非线性动力学方程转化为非线性代数方程组，同时将整体结构非线性动力学响应的求解转化到仅与连接非线性相关的自由度上进行求解，减小迭代过程的计算规模，提高计算效率。采用Iwan模型描述连接界面的局部迟滞非线性力学行为，利用三自由度质量弹簧振子和连接梁结构验证本文的方法。

1 时域非线性动力学降阶算法

1.1 非线性动力学方程代数化

考虑连接界面上局部非线性行为的动力学微分方程为

式中  分别为质量、阻尼、刚度矩阵;为动力学响应;为局部非线性恢复力;为非线性模型的参数;为外激励力。

采用Newmark法将非线性动力学微分方程转化为非线性代数方程组。采用逐步积分法建立由时刻到时刻的状态向量的递推关系。时刻的未知量，和满足

式（9）中，时刻的位移响应与非线性恢复力是耦合的，需要进行迭代求解。采用Newton's迭代方法对子步的非线性响应进行求解，定义残差方程为

1.2 动力学降阶方法

在机械结构中，由连接引起的非线性往往是局部的，与非线性相关的自由度数目远远小于整体结构的自由度。利用一位置转换矩阵提取局部非线性动力学响应。

利用局部坐标中非线性恢复力构造原物理坐标系下的非线性恢复力。

式中  为连接界面局部迟滞非线性恢复力。

忽略非线性自由度之间耦合特性，局部非线性动力学响应的迭代格式定义为

式中  为连接非线性相关的自由度序号。

由于Jacobian矩阵的病态奇异性，非线性迭代过程往往是条件收敛的。本文采用松弛迭代原理提高迭代的收敛性能，式（16）中局部非线性动力学响应不完全更新。

式中  为松弛因子，介于0?1之间，越小，非线性响应更新的速率越慢，但收敛性越好。

将式（16）和（17）中迭代收敛的解代回式（9）和（10）可以获得整体结构的非线性动力学响应。

1.3 算  例

如图1所示，利用三自由度质量弹簧振子系统验证本文的时域非线性动力学降阶方法。利用Iwan模型描述连接界面非线性恢复力

如图2所示，基于Masing映射准则[34?36]，周期性激励载荷下Iwan模型的迟滞非线性恢复力为

基于文献[37]中螺栓连接结构的实验结果进行仿真，参数选取为：m1=5.28 kg，m2=0.55 kg，m3=5.21 kg，k1=1.09×107 N/m，k2=1.9×107 N/m，c1=0，c2=200 N?s/m， fext=Asin（2πft）。外激励载荷幅值，激励频率，计算瞬态动力学响应如图3所示，各激励频率稳态非线性动力学响应如图4所示。

由图3和4可知，本文方法预测的非线性动力学响应与降阶前的结果吻合较好，验证了本文方法的有效性。共振峰附近频响曲线变平是由迟滞非线性引起的能量耗散造成的。降阶前与降阶之后计算耗费对比结果如表1所示，数值仿真均在4?core Intel（R） Core（TM） i5?3470 3.2 GHz CPU上开展。结果表明，降阶之后，每个子步迭代矩阵的维数有所降低，非线性动力学响应的计算耗时也大大减少，计算效率约提高30%。

2 频域非线性动力学降阶算法

2.1 非线性动力学方程代数化

利用多谐波平衡法将连接界面局部非线性恢复力、外激励载荷与非线性响应进行谐波级数展开：

式中  为谐波的阶数;为外载荷的激励频率;为稳态非线性响应的谐波系数;和分别为激励载荷和非线性恢复力的谐波系数;为取实部算子;i=。

将式（23）?（25）代入式（1），非线性动力学微分方程转化为非线性代数方程组

式中  为动刚度的逆矩阵或传递函数，定义为

Newton's迭代的残差函数为

整体结构非线性动力学响应的迭代格式定义为

式（29）中，每个连接非线性自由度上动力学响应的待求系数为项，包括零阶常数项、每阶谐波系数的余弦项系数和正弦项系数。

2.2 动力学降阶方法

与时域动力学降阶方法相似，连接界面上局部非线性动力学响应定义为

Newton's迭代的殘差函数重新定义为

上式仅考虑与连续非线性相关的自由度上的响应与外激励和非线性恢复力之间的传递关系。

忽略各非线性自由度之间的耦合，局部非线性动力学响应的迭代格式为

不同于式（16）和（32），迭代矩阵的构造与非线性恢复力和响应的谐波系数之间的导数相关，难以直接写出解析的 ，本文采用普适的中心差分方法计算迭代矩阵的每一列[12]

式（33）表示局部非线性动力学响应的每一阶谐波系数的微元变化引起的非线性恢复力的微元变化。

2.3 算  例

如图5所示，利用一连接梁结构验证本文的频域非线性动力学降阶方法。有限元模型含有8个线性Euler梁单元，1个非线性连接单元。两个方向的Iwan模型用来描述连接界面上局部非线性恢复力。每个节点含有横向平动u、垂向平动v、转动θ三个自由度，共计30个自由度。线性梁单元的刚度矩阵、质量矩阵、非线性连接单元的恢复力表达式可以参考文献[9，12，38?39]。

局部非线性动力学响应的转换矩阵为

式中  为连接梁单元的长度;为宽度。

图5中，为残余刚度系数;两个Iwan模型的临界滑移力，黏着刚度分别定义为

式中  为梁单元的弹性模量;为截面转动惯量。

连接梁结构的材料为钢，密度ρ=7.85103 kg/m3，弹性模量为。激励幅值为，激励频率范围根据第一阶和第三阶弯曲模态频率进行选取，分别为和，谐波阶数取10（H=10），计算结果如图6所示。

由图6可知，本文方法预测的非线性动力学响应与降阶前的结果吻合较好，频响函数能够较好地反映非线性力学行为的影响，尤其是在共振峰附近的区域。由表2可知，降阶之后，Jacobian矩阵的维数大大降低，非线性动力学响应的计算耗时也大大减少，计算效率约提高70倍。

3 结  论

本文提出一种基于局部非线性转化的降阶非线性动力学算法。采用Newmark和多谐波平衡法将非线性动力学微分方程转化为非线性代数方程组，将整体结构的动力学响应转化到仅与连接非线性相关的自由度上进行求解，减小迭代过程的计算耗费。利用含有Iwan模型的三自由度质量弹簧振子和连接梁结构验证本文的非线性动力学降阶方法。

本文方法预测的非线性动力学响应与未降阶的Newmark和多谐波平衡法的结果吻合较好。频响函数能够较好地反映迟滞非线性引起的能量耗散的影响，尤其是在共振峰附近。降阶之后，非线性迭代过程的计算效率明显提高，时域降阶提高约30%，而频域降阶提高近70倍。

参考文献：

[1] Segalman D J， Gregory D L， Starr M J， et al. Handbook on dynamics of jointed structures[R]. Technical Report， Sandia National Laboratories， NM， USA， 2009.

[2] Brake M R W. The Mechanics of Jointed Structures： Recent Research and Open Challenges for Developing Predictive Models for Structural Dynamics[M]. Springer， 2018.

[3] 李玲， 蔡安江， 阮晓光. 栓接结合部非线性等效线性化方法[J]. 振动工程学报， 2015， 28（4）： 560-566.

Li L， Cai A J， Ruan X G. Equivalent linear method for hysteresis nonlinear of bolted joints[J]. Journal of Vibration Engineering， 2015， 28（4）： 560-566.

[4] Nayfeh A H， Mook D T. Nonlinear Oscillations[M].  Wiley， 1995： 49-65.

[5] 張相盟. 摩擦连接的结构非线性动力学研究[D]. 哈尔滨：哈尔滨工业大学， 2013.

Zhang X M. Research on structural nonlinear dynamics of frictional joints?[D]. Harbin： Harbin Institute Technology， 2013.

[6] Zucca S， Firrone C M. Nonlinear dynamics of mechanical systems with friction contacts： Coupled static and dynamic multi-harmonic balance method and multiple solutions[J]. Journal of Sound and Vibration， 2014， 333（3）： 916-926.

[7] Oldfield M， Ouyang H， Mottershead J E. Simplified models of bolted joints under harmonic loading?[J]. Computers and Structures， 2005， 84（1）： 25-33.

[8] Song Y， Hartwigsen C J， Bergman L A， et al. A three-dimensional nonlinear reduced-order predictive joint model[J]. Earthquake Engineering and Engineering Vibration， 2003， 2（1）： 59-73.

[9] Song Y， Hartwigsen C J， Mcfarland D M. Simulation of dynamics of beam structures with bolted joints using adjusted Iwan beam elements[J]. Journal of Sound and Vibration， 2004， 273（1）： 249-276.

[10] Gaul L， Lenz J. Nonlinear dynamics of structures assembled by bolted joints[J]. Acta Mechanica， 1997， 125（1-4）： 169-181.

[11] Matthew S， Allen H S， Epp D S. Piecewise-linear restoring force surfaces for semi-nonparametric identification of nonlinear systems?[J]. Nonlinear Dynamics， 2008， 54（1-2）： 123-135.

[12] Wang D， Zhang Z. High-efficiency nonlinear dynamic analysis for joint interfaces with Newton-Raphson iteration process[J]. Nonlinear Dynamics， 2020，100（1）： 1-17.

[13] 刘南， 白俊强， 华俊，等. 高阶谐波平衡方法中非物理解来源分析及改进方法研究[J]. 力学学报， 2016， 48（4）： 897-906.

Liu N， Bai J Q， Hua J， et al. Investigation of the source and improvement of non-physical solutions in high-order harmonic balance?[J]. Chinese Journal of Theoretical and Applied Mechanics，2016，48（4）?：897-906.

[14] David J W， Mitchell L D， Daws J W. Using transfer matrices for parametric system forced response[J]. Journal of Vibration and Acoustics， 1987，109（4）?： 356-360.

[15] Lacayo R， Pesaresi L， Gro? J， et al. Nonlinear modeling of structures with bolted joints： A comparison of two approaches based on a time-domain and frequency-domain solver[J]. Mechanical Systems and Signal Processing， 2019， 114： 413-438.

[16] Zhou B， Thouverez F， Lenoir D. A variable-coefficient harmonic balance method for the prediction of quasi-periodic response in nonlinear systems[J]. Mechanical Systems and Signal Processing， 2015， 64-65： 233-244.

[17] Ren Y. Identification of properties of nonlinear joints using dynamic test data?[J]. Journal of Vibration and Acoustics， 1998， 120（2）： 324-330.

[18] Groll G V， Ewins D J. The harmonic balance method with arc-length continuation in rotor/stator contact problems[J]. Journal of Sound and Vibration， 2001， 241（2）： 223-233.

[19] Cameron T M， Griffin J H. An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems[J]. Journal of Applied Mechanics， 1989， 56（1）： 149-154.

[20] Ling F H. Discussion： “An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems” （Cameron， T. M.， and Griffin， J. H.， 1989， ASME J. Appl Mech.， 56，PP.149-154）[J]. Journal of Applied Mechanics， 1990， 57（1）： 251.

[21] Zhang Z Y， Chen Y S. Harmonic balance method with alternating frequency/time domain technique for nonlinear dynamical system with fractional exponential[J]. Applied Mathematics and Mechanics， 2014， 35（4）： 423-436.

[22] Liu Y， Bo S， Xu Z. Improved hybrid frequency-time domain method for nonlinear analysis of frictionally damped blade systems[C]. Turbomachinery Technical Conference and Exposition， Seoul， South Korea， 2016.

[23] Petrov E P. A high-accuracy model reduction for analysis of nonlinear vibrations in structures with contact interfaces[J]. Journal of Engineering for Gas Turbines and Power， 2011， 133（10）： 102503.

[24] Sever I A， Petrov E P， Ewins D J. Experimental and numerical investigation of rotating bladed disk forced response using underplatform friction dampers[J]. Journal of Engineering for Gas Turbines and Power， 2008， 130（4）： 042503.

[25] Petrov E P， Ewins D J. Analytical formulation of friction interface elements for analysis of nonlinear multi-harmonic vibrations of bladed discs[J]. Journal of Turbomachinery， 2003， 125（2）： 364-371.

[26] Petrov E P， Ewins D J. State-of-the-art dynamic analysis for non-linear gas turbine structures[J]. Proceedings of the Institution of Mechanical Engineers， Part G： Journal of Aerospace Engineering， 2004， 218（3）： 199-211.

[27] Krack M， Gross J. Harmonic Balance for Nonlinear Vibration Problems[M].  Springer， 2019.

[28] Krack M， Salles L， Thouverez F. Vibration prediction of bladed disks coupled by friction joints[J]. Archives of Computational Methods in Engineering， 2017， 24（3）： 589-636.

[29] Amata S， Argyrosb I K， Busquiera S. Newton-type methods on Riemannian manifolds under Kantorovich-type conditions[J]. Applied Mathematics and Computation， 2014， 227（1）： 762-787.

[30] Argyros I K. On Newton's method under mild differentiability conditions and applications[J]. Applied Mathematics and Computation， 1999， 102（2-3）： 177-183.

[31] Afzal M， Arteaga I L， Kari L. An analytical calculation of the Jacobian matrix for 3D friction contact model applied to turbine blade shroud contact[J]. Computers and Structures， 2016， 177（1）： 204-217.

[32] Sombroek C S M， Tiso P， Renson L， et al. Numerical computation of nonlinear normal modes in a modal derivative subspace[J]. Computers and Structures， 2018， 195（1）： 34-46.

[33] Peletan L， Baguet S， Torkhani M， et al. A comparison of stability computational methods for periodic solution of nonlinear problems with application to rotordynamics[J]. Nonlinear Dynamics， 2013， 72（3）： 671-682.

[34] Sü? D， Willner K. Investigation of a jointed friction oscillator using the multiharmonic balance method[J]. Mechanical Systems and Signal Processing， 2015， 52-53（1）： 73-87.

[35] Rajaei M， Ahmadian H. Development of generalized Iwan model to simulate frictional contacts with variable normal loads?[J]. Applied Mathematical Modelling， 2014， 38（15-16）： 4006-4018.

[36] Argatov I I， Butcher E A. On the Iwan models for lap-type bolted joints[J]. International Journal of Non-Linear Mechanics， 2011， 46（2）： 347-356.

[37] Segalman D J， Starr M J. Inversion of Masing models via continuous Iwan systems[J]. International Journal of Non-Linear Mechanics， 2008， 43（1）： 74-80.

[38] Yuan P P， Ren W X， Zhang J. Dynamic tests and model updating of nonlinear beam structures with bolted joints[J]. Mechanical Systems and Signal Processing， 2019， 126（1）： 193-210.

[39] Ahmadian H， Rajaei M. Identification of Iwan distribution density function in frictional contacts[J]. Journal of Sound and Vibration， 2014， 333（15）： 3382-3393.

作者简介： 王  东（1988-），男，助理研究员。电话：（0816）2485436; E-mail：king_east@sina.cn

通讯作者： 张周锁（1963-），男，教授。电话：（029）82663689; E-mail：zzs@mail.xjtu.edu.cn