异形脊波导传输特性分析

题 目 异形脊波导传输特性分析 学生姓名 杨璐 学号 1213014118

所在学院 物理与电信工程学院 专业班级 电子1204班 指导教师 帅春江 完成地点 陕西理工学院

2016年 6月 5日

陕西理工学院毕业设计

异形脊波导传输特性分析

作者：杨璐

（陕西理工学院 ,物理与电信工程学院, 电子信息工程系2012级4班，陕西 汉中 723000）

指导老师：帅春江

[摘要]在现代微波工程中，为了满足微波传输系统性能的某些需求，需要不断探索和研究具有特殊截面形状的

各类新型波导，多脊波导的理论研究对解决实际问题有很好的参考价值。由于其能实现微波小型化及具有良好的宽带特性等优点，具有很好的工程应用价值。有限差分法突破了解析法只能分析结构简单的波导器件，并且应用其处理脊形波导结构具有不需要复杂的物理模型等许多优点。本文从泰勒级数出发推导出用有限差分法计算脊波导的截止频率，以单脊波导为例，编程计算及验证了单脊波导内TM波的截止频率，并分析计算了新型脊波导-梯形脊波导的TM波的截止频率，所得结果为多脊波导的设计与制作提供理论参考。

[关键词]有限差分；单脊波导；梯形脊波导；截止频率；TM波

陕西理工学院毕业设计

The analysis of transmission characteristics about the

shaped trapezium ridge waveguide

Author： LuYang

（Grade 11, Class 1, Major electronics and information engineering, School of Physics and Electronic

Information Engineering, Shaanxi University of Technology, Hanzhong 723000, Shaanxi）

Tutor： Chunjiang Shuai

Abstract:In modern microwave engineering, in order to meet some of the requirements of microwave transmission

system performance, we need to continue to explore and research new types of waveguide having a special cross-sectional shape, more ridge waveguide theory to solve practical problems for a good reference value. Due to its size and to achieve microwave has good broadband characteristics, etc., and has good engineering application value. Finite difference method to break the analytical method can only analyze a simple structure waveguide devices, and applications which deal with ridge waveguide structure has many advantages not require complex physical models. From the Taylor series starting deduced ridge waveguide by finite difference method cutoff frequency to a single ridge waveguide, for example, calculation and verification program within a single ridge waveguide cutoff frequency field intensity distribution and TM waves, and analysis intensity distribution shaped trapezium ridge waveguide TM wave and the cut-off frequency, the results provide a theoretical reference for the design and production of multi-ridge waveguide.

Key words: Finite-Difference;single ridge waveguide; trapezium ridge waveguide;Cutoff frequency;TM wave

陕西理工学院毕业设计

目录

1 引言 ........................................................................................................................................................ 1

1.1 脊波导的发展 ........................................................................................................................ 1 1.2 课题意义 ................................................................................................................................... 2 1.3 数值法简介.............................................................................................................................. 2

2 有限差分基本原理 ...................................................................................................................... 3

2.1 基本差分公式 ........................................................................................................................ 3 2.2差分方程的求解过程 .......................................................................................................... 4 2.3 波导中的电磁场方程 ........................................................................................................ 6 2.4 亥姆霍兹方程的差分表达式 ........................................................................................ 8

3 FORTRAN语言 ............................................................................................................................ 12

3.1 简介 ........................................................................................................................................... 12 3.2 基本程序结构 ...................................................................................................................... 12

4 脊波导TM波传输特性分析 ............................................................................................. 13

4.1 单脊波导传输特性分析 ................................................................................................. 13 4.2 梯形脊波导传输特性分析 ............................................................................................ 15

5 总结 ...................................................................................................................................................... 23 致谢 ............................................................................................................................................................ 24 参考文献 ................................................................................................................................................. 20 附录 A 外文文献 ......................................................................................................................... 21 附录 B 中文翻译 ......................................................................................................................... 30 附录 C 程序设计 ......................................................................................................................... 29

陕西理工学院毕业设计

1 引言

电磁波在有界空间的传播，就被称为导波系统中的电磁波。所谓导波系统是指引导电磁波沿一定方向传播的装置，它不仅是能量和信息的载体和传递工具，而且是构成各种高频、微波元件和仪

[1]

器的基础。脊波导是其中的一种常见的导波结构，它是矩形波导的一种变形，又称凸缘波导。

近年来，微波因其本身的特点，在信息获取、传输中所占地位愈显重要，在这样的大环境下，各种高质量的微波器件（包括用脊波导制作的各种器件）会有很大的实际需求，在雷达及卫星通讯设备中，最常用的是矩形波导和圆波导，在现代微波工程中矩形脊波导具有较低的截止频率,较宽的

[2]

工作带宽,低特性阻抗等优点,使得矩形脊波导在微波和毫米波器件中被广泛应用，如：宽带脊波导滤波器、宽带定向耦合器、双工器、变频器、移相器、混频器、低噪声放大器、脊波导缝隙天线阵等。

尽管脊波导的应用越来越广泛,然而有关它的理论分析计算的资料却较少，因此如何计算脊形波导中传输的电磁波的截止频率是一个值得研究讨论的问题。经过对论文的研究发现，对非规则脊波导研究的文献很少，在本文中将采用有限差分法这一经典的数值方法结合计算软件Visual fortran来分析梯形脊波导的传输特性。 1.1 脊波导的发展

早在1947年，Cohn Seymour. B在论文“Properties of ridged waveguide” 中对脊波导特性进行了初步研究，利用横向谐振技术给出了求解矩形脊波导的截止波长和特性阻抗的方程和特性图，从此脊波导便一直受到广泛的关注。

1955年，Hopfer 在论文“The Design of Ridged waveguides”中给出了一些不同外观比例的脊波导各个参数的计算方法，还给出了一些设计曲线。到了60年代，Pyle用准静态法得到了任意外形比例的脊波导的特征值，其结果比Hopfer的更精确。以上文章对脊波导的研究属于开创性的，但都只能解决脊波导的TEn0模的特征值问题，都对其它高次模的解无能为力，特别是当脊很薄和脊沟很窄时，其精度受很大限制。70年代初Konishi等把用脊波导制作的12GHz低噪声变频器等器件成功地用于卫星通信，并对矩形脊波导的不连续性作了分析计算，使得矩形脊波导在微波和毫米波器件中被广泛应用。80年代，Monsour用模式分析技术得到了双脊波导的TE模的近似闭式表达式，其中用到了横向谐振的概念。90年代，J.Helszajn用有限元法分析了脊波导的特性阻抗，给出了精确的曲线图，并提出了梯形脊波导，并分析了其阻抗特性。2000年Javis给出了任意高度的双脊波导设计

[3]

曲线，他们所用的方法都属于与上述文章类似的等效电路法。

国内在脊波导的研究领域起步较晚，但发展迅速。1997年黄彩华在文献“矩形变形脊波导主模截止波长和特性计算”中，用等效电路法得到了矩形变形脊波导的主模和特性阻抗。1999年李永明等在论文“有限元外推法在波导本征值问题中的应用”中用有限元法对矩形双脊波导进行了研究，得到了

[4]

主模截止波长相关数据。2004年任列辉等在“用电磁场算子理论分析脊波导的传输特性”中，用电磁场算子理论求解脊波导的本征值，在此基础上讨论了脊波导的的传输特性。王萍在“脊波导各种参数的计算”一文中，用等效电路法对脊波导进行了分析，得到了截止波长，特性阻抗，功率容量等数据。2006年薛德强在论文“下限单脊矩形波导传输特性研究”中，用有限差分法对下限单脊矩形波导进行了研究。2007年陈小强教授等在论文“两种新型双脊波导传输特性的研究”中提出了三角形和倒

[4]

梯形两种新型对称双脊波导，并计算了TE模式的传输特性。

