Đề tài Dao động uốn của tấm mỏng hình chữ nhật với điều kiện biên không đối xứng

kính cong : Thay vào (1.25), ta có : (1.27) Biểu thức (1.27) cho ta biết độ biến đổi hình dáng của tấm nếu đã Mx và My đã biết. Nếu Mx = 0 thì , hoặc My = 0 thì , do đó rx và ry trái dấu nhau, ta gọi đó là uốn lệch pha ( anticlastic). Nếu Mx = My = M thì , sự uốn theo mặt cầu gọi là uốn đồng pha ( syclastic) và . Mx My Hình 1.8 – Uốn lệch pha 1.2.5.2 Tấm chịu uốn và xoắn đồng thời. a) Mômen chính và độ cong chính. Thông thường, mômen uốn tác dụng lên tấm sẽ không nằm trong mặt phẳng vuông góc với các cạnh. Vì thế, mômen uốn có thể tách thành các thành phần đơn giản gồm thành phần tiếp tuyến và thành phần pháp tuyến. Thành phần pháp tuyến Mx và My đã đề cập ở phần trên, còn các thành phần tiếp tuyến Mxy và Myx ( cũng là mômen đơn vị ) gây ra sự xoắn của tấm đối với các trục song song với trục x và y. Chỉ số thứ nhất của mômen xoắn chỉ phương pháp tuyến của mặt đang xét, chỉ số thứ hai của mômen xoắn chỉ độ dài cạnh song song trục nào. Quy ước dấu : Mx , My > 0 nếu làm căng thớ ở z > 0; Mxy ,Myx > 0 nếu quay ngược chiều kim đồng hồ khi nhìn dọc trục của nó theo hướng song song với hướng dương của trục tương ứng x hoặc y. Trong hình vẽ tất cả các cường độ mômen đều dương. Mxy Myx( = - Mxy ) My Mx Mx x Mxy y Myx My z Hình 1.9 – Uốn và xoắn đồng thời tấm phẳng D Mxy Myx My Mx a Mx x Mxy y F Myx My z Mt Mn A a Mx C Mxy Mxy My B Hình 1.10 – Mômen trong mặt phẳng bất kì Vì Mxy và Myx sinh ra ứng suất tiếp txy = - tyx nên Mxy = - Myx. Trên mặt phẳng chéo FD có các mômen tác dụng là mômen tiếp Mt và mômen pháp Mn . Chúng ta có thể biểu diễn các mômen này qua các Mx, My, Mxy nhờ vào các phương trình cân bằng phân tố tam giác ABC ( hình 1.10 ). Trong mặt phẳng vuông góc với AC có : ị Mn = Mx cos2a + Mysin2a - Mxysin2a (1.28) Tương tự cho cân bằng trong mặt phẳng song song với AC, cũng có : (1.29) Từ (1.28) và (1.29) ta nhận thấy có hai giá trị của a, chênh nhau 900 thoả mãn và : (1.30) Thay (1.30) vào (1.28), Mn có hai giá trị max và min, ta có : (1.31) Các mômen chính gây ra các độ cong chính tương ứng. Trong trường hợp tấm chịu uốn và xoắn thuần tuý trong đó Mx, My và Mxy không đổi trong suốt chiều dài các cạnh của tấm thì các mômen chính là đại lượng lớn nhất và nhỏ nhất trong tấm và được tính theo công thức (1.31). Từ đó ta thấy không có ứng suất cắt trên mặt phẳng đó và các ứng suất pháp tương ứng ( phụ thuộc vào z và Mx , My , Mxy) là ứng suất lớn nhất và nhỏ nhất trong tấm. b) Liên hệ giữa mômen xoắn và chuyển vị w. Trở lại tấm chịu tải như hình 1.10 a), chúng ta đã thành lập được mối liên hệ giữa cường độ mômen uốn Mx và My với chuyển vị w của tấm cho bởi các công thức (1.27). Tiếp theo chúng ta sẽ tìm mối liên hệ giữa mômen xoắn Mxy và chuyển vị w. Từ định lý siêu định vị chúng ta có thể xét riêng Mxy tác dụng bỏ qua Mx và My. Như ta đã biết Mxy bị chống lại bởi hệ các cặp ứng suất bù nhau trong mặt phẳng thẳng đứng của phân tố được lấy ra suốt chiều dày của tấm và song song với các trục x và y. B A E z Xét một phân tố thuộc phân tố đang xét ( hình 1.11 ). Các ứng suất tiếp bù nhau trên lân cận của phân tố cách một khoảng z phía dưới mặt trung hoà là txy. Cho nên trên mặt ABCD có : dx dy và trên mặt ADEF có : Mxy Mxy dz txy txy F C Hình 1.11 – ứng suất tiếp do mômen xoắn gây ra (1.32) Trong đó : G : mô đun đàn hồi. gxy: độ dãn dài góc. Chúng ta có mối quan hệ giữa biến dạng góc gxy và các chuyển vị cho bởi công thức : Ta cần biểu diễn gxy qua chuyển vị w của tấm. Lấy một phân tố có chiều dày của tấm sẽ bị quay một góc bằng và trong các mặt phẳng xz và yz. Xét sự quay của một phân tố trong mặt phẳng xz của một điểm ở phía dưới mặt trung hoà một khoảng z là : Tương tự, chuyển vị theo phương y là : z O x h dz z -u Hình 1.12 – Góc xoay của phân tố Thay u, v vào biểu thức gxy ta có : (1.33) Từ (1.32) : Thay cho ta : Nhân cả tử và mẫu của biểu thức với (1- m) thì : (1.34) Biểu thức (1.27) và (1.34) cho ta mối quan hệ giữa mômen uốn và xoắn với chuyển vị của tấm và chúng cũng tương đương với mối quan hệ giữa mômen uốn và độ cong trong dầm đơn. 1.2.5.3 Tấm chịu tải trọng phân bố vuông góc với mặt tấm. q y x z Hình 1.13 – Tấm chịu tải trọng phân bố vuông góc Mối quan hệ giữa mômen uốn và xoắn với chuyển vị của tấm sẽ được sử dụng trong việc thiết lập phương trình vi phân tổng quát cho bài toán tấm hình chữ nhật mỏng chịu tải trọng phân bố vuông góc với mặt tấm (hình 1.13) Trong trường hợp tổng quát tải trọng phân bố có thể thay đổi theo một hàm phụ thuộc vào x và y. Chúng ta vẫn giả thiết mặt trung gian trùng với mặt trung hoà, mặt cắt ngang vẫn phẳng trong suốt quá trình biến dạng. Giả thiết sau chỉ rõ sự mâu thuẫn trong lý thuyết, đó là nếu mặt cắt ngang vẫn phẳng trong quá trình biến dạng thì các biến dạng góc gxz và gyz phải bằng 0, trong khi tải trọng phân bố gây các lực cắt vuông góc mặt tấm vẫn tồn tại ( và là ứng suất đã biết ). Vì thế chúng ta giả thiết rằng mặc dù và là không đáng kể nhưng các lực cắt tương ứng lại giống với độ lớn của tải trọng phân bố q và các mômen Mx , My và Mxy. Giả thiết này tương tự trong lý thuyết dầm đơn mà trong đó biến dạng góc được bỏ qua. Xét một phân tố của tấm như hình 1.14, chịu mômen uốn và xoắn như mô tả lần trước, và các lực Qx và Qy tác dụng lên một đơn vị chiều dài trong từng mặt phẳng vuông góc với trục x và y. Thay đổi của ứng suất cắt txz và tyz theo các cạnh nhỏ dx và dy của phân tố là không đáng kể và hợp lực của lực cắt Qxdy và Qydx đặt ở trọng tâm các mặt phẳng của phân tố. Ta có : và (1.35) Phương trình cân bằng lực của phân tố theo trục Oz và giả thiết là trọng lượng riêng của tấm được kể trong q, có dạng : (1.36) Lấy mômen với trục x : Đơn giản phương trên và bỏ qua các đại lượng vô cùng bé bậc cao, ta được: (1.37) Tương tự lấy mômen đối với trục y, ta có : (1.38) 1.2.5.4 Tấm chịu uốn đồng thời với lực tác dụng trong mặt phẳng trung gian của tấm. Trong các trường hợp đã xét chúng ta giả thiết rằng các mặt trung hoà của tấm không có ứng suất. Nếu trong mặt phẳng của tấm thêm vào tải trọng kéo, nén hoặc tiếp thì sẽ gây ra ứng suất trong mặt trung gian và nếu như nó đủ lớn sẽ ảnh hưởng đến sự uốn của tấm. Xét một phân tố nhỏ dxdy thuộc mặt phẳng trung gian của tấm mỏng, vị trí chuyển vị của nó như hình 1.15. Chiều và độ lớn của lực trên một đơn vị dài được tạo bởi tải trọng trong mặt phẳng tấm là Nx, Ny , và Nxy và quy ước dấu của nó như đối với ứng suất pháp và ứng suất tiếp. Ngoại lực theo phương vuông góc với tấm không tham gia vào phương trình cân bằng chiếu theo phương x và y, do đó : dx Mxy dy q Qy My Mxy Qx Mx Hình 1.14 – Trạng thái chịu lực của phần tử tấm Theo phương x : do w nhỏ nên nhỏ, suy ra : cos ằ 1, cos ằ 1 Vậy ta có : = 0 (1.39) O dx x Nx- O Nyx Ny x Nxy y Nx Hình 1.15 – Lực tác dụng trong mặt phẳng tấm Tương tự theo phương y : = 0 (1.40) Những phương trình này hoàn toàn độc lập với ba phương trình cân bằng ( đã nêu trong 2.5.3) và vì thế ta có thể sử dụng chúng một cách độc lập. Khi xét hình chiếu của các lực lên trục z, ta phải kể đến ảnh hưởng của các lực Nx, Ny, Nxy và Nyx do khi tấm bị uốn mà có những góc nhỏ giữa các phương của các lực này. Thành phần hình chiếu lên trục z do Nxy là : ( bỏ qua vô cùng bé bậc cao ) Tương tự ta có thành phần hình chiếu lên trục z do Nyz là : Thành phần hình chiếu lên trục z do Nx là : Tương tự do Ny là : Do đó tổng hợp lực do Nx, Ny, Nxy,( Nxy = Nyx ) chiếu lên phương z là : Sử dụng (1.39) và (1.40) ta có: Cộng hợp lực này với tải trọng qdxdy tác động vào phân tố và dùng phương trình (1.36) ta có phương trình cân bằng dưới đây: (1.41) Chương 2 Dao động uốn của tấm mỏng 2.1 Thiết lập phương trình uốn của tấm mỏng. Xét dao động uốn của tấm mỏng đồng chất bề dày h, mật độ khối r không đổi. Để thiết lập các phương trình dao động của tấm, ta thừa nhận các giả thiết của lý thuyết tấm mỏng cổ điển như sau : Biến dạng uốn của tấm khi dao động là những biến dạng nhỏ tuân theo định luật Hooke. Trong tấm luôn tồn tại một lớp trung hoà mà khoảng cách giữa các điểm của nó không thay đổi. Khi tấm đồng chất bị uốn ít, lớp trung hoà trùng với mặt cong trung bình chia đôi bề dày của tấm. Ta gọi mặt này là mặt trung hoà. Các phần tử của tấm nằm trên đường thẳng vuông góc với mặt trung hoà, khi tấm bị uốn vẫn nằm trên đường thẳng đó và đường thẳng đó vẫn vuông góc với mặt trung hoà. Không xét đến các ứng suất vuông góc với mặt trung hoà. Bỏ qua ảnh hưởng của biến dạng trượt và quán tính quay. Ta chọn mặt phẳng trùng với mặt trung hòa trạng thái chưa biến dạng làm mặt phẳng toạ độ xy, trục z được chọn vuông gócvới mặt phẳng xy và hướng về phía dưới. Ta ký hiệu u, v, w là các thành phần dịch chuyển của điểm M(x,y,z) của tấm tương ứng theo các trục x, y và z. Ký hiệu u0,v0,w0 là các thành phần dịch chuyển tương ứng của điểm A thuộc mặt trung hoà mà MA vuông góc với mặt trung hoà. Từ các giả thiết trên ta có: u0= v0 = 0 , w0= w(x,y,t) u = -z , v = -z (2.1) Với các lực và mômen đã vẽ như trên, áp dụng nguyên lý d’Alembert ta nhận được các phương trình sau: (2.2) (2.3) (2.4) Thế các biểu thức (2.3) và (2.4) vào phương trình (2.2) ta được : (2.5) Trong các phương trình trên ta sử dụng các ký hiệu: Qx , Qy : Lực cắt trên một đơn vị chiều dài ở các mặt cắt x = const, y = const theo hướng z. x y z Mx , My : Mômen uốn trên một đơn vị chiều dài, vuông góc với mặt cắt x = const, y = const, (Mxy = Myx). p(x,y,t) : Tải trọng ngoài trên một đơn vị diện tích, vuông góc với mặt trung hoà. w(x,y,t) : Dịch chuyển của các điểm thuộc mặt trung hoà theo phương z. Từ các công thức Cauchy quen thuộc trong lý thuyết đàn hồi tuyến tính ta có các quan hệ hình học sau [ ] : , , (2.6) Từ các định luật Hooke đối với trạng thái suất phẳng [ ] ta có : (2.7) Từ đó ta tính được các mômen mặt cắt [ ] : Mx = - D My = - D (2.8) Mxy = - D Trong đó : D = (2.9) Gọi là độ cứng trụ của tấm, còn m gọi là hệ số Poisson. Thế các biểu thức (2.8) vào phương trình (2.5) ta nhận được phương trình dao động uốn của tấm mỏng theo lý thuyết Kirchhoff. D (2.10) Phương trình (2.10) là phương trình dao động uốn của tấm mỏng. Nếu ta sử dụng toán tử (2.11) Thì phương trình (2.10) có dạng : D (2.12) Nếu chỉ xét dao động uốn tự do của tấm, ta lấy p(x,y,t) = 0, từ phương trình (2.12) ta suy ra : D (2.13) 2.2 Các điều kiện biên. Điều kiện biên là những điều kiện trên mặt ngoài của tấm mà ta cần cho trước để nghiệm phương trình (2.10) tương ứng với từng bài toán cụ thể. Trong các điều kiện biên có cả tải trọng q(x,y) tác động ở mặt trên và mặt dưới của tấm. Khi đặt bài toán tổng quát của tấm đã tính đến nó và nó đã có mặt ở một số hạng tự do của phương trình (2.10). Do đó ta chỉ còn điều kiện trên các cạnh của tấm. 2.2.1 Tấm tựa tự do ( gối tựa). Độ võng bằng 0, mômen bằng 0. VD : x = 0 có wùx = 0 = 0, Mxùx = 0 = = 0 Mặt khác, khi wùx = 0 = 0 ị = 0 Vậy điều kiện biên là : wùx = 0 = 0 , = 0 2.2.2 Tấm bị ngàm. Độ võng và góc xoay bằng 0. VD : x = 0 có wùx = 0 = 0 , = 0 2.2.3 Biên tự do. Lực cắt, mômen uốn và mômen xoắn đều bắng 0. VD : x = 0 có Hai điều kiện cuối được sát nhập theo điều kiện Kirchoff. Như vậy ta có điều kiện biên đối với cạnh tự do x = 0 là : 2.3 Dao động uốn của tấm mỏng hình chữ nhật. 2.3.1 Dao động tự do. Ta tìm nghiệm phương trình (2.13) dưới dạng : w(x,y,t) = W(x,y)sin(wt + d) (2.14) Thế (2.14) vào phương trình (2.13) ta có : (2.15) Trong đó ta ký hiệu (2.16) Phương trình (2.15) có thể viết dưới dạng : (2.17) Trong đó toán tử có dạng Nghiệm của phương trình (2.17) có thể tìm dưới dạng W(x,y) = W1(x,y) + W2(x,y) Trong đó W1, W2 tương ứng là nghiệm của các phương trình (2.18a) (2.18b) Nếu xét bài toán dao động tự do của tấm trên nền đàn hồi, thì phương trình dao động tự do của tấm có dạng : Phương trình này có thể đưa về dạng (2.15), nếu ta ký hiệu : Phương trình vi phân đối với các dạng dao động riêng của tấm: (2.19) Trong hệ toạ độ Đề các vuông góc, phương trình (2.19) có dạng: (2.20) Ta có nghiệm phương trình trên dưới dạng W = W1 + W2. Trong đó W1 , W2 là nghiệm của phương trình : (2.21) (2.22) Nghiệm của phương trình (2.21) có dạng : (2.23) Nghiệm của phương trình (2.22) có dạng : W2(x,y) = X(x).Y(y) (2.24) Thế biểu thức (2.24) vào phương trình (2.22) ta có : (2.25) Từ đó suy ra : (2.26) (2.27) (2.28) Nghiệm của (2.26),(2.27) có dạng : (2.29a) (2.29b) Nghiệm của phương trình (2.19) bây giờ có dạng : Trong đó Các giá trị và do đó b và w được xác định từ điều kiện biên. 2.3.2 Dao động cưỡng bức. Dao động uốn cưỡng bức không cản của tấm hình chữ nhật được biểu diễn bởi phương trình vi phân : (2.30) Để tìm nghiệm riêng của phương trình (2.30) ta có nhiều phương pháp khác nhau. Nếu đã biết được các hàm riêng của bài toán dao động tự do Wm,n(x,y), ta có thể tìm một nghiệm riêng của phương trình (2.30) dưới dạng: (2.31) Trong đó qm,n(t) là các hàm cần tìm. Dựa trên tính chất trực giao của các hàm riêng, ta có thể tìm phương trình vi phân đối với hàm qm,n(t) tương tự như trong tính toán dao động uón của dầm. Thế biểu thức (2.31) vào phương trình (2.30) ta được : (2.32) Do Wm,n(x,y) là các hàm riêng nên theo phương trình (2.15) ta có : (2.33) Thế biểu thức (2.33) vào phương trình (2.32) ta suy ra : (2.34) Nhân cả hai vế của phương trình (2.34) với hàm riêng rồi lấy tích phân trên diện tích mặt tấm, chú ý đến tính chất trực giao của hàm riêng . ta nhận được hệ phương trình vi phân thường đối với qm,n(t) (2.35) với (2.36) Nếu không biết trước các hàm riêng Wm,n(x,y) ta áp dụng phương pháp Rit – Galerkin tìm nghiệm dưới dạng : (2.37) trong đó được chọn sao cho thoả mãn điều kiện biên. 2.4 Một Số Bài Toán Ví Dụ 2.4.1 Bài toán 1. Xét dao động cảu tấm mỏng hình chữ nhật chịu điều kiện biên như hình vẽ, chịu tác động của lực phân bố đều q. a x y Giải : Phương trình dao động uốn của tấm mỏng : (1) Với Điều kiện biên bài của toán : x = 0, x = a : y = 0 : y = : Ta tìm nghiệm của (1) dưới dạng : w = w1 + w2 (2) Trong đó : (3) (4) Với w1 : là nghiệm của bài toán chịu mômen phân bố đều trên một đơn vị dài của cạnh ( tương đương với cạnh chịu ngàm ). w2 : là nghiệm của bài toán tấm có 4 cạnh tựa tự do. Khi tấm có 4 cạnh tựa tự do thì độ dốc của mặt võng tại cạnh là : (5) Tại cạnh ngàm chặt nên độ võng do mômen gây ra trên cạnh này sẽ là: (6) Hai biểu thức (5) và (6) bằng nhau nên ta có : (7) Ta có : (8) (9) (10) (11) Mômen uốn Mx được tính như sau : Mx = Mx1 + Mx2 (12) Trong đó : (13) (14) Mômen uốn My được tính như sau : My = My1 + My2 (15) Trong đó : (16) (17) 2.4.2 Bài toán 2. b a Xét dao động của tấm mỏng hình chữ nhật chịu liên kết như hình vẽ, chịu tác dụng của lực phân bố đều trên toàn mặt tấm. A B x y Giải Có thể coi tấm chữ nhật ABCD có hai cạnh kề nhau, x = 0, y = 0 tựa tự do, hai cạnh kia bị ngàm như một phần của tấm bị ngàm tại toàn bộ chu vi x = ± a, y =± b. b a y Tải trọng phân bố đều trên diện tích của tấm cho trước. Lúc này, việc chia ô trên diện tích 2ax2b như hình vẽ sẽ xác định điều kiện tựa tự do tại x = 0 và y = 0. Như vậy, bài toán này được đưa về bài toán tấm ngàm tại toàn bộ chu vi. A B x D C Điều kiện biên của bài toán : x = 0 : y = 0 : x = a : x = b : Ta tìm nghiệm bài toán dưới dạng : w = w1 + w2 (1) Trong đó : w2 – nghiệm của bài toán tấm có 4 cạnh tựa tự do. (2) w1 – nghiệm của bài toán tấm chịu mômen phân bố đều trên một đơn vị dài. (3) ở đây : Góc xoay của mặt phẳng tấm hay độ dốc ở cạnh y = b là : (4) Độ dốc tương ứng với độ võng w1 ở cạnh y = b là : (5) Hai biểu thức (4) và (5) bằng nhau về trị số nhưng ngược dấu nên ta có : (6) Thay (6) vào (3) ta được : Mômen uốn Mx tương ứng với độ võng w1 là : (7) Đạo hàm (3) hai lần theo x và hai lần theo y, sau đó thay vào (7) ta được : (8) (9) (10) Với Em được tính trong công thức (6). Mômen uốn Mx tương ứng với độ võng w2 là : (11) Đạo hàm (2) hai lần theo x và hai lần theo y ta có : (12) (13) Thay (12) và (13) vào (11) ta có : (14) Mômen uốn My tương ứng w1 là : (15) Thay (8) và (9) vào (15) ta có : (16) Mômen uốn My tương ứng w2 là : (17) Thay (12) và (13) vào (17) ta được : (18) Như vậy, mômen uốn Mx sẽ được tính là : Mx = M1x + M2x Với M1x được tính trong (10), M2x được tính trong (14). Mômen uốn My sẽ được tính là : My = M1y + M2y Với M1y được tính trong (16), M2y được tính trong (18). 2.4.3 Bài toán 3. Xét dao động của tấm mỏng hình chữ nhật chịu điều kiện biên như hình vẽ, dưới tác dụng của lực phân bố đều q. x y a b Giải Điều kiện biên của bài toán : x = 0, x = a : (1) y = 0 : (2) y = b : (3) Ta tìm nghiệm của bài toán dưới dạng : w = w1 + w2 Trong đó : w1 – là độ võng của dải tựa tự do có chiều dài a chịu tải trọng phân bố đều và được biểu thị bằng chuỗi : (4) Còn w2 có thể xác định theo chuỗi : (5) Với (6) Những chuỗi (4) và (5) thoả mãn điều kiện biên (1), bốn hằng số trong biểu thức (6) được chọn sao cho điều kiện (2) và (3) được nghiệm đúng. Từ điều kiện (2) ta được : (7) Từ hai điều kiện còn lại, tìm thấy : (8) Trong đó : Sau khi thay các hằng số (7) và (8) vào phương trình (6) và dùng các chuỗi (4) và (5) ta rút ra biểu thức của mặt võng. Đạo hàm hai lần theo x biểu thức (4) ta có : (9) Đạo hàm hai lần theo x biểu thức (5) ta có : (10) Đạo hàm hai lần theo y biểu thức (5) ta có : (11) Mômen uốn Mx tương ứng với độ võng w1 là : (12) Mômen uốn Mx tương ứng với độ võng w2 là : (13) Mômen uốn My tương ứng w1 là : (14) Mômen uốn My tương ứng w2 là : (15) Như vậy, mômen uốn Mx sẽ được tính là : Mx = M1x + M2x Với M1x được tính trong (12), M2x được tính trong (13). Mômen uốn My sẽ được tính là : My = M1y + M2y Với M1y được tính trong (14), M2y được tính trong (15). Chương 3 Phương pháp phần tử hữu hạn 3.1 Giới thiệu chung. Phương pháp phần tử hữu hạn là một phương pháp rất tổng quát và hữu hiệu cho lời giải số nhiều lớp bài toán kỹ thuật khác nhau. Bằng cách thay thế kết cấu thực tế bằng một mô hình dùng để tính toán bao gồm một số hữu hạn phần tử riêng lẻ liên kết với nhau chỉ ở một số hữu hạn điểm nút, tại các điểm nút tồn tại các lực tương tác biểu thị tác dụng qua lại giữa các phần tử kề nhau, nhằm đơn giản hoá tính toán mà vẫn đảm bảo đủ mức chính xác yêu cầu. 3.2 Xấp xỉ bằng phần tử hữu hạn. Giả sử V là miền xác định của một đại lượng cần khảo sát nào đó ( chuyển vị, ứng suất, biến dạng, nhiệt độ…). Ta chia V ra nhiều miền con ve có kích thước và bậc tự do hữu hạn. Đại lượng xấp xỉ của đại lượng trên sẽ được tính trong tập hợp các miền ve. Phương pháp xấp xỉ nhờ các miền con ve được gọi là phương pháp xấp xỉ bằng phương pháp phần tử hữu hạn, nó có một số đặc điểm sau : Xấp xỉ nút trên mỗi miền con ve chỉ liên quan đến những biến nút gắn vào nút của ve và biên của nó. Các hàm xấp xỉ trong mỗi miền con ve được xây dựng sao cho chúng liên tục trên ve và phải thoả mãn các điều kiện liên tục giữa các miền con khác nhau. Các miền con ve được gọi là các phần tử hữu hạn. 3.3 Định nghĩa hình học các phần tử hữu hạn. 3.3.1 Nút hình học. Nút hình học là tập hợp n điểm trên miền V để xác định hình học các phần tử hữu hạn. Chia miền V theo các nút trên, rồi thay miền V bằng một tập hợp các phần tử ve có dạng đơn giản hơn. Mỗi phần tử ve cần chọn sao cho nó được xác định giải tích duy nhất theo các toạ độ nút hình học của phần tử đó, có nghĩa là các toạ độ nằm trong ve hoặc trên biên của nó. 3.3.2 Quy tắc chia miền thành các phần tử. Việc chia miền V thành các phần tử ve phải thoả mãn hai quy tắc sau : Hai phần tử khác nhau chỉ có thể có những điểm chung nằm trên biên của chúng. Điều này loại trừ khả năng giao nhau giữa hai phần tử. Biên giới giữa các phần tử có thể là điểm, đường hay mặt. Biên giới Biên giới Biên giới Hình 3.1 – Biên giới của phần tử hữu hạn Tập hợp tất cả các phần tử ve phải tạo thành một miền càng gần với miền V cho trước càng tốt. Tránh không được tạo lỗ hổng giữa các phần tử. 3.3.3 Các dạng phần tử hữu hạn. Có nhiều dạng phần tử hữu hạn: phần tử một chiều, hai chiều và ba chiều. Trong mỗi dạng đó đại lượng khảo sát có thể biến thiên bậc nhất ( gọi là phần tử bậc nhất), bậc hai hoặc bậc ba…Dưới đây, chúng ta làm quen với một số dạng phần tử thường gặp. a. Phần tử một chiều. PT bậc nhất PT bậc hai PT bậc ba b. Phần tử hai chiều. PT bậc nhất PT bậc hai PT bậc ba Hình 3.3 – Phần tử hai chiều c. Phần tử ba chiều. PT bậc ba PT bậc hai PT bậc nhất Hình 3.4 – Phần tử ba chiều 3.3.4 Lực, chuyển vị, biến dạng và ứng suất. Có thể chia lực tác dạng ra làm ba loại và ta biểu diễn chúng dưới dạng véctơ cột : Lực thể tích f: f = f [fx, fy, fz]T. Lực diện tích T: f = T [Tx, Ty, Tz]T. Lực tập trung Pi: fi = Pi [Tx, Ty, Tz]T. Chuyển vị của một điểm thuộc vật được kí hiệu bởi : u = [u, v, w]T Các thành phần của tenxơ biến dạng được kí hiệu bởi ma trận cột : e = [ex, ey, ez, gyz, gxz, gxy]T Trường hợp biến dạng bé : Các thành phần của tenxơ ứng suất được kí hiệu bởi ma trận: s = [sx, sy, sz, syz, sxz, sxy]T Với vật liệu đàn hồi tuyến tính và đẳng hướng, ta có quan hệ giữa ứng suất với biến dạng : s = De Trong đó : E là môđun đàn hồi, m là hệ số Poisson của vật liệu. 3.3.5 Nguyên lý cực tiểu hoá thế năng toàn phần. Thế năng toàn phần P của một vật thể đàn hồi là tổng của năng lượng biến dạng U và công của ngoại lực tác dụng W : Với vật thể đàn hồi tuyến tính thì năng lượng biến dạng trên một đơn vị thể tích được xác định bởi : sTe Do đó năng lượng biến dạng toàn phần : Công của ngoại lực được xác định bởi : Thế năng toàn phần của vật thể đàn hồi sẽ là : Trong đó u là véctơ chuyển vị và Pi là lực tập trung tại i có chuyển vị là ui. áp dụng nguyên lý cực tiểu hoá thế năng : Đối với một hệ bảo toàn, trong tất cả các di chuyển khả dĩ, di chuyển thực ứng với trạng thái cân bằng sẽ làm cho thế năng đạt cực trị. Khi thế năng đạt giá trị cực tiểu thì vật (hệ) ở trạng thái cân bằng ổn định. 3.4 áp dụng phương pháp phần tử hữu hạn để giải bài toán uốn tấm. 3.4.1 Phần tử không tương thích dạng tam giác. Trong đồ án này chúng ta sẽ xét cụ thể phần tử tam giác và xây dựng chúng theo lý thuyết Kirchoff, việc tính toán tấm mỏng chịu uốn không kể tới biến dạng trượt theo phương pháp phần tử hữu hạn. Có thể thấy rằng trong lý thuyết Kirchoff, khi bỏ qua biến dạng trượt ( gxy =gyz = 0 ), trạng thái ứng suất và biến dạng của tấm hoàn toàn được biểu diễn qua hàm độ võng w(x, y) của mặt trung gian của tấm. Vậy theo phương pháp phần tử hữu hạn, nếu hàm chuyển vị được xấp xỉ hoá và cần tìm là hàm độ võng w(x,y) thì đòi hỏi xấp xỉ của nó không chỉ là liên tục mà còn đảm bảo tính liên tục của các đạo hàm của nó giữa các phần tử. Ngoài ra theo các yêu cầu hội tụ, đa thức xấp xỉ của hàm độ võng w phải có khả năng mô tả các trạng thái biến dạng mà có độ cong không đổi, tức đảm bảo và là hằng số. Và tất nhiên đa thức xấp xỉ cũng đảm bảo tính đẳng hướng hình học. Điều này rõ ràng làm việc chọn xấp xỉ thoả mãn tất cả các yêu cầu trên lại càng phức tạp hơn. Để đơn giản chúng ta chỉ cần đảm bảo sự liên tục của chuyển vị và các góc xoay tại các nút góc theo phương pháp tuyến của các cạnh biên. Những phần tử được xây dựng với hàm xấp xỉ thoả mãn các điều kiện này được gọi là các phần tử tương thích ( compatible hay conformal ). Tuy nhiên việc tìm ma trận độ cứng với các phần tử tương thích là rất phức tạp. Do đó người ta sử dụng một loại phần tử tấm chỉ đòi hỏi sự liên tục của chuyển vị và góc xoay tại các nút góc. Khi đó, các bậc tự do của mỗi nút sẽ là giá trị của và tại nút đó. Phần tử loại này được gọi là phần tử không tương thích ( incompatible hay non- conformal ). Và thực tế cho thấy các phần tử không tương thích vẫn cho kết quả tốt và thường được sử dụng để tính toán tấm chịu uốn. Như vậy có thể thấy rằng các bước phân tích và xây dựng cho các phần tử tấm dựa trên cơ sở lý thuyết Kirchoff vừa dẫn đến việc phải sử dụng các phần tử không tương thích vừa kéo theo các thủ tục phức tạp trong phân tích và lập trình. Trong thời gian gần đây người ta đã xây dựng phần tử tấm chịu uốn dẳng tham số đơn giản và hiệu quả với việc sử dụng các đại lượng cần tìm chỉ đòi hỏi liên tục về giá trị. Các phần tử này được phát triển trên cơ sở lý thuyết của Mindlin và lý thuyết có kể đến biến dạng trượt của tấm. Xét một phần tử tấm mỏng dạng tam giác chịu uốn cùng hệ trục toạ độ địa phương xyz như hình vẽ. z y q7 q9 k q8 q1 q4 i q2 j q5 x q3 q6 Hình 3.5 - Các thành phần chuyển vị y k(0,b) b i(0,0) a j(a,0) x Tại mỗi nút phần tử, độ võng w và các góc xoay qx, qy theo các trục x, y được xem là các bậc tự do của nút. Bỏ qua ảnh hưởng của biến dạng trượt ( theo lý thuyết Kirchoff ) nên các góc xoay chính là đạo hàm của hàm độ võng w theo x, y. Cụ thể các bậc tự do của nút hay còn gọi là các chuyển vị của mỗi nút là : và 3.4.2 Chọn hàm dạng cho phần tử tam giác không tương thích. Theo mô hình tương thích ( tính toán theo chuyển vị ) ta thấy rằng hàm độ võng w(x,y) là đại lượng cần tìm trong bài toán tấm chịu uốn theo lý thuyết Kirchoff. Vì phần tử tam giác như trên có 9 bậc tự do nên hàm độ võng w(x,y) được xấp xỉ hoá bằng một đa thức 9 tham số. Để đảm bảo tính đẳng hướng hình học hàm đa thức xấp xỉ của độ võng như sau : w(x,y) = a1 + a2x + a3y + a4x2 + a5xy + a6 y2 + a7x3 + a8(x2y + xy2) + a9y3 Đến đây ta thử tính tương thích của phần tử tam giác đang xét với hàm độ võng đang xét. Điều này có nghĩa là cần xem rằng chuyển vị và các đạo hàm của nó trên các ( mặt ) biên có xác định một cách duy nhất theo các chuyển vị ( bậc tự do ) không. B A i j i j Hình 3.6 - Xét điều kiện tương thích dọc theo cạnh ij Xét cạnh biên ij là cạnh biên chung giữa hai phần tử A và B kề sát nhau. Trong hệ toạ độ địa phương cạnh ij có phương trình y = 0. Chuyển vị theo phương z hay độ võng dọc theo cạnh này là : wij = (w)y=0 = a1 + a2x + a4x2 + a7x3 (3.1) Độ dốc theo phương x hay dọc theo cạnh này là : a2+ 2a4x + 3 a7x2 (3.2) Vì nút i và nút j là chung của hai phần tử A và B nên các bậc tự do hay các chuyển vị nút này cũng là chung đối với hai phần tử có chung cạnh biên ij này. Tuy nhiên nếu để ý đến bốn chuyển vị nút wi ,, wj ,ta có thể thấy rằng bằng phép đồng nhất bốn bậc tự do này với giá trị hàm w và giá trị đạo hàm tại hai điểm nút i( x = 0 ) và j( x = a ) tức là chúng ta có bốn điều kiện biên sau : (3.3) thì chúng ta có thể hoàn toàn xác định được bốn tham số a1, a2, a4 và a7 một cách duy nhất theo bốn bậc tự do trên. Do đó hoàn toàn được xác định trên cạnh biên ij và tính liên tục của chuyển vị và độ dốc dọc theo cạnh này được đảm bảo khi chuyển từ phần tử A sang phần tử B. Tuy nhiên kết luận trên không đúng đối với độ dốc theo phương y, tức dọc theo cạnh biên chung ij. Thật vậy, trên cạnh biên ij ( y = 0 ) ta có a3 + a5x + a8x2 (3.4) Để xác định được dọc theo cạnh biên ij ta cần xác định được ba tham số a3 ,a5 và a8 . Nhưng chúng ta chỉ có hai phương trình : (3.5) Do đó độ dốc vuông góc với biên ij là không xác định. Độ dốc này tại các nút i và j là như nhau đối với hai phần tử, nhưng có thể là khác nhau tại các điểm khác dọc theo cạnh biên chung ij. Do vậy phần tử này sẽ không tương thích. Sau đay ta sẽ xây dựng phương pháp chuyển vị với phần tử tam giác theo lý thuyết Kirchoff. Theo cách biểu diễn hàm độ võng dưới dạng ma trận như sau: w(x,y) = [P(x,y)] {a} Trong đó : [P(x,y)] = [ 1 x y x2 xy y2 x3 ( x2y + xy2 ) y3 ] { a} = { a1 a2 a3 ... a9 }T Véc tơ chuyển vị nút của phần tử {q}e gồm 9 thành phần là chuyển vị của ba nút i, j, k. Ta có : {q}e= Các tham số [a] đợc xác định bằng cách đồng nhất bậc tự do của nphânf tử với giá trị nhận được tại các nút trong hệ toạ độ địa phương. Với các toạ độ nút i( 0, 0); j( a, 0) và k( 0, b) ta có : Hay ở dạng ma trận : {q}e = [A] {a} (3.6) Trong đó [A] là ma trận : [A] = (3.7) Nghịch đảo của [A] là : [A]-1 = (3.8) Trong đó c = b – a Hàm độ võng được biểu diễn như tổ hợp tuyến tính các hàm dạng : (3.9) [N] = [P(x,y)].[A]-1 = [N1 N2 ... N9] Với các hàm dạng Ni như sau : 3.4.3 Xây dựng ma trận độ cứng phần tử. Các biến dạng trong tấm của phần tử khảo sát. (3.10) Vì w(x,y) = [N] {q}e = [P(x,y)] [A]-1{q}e (3.11) Nên biến dạng được biểu diễn theo {q}e như sau : (3.12) Hay gọn lại : {e}e = [B]{q}e (3.13) Trong đó ma trận biến dạng [B] là : (3.14) ở đây : (3.15) Ma trận độ cứng phần tử tấm dạng tam giác 3 nút không tương thích : (3.16) Bằng cách thay [B] xác định theo (3.14) vào (3.16), ta có : (3.17) Vì [A]-1 chỉ chứa các hằng số nên đưa ra khỏi dấu tích phân, và ta được : (3.18) ở đây S là diện tích phần tử. Hay ở dạng gọn hơn : (3.19) Trong đó : (3.20) ở đây [D]t là ma trận các hệ số đàn hồi của tấm chịu uốn và xác định như sau: [D]t = , với D = là độ cứng trụ của tấm. Với trong (3.15), ta tính sẵn tích các ma trận là (3.21) Trong đó : Để tính tích phân ta chỉ cần lưu ý rằng : là diện tích tam giác phần tử. như mômen tĩnh của tam giác phần tử với trục y. như mômen quán tính của tam giác phần tử với trục y. như mômen quán tính ly tâm của phần tử với hệ trục xy. Thực hiện các tích phân (3.20), ta xác định được ma trận [I] như sau : (3.22) Trong đó : Cuối cùng thực hiện nhân ma trận như (3.19) với [A]-1 đã có ở (3.8), ta xác định được ma trận độ cứng của phần tử tam giác ba nút. [K]e =([A]-1)T[I] [A]-1 Ma trận độ cứng phần tử trong hệ toạ độ tổng thể x’y’z’ được xác định bằng công thức : [K’]e = [T]eT [K]e [T]e (3.23) Trong đó [T]e là ma trận chuyển đổi hệ trục toạ độ. Nếu hệ trục toạ độ tổng thể x’y’z’ có mặt phẳng toạ độ (x’y’) là trùng với mặt phẳng (xy) của hệ toạ độ địa phương, tức là trùng với mặt trung gian tấm thì ma trận [T]e được cho bởi : (3.24) Trong đó : (lx, mx) và (ly, my) là các cosin chỉ phương của các trục x và y trong mặt phẳng (x’,y’) với các trục x’, y’. 3.4.4 Xây dựng ma trận khối lượng phần tử. Với phần tử tấm mỏng chịu uốn dạng tam giác thì hàm độ võng w(x,y) được xấp xỉ bởi đa thức sau : w(x,y) = [P(x,y)] {a} (3.25) với ma trận [P(x,y)] cho bởi : [P(x,y)] = [ 1 x y x2 xy y2 x3 ( x2y + xy2 ) y3 ] {a} = { a1 a2 a3 ... a9 }T và được biểu diễn qua các bậc tự do của phần tử {q}e bởi phương trình : {a} = [A]-1{q}e (3.26) với [A]-1 được cho trong (3.8). Vậy hàm độ võng được biểu diễn theo {q}e là : w(x,y) = ( [P(x,y)] [A]-1) {q}e (3.27) Ngoài ra, các điểm nằm cách mặt trung gian tấm cũng có các chuyển vị theo các phương trục x, y là : (3.28) Vậy cả ba chuyển vị thẳng có thể được biểu diễn theo véctơ {q}e các bậc tự do của phần tử, với chú ý sử dụng (3.27) và (3.28) sẽ là : (3.29) Hay (3.30) ở đây : (3.31) và (3.32) Như vậy, ma trận khối lượng tương thích của phần tử tấm dạng tam giác là (3.33) Ma trận khối lượng tính theo (3.33) trên cơ sở ma trận [N] cho bởi (3.32) là kể tới cả quán tính quay của phần tử do các đoạn thẳng vuông góc mặt trung bình quay xung quanh trục x và y trong quá trình tấm bị uốn. Nếu bỏ qua quán tính quay như trong thực tế vẫn làm thì cũng có nghĩa là bỏ qua các chuyển vị u và v. Khi đó ma trận chính là [P(x,y)]. Và ma trận khối lượng tương thích của phần tử tấm chịu uốn dạng tam giác trong trường hợp này là : (3.34) Trong đó : Sau khi nhân [P(x,y)]T với [P(x,y)] rồi lấy tích phân ta sẽ được ma trận [H]. Trong đó : Sau khi tìm được ma trận [H], thay vào phương trình (3.34) ta có ma trận khối lượng tương thích của phần tử tấm dạng tam giác [M]e . Với phần tử lẻ thì ma trận khối lượng tổng thể sẽ là [M]e’ = [M]e, còn với phần tử chẵn thì ma trận khối lượng tổng thể sẽ được tính như sau: [M]e’ = [T]eT[M]e [T]e 3.4.5 Xây dựng véc tơ tải trọng. a. Trường hợp tấm chịu tải trọng phân bố p(x,y) vuông góc mặt phẳng tấm. Khi đó công dA của di chuyển khả dĩ dw(x,y) bằng : (3.36) Ta có : w(x,y) = [N] {q}e, nên dw(x,y) = [N] {dq}e (3.37) Thay (3.37) vào (3.36) ta được : (3.38) Thế năng biến dạng : dw = (3.39) Ta có : {s} = [D]t{e} ị {s}T = {e}T[D]tT = {e}T[D]t (3.40) {e} = [B] {q}e ị {de}e = {dq}e[B]T (3.41) {e} = {q}eT [B]T (3.42) Trong đó ma trận [B] được tình theo công thức (3.14). Thay các công thức (3.40), (3.41), (3.42) vào công thức (3.39) ta được : dw = (3.43) Từ phương trình cân bằng dw = dA ta có : (3.44) Rút gọn {q}e ở hai vế ta được : (3.45) Trường hợp tấm chịu tải trọng phân bố đều p(x,y) = p0 {P}e = p0 (3.46) Tính tích phân trên ta được : b. Lực tập trung và mômen tập trung. Ta có véctơ tải tổng thể : (3.47) Trong đó : là véctơ tải trọng tập trung đặt tại các nút theo phương tương ứng của các thành phần trong véctơ chuyển vị nút kết cấu z k P i j Hình 3.7 - Tải trọng tập trung Ví dụ : Trường hợp tải trọng tập trung P đặt tại nút j của phần tử ( như hình 3.7). Ta có véctơ tải trọng nút như sau : ={ 0 0 0 P 0 0 0 0 0 }T Thay các phương trình (3.23), (3.34) và (3.47) vào phương trình Lagrange loại 2 ta có được : (3.48) Trong đó : là vectơ gia tốc nút tổng thể . là ma trận cản tổng thể. Sau khi áp đặt các điều kiện động học, chúng ta nhận được phương trình động lực học của hệ có dạng : (3.49) Phương trình (3.49) còn được gọi là phương trình dao động cưỡng bức của kết cấu hay vật thể theo mô hình tương thích của phương pháp phần tử hữu hạn. Nếu xem rằng trong quá trình dao động, hệ không hề mất mát năng lượng do các nguyên nhân như ma sát trong, ma sát của các liên kết , lực cản của môi trường,...tức là nếu bỏ qua sự cản thì phương trình dao động của kết cấu hay vật thể là : (3.50) Phương trình (3.50) được gọi là phương trình dao động cưỡng bức không cản. Giải hệ phương trình tuyến tính (3.50) ta tìm được {q} và từ đó tìm được hàm độ võng theo công thức (3.9). Chương 4 giải thuật và chương trình 4.1 Sơ đồ thuật giải. Chương trình được viết trên ngôn ngữ Matlab 6.5, với ưu điểm: giao diện thiết kế nhanh, gọn, đẹp và dễ sử dụng. Sơ đồ thuật giả của của chương trình như mô tả trong sơ đồ. 4.2 Các bước tính toán. Một số kỹ thuật chính của chương trình : Vẽ cách chia phần tử trên tấm và vẽ chỉ số nút. Đánh chỉ số nút và phần tử. Xây dựng ma trận độ cứng phần tử và ghép nối. Xây dựng ma trận độ khối lượng phần tử và ghép nối. Xây dựng véc tơ tải trọng phần tử và ghép nối. áp đặt điều kiện biên. Giải phương trình chung tìm véc tơ chuyển vị tổng thể. 4.2.1 Vẽ cách chia phần tử trên tấm và chỉ số nút. Kích thước tấm và số đoạn chia theo từng chiều được nhập vào từ bàn phím thông qua FemPlan. 4.2.2 Đánh chỉ số nút và phần tử. 2 Việc đánh số nút và số phần tử tự động, mỗi phần tử với các nút của nó sẽ được lưu trong mảng các phần tử. 3 Chẵn Lẻ 1 Hình 5.1 – Chia phần tử và chiều đánh số nút Nhập a, b, t, E, n, số đoạn chia theo các cạnh, lực P, Mx, My Đánh số phần tử và số nút Xây dựng ma trận độ cứng phần tử [K]eđ [K’]e Xây dựng ma trận khối lượng phần tử [M]eđ [M’]e Véc tơ tải trọng [P]e đ [P’]e Gép nối : Ma trận độ cứng tổng thể [K’] Ma trận khối lượng tổng thể [M’] Véc tơ tải tổng thể [P’] ị Thu gọn ma trận [K’], ma trận [M’], véc tơ tải trọng [P’] Giải hệ phương trình Kết quả Điều kiện biên End Begin Việc đánh số này có ý nghĩa quan trọng khi ghép nối [K’], [M’] và [P’]. Chú ý rằng việc chia phần tử (như hình vẽ) tồn tại hai loại phần tử với cách đánh số khác nhau và chúng có ma trận chuyển toạ độ khác nhau, tạm gọi chúng là phần tử lẻ và phần tử chẵn. for Hang = 1: m_SoHang for Cot = 1: m_SoCot So_Phan_Tu = ((Hang - 1)*m_SoCot + Cot)*2 - 1; Nuti = (Hang - 1)*(m_SoCot+1) + Cot ; Nutj = Hang*(m_SoCot+1) + Cot ; %So phan tu le*/ PTLe(So_Phan_Tu,1) = Nuti ; PTLe(So_Phan_Tu,2) = Nuti + 1; PTLe(So_Phan_Tu,3) = Nutj ; %/*So phan tu chan*/ PTChan(So_Phan_Tu,1) = Nutj+1; PTChan(So_Phan_Tu,2) = Nutj ; PTChan(So_Phan_Tu,3) = Nuti + 1; end end 4.2.3 Xây dựng ma trận độ cứng phần tử và ghép nối. Nhập các thông số của vật liệu gồm modun đàn hồi E, hệ số Poatxong m. Từ đây xây dựng ma trận độ cứng phần tử [K]e với phần tử chẵn và phần tử lẻ. Ma trận độ cứng phần tử lẻ khi chuyển về hệ toạ độ tổng thể không thay đổi ( vì ma trận chuyển Te là ma trận đơn vị ). - Phần tử lẻ: MaTranKeLe( ). Xây dựng theo công thức (3.23) ở đây có một chú ý nhỏ : Do phần tử lẻ phương của các thành phần trong vectơ chuyển vị nút phần tử và tổng thể trùng nhau nên ma trận độ cứng phần tử lẻ trong hai hệ toạ độ địa phương và tổng thể như nhau. - Phần tử chẵn: MaTranKeChan( ). Khi xây dựng được MaTranKeLe( ) thì ta có : MaTranKeChan( ) = MatranTeT( ).MaTranKeLe( ).MatranTe( ) Với cách đánh số nút và phần tử như ở mục 4.2.2 ta có thể ghép các phần tử như sau : -------------------------------------------------------------- ‘Gán giá trị cho biến’ a = Chieudai/sodoan b = Chieurong/sodoan c = b – a D = E*t^3/ (12*(1- n^2) ) -------------------------------------------------------------- ‘Gán giá trị cho ma trận A-1: MaTranA_1(m_a, m_b, m_c)’ MatrixA_1(1,1) = 1.0 ; MatrixA_1(2,2) = 1.0 ; MatrixA_1(3,3) = -1.0 ; MatrixA_1(4,1) = -3/(m_a^2) ; MatrixA_1(4,2) = -2/m_a ; MatrixA_1(4,4) = 3/(m_a^2) ; MatrixA_1(4,5) = -1/m_a ; MatrixA_1(5,2) = m_a/(m_b*m_c) ; MatrixA_1(5,3) = m_b/(m_a*m_c) ; MatrixA_1(5,6) = -m_b/(m_a*m_c) ; MatrixA_1(5,8) = -m_a/(m_b*m_c) ; MatrixA_1(6,1) = -3/(m_b^2) ; MatrixA_1(6,3) = 2/m_b ; MatrixA_1(6,7) = 3/(m_b^2) ; MatrixA_1(6,9) = 1/m_b ; MatrixA_1(7,1) = 2/((m_a^3)) ; MatrixA_1(7,2) = 1/((m_a^2)) ; MatrixA_1(7,4) = -2/((m_a^3)) ; MatrixA_1(7,5) = 1/((m_a^2)) ; MatrixA_1(8,2) = -1/(m_b*m_c) ; MatrixA_1(8,3) = -1/(m_a*m_c) ; MatrixA_1(8,6) = 1/(m_a*m_c) ; MatrixA_1(8,8) = 1/(m_b*m_c) ; MatrixA_1(9,1) = 2/((m_b^3)) ; MatrixA_1(9,3) = -1/((m_b^2)) ; MatrixA_1(9,7) = -2/((m_b^2)) ; MatrixA_1(9,9) = -1/((m_b^2)) ; -------------------------------------------------------------- ‘ Gán giá trị cho ma trận I: MaTranI(m_a, m_b, m_D, m_v)’ MatrixI(4,4) = 2*m_D*m_a*m_b ; MatrixI(5,5) = m_D*(1-m_v)*m_a*m_b ; MatrixI(6,4) = 2*m_D*m_v*m_a*m_b ; MatrixI(6,6)= 2*m_D*m_a*m_b; MatrixI(7,4) = 2*m_D*(m_a^2)*m_b; MatrixI(7,6) = 2*m_D*m_v*(m_a^2)*m_b; MatrixI(7,8) = 3*m_D*(m_a^3)*m_b; MatrixI(8,4) = 2*m_D*(m_a*(m_b^2)+m_v*(m_a^2)*m_b)/3; MatrixI(8,5) = 2*m_D*(1-m_v)*(m_a*(m_b^2) + (m_a^2)*m_b)/3; MatrixI(8,6) = 2*m_D*(m_v*m_a*(m_b^2) + (m_a^2)*m_b)/3; MatrixI(8,7) = m_D*(9*(m_a^2)*(m_b^2)+6*m_v*(m_a^3)*m_b)/6; MatrixI(8,8) = m_D*(((m_a^3)*m_b+(m_b^3)*m_a)*(3- 2*m_v)+(m_a^2)*(m_b^2)*(2-m_v))/3 ; MatrixI(9,4) = 2*m_D*m_v*m_a*(m_b^2) ; MatrixI(9,6) = 2*m_D*m_a*(m_b^2) ; MatrixI(9,7) = 27*m_D*m_v*(m_a^2)*(m_b^2)/6 ; MatrixI(9,8) = m_D*(6*m_v*(m_a^3)*m_b+3*(m_a^2)*(m_b^2))/6; MatrixI(9,9) = 3*m_D*m_a*(m_b^3) ; for i =1:9 for j =1:9 MatrixI(i,j) = MatrixI(i,j) ; MatrixI(j,i) = MatrixI(i,j) ; end end -------------------------------------------------------------- ‘ Gán giá trị cho ma trận chuyển đổi hệ toạ độ: MatranTe( )’ MatrixTe(1,1) = 1 ; MatrixTe(2,2) = -1 ; MatrixTe(3,3) = -1 ; MatrixTe(4,4) = 1 ; MatrixTe(5,5) = -1 ; MatrixTe(6,6) = -1 ; MatrixTe(7,7) = 1 ; MatrixTe(8,8) = -1 ; MatrixTe(9,9) = -1 ; -------------------------------------------------------------- ‘ Xây dựng ma trận độ cứng phần tử [K]e’ [K]e = MatrixA_1t * MatrixI * MatrixA_1 ‘ Tính [K]e trong hệ toạ độ tổng thể bằng cách nhân với ma trận chuyển toạ độ MatrixTe’ Với phần tử lẻ mảng lưu ma trận độ cứng là MaTranKeLe( ) for i = 1: 2 : SoPhanTu ChiSo(1) = 3*PhanTuLe(i,1)-2 ; ChiSo(2) = 3*PhanTuLe(i,1)-1 ; ChiSo(3) = 3*PhanTuLe(i,1) ; ChiSo(4) = 3*PhanTuLe(i,2)-2 ; ChiSo(5) = 3*PhanTuLe(i,2)-1 ; ChiSo(6) = 3*PhanTuLe(i,2) ; ChiSo(7) = 3*PhanTuLe(i,3)-2 ; ChiSo(8) = 3*PhanTuLe(i,3)-1 ; ChiSo(9) = 3*PhanTuLe(i,3) ; for j = 1: 9 for k = 1: 9 KChung(ChiSo(j), ChiSo(k)) = KChung(ChiSo(j), ChiSo(k)) + KLe(j, k) ; end end end Với phần tử lẻ mảng lưu ma trận độ cứng là MaTranKeChan( ) for i = 2: 2: SoPhanTu ChiSo(1) = 3*PhanTuChan(i-1,1)-2 ; ChiSo(2) = 3*PhanTuChan(i-1,1)-1 ; ChiSo(3) = 3*PhanTuChan(i-1,1) ; ChiSo(4) = 3*PhanTuChan(i-1,2)-2 ; ChiSo(5) = 3*PhanTuChan(i-1,2)-1 ; ChiSo(6) = 3*PhanTuChan(i-1,2) ; ChiSo(7) = 3*PhanTuChan(i-1,3)-2 ; ChiSo(8) = 3*PhanTuChan(i-1,3)-1 ; ChiSo(9) = 3*PhanTuChan(i-1,3) ; for j = 1: 9 for k = 1: 9 KChung(ChiSo(j),ChiSo(k)) = KChung(ChiSo(j),ChiSo(k)) + KChan(j, k) ; end end end Ta thấy mảng Chiso(1 to 9) theo cách gán trên chính là chỉ số của phần tử trong ma trận tổng thể. 4.2.4 Xây dựng ma trận khối lượng và ghép nối. Tương tự như xây dựng ma trận độ cứng phần tử [K]e, do có hai loại phần tử là phần tử chẵn và phần tử lẻ nên khi xây dựng ma trận khối lượng phần tử [M]e cũng có hai mảng lưu ma trận khối lượng phần tử đó là: MaTranMeLe() và MaTranMehan( ).Với hai mảng lưu ma trận khối lượng phần tử như trên ta có thể ghép ma trận khối lượng như sau. -------------------------------------------------------------- ‘ Gán biến cho ma trận H: MaTranH(a,b)’ H81 = a^3*b^2/60 + a^2*b^3/60 ; H82 = a^4*b^2/120 + a^3*b^3/180 ; H83 = a^3*b^3/180 + a^2*b^4/120 ; H84 = a^5*b^2/210 + a^4*b^3/420 ; H85 = a^4*b^3/420 + a^3*b^4/420 ; H86 = a^3*b^4/420 + a^2*b^5/210 ; H87 = a^6*b^2/336 + a^5*b^3/840 ; H88 = a^5*b^3/840 + a^4*b^4/560 ; H89 = a^3*b^5/840 + a^2*b^6/336 ; H18 = H81; H28 = H82; H38 = H83 ; H48 = H84; H58 = H85; H68 = H86 ; H78 = H87; H98 = H89 ; for i = 1: 9 for j = 1: 9 H(i,j) = H(i,j) ; H(j,i) = H(i,j) ; end end -------------------------------------------------------------- ‘Tính [M]e trong hệ toạ độ tổng thể bằng cách nhân với ma trận chuyển toạ độ MatrixTe’ Với phần tử lẻ mảng lưu ma trận độ cứng là MaTranMeLe( ) for i = 1: 2: SoPhanTu ChiSo(1) = 3*PhanTuLe(i,1)-2 ; ChiSo(2) = 3*PhanTuLe(i,1)-1 ; ChiSo(3) = 3*PhanTuLe(i,1) ; ChiSo(4) = 3*PhanTuLe(i,2)-2 ; ChiSo(5) = 3*PhanTuLe(i,2)-1 ; ChiSo(6) = 3*PhanTuLe(i,2) ; ChiSo(7) = 3*PhanTuLe(i,3)-2 ; ChiSo(8) = 3*PhanTuLe(i,3)-1 ; ChiSo(9) = 3*PhanTuLe(i,3) ; for j = 1: 9 for k = 1:9 MChung(ChiSo(j), ChiSo(k)) = MChung(ChiSo(j), ChiSo(k)) + MLe(j, k) ; end end end Với phần tử lẻ mảng lưu ma trận độ cứng là MaTranMeChan( ) for i = 2: 2: SoPhanTu ChiSo(1) = 3*PhanTuChan(i-1,1)-2 ; ChiSo(2) = 3*PhanTuChan(i-1,1)-1 ; ChiSo(3) = 3*PhanTuChan(i-1,1) ; ChiSo(4) = 3*PhanTuChan(i-1,2)-2 ; ChiSo(5) = 3*PhanTuChan(i-1,2)-1 ; ChiSo(6) = 3*PhanTuChan(i-1,2) ; ChiSo(7) = 3*PhanTuChan(i-1,3)-2 ; ChiSo(8) = 3*PhanTuChan(i-1,3)-1 ; ChiSo(9) = 3*PhanTuChan(i-1,3) ; for j = 1: 9 for k = 1: 9 MChung(ChiSo(j),ChiSo(k)) = MChung(ChiSo(j),ChiSo(k)) + MChan(j, k) ; end end end 4.2.5 Xây dựng véc tơ tải tổng thể. Nhập tải trọng tác dụng lên tấm bao gồm lực phân bố, lực tập trung, mômen uốn tập trung. Việc xử lý lực nhập như sau : a. Lực tập trung ‘ Gán véc tơ tải phần tử tập trung : LucTapTrung(j,P)’ for i = 1: 9 Pp(i) = 0 ; Pp(j) = P ; end b. Lực phân bố ‘ Gán véc tơ tải phần tử phân bố: LucPhanBo(a,b,Po)’ c = b-a ; if c == 0 'Chon lai a,b' else pp = [] ; pp(1) = a*b/5 ; pp(2) = a^2*b^2/(60*c) + a*b^2/20 - a*b^3/(40*c) ; pp(3) = a^2*b^2/(60*c) - a^2*b/20 - a^3*b/(40*c) ; pp(4) = 3*a*b/20 ; pp(5) = a*b^3/(40*c) - a^2*b^2/(60*c) ; pp(6) = a^2*b/12 - a*b/20 ; pp(7) = 3*a*b/20 ; pp(8) = - a*b^2/30 ; pp(9) = a^3*b/(40*c) - a^2*b^2/(60*c) ; Pphanbo = P0*pp' ; end ‘ Ghép véc tơ tải trọng’ Với phần tử chẵn ta có: for i = 2: 2: Sopt ChiSo(1) = 3*PhanTuChan(i-1,1)-2 ; ChiSo(2) = 3*PhanTuChan(i-1,1)-1 ; ChiSo(3) = 3*PhanTuChan(i-1,1) ; ChiSo(4) = 3*PhanTuChan(i-1,2)-2 ; ChiSo(5) = 3*PhanTuChan(i-1,2)-1 ; ChiSo(6) = 3*PhanTuChan(i-1,2) ; ChiSo(7) = 3*PhanTuChan(i-1,3)-2 ; ChiSo(8) = 3*PhanTuChan(i-1,3)-1 ; ChiSo(9) = 3*PhanTuChan(i-1,3) ; for j = 1:9 P(ChiSo(j)) = P(ChiSo(j)) + Pe(j); end end Với phần tử lẻ ta có: for i = 1: 2: Sopt ChiSo(1) = 3*PhanTuLe(i,1)-2 ; ChiSo(2) = 3*PhanTuLe(i,1)-1 ; ChiSo(3) = 3*PhanTuLe(i,1) ; ChiSo(4) = 3*PhanTuLe(i,2)-2 ; ChiSo(5) = 3*PhanTuLe(i,2)-1 ; ChiSo(6) = 3*PhanTuLe(i,2) ; ChiSo(7) = 3*PhanTuLe(i,3)-2 ; ChiSo(8) = 3*PhanTuLe(i,3)-1 ; ChiSo(9) = 3*PhanTuLe(i,3) ; for j = 1: 9 P(ChiSo(j)) = P(ChiSo(j)) + Pe(j) end end 4.2.6 áp đặt điều kiện biên Kiểu liên kiết được chọn từ Dialog “Nhập liên kết ”. Với các kiểu liên kết sau : Ngam : Hạn chế cả 3 chuyển vị của nút Goi : Hạn chế 1 chuyển vị là (w) Goix : Hạn chế 2 chuyển vị là (w) và qx Goiy : Hạn chế 2 chuyển vị là (w) và qy Thuật toán để áp đặt điều kiện biên. Chuyenvi =[] ; KieuLienKet = [] ; SoNut = (SoHang+1)*(SoCot+1) ; for i = 1:SoNut*3 Chuyenvi(i) = 1; KieuLienKet(i)=0 ; KieuLienKet(1) = 2 ; KieuLienKet(2) = 3 ; end Dem = 3*SoNut ; for i = 1:SoNut if KieuLienKet(i) == 0 %Tu do Chuyenvi(3*i-2) = 1; Chuyenvi(3*i-1) = 1; Chuyenvi(3*i) = 1; end if KieuLienKet(i) == 1 %ngam Chuyenvi(3*i-2) = 0; Chuyenvi(3*i-1) = 0; Chuyenvi(3*i) = 0; Dem = Dem - 3; end ; if KieuLienKet(i) == 2 %Goi do Chuyenvi(3*i-2) = 0; Dem = Dem - 1; end if KieuLienKet(i) == 3 %Goi x Chuyenvi(3*i-2) = 0; Chuyenvi(3*i-1) = 0; Dem = Dem - 2; end if KieuLienKet(i) == 4 %Goi y Chuyenvi(3*i-2) = 0; Chuyenvi(3*i) = 0; Dem = Dem - 2; end end 4.2.7 Thu gọn ma trận khối lượng tổng thể, ma trận độ cứng tổng thể và ma trận tải trọng tổng thể. Kthugon = [] ; Mthugon = [] ; Pthugon = [] ; CVRutgon = [] ; Row = 1; for i= 1:SoNut*3 col = 1; if (Chuyenvi(i) == 1) for j =1: SoNut*3 if (Chuyenvi(j) == 1) Kthugon(Row,col) = K(i,j) ; Mthugon(Row,col) = M(i,j) ; col = col+1; end end Pthugon(Row) = P(i) ; Row = Row + 1 ; end end 4.2.8 Giải hệ phương trình tìm véc tơ chuyển vị nút. Sau khi áp đặt điều kiện biên ta có hệ phương trình thu gọn : Sau khi giải hệ ta tìm được véc tơ chuyển vị nút thu gọn, chương trình sẽ sắp xếp lại véc tơ chuyển vị tổng thể như sau : for i = 1: 3*Sonut if Chuyenvi(i) = 1 Chuyenvi(i) = q(k) ; k = k + 1 ; end end. 4.3 Kiểm nghiệm kết quả tính toán thông qua ví dụ. Trở lại bài toán ví dụ 3 với các thông số : a = 4 (cm) b = 2 (cm) h = 0.5 (cm) E = 2000 (kN/cm2) m = 0.3 q = 100 (kN) Nhập các thông số bài toán qua Dialog “ FemPlan.m ”. Nhập các điều kiện biên tại các nút sau đó nhấn nút Run ta có kết quả : TKieuLienKet = 1 2 0 2 SoDem = 7 SoChuyenVi = 0 0 0 0 1 1 1 1 1 0 1 1 Krutgon = 386.1416 -68.6813 -68.6813 -274.7253 149.5726 -500.6105 -149.5726 -68.6813 295.7112 0 149.5726 0 0 320.5128 -68.6813 0 768.3723 17.1703 377.7473 103.0220 34.3407 -274.7253 149.5726 17.1703 386.1416 -68.6813 183.1502 68.6813 149.5726 0 377.7473 -68.6813 295.7112 -218.2540 -66.7735 -500.6105 0 103.0220 183.1502 -218.2540 729.5482 218.2540 -149.5726 320.5128 34.3407 68.6813 -66.7735 218.2540 570.4365 Mrutgon = 3.4844 1.3067 -2.0533 -2.9867 -1.1511 1.2444 -0.5600 1.3067 0.4978 -0.6067 -1.1511 -0.1867 0.1556 -0.4667 -2.0533 -0.6067 3.9200 3.3600 0.8400 -0.7467 0.7000 -2.9867 -1.1511 3.3600 3.4844 1.3067 -0.7467 0.7156 -1.1511 -0.1867 0.8400 1.3067 0.4978 0.3733 0.2178 1.2444 0.1556 -0.7467 -0.7467 0.3733 -0.0000 -0.8400 -0.5600 -0.4667 0.7000 0.7156 0.2178 -0.8400 0.5600 Tài liệu tham khảo Nguyễn Văn Vượng. Lý thuyết đàn hồi ứng dụng. Nhà xuất bản giáo dục – 1999. Nhữ Phương Mai, Nguyễn Nhật Thăng. Bài tập đàn hồi ứng dụng. Nhà xuất bản giáo dục – 2003. Đặng Việt Cương, Nguyễn Nhật Thăng, Nhữ Phương Mai. Sức bền vật liệu, tập 1 và 2. Nhà xuất bản khoa học và kỹ thuật – 2002. Nguyễn Văn Khang. Dao động kỹ thuật. Nhà xuất bản khoa học và kỹ thuật – 2003. Timoshenko. Những vấn đề dao động trong kỹ thuật. Nhà xuất bản khoa học và kỹ thuật – 2000 Timoshenko, Woinowsky Krieger. Theory of plates and shells. Mc Graw – Hill, New York, 1959. Bùi Trọng Lựu, Nguyễn Văn Vượng. Bài tập sức bền vật liệu. Nhà xuất bản giáo dục – 1996. Đặng Việt Cương, Lê Quang Minh. Lý thuyết dẻo ứng dụng. Nhà xuất bản khoa học và kỹ thuật - 2003. Chu Quốc Thắng. Phương pháp phần tử hữu hạn. Nhà xuất bản khoa học và kỹ thuật – 1997. Trần ích Thịnh, Trần Đức Trung, Nguyễn Việt Hùng. Phương pháp phần tử hữu hạn trong kỹ thuật. Nhà xuất bản khoa học kỹ thuật - 2000. Phạm Quang Hân, Phạm Quang Huy, Hồ Xuân Phương. Giáo trình thực hành tính toán kết cấu công trình trong xây dựng. Nhà xuất bản thống kê – 2001. Nguyễn Phùng Quang. Matlab và Simulink. Nhà xuất bản khoa học và kỹ thuật – 2004. Phạm Văn Điển. Tính toán, lập trình và giảng dạy toán học trên Maple. Nhà xuất bản khoa học và kỹ thuật – 2002. Mục lục Tài liệu tham khảo Nguyễn Văn Vượng. Lý thuyết đàn hồi ứng dụng. Nhà xuất bản giáo dục – 1999. Nhữ Phương Mai, Nguyễn Nhật Thăng. Bài tập đàn hồi ứng dụng. Nhà xuất bản giáo dục – 2003. Đặng Việt Cương, Nguyễn Nhật Thăng, Nhữ Phương Mai. Sức bền vật liệu, tập 1 và 2. Nhà xuất bản khoa học và kỹ thuật – 2002. Nguyễn Văn Khang. Dao động kỹ thuật. Nhà xuất bản khoa học và kỹ thuật – 2003. Timoshenko. Những vấn đề dao động trong kỹ thuật. Nhà xuất bản khoa học và kỹ thuật – 2000 Timoshenko, Woinowsky Krieger. Theory of plates and shells. Mc Graw – Hill, New York, 1959. Bùi Trọng Lựu, Nguyễn Văn Vượng. Bài tập sức bền vật liệu. Nhà xuất bản giáo dục – 1996. Đặng Việt Cương, Lê Quang Minh. Lý thuyết dẻo ứng dụng. Nhà xuất bản khoa học và kỹ thuật - 2003. Chu Quốc Thắng. Phương pháp phần tử hữu hạn. Nhà xuất bản khoa học và kỹ thuật – 1997. Trần ích Thịnh, Trần Đức Trung, Nguyễn Việt Hùng. Phương pháp phần tử hữu hạn trong kỹ thuật. Nhà xuất bản khoa học kỹ thuật - 2000.