基于分形理論的結合面法向接觸剛度模型

陳虹旭，董冠華，殷 勤，譚 峰，殷國富

(四川大學 制造科學與工程學院，成都 610000)

動態特性是制約高端機床發展的重要因素，目前對機床動態特性的研究主要采用實驗和有限元相結合的方法[1]。研究表明，機床結合面提供了機床60%～80%的柔度特性，如何對結合面準確的建模是建立機床有限元模型的關鍵問題[2]。

目前結合面法向接觸剛度模型大多基于M-B分形理論，該理論認為微凸體變形前的頂端曲率半徑R是一個隨微凸體變形量δ變化的值[3-6]，這導致對變形量δ求導時不能把曲率半徑R視為常數[7]。然而，從微凸體變形前的頂端曲率半徑的實際含義出發，該參數應僅受微凸體尺度的影響，不應該與微凸體的接觸變形量有關。之所以出現該矛盾，是因為M-B模型在推導過程中忽略了微凸體的壓縮過程。Morag等[8]對該問題提出了一個修正的分形理論，遺憾的是在他們在推導過程中存在錯誤，同時也未能給出最終的力學表達式。

本文基于修正的分形理論模型，推導出了結合面法向接觸剛度模型。該模型考慮了微凸體的接觸變形過程，解釋了微凸體變形前的頂端曲率半徑R與接觸變形量的關系。基于該模型，計算出了結合面法向接觸剛度值并錄入有限元模型進行仿真，進而與實驗結果進行對比。對比結果表明，仿真結果與實驗結果基本一致，該模型可有效的進行結合面法向接觸剛度值計算。

1 粗糙表面接觸模型

1.1 W-M函數

對于具有連續、自仿射、不可微的分形特征的表面輪廓，可用W-M函數描述[9]

(1)

式中：z(x)為粗糙表面輪廓的高度；D為粗糙表面輪廓的分形維度，1

W-M函數由無數個余弦函數組成，設截止長度為lts=γ-(nh+1)，nh為與最高頻率γnh對應的系數，波長小于lts的分形細節將被忽略。

1.2 微凸體接觸變形模型

