Bài giảng Phương pháp số (phương pháp phần tử hữu hạn)

     2 2 3 1 1 2EFK 1 1 3a       => 1 2 3 1 1 0 1 EFK 1 1 1 1 2 a 0 1 1 3                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 theo công thức (III.4) ta có :   01 1 1p .ap 1 22        ;   02 1 2p .ap 1 32        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  np sẽ là :  n t R p 0 p            , khi đó véc tơ tải tổng thể sẽ là :   0 0 0 t 0 t p a R1 R 2p aP 1 1 0 p a 2 1 p p a p 2                            => có phương trình    K q P    như sau : Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 47 0 1 2 0 3 0 t p a R1 1 0 q 2EFK 1 2 1 q p a a 0 1 1 q p a p 2                                 (III.6) 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 = 0, do vậy trong hệ thống phương trình (III.6) ta “tạm bỏ đi” hàng và cột ứng với 1q trong ma trận [K], tức là hàng 1 và cột 1 => dẫn về hệ hai phương trình xác định 2q và 3q : 0 2 0 3 t p aq2 1EF p aq1 1a p 2                  5. Giải hệ phương trình trên dẫn đến 2 2 2 0 t 3 0 t 3a a 2a 2aq p p ; q p p 2EF EF EF EF     do đó véc tơ chuyển vị nút tổng thể sẽ là :   2 0 t 2 0 t 0 3a aq p p 2EF EF 2a 2ap p EF EF                    (III.7) Thay vào (III.6) xác định được phản lực liên kết : R = - 2q - 2 0p a 2EF = 2 0 t p a a2 p EF EF   và các véc tơ chuyển vị nút phần tử sẽ là :   1 21 2 0 t 0 q q 3a aq p p 2EF EF                 ;   2 0 t 2 2 2 3 0 t 3a ap pq 2EF EFq q 2a 2ap p EF EF                  Có hàm dạng của chuyển vị là [N]1 = [N]2 = x x1 a a        do vậy : chuyển vị của mỗi phần tử : u1(x) = [N]1 {q}1 = x x1 a a        2 0 t 0 3a ap p 2EF EF          = 0 t 3a 1p p x 2EF EF      Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 48 u2(x) = [N]2 {q}2 = x x1 a a        2 0 t 2 0 t 3a ap p 2EF EF 2a 2ap p EF EF            = (III.8) => u2(x) = 2 0 t 3a ap p 2EF EF  + 2 0 t x a ap p a 2EF EF       Chú ý : biến x trong biểu thức của u1(x) và u2(x) là biến địa phương : 0 x a  Ma trận tính biến dạng [B] =  1 1 1 1 1 a a a       Lực dọc trên mỗi phần tử Ne = EF. xε = EF.[B].{q}e  N1 = EF.  1 1 1a  2 0 t 0 3a ap p 2EF EF          = 0 t 3a p p 2   N2 = EF.  1 1 1a  2 0 t 2 0 t 3a ap p 2EF EF 2a 2ap p EF EF            = 0 t a p p 2  (III.9) Để 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 : d du(x)EF p(x) 0 dx dx        , trong đó p(x) là cường độ tải trọng phân 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 : 2 02 d u(x)EF p 0 ; 0 x 2a dx     (III.10) và các điều kiện biên : x 0 tx 2a x 2a u(x) 0 du(x)N EF p dx       (III.11) Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 49 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 của bài toán : u(x) =   2t 0 0 1 2p 4p a x p x 2EF     và N(x) = t 0 0p 2p a p 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’ Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 50 hàm chuyển vị u(x) = [N]e . {q}e = x1 L      q1 + x L q2 (III.13) như vậy hàm dạng của chuyển vị là [N]e = x x1 L L        (III.14) Nếu trục phần tử thanh lập với 0x’ góc α thì gọi cosα ; m sin α  ta sẽ được ma trận chuyển tọa độ như của (II.27)’ là : [T]e = m 0 0 0 0 m         (III.15) 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 hệ tông thể là (xi , yi) và ( xj , yj) thì ta cùng có : j i j i x x y y ; m L L     với    2 2j i j iL x x y y    - độ 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}e = [T]e . {q’}e => 1 1 2 2 3 4 q ' q q 'm 0 0 q q '0 0 m q '                       (III.16) và {P}e = [T]e {P’}e => 1 1 2 2 3 4 p ' p p 'm 0 0 p p '0 0 m p'                       (III.17) ma trận cứng phần tử trong tọa độ địa phương sẽ là : e[K] = 1 1E.F 1 1L      ma trận cứng phần tử trong tọa độ tổng thể sẽ là Te e e e[K '] [T] [K] [T] hay như (II.28) : Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 51 Ma trận cứng của phần tử thanh chịu kéo, nén dọc trục : [K’]e = 2 2 2 2 2 2 2 2 m m m m m mE.F L m m m m m m                              (III.19) Ma trận tính biến dạng : [B] =  1 1 1 1 1 L L L       (III.20) Có        0 0e e e eσ [D].[B] q D ε σ   , trong trường hợp này ứng suất chỉ còn một thành phần x xσ E.ε [D].{ε}  => [D] = E Lực dọc thanh : Ne = F. x x eσ E.Fε EF.[B]{q}  do đó lực dọc thanh Ne trong hệ tọa độ tổng thể : Ne = E.F.[B]. [T]e . {q’}e = [S’]e. {q’}e (III.21) với [S’]e =  EF 1 1L  m 0 0 0 0 m         =  EF m m L    (III.22) 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: Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 52 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  Loại các hàng, cột trong    K' q ' p '    ’ ( liên quan đến chuyển vị đã biết : ở bài này là các chuyển vị = 0) => nhận được hệ phương trình xác định các chuyển vị chưa biết còn lại :    ' ' '* * *K q p     Giải hệ phương trình : Tính được các chuyển vị 'kq  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) Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 53 1). Thứ tự các thanh và các nút đầu , cuối : Bảng thứ tự thanh – nút Bảng thứ tự thanh - bậc tự do ( ma trận chỉ số [b] ) Nút đầu - i Nút cuối - j TT Thanh Nút đầu ( i ) Nút cuối ( j ) TT Thanh q2i-1 q2i q2j-1 q2j (1) 1 2 (1) 1 2 3 4 (2) 2 3 (2) 3 4 5 6 (3) 3 4 (3) 5 6 7 8 (4) 4 5 (4) 7 8 9 10 (5) 5 6 (5) 9 10 11 12 (6) 6 2 (6) 11 12 3 4 (7) 2 4 (7) 3 4 7 8 (8) 2 5 (8) 3 4 9 10 2). Tính độ dài và góc α của thanh thứ i : ký hiệu L(i) và α (i) : L(1) = L(2) = L(3) = L(7) = a ; L(5) = L(6) = L(8) = 2a 3 ; L(4) = a 3 α (1) = 300 ; α (2) = 1500 ; α (3) = 300 ; α (4) = 00 ; α (5) = - 600 ; α (6) = 1800 α (7) = 900 ; α (8) = 600 ; 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 ; Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 54 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ó :   (1) 1 (1) 2 (3) 1 (3) 2 0 R R 0 0 R Rp ' 0 0 0 0 0 p                                          Đ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ể : Phần khai báo: P : array [1..2*R] of real { biến toàn cục- R số nút của hệ } Procedure VecToTai 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á 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 kiện biên để từ hệ    K' . q ' p '    ta dẫn về hệ    ' ' '* * *K q p    end; Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 55 6). Giải hệ phương trình : Tính được các chuyển vị 'kq 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à :   i 3i 2 i 3i 1 i 3i e j 3 j 2 j 3 j 1 j 3j u ' q ' v ' q ' w ' q ' q ' u ' q ' v ' q ' w ' 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 cuối là j ) vẫn là :   ie j u q u        trong đó ui và uj là chuyển vị nút theo phương trục thanh ( trục x). Quan hệ giữa {q’}e và {q}e vẫn là {q}e = [T]e. {q’}e với  e m n 0 0 0 T 0 0 0 m n          (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 tính theo công thức : j i j i j ix ' x ' y ' y ' z ' z '; m ; n L L L       (III.24) với      2 2 2j i j i j iL x ' x ' y ' y ' z ' z '      ma trận cứng phần tử trong tọa độ địa phương sẽ là : e[K] = 1 1E.F 1 1L      ma trận cứng phần tử trong tọa độ tổng thể sẽ là Te e e e[K '] [T] [K] [T] tức là : Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 56   2 2 2 2 2 2 e 2 2 2 m n m n m mn m m mn n n nm nEFK ' L m n m mn n                                       (III.25) 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 – nút Bảng thứ tự thanh - bậc tự do ( ma trận chỉ số [b] ) Đối xứng Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 57 Nút Nút đầu – i Nút cuối - j TT thanh đầu i cuối j TT thanh q3i-2 q3i-1 q3i q3j-2 q3j-1 q3j (1) 1 2 (1) 1 2 3 4 5 6 (2) 2 3 (2) 4 5 6 7 8 9 (3) 4 3 (3) 10 11 12 7 8 9 (4) 1 4 (4) 1 2 3 10 11 12 (5) 2 8 (5) 4 5 6 22 23 24 (6) 3 7 (6) 7 8 9 19 20 21 (7) 1 5 (7) 1 2 3 13 14 15 (8) 4 6 (8) 10 11 12 16 17 18 (9) 5 8 (9) 13 14 15 22 23 24 (10) 8 7 (10) 22 23 24 19 20 21 (11) 6 7 (11) 16 17 18 19 20 21 (12) 5 6 (12) 13 14 15 16 17 18 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 thanh L cosin góc lập 0x với trục thanh:  cosin góc lập 0y với trục thanh : m cosin góc lập 0z với trục thanh : n (1) 4a 1 0 0 (2) 4a 0 1 0 (3) 4a 1 0 0 (4) 4a 0 1 0 (5) 3a -1/3 1/3 7 / 3 (6) 3a -1/3 -1/3 7 / 3 (7) 3a 1/3 -1/3 7 / 3 (8) 3a 1/3 1/3 7 / 3 (9) 2a 1 0 0 (10) 2a 0 1 0 (11) 2a 1 0 0 (12) 2a 0 1 0 3). Ma trận cứng phần tử trong tọa độ tổng thể: Áp dụng công thức (III.25)     3 1 (10 11 12 7 8 9) 1 2 3 4 5 6 1 0 0 1 0 0 1 10 0 0 0 0 0 0 2 11 0 0 0 0 0 0 3 12EFK ' ; K ' 1 0 0 1 0 0 4 74a 0 0 0 0 0 0 5 8 0 0 0 0 0 0 6 9                    Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 58     4 2 (1 2 3 10 11 12) 4 5 6 7 8 9 0 0 0 0 0 0 4 1 0 1 0 0 1 0 5 2 0 0 0 0 0 0 6 3EFK ' ; K ' 0 0 0 0 0 0 7 104a 0 1 0 0 1 0 8 11 0 0 0 0 0 0 9 12                     5 4 5 6 22 23 24 1 1 7 1 1 7 4 1 1 7 1 1 7 5 7 7 7 7 7 7 6EFK ' 2227a 1 1 7 1 1 7 231 1 7 1 1 7 247 7 7 7 7 7                                  6 7 8 9 19 20 21 1 1 7 1 1 7 7 1 1 7 1 1 7 8 7 7 7 7 7 7 9EFK ' 1927a 1 1 7 1 1 7 201 1 7 1 1 7 217 7 7 7 7 7                                    7 1 2 3 13 14 15 1 1 7 1 1 7 1 1 1 7 1 1 7 2 7 7 7 7 7 7 3EFK ' 1327a 1 1 7 1 1 7 141 1 7 1 1 7 157 7 7 7 7 7                                 Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 59  8 10 11 12 16 17 18 1 1 7 1 1 7 10 1 1 7 1 1 7 11 7 7 7 7 7 7 12EFK ' 1627a 1 1 7 1 1 7 171 1 7 1 1 7 187 7 7 7 7 7                                     9 11 (16 17 18 19 20 21) 13 14 15 22 23 24 1 0 0 1 0 0 13 16 0 0 0 0 0 0 14 17 0 0 0 0 0 0 15 18EFK ' ; K ' 1 0 0 1 0 0 22 192a 0 0 0 0 0 0 23 20 0 0 0 0 0 0 24 21                        10 12 ( 13 14 15 16 17 18 ) 22 23 24 19 20 21 0 0 0 0 0 0 22 13 0 1 0 0 1 0 23 14 0 0 0 0 0 0 24 15EFK ' ; K ' 0 0 0 0 0 0 19 162a 0 1 0 0 1 0 20 17 0 0 0 0 0 0 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 '   : Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 60 1,3 2 ,3 3 ,3 4 ,3 0 0 0 0 0 0 R 0 0 R 0 0 R 0 0 R P ' 0 0 p 0 0 p 0 0 p 0 0 p                                                                               Chú ý Vì có các chuyển vị : q '3 = q '6 = q '9 = q '12 = 0 do đó trong hệ phương trình K'   . q ' = P '   ta tạm loại đi các hàng và cột : 3 , 6 , 9 , 12 như vậy sẽ dẫn về hệ phương trình '*K   . ' *q   = ' *P   ( còn lại 20 phương trình , xác định 20 ẩn q 'k ) Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 61 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 hay với [P(x)] = [ 1 x x2 x3 ] và {a} = 1 2 3 4 a a a 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  = a2 + 2a3x + 2a4x2 = [ 0 1 2x 3x2] . {a} (III.28) 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 ) là 4 , tức là ta có :   1 1 1 2 e 2 3 2 4 v q θ q q v q θ q                            , do có (III .27) và (III.28) tính tại x = 0 và x = L ta có được : Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 62   1 1 1 1 2 2 2 3e 2 3 3 2 2 4 4 1 0 0 0v q a 0 1 0 0θ q a q v q a1 L L L θ q a0 1 2L 3L                                                  hay {q}e = [A] . {a} có [A]-1 = 2 2 3 3 2 1 0 0 0 0 1 0 0 3 / L 2 / L 3 / L 1/ L 2 / L 1/ L 2 / L 1/ L                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 ] với 2 3 2 3 v v 1 22 3 2 2 3 2 3 v v 3 42 3 2 3x 2x 2x xN 1 ; N x LL L L 3x 2x x xN ; N LL L L            Khi đó hàm chuyển vị v(x) của phần tử dầm chịu uốn là : v(x) = [Nv(x)] . {q}e = 4 v i i i 1 N .q   (III.29) Còn ma trận hàm dạng của góc xoay Ө(x) sẽ là : [Nө] = [Pө(x)] . [A-1] [Nө] = [0 1 2x 3x2 ] . 2 2 3 3 2 3 2 1 0 0 0 0 1 0 0 3 / L 2 / L 3 / L 1/ L 2 / L 1/ L 2 / L 1/ L                = [Nө1 Nө2 Nө3 Nө4 ] Ө(x) = [Nө(x)] . {q}e = 4 θ i i i 1 N .q   (III.30) với Nө1 = 2 2 3 6x 6x L L   ; Nө2 = 2 2 4x 3x1 L L   ; Nө3 = 2 2 3 6x 6x L L  ; Nө4 = 2 3 2 2x 3x L L   Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 63 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) và độ võng v(x) có quan hệ : dv(x)u(x) y dx   , do đó biến dạng dọc trục 2 2 x e e2 2 du d v(x) d [N] ε y y {q} [B].{q} dx dx dx       với 2 2 d [N][B] y dx   như vậy [B] là ma trận hàng 2 3 2 2 3 2 6 12x 4 6x 6 12x 2 6x[B] y L LL L L L L L                                  (III.31)  Ứng suất tại mọi điểm của dầm chịu uốn :    xσ [D]. ε với [ D] = E  Ma trận cứng của phần dầm chịu uốn :   e T T T e V L F L K [B] [D].[B] dv E [B] .[B]dF .dx = E.F. [B] [B]dx     = do đó :   2 2 3e 2 2 12 6L 12 6L 6L 4L 6L 2LEJK 12 6L 12 6LL 6L 2L 6L 4L                (III.32) với 2 F J y dF  - mô men quán tính của mặt cắt ngang đối với đường trung hòa. Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 64 Nếu thiết diện dầm là chữ nhật : rộng - b ; cao - h thì J = 3bh 12  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 lực) ta có được công thức tính véc tơ tải phần tử sau : q M Tn nTT e Qi i Mi i i 1 i 1L dN{P} {N} p(x)dx N(x ) Q (x ) M dx             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  Trường hợp tải phân bố đều với cường độ p0 :   2 3 0 2 3 v 2 3 211 0 v 2 2 2 0 0e v 2 3 3 0L L3 2 3v4 4 22 3 0 2 3x 2x Lp1 L L 2 Np 2x x L pxp N L L 12p .p dx = .p dx = p LpN 3x 2x 2p L LN L px x 12L L                                                                                (III.33) 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 :   2 3 2 3 2 3 1 2Tv2 e 2 3 3 2 3 4 2 3 2 3a 2a1 L L p 2a aap L Lp N (a) .Q Q p 3a 2a p L L a a L 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ó Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 65     1 T2 e 3 4 1 2 p L p 8p N(a) .Q Q p 1 2p L 8                                  (III.35)  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ì: 2 2 e e2 2 d v dM(x) EJ EJ [N]{q} EJ [N ] {q} dx dx    (III.36) với        v v v v1 2 3 4[N ] N N N N          Để ý rằng các hàm viN (x) là các đa thức bậc 3, nên  viN (x)  là các hàm bậc nhất dẫn đế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) Gọi {M}e = M(nút1) M(nút 2)       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, =>        1e e ee 2 M N (x 0) M EJ. q S q M N (x L)               (III.37) hay là      1e ee 2 M M S q M        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) Với   2 2 3 2 2 6 4 6 2 6 2 6 4e L L L LEJ S L L L L L           (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 Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 66 Các kết quả cần nhớ :  Phần tử dầm độ dài L, chịu tải trọng đều với cường độ q thì ta có : 2 2 0 2 2 12 q qL qL M (x) x x    (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ó : 0 1 2 a a M (x) p x L              với 0 x a  (III. 40) 0 2 pa L a M (x) x L       với a x L  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 Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 67 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 vuông góc trục dầm, góc xoay. Nút đầu (i) Nút cuối (j) Chỉ số đphương Phần tử (1) 1 2 3 4 (2) 3 4 5 6 2) Thiết lập ma trận cứng phần tử Sử dụng công thức (III.32)   2 2 31 2 2 1 2 3 4 12 6a 12 6a 1 6a 4a 6a 2a 2EJK 12 6a 12 6a 3a 46a 2a 6a 4a                  2 2 32 2 2 3 4 5 6 12 6a 12 6a 3 6a 4a 6a 2a 4EJK 12 6a 12 6a 5a 66a 2a 6a 4a                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 độ. Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 68 * Để 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 trình dẫn đến ma trận cứng phần tử [ *K ] 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 dẫn đến ma trận [ *K ] còn lại : 2 2 * 3 2 2 4 6 8a 2a 4EJK 6a 2a 4a            3) véc tơ tải phần tử   2 1 2 aq 1 2 a q 2 12p aq 32 a q 412                        ;   02 1 3 2 a 4 8p p 1 5 2 a 68                       => tải tổng thể   2 0 2 0 0 0 aq 2 a q 12 paq 2 2p apa q 12 8 p 2 ap 8                                  chú ý rằng trong véc tơ chyển vị nút chỉ còn thành phần q4 , q6 =>   4 6 * qq q        và véc tơ tải   2 0 * 0 apa q 12 8p ap 8              4) Giải hệ phương trình : Phương trình    * * *K q p     2 2 3 2 2 8a 2aEJ a 2a 4a        4 6 q q       = 2 0 0 apa q 12 8 ap 8             Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 69 => 0 2 4 6 0 paq q4 1 a 3 2 q1 2 p8EJ 2                      giải ra ta được : 2 0 4 3pa 2aqq 56EJ 3 2       và 2 0 6 5pa aqq 56EJ 3 2        5) Tính các đại lượng còn lại : M(x) ; v(x) Sử dụng (III.37) và (III.38) ta có có     2 0 4 1 2 2 0 2 64 0 0 0 00 320 0 56 3 2 0 00 0 32 556 3 2 3 2 pa aq q EJ q ; q pa aq qq pa aq EJ 56EJ                                                                           ;         1 1 2 0 21 2 2 0 0 0 3 20 3 2 3 2 30 283 3 2 6 4 32 2 3 56 3 2 e M S q ap a q a a2EJ 1 M a aa ap a q pa aq EJ                                         2 0 20 22 0 2 0 0 32 56 3 23 2 3 2 283 3 2 0 7 25 3 2 pa aq ap a qEJa a2EJ 1 M a a apa pa aq 56EJ                                                Nếu p0 = qa khi đó ta có Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 70   2 0 2 1 2 0 3 2 12 3 28 84 56 4 2 3 ap a q 1 a q M ap a q                  ;   20 2 2 0 12 28 567 7 2 ap a q 1 a q M ap                  2 2 01 2 2 12 q qa qa M (x) x x    ; 02 2 4 qa a M (x) x       với 0 2 a x  ; 02 3 2 4 qa a M (x) x       với 2 a x a  ( x tọa độ địa phương trên phần tử) Biểu đồ mô men uốn : Giá trị của các M(x) được nhân với qa2 Kết quả cuối cùng ta được các mô men uốn : 2 2 1 3 0 2 7 14 qx qax qa M (x) ; x a      Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 71 2 2 2 9 0 14 7 2 5 5 28 28 2 qa qa a x khi x M (x) qa qa a x khi x a            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ó 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ử:   1 1 1 2 2 3 2 4 e v q q q v q 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 :   11 21 31 42 52 62 e qu qv q q qu qv q                                         và {q}e = [T]e . {q'}e với   0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 y y e y y m T m                (III.41) các y và ym là các cosin chỉ phương của 0y với các trục 0x' , 0y'. nếu gọi  là góc lập bởi trục phần tử dầm với 0x' thì y = - sin ; my = cos , Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 72 do vậy ta vẫn gọi  = cos ; m = sin và sẽ dẫn đến:   0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 e m T m               (III.42) sử dụng (III.32) và [K']e = [T]eT [K]e . [T]e ta được :  Ma trận cứng của phần tử dầm chịu uốn trong khung phẳng:   2 2 2 2 2 2 3 2 2 2 12 12 6 12 12 6 12 6 12 12 6 4 6 6 2 12 12 6 12 6 4 e m m am m m am a m a a am a aEJ K a m m am a a                                    (III.43)  Mô men uốn tính theo trong tọa độ tổng thể : Theo (III.37) và (III.38) ta có      1e ee 2 M M S q M        với   2 2 3 2 2 6 4 6 2 6 2 6 4e a a a aEJ S a a a a a           Thay {q}e = [T]e . {q'}e dẫn đến :                e e e ee e e eM S q S T q S q     với      e e eS S T    2 2 3 2 2 6 6 4 6 6 2 6 6 2 6 6 4e am a a am a aEJ S a am a a am a 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à : đối xứng Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 73   1 1 1 2 2 2 e u v q u v                      và trong tọa độ tổng thể sẽ là :   11 21 31 42 52 62 e qu qv q q qu qv q                                         (III.45) 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ể : (III.46)   2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 12 12 6 12 12 6 12 6 12 12 6 4 6 6 2 12 12 6 12 6 4 e F. m H F m mH amH F. m H F m mH amH Fm H a H F m mH Fm H a H a H amH a H a HE K a F. m H F m mH amH Fm H a H a H                                                     trong đó :  = cos  ; m = sin ;  - góc lập bởi trục phần tử dầm và 0x' a – độ dài thanh dầm F – diện tích thiết diện ; H = 2 J a 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) đối xứng Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 74 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 :  Với phần tử (1) : 0 1 2 ; m       2 4 5 6 12 0 6 4 0 0 5 66 0 4 * e H aH E K F a aH a H               Với phần tử (2) : 0 1 0; m     2 4 5 6 0 0 4 0 12 6 5 60 6 4 * e F E K H aH a aH a H              Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 75  Ma trận cứng tổng thể *K   ( đã tính đến các điều kiện biên) 2 12 0 6 0 12 6 6 6 * H F aH E K H F aH a aH aH 8a H              ; H = 2 J a 3) Véc tơ tải tổng thể :  Có   1 1 1 1 1 2 1 1 2 30 0 4 5 60 0 y y x x R R R R p Q qa Q qa                                              Véc tơ tải trên phần tử dầm (2) do tải trọng đều q và trọng tải nút P0 tác động. Theo (III.13) ta có   2 2 2 2 12 2 12 q aq a q p aq a q                        mà    12 2 q qTp [T] p  áp dụng (III.42) với : 0  nên 1 0  ; m , do đó ta được Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 76   2 2 2 2 2 2 1 0 0 0 02 0 0 0 0 0 1 0 0 12 12 0 0 1 0 2 20 0 0 0 00 0 0 1 12 12 q aq aq a q a q p aq aq a q a q                                                                 ; véc tơ tải trọng nút   2 1 22 2 0 0 nut ( ) x ( ) y Q Q p' R R                    Do vậy véc tơ tải phần tử (2) trong tọa độ tổng thể là :       2 2 1 2 22 1 22 2 2 2 2 2 2 2 2 2 2 3 2 22 0 0 12 1212 2 22 0 0 12 12                                                                                    nut q ( ) ( ) ( )x x x ( ) y ( ) ( ) y y aq aqaq Q Q Q aq Q a q a qa q p p p aq aqRaq R R R R R a q a q a 4 5 6 7 8 912                               q Vậy có véc tơ tải tổng thể :   1 1 2 2 2 2 1 2 30 5 42 52 612 72 8 912                                            ( ) y ( ) x ( ) x ( ) y R R aq qa p ' a q aq R R a q =>   2 5 2 2 12 * aq p ' qa a q                  Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 77 và ta có hệ phương trình :   4 5 2 2 6 5 12 0 6 2 0 12 6 2 6 6 12 * * aq H F aH q E K q H F aH q qa a qaH aH 8a H a q                                       với H = 2 J a Giải hệ phương trình trên ta tính được các giá trị 4 5 6q , q ,q 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 2 5 6 12 0 6 30 0 12 6 24 12 6 6 1 H F H q a q H F H q E H H 8H aq                             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 Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 78 Trong đó : q1 và q7 chuyển vị dọc trục dầm gây ra biến dạng dọc trục dầm q2 và q8 chuyển vị thẳng theo trục y q6 và q12 : góc xoay trong mặt phẳng xy gây ra biến dạng uốn trong mặt phẳng xy q3 và q9 chuyển vị thẳng theo trục z q5 và q11 : góc xoay trong mặt phẳng xz gây ra biến dạng uốn trong mặt phẳng xz q4 và q10 góc xoắn quanh trục x gây ra biến dạng xoắn của thanh 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 : a) Biến dạng dọc trục ( q1 và q7) 1 7 1 (e) 7 q q q1 1E.F[K ] q1 1L       (III.47) Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 79 b) Biến dạng xoắn theo trục x ( q4 và q10) Thấy rằng đa thức xấp xỉ của hàm góc xoắn θ (x) là hàm bậc nhất : θ (x) = a1 + a2x ( 0 ≤ x ≤ L ) hay : 1 2 a θ(x) [1 x ] . a        = [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 = x1 L      q4 + x L q10 (III.48) ở đây [N]e = x x1 L L        Nếu là thanh tròn thì biến dạng xoắn sẽ là xyz dr dx    và ứng suất tiếp xy xyG.   ( 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) áp dụng  ε xoắn =   θ xoắn =  [N]{q} xoắn = [B] {q} xoắn  σ xoắn = { yzτ }=  [T] q xoắn (III.49) trong đó : [T] = [D]. [B] là ma trận tính ứng suất gọi [Kxoắn ]e - Ma trận cứng phần tử [Kxoắn ]e e T V [B] [D][B]dv  (III.50) ở đây ta có { } = { yz } và { } = [ B ] . {q }(xoắn) với B = r r L L      ; [D ] = [ G ] Do đó ta được : Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 80 [Kxoắn ]e = L 2 0 F 1 1 1LG dx r dF. . 1 L L L                (III.51) => [Kxoắn ]e = 4 10 4x 10 q q q1 1GJ q1 1L      (III.52) ở đây Jx là mô men quán tính độc cực của mặt cắt ngang: 4 x RJ 2   ; R – bán kính dầm 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) a/b 1 1,5 2 3 5 10 c 0,141 0,196 0,229 0,263 0,291 0,312 c) Biến dạng uốn trong mặt phẳng xz ( do { q}xz = { q3 , q5 , q9 , q11}T ) áp dụng   3 5 9 11 3 2 2 y 5 xz 3e 9 2 2 11 q q q q 12 6L 12 6L q EJ q6L 4L 6L 2L K q12 6L 12 6LL q6L 2L 6L 4L                (III.53) với 2y F J z dF  - mô men quán tính của mặt cắt ngang đối với đường trung hòa. Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 81  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 = 3ab 12  Nếu thiết diện tròn bán kính R thì : 4 y RJ 4   d) Biến dạng uốn trong mặt phẳng xy ( do { q}xz = { q2 , q6 , q8 , q12}T ) áp dụng 2 6 8 12 2 2 2 6z xy 3e 8 2 2 12 q q q q 12 6L 12 6L q q6L 4L 6L 2LEJK q12 6L 12 6LL q6L 2L 6L 4L                   (III.54) với 2z F J y dF  - mô men quán tính của mặt cắt ngang đối với đường trung hòa.  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 = 3a b 12  Nếu thiết diện tròn bán kính R thì : 4 z RJ 4   Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 82  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) trong hệ tọa độ tổng thể Ta có [ K' ](e) =      T(e) (e) (e)T K T 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 :  (e) [n] [0] [0] [0] [0] [n] [0] [0] T [0] [0] [n] [0] [0] [0] [0] [n]             (III.55) Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 83 với : [0] = 0 0 0 0 0 0 0 0 0          và [n] = x x x y y y y z z m n m n m 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' khi đó : x x ' y [n]. y ' z z '                   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 Nút T.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 Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 84 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 thanh L E G Jx Jy Jz x mx nx y my ny z mz nz 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 : Nút đầu : i Nút cuối : J T.T 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 : Ke(1 to 12, 1 to 12) as double Bước 6. Tính [K'](e) =  T(e)T . [K](e) . [T](e) được ghi vào mảng : Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 85 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. Bước 10. Giải hệ phương trình tính các chuyển vị kq ( giá trị các bậc tự do theo hệ tọa độ 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 Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 86 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 Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 87 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 Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 88 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 tg = tg + Ke(i,s) * T(s,j) next MaTranTG(i,j) = tg next Next For i =1 to 12 for j = 1 to 12 tg = 0 for s = 1 to 12 tg = tg + T(s,i) * MaTranTG(s,j) next KeT(i,j) = tg next Next End sub Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 89 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 Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012 90 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