基于Iwan 模型的接合面切向響應建模1)

占旺龍 李 衛 黃 平

*(深圳技術大學中德智能制造學院,深圳 518118)

?(華南理工大學機械與汽車工程學院,廣州 510640)

引言

機械工程結構中廣泛存在著以螺紋連接、法蘭盤和鉚接為代表的搭接接頭形式,如圖1 所示.搭接連接的典型特征是接合面在法向預緊力作用下發生接觸,由于表面粗糙度的存在使得接觸界面上承受著非均勻的壓力分布.當連接件受到切向力作用時,接觸面間的摩擦力會阻礙界面間的宏觀相對運動.假設接觸界面上摩擦系數處處相等,則顯然接觸界面上各處抵抗切向相對運動的能力各不相同,法向接觸壓力大的區域抵抗切向相對運動的能力強,而法向接觸壓力小的區域則抵抗切向相對運動的能力弱.由于界面接觸壓力大小與粗糙峰高度正相關,因此在切向力作用下微凸體高度小的先發生滑移,而高度大的粗糙峰后發生滑移,當切向力增大到一定程度時,滑移區將覆蓋整個接觸界面從而導致連接失效.在往復振蕩切向載荷作用下,接觸界面上的黏著區與滑移區會交替變化,力與位移關系在笛卡爾坐標系中形成遲滯曲線,遲滯曲線所圍成的面積即為振蕩加載過程中的能量耗散[1-2].

圖1 螺栓搭接接頭Fig.1 Bolted lap joint

在傳統的機械動力學研究中,一般都采用等效的彈簧阻尼元件來模擬連接界面間的接觸力學行為而忽略接合面的非線性特征[3],模型中的參數需要通過實驗進行標定.然而,由于連接界面間的非線性特征,在某種條件下標定得到的模型參數并不一定適用于其他實驗條件[4].有限元方法也被廣泛應用于簡單裝配結構的建模分析工作中,但是若想要準確地了解接合面的滑移行為則需要更加細密的網格,這將大大增大計算時間.與此同時,對于由隨機粗糙面組成的接合面,要建立精確的有限元模型還存在很大的困難,因此開展搭接接頭非線性特性研究并建立基于實際可測物理參數的理論模型具有重要的意義[5-7].

目前用于描述接合面非線性力學行為的代表模型有:Iwan 模型[8]、LuGre 模型[9]和Valanis 模型[10]等.但是這些模型大部分都是唯象的,模型中的參數靠實驗標定而往往缺少明確的物理含義,因此也就無法反應出接合面粗糙程度對宏觀力學行為的影響.Iwan 模型[11]最早用來描述土木結構的彈塑性遲滯行為,而后Segalman[8]通過定義連接界面臨界滑移力密度分布函數,將該模型推廣至連接界面黏滑行為建模.在小幅振蕩周期性載荷作用下,單個加載周期內的能量耗散值D與切應力幅值F0之間的關系可以表示為

式中,ν 和χ 通過實驗數據擬合得到,χ 范圍在-1 到0 之間[12].利用該關系,Segalman 定義的臨界滑移力分布函數為

式中,H(·)為Heaviside 單位階躍函數;φmax為發生宏觀滑移時的位移;S為界面開始發生宏觀滑移時力-位移曲線斜率不連續程度,R和S的具體計算方法見文獻[8].

從式(2)可以看出,當φ 增大時,ρ(φ)會逐漸減小;當φ=φmax時,ρ(φ)會發生突變;而當φ ＞φmax時,ρ(φ)=0,其對應的圖像形式如圖2(a)所示,其為含截斷冪律分布和單脈沖的四參數非均勻密度函數[8].Song 等[13]和張相盟等[14]基于均勻密度函數建立起摩擦接頭Iwan 模型(如圖2(b)所示).其他一些比較常見的臨界滑移力分布函數如圖2(c)和圖2(d)所示,具體解釋可參閱文獻[15].而后Argatov 等[16]、王東等[17]利用Kragelsky-Demkin 粗糙面接觸理論求解了螺栓連接接合面間的非線性力學行為.Brake[15]還將Iwan 模型擴展到連接界面發生宏觀滑移后的階段.另外李一堃等[18-19]提出一種含截斷冪律分布和雙脈沖的非均勻密度函數的六參數Iwan 模型.

