連接局部遲滯非線性的時頻域動力學降階方法

王東 張周鎖

摘要： 針對連接結構時域和頻域的非線性振動問題，提出一種基于局部非線性轉化的動力學降階方法。采用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