接觸問題的非線性有限元分析

1、第五章 接觸問題的非線性有限元分析5.1引言在工程結(jié)構(gòu)中，經(jīng)常會(huì)遇到大量的接觸問題。火車車輪與鋼軌之間，齒輪的嚙合是典型的接觸問題。在水利和土木工程中，建筑物基礎(chǔ)與地基，混凝土壩分縫兩側(cè)，地下洞室襯砌與圍巖之間，巖體結(jié)構(gòu)面兩側(cè)都存在接觸問題。對(duì)于具有接觸面的結(jié)構(gòu)，在承受荷載的過(guò)程中，接觸面的狀態(tài)通常是變化的，這將影響接觸體的應(yīng)力場(chǎng)。而應(yīng)力場(chǎng)的改變反過(guò)來(lái)又影響接觸狀態(tài)，這是一個(gè)非線性的過(guò)程。由于接觸問題對(duì)工程實(shí)踐的重要性，本章將作為專門問題進(jìn)行研究。最早對(duì)接觸問題進(jìn)行系統(tǒng)研究的是H. Hertz，他在1882年發(fā)表了彈性接觸問題一書中，提出經(jīng)典的Hertz彈性接觸理論。后來(lái)Boussinesg等

2、其他學(xué)者又進(jìn)一步發(fā)展了這個(gè)理論。但他們都是采用一些簡(jiǎn)單的數(shù)學(xué)公式來(lái)研究接觸問題，因而只能解決形狀簡(jiǎn)單（如半無(wú)限大體）、接觸狀態(tài)不復(fù)雜的接觸問題。二十世紀(jì)六十年代以后，隨著計(jì)算機(jī)和計(jì)算技術(shù)的發(fā)展，使應(yīng)用數(shù)值方法解決復(fù)雜接觸問題成為可能。目前，分析接觸問題的數(shù)值方法大致可分為三類：有限元法、邊界元法和數(shù)學(xué)規(guī)劃法。數(shù)學(xué)規(guī)劃法是一種優(yōu)化方法，求解接觸問題時(shí)，根據(jù)接觸準(zhǔn)則或變分不等式建立數(shù)學(xué)模型，然后采用二次規(guī)劃或罰函數(shù)方法給出解答。邊界元方法也被用來(lái)求解接觸問題，1980年和1981年，Anderson先后發(fā)表兩篇文章，用于求解無(wú)摩擦彈性接觸和有摩擦彈性接觸問題。近年來(lái)雖有所發(fā)展，但仍主要用于解決彈性

3、接觸問題。就目前的發(fā)展水平來(lái)看，數(shù)學(xué)規(guī)劃法和邊界元法只適合于解決比較簡(jiǎn)單的彈性接觸問題。對(duì)于相對(duì)復(fù)雜的接觸非線性問題，如大變形、彈塑性接觸問題，還是有限元方法比較成熟、比較有效。早在1970年，Wilson和Parsons提出一種位移有限元方法求解接觸問題。Chan和Tuba，Ohte等進(jìn)一步發(fā)展了這類方法。它的基本思想是假定接觸狀態(tài)，求出接觸力，檢驗(yàn)接觸條件，若與假定的接觸狀態(tài)不符，則重新假定接觸狀態(tài)，直至迭代計(jì)算得到的接觸狀態(tài)與假定狀態(tài)一致為止。具體做法是：對(duì)于彈性接觸的兩個(gè)物體，通過(guò)有限元離散，建立支配方程 (5.1)式中，為初始的整體勁度矩陣，它與接觸狀態(tài)有關(guān)，通常根據(jù)經(jīng)驗(yàn)和實(shí)際情況假

