空間問題的四面體單元

1、圖3-7空間四面體單元第三章 軸對稱、三維和高次單元 3-2空間問題的四面體單元空間問題的有限單元法，和平面問題及軸對稱問題的有限單元法的原理和分析過程完 全相同。由于空間問題應(yīng)采用三維坐標(biāo)系，因此單元的自由度、剛度矩陣的元素個數(shù)，方 程組內(nèi)方程個數(shù)等要較平面問題和軸對稱問題多，所以空間問題的規(guī)模一般比軸對稱問題 和平面問題大得多。它要求計算機的內(nèi)存大，且計算時間長，費用高。這些問題都給三維 有限單元法的具體運用帶來許多困難。和平面問題一樣，空間有限單元法采用單元 也是多種多樣的，其中最簡單的是四節(jié)點四面體 單元。采用四面體單元和線性位移模式來處理空 間問題，可以看作平面問題中三角形單元的推廣

2、。在采用四面體單元離散化后的空間結(jié)構(gòu)物 中，一系列不相互重疊的四面體之間僅在節(jié)點處 以空間鉸相互連接。四節(jié)點四面體單元僅在四個 頂點處取為節(jié)點，其編號為i,j,m,p 。每個單元的 計算簡圖如圖3-7所示。在位移法中，取節(jié)點位移為基本未知量，四 節(jié)點四面體單元共有十二個自由度（位移分量），其節(jié)點位移列陣為VUiViwiUj VjwjUmTWm Up Vp Wp其子矩陣Ui Vi Wi (i,j,m)相應(yīng)的節(jié)點力列陣為FiFjFmFp其子矩陣FiUiViw一、單元法位移函數(shù)結(jié)構(gòu)中各點的位移是坐標(biāo)X、y、z的函數(shù)。當(dāng)單元足夠小時，單元內(nèi)各點的位移可用簡單的線性多項式來近似描述，即u1 2 X3y4

3、Zv56 X7y8Z(3-49)w0 10Xny12Z曰2，12是卜二個待定系數(shù)，它們可由單元的節(jié)點位移和坐標(biāo)確定。假定節(jié)點 i,j,m,p的坐標(biāo)分別為(xi yiZi )、(xjyj zj)、(Xm將它們代入(3-49)式的第一式可得各個節(jié)點在X方向的位移Ui12Xi3Yi4Zuj12Xj3Yj4ZjUm12Xm3Ym4 ZmUp12Xp3Yp4 Zp解上述線性方程組，可得到1 ,2 ,3 ,4 ,再代入U 6V(ai bXcydiZ)Ui(ajbjx(am bmX Cmydmz)Um(ap bpXCp(3-50)式，得ymZm )、(Xp y p Zp )，5y 3)5(3-51)式中1

4、,其中V為四面體ijmp的體積，a,bi,，cp,d p為系數(shù)。(3-50)y dpZ)Up1XiYiZi1XjyjZj1XmymZm1XpYpZp(3-52)XjyjZj1 yj召Xmy mZmb1 ymZmXpypZp1 ypzpaiXj 1 ZjXj yj 1Xm 1 ZmdiXmym 1Xp 1 ZpXp yp 1(i,j,m,p)(3-53)Cii,j,m,pi tj t m的轉(zhuǎn)向轉(zhuǎn)動時,必須按照一定的順序：向 p的方向前進，象圖為了使四面體的體積 v不致為負(fù)值，單元四個節(jié)點的標(biāo)號 在右手坐標(biāo)系中，要使得右手螺旋在按照3-1中單元那樣。用同樣方法,可以得出其余二個位移分量:16V(a

5、m16V(am(aibixcy diZ)v (ajbjX Cjy djZ)VjbmXCmYdmZ)Vm (ap bpX CpydpZ)Vp(aibixqy djZ)w (aj bjX ciY djZ)WjbmXCmYdmZ)Wm (ap bpX CpydpZ)Wp(3-54)(3-55)綜合表達式(3-51)、(3-54)及(3-55)，可以將位移分量表示成為Tefu v w N INi INj INm IN pe(3-56)其中1是三階的單位矩陣，N為形函數(shù)矩陣，而各個形函數(shù)為Ni(ai bix ciy diz)/6V(i,m)(3-57)Nj(a bix c y diz)/6V(j, p)

