范立峰，趙 璐，聶 雯，劉小明

(1. 北京工業(yè)大學建筑工程學院，北京 100124；2. 河北工業(yè)大學土木與交通學院，天津 300401；3. 中國科學院力學研究所，北京 100084)

粗糙界面接觸剛度是影響機械力學性能的重要因素之一，復雜機械的界面接觸剛度可達整體靜態(tài)剛度的60%～80%[1]。界面接觸剛度受接觸荷載、表面粗糙度和基體力學性能的顯著影響[2?4]，研究界面接觸剛度與這些因素間的關系，對優(yōu)化機械界面設計、保障復雜機械工作性能具有重要的工程意義。

界面建模是研究界面接觸剛度的基礎，通常采用半球形[5 ? 9]、橢球形[10 ? 11]和正弦形[12 ? 20]對粗糙表面簡化建模。基于Hertz的球體彈性接觸解[21]和Hisakado[22]的表面形貌測量，Greenwood和Williamson[5]采用孤立彈性球形微凸體描述了粗糙表面的微觀形貌，首次建立了粗糙表面的接觸理論模型(GW模型)。結合微凸體塑性變形體積恒定假設，Chang等[6]推導出球形微凸體的彈-塑性接觸規(guī)律，將GW彈性接觸模型推廣至彈-塑性接觸統(tǒng)計模型(CEB模型)。Zhao等[7]采用樣板曲線描述了微凸體從彈性向完全塑性轉化的過程，進一步構建了粗糙表面的彈塑性接觸模型(ZMC模型)。Horng[10]將CEB球形彈-塑性接觸模型推廣至橢球形彈-塑性接觸模型。Jeng和Wang[11]將ZMC球形彈塑性接觸模型推廣至橢球形彈塑性接觸模型，研究了塑性指數對粗糙面平均間距和真實接觸面積的影響。此類理論模型被廣泛地應用于界面接觸性能預測。

另一方面，數值計算由于高效、便捷的特點被廣泛地應用于界面接觸剛度的研究。Kogut和Etsion[8]采用有限元方法模擬了球形微凸體的彈塑性變形演化過程，建立了微凸體的彈塑性變形模型(KE模型)。Zhao等[9]在彈塑性球形微凸體接觸有限元模型中考慮了冪律硬化和肩-肩接觸影響，建立了法向接觸剛度統(tǒng)計模型。Liu和Proudhon[12]將正弦微凸體變形演化劃分為彈性變形和塑性變形階段。隨后，Liu和Proudhon[13]又將正弦微凸體變形劃分為五個階段來進一步精確描述正弦微凸體變形演化過程。Song等[14]、Liu[15]和Yastrebov[16]在有限元模型中分別采用4.0倍、6.6倍和7.2倍的基體厚度-微凸體寬度比研究界面的接觸行為。為了考慮微凸體相互作用和基體力學性能的影響，姜英杰等[23]研究了高度和寬度隨機分布的粗糙界面接觸剛度，結果表明界面接觸特性受微凸體分布影響顯著。Gao等[17]研究表明：波狀表面基體厚度大于3倍波長寬度時，基體可近似為半空間體。因此，合理簡化粗糙表面模型、揭示力學機理的同時有效提高計算效率，對粗糙界面接觸行為的研究具有重要意義。

本文基于最小二乘法采用三角函數疊加來表征粗糙表面，建立具有不同特征尺度參數和分形維數的界面接觸有限元模型，確定基體最優(yōu)建模厚度。通過彈塑性有限元計算研究法向接觸壓力、法向接觸變形和法向接觸剛度的關系，揭示分形維數和特征尺度參數對法向接觸剛度的影響。

1 粗糙表面建模

1.1 三角函數疊加建模

粗糙表面測量表明正弦形接觸比球形接觸能更好地模擬粗糙表面微凸體接觸[18 ? 20]，三角函數疊加建模比球形建模更接近真實情況[24]。此外，由于正弦形微凸體底部與基體平滑連接，可以有效避免應力集中，因此三角函數疊加建模在微凸體塑性變形有限元模擬中具有更好的適應性[20]。Zhang等[20]以輪廓高度0位線為參考，將0位線以上的粗糙輪廓簡化為三角函數，如圖1所示。數學方程為：

式中：x1和xn為三角函數曲線與參考線交點；h為三角函數曲線高度；λ為三角函數曲線波長。

圖 1 粗糙表面和簡化的正弦表面Fig.1 Rough profile and simplified sinusoidal profile

