1D FDTD plane wave propagation simulation
42 views (last 30 days)
Show older comments
I'm currently trying to simulate a simple case of wave propagation in free space before adding in more complexities, and already I'm stumped. Most examples I see online set the Courant stability condition to 1, i.e. the wave travels 1 spatial step in 1 time step. However, I set mine to 2 with the intention of increasing it further should I require more resolution.
Problem: My source is a unity amplitude sine wave (in exponential form), and when I plot it the amplitude is around 2 with increasing ringing with increasing distance.
I'm using these slides as my reference (https://empossible.net/wp-content/uploads/2020/01/Lecture-Formulation-of-1D-FDTD.pdf). The relevant slides are on p.12 and 20.
c0=299792458; %speed of light
u0 = 1.256637e-6;
e0 = 8.85418782e-12;
n1 = 1.; er1=n1^2; %refractive index and relative permittivity
fmax = 1e6; %frequency
lambda_min = c0/(fmax*n1); %minimum wavelength
z_step = lambda_min/(100); %step size in space
t_step = z_step/(2*c0); %Courant stability condition
t_total = t_step*4000; %total runtime
T = t_total/t_step; %total number of timesteps
Z = 400; %total number of spatial steps
Zarray = linspace(1,Z,Z);
Hx = zeros(Z,1); Ey = zeros(Z,1);
C_Ey = ones(Z,1).*c0*t_step/z_step/er1;
C_Hx = ones(Z,1).*c0*t_step/z_step;
H1=0;H2=0;E1=0;E2=0;
for t_count = 1:T
disp(t_count)
H2=H1; H1=Hx(1);
for z_count = 1:Z-1
Hx(z_count) = Hx(z_count) + C_Hx(z_count)*(Ey(z_count+1)-Ey(z_count));
end
Hx(Z) = Hx(Z) + C_Hx(Z)*(E1-Ey(Z));
E2=E1; E1=Ey(Z);
Ey(1) = Ey(1) + C_Ey(1)*(Hx(1)-H1);
for z_count = 2:Z
Ey(z_count) = Ey(z_count) + C_Ey(z_count)*(Hx(z_count)-Hx(z_count-1));
end
Ey(50) = Ey(50)+exp(1i*2*pi*fmax*t_count*t_step); %soft source
plot(Zarray,real(Ey))
drawnow
end
0 Comments
Answers (1)
Romain
on 12 Jul 2024
Edited: Romain
on 12 Jul 2024
Hi Jerry,
At line 33, change:
Ey(50) = Ey(50)+exp(1i*2*pi*fmax*t_count*t_step); %soft source
To:
Ey(1) = Ey(1)+exp(1i*2*pi*fmax*t_count*t_step); %soft source
And you get a much smoother wave:
Additionally, a small error is present in the definition of lambda_min: replace n2 with n1 or define n2 before.
If you sorely need the wave to generate at abscissa 50, do not perform the above modification but change the for range at lines 22 and 29:
for z_count = 50:Z-1
Hx(z_count) = Hx(z_count) + C_Hx(z_count)*(Ey(z_count+1)-Ey(z_count));
end
for z_count = 51:Z
Ey(z_count) = Ey(z_count) + C_Ey(z_count)*(Hx(z_count)-Hx(z_count-1));
end
And you will get:
2 Comments
Romain
on 15 Jul 2024
Hi Jerry,
Is there a reason why you generate your sine wave in exponential form ?
How do you deal with boundaries Z=50 and Z=400 ?
Hx(Z) = Hx(Z) + C_Hx(Z)*(E1-Ey(Z));
But you set a each iteration
E1 = Ey(Z)
So Hx(Z) is never modified. Same goes for Ey(50), where you set (I deliberately changed indexes from 1 to 50):
H1 = Hx(50)
And then
Ey(50) = Ey(50) + C_Ey(50)*(Hx(50)-H1)
Looking at the last slide of the pdf you mentioned:
Should not Ey and Hx be equal to 0 outside of the boundaries ? So E1=0 and H1=0 at any time ?
If you set it to 0, I get the reflecting behavior you asked for
See Also
Categories
Find more on General Physics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!