王 豪,黃為彬,方 兵*
(1.福建農林大學 機電工程學院,福建 福州 350002;2.福建省農業信息感知技術重點實驗室,福建 福州 350002)
目前,基于計算機輔助設計的理念已經廣泛地應用于精密數控機床產品的開發中。設計者期望在設計階段準確掌握產品的剛度、熱特性、動態響應以及加工精度等性能。但由于機床零部件數量眾多、結構復雜,存在多種類型的結合面;結合面不連續的特點對機床剛度、熱傳導特性的影響顯著,且結合面自身的影響因素較多、規律不明顯。因此,結合面的建模準確與否在較大程度上影響整體分析的結果[1-3]。
在處理結合面問題時,“彈簧-阻尼”法是常用的方法之一。該方法用若干個“彈簧-阻尼”單元來模擬結合面的靜、動態特性。如朱堅民等人[4]使用了“彈簧-阻尼”單元模擬主軸刀柄結合面,并進行了動力學參數識別,結果表明實驗值與計算值相對誤差在4%以內,具有較高精度。李院生等人[5]采用了4節點的“彈簧-阻尼”單元模型來等效結合面,解決了螺栓固定結合面參數識別問題。蘭國生等人[6]采用了“彈簧-阻尼”法,建立了固定結合面法向接觸阻尼模型及結合面間阻尼損耗因子模型。
然而“彈簧-阻尼”單元法存在一些缺點。該方法是一種簡化自由度的方法,需要根據具體的問題,選擇單元的數量,確定單元的位置,并且忽略了各單元之間的相互影響。
為了解決上述問題,近年來研究人員已經開始采用虛擬材料法建立固定結合面模型。
虛擬材料法是一種根據結合面接觸部位微觀結構,采用具備彈性模量、泊松比以及密度的材料來模擬結合面靜力學和動力學特性的方法。比如,張學良等人[7]提出了一種利用宏觀各向同性虛擬材料等效的固定結合部動力學參數化建模方法,即采用虛擬材料的彈性模量、泊松比等來模擬結合面靜力學和動力學特性,結果顯示,計算值與實驗值的絕對誤差在10%以內。孫偉等人[8-9]采用了一層虛擬材料來對含螺栓復合材料結構進行動力學分析,均達到了較高的精度。
然而,上述研究并未考慮彈塑性變形對虛擬材料屬性的影響。
因此,筆者基于Greenwood-Williamson(G-W)接觸分析模型,考慮彈塑性變形的影響,依據粗糙表面形貌統計參數,提出一種新的考慮摩擦力影響的剛度模型;根據結合面的剛度特性,將結合面等效成一層虛擬材料,建立固定結合面的靜力學分析模型,將其應用于整機的有限元分析,并且開展機床靜剛度測試實驗,以驗證上述分析結果的正確性。
結合面實質上是由兩個相互接觸的粗糙表面組成,可以將其簡化為一個粗糙表面與一個剛性平面接觸。
GREENWOOD J A等人[10]早在20世紀60年代就提出了基于統計理論的G-W接觸分析模型。該模型假定凸峰的高度服從某種特定的分布(如高斯分布),凸峰具有相同的曲率半徑。當凸峰的變形量ω大于彈性極限值ωe時,凸峰開始呈現塑性變形的特點。然而,KOGUT通過有限元分析的方法發現,凸峰由彈性變形階段開始進入塑性變形階段時,先要經歷一個彈塑性混合變形階段;然后,直到其變形量大于等于110×ωe時,凸峰才呈現完全塑性變形的特點[11]。因此,單個凸峰的接觸過程可由彈性、彈塑性和塑性三個階段來描述。
根據Hertz接觸理論,可建立彈性接觸面積ae、彈性接觸載荷pe,以及彈性變形量ω間的關系如下:
ae=πRω
(1)
(2)

當結合面間存在摩擦力時[12],由凸峰發生屈服的臨界平均接觸壓力pave=kυkμH可得到彈性變形極限,即:
(3)
式中:kυ為平均接觸壓力系數,由材料泊松比υ確定[13],kυ=0.464 5+0.314 1υ+0.194 3υ2(若取υ=0.3,則kυ≈0.6);kμ為摩擦力影響系數,其值與摩擦系數μ有關[9];H為材料硬度。
kμ的計算公式如下:
(4)
式中:μ為摩擦系數。
(5)
根據文獻[14]可知,單個凸峰的切向剛度為:
(6)

