螺栓結合面微觀接觸特性的不確定性量化?

李 玲，呂國豪，林 紅，王晶晶，蔡安江

（西安建筑科技大學機電工程學院 西安，710055）

引言

機械設備中存在大量螺栓結合面，螺栓結合面接觸特性的研究受諸多不確定性因素影響，如接觸剛度和阻尼、表面粗糙度、裝配工藝和載荷等。這些因素使得機械系統的動態特性具有一定的不確定性，導致其動力學特性極為復雜，難以建立完整的理論模型［1-2］。這些不確定性因素的測量、建模和量化能夠促進和優化工程結構設計［3］。因此，研究結合面微觀接觸特性的不確定性，可以提高機械系統設計結果的可信度。

Greenwood 等［4］在Hertz 接觸理論的基礎上利用概率統計方法建立了粗糙表面之間的彈性接觸（Greenwood-Williamson，簡稱GW）模型。Chang等［5］利用微凸體體積守恒原理考慮了微凸體的塑性變形，并將GW 模型推廣到塑性（Chang-Etsion-Bogy，簡稱CEB）模型。Zhao 等［6］基于接觸力學理論，將粗糙表面接觸時微凸體發生的彈性變形和塑性變形聯系起來，提出了彈塑性微觀接觸（Zhao-Maietta-Chang，簡稱ZMC）模型。文獻［7-8］用有限元方法求解了半球形接觸問題，得到了彈性、彈塑性和全塑性狀態下的接觸本構規律，并將其推廣到整個結合面（Kougut-Etsion，簡稱KE）模型。Li 等［9］基于GW 模型，利用不動點迭代法獲得了基體變形影響下的微凸體變形，得到了更加精確的彈性接觸模型。上述概率統計模型受測量儀器和數學方法限制，使得模型有一定局限性［10］。Sayles等［11］通過對機械加工表面進行研究，發現機加工表面符合分形特征，并提出了一種基于分形理論模擬粗糙表面的方法。但是，無論是概率統計模型還是分形模型，都未能考慮表面形貌參數的不確定性對接觸特性的影響。

近年來，有學者開始考慮不確定性因素對螺栓連接接觸特性的影響。Silva 等［12］采用概率法，從物理影響和幾何參數兩方面分析了螺栓連接圓柱殼結構非線性振動和穩定性的影響。Hakula 等［13］在概率法的基礎上，求解出具有不確定性參數的復雜殼體的頻響函數，并采用蒙特卡洛法驗證了該方法的正確性。然而，在不確定性的分析中，采用概率法難以獲取不確定性參數的概率分布［14-15］，因此學者們提出了一種區間分析法，并將其應用于轉子［16-19］、齒輪［20-21］和螺栓［22-23］等的動力學不確定性分析中。Castelluccio 等［24］采用有限元模擬，研究了螺紋緊固件的模型輸入和形狀不確定性，并將緊固件建模中的不確定性和可變性的來源描述為幾何、材料、力學和方法論中的不確定性。

筆者考慮了螺栓結合面表面形貌參數的不確定性，通過分形理論表征了微凸體表面輪廓高度，采用矩譜法求解了螺栓結合面表面形貌參數區間，并根據中心極限定理將表面形貌參數區間轉變為符合微凸體高度分布的高斯分布函數。在此基礎上，通過蒙特卡洛法建立了GW 模型、改進的GW 模型和ZMC 模型對考慮結合面表面形貌參數不確定性的區間估計。研究結果可為螺栓結合面的動態性能預測與優化提供更加準確的理論依據和實驗參考。

1 粗糙表面形貌參數表征

1.1 表面輪廓高度表征

在室溫條件下，通過實驗測量了結合面的分形參數，試樣及其測試位置如圖1 所示。實驗儀器為放大倍數×1 000 的Talysurf 表面輪廓儀，表面測量儀的最小采樣間距為0.25 μm，采樣長度為3 mm；待測表面的尺寸為20 mm×20 mm，材料為20CrMo；測量選擇的樣品粗糙度Ra為6.529 6 μm。圖中紅色1，2，3 線段表示測量位置，由下向上進行測量。

