Đề tài Ứng dụng phương pháp phần tử hữu hạn mở rộng trong việc tính hệ cường độ ứng suất

TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 13, SỐ K5 - 2010 Bản quyền thuộc ĐHQG-HCM Trang 51 ỨNG DỤNG PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN MỞ RỘNG TRONG VIỆC TÍNH HỆ CƯỜNG ĐỘ ỨNG SUẤT Vũ Công Hòa, Nguyễn Công Đạt 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 18 tháng 10 năm 2010) TÓM TẮT: Phương pháp phần tử hữu hạn là một phương pháp số ñược ứng dụng rất hữu hiệu trong cơ học khi dự ñoán và mô hình hóa ứng xử cơ học của vật liệu và của kết cấu. Tuy nhiên trong một số trường hợp phương pháp phần tử hữu hạn trở nên phức tạp như việc mô phỏng sự di chuyển của những miền không liên tục, dẫn ñến việc chia lại lưới phần tử. Phương pháp phần tử hữu hạn mở rộng (PP-PTHHMR) cho ta một cách thức mới trong việc mô hình hóa vết nứt trên nền tảng phương pháp phần tử hữu hạn. Phương pháp này cho phép vết nứt ñược thể hiện một cách ñộc lập với lưới phần tử, do ñó không cần phải chia lại lưới phần tử khi mô hình vết nứt lan truyền. Bài báo này ñề cập tới việc hiện thực hóa phương pháp phần tử hữu hạn mở rộng trong tính toán hệ số mật ñộ ứng suất, một tham số quan trọng trong việc dự ñoán ñược hướng của vết nứt ngay khi vết nứt không còn phát triển. Từ khóa: Phương pháp phần tử hữu hạn, Phương pháp phần tử hữu hạn mở rộng, hệ số cường ñộ ứng suất, Abaqus. 1. GIỚI THIỆU Trong những năm gần ñây PP-PTHHMR xuất hiện như một kỹ thuật hiệu quả trong việc phân tích những vấn ñề của vết nứt. Nó ngày càng ñược sử dụng rộng rãi như một phương pháp khả thi trong mô hình vết nứt phát triển dưới giả thuyết của cơ học rạn nứt ñàn hồi tuyến tính [1, 2, 3, 4]. Nguyên tắc của PP- PTHHMR ở chỗ kết hợp những hàm mở rộng vào những phần tử suy biến ñể tính chuyển vị ở gẩn ñỉnh vết nứt. So sánh với PP-PTHH cổ ñiển, PP- PTHHMR cung cấp những thuận lợi trong việc mô phỏng sự lan truyền của vết nứt. Phương pháp này dựa trên sự mở rộng của bậc tự do của những nút bị chia cắt bởi vết nứt. 2. PHƯƠNG TRÌNH CƠ BẢN Khảo sát một miền Ω có biên là Г bao gồm Гu , Гt , Гc với Г = Г u U Гt U Гc Hình 1. Trạng thái cân bằng của vật có vết nứt Với: Г u là biên của chuyển vị, Гt là biên của ngoại lực, Гc là bề mặt kéo tự do (vết nứt), t là thời gian. Khi ñó phương trình cân bằng ñược viết : 0bfsÑ + = trong miền Ω (1) Điều kiện biên như sau: . t n fs = trên biên tG tΓ uΓ Γ n tf t 0= cΓ Ω Science & Technology Development, Vol 13, No.K5- 2010 Trang 52 Bản quyền thuộc ĐHQG.HCM =u u là trường chuyển vị trên biên u G . 0ns = trên biên c G với σ là tensor ứng suất , bf là lực khối , tf là ngoại lực, n là pháp vector ñơn vị. 3. XẤP XỈ TRONG PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN MỞ RỘNG (PP-PTHHMR) Ý tưởng cơ bản của PP-PTHHMR là mở rộng không gian hữu hạn phần tử bằng cách cộng thêm những hàm mở rộng. Khảo sát một ñiểm x thuộc miền phần tử, xấp xỉ chuyển vị tại ñiểm x ñược tính như sau [1]: ( ) ( ). ( ). ( ). enr e h fem enr I I J J I N J N u x u u N x u N x x aψ ∈ ∈ = + = +∑ ∑ (2) Trong (2): ( ). ; e fem I I I N u N x u ∈ = ∑ ( ). ( ). enr enr J J J N u N x x aψ ∈ = ∑ Với: ( )hu x là xấp xỉ chuyển vị tại ñiểm x; Iu là chuyển vị nút liên tục; Ja là chuyển vị nút không liên tục; ( )IN x và ( )JN x là các các hàm dạng tương ứng; ( )xψ là hàm mở rộng tại các nút không liên tục; eN là tập các nút của phần tử; enrN là tập các nút bị mở rộng. Trong trường hợp mô hình vết nứt phẳng ta có ñược xấp xỉ [2]: 4 1 ( ) ( ). ( ). . . . dis asympt h I I J j J K K K I N J N K N u x N x u N x H a N F bα α α∈ =∈ ∈ = + +∑ ∑ ∑ ∑ (3) Ở ñây: N là tập các nút không mở rộng; disN tập các nút bị chia cắt bởi vết nứt; asymptN là tập các nút chứa ñỉnh vết nứt; Kbα là bậc tự do mở rộng dưới ảnh hưởng của hàm KFα tại nút K ñược ñịnh nghĩa như sau: ( ) ( )= −K KF F x F xα α α (4) Khi ñó trường chuyển vị .      = =        Txh fem enr fem enr I I I I y u u N N u u u (5) Điều này ñược thể hiện rõ hơn thông qua phần tử tứ giác với hàm dạng tuyến tính. Đây là một phần tử thông dụng trong PP-PTHHMR vì việc tính toán dựa trên phần tử này không quá phức tạp và chính xác hơn so với phần tử tam giác nói chung. Xét một phần tử tứ giác TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 13, SỐ K5 - 2010 Bản quyền thuộc ĐHQG-HCM Trang 53 Hình 2. Phần tử tứ giác trong hệ tọa ñộ tổng thể và hệ tọa ñộ ñịa phương. Ta có các hàm dạng trong tọa ñộ phần tử tương ứng: 1 2 3 4 1 1( , ) (1 )(1 ) ; ( , ) (1 ) (1 ) 4 4 1 1( , ) (1 )(1 ) ; ( , ) (1 ) (1 ) 4 4 = − − = + − = + + = − + N N N N ξ η ξ η ξ η ξ η ξ η ξ η ξ η ξ η (6) Lúc này trường chuyển vị: = +h fem enri i iu u u (7)      = =        I I I I Txh fem enr fem enr i y u u u N N u u (8) 1 2 3 4 1 2 3 4 0 0 0 0 0 0 0 0   =     I fem N N N N N N N N N (9) 1 2 3 4 1 2 3 4 0 0 0 0 0 0 0 0   =      I enr N N N N N N N N N (10) Với: , ( ) ( ) ( )= = −I I I I i i iIN N x x xψ ψ ψ ψ (11) Khi ñó: 1 2 3 4 1 2 3 4 =  I Tfem x x x x y y y yu u u u u u u uu (12) . 1 2 3 4 1 2 3 4 =  I Tenr x x x x y y y ya a a a a a a au (13) 4. RỜI RẠC HÓA PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN MỞ RỘNG Theo thuyết cân bằng năng lượng [1]: = in extW W (14) Tương ñương: Ω Ω Γ Ω = Ω + Γ∫ ∫ ∫ b td d f u d f u dσ ε δ δ (15) Việc hiện thực hóa phương trình trên sử dụng PP-PTHHMR thu ñược phương trình sau: = hK u f (16) Với K là ma trận cứng tổng thể; hu là vector bậc tự do nút bao gồm bậc tự do mở rộng và f là vector ngoại lực. Trong PP-PTHHMR thì K , f ñược ñịnh nghĩa như sau:   =     fe- fe fe-en en- fe en-en K KK K K (17) Science & Technology Development, Vol 13, No.K5- 2010 Trang 54 Bản quyền thuộc ĐHQG.HCM Với Ω Ω − Ω Ω = Ω = Ω = Ω = Ω = ∫ ∫ ∫ ∫ T T T T fe enfe- fe fe en-en en fe enfe-en en fe en fe d d d d K B C B K B C B K B C B B C B K (18) 31 2 4[ ] ( 1,2,3,4) Γ Ω Γ Ω Γ Ω = = Γ + Ω = Γ + Ω = Γ + Ω = ∫ ∫ ∫ ∫ ∫ ∫ t t t bb b bu a T i i i i i i u t b i i i a t b i i i b t b i i i f f f f f f f f N f d N f d f N H f d N H f d f N F f d N F f dα α α α (19) Với B là ma trận ñạo hàm của hàm dạng: , , , , , , , , , , , , ( ) ( )0 0 0 0 0 ( ) 0 ( ) ( ) ( )( ) ( ) ( 1 4)             = = =                  = ÷ i x i x i x u a i i y i i y i i y i y i y i yi x i x i x N N H N F B N B N H B N F N N H N FN N H N F α α α α α α (20) Xét trong trưởng hợp phần tử tứ giác. Ta có tensor biến dạng : ( ) 2.     = =      xx h yy i xy Du x ε ε ε ε (21) Với D là toán tử ñạo hàm, khi ñó: = = h h I i iD N u Buε (22) Kết hợp hai trường hợp cơ bản và mở rộng ta có ñược :  =  I I fem enrB B B (23) 1, 2, 3, 4,1, 2, 3, 4, . 1, 2, 3, 4,1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4,1, 2, 3, 4, 1, 2, 3, 4, 0 0 0 00 0 0 0 0 0 0 0 ; 0 0 0 0I I x x x x x x x x fem enr y y y yy y y y y y y y x x x xy y y y x x x x N N N NN N N N N N N N N N N N N N N N N N N N N N N N N N N N       = =            B B (24) Công việc còn lại của việc tính toán là ñịnh nghĩa những hàm mở rộng Ψ [5]. 4.1 Hàm Heaviside ( ) ( )=x Hψ ξ Hàm của xấp xỉ ( )hu x ñược viết lại ở dạng: . ( . ( ) . ( )) ∈ ∈ = + −∑ ∑ enr h I I J J i J I N J N u N u N H N H aξ ξ (25) TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 13, SỐ K5 - 2010 Bản quyền thuộc ĐHQG-HCM Trang 55 , 1=iH tại vết nứt, , 0iH = tại những nơi khác. Suy ra công thức (20) ñược viết lại như sau: , , , , 0 0 i x a i i y i y i x N H B N H N H N H é ù ê ú ê ú = ê ú ê ú ê úë û (26) 4.2 Hàm dốc ( ) ( )=x xψ φ Đạo hàm của ( )xψ ñược tính ( ) ,, ( ) ( ( )) ( )= iix sign x xψ φ φ Đạo hàm của ( )xφ theo hai biến ,x y ñược tính như sau: 1 2 , 1, 2, 3, 4, 3 4 ( ) [ ] ( , )       = =       i i i i ix N N N N i x y φ φφ φ φ (27) 4.3 Hàm mở rộng gần ñỉnh vêt nứt ( ) ( , )x F ray q= Hàm mở rộng tại ñỉnh vết nứt ñược ñịnh nghĩa ở dạng hệ trục tọa ñộ ( , )r q gắn với ñỉnh vết nứt. Hình 3. Hệ tọa ñộ tổng thể và hệ tọa ñộ ñịa phương. ( , ) sin , cos , sin sin , sin cos 2 2 2 2 F r r r r ra q q q qq q q ì üï ïï ï = í ýï ïï ïî þ (28) Đạo hàm của ( , )F ra q trong hệ trục ( , )r q 1, 1, 2, 2, 3, 3, 4, 4, 1 1 sin cos cos sin 2 2 2 2 2 22 2 1 1 sin sin ( cos sin sin cos ) 2 2 2 22 1 1 cos sin ( sin cos cos cos ) 2 2 2 22 r r r r r rF F F F r r F F r r F F r r q q q q q q q q q q qq q q q q qq q q = = = = - = = + = = - + (29) Trong hệ trục tọa ñộ 1 2( , )x x y x X1 X2 α α Vết nứt Science & Technology Development, Vol 13, No.K5- 2010 Trang 56 Bản quyền thuộc ĐHQG.HCM 1 2 1 2 1 2 1 2 1, 1, 2, 2, 3, 3, 4, 4, 1 1 1 1 sin cos cos sin 2 2 2 22 2 2 2 1 3 1 3 1 3 sin sin (sin sin cos ) cos sin 2 2 2 22 2 2 1 3(cos cos cos ) 2 22 x x x x x x x x F F F F r r r r F F F r r r F r q q q q q q q qq q q q q q = = = = = = + = = + (30) Cuối cùng trong hệ trục tổng thể ta thu ñược 1 2 1 2, , , , , , cos sin sin cosx x x y x xF F F F F Fa a a a a aa a a a= - = + (31) Trong (31): a là góc hợp bởi vết nứt và trục x 5. TÍNH HỆ SỐ MẬT ĐỘ ỨNG SUẤT DỰA TRÊN PHƯƠNG PHÁP TÍCH PHÂN TƯƠNG TÁC Hệ số mật ñộ ứng suất là một tham số quan trong việc phân tích vết nứt phát triển. Hệ số mật ñộ ứng suất ñược ño bằng sự thay ñổi ứng suất tại vùng lân cận ñỉnh vết nứt. Vì vậy hệ số mật ñộ ứng suất có vai trò quan trọng trong việc biết ñược hướng của vết nứt ngay khi vết nứt không còn lan truyền. Phương pháp tích phân tương tác là một kỹ thuật rất hữu hiệu trong việc lập trình ñể tính hệ số mật ñộ ứng suất. Xét một vết nứt trong tọa ñộ ñề-các, với Γ là chu tuyến bao quanh ñỉnh vết nứt. Hình 4. Tích phân J xung quanh ñỉnh vết nứt. Tích phân J theo chu tuyến Γ ñược ñịnh nghĩa như sau [3]: 2 1 1 1 Γ Γ    ∂ ∂ = − Γ = − Γ   ∂ ∂   ∫ ∫ i i i i ij j u uJ Wdx T d W n d x x δ σ (32) Với .=i ij jT nσ là lực kéo trên chu tuyến ;Γ jn pháp vector ngoài của ;Γ W là mật ñộ năng lượng biến dạng. Trong phương pháp tích phân tương tác, một trường bổ trợ ñược ñặt thêm vào ñối tượng chứa vết nứt cùng với trường hiện có. Lúc này tích phân J là tổng của hai trường này. 2 1 1 ( )1 ( )( ) ( ) 2 aux aux aux aux i ij ij ij ij j ij ij j u uJ n d x σ σ ε ε δ σ σ Γ  ∂ + = + + − + Γ ∂ ∫ (33) Co n n e1 e2 m C+ C- TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 13, SỐ K5 - 2010 Bản quyền thuộc ĐHQG-HCM Trang 57 Trong (33): ijσ là các thành phần ứng suất; auxijσ là các thành phần ứng suất bổ trợ; ijε là các thành phần biến dạng; aux ijε à các thành phần biến dạng bổ trợ. Dạng rút gọn của công thức (33) là: = + +act auxJ J J M (34) 1 1 1 aux aux Mi i ij ij j j u u M W n d x x σ σ σ Γ  ∂ ∂ = + − Γ ∂ ∂ ∫ (35) và M aux auxij ij ij ijW σ ε σ ε= = (36) ( actJ là tích phân J thực; auxJ là tích phân J bổ trợ; M ñược gọi là tích phân tương tác và MW ñược gọi là năng lượng tương tác) Từ mối quan hệ của tích phân J và hệ số mật ñộ ứng suất, ta có: 2 2 ' ' 1 ( ) ; 2 ( . . ) I II aux aux I I II II J K K E M K K K K E = + = + (37) Ở ñây: IK và aux IK là hệ số cường ñộ ứng suất và hệ số cường ñộ ứng suất trạng thái bổ trợ theo dạng nứt mode I; IIK và aux IIK là hệ số cường ñộ ứng suất và hệ số cường ñộ ứng suất trạng thái bổ trợ theo dạng nứt mode II. Đối với từng dạng vết nứt ta chọn 1;auxIK = 0= aux IIK cho dạng nứt thứ nhất (mode I), và ngược lai cho dạng nứt thứ hai (mode II). Khi ñó: ' 2 = EK M (38) Trong công thức (38), E E′ = khi có trạng thái ứng suất phẳng và 1 EE ν ′ = − khi có trạng thái biến dạng phẳng. Với E là mô-ñun ñàn hồi dọc (Young’s modulus) và ν là hệ số poisson. 6. KẾT QUẢ VÀ SO SÁNH Trong phần này giới thiệu việc sử dụng phương pháp phần tử hữu hạn mở rộng ñể mô phỏng một ví dụ ñiển hình của cơ học nứt. Một tấm có kích thước (0.08 x 0.04) m như hình 5, chịu kéo với ứng suất σ = 1000 MPa cạnh dưới cố ñịnh theo phương y, các tham số vật liệu là E = 117.103 MPa, ν = 0.34. Hình 5. Tấm hình chữ nhật chịu kéo a W H σ Science & Technology Development, Vol 13, No.K5- 2010 Trang 58 Bản quyền thuộc ĐHQG.HCM Theo lý thuyết hệ số mật ñộ ứng suất KI ñược tính như sau [1]: IK C as p= (39) Trong ñó a là chiều dài vết nứt, W là chiều rộng tấm, và C là hệ số thực nghiệm [1] 2 3 4 1.12 0.231 10.55 21.72 30.39 a a a a aC W W W W W æ ö æ ö æ ö æ ö æ ö÷ ÷ ÷ ÷ ÷ç ç ç ç ç= - + - +÷ ÷ ÷ ÷ ÷ç ç ç ç ç÷ ÷ ÷ ÷ ÷ç ç ç ç çè ø è ø è ø è ø è ø (40) khi 0.015a m= thì hệ số mật ñộ ứng suất theo lý thuyết bẳng 428.3 Mpa m . Để tiện cho việc so sánh một tấm hình chữ nhật chịu kéo ñược chia với nhiều lưới phần tử khác nhau bởi phần tử tứ giác, ñể tiếp cận với kết quả chính xác. Trên hình 6 chuyển vị theo phương Y của tấm ñược giải bằng Abaqus và XFEM, một phần mềm khá mạnh trong lĩnh vực cơ học phi tuyến trên nền tảng phương pháp phần tử hữu hạn. Bảng 1 và hình 7 so sánh chuyển vị lớn nhất theo phương y (UYmax) của tấm với nhiều lưới phần tử khác nhau khi giải bằng Abaqus và XFEM. a) b) Hình 6. Chuyển vị theo phương Y của tấm giải bằng a) Abaqus và b) XFEM Bảng 1.Chuyển vị lớn nhất theo phương y Số phần tử UY max (Abaqus) UY max (PTHHMR) Sai Số (%) 171 1,801.10-3 1,5360.10-3 14,714 361 1,801.10-3 1,5643.10-3 13,143 741 1,801.10-3 1,5921.10-3 11,599 1521 1,801.10-3 1,5924.10-3 11,582 3081 1,801.10-3 1,6082.10-3 10,705 4661 1,801.10-3 1,6140.10-3 10,381 TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 13, SỐ K5 - 2010 Bản quyền thuộc ĐHQG-HCM Trang 59 Bảng 2. Ứng suất lớn nhất theo phương y Số phần tử σYY max (Abaqus) σYY max (PTHHMR) 171 9,888.103 2,8355.103 361 9,888.103 4,7754.103 741 9,888.103 6,7139.103 1521 9,888.103 7,4157.103 3081 9,888.103 8,0143.103 4661 9,888.103 10,033.103 Hình 7. Quan hệ giữa tổng số phần tử và chuyển vị lớn nhất UYmax Bảng 2 và hình 9 so sánh ứng suất lớn nhất theo phương Y (σYY max) của tấm với nhiều lưới phần tử khác nhau khi giải bằng Abaqus và XFEM. a) (b) Hình 8. Ứng suất trong tấm theo phương Y giải bằng a) Abaqus và b) XFEM Science & Technology Development, Vol 13, No.K5- 2010 Trang 60 Bản quyền thuộc ĐHQG.HCM Hình 9. Đồ thị quan hệ giữa tổng số phần tử và ứng suất lớn nhất theo trục Y (σYY max). Bảng 3 và hình 10 so sánh hệ số mật ñộ ứng suất với các lưới phần tử khác nhau và hệ số mật ñộ ứng suất theo lý thuyết Bảng 3.Hệ số mật ñộ ứng suất KI Số nút KI (lý thuyết) KI (xấp xỉ) Sai số (%) 50 4.283.102 3.942.102 7.966 200 4.283.102 4.089.102 4.532 800 4.283.102 4.217.102 1.547 1600 4.283.102 4.218.102 1.511 3200 4.283.102 4.253.102 0.695 4800 4.283.102 4.265.102 0.4144 Hình 10. Đồ thị ñánh giá sai số % của hệ số mật ñộ ứng suất KI so với lý thuyết dựa trên tổng số nút N. Tổng số phần tử σ Y Y m a x Tổng số nút N Sa i s ố % TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 13, SỐ K5 - 2010 Bản quyền thuộc ĐHQG-HCM Trang 61 Lưới (20 x 40) Lưới (20 x 20) Hình 11. Chu tuyến tích phân J Bảng 4: Quan hệ giữa bán kính chu tuyến và tổng số nút Tổng số nút Bán kính chu tuyến ( )%I I XFEM I K K K − 400 800 1600 1800 2400 3200 10.4 x10-3 7.3x10-3 5.1 x10-3 4.8 x10-3 4.1 x10-3 3.6 x10-3 2.8229 1.6036 3.1759 3.1397 2.5076 1.7101 Tich phân J Tổng số nút Giải tích XFEM Sai số (%) 400 800 1600 1800 2400 3200 1.380 1.380 1.380 1.380 1.380 1.380 1.303 1.336 1.293 1.294 1.311 1.333 5.580 3.188 6.304 6.232 5.000 3.406 Science & Technology Development, Vol 13, No.K5- 2010 Trang 62 Bản quyền thuộc ĐHQG.HCM Tông sô Node 0 500 1000 1500 2000 2500 3000 3500 Ba n ki n h 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.010 0.011 Hình 12. Quan hệ giữa tổng số nút N và bán kính chu tuyến R. 7. KẾT LUẬN Việc áp dụng phương pháp số ñể giải quyết các vấn ñề của cơ học rạn nứt là cần thiết trong thực tế. Thông qua sự hỗ trợ của máy tính và PP – PTHHMR, những mô hình vết nứt ñược giải quyết một cách thuận lợi, nhanh chóng. Ví dụ như hệ số mật ñộ ứng suất KIC trước ñây ñược tính thông qua thực nghiệm, nhưng ñiều này ñã ñược khắc phục thông qua phương pháp phần tử hữu hạn mở rộng bằng chương trình mô phỏng trên máy tính. Bài báo này dừng lại ở chỗ chỉ mô phỏng một ví dụ cơ bản của cơ học nứt. Trong các nghiên cứu kế tiếp theo sự tính toán và mô phỏng sự lan truyền các dạng mô hình vết nứt sẽ ñược tiếp tục phát triển. APPLYING OF EXTENDED FINITE ELEMENT METHOD FOR CALCULATING STRESS INTENSITY FACTOR Vu Cong Hoa, Nguyen Cong Dat University of Technology, VNU – HCM ABSTRACT: Finite element method is a very powerful numerical method to predict and model mechanical behavior of material and structure. However, in some cases finite element method is more complicated like the modeling of moving discontinuities, hence the need to update the mesh to match the geometry of discontinuity. Extended finite element method (XFEM) allows us a new technique to TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 13, SỐ K5 - 2010 Bản quyền thuộc ĐHQG-HCM Trang 63 modeling crack independently of the mesh; hence it is no need to remesh during propagation of the crack. In this paper, an extended finite element method is used to calculate stress intensity factor. It’s important parameter when we predict the direction of crack in the event of crack stops propagation. Keywords: finite element method, extended finite element method, stress intensity factor, Abaqus. TÀI LIỆU THAM KHẢO [1].Mohammadi S., Extended Finite Element Method for Fracture Analysis of structure, Blackwell Publishing, 24-115 (2008). [2].Moës N., Sukumar N., Moran B., Belytschko T., An Extended Finite Element Method for Two and Three Dimentional Crack Modeling, in ECCOMAS 2000, Barcelona, Spain, (9/2000). [3].Carlos Cueto-Felgueroso, Implementation Domain Integral Approach for J Integral Evaluations, Transactions, SMiRT 16, Washington DC, 1355, (8/2001). [4].Vinh P. N., Timon R., Stéphane B., Meshless methods: A review and computer implementation aspects, Mathematics and Computers in Simulation 79 (3), 763-813 (2008). [5].Bodas S., Extended Finite Element Method and Level Set Method with Applications to Growth of Cracks and Biofilms, Ph.D thesis, Northwestern University (2003). [6].Ferreira A. J. M., Matlab Codes for Finite Element Analysis Solids and Structure, Springer Publisher, 143-160 (2008). [7].Chessa J., Programing the Finite Element Method with Matlab, Northwestern University (2002).