條(環)狀干摩擦阻尼器的微滑移數值模型*

趙寧，藺彥虎(西北工業大學機電學院，陜西西安710072)

條(環)狀干摩擦阻尼器的微滑移數值模型*

趙寧，藺彥虎
(西北工業大學機電學院，陜西西安710072)

提出一種求解彈性條(環)狀阻尼器微滑移接觸運動的數值方法。將阻尼器和外部激勵歷程在空間和時間上離散，將相同數量的干摩擦觸點布置于離散阻尼器上；把接觸運動判據應用到各離散接觸點，確定其運動狀態并修正剛度矩陣，求解整個阻尼器的平衡方程。該方法避免了有限元軟件求解含摩擦接觸問題的迭代過程，從而保證了求解的可執行性。同時，克服了微滑移模型理論解法對法向載荷分布規律及載荷時變性的限制，為求解具有局部性以及時變性的法向載荷的結構動態響應提供了更為精確的邊界條件，從而可提高結構頻響分析的準確性。應用多諧波平衡法分別計算宏滑移和微滑移阻尼器約束下的結構動態響應，發現在結構減振中，微滑移模型能夠適應更寬范圍的法向力。

干摩擦阻尼器；微滑移；接觸運動；多諧波平衡法；動響應

由于結構簡單、在高溫及復雜環境的穩定性及優良的減振性能，干摩擦阻尼器被廣泛應用于工程中。如為了降低筒形結構的航空發動機的篦齒封嚴裝置振動幅值，在其定子或轉子上安裝有阻尼環或阻尼套筒；航空發動機的葉-盤系統中，葉片設計有凸臺，通過振動時凸臺間的擠壓和相對運動來降低葉片振動幅值；航空薄壁齒輪上安裝的阻尼環等結構［1］，均是應用干摩擦進行結構減振的實例。

葉片減振結構中，凸臺面積以及凸臺間的法向力較小，可將其簡化為單對觸點，進行非線性結構的動響應求解［2］。然而，對于安裝阻尼環、阻尼套筒的結構，由于接觸面積大，加上筒狀結構進行周波型振動或者盤狀結構進行節徑型振動時，將引起不同接觸部位的法向力變化，局部法向力可能會變得較大，這時還將阻尼器與結構間的作用通過單對接觸點來表示，將不符合阻尼器與結構間的實際。

微滑移模型用多對觸點來表達其與結構間的局部摩擦作用，Menq［3-4］建立了考慮微滑移效應的矩形桿模型，桿與結構間分布有厚度可忽略的剪切層，用該剪切層模擬阻尼器和結構表面微觀凸起，給出了接觸運動以及干摩擦力的理論解［3］，并將該理論解應用于結構動態響應計算，并與實驗結果［4］對比，驗證了該模型的準確性。徐自力［5］將該模型引入葉片減振分析，并分析了最優正壓力的影響因素。Csaba［6］的微滑移模型忽略了Menq模型中的剪切層，用單元節點切向剛度模擬阻尼器的切向剛度。Menq與Csaba模型均采用解析法求解，法向力分布為特定形式，不可時變，激勵力或激勵位移也被設定為諧波形式。這些假設往往不符合應用中的載荷分布以及變化特性，具有一定的局限性。

商用有限元軟件應用迭代法求解含摩擦的接觸問題，迭代過程實際上是確定觸點運動狀態的過程。有限元軟件能較好地模擬含摩擦的局部接觸作用，并對結構做時域內的響應分析，然而，多數商業有限元軟件，如Ansys、Abaqus等均不能用來做非線性頻響分析，而非線性頻響求解又是薄壁結構減振設計的重要內容。馬曉秋［7］用諧波法描述干摩擦力，將干摩擦力等效為結構內阻，用Ansys軟件計算了葉片頻響。然而一階諧波法是一種近似方法，不能準確表達接觸運動的粘-滑效應，等效方法對求解的準確性也有一定影響。因此有必要針對阻尼器具體結構，發展出能夠方便求解接觸運動的數值模型，將其應用于非線性頻響分析，擺脫應用有限元軟件求解接觸問題的束縛。

