Luận văn Phương pháp phần tử hữu hạn cho tấm composite lớp trên nền đàn hồi

đó: (1.1) Còn là chuyển vị của phương ngang, phương dọc và độ võng của các điểm giữa thuộc mặt phẳng của tấm; là các biến dạng tại mặt giữa; là các biến thiên độ cong của tấm. Chúng thỏa mãn phương trình tương thích biến dạng: (1.2) 1.1.2. Quan hệ ứng suất biến dạng của tấm composite lớp Sử dụng giả thiết Kirchhoff có thể bỏ qua thành phần ứng suất vuông góc với mặt giữa: Liên hệ ứng suất biến dạng của lớp composite thứ k của tấm là [5]: (1.3) Trong đó ký hiệu các thành phần biến dạng trong mặt phẳng lớp thứ k: , , Trường hợp phương của sợi lệch một góc với trục của tấm, thay ma trận bằng ma trận . Trong đó tính qua theo công thức [4]: (1.4) Biểu thức các hằng số độ cứng qua các mô đun đàn hồi trong hệ trục chính như sau: (1.5) , , trong đó: ,là các môđun đàn hồi của tấm theo phương trục chính của lớp vật liệu composite; là hệ số Poisson của vật liệu, là môđun trượt trong hệ trục chính của lớp vật liệu. Lực pháp, lực tiếp, mômen uốn, mômen xoắn được xác định theo công thức: (1.6) Ở đây: Thay (1.3) vào (1.6) ta được: Với: (i, j = 1, 2, 6) (1.7) Mối quan hệ giữa lực pháp, lực tiếp, momen uốn, momen xoắn và biến dạng, biến thiên độ cong được viết dưới dạng ma trận như sau: (1.8) với: N = 4 là số lớp của tấm; lần lượt là ma trận độ cứng dãn nén, độ cứng tương tác dãn – uốn – xoắn và độ cứng uốn của tấm composite lớp. Giả thiết tấm xếp lớp đối xứng qua mặt giữa ta có và xem các đại lượng là nhỏ có thể bỏ qua. Khai triển (1.8) ta được biểu thức lực màng của tấm composite lớp: (1.9) Giải ngược lại suy ra: , , Trong đó: , , Và mômen trong của tấm composite lớp được tính theo công thức: (1.10) Trong đó: , , , , 1.1.3. Phương trình chuyển động của tấm composite lớp trên nền đàn hồi Phương trình chuyển động của tấm composite lớp trên nền đàn hồi với hệ số nền K theo mô hình Vinkler được viết như sau [5]: (1.11) Trong đó được xác định theo công thức: là mật độ khối lượng của lớp thứ k, là lực phân bố 1.2. Nghiệm của bài toán Giả thiết lực ngang phân bố đều và mật độ khối lượng của lớp thứ k là hằng số. Khi đó ta có: Theo Volmir [8] thì các số hạng quán tính trong hai phương trình đầu của (1.11) có thể bỏ qua. Do vậy phương trình chuyển động có dạng: (1.12) (1.13) (1.14) Phương trình (1.12), (1.13) thỏa mãn khi đưa vào hàm ứng suất φ dạng: (1.15) Thay (1.9) và (1.15) vào phương trình tương thích biến dạng ta được: Suy ra: Vì nên Do đó phương trình tương thích biến dạng trở thành: (1.16) Thay (1.10) và (1.15) vào (1.14) ta được: Phương trình (1.14) trở thành: hay Do đó (1.14) được viết dưới dạng (1.17) Cuối cùng bài toán bao gồm hai phương trình đạo hàm riêng phi tuyến (1.16), (1.17). Điều kiện biên tựa bản lề thỏa mãn nếu chọn hàm độ võng dạng: (1.18) trong đó a và b là độ dài các cạnh của tấm, là độ võng cực đại. Thế (1.18) vào vế phải của phương trình (1.16) ta có Suy ra Vậy (1.16) trở thành: Khi đó nghiệm của phương trình sẽ là: (1.19) Thế (1.19) vào (1.16) ta được: Suy ra Do đó (1.20) Thế (1.18), (1.19) vào (1.17) ta được: (1.21) Đây là phương trình cân bằng của tấm composite lớp trên nền đàn hồi. 1.2.1. Bài toán tĩnh Ta có , lực phân bố đều = const Khi đó phương trình cân bằng của tấm trên nền đàn hồi có dạng: Áp dụng phương pháp Bubnov – Galerkin: ta được Sau một số phép biến đổi và sử dụng (1.20) khi đó phương trình (1.22) trở thành: (1.23) trong đó , Đây là phương trình cân bằng xác định độ võng cực đại theo phương pháp lý thuyết bản vỏ mỏng trong trường hợp không có dao động. 1.2.2. Bài toán động Ta có , lực phân bố đều Khi đó phương trình cân bằng của tấm trên nền đàn hồi có dạng: Áp dụng phương pháp Bubnov – Galerkin: ta được (1.24) Sau một số phép biến đổi và sử dụng (1.20) khi đó phương trình (1.24) trở thành: (1.25) trong đó: , , Đây là phương trình dao động của tấm composite trên nền đàn hồi xác định độ võng cực đại theo phương pháp lý thuyết bản vỏ mỏng. CHƯƠNG 2 PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN CHO TẤM COMPOSITE LỚP TRÊN NỀN ĐÀN HỒI 2.1. Khái niệm về phương pháp phần tử hữu hạn Phương pháp phần tử hữu hạn (PP PTHH) là một phương pháp số đặc biệt có hiệu quả để tìm dạng gần đúng của một hàm chưa biết trong miền xác định V của nó. Tuy nhiên PP PTHH không tìm dạng xấp xỉ của hàm cần tìm trên toàn miền V mà chỉ trong từng miền con Ve (phần tử) thuộc miền xác định V. Do đó phương pháp này rất thích hợp với hàng loạt bài toán vật lý và kỹ thuật trong đó hàm cần tìm được xác định trên những miền phức tạp gồm nhiều vùng nhỏ có đặc tính hình học, vật lý khác nhau chịu những điều kiện biên khác nhau. Trong PP PTHH miền V được chia thành một số hữu hạn các miền con, gọi là phần tử. Các phần tử này được nối kết với nhau tại các điểm định trước trên biên phần tử, gọi là nút. Trong phạm vi mỗi phần tử đại lượng cần tìm được lấy xấp xỉ trong một hàm đơn giản gọi là các hàm xấp xỉ. Và các hàm xấp xỉ này được biểu diễn qua các giá trị của hàm (và có khi cả các giá trị của đạo hàm của nó) tại các điểm nút trên phần tử. Các giá trị này được gọi là các bậc tự do của phần tử và được xem là ẩn số cần tìm của bài toán. 2.2. Tính toán tấm composite lớp trên nền đàn hồi theo phương pháp phần tử hữu hạn Bài toán: Xét tấm composite chữ nhật N lớp trên nền đàn hồi có kích thước (a x b), độ dày h, chịu lực phân bố đều q(x,y) = q = const. Bước 1: Rời rạc hóa miền vật thể. Chia tấm thành Ne phần tử hình chữ nhật Ve (e = 1Ne) có kích thước bằng nhau. Chọn hệ tọa độ Đề-các có x, y là các trục tọa độ nằm trong mặt phẳng giữa theo các cạnh, trục z hướng theo phương pháp tuyến với mặt giữa. Giả thiết tấm có các lớp được xếp đối xứng qua mặt trung bình. Khi đó xét điểm M(x,y) nằm cách mặt trung hòa khoảng z. Khi đó chuyển vị tại M theo phương Ox, Oy tương ứng là: , Tại mỗi điểm M(x,y) trong phần tử e giả sử có các thành phần chuyển vị là {w(x, y), θx(x, y), θy(x, y)}. Mỗi phần tử hình chữ nhật có bốn nút, mỗi nút có ba bậc tự do. Do đó phần tử có mười hai bậc tự do nên việc chọn hàm xấp xỉ phải đảm bảo có mười hai tham số. Bước 2: Chọn hàm xấp xỉ Chọn hàm độ võng dưới dạng đa thức gồm mười hai số hạng [3]: Góc xoay của tấm tại M trên mặt giữa theo các phương Ox, Oy tương ứng: Chuyển dịch của phần tử Thay tọa độ các nút i, j, k, l đã biết vào phương trình (2.1) ta được véctơ chuyển vị nút : với Suy ra : (2.2) Từ (2.1) và (2.2) có : trong đó: [A] là ma trận tọa độ nút phần tử, [N(x,y)] là ma trận hàm dạng Bước 3 : Biểu diễn biến dạng và ứng suất qua chuyển vị nút Trạng thái biến dạng phần tử tính theo hệ thức Cauchy: (2.3) Liên hệ ma trận biến dạng và ma trận độ cong như sau: với (2.4) Liên hệ ứng suất - biến dạng được biểu diễn theo định luật Hooke: Công thức xác định mômen uốn, mômen xoắn là: , , (2.5) Suy ra (2.6) (do ) trong đó là ma trận độ cứng uốn của tấm composite lớp: Các phần tử Dij được tính như sau: Công thức tính tại từng lớp được tính như sau: Lớp 1 và lớp 4 đối xứng nhau (góc xoay ) nên . Cụ thể: Lớp 2 và lớp 3 đối xứng nhau (góc xoay) nên . Cụ thể: trong đó được xác định theo công thức (1.5). Ký hiệu :. Từ (2.6) ta được: (2.7) Bước 4: Thiết lập thế năng toàn phần, ma trận độ cứng, véctơ tải của phần tử. Thế năng biến dạng của phần tử tấm: Thay (2.3) vào hệ thức trên ta được : Theo (2.5) ta được : Từ (2.4) và (2.7) suy ra trong đó được gọi là ma trận độ cứng phần tử Thế năng của nền đàn hồi với hệ số nền K là : với Vậy thế năng của tấm composite lớp trên nền đàn hồi sẽ là: với được gọi là ma trận độ cứng phần tử tấm composite lớp trên nền đàn hồi. Công của ngoại lực do q gây ra trên độ võng của phần tử e là: trong đó là véctơ tải phần tử của tấm composite trên nền đàn hồi. Bước 5: Ghép nối phần tử, tìm ma trận độ cứng tổng thể và véctơ tải tổng thể. Giả sử là véctơ chuyển vị nút tổng thể , còn (e = 1...Ne) là chuyển vị nút của một phần tử tương ứng. Vậy giữa hai đại lượng này phải có liên hệ: , trong đó gọi là ma trận định vị của phần tử có các phần tử cho bởi quy tắc sau: Thế năng toàn phần của tấm composite lớp trên nền đàn hồi là: Trong đó : tương ứng là ma trận độ cứng tổng thể, véctơ tải tổng thể của tấm composite lớp trên nền đàn hồi. Bước 6: Áp điều kiện biên tìm chuyển vị nút chưa biết Điều kiện cực trị thế năng toàn phần là Hệ trên trở thành : (2.8) Đây chính là hệ phương trình xác định chuyển vị nút. Áp điều kiện biên : Tấm khớp bản lề tại bốn cạnh: Trên biên có: W = 0 Trên cạnh y = 0, y = b thì θx = 0 (2.9) Trên cạnh x = 0, x = a thì θy = 0 Tấm có hai cạnh ngàm, hai cạnh khớp bản lề Trên biên có W = 0 Tại hai cạnh ngàm: θx = θy = 0 (2.10) Tại hai cạnh khớp bản lề: θx = 0 hoặc θy = 0 Khi đó hệ phương trình xác định chuyển vị nút được viết dưới dạng: (2.11) trong đó thu được bằng cách bỏ đi số hàng và số cột của ma trận độ cứng tổng thể mà ở các vị trí đó các thành phần của véctơ bằng không,cũng thu được từ bằng cách bỏ đi các thành phần tương ứng. Nghiệm chuyển vị nút cần tìm được tính theo công thức: (2.12) CHƯƠNG 3 TÍNH TOÁN SỐ CHO TẤM COMPOSITE LỚP TRÊN NỀN ĐÀN HỒI Xét tấm composite lớp với các thông số hình học và cấu tạo như sau : - Kích thước tấm: a = b = 0,4m. Chiều dày tấm h = 0,01m - Kiểu sắp xếp lớp: [45/-45/-45/45]. - Vật liệu có các hằng số đàn hồi của mỗi lớp : E1 = 127,4 GPa; E2 = 13,0 GPa; ν = 0,38; G12 = 6,4 GPa, kg/m3 Lực phân bố đều : q0 = 1000 (N/m2), Hệ số nền K = 108 (N/m3) Điều kiện biên: Bản lề bốn cạnh 3.1. Tính toán số bài toán tĩnh của tấm composite lớp trên nền đàn hồi Chia bản composite thành 16 miền chữ nhật như hình vẽ với 25 nút. Các phần tử là các hình chữ nhật giống nhau kích thước . y 25 24 23 22 21 20 19 18 17 16 15 14 12 13 11 6 10 9 8 7 5 4 3 2 1 Hình 2: Bản composite được chia thành 16 phần tử với 25 nút x Phương trình (1.23) có nghiệm độ võng cực đại: f = 0,99991.10-5 m Kết quả tính toán theo phương pháp PTHH cho chuyển vị tại các nút chưa biết: (2) = Q6 = 0,78307.10-4 (3 ) = Q9 = 0,75017.10-4 (4) = Q12 = 0,39923.10-4 (6) = Q17 = 0,78307.10-4 w(7) = Q19 = 0,52654.10-5m (7) = Q20 = 0,30816.10-4 (7) = Q21 = 0,30816.10-4 w(8) = Q22 = 0,64841.10-5m (8) = Q23 = -0,43885.10-5 (8) = Q24 = 0,43322.10-4 w(9) = Q25 = 0,44177.10-5m (9) = Q26 = -0.37644.10-4 (9) = Q27 = 0.37644.10-4 (10) = Q29 = -0.39923.10-4 (10) = Q32 = 0.75017.10-4 w(12) = Q34 = 0.64841.10-5m (12) = Q35 = 0.43322.10-4 (12) = Q36 = -0.43885.10-5 w(13) = Q37 = 0.86056.10-5m (13) = Q38 = 0,5014.10-34 (13) = Q39 = -0,95200.10-35 w(14) = Q40 = 0,64841.10-5m (14) = Q41 = -0,43322.10-4 (14) = Q42 = 0,43885.10-5 (15) = Q44 = -0,75017.10-4 (16) = Q47 = 0,39923.10-4 w(17) = Q49 = 0,44177.10-5m (17) = Q50 = 0,37644.10-4 (17) = Q51 = -0,37644.10-4 w(18) = Q52 = 0,64841.10-5m (18) = Q53 = 0.43885.10-5 (18) = Q54 = -0.43322.10-4 w(19) = Q55 = 0 .52654.10-4m (19) = Q56 = -0.30816.10-4 (19) = Q57 = -0.30816.10-4 (20) = Q59 = -0.78307.10-4 (22) = Q66 = -0.39923.10-4 (23) = Q69 = -0.75017.10-4 (24) = Q72 = -0.78307.10-4 Sai số hai phương pháp giải tích và phần tử hữu hạn là: Ta xét thêm một số trường hợp sau: TH1: a = b = 0,4 m, ta có , TH2: a = b = 0,4 m, ta có TH3: a = 0,4 m, b = 0,35m, N/m3 ta có TH4: a = 0,4 m, b = 0,35m, N/m3 ta có TH5: a = 0,4 m, b = 0,35m, N/m3 ta có Đồ thị Khi tấm có kích thước bằng nhau a = b = 0,4m chịu lực phân bố đều q0 = 1000N/m2 và hệ số nền K thay đổi bằng (0,5.108; 108; 2.108)N/m3 thì độ võng của tấm tính theo phương pháp lý thuyết vỏ mỏng và phương pháp phần tử hữu hạn được thể hiện trong các hình (3, 4, 5) sau: Hình 3. Đồ thị biểu diễn độ võng của tấm vuông trên nền có hệ số K=0,5.108N/m3 theo PP lý thuyết bản vỏ mỏng và PP phần tử hữu hạn Hình 4. Đồ thị biểu diễn độ võng của tấm vuông trên nền có hệ số K = 108N/m3 theo hai phương pháp lý thuyết bản vỏ mỏng và phần tử hữu hạn Hình 5. Đồ thị biểu diễn độ võng của tấm vuông trên nền có hệ số K = 2.108N/m3 theo hai phương pháp lý thuyết bản vỏ mỏng và phần tử hữu hạn Khi tấm có kích thước a = 0,4m; b = 0,35m chịu lực phân bố đều q0 = 1000N/m2 và hệ số nền K thay đổi bằng (0,5.108; 108; 2.108)N/m3 thì độ võng của tấm tính theo phương pháp lý thuyết vỏ mỏng và phương pháp phần tử hữu hạn được thể hiện trong các hình (6, 7, 8) sau: Hình 6. Đồ thị biểu diễn độ võng của tấm chữ nhật với hệ số nền K = 0,5.108N/m3 theo hai phương pháp lý thuyết bản vỏ mỏng và phần tử hữu hạn Hình 7. Đồ thị biểu diễn độ võng của tấm chữ nhật với hệ số nền K = 108N/m3 theo hai phương pháp lý thuyết bản vỏ mỏng và phần tử hữu hạn Hình 8. Đồ thị biểu diễn độ võng của tấm chữ nhật trên nền có hệ số K = 2.108N/m3 theo hai phương pháp lý thuyết bản vỏ mỏng và phần tử hữu hạn Nhận xét: - Trong điều kiện tải trọng như nhau tác động lên các tấm composite thì tấm đặt trên nền có hệ số lớn hơn sẽ có độ võng nhỏ hơn. - Sai số hai phương pháp sẽ nhỏ hơn khi áp dụng cho tấm composite xếp lớp có kích thước hai cạnh gần nhau hơn 3.2. Tính toán số bài toán động của tấm composite lớp trên nền đàn hồi 3.2.1. Phương pháp Newmark Phương trình dao động phi tuyến có dạng: hay Đặt và Phương trình trên trở thành: (3.1) Giải phương trình (3.1) bằng phương pháp bước lặp với sơ đồ Newmark. Chia quá trình tính toán thành n giai đoạn, phương trình vi phân (3.1) chuyển về dạng phương trình đại số trong từng giai đoạn để xác định khi đã biết ,, theo công thức sau: (3.2) trong đó: (3.3) Thay (3.3) theo giai đoạn n đã biết được phương trình (3.2). Do hệ số là phi tuyến nên ta sử dụng phương pháp bước lặp để giải phương trình: ( k là số bước lặp) Sau khi giải (3.2) tìm được . Ta cần tính tiếp theo công thức sau: (3.4) Sau khi giải (3.2) và (3.4) ta tính được: , Lại lặp lại quá trình trên tìm được , với điều kiện đầu: , , Tại mỗi bước thứ k+1 thì hệ số được tính từ bước thứ k trước đó. Quá trình tính toán dừng lại khi thỏa mãn điều kiện hội tụ : 3.2.2. Kết quả giải số Trường hợp tấm có kích thước a = b = 0.4m Tần số dao động tuyến tính tự do trong phương trình (1.23) là rad/s tương ứng hệ số nền K thay đổi bằng (0,5.108; 108; 2.108)N/m3 ta chọn tần số lực ngoài ω0 = 2000rad/s, với biên độ lực ngoài q0 = 1000N/m2, thì đáp ứng phi tuyến được thể hiện trong các hình (9, 10, 11) sau: Hình 9. Đồ thị biểu diễn đáp ứng phi tuyến của tấm vuông trên nền có hệ số K = 0,5.108N/m3 Hình 10. Đồ thị biểu diễn đáp ứng phi tuyến của tấm vuông trên nền có hệ số K = 108N/m3 Hình 11. Đồ thị biểu diễn đáp ứng phi tuyến của tấm vuông với hệ số nền K = 2.108N/m3 Chọn biên độ lực ngoài q0 = 1000N/m2. Với hệ số nền K = 108N/m3, tần số dao động tuyến tính tự do rad/s, khi tần số lực ngoài thay đổi thì đáp ứng phi tuyến của tấm được thể hiện như hình (12, 13, 14). Còn với các hệ số nền K = (0,5.108; 2.108)N/m3, tần số rad/s, ta chọn tần số lực ngoài rad/s thì đáp ứng phi tuyến của tấm được thể hiện như hình (15, 16). Hình 12. Đồ thị biểu diễn đáp ứng phi tuyến của tấm vuông với tần số lực ngoài ω0 = 2000rad/s Hình 13. Đồ thị biểu diễn đáp ứng phi tuyến của tấm vuông với tần số lực ngoài ω0 = 2800rad/s Hình 14. Đồ thị biểu diễn đáp ứng phi tuyến của tấm vuông trên nền có hệ số K = 108 N/m3 với tần số lực ngoài ω0 xấp xỉ tần số dao động tuyến tính tự do ω1. Hình 15. Đồ thị biểu diễn đáp ứng phi tuyến của tấm vuông trên nền có hệ số K = 0,5.108N/m3 khi tần số lực ngoài ω0 xấp xỉ tần số dao động tuyến tính tự do ω1. Hình 16. Đồ thị biểu diễn đáp ứng phi tuyến của tấm vuông trên nền có hệ số K = 2.108 N/m3 khi tần số lực ngoài ω0 xấp xỉ tần số dao động tuyến tính tự do ω1. Với hệ số nền K = 108N/m3, tần số dao động tuyến tính tự do rad/s, ta chọn tần số lực ngoài ω0 = 2000rad/s, khi biên độ ngoại lực thay đổi thì đáp ứng phi tuyến thay đổi như hình (17, 18, 19) như sau: Hình 17. Đồ thị biểu diễn đáp ứng phi tuyến của tấm vuông khi cho biên độ ngoại lực thay đổi q0 = 500N/m2 Hình 18. Đồ thị biểu diễn đáp ứng phi tuyến của tấm vuông khi cho biên độ ngoại lực thay đổi q0 = 1000N/m2 Hình 19. Đồ thị biểu diễn đáp ứng phi tuyến của tấm vuông khi cho biên độ ngoại lực thay đổi q0 = 2500N/m2 Trường hợp tấm có kích thước a = 0.4m, b = 0.2m Với biên độ lực ngoài q0 = 1000N/m2, hệ số nền K = 108 N/m3, tần số dao động tuyến tính tự do rad/s ta chọn tần số lực ngoài ω0 = 2000rad/s và ω0 = 3130rad/s thì đáp ứng phi tuyến của tấm được thể hiện như hình (20, 21). Hình 20. Đồ thị biểu diễn đáp ứng phi tuyến của tấm chữ nhật khi cho tần số lực ngoài thay đổi ω0 = 2000rad/s Hình 21. Đồ thị biểu diễn đáp ứng phi tuyến của tấm vuông khi tần số lực ngoài xấp xỉ tần số dao động tuyến tính tự do Nhận xét Với tấm có cùng kích thước: - Trong điều kiện biên độ tải trọng và tần số lực ngoài như nhau thì tấm nào đặt trên nền có hệ số lớn hơn sẽ có biên độ đáp ứng phi tuyến nhỏ hơn. - Trong điều kiện các tấm chịu tải trọng như nhau và được đặt trên các nền có cùng hệ số, khi tần số dao dộng lực ngoài ω0 và tần số dao động tuyến tính tự do ω1 càng gần nhau thì biên độ đáp ứng phi tuyến của tấm càng cao. - Trong điều kiện biên độ tải trọng khác nhau thì đáp ứng phi tuyến của tấm là khác nhau: biên độ tải trọng càng tăng thì biên độ đáp ứng phi tuyến của tấm càng cao. Với các tấm có cùng hệ số nền và lực ngoài như nhau nhưng kích thước khác nhau thì tấm có kích thước nhỏ hơn sẽ cho biên độ đáp ứng phi tuyến cũng nhỏ hơn và ngược lại (Hình 14, 21). KẾT LUẬN Luận văn đã phân tích và đưa ra các hệ thức cho bài toán tĩnh và động lực của tấm composite lớp trên nền đàn hồi dựa trên lý thuyết vỏ mỏng có tính đến phi tuyến hình học. Dùng phương pháp hàm ứng suất và Bubnov-Galerkin nhận được phương trình dao động phi tuyến của tấm. Đã khảo sát ảnh hưởng của nền, tần số dao động lực ngoài, biên độ lực, kích thước tấm tới độ võng và đáp ứng phi tuyến tấm composite lớp. Lời giải số bài toán dao động tìm được theo phương pháp bước lặp và sử dụng sơ đồ tính toán Newmark. Luận văn đã sử dụng phương pháp phần tử hữu hạn với mô hình phần tử tứ giác có 12 bậc tự do để giải bài toán tĩnh của tấm composite trên nền đàn hồi. Đã so sánh kết quả độ võng theo các phương pháp giải tích và phần tử hữu hạn. TÀI LIỆU THAM KHẢO Vũ Đỗ Long (2006). Tính toán phi tuyến vỏ thoải composite lớp. Tuyển tập Hội nghị khoa học toàn quốc Cơ học Vật Rắn Biến Dạng, lần thứ 8, NXB Khoa Học Công Nghệ Việt Nam, 463-472 Vũ Đỗ Long, Phạm Văn Khoa (2009). Phân tích phi tuyến tĩnh và động tấm composite lớp có gân gia cường. Tuyển tập Hội nghị Cơ học toàn quốc lần thứ 8, Hà nội. Chu Quốc Thắng (1997). Phương pháp Phần tử hữu hạn, NXB Khoa học và kỹ thuật. Trần Ích Thịnh (1994). Vật liệu composite. Cơ học và tính toán kết cấu, NXB Giáo dục. Dao Huy Bich (2006). Nonlinear analysis of laminated composite doubly curved shallow shells. Vietnam Journal of Mechanics, VAST, Vol 28, No 1, 43-55. Dao Huy Bich, Tran Thanh Tuan (2005). Dynamical analysis of multilayered reinforced composite cylindrical shells. Vietnam Journal of Mechanics, VAST, Vol 27, No 4, 220-228. Dao Huy Bich, Vu Do Long (2007). Nonlinear dynamic analysis of laminated reinforced composite doubly curved shallow shells. Vietnam Journal of Mechanics, VAST, Vol 29, No 3, 257-269. Volmir A.S (1972). Non-linear dynamics of plates and shells. Science Edition, Moscow. PHỤ LỤC Chương trình dùng Matlap tính toán trong trường hợp bài toán tĩnh function uonbancomposite1 syms x y z f a=0.4;b=0.4;h=0.01; q0=1000; k=10^8; E1=12.74*10^10; E2=13*10^9; v12=0.38; G12=6.4*10^9; Q11=E1/(1-E2*(v12^2)/E1); Q12=E2*v12/(1-E2*(v12^2)/E1); Q22=E2/(1-E2*(v12^2)/E1); Q66=G12; % Lop 1 và 4 doi xung nhau t1=pi/4; Q11n=Q11*cos(t1)^4+ Q22*sin(t1)^4+2*(Q12+2*Q66)*(sin(t1)^2)*(cos(t1)^2); Q12n=(Q11+Q22-4*Q66)*(sin(t1)^2)*(cos(t1)^2)+Q12*(sin(t1)^4+cos(t1)^4); Q16n=(Q11-Q12-2*Q66)*sin(t1)*(cos(t1)^3)+(Q12-Q22+2*Q66)*(sin(t1)^3)*cos(t1); Q22n=Q11*sin(t1)^4+ Q22*cos(t1)^4+2*(Q12+2*Q66)*(sin(t1)^2)*(cos(t1)^2); Q26n=(Q11-Q12-2*Q66)*(sin(t1)^3)*cos(t1)+(Q12-Q22+2*Q66)*sin(t1)*(cos(t1)^3); Q66n=(Q11+Q22-2*Q12-2*Q66)*(sin(t1)^2)*(cos(t1)^2)+ Q66*(sin(t1)^4+cos(t1)^4); % Lop 2 và lop 3 doi xung nhau t2=-pi/4; Q11n2=Q11*cos(t2)^4+ Q22*sin(t2)^4+2*(Q12+2*Q66)*(sin(t2)^2)*(cos(t2)^2); Q12n2=(Q11+Q22-4*Q66)*(sin(t2)^2)*(cos(t2)^2)+Q12*(sin(t2)^4+cos(t2)^4); Q16n2=(Q11-Q12-2*Q66)*sin(t2)*(cos(t2)^3)+(Q12-Q22+2*Q66)*(sin(t2)^3)*cos(t2); Q22n2=Q11*sin(t2)^4+ Q22*cos(t2)^4+2*(Q12+2*Q66)*(sin(t2)^2)*(cos(t2)^2); Q26n2=(Q11-Q12-2*Q66)*(sin(t2)^3)*cos(t2)+(Q12-Q22+2*Q66)*sin(t2)*(cos(t2)^3); Q66n2=(Q11+Q22-2*Q12-2*Q66)*(sin(t2)^2)*(cos(t2)^2)+ Q66*(sin(t2)^4+cos(t2)^4); % Tính ma tran Aij A11=int(Q11n,z,h/4,h/2)+int(Q11n2,z,0,h/4)+int(Q11n2,z,-h/4,0)+int(Q11n,z,-h/2,-h/4); A12=int(Q12n,z,h/4,h/2)+int(Q12n2,z,0,h/4)+int(Q12n2,z,-h/4,0)+int(Q12n,z,-h/2,-h/4); A16=int(Q16n,z,h/4,h/2)+int(Q16n2,z,0,h/4)+int(Q16n2,z,-h/4,0)+int(Q16n,z,-h/2,-h/4) A22=int(Q22n,z,h/4,h/2)+int(Q22n2,z,0,h/4)+int(Q22n2,z,-h/4,0)+int(Q22n,z,-h/2,-h/4); A26=int(Q26n,z,h/4,h/2)+int(Q26n2,z,0,h/4)+int(Q26n2,z,-h/4,0)+int(Q26n,z,-h/2,-h/4); A66=int(Q66n,z,h/4,h/2)+int(Q66n2,z,0,h/4)+int(Q66n2,z,-h/4,0)+int(Q66n,z,-h/2,-h/4); A11=vpa(A11,5); A12=vpa(A12,5); A22=vpa(A22,5); A66=vpa(A66,5); % Tính ma tran Bij B11=int(Q11n*z,z,h/4,h/2)+int(Q11n2*z,z,0,h/4)+int(Q11n2*z,z,-h/4,0)+int(Q11n*z,z,-h/2,-h/4); B12=int(Q12n*z,z,h/4,h/2)+int(Q12n2*z,z,0,h/4)+int(Q12n2*z,z,-h/4,0)+int(Q12n*z,z,-h/2,-h/4); B16=int(Q16n*z,z,h/4,h/2)+int(Q16n2*z,z,0,h/4)+int(Q16n2*z,z,-h/4,0)+int(Q16n*z,z,-h/2,-h/4); B22=int(Q22n*z,z,h/4,h/2)+int(Q22n2*z,z,0,h/4)+int(Q22n2*z,z,-h/4,0)+int(Q22n*z,z,-h/2,-h/4); B26=int(Q26n*z,z,h/4,h/2)+int(Q26n2*z,z,0,h/4)+int(Q26n2*z,z,-h/4,0)+int(Q26n*z,z,-h/2,-h/4); B66=int(Q66n*z,z,h/4,h/2)+int(Q66n2*z,z,0,h/4)+int(Q66n2*z,z,-h/4,0)+int(Q66n*z,z,-h/2,-h/4); % Tính ma tran Dij D11=int(Q11n*z^2,z,h/4,h/2)+int(Q11n2*z^2,z,0,h/4)+int(Q11n2*z^2,z,-h/4,0)+int(Q11n*z^2,z,-h/2,-h/4); D12=int(Q12n*z^2,z,h/4,h/2)+int(Q12n2*z^2,z,0,h/4)+int(Q12n2*z^2,z,-h/4,0)+int(Q12n*z^2,z,-h/2,-h/4); D16=int(Q16n*z^2,z,h/4,h/2)+int(Q16n2*z^2,z,0,h/4)+int(Q16n2*z^2,z,-h/4,0)+int(Q16n*z^2,z,-h/2,-h/4); D22=int(Q22n*z^2,z,h/4,h/2)+int(Q22n2*z^2,z,0,h/4)+int(Q22n2*z^2,z,-h/4,0)+int(Q22n*z^2,z,-h/2,-h/4); D26=int(Q26n*z^2,z,h/4,h/2)+int(Q26n2*z^2,z,0,h/4)+int(Q26n2*z^2,z,-h/4,0)+int(Q26n*z^2,z,-h/2,-h/4); D66=int(Q66n*z^2,z,h/4,h/2)+int(Q66n2*z^2,z,0,h/4)+int(Q66n2*z^2,z,-h/4,0)+int(Q66n*z^2,z,-h/2,-h/4); D11=vpa(D11,5); D12=vpa(D12,5); D16=vpa(D16,5); D22=vpa(D22,5); D26=vpa(D26,5); D66=vpa(D66,5); D1=D11; D2=-D22; D3=2*D66+D12; Ds=[D11 D12 D16;D12 D22 D26;D16 D26 D66] E1s=(A11*A22-A12^2)/A22; E1s=vpa(E1s,5); E2s=(A11*A22-A12^2)/A11; E2s=vpa(E2s,5); Gs=A66; m1=(pi^4)*(D1/a^4+2*D3/(a^2*b^2)+D2/b^4); m1=vpa(m1,5) m3=((pi^4)/16)*(E2s/b^4+E1s/a^4); m3=vpa(m3,5) % Bai toan tinh f=solve((k+m1)*f+m3*f^3-16*q0/pi^2); f=vpa(f,5) % Ham do vong w=[1 x y x^2 x*y y^2 x^3 x^2*y x*y^2 y^3 x^3*y x*y^3]; tx=diff(w,x); ty=diff(w,y); p=[w;tx;ty]; p1=subs(p,[x,y],[0,0]); p2=subs(p,[x,y],[a/4,0]); p3=subs(p,[x,y],[a/4,b/4]); p4=subs(p,[x,y],[0,b/4]); A1=[p1;p2;p3;p4]; N=p*inv(A1); N1=N(1,:); B=[diff(diff(N1,x),x);diff(diff(N1,y),y);2*diff(diff(N1,x),y)]; % Ma tran do cung phan tu cua tam composite K1=int(int(B'*Ds*B,x,0,a/4),y,0,b/4); K1=vpa(K1,5); %Ma tran do cung phan tu cua nen dàn hoi K2=k*int(int(N1'*N1,x,0,a/4),y,0,b/4); K2=vpa(K2,5); % Ma tran do cung phan tu cua ban composite tren nen dan hoi Ke=K1+K2; %Vecto tai phan tu cua tam P=int(int(N1'*q0,x,0,a/4),y,0,b/4); P=vpa(P,5); % Ma tran dinh vi I=eye(75,75); L1=[I(1:6,:);I(19:21,:);I(16:18,:)]; L2=[I(4:9,:);I(22:24,:);I(19:21,:)]; L3=[I(7:12,:);I(25:27,:);I(22:24,:)]; L4=[I(10:15,:);I(28:30,:);I(25:27,:)]; L5=[I(16:21,:);I(34:36,:);I(31:33,:)]; L6=[I(19:24,:);I(37:39,:);I(34:36,:)]; L7=[I(22:27,:);I(40:42,:);I(37:39,:)]; L8=[I(25:30,:);I(43:45,:);I(40:42,:)]; L9=[I(31:36,:);I(49:51,:);I(46:48,:)]; L10=[I(34:39,:);I(52:54,:);I(49:51,:)]; L11=[I(37:42,:);I(55:57,:);I(52:54,:)]; L12=[I(40:45,:);I(58:60,:);I(55:57,:)]; L13=[I(46:51,:);I(64:66,:);I(61:63,:)]; L14=[I(49:54,:);I(67:69,:);I(64:66,:)]; L15=[I(52:57,:);I(70:72,:);I(67:69,:)]; L16=[I(55:60,:);I(73:75,:);I(70:72,:)]; %Ma tran do cung tong the cua ban composite tren nen dan hoi Kt=L1'*Ke*L1+L2'*Ke*L2+L3'*Ke*L3+L4'*Ke*L4+L5'*Ke*L5+L6'*Ke*L6+L7'*Ke*L7+L8'*Ke*L8+L9'*Ke*L9+L10'*Ke*L10+L11'*Ke*L11+L12'*Ke*L12+L13'*Ke*L13+L14'*Ke*L14+L15'*Ke*L15+L16'*Ke*L16; Kt=vpa(Kt,5); %vecto tai tong the cua ban Pt=L1'*P+L2'*P+L3'*P+L4'*P+L5'*P+L6'*P+L7'*P+L8'*P+L9'*P+L10'*P+L11'*P+L12'*P+L13'*P+L14'*P+L15'*P+L16'*P; Pt=vpa(Pt,5); % Ban le bon canh Kt([75,74,73,71,70,68,67,65,64,63,62,61,60,58,48,46,45,43,33,31,30,28,18,16,15,14,13,11,10,8,7,5,4,3,2,1],:)=[]; Kt(:,[75,74,73,71,70,68,67,65,64,63,62,61,60,58,48,46,45,43,33,31,30,28,18,16,15,14,13,11,10,8,7,5,4,3,2,1])=[]; Pt([75,74,73,71,70,68,67,65,64,63,62,61,60,58,48,46,45,43,33,31,30,28,18,16,15,14,13,11,10,8,7,5,4,3,2,1],:)=[]; disp('chuyen vi tai cac nut chua biet la') Q=inv(Kt)*Pt; Q=vpa(Q,5) Wmax=Q(19,:) %th