CN107329934A - 一种凸面型压头与弹性体接触的分析方法 - Google Patents

Abstract

本发明公开一种凸面型压头与弹性体接触问题的求解方法，属于弹性力学领域。该方法针对的弹性体主要是正交各向材料的弹性体，不过此法对各向同性材料的弹性体也适用。本发明主要研究平面问题的接触模型，根据Eshelby‑Reid‑Shockley公式、线性材料叠加原理，求解凸面型压头与弹性体的接触问题。相对于一般的求解方法，该方法具有可以求解有界弹性体问题、凸面型的压头可以复杂化等特点。已知材料参数的情况下，给定接触面积可以求解出反作用力、应力分布及弹性体各点的位移等。因此本发明所提出的求解方法能为一些实际接触问题提供指导依据，进而通过触觉辨识所接触的组织。

Description

一种凸面型压头与弹性体接触的分析方法

技术领域

本发明涉及弹性力学领域，阐述了一种凸面型压头与弹性体接触问题的求解方法。

背景技术

随着临床医疗技术向高、准、细方向发展，微创手术在临床治疗过程中扮演着越来越重要的角色，而触觉传感对于微创手术的成功实施作用明显。触觉传感一般要通过触觉传感器来实现，通过触觉传感器测量力、力的分布、传感器压入深度等数据计算得到弹性材料的弹性性能、变形难易程度。触觉信息求解的过程便涉及到了弹性材料的接触问题。

目前，国内外学者根据弹性组织的特点和相关边界条件，对圆锥体、圆柱体、球体等相关规则体的接触受力进行了详细分析，建立相关规则体的弹性接触模型。由于规则体对组织接触作用所造成的变形呈现一定的规律性，这为相关规则体的接触建模分析提供了依据。一般研究较多的是规则体的接触，而且其中大部分都局限于弹性半空间体，故最后得到的结论准确度不高，无法直观的得到弹性材料的接触状态。

发明内容

本发明所要解决的技术问题是凸面型压头与弹性体接触问题。本发明提供一种新的技术方案：一种凸面型压头与弹性体接触的分析方法，其特征在于，包括以下步骤：

1)将弹性体按是否与压头发生接触进行划分为若干部分；

2)根据初始的边界条件给定划分的若干部分弹性体对应的边界条件，同时划分的交接处要满足连续性条件，所述若干部分包括与压头接触的接触区域及与压头不接触的非接触区域；

3)根据边界条件和连续性条件使用“D-S计算”求解弹性体划分的各个部分的位移函数和应力函数；

4)将若干部分弹性体顶部的位移函数合并成一个整体的位移函数；

5)将步骤4)中合并得到位移函数作为整个弹性体的边界条件进行“D-S计算”，得到整个弹性体的位移函数和应力函数；

6)根据步骤5)的应力函数判断非接触区域的顶部正应力和切应力是否接近于0；若是，执行步骤10；若不是，执行步骤7；

7)将非接触区域顶部应力分别取反后作为各自的边界条件进行“D-S计算”，得到各个非接触弹性体的位移函数和应力函数；

8)将步骤5)与步骤7)非接触区域顶部的位移函数叠加作为非接触区域顶部的位移，并计算非接触区域与接触区域交接点的位移，将该交接点的位移添加到压头的外廓函数作为接触区域的位移；

9)回到步骤4)执行算法；

10)根据位移函数和应力函数求需要的量。

进一步的，上述的分析方法应用于触觉识别。

有益效果：在一般的压痕试验中，常常选用的压头为圆柱型、平头等较为规则的几何外廓，基于本专利的技术，在压痕试验中可以采用更多样式的压头，如抛物线型、楔形、倒圆角的楔形以及其他组合形式的凸型压头。本专利的技术可以应用于触觉识别：获取压头的接触宽度及其所受的压力，可以辨别弹性体的变形的难易程度。在已知弹性体的材料参数的情况下，同时测得压头与弹性体之间的接触面积，可以推算出压头所受的作用力、应力分布以及弹性体上个点的位移。本专利的技术适用于复杂边界条件问题的求解，能求解有界弹性体与压头接触的问题。一般的接触问题求解是将弹性体考虑成弹性半空间体进行计算，而本专利的求解方法将弹性体作为有界弹性体进行计算，其更接近实际情况。

进一步的，所述步骤3)中的D-S计算包括以下步骤：

3.1)分别设定所述若干部分的f(z)函数；

3.2)依据位移和应力与f(z)函数的关系确立位移函数和位移函数关于z的表达式；