采用最小二乘法基于最小均方差獲取最優(yōu)擬合函數，以此建立粗糙表面的最優(yōu)簡化模型[25]。三角函數曲線與粗糙表面間的最小均方差可表示為：

式(2)中，微凸體的高度h由?Er/?h=0求得：

通常采用分割三角函數的方法提高擬合精度。將三角函數谷值坐標zv與峰值坐標zp的比值定義為谷峰比，即zv/p=zv/zp。當谷峰比小于臨界值zc時，將原有三角函數劃分為兩個三角函數，以此提高擬合精度。現有研究結果顯示，當zc=0.30時，使用雙三角函數比單三角函數具有更高的擬合精度[26]。

1.2 粗糙表面簡化流程圖

采用谷峰比的臨界比值zc作為將單個微凸體原始粗糙表面簡化為兩個三角函數單元的判據。第一步，根據參考線與輪廓上離散點之間的關系，將分布在參考線同側的離散點作為一個微凸體，確定每個微凸體的起始點和結束點；第二步，基于第一步確定的微凸體單元，將滿足谷峰比zv/p≤0.30的微凸體劃分為兩個微凸體。此時，當離散點高度z(i)>z(i+1)且z(i)>z(i?1)時，則z(i)為峰值點；當離散點高度z(i)

圖 2 粗糙表面簡化流程圖Fig.2 Flow chart of rough profile simplified

1.3 粗糙表面的獲取與簡化

根據粗糙表面的分形特征，Weierstrass-Mandelbrot(WM)分形函數[27? 28]由于其處處連續(xù)但不可導和自放射性等特征被廣泛應用于擬合粗糙表面[27, 29 ? 31]。WM分形函數為：

式中：z(x)為粗糙表面的輪廓高度；x為表面輪廓的測量位置；G為粗糙表面的特征尺度參數，與輪廓曲線粗糙度的幅值有關；D為二維粗糙表面的分形維數，與輪廓曲線粗糙度的橫向密度有關；γ為控制表面起伏密度的參量，對于服從正態(tài)分布的表面，通常取γ=1.50(用來提供相位的隨機性和高頻譜密度的特性)；γn為隨機輪廓曲線的空間頻率，其大小對應于粗糙表面上單個微凸體波長的倒數γn=1/ln；nmin為與曲線結構的最低截止頻率指數；nmax為最高截止頻率指數。在實際應用中，表面輪廓的最高截止頻率也是有限的[32]。

根據現有研究，本文選取特征尺度參數G=10?12m～10?9m和頻率指數n=20～40來研究機械界面 接 觸 剛 度[28,32 ?33]。根 據 最 高 截 止 頻 率 指 數nmax和采樣間隔l的關系[34]，可以確定采樣間隔為0.05 μm。如圖2所示粗糙表面簡化流程圖以及式(4)，可以得到典型的粗糙表面三角函數擬合結果如圖3所示(分形維數D=1.40和特征尺度參數G=1.0×10?9m)。根據文獻[27, 29, 35]可知，采用三角函數疊加構建的粗糙表面輪廓具有顯著的分形特征，仍然為分形表面。本文對于構建得到的分形表面采用原始的分形維數D和特征尺度參數G表示。

圖 3 分形維數D=1.40和特征尺度參數G=1.0×10?9 m的粗糙表面和簡化的正弦曲線表面Fig.3 Rough profile with fractal dimension D=1.40 and characteristic length scale G=1.0×10?9 m and simplified sinusoidal profile

2 粗糙表面的有限元模擬方法

2.1 材料參數

金屬材料一般具有明顯的彈塑性力學行為，通常采用冪率硬化定律描述金屬的塑性行為[36]，其應力-應變本構關系如下[37]：

式中：E為材料彈性模量；σs為屈服強度；n為硬化指數。通常硬化指數0≤n≤1。當n=0時，表示材料為理想彈塑性；當n=1時，表示材料為彈性。本研究中，粗糙表面的材料為18CrMo4鋼。其中，彈性模量E=197.60 GPa；屈服強度σs=469 MPa；泊松比ν=0.29和硬化指數n=0.15[9]。

2.2 粗糙界面接觸的有限元模型