圖2 臨界滑移力密度分布函數圖[20]Fig.2 Density distribution function of critical slip force[20]

本文針對前人研究中臨界滑移力分布函數需要進行參數標定的缺陷,從接觸面統計接觸模型出發,推導出新的臨界滑移力分布函數,在此基礎上對典型工程摩擦副表面進行數值計算并對結果進行分析,最后還將推導出的函數用于與已發表的實驗結果進行對比,結果表明理論推導出的分布函數與實驗結果具有很好的一致性,該研究成果可為搭接接頭的切向響應研究提供參考.

1 Iwan 模型

組成Iwan 模型的兩種基本單元如圖3 所示,每個基本單元由一個線性彈簧元件和一個庫侖摩阻片按照串聯或者并聯方式組成.并聯單元允許在位移不發生改變的情況下改變力的大小;而串聯模型允許位移發生改變而力不發生變化.

一般可以將隨機粗糙面的接觸問題簡化成剛性平面與球體高度隨機分布的彈性體之間的接觸,如圖4(a)所示.由于球體高度隨機分布,那么每個球體所受到的法向力則會是球體高度的函數.根據庫侖摩擦定律,即摩擦副間的摩擦力與法向力成正比,這樣由于每個微凸體所承受的法向載荷不同而使得其所能承受的最大切應力也各不相同.與Iwan 模型對比發現[16],可以將摩擦面間的切向接觸問題用串-并聯單元來描述,如圖4(b)所示.該模型由若干個Jenkins 單元并聯組成,而每個Jenkins 單元又是由彈性常數為k的彈簧和閾值為qi的庫侖摩阻片串聯而成.通過類比可以將圖4(a)中每個微凸體的接觸行為用圖4(b)中的一個Jenkins 單元進行描述.

圖3 Iwan 模型基本單元Fig.3 Basic element of Iwan model

圖4 (a)粗糙面接觸模型;(b)簡化并-串聯Iwan 模型Fig.4 (a)Rough surface contact model;(b)Simplified parallel-series Iwan model

Iwan 模型的工作原理如下:當單個Jenkins 單元所受到的切向力Fi小于摩阻片滑動閾值qi時,彈簧將會發生彈性形變從而承受切向力作用;而當切向力Fi大于摩阻片滑動閾值qi時,摩阻片將會發生滑動,此時該單元所能承受的最大切向力為qi.假設在力F作用下,上表面相對于下表面切向相對位移為x,根據上述描述的Jenkins 切向滑動行為,可以用數學關系表示為

由于微凸體高度隨機分布,根據類比可以假定摩阻片滑動力閾值也隨機分布.假設其摩阻片臨界滑移力密度分布函數為ρ(q),這樣當＞0 時,整個體系切向力與相對位移之間的關系可以表示為

式中,ρ(q)dq為滑動力閾值在(q,q+dq)之間的Jenkins 單元數.式(4)中等號右邊第一項為處于滑動狀態下Jenkins 單元所產生的總切向力,此時所對應的微凸體發生滑動;第二項為處于彈性變形狀態下Jenkins 單元所產生的總切向力,此時微凸體處于黏著狀態.若式(4)等號右端第二項為零,此時界面發生宏觀滑移;而若該項不為零,此時接觸界面部分微凸體處于滑移狀態.因此摩擦副界面間的運動可以分為兩個階段:

(1)微滑移階段:當外界施加的切向載荷較小時,接觸面承受法向載荷較小的局部會發生滑動,而其他部分處于黏著狀態;

(2)宏觀滑移階段:當外界施加的切向載荷大到一定值的時候,接觸面間會發生明顯的宏觀滑動,此時整個接觸面都發生相對滑動.

利用以下變量關系式對式(4)進行變換

式中,σ 為微凸體高度分布方均根.經過上述變量代換后,式(4)中變量的量綱將發生變化,由于σ 量綱為長度,因此和φ 為無量綱量,同理的量綱為力.代入式(4)后可得

對上式求一階導數,可以得到摩擦副間的切向剛度kt為

