Tìm nghiệm của phương trình Poisson ba chiều bằng phương pháp CGS

TÌM NGHIỆM CỦA PHƯƠNG TRÌNH POISSON BA CHIỀU<br /> BẰNG PHƯƠNG PHÁP CGS<br /> LƯƠNG VĨNH GIANG - ĐINH NHƯ THẢO<br /> Trường Đại học Sư phạm – Đại học Huế<br /> Tóm tắt: Trong bài báo này chúng tôi trình bày việc sử dụng thuật toán CGS<br /> để giải phương trình Poisson ba chiều trong chương trình mô phỏng linh<br /> kiện đi-ốt p-i-n bán dẫn GaAs bằng phương pháp Monte – Carlo tập hợp tự<br /> hợp. Các kết quả thu được cho thấy chương trình giải phương trình Poisson<br /> dựa trên thuật toán CGS tuy có độ chính xác không cao nhưng thời gian mô<br /> phỏng được cải thiện đáng kể so với thuật toán BICGSTAB(3).<br /> <br /> 1. GIỚI THIỆU<br /> Trong nghiên cứu khoa học, việc nghiên cứu thực nghiệm và lý thuyết cùng phát triển<br /> song song với nhau. Trong một số trường hợp, việc nghiên cứu lý thuyết có phần thuận<br /> lợi hơn, chẳng hạn trong vấn đề nghiên cứu các linh kiện nano, bước đầu nghiên cứu<br /> thực nghiệm sẽ gặp nhiều khó khăn và tốn kém nên người ta chọn phương pháp nghiên<br /> cứu lý thuyết [1], [2], [3]. Một trong số các phương pháp được chọn ở đây là phương<br /> pháp mô phỏng Monte – Carlo tập hợp tự hợp vì nó có nhiều ưu điểm nổi trội, đặc biệt<br /> là tính chính xác và tính ổn định số. Đây là phương pháp bán cổ điển với tốc độ tán xạ<br /> được tính toán dựa trên quy tắc vàng Fermi và việc khảo sát động lực học của hạt tải<br /> dựa trên các phương trình chuyển động Newton [3].<br /> Một trong những bài toán cơ bản cần giải quyết là bài toán hạt tải chuyển động trong<br /> điện trường. Để biết điện trường ta cần biết điện thế, phân bố điện thế trong linh kiện<br /> được tìm bằng cách giải phương trình Poisson. Đây là một hệ phương trình tuyến tính<br /> rất lớn, do đó cần dùng các giải thuật đặc biệt để giải [3].<br /> Đến nay đã có nhiều phương pháp được xây dựng như: Jacobi, Gauss – Seidel, SOR, đa<br /> ô lưới (multigrid), iLU [4]. Các phương pháp này có thể cho các kết quả chính xác tuy<br /> nhiên độ ổn định số không cao và độ hội tụ thấp. Gần đây, các phương pháp không gian<br /> con Krylov như CG, GMRES, CGS, QMR, BiCG, BICGSTAB, BICGSTAB(3),<br /> BICGSTAB tiền điều kiện đã được phát triển và sử dụng hiệu quả trong việc giải các hệ<br /> phương trình tuyến tính thưa loại lớn [4], [5]. Một số tác giả đã sử dụng các phương<br /> pháp BICGSTAB, BICGSTAB(3), BICGSTAB tiền điều kiện để giải phương trình<br /> Poisson và đã thu được các kết quả chính xác với thời gian tính toán được rút ngắn<br /> nhiều lần [6], [7], [8], [9]. Để khảo sát mở rộng hướng nghiên cứu trên chúng tôi tiến<br /> hành tìm nghiệm của phương trình Poisson bằng phương pháp CGS với mục đích tìm ra<br /> những phương pháp tối ưu, hoạt động ổn định, cho kết quả nhanh hơn.<br /> Chúng tôi đã xây dựng chương trình mô phỏng mới và thực hiện tính toán trên máy tính<br /> có cấu hình Intel(R) Pentium(R) D CPU 2.80 GHz, DDRII 1GB. Kết quả chỉ ra rằng<br /> chương trình mô phỏng sử dụng thuật toán CGS tối ưu hơn.<br /> Tạp chí Khoa học và Giáo dục, Trường Đại học Sư phạm Huế<br /> ISSN 1859-1612, Số 01(21)/2012: tr. 5-11<br /> <br /> 6<br /> <br /> LƯƠNG VĨNH GIANG - ĐINH NHƯ THẢO<br /> <br /> 2. GIẢI PHƯƠNG TRÌNH POISSON BA CHIỀU BẰNG THUẬT TOÁN CGS<br /> Phương trình Poisson ba chiều trong trường hợp vật liệu đồng nhất có dạng<br /> <br /> ∂ 2ϕ ∂ 2ϕ ∂ 2ϕ<br /> ρ<br /> + 2 + 2 =− ,<br /> 2<br /> ∂x<br /> ∂y<br /> ∂z<br /> εS<br /> <br /> (1)<br /> <br /> ở đây ϕ là điện thế, ρ là mật độ điện tích, ε S là hằng số điện môi tĩnh trong linh kiện;<br /> x , y , z là ba biến không gian. Thực hiện phép sai phân hữu hạn và chia mô hình linh<br /> kiện thành các ô lưới bằng nhau Δx = Δy = Δz thì phương trình (1) được viết lại như sau<br /> <br /> ϕi −1, j ,k + ϕi , j −1,k + ϕi , j ,k −1 − 6ϕi , j ,k + ϕi +1, j ,k + ϕi , j +1,k + ϕi , j ,k +1 = −<br /> <br /> ρi , j , k 2<br /> Δx ,<br /> εS<br /> <br /> (2)<br /> <br /> ở đây i = 1, N x , j = 1, N y , k = 1, N z với N x , N y , N z lần lượt là số nút lưới theo các<br /> chiều không gian Ox , Oy , Oz . Đây chính là hệ phương trình tuyến tính thưa cực lớn<br /> mà ta cần giải. Thuật toán CGS để tìm nghiệm phương trình Poisson được khai triển<br /> trong Bảng 1.<br /> Bảng 1. Thuật toán CGS để tìm nghiệm của phương trình Poisson<br /> <br /> Thuật toán CGS<br /> 1<br /> <br /> Chọn x0 , rˆ0 bất kỳ<br /> <br /> 11<br /> <br /> v = Api<br /> <br /> 2<br /> <br /> Tính r0 = b -­‐ Ax0<br /> <br /> 12<br /> <br /> a = r i / (rˆ0 , v)<br /> <br /> 3<br /> <br /> Lấy rˆ0 = r0<br /> <br /> 13<br /> <br /> qi = u -­‐ a v<br /> <br /> 14<br /> <br /> uˆ = u + qi<br /> <br /> 15<br /> <br /> xi = xi-­‐ 1 + a uˆ<br /> <br /> 4<br /> 5<br /> <br /> r0 = 1<br /> p0 = q0 = 0<br /> <br /> 6<br /> 7<br /> <br /> Bắt đầu vòng lặp chính<br /> ri = ri-­‐ 1 -­‐ a Auˆ r i = (rˆ0 , ri-­‐ 1 )<br /> <br /> 16<br /> <br /> Nếu xi đủ chính xác thì thoát<br /> <br /> 8<br /> <br /> 17<br /> <br /> ri = ri-­‐ 1 -­‐ a Auˆ<br /> <br /> 9<br /> <br /> b = r i / r i-­‐ 1<br /> u = ri-­‐ 1 + b qi-­‐ 1<br /> <br /> 18<br /> <br /> Kết thúc vòng lặp chính<br /> <br /> 10<br /> <br /> pi = u + b (qi-­‐ 1 + b pi-­‐ 1 )<br /> <br /> 19<br /> <br /> Kết thúc thuật toán<br /> <br /> TÌM NGHIỆM CỦA PHƯƠNG TRÌNH POISSON BA CHIỀU BẰNG THUẬT TOÁN CGS<br /> <br /> 7<br /> <br /> 3. KẾT QUẢ MÔ PHỎNG VÀ THẢO LUẬN<br /> Mô hình cấu trúc của đi-ốt p-i-n bán dẫn GaAs gồm một lớp bán dẫn thuần i kẹp giữa<br /> hai lớp bán dẫn pha tạp loại p và loại n như được chỉ ra trong Hình 1, trong đó mỗi lớp<br /> có độ dày tương ứng là di , d p và d n . Mật độ pha tạp acceptor và donor tương ứng là<br /> <br /> N A và N D , các tạp được phân bố từ bề mặt của các lớp p và n vào sâu bên trong linh<br /> kiện theo hàm phân bố Gauss. Trạng thái cân bằng nhiệt của linh kiện được xác lập<br /> bằng mô phỏng thời gian thực trước khi chiếu xung laser vào linh kiện.<br /> Chúng tôi đã sử dụng phương pháp Monte – Carlo tập hợp tự hợp để mô phỏng động<br /> lực học của hạt tải trong linh<br /> kiện trong trường hợp chiếu<br /> một xung laser với chiều dài<br /> của xung là 12 fs và năng<br /> lượng photon là 1.49 eV .<br /> Các tham số cấu trúc vùng<br /> năng lượng được sử dụng<br /> Γ<br /> như sau: Egap<br /> = 1.42 eV ,<br /> <br /> mΓ* e = 0.063 m0 ,<br /> mΓ* h = 0.45 m0 ,<br /> <br /> Hình 1. Mô hình đi-ốt p-i-n GaAs<br /> <br /> *<br /> Le<br /> <br /> m = 0.222 m0 , và độ chênh<br /> lệch năng lượng giữa Γ và L là ΔEΓL = 0.29 eV . Chúng tôi giả thiết rằng<br /> d p = d n = 50 nm ,<br /> <br /> di = 340 nm<br /> <br /> và<br /> <br /> N A = 0.5 ×1017 cm −3 ,<br /> <br /> N D = 2.5 × 1017 cm −3 ,<br /> <br /> N ex = 5 × 1016 cm −3 sau thời gian 1 ps . Kích thước theo ba chiều không gian của đi-ốt là<br /> Lx × Ly × Lz = 440 nm ×100 nm ×100 nm , giả sử đi-ốt được nuôi cấy theo phương Ox .<br /> Mô hình linh kiện được chia thành các ô lưới không gian với khoảng cách giữa các nút<br /> lưới là Δx = Δy = Δz = 50 ×10−10 m . Như vậy theo phương Ox ta sẽ có N x = 89 nút<br /> lưới, theo phương Oy và Oz có N y = N z = 21 nút lưới. Điện trường ngoài được đặt<br /> vào linh kiện dọc theo phương Ox và đi-ốt được phân cực nghịch, xem Hình 1.<br /> Hình 2 mô tả sự thay đổi vận tốc trôi dạt toàn phần và vận tốc trôi dạt của điện tử theo<br /> các phương Ox , Oy , Oz ứng với điện trường ngoài Eext = 100 kV cm , được tính toán<br /> với lời giải phương trình Poisson bằng thuật toán CGS. Từ đồ thị ta thấy rằng điện tử<br /> chủ yếu chuyển động trôi dạt theo phương Ox còn các thành phần vận tốc theo phương<br /> Oy , Oz đóng góp không đáng kể. Đặc biệt, sau khi chiếu xung laser, vận tốc trôi dạt<br /> toàn phần của điện tử (vận tốc trôi dạt của điện tử theo phương Ox ) tăng nhanh vượt xa<br /> giá trị bão hòa sau đó giảm nhanh về giá trị bão hòa, hiện tượng này được gọi là sự vượt<br /> quá vận tốc [1], [3].<br /> <br /> 8<br /> <br /> LƯƠNG VĨNH GIANG - ĐINH NHƯ THẢO<br /> <br /> Hình 2. Vận tốc trôi dạt toàn phần của điện<br /> tử và vận tốc trôi dạt theo các phương khác<br /> nhau như là hàm của thời gian ứng với<br /> <br /> Hình 3. Vận tốc trôi dạt toàn phần của điện tử<br /> như là hàm của thời gian ứng với các điện<br /> trường ngoài khác nhau<br /> <br /> Eext = 100 kV cm<br /> Hình 3 mô tả sự phụ thuộc của vận tốc trôi dạt của điện tử theo thời gian ứng với các giá<br /> trị điện trường ngoài Eext = 70kV / cm và Eext = 100 kV cm . Kết quả cho thấy điện<br /> trường ngoài càng cao thì sự vượt quá vận tốc xảy ra càng sớm, đỉnh vượt quá vận tốc<br /> càng cao và vận tốc càng nhanh chóng về tiệm cận giá trị bão hòa. Điều này được lý<br /> giải, khi điện trường ngoài càng cao thì số điện tử nằm trong các trạng thái có thể tham<br /> gia vào quá trình tán xạ liên thung lũng càng lớn. Khi điện tử bị tán xạ từ thung lũng Γ<br /> sang thung lũng L vận tốc của điện tử bị giảm nhiều do khối lượng hiệu dụng của điện<br /> tử trong thung lũng L lớn hơn nhiều lần khối lượng hiệu dụng của điện tử trong thung<br /> lũng Γ. Hệ quả là vận tốc của điện tử càng giảm nhanh về giá trị bão hòa.<br /> Hình 4 cho kết quả so sánh vận tốc của điện tử thu được bằng hai chương trình mô<br /> phỏng ba chiều sử dụng thuật toán BICGSTAB(3) và thuật toán CGS [7]. Hai đồ thị gần<br /> như trùng nhau hoàn toàn. Điều đó nói lên sự thành công trong việc sử dụng thuật toán<br /> đơn giản CGS để giải bài toán chuyển động của hạt tải trong linh kiện đi-ốt p-i-n bán<br /> dẫn GaAs (chương trình sử dụng thuật toán BICGSTAB(3) đã được chứng minh là cho<br /> kết quả phù hợp với thực nghiệm [7]).<br /> Hình 5 mô tả sự phân bố điện thế không gian trong đi-ốt p-i-n bán dẫn GaAs theo các<br /> phương Ox và Oy tại mặt cắt z = 10nm ứng với điện trường ngoài Eext = 100 kV cm .<br /> Đồ thị cho thấy điện thế phân bố đều trong linh kiện và chỉ biến thiên nhỏ xung quanh<br /> một số các điểm lưới theo phương Ox và gần như không đổi theo hai phương Oy , Oz .<br /> Kết quả hoàn toàn hợp lý do điện trường ngoài được đặt vào linh kiện theo phương Ox<br /> và cũng phù hợp với các kết quả đã được công bố trước đây [1], [7].<br /> <br /> TÌM NGHIỆM CỦA PHƯƠNG TRÌNH POISSON BA CHIỀU BẰNG THUẬT TOÁN CGS<br /> <br /> Hình 4. Vận tốc trôi dạt toàn phần của điện<br /> tử theo thời gian ứng với hai thuật toán<br /> BICGSTAB(3) và CGS dưới điện trường<br /> ngoài Eext = 70 kV cm<br /> <br /> 9<br /> <br /> Hình 5. Phân bố điện thế không gian trong điốt p-i-n bán dẫn GaAs theo hai phương Ox, Oy<br /> tại mặt cắt z = 10nm ứng với<br /> <br /> Eext = 100 kV cm<br /> <br /> Hình 6 chỉ ra việc so sánh sự phân bố điện thế trong không gian theo trục Ox tại y = z =<br /> 10 nm ứng với hai thuật toán BICGSTAB(3) và CGS dưới điện trường ngoài<br /> Eext = 70 kV cm . Ta thấy hai đường này trùng khít nhau, qua đó một lần nữa nói lên sự<br /> thành công trong việc sử dụng thuật toán CGS để giải phương trình Poisson ba chiều.<br /> Các kết quả thu được nói lên sự thành công trong việc sử dụng thuật toán CGS, đến đây<br /> ta nhận thấy với bài toán tìm nghiệm của phương trình Poisson ba chiều thì việc sử<br /> dụng thuật toán CGS có ưu điểm hơn các thuật toán khác bởi tính đơn giản của nó. Tuy<br /> nhiên, một vấn đề không kém phần quan trọng cần được thảo luận là độ hội tụ, tính ổn<br /> định số và chi phí thời gian tính toán.<br /> Hình 7 so sánh độ hội tụ của thuật toán BICGSTAB(3) và thuật toán CGS. Kết quả cho<br /> thấy tuy độ hội tụ của thuật toán BICGSTAB(3) nhanh hơn nhưng đồ thị tương ứng với<br /> thuật toán CGS trơn hơn điều này nói lên tính ổn định số của thuật toán CGS. Khi tiến<br /> hành chạy chương trình mô phỏng trên máy tính có cấu hình Intel(R) Pentium(R) D<br /> CPU 2.80 GHz, DDRII 1 GB với sai số 10-12 thì chương trình mô phỏng với thuật toán<br /> BICGSTAB(3) chạy trong khoảng thời gian 61 phút trong khi chương trình mô phỏng<br /> với thuật toán CGS chạy trong khoảng thời gian 42 phút. Qua đó chúng ta thấy việc sử<br /> dụng thuật toán CGS đối với bài toán tìm nghiệm của phương trình Poisson ba chiều là<br /> tối ưu hơn.<br /> <br />