通常將三維表面進行降維處理來研究機械表面形貌[38?39]。此時，粗糙界面接觸可簡化為平面應變問題[40]。本文采用有限元軟件ANSYS建立粗糙表面有限元模型。如圖4所示為分形維數D=1.40和特征尺度參數G=1.0×10?9m的粗糙表面與剛性光滑平板接觸有限元模型。粗糙表面長度為L=3.0 mm，粗糙面以下基體厚度為Hs=5.0 mm。采用八節(jié)點PLANE183平面應變單元對粗糙表面和基體進行離散化，接觸單元(CONTA172)和剛性二節(jié)點目標單元(TARGE169)模擬接觸。

圖 4 分形維數D=1.40和特征尺度參數G=1.0×10?9 m的粗糙表面與剛性光滑平板接觸的有限元模型Fig.4 Finite element model of contact between a rigid smooth plate and a rough surface with fractal dimension D=1.40 and characteristic length scale G=1.0×10?9 m

本模型接觸算法采用增廣拉格朗日法計算無摩擦、無粘附接觸問題。粗糙界面接觸有限元模型邊界條件如下：1) 約束基體底部所有節(jié)點x方向位移和z方向位移，即Ux(x, 0)=Uz(x, 0)=0；2) 粗糙表面兩邊x方向位移和z方向位移自由，即Ux(0,z)和Ux(3,z)自由；3) 使用接觸向導，在剛性光滑平板最右端生成控制節(jié)點(Pilot節(jié)點)控制剛性光滑平板自由度；4) 在Pilot節(jié)點上施加將豎向位移δn。

為了驗證網格敏感性和邊界條件正確性，定義相對誤差Ep=(|PTH?PFE|/PTH)×100%，其中，PTH和PFE分別是理論計算荷載和有限元計算荷載，PTH=(δs/Hs)×(E/(1?ν2))。光滑表面上承受屈服荷載時，有限元模型計算的荷載和理論計算荷載的最大相對誤差為1.46%。結果表明：本研究邊界條件合理、網格尺寸滿足計算精度要求。

3 結果與討論

細觀微凸體廣泛存在于實際機械接觸表面，如圖5所示。機械界面接觸性能不僅與微凸體有關，也與基體的力學行為有關。因此，確定最優(yōu)的基體分析厚度，不僅能精確地反映接觸面力學行為，更能提高計算效率(在保證計算精度時，減少基體劃分時網格數量)，在界面接觸有限元數值模擬中尤為重要。本研究將參考線以下基體Hs分為基體最優(yōu)建模厚度Hs,1和基體厚度Hs,2，如圖5(b)。以18CrMo4鋼表面承受屈服荷載(即法向接觸壓力Pn=469 MPa)時，Von Mises應力分布規(guī)律確定基體最優(yōu)建模厚度Hs,1。

圖 5 機械界面Fig.5 Mechanical interface

3.1 特征尺度參數G對基體最優(yōu)建模厚度Hs,1的影響

圖6所示為不同特征尺度參數G的表面內Von Mises應力分布圖。從圖6(a)～圖6(d)可知，當分形維數D=1.40且表面法向接觸壓力Pn=469 MPa時，受微凸體隨機幾何分布特性影響，Von Mises應力在接觸表面附近區(qū)域整體復雜分布。對于單個微凸體，Von Mises應力從與剛性光滑平板接觸的峰值點向微凸體底部逐漸減小，直到與基體內Von Mises應力相等。圖6進一步顯示，未接觸區(qū)域基體內也會出現應力分布，這是由于相鄰微凸體相互擠壓，底部基體內應力相互影響形成的。

從圖6可知，隨著特征尺度參數G的增加，微凸體底部的應力從相鄰微凸體相互影響變?yōu)橄噜彾鄠€微凸體相互影響，如從圖6(b)～圖6(c)。隨著特征尺度參數G的增加，表面微凸體底部應力相互影響范圍增加，如圖6(c)和圖6(d)。在接觸區(qū)域附近應力分布不均勻，基體大部分區(qū)間內應力等于469 MPa。這說明特征尺度參數G對應力分布的均勻性有影響，表明表面下基體存在兩個明顯的區(qū)域，即應力均勻區(qū)域和不均勻區(qū)域。根據應力分布的均勻情況，確定應力不均勻區(qū)域為基體最優(yōu)建模厚度Hs,1。因此，選取粗糙表面的參考線與469 MPa的等值線間最小厚度作為基體最優(yōu)建模厚度Hs,1。從圖6可知，當分形維數D一定時，隨著特征尺度參數G的增加，表面下應力不均勻分布的深度增加，基體最優(yōu)建模厚度Hs,1增加。在本研究中，對于特征尺度參數分別為G=1.0×10?12m、1.0×10?11m、1.0×10?10m和1.0×10?9m的粗糙表面，基體最優(yōu)建模厚度Hs,1分別為0.14 mm、0.18 mm、0.45 mm和0.74 mm。

