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

Tạp chí Khoa học và Giáo dục, Trường Đại học Sư phạm Huế ISSN 1859-1612, Số 01(21)/2012: tr. 5-11 TÌM NGHIỆM CỦA PHƯƠNG TRÌNH POISSON BA CHIỀU BẰNG PHƯƠNG PHÁP CGS LƯƠNG VĨNH GIANG - ĐINH NHƯ THẢO Trường Đại học Sư phạm – Đại học Huế 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 để giải phương trình Poisson ba chiều trong chương trình mô phỏng linh 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ự hợp. Các kết quả thu được cho thấy chương trình giải phương trình Poisson dựa trên thuật toán CGS tuy có độ chính xác không cao nhưng thời gian mô phỏng được cải thiện đáng kể so với thuật toán BICGSTAB(3). 1. GIỚI THIỆU 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 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 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 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 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 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 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ạ đượ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 dựa trên các phương trình chuyển động Newton [3]. 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 đ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 đượ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 rất lớn, do đó cần dùng các giải thuật đặc biệt để giải [3]. Đến nay đã có nhiều phương pháp được xây dựng như: Jacobi, Gauss – Seidel, SOR, đa ô 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 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 con Krylov như CG, GMRES, CGS, QMR, BiCG, BICGSTAB, BICGSTAB(3), 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ệ 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 pháp BICGSTAB, BICGSTAB(3), BICGSTAB tiền điều kiện để giải phương trình 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 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 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 những phương pháp tối ưu, hoạt động ổn định, cho kết quả nhanh hơn. 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 có cấu hình Intel(R) Pentium(R) D CPU 2.80 GHz, DDRII 1GB. Kết quả chỉ ra rằng chương trình mô phỏng sử dụng thuật toán CGS tối ưu hơn. LƯƠNG VĨNH GIANG - ĐINH NHƯ THẢO 6 2. GIẢI PHƯƠNG TRÌNH POISSON BA CHIỀU BẰNG THUẬT TOÁN CGS Phương trình Poisson ba chiều trong trường hợp vật liệu đồng nhất có dạng 2 2 2 2 2 2 , Sx y z ϕ ϕ ϕ ρ ε ∂ ∂ ∂ + + = − ∂ ∂ ∂ ở đâ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; 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 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 , , 2 1, , , 1, , , 1 , , 1, , , 1, , , 16 , i j k i j k i j k i j k i j k i j k i j k i j k S x ρ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ε− − − + + + + + − + + + = − Δ ở đây 1, xi N= , 1, yj N= , 1, zk N= với xN , yN , zN lần lượt là số nút lưới theo các 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 mà ta cần giải. Thuật toán CGS để tìm nghiệm phương trình Poisson được khai triển trong Bảng 1. Bảng 1. Thuật toán CGS để tìm nghiệm của phương trình Poisson Thuật toán CGS 1 2 3 4 5 6 7 8 9 10 Chọn 0x , 0ˆr bất kỳ Tính 0 0r b Ax= -­‐ Lấy 0 0rˆ r= 0 1r = 0 0 0p q= = Bắt đầu vòng lặp chính 1 ˆi ir r Aua-­‐= -­‐ ( )0 1ˆ ,i ir rr -­‐= 1/i ib r r -­‐= 1 1i iu r qb-­‐ -­‐= + ( )1 1i i ip u q pb b-­‐ -­‐= + + 11 12 13 14 15 16 17 18 19 iv Ap= ( )0ˆ/ ,i r va r= iq u va= -­‐ ˆ iu u q= + 1 ˆi ix x ua-­‐= + Nếu ix đủ chính xác thì thoát 1 ˆi ir r Aua-­‐= -­‐ Kết thúc vòng lặp chính Kết thúc thuật toán (1) (2) TÌM NGHIỆM CỦA PHƯƠNG TRÌNH POISSON BA CHIỀU BẰNG THUẬT TOÁN CGS 7 3. KẾT QUẢ MÔ PHỎNG VÀ THẢO LUẬN 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 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 có độ dày tương ứng là id , pd và nd . Mật độ pha tạp acceptor và donor tương ứng là AN và DN , 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 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 bằng mô phỏng thời gian thực trước khi chiếu xung laser vào linh kiện. Chúng tôi đã sử dụng phương pháp Monte – Carlo tập hợp tự hợp để mô phỏng động lực học của hạt tải trong linh kiện trong trường hợp chiếu một xung laser với chiều dài của xung là 12 sf và năng lượng photon là 1.49 eV . Các tham số cấu trúc vùng năng lượng được sử dụng như sau: 1.42gapE eV Γ = , * 00.063em mΓ = , * 00.45hm mΓ = , * 00.222Lem m= , và độ chênh lệch năng lượng giữa Γ và L là 0.29LE eVΓΔ = . Chúng tôi giả thiết rằng 50p nd d nm= = , 340id nm= và 17 30.5 10AN cm −= × , 17 32.5 10DN cm −= × , 16 35 10exN cm −= × sau thời gian 1 ps . Kích thước theo ba chiều không gian của đi-ốt là 440 100 100x y zL L L nm nm nm× × = × × , giả sử đi-ốt được nuôi cấy theo phương Ox . 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 lưới là 1050 10x y z m−Δ = Δ = Δ = × . Như vậy theo phương Ox ta sẽ có 89xN = nút lưới, theo phương Oy và Oz có 21y zN N= = nút lưới. Điện trường ngoài được đặt 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. 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 các phương Ox , Oy , Oz ứng với điện trường ngoài 100extE kV cm= , được tính toán 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ử 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 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 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 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 quá vận tốc [1], [3]. Hình 1. Mô hình đi-ốt p-i-n GaAs LƯƠNG VĨNH GIANG - ĐINH NHƯ THẢO 8 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á trị điện trường ngoài ex 70 /tE kV cm= và 100extE kV cm= . Kết quả cho thấy điện 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 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ý 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 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 Γ 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 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 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. 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ô 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 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 đơ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 dẫn GaAs (chương trình sử dụng thuật toán BICGSTAB(3) đã được chứng minh là cho kết quả phù hợp với thực nghiệm [7]). 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 phương Ox và Oy tại mặt cắt 10z nm= ứng với điện trường ngoài 100extE kV cm= . Đồ thị cho thấy điện thế phân bố đều trong linh kiện và chỉ biến thiên nhỏ xung quanh 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 . 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 và cũng phù hợp với các kết quả đã được công bố trước đây [1], [7]. Hình 2. Vận tốc trôi dạt toàn phần của điện tử và vận tốc trôi dạt theo các phương khác nhau như là hàm của thời gian ứng với 100extE kV cm= Hình 3. Vận tốc trôi dạt toàn phần của điện tử như là hàm của thời gian ứng với các điện trường ngoài khác nhau TÌM NGHIỆM CỦA PHƯƠNG TRÌNH POISSON BA CHIỀU BẰNG THUẬT TOÁN CGS 9 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 = 10 nm ứng với hai thuật toán BICGSTAB(3) và CGS dưới điện trường ngoài 70extE 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ự 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. 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 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ử 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 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 định số và chi phí thời gian tính toán. 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 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 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 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 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 BICGSTAB(3) chạy trong khoảng thời gian 61 phút trong khi chương trình mô phỏng 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ử 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à tối ưu hơn. Hình 4. Vận tốc trôi dạt toàn phần của điện tử theo thời gian ứng với hai thuật toán BICGSTAB(3) và CGS dưới điện trường ngoài 70extE kV cm= 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 tại mặt cắt 10z nm= ứng với 100extE kV cm= LƯƠNG VĨNH GIANG - ĐINH NHƯ THẢO 10 4. KẾT LUẬN Chúng tôi đã 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 trong mô phỏng linh kiện nano bán dẫn bằng phương pháp Monte – Carlo tập hợp tự hợp. Để khảo sát các đặc trưng của phương pháp chúng tôi đã áp dụng để mô phỏng động lực học ba chiều của hạt tải trong đi-ốt p-i-n bán dẫn GaAs và so sánh với các kết quả mô phỏng đã được công bố trước đây. Kết quả cho thấy chương trình giải phương trình Poisson sử dụng thuật toán CGS có độ hội tụ chậm hơn chương trình sử dụng thuật toán BICGSTAB(3) nhưng cho kết quả ổn định hơn và thời gian mô phỏng được cải thiện đáng kể so với chương trình sử dụng thuật toán BICGSTAB(3). LỜI CẢM ƠN Các tác giả xin gửi lời cảm ơn chân thành tới PGS. TS. Nguyễn Hồng Quang, Ban Hợp tác Quốc tế, Viện Khoa học và Công nghệ Việt Nam, chủ nhiệm đề tài NAFOSTED (103.02.109.09) đã hỗ trợ mọi mặt cho các tác giả trong quá trình hoàn thành bài báo này. TÀI LIỆU THAM KHẢO [1] D. N. Thao, S. Katayama, and K. Tomizawa (2004). Numerical simulation of THz radiation by coherent LO phonons in GaAs p-i-n diodes under high electric fields. Journal of the Physical Society of Japan 73, 3177 – 3181. [2] G. Klatt et al. (2011). Photo-Dember terahertz emitter excited with an Er: fiber laser. Appl. Phys. Lett. 98, 021114 – 021114-3. [3] K. Tomizawa (1993). Numerical simulation of submicron semiconductor devices. Artech House, Boston London. Hình 6. Phân bố điện thế không gian trong đi-ốt p-i-n bán dẫn GaAs theo phương Ox tại mặt cắt y = 10z nm= ứng với hai thuật toán BICGSTAB(3) và CGS dưới điện trường ngoài 70extE kV cm= Hình 7. Sự phụ thuộc của chuẩn Euclid của vectơ thặng dư vào số vòng lặp của chương trình con Poisson ứng với 100extE kV cm= TÌM NGHIỆM CỦA PHƯƠNG TRÌNH POISSON BA CHIỀU BẰNG THUẬT TOÁN CGS 11 [4] H. A. Vorst (2003). Iterative Krylov methods for large linear systems. Cambridge University. [5] Y. Saad (2000). Iterative Methods for Sparse Linear Systems. SLAM. [6] G. Speyer, D. Vasileska and S. M. Goodnick (2001). Efficient Poisson equation solvers for large scale 3D simulations. Technical Proceedings of the 2001 International Conference on Modeling and Simulation of Microsystems, Nanotech 2001, Vol. 1, 23-26. [7] D. N. Thao and L. H. Hai (2010). 3D simulation of semiconductor devices using BICGSTAB (3) for the solution of Poisson’s equation. Journal of Science and Education 15, 19-26. [8] D. N. Thao and N. T. Ngoc (2010). 3D simulation of semiconductor devices using preconditioned BICGSTAB algorithm with Jacobi preconditioner for the solution of Poisson’s equation. Journal of Science and Education, 16, 34-41. [9] D. N. Thao, D. T. D. My, N. C. P. Thi and N. T. Thuy (2011)., Three-dimensional simulation of nano semiconductor devices using GPBICG algorithm for the solution of the Poisson's equation. Journal of Science, 65, 215-223. Title: OBTAINING THE SOLUTION OF THE 3D POISSON'S EQUATION BY MEANS OF THE CGS METHOD Abstract: In this article, we present the CGS algorithm used to solve Poisson’s equation in three-dimensional simulation in GaAs p-i-n diodes semiconductor by the self-consitent ensemble Monte Carlo method. These results suggest that CGS algorithm solve 3D Poisson’s equation less accuracy than BICGSTAB(3) algorithm but the simulation time is significantly improved. TS. ĐINH NHƯ THẢO Khoa Vật lý, Trường Đại học Sư phạm – Đại học Huế ĐT: 099-686-7668. Email: dnthao@gmail.com LƯƠNG VĨNH GIANG Học viên Cao học, Trường Đại học Sư phạm – Đại học Huế ĐT: 0974.193.736