圖1 試樣及其測試位置（單位：mm）Fig.1 Sample and its test location (unit:mm)

根據測得的實驗數據，結合結構函數法［24］得到結合面表面分形參數，Ra為6.529 6 μm 時分形維數D與分形粗糙度G的值如表1 所示。由表可知，同一被測表面多次測量的分形參數值是在一定范圍內波動的。D和G的區間范圍分別為［1.172 7～1.178 3］和［0.628 2～1.648 6］。

表1 Ra=6.529 6 μm 時分形維數D 與分形粗糙度G 的值Tab.1 Values of fractal dimension D and fractal roughness G when Ra=6.529 6 μm

為表征結合面表面輪廓高度，對分形參數進行模擬。根據分形理論［25］，二維輪廓高度曲線表達式為

其中：z（x）為微凸體高度；x為橫向距離；L為樣品長度；n為空間頻率指數；γ為標度參數，通常γ=1.5；φ1，n為隨機相位；nmax為最大頻率指數。

隨機選取了第3，5，6 組的數據，并代入式（1）中進行模擬。同一粗糙度下微凸體輪廓高度如圖2 所示，在粗糙度為6.529 6 μm 時微凸體輪廓高度峰值和谷值存在明顯差異，輪廓高度的數值隨著D和G的增大而增大。由D和G的值對微凸體輪廓高度分布的表征可以發現，即使分形維數與分形粗糙度系數的變化很小，對輪廓高度的影響依然顯著。因此，即使是同一被測表面，其表面形貌參數也具有不確定性。

圖2 同一粗糙度下微凸體輪廓高度Fig.2 The profile height of the micro convex under the same roughness

1.2 表面形貌參數的確定

為求解受分形參數影響下的表面形貌參數區間（粗糙密度η、粗糙半徑R和表面高度標準差σs），引入矩普法［25］。表面形貌參數計算式為

其中：m0，m2和m4分別為零階、2 階和4 階譜矩。

其中：Φ（δ）為功率譜函數；wh和wl分別為頻譜帶寬的上限和下限，wl=1/L，wh=1/T；L為樣本長度；T為測量儀器的最小采樣間距。

由式（2）～（7）可知，表面形貌參數的值與分形參數密切相關，并且表面形貌參數的變化直接影響著微凸體形態的變化。表面形貌參數區間如表2 所示，可通過實驗測得的分形參數區間來計算。

表2 表面形貌參數區間Tab.2 Surface topography parameter interval

2 結合面接觸參數的不確定性量化

2.1 形貌參數的不確定性表征

由于實驗次數及測量方法的限制，難以獲得完備數據。為此，采用蒙特卡洛法對表面形貌參數區間隨機抽樣，以彌補實驗次數及測量方法的不足。利用中心極限定理，將表面形貌參數區間轉變為與微凸體高度分布［4］一致的高斯分布函數，從而避免隨機抽樣誤差的累加造成的置信水平降低。

根據中心極限定理，設隨機變量x1，x2，…，xn相互獨立，服從同一分布且有有限的數學期望μ和方差σ2>0，對于任意x滿足

其中：Yn為隨機變量xi（i=1，2，…，N）之和；Fn（x）為Yn標準化的分布函數。

標準化后Var（xi）會存在微小的誤差，此時Yn是微小誤差的累加，即每一次抽取的隨機數都會產生誤差，經過至少1 000 次隨機模擬后產生的誤差將十分顯著。為減小抽取隨機數標準化后存在的誤差，假設隨機變量序列Xi相互獨立，且具有相同的期望和方差，即E（Xi）=μ，D（Xi）=σ2，令Yn=X1+X2+…+Xn，可得表面形貌參數服從高斯分布的表達式為

其中：Fn為表面形貌參數的分布律。

根據Fn→N（μ，σ2），能夠求解符合微凸體高斯分布的表面形貌參數區間，如表3 所示。

