王東 張周鎖

摘要： 針對連接結(jié)構(gòu)時域和頻域的非線性振動問題，提出一種基于局部非線性轉(zhuǎn)化的動力學(xué)降階方法。采用Newmark法和多諧波平衡法分別將時域和頻域的非線性動力學(xué)方程轉(zhuǎn)化為非線性代數(shù)方程組，將整體結(jié)構(gòu)的非線性動力學(xué)響應(yīng)降階到僅與連接非線性相關(guān)的自由度上進行求解，通過減小迭代過程中Jacobian矩陣的維數(shù)來提高計算效率。采用Iwan模型描述連接界面的局部遲滯非線性，利用三自由度質(zhì)量彈簧振子和連接梁結(jié)構(gòu)研究非線性動力學(xué)分析方法的計算精度和效率。結(jié)果表明，本文建立的降階方法預(yù)測的非線性動力學(xué)響應(yīng)與現(xiàn)有的未降價方法吻合較好，提出的非線性動力學(xué)降階方法能夠有效地減少迭代過程的計算耗費，提高計算效率，時域方法約提高30%，頻域方法提高近70倍。

關(guān)鍵詞： 遲滯非線性; 連接界面; 動力學(xué)降階; 多諧波平衡法; Newmark法

中圖分類號： O322? ? 文獻標志碼： A? ? 文章編號： 1004-4523（2021）03-0559-08

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

引? 言

連接是復(fù)雜機械裝配結(jié)構(gòu)中不可或缺的重要組成單元，廣泛地應(yīng)用于部組件之間載荷和能量傳遞[1]。連接界面的存在引起復(fù)雜的動力學(xué)特征，也是裝配結(jié)構(gòu)中的薄弱環(huán)節(jié)，嚴重地影響機械設(shè)備的可靠性、穩(wěn)定性和實用性能。動力學(xué)響應(yīng)的預(yù)測對機械結(jié)構(gòu)的設(shè)計、優(yōu)化、控制、健康監(jiān)測有重要的作用[2?3]。

非線性動力學(xué)響應(yīng)的預(yù)測方法主要分為時域和頻域兩大類方法[4]。時域方法又分為直接數(shù)值積分和半解析方法。直接時程積分方法包括中心差分法、Runge?Kutta法、Newmark法、Wilson?θ法等[5]。采用時程積分方法求解非線性動力學(xué)問題需要耗費巨大的計算資源[6]。Oldfield等采用Iwan、Bouc?Wen模型描述單螺栓連接結(jié)構(gòu)的遲滯非線性行為，利用四階Runge?Kutta法對非線性動力學(xué)方程進行求解[7]。Song等[8?9]和Gaul等[10]也采用時程積分法對螺栓連接結(jié)構(gòu)的瞬態(tài)動力學(xué)響應(yīng)進行求解。半解析法可以直接獲得結(jié)構(gòu)的動力學(xué)響應(yīng)，將非線性系統(tǒng)截斷為多個考慮不同初值的線性系統(tǒng)，但需要在每個時程截斷點求解過渡方程，僅適用于特殊的非線性系統(tǒng)，如分段線性化模型[11]。在求解小阻尼結(jié)構(gòu)穩(wěn)態(tài)響應(yīng)時，計算瞬態(tài)響應(yīng)仍需要花費大量的計算資源，尤其是高頻激勵載荷。因此，研究者傾向于在頻域求解非線性動力學(xué)方程[12]。

