王世軍, 劉鑫, 范凌松, 崔圣奇, 李鵬陽

(西安理工大學(xué)機械與精密儀器工程學(xué)院, 710048, 西安)

機械結(jié)構(gòu)中兩個相互接觸的表面稱為結(jié)合面。宏觀上是光滑平整的零件表面,微觀上總是粗糙不平的。兩個零件表面的接觸,微觀上都是兩個粗糙表面的接觸,接觸性質(zhì)并不是兩個光滑表面的接觸性質(zhì)。結(jié)合面的形貌特征和材料特性決定了結(jié)合面的接觸特性。如何從接觸面的微觀形貌入手構(gòu)建準確的接觸模型,對于研究結(jié)合面的接觸特性[1]乃至分析機械整機的靜態(tài)和動態(tài)性能,都有重要意義。

在結(jié)合面的接觸問題研究中,學(xué)者大多將結(jié)合面簡化成粗糙表面和剛性平面的接觸,相當于峰-峰接觸,此時兩微凸體的接觸角度為0°。雙粗糙表面微觀上微凸體之間除了正接觸,大多為肩-肩接觸,接觸角度范圍在0°～90°之間[2]。

針對側(cè)接觸問題,Sepehri等[3]將微凸體定義為橢球體,研究了傾斜接觸力在法線和切線方向上的分量,建立了SF統(tǒng)計模型。Gorbatikh等[4]提出了一個接觸角度分布函數(shù),并用統(tǒng)計方法研究了微凸體之間的側(cè)向接觸性質(zhì),但是沒有給出這個角度分布函數(shù)的理論依據(jù)。Misra等[5]通過引入一個方位函數(shù)來模擬接觸角和各向異性表面,利用微力學(xué)模型來研究傾斜接觸微凸體的滑動行為,并將剪切方向與法向的接觸剛度耦合起來,建立了考慮接觸方向分布和微滑移特性的側(cè)接觸統(tǒng)計模型。高志強等[6]引入文獻 [4]的接觸角度函數(shù),建立了考慮微凸體斜接觸和相互作用的結(jié)合面接觸剛度及阻尼模型。Jamshidi等[7]根據(jù)微凸體接觸時的穿透量來界定接觸狀態(tài),同時對Misra的接觸角度函數(shù)進行修正,由此建立了接觸界面法向載荷與剪切力耦合的HJ模型。文獻 [4-7]在研究斜接觸問題時雖然給出了接觸角的分布函數(shù),但對其分布規(guī)律的來源和依據(jù)沒有明確說明。范凌松等[8-9]研究了粗糙表面上相鄰微凸體的水平距離,發(fā)現(xiàn)水平距離遵循高斯分布規(guī)律,建立了包含法向和切向接觸特性的FLS模型,研究了相關(guān)輪廓統(tǒng)計參數(shù)及相互作用對結(jié)合面接觸剛度的影響。

利用分形模型研究側(cè)接觸問題時,孫見君等[10]根據(jù)微凸體斜接觸時上下兩輪廓的空間坐標與接觸角的聯(lián)系,建立了雙粗糙面的側(cè)接觸分形模型,分析了表面形貌參數(shù)、接觸壓力對孔隙度和實際接觸面積的影響,該模型采用的接觸角函數(shù)是由空間幾何關(guān)系獲得,與水平距離無關(guān)。Zhang等[11]基于MB模型和側(cè)接觸理論建立了考慮接觸角分布的ZK改進分形模型,該模型假定接觸角的范圍在0°～45°之間,通過給定接觸面積與接觸角的函數(shù)關(guān)系來研究接觸角對接觸剛度的影響,該模型忽略了水平距離的存在,且假定的角度范圍與實際不符。以上模型都未對接觸角的理論來源加以研究。