本文將Yang［8］發展的宏滑移模型黏滯—滑移—分離判據應用于微滑移模型，形成了求解微滑移接觸運動的數值方法。相比微滑移模型理論解法，該方法對微滑移模型法向載荷分布規律及其時變性沒有限制，同時取消了Menq［3］模型中防止剛體位移的限位彈簧，可對剛體位移進行準確求解。求解過程無須迭代，從而保證求解的可執行性。最后，應用多諧波平衡法(Multiple Harmonic Balance Method，MHBM)計算微滑移數值模型約束下的結構動響應，以證實所發展模型的實用性，并分析了宏滑移模型和微滑移模型對結構動響應影響的異同。

1 Yang模型及算例

1.1 Yang宏滑移模型

相比于微滑移模型，對宏滑移模型的研究較多，也較為成熟。Sanliturk和Ewins［9］發展了一種求解宏滑移模型的“軌跡跟蹤法”，該方法將求解周期離散為若干時間點，對每一時間步的運動狀態進行判斷，確定運動狀態后，計算出相應狀態下的位移量及干摩擦力。重復進行以上步驟，即可得到整個周期內的觸點運動狀態及干摩擦力，一般地，跟蹤過程只需維持兩個周期，即可得到穩定的運動狀態。單穎春［10］應用類似方法計算葉片動響應，并與實驗值對比，取得了很好的效果。

Yang于1998年發表文章［8］詳細討論了宏滑移模型的運動狀態確定以及狀態判定依據，并假設激勵為諧波形式時運動狀態的轉換角。所建立的滑移模型如圖1所示。

圖1中，body1和body2表示兩接觸物體，通過觸點表示兩者的作用力，body2可運動也可固定。n0表示初始法向力；ku，kv分別為切向及法向剛度；u，v分別為切、法向運動；f表示干摩擦力；μ為干摩擦系數；w表示觸點切向運動；n表示法向力。兩接觸體可能存在3種運動狀態，黏滯、滑移或分離狀態，在接觸運動周期內，這幾種狀態可能會相互轉化，如表1所示。

圖1 Yang接觸運動模型Fig.1 Yang’s contact kinematicmodel

表1中，“E”表示黏滯狀態，“P”表示正向滑動，“N”表示反向滑移，“S”表示分離狀態，“當前”表示當前時刻觸點的運動狀態，“后繼”表示運動狀態從當前到下一時刻的轉化。

表1 運動狀態轉化條件Tab.1 Translation criterion of kinematic state

干摩擦力遵循Coulomb摩擦定律，即

式中，sgn為符號函數，abs為絕對值符號，t為時間。

1.2 宏滑移接觸運動模型算例

選取圖1中“body1”的運動參數如表2所示，“body2”固定不動，求解該模型。求解時將周期離散，按表1中所示的“當前狀態”，計算運動狀態參數及力值，然后依據表1中的不等式，判定是否發生了運動狀態轉化，如果轉化，則按照新的狀態求解各參量，如果未發生狀態轉化，則計算下一個時刻各參數。計算取摩擦系數μ為0.5，t代表離散時間點。

結果如圖2所示，圖中每列為表2中計算參數的運動狀態、力及滯回環。可以看出應用Yang發展的判據可求解復雜的運動狀態。從圖2第二行可以看出在運動過程中，兩物體遍歷了分離—黏滯—滑移運動，所以該行的第三列所示的滯回環中有相當一部分值為0。第三行表示在平面內的運動頻率不同時的運動狀態，類似于李薩如的運動形式，這時摩擦力及滯回曲線顯得更為復雜，平面內y向運動頻率為x向的3倍，所以滯回環也呈現出這種關系。

表2 宏滑移模型計算參數Tab.2 Parameters ofmacro-slip model

圖2 運動狀態、力及滯回曲線Fig.2 Curve of kinematic state，force and hysteresis

2 微滑移數值模型及求解方法

2.1 條(環)狀阻尼器的離散

工程中經常用到的阻尼器為環狀或條狀結構，如阻尼環、阻尼塊。干摩擦模型的建立是為了研究法向力存在時的切向運動及受力，故適合用桿模型等效條狀阻尼器［3］。圖3(a)為Menq［3］建立的微滑移模型，在阻尼塊與結構的接觸面上布有剪切層，用來模擬接觸面的切、法向力學行為。該模型認為，當接觸體表面粗糙度足夠小時，將接觸表面均勻劃分為若干塊，盡管各分塊的凸起和凹陷形貌不盡相同，然而，從統計學的角度看，他們表現出來的力學行為是相同的。因此剪切層的不同部位具有相同的切、法向剛度。

