本發(fā)明涉及彈性力學(xué)領(lǐng)域，闡述了一種凸面型壓頭與彈性體接觸問題的求解方法。

背景技術(shù)：

隨著臨床醫(yī)療技術(shù)向高、準(zhǔn)、細(xì)方向發(fā)展，微創(chuàng)手術(shù)在臨床治療過程中扮演著越來越重要的角色，而觸覺傳感對于微創(chuàng)手術(shù)的成功實(shí)施作用明顯。觸覺傳感一般要通過觸覺傳感器來實(shí)現(xiàn)，通過觸覺傳感器測量力、力的分布、傳感器壓入深度等數(shù)據(jù)計算得到彈性材料的彈性性能、變形難易程度。觸覺信息求解的過程便涉及到了彈性材料的接觸問題。

目前，國內(nèi)外學(xué)者根據(jù)彈性組織的特點(diǎn)和相關(guān)邊界條件，對圓錐體、圓柱體、球體等相關(guān)規(guī)則體的接觸受力進(jìn)行了詳細(xì)分析，建立相關(guān)規(guī)則體的彈性接觸模型。由于規(guī)則體對組織接觸作用所造成的變形呈現(xiàn)一定的規(guī)律性，這為相關(guān)規(guī)則體的接觸建模分析提供了依據(jù)。一般研究較多的是規(guī)則體的接觸，而且其中大部分都局限于彈性半空間體，故最后得到的結(jié)論準(zhǔn)確度不高，無法直觀的得到彈性材料的接觸狀態(tài)。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明所要解決的技術(shù)問題是凸面型壓頭與彈性體接觸問題。本發(fā)明提供一種新的技術(shù)方案：一種凸面型壓頭與彈性體接觸的分析方法，其特征在于，包括以下步驟：

1)將彈性體按是否與壓頭發(fā)生接觸進(jìn)行劃分為若干部分；

2)根據(jù)初始的邊界條件給定劃分的若干部分彈性體對應(yīng)的邊界條件，同時劃分的交接處要滿足連續(xù)性條件，所述若干部分包括與壓頭接觸的接觸區(qū)域及與壓頭不接觸的非接觸區(qū)域；

3)根據(jù)邊界條件和連續(xù)性條件使用“d-s計算”求解彈性體劃分的各個部分的位移函數(shù)和應(yīng)力函數(shù)；

4)將若干部分彈性體頂部的位移函數(shù)合并成一個整體的位移函數(shù)；

5)將步驟4)中合并得到位移函數(shù)作為整個彈性體的邊界條件進(jìn)行“d-s計算”，得到整個彈性體的位移函數(shù)和應(yīng)力函數(shù)；

6)根據(jù)步驟5)的應(yīng)力函數(shù)判斷非接觸區(qū)域的頂部正應(yīng)力和切應(yīng)力是否接近于0；若是，執(zhí)行步驟10；若不是，執(zhí)行步驟7；

7)將非接觸區(qū)域頂部應(yīng)力分別取反后作為各自的邊界條件進(jìn)行“d-s計算”，得到各個非接觸彈性體的位移函數(shù)和應(yīng)力函數(shù)；

8)將步驟5)與步驟7)非接觸區(qū)域頂部的位移函數(shù)疊加作為非接觸區(qū)域頂部的位移，并計算非接觸區(qū)域與接觸區(qū)域交接點(diǎn)的位移，將該交接點(diǎn)的位移添加到壓頭的外廓函數(shù)作為接觸區(qū)域的位移；

9)回到步驟4)執(zhí)行算法；

10)根據(jù)位移函數(shù)和應(yīng)力函數(shù)求需要的量。

進(jìn)一步的，上述的分析方法應(yīng)用于觸覺識別。

有益效果：在一般的壓痕試驗(yàn)中，常常選用的壓頭為圓柱型、平頭等較為規(guī)則的幾何外廓，基于本專利的技術(shù)，在壓痕試驗(yàn)中可以采用更多樣式的壓頭，如拋物線型、楔形、倒圓角的楔形以及其他組合形式的凸型壓頭。本專利的技術(shù)可以應(yīng)用于觸覺識別：獲取壓頭的接觸寬度及其所受的壓力，可以辨別彈性體的變形的難易程度。在已知彈性體的材料參數(shù)的情況下，同時測得壓頭與彈性體之間的接觸面積，可以推算出壓頭所受的作用力、應(yīng)力分布以及彈性體上個點(diǎn)的位移。本專利的技術(shù)適用于復(fù)雜邊界條件問題的求解，能求解有界彈性體與壓頭接觸的問題。一般的接觸問題求解是將彈性體考慮成彈性半空間體進(jìn)行計算，而本專利的求解方法將彈性體作為有界彈性體進(jìn)行計算，其更接近實(shí)際情況。

