1130阅读
0回复

[求助]为什么高斯波（软源）在某个步数后在加载点变为恒值（不为0）

UID ：58939

注册：2010-05-07
登录：2012-06-17
发帖：31
等级：仿真新人

0楼发表于: 2010-07-14 21:29:06

麻烦问下，我在改Susan C. Hagness的3-D FDTD code with PEC boundaries，将所加的激励源改了一下后，怎么在60多步后加载点变为恒值，电场不在变化，不为0.
印象中一定步数后高斯波就应该没了呀，怎么一直在呀？
程序如下：
clear

%***********************************************************************
%     Fundamental constants
%***********************************************************************

cc=2.99792458e8;            %speed of light in free space
muz=4.0*pi*1.0e-7;          %permeability of free space
epsz=1.0/(cc*cc*muz);       %permittivity of free space

%***********************************************************************
%     Grid parameters
%***********************************************************************

ie=50;    %number of grid cells and Ex samples along x-direction
je=24;    %number of grid cells and Ey samples along y-direction
ke=10;    %number of grid cells and Ez samples along z-direction

ib=ie+1;    %grid boundary coordinates
jb=je+1;
kb=ke+1;

is=26;       %location of z-directed current source
js=13;       %location of z-directed current source

%kobs=5;
ks=6;

dx=0.002;          %space increment of cubic lattice
dt=dx/(2.0*cc);    %time step

nmax=200;          %total number of time steps

%***********************************************************************
%     Differentiated Gaussian pulse excitation
%***********************************************************************

rtau=50.0e-12;
tau=rtau/dt;
ndelay=3*tau;
J0=-1.0;

%***********************************************************************
%     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(ie,jb,kb);
ey=zeros(ib,je,kb);
ez=zeros(ib,jb,ke);
hx=zeros(ib,je,ke);
hy=zeros(ie,jb,ke);
hz=zeros(ie,je,kb);

%figure
%set(gcf,'DoubleBuffer','on')

%***********************************************************************
%     BEGIN TIME-STEPPING LOOP
%***********************************************************************

for n=1:nmax

%***********************************************************************
%     Update electric fields
%***********************************************************************

ex(1:ie,2:je,2:ke)=ca*ex(1:ie,2:je,2:ke)+...
                   cb*(hz(1:ie,2:je,2:ke)-hz(1:ie,1:je-1,2:ke)+...
                       hy(1:ie,2:je,1:ke-1)-hy(1:ie,2:je,2:ke));

ey(2:ie,1:je,2:ke)=ca*ey(2:ie,1:je,2:ke)+...
                   cb*(hx(2:ie,1:je,2:ke)-hx(2:ie,1:je,1:ke-1)+...
                       hz(1:ie-1,1:je,2:ke)-hz(2:ie,1:je,2:ke));

ez(2:ie,2:je,1:ke)=ca*ez(2:ie,2:je,1:ke)+...
                   cb*(hx(2:ie,1:je-1,1:ke)-hx(2:ie,2:je,1:ke)+...
                       hy(2:ie,2:je,1:ke)-hy(1:ie-1,2:je,1:ke));

%ez(is,js,ks)=ez(is,js,ks)+...
              % J0*exp(-((n-ndelay)^2/tau^2));

ez(is,js,ks)=ez(is,js,ks)+J0*exp(-0.5*((n-25)^2/20^2));
%***********************************************************************
%     Update magnetic fields
%***********************************************************************

hx(2:ie,1:je,1:ke)=hx(2:ie,1:je,1:ke)+...
                   db*(ey(2:ie,1:je,2:kb)-ey(2:ie,1:je,1:ke)+...
                       ez(2:ie,1:je,1:ke)-ez(2:ie,2:jb,1:ke));

hy(1:ie,2:je,1:ke)=hy(1:ie,2:je,1:ke)+...
                   db*(ex(1:ie,2:je,1:ke)-ex(1:ie,2:je,2:kb)+...
                       ez(2:ib,2:je,1:ke)-ez(1:ie,2:je,1:ke));

hz(1:ie,1:je,2:ke)=hz(1:ie,1:je,2:ke)+...
                   db*(ex(1:ie,2:jb,2:ke)-ex(1:ie,1:je,2:ke)+...
                       ey(1:ie,1:je,2:ke)-ey(2:ib,1:je,2:ke));

%***********************************************************************
%     Visualize fields
%***********************************************************************

timestep=int2str(n);
tview(:,:)=ez(:,:,ks);
sview(:,:)=ez(:,js,:);

subplot('position',[0.15 0.45 0.7 0.45]),imagesc(tview');
shading flat;
caxis([-1.0 1.0]);
colorbar;
axis image; axis xy;
title(['Ez(i,j,k=5), ti ..

未注册仅能浏览部分内容，查看全部内容及附件请先登录或注册

发帖回复