Phân tích động lực học vết nứt trong vật liệu lẫn hạt cứng và lỗ rỗng bằng phương pháp phần tử hữu hạn nội suy liên tiếp mở rộng

2261(8) 8.2019
Khoa học Tự nhiên
Giới thiệu 
Độ bền của cấu trúc vật liệu pha hạt cứng phụ thuộc rất 
nhiều vào sự xuất hiện của các biên bất liên tục. Đối với việc 
xấp xỉ những lời giải không liên tục, phương pháp phần tử 
hữu hạn truyền thống sử dụng không gian xấp xỉ đa thức và 
phụ thuộc rất nhiều vào lưới để đảm bảo độ chính xác của 
các kết quả gần miền suy biến hay những vùng có gradient 
cao. Việc mô phỏng các biên bất liên tục như vết nứt, lỗ 
trống bằng phương pháp phần tử hữu hạn truyền thống đòi 
hỏi mật độ lưới rất lớn. Và việc làm mịn lưới đòi hỏi một 
lượng tài nguyên máy tính khá lớn. Hơn nữa, việc làm mịn 
lưới thường khó có thể tiến hành một cách tự động mà đòi 
hỏi phải có sự can thiệp thủ công của người dùng. Phương 
pháp phần tử hữu hạn mở rộng (extended finite element 
method - XFEM) được giới thiệu bởi Belytschko và Black 
[1], Moës và các cộng sự [2]. Phương này thừa hưởng nền 
tảng lý thuyết vững chắc của phương pháp phần tử hữu hạn 
truyền thống và hạn chế được sự khó khăn trong quá trình 
làm mịn lưới và chia lưới lại. Gần đây nhất, tác giả Bui 
và các cộng sự [3] đã thành công trong việc thiết lập phần 
tử 4 nút nội suy liên tiếp (consecutive-interpolation 4-node 
quadrilateral element - CQ4), dựa trên ý tưởng theo [4, 5]. 
Các hàm cơ bản của CQ4 được xây dựng với hai lần nội 
suy. Lần nội suy thứ nhất hoàn toàn giống với phương pháp 
phần tử hữu hạn tiêu chuẩn. Lần nội suy thứ hai, hàm xấp xỉ 
được nội suy thông qua chuyển vị nút và trung bình đạo hàm 
chuyển vị tại nút. Do đó, hàm dạng nhận được liên tục và có 
đa thức bậc cao hơn mà không làm tăng thêm tổng số bậc tự 
do. Trường ứng suất trở nên liên tục mà không cần những 
biện pháp xử lý phức tạp, kết quả tính toán có sự hội tụ tăng 
đáng kể so với phương pháp phần tử hữu hạn truyền thống.
Bài báo này phát triển phần tử làm giàu cho những 
đối tượng lỗ trống và hạt cứng trong vật liệu dạng hạt dựa 
trên ý tưởng nội suy liên tiếp mở rộng (extended twice-
interpolation finite element method - XTFEM). Điều này 
sẽ tận dụng được những ưu điểm của ý tưởng phương pháp 
phần tử hữu hạn nội suy liên tiếp trong phần tử tứ giác bốn 
nút và ý tưởng làm giàu của XFEM cho những bài toán bất 
liên tục. Các hàm xấp xỉ của phần tử nội suy liên tiếp sẽ 
được mở rộng bằng cách thêm các hàm làm giàu mô tả sự 
bất liên tục của vết nứt, lỗ trống và tạp chất. Kết quả tính 
toán được so sánh với kết quả đã công bố của bài báo khoa 
học quốc tế uy tín.
Cơ sở lý thuyết
Phần tử tứ giác, bốn nút nội suy liên tiếp
Theo [3], quá trình nội suy liên tiếp phần tử tứ giác bốn 
Phân tích động lực học vết nứt trong vật liệu 
lẫn hạt cứng và lỗ rỗng bằng phương pháp 
phần tử hữu hạn nội suy liên tiếp mở rộng
Trương Tích Thiện1*, Trần Kim Bằng1, Phan Ngọc Nhân1, Bùi Quốc Tính2 
1Trường Đại học Bách khoa, Đại học Quốc gia TP Hồ Chí Minh
2Khoa Xây dựng Dân dụng và Môi trường, Viện Công nghệ Tokyo, Nhật Bản
Ngày nhận bài 25/1/2019; ngày chuyển phản biện 31/1/2019; ngày nhận phản biện 25/2/2019; ngày chấp nhận đăng 22/3/2019
Tóm tắt:
Vật liệu có lẫn những hạt cứng là một trong những loại vật liệu được sử dụng phổ biến trong nền công nghiệp hiện 
đại. Vết nứt và khuyết tật xuất hiện sẽ gây ra hiện tượng tập trung ứng suất và làm ảnh hưởng lớn đến độ bền của 
kết cấu. Các khuyết tật trong vật liệu có thể được mô tả dưới dạng các lỗ trống. Ứng xử của vết nứt trong miền 
xuất hiện lỗ trống và các hạt cứng sẽ phức tạp hơn dưới tác dụng của tải trọng động. Trong bài báo này, nhóm tác 
giả phát triển ma trận độ cứng và khối lượng cho các phần tử mô tả vết nứt, lỗ trống và hạt cứng trong vật liệu nền 
bằng phương pháp phần tử hữu hạn nội suy liên tiếp mở rộng (extended twice-interpolation finite element method - 
XTFEM) cho bài toán động lực học, tính toán hệ số cường độ ứng suất động theo thời gian, khảo sát sự ảnh hưởng 
của lỗ trống, hạt cứng gần vết nứt. Các kết quả tính toán hệ số cường độ ứng suất động tại đỉnh vết nứt bằng XTFEM 
sẽ được so sánh với kết quả đã được công bố trên tạp chí khoa học quốc tế uy tín để kiểm chứng độ tin cậy.
Từ khóa: hạt cứng, lỗ trống, mở rộng, nội suy liên tiếp, tải trọng động, vết nứt, XTFEM.
Chỉ số phân loại: 1.9
*Tác giả liên hệ: Email: tttruong@hcmut.edu.vn
2361(8) 8.2019
Khoa học Tự nhiên
nút bao gồm hai giai đoạn:
Giai đoạn nội suy thứ nhất:
Trong phần tử hữu hạn truyền thống, chuyển vị xấp xỉ tại 
điểm x được tính như sau:
[ ] [ ] [ ] [ ]( ) i j k mi j k mu x N u N u N u N u= + + + (1)
1 1
(1 )(1 ), (1 )(1 ),
4 4
1 1
(1 )(1 ), (1 )(1 )
4 4
i j
k m
N r s N r s
N r s N r s
= - - = + -
= + + = - +
(2)
Với u[i], u[j], u[k], u[m] là chuyển vị tại các nút i, j, k, m của 
phần tử và r, s là tọa độ trong hệ tọa độ tự nhiên của phần 
tử tứ giác 4 nút.
Miền hỗ trợ nút i, S
i
 chứa tất cả các phần tử có liên quan 