圖 6 在表面法向接觸壓力Pn=469 MPa時，具有相同分形維數D=1.40和不同特征尺度參數G的表面下Von Mises應力等值線云圖Fig.6 The Von Mises stress contour map of surface with the same fractal dimension D=1.40 and different characteristic length scale G under normal contact pressure Pn=469 MPa

3.2 分形維數D對基體最優(yōu)建模厚度Hs,1的影響

圖7所示為不同分形維數D的表面內Von Mises應力分布圖。從圖7可知，對于特征尺度參數G=1.0×10?9m和不同分形維數D的表面，當表面法向接觸壓力Pn=469 MPa時，隨著分形維數D的增加，微凸體底部的應力從相鄰多個微凸體相互影響變?yōu)橄噜徫⑼贵w相互影響，如從圖7(b)～圖7(c)。隨著分形維數D的增加，表面微凸體底部應力相互影響范圍減小，并且應力不向表面中心下方聚集，如圖7(c)和圖7(d)。在接觸區(qū)域附近應力分布不均勻，基體大部分區(qū)間內應力等于469 MPa。這說明分形維數D對應力分布的均勻性有影響，表明表面下基體存在兩個明顯的區(qū)域，即應力均勻區(qū)域和不均勻區(qū)域。從圖7可知，當特征尺度參數G一定時，隨著分形維數D的增加，表面下應力不均勻分布的深度減小，基體最優(yōu)建模厚度Hs,1減小。采用類似的定義方法可知，對于具有相同特征尺度參數G=1.0×10?9m和分形維數分別為D=1.40、1.50、1.60和1.70的粗糙表面，基體最優(yōu)建模厚度Hs,1分別為0.74 mm、0.37 mm、0.15 mm和0.09 mm。

圖 7 在法向接觸壓力Pn=469 MPa時，具有相同特征尺度參數G=1×10?9 m和不同分形維數D的表面下Von Mises應力等值線云圖Fig.7 The Von Mises stress contour map of surface with the same characteristic length scale G=1×10?9 m and different fractal dimension D under normal contact pressure Pn=469 MPa

3.3 特征尺度參數G對界面接觸行為的影響

3.1節(jié)討論了特征尺度參數G對最優(yōu)建模厚度Hs,1的影響，本部分在表面接觸行為中考慮基體變形的影響。具有相同分形維數D=1.40和不同特征尺度參數G的表面接觸行為如圖8所示。從圖中可知，法向接觸壓力Pn隨著法向接觸變形δn增加而非線性增加。在相同的法向接觸變形δn下，法向接觸壓力Pn隨著特征尺度參數G而減小。圖8(b)所示為法向接觸剛度Kn和法向接觸壓力Pn的關系。從圖中可知，法向接觸剛度Kn隨著法向接觸壓力Pn增加而非線性增加。在相同法向接觸壓力Pn下，法向接觸剛度Kn隨著特征尺度參數G的增加而減小，這一規(guī)律與前人研究[33, 41 ? 42]一致。這是因為特征尺度參數G的增加意味著表面粗糙度Ra增加[35,43]，Ra是輪廓曲線的算術平均值以上微凸體高度的平均值[28]。對于相同分形維數D和不同特征尺度參數G的表面，隨著特征尺度參數G增加，粗糙輪廓形狀不變，但整個粗糙輪廓的粗糙幅值增加[28 ? 29]。

圖8(c)所示為在不同法向接觸壓力Pn和相同分形維數D下，特征尺度參數G對表面法向接觸剛度Kn的影響。法向接觸壓力Pn分別取0.2σy、0.4σy、0.6σy、0.8σy和1.0σy，σy為材料的屈服強度。從圖中可知，法向接觸剛度Kn隨著特征尺度參數G的增加而非線性減小。在相同特征尺度參數G下，法向接觸剛度Kn隨著法向接觸壓力Pn的增加而增加。

3.4 分形維數D對界面接觸行為的影響

