Ứ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