基于Matlab的光子晶体波导仿真研究

仿真技术

文章编号：１００８－０５７０（２００７）０１－１－０２６４－０２

中文核心期刊《微计算机信息》（测控自动化）２００７年第２３卷第１－１期

基于Ｍａｔｌａｂ的光子晶体波导仿真研究

ＦＤＴＤＡｒｉｔｈｍｅｔｉｃＰｒｏｇｒａｍｍｉｎｇｉｎｐｈｏｔｏｎｉｃｃｒｙｓｔａｌｓｗａｖｅｇｕｉｄｅｏｎｔｈｅＢａｓｉｓｏｆＭａｔｌａｂＬａｎｇｕａｇｅ

（江苏大学）吴炳坚

沈廷根

ＷＵＢＩＮＧＪＩＡＮＳＨＥＮＴＩＮＧＧＥＮ

摘要：本文首先介绍了光子晶体波导的原理，然后分析时域有限差分法微分方程及边界条件，最后运用ｍａｔｌａｂ语言实现二维

光子晶体波导的仿真。

关键词：光子晶体波导；时域有限差分法；ｍａｔｌａｂ语言中图分类号：ＴＰ３９１文献标识码：Ａ

Ａｂｓｔｒａｃｔ：ｗｅｉｎｔｒｏｄｕｃｅｔｈｅｔｈｅｏｒｙｏｆｐｈｏｔｏｎｉｃｃｒｙｓｔａｌｓｗａｖｅｇｕｉｄｅ，ｔｈｅｎａｎａｌｙｚｅＦＤＴＤａｒｉｔｈｍｅｔｉｃａｎｄｂｏｕｎｄａｒｙｃｏｎｄｉｔ．Ｆｉｎａｌｌｙ，ｗｅｕｓｅＭａｔｌａｂｌａｎｇｕａｇｅｔｏｐｒｏｇｒａｍｍｅＦＤＴＤａｒｉｔｈｍｅｔｉｃｉｓｃａｒｒｉｅｄｏｕｔ，ｔｈｅｒｅｓｕｌｔｓｐｒｅｓｅｎｔａｔｈｅｏｒｅｔｉｃａｌｂａｓｅｆｏｒｐｒｏｃｅｓｓｉｎｇｐｈｏｔｏｎｉｃｄｅｖｉｃｅ．Ｋｅｙｗｏｒｄｓ：ｐｈｏｔｏｎｉｃｃｒｙｓｔａｌｓｗａｖｅｇｕｉｄｅ，ＦＤＴＤ，Ｍａｔｌａｂｌａｎｇｕａｇｅ

技术创新

１引言

光子晶体最早是在１９８７年由Ｅ．Ｙａｂｌｏｎｏｖｉｔｃｈ与Ｓ．Ｊｏｈｎ提出，是由不同介电系数的物质周期性排列所组成。和在半导体晶体中电子会形成能带结构、带与带间有能隙的情形类似，在光子晶体中传播的光波色散曲线，因介质介电系数周期性排列的缘故，也形成带状结构即“光子禁带”。特定频率的光因为落在此禁带内而被禁止传播。

我们通过利用这一性质做成光子晶体波导，为整个回路信号的传输提供一个快速通道，能为将来的超大规模全光或光电子集成回路提供一个优良的基础物质平台。

时域有限差分（ＦＤＴＤ）方法自Ｙｅｅ（１９６６年）提出以来就得到迅速发展，其主要思想是以Ｙｅｅ元胞为空间电磁场离散单元，将麦克斯韦旋度方程转化为差分方程，在时间轴和空间轴上逐步推进地求解，最终求出空间场的分布。本文主要结合

分别是媒质的σ是电导率（ｓ／ｍ），ｓ是磁电阻率（Ω／ｍ）μ和ε磁导率和电容率。

图１Ｙｅｅ交错晶格与垂直的Ｘ－Ｙ平面

ＦＤＴＤ算法边界条件特点，用Ｍａｔｌａｂ语言进行编程模拟实现

光子晶体波导传播。