本文基于分形理論研究粗糙表面上相鄰微凸體之間的高度差和水平距離的分布規(guī)律,采用分形方法描述相鄰微凸體之間的距離分布,結(jié)合相鄰微凸體之間高度差的分形規(guī)律,遵循微凸體峰高與水平距離的統(tǒng)一性,得到粗糙表面上微凸體側(cè)向接觸時接觸角的分布函數(shù),同時考慮彈塑性邊界條件的連續(xù)性及域擴展因子對真實接觸面積的影響,以此為根據(jù)建立考慮微凸體水平距離和接觸角分布的側(cè)接觸分形模型。通過仿真分析探究了相關(guān)接觸參數(shù)對結(jié)合面接觸剛度的影響。

1 微凸體中心水平距離分布規(guī)律

針對45鋼磨削加工表面,利用DCM3D萊卡顯微鏡分別沿著輪廓水平方向按1 μm分辨率和高度方向上按0.25 μm分辨率進行采樣,根據(jù)采樣數(shù)據(jù)繪制的三維表面輪廓圖如圖1所示。

圖1 磨削表面三維輪廓

本文統(tǒng)一按“三點峰”的定義來獲取一組輪廓上n個峰值點的水平坐標xi,i取1～n,以第1個峰為基準,之后每個峰相對上一個峰的坐標差定義為水平距離si,同理把后續(xù)每個峰與前一個峰的高度差的絕對值定義為峰的高度差hi,如圖2所示。

圖2 微凸體水平和峰高距離定義

針對磨削試樣按照順紋理方向提取多組斷面的二維輪廓數(shù)據(jù),分別統(tǒng)計各組輪廓不同間隔下的水平距離數(shù)據(jù),擬合出的輪廓曲線如圖3所示。

圖3 不同間隔水平距離輪廓

已有研究表明,磨削加工得到的機械表面在微觀結(jié)構(gòu)上微凸體的接觸類型主要為側(cè)向接觸[10]。本文采用PSD法來獲取上述輪廓的分形參數(shù),首先利用MATLAB編寫程序,得到對數(shù)坐標下的功率譜密度函數(shù),利用最小二乘法線性擬合求出其分形參數(shù),從而驗證水平距離是否具有分形規(guī)律。

對文獻 [12]中的PSD函數(shù)左右各取對數(shù)

lnP(ω)=(2D-5)lnω+(2D-2)lnG-ln(2lnγ)

(1)

式中:D為分形維數(shù);G為分形尺度參數(shù);ω為空間截止頻率;γ為尺度參數(shù),通常取1.5。PSD曲線斜率為2D-5,截距為(2D-2)lnG-ln(2lnγ)。按照該方法對圖3中的水平距離輪廓數(shù)據(jù)進行計算。間隔1 μm下水平距離輪廓的PSD曲線和最小二乘法擬合的曲線如圖4所示。

圖4 1 μm間隔水平距離功率譜密度曲線

按照相同方法分析相鄰微凸體峰的高度差輪廓,峰的高度差滿足分形規(guī)律。根據(jù)三點峰兩組表面輪廓數(shù)據(jù)計算得到的分形參數(shù)如表1所示,分析表1中數(shù)據(jù)可以看出,不同采樣間隔下三點峰水平輪廓和峰的高度差輪廓的斜率k在 [-3,-1]之間[13],證明本文按三點峰定義的兩種輪廓滿足分形規(guī)律,因此可以利用二者的分形規(guī)律建立法向和切向剛度模型。

表1 輪廓分形參數(shù)統(tǒng)計

2 單對微凸體側(cè)接觸模型

微凸體在法向力P和切向力T作用下發(fā)生側(cè)接觸時,其接觸過程如圖5所示,S、R分別為微凸體的高度、曲率半徑,下標1、2代表上、下兩微凸體,下標n、τ代表法線和切線方向,結(jié)合面的平均間距用d表示,微凸體中心水平距離和接觸角度分別用r、θ表示。在Z′OX′坐標系中,法向接觸力P沿上下微凸體接觸面的法線和切線方向可分解為Pn、Pτ,同理切向力T可分解為Tn、Tτ,上下微凸體在Z′軸方向的形變量為ωZ′,X′軸方向位移量為ζX′。在ZOX坐標系中,Z軸方向變形為ωZ,X軸方向位移為ζ。

