求解帶摩擦接觸問題的擬高斯迭代法

1、求解帶摩擦接觸問題的擬高斯迭代法（標(biāo)題2、三號(hào)、黑體、加粗）金朝 作者2 （宋體、五號(hào)字） 作 者 （宋體、五號(hào)字）中國鐵道科學(xué)院 （宋體、五號(hào)字） 單 位2 （宋體、五號(hào)字）（正文字體為宋體、五號(hào)字，行間距為18磅）摘 要：（小四號(hào)字體加粗）本文提出了一種新的求解帶摩擦的接觸問題的數(shù)值算法，即擬高斯迭代法。它對(duì)法向接觸力和切向接觸力進(jìn)行交替迭代，并利用高斯迭代法求解法向接觸力，利用分塊高斯迭代法求解切向接觸力。同其它的數(shù)值算法相比，該算法保留了關(guān)于接觸力的柔度矩陣的稀疏性和對(duì)稱性，利用矩陣乘向量可以分步進(jìn)行的技巧，該算法只存儲(chǔ)關(guān)于接觸力的柔度矩陣的下三角形矩陣的非零元和對(duì)角矩陣。根據(jù)可能接觸

2、邊界的分布特點(diǎn)，將區(qū)域分解成不同的子區(qū)域，引進(jìn)拉格朗日乘子表示接觸力，保證了各子區(qū)域的網(wǎng)格剖分和位移求解是完全獨(dú)立的。基于上述算法和有限元程序自動(dòng)生成系統(tǒng)開發(fā)了相應(yīng)的求解帶摩擦接觸問題的軟件，數(shù)值實(shí)驗(yàn)表明，程序是正確的，算法是高效的。（字體楷體_GB2312、五號(hào)字）關(guān)鍵詞：摩擦，接觸問題，擬高斯迭代法，拉格朗日乘子，區(qū)域分解（請(qǐng)選出4至5個(gè)關(guān)鍵字，如：有限元、接觸、FEPG、并行）Abstract（小四號(hào)字體加粗） A numerical code RFPA3D based on linear elastic damage mechanics on mesoscopic scale and

3、FEM are developed to simulate the complete progressive 3-D failure and macroscopic mechanical behaviours of rock-like materials. Rock properties such as elastic constants, peak strength, and Poisson ratio are randomly distributed to reflect the initial random distributed weakness in mesoscopic scale

4、. By introducing heterogeneity of rock properties into the model, the RFPA3D code can simulate non-linear deformability of a quasi-brittle behaviour with an ideal brittle constitutive law. Each element in rock samples is considered to be homogenous and obeys a simple constituent law. A self-adaptive

5、 loading method is adopted to control the number of elements as less as possible that reach the failure threshold. A specimen was subjected to uniaxial tension. Complete stress-strain curves are obtained. Different from 2-D failure, 3-D samples show more ductile features. Crack initiation, propagati

6、on, coalescence and stress distribution can all be observed obviously during progressive failure process. （字體楷體_GB2312、五號(hào)字）Key words: （小四號(hào)字體加粗） impact，thin-walled tube，dynamic progressive buckling，non-linear FEM（字體楷體_GB2312、五號(hào)字）（注意：本模板僅作為對(duì)論文字體、段落、圖片、參考文獻(xiàn)和標(biāo)題等的示范。除中英文的摘要和關(guān)鍵字必須編寫之外，論文文章結(jié)構(gòu)由作者自行決定。本文的文章結(jié)

