SA13006060_左元_光子晶体光纤中导模的仿真

光子晶体光纤中导模的仿真
左元SA13006060
June29,2014
1背景简介
1.1光子晶体光纤
光子晶体是目前一个热门的研究方向，吸引着越来越多的研究人员的关注。

光子晶体是一种周期结构，这种结构的折射率在空间分布上存在着周期性。

研究人员希望通过光子晶体这种材料控制光的光学性质，利用光子晶体，可以让特定频率的光实现完美的反射，或者让它们只在某个特定的方向传播。

光子晶体这种材料的这些性质，显示出它在激光工程、高速通信和计算等领域的潜在价值[1]。

图1显示的分别是一到三维的光子晶体示意图，从图中可以看到光子晶体材料的折射率在空间分布存在着固定的周期，这也是被称为“晶体”的一个原因。

类比于常规的晶体，晶格对波的散射性质可以知道，光子晶体对光也会有类似的性质。

类似于晶体的能量禁带概念，光子晶体也有光子带隙的概念。

光子带隙会阻止特定方向传播的特定频率的光[1]。

图1:光子晶体
利用光子晶体的特性，可以制作出光子晶体光纤，也叫微结构光纤。

图2是三种不同类型的光子晶体光纤，分别是布拉格光纤（一维光子带隙光纤）、二维光子带隙光纤、Holey光纤。

应用中用得最多的是后面两种，光子带隙光纤利用光子带隙对光进行约束，让光在纤芯的低折射率的孔洞中传播。

通常孔洞中导光材料是空气，因此可以有效的减少损耗、不希望的非线性特性以及其他不希望的特性。

另一种是折射率导光光子晶体光纤(第三种)，这种结构的光纤不是
利用光子带隙，而是利用这种周期结构形成低有效折射率，而纤芯是高折射率材料，从而可以利用全内反射进行导光，将光约束在纤芯中[1]。

图2:光子晶体光纤
1.2时域有限差分(FDTD)方法
1966年，Yee首先提出麦克斯韦方程的时域有限差分求解方法，用来处理电磁场的传播和反射问题[2]。

通过将微分方程离散化，利用数值方法求解方程的数值解。

此后该方法得到进一步发展，1981年，Mur提出了在计算区域界断边界处的一阶和二阶吸收边界条件(ABC)[3]。

后来，Berenger进一步提出了用完全匹配层(PML)做吸收边界条件[4]。

这两种吸收边界条件被广泛地利用到仿真当中。

由于FDTD方法是直接从时域求解麦克斯韦方程，随时间的推进可以方便的给出电磁场随时间的演化过程。

在计算机上利用伪彩色显示这种动态过程，这种可视化清楚的显示了物理过程，也有助于对物理问题的分析和理解。

本文在仿真过程中也采用了这种可视化方法，得到了电场随时间演化的视频，作为附件放在邮件之中。

2理论分析
2.1光子晶体光纤的数值模型
本文将采用FDTD方法对光子晶体光纤的导模场进行仿真。

光子晶体光纤的横截面介电常数分布如图3所示，深色的表示空气，光在纤芯的空气层中传导。

光在介质中的传播行为可以用介质中的麦克斯韦方程来描述，其中两个旋度方程为
−σm H−µ∂H
∂t
=∇×E(1)
σE+ϵ∂E
∂t
=∇×H(2)
其中µ,ϵ,σ,σm分别是磁导率、介电常数、电导率、磁导系数，他们与位置有关。

旋度方程离散化后可以得到6个标量一阶微分方程。

对于光子晶体光纤仿真问题，假设场沿z轴传播，具有传播常数β。

即场分量具有形式ϕ(x,y,z)=ϕ(x,y)e iβz。

因此，对z的偏导可以用常数iβ替代，并且可以将3维的仿真问题转换为2维的FDTD仿真问题。

进一步，为了避免复数的出现，让E z,H x,H y具有cos(βz)成分，而E x,E y,H z具有sin(βz)成分。