nút i. Hàm trọng eω được tính như sau: 
' 'i
e
e
e S e
ω
∈
=
∑


 (3)
Với Δ
e
 là diện tích của phần tử e.
Δ
e’
 là diện tích của phần tử e’ trong S
i
.
Hình 1. Điểm cần nội suy và miền hỗ trợ của nút có tọa độ x.
Đạo hàm trung bình tại nút i có thể viết như sau:
[ ][ ] [ ][ ]( )[ ], , ,
1
s
i i
n
i e i ei
x e x e l x l
e l e
u u N uω ω
∈ = ∈
 
= =  
 
∑ ∑ ∑
S S
 (4)
[ ] [ ][ ]( ), ,
i
i i e
l x e l x
e
N Nω
∈
= ∑
S
 (5)
Giai đoạn nội suy thứ hai:
Trong giai đoạn nội suy lần hai, giá trị nội suy trên điểm 
x được tính như sau:
[ ] [ ] [ ] [ ] [ ] [ ]
[ ] [ ] [ ] [ ] [ ] [ ]
, ,
, ,
, ,
, ,
ˆ( )
i j
x x
k m
x x
i i j j
i ix iy y j jx iy y
k k m m
k kx ky y m mx my y
u u u u u u u
u u u u u u
φ φ φ φ φ φ
φ φ φ φ φ φ
= + + + + +
+ + + + + +
x 
(6)
An extended twice-interpolation 
finite element method applied 
to simulate dynamic crack 
behaviour in matrix inclusion 
materials with random holes
Tich Thien Truong1*, Kim Bang Tran1, 
Ngoc Nhan Phan1, Quoc Tinh Bui 2
1VNUHCM - Ho Chi Minh City University of Technology 
2Department of Ci ... ác nhau đó được thể hiện như hình 3. 
 (A) (B) 
