Bài giảng Phương pháp số (Phương pháp phần tử hữu hạn) - Vũ Khắc Bảy (ĐH Lâm nghiệp)

Tr­êng ®¹i häc l©m nghiÖp

Bé m«n To¸n

Vò Kh¾c B¶y

Bµi gi¶ng

ph­¬ng ph¸p sè

(ph­¬ng ph¸p phÇn tö h÷u h¹n)

Hµ néi - N¨m 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Để giải và tính toán các bài toán về kêt cấu cơ học, ngoài các phương pháp giải tích

LỜI NÓI ĐẦU

ta còn có các phương pháp số. Do các bài toán cơ học thường dẫn đến việc giải các

phương trình vi phân với các điều kiện biên xác định nào đó. Vì vậy thời kỳ đầu của các

phương pháp số là : các phương pháp tích phân số và phương pháp sai phân hữu hạn. Cùng

với sự phát triển của máy tính điện tử, phương pháp phần tử hữu hạn ra đời và phát triển

rất mạnh mẽ và là một phương pháp được dùng rất phổ biến hiện nay khi tính toán các

bài toán cơ học. Nó cũng đã được áp dụng để có được nhiều chương trình tính cho các

dạng bài toán cơ học khác nhau: Tính cho dàn thanh, khung không gian, các kết cấu

dạng tấm , vỏ ,...

Phương pháp phần tử hữu hạn là môn học cơ sở của các ngành kỹ thuật liên quan

đến tính toán các kết cấu và hiện cũng là một môn học của ngành Xây dựng và Kỹ thuật

công trình thuộc trường ĐHLN.

Trong năm trước đây chúng tôi có biên soạn nội dung bài giảng : Phương pháp

phần tử hữu hạn để phục vụ cho công tác giảng dạy môn học : Phương pháp số. Vẫn biết

rằng tài liệu viết về môn học này đã có rất nhiều trên các dạng : sách , bài giảng và trên

mạng, song thiết nghĩ thì việc biên soạn một tài liệu dạng bài giảng về phương pháp

phần tử hữu hạn với thời lượng 2 tín chỉ cũng là điều cần thiết để các em sinh viên ( và

cả các độc giả lần đầu biết về phương pháp này) tiếp cận với môn học này thuận lợi hơn.

Tài liệu mới chỉ tiếp cận đến một số nội dung và khái niệm cơ bản của phương

pháp phần tử hữu hạn. Các vấn đề trình bày mới dừng đến việc tính toán cho dàn, khung

không gian.

Tài liệu cũng đã đưa ra một số thủ tục cơ bản trong lập trình tính toán, các thủ tục

này được viết trong Visual Basic, độc giả có thể chuyển đổi dễ dàng sang các môi trường

lập trình khác.

Mong rằng với ý muốn như thế sẽ giúp ích được phần nào cho quá trình học tập

môn học này của các em sinh viên, và tất nhiên rất mong được các đóng góp của độc giả

về các vấn đề trình bài trong tài liệu.

Tác giả

1 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Chương I MỘT SỐ VẤN ĐỀ CƠ BẢN TRONG CƠ HỌC VẬT RẮN BIẾN DẠNG & CÁC PHƯƠNG PHÁP TÍNH TRONG CƠ HỌC

I.1 CÁC VẤN ĐỀ CƠ BẢN TRONG CƠ HỌC VẬT RẮN BIẾN DẠNG I.1.1 Ten xơ ứng suất.

Dưới tác dụng của lực ngoài, vật thể chịu lực bị biến dạng và bên trong nó sẽ

xuất hiện ứng suất. Ứng suất tại mỗi điểm khác nhau là khác nhau, véc tơ ứng suất không

những phụ thuộc vào điểm mà còn phụ thuộc vào hướng của thiết diện qua nó mà được

 nT

. Như vậy tập hợp cặp véc tơ ứng suất và véc tơ

 xác định bởi pháp tuyến có hướng n  n tại điểm P sẽ xác định trạng thái ứng suất tại điểm đó. Trạng thái ứng suất tại điểm

hoàn toàn được xác định qua ten-xơ ứng suất – là một ten xơ đối xứng hạng hai, nên nó

σ

σ

σ

11

12

13



σ

có 6 thành phần độc lập:

σ

với σ

ji

ij

ij

21

22

23

σ σ

σ σ

σ σ

31

32

33

     

    

; τ

; τ

Trong hệ tọa độ De-cac các thành phân của ten xơ ứng suất được ký hiệu là :

σ ;σ ; σ ; τ y

x

z

xy

xz

yz

I.1.2 Phương trình cân bằng

Tách phần thể tích V tùy ý giới hạn bởi mặt S của môi trường liên tục ở hình thái

biến dạng, xét sự cân bằng các lực tác dụng lên thể tích đó ( không kể lực quán tính) ta

được :

 dVK

 dVK



0



0

 dST n

 dSnT i i



 



 

S

S

V

V

hay là

dV



 dSnT i i





 T  i x 

V

S

i

( công thức Gaoxơ - Ôtrôgratxki) nên ta có : Do





 T  i x 

V

i

  

   dVK  

  ij

, vì V là thể tích tùy ý nên biểu thức dưới dấu tích phân bằng không

hay

lµ

 K 





0



K 



0

j

x

 T  i x 



i

i

=> ta được : (I.1)

2 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

 Phương trình (I.1) được gọi là phương trình cân bằng, ρ là mật độ khối lượng , K

lực

khối.

 τ







ρ

Viết (I.1) dưới dạng tường minh ta được :

K



0

x

xy  y  σ

 τ xz  z  τ







ρ

(I.2)

K



0

y

yz  z

y  y  τ







ρ

K



0

z

τ  zx  x

zy  y

σ  z  z

  σ x   x   τ  yx   x    

 τ





ρ

Với bài toán hai chiều ( Tấm , vỏ ), phương trình cân bằng có dạng :

K



0

x

(I.3)

xy  y  σ





ρ

K



0

y

y y



 σ x   x    τ  yx   x 

Còn trong bài toán một chiều phương trình cân bằng sẽ là :





0

(I.4)

xK ρ

 σ x  x

I.1.3 Quan hệ giữa biến dạng – chuyển vị ( hệ thức Cô-si)

 đổi kích thước dài đoạn vô cùng nhỏ bất kỳ Xd

 lấy từ điểm X

Khi xây dựng hệ thức quan hệ giữa biến dạng và chuyển vị, xuất phát từ sự thể hiện thay

và thay đổi góc giữa hai

đoạn vô cùng nhỏ bất kỳ lấy từ điểm đó người ta đã dẫn đến ten-xơ biến dạng hữu hạn

 U

j

i

k

viết trong hệ tọa độ Đề các : Grin và Anmăngxi

γ







ij

 X

j

i

U U   k   X X i

j

 1 U    2 X  

   

 u

j

i

k

(Grin) (I.5)

 γ







ij

1 2

 u  x

 x

 u k  x

u   x

j

i

i

j

   

   

(Anmăngxi) (I.6)

Trong đó Xn - biến Lagrăng, xk – biến Ơ le , Um và un là các thành phần chuyển vị theo

biến Lagrăng và Ơle.

3 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Trường hợp biến dạng nhỏ ( t.là khi bỏ qua VCB bậc cao – là các thành phần phi tuyến

trong (1.5) và (1.6) )) khi đó hai ten xơ này xấp xỉ bằng nhau và các thành phần của ten-

11

3

2

xơ biến dạng nhỏ sẽ có dạng:



, ε



, ε



,



ij

21

ε 11

33

22

 

u x

 u  x

 u 1  x 1

3

2

  12 13  23 22   33 32

31

    

    

2

3

1

3





, ε





, ε





ε 12

23

13

1 2

 u 1  x

1 2

u   x

 u  x

1 2

 u  x

u   x

2

u  2  x 1

3

2

3

1

  

  

  

  

  

  

với các thành phần :

Để thuận lợi cho các công thức sau này trong tính toán theo phương pháp PTHH, người

ta ký hiệu :

- Các thành phần chuyển vị : u , v , w

ε



; ε



; ε



;

x

y

z

u   x

v   y

w   z

- Các thành phần của ten-xơ biến dạng ε x ; ε y ; ε z ; γ xy ; γ xz ; γ yz với :

γ





; γ





; γ





xy

xz

yz

u   y

 v  x

u   z

 w  x

v   z

w   y

(I.7)

0

0



x

0

 y

y



0

0

 z

z

Hay viết dưới dạng ma trận :





xy

0

 y

 x

u   . v   w 

    

ε ε ε γ γ

yz





0

 z



y

γ

xz

        

        





0

  x  0            

 z

 x

             

(I.8)

I.1.4 Phương trình liên tục

Hệ thức Cô-si (I.8) cho sự liên hệ giữa 6 thành phần biến dạng xác định duy nhất qua 3

thành phần chuyển vị cho trước. Như vậy với 6 thành phần biến dạng cho trước chỉ từ

quan hệ (I.8) sẽ không cho duy nhất 3 thành phần chuyển vị, do đó giữa các thành phần

4 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

biến dạng sẽ có 6 phương trình tương thích biến dạng ( còn gọi là các phương trình tương

thích biến dạng Xanhvơnăng), các phương trình này đảm bảo cho sự biến dạng liên tục

2

2

2

2

 γ

 γ



γ









2





trong môi trường.

ε y 2

  x

 γ zx  y

yx  z

zy  x

 ε x   y z

xy   x y

 x

  

  

2

 ε x 2  y 2

2

2

 γ



γ

ε





γ











2





(I.9)

ε y 2

ε z 2

  y

zy  x

yx  z

 γ zx  y

y   x z

yz   z y

 z

  

  

 y 2

2

2

2

 γ

 γ













2

ε z 2

γ  xz   x z

  z

 γ zx  y

yz  x

xy  z

ε  z   y x

 

ε x 2 z

 x

  

  

2

2

2



γ







 Trong bài toán 2 chiều : (I.9) còn 1 phương trình

ε y 2

xy   x y

 ε x 2  y

 x

 Trong bài toán 1 chiều : các phương trình trên đều thỏa mãn.

I.1.5 Điều kiện biên

Điều kiện biên liên quan đến chuyển vị hoặc ứng suất

 Điều kiện biên liên quan đến chuyển vị ( gọi là điều kiện biên động học) thường

được cho trước chuyển vị của một điểm hoặc một phần mặt biên nào đó

 Điều kiện biên liên quan đến ứng suất ( gọi là điều kiện biên tĩnh học) đòi hỏi sự

cân bằng giữa ứng suất trên mặt biên với ngoại lực đặt lên đó.

Ví dụ

Một thanh chiều dài  , chiều dầy h,

bị ngàm chặt một đầu, một đầu tự

do, chịu tác dụng của lực phân bố

đều có cường độ q như hình vẽ.

Chọn hệ tọa độ : 0x theo chiều dài,

sát mặt dưới, 0y hướng lên trên

Hình 1.1

Đây là bài toán hai chiều và các điều kiện biên được đưa về các hệ thức sau :



0

 v(0, y)  x

- Tại x = 0 : u(0,y) = v(0,y) = 0 ;

5 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

τ

(x, h)



  q

0 ; σ (x, h) y

τ

(x,0)





0

- Tại mặt trên ( y = h) : xy

0 ; σ (x, 0) y



τ

 ( , y)





0

- Tại mặt dưới ( y = 0) : xy

0 ; σ ( , y) x

- Tại đầu B ( x =  ) : xy

I.1.6 Phương trình vật lý (phương trình trạng thái)–Quan hệ ứng suất và biến dạng.

Trong giáo trình này chỉ xét giai đoạn làm việc của vật liệu ở giai đoạn đàn hồi,

biến dạng là nhỏ và đàn hồi là tuyến tính. Như vậy quan hệ giữa ứng suất và biến dạng ở

đây được áp dụng bởi định luật Hooke.

Xét với vật liệu đẳng hướng :

ε



σ



ν(σ



γ



τ



τ

x

x

y

xy

xy

xy

 

 σ ) ,  z

I.1.6.1 Bài toán 3 chiều : định luật Hooke có dạng :

ε



σ



ν(σ



γ



τ



τ

y

y

x

yz

yz

yz

 

 σ ) ,  z

ε



σ



ν(σ



γ



τ



τ

z

z

x

xz

xz

xz

 

 σ ) ,  y

1 G 1 G 1 G

 2(1 ν) E  2(1 ν) E  2(1 ν) E

1 E 1 E 1 E

(I.10)





Chú ý rằng các thành phần của ten-xơ biến dạng được tính theo chuyển vị qua (I.7)

  ε

   C . σ



(I.11) Nếu viết dưới dạng ma trận và có kể đến biến dạng ban đầu thì (I.10) có dạng:  0 ε

trong đó:



  ε

xy

yz

y

z

zx

 ε ,ε ,ε , γ , γ , γ x

T

- là véc tơ biến dạng



, τ

, τ

x

z

yz

zx

xy

 σ ,σ ,σ , τ y

T

- là véc tơ ứng suất   σ



,ε

,ε

, γ

, γ

, γ

  ε

0

0x

0y

0z

0xy

0yz

0zx

 ε

T

- là véc tơ biến dạng ban đầu

( Chữ T – ký hiệu chuyển vị ma trận)

[C] – ma trận các hệ số đàn hồi ,

6 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .



ν



ν

0

0

0

ν

ν

ν ν

C







1 E

1  0 0 0

 1 0 0 0

0 0  2(1 ν) 0 0

0 0 0  2(1 ν) 0

0 0 0 0  2(1 ν)

  0 0 0

        

1         

(I.12)

0

với E – mô đun đàn hồi Young , G – mô đun trượt , ν - hệ số Poát-xông của vật liệu.



αT 1 , 1 , 1, 0, 0, 0



T

σ



, trong Trường hợp biến dạng ban đầu do nhiệt độ thì   0ε

    D ε

0   ε

đó α - hệ số dãn nở vì nhiệt, T0 – độ biến thiên của nhiệt độ. Biểu diễn ứng suất qua các thành phần biến dạng ta sẽ có :    

0



    σ

  D ε

E αT  1 2ν

1     1     1   0     0   0  

hay là :

ν  1 ν ν

ν ν  1 ν

0 0 0

0 0 0

0

0

0

0

0

với ma trận hệ số D là :

D







0 0 0  1 2ν 2

E (1 ν).(1 2ν)





0

0

0

0

0

 1 2ν 2

0

0

0

0

0

 1 ν  ν   ν         

 1 2ν 2

            

(I.13)

I.1.6.2 Bài toán 2 chiều :

 Bài toán ứng suất phẳng : ví dụ như xét các bài toán về tấm, vỏ với tải

trọng nằm trong mặt phẳng giữa tấm, phân bố đều theo bề dầy của tấm khi đó chọn trục z

vuông góc với mặt phẳng tấm, sẽ dẫn đến thể giả thiết rằng :

7 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

σ



τ



τ

 0

z

xz

yz

ε





  C σ

 0 ε

 ứng suất không đổi theo chiều dầy của tấm. với giả thiết này các biểu thức của định luật Hooke có dạng :    

x

x

0x

0











  ε

  ; σ

  ; ε





0

y

y

0y

1 E

ε ε γ

σ σ τ

ε ε γ

1  ν 0

 ν 1 0

xy

xy

0xy

1     αT 1 ; C     0  

    

