CODEWITHSUNDEEP
matlabnumerical-methods
Sean James Jamieson
Sean James Jamieson

Reputation: 101

Finite Difference Time Domain (FTDT) method for 1D EM Wave

I have attempted to write a code in order to solve the following coupled partial differential EM wave equations:

enter image description here

The code employs finite difference time domain using the Yee algorithm which can be read about in the following two online documents:

http://www.eecs.wsu.edu/~schneidj/ufdtd/ufdtd.pdf

http://www.eecs.wsu.edu/~schneidj/ufdtd/chap3.pdf

I want my source at the left hand side boundary to be a sinusoidal wave with parameters as such: 

Ex(0,t) = E0 sin(2πft)      for 0 ≤ t ≤ 1/f

The code updates the electric and magnetic field properties of the wave with each loop.

My initial code is as follows:

%FDTD Yee algorithm to solve coupled EM wave equations
clear
clc


G=50;           %Specify size of the grid
f=10^3;         %choose initial frequency of wave in hertz
e=1;            %specify permitivity and permeability (normalised condition)
u=1;            
Nt=150;         %time steps
E0=1;           %Electric Field initial amplitude

%Specify step sizes using corruant condition
c=3*10^8;
dx=0.01;
dt=dx/2*c;

%make constant terms
c1=-dt/(dx*e);
c2=-dt/(dx*u);

%create vgector place holders
Ex=zeros(1,G);
Hy=zeros(1,G);

%create updating loop
M=moviein(Nt);
for t=1:Nt

    % Spatial Ex

    for k=2:G-1
        Ex(k)=Ex(k)+c1*(Hy(k)-Hy(k-1));
    end
    Ex(G)=0; %PEC boundary condition
    %E Source at LHS boundary 
    Ex(1)=E0*sin(2*pi*f*t);
    %Spatial Hy
    for n=1:G-1
        Hy(n)=Hy(n)+c2*(Ex(n)-Ex(n+1));
    end
    Hy(G)=0; %PMC boundary condition

plot(Ex);
M(:,t) = getframe;
end
movie(M,1);

Basically I want to create an updating movie which shows the sinusoidal wave propagating to the right hand side boundary coded as a perfect electrical conductor; therefore reflecting the wave, and then propagating back to the left hand side boundary coded as a perfect insulator; absorbing the wave.

The problems I have are as follows:

1) I'm not sure how to properly implement my desired source. It don't appear to be purely sinusoidal.

2) The wave I've coded begins to propagate but it quickly disappears for the majority of the simulation. I do not know why this is occurring

3) I do not know much about running a movie simulation and the plot oscillates as the solution is being run. How can I stop this?

Upvotes: 1

Views: 3365

Answers (2)

Paul Mackenzie
Paul Mackenzie

Reputation: 11

You should set the Courant number closer to 1 say 0.995. Thus delta_t = 0.995*delta_x/c. Also assuming delta_x is in METRIC units then e and u should be in metric units. I do not know about the specific coding language used but in c or c++ there is no need for intermediate variable Ey1 etc. Also there should be at least 10 samples per wavelength for accuracy ( preferably 60). Thus wavelength = 60*delta_x and thus the frequency equals roughly of the order 10 to power of 9. Also, I think the sinesoidal source should be E0 * sin(2* pi * f * t * delta_t). You need to adjust your constants, and try it again

Upvotes: 0

user2271770
user2271770

Reputation:

Your wave attenuates because the diference equations are not correctly implemented; instead:

Ex(k)=Ex(k)+c1*(Hy(k)-Hy(k-1));

you should use

Ex1(k)=Ex(k)+c1*(Hy(k)-Hy(k-1));

and instead of:

Hy(n)=Hy(n)+c2*(Ex(n)-Ex(n+1));

you should use:

Hy1(n)=Hy(n)+c2*(Ex(n)-Ex(n+1));

and, in the end of the loop update the current "dataframe":

Hy = Hy1;
Ey = Ey1;

(you should take care also the boundary conditions).

If you want a fixed plot frame that doesn't change when your data changes, create first a axis where you can plot into, with a defined xmin/max and ymin/max, see http://www.mathworks.com/help/matlab/ref/axis.html

Upvotes: 1

Related Questions