时域有限差分法的Matlab仿真

张通　孙晶

【摘 要】文章介绍了时域有限差分法的基本原理，利用matlab仿真，实现了用时域有限差分程序来计算二维问题空间中的电场分布。

【关键词】时域有限差分法；Matlab；电场分布

Simulation of Finite Difference Time Domain Method Using Matlab

ZHANG Tong SUN Jing

（College of Physics，Mechanical and Electrical Engineering，Jishou University， JiShou Hunan 416000，China）

【Abstract】The basic principle of finite difference time domain is introduced in this paper.With two-dimensional finite difference time domain program to calculate the electric field distribution of the problem space is implemented using Matlab.

【Key words】FDTD；Matlab；Electric field distribution

0 引言

時域有限差分（Finite Difference Time Domain，FDTD）法是K.S.Yee在1966年给出的利用有限差分式把麦克斯韦（Maxwell）旋度方程替换为一组差分方程[1]，并提供所解问题中电磁特性物理意义的算法，可直接在时域中求解。

Matlab是一种功能强大、高效的高级技术计算语言和交互式环境[2]，在科学和工程领域中赢得了极为广泛的应用，将其用于FDTD法的数值计算及仿真，不仅可以简化程序设计、操作方便，另外运算结果也更简洁。

因此，本文将结合Matlab强大的数组运算和绘图功能，通过对FDTD法编程来模拟出二维问题空间中电场分布。

1 FDTD法的基本原理

FDTD算法将问题空间离散为电场和磁场分量在其位置上交叉放置的空间网格点，并以中心差分的方式近似Maxwell方程中关于空间和时间的导数，通过时间向前推进的差分方程模拟出电磁场在时域的进程。空间网格中，电场分量位于Yee元胞网格单元每条棱的中心，磁场分量位于网格单元每个面的中心[3]，如图1所示。

1.1 Maxwell方程的差分形式

Maxwell旋度方程为：

？荦×H=+J；？荦×E=--Jm（1）

已知本构关系表达式为：

D=？着E；B=？滋H；J=？啄E；Jm=？啄mH

在直角坐标系中，根据本构关系把（1）式写为：

（2）

下面我们求解（2）式的中心差分，令f（x，y，z，t）表示E或H某一分量，离散形式写为：

f（x，y，z，t）=f（i？驻x，j？驻y，k？驻z，n？驻t）=fn（i，j，k）（3）

在二维问题空间中，假定任意的电磁场分量只与x，y坐标有关，与z坐标无关，即？坠/？坠z=0，以TE波为例，Hx=Hy=Ez=0，由（2）式可得

（4）

用中心差分式来近似（4）式中的导数，根据场分量的位置，并采取？驻x=？驻y=？驻z=？啄离散方式，得到了关于TE波的FDTD公式为：

Ex=CA·Ex+CB'·Hz-Hz（5）

Ey=CA·Ey+CB'·Hz-Hz（6）

（7）

式中系数CA，CB'，CQ'的定义为：

CA=1-/1+；CB'=/1+

CP=1-/1+；CQ'=/1+

为了统一TE波、TM波两者方程的离散形式，分别将（5）、（6）、（7）式中的空间位置标号移动1/2，时间移动？驻t/2，以上式子分别写为：

Ex=CA·Ex+CB'·Hz-Hz（8）

Ey=CA·Ey+CB'·Hz-Hz（9）

（10）

利用TE波与TM波之间的对偶关系，写出通用于计算求解二维问题空间中TE波与TM波的FDTD程序。

1.2 数值色散及稳定性条件

FDTD方法为场的行为提供了一种解，连续函数的导数有限差分近似给解引入了误差，我们把用FDTD数值方法得到的相速与实际的相速之间的差别称为数值色散。为保证结果准确性，空间网格大小应满足？姿min≥10？驻，？驻=min（？驻x，？驻y，？驻z），？姿min表示媒质空间中最小波长值，减小网格大小虽然会减小数值色散，但在计算中将会占用更多的内存。

为保证数值计算稳定性，根据Cournant稳定条件，算法中的时间步长应满足：c？驻t≤

其中c=1/为介质中的光速，一般选取？驻t=δ/（2c）。

2 问题空间中电场分布的模拟

以图2所示二维问题为例，几何图形中包含半径为0.2m、介电常数为4的圆柱，激励信号频率为1GHz。此问题空间由边长为5mm的正方形网格构成，端接为8mm的PML，圆柱与PML边界问题的空气隙在xn、yn和yp方向上为30个网格，在xp方向上为80个网格，其激励源为一正弦波形的外加电流密度。

其中，J1为激励线源，E1为电场取样点。此例中，定义了两种输出类型，包括瞬态电场分布和某频率下的电场分布。通过Matlab编程来实现问题的定义与模拟，其主要编程如下：

（1）定义几何体。设置圆柱的位置范围及材料类型，包括中心坐标、半径及电介质。

（2）定义激励源。设置激励源波形及参量，包括初始值矢量坐标、波形类型及量级。

（3）定义并初始化输出参量。设置频率边界，定义并初始化取样电场、磁场及瞬态电场、某频率下的电场，取更新频率为10个时间步。

（4）显示取样参数，计算并获取在节点上的电场。

（5）获取取样电场给定频率的频域响应，规定此程序中场量在6000个时间步以后获取。

（6）显示频率边界输出，计算得到的响应曲线。

使用Matlab运行仿真后，仿真结果如下图所示：

假定正弦激励的时域响应在6000时间步达到稳定后获得稳态场的幅度，所以用给定的程序能够获得场频域响应的幅度，而且场的幅度响应仅属于此单一的激励频率。由图5发现，某些频率的取样电场幅度和相位响应并不能同时求解，只能显示1GHz单一频率幅度。为解决此问题，我们将采用更有效地离散傅里叶变换（DFT）[4]。

如前面所讨论的，我们将求解多频率下电场幅度，如果激励信号是具有一定频谱宽度的波形，那么使用DFT就可以得到多个频率的结果。首先，对于激励波形，其频谱应包含所想要的频率分量。其次，为计算多个频率的场分布，必须实施实时DFT。使用DFT技术同时求解出多频率取样电场的幅度和相位响应如图8所示。另外，如图9所示求解出的在1GHz频率下的电场幅度分布与图6结果相同，同时也求解出在2GHz频率下的电场幅度分布，如图10所示，验证了此方法的可行性。