1、    含節理單元的三維P型自適應有限元解法            作者：費文平，陳勝宏時間：2007-11-25 12:11:00                     摘要：本文提出了含三維無厚度節理單元、等厚度節理單元和變厚度節理單

2、元的p型自適應有限元模型，給出了三維節理單元升階譜有限元法的解題步驟，通過具體算例，驗證了p型升階譜有限元法在求解含三維節理單元的有限元問題時的可行性及優越性。 關鍵詞：節理單元 升階譜有限元 p型有限元     有限元法是求解微分方程數值解的一種重要方法，對于一個給定的問題，為改善其有限元解的精度，可以采用以下3種方法。(1)h型有限元法1，這種方法通過減小單元尺寸來提高有限元解的精度。(2)p型有限元法2，這種方法通過增加基底函數的階次來提高有限元解的精度。(3)hp型有限元法3，這種方法是以上兩種方法的綜合，它既減小單元尺寸，又增加基底函數的階次。

3、作者所在的研究小組從1995年開始研究水工結構的h型彈粘塑性有限單元方法，目前已建立了實用的二維分析軟件體系4,5，并在三維分析方面取得了進展6。從1999年開始，在水工結構的p型自適應分析方面也有所突破，1999年，程昭7等人針對水工結構分析問題提出了三維升階譜有限元分析方法。2001年，陳勝宏8等人進一步提出了二維問題的p型自適應分析策略，并將自適應有限元方法歸類為全域升階方法、單元升階方法和自由度升階方法等三類。之后，費文平9等人將p型自適應有限元分析方法推廣到三維彈粘塑性領域。但是，以上有關p型有限元的研究成果中均未涉及到斷層、節理這一類特殊單元。   

4、大壩壩基、壩肩和巖石高邊坡等部位總是存在斷層、節理和軟弱夾層等大規模的不連續面，且對結構的變形和穩定影響巨大，故在有限元分析中應給予高度的重視。古德曼(Goodman)最初運用有限元技術模擬巖體工程中的非線性不連續面問題，并提出了無厚度的節理元的概念10。隨后，朱伯芳于1979年提出了等厚度節理元模型，并將其與無厚度節理元模型形成統一的計算公式11，在此基礎上，王鴻儒等人提出了變厚度的節理單元的彈塑性模型并將其應用到工程實踐中12。目前，國內外在有關p型自適應有限元分析的研究中，尚未涉及到這類特殊單元的處理問題，從而使研究成果的工程應用受到一定程度的限制。   

5、對于常規塊體單元的三維p型有限元模型，作者在文獻9中已有詳盡的論述，本文主要給出三維無厚度節理單元、等厚度節理單元和變厚度節理單元的p型有限元模型，并給了升階譜的計算格式。實例分析結果表明，用p型有限元法來求解含三維節理單元的有限元問題具有收斂速度快、計算精度高的優點。1  三維p型無厚度節理單元和等厚度節理單元模型    對厚度很小和厚度變化不大的節理，可以分別采用無厚度的節理單元和等厚度的節理單元進行模擬，可取如圖1所示的六面體節理單元，進行單元的網格劃分。1.1  形函數及位移函數  對節理單元上下面的相應點、棱和面可取相同的點

6、基函數、棱基函數和面基函數，對無厚度節理單元或等厚度節理單元，不存在連續上、下兩面的棱基函數和面基函數，也不存在體基函數。基底函數的具體形式如下13：    點基函數(p1)(1)式中：，。    棱基函數(p2)，(2)    面基函數(p4)（i+j=p,  i,j2）(3)式中：而為Legender多項式。令位移函數為(4)(5)(6)    同理可以寫出V下，V上，W下，W上及v，w的具體表達式。    將基函數Ni,pEi,p