進(jìn)一步的，所述步驟3)中的d-s計算包括以下步驟：

3.1)分別設(shè)定所述若干部分的f(z)函數(shù)；

3.2)依據(jù)位移和應(yīng)力與f(z)函數(shù)的關(guān)系確立位移函數(shù)和位移函數(shù)關(guān)于z的表達(dá)式；

3.3)根據(jù)邊界條件和連續(xù)性條件獲得方程式；

3.4)對方程式進(jìn)行傅里葉積分并進(jìn)行求解；

3.5)將求解得到的參數(shù)代回位移、應(yīng)力關(guān)于z的關(guān)系式得到位移函數(shù)和應(yīng)力函數(shù)。

進(jìn)一步的，所述步驟7)中的d-s計算包括以下步驟：

7.1)設(shè)定彈性體整體的f(z)函數(shù)；

7.2)依據(jù)位移和應(yīng)力與f(z)函數(shù)的關(guān)系確立位移函數(shù)和位移函數(shù)關(guān)于z的表達(dá)式；

7.3)根據(jù)邊界條件獲得方程式；

7.4)對方程式進(jìn)行傅里葉積分并進(jìn)行求解；

7.5)將求解得到的參數(shù)代回位移、應(yīng)力關(guān)于z的關(guān)系式得到位移函數(shù)和應(yīng)力函數(shù)。

附圖說明

圖1是本發(fā)明的流程示意圖；

圖2是本發(fā)明的d-s計算流程示意圖；

圖3是本發(fā)明的整體模型圖；

圖4是本發(fā)明的彈性體拆分成三部分時的邊界條件圖；

圖5是本發(fā)明的三部分彈性體重新看作整體時的邊界條件圖；

圖6是本發(fā)明的在修正結(jié)果未收斂時part1的邊界條件圖；

圖7是本發(fā)明的在修正結(jié)果未收斂時part3的邊界條件圖；

圖8是本發(fā)明的整體計算結(jié)果與part1和part3單獨(dú)計算結(jié)果的疊加圖。

具體實(shí)施方式

下面結(jié)合附圖對所發(fā)明的技術(shù)方案進(jìn)行詳細(xì)說明：

本發(fā)明的流程示意圖如圖1所示，步驟如下：

1)將彈性體按是否與壓頭發(fā)生接觸進(jìn)行劃分為若干部分；

2)根據(jù)初始的邊界條件給定劃分的若干部分彈性體對應(yīng)的邊界條件，同時劃分的交接處要滿足連續(xù)性條件，所述若干部分包括與壓頭接觸的接觸區(qū)域及與壓頭不接觸的非接觸區(qū)域；

3)根據(jù)邊界條件和連續(xù)性條件使用“d-s計算”求解彈性體劃分的各個部分的位移函數(shù)和應(yīng)力函數(shù)；

4)將若干部分彈性體頂部的位移函數(shù)合并成一個整體的位移函數(shù)；

5)將步驟4)中合并得到位移函數(shù)作為整個彈性體的邊界條件進(jìn)行“d-s計算”，得到整個彈性體的位移函數(shù)和應(yīng)力函數(shù)；

6)根據(jù)步驟5)的應(yīng)力函數(shù)判斷非接觸區(qū)域的頂部正應(yīng)力和切應(yīng)力是否接近于0；若是，執(zhí)行步驟10；若不是，執(zhí)行步驟7；

7)將非接觸區(qū)域頂部應(yīng)力分別取反后作為各自的邊界條件進(jìn)行“d-s計算”，得到各個非接觸彈性體的位移函數(shù)和應(yīng)力函數(shù)；

8)將步驟5)與步驟7)非接觸區(qū)域頂部的位移函數(shù)疊加作為非接觸區(qū)域頂部的位移，并計算非接觸區(qū)域與接觸區(qū)域交接點(diǎn)的位移，將該交接點(diǎn)的位移添加到壓頭的外廓函數(shù)作為接觸區(qū)域的位移；

9)回到步驟4)執(zhí)行算法；

10)根據(jù)位移函數(shù)和應(yīng)力函數(shù)求需要的量，如壓痕深度、最大應(yīng)力等，也可通過積分求得作用力。

在圖1的流程示意圖和上述步驟中提到的“d-s計算”如圖2所示。為了便于說明這個計算方法，以下做些基本說明和公式推導(dǎo)。

