王 民，孫鐵偉，董朝陽，孔德順，高相勝

(1.北京工業(yè)大學材料與制造學部先進制造技術重點實驗室，北京100124；2.電火花加工技術北京市重點實驗室，北京100191；3.中國鐵道科學研究院標準計量研究所，北京100081)

滾珠絲杠副具有傳動精度高、同步性能好、傳動效率高等優(yōu)良性能，已成為航空、航天、石油鉆井和科學測量儀器等行業(yè)所必需的部件[1]。在滾珠絲杠副的運行過程中，滾珠與滾道的摩擦和磨損，會降低其精度的穩(wěn)定性和使用壽命。

滾珠絲杠副中滾珠與螺母或絲杠滾道之間產生的接觸應力分析、表面變形以及微觀接觸特性是滾珠絲杠副摩擦學研究的基本問題。這些問題對滾珠絲杠副的磨損分析、疲勞壽命分析以及滾珠運動狀態(tài)分析具有重要意義。國內外許多學者對滾珠絲杠副的彈性變形和接觸應力進行了深入研究。劉暢等[2]對滾珠絲杠副的接觸變形進行了原理性闡述，建立了兩個粗糙的曲面之間的接觸多尺度力學模型。Liu等[3]建立了考慮彈性、彈塑性、塑性變形的表面微凸體接觸等接觸參數來計算總實際接觸面積和接觸載荷。Zhao等[4]在考慮滾珠幾何誤差的情況下，使用運動分析方法計算了載荷分布和軸向彈性變形。胡建忠等[5]采用赫茲接觸理論建立了滾珠絲杠副軸向接觸剛度計算模型。Wei[6]、Verl[7]等基于描述滾珠絲杠副接觸狀態(tài)的赫茲理論，對高速滾珠絲杠副進行磨損分析。Chen等[8]研究了影響滾珠絲杠副軸向接觸剛度的各個要素，分別對單螺母彈性變形和軸向剛度進行分析。

以上學者普遍采用赫茲接觸理論來對各種接觸情況進行近似分析。但是，在滾珠絲杠副中，接觸區(qū)域是由公稱半徑、滾道形狀、螺旋升角、滾珠半徑等共同決定的。由于絲杠滾道和螺母滾道具有螺旋升角，從而造成滾道截面形狀的非對稱性，且螺旋升角越大導致的非對稱性越嚴重，因此采用Hertz接觸理論會造成較大的誤差。

為了便于對一般3維型面兩接觸體之間進行精確的摩擦磨損分析和疲勞壽命分析，國內外許多專家根據最小余能原理和共軛梯度法（conjugate gradient method，CGM）[9-10]，采用半解析法對一般接觸問題進行了深入研究。Kubo[11]、Francis[12]等提出利用最小余能原理計算彈性接觸中的壓力分布和應力分布，這種方法的優(yōu)點在于可適用于不規(guī)則接觸表面的接觸力學分析。Wang等[13]基于最小余能方程求解干接觸問題，接觸壓力和切向應力使用共軛梯度法求解，通過解析法得到影響系數，由此計算表面變形或位移，并采用快速傅里葉變換（fast Fourier transform，F(xiàn)FT）加速表面變形的計算[14]。Ai等[15]采用半解析法進行粗糙表面的接觸問題分析。因此，作者基于最小余能原理的半解析法來求解高速滾珠絲杠副滾珠與滾道之間接觸應力分布，使計算結果更加精確。

1 接觸幾何描述

滾珠絲杠副的接觸區(qū)域的幾何形狀是由公稱半徑、滾道截面圓弧半徑、螺旋升角、接觸角以及滾珠直徑等參數共同決定的。計算中所用大導程滾珠絲杠副來源于NSK的產品手冊，其各項參數如表1所示。

表1 滾珠絲杠副結構參數Tab.1 Structural parameters of ball screw mechanism

為了對滾珠絲杠副進行運動學分析，首先需要建立滾珠絲杠副系統(tǒng)的坐標系[15]，如圖1所示。分別在絲杠滾珠接觸點A和螺母滾珠接觸點B建立坐標系AXAYAZA和BXBYBZB，原點A和B位于剛接觸上時（即沒有正壓力時）的接觸點， ZA軸和 ZB軸分別垂直于與接觸點處的共切平面，XA軸和XB軸平行于滾珠中心軌跡螺旋線的切向方向t，YA軸和YB軸由右手定則確定。

