基于Eshelby-Stroh公式各向异性弹性体接触问题研究

颜灯灯，李成刚，*，申景金，王艳，王春明，宋伟山

(1.南京航空航天大学 机电学院，南京210016; 2.南京邮电大学 自动化学院，南京210023)

压头与弹性体之间平面问题的研究在工程力学中有着广泛的应用，通过建立载荷与压痕深度或变形曲线等参数之间的数学模型，再依据压痕试验便可以测量材料的力学性能。随着对压头和弹性体接触问题的研究不断深入，压痕试验被广泛应用于软骨组织的弹性模量测量［1］、生物组织的性能评估［2］和力学性能测量［3］、微创手术用的触诊探针［4］、弱锤冲击初始损伤试验［5］等方面。近年来，使用原子力显微法的压痕试验技术成为量化软质材料和生物材料力学性能的重要技术［6-7］。

在许多应用中，压头的弹性模量要比被压弹性体的大很多，因此压头一般都被视为刚体。被压的弹性体视为固定在一个刚体的基座上，即弹性体的其中一面完全固定。当压痕深度相对于弹性体的厚度、宽度以及压头某些参数来说比较小时，弹性体的变化可以视为是线性变化的。

平面压痕问题是混合边界值问题，多年来有着众多的研究。Muskhelishvili和England等使用解析延拓的方法求解典型结构压头的接触问题［8-9］，Gladwell和 England 使用特定的正交多项式展开的方法解决混合边界值接触问题［10］，Okumura等利用有理映射函数和复应力函数分析平面压痕问题［11］。上述研究面向各向同性的弹性材料，若将对象拓展到各向异性材料，求解方法变得复杂。结合Muskhelishvili［8］的解析延拓方法，使用Eshelby-Stroh公式可以求解任意轮廓的平面压头与各向异性弹性半空间体的接触问题［12-13］。Hwu和Fan基于Eshelby-Stroh公式和解析延拓，并使用共轭梯度法研究压头与各向异性弹性半空间体的接触问题［14-15］，可以考虑压头与弹性体之间的摩擦。Batra和Jiang基于Eshelby-Stroh公式设定级数函数表示位移函数和应力函数，研究了圆柱压头和平头的接触问题，根据接触宽度求解各处的位移和应力［16］。Jiang和Batra依据Eshelby-Stroh公式，改变了级数函数的形式，研究压头与内含缺陷的弹性体之间的接触问题［17］。

基于Eshelby-Stroh公式，可以求解边界条件较为复杂的接触问题，但是应力往往要比位移收敛得慢［18］，在级数项数较少时求解载荷值误差较大。为了获取较为准确的计算结果，计算过程中要设定较大的K值(函数中级数的项数)。本文提出使用整体位移约束法和线性迭代，提高求解结果的收敛性，在选择较小的K值时求得更加准确的计算结果。

1 接触问题分析

如图1所示，本文要研究的是各向异性的弹性体与凸面型压头接触作用的平面问题。各向异性材料的弹性体在x1轴上长L，在x2轴上高h。假设凸面型压头和弹性体在x3轴方向上(垂直于纸面向外)的长度相比于纸面上的参数(如长L、高h)足够大，凸面型压头受到载荷P，沿x2轴负方向垂直压入弹性体，记压头压入的深度为D，接触宽度为2c，压头与弹性体发生接触的高度为d，则有非接触高度u0=D－d。图中的凸面型压头是一个倒圆角的楔形。其中，楔形的斜度为θ，圆角半径为R。根据刚体压头有无与弹性体发生接触，可将弹性体按“竖直虚线”划分为Ⅰ、Ⅱ和Ⅲ三部分。

在直角坐标系下，弹性体各处的位移、应力和应变可以表示为关于x1和x2的函数。在不考虑弹性体自身重力时，存在以下关系:

图1 压头与弹性体的接触模型Fig.1 Contact model of indenter and elastic body

式中:应力 σij= σji，σij，j= σij/ xj;εks为应变;uk为在 xk方向上的位移;uk，s为 uk对 xs求导;Cijks为弹性体的材料参数，满足全对称性:

将式(3)代入式(2)并照式(1)求导可得

参照图1中的接触模型，弹性体的底部固定，左右两侧自由，顶部在中间长度2c区域发生接触。由此可得边界条件为

式中:μ为刚体压头与弹性体之间的摩擦系数; 为刚体压头与弹性体接触的切向角，该公式由切应力和正应力的平衡关系得到。若将刚体压头与弹性体之间的摩擦视为0，则有

对于图1中倒圆角的楔形，以区域Ⅰ和区域Ⅱ顶部的交接处为原点，x1和x2方向与图1中的坐标方向相同，则有外廓函数为

