第1页 / 共12页
第2页 / 共12页
第3页 / 共12页
第4页 / 共12页
第5页 / 共12页
第6页 / 共12页
第7页 / 共12页
第8页 / 共12页
FTDT算法 二维三维程序，无错可直接运行.doc
3-D FDTD 二维和三维
Program author: Susan C. Hagness
Date of this version: February 2000
This MATLAB M-file implements the finite-difference time-domain
solution of Maxwell's curl equations over a three-dimensional
Cartesian space lattice comprised of uniform cubic grid cells.
To illustrate the algorithm, an air-filled rectangular cavity
resonator is modeled. The length, width, and height of the
cavity are 10.0 cm (x-direction), 4.8 cm (y-direction), and
2.0 cm (z-direction), respectively.【矩形谐振腔 10 * 4.8 * 2.0 cm】
Department of Electrical and Computer Engineering
University of Wisconsin-Madison
1415 Engineering Drive
Madison, WI 53706-1691
608-265-5739
hagness@engr.wisc.edu
%***********************************************************************
%
%***********************************************************************
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
To execute this M-file, type "fdtd3D" at the MATLAB prompt.
This M-file displays the FDTD-computed Ez fields at every other
time step, and records those frames in a movie matrix, M, which
is played at the end of the simulation using the "movie" command.
【显示 Ez 的场分布】
where tau=50 ps. The FWHM spectral bandwidth of this zero-dc-
content pulse is approximately 7 GHz. The grid resolution
(dx = 2 mm) was chosen to provide at least 10 samples per
wavelength up through 15 GHz.【为了一个波长 10 个采样点，网格分割 2mm】
The computational domain is truncated using PEC boundary
conditions:
ex(i,j,k)=0 on the j=1, j=YB, k=1, and k=ZB planes
ey(i,j,k)=0 on the i=1, i=XB, k=1, and k=ZB planes
ez(i,j,k)=0 on the i=1, i=XB, j=1, and j=YB planes
These PEC boundaries form the outer lossless walls of the cavity.
The cavity is excited by an additive current source oriented
along the z-direction. The source waveform is a differentiated
Gaussian pulse given by 【在腔体的一端加入高斯激励源】
J(t)=-J0*(t-t0)*exp(-(t-t0)^2/tau^2),
%***********************************************************************
clear
%***********************************************************************
%
%***********************************************************************
Fundamental constants
cc=2.99792458e8;
muz=4.0*pi*1.0e-7;
epsz=1.0/(cc*cc*muz);
%speed of light in free space【光速】
%permeability of free space【磁导率】
%permittivity of free space【介电常数】
%***********************************************************************
%
%***********************************************************************
Grid parameters
XGridNum=50;
YGridNum=24;
ZGridNum=10;
%number of grid cells in x-direction【x 方向单元格数目】
%number of grid cells in y-direction【y 方向单元格数目】
%number of grid cells in z-direction【z 方向单元格数目】
XB=XGridNum+1;
YB=YGridNum+1;
ZB=ZGridNum+1;
%location of z-directed current source 【z 方向电流激励源的位置 x】
%location of z-directed current source 【z 方向电流激励源的位置 y】
XS=26;
YS=13;
kobs=5;
dx=0.002;
dt=dx/(2.0*cc);
%space increment of cubic lattice【空间立方格步进】
%time step【时间步进】
nmax=500;
%total number of time steps【总步数】
Differentiated Gaussian pulse excitation【不同的高斯脉冲激励】
%***********************************************************************
%
%
%***********************************************************************
J(t)=-J0*(t-t0)*exp(-(t-t0)^2/tau^2)
rtau=50.0e-12;
tau=rtau/dt;
ndelay=3*tau;
srcconst=-dt*3.0e+11;
%***********************************************************************
%
%***********************************************************************
Material parameters【材料参数】
eps=1.0;
sig=0.0;
%***********************************************************************
%
%***********************************************************************
Updating coefficients【更新系数】
ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps));
cb=(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps));
da=1.0;
db=dt/muz/dx;
%***********************************************************************
%
%***********************************************************************
Field arrays
ex=zeros(XGridNum,YB,ZB);
ey=zeros(XB,YGridNum,ZB);
ez=zeros(XB,YB,ZGridNum);
hx=zeros(XB,YGridNum,ZGridNum);
hy=zeros(XGridNum,YB,ZGridNum);
hz=zeros(XGridNum,YGridNum,ZB);
%***********************************************************************
%
%***********************************************************************
Movie initialization
tview(:,:)=ez(:,:,kobs);
sview(:,:)=ez(:,YS,:);
subplot('position',[0.15 0.45 0.7 0.45]),pcolor(tview');
shading flat;
caxis([-1.0 1.0]);
colorbar;%【显示右侧颜色示意条】
axis image;%【添加图形坐标】？
title(['Ez(i,j,k=5), time step = 0']);%【添加图形标题】
xlabel('i coordinate');%【添加 i 坐标显示】
ylabel('j coordinate');%【添加 j 坐标显示】
subplot('position',[0.15 0.10 0.7 0.25]),pcolor(sview');
shading flat;%【删除网格】
caxis([-1.0 1.0]);
colorbar;
axis image;
title(['Ez(i,j=13,k), time step = 0']);
xlabel('i coordinate');
ylabel('k coordinate');
rect=get(gcf,'Position');
rect(1:2)=[0 0];
M=moviein(nmax/2,gcf,rect);
%***********************************************************************
%
%***********************************************************************
BEGIN TIME-STEPPING LOOP【时域循环】
for n=1:nmax
%***********************************************************************
%
%***********************************************************************
Update electric fields【更新电场】
ex(1:XGridNum,2:YGridNum,2:ZGridNum)=ca*ex(1:XGridNum,2:YGridNum,2:ZGridNum)+...
cb*(hz(1:XGridNum,2:YGridNum,2:ZGridNum)-hz(1:XGridNum,1:YGridNum-1,2:ZGridNum)+.
..
hy(1:XGridNum,2:YGridNum,1:ZGridNum-1)-hy(1:XGridNum,2:YGridNum,2:ZGridNum));
ey(2:XGridNum,1:YGridNum,2:ZGridNum)=ca*ey(2:XGridNum,1:YGridNum,2:ZGridNum)+...
cb*(hx(2:XGridNum,1:YGridNum,2:ZGridNum)-hx(2:XGridNum,1:YGridNum,1:ZGridNum-1)
+...
hz(1:XGridNum-1,1:YGridNum,2:ZGridNum)-hz(2:XGridNum,1:YGridNum,2:ZGridNum));
ez(2:XGridNum,2:YGridNum,1:ZGridNum)=ca*ez(2:XGridNum,2:YGridNum,1:ZGridNum)+...
cb*(hx(2:XGridNum,1:YGridNum-1,1:ZGridNum)-hx(2:XGridNum,2:YGridNum,1:ZGridNum)
+...
hy(2:XGridNum,2:YGridNum,1:ZGridNum)-hy(1:XGridNum-1,2:YGridNum,1:ZGridNum));
ez(XS,YS,1:ZGridNum)=ez(XS,YS,1:ZGridNum)+...
srcconst*(n-ndelay)*exp(-((n-ndelay)^2/tau^2));
fprintf('\n\n sourse = %d \n',srcconst*(n-ndelay)*exp(-((n-ndelay)^2/tau^2)));
for tt=1:ZGridNum
fprintf('%d \t',ez(24,10,tt));
end
%***********************************************************************
%
%***********************************************************************
Update magnetic fields【更新磁场】
hx(2:XGridNum,1:YGridNum,1:ZGridNum)=hx(2:XGridNum,1:YGridNum,1:ZGridNum)+...
db*(ey(2:XGridNum,1:YGridNum,2:ZB)-ey(2:XGridNum,1:YGridNum,1:ZGridNum)+...
ez(2:XGridNum,1:YGridNum,1:ZGridNum)-ez(2:XGridNum,2:YB,1:ZGridNum));
hy(1:XGridNum,2:YGridNum,1:ZGridNum)=hy(1:XGridNum,2:YGridNum,1:ZGridNum)+...
db*(ex(1:XGridNum,2:YGridNum,1:ZGridNum)-ex(1:XGridNum,2:YGridNum,2:ZB)+...
ez(2:XB,2:YGridNum,1:ZGridNum)-ez(1:XGridNum,2:YGridNum,1:ZGridNum));
hz(1:XGridNum,1:YGridNum,2:ZGridNum)=hz(1:XGridNum,1:YGridNum,2:ZGridNum)+...
db*(ex(1:XGridNum,2:YB,2:ZGridNum)-ex(1:XGridNum,1:YGridNum,2:ZGridNum)+...
ey(1:XGridNum,1:YGridNum,2:ZGridNum)-ey(2:XB,1:YGridNum,2:ZGridNum));
%【注意】本程序没有吸收边界条件；没有数值稳定性分析；
%***********************************************************************
%
%***********************************************************************
Visualize fields【观察场分布】
if mod(n,2)==0;
timestep=int2str(n);
tview(:,:)=ez(:,:,kobs);
sview(:,:)=ez(:,YS,:);
TT=ez(5,5,5);
subplot('position',[0.15 0.45 0.7 0.45]),pcolor(tview');
shading flat;
caxis([-1.0 1.0]);
colorbar;
axis image;
title(['Ez(i,j,k=5), time step = ',timestep]);
xlabel('i coordinate');
ylabel('j coordinate');
subplot('position',[0.15 0.10 0.7 0.25]),pcolor(sview');
shading flat;
caxis([-1.0 1.0]);
colorbar;
axis image;
title(['Ez(i,j=13,k), time step = ',timestep]);
xlabel('i coordinate');
ylabel('k coordinate');
nn=n/2;
M(:,nn)=getframe(gcf,rect);
end;
%***********************************************************************
%
%***********************************************************************
END TIME-STEPPING LOOP
end
3-D FDTD code with PEC boundaries
Program author: Susan C. Hagness
movie(gcf,M,0,10,rect);
%***********************************************************************
%
%***********************************************************************
%
%
%
%
%
%
Department of Electrical and Computer Engineering
University of Wisconsin-Madison
1415 Engineering Drive
Madison, WI 53706-1691
608-265-5739
hagness@engr.wisc.edu
Date of this version: February 2000
This MATLAB M-file implements the finite-difference time-domain
solution of Maxwell's curl equations over a three-dimensional
Cartesian space lattice comprised of uniform cubic grid cells.
The computational domain is truncated using PEC boundary
conditions:
ex(i,j,k)=0 on the j=1, j=YB, k=1, and k=ZB planes
ey(i,j,k)=0 on the i=1, i=XB, k=1, and k=ZB planes
ez(i,j,k)=0 on the i=1, i=XB, j=1, and j=YB planes
These PEC boundaries form the outer lossless walls of the cavity.
To illustrate the algorithm, an air-filled rectangular cavity
resonator is modeled. The length, width, and height of the
cavity are 10.0 cm (x-direction), 4.8 cm (y-direction), and
2.0 cm (z-direction), respectively.【矩形谐振腔 10 * 4.8 * 2.0 cm】
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%***********************************************************************
To execute this M-file, type "fdtd3D" at the MATLAB prompt.
This M-file displays the FDTD-computed Ez fields at every other
time step, and records those frames in a movie matrix, M, which
is played at the end of the simulation using the "movie" command.
【显示 Ez 的场分布】
where tau=50 ps. The FWHM spectral bandwidth of this zero-dc-
content pulse is approximately 7 GHz. The grid resolution
(dx = 2 mm) was chosen to provide at least 10 samples per
wavelength up through 15 GHz.【为了一个波长 10 个采样点，网格分割 2mm】
The cavity is excited by an additive current source oriented
along the z-direction. The source waveform is a differentiated
Gaussian pulse given by 【在腔体的一端加入高斯激励源】
J(t)=-J0*(t-t0)*exp(-(t-t0)^2/tau^2),
clear
%***********************************************************************
%
%***********************************************************************
Fundamental constants
cc=2.99792458e8;
%speed of light in free space【光速】
muz=4.0*pi*1.0e-7;
epsz=1.0/(cc*cc*muz);
%permeability of free space【磁导率】
%permittivity of free space【介电常数】
%***********************************************************************
%
%***********************************************************************
Grid parameters
XGridNum=50;
YGridNum=24;
ZGridNum=10;
%number of grid cells in x-direction【x 方向单元格数目】
%number of grid cells in y-direction【y 方向单元格数目】
%number of grid cells in z-direction【z 方向单元格数目】
XB=XGridNum+1;
YB=YGridNum+1;
ZB=ZGridNum+1;
%location of z-directed current source 【z 方向电流激励源的位置 x】
%location of z-directed current source 【z 方向电流激励源的位置 y】
XS=26;
YS=13;
kobs=5;
dx=0.002;
dt=dx/(2.0*cc);
%space increment of cubic lattice【空间立方格步进】
%time step【时间步进】
nmax=500;
%total number of time steps【总步数】
Differentiated Gaussian pulse excitation【不同的高斯脉冲激励】
%***********************************************************************
%
%
%***********************************************************************
J(t)=-J0*(t-t0)*exp(-(t-t0)^2/tau^2)
rtau=50.0e-12;
tau=rtau/dt;
ndelay=3*tau;
srcconst=-dt*3.0e+11;
%***********************************************************************
%
%***********************************************************************
Material parameters【材料参数】
eps=1.0;
sig=0.0;
%***********************************************************************
%
Updating coefficients【更新系数】
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
