用c语言模拟fdtd方法.docx

用c语言模拟fdtd方法.docx

《用c语言模拟fdtd方法.docx》由会员分享，可在线阅读，更多相关《用c语言模拟fdtd方法.docx（33页珍藏版）》请在冰豆网上搜索。

用c语言模拟fdtd方法

学号：

2009040470

哈尔滨师范大学学士学位论文

题

目

用c语言模拟三维FDTD方法研究

学

生

杨明川

指导教师

荣宪伟副教授

年

级

2009级

专

业

光电信息工程

系

别

光电信息工程系

学

院

物理与电子工程学院

哈尔滨师范大学

学士学位论文开题报告

论文题目用c语言模拟三维FDTD方法研究

学生姓名杨明川

指导教师荣宪伟副教授

年级09级

专业光电信息工程

2013年4月

本表需在指导教师和有关领导审查批准的情况下，要求学生认真填写。

说明课题的来源（自拟题目或指导教师承担的科研任务）、课题研究的目的和意义、课题在国内外研究现状和发展趋势。

若课题因故变动时，应向指导教师提出申请，提交题目变动论证报告。

课题来源：

与指导老师商讨命题

课题研究的目的和意义：

FDTD方法，即时域限差分方法（Finite-DifferenceTune-Domain）是在电磁场计算领域常用的方法。

在FDTD方法提出之后，随着计算技术，特别是电子计算机技术的发展，FDTD方法得到了长足的发展，在电磁学，电子学，光学等领域都得到了广泛的应用。

例如：

辐射天线的分析，微波器件结构的研丸，散射和雷达反射截面的计算，周期结构分析，电子封装，电磁兼容分析，微光学器件中比的传播和衍射特性硏究，核电磁脉冲的传播和散射研究。

国内外同类课题研究现状及发展趋势：

1966年，Yee首先提出麦克斯韦方程的差分离散方式，并用来处理电磁脉冲的传播和反射问题。

1969年Taylor等人用时域有限差分方法分析非均匀介质的电磁辐射，提出用吸收边界來吸收外向行波，吸收边界采用的是简单的插值方法

1971年Merewether用时域有限差分方法计算旋转体上由入射脉冲引发的感生电流，采用了辐射边界条件。

1975年，Taflove等人用时域有限差分方法计算非均匀介质在正弦波入射时的时谐场电磁散射，讨论了和谐场情况的近一一远场外推，以及数值稳定性条件

1981年，Mur提出了在计算区域截断边界处的一阶和二阶吸收边界条件及其在其时域有限差分方法中的离散形式。

这是时域有限差分方法中的一种十分有效的吸收边界条件，获得广泛应用。

目前随着电子计算机技术的飞速发展，FDTD方法得到了长足的发展，在电磁学，电子学，光学等领域都得到了广泛的应用。

课题研究的主要内容和方法，研究过程中的主要问题和解决办法：

搜集、整理课题研究材料，筛选有用的，剔除无用的，然后把有用的材料集中起來进行总结、提炼。

用调査法抽象，总结，它是有目的、有计划、有系统地搜集有关研究对象现实状况或历史状况的材料的方法。

调査方法是科学研究中常用的基本研究方法，对调查搜集到的大量资料进行分析、综合、比较、归纳，从而为人们提供规律性的知识。

研究过程中的主要问题是用什么语言编写程序可以使系统更加稳定，和如何研究需求分析，根据口C的经验和需求用WORD或相关的工具再做出-份文档系统的功能需求文档。

课题研究起止时间和进度安排:

2013年3J],进行材料收集和整理，修改开题报告初稿;

2013年4月，收集资料形成写作思路，撰写初稿；

2013年5月，初稿修改并撰写论文第二稿2013年6月，最后定稿，完成论文

课题研究所需主要设备.仪器及药品:

电脑相应的C语言计算模拟软件

外出调研主要单位，访问学者姓名:

北京航空航天大学怀进鹏院士

指导教师市査意见:

指导教师

（签字）年月

教研室（研究室）评审意见:

教研室（研究室）主任

（签字）年月

院（系）审査意见:

院（系）主任

（签字）年月

学士学位论文

题目用c语言模拟三维FDTD方法研究

学生杨明川

指导教师荣宪伟副教授

年级2009级

专业光电信息工程

系别光电信息工程系

学院物理与电子工程学院

哈尔滨师范大学

2013年4月

利用c语言模拟三维FDTD方法硏究

杨明川

摘要：

FDTD方法在当今的电磁场计算中被广泛的使用，并发挥了重要的作用.本文简要的介绍了FDTD方法的发展简史及目前的发展应用现状。

本文用一定的篇幅介绍了FDTD方法的推导过程，并逐一由一维到三维的将该方法的公式推导出来，随后对吸收边界及数值稳定性条件进行了简单讨论.最后将c语言与FDTD方法相结合，并结合实例用c语言进行计算机模拟脸证.

关键词：

麦克斯韦方程组，FDTD方法，c语言模拟

一、绪论

11引言

FDTD算法直接对麦克斯韦方程作差分处理、來解决电磁脉冲在电磁介质中传播和反射问题的算法。

基木思想是：

FDTD计算域空间节点采用Yee元胞的方法，将并离散点上的数值解来逼近连续场域内的真实解，经过时间演算，町以计算出电磁波随时间演化的规律同时电场和磁场节点空间与时间上都采用交错抽样;也即把整个计算域划分成包括散射体的总场区以及只有反射波的散射场区，这两个区域是以连接边界相连接，最外边是采用特殊的吸收边界，同时在这两个边界之间有个输出边界，用于近、远场转换；在连接边界上采用连接边界条件加入入射波，从而使得入射波限制在总场区域；在吸收边界匕釆用吸收边界条件，尽量消除反射波在吸收边界上的非物理性反射波。

FDTD算法，其空间节点采用Yee元胞的方法，电场和磁场节点空间与时间上都采用交错抽样，因而使得麦克斯韦旋度方程离散后构成显式差分方程，相比较与前面的波动方程求解，计算得到人人简化。

由于FDTD采用吸收边界条件的方法，使得计算口J■以在有限的空间范開内进行，这样就町以降低程序对计算机破件的要求。

这是是-种有效的求解电磁波的方法，可适用于任何复杂的界面，可以得到任何形状结构内部的电磁波。

1.2FDTD方法的研究现状

FDTD方法在现代各个门类的科研与技术实践中得到了广泛应用，其突出的应用表现为以下几个方面：

1.辐射天线的分析，例如柱状和锥状天线，接地导体附近的天线，喇叭天线，微带天线，手机天线，以及天线阵列。

2.微波器件的研究例如波导，微带传输，波导中的孔缝耦合加载谐振腔等

3.散射和雷达反射截面的计算，例如导体介质物体和具冇复杂结构及形状物体（导弹,飞机）的雷达截面，人体对电磁波的吸收，地下物体散射等

4.周期结构分析，例如频率选择表面，光栅传输特性，周期阵列天线，光子带隙结构,以及随机粗糙物体表面等

5.电子対装，电磁兼容分析，例如多线时装及高密度时的数字信兮传输，分析环境和结构对元器件和系统电磁参数的影响。

6.核电磁脉冲的传播和散射，在地面的反射及对电缆传输的干扰。

7.微光学器件中光的传播和衍射特性。

8.双负介质中电磁波的传播特性。

1.2.1FDTD技术的发展

11966年，YeeS先提出麦克斯韦方程的差分离散方式，并用来处理电磁脉冲的传播和反射问题。

21969年.Taylor等人用时域有限差分方法分析非均匀介质的电磁辐射，提出用吸收边界来吸收外向行波，吸收边界采用的是简单的插值方法

