《电磁场问题数值算法》经典习题解答.pdf

电 磁 场 问 题 数 值 算 法 学 习 经 典 习 题吕波矩 量 法问 题 ： 求 解 金 属 圆 柱 在 水 平 入 射 均 匀 平 面 波 照 射 下 ， 面 电 流 分 布 及 由 此 产 生 的 远区 散 射 场 分 布 。分 析 ： 可 将 面 电 流 密 度 用 脉 冲 函 数 展 开 ， 即 ：1 ( )N isz i iiJ r r （ 1-1）它 表 示 把 导 体 横 截 面 的 周 界 分 成 N个 小 段 ， 且 设 每 一 个 小 段 中 的 面 电 流 密 度 i 作均 匀 分 布 。 这 样 算 子 方 程 为 ： 201 1( ) ( ) ( | |)4 iN Ni isz i i i li i kL J L r r H k r r dl 若 以 狄 拉 克 函 数 作 为 权 函 数 切 令 它 定 义 在 各 分 段 的 中 点 ( , )j jx y 上 ， 即 令( , )j j jW x x y y （ 1-2）则 积 分 方 程 可 以 转 化 为 代 数 方 程 组 l g （ 1-3）式 中 ( , ), ( , ) i ij j j z z j jg x x y y E E x y （ 1-4）(2) 2 20( , ), ( ) ( ( ) ( )4 iiij ji j j i i j i j ilkl l x x y y L r r H k x x y y dl （ 1-5）当 i j 时 ， 采 用 近 似 式 (2) 2 20 ( ( ) ( ) )4ij ji i j i j ikl l H k x x y y l （ 1-6）当 i j 时 ， 21 (ln 1)4 4 iii i Dk lkl l j （ 1-7） g中 个 元 素 可 按 下 面 的 分 析 定 之 。 设 自 由 空 间 中 有 一 均 匀 平 面 波 沿 r 方 向 传 播 ，r 和 x轴 之 间 的 夹 角 为 i， 则 在 空 间 任 意 点 （ x,y） 处 的 电 场 强 度 可 表 示 为 ：( cos sin )( ) i ijk x yizE e 导 体 表 面 上 匹 配 点 处 的 电 场 强 度 应 为( cos sin )( )( , ) j i j ijk x yiz j jE x y e （ j=1,2, ,N） （ 1-8）Matlab源 程 序 如 下 ： %矩 量 法 计 算 圆 柱 散 射%copyright by lvbo% 2003,6,24clear;yita=120*pi%自 由 空 间 波 阻 抗freq=5*108%频 率 为 500MHzk=(2*pi*freq)/(3*108)%该 频 率 的 传 播 常 数N=100%将 圆 柱 导 体 横 截 面 的 周 界 l分 成 N个 小 段r=0.1%圆 柱 导 体 的 半 径 为 0.1detaL=2*pi*r/N%每 个 小 段 的 长 度 %计 算 l矩 阵 各 个 元 素D=1.781;for wx=1:N % this is for fufill the matrixfor wy=1:Nx(wx)=r*cos(pi/N+(wx-1)*2*pi)/N);x(wy)=r*cos(pi/N+(wy-1)*2*pi)/N);y(wx)=r*sin(pi/N+(wx-1)*2*pi)/N);y(wy)=r*sin(pi/N+(wy-1)*2*pi)/N);if wx=wyl(wx,wy)=(k*yita*detaL/4)*(1-j*2*(log(D*k*detaL/4)-1);else l(wx,wy)=(k*yita/4)*besselh(0,2,k*sqrt(x(wx)-x(wy)2+(y(wx)-y(wy)2)*detaL;l(wy,wx)=(k*yita/4)*besselh(0,2,k*sqrt(x(wx)-x(wy)2+(y(wx)-y(wy)2)*detaL;endendend;%给 g矩 阵 赋 值%入 射 波 以 0度 (水 平 )入 射 圆 柱 导 体 for px=1:Ng(px)=exp(j*k*(x(px)*cos(0)+y(px)*sin(0);end%求 解 电 流 展 开 系 数afa=g*inv(l);%求 解 面 电 流 密 度for py=1:N J(py)=afa(py);magJ(py)=abs(J(py);phasJ(py)=angle(J(py);endsubplot(2,2,1);polar(pi/N:2*pi/N:2*pi,magJ);title(面 电 流 密 度 幅 度 分 布 );subplot(2,2,2);plot(phasJ);title(面 电 流 密 度 相 位 分 布 ); %求 解 远 区 散 射 场fai=pi/N:2*pi/N:2*pi;p=1000%所 求 远 区 场 与 柱 体 的 距 离Es=0;for pj=1:NEs=Es+(k*yita*exp(-j*k*p+0.75*pi)/sqrt(j*8*k*p)*afa(pj)*detaL*exp(-j*k*(x(pj)*cos(fai)+y(pj)*sin(fai);endsubplot(2,2,3); polar(fai,abs(Es)/max(abs(Es);title(电 场 幅 度 随 离 开 柱 体 距 离 的 变 化 趋 势 );subplot(2,2,4);plot(1:N,angle(Es);title(电 场 相 位 随 角 度 变 化 趋 势 );程 序 结 果 ： 1、 时 域 有 限 差 分 法问 题 ： 在 TM模 式 下 ， 编 写 matlab 程 序 ， 使 用 FDTD方 法 分 析 点 脉 冲 源 激 励 情况 下 ， 正 方 形 MUR边 界 的 性 能 。分 析 ： 对 TM波 ， FDTD公 式 为 ：1 12 2 ( , 1) ( , )1 1( , ) ( ) ( , ) ( )2 2 n nn n z zx x E i j E i jH i j cp m H i j cq m y 1 12 2 ( 1, ) ( , )1 1( , ) ( ) ( , ) ( )2 2 n nn n z zy y E i j E i jH i j cp m H i j cq m x 1 1 2 21 1 1( , ) ( , )2 2( , ) ( ) ( , ) ( ) n ny yn nz z H i j H i jE i j ca m E i j cb m x 1 12 21 1( , ) ( , )2 2n nx xH i j H i jy 其 中 ， 真 空 中 ca(m)=cp(m)=1; 0( ) tcb m ; 0( ) tcq m 二 维 Mur吸 收 边 界 条 件 的 FDTD形 式 如 下 （ 一 阶 ） ：1 1( , ) ( 1, ) ( 1, ) ( , )n n n nz z z zc tE i j E i j E i j E i jc t （ 左 边 界 ）1 1( , ) ( 1, ) ( 1, ) ( , )n n n nz z z zc tE i j E i j E i j E i jc t （ 右 边 界 ）1 1( , ) ( , 1) ( , 1) ( , )n n n nz z z zc tE i j E i j E i j E i jc t （ 上 边 界 ）1 1( , ) ( , 1) ( , 1) ( , )n n n nz z z zc tE i j E i j E i j E i jc t （ 下 边 界 ）二 维 Mur吸 收 边 界 角 点 的 处 理 如 下 ：1 12( , ) ( 1, 1) ( 1, 1) ( , )2n n n n z z z zc tE i j E i j E i j E i jc t （ 左 下 角 ）1 12( , ) ( 1, 1) ( 1, 1) ( , )2n n n nz z z zc tE i j E i j E i j E i jc t （ 右 上 角 ）1 12( , ) ( 1, 1) ( 1, 1) ( , )2n n n nz z z zc tE i j E i j E i j E i jc t （ 左 上 角 ）1 12( , ) ( 1, 1) ( 1, 1) ( , )2n n n nz z z zc tE i j E i j E i j E i jc t （ 右 下 角 ）Matlab 程 序 如 下 ： yps=8.85*10-12;%真 空 中 介 电 常 数u=4*pi*10-7;c=3*108;%真 空 中 光 速L=0.1;%TM波 约 定 ：%Hx(i,j)=Hx(i*deltaX,(j+1/2)*deltaY);%Hy(i,j)=Hy(i+1/2)*deltaX,j*deltaY);%Ez(i,j)=Ez(i*deltaX,j*deltaY);%区 域 设 定 N=10%将 正 方 形 区 域 均 分 成 为 N*N 个 小 正 方 区 域%初 始 电 磁 场 值%旧 时 刻 场 量 值Ez1=zeros(N+1,N+1); Hx1=zeros(N+1,N);Hy1=zeros(N,N+1);%迭 代 后 新 时 刻 场 量 值Ez2=zeros(N+1,N+1);Hx2=zeros(N+1,N);Hy2=zeros(N,N+1);%设 置 脉 冲 源tao=2*10-10;%高 斯 脉 冲 的 宽 度t0=tao/2%高 斯 脉 冲 峰 值 出 现 的 时 刻 deltaL=L/N;%空 间 间 隔 要 求deltaT=tao/(4*N);%时 间 间 隔 要 求t=0;%可 以 修 改 t的 值 一 观 察 不 同 时 刻 的 电 磁 场%迭 代 过 程fordc=1:1000%迭 代 的 次 数ca=1;%迭 代 系 数cb=deltaT/yps;cp=1;cq=deltaT/u; %激 励 源 加 入It=exp(-4*pi*(t-t0).2/(tao2);%非 边 界 场 值 迭 代clearij;fori=1:(N+1)forj=1:NHx2(i,j)=cp*Hx1(i,j)-cq*(Ez1(i,j+1)-Ez1(i,j)/deltaL;endend clearij;fori=1:Nforj=1:(N+1)Hy2(i,j)=cp*Hy1(i,j)+cq*(Ez1(i+1,j)-Ez1(i,j)/deltaL;endendclearij;fori=2:Nforj=2:N Ez2(i,j)=ca*Ez1(i,j)+cb*(Hy2(i,j)-Hy2(i-1,j)/deltaL-(Hx2(i,j)-Hx2(i,j-1)/deltaL);endendclearij;i=8;j=8;%可 以 设 置 激 励 源 的 位 置Ez2(i,j)=Ez1(i,j)+It;%激 励 加 入%边 界 非 角 点 场 值 迭 代%左 边 界 的 迭 代 clearij;i=1;forj=2:NEz2(i,j)=Ez1(i+1,j)+(c*deltaT-deltaL)*(Ez2(i+1,j)-Ez1(i,j)/(c*deltaT+deltaL);end%右 边 界 的 迭 代clearij;i=N+1;forj=2:NEz2(i,j)=Ez1(i-1,j)+(c*deltaT-deltaL)*(Ez2(i-1,j)-Ez1(i,j)/(c*deltaT+deltaL); end%上 边 界 的 迭 代clearij;j=N+1;fori=2:NEz2(i,j)=Ez1(i,j-1)+(c*deltaT-deltaL)*(Ez2(i,j-1)-Ez1(i,j)/(c*deltaT+deltaL);end%下 边 界 的 迭 代clearij; j=1;fori=2:NEz2(i,j)=Ez1(i,j+1)+(c*deltaT-deltaL)*(Ez2(i,j+1)-Ez1(i,j)/(c*deltaT+deltaL);end%角 点 的 迭 代clearij;Ez2(1,1)=Ez1(2,2)+(c*deltaT-sqrt(2)*deltaL)*(Ez2(2,2)-Ez1(1,1)/(c*deltaT+sqrt(2)*deltaL);Ez2(N+1,N+1)=Ez1(N,N)+(c*deltaT-sqrt(2)*deltaL)*(Ez2(N,N)-Ez1(N+1,N+1)/(c*deltaT+sqrt(2)*deltaL);Ez2(1,N+1)=Ez1(2,N)+(c*deltaT-sqrt(2)*deltaL)*(Ez2(2,N)-Ez1(1,N+1)/(c*deltaT+sqrt(2)*deltaL); Ez2(N+1,1)=Ez1(N,2)+(c*deltaT-sqrt(2)*deltaL)*(Ez2(N,2)-Ez1(N+1,1)/(c*deltaT+sqrt(2)*deltaL);%一 次 迭 代 结 束 ， 新 得 到 的 值 变 为 下 次 迭 代 的 旧 值clearij;fori=1:(N+1);forj=1:(N+1)Ez1(i,j)=Ez2(i,j);endendclearij; fori=1:(N+1)forj=1:NHx1(i,j)=Hx2(i,j);endendclearij;fori=1:Nforj=1:(N+1)Hy1(i,j)=Hy2(i,j);end end%一 次 迭 代 结 束end%总 的 迭 代 结 束clearij;i=1:(N+1);j=1:(N+1); figure(1)contour(i,j,Ez2(i,j);结 果 如 下 :