含節(jié)理單元的三維P型自適應(yīng)有限元解法

1、含節(jié)理單元的三維P型自適應(yīng)有限元解法             摘要本文提出了含三維無厚度節(jié)理單元、等厚度節(jié)理單元和變厚度節(jié)理單元的p型自適應(yīng)有限元模型，給出了三維節(jié)理單元升階譜有限元法的解題步驟，通過具體算例，驗(yàn)證了p型升階譜有限元法在求解含三維節(jié)理單元的有限元問題時(shí)的可行性及優(yōu)越性。 關(guān)鍵詞節(jié)理單元 升階譜有限元 p型有限元   有限元法是求解微分方程數(shù)值解的一種重要方法，對于一個(gè)給定的問題，為改善其有限元解的精度，可以采用以下3種方法。(1)h型有限元

2、法1，這種方法通過減小單元尺寸來提高有限元解的精度。(2)p型有限元法，這種方法通過增加基底函數(shù)的階次來提高有限元解的精度。(3)hp型有限元法，這種方法是以上兩種方法的綜合，它既減小單元尺寸，又增加基底函數(shù)的階次。作者所在的研究小組從1995年開始研究水工結(jié)構(gòu)的h型彈粘塑性有限單元方法，目前已建立了實(shí)用的二維分析軟件體系4,5，并在三維分析方面取得了進(jìn)展。從1999年開始，在水工結(jié)構(gòu)的p型自適應(yīng)分析方面也有所突破，1999年，程昭等人針對水工結(jié)構(gòu)分析問題提出了三維升階譜有限元分析方法。2001年，陳勝宏等人進(jìn)一步提出了二維問題的p型自適應(yīng)分析策略，并將自適應(yīng)有限元方法歸類為全域升階方法、單元

3、升階方法和自由度升階方法等三類。之后，費(fèi)文平等人將p型自適應(yīng)有限元分析方法推廣到三維彈粘塑性領(lǐng)域。但是，以上有關(guān)p型有限元的研究成果中均未涉及到斷層、節(jié)理這一類特殊單元。   大壩壩基、壩肩和巖石高邊坡等部位總是存在斷層、節(jié)理和軟弱夾層等大規(guī)模的不連續(xù)面，且對結(jié)構(gòu)的變形和穩(wěn)定影響巨大，故在有限元分析中應(yīng)給予高度的重視。古德曼(Goodman)最初運(yùn)用有限元技術(shù)模擬巖體工程中的非線性不連續(xù)面問題，并提出了無厚度的節(jié)理元的概念10。隨后，朱伯芳于1979年提出了等厚度節(jié)理元模型，并將其與無厚度節(jié)理元模型形成統(tǒng)一的計(jì)算公式11，在此基礎(chǔ)上，王鴻儒等人提出了變厚度的節(jié)理單元的彈塑性

4、模型并將其應(yīng)用到工程實(shí)踐中12。目前，國內(nèi)外在有關(guān)p型自適應(yīng)有限元分析的研究中，尚未涉及到這類特殊單元的處理問題，從而使研究成果的工程應(yīng)用受到一定程度的限制。   對于常規(guī)塊體單元的三維p型有限元模型，作者在文獻(xiàn)中已有詳盡的論述，本文主要給出三維無厚度節(jié)理單元、等厚度節(jié)理單元和變厚度節(jié)理單元的p型有限元模型，并給了升階譜的計(jì)算格式。實(shí)例分析結(jié)果表明，用p型有限元法來求解含三維節(jié)理單元的有限元問題具有收斂速度快、計(jì)算精度高的優(yōu)點(diǎn)。1  三維p型無厚度節(jié)理單元和等厚度節(jié)理單元模型   對厚度很小和厚度變化不大的節(jié)理，可以分別采用無厚度的節(jié)理單元和等

5、厚度的節(jié)理單元進(jìn)行模擬，可取如圖1所示的六面體節(jié)理單元，進(jìn)行單元的網(wǎng)格劃分。1.1  形函數(shù)及位移函數(shù)  對節(jié)理單元上下面的相應(yīng)點(diǎn)、棱和面可取相同的點(diǎn)基函數(shù)、棱基函數(shù)和面基函數(shù)，對無厚度節(jié)理單元或等厚度節(jié)理單元，不存在連續(xù)上、下兩面的棱基函數(shù)和面基函數(shù)，也不存在體基函數(shù)?；缀瘮?shù)的具體形式如下13：   點(diǎn)基函數(shù)(p1)式中：，。   棱基函數(shù)(p2)   面基函數(shù)(p4)式中：而為Legender多項(xiàng)式。令位移函數(shù)為   同理可以寫出V下，V上，W下，W上及v，w的具體表達(dá)式。 &#

