时域有限差分法 IIPPT课件

1稳定性对时间间隔的要求稳定性对时间间隔的要求2121nnnnffffq1q要求12tTt 即或定义数值增长因子稳定性对时间间隔的限制nnnfjtff2121012tqjq2212ttjq第1页/共80页23.2 Courant 稳定性条件稳定性条件齐次波动方程 022222222fczfyfxf平面波的解 tzkykxkjftzyxfzyxexp,0第2页/共80页3齐次波动方程差分近似齐次波动方程差分近似fxxkfxxjkxjkxfxxx2222222sinexp2exp022sin22sin22sin22222222czzkyykxxkzyx可得 2222xxxfxfxxfxf第3页/共80页4CourantCourant条件条件12t1222sin22sin22sin222222222 tzzkyykxxktczyx上式成立的充分条件是 11112222zyxtc时间间隔限制1q第4页/共80页5CourantCourant稳定性条件稳定性条件 2221111zyxtc上式即 确定空间间隔后对时间步长的限制第5页/共80页6CourantCourant条件条件几个特殊情况几个特殊情况 zyx设三维 3tc二维 2tc一维tc第6页/共80页73.3 3.3 数值色散对空间离散间隔的要求数值色散对空间离散间隔的要求考虑一维情形下波动方程 ck02222fcxf对于平面波tkxjftxfexp,0色散关系（解析式）为fxxkfxxjkxjkxfxxx2222222sinexp2exp第7页/共80页8数值色散关系数值色散关系对空间离散间隔的限制对空间离散间隔的限制离散后得到数值色散关系022sin2222 cxxk122xk要求数值色散关系和解析式相同，必须12x或对离散空间间隔的限制根据三角函数，当 时，/12sinck第8页/共80页93.4 3.4 差分近似后的差分近似后的各向异性各向异性特性特性波动方程01222222222tfczfyfxf对于平面波，差分近似后 022sin122sin22sin22sin222222222ttczzkyykxxkzyx第9页/共80页10二维情况二维情况sincoskkkkyx2222222sinsinsinsincoscoscossincv前式变为上式表明相速与平面波传播方向有关：数值各向异性 设kv第10页/共80页11传播速度空间各向异性的图示传播速度空间各向异性的图示0.00.51.0045901351802252703150.00.51.0 82043各向异性可以忽略的条件:8对空间间隔的限制图3-4-1给出了以/为参变量时相速与光速之比 v/c和平面波传播方向角之间的关系，即差分离散所带来的各向异性。第11页/共80页12 单向行波算子（Mur）吸收边界 Berenger完全匹配层（PML）吸收边界 单轴各向异性介质完全匹配层（UPML）吸收边界第四章第四章 吸收边界条件吸收边界条件第12页/共80页13开域问题：散射问题开域问题：散射问题1 截断边界使散射计算区域成为有限域截断边界散射体入射波散射波第13页/共80页14开域问题吸收边界的必要性开域问题吸收边界的必要性 从物理观点：只有在实验室墙壁上敷以吸波材料使波在此界面无反射，形成微波暗室。 从计算观点：在截断边界上FDTD的E、H分量计算公式中，至少有一个环绕分量位于截断边界之外，需要特殊处理。第14页/共80页154.1 Engquist-Majda4.1 Engquist-Majda吸收边界条件吸收边界条件波动方程012222222tfcyfxfykxktjAtyxfyxexp,平面波的解： 22222ckkkyx其中第15页/共80页16左行波和右行波左行波和右行波 ykxkktjAfyy22expykxkktjAfyy22expxf-f+入射波反射波yx0 区域界面左行波右行波 第16页/共80页17波动微分算子波动微分算子02222fkkxfy改写波动方程为定义微分算子L2222ykkxL因式分解 LLkkjxkkjxLyy2222第17页/共80页18左行波和右行波算子左行波和右行波算子00fLfL为了使截断界面处右行波，即反射波成分等于零，在截断边界处设置条件 00 xfL0022xyfkkjx对于左行波和右行波有即第18页/共80页19单向行波吸收边界条件单向行波吸收边界条件yjktcjky,1频域到时域算子过渡：01022222xfytcx0122222bxfytcx以及左侧边界右侧边界得0022xyfkkjx第19页/共80页204.2 4.2 一阶近似吸收边界条件一阶近似吸收边界条件 左行波算子21kkjkxLy利用Taylor级数展开（一阶） jkxL010 xftcx过渡导到时域第20页/共80页212211kkjkxLy02102xyfjkkjkx021110222222xfytctxc第21页/共80页22残留的反射波与入射波之比残留的反射波与入射波之比( (反射系数反射系数) ) 03060900.00.20.40.60.81.0 f+/f- 一阶近似 二阶近似接近0度（垂直入射）时反射最小第22页/共80页234.5 4.5 三维吸收边界条件三维吸收边界条件二阶近似 02111022222222xfzytctxc 432143212000101010122QfQfQfQfPfPfPfPftctcPfQftcPfQftctcQfPfnnnnnnnnnnnnnn离散后第23页/共80页24二阶吸收边界条件所涉及的二阶吸收边界条件所涉及的1010个节点（三维）个节点（三维）截断边界P0P3P2P1P4Q0Q4Q3Q2Q1第24页/共80页25一阶吸收边界条件涉及的一阶吸收边界条件涉及的2 2个节点个节点001001PfQftctcQfPfnnnn一阶近似可以应用于三维、二维、一维问题第25页/共80页26三维三维FDTD区的角顶区区的角顶区zxy三维长方体区域有6个平面UPML区，12个棱边区和8个角顶区第26页/共80页27三维截断边界上的场分量节点三维截断边界上的场分量节点Ezxyz(i0,j0,k0)(i0+1,j0+1,k0)EzExEyEz(ia,jb,kc)第27页/共80页28三维截断边界区分三种情况三维截断边界区分三种情况 6个截断边界面：面上有E切向分量和H法向分量节点；H法向分量节点不需要应用吸收边界条件 12条棱边：边上只有E的切向分量节点 8个角顶点：无电磁场分量节点Ezxyz(i0,j0,k0)(i0+1,j0+1,k0)EzExEyEz(ia,jb,kc)第28页/共80页29第五章第五章 FDTDFDTD中常用激励源中常用激励源 几种随时间变化的源 入射波的加入总场边界条件 第29页/共80页305.1 5.1 几种随时间变化的源几种随时间变化的源 时谐场源 0 sin0 0)(0ttEttEi当当 ttUEtEisin)(0或 开关函数0100)(tttU或ttttttttU00010/cos15 . 000)(（升余弦函数）（升余弦函数） 考虑到建立过程，在 (5-1-1)式所示激励源情况下达到时谐场的稳态，通常需要35个周期。当然，对于散射问题所需的周期数还与散射体大小及形状有关。例如，对于具有凹腔结构物体，Taflove等指出，达到稳定状态所需经过的周期数大约等于所模拟散射结构的Q值。为了缩短稳态建立时间，减小冲击效应，可以引入开关函数，例如采用升余弦函数 第30页/共80页31脉冲源脉冲源 高斯脉冲 微分高斯脉冲 调制高斯脉冲 升余弦脉冲 截断三余弦脉冲 截断三正弦脉冲 双指数脉冲 第31页/共80页32脉冲源脉冲源: :高斯脉冲高斯脉冲 220)(4exp)(tttEi42exp2)(220ftfjfEi时域形式 频谱 0tt 其中为脉冲峰值出现时间 决定高斯脉冲的宽度 0.10.20.30.40.50.60.700.10.20.30.40.50.60.70.80.91t / nsE(t) / v第32页/共80页33n=1024; %逆傅立叶变换采样点df=10.24/n; %频率间隔10e-11dt=1./n/df; %逆傅立叶变换的dtEit(:,2)=ifft(Ei1(:,2);第33页/共80页34不同频率时高斯脉冲的频谱值与最大值之比不同频率时高斯脉冲的频谱值与最大值之比 f0.0432 0.0948 0.4562/13 maxfEfE通常可取 为高斯脉冲的频宽 2第34页/共80页35高斯脉冲的时域波形高斯脉冲的时域波形( (带宽为带宽为0 08G) 8G) GHzf86 . 1maxnsf2 . 06 . 1max0.10.20.30.40.50.60.700.10.20.30.40.50.60.70.80.91t / nsE(t) / v第35页/共80页36高斯脉冲的频谱高斯脉冲的频谱( (带宽为带宽为0 08G)8G)0246810121416012345678910 x 10-11f / GHz|E(f)| / (v.s)FFT结果解析 结 果 GHzf86 . 1maxnsf2 . 06 . 1max第36页/共80页37脉冲源脉冲源: :微分微分高斯脉冲高斯脉冲 时域形式 频谱 42exp8)(2202ftfjfjfEi2200)(4exp)()(tttttEi第37页/共80页38微分微分高斯脉冲高斯脉冲( (频率上限为频率上限为6G) 6G) 00.511.5-0.2-0.15-0.1-0.0500.050.1t / nsE(t) / vGHzf64 . 2maxnsf4 . 04 . 2max时域波形第38页/共80页39微分微分高斯脉冲高斯脉冲( (频率上限为频率上限为6G)6G)0123456700.511.52x 10-11f / GHz|E(f)| / (v.s)FFT结果解析 结 果 频域波形特点：无零频直流分量GHzf64 . 2maxnsf4 . 04 . 2max第39页/共80页40脉冲源脉冲源: :调制调制高斯脉冲高斯脉冲 2204expcosttttEi时域形式 频谱 00220002202exp4exp4 2exp4exp4)(tffjfftffjfffEi第40页/共80页41调制调制高斯脉冲（中心频率为高斯脉冲（中心频率为6.5G) 6.5G) 0123456-1-0.8-0.6-0.4-0.200.20.40.60.81t / nsE(t) / v时域波形第41页/共80页42调制调制高斯脉冲（中心频率为高斯脉冲（中心频率为6.5G)6.5G)024681012140246810121416x 10-11f / GHz|E(f)| / (v.s)解析 结 果 FFT结果频域波形特点：有中心频率和带宽第42页/共80页435.6 5.6 入射波的加入入射波的加入总场边界条件总场边界条件 散射体截断边界散射场区总场区总场边界入射波将计算区域划分为总场区和散射场区 sincsincHHHEEE第43页/共80页44应用等效原理设置入射波应用等效原理设置入射波 AincHnJAincmEnJ 为了使入射波限制在总场区内的空间有限区域，根据等效原理，在区域界面A上设置等效面电磁流，并设A面外的场为零 在总场-散射场区的分界面上设置入射波电磁场的切向分量便可将入射波只引入到总场区。 incEincHn 零场总场区散射场区第44页/共80页45三维情况总场三维情况总场- -散射场边界散射场边界 ZXY 总场边界总场外边界a b00,kjia00,kjibckji,00000kji,21,21,2100kjia21,21,210kjiba21,21,210cakji21,21,2100ckji总场边界：6个面，12条棱边，角顶点处无节点第45页/共80页46MaxwellMaxwell旋度方程对于总场或散射场均适用旋度方程对于总场或散射场均适用xxyzEtEzHyHxmxyzHtHzEyE电场 x 分量磁场 x 分量第46页/共80页47MaxwellMaxwell旋度方程旋度方程的的FDTDFDTD对于总场或散射场均适用对于总场或散射场均适用zkjiHkjiHykjiHkjiHmCBkjiEmCAkjiEnynynznznxnx21,2121,21,21,21,21,21)(,21)(,21212121211要求：FDTD公式中涉及的所有场分量（电场、磁场）或者属于总场，或者属于散射场第47页/共80页48总场总场- -散射场边界附近的场分量散射场边界附近的场分量 EzEzHxHy(i0,j)(i0,j0)(ia,j0)(ia,jb)(i0,jb)(i0-1/2,j0-1/2)(ia+1/2,j0-1/2)(i0-1/2,jb+1/2)(ia+1/2,jb+1/2)总场区总场边界散射场区总场外边界HxHxHxHyHyHyEzEzEz/2/2xyz(i,jb)(ia,j)(i,j0) 在为总场边界上属于总场区。 距离总场边界半个网格处为总场外边界，其上磁场分量属于散射场区。 第48页/共80页49总场边界上的总场边界上的E E分量（属于总场）被分量（属于总场）被4 4个个H H分量环绕，其中位于散分量环绕，其中位于散射场区的分量属于散射场射场区的分量属于散射场 (i,j,k)元胞中心元胞中心元胞中心元胞中心(i-1/2,j+1/2,k+1/2)EzHyHxHxHyxzyxy第49页/共80页50总场外边界上的总场外边界上的H H分量（属于散射场）被分量（属于散射场）被4 4个个E E分量环绕，其中位分量环绕，其中位于总场区的分量属于总场于总场区的分量属于总场(i+1,j,k)元胞中心(i+1/2,j+1/2,k-1/2)HzEyExExEyxzy元胞中心(i+1/2,j+1/2,k+1/2)(i,j+1,k)xy第50页/共80页51总场边界处总场边界处FDTDFDTD公式的修改公式的修改 在总场散射场边界处应用FDTD公式时，同一公式中可能既有总场节点，也有散射场节点。因此，原来FDTD公式需要修改。 认定FDTD公式中全部节点为总场，则要在散射场节点加上入射场；认定全部节点为散射场，则要在总场节点减去入射场。 修改后结果正好与等效原理一致，在总场区加入了入射波。第51页/共80页52二维二维TMTM情况总场情况总场- -散射场边界散射场边界 EzEzHxHy(i0,j)(i0,j0)(ia,j0)(ia,jb)(i0,jb)(i0-1/2,j0-1/2)(ia+1/2,j0-1/2)(i0-1/2,jb+1/2)(ia+1/2,jb+1/2)总场区总场边界散射场区总场外边界HxHxHxHyHyHyEzEzEz/2/2xyz(i,jb)(ia,j)(i,j0) 在为总场边界上属于总场区。距离总场边界半个网格处为总场外边界，其上磁场分量属于散射场区。 zEyxHH ,第52页/共80页53总场边界附近元胞（二维总场边界附近元胞（二维TMTM） (i-1/2,j0)(i,j0+1/2)HyHxyxHx(i,j0)Ez(i,j0-1/2)(i+1/2,j0)Hy总场区总场边界散射场区下边界第53页/共80页54总场边界附近总场边界附近FDTDFDTD公式需改写公式需改写 yjiHtjiEjiEnincxFDTDnznz)21,(),(),(021,001入射场总场yjiHjiHxjiHjiHmCBjiEmCAjiEnxnxnynynznz21,21,21,21)(),()(),(212121211第54页/共80页55总场边界附近总场边界附近FDTDFDTD公式需改写公式需改写 yjiEtjiHjiHninczFDTDnxnx),(21,)21,(0,021021入射场散射场112200,00011( ,)( ,)22( ,)( ,)( ,1)nnxxnnnz inczzHi jHi jEi jEi jEi jttyy第55页/共80页56总场边界附近总场边界附近FDTDFDTD公式需改写公式需改写 FDTDnynyjiHjiH021021,21),21(总场总场xjiEjiEmCQjiHmCPjiHnznznyny),(), 1()(,21)(,212121第56页/共80页57二维二维TMTM波总场区四个角点波总场区四个角点 zExjiHyjiHtjiEjiEnincynincxFDTDnznz),21()21,(),(),(0021,0021,00001以左下角点为例 基本思想：认定公式全部节点属于总场（散射场），则应当在散射场（总场）节点处加上（减去）入射场值。入射场总场yjiHjiHxjiHjiHmCBjiEmCAjiEnxnxnynynznz21,21,21,21)(),()(),(212121211第57页/共80页58二维二维TMTM波总场区四个角点波总场区四个角点 zEEzEzHxHy(i0,j)(i0,j0)(ia,j0)(ia,jb)(i0,jb)(i0-1/2,j0-1/2)(ia+1/2,j0-1/2)(i0-1/2,jb+1/2)(ia+1/2,jb+1/2)总场区总场边界散射场区总场外边界HxHxHxHyHyHyEzEzEz/2/2xyz(i,jb)(ia,j)(i,j0)图5-6-3 二维总场-散射场边界第58页/共80页59二维二维TMTM波入射场设置算例波入射场设置算例（幅值分布）（幅值分布） nm35140775. 8nm总场边界：80，80；80，80 第59页/共80页60二维二维TMTM波入射场设置算例波入射场设置算例（相位分布）（相位分布）第60页/共80页61一维情况总场一维情况总场- -散射场区散射场区散射场 区总场 区 d00,0r0划分总场区和散射场区后，介质板置于总场区；反射波和透射波便于提取第61页/共80页62一维总场一维总场- -散射场边界散射场边界xEyHyHz散射场区总场区总场边界注意：总场边界上的E节点属于总场； 总场外边界处的H节点位于散射场区，属于散射场。第62页/共80页63一维总场边界上入射波的加入（电场公式）一维总场边界上入射波的加入（电场公式） 21021,0101kHztkEkEniyFDTDnxnx总场边界上的电场节点属总场，而总场外边界处磁场节点属散射场，应当加上入射场。FDTD公式修改为 212121211kHkHztkEkEnynynxnx通常FDTD电场公式应用于总场边界时总场 总场 散射场入射场总场第63页/共80页64一维总场边界上入射波的加入（磁场公式）一维总场边界上入射波的加入（磁场公式） 121212121kEkEztkHkHnxnxnyny 0,0210212121kEztkHkHnixFDTDnyny总场外边界处磁场节点属散射场，而总场边界上的电场节点属总场，应当减去入射场。FDTD公式修改为 通常FDTD磁场公式应用于总外场边界时散射场 总场 散射场入射场散射场第64页/共80页65一维总场边界入射波加入算例：一维总场边界入射波加入算例：400dt400dt时的波形时的波形 m1cdzt2FDTD区为1,1000。 吸收边界为一阶MUR。总场边界在500。在总场边界加入右行波。02004006008001000-1.0-0.50.00.51.0Ext/dt B第65页/共80页66一维总场边界入射波加入算例：一维总场边界入射波加入算例：1200dt1200dt时的波形时的波形 02004006008001000-1.0-0.50.00.51.0Ezt/dt B第66页/共80页67通过总场边界加入入射波通过总场边界加入入射波 确定总场边界（三维，二维，一维）； 在总场边界处，修改通常FDTD公式，引进入射波（切向分量）； 总场边界上的等效面电磁流将只在总场区内产生入射波； 检验：平面波入射（总场区内无散射体）第67页/共80页68第六章第六章 近近远场外推远场外推电磁场的远场和近场划分(转) 电磁辐射源产生的交变电磁场可分为性质不同的两个部分，其中一部分电磁场能量在辐射源周围空间及辐射源之间周期性地来回流动，不向外发射，称为感应场；另一部分电磁场能量脱离辐射体，以电磁波的形式向外发射，称为辐射场。 一般情况下，电磁辐射场根据感应场和辐射场的不同而区分为近区场（感应场）和远区场（辐射场）。由于远场和近场的划分相对复杂，要具体根据不同的工作环境和测量目的进行划分，一般而言，以场源为中心，在三个波长范围内的区域，通常称为近区场，也可称为感应场；在以场源为中心，半径为三个波长之外的空间范围称为远区场，也可称为辐射场。 第68页/共80页69第六章第六章 近近远场外推远场外推 近区场通常具有如下特点： 近区场内，电场强度与磁场强度的大小没有确定的比例关系。一般情况下，对于电压高电流小的场源(如发射天线、馈线等)，电场要比磁场强得多，对于电压低电流大的场源(如某些感应加热设备的模具)，磁场要比电场大得多。 近区场的电磁场强度比远区场大得多。从这个角度上说，电磁防护的重点应该在近区场。 近区场的电磁场强度随距离的变化比较快，在此空间内的不均匀度较大。远区场的主要特点如下： 在远区场中，所有的电磁能量基本上均以电磁波形式辐射传播，这种场辐射强度的衰减要比感应场慢得多。 在远区场，电场强度与磁场强度有如下关系：在国际单位制中，E=377H，电场与磁场的运行方向互相垂直，并都垂直于电磁波的传播方向。 远区场为弱场，其电磁场强度均较小第69页/共80页706.1 6.1 等效原理等效原理目标E, HE, H自由空间界面A原问题：例如辐射或散射第70页/共80页71等效原理：自由空间辐射等效原理：自由空间辐射E, H零场自由空间界面A等效问题：n 区域内为零场；n 界面A上有面电流和面磁流；n 面电磁流在自由空间产生辐射场。自由空间HeJn EeJnmne 第71页/共80页72外推边界（输出边界）外推边界（输出边界）散射场区散射场区总场区目 标连接边 界输出边 界吸收边 界外推输出边界在散射场区内第72页/共80页73等效面电磁流辐射的计算公式等效面电磁流辐射的计算公式 HnJ EnJmVdrrGrJrFVdrrGrJrAm),()()(),()()(面电磁流辐射场的势函数其中G 为Green函数（远场及近场）第73页/共80页74由势函数到电磁场的公式由势函数到电磁场的公式AjAjFAjFE11FjFjAFjAH11VdrrGrJrFVdrrGrJrAm),()()(),()()(其中第74页/共80页75Visualization of plane-wave penetration and scattering of a missile radome containing a horn antenna (from Taflove, 2000)The impinging plane wave propagates from right to leftand is obliquely incident at 15 from boresight 第75页/共80页76Visualization of the FDTD-computed specific absorption rate (SAR) distribution of a head model (from Taflove, 2000)An one-quarter wavelength whip antenna is in anear-level plane for a 1,900-MHz cellular telephone held vertically against a tilted-head model.第76页/共80页77Surface currents on “Rund” aircraft model for vertical polarization (from Andersson 2001) Rund 是瑞典空军研究所一种飞机模型，机身长1米，翼展1米，高0.5米第77页/共80页78Surface currents of the Saab 2000 aircraft 125 ns (1500 time steps) after a lightning stroke the nose (from Andersson 2001)Also the magnitude of magnetic field is shown in a cutting plane across the wings perpendicular to the fuselage. 第78页/共80页79Surface currents on the interior of the Saab 2000 (from Andersson 2001)The view is from the center of aircraft towards to the cockpit. 第79页/共80页80感谢您的观看。第80页/共80页