Phương pháp đa cực nhanh trong tính toán trường điện từ

Sinh viên: Nguyễn Trung Thành – Lớp KSTN – ĐTVT – K52 1

Lời nói đầu

Sự phát triển nhanh chóng của công nghệ thông tin, đặc biệt là sự xuất hiện của các hệ
thống siêu máy tính có tốc độ tính toán nhanh đã mở ra một phƣơng pháp mới trong
nghiên cứu khoa học. Các kỹ thuật sử dụng máy vi tính nhằm giải quyết các vấn đề vật lý
thay cho các phƣơng pháp cũ đã rất phát triển. Trong lĩnh vực nghiên cứu điện từ trƣờng,
sử dụng phƣơng pháp đa cực nhanh (FMM) trên máy tính nhằm giải quyết nhanh các bài
toán dữ liệu đầu vào khổng lồ đã đóng một vai trò vƣợt trội hơn so với các phƣơng pháp
khác.
Phƣơng pháp đa cực nhanh là một kỹ thuật toán học đƣợc phát triển để tăng tốc độ tính
toán trong vấn đề N-body. Thuật toán thực hiện bằng cách mở rộng hàm Green sử dụng
sự mở rộng đa cực, trong đó cho phép một nhóm các nguồn nằm gần nhau và đối xử với
họ nhƣ thể họ là một nguồn duy nhất.
FMM cũng đã đƣợc áp dụng trong việc tăng tốc độ tính toán của phƣơng pháp Moment
(MOM) cũng nhƣ áp dụng cho vấn đề tính toán trƣờng điện từ. Năm 1985, FMM lần đầu
tiên đƣợc giới thiệu bởi Greengard và Rokhlin, dựa trên việc mở rộng đa cực của các
phƣơng trình Helmholtz dạng vector. Nếu FMM đƣợc áp dụng một cách có thứ bậc, nó có
thể cải thiện sự phức tạp của tích vector ma trận từ O(N
2
) xuống O(N). FMM cũng đã
đƣợc áp dụng một cách hiệu quả trong các bài toán tƣơng tác Coulomb hay động lực học
phân tử. FMM đƣợc coi là một trong những thuật toán hàng đầu của thế kỷ 20
Mục đích chính của đồ án này là nghiên cứu các bƣớc thực hiện phƣơng pháp đa cực
nhanh, từ đó sử dụng phần mềm Matlab tính toán các thông số về thời gian, tốc độ tính
toán trong bài toán tƣơng tác trƣờng điện từ.
Em xin cảm ơn các thầy cô giáo trong khoa Điện tử - Viễn thông, cũng nhƣ các thầy cô
giáo tại trƣờng Đại học Bách Khoa Hà Nội đã truyền đạt những kiến thức quý báu, cần

Phương pháp đa cực nhanh trong tính toán trường điện từ

Sinh viên: Nguyễn Trung Thành – Lớp KSTN – ĐTVT – K52 4

Abstract

The purpose of my thesis is research about Fast Multipole Method (FMM) in calculating
electromagnetics field. It present programmes to compute parameters about the execution
time, calculating speed…with Matlab. Namely, this thesis introduces a brief review of
electromagnetics (Maxwell equations, Green’s function…), the steps for calculating
electromagnetics by FMM. At the end, it presents programmes to compute and model: the
potential, the gradient of the potential, the execution time, calculating speed …with
Matlab and analyses the results by using different methods: FMM and classic method.
My thesis includes 3 chapters, in which:
Chapter 1: Basical theory about Electromagnetics
Chapter 2: Theory about Fast Multipole Method
Chapter 3: Computing and modeling with Matlab


3.1.2: Các thông số đầu ra 51
3.2 Xây dựng và tính toán cho mô hình điểm nguồn và đích nằm trong hình tròn 51
3.2.1 Giá trị các tham số đầu vào. 51
3.2.1Giá trị các tham số đầu ra 53
Phương pháp đa cực nhanh trong tính toán trường điện từ

Sinh viên: Nguyễn Trung Thành – Lớp KSTN – ĐTVT – K52 6

