时域有限元方法

2024-07-02

时域有限元方法（精选7篇）

时域有限元方法篇1

时域有限差分方法 (FDTD) 1966年由Yee提出后, 经过许多年的发展, 在电磁波散射方面得到了广泛的应用。现在, 在光子晶体能隙的计算方面也得到了很多应用。笔者在一维、二维的FDTD编程计算中, 发现时间步长的选择不同, 会导致结果的不同, 下面进行讨论。

一、一维FDTD差分方程表达形式

一维波动方程, 在一维无限长介质中, 其解为u=f (x-vt) , 是以速度v向x轴正向传播的行波。用差分形式, 把波动方程进行改写:

二、初始、边界条件的确定以及编程计算的结果

边界条件:u (t, 1) =0, u (t, Nx) =0, 对t=1到Nt。

初始条件:由于在差分形式的递推公式 (1) 中, 要计算t=3的值, 需要知道t=1, 2时的值。本文规定初始条件为u (t, i) =0, 对t=1, 2.

下面是两组不同的参数, 用Matlab编程计算的结果:

第1组:Lx=1×10-4, Δx=1×10-6, Nx=100, Nt=100, v=c=3×108,

信号源, 观察点在i=80, 也就是x=80Δx处。

第2组:Lx=1×10-4, Δx=1×10-6, Nx=100, Nt=1000, v=c=3×108,

信号源, 观察点在i=80, 也就是x=80Δx处。参数中时间长度Nt与第1组不同, 其他参数相同。

第3组:Lx=1×10-4, Δx=1×10-7, Nx=1000, Nt=3000, v=c=3×108,

信号源, 观察点在i=500, 也就是x=500Δx处。

从图中可以看出:

t=1到600, 在观察点无信号, 信号还没有到达;

t=600到900, 是过渡阶段, 与数值化的计算有关, 在观察点信号尚不稳定;

t=900到2400, 观察点信号与用连续方程得到的结果一致;

t=2400以后, 观察点的信号的幅度比源的信号还大, 是因为文中取一维两端边界为零, 反射波到达观察点, 与正向波形成干涉。

用快速傅里叶变换得出的频谱透射率, 由于已经包含了开始的过渡阶段和以后的反射波的干涉的影响, 在0.8左右, 与波的无衰减传播结果接近。

三、结论

在一维的FDTD计算中, 直接使用整数时刻和整数空间坐标, 而不需要使用半整数时间和半整数空间坐标, 简化了坐标的标记形式;时间步长取, 时间总长Nt的取值, 要求在这个时间信号源的信号可以到达观察点, 但又没有界面反射波到达观察点, 这样的参数比较理想, 可以得到与连续方程相同的结果。

摘要：在用一维时域有限差分方法计算波的传播时, 为了防止出现发散的结果, 要求时间步长Δt<Δx/v, 但具体取多大的值, 没有定论。通过编程试验, 时间步长取Δt<0.75Δx/c, 而时间步数Nt大于波从波源传到观察点所需要的时间, 又小于反射波到达观察点的时间。而且计算的时间格点和位置格点全部都可以取整数。

关键词：FDTD,有限时域差分方法,参数

参考文献

[1]葛德彪, 阎玉波.电磁波时域有限差分方法[M].第2版.西安:西安电子科技大学出版社, 2005.

[2]曾辉, 杨亚培.FDTD法与平面波展开法在光子禁带计算中的差异分析[J].电子科技大学学报, 2005, 34 (6) :901-904.

[3]张国华, 袁乃昌, 付云起.FDTD方法分析光子带隙能带结构[J].微波学报, 2001, 17 (4) :14-17.

[4]张亮, 寇晓艳.FDTD的Matlab语言的实现[J].延安大学学报:自然科学版, 2009, 28 (2) :57-59.

[5]朱章虎, 卢万铮, 冯奎胜.FDTD计算中PML的简化应用及编程实现[J].空军工程大学学报 (自然科学版) , 2006, 7 (2) :55-57.

时域有限元方法篇2

时域有限差分方法[finite-difference time-domain (FDTD) method]自1966年由Yee[1]提出来后, 一直是计算电磁学常用的方法之一, 并广泛应用于电磁散射、天线的分析与设计、雷达截口的计算等电子工业与国防工业。但是, 由于FDTD方法是条件稳定的, 即在二维情形下, 时间和空间步长分别为Δt, Δx, Δy;必须满足Courant-Friedrichs-Lewy (CFL) 稳定条件: cΔt≤[1/ (Δx) 2+1/ (Δy) 2]-1/2。其中c是在介质中的光速, ε, μ是介电常数和磁导率。为了解决上述问题, 由文献[2,3]提出了交替方向隐式时域有限差分方法 (ADI-FDTD) 方法, 并用Fourier方法证明了这种格式是无条件稳定的。近年来[4], 将分裂算子方法与FDTD方法结合提出新颖而且简单的分裂算子的时域有限差分方法 (S-FDTD) , 与二维ADI-FDTD相比都是无条件稳定的二阶格式, 但是S-FDTD格式采用算子分裂技术, 所以格式简单, 计算时间短, 并且在模拟一种散射问题时, S-FDTD比ADI-FDTD更精确。麦克斯韦方程的高阶时域有限差分方法已有较多的研究工作[4,5,6,7,8,9], 但是这些方法都没有使用分裂算子技术。

本文将高阶差分与S-FDTD方法相结合提出一种新的方法, 即在分裂后的方程基础上对空间采取四阶中心差分, 从而得到分裂 (splitting) 的高阶 (high-order) 的时域有限差分 (FDTD) 格式SHO-FDTDⅠ;并在此格式的基础上通过增加扰动项减小分裂误差给出它的修正格式SHO-FDTDⅡ, 给出了具体的计算步骤, 然后用Fourier方法证明了这两种格式是无条件稳定的, 并给出这两种格式的数值弥散关系式和数值弥散误差的计算。与条件稳定的高阶FDTD相比, 新方法是无条件稳定的;与高阶的ADI-FDTD相比, 新格式更简单, 易于编程和实现。

1 二维麦克斯韦方程的SHO-FDTD格式

考虑以下二维横向电波:

$\frac{\partial E_{x}}{\partial t} = \frac{1}{ε} \frac{\partial Η_{z}}{\partial y}$ (1.1)

$\frac{\partial E_{y}}{\partial t} = - \frac{1}{ε} \frac{\partial Η_{z}}{\partial x}$ (1.2)

$\frac{\partial Η_{z}}{\partial z} = \frac{1}{μ} (\frac{\partial E_{x}}{\partial y} - \frac{\partial E_{y}}{\partial x})$ (1.3)

其中 $\vec{E} = (E_{x} (x, y, t), E_{y} (x, y, t))$ 表示电场, $Η_{z} = Η_{z} (x, y, t)$ 表示磁场, ε和μ分别是介电常数和磁导率, Ω=[0, a]×[0, b], t ∈ (0, T]。其边界满足理想导体边界条件: $E_{x} (x, 0, t) = E_{x} (x, b, t) = E_{y} (0, y, t) = E_{y} (a, y, t) = 0$ 。初始条件: Ex0=Ex0 (x, y) , Ey0=Ey0 (x, y) , Hz0=Hz0 (x, y) 。

为了书写的方便, 我们仅讨论常系数的情形, 这里所用的方法容易推广到变系数的情形。对空间区域Ψ以及时间区域[0, T]采取如下剖分:

其中Δx和Δy分别是沿x轴方向和沿y轴方向的空间离散步长, Δt是时间步长, I, J, N为整数。

对于任意一个网格函数F (t, x, y) , 引入以下记号:

F $_{α, β}^{m}$ =F (mΔt, αΔx, βΔy) , $δ_{t} F_{α, β}^{m} = \frac{F_{α, β}^{m + \frac{1}{2}} - F_{α, β}^{m - \frac{1}{2}}}{Δ t},$

$δ_{x}^{2} F_{α, β}^{m} = \frac{1}{Δ x} (\frac{1}{24} F_{α - \frac{3}{2}, β}^{m} - \frac{9}{8} F_{α - \frac{1}{2}, β}^{m} + \frac{9}{8} F_{α + \frac{1}{2}, β}^{m} - \frac{1}{24} F_{α + \frac{3}{2}, β}^{m})$ ,

$δ_{y}^{2} F_{α, β}^{m} = \frac{1}{Δ y} (\frac{1}{24} F_{α, β - \frac{3}{2}}^{m} - \frac{9}{8} F_{α, β - \frac{1}{2}}^{m} + \frac{9}{8} F_{α, β + \frac{1}{2}}^{m} - \frac{1}{24} F_{α, β + \frac{3}{2}}^{m})$ 。

用E $_{v_{α, β}}^{m}$ (v=x, y) 表示电场Ev (tm, xα, yβ) 的近似值。H $_{z_{α, β}}^{m}$ 表示磁场Hz (tm, xα, yβ) 的近似值。将分裂算子的时域有限差分方法[9]与高阶差分方法[4]相结合, 提出如下格式:

第一步:

$\frac{E_{y_{i, j + \frac{1}{2}}}^{n + 1} - E_{y_{i, j + \frac{1}{2}}}^{n}}{Δ t} = \frac{1}{2 ε} δ_{x}^{2} = (Η_{z_{i, j + \frac{1}{2}}}^{n + \frac{1}{2}} + Η_{z_{i, j + \frac{1}{2}}}^{n})$ (1.4a)

$\frac{Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}} - Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n}}{Δ t} = - \frac{1}{2 μ} δ_{x}^{2} (E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n})$ (1.4b)

第二步:

$\frac{E_{x_{i + \frac{1}{2}, j}}^{n + 1} - E_{x_{i + \frac{1}{2}, j}}^{n}}{Δ t} = \frac{1}{2 ε} δ_{x}^{2} = (Η_{z_{i + \frac{1}{2}, j}}^{n + \frac{1}{2}} + Η_{z_{i + \frac{1}{2}, j}}^{n})$ (1.5a)

$\frac{Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} - Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}}}{Δ t} = \frac{1}{2 μ} δ_{y}^{2} (E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n})$ (1.5b)

此格式称为分裂算子 (splitting) 的高阶 (high-order) 时域有限差分方法, 记为SHO-FDTDⅠ。其中, 边界条件为: $E_{x_{i + \frac{1}{2}, 0}}^{n} = E_{x_{i + \frac{1}{2}, j}}^{n} = E_{y_{0, j + \frac{1}{2}}}^{n} = E_{y_{i, j + \frac{1}{2}}}^{n} = 0$ 。初始值:

求解步骤为:先解第一步, 由式 (1.4b) 得:

$Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}} = Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n} - \frac{Δ t}{2 μ} δ_{x}^{2} (E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n})$ (1.6)

将式 (1.6) 代入式 (1.4a) 整理得:

显然这是一个七对角的方程组, 系数矩阵元素都是常数而且在每一个时间层上都不变。因此可由一些求解线性方程组的方法如:超松弛迭代法, 共轭梯度法等, 求出 ${E_{y_{i, j + \frac{1}{2}}}^{n + 1}}$ 。然后, 将 $E_{y_{i, j + \frac{1}{2}}}^{n + 1}$ 代入式 (1.6) 显式求出 $Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}}$ 。式 (1.5) 的计算与第一步相似, 但要用到第一步算出的 $Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}}$ 。方程 (1.7) 中的第一个方程和最后一个方程的系数需要调整, 可根据文献[5]中的方法进行, 这里略去。

检查格式SHO-FDTDⅠ的截断误差发现关于时间它的精度不高。为了提高精度, 引入一个扰动项, 得到修正格式SHO-FDTDⅡ如下:

第一步:

第二步与SHO-FDTDⅠ中的第二步相同.

这种格式的初边值条件与式 (1.4) 、式 (1.5) 的初边值条件完全相同, 求解方法也与前一格式相同。为了了解格式SHO-FDTDⅡ的逼近精度, 由式 (1.5a) 、式 (1.5b) 、式 (1.8a) 、式 (1.8b) 消去中间项 $Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}}$ , 可得到SHO-FDTDⅡ的等价格式:

$δ_{t} E_{x_{i + \frac{1}{2}, j}}^{n + \frac{1}{2}} = \frac{1}{2 ε} δ_{y}^{2} (Η_{z_{i + \frac{1}{2}, j}}^{n + 1} + Η_{z_{i + \frac{1}{2}, j}}^{n}) (1.9 a)$

$\begin{array}{l} δ_{t} E_{y_{i, j + \frac{1}{2}}}^{n + \frac{1}{2}} = - \frac{1}{2 ε} δ_{x}^{2} (Η_{z_{i, j + \frac{1}{2}}}^{n + \frac{1}{2}} + Η_{z_{i, j + \frac{1}{2}}}^{n}) + \\ \frac{Δ t}{4 μ ε} δ_{x}^{2} δ_{y}^{2} (E_{x_{i, j + \frac{1}{2}}}^{n + 1} - E_{x_{i, j + \frac{1}{2}}}^{n}) (1.9 b) \end{array}$

$\begin{array}{l} δ_{t} Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + \frac{1}{2}} = \frac{1}{2 μ} [δ_{y}^{2} (E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n}) - \\ δ_{x}^{2} (E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n})] (1.10) \end{array}$

格式SHO-FDTDⅠ的等价格式与式 (1.9a) 、式 (1.9b) 、式 (1.10) 相似, 只要将式 (1.9) 最右端后一项中的“-”改为“+”。由此可以看出, SHO-FDTDⅡ的摄动误差 (二阶) 比SHO-FDTDⅠ的摄动误差 (二阶) 要高一阶, SHO-FDTDⅡ是关于时间2阶, 关于空间是4阶的式 (2.4) 格式, SHO-FDTDⅠ则是式 (1.4) 格式。

2 稳定性分析与数值弥散分析

用Fourier方法分析这两种格式的稳定性。假定格式的差分解具有下列形式。

$\vec{E}_{α, β}^{n} = \vec{E}_{0} ξ^{n} e^{- i (k_{x} α Δ x + k_{y} β Δ y)}$ ;

H $_{z_{α, β}}^{n}$ =Hz0ξne-i (kxαΔx+kyβΔy) 。

其中 $i = \sqrt{- 1}$ , $\vec{E}_{0} = (E_{x_{0}}, E_{y_{0}})^{Τ}$ 和Hz0是振幅, $\vec{k} = (k_{x}, k_{y})$ 是波矢量, $| k | = \sqrt{k_{x}^{2} + k_{y}^{2}}$ , ξ是增长因子。下面将E $_{α, β}^{n}$ , H $_{z_{α, β}}^{n}$ 代入式 (1.5a) 、式 (1.9) 、式 (1.10) , 化简得:

$(ξ - 1) E_{x_{0}} - i \frac{Δ t}{ε} (ξ + 1) b_{y} Η_{z_{0}} = 0$ (2.1)

$\begin{array}{l} (ξ - 1) E_{y_{0}} + \frac{(Δ t)^{2}}{ε μ} (ξ - 1) a_{x} b_{y} E_{x_{0}} + \\ i \frac{Δ t}{ε} (ξ + 1) a_{x} Η_{z_{0}} = 0 (2.2) \end{array}$

$\begin{array}{l} (ξ - 1) Η_{z_{0}} - i \frac{Δ t}{μ} (ξ + 1) b_{y} E_{x_{0}} + \\ i \frac{Δ t}{μ} (ξ + 1) a_{x} E_{y_{0}} = 0 (2.3) \end{array}$

其中 $a_{x} = \frac{(\frac{1}{12} \sin \frac{3}{2} k_{x} Δ x - \frac{9}{4} \sin \frac{1}{2} k_{x} Δ x)}{2 Δ x}, b_{y} = \frac{(\frac{1}{12} \sin \frac{3}{2} k_{y} Δ y - \frac{9}{4} \sin \frac{1}{2} k_{y} Δ y)}{2 Δ y}$ 。由于向量 (Ex0, Ey0, Hz0) 是非零向量, 所以关于Ex0, Ey0, Hz0的方程组的系数矩阵的行列式的值为零。整理得:

(ξ-1) (d0ξ2+2d1ξ+d0) 2=0 (2.4)

式 (2.4) 中

$d_{0} = 1 + \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2};$

$d_{1} = - 1 + \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2}$ 。

解方程式 (2.4) 得: $ξ_{1} = 1, ξ_{2} = d_{0}^{- 1} (- d_{1} + i \sqrt{d_{0}^{2} - d_{1}^{2}}) 。 ξ_{3} = d_{0}^{- 1} (- d_{1} - i \sqrt{d_{0}^{2} - d_{1}^{2}})$ 。显然方程三个根的模都是1, 所以格式SHO-FDTDⅡ是无条件稳定的。对于格式SHO-FDTDⅠ, 同理可得:

c3ξ3+c2ξ2+c1ξ+c0=0 (2.5)

式 (2.5) 中

$\begin{array}{l} c_{3} = 1 + \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2}, \\ c_{2} = - 3 + \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + 3 \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2}; \end{array}$

$\begin{array}{l} c_{1} = 3 - \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + 3 \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2}, \\ c_{0} = - 1 - \frac{(Δ t)^{2}}{ε μ} (a_{x}^{2} + b_{y}^{2}) + \frac{(Δ t)^{4}}{ε^{2} μ^{2}} (a_{x} b_{y})^{2} 。 \end{array}$

