FDTD

摘要 同轴线馈电是超宽带天线常用的馈电方式。文中运用时域有限差分法

（FDTD），对同轴馈电的单极子天线的电磁特性进行时域仿真，基本思路是把同轴馈电的天线计算分为激励网格和天线网格，并通过反射场/总场的分离边界把入射场加入到天线网格中。首先对算法的原理介绍及其在仿真计算中的具体实现作介绍，同时给出了阶梯近似对同轴线馈电和单极振子天线建模的具体算法，然后实例分析了同轴馈电单极子天线的若干性能，对比粗网格模拟和细网格模拟单极振子天线所产生的不同。用FDTD法对天线的仿真分析，可用作天线设计的一种快速、

经济的辅助手段。

关键词 时域有限差分法 阶梯近似网格剖分 二阶MUR吸收边界条件 单极子天线

MODELING AND SIMULATION

OF COAXIAL FEED LINE STRUCTURES FOR

ULTRA-WIDE BAND ANTENNA

Abstract:Coaxial feed structures are widely used in ultra-wide band antennas .

This paper modeled the characteristic of the monopole antenna feeded by coaxial line by FDTD in the time-domiain,which showes that . Firstly, it introduced the theory of the arithmetic and the particularly realization in the calculation; then it described the use in the time-domain; finally it analysed several characteristics of the monopole antenna. The arithmetic used in the microstrip antenna is also a quick and economical way to design the antenna.

Key words:FDTD;2ndorder MUR absorbing-boundary; staircasing technique; monopole antenna

1.引言

超宽带天线是针对天线收发信号的相对带宽而言的：当信号的带宽与中心频率

之比小于1%时称为窄带，带宽与中心频率之比在1%与25%之间称为宽带，带宽与中心频率之比大于25%时称为超宽带。超宽带天线就是用来发射和接收超宽带信号的天线装置。而天线的辐射场要靠源来激发，如何设置符合实际的激励源，是计算天线辐射特性的关键之一。源的设置方式要与天线的馈电方式相一致，以保证天线的辐射，图形与实际一致。不同天线有不同的馈电方式，相应的激励源的设置也有所不同。同轴线是以TEM模为主模的传输线，可以传输超宽频带信号，由于其本身的一些优点，同轴线馈电是超宽带天线和高功率微波常用的馈电方式。

1966年K．S．Yee创立了计算电磁场的时域有限差分法(FDTD)[1]，它是一种计算时变电磁场有效的数值方法。它把含时间变量的Maxwell旋度方程转化为差分形式，并在时间轴上依次推进的求解空间电磁场的数值计算方法。其已在计算电磁场的各个领域得到了广泛应用，尤其是在天线分析。作为一种电磁场的数值计算方法，时域有限差分法在计算天线的特性上具有一些很突出的优点：①

作为一种瞬态方法，在脉冲波的激励下，FDTD方法的一次计算结果经Fourier变换后可获得丰富的频域信息；②适合模拟各种复杂电磁结构，用FDTD的离散空间网格可以较精确的模拟天线(阵)的实际结构；③易于得到计算空间场的暂态分布情况，这既便于定性理解其工作的物理过程，又便于得到供定量分析的有关电参量；④它所需要的计算机内存和CUP时间与网格单元成正比，并且不需要矩阵求逆，明显优于传统的矩量法。本文首先从FDTD几个关键步骤入手，从理论上加以分析说明。并用这些措施对单极子天线进行分析计算。给出电磁计算的仿真结果。因此FDTD为天线的设计提供一种快速、经济的辅助手段。下面我们将用时域有限差分法对超宽带天线同轴馈电结构实例进行研究，即对同轴馈电的单极子天线进行建模与模拟研究,天线结构如图1所示。

图 1 金属平板上圆柱天线模型

2.天线建模与仿真

本文对图1中天线进行建模和模拟研究的基本思路是把同轴馈电的天线计算

分为激励网格和天线网格，并通过反射场/总场的分离边界把入射场加入到天线网格中。这种思想将激励设置划分出来成为一个单独的网格空间(激励空间)，在此网格空间只存在与天线馈线相同的传输线，无任何其它结构，而所研究的天线结构处于另一个网格空间之内(天线结构空间)，见图2。激励空间的作用是迭代产生天线馈线入射波场，为使激励空间仅有入射波，传输线两端用吸收边界端接，在本文中吸收边界采二阶MUR吸收边界。