用數值方法對微滑移理論模型求解的第一步是對阻尼器離散。根據求解問題，在Menq模型底部預設一定數量觸點來替代剪切層，將模型均勻離散，離散數量與觸點數相同，如圖3(b)所示。由于均勻離散，各觸點的力學特征相同，所以觸點切向剛度Ku相等。條(環)狀阻尼器法向尺寸遠小于接觸平面尺寸，法向力引起的彎曲效應較小，可將法向力作用直接施加在觸點上。若阻尼器為開口環或條狀阻尼器，則連接第n塊與第1塊的彈簧將不存在，對于整體環，則彈簧存在。圖中fn，x(t)(x=1，2，…，n)表示法向載荷，該載荷可時變。k為離散單元間的剛度。

圖3 Menq微滑移模型及其離散Fig.3 Menq’smicro-slipmodel and its discretization

2.2 求解流程

與宏滑移模型相似，求解微滑移模型也需要在求解周期內離散。

如圖3(b)所示，假設在第1號離散體上作用位移載荷，則該位移載荷也可被離散為有限多個值，用uj和d uj表示j時刻的位移以及j+1時刻與j時刻的位移變化量。用Fi，j表示觸點i在j時刻的受力，d Fi，j表示第j+1時刻與j時刻觸點i受力變化量。類似地，用wi，j和d wi，j表示觸點位移量及其變化量，用xi，j和d xi，j表示離散段的位移量及其變化量。用s表示u，F，x，w則有

d x為一向量，表示離散位移增量，長度等于離散體個數。k(t)為離散體剛度矩陣，呈帶狀，是時變量:黏附于離散體的觸點(如第q個觸點)處于黏滯時，剛度矩陣對應的k(q，q)值為k(q，q)= 2k+ku，ku為剪切剛度，如圖3(b)所示。若處于滑移或分離狀態時，對應節點的剛度值為k(q，q) =2k。需要注意的是，由于激勵形式為位移激勵，因此在任何時間點上任何離散塊體都處于受力平衡狀態，即每個離散塊體所受合力為0，故f為0向量，注意式(3)中的f是各離散體所受合力，不是摩擦觸點的受力。求解結果有可能出現阻尼器的剛體位移，這取決于位移載荷的大小和法向力值。圖4給出求解流程，圖中i表示干摩擦觸點，j表示離散時間點。

圖4 微滑移模型計算流程Fig.4 Flowchart of solving micro-slip model

3 算例

3.1 算例1均布載荷的整體環狀阻尼器

圖5所示的整體環狀阻尼器，與圖3(a)相比，除了結構與載荷變化外，做了兩個調整:一是激勵形式變為位移形式，帶來的好處是可以求解剛體位移，既可求解出整體滑移也可求解出局部滑移。在阻尼器實際工作中，也屬位移激勵的情形，摩擦力伴隨接觸運動產生。二是解除了防止阻尼器產生剛體位移的彈簧β，對原模型的假設限制進一步減少。

圖5 均布載荷的環狀阻尼器Fig.5 Circular damper with uniform load

如圖5所示，阻尼環橫截面為矩形，彈性模量E=2.07×105MPa，寬b=30mm，阻尼環內徑r= 140mm，外徑R=150mm，μ=0.5，切向剛度kd按Mindlin［11］理論計算。激勵取x(t)=0.31sin(2πt)，方向沿環的切向，法向力q(t)=137N，t為離散時間點，將結構均勻離散，得到圖6所示的滯回環。

如圖6所示，離散數取n=1時，阻尼器模型為宏滑移模型，由于未考慮阻尼器本身的彈性變形，其滯回環的斜率最大。其求解結果與跟蹤求解［9］方法一致，當n逐步增大時，求解結果越逼近微滑移模型理論解，當離散數n=15時，數值模型的滯回曲線與理論解重合，已經足夠接近理論解［3］了。在激勵歷程中，微滑移觸點逐次滑移，最終達到整體滑移，形成剛體運動。

