接觸問題的三角形載荷離散FFT加速算法

陳楠 朱凱 蔣志楨 龔詩雨 李璞 金曉清

摘要：接觸問題控制方程的有效求解，往往涉及到復雜的數學理論知識，而在實際工程應用中，接觸應力分布又具有高度隨機性。為高效快速求解任意載荷分布下固體的接觸響應，基于三角形載荷離散單元，嵌入離散卷積快速傅里葉變換（DC-FFT）算法，提供了一種高精度、高可靠度的計算方法。相比于通常采用的分段均布載荷離散方法，三角形單元的解析求解略顯復雜，但能更好地模擬接觸載荷任意分布的特性，對于接觸邊緣處載荷由零遞增或遞減為零的情況，也可以予以充分考慮。為優化三角形載荷離散單元的求解方法，根據接觸影響系數矩陣的“激勵-響應”特性，推導了三角形載荷單元和均布載荷單元作用下的應力分量解析解。通過構造包含影響系數矩陣的離散卷積形式應力解，將某一目標節點在所有載荷單元作用下，重復度極高的矩陣運算疊加過程，采用DC-FFT算法來簡化加速計算。通過程序編程計算，分析驗證了所提出算法的精確度和高效性。

關鍵詞：三角形單元；接觸應力；DC-FFT；數值解；應力場

中圖分類號：TH123????????? 文獻標志碼：A????? 文章編號：1000-582X（2024）02-095-11

FFT acceleration algorithm for contact problems based on triangular element discretization

CHEN Nan1a， ZHU Kai1a， JIANG Zhizhen1a， GONG Shiyu1a， LI Pu2， JIN Xiaoqing1a，1b

（1a. College of Aerospace Engineering; 1b. State Key Laboratory of Mechanical Transmissions， Chongqing University， Chongqing 400044， P. R. China; 2. School of Science， Harbin Institute of Technology， Shenzhen 518055， Guangdong， P. R. China）

Abstract： Effectively solving the governing equations for contact problems often involves complex mathematical theory， while the distribution of contact stress is highly random in practical engineering applications. This study proposes a novel algorithm based on the triangular load discrete element and the discrete convolution fast Fourier transform （DC-FFT） algorithm. This algorithm provides a high-precision and reliable method for efficiently solving the contact response of a solid under any load distribution. Compared to the commonly used uniform load element discrete method， the analytical solution of the triangular element is more complex. However， it better simulates the characteristics of contact load distribution， accounting for situations where the load at the contact edge increases from zero or decreases to zero. The stress component under the action of the triangular and uniform load elements is derived based on the “excitation-response” characteristics of the contact influence coefficient matrix. This information is used to optimize the solution method of the triangular load discrete element. By constructing the stress solution in the form of a discrete convolution， including the influence coefficient matrix， the stress superposition effect of a target node under the action of all elements can be further simplified and accelerated by using the DC-FFT algorithm for highly repetitive matrix calculations. Programming and calculation analysis show that the proposed algorithm based on the triangular load element is accurate and efficient.

Keywords： triangular element; contact stress; DC-FFT; numerical solution; stress field

隨著機械裝備在高速度、高頻率和高精度等方面的需求越發苛刻，由局部壓力波動和微觀接觸形貌所造成的影響變得越發不可忽視[1]，由此引起的表界面接觸應力，進一步呈現出更加隨機和任意的分布特征。對上述接觸應力分布的高精度捕捉和分析，是揭示機械零部件失效和破壞的重要理論依據。然而，由于任意載荷分布情況下的封閉解析解很難獲得，數值計算方法成為求解此類問題的一種有效而通用的手段。早期對這類問題的研究是通過應用已知函數的無窮級數代替表面接觸壓力分布來求解的[2]，但是這種方法必須選擇正確的已知函數，否則會導致極大的誤差，而大多數情況下這些函數卻很難得到。Nowell等[3]利用三角形離散單元求解平面薄層彈性接觸問題的方法，顯示出了較高的計算精度，為任意分布載荷的接觸問題分析提供了一種有效思路，然而該方法消耗的時間較長。

