第6章 平面問題高階單元(簡)

第6章平面問題高階單元6.1位移模式階次的選擇在前面兩章中討論了平面問題三結(jié)點(diǎn)三角形單元，其位移模式的最高階是坐標(biāo)x、y的一次項(xiàng)。這種位移模式導(dǎo)致單元常應(yīng)變、常應(yīng)力特性，單元應(yīng)變矩陣、應(yīng)力矩陣、剛度矩陣均為常數(shù)矩陣，因此計(jì)算非常簡單。但這種單元難以反映應(yīng)力梯度的迅速變化。2021/5/91

要想提高計(jì)算精度，必須細(xì)分網(wǎng)格，增加單元數(shù)和點(diǎn)數(shù)，因而加大輸入數(shù)據(jù)的工作量。

提高計(jì)算精度的另一條有效途徑是采用高階單元。由于高階單元的應(yīng)變、應(yīng)力不再是常數(shù)，因此采用少量單元就可能達(dá)到較高的精度。圖7-1懸臂梁分別采用高、低階單元計(jì)算就是一個典型的例子。2021/5/92h4hPAB懸臂深梁解析解：

A=1.0

B=1.0常應(yīng)變單元：

A=0.866

B=0.619高階單元：

A=0.99

B=0.99??????圖6-12021/5/93選擇位移模式時(shí)，第2章提到要考慮解的收斂性，即要考慮到位移模式的完備性和協(xié)調(diào)性。實(shí)際操作中，一般應(yīng)考慮位移模式的對稱性。這是因?yàn)?，有限元位移模式的選擇實(shí)際是以帕斯卡（Pascal）三角形基礎(chǔ)上的(如圖6-2所示)，由低價(jià)至高階，順序選取，組成多項(xiàng)式。多項(xiàng)式中的項(xiàng)數(shù)等于單元節(jié)點(diǎn)自由度數(shù)。如三節(jié)點(diǎn)三角形單元，位移模式取完全一次式，共3項(xiàng)。六節(jié)點(diǎn)三角形單元，位移模式取完全二次式共6項(xiàng)。如果某一階次不能全取，則應(yīng)按對稱性原則適當(dāng)選取。

1xyx2xyy2x3x2yxy2y3

x4x3yx2y2xy3y4

圖6-2多項(xiàng)式選擇的怕斯卡三角形2021/5/94

1xyx2xyy2x3x2yxy2y3

x4x3yx2y2xy3y4

圖6-2多項(xiàng)式選擇的怕斯卡三角形例如在下節(jié)將要討論的四結(jié)點(diǎn)矩形單元中，位移模式不能取1，x，y，x2四項(xiàng)，也不能取1，x，y，y2四項(xiàng)，而應(yīng)取1，x，y，xy四項(xiàng)。6.2四節(jié)點(diǎn)矩形單元圖6-3示出的矩形單元，邊長分別為2a和2b。取4個角點(diǎn)為節(jié)點(diǎn),編號為i，j，l，m。將x軸和y軸置于單元的對稱軸上單元的位移函數(shù)可取為：1、位移函數(shù)2021/5/95

在上式表示的位移模式中，a1,a2,a3,a5,a6,a7，a8反映了單元的剛體位移和常應(yīng)變。在單元的邊界(x=±a或y=±a)上（或），位移是按線性分布的。因此，相鄰單元在公共邊上的位移是連續(xù)的。這樣，位移模式滿足了解答收斂性的充分條件。ijlmxyaabb圖6-3在式（6-1）中代入節(jié)點(diǎn)位移和節(jié)點(diǎn)坐標(biāo)后，可解出（6-1）2021/5/96式中形函數(shù)為：（6-3）（6-2）各待定系數(shù)(a1…a8)。將這些系數(shù)再代入式(6-1),可得：2021/5/97則式（6-3）可簡寫為（6-4）將位移函數(shù)寫成矩陣形式，即有與式（2-20）相同的形式（6-5）式中（6-6）令在節(jié)點(diǎn)上的值為：2021/5/98（6-7）其中，I為二階單位矩陣。2、應(yīng)變矩陣根據(jù)幾何方程，可得與式（2-25）同樣的形式（6-8）把應(yīng)變矩陣[B]寫成子矩陣形式（6-9）

其中（6-10）2021/5/99

由此可見，[B]是

、