6、160; 將基函數(shù)Ni,pEi,pF統(tǒng)一記為i，位移差(uN,i+4-uNi)，(uE,i+4-uEi),(uF2-uF1)統(tǒng)一記為廣義結(jié)點(diǎn)位移差ui，設(shè)單元基底函數(shù)個(gè)數(shù)為fe(p)，上式簡化為，同理有,1.2  三維等厚度節(jié)理單元或無厚度節(jié)理單元升階過程  當(dāng)p=1時(shí)：N1-1I-2I-3I-4I1I2I3I4I,I為3×3階的單位陣。當(dāng)p=2時(shí)，N2可在N1的基礎(chǔ)上進(jìn)行擴(kuò)充，擴(kuò)充矩陣Np=2為Np=2=-5I-6I-7I-8I5I6I7I8I。依此類推，可得最終的形函數(shù)矩陣為1.3  坐標(biāo)插值及坐標(biāo)變換  對于節(jié)理上、下面坐標(biāo)的插值仍采用各

7、面上的四個(gè)節(jié)點(diǎn)進(jìn)行插值，即式中：(xi,yi,zi)i=18為節(jié)理間面體單元8個(gè)頂點(diǎn)的整體坐標(biāo)。   定義整體坐標(biāo)系的x軸朝北，y軸朝西，z軸朝上，定義等厚度的節(jié)理單元或無厚度的節(jié)理的局部坐標(biāo)系的z為中面的法線朝上方向，y指向節(jié)理面的傾向，x軸由右手法則確定，并設(shè)等厚度節(jié)理的傾角為，傾向?yàn)椤?#160;  三維等厚節(jié)理或無厚節(jié)理單元局部坐標(biāo)與整體坐標(biāo)的轉(zhuǎn)換矩陣為1.4  三維等厚度節(jié)理的單元?jiǎng)偠染仃囀街校簭椥跃仃?；單元?yīng)變矩陣BLB1/eLNp根據(jù)虛功原理，單元?jiǎng)偠染仃嚍?.5  三維無厚度節(jié)理的單元?jiǎng)偠染仃?#160; 當(dāng)?shù)群穸裙?jié)理單元的厚度

8、e0時(shí)，即形成無厚度的節(jié)理單元。此時(shí)，可假定單元內(nèi)應(yīng)力分量與位移差成正比，同理可得單元?jiǎng)偠染仃嚍槭街校簽閱卧獎(jiǎng)哦染仃嚒?#160;        2  三維p型變厚度的節(jié)理單元模型    當(dāng)節(jié)理單元的厚度變化較大時(shí)，應(yīng)將等厚度節(jié)理單元推廣得到變厚度的節(jié)理單元，如圖2所示建立局部坐標(biāo)系。2.1  形函數(shù)及位移函數(shù)  六面體變厚度的節(jié)理單元的基底函數(shù)，由點(diǎn)基函數(shù)、棱基函數(shù)、面基函數(shù)和體基函數(shù)組成，各基底函數(shù)的具體形式點(diǎn)基函數(shù)(p1)式中：i=(-1)i,i=(-1)i/2+0.5

9、，i=(-1)i/4+0.75。   棱基函數(shù)(p2)式中：2i=1-i/12+0.65，2i12(1-(-1)i-1/4)；2i=i-18，將1i，1i,1i用向量的形式表達(dá)為1=0000-1-111-11-11T，1-11-110000-1-11T，1-1-111-11-110000T。     面基函數(shù)(p4)式中：i,j2，i+j=p。   體基函數(shù)(p6)令位移函數(shù)   將基底函數(shù)Ni,pEi,pFi,pB統(tǒng)一記為i，位移uNi,uEi,uFi,uB統(tǒng)一記為廣義結(jié)點(diǎn)位移ui，設(shè)單元基底