0   0   2(1 ν 

    

    

    

    

    

    

0



ε





với :

 σ D ε



  D ε

0

      



E αT  1 ν

1     1      0

D



1 ν

ν 1

Hay biểu diễn ngược lại:  



2

E  1 ν

0

0

      

  0  0   1 ν   2

0

với 

ε



σ



σ



αT

z

x

y





 ν E

biến dạng theo phương z vẫn tồn tại và

 Bài toán biến dạng phẳng : Khi xét vật thể hình lăng trụ dài có mặt cắt

ngang không đổi theo chiều dài ( theo chiều trục 0z) , chịu tải trọng đều



0

ε



z

 w  z

vuông góc với 0z , khi đó ta có : w = 0 ; ; các đại lượng





  ε



   C . σ

 0 ε

0x









0  1 ν αT 1 ; C



  ε





0

0y

 1 ν E

 ν  1 ν 0

0 0 2

ε ε γ

0xy

1         0  

1 ν     ν   0 

    

    

    

0



ε





ứng suất và biến dạng chỉ phụ thuộc vào các biến x và y.

 σ D ε



  D ε

0

      



E αT  1 2ν

1     1      0

Hay biểu diễn ngược lại:  

8 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

D



ν  1 ν



E (1 ν).(1 2ν)





0

0

  1 ν   ν    

      

với 



σ



σ



0 E αT ; τ



τ



0

z

y

x

yz

xz

0 0  1 2ν 2  ν σ



ứng suất theo phương z vẫn tồn tại và

0

0

I.1.6.3 Bài toán 1 chiều :

ε



σ



αT

σ



Eε



EεT

x

x

x

x

1 E

hoặc ( => D = E – mô đun đàn hồi)

I.1.7 Đặt bài toán đàn hồi : Thiết lập bài toán đàn hồi bao gồm việc thiết lập các

phương trình và các điều kiện biên, chúng phải lập thành một hệ kín để có thể giải ra

được các ẩn cần tìm đó là các giá trị của các thành phần ten – xơ biến dạng, ứng suất ,

véc tơ chuyển vị. Các phương trình gồm có :

 Phương trình cân bằng

 Hệ thức Cô-si ( liên hệ chuyển vị và biến dạng)

 Phương trình trạng thái ( liên hệ ứng suất và biến dạng : định luật Hooke)

và các điều kiện biên động học hoặc tĩnh học.

Người ta đã chứng minh được sự tồn tại và duy nhất nghiệm của bài toán đàn hồi.

I.2 CÁC PHƯƠNG PHÁP GIẢI BÀI TOÁN ĐÀN HỒI

Như phần trên đã trình bày, ta thấy rằng để giải bài toán đàn hồi tuyến tính : 3

chiều ta có tới 15 phương trình cùng các điều kiện biên để tìm ra giá trị của 15 ẩn : 3

thành phần chuyển vị, 6 thành phần biến dạng và 6 thành phần ứng suất. Cho đến nay

cũng đã có được các phương pháp giải đúng và gần đúng. Có những phương pháp có thể

áp dụng tốt cho một lớp các dạng bài toán cơ học biến dạng nhưng lại khó khăn áp dụng

nó cho dạng khác. Có thể tổng kết ra đây theo sơ đồ sau

9 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Các phương pháp giải

Các phương pháp số

Các phương pháp giải tích

Các phương pháp chính xác

Các phương pháp gần đúng ( biến phân)

Các phương pháp số giải các phương trình vi phân

Phương pháp phần tử hữu hạn ( PTHH )

Mô hình cân bằng

Mô hình hỗn hợp

Mô hình tương thích

Các phương pháp tích phân số

Phương pháp sai phân hữu hạn

I.2.1 Phương pháp chính xác

Phương pháp giải tích giải các bài toán đàn hồi : Có thể giải theo chuyển vị , hoặc

giải theo ứng suất.

ví dụ : xét bài toán uốn dầm như hình 1.2 : dầm chịu tải trọng đều với cường độ q , tựa

gối khớp tại hai đầu A , B

Lời giải bài toán trên đã có trong sức bền 0  x

EJ

 

q 0 ; ( 0 x

   (1.14) )

4 d w 4 dx

vật liệu :

z





với điều kiện biên : hình 1.2

w(0) w( ) 0 

 



w (0) w ( ) 0





  

4

2

Tích phân liên tiếp (1.14) ta được chuyển vị theo

w(x)





3  C x C x



 C x C

2

1

3

4

 x q  EJ 24 

  

phương z : , từ các điều kiện biên ta xác định

10 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

3

 

; C



C 1

3

 12

 24

4

3



. dẫn đến nghiệm bài toán được các hằng số tích phân : C2= C4 = 0 ;

w(x)



x



3  2 x

x





q 24EJ

2

2

theo chuyển vị => tìm được mô men uốn:

M(x)

 

EJ



x

 

x

/ 2





45q  384EJ

d w q 2 2 dx

; độ võng lớn nhất tại x= là wmax =

I.2.2 Các phương pháp biến phân :

Các bài toán cơ học nói chung có thể biểu diễn dưới dạng tổng quát :

L(u)

  g

0

C(u)



p

trong miền V

  

trên biên S (1.15)

trong đó L , C là các toán tử vi phân , u = u (x,y,z) đại lượng cần tìm ,

g = g(x,y,z) và p=p(x,y,z) là các hàm cho trước.

Phương pháp biến phân là phương pháp gần đúng nhằm tìm một nghiệm xấp xỉ dựa trên

một tiêu chuẩn nào đó. Các phương pháp biến phân thường gặp đó là phương pháp Ritz,

phương pháp Galerkin, ...

1. Phương pháp Ritz :

Trong một số các bài toán đàn hồi thường tồn tại một phiếm hàm I dạng :

I

 F(x, y, u, u , u , u , u ..)dxdy xx

 yy

 y

 x

 

V

và từ điều kiện dừng của phiếm hàm này sẽ dẫn ra

các phương trình vi phân của bài toán.

Phương pháp Ritz tìm nghiệm xấp xỉ gần đúng trong một dạng tổ hợp tuyến tính các hàm

N

biết trước, chẳng hạn như trong bài toán phẳng , theo Ritz thì chuyển vị có dạng :

 u(x, y) φ (x, y)

0

C φ (x, y) i

i

 

 i 1

(1.16)

trong đó :

 φ 0(x,y) là hàm được chọn sao cho thỏa mãn các điều kiện biên không thuần

nhất.

 φ i(x,y) là các hàm thỏa mãn :

 Khả vi đến cấp cần thiết

11 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

 Thỏa mãn các điều kiện biên thuần nhất

i

là độc lập tuyến tính và đầy đủ  Hệ hàm  N iφ

Với các điều kiện trên sẽ làm cho nghiệm của Ritz hội tụ đến nghiệm chính

xác ( khi N   )

Ci là các hằng số cần tìm, chúng được xác định khi thay u(x,y) vào phiếm hàm I,

lấy tích phân trên V, khi đó I sẽ là hàm của các Ci, cho thỏa mãn điều kiện dừng của I

(tức là các đạo hàm riêng của I theo các Ci bằng 0) sẽ dẫn đến hệ phương trình đại số xác

định các hằng số Ci.

2. Phương pháp phần dư có trọng (Weighted Residual Method)

Gọi R(u) = L(u) + g được gọi là hàm phần dư, như vậy nghiệm của bài toán

N

L(u) + g = 0 sẽ được chuyển thành tìm u là nghiệm của bài toán R(u) = 0. Người ta cũng

u



φ

0

N

C φ i

i

 

 i 1

tìm nghiệm dưới dạng , trong đó :

 φ 0 : được chọn sao cho thỏa mãn các điều khiện biên không thuần nhất

 φ i : khả vi đến số lần cần thiết , thoả mãn các điều kiện biên thuần nhất và

hệ hàm { φ i} là độc lập tuyến tính.

khi đó R(uN) = R(x,y,z,C1, C2, ...,CN).

Các hằng số Ci được tìm bằng phương pháp : lấy tích phân trên V của tích các hàm ψ k

với R(uN) , tức là :

ψ R(u ) dv



 0 ; k 1, 2,..., N

k

N



V

(1.17)

các hàm ψ k được gọi là các hàm trọng số (Weighted function)

Chú ý rằng toán tử vi phân L trong (1.15) là tuyến tính nên thay R(uN) vào (1.17)

N

ta sẽ dẫn đến :

ψ L(φ )dv





0

ψ R(u ) dv



C



 g dv

k

i

 ψ L(φ ) 0

k

N

k

i









 i 1

V

V

V

  

  

; k =1,2,...,N (1.18)

ψ L(φ )dv



 g dv

k

i

 ψ L(φ ) 0

k





V

V

Nếu đặt = Aki ; = Bk khi đo ta sẽ có hệ

phương trình đại số tuyến tính : {A}{C} = {B} xác định N hằng số Ci , thay vào uN ta

nhận được nghiệm gần đúng. 12 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

ψ k  φ k . Như vậy nếu phương pháp Galerkin chọn hệ hàm { φ k} là hệ trực giao thì sẽ

3. Phương pháp Galerkin : Phương pháp Galerkin là phương pháp dư có trọng lấy

rất thuận lợi dẫn đến (1.18) .

EJ

 

q 0 ; ( 0 x

   )

4 d w 4 dx

Ví dụ: Trong ví dụ trên phần 1.2.1 ta co phương trình vi phân của độ võng

N

N

nghiệm gần đúng theo phương pháp Galerkin sẽ được chọn :

w





x

N

C w (x) i

i

C sin i





πi 

 i 1

 i 1





,

w(0) w( ) 0 

 



w (0) w ( ) 0





  

N

như vậy wN thỏa mãn các điều kiện biên :

x





0

C sin i



πi 

q EJ

4 4 π i 4 

 i 1

Thay vào phương trình vi phân của độ võng ta được : ,

sin

x

πk 





N

x sin



sin

xdx

C sin i

 



πk 

πk 

πi 

q EJ

4 4 π i 4 

 i 1 0

0

  

 x dx  

Nhân cả hai vế với và lấy tích phân hai vế ( k = 1, 2,3,..., N) dẫn đến :

 sin  

πk   

chú ý rằng hệ là một hệ trực giao trong không gian có tích vô hướng là



k



tích phân xác định trên [0 ,  ] nên ta luôn có :

sin

x sin

x dx



πi 

πk 

khi i



k

0

    

0 khi i  2





sin

xdx



;

   1 ( 1)

k



 khi k 2m 1



πk 

 kπ

0

0 khi k 2 m  2  (2m 1)π

     

4

Bk =

4 π k 3 2

; như vậy ta có ma trận (A) là ma trận chéo với A kk =

 2 q  (2m 1)πEJ

B2m = 0 ; B2m -1 =

13 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

3

 2



4

4

5

 2 q 

(2m 1)πEJ π (2m 1)



4  4 q 5  (2m 1) π EJ

N

M



C2m = 0 ; C2m - 1 =

w





sin

x

N

C w (x) i

i

5





π(2m 1) 

4  4 q 5  (2m 1) π EJ

 i 1

 m 1

vậy nghiệm gần đúng

/ 2

4

4

M







w



sin







max

5



π(2m 1) 2

0,013054 q EJ

4  5 q 384 EJ

0,013021 q EJ

4  4 q 5  (2m 1) π EJ

 m 1

Độ võng lớn nhất tại x = => Nếu lấy M =10 ta được

sai số 0,25%

4. Phương pháp sai phân : Phương pháp sai phân là phương pháp biểu diễn gần đúng

các giá trị của đạo hàm theo các giá trị của hàm số tại các điểm lân cận ( xuất phát từ khai

 f (x Δx)

 f (x Δx)



df dx

 2Δx

x

 f (x Δx)







 2

2 d f 2 dx

f (x Δx) 2f (x) Δx

x

f (x 2Δx) 3f (x Δx) 3f (x)









 f (x Δx)



 3

3 d f 3 dx

Δx

x

f (x 2Δx) 4f (x Δx) 6f (x) 4f (x Δx)















 f (x 2Δx)



4

4 d f 4 dx

Δx

x

triển Tay-lo của hàm số). Chẳng hạn đối với hàm một biến f(x) ta có được :

khi đó việc giải các phương trình vi phân được đưa về việc tìm các giá trị hàm số tại các

điểm nút lưới ( mà khoảng cách giữa các điểm chính là Δx ).

ví dụ : xét bài toán uốn dầm như hình 1.2 : dầm chịu tải trọng đều với cường độ q , tựa

gối khớp tại hai đầu A , B

EJ

 

q 0 ; ( 0 x

  

)

4 d w 4 dx

Có phương trình vi phân của độ võng 0  x





w(0) w( ) 0 

 



w (0) w ( ) 0





  

z với điều kiện biên :

14 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Chia lưới sai phân như hình vẽ : 4 đoạn với 5 nút (từ 0 đến 4, với nút 0 tại x = 0 ; nút 4



 4

2

w

2w



 w 0.Δx

0

tại x =  ), có Δx

 => w-1 = - w1

  1

0

1

4

tại nút 0 : w0 = 0

w

4w +6w 4w





w



  1

1

0

3

2

pΔx EJ

4

tại nút 1:

w 4w +6w









0

1

4w w 3

4

2

4

tại nút 2:

w 4w +6w 4w







w



2

3

1

5

4

pΔx EJ pΔx EJ

2

tại nút 3:

 w 2w



w



0.Δx

 0

3

4

5

4

4

tại nút 4: w4 = 0



pΔx EJ

 0,013672p EJ

sai số tương đối : 0.049987 => w2 =3,5.

I.3 CÁC NGUYÊN LÝ BIẾN PHÂN

I.3.1 Nguyên lý thế năng toàn phần dừng ( nguyên lý biến phân về chuyển vị)

Thế năng toàn phần của một hệ đàn hồi  được xác định là :  = U - A trong đó :

U - thế năng biến dạng của vật thể đàn hồi được tích lũy trong quá trình biến dạng

A – công của ngoại lực sinh ra trên các chuyển dời của ngoại lực do vật thể bị biến

dạng.

Nguyên lý phát biểu rằng : Trong tất cả các trường chuyển vị ( các trạng thái chuyển

vị ) khả dĩ động ( tức là thỏa mãn điều kiện tương thích và điều kiên biên động học) thì

trường chuyển vị thực ( tương ứng với sự cân bằng của vật thể) sẽ làm cho thế năng

toàn phần  đạt giá trị dừng (nhỏ nhất).

Tức là khi đó: δ  ({u}) = δ U({u}) - δ A({u}) = 0.

T

T

U



{ε} {σ}dv hay U



{ε} [D]{ε}dv





1 2

1 2

V

V

Thế năng biến dạng U được tính theo công thức :

15 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

T

T

U



{ε} [D]{ε}dv



{ε} [D]{ε }dv

0





1 2

V

V

nếu có biến dạng đầu thì :

T

T

A



{u} {g}dv



{u} {p}ds





V

S

t

Còn công của ngoại lực ( gồm lực khối {g} , lực mặt {p}) trên chyển dời {u} sẽ là :

T

T

T

khi đó thế năng toàn phần sẽ là :



 ({u}) U A





{ε} [D] {ε}





{u} {g}dv



{u} {p}ds



 2{ε } dv 0







1 2

V

V

S

t

(I.19)

Ví dụ :

Xét dầm có chiều dài  , tựa khớp hai đầu và chịu tải phân bố đều q,



0  x Thế năng biến dạng

M(x)

 

EJ

dx



2 1 M (x) 2

EJ

2 d w 2 dx

0

2

U = mà





z

U



dx

q.w dx





EJ 2

2 d w 2 dx

0

0

  

  

2





=> ; A =

 

dx



q.w dx





EJ 2

2 d w 2 dx

0

0

  

  

khi đó thế năng toàn phần : .

 

 δU δA

 , tức

0

Theo nguyên lý thế năng toàn phần dừng : w(x) là thực nếu : δ



2

2

2

2

2

3







 δU EJ

dx



EJ

dx



EJ



EJ

dx

2

2







dx

d w d δw 2 dx

dx

d w dδw 2 dx dx

dδw d w . 3 dx dx

0

0

0

 d w d w δ  2 dx 

  

  

  

0



3

2



là:

δU



EJ



δw



δw.

dx

2



d w dδw d w 3 dx dx dx

4 d w 4 dx

0

  

  

0

    

    



3

2





=>

δ

 

EJ.



q δw dx EJ





δw

δA

q.δw dx

2



 

4 d w 4 dx

d w dδw d w 3 dx dx dx

0

0

  

  

  

  

0





còn ; vậy

w(0) w( ) 0 

 



w (0) w ( ) 0





  

do điều kiện biên :

16 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .



δ

 

EJ.



q δw dx



2 d w 2 dx

4 d w 4 dx

0

  

  

EJ.

 = 0 , đây chính là phương trình

q

do đó tại x = 0 và x =  ta có : = 0 ; δw =0 =>

4 d w 4 dx

để thỏa mãn với mọi δw thì cần phải có

cân bằng của dầm chịu uốn viết cho chuyển vị.

I.3.2 Nguyên lý cực tiểu của năng lượng bù toàn phần ( nguyên lý biến phân theo

ứng suất)

* được định nghĩa là :

* = U* - A* , trong đó U* là năng lượng bù của biến dạng.

T

Năng lượng bù toàn phần của vật thể đàn hồi

{σ} {ε}dv

1 2 

V

*

T

Trong giai đoạn đàn hồi tuyến tính thì U* = U = , nếu có biến

U



{σ} [C]{σ} 2{ε } dv







0

0{ε } thì



1 2

V

T

dạng đầu

{p} {u}ds



S

t

T

T

A* là công bù của ngoại lực : A* = trong đó st là phần biên đã biết chuyển

*  

{σ} [C]{σ} 2{ε } dv





{p} {u}ds





0





1 2

S

V

t

vị {u}, vì vậy

Nguyên lý cực trị năng lượng bù toàn phần được phát biểu như sau:

Trong tất cả các trường ứng suất khả dĩ tĩnh ( tức là thỏa mãn các điều kiện cân

bằng và điều kiện biên tĩnh học trên St ) thì trường ứng suất thực ( tương ứng thỏa mãn

* đạt giá trị dừng.

*

δU



* δU ({σ})



* δA ({σ})

 0

điều kiện tương thích) sẽ làm cho năng lượng bù toàn phần

17 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Chương II

CƠ SỞ VÀ CÁC NỘI DUNG CƠ BẢN CỦA

PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN

II.1 Khái niệm về phương pháp phần tử hữu hạn (PTHH)

Phương pháp phần tử hữu hạn (PP PTHH) là một phương pháp số để tìm nghiệm

gần đúng của một hàm chưa biết trong miền xác định V. Tuy nhiên PP PTHH không tìm

dạng xấp xỉ của hàm cần tìm trên toàn miền V mà chỉ trong từng miền con Ve ( phần tử)

thuộc miền xác định V. Chính vì lẽ đó nên phương pháp này rất thích hợp để tìm nghiệm

gần đúng cho các bài toán vật lý, kỹ thuật khi mà hàm cần tìm được xác định trên những

miền phức tạp là những vùng nhỏ có các đặc trưng hình học, vật lý khác nhau, chịu các

điều kiện biên khác nhau. Phương pháp được phát biểu một cách tổng quát chặt chẽ như

một phương pháp biến phân hay phương pháp dư có trọng số trên mỗi phần tử.

Trong PP PTHH , miền V được chia thành một số hữu hạn các miền con được

gọi là các phần tử. Các phần tử này được kết nối với nhau tại các điểm trên biên được gọi

là các nút. Trong phạm vi mỗi phần tử, đại lượng cần tìm ( chẳng hạn đó là các biến

dạng, dịch chuyển, ứng suất ,…) được lấy xấp xỉ trong một dạng hàm đơn giản – được

gọi là các hàm xấp xỉ ( approximation function). Các hàm xấp xỉ này được được tính

thông qua các giá trị của nó ( đôi khi qua các giá trị đạo hàm) tại các điểm nút trên phần

tử và các giá trị này được gọi các bậc tự do của phần tử mà ta xem như là các ẩn cần tìm

của bài toán.

Trong bài toán cơ học vật rắn biến dạng và cơ kết cấu tùy theo ý nghĩa vật lý của

các hàm xấp xỉ ta có thể áp dụng bài toán theo ba loại mô hình sau:

1. Mô hình tương thích : Xem chuyển vị là đại lượng cần tìm trước và hàm xấp xỉ

biểu diễn gần đúng dạng phân bố của chuyển vị trong phần tử. Các ẩn số được xác

định từ hệ phương trình thiết lập trên cơ sở nguyên lý thế năng toàn phần ( hay

nguyên lý biến phân Lagrange).

2. Mô hình cân bằng : Hàm xấp xỉ biểu diễn gần đúng dạng phân bố của ứng suất

hay nội lực trong phần tử. Các ẩn số được xác định từ hệ phương trình thiết lập

18 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

trên cơ sở nguyên lý năng lượng hệ toàn phần dừng ( hay nguyên lý biến phân về

ứng suất – nguyên lý Castigliano)

3. Mô hình hỗn hợp : Xem các đại lượng chuyển vị và ứng suất là hai hai yếu tố

độc lập. Các hàm xấp xỉ biểu diễn gần đúng dạng phân bố của chuyển vị và ứng

suất trong phần tử. Các ẩn cần tìm được xác định từ hệ phương trình thiết lập trên

cơ sở nguyên lý biến phân Reisner.

Sau khi tìm được giá trị các ẩn số ( bằng việc giải một hệ phương trình đại số), như

vậy ta đã tìm được xấp xỉ các đại lượng cần tìm, từ đó tìm được giá trị của các đại

lượng còn lại.

Mô hình tương thích được áp dụng rộng rãi. Trong giáo trình này chủ yếu các bài

toán được giải theo mô hình tương thích.

II.2 Trình tự các bước phân tích bài toán theo phương pháp PTHH

 Bước 1. Rời rạc hóa miền khảo sát

Miền khảo sát V được chia thành các miền con Ve ( phần tử) có dạng hình học thích hợp.

Với bài toán cụ thể thì số phần tử, hình dạng hình học của phần tử và kích thước các

phần tử phải được xác định cụ thể. Số điểm nút mỗi phần tử không được lấy tùy tiện mà

phải phụ thuộc vào dạng hàm xấp xỉ định chọn.

Các phần tử có các dạng hình học đơn giản :

Phần tử một chiều :

Phần tử bậc 1 Phần tử bậc 2 Phần tử bậc 3

Phần tử hai chiều:

Phần tử bậc 1 Phần tử bậc 2 Phần tử bậc 3

19 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Phần tử ba chiều :

Dạng tứ diện

Phần tử bậc 1 Phần tử bậc 2 Phần tử bậc 3

Dạng lăng trụ

Phần tử bậc 1 Phần tử bậc 2 Phần tử bậc 3

 Bước 2 Chọn hàm xấp xỉ thích hợp : Chọn dạng hàm xấp xỉ sao cho đơn

giản đối với tính toán, nhưng vẫn đảm bảo các tiêu chuẩn hội tụ . Thường chọn các

hàm này có dạng đa thức.

sau khi chọn dạng hàm xấp xỉ ta biểu diễn các hàm này (kể cả đạo hàm của

nó) theo tập hợp các giá trị tại các nút của phần tử {q}e

 Bước 3 Xây dựng phương trình phần tử, tức là thiết lập ma trận độ cứng

phần tử [K]e và véc tơ tải phần tử {P}e .

Kết quả nhận được phương trình có dạng : [K]e {q}e = {P}e

 Bước 4 Ghép nối các phần tử trên cơ sở mô hình mà kết quả là hệ thống



{P}

[K] là ma trận độ cứng tổng thể ( toàn miền V)

{q} là véc tơ tập hợp các giá trị đại lượng cần tìm tại tất cả các nút ( tức là véc

trong đó : phương trình : [K] {q}

{P} là véc tơ số hạng tự do tổng thể ( tức là véc tơ tải tổng thể )

tơ chuyển vị nút tổng thể)

*

sau đó sử dụng điều kiện biên của bài toán sẽ nhận được hệ phương trình :

* [K ] {q }



* {P }

- là hệ phương trình hệ thống hay còn gọi là hệ phương trình để giải

20 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

*

* [K ] {q }



* {P }

*

 Bước 5 Giải hệ phương trình đại số : , tìm được

* [K ] {q }



* {P }

chuyển vị của các nút. Việc giải hệ phương trình đối với bài toán

tuyến tính không gặp khó khăn, nhưng với bài toán phi tuyến thì sẽ dùng phương pháp

lặp ( mà được tuyến tính hóa , chẳng hạn như phương pháp Newton – Raphson) mà ở

* [K ] và

*{P } sẽ thay đổi.

mỗi bước lặp ma trận độ cứng

 Bước 6 Hoàn thiện: tính giá trị các đại lượng còn lại: ứng suất, biến dạng

II.3 Hàm xấp xỉ - đa thức xấp xỉ . Phép nội suy

1. Hàm xấp xỉ :

Tư tưởng chính của PP PTHH là xấp xỉ hóa đại lượng cần tìm trong mỗi miền con

– phần tử Ve . Do đó đầu tiên phải chọn hàm số đơn giản mô tả gần đúng đại lượng cần

tìm trong mỗi phần tử. Hàm số đơn giản hay được chọn có dạng đa thức, vì:

 Đa thức được xem như tổ hợp tuyến tính các đơn thức , các đơn thức này

thỏa mãn yêu cầu của Ritz , Galerkin.

 Hàm xấp xỉ dạng đa thức thường dễ tính toán, dễ thiết lập công thức khi

xây dựng các phương trình PP PTHH, dễ đạo hàm, dễ lấy tích phân.

 Có khả năng tăng độ chính xác ( bằng cách tăng số bậc của đa thức), tuy

nhiên trong thực tế thường chỉ lấy bậc thấp mà thôi.

2. Phép nội suy

Trong PP PTHH , các hệ số trong các hàm đa thức xấp xỉ được biểu diễn qua các

giá trị của nó ( cả những giá trị đạo hàm) tại các điểm nút được định trước trên

mỗi phần tử.

Ví dụ :

u(x) u(x) u(x)

x x x

 b a a b 2

a b a b

21 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .



u



u(

)

Nội suy xấp xỉ hằng số Nội suy tuyến tính Nội suy bậc hai

0

 a b 2



u





x

Nội suy hằng số (u)x

0

 b.u(a) a.u(b)  b a

Nội suy tuyến tính : (u)x

 u(b) u(a) b a  = a1 + a2x + a3 x2

0u

Nội suy bậc hai : (u)x

3. Dạng đa thức xấp xỉ :

Hàm xấp xỉ được chọn dưới dạng đa thức đơn giản :

Bài toán 1 –D : ( xấp xỉ tuyến tính)

( xấp xỉ bậc 2)

 n 1

 i 1

( xấp xỉ bậc 3) u(x) = a1 + a2x u(x) = a1 + a2x + a3 x2 u(x) = a1 + a2x + a3 x2 + a4x3

a x i



1

như vậy nếu u(x) xấp xỉ bậc n thì u(x) = hay là:

a a 2  a  n 1

i 1        

      

u(x) = [1 x x2 x3 ….. xn] . hay u(x) = [P(x)] . {a}

trong đó [P(x)] – ma trận các đơn thức {a} - véc tơ tọa độ tổng quát (hay véc tơ các

tham số)

Bài toán 2-D

Nếu chọn xấp xỉ bậc hai thì khi đó u(x,y) = a1 + a2x + a3y + a4x2 + a5y2 + a6xy

6

a   1   a   2         a

hay u(x,y) = [ 1 x y x2 y2 xy] . => u(x,y) = [P(x,y)] . {a}

Bài toán 3 – D có u(x,y,z) = [P(x,y,z)] . {a} (II.1)

1





  u

e

a a 2 

u(x, y, z) v(x, y, z) w(x, y, z)

y z 0 0 0 0 0 0 0 0 0 0 1 x 0 0 0 0 0 0 1

y z 0 0 0 0 . x

y z

    

    

1 x   0 0   0 0 

e

a

12

           

      

Trường hợp xấp xỉ tuyến tính có [P(x,y,z)] = [ 1 x y z ] khi đó

22 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

4. Chọn bậc của đa thức xấp xỉ.

Các đa thức xấp xỉ cần thỏa mãn được các yêu cầu sau:

 Các đa thức xấp xỉ phải thỏa mãn điều kiện hội tụ:

Do PP PTHH là một phương pháp số nên phải đảm bảo được rằng khi kích thước

các phân tử giảm đi thì kết quả tính phải hội tụ đến giá trị chính xác. Để có được

điều này thì các đa thức xấp xỉ ue phải thỏa mãn 3 điều kiện sau:

- Liên tục trong phần tử Ve . Điều này được thỏa mãn vì xấp xỉ là đa thức.

- Bảo đảm trong phần tử có trạng thái đơn vị ( hằng số) và có các đạo hàm riêng

(m)

của nó đến bậc cao nhất mà phiếm hàm I(u) đòi hỏi.

I(u)

F(x, u, u ', u",..., u

) dv

 

V

- x là tọa độ điểm, u là hàm xấp xỉ

- Trên biên phần tử , u và các đạo hàm của nó có đến cấp m -1 liên tục.

Chẳng hạn như khi u là chuyển vị thì muốn đảm bảo trạng thái đơn vị và dịch

chuyển cứng thì trong đa thức xấp xỉ không được bỏ qua số hạng a1 , hay không được

bỏ qua thành phần 1 trong [P(x,y,z)].

Khi làm mịn lưới các phần tử cần tuân theo các quy tắc sau :

+ Lưới sau mịn hơn trên cơ sở lưới trước, các điểm nút lưới trước cũng có trong

tập các nút lưới sau.

+ Các phần tử có kích thước nhỏ hơn trước, nhưng dạng hình học vẫn phải được

giữ như trước.

+ Dạng đa thức không đổi trong quá trình mịn hóa lưới phần tử.

 Các đa thức xấp xỉ được chọn sao cho không làm mất tính đẳng hướng hình

học. Dạng đa thức được chọn từ tam giác Passcal ( cho bài toán 2 chiều), tháp

Passcal cho bài toán 3 chiều.

23 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Số tham số

Bậc đa thức

Hằng số

1

Tam giác Passcal

1

tuyến tính

3

x y

bậc 2

6

bậc 3

10

x2 xy y2

x3 x2y xy2 y3

 Số các phần tử của {a} – là các tham số của đa thức xấp xỉ phải bằng số bậc

tự do của phần tử {q}e , khi đó có thể nội suy đa thức xấp xỉ thep giá trị đại

lượng cần tìm tại nút phần tử.

5. Biểu diễn đa thức xấp xỉ theo véc tơ các bậc tự do của phần tử. Ma trận hàm dạng

 Bậc tự do của một nút ( nodal degree of freedom ) là giá trị ( có thể có cả giá

trị đạo hàm ) của hàm ( đa thức ) xấp xỉ tại nút.

 Tập hợp tất cả các bậc tự do của các nút trên phần tử được gọi là véc tơ bậc tự

do của phần tử, ký hiệu là { q}e - còn được gọi là véc tơ chuyển vị nút phần

tử . Các bậc tự do này là ẩn số của bài toán PP PTHH

Ví dụ : Trong bài toán phẳng đàn hồi , khi dùng phần tử tam giác có 3 điểm nút, mỗi nút

có 2 bậc tự do : đó là các chuyển vị theo phương x và y , do đó tập hợp chuyển vị ở 3 nút

là véc tơ chuyển vị của phần tử :

T {q}e = { ui , vi , uj , vj , uk ,vk }e T

 { q1 , q2 , q3 , q4 , q5 , q6}e

Như vậy nếu phần tử e có r nút , mỗi nút có s bậc tự do thì véc tơ chuyển

vị phần tử {q}e có số thành phần ne = s × r

Các đa thức xấp xỉ được biểu diễn theo các bậc tự do của phần tử {q}e , và phải thỏa

mãn điều kiện : Các giá trị của đa thức xấp xỉ ( có thể cả đạo hàm của nó) tại các điểm

nút thuộc phần tử phải đồng nhất bằng giá trị các bậc tự do của phần tử

24 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

1

1

1

1

1

2

2

2

2

2

Tức là :



. {a}



{q} e

u(nút1) u(nút 2) 

u(x , y , z ) 1 u(x , y , z ) 2 

u(nút1) u(nút 2) 

[P(x , y , z )]

u(nút r)

u(nút r)

u(x , y , z ) r

r

r

r

r

r

      

      

      

      

      

[P(x , y , z )]     [P(x , y , z )]           

     

1

1

1

2

2

2

mà

 . {a} {q} e

 {a} {q} e

[P(x , y , z )]

r

r

r

[P(x , y , z )]   [P(x , y , z )]     

     

do vậy , ta sẽ có : hay là [A]. (II.2)

Chú ý rằng ma trận [A] là ma trận vuông có ne × ne phần tử và chỉ chứa tọa độ các điểm nút của phần tử => {a} = A-1.{q}e . Sau khi tính được {a} ta thay vào (II.1) tức là :

u(x,y,z) = [P(x,y,z)].{a}

=> (II.3) u(x,y,z) = [P(x,y,z)]. A-1.{q}e hay u(x,y,z) = [N].{q}e

với [N] = [P(x,y,z)]. A-1 (II.4)

được gọi là ma trận các hàm nội suy hay là ma trận hàm dạng

Ví dụ 1: Tìm ma trận hàm dạng của phần tử thanh lăng trụ chịu kéo, nén dọc trục ( Biết

thanh chiều dài L , diện tích mặt cắt ngang : F , mô đun đàn hồi E)

ε



x

du dx

Giải Tại mọi điểm chỉ tồn tại chuyển vị và biến dạng dọc trục : u(x) và

U



T {ε} σ dv





ε σ dx x

x

ε E.ε dx x

x







F 2

F 2

1 2

L

L

V

do đó phiếm hàm thế năng đàn hồi :

U



2 (ε ) dx

(

2 ) dx

x





EF 2

EF 2

du dx

L

L

=> =

thấy rằng đa thức xấp xỉ u(x) chỉ đòi hỏi tồn tại đạo hàm bậc nhất, do đó có thể lấy u(x) là

1

hàm xấp xỉ tuyến tính : u(x) = a1 + a2x ( 0 ≤ x ≤ L ) hay :

u(x)



[1 x ] .

a 1 a

a a

2

2

  

  

  

  

= [P(x)] . {a} => [P(x)] = [ 1 x ] , còn {a} =

Như vậy {a} có 2 tham số, nên véc tơ chuyển vị nút {q}e của phần tử cũng chỉ có 2 bậc tự

do - là 2 chuyển vị dọc trục x tại hai đầu nút : 25 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

T T ≡ {u1 u2 }e

u(nút1)

u(x



0)

{q}e = {q1 q2}e









a

a 1 a

q 1 q

u(nút 2)

 u(x L)

1

a 1  a L 2

2

2

  

  

  

  

1 0   1 L 

  

  

  

  

  

  

  

ta có :



 P(x ) 1  P(x ) 2

1 0   1 L 

  

    

  

1    1  L

0 1 L

   



[N (x) N (x)]

1



vậy [A] = => A-1 = =>

2

1

x L

x L

  

  

  

  

   

1  1 L

0 1 L

    

hàm dạng : [N]e = [P(x)]. [A]-1 = [1 x ] .

x L

x L

  1 

  

  

  

- hàm dạng của phần tử thanh chịu kéo , nén dọc trục [N]e =

Thay vào (II.3) ta được biểu diễn xấp xỉ chuyển vị u(x) theo các chuyển vị tại các nút:

x L

x L

  1 

  

u(x) = [N]e . {q}e = q1 + q2

Các hàm ui(x) còn được gọi là các hàm nội suy Lagrange bậc 1.

Ví dụ 2. Chọn đa thức xấp xỉ và tìm ma trận hàm dạng phần tử dầm 2 điểm nút chịu uốn.

Giải

Bỏ qua biến dạng dọc trục, khi đó trạng

thái chuyển vị của dầm được đặc trưng bởi

chuyển vị v(x) theo phương vuông góc với

trục dầm ( trục 0x)

Do sự tương thích biến dạng, nên tại biên giữa các phần tử cần đảm bảo sự liên tục về

θ



dv dx

chuyển vị v và góc xoay => số bậc của phần tử là 4 : đó là các chuyển vị và

góc xoay tại hai đầu:

{q}e = {v1 θ 1 v2 θ 2 }e

T T

(a) ≡ {q1 q2 q3 q4 }e

Như vậy để nội suy hàm xấp xỉ v(x) thì véc tơ tham số {q}e cũng phải có 4 tham số =>

đa thức xấp xỉ của v(x) là đa thức bậc 3 ( vì đa thức bậc 3 có 4 hệ số)

26 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

1

2

3

a a a a

4

     

     

1

2

v(x) = a1 + a2x + a3x2 + a4x3 = [ 1 x x2 x3 ] . = [Pv(x)]. {a}

dv dx

3

a a a a

4

     

     

Góc xoay θ (x) = = a2 + 2a3x + 3a4x2 = [ 0 1 2x 3x2 ] . = [Pө(x)]. {a}

Từ (a) ta có :

a

1

x 0v

 



a

hay ở dạng ma trận : q1 ≡ v1 =

2

0 1 0 1

0 0

0 0

2

dv dx  x 0



2 1 L L

3 L

3

q2 ≡ θ 1 =

x Lv 

q 1 q q q

4

0 1

2 2L 3L

      

      

a   1   a   2   a   3   a   4

      

      



a





q3 ≡ v2 = = a1 + a2L + a3L2 + a4L3

2

2a L 3

2 3a L 4

dv dx  x L

q4 ≡ θ 2 =  {q}e = [A]. {a}

2

Chú ý rằng u(nút k ) và P(x,y,z) trong (II.2) ở đây sẽ là :

v θ

v 1 θ 1

2

  

  

  

  

2

3

1 x x

x

u(nút 1) = ; u(nút 2) =



2

P (x) v P (x) θ

0 1 2x 3x

  

  

    

   

=> [P(x,y,z)] =

1 0 0 1

0 0

0 0

1 0

0 1

0 0

0 0

từ đó ta cũng có thể thấy ngay :

[A]





2 1 L L

3 L

2 3 / L



2 2 / L 3 / L



3 1/ L



0 1

2 2L 3L

3 2 / L

2 1/ L

3 2 / L



2 1 / L

[P (0)]  v  [P (0)]  θ  [P (L)] v  [P (L)]  θ

     

      

      

      

      

=> A-1 =

27 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

1 0

0 1

0 0

0 0

Do đó ta có ma trận hàm dạng của chuyển vị v(x) sẽ là : [Nv] = [Pv(x)] . [A-1], tức là :

1 Nv

2 Nv

3 Nv

4 ]

