考慮彈塑性變形機制的結合面法向接觸剛度建模＊

張學良，陳永會，溫淑花，蘭國生，丁紅欽，王南山，張宗陽

（1.太原科技大學機械工程學院，山西 太原030024；2.東南大學機械工程學院，江蘇 南京211100；

3.西安交通大學機械工程學院，陜西 西安710049；4.徐州徐工基礎工程機械有限公司，江蘇 徐州221004）

引 言

機械結構都是組合結構，其中包含了各種各樣的結合面，結合面的存在使得機械結構不再具有連續性，導致了其結構靜動態特性的復雜性。結合面的接觸剛度常常是機械結構整體剛度的重要組成部分，如對機床結構而言，其結合面的接觸剛度約占機床總剛度的60%～80%［1］。關于結合面接觸剛度的研究［3～6，9～10，12～14］，往 往 僅 考 慮 了 結 合 面 微 凸 體彈性接觸變形和塑性接觸變形兩個階段，而忽略了微凸體的彈塑性接觸變形。文獻［2，18］分別首先提出了粗糙表面接觸的三維和二維接觸分形模型；文獻［3，4］最早建立了結合面法向接觸剛度分形模型，在此基礎上，文獻［4］研究建立了考慮微接觸域擴展因子影響的結合面法向接觸剛度分形模型，之后，文獻［7，9，10，12～14］再次對結合面法向接觸剛度模型進行了改進，但都沒有考慮彈塑性接觸變形，而僅僅考慮了彈性接觸變形和塑性接觸變形。文獻［15］提出了包含彈性變形、彈塑性變形和塑性變形情況下的接觸模型，但該模型對于彈塑性變形區間的預估不夠精確，與文獻［16］有限元分析結果存在較大誤差。文獻［7，8，11］基于文獻［16］關于彈塑性變形階段的分析結果，在考慮結合面微凸體彈性變形、彈塑性變形和塑性變形的基礎上，分別建立了結合面法向接觸剛度和靜摩擦因數的統計模型，并通過仿真分析揭示了結合面微觀形貌對結合面法向接觸剛度和靜摩擦因數的影響規律。但是，法向接觸剛度統計模型中的粗糙表面統計參數與取樣長度和測試儀器的分辨率有關，不具有客觀唯一確定性。鑒于此，本文基于接觸分形理論和MB模型及文獻［16］的彈塑性接觸變形機制，綜合考慮結合面微凸體的彈性接觸變形、塑性接觸變形和彈塑性接觸變形的兩個不同階段，建立了結合面法向接觸剛度的分形模型，并通過仿真揭示了彈塑性接觸變形機制及相關參數對結合面法向接觸剛度的影響，進而應用實驗實例驗證了該模型的正確可信性。該模型的特點是具有尺度獨立性，綜合考慮了結合面微凸體的彈性接觸變形、塑性接觸變形和彈塑性接觸變形的兩個不同階段，因此更加符合結合面實際接觸情況，具有較高的計算精度。

1 結合面微凸體接觸變形機制

結合面是由兩個粗糙表面接觸形成，通常將兩粗糙表面的接觸簡化為一粗糙表面和一剛性平面的接觸，對于粗糙表面的單個微凸體，將其等效為半球體，其等效曲率半徑為R，如圖1（a）所示。當在法向載荷p作用下，接觸區域半徑為r，法向接觸變形為δ，如圖1（b）所示。

圖1 等效微凸球體和剛性平面接觸前后狀態Fig.1 The two states of the equivalent microcontact sphere against a rigid flat plane before and after contact

1.1 微凸體彈性接觸變形階段

當微凸體的接觸變形在完全彈性變形范圍時，根據Hertz理論，球體接觸變形量為δ時，單個微凸體的圓形接觸面積ae和法向接觸載荷Fe分別為

式中E為復合彈性模量，，E2為兩接觸材料的彈性模量；ν1，ν2為兩接觸材料的泊松比。

1.2 微凸體彈塑性變形階段

文獻［16］利用有限元法對微凸體彈塑性變形階段進行了分析，當微凸體變形量大于彈性臨界變形量δc時，開始發生屈服現象，微凸體即進入彈塑性變形區域，彈性臨界變形量δc可表示為［17］

