偶应力弹性半空间轴对称完全粘着接触

Axisymmetric full adhesive contact of couple-stress elastic half-space

摘要

随着高新技术产业的发展，材料和结构在微纳尺度下的力学性能受到研究者的广泛关注。当构件的接触尺寸由宏观缩小到微观时，研究尺度效应对接触力学行为的影响将变得非常重要。通过理论分析的方法，研究了偶应力弹性半空间在球形压头作用下的完全粘着接触问题。通过引入材料特征长度参数来描述尺度效应，借助Hankel积分变换将力学模型转化为奇异积分方程，并利用数值方法对其进行求解。结果表明，材料特征长度参数对完全粘着接触时的应力分布、压力-压痕曲线有显著影响。可通过改变材料的材料特征长度参数来调节材料的力学性能，为材料的优化设计和抑制接触损伤提供可靠的理论依据。

Abstract

With the development of high-tech industry, the mechanical properties of materials and structures at micro-nano scale have been widely concerned. When the contact size of the component is reduced from macroscopic to microscopic, the study of the size effect on the contact mechanical behavior becomes crutia. This paper investigates the full adhesive contact of couple-stress elastic half-space under spherical indenter by means of theoretical analysis. The material characteristic length parameter was introduced to describe the size effect. The adhesive contact model was transformed into a singular integral equation which was solved by numerical method by using the Hankel integral transformation. The results show that the material characteristic length parameters have significant effects on the stress distribution and pressure-indentation curve of full adhesive contact. The mechanical properties of the material can be adjusted by changing the material characteristic length parameter, which provides a reliable theoretical basis for the optimal design of the material and the suppression of contact damage.

Graphical abstract

关键词

偶应力弹性体 / 完全粘着 / 尺度效应 / Hankel积分变换 / 奇异积分方程

Key words

couple-stress elastomer / full adhesion / size effect / Hankel integral transformation / singular integral equation

随着高新技术产业的发展，微电机系统受到了研究者们的关注。由于微机电系统下的材料和结构尺寸达到了微米或纳米级，因此，研究构件在微观尺度下的力学性能变得非常重要。接触问题主要研究外载荷作用下的两个接触体相互作用所产生的应力和应变，最早是由Hertz^[1]提出了两个宏观弹性体之间的接触Hertz解，后来Hertz解被广泛应用于岩土、机械、材料等领域，以此来获得表面接触应力分布，但其只适用于宏观弹性体之间的无摩擦接触问题。随着材料尺寸减小到了微米甚至纳米级，经典的Hertz解无法描述尺度效应。后来Mindlin^[2]建立了偶应力理论，并引入了表征材料微观结构尺寸的材料特征长度参数来描述尺度效应。在此基础上，有大量学者利用偶应力理论来研究弹性体接触的尺度依赖接触问题^[3]。Zisis等^[4]基于偶应力理论得出了二维接触问题的一般解。Song等^[5]分析了基于偶应力理论的二维均匀半平面滑动摩擦接触问题，分析了尺度效应和摩擦系数对接触力学行为的影响，其结果表明，在微纳尺度下材料表现出很强的尺度效应。Li等^[6-8]分析了二维均匀半平面黏附接触问题、二维梯度纳米涂层均匀半平面黏附接触问题以及二维梯度纳米涂层完全黏附接触问题，其结果都表明，材料特征长度对材料的力学性能有显著影响。由于二维问题存在边界效应且无法满足实验要求，一些学者将微纳尺度下的接触问题扩展到了三维。Gourgiotis等^[9]证明了引入材料特征长度参数的偶应力理论对于复杂微结构材料的建模是非常有效且相当简单的，并通过积分变换和奇异积分方程，得到了轴对称接触问题的精确解，还分析了偶应力三维轴对称均匀半空间无摩擦接触问题，证明偶应力理论可以在微观结构中捕捉到并实验观察到压痕响应和相关尺度效应。Li等^[10]考虑了偶应力下的轴对称均匀半空间无摩擦黏附接触问题，其结果表明，在考虑偶应力时，材料表现出很强的尺度效应。之前的这些研究结果，关于偶应力弹性体的完全粘着接触问题还未考虑到，因此，本文将基于偶应力理论，建立偶应力弹性半空间在球形压头作用下的完全粘着接触力学模型，借助Hankel积分变换^[11]，将偶应力弹性半空间轴对称完全粘着接触力学模型转化为奇异积分方程并进行数值求解，对材料特征长度参数完全粘着接触下的力学响应进行分析。

1 偶应力弹性半空间在球形压头作用下的完全粘着接触力学模型

如图1所示，在法向载荷P作用下，半径为R的球形压头与偶应力弹性半空间发生接触，形成半径为a的接触区。建立柱坐标系，其中，r坐标沿水平方向向右，z坐标竖直向下。假设球形压头与弹性半空间接触区是完全粘着的，即不发生相对滑动。

由文献[9]可知，基于偶应力理论，柱坐标下的弹性体几何方程如下：

{\epsilon }_{rr}=\frac{\partial u_r}{\partial r},{\epsilon }_{\theta \theta }=\frac{u_r}{r},{\epsilon }_{zz}=\frac{\partial u_z}{\partial z}(1)

{\epsilon }_{zr}=\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\partial u_z}{\partial r}+\frac{\partial u_r}{\partial z}\right)={\epsilon }_{rz}(2)

{\omega }_{\theta }=\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\partial u_r}{\partial z}-\frac{\partial u_z}{r}\right)(3)

{\kappa }_{\theta r}=-\frac{{\omega }_{\theta }}{r},{\kappa }_{r\theta }=\frac{\partial {\omega }_{\theta }}{\partial r},{\kappa }_{z\theta }=\frac{\partial {\omega }_{\theta }}{\partial z}(4)

其中，{\epsilon }_{rr}\mathrm{、}{\epsilon }_{\theta \theta }\mathrm{、}{\epsilon }_{zz}\mathrm{、}{\epsilon }_{zr}\mathrm{、}{\epsilon }_{rz}为应变分量，u_r\mathrm{、}{u}_z为位移分量，{\omega }_{\theta }为旋转矢量，{\kappa }_{\theta r}\mathrm{、}{\kappa }_{r\theta }\mathrm{、}{\kappa }_{z\theta }为曲率张量。

平衡方程如下：

\frac{\partial {\sigma }_{rr}}{\partial r}+\frac{\partial {\sigma }_{zr}}{\partial z}+\frac{{\sigma }_{rr}-{\sigma }_{\theta \theta }}{r}=\mathrm{0}(5)

\frac{\partial {\sigma }_{rz}}{\partial r}+\frac{\partial {\sigma }_{zz}}{\partial z}+\frac{{\sigma }_{rz}}{r}=\mathrm{0}(6)

\frac{\partial m_{r\theta }}{\partial r}+\frac{\partial m_{z\theta }}{\partial z}+\frac{m_{r\theta }+m_{\theta r}}{r}+{\sigma }_{zr}-{\sigma }_{rz}=\mathrm{0}(7)

其中，{\sigma }_{rr}\mathrm{、}{\sigma }_{\theta \theta }\mathrm{、}{\sigma }_{zz}\mathrm{、}{\sigma }_{zr}\mathrm{、}{\sigma }_{rz}为非对称的应力分量，m_{r\theta }\mathrm{、}{m}_{z\theta }\mathrm{、}{m}_{\theta r}为非对称偶应力分量。

物理方程如下：

{\sigma }_{rr}=\frac{\mathrm{2}\mu v}{\mathrm{1}-\mathrm{2}v}\left({\epsilon }_{rr}+{\epsilon }_{\theta \theta }+{\epsilon }_{zz}\right)+\mathrm{2}\mu {\epsilon }_{rr}(8)

{\sigma }_{\theta \theta }=\frac{\mathrm{2}\mu v}{\mathrm{1}-\mathrm{2}v}\left({\epsilon }_{rr}+{\epsilon }_{\theta \theta }+{\epsilon }_{zz}\right)+\mathrm{2}\mu {\epsilon }_{\theta \theta }(9)