3.1971年Merewether用时域付限差分方法计算旋转体上由入射脉冲引发的感生电流，采用了辐射边界条件。

41975年，Taflove等人用时域有限差分方法计算非均匀介质在正弦波入射时的时谐场电磁散射，讨论了时谐场情况的近一一远场外推，以及数值稳定性条件

5.1981年,Mur提出了在计算区域截断边界处的一阶和二阶吸收边界条件及其在其时域有限差分方法中的离散形式。

这是时域有限差分方法中的一种十分有效的吸收边界条件，获得广泛应用。

61982年，Umsshankar和Taflove用时域优先差分方法计算目标雷达散射截面，提出将计算区域分为总场区和散射场区。

并提出连接边界条件，是散射计算中入射波设豐的一种简便有效的方法。

71982年Mei等人提出的共形网格技术

8.1986年，Choi和Hoefer用时域何限差分法分析波导腔体的谐振问题，计算其谐振频率。

91987年，Kasher和Yee提出的亚网格技术

101990年，Maloney等人用柱坐标系卜的时域有限差分方法分析了柱状和锥状天线位于理想导体平面上的辐射。

得到宽带天线的输入阻抗是瞬态辐射场的直观町视化報示。

111991年，丫庆等人和Luebbei*等人分别提出了三维时域有限差分方法时域近一—远场外推方法。

12.1992年，Maloney和Smith提出将阻抗边界条件应用于时域有限差分方法

ISBerenger提出（1994,1996年）提出将麦克斯韦方程扩展为场分量分裂形式，并构成完全匹配层（PML）,这是一种全新的吸收边界。

14Sacks等人（1995年）和Gedney（1996年）提出各项异性介质的PML,其支配方程是各项异性介质的麦克斯韦方程。

151999年，Prather和Shi分析轴对称衍射透镜，给出波长为微米的平面波和高斯波入射时，直径为10247微米的衍射透镜的光波传播特性。

1.3课题的研究意义和主要研究内容

13.1课题的研究意义

FDTD方法是现代电磁学理论中比较常用也比较重要的方法。

近些年FDTD发展的十分迅速，在许多方面都有重要应用，包括天线设计，微波电路设il•，电磁兼容分析，电磁散射计算,光子学应用等等。

FDTD将空间和时间进行离散化处理，将齐离散点上的数值解來逼近连续场域内的真实解，经过时间演算，町以计算出电磁波随时间演化的规律，给出精确的结构。

是一种有效的并且较为准确的求解电磁波方程的方法。

1.3.2课题的主要研究内容

本文主要研究FDTD方法的基本原理，同时对麦克斯韦方程组与FDTD方法之间的联系加以阐述并且由一维EDTD力汉、逐步深入至三维。

最后用c语言将所建立的模型加以模拟。

二、Yee元胞，麦克斯韦方程组与FDTD方程的推演

2.1FDTD方法的主要思想

FDTD算法直接对麦克斯韦方程作差分处理、來解决电磁脉冲在电磁介质中传播和反射问题的算法。

基本思想是：

TOTD计算域空间节点采用Yee元胞的方法，将各离散点上的数值解来逼近反射波的散射场区，这两个区域是以连接边界柑连接，最外边是采用特殊的吸收边界，同时在这两个近连续场域内的真实解，经过时间演算，可以计算出电磁波随时间演化的规律同时电场和磁场节点空间与时间上都采用交错抽样；也即把整个计算域划分成包括散射体的总场区以及只有边界之间有个输出边界，用于近、远场转换；在连接边界上采用连接边界条件加入入射波，从而使得入射波限制在总场区域：

在吸收边界上釆用吸收边界条件,尽量消除反射波在吸收边界上的菲物理性反射波。

2.2FDTD方程的推导

FDTD算法，其空间节点采用Yee元胞的方法，电场和磁场节点空间与时间上都采用交错抽样，冈而使得麦克斯韦旋度方程离散后构成显式差分方程，相比较与前面的波动方程求解，计算得到大人简化。