3.3 Xây dựng và tính toán cho mô hình điểm nguồn và đích nằm trong hình vuông. 60
3.3.1Giá trị các tham số đầu vào. 60
3.3.2Giá trị các tham số đầu ra. 60
Kết luận 67
Tài liệu tham khảo: 68

Danh mục các hình vẽ
Chƣơng 1
Hình 1.1 Điều kiện biên. 16
Hình 1.2 Miền khối bao quanh bởi mặt kín. 19
Hình 1.3 Hai loại hàm giới hạn. 26
Hình 1.4 Trục toạ độ của hàm giới hạn cho không gian 2 chiều. 29

Chƣơng 2
Hình 2. 1 Hai tập hợp hạt đủ xa trên mặt phẳng 35
Hình 2. 2 Dịch chuyển tâm của khai triển đa cực. 38
Hình 2. 3 Ý tƣởng tính lực xấp xỉ trong FMM 40
Hình 2. 4 Phân mảnh không gian và chỉ số thứ tự Morton, L=2 41
Hình 2. 5 Ví dụ về phân mảnh không gian và đánh số hộp, L=3 41
Hình 2. 6 Một vài mức phân chia trong FMM 42
Hình 2. 7 Các miền xếp theo thứ bậc không gian 44
Hình 2. 8 Pha M2M trong thuật toán FMM. 45

Bảng 3. 9: Giá trị điện thế tại 1 số điểm nguồn trong mô hình 2 61
Bảng 3. 10: Giá trị điện thế tại 1 số điểm đích trong mô hình 2 63
Bảng 3. 11: Thời gian tính toán khi thay đổi số điểm nguồn trong mô hình 2 64
Bảng 3. 12: Tốc độ tính toán khi thay đồi số điểm nguồn trong mô hình 2 65
Phương pháp đa cực nhanh trong tính toán trường điện từ

Sinh viên: Nguyễn Trung Thành – Lớp KSTN – ĐTVT – K52 8

Danh sách các từ viết tắt

FMM
Fast multipole method
Phƣơng pháp đa cực nhanh
SP
Space partition
Phân mảnh không gian
UP
Upward pass
Pha lên
DP
Downward pass
Pha xuống
M2L
Multipole to local
Chuyển đồi đa cực- địa phƣơng
M2M
Multipole to multipole
Chuyển đồi đa cực- đa cực
L2L
Local to local

0 B

(1.3)
ρ D

(1.4)
HB
0



(1.5)
ED
0



(1.6)
Đối với các vật liệu dẫn điện, định luật bảo toàn điện tích đƣợc biểu diễn bởi quan hệ:

(1.7)
Mật độ dòng
J
và cƣờng độ điện trƣờng
E
liên hệ với nhau bởi định luật Ôm :
EJ




J
là mật độ dòng
điện dẫn, ρ là mật độ điện tích. Cuối cùng, ε
0
, μ
0
là hệ số điện môi và hệ số từ thẩm trong
không gian tự do,γ là hệ số phụ thuộc tính dẫn điện của môi trƣờng. Đối với vật liệu điện môi và vật liệu từ thì véc tơ
phân cực
P
và véc tơ từ hoá
M
đƣợc định nghĩa:
HBM
EDP
0


0
1



(1.10)
Các phƣơng trình Maxwell đƣợc biểu diễn lại nhƣ sau:
t


4 .10 [H/m]

















 M
P
J
E
B
tt
000

Phương pháp đa cực nhanh trong tính toán trường điện từ

Sinh viên: Nguyễn Trung Thành – Lớp KSTN – ĐTVT – K52 11




(1.15)
Đây là phƣơng trình sóng của
E
. Tƣơng tự, ta thu đƣợc phƣơng trình sóng của
H
từ
phƣơng trình (1.1):
0
2
2
2




t
H
H


(1.16)
Khi giải phƣơng trình Maxwell để thuận tiện ta định nghĩa hàm vô hƣớng φ và hàm vector
thế
A
thoả mãn :




A





2
2
2
t

(1.21)
Do
E
và φ là các hàm tuỳ ý nên ta có thể chọn chúng sao cho :
t