下面首先介绍时域有限差分基本原理，然后讨论激励网格和天线网格的建模，同时给出同轴线的激励源设置方法，最后实例分析了同轴馈电单极子天线的若干性能，对比粗网格模拟和细网格模拟单极振子天线所产生的不同点。

图2.1 天线网格与入射场网格的连接示意图

2.1时域有限差分法基本原理

2.1.1麦克斯韦方程及其差分格式

麦克斯韦方程的积分方程

lHdlDJdSst （2-1-1）

lBEdldSst （2-1-2）

BdSl0 （2-1-3） q （2-1-4）

DdSl它们相应的微分方程为

HDtBE-JtJ （2-1-5）

m （2-1-6）

B0D （2-1-7） （2-1-8）

其中E为电场强度，单位为伏特／米(V/m)； D为电通量密度，单位为库仑/米2 (C/m2)； H为磁场强度，单位为安培／米(A/m)； B为磁通量密度，单位为韦伯/米2 (Wb/m2)； J为电流密度，单位为安培/米2(A/m2)； 2Jm 为磁流密度，单位为伏特/米 (V/m)。 2各向同性线性介质中的本构关系为

D=E B=μH J=σE JmmH （2-1-9）

其中表示介质介电系数，单位为法拉／米(F／m)；μ表示磁导系数，单位为亨利／米(H／m)；σ表示电导率，单位为西门子／米(S／m)；m表示导磁率，单位为欧姆／米(Ω／m)。σ和m分别为介质的电损耗和磁损耗。真空中σ＝0，m＝0，以及

=m=8.851012F/m

μ=m4107H/m

在直角坐标系中，麦克斯韦旋度方程（2-1-5）、（2-1-6）式写为

HzyHzyExtEx （2-1-10）

HxzHxEzyyHzxHxyEyzEytEztEy （2-1-11）

Ez （2-1-12）

HxtmHx （2-1-13）

ExzEyxEzxExyHytEztmHy （2-1-14）

mEz （2-1-15）

以上六个方程是FDTD离散差分的基础。

首先，将问题空间沿三个坐标轴向分成很多网格单元，用公、却和酝分别 表示在x、y和Z坐标方向的网格空间步长，用△t表示时间步长，任意一个空间和时间的函数可表示为:

f(x,y,z,t)f(ix,jy,kz,nt)f(i,j,k) （2-1-16）

n然后用中心差分式来表示函数对空间和时间的偏导数，这种差分式实质上是 一种蛙跳法，具有二阶精度。

为了实现空间坐标的差分计算，并考虑到电磁场在空间互相正交和铰链的关系，在FDTD离散中电场和磁场各节点空间排布如图2.2，这就是著名的Yee元胞。

图 2.2 FDTD离散中的Yee元胞

Yee元胞中E、H个分量空间节点与时间步取值的整数和半整数约定如表2-1所示。

表2-1 Yee元胞中E、H各分量节点位置 空间分量取样 时间轴t取样 电磁场分量 X坐标 Y坐标 Z坐标

1 j k

Ex i

E节点 n 2i 1k

Ey j

2Ez

i j

k12

H节点

Hx Hyi

ii1212j12

kk1212

12

j

j12 n

Hz

k

下面直接给出直角坐标中三维情形下的FDTD形式：

设观察点（x,y,z）为Ex的节点，即(i1/2,j,k)，以及时刻

t2-1-10）离散为 n1/2t，于是（

1 n1Ex(i,j,k) 2 11nCA(i,j,k)Ex(i,j,k) 221 nn1111122H(i,j,k)H(i,j,k)zz 12222CB(i,j,k) 2y  11nn11112 Hy(i,j,k)Hy2(i,j,k)2222

z



（2-1-17） 式中