在數值求解接觸問題時，由于在目標域邊界處容易產生周期性誤差，為了獲得更準確的結果，常常要以計算效率為代價，選擇遠大于目標尺寸的計算域，來彌補計算結果的誤差。基于快速傅里葉變換（fast Fourier transform，FFT）的接觸問題分析方法能大大減小計算負擔，是解決大規模復雜問題的有力工具。Ju等[4]首次將FFT引入到計算接觸問題的線性卷積過程中，他們利用傅里葉變換將接觸的時域問題轉換為傅里葉空間下的頻域問題，通過結合頻譜分析和FFT，為粗糙表面接觸的計算提供一種高效的技術。相關研究表明[5?6]，離散卷積快速傅里葉變換（discrete convolution fast Fourier transform，DC-FFT）算法在處理復雜的平面接觸問題時更便捷、更有效。Wang等[7]總結了涉及快速傅里葉變換（FFT）的不同算法，對不同的接觸問題、誤差控制，以及不同幾何形狀、材料和物理問題的求解方案進行了分析。Li等[8]利用FFT研究了無黏著接觸區內外的表面張力。基于FFT的數值方法也適用研究粗糙表面的微動接觸問題[9]。快速傅里葉變換也可以通過和Westergaard基本解結合，來有效地解決約束最小化問題，嚴格驗證接觸正交性[10]。Yu等[11]將FFT與共軛梯度法相結合，用于研究三維熱彈性滾動接觸的穩態模型。在研究微動接觸和微滑接觸時[12]，共軛梯度法和FFT的結合顯示出來更高的計算效率和準確性。在接觸球軸承的研究中[13]，FFT可以被用來求解橢圓接觸下，有油膜潤滑的軸承內部接觸載荷的分布。黏彈性問題的瞬態及穩態響應也可以通過FFT來提高計算效率[14]。

在復雜接觸應力情況下，為達到計算精度和計算消耗時間的優化，筆者提出用重疊的三角形單元離散任意函數形式分布的表面接觸應力，來求解受載固體內任意位置應力的方法。推導均布載荷單元和三角形載荷下的單元解，并進一步獲得卷積形式的應力公式，通過借助DC-FFT方法，實現加速計算。通過考慮赫茲分布形式的接觸應力，數值解與解析解顯示出較好的吻合度，驗證了所提出的數值計算方案的有效性。此外，通過控制單元寬度和單元數量，對比均布載荷單元離散與三角形單元離散計算出的應力場，通過對比分析2種離散方法的相對誤差，表明所提出的基于三角形載荷單元的快速傅立葉計算方法的優越性。最后，給出單元直接疊加方法與本文算法的效率對比，分析不同單元數量對CPU時間的影響。

1 任意載荷作用下彈性半平面模型

在彈性半平面接觸問題中，由于接觸問題的復雜性，可能會導致接觸表面載荷分布的隨意性，表面接觸區域（l，c）上的任意載荷可分解為法向載荷p（t）及切向載荷q（t），如圖1所示。

在表面距離O點t處取微元寬度dt，則點（x，z）處的應力分量可以表示為