4、定。是結(jié)點(diǎn)位移列陣，為結(jié)點(diǎn)荷載列陣。求解式（5.1），得到結(jié)點(diǎn)位移，再計(jì)算接觸點(diǎn)的接觸力，將和代入與假定接觸狀態(tài)相應(yīng)的接觸條件，如果不滿足接觸條件，就要修改接觸狀態(tài)。根據(jù)修改后新的接觸狀態(tài)，建立新的勁度矩陣和支配方程 (5.2)再由式（5.2）解得，進(jìn)一步計(jì)算接觸力，將和代入接觸條件，驗(yàn)算接觸條件是否滿足。這樣不斷的迭代循環(huán)，直至和滿足接觸條件為止，此時(shí)得到的解答就是真實(shí)接觸狀態(tài)下的解答。在以上的研究中，沒有考慮接觸面的摩擦力。不考慮摩擦力的接觸過(guò)程是一種可逆的過(guò)程，即最終結(jié)果與加載途徑無(wú)關(guān)。此時(shí)，只需要進(jìn)行一次加載，就能得到最終穩(wěn)定的解。如果考慮接觸面的摩擦力，接觸過(guò)程就是不可逆的，必須采用

5、增量加載的方法進(jìn)行接觸分析。1973年，Tusta和Yamaji的文章詳細(xì)討論了接觸過(guò)程的可逆性和不可逆性。從Wilson和Parsons的方法可看出，每一次接觸狀態(tài)的改變，都要重新形成整體勁度矩陣，求解全部的支配方程，既占內(nèi)存，又費(fèi)機(jī)時(shí)。實(shí)際上，接觸狀態(tài)的改變是局部的，只有與接觸區(qū)域有關(guān)的一小部分需要變動(dòng)，為此又提出一些改進(jìn)的方法。1975年，F(xiàn)rancavilla和Zienkiewicz 提出相對(duì)簡(jiǎn)單的柔度法。圖5.1示出兩個(gè)相互接觸的物體A和B，假定A上有外力R作用，B有固定邊界。接觸面作用在A上的接觸力是，作用在B上的接觸力是，對(duì)于二維問題， (5.3)這些接觸力是未知的，假定有m個(gè)接

6、觸點(diǎn)對(duì)，則增加了4m個(gè)未知量，為此需要補(bǔ)充4m個(gè)方程?，F(xiàn)列出接觸點(diǎn)的柔度方程 (5.4)其中，和分別是物體A和B在接觸點(diǎn)i處的位移，和分別表示物體A和B因j點(diǎn)作用單位力時(shí)在i點(diǎn)引起的位移（即柔度系數(shù)）所組成的柔度子矩陣，m1是外荷載作用的點(diǎn)數(shù)，為第k個(gè)荷載作用點(diǎn)上的荷載向量。如果物體A和B之間的接觸屬于連續(xù)接觸，則接觸條件為 (5.5) (5.6)（5.5）和（5.6）是4m個(gè)補(bǔ)充方程，式中，是第i個(gè)接觸點(diǎn)對(duì)的初始間隙向量。由于式（5.6）的存在，令，未知量數(shù)目減少，增加的未知量剩下2m個(gè)。將式（5.4）和（5.6）代入（5.5）得 (5.7)式（5.7）共有2m個(gè)補(bǔ)充方程。對(duì)于滑動(dòng)接觸和不接

7、觸的自由邊界，同樣可根據(jù)相應(yīng)的接觸條件列出與式（5.7）類似的補(bǔ)充方程求解。引入接觸條件后，接觸狀態(tài)變化時(shí)，計(jì)算對(duì)象的整體勁度矩陣不再改變，出現(xiàn)的問題是增加了未知量數(shù)，需要建立補(bǔ)充方程。但由于補(bǔ)充方程（5.7）中，、和不隨接觸狀態(tài)的改變而變化，而且接觸點(diǎn)的數(shù)目遠(yuǎn)小于整體的結(jié)點(diǎn)數(shù)，因而可大大節(jié)約計(jì)算時(shí)間，提高了求解接觸問題的效率。另外一種提高接觸問題計(jì)算效率的方法是把接觸點(diǎn)對(duì)作為“單元”考慮。1979年，Okamoto和Nakazawa提出“接觸單元”，它是根據(jù)接觸點(diǎn)對(duì)位移與力之間的接觸條件建立的。接觸單元和普通單元一樣，可以直接組裝到整體勁度矩陣中去。然后對(duì)支配方程進(jìn)行“靜力凝聚”，保留接觸面