由T-FDTD采用吸收边界条件的方法，使得计算町以在有限的空间范用内进行，这样就町以降低程序对计算机硬件的要求。

这是是一种有效的求解电磁波的方法，町适用于任何复杂的界面，町以得到任何形状结构内部的电磁波。

其算法的思想具体如下：

对于麦克斯韦方程的旋度方程（方程①、②）：

VxH=—+J；

at

VxE=-—-Jm；

m

莫中Eg：

电场强度，D是电通屋密度，H是磁场强度，B是磁通量密度，J是电流密度，Jm是磁流密度，单位均为国际单位制。

在各项同性线性介质中的本构关系为方程组③:

厂D=$E

B=“H

J=aE

其中：

£是介质的介电常数.“是磁导系数.在真空中＜7=0.am=0,即两者分别

为介质的电损耗和磁损耗。

在直角坐标系中町以转化为如下一组方程即方程组④:

®2_

dz

叽

dz

dx.

5HV

dx.

驾

dz.

理.

dz

dx

%.

4

a

基+o-E

a7

=一〃讐一％H

dt

这时引入一个时空变换前数f=（x,y,z.t）代表E或H在直角坐标系中的某一分帚〉既是方程⑤：

f（凡y,乙t）=f（iAxJ4y,k&nAt）=垃⑤

对于上式怖数在时间和空间上分别取一阶偏导，并H取中心差分近似•有方程

贲（x,y,z,t）

dx

Ax

/（兀y,乙t）dz

加（x,y,乙t）

6

5f（x,y,z,t）

«

At

⑥

⑦

⑧

⑨

由上述近似方程即可以得到著名的Yee元胞如图1：

（ijk）Ey

Yee元胞

其中由图片和公式可以得到Yee尤胞中E节点和H芳点的位置关系如卜表1：

表1Yee元胞中E,H齐个分量节点位置

电磁场分最

空间分量取样

时间t取样

X坐标

y坐标

z坐标

E节点

氏

•1

1+—

2

•

J

k

n

•1

J+I

k

Ez

•

1

•

J

k+g

H节点

Hx

.1

J+7

k+-

2

1n+—

2

Hy

.1

1+—

2

■

J

k+-

2

Hz

•1

1+—

2

•1

J+I

k

可见电场和磁场在时间顺序上交替取值，时间间隔彼此相差半个时间步长■，从而可以在给定相应电磁问题的初始条件和边界条件后，利用时域冇限差分法町以逐步推进地求出以后各个时刻空间电磁场的分布。

22.1FDTD一维方程的推导

对于一维情形，TEM波只是沿着z轴传播，所有物理最均于x,y无关，故方程町以变为方程组⑩：

dz

其FDTD方程为方程丨和方程II：

1-

2耳

△匕

1+

X

（Ex）r

Az

1&（%）

杲2/M

叫=飞兀

1+二

2忖

in——（Hy）v

yk+_

At

/z1

4（EJ鴛-（EJ；：

（II）

Az

方程变为方程组⑪：

—^=-/z—-

dy&…

dEcHy—

旦一旦"生+叫

22.1FDTD二维方程的推导

对于二维情形，假设所有物理量与Z轴无关,

5HdE

——=£—+<7E.

dy&x

cH.5Ey「

-——=£—+

dx.dt

西ydExdHz口

=-//—一久Hn

axdyara2

二维情况卜电磁场的直角分呈町以划分为独立的两组：

TE波：

Ex,Ey,Ez

TM波：

Hx,Hy>Hz

TM波和TE波的Yee尤胞中E,H各个分量的位置如表2：

表2Yee元胞中E,H各个分量节点位置

电磁场分量

空间分量取样

时间t取样

龙坐标

y坐标

TE波

Ex

・1

1+—

2