在研究小变形接触问题时， 的取值较小时便可以忽略，因此轴向载荷P可表示为

若考虑接触点处的斜率，轴向载荷P可表示为

2 接触问题的解析解

2.1 根据Eshelby-Stroh公式求解边界值问题

假设弹性体的位移、应力和应变是关于x1和x2的函数，根据Eshelby-Stroh公式可以求解上述边界值问题［19］。

首先假设

或

式中:f(z)为关于z的函数;ai可以根据材料参数求解;u和a分别为位移和位移参数。

其中:p为特征值。

将式(14)对xs求导可得

式中:δsi为克罗内克函数。再将式(17)对xj求导，后依据式(5)可得

即

将式(19)写成矩阵形式为

式中:Q、R、T都为3×3的矩阵，满足

求解式(20)可以引入

式中:ζ为特征向量。

p为复数，6个p中，有3对是共轭复数。假定

根据弹性体的材料参数，便可确定p和ai的值。根据式(2)和式(17)可得

因此，可以得位移和应力表达式分别为

式中:σ1、σ2为应力。u、σ1、σ2都是 3 ×1 的列向量，且

在区域Ⅰ与区域Ⅱ的交接处和区域Ⅱ与区域Ⅲ的交接处满足连续性条件:交接处左边的位移和应力与交接处右边位移和应力相等。将弹性体分成3个部分后，每一部分弹性体以其左下角为原点，不改变轴x1和轴x2的方向，则有

式中:u、σ1和 σ2右上角括号内的上标1、2、3对应弹性体Ⅰ、Ⅱ、Ⅲ。

假设每一部分组织的弹性体的f(z)函数为

则根据式(27)～式(29)可以得到各部分弹性体的位移函数和应力函数，如位移函数可以表示为

式中:

其中:、、都是3×1的列向量，其元素都为复数;conjugate表示式子中conjugate前面部分的共轭。

根据式(6)～式(9)给定的边界条件，则有

式中:、、和分别为弹性体上边、下边、左边和右边提取位移约束的3×3矩阵;、和分别为弹性体上边、下边、左边和右边提取应力约束的3×3矩阵。

为了满足式(35)的边界条件和式(30)的连续性条件，采用傅里叶积分的方式来进行计算:

式中:j=0，1，2，…;n=1，2，3。

通常情况下，非接触高度u0是未知量，可以通过求导的方式消去u0，区域Ⅱ顶部的边界条件式可以转化为

式中:

式(45)代入式(43)可得

为了保证函数在各边界、交接处的收敛效果相近，本文令各部分的k和m的上界满足:

式中:Ceil(*)表示大于或等于*的最小整数;K为所有K(n)相加之和。根据k和m的取值，选取合适的j值可求解()函数中的各个参数，便可分别求得三部分弹性体的位移和应力函数表达式。在实际计算过程中，k和m的取值不可能无限大，且求解到的位移的收敛性要比应力的收敛性好。并且在区域Ⅰ与区域Ⅱ交接处和区域Ⅱ与区域Ⅲ交接处往往存在应力突变，同时，非接触区域存在应力不近似于零的问题。这样求解到的函数结果不能很好地吻合最初设定的边界条件。为了解决这些问题，以下采取整体位移约束法和线性叠加原理，通过迭代方式获取更匹配最初边界条件的位移函数和应力函数。

2.2 整体位移约束法解决交接处应力突变问题

在划分成三部分弹性体计算时，在x2方向上，区域Ⅰ和区域Ⅲ顶部给定的是应力约束条件σ22=0，而区域Ⅱ给定的是位移约束条件 u2=g(x1)+u0。可见弹性体顶部存在不同形式的约束，导致交接处存在应力突变问题。为此，将三部分弹性体再视为整体进行考虑，指定以下约束条件:

并假设

此时的边界条件可以表示为

式中:

按照 2.1节中的傅里叶积分的方式对式(52)求解可以求解fα(zα)函数中的各个参数，便可求得弹性体的位移和应力函数表达式u、σ1、σ2。

2.3 线性叠加法解决非接触区域应力不为零的问题

通过整体法得到位移和应力函数后，在非接触区域仍然存在着正应力不近似于零的问题。此时分别提取区域Ⅰ和区域Ⅲ部分进行单独分析:先通过整体位移约束法的应力函数计算得到弹性体区域Ⅰ和区域Ⅲ顶部的应力σ12和σ22，之后将得到的应力σ12和σ22取反后分别作为单独计算时区域Ⅰ和区域Ⅲ顶部的应力约束;区域Ⅰ和区域Ⅲ底部固定，在x1、x2和x3方向上给定位移约束;对于区域Ⅰ和区域Ⅲ左右两侧，在x1和x2方向上给定应力约束，在x3方向上给定位移约束。因此，此时弹性体区域Ⅰ和区域Ⅲ位移、应力函数计算的边界条件为