方程 (2.5) 的根的表达式非常复杂 (为了简单, 这里省略) , 通过对根的表达式的研究发现 $| ξ | = 1 + Ο (Δ t)$ , 因此格式SHO-FDTDⅠ是耗散的无条件稳定的。为了直观的感受 $| ξ |$ , 下面给出它在不同情况下的变化曲线图。设kx=kcos φ。ky=ksin φ。其中k, φ。是向量k的圆柱坐标, 再设 $Ν_{λ} = \frac{λ}{h}, ω = c k, ω$ 是频率, λ是波长, Δx=Δy=h, 则Nλ是一个波长内的节点个数。

令 $S = \frac{c Δ t}{h}$ , 由CFL的定义可知

$c Δ t \sqrt{\frac{1}{Δ x^{2}} + \frac{1}{Δ y^{2}}} = \frac{c Δ t}{h} \sqrt{2} = \sqrt{2} S$ 。因此, 不妨称S为CFL数。由方程 (2.4) 或式 (2.5) 及其系数的表达式可知, 增长因子ξ是关于S, φ, Nλ的函数。例如, 可推出 $\sin \frac{3}{2} k_{x} Δ x = \sin (\frac{3 π}{Ν_{λ}} \cos φ)$ 。方程 (2.5) 的根的表达式比较复杂。下面我们用Matlab求出它的近似解, 并画出不同情况下它的根的变化曲线。

图1是在S=0.35, Nλ=10情况下给出了根 $| ξ |$ 随角度φ的变化曲线。上边的曲线是方程复数根的模, 下边的曲线是实数根的绝对值。很容易看出复数根的模比1大, 但是 $| ξ | = 1 + Ο (Δ t)$ , 其中, $Δ t = S \times \frac{1}{Ν_{λ}} = 0.035$ 。从图1中显然可以看出曲线位于: y=1±0.035两条直线之间, 这就表明SHO-FDTDⅠ格式是无条件稳定的。

图2和图3分别是在S=1.5, φ=35°和Nλ=10, φ=65°情况下, 方程 (2.5) 的根的模分别随着Nλ和S的变化情况, 从这两个图中同样可以看出SHO-FDTD格式是无条件稳定的。

2.2 数值弥散关系

假设麦克斯韦方程的差分解为: $\vec{E}_{α, β}^{n} =$ $\vec{E}_{0} e^{i (k_{x} α Δ x + k_{y} β Δ y - ω n Δ t)}$ , H $_{z_{α, β}}^{n}$ =Hz0ei (kxαΔx+kyβΔy-ωnΔt) , 其中, ω, kx, ky, $\vec{E}_{0}$ 和Hz0与2.1节的符号相同。将E $_{α, β}^{n}$ , H $_{z_{α, β}}^{n}$ 代入SH-FDTDⅡ的等价格式 (1.5a) 、式 (1.9) 、式 (1.10) , 类似于2.1节对增长因子ξ的推导, 可得:

$\begin{array}{l} \sin^{2} (\frac{1}{2} ω Δ t) = (c Δ t)^{2} [a_{x}^{2} + b_{y}^{2} + (c Δ t)^{2} a_{x}^{2} b_{y}^{2}] \\ \cos^{2} (\frac{1}{2} ω Δ t) (2.6) \end{array}$

式 (2.6) 中 $c^{2} = \frac{1}{ε μ}$ , 式 (2.6) 就是SHO-FDTDⅡ的数值弥散关系式。

同理可得SHO-FDTDⅠ格式的数值弥散关系式:

$\begin{array}{l} \sin^{2} (\frac{1}{2} ω Δ t) = (c Δ t)^{2} [a_{x}^{2} + b_{y}^{2} + (c Δ t)^{2} a_{x}^{2} b_{y}^{2} \times \\ \cos (\frac{1}{2} ω Δ t) \sin^{- 1} (\frac{1}{2} ω Δ t)] \times \\ \cos^{2} (\frac{1}{2} ω Δ t) (2.7) \end{array}$

由于 $\lim_{x \to 0} a_{x} = \lim_{x \to 0} a_{x}^{2} = \frac{k_{x}^{2}}{4}, \lim_{x \to 0} b_{y} = \lim_{x \to 0} b_{y}^{2} = \frac{k_{y}^{2}}{4}$ ;所以当Δx, Δy, Δt趋于0时, 数值弥散关系式 (2.6) 、式 (2.7) 趋向于理想的数值弥散关系: ω2=c2[k $_{x}^{2}$ +k $_{y}^{2}$ ]。此外, 比较式 (2.6) 、式 (2.7) 容易看出SHO-FDTDⅠ和SHO-FDTDⅡ的格式在时间上的精度分别是一阶的和二阶的, 并且SHO-FDTDⅡ比SHO-FDTDⅠ的数值弥散误差要小得多。

2.3 数值弥散误差

下面通过试验求出格式SHO-FDTDⅠ和格式SHO-FDTDⅡ的数值弥散误差, 并与理论结果进行对比分析。令ξ=eiωΔt是方程式 (2.1) 、式 (2.4) 的根。ω=ωR+iωI, 则ξ=e-ωIΔt[cos (ωRΔt) +isin (ωRΔt) ]从而, $\tan (ω_{R} Δ t) = \frac{Ι_{m} (ξ)}{R_{e} (ξ)}$ ; 其中Im (ξ) , Re (ξ) 分别表示ξ的虚部与实部。波的数值相速vp与光速c的比的数值弥散误差: $\frac{v_{p}}{c} = \frac{ω_{R}}{k c} = \frac{1}{c k Δ t} \arctan \frac{Ι_{m} (ξ)}{R_{e} (ξ)} = \frac{Ν_{λ}}{2 π S} \arctan \frac{Ι_{m} (ξ)}{R_{e} (ξ)}$ , 其中S是CFL数, Nλ是一个波长内的节点数, k是波长值。

图4—图6给出了Vp/c在不同情况下的变化曲线。图4给出了SHO-FDTDⅠ与SHO-FDTDⅡ的Vp/c在S=0.35, Nλ=10的条件下随着角度φ的变化图, 从图4中可以看出SHO-FDTDⅡ格式的Vp/c更接近于1。图5给出了两种格式在S=2.4, φ=120°的条件下Vp/c随着Nλ的变化图, 从图5中可以看出SHO-FDTDⅡ格式的数值弥散误差小于SHO-FDTDⅠ的误差。图6是两种格式在φ=65°, Nλ=10的条件下Vp/c随着S的变化图, 从图6中可以看出随着S的增大, 数值弥散误差变得越来越大, 但是在所有的情况下, SHO-FDTDⅡ格式的数值弥散误差比SHO-FDTDⅠ格式的误差要小得多。

参考文献

[1] Yee K S.Numerical solution of initial boundary value problems invol-ving Maxwell’s equations in isotropic media.IEEE Trans Antennasand Propagation, 1966;14:302—307

[2] Zheng F, Chen Z, Zhang J.Toward the development of a three-dimen-sional unconditionally stable finite—difference time—domain method.IEEE Trans Microwave Theory Tech, 2000;48:1550—1558

[3] Namiki T.AnewFDTD algorithmbased on alternatingdirection implicitmethod.IEEE Trans Microwave Theory Tech, 1999;47:2003—2007

[4] Turkel E, Yefet A.Fourth order compact method for the Maxwell e-quations with discontinuous coefficients.Applied Numerical Mathe-matics, 2000;33:125—134

[5] Shang J S.Higher-order compact difference schemes for time depend-ent Maxwell’s equations.Journal of Computational Physics, 2003;153:312—333

[6] Xu L J, Yuan N C.高阶ADI-FDTD算法的数值色散分析.电子信息学报, 2005;27 (3) :1662—1665

[7] Xie Z Q, Zhang B, Chan C H.An explicit fourth-order staggered finitedifference time-domain method for Maxwell’s equations.Journal ofComputational and Applied Mathematics, 2002;147:75—98

[8] Manry C W, Broschat S L, Scneider J B.Higher-order FDTD meth-ods for large problems.J Applied Computational Electromagnetics So-ciety, 1995;10:17—29

时域有限元方法篇3

时域有限差分方法 (finite-difference time-domain (FDTD) method是由Yee在1966年[1]提出) , 在计算电磁学领域是一种非常常用、有效的数值计算方法[2]。在实际应用FDTD计算一种二维问题时, 时间和空间剖分步长Δt, Δx和Δy要满足Courant-Friedrichs-Lewy (CFL) 稳定性条件即: $(c Δ t)^{2} < \frac{1}{Δ x^{2}} + \frac{1}{Δ y^{2}}$ , 其中c是光的速度的限制, 因此当要模拟的问题具有细微结构时, 为了准确模拟其电磁特性, 时间步长要相应的取很小, 这时, 计算时间大大加长, 内存量大量增加, 这给实际计算带来很大困难。