表3 高斯分布的表面形貌參數區間Tab.3 Gaussian distribution of surface topography parameter interval

2.2 確定性的結合面接觸參數

由圖（1）可知，結合面分形參數存在不確定性，其不確定性首先通過矩譜法傳遞至表面形貌參數。為分析表面形貌參數的不確定性對結合面接觸特性的影響程度，給出文獻［9］模型中接觸載荷、接觸面積和法向接觸剛度的表達式為

其中：下標※表示不確定性變量；f（zn－dn）為單個微凸體所受載荷；E為等效彈性模量；dn*為無量綱接觸間隙，dn*=sqrt（6）－2sqrt（3/N）；N=ηAn，An為名義接觸面積；R*為不確定性微凸體半徑；v為較軟材料的泊松比；λ為法向接觸剛度與切向接觸剛度的比值系數，一般取為2。

由于微凸體半徑R和粗糙密度η的變化，使得結合面實際接觸面積發生變化，從而導致結合面接觸載荷產生顯著變化。表面高度標準差σs直接決定了結合面微凸體分布特征，其很大程度影響著結合面接觸剛度。因此，結合面表面形貌參數的變化必將影響結合面接觸特性。

2.3 接觸參數不確定性量化

通過對實驗獲得的D和G的不確定性區間進行隨機抽樣，并根據中心極限定理得到服從高斯分布的表面形貌參數樣本Xi（i=1，2，…，nSE），nSE為樣本點的容量。表面形貌參數的不確定性按照前一次抽取的每個樣本點Xi進行傳遞。在給定樣本點Xi的情況下，一個樣本容量為nSA的隨機樣本aij（j=1，2，…，nSA）對應的結合面接觸特性x的概率分布函數［26］可以表示為

結合面接觸特性的最小值xm和最大值xM為

結合面接觸特性x為公式中的形式參數，實際參數分別為結合面接觸載荷F、接觸面積A、法向接觸剛度Kn與切向接觸剛度Kt，其概率密度函數為。

為量化結合面接觸特性的不確定性，采用圖3所示的雙循環的蒙特卡洛法抽樣技術。

圖3 雙循環的蒙特卡洛法抽樣技術Fig.3 Monte Carlo sampling technique with double loop

3 討論與分析

3.1 模型對比驗證

為驗證本研究考慮表面形貌參數不確定性的有效性，選擇GW 模型［4］、ZMC 模型［6］和本研究模型，將其法向和切向接觸剛度隨法向載荷的變化進行對比，并建立了3 種模型的不確定性帶寬，以驗證本研究模型的正確性與適用性。選用表3 和表4 的參數值作為仿真參數。

表4 仿真參數Tab.4 Simulation parameters

由于表面形貌參數存在不確定性，使得所建立的結合面接觸剛度存在不確定性帶寬，相較于求解結合面接觸特性的確定性模型，采用蒙特卡洛法能夠建立完整的結合面不確定性接觸模型。這3 種基于表面形貌參數不確定性模型的變化趨勢分別與確定性模型曲線一致，其中接觸剛度的不確定性帶寬隨著法向接觸載荷的增大而逐漸變寬。這是因為隨著結合面接觸更加充分，接觸的微凸體個數越來越多，從而導致表面形貌參數傳遞至結合面接觸載荷與接觸剛度的不確定性逐漸累加。

剛度模型對比驗證如圖4 所示。由圖4（a）可知，本研究模型的不確定性帶寬與GW 模型和ZMC模型的帶寬都有一定接觸，且位于二者的不確定性帶寬中間。同時可以發現，GW 模型和ZMC 模型的接觸剛度與法向接觸載荷的線性關系更為相似。這是由于本研究模型與GW 模型考慮了微凸體的彈性變形，而ZMC 模型考慮了微凸體的彈塑性變形，本研究模型在GW 模型對于微凸體彈性變形的基礎上，考慮了結合面間微凸體接觸對基體變形的影響，屬于對GW 模型的修正。因此，本研究模型的不確定性帶寬與GW 模型的不確定性帶寬更為接近。

