陳楠 朱凱 蔣志楨 龔詩(shī)雨 李璞 金曉清

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

關(guān)鍵詞：三角形單元；接觸應(yīng)力；DC-FFT；數(shù)值解；應(yīng)力場(chǎng)

中圖分類(lèi)號(hào)：TH123????????? 文獻(xiàn)標(biāo)志碼：A????? 文章編號(hào)：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

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

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

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

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

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

在表面距離O點(diǎn)t處取微元寬度dt，則點(diǎn)（x，z）處的應(yīng)力分量可以表示為

{（σ_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）中包含奇異積分項(xiàng)。在載荷p（t），q（t）較復(fù)雜的情況下，各應(yīng)力分量的解析求解較為困難，因此，通常采用數(shù)值方法進(jìn)行離散求解。

2 數(shù)值方法

2.1 均布載荷單元的單元解

現(xiàn)行數(shù)值計(jì)算方法，通常采用分段均布載荷單元對(duì)載荷進(jìn)行離散，當(dāng)單元?jiǎng)澐值米銐蚓?xì)時(shí)也可以較好地接近真實(shí)解。均布載荷單元載荷如圖2所示，載荷分布區(qū)域半寬為b，常數(shù)p、q分別為法向和切向載荷的大小，r_1 、r_2 、r分別為表面上（b，0）（-b，0）、O到響應(yīng)點(diǎn)（x，z）的距離。

將載荷p代入公式（1），推導(dǎo)得到法向均布載荷單元作用下響應(yīng)點(diǎn)的各應(yīng)力分量為

