ĐẠI HỌC THÁI NGUYÊN
TRƯỜNG ĐẠI HỌC SƯ PHẠM
PHẠM XUÂN TRƯỜNG
NGHIÊN CỨU CẤU TRÚC
VÀ CƠ CHẾ KHUẾCH TÁN TRONG SiO2 LỎNG
BẰNG PHƯƠNG PHÁP MÔ PHỎNG
LUẬN VĂN THẠC SĨ KHOA HỌC VẬT CHẤT
THÁI NGUYÊN - 2017
ĐẠI HỌC THÁI NGUYÊN
TRƯỜNG ĐẠI HỌC SƯ PHẠM
PHẠM XUÂN TRƯỜNG
NGHIÊN CỨU CẤU TRÚC
VÀ CƠ CHẾ KHUẾCH TÁN TRONG SiO2 LỎNG
BẰNG PHƯƠNG PHÁP MÔ PHỎNG
Chuyên ngành: Vật lí chất rắn
Mã số: 60.44.01.04
LUẬN VĂN THẠC SĨ KHOA HỌC VẬT CHẤT
Người hướng dẫn khoa học: TS. Phạm Hữu Kiên
THÁI NGUYÊN - 2017
LỜI CAM ĐOAN
Tôi xin cam đoan đây là đề tài của riêng tôi, do chính tôi thực hiện dưới
sự hướng dẫn của TS. Phạm Hữu Kiên và trên cơ sở nghiên cứu các tài liệu
tham khảo. Nó không trùng kết quả với bất kì tác giả nào từng công bố. Nếu sai
tôi xin chịu trách nhiệm trước hội đồng.
Thái Nguyên, ngày tháng 4 năm 2017
Tác giả luận văn
Phạm Xuân Trường
i
LỜI CẢM ƠN
Tôi xin bày tỏ lòng biết ơn sâu sắc đến TS Phạm Hữu Kiên, thầy đã nhiệt
tình hướng dẫn và giúp đỡ cho tôi trong quá trình hoàn thành luận văn.
Tôi xin chân thành cảm ơn các thầy, cô giáo Khoa Vật lý, Trường Đại
học Sư phạm Thái Nguyên đã tạo điều kiện cho tôi được học tập và làm việc
trong quá trình thực hiện luận văn.
Xin chân thành cảm ơn Phòng đào tạo, Trường Đại học Sư phạm Thái
Nguyên đã tạo điều kiện thuận lợi cho tôi trong quá trình thực hiện luận văn.
Cuối cùng tôi xin cảm ơn gia đình, đồng nghiệp, anh chị em lớp Cao học
Vật lý chất rắn K23 đã dành nhiều tình cảm, động viên, giúp đỡ tôi vượt qua
những khó khăn để hoàn thành luận văn.
Thái Nguyên, ngày tháng 4 năm 2017
Tác giả luận văn
Phạm Xuân Trường
ii
MỤC LỤC
LỜI CAM ĐOAN ................................................................................................ i
LỜI CẢM ƠN .................................................................................................... ii
MỤC LỤC ......................................................................................................... iii
DANH MỤC CÁC KÝ HIỆU VÀ CHỮ VIẾT TẮT ..................................... iv
DANH MỤC BẢNG .......................................................................................... v
DANH MỤC HÌNH .......................................................................................... vi
MỞ ĐẦU ............................................................................................................. 1
1. Lý do chọn đề tài ......................................................................................... 1
3. Phương pháp nghiên cứu của đề tài ............................................................. 2
4. Đối tượng và phạm vi nghiên cứu của đề tài ............................................... 2
5. Dự kiến đóng góp của đề tài ........................................................................ 2
6. Cấu trúc của đề tài ....................................................................................... 2
Chương 1 TỔNG QUAN ................................................................................... 4
1.1. Cấu trúc và các tính chất về vật liệu SiO2 ở trạng thái lỏng..................... 4
1.2. Một số phương pháp mô phỏng ................................................................ 7
1.2.1.Tổng quan về các phương pháp mô phỏng ........................................ 7
1.2.2. Các phương pháp mô phỏng .............................................................. 9
1.3. Mô phỏng cơ chế khuếch tán .................................................................. 12
1.3.1. Các định luật khuếch tán ................................................................. 12
1.3.2. Cơ chế khuếch tán ........................................................................... 13
Chương 2 PHƯƠNG PHÁP TÍNH TOÁN.................................................... 20
2.1. Phương pháp động lực học phân tử ........................................................ 20
2.2. Thế tương tác .......................................................................................... 24
2.3. Gần đúng Ewald-Hansen ........................................................................ 25
2.4. Xác định các đặc trưng vi cấu trúc và tính chất của mô hình ................ 26
2.4.1. Hàm phân bố xuyên tâm .................................................................. 27
2.4.2. Xác định số phối trí và độ dài liên kết ............................................. 30
iii
2.4.3. Xác định phân bố góc ...................................................................... 31
2.4.4. Xác định cơ chế khuếch tán ............................................................. 32
Chương 3 KẾT QUẢ VÀ THẢO LUẬN ....................................................... 35
3.1. Khảo sát cấu trúc của SiO2 lỏng theo áp suất ......................................... 35
3.2. Khảo sát tính đa thù hình của SiO2 lỏng ................................................ 41
3.3. Khảo sát cơ chế khuếch tán trong SiO2 lỏng theo nhiệt độ .................... 46
KẾT LUẬN ....................................................................................................... 57
CÔNG TRÌNH ĐÃ CÔNG BỐ LIÊN QUAN ĐẾN LUẬN VĂN ............... 58
TÀ I LIỆU THAM KHẢ O ............................................................................... 59
iv
DANH MỤC CÁC KÝ HIỆU VÀ CHỮ VIẾT TẮT
VĐH : Vô định hình
HPBXT : Hàm phân bố xuyên tâm
ĐLHPT : Động lực học phân tử
SPTTB : Số phối trí trung bình
SPT : Số phối trí
BKS : van Beest, Kramer and van Santen
PBGLK : Phân bố góc liên kết
iv
DANH MỤC BẢNG
Bảng 2.1. Các hệ số của thế BKS đối với hệ SiO2 ........................................ 25
Bảng 3.1. Các đặc tính cấu trúc của SiO2 lỏng. .............................................. 37
Bảng 3.2. Phân bố liên kết cầu O giữa hai đơn vị cấu trúc SiOx với mlà số
nguyên tử O tham gia liên kết cầu giữa hai đơn vị cấu trúc SiOx lân
cận. Các cột tiếp theo chỉ ra tỷ lệ phần trăm liên kết cầu tương ứng
với m. Ví dụ, 16,20% số liên kết giữa hai đơn vị cấu trúc lân cận có
hai nguyên tử O tham gia cầu liên kết ở áp suất nén 9,83 GPa. ..... 41
Bảng 3.3. Tỷ lệ các phản ứng SiOx SiOx' và OSiy OSiy' trong SiO2 lỏng theo
nhiệt độ 2600K, 3000K, 3200K và 3500K. .................................... 48
v
DANH MỤC HÌNH
Hình 1.1. Sơ đồ nghiên cứu bằng phương pháp mô phỏng ............................... 8
Hình 1.2. Cơ chế khuếch tán xen kẽ................................................................ 14
Hình 1.3. Cơ chế khuếch tán qua nút khuyết .................................................. 16
Hình 1.4. Cơ chế khuếch tán tập thể ............................................................... 17
Hình 2.1. Mô hình tính toán gần đúng Ewald –Hansen trong không gian 2 chiều,
mạng tuần hoàn 3x3 được dựng lên từ ô cơ sở có tâm n (0,0). .............. 26
Hình 2.2. Mô hình hóa các loại phản ứng trong SiO2 lỏng ............................. 33
Hình 3.1. Hàm phân bố xuyên tâm thành phần của SiO2 lỏng ở các áp suất
khác nhau và T=3200K. .................................................................. 36
Hình 3.3. Các đơn vị cấu trúc ................................................................................. 39
Hình 3.2. Mạng cấu trúc của mẫu SiO2 lỏng ở 3200 K................................... 38
Hình 3.4. Sự phụ thuộc của tỷ lệ các đơn vị cấu trúc SiO4, SiO5 và SiO6 vào
áp suất của mẫu SiO2 ở nhiệt độ 3200K. ......................................... 40
Hình 3.5. Phân bố khoảng cách liên kết trong: a) SiO4, b) SiO5 và c) SiO6. .. 42
Hình 3.6. Phân bố góc liên kết O-Si-O trong các đơn vị cấu trúc: ................. 43
Hình 3.7. Phân bố góc liên kết Si-O-Si giữa các đơn vị cấu trúc ................... 44
Hình 3.8. Phân bố góc Si-O-Si trong các đơn vị cấu trúc ............................... 44
Hình 3.9. Sự phụ thuộc của mật độ vào áp suất của mẫu SiO2 lỏng ở nhiệt độ
3200K .............................................................................................. 46
Hình 3.10. Hình minh họa các đơn vị cấu trúc ................................................. 47
Hình 3.11. Sự phụ thuộc của các loại phản ứng khác nhau theo bước ĐLHPT.
Nhiệt độ là 2600 K. ......................................................................... 49
Hình 3.12. Số lượng đám và số lượng nguyên tử trong đám lớn nhất đối với hai
trường hợp: 1/ Các ô phối trí được chọn ngẫu nhiên trong cấu trúc
mạng; 2/ Các ô sai hỏng phát hiện trong hệ. Nhiệt độ là 2600 K. .. 49
Hình 3.13. Sự phụ thuộc số lượng các nguyên tử trong đám lớn nhất theo bước
ĐLHPT. Nhiệt độ là 2600 K. .......................................................... 50
vi
Hình 3.14. Số đám không phản ứng và số lượng nguyên tử trong đám lớn nhất
cho các khoảng thời gian khác nhau. Nhiệt độ là 2600 K. ............. 52
Hình 3.15. Số đám phản ứng và số lượng nguyên tử trong đám lớn nhất theo
các khoảng thời gian khác nhau. Nhiệt độ là 2600 K. .................... 52
Hình 3.16. Số lượng nguyên tử trong đám lớn nhất và khoảng cách dịch chuyển
bình phương trung bình theo khoảng thời gian khác nhau, mỗi
khoảng thời gian là 2000 bước. Nhiệt độ là 2600 K ....................... 53
Hình 3.17. Số lượng nguyên tử trong đám lớn nhất và khoảng cách dịch chuyển
bình phương trung bình theo khoảng thời gian khác nhau, mỗi
khoảng thời gian là 5000 bước. Nhiệt độ là 2600 K. ...................... 53
Hình 3.18. Số lượng nguyên tử phản ứng và số nguyên tử trong đám không
phản ứng lớn nhất. ........................................................................... 55
Hình 3.19. Hình minh họa sự phản ứng cho các trường hợp ............................ 56
vii
MỞ ĐẦU
1. Lý do chọn đề tài
Hiện nay, nghiên cứu cấu trúc và động học trong các chất lỏng SiO2,
MgO và Al2O3 có cấu trúc mạng đã nhận được sự quan tâm đặc biệt của các
nhà khoa học. Khi thay đổi áp suất và nhiệt độ của các chất lỏng có cấu trúc
mạng thể hiện nhiều tính chất lý thú chẳng hạn như: tính đa thù hình và tính
không đồng nhất động học [1, tr.7]. Hiện tượng không đồng nhất được quan
sát trực tiếp trong chất keo, chất lỏng được làm nguội nhanh bằng phương
pháp thực nghiệm. Các hiện tượng không đồng nhất cấu trúc có thể được giải
thích bằng phương pháp mô phỏng và lý thuyết. Trong đó, ôxít silic (SiO2) là
một trong những đối tượng được nhiều nhà khoa học tập trung nghiên cứu
[18, tr.12].
Vật liệu SiO2 có nhiều ứng dụng trong việc chế tạo các loại linh kiện
và vật liệu. Sự hiểu biết về cấu trúc, các tính chất vật lý đặc trưng và cơ chế
động học ở mức độ nguyên tử của loại vật liệu này dưới ảnh hưởng áp suất
và nhiệt độ là rất cần thiết. Tuy nhiên, hiểu biết chi tiết về cấu trúc vi mô
của vật liệu ôxít SiO2 cũng như cơ chế khuếch tán, mối liên hệ giữa các đặc
trưng cấu trúc và một số hiện tượng động học trong các vật liệu ô xít SiO2
vẫn còn bỏ ngỏ. Đặc biệt, sự thay đổi cấu trúc của SiO2 ở nhiệt độ 3200 K
trong dải áp suất 0-25 GPa vẫn đang là nội dung hấp dẫn. Do đó, chúng tôi
lựa chọn đề tài “Nghiên cứu cấu trúc và cơ chế khuếch tán trong SiO2
lỏng bằng phương pháp mô phỏng” để bổ sung thêm những hiểu biết về
cấu trúc, tính chất vật lí cũng như các tính chất khuếch tán trong SiO2 lỏng.
2. Mục đích của đề tài
Xây dựng mẫu ôxít SiO2 lỏng với kích thước 1998 nguyên tử (666 Si
và 1332 O) ở nhiệt độ 3200 K. Các đặc trưng cấu trúc của các mẫu vật liệu
được khảo sát thông qua hàm phân bố xuyên tâm (HPBXT) thành phần, phân
bố số phối trí (SPT) và phân bố góc liên kết (PBGLK) trong giải áp suất từ 0
đến 25GPa.
1
Khảo sát và làm sáng tỏ tính đa thù hình của ôxít SiO2 ở các áp suất khác
nhau từ 0-25GPa. Từ đó xác định được mật độ của mẫu SiO2 lỏng ở một trạng
thái nào đó thông qua tỷ lệ các đơn vị cấu trúc SiOx.
Khảo sát cơ chế khuếch tán của nguyên tử O và Si theo nhiệt độ theo cơ
chế đám.
3. Phương pháp nghiên cứu của đề tài
Luận văn sử dụng phương pháp mô phỏng động lực học phân tử, phương
pháp phân tích vi cấu trúc để xây dựng, phân tích và tính toán các đặc trưng cấu
trúc và cơ chế khuếch tán của vật liệu SiO2 lỏng.
4. Đối tượng và phạm vi nghiên cứu của đề tài
Mẫu vật liệu SiO2 lỏng ở nhiệt độ 3200 K trong giải áp suất 0-25 GPa để
nghiên cứu cấu trúc và tính đa thù hình.
Mẫu vật liệu SiO2 lỏng có nhiệt độ lần lượt là 2600 K, 3000 K, 3200 K
và 3500 K để nghiên cứu cơ chế khuếch tán theo nhiệt độ.
5. Dự kiến đóng góp của đề tài
Kết quả của luận văn có những đóng góp:
i/ Cho thấy các thông tin về đặc trưng vi cấu trúc và tính đa thù hình của
ôxít SiO2 lỏng ở nhiệt độ 3200 K, trong dải áp suất 0-25 GPa.
ii/ Cung cấp những hiểu biết về cơ chế khuếch tán theo nhiệt độ cho các
nhà nghiên cứu.
6. Cấu trúc của đề tài
Luận văn gồm phần mở đầu và kết luận, nội dung của luận văn gồm 3
chương, trong đó:
Chương 1, trình bày tổng quan về tình hình nghiên cứu cấu trúc và tính
chất của vật liệu SiO2. Sau đó trình bày tổng quan về một số phương pháp mô
phỏng và mô phỏng cơ chế khuếch tán.
2
Chương 2, trình bày kỹ thuật mô phỏng, phương pháp tính toán hệ ôxít
SiO2 ở trạng lỏng. Trước tiên, chúng tôi trình bày phương pháp mô phỏng động
lực học phân tử, thế tương tác cặp BKS, gần đúng Ewald-Hansen. Tiếp theo
chúng tôi trình bày các kỹ thuật phân tích cấu trúc như: Hàm phân bố xuyên
tâm, phân bố số phối trí, phân bố góc liên kết. Cuối cùng, chúng tôi trình bày
cơ chế khuếch tán thông qua độ dịch chuyển bình phương trung bình.
Chương 3, trình bày về các kết quả nghiên cứu đặc trưng cấu trúc, đa thù
hình và cơ chế khuếch tán của nguyên tử Si và O trong SiO2 lỏng.
3
Chương 1
TỔNG QUAN
Trong chương này, chúng tôi trình bày: tổng quan một số kết quả nghiên
cứu cấu trúc và các tính chất về vật liệu SiO2 ở trạng thái lỏng bằng cả thực
nghiệm, lý thuyết lẫn mô phỏng máy tính. Một số cơ chế khuếch tán trong vật
liệu cũng được đề cập trong chương này.
1.1. Cấu trúc và các tính chất về vật liệu SiO2 ở trạng thái lỏng
Trong nhưng năm gần đây, nhiều nghiên cứu thực nghiệm và lý thuyết
trong lĩnh vực khoa học vật liệu đã chứng minh rằng cách xắp xếp nguyên tử và
quá trình hình thành cấu trúc có ảnh hưởng đến các tính chất vật lý của SiO2.
Tuy nhiên, nghiên cứu sự chuyển pha cấu trúc của SiO2 rất khó đo được bằng
thực nghiệm. Sự hiểu biết biết về hiện tượng này vẫn chưa được thỏa đáng và
vẫn đang còn nhiều vấn đề cần thảo luận, đặc biệt là sự thay đổi cấu trúc ở
nhiệt độ 3200K trong một dải áp suất vẫn đang là một đề tài nóng.
Trong lĩnh vực khoa học trái đất, SiO2 lỏng trong lòng trái đất lại tồn tại
ở nhịêt độ 3200K nên được rất nhiều nhà khoa học tập trung nghiên cứu bằng
nhiều kĩ thuật thực nghiệm và tính toán lý thuyết như nhiễu xạ tia X, nhiễu xạ
nơtron, phổ Raman, phổ hấp thụ tia X, cộng hưởng từ hạt nhân. Còn đối với
phương pháp mô phỏng, việc nghiên cứu sẽ đem lại nhiều kết quả có giá trị
giúp so sánh dự báo các tính chất mới lạ của vật liệu. Trong công trình [23]
trạng thái của SiO2 lỏng bao gồm các đơn vị cấu trúc cơ bản liên kết nhau trong
một mạng liên tục trong không gian ba chiều hữu hạn và không có trật tự xa.
Mỗi đơn vị cấu trúc cơ bản là một khối tứ diện SiO4 với 4 nguyên tử ôxi(O)
nằm ở đỉnh và tâm là nguyên tử silic(Si). Liên kết giữa hai tứ diện đòi hỏi một
liên kết góc Si-O-Si và hai góc nhị diện. Sự biến đổi hai góc này được xem như
nguyên nhân chính làm thay đổi cấu trúc của SiO2 lỏng. Trong công trình [20]
bằng kĩ thuật nhiễu xạ tia X đã cho kết quả: trong mỗi đơn vị cấu trúc SiO4 các
4
thông tin cấu trúc được xác định bởi: số phối trí trung bình (SPTTB); độ dài
liên kết Si-Si, O-O, Si-O; góc liên kết Si-O-Si và góc liên kết O-Si-O. Trong
thực nghiệm từ đường cong tán xạ tia X hoặc tán xạ nơtron ta có thể xác định
được thông số quan trọng xác định cấu trúc vật liệu. Thừa số cấu trúc cho phép
xác định số lượng trung bình các nguyên tử ở khoảng cách bất kì tính từ
nguyên tử đang xét. Khi phân tích Phu-ri-ê thừa số cấu trúc ta còn thu được
HPBXT, một thông số dùng để xác định trật tự gần của các vật liệu có cấu trúc
mất trật tự.
Thực nghiệm đã chứng tỏ rằng cấu trúc của SiO2 lỏng phụ thuộc khá
mạnh vào áp suất và ít thay đổi theo nhiệt độ. Phân tích thừa số cấu trúc nhiễu
xạ tia X cho thấy khi tăng áp suất có sự thay đổi vị trí đỉnh và cường độ đỉnh
nhiễu xạ thứ nhất. Khi áp suất tăng từ 0,1MPa đến 8GPa thì vị trí đỉnh nhiễu xạ
thứ nhất dịch từ Q~ 1,55 đến 1,92 Å-1 trong khi cường độ của nó hấu như
không thay đổi. Sự thay đổi của thừa số cấu trúc xảy ra mạnh nhất trong vùng
áp suất từ 8 GPa đến 28 GPa, cường độ đỉnh nhiễu xạ thứ nhất giảm đi 50%
trong khi vị trí đỉnh nhiễu xạ thứ nhất dịch từ Q~ 1,92 đến 2,29 Å-1 và xuất hiện
thêm đỉnh nhiễu xạ mới ở vị trí 3,18 Å-1. Nếu tiếp tục tăng áp suất đến 42 GPa
thì thừa số cấu trúc gần như không thay đổi nữa. Tiến hành phân tích chuỗi
Phu-ri-ê thừa số cấu trúc thu được HPBXT g(r), từ đó xác định độ dài liên kết
trung bình Si-O, O-O, Si-Si trong SiO2 ở áp suất thường tương ứng bằng 1,59;
2,61; và 3,07 với sai số 0,01 Å-1. Phân tích các giá trị đỉnh độ dài liên kết Si-O
và O-O còn cho thấy góc liên kết O-Si-O bên trong đa diện có giá trị khoảng
960 và ở áp suất 42 GPa. Giá trị này nằm giữa hai giá trị 109,50 và 900 tương
ứng với cấu trúc tứ diện và bát diện. Khi tăng áp suất, SPTTB của Si-O chuyển
dần từ 4 đến 6 [6]. Nhiều công trình nghiên cứu cấu trúc SiO2 về PBGLK Si-O-
Si. công trình Mozzi và Warre [27], công bố kết quả xác xuất PBGLK của Si-
O-Si trong thủy tinh cho kết quả gần đúng với giá trị 1440 với độ rộng ở vị trí
nửa cực đại khoảng 360 và bị lệch về phía góc nhỏ hơn. Sau đó, năm 1995,
5
Poulsen và các cộng sự đã tiến hành phân tích nhiễu xạ tia X năng lượng cao
đối với SiO2 VĐH đã thu được PBGLK của Si-O-Si rất gần với Mozzi và
Warre tương ứng với 1470 và 350[27]. Gần đây, năm 2008, Trong công trình
[21] đã phân tích kết quả phổ O17 NMR của SiO2 thủy tinh với góc liên kết Si-
O-Si trung bình là 1500 với độ rộng phổ ở vị trí cực đại rất hẹp khoảng 160. Sau
đó, trong công trình [25] thực hiện thay đổi cấu trúc hình học của vật liệu SiO2
dưới áp suất khác nhau. kết quả cho thấy độ cao của đỉnh hàm phân bố cặp
giảm đi khi áp suất nén tăng lên vượt ngưỡng 15GPa và điều này chỉ ra các cấu
trúc đơn vị tứ diện co lại dưới tác dụng của áp suất. Khi áp suất vượt qua
ngưỡng trên, chiều dài liên kết trung bình Si-O và số phối trí tăng tuyến tính
với áp suất.
Cấu trúc của vật liệu SiO2 lỏng cũng được nghiên cứu bằng các phương
pháp mô phỏng khác nhau nhưng phương pháp ĐLHPT được sử dụng rộng
rãi nhất. Mô hình SiO2 được xây dựng đầu tiên bằng phương pháp động lực
học sử dụng thế Born-Mayer. Mô hình gồm 162 nguyên tử chứa trong hình
hộp lập phương mô phỏng và sử dụng điều kiện biên tuần hoàn [6]. Gần
đúng Eward đã sử dụng để tính tương tác Cu-Lông ở khoảng cách xa. Bằng
phương pháp này, người ta đã tạo ra mô hình SiO2 thủy tinh ở nhiệt độ
300K. Mô hình thu được HPBXT khá phù hợp với thực nghiệm. Sau đó, mô
hình SiO2 ở nhiệt độ 1500K gồm 375 nguyên tử được xây dựng bằng
phương pháp ĐLHPT. Kết quả thu được là trạng thái thủy tinh có cấu trúc
tứ diện với các đỉnh là nguyên tử O. HPBXT thành phần từ mô hình này
cho kết quả tương đối phù hợp với thực nghiệm nhiễu xạ tia X ở vị trí các
đỉnh thứ nhất và thứ hai, PBGLK của O-Si-O tính toán được có vị trí
khoảng 109,50 và góc Si-O-Si có đỉnh ở vị trí khoảng 1510. Kết quả còn chỉ
ra rằng, SiO2 có cấu trúc xốp và chứa nhiều quả cầu lỗ hổng [31]. Ở áp suất
35 GPa, nhiệt độ 6000 K mô hình SiO2 được nghiên cứu bằng phương pháp
ĐLHPT cho thấy số phối trí tăng lên theo áp suất.
6
Gần đây, bằng phương pháp mô phỏng ĐLHPT cho SiO2 lỏng, P.K.
Hung và các cộng sự đã chỉ ra rằng ở trạng thái lỏng, SiO2 được cấu tạo từ các
đơn vị cấu trúc SiO4, SiO5 và SiO6. Thể tích của không gian mô phỏng SiO2 là
một hàm tuyến tính của nồng độ các đơn vị cấu trúc đó. Khi nén SiO2 lỏng ở
các áp suất khác nhau thì mật độ của SiO2 thay đổi là do sự thay đổi của nồng
độ các đơn vị cấu trúc SiO4, SiO5và SiO6. Tuy nhiên, khi áp suất thay đổi thì độ
dài liên kết trung bình Si–O, O–O và Si–Si và góc liên kết trung bình Si–O–Si
trong mỗi đơn vị cấu trúc hầu như không thay đổi. Kết quả mô phỏng ở nhiệt
độ 6000 K và áp suất 20 GPa cho kết quả độ dài liên kết trung bình Si–O, O–O
và Si–Si tương ứng bằng 1,60; 2,54 và 3,12 Å với sai số 0,02 Å và góc liên kết
trung bình Si–O–Si khoảng 145 Å [32, 10]. Sự ảnh hưởng của kích thước mô
hình lên các PBGLK của Si–O–Si và O–Si–O, cũng như là độ dài các liên kết
Si–O, O–O và Si–Si được nghiên cứu bằng mô hình ĐLHPT [26]. Kết quả
nghiên cứu chỉ ra rằng không thấy bất cứ hiệu ứng kích thước lên tính chất vật
lý và hóa học của vật liệu SiO2 lỏng.
Cho tới nay đã có một lượng lớn công trình nghiên cứu về sự chuyển pha
cấu trúc của SiO2 [1, tr.12]. Tuy nhiên, sự hiểu biết đầy đủ về hiện tượng này
vẫn chưa thoả đáng và còn nhiều vấn đề đang được thảo luận. Đặc biệt, sự thay
đổi cấu trúc ở nhiệt độ 3200 K trong một dải áp suất vẫn đang là đề tài nóng (vì
đây là nhiệt độ trong lòng Trái đất nơi tồn tại SiO2). Vì vậy, trong nghiên cứu
này chúng tôi muốn cung cấp thêm một vài thông tin về vi cấu trúc cũng như sự
chuyển pha cấu trúc trong vật liệu SiO2 lỏng khi nén mô hình.
1.2. Một số phương pháp mô phỏng
1.2.1.Tổng quan về các phương pháp mô phỏng
Mô phỏng là việc nghiên cứu trạng thái của mô hình để qua đó hiểu được
hệ thống thực, mô phỏng là tiến hành thử nghiệm trên mô hình. Đó là quá trình
tiến hành nghiên cứu trên vật nhân tạo, tái tạo hiện tượng mà người nghiên cứu
cần để quan sát và làm thực nghiệm, từ đó rút ra kết luận tương tự vật thật.
7
Ta có thể thực hiện việc mô phỏng từ những phương tiện đơn giản như
giấy, bút đến các nguyên vật liệu tái tạo lại nguyên mẫu (mô hình bằng gỗ,
gạch, sắt…) hay hiện đại hơn là dùng máy tính điện tử.
Tất cả các kỹ thuật sử dụng máy tính để nghiên cứu, khảo sát các đối
tượng, quá trình vật lý xảy ra được gọi là mô phỏng hay mô hình hóa trong vật
lý. Các đối tượng và các quá trình mà chúng ta quan tâm được gọi là các hệ vật
lý. Khi mô phỏng chúng ta phải xây dựng một tập hợp các giả thiết để mô tả
hoạt động của hệ thống. Các giả thiết này bao gồm các mối quan hệ logic, các
công thức toán học. Chúng cho phép xây dựng nên các mô hình trợ giúp cho
việc khảo sát hệ thống và các quá trình vật lý xảy ra trên nó. Nếu mô hình phức
tạp chúng ta giải quyết vấn đề với sự trợ giúp của thí nghiệm số hay phương
pháp mô phỏng.
Quá trình nghiên cứu bằng phương pháp mô phỏng được thể hiện trong
sơ đồ sau:
Hình 1.1. Sơ đồ nghiên cứu bằng phương pháp mô phỏng
Nhìn vào hình 1.1 thấy rằng để nghiên cứu hệ thống thực ta phải tiến
hành mô hình hóa tức là xây dựng mô hình mô phỏng. Khi có mô hình mô
phỏng sẽ tiến hành làm các thực nghiệm trên mô hình để thu được các kết quả
mô phỏng. Thông thường kết quả mô phỏng có tính trừu tượng của toán học
nên phải thông qua xử lý mới thu được các thông tin kết luận về hệ thống thực.
Sau đó dùng các thông tin và kết luận trên để hiệu chỉnh hệ thực theo mục đích
nghiên cứu đã đề ra.
8
Các dạng mô phỏng bao gồm: Mô phỏng động (thời gian đóng vai trò
quan trọng đối với thực nghiệm mô phỏng); Mô phỏng tĩnh (không có biến thời
gian); Mô phỏng xác định (các sự kiện xảy ra trong thực nghiệm mô phỏng
theo một quy luật xác định, chính xác không có yếu tố ngẫu nhiên); Mô phỏng
ngẫu nhiên (có yếu tố ngẫu nhiên); Mô phỏng liên tục (các sự kiện xảy ra trong
thời gian liên tục); Mô phỏng gián đoạn (số lượng thời gian xác định).
1.2.2. Các phương pháp mô phỏng
Mô hình hoá là sự gắn kết giữa lí thuyết và thực nghiệm. Các vật liệu
ôxít đã và đang được nghiên cứu bằng cả phương pháp thực nghiệm và mô
phỏng. Mô phỏng cho phép xây dựng mẫu vật liệu ở dạng mô hình và khảo sát
tính chất vật lí của chúng. Thực chất là quá trình mô phỏng lại kết quả nghiên
cứu ở phòng thí nghiệm, còn được gọi là phương pháp thực nghiệm mô hình
hay phương pháp thực nghiệm máy tính.
Gần đây phương pháp mô phỏng đã cho một số lượng lớn các công trình
nghiên cứu và được ứng dụng tốt. Các phương pháp mô phỏng thường được sử
dụng như: Nguyên lí ban đầu, Monte-Carlo, liên kết chặt và động lực học phân
tử. Trong đó:
Phương pháp nguyên lí ban đầu dựa trên việc giải hệ phương trình
Schrodinger cho hệ nhiều điện tử và không sử dụng bất cứ thông số thực
nghiệm nào. Phương pháp này có nhiều ưu điểm và được ứng dụng rộng rãi tuy
nhiên mặt hạn chế là áp dụng cho các hệ nhỏ từ vài chục đến vài trăm nguyên
tử. Đối với phương pháp Monte-Carlo tính toán chuyển đổi cùng một vị trí của
nguyên tử tuân theo thống kê Boltzmann. Phương pháp liên kết chặt tính toán
Hamintonien và các ma trận cơ sở dựa trên một số dữ liệu thực nghiệm và xét
đến ảnh hưởng của các hiệu ứng lượng tử. Phương pháp này có thể áp dụng cho
những hệ lớn nhiều nguyên tử, được sử dụng đề tài nghiên cứu cấu trúc điện tử.
Đối với phương pháp ĐLHPT các tính toán được thực hiện trên cơ sở phương
trình chuyển động Newton cho các nguyên tử. Phương pháp này cho phép theo
9
dõi chuyển động của một tập hợp các nguyên tử theo thời gian và có thể xác
định ảnh hưởng của nhiệt độ, áp suất đến các tính chất hoá lý của chúng. Một
số tính chất vật lý như cấu trúc địa phương, các tính chất nhiệt động, tính chất
khuếch tán…có thể được khảo sát bằng phương pháp ĐLHPT. Phương pháp
ĐLHPT cổ điển với thế tương tác cặp có thế mạnh khi mô tả vi cấu trúc nhưng
nó không thể mô tả đúng đắn tất cả các tính chất vật lý của vật liệu. Điều này
liên quan trực tiếp đến mức độ tin cậy của mô hình ĐLHPT cổ điển và cũng có
thể nhận thấy các mô hình với thế tương tác khác nhau sẽ cho các số liệu khác
nhau. Tuy nhiên giá trị của các mô hình này là dự báo nhiều hiện tượng thú vị,
có tính chất định hướng và dẫn đến nhiều nghiên cứu bằng các phương pháp
khác, chẳng hạn khả năng thay đổi cơ chế khuyếch tán khi chuyển từ pha mật
độ thấp sang pha mật độ cao.
Gần đây, một số tác giả sử dụng đồng thời cả hai phương pháp thực
nghiệm và mô phỏng cho nghiên cứu của mình [5, 7,27], ví dụ công trình [9] đã
nghiên cứu vi cấu trúc của Al2O3 bằng phương pháp cả nhiễu xạ tia X, nhiễu xạ
nơtron và Monte-Carlo đảo. Một số tác giả khác kết hợp các ưu điểm của nhiều
phương pháp mô phỏng với nhau để thực hiện nghiên cứu, ví dụ trong công
trình [29] các tác giả đã kết hợp phương pháp nguyên lý ban đầu và ĐLHPT để
nghiên cứu tính chất cấu trúc và cơ chế khuếch tán Hiđrô trong SiO2 nóng chảy.
Khi thực hiện nghiên cứu mô phỏng những vấn đề chính ảnh hưởng đến
độ tin cậy của kết quả thu được là: cách chọn thế tương tác, điều kiện biên, kích
thước mô hình .
Thứ nhất về thế tương tác: Xét về khía cạnh vật lí, thế tương tác giữa các
nguyên tử được xác định bởi tương tác giữa các nguyên tử. Thế tương tác được
sử dụng trong mô phỏng là thế van Beest, Kramer and van Santen (BKS), được
trình bày rõ ở chương 2.
10
Thứ hai về cách chọn điều kiện biên: Có các loại điều kiện biên chính
như biên tự do, biên cứng, biên mềm, biên tuần hoàn.
Biên tự do: là khoảng cách chân không tuyệt đối bao quanh không gian
tính toán, thường dùng mô phỏng các hiện tượng bề mặt. Sử dụng biên tự do
tuy đơn giản song kém chính xác.
Biên cứng: là lớp nguyên tử biên đứng yên bao bọc quanh không gian
tính toán mà chiều dày của lớp này lớn hơn khoảng tương tác giữa các nguyên
tử, các nguyên tử trong lớp bao bọc này có thể tương tác với các nguyên tử
trong không gian tính toán. Biên cứng thường sử dụng để nghiên cứu các
khuyết tật điểm.
Biên mềm: gồm hai lớp nguyên tử biên xung quanh không gian tính toán.
Lớp ngoài tương tự như biên cứng, lớp thứ hai dịch chuyển theo biên độ cho
trước và tắt dần từ trong ra ngoài, có tác dụng giảm hiệu ứng bề mặt tốt hơn
biên cứng, thường được sử dụng cho mô hình hoá các khuyết tật kéo dài.
Biên tuần hoàn: là không gian tương tác lặp lại tuần hoàn trong không
gian. Ở điều kiện tuần hoàn các nguyên tử trong không gian tính toán của cực
phải tương tác với cực trái, ở đỉnh tương tác với đáy, phía trước tương tác với
phía sau. Đường kính không gian tính toán lớn hơn hai lần khoảng cách tương
tác giữa hai nguyên tử riêng biệt, không tồn tại các nguyên tử bề mặt. Biên tuần
hoàn được sử dụng để mô hình hoá những hệ lớn.
Ngày nay, với sự trợ giúp của kĩ thuật tính toán hiện đại (tính toán song
song, tính toán phân tán) cùng với sự ra đời của các loại máy tính có tốc độ cao
thì kích thước mô hình vật liệu đã được tăng lên đáng kể. Nếu như trước đây
mô hình tính toán cho vật liệu vào khoảng vài trăm nguyên tử thì ngày nay lên
đến hàng triệu nguyên tử. Tuy nhiên nghiên cứu vi mô của một số ôxít về cơ
chế khuyếch tán, tính đa thù hình, chuyển pha và sự ảnh hưởng của áp suất,
nhiệt độ đến các tính chất hệ khi khảo sát kích thước mô hình tính toán với 500;
2000; 3000 là gần như nhau.
11
Trên những cơ sở đó luận văn này tôi chọn phương pháp ĐLHPT với
điều kiện biên tuần hoàn, thế tương tác BKS.
1.3. Mô phỏng cơ chế khuếch tán
Khuếch tán là sự di chuyển của các phân tử, nguyên tử hay ion từ vị trí
này đến vị trí khác trong vật liệu. Khuếch tán là quá trình phụ thuộc thời gian.
Tốc độ khuếch tán được đánh giá qua khái niệm dòng khuếch tán.
Dòng khuếch tán j là lượng chất M khuếch tán qua đơn vị diện tích của
chất rắn vuông góc với phương khuếch tán trong 1 đơn vị thời gian.
1.3.1. Các định luật khuếch tán
Các quy luật khuếch tán được thể hiện qua định luật Fick, được đề xuất
bởi Adolf Fick vào năm 1855. Ngoài ra, nhiều trường hợp khuếch tán tuân theo
quy luật Arrhenius đơn giản hơn.
+ Định luật Fick 1
Xem xét dòng các hạt khuếch tán trong trường hợp 1 chiều. Các hạt có
thể là các nguyên tử, phân tử hoặc ion. Định luật Fick thứ nhất đối với môi
trường đẳng hướng có thể viết như sau:
(1.1)
ở đây Jx là dòng hạt khuếch tán và C là nồng độ. Dấu trừ trong phương
trình (1.1) cho thấy dòng khuếch tán theo chiều giảm của nồng độ.
Tổng quát trong trường hợp không gian ba chiều, định luật Fick có thể
viết như sau:
(1.2)
trong đó D là hệ số khuếch tán, có đơn vị là cm2 s-1.
+ Định luật Fick 2
Nếu nồng độ C không những là hàm của x mà còn phụ thuộc vào thời
gian t thì người ta thường sử dụng định luật Fick 2. Phương trình của định luật
Fick 2 như sau:
12
(1.3)
Trong trường hợp hệ số khuếch tán D không phụ thuộc vào nồng độ thì
phương trình Fick 2 có thể viết như sau:
(1.4)
Nghiệm của phương trình (1.4) trong trường hợp khuếch tán một chất
trên bề mặt vào bên trong mẫu với nồng độ ban đầu Co
có nồng độ Cs
(Cs>C0) có dạng:
(1.5)
trong đó erf(L) là hàm sai số của đại lượng L.
+ Quy luật Arrhenius
Trong nhiều trường hợp thì hệ số khuếch tán phụ thuộc vào nhiệt độ theo
quy luật Arrhenius:
(1.6)
Ở đây D0 là một hằng số, E là năng lượng kích hoạt, k là hằng số
Boltzomann.
1.3.2. Cơ chế khuếch tán
Khuếch tán nói chung, đặc biệt là khuếch tán trong vật liệu ô xít dạng
lỏng luôn là đề tài được nghiên cứu bằng cả lý thuyết và thực nghiệm. Tuy
nhiên, để hiểu rõ được cơ chế khuếch tán không phải là một vấn đề đơn giản.
Nhiều nghiên cứu đã đưa ra nhiều cách giải thích cho cơ chế khuếch tán
trong vật liệu ô xít như: cơ chế vacancy (cơ chế nút khuyết), cơ chế khuếch
tán tập thể, cơ chế khuếch tán theo các khe trống trong mạng tinh thể
(khuếch tán xen kẽ)…
13
+ Cơ chế khuếch tán xen kẽ
Các nguyên tử khuếch tán (có kích thước rất nhỏ hơn các nguyên tử nền)
chuyển động qua khe của các nguyên tử nền và định sứ ở lỗ hổng tạo ra bởi các
nguyên tử nền (hình 1.2).
Nguyên tử nền
Nguyên tử khuếch tán
Hình 1.2. Cơ chế khuếch tán xen kẽ
Phạm Khắc Hùng và đồng nghiệp đã nghiên cứu quá trình tự khuếch tán
trong mô hình hợp kim vô định hình một chiều ở OK bằng phương pháp thống
kê hồi phục tiếp diễn. Trong hệ cấu trúc chuẩn (canonical structures) sử dụng
thế tương tác cặp dạng:
(1.7)
với m = 6, 18, 30 (mô hình gồm 686 hạt chứa trong khối lập phương, r0
= 0,1 nm) và mô hình Fe VĐH (a - Fe) sử dụng thế tương tác Paka - Doyama.
Các tác giả đã khảo sát khuếch tán của các nguyên tử theo cơ chế điền kẽ, bởi
vì thế tương tác giữa các nguyên tử điền kẽ với các nguyên tử khác tương tự
như tương tác giữa các nguyên tử trong vật rắn vô định hình.
14
Trước hết xác định phân bố của năng lượng ở trạng thái vị trí và trạng
thái chuyển tiếp. Sau đó tính toán khoảng cách từ một điền kẽ tới các điền kẽ
lân cận khi thế năng tương tác đạt giá trị cực tiểu. Cuối cùng đã xác định hệ số
khuếch tán cho các nguyên tử điền kẽ bằng cách sử dụng phương pháp xác định
tốc độ dịch chuyển khỏi vị trí ban đầu dưới tác dụng của trường ngoài bằng
phương trình:
(1.8)
Kết quả xác định được sự phụ thuộc của hệ số tự khuếch tán có dạng của
phương trình Arrhenius:
(1.9)
Giá trị của hệ số B đối với cấu trúc vô định hình nằm giữa giá trị của cấu
trúc mạng lập phương tâm khối BCC và mạng lập phương tâm mặt FCC. Cơ
chế khuếch tán theo khe trống trong mạng tinh thể có thể được sử dụng khi
khảo sát sự khuếch tán trong hợp kim vô định hình. Các tác giả sử dụng
phương pháp hồi phục kích hoạt để mô phỏng cơ chế.
Roth đã tìm thấy kết quả cho nồng độ thấp, khi năng lượng E = 2,70V và
với thừa số trước hàm mũ là D0 = 2,2.104 cm2s-1. Giá trị cao của D0 không phù
hợp cho mẫu bẫy đơn. Điều này càng trở nên rắc rối nếu quá trình ủ đó dùng
cấu trúc hồi phục. Trong trường hợp này hàm phụ thuộc thời gian của hệ số
khuếch tán có thể bị phá vỡ.
+ Cơ chế khuếch tán Vacancy (Nút khuyết)
Các nguyên tử nhẩy từ vị trí cân bằng của nó sang vị trí nút khuyết lân
cận (hình 1.3).
Theo cơ chế khuếch tán vacancy, người ta cho rằng quá trình khuếch tán
của nguyên tử là sự trao đổi vị trí giữa các nguyên tử nằm tại nút mạng với các
vacancy bên cạnh và hệ số khuếch tán sẽ tỷ lệ với nồng độ cân bằng vacancy.
15
Nút khuyết (vacancy)
Hình 1.3. Cơ chế khuếch tán qua nút khuyết
Năng lượng kích hoạt bằng tổng năng lượng tạo vacancy và năng lượng
dịch chuyển vacancy.
Phạm Khắc Hùng và các tác giả đã chỉ rõ rằng một vacancy trong pha vô
định hình giống như một lỗ trống có thể lặp lại nhiều lần sự thay đổi vị trí với
các nguyên tử bên cạnh. Trong mẫu phân bố thông thường, chỉ có một lỗ trống
được tìm ra. Nếu áp dụng PPHPT sau khi nguyên tử và lỗ trống trao đổi vị trí,
hai hạt còn lại ở vị trí mới của chúng. Trong các trường hợp R < 80 pm sự đổi
chỗ của nguyên tử sau nhiều lần sẽ trở lại vị trí ban đầu. Trong kim loại vô định
hình, lỗ trống có bán kính R < 80 pm không đóng vai trò quan trọng. Xác suất
để một vacancy biến mất sau n bước . Số bước trung bình mà
một vacancy thực hiện trước khi biến mất là:
. (1.10)
Trong trường hợp ổn định nhất của mẫu ba chiều, một vacancy có thể
thực hiện trung bình 8 bước trước khi biến mất. Phạm Khắc Hùng đưa ra hình
thức khuếch tán trong vô định hình dẻo bằng cơ chế vacancy. Trước hết do kết
quả của dao động nhiệt, một lỗ trống xuất hiện. Sau đó lỗ trống di chuyển, trao
đổi vị trí với các nguyên tử bên cạnh và thực hiện trung bình 1/α bước nhảy.
Tiếp đó vacancy rơi vào một bẫy và trở thành một lỗ trống tương đối nhỏ hơn
và không thể tham gia vào quá trình thay đổi vị trí.
16
+ Cơ chế khuếch tán tập thể
Trong dung dịch thay thế, kích thước nguyên tử của chất tan tương tự
kích thước nguyên tử dung môi. Cơ chế khuếch tán của các nguyên tử chất tan
và các nguyên tử dung môi là sự đổi chỗ trực tiếp của các nguyên tử lân cận
hoặc có thể đổi chỗ theo đường vòng (hình 1.4).
Nguyên tử dung môi
Nguyên tử chất tan
Hình 1.4. Cơ chế khuếch tán tập thể
Đối với hệ phi tinh thể, người ta sử dụng phương pháp động lực học
phân tử để tính hệ số tự khuếch tán. Trong quá trình hồi phục của phương pháp
động lực học phân tử, thời gian mỗi bước nhảy thực hiện độ dịch chuyển nhỏ
có liên quan đến sự dịch chuyển của các nguyên tử lân cận. Điều này dẫn đến
sự xuất hiện cơ chế khuếch tán tập thể. Bình phương trung bình của độ dịch
chuyển của nguyên tử
Sự có mặt của số hạng thứ nhất là do chuyển động dao động của nguyên
tử, khi thời gian lớn thì thành phần này không thay đổi và sự phụ thuộc vào
thời gian của
thẳng này ta tính được hệ số khuếch tán (D).
Đối với chất lỏng đơn giản, hệ số khuếch tán vào cỡ 10-5 cm2s-1, chỉ cần
vài trăm bước là đủ để nhận được kết quả tin cậy. Với mẫu chất lỏng nhớt
17
(SiO2, Al2O3…), hệ số tự khuếch tán gần điểm nóng chảy (Tm) khoảng 10-6
cm2s-1, số bước phải tăng lên đến cỡ hàng trăm ngàn mới đạt được kết quả tốt.
Do đó thời gian
vào cỡ 10-7 cm2s-1 hoặc nhỏ hơn thì về mặt thực hành không thể xác định hệ số
này bằng phương pháp động lực học phân tử trong giới hạn tốc độ máy tính
hiện nay.
Tuy nhiên, phương pháp động lực học phân tử lại cho khả năng nghiên
cứu ảnh hưởng của áp suất, nhiệt độ và khối lượng nguyên tử lên hệ số khuếch
tán. Trong trường hợp đơn giản nhất, mô hình chất lỏng có dạng quả cầu cứng,
độ linh động của nguyên tử được xác định bởi hệ số kích cỡ D(m/kT)1/2/
(trong đó m là khối lượng hạt, là đường kính hạt, T là nhiệt độ tuyệt đối, k là
hằng số Boltzmann). Ở nhiệt độ không đổi, hệ số này được tính toán bởi
phương pháp động lực học phân tử, nó giảm dần với thể tích. Vì thế hệ số tự
khuếch tán gần như không thể xác định được cho trường hợp V/Vo < 1.45, với
Vo là thể tích tinh thể xếp chặt tương ứng ở nhiệt độ 0K. Nếu để mô tả sự thay
đổi của hệ số khuếch tán theo nhiệt độ theo phương trình Arrhenius thì kết quả
thu được bị lệch so với thực tế. Lí do là trên thực tế cơ chế khuếch tán trong
chất lỏng không phụ thuộc kích hoạt tự nhiên, mà phụ thuộc vào nhiệt độ theo
công thức .
Khi mô phỏng hệ phi tinh thể với thế tương tác gần, sự gia tăng áp suất
làm giảm hệ số tự khuếch tán. Trong trường hợp sử dụng thế cặp Gauss:
(1.12)
Stilling và Weber đã phát hiện sự gia tăng bất bình thường của hệ số tự
khuếch tán phụ thuộc áp suất của hệ. Hiệu ứng này xuất hiện vì khi khoảng
cách giảm thì hàm lực –du(r)/dr vượt qua giá trị cực đại và sau đó giảm theo.
Đối với chất lỏng làm lạnh nhanh thì giá trị thu được thường lớn hơn thực tế.
Ví dụ trong mẫu của hệ thống sử dụng thế Lennad - Jones, giá trị của hệ số
18
khuếch tán tại nhiệt độ 0,1Tm là 10-7 cm2s-1, nghĩa là quá lớn so với thực tế. Với
mẫu KCl và SiO2 cũng có kết quả tương tự. Sự sai lệch này có thể do kích
thước quá nhỏ của mẫu hoặc do hiệu ứng lực tương tác xa, thực tế tương tác
này là rất nhỏ rất khó xác định nên thường bỏ qua. Trong tương lai, với việc
ứng dụng những thuật toán mới và máy tính hiện đại chúng ta hy vọng có thể
nghiên cứu hệ với hàng triệu hạt.
Các mẫu vật liệu khác nhau và các phương pháp phân tích khác nhau đã
đạt được thành công về việc mô tả cơ chế khuếch tán tập thể trong chất lỏng
siêu lạnh. Người ta dự đoán rằng các nguyên tử có thể chia làm hai loại: một
loại như chất rắn không thể dịch chuyển, loại còn lại như chất lỏng có thể tham
gia vào quá trình di chuyển khuếch tán. Thể tích hiệu dụng của ô lưới mà trong
đó hạt đang lưu trú và tỷ lệ giữa thể tích này và thể tích tới hạn cho phép xác
định nguyên tử thuộc loại nào. Ví dụ trong hệ với tương tác cặp hoàn toàn là
đẩy thuộc dạng: , trong đó r là khoảng cách nguyên tử, và ro
là các hằng số, vai trò quyết định thuộc về hệ số kích cỡ:
(1.13)
Với N là số hạt trong thể tích V. Khi * của vật liệu phi tinh thể tăng từ
1,1 đến 1,5 thì thành phần hạt như trong chất lỏng giảm từ 1 đến 0, nghĩa là ảnh
hưởng lên khuếch tán. Tiêu chuẩn này không đơn thuần là thể hiện tính cấu trúc
mà còn thể hiện ý tưởng về mối tương quan giữa độ linh động của hạt và thể
tích tự do.
Trong vùng ổn định của trạng thái lỏng, phương pháp động lực học phân
tử khiến ta dễ dàng và chắc chắn trong việc tính hệ số khuếch tán nếu thế tương
tác giữa các hạt được biết tới. Mặc dù vậy, trong nhiều năm gần đây, dường
như không có nhiều sự phát triển thực sự trong việc phân tích cơ chế khuếch
tán tập thể cho chất lỏng ở vùng lạnh sâu.
Quá trình khuếch tán của Si và O thông qua cơ chế chuyển động tập thể
của cả khối SiOx (x=4, 5, 6).
19
Chương 2
PHƯƠNG PHÁP TÍNH TOÁN
Trong chương này, kỹ thuật mô phỏng hệ ôxít SiO2 ở trạng lỏng sẽ được
trình bày chi tiết. Trước tiên, chúng tôi trình bày phương pháp mô phỏng động
lực học phân tử, thế Born-Mayer-Huggins, gần đúng Ewald-Hansen. Tiếp theo
chúng tôi trình bày các kỹ thuật phân tích cấu trúc như : Hàm phân bố xuyên
tâm, phân bố số phối trí, phân bố góc liên kết. Cuối cùng, chúng tôi trình bày
cơ chế khuếch tán thông qua độ dịch chuyển bình phương trung bình và các
loại phản ứng trong đám thông qua phản ứng SiOxSiOx' với x, x'=4, 5 và
OSiy OSiy' với y, y'=2.
2.1. Phương pháp động lực học phân tử
Phương pháp ĐLHPT là một công cụ cho phép chúng ta xây dựng mô
hình vật liệu dựa trên hệ phương trình chuyển động của Newton. Phương trình
chuyển động được khảo sát với vận tốc chuyển động của hạt tính bằng thuật
toán Verlet theo bước thời gian dt.
Xét một hệ gồm N nguyên tử được gieo vào khối hình lập phương cạnh
L. Tọa độ ban đầu của các nguyên tử có thể lấy ngẫu nhiên nhưng phải thoả
mãn điều kiện không có bất kỳ hai nguyên tử nào quá gần nhau. Dưới tác dụng
của lực tương tác, các nguyên tử sẽ dịch chuyển dần đến vị trí cân bằng. Trạng
thái cân bằng của mô hình được xác định bởi nhiệt độ và áp suất. Chuyển động
của các nguyên tử trong mô hình tuân theo định luật cơ học cổ điển Newton.
Đối với hệ gồm N hạt, phương trình chuyển động của định luật hai Newton có
thể viết như sau:
Phương pháp ĐLHPT dựa trên phương trình chuyển động Newton:
(2.1)
(2.2)
20
trong đó, Fi là lực tổng hợp tác dụng lên nguyên tử thứ i từ các nguyên tử
còn lại; mi và ai lần lượt là khối lượng và gia tốc của nguyên tử thứ i. Lực Fi
được xác định theo công thức:
(2.3)
trong đó, là thế tương tác giữa nguyên tử thứ i và nguyên tử thứ j và
rij là khoảng cách giữa chúng.
Trong mô phỏng ĐLHPT, ta sử dụng thuật toán Verlet để giải hệ phương
trình chuyển động của các nguyên tử theo định luật hai Newton. Trong thuật
toán này, toạ độ của nguyên tử i ở thời điểm (t + dt) được xác định thông qua
tọa độ của nó ở hai thời điểm t và (t - dt) bằng biểu thức:
(2.4)
Vận tốc ở thời điểm t được xác định thông qua tọa độ ở thời điểm (t - dt)
và (t + dt) theo biểu thức:
(2.5)
Lực Fi(t) được phân tích theo ba thành phần tương ứng với các phương
Ox, Oy và Oz của hệ tọa độ Đề các:
(2.6)
trong đó được xác định như sau:
(2.7)
với là véctơ đơn vị của trục Ox. Các thành phần được xác
định tương tự như phương trình (2.7).
Khi nghiên cứu các mô hình vật liệu bằng phương pháp ĐLHPT, tuỳ
theo mục đích cần nghiên cứu mà người ta thường chọn một trong các mô hình
sau đây: mô hình NVE, NVT, NPH, NTP, TV và TP. Trong đó: N, E, V, T,
21
P, H và lần lượt là số nguyên tử, năng lượng toàn phần, thể tích, nhiệt độ, áp
suất, entanpy và thế hoá học. Đối với mô hình NVE thì các đại lượng N, V và E
không đổi trong suốt thời gian mô phỏng. Còn đối với các mô hình khác sẽ có
các đại lượng tương ứng không thay đổi.
Trong quá trình mô phỏng ĐLHPT, U và K lần lượt là thế năng và động
năng của hệ và được tính theo biểu thức sau:
(2.8)
(2.9)
Năng lượng E của hệ có thể tính theo công thức:
(2.10)
Nhiệt độ của mô hình ĐLHPT có thể được xác định thông qua động năng
của hệ theo công thức:
(2.11)
trong đó kB là hằng số Boltzman.
Động năng của hệ được xác định thông qua vận tốc của các nguyên tử
theo công thức:
(2.12)
Trong mô hình NVT, để giữ nhiệt độ có giá trị không đổi người ta
thường sử dụng kỹ thuật điều chỉnh nhiệt độ (Temperature Scaling). Ý tưởng
của thuật toán này là điều chỉnh vận tốc của tất cả các hạt bởi một thừa số được
xác định bởi tỷ số giữa nhiệt độ mong muốn và nhiệt độ hiện tại được xác định
từ phương trình (2.11). Giả sử nhiệt độ được tính từ phương trình (2.10) là T,
của tất cả các nhiệt độ mong muốn của hệ đạt được là T0, điều chỉnh vận tốc
nguyên tử theo phương trình sau:
22
(2.13)
chúng ta sẽ thu được:
(2.14)
Chọn áp suất của mô hình ĐLHPT có thể được điều chỉnh thông qua
kích thước của mô hình. Mô hình NPT sẽ điều chỉnh áp suất P thông qua
việc nhân tọa độ của tất cả các nguyên tử với thừa số điều chỉnh λ. Khi áp
suất của hệ nhỏ hơn giá trị cho phép, ta sẽ λ > 1, và ngược lại nếu áp suất lớn
hơn giá trị cho trước ta chọn λ < 1. Trong chương trình, áp suất được điều
chỉnh như sau: Nhập giá trị áp suất Pmới, nếu Pmới >Phệ thì λ = 1- dP, ngược
lại λ = 1+ dP với giá trị dP được chọn là 10-4. Do vậy, toạ độ mới của các
nguyên tử được xác định:
(2.15)
Khi đó, kích thước mô hình sẽ có giá trị L’ = Lλ
Khi xây dựng mô hình ĐLHPT, các thông số nhiệt độ và áp suất ở thời
điểm t được xác định như sau:
(2.16)
(2.17)
(2.18)
Các mô hình mô phỏng NVE, NPT, NVT được sử dụng trong mô phỏng .
Mô hình NVE cô lập với môi trường bên ngoài do vậy hầu như không chịu tác
động của ngoại lực. Đây là mô hình có thể sử dụng để khảo sát sự dịch chuyển
của các nguyên tử mô hình và từ đó có thể tính được hệ số tự khuếch tán của
các nguyên tử. Nhược điểm của mô hình NVE là để khảo sát ở nhiệt độ T và áp
23
suất P cho trước ta phải thực hiện một số rất lớn các bước lặp ĐLHPT, do đó
thời gian mô phỏng sẽ kéo dài. Để khắc phục nhược điểm trên, ban đầu chúng
ta mô phỏng theo mô hình NPT hoặc NVT để đạt được các thông số T và P đã
cho. Bước tiếp đến, thực hiện mô phỏng theo mô hình NVE, do đó thời gian
mô phỏng sẽ được giảm đi rất nhiều.
2.2. Thế tương tác
Cho đến nay, đã có một một số hàm thế tương tác được phát triển để mô
phỏng ô xít SiO2 ở cả trạng thái lỏng, tinh thể và VĐH [2, tr.5]. Tuy nhiên, yêu
cầu cần thiết đối với các hàm thế đó là phải bảo đảm sự ổn định của các tứ diện
SiO4 và phải thể hiện được bản chất của liên kết Si-O. Một phần của liên kết
này có bản chất cộng hoá trị và một phần có bản chất ion. Thêm nữa, hàm thế
cũng phải có giá trị thích hợp khi mô phỏng các dạng thù hình khác nhau của ô
xít SiO2. Một số hàm thế thoả mãn điều kiện trên là thế tương tác cặp, một số
hàm thế khác kết hợp cả thế tương tác cặp và thế ba thành phần. Nhiều hàm thế
được rút ra từ dạng thế Born-Mayer-Huggins như sau:
(2.19)
Trong phương trình (2.19), là thế tương tác giữa nguyên tử thứ i và
nguyên tử thứ j; là khoảng cách giữa chúng; qi, qj là điện tích của các
nguyên tử i, j; Aij, Bij, Cij và Dij là các thông số tương ứng của hàm thế. Hai
thành phần sau cùng ở vế phải của phương trình (2.19) liên quan đến tính đa
cực của các nguyên tử và có thể bỏ qua nếu như các nguyên tử không có sự
phân cực.
Để nghiên cứu cấu trúc và cơ chế khuếch tán của ô xít SiO2, chúng tôi đã
sử dụng thế tương tác cặp BKS. Tuy loại thế này khá đơn giản nhưng thực tế
cho thấy việc sử dụng nó đã cho phép mô phỏng được nhiều tính chất của SiO2
ở trạng lỏng (như cấu trúc, mật độ và sự giãn nở nhiệt, v.v... ) Các kết quả mô
phỏng phù hợp tốt với thực nghiệm. Thế tương tác này có dạng:
24
(2.20)
trong đó, rij là khoảng cách giữa hai nguyên tử thứ i và j. Các hệ số của
hàm thế được cho trong bảng sau:
Bảng 2.1. Các hệ số của thế BKS đối với hệ SiO2
Cặp Aij Bij Cij Điện tích (e) nguyên tử (eV) (Å-1 ) (eV Å6)
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
2.3. Gần đúng Ewald-Hansen
Các hàm thế tương tác sử dụng trong mô phỏng thường có hai phần:
tương tác gần và tương tác xa. Thế tương tác U(r) ~ được xem là tương tác
gần khi s >3. Khi s 3 thế V(r) được xem như là thế tương tác xa. Xét một hệ
gồm các nguyên tử được chứa trong một hình cầu bán kính R, năng lượng
tương tác giữa các nguyên tử tỷ lệ với ( r là khoảng cách giữa hai nguyên
tử). Thế năng tác dụng nên nguyên tử thứ i là:
(2.21)
trong đó N là số nguyên tử trong hình cầu. Giả sử hệ là đẳng hướng với
mật độ , phương trình (2.21) có thể viết lại như sau:
(2.22)
Trong luận văn này, chúng tôi đã sử dụng kỹ thuật gần đúng Ewald đối
với tương tác xa (tương tác coulomb) và kỹ thuật ngắt tương tác đối với tương
tác gần.
25
Xét một hệ gồm N nguyên tử đặt trong không gian tính toán là một hình
lập phương có kích thước L. Thế tương tác của các nguyên tử trong hệ và các
ảnh của nó với điều kiện biên tuần hoàn được xác định như sau:
(2.23)
Trong đó qi , qj lần lượt là điện tích của nguyên tử thứ i và j; n là vectơ
toạ độ tâm của các hình hộp ảnh và n = (n1, n2, n3) = n1.Lx + n2.Ly +n3.Lz với x,
y, z là vectơ đơn vị trong hệ toạ độ Đề các. Tâm hình hộp chứa hệ N nguyên tử
mô phỏng có tọa độ n = (0, 0, 0) còn các hình hộp ảnh có tâm ở toạ độ L.n theo
ba chiều với n tiến đến vô cùng. Hình 2.1 mô tả mô hình tính toán gần đúng
Ewald –Hansen trong không gian 2 chiều, mạng tuần hoàn 3x3 được dựng lên
từ ô cơ sở có tâm n(0,0)
n(-1,1) n(0,1) n(1,1)
n(-1,0) n(0,0) n(1,0)
n(-1,-1) n(0,-1) n(1,-1)
Hình 2.1. Mô hình tính toán gần đúng Ewald –Hansen trong không gian 2 chiều, mạng tuần hoàn 3x3 được dựng lên từ ô cơ sở có tâm n(0,0).
Trong biểu thức, tổng đầu tiên chỉ xét điều kiện i j với n = 0; rij,n là
khoảng cách giữa một nguyên tử với các nguyên tử khác trong hệ hoặc trong
các hình hộp ảnh.
(2.24)
Trong mô phỏng ĐLHPT, việc xác định thế tương tác Coulomb chiếm
phần lớn thời gian mô phỏng. Tuy nhiên, sử dụng gần đúng Ewald-Hansen cho
phép giảm đáng kể thời gian tính toán.
2.4. Xác định các đặc trưng vi cấu trúc và tính chất của mô hình
Nội dung của phần này tập trung vào việc trình bày lần lượt cách xác
định các đại lượng đặc trưng cho vi cấu trúc của vật liệu như hàm phân bố
xuyên tâm, phân bố số phối trí, phân bố góc. Ngoài ra, cách xác định hệ số
khuếch tán, một đại lượng quan trọng liên quan mật thiết đến vi cấu trúc SiO2
lỏng cũng sẽ được trình bày chi tiết.
26
2.4.1. Hàm phân bố xuyên tâm
Trong mô phỏng vật liệu ở trạng thái lỏng một đại lượng tuân theo quy
tắc thống kê được sử dụng để xác định vi cấu trúc của vật liệu ở mức nguyên tử,
đó là hàm phân bố xuyên tâm (HPBXT). HPBXT cũng có thể được xác định từ
thực nghiệm thông qua thừa số cấu trúc. Thông qua HPBXT các phân bố về số
phối trí trung bình, khoảng cách liên kết trung bình và góc liên kết trung bình
cũng sẽ được xác định.
Theo cơ học thống kê, hàm tương quan cặp hay HPBXT g(r) được xác
định như sau:
(2.25)
trong đó V là thể tích của mẫu vật liệu và N chính là số nguyên tử chứa
trong thể tích V đó. Phương trình (2.25) có thể viết lại một cách tường minh
hơn như sau:
(2.26)
ở đây rij = ri - rj và ri, rj là véc tơ toạ độ của các hạt thứ i và thứ j. Véc tơ r
là một thông số xuất hiện như một biến thực ở vế trái của phương trình (giá trị
của r do chúng ta chọn). Hàm g(r) có thể hiểu là tỷ lệ thuận với xác suất tìm
thấy nguyên tử cách nguyên tử trung tâm một véc tơ r. Đối với hệ đẳng hướng,
g(r) chỉ phụ thuộc vào độ dài của véc tơ r. Lấy tích phân qua thể tích V(r, r)
giữa r và r+dr và giả sử rằng lớp vỏ hình cầu là đủ mỏng chúng ta sẽ thu được:
(2.27)
Thay phương trình (2.26) vào (2.27) chúng ta có:
(2.28)
27
Tích phân hàm delta, chúng ta sẽ tính được số hạt trong lớp hình cầu là
ni(r,r) và
(2.29)
Thay (2.29) vào (2.28) chúng ta tìm được:
(2.30)
Phương trình (2.30) có thể viết lại một cách đơn giản như sau:
(2.31)
với 0 chính là mật độ nguyên tử trung bình trong thể tích V của mẫu vật
liệu và (r) là mật độ nguyên tử ở khoảng cách r tính từ nguyên tử trung tâm.
(2.32)
HPBXT cũng có thể được xác định từ thực nghiệm. Đại lượng có thể
đo được trực tiếp từ thực nghiệm nhiễu xạ là cường độ nhiễu xạ I(). Trong
đó, là góc giữa tia tới và tia tán xạ. Gọi kin và kout tương ứng là véc tơ sóng
tới và véc tơ sóng tán xạ. Bởi vì tán xạ là đàn hồi, |kin|=|kout|, với k=kin- kout
chúng ta có:
(2.33)
Cường độ tán xạ có thể được tách thành hai phần: thừa số dạng nguyên
tử f(K) và thừa số cấu trúc S(K) như sau:
(2.34)
28
Thừa số hình dạng đặc trưng cho loại nguyên tử và phụ thuộc vào việc
hiệu chỉnh thiết bị đo. Thừa số cấu trúc được xác định bởi phương trình (2.35)
và chứa tất cả các thông tin về vị trí của các nguyên tử.
(2.35)
Liên hệ giữa thừa số cấu trúc với hàm HPBXT, chúng ta dùng định nghĩa
chuẩn về hàm phân bố xuyên tâm (2.25) và biểu diễn chuyển đổi Fourier của
hàm Dirac delta như trong phương trình sau:
(2.36)
Tách tổng trong phương trình (2.35) thành 2 phần ứng với l=m và lm
(2.37)
Áp dụng chuyển đổi fourier trong không gian ba chiều đối với (2.35)
chúng ta được:
(2.38)
Thay phương trình (2.25) vào (2.38) chúng ta được:
(2.39)
Giản ước hàm delta trong phương trình (2.39) chúng ta được:
(2.40)
Từ phương trình (2.40) chúng ta thấy HPBXT, g(r) có thể được xác định
từ thực nghiệm thông qua thừa số cấu trúc. Trong luận văn này các HPBXT
thành phần của hệ SiO2 lỏng được tính trong chương trình mô phỏng như sau:
29
(2.41)
Ở đây, N là tổng số nguyên tử trong mô hình, N và N lần lượt là số
nguyên tử loại và loại ; 0 là mật độ nguyên tử trung bình trong thể tích V.
Hàm phân bố xuyên tâm tổng cộng được xác định như sau:
(2.42)
Trong đó ci, cj là nồng độ của nguyên tử loại α và β, bi và bj là hệ số tán
xạ cho nguyên tử loại α và β.
2.4.2. Xác định số phối trí và độ dài liên kết
Số phối trí trung bình Z được xác định bằng biểu thức tích phân píc thứ
nhất HPBXT tương ứng:
(2.43)
trong đó, rc là bán kính ngắt, thường được chọn là vị trí cực tiểu ngay sau
píc thứ nhất của HPBXT g(r). Giá trị của Z cho ta biết trong hình cầu có
tâm ở vị trí của một nguyên tử loại và bán kính là rc, có bao nhiêu nguyên tử
loại . Từ vị trí của đỉnh thứ nhất trong các hàm phân bố xuyên tâm thành phần
cho phép ta xác định được độ dài liên kết giữa các cặp nguyên tử. Cụ thể, từ vị
trí đỉnh thứ nhất của HPBXT thành phần gSi-Si(r) ta suy ra được khoảng cách lân
cận gần nhất giữa hai nguyên tử Si-Si. Tương tự, ta có thể tính được độ dài liên
kết giữa các cặp nguyên tử O-Si và O-O từ vị trí đỉnh thứ nhất của HPBXT
thành phần gO-Si(r) và gO-O(r).
30
2.4.3. Xác định phân bố góc
Như chúng ta đã biết, cấu trúc của ôxít SiO2 ở trạng thái lỏng được tạo
thành từ các đơn vị cấu trúc cơ bản SiO4, SiO5, SiO6. Phân bố góc Si-O-Si sẽ
cho biết liên kết của các nguyên tử trong một đơn vị cấu trúc cơ bản SiOx
(thông tin về trật tự gần). Còn phân bố góc O-Si-O cho biết liên kết giữa các
đơn vị cấu trúc cơ bản. Như vậy, từ kết quả tính phân bố góc có thể xác định
được sự thay đổi trật tự gần trong các đơn vị cấu trúc cũng như sự thay đổi trật
tự trong khoảng trung, liên quan đến sự thay đổi liên kết giữa các đơn vị cấu
trúc. Trong các công trình [22,9] các kết quả tính toán phân bố góc chưa cung
cấp được thông tin chi tiết về trật tự gần của các đơn vị cấu trúc cơ bản. Phân
bố góc O-Si-O trong các công trình trên nhận được bằng cách tính trung bình
cho tất cả các đơn vị cấu trúc trong hệ. Do vậy, trong luận văn này, việc tính
toán phân bố góc được thiết lập cho từng đơn vị cấu trúc SiOx một cách riêng
biệt. Vì thế so với kết quả trong các công trình trên, kết quả tính toán trong luận
văn này cung cấp nhiều thông tin chi tiết hơn về vi cấu trúc.
Thuật toán sử dụng xác định phân bố góc trong luận văn này được thiết
lập như sau: Trước tiên, chúng tôi tiến hành xác định tập hợp các đơn vị cấu
trúc cơ sở SiOx. Sau đó, ứng với mỗi loại đơn vị cấu trúc, chúng tôi xác định
phân bố góc O-Si-O. Tuy nhiên, đối với phân bố góc Si-O-Si, chúng tôi cũng
tiến hành xác định tất cả các đơn vị cấu trúc OSiy (y=2, 3, 4) trước, sau đó xác
định phân bố góc Si-O-Si.
Để xác định góc O-Si-O hoặc Si-O-Si, khi biết toạ độ của các nguyên tử
tương ứng, chúng tôi làm như sau: giả sử chúng ta xét một tập gồm ba nguyên
tử với toạ độ tương ứng: O1(x1,y1,z1); Si(x2,y2,z2); O2(x3,y3,z3). Góc O-Si-O
được xác định bằng biểu thức:
31
(2.44)
trong đó:
(2.45)
Ở đây, trật tự của các toạ độ x1, x2, x3 (tương tự với toạ độ y và z) đóng
vai trò quan trọng trong việc xác định sự phân bố góc. Góc Si-O-Si được xác
định hoàn toàn tương tự như đối với góc O-Si-O.
2.4.4. Xác định cơ chế khuếch tán
Các mẫu vật liệu khác nhau và các phương pháp phân tích khác nhau
đã đạt được thành công về việc mô tả cơ chế khuếch tán tập thể trong chất
lỏng siêu lạnh. Quá trình khuếch tán trong vật liệu ô xít SiO 2 ở trạng thái
lỏng là chủ đề của những nghiên cứu chuyên sâu trong các thập kỷ qua. Kể
từ năm 1965, Adam và Gibbs [3] đề xuất một lý thuyết để mô tả sự chuyển
động của chất lỏng siêu lạnh. Lý thuyết này cho rằng các chất lỏng thay đổi
cấu trúc của nó theo các vùng độc lập. Các vùng có chứa các nguyên tử lại
có xu hướng hình thành các đám, từ đó phát triển các khái niệm về động
lực học không gian không đồng nhất. Đối với các chất lỏng động không
đồng nhất cũng được xây dựng bởi các thí nghiệm và mô phỏng. Kết quả
thấy rằng hầu hết các nguyên tử linh động và không linh động đã hình
thành các đi chuyển chậm và nhanh trong không gian. Chúng tôi đã nghiên
cứu sự khuếch tán trong SiO2 lỏng sau một cách tiếp cận rằng sự di chuyển
các hạt được đánh giá thông qua phản ứng SiOxSiOx' với x, x'=4, 5 và
OSiy OSiy' với y, y'=2.
32
A)
4 4 4 2 2 3 12 12 12 3 1 3 1 6 1 6
B)
4 4 4 2 2 2
10 10 10 3 3 3 1 5 1 1 5
Hình 2.2. Mô hình hóa các loại phản ứng trong SiO2 lỏng
A) là phản ứng SiO4SiO5, SiO5SiO4 và liên kết Si12- O2 được thay thế bởi liên kết Si12- O6;
B) là phản ứng SiO5SiO4, SiO4SiO5 và các liên kết Si10- O5 bị hỏng, sau đó phục hồi. Ở đây, hình cầu
màu đen và màu xám vẽ nguyên tử O và nguyên tử Si; đường nối các nguyên tử Si-O là liên kết Si-O.
Hình 2.2 minh họa việc sắp xếp các nguyên tử trong các ô phối trí khi hai
phản ứng xảy ra: ở hình 2.2(a) là phản ứng SiO4SiO5, SiO5SiO4 xảy ra và
liên kết Si12- O2 được thay thế bởi Si12- O6; hình 2.2(b) là phản ứng
SiO5SiO4, SiO4SiO5 xảy ra khi các liên kết Si10- O5 bị hỏng, sau đó phục
hồi. Khi một phản ứng xảy ra, liên kết hiện tại bị phá vỡ hoặc liên kết mới được
tạo ra: trong trường hợp thứ nhất một liên kết được thay thế bằng một liên kết
mới, trong trường hợp thứ hai một liên kết bị phá vỡ sau đó hồi phục, các
nguyên tử trong trường hợp thứ hai gần như chỉ dao động xung quanh vị trí cân
bằng cố định. Phản ứng trong trường hợp đầu tiên gọi là phản ứng có hiệu quả
và nó chiếm ưu thế. Chúng tôi nghiên cứu tỉ lệ phản ứng có hiệu quả để chỉ ra
sự khuếch tán trong chất lỏng SiO2. Hằng số khuếch tán có thể được mô tả bởi:
33
(2.46)
Trong đó
tử cho mỗi phản ứng có hiệu quả; D là tốc độ phản ứng; là tỉ lệ phản ứng có
hiệu quả.
Một đám được định nghĩa là một tập các nguyên tử mà mỗi nguyên tử có
thể kết nối với nhau thông qua một đường dẫn bao gồm các liên kết. Việc xác
định các đám từ các nguyên tử N được thực hiện như sau. Thứ nhất, tất cả các
liên kết được xác định hình thành từ N các nguyên tử. Sau đó, mỗi nguyên tử
được gán cho một nhãn đám i; i = 1, 2, .. N. Tiếp theo, nguyên tử k1 sẽ được
gán lại cho nhãn mới, nếu nó tạo thành một liên kết với một nguyên tử k2 và
k1≠ k2. Cả hai nguyên tử được gán lại cho nhãn k1 nếu k1 Thủ tục gán nhãn đám đã được thực hiện cho đến khi tất cả các cặp nguyên tử tạo thành một liên kết có cùng một nhãn. Mô phỏng cho thấy trong SiO2 chứa một lượng nhỏ các khuyết tật (sai hỏng mạng). Khuyết tật phân bố không đồng nhất trong chất lỏng mà chúng có xu hướng di chuyển lại gần nhau và tạo thành đám. Quá trình khuếch tán thực hiện chủ yếu thông qua phản ứng SiOxSiOx' với x, x'=4, 5 và OSiy OSiy' với y, y'=2, 3. Phản ứng SiOxSiOx' xảy ra không phải ngẫu nhiên mà có xu hướng tạo thành đám. Chúng tôi phát hiện trong SiO2 có hai loại đám. Đám thứ nhất gồm các nguyên tử mà không có phản ứng SiOxSiOx' nào xảy ra. Đám thứ hai có nguyên tử ở đó phản ứng SiOxSiOx' xảy ra tuần hoàn. Chúng tôi xác định độ linh động của nguyên tử trong đám không phản ứng so với đám phản ứng và xác định kích thước của hai đám mô tả phụ thuộc vàhời gian mô phỏng và nhiệt độ. Từ đó xác định được cơ chế khuếch tán trong SiO2 lỏng. 34 Chương 3 KẾT QUẢ VÀ THẢO LUẬN trong khoảng áp suất từ 0 đến 25 GPa ở nhiệt độ 3200 K. Tiếp theo chúng tôi khảo sát tính đa thù hình của SiO2 lỏng. Cuối cùng, chúng tôi khảo sát cơ chế khuếch tán của nguyên tử Si và nguyên tử O theo quan điểm khuếch tán theo đám địa phương. 3.1. Khảo sát cấu trúc của SiO2 lỏng theo áp suất Để nghiên cứu cấu trúc của SiO2 lỏng theo áp suất ở nhiệt độ 3200 K, chúng tôi tạo ra sáu mẫu SiO2 chứa 1998 nguyên tử, trong đó có 666 nguyên tử
Si và 1332 nguyên tử O, mật độ nằm trong khoảng từ 2,634 đến 3,991 g/cm3, các mẫu này được tạo ra bằng cách nén từ mẫu có mật độ thấp. Ở đây, chúng tôi sử dụng thế tương tác BKS và thuật toán Verlet để lấy tính phân phương trình chuyển động. Chi tiết về thế BKS và thuật toán Verlet được trình bày chi tiết trong chương 2. Bước thời gian trong mô phỏng là 0,47 fs. Cấu hình đầu tiên của mẫu được tạo ra bằng cách gieo ngẫu nhiên các nguyên tử trong không gian mô phỏng. Cấu hình này được nung nóng tới nhiệt độ 5000 K để phá vỡ trạng thái nhớ ban đầu và được giữ ở nhiệt độ này trong
hơn 5×104 bước thời gian mô phỏng. Sau đó, mẫu vật liệu được làm nguội xuống các nhiệt độ 4500 K, 4000 K và cuối cùng là 3200 K. Ở các nhiệt độ
4500K và 4000 K, mẫu được hồi phục trong khoảng 2×105 bước thời gian. Ở
nhiệt độ 3200 K, mẫu được hồi phục trong 106 bước thời gian. Sau khi nhận được mẫu SiO2 lỏng ở nhiệt độ 3200 K và ở áp 0 GPa, chúng tôi đã tạo ra sáu mẫu SiO2 ở các áp suất nén khác nhau (mật độ khác nhau). Chúng tôi sử dụng mô hình NPT (số hạt, áp suất, nhiệt độ không đổi) để tạo ra các mẫu có áp suất khác nhau. Sau đó, các mẫu có áp suất khác nhau được hồi phục ở thể tích không đổi bằng cách sử dụng mô hình NVE (số hạt, thể tích, năng lượng không
đổi). Mẫu đạt được trạng thái cân bằng sau khoảng 106 bước thời gian mô phỏng. Khi mẫu đạt được trạng thái ổn định, chúng tôi đã tiến hành xác định các đặc trưng cấu trúc và tính đa thù hình của mẫu. 35 )
r
(
i
S
-
i
S
g 36 Để đảm bảo độ tin cậy, chúng tôi đã tính toán các đặc trưng cấu trúc cơ bản như HPBXT, SPT cặp trung bình và so sánh kết quả mô phỏng với kết quả thực nghiệm cũng như kết quả mô phỏng của các tác giả khác. Trên hình 3.1 vẽ HPBXT thành phần Si-Si, Si-O, O-O của SiO2 lỏng ở áp suất khác nhau và ở nhiệt độ 3200 K. Như có thể thấy, HPBXT trên 3.1 có hình dạng, vị trí và độ cao của các cực đại tương tự với kết quả tính toán của các tác giả trong các công trình [19,29]. Chi tiết về các đặc trưng của HPBXT và SPT trung bình được liệt kê trong bảng 3.1. Như có thể thấy trong bảng, độ dài liên kết giữa các cặp nguyên tử Si-Si, Si-O, O-O cũng phù hợp tốt với kết quả thực nghiệm (ở áp suất 0 GPa) đã được công bố trong các trong công trình [24]. Kết quả đo thực nghiệm HPBXT của SiO2 lỏng của một số công trình được báo cáo trong [11] đều cho thấy: khoảng cách liên kết Si-O nằm trong khoảng từ 1,60 Å đến 1,62 Å; khoảng cách liên kết O-O nằm trong khoảng từ 2,62 Å đến 2,65 Å; khoảng cách liên kết Si-Si là từ 3,10 Å đến 3,13 Å. gij rij(Å) Zij Áp suất, rij GPa 1-2 1-1 1-2 2-2 1-1 2-2 1-1 1-2 2-1 2-2 (Å) 9,1 0,10 3,10 1,60 2,60 0,01 2,89 2,75 4,49 4,07 2,03 8,17 4,87 3,08 1,60 2,56 0,01 2,57 7,22 2,48 5,71 4,40 2,2 11,13 9,83 3,08 1,60 2,50 0,01 2,42 6,14 2,40 6,90 4,78 2,39 12,95 15,73 3,08 1,62 2,50 0,01 2,38 5,67 2,41 7,96 5,08 2,54 13,91 20,15 3,08 1,62 2,46 0,01 2,36 5,42 2,43 8,42 5,31 2,65 14,61 25,20 3,08 1,64 2,44 0,01 2,35 5,29 2,46 8,98 5,50 2,75 15,20 [13] 3,12 1,62 2,65 rij, gij là vị trí, độ cao của đỉnh thứ nhất của các hàm phân bố xuyên tâm thành phần; rij là sai số của rij; Zij là số phối trí cặp trung bình. Trong đó, 1-1 là cặp Si-Si; 1-2 là cặp Si-O; 2-1 là cặp O-Si và 2-2 là cặp O-O Kết quả tính toán SPT cặp trung bình trên bảng 3.1 chỉ ra rằng ở áp suất 0,10 GPa phần lớn các nguyên tử Si được bao quanh bởi 4 nguyên tử O, và 37 phần lớn các nguyên tử O liên kết với hai nguyên tử Si. Khi áp suất tăng, SPT của tất cả các cặp Si-O, O-O và Si-Si đều tăng. Kết quả tính toán SPT trung bình một lần nữa cho thấy sự phù hợp khá tốt với số liệu tính toán đã được công bố trong các các công trình [27,29]. (a) áp suất 0,10 GPa; (b) 25,20 GPa; O là quả cầu mầu đen, Si là trắng. 38 Hình 3.2 cho thấy ảnh chụp cấu trúc mạng của SiO2 lỏng ở áp suất 0 GPa (Hình 3.2 a) và ở áp suất 25 GPa (Hình 3.2 b). Cấu trúc được vẽ trong hình lập phương có các cạnh là 20 Å. Như có thể thấy, sự thay đổi cấu trúc mạng theo áp suất của mẫu SiO2 lỏng có thể quan sát một cách trực quan trên cấu trúc mạng. Kết qủa cho thấy, hầu như không quan sát thấy SPT của cặp Si-O có giá trị 5, 6 hoặc SPT của cặp O-Si có giá trị 3, 4. Tiếp theo chúng tôi phân tích đặc tính của các đơn vị cấu trúc. Kết quả tính toán SPT cặp trung bình cho thấy cấu trúc của SiO2 được tạo thành từ các đơn vị cấu trúc cơ bản SiO4, SiO5, SiO6 (Hình 3.3). (a) là cấu trúc SiO4; (b) là cấu trúc SiO5; (c) là cấu trúc SiO6 và (d) là liên kết giữa hai đơn vị cấu trúc. Quả cầu mầu đen là nguyên tử O, mầu xám là Si Hình 3.4 vẽ sự phụ thuộc của tỷ lệ các đơn vị cấu trúc SiOx (x = 4, 5, 6) theo áp suất. Có thể thấy, khi áp suất tăng từ 0-25 GPa thì tỷ lệ các đơn vị cấu trúc SiO4 giảm nhanh ngược lại tỷ lệ các đơn vị cấu trúc SiO5 và SiO6 tăng. 39 Đặc biệt, trong khoảng áp suất từ 12 GPa đến 15 GPa thì tỷ lệ đơn vị cấu trúc SiO5 đạt cực đại. Khi áp suất trên 15 GPa thì tỷ lệ đơn vị cấu trúc SiO5 giảm trong khi đó tỷ lệ các đơn vị cấu trúc SiO6 tăng. Ở áp suất 0,10 GPa, tỷ lệ các đơn vị cấu trúc SiO4 chiếm 92%, tỷ lệ các đơn vị cấu trúc SiO5 khoảng 8% còn tỷ lệ các đơn vị cấu trúc SiO6 hầu như không đáng kể. Ở áp suất 25,20 GPa tỷ lệ SiO4 khoảng 5%, SiO5 khoảng 36% và SiO6 khoảng 58%. Ở áp suất cao, các đơn vị cấu trúc SiO3 và SiO7 chiếm tỷ lệ nhỏ hơn 1%. Kết quả phân tích trên chứng tỏ rằng, khi nén áp suất từ 0 lên 25 GPa, trong ôxít SiO2 lỏng có sự chuyển pha cấu trúc từ cấu trúc tứ diện SiO4 sang cấu trúc bát diện SiO6 thông qua các đơn vị cấu trúc SiO5. Như chúng ta đã biết, các đơn vị cấu trúc SiO4, SiO5, SiO6 liên kết với nhau thông qua các nguyên tử cầu O ở đỉnh của các đa diện SiOx tạo thành mạng ngẫu nhiên trong không gian ba chiều. Hai đơn vị cấu trúc SiOx lân cận có thể kết nối với nhau thông qua một, hai, hoặc ba nguyên tử cầu O chung. 40 Các nguyên tử ôxy liên kết giữa các đơn vị cấu trúc gọi là liên kết cầu ôxy. Bảng 3.2 liệt kê sự phụ thuộc của phân bố liên kết cầu ôxy vào áp suất của SiO2. Kết quả cho thấy, ở áp suất thấp hầu hết các liên kết giữa hai đơn vị cấu trúc SiOx thông qua một nguyên tử O dùng chung. Ở áp suất cao, số liên kết thông qua một nguyên tử O dùng chung này giảm, đồng thời trong cấu trúc xuất hiện các liên kết thông qua hai và ba nguyên tử O dùng chung. Khi áp suất tăng, số liên kết thông qua hai nguyên tử O dùng chung tăng đáng kể (khoảng 20,5 % ở áp suất 25,20 GPa). Tuy nhiên, các liên kết thông qua ba nguyên tử O dùng chung chiếm tỷ lệ nhỏ không đáng kể (nhỏ hơn 2% ở tất cả các áp suất nghiên cứu). 1 97,88 91,17 82,52 80,04 77,46 76,78 2 2,05 8,43 16,20 18,01 20,57 21,51 3 0,07 0,40 1,27 1,95 1,96 1,71 3.2. Khảo sát tính đa thù hình của SiO2 lỏng Như đã trình bày trong mục 3.1, chúng ta đã biết cấu trúc của mẫu SiO2 lỏng được tạo thành từ các đơn vị cấu trúc cơ bản SiO4, SiO5 và SiO6. Tỷ lệ của các đơn vị cấu trúc này phụ thuộc mạnh vào mật độ khối lượng hay nói cách khác là phụ thuộc mạnh vào áp suất nén. Nhiều tính chất vật lý liên quan đến tỷ lệ của các đơn vị cấu trúc SiOx (x = 4, 5, 6) trong mẫu. Do đó, trong mục này chúng tôi khảo sát tính đa thù hình của SiO2 lỏng thông qua các đặc tính của các đơn vị cấu trúc SiOx. 41 Thông tin khá thú vị liên quan đến tính đa thù hình của SiO2 lỏng được đặc trưng bởi tính đa thù hình của SiO2 lỏng trong các đơn vị cấu trúc. Hình 3.5 vẽ phân bố khoảng cách liên kết của cặp Si-O trong các đơn vị cấu trúc SiO4, SiO5 và SiO6. Có thể thấy, phân bố khoảng cách liên kết Si-O trong các đơn vị cấu trúc SiOx không phụ thuộc vào áp suất (hay không phụ thuộc vào mật độ). Phân bố khoảng cách liên kết Si-O trong các đơn vị cấu trúc SiO4, SiO5 và SiO6 có các đỉnh cực đại %
m
ă
r
t
n
ầ
h
p
ệ
l
ỉ T Khoảng cách, Å ở khoảng cách liên kết lần lượt khoảng 1,60 Å, 1,64 Å và 1,68 Å. a) SiO4, b) SiO5 và c) SiO6. Thông tin nữa cũng liên quan đến tính đa thù hình của SiO2 lỏng có thể nhận được từ phân bố góc. Hình 3.6 vẽ phân bố góc liên kết trung bình O-Si-O trong các đơn vị cấu trúc SiO4, SiO5 và SiO6. Hình vẽ cho thấy phân bố góc O-Si-O trong các đơn vị cấu trúc phụ thuộc không đáng kể vào áp suất nén (không phụ thuộc vào mật độ). Phân bố góc liên kết trung bình O-Si-O trong các đơn vị cấu
trúc SiO4 đạt cực đại ở góc khoảng 105o. Mặt khác chúng ta biết góc O-Si-O
trong tứ diện đều SiO4 có giá trị là 1090. So sánh với kết quả phân tích trên về 42 trong SiO2 lỏng phân bố góc, chúng ta có thể kết luận các đơn vị cấu trúc SiO4 ) ( %
m
ă
r
t
n
ầ
h
p
ệ
l
ỷ
T Góc (độ) là các tứ diện bị biến dạng (méo). Trên hình 3.6 còn cho thấy phân bố góc liên
kết O-Si-O trong các đơn vị cấu trúc SiO5 xuất hiện hai cực đại đặt tại khoảng 90o
và 155o. Trong khi đó, trong các đơn vị cấu trúc SiO6 có hai cực đại đặt ở khoảng
85o và 160o. Hình 3.7 vẽ phân bố góc liên kết Si-O-Si giữa hai đơn vị cấu trúc SiOx lân cận. Hình vẽ cho thấy phân bố góc liên kết trung bình Si-O-Si giữa các đơn vị cấu trúc SiOx phụ thuộc mạnh theo áp suất. Cụ thể: áp suất thấp 0,10 GPa phân bố góc liên kết Si-O-Si chỉ có một đỉnh cực đại ở góc khoảng 144o, còn áp suất cao phân bố góc Si-O-Si tách thành hai cực đại đặt tại khoảng 95o và 125o. Hình 3.8 vẽ phân bố góc liên kết Si-O-Si trong các đơn vị cấu trúc OSi2, OSi3 và OSi4. Kết quả cho thấy phân bố góc liên kết Si-O-Si trong các đơn vị cấu trúc OSi2 có một cực đại khoảng 145o. Phân bố góc liên kết Si-O-Si trong 43 các đơn vị cấu trúc OSi3 có hai cực đại đặt tại khoảng 100o và khoảng 120o. Phân bố góc liên kết Si-O-Si trong các đơn vị cấu trúc OSi4 có một cực đại đặt tại ) % ( m
ă
r
t
n
ầ
h
p
ệ
l
ỷ
T Góc(độ) khoảng 90o. ) ( %
m
ă
r
t
n
ầ
h
p
ệ
l
ỷ
T Góc(độ) (a) là cấu trúc OSi2, (b) là cấu trúc OSi3 và (c) là cấu trúc OSi4 44 Hình 3.8 còn cho thấy phân bố góc Si-O-Si không phụ thuộc vào áp suất nén. Như vậy sự thay đổi góc liên kết trung bình Si-O-Si trong là do sự thay đổi tỷ lệ các đơn vị cấu trúc OSiy (y= 2, 3, 4). Ở áp suất thấp các đơn vị cấu trúc OSi2 chiếm đa số, khi áp suất tăng thì tỷ lệ các đơn vị cấu trúc OSi2 giảm, tỷ lệ các đơn vị cấu trúc OSi3 và OSi4 tăng. Tỷ lệ của các đơn vị cấu trúc OSiy thay đổi là do sự thay đổi SPT trung bình của cặp O-Si như đã mô tả trong bảng 3.1. Như vậy, các kết quả trên đã cho thấy cấu trúc của các đa diện SiO4, SiO5, và SiO6 không phụ thuộc vào áp suất. Điều này có nghĩa là cấu trúc của chúng là giống nhau trong các mẫu SiO2 lỏng với mật độ khác nhau, chỉ có tỷ lệ của các đơn vị cấu trúc SiO4, SiO5, và SiO6 là thay đổi. Từ những phân tích trên, có thể khẳng định cấu trúc của mẫu SiO2 lỏng có thể coi được tạo thành từ các đa diện SiO4, SiO5 và SiO6, các đa diện này có cấu trúc không thay đổi. Tuy nhiên, tỷ lệ của các đa diện này phụ thuộc mạnh vào áp suất nén. Do đó một số tính chất của mẫu SiO2 lỏng có thể được biểu diễn như là hàm của tỷ lệ các đơn vị cấu trúc SiO4, SiO5 và SiO6, hàm này được biểu diễn như sau: (3.1) Ở đây CSiO4, CSiO5 và CSiO6 lần lượt là tỷ lệ các đa diện SiO4, SiO5 và SiO6, k, m, l là các hệ số fit, các hệ số được fit từ phương trình (3.3). Từ đó suy ra mật độ của SiO2 lỏng được xác định bởi: (3.2) Ở đây, SiO2 là mật độ của mẫu, đặt SiO4 = 2,564 g/cm3 là mật độ trong vùng mà toàn bộ các đơn vị cấu trúc cơ bản là SiO4, SiO5 = 3,666 g/cm3là mật độ trong vùng mà toàn bộ các đơn vị cấu trúc cơ bản là SiO5 và SiO6 = 4,614 g/cm3 là mật độ trong các vùng mà toàn bộ các đơn vị cấu trúc là SiO6. Như vậy, mật độ của mẫu SiO2 lỏng ở một trạng thái nào đó có thể được xác định thông qua tỷ lệ các đơn vị cấu trúc SiO4, SiO5 và SiO6. 45 Hình 3.9 vẽ sự phụ thuộc vào mật độ của mẫu SiO2 lỏng theo áp suất nén từ 0 đến 25 GPa theo dữ liệu mô phỏng và công thức (3.2). Từ hình vẽ ta thấy, kết quả tính mật độ tính theo công thức (3.2) tương tự với dữ liệu tính được từ mô phỏng. 3.3. Khảo sát cơ chế khuếch tán trong SiO2 lỏng theo nhiệt độ Trong mục này, để nghiên cứu cơ chế khuếch tán theo nhiệt của SiO2 lỏng, chúng tôi đã tạo ra bốn mẫu SiO2 lỏng chứa 1998 nguyên tử, trong đó có 666 nguyên tử Si và 1332 nguyên tử O, các mẫu có nhiệt độ lần lượt là 2600 K, 3000 K, 3200 K và 3500 K. Quá trình chế khuếch tán được chúng tôi khảo sát thông qua hai loại đám đặc biệt. Đám thứ nhất gồm có các nguyên tử mà ở đó không có phản ứng SiOxSiOx' nào xảy ra. Đám thứ hai có các nguyên tử ở đó phản ứng SiOxSiOx' xảy ra liên tục. Hình 3.10 chúng tôi vẽ minh họa các đơn vị cấu trúc SiOx, OSiy, các đám và hai phản ứng SiOx SiOx' và OSiy OSiy'. Như ví dụ trong hình 3.10 d và e phản ứng xảy ra khi cấu trúc SiO5 chuyển thành cấu trúc SiO4, cấu trúc OSi2 chuyển thành cấu trúc OSi3. Sự khuếch tán 46 được thực hiện thông qua hai loại phản ứng: SiOx → SiOx’ và OSiy → OSiy’. Ở đây hầu hết các phản ứng phải thoả mãn điều kiện | x - x '| = 1; | y - y '| = 1. Còn các phản ứng SiO4→SiO6 hoặc SiO6→SiO4 hầu như không xảy ra. Khi phản ứng SiOx → SiOx’ xảy ra, thì một phản ứng OSiy → OSiy’ cũng diễn ra. a) là cấu trúc SiOx; b) là cấu trúc OSiy, c) là đám và d), e) là các phản ứng SiOx SiOx' và OSiy OSiy'. Quả cầu đen và xanh tương ứng với nguyên tử O và Si . Đường mầu đen chỉ liên kết giữa hai nguyên tử. 47 Tiếp theo, chúng tôi khảo sát các đặc tính của các loại phản ứng và các đám trong các mẫu SiO2 lỏng. Hình 3.11 cho thấy số các phản ứng loại khác nhau như là hàm của bước mô phỏng ĐLHPT. Tất cả các điểm nằm phù hợp trên đường thẳng. Độ dốc của các đường này được sử dụng để xác định tốc độ tương tác. Như thấy trong bảng 3.3, tốc độ SiOx SiOx' hầu như bằng tốc độ của OSiy OSiy'. Nó nghĩa là hầu hết tốc độ phản ứng của các loại thỏa mãn điều kiện | x - x' | =1 hoặc | y - y' | =1. Do đó, tốc độ SiOx SiOx' và OSiy OSiy' lớn hơn rất ít tốc độ đó với x, x'=4, 5 và y, y'=2, 3. Điều này dẫn đến hầu hết các ô sai hỏng là SiO5 và OSi3. Tỷ lệ phản ứng 0,2502×10-3 0,2902×10-3 0,3097×10-3 0,5278×10-3 SiOx SiOx' Tỷ lệ phản ứng 0,2502×10-3 0,2906×10-3 0,3099×10-3 0,5276×10-3 OSiy OSiy' Tỷ lệ phản ứng 0,2206×10-3 0,2367×10-3 0,2727×10-3 0,4086×10-3 SiOxSiOx' với x, x'=4, 5 Tỷ lệ phản ứng 0,2225×10-3 0,2375×10-3 0,2641×10-3 0,4100×10-3 OSiy OSiy' với y, y'=2, 3 48 Các phản ứng SiOx và OSiy
Các phản ứng SiOx
Các phản ứng SiO4 và SiO5 g
n
ứ
n
ả
h
p
c
á
c
g
n
ợ
ư
l
ố
S Bước ĐLHPT Hình 3.11. Sự phụ thuộc của các loại phản ứng khác nhau theo bước Các ô sai hỏng
Ô phối trí ngẫu nhiên m
á
đ
ố
S m
á
đ
ố
S t
ấ
h
n
n
ớ
l m
á
đ g
n
o
r
t
ử
t
n
ê
y
u
g
n
ố
S Bước ĐLHPT ĐLHPT. Nhiệt độ là 2600 K. 49 Chất lỏng gồm chủ yếu các ô phối vị SiO4 và OSi2. Một lượng nhỏ các ô sai hỏng được phân bố thông qua mẫu cấu trúc. Đặt Nd(n) là số các ô sai hỏng xác định ở bước n. Nd(n) thay đổi trong khoảng rộng từ 65 đến 110. Các ô sai hỏng tạo thành số đám NC(n). Đại lượng quan tâm ở đây là số nguyên tử bên trong các đám lớn NLC(n). Để làm rõ đặc trưng về phân bố riêng của các ô sai hỏng, chúng tôi lấy ngẫu nhiên Nd(n) nguyên tử trong mẫu cấu trúc. Sau đó chúng tôi tìm các đám tạo bởi các ô phối vị của các nguyên tử này (ô phối trí ngẫu nhiên). Các đại lượng NC(n), NLC(n) được xác định đối với các đám tìm thấy. Kết quả được thấy trên hình 3.12. So với trường hợp các ô phối trí ngẫu nhiên, số đám chứa ô sai hỏng NC(n) giảm đáng kể và số nguyên tử bên trong các đám lớn chứa ô sai hỏng NLC(n) tăng ngược lại. Có nghĩa là các ô sai hỏng phân bố không đồng nhất trong chất lỏng, nhưng chúng có xu hướng gần nhau và tạo thành các đám. Điều này cho thấy, chất lỏng gồm vùng nhiều và vùng ít t
ấ
h
n
n
ớ
l Các đám nguyên tử phản ứng
Các đám nguyên tử không phản ứng m
á
đ
g
n
o
r
t
ử
t
n
ê
y
u
g
n
c
á
c
g
n
ợ
ư
l
ố
S Bước ĐLHPT sai hỏng. 50 Thêm nữa, chúng tôi tìm đám lớn nhất và xác định kích thước của nó. Hình 3.13 cho thấy số nguyên tử trong đám lớn nhất tìm thấy đối với các nguyên tử phản ứng và không phản ứng. Khi thời gian tăng, kích thước của đám không phản ứng giảm, ngược lại kích thước của đám phản ứng tăng. Hơn nữa, kích thước của cả hai đám mô tả trên là bằng nhau tại bước n = 12510. Các kết quả trên cho thấy phản ứng mở rộng trên nhiều nguyên tử hơn trong chất lỏng khi tăng thời gian. Tiếp theo, chúng tôi xác định các phản ứng xảy ra với thời gian định rõ td. Nếu td không đủ lớn, thì các phản ứng trải rộng trên tất cả các nguyên tử và mẫu gồm chỉ một đám tương tác. Trong trường hợp td nhỏ, số các phản ứng và không phản ứng xuất hiện. Đặt Nre là số phản ứng xảy ra với thời gian td. Ở cuối của gian đoạn thời gian td chúng tôi xác định số đám. Sau đó chúng tôi xác định số nguyên tử trong đám lớn nhất trong các đám tìm thấy. Chúng tôi cũng xét trường hợp khi các phản ứng Nre được phân bố ngẫu nhiên thông qua các nguyên tử trong mẫu. Sau đó chúng tôi xác định số đám phản ứng ngẫu nhiên và không ngẫu nhiên. Chúng tôi đã xác định số các phản ứng đối với 500 bước và các khoảng thời gian tiếp theo [0, 500], [500, 1000], ... [1.95×104, 2×104]. Các kết quả được biểu diễn trong hình 3.14 và 3.15 cho thấy so với trường hợp phân bố ngẫu nhiên, số các đám là nhỏ hơn nhưng kích thước đám lớn nhất lại lớn hơn. Điều này có nghĩa các phản ứng phân bố không đồng nhất qua các nguyên tử trong chất lỏng. 51 t
ấ
h
n
n
ớ
l n
ê
y
u
g
n
c
á
c
g
n
ợ
ư
l
ố
S m
á
đ
g
n
o
r
t
ử
t Đám không phản ứng
Đám không phản ứng ngẫu nhiên m
á
đ
g
n
ợ
ư
l
ố
S Khoảng thời gian t
ấ
h
n
n
ớ
l m
á
đ
g
n
o
r
t ử
t
n
ê
y
u
g
n
c
á
c
g
n
ợ
ư
l
ó
S Đám phản ứng
Đám phản ứng ngẫu nhiên m
á
đ
c
á
c
g
n
ợ
ư
l
ố
S Khoảng thời gian 52 . n
ớ
l ử
t
n
ê
y
u
g
n
g
n
ợ
ư
l
ố
S m
á
đ
g
n
o
r
t 2 Å Đám phản ứng
Đám không phản ứng h
n
ì
b
h
c
á
c
g
n
ả
o
h
K ,
h
n
ì
b
g
n
u
r
t
g
n
ơ
ư
h
p Khoảng thời gian n
ớ
l ử
t
n
ê
y
u
g
n
g
n
ợ
ư
l
ố
S m
á
đ
g
n
o
r
t 2 Å Đám phản ứng
Đám không phản ứng h
n
ì
b
h
c
á
c
g
n
ả
o
h
K ,
h
n
ì
b
g
n
u
r
t
g
n
ơ
ư
h
p Khoảng thời gian 53 Trong mục tiếp theo chúng tôi xác định các đám phản ứng và không phản ứng đối với hai trường hợp td bằng 2000 và 5000 bước. Chúng tôi đã tính độ dịch chuyển bình phương trung bình trên nguyên tử trong các trường hợp này. Kết quả tính toán được vẽ trên hình 3.16 và 3.17. Chúng ta có thể thấy rằng kích thước của đám không phản ứng lớn nhất thay đổi trong khoảng 850 đến 1110 đối với trường hợp thứ nhất, và từ 1340 đến 1520 đối với trường hợp hai. Khoảng tương ứng này đối với đám phản ứng là 50-150 và 20-60 đối với trường hợp hai. Số các phản ứng xảy ra ở trung bình trên một nguyên tử trong các đám phản ứng biến đổi từ 46 đến 211 và từ 365 đến 1221 đối với trường hợp td bằng 2000 và 5000 bước. Nó là rất thú vị nhớ rằng độ dịch chuyển bình phương trung bình trên nguyên tử với đám không phản ứng là nhỏ hơn đáng kể so với đám phản ứng. Sự khác nhau về độ dịch chuyển bình phương trung bình trên nguyên tử giữa hai đám là lớn hơn với trường hợp td là 2000 bước. Nghĩa là trong suốt thời gian td chúng tôi tìm thấy hai đám riêng biệt ở đó độ linh động của nguyên tử khác nhau đáng kể với mỗi loại đám. Kích thước của các đám này cũng khác nhau về độ dịch chuyển bình phương trung bình trên nguyên tử giữa các đám phụ thuộc vào thời gian td. Như vậy, khuếch tán diễn ra như sau: Ban đầu, các phản ứng xảy ra trên các nguyên tử tạo thành các đám riêng biệt. Số các đám phản ứng khác nhau rất lớn đối với trường hợp khi các phản ứng được phân bố ngẫu nhiên qua nguyên tử (Hình 3.14 và 3.15). Khi số phản ứng tăng, các đám này phát triển và sau đó nhập thành các đám lớn. Sau đó phản ứng tạo thành đám thấm qua toàn mẫu. Do đó, khuếch tán phụ thuộc không chỉ vào tốc độ phản ứng mà còn phụ thuộc vào các phản ứng được phân bố không gian như thế nào. Do đó, nó là rất thú vị để làm rõ ảnh hưởng của nhiệt độ đến phân bố không gian của các tương tác. 54 Trong hình 3.18 chúng tôi vẽ số nguyên tử (nguyên tử tương tác) ở đó phản ứng xảy ra và số nguyên tử trong đám không phản ứng lớn nhất đối với hai trường hợp, khi nhiệt độ là 2600 và 3000 K. Chúng tôi thấy rằng đối với số các phản ứng giống nhau (số các nguyên tử phản ứng) lớn hơn trong trường hợp 3000 K, cụ thể phản ứng mở rộng trên nhiều nguyên tử hơn đối với trường hợp nhiệt độ cao. Hơn nữa, kích thước của đám không phản ứng lớn hơn đối với trường hợp nhiệt độ thấp lớn hơn (Hình 3.18). Nghĩa là khi nhiệt độ giảm, các phản ứng có xu hướng tập trung vào trong các vùng nhỏ hơn. Xu hướng này được mô tả trên hình 3.19. Mặt khác, các đám không phản ứng trở lên lớn hơn. Như một kết quả, đám không phản ứng ngăn chặn sự khuếch tán thâm nhập vào ử
t
n
ê
y
u
g
n
g
n
ợ
ư
l
ố
S g
n
ứ
n
ả
h
p ử
t
n
ê
y
u
g
n
g
n
ợ
ư
l
ố
S g
n
ô
h
k
m
á
đ
g
n
o
r
t t
ấ
h
n
n
ớ
l
g
n
ứ
n
ả
h
p toàn mẫu với sự giảm nhiệt độ. Số phản ứng Hình 3.18. Số lượng nguyên tử phản ứng và số nguyên tử trong đám không phản ứng lớn nhất. 55 a) b) a) nhiệt độ cao; b) nhiệt độ thấp. Vòng tròn mầu đen biểu diễn nguyên tử mà không có bất kỳ phản ứng xảy ra; vòng tròn màu xanh và màu đỏ biểu diễn nguyên tử mà ở đó phản ứng xảy ra với tần số thấp và tần số cao. Khi nhiệt độ tăng, các phản ứng mở rộng trên nhiều nguyên tử và phân bố đồng nhất hơn. Hình 3.19. Hình minh họa sự phản ứng cho các trường hợp 56 KẾT LUẬN Luận văn đã đạt được một số kết quả chính như sau - Đã tạo ra và khảo sát cấu trúc của ôxít SiO2 lỏng ở nhiệt độ 3200K trong dải áp suất 0-25 GPa. Cấu trúc của SiO2 lỏng được tạo thành từ các khối đa diện SiOx (x = 4 ,5, 6) liên kết với nhau thông qua các nguyên tử cầu ôxy. Ở áp suất thấp, phần lớn các khối đa diện là tứ diện SiO4, liên kết giữa các khối đa diện chủ yếu bởi một nguyên tử cầu ôxy. Khi áp suất tăng, các đa diện liên kết với nhau bởi hai, ba nguyên tử cầu ôxy tăng đáng kể. Nghĩa là khi áp suất tăng, SiO2 lỏng chuyển từ cấu trúc tứ diện sang bát diện. - Khảo sát tính đa thù hình của SiO2 lỏng chúng tôi phát hiện các đơn vị cấu trúc SiOx (x = 4 ,5, 6) không bị thay đổi, nhưng tỉ lệ các đơn vị SiOx (x = 4 ,5, 6) thay đổi mạnh theo áp suất tăng. Kết quả này cho phép xác định mật độ của mẫu SiO2 lỏng ở một trạng thái nào thông qua các đơn vị SiOx (x = 4 ,5, 6). - Mô phỏng cho thấy trong SiO2 chứa một lượng nhỏ các khuyết tật (sai hỏng mạng). Khuyết tật phân bố không đồng nhất trong chất lỏng mà chúng có xu hướng di chuyển lại gần nhau và tạo thành đám. Quá trình khuếch tán thực hiện chủ yếu thông qua phản ứng SiOxSiOx' với x, x'=4, 5 và OSiy OSiy' với y, y'=2, 3. Phản ứng SiOxSiOx' xảy ra không phải ngẫu nhiên mà có xu hướng xảy ra ở vùng có nguyên tử mà tạo thành đám. Chúng tôi phát hiện trong SiO2 có hai loại đám. Đám thứ nhất gồm các nguyên tử mà không có phản ứng SiOxSiOx' nào xảy ra. Đám thứ hai có nguyên tử ở đó phản ứng SiOxSiOx' xảy ra tuần hoàn. Độ linh động của nguyên tử trong đám không phản ứng nhỏ hơn đáng kể so với đám phản ứng. Kích thước của hai đám và độ linh động của chúng phụ thuộc mạnh vào thời gian mô phỏng và nhiệt độ. Do đó cơ chế khuếch tán trong SiO2 lỏng chủ yếu ở vùng có chứa loại đám thứ 2. 57 CÔNG TRÌNH ĐÃ CÔNG BỐ LIÊN QUAN ĐẾN LUẬN VĂN Giap Thi Thuy Trang, Pham Duc Linh, Pham Xuan Truong and Pham Huu Kien, About Structural Transition of Nickel Metal, International Journal of Innovative Science, Engineering & Technology, Vol. 2 Issue 11, November 2015, trang 132-138. 58 TÀ I LIỆU THAM KHẢ O Tiếng Việt 1. Nguyễn Văn Hồng (2010), Mô phỏng ô xít hai nguyên ở trạng thái vô định hình vào lỏng, Luận án Tiến sỹ, Trường Đại học Bách Khoa Hà Nội. 2. Bù i Thi ̣ Quỳnh Trang (2010), Mô hình mới đối với sự khuếch tán của
nguyên tử trong chất rắn vô định hình, Luận văn thạc sỹ, Trường Đại học
Sư phạm Hà Nội 1. Tiếng Anh 3. Adam, G. and Gibbs, J. H. J. (1965), Chemical Physics 43(1), pp.139 -146. 4. Andrea Trave et al (2002), Pressure-Induced Structural Changes in Liquid SiO2 from Ab Initio Simulations, Physical Review Letters, 89(24), 245504. 5. B.T. Poe, P.F. McMillan, D.C. Rubie, S. Chakraborty, J. Yarger,J.Diefenbacher (1997), Silicon and oxygen self-diffusivities in silicate liquids measured to 15 gigapascals and 2800 Kelvin, Science 276, pp.1245-1248. 6. Brown I. D. (1992), Chemical and steric constraints in inorganic solids, Acta Crystallographica B48, pp. 553-572. 7. C.A. Angell, P.A. Cheeseman, S. Tamaddon (1982), Pressure enhancement of ion mobilities in liquid silicates from computer simulation studies to 800 kilobars, Science 218, pp.883-884. 8. Diane M. Christie and James R. Chelikowsky (2000), Electronic and structural properties of germania polymorphs, Physical Review B62,14703. 9. G.A. Appignanesi, J.A. Rodriguez Fris (2009), Journal of Physics: Condensed. Matter 21 203103. 10. Gernot Kostorz (2001), Phase Transformations in Materials, Wiley-VCH Verlag GmbH, D-69469 Weinheim (Federal Republic of Germany. 59 11. Grimley D. I., Wright A. C. and Sinclair R. N. (1990), Neutron scattering from vitreous silica IV. Time-of-flight diffraction, Journal of Non- Crystalline Solids 119, pp. 49. 12. H. Sillescu, R. Bohmer, G. Diezemann, G. Hinze (2002), Journal on Crystalline Solids 307-310, pp.16-23. 13. Horbach J. (2008), Molecular dynamics computer simulation of amorphous silica under high pressure, Journal of Physics: Condensed. Matter 20 244118. 14. Horbach J. and Kob W. (1999), "Static and dynamic properties of a viscous silica melt", Physical Review Section B, Vol 60 (5), pp.3169-3181 15. J.S. Langer (2014), "Theories of glass formation and the glass transition” Report on Progress in Physics 77, number 4. 16. James Badro, (1997), Theoretical study of a five-coordinated silica polymorph, Physical Review SectionB 56, pp. 5797 - 5806. 17. K. Schmidt-Rohr, H.W. Spiess (1991), Physical Review. Lettets 66 3020. 18. K.L. Ngai (2000), Journal on Crystalline Solids 275, pp.7-51. 19. Konnert J. H. and Karle J (1973), The computation of radial distribution functions for glassy materials, Acta Crystallographica Section A 29, pp.702-710. 20. Kubicki J. D. and Toplis M. J. (2002), Molecular orbital calculations on aluminosilicate tricluster molecules: implications for the structure of aluminosilicate glasses, American Mineralogist 87, pp. 668–678. 21. L.I. Tatarinova (Nauka, Moscow, 1983), in The Structure of Solid Amorphous and Liquid Substances, pp.93. 22. M. .S. Shell, G.D. Pablo, Z.P. Athanassios (2002), Molecular structural order and anomalies in liquid silica, Physical Review E 66, 011202. 23. M. Micoulaut, X. Yuan, L.W. Hobbs (2007), Co-ordination and intermediate-range order alterations in densified germania, Journal of Non-Crystalline Solids 353, pp. 1961-1965. 60 24. M.T. Cicerone, F.R. Blackburn, M.D. Ediger (1995), Journal of Chemical Physics102, pp. 471-479. 25. Micoulaut M, Cormier L and Henderson G. S. (2006), The structure of amorphous, crystalline and liquid GeO2, Journal of Physics: Condensed Matter 18, pp.753-784. 26. Mishima, O. H., Calvert, L. D., Whalley, E. (1984), ‘Melting ice’ I at 77 K and 10 Kbar: a new method of making amorphous solids, Nature 310, pp.393-395. 27. Mozzi R. L. and Warren B. E. (1969), The structure of vitreous silica, Journal of Applied Crystallography 2, pp.164-172. 28. P K Hung, N V Hong, and L T Vinh (2007), Journal of Physics: Condensed Matter 19, 466103. 29. Robert L. Erikson and Charles J. Hostetler (1987), "Application of empirical ionic models to SiO2 liquid: Potential model approximations and integration of SiO2 polymorph data", Geochimica et Cosmochimica Acta, Vol. 51, pp.1209-1218 30. S. Tsuneyuki, Y. Matsui (1995), First-Principles Interatomic Potential of Silica Applied to Molecular Dynamics, Physical Review Letters 74, pp.3197. 31. Warren B. E. (1932), Detemination of the structure of glass, the Journal of American Chemical Society 54, pp. 3841. 32. Wolf, G.H, Durben, D.J, McMillan, P.F (1990), High-pressure spectroscopic study of sodium tetrasilicate (Na2Si4O9) glass, Journal of Chemical Physics 93, pp. 2280–2288. 61Trong chương này, trước tiên chúng tôi khảo sát cấu trúc của SiO2 lỏng
Hình 3.1. Hàm phân bố xuyên tâm thành phần của SiO2 lỏng
ở các áp suất khác nhau và T=3200K
Bảng 3.1. Các đặc tính cấu trúc của SiO2 lỏng
(a)
(b)
Hình 3.2. Mạng cấu trúc của mẫu SiO2 lỏng ở 3200 K
Hình 3.3. Các đơn vị cấu trúc
Hình 3.4. Sự phụ thuộc của tỷ lệ các đơn vị cấu trúc SiO4, SiO5 và SiO6 vào áp suất
của mẫu SiO2 ở nhiệt độ 3200K.
Bảng 3.2. Phân bố liên kết cầu O giữa hai đơn vị cấu trúc SiOx với m
là số nguyên tử O tham gia liên kết cầu giữa hai đơn vị cấu trúc SiOx lân cận.
Các cột tiếp theo chỉ ra tỷ lệ phần trăm liên kết cầu tương ứng với m.
Ví dụ, 16,20% số liên kết giữa hai đơn vị cấu trúc lân cận có hai nguyên tử O
tham gia cầu liên kết ở áp suất nén 9,83 GPa.
Áp suất (GPa)
m
0,10
4,87
9,83
15,73
20,15
25,20
Hình 3.5. Phân bố khoảng cách liên kết trong
Hình 3.6. Phân bố góc liên kết O-Si-O trong các đơn vị cấu trúc
(a) là cấu trúc SiO4, (b) là cấu trúc SiO5, (c) là cấu trúc SiO6
Hình 3.7. Phân bố góc liên kết Si-O-Si giữa các đơn vị cấu trúc
Hình 3.8. Phân bố góc Si-O-Si trong các đơn vị cấu trúc
Hình 3.9. Sự phụ thuộc của mật độ vào áp suất của mẫu SiO2 lỏng
ở nhiệt độ 3200K
a)
b)
c)
d)
e)
Hình 3.10. Hình minh họa các đơn vị cấu trúc
Bảng 3.3. Tỷ lệ các phản ứng SiOx SiOx' và OSiy OSiy' trong SiO2 lỏng
theo nhiệt độ 2600K, 3000K, 3200K và 3500K.
Nhiệt độ
2600 K
3000 K
3200 K
3500 K
Hình 3.12. Số lượng đám và số lượng nguyên tử trong đám lớn nhất đối với hai
trường hợp: 1/ Các ô phối trí được chọn ngẫu nhiên trong cấu trúc mạng;
2/ Các ô sai hỏng phát hiện trong hệ. Nhiệt độ là 2600 K.
Hình 3.13. Sự phụ thuộc số lượng các nguyên tử trong đám lớn nhất
theo bước ĐLHPT. Nhiệt độ là 2600 K.
Hình 3.14. Số đám không phản ứng và số lượng nguyên tử trong đám
lớn nhất cho các khoảng thời gian khác nhau. Nhiệt độ là 2600 K.
Hình 3.15. Số đám phản ứng và số lượng nguyên tử trong đám lớn nhất
theo các khoảng thời gian khác nhau. Nhiệt độ là 2600 K.
Hình 3.16. Số lượng nguyên tử trong đám lớn nhất và khoảng cách dịch chuyển
bình phương trung bình theo khoảng thời gian khác nhau, mỗi khoảng thời gian
là 2000 bước. Nhiệt độ là 2600 K
Hình 3.17. Số lượng nguyên tử trong đám lớn nhất và khoảng cách dịch chuyển
bình phương trung bình theo khoảng thời gian khác nhau, mỗi
khoảng thời gian là 5000 bước. Nhiệt độ là 2600 K.