簡諧激勵下一類干摩擦振子的響應計算

王本利，張相盟，劉 源，馬 平

（1.哈爾濱工業大學 衛星技術研究所，150080哈爾濱；2.中國人民解放軍91216部隊，125106遼寧興城）

裝配結構中，大量存在著各種機械連接，諸如螺接、鉚接、法蘭式連接或其他緊固式連接，這些連接也被稱為“硬連接”［1］.實際上，此類連接并不完全緊固，當外激勵量級較大時，其連接接觸面常會產生微滑移（Microslip）甚至宏滑移（Macroslip），與之伴隨的干摩擦是結構阻尼的一個重要來源，這種干摩擦阻尼甚至可占到整個結構阻尼的90%［2］.由于這種連接還會改變連接處的局部剛度，從而會影響整個結構的動力學特性及動態響應［3］.因此，在進行整體裝配結構計算時，需要對連接處特別考慮.

庫倫模型是一種簡單常用的摩擦模型，它已被廣泛應用于描述摩擦阻尼或摩擦連接處的建模和計算中.這方面的工作可以追溯到Den Hartog［4］，他首次計算了同時含有庫倫阻尼器和粘性阻尼器的振子系統的穩態響應.Lee等［5］計算了含有一個摩擦連接的梁在正弦激勵下的穩態解，其連接處摩擦力用庫倫模型來描述，得到的數值計算結果與實驗結果吻合.Ding等［6］提出了用于計算含有庫倫摩擦阻尼器的單自由度振子系統響應的解析計算方法.

庫倫模型只能用來描述接觸面的粘滯狀態和宏滑移狀態，而不能用于描述微滑移狀態.而實際上，在連接處若產生宏滑移，常會引起連接失效，其并不多見；而微滑移現象更為普遍，它是結構阻尼的重要來源且不會引起連接失效［7-8］，因此其更具研究價值.目前，國內外學者已發展了多種可以描述接觸面宏滑移和微滑移的遲滯模型，如Iwan模型［9］、Valanis模型［10］以及Bouc-Wen模型［11-12］等.其中，Iwan模型由于其構型簡單直觀，近年來頗受關注，已被廣泛應用于研究含摩擦連接結構的力學行為和阻尼問題［13-16］.圖1（a）和（b）分別是經典Iwan模型和改進Iwan模型構型示意圖.經典Iwan模型是由N個Jenkins單元并聯而成，每個Jenkins單元（也稱為雙線性遲滯模型）是由剛度為ki的線性彈簧和屈服力（即臨界摩擦力）為fi?的庫倫摩阻片串聯構成.修正Iwan模型是在經典Iwan模型的基礎上并聯一個彈簧單元得到的.由于每個Jenkins單元的屈服位移xi?（xi?=fi?／ki）不同，使得Iwan模型可以用來描述連接處的宏滑移行為（所有Jenkins單元產生滑動）和微滑移行為（部分Jenkins單元出現滑動）.

圖1 Iwan模型示意圖

本文研究的是由有限個Jenkins單元構成的Iwan模型的摩擦振子在簡諧激勵下的響應求解問題.針對此類分段線性的振子系統，在文獻［6］中方法的基礎上，給出了簡諧激勵下振子系統響應的精確解析求解方法，并分析系統在不同量級的外激勵下的動力學特征和阻尼特性.

1 摩擦振子模型

本文研究的振子模型如圖2所示，模型為剛度kα，粘性阻尼c，質量塊m的振子系統.在外部簡諧激勵Fesin（ωt）的作用下，質量塊將沿摩擦面上運動，位移為x.基于Iwan模型描述的摩擦接觸模型見圖3：圖中，ui為Iwan模型中任意一個摩阻片的位移.系統的運動方程為

其中，｛ui｝是各摩阻片位移的集合，f（x，x，｛ui｝）是Iwan模型描述的接觸面摩擦力函數，如下

圖2 振子模型示意圖

圖3 含Iwan摩擦模型的振子模型

當處于粘滯狀態時（即所有Jenkins單元均未滑動），系統的總剛度為

引入下面無量綱化參數

得到無量綱化后的方程（1）為

其中y＇和y″分別表示y對τ的一、二階導數，F（y，y′，｛zi｝）為無量綱摩擦力函數式如下所示.

2 振子響應計算