8、各點(diǎn)的自由度，得到在接觸點(diǎn)凝聚的支配方程。由于接觸點(diǎn)數(shù)遠(yuǎn)小于結(jié)點(diǎn)數(shù)，凝聚后的方程階數(shù)比未凝聚時(shí)方程階數(shù)低得多。當(dāng)接觸狀態(tài)改變時(shí)，只需對(duì)凝聚的支配方程進(jìn)行修正和求解，因而可節(jié)約計(jì)算時(shí)間。1975年，Schafer根據(jù)虛功原理推導(dǎo)了“連接單元”，也可以象普通單元一樣地形成和組裝到整體勁度矩陣中。連接單元包含有接觸面的接觸特性，通過(guò)改變形成單元的某些參數(shù)，來(lái)反映不同的接觸狀態(tài)。1977年，J. T. Stadter和R. O. Weiss提出間隙元方法。“間隙元”是一種虛設(shè)的具有一定物理性質(zhì)的特殊接觸單元，其內(nèi)部的應(yīng)力應(yīng)變反映了接觸狀態(tài)，并利用塑性力學(xué)中的“應(yīng)力不變”準(zhǔn)則來(lái)模擬接觸過(guò)程。目前的接觸研

9、究主要集中在彈性接觸問題，關(guān)于彈塑性接觸問題的研究也有了相當(dāng)?shù)倪M(jìn)展，但有關(guān)大變形彈塑性接觸的研究成果還很少。5.2 彈性接觸問題基本假定在分析彈性接觸問題時(shí)，有如下的基本假定：(1) 接觸物體的材料是線彈性的，位移和變形是微小的；(2) 作用在接觸面上的摩擦力服從Mohr-Coulomb準(zhǔn)則；(3) 接觸面連續(xù)平滑。所謂接觸條件，是指接觸面上接觸點(diǎn)處的位移和力的條件。利用接觸條件，可以判斷接觸物體之間的接觸狀態(tài)。接觸狀態(tài)可分為三類：連續(xù)接觸，滑動(dòng)接觸和自由邊界。為了更方便地表示接觸條件，需要在接觸面上建立局部坐標(biāo)系，如圖5.2所示。由于一般情況下，A、B兩個(gè)物體在接觸點(diǎn)處無(wú)公共切面和公共法線，

10、因此，局部坐標(biāo)系的軸只能盡可能地接近公法線方向，平面盡可能地接近公切面。令和分別是第j個(gè)接觸物體（j=A, B）沿第i個(gè)局部坐標(biāo)（i=）的位移和接觸力，則三類接觸條件可表示為：(1) 連續(xù)接觸條件 （i=） (5.8) （i=） (5.9)同時(shí)要滿足沿接觸面的切平面方向不滑動(dòng)的條件：和 (5.10)以上式中，是接觸面在方向的初始間隙，f是接觸面之間的滑動(dòng)摩擦系數(shù)。(2) 滑動(dòng)接觸條件 (5.11) （i=）或者表示為和 (5.12)其中，(3) 自由邊界條件 （i=） (5.13) (5.14)以上接觸條件中出現(xiàn)的位移和接觸力通常都是未知量，因此需要采用迭代算法，即首先假定接觸狀態(tài)，根據(jù)假定的