其中:n=1，3。

同样依照前面的计算方式求解出()函数中的各个参数，得到⌒式(54)边界条件下位移、应力的函数表达式为、。然后将其与之前整体位移约束法计算的结果、、行叠加，并将叠加的结果作为边界条件的一部分进行再一次的整体计算，与2.2节中边界条件相同，但是顶部位移约束条件的数值发生改变:

依据边界条件和上述的位移约束条件函数求解出位移、应力表达式u、σ1和σ2。到此，便完成了线性叠加。

2.4 迭代法逼近理想解

通过2.2节中的整体位移约束法计算和2.3节中的线性叠加法计算，非接触区域的应力在一定程度上变小，但还是不为零。由于应力函数给定的是傅里叶级数的函数形式，且k和m的取值有限，因此不能保证弹性体顶部非接触区域每个点的应力满足σ12=0和σ22=0，只能使其值近似于零。

对于σ12和σ22不近似于零时，需要将2.3节计算出来的u、σ1和σ2记为，然后重复2.3节中的线性叠加法，在区域Ⅰ和区域Ⅲ单独计算出叠加后作为整体计算时的边界条件，再次得到u、σ1和σ2。最后再次判断是否σ12≈0和σ22≈0，若不满足则再次重复上述操作，进行迭代。在计算出线性叠加后的位移函数和应力函数后，可通过每次迭代后对应力改变量是否小于一个值来判定迭代是否结束。

3 算法验证

3.1 圆柱压头接触问题验证

为了验证算法的正确性以及优越性，参考文献［16］中的数据给定参数:

式中:E1、E2、E3分别为坐标轴3个方向的弹性模量;G23、G31、G12分别为坐标轴3个方向的剪切模量;ν23、ν31、ν12分别为坐标轴 3 个方向的泊松比。

将刚体压头的外廓修改成圆柱形，则外廓函数为

在MATLAB中编写程序进行计算，再与文献［14-15］的解析解(弹性半空间下的求解)做比较。图2给出了K=500时顶部应力分布与文献［14-15］弹性半空间法求解的应力分布。可以看到，在非接触区域，当前解与弹性半空间解一致;在接触区域，当前解也与弹性板空间解较吻合。图3给出了K=500时接触区域的相对误差εre，进一步说明它的吻合程度。在图3中，竖直2条虚线的中间区域代表了弹性体顶部

x1－L /2 ＜0.9c区域，在此区域内的最大相对误差为2.6%。

图2 当前解与弹性半空间解的比较Fig.2 Comparison between present solution andsolution of elastic half space

图3 K=500时接触区域应力的相对误差Fig.3 Relative error of stress in contact zone at K=500

表1 K=400时的迭代结果Table 1 Iterative results at K=400

得到应力分布函数后，通过积分求得载荷值，表1为K取400时的迭代结果。随着迭代次数增加，中间点的位移与应力、载荷的变化越来越小。经过15次迭代后载荷值为5305.60N。用文献［14-15］的方法计算得到的载荷 Portho为5278.14 N，采用式(59)计算相对误差，并与文献［16］的计算结果进行比较。在文献［16］的计算结果中，K=400时相对误差为 1.9%，K=1000时相对误差为0.8%，而本文在K=400时计算得到的相对误差仅为0.52%。可见，本文计算结果的相对误差更小，且在K取较小值时便能获得较优的结果。

3.2 圆柱压头与有界弹性体接触问题验证

该求解方法与弹性半空间平面接触问题求解不同，可以适用于求解有界弹性体的接触问题。

给定参数如下:

取K=400，求解得到刚体压头的压痕深度为4.788mm(实际收敛结果为4.78779mm)，载荷值为5 832.3 N。在ABAQUS软件中依据上述的参数建立模型，由于ABAQUS中依据接触宽度作为边界条件比较困难，因此采用给定压痕深度，反过来求接触宽度。

在ABAQUS中给予刚体压头－4.788 mm的位移约束，仿真得到的接触应力如图4所示，求得接触宽度为100.33 mm，载荷值为5 793.6 N。容易得到，计算的接触宽度对于仿真结果的相对误差为0.33%，计算得到的载荷值相对于仿真结果的相对误差为0.67%。图5显示了弹性体接触区域应力分布的计算结果和仿真结果，定义绝对误差εab为计算值减仿真值，则最大绝对误差为3.648 MPa;除去右边靠近交接处的位置，其他位置的绝对误差值都小于3 MPa。可见，本文求解方法在解决有界弹性体接触问题上具有可行性。