針對非線性動力學(xué)方程的穩(wěn)態(tài)響應(yīng)，研究者提出了多種半解析的頻域近似方法，如攝動法、平均法、漸近法、多尺度法、諧波平衡法、增量諧波平衡法等[4，13?16]。其中，諧波平衡法求解過程簡單，已經(jīng)廣泛地被用于求解非線性系統(tǒng)的周期穩(wěn)態(tài)響應(yīng)。諧波平衡法將動力學(xué)響應(yīng)展開成一系列含未知系數(shù)的傅里葉級數(shù)，將非線性偏微分方程轉(zhuǎn)化為非線性代數(shù)方程，通過匹配方程的諧波系數(shù)獲得系統(tǒng)的穩(wěn)態(tài)非線性響應(yīng)[10，17?18]。利用諧波平衡法求解非線性動力學(xué)問題時，一般將研究對象簡化為單自由度或低自由度的等效模型，截取的諧波階數(shù)對計算精度有重要的影響，以上的研究在求解過程中普遍保留一階或較少階的諧波項。隨著諧波階次的增加，計算精度變高但計算量迅速增大，計算效率降低。

對非線性動力學(xué)方程的非線性項進行傅里葉變換時，可采用混合時頻域方法對結(jié)果進行求解。Cameron?Griffin首先將該方法應(yīng)用于求解含遲滯非線性的單自由度系統(tǒng)，基于時頻交替變換和諧波平衡法求解穩(wěn)態(tài)非線性動力學(xué)響應(yīng)。這種方法融合了頻域求解振動方程的高效性和時域判斷非線性力的便捷性，通過快速離散傅里葉正?逆變換，反復(fù)迭代獲得系統(tǒng)的穩(wěn)態(tài)非線性動力學(xué)響應(yīng)[6，18?22]。Petrov等開發(fā)的諧波平衡法求解器FORSE （Force Response Suite）已被廣泛地應(yīng)用于摩擦系統(tǒng)的非線性動力學(xué)響應(yīng)求解，如葉盤、發(fā)動機系統(tǒng)[23?26]。Krack等研發(fā)的開源程序NLVIB （Nonlinear Vibration）也可以用來求解非線性方程組[27?28]。兩種求解器均采用Newton's迭代法搜尋匹配的諧波系數(shù)[29?30]。由于Jacobian矩陣的病態(tài)（奇異矩陣），求解非線性響應(yīng)的迭代過程往往是條件收斂的[31?33]。在迭代過程中，Jacobian矩陣的構(gòu)造與求逆（或特征值求解）耗費了大量的計算資源。因此，發(fā)展有效的非線性動力學(xué)降階方法，減小Jacobian矩陣的規(guī)模，提高迭代過程的計算效率是非常必要的。

本文提出一種基于局部非線性轉(zhuǎn)化的動力學(xué)降階方法。針對連接結(jié)構(gòu)時域和頻域非線性振動問題，分別采用Newmark法和多諧波平衡法將非線性動力學(xué)方程轉(zhuǎn)化為非線性代數(shù)方程組，同時將整體結(jié)構(gòu)非線性動力學(xué)響應(yīng)的求解轉(zhuǎn)化到僅與連接非線性相關(guān)的自由度上進行求解，減小迭代過程的計算規(guī)模，提高計算效率。采用Iwan模型描述連接界面的局部遲滯非線性力學(xué)行為，利用三自由度質(zhì)量彈簧振子和連接梁結(jié)構(gòu)驗證本文的方法。

1 時域非線性動力學(xué)降階算法

1.1 非線性動力學(xué)方程代數(shù)化

考慮連接界面上局部非線性行為的動力學(xué)微分方程為

式中? 分別為質(zhì)量、阻尼、剛度矩陣;為動力學(xué)響應(yīng);為局部非線性恢復(fù)力;為非線性模型的參數(shù);為外激勵力。

采用Newmark法將非線性動力學(xué)微分方程轉(zhuǎn)化為非線性代數(shù)方程組。采用逐步積分法建立由時刻到時刻的狀態(tài)向量的遞推關(guān)系。時刻的未知量，和滿足

式（9）中，時刻的位移響應(yīng)與非線性恢復(fù)力是耦合的，需要進行迭代求解。采用Newton's迭代方法對子步的非線性響應(yīng)進行求解，定義殘差方程為