6、和平面問題相似，(3-49)式中的系數(shù)12代表常量的正應(yīng)變；其余1 ,5 ,6代表剛性移動U0 , V0 , Wo ;系6個系數(shù)反映了剛性轉(zhuǎn)動Wx , Wy , Wz和常量剪應(yīng)變。這就是說，12個系數(shù)充分反映了單元的剛體位移和常量應(yīng)變。同時，可以證明： 由于位移模式是線性的， 兩個相鄰單元的共同邊界在變形過程中，始終是相互貼合的，使 得離散的模型變形中保持為連續(xù)體。這樣，選用的位移函數(shù)滿足收斂的充分必要條件，保 證了有限單元法解答收斂于精確解。二、載荷移置空間問題的單元載荷移置和平面問題一樣，也是根據(jù)靜力等效原則，將不作用在節(jié)點 上的集中力、體力、面力移置成作用在節(jié)點上的等效節(jié)點載荷。其通用公

7、式的形式和平面 問題也是一樣的，只不過多出一維空間分量。1.集中力T設(shè)單元上某點(x,y,z)作用有集中力PPxPyP,則仍然得到等效節(jié)點載荷R NT P(3-58)這里eRXi Y 乙Xj YjZ jX mYmZ mXp Yp Zp2.分布體力單元上作用有分布體力P XY Zt，則R eNT PdV(3-59)其中dV是單元中的微分體積，對于直角坐標(biāo)糸上式為ReNt p dxdydz(3-60)e3.分布面力單元的某一邊界面s上作用有一分布面力P X Y Z TR e Nt P dA其中dA是邊界面S上的微分面積。4.常見載荷的移置上列公式是空間問題載荷移置的通用公式。 對于四節(jié)點四 面體單

8、元，由于其采用線性位移模式， 采用直接計算虛功的方 法求出節(jié)點載荷比較簡單。下面介紹常見的二種載荷的移置。重力四面體單元的自重為 W作用在質(zhì)心C處(如圖3-8)。為 求得節(jié)點載荷 X,Yi ,Zi,可分別假想發(fā)生 u* 1 , V* 1或 w；1的虛位移。在Ui* 1或V；1時，整個單元上各點的均沒有 z方向上圖3-8重力移置當(dāng)Wi 1時，jmp面上各點的虛位移為零，即* 1Wb0，又因bc bi ,所以有4的虛位移，重力 W不做功，所以 X=Yi=O。*1“Wwc-，Zi44對于其余三個節(jié)點可得同樣結(jié)論，于是有(3-61)eRi00TW(i,j,m,p)41移置到每個節(jié)點上即可。4P,共在三

9、個節(jié)點上的強度分別為即，對于四節(jié)點四面體單元承受的重力載荷，只需要把共(2)界面壓力設(shè)四面體的一個邊界面ijm上受有一線性分布的壓力作用于1/4。于ijm面上的d點， 二是可得d點到ij邊和im邊的距離分別為m 至U ij及j到im邊的距離的TTePPP c1,11Ri-0q ijm1 -0(3-62)2446 j22qi,0,0。很容易看出，該力向p點移置的等效節(jié)點力為零。由水力學(xué)知，總壓力P 1 qi ijm ,3所得各節(jié)點載荷的方向和分布力的方向相同，要求各節(jié)點載荷分量還需乘上相應(yīng)的方向余 弦。由上述面力移置結(jié)果，可求出任意線性分布面的等效節(jié)點載荷。如在ijm面受有線性分布面力在各點強度

10、分別為qi, qj，qm,時，在i節(jié)點的等效載荷為111P (qiqjqm) jm (i ，j， m) (3-63)6 2 2三、應(yīng)力應(yīng)變矩陣空間問題幾何方程為yTzxyuvwxyzuvyxu zwywuxzz將四面體單兀之位移表達式(3-52)、(3-54)和(3-55)代入幾何方程，即得單兀應(yīng)變。用節(jié)點位移可表示為eBE3iBjBmBpe(3-64)式中應(yīng)變矩陣子矩陣為 6X 3矩陣:bi000Ci01 00diBic (i，j，m,p)(3-65)6V cibi00diqdi0bi由上式可以看出，每一個單元的應(yīng)變矩陣是一個常量矩陣；因此，采用線性位移模式 的四面體單元是常應(yīng)變單元。這與平

11、面問題中的三角形單元是一樣的。而與平面問題的不 同之處僅在于應(yīng)變矩陣的階數(shù)不同。將表達式(3-16)代入空間問題的物理方程，即可得出用單元節(jié)點位移表示的單元應(yīng)力:eeD DB S(3-66)式中彈性矩陣D為應(yīng)力矩陣D1100122(1)000000稱1 22(1 )01 22(1 )SA1Si Sj Sm Sp(3-67)A2biAbSiE(1)Ab6(1)(1 2 )V A2 GA2diAACiAAdi(i,j,m,p)(3-68)A2bi0A2diA2G0A2biS是常量矩陣，所以，四面體單元是顯然，式(3-68)中各元素均為常量，應(yīng)力矩陣 常應(yīng)力單元。四、單元剛度矩陣空間問題的單元剛度由