图4 ABAQUS仿真接触应力Fig.4 Contact stress of ABAQUS simulation

图5 圆柱压头与有界弹性体接触应力分布Fig.5 Distribution of contact stress between cylindrical indenter and bounded elastic body

3.3 倒圆角楔形压头接触问题验证

在图1所示模型中，刚体压头是倒圆角的楔形，式(11)是它的外廓函数。该外廓函数对圆角半径R=0时也适用。基于式(60)中弹性体的参数，给定楔形的参数:楔形斜度θ和圆角半径R。在 θ取 10°时，R 分别取 0、50、100、150、200、250 mm进行计算，得到压痕深度后作为位移约束在ABAQUS中进行对应的有限元仿真，得到表2、表3和图6。

表2 不同R值下的压痕深度与仿真接触宽度Table 2 Indentation depth and simulation contact width at different R values

表2和表3表示在K=500时不同R值下的计算结果和仿真结果，并给出载荷值和接触宽度的相对误差。如表2所示，在有圆角的情况下，随着R值的增加，达到相同的接触宽度(2c=100 mm)，需要的压痕深度值逐渐减少;如表3所示，随着R值的增加，载荷值也逐渐减少。载荷值的相对误差都小于2%，接触宽度2c相对误差都小于1%。可见，K=500时计算值和其对应的仿真值是比较接近的。

表3 不同R值下的计算载荷与仿真载荷Table 3 Calculation load and simulation load at different R values

图6 在θ=10°，R=0，50，100，150，200，250 mm时的计算与仿真应力分布及两者的绝对误差Fig.6 Computed and simulated stress distribution and absolute error between them at θ=10°and R=0，50，100，150，200，250 mm

图6 为 R 取0、50、100、150、200、250mm 时计算得到的应力分布和仿真得到的接触应力分布以及两者的绝对误差。在图6(a)～图6(f)中，坐标左右两端的误差比较大，一个原因在于左右两端是接触区域和非接触区域的交接处，交接处的应力计算时不容易收敛且此处附近的点在x1轴方向上的位移变化相对中间区域的点要大。实际上，算法中给顶部的位移约束不是刚体压头的外廓函数减间隙高度，因为接触点处不光发生了x2轴方向的位移，也在x1轴方向上也发生位移，因此实际给的外廓函数在x1轴方向上发生了形变。接触区域各点在x1方向发生的移动指向接触中心处，计算用的刚体压头相对于仿真用的刚体压头要略微“偏瘦”。另一个原因在于原本设定的级数函数就是为了逼近实际的位移函数和应力函数，求得的函数在理想的函数上就有一些波动。算法本身便会导致误差，在斜率较大或突变的地方会更容易造成较大的误差。

另外，在切点(圆弧和斜线的过度处)附近的绝对误差值也要比正常地方的绝对误差值大。例如，在 R=150 mm 时，Rsin 10°=26.04 mm，在图6(d)中可以看出在x1－L/2=±26.04 mm附近绝对误差值较大;在R=200 mm时，Rsin 10°=34.73 mm，在图 6(e)中可以看出 x1－L/2=±34.73 mm附近绝对误差也较大。因此，随着R值的增加，绝对误差较大处逐渐往外移动。

对于有圆角的楔形，从图6(b)～图6(f)可以看出，在接触区域，计算值与仿真值的绝对误差值小于10 MPa，除去交接点附近和切点附近，绝对误差值小于5 MPa。

虽然接触应力的绝对误差并不小，但计算得到的应力函数在仿真得到的应力分布曲线上波动，由于误差有正有负，能够得到消除，使得最后得到的载荷值误差较小。

4 结论

1)针对线性各向异性弹性体小变形接触问题，本文基于Eshelby-Stroh公式，运用线性材料叠加原理和整体位移约束法，解决了圆柱型压头和倒圆角楔形压头与弹性体的接触问题。

2)针对正交各向异性弹性体接触问题，本文以圆柱压头为例，求得的结果在载荷值比文献［16］求解结果更接近于文献［14-15］的求解结果，且在K取较小值时便能算出相对误差较小的载荷值:K=400时载荷计算值相对于文献［14-15］弹性半空间计算结果的相对误差为0.52%。本文算法更佳。

3)本文算法解决了圆柱压头与有界线性正交各向异性弹性体的接触问题，并用ABAQUS仿真验证:计算的接触宽度和载荷对于仿真结果的相对误差分别为0.33%、0.67%。

