Tài liệu Mô hình hóa và mô phỏng ứng dụng xử cơ học của ống và tấm mỏng có kích cơ nano mét

BỘ GIÁO DỤC VÀ ĐÀO TẠO TRƯỜNG ĐẠI HỌC BÁCH KHOA HÀ NỘI Nguyễn Danh Trường MÔ HÌNH HÓA VÀ MÔ PHỎNG ỨNG XỬ CƠ HỌC CỦA ỐNG VÀ TẤM MỎNG CÓ KÍCH CỠ NANO MÉT Chuyên ngành : Cơ kỹ thuật Mã số : 62520101 TÓM TẮT LUẬN ÁN TIẾN SĨ CƠ HỌC Hà Nội - 2015 Công trình được hoàn thành tại: TRƯỜNG ĐẠI HỌC BÁCH KHOA HÀ NỘI Người hướng dẫn khoa học: PGS.TS. Lê Minh Quý Phản biện 1: ………………………………………………… Phản biện 2: ………………………………………………… Phản biện 3: ………………………………………………… Luận án sẽ được bảo vệ trước Hội đồng chấm luận án tiến sĩ cấp Trường họp tại Trường Đại học Bách Khoa Hà Nội Vào hồi …….. giờ, ngày ….... tháng ….... năm ………... Có thể tìm hiểu luận án tại thư viện: 1. Thư viện Tạ Quang Bửu - Trường ĐH Bách Khoa Hà Nội. 2. Thư viện Quốc gia Việt Nam. 1 MỞ ĐẦU Lý do chọn đề tài: Nghiên cứu về vật liệu nano là một trong những lĩnh vực nghiên cứu sôi động nhất trong khoảng hai thập niên trở lại đây. Điểm nhấn quan trọng là vào năm 1991, khi Sumio Iijima trình bày công trình khoa học của ông về quá trình tổng hợp tạo ra ống cácbon nano đa lớp (MWCNTs) đầu tiên (Iijima 1991). Đến năm 1993, ông và cộng sự công bố tiếp việc tổng hợp được ống cácbon nano đơn lớp (SWCNT). Kể từ đó tới nay, hàng loạt các nghiên cứu về ống cácbon nano (CNT) được thực hiện cho thấy những đặc tính rất ưu việt, hứa hẹn CNT sẽ có nhiều ứng dụng quan trọng trong nhiều lĩnh vực khác nhau. Như làm chất gia cường cho vật liệu composite, linh kiện điện tử, pin Lithium ion, siêu tụ điện, bộ cảm ứng ánh sáng, nhiệt, sóng điện từ hoặc nhận biết những hóa chất độc hại với độ nhạy rất cao. Những ứng dụng thú vị vừa nêu mới chỉ dừng ở quy mô phòng thí nghiệm, để có thể sản xuất đại trà đòi hỏi nhiều nghiên cứu hơn nữa. Ngoài CNT, cho đến nay, các nhà khoa học đã phát hiện có thêm nhiều loại vật liệu nano dạng tấm và ống có cấu trúc dạng lục giác tương tự, như: graphene, boron nitride (BN), silicon carbide (SiC), silicene (Si),... Sau khi được phát hiện, việc tìm hiểu, dự đoán ứng xử cơ học của các vật liệu trên là rất cần thiết trước khi sản xuất, ứng dụng chúng. Do đó, nghiên cứu sinh đã chọn đề tài cho luận án tiến sĩ của mình là: “Mô hình hóa và mô phỏng ứng xử cơ học của ống và tấm mỏng có kích cỡ nano mét”. Mục đích, đối tượng và phạm vi nghiên cứu: Mô hình hóa và mô phỏng số tìm ứng xử cơ học như mô đun đàn hồi, mô đun đàn hồi trượt, hệ số Poisson, đường cong ứng suất biến dạng,... của các ống và tấm vật liệu nano đơn lớp có cấu trúc dạng lưới lục giác thông qua mô phỏng các thí nghiệm kéo, trượt và uốn. Một số vật liệu nano được chọn để mô phỏng là: graphene, BN, SiC, BSb, Si. Bên cạnh mô hình lý tưởng, luận án cũng xét tới khuyết tật mất hai nguyên tử liền kề và khuyết tật Stone-Wales xảy ra riêng lẻ tại trung tâm của tấm hoặc ống. Phương pháp nghiên cứu: Cho tới nay, việc tiến hành thực nghiệm trên các vật liệu nano là rất khó khăn và phức tạp. Do đó phương pháp mô phỏng số trên máy tính ngày càng được coi trọng. Mô hình hóa và mô phỏng số các loại vật liệu nano hiện nay thường dùng các phương pháp trong cơ học lượng tử cho kết quả chính xác cao như lý thuyết hàm mật độ (DFT) hay mô phỏng ở cấp độ nguyên tử như là phương pháp động lực phân tử (MD). Ngoài ra, phương pháp phần tử hữu hạn nguyên tử (AFEM) được đề xuất khoảng chục năm trở lại đây cũng cho thấy nhiều ưu điểm đáng chú ý. Trong luận án này, nghiên cứu sinh cùng thầy hướng dẫn đã phát triển phương pháp AFEM để mô phỏng ứng xử cơ học của các vật liệu nano. Kết quả thu được sẽ được kiểm chứng bằng cách so sánh với các phương pháp MD, DFT và nhiều phương pháp tin cậy khác. Quá trình xây dựng mô hình cũng như mô phỏng số tính toán kết quả được nghiên cứu sinh lập trình trên phần mềm Matlab. Ý nghĩa khoa học và thực tiễn của đề tài: Ở kích thước cỡ nano mét, việc tiến hành thực nghiệm trên các vật liệu nano nêu trên là rất khó khăn. Thậm chí có những vật liệu mới được phát hiện tồn tại trên lý thuyết, chưa chế tạo được trên thực tế thì việc thực nghiệm là không thể. Do đó việc sử dụng phương pháp mô hình hóa và mô phỏng giúp ta chuẩn đoán chính xác và nhanh chóng đặc tính của vật liệu nano trước khi triển khai sản xuất chúng. Nó giúp tiết kiệm chi phí cho quá trình thiết kế, sản xuất thử và sản xuất hàng loạt các vật liệu này. Bố cục của luận án: Luận án bao gồm phần mở đầu, 4 chương, kết luận và hướng phát triển của luận án, tài liệu tham khảo và phụ lục. 2 CHƯƠNG 1 TỔNG QUAN 1.1 Vật liệu nano cấu trúc dạng lục giác Năm 1991, Sumio Iijima phát hiện ra ống cácbon nano đa lớp (MWCNTs) trong khi đang tiến hành khảo sát fullerene C60 (Iijima 1991). Đến năm 1993 Iijima và cộng sự tiếp tục báo cáo việc tổng hợp được ống cácbon nano đơn lớp (SWCNT) với đường kính 1nm. Năm 2004, graphene, tấm graphite đơn lớp đầu tiên được bóc tách bởi hai nhà khoa học Kostya Novoselov và Andre Geim. Đến 2010 họ đã được trao giải Nobel Vật lý cho những đóng góp của họ trong việc tạo ra và tiến hành thực nghiệm trên tấm graphene. Ngoài graphene và CNT, cho đến nay đã có thêm nhiều vật liệu nano có cấu trúc dạng lục giác tương tự đã được dự đoán tồn tại trên lý thuyết. Trong đó có tấm và ống BN đã được tổng hợp trên thực tế. Mới đây, năm 2012 tấm SiC với độ dày 0,5–1,5 nm cũng được tạo ra ở quy mô phòng thí nghiệm. Bên cạnh đó Si cũng được dự đoán tồn tại trên lý thuyết ở dạng cấu trúc low-buckled (các nguyên tử Si nằm trên hai mặt phẳng song song với nhau). Si ở dạng tấm và dạng dải hẹp đã được tổng hợp trên nền bạc vào năm 2010 và năm 2012, trên nền Zirconium diboride (ZrB2) năm 2012, trên nền Iridium (Ir) năm 2013. Sau khi được phát hiện, các vật liệu trên cần được tìm hiểu, dự đoán các đặc trưng cơ, lý, hóa để phục vụ cho sản xuất và ứng dụng chúng. 1.2 Một số phương pháp mô phỏng vật liệu nano Những phương pháp thường được dùng để mô phỏng, tính toán đặc trưng của các vật liệu nano có thể được phân chia làm hai nhóm: nhóm tính toán ở cấp độ electron và nhóm tính toán ở cấp độ nguyên tử. Trong đó, ở cấp độ electron thì có phương pháp lý thuyết hàm mật độ (DFT), ở cấp độ nguyên tử thì có phương pháp mô phỏng động lực phân tử (MD) là hai phương pháp chuẩn mực được sử dụng nhiều trong các nghiên cứu lý thuyết hoá học, vật lý và khoa học vật liệu. Nhiều nghiên cứu, phương pháp ra đời sau này thường lấy DFT và MD làm chuẩn so sánh. Và ở luận án này, tác giả đã phát triển phương pháp phần tử hữu hạn nguyên tử (AFEM) để sử dụng trong nghiên cứu. Kết quả thu được từ AFEM được so sánh chủ yếu với kết quả từ DFT, MD để đánh giá độ tin cậy của phương pháp cũng như mô hình. AFEM đã được phát triển để phân tích các vật liệu cấu trúc nano một cách hiệu quả (Liu và cs 2004, Nasdala và cs 2005, Wang và cs 2006, Wackerfuß 2009, Nasdala và cs 2010). Trong AFEM, nguyên tử được coi như là nút và các chuyển vị của nó được xem như là các chuyển vị nút. Mỗi phần tử trong AFEM được xây dựng để mô tả cho hàm thế năng. Ma trận độ cứng của các phần tử đó được thiết lập dựa trên hàm thế năng của chúng. Cũng giống như trong FEM, ma trận độ cứng tổng thể trong AFEM được ghép nối từ các ma trận độ cứng phần tử. Mỗi liên hệ giữa chuyển vị nút (chuyển vị của các nguyên tử) và ngoại lực tác dụng cũng được mô tả bởi một hệ phương trình. Về cơ bản, AFEM là sự kết hợp của phương pháp mô phỏng cấp độ nguyên tử và phương pháp phần tử hữu hạn. AFEM cho ta kết quả chính xác như các phương pháp mô phỏng ở thang nguyên tử như MD đồng thời lại cho ra tốc độ hội tụ nhanh hơn do AFEM sử dụng cả đạo hàm bậc nhất và đạo hàm bậc hai của hàm thế khi tính toán cực tiểu hóa năng lượng tìm vị trí cân bằng của hệ, trong khi MD chỉ sử dụng đạo hàm bậc nhất. Tác giả cùng người hướng dẫn nhận định đây là phương pháp mới, còn nhiều điều có thể phát triển để áp dụng mô phỏng cho các vật liệu nano mới tìm ra. Do đó, AFEM được tác giả sử dụng để mô hình hóa và mô phỏng các vật liệu nano trong luận án này. CHƯƠNG 2 CẤU TRÚC NGUYÊN TỬ VÀ THẾ NĂNG TƯƠNG TÁC 2.1 Cấu trúc hình học tấm và ống vật liệu nano dạng lục giác Tấm graphene cũng như các tấm vật liệu nano có cấu trúc dạng lục giác khác chứa các nguyên tử nằm tại đỉnh của các hình lục giác xếp khít với nhau tạo nên dạng lưới như hình tổ ong (hình 2.1). Có hai dạng tấm là tấm phẳng (các nguyên tử cùng nằm trên một mặt phẳng, hình 2.1a) và tấm low-buckled (các nguyên tử nằm trên hai mặt phẳng song song với 3 nhau, hình 2.1b). Như vậy kích thước hình học một tấm vật liệu nano dạng lục giác được xác định bởi hai thông số độc lập là chiều dài liên kết ban đầu lo và khoảng  Low-buckled . Từ đó   rút ra được các thông số góc liên kết θ và độ dài véc tơ đơn vị a1 , a2 như sau:   a  a1  a2  3 lo2   2Low-buckled    3 2Low-buckled  lo2   2lo2   Hình chiếu bằng (2.1) (2.2)   arccos  a1 a2 Tấm phẳng có Low-buckled  0 , nên theo công  tơ Ch như sau:     C h  OA  na1  ma2 lo 120o (a) a1 a2 y, Armchair thức (2.2), góc liên kết θ luôn bằng 120o. Còn góc liên kết của tấm low-buckled sẽ luôn nhỏ hơn 120o. Năm 1995, Dresselhaus và cộng sự đưa ra một hệ thống các tham số biểu diễn hình học để mô tả cấu trúc tấm graphene trước khi được cuộn thành ống SWCNT và sau này được dùng cho các vật liệu nano cấu trúc lưới lục giác khác (hình 2.2). Khi đó, một ống nano xem như được tạo thành bằng cách cuộn một tấm vật liệu nano theo phương véc Hình chiếu cạnh lo θo (2.3) Trong đó cặp chỉ số  n, m là số bước (b) x, Zigzag ∆Low-buckled Hình 2.1 Hình chiếu bằng và hình chiếu cạnh  của tấm vật liệu nano cấu trúc lục giác: a) tấm   và a1 , a2 là các véc tơ đơn vị. Véctơ C h xác phẳng với góc liên kết luôn là θ=120o; b) tấm  định độ xoắn của ống. Khi Ch trùng với low-buckled với góc liên kết θ<120o. dọc theo liên kết dích dắc của lưới lục giác phương zigzag ta có ống dạng zigzag  n, 0  , (n,0) Zigzag tube  Ch trùng với phương armchair ta có ống dạng armchair  n, n  , còn lại là các ống dạng chiral. Xét về phương diện toán học, một tấm vật liệu nano gồm N nguyên tử có (xi, yi, zi) là tọa độ của nguyên tử thứ i , trong đó phương trục z vuông góc với mặt phẳng tấm. Vị trí của nguyên tử i sau khi cuộn  Xi ,Yi , Zi  Ch na1 sẽ được xác định như sau:  Khi Ch trùng với phương zigzag (phương x)  X i  R sin  xi / R   Yi  yi   Z i  R  cos  xi / R   1 y (2.4) Armchair tấm thành ống là (n,n) Armchair tube φ a1 ma2 a2 Zigzag x Hình 2.2 Thông số hình học tấm vật liệu nano cấu trúc hình lục giác. 4  Khi Ch trùng với phương armchair (phương y): X  x i  i Yi  R sin  yi / R    Z i  R cos  yi / R   1 (2.5) (a) Trong đó, R là bán kính ống: R  Ch 2  a 2 n 2  mn  m2  . (2.6) Như vậy, một ống nano được xác định về kích thước khi biết các thông số là cặp chỉ số  n, m (b) (c) y, Armchair và chiều dài ống L. Bên cạnh cấu trúc lục giác lý tưởng như trên, tấm và ống vật liệu sau khi tổng hợp được cũng thường xuất hiện các khuyết tật là khuyết tật Stone-Wales (SW) (có hai loại: SW1 hình 2.5.a và SW2 hình 2.5.b) và khuyết tật mất đi hai nguyên tử liền kề (hình 2.5.c). Những khuyết tật này ảnh hưởng đáng kể tới các đặc trưng cơ học của các vật liệu nano. Hình 2.5 Mô hình tấm vật liệu có cấu trúc lưới lục giác cấu tạo từ 2 loại nguyên tử. (a) khuyết tật SW1 (một liên kết song song phương amchair quay 90o); b) khuyết tật SW2 (một liên kết nghiêng 60o so với phương armchair quay 90o); (c) khuyết tật mất đi 2 nguyên tử liền kề. x, Zigzag 2.2 Thế năng tương tác giữa các nguyên tử Trong luận án này, hàm thế điều hòa và hàm thế Tersoff được chọn cho mô phỏng tương tác giữa các nguyên tử. Ở biến dạng bé, thế năng biến dạng dài và biến dạng góc giữa các nguyên tử có thể được biểu diễn bởi dạng hàm điều hòa như sau: 2 1 (2.8) k r  lij  lo  2 2 1 (2.12) Eijk  k  ijk   o  2 lij là khoảng cách giữa hai nguyên tử i , j . lij  lo là khoảng cách giữa hai nguyên tử i , j Eij  khi hệ ở trạng thái cân bằng. và jk . ijk là góc có đỉnh là nguyên tử j tạo bởi hai liên kết thẳng ji ijk  o khi hệ ở vị trí cân bằng. Hàm thế điều hòa chỉ phù hợp cho biến dạng bé. Hàm thế Tersoff được đề xuất lần đầu để mô tả năng lượng tương tác giữa các nguyên tố Si, Ge và C (Tersoff 1989). Và sau này được đề xuất cho các nguyên tố khác như: Al, B, N, Ga,... Hàm Tersoff có dạng như sau: E 1 Vij Vij   2 i j j i (2.14a) Trong đó: Vij  f C  rij   f R  rij   bij f A  rij   (2.14b) 5 rij  Rij ; 1,   rij  Rij  1 1 fC  rij     cos   . , R  r  Sij ;  Sij  Rij  ij ij 2 2   0, rij  Sij  f R  rij   Aij exp  ijI rij  ; f A  rij    Bij exp  ijII rij  ,  bij  ij 1  ini  ijni  1 2 ni ;  ij   (2.14c) (2.14d) fC  rik  ik g ijk ; g ijk   1  k i , j ci2 ci2  , di2 d 2   h  cos  2 i i ijk ijI   iI   jI  / 2; ijII   iII   jII  / 2; Aij  Ai Aj ; Bij  Bi B j ; Rij  Ri R j ; Sij  Si S j Bằng việc đặt: Rij0  Sij  Rij  2 ; Dij  Sij  Rij  2 . (2.15) một số tài liệu viết lại hàm ngắt ở phương trình (2.14c) ở dạng sau: 1, rij  Rij0  Dij ;   1 1  r  Rij0  0 0 fC  rij     sin   . ij  , Rij  Dij  rij  Rij  Dij ;  2 2 2 D ij     rij  Rij0  Dij 0, (2.16) CHƯƠNG 3 MÔ HÌNH PHẦN TỬ HỮU HẠN NGUYÊN TỬ 3.1 Cở sở lý thuyết phương pháp phần tử hữu hạn nguyên tử 3.1.1 Thiết lập và giải phương trình trong AFEM  o Xét một hệ N nguyên tử. Gọi ri , ri tương ứng là véc tơ vị trí của nguyên tử i ban đầu và sau khi chuyển vị. ui , fi tương ứng là điều kiện biên chuyển vị và ngoại lực tác dụng lên nguyên tử i . Ta có tổng năng lượng của hệ là (Wackerfuß 2009): ET  r   E  r   N  f .r i i (3.1) i 1 Trong đó: r   r1 , r2 ,..., rN T  là véc tơ vị trí của hệ và E  E  r  là nội năng hay thế năng tương tác của hệ có N nguyên tử. Khi hệ ở trạng thái cân bằng, năng lượng đạt giá trị cực tiểu nên đạo hàm bậc nhất của ET (3.2) 0 r Khai triển chuỗi Taylor hàm năng lượng ET quanh vị trí cân bằng ban đầu r  o  ta có (Liu tổng năng lượng của hệ khi đó phải bằng không. Ta có:   và cs 2004): ET  r   ET r o  ET r 1  2 ET .u  uT 2 rr r r o  .u   , (3.3) r  r o Trong đó u  r  r  o  là chuyển vị đủ nhỏ quanh vị trí cân bằng r  o  . u càng nhỏ khai triển gần đúng (3.3) sẽ càng chính xác. Thế (3.3) vào (3.2) và bỏ qua vô cùng bé bậc cao ta có: ET ET  r r  r r o   2 ET rr Biến đổi phương trình (3.4) ta được: Trong đó: P ET r f  r  r o  (3.4) .u  0 r r  o  E r Ku  P , (3.5) (3.6) o r r   6 là lực của hệ ở trạng thái chưa cân bằng. P sẽ tiến tới 0 khi hệ tiến tới vị trí cân bằng mới. K  2 ET rr  r r o  2 E rr (3.7) o r r   là ma trận độ cứng của hệ ở trạng thái cân bằng ban đầu. Trường hợp hệ tuyến tính dấu bằng ở biểu thức (3.3) và (3.4) sẽ xảy ra, dẫn tới phương trình (3.5) chỉ cần giải một lần cho ra ứng xử của hệ ngay tại gốc cân bằng ban đầu. Trường hợp hệ phi tuyến, phương trình (3.5) cần được giải lặp lại cho đến khi lực của hệ ở trạng thái không cân bằng hoặc chuyển vị của hệ ở bước lặp sau đó tiến tới không. Trong thực hành người ta thường tính lặp cho tới khi chúng nhỏ thua một số dương δ bé cho trước   0 , tức là:  , (3.8) P k    . (3.9) u k hoặc k Ở bước lặp thứ k sẽ tính ra được u   . Sau đó véc tơ vị trí của hệ ở bước k  1 sẽ được tính từ bước thứ k trước đó như sau: (3.10) r  k 1  r  k   u  k  r P  k 1  k 1 , K sẽ được dùng để tính P k 1  k 1 ,K  k 1 thông qua biểu thức (3.6) và (3.7). Sau đó sẽ được thế vào (3.5) để giải ra u  k 1 . Vòng lặp cứ thế lặp lại cho đến khi một trong hai điều kiện (3.8) hoặc (3.9) thỏa mãn. Mặt khác, nếu ở phương trình (3.5) ta đặt vào một điều kiện biên chuyển vị u đủ nhỏ (càng nhỏ càng chính xác) để thỏa mãn ngay hoặc chỉ sau một bước lặp là thỏa mãn điều kiện (3.8) thì ta sẽ thu được kết quả gần đúng đủ chính xác từ phương trình (3.5) ngay sau một bước lặp đầu tiên. Điều này tiết kiệm được thời gian tính toán rất nhiều vì không phải tính lặp lại nhiều lần. Toàn bộ quá trình tính toán thiết lập và giải phương trình đặc trưng nêu trên được tác giả lập trình trên phần mềm Matlab. 3.2 Mô hình phần tử hữu hạn nguyên tử với hàm thế điều hòa 3.2.1 Thiết lập ma trận độ cứng phần tử Bỏ qua các thành phần thế năng tương tác đến từ biến dạng xoắn, đảo liên kết, lực van der Waals và tương tác điện từ, thế năng của hệ nguyên tử có thể được biểu diễn như sau: E  Er  E (3.12c) Trong đó Er là thế năng biến dạng dài và E là thế năng biến dạng góc. Ở biến dạng nhỏ, chúng thường được mô tả bằng các hàm điều hòa như sau: 1 M e e2 (3.13) Cr lij , 2 e 1 2 1 N e E   Cijk ijke . (3.14) 2 e 1 e e Trong đó, lij là biến dạng dài của liên kết thẳng thứ e giữa hai nguyên tử i và j ,  ijk là biến dạng góc giữa hai liên kết thẳng ji và jk (hình 3.8). M và N tương ứng là tổng số biến Er      dạng dài và biến dạng góc. Cre và Cijke là các hằng số lực tương ứng cho biến dạng dài và biến dạng góc. Với những tấm vật liệu cấu tạo bởi hai loại nguyên tử sẽ có hai giá trị Cijk tương e ứng phụ thuộc vào kiểu nguyên tử ở đỉnh của góc có biến dạng. Trong hình 3.8b, nguyên tử i và k là cùng loại. 7 Hình 3.8 Mô hình hai kiểu j i j i phần tử: a) phần tử biến dạng dài lij và biến dạng  ijk  lij lij dài  lij ; b) phần tử biến dạng góc ijk và biến dạng  ijk a) góc ijk . (b) k Khi đó thế năng biến dạng của một phần tử hai nút là: 2 1 Ere  Cre  l eji  . 2 (3.15) Khoảng cách giữa hai nguyên tử l ji có thể được biểu diễn thông qua tọa độ e nguyên tử i và tọa độ x j , y j , z j  của nguyên tử l    x  x    y  y    z  z  .  e 2 ji 2 i j 2 i j của j như sau: 2 i  xi , yi , zi  (3.16) j Vi phân hai vế phương trình (3.16), ta có: l eji l eji   xi  x j  xi   x j    yi  y j  yi   y j    zi  z j  zi   z j  . (3.17) Đặt  xi  q1 ,  yi  q2 ,  zi  q3 và  x j  q4 ,  y j  q5 ,  z j  q6 là các chuyển vị của nguyên tử i và j theo các trục x, y và z tương ứng. Phương trình (3.17) trở thành:  l eji  1  xij  q1  q4   yij  q2  q5   zij  q3  q6   . l eji  Trong đó: xij  xi  x j , (3.18) yij  yi  y j , zij  zi  z j . (3.19) e Ở trạng thái cân bằng ban đầu, ta có l ji  l0 . Thay phương trình (3.18) vào phương trình (3.15) ta có: 2 C re  x q  q4   yij  q2  q5   z ij  q3  q6   . E re  2  ij  1 2  l0  (3.20) Đạo hàm bậc hai hàm thế năng được xác định từ phương trình (3.20) theo chuyển vị nút cho ta được ma trận độ cứng của phần tử hai nút (Zienkiewicz và cs 2005): k r mp  2 Ere  , m, p  1, 2...6. qm q p (3.21) Lập luận tương tự, thế năng của một phần tử ba nút (hình 3.8b) có dạng: 2 1 e (3.22) Cijk  ijke . 2 e Trong đó góc  ijk giữa hai véc tơ l ji (véc tơ nối nguyên tử j tới i ) và l jk (véc tơ nối Ee    nguyên tử j tới k ) có thể được tính từ biểu thức sau: l ji .l jk cos ijke  xij xkj  yij ykj  zij zkj , (3.23) với xij , yij , zij và xkj , ykj , zkj có ký hiệu tương tự như (3.19). Vi phân hai vế của phương trình (3.23) và biến đổi ta có: l  x x l  y y z z  ijke sin  ijke   ji  jk  cos  ijke   ij kj  ij kj  ij kj  l  l l l jk  l ji l jk l ji l jk  ji  ji jk Chú ý rằng:   xij  x kj yij  y kj z ij  z kj        l l l l l ji l jk  ji jk ji jk   (3.24)  xij  q1  q4 ;  yij  q2  q5 ;  zij  q3  q6 ;  xkj  q7  q4 ;  ykj  q8  q5 ; zkj  q9  q6 . (3.25) 8 Trong đó q7 , q8 và q9 là chuyển vị của nguyên tử k theo phương x, y và z tương ứng. Đối với cấu trúc tấm hình lục giác ở điều kiện cân bằng ban đầu, chưa có biến dạng, ta có l ji  l jk  lo và ijke  0  1200 . Chú ý đến phương trình (3.18) và thay (3.25) vào phương trình (3.24) ta có: ijke  1  xij cos 0  x jk   q1  q4    yij cos 0  y jk   q2  q5    zij cos  0  z jk   q3  q6     l02 sin 0    x jk cos  0  xij   q4  q7    y jk cos  0  yij   q5  q8    z jk cos  0  zij   q6  q9     (3.26) Thay phương trình (3.26) vào phương trình (3.22) ta có: Ee  Cijke 2  l02 sin 0  2  xij cos 0  x jk   q1  q4    yij cos 0  y jk   q2  q5    zij cos 0  z jk   q3  q6        x jk cos  0  xij   q4  q7    y jk cos  0  yij   q5  q8    z jk cos 0  zij   q6  q9     2 (3.27) Ma trận độ cứng k của phần tử ba nút biến dạng góc được xác định như sau:   kmp   2 Ee , m, p  1, 2...9. qm q p (3.28) Như vậy, ma trận của phần tử hai nút và ba nút là k r và k  đã hoàn toàn xác định và được đưa ra ở dạng hiển qua các biểu thức (3.21) và (3.28). Chúng chỉ phụ thuộc vào tọa độ của các nguyên tử ở vị trí cân bằng ban đầu. 3.2.2 Thông số hàm thế điều hòa Bảng 3.1 Các hằng số lực và chiều dài liên kết ban đầu của graphene, BN, SiC và BSb ở nhiệt độ 0K. Vật liệu Nguyên tử 1 Nguyên tử 2 Graphene BN SiC BSb C B Si B C N C Sb Cre e C121 e C212 kcal/(mol.Å2) 1051 856 447 331 kcal/(mol.rad2) 195 95 155 52 kcal/(mol.rad2) 195 191 105 123 l0 (Å) 1,42 1,45 1,77 2,12 Trong đó chiều dài liên kết ban đầu được lấy từ tài liệu tham khảo (Sahin và cs 2009); các hằng số lực của mô hình vật liệu graphene được tham khảo từ (Berinskii và cs 2013), của vật liệu BN tham khảo từ (Jiang và cs 2011), và của vật liệu BSb tham khảo từ (Rappe và cs 1992); hằng số lực biến dạng góc trong mô hình vật liệu SiC được tham khảo từ tài liệu (Rappe và cs 1992) còn hằng số lực biến dạng dài của SiC được triển khai từ hàm thế năng Tersoff (Erhart và cs 2005). 3.2.3 Kích thước tấm nano Mô hình tấm vật liệu được xây dựng trong luận án có dạng hình vuông gồm 4032 nguyên tử. Chiều dài liên kết ban đầu của các vật liệu được đưa ra trong bảng 3.1. Mô hình ống được xây dựng bằng cách cuộn từ tấm tương ứng, với chiều dài ống sẽ được thiết lập tùy vào trường hợp cụ thể khi mô phỏng. 3.3 Mô hình phần tử hữu hạn nguyên tử với hàm thế Tersoff 3.3.1 Phần tử và thông số hàm thế Tersoff Hai kiểu phần tử được xây dựng như trong hình 3.9 để mô tả từng thành phần tương ứng trong hàm thế Tersoff. 9 k k i i Hình 3.9 Mô hình hai kiểu phần tử mô tả hàm thế Tersoff: (a) phần tử ba nút; (b) phần tử bốn nút. j j a) b) Bảng 3.3 Bảng thông số hàm thế Tersoff cho các vật liệu graphene, BN, SiC, Si. Vật liệu Thông số A (eV)  I (Å-1) B (eV)  II (Å-1)  β n  c d h Graphene Tham khảo của (Lindsay và cs 2010) 1393,6 3,4879 430 2,2119 1,0 1,0 0,72751 1,25724.10-7 38049 4,3484 -0,93 BN Tham khảo của (Sevik và cs 2011) 1380 3,568 340 2,199 1,0 1,0 0,72751 1,25724.10-7 25000 4,3484 -0,89 SiC Tham khảo của (Erhart và cs 2005) 1779 3,2656 225,1895 1,7681 1,0 1,0 1,0 0,011877 273987 180,314 -0,68 l Si Tham khảo của (Tersoff 1989) 1830 2,4799 471,18 1,7322 1,0 1,0 0,78734 1,1.10-6 1,0039.105 16,217 -0,59825 50 0 40 -1 30 Sử dụng hàm ngắt Gỡ bỏ hàm ngắt 20 Năng lượng, eV Lực kéo, eV/Å 3.3.2 Loại bỏ hàm ngắt -2 Sự dụng hàm ngắt -3 -4 10 -5 0 -6 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Biến dạng Gỡ bỏ hàm ngắt -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Biến dạng (a) (b) Hình 3.10 Mô phỏng tương tác của liên kết B-N: a) năng lượng tương tác; b) lực chống lại biến dạng kéo khi sử dụng (nét liền) và không sử dụng (nét đứt) hàm ngắt. Để giảm tải khối lượng tính toán, hàm ngắt có dạng như phương trình (2.14c) được sử dụng nhiều trong mô phỏng MD với mục đích hạn chế số cặp liên kết. Quan sát hình 3.10 ta thấy, trường hợp sử dụng hàm ngắt khi biến dạng dài của liên kết giữa hai nguyên tử B và N đạt ~32% so với ban đầu thì lực liên kết có xu hướng tăng đột ngột và đạt đỉnh lớn nhất ở ~35%. Sự thay đổi lạ này của lực liên kết cũng được nhắc tới trong nghiên cứu của Belytschko và cộng sự (Belytschko và cs 2002) ở liên kết C-C. Vì hiện tượng này mà một số nghiên cứu sử dụng mô phỏng MD với hàm thế Tersoff sử dụng hàm ngắt đã tính toán giá trị 10 ứng suất, biến dạng phá hủy của BN khá lớn và không hợp lý (Mortazavi và cs 2012, Zhao và cs 2013, Han và cs 2014). Để tránh điều này, một số nghiên cứu đã coi Rij=Sij trong phương trình (2.14c). Tuy nhiên, nghiên cứu Belytschko (Belytschko và cs 2002) cho thấy việc dùng hàm ngắt vẫn ảnh hưởng tới ứng xử phá hủy của vật liệu cho dù giá trị hàm ngắt đã lấy thay đổi ở biến dạng 100%. Để loại bỏ hoàn toàn ảnh hưởng của hàm ngắt tới ứng xử phá hủy của vật liệu, luận án này đã gỡ bỏ hàm ngắt f C rij   khi sử dụng hàm thế Tersoff. 3.3.3 Kích thước tấm nano Thông số chiều dài liên kết ban đầu lo trong tấm phẳng được tham khảo từ nghiên cứu (Le Minh Quy 2014) bằng cách cân bằng hệ ở 0K sử dụng mô phỏng MD. Bảng 3.4 Thông số kích thước tấm graphene, BN và SiC gồm 4032 nguyên tử. Tấm vật liệu Chiều dài theo Chiều dài theo Chiều dài liên kết ban đầu nano phương zigzag phương armchair Luận án Tài liệu (Sahin và cs 2009) Graphene 103,38 101,79 1,441 1,42 BN 103,60 102,17 1,444 1,45 SiC 130,46 128,79 1,816 1,77 Thông số tấm low-buckled Si được liệt kê trong bảng 3.5. Trong đó góc liên kết 0=116,4o được tham khảo từ tài liệu (Sahin và cs 2009). Sau đó sử dụng cực tiểu hóa năng lượng liên kết của một liên kết bằng cách giữ góc 0=116,4o không đổi ta thu được chiều dài liên kết lo=2,313 Å theo hàm thế Tersoff. Bảng 3.5 Thông số kích thước tấm Si. Nguồn trích dẫn 0, (o) a, Å l 0, Å 0, Å AFEM, luận án 116,4 3,932 2,313 0,4444 DFT (Sahin và cs 2009) 116,4 3,83 2,25 0,44 DFT (Qin và cs 2012) 3,83 2,25 0,42 DFT (Zhao 2012) 116,25 2,28 0,45 DFT (Peng và cs 2013) 116,134 3,865 2,277 0,454 MD với hàm thế MEAM (Pei và cs 2014) 4,05 2,5 0,85 MD với hàm thế Reax (Roman và cs 2014) 3,76 2,11 0,69 Mô hình các tấm và ống vật liệu BN, SiC, Si chứa các khuyết tật hay gặp trong vật liệu nano dạng lưới lục giác như khuyết tật Stone-Wales (SW) và khuyết tật mất đi hai nguyên tử liền kề cũng được xây dựng khi cần, để khảo sát ảnh hưởng của chúng tới các đặc trưng cơ học, ứng xử phá hủy của các vật liệu. CHƯƠNG 4 KẾT QUẢ VÀ BÀN LUẬN 4.1 Giới thiệu Trong luận án này, mô đun đàn hồi hai chiều Ys (N/m), mô đun đàn hồi trượt hai chiều Gs (N/m) và ứng suất hai chiều  t (N/m) được tính toán chủ yếu. Khi cần thiết tính mô đun đàn hồi Y  Ys / t và mô đun đàn hồi trượt G  Gs / t để so sánh, độ dày tấm được lấy t=0,335nm. 4.2 Kết quả và bàn luận của mô hình sử dụng hàm thế điều hòa 4.2.1 Kéo và trượt thuần túy tấm graphene, BN, SiC và BSb Trong mục này, bài toán kéo đơn trục theo hai phương armchair, zigzag và bài toán trượt thuần túy được mô phỏng tìm đặc trưng đàn hồi của bốn tấm vật liệu: graphene, BN, SiC, BSb. Kết quả được đưa ra trong bảng 4.1 cho thấy rằng, sử dụng cùng một hàm thế, mô đun đàn hồi và mô đun đàn hồi trượt tính bởi phương pháp AFEM sai lệch <5% so với kết quả tính bởi phương pháp MD lấy từ nghiên cứu (Le Minh Quy và cs 2015). Điều này cho thấy tính đúng đắn của phương pháp AFEM được đề xuất trong nghiên cứu này. 11 Bảng 4.1 Đặc trưng đàn hồi của tấm graphene, BN, SiC, BSb được nghiên cứu bởi luận án và so sánh với kết quả của những nghiên cứu khác. Vật liệu Phương pháp nghiên cứu, tác giả Graphene AFEM, luận án MD, (Le Minh Quy và cs 2015) DFT mô hình vỏ,(Kudin và cs 2001) Mô hình Hyperelastic & DFT, (Xu và cs 2012) DFT, (Liu và cs 2007) DFT,(Sahin và cs 2009) Mô hình phi tuyến liên tục & DFT, (Wei và cs 2009) Mô hình vỏ với hàm thế TersoffBrenner, (Yakobson và cs 1996) Mô hình động lực tinh thể, (Popov và cs 2000) Thí nghiệm tĩnh và cộng hưởng âm, siêu âm, (Blakslee và cs 1970) Tán xạ tia X, (Bosak và cs 2007) Thực nghiệm, Lee và cộng sự (2008) (Lee và cs 2008) BN SiC BSb Mô đun đàn hồi Ys, N/m 353 (zig) 347 (arm) 364 (zig) 352 (arm) 345 350 351 335 348 158 Hệ số Poisson 0,179 (zig) 0,175 (arm) 0,15-0,17 150 143 0,149 0,221 149 0,186 0,16 0,169 360 0,19 Y=1002 GPa G=414 GPa 0,21 Y=1029 GPa G=440 GPa 0,17 Y=1092 GPa G=485 GPa 340 ± 50 271 (zig) 266 (arm) 280 (zig) MD, (Le Minh Quy và cs 2015) 271 (arm) DFT,mô hình vỏ, (Kudin và cs 2001) 271 DFT, (Sahin và cs 2009) 267 Mô hình phi tuyến liên tục & 278 DFT, (Peng và cs 2012) DFT & phương trình cân bằng, 275,8 (Andrew và cs 2012) MD, (Zhang và cs 2011) 267 Tán xạ tia X, (Bosak và cs 2006) Y=776 GPa 151 (zig) AFEM, luận án 148 (arm) 156 (zig) MD, (Le Minh Quy và cs 2015) 150 (arm) DFT, (Sahin và cs 2009) 166 DFT & phương trình cân bằng, 163,5 (Andrew và cs 2012) 92 (zig) AFEM, luận án 90 (arm) 95 (zig) MD, (Le Minh Quy và cs 2015) 92 (arm) DFT, (Sahin và cs 2009) 91 AFEM, luận án Mô đun trượt Gs, N/m 152 114 117 112 114 116 G=321 GPa 65 68 0,226 (zig) 0,222 (arm) 0,19-0,22 0,211 0,21 0,225 0.15 0,21 0,176 (zig) 0,172 (arm) 0,16-0,17 0,29 36 37 0,327 (zig) 0,320 (arm) 0,18-0,27 0,34 12 Bảng 4.1 cho thấy, bằng việc sử dụng phương pháp AFEM, đặc trưng đàn hồi của các tấm graphene, BN, SiC, BSb tính ra trong luận án này rất phù hợp với các kết quả của những công trình sử dụng những phương pháp chuẩn mực khác như MD, DFT. Điều này cho thấy rằng mô hình phần tử hữu hạn, ma trận độ cứng mà luận án đã xây dựng được là hoàn toàn đáng tin cậy. Kết quả này đóng góp thêm một phương pháp tính toán đặc trưng đàn hồi của các vật liệu nano một cách nhanh chóng mà vẫn đảm bảo độ chính xác như những phương pháp chuẩn mực khác như DFT và MD. 4.2.2 Kéo ống C, BN, SiC, BSb Đầu tiên, để thấy được sự phụ thuộc của giá trị mô đun đàn hồi ống vào chiều dài ống, hai ống SWCNT (10,0) và (10,10) với các chiều dài thay đổi được chọn để tính toán. Bảng 4.2 Mô đun đàn hồi hai chiều Ys (N/m) của SWCNT phụ thuộc vào tỷ lệ chiều dài (L) trên đường kính ống (D). Kích thước (nm) (10,0), D=0,783 (10,10), D=1,356 L = 5D 337,7 347,6 L = 10D 339,5 348,2 L = 15D 340,1 348,4 L = 20D 340,4 348,5 L = 22D 340,5 348,5 Từ bảng 4.2 ta thấy mô đun đàn hồi của cả ống armchair (10,10) và ống zigzag (10,0) đều tăng nhẹ khi chiều dài của ống tăng. Tuy nhiên, khi chiều dài tăng đủ lớn thì giá trị mô đun đàn hồi của ống dần ổn định. Cụ thể từ giá trị L=15D trở đi, mô đun đàn hồi không còn thay đổi nhiều. Do đó tỷ lệ L=15D được chọn cố định cho các ống có đường kính khác nhau để tìm sự ảnh hưởng của đường kính tới giá trị mô đun đàn hồi của ống. Mô đun đàn hồi hai chiều Ys, N/m 400 CNT(n,n) CNT(n,0) BN(n,n) BN(n,0) SiC(n,n) SiC(n,0) BSb(n,n) BSb(n,0) 350 SWCNT 300 250 Ống BN 200 150 Ống SiC 100 Ống BSb 50 0 5 10 15 20 25 Đường kính ống, Å Hình 4.2: Mô đun đàn hồi hai chiều của 4 loại ống vật liệu theo đường kính ống. (Các đường nằm ngang là mô đun đàn hồi của các tấm zigzag (nét liền) và tấm armchair (nét đứt) của vật liệu tương ứng). Hình 4.2 cho thấy, ở cùng một giá trị đường kính, mô đun đàn hồi của ống armchair lớn hơn so với ống zigzag. Tuy vậy cả hai giá trị đều tăng khi đường kính ống tăng, và tiến sát giá trị của vật liệu dạng tấm tương ứng khi đường kính ống đủ lớn. 4.3 Kết quả và bàn luận của mô hình sử dụng hàm thế Tersoff 4.3.1 Kéo tấm graphene, BN, và SiC Kết quả thu được trong mục này được so sánh với các kết quả thu được bằng mô phỏng MD sử dụng cùng một bộ thông số hàm thế bởi tác giả Lê Minh Quý [Le]. 13 AFEM, MDFEM,zig zig. MDFEM,arm arm. AFEM, MD,[Le], zig. zig MD MD,arm. MD [Le], arm Ứng suất hai chiều σt, N/m 40 35 30 25 20 15 40 (a) Ứng suất hai chiều σt, N/m 45 Graphene 10 5 (b) MDFEM,zig zig. AFEM, MDFEM, arm. AFEM, arm MD[Le], [21], zig. MD zig [21], arm arm. MDMD [Le], 35 30 25 20 15 BN 10 5 0 0 0 0.05 0.1 0.15 0.2 Biến dạng kéo 0 0.25 Hình 4.3 So sánh đường cong ứng suất biến 25 Ứng suất hai chiều σt, N/m dạng kéo tính bởi AFEM và MD khi kéo các tấm nguyên: (a) graphene; (b) tấm BN; (c) tấm SiC. 0.05 20 0.1 0.15 0.2 0.25 Biến dạng kéo MDFEM,zig zig. AFEM, MDFEM, arm. AFEM, arm MD zig MD[Le], [21], zig. [21], arm. MDMD[Le], arm 0.3 (c) Quan sát kết quả đường cong ứng suất15 biến dạng của các tấm vật liệu nguyên khi chịu kéo đơn trục trên hình 4.3 ta thấy, kết quả khi 10 tính bằng MD và AFEM gần như trùng khít SiC nhau trước khi tới điểm phá hủy. Bên cạnh đó, mô đun đàn hồi thu được giữa hai phương 5 pháp MD và AFEM sai khác không quá 5% (xem bảng 4.4). Nhưng do mô phỏng MD vẫn 0 sử dụng hàm ngắt nên tương tác giữa các 0 0.05 0.1 0.15 0.2 0.25 0.3 nguyên tử sẽ bị mất đi khi khoảng cách giữa Biến dạng kéo chúng vượt qua giá trị Sij trong hàm ngắt. Trong khi thực tế tương tác giữa chúng có thể vẫn còn. Điều này cho thấy sự hợp lý của AFEM bằng việc không dùng hàm ngắt sẽ cho ra ứng suất và biến dạng phá hủy cao hơn so với khi tính bằng MD. Bảng 4.4 liệt kê chi tiết các kết quả có được của luận án cùng với các kết quả của các nghiên cứu trước đây để so sánh. Bảng 4.4 Đặc trưng cơ học khi kéo của tấm graphene, BN và SiC nguyên tính bởi luận án và các tác giả khác. Vật liệu Phương pháp nghiên cứu, tác giả Graphene AFEM, luận án MD sử dụng hàm thế Tersoff, (Le Minh Quy và cs 2015) DFT, (Liu và cs 2007) BN Mô hình Hyperelastic và DFT,(Xu và cs 2012) Thực nghiệm, (Lee và cs 2008) AFEM, luận án Mô đun đàn hồi Ys, N/m 358 (zig) 350 (arm) 373 (zig) 350 (arm) 351 350 340 ± 50 258 (zig) 251 (arm) Ứng suất phá hủy σt, N/m 40,3 (zig) 37,0 (arm) 41,1 (zig) 28,5 (arm) 40,4 (zig) 36,7 (arm) 40,0 (zig) 36,4 (arm) 42 ± 4 37,7 (zig) 35,5 (arm) Biến dạng phá hủy, % 20,0 % (zig) 18,1% (arm) 19,6% (zig) 11,4% (arm) 26,6% (zig) 19,4% (arm) 24% (zig) 19% (arm) 25% 25,7 % (zig) 26,4% (arm) 14 MD sử dụng hàm thế Tersoff, (Le Minh Quy 2014) DFT, (Kudin và cs 2001) Tán xạ tia X, (Bosak và cs 2006) DFT, (Sahin và cs 2009) 263,4 (zig) 253,3 (arm) 271 Y=776 GPa 267 DFT, (Topsakal và cs 2010) 258 (dải armchair) 267 278 275,8 MD, (Zhang và cs 2011) DFT, (Peng và cs 2012) DFT và phương trình trạng thái, Andrew và cộng sự (2012). AFEM, luận án SiC 174 (zig) 171 (arm) 179,6 (zig) 173,4 (arm) 166 163,5 36,1 (zig) 29,7 (arm) 22,7% (zig) 17,7% (arm) 20,7 17,9 20,5 17,6 24,8 % (zig) 18,5% (arm) 22,8% (zig) 17,4% (arm) (zig) (arm) (zig) (arm) SW1bị SW1 Tấm SW2bị SW2 Tấm 30 Ứng suất hai chiều σt, N/m Ứng suất hai chiều σt, N/m MD sử dụng hàm thế Tersoff (Le Minh Quy 2014) DFT, Sahin và cộng sự (2009) DFT và phương trình trạng thái, Andrew và cộng sự (2012). Mô đun đàn hồi, ứng suất và biến dạng phá hủy của các tấm nguyên tính được trong luận án này rất phù hợp với kết quả khi tính bởi phương pháp mô phỏng DFT, MD và thực nghiệm. Luận án này chỉ xét tới từng khuyết tật riêng lẻ xảy ra ở trung tâm của tấm là: tấm bị khuyết tật SW1, tấm bị khuyết tật SW2 và tấm bị khuyết tật mất đi hai nguyên tử liền kề. 40 40 Pristine Pristine Tấm nguyên Tấm nguyên (b) (a) 2-atom Tấm mấtvacancy 2 nguyên tử Tấm mấtvacancy 2 nguyên tử 2-atom 35 35 25 20 15 Kéo tấm BN phương zigzag 10 5 0 Tấm SW1bị SW1 SW2bị SW2 Tấm 30 25 20 15 Kéo tấm BN phương armchair 10 5 0 0 0.05 0.1 0.15 0.2 Biến dạng kéo 0.25 0.3 0 0.05 0.1 0.15 0.2 0.25 Biến dạng kéo 0.3 Hình 4.4 So sánh đường cong ứng suất biến dạng khi kéo tấm BN trong hai trường hợp tấm nguyên và tấm bị khuyết tật: (a) phương zigzag ; (b) phương armchair. Quan sát hình 4.4 ta thấy đường cong ứng suất biến dạng của tấm BN nguyên và tấm BN bị khuyết tật rất khớp nhau trước khi tới điểm phá hủy. Kết quả của SiC cũng tương tự. Điều này cho thấy mô đun đàn hồi của tấm nguyên và tấm bị khuyết tật gần như tương đương nhau, trong khi tấm bị khuyết tật bị phá hủy sớm hơn so với tấm nguyên. Vị trí xảy ra phá hủy được tổng hợp chi tiết trong bảng 4.5. Các trường hợp tấm nguyên, vị trí xảy ra phá hủy là ở biên kéo (hình 4.5). Trong khi trường hợp tấm bị khuyết tật mất đi hai nguyên tử liền kề thì vị trí phá hủy luôn xảy ra tại vị trí khuyết tật (hình 4.6). Điều này có thể dễ dàng giải thích bởi việc mất đi nguyên tử làm tấm bị yếu đi cục bộ tại vị trí khuyết tật. Trong khi đó, trường hợp khuyết tật SW1 và SW2 số lượng nguyên tử và liên kết vẫn giữ nguyên nên khuyết tật SW1 và SW2 không làm tấm yếu đi nhiều như trường hợp khuyết tật 15 mất đi hai nguyên tử. Bên cạnh đó, liên kết sau khi bị khuyết tật SW1 trở nên song song với phương zigzag còn đóng góp vào việc chống lại lực kéo và khả năng bị phá hủy xảy ra tại vị trí khuyết tật theo phương zigzag. Ngược lại, liên kết sau khi bị khuyết tật SW2 tạo góc 30o so với phương armchair thì đóng góp vào việc chống lại lực kéo và khả năng bị phá hủy xảy ra tại vị trí khuyết tật theo phương armchair. Mặt khác, ở biến dạng nhỏ hệ số Poisson tính được trong nghiên cứu này cho tấm BN là 0,29 trong khi của SiC là 0,18, tức là khi chịu kéo thì phương vuông góc trong tấm BN sẽ bị co lại nhiều hơn trong tấm SiC. Điều này dẫn tới, tấm BN khuyết tật SW1 khi kéo theo phương zigzag và tấm BN khuyết tật SW2 khi kéo theo phương armchair vẫn bị phá hủy ở biên kéo trong khi tấm SiC thì không. Bảng 4.5 So sánh đặc trưng cơ học khi kéo tấm BN, SiC nguyên và bị khuyết tật. Phương kéo, Dạng Mô đun đàn Ứng suất Biến dạng Vị trí bị phá hủy vật liệu khuyết tật hồi Ys , phá hủy σt, phá hủy, N/m N/m % Zigzag Tấm nguyên 257,9 37,7 25,7 Biên kéo BN Mất hai 256,6 31,0 16,9 Vị trí khuyết tật nguyên tử SW1 257,2 35,1 21,1 Biên kéo SW2 257,9 31,4 17,0 Vị trí khuyết tật Armchair Tấm nguyên 251,1 35,5 26,4 Biên kéo BN Mất hai 248,4 30,1 17,9 Vị trí khuyết tật nguyên tử SW1 250,9 31,8 19,9 Vị trí khuyết tật SW2 250,6 33,2 22 Biên kéo Zigzag Tấm nguyên 174,1 20,7 24,8 Biên kéo SiC Mất hai 173,3 15,6 13,0 Vị trí khuyết tật nguyên tử SW1 173,5 18,6 17,8 Vị trí khuyết tật SW2 174,0 14,9 12,0 Vị trí khuyết tật Armchair Tấm nguyên 171,3 17,9 18,5 Biên kéo SiC Mất hai 171,0 15,1 12,5 Vị trí khuyết tật nguyên tử SW1 171,0 14,3 11,3 Vị trí khuyết tật SW2 171,0 17,0 15,9 Vị trí khuyết tật Hình ảnh tấm bị phá hủy được minh họa bởi hình 4.5 và 4.6, ở đó ta quy ước một liên kết vẽ bằng đường màu đỏ nếu khoảng cách liên kết lớn hơn giá trị Sij trong hàm ngắt ở phương trình 2.14c và không được vẽ khi chiều dài lớn hơn hai lần chiều dài liên kết ban đầu. (a) 24,8% (b) 25,0% Hình 4.5 Hình dạng tấm SiC nguyên bị kéo theo phương zigzag: a) tấm chưa xuất hiện phá hủy ở =24,8%; b) tấm bị phá hủy ở =25,0%. 16 a)=13,0%; b)=13,4%; c) =13,7% Hình 4.6 Hình dạng tấm SiC khuyết mất hai nguyên tử ở trung tâm bị kéo đơn trục theo phương zigzag: a) tấm ở biến dạng =13,0% ứng với ứng suất đạt giá trị lớn nhất; b) tấm ở biến dạng =13,4% và c) tấm ở biến dạng =13,7%. Qua đây ta thấy, tấm BN và SiC có cơ chế phá hủy giòn, việc phá hủy xảy ra rất nhanh và đường cong ứng suất biến dạng cũng rơi đột ngột trên đồ thị. Bảng 4.5 cho thấy, mô đun đàn hồi của tấm bị khuyết tật giảm chưa tới 1,5% so với tấm nguyên. Trong khi đó ứng suất và biến dạng phá hủy của tấm bị khuyết tật lại giảm đáng kể. Cụ thể, ứng suất phá hủy của tấm BN bị mất hai nguyên tử giảm từ 15-18% và biến dạng phá hủy giảm 32-34%. Với tấm SiC, khuyết tật mất đi hai nguyên tử làm giảm 16-25% ứng suất phá hủy và 32-48% biến dạng phá hủy. Còn khuyết tật SW làm giảm từ 7-17% ứng suất phá hủy trong tấm BN và từ 5-28% trong tấm SiC. Việc ứng suất và biến dạng phá hủy giảm mạnh trong khi mô đun đàn hồi giảm ít với tấm có khuyết tật nhỏ cũng từng được trình bày trong các nghiên cứu (Zhang và cs 2005, Khare và cs 2007) tính bằng cơ học phân tử (MM) và các nghiên cứu (Chowdhury và cs 2012, Xu và cs 2013) tính bằng động lực phân tử (MD) cho vật liệu graphene và CNT. 4.3.2 Kéo tấm Si Bảng 4.6 Đặc trưng cơ học khi chịu kéo của tấm Si nguyên. Phương pháp nghiên cứu, tác giả Mô đun đàn hồi Ys , N/m 60,0 (zig) 59,7 (arm) 62 63 63,8 62 63,5 (zig) 60,1 (arm) 62,4 (zig) 59,1 (arm) Hệ số Poisson MD sử dụng hàm thế EDIP (Jing và cs 2013) 65,0 (zig) 64,6 (arm) - MD sử dụng tham số lực Reax (Roman và cs 2014) 50,4 (zig) 62,3 (arm) - AFEM, luận án DFT, (Sahin và cs 2009) DFT, (Qin và cs 2012) DFT,(Pei và cs 2014) DFT,(Topsakal và cs 2010) DFT, (Zhao 2012) DFT,(Jing và cs 2013) 0,28 (zig) 0,25 (arm) 0,30 0,31 0,325 0,37 (zig) 0,41 (arm) - Ứng suất phá hủy σt, N/m 8,2 (zig) 7,9 (arm) 7,07 (zig) 5,66 (arm) - Biến dạng phá hủy, % 23,2% (zig) 21,2% (arm) ~ 20% ~23% dải armchair 18% (zig) 14% (arm) - 7,6 (zig) 6,5 (arm) 19,5% (zig) 15,5% (arm) 5,85 (zig) 4,78 (arm) ~18% (zig) 8,86% (arm) Tấm vuông 21nm x 21nm. Tấm vuông 10nm x 10nm. 17 9 dạng khi kéo của tấm Si nguyên. Tính ở biến dạng nhỏ, mô đun đàn hồi hai chiều của tấm Si nguyên theo hai phương khá gần nhau là 60,0 và 59,7 N/m, hệ số Poisson của tấm Si nguyên là 0,28 theo phương zigzag và 0,25 theo phương armchair. Các giá trị nêu trên rất gần với giá trị tính được ở các công trình khác như đã liệt kê trong bảng 4.6. Zigzag 8 Ứng suất hai chiều σt, N/m Hình 4.8 Đồ thị đường cong ứng suất-biến Armchair 7 6 5 4 3 Si, tấm nguyên 2 1 0 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 Biến dạng kéo 9 7 6 5 4 3 Si, zigzag 2 1 Tấm Pristine nguyên Tấm mất vacancy 2 nguyên tử 2-atom Tấm SW1bị SW1 Tấm SW2bị SW2 8 2-atom vacancy Tấm mất 2 nguyên tử SW1bị SW1 Tấm Tấm SW2bị SW2 (a) Ứng suất hai chiều σt, N/m 8 Ứng suất hai chiều σt, N/m 9 Pristine Tấm nguyên 0 7 6 5 4 3 Si, armchair 2 1 (b) 0 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 Biến dạng kéo 0 Hình 4.9 So sánh đường cong ứng suất-biến 0.40 Hệ số Poisson dạng khi kéo tấm Si trường hợp nguyên và bị khuyết tật: a) phương zigzag; b) phương armchair. 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 Biến dạng kéo Tấm nguyên Tấm mất hai nguyên tử Tấm bị SW1 Tấm bị SW2 0.35 Hình 4.9 cho thấy đường cong ứng suất biến dạng của tấm Si khi bị khuyết tật so với tấm Si nguyên gần như giống hệt nhau Si, zigzag trước khi tới điểm phá hủy. Tuy nhiên, 0.30 cũng giống như BN và SiC, giá trị ứng suất và biến dạng phá hủy của tấm Si bị khuyết tật giảm đi đáng kể so với tấm Si nguyên. Si, armchair Bảng 4.7 cho thấy khuyết tật mất đi hai 0.25 nguyên tử vẫn làm tấm Si bị yếu đi nhiều 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 Biến dạng kéo hơn là khuyết tật SW. Mặt khác, như đã giải Hình 4.10 So sánh hệ số Poisson khi kéo tấm thích trên tấm BN và SiC, kiểu SW1 làm tấm Si bị yếu đi theo phương armchair Si trường hợp nguyên và bị khuyết tật. nhiều hơn là theo phương zigzag trong khi SW2 làm tấm Si bị yếu đi theo phương zigzag nhiều hơn là theo phương armchair. Tóm lại, tấm bị khuyết tật mất hai nguyên tử liền kề ở trung tâm tấm làm cho ứng suất phá hủy của tấm giảm 18-20%, biến dạng phá hủy giảm 3335% và tấm bị khuyết tật SW làm ứng suất phá hủy của tấm giảm 5-20% và biến dạng phá hủy giảm 11-38% so với tấm nguyên. Trong khi mô đun đàn hồi của tấm Si khi bị khuyết tật chỉ giảm khoảng 0,5% so với tấm Si nguyên. 18 Bảng 4.7 So sánh đặc trưng khi kéo của tấm Si trường hợp nguyên và bị khuyết tật. Phương kéo Dạng khuyết tật Zigzag Tấm nguyên Mất hai nguyên tử SW1 SW2 Tấm nguyên Mất hai nguyên tử SW1 SW2 Armchair Mô đun đàn hồi Ys , N/m 60,0 59,7 Ứng suất phá hủy σt, N/m 8,23 6,74 Biến dạng phá hủy, % 23,2 15,0 Vị trí bị phá hủy 59,7 59,9 59,7 59,6 6,85 6,61 7,95 6,35 15,2 14,3 21,2 14,1 Biên kéo Biên kéo Biên kéo Vị trí khuyết tật 59,6 59,5 6,79 7,51 15,8 18,9 Vị trí khuyết tật Vị trí khuyết tật Biên kéo Vị trí khuyết tật Về vị trí trên tấm bị phá hủy, nhìn vào sự liệt kê trong bảng 4.7 cho ta thấy, tấm Si nguyên khi kéo đều bị phá hủy ở biên kéo còn tấm Si bị khuyết tật mất đi hai nguyên tử liền kề thì đều bị phá hủy ở vị trí khuyết tật ở cả hai phương kéo zigzag và armchair giống như tấm BN và SiC. Với khuyết tật SW thì khi kéo theo phương zigzag tấm vẫn vị phá hủy tại biên kéo trong khi kéo theo phương armchair thì tấm lại bị phá hủy tại vị trí khuyết tật. Điều này có thể giải thích là do hệ số Poisson theo phương zigzag luôn lớn hơn theo phương armchair (hình 4.10) nên khi kéo theo phương zigzag, phương vuông góc với phương kéo bị co lại nhiều hơn là khi kéo theo phương armchair. 4.3.3 Kéo ống BN Trước hết, để khảo sát sự phụ thuộc của chiều dài ống tới ứng xử khi chịu kéo của ống, ống BN (8,8) và (14,0) với các chiều dài L=10D, 15D và 20D đã được chọn để tính toán. 40 Ứng suất hai chiều σt, N/m 35 30 0.50 BN(8,8) 0.40 25 20 BN(14,0) 15 L=10D L=15D L=20D 0.45 Hệ số Poisson L=10D L=15D L=20D 0.30 0.25 10 0.20 5 0.15 0 BN(14,0) 0.35 BN(8,8) 0.10 0 0.05 0.1 0.15 0.2 0.25 Biến dạng kéo 0.3 0.35 Hình 4.14 Đường cong ứng suất-biến dạng khi kéo đúng tâm của ống BN(8,8) và BN(14,0) ở các chiều dài ống L=10D, 15D và 20D. 0 0.05 0.1 0.15 0.2 0.25 Biến dạng kéo 0.3 0.35 Hình 4.15 Sự phụ thuộc của hệ số Poisson theo biến dạng kéo đúng tâm của ống BN(8,8) và BN(14,0) ở các chiều dài ống L=10D, 15D và 20D. Nhận xét đầu tiên nhìn trên hình 4.14 là đường cong ứng suất-biến dạng ống zigzag (14,0) ở các chiều dài ống L=10D, 15D, 20D gần như trùng khít nhau từ đầu cho tới khi ống bị mất ổn định (điểm tới hạn) và phá hủy, điều này cũng xảy ra tương tự với ống (8,8). Cụ thể ống BN(14,0) với các chiều dài khác nhau đều có biến dạng tới hạn 25,0-25,1% và ứng suất