时域有限差分法

Problem 5.1

In this illustrative solution, the electric-field hard source condition of (5.1) is implemented at the far-left grid boundary. The source time function has an amplitude of 1.0 V/m and a frequency of 10 GHz. The reflecting barrier (PEC) is implemented at the far-right grid boundary. The computational domain represents a physical length of 15 cm.

Matlab code:

%***********************************************************************

% 1D FINITE-DIFFERENCE TIME-DOMAIN SOLUTION: PLANE WAVE PROPAGATION

%***********************************************************************

%

% Program author:

% Prof. Susan C. Hagness

% Department of Electrical and Computer Engineering

% University of Wisconsin-Madison

% 1415 Engineering Drive

% Madison, WI 53706-1691

% hagness@

%

%***********************************************************************

clear;

%..........Material Parameters............

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

eps=[1.0]; %relative permittivity

sig=[0.0]; %electric conductivity

mur=[1.0]; %relative permeability

sim=[0.0]; %magnetic loss

media=length(eps);

%..........Space, Time, and Source Parameters...

S=1.0;

freq=10e9; %frequency of sinusoidal excitation = 10 GHz

E0=1.0; %amplitude of sinusoidal excitation = 1.0 V/m

lambda=cc/freq;

length=0.15; %physical length of grid (in units of m)

dx=lambda/20; %grid resolution of 20 cells per wavelength

dt=S*dx/cc;

ie=round(length/dx)+1; %number of Ez samples in grid

ih=ie-1; %number of Hy samples in grid

nmax=3*round(ie*S);

source(1:nmax)=E0*sin(2*pi*freq*(1:nmax)*dt);

%..........Initial Conditions...........

ez(1:ie)=0.0;

hy(1:ih)=0.0;

Problem 5.1 (cont.)

%..........Update Coefficients.........

for i=1:media

eprop=sig(i)*dt/(2.0*epsz*eps(i));

ca(i)=(1.0-eprop)/(1.0+eprop);

cb(i)=dt/(epsz*eps(i)*dx)/(1.0+eprop);

hprop=sim(i)*dt/(2.0*muz*mur(i));

da(i)=(1.0-hprop)/(1.0+hprop);

db(i)=dt/(muz*mur(i)*dx)/(1.0+hprop);

end

mediaez(1:ie)=1; %media pointers

mediahy(1:ih)=1;

%..........Time-Stepping..............

for n=1:nmax

%%%%%%%%%%% EZ UPDATES %%%%%%%%%%%

ez(2:ih)=ca(mediaez(2:ih)).*ez(2:ih)+cb(mediaez(2:ih)).*(hy(2:ih)-hy(1:ih-1)); ez(1)=source(n); % incident wave source condition (hard source) %%%%%%%%%%% HY UPDATES %%%%%%%%%%%

hy(1:ih)=da(mediahy(1:ih)).*hy(1:ih)+db(mediahy(1:ih)).*(ez(2:ie)-ez(1:ih));

%%%%%%%%%%% DATA OUTPUT %%%%%%%%%%%

subplot(2,1,1),plot(ez,'r'),axis([1 ie -3 3]);

ylabel('Ez'); title(['n=',num2str(n)]);

subplot(2,1,2),plot(hy,'b'),axis([1 ih -8e-3 8e-3]);

ylabel('Hy'); xlabel('grid coordinate')

pause(0.05)

end