在圖1中，坐標系OWXWYWZW為慣性參考坐標系，使 ZW軸與絲杠的中心軸線重合，但不隨絲杠旋轉；OXY Z為與絲杠固聯(lián)并隨絲杠繞軸線旋轉的坐標系，且 Z軸與絲杠軸線重合，使X軸與滾珠球心軌跡形成的螺旋線相交；OHtnb為滾珠球心軌跡形成的螺旋線上任意一點的Frenet-Serret坐標系，各個坐標軸的方向隨著原點在曲線上的位置不同而發(fā)生變化。原點OH位于螺旋線上，且與螺旋曲線（即絲杠）的相對位置由底角θ確定；t為螺旋線的切線方向且沿著螺旋線的上升方向；n為螺旋線的主法線方向；b為螺旋線的副法線方向。

圖1 滾珠絲杠副系統(tǒng)的坐標系Fig.1 Coordinate system of ball screw mechanism

為了準確地描述滾道接觸點處的幾何形狀、接觸變形和接觸應力，需要建立滾道接觸點處的坐標系，如圖2所示。

圖2 滾珠絲杠副接觸坐標系Fig.2 Contact coordinate system of ball screw mechanism

由于滾珠的對稱性，在兩個接觸點坐標系下滾珠表面形貌的參數方程均相同，可表示為[16-17]：

式中，A和B分別為絲杠與滾珠、螺母與滾珠的接觸點。滾道表面幾何形貌在接觸點A和B處所建立的坐標系下的參數方程可通過坐標轉換的方法獲得。Frenet-Serret坐標系OHtnb中，接觸點A附近絲杠滾道截面圓弧的參數表達式為：

式中，rS為絲杠滾道截面圓弧半徑，rb為滾珠半徑，βA為滾珠與絲杠滾道之間的接觸角，βAS為所分析絲杠滾道區(qū)域截面圓弧段對應的圓心角，βA0用于定義絲杠滾道圓心角的范圍。

同理，可得接觸點B附近螺母滾道截面圓弧的參數表達式。

由滾道圓弧參數方程（3）和（4）可以得到滾道曲面形貌即絲杠滾道和螺母滾道接觸點附近區(qū)域在接觸點坐標系中的分布情況。但是，3維表面在接觸點坐標系的平面AXAYA和BXBYB上的投影并非規(guī)則的矩形區(qū)域，因此，需要通過2維插值算法，重新進行網格劃分，以便于后續(xù)的接觸力學分析。滾道曲面插值結果如圖3所示。

圖3 滾道曲面插值結果Fig.3 Interpolation resultsof raceway surface

2 接觸區(qū)域的法向應力分析

2.1 接觸區(qū)域的離散化

當一個球和一個彈性半無限體的接觸，在求解接觸應力分布的數值解時，首先要對接觸區(qū)域和控制方程進行離散。為計算方便，不考慮摩擦力對接觸應力計算的影響，給出的方程都是離散化形式，其中M和N分別為x與y方向的網格點數。兩表面的接觸間隙和接觸壓力滿足如下方程[14]：

式中，ci，j為接觸表面間隙矩陣c的元素，d0為剛體接近距離，cgi，j為接觸體幾何形狀矩陣cg的元素，cri，j為接觸體表面粗糙度矩陣cr的元素，Vi，j為表面彈性變形矩陣V的元素。

表面彈性變形的離散形式的計算方程為：

式中，pi，j為接觸壓力矩陣P的元素，ki-k，j-l為變形影響系數矩陣K的元素。

同時接觸間隙和接觸壓力要滿足的約束條件為：

在接觸區(qū)內，ci，j=0，pi，j≥0；在接觸區(qū)外，ci，j＞0，pi，j=0。顯然，ci，j和pi，j不能同時為0，即：

2.2 接觸變形的計算方法

對于兩彈性體接觸表面產生的彈性變形的計算，傳統(tǒng)的計算方法是在單元網格上用多項式近似壓力分布得到影響系數后采用直接矩陣相乘的方法來計算彈性變形。但為了達到較高的精度，需要對網格進行細分，這樣會顯著增加計算量，效率較低。Brandt等[18]提出了多重網格積分法，在求解帶有粗糙表面的接觸問題時可以得到較為精確的結果，且計算效率較高。

由于接觸表面變形的計算相當于求影響系數矩陣與法向應力之間的卷積，因此，在數值求解過程中，采用快速傅里葉變換（FFT）方法可以方便快速地計算彈性變形。通過該方法求解接觸變形量分布的計算步驟可參考文獻[14]。計算結果表明采用快速傅里葉變換計算彈性變形能夠達到所需精度，而且計算效率得到明顯提高。