式中H為較軟材料的硬度，對于大部分延性材料，硬度和屈服強度Y之間關系為H＝2.8Y；根據文獻［17］，硬度系數K和較軟材料的泊松比ν之間的關系為K＝0.454＋0.41ν。

根據文獻［16］的研究結果，彈塑性變形區域的范圍介于δc≤δ≤110δc，并將彈塑性區域進一步劃分為兩個區域，變形量在δc≤δ≤6δc范圍內，屈服區域都發生在接觸面下方，而當變形量隨載荷逐漸增加時，接觸面下方的屈服區域逐漸擴大。當變形量為δ＝6δc時，屈服區域由接觸面下方擴大至接觸表面。當變形量為δ＝68δc時，屈服區域擴大至整個接觸表面。但當平均接觸壓力仍小于材料的硬度，當變形量增加至δ＝110δc時，平均接觸壓力才等于硬度。所以當變形量達到彈性臨界變形量δc的110倍，即δp＝110δc時，微凸體才進入完全塑性變形區域。有限元分析結果，得到微凸體彈塑性變形階段的接觸載荷、接觸面積和變形量之間的關系為［16］：

當δc≤δ≤6δc時，

當6δc≤δ≤110δc時，

1.3 微凸體塑性變形階段

當微凸體的變形量δ＞δp＝110δc時，進入完全塑性變形階段，微凸體的接觸面積和法向接觸載荷可以分別表示為：

2 結合面微凸體彈塑性接觸變形模型

2.1 結合面接觸分形基本理論

結合面微接觸點的面積分布密度函數為［18］

式中al為最大微接觸點的面積，D為粗糙表面的分形維數。

根據文獻［2，12］，結合面微凸體的變形量為

式中G為粗糙表面的分形粗糙度參數，α為大于1的常數，對于服從正態分布的隨機表面，α＝1.5較符合高頻密度和相位隨機的情況。

根據微凸體變形量與曲率半徑之關系，微凸體曲率半徑為［2，12］

當δ＝δc時，相應的接觸面積即為臨界接觸面積ac，且

2.2 結合面微凸體接觸載荷與接觸剛度模型

將式（11）和（12）代入式（2）中，得微凸體在彈性接觸變形階段的接觸載荷為

因此，單個微凸體在彈性接觸變形階段的法向接觸剛度為

將式（3），（11）和（12）代入式（5）中，得微凸體在彈塑性接觸變形第一階段的接觸載荷為

因此，單個微凸體在彈塑性接觸變形第一階段的法向接觸剛度為

將式（3），（11）和（12）代入式（7）中，得微凸體在彈塑性接觸變形第二階段的接觸載荷為

因此，單個微凸體在彈塑性接觸變形第二階段的法向接觸剛度為

3 考慮彈塑性接觸變形機制的結合面法向接觸剛度分形模型

結合面的真實接觸面積包括處于完全彈性接觸真實接觸面積、處于完全塑形接觸真實接觸面積和處于彈塑性接觸變形第一階段的真實接觸面積及處于彈塑性接觸變形第二階段的真實接觸面積，即

將式（10）代入（20），則

當結合面上的微凸體的法向接觸變形量δ≥110δc時，將處于完全塑性變形范圍，即當接觸變形處于彈性和彈塑性階段的微凸體才存在法向接觸剛度。

因此，結合面的法向接觸剛度為

將式（10），（15），（17）和（19）代入式（22）中，得

結合面法向總載荷與結合面真實接觸面積之間的關系為

當1＜D＜2，且D≠1.5時

當D＝1.5時

將式（23），（25）和（26）無量綱化，得

當D≠1.5時

當D＝1.5時

以上即為綜合考慮彈塑性接觸變形機制（即彈性接觸變形、彈塑性接觸變形第一階段和第二階段、塑性接觸變形）的結合面法向接觸載荷和法向接觸剛度分形模型，且法向接觸剛度與法向接觸載荷之間存在復雜的非線性關系。

需要指出的是，本文是將結合面（兩粗糙表面接觸）等效為一粗糙表面與一剛性平面的接觸，因此上述相關式中涉及到的參數D和G實質上是等效粗糙表面的參數，其確定方法參見文獻［5，13］。