于是，可以得到
图3:光子晶体光纤横截面介电常数分布
FDTD 的更新方程
H x |n +12
i,j +12=
1−σm ∆t 2µ1+σm ∆t 2µH x |n −12i,j +12−11+σm ∆t 2µ∆t µ∆y (E z |n i,j +1−E z |n i,j )+11+σm ∆t 2µβ∆t µE y |n i,j +12H y |n +12
i +12,j =
1−σm ∆t 2µ1+σm ∆t 2µH x |n −12i +12,j +11+σm ∆t 2µ∆t µ∆y (E z |n i +1,j −E z |n i,j )−11+σm ∆t 2µβ∆t µE y |n i +12,j H z |n +12i +12,j +1
2=1−σm ∆t 2µ1+σm ∆t 2µH z |n −12i,j +12−11+σm ∆t 2µ
∆t µ∆x (E y |n i +1,j +12−E y |n i,j +12)+11+σm ∆t 2µ
∆t µ∆y (E x |n i +12,j +1−E x |n i +12,j )E x |n +1i +12,j =1−σ∆t 2ϵ1+σ∆t 2ϵE x |n i +12,j +11+σ∆t 2ϵ∆t ϵ∆y
(H z |n +12i +12,j +12−H z |n +12
i +12,j −12)+11+σ∆t 2ϵβ∆t ϵH y |n +12i +12,j E y |n +1i,j +12=1−σ∆t 2ϵ1+σ∆t 2ϵE x |n i,j +12−11+σ∆t 2ϵ∆t ϵ∆y
(H z |n +12i +12,j +12−H z |n +12
i −12,j +12)−11+σ∆t 2ϵ
β∆t ϵH x |n +12i,j +12E z |n +1i,j =1−σ∆t 2ϵ1+σ∆t 2ϵE x |n i,j +12+11+σ∆t 2ϵ∆t ϵ∆x (H y |n +12i +12,j −H y |n +12
i −12,j )−11+σ∆t 2ϵβ∆t ϵ(H x |n +1
2i,j +12−H x |n +12i,j −12)(3)
其中Yee 元包中节点的选取如表1所示，电厂和磁场在时间和空间上交替采样。

本文仿真的光子晶体光纤是有缺陷的Holey 光纤，它的横截面介电常数分布如图3所示。

为了求解导模的场分布，可以通过初始化所有模场分量都有的一个初始场，通过时间的筛选，最终达到稳定时剩下的就是导模[5,6]。

从仿真过程来看，不是导模的成分就会被完美匹配层吸收掉，物理上来看，则是辐射出去了，只有导模才会被光子晶体光纤的周期结构束缚在光纤中。

表1:Yee元包场节点选取
t x y
H x n+1/2i j+1/2
H y n+1/2i+1/2j
H z n+1/2i+1/2j+1/2
E x n i+1/2j
E y n i+1/2j+1/2
E z n i j
横向场H⊥满足本征方程
∇⊥×(1
ϵr ∇⊥×H⊥)−1
ϵr
∇⊥(∇⊥·H⊥)=(µϵ0ω2−k
2
z
ϵr
)H⊥(4)
其中ϵr=ϵ/ϵ0。

横向本振场求解出来后，可以利用散度方程∇·H=0得到H z，电场可以利
用它与磁场的关系
E=
1
iωϵ
∇×H(5)
当得到所有的本征场后，初始场可以改写为
H=Σn(a n e iωn t+b n e−iωt)H n(6)
E=Σn(a n e iωn t−b n e−iωt)E n(7)为了方便计，事实上只要初始场各个分量都有即可，所以我在初始化时就直接初始化为高斯分布的场。

2.2完美匹配层(PML)
完美匹配层(PML)首先由Berenger提出，通过在FDTD仿真区域截断边界处设置一种特殊介质层，该层的波阻抗与相邻介质波阻抗完全匹配，并且对场有吸收，因而可以没有反射地吸收仿真区域外的场[4]。

完美匹配层的电磁参数需要满足阻抗匹配条件
σϵ=
σm
µ
(8)
同时还需要满足在分界面的两侧，横向电导率相同。

所以，可以得到PML的电导率(σx,σm x,σy,σm y)的分布图如图4所示。

而在实际仿真中，为了避免电导率突变所带来的反射，通常让电导率从分界面到最外侧逐渐递增。