2.3 共軛梯度法求解接觸壓力分布

接觸問題的核心就在于求解壓力分布和接觸體表面幾何形貌之間的關系，式（5）～（8）所描述的接觸問題實際上構成一個線性補余問題，可用矩陣形式表示[14]：

在約束條件ci，j≥0，pi，j≥0，ci，jpi，j=0下的極小值。

通過變分原理將接觸問題轉化為式（11）所描述的極值問題后，可以采用共扼梯度法進行循環(huán)迭代來實現(xiàn)快速收斂。作者所使用的共軛梯度法收斂性的詳細證明和迭代求解過程可參考文獻[14,19]。

3 計算結果驗證與分析

3.1 計算結果驗證

在研究滾珠絲杠副中滾珠與滾道之間的接觸狀態(tài)時，國內外許多學者均采用赫茲接觸理論來進行近似分析。但是由于滾道曲面形狀的復雜性和非對稱性，因此直接采用赫茲接觸模型計算是不精確的。為了驗證基于最小余能原理接觸應力分布的數值解即非赫茲解的正確性，將該原理與赫茲接觸理論分別用于求解光滑的球與平面之間的彈性接觸應力分布，并將兩種方法的計算結果進行對比分析[17]。

赫茲接觸理論[20]提出了針對兩個不同球體之間彈性接觸問題的經典解。根據該理論，使其中一個球體的半徑趨于無窮大，即可得到光滑圓球與光滑平面接觸時以半徑為a的圓形接觸區(qū)域內的赫茲接觸壓力分布的解析解即赫茲解。

式中，E1、E2和ν1、ν2分別為兩接觸體的彈性模量和泊松比，R1為球體的半徑，R2→+∞表示平面。

接觸區(qū)域半徑與法向趨近量之間的關系為：

根據式（12）～（15）即可求出光滑圓球與光滑平面接觸應力分布的情況，兩種方法的計算結果如圖4所示，其中實線表示接觸應力分布的赫茲解。從圖4中，可以看出非赫茲解與赫茲解的重合度很好，因此驗證了前述非赫茲解的正確性。

圖4 球與平面彈性接觸壓力分布的非赫茲解Fig.4 Comparison of the non-Hertzian and Hertzian solutionsfor the pressuredistribution of the elastic contact between a ball and a plane

3.2 接觸應力計算

為了更加準確、全面描述滾珠絲杠副的接觸應力分布情況，分別基于最小余能原理和采用赫茲接觸理論對滾珠絲杠副中滾珠與滾道的接觸問題進行詳細計算，并對比分析兩種方法的計算結果。

在求解過程中，法向正壓力取170 N并綜合考慮了滾道的實際3維曲面形貌，圖5和6分別給出了螺旋升角為32.48°時滾珠絲杠副接觸點A和接觸點B處接觸區(qū)域上應力分布的赫茲解、非赫茲解。由圖5和6可知，兩種求解方法所得接觸應力分布存在一定差異，這是由于其接觸區(qū)域的非對稱性變化造成的。同時，可以看到絲杠滾珠接觸點A處接觸應力比螺母滾珠接觸點B處接觸應力稍大一些，這是因為兩接觸點處的滾道3維表面形貌不同。

為了對比兩種方法計算結果的不同，圖7給出了螺旋升角為32.48°時滾珠絲杠副接觸點A和接觸點B處應力的非赫茲解與赫茲解之差的等高線圖。由圖7可知，兩種解法所得結果存在較大的不同，尤其是接觸區(qū)域的邊緣處應力差值最大，距離接觸中心越近差值越小。因此，采用赫茲接觸理論求解滾珠絲杠副的接觸應力將產生較大的誤差，而采用非赫茲解來分析滾道的接觸應力分布情況計算精度更高，且方法簡便。

圖5 螺旋升角為32.48°時接觸區(qū)域應力分布赫茲解Fig.5 Hertzian solution of stress distribution in contact area with helix rise angle of 32.48°

圖6 螺旋升角為32.48°時接觸區(qū)域應力分布非赫茲解Fig.6 Non-Hertzian solution of stress distribution in contact area with helix rise angle of 32.48°

3.3 接觸區(qū)域分析

