介质波导的数值模拟

电磁场与电磁波报告 2002.12～2003.1混合 009 陈达如 3001144007 混合 003 胡睿 300300105911介质波导的数值模拟摘要：我们从麦氏方程出发，采用标量法对波导中电磁波满足的微分方程进行数值求解，并对其结果给出形象的图形表示。对于求解微分方程，我们采用有限差分法，将求解微分方程转化为等价的求解矩阵特征值和相应特征向量的问题，利用 MATLAB 在处理矩阵运算方面的强大功能，求解出电场在介质波导截面的分布。此外，我们对不同结构、不同折射率分布的波导进行了模拟分析，作为扩展给出了金属波导和光纤的近似数值模拟方法。本文深入浅出的介绍了这种介质波导数值模拟方法，学过电磁场理论并对 MATLAB 有一定了解的同学只要参考这篇文章便能很快学会这种方法。第一部分 原理一、电磁波方程的导出在无源简单介质中，麦克斯韦方程写成： HEj0从两个旋度方程消去 E 或 H 得到 2Ek其中 2设 ，取截面的二维矢量波方程为0zt其中kkzzt n2202即得如果我仅仅考虑电场的标量方程，则电场大小 E(Ex 或 Ey)满足：令 得电磁场与电磁波报告 2002.12～2003.1混合 009 陈达如 3001144007 混合 003 胡睿 300300105922----------------------------------------------------------(*) 这便是我们进行数值求解的微分方程。二、有限差分法有限差分法的基本思想是利用网格线将定解区域(场域)离散化为网格离散节点的集合，然后，基于差分原理的应用，以离散点上函数的差商来近似替代该点的偏导数，这样，待求的偏微分方程定解问题可转化为一组相应的差分方程组的问题。根据差分方程组(代数方程组) ，解出各个离散点上的待求函数值，即为所求定解问题的离散解，再应用插值方法便可从离散解得到定解问题在中个场域上的近似解。对于包括电磁场在内的物理场，应用有限差分法进行数值计算的步骤如下：(1) 采用一定的网格划分方式离散化场域；(2) 基于差分原理的应用，对场域内偏微分方程以及定解条件进行差分离散化处理(构造差分格式) ；(3) 对所建立的差分格式(即与偏微分方程定解对应的代数方程组) ，进行计算机求解。有关这方面的问题可参考偏微分方程数值解的书，下面结合我们的问题，介绍应用有限差分法求解方程(*) 。三、问题的求解我们考虑的介质波导结构如图 1 所示方型波导生长在 SiO2 衬底上，芯层折射率大于包层折射率(如图中所示 n1>n2) 。显然，如果芯层折射率比包层折射率大的多，电磁波将被限制在芯层中传播，在包层介质 n2 中，电磁波已经很弱，因此我们将包层介质与空气及衬底边界的电场设为零，在这样的假设下，求解介质波导截面电场分布就转化成下面的微分方程求解问题：由于波导形状规则，我们很容易将其作网格划分，网格线交点(节点)处的电场大小就是我们要求解的电场离散解(如图 2 所示) 。每一个节点的电场大小都是未知数(除了边界点外) ，我们求解的是(Nx-2) *(Ny-2 )个未知数 E(i,j) (i=2,3,4…Nx-1; j=2,3,4…Ny-1) ，电磁场与电磁波报告 2002.12～2003.1混合 009 陈达如 3001144007 混合 003 胡睿 300300105933下面我们从偏微分方程(*)中“提取”信息,构造求离散解所需的方程组。我们分析偏微分方程(*) 。首先将偏导数差分化，考虑函数 f(x),取小量△x=h,则同样的二阶微分,同样对 E(x,y) 的二阶偏导有,_代入偏微分方程(*)得：如图 2 所示，我们采用 E(i,j)表示节点处的电场，则有：如果记图中 A 点的电场为为 E(i,j) ，则上式给出了节点 B、C、D、E 处电场和节点 A 处电场的关系，即所谓的五点差分格式。对所有节点列出这种关系式，并将其写成矩阵的形式，得到：其中 X 是由各节点电场 E(1 ，1) ，E(1，2)…组成的Nx*Ny 个元素的列向量，A 是 Nx*Ny 行、Nx*Ny 列的矩阵，其每一行对应一个节点的五点差分格式方程。作为例子，我们给出如图 3 所示网格(节点处折射率均为 n)的矩阵方程：电磁场与电磁波报告 2002.12～2003.1混合 009 陈达如 3001144007 混合 003 胡睿 300300105944从简单的例子中可以看出矩阵 A 是个数字分布有规律的对称的庞大的稀疏矩阵，在我们的程序中给出了一种巧妙的方法来构造这个矩阵。至此，我们的问题就转化为求解矩阵 A 的特征值以及相应的特征向量，从电磁波理论上讲，这里的一个特征向量对应一种电磁场在波导中的模式。而求解求矩阵 A 的特征值以及相应的特征向量，在 MATLAB 中，只要用一个语句就行了(用 eigs()函数或者 eig()) 。求出特征向量后，进行归一化，就得到离散节点处电场相对强弱的分布，再用 MATLAB中的作图函数直观地显示出来。值得注意的一点是我们似乎没有对芯层折射率进行考虑，其实不然，在稀疏矩阵 A 中已经包含每个节点对应的折射率 n(i,j)的信息。正是基于这一点，我们可以考虑多种介质复杂结构的波导的数值模拟，在程序中使用者可以自由地设计波导结构和各个部分的折射率。当然，由于这种模拟方法的局限性，不能对超出矩形和其组合的波导，如三角形，圆形进行数值模拟。自由是相对的。四、金属波导和光纤的模拟我们这里不是从金属波导或光纤中电磁波所满足的电磁波方程来给出其数值模拟，而是仅仅从改变上述介质波导的结构和边界条件，仿造金属波导和光纤的结构，然后给出其近似的数值模拟结果。在金属波导中我们考虑先求出 x 方向、y 方向电场分量，再将其合成截面的总电场。x 方向、y 方向电场分量的求法是相似的，这里是讨论 Ey 的求解。我们仍然采用图 2 作说明。对金属波导来说，波导折射率为 1，因此将图 2 中芯层和包层的折射率均设为 1；对求解来说，边界条件是在 x=a,x=b 处的电场为零，在 y=c,y=b 处应该没有限制，我们采用近似的处理方法，将 y=b 趋向正无穷远处，y=c 趋向负无穷远处(实际上给出 y 方向充分大的线度即可，我们在程序中默认用 3 倍于中心介质层模拟) ，然后将 y=c, y=b 的边界置为零就不会影响结果了。然而，这种方法只能应用到 10 和 01 模，对于多模情况无效，这是我们通过对比解析解得出的结论，其原因有待进一步的研究。对于光纤来说，由于光场集中分布在纤芯中，在离纤芯较远的边界设将光场设为零的条件是可以接受的。如图 4 所示，芯层、包层、外层的折射率递减，外层的边界电场设为零。根据折射率设置的不同，我们可以对阶跃光纤和渐变光纤进行数值模拟。电磁场与电磁波报告 2002.12～2003.1混合 009 陈达如 3001144007 混合 003 胡睿 300300105955第二部分 程序一、概述在这一部分里我们将结合附