Phân tích bài toán nứt phẳng đàn hồi tuyến tính bằng phần tử hữu hạn mở rộng tứ giác nội suy kép (xcq4)

15 Tạp chí Khoa học Lạc Hồng Số 04 Journal of Science of Lac Hong University Vol.4 (12/2015), pp. 15-20 Tạp chí Khoa học Lạc Hồng Số 4 (12/2015), trang 15-20 PHÂN TÍCH BÀI TOÁN NỨT PHẲNG ĐÀN HỒI TUYẾN TÍNH BẰNG PHẦN TỬ HỮU HẠN MỞ RỘNG TỨ GIÁC NỘI SUY KÉP (XCQ4) Optimization of linear elastic fracture mechanics by using an extended consecutive – Interpolation quadrilateral element (XCQ4) method Nguyễn Đình Dư1, Nguyễn Duy Phích2 1dinhdu85@gmail.com Khoa Kỹ Thuật Công trình Trường Đại học Lạc Hồng, Đồng Nai, Việt Nam Đến tòa soạn 10/1/2015; Chấp nhận đăng: 1/2/2015 Tóm tắt. Bài báo trình bày một phần tử hữu hạn mở rộng tứ giác (XCQ4) dựa trên thủ tục nội suy kép với ứng suất liên tục tại nút để mô phỏng trường ứng suất vùng lân cận đỉnh vết nứt hai chiều. Khác với phương pháp truyền thống, hàm xấp xỉ trong nghiên cứu này bao gồm giá trị tại nút và trung bình cộng giá trị đạo hàm của nó từ bước đầu tiên. Mục tiêu chính của bài viết này là nhằm giới thiệu một sự phát triển của phần tử CQ4 được công bố gần đây với kỹ thuật làm giàu nhằm tính chính xác hệ số cường độ ứng suất tại đỉnh nứt (SIFs). Sự chính xác của phần tử XCQ4 trong nghiên cứu này được chứng minh thông qua các ví dụ số có hình học từ đơn giản đến phức tạp. Độ chính xác cũng như tốc độ hội tụ thu được SIFS từ XCQ4 là cao hơn XQ4 truyền thống với điều kiện như nhau. Phần tử mới XCQ4 có thể được mở rộng và áp dụng cho các bài toán phức tạp hơn trong thực tế. Từ khoá: Phần tử hữu hạn mở rộng; Cơ học rạn nứt; Hệ số cường độ ứng suất; Kỹ thuật làm giàu; Nội suy kép Abstract. This work presents a novel extended 4-node quadrilateral finite element (XCQ4) method based on the consecutive- interpolation procedure (CIP) with continuous nodal stress in oder toaccurately model singular stress fields near crack tips of two-dimensional (2D) elastic cracked solids. Incontrast with conventional FEM methods, the approximation functions constructed based on CIP involve both nodal values and averaged nodal gradients as interpolation conditions. As a pioneering extension of a recently developed CQ4 element associated with enrichment method, the proposed XCQ4 extracts the stress intensity factors (SIFs) at the crack tipsprecisely. Accuracy and convergence of the SIFs results obtained by the proposed method are as high as those of the standard XQ4 solutions. Importantly, the proposed XCQ4 elementmethod is highly promising for use in other complex engineering problems. Keywords: Fracture; Extended finite element method; Stress intensity factors; Smooth nodal stress; Enrichment; Consecutive- interpolation 1. GIỚI THIỆU Mô hình hóa chính xác trường chuyển vị và biến dạng đỉnh nứt vẫn đang là một vấn đề thách thức trong cơ học rạn nứt. Việc dự đoán chính xác trường ứng suất cạnh đỉnh nứt đóng một vai trò quan trọng trong việc bảo trì, dự đoán tuổi thọ và đánh giá sự an toàn của vật liệu cũng như kết cấu trong tương lai. Trong một vài thập kỷ qua, phương pháp phần tử hữu hạn (FEM) được xem là một công cụ số hiệu quả giải quyết nhiều bài toán kỹ thuật phức tạp mà lời giải giải tích không đáp ứng được. Tuy nhiên, khi giải quyết bài toán nứt gặp nhiều khó khăn do sự liên kết cấu trúc phần tử phải được cập nhập một cách liên tục trong suốt quá trình tái chia lưới. Do đó, có nhiều nỗ lực nhằm cải tiến phương pháp FEM để phù hợp với yêu cầu kỹ thuật, một trong những phương pháp đó chính là phần tử hữu hạn mở rộng (XFEM). Trên cơ sở đó, trong phạm vi nghiên cứu này, một nghiên cứu mới cho bài toán nứt của kết cấu đàn hồi hai chiều được sử dụng bởi phần tử hữu hạn mở rộng tứ giác nội suy kép (XCQ4) được trình bày. Kết quả số thu được từ phương pháp được so sánh với XFEM, các phương pháp số khác cũng như lời giải chính xác nhằm kiểm chứng hiệu suất của XCQ4. 2. XÂY DỰNG PHẦN TỬ XCQ4 PHÂN TÍCH NỨT Về cơ bản mà nói, XCQ4 được nghiên cứu trong bài báo này là một phiên bản cải tiến của chúng tôi về CQ4 được tích hợp các chức năng làm giàu để mô phỏng chính xác trường biến dạng của vết nứt. 2.1 Phần tử CQ4 và các thuộc tính Mô tả chi tiết phần tử CQ4 có thể tìm thấy trong [1]. Tuy nhiên, để tiện theo dõi, tác giả xin trình bày ngắn gọn CQ4 được trình bày trong bài báo này. Một điểm cần nội suy có tọa độ x(x, y) trong phần tử tứ giác có bốn nút lần lượt i, j, k, m được minh họa trong Hình 1. Như mô tả trong hình vẽ, chúng ta có các miền phần tử Si, Sj, Sk và Sm là các miền chứa tất cả các phần tử có liên quan lần lượt đến các nút i, j, k và m. Như vậy, những nút hỗ trợ cho điểm x trong phần tử CQ4 bao gồm tất cả các nút của miền phần tử Si, Sj, Sk và Sm. Nguyên Đình Dư, Nguyên Duy Phích 16 Tạp chí Khoa học Lạc Hồng Số 04 Hình 1. Hình minh họa CQ4 trong 2D Theo đó, miền hỗ trợ CQ4 cho điểm x thì rộng hơn miền hỗ trợ Q4 trong FEM chuẩn. Lời giải gần đúng tại điểm x có thể viết như sau: (1) Trong phương trình (1), hàm dạng nội suy kép được xác định như sau [1], [2]: (2) Trong đó, dl là chuyển vị tại nút, còn là hàm dạng của nút i và ns là tổng số nút hỗ trợ có liên quan đến điểm x. Trong công thức nội suy, việc xây dựng đạo hàm trung bình tại nút i (những nút khác tương tự) được thực hiện như [1] [2], có thể viết như sau: (3) Trong phương trình (3), thành phần là đạo hàm của được tính trong phần tử e, và we là hàm trọng số của phần tử e ÎSi, nó được định nghĩa như sau: , với e ÎSi (4) và De là diện tích của phần tử e. Trong phương trình (2), các hàm số i, ix và iylà các đa thức cơ sở có liên quan đến nút i phải thỏa mãn điều kiện sau: i (xl) = il, i,x (xl) = 0 , i,y (xl) = 0 ix (xl) = 0 , ix,x (xl) = il, ix,y (xl) = 0 (5) (13) iy (xl) = 0 , iy,x (xl) = 0 , iy,y (xl) = il Trong đó, l lần lượt là một trong các chỉ số i, j, k và m, và (6) Cần lưu ý rằng, các điều kiện trên cũng cần áp dụng tương tự cho cho các hàm số khác như j , jx , jy; k , kx , kyvà m , mx , my. Sau cùng, các hàm đa thức cơ sở trong CQ4 được định nghĩa như sau: (7) Trong phương trình (7), p = 1/2 và các hàm số j , jx , jy; k , kx , ky ; m , mx , my cũng được thực hiện tương tự bằng cách xoay vòng các chỉ số i, j, k và m. Theo đó, Li, Lj, Lk và Lm là hệ trục tọa độ diện tích của điểm cần nội suy x trong phần tử tứ giác i, j, k và m. Tất cả các tính chất của phần tử CQ4 có thể được tìm thấy trong [1] và một điểm quan trọng cần nhắc đến ở đây đó là hàm dạng của CQ4 thì trơn và mịn hơn Q4 chuẩn được thể hiện trong Hình 2. Hình 2. Hàm dạng của Q4 (a), CQ4 (b) trong 2D. 2.2Mở rộng xấp xỉ CQ4 cho nứt Cũng giống như ý tưởng chủ đạo của XFEM [3] là sử dụng một chuyển vị xấp xỉ có khả năng mô tả chính xác sự bất liên tục và kỳ dị của vùng lân cận đỉnh nứt. Do đó, xấp xỉ mở rộng CQ4 của chuyển vị cho vết nứt có thể viết như sau: (a) (b) Phân tích bài toán nứt phẳng đàn hồi tuyến tính bằng phần tử hữu hạn mở rộng tứ giác nội suy kép (XCQ4) 17 Tạp chí Khoa học Lạc Hồng Số 04 (8) Trong đó, đại diện cho hàm dạng CQ4 liên quan đến nút i để nội suy trường chuyển vị tiêu chuẩn, Jcut là tập hợp các nút mở rộng thuộc phần tử bị vết nứt cắt qua, Ktip là tập hợp các nút mở rộng thuộc phần tử chứa đầu vết nứt, ui là vectơ chuyển vị của nút i. H(x) là hàm Heaviside có giá trị +1 trên đường nứt và – 1 dưới đường nứt, aj là chuyển vị bậc tự do mở rộng thuộc Jcut. Cuối cùng, là bậc tự do chuyển vị tại nút K thuộc Ktip phù hợp với hàm mở rộng tiệm cận . Đối với bài toán đàn hồi tuyến tính đẳng hướng 2D, hàm được định nghĩa như sau: (9) Trong đó, (r, ) là hệ tọa độ cực địa phương có gốc tọa độ là đỉnh nứt. 2.3 Phương trình rời rạc Các hệ thống rời rạc của phương trình cân bằng tuyến tính trong điều kiện đàn hồi biến dạng nhỏ có thể được thể hiện như sau: Kd = F (10) Trong đó, K là ma trận độ cứng tổng thể, d là vectơ chuyển vị nút tổng thể bao gồm nút tiêu chuẩn và nút mở rộng. Đối với phần tử được làm giàu, ma trận độ cứng phần tử thu được: (11) Trong khi đối với phần tử tiêu chuẩn thì (12) với ; r, s@u, a, b (13) Trong đó, B ma trận tính biến dạng, là đạo hàm đối xứng rời rạc của ma trận hàm dạng mở rộng. Các thành phần của B được cho như sau: ; ; (14) Trong phương trình (10), F là vectơ lực nút tác động từ bên ngoài và bản thân của phần tử dẫn đến vectơ lực tổng thể phần tử được viết như sau: (15) Đối với phần tử tiêu chuẩn thì (16) Các thành khác là của phần tử mở rộng, chúng được định nghĩa chi tiết như sau: (17) (18) (19) với và lần lượt là lực bản thân trên một đơn vị thể tích và lực tác dụng lên biên. 3. TÍCH PHÂN TƯƠNG TÁC VÀ CÁCH TÍNH SIFS Để thu được SIFs, dạng miền của tích phân tương tác được chọn trong bài viết này [4]. Hai trạng thái của tấm chứa vết nứt được xem xét. Trạng thái #1 là trạng thái thực tương ứng với sự làm việc hiện tại, còn trạng thái #2 là trạng thái ảo, trạng thái xấp xỉ, là trạng thái được chọn như là trường tiệm cận khu vực đầu đỉnh nứt cho mode I và mode II. Mối quan hệ giữa tích phân tương tác và SIFs đa mode như sau: (20) Trong đó, M(1,2) là tích phân tương tác và được định nghĩa như sau: (21) Trong (20), là tham số vật liệu và được định nghĩa như sau: (22) và trong (21), . Nếu trạng thái #2 được giả định là trường tiệm cận của mode I với KI(2) = 1 và KII(2) = 0, sau đó SIF mode I cho trạng thái #1 thu được như sau: (23) Ứng suất phẳng Biến dạng phẳng Nguyên Đình Dư, Nguyên Duy Phích 18 Tạp chí Khoa học Lạc Hồng Số 04 Tương tự cách tính trên, SIF mode II cho trạng thái #1 thu được như sau: (24) Tích phân M(1,2) trong (21) không phù hợp lắm với phần tử hữu hạn, gây rất nhiều khó khăn trong tính toán. Vì thế, nó được chuyển sang tích phân miền bằng cách nhân thêm một hàm trọng số biên q(x), nó được minh họa như Hình 3, giá trị của nó thay đổi từ 1 khi ở gần đỉnh nứt và bằng 0 trên đường viền bên ngoài theo quy định. Cuối cùng, tích phân tương tác trong hình thể mới được xác định bởi: (25) Một vấn đề cần lưu ý thêm khi thực hiện nội suy kép trong hàm xấp xỉ, đó là miền hỗ trợ cho các phần tử thuộc J-domain trong XCQ4 là lớn hơn XQ4 chuẩn. J-domain được thể hiện trong Hình 4a cho XCQ4 và 4b cho XQ4. Do đó hệ số SIFS thu được từ XCQ4 có độ chính xác cao hơn XQ4 và được kiểm chứng trong các ví dụ số. Hình 3. Hình ảnh minh họa hàm trọng số q(x) dùng trong tích phân tương tác: (a) mặt bằng, (b) mặt đứng (a) (b) Hình 4. Hình ảnh minh họa J-domain giữa Q4 (a) và XCQ4 (b) 4. VÍ DỤ SỐ Hai ví dụ về bài toán cơ học rạn nứt đàn hồi tuyến tính hai chiều bao gồm mode đơn và hỗn hợp được xem xét. Hệ số cường độ ứng suất (SIFs) được tính bằng tích phân tương tác thu được từ XCQ4 sau đó so sánh với XQ4 tiêu chuẩn hoặc các phương pháp số hiện có khác. Các phân tích về tổng số bậc tự do, cách chia lưới quy tắc hay bất quy tắc, bán kính tích phân tương tác nhằm đánh giá mức độ hội tụ cũng như độ chính xác của phương pháp. Trong các ví dụ số, điều kiện biến dạng phẳng được áp dụng. 4.1 Bài toán 1 Trong ví dụ này, một tấm hữu hạn chứa vết nứt cạnh chịu lực cắt phân bố phân bố đều tác dụng vào cạnh biên trên như Hình 5.a. Tấm có chiều rộng W = 7, chiều cao L = 16 và chiều dài vết nứt a = 3.5. Lực cắt có cường độ t =1 được áp vào biên trên của tấm. Dưới tác dụng của lực này, tấm sẽ bị phá hoại theo hai trạng thái đó là mode 1 và mode 2. Lời giải chính xác [1] cho kết quả KI = 34.0, KII = 4.55, là điều kiện dùng để so sánh với XCQ4. (a) (b) Hình 5. Tấm hữu hạn chứa vết nứt cạnh chịu lực cắt (a) và lực kéo (b) Bảng 1. Độ hội tụ của SIFs với mật độ lưới của tấm hữu hạn chịu lực cắt R Phương Pháp Mật độ lưới 9´19 19´39 29´59 39´79 2 K1 XQ4 31.3742 33.0055 33.4473 33.6480 XCQ4 31.7094 33.1934 33.6175 33.8123 K2 XQ4 4.4075 4.4632 4.4775 4.4839 XCQ4 4.4987 4.5301 4.5635 4.5553 2.5 K1 XQ4 31.1810 32.8235 33.2693 33.4731 XCQ4 31.5365 33.0080 33.4352 33.6361 K2 XQ4 4.3953 4.4522 4.4698 4.4774 XCQ4 4.4467 4.4764 4.4838 4.4911 3 K1 XQ4 31.2015 32.8516 33.3274 33.5318 XCQ4 31.3362 33.0001 33.4795 33.6808 K2 XQ4 4.3912 4.4477 4.4664 4.4737 XCQ4 4.4265 4.4563 4.4709 4.4989 Sự ảnh hưởng của tích phân J phụ thuộc vào bán kính R được mô tả như Hình 6. Trong nghiên cứu này, bán kính J-domain R được khảo sát ở ba giá trị khác nhua bao gồm R = 2, R = 2.5, R = 3. Mật độ chia lưới hay tổng số bậc tự do là một trong những yếu tố ảnh hưởng rất lớn đến kết quả tính toán. Theo đó, bốn cách chia lưới có quy tắc 9´19; 19´39; 29´59 và 39´79 được xem xét. Kết quả tính toán của mỗi miến tích phân cùng với cách chia lưới được thể hiện trong Bảng 1. Kết quả số thu được là rất tốt, sự chính xác của hệ số SIFs càng tăng gần với lời giải chính xác khi số phần tử tăng lên. Đồng thời, giá trị trong bảng 1 cũng cho thấy rằng kết quả SIFs thu được từ XCQ4 là tốt hơn XQ4 truyền thống. Tầm ảnh hưởng của bán kính J-domain khi tính hệ số SIFs ở chế độ đa mode là nhẹ nhàng và không đáng kể. Tuy nhiên, một bán kính miền không quá lớn cũng như không quá nhỏ được chọn nhằm đảm bảo tính chính xác của phương pháp. Do đó, R = 2; 2.5 hoặc 3 có thể được dùng. Sai số tương đối của hệ số SIFs khi tính bằng XCQ4 và XQ4 ở mode I và mode II so với lời giải chính xác được phác họa ở Hình 6 cho các bán kính miền J khác nhau. Dễ dàng nhận thấy tốc độ hội tụ của cả hai mode là tương tự nhau. Điều đó có nghĩa rằng tốc độ hội tụ của XCQ4 cho các trường hợp R khác nhau thì tốt hơn XQ4 cho cả hai mode. Do đó, khi thực hiện với XCQ4 thì ít suy biến hơn XQ4 truyền thống. (a) (b) Phân tích bài toán nứt phẳng đàn hồi tuyến tính bằng phần tử hữu hạn mở rộng tứ giác nội suy kép (XCQ4) 19 Tạp chí Khoa học Lạc Hồng Số 04 Hình 6. Kết quả hội tụ của mode I và mode II thu được từ XCQ4 và XQ4 với R khác nhau 4.2 Bài toán 2 Tấm hữu hạn chịu lực kéo như sơ đồ tính thể hiện trong Hình 5.b. Lực kéo có cường độ =1 đặt tại biên trên của tấm, với điều kiện của bài toán như vậy thì chỉ có mode I bị phá hoại. Tương tự như lực cắt ở trên, ta cũng đi xét sự ảnh hưởng của bán kính miền tích phân J nhưng chỉ hai giá trị R = 2 và R = 2.5 được xem xét. Độ hội tụ của hệ số SIF vào mật độ lưới cũng được phân tích. Kết quả số hệ số SIF mode I được kiểm chứng với lời giải chính xác được cho bởi Ewalds và Wanhill [5]. (26) Hình 7. Kết quả hội tụ của mode I thu được từ XCQ4 và XQ4 với hai giá trị R Trong đó, hệ số điều chỉnh C được tính như sau: (27) Kết quả tính toán và giá trị lời giải chính xác được thể hiện trong Bảng 2. Như mong đợi, dễ dàng quan sát được sự chính sát của kết quả số thu được khi tính hệ số SIF mode I bởi XCQ4 là tốt hơn XQ4 thông thường. Đường biểu diễn kết quả số SIF của nghiên cứu này tiệm cận với lời giải chính xác khi tăng số lượng phần tử được thể hiện ở Hình 7. Cũng giống như ví dụ ở trước, sự ảnh hưởng của bán kính miền tích phân J cũng tìm thấy trong ví dụ này. Bảng 2. Độ hội tụ của KI phụ thuộc vào mật độ lưới R Phương Pháp Mật độ lưới Chính xác 9´19 19´39 29´59 39´79 2 XQ4 8.7292 9.1050 9.2137 9.2649 9.3721 XCQ4 8.7808 9.1499 9.2591 9.3105 2.5 XQ4 8.6869 9.0602 9.1681 9.2192 XCQ4 8.7362 9.1012 9.2108 9.2632 Hơn nữa, sai số tương đối của hệ số SIF mode I được rời rạc hóa theo bán kính miền R của tích phân J được minh họa trong Hình 8 rõ ràng cho thấy thêm lần nữa là XCQ4 tốt hơn XQ4. Nguyên Đình Dư, Nguyên Duy Phích 20 Tạp chí Khoa học Lạc Hồng Số 04 Hình 8. Trường ứng suất von Mises giữa XQ4 (a) và XCQ4 (b) Để tiếp tục giải thích biểu hiện sự hiệu quả của phương pháp XCQ4 trong vấn đề phân bố ứng suất. Ứng suất von Mises trong ví dụ này với kích thước chia lưới có quy tắc 39´79 được thể hiện trong Hình 8. Kết quả cho thấy rằng việc sử dụng phần tử XCQ4 với thủ tục nội suy kép được đề xuất trong nghiên cứu này cho trường ứng suất mịn màng hơn XQ4 truyền thống. 5. KẾT LUẬN Một phương pháp số mở rộng bằng cách sử dụng các hàm số cơ bản của CQ4 để mô hình hóa chính xác các vết nứt đã được phát triển bởi kỹ thuật làm giàu. Nó bao gồm hàm bất liên tục Heaviside và một trường tiệm cận gần đỉnh nứt. Nghiên cứu đã sử dụng các mở rộng của phần tử CQ4 (XCQ4), nó đặc biệt phù hợp cho mô hình vết nứt, để phân tích hệ số cường độ ứng suất và một số vấn đề trong bài toán cơ rạn nứt đàn hồi tuyến tính 2D. Trong mỗi trường hợp nghiên cứu, các SIFs và độ hội tụ của nó được tìm thấy khớp với các phân tích tương ứng như XQ4 nói chung và các phương pháp số khác nói riêng. Nói chung, các kết quả số thu được từ XCQ4 cho thấy đạt một hiệu suất cao vì nó mang lại sự chính xác cao hơn phần tử thông thường. Mặt khác, XCQ4 còn thể hiện được trường ứng suất mịn màng ở đỉnh nứt mà không thể dễ dàng thu được tương tự ở phần tử thông thường. 6. TÀI LIỆU THAM KHẢO [1] Q.T. Bui, Q.D.Vo, C.H. Zhang, D.D. Nguyen, “A consecutive-interpolation quadrilateral element (CQ4): Formulation and applications,” Finite Elem Anal Des, vol. 84, pp. 14-31, February 2014. [2] Y.T. Yang, X.H.Tang, H. Zheng, “A three-node triangular element with continuous nodal stress,” Comput Struct, vol. 141, pp. 46-58, 2014. [3] N. Moes, J. Dolbow, T. Belytschko, “A finite element method for crack growth without remeshing,” Int. J. Numer. Meth. Engng, vol. 46, pp. 131-150, 1999. [4] Q.T. Bui,C.H. Zhang, “Extended finite element simulation of stationary dynamic cracks in piezoelectric solids under impact loading,” Comput Mater Sci, vol. 62, pp. 243-257, 2012. [5] H. Ewalds andR. Wanhill, Fracture Mechanics, Edward Arnold: New York, 1993. TIỂU SỬ TÁC GIẢ Nguyễn Đình Dư Năm sinh 1985, Tuy Phước, Bình Định. Tốt nghiệp Đại học và Thạc sĩ tại Trường Đại học Bách Khoa Thành phố Hồ Chí Minh năm 2009 và 2013. Hiện anh đang là Giảng viên cơ hữu tại Khoa Kỹ thuật Công trình. Lĩnh vực nghiên cứu: Phương pháp số, Phần tử hữu hạn, Cơ học vật rắn biến dạng.