如圖3所示，這是一個楔形凸面與彈性體接觸模型，研究的彈性體是正交各向異性材料且彈性體底面固定，同時不考慮摩擦力的因素且楔形的斜度較小。以紙面的橫向?yàn)閤1軸方向，紙面的縱向?yàn)閤2軸方向，垂直于紙面向外方向?yàn)閤3軸方向，并假設(shè)彈性體在x3軸方向上足夠長，從而可將接觸問題視為x1-x2平面問題。將彈性體分成三個部分：part1(圖中i區(qū)域)、part2(圖中ii區(qū)域)和part3(圖中iii區(qū)域)，其中只有part2與楔形凸面發(fā)生接觸。其中，彈性體的寬度為l，高度為h，剛體壓頭與彈性體發(fā)生的接觸長度為2c，接觸高度為d，壓頭從剛接觸到最后的穩(wěn)定狀態(tài)在x2方向的移動距離為d，接觸間隙高度u0＝d-d。不考慮彈性體的自身重力，彈性體滿足以下條件：

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

σij＝cijksεks，(k，s＝1，2，3)(2)

其中，應(yīng)力σij＝σji，εks表示彈性體應(yīng)變，uk、us分別表示彈性體在xk、xs方向上的位移，cijks是彈性體的材料參數(shù)，并滿足全對稱性：

cijks＝cjiks＝cksij＝cijsk(4)

將式(3)代入式(2)后求導(dǎo)可得

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

根據(jù)eshelby-reid-shockley公式，設(shè)定ui＝aif(z)或者u＝af(z)，其中z＝x1+px2，則有

uk，s＝(δs1+pδs2)akf′(z)(6)將其代入到式(5)可得

cijks(δj1+pδj2)(δs1+pδs2)ak＝0

(7)

或{ci1k1+p(ci1k2+ci2k1)+p²ci2k2}ak＝0

可以寫成

{q+p(r+r^t)+p²t}a＝0(8)

其中，q、r、t都是3x3的矩陣：qik＝ci1k1，rik＝ci1k2，tik＝ci2k2。對于正交各向異件材料，則有：

可以看出在應(yīng)變能為正的時候，q和t是對稱矩陣且是正定的。想要式子中的a有非零解，必須要滿足

|q+p(r+r^t)+p²t|＝0(12)

定義

對于這個式子的求解，可以引入矩陣n使得

nζ＝pζ(14)

其中，

p是特征值，ζ是特征向量。p是復(fù)數(shù)，在6個特征值p中有三對是共軛復(fù)數(shù)。假定im(pα)＞0，是pα的共軛復(fù)數(shù)。便可獲得位移u的通解：

同樣地，應(yīng)力的通解可以寫成：

即

其中，(σ1)i＝σi1，(σ2)i＝σi2。

假定函數(shù)

其中，

則位移和應(yīng)力可表示成：

其中，a＝[a1a2a3]，b＝[b1b2b3]，<φ*ψ*χ*>＝diag[φ1ψ1χ1，φ2ψ2χ2，φ3ψ3χ3]，

至此，三個彈性體的應(yīng)力和位移可以表示為：

在有了以上的基礎(chǔ)后，便可以計算彈性體分3部分時各自的位移和應(yīng)力。將彈性體分成3部分之后的邊界條件如圖4所示，由于將接觸問題視為平面問題，故始終存在

u3＝0(25)

又由于part1和part3的頂部未與壓頭發(fā)生接觸，則這兩部分滿足

σ12＝0，σ22＝0(26)

彈性體part2頂部與壓頭發(fā)生接觸，忽略摩擦力的因素，則有

σ12＝0，u2＝u2(x1)＝g(x1)+u0(27)

其中，g(x1)是壓頭的外廓函數(shù)。在本發(fā)明中是以凸型斜面為例，若只針對part2分析，外廓函數(shù)為

若針對整個彈性體分析，外廓函數(shù)為

part1、part2和part3的底面是固定的，故滿足邊界條件：

u1＝0，u2＝0(30)

part1左側(cè)和part3右側(cè)也不發(fā)生接觸，故滿足邊界條件：

σ11＝0，σ21＝0(31)

以上的邊界條件總結(jié)為

part1右側(cè)與part2左側(cè)相接，part2右側(cè)與part3左側(cè)相接，因此在這要滿足位移和應(yīng)力的連續(xù)性，固有約束條件：

轉(zhuǎn)換到各個部分，可表示為

將邊界條件和連續(xù)性條件的式子進(jìn)行傅里葉積分，得到

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

