1D FDTD plane wave propagation simulation

57 次查看（过去 30 天）
Jerry Yeung
Jerry Yeung 2024-7-12
评论： Romain 2024-7-15
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

请先登录，再回答此问题。

回答（1 个）

Romain
Romain 2024-7-12
编辑：Romain 2024-7-12
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 个评论
Jerry Yeung
Jerry Yeung 2024-7-12
编辑：Jerry Yeung 2024-7-12
Thanks for the reply. Unfortunately this does not solve my issue. In your figures the wave has a peak magnitude of >20, whereas my source has a magnitude of unity. This is due to reflection from the left boundary, but you move the source anywhere within the grid and the rightward propagating wave still has a magnitude of 2. Furthermore, I require 1) my source to not be next to the boundaries and 2) simulation of the back-reflected (leftward-travelling) field later down the line. I can't simply not simulate anything before abscissa 50.
Romain
Romain 2024-7-15
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 General Physics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by