A

(1.22)
Trong môi trƣờng suy hao thì phƣơng trình sóng sử dụng dạng sau :













t
t
2
2
2
(1.25)
Phƣơng trình (1.23) và (1.24) đƣợc dùng để tính toán sóng bức xạ của anten, trƣờng tán
xạ của vật liệu và sự truyền sóng trong ống dẫn sóng hay các thiết bị điện tử khác.
Trường cân bằng:
Khi bài toán đƣợc xét trong điều kiện trƣờng biến đổi theo thời gian rất chậm thì trạng
thái cân bằng xấp xỉ đƣợc sử dụng. Tiêu chuẩn đƣợc gọi là chậm nếu nó thoả mãn điều
kiện sau:




ω là tần số của tín hiệu hình sin.
Tiêu chuẩn này có nghĩa rằng dòng dẫn chiếm ƣu thế và dòng dịch có thể đƣợc bỏ qua.
Do đó, từ trƣờng xoáy sinh ra bởi điện trƣờng không tồn tại. Không có mối liên hệ giữa
sự thay đổi vị trí và biến đổi theo thời gian của trƣờng. Vì vậy không có sự truyền sóng.
Thông thƣờng, trong các bài toán trƣờng cân bằng đại lƣợng

E
và
H
tuân theo phƣơng trình truyền parabol:
EE
HH


j
j


2
2

(1.30)
Trong trƣờng hợp nhƣ vậy để thuận tiện ta giả thiết sự tồn tại của véc tơ từ thế
A
và véc
tơ điện thế
T
. Việc xác định
A
và
T
xuất phát trực tiếp từ hệ phƣơng trình Maxwell
0 B
và
0 J
:












1

(1.33)
0
1












jj TT


0J 0B JH 
(1.36)
Dựa vào phƣơng trình
00  B ,E
cả điện thế, từ thế vô hƣớng φ, φ
m
và véc tơ từ
thế
T
đƣợc biểu diễn dƣới dạng :
m





H
AB
E
(1.37)
Từ phƣơng trình (1.35), (1.36), và (1.37) Ta thu đƣợc hệ phƣơng trình Poisson và Laplace
:
Phương pháp đa cực nhanh trong tính toán trường điện từ

Sinh viên: Nguyễn Trung Thành – Lớp KSTN – ĐTVT – K52 15

0
2
2
2

n









21
21
21
21
21
0
0
JJ .
KHH .
EE
DD .
BB .
(1.39)
Trong đó
n
là véc tơ chuẩn đơn vị của bề mặt trong hình 1.1,
,
1
E
,


1 1 1
,,
  

2 2 2
,,
  

+
-


n
K
Phương pháp đa cực nhanh trong tính toán trường điện từ

Sinh viên: Nguyễn Trung Thành – Lớp KSTN – ĐTVT – K52 16

Hình 1.1 Điều kiện biên.

Nếu véc tơ thế điện vô hƣớng đƣợc coi nhƣ là một biến thì điều kiện biên giữa 2 mặt là :









là 3 thành phần của
A
,
t
và
s
là 2 véc tơ đơn vị trực giao với véc
tơ chuẩn hƣớng
n
. Với tiêu chuẩn Coulomb thành phần chuẩn
n
A
thoả mãn:
nn 21
AA 

Tính liên tục của thành phần tiếp tuyến của cƣờng độ từ trƣờng
H
đƣợc biểu diễn bằng
biểu thức sau:

































 AAA
1
2
1
2
2
1
1









 AAA
1
2
1
2
2
1
1




(1.44)










trực tiếp từ lý thuyết phân rã.


 S
Sdd AA
(1.45)
Trong đó Ω là miền bao quanh bởi mặt kín S.
A
là hàm véc tơ của vị trí. Giả thiết u và v
là 2 hàm vô hƣớng bất kỳ của vị trí. Nếu u, v cùng vi phân bậc nhất, bậc hai của chúng
liên tục trong không gian Ω và trên bề mặt S thì ta có đƣợc biến đổi sau theo định lý
Divergence:


 S