如圖1所示，用k(k=0,1,2,…,nh)定義與剛性平面接觸的微凸體的尺度，則該尺度下的微凸體名義接觸長度ltk(lsk-1

(2)

式中：int[…]為向下取整；ltk為k尺度下的微凸體名義接觸長度，lsk-1

(3)

微凸體變形前的頂端曲率半徑可由式(4)求得

(4)

微凸體高度為δk=GD-1(γk-nh)2-D，微凸體的接近量wk為

(5)

圖1 k尺度下的微凸體Fig.1 Asperity of scale k

該模型將k尺度下的微凸體高度δk和微凸體接觸變形量wk進行區分，微凸體變形前的頂端曲率半徑不再是接觸變形量的函數。

1.3 微凸體臨界接觸變形量

第k尺度臨界接觸變形量可表達為[10]

(6)

當wk

(7)

(8)

2 結合面法向接觸剛度模型

2.1 微凸體彈性變形模型的誤差補償

在彈性變形區(ke≤k≤nh-nl)內，用剛性平面所在的那一尺度微凸體進行受力計算。由于接觸變形是剛性平面從較小尺度微凸體壓縮至較大尺度微凸體的一個連續過程，如圖2所示。該過程不同于赫茲接觸過程中的單個微凸體變形。為了使第k尺度微凸體和第k+1尺度微凸體在兩尺度交界位置的計算值相等，進行如下分析計算

圖2 誤差補償示意圖Fig.2 Error compensation

(9)

式中：hk+1為剛性平面進入k+1尺度時第k+1尺度的預壓縮量。

為便于計算，代入常用值γ=1.5(同理也可做一般性推論)，則有

(10)

根據赫茲接觸理論[11]，兩尺度交界處受力的差值為

pkm-pk+1h=pkm·Δ

(11)

式中：pkm為第k尺度微凸體處于最大壓縮量時的受力；pk+1h為第k+1尺度微凸體壓縮量為hk+1時的受力；Δ為一個常數，可由式(12)求得

(12)

由式(11)可知，彈性變形從某一尺度進入后一尺度時，用該模型計算會丟失一部分力。考慮通過誤差補償進行修正，若k0尺度的受力遠小于所關心的尺度k0+t(t=1,2,3,…)，從k0尺度開始修正，有

(13)

由式(13)可求得k0+t尺度修正后的彈力為

pk0+1m·Δt-1+pk0m·Δt

(14)

考慮到0.204 5<Δ<0.469 7，近似有

(15)

(16)

式中：pk0+th為k0+t(t=1,2,3,…)尺度修正前微凸體壓縮量為hk0+t時的受力。

(17)

(18)

(19)

因此，式(17)還可寫為

(20)

2.2 結合面法向載荷

采用M-B模型中的面積分布函數[12-14]，假設面積大于a′的接觸點總數為N，N滿足

(21)

(22)

(23)

粗糙表面塑性變形總載荷

(24)

式中：Pp為粗糙表面的塑性變形總載荷。粗糙表面處于彈性變形區(ke≤k≤nh-nl)第k尺度的微凸體受力

(25)

粗糙表面彈性變形總載荷為

(26)

(27)

(28)

粗糙表面總載荷為

P=Pp+Pe

(29)

2.3 結合面法向接觸剛度

根據球體間赫茲接觸推導得出的彈性接觸的面積公式，單個微凸體發生彈性變形的真實接觸面積為

(30)

采用和式(26)相似過程對式(30)積分，可得結合面彈性變形的真實接觸面積為

(31)

式中：f3，f4分別由式(32)、式(33)求得

(32)

(33)

式(17)兩端對接觸變形量wk求導，可得單個微凸體的法向接觸剛度為

(34)

考慮整個名義面積內所有微凸體的接觸剛度，采用和式(26)相似過程對式(34)進行積分，可得結合面的法向接觸剛度表達式為

(35)

式中：f5,f6分別由式(36)、式(37)求得

(36)

(37)

3 法向接觸剛度模型實驗驗證

3.1 分形參數的測定

實驗結合面材料為45號鋼，表面經磨削加工，彈性模量為210 GPa，表面硬度為300 MPa，采用M8×35的螺栓連接，各螺栓連接力矩為20 N·m，其示意圖如圖3所示。

圖3 結合面示意圖Fig.3 Sketch of joint surface

加工過程中，使用一塊300 mm長的毛坯進行磨削加工，之后采用線切割切下10 mm長的表面制作檢測樣品，保證樣品和結合面的一致性。采用布魯克原子力顯微鏡對樣品表面情況進行檢測，檢測結果如圖4所示。

圖4 表面檢測結果Fig.4 Surface test results

采用結構函數法提取其分形參數，結構函數可寫為

S(τ)=CG2(D-1)τ(4-2D)

(38)

式中：C可由式(39)求出

(39)

在雙對數坐標系下，結構函數S(τ)和τ呈直線關系。當斜率kτ滿足0

D=(4-kτ)/2

(40)

再由直線在S(τ)軸上的截距可計算出G值。實測輪廓數據為離散值，記采樣間隔為Δt，共采樣N組數據，則可用nΔt(n=0,1,…)代替τ，離散的結構函數為

(41)

受限于實驗設備的測量精度和數據的舍入誤差，結構函數并非絕對線性，因此可采用最小二乘法擬合一條直線。將粗糙表面輪廓數據導入Matlab中進行處理，得到結構函數和最小二乘擬合直線，如圖5所示。對10 μm×10 μm和20 μm×20 μm檢測表面分別取3組輪廓數據，計算各組數據的分形維度和特征尺度參數，并按式(29)、式(35)計算法向接觸剛度，計算結果見表1。

圖5 結構函數與擬合直線Fig.5 Structure function and fitting line

表1 分形參數及法向接觸剛度計算值Tab.1 The calculated values of fractal parameters and normal contact stiffness

3.2 模態實驗

為研究螺栓結合面的動力學特性，對比理論建模與實驗數據之間的差距，本文開展了基于m+p噪聲振動測試系統的模態實驗，實驗中結合面建模如圖6所示。

螺栓結合面采用兩柔性繩懸吊，模擬自由邊界狀態，以移動力錘法對各測點的法向和切向進行激振，在測點1處設置三向傳感器進行數據采集。測試范圍為3 000 Hz，測得的數據通過以太網上傳至計算機，經SO Analyzer后處理，得到結合面前8階模態。模態判據MAC圖如圖7所示。可見實驗結果中各階模態已充分解耦，具有很高的置信度。

圖6 結合面幾何建模Fig.6 Geometric modeling of joint surface

圖7 模態判據MAC圖Fig.7 Diagram form of modal criterion(MAC)

3.3 法向接觸剛度與分形維度的關系

取表1中10 μm×10 μm檢測表面中3組法向接觸剛度平均值作為第一組結合面法向接觸剛度K1，可得K1=2.27×107N/m；取20 μm×20 μm檢測表面中3組法向接觸剛度平均值作為第二組結合面法向接觸剛度K2，可得K2=7.93×106N/m。

之所以分別計算兩組法向接觸剛度值，是因為在測量不同尺度的檢測表面時，受檢測儀器精度的影響，測得的分形維度值存在一定的誤差。注意到特征尺度參數處于同一量級，而計算出的法向接觸剛度值有較大差異，為探究法向接觸剛度與分形維度之間的關系，取結合面法向載荷為50 000 N，特征尺度參數G為5×10-11，D=1.25,1.3,…,1.9。以分形維度D作為橫坐標，法向接觸剛度的對數lgK作為縱坐標，繪制關系曲線如圖8所示。

圖8 法向接觸剛度和分形維度的關系Fig.8 The relation between normal contact stiffness and fractal dimension

由圖8可知，當分形維度D<1.5時，分形維度每提高0.05，法向接觸剛度就增加一個數量級，可見在一定范圍內分形維度對法向剛度的影響是很大的。當分形維度D>1.5時，隨著分形維度的增高，法向接觸剛度的變化率逐漸減小。為進一步說明，定義法向接觸剛度對分形維度的靈敏度S(后面簡稱靈敏度S)，記

(42)

式中：KD0為分形維度為D0時所求得的法向接觸剛度；KD0+ΔD為分形維度為D0+ΔD時所求得的法向接觸剛度。

S表征了在分形維度D0附近出現ΔD的變化時，法向剛度計算值對分形維度變化的敏感程度，S越接近于1，計算出的法向接觸剛度對分形維度的變化就越敏感。取D0=1.25,1.3,…,1.85，ΔD為0.05，其余參數與圖8相同，繪制S與D0的關系如圖9所示。

圖9 靈敏度S與分形維度D0的關系Fig.9 The relation between sensitivity S and fractal dimension D0

由圖9可知，當分形維度小于1.5時，靈敏度S>0.8，隨分形維度的增加變化緩慢，法向接觸剛度受分形維度變化的影響較大；當分形維度大于1.5時，靈敏度S迅速下降，法向接觸剛度受分形維度的影響減小。

3.4 理論模態與實驗模態對比

結合面切向剛度采用參數識別的方法獲取[16]，并通過有限元仿真與實驗結果的對比來保證其準確性，所得剛度值為Kt=7.5×109N/m。根據所求得的剛度值，基于吉村允孝積分法的思想，在ANSYS中采用COMBIN14單元錄入結合面剛度信息，建立有限元仿真模型。仿真結果表明，前8階理論振型與實驗振型一致，限于篇幅僅列出第二組數據固有頻率誤差最大的第4階振型對比和誤差最小的第6階振型對比，如圖10所示。理論固有頻率與實驗固有頻率對比見表2。

表2 理論固有頻率與實驗固有頻率對比Tab.2 Comparison between theoretical natural frequency and experimental natural frequency

圖10 理論振型與實驗振型對比Fig.10 Comparison between theoretical modal shapes and experimental modal shapes

由表2可知，固有頻率計算值均小于實驗值，說明所求得的法向剛度值偏小。其中第4階為結合面法向二階彎曲，誤差最大，第6階為結合面切向一階彎曲，誤差最小。所求得的法向剛度值偏小的原因是因為該模型在計算受力時考慮了塑性變形，而在計算剛度時忽略了塑性變形和彈塑性變形的影響，因此計算出的剛度值是偏小的。

4 結 論

(1)建立了結合面法向接觸剛度模型并進行了實驗驗證，結果表明理論模型和實驗模型前8階模態的振型一致，固有頻率的誤差為-10.2 %～-1%，間接論證了模型的正確性。

(2)本文模型在計算受力時考慮了塑性變形，而在計算剛度時忽略了塑性變形和彈塑性變形的影響，因此計算出的結合面法向接觸剛度值是偏小的。

(3)當分形維度小于1.5時，法向接觸剛度對分形維度的靈敏度S>0.8且變化緩慢，法向接觸剛度受分形維度的影響較大；當分形維度大于1.5時，靈敏度S隨分形維度的增大而快速減小，法向接觸剛度受分形維度的影響減小。因此，當結合面的分形維度小于1.5時，提高分形維度可以有效提升結合面的法向接觸剛度。