２ＦＤＴＤ算法和边界处理

２．１时域有限差分法

时域有限差分法是将麦克斯韦方程式之中对空间与时间的微分，改以差分来代替。在此结构下，每一个电场分量由磁场分量围绕。由安培定律可根据周围磁场，决定中间电场的变化量。每一个磁场分量由电场分量围绕。由法拉第定律，可根据周围电场决定中间磁场的变化量。如图（１，２）所示对于二维横磁模（ＴＭ）

图２ｘ＝ｉ、ｙ＝ｊ的平面晶格

２．２边界条件

在图１中，最边缘的Ｈｘ、Ｈｙ与Ｅｚ由于没有更边缘的电磁

场值可供计算，故无法确定其值。所以，如以零为默认值，则波传到边缘会好像碰到金属，波会反弹回解析空间。我们希望波在未到边缘前，其电磁场值，便可衰减为零，如此便可与边缘的预设零值相接。此处我们以Ｂｅｒｅｎｇｅｒ的ＰＭＬ边界理论（如图３所示），来讨论程序的写法。

，Ｈｚ＝０。我们得到ＦＤＴＤ的基本方程：

吴炳坚：硕士研究生

基金项目：江苏省自然科学基金（ＢＫ２００４０５９）

－２６４－

＊图３：ＰＭＬ内，各区域的导电率σ与导磁率σ分布示意图，＊＊

在σ由内层至外层逐渐增加。、＝０及σ≠的区域中，σσ

３６０元／年邮局订阅号：８２－９４６现场总线技术应用２００例》

您的论文得到两院院士关注

＊

理论上σσＸ愈大，则波衰减愈快，但为解析空间的精准Ｘ、

仿真技术

Ｓｉｇｍａｘ＿ｙ１＝ｒｅｐｍａｔ（Ｓｉｇｍａｘ（１：２：ＮＰＭＬ＊２－１）＇，１，ＮＰＭＬ）；

Ｓｉｇｍａｙ＿ｚ１＝ｆｌｉｐｌｒ（ｒｅｐｍａｔ（Ｓｉｇｍａｘ（２：２：ＮＰＭＬ＊２），ＮＰＭＬ，１））；Ｓｉｇｍａｙ＿ｘ１＝ｆｌｉｐｌｒ（ｒｅｐｍａｔ（Ｓｉｇｍａｘ（１：２：ＮＰＭＬ＊２－１），ＮＰＭＬ，１））；Ｓｉｇｍａｙ＿ｙ１＝ｆｌｉｐｌｒ（ｒｅｐｍａｔ（Ｓｉｇｍａｘ（２：２：ＮＰＭＬ＊２），ＮＰＭＬ，１））；％这里只设置一个区域的ＰＭＬ层，其余类同。

２．３．４程序模拟效果图

度，则交界处波的连续性要良好，因此σ通常σ也不能太大。Ｘ、

＊

Ｘ

＊

σσσＸ在ＰＭＬ内的安排是采渐进式，最外层，σＸ、Ｘ（Ｌ）、Ｘ（ｍａｘ）

且：

其中Ｒ（０）是入射角零度的反射率，Ｍ为一个介于２～４的系数，Ｃ为真空中的光速，Ｌ为ＰＭＬ的总层数，Δｘ为每一层的宽度。

ＦＤＴＤ网格形式的划分，必须在截断处设置适当的吸收边界条

件，以便用有限网格空间模拟开放的无限空间或无限长的传输结构。在网格的划分过程中，时间步长取为

其中ｃ为自由空间光速。最后进行波源的设置。

，

２．３Ｍａｔｌａｂ语言部分源代码

２．３．１初始化设置

Ｗ＝０．３６；％设定光子晶体波源频率，需要注意的是只有在

光子晶体禁带范围的频率才能产生波导

图４二维光子晶体波导结构图

ＭＬａｔｘ＝１１；％ｘ方向的晶格数ＭＬａｔｙ＝１１；％ｙ方向的晶格数ＮＭｌａｔ＝２１；％每个晶格的网格数