为了克服这一限制条件 (CFL) , Zheng, Chen and Zhang和Namiki (没有利用分裂方程) 提出了三维和二维麦克斯韦方程的交替方向时域有限差分方法 (ADI-FDTD[3,4]) 被证明是无条件稳定, 并得到电磁学界的高度关注和广泛研究[2]。与ADI-FDTD不同, 最近, 文献[5,6]运用分裂算子 (或方程) 的方法提出了电导率为零的麦克斯韦方程的另外一种无条件稳定的有限差分方法, 称为分裂时域有限差分方法 (S-FDTD) 。S-FDTD包括两种格式, S-FDTDI和S-FDTDII, 其中前者关于时间是一阶格式, 后者利用增加摄动项的方法降低分裂误差, 是二阶格式。在分裂时域有限差分方法S-FDTDI的基础上, Chen, Li and Liang[7]利用电磁和磁场的对称性, 提出了对称的能量守恒的分裂时域有限差分 (简记为SS-FDTD ) 方法。本文将二维的S-FDTD方法和二维SS-FDTD方法推广到带有非零电导率的二维麦克斯韦方程中, 考虑这种方法对一类波导问题的计算和模拟, 并给出与ADI-FDTD的比较。

1 Maxwell方程及分裂的FDTD方法

考虑如下的Maxwell方程:

$\begin{array}{l} \frac{\partial E_{x}}{\partial t} = \frac{1}{ε} (\frac{\partial Η_{z}}{\partial y} - σ E_{x}) (1) \\ \frac{\partial E_{y}}{\partial t} = \frac{1}{ε} (- \frac{\partial Η_{z}}{\partial x} - σ E_{y}) (2) \\ \frac{\partial Η_{z}}{\partial t} = \frac{1}{μ} (\frac{\partial E_{y}}{\partial x} - \frac{\partial E_{x}}{\partial y}) (3) \end{array}$

其中E= (Ex (x, y, t) , Ey (x, y, t) ) , Hz=Hz (x, y, t) , (x, y) ∈Ω=[0, a]×[0, b], t∈[0, T], ∂Ω为Ω的边界, ε为介电常数, μ为媒质磁导率, σ为电导率。假设介质为均匀的, 且充满在矩形的边界为理想导体区域中, 则电场满足如下的边界条件

$(E, 0) \times (\vec{n}, 0) = 0 ‚ o n (0, Τ) \times \partial Ω (4)$

假设初始条件为

E (x, y, 0) =E0 (x, y) = (Ex (x, y, 0) , Ey (x, y, 0) ) ,

Hz (x, y, 0) = Hz0 (5)

为了简单起见, 只考虑ε, μ, σ为常数时的情况。对Ω采用与Yee相同的离散方法[1], 即交错网格剖分。设x方向上的步长为Δx, y方向上的步长为Δy, 令 $x_{α} = α Δ x, y_{β} = β Δ y, α = i, i + \frac{1}{2}, β = j, j + \frac{1}{2}, i = 0, 1, \dots, Ι - 1 ‚ j = 0, 1, \dots, J - 1 ‚ y_{J} = J Δ y = b$ , 对[0, T]进行等距剖分, Δt为时间步长, $t^{n} = n Δ t ‚ t^{n + \frac{1}{2}} = (n + \frac{1}{2}) Δ t, n = 0, 1, \dots, Ν - 1$ 。

为书写简单, 对函数F (x, y, t) , 定义一些差分算子, 对于m∈N+, 令

$\begin{array}{l} F_{α, β}^{m} = F (m Δ t, α Δ x, β Δ y), δ_{t} F_{α, β}^{m} = \frac{F_{α, β}^{m + \frac{1}{2}} - F_{α, β}^{m - \frac{1}{2}}}{Δ t}, \\ δ_{x} F_{α, β}^{m} = \frac{F_{α + \frac{1}{2}, β}^{m} - F_{α - \frac{1}{2}, β}^{m}}{Δ x}, δ_{y} F_{α, β}^{m} = \frac{F_{α, β + \frac{1}{2}}^{m} - F_{α, β - \frac{1}{2}}^{m}}{Δ y} \end{array}$

对于网格函数定义在Yee的交错网络点上的{Uα, β}, 在本文中令Δv=ΔxΔy, 定义

$\begin{array}{l} | U |_{E_{x}}^{2} = \sum_{i = 0}^{Ι - 1} \sum_{j = 1}^{J - 1} (U_{i + \frac{1}{2}, j})^{2} Δ v ‚ \\ | U |_{E_{y}}^{2} = \sum_{i = 1}^{Ι - 1} \sum_{j = 0}^{J - 1} (U_{i, j + \frac{1}{2}})^{2} Δ v ‚ \\ | U |_{Η}^{2} = \sum_{i = 0}^{Ι - 1} \sum_{j = 0}^{J - 1} (U_{i + \frac{1}{2}, j + \frac{1}{2}})^{2} Δ v 。 \end{array}$

根据文献[5]中的方法, 定义如下的两种格式, 分别对空间和时间进行离散得到S-FDTDII。

Stage 1:

${\begin{cases} \frac{E_{y_{i, j + \frac{1}{2}}}^{n + 1} - E_{y_{i, j + \frac{1}{2}}}^{n}}{Δ t} = - \frac{1}{2 ε} δ_{x} (Η_{z_{i, j + \frac{1}{2}}}^{*} + Η_{z_{i, j + \frac{1}{2}}}^{n}) - \\ \frac{δ}{2 ε} (E_{y_{i, j + \frac{1}{2}}}^{n + 1} + E_{x_{i, j + \frac{1}{2}}}^{n}) - \frac{Δ t}{2 μ ε} δ_{x} δ_{y} E_{x_{i, j + \frac{1}{2}}}^{n}, \\ \frac{Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{*} - Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n}}{Δ t} = - \frac{1}{2 μ} δ_{x} (E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n}) \end{cases}$

(6)

Stage 2:

${\begin{cases} \frac{E_{x_{i + \frac{1}{2}, j}}^{n + 1} - E_{x_{i + \frac{1}{2}, j}}^{n}}{Δ t} = \frac{1}{2 ε} δ_{y} (Η_{z_{i + \frac{1}{2}, j}}^{n + 1} + Η_{z_{i + \frac{1}{2}, j}}^{*}) - \\ \frac{δ}{2 ε} (E_{x_{i + \frac{1}{2}, j}}^{n + 1} + E_{x_{i + \frac{1}{2}, j}}^{n}), \\ \frac{Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} - Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{*}}{Δ t} = \frac{1}{2 μ} δ_{y} (E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n + 1} + E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n}) \end{cases} (7)$

根据文献[7]中的方法, 提出如下格式。

1) 在奇数时间层上

Stage 1:

${\begin{cases} \frac{E_{y_{i, j + \frac{1}{2}}}^{2 k + 1} - E_{y_{i, j + \frac{1}{2}}}^{2 k}}{Δ t} = - \frac{1}{2 ε} δ_{x} (Η_{z_{i, j + \frac{1}{2}}}^{*} + Η_{z_{i, j + \frac{1}{2}}}^{2 k}) - \\ \frac{δ}{2 ε} (E_{y_{i, j + \frac{1}{2}}}^{2 k + 1} + E_{x_{i, j + \frac{1}{2}}}^{2 k}), \\ \frac{Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{*} - Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{2 k}}{Δ t} = - \frac{1}{2 μ} δ_{x} (E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{2 k + 1} + E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{2 k}) \end{cases} (8)$

Stage 2:

${\begin{cases} \frac{E_{x_{i + \frac{1}{2}, j}}^{2 k + 1} - E_{x_{i + \frac{1}{2}, j}}^{2 k}}{Δ t} = \frac{1}{2 ε} δ_{y} (Η_{z_{i + \frac{1}{2}, j}}^{2 k + 1} + Η_{z_{i + \frac{1}{2}, j}}^{*}) - \\ \frac{δ}{2 ε} (E_{x_{i + \frac{1}{2}, j}}^{2 k + 1} + E_{x_{i + \frac{1}{2}, j}}^{2 k}), \\ \frac{Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{2 k + 1} - Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{*}}{Δ t} = \frac{1}{2 μ} δ_{y} (E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{2 k + 1} + E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{2 k}) \end{cases} (9)$

2) 在偶数时间层上

Stage 1:

${\begin{cases} \frac{E_{x_{i + \frac{1}{2}, j}}^{2 k + 2} - E_{x_{i + \frac{1}{2}, j}}^{2 k + 1}}{Δ t} = \frac{1}{2 ε} δ_{y} (Η_{z_{i + \frac{1}{2}, j}}^{* *} + Η_{z_{i + \frac{1}{2}, j}}^{2 k + 1}) - \\ \frac{δ}{2 ε} (E_{x_{i + \frac{1}{2}, j}}^{2 k + 2} + E_{x_{i + \frac{1}{2}, j}}^{2 k + 1}), \\ \frac{Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{* *} - Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{2 k + 1}}{Δ t} = \frac{1}{2 μ} δ_{y} (E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{2 k + 2} + E_{x_{i + \frac{1}{2}, j + \frac{1}{2}}}^{2 k + 1}) \end{cases}$

(10)

Stage 2:

${\begin{cases} \frac{E_{y_{i, j + \frac{1}{2}}}^{2 k + 2} - E_{y_{i, j + \frac{1}{2}}}^{2 k + 1}}{Δ t} = - \frac{1}{2 ε} δ_{x} (Η_{z_{i, j + \frac{1}{2}}}^{2 k + 2} + Η_{z_{i, j + \frac{1}{2}}}^{* *}) - \\ \frac{δ}{2 ε} (E_{y_{i, j + \frac{1}{2}}}^{2 k + 2} + E_{x_{i, j + \frac{1}{2}}}^{2 k + 1}), \\ \frac{Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{2 k + 2} - Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{* *}}{Δ t} = - \frac{1}{2 μ} δ_{x} (E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{2 k + 2} + E_{y_{i + \frac{1}{2}, j + \frac{1}{2}}}^{2 k + 1}) \end{cases}$

(11)

2 能量恒等式和稳定性分析

在本节中推导SS-FDTD的能量恒等式, 由此给出稳定性分析。

定理1 设 $E_{x_{i + \frac{1}{2}, j}}^{2 n}, E_{y_{i, j + \frac{1}{2}}}^{2 n} ‚ Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{2 n}$ 是格式 (8) —式 (11) 的解, 则下列恒等式成立。

$\begin{array}{l} ε | E_{y}^{2 Ν + 2} |_{E_{y}}^{2} + ε | E_{x}^{2 Ν + 2} |_{E_{x}}^{2} + μ | Η_{z}^{2 Ν + 2} |_{Η}^{2} + \\ \frac{σ}{2} Δ t \sum_{k = 0}^{Ν} {| E_{y}^{2 k + 2} + E_{y}^{2 k + 1} |_{E_{y}}^{2} + | E_{y}^{2 k + 1} + E_{y}^{2 k} |_{E_{y}}^{2} + \\ | E_{x}^{2 k + 2} + E_{x}^{2 k + 1} |_{E_{x}}^{2} + | E_{x}^{2 k + 1} + E_{x}^{2 k} |_{E_{x}}^{2}} = \\ ε | E_{y}^{0} |_{E_{y}}^{2} + ε | E_{x}^{0} |_{E_{x}}^{2} + μ | Η_{z}^{0} |_{Η}^{2} 。 \end{array}$

证明式 (8) 的第二个式子的两端同乘以 $(Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{*} + Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{2 k}) μ Δ v Δ t$ , 关于i, j求和, 得

$\begin{array}{l} μ | Η_{z}^{*} |_{Η}^{2} - μ | Η_{z}^{2 k} |_{Η}^{2} = \\ \frac{Δ t}{2} \sum_{i = 1}^{Ι - 1} \sum_{j = 0}^{J - 1} {δ_{x} Η_{z_{i, j + \frac{1}{2}}}^{*} E_{y_{i, j + \frac{1}{2}}}^{2 k + 1} + δ_{x} Η_{z_{i, j + \frac{1}{2}}}^{2 k} E_{y_{i, j + \frac{1}{2}}}^{2 k + 1} + \\ δ_{x} Η_{z_{i, j + \frac{1}{2}}}^{*} E_{y_{i, j + \frac{1}{2}}}^{2 k} + δ_{x} Η_{z_{i, j + \frac{1}{2}}}^{2 k} E_{y_{i, j + \frac{1}{2}}}^{2 k}} Δ v 。 \end{array}$

式 (8) 的第一个式子的两端同乘以 $(E_{y_{i, j + \frac{1}{2}}}^{2 k + 1} + E_{y_{i, j + \frac{1}{2}}}^{2 k}) ε Δ v Δ t$ , 关于i, j求和, 得

$\begin{array}{l} ε | E_{y}^{2 k + 1} |_{E_{y}}^{2} - ε | E_{y}^{2 k} |_{E_{y}}^{2} + \frac{σ}{2} Δ t Δ v \sum_{i = 1}^{Ι - 1} \sum_{j = 0}^{J - 1} (E_{y_{i, j + \frac{1}{2}}}^{2 k + 1} + \\ E_{y_{i, j + \frac{1}{2}}}^{2 k})^{2} = - \frac{Δ t}{2} \sum_{i = 1}^{Ι - 1} \sum_{j = 0}^{J - 1} {δ_{x} Η_{z_{i, j + \frac{1}{2}}}^{*} E_{y_{i, j + \frac{1}{2}}}^{2 k + 1} + \\ δ_{x} Η_{z_{i, j + \frac{1}{2}}}^{2 k} E_{y_{i, j + \frac{1}{2}}}^{2 k + 1} + δ_{x} Η_{z_{i, j + \frac{1}{2}}}^{*} E_{y_{i, j + \frac{1}{2}}}^{2 k} + \\ δ_{x} Η_{z_{i, j + \frac{1}{2}}}^{2 k} E_{y_{i, j + \frac{1}{2}}}^{2 k}} Δ v 。 \end{array}$

将上面两个式子相加, 得

$\begin{array}{l} ε | E_{y}^{2 k + 1} |_{E_{y}}^{2} + μ | Η_{z}^{*} |_{Η}^{2} + \frac{σ}{2} Δ t Δ v \sum_{i = 1}^{Ι - 1} \sum_{j = 0}^{J - 1} (E_{y_{i, j + \frac{1}{2}}}^{2 k + 1} + \\ E_{y_{i, j + \frac{1}{2}}}^{2 k})^{2} = ε | E_{y}^{2 k} |_{E_{y}}^{2} + μ | Η_{z}^{2 k} |_{Η}^{2} 。 \end{array}$

类似地, 由式 (9) , 也可以得到如下结果

$\begin{array}{l} ε | E_{x}^{2 k + 1} |_{E_{x}}^{2} + μ | Η_{z}^{2 k + 1} |_{Η}^{2} + \frac{σ}{2} Δ t Δ v \sum_{i = 0}^{Ι - 1} \sum_{j = 1}^{J - 1} (E_{x_{i + \frac{1}{2}, j}}^{2 k + 1} + \\ E_{x_{i + \frac{1}{2}, j}}^{2 k})^{2} = ε | E_{x}^{2 k} |_{E_{x}}^{2} + μ | Η_{z}^{*} |_{Η}^{2} 。 \end{array}$

从而由上面两个式子相加, 得到奇数时间层上有

$\begin{array}{l} ε | E_{y}^{2 k + 1} |_{E_{y}}^{2} + ε | E_{x}^{2 k + 1} |_{E_{x}}^{2} + μ | Η_{z}^{2 k + 1} |_{Η}^{2} + \\ \frac{σ Δ t Δ v}{2} {\sum_{i = 1}^{Ι - 1} \sum_{j = 0}^{J - 1} (E_{y_{i, j + \frac{1}{2}}}^{2 k + 1} + E_{y_{i, j + \frac{1}{2}}}^{2 k})^{2} + \\ \sum_{i = 0}^{Ι - 1} \sum_{j = 1}^{J - 1} (E_{x_{i + \frac{1}{2}, j}}^{2 k + 1} + E_{x_{i + \frac{1}{2}, j}}^{2 k})^{2}} = \\ ε | E_{y}^{2 k} |_{E_{y}}^{2} + ε | E_{x}^{2 k} |_{E_{x}}^{2} + μ | Η_{z}^{2 k} |_{Η}^{2} 。 \end{array}$

在偶数时间层上, 可以得到类似的结果

$\begin{array}{l} ε | E_{y}^{2 k + 2} |_{E_{y}}^{2} + ε | E_{x}^{2 k + 2} |_{E_{x}}^{2} + μ | Η_{z}^{2 k + 2} |_{Η}^{2} + \\ \frac{σ Δ t Δ v}{2} {\sum_{i = 1}^{Ι - 1} \sum_{j = 0}^{J - 1} (E_{y_{i, j + \frac{1}{2}}}^{2 k + 2} + E_{y_{i, j + \frac{1}{2}}}^{2 k + 1})^{2} + \\ \sum_{i = 0}^{Ι - 1} \sum_{j = 1}^{J - 1} (E_{x_{i + \frac{1}{2}, j}}^{2 k + 2} + E_{x_{i + \frac{1}{2}, j}}^{2 k + 1})^{2}} = \end{array}$

$ε | E_{y}^{2 k + 1} |_{E_{y}}^{2} + ε | E_{x}^{2 k + 1} |_{E_{x}}^{2} + μ | Η_{z}^{2 k + 1} |_{Η}^{2}$ 。

对于以上得到的两个式子, 关于时间求和, 有

$\begin{array}{l} ε | E_{y}^{2 Ν + 2} |_{E_{y}}^{2} + ε | E_{x}^{2 Ν + 2} |_{E_{x}}^{2} + μ | Η_{z}^{2 Ν + 2} |_{Η}^{2} + \\ \frac{σ Δ t Δ v}{2} \sum_{k = 0}^{Ν} \sum_{i = 1}^{Ι - 1} \sum_{j = 0}^{J - 1} {(E_{y_{i, j + \frac{1}{2}}}^{2 k + 2} + E_{y_{i, j + \frac{1}{2}}}^{2 k + 1})^{2} + \\ (E_{y_{i, j + \frac{1}{2}}}^{2 k + 1} + E_{y_{i, j + \frac{1}{2}}}^{2 k})^{2} \end{array}$

$+ \frac{σ Δ t Δ v}{2} \sum_{k = 0}^{Ν} \sum_{i = 0}^{Ι - 1} \sum_{j = 1}^{J - 1}$

$\begin{array}{l} (E_{x_{i + \frac{1}{2}, j}}^{2 k + 2} + \\ E_{x_{i + \frac{1}{2}, j}}^{2 k + 1})^{2} + E_{x_{i + \frac{1}{2}, j}}^{2 k + 1} + E_{x_{i + \frac{1}{2}, j}}^{2 k})^{2}} = \end{array}$

$ε | E_{y}^{0} |_{E_{y}}^{2} + ε | E_{x}^{0} |_{E_{x}}^{2} + μ | Η_{z}^{0} |_{Η}^{2}$ 。

由上面的式子可以看出, 等式右端项是一个固定值, 而等式左端的每一项都是非负的, 而且, 与CFL条件无关, 因此, 该对称格式是无条件稳定的。

3 数值验证

在本部分我们求解一个具体问题式 (1) —式 (3) , 其中, T=1, ε=1, μ=1, σ=3π, 该问题的真解为:Ex=e-πtcosπxsinπy, Ey=-e-πtsinπxcosπy, Hz=-2e-πtcosπxcosπy。设 $E_{x} (t^{n}, x_{i + \frac{1}{2}}, y_{j}) ‚ E_{y} (t^{n}, x_{i}, y_{j + \frac{1}{2}}) ‚ Η_{z} (t^{n}, x_{i + \frac{1}{2}}, y_{j + \frac{1}{2}})$ 表示真解在相应点上的值, ${E_{x_{i + \frac{1}{2}, j}}^{n}} ‚ {E_{y_{i, j + \frac{1}{2}}}^{n}} ‚ {Η_{z_{i + \frac{1}{2}, j + \frac{1}{2}}}^{n}}$ 表示差分解。令 $E r r - E = | E (t^{n}) - E^{n} |_{E} ‚ E r r - Η = | Η_{z} (t^{n}) - Η^{n} |_{Η} ‚ R a - E = \frac{| E (t^{n}) - E^{n} |_{E}}{| E (t^{n}) |} ‚ R a - Η = \frac{| Η_{z} (t^{n}) - Η^{n} |_{Η}}{| Η_{z} (t^{n}) |}$ 分别表示差分解的绝对误差和相对误差。Ra-E , Ra-H分别表示:E, Hz的收敛阶。

从表1可以看出SS-FDTD在空间上和时间上的收敛阶都为二阶。

表2给出了剖分步长Δx=Δy=Δt=0.005下, 时间长度分别为T=1, 4, 8三种情况下, 由格式SS-FDTD计算出的问题的解。由此可以看出分裂对称格式SS-FDTD是无条件稳定的 (此时, CFL数是 $\sqrt{2} > 1)$ , 而且对于长时间的计算结果也是很好的。

由此可以看出, 与 ADI-FDTD相比, SS-FDTD计算结果更好一些。

参考文献

[1] Chen W, Li X, Liang D. Symmetric energy-conserved splitting FDTD schemes for Maxwell’s equations Comm. Comput Phys, 2009;6:804—825

[2] Gao L, Zhang B, Liang D. The splitting finite-difference time-domain methods for Maxwell’s equation in two dimensions. J Comput Appl Math, 2007;205:207—230

[3] Gao L, Zhang B, Liang D.Splitting finite difference methods on stag-gered grids for the three dimensional time dependent Maxwell’s equa-tions.Commun Comput Phys, 2008;4:405—432

[4] Namiki T.A new FDTD algorithm based on alternating-direction im-plicit method.IEEE Trans Microwave Theory Tech, 1999;47 (10) :2003—2007

[5] Taflove A, Hagness S.Computational electrodynamics:the finite-differ-ence time-domain method, (second ed) .Boston, MA:ArtechHouse, 2000

[6] Yee K S.Numerical solution of initial boundary vaule problem invol-ing Maxwell’s equation in isotropic media.IEEE Trans Antcn Propa-gat, 1966;AP-14:302—307

时域有限元方法篇4

目前, 电磁场的时域计算方法越来越引人注目。时域有限差分法 (Finite Difference Time Domain, FDTD) 作为一种主要的电磁场时域计算方法, 最早是在1966年由K.S.Yee提出的。这种方法通过将Maxwell旋度方程转化为有限差分式而直接在时域求解, 通过建立时间离散的递进序列, 在相互交织的网格空间中交替计算电场和磁场。经过三十多年的发展, 这种方法已经广泛应用到各种电磁问题的分析之中。

Matlab作为一种工程仿真工具得到了广泛应用。用于时域有限差分法, 可以简化编程, 使研究者的研究重心放在FDTD法本身上, 而不必在编程上花费过多的时间。

下面将采用FDTD法, 结合Matlab的图形功能来模拟自由空间中目标的散射近场, 说明了将二者结合起来的优越性。

(二) FDTD方法的基本原理

时域有限差分法的主要思想是把Maxwell方程在空间、时间上离散化, 用差分方程代替一阶偏微分方程, 求解差分方程组, 从而得出各网格单元的场值。FDTD空间网格单元上电场和磁场各分量的分布如图1所示。

电场和磁场被交叉放置, 电场分量位于网格单元每条棱的中心, 磁场分量位于网格单元每个面的中心, 每个磁场 (电场) 分量都有4个电场 (磁场) 分量环绕。这样不仅保证了介质分界面上切向场分量的连续性条件得到自然满足, 而且还允许旋度方程在空间上进行中心差分运算, 同时也满足了法拉第电磁感应定律和安培环路积分定律, 也可以很恰当地模拟电磁波的实际传播过程。

1. Maxwell方程的差分形式

Maxwell旋度方程为:

令f (x, y, z, t) 代表E或H在直角坐标系中某一分量, 在时间和空间域中的离散形式取以下符号表示:

其中:i, j, k和n为整数。对于二维问题, 设所有物理量均为与z坐标无关, 即∂/∂z=0, 于是将 (1) 式利用二阶精度的中心有限差分式来表示函数对空间和时间的偏导数, 下面以TE波为例, 即可得到如下FDTD基本差分式:

上式中的参数Ca、bC、aD、bD的定义如下:

其中:ε, μ, σ, σm是各点处的电磁参数。在 (3) 式的推导过程中假设时间步进Δx=Δy=δ。由于计算机存储时数组下标是整数, 所以对 (3) 式进行了修改, 修改后的更新方程如下:

利用TM波和TE波之间的对偶关系, 即可编写统一适用于TE波和TM波情况的二维FDTD计算程序。

2. 数值色散及稳定性条件

为了减小数值色散, 在选取空间网格尺寸δ时, 应满足λmin≥10δ, δ=min (Δx, Δy, Δz) , 其中λmin是被研究媒质空间的最小波长值。由此可以看出:减小网格尺寸可以减小数值色散, 但是会引起计算存储量的增大, 因此需综合考虑, 权衡处理。

为了使数值计算稳定, 时间步长的选择应满足:

(三) 自由空间中目标散射近场的可视化模拟

1. 自由空间中未加散射体时的散射近场

FDTD模拟参数为:计算区域大小150×150 (网格数) , 正中的130×130为总场区, 其外是散射场区。取空间网格大小δ=λ/40, 时间步进Δt=δ/2c。平面波从左边入射, 入射角度为0°。图2是自由空间中未加散射体时的散射近场。

从图中我们得到:在自由空间内没有散射体的情况下, 整个空间就没有散射场。入射波只存在于总场区内, 散射场区域内没有电磁波。即入射波完全由总场散射场边界产生, 又由总场散射场边界吸收, 且在总场区域内也保持了平面波传播的特性, 这与事实相符。由此可见, 该方法的有效性。

2. 自由空间中金属圆柱的散射近场

FDTD模拟参数为:计算区域大小150×150 (网格数) , 正中的130×130为总场区, 其外是散射场区。圆形金属散射体的圆心位于网格正中心, 圆柱半径r=λ, λ=1×10-2m, 取空间网格大小δ=λ/40, 时间步进Δt=δ/2c, 可见圆柱半径为40个网格。平面波从左边入射到目标, 入射角度为0°。图3是位于自由空间中金属圆柱在不同时刻的散射近场可视化模拟结果。

3. 自由空间中金属方柱的散射近场

FDTD模拟参数与上述第1点相同, 只是将散射体从金属圆柱换成金属方柱。图4是位于自由空间中金属方柱在不同时刻的散射近场可视化模拟结果。

由于散射体是金属, 其散射性很强, 从图3、图4中可以看出圆柱和方柱散射体内部均没有场, 遵循了金属散射体的电磁规律。图中总场边界 (入射波在此边界引入) 处场为不连续, 这是因为在总场边界以内为总场, 以外为散射场, 目标位于总场边界内部。通过这两个算例可以看出FDTD在可视化近场时非常实用。

(四) 自由空间中目标的远场RCS

1. 自由空间中金属圆柱的远场RCS

2. 自由空间中金属方柱的远场RCS

图5和图6给出了不同形状散射体的双站RCS (用波长归一化) , 包括TE波和TM波情况。所得结果与参考文献[1]作了对比, 可以看出一致性很好。

(五) 结语

以上结合FDTD和Matlab对自由空间中目标的近场散射做了仿真分析, 所编Matlab程序简洁明了, 运行效率也较高。FDTD法在电磁场数值分析方面有很大的优越性, 而Matlab具有强大的数据处理和图形处理功能, 可以快速地编出高效高质量的程序。将二者的优势有效地结合起来, 可以将算法迅速程序化, 并获得很好的数据处理结果, 使研究者可以集中精力在FDTD方法和研究对象本身上, 而只需花费少量的时间在程序的实现上。

参考文献

[1]葛德彪, 闫玉波.电磁场时域有限差分方法[M].西安:西安电子科技大学出版社, 2005.

[2]Yee K S.Numerical solution of initial boundary value problems involving Maxwell equations in isotropic media[J].IEEE Trans, Antennas Propagat, May1996, AP-14 (3) :302-307.

[3] (美) Duane Hanselman, Bruce Littlefield著.精通Matlab7[M].朱仁峰, 译.北京:清华大学出版社, 2006.5.

时域有限元方法篇5

1 天线结构与理论计算

同轴线馈电的微带天线建模需要考虑同轴线的建模问题。对于此问题的处理可以采用阶梯近似法和一维传输线模型法[4—8]。阶梯近似法首先将同轴线分解成内、外导体及中间的介质三部分。然后对各个部件分别建立几何外形及尺寸的描述文件并进行FDTD剖分。因此接地板下方需要增加同轴线的模型空间, 同轴线的建模无疑增加了计算量。一维传输线模型法虽然不用对同轴线进行建模, 但也需要同轴线进行处理。需要将同轴线划分为接地板下的同轴线部分和馈电探针部分。同轴线部分采用一维电压V和电流I进行迭代馈电提供激励源[9]。接地板上方的馈电探针部分采用细导线模型。但是, 在同轴线与接地板的连接处需要考虑同轴线与接地板口径处的耦合问题[10]。一维传输线模型法虽然减少了计算空间, 但需要对部分空间进行特殊处理。

本文针对一维传输线模型法进行改进, 将馈电探针与贴片相连接, 探针下方与接地板之间留有一个网格作为激励源的位置, 采用点源作为激励。这样接地板与探针之间就不必采用耦合方式馈电, 于是省略了一维传输线模型的迭代。新方法建模如图1所示。

1.1 同轴探针处理

馈电探针处的处理采用细导线模型, 设细导线在 (ia, ja, k) 处。图2给出了在FDTD网格中半径为a (a<0.5Δx) 轴线与Ez分量重合的细导线。对于细导线附近的场量Hy、Ex, Hx、Ey在细导线附近有1/r的依赖关系, 以Hy、Ex为例:

对环绕细导线区域应用法拉第电磁感应定律:, 将上式带入, 得Hy的表达式为

同样的方法可以求出其他场分量。由于模型中简化了同轴线, 所以不必考虑探针处同轴线端口场量的耦合。

1.2 激励源

为了得到天线的宽频带信息激励源Vs选择高斯脉冲。激励源在细导线与接地板之间的一个网格处以电压源形式加入。电压源内阻R为50Ω, 设介质无损耗, 则激励源处的迭代式为:

1.3 电压电流取样

为了求得同轴线内电压、电流的波形, 可在探针处进行电压与电流的取样。设探针与接地板之间的电压V, 电流I;

式 (4) 中, k1, k2为接地板到探针之间的k坐标。因为探针与接地板间只有一个网格, 则上式简化为

对磁场应用安培环路定理:

值得注意的是以上的取样信息是在时间上相差半个时间步。通过得到的V、I的时域波形可以进一步求得微带天线的反射系数、输入阻抗等频域信息。

2 算例

采用以上方法编程分析了图1所示同轴探针馈电的微带天线。天线结构参数如下:贴片lx×ly=36.6 mm×26 mm, 介质厚度h=1.58 mm, 相对介电常数εr=2.17, 同轴探针半径r=0.24 mm, 馈电点位置在中线上lx, ds=6.5 mm, FDTD网格尺寸Δx=1.83 mm, Δy=1.625 mm, Δz=0.79 mm, 时间步长Δt=Δz/2c, 其中c为光速。吸收边界采用10层卷积完全匹配层 (CPML) 。图3为该微带天线的S11曲线。从图中可以看出计算结果与参考文献结果吻合很好, 说明此方法简化分析同轴馈电微带天线的可行性。图4、图5分别为微带天线的输入阻抗和方向图, 进一步说明在对于同轴线的简化问题上并不影响天线的其他参数的计算。

3 结论

简化同轴馈电天线模型同时结合同轴线的细导线分析方法分析了探针馈电的微带天线。数值结果和参考文献吻合很好, 说明这种同轴线的简化建模方式是可行的。建模过程中由于省略了对同轴线一维模型或阶梯近似方法的引入, 使得分析过程更简洁, 编程更容易。而且简化同轴线的思想并不影响天线的阻抗及辐射特性, 使得此方法更具实用性。

摘要：提出了一种简化的同轴馈电模型用于分析同轴探针馈电的微带天线。在时域有限差分方法 (FDTD) 中对于同轴馈电模型一般采用一维电压电流递推方式引入入射波形;同时结合细导线模型模拟同轴线馈电。需要对接地板同轴线接口处的电场值进行特别处理, 过程较复杂。简化方法仅采用细导线模型作为馈电端口, 在接地板与细导线探针间加电压源激励产生入射波形。省略了同轴线的建模及接口处电场值的处理过程, 使得模型更加紧凑。数值仿真验证了该算法的有效性和准确性。

关键词：时域有限差分法 (FDTD) ,同轴馈电,细导线,微带天线

参考文献

[1] 李增瑞, 王均宏, 蒋开波, 等.采用FDTD法研究微带天线雷达散射截面的减缩.电子学报, 2007;35 (6) :1065—1068

[2] 何芒, 徐晓文.探针馈电圆柱共形微带天线的FDTD法分析.微波学报, 2004;20 (3) :1—5

[3] Wong K L, Design of nonplanar microstrip antennas and transmission line.John Wiley&Sons, Inc., 1999

[4] 陈文俊, 黎滨洪, 谢涛.FDTD法分析探针加载微带天线的简便方法.电波科学学报, 2006;21 (1) :146—149

[5] 许锋, 洪伟, 冯祖伟.同轴馈电耦合微带贴片天线的时域有限差分法分析.电波科学学报, 2001;16 (1) :34—38

[6] 袁峰, 卢春兰.基于FDTD法分析计算同轴馈电贴片天线.军事通信技术, 2007;28 (2) :12—15

[7] 姜永金, 潘谊春, 傅文斌.同轴馈电耦合微带天线的MPSTD—FDTD分析.微波学报, 2011;27 (2) :19—23

[8] Supriyo D, Mitrrra R.A locally conformal finite-difference time-domain (FDTD) algorithm for modeling three-dimensional perfectly conducting objects.IEEE Trans AP, 1997;45 (7) :273—276

[9] Maloney J G, Shlager K L, Smith G S.A simple FDTD model for transient excitation of antennas by transmission lines.IEEE Trans on AP, 1994;42 (2) :289—292

时域有限元方法篇6

自1987年E.Yablonovitch提出光子带隙结构(Photonic Bandgap Structure)概念以来,其应用已涵盖光通信领域和微波领域等各方面。光子带隙结构是由周期性结构组成的,具有对特定频段的电磁波产生阻带特性,可明显提高微波集成电路、微波印刷天线等性能。研究电磁波在PBG结构中的传播特性,可发掘其潜在的应用价值:应用PBG结构于微波领域中可有效提高微波元件与系统的整体性能。尤其近年来,光子带隙结构在微波电路和天线中应用的研究受到诸多学者的关注[1,2,3,4,5]。

自20世纪80年代末人们开始关注EBG结构以来,就发展了一些用来研究EBG结构的带阻特性和色散特性的方法,如有:平面波展开法(Plane Wave Expansion,简称PWE)、转移矩阵法(Transfe Matrix Method,简称TMM)、时域有限差分法(Finite Difference Time Domain,简称FDTD)、有限元法(Finite Element Method,简称FEM)、多散射矩阵法等等[6,7,8]。时域有限差分法作为一种全面而比较精确的时域电磁场数值计算方法,几乎可以分析所有电磁问题。它的基本思想虽然简单明了,其实包含了许多数学、物理问题,而在实际应用中又有诸多的变种和改进形式,其本身已经成了一个庞大的研究方向。本文对基于时域有限差分法的二维磁带隙结构的能带研究的基本原理、主要的数值理论进行研究。

1 基本原理

1.1 格矢与倒格矢

光子带隙结构是一种周期结构,如果把每一个基本单元抽象成一个点,整个周期结构就可看成在空间里周期排列的点,这些点就被称为格点。而格点所代替的一个周期单元就是原胞,任意两个格点之间的方向矢量就是格矢,可表示为:

式(1)中,l、m、n为整数,为基矢。

与格矢相对应的倒易空间的矢量称为倒格矢,表示为:

式(2)中,l′、m′、n′为整数,为倒格基矢。以为例,满足:

所以ai·bi=2πδij(i,j=1,2,3)。

以一个倒格点为中心,作出指向邻近的各倒格点的矢量,由垂直平分这些倒格矢的平面包围的区域就是第一Brillouin区。图一即二维三角光子带隙结构的基矢、倒格基矢及第一Brillouin区的示意图。

1.2 FDTD对二维光子晶体的TM模的分析

对于无限周期结构,FDTD算法需要设定满足周期性边界条件以及较为特殊的初始场,算法的编制变得较为复杂,计算难度大。从无源、线性各向同性媒质中的Maxwell方程出发有:

必须注意到,如果计算区域空间结构是周期变化的,则电场、磁场也是伪周期的,即有:

式(8)和式(9)中,k是波矢量,E'、H'是周期函数。

通过计算可以得到:

考虑在二维的情况下,将场分解为TM波和TE波,令为二维矢波,则对TM波有:

对于TE波有:

2 数值计算

本文采用时域有限差分法计算了多种PBG结构的色散曲线,并进行了比较。以三角点阵的光子带隙结构为例,晶格常数a=6mm,介质柱的介电常数为12.29,半径r=2mm,衬底材料为有机玻璃,介电常数εb=2.56。TM模的能带图如图二所示,其中阴影部分即为全带隙,此频段的电磁波不能在光子晶体里任意方向传播。

本文对二维正方点阵结构的周期单元采用网格点总数为20×20的剖分,该正方周期的边长为12mm,相对介电常数为10.2,介质填充比为0.25。其计算结果如图三所示,与文献[8]给出的结果吻合得很好。

3 结束语

时域有限差分法作为一种全面而比较精确的时域电磁场数值计算方法,本文利用时域有限差分法是分析二维光子带隙结构的能带计算的有效方法。通过对二维光子带隙结构分析计算,和参考文献的结论吻合得很好,证明了本文对基于时域有限差分法的二维磁带隙结构的能带研究的基本原理、主要的数值理论进行研究的正确性。

参考文献

[1]C.T.Chan,W.Y.Zhang,z.L.Wang,X.Y.Lei,Da gui Zheng,W.Y.Tam,ESheng.Photonic Band Gapsfrom Metallo—Dielectric Spheres.Physica B,2792000:150-154.

[2]HanhuiFanmPierreR.Villenenve,J.D.Joanno-poulos.Large Omnidirectional Band Gaps in Metal-lodielectric Photonic Crystals.Physical Review B,1996,54(16):1245-1251.

[3]C.Y.Chang,W.C.Hsu.Photonic Bandgap Di-electric Waveguide Filter,”IEEE Microwave and Let-ters,2002,12(04).

[4]R.A.English,J.M.Pitarlce,J.B.Pendry.Order-N Effective Response of Two-dimensional Metallic Structures.SurfaceScience454-456(2000):1091-1093.

[5]S.K.Padhi,and M.E.Bialkowski.Investigations of an Aperture Coupled Microstrip Yagi Antenna Us-ing PBG Structure.IEEE AP-S2002:752-755.

[6]M.E.Zoorob,D.B.Charlton,G.J.Parker,el al.Complete and Absolute Photonic Bandgaps in Highly Symmetric Photonie QasicrystaIs Embedded in Low Refractive Index Materials.Materials Science and En-gineeringB74(2000):168-174.

[7]C.A.Kyriazidou,H.F.Contopanagos,M.Willia-am,N.G.Merrill.Artificial wersus Natural Custals:Effective Wave Impedance of Printed Photonic Bandgap Materials.IEEE Transaction on Antennas and Propagation,2000,78(01):95-105.

时域有限元方法篇7

1 导频和数据的设计

训练序列 (PN序列) 在通信系统中获得了广泛的应用。而其中的最大长度线性移位寄存器序列 (m序列) 有非常优越的相关特性, 本文利用的就是这种序列。

长度为N的m序列有如下循环自相关性:

长度为63的m序列的循环自相关结果:

本文中, 导频序列为重复两次的一个m序列, 并在前后各增加1位, 总长度为63*2+2=128个符号。在导频序列之后紧跟长度为128*19位数据, 使导频序列的开销为5%。高速数据通道的数据帧格式如图2所示。

2 基于PN序列的时域信道估计的方法

设信道冲击响应为h (n) , 对于两径信道模型

信号x (n) 经过信道, 在接收端的信号为

所以对于导频段的m序列, 在接收端的信号为

若记m序列的自相关结果为R (n) , 经过相关器后, 可以得到

当n=0时, , 因为, 故

同理

由此, 只要经过后续的相关峰检测和参数提取, 就可以得到相应的信道参数了。

3 主要模块的FPGA实现方法

3.1 相关器

相关器是信道估计中的核心模块, 后续的所有操作都是基于相关结果是否准确的基础上的, 在本文中, 采用的是基于流水线加法器的数字相关器。

数字相关器可以看做是滤波器, 即输入数据与本地码在同步时钟的驱动下, 逐级相乘并累加, 其过程如下图所示。其中, 累加的过程是一个流水线过程。

因为数据与本地码相乘改变的只是符号, 所以用符号判决来代替乘法器, 这样就大大节省了乘法器资源, 并且达到了目的。下图是用Chipscope得到的FPGA中相关器的结果。

从图中可以看出, 相关的结果是四个相关峰, 其中第一个相关峰由于有多径的缘故, 其结果受到部分数据的影响, 准确性不如第三个相关峰, 同样的, 第四个相关峰准确性不如第二个相关峰, 所以, 在峰值搜索的时候, 搜索的是第二和第三个相关峰。

3.2 相关峰搜索和I, Q路信道参数检测

相关峰搜索的目的是得到各个相关峰的位置和时延参数, 在实际的实现中, 其实就是一个状态机。

具体过程如下:

(1) 搜索态 (复位时也在此状态) , 此时, 同时对两路输入数据进行搜索, 其中第二路数据时第一路数据延时一个m序列长度产生的。注意, 在这种情况下, 相关峰是成对出现的, 所以, 当搜索到一对相关峰 (无多径时) 或两对相关峰 (有多径时) , 转入下一个状态, 否则, 一直处于搜索态。

(2) 保护态, 此状态主要是防止偶然的数据错误造成判断错误, 在此状态下, 所进行的工作和同步态一样, 但不输出同步标志和其他相关的参数, 若连续不能搜索到相关峰, 则返回搜索态。

(3) 同步态, 此状态下不再是搜索相关峰, 而是周期性检测预期位置相关峰, 若存在, 则输出同步标志, 并输出相关参数, 否则转入保护态。

当同步后, 就可以得到各个峰值的相对位置, 并且, 在同步的同时, 将各路相关的结果存入一个双口RAM, 我们只要根据各个峰值的相对位置, 就可以得到两路的信道参数了。

4 结语

本文提出了一种利用PN序列实现时域信道估计的方法, 主要介绍了如何在FPGA中实现。由于利用了PN序列的良好的自相关性, 使得时域信号估计具有更高的精度这中方法简单有效准确[3], 另外, PN序列在进行信道估计的同时, 还可以用于频偏估计和定时同步等[4,5]。这种方法简单准确, 而且在实现的过程中占用的资源比较少, 达到了所要求的目的。

参考文献

[1]杨霄鹏, 姚昆, 史浩山等.航空信道仿真研究[J].空军工程大学学报, 2006, 7 (3) :16-17.

[2]Simon Haykin.通信系统 (第4版) [M].北京:电子工业出版社:2003.

[3]姜永权, 马逸新.一种基于PN序列的OFDM时域信道估计方法[J].汕头大学学报, 2006 (11) .

[4]Tufvesson F, Edfors O, Faulkner M.Time and frequencysynchronization for OFDM using PN-sequence preambles[C].Amsterdam:Proceedings of IEEE Vehicular Technol-ogy Conference, 1999:2203-2207.

【时域有限元方法】推荐阅读：

时域方法06-11

时域反射仪02-24

时域门技术02-27