圖6 微滑移模型滯回環Fig.6 Hysteresis curve ofmicro-slip model

3.2 算例2非均布載荷

如圖7所示，阻尼器橫截面為矩形，楊氏模量與寬度，即干摩擦系數同算例1，高h=3mm，激勵x(t) =0.32sin(2πt)，長L=80mm；載荷呈拋物線狀，最大值分布于兩端，均為34.5N，最小值為5N。

圖7 非均布載荷模型Fig.7 Modelwith non-uniform load

求解阻尼器離散為15段，得到如圖8所示的滯回環。從滯回環可以看出，該模型在激勵下并未發生剛體滑動，只是產生了局部滑移。

圖8 非均布載荷下的滯回環曲線Fig.8 Hysteresis curve of non-uniform load

圖9滑移等值線圖Fig.9 Contourmap of slip

圖9 為各摩擦觸點在一個運動周期內的滑移量等值線圖，可以看出8號觸點的滑移量最大，成為滑移核心，結合圖7，可以發現這是由于法向載荷較小造成的；等值線圖不是關于點8對稱而是偏向于編號小的觸點，這是由于計算中考慮了阻尼器的變形以及切向載荷的位置共同造成的。1～4號觸點及11～15號觸點一直處于黏滯狀態，所以模型呈部分滑移狀態，與圖8的結論是一致的。

該計算流程還可以求解法向力時變條件下的微滑移接觸運動及干摩擦力，也可計算有分離運動的實例。限于篇幅，這里不再給出算例。

4 MHBM法計算計及微滑移的結構動響應

具有干摩擦約束的結構呈現出強邊界非線性，傳統的數值積分方法由于計算耗時太長而不能滿足結構設計過程所要進行的頻響分析。因而人們發展出了攝動法，諧波平衡法，加遼金法等。經研究［12］，應用MHBM法結合Fourier變換及其逆變換法求解含干摩擦的動響應問題，具有足夠的準確度和求解速度。本文用該方法結合快速Fourier變換及快速逆Fourier變換法，計算計及微滑移的結構動態響應。

4.1 多諧波平衡法

圖10為一裝有微滑移阻尼器的振動模型，該模型可用來表示航空發動機渦輪葉盤B-D型減振器。M，C，K分別表示質量、阻尼和剛度。fex(t)為激振力，q為阻尼器法向載荷。其控制方程為

式中fnl［x，˙x，t］為由阻尼器產生的非線性力。諧波平衡假設認為，在周期性載荷作用下，系統的響應也是周期的，因此可用傅里葉序列表示

圖10 裝有微滑移阻尼器的振子模型Fig.10 A vibratormodelwith a micro-slip damper

式中Λk=-(kω)2M+j kωC+K，為動剛度矩陣。

求解時，首先對位移序列傅里葉系數進行傅里葉逆變換，得出時域位移值，用該位移序列求解干摩擦力，之后對干摩擦力進行傅里葉變換，得到干摩擦力在頻域內的諧波系數值。最后對位移序列的傅里葉系數值迭代，直到滿足式(8)為止。式(8)為一非線性方程，需用迭代方法求解，Newton法及其衍生格式是求解該方程的基礎。本文采用Broyden法求解。

4.2 計算結果

選取表3所列的計算參數，分別應用微滑移模型和宏滑移模型計算結構動響應。

表3 計算參數表Tab.1 Parameters used in the calculation

微滑移模型計算的頻幅響應如圖11所示，圖中Fn為法向載荷，即圖10中的q，計算微滑移模型離散為15段，切向剛度求解用Mindlin［11］理論。

圖11 微滑移模型結構頻響圖Fig.11 Frequency response of structure withmicro-slip model

