Đăng ký Đăng nhập
Trang chủ Giáo dục - Đào tạo Cao đẳng - Đại học Development of new finite elements based on consecutive interpolation for 2d and...

Tài liệu Development of new finite elements based on consecutive interpolation for 2d and 3d thermal mechanical problems

.PDF
31
169
109

Mô tả:

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 eS 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 -

Tài liệu liên quan