Hình 3. (A) Tọa độ địa phương tại đỉnh vết nứt và (B) miền hỗ trợ chứa các phần 
tử được dùng để tính tích phân J trong XFEM và XTFEM 
Mô hình mô phỏng 
Mô hình 1 
Xét bài toán tấm phẳng với chiều rộng W = 20 mm, chiều cao H = 40 mm, vết nứt 
nằm ở giữa tấm có kích thước 2a = 4,8 mm. Lực phân bố đều f0 = 1 N/m ở cạnh trên 
và dưới của tấm. Mô-đun đàn hồi Young: E = 199,992 GPa, hệ số Poisson: v = 0,3, 
khối lượng riêng ρ = 5000 kg/m3. Bước thời gian là Δt = 0,05 μs và tải tác dụng từ 0 
đến 13,62 μs như hình 4. Hạt cứng trong vật liệu nền có bán kính r = 2 mm cách tâm 
vết nứt với khoảng cách d = 6 mm theo phương x. Hạt cứng có mô-đun đàn hồi 
Young: E = 199,992 x 103 GPa, hệ số poisson v của hạt cứng và vật liệu nền đều 
bằng 0,3. Trường hợp đang xét là biến dạng phẳng. Vận tốc lan truyền sóng dọc 
Cd = 7,34 mm μs
-1. Kết quả tính toán từ XTFEM sẽ được so sánh với kết quả tính theo 
XFEM, đã được công bố trong tài liệu [9]. Thời gian sẽ được chuẩn hóa như sau: 
_ 2 /chuan hoa dt C t h với h = H là chiều cao tấm 
Hệ số cường độ ứng suất động Mode I sẽ được chuẩn hóa như sau: 
*
0 0
0
; II
K
K K f a
K
Hình 4. Tấm phẳng với vết nứt ở giữa gầ hạt cứng hình tròn. 
2761(8) 8.2019
Khoa học Tự nhiên
Sự so sánh các kết quả tính toán giữa hai phương pháp 
XFEM [9] và XTFEM được thể hiện trong hình 5. Các kết 
quả hệ số cường độ ứng suất theo thời gian thu được từ 
XTFEM khá tương đồng với các kết quả được tham khảo 
từ [9].
Mô hình 2
Trong ví dụ tiếp theo, xét mô hình tấm phẳng bị nứt có 
chứa đồng thời lỗ trống và hạt cứng. Các thông số về kích 
thước, vật liệu và tải trọng được cho tương tự như mô hình 
1. Tấm chứa 2 lỗ trống với tọa độ tâm và bán kính lần lượt 
là: O
1
(1,25; 3,75) m, r
1 
= 0,6 m và O
2
(3,75; 1,25) m, r
2 
= 
0,7 m. Tấm chứa 2 hạt cứng với tọa độ tâm và bán kính lần 
lượt là: O3(1,25; 1,25) m, r3 = 0,55 m và O4(3,75; 3,75) m, 
r
4 
= 0,8 m. Tấm có vết nứt bắt đầu từ giữa cạnh trái đi về 
tâm của tấm, với chiều dài vết nứt là a = 2,5 m. Lực phân bố 
cạnh trên theo thời gian như hình 6. Cạnh dưới được ngàm. 
Trường hợp đang xét là biến dạng phẳng.
Hình 6. Tấm phẳng chứa vết nứt, lỗ trống và hạt cứng chịu tải 
theo thời gian.
Đồ thị hệ số cường độ ứng suất động Mode I của tấm 
theo thời gian và chuyển vị theo phương x, y, trường ứng 
suất theo hai phương x, y cũng được minh họa qua hình 7 
và 8.
Hình 8. (A) Chuyển vị theo phương x; (B) Chuyển vị theo phương 
y; (C) Ứng suất theo phương x; (D) Ứng suất theo phương y tại 
thời điểm t = 0,15s.
Tại thời điểm t = 0,15s, lúc tải có giá trị lớn nhất, thì ứng 
suất chủ yếu tập trung tại đỉnh vết nứt, còn cạnh vết nứt ứng 
suất gần như bằng 0, từ đó cho thấy trong trường hợp này 
vết nứt chính là biên bất liên tục dạng mạnh nhất. Ứng suất 
cũng tập trung ở hạt cứng phía trên vết nứt, gây ảnh hưởng 
rõ rệt đến hệ số cường độ ứng suất tại đỉnh vết nứt. Với việc 
sử dụng XTFEM, trường ứng suất khá mịn. Xét về chuyển 
vị theo phương y, chuyển vị trên và dưới vết nứt có sự đối 
lập nhau rõ rệt, phần phía trên vết nứt bị kéo ra nên chuyển 
vị rất lớn.
Kết luận
Trong bài báo này, phương pháp phần tử hữu hạn nội suy 
liên tiếp mở rộng (XTFEM) đã được áp dụng để xây dựng 
chương trình mô phỏng bài toán động lực học vết nứt với sự 
ảnh hưởng của các khuyết tật lỗ rỗng và các hạt cứng phân 
bố trong vật thể dựa trên ngôn ngữ lập trình Matlab. Phương 
pháp này có ưu điểm là dễ dàng mô tả được tính chất vật 
lý của vết nứt, lỗ rỗng và hạt cứng thông qua hàm làm giàu 
mà không cần phụ thuộc vào lưới mô hình. Đồng thời có 
thể đem lại một trường ứng suất và biến dạng trơn, liên tục. 
Điều này sẽ ảnh hưởng đến kết quả tính hệ số cường độ ứng 
suất tại đỉnh vết nứt. Một vài ví dụ mô phỏng số đã được 
thực hiện và so sánh để chứng tỏ tính đúng đắn của chương 
trình. Kết quả thu được từ XTFEM khá tương đồng với kết 
Hình 5. So sánh kết quả hệ số cường độ ứng suất động giữa hai phương pháp. 
Sự so sánh các kết quả tính toán giữa hai phương pháp XFEM [9] và XTFEM được 
thể hiện trong hình 5. Các kết quả hệ số cường độ ứng suất theo thời gian thu được từ 
XTFEM k á tương đồng với các kết quả được tham khảo từ [9]. 
Mô hình 2 
Trong ví dụ tiếp theo, xét mô hình tấm phẳng bị nứt có chứa đồng thời lỗ trống và 
hạt cứng. Các thông số về kích thước, vật liệu và tải trọng được cho tương tự như mô 
hình 1. Tấm chứa 2 lỗ trống với tọa độ tâm và bán kính lần lượt là: O1(1,25; 3,75) m, 
r1 = 0,6 m và O2(3,75; 1,25) m, r2 = 0,7 m. Tấm chứa 2 hạt cứng với tọa độ tâm và bán 
kính lần lượt là: O3(1,25; 1,25) m, r3 = 0,55 m và O4(3,75; 3,75) m, r4 = 0,8 m. Tấm có 
vết nứt bắt đầu từ giữa cạnh trái đi về tâm của tấm, với chiều dài vết nứt là a = 2,5 m. 
Lực phân bố cạnh trên theo thời gian như hình 6. Cạnh dưới được ngàm. Trường hợp 
đang xét là biến dạng phẳng. 
Hình 6. Tấm phẳng chứa vết nứt, lỗ trống và hạt cứng chịu tải theo thời gian. 
Đồ thị hệ số cường độ ứng suất động Mode I của tấm theo thời gian và chuyển vị 
theo phương x, y, trường ứng suất theo hai phương x, y cũng được minh họa qua hình 7 
và 8. 
Hình 7. Hệ số cường độ ứng suất động mode I của tấm theo thời gian. 
Hình 5. So sánh kết quả hệ số cường độ ứng suất động giữa hai phương pháp. 
Sự so sánh các kết quả tính toán giữa hai phương pháp XFEM [9] và XTFEM được 
thể hiện trong hình 5. Các kết quả hệ số cường độ ứng suất theo thời gian thu được từ 
XTFEM khá tương đồng với các kết quả được tham khảo từ [9]. 
Mô hình 2 
Trong ví dụ tiếp theo, xét mô hình tấm phẳng bị nứt có chứa đồng thời lỗ trống và 
hạt cứng. Các thông số về kích thước, vật liệu và tải trọng được cho tương tự như mô 
hình 1. Tấm chứa 2 lỗ trống với tọa độ tâm và bán kính lần lượt là: O1(1,25; 3,75) m, 
r1 = 0,6 m và O2(3,75; 1,25) m, r2 = 0,7 m. Tấm chứa 2 hạt cứng với tọa độ tâm và bán 
kính lần lượt là: O3(1,25; 1,25) m, r3 = 0,55 m và O4(3,75; 3,75) m, r4 = 0,8 m. Tấm có 
vết nứt bắt đầu từ giữa cạnh trái đi về tâm của tấm, với chiều dài vết nứt là a = 2,5 m. 
Lực phân bố cạnh trên theo thời gian như hình 6. Cạnh dưới được ngàm. Trường hợp 
đa g xé là biến dạng phẳn . 
Hình 6. Tấm phẳng chứa vết nứt, lỗ trống và hạt cứng chịu tải theo thời gian. 
Đồ thị hệ số cường độ ứng suất động Mode I của tấm theo thời gian và chuyển vị 
theo phương x, y, trường ứng suất theo hai phương x, y cũng được minh họa qua hình 7 
và 8. 
Hình 7. Hệ số cường độ ứng suất động mode I của tấm theo thời gian. 
 (A) (B) (C) (D) 
