Chapter 2 馬克斯威爾方程組
2.4 卷積完美匹配吸收層
了高吸收效能的完美匹配層 (Perfectly Matched Layer, PML) 邊界條件 [7]，才使得 FDTD 方法煥發出新的生命力。目前，已發展出數種 PML 方法，即 Gedney 所提 出的軸向完美匹配層 (uniaxial PML, UPML) [23]，Abarbanel 所提出的各向異性完 美匹配層 (anisotropic PML, APML) [24] 等方法。本文所使用的吸收邊界稱為卷積 完美匹配層 (convolution PML, CPML) [9]，其優點在於它對於消逝波 (evanescent wave) 具有很好的吸收能力，此外它尚能大幅度減少記憶體的占用。除了吸收效 果較好之外，亦具普適性，CPML 吸收邊界完全獨立於 FDTD 計算域內的介質，
可以不做任何修改地應用到各向同性、各向異性、吸收、色散或是非線性介質的 計算中。因此，CMPL 相對于 UPML 吸收邊界而言，更適用於非自由空間中目標 散射特性的研究。在此，我們給出了其張量係數 (tensor coefficient)sw 的形式 [25]
sw = κw+ σw
(其 中 u(t) 和 δ(t) 分 別 是 單 位 步 階 函 數 (Unit step function) 和 脈 衝 函 數 (Delta 卷積項的計算，將會消耗大量的記憶體和計算時間。透過 Luebbers 和 Hunsberger [26] 所採用的遞推卷積 (RC) 技術，可以有效的解決此一問題。
際的計算中將嚴重地影響計算速度。Luebbers 和 Hunsberger 將式 (2.29) 寫為以下
其中
其中
σw = σw,max
(d− w d
)m
, aw = aw,max
(d− w d
)ma
; 0≤ w ≤ d, kw = 1 + (kw,max− 1) ·
(d− w d
)m
. (2.36)
式 (2.36) 中的 d 代表 CPML 厚度，而 σmax 之大小可以下式計算之
σw,max = 0.8(m + 1)
η0∆√ϵrµr, (2.37) 其中 η0 =√
µ0/ϵ0為真空中的波阻抗。
圖2.1為 CPML 之示意圖，利用圖2.2說明了以理想馬克斯威爾方程為核心概 念，並延伸至具不同材料色散介質特性的應用。CPML 的概念圖如下。
Fig 2.1: CPML 吸收邊界之示意圖。
Fig 2.2: 求解電磁波傳遞問題的方程及計算空間之示意圖。
第 三 章 數值方法
求解電磁波方程的數值方法主要可分為時域和頻域兩大類。頻域法包括有矩 量法 (Moment method, MoM)、有限元素法 (Finite element method, FEM) 等；時域 法包括有時域有限差分法 (Finite difference time domain, FDTD)、時域積分方程法 (Time domain integral equation, TDIE) 等。本論文根據差分理論 (difference theory)，
利用泰勒級數展開式 (Taylor series) 將微分方程轉換成差分方程。計算區域將劃分 成有限個網格點，然後於網格點上對統御方程進行離散與計算。
由式 (3.1) 可知，馬克斯威爾方程式可區分為橫向電場極化 (transverse electric
同樣地，由 (3.1) 之方程組可得橫向磁場極化 (transverse magnetic polarization，TM mode)，其方程由 (Hx, Hy, Ez) 所組成 離存放的作法 (staggered grid)；本文是將電場與磁場存放在同一點 (non-staggered grid)，圖3.2表示非交錯式網格系統。以非交錯網格系統處理空間離散時，必須特 別注意是否會引起不必要之非物理性的誤差 (奇偶次項震盪所造成的問題)。以
∂Hx/∂y 於 (i, j, k) 點離散展開為例，需將參考點 ∂Hx/∂y|(i,j,k) 本身也納入考慮，
以避免奇偶震盪問題 (even-odd decoupling problem)；為了提升數值穩定性，吾人
引入了緊緻格式 (compact scheme) 的方式 [? ? ]。參閱圖3.3，∂H∂xy|i,j,k和 ∂H∂yx|i,j,k在
並透過第二類型之修正方程 (Modified equation of second kind) 以及色散關係保持 (Dispersion relation preserving, DRP) 的觀念，以求得具最佳色散關係保持性的離散 方程式。
3.2 具辛結構之 PRK 時間離散方法
馬克斯威爾方程敘述了電磁波於空間和時間之傳播，它們必須滿足漢彌爾頓 結構 (Hamiltonian structure) 的性質 [27, 28]。於 [29, 30] 等文章中，說明了若要得 到具漢彌爾頓的性質，需引入具辛結構 (symplectic) 之時間格式，方能比未引入 之時間格式得到較佳之理論守恆性質。本論文將使用具辛結構性質之時間離散格
其中，漢彌爾頓系統 H 可寫成如下的分離格式 具辛結構的顯式時間格式 (explicit symplectic time integrator) [31]。本文的時間離 散使用二階段 (two-stage) 二階準確 (second-order accurate) 之具辛結構之顯式分離 Runge-Kutta 方法 (explicit symplectic partitioned Runge-Kutta method)，即
Q0 = qn
(Modified equation of second kind) 式。
3.3 空間離散方程之推導
本章節將在空間離散格式 (3.6-3.7) 以及具辛結構之時間顯示離散格式 (ex-plicit symplectic partitioned Runge-Kutta) 的基礎下，利用第二類修正方程之分析 (Modified equation of second kind analysis)，推導空間之離散方程，
首先，將式 (3.7c)，回推一個時間步，可得 Hn = Hn−1/2− 2µdt∇ × En。將上
Ez|n±1/2i,j,k = Ez|ni,j,k± 1 威爾方程組推導出的電場或磁場波動方程式 (wave equation)，並將其使用泰勒級 數展開，將上述之高階時間微分項轉為空間微分項，所得的代換結果如下
接著，強迫式 (3.3.4) 滿足原方程 ∂E∂tz = 1ϵ (∂Hy
∂x −∂H∂yx)
，即可推導出在三維情況下 的兩條代數方程
3a1+ 2a2+ a3 = 1
2 (3.3.5)
9a1+8 3a2+1
3a3− Cr2x
12 (3a1+ 2a2+ a3) = 0 or
9a1+8 3a2+1
3a3− Cr2y
12 (3a1+ 2a2+ a3) = 0 (3.3.6) 其中 Crx ≡ c∆t∆x，Cry ≡ c∆t∆y 為主導方程式之無因次參數。透過其他方向分量之安 培定律式 ((3.3.1b)-(3.3.1c))，同樣能推導出如上述之代數方程組。此時，仍缺少 一條代數方程式，方得以求解三個未知數。最後一條方程式的推導，係利用波動 方程 (wave equation) 離散式，引入平面波的概念，以推導出數值色散關係方程式 (Dispersion relation equation)，將計算域的時間、空間同步地導入頻域中的角頻率 與波數，以提高格式對波數的保持能力，參閱圖3.4。關於波動方程離散式以及色 散關係分析將於下一章節進行討論。
Fig 3.1: 交錯式之 Yee 網格系統。
( , )
z
x
y
E H
(E,H) (E,H)
(E,H)
(E,H) (E,H) (E,H)
(E,H)
(E,H)
(E,H) (E,H) (E,H)
Fig 3.2: 本文所採用之非交錯式網格系統。
Fig 3.3: 非交錯網格內部節點之示意圖及編號。
t
x
(a) 時間與空間座標系之示意圖。
ω
κ
(b) 角頻率與波數座標系之示意圖。
Fig 3.4: 座標系示意圖。
第 四 章 具色散關係式保持特性之離 散方法及其分析
本文將在本章節中引入于第三章所推導之格式 (3.6-3.7)，針對波動方程式 (Wave equation) 整理出差分離散式。對於具色散關係保持性的離散方程 (Dispersion relation equation, DRE)，引入平面波的概念，將原本時間 (time) 對應空間 (space) 之座標系統轉換為角頻率 (Angular frequency) 對應波數 (Wavenumber) 的座標系 統，使得時間與空間在計算上得以取得更正確的聯結。為了避免數值色散誤差 (dispersion error) 的累積，本文將導入可以精確描述相速度的概念，進而求得具最 佳色散關係保持性的離散方程，以提高空間離散格式與時間格式的合宜連結。此
的時間離散格式∂ϕ∂t|ni = ϕ| (Dispersion relation equation)
1
方位角 (Azimuth angle)，參考圖4.1。式 (4.1.6) 經整理後，可得
為了盡量提高並保持準確的色散關係式，首先定義誤差方程 (Error function) 為[
常需伴隨著混疊誤差 (aliasing error)；為了能描述一完整波形，整個積分範圍將 盡量地保持或貼近 π。圖4.2為固定 Cr = 0.2 時，積分域 m 為 12、3
7、2
5、3
8，以
及 13 與實解之色散關係之比較圖，橫軸及縱軸分別代表了修正的波數 (modified wavenumber) kh 和角頻率 (angular frequency)。
吾人的目標是儘可能地使數值的色散關係貼近實解的色散關係，以提高相位 的正確性。當積分域越趨近完整波形 π 時，越能保持離散方程的色散關係。
4.1.2 Cr 數之影響
消去泰勒級數 (Taylor series) 中的低階誤差項次並不能將時域與頻域作一物理的聯 結。因此，配合顯式具辛結構時間格式 (explicit symplectic time integrator)[31? ] ， 可推導出具色散關係之離散方程 (Dispersion relation equation，DRE)。以此方程再 將數值相速度與實解相速度之間誤差作一優化，可以求得具色散關係式保持的第4.2 數值分析之結果與討論
本文對於空間的一次微分項，在緊緻格式的架構下 [32, 33]，配合顯式具辛結 構時間格式 (explicit symplectic time integrator) [31]，利用降低數值相速度與實解相 速度之間誤差的理念，求得了具色散關係式保持性質的係數 a1 ∼ a3。本文所開發 的具色散關係保持性質之空間離散系統，特別針對在不同天頂角 (Zenith angle) 和 方位角 (azimuth angle) 的情況下，開發各天頂角 (Zenith angle) 和方位角 (azimuth angle) 所對應的係數，以確保波在不同的傳播角度時皆能得到色散關係得以完美 的保持，以降低因各向異性 (anisotropy) 所造成的誤差。
對於因差分離散格式所引入的非物理性數值色散現象，透過前兩章節的分析，
可知本文所發展之具色散關係保持性質的數值方法，對於因數值色散所引起的色 散誤差以及各向異性皆具良好的色散保持效果。針對不同傳播角度所求得相對應 具色散關係保持的係數，能夠有效地減少因各向異性所造成的誤差。
Cr EzL2−error norm CPU TIME(s) 0.2 2.60336× 10−4 46.34 0.05 3.21679× 10−5 192.78
Table 4.1: 經長時間 (即 T 為 30(s)) 計算後，在 Cr = 0.2 及 Cr = 0.05 情況下計算 誤差以及所需 CPU TIME (s) 之比較。
Table 4.2: 取天頂角 (Zenith angle) θ 為 0◦ 和 90◦、30◦ 和 60◦ 以及 45◦，配合方位角 (Azimuth angle)ϕ 為 0◦和 90◦、6◦和 84◦、9◦ 和 81◦、12◦ 和 78◦、22.5◦ 以及 67.5◦、 30◦ 和 60◦、36◦和 54◦、6◦ 和 84◦ 以及 45◦ 所求係數 a1 ∼ a3 的分布情況。當天頂 角 (Zenith angle)θ 為 90◦ 時，可得到二維情況下的係數；當天頂角 (Zenith angle)θ 為 90◦和方位角 (Azimuth angle)ϕ 為 0◦ 或 90◦ 時，可得到一維情況下的係數。
present Yee’s method 1D ∆t≤ 1.167132hc ∆t≤ hc 2D ∆t≤ 0.825287hc ∆t≤ 0.707107hc 3D ∆t≤ 0.673844hc ∆t≤ 0.57735hc
Table 4.3: 在三維情況下，本文所提出的 FDTD 數值方法與 Yee 方法的穩定性範圍 比較。
Fig 4.1: 平面波在三維情況下傳遞時，定義出天頂角 (Zenith angle) θ，方位角 (Azimuth angle) ϕ。
κx
ω
0 0.5 1 1.5 2 2.5 3
0 0.5 1 1.5 2 2.5 3
Exact m = 1/3 m = 3/8 m = 2/5 m = 3/7 m = 1/2
Fig 4.2: 固定 Cr 為 0.2 時，在不同積分範圍下之色散關係之實解與數值解的比較。
× × × × × × × × × × × × × × × × × × × × kh
ω
0.5 1 1.5 2 2.5 3
0.5 1 1.5 2 2.5 3
Exact Cr = 0.05 Cr = 0.2 Cr = 0.3 Cr = 0.4 Cr = 0.5 Cr = 0.6
×
Fig 4.3: 固定積分範圍 −π2 ∼ π2 時，在不同 Cr 值情況下之色散關係之實解與數值 解的比較。
第 五 章 數值方法之驗證
在進行電磁波問題的模擬前，為了確保分析結果之正確性，必須先對數值程 式作一驗證。本章節將以具解析解之三維馬克斯威爾方程來驗證本文所提出之顯 式 (explicit) 具辛結構 (symplectic) 的時間離散與具 DRP 性質之空間離散的格式；
本論文所使用之數值方法為直接求解法拉第定律與安培定律以得到電場與磁場，
電場、磁場與時間收斂階數之比較結果如表5.1所示。並針對吾人所開發之方法，
對於空間離散之實際收斂階數與理論上吾人所推導出之空間收斂階數進行比較，
吾人所選取之 h、電場、磁場與空間收斂階數之結果如表5.2所示。
Scheme dt E-maxError H-maxError E-T.r.o.c H.T.r.o.c 0.5h 6.5922E-06 7.8720E-05 - -PRK-FDTD 0.25h 1.6816E-06 2.0362E-05 1.9709 1.9509
0.125h 5.1574E-07 5.7744E-06 1.7051 1.8181 0.5h 1.673E-04 1.511E-04 - -ADI-FDTD 0.25h 1.135E-04 1.015E-04 0.5595 0.5736
0.125h 1.001E-04 8.939E-05 0.1818 0.1839
Table 5.1: 比較本文所使用的 PRK-DRP FDTD 與 ADI-FDTD [1] 在 h = 0.01 情況下 計算電場 E 與磁場 H 的最大誤差和時間收斂斜率。
Scheme h E-maxError H-maxError E-S.r.o.c H.S.r.o.c 1/10 3.0899E-05 1.0312E-04 - - PRK-FDTD 1/20 2.5787E-06 1.1957E-05 3.5828 3.1084
1/30 5.2522E-07 3.4710E-06 3.9244 3.0505 1/40 2.3758E-07 2.1476E-06 2.7576 1.6688
Table 5.2: 本文所使用的 PRK-DRP FDTD 當 t = 1，並選取 CFL number =0.2 時，
計算電場 E 與磁場 H 的最大誤差和空間收斂斜率。
5.2 驗證結果與討論
藉由實解驗證的結果，可以確保本文所開發之具辛結構 (symplectic) 的時間離 散與具 DRP 性質之空間離散的數值方法，在求解馬克斯威爾方程問題時，能量密 度經長時間後仍得以保持漢彌爾頓之守恆性質。此外，利用磁場的散度隨時間的 變化，可以證明此數值方法在僅求解法拉第定律和安培定律的情況下，仍可滿足 馬克斯威爾方程組中的高斯定律。與傳統之 Yee 離散方法做一比較後可知，此方 法具有更準確且更高之空間收斂階數。
透過實解驗證，可以確保本論文開發之具辛結構的時間離散格式，以及具備 最佳數值色散關係之空間離散格式的數值方法，在非交錯網格下求解馬克斯威 爾方程組問題時，具有良好的效果。由結果可以得知，吾人所開發之 PRK-DRP FDTD，在時間上具有二階準確，在空間上更具有四階準確。此外，在和 ADI-FDTD [1] 結果比較後，可知 PRK-DRP ADI-FDTD 具備更佳的收斂性以及準確度。在 接下來的章節，將使用本文所開發之 PRK-DRP FDTD 進行實際問題的模擬應用。
第 六 章 人體電磁比吸收率之分析
6.1 實際物理問題之描述
隨著科技日新月異，目前手機在我們的日常生活中已被廣泛的使用，然而手 機因發射與接收訊號時的輻射可能會引起健康問題的疑慮，是一項相當值得重視 的議題。由於使用手機時，人體將暴露在近場電磁場之輻射中，其中是否將引起
隨著科技日新月異，目前手機在我們的日常生活中已被廣泛的使用，然而手 機因發射與接收訊號時的輻射可能會引起健康問題的疑慮，是一項相當值得重視 的議題。由於使用手機時，人體將暴露在近場電磁場之輻射中，其中是否將引起