尽管波导历史悠久,但时至今日仍有各种新的导波结构不断涌现。如今种类繁多的波导广泛应用于现代通讯,雷达及微波,毫米波,光电系统中。其种类包括了传输线类波导;金属类空心波导(如圆波导,矩形波导,椭圆波导,三角形波导及近年来新兴的五边形波导等);微带线类波导(如屏蔽微带线,带状线

[5]

等);介质类波导(如阶跃型光纤,渐变折射率型光纤,椭圆光纤,蝴蝶结光纤,光纤光栅等)。各种新材料波导更是层出不穷,等离子体金属波导,等离子体介质波导,非线性波导,光子晶体光纤等成为波导研究的重要领域。纵观今日的各类波导,工作频率已从射频贯穿到光频,波导横截面尺寸从米量级直到亚微米级。

第 1 页 共 35 页

陕西理工学院毕业设计

1.2 课题意义

在现代微波工程中脊波导具有较低的截止频率，较宽的工作带宽，低特性阻抗等优点，使得脊波导被广泛应用。它常用作宽频带传输系统及宽频带测试系统，比如宽带脊波导滤波器、宽带定向耦合器、双工器等等；也常作为微波管的输出波导，比如：变频器、移相器、脊波导缝隙天线阵等等。为了满足微波传输系统性能的某些需求，我们还需要不断探索和研究具有特殊截面形状的各种新型波导。最近几十年来各种结构形状的脊波导应运而生，它们可以看作是矩形波导的变形，本文所研究的梯形脊波导就是其中一种情形。

长期以来分析波导问题的解析方法和数值方法有很多，有限差分法和有限元法对复杂问题具有较强的适应性，近年来,为了解决大规模电磁场数值计算问题，基于有限差分和有限元的区域分解法已受到人们的广泛重视。

本文就是采用有限差分法对梯形脊波导的传输特性进行了分析，得到了它们的截止频率的一些变化规律，给脊波导在今后工程应用中提供了一定的数据参考价值。 1.3 数值法的简介

所使用的解析法和数值法都是基于实际问题的解析模型进行求解，其过程是从基本的麦克斯韦方程组或波动方程出发，加上特定的边界条件和本构参数，推导出微分方程或积分方程的体系，然

[6]

后尽可能地采用解析手段具体问题具体分析，最后则可以编制出一个紧凑而且高效的计算程序。

数值法是指直接将待求解的数学方程进行离散化处理，将无限维的连续性问题化为有限维的离散问题，将解析方程的求解问题化为代数方程的计算问题的一类方法。用高性能的计算机，可以直接以数值的、程序的形式代替解析形式来描述电磁场问题。在纯数值法中，通常以差分代替微分、用有限求和代替积分，将问题化为求解差分方程或代数方程问题。其优点是原则上它可适用于任何复杂的边界问题的求解，既可用于求解第一、二类可解的问题，又可用于求解因边界复杂由它们无法解决的问题，原则上可以求得所需要的任意精度，仅受限于计算机容量、速度和舍入误差、但随着大容量、高速度巨型电子计算机的发展和数值计算方法本身的发展，可以预料，许多过去只能靠

[7]

实验来分析和设计的电磁场边值问题都可以靠数值法来实现，因而数值法有着广阔的应用前景。

数值法，包括有限差分法、矩量法、边界元法及有限元法等。矩量法及边界元法的特点是计算精度高，但对于复杂边界形状不易导出特定的积分，不能模拟任意的媒质分布参数问题。有限元法适应性强，可以用于各种复杂的媒质分布等问题。有限差分法由于其直接从麦克斯韦方程组出发，不需要任何导出方程，这就避免了使用更多的数学工具，使其成为所有电磁场计算方法中最简单的一种。还有一种数值方法就是时域有限差分法,它是在有限差分法的基础上提出的，也是当今重要的电磁场数值计算方法之一，但由于其通常对计算机的内存配置和CPU的运算能力要求较高，其应用受到一定限制。

第 2 页 共 35 页

陕西理工学院毕业设计

2 有限差分基本原理

2.1 基本差分公式

所我们要求的电位函数u，它在区域D内满足下面的拉普拉斯方程

2u2u20 2xy (2-1)

在边界上S，它服从以下条件：

uSfp (2-2)

式中fp为边界点p的函数。这类问题一般称为第一类边值问题或称狄里赫利问题。

为了用差分方法求解电位分布，先在xy平面分别作两族平行于x轴和y轴的直线,线间的距离为h，于是各直线的x和y坐标分别为：

xiih; yjjh

式中i,j为正整数，取值1、2、……。这样区域D就被许多边长为h的正方形所覆盖，在图2-1中示出了这种情况。各正方形的顶点被称为网格的节点，从图可以看到，各节点所处位置有所不同。一些节点（例如a节点）恰落在边界上S，我们把它叫做边界节点。有些节点到边界的距离不足h（例如节点b）,这些节点叫做不规则节点。但是大部分节点到边界的距离大于h，例如图上的0点，它们属于规则节点。差分法就是求这些离散节点处u的近似值。以后将用符号(i,j)来表示节点(xi,yj)，而将节点处(i,j)的值u用uij代表[8]。

y h

j=1,2,3…

2 D b （x0，y0） 0 a 3 h3 h4 h2 0 h1 1 h x

i=1,2,3… 图2-1 将场域划分成网格

4 图2-2 不等距网格

当用差商来替代偏导数，可以从(2-1)式导出它的差分表达式。在直角坐标系如果采用图2-2所示的不等距网格时，求出的差分表达式如下

1u1u2u3u41u00 (2-3) h1h1h3h2h2h4h3h1h3h4h4h2h1h3h2h4 对于常用的等距网格情况，h1h2h3h4h。根据式（2-3）就得到

u1u2u3u44u00 (2-4)

在轴对称情况，拉普拉斯方程就成为下面的形式

第 3 页 共 35 页

陕西理工学院毕业设计

2u1u2u20 (2-5) 2rrzr它的差分表达式对不等矩网格为

r ２

h2 ３ h3 h1

１

r0 轴线 ４

Z

图2-3 轴对称情况的不等距网格

2r0h422u1u3u4h1h1h3h3h1h3r0h2h2h422r0h22r0h4h2u4hhr0h4h2h4h2h4r013u0 (2-6)

式中各量见图2-3。

在等距网格情况，h1h2h3h4h。如果令轴线处j1，而0点落j行并且j1，则r0j1 h。根据（2-6）式就得出

11u1u31 u122j1 u44 u00 (2-7)

2j1假使0点落在j1的轴上，可以求出差分格式为

u1u34 u26 u00 (2-8)

2.2差分方程的求解过程

下面用一个计算金属槽内电位分布的例子来说明差分方程的求解过程。设金属槽上盖板电位等

于100V，侧壁与底壁的电位为0V。将场域用正方形网格划分，要求的是场域内部九个节点的电位值。对网格节点1，根据式（2-4）的等距差分公式可以得到

100 100

100

4u1u2u41000

也可以扩展成以下的形式

0 0 0

0 0 0 0

0

0

4u1u20u3u40u50u60u70u80u9100

对其它标号为2~9的各个节点也可写出类似的扩展方程，将它们组合起来，就得到下面的依矩阵形式表示的线性代数方程组

图2-4 金属槽结构

第 4 页 共 35 页

陕西理工学院毕业设计

0100000u110041141010000u210014001000u31000u10041010004010141010 u50 u00101400106000100410u70u00001014108000010140u90这是一个包含九个变量u1~u9的方程组，可采用线性代数中塞德尔迭代法求解。但是采用这种方法，

必须先给各个变量赋以初值，这样就可从第一至第九个方程依次求出u1~u9的改善值。不断重复种计算，直至相邻两次算出的各个u值收敛到我们需要的精度为止。

在实际的计算过程中，当然用不着先列出上面所述的方程组，它仅在网格节点很少时候才有可能。注意到前面列出的各个差分公式中，网格中某节点的电位只和它四邻节点电位有关这一特点，例如图1-6中节点5的电位仅和它四周相邻点2、4、6、8的电位有关。就可以从左至右，从上至下

