ĐẠ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ử được mô tả bởi công thức:

= a(t) + 6D.t (1.11)

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 sẽ tiệm cận với đường thẳng. Dựa vào độ dốc của đường

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 đạt tới đường tiệm cận sẽ tăng lên. Nếu hệ số khuếch tán

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 SiOxSiOx' 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à lm

(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 SiOxSiOx' 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 SiO4SiO5, SiO5SiO4 và liên kết Si12- O2 được thay thế bởi liên kết Si12- O6;

B) là phản ứng SiO5SiO4, SiO4SiO5 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 SiO4SiO5, SiO5SiO4 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

SiO5SiO4, SiO4SiO5 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 đó là độ dịch chuyển bình phương trung bình của nguyên

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 SiOxSiOx' với x, x'=4, 5 và OSiy OSiy' với

y, y'=2, 3. Phản ứng SiOxSiOx' 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 SiOxSiOx' nào xảy ra. Đám thứ

hai có nguyên tử ở đó phản ứng SiOxSiOx' 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 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

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

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

Để đả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 Å.

Bảng 3.1. Các đặc tính cấu trúc của SiO2 lỏng

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)

(b)

Hình 3.2. Mạng cấu trúc của mẫu SiO2 lỏng ở 3200 K

(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).

Hình 3.3. 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 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.

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.

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).

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

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 Å.

Hình 3.5. Phân bố khoảng cách liên kết trong

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.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 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.

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

)

(

% m ă r t n ầ h p ệ l ỷ T

Góc(độ)

Hình 3.8. Phân bố góc Si-O-Si trong các đơn vị cấu trú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. Sự phụ thuộc của mật độ vào áp suất của mẫu SiO2 lỏng

ở nhiệt độ 3200K

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 SiOxSiOx' nào xảy ra. Đám thứ hai có các nguyên tử ở đó

phản ứng SiOxSiOx' 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)

b)

c)

d)

e)

Hình 3.10. Hình minh họa các đơn vị cấu trúc

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.

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

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

SiOxSiOx' 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.

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

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.

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

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

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.

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

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

.

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

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

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

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

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 SiOxSiOx' với x, x'=4, 5 và OSiy OSiy' với

y, y'=2, 3. Phản ứng SiOxSiOx' 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

SiOxSiOx' nào xảy ra. Đám thứ hai có nguyên tử ở đó phản ứng SiOxSiOx'

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.

61