的函數(shù)，即是x、y的函數(shù)。因此單元中的應(yīng)變不再是常數(shù)。3、應(yīng)力矩陣根據(jù)應(yīng)力-應(yīng)變關(guān)系，可以計(jì)算單元中的應(yīng)力，得到式（2-28）同樣形式（6-11）應(yīng)力矩陣[S]具有與式（2-29）同樣形式（6-12）將[S]寫成子矩陣形式（6-13）2021/5/910其中（6-14）上式對應(yīng)平面應(yīng)力情形。對于平面應(yīng)變情形，只需將其中的E，

作相應(yīng)的改變即可。4、單元剛度矩陣

單元剛度矩陣可采用式（2-33a）進(jìn)行計(jì)算（2-33a）2021/5/911在四節(jié)點(diǎn)矩形單元中，[k]是一個8×8的矩陣。將[k]寫成分塊形式：（6-16）其中的子矩陣[krs]2×2可由下式計(jì)算2021/5/912（6-17）上式對應(yīng)平面應(yīng)力情形。對于平面應(yīng)變情形，只須將上式中的E、

作相應(yīng)的改變。5、等價(jià)節(jié)點(diǎn)力

單元體積力和表面力引起的節(jié)點(diǎn)力仍可用式（2-45）和（2-46）進(jìn)行計(jì)算。2021/5/913

對本問題給定的位移函數(shù)，若體積力是重力的情形（設(shè)重度為

），單元等價(jià)節(jié)點(diǎn)載荷列陣為：（2-45）（2-46）（6-18）

有了對單元的上述結(jié)果，便可應(yīng)用第5章的方法組集結(jié)構(gòu)剛度矩陣和節(jié)點(diǎn)荷載向量；求解節(jié)點(diǎn)位移；計(jì)算內(nèi)力和應(yīng)力。

四節(jié)點(diǎn)矩形單元采用較高階的位移模式，具有比三節(jié)點(diǎn)三角形單元較高的計(jì)算精度。但矩形單元也有缺點(diǎn)，2021/5/914

在三角形單元i,j,m的各邊中點(diǎn)增設(shè)一個節(jié)點(diǎn)，使每個單元具有6個節(jié)點(diǎn)，得到圖6-4所示的六節(jié)點(diǎn)三角形單元。這種單元具有12個自由度，可以采用完全二次多項(xiàng)式的位移模式：

一是不能適應(yīng)斜線及曲線邊界，二是不便于采用大小不同的單元。6.3六節(jié)點(diǎn)三角形單元1、位移模式

???ijmi

j

m

xy圖6-42021/5/915（6-20）所取位移模式反映了單元的剛體位移和常應(yīng)變；單元內(nèi)部是連續(xù)的；在單元邊界上位移分量按拋物線變化，而每條公共邊界上有3個公共結(jié)點(diǎn)，可以保證相鄰兩單元位移的連續(xù)性。因此，上述位移模式滿足收斂的必要和充分條件。

上述位移模式確定之后，可以用分析三節(jié)點(diǎn)三角形單元和四節(jié)點(diǎn)矩形單元相同的方法進(jìn)行分析。得到形函數(shù)、應(yīng)變矩陣、應(yīng)力矩陣、單元剛度矩陣、等價(jià)節(jié)點(diǎn)力向量。但其過程十分繁復(fù)，采用面積坐標(biāo)可以大大簡化計(jì)算。2021/5/9162、面積坐標(biāo)

對于一個三角形ijm（圖6-5），三角形內(nèi)任一點(diǎn)P(x,y)的位置，可以用如下的三個比值來確定：ijmxy圖6-5·P（6-21）AiAjAm（1）定義2021/5/917其中A為三角形ijm的面積，Ai,Aj,Am分別為三角形的Pjm,Pmi,Pijd的面積。這三個比值Li,Lj,Lm稱為P點(diǎn)的面積坐標(biāo)。由于則（6-22）由此可見，P點(diǎn)的三個面積坐標(biāo)不是獨(dú)立的。同時(shí)，面積坐標(biāo)只是用以確定三角形內(nèi)部某點(diǎn)的位置，因而是一種局部坐標(biāo)。下面進(jìn)一步給出面積坐標(biāo)的幾個性質(zhì)。（2）面積坐標(biāo)與直角坐標(biāo)的關(guān)系在圖7-5中，三角形Pjm的面積為2021/5/918（6-23）由式（6-23），式（6-21）化為（6-24）將式（6-24）（6-23a）和式(2-18)(2-17)對比，可知，面積坐標(biāo)就是三節(jié)點(diǎn)三角形單元的形函數(shù)（6-23a）Ni、Nj、Nm。2021/5/919