{（σ_x=-2z/π ∫_l^c?〖（p（t）〖（x-t）〗^2）/（[〖（x-t）〗^2+z^2 ]^2 ） dt-〗? 2/π ∫_l^c?〖（q（t）〖（x-t）〗^3）/（[〖（x-t）〗^2+z^2 ]^2 ） dt ，〗@σ_z=-（2z^3）/π ∫_l^c?〖（p（t））/（[〖（x-t）〗^2+z^2 ]^2 ） dt-〗? （2z^2）/π ∫_l^c?〖（q（t）（x-t））/（[〖（x-t）〗^2+z^2 ]^2 ） dt 〗，@τ_xz=-（2z^2）/π ∫_l^c?〖（p（t）（x-t））/（[〖（x-t）〗^2+z^2 ]^2 ） dt-〗? 2z/π ∫_l^c?〖（q（t）〖（x-t）〗^2）/（[〖（x-t）〗^2+z^2 ]^2 ） dt 〗 。）┤? （1）

式（1）中包含奇異積分項。在載荷p（t），q（t）較復雜的情況下，各應力分量的解析求解較為困難，因此，通常采用數值方法進行離散求解。

2 數值方法

2.1 均布載荷單元的單元解

現行數值計算方法，通常采用分段均布載荷單元對載荷進行離散，當單元劃分得足夠精細時也可以較好地接近真實解。均布載荷單元載荷如圖2所示，載荷分布區域半寬為b，常數p、q分別為法向和切向載荷的大小，r_1 、r_2 、r分別為表面上（b，0）（-b，0）、O到響應點（x，z）的距離。

將載荷p代入公式（1），推導得到法向均布載荷單元作用下響應點的各應力分量為

{（σ_x=-p/π [arctan （b-x）/z+arctan （b+x）/z ├ -（z（b+x））/（r_1^2 ）+（z（x-b））/（r_2^2 ）]? ，┤@σ_z=-p/π [arctan （b-x）/z+arctan （b+x）/z ├ +（z（b+x））/（r_1^2 ）-（z（x-b））/（r_2^2 ）]? ，┤@τ_xz=-（z^2 p）/π （1/（r_1^2 ）-1/（r_2^2 ））? 。）┤??? （2）

其中，

{（r_1^2=〖（b-x）〗^2+z^2? ，@r_2^2=〖（b+x）〗^2+z^2? ，@r^2=x^2+z^2? 。）┤???? （3）

同理，代入q，得到切向均布單元載荷作用下響應點的各應力分量為

{（σ_x=-q/π （z^2/（r_2^2 ）-z^2/（r_1^2 ）+2ln r_2/r_1 ）? ，@σ_z=-（z^2 q）/π （1/（r_1^2 ）-1/（r_2^2 ））? ，@τ_xz=-q/π [arctan （b-x）/z+arctan （b+x）/z ├ -（z（b+x））/（r_1^2 ）+（z（x-b））/（r_2^2 ）]? 。┤ ）┤? （4）

利用均布單元離散任意載荷分布的接觸問題時，由于相鄰單元的載荷大小不盡相同，導致各個單元之間的連接是跳躍間斷的。應用該方法求解位移時，會出現單元交接處的位移梯度無窮大，進而得到的結果也不夠光滑連續，尤其是當接觸載荷從零開始變化的時候，任意載荷大小的均布單元也無法捕捉到接觸邊界上載荷為零的情況。

2.2 三角形載荷的單元解

為克服上述均布載荷單元的不足，將受載面上任意法向和切向載荷用一系列重疊等寬的三角形單元進行離散。相應的單個三角形分布的載荷單元如圖3所示。

由圖3可知，三角形單元的法向和切向載荷在-a～a的范圍內，從0線性增加至最大值，再線性減小為0。該三角形單元的半寬為a，單元載荷的公式為

（p（t）=p_0/a（a-|t|），|t|≤a? ，@q（t）=q_0/a（a-|t|），|t|≤a? 。） ???? （5）

用x表示半平面內響應點（x，z）與三角形單元底邊中點在t軸方向的相對距離，z表示響應點在z方向的深度，將表面上（-a，0）（a，0）、O到（x，z）的距離用d_1 、d_2 、d表示

{（d_1^2=〖（a-x）〗^2+z^2?? ， @d_2^2=〖（a+x）〗^2+z^2?? ，@d^2=x^2+z^2?? 。）┤??? （6）

將式（5）代入式（1），經過一系列積分運算與簡化，得三角形分布法向載荷p_0作用下各應力分量為

{（σ_x=-p_0/πa [（a-x）arctan （a-x）/z+（x+a）arctan （a+x）/z ├ -2xarctan x/z+2zln d^2/（d_1 d_2 ）]?? ，┤@σ_z=-p_0/πa [（x-a）arctan （x-a）/z┤+（x+a）arctan （a+x）/z ├ -2xarctan x/z]?? ，@τ_xz=（zp_0）/πa [arctan （a+x）/z+arctan （x-a）/z-2arctan x/z]?? 。）┤?? （7）

同理，三角形分布切向載荷q0作用下半平面內任意一點H（x，z）處各應力分量為

{（σ_x=q_0/π {2ln d_1/d_2 +2x/a ln d^2/（d_1 d_2 ）+3z/a （2arctan x/z ├ ├ -arctan （a+x）/z-arctan （a-x）/z） }?? ，┤ ┤@σ_z=（zq_0）/πa （arctan （a+x）/z+arctan （x-a）/z-2arctan x/z）?? ，@τ_xz=-q_0/πa [（a-x）arctan （a-x）/z+（x+a）arctan （a+x）/z ├ -2xarctan x/z+2zln d^2/（d_1 d_2 ）] ┤?? 。）┤? （8）

需要指出的是，式（7）～（8）與文獻[15]相應的公式不同，本文公式更加清晰地顯示出了場點和源點之間的作用關系，采用本文的形式可以更好地利用離散卷積的性質，便于后續的加速計算。

3 數值求解

3.1 三角形、均布載荷離散方法

對任意分布的載荷用一系列重疊等寬的三角形單元進行離散，如圖4所示。

將載荷劃分為n個單元，從左到右編號依次從1遞增到n，其第j個單元底邊中點位置為（x_j，0），用x_ij表示任意深度z=k的第j個三角形單元與響應點（x_i，z_k）在x方向的相對距離|x_i-x_j |，i表示計算區域內的第i個單元。在該法向單元載荷作用下，響應點處x方向的應力為

{σ_x }_ij=-p_j/πa [（a-x_ij）arctan （a-x_ij）/z_k ┤+（x_ij+a）arctan（a+x_ij）/z_k? ├ -2xarctan x_ij/z_k +2z_k ln d^2/（d_1 d_2 ）]，?? （9）

式中，

{（d_1^2=（a-x_ij ）^2+z_k^2?? ，@d_2^2=（a+x_ij ）^2+z_k^2?? ，@d^2=x_ij^2+z_k^2?? 。）┤? （10）

取應力函數的影響系數T_（j-i）為

T_（j-i）=-1/πa [（a-x_ij）arctan （a-x_ij）/z_k +（x_ij+a）arctan （a+x_ij）/z_k ┤ ├ -2xarctan x_ij/z_k +2z_k ln d^2/（d_1 d_2 ）] 。??? （11）

受載半平面內任意位置點（x_i，z_k）在n個三角形載荷單元的疊加作用下，其x方向的應力可表示為：

σ_（x_i ）=∑_（j=1）^n?T_（j-i）? p_j 。? （12）

同理，其他方向以及切向力作用下的應力均可表示為相同的形式。

均布載荷單元離散如圖5所示，疊加思想過程與三角形離散過程相同。

3.2 離散卷積快速傅里葉變換求解

本文研究的問題可以利用“激勵-響應”機制求解，彈性基體次表層的各應力分布，都是表面上n個三角形單元疊加作用下的響應，可以寫成如下卷積形式，其中*表示卷積

（f *g）（x）=f（x）*g（x）=∫_（-∞）^∞?〖f（α）〗 g（x-α）dα? 。???? （13）

對于所有變量都僅與相對距離相關的連續系統，可以利用快速傅里葉變換來解決計算規模龐大的復雜問題。為方便編程計算，將半平面離散為合適尺寸的均布載荷單元，如圖5所示。在載荷作用范圍內，深度為z的平面上取相同數量n個均布載荷區域的節點，利用式（12），該n個節點的應力值可表示為如下應力影響矩陣T_（j-i）與載荷矩陣P_j的乘積，即

[（σ_（x_1 ）@σ_（x_2 ）@?@σ_（x_n ） ）]=[（T_0&T_（-1）&…&T_（1-n）@T_1&T_0&…&T_（2-n）@?&?&?&?@T_（n-1）&T_（n-2）&…&T_0 ）][（p_1@p_2@?@p_（n-1）@p_n ）] 。????? （14）

式（14）影響矩陣為n階Toeplitz矩陣，為引入FFT算法，需將影響矩陣擴充為循環矩陣。現取其首列、首行來構造包卷循環a_2n，

a_2n=[（T_0&T_1&…&T_（n-1）&（0&T_（1-n）&…）&T_（-1） ）]^T 。???? （15）

將這個包卷循環作為循環矩陣C_2n的第一列，利用循環矩陣的性質來將矩陣補充完整。載荷矩陣P_n通過在第n項后補n個0，將其擴充成為2n項，得到P_2n。則擴充后的應力分量表示為

σ_2n=C_2n P_2n=C_2n [（P_n@0_n ）]=[（T_（j-i） P_n@Θ）]=[（σ_n@Θ）] ，? （16）

式中，Θ表示多余項，與求解無關，原本需要求解的應力矩陣就是擴充運算后σ_2n的前n項。利用循環矩陣如下性質

C_2n=F_2n^（-1） diag（F_2n a_2n）F_2n 。????? （17）

式中，F_2n是離散傅里葉矩陣，diag函數用于構造一個對角矩陣，式（17）代入式（16）得

σ_2n=C_2n P_2n=F_2n^（-1） diag（F_2n a_2n）F_2n P_2n 。???? （18）

式（18）兩邊左乘離散傅里葉矩陣得

F_2n σ_2n=F_2n F_2n^（-1） diag（F_2n a_2n）F_2n P_2n=diag（F_2n a_2n）F_2n P_2n? 。????? （19）

由式（19）可知，求解擴展后的應力矩陣時只需對式（15）以及載荷矩陣P_2n做傅里葉變換后對應項相乘，再對整體做一次逆變換即可：

σ_2n=F_2n^（-1） [diag（F_2n a_2n）F_2n P_2n] 。? （20）

最終取σ_2n的前n項即所求節點的應力矩陣。引入離散傅里葉矩陣可以細化運算結構，把原始高階矩陣運算，依次分解成一系列低階矩陣運算。充分利用離散傅里葉變換所具有的對稱性質和周期性質，并進行適當組合，就可以去除重復計算，減少計算過程中的乘法運算，讓整體的運算結構更清晰簡便，這就是快速傅里葉變換應用在此類“激勵-響應”系統中起到加速效果的機理。而且在具體運用中涉及到的計算規模越大，FFT算法在計算效率上的優越性就越顯著。

4 結果與討論

4.1 三角形單元離散數值解驗證

為驗證三角形單元離散數值解法的有效性，現引入赫茲型的接觸載荷，如式（21）所示，

P（x）=P_0 √（1-x^2/r^2 ） 。? （21）

切向載荷大小是法向載荷的0.3倍，載荷半徑r為5 mm，中心軸為z軸，P_0為100 MPa。第j個三角形離散單元的頂點處對應的法向載荷為

P（j）=P_0 √（1-〖x_j〗^2/r^2 ） 。???? （22）

將接觸應力分為40個離散單元，每個單元半寬a=0.25 mm，將受載半平面用0.25 mm×0.25 mm的正方形單元進行離散，如圖6所示。

利用Fortran編程計算，將赫茲型載荷解析解得到的結果與相同接觸應力作用下三角形單元離散方法得到的各應力分量進行對比驗證。

圖7是法向力作用下深度為1 mm時，x從-4 mm到4 mm各節點上三角形離散數值解與解析解各應力分量對比。圖7中，3條曲線分別是通過解析法計算得出的σ_xx 、σ_zz 、σ_xz值，相同顏色散點是三角形單元離散數值方法計算得出的對應的應力值。圖8是切向力作用下2種方法得到的各應力值對比。觀察可知，此種方法得到的數值解與解析解無論在法向還是切向接觸應力的作用下，深度為1 mm處各應力分量都能保持較高的吻合度。

為驗證本文數值算法的有效性，對比驗證了受法向與切向載荷共同作用且接觸平面下深度z分別為0.25、1.25、2.25 mm時，解析解和數值解得到的各節點應力響應，如圖9～11所示。

在法向與切向載荷共同作用下，三角形離散數值方法求得的不同深度各應力分量（散點）數值解也可以與解析解（曲線）保持很好的吻合度。圖10中z方向應力值在z=1.25時，偏差相對較大，經計算，該組數據中數值解與解析解的平均相對誤差實際只有0.123%，這并不影響本數值方法的有效性。經計算，在總寬10 mm，單元寬度為0.25 mm的情況下，受載半平面內各點的應力分量數值解和解析解的相對誤差平均值僅為0.048 7%。因此，利用一系列重疊等寬的三角形單元離散表面接觸應力，進行數值計算是可行的。

4.2 三角形單元離散與均布載荷單元離散數值解對比

利用3.1節中提到的數值計算方法，用均布載荷單元離散法向和切向接觸應力，通過Fortran編程計算，得到受載半平面內各節點應力場數值解。將接觸平面深度z=0.25 mm處的各節點均布載荷單元離散數值解、三角形單元離散數值解與赫茲解析解的相對誤差作對比，如圖12～13所示。

由圖12～13觀察可知，在相同的單元寬度和單元數量下，用三角形單元進行離散比用均布載荷單元離散得到的各應力分量數值解與解析解的相對誤差更小。為更全面地觀察2種離散方法的準確性，用各應力分量換算得到Mises屈服應力，將數值方法得到的在接觸面下z=1 mm時，x從-3 mm到3 mm的一系列節點上Mises屈服應力相對解析解的誤差值進行對比，得到結果如表1所示。

由表1分析可知，均布載荷離散數值解與相同位置Mises應力解析解的相對誤差比三角形離散數值解與相同位置Mises應力解析解的相對誤差大一倍。這驗證了三角形單元離散的方法由于各個單元之間不存在跳躍式的間斷，而是相互重疊，使單元與單元之間過渡更加平穩，與原本光滑連續的接觸應力實際作用情況更為接近，所以在計算精度上要優于普通的均布載荷離散方法。

4.3 DC-FFT與直接計算效率對比

在數值離散方法中，離散單元越精細，數量越多，就越會得到更接近真實解的結果，但這也會帶來極大的計算量，所以采用更為高效的計算方法來提高運算效率非常必要。DC-FFT方法可以將原本N階離散傅里葉變換的時間復雜度從O（N^2）降到O（NlogN），其中О表示算法復雜度的階次范疇，運算的復雜度越高，DC-FFT方法節約的運算成本就越大。為了對比本文算法和直接計算方法在計算效率上的具體差距，分別用這2種方法求解法向載荷作用下半平面深度為1 mm的一系列節點的x方向應力，通過改變計算的節點數量，對比2種計算方法所需CPU時間，得到結果如表2所示。

由表2可知，計算2^13個節點的應力值時，2種方法計算耗時相差184倍。隨著計算量增加，可以看到，節點數從2^13個增加到2^18個的過程中，直接計算所需要的CPU時間急劇增加，而用FFT計算所需的時間卻始終不到0.1 s，在計算量增加到2^18個節點時，DC-FFT方法所需時間剛好超過0.1 s，而直接計算消耗的計算時間已經接近1 h，從比率上看兩者時間更是相差34 957倍。以這樣的趨勢和DC-FFT方法的特性上看，在需要精確計算基體更大區域范圍內各應力分量時，由于計算規模急劇增加，2種方法所消耗的計算資源差距也會更加明顯。因此，DC-FFT方法是能穩定高效地得到相當高精度數值解的優質算法。

5 結? 論

提出了一種高效的數值化計算方法，用以快速求解復雜接觸應力情況下基體內部應力場。主要得到如下結論：

1）分別獲得了在三角形和分段均布接觸載荷分布下，基體內任意點應力的顯式解析解；

2）通過與快速傅里葉變換算法相結合，分別構建了基于三角形和分段均布接觸載荷單元的高效計算方案；

3）相比于分段均布載荷的矩形離散算法，在相同單元數和單元寬度下，基于三角形接觸載荷單元的離散算法得到的數值解更精確。

4）基于三角形接觸載荷單元和離散卷積-快速傅里葉變換算法（DC-FFT）所構建的計算方案，具有優異的魯棒性、計算效率和計算精度。

參考文獻

［1］? Li Y Y， Yang Y， Li M， et al. Dynamics analysis and wear prediction of rigid-flexible coupling deployable solar array system with clearance joints considering solid lubrication[J]. Mechanical Systems and Signal Processing， 2022， 162： 108059.

［2］? Brebbia C A. The boundary element method in engineering practice[J]. Engineering Analysis， 1984， 1（1）： 3-12.

［3］? Nowell D， Hills D A. Tractive rolling of tyred cylinders[J]. International Journal of Mechanical Sciences， 1988， 30（12）： 945-957.

［4］? Ju Y， Farris T N. Spectral analysis of two-dimensional contact problems[J]. Journal of Tribology， 1996， 118（2）： 320-328.

［5］? Sun L L， Wang Q J， Zhao N， et al. Discrete convolution and FFT modified with double influence-coefficient superpositions （DCSS-FFT） for contact of nominally flat heterogeneous materials involving elastoplasticity[J]. Computational Mechanics， 2021， 67（3）： 989-1007.

［6］? Sun L L， Wang Q J， Zhang M Q， et al. Discrete convolution and FFT method with summation of influence coefficients （DCS-FFT） for three-dimensional contact of inhomogeneous materials[J]. Computational Mechanics， 2020， 65（6）： 1509-1529.

［7］? Wang Q J， Sun L L， Zhang X， et al. FFT-based methods for computational contact mechanics[J]. Frontiers in Mechanical Engineering， 2020， 6： 61.

［8］? Li Q A， Popov V L. Non-adhesive contacts with different surface tension inside and outside the contact area[J]. Frontiers in Mechanical Engineering， 2020， 6： 63.

［9］? 牛睿， 萬強， 靳凡， 等. 基于共軛梯度法和快速傅立葉變換的粗糙表面微動接觸數值研究[J]. 計算力學學報， 2017， 34（3）： 312-321.

Niu R， Wan Q， Jin F， et al. A numerical analysis of fretting contact with rough surface based on conjugate gradient method and fast Fourier transform[J]. Chinese Journal of Computational Mechanics， 2017， 34（3）： 312-321.（in Chinese）

［10］? Rey V， Anciaux G， Molinari J F. Normal adhesive contact on rough surfaces： efficient algorithm for FFT-based BEM resolution[J]. Computational Mechanics， 2017， 60（1）： 69-81.

［11］? Yu Y H， Suh J. Numerical analysis of three-dimensional thermo-elastic rolling contact under steady-state conditions[J]. Friction， 2022， 10（4）： 630-644.

［12］? 吳桐. 粗糙表面微動接觸分析[D]. 武漢： 武漢科技大學， 2019.

Wu T. Fretting contact analysis of rough surfaces[D]. Wuhan： Wuhan University of Science and Technology， 2019.（in Chinese）

［13］? 帥琪琪， 陳曉陽， 陳世金， 等. 預緊方式對彈流潤滑下角接觸球軸承內部力學特性的影響[J]. 摩擦學學報， 2022， 42（1）： 85-94.

Shuai Q Q， Chen X Y， Chen S J， et al. Influences of preload methods on internal mechanical characteristics of angular contact ball bearings under elastohydrodynamic lubrication[J]. Tribology， 2022， 42（1）： 85-94.（in Chinese）

［14］? Zhang X， Wang Q J， He T. Transient and steady-state viscoelastic contact responses of layer-substrate systems with interfacial imperfections[J]. Journal of the Mechanics and Physics of Solids， 2020， 145： 104170.

［15］? Johnson K L. Contact mechanics[M]. Cambridge： Cambridge University Press， 1987.

（編輯? 鄭潔）