Phân tích giới hạn tấm dày 5 bậc tự do sử dụng phương pháp không lưới Element Free Galerkin (EFG)

ĐỊA KỸ THUẬT SỐ 3-2017 3 PHÂN TÍCH GIỚI HẠN TẤM DÀY 5 BẬC TỰ DO SỬ DỤNG PHƯƠNG PHÁP KHƠNG LƯỚI ELEMENT FREE GALERKIN (EFG) NGUYỄN NGỌC PHÚC*, NGUYỄN HỒNG PHƢƠNG **, HỒ THỊ ĐOAN TRANG * The critical state of thick-plate using stabilized mesh-free method as a five-node plate bending element based on mindlin/reissener plate theory- element free galerkin (efg) Abstract: This papers concern about critical state of thick-plate using Stabilized Mesh-Free Method as a five-node plate bending element based on Mindlin/Reissener plate theory. Two case studies rectangular plates with lean and rigid boundary condition were considered. The result shows main failure forms by critical state which we want to know when using as mat, diaphragm wall, retaining wall, Keywords: Thick-plate; Stabilized Mesh-Free Method; Mindlin/Reissener plate theory; phần tử 5 bậc tự do; định lý cận trên. 1. GIỚI THIỆU * Trong những thập niên gần đây việc xác định tải trọng giới hạn của cơng trình ngày càng đƣợc quan tâm L thuyết phân tích giới hạn ngày càng phát triển phù hợp với các l thuyết tấm khác nhau Trong phân tích giới hạn trƣờng chuyển vị hoặc trƣờng ứng suất sẽ đƣợc rời rạc sau đĩ định l cận trên hoặc định l cận dƣới đƣợc áp dụng để xác định tải trọng giới hạn Bên cạnh đĩ các phƣơng pháp số cũng khơng ngừng đƣợc phát triển và là cơng cụ đắc lực để nâng cao hiệu quả tính tốn Một lớp phƣơng pháp số mới đƣợc phát triển trong thời gian gần đây là phƣơng pháp khơng lƣới (meshfree hay meshless) Gần đây nhiều phƣơng pháp khơng lƣới đƣợc phát triển nhƣ phƣơng pháp khơng lƣới Element Free Galerkin (EFG) khơng lƣới Local Petrov Galerkin (MLPG) khơng lƣới Radial Point Interpolation Method (RPIM) khơng lƣới Local Radial Point Interpolation Method (LRPIM) khơng lƣới * Khoa Xây dựng, Cao đẳng Xây dựng 2 ** Khoa Kiến trúc-Xây Dựng - Mỹ Thuật Ứng dụng, Đại Học Nguyễn Tất Thành Moving Kriging (MGK) Khác nhau cơ bản giữa các phƣơng pháp này là kỹ thuật nội suy cĩ nhiều kỹ thuật nội suy đƣợc áp dụng nhƣ Kernel Partical Method, Moving Least Square Approximate, Partition of Unity, Kringing Interpolation Phƣơng pháp EFG là một phƣơng pháp khơng lƣới đƣợc phát triển bởi Belytchko et al., 1994 Trong phƣơng pháp EFG xấp xỉ bình phƣơng cực tiểu MLS (Moving Least Square) đƣợc sử dụng để xây dựng hàm dạng phƣơng trình hệ thống đƣợc xây dựng thơng qua dạng yếu Galerkin những ơ nền đƣợc yêu cầu cho việc tính tích phân từng phần Khi sử dụng phƣơng pháp EFG để rời rạc trƣờng chuyển vị cho bài tốn cận trên số lƣợng biến trong bài tốn ít hơn nhiều so với khi rời rạc bằng FEM vì phƣơng pháp EFG chỉ yêu cầu một bậc tự do tại mỗi nút thay vì các bậc tự do của điểm Guass trong FEM Để đảm bảo tính chính xác cho lời giải khi sử dụng phƣơng pháp EFG tích phân nút ổn định (Stablised Confroming Nodal Intergration (SCNI)) đƣợc áp dụng để làm trơn biến dạng Khi đĩ tích phân đƣợc tính trực tiếp ĐỊA KỸ THUẬT SỐ 3-2017 4 tại các nút mà khơng sử dụng những điểm Gauss giúp giảm nhẹ chi phí tính tốn Trong bài báo này phần tử EFG cho nút 5 bậc tự do cho một nút đƣợc sử dụng Kết quả đạt đƣợc cĩ cải thiện so với các phƣơng pháp khác sử dụng cho nút 3 bậc tự do 2. CƠ SỞ LÝ THUYẾT 2.1. Lý thuyết tấm Mindlin 5 bậc tự do L thuyết cắt bậc nhất (L thuyết tấm dày Mindlin-Reissner): Cùng xem xét miền thể tích trong 2 R với mặt phẳng giữa tấm (mặt trung bình) Trƣờng chuyển vị theo l thuyết FSDT [5] gồm 5 bậc tự do nhƣ sau:                 0 0 0        , , , , , , , , , , , x y u x y z u x y z x y v x y z v x y z x y w x y z w x y ( 1) Biến dạng trong mặt phẳng đƣợc biểu hiện theo cơng thức 0        ε κ T xx yy xy z (2) Với biến dạng màng và biến dạng cong 0 0 ε us   1 2   κ β βT (3) (4) Biến dạng cắt  ε βs w (5) 2.2. Phƣơng pháp EFG Phƣơng pháp làm trơn biến dạng lần đầu tiên đƣợc trình bày bởi Chen et al (2000) và đƣợc hiệu chỉnh bằng cách sử dụng phép tích phân nút bởi Chen et al (2001):      φ , dΩ    h h ij J ij J J ε x ε x x x x ( 6) trong đĩ hijε là giá trị đƣợc làm trơn của h ijε tại nút J và υ là hàm phân phối (hàm trơn) và phải thỏa mãn những điểm sau (Chen et al 2000; You et al., 2004): φ 0 φdΩ 1    J và (7) Để đơn giản hàm υ đƣợc giả sử là những hàm nhỏ khơng đổi:   1 , aφ , 0,        J JJ J x x x x x (8) trong đĩ aJ là diện tích miền đại diện của nút J Thay phƣơng trình (8) vào (6) và áp dụng định l phân kì ta đƣợc:       , , Ω Γ 1 1 dΩ 2 1 = d 2 J J h h h ij J i j j i J h h i j j i J u u a u n u n a       ε x (9) trong đĩ ΓJ là biên của miền Ω Với xấp xỉ bình phƣơng cực tiểu của trƣờng chuyển vị dạng trơn của biến dạng cĩ thể đƣa ra nhƣ sau:             ε x x ε x x ε x h xx J h h h J yy J m b J s h xy J z                ε B B d; B d (10) trong đĩ: T 1 1 1 1 1,..., , ,..., , ,..., , ,..., , ,..., n n n x xn y ynu u v v w w       d ĐỊA KỸ THUẬT SỐ 3-2017 5 1, 2, , 1, 2, , 1, 2, , 1, 2, , 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 x x n x m y y n y y y n y x x n x                    B 1, 2, , 1, 2, , 1, 2, , 1, 2, , 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 x x n x b y y n y y y n y x x n x                    B 1, , 1 11, 2, 0 ... 0...0 ... 0 ... ......0 ... 0 0 ... 0 x n x n s ny y              B trong đĩ  là dạng trơn của  ns là số đoạn tạo nên miền nút ΩJ Bài tốn phân tích giới hạn cho tấm dày 5 bậc tự do Tiêu chuẩn phá hoại von Mises đƣợc sử dụng nhƣ sau:   T Tφ , 0b s P   σ τ σ P σ τ P τ (11) Với P là ứng suất chảy dẻo của vật liệu và các thơng số vật liệu của tiêu chuẩn von Mises trong bài tốn ứng suất phẳng đƣợc thể hiện nhƣ sau 2 1 0 3 01 1 1 2 0 , 0 32 2 0 0 6 b s               P P (12) Hàm năng lƣợng tiêu tán dẻo trên một đơn vị diện tích theo tiêu chuẩn von Mises   T TD ,P P b s   ε Q ε γ Q γ (13) Với 1 4 2 0 1 2 4 0 3 0 0 1             bb Q P và 1 1 01 3 0 1         ss Q P       2 0 0 2            ε κ Q ε κ γ Q γ / T T / , d d t P P b s t D z z z (14) Thực hiện phép đổi biến 2 t z   ta đƣợc biểu thức sau   1 0 0 1 1 2 2 0 0 1 2 2 2 2 4 2 4 2 2                                                             ε κ Q ε κ γ Q γ ε κ Q ε κ γ Q γ = T T T T , d d d d P P b s P b s t t t D t t t t t t (15) Sử dụng tích phân Guass để tính tích phân theo chiều dày tấm:   2 2 0 0 2 4 2 4 2 2                                     ε κ Q ε κ γ Q γ T T , ngnno g gP k P g k k b k k k s k k g t tt t t t D A (16) Thực hiện phép đổi biến để đƣa về bài tốn SOCP cĩ dạng sau:   2 2      ε γ, ngnno P T T k P g bg k sg k k g D A C C (17) ĐỊA KỸ THUẬT SỐ 3-2017 6 Sử dụng chéo hĩa ma trận Qb ta cĩ Q = VDV T b 2 2 4 2          C B B V D C B Q T T, gI I I bg m b sg s s tt t (18) Để thuận tiện biến thêm vào để đƣợc định nghĩa nhƣ sau: 1 4 2 5 3                      C C T T; gk gk gk bg k sg kgk gk (19) Năng lƣợng tiêu tán đƣợc viết lại:      ρ, ngnno P k P g gk k g D A (20) Bài tốn tối ƣu đƣợc phát biểu dƣới dạng hình nĩn bậc hai: s.t i=1,2,...,nno×ng        A d b ρ min , ngnno k P g i k g eq eq ii A t t (21) Với 0 0 0 W w 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 x y eq u eq v eq eq eq eq eq                       A A A A A A A ; w T eq 1 0 0 ... 0 0 0 ... 0 0 0 ... 0 0 0 ... 0 0 0 ... 0 yu v x ddd d d          b 3. VÍ DỤ SỐ Ví dụ số đƣợc thực hiện bằng ngơn ngữ lập trình Matlab Các thơng số vật liệu đƣợc xem là tấm đồng nhất với ứng suất dẻo là 250( Mpa) Ta xét tấm sàn hình vuơng với các điều kiện biên: 4 biên tựa đơn và 4 biên ngàm Kích thƣớc hình học của tấm: a = 10m (hình 1) Để xét sự hội tụ của lời giải hệ nút đƣợc bố trí đều theo hai phƣơng tăng dần nhƣ sau: 21×21(441 nút), 27×27(729 nút). Kết quả khảo sát sự ảnh hƣởng của việc thay đổi chiều dày tấm đƣợc thể hiện thơng qua bảng 1. Kết quả đã thể hiện nhƣ mong đợi khi hội tụ dần về với sự chênh lệch giữa các trƣờng hợp độ mảnh giảm dần Điều này thể hiện đƣợc sự hội tụ của phƣơng pháp 4 biên tựa đơn 4 biên ngàm Hình 1. Bài tốn tấm hình vuơng chịu tải phân bố đều Bảng 1. Hệ số tải trọng giới hạn của tấm vuơng 4 biên ngàm Hệ số độ mảnh 21x21 27x27 Nghiên cứu này Tham khảo [18] Nghiên cứu này Tham khảo [18] 1 9,03 8,93 8,87 8,86 2 17,94 17,77 17,67 17,68 4 30,77 30,75 30,18 30,6 8 39,24 40,76 39,05 40,45 ĐỊA KỸ THUẬT SỐ 3-2017 7 10 40,18 42,61 40,07 42,25 20 40,81 45,46 40,52 44,97 40 40,91 46,19 40,55 45,67 100 41,05 46,39 40,56 45,87 Bảng 1. So sánh hệ số tải trọng giới hạn bài tốn tấm vuơng 4 biên ngàm Tác giả Số BTD Phƣơng pháp 100 L t  1 L t  Nghiên cứu này 5 EFG 40,56 8,87 C.V. Le et al. (2016) [18] 3 MeshFree 45,87 8,86 C.V. Le (2013) [20] 3 ES-DSG 46,84 9,02 C.V. Le et al. (2010) [19] 3 HTC&EM 45,12 - Capsoni and Corradi (1999)[5] - FEM 46,18 - Kết quả khi xem xét tấm dày cĩ độ chênh lệch rất thấp (0 14% so với C V Le (2016)). Điều này thể hiện sự tƣơng đồng khi tính tốn l thuyết tấm dày Kết quả chênh lệch khi xem xét tấm mỏng dần (11,57 % so với C V Le et al (2016)[18]; 13,40 % C.V. Le (2013) [20]; 10,10 % C.V. Le et al. (2010) [19]; 12,17 % Capsoni and Corradi (1999)[5]). Ta thấy giá trị của phƣơng pháp này chênh lệch so với các nghiên cứu trƣớc đây là tƣơng đƣơng nhau cho thấy phƣơng pháp này cho kết quả đáng tin cậy Hình 2. Cơ cấu phá hoại tấm hình vuơng 4 biên ngàm ở trạng thái giới hạn Cơ cấu này phù hợp với phá hoại thực tế của tấm sàn hình vuơng biên ngàm chịu tải phân bố đều Cơ cấu phá hoại đƣợc thể hiện thơng qua các vị trí tập trung mật độ lớn năng lƣợng tiêu tán dẻo (ở vị trí nào cĩ năng lƣợng tiêu tán dẻo tập trung lớn thì sự phá hoại xuất hiện ở đĩ) Từ cơ cấu phá hoại ta cĩ thể thấy ở trạng thái giới hạn tấm sàn 4 biên ngàm cĩ xu hƣớng tập trung dọc theo biên, rẽ quạt xuất phát từ tâm tấm và mở rộng dần ra phía ngồi hƣớng về 4 gĩc, một điểm đáng lƣu là ở 4 gĩc khơng bị chảy dẻo, đƣờng chảy dẻo tạo thành một đƣờng cong bo trịn tại các gĩc. Kết quả khảo sát bài tốn tấm hình vuơng với 4 biên tựa chịu tải phân bố đều đƣợc trình bày trong bảng 3 Sự hội tụ của phƣơng pháp cũng đƣợc thể hiện khi việc sai số giảm dần giữa các trƣờng hợp khi độ mảnh tăng dần Bảng 2. Hệ số tải trọng giới hạn của tấm vuơng 4 biên tựa khi tấm mỏng dần Hệ số 21x21 27x27 ĐỊA KỸ THUẬT SỐ 3-2017 8 độ mảnh Nghiên cứu này Nghiên cứu này 1 8,98 8,86 2 16,38 16,31 4 20,05 20,12 8 20,72 20,78 10 20,84 20,92 20 21,04 21,14 40 21,10 21,26 100 21,18 21,34 Bảng 3. So sánh hệ số tải trọng giới hạn bài tốn tấm vuơng 4 biên tựa Tác giả Số BTD Phƣơng pháp 100 L t  1 L t  Nghiên cứu này 5 EFG 21,34 8,86 C.V. Le et al. (2016)[18] 3 MeshFree 25,22 8,86 C.V. Le (2013)[20] 3 ES-DSG 25,05 9,03 C.V. Le et al. (2010)[19] 3 HTC-EM 25,02 - C.V. Le et al. (2009)[15] 3 EFG 25,01 - Kết quả khi xem xét tấm dày cĩ độ chênh lệch rất thấp (0 % so với C V Le (2016)) Điều này thể hiện sự tƣơng đồng khi tính tốn l thuyết tấm dày Kết quả xét sự chênh lệch khi xem xét tấm mỏng dần (15 38 % so với C V Le (2016)[18]; 14 81 % so với C.V. Le (2013)[20]; 14 70 % so với C.V. Le et al. (2010)[19]; 14 67 % so với C.V. Le et al. (2009)[15]) Ta thấy giá trị của phƣơng pháp này chênh lệch so với các nghiên cứu trƣớc đây là tƣơng đƣơng nhau cho thấy phƣơng pháp này cho kết quả đáng tin cậy Do năng lƣợng dẻo tiêu tán hấp thụ do các chuyển vị trong mặt phẳng đƣợc xét đến và vì vậy hệ số tải trọng sẽ giảm so với mơ hình tấm Mindlin 3 bậc tự do Đối với tấm hình vuơng 4 biên tựa chịu tải phân bố đều ta thấy đƣờng chảy dẻo cĩ xu hƣớng hình thành gĩc 45 độ từ gĩc tấm hình thể hiện trong hình 2 và 3 Tuy nhiên điểm khác biệt khi trƣờng hợp tấm đủ dày thì cơ cấu phá hoại vẫn dọc theo biên Điều này cĩ thể giải thích là tấm dày cĩ xu hƣớng phá hoại cục bộ dọc biên nên khác đƣờng chảy dẻo của tấm mỏng ĐỊA KỸ THUẬT SỐ 3-2017 9 (a) Tấm dày (L/t=1) (b) Tấm Mỏng (L/t=100) Hình 3. Cơ cấu phá hoại của tấm hình vuơng 4 biên tựa ở trạng thái giới hạn 4. KẾT LUẬN Nghiên cứu đã trình bày phƣơng pháp phân tích giới hạn tấm dày 5 bậc tự do sử dụng phần tử EFG và kỹ thuật tích phân nút ổn định SNCI với các hình dạng tấm khác nhau (tấm hình vuơng tấm hình chữ nhật tấm hình trịn và tấm hình chữ L) và các điều kiện biên khác nhau (điều kiện biên ngàm chu vi và điều kiện biên tựa chu vi) Khi xem xét hiện tƣợng tấm mỏng dần phƣơng pháp dần hội tụ về giá trị tấm mỏng Sai số đối với các kết quả tham khảo tƣơng đối nhỏ Qua các kết quả đạt đƣợc và so sánh với các nghiên cứu của các tác giả khác phƣơng pháp đạt đƣợc độ tin cậy cao Bên cạnh đĩ sự tập trung của năng lƣợng tiêu tán dẻo phần nào giúp dự đốn đƣợc cơ cấu phá hoại của tấm dày Qua đĩ chúng ta cĩ thể cĩ các biện pháp gia cƣờng hợp lí khi phân tích kết cấu chịu lực TÀI LIỆU THAM KHẢO 1. T. T. M. Doan, C. V. Le and T. Q. Chu, "Limit load computation of Mindlin-Reissner plates using the ES-DSG method and second- order cone programming," The International Conference on Advances in Computational Mechanics (ACOME), Ho Chi Minh City, Vietnam, August 2012. 2. N. T. Nguyen, C. V. Le, T. Q. Chu, N. T. Tran and H. N. Pham,, "A locking-free stabilized meshfree method for computation of limit load of Mindlin-Reissner plates," International Conference on Green Technology and Sustainable Development, HCMC, September 2012. 3. Nguyen, Danh An; Bui, Thanh Cong; Nguyen, Hung Dang, "A recursive approach for limit analysis of frame," Proceedings of the Sixth National Conference on Solid Mechanics, Hanoi, 11 1999. 4. Nguyen, Hung Dang; Yan Ai-Min; Bui, Thanh Cong, "On the Limit and Shakedown Analysis of Plastified and Cracked Structures," Proceedings of The First Vietnam-Japan Symposium in Advances in Applied Electromagnetics and Mechanics HoChiMinh City, Vietnam, 19-21 January 1998 5. Capsoni A Corradi L “Limit analysis of plates - a finite element formulation” Structural Engineering and Mechanics 1999; 8:325–341 6. Zienkiewicz OC, Taylor RL, Too JM. ĐỊA KỸ THUẬT SỐ 3-2017 10 “Reduced integration technique in general analysis of plates and shells: simple and efficient element for plate bending ” International Journal for Numerical Methods in Engineering 1971; 3:275–290. 7. Hughes TJR, Taylor RL, Kanoknukulchai W “Simple and efficient element for plate bending” International Journal for Numerical Methods in Engineering 1977; 11:1529–1543. 8. Zienkiewicz OC, Lefebvre D “A robust triangular plate bending element of the Reissner– Mindlin type” International Journal for Numerical Methods in Engineering 1988; 26:1169–1184. 9. Lee SW Wong C “Mixed formulation finite elements for mindlin theory plate bending” International Journal for Numerical Methods in Engineering 1982; 18:1297–1311. 10. Simo JC Rifai MS “A class of mixed assumed strain methods and the method of incompatible modes” International Journal for Numerical Methods in Engineering 1990; 29:1595–1638. 11. Bathe KJ Dvorkin EN “A four-node plate bending element based on Mindlin/Reissener plate theory and a mixed interpolation” International Journal for Numerical Methods in Engineering 1985; 21:367–383. 12. Bletzinger KU BischoffM Ramm E “A unified approach for shear-locking free triangular and rectangular shell finite elements” Computers and Structures 2000; 75:321–334. 13. Le CV, Nguyen-Xuan H, Askes H, Rabczuk T, Nguyen-Thoi T “Computation of limit load using edge-based smoothed finite element method and second-order cone programming” International Journal of Computational Methods 2013; 10(1):1340004 14. Le CV, Nguyen-Xuan H, Askes H, Bordas S, Rabczuk T, Nguyen-Vinh H “A cellbased smoothed finite element method for kinematic limit analysis” International Journal for Numerical Methods in Engineering 2010;83:1651–74. 15. Le, C.V. and Gilbert, M. and Askes, H., "Limit analysis of plates using the EFG method and second-order cone programming," International journal for numerical methods in engineering, vol. 78, no. 13, pp. 1532--1552, 2009. 16. Belytschko, T. and Lu, Y.Y. and Gu, L., "Element-free Galerkin methods," International journal for numerical methods in engineering, vol. 37, no. 2, pp. 229--256, 2005 17. Hodge, Philip Gibson and Belytschko, Ted, "Numerical methods for the limit analysis of plates," Journal of Applied Mechanics, vol. 35, p. 796, 1968. 18. Le C V & Chu T Q Plastic “Collapse Analysis of Mindlin-Reissner Plates Using a Stabilized Mesh-Free Method” International Journal of Computational Methods, vol 13, 1650004, 2016. 19. Le, C.V. and Nguyen-Xuan, H. and Nguyen-Dang, H., "Upper and lower bound limit analysis of plates using FEM and second-order cone programming," Computers \& structures, vol. 88, no. 1, pp. 65--73, 2010. 20. Le C V “A stabilized discrete shear gap finite element for adaptive limit analysis of Mindlin–Reissner plates” International Journal for numerical methods in engineering Int. J. Numer. Meth. Engng 2013; 96:231–246. 21. Hopkins, Harry Geoffrey and Wang, Alexander Jen “Load-carrying capacities for circular plates of perfectly-plastic material with arbitrary yield condition” Journal of the Mechanics and Physics of Solids, vol. 3, no. 2, pp. 117--129,1955. ĐỊA KỸ THUẬT SỐ 3-2017 11 Người phản biện: TS NGUYỄN VIỆT TUẤN