Ứng dụng phương pháp FDTD 2 chiều trong mô phỏng trường điện từ

Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số 31, 06 - 2014 109 øng dông ph­¬ng ph¸p fdtd 2 chiÒu trong m« pháng tr­êng ®iÖn tõ NGUYỄN HUY HOÀNG, NGUYỄN VĂN TRUNG, NGUYỄN THÙY LINH Tóm tắt: Bài báo giới thiệu tóm tắt việc ứng dụng phương pháp sai phân hữu hạn trong miền thời gian (Finite difference Time Domain - FDTD) hai chiều trong mô phỏng trường điện từ với các nội dung chính: Trình bày tóm tắt các vấn đề về rời rạc hóa các phương trình Macxoen bằng phương pháp FDTD và điều kiện biên hấp thụ trong mô phỏng 2 chiều hay còn gọi là lớp hấp thụ (Perfect Matched Layer - PML); trên cơ sở đó tiến hành mô phỏng 2 chiều với mô hình sóng điện từ phẳng và đưa ra các nhận xét từ kết quả mô phỏng. Các chương trình mô phỏng được thực hiện trên phần mềm Matlab và kết quả mô phỏng thu được phù hợp với lý thuyết trường điện từ. Từ khóa: FDTD, Trường điện từ, Phương trình Macxoen, PML, Mô phỏng. 1. MỞ ĐẦU Phương pháp FDTD được Kane Yee người Nhật giới thiệu vào năm 1966, phương pháp này được đưa ra nhằm mục đích giải trực tiếp bằng số các phương trình Macxoen trong các môi trường và các miền không gian khác nhau trong miền thời gian. Trong phương pháp này, điện trường và từ trường được rời rạc hóa trong phép lấy vi phân các phương trình Macxoen theo phương pháp sai phân trung tâm và sau đó các giá trị rời rạc của chúng sẽ được tính toán bằng máy tính. Trong các phương pháp được sử dụng để tính toán số và mô phỏng trường điện từ như phương pháp mô men, phương pháp phần tử hữu hạn, phương pháp FDTD thì phương pháp FDTD được sử dụng phổ biến hơn cả vì nó cho phép giải quyết được số lượng lớn các bài toán điện từ, đặc biệt là các bài toán liên quan đến các vật thể có cấu trúc phức tạp (2D và 3D) hay các bài toán có liên quan đến các vật thể có kích thước so sánh được với bước sóng cũng như các bài toán yêu cầu miền tần số cần khảo sát lớn. Với các ưu điểm như vậy, phương pháp FDTD hiện là một công cụ rất mạnh mẽ và hữu hiệu được ứng dụng rộng rải khi giải các bài toán phức tạp liên quan đến điện từ trường trong nhiều lĩnh vực như thiết kế anten, kỹ thuật siêu cao tần, radar... Cơ sở lý thuyết và nội dung của phương pháp FDTD đã được trình bày chi tiết trong [4, 7], bài báo này tập trung vào nghiên cứu ứng dụng phương pháp FDTD để mô phỏng 2 chiều sóng điện từ phẳng trong cả 2 trường hợp không và có thiết lập lớp hấp thụ PML. 2. RỜI RẠC HÓA PHƯƠNG TRÌNH MACXOEN VÀ LỚP HẤP THỤ 2.1 Rời rạc hóa phương trình Macxoen z Trong mô phỏng hai chiều ta chọn một trong hai nhóm gồm 3 vectơ sau để mô phỏng [6]: Trường từ ngang (TM), gồm các thành phần yxz HHE ,, ~ hoặc trường điện ngang (TE) gồm các thành phần yxz EEH ~ , ~ , . Ta sẽ làm việc với trường TM. Với trường TM, phương trình Macxoen [1] được tiến hành rời rạc hóa tương tự như trường hợp mô phỏng một chiều bằng phương pháp FDTD như sau: 1/2 1/2 0 0 ( 1/ 2, ) ( 1/ 2, ) ( , 1/ 2) ( , 1/ 2)( , ) ( , ) 1 (1) n n n nn n y y x xz z H i j H i j H i j H i jD i j D i j t x y                  1 1/2 1/2 0 0 ( , 1/ 2) ( , 1/ 2) ( , 1) ( , )1 (2) n n n n x x z zH i j H i j E i j E i j t y             Khoa hoc máy tính N. H. Hoµng, N. V. Trung, N. T. Linh,“øng dông ph­¬ng ph¸p FDTD tr­êng ®iÖn tõ.” 110 1 1/2 1/2 0 0 ( 1/ 2, ) ( 1/ 2, ) ( 1, ) ( , )1 (3) n n n n y y z z H i j H i j E i j E i j t x            ở đây, ta vẫn sử dụng điều kiện ổn định nghiệm như trong mô phỏng một chiều [4, 7], đó là: 0.2/ cxt  , với c0 là vận tốc ánh sáng trong chân không và để đơn giản ta chọn xy  . Các phương trình rời rạc hóa (1), (2), (3) sẽ được sử dụng để viết mã mô phỏng FDTD 2 chiều. 2.2. Lớp hấp thụ Giả sử ta đang mô phỏng một sóng điện từ phát ra từ một nguồn điểm truyền trong không gian tự do. Sóng sẽ lan truyền ra xung quanh, đến biên của vùng không gian mô phỏng và bị phản xạ lại không như ý tại biên, lúc đó ta sẽ không xác định được đâu là sóng thực và đâu là sóng phản xạ. Đây chính là lý do phải xác lập điều kiện biên hấp thụ trong phương pháp FDTD. Một trong những điều kiện biên hấp thụ hiệu quả và linh hoạt nhất trong mô phỏng FDTD 2 chiều đó là tầng phối hợp trở kháng hay lớp hấp thụ (PML) do Berenger xây dựng. Ý tưởng cơ bản để xây dựng nên PML là, khi sóng truyền từ môi trường A sang môi trường B thì hệ số phản xạ của sóng sẽ phụ thuộc vào trở sóng của hai môi trường đó [2, 3, 6]: CA CB CA CB Z Z R Z Z    Trở sóng của môi trường lại phụ thuộc vào hằng số điện môi và độ từ thẩm của môi trường theo công thức: CZ   Nếu ta cho  thay đổi theo  sao cho trở sóng của hai môi trường luôn bằng nhau thì sẽ không xảy ra hiện tượng phản xạ tại mặt phân cách, hay hệ số phản xạ bằng 0. Tiếp theo ta sẽ tạo ra một môi trường tổn hao sao cho khi sóng truyền vào sẽ bị suy giảm hoàn toàn trước khi tới được biên. Điều này được thực hiện khi chúng ta coi  và  là các số phức vì phần ảo của chúng thể hiện sự suy giảm như đã biết từ lý thuyết trường điện từ [2]. Sau một vài phép biến đổi bằng cách sử dụng thêm các hằng số  và  ảo, ký hiệu là *** ,, FyFxFz  , hệ phương trình Macxoen [1] sẽ có dạng như sau: * * 0 . ( ). ( ) . (4) ( ) '*( ). ( ) y x z Fz Fz z z H H i D x y c x y D E                 (5) * * 0 . ( ). ( ) . (6) z x Fx Fx E i H x y c y        * * 0 . ( ). ( ) . (7) z y Fy Fy E i H x y c x       Các hằng số ảo *F và * F trong các công thức trên phụ thuộc vào cả hai hướng x, y và chúng không ảnh hưởng đến các hằng số thực của môi trường. Đây là các đại lượng phức và chúng có dạng như sau: * * 0 0 ;Dm HmFm Fm Fm Fm i i             Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số 31, 06 - 2014 111 trong đó, Dm và Hm là các tham số thỏa mãn điều kiện [6]: 0 0 0 Hm Dm D        với D là độ dẫn điện của môi trường đang xét; m = x hoặc y. Ta xây dựng PML theo phương x. Sử dụng điều kiện ổn định nghiệm 2/1/).( 0  xct , phương trình (4) được biến đổi thành: 1/ 2 1/ 2 ( , ) 3( ). ( , ) 2( ).0.5. ( 1 / 2, ) ( 1 / 2, ) ( , 1 / 2) ( , 1 / 2) (8) n n z z n n n n y y x y D i j gi i D i j gi i H i j H i j H i j H i j              Các tham số gi2 và gi3 có dạng sau: 0 0 0 1 2( ) (9) 1 ( ). / (2 . ) 1 ( ). / (2. ) 3( ) (10) 1 ( ). / (2 . ) D D D gi i i t i t gi i i t               Ta biến đổi tương tự đối với phương trình (7): 1/ 2 1/ 2 ( 1 / 2, ) 3( 1 / 2). ( 1 / 2, ) 2( 1 / 2).0.5. ( 1, ) ( , ) (11) n n y y n n z z H i j fi i H i j fi i E i j E i j              trong đó, 0 0 0 1 2( 1 / 2) (12) 1 ( 1 / 2). / (2. ) 1 ( 1 / 2). / (2. ) 3( 1 / 2) (13) 1 ( 1 / 2). / (2. ) D D D fi i i t i t fi i i t                    Phương trình (6) được biến đổi thành: )2/1,( 2 ).(. )2/1,( )2/1,( ).( . .. )2/1,()2/1,( 2/1 0 0 2/1 0 001                 jiI tx rotE x tc jiH jiI tx x tc rotE x tc jiHjiH n Hx Dn x n Hx Dn x n x     Ở đây, ta ký hiệu      T n n Hx x rotE jiI 0 2/1 )2/1,( . Do vậy khi viết chương trình mô phỏng phương trình (6) sẽ được thực hiện thông qua các phương trình dưới đây: 1/2 1/2 1/2 1/2 1 1/2 ( , ) ( , 1) (14) ( , 1/ 2) ( , 1 / 2) (15) ( , 1 / 2) ( , 1 / 2) 0.5. 1( ). ( , 1 / 2) (16) n n z z n n Hx Hx n n n x x Hx rotE E i j E i j I i j I i j rotE H i j H i j rotE fi i I i j                    trong đó, 0 ( ). 1( ) 2 D i tfi i     Giá trị biến thiên của các tham số như sau [6]: fi1(i) : (0 ÷ 0.333); gi2(i) : (1 ÷ 0.75); gi3(i) : (1 ÷ 0.5) Ở trên ta đã xây dựng PML theo phương x. Một cách tương tự ta hoàn toàn có thể xây dựng theo phương y. Ta có phương trình cho Dz: Khoa hoc máy tính N. H. Hoµng, N. V. Trung, N. T. Linh,“øng dông ph­¬ng ph¸p FDTD tr­êng ®iÖn tõ.” 112 1/2 1/2( , ) 3( ). 3( ) ( ) 2( ). 2( ).0.5. ( 1 / 2, ) ( 1 / 2, ) ( , 1 / 2) ( , 1 / 2) (17) n n z z n n n n y y x x D i j gi i gj j D i gi i gj j H i j H i j H i j H i j              và 1/2 1/2 1/2 1/2 ( 1, ) ( , ) (18) ( 1/ 2, ) ( 1/ 2, ) (19) n n z z n n Hy Hy rotE E i j E i j I i j I i j rotE            1 1/2 ( 1 / 2, ) 3( 1 / 2). ( 1 / 2, ) 2( 1 / 2).0.5. 1( ). ( 1 / 2, ) (20) n n y y n Hy H i j fi i H i j fi i rotE fi j I i j            cuối cùng, Hx theo phương x có dạng như sau: 1/2 1/2 1/2 1/2 1 1/2 ( , ) ( , 1) ( , 1/ 2) ( , 1/ 2) ( , 1/ 2) 3( 1/ 2). ( , 1/ 2) 2( 1/ 2).0.5. 1( ). ( , 1/ 2) n n z z n n Hx Hx n n x x n Hx rotE E i j E i j I i j I i j rotE H i j fj j H i j fj j rotE fi i I i j                       Kết hợp tính toán theo cả 2 phương x và y, ta đã có đủ các tham số của môi trường PML với các giá trị biến thiên của các tham số như sau [6]: fi1(i) và fj1(j): (0 ÷ 0.333); fi2(i), gi2(i), fj2(j) và gj2(j) : (1 ÷ 0.75) fi3(i), gi3(i), fj3(j) và gj(3) : (1 ÷ 0.5) Lưu ý rằng ta có thể “tắt” môi trường PML trong miền không gian khảo sát bằng cách cho hai tham số fi1 và fj1 bằng 0 còn các tham số khác bằng 1. 3. MÔ PHỎNG SÓNG ĐIỆN TỪ PHẲNG 3.1. Mô hình sóng điện từ phẳng Mô hình sóng phẳng thường được sử dụng trong mô phỏng trường điện từ [4, 7]. Nhiều bài toán thực tế, ví dụ như tính toán diện tích phản xạ hiệu dụng trong radar..., sử dụng mô hình sóng phẳng. Hơn nữa, tại khoảng cách xa cỡ chục lần bước sóng thì trường bức xạ của hầu hết anten có thể coi là sóng phẳng [2, 3]. Để mô phỏng sóng phẳng trong chương trình FDTD 2 chiều, không gian mô phỏng được chia làm 2 vùng: Vùng trường tổng hợp và vùng trường tán xạ [6]. Vùng trường tổng hợp là vùng mà trong đó giá trị của trường là sự tổng hợp của sóng tới và sóng tán xạ đi ra từ vật thể đặt trong vùng này. Vùng trường tán xạ là vùng mà giá trị của trường chính là của sóng tán xạ ra từ vật thể đặt trong vùng trường tổng hợp. Nếu trong vùng trường tổng hợp không đặt vật thể nào thì trường tổng hợp chính là sóng tới, còn trường tán xạ khi đó không tồn tại. Hình 1. Vùng trường tổng hợp và vùng trường tán xạ trong mô phỏng sóng phẳng. Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số 31, 06 - 2014 113 Hình 1 minh họa cho điều này. Sẽ có một ma trận sóng tới một chiều tạo ra sóng phẳng, một điểm trong ma trận đóng vai trò là nguồn và trường Ez sẽ được tạo ra tại điểm đó. Khi đó sóng sẽ lan truyền về cả hai phía không gian. Vì PML thực tế không lý tưởng, nên trước khi sóng phẳng truyền đến lớp này sẽ bị làm suy giảm để sao cho khi đi vào PML hoàn toàn không xảy ra hiện tượng phản xạ. Hình 2. Vùng trường tổng hợp và vùng trường tán xạ trong mô phỏng sóng phẳng. Minh họa trên hình 2 cho thấy, mỗi điểm trong không gian mô phỏng 2 chiều chỉ có thể nằm trong vùng trường tổng hợp hoặc ngoài vùng đó, không có điểm nào nằm trên biên cả. Do đó, nếu có một điểm trong vùng trường tổng hợp sử dụng các điểm ở ngoài để tính giá trị của nó thì điểm đó phải được tính toán lại và ngược lại. Đây chính là nguyên nhân phải có ma trận sóng tới, nó chứa các giá trị cần thiết để tạo ra sự thay đổi này. Có ba vị trí cần phải tính toán lại đó là: - Giá trị Dz tại j = ja hoặc j = jb: Dz(i, ja) = Dz(i, ja) + 0.5.Hxt( ja – 1/2);Dz(i, jb) = Dz(i, jb) - 0.5.Hxt( jb + 1/2) - Giá trị Hx ngay ở ngoài j = ja và j = jb: Hx(i, ja – 1/2) = Hx(i, ja – 1/2) + 0.5.Ezt( ja);Hx(i, jb + 1/2) = Hx(i, jb + 1/2) - 0.5.Ezt( jb) - Giá trị Hy ngay ở ngoài i = ia và i = ib: Hy(i – 1/2, j) = Hx(i – 1/2, j) - 0.5.Ezt( j); Hy(i + 1/2, j) = Hx(i + 1/2, j) + 0.5.Ezt( j) Với ja, jb và ia, ib lần lượt là giới hạn vùng không gian mô phỏng theo các trục y và x; còn Et và Ht là ký hiệu các thành phần của sóng tới. 3.2. Mô phỏng 2 chiều FDTD khi chưa thiết lập PML PhÇn nµy tr×nh bµy kÕt qu¶ m« pháng 2 chiÒu mét xung Gauss truyÒn trong m«i tr­êng kh«ng khÝ b»ng ph­¬ng ph¸p FDTD khi ch­a thiÕt lËp PML. Hình 3. Xung Gauss được tạo ra tại tâm của không gian mô phỏng tại T=10. Hình 4. Xung lan truyền ra xung quanh tại T=50. Khoa hoc máy tính N. H. Hoµng, N. V. Trung, N. T. Linh,“øng dông ph­¬ng ph¸p FDTD tr­êng ®iÖn tõ.” 114 C¸c h×nh 3 ®Õn h×nh 6 thÓ hiÖn kÕt qu¶ víi c¸c chØ sè b­íc thêi gian T kh¸c nhau, ta thÊy r»ng xung bÞ ph¶n x¹ l¹i vµ giao thoa víi xung tíi khi ®i ra biªn. 3.3. Mô phỏng 2 chiều FDTD có thiết lập PML Phần này trình bày kết quả mô phỏng 2 chiều một xung hình sin có dạng: )1015002sin( 6 Tdt , chiều dài của lớp PML được thiết lập bằng 20 lần kích thước cell, môi trường là không khí. Kết quả mô phỏng trong các hình 7 đến hình 10 với các chỉ số bước thời gian T khác nhau chỉ ra rằng, xung không bị phản xạ lại khi đi ra biên mà bị hấp thụ và suy giảm khi đi vào lớp PML và cuối cùng bị triệt tiêu tại biên. Hình 5. Xung bắt đầu truyền đến biên tại T=64. Hình 6. Hình ảnh xung bị giao thoa do phản xạ nếu không thiết lập PML tại T=100. Hình 9. Xung truyền vào PML và không bị phản xạ lại tại T=70. Hình 10. Xung tiếp tục truyền vào PML và bị hấp thụ tại T=100. Hình 7. Xung được tạo ra tại tâm của vùng không gian mô phỏng tại T=10. Hình 8. Xung truyền ra biên có thiết lập PML tại T=50. Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số 31, 06 - 2014 115 3.4. Mô phỏng 2 chiều FDTD sóng phẳng có thiết lập PML Sóng phẳng 2 chiều được mô phỏng trong phần này có dạng: )104002sin( 6 Tdt , truyền trong không khí, lớp hấp thụ PML cũng được thiết lập với chiều dài bằng 20 lần kích thước cell. Kết quả chỉ ra trong các hình 11 đến hình 14 với các chỉ số bước thời gian T khác nhau, cho thấy sóng bắt đầu hình thành từ 1 biên và lan truyền đến biên kia, khi ra biên không bị phản xạ lại mà bị hấp thụ và suy giảm khi đi vào lớp PML và cuối cùng bị triệt tiêu. 4. KẾT LUẬN Bài báo đã giới thiệu tóm tắt về cách rời rạc hóa các phương trình Macxoen cho mô phỏng FDTD 2 chiều, phương pháp thiết lập PML và mô hình mô phỏng sóng phẳng đơn giản. Việc thiết lập PML chủ yếu dựa vào các tham số f và g, và cách lựa chọn giá trị của chúng phụ thuộc vào từng bài toán mô phỏng cụ thể, việc thiết lập chiều dài của lớp PML cũng đóng một vai trò quan trọng để sao cho khi sóng truyền ra biên không bị phản xạ lại. Mô hình sóng phẳng trong phương pháp FDTD 2 chiều cũng được thực hiện dựa trên sự thiết lập PML và việc phân chia không gian mô phỏng làm 2 vùng là vùng trường tổng hợp và vùng trường tán xạ. Để tạo ra sóng phẳng trong không gian 2 chiều ta chỉ đơn giản sử dụng một ma trận sóng tới 1 chiều. Tất cả kết quả mô phỏng đều được thực hiện trên phần mềm Matlab minh họa cho việc trường điện từ lan truyền trong không gian 2 chiều Hình 11. Sóng bắt đầu hình thành từ một biên tại T=20. Hình 12. Sóng tiếp tục lan truyền đến biên kia tại T=50. Hình 13. Sóng bị triệt tiêu khi đến biên tại T=120. Hình 14. Sóng tiếp tục bị triệt tiêu khi đến biên tại T=150. Khoa hoc máy tính N. H. Hoµng, N. V. Trung, N. T. Linh,“øng dông ph­¬ng ph¸p FDTD tr­êng ®iÖn tõ.” 116 khi có và không có PML, cũng như mô phỏng sự lan truyền sóng phẳng trong không gian 2 chiều bằng phương pháp FDTD. TÀI LIỆU THAM KHẢO [1]. J. C. Maxwell, “A Treatise on Electricity and Magnetism”, Clarendon Press Series (1873). [2]. J. A. Straton, “Electromagnetic Theory”, McGraw-Hill Book Company, Inc. (1941). [3]. C. A.Balanis, “Advanced Engineering Electromagnetics”, John Wiley & Sons Inc. (1989). [4]. C. K. S. Kunz, R. J. Luebbers, “Finite Difference Time Domain Method for Electromagnetics”, CRC Press (1993). [5]. A. Taflove, S. C.Hagness, “Computational Electrodynamics The Finite-Difference Time-Domain Method”, Artech House (2005). [6]. D. M. Sullivan, “Electromagnetic Simulation Using The FDTD Method”, IEEE Press Series on RF and Microwave Technology (2000). [7]. S. D.Gedney, “Introduction to the Finite-Difference Time-Domain (FDTD) Method for Electromagnetic”, Morgan & Claypool Publishers (2001). ABSTRACT APPLICATION OF 2D FINITE-DIFFERENCE TIME DOMAIN (FDTD) METHOD TO ELECTROMAGNETIC FIELD SIMULATION. The paper introduces briefly the application of 2D FDTD to simulating Electromagnetic Field with several main contents: Present concisely the process of discreting Maxwell field equations by using FDTD method and the absorbing bound condition in 2D simulation, also called Perfect Matched Layer (PML); Carry out the 2D simulation with simple plane wave and give some remarks based on the simulated results. Keywords: FDTD, Electromagnetic field, Field Maxwell equation, PML, Simulation. Nhận bài ngày 19 tháng 12 năm 2013 Hoàn thiện ngày 20 tháng 5 năm 2014 Chấp nhận đăng ngày 28 tháng 5 năm 2014 Địa chỉ: Học viện Kỹ thuật quân sự Email: hoangnh@mta.edu.vn, trungcs2@mta.edu.vn, truonglinh0808@gmail.com Di động: 0902146368