11、接觸狀態(tài)建立有限元求解的支配方程，求解方程得到接觸面的位移和接觸力，并校核接觸條件是否與原來(lái)假定的接觸狀態(tài)相符。若不同，就要修正接觸狀態(tài)，這樣不斷地循環(huán)，直到接觸狀態(tài)穩(wěn)定為止。實(shí)際上，這是一個(gè)局部的幾何非線性問題。對(duì)于接觸問題，存在可逆和不可逆兩種接觸狀態(tài)。所謂“可逆”，是指沿不同的加載途徑，其最終的結(jié)果是相同的。“不可逆”則是指對(duì)于不同的加載途徑，最后的結(jié)果不同。發(fā)生不可逆過(guò)程的原因是由于接觸面出現(xiàn)了滑動(dòng)摩擦，下面的例子可以說(shuō)明這一點(diǎn)。qt1t2t3-fP0fP0Q0t4QQPst0圖5.3示出一個(gè)由A、B、C三個(gè)物體組成的接觸問題。物體A上面作用有勻布荷載R，左邊為鉸支座。物體C左邊和下面

12、均有鉸支座。物體B在勻布荷載Q的作用下，可以沿著上下兩個(gè)接觸面滑移。加載分三步：(1) 施加荷載R，并保持R0不變；(2) 施加荷載Q，從0增加到Q0；(3) 逐漸減小荷載R，回到0?，F(xiàn)考察A、B接觸面上某一點(diǎn)s的切向力q隨荷載Q的變化過(guò)程。當(dāng)Q0時(shí)，假定q0。隨著Q的增大，q也增大，直到物體B發(fā)生滑動(dòng)，此時(shí)qfP0，見圖5.3中的t1點(diǎn)。q保持此值，直至QQ0為止，此時(shí)，相應(yīng)的點(diǎn)為t2。接著，Q開始減小，從平衡的角度，q也減小，逐漸到零。由于這時(shí)Q還未減小到零，q會(huì)繼續(xù)減小，實(shí)際上是改變符號(hào)，向相反方向增加，直到qfP0（t3點(diǎn)），當(dāng)Q減小到零時(shí)，回到點(diǎn)t4?？梢钥闯?，由于接觸面滑動(dòng)摩擦的存

13、在，最終狀態(tài)t4與初始狀態(tài)t0是不同的，說(shuō)明切向接觸力是不可逆的。因此，凡是考慮接觸面切向摩擦力的接觸問題，都應(yīng)當(dāng)按復(fù)雜加載過(guò)程來(lái)研究，即通過(guò)增量的方式求解。對(duì)于不考慮摩擦的可逆過(guò)程，是一種簡(jiǎn)單加載過(guò)程，可以一步加載完成求解。5.3 彈性接觸問題有限元基本方程和柔度法求解假設(shè)A、B是相互接觸的兩個(gè)物體，為了研究的方便，將它們分開，代之以接觸力PA和PB，如圖5.4所示。然后建立各自的有限元支配方程： (5.15)式中，KA、和RA分別是物體A的整體勁度矩陣、結(jié)點(diǎn)位移列陣和外荷載，KB、和RB分別是物體B的整體勁度矩陣、結(jié)點(diǎn)位移列陣和外荷載。顯然，接觸力PA和PB都是增加的未知量，無(wú)法由式（5.

14、15）求出，必須根據(jù)接觸面上接觸點(diǎn)對(duì)的相容條件確定。設(shè)A、B上的接觸點(diǎn)對(duì)為iA和iB（i=1, 2, , m），假定勁度矩陣KA和KB非奇異，可求逆，則由式（5.15）得到接觸點(diǎn)的柔度方程 (5.16)式中，i、j=1, 2, , m表示結(jié)點(diǎn)號(hào)，m是接觸點(diǎn)對(duì)數(shù)目，nA、nB分別為作用在物體A和B上外荷載的作用點(diǎn)數(shù)，和表示物體A和B上接觸點(diǎn)i的位移、是A和B上接觸點(diǎn)j的接觸力、為A和B上結(jié)點(diǎn)k的外荷載、表示物體A和B上，由j點(diǎn)的單位力引起的i點(diǎn)在x、y、z三個(gè)方向的位移，是一個(gè)33階的柔度矩陣。在列出相容條件，求解接觸問題之前，有兩個(gè)問題需要解決。首先是消除剛體位移的問題。因?yàn)榈玫椒匠蹋?.16