Hình 8. (A) Chuyển vị theo phương x; (B) Chuyển vị theo phương y; (C) Ứng suất 
theo phương x; (D) Ứng suất theo phương y tại thời điểm t = 0,15s. 
Tại thời điểm t = 0,15s, lúc tải có giá trị lớn nhất, thì ứng suất chủ yếu tập trung tại 
đỉnh vết nứt, còn cạnh vết nứt ứng suất gần như bằng 0, từ đó cho thấy trong trường 
hợp này vết nứt chính là biên bất liên tục dạng mạnh nhất. Ứng suất cũng tập trung ở 
hạt cứng phía trên vết nứt, gây ảnh hưởng rõ rệt đến hệ số cường độ ứng suất tại đỉnh 
vết nứt. Với việc sử dụng XTFEM, trường ứng suất khá mịn. Xét về chuyển vị theo 
phương y, chuyển vị trên và dưới vết nứt có sự đối lập nhau rõ rệt, phần phía trên vết 
nứt bị kéo ra nên chuyển vị rất lớn. 
Kết luận 
Trong bài báo này, phương pháp phần tử hữu hạn nội suy liên tiếp mở rộng 
(XTFEM) đã được áp dụng để xây dựng chương trình mô phỏng bài toán động lực học 
vết nứt với sự ả h hưởng của các khuyết tật lỗ rỗng và các hạt cứng phân bố trong vật 
thể dựa trê ngôn n ữ lập trình Matlab. Phương pháp này có ưu điểm là dễ dàng mô tả 
được tính chất vật lý của vết nứt, lỗ rỗng và hạt cứng thông qua hàm làm giàu mà 
không cần phụ thuộc vào lưới mô hình. Đồng thời có thể đem lại một trường ứng suất 
và biến dạng trơn, liên tục. Điều này sẽ ảnh hưởng đến kết quả tính hệ số cường độ 
ứng suất tại đỉnh vết nứt. Một vài ví dụ mô phỏng số đã được thực hiện và so sánh để 
chứng tỏ tính đúng đắn của chương trình. Kết quả thu được từ XTFEM khá tương 
đồng với kết quả tham khảo từ tài liệu uy tín [9]. Đề tài này có thể phát triển theo 
nhiều hướng liên quan đến bài toán nhiều biên bất liên tục với các loại vật liệu khác. 
LỜI CẢM ƠN 
 ghiên cứu được tài trợ bởi Trường Đại học Bách khoa - Đại học Quốc gia TP Hồ 