圖5 微凸體肩對肩接觸

單對微凸體在發(fā)生側(cè)接觸時存在以下幾何關(guān)系

(2)

(3)

δ=S1+S2-d-r2/2Rs

(4)

式中Rs=R1+R2[14]。

根據(jù)側(cè)接觸受力關(guān)系可得

P=PZ′cosθ-TX′sinθ

(5)

T=PZ′sinθ+TX′cosθ

(6)

式中:PZ′=TX′tan(α+θ),其中PZ′、TX′分別代表Z′、X′方向所受的合力。

2.1 彈性接觸變形階段

單對微凸體發(fā)生側(cè)接觸時,實質(zhì)上為兩個半球體在側(cè)臂上的接觸,此時兩球體存在一個接觸角θ,定義單對微凸體接觸時Z′方向上的變形量為ωZ′,由Weierstrass-Mandelbrot分形函數(shù)可得單個微凸體的變形及等效曲率半徑[15]

ωZ′=2GD-1(lnγ)0.5(2r′)2-D

(7)

R=r′D/[24-DGD-1(lnγ)0.5]

(8)

式中:r′為微凸體的截斷半徑。

根據(jù)圖中幾何關(guān)系可以得到微凸體在Z′方向上的變形

ωZ′=ωZcosθ+ζsinθ

(9)

當微凸體達到臨界變形接觸狀態(tài)時,在Z′方向上變形量為[16]

δZ′ec=(πKH/2E)2R

(10)

聯(lián)立式(7)、(8)和(10)可得微凸體在Z′方向上的臨界接觸面積

(11)

微接觸點的截面積a的面積分布函數(shù)[18]為

(12)

式中:al為最大微接觸點的截斷面積;ψ為與分形維數(shù)D有關(guān)的域擴展因子,且近似滿足ψ=1.549D-1.253+1.069。

當微凸體在Z′方向上變形量ωZ′<ωZ′ec時,此時微凸體為完全彈性接觸狀態(tài),依據(jù)赫茲接觸理論及文獻 [19],可得微凸體彈性接觸階段的法向合力,切向合力與變形量ωZ′關(guān)系[19]如下

pe=4ER1/2b3/2c1/3

(13)

te=4ER1/2b3/2c2/3

(14)

b=ωZcosθ+ζsinθ

(15)

c1=[cosθ-sinθ/tan(θ+α)]

(16)

c2=[sinθ+cosθ/tan(θ+α)]

(17)

彈性階段,兩微凸體之間的法向接觸剛度

kne=dpe/dω=2ER0.5b0.5c1cosθ

(18)

切向接觸剛度為

kτe=dte/dζ=2ER0.5b0.5c2sinθ

(19)

2.2 彈塑性變形階段

上下微凸體發(fā)生側(cè)接觸時在彈塑性階段的臨界變形范圍為ωZ′ec≤ωZ′≤110ωZ′ec。根據(jù)文獻 [20],在ZOX坐標系中,微凸體在彈塑性區(qū)法向合力、切向合力與變形量ωZ′之間的關(guān)系分別為

(20)

(21)

在ZOX坐標系中,微凸體的法向接觸剛度、切向接觸剛度可表示為

(22)

(23)

式中:m=ln(330/K)/ln(110)。

2.3 塑性變形階段

當微凸體在Z′方向上變形量ωZ′≥110ωZ′ec時,微凸體進入完全塑性接觸狀態(tài),此時微凸體的平均接觸載荷與兩接觸材料中較軟一方的硬度相等。完全塑性變形階段Z方向上的法向合力與變形量ωZ′的關(guān)系[21]如下

pp=2πHRbc1

(24)

此時法向接觸剛度

(25)

3 結(jié)合面的接觸剛度分形模型

