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

（Ｈｘ（１：ＮＴｘ－ １，２：ＮＴｙ）－ Ｈｘ（１：ＮＴｘ－ １，１：ＮＴｙ－ １））／Ｄｙ）．／Ｅｐ
２．３．３ 边界设置 ｎ＝４； ％ＰＭＬ 层系数 Ｒ＝１ｅ－ １０； Ｄｅｌｔａ＝ＮＰＭＬ＊Ｄｘ； ＳｉｇｍａＭａｘ＝－ （ｎ＋１）＊ｅ０＊ｃ＊ｌｏｇ（Ｒ）／（Ｄｅｌｔａ＊２）； ＮＵＭ＝ＮＰＭＬ＊２：－ １：１； Ｓｉｇｍａｘ＝ＳｉｇｍａＭａｘ＊ （（ＮＵＭ＊Ｄｘ／２ ＋Ｄｘ／２）．＾（ｎ ＋１）－ （ＮＵＭ＊Ｄｘ／２ － Ｄｘ／２）．＾（ｎ＋１））／（Ｄｅｌｔａ＾ｎ＊Ｄｘ＊（ｎ＋１））； Ｓｉｇｍａｙ＝Ｓｉｇｍａｘ； ＳｉｇｍａＢｏｕｎｄ＝ＳｉｇｍａＭａｘ＊（Ｄｘ／２）．＾（ｎ＋１）／（Ｄｅｌｔａ＾ｎ＊Ｄｘ＊（ｎ＋１））； ＥｚｘＰＭＬ１＝ｚｅｒｏｓ（ＮＰＭＬ，ＮＰＭＬ）； ＥｚｙＰＭＬ１＝ｚｅｒｏｓ（ＮＰＭＬ，ＮＰＭＬ）； ＨｘＰＭＬ１＝ｚｅｒｏｓ（ＮＰＭＬ，ＮＰＭＬ）； ＨｙＰＭＬ１＝ｚｅｒｏｓ（ＮＰＭＬ，ＮＰＭＬ）； Ｓｉｇｍａｘ＿ｚ１＝ｒｅｐｍａｔ（Ｓｉｇｍａｘ（２：２：ＮＰＭＬ＊２）＇ ，１，ＮＰＭＬ）；
吴炳坚：硕士研究生 基 金 项 目 ： 江 苏 省 自 然 科 学 基 金 （ＢＫ ２００４０５９）
－２６４－ ３６０元 ／ 年 邮局订阅号： ８２－９４６
图 ３：ＰＭＬ 内 ， 各 区域 的 导 电率 σ与 导 磁率 σ＊ 分 布 示 意 图 ， 在 σ＝０ 及 σ＊≠的区域中， σ、σ＊ 由内层至外层逐渐增加。
图 ７ 取频率 ｗ＝０．５ 时出现的光波图
３ 结论
本文 通 过 运 用 时 域 有 限 差 分 法 及 ＰＭＬ 边 界 条 件 来 仿 真 光 子 晶体 波 导 克服 了 以 往用 平 面 波展 开 法 带来 的 计 算复 杂 、精 确 度 不 高 、占用 计 算 机内 存 空 间大 、计 算 时间 周 期 长等 缺 点， 同 时 借 助 Ｍａｔｌａｂ 工 具 在 数 值 计 算 和 科 学 绘 图 等 方 面 显 示 出 比 其 他 编程工具所难以比拟的优越性， 只要经过很短的时间就能得出 精确的结果。（下转第 ３００ 页）
９ｇ￣１１ｇ。
有关设计信息如下： ①按前文方法确定属性集
＝｛冲击加速度， 振动频率， 振动加速度｝； ②根据式（２），
把 ｚ３ 的属 性 值 处理 为 ，
式 中 常 数 项 θ取
为 ２．０； ｚ１ 的 属 性值 表 达 为模 糊 集
；
③根 据 专 家 经 验 ， 确 定 权 重 集 Ｗ＊＝｛０．４０， ０．３０， ０．３０｝； ④从 知 识 库 中 调 出 有 关 特 征 模 型 的 数 据 ： Ｍ＝｛旋 转 式 ， 永 磁 式 ， 拍 合 式｝
技 ＝
术
创
按照前文所述的方法， 算得
新 ０．８８１＋０．３０×０．９４３ ＝ ０．９１； 同理， 算得
＝ ０．４０×０．９０３＋０．３０×
＝０．６２；
＝
０．５８。因此， 本产品应选择旋转式结构。
为设计方案的可
信度。这里， ｜Ｚ＇ ｜为满足“约定④”前提条件的属性个数； ｜Ｚ｜＝ｎ， 为
属性集的元素个数。｜Ｚ＇ ｜／｜Ｚ｜偏小说明设计信息不完整， ＣＦ 偏
小则综合体现了信息的不确定。此时，Байду номын сангаас应终止“模糊识别”过程， 启
动 ＣＢＲ。
３．３．４ 采用模糊识别方法求解结构设计问题的应用实例
设： 欲 设 计一 种 电 磁继 电 器 ， 要求 是 旋 转式 、永 磁 式或 者 拍
Ｅｚ（１，Ｎｐｙ－ （ＮＭｌａｔ－ １）／２：Ｎｐｙ＋（ＮＭｌａｔ－ １）／２）＝Ｅｚ（１，Ｎｐｙ－ （ＮＭｌａｔ－
１）／２：Ｎｐｙ＋（ＮＭｌａｔ－ １）／２）＋ｓｉｎ（Ｗ＊ｍ＊Ｄｔ）；
％模拟连续平面波源 ｅｘｐ（－ （ｍ＊Ｗ＊Ｄｔ－ ３）＾２）
２．３．２ 数学模拟 ＦＤＴＤ 方程
Ｈｘ（：，２：ＮＴｙ－ １）＝Ｈｘ（：，２：ＮＴｙ－ １）－ Ｄｔ＊（Ｅｚ（：，２：ＮＴｙ－ １）－ Ｅｚ（：，１：ＮＴｙ－
在此禁带内而被禁止传播。
我们通过利用这一性质做成光子晶体波导， 为整个回路信
号的传输提供一个快速通道， 能为将来的超大规模全光或光电
子集成回路提供一个优良的基础物质平台。
时域 有 限差 分（ＦＤＴＤ） 方 法自 Ｙｅｅ （１９６６ 年） 提 出以 来 就 得
到 迅 速 发 展， 其 主 要 思 想 是 以 Ｙｅｅ 元 胞 为 空 间 电 磁 场 离 散 单
合 式 ， 动 作 电 压 １８Ｖ， 额 定 电 压 ２７Ｖ， 输 入 回 路 电 流 大 于 ４Ａ 时
动 作 ， 工 作 环 境 ： 冲 击 加 速 度 为 １０ｇ、２０ｇ、３０ｇ、４０ｇ、５０ｇ 的 可 能
性 分 别 为 ０．２、０．７、１．０、０．６、０．１； 振 动 频 率 １０００Ｈｚ； 振 动 加 速 度
２））／（Ｄｙ＊ｍｕ０）；％我 们 考 虑 的 是 无 电 或 磁 损 耗
Ｈｙ（２：ＮＴｘ－ １，：）＝Ｈｙ（２：ＮＴｘ－ １，：）＋Ｄｔ＊（Ｅｚ（２：ＮＴｘ－ １，：）－ Ｅｚ（１：ＮＴｘ－
２，：））／（Ｄｘ＊ｍｕ０）；
的媒质 ｓ＝０．
Ｅｚ＝Ｅｚ＋Ｄｔ＊（（Ｈｙ（２：ＮＴｘ，１：ＮＴｙ－ １）－ Ｈｙ（１：ＮＴｘ－ １，１：ＮＴｙ－ １））／Ｄｘ－
技 １ 引言
术
光 子 晶 体 最 早 是 在 １９８７ 年 由 Ｅ．Ｙａｂｌｏｎｏｖｉｔｃｈ 与 Ｓ．Ｊｏｈｎ 提
出， 是由不同介电系数的物质周期性排列所组成。和在半导体
创 晶体中电子会形成能带结构、带与带间有能隙的情形类似， 在
光子晶体中传播的光波色散曲线， 因介质介电系数周期性排列
新 的 缘 故 ， 也形 成 带 状结 构 即“光子 禁 带 ”。特定 频 率 的光 因 为 落
仿真技术 文章编号：１００８－ ０５７０（２００７）０１－ １－ ０２６４－ ０２
中 文 核 心 期 刊《 微 计 算 机 信 息 》（ 测 控 自 动 化 ）２００７ 年 第 ２３ 卷 第 １－１ 期
基于 Ｍａ ｔｌａ ｂ 的光子晶体波导仿真研究
ＦＤＴＤ Ａｒｉｔｈ ｍ ｅ ｔｉｃ Ｐｒｏ ｇ ｒａ ｍ ｍ ｉｎ ｇ ｉｎ ｐ ｈ ｏ ｔｏ ｎ ｉｃ ｃｒｙｓ ｔａ ｌｓ ｗ ａ ｖｅ ｇ ｕ ｉｄ ｅ ｏ ｎ ｔｈ ｅ Ｂａ ｓ ｉｓ ｏ ｆ Ｍａ ｔｌａ ｂ Ｌａ ｎ ｇ ｕ ａ ｇ ｅ
《 Ｐ ＬＣ 技术应用 ２００ 例》
邮局订阅号： ８２－９４６ ３６０ 元 ／ 年 － ２６５ －
软件天地
中 文 核 心 期 刊《 微 计 算 机 信 息 》（ 测 控 自 动 化 ）２００７ 年 第 ２３ 卷 第 １－１ 期
补充说明： 上文中“约定④”所述情况， 说明用户提供的设计信息
不完整。为此， 定义ＣＦ＝
图 ４ 二维光子晶体波导结构图
技
图 ５ 当 ＦＤＴＤ 时间模拟 ３８０ 步时图形 我们发现处在光子晶体禁带的光频率只在我们预先设计 好的通道中通过， 而不处于禁带的光频率通过波导时会产生散 射现象， 程序结果较好的验证了波导理论， 为光子晶体集成化， 提供实验依据。
术 创 新
图 ６ 当 ＦＤＴＤ 模拟 ２５００ 步时图形
且：
其中 Ｒ（０） 是入射 角 零 度的 反 射 率， Ｍ 为 一 个 介于 ２ ～４ 的 系 数 ， Ｃ 为 真 空 中 的 光 速 ， Ｌ 为 ＰＭＬ 的 总 层 数 ， Δｘ 为 每 一 层 的 宽度。
ＦＤＴＤ 网格形式的划分， 必须在截断处设置适当的吸收边界条
件， 以便用有限网格空间模拟开放的无限空间或无限长的传输结
元， 将麦克斯韦旋度方程转化为差分方程， 在时间轴和空间轴
上逐步推进地求解， 最终求出空间场的分布。本文主要结合
ＦＤＴＤ 算法边界条件特点， 用 Ｍ ａｔ ｌａｂ 语言进行 编 程 模拟 实 现
光子晶体波导传播。
２ ＦＤＴＤ 算法和边 界处理
２．１ 时域有限差分法
时域有限差分法是将麦克斯韦方程式之中对空间与时间
图 １ Ｙｅｅ 交错晶格与垂直的 Ｘ－ Ｙ 平面
图 ２ ｘ＝ ｉ、ｙ＝ ｊ 的平面晶格 ２．２ 边界条件 在图 １ 中 ， 最边 缘 的 Ｈｘ、Ｈｙ 与 Ｅｚ 由 于 没 有 更 边 缘 的 电 磁 场值可供计算， 故无法确定其值。所以， 如以零为默认值， 则波 传到边缘会好像碰到金属， 波会反弹回解析空间。我们希望波 在未到边缘前， 其电磁场值， 便可衰减为零， 如此便可与边缘的 预设零值相接。此处我们以 Ｂｅｒｅｎｇｅｒ 的 ＰＭＬ 边界理论 （ 如 图 ３ 所示） ， 来讨论程序的写法。
（江苏大学）吴 炳 坚 沈 廷 根
ＷＵ ＢＩＮＧＪＩＡＮ ＳＨＥＮ ＴＩＮＧＧＥＮ
摘要：本 文 首 先 介 绍 了 光 子 晶 体 波 导 的 原 理 ， 然 后 分 析 时 域 有 限 差 分 法 微 分 方 程 及 边 界 条 件 ， 最 后 运 用 ｍａｔｌａｂ 语 言 实 现 二 维
光子晶体波导的仿真。
关键词：光子晶体波导；时域有限差分法；ｍａｔｌａｂ 语言
中 图 分 类 号 ： Ｔ Ｐ ３９１
文献标识码：Ａ
Ａｂｓｔｒ ａｃｔ：ｗｅ ｉｎｔｒｏｄｕｃｅ ｔｈｅ ｔｈｅｏｒｙ ｏｆ ｐｈｏｔｏｎｉｃ ｃｒｙｓｔａｌｓ ｗａｖｅｇｕｉｄｅ， ｔｈｅｎ ａｎａｌｙｚｅ ＦＤＴＤ ａｒｉｔｈｍｅｔｉｃ ａｎｄ ｂｏｕｎｄａｒｙ ｃｏｎｄｉｔ．Ｆｉｎａｌｌｙ， ｗｅ ｕｓｅ Ｍａｔｌａｂ ｌａｎｇｕａｇｅ ｔｏ ｐｒｏｇｒａｍｍｅ ＦＤＴＤ ａｒｉｔｈｍｅｔｉｃ ｉｓ ｃａｒｒｉｅｄ ｏｕｔ，ｔｈｅ ｒｅｓｕｌｔｓ ｐｒｅｓｅｎｔ ａ ｔｈｅｏｒｅｔｉｃａｌ ｂａｓｅ ｆｏｒ ｐｒｏｃｅｓｓｉｎｇ ｐｈｏｔｏｎｉｃ ｄｅｖｉｃｅ． Ｋｅｙ ｗｏｒ ｄｓ：ｐｈｏｔｏｎｉｃ ｃｒ ｙｓｔａｌｓ ｗａｖｅｇｕｉｄｅ，ＦＤＴＤ，Ｍａｔｌａｂ ｌａｎｇｕａｇｅ
构。在网格的划分过程中， 时间步长取为
，
其中 ｃ 为自由空间光速。最后进行波源的设置。 ２．３ Ｍａｔｌａｂ 语言部分源代码 ２．３．１ 初始化设置 Ｗ＝０．３６； ％设 定 光 子 晶 体 波 源 频 率 ， 需 要 注 意 的 是 只 有 在
光子晶体禁带范围的频率才能产生波导
ＭＬａｔｘ＝１１； ％ｘ 方向的晶格数 ＭＬａｔｙ＝１１； ％ｙ 方向的晶格数 ＮＭｌａｔ＝２１； ％每个晶格的网格数 ＮＴｘ＝ＭＬａｔｘ＊ＮＭｌａｔ＋１； ％ｘ 方向的总网格数 ＮＴｙ＝ＭＬａｔｙ＊ＮＭｌａｔ＋１； ％ｙ 方向的总网格数 Ｒ＝０．２； ％晶格半径 ｅａ＝１１．４； ％介电常数 ｍｕ０＝４＊ｐｉ＊１．０ｅ－ ７； ％真空中的电容率 ｅ０＝８．８５＊１ｅ－ １２； ％真空中的介电常数 ｃ＝１／ｓｑｒｔ（ｍｕ０＊ｅ０）； ％光速
的微分 ， 改以差分来代替。在此结构下， 每一个电场分量由磁场
分量围绕。由安培定律可根据周围磁场， 决定中间电场的变化
量。每一个磁场分量由电场分量围绕。由法拉第定律， 可根据周
围电场决定中间磁场的变化量。如图（ １， ２） 所示对于二维横磁
模（ ＴＭ）
，Ｈｚ＝０。我们得到 ＦＤＴＤ 的基本方程：
σ是电导率（ ｓ／ｍ） ，ｓ 是磁电阻率 （ Ω／ｍ） μ和 ε分 别是 媒 质 的 磁导率和电容率。
《 现场总线技术应用 ２００ 例》
您的论文得到两院院士关注
仿真技术
理 论 上 σＸ、σ＊Ｘ 愈 大 ， 则 波 衰 减 愈 快 ， 但 为 解 析 空 间 的 精 准 度， 则交界处波的连续性要良好， 因此 σＸ、σ＊Ｘ 也不能太大。通常 σＸ、σ＊Ｘ 在 ＰＭＬ 内 的 安排 是 采 渐进 式 ， 最 外层 ， σＸ（Ｌ）、σＸ（ｍａｘ）
Ｓｉｇｍａｘ＿ｘ１＝ｒｅｐｍａｔ（Ｓｉｇｍａｘ（２：２：ＮＰＭＬ＊２）＇ ，１，ＮＰＭＬ）；
Ｓｉｇｍａｘ＿ｙ１＝ｒｅｐｍａｔ（Ｓｉｇｍａｘ（１：２：ＮＰＭＬ＊２－ １）＇ ，１，ＮＰＭＬ）； Ｓｉｇｍａｙ＿ｚ１＝ｆｌｉｐｌｒ（ｒｅｐｍａｔ（Ｓｉｇｍａｘ（２：２：ＮＰＭＬ＊２），ＮＰＭＬ，１））； Ｓｉｇｍａｙ＿ｘ１＝ｆｌｉｐｌｒ（ｒｅｐｍａｔ（Ｓｉｇｍａｘ（１：２：ＮＰＭＬ＊２－ １），ＮＰＭＬ，１））； Ｓｉｇｍａｙ＿ｙ１＝ｆｌｉｐｌｒ（ｒｅｐｍａｔ（Ｓｉｇｍａｘ（２：２：ＮＰＭＬ＊２），ＮＰＭＬ，１））；％这 里 只设置一个区域的 ＰＭＬ 层， 其余类同。 ２．３．４ 程序模拟效果图