Luận văn Mô phỏng cấu trúc và động học của hệ silica lỏng với mô hình kích thước lớn

y, nhưng áp suất cao cấu trúc mạng bao gồm các đơn vị cấu trúc vị SiO4, SiO5 và SiO6. Tính dẻo của đơn vị cấu trúc SiO5 là nguyên nhân tăng cơ tính của silica [49]. Sự giòn, dễ gãy vỡ của SiO2 đã được nghiên cứu bằng cả thực nghiệm [50] và MD [49,50]. Bằng cách sử dụng kính hiển vi nguyên tử và phương pháp MD, những nghiên cứu này cho thấy sự tồn tại của dòng chảy plastic trước khi xảy ra hiện tượng bề mặt tự do bị nứt. Các vết nứt sẽ lan truyền, phát triển và kết tụ tạo thành các hốc sai hỏng ở kích thước nano. Ngược lại, có một số công trình vẫn còn tranh cãi chưa thống nhất được về quan điểm hình thành vết nứt trong hệ. Hơn nữa, bằng mô phỏng người ta đã quan sát thấy trong hệ không đồng nhất động học tồn tại các nguyên tử chuyển động nhanh và các nguyên tử chuyển động chậm, các nguyên tử "chuyển động nhanh" hoặc "chuyển động chậm" tạo thành các đám chuyển động trong không gian [18]. Tính không đồng nhất động học (DH) đã được phát hiện trong các mô hình MD bởi các hàm tương quan hai và bốn Số hóa bởi Trung tâm Học liệu và Công nghệ thông tin – ĐHTN điểm, kỹ thuật trực quan hóa và phân tích đám [19]. Tuy nhiên, các hàm và phương pháp trực quan này dường như không thể nắm bắt được sự khác nhau của vùng “nhanh” và “chậm”. Do đó, các phương pháp làm rõ DH vẫn cần phải được cải thiện để làm rõ cơ chế mức nguyên tử của hiện tượng không đồng nhất động học, mối liên hệ giữa cấu trúc và động học không đồng nhất. Số hóa bởi Trung tâm Học liệu và Công nghệ thông tin – ĐHTN Chương 2 PHƯƠNG PHÁP TÍNH TOÁN Nội dung chương này, chúng tôi trình bày các kỹ thuật mô phỏng hệ ôxít SiO2 lỏng tại các áp suất khác nhau. Cụ thể như sau: các mô hình SiO2 lỏng được tạo ra bằng cách gieo ngẫu nhiên các nguyên tử O và Si (13332 nguyên tử O và 6666 nguyên tử Si) vào khối hình lập phương; kích thước của khối hình lập phương được xác định dựa trên mật độ thực của mẫu SiO2 lỏng. Chúng tôi sử dụng kỹ thuật mô phỏng ĐLHPT để tạo ra các mô hình SiO2 lỏng ở nhiệt độ 3500 K với áp xuất trong dải từ 0 GPa đến 45 GPa. Trên cơ sở các mô hình SiO2 lỏng, các đặc trưng về cấu trúc cũng như các tính chất động học của SiO2 lỏng sẽ được xác định và phân tích. 2.1. Xây dựng mô hình động lực học phân tử 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. Các nguyên tử có tọa độ ban đầu 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. Do chịu 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 các yếu tố nhiệt độ và áp suất. Các nguyên tử trong mô hình chuyển động tuân theo định luật cơ học cổ điển Newton. Đối với hệ này, phương trình chuyển động của định luật hai Newton có thể viết như sau: 2 i i i i i 1 N2 d r m a =m =F(r ...,r ) dt (2.1) 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: 1 N ij i j ij U F r      (2.2) Số hóa bởi Trung tâm Học liệu và Công nghệ thông tin – ĐHTN trong đó, ijU 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. Để tính toán tương tác xa, gần đúng Ewald-Hansen đã được sử dụng. Trong mô phỏng ĐLHPT, có nhiều thuật toán được sử dụng để giải hệ phương trình chuyển động của các nguyên tử theo định luật hai Newton. Trong số các thuật toán được sử dụng hiện nay, thuật toán Verlet được sử dụng rộng rãi hơn do tính đơn giản của nó. 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 i i i i i F(t) r (t + dt) = 2r (t) - r (t - dt) + (dt) m (2.3) Vận tốc của nguyên tử ở 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: i i i r (t+dt)-r (t-dt) v (t) = 2dt (2.4) 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 như sau: i i i ij ij ij y zi x x y z j j j F ( t ) F F F F F F        (2.5) trong đó ijx j F được xác định như sau: 0ij ij i j x j ij ij U( r ) x x F x . . r r           (2.6) với 0x là véctơ đơn vị của trục Ox. Các thành phần ij ijy zF , F được xác định tương tự như ijx F ở phương trình (2.6). Trong quá trình mô phỏng, năng lượng E của hệ có thể tính theo công thức: E=K+U (2.7) trong đó U và K lần lượt là thế năng và động năng của hệ và được tính theo Số hóa bởi Trung tâm Học liệu và Công nghệ thông tin – ĐHTN biểu thức sau: ij ij i>j U= U (r ) (2.8) 2N i i i=1 m v K= 2  (2.9) 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 lựa 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, 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. Đối với các mô hình khác sẽ có các đại lượng tương ứng không thay đổi. 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 3 B T=K Nk (2.10) 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 2 2            2N N i i i i i i=1 i m v m r( t dt ) r ( t dt ) K= dt (2.11) Trong mô hình NVT, người ta thường sử dụng kỹ thuật điều chỉnh nhiệt độ để giữ nhiệt độ có giá trị không đổi. Ý 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.10). Điều chỉnh vận tốc iv của tất cả các nguyên tử theo phương trình sau: Số hóa bởi Trung tâm Học liệu và Công nghệ thông tin – ĐHTN ' 0 i i T v = v T (2.12) Thay vận tốc mới vào công thức (2.10) và (2.11) chúng ta sẽ thu được: N N ' 2 20 i i i i 0 i=1 i=1B B T1 1 T'= m (v ) = m v =T 3k N 3k N T   (2.13) Á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ẽ chọn λ > 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: a a a a a a b b b b b b x' [i] = x [i]. ; y' [i] = y [i]. ; z' [i] = z [i]. x' [i] = x [i]. ; y' [i] = y [i]. ; z' [i] = z [i].       (2.14) 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 N i i i=1B B 2 K(t) 1 T(t)= = m (v (t) ) 3 k N 3Nk  (2.15) 1 3 ij ij i j N P( t ) kT( t ) r ( t )F ( t ) V V     (2.16) 2.2. Thế tương tác đối với hệ SiO2 Đề tài sử dụng thế tương tác cặp Van Beets-Kramer-Van Santen (BKS) để nghiên cứu cấu trúc, tính chất động học của hệ SiO2 lỏng theo áp suất. Thế tương tác này có dạng: 2 6i j ij ij ij ij ij ij ij q q e U( r ) A exp( B r ) C r r     (2.17) Số hóa bởi Trung tâm Học liệu và Công nghệ thông tin – ĐHTN trong đó, rij là khoảng cách giữa hai nguyên tử thứ i và j. Các hệ số thế được liệt kê trong bảng 2.1. Bảng 2.1. Các hệ số trong thế BKS đối với hệ SiO2 Cặp nguyên tử Aij (eV) Bij (Å-1 ) Cij (eV×Å6) Điện tích (e) O-O 1388,773 2,760 175,000 qO = − 1,2 Si-O 18003,757 4,873 33,538 qSi = + 2,4 Si-Si 0,0 0,0 0,0 Thế BKS có một đặc trưng phi vật lý, đó là năng lượng tương tác của các cặp O-O, Si-O tiến tới âm vô cùng khi khoảng cách tương tác giữa các nguyên tử tiến tới 0. Để tránh hiệu ứng phi vật lý ở khoảng cách gần này, chúng tôi đã hiệu chỉnh thế BKS bằng cách cộng thêm thế tương tác Lennard-Jones vào thế BKS. Việc hiệu chỉnh này không làm thay đổi dạng của thế BKS ở khoảng cách xa. Thế BKS ban đầu và thế BKS đã hiệu chỉnh được trình bày trên hình 2.1. Kết quả cho thấy hai thế này hầu như trùng khít nhau ở khoảng cách r > 1,3 Å đối với tương tác Si-O và r > 1,7 Å đối với tương tác O-O. . Hình 2.1. Dạng đồ thị của thế BKS: (1) tương tác O-O; (2) dạng hiệu chỉnh của tương tác O- O; (3) tương tác Si-O và (4) dạng hiệu chỉnh đối với tương tác Si-O 1 2 3 4 5 -5 0 5 10 15 20 4 3 1 2 U (r )( eV ) r (Å) Số hóa bởi Trung tâm Học liệu và Công nghệ thông tin – ĐHTN 2.3. Phương pháp gần đúng Ewald-Hansen Trong mô phỏng, với thế tương tác xa ta thường sử dụng kỹ thuật gần đúng Ewald. 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: 0 1 1 1 2 N N i j n i j ij ,n q q U r       (2.18) 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) = n1Lx + n2Ly +n3Lz 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. Bảng 2.2 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). Trong biểu thức (2.18), 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. ij ,n jn i i jr r r r r nL     (2.19) Bảng 2.2. Biểu diễn cách tính toán gần đúng Ewald –Hansen trong không gian 2 chiều 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) Theo gần đúng Ewald, thế tương tác (2.17) được tách thành ba thành phần trong đó có hai thành phần hội tụ nhanh và một số hạng không đổi như sau: Ewald r m oU U U U   (2.20) Số hóa bởi Trung tâm Học liệu và Công nghệ thông tin – ĐHTN Trong đó, Ur là tổng trong không gian thực, Um là tổng trong không gian đảo và Uo là số hạng không đổi. Các thành phần này có dạng như sau: 1 2 N ij ,n r i j i , j n ij ,n erfc( r ) U q q r    (2.21) 2 2 0 2 1 2 i j N m i j i , j m πm exp πim.( r r ) U q q πV m                (2.22) 2 1 N o i i U q π      (2.23) Trong đó V là thể tích của hình hộp mô phỏng; m = (l,j,k) là véc tơ mạng đảo; n = n1.L.x + n2.L.y+n3.L.z = (n1,n2,n3). Trong biểu thức tổng (2.22) khi n = 0 chỉ xét i j. Hàm sai bù của biến x là erfc(x) được xác định theo phương trình sau:     2 0 2 x uerfc x 1 – erf x 1 e du π     (2.24) Từ (2.24) ta có thể xác định lực tác dụng lên hạt thứ i: FP(i) = FPm(i) + FPr(i) (2.25) với   23 1 2N ij ,n P r i j ij ,n ij ,n ij ,n j , j i n ij ,n ( r ) FP i q q erfc( r ) r exp( ( r ) ) r π                (2.26)   2 2 1 1 0 2 2Ni P m ij j , j m q m m π FP i exp sin m.r L m L L                       (2.27) trong đó P = x,y,z Trong gần đúng Ewald-Hansen, thế tương tác dạng (2.20) được viết dưới dạng:  o E-H r mU U U U   (2.28) trong đó:   1 N o i i o o i Z .Z U U .e L  (2.29) Số hóa bởi Trung tâm Học liệu và Công nghệ thông tin – ĐHTN 1 ij N r i j o i j ij .r erfc L U Z .Z .e r L                      (2.30) 2 4 4 6 6 8 8 N i j m o ij o o i j Z Z U a exp( π )(V V S V V V S .e L           (2.31) với ij ij r L    04 2 4 6o 4 6 8 1V a a . a . a        0 2 4 6 4 4 6 8 1V b b . b . b       0 2 4 6 6 4 8 1 12V C C . C . C      và 4 4 44 x y zS f     2 2 2 6 x y zS      8 8 8 8 x y z S      (2.32) Lực tác dụng lên hạt thứ i gồm hai số hạng được xác định như sau: 2 2 1 1 2 ijij r N ijr (i) i j o j , ji ij ij πrr experfc π r .LLU F Z Z e r r r .L                    (2.33) 2 4 4 6 6 8 82 1 2 r N i j om ( i ) o ji Z Z eU F exp( ).( )(V V S V S V S ) r L               + exp(- ρ2 8 6 6 84 4 4 8 6 4 6 8 oV V V S SV .SS S S V V V                               (2.34) trong đó: Uo, eo, ao, a4, a6, a8, a10, b4, b6, b8, b10, C6, C8, C10, C12, d8, d10, d12 là các hệ số; Zi và Zj lần lượt là điện tích của ion i và j; L là kích thước hình hộp mô phỏng và rij là khoảng cách giữa ion thứ i và ion thứ j. 2.4. Xác định các đặc trưng vi cấu trúc Tiếp theo, chúng tôi trình bày các đại lượng đặc trưng vi cấu trúc của SiO2 lỏng đó là: HPBXT, SPT, phân bố góc liên kết, phân bố độ dài liên kết và trực quan hóa các đơn vị cấu trúc (ĐVCT). Số hóa bởi Trung tâm Học liệu và Công nghệ thông tin – ĐHTN 2.4.1. Hàm phân bố xuyên tâm Trong mô hình vật liệu ở trạng thái lỏng, HPBXT là 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ử. Trong đó HPBXT g(r) được xác định như sau: 2 ij i , j i V g( r ) ( r r ) N     (2.35) 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.35) có thể viết lại một cách tường minh hơn như sau: 1 22 N N N N ij i , j j V g( r ) drdr ...dr P ( r ) ( r r ) N      (2.36) ở đâ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: 24 V ( r , r ) dr g( r ) r rg( r. )    (2.37) Thay phương trình (2.37) vào (2.36) chúng ta có: 1 22 2 22 2 4 4 N N N N ij i , j iV ( r , r ) N N N N ij i , j i V ( r , r ) V g( r ) dr dr dr ...dr p ( r ) ( r r ) r rN V dr ...dr p ( r ) dr ( r r ) r rN                  (2.38) 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à i ij i j V ( r , r ) n ( r, r ) dr ( r r )       (2.39) Thay (2.39) vào (2.37) chúng ta tìm được: 1 22 2 2 2 4 4 N N N i i i i V g( r ) dr dr ...dr P ( r ) n ( r, r ) r rN V n ( r, r ) r rN           (2.40) Số hóa bởi Trung tâm Học liệu và Công nghệ thông tin – ĐHTN Phương trình (2.40) có thể viết lại một cách đơn giản như sau: 0 ( r ) g( r )    (2.41) 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. 0 2 1 4 i i N V n ( r, r ) N ( r ) r r             (2.42) 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ạ. Do tán xạ là đàn hồi, |kin|=|kout|, với k=kin- kout chúng ta có: 4 2 in K= k = sin( / )    (2.43) 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: (θ) (K) (K)I =f NS (2.44) 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.45) và chứa tất cả các thông tin về vị trí của các nguyên tử. N l m l,m 1 S(K)= exp[ iK(r -r )] N  (2.45) Liên hệ giữa thừa số cấu trúc với HPBXT, chúng ta dùng định nghĩa chuẩn về HPBXT (2.35) và biểu diễn chuyển đổi Fourier của hàm Dirac delta như trong phương trình sau: 1 2 iKx( x ) dK e      (2.46) Số hóa bởi Trung tâm Học liệu và Công nghệ thông tin – ĐHTN Tách tổng trong phương trình (2.45) thành 2 phần ứng với l=m và lm 1 1 N N l m l l m N lm l m 1 S(K)= exp[iK.0] exp[iK.(r -r )] N N 1 exp[iK.(r )] N         (2.47) Áp dụng chuyển đổi Fourier trong không gian ba chiều đối với (2.47) chúng ta được: 1 NiK.r lm3 l m 1 dK e S( K ) ( r ) ( r r ) (2π) N       (2.48) Thay phương trình (2.46) vào (2.48) ta được: iK.r 3 1 N dK e S(K)=δ(r)+ g(r) (2π) V (2.49) Giản ước hàm delta trong phương trình (2.49) ta được: iK.r 3 1 N dK e [S(K)-1]= g(r) (2π) V (2.50) Từ phương trình (2.50) 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 được tính trong chương trình mô phỏng như sau:   0 0 1 βα NN , αβ ij i j αβ α β g ( r ) N δ(r -r) α, β Si, O N α=β N ( N ) N N α β ρ N N                (2.51) Ở đâ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 i i j j , ij i i i c b c b g ( r ) g( r ) c b            (2.52) Số hóa bởi Trung tâm Học liệu và Công nghệ thông tin – ĐHTN 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 đỉnh thứ nhất của HPBXT tương ứng: 2 0 4 cr jZ g ( r )r dr   (2.53) trong đó, rc là bán kính ngắt, thường được chọn là vị trí cực tiểu ngay sau đỉnh 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 HPBXT 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). 2.4.3. Xác định phân bố góc Để xác định góc O-T- O hoặc T- O -T, khi biết toạ độ của các nguyên tử tương ứng, chúng tôi làm như sau: Giả sử xét một tập hợp gồm ba nguyên tử với toạ độ tương ứng: O1(x1,y1,z1); T(x2,y2,z2); O2(x3,y3,z3). Góc O-T-O được xác định bằng biểu thức: 1 2 1 2 1 2 2 2 2 2 2 2 1 1 1 2 2 2 l .l m .m n .n arccos l m n . l m n             (2.54) trong đó: 1 2 1 1 2 1 1 2 1 2 2 3 2 2 3 2 2 3 l x x ; m y y ; n z z ; l x x ; m y y ; n z z ;             (2.55) Ở đây, cần chú ý rằng 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 T-O-T được xác định hoàn toàn tương tự như đối với góc O-T-O. Số hóa bởi Trung tâm Học liệu và Công nghệ thông tin – ĐHTN 2.4.4. Trực quan hóa dữ liệu các đơn vị cấu trúc Cấu trúc của Silica lỏng được phân tích với các mẫu SiO2 ở nhiệt độ 3500 K và áp suất 0÷45 GPa. Trong dải áp suất được khảo sát, mạng SiO2 được hình thành từ các ĐVPT SiOx (x = 4, 5, 6) và OSiy (y = 2, 3) (Hình 2.3). Các ĐVPT cơ bản SiOx (x = 4, 5, 6) nghĩa là nguyên tử Si được bao quanh bởi x nguyên tử O ở khoảng cách lân cận gần nhất. Tương tự, OSiy (y = 2, 3) nghĩa là nguyên tử O được bao quanh bởi y nguyên tử Si ở khoảng cách lân cận gần nhất. Các ĐVPT khác chiếm tỉ phần nhỏ và không xét tới ở đây, liên kết Si-O được hình thành bởi Si và O trong khoảng cách nhỏ hơn so với rSi-O. Cấu trúc mạng Si-O bao gồm các nguyên tử được kết nối với nhau thông qua liên kết Si-O. Các đơn vị cấu trúc mô tả sơ đồ trong Hình 2.3. Trong luận văn này, chúng tôi sử dụng các chương trình Matlab để trực quan hóa ảnh 3 chiều các ĐVCT của hệ SiO2 lỏng ở nhiệt độ 3500 K theo áp suất, phân tích chi tiết có thể thấy trong chương 3. Hình 2.2. Mô tả 5 ĐVPT cấu thành mạng SiO2 lỏng SiO4 (a); SiO5 (b); SiO6 (c); OSi2 (d); OSi3(e) và đám các ĐVCT (f). Ở đây, quả cầu nhỏ là O và quả cầu lớn là Si. (f) Số hóa bởi Trung tâm Học liệu và Công nghệ thông tin – ĐHTN 2.5. Phương pháp khảo sát động học trong hệ SiO2 lỏng Tính chất động học của hệ SiO2 lỏng là do quá trình khuếch tán của các nguyên tử Si và O trong hệ gây ra. Một đặc trưng khá quan trọng của quá trình khuếch tán, đó là sự phụ thuộc độ dịch chuyển bình phương trung bình của nguyên tử vào số bước dịch chuyển của mô hình. Trên thực tế, đây chính là sự phụ thuộc của vào thời gian. Trên cơ sở sự phụ thuộc này, chúng ta có thể xác định được hệ số khuếch tán D. Theo lý thuyết khuếch tán, trong mô hình ĐLHPT giá trị trung bình của 2r ( t )  phụ thuộc vào thời gian t theo biểu thức: 2 6r ( t ) Dt  (2.56) Như vậy, nếu xây dựng đồ thị phụ thuộc giữa 2r ( t )  vào t chúng ta sẽ nhận được sự phụ thuộc tuyến tính. Tuy nhiên, để xác định 2r ( t )  sẽ đòi hỏi thực hiện một khối lượng lớn các thực nghiệm mô phỏng và tiêu tốn rất nhiều thời gian. Một phương pháp khác hiệu quả hơn đó là chúng ta tiến hành xác định hàm nội suy tuyến tính: 2r ( t ) at b   (2.57) Theo phương pháp bình phương tối thiểu, chúng ta sẽ xác định được các hệ số a và b. Hệ số a trong biểu thức (2.57) sẽ là hệ số khuếch tán D cần tìm. Số hóa bởi Trung tâm Học liệu và Công nghệ thông tin – ĐHTN Chương 3 KẾT QUẢ VÀ THẢO LUẬN Nội dung 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 45 GPa ở nhiệt độ 3500 K bằng cách phân tích các mẫu SiO2 lỏng kích thước lớn, cụ thể mẫu chứa 19998 nguyên tử. Tiếp theo, 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 của hệ SiO2 lỏng trong khoảng áp suất từ 0 đến 45 GPa ở nhiệt độ 3500 K bằng phương pháp pháp tích độ dịch chuyển bình phương trung bình (DCBPTB) của các nguyên tử. 3.1. Khảo sát cấu trúc SiO2 lỏng theo áp suất Bằng các phương pháp xây dựng mẫu mô phỏng như trình bày ở chương 2, chúng tôi tạo ra các mẫu SiO2 lỏng chứa 19998 nguyên tử, có mật độ lấy từ dữ liệu thực nghiệm, ở nhiệt 3500 K với thế tương tác cặp BKS. Hình 3.1 là hình ảnh ba chiều mô tả sự sắp xếp các nguyên tử Si và O trong mẫu SiO2 lỏng ở áp suất 0 GPa và nhiệt độ 3500 K. Có thể nhận thấy mẫu SiO2 lỏng có cấu trúc trật tự xa, các nguyên tử Si và O phân bố có trật tự trong mẫu. Để khảo sát các đặc trưng cấu trúc cơ bản của SiO2 lỏng ở nhiệt độ 3500 K và áp suất trong khoảng 0 đến 45 GPa, đầu tiên chúng tôi sử dụng phương pháp phân tích các đặc điểm của các HPBXT. Trật tự tầm trung được phân tích thông qua HPBXT gSi-Si(r). Trên hình 3.2 cho thấy, độ cao đỉnh thứ nhất HPBXT gSi- Si(r) giảm đáng kể. Cụ thể là độ cao giảm từ 3,09 đến 2,27 Ǻ trong khoảng áp suất từ 0 đến 15 GPa, và giảm rất ít trong khoảng áp suất từ 15 đến 45 GPa. Vị trí đỉnh thứ nhất HPBXT cặp Si-Si thay đổi từ 3,10 Ǻ ở áp suất 0 GPa đến 3,02 Ǻ ở áp suất 5 GPa và không thay đổi trong khoảng áp suất từ 5 đến 30 GPa và ổn định trong khoảng áp suất từ 30 đến 45 GPa, do đó khoảng cách liên kết trung bình của cặp Si-Si thay đổi không đáng kể khi áp suất tăng từ 0 đến 45 GPa. Số hóa bởi Trung tâm Học liệu và Công nghệ thông tin – ĐHTN Hình 3.1. Ảnh 3D mô tả cấu trúc của hệ SiO2 lỏng ở 3500 K, áp suất phòng. Ở đây, quả cầu màu đỏ là nguyên tử Si; quả cầu màu xám là nguyên tử O Trên hình 3.3 trình bày HPBXT gSi-O của SiO2 lỏng ở nhiệt độ 3500 K và các áp suất khác nhau. Đối với cặp Si-O ở áp suất không, đỉnh thứ nhất HPBXT khá nhọn, chứng tỏ cấu trúc trật tự tầm gần (cấu trúc địa phương) trong SiO2 lỏng