簡諧激勵(lì)下一類干摩擦振子的響應(yīng)計(jì)算

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

（1.哈爾濱工業(yè)大學(xué) 衛(wèi)星技術(shù)研究所，150080哈爾濱；2.中國人民解放軍91216部隊(duì)，125106遼寧興城）

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

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

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

圖1 Iwan模型示意圖

本文研究的是由有限個(gè)Jenkins單元構(gòu)成的Iwan模型的摩擦振子在簡諧激勵(lì)下的響應(yīng)求解問題.針對此類分段線性的振子系統(tǒng)，在文獻(xiàn)［6］中方法的基礎(chǔ)上，給出了簡諧激勵(lì)下振子系統(tǒng)響應(yīng)的精確解析求解方法，并分析系統(tǒng)在不同量級的外激勵(lì)下的動力學(xué)特征和阻尼特性.

1 摩擦振子模型

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

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

圖2 振子模型示意圖

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

當(dāng)處于粘滯狀態(tài)時(shí)（即所有Jenkins單元均未滑動），系統(tǒng)的總剛度為

引入下面無量綱化參數(shù)

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

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

2 振子響應(yīng)計(jì)算

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

2.1 各階段的運(yùn)動方程

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

接觸面摩擦力可寫成

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

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

2.2 方程解的一般表達(dá)式

假設(shè)在τa時(shí)刻，運(yùn)動進(jìn)入階段Si，則方程（8）的解可寫成下面一般形式

其中ηh為瞬態(tài)部分，ηp為穩(wěn)態(tài)部分，它們各自為

式中

很明顯，當(dāng)λi=1時(shí)，Bi達(dá)到最大值1／（2~），亦即1／（ξΩ）.

當(dāng)初值ηi，0與η′i，0已知時(shí)，便可求得該運(yùn)動階段內(nèi)任意時(shí)刻的η（τ）值，代入式（7）便可進(jìn)一步求得位移y（τ）.假設(shè)Si的前一個(gè)運(yùn)動階段的系統(tǒng)響應(yīng)已知，由速度和位移在τa處的連續(xù)性，容易求得初值ηi，0與η′i，0分別為

2.3 摩阻片位移值計(jì)算

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

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

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

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

2.4 階段轉(zhuǎn)換時(shí)刻的確定

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

運(yùn)動進(jìn)入階段Sj，如果速度為零

則所有Jenkins單元全部停止滑動，運(yùn)動轉(zhuǎn)入S0.為了確定狀態(tài)轉(zhuǎn)換時(shí)間，我們構(gòu)造下面函數(shù)

以及

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

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

3 仿真算例及結(jié)果討論

假設(shè)Iwan模型中Jenkins單元的剛度和屈服位移滿足下面關(guān)系：

模型參數(shù)見表1.

表1 模型參數(shù)

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

圖4 振子無量綱位移、速度和恢復(fù)力時(shí)間歷程圖

圖5 振子穩(wěn)態(tài)階段一個(gè)周期內(nèi)位移、速度和摩擦力的時(shí)域曲線

圖6 系統(tǒng)穩(wěn)態(tài)段遲滯環(huán)

圖7是不同激勵(lì)幅值下，激勵(lì)頻率在［0.1，1.5］區(qū)間內(nèi)振子系統(tǒng)響應(yīng)的幅頻曲線.圖中，γ=Fe／1.579 1.結(jié)果顯示，當(dāng)γ=0.2時(shí)，其幅頻曲線峰值出現(xiàn)在Ω=1處，其值為20，正好等于Bi的極值1／（ξΩ），表明在此激勵(lì)下，整個(gè)振動周期內(nèi)F（y′，｛zi｝）均為零，表明振子運(yùn)動一直處于S0階段，因此振子運(yùn)動是完全線性的.故該幅頻曲線也是線性的.

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

對A求偏導(dǎo)，得到

可見，在區(qū)間（xi?，+∞）上，ceq，i隨A先增后減，當(dāng)A的值為2xi?時(shí)，ceq，i有最大值為

圖7 不同量級外激勵(lì)下的振子響應(yīng)的幅頻曲線

圖8 單個(gè)Jenkins單元的遲滯回線示意圖

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

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

圖9 當(dāng)γ＝15，Ω＝0.4和Ω＝0.6時(shí)分別對應(yīng)的遲滯環(huán)

4 結(jié) 論

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

［1］肖世富，陳濱，杜強(qiáng).寬帶隨機(jī)激勵(lì)下非線性連接結(jié)構(gòu)振動響應(yīng)分析［C］／／第十三屆全國結(jié)構(gòu)工程學(xué)術(shù)會議論文集（第Ⅰ冊）.井岡山：中國力學(xué)學(xué)會結(jié)構(gòu)工程專業(yè)委員會，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，F(xiàn)ENG 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.