[9]

一行一行地扫描逐一地求出各网格节点的电位。

为判别计算的收敛程度，一般采用以下两种方法：

（一） 对所有网格点，找出相继的两次迭代的最大相对误差

erukk1uiimax1k (2-9)

uii式中k表示第k次迭代。

（二） 找出相继两次迭代形成的各网格节点的平均相对误差erc

uuercukijki,jkiji,jk1ij (2-10)

对同一问题，由第二种方法求出的收敛误差比第一种方法小1~2个数量级，所以如果用第一种

56[6]

方法是收敛精度误差控制在104，则用第二种方法时就应控制在10~10。

根据以上所说，可以写出用差分方法计算静电电位问题的框图如图2-5。

第 5 页 共 35 页

陕西理工学院毕业设计

启动 给定边界给定边界电位和给其它节点赋初值 位置信息 迭代次数k＝0 k＝k +1 j＝1 i＝1 规则点乎？ 是 用规则点的五点 差分公式计算 否 用不规则点差分公式计算 第j行扫描结束乎？ 否 是 否 j＝j +1 全区扫描结束乎？ i＝i +1 是 否 计算收敛到指定精度乎？ 是 打印计算结果 图2－5 差分法计算框图

第 6 页 共 35 页

陕西理工学院毕业设计

2.3 波导中的电磁场方程

当波在任意截面的波导中传播时，可以将其中传播的电磁波分成TE波和TM波别把各个场分量关系总括如下：

对TM波，如令Ezu，则其它各个场分量为：

[10]

。现在将分

Hz0

ExjEz2，EyjEz2 (2-11)

kkpxkkpyHjEzxk2y，EjEzk kpx2kpx以上诸式中的横向电磁场分量可以归并成以下简洁的形式：

EiEzEztjk2kpxxiyy jk2kptEzjk2kptuHjtk2iEzEzkpxyiyxjjk2iztEzk2iztukpkp对TE波，可令Hzu，则其它各个场分量为

Ez0

ExjHzjHzk2，Ek2 kpyykpxHjHzjxk2，HHzkpxyk2 kpy各个横向分量也可合并如下：

HjHzjHztixk2ikpxyk2kpyjjk2tHzkpk2kptuEHztjk2kpiHzxyiyx jj k2ikpztHzk2ikpztu在上面各式中，为相位常数，它等于

第 7 页 共 35 页

(2-12) (2-13) (2-14) (2-15) (2-16) (2-17) (2-18)

陕西理工学院毕业设计

2k2kp (2-19)

kkp被称为截止波数，他由下式决定

kkp2fkp/vc (2-20)

式中，fkp——波导的截止频率； vc——光速。

2.4 亥姆霍兹方程的差分表达式

无论是TE波还是TM波，上面定义的u函数均满足亥姆霍兹方程

2u2u22kkpu0 (2-21) 2xy但是对于不同波型有不同的边界条件。对TM波，ul0；对TE 波,面的边界线。

用求拉普拉斯方程的差分方程一样的方法，不难求出在不等距网格情况下，亥姆霍兹方程的差

分表达式为

u0，l代表波导截nlu1u31u2u41222u0kkphu00 (2-22) pprqqsrprsqsprqs式中的p、q、r、s定义见图2-6。

对等距网格情况，pqrs1，就得到

22u1u2u3u44kkphu00 (2-23)

22当用式（2-21）或（2-22）求解亥姆霍兹方程时，和解拉普拉斯方程有所不同。前者需要给定kkph值，但这个值在计算前是不知道的。为消除这个困难，可采用下面的措施，就是在计算开始前，除

22了给内部网格节点赋初值外，也给kkph一个初值，有了这个初值，就可进行u值迭代计算而将各节

点u值改善到一个较好的精度，如果松弛系数选得适当，只要几次迭代就行了。然后用各节点的这

2222些u值用下面讲的方法求出改进的kkph值。再利用这个kkph值去求u，这种过程循环重复，直到22kkph和各节点u值均收敛到应用的精度为止[11]。

为了从各网格节点u值求kkph，可采用以下方法。用u乘式（2-20）后并将整个等式对波导横截面S进行面积分，得到：

222u2u22udSkudS0 kp22ssyx或者写成

第 8 页 共 35 页

陕西理工学院毕业设计

uudSkss2t22kpudS0

2式中t为对横向坐标的拉普拉斯算子。由上式解出kkp

2kkput2udSsus2dS (2-24)

利用这个公式求kkp，直观上也可看出比直接用式（2-21）来求更合理，因为它意味着在整个场域内的平均，这就可能使计算形成的误差部分地被抵消。

（2-23）式可用下面的差分近似式来表示

22kkph2uijui,j1ui,j1ui1,jui1,j4uijSiji,j2uijSiji,j (2-25)

上式中的代表对所有网格节点求和，Sij代表第（i，j）节点所占有的面积则点，例如图2-7中的g点，Sijh。但对于面积上的W点，Sij

2 3 2[15]

。对于区域内的规

12h； 2z

rh sh qh 0 p,q,r,s1 1 ph 4

g

w

图2-6 不等距网格

对于边界上的角点z，Sij12h，这些关系由图上的阴影区容易看出。在TM波情况，由于边界4图2-7 求边界节点的Sij

2上uEz0，所以沿边界的乘积uijSij0，这是在式（2-24）中所有其他内部点Sij均等于h而可约去，式（2-24）就简化成

khkp2uijui,j1ui,j1ui1,jui1,j4uiji,juij2i,j (2-26)

在以后对波导计算时，编制的计算机程序所采用的方法同于计算电位的超松弛法，但稍作改变如下：

将式（2-22）改写成下面的形式

22u1u2u3u4kkphu0,old4u0,oldR0 (2-27)

式中R0叫做余量，u0,old表示0点的第k1次迭代求出的u值，各点1，2，…等位置见图2-2。再

第 9 页 共 35 页

陕西理工学院毕业设计

,new，它应使得式（2-26）的余量为零，即 来求0点的新u值u022,new0 (2-28) u1u2u3u44kkphu0为了加速收敛，乃令第k次迭代求出的u0值为u0,new，它等于

,newu0,old (2-29) u0,newu0,oldu0式中为松弛系数。用式（2-26）减去式（2-25）可得

4khu22kp224k0,newkphu0,oldR0

或者

,newu0,oldu0将式（2-27）代入式（2-28）就有

R0 (2-30) 224kkphu0,newu0,old4kR02kph2 (2-31)

22计算过程是在给定了kkph后，先由（2-29）式计算各节点余量R0，再根据式（2-30）求出改

善后的u0,new，用这种方法算出场区各点的u值。重复这样的计算，经对指定的迭代次数后，再由式

222222（2-24）或（2-25）求出较精确的kkph值。利用这个新的kkph值，重复上述过程直至kkph值收敛

为止

[16]

。

22在求出kkph后，就可以算出截至频率fkp。因为

2fkp222kkphvh

c由上面的公式可解出

2fkp1vc2hkhkp2 (2-32)

根据以上所述，整个用差分法求解波导问题的计算框图示如图2-8。

第 10 页 共 35 页

陕西理工学院毕业设计

给边界点赋值和给内部节点赋初值 输入khkp2的初始估计值 进行指定次数的迭代以改善场量u的值 计算新的khkp2值 khkp2值是否收敛 否 总迭代次数大于指定次数乎？ 否

是 给边界点和内部点赋值 给边界点赋值 是 图2-8求解波导问题的计算框图

第 11 页 共 35 页

陕西理工学院毕业设计

3 FORTRAN语言

FORTRAN语言，是第一个面向过程的高级程序设计语言，主要用于科学计算，也可用于数据处理和仿真。FORTRAN是英文FORmula TRANslator的缩写，原意是公式翻译。FORTRAN语言可使程序员用一种非常接近于常用数学表达式和英语自然语言的方式编制计算机程序。自1956年开始使用以来，一直在国际上广泛流行，是使用最广泛的程序设计语言之一。 3.1 简介