{\sigma }_{zz}=\frac{\mathrm{2}\mu v}{\mathrm{1}-\mathrm{2}v}\left({\epsilon }_{rr}+{\epsilon }_{\theta \theta }+{\epsilon }_{zz}\right)+\mathrm{2}\mu {\epsilon }_{zz}(10)

\left\{\begin{array}{c}{\sigma }_{zr}={\tau }_{zr}+{\alpha }_{zr}\\ {\sigma }_{rz}={\tau }_{rz}+{\alpha }_{rz}\\ {\sigma }_{rz}={\tau }_{rz}+{\alpha }_{rz}\end{array}\right.(11)

{\alpha }_{zr}=-\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\partial m_{r\theta }}{\partial r}+\frac{m_{r\theta }}{r}+\frac{m_{\theta r}}{r}+\frac{\partial m_{z\theta }}{\partial z}\right)=-{\alpha }_{rz}(12)

m_{r\theta }=\mathrm{4}\mu l^{\mathrm{2}}\left({\kappa }_{r\theta }+\beta {\kappa }_{\theta r}\right)(13)

m_{\theta r}=\mathrm{4}\mu l^{\mathrm{2}}\left({\kappa }_{\theta r}+\beta {\kappa }_{r\theta }\right)(14)

m_{z\theta }=\mathrm{4}\mu l^{\mathrm{2}}{\kappa }_{z\theta },m_{\theta z}=\mathrm{4}\mu l^{\mathrm{2}}\beta {\kappa }_{z\theta }(15)

其中，{\tau }_{zr}\mathrm{、}{\tau }_{rz}\mathrm{、}{\alpha }_{zr}\mathrm{、}{\alpha }_{rz}为切应力分量的对称和反对称的部分，\mu \mathrm{、}\nu分别为剪切模量和泊松比，\beta为一个无量纲的参数，，l为偶应力理论中表征材料尺度效应的材料特征长度参数。

将几何方程(1)~(4)和物理方程(8)~(15)代入平衡方程(5)~(7)，得到用位移分量表示的控制方程如下：

\begin{array}{l}\left({\nabla }^{\mathrm{2}}-\frac{\mathrm{1}}{r^{\mathrm{2}}}\right)u_r+\left(\lambda +\mu \right)\frac{\partial }{\partial r}\left(\frac{\partial u_r}{\partial r}+\frac{\partial u_z}{\partial z}+\frac{u_r}{r}\right)\\ -\mu l^{\mathrm{2}}\left({\nabla }^{\mathrm{2}}-\frac{\mathrm{1}}{r^{\mathrm{2}}}\right)\left(\frac{{\partial }^{\mathrm{2}}{u}_r}{\partial z^{\mathrm{2}}}-\frac{{\partial }^{\mathrm{2}}{u}_z}{\partial z\partial r}\right)=\mathrm{0}\end{array}(16)

\begin{array}{c}{\nabla }^{\mathrm{2}}{u}_z+\left(\lambda +\mu \right)\frac{\partial }{\partial z}\left(\frac{\partial u_r}{\partial r}+\frac{\partial u_z}{\partial z}+\frac{u_r}{r}\right)\\ -\mu l^{\mathrm{2}}{\nabla }^{\mathrm{2}}\left(\frac{{\partial }^{\mathrm{2}}{u}_z}{\partial r^{\mathrm{2}}}+\frac{\mathrm{1}}{r}\left(\frac{\partial u_z}{\partial r}-\frac{\partial u_r}{\partial z}\right)-\frac{{\partial }^{\mathrm{2}}{u}_z}{\partial z\partial r}\right)=\mathrm{0}\end{array}(17)

其中，{\nabla }^{\mathrm{2}}为拉普拉斯算子。

Mindlin等^[12]给出了上述方程组的通解如下：

u_r=-\frac{\partial }{\partial r}\left\{l^{\mathrm{2}}\frac{\partial \Psi }{\partial z}+\alpha \left[z\left(\mathrm{1}-l^{\mathrm{2}}{\nabla }^{\mathrm{2}}\right)\Psi +\Phi \right]\right\}(18)

u_z=\Psi -\frac{\partial }{\partial z}\left\{l^{\mathrm{2}}\frac{\partial \Psi }{\partial z}+\alpha \left[z\left(\mathrm{1}-l^{\mathrm{2}}{\nabla }^{\mathrm{2}}\right)\Psi +\Phi \right]\right\}(19)

\alpha =\frac{\lambda +\mu }{\mathrm{2}\left(\lambda +\mathrm{2}\mu \right)}=\frac{\mathrm{1}}{\mathrm{4}\left(\mathrm{1}-v\right)}(20)

其中，\lambda为第一拉梅常数。

\Psi和\Phi是r和z的函数，且满足如下方程：

{\nabla }^{\mathrm{2}}\Phi =\mathrm{0}(21)

\left(\mathrm{1}-l^{\mathrm{2}}{\nabla }^{\mathrm{2}}\right){\nabla }^{\mathrm{2}}\Psi =\mathrm{0}(22)

式(21)和式(22)可利用Hankel变换将其变形。

Hankel变换和Hankel逆变换^[11]公式如下：

{〈F\left(s\right)〉}_n={\int }_{\mathrm{0}}^{\mathrm{\infty }}rF\left(x\right)J_n\left(sr\right)\mathrm{d}r(23)

F\left(x\right)={\int }_{\mathrm{0}}^{\mathrm{\infty }}s{〈F〉}_n{J}_n\left(sr\right)\mathrm{d}s(24)

其中，J_n(\cdot )表示n阶第一类贝塞尔函数，{〈\cdot 〉}_n表示n阶Hankel积分变换，s为积分变量。

对式(21)和式(22)进行Hankel积分变换，得到方程如下：

\left(\frac{{\mathrm{d}}^{\mathrm{2}}}{\mathrm{d}{z}^{\mathrm{2}}}-s^{\mathrm{2}}\right){〈\Phi 〉}_{\mathrm{0}}=\mathrm{0}(25)

\left[\mathrm{1}-l^{\mathrm{2}}\left(\frac{{\mathrm{d}}^{\mathrm{2}}}{\mathrm{d}{z}^{\mathrm{2}}}-s^{\mathrm{2}}\right)\right]\left(\frac{{\mathrm{d}}^{\mathrm{2}}}{\mathrm{d}{z}^{\mathrm{2}}}-s^{\mathrm{2}}\right){〈\Psi 〉}_{\mathrm{0}}=\mathrm{0}(26)

求解方程(25)、(26)得到通解如下：

\overline{\Phi }\left(s,z\right)=s^{-\mathrm{1}}{A}_{\mathrm{1}}\left(s\right){\mathrm{e}}^{-sz}(27)

\overline{\Psi }\left(s,z\right)=s^{-\mathrm{1}}{A}_{\mathrm{2}}\left(s\right){\mathrm{e}}^{-sz}+s^{-\mathrm{1}}{A}_{\mathrm{3}}\left(s\right){\mathrm{e}}^{-\zeta z/l}(28)

其中，\zeta =\sqrt[]{\mathrm{1}+l^{\mathrm{2}}{s}^{\mathrm{2}}}，A_{\mathrm{1}}\mathrm{、}{A}_{\mathrm{2}}\mathrm{、}{A}_{\mathrm{3}}为边界条件确定的待定系数。

同理，将旋转矢量方程(3)和应力分量方程(10)、方程(11)、方程(15)以及位移分量(18)、(19)对r进行Hankel积分变换得到如下方程：

{〈u_r〉}_{\mathrm{1}}=l^{\mathrm{2}}s\frac{\mathrm{d}}{\mathrm{d}z}{〈\Psi 〉}_{\mathrm{0}}+\alpha sz\left(\mathrm{1}-l^{\mathrm{2}}\left(\frac{{\mathrm{d}}^{\mathrm{2}}}{\mathrm{d}{z}^{\mathrm{2}}}-s^{\mathrm{2}}\right)\right){〈\Psi 〉}_{\mathrm{0}}+\alpha s{〈\Phi 〉}_{\mathrm{0}}(29)