7、構(gòu)不屬于我們對(duì)論文編寫的規(guī)定之一）1、引言（大標(biāo)題頂格、標(biāo)題三、四號(hào)字、加粗）接觸問題在工程和自然界中大量存在，如高速火車的弓網(wǎng)接觸問題，帶斷層或裂縫的拱壩穩(wěn)定性問題，地震震源破裂和滑坡問題等。其特點(diǎn)和難點(diǎn)是接觸邊界和接觸力事先是未知的，當(dāng)考慮摩擦滑動(dòng)時(shí)，能量變分原理不再成立1。目前求解接觸問題的方法有解析法23，半解析法45，和數(shù)值算法，如實(shí)驗(yàn)誤差法67，線性互補(bǔ)法8，非線性互補(bǔ)法1，罰函數(shù)法910和拉格朗日乘子法11。隨著計(jì)算機(jī)的發(fā)展，基于有限元的數(shù)值算法成了求解接觸問題的主要方法。最早采用而且最常用的數(shù)值算法是實(shí)驗(yàn)誤差法，即對(duì)狀態(tài)進(jìn)行開閉迭代67，在給定狀態(tài)下，需要求解的是關(guān)于位移和接觸

8、力的等式方程。其優(yōu)點(diǎn)是力學(xué)意義清晰，程序易實(shí)現(xiàn)。其缺點(diǎn)是，當(dāng)考慮滑動(dòng)時(shí)，關(guān)于位移和接觸力系數(shù)矩陣是非對(duì)稱不定的。當(dāng)接觸邊界較大，接觸狀態(tài)變化較快時(shí)，已有的實(shí)驗(yàn)誤差法一般是不收斂的。本文提出了一種新的求解帶摩擦接觸問題的數(shù)值算法，即擬高斯迭代法。它既具有實(shí)驗(yàn)誤差法的優(yōu)點(diǎn)，又避免了實(shí)驗(yàn)誤差法的缺點(diǎn)。它對(duì)法向接觸力和切向接觸力進(jìn)行交替迭代，并利用高斯迭代法求解法向接觸力，利用分塊高斯迭代法求解切向接觸力。同其它數(shù)值算法相比，該算法保留了關(guān)于接觸力的柔度矩陣的稀疏性和對(duì)稱性，利用矩陣乘向量可以分步進(jìn)行的技巧，該算法只存儲(chǔ)關(guān)于接觸力的柔度矩陣的下三角形矩陣的非零元和對(duì)角矩陣。從子結(jié)構(gòu)的基礎(chǔ)上發(fā)展起來的區(qū)

9、域分解算法12是求解大規(guī)模問題最有效的方法。本文的算法根據(jù)可能接觸邊界的分布特點(diǎn)，將區(qū)域分解成不同的子區(qū)域，引進(jìn)拉格朗日乘子表示接觸力，保證了各子區(qū)域的網(wǎng)格剖分和位移求解是完全獨(dú)立的。對(duì)于罰函數(shù)法910來說，在不引進(jìn)新的邊界變量的情況下，整個(gè)區(qū)域的位移場必須同時(shí)求解，而且關(guān)于整個(gè)區(qū)域的位移場的剛度矩陣是變化的。線性互補(bǔ)8和非線性互補(bǔ)法1在一定的條件下彌補(bǔ)了實(shí)驗(yàn)誤差法收斂性很難保證的缺陷，但是它們破壞了關(guān)于接觸力的柔度矩陣的稀疏性。當(dāng)接觸區(qū)域很大或接觸邊界的網(wǎng)格剖分很密時(shí)，關(guān)于接觸力的柔度矩陣的存儲(chǔ)量和計(jì)算量會(huì)大得讓人難以接受。本文的算法保留了關(guān)于接觸力的柔度矩陣的稀疏性和對(duì)稱性?；诒疚牡乃惴?/p>