4)对于倒圆角楔形压头与线性正交各向异性弹性体小变形接触问题，在接触区域和非接触区域交接处的绝对误差相对普通接触点较大。同时在刚体压头外廓上，圆弧和斜线过度处也会发生应力绝对误差较大的现象，并随着R值的增加，绝对误差较大处逐渐往外移动。至于普通接触点的接触应力，计算值和仿真值比较接近。虽然某些地方接触应力的绝对误差并不小，但通过积分得到载荷值的相对误差都小于2%。因此，已知倒圆角楔形压头与线性正交各向异性弹性体的接触宽度，通过本算法求解压头的载荷是可行的。

［1］ T YR S J，LYYRA-LAITINEN T，NIINIM K M，et al.Estimation of the Young’s modulus of articular cartilage using an arthroscopic indentation instrument and ultrasonic measurement of tissue thickness［J］.Journal of Biomechanics，2001，34(2):251-256.

［2］ KORHONEN R K，SAARAKKALA S，TOYRAS J，et al.Experimental and numerical validation for the novel configuration of an arthroscopic indentation instrument［J］.Physics in Medicine ＆Biology，2003，48(11):1565-1576.

［3］ DIMITRIADIS E K，HORKAY F，MARESCA J，et al.Determination of elastic moduli of thin layers of soft material using the atomic force microscope［J］.Biophysical Journal，2002，82(5):2798-2810.

［4］ WANNINAYAKE I B，DASGUPTA P，SENEVIRATNE L D，et al.Air-float palpation probe for tissue abnormality identification during minimally invasive surgery［J］.IEEE Transactions on Biomedical Engineering，2013，60(10):2735-2744.

［5］郭渊，关志东，刘德博，等.复合材料静压痕与落锤冲击初始损伤对比试验［J］.北京航空航天大学学报，2009，35(8):1018-1021.GUO Y，GUANG Z D，LIU D B，et al.Comparison between quasi-static indentation testing and drop-weight impact testing on delamination on set damage［J］.Journal of Beijing University of Aeronautics and Astronautics，2009，35(8):1018-1021(in Chinese).

［6］ STOLZ M，GOTTARDI R，RAITERI R，et al.Early detection of aging cartilage and osteoarthritis in mice and patient samples using atomic force microscopy［J］.Nature Nanotechnology，2009，4(3):186-192.

［7］ LIAO Q，HUANG J，ZHU T，et al.A hybrid model to determine mechanical properties of soft polymers by nanoindentation［J］.Mechanics of Materials，2010，42(12):1043-1047.

［8］ MUSKHELISHVILI N I.Some basic problems of the mathematical theory of elasticity［J］.Mathematical Gazette，1953，48(365):351.

［9］ ENGLAND A H，SIH G C.Complex variable methods in elasticity［M］.London:Wiley-Interscience，1971:318.

［10］ GLADWELL G M L，ENGLAND A H.Orthogonal polynomial solutions to some mixed boundary-value problems in elasticity theory［J］.Quarterly Journal of Mechanics ＆ Applied Mathematics，1977，30(2):175-185.

［11］ OKUMURA M，HASEBE N，NAKAMURA T.Crack due to wedge-shaped punch with friction［J］.Journal of Engineering Mechanics，1990，116(10):2173-2185.

［12］ STROH A N.Dislocations and cracks in anisotropic elasticity［J］.Philosophical Magazine，1958，3(30):625-646.

［13］ STROH A N.Steady state problems in anisotropic elasticity［J］.Studies in Applied Mathematics，1962，41(1):77-103.

［14］ FAN C W，HWU C.Punch problems for an anisotropic elastic half-plane［J］.Journal of Applied Mechanics，1996，63(1):69-76.

［15］ HWU C，FAN C W.Sliding punches with or without friction along the surface of an anisotropic elastic half-plane［J］.Quarterly Journal of Mechanics ＆ Applied Mathematics，1998，51(1):159-177.

［16］ BATRA R C，JIANG W.Analytical solution of the contact problem of a rigid indenter and an anisotropic linear elastic layer［J］.International Journal of Solids and Structures，2008，45(22-23):5814-5830.

［17］ JIANG W，BATRA R C.Indentation of a laminated composite plate with an interlayer rectangular void［J］.Composites Science and Technology，2010，70(6):1023-1030.

［18］ VEL S S，BATRA R C.The generalized plane strain deformations of thick anisotropic composite laminated plates［J］.International Journal of Solids and Structures，2000，37(5):715-733.

［19］ TINGT C T.Anisotropic elasticity theory and applications［M］.Oxford:Oxford University Press，1996:134-142.