若將u0作為一個未知量來考慮，需要通過求導(dǎo)來消去u0，part2的頂部的邊界條件式子轉(zhuǎn)化為

其中，

根據(jù)之前的k和m值選取合適的j值(保證線性方程組35～40有唯一解)求解出式中的各個參數(shù)，便可得到三部分彈性體的位移、應(yīng)力表達(dá)式。從“設(shè)定”到求解出彈性體的位移應(yīng)力的這一過程，在本發(fā)明中稱為“d-s計算”，流程如圖2所示，上述的計算是將彈性體劃分為三個區(qū)域的d-s計算。

由于在實(shí)際計算過程中，k和m值不可能無限大，求解得到的結(jié)果往往位移收斂性較好，而應(yīng)力并沒有很好地收斂。在part1與part2交接處和part2與part3交接處存在著應(yīng)力突變以及應(yīng)力不近似于零等問題，為了解決這些問題，依據(jù)線性疊加原理使用迭代法得到更好的函數(shù)。

對于彈性體的頂部，在x2方向上，part1和part3給定的是位移約束σ22＝0，而part2給定的是位移約束u2＝u2(x1)＝g(x1)+u0。為了將三個部分的彈性體重新作為一個彈性體考慮，彈性體頂部需要給定一樣的位移約束，具體的邊界條件如圖5所示。其中，

假設(shè)整體分析時

其中，zα＝x1+pαx2，

與三部分相同，得到整體的位移、應(yīng)變表達(dá)式：

u(0～l，0～h)、σ1(0～l，0～h)、σ2(0～l，0～h)。

邊界條件為

其中，

對這四個邊界條件式子進(jìn)行傅里葉積分，得到

其中，g^t(x1)＝[0，u2(x1)，0]^t，g^b(x1)＝g^l(x2)＝g^r(x2)＝[0，0，0]^t。

同樣根據(jù)k和m的值設(shè)定合適的j值來求解方程組，求得fα(zα)各項(xiàng)參數(shù)，便可得到整個彈性體位移、應(yīng)力的函數(shù)表達(dá)式在交接處的位移和應(yīng)力的連續(xù)性得到了很大改善，但仍會存在非接觸區(qū)域應(yīng)力不近似于零的現(xiàn)象。這時，需要提取出part1和part3進(jìn)行處理，將整體計算得到的σ12、σ22反向后最為part1和part3的應(yīng)力約束，邊界條件如圖6和圖7所示，特別需要注意的是σ12＝-σ12和σ22＝-σ22是賦值語句，是將整體計算得到的結(jié)果取反后作為邊界條件。

其中，n＝1，3，和表示整體計算時得到的應(yīng)力。對上述式子進(jìn)行傅里葉積分得到

其中，同樣根據(jù)k和m的值設(shè)定合適的j值來求解方程組，求得各項(xiàng)參數(shù)，便可得到圖6和圖7邊界條件下位移、應(yīng)力的函數(shù)表達(dá)式

然后將part1和part3的計算結(jié)果與之前整體計算的結(jié)果進(jìn)行疊加，疊加之后的邊界約束如圖8所示。

然后再進(jìn)行一次整體計算，邊界條件與之前整體計算的邊界條件相同，如圖5所示，但頂面的位移約束條件發(fā)生了改變：

完成本次整體計算之后判斷非接觸區(qū)域(part1和part3)頂部各個部位是否滿足σ12≈0、σ22≈0。由于σ2是傅里葉級數(shù)的形式，且k和m的取值有限，很難保證頂部的每個點(diǎn)都滿足σ12＝0、σ22＝0，因此需要設(shè)定合適的閾值判斷其是否近似為0。通常情況下，在part1、part3與part2相接處附近的各點(diǎn)(不包括part2上的點(diǎn))滿足σ12≈0和σ22≈0時，便可認(rèn)為得到的fα(zα)是滿足我們最初邊界條件的函數(shù)。假如不滿足這兩個條件，需要進(jìn)行迭代運(yùn)算，重復(fù)“分別提取part1和part3并給予反向應(yīng)力計算”、“單獨(dú)計算結(jié)果疊加到整體計算結(jié)果后作為頂部位移約束進(jìn)行整體計算”。

在得到滿足要求的fα(zα)函數(shù)后，將fα(zα)代入式(16)、式(19)和式(20)，便得到了彈性體各點(diǎn)的位移函數(shù)和應(yīng)力函數(shù)。再對part2頂部的σ22進(jìn)行積分便能得到作用力p：

若考慮壓頭的外廓，則有

其中，θ是接觸點(diǎn)切線與坐標(biāo)軸x1的夾角。