從式（3）可以看出，可根據Iwan模型中發生屈服（滑動）的Jenkins單元數量，將振子每次加載（振子速度為正）／卸載（振子速度為負）的過程分為N+1個階段，分別命名為S0，S1，S2，…，SN，其中下標i（i=0，1，2，…，N）對應為在該運動階段已屈服的Jenkins單元的數量.在每個階段內，振子的剛度不變，其恢復力-位移關系是線性的.在周期外載荷激勵下，振子反復經歷加載和卸載過程，在任意一次加載／卸載過程中，如果接觸面存在滑動，系統將歷經多個運動階段.下面任取一個加載／卸載過程中的一個階段的運動進行分析.

2.1 各階段的運動方程

將所有Jenkins單元的屈服位移按從小到大排列，則對于任意運動階段Si，無量綱位移滿足

接觸面摩擦力可寫成

式中F（y＇，｛zi｝）為常數項，Kires為階段Si中Iwan模型的剩余剛度，兩者分別為

將式（5）代入式（2），得到階段Si的運動方程為

2.2 方程解的一般表達式

假設在τa時刻，運動進入階段Si，則方程（8）的解可寫成下面一般形式

其中ηh為瞬態部分，ηp為穩態部分，它們各自為

式中

很明顯，當λi=1時，Bi達到最大值1／（2~），亦即1／（ξΩ）.

當初值ηi，0與η′i，0已知時，便可求得該運動階段內任意時刻的η（τ）值，代入式（7）便可進一步求得位移y（τ）.假設Si的前一個運動階段的系統響應已知，由速度和位移在τa處的連續性，容易求得初值ηi，0與η′i，0分別為

2.3 摩阻片位移值計算

式（10）中變量｛zi｝尚未求得.將構成Iwan模型的Jenkins單元按屈服位移從小到大依次命名為J1、J2、…、JN，并將Si所處的加載／卸載過程的初始時刻和末端時刻分別標記為τ01和τ02.則在運動階段Si中，Jenkins單元J1至Ji屈服，其余均處于粘滯狀態.對于處于粘滯狀態的單元，其摩阻片在該階段任意時刻的位移值與加載／卸載初始時刻的位移值相同：

而對于已滑動的單元，在Si階段滑動的位移與振子在該階段振子運動位移相同，因此

式（15）表明，每次加載／卸載過程中，初始時刻各摩阻片的位移值正是式（10）中待求的｛zi｝，因此只需求得這些時刻已屈服單元中摩阻片的位移值即可.假設當前加載／卸載過程中，振子的運動經歷了S0到Sn（n≤N）的階段，且zj（τ01）已知，結合式（13）不難求出τ02時刻（下一個卸載／加載初始時刻）各摩阻片的位移值為

式（14）也可視為每次加載／卸載完成后，庫倫摩阻片位移值的更新.由于零時刻點摩阻片的位移值已知，則第一個加載／卸載過程內振子的位移值便可求得，下一個卸載／加載過程初始時刻摩阻片的位移值可通過式（14）～（15）求得，依次遞推求解，便可獲得整個時域內振子響應和每次加載／卸載初始時刻各個摩阻片的位移值.

2.4 階段轉換時刻的確定

確定階段轉換時刻，就是確定各個運動階段的初始時刻和末端時刻，如式（9）～（11）中的τa，該值顯然會影響解的精度.當不等式（4）不再滿足時，運動跳出階段Si，轉入其他運動階段.這里存在兩種情況，如果

運動進入階段Sj，如果速度為零

則所有Jenkins單元全部停止滑動，運動轉入S0.為了確定狀態轉換時間，我們構造下面函數

以及

不難看出，確定狀態轉換時刻實際上就是尋找g1（τ）或g2（τ）的零值點所對應的時刻值，亦即為求方程g1（τ）=0或g2（τ）=0的根.在每一步的計算過程中，應首先判斷g2（τ）是否經過零值點：如果過零，則狀態轉入S0；如果未過零，進一步判斷g1（τ）：如果過零，則轉入更高階段Sj，如果未過零，則進行下一步計算.為了獲得更準確的兩函數零值點所對應的時間值，當g1（τ）或g2（τ）處于零值附近時，應減小步長，可用二分法確定零值點所對應的時刻值.