3.2節(jié)討論了分形維數D對最優(yōu)建模厚度Hs,1的影響，本部分在表面接觸行為中考慮基體變形的影響。具有相同特征尺度參數G=1×10?9m和不同分形維數D的粗糙界面接觸行為如圖9所示。圖9(a)所示為法向接觸壓力Pn和法向接觸變形δn的關系。從圖中可知，法向接觸壓力Pn隨著法向接觸變形δn的增加而非線性增加。在相同法向接觸變形δn下，法向接觸壓力Pn隨著分形維數D的增加而增加。圖9(b)所示為法向接觸剛度Kn和法向接觸壓力Pn的關系。從圖中可知，法向接觸剛度Kn隨著法向接觸壓力Pn的增加而非線性增加。在相同法向接觸壓力Pn下，法向接觸剛度Kn隨著分形維數D的增加而增加，這一規(guī)律與前人研究[33, 41, 44]一致。這是因為分形維數D的增加意味著表面粗糙度Ra減小[35]。對于相同特征尺度參數G和不同分形維數D的表面，隨著分形維數D的增加，粗糙表面輪廓沿橫向擠壓，輪廓變得越來越密集和復雜，但整個表面的輪廓粗糙幅值減小，表面變得越來越光滑[28 ? 29, 45]。

圖 8 相同分形維數D=1.40下，特征尺度參數G對表面接觸行為的影響Fig.8 Influence of characteristic length scale G on contact behavior of surface with the same fractal dimension D=1.40

圖 9 相同特征尺度參數G=1×10?9 m，分形維數D對表面接觸行為的影響Fig.9 Influence of fractal dimension D on contact behavior of surface with the same characteristic length scale G=1×10?9 m

圖9(c)所示為不同法向接觸壓力Pn和相同特征尺度參數G下，分形維數D對表面法向接觸剛度Kn的影響。法向接觸壓力Pn分別取0.2σy、0.4σy、0.6σy、0.8σy和1.0σy，σy為材料的屈服強度。從圖中可知，法向接觸剛度Kn隨著分形維數D的增加而非線性增加。在相同的分形維數D下，法向接觸剛度Kn隨著法向接觸壓力Pn的增加而增加。

3.5 有限元結果和理論模型對比

將本文有限元結果與文獻[44]中建立的不考慮微凸體相互作用和基體變形的彈塑性法向接觸剛度分形理論模型(簡稱分形理論模型)進行了對比，如圖10所示。從圖10可知，對于具有相同分形維數D和不同特征尺度參數G的表面，在相同的法向接觸壓力Pn下，分形理論模型計算的法向接觸剛度Kn比有限元計算的法向接觸剛度Kn大。造成此現象的原因是分形理論模型沒有考慮微凸體相互作用和基體變形的影響[33,46]。同時，特征尺度參數G=1×10?11m的表面其有限元結果與分形理論模型結果比特征尺度參數G=1×10?9m的表面更接近。原因是特征尺度參數G=1×10?9m的表面受到微凸體相互作用和基體變形影響比特征尺度參數G=1×10?11m的表面更大，而分形理論模型沒有考慮微凸體相互作用和基體變形。另外，Wang等[33]的研究也表明不考慮微凸體相互作用和基體變形的接觸剛度比試驗結果大，考慮微凸體相互作用的接觸剛度模型與試驗吻合性更好，說明本文考慮的基體變形和微凸體相互作用是不可忽略的。

圖 10 具有分形維數D=1.40，特征尺度參數分別為G=1×10?11 m和G=1×10?9 m表面，有限元模型結果和分形理論模型結果對比Fig.10 Comparison between finite element model and fractal theoretical model for surface with fractal dimension D=1.40,characteristic scale parameter G=1×10?11 m and G=1×10?9 m respectively

4 結論

本文建立了粗糙表面與剛性光滑平板接觸有限元模型，基于接觸面下Von Mises應力分布特征，提出了基體最優(yōu)建模厚度Hs,1?；诨w最優(yōu)建模厚度Hs,1研究了不同分形維數和特征尺度參數條件下粗糙表面的接觸行為。進一步討論了分形維數和特征尺度參數對法向接觸剛度的影響。主要得到以下幾個結論：

(1)對于具有相同分形維數D和不同特征尺度參數G的粗糙表面，或者相同特征尺度參數G和不同分形維數D的粗糙表面，均存在基體最優(yōu)建模厚度Hs,1?；w最優(yōu)建模厚度可以在相同計算精度下有效提高計算效率。

(2)在相同法向接觸壓力Pn下，表面的法向接觸剛度Kn隨著特征尺度參數G的增加而非線性減少。

(3)在相同法向接觸壓力Pn下，表面的法向接觸剛度Kn隨著分形維數D的增加而非線性增加。