由于水平距離和峰的高度差都具有分形規(guī)律,因此可用W-M函數(shù)來表征相鄰微凸體的輪廓

(26)

(27)

式中:Dh,s、Gh,s分別為分形維數(shù)、分形尺度參數(shù);γn為表征粗糙表面的頻譜;n為頻率指數(shù);xs為水平距離;xh為峰-峰相對距離。

在坐標系中,接觸角θ與兩個輪廓的正切值有關(guān),因此接觸角的分布函數(shù)為

θ(xs,xh)=arctan[H(xh)/S(xs)]

(28)

根據(jù)接觸角幾何關(guān)系及式(28)可得微凸體接觸角的概率密度函數(shù)

f(θ)=

(29)

式中:θ為肩對肩接觸角,范圍為0°～90°;L為接觸表面加載次數(shù)。

根據(jù)式(29)擬合出分布函數(shù)曲線,分布規(guī)律如圖6所示,可以發(fā)現(xiàn):接觸角度函數(shù)值隨著接觸角θ的增加而減小;隨著加載次數(shù)增加,接觸面逐漸壓平,接觸角度函數(shù)值變化逐漸緩慢。

圖6 側(cè)接觸角度分布規(guī)律

微凸體發(fā)生彈性、彈塑性和全塑性變形的3個區(qū)域面積之和構(gòu)成了結(jié)合面的真實接觸面積[20]

(30)

將式(12)代入式(30),可得

amax=Ar(2-D)Ar/[Dψ(1-0.5D)]

(31)

根據(jù)文獻 [21]及式(7)、(8)、(10)、(12)和(31),結(jié)合面無量綱法向接觸載荷P*可表示為

((2-D)/D)0.5Df(θ)dadθ

(32)

無量綱切向接觸載荷T*為

a[(0.5D-1)+(1-D)m]c2f(θ)dadθ

(33)

(34)

(35)

4 結(jié)果與分析

(a)D對的影響

對上述仿真曲線進行分析可得如下結(jié)果。

圖8 材料參數(shù)對切向接觸剛度的影響

(3)接觸參數(shù)對真實接觸面積的影響。真實接觸面積隨法向載荷的變化曲線如圖9所示,可知D、塑性指數(shù)φ的增大會使真實接觸面積不斷變大,D越大,表面輪廓越趨近于光滑,而φ越大微凸體越不易發(fā)生塑性變形,從而使臨界的接觸面積不斷變小,而真實接觸面積的占比不斷變大。因此,改善分形維數(shù)和塑性指數(shù)有助于提高真實接觸面積。

圖9 接觸參數(shù)對真實接觸面積的影響

(4)法向載荷系數(shù)的影響。圖10為結(jié)合面接觸剛度與法向載荷系數(shù)的變化關(guān)系。h分別取2～50時,無量綱法向和切向接觸剛度均隨著法向載荷系數(shù)h的增大而增大,隨著分形尺度參數(shù)G*的增大而減小。固定分形尺度參數(shù)G*時,保持結(jié)合面法向載荷P不變,h的增大意味著結(jié)合面切向接觸載荷T減小,則接觸面抵抗切向變形的能力增強,從而結(jié)合面法向和切向接觸剛度都會增大。

與h的變化曲線

(5)接觸方式的影響。圖11為兩種模型的接觸剛度和真實面積隨法向載荷變化的曲線。其中LZT模型[16]的接觸類型為正接觸,本文模型則考慮了側(cè)接觸,給定相同的無量綱法向載荷,新模型的接觸剛度和真實接觸面積小于LZT模型,根據(jù)圖6及文獻 [6]的研究可知,正接觸的接觸角度函數(shù)值相比側(cè)接觸較大,而LZT模型忽略了側(cè)接觸的影響,導(dǎo)致該模型的剛度偏高。

(a)真實接觸面積

5 有限元仿真與試驗驗證

