王新柯07B911020 物理系时域有限差分算法对于生物电磁学的应用博士研究生：王新柯学号：07B911020导师：张岩所在单位：物理系- I -王新柯07B911020 物理系目录第1章时域有限差分算法简介 (3)第2章时域有限差分方法的基本理论 (7)2.1 麦克斯韦方程和Yee元胞 (7)2.2 直角坐标中的三维FDTD (10)第3章FDTD对于生物电磁学的应用 (19)3.1 生物组织的电磁特性和人体的电磁模型 (19)3.2 平面波照射下人体内的电磁效应 (23)3.3 动力电的人体效应 (30)3.4 高压EMP的生物效应 (34)第4章结论 (38)参考文献 (39)- II -王新柯07B911020 物理系第1章时域有限差分算法简介自1873年麦克斯韦建立电磁场基本方程以来，电磁波理论和应用的发展已经有一百多年的历史。

目前，电磁波的研究已深入到各个领域，应用十分广泛，例如无线电波传播、光纤通信和移动通信、雷达技术、微波、天线、电磁成像、地下电磁探测、电磁兼容，等等。

电磁波在实际环境中的传播过程十分复杂，例如各种复杂目标的散射，复杂结构天线的影响，等等。

具体实际地研究电磁波的特性有着十分重要的意义。

实验和理论分析技术是相辅相成的重要手段。

分析计算途径需要结合实际环境电磁参数求解麦克斯韦方程边值问题，通常只有一些经典问题有解析解。

应当说，解析解具有重要的指导意义。

然而，由于实际环境的复杂性，往往需要通过数值解得到具体环境下的电磁波特性。

随着计算机技术的发展，已经提出求解麦克斯韦方程的很多有意义的数值解方法，例如矩量法、有限元法、边界元法已经时域有限差分方法等等。

并且，随着电磁波应用的广泛和计算机技术的发展，各种方法的研究也更加深入。

1966年K. S. Yee首次提出了一种电磁场数值计算的新方法——时域有限差分（Finite Difference Time Domain, FDTD）方法[1]。

对电磁场E、H分量在空间和时间上采用交替抽样的离散方式，每一个E（或H）场分量周围有四个H（或E）分量环绕，应用这种离散方式将含有时间变量的麦克斯韦旋度方程转化为一组差分方程，并在时间轴上逐步推进地求解空间电磁场。

Yee提出的这种抽样方式后来被称为Yee元胞。

FDTD方法是求解麦克斯韦微分方程的直接时域方法。

在计算中将空间某一样本点的电场或磁场与- 3 -王新柯07B911020 物理系周围格点的磁场或电场直接相关联，且介质参数已赋值给空间每一个元胞，因此这一方法可以处理复杂形状目标和非均匀介质物体的电磁散射、辐射等问题。

同时，FDTD的随时间推进可以方便地给出电磁场的时间演化过程，在计算机上以伪彩色方式显示，这种电磁场可视化结果清楚地显示了物理过程，便于分析和设计。

FDTD方法是求解麦克斯韦方式的直接时域方法，经过三十多年的发展已成为一种成熟的数值方法，目前已被广泛研究和应用。

以下只是简单回顾FDTD的发展：1966年，Yee首次提出麦克斯韦方程的差分离散方式，并用来处理电磁脉冲的传播和反射问题[1]。

1969年，Taylor等用FDTD分析非均匀介质的电磁散射，提出用吸收边界来吸收外向行波，吸收边界采用的是最简单插值方法[2]。

1971年，Merewether用FDTD计算旋转体上有入射脉冲所引起的感生电流，采用了辐射边界条件[3]。

1975年，Taflove等用FDTD计算非均匀介质在正弦波入射的时谐场（稳态）电磁散射，讨论了时谐场情况的近—远场外推，以及数值稳定性条件[4]。

Holland在1977年和Kunz在1978年，分别独立采用FDTD计算F117飞机这种复杂目标的电磁波散射[5，6]。

1981年，Mur提出在计算区域截断边界处的一阶和二阶吸收边界条件及其在FDTD的离散形式。