nSdvudvu )()(
(1.46)
n
là chuẩn véc tơ đơn vị ngoài của bề mặt S. Với S = Sn, Sử dụng phƣơng trình véc tơ:

vuvuvu 
2
)(
(1.47)
và chú ý rằng:
Phương pháp đa cực nhanh trong tính toán trường điện từ

Sinh viên: Nguyễn Trung Thành – Lớp KSTN – ĐTVT – K52 18

n

(1.50)
Phƣơng trình này là dạng đối xứng của (1.49) Nếu trừ vế với vế của (1.49) và (1.50) ta
đƣợc phƣơng trình sau :













 S
dS
n
u
v
n
v
uduvvu )(
22
(1.51)
(Đây là định lý Green dạng thứ 2 thƣờng đƣợc biết đến nhƣ định lý Green).
Nó là định lý tích phân bao gồm Gradient của hàm bị tích. Ý nghĩa của định lý
Divergence là chuyển từ tích phân khối thành tích phân mặt.
Trong trƣờng hợp đặc biệt nếu để u = v và u là nghiệm của phƣơng trình Laplace thì





 S
dS
RnnR
d
R
r
r
11
4
1)(
4
1
)(
,






(1.53)
Đây là phƣơng trình tích phân của thế φ(r) nó không đƣa ra nghiệm của thế. Trong
phƣơng trình này ρ(r’) là mật độ điện tích khối, R = |r– r’|, Ω là khối đƣợc bao quanh bởi
mặt kín S chỉ ra trong hình 1.2.
Phương pháp đa cực nhanh trong tính toán trường điện từ


dS
RnnR
r
11
4
1
)(




(1.54)
Nhƣ vậy tích phân bề mặt trong phƣơng trình (1.53) biểu diễn sự phân bố của nguồn
ngoài mặt S. Hay nói cách khác điều kiện biên biểu diễn sự phân bố của nguồn ngoài mặt
S. Kết luận này ngụ ý rằng điều kiện biên có thể được biểu diễn bởi một nguồn ngoài
tương đương. Nếu không có nguồn ngoài thì tích phân mặt sẽ biến mất. Kiểm tra lại
phƣơng trình (1.51) chỉ ra rằng nếu hàm u là một hàm điều hoà (ví dụ là nghiệm của
phƣơng trình Laplace ) và v = 1 thì dạng 2 của định lý Green sẽ rút gọn thành :

0



S
dS
n
u
(1.55)
z
x

1.4.2 Sự tương đương vector định lý Green

Trong mục này ta đƣa ra phƣơng trình tích phân trong đó véc tơ thế
A
đƣợc coi nhƣ chƣa
biết. Giả thiết
P
và
Q
là hàm véc tơ liên tục của vị trí trong khối kín Ω bao quanh bởi
mặt S cả
P
,
Q
và vi phân từng phần bậc nhất, bậc hai của nó đều tồn tại trên mặt S và
khối Ω. Sử dụng định lý Divergence ta có:
   


 S
dSd QPQP
(1.57)
Phân tích tích phân khối sử dụng biến đổi sau:
 
BAABBA 
(1.58)
Ta thu đƣợc:
   

























11
4
1)'(
4
)(
0
AA
AJ




(1.62)
Do đó:


 drrrfrf )',()()'(

(1.63)
trong đó
)'()'()'()',( zzyyxxrr 


Với phƣơng trình toán tử bất kỳ
(1.64)
Hàm u có thể đƣợc biểu diễn nhƣ sau:


 drrrfru )',()'()(

-1
L
(1.65)
Trong đó L biểu diễn một toán tử bất kỳ (ví dụ với phƣơng trình Poisson
2
L
) và L
-1
là

L
(1.67)
Nếu
L
thì

2 , ,
( , ) ( , )G r r r r

  
(1.68)
Lấy tích phân cả hai vế của phƣơng trình (1.66) thì đƣợc một phƣơng trình biểu diễn toán
tử ngƣợc L
-1
nhƣ sau :

1 1 , ,
( , ) ( , )r r d G r r d



    

LL
(1.69)
Nhƣ vậy toán tử ngƣợc của toán tử vi phân là một toán tử tích phân trong đó hàm Green
là trung tâm. Tuy nhiên hàm Green không đƣợcđịnh nghĩa phù hợp với phƣơng trình
(1.67). Nếu có một hàmg(r) bất kỳ nào thoả mãn Lg = 0 cộng vớiG(r,r’)
~
,,

,
, , , ,
( , ) ( )
( ) ( , ) ( , ) ( ) ( ) ( , )
S
G r r u r
u r r r G r r f r d u r G r r dS
nn





     





(1.69)
Kết hợp (1.63) với (1.69) ta đƣợc phƣơng trình sau:
,
, , ,
( , ) ( )
( ) ( , ) ( ) ( ) ( , )
S
G r r u r
u r G r r f r d u r G r r dS
nn


Phƣơng trình này làđịnh lý Green dạng thứ ba. Trong phƣơng trình (1.71)bao gồm cả giá
trị biên củau và
/un
trên bề mặtS.
Với điều kiện Dirichlet cho
,
( , ) | 0
S
G r r 
thì phƣơng trình (1.71) rút gọn còn:
, , ,
( ) ( , ) ( ) ( ) ( , )
S
u r G r r f r d u r G r r dS
n


  


(1.72)
Với điều kiện Neumann,
,
( , )/ | 0
S
G r r n  
phƣơng trình rút gọn là:

, , ,
( ) ( , ) ( ) ( , )

12
( , )
( ) ( ) ( , ) 0
S
G r r
f r f r G r r
n







Thế phƣơng trình (1.74) vào phƣơng trình (1.72):
, , ,
3
1
1
( ) ( , ) ( ) ( ) ( , )
S
u r G r r f r d f r G r r dS
f

  

(1.75)
hay
,
,,

r G r r d






(1.77)
Phƣơng trình này cho thấy nếu hàm Green của toán tửđƣa ra làđã biết thì nghiệm của bất
kỳ nguồn phân bố nào cũng có thểđƣợc tính toán từ phƣơng trình (1.77). Do đó hàm
Green là công cụ cơ bảnđể phân tích các vấnđề toán lý khác nhau.
Hàm Green của phƣơng trình Poisson trong không gian 3 chiều xét với không gian tự do
là

,
,
0
1
( , )
4
G r r
rr



(1.78)
Đây là nghiệm của phƣơng trình Poisson với nguồn xung đơn vị nếu xétđếnảnh hƣởng
của mặtđất thì

,

là toán tử bất kỳ. Chúý rằng sự khác nhau giữa hàm Green và nghiệm cơ
bảnđó là hàm Green thì liên quan đếnđiều kiện biên nhƣng nghiệm cơ bản thìđƣợc
xácđinh trong không gian không tự do không có biên. Có thể thay thế hàm Green trong
không gian tự do là nghiệm cơ bản của cùng một toán tử. Nghiệm cơ bản của phƣơng
trìnhLaplace trong không gian 2 chiều và 3 chiều có thểđƣợc giải nhƣ sau:
Trong hệ toạđộ cực 2 chiều phƣơng trìnhLaplaceđƣợc khai triển thành:

1
0
d du
r
r dr dr




(1.81)
Phương pháp đa cực nhanh trong tính toán trường điện từ

Sinh viên: Nguyễn Trung Thành – Lớp KSTN – ĐTVT – K52 25

thì nghiệm là

12
lnu C r C
(1.82)
Trong không gian 3 chiều nghiệm của
2
0
d du

R


(1.84)
Phƣong trình cơ bản
Nghiệm cơ bản
2-D
3-D
Phƣơng trình Laplace
2
0F

  

,
11
ln
2
F
rr




,
11
4
F
rr


1
exp -
4
F jk r r
rr




Phƣơng trình khuyếch tán
2
1
( ) ( ) 0
1
F
F r t
kt
k



   



2
,
1
exp
44



Phƣơng trình sóng
2
22
2
2
( ) ( ) 0
1
F
v F r t
t
v



   





,
2
3/ 2 ,
(
2
H vt r r
F
v v r r