3.3)根据边界条件和连续性条件获得方程式；

3.4)对方程式进行傅里叶积分并进行求解；

3.5)将求解得到的参数代回位移、应力关于z的关系式得到位移函数和应力函数。

进一步的，所述步骤7)中的D-S计算包括以下步骤：

7.1)设定弹性体整体的f(z)函数；

7.2)依据位移和应力与f(z)函数的关系确立位移函数和位移函数关于z的表达式；

7.3)根据边界条件获得方程式；

7.4)对方程式进行傅里叶积分并进行求解；

7.5)将求解得到的参数代回位移、应力关于z的关系式得到位移函数和应力函数。

附图说明

图1是本发明的流程示意图；

图2是本发明的D-S计算流程示意图；

图3是本发明的整体模型图；

图4是本发明的弹性体拆分成三部分时的边界条件图；

图5是本发明的三部分弹性体重新看作整体时的边界条件图；

图6是本发明的在修正结果未收敛时part 1的边界条件图；

图7是本发明的在修正结果未收敛时part 3的边界条件图；

图8是本发明的整体计算结果与part 1和part 3单独计算结果的叠加图。

具体实施方式

下面结合附图对所发明的技术方案进行详细说明：

本发明的流程示意图如图1所示，步骤如下：

1)将弹性体按是否与压头发生接触进行划分为若干部分；

2)根据初始的边界条件给定划分的若干部分弹性体对应的边界条件，同时划分的交接处要满足连续性条件，所述若干部分包括与压头接触的接触区域及与压头不接触的非接触区域；

3)根据边界条件和连续性条件使用“D-S计算”求解弹性体划分的各个部分的位移函数和应力函数；

4)将若干部分弹性体顶部的位移函数合并成一个整体的位移函数；

5)将步骤4)中合并得到位移函数作为整个弹性体的边界条件进行“D-S计算”，得到整个弹性体的位移函数和应力函数；

6)根据步骤5)的应力函数判断非接触区域的顶部正应力和切应力是否接近于0；若是，执行步骤10；若不是，执行步骤7；

7)将非接触区域顶部应力分别取反后作为各自的边界条件进行“D-S计算”，得到各个非接触弹性体的位移函数和应力函数；

8)将步骤5)与步骤7)非接触区域顶部的位移函数叠加作为非接触区域顶部的位移，并计算非接触区域与接触区域交接点的位移，将该交接点的位移添加到压头的外廓函数作为接触区域的位移；

9)回到步骤4)执行算法；

10)根据位移函数和应力函数求需要的量，如压痕深度、最大应力等，也可通过积分求得作用力。

在图1的流程示意图和上述步骤中提到的“D-S计算”如图2所示。为了便于说明这个计算方法，以下做些基本说明和公式推导。

如图3所示，这是一个楔形凸面与弹性体接触模型，研究的弹性体是正交各向异性材料且弹性体底面固定，同时不考虑摩擦力的因素且楔形的斜度较小。以纸面的横向为x₁轴方向，纸面的纵向为x₂轴方向，垂直于纸面向外方向为x₃轴方向，并假设弹性体在x₃轴方向上足够长，从而可将接触问题视为x_1-x2平面问题。将弹性体分成三个部分：part 1(图中I区域)、part 2(图中II区域)和part3(图中III区域)，其中只有part 2与楔形凸面发生接触。其中，弹性体的宽度为L，高度为h，刚体压头与弹性体发生的接触长度为2c，接触高度为d，压头从刚接触到最后的稳定状态在x₂方向的移动距离为D，接触间隙高度u₀＝D-d。不考虑弹性体的自身重力，弹性体满足以下条件：

σ_ij，j＝0，(i，j＝1，2，3) (1)

σ_ij＝C_ijksε_ks，(k，s＝1，2，3) (2)

其中，应力σ_ij＝σ_ji，ε_ks表示弹性体应变，u_k、u_s分别表示弹性体在x_k、x_s方向上的位移，C_ijks是弹性体的材料参数，并满足全对称性：

C_ijks＝C_jiks＝C_ksij＝C_ijsk (4)

将式(3)代入式(2)后求导可得

C_ijksu_k，sj＝0，(i，j＝1，2，3；k，s＝1，2，3) (5)

根据Eshelby-Reid-Shockley公式，设定u_i＝a_if(z)或者u＝af(z)，其中z＝x₁+px₂，则有

u_k，s＝(δ_s1+pδ_s2)a_kf′(z) (6)将其代入到式(5)可得

C_ijks(δ_j1+pδ_j2)(δ_s1+pδ_s2)a_k＝0