從前面分析過程不難看出，本文提出的方法是一種“分段-解析”法：先對分段線性運動方程（2）進行分段分析，由于各運動階段內的運動方程（6）都是線性的，通過變量代換（式（7）），便獲得了代換后方程的精確解析解（9），這樣順次求得相應時域內各個階段運動方程的解析解，自然獲得了所求時域上振子系統的響應.由于式（9）為線性方程（8）完全精確的解析解，而方程（8）是通過對方程（6）進行變量代換得到的，變量代換顯然不會引入誤差，因而將式（9）代入代換式（7）所得的解便是方程（6）的精確解析解.由此可見，如果階段轉換時刻能精確確定，由各階段的精確解析解連接起來所形成的方程（2）的解便完全精確.但由階段轉換時刻只能通過數值方法確定，這個過程不可避免會產生數值誤差，其誤差大小由事先設置的誤差容限（迭代終止條件）決定.因此，這里的“分段-解析”法的誤差僅來源于確定階段轉換時刻的數值誤差，這種誤差顯然要遠小于求解非線性振動方程的近似解析方法（如諧波平衡法、小參數攝動法等）由近似所致的誤差.與純數值方法相比（如用四階龍格庫塔法直接求解方程（2）），該“分段-解析”法可給出每個運動階段所對應方程的精確解析解，因此比數值方法更精確和高效.需要說明的是，當Jenkins單元數量非常大而趨于無窮時，其振子系統不再是分段線性系統，這種情況下，本文方法不再適用.

3 仿真算例及結果討論

假設Iwan模型中Jenkins單元的剛度和屈服位移滿足下面關系：

模型參數見表1.

表1 模型參數

取方程參數為表1中的數據，由上節介紹的方法求得的大約30個周期的無量綱位移、速度和摩擦力的時間歷程如圖4所示.從圖中可以看出，由于摩擦阻尼的影響，使得系統運動很快進入穩態.圖5為圖4中位移和摩擦力穩態階段一個周期時域曲線的放大圖，曲線上的圓圈為階段轉換點.可以看出，在該量級的激勵下，在一個加載／卸載過程中依次經歷了S0到S4的5個運動階段.從圖還可看出，每個轉換點前后位移曲線光滑性要好于恢復力曲線，這是由于恢復力變化量是剛度與位移變化量的乘積，因此階段轉換引起的剛度變化對其影響較大.以圖5中位移和摩擦力分別為x軸和y軸，便得到了如圖6所示的遲滯環曲線，曲線的斜率值表示了每個運動階段內系統的剛度.很明顯，在S4階段，Iwan模型恢復力值為定值，Iwan模型剛度為零，表明Iwan模型中所有Jenkins單元均屈服，模型處于宏滑移階段.

圖4 振子無量綱位移、速度和恢復力時間歷程圖

圖5 振子穩態階段一個周期內位移、速度和摩擦力的時域曲線

圖6 系統穩態段遲滯環

圖7是不同激勵幅值下，激勵頻率在［0.1，1.5］區間內振子系統響應的幅頻曲線.圖中，γ=Fe／1.579 1.結果顯示，當γ=0.2時，其幅頻曲線峰值出現在Ω=1處，其值為20，正好等于Bi的極值1／（ξΩ），表明在此激勵下，整個振動周期內F（y′，｛zi｝）均為零，表明振子運動一直處于S0階段，因此振子運動是完全線性的.故該幅頻曲線也是線性的.

比較圖中的曲線可發現：隨著激勵力幅值增大，曲線峰值明顯左移，表現出剛度軟化特征，其原因顯而易見：激勵力幅值越大，產生屈服Jenkins單元越多，系統剛度遞減量就越大.從圖中還可發現，隨著激勵力幅值的增大，幅頻曲線的峰值先減小后增大，表明系統阻尼也隨之先減小后增大.文獻［17］中的實驗結果也出現了此現象，但作者并未給出解釋.顯然這種阻尼特征是由Iwan模型決定的，由于Iwan模型是由一系列Jenkins單元并聯而成，故其阻尼特性與單個Jenkins單元的阻尼特性密切相關.對于任意一個Jenkins單元，在簡諧激勵下發生屈服時的遲滯環的形狀如圖8所示，圖中A為振幅.當振動頻率為ω時，Jenkins單元的等效粘性阻尼為

對A求偏導，得到

可見，在區間（xi?，+∞）上，ceq，i隨A先增后減，當A的值為2xi?時，ceq，i有最大值為

圖7 不同量級外激勵下的振子響應的幅頻曲線

圖8 單個Jenkins單元的遲滯回線示意圖

而Iwan模型的等效粘性阻尼為

