Chapter 1 序論
1.5 論文大綱
本文本諸理想的馬克斯威爾方程與其色散關係方程式的核心概念，發展出一 具最佳數值色散關係式之非交錯格式的空間離散方法。首先，於第一章中介紹計 算電磁波的相關背景與發展之文獻回顧，以及本研究的動機與目標。本文的第二 章介紹三維的馬克斯威爾方程組，以及在時域中，開發時間與空間具物理相關特 徵的顯式有限差分離散格式。第三章中，為了有效的執行有限差分法，此章節簡 要的描述了非交錯網格中所使用的數值模型。然而，所提出的 FDTD (時域有限差 分法) 數值模型滿足了具辛結構以及數值色散關係方程的性質。此外，也針對被 卷積完全匹配層 (CPML) [9] 所包圍的三維計算空間中進行了簡要的描述。第四章 中，詳細敘述了吾人所開發之數值方法，將時間與空間之離散的方法進行推導。
在第五章，透過實解驗證，並且與 ADI-FDTD [1] 交替方向隱式的時域有限差分方 法進行時間及空間上收斂斜率之比較，同時驗證了本文所開發之方法準確性，以 及程式的可信度。第六章中，詳細介紹了正在研究的問題，其中包括貼置在模擬 人頭部耳上的 Apple iPhone4-like phone 的模型建模，以及建模方法。並且回顧了 大腦組織的主要組成。在第七章中，詳細討論了透過求解馬克斯威爾方程組模擬 手機對人腦之輻射問題的模擬結果。最後，我們將在第八章中給出結論以及未來 的展望。
第 二 章 馬克斯威爾方程組
2.1 法拉第/安培/高斯方程組及其推導
馬克斯威爾方程組係由法拉第定律、安培定律及高斯定律所組成。在時域中，
分別可以表示為如下的安培定律及法拉第定律表示式，
∇ × ⃗E = −∂ ⃗B
∂t − ⃗Jm, (2.1)
∇ × ⃗H = ∂ ⃗D
∂t + ⃗J , (2.2)
以及高斯定律所組成的限制方程式
∇ · ⃗B = 0, (2.3)
∇ · ⃗D = ρ. (2.4)
上述方程組中， ⃗H 為磁場強度 (Magnetic field intensity)、 ⃗E 為電場強度 (Electrical field intensity)、 ⃗B、磁位移 (Magnetic displacement) 和 ⃗D 電位移 (Electrical displace-ment)。式 (2.2) 與 (2.1) 中之 ⃗J 和 ⃗Jm分別代表電流密度 (Electric current density) 和 磁流密度 (Equivalent magnetic current density)。于 (2.4) 之 ρ 代表電荷密度 (Electri-cal charge density)。在各向同性 (Isotropic) 的線性介質中需滿足如下的本構關係 式
D = ϵ ⃗⃗ E, ⃗B = µ ⃗H, ⃗J = σ ⃗E, ⃗Jm = σmH.⃗ (2.5)
上 式 中，ϵ 和 µ 分 別 表 示 介 電 係 數 (Electric permittivity) 和 磁 導 係 數 (Magnetic permeability)，σ 和 σm 分別為電導率 (Electric conductivity) 和磁導率 (Equivalent magnetic loss)，兩者分別代表介質的損耗和磁損耗。本文將所推導的馬克斯威爾 方程之離散方法係居于顯式的格式，用來模擬電磁波於簡單介質 (Simple medium) 中的傳播行為。主要探討之材料性質係數中的導磁係數 (Magnetic permeability，µ) 以及介電係數 (Electric conductivity，ϵ) 皆視為常數，此二材料性質係數主要是在
描述材料之介電性質以及導磁特性且皆為正其中，ϵ = ϵ0· ϵr 而 µ = µ0 · µr；兩材 到瞬時之電場與磁場。對於式 (2.8-2.9) 之無散度條件 (divergence-free constraint conditions) 之限制，吾人將於程式驗證中進行分析。
馬克斯威爾方程組可以利用漢彌爾頓動力系統 (Hamiltonian dynamical system) 表示之，即
其中，漢彌爾頓方程 (Hamiltonian function，H) 可以下式表示之 [21]
H( ⃗H, ⃗E) =
對稱矩陣 (skew-symmetric matrix) 表式之 [22]
G =
( 0 −µ−1/2∇ × (ϵ−1/2) ϵ−1/2∇ × (µ−1/2) 0
)
. (2.12)
馬克斯威爾方程組可由 ∂ ⃗ψ(t) / ∂t = G ⃗ψ(t) 表示之，如此吾人可得該方程之通解 ψ(t) = exp⃗
(
t G ⃗ψ (t = 0)
)，其中 exp(t G) 顯示了馬克斯威爾方程組之解係隨時間 演進呈指數的變化。將向量 ⃗ψ 特徵化為∫
Ωϵ ⃗E· ⃗E + µ ⃗H· ⃗HdΩ，它與電磁場之能 量密度 (energy density，w(t)) 存在直接的理論關聯，即
w(t) =
∫
Ω
ϵ ⃗E· ⃗E + µ ⃗H· ⃗H dΩ. (2.13) 式 (2.13) 之能量密度不會隨時間而改變。對於式 (2.11) 及 (2.13) 隨時間皆不變之性 質，吾人可將其作為程式驗證之用。
2.3 色散介質
常 見 的 色 散 介 質 模 型 包 括 德 拜 模 型 (Debye model)、 洛 侖 茲 模 型 (Lorentz model) 和德魯模型 (Drude)，它們所對應的單極介電係數表達式如下:
(1) Debye 介質 其中 ⃗Jd項為介質的色散極化所引起的極化電流項 (Polarization current density)，對 於不同類型的色散介質，它們分別滿足如下的不同輔助偏微分方程
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) 所造成的誤差。
對於因差分離散格式所引入的非物理性數值色散現象，透過前兩章節的分析，
對於因差分離散格式所引入的非物理性數值色散現象，透過前兩章節的分析，