1.2 動力學(xué)降階方法

在機械結(jié)構(gòu)中，由連接引起的非線性往往是局部的，與非線性相關(guān)的自由度數(shù)目遠遠小于整體結(jié)構(gòu)的自由度。利用一位置轉(zhuǎn)換矩陣提取局部非線性動力學(xué)響應(yīng)。

利用局部坐標中非線性恢復(fù)力構(gòu)造原物理坐標系下的非線性恢復(fù)力。

式中? 為連接界面局部遲滯非線性恢復(fù)力。

忽略非線性自由度之間耦合特性，局部非線性動力學(xué)響應(yīng)的迭代格式定義為

式中? 為連接非線性相關(guān)的自由度序號。

由于Jacobian矩陣的病態(tài)奇異性，非線性迭代過程往往是條件收斂的。本文采用松弛迭代原理提高迭代的收斂性能，式（16）中局部非線性動力學(xué)響應(yīng)不完全更新。

式中? 為松弛因子，介于0?1之間，越小，非線性響應(yīng)更新的速率越慢，但收斂性越好。

將式（16）和（17）中迭代收斂的解代回式（9）和（10）可以獲得整體結(jié)構(gòu)的非線性動力學(xué)響應(yīng)。

1.3 算? 例

如圖1所示，利用三自由度質(zhì)量彈簧振子系統(tǒng)驗證本文的時域非線性動力學(xué)降階方法。利用Iwan模型描述連接界面非線性恢復(fù)力

如圖2所示，基于Masing映射準則[34?36]，周期性激勵載荷下Iwan模型的遲滯非線性恢復(fù)力為

基于文獻[37]中螺栓連接結(jié)構(gòu)的實驗結(jié)果進行仿真，參數(shù)選取為：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）。外激勵載荷幅值，激勵頻率，計算瞬態(tài)動力學(xué)響應(yīng)如圖3所示，各激勵頻率穩(wěn)態(tài)非線性動力學(xué)響應(yīng)如圖4所示。

由圖3和4可知，本文方法預(yù)測的非線性動力學(xué)響應(yīng)與降階前的結(jié)果吻合較好，驗證了本文方法的有效性。共振峰附近頻響曲線變平是由遲滯非線性引起的能量耗散造成的。降階前與降階之后計算耗費對比結(jié)果如表1所示，數(shù)值仿真均在4?core Intel（R） Core（TM） i5?3470 3.2 GHz CPU上開展。結(jié)果表明，降階之后，每個子步迭代矩陣的維數(shù)有所降低，非線性動力學(xué)響應(yīng)的計算耗時也大大減少，計算效率約提高30%。

2 頻域非線性動力學(xué)降階算法

2.1 非線性動力學(xué)方程代數(shù)化

利用多諧波平衡法將連接界面局部非線性恢復(fù)力、外激勵載荷與非線性響應(yīng)進行諧波級數(shù)展開：

式中? 為諧波的階數(shù);為外載荷的激勵頻率;為穩(wěn)態(tài)非線性響應(yīng)的諧波系數(shù);和分別為激勵載荷和非線性恢復(fù)力的諧波系數(shù);為取實部算子;i=。

將式（23）?（25）代入式（1），非線性動力學(xué)微分方程轉(zhuǎn)化為非線性代數(shù)方程組

式中? 為動剛度的逆矩陣或傳遞函數(shù)，定義為

Newton's迭代的殘差函數(shù)為

整體結(jié)構(gòu)非線性動力學(xué)響應(yīng)的迭代格式定義為

式（29）中，每個連接非線性自由度上動力學(xué)響應(yīng)的待求系數(shù)為項，包括零階常數(shù)項、每階諧波系數(shù)的余弦項系數(shù)和正弦項系數(shù)。

2.2 動力學(xué)降階方法

與時域動力學(xué)降階方法相似，連接界面上局部非線性動力學(xué)響應(yīng)定義為