求出宏滑移模型頻響后，即可對比兩者減振性能的差異。如圖12所示，法向力較小時，微滑移模型與宏滑移模型的效果相差無幾，微滑移模型對振動的抑制效果稍差于宏滑移模型，這主要是微滑移模型考慮了阻尼器的彈性形變而顯示出切向剛性不足，不能給結構提供足夠的干摩擦力約束振動幅值。隨著法向載荷的增加，則顯示出大的差異:采用宏滑移模型時，阻尼器很快處于黏滯狀態，這時可等效為一彈簧，因此隨著法向力的進一步增加，結構振動幅值基本不再增加；而微滑移模型考慮了結構的局部滑移，阻尼器仍未完全黏滯，因此可以提供不同幅值的干摩擦約束力，從圖6也可以觀察出這一點。直到法向載荷增加至1007N，微滑移模型所有觸點才完全黏滯，達到與宏滑移模型相同的振動幅值。由于在較大法向載荷下仍有觸點可以滑動，因此微滑移模型在較大法向載荷時，仍具有較明顯的減振效果。

圖12 宏、微滑移模型減振效果對比Fig.12 Contrast between themacro-slip model and themicro-slip model

Griffin［13］的用宏滑移遲滯模型為航空發動機葉片B-D型阻尼器建模，計算了葉片的頻響及動應力，并搭建實驗臺驗證計算結果的可靠性。計算及實驗結果如圖13所示［13］。圖13為兩種不同剛度的阻尼器在不同法向力下葉片的動應力。文獻［13］未給出葉片及阻尼器的幾何及力學參數，因此無法用本文發展的微滑移模型計算其動響應。但觀察其計算結果與實驗結果可發現:法向力值較低時，實驗結果和計算結果相符合；但隨著法向力的增加，實驗值與計算值出現偏離，尤其是在阻尼器剛度較小的條件下，并且法向力較大時，動應力的計算值總是大于實驗值。

切向剛度較大時，阻尼器更像是一個剛體，因而可用宏滑移模型描述，其計算結果與實驗值相差不大；當切向剛度較小時，阻尼器在整體滑動之前更容易產生局部滑移，這時如果還用宏滑移模型來描述阻尼器，則會在法向力較大時無法計及局部滑移而使得計算值偏大。對比圖12，當法向力增大時，使用宏滑移模型計算的結構頻響值會很快上升，而使用微滑移模型的頻響值上升緩慢。可見用微滑移模型，動應力計算結果會緩慢上升，因而更加接近于實驗值。這也說明使用微滑移模型可為結構提供更加精確的邊界條件，使結算結果更符合實際；微滑移模型能夠在更寬的法向力范圍為結構提供振動幅值衰減，從而使動應力計算值降低。

圖13 Griffin計算結果與實驗結果對比Fig.13 Contrast between calculation result and test result of Griffin

5 結論

1)依據Yang發展的單觸點分離—接觸—滑移判據，求解了不同法向力及不同牽引切向運動時，干摩擦觸點的運動狀態及干摩擦力。

2)將Yang判據應用于條狀或環狀微滑移模型，形成了求解微滑移觸點運動狀態及受力的數值方法，通過與理論模型對比，檢測了算法的準確性。該方法克服了理論模型對于方向載荷分布及激勵條件的諸多假設，拓寬了微滑移模型的應用范圍，可為結構減振計算提供更加真實的邊界條件，同時該求解方法無須迭代，從而避免了有限元軟件求解含摩擦接觸的收斂問題。

3)用MHBM法分別求解了宏滑移及微滑移阻尼器約束的結構動響應。結果表明，在法向力較小條件下，兩者相差不大，隨著法向力的增加，兩者減振效果出現明顯差別，微滑移模型對法向力的適應能力更強，可在較寬法向載荷范圍內達到結構減振的目的。

References)

［1］晏礪堂，朱梓根，李其漢，等.高速旋轉機械振動［M］.北京:國防工業出版社，1994.YAN Litang，ZHU Zigen，LI Qihan，et al.Mechanical vibration with high rotational speed［M］.Beijing:National Defense Industry Press，1994.(in Chinese)

［2］單穎春，郝燕平，朱梓根，等.干摩擦阻尼塊在葉片減振方面的應用與發展［J］.航空動力學報，2001，16(3): 218-223.SHAN Yingchun，HAO Yanping，ZHU Zigen，et al.Application and development of platform friction damper for depressing resonant vibration of blades［J］.Journal of Aerospace Power，2001，16(3):218-223.(in Chinese)

［3］Menq C H，Bielak J，Griffin JH.The influence ofmicroslip on vibratory response，part I:a new microslip model［J］.Journal of Sound and Vibration，1986，107(2):279-293.

