Luận án Mô phỏng ứng xử cơ học của ống Nano Phốt pho đen bằng phương pháp phần tử hữu hạn nguyên tử

( ) ( )i j i j i j i j i j = + C R AV f r f r b f r (2.7) ( ) i j i j i j i j i j i j i j i j i j i j i j i j 1, 1 1 c o s . , 2 2 0 ,     − = +     −    C r R r R f r R r S S R r S (2.8) ( ) ( ) ( ) ( )ij ij ij ij ij ij ij ijex p ; ex p = − = − − I II R A f r A r f r B r (2.9) 35 ( ) ( ) ( ) ( ) ( ) 1 / 2 ij i j i j i j i j , 2 2 ij 22 2 ij 1 ; 1 o s i i i n n n i C ik ik k k i j i i k i i i k b f r g c c g d d h c         −  = + = = + − + −  ( ) ( )i j i j i j i j i j i j i j i j ; ; ; 2 2 ; I I II II I II i j i j i j i j A A A B B B R R R S S S       + + = = = = = = Nếu đặt: ( ) ( )ij ij ij ij0 ij ij ; 2 2 + − = = S R S R R D (2.10) Hàm ngắt ở phương trình (2.8) sẽ có dạng sau: ( ) 0 ij i j i j 0 ij i j 0 0 ij i j i j i j i j i j i j 0 ij i j i j 1, 1 1 s in . , 2 2 2 0 ,    −   − = + −   +        +  C r R D r R f r R D r R D D r R D (2.11) trong biểu thức (2.7), (2.8) và (2.9): + i, j và k là ký hiệu cho nguyên tử i, j và k; + rij là khoảng cách của liên kết i và j; + ijk là góc giữa 2 liên kết ij và ik; + ( ) ( )ij ij,A Rf r f r tương tứng là các thành phần tương tác cặp hút và đẩy giữa 2 nguyên tử i và j; + i j b là thành phần xét tới ảnh hưởng của các nguyên tử liền kề với nguyên tử i ngoài j; Tuy nhiên, các thành phần tương tác đẩy và hút được định nghĩa như (2.9) với dạng hàm số i j − r e , do đó hàm Tersoff là hàm thế năng tương tác gần và lúc này hàm ngắt được đề xuất như (2.11) để giới hạn số lượng nguyên tử liền kề có tương tác với nhau, việc này giúp giảm thiểu đáng kể khối lượng tính toán khi mô phỏng. 36 2.2.3 Hàm thế Stillinger-Weber Hàm thế Stillinger-Weber được đưa ra vào năm 1989 bởi 2 nhà khoa học Frank H. Stillinger và Thomas A. Weber cho tương tác năng lượng giữa các nguyên tử Si-F [111]. Do đó, hàm thế này được đặt tên theo tên của 2 nhà khoa học đã tìm ra nó và gọi là hàm thế Stillinger-Weber. Đây là hàm thế phù hợp để mô tả năng lượng tương tác giữa các nguyên P-P, Mo-S và Mo-Se [40]. Thế năng E của cấu trúc nguyên tử là tổng năng lượng kéo dãn liên kết Er và năng lượng uốn liên kết E:  = + r E E E (2.12) 2 3 1 1 ,  = = = =  M N r e e E V E V (2.13) ( ) ( ) m a x 4 2 1  −  = − i j i j r r ij V A e B r (2.14) ( ) ( ) ( ) m a x m a x 2 3 c o s c o s      − + −  = − i j i j i j ik ik ik r r r r ijk o V K e (2.15) trong đó: + V2 tương ứng là thế năng kéo dãn liên kết. + V3 là thế năng uốn liên kết. + M và N ký hiệu tổng phần tử kéo dãn liên kết và uốn liên kết. + Bán kính ngắt của hàm ngắt rmaxij, rmaxik được xác định dựa trên cấu trúc vật liệu. + A, K là hệ số. ρ, B, ρij, ρik, θo là 5 thông số hình học. + rij, rik là chiều dài của liên kết ij và ik. + ijk là góc giữa liên kết ij và ik. Chi tiết hơn về hàm thế Stillinger-Weber được đưa ra trong tài liệu tham khảo [40]. Hàm thế Stillinger-Weber cho thấy ưu điểm khi dùng để mô tả tương tác giữa các nguyên tử của cấu trúc vật liệu nano như MoS2, MoSe2, phốt pho đen [38, 83, 102, 133]. Điều này cho thấy việc sử dụng hàm thế Stillinger-Weber để mô hình hóa tương tác giữa các nguyên tử của cấu trúc ống vật liệu nano phốt pho đen là hoàn toàn phù hợp. 37 2.3 Phương pháp phần tử hữu hạn nguyên tử 2.3.1 Giới thiệu về phương pháp phần tử hữu hạn nguyên tử Trong khi việc thí nghiệm ở kích cỡ nanomét vẫn còn nhiều khó khăn. Các phương pháp tính toán mô phỏng bằng máy tính chiếm vai trò quan trọng trong việc nghiên cứu các thuộc tính của vật liệu cấu trúc nano. So với phương pháp mô phỏng MD thì kết quả của AFEM có độ chính xác khá cao, sai số có thể bỏ qua khi so với MD tính ở xấp xỉ 0ºK. Tuy nhiên, MD có tính đến ảnh hưởng của nhiệt độ, còn AFEM thì không. Trong khi đó, AFEM tính nhanh hơn MD vì AFEM dùng đạo hàm bậc hai của năng lượng, tìm chuyển vị bằng cách giải hệ phương trình tuyến tính như FEM, còn MD thì dùng đạo hàm bậc nhất của năng lượng và dùng giải thuật tối ưu hóa để tính chuyển vị của nguyên tử. Phương pháp phần tử hữu hạn nguyên tử (AFEM) đã được phát triển để phân tích cấu trúc vật liệu nano một cách có hiệu quả [67, 80, 81, 121, 126]. Trong AFEM, các nguyên tử và chuyển vị của các nguyên tử được coi như là các nút và chuyển vị nút. Ma trận độ cứng của các phần tử được lắp ráp dựa trên hàm thế năng tương tác giữa các nguyên tử. Giống như phương pháp phần tử hữu hạn (FEM), ma trận độ cứng tổng thể có được bằng cách lắp rắp các ma trận độ cứng phần tử. Do đó, mối quan hệ giữa chuyển vị và lực có được bằng cách giải hệ phương trình tuyến tính. Thủ tục của phương pháp này có một số điểm giống FEM nên gọi là AFEM. 2.3.2 Cơ sở lý thuyết của phương pháp phần tử hữu hạn nguyên tử Phương pháp AFEM đã được Lê Minh Quý và Nguyễn Danh Trường ứng dụng thành công cho nhóm vật liệu BN, SiC, Si với hàm thế Tersoff [1, 51-54]. Nội dung cụ thể của phương pháp được trình bày dưới đây. Giả sử, cấu trúc nguyên tử (hệ) gồm N nguyên tử. Gọi ( 0 ) i r là véc tơ vị trí của nguyên tử i ở thời điểm ban đầu. ri là véc tơ vị trí của nguyên tử i sau khi dịch chuyển trong không gian. Khi đó = − ( 0 ) i i u r r là chuyển vị của nguyên tử i. = − i j i j r r r là véc tơ nối giữa hai nguyên tử i và j. i u là điều kiện biên chuyển vị tác dụng lên nguyên tử thứ i. i f là điệu kiện biên ngoại lực tác dụng lên nguyên tử thứ i. 38 Khi đó tổng năng lượng của hệ sẽ bằng nội năng của hệ trừ đi lượng cân bằng với công của ngoại lực [1]: ( ) 1 ( ) . N i T i i E E = = −  fr r r (2.16) trong đó: -  = 1 2, ,..., T N r r r r là véc tơ vị trí của hệ. - ( )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 tổng năng lượng của hệ khi đó phải bằng không:  =  0 T E r (2.17) Khiển triển Taylor hàm năng lương ET tại ví trí cân bằng ban đầu  =   ( 0 ) ( 0 ) ( 0 ) ( 0 ) 1 2 , , . . . , T N r r r r ta có: ( ) ( )( ) ( ) ( )0 0 2 0 1 . . . 2 TT T T T E E E E = =    + +    r r r r u r r r ur r u (2.18) trong đó, = − ( 0 ) u r r là chuyển vị đủ nhỏ quanh vị trí cân bằng ( 0 )r . u càng nhỏ thì khai triển Taylor của phương trình (2.18) càng chính xác. Thay phương trình (2.18) vào (2.17) và bỏ qua các vô cùng bé bậc cao, ta có: ( ) ( )00 2 . T E E f ==   = −    r rr r r u r r (2.19) Biến đổi (2.19) ta có: =K u P (2.20) trong đó, ta dặt: ( ) ( )00 2 2 ==   = =     T E E r rr r K r r r r , là ma trận độ cứng của hệ ở trạng thái cân bằng (2.21) ( )0 =  = −  E f r r P r , là véc tơ lực của hệ ở trạng thái cần bằng (2.22) Phương trình (2.20) là phương trình đặc trưng của AFEM. Trong trường hợp tổng quát (2.20) là hệ phương trình phi tuyến. Do đó để giải (2.20) ta dùng các 39 phương pháp lặp cho đến khi lực của hệ ở trạng thái cân bằng nhỏ thua  là một số dương đủ nhỏ cho trước. Cụ thể:  ( ) ( ) ( ) 1 2 m a x( , , ..., ) k k k M P P P (2.23) trong đó: ( )k i P là véc tơ lực của nguyên tử thứ i tại bước thứ k. Ở bước lặp thứ k sẽ tính ra được ( )k 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: + = + ( 1) ( ) ( )k k k r r u (2.24) +( 1)k r sẽ được dùng để tính +( 1)kP , +( 1)kK thông qua biểu thức (2.21) và (2.22). Sau đó +( 1)k P , +( 1)kK sẽ được thế vào (2.20) giải ra +( 1)k u . Vòng lặp cứ thế lặp lại cho đến khi điều kiện (2.23) thỏa mãn (xem hình 2.4), giống với giải thuật Newton. Tuy nhiên, nếu ở phương trình (2.20) ta đặt vào một điều kiện u đủ nhỏ (càng nhỏ càng chính xác) để thỏa mãn ngay hoặc chỉ sau một số bước lặp là thỏa mãn điều kiện (2.23) thì ta sẽ thu được kết quả gần đúng đủ chính xác từ phương trình (2.20) 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 (xem hình 2.5). Chương trình tính bằng Matlab trong luận án sử dụng sơ đồ giải thuật lặp này. Điều kiện biên chuyển vị được chọn rất nhỏ (rất mịn) để đảm bảo là không cần có chương trình lặp con trong mỗi bước tính mà vẫn đảm bảo đạt độ chính xác cần thiết của kết quả tính. Cụ thể, điều kiện biên chuyển vị trong luận án chọn bằng 0.001 Å. Tóm lại, nếu ta có một hệ các nguyên tử đang cân bằng ở vị trí ban đầu, ta cần tìm vị trí cân bằng mới dưới tác dụng của điều kiện biên bất kỳ thì sơ đồ khối giải lặp hình 2.4 được sử dụng để cho kết quả nhanh và chính xác. Nếu ta cần xác định một vị trí cân bằng mới ngay sát vị trí cân bằng trước đó thì chỉ cần sử dụng một bước lặp với sơ đồ khối giải lặp với điều kiện biên đủ nhỏ như hình 2.5. Sơ đồ hình 2.5 này rất phù hợp cho việc tìm đường cong quan hệ ứng suất-biến dạng của hệ các nguyên tử. 40 Hình 2.4 Sơ đồ khối chương trình giải lặp với điều kiện biên tùy ý [1] Đúng Xác định các đại lượng: Tính ma trận K, P theo biểu thức (2.21) và (2.22) Bước lặp k=0 Giải hệ phương trình Ku=P, xác định véc tơ chuyển vị u max(|P1(k)|, |P2(k)|,, |PM(k)|) Kết luận vị trí cân bằng mới của hệ: Tính và in các kết quả Vào thông số mô hình vật liệu, hàm thế và vị trí cân bằng ban đầu. Điều kiện biên. Điều kiện dừng vòng lặp. Bước lặp k=k+1 Xác định vị trí mới của hệ Sai 41 Hình 2.5 Sơ đồ khối chương trình giải lặp với điều kiện biên là một chuyển vị đủ nhỏ [1] Xác định các đại lượng: E, r (0) Tính ma trận K, P theo biểu thức (2.21) và (2.22) Giải hệ phương trình Ku=P, xác định véc tơ chuyển vị u Kết luận vị trí cân bằng mới của hệ: Tính và in các kết quả Vào thông số mô hình vật liệu, hàm thế và vị trí cân bằng ban đầu, điều kiện dừng. Vào điều kiện biên đủ nhỏ: đủ nhỏ và chia thành N bước max(|P1(k)|, |P2(k)|,, |PM(k)|) Đúng Sai 42 Để giải thích rõ hơn sự khác nhau giữa hai sơ đồ giải thuật trên, giả sử có hệ F=f(x) đang cần bằng tại vị trí x(0)=0 (hình 2.6 và 2.7). Nếu cần xác định vị trí cân bằng tại điểm A, ta áp dụng sơ đồ khối hình 2.4. Chỉ sau m=4 bước lặp cho ta kết quả tốt như hình 2.6. Mặt khác nếu sử dụng sơ đồ khối hình 2.5 ta sẽ chia khoảng từ x(A) tới x(0) thành n khoảng, lúc này điều kiện biên sẽ là: ( ) ( 0 )A x x n − =u (2.25) Lặp sau n bước sẽ đến được điểm cân bằng A. Rõ ràng n càng lớn cho kết quả càng chính xác. Như trong 2.7 là 6 bước tính. Trong trường hợp này số bước tính theo sơ đồ 2.4 nhỏ hơn số bước tính sử dụng giải thuật lặp với điều kiện biên nhỏ hình 2.5. Nghĩa là lúc này sơ đồ giải thuật hình 2.4 có lợi hơn vì số bước lặp ít hơn. Một điểm lợi thế khác của sơ đồ giải lặp với điều kiện biên là một chuyểm vị đủ nhỏ (hình 2.5) là nó có khả năng xác định được các vị trí cân bằng tại điểm B và C như hình 2.7. Trong khi đó, sơ đồ hình 2.4 lại không thể tính do quá trình lặp không hội tụ. Hình 2.6 Mô phỏng quá trình tìm nghiệm sử dụng sơ đồ lặp với điều kiện biên tùy ý [1] F F 1 O A B C x ( )1 1 x ( )2 1 x ( )4 1 x ( )3 1 x ( ) ( )4 3 3 u = x - x < δ 43 Hình 2.7 Mô phỏng quá trình tìm nghiệm sử dụng sơ đồ lặp với điều kiện biên là một chuyển vị đủ nhỏ [1] Ngược lại, để thiết lập được cả đường cong f(x) tới điểm cân bằng A nếu sử dụng sơ đồ hình 2.4 thì cần phải chia từ điểm cân bằng x(A) tới x(0) thành n khoảng, mỗi giá trị n ta cần có m bước lặp, tức là ta cần phải giải m.n lần để xây dựng được đường cong f(x) gần đúng cần thiết. Trong khi đó, nếu sử dụng sơ đồ tính hình 2.5 ta vẫn chỉ cần n=6 lần như trên. Khi đó, rõ ràng sử dụng sơ đồ hình 2.5 có lợi hơn. Với những ưu điểm của phương pháp giải thuật lặp với điều kiện biên là đủ nhỏ như hình 2.5 như trên, nghiên cứu sinh tiếp tục ứng dụng phương pháp giải thuật lặp với điều kiện biên là một chuyển vị đủ nhỏ để mô phỏng tính toán ứng xử cơ học của ống vật liệu nano phốt pho đen. 2.3.3 Kiểu phần tử trong AFEM Hệ gồm N nguyên tử tương tác với nhau thì một nguyên tử chỉ tác động qua lại với các nguyên tử lân cận chúng, tạo thành nhóm nguyên tử tương tác lân cận. Do đó, hệ N nguyên tử được coi là bao gồm vô số các nhóm nguyên tử tương tác lân cận. Khi đó ta coi mỗi nhóm nguyên tử tương tác lân cận là một phần tử giống như trong FEM và các nguyên tử được coi là các nút. Có hai kiểu xây dựng phần tử trong AFEM là xây dựng phần tử dựa trên hàm thế và xây dựng phần tử dựa trên cấu trúc nguyên tử. Kiểu xây dựng phần tử dựa trên hàm thế: người ta sẽ sử dụng một phần tử m nút để mô tả cho các thành phần của hàm thế. Việc xây dựng phần tử dựa trên hàm thế F x A B O O x1 x2 x3 x4 x5 x6 F1 F1 F2 F3 F4 F5 F6 44 sẽ phụ thuộc trực tiếp vào hàm thế được chọn. Do đó, các vật liệu khác nhau nhưng sử dụng chung một hàm thế thì ta cũng chỉ cần dùng chung một kiểu phần tử để mô tả cho chúng. Kiểu xây dựng phần tử dựa trên cấu trúc nguyên tử: phần tử hữu hạn trong trường hợp này bao gồm một nguyên tử trung tâm cùng với tất cả các nguyên tử khác nằm trong phạm vi tương tác với nguyên tử trung tâm đó. Trong cách xây dựng phần tử này, ta thấy nó không liên quan đến một dạng hàm thế cụ thể nào. Nghĩa là, với cùng một cấu trúc của một vật liệu ta chỉ có thể chọn một kiểu phần tử duy nhất để mô phỏng cho một hàm thế bất kỳ. Phần tử trong trường hợp này sẽ có khá nhiều nút. Các phần tử trùng nhau nhiều nút. Điều này sẽ dẫn đến tính toán ma trận phần tử và ghép nối ma trận tổng thể khó khăn. Trong luận án này, tác giả dùng kiểu xây dựng phần tử dựa trên hàm thế. Ưu điểm là thuận lợi trong việc tính thế năng và ma trận độ cứng của phần tử. 2.3.4 Mô hình phần tử hữu hạn nguyên tử với hàm thế Stillinger-Weber Trong nghiên cứu này, hàm thế SW được sử dụng để mô hình hóa tương tác giữa các nguyên tử P-P [40]. Thế năng E của cấu trúc nguyên tử là tổng năng lượng kéo dãn liên kết r E và năng lượng làm uốn liên kết E  :  = + r E E E (2.26) 2 3 1 1 ,  = = = =  M N r e e E V E V (2.27) ( ) ( ) m a x 4 2 1  −  = − i j i j r r ij V A e B r (2.28) ( ) ( ) ( ) m a x m a x 2 3 c o s c o s      − + −  = − i j i j i j ik ik ik r r r r ijk o V K e (2.29) trong đó V2 tương ứng là thế năng kéo dãn liên kết và V3 là thế năng uốn liên kết. M và N ký hiệu tổng phần tử kéo dãn liên kết và uốn liên kết trên hình 2.8. Bán kính ngắt của hàm ngắt rmaxij, rmaxik được xác định dựa trên cấu trúc vật liệu. A, K là hệ số. ρ, B, ρij, ρik, θo là 5 thông số hình học. rij, rik là chiều dài của bond ij và ik. ijk là góc giữa bond ij và ik. Chúng được xác định như sau: ( ) ( ) ( ) 2 2 2 = − + − + − i j j i j i j i r x x y y z z (2.30) 45 ( ) ( ) ( ) ( ) ( ) ( ) co s .  − − + − − + − − = j i k i j i k i j i k i ijk ij ik x x x x y y y y z z z z r r (2.31) trong đó (xi, yi, zi ), (xj, yj, zj ) và (xk, yk, zk ) thứ tự là các tọa độ của các nguyên tử i, j và k. Có hai kiểu phần tử được sử dụng để mô tả thế năng tương tác giữa các nguyên tử với hàm thế Stillinger-Weber là phần tử kéo dãn liên kết và phần tử uốn liên kết (xem hình 2.8). Hình 2.8 Hai kiểu phần tử sử dụng hàm thế Stillinger-Weber a) Hai nguyên tử (kéo dãn liên kết) b) Ba nguyên tử (uốn liên kết) Hình 2.9 thể hiện hình ảnh phóng to của một nhóm gồm 6 nguyên tử tách ra từ cấu trúc ống nano phốt pho đen. Khi đó sẽ có 5 phần tử kéo dãn liên kết được đánh số thứ tự 1, 2, 3, 4, 5 và 6 phần tử uốn liên kết được đánh số thứ tự là g1, g2, g3, g4, g5, g6 (hình 2.9).  ijk r ij i j j k i a) b) 46 Hình 2.9 Hình ảnh phóng to nhóm gồm 6 nguyên tử tách ra từ cấu trúc ống phốt pho đen Bộ thông số hàm thế SW được lấy từ [40] cho tương tác giữa P-P của phốt pho đen và được trình bày trong bảng 2.2 và 2.3. Bộ thông số này được xác định bằng phương pháp tính toán theo nguyên lý ban đầu với lý thuyết hàm mật độ. Bảng 2.2 Thông số hàm thế Stillinger-Weber cho tương tác giữa 2 nguyên tử A (nN.Å) ρ (Å) B (Å4) ro (Å) rmaxij (Å) P – P 5,7962 0,809 14,287 2,224 2,790 Bảng 2.3 Thông số hàm thế Stillinger-Weber cho tương tác giữa 3 nguyên tử (uốn liên kết). Pb và Pt là các nguyên tử thuộc đường kính trong và đường kính ngoài của ống Loại K (nN.Å) 𝜃o (độ) ρij (Å) ρik (Å) rmaxij (Å) rmaxik (Å) rmaxjk (Å) Pt – Pt – Pt 57,068 96,359 0,809 0,809 2,790 2,790 3,89 Pb – Pb – Pb 57,068 96,359 0,809 0,809 2,790 2,790 3,89 Pt – Pt – Pb 51,162 102,094 0,809 0,809 2,790 2,790 3,89 Pb – Pb – Pt 51,162 102,094 0,809 0,809 2,790 2,790 3,89 47 Phương pháp AFEM với hàm thế Tersoff đã được Lê Minh Quý và Nguyễn Danh Trường áp dụng thành công trong việc kéo ống BN, SiC và BSb, kéo tấm phẳng BN và SiC [1, 51-54]. Kết quả tính AFEM được so sánh với kết quả tính bằng phương pháp DFT và MD là phù hợp. Do vậy, nghiên cứu sinh chọn và tiếp tục ứng dụng phương pháp AFEM để mô phỏng cho các vật liệu cấu trúc nano mới tìm ra là phốt pho đen và kết hợp với hàm thế mới là hàm thế Stillinger-Weber. 2.4. Kết luận Trong chương này, luận án đã trình bày cấu trúc hình học của ống nano phốt pho đen. Trên cơ sở đó, nghiên cứu sinh sẽ lựa chọn hàm thế và xây dựng chương trình tính mô phỏng thí nghiệm kéo và nén ống nano phốt pho đen trong môi trường Matlab để xác định các đặc trưng cơ học của ống nano phốt pho đen. Thêm vào đó, nghiên cứu sinh đã đưa ra và phân tích các dạng hàm thế năng hay dùng để mô hình hóa thế năng tương tác giữa các nguyên tử là hàm thế Tersoff và hàm thế Stillinger-Weber. Ưu nhược điểm của từng hàm thế đã được chỉ rõ. Trên cơ sở đó, nghiên cứu sinh đã chọn hàm thế Stillinger-Weber dùng để mô hình hóa thế năng tương tác giữa các nguyên tử cho cấu trúc vật liệu ống nano phốt pho đen trong luận án của mình. Bên cạnh đó, mô hình và thuật toán bằng phương pháp phần tử hữu hạn nguyên tử đã được trình bày để giải quyết bài toán xác định các đặc trưng cơ học của ống nano phốt pho đen. Trên cơ sở đó, nghiên cứu sinh áp dụng nó để mô phỏng các thí nghiệm kéo và nén ống vật liệu nano phốt pho đen và qua đó xác định các đặc trưng cơ học của vật liệu này. 48 Chương 3 KẾT QUẢ KÉO ỐNG NANO PHỐT PHO ĐEN 3.1 Đánh giá độ tin cậy của chương trình tính kéo ống nano phốt pho đen Để đánh giá độ tin cậy của chương trình tính AFEM, đầu tiên nghiên cứu sinh thực hiện thí nghiệm mô phỏng kéo ống nano phốt pho đen armchair (0, 20) bằng AFEM so sánh với kết quả tính bằng MD của Liao [63] cho ống armchair (0, 20) có cùng đường kính D=29 Å và chiều dài L=165 Å sử dụng cùng dạng hàm thế Stillinger-Weber. Hình 3.1 thể hiện đường cong ứng suất-biến dạng tính bằng phương pháp AFEM so sánh với kết quả tính bằng MD của Liao [63] với cùng hàm thế là Stillinger-Weber. Kết quả cho thấy, đồ thị đường cong ứng suất-biến dạng tính bằng MD [63] và AFEM trùng khít nhau. Sự trùng khớp này khẳng định độ tin cậy của chương trình tính AFEM. Do đó, việc sử dụng chương trình tính AFEM này để xác định cơ tính của ống vật liệu nano phốt pho đen là đáng tin cậy và chính xác. Hình 3.1 Đường cong ứng suất-biến dạng tính bằng AFEM và MD khi kéo ống nano phốt pho đen 49 3.2 Khảo sát ảnh hưởng của chiều dài ống đến đặc trưng cơ học của ống nano phốt pho đen 3.2.1 Ảnh hưởng của chiều dài ống đến đường cong ứng suất-biến dạng của ống nano phốt pho đen chịu kéo Với độ chính xác của chương trình tính matlab AFEM được đề cập ở phần trên, trong mục này nghiên cứu sinh sẽ làm các thí nghiệm mô phỏng kéo một cặp ống nano phốt pho đen có đường kính xấp xỉ nhau là armchair (0, 8) và zigzag (10, 0) với tỷ số chiều dài/đường kính (L/D) thay đổi từ 6-20 để xét ảnh hưởng của chiều dài ống đến đường cong ứng suất-biến dạng. Kết quả trên hình 3.2 và 3.3 chỉ ra rằng, đường cong ứng suất-biến dạng của cả hai kiểu ống armchair (0, 8) và zigzag (10, 0) đều trùng khít nhau khi tỷ số chiều dài/đường kính (L/D) lần lượt là L/D=6, 8, 12, 16 và 20. Điều này, khẳng định rằng đặc trưng cơ học của ống nano phốt pho đen khi chịu kéo dọc trục sẽ không phụ thuộc vào chiều dài ống. Hình 3.2 Đường cong ứng suất-biến dạng của ống armchair (0, 8) khi tỷ số chiều dài/đường kính, L/D=6, 8, 12, 16, 20 50 Hình 3.3 Đường cong ứng suất-biến dạng của ống zigzag (10, 0) khi tỷ số chiều dài/đường kính, L/D=6, 8, 12, 16, 20 3.2.2 Khảo sát ảnh hưởng của chiều dài ống đến giá trị mô đun đàn hồi của ống nano phốt pho đen chịu kéo Trong phần này, thí nghiệm kéo một cặp ống nano phốt pho đen có đường kính xấp xỉ nhau là armchair (0, 8) và zigzag (10, 0) với chiều dài thay đổi L/D=6- 20 được thực hiện. Mục đích của thí nghiệm này là để xác định sự ảnh hưởng của chiều dài ống đến giá trị mô đun đàn hồi ống. Kết quả được chỉ ra trên hình 3.4 và 3.5. Ta thấy, mô đun đàn hồi của cả hai kiểu ống armchair (0, 8) và zigzag (10, 0) thay đổi rất ít (1,3 %) khi tỷ số chiều dài/đường kính (L/D) tăng từ 6-20. Khi đó, ta có thể coi mô đun đàn hồi là không thay đổi khi tỷ số chiều dài ống thay đổi với ống nano phốt pho đen armchair và zigzag chịu kéo dọc trục. 51 Hình 3.4 Mô đun đàn hồi của ống armchair (0, 8) khi tỷ số chiều dài/đường kính (L/D) thay đổi từ 6 đến 20 Hình 3.5 Mô đun đàn hồi của ống zigzag (10, 0) khi tỷ số chiều dài/đường kính (L/D) thay đổi từ 6 đến 20 52 3.2.3 Khảo sát ảnh hưởng của chiều dài ống đến giá trị ứng suất phá hủy của ống nano phốt pho đen chịu kéo Hình 3.6 và 3.7 thể hiện sự thay đổi giá trị của ứng suất phá hủy của ống armchair (0, 8) và zigzag (10, 0) nano phốt pho đen khi tỷ số chiều dài trên đường kính (L/D) thay đổi từ 6-20. Kết quả chỉ ra rằng, ứng suất phá hủy của các hai kiểu ống armchair (0, 8) và zigzag (10, 0) hầu hư không thay đổi khi tỷ số chiều dài/đường kính (L/D) thay đổi từ 2-20. Do đó, ta có thể khẳng định rằng giá trị ứng suất phá hủy của ống nano phốt pho đen không thay đổi khi chiều dài ống thay đổi. Hình 3.6 Ứng suất phá hủy của ống armchair (0, 8) khi tỷ số chiều dài/đường kính (L/D) thay đổi từ 6 đến 20 53 Hình 3.7 Ứng suất phá hủy của ống zigzag (10, 0) khi tỷ số chiều dài/đường kính (L/D) thay đổi từ 6 đến 20 3.2.4 Khảo sát ảnh hưởng của chiều dài ống đến giá trị biến dạng phá hủy của ống nano phốt pho đen chịu kéo Hình 3.8 và 3.9 thể hiện sự thay đổi của giá trị biến dạng phá hủy của một cặp ống có đường kính xấp xỉ nhau là armchair (0, 8) và zigzag (10, 0) khi tỷ số chiều dài/đường kính ống (L/D) thay đổi từ 2-20. Kết quả cho thấy, biến dạng phá hủy của hai kiểu ống không thay đổi khi tỷ số chiều dài/đường kính (L/D) ống thay đổi. Từ các kết quả trên, ta có thể khẳng định rằng cơ tính của ống nano phốt pho armchair và zigzag khi chịu kéo dọc trục là không phụ thuộc vào sự thay đổi của chiều dài ống. Do đó, từ các thí nghiệm về sau nghiên cứu sinh chọn tỷ số chiều dài/đường kính của tất cả các ống nghiên cứu đều là L/D=8. Khi đó, ứng xử cơ học của ống nano phốt pho đen chỉ bị ảnh hưởng của thông số hình học là đường kính ống. Những ảnh hưởng cụ thể của đường kính đến cơ tính của ống nano phốt pho đen sẽ được nghiên cứu sinh trình bày rõ trong phần dưới đây. 54 Hình 3.8 Biến dạng phá hủy của ống armchair (0, 8) khi tỷ số chiều dài/đường kính (L/D) thay đổi từ 6 đến 20 Hình 3.9 Biến dạng phá hủy của ống zigzag (10, 0) khi tỷ số chiều dài/đường kính (L/D) thay đổi từ 6 đến 20 55 3.3 Kết quả kéo ống nano phốt pho đen bằng phương pháp AFEM 3.3.1 Đường cong ứng suất-biến dạng Trong mục này, thí nghiệm kéo 12 cặp ống armchair và zigzag nano phốt pho đen được mô phỏng để xác định các đặc trưng cơ học của cả hai loại ống armchair và zigzag nano phốt pho đen. Tỷ số chiều dài/trên đường kính (L/D) của tất cả các ống đều là L/D=8. σ và ε ký hiệu cho ứng suất kéo và biến dạng kéo. Mô đun đàn hồi Y được xác định từ đạo hàm bậc nhất của đường cong ứng xuất-biến dạng tại biến dạng =0. Cụ thể được xác định bằng công thức: Y   = (3.1) trong đó: σ và ε lần lượt là ứng suất pháp và biến dạng dài trong thí nghiệm kéo đúng tâm. Do việc xác định chiều dày “t” của các tấm và ống vật liệu 2D là khó khăn. Nên dùng Yt và t ký hiệu cho mô đun đàn hồi 2D và ứng suất 2D, (t là chiều dày của ống). Với: + Mô đun đàn hồi hai chiều được xác định bằng công thức: Yt=Y*t (3.2) + Ứng suất hai chiều được xác định bằng công thức: t=*t (3.3) trong đó: t là chiều dày của tấm hoặc ống vật liệu nano. Thông số hình học của 12 ống armchair có cấu hình lần lượt là (0, 8), (0, 10), (0, 12), (0, 14), (0, 15), (0, 16), (0, 18), (0, 19), (0, 20), (0, 21), (0, 22), (0, 24) và 12 ống zigzag nano phốt pho đen có cấu hình lần lượt là (10, 0), (13, 0), (16, 0), (18, 0), (20, 0), (21, 0), (23, 0), (26, 0), (27, 0), (28, 0), (29, 0), (31, 0) được thể hiện trên