计算电磁学之FDTD算法的MATLAB语言实现

South China Normal University 课课程程设设计计实实验验报报告告 课程名称 课程名称 计算电磁学计算电磁学 指导老师 指导老师 专业班级 专业班级 20142014 级级 电路与系统电路与系统 姓姓 名 名 学学 号 号 FDTDFDTD 算法的算法的 MATLABMATLAB 语言实现语言实现 摘要 摘要 时域有限差分 FDTD 算法是 K S Yee 于 1966 年提出的直接对麦克斯韦 方程作差分处理 用来解决电磁脉冲在电磁介质中传播和反射问题的算法 其基 本思想是 FDTD 计算域空间节点采用 Yee 元胞的方法 同时电场和磁场节点空间 与时间上都采用交错抽样 把整个计算域划分成包括散射体的总场区以及只有 反射波的散射场区 这两个区域是以连接边界相连接 最外边是采用特殊的吸收 边界 同时在这两个边界之间有个输出边界 用于近 远场转换 在连接边界上采 用连接边界条件加入入射波 从而使得入射波限制在总场区域 在吸收边界上采 用吸收边界条件 尽量消除反射波在吸收边界上的非物理性反射波 本文主要结合 FDTD 算法边界条件特点 在特定的参数设置下 用 MATLAB 语 言进行编程 在二维自由空间 TEz 网格中 实现脉冲平面波 关键词 关键词 FDTD MATLAB 算法 1 1 绪论绪论 1 11 1 课程设计背景与意义课程设计背景与意义 20 世纪 60 年代以来 随着计算机技术的发展 一些电磁场的数值计算方 法逐步发展起来 并得到广泛应用 其中主要有 属于频域技术的有限元法 FEM 矩量法 MM 和单矩法等 属于时域技术方面的时域有限差分法 FDTD 传输线矩阵法 TLM 和时域积分方程法等 其中 FDTD 是一种已经获 得广泛应用并且有很大发展前景的时域数值计算方法 时域有限差分 FDTD 方法于 1966 年由 K S Yee 提出并迅速发展 且获得广泛应用 K S Yee 用后来被称作 Yee 氏网格的空间离散方式 把含时间变量的 Maxwell 旋度方程转化为差分方程 并成功地模拟了电磁脉冲与理想导体作用的时域响 应 但是由于当时理论的不成熟和计算机软硬件条件的限制 该方法并未得到 相应的发展 20 世纪 80 年代中期以后 随着上述两个条件限制的逐步解除 FDTD 便凭借其特有的优势得以迅速发展 它能方便 精确地预测实际工程中 的大量复杂电磁问题 应用范围几乎涉及所有电磁领域 成为电磁工程界和理 论界研究的一个热点 目前 FDTD 日趋成熟 并成为分析大部分实际电磁问 题的首选方法 1 21 2 时域有限差分法的发展与应用时域有限差分法的发展与应用 经过四十多年的发展 FDTD 已发展成为一种成熟的数值计算方法 在发展 过程中 几乎都是围绕几个重要问题展开的 即数值稳定性 计算精度 数值 色散 激励源技术以及开域电磁问题的吸收边界条件等 数值稳定和计算精度 对任何一种数值计算方法都是至关重要的 其中激励源的设计和引入也是 FDTD 的一个重要任务 目前 应用最广泛的激励源引入技术是总场 散射场体系 对 于散射问题 通常在 FDTD 计算空间中引入连接边界 它将整个计算空间划分为 内部的总场区和外部的散射场区 如图 1 1 利用 Huygens 原理 可以在连接 边界处引入入射场 使入射场的加入变得简单易行 图 1 1 总场 散射场区 2 2 FDTDFDTD 法基本原理法基本原理 2 12 1 MaxwellMaxwell 方程和方程和 YeeYee 氏算法氏算法 根据电磁场基本方程组的微分形式 若在无源空间 其空间中的煤质是各 向同性 线性和均匀的 即煤质的参数不随时间变化且各向同性 则 Maxwell 旋度方程可写成 式中 E 是电场强度 单位为伏 米 V m H 是磁场强度 单位为安 米 A m 表示介质介电系数 单位为法拉 米 F m 表示磁导系数 单位为亨利 米 H m 表示介质电导率 单位为西门子 米 S m m表示导磁率 单 位为欧姆 米 m K S Yee 在 1966 年建立了如图 2 1 所示的空间网格 这就是著名的 Yee 氏 元胞网格 图 2 1 Yee 氏网格及其电磁场分量分布 电场和磁场被交叉放置 电场分量位于网格单元每条棱的中心 磁场分量位 于网格单元每个面的中心 每个磁场 电场 分量都有 4 个电场 磁场 分量环绕 这样不仅保证了介质分界面上切向场分量的连续性条件得到自然满足 而且还允 许旋度方程在空间上进行中心差分运算 同时也满足了法拉第电磁感应定律和安 培环路积分定律 也可以很恰当地模拟电磁波的实际传播过程 2 22 2 时域有限差分法相关技术时域有限差分法相关技术 2 2 12 2 1 数值稳定性问题数值稳定性问题 FDTD 方程是一种显式差分方程 在执行时 存在一个重要的问题 即算 法的稳定性问题 这种不稳定性表现为在解显式方程时 随着时间步数的继续 增加 计算结果也将无限制地增加 Taflove 于 1975 年对 Yee 氏差分格式的稳 定性进行了讨论 并导出了对时间步长的限制条件 数值解是否稳定主要取决 于时间步长 t 与空间步长 x y z 的关系 对于非均匀媒质构成的计 算空间选用如下的稳定性条件 上式是空间和时间离散之间应当满足的关系 又称为 Courant 稳定性条件 若 采用均匀立方体网格 x y z s 其中 v 为计算空间中的电磁波的最大速度 2 2 22 2 2 数值色散数值色散 FDTD 方程组是对 Maxwell 旋度方程进行差分近似 在进行数值计算时 将会在计算网格中引起数字波模的色散 即在 FDTD 网格中 电磁波的相速与 频率有关 电磁波的相速度随波长 传播方向及变量离散化的情况不同而改变 这种关系由非物理因素引起 且色散将导致非物理因素引起的脉冲波形畸变 人为的各向异性和虚假折射等现象 为获得理想的色散关系 问题空间分割应 按照小于正常网格的原则进行 一般选取的最大空间步长为 max min 20 min为所研究范围内电磁波的最小波长 由上分析说明 数值 色散在用 FDTD 法分析电磁场传播中的影响是不可能避免的 但我们可以尽可 能的减小数值色散的影响 2 2 32 2 3 离散网格的确定离散网格的确定 无论是简单目标还是复杂目标 在进行 FDTD 离散时网格尺寸的确定 除了受计算资源的限制不可能取得很小外 还需要考虑以下几个因素 1 目标离散精确度的要求 网格应当足够小以便能精确模拟目标几何形状 和电磁参数 2 FDTD 方法本身的要求 主要是考虑色散误差的影响 设网格为立方体 x x y z 所关心频段的频率上限为 f max 对应波长为 f min 则考虑 FDTD 的数值色散要求 min N 通常 N 10 上式是根据已知所关心频率上限情况下来确定 FDTD 网格尺寸 的 反之 若给定 则 FDTD 计算结果可用的上限频率也随之确定 2 32 3 吸收边界条件吸收边界条件 由时域有限差分法的基本原理可知 在利用时域有限差分法研究电磁场时 需在全部问题空间建立 Yee 氏网格空间 并存储每个单元网格上任一时间步的 六个场分量用于下一时间步的计算 而在对于辐射 散射这类开放系统的实际 研究中 不可能有无限大的存储空间 因此 必须在某处将网格空间截断 且 在截断边界网格点处运用特殊的场分量计算方法 使得向边界面行进的波在边 界处保持 外向行进 特性 无明显的反射现象 并且不会使内部空间的场产 生畸变 从而用有限网格空间模拟电磁波在无界空间中传播的情况 具有这种 功能的边界条件称之为吸收边界条件 或辐射边界条件 或网格截断条件 如 图 2 2 所示 图 2 2 附加截断边界使计算区域变为有限域 从 FDTD 的基本差分方程组可以看出 在截断边界面上切向场分量的计算需 要利用计算空间以外的电磁场分量 因此 FDTD 基本差分方程对这些截断边界面 上的场分量失效 如何处理截断边界上的场分量 使之与需要考虑的无限空间 有尽量小的差异 是 FDTD 中必须很好解决的一个重要问题 实际上 这是要求 在误差可容忍的范围内 计算空间中的外向波能够顺利通过截断边界面而不引 起波的明显反射 使有限计算空间的数值模拟与实际情况趋于一致 对外向波 而言 就像在无限大空间中传播一样 所以 需要一种截断边界网格处的特殊 计算方法 它不仅要保证边界场计算的必要精度 而且还要大大消除非物理因 素引起的波反射 使得用有限的网格空间就能模拟电磁波在无限空间中的传播 但是如果处理不当 截断边界面可能造成较大反射 构成数值模拟误差的一部 分 甚至可能造成算法不稳定 加于截断边界场分量符合上述要求的算法就称为吸收边界条件 Absorbing Boundary Conditions 2 42 4 FDTDFDTD 计算所需时间步的估计计算所需时间步的估计 为了使计算达到稳定 通常计算所需要时间步按照电磁波往返穿越 FDTD 计 算区对角线 3 5 次来估计 若 FDTD 计算区总元胞数为 N 则对角线上元胞约 为 三维 按照 Courant 稳定条件 设计算中心区 t 2c2c 即穿 越对角线一次需要时间步为 总计算时间步约为 需 对于二维 情况则约为 或者说 计算时间步大约等于 FDTD 计算区对角线 上元胞数目的 12 20 倍 实际上 计算所需时间步还与散射体具体形状 结构 有关 图 2 3 给出了应用 FDTD 分析电磁场问题时的程序流程图 图 2 3 FDTD 程序流程图 3 3 课程设计要求及课程设计要求及 MATLABMATLAB 语言实现语言实现 3 13 1 课程设计要求课程设计要求 采用总场 散射场技术 在二维自由空间 TEz 网格中 实现脉冲平面波 总场区域 100 100 个网格 周边散射场区域 20 个网格 散射场区域外采用 PML 吸收边界条件 平面波 0 度角入射 使用高斯脉冲 中心频率为 1GHz 带宽为 1 e 采用正方形 Yee 网格 离散步长 1 5cm 时间离散步长为 2c 3 23 2 MATLABMATLAB 程序程序 此程序实现自由空间的平面波 此程序采用总场 散射场技术 此程序使用PML设置吸收边界条件 KE 128 真空区域网格数 IE KE JE KE T 0 NSTEP 300 迭代次数 e 2 71828 初始化 dz zeros KE KE ez zeros KE KE hx zeros KE KE hy zeros KE KE ihx zeros KE KE ihy zeros KE KE ga ones KE KE ez inc zeros 1 JE hx inc zeros 1 JE PML 参数 gi2 ones 1 IE gi3 gi2 fi1 zeros 1 IE fi2 ones 1 IE fi3 fi2 gj2 ones 1 JE gj3 gj2 fj1 zeros 1 JE fj2 ones 1 JE fj3 fj2 npml 8 for i 1 8 设置PML层中的参数 xnum npml i xd npml xxn xnum xd xn 0 33 xxn 3 gi2 i 1 0 1 0 xn gi2 IE 1 i 1 0 1 0 xn gi3 i 1 0 xn 1 0 xn gi3 IE 1 i 1 0 xn 1 0 xn xxn xnum 0 5 xd xn 0 25 xxn 3 fi1 i xn fi1 IE 2 i xn fi2 i 1 0 1 0 xn fi2 IE 2 i 1 0 1 0 xn fi3 i 1 0 xn 1 0 xn fi3 IE 2 i 1 0 xn 1 0 xn end for j 1 8 xnum npml j xd npml xxn xnum xd xn 0 33 xxn 3 gj2 j 1 0 1 0 xn gj2 JE 1 j 1 0 1 0 xn gj3 j 1 0 xn 1 0 xn gj3 JE 1 j 1 0 xn 1 0 xn xxn xnum 0 5 xd xn 0 25 xxn 3 fj1 j xn fj1 JE 2 j xn fj2 j 1 0 1 0 xn fj2 JE 2 j 1 0 1 0 xn fj3 j 1 xn 1 0 xn fj3 JE 2 j 1 0 xn 1 0 xn end 总场 散射场边界 ia 28 ib IE ia 1 ja 28 jb JE ja 1 ez inc low m1 0 ez inc low m2 0 ez inc high m1 0 ez inc high m2 0 信号源 t0 80 脉冲高度 spread 12 脉冲宽度 网格 ddx 0 015 离散步长 dt ddx 2 3e8 计算时间离散步长 epsz 8 851 1e 12 自由空间的介电常数 迭代求解电场与磁场 for n 1 NSTEP a 1 T T 1 for j 2 JE ez inc j ez inc j 0 5 hx inc j 1 hx inc j end 入射波缓冲区 ez inc 1 ez inc low m2 ez inc low m2 ez inc low m1 ez inc low m1 ez inc 2 ez inc JE 1 ez inc high m2 ez inc high m2 ez inc high m1 ez inc high m1 ez inc JE 2 计算 dx 区域 for i 2 KE for j 2 KE dz i j gi3 i gj3 j dz i j gi2 i gj2 j 0 5 hy i j hy i 1 j hx i j hx i j 1 end end 输入高斯脉冲信号源 pulse exp 0 5 t0 T 2 spread 2 ez inc 3 pulse 入射波 DZ 参数 for i ia ib dz i ja dz i ja 0 5 hx inc ja 1 dz i jb dz i jb 0 5 hx inc jb end 计算电场dx区域 for i 1 KE for j 1 KE ez i j ga i j dz i j end end 设置电场 EZ 边缘为0 作为PML的参数 for j 1 JE 1 ez 1 j 0 ez IE j 0 end for i 1 IE 1 ez i 1 0 ez i JE 0 end for j 1 JE 1 hx inc j hx inc j 0 5 ez inc j ez inc j 1 end 计算磁场hx区域 for i 1 KE for j 1 KE 1 curl e ez i j ez i j 1 ihx i j ihx i j fi1 i curl e hx i j fj3 j hx i j fj2 j 0 5 curl e ihx i j end end 入射磁场Hx区域参数 for i ia ib hx i ja 1 hx i ja 1 0 5 ez inc ja hx i jb hx i jb 0 5 ez inc jb end 计算磁场hy区域 for i 1 KE 1 for j 1 KE curl e ez i 1 j ez i j ihy i j ihy i j fj1 j curl e hy i j fi3 i hy i j fi2 i 0 5 curl e ihy i j end end 入射磁场hy区域参数 for j ja jb hy ia 1 j hy ia 1 j 0 5 ez inc j hy ib j hy ib j 0 5 ez inc j end 图形展示 m moviein 500 x 1 KE y 1 KE X Y meshgrid x y surf X Y ez axis 0 KE 0 KE 1 1 set gcf Color white Number off Name sprintf Simulation FDTD 2D Iteration i n title sprintf t