\begin{array}{l}{〈u_z〉}_{\mathrm{0}}=\left(\mathrm{1}-l^{\mathrm{2}}\frac{{\mathrm{d}}^{\mathrm{2}}}{\mathrm{d}{z}^{\mathrm{2}}}\right){〈\Psi 〉}_{\mathrm{0}}\\ -\alpha \frac{\mathrm{d}}{\mathrm{d}z}\left(z\left(\mathrm{1}-l^{\mathrm{2}}\left(\frac{{\mathrm{d}}^{\mathrm{2}}}{\mathrm{d}{z}^{\mathrm{2}}}-s^{\mathrm{2}}\right)\right){〈\Psi 〉}_{\mathrm{0}}+{〈\Phi 〉}_{\mathrm{0}}\right)\end{array}(30)

{〈{\omega }_{\theta }〉}_{\mathrm{1}}=\frac{\mathrm{1}}{\mathrm{2}}s{〈\Psi 〉}_{\mathrm{0}}(31)

\begin{array}{l}{〈{\sigma }_{zr}〉}_{\mathrm{1}}=\mu \left[\frac{\mathrm{d}{〈u_r〉}_{\mathrm{1}}}{\mathrm{d}z}-s{〈u_z〉}_{\mathrm{0}}\right.\\ \left.-l^{\mathrm{2}}\left(\frac{{\mathrm{d}}^{\mathrm{2}}}{\mathrm{d}{z}^{\mathrm{2}}}-s^{\mathrm{2}}\right)\left(\frac{\mathrm{d}{〈u_r〉}_{\mathrm{1}}}{\mathrm{d}z}+s{〈u_z〉}_{\mathrm{0}}\right)\right]\end{array}(32)

{〈{\sigma }_{zz}〉}_{\mathrm{0}}=\lambda s{〈u_r〉}_{\mathrm{1}}+\left(\lambda +\mathrm{2}\mu \right)\frac{\mathrm{d}}{\mathrm{d}z}{〈u_z〉}_{\mathrm{0}}(33)

{〈m_{z\theta }〉}_{\mathrm{1}}=\mathrm{4}\mu l^{\mathrm{2}}\frac{\mathrm{d}}{\mathrm{d}z}{〈{\omega }_{\theta }〉}_{\mathrm{1}}(34)

将通解式(27)和(28)代入方程(29)~(34)并用矩阵表示结果如下：

\left\{S\right\}=\left[T\left(z\right)\right]\left\{A\right\}(35)

\left\{S\right\}={\left[{〈u_r〉}_{\mathrm{1}},{〈u_z〉}_{\mathrm{0}},{〈{\omega }_{\theta }〉}_{\mathrm{1}},{〈{\sigma }_{zz}〉}_{\mathrm{0}},{〈{\sigma }_{zr}〉}_{\mathrm{1}},{〈m_{z\theta }〉}_{\mathrm{1}}\right]}^{\mathrm{T}}(36)

\left[T(z)\right]=\left[\begin{array}{ccc}\alpha {\mathrm{e}}^{-sz}& \left(\alpha z-s{l}^{\mathrm{2}}\right){\mathrm{e}}^{-sz}& -l\zeta {\mathrm{e}}^{-\zeta z/l}\\ \alpha {\mathrm{e}}^{-sz}& \frac{\left[\mathrm{1}-s^{\mathrm{2}}{l}^{\mathrm{2}}-\alpha \left(\mathrm{1}-sz\right)\right]{\mathrm{e}}^{-sz}}{s}& -s{l}^{\mathrm{2}}{\mathrm{e}}^{-\zeta z/l}\\ \mathrm{0}& {\mathrm{e}}^{-sz}/\mathrm{2}& {\mathrm{e}}^{-\zeta z/l}/\mathrm{2}\\ -\mathrm{2}\mu \alpha s{\mathrm{e}}^{-sz}& \left[\lambda \left(\mathrm{2}\alpha -\mathrm{1}\right)-\mathrm{2}\mu \left(\mathrm{1}-s^{\mathrm{2}}{l}^{\mathrm{2}}-\alpha \left(\mathrm{2}-sz\right)\right)\right]{\mathrm{e}}^{-sz}& \mathrm{2}\mu sl\zeta {\mathrm{e}}^{-\zeta z/l}\\ -\mathrm{2}\mu \alpha s{\mathrm{e}}^{-sz}& \mu \left[\mathrm{2}{s}^{\mathrm{2}}{l}^{\mathrm{2}}+\mathrm{2}\alpha \left(\mathrm{1}-sz\right)-\mathrm{1}\right]{\mathrm{e}}^{-sz}& \mathrm{2}\mu s^{\mathrm{2}}{l}^{\mathrm{2}}{\mathrm{e}}^{-\zeta z/l}\\ \mathrm{0}& {\mathrm{e}}^{-sz}/\mathrm{2}& -\mathrm{2}\mu l\zeta {\mathrm{e}}^{-\zeta z/l}\end{array}\right](37)

\left\{A\right\}={\left[A_{\mathrm{1}},A_{\mathrm{2}},A_{\mathrm{3}}\right]}^{\mathrm{{\rm T}}}(38)

给出边界条件如下：

{\sigma }_{zz}\left(r,\mathrm{0}\right)=-p\left(r\right), \mathrm{0}\le r\le a(39)

(40)

(41)

(42)

(43)

将边界条件式(39)~(43)代入式(35)求出待定系数

\begin{array}{l}{A}_{\mathrm{1}}=\frac{\mathrm{2}\left(\mathrm{1}-v\right)}{\mu sN\left(s\right)}\left[N\left(s\right)-\mathrm{2}\left(\mathrm{1}-v\right)\zeta \right]{〈p\left(s\right)〉}_{\mathrm{0}}\\ -\frac{\mathrm{4}\zeta {\left(\mathrm{1}-v\right)}^{\mathrm{2}}}{s\mu N\left(s\right)}{〈q\left(s\right)〉}_{\mathrm{1}}\end{array}(44)

A_{\mathrm{2}}=\frac{\mathrm{2}\left(\mathrm{1}-v\right)\zeta }{\mu N\left(s\right)}{〈p\left(s\right)〉}_{\mathrm{0}}+\frac{\mathrm{2}\left(\mathrm{1}-v\right)\zeta }{\mu N\left(s\right)}{〈q\left(s\right)〉}_{\mathrm{1}}(45)

A_{\mathrm{3}}=-\frac{\mathrm{2}\left(\mathrm{1}-v\right)ls}{\mu N\left(s\right)}{〈p\left(s\right)〉}_{\mathrm{0}}-\frac{\mathrm{2}\left(\mathrm{1}-v\right)ls}{\mu N\left(s\right)}{〈q\left(s\right)〉}_{\mathrm{1}}(46)

其中：

N\left(s\right)=\zeta +\mathrm{4}\left(\mathrm{1}-v\right)l^{\mathrm{2}}{s}^{\mathrm{2}}\left(\zeta -ls\right)(47)

{〈p\left(s\right)〉}_{\mathrm{0}}={\int }_{\mathrm{0}}^{\mathrm{\infty }}rp\left(r\right)J_{\mathrm{0}}\left(sr\right)\mathrm{d}r(48)

{〈q\left(s\right)〉}_{\mathrm{1}}={\int }_{\mathrm{0}}^{\mathrm{\infty }}rq\left(r\right)J_{\mathrm{1}}\left(sr\right)\mathrm{d}r(49)

将式(44)~(46)代入式(35)并取z=\mathrm{0}就可得到表面法向位移{〈u_z\left(s\right)〉}_{\mathrm{0}}和粘着区径向位移{〈u_r\left(s\right)〉}_{\mathrm{1}}：