这是FDTD的一种十分有效的吸收边界条件，获得广泛应用[7]。

- 4 -王新柯07B911020 物理系1982年，Umashankar和Taflove用FDTD计算目标雷达散射截面，提出将FDTD区划分成总场区和散射场区，并提出连接边界条件，是散射计算中入射波设置的一种简便有效方法[8]。

Umashankar和Taflove等在1987年和1988年，利用FDTD分析了自由空间及腔体中导线上的感应电流，讨论了FDTD中细导线的处理方法[9，10]。

Choi和Hoefer于1986年，用FDTD分析了波导腔体的谐振问题，计算其谐振频率[11]。

Kasher和Yee在1987年，提出亚网格技术[12]，Mei等人在1984年提出共形网络技术[13]。

Zhang和Mei在1988年，Liang等在1989年，Gwarek在1988年，Sheen和Kong等在1990年，用FDTD分析计算了波导、同轴线、微带天线及微带不连续性问题，得到相应的阻抗、传播常数及S参数[14，15，16，17]。

Maloney等在1990年，用圆柱坐标下FDTD分析了柱状和锥状天线位于理想导体平面上的辐射，得到宽带天线的输入阻抗及瞬态辐射场的直观可视化显示[18]。

Sullivan在1990年，用FDTD计算了60-70 MHz电磁波照射下透入到人体内部的电磁场，研究了生物电磁学问题[19]。

Britt在1989年，首次给出了时域远场结果[20]；Yee等和Luebbers等分别在1991年提出了三维FDTD时域近—远场外推方法[21，22]；随后，Luebbers等在1992年提出了二维FDTD时域近—远场外推方法[23]。

Larson, Perlik和Taflove等分别在1989年提出了研究适用于FDTD的专用计算机，以便计算电磁波与大尺寸物体的相互作用[24，25]。

- 5 -王新柯07B911020 物理系Luebbers和Hunsberger等在1990年，研究了色散介质在FDTD中的处理方法[26]。

Maloney和Smith在1992年，提出将阻抗边界条件应用于FDTD[27]。

Sui等于1992年，提出用二维FDTD计算有集中参数元件的数字和微波电路模型，包括电阻、电容、电感、二极管、晶体管等元件[28]。

Berenger在1994年，提出将麦克斯韦方程扩展为场分量分裂形式，并构成完全匹配层（PML），这是一种全新的吸收边界[29]；Sacks等人在1995年和Gedney在1996年分别提出各向异性介质的PML，其支配方程是各向异性介质麦克斯韦方程[30，31]；在FDTD计算中这两种PML作为吸收边界已得到广泛应用。

Prather和Shi在1999年，分析轴对称衍射透镜，给出波长为平面波和高斯波入射时，直径衍射透镜的光波传播特性[32]。

- 6 -王新柯 07B911020 物理系- 7 -第2章 时域有限差分方法的基本理论2.1 麦克斯韦方程和Yee 元胞麦克斯韦方程组是支配宏观电磁现象的一组基本方程。

这组方程即可以写成微分形式，又可以写成积分形式。

FDTD 方法是由微分形式的麦克斯韦方程出发进行差分离散。

众所周知，麦克斯韦为旋度方程可写为，D H J t ∂∇⨯=+∂ （1） B E J t∂∇⨯=--∂ （2） 其中，E 为电场强度，D 为电通量密度，H 为磁场强度，B 为磁通量密度，J 为电流密度，m J 为磁流密度。

各向同性线性介质中的本构关系为，,,,.m m D E B H J E J H εμσσ==== （3） 其中ε表示介质介电系数，μ表示磁导系数，σ表示电导率，m σ表示磁导率。

σ和m σ分别为介质的电损耗和磁损耗。

由于真空中，有0,0m σσ==，以及1208.8510/F m εε-==⨯和70410/H m μμπ-==⨯。