10、和有限元程序自動(dòng)生成系統(tǒng)13開發(fā)了相應(yīng)的軟件，數(shù)值實(shí)驗(yàn)表明，程序是正確的，算法是高效的。2、接觸問題的基本方程及其離散格式根據(jù)可能接觸邊界的分布特點(diǎn)，將區(qū)域分解成不同的子區(qū)域?？紤]由 M 個(gè)子區(qū)域 組成的區(qū)域。對(duì)于子區(qū)域 ，其動(dòng)力學(xué)平衡方程可表示為 (2.1)其中 表示密度， 表示阻尼， 表示應(yīng)力張量， 表示位移， 表示速度， 表示加速度， 表示體力。子區(qū)域 的初邊值條件可表示為： 其中和分別表示位移邊界和應(yīng)力邊界。引進(jìn)拉格朗日乘子 表示接觸力，不同子區(qū)域之間的可能接觸邊界條件為：a)分離狀態(tài) (2.3a).b)粘著狀態(tài) (2.3b)c)滑動(dòng)狀態(tài) (2.3c)其中 分別表示法向和切向分量，表示

11、初始間隙。分別表示法向和切向接觸力， 為摩擦系數(shù)。3、擬高斯迭代法及其矩陣表示下面我們將利用擬高斯迭代法求解第N個(gè)時(shí)間步滿足關(guān)系式（2.7）(2.9)的位移增量與接觸力增量。為了敘述方便，這里我們將關(guān)于時(shí)間步N的下標(biāo)省略。將式（2.7）寫成整體的形式，有（3.1）其中其中L 表示拉格朗日乘子的個(gè)數(shù)，為與第 個(gè)拉格朗日乘子相關(guān)的子區(qū)域號(hào)。 5、數(shù)值實(shí)驗(yàn)5.1 算例1（小標(biāo)題頂格、小四、加粗） 將一根無限長圓柱放在非常光滑的平面上，并在它的頂端施加單位載荷 F ，如圖1所示。RFA （圖片/文本框的版式應(yīng)選為：“嵌入型”）做法：按鼠標(biāo)右鍵，選擇設(shè)置圖片格式。在版式頁中選擇“嵌入型”圖1 模型（圖注

12、應(yīng)在圖片正下方）這是一個(gè)經(jīng)典的 Herz接觸問題，最大接觸力發(fā)生在接觸點(diǎn) A 處，由 Goldsmith17 可知，最大接觸壓力 和接觸半徑 b 的理論解為 (5.1) 其中 彈模， 泊松比。這是一個(gè)平面應(yīng)變問題，考慮到對(duì)稱性，取圓柱截面四分之一進(jìn)行計(jì)算。為了進(jìn)行接觸判定，在圓柱截面下面加了一個(gè)不變形的剛性體。雖然進(jìn)行了網(wǎng)格剖分，但每個(gè)節(jié)點(diǎn)的自由度為固定約束。有限元模型及其網(wǎng)格如圖2所示。 取 A 點(diǎn)垂直方向的正應(yīng)力同最大接觸壓力進(jìn)行比較。對(duì)于不同的載荷 F 和不同的彈模 E ，接觸半徑和最大接觸力與理論值的比較見表 1，對(duì)應(yīng)于表1的最后一種工況，圓柱體受載后的變形情況如圖3所示。表1 接觸半

13、徑和最大接觸力與理論值的比較 （序號(hào)與文字中間空一格，整體在表格之上并居中） F E(彈模) 1600200027232672374536342284%16004000192519225289952060159%16006000157215886478864067111%8002000192519222645026030158%80040001361137237405370790872%800600011111146458124555104569% 圖 2有限元模型及其網(wǎng)格 圖 3：變形圖 由表1可知，接觸壓力的誤差在 之內(nèi)，說明計(jì)算結(jié)果是正確的；接觸壓力的誤差隨著接觸半徑變小而變小，這與He