10、函數(shù)的個(gè)數(shù)為fe(p)，上式簡化為，同理。2.2  三維變厚度節(jié)理單元升階過程  三維變厚度節(jié)理單元可按表1的方式進(jìn)行升階。2.3  坐標(biāo)插值與坐標(biāo)變換  對于單元任一點(diǎn)的坐標(biāo)，仍采用單元的8個(gè)實(shí)結(jié)點(diǎn)進(jìn)行插值，即，此處i=Ni(i=1,2 ,8),(xi,yi,zi)為立方體單元的8個(gè)頂點(diǎn)坐標(biāo)。參見圖2，取局部坐標(biāo)系的xoy平面與變厚度六面體單元的中面重合，z軸垂直于中面。根據(jù)變厚度節(jié)理單元8個(gè)實(shí)結(jié)點(diǎn)的整體坐標(biāo)，可按如下方法建立坐標(biāo)轉(zhuǎn)換矩陣及局部坐標(biāo)系。   先求出變厚度六面體單元中面法向向量nx,y,z，其中x(y2+y6-y1-y

11、5)(z3+z7-z1-z5)-(z2+z6-z1-z5）（y3+y7-y1-y5),y(z2+z6-z1-z5)(x3+x7-x1-x5)-(x2+x6-x1-x5)(z3+z7-z1-z5)，z(x2+x6-x1-x5)(y3+y7-y1-y5)-（y2+y6-y1-y5)(x3+x7-x1-x5)。   轉(zhuǎn)換矩陣L的形成，  令    向量l2 m2 n2T可根據(jù)右手法則，由向量l1 m1 n1T和l3 m3 n3T確定，即2.4  三維變厚度節(jié)理的單元?jiǎng)偠染仃?#160; 變厚度節(jié)理單元的單元應(yīng)變和應(yīng)力分別為

12、0;  變厚度節(jié)理單元的單元?jiǎng)偠染仃嚍?#160;  若用分塊矩陣的形式來表示單元?jiǎng)偠染仃嚕瑒t3  含節(jié)理單元的三維p型自適應(yīng)有限元算法3.1  控制方程與離散  根據(jù)實(shí)際情況分類進(jìn)行單元的網(wǎng)格剖分，形成三種特殊單元及常規(guī)塊體單元的單元?jiǎng)偠染仃?，并組合成總體剛度矩陣。形成荷載列向量，即   式中右端三項(xiàng)分別為體力、面力和集中力，并假設(shè)面荷載作用在1面上(余類推)，而Ex2,+y2,+z2,Ex2,+y2,+x2,Ex,x,+y,y,+z,z,可得到控制方程3.2  升階方法與誤差估計(jì)  常用的升階分析方法

13、有3種，即全域升階法、自由度升階法和單元升階法。此處選用單元升階法進(jìn)行升階計(jì)算，即對精度不滿足要求的單元進(jìn)行升階。在誤差估計(jì)分析中，取高階的應(yīng)力和應(yīng)變?yōu)閼?yīng)力和應(yīng)變精確值的“最佳估計(jì)”。定義誤差能范為總能范數(shù)為3.3  含節(jié)理單元的三維p型自適應(yīng)有限元分析過程  設(shè)相對控制誤差能范為et，對每一個(gè)單元進(jìn)行誤差估計(jì)，按下式計(jì)算單元的誤差參數(shù)   若i1，則將該單元進(jìn)行升階，誤差估計(jì)結(jié)束后，重新進(jìn)行彈粘塑性的有限元分析，如此反復(fù)，直到所有的單元均滿足精度要求。         4&

14、#160; 算 例4.1  算例1  為驗(yàn)證含節(jié)理單元的三維p型自適應(yīng)有限元理論及程序的可行性和優(yōu)越性，考察具有代表性的三維變厚度節(jié)理單元，取如圖3所示的計(jì)算試件，對程序進(jìn)行考核，單元及結(jié)點(diǎn)編號如圖3。其中，單元、為塊體單元，單元為變厚度節(jié)理單元，試件尺寸1m×1m×2m試件頂部作用一面力荷載-10.MPa，底部4結(jié)點(diǎn)固定結(jié)束，相對控制誤差取3.0%，材料參數(shù)見表2。同時(shí)，將該試件進(jìn)行網(wǎng)格細(xì)分，得到如圖4所示的計(jì)算試件細(xì)分網(wǎng)格模型，進(jìn)行常規(guī)有限元計(jì)算，并將p型有限元和各階計(jì)算結(jié)果與細(xì)網(wǎng)格的常規(guī)有限元計(jì)算結(jié)點(diǎn)位移值進(jìn)行比較，見表3。計(jì)算結(jié)果表明，在p=4時(shí)

