彈性力學平面問題的有限元法

一、    平面問題的基本方程

平面問題是指彈性體內一點的應力、應變或位移狀態只與兩個坐標方向的變量有關的二維問題。

1.  平面應力和平面應變

應力分量：

 式3-1

應變分量：

   式3-2

應力分量與應變分量間的關系：

   式3-3

或

   式3-4

式中：[D]—彈性矩陣

對于平面應力問題

   式3-5a

對于平面應變問題（E換成，μ換成）

   式3-5b

2.  平面問題的基本方程

平面問題的總位能表達式

當δΠ=0，可以得到用位移u和v表示的基本方程

   式3-6

如采用應力函數，可得到平面問題的基本微分方程（雙調和方程）

   式3-7

用有限元法求解平面問題的思路：

①   剖分和插值

把整個平面區域S用三角形板單元或矩形板單元等進行剖分并在單元內進行位移函數（形狀函數）的插值。

   式3-8

②   單元分析

把形狀函數代入位能泛函式Π_i，并按單元進行計算。

   式3-9

③   單元組集

把各單元重新組集起來。

   式3-10

說明：{q}是單元各節點的位移列陣；[K]是單元剛陣；{F}是單元的廣義載荷列陣；q是整個區域上各單元節點的位移總和的列陣；K是總剛陣；F是整個區域上的廣義載荷列陣，S是單元總數。

二、    位移函數

①   設定的位移函數是泛函的極限條件，即控制方程的近似解。

②   選擇位移函數的階次應考慮下列因素：

ⅰ、滿足完備性和協調性。一般采用一個由低階算起完全的多項式表示位移函數。如u=a₁+a₂x+a₃y+a₄x²+a₅xy+a₆y²……

ⅱ、對稱性即該多項式位移函數與局部坐標系的方位無關。

ⅲ、多項式的項數與節點自由度相等。

ⅳ、收斂性

③   設定位移函數時應符合

ⅰ、在單元內部和邊界上（包括節點處）處處都能滿足力的平衡條件和變形協調條件。

ⅱ、在單元內部要求應變或應力最少應是常值（或線性變化的）。

ⅲ、包含有代表剛體運動的項。

一個單元內各點的位移實際上包含著兩部分，一是單元本身變形引起的部分；二是剛體位移部分。位移函數必須能反映這兩種位移。

三、    三角形單元分析

1.  位移函數d9b9f79xv7n

圖3-1 三角形板單元

設三角形板單元單元內某一點的位移函數為

代數形式：    式3-11a

矩陣形式：    式3-11b

那么單元的三個節點i、j、k可以寫成

    式3-12

由式3-12解得，

    式3-13

式中：

為了不使A為負值，i、j、k的順序必須是逆時針方向，如圖3-1。

將代人中，整理得

    式3-14

式中：[N]形狀函數

    式3-15

2.  應變{ε}

由，可得

    式3-16

式中：[B] 應變矩陣

    式3-17

3.  應力{σ}

根據虎克定律，得

    式3-18

式中：[R] 應力矩陣

對于平面應力問題

  式3-19a

對于平面應變問題

 式3-19b

4.  剛陣[K]

根據虛位移原理，可推得

    式3-20

式中：[K] 剛度矩陣

    式3-21

式中：

對于平面應力問題

對于平面應變問題

四、    矩形單元分析

1.  位移函數d9b9f79xv7n

圖3-2 矩形板單元

設矩形板單元單元內某一點的位移函數為

   式3-22

根據矩形板四個坐標值，求得α₁、α₂、α₃、α₄、α₅、α₆、α₇、α₈，并回帶式4-1中，可以得到新的位移函數

   式3-23

2.  應變{ε}

   式3-24

式中：

3.  應力{σ}

   式3-25

式中：對于平面應力問題

對于平面應變問題

4.  剛陣[K]

對于平面應力問題

  式3-26a

對于平面應變問題

  式3-26b

五、    形狀函數

對于平面梁位移，可以用下述形狀函數來表示

v(x)=[N]{q}=N₁v₁+N₂v₂+N₃v₃+N₄v₄

對于平面問題三角形單元，可以用下述形狀函數來表示

u=N_iu_i+N_ju_j+N_ku_k+N_lu_l

v=N_iv_i+N_jv_j+N_kv_k+N_lv_l

形狀函數的幾何意義反映了單元體的變形情況，即位移分布狀態。

形狀函數的兩個重要性質：

1、

2、在單元任一點上三個形狀函數之和等于1。

六、    三角形面積坐標

面積坐標就是用面積的比例關系來表示三角形單元中任意一點P(x,y)的位置。

其中：

圖3-3

當采用面積坐標時，三角形單元內某點的位移函數為

  式3-27

面積坐標與直角坐標之間的關系

  式3-28a

  式3-28b

七、    載荷的轉移

1、三角形單元上載荷的轉移