ＮＴｘ＝ＭＬａｔｘ＊ＮＭｌａｔ＋１；％ｘ方向的总网格数ＮＴｙ＝ＭＬａｔｙ＊ＮＭｌａｔ＋１；％ｙ方向的总网格数Ｒ＝０．２；％晶格半径ｅａ＝１１．４；％介电常数

ｍｕ０＝４＊ｐｉ＊１．０ｅ－７；％真空中的电容率ｅ０＝８．８５＊１ｅ－１２；％真空中的介电常数ｃ＝１／ｓｑｒｔ（ｍｕ０＊ｅ０）；％光速

Ｅｚ（１，Ｎｐｙ－（ＮＭｌａｔ－１）／２：Ｎｐｙ＋（ＮＭｌａｔ－１）／２）＝Ｅｚ（１，Ｎｐｙ－（ＮＭｌａｔ－１）／２：Ｎｐｙ＋（ＮＭｌａｔ－１）／２）＋ｓｉｎ（Ｗ＊ｍ＊Ｄｔ）；

％模拟连续平面波源ｅｘｐ（－（ｍ＊Ｗ＊Ｄｔ－３）＾２）２．３．２数学模拟ＦＤＴＤ方程

Ｈｘ（：，２：ＮＴｙ－１）＝Ｈｘ（：，２：ＮＴｙ－１）－Ｄｔ＊（Ｅｚ（：，２：ＮＴｙ－１）－Ｅｚ（：，１：ＮＴｙ－２））／（Ｄｙ＊ｍｕ０）；％我们考虑的是无电或磁损耗

Ｈｙ（２：ＮＴｘ－１，：）＝Ｈｙ（２：ＮＴｘ－１，：）＋Ｄｔ＊（Ｅｚ（２：ＮＴｘ－１，：）－Ｅｚ（１：ＮＴｘ－

的媒质ｓ＝０．２，：））／（Ｄｘ＊ｍｕ０）；

Ｅｚ＝Ｅｚ＋Ｄｔ＊（（Ｈｙ（２：ＮＴｘ，１：ＮＴｙ－１）－Ｈｙ（１：ＮＴｘ－１，１：ＮＴｙ－１））／Ｄｘ－（Ｈｘ（１：ＮＴｘ－１，２：ＮＴｙ）－Ｈｘ（１：ＮＴｘ－１，１：ＮＴｙ－１））／Ｄｙ）．／Ｅｐ

２．３．３边界设置ｎ＝４；％ＰＭＬ层系数Ｒ＝１ｅ－１０；

Ｄｅｌｔａ＝ＮＰＭＬ＊Ｄｘ；

ＳｉｇｍａＭａｘ＝－（ｎ＋１）＊ｅ０＊ｃ＊ｌｏｇ（Ｒ）／（Ｄｅｌｔａ＊２）；ＮＵＭ＝ＮＰＭＬ＊２：－１：１；

Ｓｉｇｍａｘ＝ＳｉｇｍａＭａｘ＊（（ＮＵＭ＊Ｄｘ／２＋Ｄｘ／２）．＾（ｎ＋１）－（ＮＵＭ＊Ｄｘ／２－Ｄｘ／２）．＾（ｎ＋１））／（Ｄｅｌｔａ＾ｎ＊Ｄｘ＊（ｎ＋１））；

Ｓｉｇｍａｙ＝Ｓｉｇｍａｘ；

ＳｉｇｍａＢｏｕｎｄ＝ＳｉｇｍａＭａｘ＊（Ｄｘ／２）．＾（ｎ＋１）／（Ｄｅｌｔａ＾ｎ＊Ｄｘ＊（ｎ＋１））；ＥｚｘＰＭＬ１＝ｚｅｒｏｓ（ＮＰＭＬ，ＮＰＭＬ）；ＥｚｙＰＭＬ１＝ｚｅｒｏｓ（ＮＰＭＬ，ＮＰＭＬ）；ＨｘＰＭＬ１＝ｚｅｒｏｓ（ＮＰＭＬ，ＮＰＭＬ）；ＨｙＰＭＬ１＝ｚｅｒｏｓ（ＮＰＭＬ，ＮＰＭＬ）；