j+丄

丿2

1n+—

2

Ey

・1

1+—

2

■

J

n

Ea

•

1

•1

J+I

n

TM波

Hx

•

1

•

J

n

Hy

■

1

•1

J+7

1n+—

2

Hz

・11+—

2

•

J

1n+—

2

对于TE波：

Hx=Hy=Ez=O

At

（in）

△t<7

△y

At

（IV）

Ater

1+

故其方程为方程组MIVVI

Ater

\n-H叫

、討Z川

（Hj..f.,i-（Hz）i+r1

「％列

At<7i

（%）：

J（E心（毗叶（哲

AxAy

（V）

2.2.1FDTD三维方程组的推导

对于三维情况，就需要结合图lYee元胞图进行具体推导，由图标1町以得到方程组

（2）

（3）（4）（6）（7）（8）：

At

■

Az

Ay

1

1⑴

1-

Ata1

i+

-^4ji

Ata.

1-

气屮

1+

△t<7i

5$

Ata]

（叫+

△t

叫+善7

Az

Ax

⑸

Ay

（H』祸At（aM）x.1

]+i_i

（Ejg（Eg

/叫出

△y

At

（Ey）：

和・（Ey）：

j+；、

I&（召）中斗

叫氐y严说严Jgx

7如

"中.、♦；

（E心站-（E」：

r（E心州・（EJ：

»「

AzAz

At

耳罠冲=严或二r（H'£冲一

%玮中

&）：

冲-（耳）：

冲

Ax

At

△y

“匸醞吐"比曉出+田」爲斗（EJ：

片山・（Ey）：

g斗（E・）：

沖斗・（£）：

»!

々冲斗AtI

匕述六个方程是时域有限差分方法中的一组非常基本的公式。

而在三维FDTD计算中常常采用Yee元胞为立方体•即是等间隔离散Ax=Ay=Az=由此町以得出

如下方程：

Vivnvnnxixxxi

At（7]

]_J中

（E」:

;[厂△:

;&）:

冲

Atai

At

Ata]

耳）九-他）：

氣一但石中2百住

vn

△t

Sc]

s专

Ata

1+

叫珀;

Atai

1一J叫

（E拓二广亦出+

1+W

LjZj

At

vni

（』冲号一At（aJxj+b+i

1+

（E«K-（E«kp+l-（Ey'片m+（EyL冲

IX

決总寸飞迸?

耳）：

鶴

1+—1―L

-如

帆）話*（E」：

A-（E・）二屮E・3]

At

「—％;冲

（a）n|冲At

⑴」讣片、■△WHJH」諾冲■」邑）纸冲_（Ey）ij中-（Ex>4^+（EJ4>-1.

1+2_—1+

"斗冲斗冲

XI

四、有关于FDTD算法中吸收边界，数值稳定性的浅析

31吸收边界

在TOTD算法中，要根据要求选择适当的吸收边界。

一般来说対角区域要求不高，计算机破件条件较差时，对以选用Mur边界条件，如波导、微带天线及微波电路等。

如果对于计算结果要求较高，且计算机硬件条件较好，则可以选用PML边界*件。

确定边界条件以后，要根据FDTD运算步骤，进行准备，一般为如下几步：

（1）求入射场在总场边界上的等效电磁切向分量。

（2）将计算出來的数学模型进行差分，将其写成差分格式，包括计算边界条件。

（3）将计算区域的网格离散化，将其写成差分格式，包括计算边界条件。

（4）根据选定的边界条件写出边界和角点上的场的递推。

（5）绘后进行对整个公式的迭代运算。

32数值稳定性

洋f（卡

首先我们以一维情形下的波动方程—+—f=0为例，对于平面波改c_

0）

f（x,t）=exp[-j（kx-at）]有_k?

+终f=0BPk=—

另一方面波的相速度为

=p対于无耗介质，应属于无色散。

