Bài giảng Cơ học môi trường liên tục - Trần Văn Liên

v u v u bcbcbc ccc bbb A x v y u y v x u e k k j j i i kkjjii kji kji e ee e e xy yy xx e = ⎪⎪ ⎪⎪ ⎭ ⎪⎪ ⎪⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ = ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪ ⎬ ⎫ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ ∂ ∂+∂ ∂ ∂ ∂ ∂ ∂ = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ = 000 000 2 1 γ ε ε ε (9.4.7) Do đó, ma trận quan hệ biến dạng – chuyển vị nút Be có dạng ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ = kkjjii kji kji e e bcbcbc ccc bbb A B 000 000 2 1 (9.4.8) Từ bảng 9.2.2, ma trận các hằng số đμn hồi của bμi toán ứng suất phẳng lμ ( )⎟ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ −− = ν ν ν ν 100 01 01 1 2 1 2 EDe (9.4.9) Ma trận các hằng số đμn hồi De của bμi toán biến dạng phẳng t−ơng tự ma trận các hằng số đμn hồi của bμi toán ứng suất phẳng (9.4.9) với ν ννν −=−= 1;1 121 EE Do ma trận quan hệ biến dạng – chuyển vị nút Be của phần tử tam giác lμ hằng số nên biến dạng vμ ứng suất tại các điểm trong phần tử lμ không đổi. Vì vậy, ma trận độ cứng (9.2.14) của phần tử tam giác có dạng ee T ee V eee T ee BDBhAdVBDBK e == ∫ (9.4.10) trong đó h lμ độ dầy của phần tử, Ae lμ diện tích phần tử tam giác vμ các nút i, j, k có thứ tự ng−ợc chiều kim đồng hồ. Việc xác định ma trận khối l−ợng (9.2.13) hay ma trận cản (9.2.15) đòi hỏi phải tích phân khá phức tạp. Đối với các vật thể đμn hồi có các thμnh phần ứng suất vμ biến dạng thay đổi nhanh thì việc chọn phần tử tam giác dẫn đến số l−ợng phần tử hữu hạn khá lớn để đảm bảo độ chính xác cần thiết. Khi đó, phần tử tứ giác th−ờng đ−ợc sử dụng thay cho phần tử tam giác. Để quy đổi tải trọng trên mặt biên về các nút, ta sử dụng hệ tọa độ phụ dọc theo biên. Xét tr−ờng hợp tải trọng theo ph−ơng x phân bố tuyến tính trên biên jk với giá Cơ học môi tr−ờng liên tục Trần Văn Liên 256 trị kj xx qq ; tại các nút j vμ k. Trong hệ tọa độ phụ với gốc tại nút j, tọa độ của các điểm trên biên lμ ( ) ( ) L syyyy L sxxxx jkjjkj −+=−+= ; (9.4.11) với L lμ chiều dμi biên jk vμ Ls ≤≤0 . Tải trọng theo ph−ơng x phân bố tuyến tính trên biên jk đ−ợc xác định theo công thức ( ) ( ) L sqqqsq jkj xxxx −+= (9.4.12) Sử dụng (9.4.11), viết lại (9.4.5) d−ới dạng (9.2.5) ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ − −=⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛= LsLs LsLs NNN NNN N kji kji e 01000 00100 000 000 (9.4.13) khi đó Ne lμ ma trận các hμm dạng của phần tử tam giác. Sử dụng công thức (9.2.17) với hμm dạng Ne xác định theo (9.4.13), ta nhận đ−ợc tải trọng nút quy đổi ( ) ( ) ( ) ( )∫ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ++== L T xxxx x T eS LqqLqq dssqsNP kjkj e 0 0 6 2 0 6 2 00 (9.4.14) Công thức (9.4.14) t−ơng tự công thức quy đổi tải trọng tuyến tính hình thang (9.3.7). Đối với tr−ờng hợp tải trọng theo ph−ơng y trên biên jk vμ các tr−ờng hợp trên biên khác, ta có kết quả t−ơng tự. Tr−ờng hợp tải trọng P0 phân bố đều theo bề dầy tại tọa độ Laas <<= 0; , sử dụng hμm Dirac ( ) ( )asPsq −= δ0 , ta có kết quả t−ơng tự (9.3.8). Do ph−ơng của hệ tọa độ địa ph−ơng trong phần tử tam giác trùng với ph−ơng của hệ tọa độ tổng thể nên ta không cần thực hiện phép chuyển ma trận độ cứng Ke vμ véc tơ tải trọng Pe về hệ tọa độ tổng thể thông qua ma trận Te . Sau khi tìm đ−ợc các chuyển vị nút, ta xác định đ−ợc biến dạng vμ ứng suất trong các phần tử theo các công thức (9.2.6) vμ (9.2.7). Ví dụ 9.4.1: Tính ứng suất trong tấm ở trạng thái ứng suất phẳng có độ dầy h không đổi với liên kết, tải trọng nh− trên hình 9.4.2. Cho biết hệ số Poisson ν=0.25. Giải: Chia tấm thμnh hai phần tử I vμ II với các nút đ−ợc đánh số thứ tự (1), (2), (3), (4) nh− hình 9.4.2. Mỗi nút gồm hai chuyển vị u, v theo các ph−ơng trục x, y vμ véc tơ các chuyển vị nút cần tìm lμ ( )Tvuvuvuvu 44332211 . Cơ học môi tr−ờng liên tục Trần Văn Liên 257 Phần tử I gồm các nút 1(0,0), 2(a,0), 4(0,2a) có thứ tự ng−ợc chiều kim đồng hồ vμ các chuyển vị nút t−ơng ứng 442211 ;;;;; vuvuvu . Diện tích phần tử I theo (9.4.4) lμ 2 200 00 111 det 2 1 a a aAI =⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ = Xác định các hệ số của ma trận BI theo (9.4.3) aaxxcyyb xxcaayyb aaxxcaayyb =−=−==−=−= =−=−==−=−= −=−=−=−=−=−= 0;000 000;202 0;220 124214 412142 241421 do đó ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ −− − − = ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ = 012021 100010 000202 2 1000 000 2 1 442211 421 421 a bcbcbc ccc bbb A B I I Ma trận các hằng số đμn hồi D không thay đổi trong cả tấm vμ khi ν=0.25 ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ = ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ = 300 082 028 15 2 375.000 0125.0 025.01 15 16 EED Hình 9.4.2. 4 v1 v2 v3 2a 1 2 3 I II q 2q a q u1 u2 u3 u4 v4 qa/2 qa/2 4qa/3 5qa/3 x y Cơ học môi tr−ờng liên tục Trần Văn Liên 258 Từ đó ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ −− −− −− = ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ −− − − ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ = 036063 800484 20016216 15 012021 100010 000202 2 1 300 082 028 15 2 a E a EDBI Ma trận độ cứng của phần tử I xác định theo (9.4.10) lμ ⎟⎟ ⎟⎟ ⎟⎟ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎜⎜ ⎜⎜ ⎝ ⎛ −− −− −− −− −−−− −−−− = ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ −− −− −− ⎟⎟ ⎟⎟ ⎟⎟ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎜⎜ ⎜⎜ ⎝ ⎛ −− −− == 800484 036063 06120126 40032432 861242010 436321035 30 036063 800484 20016216 15 010 100 200 002 210 102 2 12 Eh a E a haDBBhAK I T III Có thể chứng minh rằng, tổng các phần tử theo một hμng hay theo một cột bất kỳ của ma trận độ cứng phần tử tam giác Ke với các hμm dạng (9.4.6) luôn bằng không. Cùng với tính chất đối xứng vμ tính chất xác định d−ơng của ma trận độ cứng, các tính chất nμy rất cần thiết khi kiểm tra kết quả tính toán. Phần tử II gồm các nút 2(a,0), 3(a,2a), 4(0,2a) có thứ tự ng−ợc chiều kim đồng hồ vμ các chuyển vị nút t−ơng ứng 443322 ;;;;; vuvuvu . Diện tích phần tử II theo (9.4.4) lμ 2 220 0 111 det 2 1 a aa aaAII =⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ = Xác định các hệ số của ma trận BII theo (9.4.3) 0;220 0;202 0;022 324324 423243 342432 =−=−=−=−=−= =−=−==−=−= −=−=−==−=−= aaxxcaayyb aaxxcaayyb aaxxcaayyb do đó Cơ học môi tr−ờng liên tục Trần Văn Liên 259 ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ −− − − = ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ = 202101 001010 020200 2 1000 000 2 1 443322 432 432 a bcbcbc ccc bbb A B II II Vì vậy ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ −− −− −− = ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ −− − − ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ = 606303 048480 01621620 15 202101 001010 020200 2 1 300 082 028 15 2 a E a EDBII Ma trận độ cứng của phần tử II xác định theo (9.4.10) lμ ⎟⎟ ⎟⎟ ⎟⎟ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎜⎜ ⎜⎜ ⎝ ⎛ −− −− −−−− −−−− −− −− = ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ −− −− −− ⎟⎟ ⎟⎟ ⎟⎟ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎜⎜ ⎜⎜ ⎝ ⎛ − − − − == 12012606 03243240 124201086 632103543 048480 606303 30 606303 048480 01621620 15 200 002 210 102 010 100 2 12 Eh a E a haDBBhAK II T IIIIII Theo ph−ơng pháp độ cứng trực tiếp (9.2.21), ta xác định đ−ợc ma trận độ cứng của cả hệ nh− sau ⎟⎟ ⎟⎟ ⎟⎟ ⎟⎟ ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎜⎜ ⎜⎜ ⎜⎜ ⎜ ⎝ ⎛ −−−− −−−− −−−− −−−− −−−− −−−− −−−− −−−− = 20012601084 03543210063 12420108600 63210354300 01084200126 10063035432 86001242010 43006321035 30 EhK Tải trọng quy đổi về các nút xác định theo công thức (9.4.14) ( ) TTyxyxyxyx qaqaqaqaPPPPPPPPP ⎟⎠⎞⎜⎝⎛ −−−−== 202340350044332211 Cơ học môi tr−ờng liên tục Trần Văn Liên 260 Từ điều kiện biên, ta có 0;0;0;0 4211 ==== uvvu . Theo (9.2.23) ta cần: - Loại bỏ ẩn 4211 ;;; uvvu trong véc tơ ẩn số. - Loại bỏ các hμng vμ các cột 1, 2, 4, 7 trong ma trận độ cứng K. - Loại bỏ các hμng 1, 2, 4, 7 trong véc tơ tải trọng nút P. Từ đó ta nhận đ−ợc ph−ơng trình rút gọn ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ −= ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ ⎟⎟ ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎜ ⎝ ⎛ −− −− −− −− 3 3 8 10 6 2012610 1220106 610353 106335 30 4 3 3 2 qa v v u u Eh Giải hệ ph−ơng trình nμy, ta đ−ợc ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ −= ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ 14 16 13 17 12 4 3 3 2 Eh qa v v u u Véc tơ chuyển vị nút viết đầy đủ lμ ( ) ( )TT Eh qavuvuvuvu 140161301700 1244332211 −= Xác định ứng suất trong từng phần tử theo công thức (9.2.7): - Phần tử I có ( ) ( )TITxyyyxxI vuvuvuDB 442211== σσσσ do đó ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ −= ⎪⎪ ⎪⎪ ⎭ ⎪⎪ ⎪⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ ⎟⎠ ⎞⎜⎝ ⎛− ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ −− −− −− = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ 0 1 35 14 0 0 17 0 0 12 036063 800484 20016216 15 h q Eh qa a E xy yy xx σ σ σ - Phần tử II có ( ) ( )TIITxyyyxxII vuvuvuDB 443322== σσσσ do đó Cơ học môi tr−ờng liên tục Trần Văn Liên 261 ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ −= ⎪⎪ ⎪⎪ ⎭ ⎪⎪ ⎪⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ ⎟⎠ ⎞⎜⎝ ⎛− ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ −− −− −− = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ 0 1 34 14 0 16 13 0 17 12 606303 048480 01621620 15 h q Eh qa a E xy yy xx σ σ σ Giá trị ứng suất tại các nút đ−ợc tính theo công thức ∑= r ri n σσ 1 trong đó n lμ số phần tử cùng chung nút i, r lμ tên phần tử có nút i. Từ đó ( ) ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ −= ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ −=+== ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ −= 0 1 34 ; 0 1 23 2 1; 0 1 35 3421 h q h q h q III σσσσσσ 9.5. Ph−ơng pháp phần tử hữu hạn trong bμi toán tấm chịu uốn Giả thiết độ võng w(x,y) của tấm lμ nhỏ. Nếu giả thiết rằng độ võng nμy độc lập với chuyển vị của các điểm thuộc mặt trung bình của tấm u0(x,y) vμ v0(x,y), tức lμ độc lập với trạng thái ứng suất trong mặt trung bình, thì việc xác định độ võng nμy trở nên đơn giản hơn. Tuy nhiên khi áp dụng ph−ơng pháp phần tử hữu hạn cho bμi toán tấm có một số khó khăn: - Nếu trong các bμi toán thanh vμ bμi toán tấm phẳng chịu tải trọng nằm trong mặt phẳng trung bình, nghiệm chuyển vị d−ới dạng đa thức lμ duy nhất thì ng−ợc lại trong bμi toán tấm chịu tải trọng ngang, nghiệm chuyển vị d−ới dạng đa thức lại không phải lμ duy nhất dẫn đến ma trận độ cứng thu đ−ợc khác nhau. Khi đó nảy sinh vấn đề ma trận độ cứng nμo lμ xấp xỉ tốt nhất với nghiệm chính xác. - Yêu cầu về tính liên tục của các đạo hμm của độ võng tại các đ−ờng biên của phần tử để đảm bảo khi ghép các phần tử không xuất hiện các vết gãy đã dẫn đến những khó khăn về mặt toán học. Thực tế tính toán cho thấy, nếu ta chấp nhận điều kiện liên tục của độ võng vμ bỏ qua điều kiện liên tục của các đạo hμm hμm độ võng tại các đ−ờng biên của phần tử (đ−ợc gọi lμ các hμm độ võng “không t−ơng thích”) thì kết quả thu đ−ợc sẽ hội tụ về nghiệm chính xác khi chia tăng số phần tử hữu hạn. Để có đ−ợc độ chính xác mong muốn, nếu ta chọn hμm độ võng thỏa mãn đồng thời điều kiện liên tục của độ võng vμ điều kiện liên tục của các đạo hμm hμm độ võng tại các đ−ờng biên của phần tử Cơ học môi tr−ờng liên tục Trần Văn Liên 262 (đ−ợc gọi lμ các hμm độ võng “t−ơng thích”) thì số phần tử hữu hạn phải chia sẽ ít hơn so với việc chọn hμm độ võng không t−ơng thích. Khi tính toán tấm bằng ph−ơng pháp phần tử hữu hạn, ng−ời ta th−ờng chia tấm thμnh các phần tử hình chữ nhật vμ phần tử hình tam giác. Tại mỗi nút của phần tử có ba chuyển vị suy rộng: một chuyển vị theo trục z vμ hai góc xoay quanh các trục x vμ y (hình 9.5.1). D−ới đây lμ hai tr−ờng hợp chọn hμm dạng độ võng cho phần tử tấm hình chữ nhật: - Hμm dạng độ võng “không t−ơng thích” Ni , i=1,2,..,12 với byax == ηξ ; ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )ηξηξξξξηξηξη ηξξηηξξηξξηξηξηξ ξηξηξηηξξηηξξη ηξηξξηξηξηηη ηξξηηξξηηξηη ηξξηξξξηξ ξηηξηηξηη ξηηξηξηηξξηξηξ 3232 12 32 11 233232 10 32 9 32 8 3322 7 32 6 3232 5 332232 4 3322 3 3322 2 33322322 1 ;2 322323; ;2233 2; 223323 22 22 222332331 +−−=+−= −++−−+=−= +−=−−++−= −+−=−++−= ++−−−+= +−−++−= −++−−= −−++++−−−= aNbN NaN bNN aNbN N aN bN N (9.5.1) - Hμm dạng độ võng “t−ơng thích” Ni , i=1,2,..,12 với byax == ηξ ; ( )( ) ( )( ) ( )( ) ( ) ( ) ( )( ) ( )( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( )( ) aNbN NaN bNN aNbN NaN bNN 22 12 22 11 22 10 22 9 22 8 22 7 22 6 22 5 22 4 22 3 22 2 22 1 1211;123 12123;231 123;2323 231;1121 23121;1211 1121;121121 ηηξξηηξξ ηηξξηηξξ ηηξξηηξξ ηηξξηηξξ ηηξξηηξξ ηηξξηηξξ −+−=−−= −+−=−−= −−−=−−= −−−=−−−−= −−+=−+−−= −−+=−+−+= (9.5.2) Từ bảng 9.2.1, ma trận quan hệ biến dạng – chuyển vị nút của phần tử tấm có dạng ( )1221 BBBBe L= với các cột Bi lμ T iii i yx N y N x NB ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ∂∂ ∂ ∂ ∂ ∂ ∂= 2 2 2 2 2 2 (9.5.3) vμ Ni lμ các hμm dạng “không t−ơng thích” (9.5.1) hay lμ các hμm dạng “t−ơng thích” (9.5.2). b q8 Hình 9.5.1. y q7 q9 q10 q12 q11 q4 q6 q5 q1 q3 q2 a x z Cơ học môi tr−ờng liên tục Trần Văn Liên 263 Từ bảng 9.2.2, ma trận các hằng số đμn hồi của phần tử tấm chịu uốn có dạng ( )⎟ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ −− = ν ν ν ν 100 01 01 )1(12 2 1 2 EDe (9.5.4) Từ (9.2.13) – (9.2.14), ta xác định đ−ợc ma trận khối l−ợng Me với các thμnh phần dxdyNNhm a b jiij ∫ ∫= 0 0 ρ (9.5.5) vμ ma trận độ cứng Ke của phần tử với các thμnh phần ( ) ( ) dxdy x N y N y N x N yx N yx N y N x N y N x N Ehk a b jijiji jjii ij ∫ ∫ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ∂ ∂ ∂ ∂−∂ ∂ ∂ ∂−∂∂ ∂ ∂∂ ∂− +⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ∂ ∂+∂ ∂ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ∂ ∂+∂ ∂ −= 0 0 2 2 2 2 2 2 2 222 2 2 2 2 2 2 2 2 2 3 2 1 2 112 112 νν (9.5.6) trong đó h lμ độ dầy phần tử tấm. Các tích phân (9.5.5), (9.5.6) có thể thực hiện dễ dμng nhờ các phần mềm symbolic nh− Maple hay MatLab. Ví dụ 9.5.1: Xác định 4 tần số riêng đầu tiên của tấm kích th−ớc aìa với các cạnh tựa tự do khi sử dụng các hμm dạng độ võng “không t−ơng thích” (9.5.1). Giải: Từ ph−ơng trình tần số ( ) 0det 2 =− MK ω , ta tìm đ−ợc các tần số riêng có dạng ( )2 3 4 112 1 νραω −= Eh haii với giá trị αi theo từng ph−ơng án chọn l−ới phần tử hữu hạn khác nhau đ−ợc dẫn ra ở bảng 9.5.1. Rõ rμng lμ khi chia tăng số phần tử hữu hạn thì tần số riêng thu đ−ợc cμng gần với nghiệm chính xác hơn. Bảng 9.5.1 L−ới chia Tần số thứ nhất Tần số thứ hai Tần số thứ ba Tần số thứ t− 2ì2 20.471 53.071 56.932 127.410 4ì4 19.848 51.448 54.020 91.018 8ì8 19.754 50.913 52.721 85.314 Nghiệm chính xác 19.739 49.348 49.348 78.957 Cơ học môi tr−ờng liên tục Trần Văn Liên 264 9.6. Ph−ơng pháp ma trận độ cứng động lực Khi ứng dụng ph−ơng pháp phần tử hữu hạn vμo các bμi toán động lực học, chỉ có một chỗ duy nhất mμ ta phải xấp xỉ tr−ờng chuyển vị trong phần tử bằng tr−ờng chuyển vị tĩnh, tức lμ đã bỏ qua yếu tố động lực học của tr−ờng chuyển vị. Nếu ta chọn các hμm dạng của phần tử hữu hạn lμ tr−ờng chuyển vị động thỏa mãn ph−ơng trình cân bằng động thì ph−ơng pháp phần tử hữu hạn không còn lμ một ph−ơng pháp gần đúng mμ lμ một ph−ơng pháp chính xác. Ph−ơng pháp ma trận độ cứng động lực đã ra đời trên cơ sở ý t−ởng nμy. Để chọn đ−ợc các hμm dạng chuyển vị động, ta phải xét bμi toán cân bằng động của phần tử hữu hạn trong miền tần số, tức lμ xét chuyển động với biên độ phức phụ thuộc vμo tần số. Sau đó, các b−ớc thực hiện của ph−ơng pháp ma trận độ cứng động lực về hình thức không khác gì ph−ơng pháp phần tử hữu hạn thông th−ờng. 9.6.1. Khái niệm ma trận độ cứng động lực Trong miền tần số, ph−ơng trình cơ bản của ph−ơng pháp phần tử hữu hạn (9.1.1) có dạng ( ) ( ) ( )ωωω PUK ˆˆˆ = (9.6.1) trong đó ( ) MCiKK 2ˆ ωωω −+= (9.6.2) Ma trận )(ˆ ωK đ−ợc gọi lμ ma trận độ cứng động lực, lμ hμm của tần số vμ phụ thuộc tuyến tính vμo các ma trận khối l−ợng, cản vμ độ cứng. Véc tơ )(ˆ),(ˆ ωω PU lần l−ợt đ−ợc gọi lμ biên độ phức của chuyển vị nút vμ tải trọng quy về nút. Nh− vậy, bμi toán động lực học của hệ đμn hồi dẫn đến giải hệ ph−ơng trình đại số đối với chuyển vị trong miền tần số, nghĩa lμ cho tất cả các tần số trong một dải tần nμo đó. Công việc chính của ph−ơng pháp ma trận độ cứng động lực lμ xây dựng ma trận độ cứng động lực )(ˆ ωK vμ véc tơ tải trọng quy về nút )(ˆ ωP . Các bμi toán cơ bản trong phân tích hệ đμn hồi sử dụng ph−ơng pháp ma trận độ cứng động lực bao gồm: a) Bμi toán phân tích tĩnh có dạng )0(ˆˆ)0(ˆ 0 PUK = (9.6.3) kết quả cho ta chuyển vị tĩnh của nút 0Uˆ . b) Bμi toán dao động riêng có dạng 0)(ˆ =φωK (9.6.4) Cơ học môi tr−ờng liên tục Trần Văn Liên 265 trong đó các tần số riêng ωj đ−ợc xác định từ ph−ơng trình 0)(ˆdet =ωK (9.6.5) Các dạng riêng φj t−ơng ứng với tần số riêng ωj đ−ợc tìm từ (9.6.4) cùng với điều kiện chuẩn hoá 1=jφ . c) Bμi toán dao động c−ỡng bức với kích động điều hoμ lμ bμi toán tổng quát (9.6.1) nêu trên. d) Bμi toán ổn định theo tiêu chuẩn ổn định động lực học đ−a về bμi toán xác định tải trọng tới hạn P sao cho 0),(ˆdet =PK ω (9.6.6) Tải trọng tới hạn P đ−ợc xác định từ điều kiện tần số dao động riêng ω lμ một số phức có phần ảo âm. Khi đó hệ sẽ mất ổn định do sự tăng dần biên độ chuyển động của các chuyển động bé của hệ ở lân cận vị trí cân bằng. Khai triển ma trận độ cứng động lực ra chuỗi Taylor theo tần số ω cho ta )()(ˆ 32 ωωωω oMCiKK +−+= (9.6.7) trong đó K, M, C lμ các ma trận của ph−ơng pháp phần tử hữu hạn thông th−ờng. Nh− vậy, ph−ơng pháp ma trận độ cứng động lực bao hμm ph−ơng pháp phần tử hữu hạn thông th−ờng nh− lμ một tr−ờng hợp riêng. 9.6.2. Ph−ơng pháp ma trận độ cứng động lực trong bμi toán thanh 9.6.2.1. Phần tử thanh chịu kéo nén hay chịu xoắn Xét một phần tử thanh có diện tích tiết diện A, môđun đμn hồi E, mật độ khối l−ợng ρ, chiều dμi L chịu tải trọng phân bố dọc trục ( )txqe ,* có chiều d−ơng h−ớng theo trục x (hình 9.6.1). Ph−ơng trình dao động dọc của thanh có dạng ( ) ( ) ( ) ( ) ( )⎥⎦ ⎤⎢⎣ ⎡ ∂ ∂+∂ ∂=+⎥⎦ ⎤⎢⎣ ⎡ ∂∂ ∂+∂ ∂ t txu t txuAtxq tx txu x txuEA eeeee ,, , ,, 22 2 * 2 3 12 2 μρμ (9.6.8) trong đó μ1 lμ hệ số cản nhớt của vật liệu, μ2 lμ hệ số cản của môi tr−ờng. Bằng cách đặt ( ) ( ) ( ) ( ) 1 ,, ,, * −= = Φ= i exqtxq extxu ti ee ti ee ω ω ω ω (9.6.9) với Φe(x,ω) vμ qe(x,ω) lμ biên độ của chuyển vị j y Ue2 Re1 Re2 x Ue1 Hình 9.6.1. k Cơ học môi tr−ờng liên tục Trần Văn Liên 266 dọc trục vμ tải trọng dọc trên thanh, ta thu đ−ợc ph−ơng trình ( ) ( ) ( )ωωλω ,~,, 22 2 xqx dx xd ee e =Φ+Φ (9.6.10) trong đó ( ) ⎟⎠ ⎞⎜⎝ ⎛ −=−= ω μρωλω 222 1ˆ;ˆ ,~ i EAE xqq ee (9.6.11) λ lμ tham số động lực học; ω lμ tần số vòng (rad/giây); nếu λ=0 tức ω=0 ta có tr−ờng hợp biến dạng tĩnh; nếu ω≠0 thì Φe(x,ω) lμ biên độ chuyển vị động. Đại l−ợng ( )ωμ11ˆ iEE += (9.6.12) lμ môđun đμn hồi phức, d−ới đây để cho dễ theo dõi ta vẫn dùng ký hiệu E nh− khi không có cản. Ký hiệu véc tơ chuyển vị nút ( )Teee UUU 21ˆ = với ( ) ( ) 21 ;0 eeee ULU =Φ=Φ (9.6.13) vμ véc tơ ứng lực tại tiết diện thanh sát nút ( )Teee RRR 21ˆ = với 21 0 ; e Lx e e x e R dx dEAR dx dEA =Φ=Φ− == (9.6.14) Nghiệm tổng quát của ph−ơng trình (9.6.10) có dạng ( ) ( )[ ] ( )∫ −+−=Φ x eeee dqxREAxxUx 011 ,~sin 1sincos τωττλλλ λλ (9.6.15) Từ các hệ thức (9.6.13) – (9.6.15), ta nhận đ−ợc ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛+⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ −− −=⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ − 2 1 1 1 1 2 2 cossin sin)(cos e e e e e e q q R U LLEA LEAL R U λλλ λλλ (9.6.16) trong đó ( )[ ] ( ) ( )[ ] ( )∫∫ −=−= L eeL ee dqLEAqdqLq 0 2 0 1 ,~cos;,~sin 1 τωττλτωττλλ (9.6.17) Sử dụng (9.6.16), ta nhận đ−ợc ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ − −+⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛=⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ − − 2 1 2 1 2 1 1sincos 0sin cos1 1cos sin e e e e e e q q LLEA LEA R R U U L L L EA λλλ λλ λ λ λ λ Cơ học môi tr−ờng liên tục Trần Văn Liên 267 hay lμ ( ) eeee PRUK ˆˆˆˆ +=ω (9.6.18) Đến đây ta nhận đ−ợc ma trận độ cứng động lực eKˆ vμ véc tơ tải trọng quy về nút ePˆ cho phần tử thanh chịu kéo nén có dạng ( ) ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ − −= ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ − −=Κ 2 1 1sincos 0sinˆ cos1 1cos sin ˆ e e e e q q LLEA LEA P L L L EA λλλ λλ λ λ λ λω (9.6.19) Nếu cho tần số ω tiến đến 0, ta nhận đ−ợc ma trận độ cứng Ke vμ ma trận khối l−ợng Me của ph−ơng pháp phần tử hữu hạn thông th−ờng (9.3.5) ( ) ... cos1 1cos sin ˆ 2 +−=⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ − −=Κ eee MKL L L EA ωλ λ λ λω Khi biên độ tải trọng dọc lμ hằng số ( ) ee qxq =ω, , tải trọng quy về nút có dạng ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ +=⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛= 1 1 cos1 sin ˆ ˆˆ 2 1 L Lq P PP e e e e λ λ λ (9.6.20) Nếu cho tần số ω tiến đến 0, ta nhận đ−ợc véc tơ tải trọng quy về nút của ph−ơng pháp phần tử hữu hạn thông th−ờng (9.3.7). Biên độ lực dọc trong thanh đ−ợc xác định theo công thức t−ơng tự (9.3.13) ( ) ( ) dx xdAExN e ωω ,ˆ, Φ= (9.6.21) Tr−ờng hợp phần tử thanh chịu xoắn, ma trận độ cứng động lực vμ véc tơ tải ngoμi vẫn có dạng nh− (9.6.19) nh−ng thay cho E vμ A lμ G vμ I. 9.6.2.2. Phần tử dầm hai chiều chịu uốn Xét một phần tử dầm có diện tích tiết diện A, mômen quán tính tiết diện Iz , môđun đμn hồi E, mật độ khối l−ợng ρ, chiều dμi L chịu uốn bởi tải trọng phân bố ( )txqe ,* có chiều d−ơng h−ớng theo trục y trong mặt phẳng xy của hệ tọa độ địa ph−ơng (hình 9.6.2). Ph−ơng trình dao động uốn của dầm có dạng ( ) ( ) ( ) ( ) ),(,,,, *22 2 4 5 14 4 txq t txv t txvA tx txv x txvEI eeeeez =⎥⎦ ⎤⎢⎣ ⎡ ∂ ∂+∂ ∂+⎥⎦ ⎤⎢⎣ ⎡ ∂∂ ∂+∂ ∂ μρμ (9.6.22) với μ1 lμ hệ số cản nhớt của vật liệu, μ2 lμ hệ số cản của môi tr−ờng. Cơ học môi tr−ờng liên tục Trần Văn Liên 268 Bằng cách đặt ( ) ( ) ( ) ( ) 1 ,, ,, * −= = Φ= i exqtxq extxv ti ee ti ee ω ω ω ω (9.6.23) với Φe(x,ω) vμ qe(x,ω) lμ biên độ của chuyển vị ngang vμ tải trọng ngang, ta thu đ−ợc ph−ơng trình ( ) ( ) ( )ωωλω ,~,, 44 4 xqx dx xd ee e =Φ−Φ (9.6.24) trong đó ( ) ( ) ⎟⎠ ⎞⎜⎝ ⎛ −== ω μρωλωω 224 1ˆ;ˆ ,,~ i IE A IE xqxq zz e e (9.6.25) λ lμ tham số động lực học; ω lμ tần số vòng (rad/giây); nếu λ=0 tức ω=0 ta có tr−ờng hợp biến dạng tĩnh; nếu ω≠0 thì Φe(x,ω) lμ biên độ chuyển vị động. Eˆ lμ môđun đμn hồi phức theo (9.6.12), d−ới đây để cho dễ theo dõi ta vẫn dùng ký hiệu E nh− khi không có cản. Ký hiệu véc tơ chuyển vị nút ( )Teeeee UUUUU 4321ˆ = với ( ) ( ) ( ) ( )LULUUU eeeeeeee Φ′=Φ=Φ′=Φ= 4321 ; ; 0 ; 0 (9.6.26) vμ véc tơ ứng lực tại tiết diện thanh sát nút ( )Teeeee RRRRR 4321ˆ = với ( ) ( ) ( ) ( )LEIRLEIR EIREIR ezeeze ezeeze Φ ′′=Φ ′′′−= Φ ′′−=Φ ′′′= 43 21 ; 0; 0 (9.6.27) Nghiệm tổng quát của ph−ơng trình (9.6.24) có dạng ( ) ( ) ( ) ( ) ( ) ( )[ ] ( )∫ −+ −++=Φ x e e z e z eee dqxK R EI xKR EI xKUxKUxKx 0 43 22 3 13 4 2 2 11 ,~1 τωττλλ λ λ λ λ λ λλ (9.6.28) với L] ,0[x ∈ vμ ( ) ( ) ( ) ( )xKxKxKxK 4321 ;; ; lμ những hμm Krylov ( ) ( ) ( ) ( ) 2 sinsinh; 2 coscosh; 2 sinsinh; 2 coscosh 4321 xxxKxxxKxxxKxxxK −=−=+=+= (9.6.29) Từ các hệ thức (9.6.26) – (9.6.28), ta nhận đ−ợc j k x Re4 Ue2 Ue1 Re1 Re2 Ue4 Ue3 Re3 Hình 9.6.2. y Cơ học môi tr−ờng liên tục Trần Văn Liên 269 ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ + ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ ⎟⎟ ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎜ ⎝ ⎛ − −−− − − = ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ 4 3 2 1 2 1 2 1 1243 2 413 2 2 3 2 2 314 2 3 3 421 4 3 4 3 )()()()( )()()()( )()()()( )()()()( e e e e e e e e zz zz zz zz e e e e q q q q R R U U LKLKLKEILKEI LKLKLKEILKEI EILKEILKLKLK EILKEILKLKLK R R U U λλλλλλλ λλλλλλλ λλλλλλλ λλλλλλλ (9.6.30) trong đó ( )[ ] ( ) ( )[ ] ( ) ( )[ ] ( ) ( )[ ] ( ) ,~;,~ ,~1;,~1 0 24 0 13 0 322 0 431 ∫∫ ∫∫ −=−−= −=−= L e z e L eze L ee L ee dqLKEIqdqLKEIq dqLKqdqLKq τωττλλτωττλ τωττλλτωττλλ (9.6.31) Sử dụng (9.6.30) để biểu diễn các lực nút qua chuyển vị nút t−ơng tự (9.6.18), ta có ( ) eeeeeee PRqHRUK ˆˆˆˆˆ +=+=ω (9.6.32) trong đó ( )Teeeee qqqqq 4321= , ( )ωeKˆ lμ ma trận độ cứng động lực của phần tử dầm chịu uốn có cấp 4ì4 đối xứng với các thμnh phần ( ) ( ) ( ) ( ) ( ) LL LLEIk LL LLLLEIkk LL LLEIkk LL LLEIk LL LLEIkk LL LLLLEIkk zz zz zz λλ λλλ λλ λλλλλ λλ λλλ λλ λλλ λλ λλλ λλ λλλλλ coscosh1 sinsinhˆ; coscosh1 cossinhsincoshˆˆ coscosh1 coscoshˆˆ; coscosh1 sinsinhˆ coscosh1 sinsinhˆˆ; coscosh1 cossinhsincoshˆˆ 244422 2 2314 3 13 2 3412 3 3311 − −=− −== − −=−=− +−= −=−=− +== (9.6.33) vμ He lμ ma trận có cấp 4ì4 với các thμnh phần 1;0;ˆ;ˆ;0;1;ˆ;ˆ 0;ˆ;ˆ;0;ˆ;ˆ 444322421241343312321131 242324221421141314121311 −===−==−=−== ===−===== hhkhkhhhkhkh hhkhkhhhkhkh (9.6.34) Nếu cho tần số ω tiến đến 0, ta nhận đ−ợc ma trận độ cứng Ke vμ ma trận khối l−ợng Me của ph−ơng pháp phần tử hữu hạn thông th−ờng (9.3.18) ( ) )(ˆ 42 ωωω oMKK eee +−= Khi biên độ tải trọng ngang lμ hằng số ( ) ee qxq =ω, , tải trọng quy về nút có dạng ( )Teeeee PPPPP 4321 ˆˆˆˆˆ = với LL LLLLqPP LL LLLLLLqPP e ee e ee λλ λλλλ λ λλ λλλλλλ