可見，Iwan模型的等效粘性阻尼也會隨著振幅的增大而先增大后減小.γ=15對應的幅頻曲線上共振帶邊緣頻率點Ω=0.4和Ω=0.6處所對應的穩態段的遲滯環形狀如圖9所示.可以看出，在一個運動周期的大部分時間內，Iwan模型處在宏滑移階段（即S4階段），遲滯環形狀與雙線性遲滯模型的遲滯環形狀（圖8）相似，且振幅遠大于Iwan模型中屈服位移值最大的J4的屈服位移的2倍，因此其等效阻尼值已經出現減小.

圖9 當γ＝15，Ω＝0.4和Ω＝0.6時分別對應的遲滯環

4 結 論

針對簡諧激勵下一類含Iwan模型的分段線性振子系統的非線性振動問題，提出了用于計算系統響應的“分段-解析”方法.通過對一算例求解，分別得到了系統響應的時域曲線和幅頻曲線.在激勵量級逐漸遞增時，幅頻曲線峰值明顯向左偏移，使得Jenkins單元滑動所致的剛度軟化效應得到了體現.另一方面，摩擦阻尼對系統響應影響也很顯著：隨著激勵量級的遞增，幅頻曲線峰值先減小后增大.通過對構成Iwan模型的Jenkins單元的等效粘性阻尼分析，得到系統等效粘性阻尼隨振幅的增大而先增大后減小，從而解釋了幅頻曲線峰值產生這種變化的原因.文中方法為存在分段線性問題的系統的響應求解提供了一種有效途徑.后面還可進一步開展關于存在摩擦連接邊界的連續體的非線性振動方面的研究.

［1］肖世富，陳濱，杜強.寬帶隨機激勵下非線性連接結構振動響應分析［C］／／第十三屆全國結構工程學術會議論文集（第Ⅰ冊）.井岡山：中國力學學會結構工程專業委員會，2004：405-411.

［2］BEARDS C F.Damping in structural joints［J］.The Shock and Vibration Digest，1992，24（7）：3-7.

［3］IBRAHIM R，PETTIT C.Uncertainties and dynamic problems of bolted joints and other fasteners［J］.Journal of Sound and Vibration，2005，279（3／4／5）：857-936.

［4］DEN HARTOG J P.Forced vibrations with combined Coulomb and viscous friction［J］.Journal of Applied Mechanics，ASME，1931，53（9）：107-115.

［5］LEE Y，FENG Z C.Dynamic responses to sinusoidal excitations of beams with frictional joints［J］.Communications in Nonlinear Science and Numerical Simulation，2004，9：571-581.

［6］DING Q，CHEN Y.Analyzing resonant response of a system with dry friction damper using an analytical method［J］.Journal of Vibration and Control，2008，14（8）：1111-1123.

［7］GUTHRIE M A，KAMMER D C.A general reduced representation of one-dimensional frictional interfaces［J］.Journal of Applied Mechanics，ASME，2008，75（1）：011019.

［8］CHEN W，DENG X.Structural damping caused by micro-slip along frictional interfaces［J］.International Journal of Mechanical Sciences，2005，47（8）：1191-1211.

［9］IWAN W D.A distributed-element model for hysteresis and its steady-state dynamic response［J］.Journal of Applied Mechanics，ASME，1966，33（4）：893-900.

［10］VALANIS K C.Fundamental consequences of a new intrinsic time measure：plasticity as a limit of the endochronic theory［J］.Archiwum Mechaniki Stossowanej，1980，32（2）：171-191.

［11］BOUC R.Forced vibration of mechanical system with hysteresis［C］／／Proceedings of the Fourth Conference on Nonlinear Oscillations.Prague：［s.n.］，1967.

［12］WEN Y K.Method of random vibration of hysteretic systems［J］.Journal of the Engineering Mechanics Division，ASCE，1976，102（2）：249-263.

［13］OUYANG H，OLDFIELD M J，MOTTERSHEAD J E.Experimental and theoretical studies of a bolted joint excited by a torsional dynamic load［J］.International Journal of Mechanical Sciences，2006，48（12）：1447-1455.

［14］SEGALMAN D J.A modal approach to modeling spatially distributed vibration energy dissipation，Technical Report No.SAND2010-4763［R］.New Mexico：Sandia National Laboratories，2010.

［15］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.

［16］QUINN D D.Modal analysis of jointed structures［J］.Journal of Sound and Vibration，2012，331（1）：81-93.

［17］JALALI H，AHMADIAN H.Identification of microvibro-impacts at boundary condition of a nonlinear beam［J］.Mechanical Systems and Signal Processing，2011，25：1073-1085.