FORTRAN语言问世以来，根据需要几经发展，先后推出了不同的版本，其中最流行的是1958年出现的FORTRANⅡ和1962年出现的FORTRANⅣ。1966年美国标准化协会(ANSI)公布了两个美国标准文本：

·标准FORTRAN（X3.9-1966），大致相当于FORTRANⅣ。 ·标准基本FORTRAN（X3.10-1966），大致相当FORTRANⅡ。 1972年国际标准化组织(ISO)接受了美国标准，在稍加修改后公布了ISO FORTRAN标准，即《程序设计语言FORTRAN ISO 1539-1972》，它分为三级，即：

· 完全的(一级)FORTRAN，相当于FORTRANⅣ。 · 中间的(二级)FORTRAN，介于FORTRANⅡ和FORTRANⅣ之间。 · 基本的(三级)FORTRAN，相当于FORTRANⅡ。

FORTRANⅣ(即FORTRAN 66)流行了十几年，几乎统治了所有的数值计算领域，许多应用程序和程序库都是用FORTRANⅣ语言编写的。美国标准化协会(ANSN)在1976年对ANSI FORTRAN(X3.9-1966)进行了修订，预定在1977年通过，为了区别于FORTRAN 66，新标准定名为FORTRAN 77。实际上到1978年4月才由ANSI正式公布作为新的美国国家标准。即FORTRAN(X3.9-1978)。1980年，FORTRAN 77被接受为国际标准，即《程序设计语言FORTRAN ISO 1539-1980》，该标准分为全集和子集。

中国制订的FORTRAN标准，基本上采用了国际标准，于1983年5月公布执行，标准号为GB3057-82。

FORTRAN 77标准完成后，新版本的修订工作也在同一时间开始进行。这个版本进行了15年，最后在1992年正式由国际标准组织ISO公布，它就是FORTRAN 90。FORTRAN 90对以往的FORTRAN语言标准作了大量的改动，使之成为一种功能强大、具有现代语言特征的计算机语言。其主要特色是加入了面向对象的概念及工具、提供了指针、加强了数组的功能、改良了旧式FORTRAN语法中的编写“版面”格式。

FORTRAN 95标准在1997年同样由ISO公布，它可以视为是FORTRAN 90的修正版，主要加强了FORTRAN在并行运算方面的支持。同时一些公司纷纷推出Visual Fortran，这为工程技术界进行科学计算和编写面向对象的工程实用软件的用户提供了极大的方便。熟悉VB或VC的读者可以很容易地掌握Visual Fortran的使用，进一步开发出自己专业领域的Windows下的界面友好的工程应用

[13]

软件。 3.2 基本程序结构

FORTRAN程序由一个或几个相对独立的程序段组成，其中必须有一个主程序段(PROGRAM),还可以有（也可以没有）子程序段（subroutine）、函数段(FUNCTION)或数据段(BLOCK DATA)。

[14]

程序段由语句组成。语句分为可执行语句和非执行语句两类。可执行语句是在程序执行时能导致系统硬件作出一个实际操作的语句，包括赋值语句、控制语句、输入输出语句等。非执行语句是在程序编译时为编译程序提供有关信息的语句，包括说明语句、格式语句、数值语句、参数语句、函数定义语句和程序段语句等。语句是由常数、变量、运算符和专用定义符等按事先规定的格式书写。语句标号则提供了引用单个语句的方法。

第 12 页 共 35 页

陕西理工学院毕业设计

4 脊波导TM波传输特性分析

4.1 对称单脊波导传输特性分析

对于脊形波导中传播的TM波，Ez0,Hz0；在边界应满足Ez0的条件。此外，就基波而言，它的Ez分量称于截面中心线，所以只要研究截面的一半区域就可以了，此时中心处的边界条

2Ez[15]

0，n为垂直中心线的外法线方向。由于沿金属边界Ez0，故而计算kkph时可以n应用式（2-25）。

件为

根据图2-10框图编制的程序列在下面(见附录C)，其中ALFA变量代表松弛系数，SQKHOLD代表上一计算得到的kkph值，SQKHNEW代表新算出的kkph值，Ez为计算内部节点时常用的一个函数子程序。如图4-1为对称单脊波导的立体模型，图4-2所示为对称单脊波导的横截面图。

[16]

22

图

4-1 对称单脊波导立体模型

图4-2 对称单脊波导截面图

对于对称单脊波导，由于其对称性，我们仅仅需要研究其一半电位分布来求取截止频率。

[17，18]

待计算的对称单脊波导的一半截面示如图4-3，采用正方形网格对波导进行划分。HH为波导侧壁高度。

c a HH b d e 图4-3 对称单脊波导的网格划分图

第 13 页 共 35 页

陕西理工学院毕业设计

4.2对称单脊波导的TM波的截止频率分析

由于工程上的所需，改变脊的高度与长度，波导的TM波的截止频率会有不同的变化。如图4-3所示，取HH=5.0mm，a=2.5mm,b=2.5mm,c=2.5mm,d=1.25mm,e=1.25mm。划分网络c=10；a=5；b=5-10；d=5；e=5-10。

计算所得截止频率数与表4-1理论近似，计算结果正确，程序无误。

表4-1计算值与参考值比较f（GHz）

[19]

本文 参考值

[20]3.2226

3.2223

(1) 取d=5不变时， 改变b/HH的比值所得到的TM波的截止频率如表4-2所示，据此绘制出如图4-4所示的对称单脊波导TM波截止频率随脊的高度比变化曲线。

表4-2对称单脊波导的TM波随高度比变化的截止频率（GHz）

b/HH

0.1

0.2

0.3

0.4

0.5

0.6

0.7

f 20.184 22.389 25.307 28.990 32.226 33.487 33.806