Ｓｉｇｍａｘ＿ｚ１＝ｒｅｐｍａｔ（Ｓｉｇｍａｘ（２：２：ＮＰＭＬ＊２）＇，１，ＮＰＭＬ）；Ｓｉｇｍａｘ＿ｘ１＝ｒｅｐｍａｔ（Ｓｉｇｍａｘ（２：２：ＮＰＭＬ＊２）＇，１，ＮＰＭＬ）；

ＰＬＣ技术应用２００例》

图５当ＦＤＴＤ时间模拟３８０步时图形

我们发现处在光子晶体禁带的光频率只在我们预先设计好的通道中通过，而不处于禁带的光频率通过波导时会产生散射现象，程序结果较好的验证了波导理论，为光子晶体集成化，提供实验依据。

技术创新

图６当ＦＤＴＤ模拟２５００步时图形

图７取频率ｗ＝０．５时出现的光波图

３结论

本文通过运用时域有限差分法及ＰＭＬ边界条件来仿真光子晶体波导克服了以往用平面波展开法带来的计算复杂、精确度不高、占用计算机内存空间大、计算时间周期长等缺点，同时借助Ｍａｔｌａｂ工具在数值计算和科学绘图等方面显示出比其他编程工具所难以比拟的优越性，只要经过很短的时间就能得出精确的结果。（下转第３００页）

邮局订阅号：８２－９４６３６０元／年－

２６５－

软件天地

中文核心期刊《微计算机信息》（测控自动化）２００７年第２３卷第１－１期

产品的智能设计系统中。文中的算法对大多数机电产品的设计都有较好的通用性，目前已应用在电器设计的专家系统中。

补充说明：上文中所述情况，说明用户提供的设计信息“约定④”不完整。为此，定义ＣＦ＝

为设计方案的可

前提条件的属性个数；｜Ｚ｜＝ｎ，为信度。这里，｜Ｚ＇｜为满足“约定④”

属性集的元素个数。｜Ｚ＇｜／｜Ｚ｜偏小说明设计信息不完整，ＣＦ偏小则综合体现了信息的不确定。此时，应终止过程，启“模糊识别”动ＣＢＲ。

３．３．４采用模糊识别方法求解结构设计问题的应用实例

永磁式或者拍设：欲设计一种电磁继电器，要求是旋转式、

合式，动作电压１８Ｖ，额定电压２７Ｖ，输入回路电流大于４Ａ时动作，工作环境：冲击加速度为１０ｇ、２０ｇ、３０ｇ、４０ｇ、５０ｇ的可能性分别为０．２、０．７、１．０、０．６、０．１；振动频率１０００Ｈｚ；振动加速度

９ｇ￣１１ｇ。

有关设计信息如下：①按前文方法确定属性集

＝｛冲击加速度，振动频率，振动加速度｝；②根据式（２），

把ｚ３的属性值处理为，式中常数项θ取

为２．０；ｚ１的属性值表达为模糊集

；

技术创新

③根据专家经验，确定权重集Ｗ＊＝｛０．４０，０．３０，０．３０｝；④从知识库中调出有关特征模型的数据：Ｍ＝｛旋转式，永磁式，拍合式｝＝

参考文献：

［１］ＴａｙＦｒａｎｃｉｓＥ．Ｈ．ａｎｄＧｕＪｉｎｘｉａｎｇ，＂Ｐｒｏｄｕｃｔｍｏｄｅｌｉｎｇｆｏｒｃｏｎ－ｃｅｐｔｕａｌｄｅｓｉｇｎｓｕｐｐｏｒｔ，＂ＣｏｍｐｕｔｅｒｓｉｎＩｎｄｕｓｔｒｙ，ｖｏｌ．４８，ｎｏ．２，ｐｐ．１４３－１５５，２００２．