(7)

或{C_i1k1+p(C_i1k2+C_i2k1)+p²C_i2k2}a_k＝0

可以写成

{Q+p(R+R^T)+p²T}a＝0 (8)

其中，Q、R、T都是3x3的矩阵：Q_ik＝C_i1k1，R_ik＝C_i1k2，T_ik＝C_i2k2。对于正交各向异件材料，则有：

可以看出在应变能为正的时候，Q和T是对称矩阵且是正定的。想要式子中的a有非零解，必须要满足

|Q+p(R+R^T)+p²T|＝0 (12)

定义

对于这个式子的求解，可以引入矩阵N使得

Nζ＝pζ (14)

其中，

p是特征值，ζ是特征向量。p是复数，在6个特征值p中有三对是共轭复数。假定Im(p_α)＞0，是p_α的共轭复数。便可获得位移u的通解：

同样地，应力的通解可以写成：

即

其中，(σ₁)_i＝σ_i1，(σ₂)_i＝σ_i2。

假定函数

其中，

则位移和应力可表示成：

其中，A＝[a₁ a₂ a₃]，B＝[b₁ b₂ b₃]， <φ*Ψ*χ*>＝diag[φ₁ψ₁χ₁，φ₂Ψ₂χ₂，φ₃Ψ₃χ₃]，

至此，三个弹性体的应力和位移可以表示为：

在有了以上的基础后，便可以计算弹性体分3部分时各自的位移和应力。将弹性体分成3部分之后的边界条件如图4所示，由于将接触问题视为平面问题，故始终存在

u₃＝0 (25)

又由于part 1和part 3的顶部未与压头发生接触，则这两部分满足

σ₁₂＝0，σ₂₂＝0 (26)

弹性体part 2顶部与压头发生接触，忽略摩擦力的因素，则有

σ₁₂＝0，u₂＝u₂(x₁)＝g(x₁)+u₀ (27)

其中，g(x₁)是压头的外廓函数。在本发明中是以凸型斜面为例，若只针对part 2分析，外廓函数为

若针对整个弹性体分析，外廓函数为

Part 1、part 2和part 3的底面是固定的，故满足边界条件：

u₁＝0，u₂＝0 (30)

Part 1左侧和part 3右侧也不发生接触，故满足边界条件：

σ₁₁＝0，σ2₁＝0 (31)

以上的边界条件总结为

Part 1右侧与part 2左侧相接，part 2右侧与part3左侧相接，因此在这要满足位移和应力的连续性，固有约束条件：

转换到各个部分，可表示为

将边界条件和连续性条件的式子进行傅里叶积分，得到

其中，l⁽¹⁾＝l⁽³⁾＝L/2-c，l⁽²⁾＝2c，

若将u₀作为一个未知量来考虑，需要通过求导来消去u₀，part 2的顶部的边界条件式子转化为

其中，

根据之前的k和m值选取合适的j值(保证线性方程组35～40有唯一解)求解出式中的各个参数，便可得到三部分弹性体的位移、应力表达式。从“设定”到求解出弹性体的位移应力的这一过程，在本发明中称为“D-S计算”，流程如图2所示，上述的计算是将弹性体划分为三个区域的D-S计算。

由于在实际计算过程中，k和m值不可能无限大，求解得到的结果往往位移收敛性较好，而应力并没有很好地收敛。在part 1与part 2交接处和part 2与part3交接处存在着应力突变以及应力不近似于零等问题，为了解决这些问题，依据线性叠加原理使用迭代法得到更好的函数。

对于弹性体的顶部，在x₂方向上，part 1和part 3给定的是位移约束σ₂₂＝0，而part 2给定的是位移约束u₂＝u₂(x₁)＝g(x₁)+u₀。为了将三个部分的弹性体重新作为一个弹性体考虑，弹性体顶部需要给定一样的位移约束，具体的边界条件如图5所示。其中，

假设整体分析时

其中，z_α＝x₁+p_αx₂，

与三部分相同，得到整体的位移、应变表达式：

u(0～L，0～h)、σ₁(0～L，0～h)、σ₂(0～L，0～h)。

边界条件为

其中，

对这四个边界条件式子进行傅里叶积分，得到

其中，g^t(x₁)＝[0，U₂(x₁)，0]^T，g^b(x₁)＝g^l(x₂)＝g^r(x₂)＝[0，0，0]^T。

