頻域有限差分法應用於二維光波導結構之探討

全文

(1)♁. 國立中山大學光電工程研究所碩士論文. 頻域有限差分法應用於二維光波導結構之探討 Analysis of 2D Optical Waveguide Structures Using Frequency-Domain Finite-Different Method. 研究生：林政衛撰指導教授：張弘文博士. 中華民國九十三年六月.

(2)

(3)

(4)

(5) 誌. 謝. 本篇論文的完成，首先要感謝我的指導教授張弘文博士，在他悉心的教導之下，使我可以在數值理論的基礎上紮根，並且對於論文理論以及程式方面給予多方的協助，並使我在對學習研究以及為人處事上獲益良多。另外還要感謝所上的老師在學業上所提供的協助，由於他們的教誨，使我獲得不少知識以及幫助。同時要感謝在實驗室中的二位博士班學長盛夢徽學長、吳祚倫學長，兩年來照顧我的學業以及生活，同時在作論文過程時更是提供了不少寶貴的建議。也感謝學弟王勝民、劉漢強、蔡昂勳、李致維、洪慶龍在口試期間的幫忙，讓我能安心準備口試。最後要感謝我父母辛苦的養育和栽培，有了他們在背後的大力支持，我才能夠完成碩士學業以及論文研究。本論文感謝中科院一所航材組資助(NSC91-2623-7-110-002)。. 林政衛. 2004,06 于高雄西子灣. I.

(6) 摘要. 本論文為運用頻域有限差分法(finite different frequency domain) 求解光入射於光被動元件，即光波導，所產生的場型。相較於時域有限差分法(finite different time domain)，頻域有限差分法需解聯立方程式，並且計算其矩陣較大，也沒有理想的近似快速解法。但頻域有限差分法應用於光波導的模擬有許多的優點，使其相較於時域有限差分法，更為有效率。而光波導為在波長分工，已知通常亦為特定頻率或波長，而頻域有限差分法有模態解析解，可合併處理以減少計算量。而我們所需探討，為所設之格點其有效折射率(index)如何取樣，以及在偏微分方程中必定遇到的邊界條件的處理，主要包括電牆 (electric Wall)、磁牆(magnetic Wall)、週期牆(periodic wall)、吸波邊界 (absortion boundary)等邊界條件。其次分別探討 TE、TM 波處理之不同。最後我們探討高次的有限差分，以及計算極大區域或極高精密度的方法；以期能利用頻域有限差分法，以最快、最準的程序以找出各種光波導的入射、反射、散射的場型。. II.

(7) Abstract. This thesis applies the method of finite-difference frequency -domain to solve the field when light propagating through the optical waveguide. Comparing with other method, method of finite-difference time-domain which required to be calculated with larger matrices and without any ideal approximate rapid solution, our method only needs to solve the joint equations. Method of finite-difference frequency -domain will be more efficient than method of finite- difference time-domain in optical waveguide simulation. However, all that we need to confer with are how to sample the effective indices for set grind points and the deal for the boundary conditions in partial differential equations, including electric wall, magnetic wall, absorption boundary and other boundary conditions, etc. Next, we confer with the difference deal between TE and TM waves respectively. Final, we confer with high-order finite-difference and the method for calculating extremely large area or extremely high precision in order to use the method of finite-difference frequency-domain to find the incident, reflective, and scattering fields by the fastest and the exactest procedures.. III.

(8) 目錄頁次. 誌謝. I. 摘要. II. 英文摘要. III. 內容目錄. IV. 圖表目錄. VI. 第一章、導論. 1. 第二章、波動方程式推導. 3. 2.1 馬克思威爾方程式與向量基本公式. 3. 2.2 一維純量波動方程式推導. 5. 2.3 二維純量波動方程式推導. 9. 第三章、以有限差分法近似波動方程式. 11. 3.1 有限差分法理論概述. 11. 3.2 基本有限差分式的推導. 13. 3.3 二維有限差分推導. 15. 3.4 四點有限差分. 18. 第四章、一維有限差分計算推導 4.1 邊界條件處理和推有限差分式. IV. 23 25.

(9) 4.2 三層介質一維有限差分法計算. 第五章、運用頻域有限差分求解波導場型. 29. 33. 5.1 計算區域有限差分式的近似. 34. 5.2 二維頻域有限分模型探討. 38. 5.3 頻域有限差分大矩陣計算簡化法. 51. 5.4 空心波導數值模擬結果. 56. 5.5 介質波導數值模擬結果. 59. 5.6 波核漸細波導數值模擬結果. 63. 第六章、結論. 68. 參考文獻. 70. 中英對照表. 72. V.

(10) 圖目錄圖 2-1 一維多層介質. 5. 圖 2-2 二維模型. 9. 圖 3-1 微分示意圖. 11. 圖 3-2 以割線近似切線. 12. 圖 3-3 h 趨近於零的情況. 12. 圖 3-4 一維有限差分關聯圖. 14. 圖 3-5 二維有限差分示意圖. 15. 圖 3-6 二維有限差分設點方式. 16. 圖 3-7 四點有限差分示意圖. 18. 圖 3-8 二維四點有限差分取點示意圖. 21. 圖 4-1 單一界面設點. 24. 圖 4-2 介面處理示意圖. 26. 圖 4-3 TE wave, Ey. 30. 圖 4-4 TE wave, Hx. 30. 圖 4-5 TM wave, Hy. 30. 圖 4-6 TM wave, Ex. 31. 圖 5-1 任意光波導結構示意圖. 34. 圖 5-2 二維有限差分關聯圖. 35. VI.

(11) 圖 5-3 ε、 r ε、 l ε、 u ε d 的等效示意圖. 37. 圖 5-4 頻域有限差分二維模型示意圖. 38. 圖 5-5 頻域有限差分取點示意圖. 39. 圖 5-6 頻域有限差分取點示意圖二. 40. 圖 5-7 空心波導上下邊界電牆場形圖. 55. 圖 5-8 空心波導上下邊界磁牆場形圖. 57. 圖 5-9、點波源空心波導上下邊界電牆場形圖. 58. 圖 5-10、點波源空心波導上下邊界磁牆場形圖. 58. 圖 5-11、點波源空心波導上下邊界理想邊界場形圖. 59. 圖 5-12 介質波導示意圖及其磁牆等效圖. 60. 圖 5-13 介質波導正解之找根圖. 62. 圖 5-14 介質波導正解之場型圖. 62. 圖 5-15 介質波導頻域有限差分解之場型圖. 63. 圖 5-16 Tapered Waveguide 示意圖. 63. 圖 5-17 Tapered Waveguide 區域劃分示意圖. 65. 圖 5-18 Tapered Waveguide 基模態場型圖. 65. 圖 5-19 Tapered Waveguide 第二高階模態場型圖. 66. 圖 5-20 Tapered Waveguide 第四高階模態場型圖. 66. 圖 5-21 吸收區示意圖. 67. VII.