［２］ＱｉｎＳ．Ｆ，ＨａｒｒｉｓｏｎＲ，ＷｅｓｔＡ．Ａ，ｅｔａｌ，＂Ａｆｒａｍｅｗｏｒｋｏｆｗｅｂ－ｂａｓｅｄｃｏｎｃｅｐｔｕａｌｄｅｓｉｇｎ，＂ＣｏｍｐｕｔｅｒｓｉｎＩｎｄｕｓｔｒｙ，ｖｏｌ．５０，ｎｏ．２，ｐｐ．１５３－１６４，２００３．

［３］崔芮华，柳杰，贾峰，等．数据库技术在ＣＡＤ软件开发中的应用．微计算机信息，２００６，２２（２－３）：８３￣８５．

作者简介：李玲玲（１９６８．２－），女，博士，副教授，硕士生导师，专业为电机与电器；执教于河北工业大学电气与自动化学院．Ｅ－ｍａｉｌ：ｈａｏｈａｏｈａｏ＠ｅｙｏｕ．ｃｏｍ。

Ｂｉｏｇｒａｐｈｙ：ＬｉＬｉｎｇｌｉｎｇ，ｂｏｒｎｉｎ１９６８．２，ｆｅｍａｌｅ，ｄｏｃｔｏｒ，ａｓｓｏｃｉａｔｅｐｒｏｆｅｓｓｏｒ，ｓｕｐｅｒｖｉｓｏｒｏｆｍａｓｔｅｒ，ｔｈｅｍａｊｏｒｉｓＥｌｅｃｔｒｉｃＭａｃｈｉｎｅｓａｎｄＥｌｅｃｔｒｉｃＡｐｐａｒａｔｕｓ．ＳｈｅｉｓｗｏｒｋｉｎｇｉｎＴｈｅＤｅｐａｒｔｍｅｎｔｏｆＥｌｅｃｔｒｉｃａｌＥｎｇｉｎｅｅｒｉｎｇａｎｄＡｕｔｏｍａｔｉｏｎ，ＨｅｂｅｉＵｎｉｖｅｒｓｉｔｙｏｆＴｅｃｈｎｏｌｏｇｙ．Ｅ－ｍａｉｌ：ｈａｏｈａｏｈａｏ＠ｅｙｏｕ．ｃｏｍ．

（３００１３０天津河北工业大学电气与自动化学院）李玲玲张惠娟

王建楠高朝晖

（ＴｈｅＤｅｐａｒｔｍｅｎｔｏｆＥｌｅｃｔｒｉｃａｌＥｎｇｉｎｅｅｒｉｎｇａｎｄＡｕｔｏｍａｔｉｏｎ，

按照前文所述的方法，算得

０．８８１＋０．３０×０．９４３＝０．９１；同理，算得

０．５８。因此，本产品应选择旋转式结构。

３．４ＣＢＲ在产品结构设计中的使用方法

＝０．４０×０．９０３＋０．３０×＝０．６２；＝

ＨｅｂｅｉＵｎｉｖｅｒｓｉｔｙｏｆＴｅｃｈｎｏｌｏｇｙ，Ｔｉａｎｊｉｎ３００１３０）ＬｉＬｉｎｇｌｉｎｇ

ＺｈａｎｇＨｕｉｊｕａｎＷａｎｇＪｉａｎｎａｎＧａｏＺｈａｏｈｕｉ

通讯地址：（３００１３０天津天津市红桥区丁字沽光荣道８号２５１信箱）李玲玲

（收稿日期：２００６．５．２５）（修稿日期：２００６．６．２８）

（上接第２６５页）

本文作者的创新点：将传统计算电磁场与电磁波传播的时域有限差分法（ＦＤＴＤ）运用于光子晶体波导提高计算光传播的精确度；ＰＭＬ边界条件的设置解决了光子晶体波导边缘的反射问题；借助Ｍａｔｌａｂ工具仿真实现光在光子晶体中传播的可视化。