21tCA(i,j,k)12(i,j,k)(i2t(i1,j,k)(i122122(i,j,k),j,k)1112,j,k)t12,j,k)2(i(i12（2-1-18）

,j,k)t12,j,k)2(it121(i11212,j,k)CB(i,j,k)(i,j,k)(i,j,k)22t21（2-1-19）

,j,k)t12,j,k)(i12(i

同样（2-1-11）离散为

111n1nEy(i,j,k)CA(i,j,k)Ey(i,j,k)

222 1nn1111122 H(i,j,k)H(i,j,k)xx12222 CB(i,j,k)2x



11nn 1111Hz2(i,j,k)Hz2(i,j,k) 2222x （2-1-20）

式中CA(i,j12,k)、CB(i,j12、（2-1-19）推得。 ,k)可由（2-1-18）

同样（2-1-12）离散为

111 n1nEz(i,j,k)CA(i,j,k)Ez(i,j,k) 2221 nn1111122H(i,j,k)H(i,j,k)yy 12222CB(i,j,k) 2x  11nn1111Hx2(i,j,k)Hx2(i,j,k)

2222

y



（2-1-21） 式中CA(i,j,k)、CB(i,j,k)可由（2-1-18）、（2-1-19）推得。

2211同样（2-1-13）离散为

11 nn111111Hx2(i,j,k)CP(i,j,k)Hx2(i,j,k) 222222 11nnE(i,j1,k)E(i,j1,k)z 11z22CQ(i,j,k) 22y  11nnE(i,j,k1)E(i,j,k)yy 22 z  （2-1-22） 式中

12t12t11221221212m(i,j))1112,k1212)t12)(i,jCP(i,j12,k12),k21))m(i,j,k2(i,j,k12(i,j,k2m(i,j,km(i,j12

,k12)t12)2(i,j,k （2-1-23）

t1212111(i,jm(i,j122,k12)1CQ(i,j,k)(i,j12t,k12)m(i,j21,k,k1221)

)t12)222(i,j,k （2-1-24） 同样（2-1-14）离散为

1 n111111nHy2(i,j,k)CP(i,j,k)Hy(i,j,k) 222222 11nnE(i,j,k1)E(i,j,k)x 11x22CQ(i,j,k) 22z  11nnE(i1,j,k)E(i,j,k)zz 22 x 

（2-1-25）

12式中CA(i12,j,k12)、CB(i12,j,k12)可由（2-1-23）、（2-1-24）推得。

同样（2-1-11）离散为

1 n111111nHz2(i,j,k)CP(i,j,k)Hz(i,j,k) 222222 11nnE(i1,j,k)E(i,j,k)y y1122CQ(i,j,k) 22z  11nnE(i,j1,k)E(i,j,k)zz 22 x 

（2-1-26）

12式中CA(i12,j12,k)、CB(i12,j12,k)可由（2-1-23）、（2-1-24）推得。

2.1.2数值稳定性

由于FDTD方程只是原Maxwell旋度方程的一种近似，在计算中存在误差。

同时，由于FDTD方法是一个迭代过程，因此它的数值稳定性至关重要。 （1）数值色散问题与空间步长

ct11(x)2 （2-1-27）

21(y)1(z)2（2）数值稳定性与时间步长

x12 （2-1-28）

2.1.3介质处理

界面上的介质参数采用介质边界算术平均值条件：

(i,j,k112)11111111((i,j,k)(i,j,k)4222222111

1((i,j,k)(i,j,k)222222 （2-1-29）

(i,j,k((i1212)121111111((i,j,k)(i,j,k)4222222,k12)(i12,j12,k1)2

,j （2-1-30）

(i12,j12,k)1111111((i,j,k)(i,j,k)2222222 （2-1-31）

m(i12,j12,k)1111111((i,j,k)(i,j,k)mm2222222 （2-1-31）

对于理想金属导体界面，其切向电场为零。

2.1.3良导体中的差分格式

在Yee的差分格式中，已经给出电导率为σ的导体中的FDTD差分格式。其中，对电流项E采用的是一般的中心差分格式。对于良导体，通常t≫2，因此

t2≫1，那么方程右边第一项系数为负数，接近于-1，这将导致时间步进趋

于不稳定。

一般来说，为使解稳定，电流项E的差分格式应介于中心平均近似和前向近似之间。Luebbers就曾对HzxEytEy用Eyn12(i)Eyn1(i)代替中

心平均近似使其差分方程对良导体稳定。此外，还可以采用指数差分格式。

下面采用这两种差分格式离散麦克斯韦方程的旋度方程： （1）前向近似差分格式

以（2-1-10）

HzyHzyExtEx为例

设观察点（x,y,z）为Ex的节点，即(i1/2,j,k)，以及时刻tn1/2t，

HHn121z(i1,j,k)Hn12z(i1,j1,k)z222y2yHHn112(i1,j,k1)Hn2i1,j,k1yz()y2222zzEn11,j,k)En1x(ix(iE,j,k)x22tt 令En12x(i12,j,k)En11x(i2,j,k)，有

En121x(i12,j,k)En1x(i2,j,k) 于是（2-1-10）离散为

En11x(i,j,k)

2 CA(m)En1x(i

2,j,k)  Hn12z(i11n1211CBm)2,j2,k)Hz(i2,j2,k)

y

 n11n1 H2,k1)H21j,k1)y(i,jy(i, 2222z 

2-1-32） 2-1-33）

2-1-34）

（（（式中

(m)CA(m)t(m)t(m)2(m)211(m)t2(m)(m)t2(m)t （2-1-35）

CB(m)1(m) （2-1-36）

(m)(m)(m)tt212(m)上式中标号m(i12,j,k)。同样（2-1-11）离散为

En11y(i,jk)CA(m)En2,y(i,j12,k)n11211n211 Hx(i,j,k)Hx(i,j,kCB(m)2222)

xHn12n12z(i1j1(i12,2,k)Hz2,j12,k)x 式中标号m(i,j12,k)。同样（2-1-12）离散为

En1z(i,j,k1)CA(m)En2z(i,j,k12)n11H2(i1n211yCB(m)2,j,k12)Hy(i2,j,k2)x

Hn12x(i,j122,k1n12)Hx(i,j112,k2)y 式中标号m(i,j,k12)。

（2）指数差分格式

仍以（2-1-10）式为例，首先将（2-1-10）式变为方程

2-1-37）（2-1-38）（ ExtEx1HzHy （2-1-39）

yz解微分方程（2-1-39）式： 令

ExtExn1(i12,j,k)Ex(i1etn12,j,k)，有

Exn1(i12,j,k)Ex(i1etn12,j,k)Ex(in12,j,k)c

1HzHy式中，c

yz整理得：

Exn1(i12,j,k)etEx(in12,j,k)c1et （2-1-40）

所以，

1nn1111122Hz(i,j,k)Hz(i,j,k)2222yExn1(i12,j,k)etEx(in12,j,k)1et1)2Hny12(i12,j,k12)Hzny12(i12,j,k

即

1nn1111122H(i,j,k)H(i,j,k)zz11n1n2222Ex(i,j,k)CA(m)Ex(i,j,k)CB(m)22yny12H(i12,j,k12)Hzny12(i12,j,k1)2(m)(m) （2-1-41） 式中CA(m)e(m)(m)tt，CB(m)1e(m)

上式中标号m(iEy(i,jn112,j,k)。同样（2-1-11）离散为

n12,k)CA(m)Ey(i,j12,k)1nn1111122H(i,j,k)H(i,j,k)xx2222 CB(m)xn12

Hz(i12,j12,k)Hzxn12(i12,j1,k)2 （2-1-42） 式中标号m(i,jEzn112,k)。同样（2-1-12）离散为

n(i,j,k12ny)CA(m)Ez(i,j,k121212)12,j,k12)HCB(m)Hxn12(i12,j,k12)Hxny(i

(i,j12,k12)Hxyn12(i,j12,k1)2 （2-1-43） 式中标号m(i,j,k)。

212.1.4吸收边界条件

用FDTD分析电磁散射、辐射等开发或者半开放性质问题时,受计算机内存容量限制,不可能直接对无限的空间进行计算,因此必须在截断处设置适当的吸收边界条件,以便用有限网格空间模拟开发的无限空间.目前对吸收边界条件的比较系统和深入的研究,主要是沿着两个方向进行的： 一是通过波动方程的因子分解获得单行波方程并取近似来建立吸收边界条件；二是在边界上引入吸收材料,电磁波在无反射地进入吸收材料后被衰减掉,如PML.下面介绍几种应用较为广泛的吸收边界条件.

（1）Mur吸收边界条件 一阶Mur吸收边界条件：

fn1(P0)=f(Q0)nfctctn1n(Q0)-f(P0) （2-1-44）

二阶Mur吸收边界条件：

fn1(P0)fn1(Q0)fctnctn1(Q0)fn1(P0)2ctnf(Q0)f(P0)2(ct)nnnnnnnn(ct)2

f(P1)f(P2)f(P3)f(P4)f(Q1)fQ2)f(Q3)fQ4) （2-1-45）

图2.3 二阶吸收边界条件所涉及的10个节点 图2.4 同轴线阶梯模拟 2.2同轴线馈电单极子天线建模实例 2.2.1同轴线的矩形网格阶梯模拟

对于具有复杂几何外形的目标，首先将其按几何外形的特点分解成几个部件，然后对各个部件分别建立几何外形及尺寸的描述文件并进行FDTD剖分。虽然同轴线的几何外形并不复杂，但这种思想仍然可取。实际上同轴线就是有几个圆柱体叠合拼接而成，因此在对其建模时可以将其拆分为对多个圆柱体的建模，然后对重叠部分别处理。

同轴线的矩形网格阶梯模拟如图2.4所示，模拟的原则是:以网格的中心点为准，如果网格中心点位于位于同轴线的内外导体之间，此网格即为空气网格，否则此网格为导体网格。具体算法见附录1。

同轴线与同轴线馈电单极子建模基本思想基本一致，这里不再赘述。 2.2.2同轴线的激励源设置

同轴线激励源的最简单设置是在内外导体之间的某个位置加入激励电场，但 由于同轴线中的主模是TEM模，如果加入的电场与TEM模不完全匹配，则在激励源附近就有相应的高次模存在，经过一段距离的传输之后，高次模逐渐衰减至可以忽略不计的程度，传输线上只剩余主模，此时的场分量才能够作为同轴线的入射波加入到所要计算的同轴馈电的天线计算区域中。如果同轴线上强迫激励平面所加入的电场与TEM模基本匹配，则激励源附近产生的高次模就很小，经过很短一段距离的传输，高次模就可衰减至忽略不计的程度，同轴线上的场就是我们需要的主模场。

激励源的时间函数选用Gauss脉冲。如图2.5轴线上TEM模电场只有ρ方向分量，而FDTD网格所用的坐标系是直角坐标系，因此强迫场只能通过激励

源平面电场Ex、Ey加入，为使强迫场与同轴线主模电场一致，Ex与Ey合成场的方向应是同轴线半径方向，并且幅度与所在位置的半径成反比，具体的计算公式如下:

E(x1,y1,t)ett0/T22a(x1x0)(y1y0)22 （2-2-1）

tan1y1y0xx01 （2-2-2）

（2-2-3） （2-2-4）

Ex(x1,y1,t)E(x1,y1,t)cosEy(x1,y1,t)E(x1,y1,t)sin

图2.5同轴线场矢量分解示意图

由此知道，激励源设置算法，见附录2。 2.2.3连接面上的迭代公式

天线结构空间的计算网格划分为总场区和反射场区。如图2.1在总场区中，所有网格点用Maxwell方程的差分格式计算总场;而在反射场区中，差分格式只用于计算反射场，入射场在总场区和反射场区的连接边界加入。

假设连接面k=ks为总场区，连接面上Ey和Ex迭代时需要用到反射场区的磁场，这些磁场就需要再加上入射磁场而成为总场。迭代公式如下:

111 Exn1(i,j,ks)CA(i,j,ks)Exn(i,j,ks)222 1nn1111122 H(i,j,ks)H(i,j,ks)zz12222 CB(i,j,ks)2y



111nnn 111111Hy2(i,j,ks)Hy2(i,j,ks)Hy,i2(i,j,kout) 222222z   （2-2-5）

111n1nEy(i,j,ks)CA(i,j,ks)Ey(i,j,ks)

222 11nnn1111111222 H(i,j,ks)H(i,j,ks)H(i,j,kout)xxx,i1222222 CB(i,j,ks)2x



1n n1111122Hz(i,j,ks)Hz(i,j,ks) 2222x  

（2-2-6）

在计算连接面外二分之一空间步长(即kks12)网格点上反射场区的磁场

分量Hx和Hy时，要用到连接面上(总场区)的电场，这些电场就需要再减去入射电场而成为反射场。迭代公式如下:

11nn1111112 Hx(i,j,ks)CP(i,j,ks)Hx2(i,j,ks)222222

11nn Ez(i,j,ks)Ez(i,j,ks)1122 CQ(i,j,ks)22y 

En(i,j1,ks)En(i,j,kout)En(i,j1,ks1)yy,iy22 z 

（2-2-7） n111111nHy2(i,j,ks)CP(i,j,ks)Hy(i,j,ks) 22222 111nnnE(i,j,ks)E(i,j,kout)E(i,j,ks1)x,ix 11x222CQ(i,j,ks) 22z  n11nE(i1,j,ks)E(i,j,ks)z z22 x  （2-2-8）

12

综上，我们可以得到FDTD同轴馈电单极子天线计算流程，见附录3。 3.同轴馈电单极子天线仿真分析