14、rz 理論解要求接觸半徑遠(yuǎn)遠(yuǎn)小于圓柱半徑 ，即 相吻合。接觸半徑的誤差相對(duì)要大些，這是由于接觸半徑與網(wǎng)格的尺寸和接觸端點(diǎn)的位置有關(guān)。5.2 算例 2我們用只有垂直方向變形的三維粘彈性塊體來檢驗(yàn)程序及算法的正確性，模型及其網(wǎng)格如下圖所示接觸面 圖4：模型及其網(wǎng)格其中塊體長、寬均為30m，高為600m，中間內(nèi)部有一傾斜接觸面（沒有接觸到柱體側(cè)面），傾角為 。從t=0時(shí)刻起在其頂面施加 (pa) 的恒定壓力，下底面固定，四個(gè)側(cè)面只有方向的位移，方向位移約束。楊氏模量為1.0e10Pa，粘性系數(shù)為1.0e10Pas，泊松比取為0.25，塊體密度和阻尼系數(shù)均定為0（不考慮加速度影響），接觸面的摩擦系數(shù)為

15、0.6。計(jì)算時(shí)間步長是1 秒，一共計(jì)算60步。柱體四個(gè)側(cè)面x，y方向位移都被固定，接觸面不會(huì)滑動(dòng)，這就相當(dāng)于一個(gè)連續(xù)體，可以與上面的理論解作比較。我們主要比較塊體上底面中心點(diǎn)z 向位移的數(shù)值解和理論解。 取 該問題的理論解18為 (5.2)在不同時(shí)間步，理論解和數(shù)值解的比較見表 2和圖5，表2 不同時(shí)間步的理論解和數(shù)值解時(shí)間 t 數(shù)值解 理論解 誤差10-0.84623E-04-0.83498E-040.11251E-0620-0.89277E-04-0.89295E-040.18387E-0730-0.89903E-04-0.89924E-040.20645E-0740-0.89987E-0

16、4-0.89992E-040.47274E-0850-0.89998E-04-0.89999E-040.11042E-0860-0.90000E-04-0.90000E-04-0.97176E-010圖5 由圖5可看出，60 秒后系統(tǒng)達(dá)到平衡，平衡后的z方向的位移以及x，y，z方向的應(yīng)力分別如圖6 (a)，(b)，(c)，(d)所示。 圖6(a) 圖6 (b) 圖6 (c) 圖6(d)從圖5可看出有限元解能同理論解很好的吻合；從圖6可看出應(yīng)力松弛平衡后，三個(gè)方向的應(yīng)力都收斂于壓力 。由上述結(jié)果可知，程序是正確的和可靠的。6、結(jié)論本文提出的求解帶摩擦的接觸問題的數(shù)值算法擬高斯迭代法具有以下特點(diǎn)：

17、對(duì)法向接觸力和切向接觸力進(jìn)行交替迭代，并利用高斯迭代法求解法向接觸力，利用分塊高斯迭代法求解切向接觸力；同其它數(shù)值算法相比，該算法保留了關(guān)于接觸力的柔度矩陣的稀疏性和對(duì)稱性，利用矩陣乘向量可以分步進(jìn)行的技巧，該算法只存儲(chǔ)關(guān)于接觸力的柔度矩陣的下三角形矩陣的非零元和對(duì)角矩陣；根據(jù)可能接觸邊界的分布特點(diǎn)，將區(qū)域分解成不同的子區(qū)域，引進(jìn)拉格朗日乘子表示接觸力，保證了各子區(qū)域的網(wǎng)格剖分和位移求解是完全獨(dú)立的，大大提高了該算法求解大規(guī)模接觸問題的能力；該算法對(duì)法向接觸力和切向接觸力進(jìn)行交替迭代，關(guān)于該迭代法的收斂性的嚴(yán)格的數(shù)學(xué)證明有待進(jìn)一步研究。參考文獻(xiàn)（宋體、四號(hào)、加粗、居中）1. 陳國慶、陳萬吉、馮恩明，三維接觸問題的非線性互補(bǔ)原理及算法，中國科學(xué)，1995年11月，25（11）， P11811190。2. K.L Johnson , Contact Mechanics, Cambrige University Press (1995)3Sijhoff & Noordhoff , (1980)4. T.F Conrey and A.Seireg, A mathematical programming method for design of elastic bodies in