［4］Menq C H，Bielak J，Griffin JH.The influence ofmicroslip on vibratory response，part II:a comparison with experimental results［J］.Journal of Sound and Vibration，1986，107(2): 295-307.

［5］徐自力，常東鋒，上官博.微滑移離散模型及在干摩擦阻尼葉片振動分析中的應用［J］.機械科學與技術，2007，26(10):1304-1307.XU Zili，CHANG Dongfeng，SHANGGUAN Bo.One-bar microslip discrete model and its application to vibration analysis of blade with dry friction damper［J］.Mechanical Science and Technology for Aerospace Engineering，2007，26(10):1304-1307.(in Chinese)

［6］Csaba G.Forced response analysis in time and frequency domains of a tuned bladed disk with friction dampers［J］.Journal of Sound and Vibration，1998，214(3):395-412.

［7］馬曉秋，王親猛，張錦，等.帶干摩擦阻尼結構葉/盤系統動力學分析［J］.航空動力學報，2002，17(1):110-114.MA Xiaoqiu，WANG Qinmeng，ZHANG Jin，et al.Dynamic analysis of bladed disc systems with vibrational dampers［J］.Journal of Aerospace Power，2002，17(1):110-114.(in Chinese)

［8］Yang B D，Chu M L，Menq C H.Stick-slip-separation analysis and non-linear stiffness and damping characterization of friction contacts having variable normal load［J］.Journal of Sound and Vibration，1998，210(4):461-481.

［9］Sanliturk K Y，Ewins D J.Modeling two-dimensional friction contact and its application using harmonic balance method［J］.Journal of Sound and Vibration，1996，193(2): 511-523

［10］單穎春，朱梓根，劉獻棟.凸肩結構對葉片的干摩擦減振研究-規律分析［J］.航空動力學報，2006，21(1):174-180.SHAN Yingchun，ZHU Zigen，LIU Xiandong.Investigation of the vibration control by frictional constraints between blade shrouds-rule analysis［J］.Journal of Aerospace Power，2006，21(1):174-180.(in Chinese)

［11］Mindlin R D.Compliance of elastic bodies in contact［J］.Journal of Applied Mechanics，1949，16:259-268.

［12］Wang JH，Chen W K.Investigation of vibration of a blade with friction damper by HBM［J］.Journal of Engineering for Gas Turbines and Power，1993，115(2):295-299.

［13］Griffin J H.Friction damping of resonant stresses in gas turbine engine airfoils［J］.Transactions of the American Society of Mechanical Engineers，Journal of Engineering for Power，1980，102:329-333.

A numericalmodel for the bar or circular m icro-slip dry friction damper

ZHAO Ning，LIN Yanhu
(College of Mechanical，Northwestern Polytechnical University，Xi’an 710072，China)

A numericalmethod to solve the contact kinematic problems of elastic bar or circularmicro-slip dry friction damper was proposed.Firstly，the damper and the time course of external excitation were discretized in the space and time domain respectively，and the same number of contact points was attached to the discrete damper.Secondly，the determination criterion of the contact kinematic state of the contact point pairswas used to each contact pair，and the stiffnessmatrix of the damper wasmodified.Thirdly，the whole balance equation with the new stiffnessmatrix was solved in each time step.Differing from the way the finite element software dealing with the contact problem with friction，the iteration was avoided，thus the feasibility of the solution was guaranteed.Meanwhile，the restriction of the distribution and the time variability of the normal load was overcome，and therefore it provided a more accurate boundary condition for solving dynamic response of structure.Lastly，the response of structures constrained by themacro-slip and micro-slip damper were computed by the MHBM(Multi Harmonic Balance Method)respectively and the distinctness of the solution was analyzed.Results show that themicro-slip model damper can have a broader normal force range for vibration reduction.

dry friction damper；micro-slip；contact kinematic；multiharmonic balancemethod；dynamic response

TH133

A

1001-2486(2015)05-067-08

10.11887/j.cn.201505011

http://journal.nudt.edu.cn

2014-10-15

國家自然科學基金資助項目(51075328)

趙寧(1958—)，男，廣東廣州人，教授，博士，博士生導師，E-mail:zhaon@nwpu.edu.cn