Nghiên cứu hiện tượng động học trong hệ SiO2 bằng mô phỏng động lực học phân tử - Luyện Thị Sâm

74 HNUE JOURNAL OF SCIENCE DOI: 10.18173/2354-1059.2018-0007 Natural Sciences 2018, Volume 63, Issue 3, pp. 74-79 This paper is available online at NGHIÊN CỨU HIỆN TƯỢNG ĐỘNG HỌC TRONG HỆ SiO2 BẰNG MÔ PHỎNG ĐỘNG LỰC HỌC PHÂN TỬ Luyện Thị San Viện Vật lí Kĩ thuật, Trường Đại học Bách Khoa Hà Nội Tóm tắt. Chúng tôi sử dụng mô phỏng động lực học phân tử (ĐLHPT) để nghiên cứu các hiện tượng động học trong hệ SiO2 lỏng. Các hiện tượng động học được phân tích thông qua cơ chế chuyển đổi giữa các ô phối trí SiOxSiOx1 và OSiyOSiy1. Kết quả nghiên cứu chỉ ra rằng các chuyển đổi chỉ xảy ra thường xuyên với một vài ô phối trí và có mối tương quan với nhau. Quá trình phá vỡ và hồi phục của các liên kết xảy ra không đồng nhất trong không gian. Sự xuất hiện các đám không chuyển đổi là nguyên nhân gây ra hiện tượng thuyên giảm động học khi nhiệt độ của hệ tiến gần tới nhiệt độ chuyển pha Tg. Từ khóa: Mô phỏng động lực học phân tử, ô phối trí, hệ SiO2 lỏng, thuyên giảm động học. 1. Mở đầu Các chất lỏng có cấu trúc mạng khi được làm lạnh nhanh xuống dưới nhiệt độ nóng chảy sẽ hình thành trạng thái thái vô định hình ở nhiệt độ Tg gọi là nhiệt độ chuyển pha. Càng gần tới nhiệt độ chuyển pha, động học của các chất lỏng này bị thuyên giảm đột ngột trong khi phạm vi nhiệt độ rất nhỏ [1-3] trong khi cấu trúc của vật liệu chỉ có những thay đổi rất nhỏ. Cơ chế gây ra hiện tượng thuyên giảm động học vẫn chưa được làm rõ. Quan sát chuyển động của các hạt keo chỉ ra rằng càng gần tới nhiệt độ chuyển pha, chuyển động của các hạt dường như có mối tương quan với nhau và các hạt keo linh động không di chuyển theo cùng hướng. Tồn tại các vùng tự sắp xếp gồm các hạt keo linh động [4]. Các nghiên cứu mô phỏng gần đây cho các hệ Lennard –Jones đã tìm thấy các bằng chứng cho hiện tượng không đồng nhất động học [5, 6]. Sử dụng trực quan hóa, các công trình [7, 8] đã quan sát thấy tập hợp các nguyên tử linh động và không linh động nhất và chúng có xu hướng kết cụm trong không gian. Trong các công trình nghiên cứu trước [9, 10], chúng tôi đã nghiên cứu cơ chế khuếch tán của hệ SiO2 lỏng thông qua các chuyển đổi SiOxSiOx1 và OSiy OSiy1. Kết quả nghiên cứu chỉ ra cả tốc độ chuyển đổi và sự phân bố trong không gian của các chuyển đổi đều ảnh hưởng tới hệ số khuếch tán của hệ. Trong bài báo này, chúng tôi sẽ thực hiện các mô phỏng số nhằm làm rõ mối liên hệ giữa các chuyển đổi và các hiện tượng động học đặc biệt như động học không đồng nhất và thuyên giảm động học. Ngày nhận bài: 21/2/2017. Ngày sửa bài: 12/4/2017. Ngày nhận đăng: 20/4/2017. Tác giả liên hệ: Luyện Thị San, e-mail: san.luyenthi@hust.edu.vn. Nghiên cứu hiện tượng động học trong hệ SiO2 bằng mô phỏng động lực học phân tử 75 2. Nội dung nghiên cứu 2.1. Phương pháp nghiên cứu Chúng tôi sử dụng thế tương tác Van Beest-Kramer-Van Santen (BKS) để xây dựng các mẫu SiO2 ở các nhiệt độ khác nhau. Loại thế này tuy đơn giản nhưng thực tế cho thấy nó mô phỏng được nhiều tính chất của các hệ SiO2 vô định hình và lỏng. Hàm thế có dạng     2 6exp i j ij ij ij ij ij ij ij q q e U r A B r C r r     (1) U(rij) là thế năng tương tác giữa hai nguyên tử thứ i và j Thế này bao gồm thành phần thế tương tác Cu-lông và phần thế liên quan tới tương tác cộng hóa trị trong phạm vi gần. Đại lượng rij là khoảng cách giữa hai nguyên tử thứ i và j. Các hệ số Aij, Bij, Cij được sử dụng cho mô phỏng hệ SiO2 được liệt kê trong Bảng 1. Mô phỏng ĐLHPT được tiến hành cho hệ SiO2 gồm 666 nguyên tử Si và 1332 nguyên tử O trong đó có sử dụng điều kiện biên tuần hoàn. Thời gian của một bước mô phỏng trong nghiên cứu này là 0,47 fs. Bảng 1. Các hệ số của hàm thế BKS đối với hệ SiO2 [11] Cặp nguyên tử Aij (eV) Bij (Å -1 ) Cij(eV Å 6 ) Điện tích (e) O-O 1388.773 2.760 175.000 qO = − 1.2 Si-O 18003.757 4.873 33.538 qSi = + 2.4 Si-Si 0.0 0.0 0.0 Hai mẫu SiO2 ở các nhiệt độ khác nhau được xây dựng theo trình tự như sau: Cấu hình ban đầu được tạo ra bằng cách gieo ngẫu nhiên tất cả các hạt trong một khối lập phương. Cấu hình ngẫu nhiên ban đầu này được nung nóng tới nhiệt độ 6000 K để hình thành một trạng thái ngẫu nhiên hoàn toàn mới của hệ và giữ trong khoảng 50000 bước mô phỏng. Tiếp theo, cấu hình được làm lạnh xuống tới các nhiệt độ 5500 K, 5000 K, 4500 K, 4000 K và cuối cùng là 3500 K. Ở nhiệt độ 3500 K, mẫu SiO2 được hồi phục trong 10 6 bước mô phỏng với áp suất không của mẫu được giữ không đổi. Từ mẫu SiO2 ở áp suất không và nhiệt độ T = 3500 K tiếp tục giảm nhiệt độ của hệ tới 2600 K bằng cách sử dụng mô hình NPT (số hạt N, áp suất P và nhiệt độ T được giữ không đổi). Cuối cùng, hai mẫu thu được sẽ được hồi phục ở thể tích và năng lượng không đổi bằng cách sử dụng mô hình NVE (số hạt N, thể tích V và năng lượng toàn phần E được giữ không đổi) trong khoảng 107 bước mô phỏng. Sau khi các mẫu vật liệu đạt trạng thái cân bằng, tiến hành xác định các tính chất cấu trúc. động học và tính đa thù hình của hệ. Để tăng độ chính xác cho các kết quả thu được, tất cả các đặc trưng về số phối trí, góc, hàm phân bố xuyên tâm được lấy trung bình từ 1000 cấu hình cuối cùng của mẫu vật liệu và cứ sau 10 bước mô phỏng lại thực hiện một phép đo. Ô phối trí là một tập hợp nguyên tử trong đó bao gồm nguyên tử trung tâm và các nguyên tử lân cận. Các nguyên tử lân cận sẽ hình thành các liên kết với nguyên tử trung tâm. Liên kết giữa nguyên tử Si và O được hình thành nếu khoảng cách giữa hai nguyên tử nhỏ hơn bán kính ngắt (rc = 2,38 Å) tương ứng với vị trí tiểu thứ nhất của hàm phân bố xuyên tâm đối với cặp Si-O. Có hai loại ô phối trí được sử dụng trong nghiên cứu này là SiOx và OSiy trong đó x và y tương ứng là số nguyên tử O và Si. Một chuyển đổi xảy ra khi một liên kết bị phá vỡ hoặc một liên kết mới được tạo thành. Các chuyển đổi giữa các ô phối trí thường xảy ra là SiOx  SiOx1 và OSiy  OSiy1 và rất ít các dạng chuyển đổi khác. Hai dạng chuyển đổi này xảy ra đồng thời với nhau. Luyện Thị San 76 2.2. Kết quả nghiên cứu và thảo luận Để đảm bảo tính tin cậy của mô hình, chúng tôi tiến hành so sánh hàm phân bố xuyên tâm của mẫu SiO2 được xây dựng trong bài báo với các công trình thực nghiệm và mô phỏng trước đây. Kết quả cho thấy hàm phân bố xuyên tâm có sự phù hợp tốt về cả vị trí đỉnh và độ rộng của hàm phân bố [12]. Ở áp suất khí quyển, số phối trí trung bình của Si và O tương ứng là 4 và 2. Như vậy ở áp suất khí quyển, SiO2 có cấu trúc mạng tứ diện. Hình 1 mô tả sự thay đổi số lượng các chuyển đổi theo thời gian mô phỏng ở hai nhiệt độ 3500 K và 2600 K. Cả hai đường đều có dạng tuyến tính. Độ dốc của đường thẳng dùng để xác định số lượng chuyển đổi. Hệ số khuếch tán của nguyên tử Si và O được xác định thông qua phương trình Einstein. Hệ số khuếch tán và tốc độ chuyển đổi cho hai mẫu SiO2 được trình bày trong Bảng 2. Tốc độ chuyển đổi giảm mạnh khi nhiệt độ giảm. Tuy nhiên so với hệ số khuếch tán, tốc độ chuyển đổi giảm chậm hơn. Điều này có nghĩa là chuyển đổi không phải là nguyên nhân duy nhất gây nên hiện tượng khuếch tán. Cách thức các chuyển đổi xảy ra và sự phân bố của chúng trong không gian cũng đóng vai trò quan trọng trong quá trình này. 0 5000 10000 15000 20000 0 10000 20000 30000 40000 T=3500K T=2600K S è l - î n g c ¸ c c h u y Ó n ® æ i Sè b-íc m« pháng Hình 1. Sự phụ thuộc của số lượng các chuyển đổi vào thời gian mô phỏng 0 5000 10000 15000 20000 0 500 1000 1500 2000 T=3500K T=2600K S è l - î n g c ¸ c n g u y ª n t ö t h a m g ia c h u y Ó n ® æ i Sè b-íc m« pháng Hình 2. Sự phụ thuộc của số lượng các nguyên tử tham gia chuyển đổi vào thời gian mô phỏng Nghiên cứu hiện tượng động học trong hệ SiO2 bằng mô phỏng động lực học phân tử 77 Khi số bước mô phỏng tăng (tương ứng với thời gian mô phỏng tăng), số lượng các nguyên tử tham gia chuyển đổi càng lớn. Ở cùng số bước mô phỏng , nhiệt độ càng cao thì số nguyên tử tham gia chuyển đổi càng nhiều (Hình 2). Quá trình chuyển đổi giữa các ô phối trí xảy ra khi có sự phá vỡ các liên kết Si-O trong mạng. Điều này có nghĩa là khi số bước mô phỏng tăng dần, số lượng các liên kết chưa từng bị phá vỡ sẽ giảm dần. Hình 3 cho thấy tốc độ giảm của các liên kết chưa từng bị phá vỡ ở nhiệt độ cao sẽ lớn hơn so với ở nhiệt độ thấp. Tại nhiệt độ T = 3500 K, tồn tại 4372 liên kết chưa từng bị phá vỡ tại n = 20000 bước mô phỏng trong khi ở nhiệt độ T = 2600 K, số liên kết này là 5282. Bảng 2. Hệ số khuếch tán và tốc độ chuyển đổi cho hai mẫu SiO2 Nhiệt độ (K) DSi , cm 2 /s DO , cm 2 /s ξ 2600 1,73×10 -7 2,30×10 -7 4,01×10 -4 3500 2,68×10 -6 3,82×10 -6 1,06×10 -3 Trong mỗi ô phối trí, mỗi nguyên tử Si hình thành 4 liên kết với 4 nguyên tử O xung quanh hoặc một nguyên tử O liên kết với 2 nguyên tử Si. Số lượng liên kết chưa từng bị phá vỡ trong từng ô phối trí ở nhiệt độ T = 3500 K với số bước mô phỏng n = 20000 được trình bày trong Bảng 3. Sau 20000 bước mô phỏng, có 44 nguyên tử Si mà 4 liên kết với nguyên tử O chưa hề bị phá vỡ và 112 nguyên tử Si có 4 liên kết đều bị phá vỡ. Độ dịch chuyển bình phương trung bình và số lượng chuyển đổi của các nguyên tử Si này có sự chênh lệch rõ ràng. Xu hướng tương tự cũng được quan sát thấy từ số liệu thu được đối với O. Như vậy, quá trình chuyển đổi tương ứng với sự phá vỡ các liên kết Si-O trong các ô phối trí xảy ra không ngẫu nhiên. Quá trình này chỉ tập trung ở một số nguyên tử xác định và những liên kết xác định. Trong hệ luôn luôn tồn tại những liên kết rất bền vững và những liên kết liên tục bị phá vỡ do tại đó không xảy chuyển đổi và xảy ra chuyển đổi. 0 5000 10000 15000 20000 2000 2500 3000 3500 4000 4500 5000 5500 T=3500K T=2600K S è l iª n k Õ t c h - a p h ¸ v ì Sè b-íc m« pháng Hình 3. Sự phụ thuộc của số liên kết chưa bị phá vỡ vào thời gian mô phỏng Tiếp theo chúng tôi khảo sát quá trình tạo đám của các nguyên tử. Đám được định nghĩa như một tập hợp các nguyên tử mà mỗi nguyên tử trong đám hình thành ít nhất một liên kết với các nguyên tử còn lại của đám. Kích thước đám được định nghĩa là số lượng nguyên tử của đám. Xét mẫu ở T = 3500 K tại thời điểm ứng với 10000 bước mô phỏng. Các nguyên tử chuyển đổi chia thành 188 đám và các nguyên tử không chuyển đổi chia thành 304 đám. Như vậy, các nguyên tử mà tại đó xảy ra chuyển đổi và không xảy ra chuyển đổi không phân bố đều trong không gian mô phỏng. Chúng có xu hướng tạo thành các đám. Trong 10000 bước mô phỏng, các nguyên tử xảy ra chuyển đổi chia làm 188 đám trong đó đám lớn nhất có 161 nguyên tử. Chúng tôi chọn đám lớn nhất này và theo dõi sự thay đổi nguyên tử của nó trong 10000 bước mô phỏng (Bảng 4). Số liệu thu được chỉ ra rằng đối với các nguyên tử xảy ra chuyển đổi, quá trình tan và hợp của các Luyện Thị San 78 đám xảy ra liên tục. Tuy nhiên, đối với các nguyên tử không xảy ra chuyển đổi, chúng tôi theo dõi đám không chuyển đổi có kích thước 59 nguyên tử trong toàn bộ 10000 bước mô phỏng. Đám này hầu như chỉ dao động tại vị trí ban đầu và vị trí tương quan giữa các nguyên tử thay đổi rất ít. Khi nhiệt độ giảm, các nguyên tử không tham gia chuyển đổi tăng và kích thước các đám không chuyển đổi tăng. Sự mở rộng của các đám không chuyển đổi trong không gian được cho là nguyên nhân gây ra hiện tượng thuyên giảm động học của hệ khi nhiệt độ tiến tới gần nhiệt độ chuyển pha Tg. Bảng 3. Số lượng liên kết chưa từng bị phá vỡ trong từng ô phối trí ở nhiệt độ T = 3500 K vố số bước mô phỏng n = 20000 Nguyên tử Số lượng liên kết được giữ nguyên từ đầu Số lượng nguyên tử Độ dịch chuyển bình phương trung bình/1 nguyên tử Số lượng chuyển đổi/1 nguyên tử Si 0 112 3.1501 53.1339 1 203 2.8086 39.6256 2 186 2.0325 24.3871 3 121 1.5536 14.4876 4 44 1.3408 6.3409 O 0 498 5.0023 25.9558 1 554 2.7307 12.0253 2 280 1.7539 3.5536 a) n = 0 b) n = 5000 c) n = 10000 Hình 4. Đám không chuyển đổi có kích thước 59 nguyên tử tại 3 thời điểm khác nhau (0, 5000 và 1000 bước mô phỏng) Bảng 4. Theo dõi đám chuyển đổi gồm 161 nguyên tử trong 10000 bước mô phỏng Bước mô phỏng 2000 4000 6000 8000 10000 Số lượng đám 17 16 18 13 1 Số nguyên tử của đám lớn nhất 102 121 139 126 161 Nghiên cứu hiện tượng động học trong hệ SiO2 bằng mô phỏng động lực học phân tử 79 3. Kết luận Theo kết quả mô phỏng, các chuyển đổi SiOxSiOx1 and OSiy OSi x1 trong hệ SiO2 ảnh hưởng tới quá trình khuếch tán. Khi nhiệt độ giảm thì hệ số khuếch tán giảm nhanh hơn tốc độ chuyển đổi. Các chuyển đổi không xảy ra ngẫu nhiên trong hệ mà chúng thường xảy ra với một vài nguyên tử và liên kết Si-O. Các nguyên tử xảy ra chuyển đổi và không xảy ra chuyển đổi có xu hướng tạo thành các đám trong không gian. Hiện tượng thuyên giảm động học có thể gây ra bởi sự liên thông của các đám không chuyển đổi trong không gian khi nhiệt độ của hệ tiến tới gần nhiệt độ chuyển pha Tg. Lời cám ơn. Nghiên cứu này được tài trợ bởi trường Đại học Bách khoa Hà Nội với mã số T-2016-PC-213. TÀI LIỆU THAM KHẢO [1] K. L. Ngai, 2000. Dynamic and thermodynamic properties of glass-forming substances. J. Non- Cryst. Solids, 275, pp. 7-51. [2] Sillescu et al., 2002. Heterogeneity at the glass transition: what do we know? J. Non-Cryst. Solids, 307-310, pp. 16-23. [3] G. A. Appignanesi and J. A. Rodriguez Fris, 2009. Space and time dynamical heterogeneity in glassy relaxation. The role of democratic clusters. J. Phys.: Condens. Matter, 21. [4] Eric R. Week et al., 2007. Short- and long-range correlated motion observed in colloidal glasses and liquids. J. Phys.: Condens. Matter, 19. [5] E. Flenner and G. Szamel, 2006. Anisotropic spatially heterogeneous dynamics in a model glass- forming binary mixture. J. Phys.: Condens. Matter, 19. [6] B. Doliwa and A. Heuer, 2000. Cooperativity and spatial correlations near the glass transition: Computer simulation results for hard spheres and disks. Phys. Rev. E, 61. [7] Asaph Widmer-Cooper, 2004. How Reproducible Are Dynamic Heterogeneities in a Supercooled Liquid? Phys. Rev. Lett., 93. [8] K. Deenamma Vargheese, 2010. Origin of dynamical heterogeneities in calcium aluminosilicate liquids. J. Chem.Phys., 132. [9] P. K. Hung and L. T. Vinh, 2006. Local microstructure of liquid and amorphous Al2O3 .J. Non- Cryst. Solids, 352, pp. 5531-5540. [10] P. K. Hung et al, 2016. The study of diffusion in network-forming liquids under pressure and temperature. Physica B, 501, pp. 18- 25. [11] B. van Beest, G. Kramer, R. van Santen, 1990. Force fields for silicas and aluminophosphates based on ab initio calculations. Phys. Rev. Lett., 64. [12] P. K. Hung, N. V. Hong and L. T. Vinh, 2007. Diffusion and structure in silica liquid: a molecular dynamics simulation. J. of Phys: Cond. Matt., 19. ABSTRACT Molecular dynamics simulation of SiO2 under temperature Luyen Thi San Shool of Engineering Physics, Hanoi University of Science and Technology We use the molecular dynamics simulation to study the diffusion in liquid silica. The diffusion is analyzed through reactions SiOxSiOx1, and OSiyOSiy1. It was revealed that they happen more frequently with some coordination cells and bonds Si-O and are strongly correlated. The breakdown and recovery of Si-O bonds are heterogeneous in space. Rigid-cluster result in the dynamic slowdown when system’s temperature closes to the transition temperature Tg. Keywords: Molecular dynamics simulation, coordination cell, liquid silica, dynamic slowdown.