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

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

2、邊界的分布特點，將區(qū)域分解成不同的子區(qū)域，引進拉格朗日乘子表示接觸力，保證了各子區(qū)域的網(wǎng)格剖分和位移求解是完全獨立的?；谏鲜鏊惴ê陀邢拊绦蜃詣由上到y(tǒng)開發(fā)了相應的求解帶摩擦接觸問題的軟件，數(shù)值實驗表明，程序是正確的，算法是高效的。（字體楷體_GB2312、五號字）關(guān)鍵詞：摩擦，接觸問題，擬高斯迭代法，拉格朗日乘子，區(qū)域分解（請選出4至5個關(guān)鍵字，如：有限元、接觸、FEPG、并行）Abstract（小四號字體加粗） 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、五號字）Key words: （小四號字體加粗） impact，thin-walled tube，dynamic progressive buckling，non-linear FEM（字體楷體_GB2312、五號字）（注意：本模板僅作為對論文字體、段落、圖片、參考文獻和標題等的示范。除中英文的摘要和關(guān)鍵字必須編寫之外，論文文章結(jié)構(gòu)由作者自行決定。本文的文章結(jié)

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

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

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

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

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

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

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

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

15、0.6。計算時間步長是1 秒，一共計算60步。柱體四個側(cè)面x，y方向位移都被固定，接觸面不會滑動，這就相當于一個連續(xù)體，可以與上面的理論解作比較。我們主要比較塊體上底面中心點z 向位移的數(shù)值解和理論解。 取 該問題的理論解18為 (5.2)在不同時間步，理論解和數(shù)值解的比較見表 2和圖5，表2 不同時間步的理論解和數(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)達到平衡，平衡后的z方向的位移以及x，y，z方向的應力分別如圖6 (a)，(b)，(c)，(d)所示。 圖6(a) 圖6 (b) 圖6 (c) 圖6(d)從圖5可看出有限元解能同理論解很好的吻合；從圖6可看出應力松弛平衡后，三個方向的應力都收斂于壓力 。由上述結(jié)果可知，程序是正確的和可靠的。6、結(jié)論本文提出的求解帶摩擦的接觸問題的數(shù)值算法擬高斯迭代法具有以下特點：

17、對法向接觸力和切向接觸力進行交替迭代，并利用高斯迭代法求解法向接觸力，利用分塊高斯迭代法求解切向接觸力；同其它數(shù)值算法相比，該算法保留了關(guān)于接觸力的柔度矩陣的稀疏性和對稱性，利用矩陣乘向量可以分步進行的技巧，該算法只存儲關(guān)于接觸力的柔度矩陣的下三角形矩陣的非零元和對角矩陣；根據(jù)可能接觸邊界的分布特點，將區(qū)域分解成不同的子區(qū)域，引進拉格朗日乘子表示接觸力，保證了各子區(qū)域的網(wǎng)格剖分和位移求解是完全獨立的，大大提高了該算法求解大規(guī)模接觸問題的能力；該算法對法向接觸力和切向接觸力進行交替迭代，關(guān)于該迭代法的收斂性的嚴格的數(shù)學證明有待進一步研究。參考文獻（宋體、四號、加粗、居中）1. 陳國慶、陳萬吉、馮恩明，三維接觸問題的非線性互補原理及算法，中國科學，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