{〈u_z\left(s\right)〉}_{\mathrm{0}}=\frac{\mathrm{2}\left(\mathrm{1}-\nu \right)\zeta -N\left(s\right)}{\mathrm{2}\mu sN\left(s\right)}{〈q\left(s\right)〉}_{\mathrm{1}}+\frac{\left(\mathrm{1}-\nu \right)\zeta }{\mu sN\left(s\right)}{〈p\left(s\right)〉}_{\mathrm{0}}(50)

{〈u_r\left(s\right)〉}_{\mathrm{1}}=\frac{N\left(s\right)-\mathrm{2}\left(\mathrm{1}-\nu \right)\zeta }{\mathrm{2}\mu sN\left(s\right)}{〈p\left(s\right)〉}_{\mathrm{0}}+\frac{\left(\nu -\mathrm{1}\right)}{\mu sN\left(s\right)}{〈q\left(s\right)〉}_{\mathrm{1}}(51)

将式(50)和式(51)通过式(24)进行Hankel逆变换，得到变换后的法向位移u_z\left(r\right)和粘着区径向位移u_r\left(r\right)方程为

u_z\left(r\right)={\int }_{\mathrm{0}}^{a}G\left(r,t\right)tp\left(t\right)dt+{\int }_{\mathrm{0}}^{a}F\left(r,t\right)tq\left(t\right)\mathrm{d}t(52)

u_r\left(r\right)={\int }_{\mathrm{0}}^{a}O\left(r,t\right)tp\left(t\right)dt+{\int }_{\mathrm{0}}^{a}Q\left(r,t\right)tq\left(t\right)\mathrm{d}t(53)

其中：