Chí Minh trong khuôn khổ đề tài mã số To-KHUD-2017-04. Các tác giả xin trân trọng 
cảm ơn. 
TÀI LIỆU THAM KHẢO 
[1] T. Belytschko, T. Black (1999), “Elastic crack growth in finite elements with 
minimal remeshing”, Int. J. Numer. Meth. Eng., 45, pp.601-620. 
[2] N. Moës, J. Dolbow, T. Belytschko (1999), “A finite element method for crack 
growth without remeshing”, Int. J. Numer. Meth. Eng., 46, pp.131-150. 
[3] Q.T. Bui, Q.D. Vo, Ch. Zhang, D.D. Nguyen (2014), “A consecutive-
interpolation quadrilateral element (CQ4)”, Formulation and Applications Finite Elem. 
Anal. Des., 84, pp.14-31. 
[4] C. Zheng, S.C. Wu, X.H. Tang, J.H. Zhang (2010), “A novel twice-interpolation 
finite element method for solid mechanics problems”, Acta Mech. Sin., 26, pp.265-
278. 
[5] S.C. Wu, W.H. Zhang, X. Peng, B.R. Miao (2012), “A twice-interpolation finite 
element method (TFEM) for crack propagation problems”, Int. J. Comput. Methods, 9, 
pp.12-55. 
 (A) (B) (C) (D) 