Qe和Pe分別表示該階段的總體切向載荷和法向載荷,其中:
(7)
式中:η為凸峰密度;A0為名義接觸面積;z為凸峰高度;d為剛性平面與粗糙表面中線的距離,且滿足關系式:ω=z-d;ωe為彈性極限值;φ(z)為高度分布的概率密度函數。
文獻[15]給出了單個凸峰所承受的切向載荷為:
(8)
式中:σy為與摩擦力方向相同的應力;a為接觸點面積;p為接觸點上的載荷。
根據接觸面積與變形的關系,可得該階段總切向載荷為:
(9)
當變形量在ωe和完全塑性臨界值ωp之間時,微凸峰處在彈塑性混合變形階段。
李玲等人[16]在接觸面積、平均接觸載荷滿足連續變化的基礎上,構造了以變形量ω為自變量的多項式表達式,其接觸面積可表示為:
aep=πRκ1(ω)
(10)
其中:
接觸載荷可表示為:
(11)
類似的可定義剛度為接觸載荷對變形量的導數,于是法向接觸剛度可以表示為:
(12)
其中:
κ3(ω)=-2g3(ω)+6(2-ω)g2(ω)+2(3ω-4)g(ω)+2。
則,切向接觸剛度可以表示為:
(13)
式中:Pep,Qep為該階段的總體切向載荷和法向載荷。
Pep和Qep表達式如下:
(14)

(15)
KOGUT L等人[11]通過有限元計算發現,當變形量ω大于110×ωe時,凸峰呈現完全塑性變形的特點。此時,接觸面積、載荷與變形量的關系可以表示為:
ap=2πRω
(16)
pp=Hap
(17)
類似的,可以求得該階段的法向接觸剛度和切向接觸剛度分別為:
(18)
(19)
Pp和Qp表示該階段的總體切向載荷和法向載荷,分別為:
(20)

(21)
綜合式(5)、式(12)和式(18),可得到總的法向接觸剛度為:

(22)
綜合式(6)、式(13)和式(19),可得到總的切向接觸剛度為:
(23)
采用虛擬材料模擬結合面的靜力學特性,需要知道虛擬材料的密度ρ(取兩接觸材料的平均密度)、彈性模量Ev、泊松比μv、虛擬材料的面積A和厚度t等參數。
根據機械加工表面形貌特征以及有限元分析的需求,參考文獻[17],可選取虛擬材料的厚度t為1 mm。
等效虛擬材料模型如圖1所示。

圖1 等效虛擬材料模型
若單位面積上的法向剛度和切向剛度分別為(量綱為[N·m-3]):
(24)
(25)
取一面積為ds的微元,在法向載荷Fn的作用下,其法向應力和應變分別為σ和ε,則有:
(26)
式中:δn為法向變形量。
由彈性模量定義及化簡式(14),得到:
(27)
式中:Ev為虛擬材料的等效彈性模量。
同理,可得虛擬材料的等效切變模量Gv為:
(28)
式中:τ為切變應力;γ為切變應變。
泊松比μv為:
(29)
某立式加工中心的床身與立柱材料均為HT200,接觸表面采用銑削方式加工。在機床X、Y、Z三個方向各布置一根絲杠分別用于驅動主軸箱、工作臺和滑板。主軸前端安裝錐度為7:24的BT50刀柄,軸承采用4+2的布置方式,即前端采用4列并列,后端采用2列并列的方式配置,型號分別為7020C和7018A。
筆者對模型進行必要的簡化,并在CATIA V5r16中建立床身、立柱、滑板、工作臺、主軸箱以及導軌等主要零部件模型,按設計尺寸裝配成整機三維模型,再將CATIA V5r16中的整機三維模型導入ANSYS Workbench中,進行網格劃分、參數設定以及邊界條件設置等,以建立有限元分析模型[18]。
得到結合面的虛擬材料屬性是建立整機靜剛度分析模型的關鍵。此處以床身立柱結合面為列,該結合面由16只螺栓固定,正壓力約為56 000 N,接觸面積為0.238 m2,平均接觸壓力為2.35×105Pa。采用輪廓儀(泰勒霍普森Form Talysurf i120)測量所研究表面的形貌,并根據文獻[19]所提供的方法,經計算得到服從高斯分布的表面凸峰的高度分布函數φ(z),平均半徑β為42.68 μm以及單位面積凸峰個數N為8.45×103個/mm2。
平均接觸壓力Pav與彈性模量Ev的關系如圖2所示。