對式(6)求二階導數可得

對于一般的摩擦副而言,其上下表面間能承受的最大切向力為摩擦力,即Fmax=μW.從上式可以看出變量φ 的期望為滑動摩擦力,并且當φ →∞時,有若將式(9)兩端同時除以摩擦力可得

從上述分析可知,若要將Iwan 模型成功地應用于粗糙表面之間的切向接觸問題研究,首先要解決的問題應該是推導出正確的分布函數該分布函數要滿足恒大于或等于零及式(7)～式(9)的條件.

2 周期性加載力-位移曲線

將求得的臨界滑移力密度分布函數代入式(6)便可得到切向力與切向位移之間的關系.一般來說,分布函數為分段函數,因此在計算加載曲線時要考慮φ 成立的區間范圍,例如,當時,切向力計算式為

對于其他區間范圍也同樣可以用上述方法求出.由于式(11)形式較為復雜,特別是當表面粗糙峰處于非高斯分布時,為得到切向力與切向位移之間的關系,可以用Simpson 數值積分方法進行求解.

當摩擦副承受位移幅值為a的周期性位移激勵時,在位移-力圖上可以得到如圖5 所示的遲滯曲線.其中脊線函數可以根據式(6)得出,而正、反向加載曲線可以通過Masing 模型[21]求出.當體系由圖5 所示的F(·)的(a,F(a))位置反向加載到(-a,F(-a))位置時,反向加載力-位移曲線方程為

圖5 力-位移遲滯曲線Fig.5 Force-displacement hysteresis curve

將式(12)和式(13)代入式(14),即可得到不同位移幅值下單個加載周期內的能量耗散值.

3 臨界滑移力分布函數推導

如圖4(a)所示,粗糙面之間的接觸問題可以等價成剛性平面與覆蓋球體但高度隨機分布的彈性面之間的接觸.一般用輪廓儀測量得到的是σ,但在現有的統計學模型中用到的是σa,兩者存在如下關系[22]

其中,β=σRη,R為微凸體曲率半徑,η 為微凸體分布密度,對一般工程表面來說,β 范圍為0.02～0.06.當微凸體高度服從正態分布時,分布函數可以寫成

對z進行無量綱化,令z*=z/σ,則分布函數可以寫成

塑性指數ψ 為[23]

式中,E為綜合彈性模量,,E1,2為摩擦副材料的彈性模量,ν1,2為摩擦副材料的泊松比,H為較軟材料硬度值,K最大應力接觸系數.

根據粗糙表面ZMC 彈塑性接觸模型[23],當微凸體法向干涉量ω ＜ω1(ω=z-d,d為摩擦面間的法向接觸間距)時,微凸體發生完全彈性變形,此時法向載荷We與微凸體法向形變量ω 的關系為

當ω ＞ω2時,微凸體發生完全塑性變形,此時法向載荷Wp與微凸體法向形變量ω 的關系為

當ω ∈[ω1,ω2]時,微凸體間既有彈性變形,又存在塑性變形,利用樣板函數插值方法可到法向載荷Wep與微凸體法向形變量ω 的關系為

式中,ka=1-2K/3,其中K為最大接觸應力系數,K=0.454+0.41ν(ν 為材料泊松比).ω1為微凸體開始發生彈塑性變形時的法向臨界干涉量,ω1=(πKH/E)2R;ω2為開始發生完全塑性變形時的法向干涉量閾值,ω2=54ω1[24].

利用統計學方法對式(19)～式(21)進行積分可以得到粗糙表面接觸時的法向載荷與變形量的關系為

式中,φ(z)為微凸體高度分布函數,An為名義接觸面積,具體推導過程可以參見文獻[23-24].由于法向載荷W與法向接觸距離d之間為單調關系,即接觸間距越小,接觸面間的法向載荷越大.因此在給定法向載荷的條件下,可以通過二分法計算得到給定條件下的法向接觸間距d.

根據庫侖摩擦定律,摩擦界面間的摩擦力為

式中所有帶星號(*)的變量都是用高度分布均方根σ 進行無量綱化后的變量,μ為滑動摩擦系數.對比式(23)與式(9),由于積分計算與積分變量無關,因此將ω*替換為φ 便可得到分布函數,為