2 3 / L



2 2 / L 3 / L



3 1/ L



3 2 / L

2 1/ L

3 2 / L



2 1 / L

      

      

[Nv(x)] = [1 x x2 x3 ] . = [Nv

2

3

3

2

  1

; N



x







v N 1

v 2

x 2 L

3x 2 L 2

2x 3 L 3

2x L 2

3





; N

 

N



v 4

v 3

x L

2x 3 L

3x 2 L

x 2 L

với :

4

Khi đó hàm chuyển vị v(x) của phần tử dầm chịu uốn là :

v N .q i

i



 i 1

v(x) = [Nv(x)] . {q}e =

1 0

0 1

0 0

0 0

Còn ma trận hàm dạng của góc xoay Ө(x) sẽ là : [Nө] = [Pө(x)] . [A-1]

1 Nө

2 Nө

3 Nө

4 ]

2 3 / L



2 2 / L 3 / L



3 1/ L



3 2 / L

2 1/ L

3 2 / L



2 1 / L

      

      

[Nө] = [0 1 2x 3x2 ] = [Nө

II.4 Các phương trình cơ bản

1. Chuyển vị , biến dạng và ứng suất trong phần tử. Ma trận cứng phần tử , véc

tơ tải phần tử.

Do ở đây ta giải bài toán theo mô hình tương thích ( phương pháp chuyển vị) , nên

đại lượng cơ bản cần xác định đầu tiên là chuyển vị. Chuyển vị được xấp xỉ hóa và nội

suy theo véc tơ chuyển vị nút các phần tử {q}e.

 Sau khi tìm được ma trận các hàm dạng [N], ta sẽ biểu diễn được trường

chuyển vị theo các bậc tự do ở các nút phần tử {q}e :

(II.5) {u}e = [N]. {q}e

28 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

 Ta có các phương trình Cô-si cho sự liên hệ giữa trạng thái biến dạng và

chuyển vị, do đó trạng thái biến dạng của các điểm trong phần tử sẽ là :

 



  ε

   u

   

[N]{q} e

[B] {q} e

e

e



(II.6)

0

0



0

 y



0

0

(II.7) trong đó [B] =  [N]





0



y

 x





0

 z

y







0

  x  0             

 z

      z        

 x

, [N] – hàm dạng của chuyển vị Chú ý :  

và [B] được gọi là ma trận tính biến dạng

 Trạng thái ứng suất tại một điểm thuộc phần tử trong trường hợp vật liệu

0

0

tuân theo Hooke sẽ là :







  σ

  [D]. ε

e

e

  ε

  σ

(II.8)





e

e

e

e

0

0

là biến dạng và ứng suất ban đầu của phần tử. với  0 ε ,  0 σ

  σ

  [D].[B] q

e

e

  σ

e

e

0

0

Sử dụng (II.6) ta được :   

  [T] q

e

e

  D ε

  D ε   σ

e

e

hay là :    (II.9)   σ

trong đó : [T] = [D]. [B] là ma trận tính ứng suất (II.10)

Như vậy nhờ có (II.5) , (II.6) , (II.9) ta có thể biểu diễn chuyển vị, biến dạng, ứng suất

trong phần tử theo véc tơ chuyển vị nút phần tử {q}e

Thế năng toàn phần của phần tử Ve theo (I.19) là :



 ({u}) U A







{u} {g} dv



{u} {p} ds

  {ε} σ dv

e

e

e

T e

e

T e

e

T e

e







1 2

V e

S e

V e

(II.11)

Sử dụng (II.5) , (II.6) và (II.9) vào (II.11) ta được :

29 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

T





{q} [B] [D][B] q dv



 

({q} ) e

e

T e

e



1 2

V e

0

T

0



{u} {g} dv



{u} {p} ds



dv



{q} [B] σ

dv

e

T e

e

T e

T e

T e

  T {q} [B] [D] ε

 









e

e

1 2

1 2

V e

S e

V e

V e

   

   

T . [N]T , thay vào ta được :

T = {q}e

Chú ý {u}e = [N]. {q}e => {u}e







  {q} [K] q

({q} ) e

e

T e

e

T {q} {P} e e

e

1 2

(II.12)

T

trong đó :[K]e - Ma trận cứng phần tử

[B] [D][B]dv

[K] e

 

V

(II.13)

T

T

T

0

0

và {P}e - véc tơ tải phần tử :



{N} {g} dv



{N} {p} ds



[B] [D] ε

dv



dv

{P} e

e

e

 

  T [B] σ









e

e

1 2

1 2

S e

V e

V e

V e

   

   

(II.14)

Nhận thấy rằng do [D] là ma trận đối xứng nên tích [B]T [D] . [B] là ma trận đối xứng

do đó [K]e là ma trận đối xứng.

Ví dụ 1: Tìm ma trận cứng và véc tơ tải của phần tử thanh lăng trụ chịu kéo, nén dọc

trục ( thanh chiều dài L , diện tích mặt cắt ngang : F , mô đun đàn hồi E)

ε



x

du dx

Giải Tại mọi điểm chỉ tồn tại chuyển vị và biến dạng dọc trục : u(x) và

thấy rằng đa thức xấp xỉ u(x) chỉ đòi hỏi

tồn tại đạo hàm bậc nhất, do đó có thể lấy

u(x) là hàm xấp xỉ tuyến tính :

u(x) = a1 + a2x ( 0 ≤ x ≤ L )

hay :

u(x)



[1 x ] .

a 1 a

2

  

  

= [P(x)] . {a} , theo kết quả của ví dụ 1 phần (II.2) ta có :

biểu diễn xấp xỉ chuyển vị u(x) theo các chuyển vị tại các nút:

30 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

x L

x L

  1 

  

u(x) = [N]e . {q}e = q1 + q2

x L

x L

  1 

  

  

  

(II.15) như vậy hàm dạng của chuyển vị là [N]e =

σ



E.ε



[D].{ε}

x

x

ứng suất chỉ còn một thành phần => [D] = E

x L

x L

d dx

x L

d x dx L

d dx

  

  

  1 

  

  1 

  

  

  

  

  

= Ma trận tính biến dạng : [B] = [  ] [N] =









 1 1

1 1 L L

1 L

  

  

L

T



[B] [D][B]dv



E

[ 1 1] .F . dx =



[K] e





1  1

1 L

1 L

E.F L

0

  

1    1 

  

1    1 

V e

[B] = , do vậy ma trận cứng phần tử theo (II.13):

[K] = e

E.F L

1

1

1   

1    

ma trận cứng của phần tử thanh chịu kéo nén dọc trục: (II.16)

( E – mô đul đàn hồi , F – thiết diện thanh , L – chiều dài thanh )

Véc tơ tải phần tử {P}e được tính theo (II.14) là :

L

T

T

x L

P =

[N] {p} ds



[N] {p(x)} dx

p(x)dx

* Do lực phân bố dọc trục :

e

e





0

S e

 1  L    x  0   L

     

P =

{P}e , nếu p(x) = p0 hằng số

p L 0 2

1     1  

(II.17) thì {P}e

0

T

0

0

i =

* Do nhiệt độ thì :

dv



E α.T dx



  [B] [D] ε



e

1 2

E α.T 2

L 1 1  2 L 0

  

1    1 

  

1    1 

V e

(II.18) {P}e

T0 độ biến thiên nhiệt độ

2. Ghép nối các phần tử - Ma trận cứng – Véc tơ tải tổng thể.

a) Ghép nối các phần tử:

31 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Vật thể V được chia thành Ne miền con hay là phần tử ( Ve) bởi R điểm nút. Nếu mỗi nút có s bậc tự do thì số bậc tự do của cả hệ sẽ là : n = R × s . Gọi  q là véc

tơ chuyển vị nút tổng thể ( còn gọi là véc tơ chuyển vị nút kết cấu) ,  q có n thành phần

và là tập hợp tất cả các bậc tự do của tất cả các nút của hệ.

Nếu mỗi phần tử có r nút khi đó số bậc tự do của mỗi phần tử sẽ là : ne = r × s.

Véc tơ chuyển vị nút phần tử {q}e có ne thành phần gồm tất cả các bậc tự do của r nút

của phần tử.

Ta có các thành phần của véc tơ {q}e sẽ nằm trong các thành phần của véc tơ  q , do vậy

ta có :

(II.19) {q}e = [L]e .  q

(ne × 1) (ne × n) (n× 1)

trong đó [L]e được gọi là ma trận định vị của phần tử có kích thước (ne × n)

Ví dụ 1: Dầm với 4 nút , mỗi nút có 2 bậc tự do như hình vẽ :

 q = { q1 q2 q3 q4 q5 q6 q7 q8}T

Véc tơ chuyển vị nút tổng thể  q sẽ là :

q 1

2







{q} 1

[L] .{q} 1

q 2 

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

3

q 1 q q q

q

0 0 0 1

0 0 0 0

4

8

1   0 1   0 0 1  

     

      

      

      

      

q

3

q 1

4







{q} 2

[L] .{q} 2

q 2 

0 0 1 0 0 0 1 0 0 0 0 1

0 0 0 0 0 0 0 0 0 0 0 0

5

q q q

q

0 0 0 0 0 1

0 0

6

8

     

     

      

      

      

      

q

5

q 1

6







{q} 3

[L] .{q} 3

q 2 

0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1

0 0 0 0 0 0

7

q q q

q

0 0 0 0 0 0 0 1

8

8

     

     

      

      

      

      

khi đó ta sẽ có :

32 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

b) Ma trận cứng - véc tơ tải tổng thể.

N

N

e

e

 







  {q} [K] q

({q} ) e

e

T e

e

T {q} {P} e e





e

1 2

 e 1

 e 1

  

  

T

T

Sử dụng (II.12) ta sẽ có thế năng toàn phần :

 

{q} [L] [K] L {q}



 

e

T e

{q} [L] {P} e

T e

e

  

eN 1    2  e 1

=> (II.20)

Áp dụng nguyên lý Lagrange ( nguyên lý thế năng toàn phần dừng) về điều kiện cân bằng



0



0

của toàn hệ tại các điểm nút :



δ

 

0

2

 0

     q



0

n

    q 1      q          q

N

N

e

e

hay dạng ma trận :



{q}





  [L] [K] L

  0

T e

e

T [L] {P} e e





e

 e 1

 e 1

     q

  

  

Ta có



  P

  K q  

 

dẫn đến phương trình (II.21)

K

  [L] [K] L

T e

e

eN 

e



  là ma trận cứng tổng thể

 

 e 1

     

  

(II.22) với : K

T [L] {P} e e

eN  

 e 1

(II.23) - và  P là véc tơ tải tổng thể   P

Nhận xét





{0}

  P

  K q  

 

1. Hệ phương trình biểu diễn sự cân bằng của vật thể tại các

điểm nút.

ijK của ma trận cứng tổng thể K



  biểu thị lực sinh ra ở nút i do

Phần tử

chuyển dịch đơn vị ở nút j khi tất cả các nút bị gắn cứng.

33 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

iP của véc tơ tải tổng thể  P là ngoại lực tác động lên các phần tử

Thành phần

(tính cả biến dạng và ứng suất ban đầu) được quy đổi về tương ứng với bậc tự do thứ i. Trường hợp hệ thanh còn phải kể thêm vào  P các ngoại lực tập trung tại tác động lên

các điểm nút theo các bậc tự do tương ứng.

2. Khi thiết lập phương trình (II.21) ta chưa đưa vào các điều kiện biên động học ( liên



  là suy biến ( tức là

quan đến dịch chuyển), vật thể lúc này là tự do, do đó ma trận K

 



  bằng 0) do vậy không tồn tại

*

*

*

. Nhưng sau khi đưa vào các điều định thức của K

1K      q

kiên biên động học ta sẽ dẫn về phương trình :  . K    



3. Việc sử dụng ma trận định vị [L]e để tính ma trận cứng K

P vào vị trí của nó trong K

  P  và  P thực chất là sắp   và  P . Tuy nhiên trong 



 và  e  eK 

xếp các phần tử

khi thực hành người ta sử dụng ma trận chỉ số tiện lợi hơn trong quá trình ghép nối.

 đối xứng nên K  eK 



  đối xứng.

4. * Do

  có dạng băng

* K 

* Bề rộng của băng tùy thuộc theo cách đánh thứ tự số nút

Ví dụ Xét khung phẳng 10 tầng ( như hình vẽ), 4 nhịp. Nếu không kể đến 5 nút ngàm

cứng thì hệ có 50 nút, mỗi nút có chuyển vị ( 3 bậc tự do) chưa biết => số ẩn của hệ là :

 có số phần tử = 1502 = 22500 thành phần. 



50 × 3 = 150 => ma trận K

a) Cách đánh nút theo chiều ngắn : B = 18 b) Cách đánh nút theo chiều dài : B = 33

34 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

0

ijK

Nhưng do ý nghĩa cơ học thấy rằng : khi i và j là các bậc tự do thuộc cùng một

K



0 ; K



0 ;K



0 ; K



K

 0

1,6

1,16

1,18

1,7

1,19

nút hoặc hai nút kề nhau, còn lại là bằng 0, chẳng hạn :

=> chiều rộng của băng B = (S+1) × s , với S sự sai khác lớn nhất của mã số nút kề

nhau.

Với cách đánh số nút như hình a) có : B = ( 5 + 1). 3 = 18

b) có B = (10 +1) .3 = 33

Như vậy với ví dụ trên cách đánh số nút theo chiều ngắn như hình a) cho bề

rộng băng B là nhỏ nhất. Như thể sẽ có lợi cho việc lưu trữ số liệu và việc giải hệ phương

trình đại số .

3. Phép chuyển trục tọa độ

Trong các phần trên ta đều tính giá trị các đại lượng : chuyển vị, biến dạng ,

ứng suất, ma trận hàm dạng [N]e , ma trận cứng phần tử [K]e , véc tơ tải phần tử {P}e đều

trong hệ tọa độ thích hợp cho mỗi phần tử - là hệ tọa độ địa phương ( local coordinate

system), do vậy phương của các bậc tự do cũng lấy theo hệ tọa độ này.

Thực tế các kết cấu mà các phần tử khác nhau thì có các hệ tọa độ địa phương

khác nhau, do đó các bậc tự do của các phần tử cũng khác nhau về phương. Chính vì lẽ

đó cần có một hệ tọa độ chung cho toàn hệ - gọi là hệ tọa độ tổng thể ( global coordinate

system).

Nếu gọi :

* hệ tọa độ địa phương là xyz và hệ tọa độ tổng thể là x’y’z’

* {q}e , {P}e , [K]e là véc tơ chuyển vị, véc tơ tải , ma trận cứng phần tử

trong hệ tọa độ địa phương xyz

* {q’}e , {P’}e , [K’]e là véc tơ chuyển vị, véc tơ tải , ma trận cứng phần tử

trong hệ tọa độ tổng thể x’y’z’

khi đó ta có mối liên hệ giữa chúng :

(II.24) {q}e = [T]e . {q’}e và {P}e = [T]e {P’}e

35 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

ma trận {T}e là ma trận biến đổi ( transformation matrix) các thành phần chuyển vị nút

từ hệ tọa độ tổng thể x’y’z’ về hệ tọa độ địa phương xyz

Chú ý :

T = [I] ma trận đơn vị, ( trường hợp nếu [T]e vuông thì nó

 Luôn có [T]e .[T]e

sẽ là ma trận trực giao, khi đó [T]T = [T]-1)

 Số thành phần của {q’}e có thể khác số thành phần của {q}e



  {q} [K] q

e sẽ là :

  e

T e

e

T {q} {P} e e

e

1 2

Do đó thế năng toàn phần của phần tử



  [K] [T] q '

  e