(12) 第一章導論隨著資訊的蓬勃發展以及網際路應用的急速增長，人們對傳輸和訊號的通道頻寬需求愈來愈高。由於光波導物理性質穩定、不受電磁干擾、頻寬大、損耗低等絕對優點，將成為網路傳輸的主要媒介。然而在設計時要事先了解光於光波導的場型分佈並不容易，這是因為光於光波導中的場型並不能單純以幾何光學(Ray Optics)、波動光學(Wave Optics)能夠求解，必須以電磁光學(Electromagnetic Optics) 求解。然而以電磁光學求解時，必然需要處理各種偏微分方程(Partial Different Equation)；一般來說，一旦邊界條件(Boundary Condition)較為複雜，偏微分方程皆難以求解，於是各種數值方法因運而生，如有限元素法、時域有限差分法、雙向多模解析法、解析連續法等等。其中時域有限差分是一個成熟的數值方法，並已應用於工程計算的各個領域。而相較於時域有限差分法，頻域有限差分法是一個新穎的數值方法；時域有限差分法其優點是計算時間正比於 N x × N z ，並且計算一次就可得到全部頻率；缺點則計算區域較大(因額外有吸波. ) µ (ω ) 必須皆為窄波區；時間階數(Time Step)約104 ~ 106 ，並且 ε (ω 、頻率常數。而頻域有限差分為較為新穎的方法，其優點是精密度較. ) µ (ω ) 可隨頻率而高，可進行色散分析，無時域有限差分近似， ε (ω 、變。缺點是一次只能計算一個頻率，計算時間正比於 N X 3 × N Z ，需求 1.

(13) 解個數為 N x × N z 的聯立方程式。本篇論文為延續張永豐學長的論文“寬頻電磁波吸波材料設計與數值模擬”，和學長的論文比較，我專門探討頻域有限差分法於光波導元件的應用，並且我的方法可計算更大的區域，以及更高的精密度，並且在解析解區可使用多層結構。因一般光波導皆注意其頻域特性，所以頻域有限差分法特別適用於光波導。首先於第二章“理論與公式之推導”找出適用於光波導並簡化的波動方程式。並於第三章“以有限差分法近似波動方程式”，簡介有限差分法，並得到波動方程式的有限差分型式。第四章“一維有限差分計算推導“討論在一維情況，有限差分的型式，以推導更複雜二維型式。第五章“二維頻域有限差分模型推導“，面臨的問題和一維類似，作法是以一維的作法延伸到二維，並使用各種計算方法以減少計算量。. 2.

(14) 第二章波動方程式推導 2-1 馬克思威爾方程式與向量基本公式. 於本篇論文，如何描述光場是一項最基本的問題，而本篇論文所用的方法是將光場以電磁波的波動方程式來加以描述，而以電磁波作描述就必然使用馬克思威爾方程式(Maxwell’s equations)及一部分基本的向量公式，所以在此為我們先對馬克思威爾方程式作一基本的描述。當存在源點(source point)，需要考慮產生電磁波的區域，馬克思威爾方程式可寫成如下形式：. ∇× E = − ∇× H =. ∂ B ∂t. (2.1-1). ∂ D+J ∂t. (2.1-2). ∇i D = ρ. (2.1-3). ∇i B = 0. (2.1-4). 其中 E 為電場強度 (electric field intensity) ， H 為磁場強度 (magnetic. field density) ， D 為電位移 (electric displacement) ， B 為磁通密度 (magnetic flux)， J 是電流密度(electric current density)， ρ 則為體電荷密度(volume charge density)。當傳播介質具有等方向性(isotropic)、均勻(homogeneous)、及非色散(nondispersive)此三項性質時，電場強度. 3.

(15) 和電位移就會有 D = ε E 的關係；同理，磁場強度和磁通密度則是有. B = µ H 的關係；其中 ε 稱為介電係數(permittivity)， µ 稱為導磁係數 (permeability)，其中 ε 、 µ 此時皆為與座標無關的常數；假若我們所處理的是一般的光波導或光纖元件，則無涉及磁性材料，所以介質的導磁係數 µ 等於在自由空間的磁導係數 µ0 。且若為無源區域(source free region)，則方程式(2.1-2)和(2.1-3)分別將會簡化成下面的形式：. ∇× H =. ∂ εE ∂t. (2.1-5). ∇iε E = 0. (2.1-6). 其中，電場和磁場均為時間諧和函數以 e jωt 來表示， µ 、 ε 亦非時間的函數，所以可得 H=. E=. ∂ = jω ，(2.1-1)和(2.1-5)將可頻域的式子如下： ∂t. 1 ∇× E − jωµ 1 jωε. (2.1-7). ∇× H. (2.1-8). 並藉由(2.1-7)、(2.1-8)兩式的互相耦合，即可求出所需波動方程式。. 4.

(16) 2.2 一維純量波動方程式推導. θr. n1. X. n2. nN −1. θi. nN. θt Z. 圖 2-1 一維多層介質使用的例子如圖 2-1，對於一維多層介質的斜角入射，雖然是行經路線為不是一維，但是 x 方向是連續量(phase matching condition)，只有 z 有變化，所以使用一維波動方程式就可表示，波動方程式分別為. TE wave 和 TM wave；由於 TE wave 電場以 y 分量為連續量，所以使用電場連續的波動方程式，TM wave 磁場 y 分量為連續量。其中由於光波導一般為非磁性材料，所以我們導磁係數(permeability)皆可設為 µ0 。. 2.2.1、TM wave 波動方程式令入射平面為 zx 平面，而 TM wave 即水平極化，磁場強度. H = H y aˆ y ，其中 H y = He − jk ( x sin θ + z cosθ ) ，H 為振幅。將式(2.1-8)等式左右兩邊取旋度 ∇× E = ∇×(. 1 jωε ( z ). ∇ × H ) = − jωµ0 H 5. (2.2-1).

(17) ∇×(. 1 ∇ × H ) = ω 2 µ0 H ε ( z). (2.2-2). 其中 ε ( z ) = ε r ε 0 ，代入式(2.2-2)得. ∇×(. ∇×(. 1. εr 1. εr. ∇ × H ) = ω 2ε 0 µ0 H. (2.2-3). ∇ × H ) = k02 H. (2.2-4). 磁場強度 y 方向為連續，微分為零，即. ax ∂ ∇× H = ∂x 0. ay ∂ ∂y Hy. ∂ =0 ∂y. az ∂H y ∂H y ∂ ax + az =− ∂z ∂z ∂x 0. (2.2-5). 由式(2.2-5)代入式(2.2-4)得. ∂H y ⎡ 1 ∂H y ⎤ ∇ × ⎢ (− ax + az ) ⎥ = k02 H ∂z ∂x ⎣εr ⎦. ∂H y ⎡ 1 ∂H y ⎤ ∇ × ⎢ (− ax + az ) ⎥ = ∂z ∂x ⎣εr ⎦ −. ax ∂ ∂x 1 ∂H y. εr ∂z. ay ∂ ∂y 0. (2.2-6). az ∂ ∂z 1 ∂H y. εr ∂x. 2 2 2 ∂ 1 ∂H y 1 ∂ Hy 1 ∂ Hy 1 ∂ Hy + ( = ax ( ) − ay ( )) + a z (− ) ε r ∂ y∂ x ε r ∂x2 ε r ∂ y∂ z ∂z εr ∂z 2 ∂ 1 ∂H y 1 ∂ Hy + = −a y ( ( )) ε r ∂x2 ∂z εr ∂z. (2.2-7). 6.

(18) 將式(2.2-7)代入式(2.2-6)可得一維 TM wave 波動方程式 2 1 ∂ H y ∂ 1 ∂H y + ( (−a y ( ))) = k02 H 2 ∂z εr ∂z εr ∂x. (2.2-8). TM wave 波動方程式只有 y 分向，如式(2.2-9) 2 1 ∂ H y ∂ ⎡ 1 ∂H y ⎤ 2 + ⎢ ⎥ + k0 H y = 0 2 ∂z ⎣εr ∂z ⎦ εr ∂x. (2.2-9). 由於 x 方向為連續量，對微分為常數，將波動方程式簡化，對 x 方向二階微分後，令 p = k sin θ 得式(2.2-10). 1. εr. (− p 2 ) H y +. ⎤ ∂ ⎡1 ∂ H y ⎥ + k02 H y = 0 ⎢ ∂ z ⎣εr ∂z ⎦. 式(2.2-10) 一維波動方程式，一階微分. 二階微分. (2.2-10). 1 ∂ H ∝ E x 為連續項， ε r ∂z y. ⎤ ∂ ⎡1 ∂ H y ⎥ ∝ H y 為連續項，可知波動方程式為連續方程 ⎢ ∂z ⎣ ε r ∂z ⎦. 式。使用此種形式的波動方程式是為了減少不連續項的產生，不連續項只由第一項和第二項貢獻，第三項為純量，以減少數值上計算的誤差。. 2.2.2、TE wave 波動方程式用推導一維 TM wave 波動方程式相同的方式可得一維 TE wave 波動方程式，TE wave 為垂直極化。. 7.

(19) 2 1 ∂ Ey 1 ∂2 + E y + k02 E y = 0 2 2 εr ∂x εr ∂z. (2.2-11). 1 ∂2 E y + k02 E y = 0 (− p ) E y + 2 εr εr ∂z. (2.2-12). 簡化後. 1. 2. 8.

(20) 2.3 二維純量波動方程式推導. 二維波動方程式使用的例子如圖 2-2，對於 zx 方向的介電常數 ε 作變化， ε = n 2 ，即在 xz 平面 n 值都可變化，使用二維波動方程式。. n2. n4. n2. n1. n3. n1. n2. n4. n2. x z. 圖 2-2 二維模型. TM wave 波動方程式 ε 隨 x 和 z 方向變化，推導方式和一維波動方程式推導類似，由式 (2.2-6)其中的旋度運算子. ∂H y ⎡ 1 ∂H y ⎤ ∇ × ⎢ (− ax + a z )⎥ = ∂z ∂x ⎣ε r ⎦ −. ax ∂ ∂x 1 ∂H y. εr ∂z. ay ∂ ∂y 0. az ∂ ∂z 1 ∂H y. εr ∂x. 2 2 ∂ 1 ∂H y ∂ 1 ∂H y 1 ∂ Hy 1 ∂ Hy = ax ( ) − ay ( ( )+ ( )) + a z (− ) ∂x εr ∂x ε r ∂ y∂ x ∂z εr ∂z ε r ∂ y∂ z. = −a y (. ∂ 1 ∂H y ∂ 1 ∂H y ( )+ ( )) ∂x εr ∂x ∂z εr ∂z. 將式(2.3-1)代回式(2.2-6)，可得二維 TM wave 波動方程式 9. (2.3-1).

(21) ⎤ ∂ ⎡1 ∂ ⎤ ∂ ⎡1 ∂ Hy⎥ + ⎢ H y ⎥ + k02 H y = 0 ⎢ ∂x ⎣εr ∂x ⎦ ∂z ⎣εr ∂z ⎦. (2.3-2). TE wave 波動方程式以相同的方法可得二維 TE wave 波動方程式 1 ∂2 1 ∂2 Ey + E y + k02 E y = 0 2 2 εr ∂x εr ∂z. 10. (2.3-3).

(22) 第三章以有限差分法近似波動方程式. 3.1 有限差分法理論概述. 有限差分法即是用離散、只含有限個未知數的差分方程式，去代替連續量的微分方程定解條件；以求出差分方程式的精確解作為求解偏微分方程式的近似解。如圖 3-1：. y h值趨近於0 f′=. dy dx. f(x). x x-h=x=x+h 圖 3-1 微分示意圖. 一階微分公式並取全極限為：. f ′( x) = lim h →0. f ( x + h) − f ( x − h) 2h. (3.1-1). 如圖 3-2，有限差分即是取有限小的 h 以近似微分，即以割線去近似切線。 mt =. f ( x + h) − f ( x − h) f ( x + h) − f ( x − h) = ( x + h) − ( x − h) 2h 11. (3.1-2).

(23) y. mt. f(x+h) f(x) f(x-h). x x-h. x. x+h. 圖 3-2、以割線近似切線. 然而近似會有誤差的，而誤差的大小正比於 ∆x 的大小。如圖 3-3，h 趨近於零，差分即等於微分。. y h值趨近於0. f(x). x x-h=x=x+h 圖 3-3、h 趨近於零的情況. 12.

(24) 3.2、基本有限差分式的推導. 有限差分的原理是差分的近似，而實際推導則以泰勒方程式進行。由泰勒級數展開式(3.2-1)及(3.2-2)聯立：. 1 (∆x) 2 f '' ( x0 ) 2! 1 1 + (∆x)3 f ''' ( x0 ) + (∆x) 4 f '''' ( x0 ) + 3! 4!. (3.2-1). 1 (∆x) 2 f '' ( x0 ) 2! 1 1 − (∆x)3 f ''' ( x0 ) + (∆x) 4 f '''' ( x0 ) + 3! 4!. (3.2-2). f ( x0 + ∆x) = f ( x0 ) + ∆x f ' ( x0 ) +. f ( x0 − ∆x) = f ( x0 ) − ∆x f ' ( x0 ) +. 推導一階微分近似，由(3.2-1)-(3.2.2)，可得(3.2-3). f ( x0 + ∆x) − f ( x0 − ∆x) = 2∆x f ' ( x0 ) +. 2 (∆x)3 f ''' ( x0 ) + 3!. (3.2-3). 並令 f ( x0 ) = u N , f ( x0 + ∆x) = u N +1 , f ( x0 − ∆x) = u N −1 ，可得(3.2-4). u N +1 − u N −1 = 2∆x u ' N +. 2 (∆x)3 u ''' N + 3!. (3.2-4). 並忽略高階項，可得(3.2-5). u'N = 其中 −. u N +1 − u N −1 1 − (∆x) 2 u ''' N 2∆x 3!. (3.2-5). 1 (∆x) 2 u ''' N 為誤差項。 3!. 推導二階微分近似，由 (1) + (2) (3.2-1)+(3.2.2)，可得(3.2-6). 13.

(25) f ( x0 + ∆x) + f ( x0 − ∆x) = 2 f ( x0 ) +. 2 (∆x) 2 f '' ( x0 ) 2!. 2 + (∆x) 4 f '''' ( x0 ) + 4!. (3.2-6). 並同上代入 f ( x0 ) = u N , f ( x0 + ∆x) = u N +1 , f ( x0 − ∆x) = u N −1 ，可得(3.2-7). u N +1 + u N −1 = 2 u N +. 2 2 (∆x) 2 u '' N + (∆x) 4 u '''' N + 2! 4!. (3.2-7). 並忽略高階項. u '' N = 其中 −. u N +1 + u N −1 − 2u N 2 − (∆x) 2 u '''' N 2 ∆x 4!. (3.2-8). 2 (∆x) 2 u '''' N 為誤差項。 4!. 由此可得有限差分一階、二階的公式，(3.2-9)、(3.2-10). u N +1 − u N −1 2 ∆x u − 2u N + u N −1 u′′N = N +1 ∆x 2 u ′N =. (3.2-9) (3.2-10). 由差分式的式子可知每個點都和左右的點有關聯，示意圖如圖. 3-4，每個位置的點都和左右的點有關聯，則必要須找到最左和最右的兩點的值？而最左和最右是自己令的，也就是解問題最重要的邊界值(B.C)，也就是在第 1 個點和第 N 個點是最為重要的，第 1 個點左邊的點和第 N 個點右邊的點是數值分析上最重要的邊界值。. u N −2. u N −1. uN. u N +1. 圖 3-4 一維有限差分關聯圖 14. uN +2.

(26) 3.3、二維有限差分推導. u(i , j +1). u(i −1, j ). u( i , j ). u(i +1, j ). u(i , j −1). 圖 3-5 二維有限差分示意圖一般，我們在如圖 3-5，在二維平面上，可得(3.3-1). ∂2 ∂2 )u( i , j ) + ∂x 2 ∂y 2 (3.3-1) u(i , j +1) u(i , j −1) u( i +1, j ) u(i −1, j ) 1 1 = + + + − 2u( i , j ) ( 2 + 2 ) ∆x 2 ∆x 2 ∆y 2 ∆y 2 ∆x ∆y. ∇ 2 u( i , j ) = (. 一般情況，我們皆可取 ∆x = ∆y = h ，則可得. ∇ 2 u( i , j ) =. u(i , j +1) + u(i , j −1) + u(i +1, j ) + u( i −1, j ) − 4u(i , j ) h2. (3.3-2). 二維有限差分設點的方式在本論文都是以行設點，在計算區域中以行為單位表示全部的場量，設點方式示意圖如圖 3-6，假設矩形為計算區域，在矩形範圍外就是邊界值，設點以行向量由左向右，以此方法表示全部的場量。. 15.

(27) u(1,1). u(1,6). u(3,1). u(3,6). 圖 3-6 二維有限差分設點方式以拉普拉斯算符(Laplacian operator). ∂2 ∂2 ( 2 + 2 )u = 0 為例，結合二維有限差分，每個點寫成一個等式，將 ∂y ∂x 每個點的等式寫在一起，寫成聯立矩陣方程式的形式，如式(3.3-3)，其中令 ∆x = ∆y = h ，令一行有 N 個點共有 M 行，大矩陣的大小為. (M×N)×(M×N)。每一行寫成一個 block，所以一個 block 為(N×N)大小的矩陣，如式(3.3-4)、(3.3-5)， u j 為第 j 行的場量，有 N 個元素，如式(3.3-6)， b j 為邊界值向量。. ⎡A ⎢I 1 ⎢ ⎢ h2 ⎢ ⎢ ⎢⎣. I A. I I. A I. ⎤ ⎡ u1 ⎤ ⎡ b1 ⎤ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥⋅⎢ ⎥=⎢ ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ I⎥ ⎢ ⎥ ⎢ ⎥ A ⎥⎦ ⎢⎣ uM ⎥⎦ ⎢⎣ bM ⎥⎦. 16. (3.3-3).

(28) ⎤ ⎡− 4 1 0 ⎥ ⎢ 1 −4 1 0 ⎥ ⎢ ⎥ A=⎢ ⎥ ⎢ 0 1 −4 1 ⎥ ⎢ ⎢⎣ 0 1 − 4⎥⎦ N× N ⎤ ⎡1 ⎢ 1 0 ⎥ ⎥ ⎢ ⎥ I=⎢ ⎥ ⎢ 1 ⎥ ⎢ 0 ⎢⎣ 1⎥⎦ N× N ⎡ c(1, j ) ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ，j=1,…,M uj = ⎢ ⎥ ⎢ ⎥ ⎢⎣c( N , j ) ⎥⎦ N ×1. 17. (3.3-4). (3.3-5). (3.3-6).

(29) 3.4 四點有限差分. 當場型變化較為快速時，則我們必須要更精密的計算，則我們可以採用取四點的有限差分。首先令 h = ∆x ，並如圖 3.7 定義取點方法. x0 − 2h. x0 − h. x0 + h. x0. x0 + 2h. 圖 3.7、四點有限差分示意圖. 並由各點的泰勒級數展開式，如下四式. 1 (h) 2 f '' ( x0 ) 2! 1 1 + (h)3 f ''' ( x0 ) + (h) 4 f '''' ( x0 ) + 3! 4!. (3.4-1). 1 (h) 2 f '' ( x0 ) 2! 1 1 − (h)3 f ''' ( x0 ) + (h) 4 f '''' ( x0 ) + 3! 4!. (3.4-2). f ( x0 + h) = f ( x0 ) + h f ' ( x0 ) +. f ( x0 − h) = f ( x0 ) − h f ' ( x0 ) +. 1 (2h) 2 f '' ( x0 ) 2! 1 1 + (2h)3 f ''' ( x0 ) + (2h)4 f '''' ( x0 ) + 3! 4!. f ( x0 + 2h) = f ( x0 ) + 2h f ' ( x0 ) +. 18. (3.4-3).

(30) 1 (2h) 2 f '' ( x0 ) 2! 1 1 − (2h)3 f ''' ( x0 ) + (2h) 4 f '''' ( x0 ) + 3! 4!. f ( x0 − 2h) = f ( x0 ) − 2h f ' ( x0 ) +. (3.4-4). 首先，推導一階微分近似；由(3.4-1)- (3.4-2)，可得(3.4-5). f ( x0 + h) − f ( x0 − h) = 2 3 ''' 2 (h) f ( x0 ) + (h)5 f ''''' ( x0 ) + 3! 5!. (3.4-5). 2 2 (2h)3 f (3) ( x0 ) + (2h)5 f (5) ( x0 ) + 3! 5!. (3.4-6). 2hf ' ( x0 ) + 由(3.4-3)- (3.4-4)，可得(3.4-6). f ( x0 + 2h) − f ( x0 − 2h) = 2(2h) f ' ( x0 ) +. 1 由(3.4-5)- (3.4-6)，可得(3.4-7) 8 8( f ( x0 + h) − f ( x0 − h)) − ( f ( x0 + 2h) − f ( x0 − 2h)) f ′( x0 ) = + O(h 4 ) 12h (3.4-7) 代入 f ( x0 ) = u N , f ( x0 + ∆x) = u N +1 , f ( x0 − ∆x) = u N −1 , f ( x0 + 2∆x) = u N + 2 ,. f ( x0 − 2∆x) = u N −2 ，可得一階微分的四點有限差分近似(3.4-8) u N′ =. 8(u N +1 − u N −1 ) − (u N + 2 − u N −2 ) 12h. (3.4-8). 推著，推導二階微分近似；由(3.4-1)- (3.4-2)，可得(3.4-9). f ( x0 + h) + f ( x0 − h) = 2 f ( x0 ) + (h)2 f '' ( x0 ) + 由(3.4-3)+(3.4-4)，可得(3.4-10). 19. 2 (h) 4 f (4) ( x0 ) + 4!. (3.4-9).

(31) f ( x0 + 2h) + f ( x0 − 2h) = 2 f ( x0 ) + (2h)2 f '' ( x0 ) + 由(3.4-9)-. 2 (2h)4 f (4) ( x0 ) + 4!. (3.4-10). 1 (3.4-10)，可得(3.4-11) 16. f ′′( x0 ) =. 4( f ( x0 + h) + f ( x0 − h)) 3h 2 f ( x0 + 2h) + f ( x0 − 2h) 5 f ( x0 ) − − + O(h 4 ) 2 2 12h 2h. (3.4-11). 代入 f ( x0 ) = u N , f ( x0 + ∆x) = u N +1 , f ( x0 − ∆x) = u N −1 , f ( x0 + 2∆x) = u N + 2 ,. f ( x0 − 2∆x) = u N −2 ，可得一階微分的四點有限差分近似(3.4-12) u N′′ =. 4(u N +1 + u N −1 ) u N + 2 + u N −2 5u N − − 2 + O(h 4 ) 2 2 3h 12h 2h. (3.4-12). 以下以四點求解有限差分來考慮：. Example:考慮如下的結構 ⎛ ∂2 2⎞ ⎜ 2 + k ⎟ Ey = 0 ⎝ ∂z ⎠. 給一維 TE wave eq.. 將 u 替代 E y ，寫成 FD 的形式，如式(3.4-13) ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 1 ⎜ 2 ⎜ ∆z ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝. ⎡ 5 ⎢− 2 ⎢ ⎢ 1 ⎢ 3 ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 0 ⎢ ⎢ ⎢ 0 ⎣. 1 3 5 − 2 −. 1 12. −. 1 12 1 3. 0 −. 1 12. 0. 5 2. 1 3. 1 12. 1 3 1 − 12. 1 3. −. 0. −. 0. −. 1 12 5 2 1 3. −. 20. ⎤ 0 ⎥ ⎥ 0 ⎥ ⎡k 2 ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥+⎢ ⎥ ⎢ ⎥ ⎢ ⎢ 1 ⎥⎥ ⎣⎢ 3 ⎥ 5⎥ − ⎥ 2⎦. ⎞ ⎟ ⎟ ⎟ ⎤⎟ ⎥⎟ ⎥⎟ ⎥⎟ ⎥⎟ = 0 ⎥⎟ ⎥⎟ ⎥⎟ k 2 ⎦⎥ ⎟ ⎟ ⎟ ⎟ ⎠.

(32) (3.4-13) 若為二維結構，則方法即為使用兩個一維有限差分，令 ∆z = ∆x = h ，如圖 3-8 取點，可化為式(3.4-14). u(i , j + 2). u(i , j +1). u(i −2, j ). u( i −1, j ). u( i , j ). u(i +1, j ). u(i + 2, j ). u(i , j −1). u(i , j −2). 圖 3.7、二維四點有限差分取點示意圖 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 1 ⎜ 2 ⎜h ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝. ⎡ ⎢A ⎢ ⎢1 I ⎢3 ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢0 ⎢ ⎢ ⎢0 ⎣. 1 I 3 A. −. 1 I 12. −. 1 I 12 1 I 3. 0 −. 1 I 3 0. −. 1 I 12. 0. A. 1 I 3. 1 I 12 0. 1 I 3 1 − I 12. −. 1 I 12. A 1 I 3. ⎤ 0⎥ ⎥ 0 ⎥ ⎡k 2 ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥+⎢ ⎥ ⎢ ⎥ ⎢ ⎢ 1 ⎥ ⎢⎣ I⎥ 3 ⎥ ⎥ A⎥ ⎦. ⎞ ⎟ ⎟ ⎟ ⎤⎟ ⎥⎟ ⎥⎟ ⎥⎟ ⎥⎟ = 0 ⎥⎟ ⎥⎟ ⎥⎟ k 2 ⎥⎦ ⎟ ⎟ ⎟ ⎟ ⎠. (3.4-14). 21.

(33) 1 ⎡ ⎢ −5 3 ⎢ ⎢1 −5 ⎢3 ⎢ ⎢ ⎢ 1 A 即為 ⎢ − 12 ⎢ ⎢ ⎢ ⎢0 ⎢ ⎢ ⎢0 ⎣. −. 1 12 1 3 1 3 0. 0 1 12. 0. −5. 1 3. −. −. 1 12 0. 1 3 1 − 12. −. 1 12. −5 1 3. 22. ⎤ 0⎥ ⎥ 0⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ 1⎥ ⎥ 3⎥ ⎥ −5⎥ ⎦.

(34) 第四章一維有限差分計算推導為了解各種複雜結構的光波導，由於結構不是固定的結構，所以我們使用的數值方法為有限差分法，在一維不連續邊界精確度探討我們解決了下面的問題：. 1、推導大對比材料有效差分式 2、邊界值處理文中使用的波方程式和一般波方程式不太一樣，式子如下. (4-1)、(4-2)，波方程式中以(4-1)為例，一階微分. 項，二階微分. ∂ E y ∝ H x 為連續 ∂z. ⎤ 1 ∂ ⎡1 ∂ E ∝ E y 為連續項，可知波方程式為連續 y ε r ∂z ⎢⎣ µ r ∂z ⎥⎦. 項；若是以有限差分的作法，不連續項的產生在等效介電常數，所以不連續項只有第一項和第二項，其中 p = k0 n sin θ 為常數。一維的介電常數和相對導磁係數隨行進方向變化 ε = ε ( z ) 。. TE wave 波方程式. 1. εr. (− p 2 ) E y +. 1 ∂2 E y + k02 E y = 0 2 ε r ∂z. (4.0-1). TM wave 波方程式. 1. εr. (− p 2 ) H y +. ⎤ ∂ ⎡1 ∂ H y ⎥ + k02 H y = 0 ⎢ ∂z ⎣ ε r ∂z ⎦ 23. (4.0-2).

(35) εr0. ε r1 x. u0. u1. u2. z. Z=0. 圖 4-1 單一界面設點一維不連續介質使用的模型如圖 4-1，設點方式是將點設在界面上，其中 u 為場量，可為電場或磁場，TE wave 為電場連續使用 TE wave 波方程式，反之 TM wave 磁場連續使用 TM wave 波方程式。. 24.

(36) 4.1 邊界條件處理和推有限差分式在邊界條件處理我們並不是以 Dirichlet’s、Neumann’s 或是. impedance boundary condition，而是以解析解的形式來逹成理想無反射邊界條件，作法是令入射波 ui = e. − jβ z. + jβz 、反射波 u r = re 和穿透. + jβz 波 ut = te 為已知的場量。. 入射區場量為入射波和反射波組成 u 1 = e − jβ1z + re + jβ1z ，穿透區場量 u 2 = te − jβ2 z ，u0 為邊界左邊 ∆z 的位置，而 u 2 為邊界右邊 ∆z 的位置，其中傳播因子 w = e − jβ∆z 。. u0 = e + jβ1∆z + re − jβ1∆z = w1−1 + rw1. (4.1-1). u 2 = te − jβ2∆z = tw2. (4.1-2). 使用邊界條件. u1 = u 1 z =0 = u 2. z =0. =1+ r = t. (4.1-3). r = u1 − 1. (4.1-4). t = u1. (4.1-5). 代回原式可以消去 r 和 t ，消去 r 和 t 代回原方程式可減少未知數。. 25.

(37) u0. u1 u−′. u′z =0. u2 u+′. 圖 4-2 介面處理示意圖在式(4-3)、(4-4)處理第二項式子時，使用差分式近似，以 TM wave 為例子，在 TM wave 式(4-3)第二項式子。在邊界各取兩邊差分乘加權常數. 1. εr. 為 u′+、u′−，為第一次差分，如式(4-8)，再取 u′+、u′− 平均得 u′z =0，. 如式 (4-9) 。再左右各差分第二次、左右乘加權常數再平均，如式. (4-10)，處理的示意圖如圖 4-2。 ×為設的場量的位置，第一次差分為△，第二次差分為□，再將□ 取平均為○，也是第二項式子的差分式。 ⎧ 1 ⎛ u2 − u1 ⎞ ⎪ u+′ = ⎜ ε 2 ⎝ ∆z ⎟⎠ z =∆z / 2 1 ∂u1 ⎪ ⎨ ε r ∂z ⎪ 1 u −u ′ = ⎛⎜ 1 0 ⎞⎟ u − ⎪ ε1 ⎝ ∆z ⎠ z =−∆z / 2 ⎩ u′z =0 =. 1 ∂u1 ε r ∂z. ≈ z =0. (4.1-6). 1 1 u −u 1 u −u ( u+′ + u−′ ) = ⎛⎜ 2 1 ⎞⎟ + ⎛⎜ 1 0 ⎞⎟ 2 2ε 2 ⎝ ∆z ⎠ 2ε 1 ⎝ ∆z ⎠. (4.1-7). 1 ⎛ u ′z =∆z / 2 − u ′z =0 ⎞ 1 ⎛ u ′z =0 − u ′z =− ∆z / 2 ⎞ 1 ∂ ⎡ 1 ∂u1 ⎤ = ⎜ ⎟+ ⎜ ⎟ (4.1-8) ⎢ ⎥ ∆z / 2 µ r ∂z ⎣ ε r ∂z ⎦ 2µ 2 ⎝ ∆z / 2 ⎠ 2µ1 ⎝ ⎠ 26.

(38) 由於界面左右為不連續的介質，所以我們在處理第二項式子時，各在邊界上左右取差分和加權平均使場量近似連續。而第三項式子為連續量。第一項式子在界面為不連續量，只有這個式子處理較為麻煩，我們是取算數平均。. TM wave 的波方程式得到的是磁場比，若是要換成電場比表示，電場比和磁場比相差一個係數。反射率 rEP = − rHP 穿透率 t EP = t HP. (4.1-9). η N +1 η1. (4.1-10). 綜合以上討論的方法可推得方程式式(7.2-13)，其中 s = n sin ϑ ⎛ 1 1 ⎞ ⎟ − w1−1 )⎜⎜ + µ1ε 1 µ1ε 2 ⎟⎠ ⎝ u1 = ⎛ ⎛ 1 ⎛ 1 1 ⎞ ⎛ 1 1 1 1 ⎞ 1 ⎞ s2 ⎟⎟ + 2∆z 2 k 02 ⎜⎜1 − ⎟⎟ + w2 ⎜⎜ ⎟⎟ + ⎜⎜ − w1 ⎜⎜ + − − − + ⎝ ε r µr ⎝ µ 2ε 1 µ 2ε 2 ⎠ ⎝ µ1ε 1 µ1ε 2 ⎠ ⎝ µ1ε 1 µ1ε 2 µ 2ε 1 µ 2ε 2 ⎠. (w. 1. ⎞ ⎟⎟ ⎠. (4.1-11) 當然這也可以推廣到多層介質，入射區場量 u 1 = e − jβ1z + re + jβ1z 和穿透區場量 u N +1 = te − jβ N +1 ( z −d ) ，d 為全部介質的寛度，不包括入射和穿透區。點 u1 和 u N 分別設在左右界面上。由左邊界 u1 = u 1 可得左邊界條件。右邊界 u N = u N +1. z =d. 樣的方法可得矩陣方程式，形式如下. 27. z =0. =1+ r ，. = t ，可得右邊界條件。依同.

(39) ( A + B + C )u = D. (4-14). 若對應波方程式，A 和 C 為 diagonal，B 為 tridiagonal，D 為向量且只有元素 1 有值，其中 B 為複數對稱矩陣. B=. A=C=. 以單一界面為例子， ε = [1,4] ， µ = [1,100] ， λ = 3e − 2 ， ϑi = 60° 與正解比較，其中 r 和 t 為反射係數和穿透係數皆為電場比。 ∆z = λ / N ，N 為點數，設 N 為 1000，列於 Table 4-1。. Table 4-1 單一界面計算值 FD. TE wave 0.4300 + 0.0068i reflection abs 0.4300 coefficients angle 0.9041 1.4300 + 0.0068i transmission abs 1.4300 coefficients angle 0.2719 正解. TE wave. 0.4290 abs 0.4290 angle 0 1.4290 transmission abs 1.4290 coefficients angle 0 reflection coefficients. 4.2 三層介質一維有限差分法計算 28. TM wave 0.8181 - 0.0100i abs 0.8182 angle -0.7006 0.9093 + 0.0500i abs 0.9107 angle 3.1487 TM wave 0.8180 abs 0.8180 angle 0 0.9099 abs 0.9099 angle 0.

(40) 以兩個界面為例子：. ε = [1,4,1] µ = [1,25,1] λ = 3 ×10 −2 ϑi = 60° d = 1.1e − 2 寛度和波長比約 1：3。其中 r 和 t 為電場比，與 ABCD method 比較，列 Table 4-2。. Table 4-2 三層介質計算值 FD（500p）. TE wave. 0.1546 + 0.1062i abs 0.1876 angle 34.4930 -0.5563 + 0.8096i transmission abs 0.9823 coefficients angle 124.4930 reflection coefficients. 正解. TE wave. 0.1523 + 0.1038i abs 0.1843 angle 34.2756 -0.5535 + 0.8122i transmission abs 0.9829 coefficients angle 124.2756 reflection coefficients. TM wave 0.8517 + 0.2419i abs 0.8853 angle 15.8536 -0.1270 + 0.4473i abs 0.4650 angle 105.8536 TM wave 0.8599 + 0.2320i abs 0.8907 angle 15.1010 -0.1185 + 0.4390i abs 0.4547 angle 105.1010. 將 TE wave 和 TM wave 在兩層界面各層的場量圖畫出，圖 4-3 和圖. 4-4 為 TE wave 的連續場量 Ey 和 Hx，圖 4-5 和圖 4-6 為 TM wave 的連續場量 Hy 和 Ex。. 29.

(41) TE wave, Ey. 場量介質寛度. 圖 4-3 TE wave, Ey TE wave, Hx. 場量. 介質寛度. 圖 4-4 TE wave, Hx TM wave, Hy. 場量. 介質寛度. 圖 4-5 TM wave, Hy. 30.

(42) TM wave, Ex. 場量介質寛度. 圖 4-6 TM wave, Ex 入射介質和穿透介質為空氣，中間層介質和空氣相差很大，折射率為 10，中間層的波長為 λ = 盪，波數為 N =. d. λ. λ0 n. = 3 × 10 −1 ，波在在中間層左右來回振. × 2 ，由此例子中間層波數約為 7。. 由以上的討論可知單一界面，一個波長放 50 點可得千分之一的誤差，如 Table 4-3，這個例子主要為探討在一個界面的誤差。三層介質，d 和真空波長 λ0 的比值約 λ0 : d = 3 : 1 ，設 500 點可得到約百分之一的誤差，如 table 7-4。. Table 4-3 二層介質誤差 TE t_abs. 0.0023. r_abs 2.4450e-004 TM. t_abs. 6.9979e-004 r_abs 8.7922e-004. 31.

(43) Table 4-4 三層介質誤差 TE t_abs t_angle. 6.1044e-004 r_abs 0.0017. r_angle. 0.0179 0.0063. TM t_abs. 0.0227. r_abs. 0.0061. t_angle. 0.0072. r_angle. 0.0498. 32.

(44) 第五章運用頻域有限差分求解波導場型. 一般二維結構較為簡單的波導，如矩形波導、圓柱波導等，模態皆有解析解；但對於結構較為複雜的波導，則可依賴頻域有限差分，求解出模態其誤差小且在可接受範圍的精確解。而因為在二維結構時，時域有很長的頻響的尾巴，所以時域和頻域在數值上無法做出精準的轉變（如 DFT, FFT），所以在二維波導結構需用頻域有限差分，方可得準確頻譜。而頻域有限差分，就是一種以離散的量等效連續量的方法，所以最佳化等效參數是必要的。而同處理一般的偏微分方程式，必須處理邊界條件，使其可以近似於無界的情況；而頻域有限差分將產生大小為計算點數的平方的大矩陣，如何去簡化計算的矩陣，亦是我們必須處理的問題。. 33.

(45) 5.1 計算區域有限差分式的近似. 首先求解二維的波動方程式的有限差分式，並定義座標如下： x. z. 首先列出波動方程式，如 TE 波式(5.1-1)，及 TM 波(5.1-2)：. TE：. 1 ∂ 1 ∂ 1 ∂ 1 ∂ [ [ Ey ] + E ] + k02 E y = 0 ε r ∂x µr ∂x ε r ∂z µr ∂z y. (5.1-1). TM：. 1 ∂ 1 ∂ 1 ∂ 1 ∂ [ [ Hy]+ H ] + k02 H y = 0 µr ∂x ε r ∂x µr ∂z ε r ∂z y. (5.1-2). 當所求為光波導時，其 µr = 1 ， ε r 隨波導形狀而變，代入(5.1-1) 及. (5.1-2)，則可得 TE：. 1 ∂2 1 ∂2 + E E y + k02 E y = 0 y 2 2 ε r ∂x ε r ∂z. (5.1-3). TM：. ∂ 1 ∂ ∂ 1 ∂ [ Hy]+ [ H y ] + k02 H y = 0 ∂x ε r ∂x ∂z ε r ∂z. (5.1-4). 則我們由此可求波動方程於光波導內的有限差分等效式。先求 TE 波，令 E y = u ，首先必須求出所需的 ε r 值，舉例如下：. n2. n1. 圖 5-1 任意光波導結構示意圖. 34.

(46) 而 ε r ( z, x) = n 2 ( z, x) 其中 k02 E y 為連續函數，所以. 1 ∂2 1 ∂2 E 、 E 亦為連 ε r ∂x 2 y ε r ∂z 2 y. 續函數。一般來說，可視 Ai B = Ai B ，即相乘後平均用平均後相乘來近似，故如. 1. εr. 、. ∂2 E y 可分別計算平均來近似波動方程式。而如上圖， ∂z 2. 此時 ε r 並非連續函數，故必須先行平均，方可計算。當取格點計算時，必須取出離散的值，其值以各種平均值為較接近原始函數的選擇。由於 ε r 於式中位在分母，所以我們使用倒數的平均，即取. 1. εr. =. 1 ， ε r ( z, x). 形式如式(5.1-5) 1. εr. ∆x 2 ∆x x0 − 2. =∫. x0 +. ∆z 2 ∆z z0 + 2. ∫. z0 +. 1 dzdx ε r ( z, x). (5.1-5). 得到 ε r 後，代入式(5.1-3)，可得(5.1-6). TE：. ∂2 ∂2 u + u + k 2u = 0，其中 k 2 = k02ε r µr = k02ε r ∂x 2 ∂z 2. (5.1-6). 並取點如下圖，則可得式(5.1-7) u(i , j +1) ∆z. 2. +. u(i , j −1) ∆z. 2. +. u(i +1, j ) ∆x. 2. +. u(i −1, j ) ∆x. 2. − 4u(i , j ) (. 1 1 + 2 ) + k 2u(i , j ) = 0 (5.1-7) 2 ∆z ∆x. u(i +1, j ). u(i , j −1). u( i , j ). u(i , j +1) x. u(i −1, j ). 圖 5-2 二維有限差分關聯圖. 可得 TE 波的有限差分式，如式(5.1-8) 35. z.

(47) 1 1 1 1 (u(i , j +1) + u(i , j −1) ) + 2 (u( i +1, j ) + u( i −1, j ) ) + (k 2 − 4( 2 + 2 ))u( i , j ) = 0 (5.1-8) 2 ∆z ∆x ∆z ∆x 2π n 2π 1 = i ，則式(5.1-8)成為式(5.1-9) 令 λ = N λ i∆x ，且 ∆x = ∆z ，又 k = N λ ∆x λ. u(i , j +1) + u(i , j −1) + u(i +1, j ) + u(i −1, j ) + ((. 2π n 2 ) − 4)u( i , j ) = 0 Nλ. (5.1-9). 接著求 TM 波：. TM：. ∂ 1 ∂ ∂ 1 ∂ [ Hy]+ [ H y ] + k02 H y = 0 ∂x ε r ∂x ∂z ε r ∂z. 其中在有限差分的計算，為簡化計算量，可取中點近似，即上下左右四點皆以 ε r 的值為值。或取格點等效，即取格點的平均 ε r 的值為值。設 H y = u 分述如下：若取中點近似，則可得式(5.1-10) ∂ 1 ∂ ∂ 1 u(i +1, j ) − u(i , j ) 1 u(i , j ) − u(i −1, j ) − [ ] Hy] = [ ∂z ε r ∂z ∂z ε r ∆z ∆z εl 1 u(i +1, j ) − u( i , j ) 1 u(i , j ) − u(i −1, j ) = − ∆z 2 ∆z 2 εr εl. (5.1-10). 其中當 ε r 和 ε l 同時出現時， ε r 代表右半格的 ε r (介電常數)值，同理可得式(5.1-11) ∂ 1 ∂ 1 u(i , j +1) − u(i , j ) 1 u( i , j ) − u(i , j −1) − [ Hy] = ∂x ε r ∂x ∆x 2 ∆x 2 εu εd. (5.1-11). 可得 TM 的有限差分式，式(5.1-12) 1 u( i +1, j ) − u(i , j ) 1 u(i , j ) − u(i −1, j ) − εr εl ∆z 2 ∆z 2 1 u(i , j +1) − u( i , j ) 1 u( i , j ) − u(i , j −1) + − + k02u( i , j ) = 0 2 2 εu εd ∆x ∆x. 1 1 1 1 u + u + u + u 2 ( i +1, j ) 2 ( i −1, j ) 2 ( i , j +1) ε r ∆z ε l ∆z ε u ∆x ε d ∆x 2 (i , j −1) 1 1 1 1 + (k − − − − )u = 0 2 2 2 ε r ∆z ε l ∆z ε u ∆x ε d ∆x 2 ( i , j ) 2 0. 36. (5.1-12). (5.1-13).

(48) 令 λ = N λ i∆x ，且 ∆x = ∆z，又 k0 = 1. εr. u(i +1, j ) +. 1. εl. u(i −1, j ) +. 1. εu. 2π. λ. =. 2π 1 i ，則式(5.1-13)可化為(5.1-15) N λ ∆x. u(i , j +1) +. 1. εd. u(i , j −1). 2π 1 1 1 1 + (( )2 − − − − )u(i , j ) = 0 Nλ εr εl εu ε d. (5.1-14). ε、 r ε、 l ε、 u ε d 的取得，可預先處理。如所取格點為 100 × 100 的矩陣，可先. 計算向上移半格的 ε 值共 101 格，則可得 ε 、 u ε d 的值皆 100 格，向左移半格的 ε 值共 101 格，則可得 ε、 r ε l 的值皆 100 格。如此則可精簡一半的計算量，如圖 5-3. ε clad. εu. εl. ε core. εr εr. εd. 圖 5-3 ε、 r ε、 l ε、 u ε d 的等效示意圖若取 ε r 的格點等效，可得式(5.1-15) ∂ 1 ∂ 1 ∂2 Hy ] = u [ ε r ∂z 2 ∂z ε r ∂z. (5.1-15). ε r 為所求點之 ε r 的格點近似，同理可得式(5.1-16) ∂ 1 ∂ 1 ∂2 Hy ] = u [ ε r ∂x 2 ∂x ε r ∂x. (5.1-16). 可得 TM 波的有限差分式式(5.1-17) ∂2 ∂2 u + u + k02ε r u = 0 ∂x 2 ∂z 2. (5.1-17). 37.

(49) 5.2 二維頻域有限分模型探討. 假設二維模型如圖 5-4 所示，令上下邊界為電牆或磁牆，作此假設是由於計算區域比波長大得多，電磁波傳遞會逐漸衰減至零；而左右邊界為吸波邊界條件，作法是一維結果的延伸，將全部空間畫分為有限差分計算區域和解析解區如圖，在有限差分計算區域界面左邊為入射場區、右邊為穿透場區，以解析解的形式表示，再由界面求得邊界值。. 解析解區. 頻域有限差分區. 解析解區. 上下邊界為電牆或磁牆或週期牆左右邊界為吸波邊界. x z. 圖 5-4 頻域有限差分二維模型示意圖而場量的型式，將和邊界條件直接相關。以下將探討在不同上下邊界時場量的假設法。二維模型求邊界值的作法和一維類似，先令解析解的形式再由邊界連續求得邊界值。其中上下邊界可採電牆、磁牆、週期牆三種型式，其中三種型式分別具有奇對稱、偶對稱、週期的性質；以下，我們分別探討，由不同的上下邊界推導出不同的左右吸收邊 38.

(50) 界；首先我們來看看幾何結構，如圖 5-5，電牆,磁牆或週期牆 ∆x / 2. u NM. u11. u1M. 吸收邊界. 吸收邊界. uN 1. ∆z. 電牆,磁牆或週期牆. 圖 5-5 頻域有限差分取點示意圖並依不同上下邊界，推導左右吸收邊界. <型一>上下邊界為電牆(EW)導出左右吸收邊界 <STEP1>以 Fourier series f ( x) = a0 + ∑ ⎛⎜ an cos N. n =1. ⎝. nπ nπ x + bn sin L L. ⎞ x⎟ ⎠. ⎛ ∂2 ⎞ ∂2 + + β 2 ⎟ u = 0 的解析解，可得 2 2 ⎝ ∂x ∂z ⎠. 求解 Maxwell eq. ⎜ 入射區場量: N. u i ( x, z ) = φ1 ( x)e − jβ1 z + ∑ rnφn ( x)e jβn z i. i. n =1. = φ1 ( x) ⎡e ⎣. j β1i z. N. − 2 j sin( β1i z ) ⎤ + ∑ rnφn ( x)e jβn z ⎦ n =1 i. N. (by e − j β1 z = e jβ1 z − 2 j sin( β1i z ) ) i. i. = −2 jφ1 ( x) sin( β1i z ) + φ1 ( x)e j β1 z + ∑ rnφn ( x)e jβn z i. i. n =1. N. ⎧1 + rn n = 1 (by an = ⎨ ) ⎩ rn n ≠ 1. = −2 jφ1 ( x) sin( β1i z ) + ∑ anφn ( x)e jβn z i. n =1. (5.2-1) u 可為 E 或 H，例如: E i ( x, z ) = Ex e− jβ z i 1. 原先假設 u i ( x, z ) 有單一入射和多重反射波的線性組合，而經過數學上 39.

(51) 的化簡，變成單一駐波和多重反射波的線性組合。穿透區場量: N. u t ( x, z ) = ∑ tnφn ( x)e − j βn ( z − Lz ). (5.2-2). t. n =1. ⎧ 2 nπ φn ( x) = sin x (n = 1… N ) ⎪ Lx Lx ⎪ 其中 ⎨ 2 2 ⎛ nπ ⎞ ⎛ nπ ⎞ ⎪ i 2 2 = − = − β ω µε k ⎜ ⎟ ⎜ ⎟ (n = 1… N ) ⎪ n Lx Lx ⎝ ⎠ ⎝ ⎠ ⎩. 當上下邊界為電牆. <STEP2>如何做有限差分 uN 1 u( N −1)0 u( N −1)1 u( N − 2)0 u( N − 2)1 u10. u11. ∆x / 2 ∆x / 2. u NM u( N −1) M. u( N −1)( M +1). u1M. u1( M +1). ∆x. ∆x / 2. u01. x z. u0M. 圖 5-6 頻域有限差分取點示意圖二. <STEP2_1>推導左吸波邊界 N. 由式(5.2-1)， u i ( x, z ) = −2 jφ1 ( x) sin( β1i z ) + ∑ anφn ( x)e jβ z ，並令 i n. n =1. ⎡ u10 ⎤ u0 = ⎢⎢ ⎥⎥ ⎢⎣u N 0 ⎥⎦ N ×1. (5.2-3). 40.

(52) ⎡ u11 ⎤ N u1 = ⎢⎢ ⎥⎥ = u i ( xk , z = 0) = ∑ anφn ( xk ) ( by k = 1, 2 n =1 ⎢⎣u N 1 ⎥⎦ N ×1 φN ( x1 ) ⎤ ⎡ φ1 ( x1 ) ⎡ a1 ⎤ ⎢ ⎥ ⎢ ⎥ = T ia =⎢ 1 ⎥ ⎢ ⎥ ⎢⎣φ1 ( xN ) ⎥ ⎢ ⎥ φN ( xN ) ⎦ N × N ⎣ aN ⎦ N ×1. N). (5.2-4). <NOTE> 在 z=0 處存在解析解和有限差分 u11 = a1φ1 ( x1 ) + a2φ2 ( x1 ) +. + aN φN ( x1 ) 位置 x1. u21 = a1φ1 ( x2 ) + a2φ2 ( x2 ) +. + aN φN ( x2 ) 位置 x2. 以此類推，並定義 T1 、 a 如下. ⎡ φ1 ( x1 ) T1 = ⎢⎢ ⎢⎣φ1 ( xN ). φN ( x1 ) ⎤. ⎥ ⎥ φN ( xN ) ⎥⎦ N × N. ⎡ 2 π sin x1 ⎢ Lx Lx ⎢ =⎢ ⎢ ⎢ 2 sin π x N ⎢⎣ Lx Lx. 2 Nπ ⎤ sin x1 ⎥ Lx Lx ⎥ ⎥ ⎥ 2 Nπ ⎥ sin xN ⎥ Lx Lx ⎦ N×N. (5.2-5) ⎡ a1 ⎤ a = ⎢⎢ ⎥⎥ ⎢⎣ aN ⎥⎦ N ×1. (5.2-6). <STEP2_1_1>求出 a an =. u i ( xk , z = 0) φn ( xk ). φn ( xk ) φn ( xk ). ⎡ a1 ⎤ ⎡ u1 φ1 ( x1 ) ⎢ ⇒ a = ⎢⎢ ⎥⎥ = ⎢ ⎢⎣ aN ⎥⎦ ⎢⎢ u1 φN ( xk ) ⎣. =. u1 φn ( xk ). φn ( xk ) φn ( xk ). ⎤ ⎡ φ (x ) ⎥ ⎢ 1 1 ⎥=⎢ ⎥ ⎢φ ( x ) ⎥⎦ ⎣ N 1. <NOTE> 符號<A|B>為 A 和 B 積分的意思 41. = u1 φn ( xk ). φ1 ( xN ) ⎤ ⎡ u11 ⎤. (5.2-7). ⎥ ⎢ ⎥ = T T iu (5.2-8) ⎥⎢ ⎥ 1 1 φN ( xN ) ⎥⎦ ⎢⎣u N 1 ⎥⎦.

(53) u1 φ1 ( x1 ) = u11φ1 ( x1 ) + u21φ1 ( x2 ) +. + u N 1φ1 ( xN ). u1 φ2 ( x1 ) = u11φ2 ( x1 ) + u21φ2 ( x2 ) +. + u N 1φ2 ( xN ). 以此類推. <STEP2_1_2>求出 u0 ⎡ u10 ⎤ ⎥ = u i ( x , z = −∆z ) u0 = ⎢⎢ k ⎥ ⎣⎢u N 0 ⎦⎥ N. = 2 jφ1 ( xk ) sin( β ∆z ) + ∑ anφn ( xk )e − jβn ∆z i 1. i. ( by. k = 1, 2. N). n =1. − j β1 ∆z ⎡ φ1 ( x1 ) ⎤ ⎡ φ1 ( x1 )e ⎢φ (x ) ⎥ ⎢ − j β1i ∆z ⎢ x e ( ) φ 1 2 i 1 2 ⎥+ = 2 j sin( β1 ∆z ) ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ i x ( ) φ ⎣ 1 N ⎦ ⎢⎣φ1 ( xN )e − j β1 ∆z = 2 j sin( β1i ∆z )iφ1 + T1 i Pi iT1T iu1 i. φ2 ( x1 )e − jβ ∆z i 2. φ2 ( x2 )e −. j β 2i ∆z. φ2 ( xN )e− j β ∆z i 2. ⎤⎡a ⎤ ⎥⎢ 1 ⎥ − ⎥ ⎢ a2 ⎥ φN ( x2 )e ⎥⎢ ⎥ ⎥⎢ ⎥ − j β Ni ∆z ⎥ a φ N ( xN )e ⎦⎣ N ⎦. φN ( x1 )e− j β. i N ∆z. j β Ni ∆z. 其中 T1 i Pi iT1T 即為 ABC(absorption boundary condition)，吸波邊界條件。. (5.2-9) <NOTE> 此處的 a 和 STEP2_1_1 一樣，等於 T1T iu1 ，並定義 Pi 、 φ1 為 ⎡e − jβ1 ∆z ⎢ Pi = ⎢ ⎢ ⎣ 0 i. ⎤ ⎥ ⎥ i ⎥ e − jβ N ∆z ⎦ 0. (5.2-10). ⎡ φ1 ( x1 ) ⎤ ⎢φ (x ) ⎥ φ1 = ⎢ 1 2 ⎥ ⎢ ⎥ ⎢ ⎥ ⎣φ1 ( xN ) ⎦. (5.2-11). <STEP2_2>推導右吸波邊界 42.

(54) N. 由(5.2-2)， u t ( x, z ) = ∑ tnφn ( x)e− jβ ( z − Lz ) 式，並定義 Pi 、 φ1 為 t n. n =1. ⎡u(1)( M +1) ⎤ ⎢ ⎥ uM +1 = ⎢ ⎥ ⎢ u N ( M +1) ⎥ ⎣ ⎦. (5.2-12). ⎡ u1M ⎤ N ⎥ = u t ( x , z = Lz ) = t φ ( x ) uM = ⎢⎢ ∑ k n n k ⎥ n =1 ⎢⎣u NM ⎥⎦ φN ( x1 ) ⎤ ⎡ t1 ⎤ ⎡ φ1 ( x1 ) ⎢ ⎥ ⎢ ⎥ = T it =⎢ ⎥⎢ ⎥ 1 ⎢⎣φ1 ( xN ) φN ( xN ) ⎥⎦ ⎢⎣t N ⎥⎦. ( by. k = 1, 2. N). (5.2-13). <STEP2_2_1>求出 t tn =. u t ( xk , z = Lz ) φn ( xk ). =. φn ( xk ) φn ( xk ). ⎡ t1 ⎤ ⎡ uM φ1 ( x1 ) ⎢ ⇒ t = ⎢⎢ ⎥⎥ = ⎢ ⎢⎣t N ⎥⎦ ⎢⎢ uM φN ( xk ) ⎣. uM φn ( xk ). φn ( xk ) φn ( xk ). ⎤ ⎡ φ (x ) ⎥ ⎢ 1 1 ⎥=⎢ ⎥ ⎢φ ( x ) ⎥⎦ ⎣ N 1. = uM φn ( xk ). (5.2-14). φ1 ( xN ) ⎤ ⎡ u1M ⎤. ⎥⎢ ⎥ = T T iu ⎥⎢ ⎥ 1 M φN ( xN ) ⎥⎦ ⎢⎣u NM ⎥⎦. <STEP2_2_2>求出 uM +1 ⎡u(1)( M +1) ⎤ N t ⎢ ⎥ t ( , ) uM +1 = ⎢ u x z Lz z tnφn ( xk )e − j βn ∆z = = + ∆ = ∑ k ⎥ n =1 ⎢ u N ( M +1) ⎥ ⎣ ⎦ ⎡ φ1 ( x1 )e − j β1 ∆z φ2 ( x1 )e− jβ2 ∆z ⎢ t t ⎢ φ1 ( x2 )e − jβ1 ∆z φ2 ( x2 )e− jβ2 ∆z =⎢ ⎢ ⎢φ ( x )e − jβ1t ∆z φ ( x )e − jβ2t ∆z N 2 ⎣ 1 N t. t. ( by. k = 1, 2. N). ⎤⎡t ⎤ ⎥ 1⎥ − j β Nt ∆z ⎢ ⎥ ⎢ t2 ⎥ φN ( x2 )e ⎥⎢ ⎥ ⎥⎢ ⎥ − j β Nt ∆z ⎥ t φ N ( x N )e ⎦⎣ N⎦. φN ( x1 )e− jβ. t N ∆z. = T1 i Pt iT1T iuM. (5.2-15) <NOTE> 此處的 t 和 STEP2_1_1 一樣，等於 T1T iuM ，並定義 Pt 、 t. 43.

(55) ⎡e − jβ1 ∆z ⎢ Pt = ⎢ ⎢ ⎣ 0. ⎤ ⎥ ⎥ − j β Nt ∆z ⎥ e ⎦. t. 0. (5.2-16). ⎡ t1 ⎤ t = ⎢⎢ ⎥⎥ ⎢⎣t N ⎥⎦. (5.2-17). <型二>上下邊界為磁牆(MW)導出左右吸收邊界 <STEP1>以 Fourier series 求解 ⎛ ∂2 ⎞ ∂2 + + β 2 ⎟ u = 0 的解析解,可得 2 2 ⎝ ∂x ∂z ⎠. Maxwell eq. ⎜ 入射區場量:. N −1. u i ( x, z ) = φ0 ( x)e − jβ0 z + ∑ rnφn ( x)e j βn z i. i. n=0. N −1. = φ0 ( x) ⎡e j β0 z − 2 j sin( β 0i z ) ⎤ + ∑ rnφn ( x)e jβn z ⎣ ⎦ n=0 i. i. N −1. (by e − jβ0 z = e j β0 z − 2 j sin( β 0i z ) ) i. i. = −2 jφ0 ( x) sin( β 0i z ) + φ0 ( x)e jβ0 z + ∑ rnφn ( x)e j βn z i. i. n=0. N −1. ⎧1 + rn n = 0 (by an = ⎨ ⎩ rn n ≠ 0. = −2 jφ0 ( x) sin( β 0i z ) + ∑ anφn ( x)e jβn z i. n=0. ). (5.2-18) 穿透區場量: N −1. u t ( x, z ) = ∑ tnφn ( x)e − j βn ( z − Lz ). (5.2-19). t. n=0. 其中. ⎧ 1 (n = 0) ⎪ Lx ⎪ φn ( x) = ⎪⎪ 2 nπ cos x (n = 1… N − 1) ⎨ Lx Lx ⎪ 2 2 ⎪ nπ ⎞ ⎛ nπ ⎞ 2 ⎪ β ni = ω 2 µε − ⎛⎜ k = − ⎟ ⎜ ⎟ (n = 0… N − 1) ⎪⎩ ⎝ Lx ⎠ ⎝ Lx ⎠. 當上下邊界為磁牆 44.

(56) <STEP2>如何做有限差分圖型和<型一>相同. <STEP2_1>推導左吸波邊界 N −1. 由式(5.2-18)， u i ( x, z ) = −2 jφ0 ( x) sin( β 0i z ) + ∑ anφn ( x)e jβ z ，並令 i n. n =0. ⎡ u10 ⎤ u0 = ⎢⎢ ⎥⎥ ⎢⎣u N 0 ⎥⎦ N ×1. (5.2-20). ⎡ u11 ⎤ N −1 u1 = ⎢⎢ ⎥⎥ = u i ( xk , z = 0) = ∑ anφn ( xk ) ( k = 1, 2 N ) n=0 ⎢⎣u N 1 ⎥⎦ N ×1 φN −1 ( x1 ) ⎤ ⎡ φ0 ( x1 ) ⎡ a0 ⎤ ⎢ ⎥ ⎢ ⎥ = T ia =⎢ 1 ⎥ ⎢ ⎥ ⎢⎣φ0 ( xN ) ⎥ ⎢ ⎥ φN −1 ( xN ) ⎦ N × N ⎣ aN −1 ⎦ N ×1. (5.2-21). <STEP2_1_1>求出 a an =. u i ( xk , z = 0) φn ( xk ). φn ( xk ) φn ( xk ). =. ⎡ a0 ⎤ ⎡ u1 φ0 ( x1 ) ⎥=⎢ ⇒ a = ⎢⎢ ⎥ ⎢ ⎢⎣ aN −1 ⎥⎦ ⎢⎢ u1 φN −1 ( xk ) ⎣. u1 φn ( xk ). φn ( xk ) φn ( xk ) ⎤ ⎡ φ (x ) ⎥ ⎢ 0 1 ⎥=⎢ ⎥ ⎢φ ( x ) ⎥⎦ ⎣ N −1 1. <STEP2_1_2>求出 u0. 45. = u1 φn ( xk ). φ0 ( xN ) ⎤ ⎡ u11 ⎤. ⎥⎢ ⎥ T ⎥ ⎢ ⎥ = T1 iu1 φN −1 ( xN ) ⎥⎦ ⎢⎣u N 1 ⎥⎦. (5.2-22).

(57) ⎡ u10 ⎤ ⎥ = u i ( x , z = −∆z ) u0 = ⎢⎢ k ⎥ ⎢⎣u N 0 ⎥⎦ N −1. = 2 jφ0 ( xk ) sin( β 0i ∆z ) + ∑ anφn ( xk )e − jβn ∆z. ( by. i. k = 1, 2. N). n=0. − j β 0 ∆z ⎡ φ0 ( x1 ) ⎤ ⎡ φ0 ( x1 )e i ⎢φ (x ) ⎥ ⎢ ⎢ φ0 ( x2 )e − j β0 ∆z 0 2 ⎥ i ⎢ = 2 j sin( β 0 ∆z ) + ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ i ⎣φ0 ( xN ) ⎦ ⎢⎣φ0 ( xN )e − jβ0 ∆z = 2 j sin( β 0i ∆z )iφ0 + T1 i Pi iT1T iu1. φ1 ( x1 )e− j β ∆z. i. i 1. φ1 ( x2 )e − j β ∆z i 1. φ1 ( xN )e− j β ∆z i 1. ⎤⎡ a ⎤ ⎥ 0 ⎥ − j β Ni −1∆z ⎢ ⎥ ⎢ a1 ⎥ φN −1 ( x2 )e ⎥⎢ ⎥ ⎥⎢ ⎥ i φN −1 ( xN )e− jβ N −1∆z ⎥⎦ ⎣ aN −1 ⎦. φN −1 ( x1 )e− jβ. i N −1∆z. (5.2-23) <NOTE> 此處的 a 和 STEP2_1_1 一樣，等於 T1T iu1 ，並令 Pi 、 φ0 ⎡ e − j β 0 ∆z ⎢ Pi = ⎢ ⎢ ⎣ 0. ⎤ ⎥ ⎥ − j β Ni −1∆z ⎥ e ⎦. i. 0. (5.2-24). ⎡ φ0 ( x1 ) ⎤ ⎢φ (x ) ⎥ φ0 = ⎢ 0 2 ⎥ ⎢ ⎥ ⎢ ⎥ ⎣φ0 ( xN ) ⎦. (5.2-25). ⎡ φ0 ( x1 ) T1 = ⎢⎢ ⎢⎣φ0 ( xN ) ⎡ 1 ⎢ ⎢ Lx =⎢ ⎢ ⎢ 1 ⎢⎣ Lx. φN −1 ( x1 ) ⎤. ⎥ ⎥ φN −1 ( xN ) ⎥⎦ N × N 2 ( N − 1)π ⎤ cos x1 ⎥ Lx Lx ⎥ ⎥ ⎥ 2 ( N − 1)π ⎥ cos xN ⎥ Lx Lx ⎦ N×N. πx 2 cos 1 Lx Lx πx 2 cos N Lx Lx. <STEP2_2>推導右吸波邊界 N −1. 由式(5.2-19)， u t ( x, z ) = ∑ tnφn ( x)e− jβ ( z − Lz ) ，令 t n. n=0. 46. (5.2-26).

(58) ⎡u(1)( M +1) ⎤ ⎢ ⎥ uM +1 = ⎢ ⎥ ⎢ u N ( M +1) ⎥ ⎣ ⎦. (5.2-27). ⎡ u1M ⎤ N −1 ⎥ = u t ( x , z = Lz ) = t φ ( x ) uM = ⎢⎢ ∑ k n n k ⎥ n =0 ⎢⎣u NM ⎥⎦ φN −1 ( x1 ) ⎤ ⎡ t0 ⎤ ⎡ φ0 ( x1 ) ⎢ ⎥⎢ ⎥ = T it =⎢ ⎥⎢ ⎥ 1 ⎢⎣φ0 ( xN ) ⎥ ⎢ φN −1 ( xN ) ⎦ ⎣t N −1 ⎥⎦. ( by. k = 1, 2. N). (5.2-28). <STEP2_2_1>求出 t tn =. u t ( xk , z = Lz ) φn ( xk ). φn ( xk ) φn ( xk ). uM φn ( xk ). =. φn ( xk ) φn ( xk ). ⎡ t0 ⎤ ⎡ uM φ0 ( x1 ) ⎢ ⇒ t = ⎢⎢ ⎥⎥ = ⎢ ⎢⎣t N −1 ⎥⎦ ⎢⎢ uM φN −1 ( xk ) ⎣. = uM φn ( xk ). ⎤ ⎡ φ (x ) ⎥ ⎢ 0 1 ⎥=⎢ ⎥ ⎢φ ( x ) ⎥⎦ ⎣ N −1 1. (5.2-29). φ0 ( xN ) ⎤ ⎡ u1M ⎤. ⎥⎢ ⎥ = T T iu ⎥⎢ ⎥ 1 M ⎥ ⎢ φN −1 ( xN ) ⎦ ⎣u NM ⎥⎦. <STEP2_2_2>求出 uM +1 uM +1. ⎡u(1)( M +1) ⎤ ⎢ ⎥ =⎢ ⎥ ⎢ u N ( M +1) ⎥ ⎣ ⎦ N −1. = u t ( xk , z = Lz + ∆z ) = ∑ tnφn ( xk )e − jβn ∆z t. ( by. k = 1, 2. N). n=0. ⎡ φ0 ( x1 )e − j β0 ∆z ⎢ t ⎢ φ0 ( x2 )e − j β0 ∆z =⎢ ⎢ ⎢φ ( x )e − j β0t ∆z ⎣ 0 N t. φ1 ( x1 )e− j β ∆z t 1. φ1 ( x2 )e− jβ ∆z t 1. φ1 ( xN )e − jβ ∆z t 1. ⎤⎡ t ⎤ 0 ⎥ ⎥ − j β Nt −1∆z ⎢ ⎥ ⎢ t1 ⎥ φN −1 ( x2 )e ⎥⎢ ⎥ ⎥⎢ ⎥ t φN −1 ( xN )e − j β N −1∆z ⎥⎦ ⎣t N −1 ⎦. φN −1 ( x1 )e− j β. t N −1∆z. = T1 i Pt iT1T iuM. <NOTE> 此處的 t 和 STEP2_1_1 一樣，等於 T1T iuM ，並定義 Pt 、 t. 47. (5.2-30).

(59) ⎡ e − jβ0 ∆z ⎢ Pt = ⎢ ⎢ ⎣ 0. ⎤ ⎥ ⎥ − j β Nt −1∆z ⎥ e ⎦. t. 0. (5.2-31). ⎡ t0 ⎤ t = ⎢⎢ ⎥⎥ ⎢⎣t N −1 ⎥⎦. (5.2-32). <型三>上下邊界為週期牆(PW)導出左右吸收邊界 <STEP1>以 Fourier series 求解 ⎛ ∂2 ⎞ ∂2 + 2 + β 2 ⎟ u = 0 的解析解，可得 2 ⎝ ∂x ∂z ⎠. Maxwell eq. ⎜ 入射區場量:. N −1. u i ( x, z ) = φ0 ( x)e − j β0 z + ∑ rnφn ( x)e jβn z i. i. n=0. = φ0 ( x) ⎡ e ⎣. j β 0i z. N −1. − 2 j sin( β 0i z ) ⎤ + ∑ rnφn ( x)e j βn z ⎦ n=0 i. N −1. (by e − jβ0 z = e j β0 z − 2 j sin( β 0i z )) i. i. = −2 jφ0 ( x) sin( β 0i z ) + φ0 ( x)e jβ0 z + ∑ rnφn ( x)e j βn z i. i. n=0. N −1. = −2 jφ0 ( x) sin( β 0i z ) + ∑ anφn ( x)e j βn z i. n =0. ⎧1 + rn n = 0 (by an = ⎨ ) ⎩ rn n ≠ 0. (5.2-33) 穿透區場量: N −1. u t ( x, z ) = ∑ tnφn ( x)e − j βn ( z − Lz ). (5.2-34). t. n=0. ⎧ 1 (n = 0) ⎪ Lx ⎪φn ( x) = ⎪ 2 2nπ 2 2nπ x+ x (n = 1… N − 1) cos sin 其中 ⎪⎨ Lx Lx Lx Lx ⎪ 2 2 ⎪ nπ ⎞ nπ ⎞ ⎛ ⎛ i 2 2 ⎪ β n = ω µε − ⎜ ⎟ = k −⎜ ⎟ (n = 0… N − 1) ⎝ Lx ⎠ ⎝ Lx ⎠ ⎩⎪. 48.

(60) 當下邊界為週期牆以上用到許多的數學等式，以下列附錄補充之。附錄一：證明 φn ( x) φn ( x) = 1 2 nπ 2 nπ x dx (by φn ( x) = x ) sin 2 sin 0 0 Lx Lx Lx Lx 2 Lx 2 nπ 1 Lx 2nπ 1 − cos 2θ x dx x dx (by sin 2 θ = = = sin 1 − cos ) ∫ ∫ 0 0 Lx Lx Lx Lx 2 Lx. φn ( x) φn ( x) = ∫ φn 2 ( x) dx = ∫. Lx 1 ⎡ 2nπ x− = sin ⎢ Lx ⎣ Lx 2nπ. Lx. x = Lx. ⎤ x⎥ =1 ⎦ x =0. 附錄二: 證明 φn ( x) φm ( x) = 0 (m ≠ n) Lx. φn ( x) φm ( x) = ∫ φn ( x)φm ( x) dx = ∫ 0. Lx. 0. 2 nπ mπ x sin x dx sin Lx Lx Lx. nπ 2 x ) sin Lx Lx 2 Lx 1 ⎡ (n + m)π (n − m)π ⎤ x − cos x ⎥ dx = − ⎢cos ∫ 0 Lx Lx Lx 2⎣ ⎦ 1 (by sin α sin β = − [ cos(α + β ) − cos(α − β ) ] ) 2 (by φn ( x) =. x = Lx. Lx 1 ⎡ Lx (n + m)π (n − m)π ⎤ x− x⎥ =− sin sin ⎢ Lx ⎣ (n + m)π Lx Lx (n − m)π ⎦ x =0 ⎧ n = 1, 2,3 ⎧ n + m = int . (n ≠ m ⎨ ) ⎨ ⎩m = 1, 2,3 ⎩ n − m = int . =0. 49.

(61) 附錄三: 證明 T1T = T1−1 :. 由 T1T iT1 = I. ⎡ φ1 ( x1 ) T iT1 = ⎢⎢ ⎢⎣φN ( x1 ). ⎡ φ1 ( x1 ) ⎢ ⎢ ⎢⎣φ1 ( xN ). T 1. N ⎡ 2 ⎢ ∑ φ1 ( xn ) ⎢ n =1 ⎢ N ∑ φ ( x )φ ( x ) = ⎢ n =1 2 n 1 n ⎢ ⎢ ⎢N ⎢ ∑ φ ( x )φ ( x ) ⎢⎣ n =1 N n 1 n. φ1 ( xN ) ⎤. ⎥ ⎥ φN ( xN ) ⎥⎦ N × N N. ∑ φ1 ( xn )φ2 ( xn ) n =1. N. ∑φ n =1. 2. 2. ( xn ). 可得 T1T = T1−1 φN ( x1 ) ⎤. ⎥ ⎥ φN ( xN ) ⎥⎦ N × N ⎤ ( xn ) ⎥ n =1 ⎥ ⎡1 N ⎥ ⎢0 φ2 ( xn )φN ( xn ) ⎥ ∑ =⎢ n =1 ⎥ ⎢ ⎥ ⎢ ⎥ ⎣0 N 2 ⎥ φN ( xn ) ∑ ⎥⎦ N × N n =1 N. ∑ φ ( x )φ 1. 50. n. N. 0 1 0. 0⎤ ⎥ ⎥ 0⎥ ⎥ 1⎦ N ×N.

(62) 5.3 頻域有限差分大矩陣計算簡化法. 首先我們用以下方法化簡，用以快速求解頻域有限所會遇到的三對角線矩陣(Solve tridiagonal matrix)。考慮如下的矩陣 ⎡ D1 ⎢ ⎢ B2 ⎢ ⎢0 ⎢ ⎢ ⎢ ⎢ ⎢0 ⎣. A1. 0. D2. A2. B3. D3. A3 BM −1. DM −1. 0. BM. 0 ⎤ ⎡ x1 ⎤ ⎡ r1 ⎤ ⎥ ⎢ ⎥ ⎢ r ⎥ ⎥ x2 ⎥ ⎢ ⎢ 2 ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ = (5.3-1) ⎢ ⎥ ⎢ ⎥ ⎥ 0 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ xM −1 ⎥ ⎢ rM −1 ⎥ AM −1 ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ ⎢⎣ xM ⎥⎦ ( N ×M )×1 ⎢⎣ rM ⎥⎦ ( N ×M )×1 ⎥ DM ⎦ ( N ×M )×( N ×M ). N：為在 x 方向取點數 M：為在 y 方向取點數. D ：為一 tridiagonal(三對角線)的矩陣，大小為 N*N A ：為一 diagonal(單對角線)的矩陣，大小為 N*N B ：為一 diagonal(單對角線)的矩陣，大小為 N*N x ：為一待求向量，大小為 N*1. r ：為一向量，大小為 N*1 <STEP1>求出 A1 和 r1 ，令 D1new = I. (5.3-2). A1new = D1−1 • A1old. (5.3-3). r1new = D1−1 • r1old. (5.3-4) 51.

(63) 可得以下矩陣 ⎡I ⎢ ⎢ B2 ⎢ ⎢0 ⎢ ⎢ ⎢ ⎢ ⎢0 ⎣. A1. 0. D2. A2. B3. D3. 0 ⎤ ⎡ x1 ⎤ ⎡ r1 ⎤ ⎥ ⎢ ⎥ ⎢ r ⎥ ⎥ ⎢ x2 ⎥ ⎢ 2 ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ =⎢ (5.3-5) ⎢ ⎥ ⎥ ⎥ 0 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ xM −1 ⎥ ⎢ rM −1 ⎥ AM −1 ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ x r ⎢ ⎥ ⎢ M M ⎣ ⎦ ⎣ ⎦⎥ ( N ×M )×1 × × N M ( ) 1 DM ⎥⎦ ( N ×M )×( N ×M ). A3 BM −1. DM −1. 0. BM. <STEP2>推導通式第一列乘 − B2 加第二列，可得原第二列式(5.3-6) ⎡0 ⎣⎢. (D − B A ) 2. (. 1. A2. 0. 0 ⎤⎥ [ x2 ] = ⎡ r2 − B2 r1 ⎤ ⎣ ⎦ ⎦. (5.3-6). ). 第二列皆除以 D − B2 A1 可得，即可得. (5.3-7). D2 new = I. (. A2 new = D − B2 A1. ). −1. (5.3-8). • A2 old. (5.3-9). r2 new = D1−1 • (r2 old − B2 r1 ). 代回式(5.3-5)可得以下矩陣，即式(5.3-10) ⎡I ⎢ ⎢0 ⎢ ⎢0 ⎢ ⎢ ⎢ ⎢ ⎢0 ⎣. A1. 0. I. A2. B3. D3. A3 BM −1. DM −1. 0. BM. 0 ⎤ ⎡ x1 ⎤ ⎡ r1 ⎤ ⎥ ⎢ ⎥ ⎢ r ⎥ ⎥ ⎢ x2 ⎥ ⎢ 2 ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ =⎢ (5.3-10) ⎢ ⎥ ⎥ ⎥ 0 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ xM −1 ⎥ ⎢ rM −1 ⎥ AM −1 ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ x r ⎢ ⎥ ⎢ ⎣ ⎦ ⎣ ⎦⎥ ( N ×M )×1 M M ( ) 1 N M × × DM ⎥⎦ ( N ×M )×( N ×M ). 以此類推,可導出通式式(5.3-11) ⎧⎪ Ai = ( Di − Bi Ai −1 ) −1 Ai i = 2,3, ⎨ −1 ⎪⎩ri = ( Di − Bi Ai −1 ) (ri − Bi ri −1 ) 52. , M −1. (5.3-11).

(64) <STEP3>求出和 rM ，即式(5.3-12) rM = ( DM − BM AM −1 ) −1 (rM − BM rM −1 ). (5.3-12). <STEP4>最後矩陣形式 ⎡I ⎢ ⎢0 ⎢ ⎢0 ⎢ ⎢ ⎢ ⎢ ⎢0 ⎣. A1. 0. I. A2. 0. I. A3 0. I. 0. 0. 0 ⎤ ⎡ x1 ⎤ ⎡ r1 ⎤ ⎥ ⎢ ⎥ ⎢ r ⎥ ⎥ ⎢ x2 ⎥ ⎢ 2 ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ =⎢ (5.3-13) ⎢ ⎥ ⎥ ⎥ 0 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ xM −1 ⎥ ⎢ rM −1 ⎥ AM −1 ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ x r ⎢ ⎥ ⎢ ⎣ ⎦ ⎣ ⎦⎥ ( N ×M )×1 M M N M × × ( ) 1 I ⎥⎦ ( N ×M )×( N ×M ). <STEP5>求解過程 xM = rM. (5.3-14). xM −1 = rM −1 − AM −1 xM. (5.3-15). x1 = r1 − A1 x2. (5.3-16). 其次，可利用陣列相乘的方式，可減少計算時間。如下：符號：.* 來源：MATLAB 的計算方式，於 MATLAB 內可直接計算運用：對角矩陣於相乘時的計算時間定義：如下首先定義. 53.

(65) ⎡ d1 ⎢0 D=⎢ ⎢ ⎢ ⎣0. 0⎤ ⎥ ⎥ 0⎥ ⎥ dn ⎦. 0 d2 0. (5.3-17). D 為一大小為 n*n 的矩陣，d 代表 diagonal。. Cnm. ⎡ c11 c12 ⎢c = ⎢ 21 ⎢ ⎢ ⎣ cn1. c1m ⎤ ⎥ ⎥ ⎥ ⎥ cnm ⎦ n×m. (5.3-18). 為一大小為 n*m 的矩陣，c 代表 constant。. Cmn. ⎡ c11 c12 ⎢c = ⎢ 21 ⎢ ⎢ ⎣cm1. c1n ⎤ ⎥ ⎥ ⎥ ⎥ cmn ⎦ m×n. (5.3-19). 為一大小為 n*m 的矩陣，c 代表 constant。當左乘對角矩陣時，可得. D iCnm. ⎡ a1c11 ⎢a c = ⎢ 2 21 ⎢ ⎢ ⎣ an cn1. a1c12. a1c1m ⎤ ⎥ ⎥ ⎥ ⎥ an cnm ⎦ n×m. (5.3-20). 觀察可知，第一列皆為 a1 ，第二列皆為 a2 ，以此類推；則我們可定義. D iCnm. a1c1m ⎤ ⎡ a1c11 a1c12 ⎢a c ⎥ 2 21 ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ an cnm ⎦ n×m ⎣ an cn1 a1 ⎤ ⎡ a1 a1 ⎡ c11 c12 ⎢a ⎥ ⎢c 2 ⎢ ⎥ .* ⎢ 21 = ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ an ⎦ n×m ⎣ cn1 ⎣ an 54. c1m ⎤ ⎥ ⎥ ⎥ ⎥ cnm ⎦ n×m. (5.3-21).

(66) 當右乘對角矩陣，同理可得 an c1m ⎤ ⎡ a1c11 a2 c12 ⎢a c ⎥ 1 21 ⎢ ⎥ Cmn i D = ⎢ ⎥ ⎢ ⎥ an cmn ⎦ m×n ⎣ a1cn1 an ⎤ ⎡ a1 a2 ⎡ c11 c12 ⎢a ⎥ ⎢c 1 ⎢ ⎥ .* ⎢ 21 = ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ an ⎦ m×n ⎣cm1 ⎣ a1. c1n ⎤ ⎥ ⎥ ⎥ ⎥ cmn ⎦ m×n. (5.3-22). 一般矩陣相乘，同一點的值須經乘 n 次，加 n-1 次才可得，無論矩陣有多少個零。若其中之一為對角矩陣，則可經由轉換而成為陣列相乘；則同一點的值只須乘 1 次即可得，可省略極大的計算量。. 55.

(67) 5.4 空心波導數值模擬結果以空心波導的模擬，來驗證吸波邊界條件是否正確以下面的兩個例子驗證：使用參數，其介質為空氣，入射波長 λ = 3 ×10−2，計算寛度(x 方向) Lx = 5λ 長度(z 方向) Lz = 5λ，入射波為垂直正向入射，計算點數大小 100x100，一個波長約 20 點，上下邊界為電牆，如圖 5-7：. 圖 5-7、空心波導上下邊界電牆場形圖. 場量顏色尺度. 其中：前五個模態反射係數絕對值 0.0012、0.0000、0.0000、0.0000、0.0000，前五個模態穿透係數絕對值 1.0000、0.0000、0.0000、0.0000、0.0000。將上下邊界改為磁牆，如圖 5-8：. 56.

(68) 圖 5-8、空心波導上下邊界磁牆場形圖. 場量顏色尺度. 其中：前五個模態反射係數絕對值 0.0019、0.0000、0.0000、0.0000、0.0000，前五個模態穿透係數絕對值 1.0000、0.0000、0.0000、0.0000、0.0000。接著，我們模擬點波源於空心波導內的場型圖：使用參數，其介質為空氣，入射波長 λ = 3 ×10−2 ，計算寛度(x 方向) Lx = 8.33λ 長度(z 方向) Lz = 8.33λ ，波源位於波導中心，計算點數大小 250x250，一個波長約 30 點，上下邊界為電牆，如圖 5-9：. 圖 5-9、點波源空心波導上下邊界電牆場形圖接下來，將上下邊界改為磁牆，如圖 5-10：. 57.

(69) 圖 5-10、點波源空心波導上下邊界磁牆場形圖並可透過將兩種解（電牆、磁牆）相加可達到另一類理想吸波效果。如圖 5-11：. 圖 5-11、點波源空心波導上下邊界理想邊界場形圖. 58.

(70) 5.5 模擬介質波導場型在完成空心波導的場型模擬後，我們進一步試著進行介質波導的模擬，並以磁牆的來模擬對稱的結構。以下我們先進行簡單波導結構的模擬，並同樣以磁牆簡化計算。考慮 TE 波:假設磁牆(MW)，存在 E y , H x , H z , 假設為如圖 5-12: X X=d 2d. MW(PEC). X=0. Z. 圖 5-12 介質波導示意圖及其磁牆等效圖利用對稱性，假設磁牆，可將波導簡化成圖 5-12 之右圖. STEP1：By Maxwell eq. 介質不隨方向改變 ⎛ ∂2 ∂2 2⎞ ⎜ 2 + 2 + k ⎟ Ey = 0 ⎝ ∂x ∂z ⎠. (5.5-1). and 1 ∂ ⎧ ⎪ H x = jωµ ∂z E y ⎪ ∇ × E = − jωµ H ⇒ ⎨ ⎪H = 1 ∂ E ⎪⎩ z − jωµ ∂x y. (5.5-2). STEP2：假設場量 ⎧⎪ E y 0 = Ae− j ( p0 x + β z ) for x > d ⎨ − j ( px + β z ) + Ce− j ( − px + β z ) for 0 < x < d ⎪⎩ E y1 = Be where p0 2 + β 2 = k0 2 p 2 + β 2 = k12. 由上可得 59. (5.5-3).

(71) 1 ⎧ − jp0 Ae− j ( p0 x + β z ) H z0 = ⎪ − jϖµ0 ⎪ ⎨ ⎪ H = 1 ( − jpBe − j ( px + β z ) + jpCe− j ( − px + β z ) ) ⎪⎩ z1 − jϖµ1. (. ). (5.5-4). STEP3：考慮邊界條件 ⎧ H z 0 = H z1 and E y 0 = E y1 at x = d ⎪ ∂ 1 ⎨ ⎪ H z1 = − jωµ ∂x E y1 = 0 at x = 0 1 ⎩. (5.5-5). 即切電切磁連續，由上可得 p0 p ⎧ − jp0 d ) = ( Be − jpd − Ce jpd ) ⎪ H z 0 = H z1 |x = d ⇒ µ ( Ae µ1 0 ⎪⎪ E y 0 = E y1 |x = d ⇒ Ae− jp0 d = Be− jpd + Ce jpd ⎨ ⎪ H z1 = 0 | x = 0 ⇒ B = C ⎪ ⎪⎩. (5.5-6). 將(5.5-6) 的第三式代入前兩式可得 p ⎧ p0 − jp0 d ) = − 2 jB sin pd ⎪ ( Ae µ1 ⎨ µ0 − jp0 d ⎪ Ae = 2 B cos pd ⎩ p p ⇒ j 0 = tan pd. µ0. (5.5-7). µ1. STEP4：假設 β > k0. p0 = − jα t. µ0 = µ1. ⎧ β 2 − α t 2 = k0 2 ⇒ α t 2 = (k12 − k0 2 ) − p 2 ⎨ 2 2 2 ⎩ β + p = k1. (5.5-8). 令 X = pd V 2 = ( k12 − k0 2 ) d 2 αt d = V 2 − X 2. 得 X tan X = V 2 − X 2 得到 X tan X = V 2 − X 2 後，代入找根公式。假設一結構如下：使用參數為波核(core)參數 n=1.5，波覆(cladding)參數 n=1.0，計. 60.

(72) 算寛度(x 方向)Lx=6.67λ，計算長度(z 方向)Lz=6.67λ；入射波長. λ=0.03，入射波為垂直正向入射，波核寬度為 0.06。代入 X tan X = V 2 − X 2 ，如圖 5-13，其中垂真線相交點為假根，水平線也是。. 圖 5-13 介質波導正解之找根圖由圖可得共五個根。代回原式，求解，並繪場型如圖 5-14：. 圖 5-14 介質波導正解之場型圖. 61.

(73) 並由以上 5.1~5.4 各節推導，繪出介質波導頻域有限差分的場型解，如圖 5-15：. 圖 5-15 介質波導頻域有限差分解之場型圖比較兩圖，可了解頻域有限差分法配合吸波邊界已可有效模擬介質波導場型圖。. 62.

(74) 5.6 波核漸細波導數值模擬結果. 接著，我們由上一節可證明吸波邊界處理和有限差分式的作法；接下來考慮特殊的波導結構，給一波導結構波核漸細波導(Tarpered. Wave)，示意圖如圖 5-16，. nclad = 1. ncore = 1.5. 圖 5-16 Tapered Waveguide 示意圖使用參數為波核(core)參數 n=1.5，波覆(cladding)參數 n=1.0，計算寛度(x 方向)Lx=13.3λ，計算長度(z 方向)Lz=23.3λ；入射波長. λ=0.03，入射波為垂直正向入射，共使用 140000 個變數。寬邊波核寬度=3.99 λ ，窄邊波核寬度=1.33 λ ，並且中間的斜率變化為常數，且左右各有長度 2.5 λ 有等寬的波核。如圖 5-17，首先可將計算區域劃分為入射解析區、穿透解析區、計算有限差分區，並且可利用對稱的原理，將計算的區域化簡成一半；如此劃分法，可有效化減計算量。. 63.

(75) 有限差分區. 入射反. 穿透解. 射解析. 析區. 區. 有限差分區圖 5-17 Tapered Waveguide 區域劃分示意圖. 模擬結果如圖 5-18、圖 5-19、圖 5-20，分別為基模態、第二高階模態、第四高階模態場型圖。圖中藍色實線，代表解析區及有限差分區的界限，以及波覆和波核的界限。並透過 5.1 節的方法，可得到波覆及波核邊界的等效介電常數。. 圖 5-17 Tapered Waveguide 基模態場型圖. 64.

(76) 圖 5-18 Tapered Waveguide 第二高階模態場型圖. 圖 5-19 Tapered Waveguide 第四高階模態場型圖以基模態為例，計算反射及穿透係數如下，可觀察到其誤差皆極小，與理論結果是符合的。反射係數前五個模態絕對值. -0.0001、0.0065、0.0122、0.0012、0.0025 穿透係數前五個模態絕對值. 1.0191、0.0201、0.0379、0.0359、0.0255. 65.

(77) 如圖 5-20，最上層電牆將會百分之百反射入射的波；而現在希望其反射波可以極少(至少小於百分之一)，如圖所示，我們採用約三十至五十層的多層結構，其中 n = 1 − jα ( x) 為各層的折射率 n。現在我們試著找出最佳的 α ( x) 公式。. 多層結構約三十至五十層. 電牆. n = 1 − jα ( x). x z 圖 5-20 吸收區示意圖. 最佳的 α ( x) 公式必須符合如下： n′( x) = const n( x ). (1 − jα ( x))′ = const ，令 const = C1 (1 − jα ( x)) − jα ′( x) = C1 (1 − jα ( x)). − jα ′( x) ∫ (1 − jα ( x)) dx = C1 x +C2 1 ∫ (1 − jα ( x)) dα ( x) = C1 x + C2. 66.

(78) ln(1 − jα ( x)) = C1 x + C2 1 − jα ( x) = C2 eC1x. α ( x ) = C2 e C x − j 1. 如圖 5-17，於上下加吸收區，並令 α ( x) = eα 0 x − j ，α 0 =-j*0.002，共一百層；結果如圖 5-20；. 圖 5-21 Tapered Waveguid 加 PML 第四高階模態場型圖. 如圖 5-21，我們可以看到，在上下邊界還有會有高階模態存在，是以吸收區尚且無法完全吸收。. 67.

(79) 第六章結論. 一般我們對光波導的特性要求皆是要求其頻域特性，例如其適用頻段、高階模態特性等等；所以頻域有限差分特別適用於光波導的模擬分析。而頻域有限差分法，其所需處理的要點包括如何推導波動方程式的頻域有限差分式，如何處理等效其介電常數，其邊界值如何處理；而頻域有限差分必須處理 ( N x × N z ) 2 的大矩陣，如何求解大矩陣，即是頻域有限差分法的重要課題。在第二章“理論與公式之推導”，我們一步步推導並簡化得到所需使用的波動方程式；並於第三章“以有限差分法近似波動方程式”介紹有限差分的概念，以及將其代入波動方程式。在第四章“一維有限差分”，正式進行模擬分析，並著重在邊界值處理和有限差分式的推導。其使用的數值是有限差分法，以單一界面為例和正解的結果比較，在波長較密區以一個波長取 50 點進行模擬計算，誤差為千分之一。使用一維模型主要在於探討不連續介質分析的精確度，找出有效的有限差分係數，滿足對大對比、大角度、還有不同極化的電磁波，對於邊界條件的處理方法也有別於一般情形。於第五章“二維頻域有限差分模型推導”，找出最適用的等效介電常數，並找到最佳的邊界條件，結果發表於 2003 年的台灣光電科技 68.

(80) 研討會；其中，傳統上，電牆為令平行邊界之電場為零，而磁牆則無法處理。而在此吸波邊界的處理可完全處理電牆及磁牆，且在對稱性的情況下可將兩種解（電牆、磁牆）相加以達到另一類理想吸波效果。其中頻域有限差分法將產生極大的矩陣，而以 block tri-diagonal 矩陣法於求解場型時可有效減少計算量。吸波邊界亦使其反射量極少，並且不需設定額外的吸波區域；這些方法使我們可精確模擬空心波導、tapered WG 等的場型。. 69.