Chapter 4 具色散關係式保持特性之離散方法及其分析
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，設定部件內部點以及部件邊界上的材料性質 (包括相對介電常數、相對磁導係數、導電率、密度等所需材料參數) 最後輸出作 為材料檔。
計算空間內，吾人所建立之模擬 iPhone4-like phone 模型如圖6.4所示。
Fig 6.4: 計算空間中，手機各個部件的建模。
6.3 波源的設置
時域有限差分法在模擬電磁波傳遞時，對於觀察場，若要求得下一個時間步 的電場或磁場值，必定是以上一個時間步所得到的答案計算出的，因此，如果沒 有源項，透過時域有限差分法計算出的任何時刻位置的電場或磁場的值必定等於 零‧因此，有效且正確的選取合適的波源，是能否完成一個正確的計算程序的重 要因素之一，以下將以常見的波源設置方法簡要介紹之。
6.3.1 硬波源 (hard-sourced)
硬波源屬於最簡易的源項處理方式，只需要在 FDTD 的計算流程中，在波源 位置處，插入欲使用的離散波函數當作源項即可，它可以是以電場或磁場作為源 項。以正弦波函數為例‧源項可以寫成：
Ez = E0sin(ω0t). (6.1) 將時間離散化後可得
Ez|nis = E0sin(2πf0n· ∆t). (6.2) 其中 is 為波源的源點，ω0、f0 為正弦波的角頻率與頻率，n 為時間迭代的步數，
∆t 為時間的間隔。硬波源雖然在程式撰寫上相當的簡單，由於波源是在固定的空 間中隨時間變化，如果遇到散射體將會導致二次反射，使得求解結果不正確。解 決的辦法為反射波到達波源前關閉源項並作更新，以確保不會造成二次的反射。
6.3.2 軟波源 (soft-sourced)
軟波源又稱附加源，在安培定律或法拉第定律中加入電流源 (Jsource) 或磁流
源 (Msource) 做為電磁或磁場的能量。當反射波傳遞到此波源區域時，會直接與電
流源或磁流源做一疊加，而不致于造成計算時的反射。
∂E
∂t = 1
ε[∇ × H − Jsource] , (6.3)
∂H
∂t =−1
µ[∇ × E − Msource] . (6.4)
6.4 波源的種類
使用時域有限差分法分析電磁問題時，針對不同問題，選擇合適的入射波作 為波源種類‧將入射波加入到 FDTD 之時間疊代過程中，可以更方便的針對不同 的電磁問題進行分析。波源主要包括兩種，一種是隨時間週期變化的時諧場源，
另一種是對時間行進過程，呈現鐘型高斯函數的波源。以下進行兩種波源的簡要 介紹
6.4.1 時諧場源 (harmonic-sourced)
以單調正弦波 (式6.2) 為例，使用此時諧場波源，一般若在計算域中欲使系統 達到穩態，通常需要經過 3 至 5 個週期，其達到穩態的所需時間，也與散射物體 所需的週期，以及散射物體的形狀和大小有關 [35] 。
6.4.2 脈衝源 (implused-sourced)
脈衝源的頻譜範圍通常具有一定的頻寬，通過了解脈衝源，以及其頻譜特性，
對於 FDTD 的計算上十分重要，以高斯脈衝函數為例，高斯脈衝在時域下的型式 為:
Ei(t) = exp(−4π(t− t0)2
τ2 ). (6.5)
其中 τ 為常數，它決定了高斯脈衝函數的頻寬。脈衝的最大值出現在 t = t0 的時 刻，時域下之高斯函數隨時間變化圖如圖6.5 所示。由於高斯函數經過傅立葉變換 後，在頻域下，其波型仍然為一個高斯脈衝，因此可利用此一性質，方便的求出 一個波源，在所關注的頻寬下的頻譜圖，進行進一步的分析。