7、F統一記為i，位移差(uN,i+4-uNi)，(uE,i+4-uEi),(uF2-uF1)統一記為廣義結點位移差ui，設單元基底函數個數為fe(p)，上式簡化為，同理有,1.2  三維等厚度節理單元或無厚度節理單元升階過程  當p=1時：N1-1I-2I-3I-4I1I2I3I4I,I為3×3階的單位陣。當p=2時，N2可在N1的基礎上進行擴充，擴充矩陣Np=2為Np=2=-5I-6I-7I-8I5I6I7I8I。依此類推，可得最終的形函數矩陣為(7)1.3  坐標插值及坐標變換  對于節理上、下面坐標的插值仍采用各面上的四個節點進行插值，即(

8、8)(9)式中：(xi,yi,zi)i=18為節理間面體單元8個頂點的整體坐標。    定義整體坐標系的x軸朝北，y軸朝西，z軸朝上，定義等厚度的節理單元或無厚度的節理的局部坐標系的z為中面的法線朝上方向，y指向節理面的傾向，x軸由右手法則確定，并設等厚度節理的傾角為，傾向為。    三維等厚節理或無厚節理單元局部坐標與整體坐標的轉換矩陣為(10)1.4  三維等厚度節理的單元剛度矩陣(11)(12)式中：彈性矩陣；單元應變矩陣B'LB1/eLNp根據虛功原理，單元剛度矩陣為(13)1.5  三維無厚度節理

9、的單元剛度矩陣  當等厚度節理單元的厚度e0時，即形成無厚度的節理單元。此時，可假定單元內應力分量與位移差成正比，同理可得單元剛度矩陣為(14)式中：為單元勁度矩陣。2  三維p型變厚度的節理單元模型     當節理單元的厚度變化較大時，應將等厚度節理單元推廣得到變厚度的節理單元，如圖2所示建立局部坐標系。2.1  形函數及位移函數  六面體變厚度的節理單元的基底函數，由點基函數、棱基函數、面基函數和體基函數組成，各基底函數的具體形式如下：點基函數(p1) Ni(,)=1/8(1i)(1+i)(1+i)  (i=

10、1，2，8)(15)式中：i=(-1)i,i=(-1)i/2+0.5，i=(-1)i/4+0.75。    棱基函數(p2)pEi=1/4(1+1i)(1+1i)(1+1i)(2ip()+2ip()+2ip()，(i,1,2,12)(16)式中：2i=1-i/12+0.65，2i12(1-(-1)i-1/4)；2i=i-18，將1i，1i,1i用向量的形式表達為1=0000-1-111-11-11T，1-11-110000-1-11T，1-1-111-11-110000T。        面基函

11、數(p4)(17)(18)(19)式中：i,j2，i+j=p。    體基函數(p6)：(i,j,k2,i+j+k=p)(20)令位移函數(21)    將基底函數Ni,pEi,pFi,pB統一記為i，位移uNi,uEi,uFi,uB統一記為廣義結點位移ui，設單元基底函數的個數為fe(p)，上式簡化為，同理。2.2  三維變厚度節理單元升階過程  三維變厚度節理單元可按表1的方式進行升階。表1  基底函數與階次的關系階次123456基底函數1-81-201-321-501-741-1052.3 

12、; 坐標插值與坐標變換  對于單元任一點的坐標，仍采用單元的8個實結點進行插值，即，此處i=Ni(i=1,2 ,8),(xi,yi,zi)為立方體單元的8個頂點坐標。參見圖2，取局部坐標系的xoy平面與變厚度六面體單元的中面重合，z軸垂直于中面。根據變厚度節理單元8個實結點的整體坐標，可按如下方法建立坐標轉換矩陣及局部坐標系。    先求出變厚度六面體單元中面法向向量nx,y,z，其中x(y2+y6-y1-y5)(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)。    轉換矩陣L的形成，  令l1=(x2+x6-x1-x5)/l,m1=(y2+y6-y1-y5)/l,n1=(z2+z6-z1-z5)/l(22) ,(23)    向量l2 m2 n2T可根據右手法則，由向量l1 m1 n1T和l3 m3 n3T確定，即l2=m3n1-m1n3,m2=n3l1-n1l3,n2=l3m1-l1m3 (24)2