在用有限差分代替波动方程的二阶倒数后，有二r=0・所以我们可以看到,差分近似后,k与血已经不是简单的线性对应c・

关系.上式所示的两者的非线性关系必然会导致相速度与频率冇关，因而出现色散。

现在这种色散与空间步长有关。

所以要想尽量避免这种色散，以一维波动方程为例.

工-,冬=0将方程中的二阶导数用差分近似代替，得到dxrlor

12

五、用c语言模拟FDTD方法

51用c语言模拟的初衷

在较为专业的FDTD方法研究中冇许多可供选择的软件，如OptiFDTD,Matlab甚至是FORTRAN等众多软件，这些软件具有运算速度快，貝有针对性，可视化好等优点，但木文之所以选用c语言是为为该语言具有广泛的群众基础，对于人多数从出软件研究及科研工作若来说c语言并不陌生，但是用c语言模拟该过程会存在计算速度较慢，可视化效果不好的缺陷，希望大家辩证來看。

52用c语言模拟三维FDTD方法的步骤

本文以微波传输带天线为例，来用c语言模拟FDTD方法。

5.2.1描述模拟问题

由图2我们对模拟对象会冇直观的了解，我们将坐标轴分度定为大约005mm.由此x轴定义为1245mm,y轴定义为16mm,对于z轴，I月为其尺度很小，如果单位尺度定义过人将严重影响准确性，所以此处采用非正方体单元比较有利于计算。

所以我们选择

Ax=0.389,Ay=0.4,Az=0.265单位均为毫米。

下面选择步长&。

我们选用尺寸最

Az

小的z轴來定义：

At==0.441picosecondso接着我们來转化为代码语肓。

2・Co

用Az=0.265作为基本尺寸,在x,y方向有如下因子:

AZ一AXAZ一Ay

0.265

0389

=0.6812

0.265

0.400

=0.6625

5.2.2为介质建立模型

在微带天线中有两种不同介质，自由介质以及金属介质对于自由介质：

ex[i][j][k]=gax[i][j][k]*dx[i][j][k]

我们设

gax[i][j][k]=l0/22

对于金属介质我们由图2町以得到图1,所以我们设定&Ey在k=0时均为零。

图3计算介质示意图

5.2.3激励源

此处采用一个平面波源：

ex[i][ja][k]^z[i][ja][k]+5*shape[i][k]*Fix_inc[ja]hxt1][ja-l][k]=hx[i][ja-l][k]+5*shape[i][k]*ez_inc[ja]

5.2.4进行c语言编程

将己推导的方程用c语言加以描述，整介，计算，最终输出结果结果如图4:

-2（）

-X>L_

0

IO

GHz

图4模拟结果示意图

六、结论

从以上分析及模拟不难看出，TOTD方法在电磁场的数值计算中是一种比较实用可行的方法，町以帮助我们更加形彖准确的理解场。

FDTD方法也是是一种计算复杂冃标电磁散射数值的有效方法,FDTD法在电磁场数值分析方面令很人的优越性，而c语言简单实用，实用灵活占用资源小的特点可以快速地编出高效高质量的程序。

将二者的优势仃效地结合起來,可以将算法迅速程序化，并获得较好的模拟结果，使研究者可以集中精力在FDTD方法和研究对彖本身上,而只需花费少量的时间在程序的实现上。

参考文献

[1]KANESYEE：

NumericalSolutionofInitialBoundaiyValueProblemsInvolvingMaxwelFs

EquationsinIsotropicMediaIEEETrans,Ap4.（1966）.302一307

[2]CCTaylor:

ElectromagneticPulseScattenngintime-varyinginhomogeneousmediaIEEETrans,AP-17,[1969],585-589

[3]DEMei^ewetlierTransientcuirentsonabodyofrevolutionbyanelectromagneticpulse[J]

IEEETransactionsonElectromagneticCompatibility,May1971>EMC-13

（2）：

41〜4