圖4 剛度模型對比驗證Fig.4 Comparison and verification of stiffness models

圖4（b）為上述3 種模型切向接觸剛度與法向接觸載荷的關系曲線，可知切向接觸剛度與法向接觸剛度的變化規律基本相同。工程實際中，結合面的表面形貌是不確定的，相較將其等效為確定值，由蒙特卡洛法建立的不確定性帶寬更能反映細微的不確定性對接觸特性產生的影響，同時也能對結合面接觸特性進行合理估計。

3.2 結合面不確定性的接觸特性分析

接觸剛度與結合面無量綱接觸間隙的關系如圖5 所示，該圖顯示結合面在同一粗糙度下，本研究模型與GW 模型考慮不確定性因素時結合面接觸剛度K、法向接觸載荷F、實際接觸面積A與無量綱表面接觸間隙d*的關系。圖中的不確定性帶寬為隨著接觸間隙的改變，對結合面接觸特性的區間估計。結合面接觸間隙和接觸剛度的曲線，相比結合面接觸載荷與接觸剛度的曲線更能直觀反映機械系統結合面動態接觸過程的性能。

圖5 接觸剛度與結合面無量綱接觸間隙的關系Fig.5 The relationship between contact stiffness and the dimensionless contact gap of the joint surface

由圖5（a，b）可知，宏觀上，本研究模型與GW模型通過蒙特卡洛法建立的不確定性帶寬差異并不明顯，且接觸程度較為相似。由式（13）可知，切向接觸剛度和法向接觸剛度與材料的泊松比相關，則切向接觸剛度的不確定性帶寬一定與法向接觸剛度的不確定性帶寬緊密相連，都隨著表面接觸間隙的減小，結合面接觸剛度的不確定性帶寬逐漸增大，并且接觸也更加充分。微觀上，由其中的放大圖可知，不確定性帶寬的上、下限曲線并不像確定性模型的曲線一樣光滑，而是存在一定程度的波動，這是采用蒙特卡洛法對接觸剛度不確定性傳遞的結果，屬于對結合面微觀接觸特性的不確定性估計。

由圖5（c，d）可知，本研究模型與GW 模型的不確定性帶寬接觸程度發生了較大變化。這是因為當接觸間隙越來越小時，結合面間相接觸的微凸體個數也越來越多，當到達一定程度時基體發生變形，結合面間實際接觸面積增大的程度逐漸放緩，而接觸載荷受影響的程度卻相對較小，接觸間隙較大時接觸特性的不確定性帶寬較小，但未完全趨于零，這是因為結合面相互接觸時，粗糙表面上峰值較大的微凸體會首先發生接觸。

由圖5 可以發現，隨著接觸間隙逐漸減小，其帶寬也逐漸增大，并且接觸更加充分，即隨著結合面接觸間隙減小，表面形貌參數的細微變化對結合面接觸特性的影響越來越顯著。這是因為隨著接觸間隙逐漸減小，接觸的微凸體個數越來越多，從而使得表面形貌參數的不確定性傳遞至結合面接觸特性時，伴隨著不確定性的累加。因此，考慮表面形貌參數的不確定性對結合面接觸特性的影響不容忽視。

4 結論

1）隨著接觸間隙逐漸減小，結合面間接觸越來越充分。當到達一定程度時基體發生變形，結合面間實際接觸面積帶寬增大的程度逐漸放緩，而接觸載荷帶寬受影響的程度卻相對較小。

2）法向接觸載荷與接觸剛度的不確定性帶寬均隨著無量綱接觸間隙的減小而顯著增加，這說明隨著接觸的微凸體個數越來越多，結合面接觸參數的不確定性也在不斷累加。

3）基于多種結合面接觸模型，分別構建了其不確定性區間，為準確量化螺栓結合面的不確定性提供了理論依據與參考。因此，該建模方法具有成本低、效率高和適用性強的特點。