采用文獻 [9]的試驗裝置來驗證本文提出的側(cè)接觸分形模型,結(jié)合部由兩半環(huán)及螺栓構(gòu)成,半環(huán)材料均為45鋼,半環(huán)接觸面為磨削表面,兩材料的厚度和半徑分別為35、275 mm,表面經(jīng)丙酮清洗,其他材料參數(shù)如表2所示。

表2 整環(huán)材料參數(shù)

截取部分二維輪廓數(shù)據(jù),經(jīng)過PSD功率譜計算可得D=1.425,G=2.928×10-10m,半環(huán)聯(lián)結(jié)面上每個螺栓能承受的最大擰緊力矩為80 N·m。本文采用錘擊法獲得整環(huán)試驗裝置的振型和固有頻率,同時建立整機的有限元接觸模型并進行模態(tài)分析,綜合比較試驗與有限元模型的結(jié)果。

由于包含結(jié)合部的機械結(jié)構(gòu)的接觸剛度具有非線性的特征,因此需要將包含接觸的非線性模型線性化,再通過有限元模態(tài)分析進行線性化求解,具體過程如下[21]:基于虛擬材料法建立整環(huán)有限元接觸模型;通過靜力分析獲得接觸層在工作靜載荷下的剛度;如果整環(huán)工作時的載荷波動遠小于靜載荷,可認為整環(huán)在工作時接觸層的剛度保持不變,即用靜力分析獲得的接觸層剛度建立整環(huán)的線性模型,將非線性模型線性化;通過整環(huán)線性模型的模態(tài)分析獲得結(jié)合部的振型及固有頻率。

依據(jù)式(34)、(35)及LZT模型計算的法向、切向接觸剛度數(shù)據(jù)可得虛擬材料層的非線性接觸參數(shù),試樣和整環(huán)有限元模型的材料參數(shù)一致。采用APDL語言將虛擬接觸層的材料參數(shù)[22]、非線性接觸參數(shù)寫入整環(huán)有限元接觸模型,整環(huán)模型共使用80 798個SOLID185單元,其中結(jié)合部虛擬接觸層單元有1 858個,厚度為1 mm[23],如圖12所示。

圖12 整環(huán)有限元模型

對上述有限元接觸模型進行靜力分析可以求得虛擬接觸層在靜載荷下的接觸剛度,再利用APDL語言將虛擬接觸層剛度代入模態(tài)分析,可得各階固有頻率及振型[22-27]。本文計算了法向面壓為0.6 MPa時本文模型與LZT模型[16]的有限元模態(tài),圖13為本文的試驗與有限元前3階振型匯總,表3為3種模型不同面壓下的模態(tài)分析結(jié)果。

(a)第一階

由表3可知LZT模型的最大相對誤差絕對值為8.12%,而本文模型為4.32%,與LZT模型相比精度提高了46.80%。利用本文研究的接觸角分布函數(shù)及側(cè)接觸建模方法來計算有限元模型的固有頻率和振型,與試驗結(jié)果相比誤差最小,本文模型的精度得到了驗證。

6 結(jié) 論

(1)磨削表面微凸體大多為肩并肩的側(cè)接觸,本文對三點峰水平距離和峰的高度差進行定義,經(jīng)過分析發(fā)現(xiàn),二者近似滿足分形規(guī)律,首次建立了考慮微凸體水平距離分布的結(jié)合面接觸剛度分形模型。

(2)本文構(gòu)造了微凸體水平距離和峰的高度差的分布函數(shù),根據(jù)側(cè)接觸相對關(guān)系確定了接觸角θ的分布函數(shù),與Gorbatikh模型的角度分布函數(shù)趨勢相同,接觸角θ和加載次數(shù)是影響接觸角度分布規(guī)律的主要原因。

(3)本文建立的有限元模型與試驗結(jié)果最大相對誤差絕對值為4.32%,說明考慮側(cè)接觸水平距離及接觸角影響的接觸剛度分形模型更加貼近磨削結(jié)合面的真實接觸特性,可為研究側(cè)接觸問題提供新的思路。