3数值仿真结果
用高斯场做初始场，初始化Ez的幅度为2V/m，仿真区域有400×400个仿真节点，网格步长∆x=∆y=∆=a/40，其中a是晶格常数，圆孔的半径R=0.375a，中间缺陷大孔的半径r=1.5R，PML的厚度B=0.5a，电导率在PML中渐变分布为σ(t)=σ0(t−B)4。

程序流程图如图5经过1000次迭代收敛后的场分布图如图6所示。

图4:PML电导率分布
从仿真结果来看，光被束缚在纤芯中，没有泄漏出来的，验证了光子晶体光纤的导光原理。

同时，由于FDTD仿真过程可以显示出场的分布变化，所以也可以看到不同时刻场分布图如图7所示。

可以看到，随着时间的演进，场向四周辐射出去，碰到晶格后，部分被反射回去，还有一部分继续前行知道碰到PML被吸收掉。

当达到稳定后，场就不在想四周辐射了，场主要被约束在纤芯中，达到稳定的状态，这个时候可以认为场分布就是导模场的分布。

最后，为了验证自己仿真结果的正确性，将仿真结果与两篇文章中的结果[5,6]进行对比，发现虽然存在一些差异，但是总体结果是一致的，场都是稳定地集中分布在纤芯中。

比较的结果如图8所示，结果基本一致。

4总结
利用FDTD方法对周期结构进行仿真是一种非常有效的方法，并且结果可以非常直观的显示出来。

本文的仿真过程我已做成一个动画，放在附件中，通过动画可以非常直观地看到场随时间演化的过程。

最初选择光子晶体光纤这个题是因为在介质导波课程上，阅读了一篇相关的综述，感觉比较有意思，在选题的时候，便顺势选择了这个题目。

后来发现之前也有人做过光子晶体的仿真作为课程设计，为了避免重复，加上调研的时候发现有人做过导模仿真的问题，所以就选取了对光子晶体光纤的导模进行仿真。

在做课程设计之前，自己用Python写了一个简单的1维FDTD仿真的例子，作为练手，并做了一个视频在报告的时候展示给大家看了。

而在做仿真的过程中，发现这个问题比之前那个1维的问题要复杂得多。

首先是设计的方程个数从2个变成了6个，如果加上PML中z分量的分裂，一共有8个，这导致编程时由于不够仔细而出错的概率大大增加。

另外，在实现PML时，总是发现反射还是不小，也就是说吸收得不够彻底，有时还不收敛，经过仔细排查才发现是电导率选取的问题，由于没有采用指数差分，导致近似误差在电导率选取过大的时候
图5:仿真流程图
图6:光子晶体光纤导模的E z 的分布图
图7:不同时刻光子晶体光纤中场E z的分布图
图8:我的结果与论文中的结果的对比
非常大，甚至导致迭代过程中，场的旧值得权值变成负数。

后来改为指数差分，并适当减小电导率，并将电导率改为渐变的，吸收效果就好得多！
最后，感谢周老师还有和我一起讨论问题的同学以及实验室师兄师姐，让我在写程序的时候得到帮助和灵感！总的来说，通过这次课程设计，还是学到很多单纯通过课程教学学不到的东西，希望这门课程的课程设计以后可以继续下去，做的越来越好！
References
[1]J.D.Joannopoulos,S.G.Johnson,J.N.Winn,and R.D.Meade,Photonic crystals:molding
theflow of light.Princeton university press,2011.
[2]K.S.Yee et al.,“Numerical solution of initial boundary value problems involving maxwell’
s equations in isotropic media,”IEEE Trans.Antennas Propag,vol.14,no.3,pp.302–307, 1966.
[3]G.Mur,“Absorbing boundary conditions for thefinite-difference approximation of the time-
domain electromagnetic-field equations,”Electromagnetic Compatibility,IEEE Transactions on,no.4,pp.377–382,1981.
[4]J.-P.Berenger,“A perfectly matched layer for the absorption of electromagnetic waves,”
Journal of computational physics,vol.114,no.2,pp.185–200,1994.
[5]B.Cowan.(2003)Fdtd modeling of photonic crystalfibers.[Online].Available:
/grp/arb/tn/arbvol4/ARDB339.pdf
[6]M.Qiu,“Analysis of guided modes in photonic crystalfibers using thefinite-difference time-
domain method,”Microwave and Optical Technology Letters,vol.30,no.5,pp.327–330, 2001.。