4 實驗結果驗證

Eriten 等[25]對不同預緊力作用下螺栓搭接接頭的切向行為進行了實驗研究,組成摩擦副的材料為420 型不銹鋼,接觸表面經銑削加工而成.材料參數為:彈性模量E=200 GPa,泊松比ν=0.24,硬度H=5.825 GPa,名義接觸面積An=156 mm2,表面粗糙度參數σ/R=0.088 8,β=0.023,σ=2.677 μm,ψ=4.725,摩擦系數μ=0.5,具體實驗裝置及過程可以參考文獻[25].需要注意的是,在文獻[25]中摩擦系數為一與法向載荷有關的參數,但在本文中設定摩擦系數為發生宏觀滑動時摩擦力與法向載荷的比值,認為其為一恒定量,本文取摩擦副發生宏觀滑動時摩擦系數為計算值,約為0.5.

圖6 給出了不同法向載荷下對摩擦副進行切向循環位移加載得到的遲滯曲線圖.其中實線為用本文中模型計算得到的遲滯曲線,點狀線的為文獻[25]中實驗得到的曲線.從圖中可以看出理論計算值與實驗結果很吻合,從而證明了本文模型的可行性.

圖7 給出了不同法向載荷下單個加載周期內的能量耗散實驗值與理論計算結果圖,其中實驗值是統計10 個加載周期內得到的統計箱式圖,黑色菱形為計算結果.從圖中可以看出計算結果與實驗結構很吻合,進一步證明了本文提出模型的有效性.

5 數值計算結果與討論

下面以鋼-鋼接觸為例,計算中所用到的材料物理參數為[23]:彈性模量E1=E2=207 GPa,泊松比ν1=ν2=0.29,硬度H1=H2=1.96 GPa,名義接觸面積An=100 mm2,摩擦系數μ=0.5.由式(24)可知滑移力密度分布函數與形貌參數σ和R都有關,根據McCool[22]的研究,曲率半徑R、粗糙峰分布密度η 及方均根σ 與測量的表面形貌高度z(x)之間為耦合關系.因此對于一般的工程應用表面,形貌參數可以用兩個參數描述,即:β和σ/R.而這兩個參數直接決定塑性指數ψ 的大小,所以下面僅討論不同塑性指數ψ 對本構關系的影響.表1 給出典型工程表面的形貌參數及對應的塑性指數值.

圖6 不同法向載荷及位移幅值下的遲滯曲線實驗與數值結果對比Fig.6 Experimental and numerical results comparision of hysteresis curves under different normal loads and displacement amplitudes

圖7 單個周期內的能量耗散值實驗與數值結果對比(箱圖為實驗值,黑色菱形為計算值)Fig.7 Experimental and numerical results comparision of energy dissipation per cycle(the box shows the experimental value,and the black diamond is the numerical value)

表1 典型工程表面參數及塑性指數Table 1 Typical engineering surface parameters and plasticity index

