Tính toán và chẩn đoán sự phát triển vết nứt của tấm bằng Fem - Wavelet

 TẠP CHÍ KHOA HỌC & CÔNG NGHỆ CÁC TRƯỜNG ĐẠI HỌC KỸ THUẬT  SỐ 74 - 2009 
78 
TÍNH TOÁN VÀ CHẨN ĐOÁN SỰ PHÁT TRIỂN VẾT NỨT 
CỦA TẤM BẰNG FEM - WAVELET 
COMPUTATION AND DIAGNOSTIC OF THE DEVELOPMENT 
OF THE CRACK IN THE PLATE BY THE FEM - WAVELET 
Nguyễn Hoài Sơn, Lâm Phát Thuận 
Trường Đại học Sư Phạm Kỹ Thuật Thành Phố Hồ Chí Minh 
TÓM TẮT 
Phát hiện sự tồn tại vết nứt và tính toán trường ứng suất kỳ dị xung quanh đỉnh vết nứt là một 
vấn đề rất quan trọng trong việc chẩn đoán và đánh giá khả năng làm việc của chi tiết máy. Có nhiều 
kỹ thuật đã được nghiên cứu và sử dụng để giải quyết vấn đề trên. Trong khuôn khổ cơ học rạn nứt 
đàn hồi tuyến tính, phương pháp phần tử hữu hạn (FEM) với việc sử dụng phần tử suy biến Barsoum 
ở đỉnh vết nứt và phép biến đổi Wavelet đã cho thấy những hiệu quả rất thiết thực. Ở bài báo này, vết 
nứt được đánh giá thông qua hệ số cường độ ứng suất (Stress Intensity Factor - SIF) được xác định 
bằng kỹ thuật tương quan chuyển vị (Displacement Correlation Technique – DCT) và trường ứng suất 
được tính toán thông qua phép ngoại suy từ kết quả phân tích FEM. Phép biến đổi Wavelet thực hiện 
cho cả trường hợp biến đổi liên tục và biến đổi rời rạc đều đạt được kết quả tốt. 
Keywords: FEM, Wavelet, Cơ học rạn nứt, Phần tử Barsoum, Hệ số cường độ ứng suất 
ABSTRACT 
Identifying cracks and computing singularity stress field surrounding the top of the cracks is a 
very important task in diagnosis and evaluation of workability of the machine details. Many techniques 
have been investigated and employed to solve the above-mentioned problem. Upon the assumptions 
of the linear elastic fracture mechanics, the Finite Element Method (FEM) together with Barsoum 
singularity element at tips of the crack and Wavelet transform reveal effectiveness in terms of 
methodology. In this paper, the cracks were evaluated through Stress Intensity Factor (SIF) 
determined by Displacement Correlation Technique (DCT) and stress field computed through 
extrapolation methods from FEM analysis. Wavelet transform carried out in the two cases - continuous 
and discrete transforms offer good results. 
Keywords FEM, Wavelet, Fracture mechanics, Barsoum element, Stress Intensity Factor 
I. GIỚI THIỆU 
Việc ứng dụng các phương pháp số tỏ ra 
rất hiệu quả trong việc giải quyết bài toán nứt. 
Và FEM là một công cụ mạnh mẽ đã được sử 
dụng từ rất sớm. Cùng với việc Barsoum đưa ra 
phần tử suy biến điểm 1/4 và Irwin với khái 
niệm hệ số cường độ ứng suất, việc đánh giá độ 
bền cũng như tính toán trường ứng suất của chi 
tiết có vết nứt có thể được xác định từ kết quả 
FEM thông qua kỹ thuật tương quan chuyển vị 
và phép ngoại suy. Phần tử Barsoum không 
những thể hiện rất tốt trường ứng suất kỳ dị gần 
đỉnh vết nứt mà còn giảm đáng kế khả năng thất 
thoát năng lượng gần đỉnh vết nứt. 
Với phép biến đổi Wavelet, sự tương 
quan giữa tín hiệu phân tích - trong bài báo này 
là trường chuyển vị của tấm có vết nứt - và hàm 
Wavelet được biểu thị thông qua độ lớn của các 
hệ số C(a,b) và hàm chi tiết Dj(k). Tại vị trí vết 
nứt, hàm chuyển vị bị mất liên tục, thể hiện tính 
tương quan cao với hàm Wavelet, đồ thị hệ số 
C(a,b) và hàm chi tiết Dj(k ) sẽ có giá trị rất lớn. 
Đây là cơ sở cho sự phát hiện vết nứt bằng 
Wavelet. 
II. HỆ SỐ CƯỜNG ĐỘ ỨNG SUẤT (SIF) 
Hệ số cường độ ứng suất được xem như 
thướt đo của trường ứng suất và biến dạng, là 
đại lượng đánh giá độ bền của tấm có vết nứt. 
Có nhiều phương pháp tính SIF từ kết quả FEM 
trong đó kỹ thuật tương quan chuyển vị là một 
trong những kỹ thuật đơn giãn nhất và được sử 
dụng sớm nhất để xác định SIF [1]. DCT sử 
dụng 4 nút gần kề đỉnh vết nứt để tính SIF, 
Hình 1. 
 TẠP CHÍ KHOA HỌC & CÔNG NGHỆ CÁC TRƯỜNG ĐẠI HỌC KỸ THUẬT  SỐ 74 - 2009 