12、虛功方程導(dǎo)出。假設(shè)該單元發(fā)生某虛位移，相應(yīng)節(jié)點虛位移為e。此時相應(yīng)的虛應(yīng)變?yōu)锽將上式及式(3-66)代入虛功方程，有* e t e* e te()F(B)DB dxdydzv通過與平面問題一樣的處理，并注意到矩陣B中的元素為常量，可以得到F e BTDBdxdydz e BTDBeV Ke e (3-69)v式中，Ke為單元剛度矩陣:KeBTDBdxdydz BTDBV(3-70)e將式(3-64)和(3-68)式代入，可以得出KHKijKimKipKe心KjjKjmKmiKmjKmmKmpKpiKpjK pmKpp(3-71)其中，Krse為3X 3階方陣:KrseE(1)36(1)(1 )

13、VbRs A2(g drds)AGbsA2brCsAdRsAbAAbrCsACrCsCrCsA2(bQs drds)AdisAcdAbrds AzdrbsACrds AzdrCsdrds A2(brbs CrCs)(r,s=i,j,m,p)(3-72)有了單元節(jié)點力和節(jié)點位移之間的關(guān)系之后，通過分析每個節(jié)點的平衡條件可得到Krse se s i,j ,m, pRre這個矩陣形式的方程實際上代表了關(guān)于r節(jié)點三個坐標(biāo)軸方向的力平衡方程式。將關(guān)于結(jié)構(gòu)物所有節(jié)點的線性方程式集合起來，可以得到K式中 代表整個結(jié)構(gòu)的節(jié)點的位移， 是所求之基本未知量； R 代表整個結(jié)構(gòu)的節(jié)點載荷； K 為整體剛度矩陣，其是

14、由每個單元剛度矩陣升階后組集得到，即NEeKKee1其為3NP階方陣。顯然，對每一個子矩陣，應(yīng)有NEK rs K rs e1和平面問題一樣，K是對稱、帶狀、稀疏矩陣，在消除剛體位移之后，它是正定的。由平衡方程組可以解出節(jié)點位移，隨后即可求得所需節(jié)點和單元應(yīng)力。五、形成四面體的對角線劃分方法在實際計算中，用一系列的四面體來組合成一個空間物體，這個形象是很難想象的。但是如果先用一系列較為直觀的六面體 （圖 3-9） 來劃分彈性體，然后由計算機來將這些六 面體及三棱柱劃分為若干個四面體，則要方便得多。同時減少許多準(zhǔn)備及輸入工作，也為 將來結(jié)果分析帶來方便?，F(xiàn)在介紹一種適合計算機進行自動劃分四面體的方

15、法對角線劃分法。1. 將六面體劃分為四面體的方法通過連接六面體上一些四邊形的對角線， 可以把一個六面體劃分為五個或六個四面體。 為敘述方便，先將六面體的八個角點進行局部編號，編號原則是先頂面后底面，對于頂面 或底面的節(jié)點來說，則是先前后后，從左到右 （ 見圖 3-9） 排列。六面體和三棱柱(1) 將一個六面體劃分為五個四面體這種方法是先過六面體的一些四邊形的對角線，從六面體的四個角上切下四個四面體，最后剩下中心的一個四面體，共得五個六面體單元。選擇被切下的角點不同，有二種不同 的劃分結(jié)果，如圖 3-10(a )和(b )所示。我們分別稱之為 A5型劃分和B5型劃分。A5型劃分所得五個四面體為

16、1246, 1347, 1467, 1567, 4678;而B5型劃分則得到1235，2348, 2358 , 2568, 3578五個四面體。以上二種劃分方法的共同特點是，六面體二對面四邊形的剖分對角線是交叉的。這就使得如果一個六面體按A5型劃分，那么與之相鄰的各個六面體必定要按B5型劃分。(a)(b)圖3-10六面體劃分為五個四面體(a) A5 型剖分； (b)B5 型剖分(2) 將一個六面體劃分成六個四面體將六面體劃分成六個四面體有很多種劃分方法。這里介紹兩種，如圖3-11所示。它們的共同特點是，六面體上兩對面四邊形的剖分對角線是“平行”的。所不同的是在A6型剖分中取大對角線36作為劃分