T T {q'} [T] e e

e

e

T {q'} [T] .[T] {P '} e e

T e

e

e

1 2

=>



  [K] [T] q '

  e

T T {q'} [T] e e

e

e

T {q'} {P '} e e

e

1 2

=> ,





  {q'} [K'] q '

  e

T e

e

T {q'} {P '} e e

[K '] e

T [T] e

[K] [T] e e

e

1 2

eN



{P '}

hay là . với (II.25)

  [K'] q '

e

   , áp dụng nguyên lý Lagrange => :

e 1

Do (II.26)

[K']

[K'] e

eN  

 e 1

Trong đó - ma trận cứng tổng thể trong hệ tọa độ tổng thể.

 q '

- véc tơ chuyển vị nút tổng thể trong hệ tọa độ tổng thể.

{P '} =



[P'] e

 {P } n

eN 

 e 1

véc tơ tải tổng thể trong hệ tọa độ tổng thể.

{P } n

với - véc tơ tải tập trung tại các nút tác dụng theo các phương tương ứng của các

- gọi là véc tơ tải trọng nút ( nodal load thành phần véc tơ chuyển vị nút kết cấu  q '

vector).

Ví dụ: Xét phần tử thanh chịu kéo nén dọc trục :

36 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

nếu lấy hệ tọa độ địa phương 0x dọc trục, 0y vuông góc trục thanh, 0x lập với 0’x’ góc α

Gọi các véc tơ cơ sở chính tắc của hệ 0xy là E = {e1 , e2} còn cơ sở chính tắc của 0’x’y’

cosα



sin α

là E’= { e’1 , e’2 } khi đó ma trận chuyển cơ sở từ E’ sang E là P với

sin α

cosα

  

  

P = ( ma trận P có các cột là tọa độ của E biểu diễn qua E’)

cosα

sin α

 1

 do đó P-1 là ma trận chuyển cơ sở từ E sang E’, do đó với véc tơ x bất kỳ ta sẽ được:

P

  1 P [x]

 [x] E

E '

sin α cosα

   

  

cosα

sin α

, có , chuyển vị q1 trong hệ tọa độ 0xy có





 [q ] 1 E

q 1 0



sin α cosα

  

  

  

 q  1    q   2

  

tọa độ là

 q1 = q’1 . cos α + q’2 . sin α , tương tự chuyển vị q2 biểu diễn

qua chuyển vị của hệ tọa độ tổng thể là

 q2 = q’3 . cos α + q’4 . sin α

0



q 1 q

cosα sin α 0 0 0

cosα sin α

2

  

  

  

4

 q  1   q    2    q  3    q 

Do vậy ta có :

cosα sin α 0

0

do đó :

0

cosα sin α

0

  

  

(II.27) [T]e =

Nếu gọi  = cos α , m = sin α khi đó ta có

37 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .



m 0 0 

m

0 0

  

  

(II.27)’ [T]e =

[K] = e

E.F L

1

1

1   

1    

Theo (II.16) ta có ma trận cứng trong hệ tọa độ địa phương là :





Áp dụng (II.25) thì ma trận cứng của phần tử trong hệ tọa độ tổng thể là



[K '] e

T [T] e

[K] [T] e e

0 

0 0 

0

m 0

0 m

m 0

E.F L

1      1  

  

0

m

     

  1      1  

2

2







 m





m

2

2



m m



 m m 

=

2

2

E.F L











m

 m

2

2





 m m 

m m

      

      

(II.28) => [K’]e =

4. Xây dựng ma trận cứng tổng thể - Véc tơ tải tổng thể.

Sử dụng ở đây hai hệ thống chỉ số để đánh số cho các bậc tự do của các nút :

1. Hệ thống chỉ số tổng thể:

Đánh số các bậc tự do của toàn hệ thống tức là thứ tự của các bậc tự do đang xét

, nó được đánh thứ tự là 1 , 2 , 3 ,…, n = R× s trong  q hoặc  q

2. Hệ thống chỉ số phần tử : Để chỉ thứ tự các bậc tự do trong phần tử , hay là thứ

tự các bậc tự do trong {q}e hoặc {q’}e , nó được đánh số 1, 2 , 3 , …, ne = r × s

trong đó R – số nút của cả hệ , r – số nút của phần tử , s – bậc tự do của một nút.

Để xác định tương ứng mỗi thành phần của {q}e trong  q ( hoặc của {q’}e trong  q )

người ta đưa ra khái niệm ma trận chỉ số [b] ( còn gọi là ma trận liên hệ Boolean) mà

mỗi giá trị của thành phần bij chính là chỉ số bậc tự do trong hệ thống chỉ số tổng thể

tương ứng với bậc tự do thứ j trong hệ thống chỉ số phần tử của phần tử thứ i.

38 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Ma trận chỉ số [b] có số hàng bằng số phần tử của hệ : Ne , có số cột bằng số bậc

tự do của một phần tử.

Khi sử dụng ma trận chỉ số [b] để xây dựng ma trận cứng tổng thể [K] và véc tơ

e

tổng thể  P ( hoặc [K ] và  P ) thì cần nhớ rằng :

ijK của ma trận cứng phần tử [K]e ( của phần tử thứ e ) sẽ

mỗi thành phần

của ma trận cứng tổng thể [K] , được cộng thêm vào thành phần mnK

e

Điều này cũng thực hiện cho :

ijK của ma trận cứng phần tử [K’]e ( của phần tử thứ e ) sẽ được

mỗi thành phần

của ma trận cứng tổng thể [K ] , cộng thêm vào thành phần mnK

với m = bei ; n = bej ( chú ý bei , bej là giá trị của thành phần hàng e cột i và cột j

trong ma trận [b]



{P '}

5. Áp đặt điều kiện biên

  [K'] q '

b

K '

K '

1

 

 

  11

  12

1



Hệ phương trình tổng thể (II.26) : sắp xếp lại dạng như sau:

b

K '

K '

2

2

   p '     p '  

    

 

 

 

 

   q '     q '  

    

21

22

   

   

( II.29)

 b q '

2

trong đó : là véc tơ chứa tất cả các bậc tự do ( chuyển vị nút) đã biết

 1 q '

véc tơ chứa các bậc tự do chưa biết ( còn lại trong  q ' )

 b p '

1

véc tơ tải với các thành phần đã biết

 2 p '

véc tơ tải chưa biết ( còn lại trong  p ' )

b

b

K '



K '



(a)

  q '

  q '

  p '

1

2

1

 

 

  11

  12

b

K '





(b)

  q '

  q '

  p '

1

2

2

 

 

 K ' 

 

21

22

b

b

K '





K '

khi đó hệ phương trình (II.29) được viết thành hai hệ như sau :

  q '

  p '

  q '

1

1

2

 

 

  11

  12

từ hệ (a) ta nhận được : (II.30)

=> véc tơ chuyển vị tổng thể được xác định Giải hệ (II.30) ta tìm được  1 q '

39 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

K '



=> véc tơ tải tổng thể được xác định Thay  1 q ' vừa tìm được vào (b) ta tính được  2 p '

  q '

 b p '

1

1

 

  11



{P '}

Chú ý : Trường hợp các chuyển vị nút bằng 0 thì (II.30) có dạng : hệ

  [K'] q '

này nhận được từ hệ bằng cách bỏ đi các hàng , các cột ứng với các

bằng 0. bậc tự do trong  q '

Ví dụ: Xét hệ thanh phẳng cùng thiết diện F và mô đul đàn hồi E , với các kích thước

và chịu lực tập trung P như hình vẽ. Hãy thiết lập ma trận cứng tổng thể và véc tơ tải

tổng thể trong hệ tọa độ tổng thể.

Giải Mô hình hóa hệ thanh bởi 4 phần tử hữu hạn, mỗi phần tử có 2 bậc tự do ( là 2

chuyển vị theo các phương trục 0x , 0y .

Áp dụng công thức (II.28) cho từng phân tử:

3 2

1 2

 Ở phần tử (1) : α = 300 =>  = cos α = , m = sin α =

1 2 3 4 <= (Vị trí cột trong ma trận cứng tổng thể)

40 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

3

3



3



3

3

1



3



1



1 [K '] e

1 2 3

E.F 4.L

 3



3

3

3

4



3

 1

3

1

      

      



=> <= Vị trí Hàng trong MT cứng tổng thể

3 2

1 2

3

4

5

6

3



3



3

3



3

1

3



1



2 [K '] e

3 4 5

E.F 4.L

 3

3

3



3

6

3

 1



3

1

      

      

 Ở phần tử (2) : α = 1500 =>  = cos α = , m = sin α =

 3 1 2 2

 3 1 2 2

5

6

7

8

5

2



3

3

2

1 2

  2 2

 1 2

2



3

3

6

1 2

2

 1 2

  2 2



3 [K '] e

E.F 4.L

3

2



3

7

  2 2

 1 2

2

1 2

3

2



3

8

 1 2

  2 2

1 2

2

            

            

 Ở phần tử (3) : α = 150 =>  = cos α = , m = sin α =

3 2

1 2

7

8

3

4

1

3

 1



3

3



3





4 [K '] e

E.F 4.L



1



3

1

3

4



3

 3

3

3

      

 3 7  8 3   3   

 Ở phần tử (4) : α = 600 =>  = cos α = , m = sin α =

 Xây dựng ma trận chỉ số [b] :

41 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Chỉ số cục bộ Nút i Nút j

Phần tử 2 1 3 4

2 (1) 1 3 4

4 (2) 3 5 6

6 (3) 5 7 8

8 (4) 7 3 4

0

3

 3



3

3

0

0

0

0

1



3

 1

0

0

0

3

6

0



3



1



3



1

2

3



3

 3

2



3

2



3

3





3







2

1 2

2

1 2

K '

 

  

E.F 4.L

2



3

2



3

1







1 2 3 4 5 6 7 8

2

1 2

2

3

2





1



3

2

1 2

2



3

2

                   

                  3  

Đối xứng

1 , px

3 , py

1 , py

 Véc tơ tải tổng thể chỉ do tải trọng nút tạo nên. Thay tại nút 1 và nút 3 bằng các 3 , và do tại các nút 1 và nút 3 : chuyển vị phản lực liên kết : px

bằng 0 ta được :

42 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

1 p x

1 2

1 p y 0 0



p

3 x

3 4 5 6 7 8

0     0     q 3    q  4   0   0     q 7    q    8

P

          3 p  y  0  

            

K '



, véc tơ chuyển vị nút tổng thể :   q ' véc tơ tải tổng thể   P '

  . q '

  P'

 

 

=>

Do véc tơ chuyển vị nút tổng thể có các thành phần : 1 , 2 , 5 , 6 bằng 0, nên trong ma

trận cứng tổng thể ta bỏ đi các hàng và cột thứ 1 , 2 , 5 , 6 và cũng bỏ đi các thành

6

0

 1



3

0

2



3

 3

0,17 0,89

2



3

1

phần 1 , 2 , 5 , 6 trong véc tơ tải tổng thể. Khi đó sẽ dẫn về :

 1



3



1



3

0,12

4LP EF

E.F 4.L

2

2

0,52

     

     

1

2



3

P

0   0   0  

     

 q  3   q   4   q 7    q  8

 q  3   q   4   q 7    q  8



3



3

 3



3

2

2

         

         

= => =

ta xác định được các véc tơ chuyển vị nút của các Sử dụng ma trận chỉ số [b] và  q '

0,17 0,89

0 0

0,12 0,52

0 0









{q '} 1

; {q '} 2

; {q '} 3

; {q '} 4

0

0,12

0,17

0,17

EF 4.L

EF 4.L

EF 4.L

EF 4.L

0

0,52

0,89

0,89

     

     

     

     

     

     

     

     

phần tử :

43 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Chương III

TÍNH TOÁN HỆ THANH

BẰNG PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN

III.1 Hệ thanh dàn

Dàn là một hệ gồm các thanh chỉ chịu kéo nén đúng tâm ( tức là chỉ chịu biến

dạng dọc trục)

1. Phần tử thanh chịu biến dạng dọc trục (bài toán một chiều) :

Xét phần tử thanh có hai điểm nút chịu biến dạng dọc trục, chịu tải trọng phân bố

dọc trụ p(x). Như vậy phần tử có 2 bậc tự do là 2 chuyển vị của nút tại hai đầu .

Do có 2 bậc tự do nên chuyển vị dọc trục u(x) của phần tử chỉ có thể là xấp xỉ tuyến tính:

x L

x L

  1 

  

u(x) = a1 + a2x , hay là u(x) = [N]e . {q}e = q1 + q2 , L độ dài thanh.

x L

x L

  1 

  

  

  

(III.1) như vậy hàm dạng của chuyển vị là [N]e =

σ



E.ε



[D].{ε}

x

x

ứng suất chỉ còn một thành phần => [D] = E

x L

x L

d dx

x L

d x dx L

d dx

  

  

  1 

  

  1 

  

  

  

  

  

= Ma trận tính biến dạng : [B] = [  ] [N] =









 1 1

1 1 L L

1 L

  

  

[B] = (III.2)

L

T



[B] [D][B]dv



E

[ 1 1] .F . dx =



[K] e





1  1

1 L

1 L

E.F L

0

  

1    1 

  

1    1 

V e

ma trận cứng phần tử theo (II.13):

[K] = e

E.F L

1

1

1   

1    

ma trận cứng của phần tử thanh chịu kéo nén dọc trục: (III.3)

( E – mô đul đàn hồi , F – thiết diện thanh , L – chiều dài thanh )

Véc tơ tải phần tử {P}e được tính theo (II.14) là :

44 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

L

T

T

x L

P =

* Nếu do lực phân bố dọc trục :

[N] {p} ds



[N] {p(x)} dx

p(x)dx

e

e





0

S e

  1 L    x  0   L

     

P =

, {P}e

p L 0 2

1     1  

(III.4) nếu p(x) = p0 hằng số thì {P}e

0

T

0

0

i =

* Nếu do nhiệt độ thì :

dv



E α.T dx



  [B] [D] ε



e

1 2

E α.T 2

L 1 1  2 L 0

  

1    1 

  

1    1 

V e

(III.5) {P}e

T0 độ biến thiên nhiệt độ

Ví dụ Giải bài toán thanh dưới đây theo PPPTHH với sơ đồ 2 phần tử. Biết chiều dài

thanh là 2a. Diện tích thiết diện F và mô đul đàn hồi E không đổi. Thanh chịu tải phân

bố đều dọc trục và có cường độ p0 , chịu tại tập trung đầu dưới là Pt .

Các bước giải:

1. Thực hiện việc rời rạc hóa vật thể khảo sát bởi việc xác định các nút, các phần

tử và đánh số thứ tự các nút, các phần tử.

1 2

Trong ví dụ này ta chia thanh thành Ne = 2 phần tử với 3 nút ( hình (b)),

2 3

  

  

Thiết lập ma trận chỉ số [b] . Ma trận [b] =

Chỉ số cục bộ Nút đầu Nút cuối Phần tử

(1) 1 2

(2) 2 3

45 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

2. Thiết lập ma trận cứng phần tử

[K]e rồi thực hiện ghép nối các phần tử để



  , ( chú

xây dựng ma trận cứng tổng thể K

ý rằng với bài toán một chiều thì x’  x ).

1

2



K





1

EF a

2

1

1    1

1 1   

2

3

1

2

3

1



1

0

1

2

Trong ví dụ này ta có [K]1 = [K]2 và sẽ có:

K



K



2

 

  

EF a

1

3

EF a

1    1

 1   

1 1 1    1 0

2 3

    

  1    1 

=>

2. Thiết lập véc tơ tải phần tử {p}e và thực hiện ghép nối véc tơ tải tổng thể  P

Do chịu lực phân bố đều có cường độ không đổi p0 và hai phần tử thanh giống nhau nên





  p

theo công thức (III.4) ta có :

1

2

p .a 0 2

p .a 0 2

1 1     1 2  

1 2     1 3  

;   p

Nút 1 có phản lực liên kết R, nút 3 có tải trọng tập trung pt nên véc tơ tải trọng nút  n p

p

R     0      p   t

R







p a 0 2 p a 0

  P

p a 0 2

R 0 p

t

1     1 1      1  

    

    



p

t

p a 0 2

       

      

, khi đó véc tơ tải tổng thể sẽ là : sẽ là :  n



=> có phương trình như sau :

  P

  K q  

 

46 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

R



K

2

p a 0 2 p a 0

 

  

EF a

1  1 0

 1 2  1

q 1 q q

3

    

0     1      1  

    



p

t

p a 0 2

       

      

(III.6)

1q = 0, do vậy

4. Áp đặt điều kiên biên : Nhận thấy do kết cấu của bài toán nên

1q trong ma trận

trong hệ thống phương trình (III.6) ta “tạm bỏ đi” hàng và cột ứng với

2q và

3q :

p a 0



2  1

EF a



p

q 2 q 3

t

  

1      1  

  

   

p a 0 2

   

2

2

[K], tức là hàng 1 và cột 1 => dẫn về hệ hai phương trình xác định

q



p



p ; q



p



p

2

0

t

3

0

t

2a EF

2a EF

3a 2EF

a EF

0

2

5. Giải hệ phương trình trên dẫn đến



p



p

0

t

p



p

0

t

        

3a 2EF 2 2a EF

        

a EF 2a EF

2

2

(III.7) do đó véc tơ chuyển vị nút tổng thể sẽ là :   q

 2



p

2q -

t

0p a 2EF

p a 0 EF

a EF

Thay vào (III.6) xác định được phản lực liên kết : R = - =

2

0

p



p

0

t

q

2









và các véc tơ chuyển vị nút phần tử sẽ là :

  q

2

1

q 1 q

p



p

2

2 q 3

  

  

  

  

0

t

p



p

    

3a 2EF

a EF

    

0

t

      

3a 2EF 2 2a EF

a EF 2a EF

      

;   q

x a

x a

  1 

  

  

  

do vậy : Có hàm dạng của chuyển vị là [N]1 = [N]2 =

0

2

chuyển vị của mỗi phần tử :

p



0

x a

x a

3a 2EF

1 EF

  

 p x  t 

p



p

  1 

  

  

  

0

t

    

3a 2EF

a EF

    

= u1(x) = [N]1 {q}1 =

47 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

2

p



p

0

t

x a

x a

  1 

  

  

  

p



p

0

t

      

3a 2EF 2 2a EF

a EF 2a EF

      

2

2

= (III.8) u2(x) = [N]2 {q}2 =

p



p

p



p

0

t

0

t

3a 2EF

a EF

x a a 2EF

a EF

  

  



 a

+ => u2(x) =

Chú ý : biến x trong biểu thức của u1(x) và u2(x) là biến địa phương : 0 x









 1 1

1 a

1 a

1 a

  

  

Ma trận tính biến dạng [B] =

xε = EF.[B].{q}e

0

2

Lực dọc trên mỗi phần tử Ne = EF.



p



p



 1 1

0

t

3a 2

1 a

p



p

0

t

    

3a 2EF

a EF

    

2

p



p

0

t

=  N1 = EF.



p



p



 1 1

0

t

1 a

a 2

p



p

0

t

      

3a 2EF 2 2a EF

a EF 2a EF

      

= (III.9)  N2 = EF.

Để nghiên cứu kỹ hơn tính chất của PPPTHH, ta xét lời giải chính xác của bài toán này:

Ta có phương trình vi phân :

EF



p(x)



0

d dx

du(x) dx

  

  

, trong đó p(x) là cường độ tải trọng phân

2

bố dọc trục thanh, ở bài toán này p(x) = p0 - hằng số , do đó ta nhận được bài toán :

EF



p



0 ; 0



x



2a

0

d u(x) 2 dx

(III.10)

u(x)



0

x



0

và các điều kiện biên :

N



EF



p

t

 x 2a

du(x) dx

 x 2a

(III.11)

48 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Giải phương trình (III.10) cùng các điều kiện biên (III.11) ta nhận được nghiệm chính xác

2

của bài toán :

p





2p







 4p a x p x

t

2p a 0

p x 0

t

0

0

 

  và N(x) =

1 2EF

u(x) = (III.12)

(a) – Biểu đồ chuyển vị u(x) (b) - Biểu đồ lực dọc trục N(x)

Trên biểu đồ chuyển vị (a):

nghiệm theo PPPTHH là đường gấp khúc ABC - đường liền

nghiệm chính xác là đường cong ABC – đường nét đứt

Trên biểu đồ lực dọc (b):

nghiệm theo PPPTHH là đường gấp khúc ABCD - đường liền

nghiệm chính xác là đường EF – đường nét đứt

Như vậy để tiệm cận đến nghiệm chính xác thì ta cần tăng số nút lên , chẳng hạn số nú là

4 thì ta có biểu đồ chuyển vị và lực dọc như sau :

2. Phần tử thanh (chịu biến dạng dọc trục) trong dàn phẳng :

Trong dàn phẳng, xem mỗi mắt dàn là một đỉnh nút và thanh chịu biến dạng dọc trục, do

đó ở mỗi nút có 2 bậc tự do đó là chuyển vị theo các phương x’ và y’

49 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

x L

x L

  1 

  

(III.13) hàm chuyển vị u(x) = [N]e . {q}e = q1 + q2

x L

x L

  1 

  

  

  



(III.14) như vậy hàm dạng của chuyển vị là [N]e =



cosα ; m sin α



ta sẽ được Nếu trục phần tử thanh lập với 0x’ góc α thì gọi



ma trận chuyển tọa độ như của (II.27)’ là :

m 0 0 

m

0 0

  

  

(III.15) [T]e =

Chú ý : Nếu phần tử thanh có nút đầu thứ i , nút cuối thứ j và có tọa độ các nút trong

2

2

x

x



y

y



j

i

j

i



hệ tông thể là (xi , yi) và ( xj , yj) thì ta cùng có :

L



x



x



y



y



; m



j

i

j

i









L

L

với - độ dài thanh.

Khi đó chuyển vị nút phần tử cũng như véc tơ tải trong tọa độ địa phương và



trong tọa độ tổng thể liên hệ theo các công thức sau :



q 1 q

m 0 0  m

0 0

2

  

  

  

q ' 1 q ' 2 q ' 3 q ' 4

        

     



(III.16) {q}e = [T]e . {q’}e =>



p 1 p

m 0 0  m

0 0

2

  

  

  

p ' 1 p ' 2 p ' 3 p '

4

        

     

và (III.17) {P}e = [T]e {P’}e =>

[K] = e

E.F L

1

1

1    

ma trận cứng phần tử trong tọa độ địa phương sẽ là :



1    T [T] e

[K '] e

[K] [T] e e

ma trận cứng phần tử trong tọa độ tổng thể sẽ là

hay như (II.28) :

50 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

2

2







 m





m

2

2



m m



 m m 

Ma trận cứng của phần tử thanh chịu kéo, nén dọc trục :

2

2

E.F L











m

 m

2

2





 m m 

m m

      

      

(III.19) [K’]e =









 1 1

1 L

  

  

0

0

Ma trận tính biến dạng : [B] = (III.20)

  [D].[B] q

e

e

  D ε

1 1 L L   σ

e

e

   , trong trường hợp này ứng suất chỉ còn một Có   σ

σ



E.ε



[D].{ε}

x

x

σ



E.Fε



thành phần => [D] = E

EF.[B]{q} e

x

x

Lực dọc thanh : Ne = F.

do đó lực dọc thanh Ne trong hệ tọa độ tổng thể :





(III.21) Ne = E.F.[B]. [T]e . {q’}e = [S’]e. {q’}e







m

m



 1 1





m 0 0 

EF L

EF L

m

0 0

  

  

= (III.22) với [S’]e =

Công thức (III.21) xác định lực dọc trong phần tử do chuyển vị nút gây ra.

Trường hợp chiều dài phần tử có các lực tác dụng dọc trục hay chịu nhiệt độ thì phải kể

thêm vào (III.21) lực dọc N0 do tác dụng của lực dọc trục hay nhiệt độ gây lên trên phần

tử khi xem các nút là cố định. Ta có biểu đồ của N0 như sau:

51 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Ví dụ . Tính toán hệ thanh cùng thiết diện F và mô đun đàn hồi E cho trên hình vẽ sau :

Các bước thực hiện :

 Lập bảng thứ tự các thanh ( phần tử) với các nút đầu , cuối :

 Lập ma trận chỉ số [b]

 Tính độ dài các thanh và các góc lập bởi trục các thanh ( phần tử) ( véc tơ

nút đầu – nút cuối thanh) với 0x’.

 Tính ma trận cứng phần tử [K’]e ( Trong hệ tọa độ tổng thể)

 Ghép nối ma trận cứng tổng thể => [ K ’]

 Viết véc tơ tải tổng thể

 Áp các điều kiện biên

     p '

 K ' q ' 

 

 Loại các hàng, cột trong ’ ( liên quan đến chuyển vị đã biết :

' '  K q  * *

chuyển vị chưa biết còn lại : ở bài này là các chuyển vị = 0) => nhận được hệ phương trình xác định các     '  p *  

' kq

 Giải hệ phương trình : Tính được các chuyển vị

 Tính các chuyển vị trên từng phần tử : các u(x)

 Tính các lực dọc trục trên từng phần tử : các N(x)

52 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

1). Thứ tự các thanh và các nút đầu , cuối :

TT

Nút đầu

Nút cuối

Nút đầu - i

TT

Thanh (1) (2) (3) (4) (5) (6) (7) (8)

( i ) 1 2 3 4 5 6 2 2

( j ) 2 3 4 5 6 2 4 5

Thanh (1) (2) (3) (4) (5) (6) (7) (8)

q2i-1 1 3 5 7 9 11 3 3

q2i 2 4 6 8 10 12 4 4

Nút cuối - j q2j 4 6 8 10 12 4 8 10

q2j-1 3 5 7 9 11 3 7 9

Bảng thứ tự thanh - bậc tự do ( ma trận chỉ số [b] ) Bảng thứ tự thanh – nút

2). Tính độ dài và góc α của thanh thứ i : ký hiệu L(i) và α (i) :

2a 3

a 3

α (1) = 300 ; α (2) = 1500 ; α (3) = 300 ; α (4) = 00 ; α (5) = - 600 ; α (6) = 1800 α (7) = 900 ; α (8) = 600 ;

L(1) = L(2) = L(3) = L(7) = a ; L(5) = L(6) = L(8) = ; L(4) =

3). Tính ma trận cứng phần tử trong tọa độ tổng thể : áp dụng (III.19)

 Ta có thể viết một đoạn thủ tục ( procedure ) để tính [K’]e - là mảng 2 chiều Ke :

Khai báo cosine , sine : real { biến toàn cục) ; i , j , e : integer ;

Ke : array[1..4; 1..4] of real ;

Procedure MaTranCungPhanTu

Begin

for i := 1 to 4 do for j:=1 to 4 do Ke[i,j] := 0;

Ke[1,1]:= cosine^2 ; Ke[1,2]:= cosine * sine ; Ke[1,3]:= - cosine^2 ;

Ke[1,4]:= - cosine * sine ; Ke[2,2]:= sine ^2 ; Ke[2,3]:= - cosine* sine ;

Ke[2,4]:= - sine ^2 ; Ke[3,3]:= cosine^2 ; Ke[3,4]:= cosine * sine ;

Ke[2,1]:= Ke[1,2] ; Ke[3,1]:= Ke[1,3] ;Ke[3,2]:= Ke[2,3] ;

Ke[4,1]:= Ke[1,4] ; Ke[4,2]:= Ke[2,4]; Ke[4,3]:= Ke[3,4]; Ke[4,4]:= sine ^2 ;

for i := 1 to 4 do for j:=1 to 4 do Ke[i,j] := Ke[i,j] * E*F/L(e) ;

end ;

53 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

4). Ghép nối ma trận cứng tổng thể :

 Ta có thể viết một thủ tục ghép nối ma trận cứng tổng thể như sau:

Khai báo R = 6 { số nút của toàn hệ }

m , n , Ne : integer {Ne – số phần tử }

K : array[1 .. 2*R; 1..2*R] of real;

Procedure GhepMaTranCung

Begin

for m :=1 to 2*R do for n:=1 to 2*R do K[i,j] : = 0 ; { rửa mảng K }

for e:= 1 to Ne do

Begin

cosin := cos( α (e)) ; m := sin( α (e)); {chú ý : phải đổi góc α (e) sang radiăng }

MaTranCungPhanTu

for i:=1 to 4 do

for j:=1 to 4 do K[b[e,i], b[e,j]] := K[b[e,i], b[e,j]] + Ke[i,j] ;

end ;

end ;

5). Véc tơ tải tổng thể : Trong ví dụ này chỉ có véc tơ tải trọng nút : tại nút (6) py = p0;

các phản lực liên kết tại nút (1) và (3) , do đó có :

R

(1) 1

 Đoạn chương trình nhập các giá trị thnàh phần của véc tơ tài

tổng thể :

(1) R 2 0 0

Phần khai báo:

R

(3) 1

P : array [1..2*R] of real { biến toàn cục- R số nút của hệ }

Procedure VecToTai

  p '

Begin for i:=1 to 2*R do p[i]:= 0 ;

p[12]:= -p0 ;

{ Ở đây ta tạm gán cho các phản lực liên kết giá trí 0 , vì các giá



  p '

' *

  K ' . q '  

 

  '  K q  *

0

          (3) R   2  0  0   0  0   0   p 

                    

kiện biên để từ hệ ta dẫn về hệ  trị tại vị trí này sẽ bị loại khỏi hệ phương trình khi ta áp các điều   ' p *  

end;

54 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

' kq

6). Giải hệ phương trình : Tính được các chuyển vị

7). Tính các chuyển vị trên từng phần tử : các u(x)

8). Tính các lực dọc trục trên từng phần tử : các N(x)

3. Phần tử thanh (chịu biến dạng dọc trục) trong dàn không gian

Xem mỗi mắt dàn như là một nút thì mỗi nút có 3 bậc tự do : đó là 3 chuyển vị thành

phần theo 3 phương x’ , y’ , z’ của hệ tọa độ tổng thể. Do đó véc tơ chuyển vị nút phần

tử với nút đầu là i nút cuối là j trong hệ tọa độ tổng thể sẽ là :

 3i 2

q ' q '

Còn trong hệ tọa độ địa phương xyz : với trục thanh là trục

x thì véc tơ chuyển vị nút của phần tử thanh (với nút đầu i





  q '

e

j

 3 j 2

  q

e

u   i    u   j

u ' i v ' i w ' i u ' v ' j w '

 3i 1 q ' 3i q ' q '  3 j 1 q '

j

3j

          

          

          

          

cuối là j ) vẫn là :

trong đó ui và uj là chuyển vị nút theo phương trục thanh ( trục x).



0

0

Quan hệ giữa {q’}e và {q}e vẫn là {q}e = [T]e. {q’}e

T



e

m n 0 

0

0

0

m n

   

  

với (III.23)

trong đó  , m , n là các côsin chỉ phương

của các trục x , y , z đối với các trục x’ , y’ , z’ trong hệ tọa độ tổng thể và chúng được

x '



y '



z '



j

x ' i

j

y ' i

j

z ' i



tính theo công thức :



; m



; n



L

L

L

2

2

2

(III.24)

L



x '





y '





z '



j

x ' i

j

y ' i

j

z ' i













với

[K] = e

E.F L

1

1

ma trận cứng phần tử trong tọa độ địa phương sẽ là :

 1   



[K '] e

[K] [T] e e

1    T [T] e

ma trận cứng phần tử trong tọa độ tổng thể sẽ là

tức là :

55 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

2

2













m

n





m



n

2

2



m



m



m



m n 2

m n 2



n





nm



n

K '







e

n 2

E F L







n

m 2

m

(III.25)

m n 2

n

          

          

Đối xứng

Ví dụ :

Tính toán chuyển vị tại các nút và lực dọc của các thanh (chịu kéo nén dọc trục)

trong hệ thanh được liên kết khớp có dạng là cạnh của một chóp cụt đều : đáy và đỉnh

chóp là hình vuông có cạnh bằng 4a và 2a, chiều dài cạnh bên bằng 3a. Thiết diện các

thanh bằng nhau và bằng F, và có mô - dul đàn hồi E. Bốn góc đỉnh chịu lực tập trung

hướng xuống dưới và đều bằng P0 . Đáy tựa trên nền cứng tuyệt đối, không ma sát.

Giải:

Chú ý : Khi thực hiện việc đánh số thứ tự các thanh ( phần tử)và các nút nên

chú ý một điểm là : nếu có các thanh đồng dạng ( về kích thước và vị trị trong không

gian) thì nên chọn thứ tự nút đầu , cuối như nhau để có các cosin chỉ hướng của trục

thanh như nhau => được ma trận cứng phần tử trong tọa độ tổng thể như nhau.

1). Thứ tự các thanh và các nút đầu , cuối :

Bảng thứ tự thanh - bậc tự do ( ma trận chỉ số [b] )

Bảng thứ tự thanh – nút

56 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

TT

TT

thanh

thanh đầu i cuối j

Nút Nút đầu – i Nút cuối - j

1

(1)

2

(1)

1

2

3

4

5

6

2

(2)

3

(2)

4

5

6

7

8

9

4

(3)

3

(3)

10

11

12

7

8

9

1

(4)

4

(4)

1

2

3

10

11

12

2

(5)

8

(5)

4

5

6

22

23

24

3

(6)

7

(6)

7

8

9

19

20

21

1

(7)

5

(7)

1

2

3

13

14

15

4

(8)

6

(8)

10

11

12

16

17

18

5

(9)

8

(9)

13

14

15

22

23

24

8

(10)

7

(10)

22

23

24

19

20

21

6

(11)

7

(11)

16

17

18

19

20

21

5

(12)

6

(12)

13

14

15

16

17

18

q3i-2 q3i-1 q3i q3j-2 q3j-1 q3j

2). Tính độ dài và góc α , β , γ của thanh thứ i : ký hiệu L(i) và α (i) ; β (i) ; γ (i)

TT thanh

Độ dài cosin góc lập 0x cosin góc lập 0y cosin góc lập 0z

1

0

0

(1) (2)

4a 4a

0

1

0

(3)

4a

1

0

0

(4)

4a

0

1

0

(5)

3a

-1/3

1/3

(6)

3a

-1/3

-1/3

(7)

3a

1/3

-1/3

(8)

3a

1/3

1/3

(9)

2a

1

0

7 / 3 7 / 3 7 / 3 7 / 3 0

(10)

2a

0

1

0

(11)

2a

1

0

0

(12)

2a

0

1

0

thanh L với trục thanh:  với trục thanh : m với trục thanh : n

8 9)

(10 11 12 7 2 3 1

4 5 6

K '













  ; K '

3

1

EF 4a

1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0

 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0

1 10 2 11 3 12 7 4 8 5 9 6

        

        

3). Ma trận cứng phần tử trong tọa độ tổng thể: Áp dụng công thức (III.25)

57 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

3 10 11 12)

2 5 6 7

8 9

(1 4

K '











   ; K '

4

2

EF 4a

0 0 0 1 0 0 0 0 0 0 0 0  1 0 0 0 0 0

0 0  1 0 0 0 0 0 1 0 0 0

1 2 3 10 11 12

4 5 6 7 8 9

0 0 0 0 0 0

        

        

6

22

23

24

5

4



1



 1

1

7

7

1

1

7

1

7

 1





1

4 5 6



7

7

7

7



7

 7





5 K '

EF 27a

1

7

1

7

 1





1

7

 1



1



7

1

1

22 23 24



7



7



7

7

7

7

          

          

19

7

8

9

20

21

 1

7

1

1



 1



7

 1

7

1

1



 1

7

7

7



7



7

7



7

7 8 9





6 K '

EF 27a

1

7



1

 1

1



7

1

7



1

 1

1



7

19 20 21



7

7

 7



7



7

7

          

          

14

15

3

2

13

1

1



7



1



1

7

1

7

1



1

 1

7



1

7



7



7

7

7

 7

1 2 3





7 K '

EF 27a

7

1



1

 1

7



1

7



1



1

1



7

1

13 14 15

7

7



7



7

7



7

          

          

58 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

10

11

12

16

17

18

1

7

1



1



1



7

1

7

1



1



1



7

7

7

7



7



7

 7

10 11 12





8 K '

EF 27a



1



7

 1

1

1

7



1



7

 1

1

1

7

16 17 18



7



7



7

7

7

7

          

          

(16 17 18 19 20 21) 13 14 15 22 23 24

K '













  ; K '

9

11

EF 2a

1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0

 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0

13 16 14 17 15 18 22 19 23 20 24 21

        

        

( 13

14 15 16 17 18 )

22 23 24 19 20 21

K '













  ; K '

10

12

EF 2a

0 0 0 0 0 0

0 0 0 1 0 0 0 0 0 0 0 0  1 0 0 0 0 0

0 0  1 0 0 0 0 0 1 0 0 0

13 22 14 23 15 24 16 19 20 17 21 18

        

        

 

 

Thực hiện ghép nối ma trận cứng tổng thể K '

4) Véc tơ tải tổng thể : Do các thanh đáy trượt trên nền cứng tuyệt đối và có trọng tải tập



  :

trung tại các đỉnh đều bằng P0 nên ta có véc tơ tải tổng thể P '

59 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

0

0

R

Chú ý Vì có các chuyển vị : q ' 3 = q ' 6 = q ' 9 = q ' 12 = 0

q ' = P '

 

 . 



  ta tạm loại

1,3 0 0

do đó trong hệ phương trình K '

R

'

'

'

đi các hàng và cột : 3 , 6 , 9 , 12 như vậy sẽ dẫn

*K 

*q 

*P 

2 ,3 0

0

về hệ phương trình   .   =  

R

3 ,3 0 0

R

P '

 

4 ,3 0

0

p

0 0 0



p

0 0

0



p

0 0 0

p

0

                                         

                                      

( còn lại 20 phương trình , xác định 20 ẩn q ' k )

60 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

III.2 Hệ khung

Do liên kết giữa các thanh là các liên kết cứng nên trong các thanh có thể xuất

hiện biến dạng dọc , biến dạng xoắn tức là khi đó trong thanh có thể có lực dọc trục ,

lực trượt và các mô men uốn

1. Phần tử dầm chịu uốn – dầm liên tục

Xét phần tử dầm chịu uốn , dầm có chiều

dài L , thiết diện không đổi F. Do dầm chịu

uốn nên có thể bỏ qua biến dạng dọc trục

và khi đó chuyển vị tại x bất kỳ được đặc

trưng bởi chuyển vị v(x) vuông góc

với trục dầm và v(x) là một đa thức xấp xỉ bậc 3 của x :

v(x) = a1 + a2x + a3x2 + a4x3

4

a   1   a   2   a   3    a

hay với [P(x)] = [ 1 x x2 x3 ] và {a} =

vậy ta có v(x) = [P(x)].{a} (III.27)

Do sự tương thích biến dạng giữa phần tử đang xét và phần tử kề bên nên tại các nút cần

có sự liên tục của chuyển vị và góc xoay θ

θ



dv(x) dx

(III.28) = a2 + 2a3x + 2a4x2 = [ 0 1 2x 3x2] . {a}

Như vậy số bậc tự do tại mỗi nút là 2 dẫn đến số bậc tựa do cho mỗi phần tử dầm ( 2 nút )

1

2

là 4 , tức là ta có :





  q

e

2

3

v 1 θ v θ

q 1 q q q

2

4

      

      

      

      

, do có (III .27) và (III.28) tính tại x = 0 và x = L ta có được :

61 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

1 0

0 1

0 0

0 0

1

2







  q

e

1

L

2 L

3 L

2

3

v 1 θ v θ

q 1 q q q

2

4

0

1

2L

2 3L

      

      

      

      

a   1   a   2   a   3   a   4

      

      

1 0

0 1

0 0

0 0

hay {q}e = [A] . {a}

2 3 / L





2 / L

2 3 / L



1 / L

3 2 / L

1/ L

3 2 / L



2 1 / L

      

      

có [A]-1 =

và { a } = [A]-1 . {q}e

Thay vào (III.27) ta có được v(x) = [P(x)]. [A]-1 . {q}e

với [N] là ma trận hàm dạng của chuyển vị :

[N] = [ P(x)] .[A]-1 = [ N1 N2 N3 N4 ]

2

3

3

2

  1



; N



x





v N 1

v 2

x 2 L

3x 2 L 2

2x 3 L 3

2x L 2

3

N





; N

 



v 3

v 4

x L

2x 3 L

3x 2 L

x 2 L

với

4

Khi đó hàm chuyển vị v(x) của phần tử dầm chịu uốn là :

v N .q i

i



 i 1

(III.29) v(x) = [Nv(x)] . {q}e =

1 0

0 0

0 0

0 1

Còn ma trận hàm dạng của góc xoay Ө(x) sẽ là : [Nө] = [Pө(x)] . [A-1]

1 Nө

2 Nө

3 Nө

4 ]

2 3 / L



2 2 / L 3 / L



3 1/ L



3 2 / L

2 1/ L

3 2 / L



2 1 / L

      

      

4

[Nө] = [0 1 2x 3x2 ] . = [Nө

θ N .q i

i



 i 1

2

2

2

2

(III.30) Ө(x) = [Nө(x)] . {q}e =



1











1 =

2 =

3 =

4 =

4x L

3x 2 L

6x 2 L

6x 3 L

2x 3x  2 3 L L

6x 2 L

6x 3 L

với Nө ; Nө ; Nө ; Nө

62 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Gọi y là khoảng cách từ điểm đang xét đến đường trung hòa thì chuyển vị dọc trục u(x)

u(x)

 

y

dv(x) dx

2

và độ võng v(x) có quan hệ : , do đó biến dạng dọc trục

ε



 

y

 

y



[B]

 

y

x

{q} e

[B].{q} e

du dx

d v(x) 2 dx

2 d [N] 2 dx

2 d [N] 2 dx

với

như vậy [B] là ma trận hàng

[B]

 

y















4 L

2 L

6 2 L

12x 3 L

6x 2 L

6 2 L

12x 3 L

6x 2 L

  

  

  

  

  

  

  

  

  

  



(III.31)



  [D]. ε

xσ

với [ D] = E  Ứng suất tại mọi điểm của dầm chịu uốn : 

T

T

T

 Ma trận cứng của phần dầm chịu uốn :

K



[B] [D].[B] dv







e



E [B] .[B]dF .dx = E.F. [B] [B]dx  



L F

L

V e

=

12

6L



12

6L

do đó :

K







e

6L  12

2 4L  6L

 6L 12

2 2L  6L

EJ 3 L

6L

2 2L



6L

2 4L

      

      

(III.32)

J

2 y dF

 

F

với - mô men quán tính của mặt cắt ngang đối với đường trung hòa.

63 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

3bh 12

Nếu thiết diện dầm là chữ nhật : rộng - b ; cao - h thì J =

 Véc tơ tải {P}e : Nếu trên chiều dài phần tử thanh có tải trọng phân bố

p(x) , có lực tập trung Qi và các mô men tập trung Mi tác dụng thì trên cơ sở công thức

(II.14) ( hoặc từ sự cân bằng về công sinh ra trên các chuyển vị nút với công của ngoại

T

n

n

q

M

T

T



{N} p(x) dx





{P} e

N(x ) Q Qi

i

Mi

i





 

 



dN dx

 i 1

 i 1

  

 (x ) M  

L

lực) ta có được công thức tính véc tơ tải phần tử sau :

1) p(x) là cường độ lực ngang trên dầm

2) mô men uốn Mi tập trung tại xMi

3) Lực tập trung Qi tại xQi

2

3

0







2

3

v N 1

p 1

0

x







N

2

v 2

3x 2 L 2x L







  p

0

0

e

2x 3 L x 2 L 3

2





p p

3

0

N

L

L

v 3





 Trường hợp tải phân bố đều với cường độ p0 :

p

4

      

      

N

v 4

      

   .p dx =    

2

3

0





            

Lp 2 2 L p 12 Lp 2 2 L p 12

            

3x 2 L x L

 1            

2x 3 L x 2 L

      .p dx =       

(III.33)

3

2





3

2

a





T

2

3a 2 L 2a L



v N (a)



 .Q Q

  p

e

2a 3 L a 2 L 3

2

 

 

p 1 p p

3



p

4

      

      

2

Trường hợp lực tập trung có giá trị Q đặt cách nút đầu phần tử khoảng cách a :





3a 2 L a L

 1            

2a 3 L 3 a 2 L

            

(III. 34)

đặc biệt nếu a = L/2 ( tức là Q đặt giữa phần tử dầm) , khi đó có

64 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

T

2





N(a)

 .Q Q

  p





e

3

p 1 p p p

4

      

      

          

1   2  L   8  1   2  L  8 

(III.35)

2

 Mô men uốn nội lực trong phần tử dầm : nếu dầm có thiết diện không đổi thì:

M(x) EJ





EJ



[N]{q} e

 EJ [N ] {q} e

2

2 d v 2 dx

d dx

(III.36)

[N ]

N

N

N

v N 1

v 2

v 3

v 4



 



 



 



 

    

  

v

với

v iN (x)

iN (x) là các đa thức bậc 3, nên 

 

Để ý rằng các hàm là các hàm bậc nhất dẫn

M(nút1)

đến M(x) là hàm bậc nhất , do vậy ta chỉ cần tính được các M(nút 1) , M(nút 2) => M(x)

M(nút 2)

  

  



0)

là véc tơ mô men uốn tại các nút 1 và nút 2 của phần tử dầm, Gọi {M}e =





EJ.



  q

    q S

e

e

e

e

 N (x 

M 1 M

N (x L)



2

  

  

  

  

(III.37) =>   M





  M

    q S

e

e

e

M 1 M

2

  

  

hay là

 6

 4

6

 2

L

L

từ đây ta có được biểu đồ mô men Mq(x) – là đoạn thẳng đi qua điểm (0, M1) và (L , M2)



  S

e

EJ 3 L

6

2

 6

4

L

2 L 2 L

L

2 L 2 L

   

   

Với (III.38)

Chú ý trong trường hợp chiều dài phần tử có lực tác dụng cần cộng thêm vào Mq(x)

mô men M0(x) do tải trọng gây ra trên phần tử khi xem tất cả các nút được ngàm

65 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Các kết quả cần nhớ :

2

 Phần tử dầm độ dài L, chịu tải trọng đều với cường độ q thì ta có :

 





x

x

M (x) 0

q 2

qL 2

2 qL 12

(III.39)

Có biểu đồ mô men M0(x) như hình bên

 Phần tử dầm độ dài L , chịu tải tập trung tại điểm cách nút đầu một

khoảng là a có độ lớn p , khi đó có :



x

 a



1





p

x

M (x) 0

a 2

a L

  

     

  

với 0





(III. 40)

x

L





x

M (x) 0

 pa L a  2 L 

  

với a

Ví dụ : Tính toán chuyển vị và nội lực trong dầm liên tục được cho trên hình vẽ : Đoạn

AB chịu tải trong phân bố đều với cường độ q ; đoạn BC chịu tải tập trung cường độ

P0 = tại điểm giữa đoạn BC . Vẽ biểu đồ mô men uốn khi p0 = q.a

66 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Giải .

1) Rời rạc hóa kết cấu : Chia dầm thành 2 phần tử : Đoạn AB và BC, đánh số phần tử ,

số nút dẫn đến ma trận chỉ số [b] , chú ý ở mỗi nút có 2 bậc tự do : chuyển vị phương

Chỉ số đphương

vuông góc trục dầm, góc xoay.

Phần tử

Nút đầu (i) Nút cuối (j)

4 (1) 1 2 3

6 (2) 3 4 5

2) Thiết lập ma trận cứng phần tử

1

2

3

4

3

4

5

6

12

6a



12

6a

12

6a



12

6a

2

2

2

2

K



K











1

2

6a  12

4a  6a

 6a 12

6a  12

4a  6a

 6a 12

EJ 3 a

EJ 3 a

2

2

2

2

6a

2a



6a

4a

6a

2a



6a

4a

      

 1  2a 2   6a 3   4  

      

 3  2a 4   6a 5   6  

Sử dụng công thức (III.32)

Chú ý :

* Hệ tọa độ tổng thể và địa phương cùng hướng (trục x) nên không có sự đổi tọa

độ.

67 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

* Để tiện lợi nên viết véc tơ tải, và các điều kiện biên để từ đó loại ra các phương

*K ]

trình dẫn đến ma trận cứng phần tử [

Xét điều kiện biên: các chuyển vị nút 1 và nút 2 và 3 bằng 0 , góc xoay tại nút 1 bằng 0

=> q1 = q2 = q3 = q5 = 0 do vậy trong ma trận [ K ] ta bỏ đi các hàng, cột 1 , 2 , 3 , 5

*K ] còn lại :

4

6

2

2

8a

2a

*

K

2

2

 

  

4 6

EJ 3 a

2a

4a

   

   

dẫn đến ma trận [



3

1



4

2







p

3) véc tơ tải phần tử

  p

0

2

1

0



5

3



aq 2 2 a q 12



6

4

          

1   2  a   8  1   2  a  8 

             

aq 2 2 a q 12 aq 2 2 a q 12

            

0

aq 2 2 a q 12 p 0 2 ap 8 p 0 2 ap 8

                 

                

q

*

4

;   p => tải tổng thể   p

q

6

   

  

0



*

2 a q 12

ap 8

chú ý rằng trong véc tơ chyển vị nút chỉ còn thành phần q4 , q6 =>   q

0

ap 8

       

      

và véc tơ tải   p

0

2

2



q

8a

2a

*

*

4

2 a q 12

ap 8

4) Giải hệ phương trình :

  p

2

2

  *  K q 

q

EJ 3 a

6

2a

4a

  

  

0

   

   

ap 8

      

      

Phương trình   =  

68 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .



2

4

aq 3

p 0 2



q q

4 1

1 2

a 8EJ

6

  

     

  

      

p 0 2

      

2

2

0

0

=>

q



q









4

6

a 56EJ

2aq 3

3p 2

a 56EJ

aq 3

5p 2

  

  

  

  

giải ra ta được : và

5) Tính các đại lượng còn lại : M(x) ; v(x)

0

2

0



0 0

2 aq 3

3 p 2

a 56 EJ

  

0 0 0

Sử dụng (III.37) và (III.38) ta có









  ; q

1

2

0 q 4 0

0

2

0

2

0 q

q



4

6

      

      

      

      

0

a 56

2 aq 3

3 p 2

EJ





  

  

        

        

aq 3

5 p 2

a 56EJ

  

  

          

             



  M

    q S

1

1

e

0



0 0 0





  M

1

 3 3

3  3

1 28

 2 a a

 a 2 a

2EJ 2 a

  

2

0





0



3 ap 2 6 ap 2

2 2 a q 3 2 4 a q 3

      

      

2 aq 3

3 p 2

a 56 EJ

  

  

           

        

0

2

0





2 a q

a 56

2 aq 3

3 p 2

EJ





  





  M

2

3 3

3 3



0

1 28

ap 0 2 7

2 a a

 a 2 a

2EJ 2 a

0

  

2

      

ap 2

      

0





aq 3

5 p 2

a 56EJ

  

  

             

             

; có   q

 Nếu p0 = qa khi đó ta có

69 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

0

0





2 a q









  M

2

1

1 28

2 a q 56

1 28

2 a q 84

0

1     7  

1    5

  

0





      

ap 2 7 ap 2

      

3 ap 2 6 ap 2

2 2 a q 3 2 4 a q 3

      

      

2

2

;   M

 





x

x

M (x) 01

q 2

qa 2

qa 12

;

x



M (x) 02

a  ; 2

qa 2

a 4

  x 

  

với 0



 ( x tọa độ địa phương trên phần tử)

x

a





x

M (x) 02

a 2

qa 2

3 a 4

  

  

với

Biểu đồ mô men uốn : Giá trị của các M(x) được nhân với qa2

2

2

 





0



;

x

 a

M (x) 1

qx 2

3 qax 7

qa 14

Kết quả cuối cùng ta được các mô men uốn :

70 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

2



0





x

khi

x

a 2

M (x) 2

2









x

khi

x

a

9 qa 14 5 qa 28

qa 7 5 qa 28

a 2

     

2. Phần tử dầm chịu uốn trong hệ khung phẳng

Xét phần tử dầm chịu uốn trong hệ khung phẳng, độ dài phần tử dầm là a , có

diện tích thiết diện và mô đun đàn hồi là F và E không đổi, trong phần tử dầm chỉ có

2





chuyển vị v(x) theo phương vuông góc với trục dầm, do đó véc tơ chuyển vị nút phần tử:

e

2

3

q 1 q q q

2

4

v  1   1  v   

      

      

      

  q

Mỗi chuyển vị vi sẽ có 2 thành phần theo phương x' và y' trong tọa độ tổng thể , do vậy

véc tơ chuyển vị nút phần tử trong tọa độ tổng thể sẽ có 6 thành phần :



0

0

0

0

m

y

y

và {q}e = [T]e . {q'}e với



T





e

0 

0 0

0 0

1 0

0 0

0 m

y

y



  q

  e

0

0

0

0

0

1

      

      

q q

 u   1    v   1    1    u  2    v 2       2

          

 q  1   q  2  q  3   q  4  5     6

(III.41)

y và

ym là các cosin chỉ phương của 0y với các trục

các

0x' , 0y'.

y = - sin  ; my = cos  ,

nếu gọi  là góc lập bởi trục phần tử dầm với 0x' thì

71 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .



m

do vậy ta vẫn gọi  = cos  ; m = sin  và sẽ dẫn đến:



T





e

0 0 

0 0

0 0 0

0 0 

0 1 0

m

0

0

1

0

0

0

  0   0  

      T [K]e . [T]e ta được :

(III.42)

sử dụng (III.32) và [K']e = [T]e

2

2





 12





12

12



m

m

6 am

m

m

6 am

2

2

 12



 12

 6 a

 12 m

 6 a

2

2



4 a

6 am

 6 a

2 a

 Ma trận cứng của phần tử dầm chịu uốn trong khung phẳng:

K



   e

2

EJ 3 a

12



 12

m

m

6 am

(III.43)

2

 12

  6 a

2

4 a

 12          

          

đối xứng

 Mô men uốn tính theo trong tọa độ tổng thể :

2

2





 6 a

4 a

6 a

2 a







Theo (III.37) và (III.38) ta có

    q S

e

e

e

e

2

2

M 1 M

EJ 3 a

2



6 a

2 a

6 a

4 a

  

  

   

   

  M với   S







 S

 S



Thay {q}e = [T]e . {q'}e dẫn đến :

    q S

    q



    q





e

e

e

   S T e

e

e

e

e

e

   S T e

e

2

2







6 am

  6 a

6 am

6 a

  M với 

S



   e

 4 a 2

 2 a 2

EJ 3 a

6 am

 6 a

2 a

6 am

  6 a

4 a

   

   

(III.44)

3. Phần tử dầm chịu uốn và kéo , nén trong khung phẳng

Chuyển vị tại mỗi nút của phần tử dầm sẽ theo 2 phương : trục dầm và vuông góc với

trục dầm . Véc tơ chuyển vị nút của phần tử dầm sẽ là :

72 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .



q

e

  e

2

2

q q

2

 u   1    v   1    1    u  2    v 2       2

          

 q  1   q  2  q  3   q  4  5     6

u   1   v   1    1    u     v      

(III.45)   q và trong tọa độ tổng thể sẽ là :  

Do tính chất cộng tác dụng nên ta sử dụng ở đây các công thức (III.19) và (III.43) ta

được:

Ma trận cứng phần tử dầm chịu uốn và kéo nén dọc trục trong tọa độ tổng thể :

2

2



12



12





12





 12

 F.

2 m H

 F m

 mH

 6 amH

 F.

2 m H

 F m

mH

 6 amH

2

2

2

2





 12









12

Fm

H

 F m

 12 mH

Fm

H



 6 a H 2 4 a H

 6 a H

 6 a H 2 2 a H

K

(III.46)



   e

2

E a



12

 12

6 amH 2 m H

 F.

 F m

6 amH

 2

mH 2





12

Fm

H

  6 a H

2 4 a H

          

          

đối xứng

trong đó :  = cos  ; m = sin  ;  - góc lập bởi trục phần tử dầm và 0x'

J 2 a

a – độ dài thanh dầm F – diện tích thiết diện ; H =

Ví dụ

Tính chuyển vị và các nội lực của các

thanh trong khung phẳng chịu kéo nén dọc

trục và chịu uốn. Các thanh đều có độ dài

a, diện tích thiết diện F, mô đun đàn hồi E,

chịu tác dụng lực ngoài và các điều kiện

biên như hình vẽ. ( lực p0 lập với các thanh góc 450)

73 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Giải:

1) Rời rạc hóa các kết cấu: đánh số phần tử và số nút => tính ma trận chỉ số [b]

Phần tử Nút đầu Nút cuối

(1) 1 2

(2) 2 3

=> ma trận chỉ số [b]

Nút đầu i Nút cuối j

Chỉ số đ.ph Phần tử 1 2 3 4 5 6

(1) 1 2 3 4 5 6

(2) 4 5 6 7 8 9

2)Thiết lập ma trận cứng tổng thể:

Chú ý rằng tại các nút 1 và 3 do bị ngàm chặt nên các chuyển vị và góc xoay bằng

0 tức là ta có q'1 = q'2 = q'3 = q'7 = q'8 = q'9 = 0 nên trong các ma trận cứng phần tử

các hàng hoặc các cột là các số 1 ; 2 ; 3 ; 7 ; 8 ; 9 thì được "loại ra" do vậy áp

dụng (III.46) ta được :

 

  

0



1

; m

 2

6

4

5

4

0

6 aH

H

0

5

0



F

K

 

*  

e

E a

6

0

2 4 a H

6 aH

    

 12    

 Với phần tử (1) :

 

0

  

1



0

; m

4

5

6

4

0

0

F

5

0

12



H

6 aH

K

 

*  

e

E a

6

0

6 aH

2 4 a H

    

    

 Với phần tử (2) :

74 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

*K   ( đã tính đến các điều kiện biên) 

 Ma trận cứng tổng thể



 H F 0

12

K

0  H F

 

*  

E a

J 2 a

6 aH

6 aH

6 aH 6 aH 2 8a H

 12    

    

; H =





1 R y

1 R y





3) Véc tơ tải tổng thể :

1

1 2 3 4

1 R x 0 Q

1 R x 0 qa





5 6

2 Q 1 0

qa 0

          

          

          

          

 Có    p

 Véc tơ tải trên phần tử dầm (2) do tải

q

trọng đều q và trọng tải nút P0 tác động.

2



             

aq 2 2 a q 12 aq 2 2 a q 12

            

q

q



p

Theo (III.13) ta có   p

  p

T [T] 1

2

2



mà   

0 



1



0

; m

áp dụng (III.42) với : nên , do đó ta được

75 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

aq 2 0





nut

q





   p

2

2

R

0 0 1 0

0 0 0  1

0 0 0 0

R

2 ( ) x 2 ( ) y

2 a q 12 aq 2 0

0 0

0 0

0 0

 1  0   0  0    

0

Q  2   Q 1  0        

          

                 0   1    

aq 2 2 a q 12 aq 2 2 a q 12

            

              

2 a q 12

              

; véc tơ tải trọng nút   p '

4



Q

2

aq 2



5

aq 2 0







Q 2  Q 1 0

Q 1 2 a q 12

3 aq 2  aq 2 a q 12

nut

q













   p

   p

   p

2

2

2

R

6 7





R

R

2 ( ) x

2 ( ) x

aq 2

aq 2

R

2 ( ) x 2 ( ) y

2 a q 12 aq 2 0

R

R

2 ( ) y

2 ( ) y

8

0

          

          

2

              

2 a q 12

              

9

              

2 a q 12

              

              

q a 12

              

Do vậy véc tơ tải phần tử (2) trong tọa độ tổng thể là :

1



R

1 ( ) y

2

R

3

1 ( ) x 0

4

5

5 aq 2  2 qa

* p '

  p '

Vậy có véc tơ tải tổng thể :

=> 





6

2 a q 12



 5 aq  2     2 qa  2 a q    12

        



R

2 ( ) x

7

aq 2

8

R

2 ( ) y

9

                     

2 a q 12

                    

76 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

0

 H F

4

và ta có hệ phương trình :

0

12

K

q

 H F

5

với H =



 *  

 

*  

E a

J 2 a

6

6 aH

6 aH

    

 12    

 6 aH q    6 aH q      q 2 8a H 



          

5 aq 2 2 qa 2 a q 12

        

q , q ,q 5 6

4

Giải hệ phương trình trên ta tính được các giá trị

0

6

30

 H F

4

12



2 a q 12 E

24   1

0 6 H

 H F 6 H

6

12     

H q     6 H q   5   8H aq  

    

    

    

Chú ý : Để giải hệ phương trình trên ta nên đưa các ẩn cần tìm về cùng một thứ nguyên:

4) Tính các nội lực: mô men uốn

Chú ý:

khi tính mô men uốn trên phần tử thanh (2) cần cộng thêm mô men M0(x) theo (III.39)

III.3 Hệ khung không gian

1. Ma trận cứng phần tử : Ta xét phần tử khung không gian là dầm thẳng có thiết diện không đổi, trên mặt cắt

ngang tồn tại lực dọc, mô men uồn tròn cả hai mặt phẳng quán tính chính và mô men

xoắn. Các bậc tự do chuyển vị đặc trưng cho trạng thái chuyển vị - biến dạng của phần tử

dầm 2 nút được biểu diễn trên hình. Hệ tọa độ địa phương xyz với trục x là trục dầm, y

và z là hai trục chính của mặt phẳng cắt ngang

Véc tơ chuyển vị phần tử hai nút là :

{q}(e) = { q1 q2 q3 q4 q5 q6 q7 q8 q9 q10 q11 q12}T

77 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Trong đó :

gây ra biến dạng dọc trục dầm q1 và q7 chuyển vị dọc trục dầm

q2 và q8 chuyển vị thẳng theo trục y gây ra biến dạng uốn trong mặt phẳng xy q6 và q12 : góc xoay trong mặt phẳng xy

q3 và q9 chuyển vị thẳng theo trục z gây ra biến dạng uốn trong mặt phẳng xz q5 và q11 : góc xoay trong mặt phẳng xz

gây ra biến dạng xoắn của thanh q4 và q10 góc xoắn quanh trục x

Như vậy trong 12 bậc tự do chuyển vị sẽ gây ra 4 nhóm biến dạng độc lập nhau. Vì vậy

ma trận cứng phần tử [ K ](e) có kích thước 12 × 12 được thành lập từ 4 ma trận con như

sau :

q

q 1

a) Biến dạng dọc trục ( q1 và q7)

[ K ]



(e)

q 1 q

E.F L

1

1

7

1   

7 1    

(III.47)

78 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

b) Biến dạng xoắn theo trục x ( q4 và q10)

θ (x) là hàm bậc nhất :

Thấy rằng đa thức xấp xỉ của hàm góc xoắn

θ (x) = a1 + a2x ( 0 ≤ x ≤ L )

1

hay :

θ(x)



[1 x ] .

a a

2

  

  

= [P(x)] . {a} , theo kết quả của ví dụ 1 phần (II.2) ta có :

biểu diễn xấp xỉ hàm góc xoắn θ (x) theo các góc xoắn tại các nút:

θ (x) = [N]e . { q}xoắn =

x L

x L

  1 

  

(III.48) q4 + q10

x L

x L

  1 

  

  

  

ở đây [N]e =

r

  yz

  xy

 xyG.

 d x dx

Nếu là thanh tròn thì biến dạng xoắn sẽ là và ứng suất tiếp

xoắn = [B] {q} xoắn

 ε xoắn =   θ

áp dụng ( theo Hooke và G là mô đun trượt , r là khoảng cách từ tâm đến điểm tính toán) xoắn =  [N]{q} 

  [T] q xoắn

 σ xoắn = { yzτ }=

(III.49)

trong đó : [T] = [D]. [B] là ma trận tính ứng suất

T

gọi [Kxoắn ]e - Ma trận cứng phần tử

[B] [D][B]dv

 

V e

(III.50) [Kxoắn ]e

ở đây ta có {  } = { yz } và {  } = [ B ] . {q }(xoắn)

r L

r L

  

  

với B = ; [D ] = [ G ]

Do đó ta được :

79 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

2 r dF.



L  G dx

(III.51) [Kxoắn ]e =



1 1 L L

  

0

F

   .      

1   L  1    L

q

4

q

4

x

GJ L

1

q 10

1    1

q 10 1    

4



(III.52) => [Kxoắn ]e =

J



x

R 2

; R – bán kính dầm ở đây Jx là mô men quán tính độc cực của mặt cắt ngang:

Trong trường hợp nếu mặt cắt ngang dầm là hình chữ nhật có các cạnh a × b thì :

Jx = c ab3 với c là hằng số được lấy theo bảng sau ( phụ thuộc vào tỷ số a/b)

1 1,5 2 3 5 10 a/b

0,141 0,196 0,229 0,263 0,291 0,312 c

q

q

q

3 12

5 6L

9 12



q 11 6L

3

c) Biến dạng uốn trong mặt phẳng xz ( do { q}xz = { q3 , q5 , q9 , q11}T )

5

K







xz

e

6L  12

2 4L  6L

 6L 12

2 2L  6L

EJ y 3 L

9

q q q q 11

6L

2 2L



6L

2 4L

      

      

áp dụng (III.53)

J

2 z dF

y

 

F

với - mô men quán tính của mặt cắt ngang đối với đường trung hòa.

80 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

3ab 12

4



 Nếu thiết diện dầm là chữ nhật có kích thước các cạnh như hình vẽ thì Jy =

J



y

R 4

 Nếu thiết diện tròn bán kính R thì :

q

q

q

2 12

6 6L

8 12



q 12 6L

2

d) Biến dạng uốn trong mặt phẳng xy ( do { q}xz = { q2 , q6 , q8 , q12}T )

6

K

xy

 

   e

6L  12

2 4L  6L

 6L 12

2 2L  6L

EJ z 3 L

8

q q q q 12

6L

2 2L



6L

2 4L

      

      

áp dụng (III.54)

J

2 y dF

z

 

F

với - mô men quán tính của mặt cắt ngang đối với đường trung hòa.

3a b 12

4



 Nếu thiết diện dầm là chữ nhật có kích thước các cạnh như hình vẽ thì Jz =

J



z

R 4

 Nếu thiết diện tròn bán kính R thì :

81 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

 Từ các kết quả trên ta xây dựng ma trận cứng phần tử [ K](e) là ma trận cỡ 12 × 12

2. Chuyển ma trận cứng phần tử [K](e) về ma trận cứng [ K' ](e) trong hệ tọa độ tổng

thể

Để ghép nối các phần tử ta phải chuyển ma trận cứng phần tử về ma trận cứng [ K' ](e)

T

T

K

trong hệ tọa độ tổng thể









(e)

(e)

 T (e)

Ta có [ K' ](e) = 

[0]

[0]

[0]

trong đó [T](e) là ma trận chuyển hệ trục tọa độ , là ma trận vuông cấp 12, nó có dạng :

T





(e)

[n] [0]

[0] [n]

[0] [0]

[0]

[0]

[n]   [0]   [0]  [0] 

     [n] 

(III.55)

82 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

0 0 0

x

x

 

với :

0 0 0

y

y



0 0 0

y

m n x m n y m n z

z

    

    

    

    

[0] = và [n] = (III.56)

trong đó  x , mx , nx là cosin chỉ phương của trục x

 y , my , ny là cosin chỉ phương của trục y

 z , mz , nz là cosin chỉ phương của trục z trong hệ tổng thể x'y'z'

x '



z '

x     y     z  

  [n]. y'   

    

khi đó :

3. Các bước thiết lập chương trình tính cho hệ khung có Ne thanh với R nút:

Bước 1. Lập bảng đánh số thanh- nút : có Ne thanh với R nút

T.T Nút

thanh Đầu : i Cuối : j

1 1 2

... ... ...

k ki kj

... ... ...

Ne Nei Nej

Ở bảng này đánh số và nhập bằng tay ( nếu số lượng ít) hoặc viết một đoạn chương trình

để tự động đánh số thứ tự các thanh và các nút. Số liệu ghi vào mảng được khai báo :

Thanh_nut( 1 to Ne, 1 to 2) as integer

83 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Bước 2. Lập bảng nhập các thông số của các thanh dầm :

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

(15)

TT

L

E

G

Jx

Jy

Jz

x mx

nx

y my

ny

z mz

nz

thanh

1

...

k

...

Ne

Số liệu ở bảng này được ghi vào mảng :

Thanh(1 to Ne, 1 to 15) as double

Chú ý : Các trục chính x , y , z ở các thanh có gốc là nút đầu và phải lập thành

tam giác thuận.

Bước 3. Lập bảng ma trận chỉ số [b] : Thanh – Nút - Bậc tự do :

T.T Nút đầu : i Nút cuối : J

thanh q6i -5 q6i -4 q6i -3 q6i -2 q6i -1 q6i q6j-5 q6j-4 q6j-3 q6j-2 q6j-1 q6j

1

...

k

...

Ne

Số liệu sẽ tự tính và ghi vào mảng :

b(1 to Ne , 1 to 12) as integer

Bước 4. Tính các giá trị của [T](e) được ghi vào mảng:

T(1 to 12 , 1 to 12) as double

Bước 5. Tính [ K](e) được ghi vào mảng :

T

Ke(1 to 12, 1 to 12) as double

T

(e)

. [K](e) . [T](e) được ghi vào mảng : Bước 6. Tính [K'](e) = 

84 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

KeT(1 to 12, 1 to 12) as double

Bước 7. Ghép [K'](e) vào ma trận cứng tổng thể [K'] được ghi vào mảng:

KT ( 1 to N, 1 to N) as double

Ở đây N = 6 . R tổng số bậc tự do

Bước 8. Tính véc tơ tải phần tử => véc tơ tải tổng thể

Bước 9. Áp các điều kiện biên

Chú ý rằng sau khi có được ma trận cứng tổng thể [K'] ( mà định thức của nó = 0) , ta sẽ

tính cho véc tơ tải tổng thể trên cơ sở tính các véc tơ tải phần tử. Nhờ các điều kiện biên

ta sẽ bỏ đi được một số phương trình và một số biến ( số lượng của chúng bằng nhau),

điều đó có nghĩa là ma trận cứng tổng thể [K'] sẽ bị bỏ đi một số hàng và một số cột và

sẽ thành ma trận cứng tổng thể [K'* ] . Ma trận [K'* ] là ma trận vuông có cấp nhỏ hơn

N và định thức của nó  0 , do vậy sẽ dẫn về một hệ Cramer để giải các biến là các bậc

tự do ma ta cần tìm.

kq ( giá trị các bậc tự do theo hệ tọa

Bước 10. Giải hệ phương trình tính các chuyển vị

độ tổng), từ đó tính {q}e = [T]e. {q’}e

Bước 11. Tính các đại lượng : biến dạng , ứng suất .

3. Một số thủ tục chương trình

Chương trình được viết với các khai báo sau :

Dim Thanh_nut( 1 to Ne, 1 to 2) as integer ' đánh số các thanh và các nút

Dim Thanh(1 to Ne, 1 to 15) as double ' các thông số về các thanh

Dim b(1 to Ne , 1 to 12) as integer 'Ma trận chỉ số

Dim T(1 to 12 , 1 to 12) as double ' Ma trận chuyển tọa độ

Dim Ke(1 to 12, 1 to 12) as double ' Ma trận cứng phần tử

Dim KeT(1 to 12, 1 to 12) as double ' Ma trận cứng phần tử trong hệ tọa độ tổng thể

Dim KT ( 1 to N, 1 to N) as double ' Ma trận cứng tổng thể

Dim i , j , k as integer

85 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Với các khai báo như trên, ta có được một số thủ tục sau:

1. Thủ tục tạo lập ma trận chỉ số : [b]

Private sub TaoChiSo()

dim kk as integer

dim s , dau , cuoi as integer

For kk =1 to Ne

dau = Thanh_nut(kk,1) : cuoi = Thanh_nut(kk,2)

for s =1 to 6

b(kk,s) = 6*(dau -1) + s

b(kk,6+s) = 6*(cuoi -1) +s

next

Next

end sub

2. Thủ tục tính Te

Private sub Te(kk as integer)

dim r , s as integer

for r =1 to 12

for s = 1 to 12

T(r,s) = 0

next

next

For r =1 to 3

for s = 1 to 3

T(r,s) = Thanh(kk, 7 + (r-1)*3 + s)

T(3+r,3+s) = T(r,s)

T(6+r,6+s) = T(r,s)

T(9+r,9+s) = T(r,s)

next

next

86 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

end sub

3. Thủ tục tính ma trận cứng phần tử :

Private Sub MaTranCungPT(kk)

dim i, j as integer

dim E , G , L, F , Jx , Jy , Jz as double

for i=1 to 12

for j=1 to 12

Ke(i,j) = 0

next

Next

L = Thanh(kk,1) : F = Thanh(kk,2) : E = Thanh(kk,3) : G = Thanh(kk,4)

Jx = Thanh(kk,5) : Jy = Thanh(kk,6) : Jz = Thanh(kk,7)

Ke(1,1)= E*F/L : Ke(1,7) = -E*F/L

Ke(2,2)= 12*E*Jz/L^3 : Ke(2,6)= 6*E*Jz/L^2 : Ke(2,8)= -Ke(2,2) : Ke(2,12)= Ke(2,6)

Ke(3,3)= 12*E*Jy/L^3 : Ke(3,5)= 6*E*Jy/L^2 : Ke(3,9)= -Ke(3,3) : Ke(3,11)= Ke(3,5)

Ke(4,4)= G*Jx/L : Ke(4,10) = - Ke(4,4)

Ke(5,5) = 4*E*Jy/L : Ke(5,9) = - 6*E * Jy/L^2 : Ke(5,11) = 2*E*Jy/L

Ke(6,6)= 4*E*Jz/L : Ke(6,8) = - 6*E * Jz/L^2 : Ke(6,12) = 2*E*Jz/L

Ke(7,7) = E*F/L

Ke(8,8) = 12*E*Jz/L^3 : Ke(8,12) = - 6*E*Jz/L^2

Ke(9,9) = 12*E*Jy/L^3 : Ke(9,11) = - 6*E*Jy/L^2

Ke(10,10) = G*Jx/L : Ke(11,11) = 4*E*Jy/L : Ke(12, 12) = 4*E*Jz/L

for i =2 to 12

for j=1 to i -1

Ke(i,j) = Ke(j,i)

next

Next

End sub

87 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

4. Thủ tục tính [K'](e) : Ma trận cứng phần tử trong hệ tọa độ tổng thể ,

[K'](e) = [T]T.[K](e).[T]

Private Sub MaTranCungPT_Tong

Dim MaTranTG(1 to 12, 1 to 12) as double

dim i, j , s as integer

dim tg as double

for i =1 to 12

for j =1 to 12

MaTranTG(i,j) = 0

KeT(i,j) = 0

next

Next

For i =1 to 12

for j = 1 to 12

tg = 0

for s = 1 to 12 next tg = tg + Ke(i,s) * T(s,j) MaTranTG(i,j) = tg

next

Next

For i =1 to 12

for j = 1 to 12

tg = 0

for s = 1 to 12 next tg = tg + T(s,i) * MaTranTG(s,j)

KeT(i,j) = tg

next

Next

End sub

88 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

5. Thủ tục ghép nối ma trận cứng phần tử vào ma trận cứng tổng thể.

Private sub GhepNoi(k as integer)

dim i , j as integer

for i = 1 to 12

for j=1 to 12

KT(b(k,i) , b(k,j)) = KT(b(k,i) , b(k,j)) + Ke(i,j)

next

Next

end sub

6. Thủ tục tính toán chính , Thủ tục này được dùng sau khi đã nhập các số liệu vào các

mảng Thanh, Thanh_nut

Private Sub Tinh

Dim K , Kj as integer

call TaoChiSo

For k = 1 to N

for kj =1 to N

KT (k, kj) = 0

next

Next

For k = 1 to N

call Te(k)

call MaTranCungPT(k)

call MaTranCungPT_Tong

call GhepNoi(k)

Next

end sub

89 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .

Tài liệu tham khảo

1. Rao S.S. The Finite Element Method in Enginieering, Second Edition, Pegamon

Press, 1989.

2. Zienkiewicz O.C and Taylor R.L. The Finite Element Method, volum 1, 2.

4th Edition McGraw – Hill Book Co. , 1991

3. Dũng N.L. Giáo trình phương pháp phần tử hữu hạn trong cơ học,Trường Đại học

Bách khoa TP Hồ Chí Minh, 1993

2. Chu Quốc Thắng, Phương pháp phần tử hữu hạn, NXB Khoa học và kỹ thuật, 1997

90 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012