15、）的前提是KA和KB非奇異可求逆，也就是說(shuō)物體A和B要有足夠的約束，不會(huì)發(fā)生剛體位移。但是有些接觸物問題中，可能會(huì)有某個(gè)物體由于約束不夠產(chǎn)生剛體位移，此時(shí)須對(duì)剛體位移進(jìn)行處理。以圖5.4中的物體A為例，假定它的約束不夠，則KA為奇異矩陣，記為。引入虛擬的約束，消除A的剛體位移，則（5.15）的第一式可改寫為 (5.17)其中，是與虛擬約束相應(yīng)的位移向量，I是單位矩陣。由上式得到 (5.18)從式（5.18）導(dǎo)出物體A上接觸點(diǎn)的柔度方程 （i=1, 2, , m） (5.19)Fi是與剛體位移相應(yīng)的柔度矩陣。第二個(gè)問題是，要將上述整體坐標(biāo)系下的量轉(zhuǎn)化到接觸面的局部坐標(biāo)系。接觸點(diǎn)位移和接觸力在不同

16、坐標(biāo)系下的表達(dá)式有以下的關(guān)系 （A、B） (5.20)式中，是結(jié)點(diǎn)i的坐標(biāo)轉(zhuǎn)換矩陣，、分別是接觸面局部坐標(biāo)系下，結(jié)點(diǎn)i的接觸力和位移。將式（5.20）代入式（5.19），得 (5.21)其中，。同樣，將式（5.20）代入式（5.16）的第二式，得 (5.22)以下將針對(duì)三類接觸條件建立相應(yīng)的相容方程。(1) 連續(xù)邊界根據(jù)前面的連續(xù)邊界條件（5.9），可以建立接觸點(diǎn)的位移相容方程 (5.23)是第i個(gè)接觸點(diǎn)對(duì)在局部坐標(biāo)系下的初始間隙。將（5.21）和（5.22）代入（5.23），并注意有，可得 (5.24)式中， (5.25) (5.26)(2) 滑動(dòng)邊界接觸面局部坐標(biāo)系方向的位移仍然滿足式（5

17、.23），但在切平面的和方向，接觸力的合力已經(jīng)達(dá)到摩擦極限，按照Mohr-Coulomb定律，則有 (5.27)(3) 自由邊界 (5.28)以上建立的相容方程，為原來(lái)的有限元支配方程增加了3m個(gè)補(bǔ)充方程，以求解3m個(gè)增加的未知接觸力Pj (j=1, 2, , m)。在建立相容方程時(shí)，必須知道接觸狀態(tài)，而接觸狀態(tài)事先也是未知的，因此這是一個(gè)迭代求解的過(guò)程。一般先假定為連續(xù)接觸狀態(tài)，按式（5.24）建立全部接觸點(diǎn)的相容方程，求出接觸力后，驗(yàn)證接觸條件是否滿足連續(xù)接觸，若是則不作修改；若為滑動(dòng)狀態(tài)，就用式（5.27）來(lái)代替這個(gè)接觸點(diǎn)在和兩個(gè)方向相應(yīng)的方程；若是自由狀態(tài)，就用式（5.28）替換這個(gè)接觸點(diǎn)的所有相應(yīng)方程。這樣通過(guò)反復(fù)迭代，就可以求得真正的接觸力和相應(yīng)的相容方程。對(duì)于具有滑動(dòng)摩擦的接觸問題，由于接觸過(guò)程的不可逆，需要采用增量方式加載、假定分級(jí)加載的次數(shù)為np，在進(jìn)行第l級(jí)加載前已經(jīng)施加的混雜為和，本級(jí)荷載增量為和，這樣式（5.24）就變成 (5.29)注意式中的各項(xiàng)有，將上述各式代回式（5.29），得 (5.30)令 (5.31)則（5.30）成為 (5.32)式（5.32）為連續(xù)接觸條件相容方程的增量形式。對(duì)于滑動(dòng)接觸條件，方向的相容方程與式（5.32）類似，和方向上相容方程的增量形式可表示為 (5.33)對(duì)于自由接觸條件，相容