圖3-3   圖3-4  圖3-5  圖3-6

載荷1：作用在單元i-j邊上的一個沿著x方向的集中載荷P（圖3-3）。

載荷2：作用在三角形單元i-j-k內的一個沿著x方向的集中載荷P（圖3-4）。

載荷3：作用在三角形單元i-j-k的i-j邊上為一個按三角形分布的分布載荷p（圖3-5）。

載荷4：作用在三角形單元i-j-k的i-j邊上為一個按梯形分布的載荷p1、p2（圖3-6）。

①單元i-j-k的任意一點受集中載荷P，則單元等效節點載荷{F}

②單元i-j-k的任意一面受體力F，則單元等效節點載荷{F}_V

③單元i-j邊上受分布表面力F，則單元等效節點載荷｛F｝V

2、矩形單元上載荷的轉移

①矩形單元自重作用

圖3-4

②矩形單元內某點作用集中力

圖3-5

用上述方法也可以計算簡支梁的等效節點載荷值。

八、    平面問題的有限元解法

平面問題的有限元解一般方法和步驟：

1、根據單元的性質和精度要求，寫出表示單元體內任意點的位移函數，其矩陣形式

  式3-29

式中：

2、對表達式（式4-8）應用節點處的邊界條件，寫出以｛α｝表示的節點位移u₁、v₁、w₁、u₂、v₂、w₂、。。。。。。，其矩陣形式

    式3-30

式中：

    式3-31

    式3-32

式中：是形狀函數，是單元體內任意點位移的關系式。

3、用表達式（式4-11）計算單元的應變

    式3-33

4、根據應力應變關系，計算單元的應力，

    式3-34

式中：是彈性矩陣，根據單元體的性質選取。

5、作用在單元上的外力對節點位移的外力位能

    式3-35

6、根據虛位移原理（或應變能），可以得到單元剛陣

    式3-36

7、由總位能泛函式的極值條件δΠ=0（或虛位移原理），可以得到力剛陣

    式3-37

8、把各單元進行組集，即

    式3-38

9、對K進行邊界條件處理，并求解方程組,得出結構各節點位移。

10、根據節點位移，應用和，計算各節點的應力。

題1：題圖1-a為一個對角受壓的正方形薄板，載荷沿厚度t均勻分布，其值為20MPa,為簡單起見取μ=0，t=1 m，求變形和應力。

解：由于對稱，可取板的四分之一來分析，如題圖1-b所示。

1、單元剖分

把計算部分（見題圖1-b）分為四個單元①、②、③、④，并對四個單元分別編出節點號碼i、j、k，如題圖1-c和題圖1-d。

2、單元分析

由于單元①與單元②、單元④具有相同的性質。所以只需單元①和單元③進行分析。

對于單元①，由于x_i=1，y_i=1；x_j=0，y_j=2；x_k=0，y_k=1，則

而

對于單元③，由于x_i=0，y_i=1；x_j=1，y_j=0；x_k=1，y_k=1，則

而

根據上列數值及μ=0，t=1，可得各單元的剛陣（此時）

3、單元組集

按節點位移序號組成全結構的總剛陣

4、邊界條件約束處理

由于u₁ = u₂ = u₄ = 0,v₄ = v₅ = v₆ =0，所以只需考慮v₁、v₂、v₃、u₃、u₅、u₆六個位移，因此縮減的剛陣是

節點力列陣是

與縮減的剛陣對應的縮減節點力列陣是

5、線性方程組建立與求解

將Fr、Kr之值代入Fr=Kr·qr中，得

解方程組，可得

6、單元應力分量的計算

由題可知，正方形薄板屬于平面應力問題。各單元①、②、③、④的應力矩陣分別為

則各單元的應力是

單元①：

單元②：

單元③：

單元④：

題2：題圖2-a為一個懸臂深梁，在右端作用著均布的剪力，其合力為P,采用圖示的簡單網格。設泊松比μ=1/3，t=1，試求結點位移量。

解：

1、單元剖分

把懸臂深梁剖分二個單元①，②，并分別編出單元節點號碼i,j。如題圖2（b）所示。

2、單元分析

對于單元①，由于x_i=2，y_i=0，x_j=0，y_j=1，x_k=0，y_k=0，則

而

對于單元②，由于x_i=0，y_i=1，x_j=2，y_j=0，x_k=2，y_k=1，則

而

根據上列數值及μ=1/3，可得各單元的剛陣（此時）

3、單元組集

按節點位移序號組成全結構的總剛陣為

4、邊界條件約束處理

由于u₁=v₁=u₄=v₄=0，所以只需考慮u₂,v₂,u₃,v₃四個位移。因此，減縮剛陣為

節點力列陣是

與縮減的剛陣對應的縮減節點力列陣是

5、線性方程組建立與求解

將Fr、Kr之值代入Fr=Kr·qr中，得

解方程組，可得