79 
Hình 1. Phần tử tam giác điểm 1/4 tại đỉnh vết 
nứt 
)]()(4[
2
1
5342 vvvv
L
G
K I 




 (1) 
Với, 






1
3 : đối với bài toán ứng suất phẳng 
iv : chuyển vị nút thứ i của phần tử tam 
giác điểm 1/4 lân cận đỉnh vết nứt (i=2,3,4,5) 
G : Mô đun đàn hồi trượt 
L : kích thước cạnh l123, l145 của phần tử 
KI : hệ số cường độ ứng suất ở mode I. 
Với bài toán tấm nứt thuần túy theo mode 
I,KI có thể được tính theo công thức lý thuyết 
như sau: 
)/(. baFaK I  (2) 
Trong đó, F(a/b) : hệ số hiệu chỉnh so 
với trường hợp tấm rộng vô hạn 
a : Chiều dài vết nứt 
b : Chiều rộng mô hình 
III. TÍNH TOÁN TRƯỜNG ỨNG SUẤT 
DỰA TRÊN PHÉP NGOẠI SUY 
Độ chính xác của ứng suất và biến dạng 
trên phần tử phụ thuộc vào vị trí điểm tính toán. 
Đối với phần tử tứ giác 8 nút, ứng suất và biến 
dạng có giá trị chính xác nhất tại điểm bốn 
điểm tích phân Gauss 2×2, hình 2, tọa độ tự 
nhiên tương ứng: 
    3/1,3/1,  
Như vậy, giá trị ứng suất và biến dạng tại 
một nút được tính toán trên các phần tử khác 
nhau chứa nút đó là khác nhau, điều này dẫn 
đến trường ứng suất và biến dạng bị mất liên 
tục. Để tạo ra trường ứng suất và biến dạng liên 
tục gồm hai bước: 
Hình 2. Điểm tích phân 2x2 và các nút góc trên 
phần tử 
Bước 1: Từ các giá trị được tính toán tại các 
điểm tích phân, sử dụng phương pháp ngoại suy 
để tìm giá trị tại các nút ở bốn đỉnh phần tử. 
Đối với, phần tử đẳng tham số bậc hai 8 
nút, công thức ngoại suy được cho như sau 
    
)()( mm
fLf  (3) 
Trong đó,  f là giá trị trên nút ở bốn góc 
phần tử cần tìm; 
)(mf lá giá trị tại điểm tích phân 2×2;  )( mL 
là ma trận ngoại suy, 
 



















2/312/12/312/1
2/12/312/12/31
2/312/12/312/1
2/12/312/12/31
)( m
L
 (4) 
Giá trị các nút giữa bằng trung bình cộng 
của hai nút trên cạnh tương ứng. 
Bước 2: Trung bình cộng các giá trị được tính 
trên các phần tử liên quan. Sau cùng trường ứng 
suất và biến dạng sẽ được nội suy theo hàm 
dạng. 
IV. PHÉP BIẾN ĐỔI WAVELET 
4.1 Biến đổi Wavelet liên tục 
Biến đổi wavelet liên tục (Continuous 
Wavelet Transform- CWT) là tổng trên toàn 
khoảng thời gian của tín hiệu nhân theo tỷ lệ 
(scale) và dịch vị (position) của hàm wavelet 
mẹ ( )x [2]. 
,
1
( , ) ( )
( ). ( )a b
x b
C a b f x dx
aa
f x x dx






 
  
 



 (5) 
 TẠP CHÍ KHOA HỌC & CÔNG NGHỆ CÁC TRƯỜNG ĐẠI HỌC KỸ THUẬT  SỐ 74 - 2009 
80 
Trong đó: a là tỷ lệ (scale) và b là dịch vị 
(position) là những số thực và a ≠ 0 
Kết quả của CWT là tạo ra các hệ số 
wavelet C(a,b) thể hiện sự tương quan của hàm 
wavelet với tín hiệu phân tích )(xf . 
Nghịch đảo CWT cho phép thu lại tín 
hiệu )(xf từ các hệ số C(a,b), được định nghĩa 
 






a b
ba
a
da
dbdxxbaC
K
xf
2,
)(),(
1
)(  (6) 
Trong đó, hằng số 

K phụ thuộc vào loại 
wavelet. 
4.2 Biến đổi Wavelet rời rạc 
Khuyết điểm của CWT là tạo ra rất nhiều 
hệ số C(a,b). Vì vậy, để rút ngắn thời gian và 
đơn giãn hóa quá trình tính toán nhưng vẫn đảm 
bảo độ chính xác cần thiết. Người ta thay việc 
sử dụng các hệ số tỷ lệ a và dịch vị b liên tục 
bằng các hệ số rời rạc, lũy thừa cơ số 2 là 
),(,2,2 Zjkkba jj  , còn gọi là mức 
dyadic [3]. Phép biến đổi như vậy gọi là phép 
biến đổi wavelet rời rạc (Discrete Wavelet 
Transform - DWT). 
/2
,
,
2 ( ). (2 )
( ). ( )
j j
j k
j k
C f x x k dx
f x x dx



 



 



 (7) 
Tín hiệu )(xf cũng có thể tái tạo từ hệ 
số 
kj
C
,
 và thuật toán này gọi là phép biến đổi 
wavelet rời rạc ngược (IDWT). Khi đó, công 
thức thu hồi tín hiệu )(xf có thể chia thành 2 
thành phần: 
 0
0
, 20
0 ,
0
1
( ) ( , ) ( )
1
( , ) ( )
a
a b
b
a b
b
da
f x C a b x dxdb
K a
D a b x db
K a








 

 

 (8) 
Từ phương trình (8) áp dụng cho trường 
hợp hệ số tỉ lệ a và dịch vị b ở mức dyadic J, ta 
thu được một dãy các hệ số 



 dxxxfkcD
kJJ
)().()(
,
 : hệ số chi tiết 
mức J (9) 
Áp dụng phương trình (9) với )(x ta thu 
được một dãy các hệ số khác 



 dxxxfkcA
kJJ
)().()(
,
 : hệ số xấp 
xỉ mức J (10) 
Từ bộ hệ số )(kcD
J
và )(kcA
J
ta có thể 
thu hồi lại tín hiệu f(x) như sau: 
 










k
kjJ
J
j k
kjJ
xkcAxkcDxf )()()()()(
,,

 (11) 
Trong đó, 
,( ) ( ) ( )j J j k
k
D x cD k x


  : hàm chi tiết ở 
mức J 




k
kjJJ
xkcAxA )()()(
,
 : hàm xấp xỉ ở mức J 
Như vậy, tín hiệu )(xf có thể biểu diễn 
như là tổng xấp xỉ của nó ở mức J cộng với 
tổng tất cả các chi tiết của nó nhỏ hơn hoặc 
bằng mức J. 



Jj
jJ
xDxAxf
)()()( (12) 
Như vậy, quá trình phân tích wavelet liên 
tục, kết quả phân tích sẽ tạo ta hệ số wavelet 
C(a,b). Giá trị của hệ số C(a,b) càng lớn thì tính 
tương quan giữa tín hiệu f(x) và hàm wavelet 
)(x
ab
 càng cao. Đối với phân tích rời rạc, kết 
quả sẽ tạo ra hàm xấp xỉ )(xA
J
ở mức J tương 
ứng với tín hiệu có tần số thấp và hàm chi tiết 
)(xD
j
 tương ứng với tín hiệu ở tần số cao. Quá 
trình này cũng tạo ra hệ số xấp xỉ )(kcA
J
 và hệ 
số chi tiết )(xcD
j
. Giá trị của các hệ số này 
cũng thể hiện tính tương quan giữa tín hiệu f(x) 
với hàm wavelet )(x
ab
 và hàm tỷ lệ )(
,
x
kj
 . 
V. ỨNG DỤNG 
5.1 Mô hình phân tích [4,5] 
Mô hình tính toán là tấm nứt cạnh, một 
đầu ngàm, một đầu chịu kéo, Hình 3. 
 TẠP CHÍ KHOA HỌC & CÔNG NGHỆ CÁC TRƯỜNG ĐẠI HỌC KỸ THUẬT  SỐ 74 - 2009 
81 
 Hình 3. Mô hình bài toán 
Với các thông số đặc trưng như sau: 
Chiều cao tấm: H=210cm 
Chiều rộng tấm: W=80cm 
Chiều dày tấm: t=1cm 
Module đàn hồi: E=2e7 N/cm2 
Hệ số poisson:  = 0.3 
Chiều dài vết nứt: a=20(cm) 
Tọa độ đỉnh vết nứt: ynut=105 
Tải trọng phân bố điều: q=50(N/cm) 
Độ bền vật liệu: KIc=140000 N/cm
3/2
5.2 Lưới PTHH 
Hình 4. Mô hình chia lưới phần tử hữu hạn sử dụng phần tử Barsoum điểm 1/4 ở đỉnh vết nứt 
5.3 Kết quả 
Hình 5. Bảng màu trường chuyển vị của tấm 
 TẠP CHÍ KHOA HỌC & CÔNG NGHỆ CÁC TRƯỜNG ĐẠI HỌC KỸ THUẬT  SỐ 74 - 2009 
82 
Bảng 1. Kết quả tính toán KI trong các trường hợp chia lưới cho bài toán nứt cạnh 
Lưới 
Số phần 
tử 
Bậc tự do Le/a KFEM Klt 
Sai số 
tương đối 
(%) 
1 164 1090 0.0833 750.6082 756.5252 0.7821 
2 332 2154 0.0625 751.6175 756.5252 0.6487 
3 440 2826 0.0625 750.3358 756.5252 0.8181 
4 560 3578 0.0500 751.9148 756.5252 0.6094 
5 620 3954 0.0500 751.8286 756.5252 0.6208 
6 700 4442 0.0500 751.4653 756.5252 0.6688 
7 740 4690 0.0417 751.8028 756.5252 0.6242 
8 818 5170 0.0417 751.5941 756.5252 0.6518 
9 856 5406 0.0417 751.5841 756.5252 0.6518 
10 900 5674 0.0417 750.9568 756.5252 0.7361 
Nhận xét 
Từ kết quả thu được ở Bảng 1 cho thấy, sai số giữa SIF tính bằng FEM và SIF tính bằng lí 
thuyết là bé. Chứng tỏ, trường chuyển vị tính toán từ FEM đạt độ chính xác cao, giải thuật FEM là hợp 
lí. 
Từ kết quả trường chuyển vị, ta có thể xác định được trường ứng suất, Hình 6. 
Hình 6. Bảng màu trường ứng suất của tấm 
với trường chuyển vị thu được, thực hiện phép biến đổi Wavelet với tín hiệu phân tích là giá trị chuyển 
vị ở tọa độ x = 0. Ta có thể nhận thấy sự tồn tại vết nứt, thông qua sự nhiễu loạn trên đồ thị hệ số 
C(a,b) trong phân tích liên tục (sử dụng họ Wavelet Daubechies bậc 5) và đồ thị hàm chi tiết Dj(k) 
trong phân tích rời rạc (sử dụng họ Wavelet Biorthogonal bậc 6.8), Hình 7. 
 TẠP CHÍ KHOA HỌC & CÔNG NGHỆ CÁC TRƯỜNG ĐẠI HỌC KỸ THUẬT  SỐ 74 - 2009 
83 
Hình 7. Sự nhiễu loạn trên đồ thị hệ số C(a,b) và Hàm chi tiết Dj(k) 
VI. KẾT LUẬN 
Kết quả trường chuyển vị và trường ứng 
suất tính toán được bằng FEM dưới những điều 
kiện của cơ học rạn nứt đàn hồi tuyến tính đã 
thể hiện được tính chất suy biến tại đáy vết nứt. 
Điều đó cho thấy, phương pháp phần tử hữu 
hạn kết hợp với việc sử dụng phần tử Barsoum 
mang lại nhiều hiệu quả cao trong giải quyết bài 
toán nứt. 
Hệ số cường độ ứng suất tính toán từ kết 
quả FEM thông qua kỹ thuật tương quan 
chuyển vị có sai số rất nhỏ so với giá trị SIF 
tính theo lí thuyết, chứng tỏ giải thuật FEM là 
hợp lí và có khả năng ứng dụng cao trong bài 
toán rạn nứt. 
Phép biến đổi Wavelet đạt hiệu quả cao 
trong việc phát hiện sự tồn tại vết nứt. Dựa trên 
sự nhiễu loạn của đồ thị hệ số C(a,b) và hàm chi 
tiết Dj(k), chúng ta có thể nhận thấy được sự 
mất liên tục trong trường chuyển vị của tấm tại 
vị trí có vết nứt. 
TÀI LIỆU THAM KHẢO 
1. Barry D. Fehl, Kevin Z. Truman; An evaluation of fracture mechanics quarter-point displacement 
techniques used for computing stress intensity factors; Engineering Structures, Vol 21, pp.406-
415 
2. G. P. Nikishkov; Introduction To The Finite Element Method; Lecture Notes. University of Aizu, 
Aizu-Wakamatsu 965-8580, Japan - 2007 
3. A.V. Ovanesova, L.E. Suarez; Applications of Wavelet transform to damage detection in frame 
structures; Engineering Structures, Vol 26, pp.39-49, (2004). 
4. Nguyễn Hoài Sơn, Lê Thanh Phong, Mai Đức Đãi; Ứng dụng phương pháp phần tử hữu hạn trong 
tính toán kỹ thuật; NXB Đại học Quốc gia Tp. HCM - 2008. 
5. Nguyễn Quận; Nghiên cứu sự phát triển vết nứt trên chi tiết cơ khí bằng FEM - Wavelet; Luận 
văn thạc sĩ - Đại học Sư Phạm Kỹ Thuật TP.HCM - 2008. 
Địa chỉ liên hệ: Nguyễn Hoài Sơn – Tel: 0903.682.514; Email: son55vn@yahoo.com 
 Trường Đại học Sư Phạm Kỹ Thuật Thành Phố Hồ Chí Minh 
 Số 1, Võ Văn Ngân, Q. Thủ Đức, Tp. Hồ Chí Minh