在真空中，式（1）和（2）可以写为差分方程，王新柯 07B911020 物理系 - 8 -y x z x y x z y y x z z H E H E y z tE H H E z x tH H E E x y tεσεσεσ∂∂∂-=+∂∂∂∂∂∂-=+∂∂∂∂∂∂-=+∂∂∂ （4） 以及y x z m x y x z m y y x z m z E H E H y z tH E E H z x tE E H H x y tμσμσμσ∂∂∂-=--∂∂∂∂∂∂-=--∂∂∂∂∂∂-=--∂∂∂ （5） 下面考虑式（4）和（5）的FDTD 的差分离散。

令(),,,f x y z t 表示E 或H 在直角坐标系中某一分量，在时间和空间域中的离散取以下符号表示：()()(),,,,,,,,n f x y z t f i x j y k z n t f i j k =∆∆∆∆= （6）对(),,,f x y z t 关于时间和空间的一阶偏导数取中心差分近似，即()()()()()()1/21/211,,,,,,,22|11,,,,,,,22|11,,,,,,,22|,,,,,,,|n n x i x n n y j y n n z k z n n t n t f i j k f i j k f x y z t x xf i j k f i j k f x y z t y y f i j k f i j k f x y z t z tf x y z t f i j k f i j k t =∆=∆=∆+-=∆⎛⎫⎛⎫+-- ⎪ ⎪∂⎝⎭⎝⎭≈∂∆⎛⎫⎛⎫+-- ⎪ ⎪∂⎝⎭⎝⎭≈∂∆⎛⎫⎛⎫+-- ⎪ ⎪∂⎝⎭⎝⎭≈∂∆∂-≈∂t∆ （7）王新柯 07B911020 物理系- 9 -在FDTD 离散中电场和磁场各节点的空间排布如图1所示，这就是著名的Yee 元胞。

由图可见每一个磁场分量由四个电场分量环绕；同样，每一个电场分量由四个磁场分量环绕。

这种电磁场分量的空间取样方式不仅符合法拉第电磁感应定律和安培环路定理得自然结构，而且这种电磁场各分量的空间相对位置也适合于麦克斯韦方程的差分计算，能够恰当地描述电磁场的传输特性。

此外，电场和磁场在时间顺序上交替抽样，抽样时间间隔彼此相差半个时间步长，使麦克斯韦旋度方程离散以后构成显示差分方程，从而可以在时间上迭代求解，而不需要进行矩阵求逆运算。

因此，由给定相应电磁问题的初始值，FDTD 方法就可以逐步推进地求得以后各个时刻空间电磁场的分布。

图1 FDTD 离散中的Yee 元胞。

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

王新柯 07B911020 物理系- 10 - 表1 Yee 元胞中,E H 各分量节点位置2.2 直角坐标中的三维FDTD现在考虑三维情况下的FDTD 形式。

先看（4）式中第一式，设观察点(),,x y z 为x E 的节点，即()1/2,,i j k +，以及时刻()1/2t n t =+∆，于是，（4）式中第一式的离散形式为，王新柯 07B911020 物理系 - 11 -111/21/21/211,,,,122,,211,,,,122 ,,221111,,,,222211,,22 n n x x n n x x n n z z n n y y E i j k E i j k i j k t E i j k E i j k i j k H i j k H i j k yH i j k H εσ++-++-⎛⎫⎛⎫+-+ ⎪ ⎪⎛⎫⎝⎭⎝⎭+ ⎪∆⎝⎭⎛⎫⎛⎫+++ ⎪ ⎪⎛⎫⎝⎭⎝⎭++ ⎪⎝⎭⎛⎫⎛⎫++-+- ⎪ ⎪⎝⎭⎝⎭=∆⎛⎫++- ⎪⎝⎭-1/211,,22i j k z⎛⎫+- ⎪⎝⎭∆ （8） 其中用了平均值近似，即11/211,,,,122,,22n n x x n x E i j k E i j k E i j k ++⎛⎫⎛⎫+++ ⎪ ⎪⎛⎫⎝⎭⎝⎭+= ⎪⎝⎭（9） 这样处理是为了在离散式中只出现表1所示的各个场分量节点。