圖7 螺旋升角為32.48°時接觸點處應力的非赫茲解與赫茲解之差的等高線圖Fig.7 Contour map of the difference between the non-Hertzian solution and the Hertzian solution of the stress at the contact point with helix rise angle of 32.48°

圖8和9為兩種方法計算得到滾珠與滾道接觸點處接觸應力分布的等高線圖。非赫茲解所得接觸橢圓相比赫茲解所得接觸橢圓存在一定的偏斜，而且通過計算可求得接觸區(qū)域的面積和長短軸，兩種解法存在一定差異。其中，接觸區(qū)域應力分布情況可以反映絲杠和螺母滾道因法向壓力引起的磨損量在滾道上形成的磨損帶的深度和寬度。此外，兩種計算方法所得接觸點B處接觸面積大于接觸點A處接觸面積，這與第3.2節(jié)中接觸點A、B處接觸應力存在不同相符合。

3.4 螺旋升角對接觸應力誤差的影響分析

為了進一步分析滾珠絲杠副螺旋升角的變化對赫茲接觸理論所得接觸應力計算誤差的影響，圖10分別給出了螺旋升角為4.55°、10.81°和17.66°時絲杠滾珠接觸點A處應力的非赫茲解與赫茲解之差的等高線圖，其中，Emax是滾珠絲杠副接觸點A處赫茲接觸應力相對于非赫茲接觸應力的最大計算誤差，由圖10可見該誤差隨著螺旋升角的增大而增大。

圖11給出了赫茲解相對于非赫茲解的接觸應力誤差隨螺旋升角的增大而逐漸增大的趨勢圖。由圖11可知，絲杠與滾珠接觸點A處的誤差始終大于螺母與滾珠接觸點B處的誤差。隨著機床的高速化、精密化的發(fā)展，螺旋升角逐漸變大，導致接觸區(qū)域的非對稱性更大。因此，當采用赫茲接觸理論計算滾珠絲杠副的接觸應力時，將會帶來更大的誤差。

圖8 赫茲解的接觸區(qū)域Fig.8 Contact area calculated by Hertzian method

圖9 非赫茲解的接觸區(qū)域Fig.9 Contact area calculated by non-Hertzian method

圖10 絲杠-滾珠接觸點A處應力的非赫茲解與赫茲解之差Fig.10 Difference between non-Hertzian solution and Hertzian solution of stress at contact point A

圖11 赫茲解相對于非赫茲解的接觸應力誤差Fig.11 Contact stress error of Hertzian solution relative to non-Hertzian solution

圖12中給出了滾珠絲杠副接觸點處非赫茲解和赫茲解的接觸應力峰值隨螺旋升角增大的變化情況。由圖12可知，滾珠絲杠副螺旋升角的變化對接觸點處赫茲解接觸應力峰值影響不大。但是，隨著螺旋升角的增加，絲杠滾道與滾珠接觸點A處非赫茲解接觸應力峰值逐漸變?。宦菽笣L道與滾珠接觸點B處非赫茲解接觸應力峰值逐漸增大，且接觸點A處應力峰值始終大于接觸點B處峰值。

圖12 接觸點處非赫茲解和赫茲解的接觸應力峰值Fig.12 Peak value of contact stressof non-Hertzian solution and Hertzian solution at the contact point

4 結 論

為了研究滾珠絲杠副滾道曲面的應力分布，并考慮接觸區(qū)域的非對稱性，克服采用赫茲接觸理論計算所帶來的誤差，基于最小余能原理對滾珠與滾道進行接觸應力分析，可以得到以下結論：

1）赫茲接觸解與最小余能非赫茲接觸精確解相比，當螺旋升角增加時，絲杠滾道和螺母滾道與滾珠接觸點處接觸應力誤差逐漸增大。而且，隨著滾珠絲杠副螺旋升角的增大，絲杠滾道與滾珠接觸點處的誤差始終大于螺母滾道與滾珠接觸點處的誤差。

2）赫茲解和非赫茲解兩種方法計算所得滾珠與滾道接觸區(qū)域的面積和長短軸有所不同。準確計算滾珠和滾道接觸區(qū)域應力分布可以更加全面和準確的計算滾道磨損帶的寬度和深度。

3）隨著數控機床高速高精化發(fā)展，具有大螺旋升角的高速滾珠絲杠副應用越來越廣泛。由于高速滾珠絲杠副接觸區(qū)域非對稱性大，有必要采用最小余能原理進行接觸應力分析，以便更準確的分析高速滾珠絲杠副的磨損機理和預測其精度保持性。