4 結合面法向接觸剛度分形模型的仿真及結果分析

4.1 仿真及結果

給定一固定的無量綱真實接觸面積A＊r，便可以計算得到無量綱法向總載荷P＊和無量綱法向接觸剛度K＊n。一般工程表面的塑性指標范圍在0.7～2.5之間，所以給定φ分別為2.5，1.5和0.7，G＊分別為10－9，10－10和10－11，D分別取1.1～1.9，仿真計算結果如圖2～5所示，其中圖2中的“elastic and plastic”是指文獻［3］中僅考慮彈性和塑性變形機制的結合面法向接觸剛度分形模型的仿真計算結果，“elastoplastic”是指本文提出的考慮彈塑性接觸變形機制的結合面法向接觸剛度分形模型的仿真計算結果，其余圖中的仿真結果數據均由本文提出的模型仿真計算獲得。

圖2 K＊n隨P＊的變化規律（G＊＝1.0×10－10，φ＝1.5）Fig.2 Change curves of K＊nwith P＊（G＊＝1.0×10－10，φ＝1.5）

圖3 D對K＊n的影響曲線（G＊＝1.0×10－10，φ＝1.5）Fig.3 Influence curves of Don K＊n（G＊＝1.0×10－10，φ＝1.5）

圖4 G＊對Kn＊ 的影響曲線（φ＝1.5）Fig.4 Influence curves of G＊on Kn＊ （φ＝1.5）

圖5 φ對Kn＊ 的影響曲線（G＊＝1.0×10－10）Fig.5 Influence curves ofφon Kn＊（G＊ ＝1.0×10－10）

4.2 仿真結果分析

（1）從圖2（a）～（i）可以看出，綜合考慮彈塑性接觸變形機制的結合面無量綱法向接觸剛度K＊n大于僅考慮彈性和塑性變形機制的結合面無量綱法向接觸剛度值。可見，彈塑性接觸變形機制對結合面無量綱法向接觸剛度影響很明顯。當分形維數D＝1.1～1.4時，結合面無量綱法向接觸剛度隨著無量綱法向總載荷的增大而增大，且表現出較強的非線性關系；而當分形維數D＝1.4～1.9時，二者之間的非線性關系逐漸弱化并接近直線關系。

（2）從圖3（a），（b）及（b），（c）可以分別看出，當分形維數D＝1.1～1.6時，結合面無量綱法向接觸剛度隨著分形維數D的增大而增大；而當分形維數D＝1.6～1.9時，結合面無量綱法向接觸剛度隨著分形維數D的增大而減小。

（3）基于上述分析結果，這里分別選取了兩個代表性的分形維數值，即D＝1.3，1.7，從相應的仿真結果圖4（a），（b）和5（a），（b）可以分別看出，結合面無量綱法向接觸剛度隨著分形粗糙度G＊的增大而減小，卻隨著塑性指數φ的增大而增大。

5 結合面法向接觸剛度分形模型的實驗實例驗證

這里借助文獻［13］的結合面實驗裝置——啞鈴模型及其實驗數據，對本文所建立的結合面法向接觸剛度分形模型進行間接驗證。

啞鈴模型結合面的配對表面分別為HT250表面和45號鋼表面，其中 HT250表面參數為［13］：彈性模量E1＝116GPa，泊松比υ1＝0.27，密度ρ1＝7 340kg／m3，屈服強度σy1＝240MPa，硬度H1＝700MPa；45號鋼表面參數為［13］：彈性模量E2＝205GPa，泊松比υ2＝0.3，密度ρ2＝7 833kg／m3，屈服強度σy2＝353MPa，硬度H2＝500MPa。結合部的2個重要分形表征參數為［13］：D＝1.4241，G＝2.9820×10－10m。

每個螺栓的擰緊力矩分別控制為［13］：30，60，90N·m。

根據上述參數，應用本文提出的結合面法向接觸剛度分形模型，可以分別計算得到三種預緊面壓下的結合面法向接觸剛度，與之相應的結合面切向接觸剛度由文獻［5］中的模型3計算得到。