同样根据k和m的值设定合适的j值来求解方程组，求得f_α(z_α)各项参数，便可得到整个弹性体位移、应力的函数表达式在交接处的位移和应力的连续性得到了很大改善，但仍会存在非接触区域应力不近似于零的现象。这时，需要提取出part 1和part3进行处理，将整体计算得到的σ₁₂、σ₂₂反向后最为part 1和part 3的应力约束，边界条件如图6和图7所示，特别需要注意的是σ₁₂＝-σ₁₂和σ₂₂＝-σ₂₂是赋值语句，是将整体计算得到的结果取反后作为边界条件。

其中，n＝1，3， 和表示整体计算时得到的应力。对上述式子进行傅里叶积分得到

其中，同样根据k和m的值设定合适的j值来求解方程组，求得各项参数，便可得到图6和图7边界条件下位移、应力的函数表达式

然后将part 1和part 3的计算结果与之前整体计算的结果进行叠加，叠加之后的边界约束如图8所示。

然后再进行一次整体计算，边界条件与之前整体计算的边界条件相同，如图5所示，但顶面的位移约束条件发生了改变：

完成本次整体计算之后判断非接触区域(part 1和part3)顶部各个部位是否满足σ₁₂≈0、σ₂₂≈0。由于σ₂是傅里叶级数的形式，且k和m的取值有限，很难保证顶部的每个点都满足σ₁₂＝0、σ₂₂＝0，因此需要设定合适的阈值判断其是否近似为0。通常情况下，在part 1、part 3与part 2相接处附近的各点(不包括part 2上的点)满足σ₁₂≈0和σ₂₂≈0时，便可认为得到的f_α(z_α)是满足我们最初边界条件的函数。假如不满足这两个条件，需要进行迭代运算，重复“分别提取part 1和part 3并给予反向应力计算”、“单独计算结果叠加到整体计算结果后作为顶部位移约束进行整体计算”。

在得到满足要求的f_α(z_α)函数后，将f_α(z_α)代入式(16)、式(19)和式(20)，便得到了弹性体各点的位移函数和应力函数。再对part 2顶部的σ₂₂进行积分便能得到作用力P：

若考虑压头的外廓，则有

其中，θ是接触点切线与坐标轴x₁的夹角。

Claims (4)

1.一种凸面型压头与弹性体接触的分析方法，其特征在于，包括以下步骤：

1)将弹性体按是否与压头发生接触进行划分为若干部分；

2)根据初始的边界条件给定划分的若干部分弹性体对应的边界条件，同时划分的交接处要满足连续性条件，所述若干部分包括与压头接触的接触区域及与压头不接触的非接触区域；

3)根据边界条件和连续性条件使用“D-S计算”求解弹性体划分的各个部分的位移函数和应力函数；

4)将若干部分弹性体顶部的位移函数合并成一个整体的位移函数；

5)将步骤4)中合并得到位移函数作为整个弹性体的边界条件进行“D-S计算”，得到整个弹性体的位移函数和应力函数；

6)根据步骤5)的应力函数判断非接触区域的顶部正应力和切应力是否接近于0；若是，执行步骤10；若不是，执行步骤7；

7)将非接触区域顶部应力分别取反后作为各自的边界条件进行“D-S计算”，得到各个非接触弹性体的位移函数和应力函数；

8)将步骤5)与步骤7)非接触区域顶部的位移函数叠加作为非接触区域顶部的位移，并计算非接触区域与接触区域交接点的位移，将该交接点的位移添加到压头的外廓函数作为接触区域的位移；

9)回到步骤4)执行算法；

10)根据位移函数和应力函数求需要的量。

2.如权利要求1所述的分析方法，其特征在于：所述步骤3)中的D-S计算包括以下步骤：

3.1)分别设定所述若干部分的f(z)函数；

3.2)依据位移和应力与f(z)函数的关系确立位移函数和位移函数关于z的表达式；

3.3)根据边界条件和连续性条件获得方程式；

3.4)对方程式进行傅里叶积分并进行求解；

3.5)将求解得到的参数代回位移、应力关于z的关系式得到位移函数和应力函数。

3.如权利要求1所述的分析方法，其特征在于：所述步骤7)中的D-S计算包括以下步骤：

7.1)设定弹性体整体的f(z)函数；

7.2)依据位移和应力与f(z)函数的关系确立位移函数和位移函数关于z的表达式；

7.3)根据边界条件获得方程式；

7.4)对方程式进行傅里叶积分并进行求解；

7.5)将求解得到的参数代回位移、应力关于z的关系式得到位移函数和应力函数。

4.一种如权利要求1至3中任一项所述的分析方法的应用，其特征在于：应用于触觉识别。