f/*104HZ 截止频率曲线变化图 3.53

图4-4 对称单脊波导TM波截止频率随脊的高度比变化曲线图

20.10.20.30.40.5 b/HH) 0.60.70.82.5 由图4-5可看出，随脊高比（b/HH）的增大，TM波的截止频率呈上升趋势，且比例为0.5时开始趋势较为平缓。

第 14 页 共 35 页

陕西理工学院毕业设计

(2) 当a=5不变时，改变d/c的比值得到的TM波的截止频率如表4-3所示，据此绘制出如图4-5所示的对称单脊波导TM波截止频率随脊的宽度比变化曲线。

表4-3对称单脊波导的TM波随宽度比变化的截止频率（GHz）

d/c

0.1 24.430

0.2 26.331

0.3 28.847

0.4 32.226

0.5 35.946

0.6 37.275

0.7 37.460

f

4截止频率曲线变化图 3.5 f/*104HZ 32.50.10.20.30.4 d/c 0.50.60.70.8 图4-5对称单脊波导TM波截止频率随脊的宽度比变化曲线

由图4-7可以看出，随脊的宽度比（d/c）的增大，TM波的截止频率呈上升趋势，比列为0.5时趋势比较平缓。

4.2.1梯形脊波导传输特性分析

图4-6梯形脊波导的实际模型 图4-7 梯形脊波导的剖面图

4.2.2梯形脊波导的TM波的截止频率分析

由于工程上的所需，改变脊的高度比，波导的TM波的截止频率会有不同的变化。如图4-8所

第 15 页 共 35 页

陕西理工学院毕业设计

示HH=10.0mm，a=5.0mm,b=5.0mm,c=28.0mm,d=5.0mm,e=8.0mm，θ=45。划分网络c=20；

。

a=10；b=10-20；d=10；e=8。

保持d不变，改变b/hh,得到的TM波的截止频率如表4-4所示，据此绘制出如图4-8所示的梯形脊波导TM波截止频率随脊的角度变化曲线。

表4-4对称梯形脊波导的TM波随高度比变化的截止频率（GHZ）

d/c

0.1 5.94

0.2 3.37

0.3 1.98

0.4 0.947

0.5 0.808

0.6 0.741

0.7 0.695

f

图4-8梯形脊波导TM波截止频率随脊的高度比变化曲线

由图4-8可以看出，随脊高比（b/HH）的增大，TM波的截止频率呈下降趋势，且比例为0.5时开始趋势较为平缓。

保持b不变，改变d/c得到的TM波的截止频率如表4-5所示，据此绘制出如图4-9所示的梯形脊波导TM波截止频率随脊的角度变化曲线。

第 16 页 共 35 页

陕西理工学院毕业设计

表4-5对称梯形脊波导的TM波随高度比变化的截止频率（GHz）

b/HH

0.1 18.56

0.2 15.24

0.3 12.93

0.4 9.47

0.5 0.6 0.7

f 8.36 8.29 8.24

图4-9梯形脊波导TM波截止频率随脊的宽度比变化曲线

由图4-9可以看出，随脊高比（b/HH）的增大，TM波的截止频率呈下降趋势，且比例为0.5时开始趋势较为平缓。

第 17 页 共 35 页

陕西理工学院毕业设计

5 总结

本文主要运用有限差分法（FDM）分析了脊波导的TM波的传输特性，用FORTRAN程序语言编制系列程序，对对称单脊与梯形脊波导中电磁场进行模拟，从而得到在这些波导中传输的电磁波的截止特性。主要做了以下工作：

（1）对有限差分法的基本概念及原理进行了综合性论述，通过求解波导截止频率给出其相关公式推导；

（2）用泰勒级数法推出直角坐标系下拉普拉斯方程的差分格式，并从直角坐标系中TE波电磁场纵横关系出发，以对称单脊波导为例，采用FDM算法，用FORTRAN程序语言编制系列程序分析计算其截止频率，并对梯形脊波导进行模拟，从而得到在梯形波导中传输的电磁波的截止特性。 （3）通过对各种脊形波导的结构变化分析，绘制出其TM截止频率随高度比与长度比变化曲线， 得到其传输特性。

第 18 页 共 35 页

陕西理工学院毕业设计

致谢

本论文的写作从开题、需求调研、搜集资料、分析设计到最后成文，历时数月。其间，得到了老师和同学及朋友的各种帮助。在此，我表达衷心地感谢。感谢导师帅春江一直以来的精心指导,从立题到完成学位论文,每一步都凝聚着他的心血，老师以其严谨求实的治学态度、高度的敬业精神、孜孜以求的工作作风和大胆创新的进取精神对我产生了重要影响。他渊博的知识、开阔的视野和敏锐的思维给了我深深的启迪，将我带入一个崭新的领域,导师对学术的严谨态度和执著追求的精神,将永远激励我努力奋斗。导师宽以待人和豁达的生活作风,将使我终生受益。深深地感谢帅老师为我创造的一切学习的条件。

感谢物电学院领导的支持和鼓励,感谢电子信息工程实验室提供的条件,感谢电子信息工程系的同学的关心和支持。

第 19 页 共 35 页

陕西理工学院毕业设计

参考文献

[1] 廖承恩.微波技术基础[M].西安:西安电子科技大学出版社,2004. [2] 李锦屏,高继森,孙春霞.电磁场与电磁波[M].兰州:兰州大学出版社,2007.

[3] 姚斌,郑勤红,帅春江,张学庆,刘玲丽.用FDTD法计算部分介质填充波导的截止频率[J].云南师范大学学 报.2006,26(4):52-55.

[4] 杨春艳,姚斌,郑勤红.波导截止特性和衰减常数的有限差分分析[J].云南师范大学学报.2005,25(3):35-40. [5] 张璟璟,周立伟,金伟其,张智诠.多重网格法求解三维静电场分布[J].电子学报:2000,28(2):93-96. [6] 盛振华.电磁场微波技术与天线[M].西安:西安电子科技大学出版社,1998. [7] 黄彩华.矩形变形脊波导主模截止波长和特性计算[J].雷达与对抗,1997,18(3):16-22.

[8] 任列辉,邢锋,徐诚等.用电磁场算子理论分析脊波导的传输特性[J].电子与信息学报,2004, 26(2):326-331 鲁加国,樊德森.非对称单脊波导的积分方程法分析[J].微波学报,1998,14(2):108-115. [9] 姚斌.用FDTD法分析介质波导和介质加载腔的电磁波模场分布[D].2006,云南师范大学 [10] 王长清.电磁场计算中的时域有限差分法[M].北京大学出版社,1994. [11] 金建铭.电磁场有限元方法[M].西安:西安电子科技大学出版社,2001.

[12] 邓素芬,杨显清,陈友臸.基于有限元法的脊波导特征值分析[J].电子对抗技术,2005, 20(1): 43-46.

[13] 赵霞.基于有限元法的波导主模特性分析[J].科技信息,2006,14(1):6-9. [14] 谭浩强等.FORTRAN语言[M].北京:清华大学出版社,1990.

[15] 喻志远.任意截面波导的模式截面场的数值分析[J].电波科学学报.2001,16(3):291-294. [16] 金永兴,徐江峰,王剑锋[J].脊形波导有限差分法分析[J].中国计量学院.2003,14(2):114-116.

[17] 胡传顺.用差分法解磁场方程求真空中场空间磁场分布[J].石油化工高等学校学报.1998,11(2):64-68. [18] QIU Zhaojie,HOU Xinyu,WEI Gao,et. A High Accuracy FEM/BEM Hybrid Method for EM Scattering from Two-Dimensional Objects[J] . Journal of Microwares.2006,22(3):1-4.

[19] Ni G Z,Yang S Y,et.Numerical caculation of electromagnetic fields project[M].Beijing:Machine Press,2004. [20] Li Z Y. Electromagnetic field boundary element method [M]. Beijing:Press of Beijing Institute of Technology,1987.

第 20 页 共 35 页

陕西理工学院毕业设计

附录A 外文文献

The FDTD Method in Electromagnetic Field

1 Introduction

FDTD was first introduced as numerical method for solving Maxwell’s equations by Yee in 1966[1]，and pioneering work on the application of this technique in the microwave domain has been done by Taflove[2].Because of its simplicity，robustness and versatility , FDTD has been applied in the optical range of frequencies, as suitable problems to be analyzed with the method have emerged .Planar optical waveguide and graded index optical waveguide are particular-1y good candidates for modeling with FDTD. Their structure is simple enough for the model to be used effectively, and due to the planarity，only a two-dimensional analysis is required for many cases[3-4].With a modem-day super-computer, structures of hundreds of wavelengths long can be analyzed in dimensions.

FDTD is, in essence，an initial-value problem, where an electromagnetic field is allowed to evolve as specified by the sources，in discrete time steps along a lattice including the structure to be analyzed .The evolution of the field is determined by the complex dielectric constants at each cell .At boundaries with differing dielectric constants, reflection，refraction and diffraction can be observed .The time-stepping is usually carried out to several complete cycles of a sinusoidal varying source, and maximum values of the field component magnitudes during a half-cycle after the last complete cycle are stored .In this way，a steady-state solution is achieved ( provided the wave is allowed to propagate through the whole model space).

In the analysis, great care must be taken in choosing the boundary conditions. Since FDTD is an initial-value time-domain method, the propagating wave will reflect at the lattice boundaries unless special conditions are imposed on the fields at the boundary cells to make them non-reflecting ( \" absorbing\").These absorbing boundary conditions absorb fields that are incident on the boundaries such that reflections do not occur. Although at this time there is no perfect absorbing boundary algorithm, many approximate boundary conditions have been developed which minimize any reflections at the lattice boundaries .First-and second-order approximations have been developed for these absorbing boundary conditions, including those by Taflove and Brodwin[5]，Enquist and Majda[6]，and Bayliss and Turkel [7]. This prevents a significant amount of energy from reflecting at the lattice boundaries，Although the absorbing boundaries extend the limits of the lattice and thus reduce space of the structure to be analyzed, FDTD presents savings in memory and execution time. Whereas other methods require storage and

23computation time on the order of 3N and 3N , respectively, where N is the number of cells in the model, FDTD requires only N for both. This is a direct consequence of the time-domain aspect of the method. 2 Formulations

Yee expressed Maxwell' s curl equations in their finite-difference form. The curl equations that are used in the Yee/FDTD algorithm are (We assume the media to be nonmagnetic, i .e.0)

EHE HE tt(1) is electric

Where  is magneto conductivity，is dielecreic constant and

conductivity .Maxwell’s equations in a rectangular coordinate system, which are

第 21 页 共 35 页

陕西理工学院毕业设计

Ex1HzHyEx

tzyEy1HxHzEy tzxEz1HyHz EztxyHx1EyEz

tzy(2a)

(2b)

(2c)

(2d)

Hy1EzEz错误！未找到引用源。 (2e) txzHz1ExEy  tyx(2f)

Fig1.Model of field component in the three

- dimension Yee cell

The components of E and H are positioned about a unit cell of the lattice (see Fig .1) and evaluated alternatively with half-time steps. The goal of the method is to model the propagation of an electromagnetic wave into a volume of space containing dielectric and/or conducting structures by time stepping, i .e.，repeatedly implementing a finite difference analog of the curl equations at each cell of the corresponding space lattice. The incident wave is tracked as it propagates to the dielectric structure and interacts with it via surface current excitation，spreading, penetration, and diffraction .Wave tracking is completed when sinusoidal steady-state behavior is completed for each lattice cell.

By expressing a continuous function of space and rime in its discretized forma function at its(nth)time step can be rewritten as

E(i,j,k)E(ix,jy,kz,Nt)

n(3)

FDTD utilizes the Yee cell for placement of the six electric-and magnetic-field components for computation at nodes of the finite-difference lattice .The dimensions of the cell are x，y and zin the x , y and z directions , respectively .By expanding the three E and H components represented in(1)in their finite-difference approximation, using the central-difference approximation, for both space and time and utilizing x=y=z=s，one can obtain The 1/2’s in(4)and(5)represent the field components' actual placement in the Yee cell, but the calculation is implemented on the computer as node or cell(i, j，k) in the computational lattice .From( 5a)-(5c) it can be noted that each E-field component for the cell ( t , j，k)at the(n+1) st time step depends only on the E-field component’s value at the previous(nth)time step and the surrounding H-field components .A similar procedure is adopted for H-field components except that the magnetic field components are calculated at a 1/2-time-step difference .This provides a routine implementation of the solution procedure through Yee’s cell.

第 22 页 共 35 页

陕西理工学院毕业设计

1111n(1/2))Hx（i,j,k）Hxn-(1/2（i,j,k）2222

11nnE（i,j,k1）-E（i,j,k）yy2211nnE（i,j,k）-E（i,j1,k）zz22t11（i,j,k）s22

(4a)

1111n(1/2)n-(1/2)Hy（i,j,k）Hy（i,j,k）222211nnE（i,j,k）-E（i,j,k1）xx2211nnE（i1,j,k）-E（i,j,k）zz22t（i,j,k）s

1212(4b)

1111))Hzn(1/2（i,j,k）Hzn-(1/2（i,j,k）2222

11nnE（i,j1,k）-E（i,j,k）xx2211nnE（i,j,k）-E（i1,j,k）yy22t（i,j,k）s1212(4c) 错误！

未找到引用源。

11（i,j,k）t11tn1n2E（Ex（i,j,k）i,j,k）x11122（i,j,k）s（i,j,k）s2221111n(1/2)n(1/2)（i,j,k）-H（i,j,k）yy22221111n(1/2)n(1/2)H（i,j,k）-H（i,j,k）zz2222H

(5a)

11（i,j,k）t11tn1n2E（Ey（i,j,k）i,j,k）y1122（i,j,k）s（i,j,k）s22(5b) 错误！

1111n(1/2)n(1/2)H（i,j,k）-H（i,j,k）xx22221111n(1/2)n(1/2)H（i,j,k）-H（i,j,k）zz2222未找到引用源。

第 23 页 共 35 页

陕西理工学院毕业设计

11（i,j,k）t11t1n2E（Ezn（i,j,k）i,j,k）z11(5c) 22（i,j,k）s（i,j,k）s221111n(1/2)n(1/2)H（i,j,k）-H（i,j,k）yy22221111n(1/2)n(1/2)H（i,j,k）-H（i,j,k）xx2222

In order to avoid possible numerical dispersion, numerical instabilities as well as coarse resolution, the

time and space discretizations must be properly chosen .In this paper, we have chosen the space step of 1/20 of the precision dimension of the problem. The time step tis determined by the cell size and must satisfy the stability condition (6)

Vmaxt111 222yzx(6)

Where Vmax is the maximum wave velocity of within the model? For the cubic cell model

x=y=z=s, the following relation is usually taken

tt 2Cmax(7)

Where Cmax is the maximum velocity of electromagnetic waves in the interaction space .We have taken Cmaxcorresponding to the velocity for electromagnetic waves in fat and bone.

In this paper, the absorption derived by Mur (8) are used .the condition for Ezat i=1, for example, can written as

Ezn1(1,j,k1/2)Ezn1(2,j,k1/2)C1tsn1Ez(2,j,k1/2)Ezn1(1,j,k1/2)C1tss(C1t)2n(Ez(1,j,k1/2)Ezn(2,j1,k1/2))2C1ts(C1ts)Ezn(1,j1,k1/2)Ezn(1,j1,k1/2)Ezn(2,j1,k1/2)(C1t)2nnnEz(2,j1,k1/2)Ez(1,j,k3/2)Ez(1,j,k1/2)2s(C1ts)nEz(2,j,k3/2)Ezn(2,j,k1/2)

错误！未找到引用源。Where Ciis the phase velocity observed at the lattice truncation. 3 Conclusion

(8)

This paper introduces the principle of finite-difference time domain method and presents the Maxwell

equations' numerical results in electromagnetic field .The FDTD algorithm can be adopted to do the full-wave analysis on six components of electromagnetic field .It can be able to be applied not only in calculation of weakly guiding optical waveguide ,but also the waveguide structure of bigger difference of refracture index, and because of its small amount of memory capacity , this method is promisingly thought as a valuable numerical method，and will be widely applied the design of optoelectronics component.

第 24 页 共 35 页

陕西理工学院毕业设计

附录B 中文翻译

FDTD法在电磁场中的应用

1 简介

FDTD法是1966由Yee氏首次引入的一种求解麦克斯韦方程组的数值方法，并由Taflove将这一技术在微波领域中做了开创性的应用。由于其简单，健全性和多功能性，FDTD已经在光频率范围内得到广泛的应用，而且适合用这种方法分析的问题也已经出现。对于平面光波导和梯度折射率光波导来说，用FDTD分析是一种尤其好的分析方法。其结构非常简单，是一种可以有效加以利用的模式，而且由于是平面的，在许多情况下[3-4]只需要两维分析即可。用超级计算机，可以分析该结构中的数以百计的波长。

时域有限差分法，实质上是一个初值问题，在电磁场中允许用指明的来源来说明，在离散时间域可以用框架结构来分析。该场的发展是由在每一个小单元中的复杂的介电常数决定的。在不同介电常数的边界上，可以观察到反射，折射和衍射。时间步长通常是几个各不相同的完整的正弦波源周期实现的，场分量的最大值出现在存储的最后一个周期的半个周期之间。通过这种方式，一个稳定的解决办法得以实现（允许波在整个模型空间传播） 。

在分析中，必须十分小心地选择相应的边界条件。由于时域有限差分法是一种初值时域方法，传播波将在格的边界上反射，除非在场的边界格上强加特殊条件，使它们不反射（ “吸收” ） 。这些吸收边界条件的吸收场发生在边界上，就如反射没有发生。虽然在这个时候有没有完善的吸收边界算法，许多近似边界条件算法已经让网格边界上的反射最小。一阶和二阶近似差分形式已经发展为包括Taflove和Brodwin [ 5 ] ，Enquist和Majda [ 6 ] ，和Bayliss和Turkel [ 7 ]所提出的边界条件 。这可以防止网格边界上的大量源反射，虽然吸收边界中的网格可以进行扩展，进行分析的结构可以

2减少空间，提出了FDTD方法节省内存和执行时间。而其他的方法需要存储和计算时间与3N和

33N成正比，其中N是网格总数，而FDTD的计算时间仅需要与N成正比。这是一个更直接有效的时域计算方法。 2公式表达

用Yee氏网格表达的Maxwell旋度方程的有限差分表达形式。用FDTD或Yee氏网格计算该旋度方程的规则如下（我们假定这个媒介是无磁性的i .e.0）

EEHHE （1）

tt这里是磁导率，是介电常数，是电导率，在直角坐标系中，Maxwell旋度方程表达式如

下：

Ex1HzHy Ex tyzEy1HxHz EytzxEz1HyHzEz txyHx1EyEz tzy(2a)

(2b)

(2c)

(2d)

图 1 三维场分量模型

第 25 页 共 35 页

陕西理工学院毕业设计

Hy1EzEz 错误！未找到引用源。 txz(2e)

Hz1ExEy  tyx(2f)

各组成部分E和H的位置结构的单位网格（见图.1 ）这一方法的目标是建立一个电磁波在介质

空间和时间空间进行传播的模型，在相应空间网格上的每个单元中多次用有限差分法模拟了旋度方程。在入射波的传播路径是它传播到电介质结构中与它所经过的介质表面发生作用激发了色散和衍射。当正弦稳态行为在每个网格中完成后使波传播完成。

E或H在直角坐标系中任一分量，在时间和空间中的离散表示为（其中n为时间步长）：

En(i,j,k)E(ix,jy,kz,Nt)

(3)

FDTD方法利用Yee网格将六个E和H的分别放置在有限差分网格的每个节点上。式中x，y和z分别是长方体网格沿x、y、z方向的空间步长。通过展开三个E和H部分来描述如图1所示的有限差分近似的方法，使用中央差分近似法，在空间和时间上，利用x=y=z=s，就可以得到如式（4）和（5）所表示的Yee氏网格中的场分量，但是计算机所执行计算是计算机网格中的节点或者是单元（i，j，k）。从式（5a）—（5c）表明了每个网格上的E分量在n+1上的时间步长，仅仅依赖于他的前一个（第n个）时间步长和他周围的H分量。H场分量采用了一个类似的程序，除了磁场分量仅仅计算在1/2时间步长处。这就提供了一个通过Yee网格解决程序的常规的执行方案。

1111n(1/2))Hx（i,j,k）Hxn-(1/2（i,j,k）222211nnE（i,j,k1）-E（i,j,k）yy2211nnE（i,j,k）-E（i,j1,k）zz22t11（i,j,k）s22

(4a)

1111n(1/2)n-(1/2)Hy（i,j,k）Hy（i,j,k）222211nnE（i,j,k）-E（i,j,k1）xx2211nnE（i1,j,k）-E（i,j,k）zz22t（i,j,k）s

1212(4b)

第 26 页 共 35 页

陕西理工学院毕业设计

1111))Hzn(1/2（i,j,k）Hzn-(1/2（i,j,k）2222

11nnE（i,j1,k）-E（i,j,k）xx2211nnE（i,j,k）-E（i1,j,k）yy22t（i,j,k）s1212(4c) 错误！

未找到引用源。

11（i,j,k）t11tn1n2E（Ex（i,j,k）i,j,k）x11122（i,j,k）s（i,j,k）s2221111n(1/2)n(1/2)H（i,j,k）-H（i,j,k）yy22221111n(1/2)n(1/2)H（i,j,k）-H（i,j,k）zz2222

(5a)

11（i,j,k）t11tn1n2E（Ey（i,j,k）i,j,k）y1122（i,j,k）s（i,j,k）s221111n(1/2)n(1/2)（i,j,k）-H（i,j,k）xx22221111n(1/2)n(1/2)H（i,j,k）-H（i,j,k）zz2222H错误！(5b)

未找到引用源。

11（i,j,k）t11t1n2E（Ezn（i,j,k）i,j,k）z1122（i,j,k）s（i,j,k）s221111n(1/2)n(1/2)H（i,j,k）-H（i,j,k）yy22221111n(1/2)n(1/2)H（i,j,k）-H（i,j,k）xx2222

(5c)

为了避免可能造成数值色散和数值不稳定性，以及误差大的计算，对于时间和空间的离散化必

须要有适当的选择。在本文中，我们选择1/20的空间步长来解决精确尺寸问题。时间不常是有网格大小来决定，而且必须满足下列稳定条件

Vmaxt111 222yzx(6)

在这个空间模型中，哪里的波速度Vmax最大？这个立体的空间模型中x=y=z=s，通常

会用到下面的关系

第 27 页 共 35 页

陕西理工学院毕业设计

tt 2Cmax(7)

在相互作用空间中，电磁波Cmax最大。我们可以得到电磁波在边界或其他空间位置中的相应的Cmax。

在本文中，我们以Ez在i=1处的条件为例，能够得到

Ezn1(1,j,k1/2)Ezn1(2,j,k1/2)C1tsn1Ez(2,j,k1/2)Ezn1(1,j,k1/2)C1tss(C1t)2nn2C1ts(C1ts)(Ez(1,j,k1/2)Ez(2,j1,k1/2))Ezn(1,j1,k1/2)Ezn(1,j1,k1/2)Ezn(2,j1,k1/2)(C1t)2nnnEz(2,j1,k1/2)Ez(1,j,k3/2)Ez(1,j,k1/2)2s(C1ts)nEz(2,j,k3/2)Ezn(2,j,k1/2)这里所说的在网格截断处所得到的是相位速度（Ci）。

(8)

3 结论

本文主要介绍了时域有限差分法的计算规则，并介绍了麦克斯韦方程组在电磁场领域中的数值计算结果。时域有限差分算法可以用来进行电磁场中的六个场分量的全波分析。它不仅能够适用于弱指导光波导的计算，而且更大的不同的波导结构指数，并且由于它所占的存储容量小，这种方法是一种非常有前景的宝贵的数值方法，它将广泛应用于光电元件的设计当中。

第 28 页 共 35 页

陕西理工学院毕业设计

附录C程序设计

C1单脊波导截止频率计算： dimension ez(18,25) common ez

data maxt,alfa,sqkhold,hh/400.0,1.3,0.25,20.0e-3/ do 1 i=1,17 do 1 j=1,21 1 ez(i,j)=0.0 do 14 j=2,10 do 14 i=2,16 14 ez(i,j)=4.0 do 12 j=11,20 do 12 i=10,16 12 ez(i,j)=4.0

alf=alfa/(4.0-sqkhold) konv=0 iter=0

2 do 4 k=1,5 do 3 j=2,10

resdl=2.0*ez(2,j)+ez(1,j-1)+ez(1,j+1) +(sqkhold-4.0)*ez(1,j) ez(1,j)=ez(1,j)+alf*resdl do 3 i=2,16

resdl=eh(i,j)+(sqkhold-4.0)*ez(i,j) 3 ez(i,j)=ez(i,j)+alf*resdl do 4 j=11,20 do 4 i=10,16

resdl=eh(i,j)+(sqkhold-4.0)*ez(i,j) 4 ez(i,j)=ez(i,j)+alf*resdl iter=iter+1

if(konv.eq.1) goto 7 rcn=0 rcd=0

do 9 j=2,10

rcn=rcn+(2.0*ez(2,j)+ez(1,j-1)+ez(1,j+1) -4.0*ez(1,j))*ez(1,j) rcd=rcd+ez(1,j)**2 do 9 i=2,16 rcn=rcn+ez(i,j)*(eh(i,j)-4.0*ez(i,j)) 9 rcd=rcd+ez(i,j)**2 do 6 j=11,20 do 6 i=10,16 rcn=rcn+ez(i,j)*(eh(i,j)-4.0*ez(i,j)) 6 rcd=rcd+ez(i,j)**2 sqkhnew=-rcn/rcd if(abs((sqkhnew-sqkhold)/sqkhnew).lt.0.5e-4) konv=1

第 29 页 共 35 页

陕西理工学院毕业设计

if(iter.ge.maxt) goto 7 sqkhold=sqkhnew alf=alfa/(4.0-sqkhold) goto 2

7 freq=sqrt(sqkhold)*20.0*300.0/6.2832/hh open(6,file='res.txt') write(6,170) freq.iter close(6)

170 format(1x,2hf=,e12.4,1x,5hiter=,i4) stop end function eh(i,j) common ez dimension ez(18,25) eh=ez(i,j-1)+ez(i,j+1)+ez(i-1,j)+ez(i+1,j) return end

C2梯形脊高变化截止频率计算： dimension ez(100,100) common ez

data maxt,alfa,sqkhold,hh/400.0,1.3,0.25,20.0e-3/ pi=3.1415926 do 1 i=1,27 do 1 j=1,21 1 ez(i,j)=0.0 do 14 j=2,20 do 14 i=2,9 14 ez(i,j)=4.0 do 13 j=10,20

do 13 i=10,(20-j+10) 13 ez(i,j)=4.0

do 10 j=2,9 do 10 i=10,21 10 ez(i,j)=4.0

do 11 j=2,9 do 11 i=22,26 11 alf=alfa/(4.0-sqkhold) konv=0 iter=0 2 do 5 k=1,5 do 3 j=2,20

resdl=2.0*ez(2,j)+ez(1,j-1)+ez(1,j+1) +(sqkhold-4.0)*ez(1,j)

第 30 页 共 35 页

陕西理工学院毕业设计

ez(1,j)=ez(1,j)+alf*resdl do 3 i=2,9

resdl=eh(i,j)+(sqkhold-4.0)*ez(i,j) 3 ez(i,j)=ez(i,j)+alf*resdl do 4 j=10,20

do 4 i=10,(20-j+10)

resdl=eh(i,j)+(sqkhold-4.0)*ez(i,j)

4 ez(i,j)=ez(i,j)+alf*resdl do 5 j=2,9

do 5 i=10,21

resdl=eh(i,j)+(sqkhold-4.0)*ez(i,j)

5 ez(i,j)=ez(i,j)+alf*resdl do 5 j=2,9

do 5 i=22,26

resdl=eh(i,j)+(sqkhold-4.0)*ez(i,j)

6 ez(i,j)=ez(i,j)+alf*resdl iter=iter+1

if(konv.eq.1) goto 7 rcn=0 rcd=0

do 12 j=2,20

rcn=rcn+(2.0*ez(2,j)+ez(1,j-1)+ez(1,j+1) -4.0*ez(1,j))*ez(1,j) rcd=rcd+ez(1,j)**2 do 12 i=2,9 rcn=rcn+ez(i,j)*(eh(i,j)-4.0*ez(i,j)) 12 rcd=rcd+ez(i,j)**2 do 8 j=10,20 do 8 i=10,(20-j+10) rcn=rcn+ez(i,j)*(eh(i,j)-4.0*ez(i,j)) 8 rcd=rcd+ez(i,j)**2 do 6 j=2,9 do 6 i=10,26

6 rcd=rcd+ez(i,j)**2 do 7 j=2,9 do 7 i=10,26

rcn=rcn+ez(i,j)*(eh(i,j)-4.0*ez(i,j)) 7 rcd=rcd+ez(i,j)**2 sqkhnew=-rcn/rcd if(abs((sqkhnew-sqkhold)/sqkhnew).lt.0.5e-4) konv=1 if(iter.ge.maxt) goto 7 sqkhold=sqkhnew alf=alfa/(4.0-sqkhold) goto 2

9 freq=sqrt(sqkhold)*20.0*300.0/6.2832/hh

第 31 页 共 35 页

陕西理工学院毕业设计

open(6,file='res.txt') write(6,170) freq,iter close(6)

170 format(1x,2hf=,e12.4,1x,5hiter=,i4) stop end function eh(i,j) common ez dimension ez(100,100) eh=ez(i,j-1)+ez(i,j+1)+ez(i-1,j)+ez(i+1,j) return end

C3梯形脊宽度变化截止频率计算： dimension ez(100,100) common ez

data maxt,alfa,sqkhold,hh/400.0,1.3,0.25,20.0e-3/ pi=3.1415926 do 1 i=1,28 do 1 j=1,21 1 ez(i,j)=0.0 do 14 j=2,20 do 14 i=2,12 14 ez(i,j)=4.0 do 13 j=10,20

do 13 i=13,(22-j+10) 13 ez(i,j)=4.0

do 10 j=2,9

do 10 i=13,(22-j+10) 10 ez(i,j)=4.0

do 15j=2,9 do 15 i=24,27 15 ez(i,j)=4.0

alf=alfa/(4.0-sqkhold) konv=0 iter=0 2 do 5 k=1,5 do 3 j=2,20

resdl=2.0*ez(2,j)+ez(1,j-1)+ez(1,j+1) +(sqkhold-4.0)*ez(1,j) ez(1,j)=ez(1,j)+alf*resdl do 3 i=2,13

resdl=eh(i,j)+(sqkhold-4.0)*ez(i,j) 3 ez(i,j)=ez(i,j)+alf*resdl do 4 j=10,20

第 32 页 共 35 页

陕西理工学院毕业设计

do 4 i=13,(22-j+10)

resdl=eh(i,j)+(sqkhold-4.0)*ez(i,j)

4 ez(i,j)=ez(i,j)+alf*resdl do 16 j=2,9

do 16 i=13,(22-j+10)

resdl=eh(i,j)+(sqkhold-4.0)*ez(i,j)

16 ez(i,j)=ez(i,j)+alf*resdl do 5 j=2,9

do 5 i=24,27

resdl=eh(i,j)+(sqkhold-4.0)*ez(i,j)

5 ez(i,j)=ez(i,j)+alf*resdl iter=iter+1

if(konv.eq.1) goto 7 rcn=0 rcd=0

do 12 j=2,20

rcn=rcn+(2.0*ez(2,j)+ez(1,j-1)+ez(1,j+1) -4.0*ez(1,j))*ez(1,j) rcd=rcd+ez(1,j)**2 do 12 i=2,13 rcn=rcn+ez(i,j)*(eh(i,j)-4.0*ez(i,j)) 12 rcd=rcd+ez(i,j)**2 do 8 j=10,20 do 8 i=13,(23-j+10) rcn=rcn+ez(i,j)*(eh(i,j)-4.0*ez(i,j)) 8 rcd=rcd+ez(i,j)**2 do 11 j=2,9 do 11 i=13,(23-j+10) 11 rcd=rcd+ez(i,j)**2 do 6 j=2,9 do 6 i=23,27

rcn=rcn+ez(i,j)*(eh(i,j)-4.0*ez(i,j)) 6 rcd=rcd+ez(i,j)**2 sqkhnew=-rcn/rcd if(abs((sqkhnew-sqkhold)/sqkhnew).lt.0.5e-4) konv=1 if(iter.ge.maxt) goto 7 sqkhold=sqkhnew alf=alfa/(4.0-sqkhold) goto 2

7 freq=sqrt(sqkhold)*20.0*300.0/6.2832/hh open(6,file='res.txt') write(6,170) freq,iter close(6)

170 format(1x,2hf=,e12.4,1x,5hiter=,i4)

第 33 页 共 35 页

陕西理工学院毕业设计

stop end function eh(i,j) common ez dimension ez(100,100) eh=ez(i,j-1)+ez(i,j+1)+ez(i-1,j)+ez(i+1,j) return end

第 34 页 共 35 页