將計算得到的法向和切向接觸剛度以彈簧單元的等效模式嵌入到啞鈴模型整體的有限元模型中。啞鈴模型整體的有限元模型總共有5 100個單元、6 360個節點，其中彈簧單元108個，如圖6所示。

圖6 啞鈴模型有限元模型及結合面單元Fig.6 Finite element modeling and joint element of dumbbell

結合面的兩接觸面上各自對應均布36個節點，各對應節點之間嵌入3個COMBIN14單元（x，y，z三個方向，其剛度值即為上述所求得的相應的結合面接觸剛度的 1／36）。根據文 獻 ［13］，實驗裝置——啞鈴模型為自由狀態，因此將有限元模型的邊界條件設置為自由狀態。對有限元模型進行自由模態分析，得到前6階模態振型及相應的固有頻率，按照振型相似，模態固有頻率定量比較的原則，表1給出了計算固有頻率與文獻［13］實驗固有頻率的比較。文獻［13］和本文的計算結果均表明，啞鈴模型的前6階振型分別為繞z軸旋轉（1階）、繞y軸旋轉（2階）、繞x軸旋轉（3階）、沿x軸平動（4階）、沿y軸平動（5階）、沿z軸平動（6階）。顯然，繞y軸旋轉（2階）、繞x軸旋轉（3階）和沿z軸平動（6階）的三個模態與結合面法向接觸剛度關系密切，而本文關于這三個模態頻率的計算結果的相對誤差的絕對值除了第3階稍大于文獻［13］的結果外，其他兩階均明顯小于文獻［13］的結果，且本文相對誤差均為負值，而文獻［13］關于這三階計算結果的相對誤差有正有負，是不太合理的。其他三階模態與結合面切向接觸剛度關系密切，本文結果與文獻［13］的比較這里不再多述。

可見，基于本文所提出的模型計算結合面的法向接觸剛度進而用于進行啞鈴模型的模態計算具有較高的計算精度，說明該模型是正確可信的，同時也具有較高的精度。

表1 計算固有頻率與實驗固有頻率的比較Tab.1 Comparison of the computed and experimental natural frequencies

6 結 論

（1）基于接觸分形理論及MB模型，綜合考慮結合面微凸體接觸變形所處的彈性接觸變形、塑性接觸變形和彈塑性接觸變形的兩個不同階段，建立了結合面法向接觸剛度分形模型。

（2）彈塑性接觸變形機制對法向接觸剛度影響明顯，且考慮彈塑性接觸變形機制的結合面法向接觸剛度大于僅考慮彈性和塑性接觸變形機制所對應的法向接觸剛度。

（3）結合面無量綱法向接觸剛度隨著無量綱法向接觸載荷的增大而增大，且分形維數取值不同而呈較強的非線性關系（D＝1.1～1.4）或近似線性關系（D＝1.4～1.9）；結合面無量綱法向接觸剛度隨著分形維數增大而增大（分形維數D＝1.1～1.6）或者減小（分形維數D＝1.6～1.9）；結合面無量綱法向接觸剛度隨著無量綱分形粗糙度的增大而減小，而隨著塑性指數的增大而增大。

（4）實驗實例驗證說明了考慮彈塑性變形機制的結合面法向接觸剛度模型是正確可信的，具有較高的計算精度。

［1］ 張學良.機械結合面動態特性及應用［M］.北京：中國科技出版社，2002.Zhang Xueliang.Dynamic Characteristics of Machine Joint Surfaces and Its Applications［M］.Press of Science and Technology of China，Beijing，2002.

［2］ Yan W，Komvopoulos K.Contact analysis of elasticplastic fractal surfaces［J］.J.Appl.Phys.，1998，（84）：3 617—3 624.

［3］ 張學良，黃玉美，傅衛平，等.粗糙表面法向接觸剛度的分形模型［J］.應用力學學報，2000，17（2）：31—35.Zhang Xueliang，Huang Yumei，Fu Weiping，et al.Fractal model of normal contact stiffness between rough surfaces［J］.Chinese Journal of Applied Mechanics，2000，17（2）：31—35.

