Chapter 3 數值方法
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 實際物理問題之描述
隨著科技日新月異，目前手機在我們的日常生活中已被廣泛的使用，然而手 機因發射與接收訊號時的輻射可能會引起健康問題的疑慮，是一項相當值得重視 的議題。由於使用手機時，人體將暴露在近場電磁場之輻射中，其中是否將引起 與導致癌症風險的問題，也引起了學界與業界大量的討論。本研究題目之目的，
是模擬手機射頻所導致的電磁場 (electromagnetic field, EMF)。輻射的程度受到許 多因素影響，包括手機的種類，參考的發射頻率，以及手機所使用的功率等等。
為了使我們的模擬能呈現出真實的情況，針對其中一款目前市售的手機進行其複 雜幾何之模擬 [34]。為此，本研究將參考 Apple iPhone4-like phone 之結構圖以及 其所發射之 EMF 進行模擬與分析。此外，使用者距離手機天線的距離、手機的使 用時間，以及手機與人體之間的介質，決定了電磁對人體之曝曬程度的另外兩個 重要參考因素，在以下的小節，我們將介紹所模擬的手機，並且詳細的描述涉及 的模擬問題之各項係數，以進行接近真實情況的模擬預測。
圖6.1描述了吾人所欲探討之問題，在一般的情況下，使用手機通話的情形，
進行趨近真實情況的模擬。吾人假設手機緊貼於頭部耳朵一側，且手機相對於人 體呈 45 度角，手機波源的選取，取工作頻段 900Mhz、600mW 之功率的正弦波，
作用於 UMTS/GSM 手機天線輻射部上進行天線輻射，以觀察手機全場輻射對於 頭部之分布。
Fig 6.1: 吾人所討論之實際問題描述。
6.2 三維複雜幾何的散射體建模
首先，吾人利用 CAD 設計圖，可以獲取複雜物體各部位的幾何結構。接著，
將複雜物體進行拆解，分別對每個部件進行建模與設置參數。並將個別部件之 CAD 設計圖轉換為.STL (STereoLithography, 立體光刻) 的 3D 列印格式檔，可以得 到以三角形網格作為結構的 3D 物體的節點資訊。最後，將這些節點資訊引入計 算空間中，可以建立在空間網格中所需的模型與其材料參數。
考慮計算空間中的一點 P，欲判斷 P 點是否位於三角形之內，可以透過 xOy 面 (或 yOz 和 xOz) 的投影來分析。首先假設三角面 ABC 和空間中一點 P 在 xOy 面內的投影分別為 A’、B’、C’和 P’。點 P’和 A’、B’、C’可分別構成向量 P’
A’、P’B’和 P’C’，最後依次做向量外積，分別為 P’A’×P’B’、P’B’×P’
C’和 P’C’×P’A’。如果三個向量外積後為同號，則 P 點位於三角形內，反之，
P 點位於三角形外。倘若三個向量外積後其中一者為零，則該點位於三角形邊界 上。
Fig 6.2: 計算空間中節點 P 與三角面三點 ABC 之關係圖。
Fig 6.3: 真實手機的各個部件。
一般在計算網格空間中判斷網格線是否與三角面具有交點時，有以下的四種 可能，其中包括
A：只有一個交點，表示網格線與目標之表面相切；
B：共有兩個交點，表示直線貫穿目標表面；
C：共有三個交點，表示直線與目標只有一個相切點
D：共有四個或以上之交點，表示目標具有多個”凸起”(鋸齒狀幾何，例如：鋸 子或梳子上之紋路)。
將不同材質的部件分別進行材質的內外點判斷，最後，將部件內部點及部件邊界 紀錄為 1，部件的外部點紀錄為 0，設定部件內部點以及部件邊界上的材料性質 (包括相對介電常數、相對磁導係數、導電率、密度等所需材料參數) 最後輸出作
將不同材質的部件分別進行材質的內外點判斷，最後，將部件內部點及部件邊界 紀錄為 1，部件的外部點紀錄為 0，設定部件內部點以及部件邊界上的材料性質 (包括相對介電常數、相對磁導係數、導電率、密度等所需材料參數) 最後輸出作