圖2 平均面壓與彈性模量間的關系
由式(7)、式(14)、式(20)和高度分布函數φ(z)可得到總載荷P與偏差d間的關系,由式(24)和式(25)可分別得到彈性模量Ev、μv和偏差d間的關系。
綜合上述兩個關系,便可得到平均接觸壓力Pav與彈性模量Ev。隨著平均接觸壓力Pav的增加,虛擬材料的彈性模量也在增加,但趨勢趨于緩慢。
所建立的有限元分析模型主要考慮了床身立柱結合面、軸承內圈與主軸的結合面、軸承外圈與箱體的結合面,以及主軸前后兩個端面;軸承、絲杠以及導軌則采用彈簧單元法進行等效。
根據實際工況,可得到主要結合面的虛擬材料參數如表1所示。

表1 主要結合面的虛擬材料屬性
仿真計算時,在主軸刀柄下端面沿著X、Y、Z三個方向施加大約1 960 N的載荷,反向載荷施加在工作臺面相應位置上,并將床身底座與地面的接觸面固定約束。
仿真結果如圖3所示。

圖3 施加1 960 N載荷的整機變形
可見,主軸在承受載荷時,變形主要體現在主軸及其周邊的零部件,其中主軸前端影響最為明顯。在相同的載荷作用下,X、Y和Z方向的位移分別為68.6 μm、44.5 μm和47.3 μm,三個方向對應的剛度分別為28.57 N/μm、44.04 N/μm和41.43 N/μm??傮w來說,Y方向的剛度最大,Z方向上的剛度次之,X向剛度最小。
為驗證整機靜剛度有限元分析結果的正確性,筆者開展了機床靜剛度測試實驗。
實驗測試原理如圖4所示。

圖4 機床靜剛度測量
筆者將機床各零部件放置在與仿真分析時對應的位置上,以98 N為間隔逐步增加載荷至1 960 N,采用精度為0.98 N的BK-2型壓力傳感器測量載荷,由二次儀表顯示測量結果,采用分辨率為2 μm的千分表測量相對變形。
實驗結果如圖5所示。

圖5 載荷-變形曲線
隨著外加載荷的增加,主軸末端的相對變形也在逐漸增加,測得的“載荷-變形”數據基本呈線性規律??蓴M合得到X、Y、Z三個方向的剛度,分別為26.60 N/μm、41.87 N/μm和40.17 N/μm。
為研究結合面對整體剛度特性的影響,筆者對未考慮結合面特性的模型進行分析計算。所謂不考慮結合面是指忽略由連接部位(也就是結合面)引起的剛度的變化。
在有限元模型中,筆者直接將連接部位黏接成一體,使得結合部與被連接部件的材料屬性一致。
實驗結果與仿真結果如表2所示。

表2 整機靜剛度仿真與實驗對比
其中,仿真的結果分為兩個模型,一個是考慮了結合面影響的模型,另一個是沒有考慮結合面影響的模型。
對比結果表明:不考慮結合面的影響時,機床工藝系統的靜剛度與實驗值的相對誤差達到17%左右;考慮結合面影響之后,其相對誤差基本在5%以內。以上結果驗證了所提出模型的準確性。
針對機床數字化設計中典型結合面建模精度誤差較大的問題,筆者提出了一種采用虛擬材料等效結合面靜力學特性的方法,建立了考慮結合面特性影響的機床整機有限元靜力學分析模型,并完成了相應的靜剛度測試實驗。
筆者基于粗糙表面接觸理論,推導了考慮摩擦力影響的固定結合面的法向剛度和切向剛度解析表達式,并確定了用于模擬結合面靜力學特性的虛擬材料的等效彈性模量和泊松比。
研究結果表明:
1)機床典型結合面等效后的虛擬材料彈性模量較基體材料小1到2個數量級。因此,對整機分析時,必須考慮結合面的影響;
2)與實驗測試結果對比,不考慮結合面影響的機床整機有限元靜力學分析模型的誤差在13%~17%之間,與實際偏差較大;
3)考慮了結合面特性影響后,得到機床整機有限元靜力學分析模型的誤差基本在5%以內,精度滿足工程設計的要求。
目前,筆者主要采用虛擬材料法研究了結合面的靜力學特性。后續工作中,筆者將在本研究的基礎上,將該模型應用于結合面的動態特性以及熱-結構耦合分析等方面。