［4］ 張學良，黃玉美，韓穎.基于接觸分形理論的機械結合面法向接觸剛度模型［J］.中國機械工程，2000，11（7）：727—729.Zhang Xueliang，Huang Yumei，Han Ying.Fractal model of the normal con tact stiffness of machine joint surfaces based on the fractal contact theory［J］.Chinese Mechanical Engineering，2000，11（7）：727—729.

［5］ 溫淑花，張宗陽，張學良，等.固定結合面剛度分形模型［J］.農業機械學報，2013，44（2）：255—260.Wen Shuhua，Zhang Zongyang，Zhang Xueliang，et al.Stiffness fractal model for fixed joint interfaces［J］.Transactions of the Chinese Society for Agricultural Machinery.2013，44（2）：255—260.

［6］ 溫淑花，張學良，武美先，等.結合面法向接觸剛度分形模型建立與仿真［J］.農業機械學報，2009，40（11）：197—202.Wen Shuhua，Zhang Xueliang，Wu Meixian，et al.Fractal model and simulation of normal contact stiffness of joint interfaces and its simulation［J］.Transactions of the Chinese Society for Agricultural Machinery，2009，40（11）：197—202.

［7］ 尤晉閩，陳天寧.結合面靜態接觸參數的統計模型研究［J］.振動與沖擊，2010，29（11）：47—50.You Jinmin，Chen Tianning.Statistical model for static contact parameters of joint surfaces［J］.Journal of Vibration and Shock.2010，29（11）：47—50.

［8］ 尤晉閩，陳天寧.結合面靜摩擦系數的統計模型［J］.振動與沖擊，2010，29（12）：26—29.You Jinmin，Chen Tianning.Statistical model of static friction coefficient between joint surfaces［J］.Journal of Vibration and Shock，2010，29（12）：26—29.

［9］ 尤晉閩，陳天寧.基于分形接觸理論的結合面法向接觸參數預估［J］.上海交通大學學報，2011，45（9）：1 275—1 280.You Jinmin，Chen Tianning.Estimation for normal parameters of joint surfaces based on fractal theory［J］.Journal of Shanghai Jiao Tong University，2011，45（9）：1 275—1 280.

［10］尤晉閩，陳天寧.結合面法向動態參數的分形模型［J］.西安交通大學學報，2009，43（9）：91—94.You Jinmin，Chen Tianning.Fractal model for normal dynamic parameters of joint surfaces［J］.Journal of Shanghai Jiao Tong University，2009，43（9）：91—94.

［11］You J M，Chen T N.Statistical model for normal and tangential contact parameters of rough surface［J］.Proc.I Mech E，225（Part C）：171—185.

［12］Shuyun Jiang，Yunjian Zheng.A contact stiffness model of machined joint surfaces［J］.Journal of Tribology，Transactions of ASME，2010，（132）：1—7.

［13］田紅亮，趙春華，朱大林，等.金屬材料結合部法切向剛度修正與實驗驗證［J］.農業機械學報，2012，43（6）：207—214.Tian Hongliang，Zhao Chunhua，Zhu Dalin，et al.Modification of normal and tangential stiffness for joint interface with metallic material and experimental validation［J］.Transactions of the Chinese Society for Agricultural Machinery，2012，43（6）：207—214.

［14］Tian H L，Li B，Liu H Q，et al.A new method of virtual material hypothesis－based dynamic modeling on fixed joint interface in machine tools［J］.International Journal of Machine Tools ＆Manufacture，2011，51（3）：239—249.

［15］Zhao Y W，Maietta D M，Chang L.An asperity microcontact model incorporating the transition from elastic deformation to fully plastic flow［J］.ASME Journal of Tribology，2000，122（1）：86—93.

［16］Kogut L，Etsion I.Elastic－plastic contact analysis of a sphere and a rigid flat［J］.ASME Journal of Applied Mechanics，2002，69（5）：657—662.

［17］Chang W R，Etsion I，Bogy D B.An elastic－plastic model for the contact of rough surfaces［J］.ASME J.Tribol，1988，101（1）：15—20.

［18］Mjumdar A，Bhushan B.Fractal model of elastic－plastic contact between rough surfaces［J］.J.Tribol.，1991，（113）：1—11.