无论是在概念设计阶段还是在详细设计阶段应用ＣＢＲ，都需检索、匹配实例。这两种情况，不同之处在于ＣＢＲ的推理起点不同。前者的推理从的终止处开始，而后者的推“模糊识别”理从对产品部件的检索开始。若采用面向对象的语言开发系统，可以为每个特征模型定义一个实例表，表中记录与该模型同属一个产品子类的实例总数和实例标识符，这些标识符都是指针类型的数据，指向各有关实例在实例库中的检索入口，系统可据此直接去实例库中检索实例。

产品设计中，由模糊性引起的不确定性信息大量存在，所以ＣＢＲ也应在模糊机制下展开。步骤为：（１）对有关实例进行模糊化处理，以正态模糊数表达其属性值；（２）以有关实例为一组典型模式，以设计问题为待识别模式，采用模糊识别的方法进行实例匹配；（３）评价相似实例是否满足设计要求；若否，则修改实例，逐步确定产品部件、零件的几何尺寸，完成详细设计。

在结构的详细设计中，文献介绍的一些关键技术可用来开发系统。

参考文献：

［１］Ｅ．Ｙａｂｌｏｎｏｖｉｔｃｈ，Ｐｈｙｓ．Ｒｅｖ．Ｌｅｔｔ．１９８７，５８：２０５９～２０６２．

２４８９．［２］Ｓ．Ｊｏｈｎ，Ｐｈｙｓ．Ｒｅｖ．Ｌｅｔｔ．１９８７，５８：２４８６～

［３］刘青，李承芳，张翔．二维光子晶体带隙结构的研究，激光杂志，２００２．

［４］高本庆．时域有限差分法．国防工业出班社．１９９５．

［５］李曦，唐琳．在ＭＡＴＬＡＢ中运行其他可执行软件的方法［Ｊ］．微计算机信息，２００６，３：２１２－２１４．

［６］葛德彪，闫玉波．电磁波时域有限差分方法．西安电子科技大学出版社，２００１．

作者简介：吴炳坚（１９８１－），男，汉，江苏启东人，计算机与通信学院，硕士研究生．Ｅ－ｍａｉｌ：ｗｂｊｆｌｙ＠ｇｍａｉｌ．ｃｏｍ；沈廷根，男，汉，教授，计

算机与通信学院。

４结论

本文创新点：（１）在建立递阶划分、“产品－部件－零件”分解的型结构基础上，提出了联合采用聚类分析、模“树”糊识别、ＣＢＲ进行产品结构设计的方法。通过把产品各技术属性按照的拓扑结构进行层次划分，采用自顶向下“树”的思想求解产品的结构设计问题；（２）针对设计信息的不完整、不确定性，提出了相应的处理方法。

这些方法适宜作为一种推理机制和问题求解策略，应用在

－３００－

Ｂｉｏｇｒａｐｈｙ：ＷｕＢｉｎｇｊｉａｎ（１９８１－），ｍａｌｅ，ｈａｎ，ｑｉｄｏｎｇ，ｊｉａｎｇｓｕ，ｓｃｈｏｏｌ

ｏｆｃｏｍｐｕｔｅｒａｎｄｃｏｍｍｕｎｉｃａｔｉｏｎ，Ｂ．Ｅ

（２１２０１３江苏镇江市江苏大学计算机科学与通信工程学院）

吴炳坚沈廷根

（ｓｃｈｏｏｌｏｆｃｏｍｐｕｔｅｒｓｃｉｅｎｃｅａｎｄｔｅｌｅｃｏｍｍｕｎｉｃａｔｉｏｎ，ｊｉａｎｇｓｕｕ－ｎｉｖｅｒｓｉｔｙ，ｚｈｅｎｊｉａｎｇ，２１２０１３）ＷｕＢｉｎｇｊｉａｎＳｈｅｎＴｉｎｇｇｅｎ

通迅地址：（２１２０１３江苏镇江市江苏省镇江市江苏大学校本部８８０信箱）吴炳坚

（收稿日期：２００６．５．２７）（修稿日期：２００６．６．２５）

３６０元／年邮局订阅号：８２－９４６现场总线技术应用２００例》