Phân tích tĩnh tấm composite có cơ tính biến thiên theo lý thuyết chuyển vị bậc 3 đầy đủ - Đỗ Văn Thơm

Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số 37, 06 - 2015 165 PHÂN TÍCH TĨNH TẤM COMPOSITE CÓ CƠ TÍNH BIẾN THIÊN THEO LÝ THUYẾT CHUYỂN VỊ BẬC 3 ĐẦY ĐỦ Đỗ Văn Thơm*, Lê Trường Sơn Tóm tắt: Bài báo phân tích tĩnh đối với tấm composite có cơ tính biến thiên (FGM) theo lý thuyết chuyển vị bậc ba đầy đủ, sử dụng phương pháp phần tử hữu hạn để đưa ra phương trình xác định chuyển vị, nội lực, ứng suất của tấm FGM chịu tải trọng tĩnh phân bố đều theo phương vuông góc với mặt phẳng trung bình của tấm. Chương trình tính toán được lập trình trên máy tính sử dụng ngôn ngữ lập trình matlab. Kết quả tính toán có so sánh để kiểm nghiệm độ tin cậy của chương trình đã lập, kết quả cho dưới dạng hình vẽ để khảo sát ảnh hưởng của một số thông số vật liệu, hình học đến chuyển vị và ứng suất của tấm FGM. Từ khóa: Tấm composite, Tấm FGM, Lý thuyết chuyển vị bậc 3. 1. ĐẶT VẤN ĐỀ Vật liệu composite có cơ tính biến thiên (FGM) được chế tạo từ hai thành phần gốm và kim loại với tỷ lệ thể tích mỗi thành phần biến đổi một cách trơn và liên tục từ mặt này sang mặt kia theo chiều dày h, tạo ra vật liệu mới có tính năng ưu việt hơn hẳn các vật liệu ban đầu. Hiện nay, vật liệu FGM được sử dụng rộng rãi trong nhiều lĩnh vực như cơ khí, điện tử, y học, hạt nhân, vũ trụ,... Việc nghiên cứu về vật liệu FGM đã thu hút nhiều sự quan tâm của các nhà khoa học trong và ngoài nước như Đào Huy Bích và các cộng sự (2011), Đào Văn Dũng (2012), Nguyễn Đình Đức (2012), Trần Minh Tú (2014), Nguyễn Đình Kiên (2014), Reddy (2000), Mantari and Guedes Soares (2013), Sidda Redddy (2014),... Các nghiên cứu trước đây đều đi sâu vào phân tích phản ứng của tấm FGM theo lý thuyết cổ điển, lý thuyết chuyển vị bậc nhất và bậc 3 không đầy đủ, và nhiều tác giả đã đi tìm lời giải bằng phương pháp giải tích, các nghiên cứu theo lý thuyết chuyển vị bậc 3 đầy đủ theo phương pháp phần tử hữu hạn vẫn còn hạn chế.Trong bài báo này, tác giả trình bày thuật toán phần tử hữu hạn (PTHH) phân tích tĩnh tấm FGM theo lý thuyết chuyển vị bậc 3 đầy đủ. Kết quả nghiên cứu có thể tham khảo khi tính toán các kết cấu làm bằng vật liệu FGM. 2. NỘI DUNG CẦN GIẢI QUYẾT 2.1. Xây dựng lý thuyết 2.1.1. Mô hình và các phương trình cơ bản Xét một tấm chữ nhật làm bằng vật liệu composite có cơ tính biên thiên (hình 1), vật liệu tấm được pha trộn bởi 2 thành phần là kim loại và gốm với tỷ lệ thay đổi liên tục theo chiều dày tấm với hàm số mũ. Hình 1. Mô hình tấm FGM. Đối với vật liệu composite có cơ tính biến thiên, theo [2, 5, 6] ta có: Cơ kỹ thuật & Kỹ thuật cơ khí động lực Đ. V. Thơm, L. Tr. Sơn, “Phân tích tĩnh tấm Composite đầy đủ.” 166 1 1; 2 n c m c z V V V h          (1a) Trong đó, h - chiều dày của tấm; n - chỉ số mũ tỷ lệ thể tích ( 0n  ); Vm , Vc- tương ứng là tỷ lệ thể tích của kim loại và gốm. Biểu thức chuyển vị của một điểm M(x, y, z) có dạng đa thức bậc 3 của z và các hệ số chỉ phụ thuộc vào (x, y) như sau [7] :   0 2 3, , ( , ) . ( , ) . ( , ) . ( , )x x xu x y z u x y z x y z x y z x y      (2)   0 2 3, , ( , ) . ( , ) . ( , ) . ( , )y y yv x y z v x y z x y z x y z x y      (3)   0 2 3, , ( , ) . ( , ) . ( , ) . ( , )z z zw x y z w x y z x y z x y z x y      (4) Với 0 0 0, ,u v w là các chuyển vị thẳng của mặt trung bình; , ,x y z   là các chuyển vị góc; , ,x y z   là các chuyển vị thẳng bậc cao và , ,x y z   là các thành phần chuyển vị góc bậc cao. Tấm chịu tải trọng có phương vuông góc với mặt trung bình của tấm. Để tính toán theo phương pháp PTHH, tấm được chia thành các phần tử chữ nhật 9 nút. Theo phương pháp PTHH, đối với bài toán tĩnh học tấm FGM, ta có phương trình cân bằng [3]:     K q F (5) Trong đó,  K - ma trận độ cứng của kết cấu;  q - véctơ chuyển vị nút của kết cấu;  F - véc tơ tải trọng tổng thể. Ma trận  K , véc tơ tải tổng thể  F của kết cấu được xác định từ ma trận độ cứng phần tử  eK , véc tơ tải trọng nút phần tử  eF của tấm FGM chịu uốn. Theo [3] ta có:    eK K ;    eF F (6) Giải phương trình (5) cũng giống như các bài toán tĩnh học thông thường, được tác giả lập trình trong môi trường Matlab. 2.1.2. Xác định các ma trận phần tử Tấm được chia thành các phần tử chữ nhật 9 nút, mỗi nút có 12 bậc tự do (hình 2). Hình 2. Mô hình phần tử tứ giác 9 nút. Ta có véc tơ chuyển vị nút của phần tử [5,7] :    , , , , , , , , , , , T i i i i xi yi zi xi yi zi xi yi ziq u v w          ; i = 1- 9 (7) Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số 37, 06 - 2015 167 Véc tơ chuyển vị tại một điểm bất kỳ của phần tử  u được xác định từ các hàm dạng Ni và các thành phần chuyển vị nút của phần tử như sau :     9 1 .i i i u N q    (8) Trường biến dạng : 0 2 3 0 2 3 ' ' 0 2 3 0 2 3 ' ' ' 0 2 3 0 2 3 ' ' ' ' ; ; w ; v w ; v w ; x x x x x x y y y y y y z z z z z z yz z y yz yz yz yz xz z x xz xz xz xz xy y x xy xy xy xy u zk z z u zk z z zk z z zk z z zk z z u v zk z z                                                          (9) Trong đó, dấu phẩy là đạo hàm theo các biến ngay sau nó. 0 0 0 0 0 0 0 0 0 0 0 0 ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ; ; ; w ; w ; ; ; ; 2 ; 2 ; 2 ; ; ; ; 3 ; 3 ; 3 ; ; x x y y z z yz y y xz x x xy y x x x x y y y z z yz z y y xz z x x xy x y y x x x x y y y z z yz z y y xz z y x xy x y y x x x u v u v                                                                     ' ' ' ' ' '; ; 0; ; ; ;x y y y z yz z y xz z x xy x y y x               (10) Các thành phần biến dạng được viết lại dưới dạng ma trận như sau :         9 0 1 1 1 1 . [ ]i i e i L u L N q B q     (11a)         9 2 2 2 1 . [ ]i i e i L u L N q B q     (11b)         9 3 3 3 1 . [ ]i i e i L u L N q B q     (11c)         9 4 4 4 1 . [ ]i i e i L u L N q B q     (11d)         1 1 1 9 0 ' ' ' 1 [ ]i i e i L u L N q B q     (11e)         2 2 9 ' ' ' ' 2 1 . [ ]i i e i L u L N q B q     (11f)         9 ' ' ' ' 3 3 3 1 . [ ]i i e i L u L N q B q     (11g)         9 ' ' ' 4 4 4 1 ' [ ]i i e i L u L N q B q     (11h) Trong đó : [ iB ] = Li 9 1 i i N   ; [ 'iB ] = 'iL 9 1 i i N   (12) Và         1 2 9, ,..., TT T T e q q q q - ma trận chuyển vị nút của phần tử. Với Li, Li’ – các ma trận toán tử vi phân và Ni là các hàm dạng được xác định theo biểu thức :      21 2 1 1 1 1 ; 1 1 ; 4 2 N r s rs N r s s       (13a) Cơ kỹ thuật & Kỹ thuật cơ khí động lực Đ. V. Thơm, L. Tr. Sơn, “Phân tích tĩnh tấm Composite đầy đủ.” 168        23 4 1 1 1 1 ; 1 1 ; 4 2 N r s rs N r s r       (13b)      25 6 1 1 1 1 ; 1 1 ; 4 2 N r s rs N r s s      (13c)       27 8 1 1 1 1 ; 1 1 ; 4 2 N r s rs N r s r        (13d)   2 29 1 1N r s   (13e) Trong đó, r và s là tọa độ trong hệ tọa độ tự nhiên Trường ứng suất [7]:  [ ] [ ]C  (14) Trong đó : [C](6x6) là ma trận đối xứng có các thành phần Cij (i, j = 1-6) được xác định như sau : 2 11 22 33 12 13 232 3 2 3 44 55 66 14 15 16 24 25 26 34 35 36 45 46 56 (1 ). ( ) .(1 ). ( ) ; ; 1 3 2 1 3 2 ( ) ; 2(1 ) 0; E z E z C C C C C C E z C C C C C C C C C C C C C C C                                     (15a)     2 2 n c m m z h E z E E E h         (15b) ,c mE E - mô đun đàn hồi của gốm và kim loại Năng lượng biến dạng đàn hồi:       1 1 [ ]. .[K ]. 2 2 e T T e e ee e V U dV q q   (16) Trong đó, ma trận độ cứng phần tử [Ke] được xác định như sau : 1 1 1 2 1 3 1 4 2 1 2 2 2 3 2 4 3 1 3 2 3 3 3 4 4 1 4 2 4 [B ] [A][B ]+[B ] [B][B ]+[B ] [D][B ]+[B ] [E][B ]+ +[B ] [B][B ]+[B ] [D][B ]+[B ] [E][B ]+[B ] [F][B ]+ +[B ] [D][B ]+[B ] [E][B ]+[B ] [F][B ]+[B ] [G][B ]+ +[B ] [E][B ]+[B ] [F][B ]+[B ] [G][B [K ]= T T T T T T T T T T T T T T T e 3 4 4 1 1 1 2 1 3 1 4 2 1 2 2 2 3 2 4 3 1 3 2 3 3 3 ]+[B ] [H][B ]+ +[B '] [A'][B '] [B '] [B'][B ']+[B '] [D'][B ']+[B '] [E'][B ']+ +[B '] [B'][B '] [B '] [D'][B '] [B '] [E'][B '] [B '] [F'][B '] [B '] [D'][B '] [B '] [E'][B ']+[B '] [F'][B '] [B '] T T T T T T T T T T T T         4 4 1 4 2 4 3 4 4 . [G'][B '] [B '] [E'][B '] [B '] [F'][B '] [B '] [G'][B '] [B '] [H'][B '] e e S T T T T T dS                                (17) Với : Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số 37, 06 - 2015 169     11 12 13 16 /2 22 23 26 33 36/2 66 2 3 4 5 6 , , , , , , ( ). 1, , , , , , h h C C C C C C C A B D E F G H E z C C sym C z z z z z z dz               (18)     /2 44 45 45 55/2 2 3 4 5 6 ', ', ', ', ', ', ' ( ). 1, , , , , , h h C C A B D E F G H E z C C z z z z z z dz          (19) Véc tơ tải phần tử do tải trọng ngoài  eF được xác định theo [3]:      . ( ) . e e e S F N P t dS  (20) trong đó, P(t) là tải trọng phân bố đều vuông góc với mặt phẳng trung bình của tấm Các tích phân trong công thức (17) và (20) được tính bằng tích phân số cầu phương Gauss như trong [3]. Sau khi tập hợp ma trận theo công thức (6) ta được ma trận độ cứng tổng thể kết cấu [K], véc tơ tải tổng thể {F}, ta xây dựng được hệ phương trình (5). Chương trình tính toán được lập trình bằng ngôn ngữ matlab có tên là Tinh_bac_cao_FGM. 2.2. Ví dụ số 2.2.1. Kiểm tra độ tin cậy Để kiểm tra tính chính xác và độ tin cậy của kết quả nghiên cứu, ta xét một ví dụ như sau: Tấm FGM hình chữ nhật có tỷ lệ các cạnh a và b với a = b, chiều dày h = a/10, liên kết tựa đơn tất cả các cạnh và có các thông số vật liệu [4]: - Nhôm (Al): mE = 70.10 9 (N/m2), m =2702 (kg/m 3); 0,3m  - Nhôm ôxit (Al2O3) : cE = 380.10 9 (N/m2), c =3800 (kg/m 3); 0,3c  - Tấm chịu tải trọng phân bố đều với cường độ q. Số mũ biến đổi thể tích giữa gốm và kim loại n = 1. Kết quả tính được độ võng lớn nhất (không thứ nguyên) tại điểm giữa tấm ( w = 10wmax . 3 4. / .cE h q a ) : Theo kết quả các tác giả : w = 0.9431. Theo tài liệu [4] (kết quả giải tích) : w = 0.9421. Sai khác : 0,11%. Như vậy, kết quả tính toán đảm bảo độ tin cậy. 2.2.2. Ảnh hưởng của tỷ lệ thể tích vật liệu đến chuyển vị lớn nhất, ứng suất tại điểm giữa của tấm FGM Xét tấm chữ nhật composite có cơ tính biên thiên: Vật liệu gồm nhôm và Ôxit nhôm, có kích thước, chịu lực phân bố đều, liên kết tựa đơn tại 4 cạnh như trên hình 3. Các đặc trưng cơ lý của các thành phần đã cho ở trên, tấm có kích thước rộng a = b = 1 m; chiều dày h = a/10. Tấm được chia thành các phần tử chữ nhật 9 nút. Tải trọng phân bố đều trên mặt tấm, vuông góc với mặt trung bình của tấm với cường độ 6,125.106 N/m2. Cơ kỹ thuật & Kỹ thuật cơ khí động lực Đ. V. Thơm, L. Tr. Sơn, “Phân tích tĩnh tấm Composite đầy đủ.” 170 Hình 3. Sơ đồ tấm FGM chịu tải trọng tĩnh phân bố đều. Sử dụng chương trình Tinh_bac_cao_FGM nhận được kết quả đồ thị chuyển vị lớn nhất và ứng suất tại điểm giữa tấm phụ thuộc vào chỉ số thể tích (n) như hình 4 và hình 5. 0 2 4 6 8 10 -3 -2.5 -2 -1.5 -1 -0.5 x 10 -3 n w (m ) Hình 4. Chuyển vị w nút giữa tấm FGM phụ thuộc vào n. 0 2 4 6 8 10 -5 -4 -3 -2 -1 x 10 8 n X ic h m a x ( N /m 2 ) 0 2 4 6 8 10 0 5 10 15 x 10 5 n X ic h m a z ( N /m 2 ) a b 0 2 4 6 8 10 -0.01 0 0.01 0.02 0.03 n T o x z ( N /m 2 ) 0 2 4 6 8 10 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 n T o x y ( N /m 2 ) c d Hình 5. Ứng suất tại điểm (a/2, b/2, h/2) của tấm FGM phụ thuộc vào n a, ứng suất x ; b, ứng suất z ; c, ứng suất xz ; d, ứng suất xy . Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số 37, 06 - 2015 171 • Nhận xét: Chuyển vị lớn nhất tại giữa tấm tăng lên khi chỉ số mũ thể tích n tăng, điều này là hoàn toàn hợp lý vì khi n tăng, có nghĩa là thành phần kim loại trong tấm tăng, thành phần gốm giảm nên “độ cứng” nói chung của tấm giảm xuống, nên độ võng sẽ tăng lên. Ứng suất pháp x tại nút giữa tấm tăng lên khi n tăng, ứng suất z lúc tăng lúc giảm khi tăng n, và các ứng suất pháp này tăng mạnh nhất khi n tăng từ 0 - 2. Ứng suất tiếp thì lúc tăng lúc giảm khi n biến đổi. 2.2.3. Ảnh hưởng của chiều dày tấm đến chuyển vị lớn nhất, ứng suất tại điểm giữa tấm FGM: Tấm FGM có kích thước a x b như trên, bây giờ ta khảo sát tấm FGM với số mũ thể tích n = 2, cho chiều dày tấm thay đổi từ a/20 đến a/5, ta nhận được đồ thị biểu diễn chuyển vị lớn nhất tại giữa tấm theo chiều dày như hình 6, đồ thị biểu diễn ứng suất pháp , x z  , và ứng suất tiếp , yz xz  ở mặt dưới của tấm tại nút giữa tấm như hình 7 dưới đây: 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 -0.015 -0.01 -0.005 0 h(m) w (m ) Hình 6. Chuyển vị tại giữa tấm FGM phụ thuộc vào chiều dày tấm. 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 -15 -10 -5 0 x 10 8 h(m) X ic h m a x ( N /m 2 ) 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 1 1.5 2 2.5 3 3.5 4 x 10 6 h (m) X ic h m a z ( N /m 2 ) a b 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0 0.05 0.1 0.15 0.2 h (m) T o x z ( N /m 2 ) 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 -0.3 -0.2 -0.1 0 0.1 0.2 h(m) T o x y ( N /m 2 ) c d Hình 7. Ứng suất tại điểm (a/2, b/2, h/2) của tấm FGM phụ thuộc vào chiều dày tấm h: a, ứng suất x ; b, ứng suất z ; c, ứng suất xz ; d, ứng suất xy . Cơ kỹ thuật & Kỹ thuật cơ khí động lực Đ. V. Thơm, L. Tr. Sơn, “Phân tích tĩnh tấm Composite đầy đủ.” 172 • Nhận xét: Chuyển vị lớn nhất tại giữa tấm giảm đi khi chiều dày tấm tăng lên, sự thay đổi rõ rệt nhất khi chiều dày tấm biến đổi trong khoảng a/20 đến a/10. Sự thay đổi nhiều của cả ứng suất pháp và ứng suất tiếp khi chiều dày tấm biến đổi trong khoảng a/20 đến a/10. 3. KẾT LUẬN Tính toán tấm FGM theo lý thuyết chuyển vị bậc cao không có hệ số hiệu chỉnh cắt nên kết quả có độ chính xác cao hơn các lý thuyết cổ điển và lý thuyết chuyển vị bậc nhất. Bài báo đã xây dựng thuật toán phần tử hữu hạn và khảo sát một số ví dụ cụ thể để xác định độ võng lớn nhất, ứng suất pháp, ứng suất tiếp của tấm FGM chịu tải trọng tĩnh. Hướng phát triển tiếp theo có thể tính toán các dạng dao động riêng và nghiên cứu động lực học tấm FGM theo lý thuyết chuyển vị bậc 3 đầy đủ. Lời cảm ơn: Nhóm tác giả cảm ơn sự giúp đỡ về ý tưởng khoa học của GS,TS Hoàng Xuân Lượng, PGS.TS Phạm Tiến Đạt. TÀI LIỆU THAM KHẢO [1]. Hoàng Xuân Lượng, Phạm Tiến Đạt, “Nghiên cứu bài toán tấm nhiều lớp composite bằng phương pháp số”, Tuyển tập CTKH Hội nghị Cơ học vật rắn biến dạng toàn quốc lần thứ V, Hà Nội, 29-30/11/1996. [2]. Vũ Hoài Nam, Nguyễn Thị Phương, “Phân tích tĩnh tấm phẳng và panel trụ composite có cơ tính biến thiên”, Tuyển tập CTKH Hội nghị Cơ học vật rắn biến dạng toàn quốc lần thứ X, Thái Nguyên, 12-13/11/2010. [3]. Chu Quốc Thắng, “Phương pháp phần tử hữu hạn”, Nxb KH&KT, Hà Nội, 1997. [4]. Tahar Hassaine Daouadji, Abdeouahed Tounsi, Lazreg Hadji, “A theoretical analysis for static and dynamic behavior of functionally graded plates”, Materials physics and Mechanics 14, 2012, pp. 110-128. [5]. Reddy J. N, “Analysis of Functionally Graded Plates”, J. Number Method Eng. Theory and Analysis, CRC, Vol. 47, 2000, pp. 663-684. [6]. Tsung-Lin W.U, K. K. Shukla, “Nonlinear Static and Dynamic Analysis of Functionally Graded Plates”, J. of Applied Mechanics and Engineering, Vol. 11, No. 3, 2006, pp. 679. [7]. Ngô Như Khoa, “Mô hình hóa và tính toán số vật liệu kết cấu tấm composite lớp”, Luận án Tiến sĩ kỹ thuật, Hà Nội. ABSTRACT STATIC ANALYSIS OF FUNCTIONNALLY GRADED COMPOSITE PLATES BASE ON THIRD-ORDER SHEAR DEFORMATION Theoretical formulation, finite element method’ solution of rectangular plates based on higher order shear deformation model are presented for static analysis of functionally graded plates (FGPs) to setting the figure, get the equations of emotion and stress for FGM plate with uniform load. Using Matlab code to solve that equations and research for dynamic with varying of mechanical properties and geometrics of this plate. Keywords: Static analysis, Higher order shear deformation, Functionally graded plates. Nhận bài ngày 09 tháng 12 năm 2014 Hoàn thiện ngày 12 tháng 02 năm 2015 Chấp nhận đăng ngày 15 tháng 6 năm 2015 Địa chỉ: Khoa Cơ khí - Học viện Kỹ thuật quân sự; *Email: promotion6699@gmail.com