Hình 8. (A) Chuyển vị theo phương x; (B) Chuyển vị theo phương y; (C) Ứng suất 
theo phương x; (D) Ứng suất theo phương y tại thời điểm t = 0,15s. 
Tại thời điểm t = 0,15s, lúc tải có giá trị lớn n ất, thì ứng suất chủ yếu tập trung tại 
đỉnh vết nứt, còn cạnh vết nứt ứng suất gần như bằng 0, từ đó cho thấy trong trường 
hợp này vết nứt chính là biên bất liên tục dạng mạnh nhất. Ứng suất cũng tập trung ở 
hạt cứng phía trê vết nứt, gây ảnh hưởng rõ rệt đế hệ số cường độ ứng suất tại đỉnh 
vết nứt. Với việc sử dụng XTFEM, trường ứng suất khá mịn. Xét về chuyể vị theo 
phươn y, c uy vị trên và dưới vết n t có sự đối lập nhau rõ rệt, phần phía trên vết 
nứt bị kéo ra nên chuyển vị rất lớn. 
Kết luận 
Trong bài báo này, phương pháp phần tử hữu hạn nội suy liên tiếp mở rộng 
(XTFEM) đã được áp dụng để xây dựng chương trình mô phỏng bài toán động lực học 
vết nứt với sự ả h hưở của ác khuyết tật lỗ rỗng và các hạt cứng phân bố trong vật 
thể dựa trên ngôn ngữ lập trình Matlab. Phương pháp này có ưu điểm là dễ dàng mô tả 
được tính chất vật lý của vết nứt, lỗ rỗng và hạt ng thông qua hàm làm giàu mà 
không cần phụ thuộc vào lưới mô hình. Đồng thời c thể đem lại một trường ứng suất 
và biến dạng trơn, liên tục. Điều này sẽ ảnh hưởng đến kết quả tí h hệ số cường độ 
ứng suất tại đỉnh vết nứt. Một vài ví dụ mô phỏng số đã được thực hiện và so sánh để 
chứng tỏ tính đúng đắn của chương trình. Kết quả thu được từ XTFEM khá tương 
đồng với kết quả tham khảo từ tài liệu uy tín [9]. Đề tài này có thể phát triển theo 
nhiều hướng liên quan đến bài toán nhiều biên bất liên tục với các loại vật liệu khá . 
LỜI CẢM ƠN 
 ghiên cứu được tài trợ bởi Trường Đại học Bách khoa - Đại học Quốc gia TP Hồ 