15、，p型有限單元法達(dá)到要求的計(jì)算精度，隨著階次的逐步升高，p型有限元的計(jì)算結(jié)果與細(xì)網(wǎng)格的常規(guī)有限元計(jì)算結(jié)果越來越接近，在p=4時(shí)兩種有限單元法的計(jì)算結(jié)果相差甚微。這說明，p型有限元單元法在處理三維節(jié)理單元時(shí)是切實(shí)可行的。 同時(shí)可以看出，與常規(guī)有限元法相比，p型有限單元法對網(wǎng)格劃分的要求很低，用較少的單元數(shù)就能同等效果地模擬較復(fù)雜的常規(guī)有限單元網(wǎng)格，且具有很快的收斂速度。  4.2  算例2  考察一混凝土重力壩，壩基中含一個(gè)變厚度的節(jié)理面節(jié)理1和一個(gè)等厚度的節(jié)理面節(jié)理2，并將壩體與壩基的膠結(jié)面視為一個(gè)無厚度的節(jié)理面節(jié)理3。壩高100m，壩頂寬10m，上游坡

16、面垂直，下游坡面斜率10.75。壩基上游取1.5倍的壩高，下游取2倍的壩高，壩基的深度取2倍的壩高。再取10m的壩段寬構(gòu)成三維六面體網(wǎng)格，單元網(wǎng)格如圖5所示(x-z平面投影)，共160個(gè)單元，384個(gè)實(shí)結(jié)點(diǎn)，各部位的材料參數(shù)見表4。    以壩踵為坐標(biāo)原點(diǎn)，x軸指向下游，z軸垂直向上，y軸由右手法則確定，兩側(cè)壩面無y向位移，壩基上下游面無x向位移，壩基底面無z向位移。假設(shè)壩基地應(yīng)力場由壩基重力場產(chǎn)生，壩體一次澆筑完成，水庫一次蓄滿，下游無水，取相對控制誤差能范為5.0%。計(jì)算結(jié)果表明，當(dāng)p=3時(shí)，計(jì)算過程收斂，升階過程的主要指標(biāo)如表5所示。  

17、 整個(gè)升階過程結(jié)束后的單元階次分布圖、屈服單元分布圖、第一主應(yīng)力、第三主應(yīng)力、最大剪應(yīng)力等值線圖(單位：MPa)、位移矢量圖和等安全系數(shù)圖分別見圖6圖12所示。 5  結(jié) 論   (1)在求解含三維節(jié)理單元的有限元問題時(shí)，p型有限單元法是一種切實(shí)有效的方法，它具有很高的收斂速度。(2)用p型有限單元法對含三維節(jié)理單元的水工結(jié)構(gòu)有限元問題進(jìn)行自 適應(yīng)分析時(shí)，不必對單元網(wǎng)格重新剖分加密，在整個(gè)計(jì)算過程中網(wǎng)格可以保持不變。(3)用p型有限單元法求解含三維節(jié)理單元的水工結(jié)構(gòu)有限元問題時(shí)，用較為簡單的網(wǎng)格就模擬常規(guī)有限單元法中較為復(fù)雜的網(wǎng)格。(4)由于采用了升階譜

18、的單元模型，使得后一步的計(jì)算可以承襲前一步計(jì)算的結(jié)果，從而大大地減小了計(jì)算工作量。         參 考 文 獻(xiàn)：1O C Zienkiewicz, D V Phillips. An automatic mesh generation scheme for plane and curved surface by isoparametric coordinatesJ. Int.J.Num. Meth. Eng.,1971,3:519-528.2O C Zieniewicz, J P DE, R. Gago，D. W

19、. Kelly. The hierarchical concept in finite element analysisJ. Computers & structures,1983, 16(1-4):53-65.3Zienkiewicz O C, Zhu J Z, Gong N G. Effective and practical hp version adaptive analysis procedures for the finite element method J.Int. J. Numer. Meth. Eng.,1989,28(1-6):879-891.4陳尚法，陳勝宏.彈粘塑性自適應(yīng)有限元在滑坡穩(wěn)定分析中的應(yīng)用J.巖石力學(xué)與工程學(xué)報(bào)，2001，20(增2)：1596-1599.5陳勝宏，王勁松，張君祿.水工結(jié)構(gòu)的彈粘塑性自適應(yīng)有限元分析J.水利學(xué)報(bào)，1996，(2)：68-75.6Xia H X, Chen S H. 3-D ad