CHƯƠNG 1
1.1
GIỚI THIỆU
Bài toán truyền nhiệt và tương tác cơ-nhiệt
Sự thay đổi nhiệt độ gây ra hiện tượng giãn nở nhiệt trong vật liệu và nếu không
được quan tâm đúng mức thì có thể gây ra ứng suất nhiệt lớn trong kết cấu và chi
tiết máy. Thêm vào đó, thuộc tính vật liệu có thể thay đổi theo nhiệt độ. Vì thế,
nghiên cứu về ứng xử cơ-nhiệt của vật rắn là một chủ đề quan trọng phục vụ các
các ngành công nghiệp.
1.2
Phương pháp phần tử hữu hạn (FEM) và hạn chế của nó
Phương pháp phần tử hữu hạn (FEM) hiện vẫn đang là phương pháp số phổ biến
nhất được dùng để giải các bài toán kỹ thuật. Tuy vậy, FEM vẫn chứa đựng nhiều
hạn chế như đã được chỉ ra bởi [1, 2]: i) lưới phần tử cần phải cập nhật nếu miền
bài toán thay đổi; ii) nghiệm thu được bị ảnh hưởng bởi sự méo dạng của phần
tử; iii) trường đạo hàm bất liên tục tại nút (không phản ánh bản chất vật lý); và
iv) hiện tượng “khóa”.
1.3
Xu hướng phát triển các phương pháp số
Hạn chế của FEM là động lực thúc đẩy các nhà nghiên cứu phát triển các phương
pháp số mới, như phương pháp Đẳng hình học (IGA) [3] và lớp các phương pháp
không lưới [4, 5, 6]. Mỗi phương pháp đều có ưu và nhược điểm riêng. Hướng
tiếp cận thứ hai là cải tiến FEM nhằm tận dụng những lợi thế vốn có của nó, đồng
thời khắc phục được những điểm yếu. Mới đây, kỹ thuật nội suy liên tiếp (CIP)
đã được đề xuất cho phần tử tam giác và phần tử tứ giác [7, 8]. So sánh với FEM
thông thường, việc tích hợp CIP mang đến cho phần tử nhiều thuộc tính tốt và do
đó, đây là một kỹ thuật đáng để đầu tư phát triển sâu hơn.
1.4
Đóng góp mới của luận án
Luận án trình bày sự phát triển của nhóm các phần tử mới dựa trên kỹ thuật nội
suy liên tiếp (CIP), trong đó tập trung vào các bài toán phân tích truyền nhiệt và
cơ-nhiệt trong vật rắn nguyên vẹn và vật rắn có vết nứt. Các đóng góp mới của
luận án như sau:
1
a) Mở rộng kỹ thuật CIP trong phần tử hai chiều để khảo sát vật rắn có vật
liệu đẳng hướng chịu tác động bởi tải cơ và tải nhiệt, với giả thiết cơ học
nứt nhiệt-đàn hồi tuyến tính.
b) Tiếp tục mở rộng các kết quả ở mục a) cho vật liệu trực hướng.
c) Phát triển kỹ thuật CIP cho phần tử ba chiều.
d) Giới thiệu một công thức tổng quát cho kỹ thuật CIP.
1.5
Ý nghĩa khoa học và thực tiễn của các đóng góp bởi luận án
- Đề xuất phương thức tính toán mới cho bài toán nứt cơ-nhiệt đàn hồi tuyến tính
hai chiều, dựa trên sự phát triển sâu phần tử XCQ4, ứng dụng trong vật liệu đẳng
hướng và vật liệu trực hướng.
- Đề xuất công thức tính toán tổng quát cho kỹ thuật CIP, áp dụng cho nhiều loại
phần tử từ 1D tới 3D.
1.6
Phương pháp nghiên cứu
Luận án tập trung vào các kỹ thuật cải tiến dựa trên CIP đối với phương pháp
Phần tử hữu hạn (FEM) và phương pháp Phần tử hữu hạn mở rộng (XFEM) cho
vật rắn không có vết nứt và vật rắn có vết nứt dưới tác dụng của tải cơ-nhiệt, với
vật liệu trực hướng và vật liệu đẳng hướng.
1.7
Đối tượng và phạm vi nghiên cứu
Phương pháp và kỹ thuật tính toán đề xuất trong luận án hướng tới bài toán cơnhiệt trong vật rắn nguyên vẹn và có vết nứt. Giả thiết đặt ra là ứng xử vật liệu
trong vùng đàn hồi tuyến tính và biến dạng nhỏ. Với vật rắn có vết nứt, chỉ xét
bài toán hai chiều. Với vật rắn không có vết nứt, biểu thức tính toán được tổng
quát hóa và áp dụng được cho bài toán một chiều, hai chiều và ba chiều.
1.8 Bố cục luận án
Luận án được tổ chức như sau:
Chương 1 đánh giá tổng quan về FEM và các nhược điểm của nó trong việc giải
các bài toán vi phân đạo hàm riêng. Chương 2 trình bày ngắn gọn về hệ phương
trình chủ đạo của bài toán cơ-nhiệt đàn hồi tuyến tính. Chương 3 mô tả chi tiết
2
kỹ thuật nội suy liên tiếp (CIP) cho bài toán một chiều và hai chiều. Chương 4
dành riêng cho đóng góp mới trong phân tích bài toán cơ học nứt nhiệt-đàn hồi
tuyến tính hai chiều với vật liệu đẳng hướng, sử dụng kỹ thuật CIP và phương
pháp hàm làm giàu. Chương 5 trình bày sự mở rộng của quy trình tính toán ở
chương 4 cho vật liệu trực hướng, trong đó tính định hướng vật liệu có vai trò
quan trọng. Chương 6 phát triển một công thức tổng quát để tích hợp kỹ thuật
CIP vào đa dạng các phần tử từ một chiều đến ba chiều. Cuối cùng, kết luận và
hướng phát triển sẽ được rút ra trong Chương 7.
CHƯƠNG 2
2.1
BÀI TOÁN CƠ-NHIỆT ĐÀN HỒI TUYẾN TÍNH
Biểu thức
Xét vật thể đẳng hướng Ω với biên Γ. Bài toán cơ-nhiệt đàn hồi tuyến tính là bài
toán cặp đôi của các vấn đề đàn hồi tuyến tính và truyền nhiệt. Dạng yếu Galerkin
của bài toán được viết như sau
cT Td T k Td Q Td
q Td h Ta T Td 0
2
3
u ud ε : C : ε ε d b ud t ud 0 .
T
(2.1)
(2.2)
t
Trong phương trình (2.1), Q kí hiệu cho nguồn thu/phát nhiệt; k là hệ số dẫn
nhiệt; ρ là khối lượng riêng; c là nhiệt dung riêng; h là hệ số đối lưu; Ta là nhiệt
độ môi trường; và T biểu diễn đạo hàm bậc nhất của nhiệt độ theo thời gian. Ở
phương trình (2.2), b là lực khối tác dụng lên vật thể; 𝒖̈ là gia tốc, tức là đạo hàm
bậc hai của chuyển vị u theo thời gian; là tensor biến dạng; và ε T là tensor biến
dạng nhiệt.
CHƯƠNG 3
KỸ THUẬT NỘI SUY LIÊN TIẾP CHO CÁC BÀI
TOÁN 1D VÀ 2D
3
Vấn đề bất liên tục phi vật lý tại nút trong phương pháp Phần tử hữu
hạn (FEM): ví dụ với phần tử thanh 2 nút
Xét miền một chiều (1D) Ω. Một hàm u(x) bất kỳ xác định trong Ω có thể xấp xỉ
3.1
bằng FEM như sau [9, 10]
n
u x u x N I x u I Nu ,
(3.1)
I 1
với NI là hàm dạng Lagrange ứng với nút I (chỉ số toàn cục); uI là giá trị của hàm
u(x) tại nút I và n là tổng số nút. Mỗi phần tử thanh hai nút (L2) có dạng một
đoạn thẳng nối liền hai nút, kí hiệu là nút i và nút j. Hàm dạng Ni và Nj tương ứng
với nút i và j, được viết trong hệ tọa độ tự nhiên như sau
1
1
Ni 1 , và N j 1
2
2
(3.2)
Đạo hàm hàm dạng theo phương x, tính tại điểm x bất kỳ nằm trong phần tử e là
dN d
1 2
1
1
Ni , x x i
, và N j , x x ,
(3.3)
d dx
2 le
le
le
với l e là chiều dài phần tử. Do đó, đạo hàm của hàm u(x) tính tại điểm x trong
phần tử e được cho bởi
1
1
1
u,[xe ] x Ni , x x N j , x x ui u j u j ui ,
le
le
le
(3.4)
Bây giờ, xét miền 1D Ω có chiều dài L = 1 được chia đều bởi hai phần tử L2, như
minh họa ở Hình 3.1. Hai phần tử được kí hiệu là e1 và e2, trong đó các nút được
đánh số toàn cục là 1, 2 và 3. Nút 2 là nút chung của phần tử e1 và phần tử e2, do
đó, đạo hàm của hàm u tại nút 2 được tính theo phương trình (3.4)
1
u,[xe1 ] x2 u2 u1 , (xem nút 2 thuộc phần tử e1)
le
u,[xe2 ] x2
1
u3 u2 , (xem nút 2 thuộc phần tử e2).
le
Hình 3.1. Miền một chiều được chia lưới bởi hai phần tử L2
4
(3.5)
(3.6)
Các giá trị tại nút u1, u2 và u3 nói chung sẽ không giống nhau. Vì vậy, đạo hàm
u,x không liên tục tại nút 2 (không đúng về vật lý). Đây là vấn đề chung của FEM
truyền thống, vì tính chất liên tục C0 của hàm dạng Lagrange. Chi tiết về liên tục
C0 được bàn luận trong [9].
3.2
Kỹ thuật nội suy liên tiếp cho phần tử thanh 2 nút: phần tử CL2
Phần tử CL2 được tạo ra bằng cách áp dụng kỹ thuật nội suy liên tiếp (CIP) vào
phần tử L2. Sử dụng kỹ thuật CIP [7, 8, 11, 12, 13], hàm u(x) ở phương trình
(3.1) được xấp xỉ bởi
u x u x I u[ I ] Ix u, x[ I ] ,
n
(3.7)
I 1
Trong đó u[I] là giá trị của u(x) tại nút I khi tính theo FEM, và u, x[ I ] là đạo hàm
trung bình tại nút I, cụ thể
n
u[ I ] u xI N I xI uI Ν[ I ]u ,
(3.8)
u, x[ I ] we u,[xe ] ( xI ) N[, xI ]u .
(3.9)
I 1
eS I
Ở đây, SI là tập hợp các phần tử có chung nút I. Trọng số we được tính theo kích
thước phần tử. Với phần tử một chiều, we là tỉ số giữa chiều dài le của phần tử e
và tổng chiều dài của tất cả các phần tử trong tập hợp SI.
Thế các phương trình (3.8 – 3.9) và phương trình (3.7), ta thu được
u x I N[ I ] Ix N, x[ I ] u Ru ,
n
(3.10)
I 1
trong đó R là vector các hàm dạng theo CIP. Các hàm bổ sung ϕI và ϕIx tùy thuộc
vào kiểu phần tử [7, 8]. Trong phạm vi phần tử thanh hai nút, các hàm bổ sung
ứng với nút cục bộ i và j được xây dựng trên cơ sở các hàm dạng Lagrange
i N i N i2 N j N 2j N i , ix le Ni2 N j
(3.11)
j N j N 2j N i N i2 N j , jx le N 2j Ni
5
(3.12)
3.2.1
Tính các hàm dạng CIP
CL2
L2
Hình 3.2. Hàm dạng ứng với nút 2
khi tính theo phần tử CL2 và phần
tử L2
Hình 3.3. Đạo hàm bậc nhất của hàm
dạng ứng với nút 2 khi tính theo phần
tử CL2 và phần tử L2
Tiếp tục ví dụ nêu ở Hình 3.1. Hàm dạng CIP ứng với nút 2, kí hiệu R2 được tính
bởi
Figure 3.3. Shape functions
associated with node
2 (global) of the Example
1
2
x
x
N
N
N
In provided
element ein
1: R2 x
(3.13)
ix
j
j
j
i ,
Figure 3.1
le computed by CL2 and L2 elements and (b) Their
first order derivatives
1
In element e2: R2 x i x jx x Ni Ni2 N j
(3.14)
le
Hình 3.2. biểu diễn hàm dạng CIP R2 và hàm dạng Lagrange N2 (FEM truyền
thống). Có thể thấy, hàm dạng CIP là hàm trơn tại nút.
3.2.2 Đạo hàm bậc nhất của hàm dạng CIP
Đạo hàm bậc nhất của hàm R2 được thể hiện ở Hình 3.3, cho thấy một đường
cong trơn. Ngược lại, đạo hàm bậc nhất của hàm N2 (FEM truyền thống) không
trơn và có bước nhảy rất rõ rệt khi đi qua nút 2. Hình ảnh này biểu hiện sự bất
liên tục tại nút của trường đạo hàm trong FEM truyền thống.
3.3
Kỹ thuật nội suy liên tiếp (CIP) cho phần tử tam giác ba nút (CT3) và
phần tử tứ giác bốn nút (CQ4)
3.3.1 Biểu thức xấp xỉ theo CIP cho miền 2 chiều
6
Hàm u(x) bất kỳ xác định trong miền hai chiều (2D) được xấp xỉ theo CIP như
sau [7, 8, 14, 15, 16]
u x u x I N[ I ] Ix N[, xI ] Iy N[, yI ] u Ru .
n
(3.15)
I 1
Phương trình (3.15) là sự mở rộng của phương trình Equation (3.10) với sự xuất
hiện của thành phần đạo hàm của u(x) theo phương y.
3.4
Kết luận
Để đảm bảo tính đầy đủ của luận án, biểu thức xấp xỉ theo CIP cho các phần tử
1D và 2D được trình bày trong chương này. Thuộc tính tốt của các phần tử tích
hợp CIP đã được bàn luận. Điểm thú vị là việc sử dụng CIP không làm thay đổi
số lượng bậc tự do và thuộc tính Kronecker vẫn được đảm bảo. CIP có nhiều tiềm
năng để tiếp tục nghiên cứu chuyên sâu.
CHƯƠNG 4
4.1
PHÂN TÍCH NỨT CƠ-NHIỆT ĐÀN HỒI TUYẾN TÍNH
HAI CHIỀU VỚI VẬT LIỆU ĐẲNG HƯỚNG, TRONG
ĐIỀU KIỆN TẢI ĐỘNG LỰC VÀ TẢI TỰA TĨNH
Mở đầu
Trong chương này, một kiểu phần tử mới, gọi là “XCQ4”, sẽ được giới thiệu và
đặc điểm của nó khi phân tích bài toán nứt cơ-nhiệt đàn hồi tuyến tính sẽ được
khảo sát. Phần tử XCQ4 là phiên bản mở rộng của phần tử CQ4 (xem Chương 3)
bằng cách tích hợp thêm các thành phần làm giàu để biểu diễn ứng xử của vật
rắn khi có vết nứt. Các trường hợp tải tĩnh và tải động lực sẽ được xem xét. Tuy
nhiên với tải động lực, chưa khảo sát trường hợp vết nứt phát triển.
4.2
Mô hình hóa vết nứt
Có hai cách tiếp cận khi xây dựng mô hình vết nứt: mô hình “vết nứt sắc cạnh”
và mô hình “vùng nứt”. Với mô hình “vết nứt sắc cạnh”, vết nứt sẽ được mô tả
là một dạng bất liên tục. Trái lại, cách tiếp cận kiểu “vùng nứt” sẽ mô hình hóa
trạng thái vật liệu bằng một tham số biểu diễn sự hư hại với giá trị trong khoảng
từ 0 (hoàn toàn nguyên vẹn) đến 1 (hư hại hoàn toàn). Nghiên cứu hiện tại giới
hạn trong việc mô tả “vết nứt sắc cạnh” thông qua các hàm làm giàu. Cụ thể,
7
phần tử tứ giác nội suy liên tiếp (CQ4) sẽ được mở rộng bằng các thành phần làm
giàu, để tạo nên phần tử XCQ4 element.
4.3
Phần tử tứ giác nội suy liên tiếp mở rộng (XCQ4)
Phần tử XCQ4 đã được khảo sát trong các tài liệu [17, 18] cho bài toán cơ học
nứt đàn hồi tuyến tính – LEFM - (chưa xét nhiệt độ). So sánh với phiên bản của
phương pháp phần tử hữu hạn mở rộng, gọi là phần tử XQ4, thì XCQ4 được báo
cáo là cho kết quả tốt hơn. Với cảm hứng từ những ưu điểm của XCQ4 trong
LEFM, phần tử tiếp tục được mở rộng để khảo sát bài toán nứt cơ-nhiệt đàn hồi
tuyến tính.
4.3.1 Biểu thức làm giàu để xấp xỉ trường chuyển vị
Biểu thức làm giàu để xấp xỉ trường chuyển vị được viết như sau
u(x) RI (x)u I RJ (x) H (x) H (x J ) a J
I W
J Wsplit
K Wtip
4
RK (x) BL (x) BL (x K ) b KL
(4.1)
L 1
trong đó, số hạng đầu tiên chính là biểu thức xấp xỉ quen thuộc theo CIP, xem
Chương 3. Bước nhảy của trường chuyển vị khi qua vết nứt được biểu diễn bởi
hàm Heaviside H(x), còn sự tập trung ứng suất quanh đỉnh vết nứt được mô tả
bởi bốn hàm nhánh được đề xuất dựa trên lời giải giải tích [19]
B(x) r sin , r cos , r sin sin , r cos sin .
2
2
2
2
(4.2)
(r, θ) là tọa độ cực với gốc đặt tại đỉnh vết nứt, trong đó r là khoảng cách từ điểm
x bất kỳ đến gốc tọa độ và θ là góc nghiêng giữa tiếp tuyến của đoạn vết nứt tại
đỉnh và vector (x – xTIP). Ở phương trình (4.1), W là tập hợp tất cả các nút; Wtip
là tập hợp các nút thuộc các phần tử chứa đỉnh vết nứt; và Wsplit là tập hợp tất cả
các nút thuộc các phần tử bị vết nứt cắt qua nhưng đồng thời không nằm trong
Wtip. aJ và bKL là các bậc tự do tăng thêm do có thêm các thành phần làm giàu.
4.3.2 Biểu thức làm giàu để xấp xỉ trường nhiệt độ
8
4.3.2.1 Trường hợp bề mặt vết nứt đoạn nhiệt
Bề mặt vết nứt đoạn nhiệt nghĩa là nhiệt độ có bước nhảy khi qua vết nứt, và
thông lượng nhiệt có hiện tượng suy biến quanh đỉnh vết nứt, tương tự như với
bước nhảy của chuyển vị và sự suy biến của ứng suất. Vì vậy, biểu thức xấp xỉ
nhiệt độ có thể viết tương tự như phương trình (4.1)
T (x) RI (x)TI RJ (x) H (x) H (x J ) cJ
I W
J Wsplit
K Wtip
RK (x) B1 (x) B1 (x K ) d K .
(4.3)
Dựa trên lời giải giải tích [20], hàm làm giàu đỉnh vết nứt được chọn là hàm
nhánh thứ nhất B1 trong phương trình (4.2).
4.3.2.2 Trường hợp bề mặt vết nứt đẳng nhiệt
Với điều kiện đẳng nhiệt, nhiệt độ trên bề mặt vết nứt được giữ cố định. Khi đó,
bước nhảy không xuất hiện ở trường nhiệt độ, mà ở trường thông lượng nhiệt.
Một dạng hàm làm giàu phù hợp cho trường hợp này theo gợi ý của [20]
(4.4)
L(x) Ri (x) | i | | (x) | .
i
Biểu thức xấp xỉ nhiệt độ vẫn có dạng như phương trình (4.3), nhưng hàm
Heaviside được thay bởi hàm L(x) ở phương trình (4.4), và hàm làm giàu đỉnh
vết nứt được chọn là hàm nhánh thứ hai của phương trình (4.2), dựa trên lời giải
giải tích tham khảo ở [20].
4.4
Tính toán hệ số cường độ ứng suất cho bài toán nứt cơ-nhiệt đàn hồi
tuyến tính
Kỹ thuật tích phân tương tác [19, 21, 22, 23] được sử dụng để trích xuất hệ số
cường độ ứng suất SIFs. Với bài toán nứt cơ-nhiệt đàn hồi tuyến tính cho vật liệu
đẳng hướng, tích phân tương tác được tính theo biểu thức
q
u (1)
u (2)
M (1,2) ij(2) i ij(1) i ik(1) ik(2)1 j
dA
x1
x1
x j
A
9
(4.11)
(1)
(1)
2ui(1) ui(2)
(2) T
(2) T
q
d
A
kk
kk x1 qdA
t 2 x1
x1
A
A0
với
E
.
1 2
Hình 4.1. Miền tính tích phân tương tác
Miền tính tích phân tương tác được thể hiện ở Hình 4.1, với C và C0 là hai đường
cong bất kỳ bao quanh đỉnh vết nứt, nhưng không chứa thêm bất kỳ dạng bất liên
tục nào khác. Trạng thái (1) là trạng thái thật, còn trạng thái (2) được chọn từ lời
giải giải tích cho trường hợp mode-I hoặc mode-II thuần túy. Hệ số SIFs được
trích xuất từ tích phân tương tác
2
M (1,2) * K I(1) K I(2) K II(1) K II(2) .
E
(4.13)
với E* được cho như
Ứng suất phẳng: E * E , hoặc Biến dạng phẳng: E *
4.5
E
1 2
(4.14)
Mô hình vết nứt phát triển
Tiêu chuẩn ứng suất pháp hướng chu vi (MTS) [24] được sử dụng để dự đoán
hướng phát triển của vết nứt. Theo tiêu chuẩn này, vết nứt sẽ trở nên bất ổn định
và phát triển từ đỉnh hiện tại theo hướng θc sao cho ứng suất pháp hướng chu vi
σθθ đạt lớn nhất.
10
4.6
Kết quả tính toán và bàn luận
4.6.4 Mô phỏng sự lan truyền nứt dưới tải tĩnh trong miền hình chữ thập từ
một vết nứt cạnh xiên góc ban đầu
Hình 4.3. Ví dụ 4.6.1: Miền
hình học và chia lưới của mẫu
chữ thập (AF = BE = HC = GD
= 3L)
Hình 4.4. Ví dụ 4.6.1: Đường
nứt dự đoán ứng với 4 trường
hợp tải cho trong Bảng 4.5, sử
dụng phần tử XCQ4
Bảng 4.1. Ví dụ 4.6.1: Bốn trường hợp tải khảo sát
TAB
TCD
TEF
TGH
σ0
Case 1
10o
0o
-10o
0o
0
Case 2
0o
0o
0o
0o
10
Case 3
10o
0o
-10o
0o
10
Case 4
20o
0o
-20o
0o
10
11
Bài toán này khảo sát hiện tượng vết nứt phát triển trong mẫu hình chữ thập chịu
tải tĩnh, xuất phát từ một đoạn nứt xiên góc có sẵn ban đầu. Hình học và điều
kiện biên được cho ở Hình 4.3, với L = 1, chiều dài đoạn nứt ban đầu a = 0.2 L
và góc nghiêng giữa vết nứt ban đầu so với phương ngang là 135o. Thông số vật
liệu [25]: Modulus đàn hồi Young E = 218400 Pa, hệ số Poisson ν = 0.3 và hệ số
giãn nở nhiệt α = 1.67 x 10-5/oC. Có tất cả bốn trường hợp tải được xem xét, như
trình bày trong Bảng 4.1. Mẫu chữ thập được chia lưới với 2000 phần tử XCQ4.
Hình 4.4 biểu diễn đường nứt dự đoán bởi mô hình đề xuất sử dụng phần tử
XCQ4 và tiêu chuẩn MTS xác định hướng vết nứt phát triển, ứng với bốn trường
hợp tải khảo sát. Các kết quả này đều phù hợp với báo cáo ghi nhận từ các tài
liệu tham khảo [25, 26, 27].
4.6.7 Phân tích hệ số cường độ ứng suất khi xét tải động lực: bài toán vết nứt
cong chịu sốc nhiệt
Ví dụ này khảo sát ứng xử của mẫu hai chiều bằng vật liệu Bismuth với một vết
nứt cong có sẵn, chịu tác động của sốc nhiệt. Mẫu Bismuth hình chữ nhật được
minh họa trong Hình 4.27, với các kích thước như sau: W = 12.0 mm, H1 = 10.0
mm, và H2 = 15.0 mm. Thông số vật liệu Bismuth được tham khảo từ tài liệu
[23]. Ban đầu, mẫu có nhiệt độ T0 = 3.5 K. Sau đó, cạnh bên trái của mẫu đột
ngột được giảm nhiệt độ với mức giảm là ΔT = Te – T0 = -0.2 K. Ngoại trừ biên
trái của mẫu, các biên còn lại đều được cách nhiệt. Ở đây bài toán được giả thiết
là biến dạng phẳng.
Cả hai phương trình cơ và nhiệt được giải với sự lựa chọn kỹ thuật tích phân theo
thời gian phù hợp, cụ thể là phương pháp Euler lùi cho bài toán truyền nhiệt và
phương pháp Newmark cho bài toán cân bằng cơ học. Đáp ứng của tấm được
khảo sát trong khoảng thời gian 4 μs. Từ các thông số vật liệu Bismuth [23], tốc
độ lan truyền sóng đàn hồi có thể tính được là Cd = 2346.43 m/s. Vì thế, sóng sẽ
lan đến đỉnh vết nứt khi t = 2.45 μs. Chú ý rằng sóng đàn hồi này được kích hoạt
từ hiện tượng sốc nhiệt. Hình 4.5 minh họa hình ảnh sóng đàn hồi, thông qua
trường chuyển vị tổng, bắt nguồn từ biên trái của tấm Bismuth, nơi xảy ra sốc
12
nhiệt, và lan truyền dần sang biên phải. Tại thời điểm t = 2.45 μs, sóng chạm đỉnh
vết, phù hợp với ước tính ở trên. Sau đó, sóng tiếp tục lan đến biên phải tấm.
Hình 4.4. Ví dụ 4.6.2: Hình học
và điều kiện biên của tấm
Bismuth hình chữ nhật với vết
nứt cong ban đầu.
4.7
Hình 4.5. Ví dụ 4.6.2: Hình ảnh làn
truyền sóng đàn hồi (thông qua
chuyển vị tổng) trong tấm Bismuth có
vết nứt cong ban đầu tại nhiều thời
điểm khác nhau: (a) t = 0.4 μs, (b) t =
1.2 μs, (c) t = 2.0 μs, (d) t = 2.45 μs,
(e) t = 3.2 μs and (f) t = 4.0 μs
Kết luận
Phần tử XCQ4 đã được mở rộng thành công cho bài toán nứt cơ-nhiệt đàn hồi
tuyến tính trong miền hai chiều. Sự tương tác giữa cơ và nhiệt được xét đến. Biểu
thức tính toán được kiểm chứng qua nhiều ví dụ, gồm cả tải tĩnh và tải động lực.
Ở tất cả các ví dụ số, phần tử XCQ4 đều cho thấy khả năng làm việc tốt hơn so
với phiên bản XFEM, tức là phần tử XQ4, nhờ vào sự nâng cấp trong biểu diễn
trường ứng suất.
13
CHƯƠNG 5
5.1
PHÂN TÍCH NỨT CƠ-NHIỆT ĐÀN HỒI HAI CHIỀU
VỚI VẬT LIỆU TRỰC HƯỚNG TRONG ĐIỀU KIỆN
TẢI TĨNH
Mở đầu
Dữ liệu thực nghiệm về vết nứt lan truyền trong vật liệu trực hướng ghi nhận sự
không nhất quán. Có hai quan sát chính: a) vết nứt luôn luôn phát triển theo một
hướng nhất định, bất chấp điều kiện tải [28, 29], và b) điều kiện tải tác động có
ảnh hưởng lớn đến hướng phát triển vết nứt [30]. Những dữ liệu này phản ánh sự
phức tạp của bài toán lan truyền vết nứt trong vật liệu trực hướng và do đó, cần
thêm nghiên cứu về lĩnh vực này. Chương 5 tiếp tục mở rộng phần tử XCQ4 ở
chương 4 cho trường hợp vật rắn trực hướng có vết nứt chịu tải cơ-nhiệt, với sự
tích hợp về tính định hướng của vật liệu.
5.2
Phần tử XCQ4 cho bài toán nứt cơ-nhiệt đàn hồi tuyến tính với vật
liệu trực hướng
Do thuộc tính vật liệu không còn đẳng hướng, nhiều chi tiết trong các biểu thức
tính toán ở chương 4 cần phải được hiệu chỉnh. Trước hết, các hàm làm giàu cho
cả trường chuyển vị và trường nhiệt độ sẽ được chọn lại. Dự đoán hướng phát
triển của vết nứt là cũng nhiệm vụ quan trọng được đặt ra. Ở đây, giả thiết rằng
hướng chính của thuộc tính nhiệt học trùng với hướng chính của thuộc tính cơ
học.
5.2.1
Phương trình đặc trưng của vật liệu trực hướng
Hình 5.1 minh họa các hệ tọa độ sử dụng trong miền vật thể trực hướng hai chiều.
Phương trình đặc trưng cho vật liệu trực hướng được đưa ra bởi Lekhnitskii [31]
S11 z 4 2S16 z 3 2S12 S66 z 2 2S26 z S 22 0 ,
(5.1)
Với Sij là các thành phần của ma trận độ mềm vật liệu, được tính từ nghịch đảo
của ma trận độ cứng vật liệu. Nghiệm phương trình (5.1) là số phức,
zk zkx izky (k = 1, 2), dưới dạng các cặp liên hợp 𝑧1 , 𝑧̅1 và 𝑧2 , 𝑧̅2 .
14
Hình 5.1. Minh họa các hệ tọa độ: hệ toàn cục X-Y; trục tọa độ địa
phương x-y đặt tại đỉnh vết nứt; và các phương chính 1-2 của vật liệu
5.2.2 Biểu thức xấp xỉ chuyển vị
Biểu thức xấp xỉ chuyển vị trong môi trường vật rắn trực hướng có vết nứt tính
theo phần tử XCQ4 sẽ có dạng tương tự như phương trình (4.1). Tuy nhiên, xét
đến lời giải giải tích của trường chuyển vị và trường ứng suất lân cận đỉnh vết
nứt [32], các hàm nhánh được viết theo tài liệu [33]
B1 r sin
1
2
1
g1 ( ) , B3 r sin
2
g 2 ( )
2
g1 ( ) , B4 r cos
(5.2)
2
g 2 ( )
(5.3)
2
2
trong đó gk và θk (k = 1, 2) là các hàm số kể đến tính định hướng của vật liệu
B2 r cos
gk
cos +z kx sin
2
z ky sin ,
2
z ky sin
cos +z kx sin
k arctan
(5.4)
(5.5)
5.2.3 Biểu thức xấp xỉ nhiệt độ
Tương tự như phần trình bày ở chương 4, hai điều kiện nhiệt học trên bề mặt vết
nứt sẽ được xét, là điều kiện đoạn nhiệt và điều kiện đẳng nhiệt. Ở đây cũng cần
chú ý đến sự định hướng của vật liệu. Điểm khác biệt so với chương 4 nằm ở sự
15
lựa chọn hàm nhánh, cụ thể là hàm nhánh thứ nhất và thứ ba sẽ được chọn cho
trường hợp đoạn nhiệt; còn hàm nhánh thứ hai và thứ tư được chọn cho trường
hợp đẳng nhiệt.
5.3
Trích xuất hệ số cường độ ứng suất từ tích phân tương tác
Sau khi tính tích phân tương tác (cách tính tương tự như ở chương 4), các hệ số
SIFs sẽ được trích xuất dựa trên quan hệ dưới đây [33]
M (1,2) 2c11 K I(1) K I(2) c12 K I(1) K II(2) K I(2) K II(1) 2c22 K II(1) K II(2) ,
(5.6)
với
S11
S22 z1 z2
Im z1 z2
Im
, c22
2
2
z1 z2
1 S11
S
c12 22 Im
Im z1 z2
2
z1 z2 2
c11
5.4
(5.7)
(5.8)
Mô hình hóa sự phát triển của vết nứt
Cahill và cộng sự [29] đề xuất một sự hiệu chỉnh cho tiêu chuẩn MTS để dùng
∗
với vật liệu trực hướng, đó là xác định hướng θc sao cho thành phần ứng suất 𝜎𝜃𝜃
đạt lớn nhất. Đây là giá trị chuẩn hóa từ ứng suất pháp hướng chu vi và được tính
bởi biểu thức
*
E
cos p 1 sin 2 p
E2
, với p current
2
(5.9)
Chú ý rằng β là góc hợp bởi phương 1 và phương X toàn cục (xem Hình 5.1),
trong khi θcurrent là góc nghiêng hiện tại của vết nứt.
Một điều cần nhấn mạnh là tiêu chuẩn ở phương trình (5.9) chỉ mới được kiểm
chứng bởi các thí nghiệm trên vât liệu liên hợp có cốt sợi gia cường theo một
phương, và chỉ chịu tác động cơ học [29]. Tuy nhiên, có kết quả thí nghiệm trên
vật liệu khác [30] cho thấy các quan sát ngược lại. Vì thế, cần thiết phải có thêm
khảo sát về tiêu chuẩn này.
16
5.5
Ví dụ tính toán số
5.5.1 Lan truyền vết nứt trong đĩa tròn bất đẳng hướng chịu tải tựa tĩnh
Thí nghiệm nén đĩa tròn kiểu Brazil do Ke và cộng sự [30, 34] thực hiện trên một
loại vật liệu có tính đẳng hướng ngang sẽ được tái lập bằng mô phỏng số sử dụng
phần tử XCQ4. Ví dụ này được thực hiện nhằm nghiên cứu về khả năng áp dụng
tiêu chuẩn dự đoán hướng phát triển của vết nứt do nhóm Cahill đề xuất (phương
trình (5.9)) trong môi trường vật liệu dị hướng có tỉ lệ E1/E2 xấp xỉ 1, vốn chưa
được tài liệu [29] nhắc đến.
Kích thước mẫu gồm đường kính 2R = 74 mm và bề dày t = 10 mm. Vết khắc
ban đầu có chiều dài 2a được tạo ngay tại tâm đĩa và nghiêng góc θ = 45o so với
phương ngang. Đĩa được nén bởi lực tập trung P như Hình 5.2, và thỏa mãn điều
kiện ứng suất phẳng.
Hình 5.2. Ví dụ 5.5.1: Hình học của thí nghiệm nén đĩa tròn có vết nứt
Thuộc tính vật liệu được lấy theo thí nghiệm [30]: E1 = 67.681 GPa, E2 = 78.302
GPa, G12 = 25.336 GPa, G23 = 25.336 GPa, ν12 = 0.185, ν23 = 0.267. Để thực hiện
mô phỏng, miền bài toán được chia thành 6197 phần tử XCQ4. Kết quả tính toán
phù hợp với dữ liệu thực nghiệm, gợi ý rằng việc sử dụng tiêu chuẩn do nhóm
Cahill là hợp lệ. Như thể hiện ở Hình 5.3, bất chấp các góc nghiêng vật liệu khác
17
nhau, vết nứt đều có xu hướng rẽ vuông góc với cạnh khắc ban đầu rồi tiến dần
về điểm đặt lực nén. Trên thực tế, quan sát này không hề tương phản với các tài
liệu [28, 29], nơi ghi nhận rằng vết nứt có xu hướng lan truyền theo phương cứng
nhất. Thay vào đó, kết quả gợi ý rằng tỉ lệ E1/E2 có sụ ảnh hưởng nhất định đến
hướng phát triển của vết nứt. Khi vật liệu có một phương chiếm ưu thế rõ rệt về
độ cứng, vết nứt sẽ lan theo phương đó. Ngược lại, nếu không phương nào chiếm
ưu thế, điều kiện chịu tải sẽ quyết định hình dáng của đường nứt.
Hình 5.3. Ví dụ 5.5.1: Đường phát triển nứt thu được từ vết khắc ban đầu
nghiêng 45o so với phương ngang, và các góc vật liệu β khác nhau
a) β = 30o, b) β = 45o, c) β = 60o
5.5.2 Lan truyền vết nứt trong mẫu hình chữ nhật dưới tác động của thông
lượng nhiệt không đổi
Với ví dụ này, bài toán vết nứt lan truyền trong mẫu hình chữ nhật dưới tác dụng
của thông lượng nhiệt không đổi sẽ được khảo sát, như trong Hình 5.4. Các kích
thước bao gồm: L/W = 4 và a/W = 0.3. Bề mặt vết nứt là đoạn nhiệt. Dòng nhiệt
không đổi song song với bề mặt vết nứt được tạo ra bằng cách cố định nhiệt độ
ở mặt đáy tấm ở mốc tham chiếu (nhiệt độ phòng), tức là ΔT = T - Tref = 0, trong
khi làm mát mặt trên với ΔT = -T0. Mục tiêu của bài toán là khảo sát tính khả
dụng của tiêu chuẩn Cahill’s criterion (Phương trình (5.9)) đối với vật thể trực
hướng chịu tải nhiệt. Vì thế, sáu bộ thuộc tính vật liệu được xem xét như sau
18
Thủy tinh/Epoxy (Mat0): E1 = 55 GPa, E2 = 21 GPa, G12 = 9.7 GPa, ν12
= 0.25, α11 = 6.30 x 10-6 1/K, and α22 = 2 x 10-5 1/K.
Vật liệu trực hướng giả định 1 (Mat1): E2 = 21 GPa, E1 = E2.
Vật liệu trực hướng giả định 2 (Mat2): E2 = 21 GPa, E1 = 1.5 E2.
Vật liệu trực hướng giả định 3 (Mat3): E2 = 21 GPa, E1 = 2 E2.
Vật liệu trực hướng giả định 4 (Mat4): E2 = 21 GPa, E1 = 10 E2.
Vật liệu đẳng hướng giả định (MatIso): E = 21 GPa, ν = 0.25, α = 6.3 x 10-6 K-1.
Hình 5.4. Ví dụ
5.5.2: Hình học và
điều kiện biên.
Hình 5.5. Ví dụ 5.5.2: Các đường nứt tương ứng với
các loại vật liệu có tỉ lệ E1/E2 khác nhau, trong trường
hợp có cùng góc định hướng β = 60o.
Toàn bộ các trường hợp vật liệu trực hướng giả định đều được thay đổi từ thông
số của vật liệu thật thủy tinh/epoxy (Mat0), bằng cách hiệu chỉnh tỉ lệ E1/E2,
trong khi thông số khác giống như Mat0. Với Mat0, tỉ lệ E1/E2 xấp xỉ 2.62. Toàn
bộ các kết quả số đều được tính ở cùng một lưới chia gồm 49 x 99 phần tử XCQ4.
Hình 5.4 hiển thị quỹ đạo vết nứt khi góc định hướng vật ltiệu là β = 60o. Hình
ảnh này tương thích với số liệu trong Bảng 5.1, và cho thấy sự thích hợp khi áp
dụng tiêu chuẩn Cahill trong dự đoán vết nứt phát triển khi có tương tác nhiệt
đàn hồi. Với trường hợp có một phương rất cứng so với các phương khác, vết nứt
sẽ lan truyền theo phương đó, bất chấp điều kiện tải. Trái lại, khi tính dị hướng
19
không rõ rệt, điều kiện tải sẽ là yếu tố chủ đạo ảnh hưởng tới hướng phát triển
vết nứt.
Bảng 5.1. Ví dụ 5.5.2: Hệ số SIFs theo Mode-I và Mode-II, và góc phát triển
vết nứt θc ứng với bước gia tải đầu tiên, khi xét các góc định hướng vật liệu
β khác nhau
β
Material
MatIso
KI
8.400
KII
0.000
θc
0.00o
Mat1 (E1/E2 = 1)
26.995
0.012
-0.09o
Mat2 (E1/E2 = 1.5)
27.040
0.016
-0.09o
Mat3 (E1/E2 = 2)
27.090
0.018
-0.09o
Mat0 (E1/E2 ≈ 2.62)
27.141
0.020
-0.09o
Mat4 (E1/E2 = 10)
27.368
0.033
-0.09o
Mat1 (E1/E2 = 1)
28.407
0.009
-0.09o
Mat2 (E1/E2 = 1.5)
33.422
1.031
25.32o
Mat3 (E1/E2 = 2)
37.093
1.913
37.57o
Mat0 (E1/E2 ≈ 2.62)
40.545
2.858
44.41o
Mat4 (E1/E2 = 10)
57.086
9.132
56.67o
β = 0o
β = 60o
5.6
Kết luận
Phần tử XCQ4 đã được mở rộng thành công để khảo sát vật rắn trực hướng có
vết nứt dưới tác động của tải cơ-nhiệt. Lần đầu tiên, tính định hướng vật liệu
được đưa vào biểu thức xấp xỉ nhiệt độ. Tính chính xác của mô hình đề xuất được
khảo sát và thể hiện qua nhiều ví dụ số, cho thấy kết quả khả quan. Nghiên cứu
hiện tại cũng lần đầu tiên đưa ra các số liệu cho thấy tiêu chuẩn MTS hiệu chỉnh
[29] phù hợp để dự đoán sự lan truyền vết nứt trong các môi trường vật liệu trực
hướng có tỉ lệ 𝐸1 /𝐸2 đa dạng.
20
- Xem thêm -