將式（6-24）的3個式子分別乘以xi,xj,xm，然后相加，并利用關(guān)系式（6-23a），有同理（6-25）（3）面積坐標(biāo)的導(dǎo)數(shù)公式

根據(jù)面積坐標(biāo)與直角坐標(biāo)的關(guān)系，由復(fù)合函數(shù)的求導(dǎo)公式，有2021/5/920（6-26）（4）面積坐標(biāo)的積分公式下面給出面積坐標(biāo)的冪函數(shù)積分公式。它們在計(jì)算單元剛度矩陣和等效結(jié)點(diǎn)載荷時(shí)有用。2021/5/921在三角形單元上進(jìn)行積分時(shí)，有（6-27）在三角形某一邊（設(shè)ij邊，邊長為l）上進(jìn)行積分時(shí)，有（6-28）3、用面積坐標(biāo)表示六節(jié)點(diǎn)三角形單元計(jì)算公式

對應(yīng)如圖6-4所示的六節(jié)點(diǎn)三角形單元，形函數(shù)可用面積坐標(biāo)表示為

（1）形函數(shù)和位移表達(dá)式2021/5/922???ijmi

j

m

xy圖6-6現(xiàn)利用形函數(shù)的性質(zhì)檢驗(yàn)式（6-29）的正確性。先考慮三角形的角點(diǎn)，例如圖6-6中的i點(diǎn)，有由式（6-21）（P16），有代入式（6-29），有（6-29）2021/5/923

再考慮三角形的邊中點(diǎn)，例如i

點(diǎn)，面積劃分如圖6-7所示。顯然有：???ijmi

j

m

xy圖6-7由式（6-21）(P16)，有代入式（6-29）(P16)，進(jìn)一步說明式（6-29）所表示的形函數(shù)的正確性。說明形函數(shù)Ni在i點(diǎn)等于1，在其它節(jié)點(diǎn)等于0，因此是正確的。2021/5/924形函數(shù)確定后，單元中任意一點(diǎn)的位移可以表示為：（6-30）其中（6-31）（6-32）其中I為二階單位陣，形函數(shù)由式（6-29）確定。（2）應(yīng)變矩陣單元中的應(yīng)變?nèi)钥杀硎緸椋海?-33）2021/5/925

式中應(yīng)變矩陣[B]為：（6-34）其中（6-35）),,;,,()(4)(4)(400)(421][23mjimjibLLbcLLccLLcbLLbABmjmjmjmjmjmjmjmji￠￠￠úúú?ùêêê?é++++=′￠2021/5/926單元中的應(yīng)力仍可表示為：（3）應(yīng)力矩陣（6-36）

式中[D]是彈性矩陣，由式（2-9）確定；應(yīng)變矩陣由式（6-34）、（6-35）確定。根據(jù)矩陣乘法，可以給出用面積坐標(biāo)表示的應(yīng)力矩陣[S]（4）單元剛度矩陣

單元剛度矩陣仍可表示為：（6-37）

根據(jù)[B]、[D]的表達(dá)式以及面積坐標(biāo)的積分公式（6-27），可以求出[k]中元素的顯式表示。由于較為繁復(fù)，這里就不列出詳細(xì)結(jié)果。

2021/5/927（5）等價(jià)節(jié)點(diǎn)力向量

由于位移模式是非線性的，因此體積力和表面力引起的節(jié)點(diǎn)力向量不能采用靜力等效原理進(jìn)行分配，而應(yīng)采用相應(yīng)公式進(jìn)行計(jì)算。單元體積力引起的等價(jià)節(jié)點(diǎn)力計(jì)算公式仍為：（6-38）

將由式（6-29）、（6-32）表示的[N]代入，并應(yīng)用積分式（6-27），可以計(jì)算

FV

e。例如對于重力引起的

FV

e

，有

它表示各邊中點(diǎn)承擔(dān)單元重力的1/3。2021/5/928單元表面力引起的結(jié)點(diǎn)力計(jì)算公式仍為：

（6-39）設(shè)在ij邊上受有x方向的均勻分布力ps，對應(yīng)的等價(jià)節(jié)點(diǎn)力向量為（圖6-8）pslh/6pslh/64pslh