Newton's迭代的殘差函數(shù)重新定義為

上式僅考慮與連續(xù)非線性相關(guān)的自由度上的響應(yīng)與外激勵和非線性恢復(fù)力之間的傳遞關(guān)系。

忽略各非線性自由度之間的耦合，局部非線性動力學(xué)響應(yīng)的迭代格式為

不同于式（16）和（32），迭代矩陣的構(gòu)造與非線性恢復(fù)力和響應(yīng)的諧波系數(shù)之間的導(dǎo)數(shù)相關(guān)，難以直接寫出解析的 ，本文采用普適的中心差分方法計算迭代矩陣的每一列[12]

式（33）表示局部非線性動力學(xué)響應(yīng)的每一階諧波系數(shù)的微元變化引起的非線性恢復(fù)力的微元變化。

2.3 算? 例

如圖5所示，利用一連接梁結(jié)構(gòu)驗證本文的頻域非線性動力學(xué)降階方法。有限元模型含有8個線性Euler梁單元，1個非線性連接單元。兩個方向的Iwan模型用來描述連接界面上局部非線性恢復(fù)力。每個節(jié)點含有橫向平動u、垂向平動v、轉(zhuǎn)動θ三個自由度，共計30個自由度。線性梁單元的剛度矩陣、質(zhì)量矩陣、非線性連接單元的恢復(fù)力表達式可以參考文獻[9，12，38?39]。

局部非線性動力學(xué)響應(yīng)的轉(zhuǎn)換矩陣為

式中? 為連接梁單元的長度;為寬度。

圖5中，為殘余剛度系數(shù);兩個Iwan模型的臨界滑移力，黏著剛度分別定義為

式中? 為梁單元的彈性模量;為截面轉(zhuǎn)動慣量。

連接梁結(jié)構(gòu)的材料為鋼，密度ρ=7.85103 kg/m3，彈性模量為。激勵幅值為，激勵頻率范圍根據(jù)第一階和第三階彎曲模態(tài)頻率進行選取，分別為和，諧波階數(shù)取10（H=10），計算結(jié)果如圖6所示。

由圖6可知，本文方法預(yù)測的非線性動力學(xué)響應(yīng)與降階前的結(jié)果吻合較好，頻響函數(shù)能夠較好地反映非線性力學(xué)行為的影響，尤其是在共振峰附近的區(qū)域。由表2可知，降階之后，Jacobian矩陣的維數(shù)大大降低，非線性動力學(xué)響應(yīng)的計算耗時也大大減少，計算效率約提高70倍。

3 結(jié)? 論

本文提出一種基于局部非線性轉(zhuǎn)化的降階非線性動力學(xué)算法。采用Newmark和多諧波平衡法將非線性動力學(xué)微分方程轉(zhuǎn)化為非線性代數(shù)方程組，將整體結(jié)構(gòu)的動力學(xué)響應(yīng)轉(zhuǎn)化到僅與連接非線性相關(guān)的自由度上進行求解，減小迭代過程的計算耗費。利用含有Iwan模型的三自由度質(zhì)量彈簧振子和連接梁結(jié)構(gòu)驗證本文的非線性動力學(xué)降階方法。

本文方法預(yù)測的非線性動力學(xué)響應(yīng)與未降階的Newmark和多諧波平衡法的結(jié)果吻合較好。頻響函數(shù)能夠較好地反映遲滯非線性引起的能量耗散的影響，尤其是在共振峰附近。降階之后，非線性迭代過程的計算效率明顯提高，時域降階提高約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] 李玲， 蔡安江， 阮曉光. 栓接結(jié)合部非線性等效線性化方法[J]. 振動工程學(xué)報， 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] 張相盟. 摩擦連接的結(jié)構(gòu)非線性動力學(xué)研究[D]. 哈爾濱：哈爾濱工業(yè)大學(xué)， 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]. 力學(xué)學(xué)報， 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