假設微凸體高度滿足正態分布,根據上述推導公式得出如圖8 所示的在不同塑性指數和不同法向載荷下的摩阻片滑動閾值分布函數和切向剛度圖.其中圖8(a)和圖8(b)分別為塑性指數為0.7 和2.5 時的臨界滑移力分布函數圖.由式(5)中第二式可知,對于給定的粗糙表面,φ 是摩阻片滑動閾值q的度量值,φ 越大即摩阻片滑動閾值q越大.另外由式(10)可知,的數學期望為1,即圖8(a)和圖8(b)圖所示分布曲線的期望值為1.對比圖8(a)和圖8(b)圖中不同法向載荷下的分布函數可以發現,隨著法向載荷的減小,出現概率最大的摩阻片滑動閾值也隨之減小,這主要是因為法向載荷越小,表面間接觸間距d變小,此時單個微凸體承擔的法向載荷減小從而使其發生相對滑動所需的切向力也隨之減小,即摩阻片滑動閾值減小.對于給定的法向載荷,雖然φ 的取值范圍為[0,+∞],但當φ ＞1 時,迅速趨于零,這是因為為了保證期望為恒定值,必須在達到一定程度后迅速減小到零,法向載荷越小,衰減速度越快.從微凸體高度分布函數也可以進行解釋,因為雖然正態分布中自變量的取值范圍為[-∞,+∞],但在±3σ范圍內的微凸體高度所占比例達到99.73%,在該范圍以外的可以忽略不計,另外實際工程表面亦不存在無窮高度的微凸體,這些因素都保證了具有水平漸近線.對比相同載荷下不同塑性指數時的分布曲線可以看出,塑性指數越大,φ 的眾數越小,眾數出現的概率越大.這主要是因為隨著塑性指數的增大,高度較高的微凸體將發生塑性變形,該狀態下法向干涉量增大但法向載荷保持不變,這樣就會導致更多其他的微凸體發生接觸但實際承擔的法向載荷較小,從而導致隨著塑性指數的增大而眾數減小,同時出現的概率增大.圖8(c)和圖8(d)為不同法向載荷下的切向剛度隨切向位移的變化關系圖,由式(7)可知圖中單條曲線為在給定載荷下臨界滑移力分布函數的尾部分布積分值.從圖中可以看出隨著φ 的增大,切向剛度迅速減小到零,此時摩擦副間將發生宏觀滑動.φ=0 所對應的量即為初始切向剛度,從圖中可以看出隨著法向載荷的增大,初始切向剛度也隨之增大,這也與文獻[26]中的結論相一致.同時塑性指數越大,初始切向接觸剛度也越大.

圖8 不同法向載荷及塑性指數下的分布函數及切向剛度Fig.8 Distribution function and tangential stiffness under different normal load and plasticity index

圖9 為在塑性指數ψ=0.7 時無量綱切向力與切向位移關系圖,從圖中可以看出加載曲線為凸函數,隨著切向位移的增大,切向力也同樣隨之增大,切向力增大速率逐漸減小并趨于零,表明切向力逐漸趨于滑動摩擦力,此時摩擦副發生宏觀滑動.法向載荷越大,發生宏觀滑動所需的切向位移越大.

圖9 無量綱化切向力-位移曲線(ψ=0.7)Fig.9 Dimensionless tangential force-displacement curve(ψ=0.7)

圖10 為當塑性指數ψ=0.7 及法向載荷W=100 N 條件下對摩擦副進行周期性加載時不同位移幅值下用Masing 準則求解的力-位移遲滯曲線,由圖中可以看出不同的位移幅值下,脊線曲線都相同.隨著位移幅值的增大,遲滯曲線所圍成的面積增大,表明單個周期內能量耗散值隨位移幅值的增大而增大.

圖10 不同切向位移幅值下的遲滯曲線(ψ=0.7,W=100 N)Fig.10 Hysteresis curves under different tangential displacement amplitudes(ψ=0.7,W=100 N)

6 總結

本文對以螺栓連接為代表的搭接接頭為例,對在法向預緊力作用下的粗糙面接觸問題進行研究.基于Iwan 模型研究界面接觸的非線性力學性質,重點推導了Iwan 模型中摩阻片臨界滑移力密度分布函數.相比于前人研究,本文推導出的分布函數不再是唯象描述,而是從材料性能參數和表面粗糙度參數推導得出.在求出臨界滑移力密度分布函數的基礎上,進而求解出接觸面間的切向接觸響應和在周期性位移加載條件下的遲滯曲線及能量耗散.主要結論如下:

(1)臨界滑移力函數開始迅速上升,到達最高點后迅速收斂到零.收斂速度和最大值出現位置與法向載荷、塑性指數及粗糙度參數有關.

(2)界面切向接觸剛度是滑移力分布函數的尾部分布.剛開始時切向剛度最大,進一步增大切向位移時,切向剛度逐漸減小到零,此時界面發生宏觀滑動.初始切向剛度與法向載荷、粗糙度參數及塑性指數有關.對于確定的接觸表面,法向力越大,初始切向剛度越大;初始切向剛度同樣也隨著塑性指數的增大而增大.

(3)使用本文推導出的函數計算得到的單位加載周期內的遲滯曲線與實驗結果很吻合,計算出的單個周期內的能量耗散也與實驗結果吻合,證明模型的正確性.