{（σ_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，得到切向均布單元載荷作用下響應(yīng)點(diǎn)的各應(yīng)力分量為

{（σ_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）

利用均布單元離散任意載荷分布的接觸問(wèn)題時(shí)，由于相鄰單元的載荷大小不盡相同，導(dǎo)致各個(gè)單元之間的連接是跳躍間斷的。應(yīng)用該方法求解位移時(shí)，會(huì)出現(xiàn)單元交接處的位移梯度無(wú)窮大，進(jìn)而得到的結(jié)果也不夠光滑連續(xù)，尤其是當(dāng)接觸載荷從零開(kāi)始變化的時(shí)候，任意載荷大小的均布單元也無(wú)法捕捉到接觸邊界上載荷為零的情況。

2.2 三角形載荷的單元解

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

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

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

用x表示半平面內(nèi)響應(yīng)點(diǎn)（x，z）與三角形單元底邊中點(diǎn)在t軸方向的相對(duì)距離，z表示響應(yīng)點(diǎn)在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），經(jīng)過(guò)一系列積分運(yùn)算與簡(jiǎn)化，得三角形分布法向載荷p_0作用下各應(yīng)力分量為

{（σ_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作用下半平面內(nèi)任意一點(diǎn)H（x，z）處各應(yīng)力分量為

{（σ_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）與文獻(xiàn)[15]相應(yīng)的公式不同，本文公式更加清晰地顯示出了場(chǎng)點(diǎn)和源點(diǎn)之間的作用關(guān)系，采用本文的形式可以更好地利用離散卷積的性質(zhì)，便于后續(xù)的加速計(jì)算。

3 數(shù)值求解

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

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

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

{σ_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）

取應(yīng)力函數(shù)的影響系數(shù)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）

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

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

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

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

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

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

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

對(duì)于所有變量都僅與相對(duì)距離相關(guān)的連續(xù)系統(tǒng)，可以利用快速傅里葉變換來(lái)解決計(jì)算規(guī)模龐大的復(fù)雜問(wèn)題。為方便編程計(jì)算，將半平面離散為合適尺寸的均布載荷單元，如圖5所示。在載荷作用范圍內(nèi)，深度為z的平面上取相同數(shù)量n個(gè)均布載荷區(qū)域的節(jié)點(diǎn)，利用式（12），該n個(gè)節(jié)點(diǎn)的應(yīng)力值可表示為如下應(yīng)力影響矩陣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算法，需將影響矩陣擴(kuò)充為循環(huán)矩陣?，F(xiàn)取其首列、首行來(lái)構(gòu)造包卷循環(huán)a_2n，

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

將這個(gè)包卷循環(huán)作為循環(huán)矩陣C_2n的第一列，利用循環(huán)矩陣的性質(zhì)來(lái)將矩陣補(bǔ)充完整。載荷矩陣P_n通過(guò)在第n項(xiàng)后補(bǔ)n個(gè)0，將其擴(kuò)充成為2n項(xiàng)，得到P_2n。則擴(kuò)充后的應(yīng)力分量表示為

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

式中，Θ表示多余項(xiàng)，與求解無(wú)關(guān)，原本需要求解的應(yīng)力矩陣就是擴(kuò)充運(yùn)算后σ_2n的前n項(xiàng)。利用循環(huán)矩陣如下性質(zhì)

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

式中，F(xiàn)_2n是離散傅里葉矩陣，diag函數(shù)用于構(gòu)造一個(gè)對(duì)角矩陣，式（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）可知，求解擴(kuò)展后的應(yīng)力矩陣時(shí)只需對(duì)式（15）以及載荷矩陣P_2n做傅里葉變換后對(duì)應(yīng)項(xiàng)相乘，再對(duì)整體做一次逆變換即可：

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

最終取σ_2n的前n項(xiàng)即所求節(jié)點(diǎn)的應(yīng)力矩陣。引入離散傅里葉矩陣可以細(xì)化運(yùn)算結(jié)構(gòu)，把原始高階矩陣運(yùn)算，依次分解成一系列低階矩陣運(yùn)算。充分利用離散傅里葉變換所具有的對(duì)稱(chēng)性質(zhì)和周期性質(zhì)，并進(jìn)行適當(dāng)組合，就可以去除重復(fù)計(jì)算，減少計(jì)算過(guò)程中的乘法運(yùn)算，讓整體的運(yùn)算結(jié)構(gòu)更清晰簡(jiǎn)便，這就是快速傅里葉變換應(yīng)用在此類(lèi)“激勵(lì)-響應(yīng)”系統(tǒng)中起到加速效果的機(jī)理。而且在具體運(yùn)用中涉及到的計(jì)算規(guī)模越大，F(xiàn)FT算法在計(jì)算效率上的優(yōu)越性就越顯著。

4 結(jié)果與討論

4.1 三角形單元離散數(shù)值解驗(yàn)證

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

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

切向載荷大小是法向載荷的0.3倍，載荷半徑r為5 mm，中心軸為z軸，P_0為100 MPa。第j個(gè)三角形離散單元的頂點(diǎn)處對(duì)應(yīng)的法向載荷為

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

將接觸應(yīng)力分為40個(gè)離散單元，每個(gè)單元半寬a=0.25 mm，將受載半平面用0.25 mm×0.25 mm的正方形單元進(jìn)行離散，如圖6所示。

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

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

為驗(yàn)證本文數(shù)值算法的有效性，對(duì)比驗(yàn)證了受法向與切向載荷共同作用且接觸平面下深度z分別為0.25、1.25、2.25 mm時(shí)，解析解和數(shù)值解得到的各節(jié)點(diǎn)應(yīng)力響應(yīng)，如圖9～11所示。

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

4.2 三角形單元離散與均布載荷單元離散數(shù)值解對(duì)比

利用3.1節(jié)中提到的數(shù)值計(jì)算方法，用均布載荷單元離散法向和切向接觸應(yīng)力，通過(guò)Fortran編程計(jì)算，得到受載半平面內(nèi)各節(jié)點(diǎn)應(yīng)力場(chǎng)數(shù)值解。將接觸平面深度z=0.25 mm處的各節(jié)點(diǎn)均布載荷單元離散數(shù)值解、三角形單元離散數(shù)值解與赫茲解析解的相對(duì)誤差作對(duì)比，如圖12～13所示。

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

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

4.3 DC-FFT與直接計(jì)算效率對(duì)比

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

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

5 結(jié)? 論

提出了一種高效的數(shù)值化計(jì)算方法，用以快速求解復(fù)雜接觸應(yīng)力情況下基體內(nèi)部應(yīng)力場(chǎng)。主要得到如下結(jié)論：

1）分別獲得了在三角形和分段均布接觸載荷分布下，基體內(nèi)任意點(diǎn)應(yīng)力的顯式解析解；

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

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

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

參考文獻(xiàn)

［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］? 牛睿， 萬(wàn)強(qiáng)， 靳凡， 等. 基于共軛梯度法和快速傅立葉變換的粗糙表面微動(dòng)接觸數(shù)值研究[J]. 計(jì)算力學(xué)學(xué)報(bào)， 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òng)接觸分析[D]. 武漢： 武漢科技大學(xué)， 2019.

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

［13］? 帥琪琪， 陳曉陽(yáng)， 陳世金， 等. 預(yù)緊方式對(duì)彈流潤(rùn)滑下角接觸球軸承內(nèi)部力學(xué)特性的影響[J]. 摩擦學(xué)學(xué)報(bào)， 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.

（編輯? 鄭潔）