Chí Minh trong khuôn khổ đề tài mã số To-KHUD-2017-04. Các tác giả xin trân trọng 
cảm ơn. 
TÀI LIỆU THAM KHẢO 
[1] T. Belytschko, T. Black (1999), “Elastic crack growt in inite elements with 
minimal remeshing”, Int. J. Numer. Meth. Eng., 45, pp.601-620. 
[2] N. Moës, J. Dolbow, T. Belytschko (1999), “A finite element method for crack 
growth without remeshing”, Int. J. Numer. Meth. Eng., 46, pp.131-150.
[3] Q.T. Bui, Q.D. Vo, Ch. Zhang, D.D. Nguyen (2014), “A consecutive-
interpolation quadrilateral element (CQ4)”, Formulation and Applications Finite Elem. 
Anal. Des., 84, pp.14-31. 
[4] C. Zheng, S.C. Wu, X.H. Tang, J.H. Zhang (2010), “A novel twice-interpolation 
finite element method for solid mechanics problems”, Acta Mech. Sin., 26, pp.265-
278. 
[5] S.C. Wu, W.H. Zhang, X. Peng, B.R. Miao (2012), “A twice-interpolation finite 
element method (TFEM) for crack propagation problems”, Int. J. Comput. Methods, 9, 
pp.12-55. 
Hình 7. Hệ số cường độ ứng suất động mode I của tấm theo thời 
gian.
 (A) (B)
 (C) (D)
2861(8) 8.2019
Khoa học Tự nhiên
quả tham khảo từ tài liệu uy tín [9]. Đề tài này có thể phát 
triển theo nhiều hướng liên quan đến bài toán nhiều biên bất 
liên tục với các loại vật liệu khác.
LỜI CẢM ƠN
Nghiên cứu được tài trợ bởi Trường Đại học Bách khoa 
- Đại học Quốc gia TP Hồ Chí Minh trong khuôn khổ đề tài 
mã số To-KHUD-2017-04. Các tác giả xin trân trọng cảm 
ơn.
TÀI LIỆU THAM KHẢO
[1] T. Belytschko, T. Black (1999), “Elastic crack growth in finite 
elements with minimal remeshing”, Int. J. Numer. Meth. Eng., 45, 
pp.601-620.
[2] N. Moës, J. Dolbow, T. Belytschko (1999), “A finite element 
method for crack growth without remeshing”, Int. J. Numer. Meth. 
Eng., 46, pp.131-150.
[3] Q.T. Bui, Q.D. Vo, Ch. Zhang, D.D. Nguyen (2014), “A 
consecutive-interpolation quadrilateral element (CQ4)”, Formulation 
and Applications Finite Elem. Anal. Des., 84, pp.14-31. 
[4] C. Zheng, S.C. Wu, X.H. Tang, J.H. Zhang (2010), “A 
novel twice-interpolation finite element method for solid mechanics 
problems”, Acta Mech. Sin., 26, pp.265-278.
[5] S.C. Wu, W.H. Zhang, X. Peng, B.R. Miao (2012), “A twice-
interpolation finite element method (TFEM) for crack propagation 
problems”, Int. J. Comput. Methods, 9, pp.12-55.
[6] Zuoyi Kang, Tinh Quoc Bui, Du Dinh Nguyen, Takahiro 
Saitoh, Sohichi Hirose (2015), “An extended consecutive-
interpolation quadrilateral element (XCQ4) applied to linear elastic 
fracture mechanics”, Acta Mech. Sin., 80, pp.17-55.
[7] S. Mohammadi (2012), XFEM fracture analysis of composites, 
John Wiley & Sons.
[8] D. Motamedi and S. Mohammadi (2010), “Dynamic crack 
propagation analysis of orthotropic media by the XFEM”, Int. J. 
Fract., 161, pp.21-39.
[9] S. Jiang, C. Du, C. Gu and X. Chen (2014), “XFEM analysis 
of the effects of voids, inclusions and other cracks on the dynamic 
stress intensity factor of a major crack”, Fatigue Fract. Eng. Mater. 
Struct., 37, pp.1-17.