Đề tài Mô phỏng sự lan truyền vết nứt trong không gian hai chiều

Science & Technology Development, Vol 13, No.K5- 2010 Trang 40 Bản quyền thuộc ĐHQG.HCM MÔ PHỎNG SỰ LAN TRUYỀN VẾT NỨT TRONG KHÔNG GIAN HAI CHIỀU Trương Tích Thiện, Trần Kim Bằng Trường Đại học Bách khoa, ĐHQG-HCM (Bài nhận ngày 28 tháng 06 năm 2010, hoàn chỉnh sửa chữa ngày 12 tháng 11 năm 2010) TÓM TẮT: Trong lĩnh vực cơ học nứt, việc dự ñoán hướng ñi của vết nứt khi xảy ra hiện tượng vết nứt lan truyền ñóng vai trò quan trọng vì việc này sẽ ñánh giá ñược vết nứt khi lan truyền liệu có xâm phạm vào những vùng quan trọng, nguy hiểm của cấu trúc hay không. Bài báo cáo này sẽ ñề cập tới ba lý thuyết dự ñoán hướng lan truyền của vết nứt là thuyết ứng suất pháp theo phương tiếp tuyến cực ñại, thuyết suất giải phóng năng lượng cực ñại và thuyết mật ñộ năng lượng biến dạng cực tiểu. Đồng thời, chương trình FRANC2D sẽ ñược sử dụng ñể mô phỏng sự lan truyền của vết nứt dựa trên cơ sở các lý thuyết trên. Từ khóa: Cơ học nứt, vết nứt, lan truyền, chương trình FRANC2D. 1. GIỚI THIỆU Hệ số cường ñộ ứng suất là thông số vô cùng quan trọng trong cơ học nứt, nói lên mức ñộ tập trung ứng suất tại ñỉnh vết nứt. Trong không gian 3 chiều, các hệ số cường ñộ ứng suất KI, KII, KIII, ñặc trưng cho 3 sự chuyển vị ñộc lập của vết nứt gồm dạng mở rộng (opening – mode I), dạng trượt (sliding – mode II) và dạng xé (tearing – mode III). Khi dự ñoán hướng lan truyền của vết nứt hai chiều, 3 phương pháp σθθmax, Smin, Gmax ñều sử dụng 2 thông số quan trọng chính là KI và KII ñể tính toán góc uốn của vết nứt. Bài báo này sẽ trình bày cơ sở lý thuyết của các phương pháp này và một số mô hình vết nứt lan truyền ñơn giản ñược tham khảo từ các tài liệu khác. 2. PHƯƠNG PHÁP DỰ ĐOÁN HƯỚNG LAN TRUYỀN CỦA VẾT NỨT 2.1. Thuyết ứng suất pháp theo phương tiếp tuyến cực ñại σθθmax Các biểu thức dạng hỗn hợp của trường ứng suất ñàn hồi quanh ñỉnh vết nứt khi ñược biểu diễn theo tọa ñộ cực như sau 21 3cos [ 1 sin sin 2 tan ] 2 2 2 22rr I II II K K K r θ θ θ σ θ pi   = + + −    (1) 21 3cos [ cos sin ] 2 2 22 I II K K r θθ θ θ σ θ pi = − (2) ( )1 cos [ sin 3cos 1 ] 22 2r I II K K r θ θ σ θ θ pi = + − (3) Trong ñó, KI, KII là hai hệ số cường ñộ ứng suất ñặc trưng cho hai dạng chuyển vị ñộc lập của vết nứt là dạng mở rộng (mode I) và dạng trượt (mode II). TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 13, SỐ K5 - 2010 Bản quyền thuộc ĐHQG-HCM Trang 41 Thuyết ứng suất pháp theo phương tiếp tuyến cực ñại σθθmax bậc nhất ñối với vật liệu ñẳng hướng khẳng ñịnh vết nứt sẽ phát triển theo hướng vuông góc với ứng suất pháp theo phương tiếp tuyến cực ñại. Thuyết này ñược Sih và Erdogan ñưa ra vào năm 1963. Hình 1. Ứng suất pháp theo phương tiếp tuyến cực ñại trong hệ tọa ñộ cực. Đạo hàm biểu thức (2) theo biến θ và gán bằng 0. 0θθσ θ ∂ = ∂ (4) Sau khi sắp xếp lại và ñặt θ = ∆θc, biểu thức (2) sẽ có dạng sau sin 3cos 1 cII I c K K θ θ − ∆ = ∆ − (5) Giải phương trình (5) theo biến ∆θc, ta sẽ tính ñược góc uốn của vết nứt. Theo sự tham khảo từ tài liệu [3], dựa trên thuyết ứng suất pháp theo phương tiếp tuyến cực ñại, góc uốn của vết nứt còn có thể ñược tính từ công thức sau 2 4 2 2 2 2 3 8 arccos 9 II I I II c I II K K K K K K θ  + +  ∆ = +   (6) Theo công thức (6), ∆θC 0. Ngoài ra, góc uốn của vết nứt còn có thể ñược tính toán theo công thức ñược tham khảo từ tài liệu [4] như sau ( )2 2 /2arctan 1 1 8 / II I c II I K K K K θ   − ∆ =  + +  (7) Theo công thức (7), nếu KII = 0 thì ∆θC = 0 ( dạng mở rộng thuần túy). Nếu KII > 0 thì góc uốn của vết nứt ∆θC < 0. Nếu KII < 0 thì góc uốn của vết nứt ∆θC > 0. 2.2. Thuyết mật ñộ năng lượng biến dạng cực tiểu Smin Thuyết này ñược Sih ñưa ra vào năm 1974. Sih ñã phát triển công thức tính mật ñộ năng lượng biến dạng S theo hệ số cường ñộ ứng suất KI và KII như sau 2 2 11 12 222I I II IIS a K a K K a K= + + (8) Với Science & Technology Development, Vol 13, No.K5- 2010 Trang 42 Bản quyền thuộc ĐHQG.HCM ( )( )11 1 1 cos cos16a θ κ θµ=  + −   (9) ( )12 1 sin 2cos 116a θ θ κµ=  − −   (10) ( )( ) ( )( )22 1 [ 1 1 cos 1 cos 3cos 1 ]16a κ θ θ θµ= + − + + − (11) E là module ñàn hồi và ν là hệ số Possion. ( )2 1 Eµ ν = + 3 4κ ν= − trong trường hợp biến dạng phẳng. 3 1 ν κ ν − = + trong trường hợp ứng suất phẳng. Vết nứt sẽ phát triển theo hướng θ = ∆θc, nơi mà mật ñộ năng lượng biến dạng ở ñó là cực tiểu. 0dS dθ = và 2 2 0 d S dθ > (12) Vết nứt bắt ñầu lan truyền khi mật ñộ năng lượng biến dạng tiến tới giá trị cực ñại S = Scr. Theo sự tham khảo từ tài liệu [5], giá trị cực ñại Scr ñược tính theo công thức sau ( )( ) 21 2 1 / 2cr ICS K Eν ν= − + (13) Với KIC là giới hạn phá hủy. 2.3. Thuyết suất giải phóng năng lượng cực ñại Gmax Thuyết này dựa trên sự tính toán của Hussain vào năm 1974. Đó là các hệ số cường ñộ ứng suất KI(θ) và KII(θ) của một vết nứt chính ban ñầu với một phần bị uốn với góc θ rất nhỏ ở ñỉnh ñược tính toán dựa theo các hệ số cường ñộ ứng suất KI và KII của vết nứt thường. ( ) ( ) 3cos sin 2I I II K g K Kθ θ θ θ = +    (14) ( ) ( ) 3cos sin 2II II I K g K Kθ θ θ θ = −    (15) ( ) 224 1 /3 cos 1 /g θ piθ piθ θ θ pi −   =   + +   (16) Hình 2. Vết nứt chính ban ñầu với một phần bị uốn với góc θ. Theo biểu thức tổng quát của Irwin, suất giải phóng năng lượng G cho vết nứt ban ñầu với một phần bị uốn với góc θ sẽ như sau ( ) ( ) ( )( )2 21 I IIG K KEθ θ θ= +′ (17) TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 13, SỐ K5 - 2010 Bản quyền thuộc ĐHQG-HCM Trang 43 Với ( )21 EE ν ′ = − cho biến dạng phẳng. E E′ = cho ứng suất phẳng. Kết hợp với các biểu thức (14), (15), (16), biểu thức (17) trở thành ( ) ( ) ( ) ( )2 2 2 2 21 [ 1 3cos 8sin cos 9 5cos ]4 I I II IIG g K K K KEθ θ θ θ θ θ= + − + −′ (18) Góc lan truyền của vết nứt ñược tìm bằng cách cực tiểu hóa G(θ). ( ) 0G θ θ ∂ = ∂ (19) Và phải thỏa mãn ñiều kiện ổn ñịnh sau ( )2 2 0 G θ θ ∂ < ∂ (20) Dạng tổng quát của biểu thức (18) có thể ñược viết gọn lại như sau: ( ) ( ) ( ) ( ) ( )2 211 22 121 [ 2 ]4 I II I IIG A K A K A K KEθ θ θ θ θ= + +′ (21) ( ) 2 11 2 12 2 22 4 3sin 2sin 2 4 5sin A A g A θ θ θ θ  −     = −      +    (22) 3. SỰ SO SÁNH GIỮA BA PHƯƠNG PHÁP DỰ ĐOÁN HƯỚNG LAN TRUYỀN CỦA VẾT NỨT Sau ñây là ñồ thị so sánh kết quả giữa thuyết ứng suất pháp theo phương tiếp tuyến cực ñại với thuyết suất giải phóng năng lượng cực ñại và thuyết mật ñộ năng lượng biến dạng cực tiểu ñược tham khảo từ tài liệu [7]. Để thuận tiện cho việc so sánh, ñặt 12 tane I II KM Kpi −   =     (23) Hình 3. Đồ thị so sánh kết quả giữa ba lý thuyết dự ñoán hướng lan truyền của vết nứt. Science & Technology Development, Vol 13, No.K5- 2010 Trang 44 Bản quyền thuộc ĐHQG.HCM Ngoài ra, bài báo cáo này xin ñược ñưa ra một kết quả so sánh khác giữa 3 phương pháp dự ñoán hướng lan truyền của vết nứt. Kết quả này ñược tham khảo từ tài liệu [2]. Hình 4. Đồ thị so sánh kết quả giữa ba lý thuyết dự ñoán hướng lan truyền của vết nứt. 4. MÔ HÌNH TÍNH TOÁN 4.1 Mô hình 1 Mô hình ñược tham khảo trong tài liệu [6] với các kích thước W = 7 (ñơn vị dài), H = 8 (ñơn vị dài), a = 3,5 (ñơn vị dài). Trường hợp ñang xét là biến dạng phẳng.Với E = 30 [(ñơn vị lực)2/(ñơn vị dài)], hệ số Poisson ν = 0,25. Ứng suất trượt τ = 1 (ñơn vị áp suất). Độ tăng trưởng vết nứt ∆a = 0,5 (ñơn vị dài). Hình 5. Tấm phẳng với một vết nứt biên và chịu ứng suất tiếp. Kết quả hình ảnh biến dạng của mô hình sau khi ñược tính toán bằng FRANC2D như sau TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 13, SỐ K5 - 2010 Bản quyền thuộc ĐHQG-HCM Trang 45 Hình 6. Kết quả biến dạng ban ñầu và sau khi vết nứt phát triển sau 7 step. So sánh kết quả biến dạng của mô hình khi vết nứt phát triển sau 7 step với kết quả tham khảo từ tài liệu [6]. Hình 7. So sánh kết quả biến dạng. So sánh kết quả tính toán hướng lan truyền của vết nứt giữa 3 thuyết σθθmax, Gmax, Smin. Science & Technology Development, Vol 13, No.K5- 2010 Trang 46 Bản quyền thuộc ĐHQG.HCM Hình 8. So sánh kết quả vết nứt lan truyền giữa ba thuyết σθθmax, Gmax, Smin. 4.2. Mô hình 2 Mô hình ñược tham khảo từ tài liệu [3]. Trường hợp ñang xét là biến dạng phẳng. Các kích thước trong hình vẽ có ñơn vị là mm. Vật liệu ñàn hồi ñẳng hướng là hợp kim nhôm 7075-T6 với E = 71,7 GPa, ν = 0,33. Chiều dài vết nứt ban ñầu a0 = 10 mm. Lực P = 20KN. Độ tăng trưởng của vết nứt ∆a = 3 mm. Hình 9. Tấm phẳng với một vết nứt biên và ba lỗ tròn. Kết quả hình ảnh biến dạng của mô hình sau khi ñược tính toán bằng FRANC2D như sau TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 13, SỐ K5 - 2010 Bản quyền thuộc ĐHQG-HCM Trang 47 Hình 10. Kết quả biến dạng ban ñầu và sau khi vết nứt phát triển sau 11 step. So sánh kết quả tính toán hướng lan truyền của vết nứt sau 11 step giữa 3 thuyết σθθmax, Gmax, Smin và các kết quả tham khảo từ tài liệu [3]. Hình 11. So sánh kết quả vết nứt lan truyền. Science & Technology Development, Vol 13, No.K5- 2010 Trang 48 Bản quyền thuộc ĐHQG.HCM Hình 12. So sánh kết quả vết nứt lan truyền giữa ba thuyết σθθmax, Gmax, Smin. Đặt Me = (2/π)tan-1(KI/KII) Bảng 1. So sánh các giá trị Me sau 11 step khi tính toán bằng ba phương pháp σθθmax, Gmax, Smin. Step σθθmax Gmax Smin 0 -0,98669 -0,98669 -0,98669 1 0,99632 0,99486 0,99632 2 0,99547 0,99614 0,99547 3 0,99238 0,99275 0,99238 4 0,98495 0,98574 0,98498 5 0,98254 0,98389 0,98254 6 0,97162 0,97303 0,97159 7 0,96291 0,96512 0,96283 8 0,95365 0,95568 0,95365 9 0,93246 0,93511 0,93235 10 0,90906 0,90128 0,90880 11 0,85618 0,85156 0,85633 TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 13, SỐ K5 - 2010 Bản quyền thuộc ĐHQG-HCM Trang 49 Hình 13. Đồ thị so sánh các giá trị Me sau 11 step khi tính toán bằng ba phương pháp σθθmax, Gmax, Smin. 5. KẾT LUẬN Các giá trị Me ñược tính toán ở mỗi step theo 3 thuyết σθθmax, Gmax, Smin ñều có giá trị xấp xỉ bằng 1 (nằm trong khoảng 0,8 – 1). Do ñó, góc uốn của vết nứt ñược tính theo 3 thuyết σθθmax, Gmax, Smin ở mỗi step có giá trị gần bằng nhau. Điều này phù hợp với ñồ thị so sánh kết quả giữa ba lý thuyết dự ñoán hướng lan truyền của vết nứt ñược tham khảo từ tài liệu [7] (Hình 3). Vì vậy, ñường ñi của vết nứt ñược mô phỏng theo 3 thuyết σθθmax, Gmax, Smin có dạng gần giống nhau. SIMULATION OF CRACK PROPAGATION IN TWO DIMENSIONAL PROBLEMS Tich Thien Truong, Kim Bang Tran University of Technology, VNU-HCM ABSTRACT: Predicting crack trajectory when crack propagation occurs plays an important role in fracture mechanics problems because this will evaluate whether important areas of structure are heavily influenced by crack propagation. This article will introduce three theories to predict crack path, including maximum tangential stress theory, maximum energy release rate theory and minimum strain energy density theory. Besides, the FRANC2D program is used to simulate the crack propagation based on three above theories. Keywords: crack trajectory, crack propagation, FRANC2D program. Science & Technology Development, Vol 13, No.K5- 2010 Trang 50 Bản quyền thuộc ĐHQG.HCM TÀI LIỆU THAM KHẢO [1]. Soheil Mohammadi, Extended Finite Element Method, Blackwell Publishing, (2008). [2]. David G. Lewicki, Crack Propagation Studies to Determine Benign or Catastrophic Failure Modes for Aerosapcas Thin-Rim Gear, Army Research Laboratory, Technical Report ARL-TR-971. [3]. E. Giner, N. Sukumar, J. E. Tarancon and F. J. Fuenmayor, An Abaqus implementation of the extended finite element method, Preprint submitted to Engineering Fracture Mechanics, (2008). [4]. N. Sukurmar and J. –H. Prevost, Modeling Quasi-Static Crack Growth with the Extended Finite Element Method. Part I: Computer Implementation, International Journal of Solids and Structures, (2003). [5]. Ali Hassan CHAHROUR and Masayasu OHTSU, Simulation of Discrete Cracking in a Concrete Gravity Dam, Vol. 16, No.2, (1994). [6]. Zhenjun Yang, Fully automatic modelling of mixed – mode crack propagation using scaled boundary finite element method, Engineering Fracture Mechanics 73, pp. 1711 – 1731, (2006). [7]. Ingraffea A. R., Lecture Notes, Cornell University, CEE 770, Fall (2007). [8]. CFG. FRANC2D Users Guide – Version 3.1, (2003).