\left\{\begin{array}{c}G\left(r,t\right)={\int }_{\mathrm{0}}^{\mathrm{\infty }}{M}_{\mathrm{1}}\left(s\right)J_{\mathrm{0}}\left(sr\right)J_{\mathrm{0}}\left(st\right)\mathrm{d}s\\ F\left(r,t\right)={\int }_{\mathrm{0}}^{\mathrm{\infty }}{M}_{\mathrm{2}}\left(s\right)J_{\mathrm{0}}\left(sr\right)J_{\mathrm{1}}\left(st\right)\mathrm{d}s\end{array}\right.(54)

\left\{\begin{array}{c}O\left(r,t\right)={\int }_{\mathrm{0}}^{\mathrm{\infty }}{M}_{\mathrm{3}}\left(s\right)J_{\mathrm{1}}\left(sr\right)J_{\mathrm{0}}\left(st\right)\mathrm{d}s\\ Q\left(r,t\right)={\int }_{\mathrm{0}}^{\mathrm{\infty }}{M}_{\mathrm{4}}\left(s\right)J_{\mathrm{1}}\left(sr\right)J_{\mathrm{1}}\left(st\right)\mathrm{d}s\end{array}\right.(55)

\left\{\begin{array}{l}{M}_{\mathrm{1}}\left(s\right)=\frac{\mathrm{1}-v}{\mu }\frac{\zeta }{N\left(s\right)}\\ M_{\mathrm{2}}\left(s\right)=\frac{\mathrm{1}}{\mathrm{2}\mu }\frac{\mathrm{2}\left(\mathrm{1}-v\right)\zeta -N\left(s\right)}{N(s)}\end{array}\right.(56)

\left\{\begin{array}{l}{M}_{\mathrm{3}}\left(s\right)=\frac{\mathrm{1}}{\mathrm{2}\mu }\frac{N\left(s\right)-\mathrm{2}\left(\mathrm{1}-v\right)\zeta }{N\left(s\right)}\\ M_{\mathrm{4}}\left(s\right)=-\frac{\mathrm{1}-v}{\mu }\frac{\zeta }{N\left(s\right)}\end{array}\right.(57)

考虑M_{\mathrm{1}}\left(s\right)\mathrm{、}{M}_{\mathrm{2}}\left(s\right)\mathrm{、}{M}_{\mathrm{3}}\left(s\right)\mathrm{、}{M}_{\mathrm{4}}\left(s\right)函数的渐近性^[13]可知

\left\{\begin{array}{c}{\beta }_{\mathrm{1}}=\underset{s\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{M}_{\mathrm{1}}\left(s\right)=\frac{\mathrm{1}-v}{\mu }\frac{\mathrm{1}}{\mathrm{3}-\mathrm{2}v}\\ {\beta }_{\mathrm{2}}=\underset{s\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{M}_{\mathrm{2}}\left(s\right)=-\frac{\mathrm{1}}{\mathrm{2}\mu }\frac{\mathrm{1}}{\mathrm{3}-\mathrm{2}v}\end{array}\right.(58)

\left\{\begin{array}{l}{\beta }_{\mathrm{3}}=\underset{s\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{M}_{\mathrm{3}}\left(s\right)=\frac{\mathrm{1}}{\mathrm{2}\mu }\frac{\mathrm{1}}{\mathrm{3}-\mathrm{2}v}\\ {\beta }_{\mathrm{4}}=\underset{s\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{M}_{\mathrm{4}}\left(s\right)=-\frac{\mathrm{1}-v}{\mu }\frac{\mathrm{1}}{\mathrm{3}-\mathrm{2}v}\end{array}\right.(59)

通过以下方法将方程改写，令

\left\{\begin{array}{c}G\left(r,t\right)={\beta }_{\mathrm{1}}{G}_{\mathrm{1}}+G_{\mathrm{2}}\\ F\left(r,t\right)={\beta }_{\mathrm{2}}{F}_{\mathrm{1}}+F_{\mathrm{2}}\end{array}\right.(60)

\left\{\begin{array}{c}O\left(r,t\right)={\beta }_{\mathrm{3}}{O}_{\mathrm{1}}+O_{\mathrm{2}}\\ Q\left(r,t\right)={\beta }_{\mathrm{4}}{Q}_{\mathrm{1}}+Q_{\mathrm{2}}\end{array}\right.(61)

(62)

(63)

(64)

(65)

\left\{\begin{array}{c}{G}_{\mathrm{2}}={\int }_{\mathrm{0}}^{\mathrm{\infty }} \left[M_{\mathrm{1}}\left(s\right)-{\beta }_{\mathrm{1}}\right]J_{\mathrm{0}}\left(sr\right)J_{\mathrm{0}}\left(st\right)\mathrm{d}s\\ F_{\mathrm{2}}={\int }_{\mathrm{0}}^{\mathrm{\infty }} \left[M_{\mathrm{2}}\left(s\right)-{\beta }_{\mathrm{2}}\right]J_{\mathrm{0}}\left(sr\right)J_{\mathrm{1}}\left(st\right)\mathrm{d}s\end{array}\right.(66)

\left\{\begin{array}{c}{O}_{\mathrm{2}}={\int }_{\mathrm{0}}^{\mathrm{\infty }}\left[M_{\mathrm{3}}\left(s\right)-{\beta }_{\mathrm{3}}\right]J_{\mathrm{1}}\left(sr\right)J_{\mathrm{0}}\left(st\right)\mathrm{d}s\\ Q_{\mathrm{2}}={\int }_{\mathrm{0}}^{\mathrm{\infty }}\left[M_{\mathrm{4}}\left(s\right)-{\beta }_{\mathrm{4}}\right]J_{\mathrm{1}}\left(sr\right)J_{\mathrm{1}}\left(st\right)\mathrm{d}s\end{array}\right.(67)

其中，K[⋅]和E[⋅]分别表示第一类和第二类完全椭圆函数^[13]。

最终得到改写后的位移方程为

\begin{array}{l}{u}_z\left(r\right)={\beta }_{\mathrm{1}}{\int }_{\mathrm{0}}^a{G}_{\mathrm{1}}\left(r,t\right)p\left(t\right)t\mathrm{d}t+{\int }_{\mathrm{0}}^a{G}_{\mathrm{2}}\left(r,t\right)p\left(t\right)t\mathrm{d}t\\ +{\beta }_{\mathrm{2}}{\int }_{\mathrm{0}}^a{F}_{\mathrm{1}}\left(r,t\right)q\left(t\right)t\mathrm{d}t+{\int }_{\mathrm{0}}^a{F}_{\mathrm{2}}\left(r,t\right)q\left(t\right)t\mathrm{d}t\end{array}(68)

\begin{array}{l}{u}_r\left(r\right)={\beta }_{\mathrm{3}}{\int }_{\mathrm{0}}^a{O}_{\mathrm{1}}\left(r,t\right)p\left(t\right)t\mathrm{d}t+{\int }_{\mathrm{0}}^a{O}_{\mathrm{2}}\left(r,t\right)p\left(t\right)t\mathrm{d}t\\ +{\beta }_{\mathrm{4}}{\int }_{\mathrm{0}}^a{Q}_{\mathrm{1}}\left(r,t\right)q\left(t\right)t\mathrm{d}t+{\int }_{\mathrm{0}}^a{Q}_{\mathrm{2}}\left(r,t\right)q\left(t\right)t\mathrm{d}t\end{array}(69)

当式(68)不考虑切向应力时，方程将退化为弹性半空间无摩擦轴对称接触问题的法向位移方程^[9]，即

u_z\left(r\right)={\beta }_{\mathrm{1}}{\int }_{\mathrm{0}}^a{G}_{\mathrm{1}}\left(r,t\right)p\left(t\right)t\mathrm{d}t+{\int }_{\mathrm{0}}^a{G}_{\mathrm{2}}\left(r,t\right)p\left(t\right)t\mathrm{d}t(70)

2 奇异积分方程的推导

将式(68)和式(69)分别对r求导，并且对p(t)进行偶延拓，对q(t)进行奇延拓，得到如下方程：

\begin{array}{l}{m}_{\mathrm{1}}\left(r\right)=\frac{\mathrm{1}}{\mathrm{2}}{\int }_{-a}^{a}p\left(t\right)\left|t\right|L_{\mathrm{1}}\left(r,t\right)\mathrm{d}t\\ +\frac{\mathrm{1}}{\mathrm{2}}{\int }_{-a}^{a}q\left(t\right)\left|t\right|L_{\mathrm{2}}\left(r,t\right)\mathrm{d}t+\frac{{\beta }_{\mathrm{1}}}{\mathrm{\pi }}{\int }_{-a}^a\frac{p\left(t\right)}{t-r}\mathrm{d}t\\ +\frac{{\beta }_{\mathrm{1}}}{\mathrm{\pi }}{\int }_{-a}^{a}p\left(t\right)H_{\mathrm{1}}\left(r,t\right)\mathrm{d}t-{\beta }_{\mathrm{2}}q(r)\end{array}(71)

\begin{array}{l}{m}_{\mathrm{2}}\left(r\right)=\frac{\mathrm{1}}{\mathrm{2}}{\int }_{-a}^{a}p\left(t\right)\left|t\right|L_{\mathrm{3}}\left(r,t\right)\mathrm{d}t\\ +\frac{\mathrm{1}}{\mathrm{2}}{\int }_{-a}^{a}q\left(t\right)\left|t\right|L_{\mathrm{4}}\left(r,t\right)\mathrm{d}t+\frac{{\beta }_{\mathrm{4}}}{\mathrm{\pi }}{\int }_{-a}^a\frac{q\left(t\right)}{t-r}\mathrm{d}t\\ +\frac{{\beta }_{\mathrm{4}}}{\mathrm{\pi }}{\int }_{-a}^{a}q\left(t\right)H_{\mathrm{2}}\left(r,t\right)\mathrm{d}t+{\beta }_{\mathrm{3}}p\left(r\right)\end{array}(72)

其中：

m_{\mathrm{1}}\left(r\right)=\frac{\partial u_{z\mathrm{0}}}{\partial r},m_{\mathrm{2}}\left(r\right)=\frac{\mathrm{1}}{r}\frac{\partial r{u}_{r\mathrm{0}}}{\partial r}(73)

H_{\mathrm{1}}\left(r,t\right)=\frac{h_{\mathrm{1}}\left(r,t\right)-\mathrm{1}}{t-r},H_{\mathrm{2}}\left(r,t\right)=\frac{h_{\mathrm{2}}\left(r,t\right)-\mathrm{1}}{t-r}(74)

(75)

(76)

L_{\mathrm{1}}\left(r,t\right)=-{\int }_{\mathrm{0}}^{\mathrm{\infty }}\left[M_{\mathrm{1}}\left(s\right)-{\beta }_{\mathrm{1}}\right]s{J}_{\mathrm{1}}\left(sr\right)J_{\mathrm{0}}\left(st\right)\mathrm{d}s(77)

L_{\mathrm{2}}\left(r,t\right)=-{\int }_{\mathrm{0}}^{\mathrm{\infty }}\left[M_{\mathrm{2}}\left(s\right)-{\beta }_{\mathrm{2}}\right]s{J}_{\mathrm{1}}\left(sr\right)J_{\mathrm{1}}\left(st\right)\mathrm{d}s(78)

L_{\mathrm{3}}\left(r,t\right)={\int }_{\mathrm{0}}^{\mathrm{\infty }}\left[M_{\mathrm{3}}\left(s\right)-{\beta }_{\mathrm{3}}\right]s{J}_{\mathrm{0}}\left(sr\right)J_{\mathrm{0}}\left(st\right)\mathrm{d}s(79)

L_{\mathrm{4}}\left(r,t\right)={\int }_{\mathrm{0}}^{\mathrm{\infty }}\left[M_{\mathrm{4}}\left(s\right)-{\beta }_{\mathrm{4}}\right]s{J}_{\mathrm{0}}\left(sr\right)J_{\mathrm{1}}\left(st\right)\mathrm{d}s(80)

对于球形压头，其法相位移满足下式^[11]

u_z\left(r\right)=\delta -\frac{r^{\mathrm{2}}}{\mathrm{2}R}(81)

其中\delta为最大压痕深度。

由式(81)可得到

m_{\mathrm{1}}\left(r\right)=\frac{\partial u_z\left(r\right)}{\partial r}=-\frac{r}{R}(82)

对于均匀材料的法相加载问题Spence^[14]提出了自相似理论，粘着区径向位移u_r(r)可表示为

u_r\left(r\right)=B{r}^{\mathrm{2}}(83)

对式(83)中r求导并延拓可得到

m_{\mathrm{2}}\left(r\right)=\frac{\mathrm{1}}{r}\frac{\partial r{u}_{r\mathrm{0}}}{\partial r}=\mathrm{3}B\left|r\right|(84)

方程满足平衡条件如下：

P=\mathrm{\pi }{\int }_{-a}^{a}p\left(t\right)\left|t\right|\mathrm{d}t(85)

\mathrm{0}={\int }_{-a}^{a}q\left(t\right)\left|t\right|\mathrm{d}t(86)

为方便利用数值方法求解方程，采用如下方法对方程无量纲

x=\frac{t}{a},y=\frac{r}{a},L=\frac{l}{a}(87)

则方程(71)和方程(72)变换为

\begin{array}{l}{m}_{\mathrm{1}}\left(y\right)=\frac{a}{\mathrm{2}}{\int }_{-\mathrm{1}}^{\mathrm{1}}p\left(x\right)\left|ax\right|L_{\mathrm{1}}\left(y,x\right)\mathrm{d}x\\ +\frac{a}{\mathrm{2}}{\int }_{-\mathrm{1}}^{\mathrm{1}}q\left(x\right)\left|ax\right|L_{\mathrm{2}}\left(y,x\right)\mathrm{d}x+\frac{{\beta }_{\mathrm{1}}}{\mathrm{\pi }}{\int }_{-\mathrm{1}}^{\mathrm{1}}\frac{p\left(x\right)}{x-y}\mathrm{d}x\\ +\frac{{\beta }_{\mathrm{1}}}{\mathrm{\pi }}{\int }_{-\mathrm{1}}^{\mathrm{1}}p\left(x\right)H_{\mathrm{1}}\left(y,x\right)\mathrm{d}x-{\beta }_{\mathrm{2}}q\left(y\right)(\mathrm{88})\end{array}

\begin{array}{l}{m}_{\mathrm{2}}\left(y\right)=\frac{a}{\mathrm{2}}{\int }_{-\mathrm{1}}^{\mathrm{1}}p\left(x\right)\left|ax\right|L_{\mathrm{3}}\left(y,x\right)\mathrm{d}x\\ +\frac{a}{\mathrm{2}}{\int }_{-\mathrm{1}}^{\mathrm{1}}q\left(x\right)\left|ax\right|L_{\mathrm{4}}\left(y,x\right)\mathrm{d}x+\frac{{\beta }_{\mathrm{4}}}{\mathrm{\pi }}{\int }_{-\mathrm{1}}^{\mathrm{1}}\frac{q\left(x\right)}{x-y}\mathrm{d}x\\ +\frac{{\beta }_{\mathrm{4}}}{\mathrm{\pi }}{\int }_{-\mathrm{1}}^{\mathrm{1}}q\left(x\right)H_{\mathrm{2}}\left(y,x\right)\mathrm{d}x+{\beta }_{\mathrm{3}}p\left(y\right)(\mathrm{89})\end{array}

平衡条件改写为

P=\mathrm{\pi }{a}^{\mathrm{2}}{\int }_{-\mathrm{1}}^{\mathrm{1}}p\left(x\right)\left|x\right|\mathrm{d}x(90)

\mathrm{0}={\int }_{-\mathrm{1}}^{\mathrm{1}}q\left(x\right)\left|x\right|\mathrm{d}x(91)

式中H_{\mathrm{1}}\left(y,x\right)和H_{\mathrm{2}}\left(y,x\right)满足

\left\{\begin{array}{l}\underset{x\to y}{\mathrm{l}\mathrm{i}\mathrm{m}}{H}_{\mathrm{1}}\left(y,x\right)=\frac{\mathrm{1}}{\mathrm{2}y}\mathrm{l}\mathrm{o}\mathrm{g}\left|x-y\right|\\ \underset{x\to y}{\mathrm{l}\mathrm{i}\mathrm{m}}{H}_{\mathrm{2}}\left(y,x\right)=-\frac{\mathrm{1}}{\mathrm{2}y}\mathrm{l}\mathrm{o}\mathrm{g}\left|x-y\right|\end{array}\right.(92)

3 完全粘着接触问题求解

用Goodman^[15]近似方法求解方程，即忽略切向应力对法向应力的影响，则方程(88)可被解耦，方程简化为

\begin{array}{l}{m}_{\mathrm{1}}\left(y\right)=\frac{a}{\mathrm{2}}{\int }_{-\mathrm{1}}^{\mathrm{1}}p\left(x\right)\left|ax\right|L_{\mathrm{1}}\left(y,x\right)\mathrm{d}x+\frac{{\beta }_{\mathrm{1}}}{\mathrm{\pi }}{\int }_{-\mathrm{1}}^{\mathrm{1}}\frac{p\left(x\right)}{x-y}\mathrm{d}x\\ +\frac{{\beta }_{\mathrm{1}}}{\mathrm{\pi }}{\int }_{-\mathrm{1}}^{\mathrm{1}}p\left(x\right)H_{\mathrm{1}}\left(y,x\right)\mathrm{d}x(\mathrm{93})\end{array}

在球形压头作用时，假设应力分布函数^[16]

p\left(x\right)=f\left(x\right)\sqrt[]{\mathrm{1}-x^{\mathrm{2}}}(94)

根据Erdogan-Gupta^[17]提出的方法，方程(93)可离散为

{\displaystyle \sum _{i=\mathrm{1}}^N}\frac{\mathrm{\pi }\left(\mathrm{1}-x_i^{\mathrm{2}}\right)f\left(x_i\right)}{N+\mathrm{1}}\left[\frac{\mathrm{1}}{\mathrm{\pi }}\frac{\mathrm{1}}{x_i-y_j}+\frac{\mathrm{1}}{\mathrm{\pi }}\frac{\mathrm{l}\mathrm{o}\mathrm{g}\left|x_i-y_j\right|}{\mathrm{2}{y}_j}+\frac{\mathrm{1}}{{\beta }_{\mathrm{1}}}{I}_{\mathrm{1}}\left(y_j,x_i\right)\right]+\frac{\varphi \left[f\left(x_{j-\mathrm{1}}\right),f\left(x_j\right)\right]}{\mathrm{2}{y}_j}=-\frac{\mathrm{1}}{{\beta }_{\mathrm{1}}}\frac{a{y}_j}{R}(95)

其中：

x_i=\mathrm{c}\mathrm{o}\mathrm{s}\left[i\mathrm{\pi }/\left(N+\mathrm{1}\right)\right],i=\left(\mathrm{1},\mathrm{2},\dots ,N\right)(96)

y_j=\mathrm{c}\mathrm{o}\mathrm{s}\left[\mathrm{\pi }\left(\mathrm{2}j-\mathrm{1}\right)/\mathrm{2}\left(N+\mathrm{1}\right)\right],\left(j=(\mathrm{1},\mathrm{2},\dots ,N,N+\mathrm{1}\right)(97)

I_{\mathrm{1}}=\frac{a^{\mathrm{2}}}{\mathrm{2}}\left|x\right|L_{\mathrm{1}}\left(y,x\right)+\frac{{\beta }_{\mathrm{1}}}{\mathrm{\pi }}\left[\frac{h_{\mathrm{1}}\left(y,x\right)-\mathrm{1}}{x-y}-\frac{\mathrm{1}}{\mathrm{2}y}\mathrm{l}\mathrm{o}\mathrm{g}\left|x-y\right|\right](98)

\begin{array}{c}\varphi \left[f\left(x_{j-\mathrm{1}}\right),f\left(x_j\right)\right]=\frac{f\left(x_{j-\mathrm{1}}\right)}{\mathrm{\pi }}\left\{\left(\mathrm{1}-y_j^{\mathrm{2}}\right)\left[\eta \mathrm{l}\mathrm{o}\mathrm{g}\left|\eta {\left(\mathrm{1}-y_j^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}}\right|-\eta \right]-\eta \left(\mathrm{1}-x_j^{\mathrm{2}}\right)\mathrm{l}\mathrm{o}\mathrm{g}\left|x_{j-\mathrm{1}}-y_j\right|\right\}\\ +\frac{f\left(x_j\right)}{\mathrm{\pi }}\left\{\left(\mathrm{1}-y_j^{\mathrm{2}}\right)\left[\eta \mathrm{l}\mathrm{o}\mathrm{g}\left|\eta {\left(\mathrm{1}-y_j^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}}\right|-\eta \right]-\eta \left(\mathrm{1}-x_j^{\mathrm{2}}\right)\mathrm{l}\mathrm{o}\mathrm{g}\left|x_j-y_j\right|\right\}\end{array}(99)

式中：N为f\left(x_i\right)在(-1,1)的离散点的总数，\eta为y_j的领域范围。

对于方程(95)，共有N+\mathrm{1}个未知数，即f\left(x_{\mathrm{1}}\right)\mathrm{、}f\left(x_{\mathrm{2}}\right)\mathrm{、}f\left(x_{\mathrm{2}}\right),\cdots ,f\left(x_N\right)和a，在实际计算中，a给定，这里N+\mathrm{1}个方程确定N个未知数，按照文献[17]给出的方法，N取偶数，并舍弃\left(j=\left(N+\mathrm{1}\right)/\mathrm{2}+\mathrm{1}\right)对应的方程，那么N个方程解N个未知数，即f\left(x_{\mathrm{1}}\right)\mathrm{、}f\left(x_{\mathrm{2}}\right)\mathrm{、}f\left(x_{\mathrm{2}}\right),\cdots ,f\left(x_N\right)。利用插值函数就可以得到每一点的应力分布p(x)。

以上是利用Goodman近似得到的近似解，可以通过迭代的方法，得到方程(88)和方程(89)的耦合解，假设

p\left(x\right)=f_{\mathrm{1}}\left(x\right)\sqrt[]{\mathrm{1}-x^{\mathrm{2}}},q\left(x\right)=f_{\mathrm{2}}\left(x\right)\sqrt[]{\mathrm{1}-x^{\mathrm{2}}}(100)

同理，将方程(87)和方程(88)离散为

{\displaystyle \sum _{i=\mathrm{1}}^N}\frac{\mathrm{\pi }\left(\mathrm{1}-x_i^{\mathrm{2}}\right)f_{\mathrm{1}}\left(x_i\right)}{N+\mathrm{1}}\left[\frac{\mathrm{1}}{\mathrm{\pi }}\frac{\mathrm{1}}{x_i-y_j}+\frac{\mathrm{1}}{\mathrm{\pi }}\frac{\mathrm{l}\mathrm{o}\mathrm{g}\left|x_i-y_j\right|}{\mathrm{2}{y}_j}+\frac{\mathrm{1}}{{\beta }_{\mathrm{1}}}{I}_{\mathrm{1}}\left(y_j,x_i\right)\right]+\frac{\varphi \left[f_{\mathrm{1}}(x_{j-\mathrm{1}}),f_{\mathrm{1}}(x_j)\right]}{\mathrm{2}{y}_j}=\frac{\mathrm{1}}{{\beta }_{\mathrm{1}}}{E}_{\mathrm{1}}\left(y_j\right)(101)

{\displaystyle \sum _{i=\mathrm{1}}^N}\frac{\mathrm{\pi }\left(\mathrm{1}-x_i^{\mathrm{2}}\right)f_{\mathrm{2}}\left(x_i\right)}{N+\mathrm{1}}\left[\frac{\mathrm{1}}{\mathrm{\pi }}\frac{\mathrm{1}}{x_i-y_j}-\frac{\mathrm{1}}{\mathrm{\pi }}\frac{\mathrm{l}\mathrm{o}\mathrm{g}\left|x_i-y_j\right|}{\mathrm{2}{y}_j}+\frac{\mathrm{1}}{{\beta }_{\mathrm{4}}}{I}_{\mathrm{4}}\left(y_j,x_i\right)\right]+\frac{\varphi \left[f_{\mathrm{2}}\left(x_{j-\mathrm{1}}\right),f_{\mathrm{2}}\left(x_j\right)\right]}{\mathrm{2}{y}_j}=\frac{\mathrm{1}}{{\beta }_{\mathrm{4}}}{E}_{\mathrm{2}}\left(y_j\right)(102)

其中

\left\{\begin{array}{l}{E}_{\mathrm{1}}\left(y\right)=-\frac{ay}{R}+{\beta }_{\mathrm{2}}q\left(y\right)-{\int }_{-\mathrm{1}}^{\mathrm{1}}q\left(x\right)I_{\mathrm{2}}\mathrm{d}x\\ E_{\mathrm{2}}\left(y\right)=\mathrm{3}Ba\left|y\right|-{\beta }_{\mathrm{3}}p\left(y\right)-{\int }_{-\mathrm{1}}^{\mathrm{1}}p\left(x\right)I_{\mathrm{3}}\mathrm{d}x\end{array}\right.(103)

x_i=\mathrm{c}\mathrm{o}\mathrm{s}\left[i\mathrm{\pi }/\left(N+\mathrm{1}\right)\right],i=\left(\mathrm{1},\mathrm{2},\dots ,N\right)(104)

y_j=\mathrm{c}\mathrm{o}\mathrm{s}\left[\mathrm{\pi }\left(\mathrm{2}j-\mathrm{1}\right)/\mathrm{2}\left(N+\mathrm{1}\right)\right],j=\left(\mathrm{1},\mathrm{2},\dots ,N,N+\mathrm{1}\right)(105)

\left\{\begin{array}{l}{I}_{\mathrm{1}}=\frac{a^{\mathrm{2}}}{\mathrm{2}{R}^{\mathrm{2}}}\left|x\right|L_{\mathrm{1}}\left(y,x\right)+\frac{{\beta }_{\mathrm{1}}}{\mathrm{\pi }}\left[\frac{h_{\mathrm{1}}\left(y,x\right)-\mathrm{1}}{x-y}-\frac{\mathrm{1}}{\mathrm{2}y}\mathrm{l}\mathrm{o}\mathrm{g}\left|x-y\right|\right]\\ I_{\mathrm{2}}=\frac{a^{\mathrm{2}}}{\mathrm{2}{R}^{\mathrm{2}}}\left|x\right|L_{\mathrm{2}}\left(y,x\right)\\ I_{\mathrm{3}}=\frac{a^{\mathrm{2}}}{\mathrm{2}{R}^{\mathrm{2}}}\left|x\right|L_{\mathrm{3}}\left(y,x\right)\\ I_{\mathrm{4}}=\frac{a^{\mathrm{2}}}{\mathrm{2}{R}^{\mathrm{2}}}\left|x\right|L_{\mathrm{4}}\left(y,x\right)+\frac{{\beta }_{\mathrm{4}}}{\mathrm{\pi }}\left[\frac{h_{\mathrm{2}}\left(y,x\right)-\mathrm{1}}{x-y}-\frac{\mathrm{1}}{\mathrm{2}y}\mathrm{l}\mathrm{o}\mathrm{g}\left|x-y\right|\right]\end{array}\right.(106)

\begin{array}{c}\varphi \left[f_k\left(x_{j-\mathrm{1}}\right),f_k\left(x_j\right)\right]=\frac{f_k\left(x_{j-\mathrm{1}}\right)}{\mathrm{\pi }}\left\{\left(\mathrm{1}-y_j^{\mathrm{2}}\right)\left[\eta \mathrm{l}\mathrm{o}\mathrm{g}\left|\eta {\left(\mathrm{1}-y_j^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}}\right|-\eta \right]-\eta \left(\mathrm{1}-x_j^{\mathrm{2}}\right)\mathrm{l}\mathrm{o}\mathrm{g}\left|x_{j-\mathrm{1}}-y_j\right|\right\}\\ +\frac{f_k\left(x_j\right)}{\mathrm{\pi }}\left\{\left(\mathrm{1}-y_j^{\mathrm{2}}\right)\left[\eta \mathrm{l}\mathrm{o}\mathrm{g}\left|\eta {\left(\mathrm{1}-y_j^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}}\right|-\eta \right]-\eta \left(\mathrm{1}-x_j^{\mathrm{2}}\right)\mathrm{l}\mathrm{o}\mathrm{g}\left|x_j-y_j\right|\right\}\left(k=\mathrm{1},\mathrm{2}\right)\end{array}(107)

将方程(101)利用Goodman近似方法求解出p_{\mathrm{0}}\left(x\right)，将p_{\mathrm{0}}\left(x\right)代入方程(102)，对于方程(102)这里有N+\mathrm{1}个方程，有N+\mathrm{1}个未知数 f_{\mathrm{2}}\left(x_{\mathrm{1}}\right)\mathrm{、}{f}_{\mathrm{2}}\left(x_{\mathrm{2}}\right)\mathrm{、}f_{\mathrm{2}}\left(x_{\mathrm{2}}\right),\cdots ,f_{\mathrm{2}}\left(x_N\right)和B。通过迭代，再利用插值函数方法，可以得到耦合解p\left(x\right)和q\left(x\right)，其具体迭代求解过程如下：

1) 给定接触区半径a_{\mathrm{0}}和材料特征长度参数l，利用Goodman近似方法求解方程(101)得到p_{\mathrm{0}}\left(x\right)。

2) 将p_{\mathrm{0}}\left(x\right)代入方程(102)，求解得到q_{\mathrm{1}}\left(x\right)和B。

3) 再将q_{\mathrm{1}}\left(x\right)代入方程(101)，求解得到p_{\mathrm{1}}\left(x\right)。

4) 循环迭代得到p_i\left(x\right)和q_i\left(x\right)，当q_i\left(x\right)满足时停止迭代过程，此时的p_i\left(x\right)和q_i\left(x\right)即为所求耦合解。

5) 计算外载荷P和最大压痕\delta。

4 数值结果与讨论

本文采用迭代方法得到了尺度依赖的偶应力弹性半空间在球形压头下的完全粘着接触问题的应力分布以及压力-压痕关系曲线。在数值计算过程中，泊松比取1/3，离散点个数N取30个，迭代误差取10^-3。为验证本文模型和计算的可靠性，忽略摩擦切向应力的影响，将模型退化为偶应力弹性半空间轴对称无摩擦接触问题。根据文献[9]给出的参数，计算过程中取无量纲接触区半径为0.01，无量纲材料特征长度参数分别取0.1、1.0。如图2所示，给出了本文的计算结果与文献[9]结果的比较，其中实线表示文献结果，离散点表示本文的计算结果。从图中可以看出，本结果与文献结果吻合较好，这表明本文的计算方法可靠，满足精度要求。

图3给出了无量纲接触半径为0.01时，不同材料特征长度参数对法向接触应力分布（图3(a)）和切向接触应力分布（图3(b)）的影响。

从图3(a)中可以明显看出，当固定接触区半径时，随着无量纲材料特征长度参数从0.1增加到1.0，在接触区中心法向接触应力的峰值迅速增加，且材料特征长度参数较小时，对应的法向接触应力趋于平缓。该结果表明，随着材料特征长度的增加，材料的刚度增加，为保持相同的接触半径，需要施加更大的载荷，进而增大了接触应力。

从图3(b)中可以看出，在接触区中心切向应力为0，其最大值在靠近接触区边缘的位置，随着材料特征长度的增加，切向应力的最大值逐渐增加，其最大值位置向接触区中心靠拢。材料特征长度参数较小时，材料特征长度参数的变化对切向应力的影响较为明显，当材料特征长度增加到一定程度时，继续增加材料特征长度对切向应力分布影响效果减小。

以上结果表明，材料的材料特征长度参数会影响材料的刚度，当材料特征长度参数较小时，对材料力学性能影响较大。因此，可通过控制材料的材料特征长度参数来调整材料表面的接触应力分布，为材料的优化设计和缓解应力集中提供理论依据。

图4中给出了不同材料特征长度下外载荷P与接触区半径a(图4(a))和最大压痕\delta(图4(b))的关系曲线。

从图4(a)中可以看出，材料特征长度一定时，随着接触半径的增加，所需外载荷逐渐增加，在固定接触区半径时，随着材料特征长度的增大，所需外载荷逐渐增大。在外载荷相同时，随着材料特征长度参数的增加，形成的接触半径减小。该结果表明，材料特征长度的增加提高了材料抵抗变形的能力，进而需要更大的外载荷才能达到相同的接触半径。在施加的外载荷相同时，特征材料长度参数越大，产生的形变越小。

从图4(b)中可以看出，在外载荷相同时，随着材料特征长度参数的增加，最大压痕深度减小。最大压痕深度一定时，随着材料特征长度参数的增加，所需的外载荷逐渐增加。该结果也表明，随着材料特征长度参数的增加，材料抵抗变形的能力增强。另外还可以看出，材料特征长度参数较小时，材料特征长度参数的变化对压力-压痕关系曲线的影响越明显，材料特征长度参数较大时，材料特征长度参数的改变对压力-压痕关系影响效果减弱。结果表明，材料特征长度参数越小时，材料表现出的尺度效应越明显。

以上结果表明，材料特征长度参数的增加，使材料的刚度增加，材料特征长度参数越小，尺度效应越明显，可通过改变材料特征长度参数来改变材料的力学性能，进而调节接触表面外载荷与接触半径和压痕的关系，为抑制压痕的产生和缓解应力集中提供帮助。

5 结论

本文求解了偶应力弹性半空间在球形压头作用下的完全粘着接触问题。利用轴对称完全粘着接触问题的基本解，将偶应力弹性半空间在球形压头作用下的完全粘着接触模型转化成奇异积分方程组并进行求解。利用迭代方法，得到了尺度依赖的弹性半空间在球形压头作用下的应力分布，分析了材料特征长度参数对其接触应力分布、外载荷与接触半径、压力-压痕关系的影响。研究结果表明：

1) 利用本文建立的模型可有效地求解偶应力弹性半空间的完全粘着接触问题。

2) 在固定接触半径时，随着材料特征长度参数的增加，其法向接触应力和切向接触应力迅速增加。材料特征长度参数较小时，对应的表面接触应力分布更平缓。在外载荷确定时，随着特征材料长度参数的增加，最大压痕深度减小。结果表明，材料特征长度参数的增加可提高材料的刚度，使材料抵抗变形的能力增强。

3) 可通过控制材料的材料特征长度参数的变化来改善材料的力学性能，为材料的优化设计和抑制接触损伤提供一个新的方向。

参考文献

[1] HERTZ H. On the contact of elastic solids[J]. Journal for Die Reine undAngewandt Mathematik, 1882, 92: 156-171.

[2] MINDLIN R D. Influence of couple-stresses on stress concentrations[J]. Experimental Mechanics, 1963, 3(1): 1-7.

[3] 王广春, 郑伟, 姜华, 等. 纯铜微镦粗过程尺寸效应的试验研究[J]. 机械工程学报, 2012, 48(14): 32-37.

[4] ZISIS T, GOURGIOTIS P A, BAXEVANAKIS K P, et al. Some basic contact problems in couple stress elasticity [J]. International Journal of Solids and Structures, 2014, 51(11/12): 2084-2095.

[5] SONG H X, KE L L, WANG Y S. Sliding frictional contact analysis of an elastic solid with couple stresses[J]. International Journal of Mechanical Sciences, 2017, 133: 804-816.

[6] LI P X, LIU T J. The size effect in adhesive contact on a gradient nanostructured coating[J]. Nanotechnology, 2021, 32(23): 235704.

[7] LI P X, LIU T J. The two-dimensional adhesive contact problem in the theory of couple stress elasticity[J]. Journal of Adhesion Science and Technology, 2020, 34(10): 1062-1082.

[8] LI P X, LIU T J. The adhesive contact problem between a graded coated half-space and a cylindrical indenter by using a Maugis model[J]. Journal of Adhesion Science and Technology, 2018, 32(22): 2494-2508.

[9] GOURGIOTIS P A, ZISIS T, GIANNAKOPOULOS A E, et al. The hertz contact problem in couple-stress elasticity[J].International Journal of Solids and Structures, 2019, 168: 228-237.

[10] LI P X, LIU T J. The axisymmetric contact in couple-stress elasticity taking into account adhesion[J]. Journal of Adhesion Science and Technology, 2023, 37(1): 50-71.

[11] 刘铁军. 功能梯度材料轴对称接触力学及微动分析[D]. 北京: 北京交通大学, 2008.

[12] MINDLIN R D, TIERSTEN H F. Effects of couple-stresses in linear elasticity[J]. Archive for Rational Mechanics and Analysis, 1962, 11: 415-448.

[13] ANDREWS G E, ASKEY R, ROY R. Special functions[M]. Cambridge: University Press, 1999.

[14] SPENCE D A. The hertz contact problem with finite friction[J]. Journal of Elasticity, 1975, 5(3): 297-319.

[15] GOODMAN L E. Contact stress analysis of normally loaded rough spheres[J]. Journal of Applied Mechanics, 1962, 29(3): 515-522.

[16] CIVELEK M B. The axisymmetric contact problem for an elastic layer on a frictionlesshalf-space[D]. Bethlehem: Lehigh University, 1972.

[17] ERDOGAN F, GUPTA G D. On the numerical solution of singular integral equations[J]. Quarterly of Applied Mathematics, 1973, 29: 525-534.

基金资助

国家自然科学基金项目(12062019)