17、線，而在 B6型中則是取大對角線 45作為劃分線。為清楚起見，可將 A6型劃分理解為先將六面體沿2367分成兩個三棱柱，再將每個三棱柱分成三個四面體，分別得到1235，2356，3567和2346，3467, 4678六個四面體(見圖3-12。當(dāng)然，A6型劃分也可看成先將六面體沿3456面剖分，得到兩個不同于前的三棱柱，但最后得到的六個四面體是相同的(圖3-13)。對于B6型劃分，六面體先以折面 2457為分界面拆分成兩個三棱柱，如圖3-14所示。于是可見，每一個“三棱柱”被劃分為三個四面體，它們分別是 1235，2345，3456和2456，4567，4678。同時，也不難證明，若以3456

18、為分界面按B6型劃分將六面體拆成的二個“三棱柱”雖與前面的不同，但是劃分成的六個四面體和前面得到的完全相同。(a)(b)圖3-11六面體劃分為六個四面體(a) A6 型剖分； (b)B6 型剖分圖3-12 A6型劃分，以折面2376為兩個“三棱柱”的分界面圖3-13 A6 型劃分，以折面3456為兩個“三棱柱”的分界面圖3-14 B6 型劃分拆成兩個“三棱柱”我們看到，A6型和B6型劃分，由于其相對四邊形的對角線“平行”，而剖分大對角線35和45并不在六面體表面上，其相鄰的六面體可以全部采用A6型劃分或B6型劃分，兩種劃分也可以交替使用。一個六面體劃分為六個四面體，各四面體的體積大小一般較為均

19、 勻；但是在相等的六面體數(shù)目下，A6型和B6型劃分所產(chǎn)生的四面體單元的總數(shù)，要比A5型和B5型產(chǎn)生的多六分之一。此外，在A6型與B6型劃分中，如果離散體的節(jié)點整體編號是按本節(jié)開頭所述，從上 到下，從左到右連續(xù)進行的；同時每個六面體八個節(jié)點的整體編號的大小次序與其局部編 號的大小次序相一致(由小到大)的話，那么在劃分中，對底面上任一節(jié)點，與它構(gòu)成四面 體的三個節(jié)點中的最小號碼，不會比其正上方那個節(jié)點的號碼更小。這是由于在A6型及B6型劃分中，注意到在連各個四邊形對角線時，不使節(jié)點編號之差較大的三個節(jié)點出現(xiàn)在 同一個三角形中的結(jié)果。例如圖3-11中，對于1357四邊形，我們連接了 35對角線而使1

20、、7兩節(jié)點分別屬于兩個三角形中。這樣的劃分能獲得一個帶寬較窄的剛度矩陣。(3) 編號推算如果將六面體的八個頂點的節(jié)點整個編號置于數(shù)組D 1:8 中，而將前述圖中的局部編號18理解為數(shù)組D 1:8 的下標(biāo)時，于是上述問題就轉(zhuǎn)化為：要在有八個元素的數(shù)組中按一定規(guī)律， 每次取四個元素構(gòu)成一個四面體單元的節(jié)點編號問題。 對于 A6 型和 A5 型劃分所得四面體頂點編號的規(guī)律性進行一些分析之后，可以導(dǎo)出下列公式，分別表示 按預(yù)定規(guī)律劃分成的四面體的各節(jié)點編號：Dm(1 3(l1)J(1 m)(lJ)Dm(3l J1)(1 m)(1l J(J 1)/2)Dm(2(1 lJ)J(l J)(1 m)(2 l

21、J(5 J)/2)(3-73)Dm(5 lJ(3J) /2 (1m)(4 l J)(l 1,2; J0,1,2;m 0,1)通過直接代入數(shù)字檢驗，知道m(xù)=0對應(yīng)著A6型劃分所得的六個四面體，m=1則對應(yīng)著A5型劃分所得的五個四面體，(此時l=J=2形成的四面體應(yīng)舍去)。例如m=0的情況，當(dāng)1=1 ， J=0， 1， 2時，得到 D1D2D3D5， D2D3D5D6， D3D5D6D7 三個四面體，而 l=2，J= 0, 1, 2 時，得到 D2D3D4D6，D3D4D6D7，D4D6D7D8三個四面體，與前述結(jié)果一致。對于 B6 型及 B5 型劃分,同樣可以導(dǎo)出一個相似的計算公式。但是也可以利用(3-73)式，只需將原來 D 1： 8中元素位置作一定更動。更動的辦法是，對于 B6型劃分，將8、 6、4、2位置的元素置于 1、2、3、4位置上,將 7、5、3、1 位置的元素置于 5, 6, 7, 8 位置