【计算电磁学】真空中电磁波传播的MATLAB一维FDTD仿真程序

模拟的空间总长1m，激励在0.5m处，两边为理想电导体边界（PEC），信号到达边界后反射，电场反向，反射系数为-1，磁场不反向，程序如下：

clc
close all
clear all

% Define initial constants
eps_0 = 8.854187817e-12;                  % permittivity of free space
mu_0  = 4*pi*1e-7;                        % permeability of free space
c     = 1/sqrt(mu_0*eps_0);               % speed of light
 
% Define problem geometry and parameters
domain_size = 1;                          % 1D problem space length in meters
dx = 1e-3;                                % cell size in meters, Δx=0.001m
dt = 3e-12;                               % duration of time step in seconds
number_of_time_steps = 2000;              % number of iterations
nx = round(domain_size/dx);               % number of cells in 1D problem space
source_position = 0.5;                    % position of the current source Jz

% Initialize field and material arrays
Ceze       = zeros(nx+1, 1);
Cezhy      = zeros(nx+1, 1);
Cezj       = zeros(nx+1, 1);
Ez         = zeros(nx+1, 1);
eps_r_z    = ones (nx+1, 1);              % free space
sigma_e_z  = zeros(nx+1, 1);              % free space
 
Chyh       = zeros(nx, 1);
Chyez      = zeros(nx, 1);
Chym       = zeros(nx, 1);
Hy         = zeros(nx, 1);
My         = zeros(nx, 1);
mu_r_y     = ones (nx, 1);                % free space
sigma_m_y  = zeros(nx, 1);                % free space
 
% Calculate FDTD updating coefficients
Ceze       = (2 * eps_r_z * eps_0 - dt * sigma_e_z) ...
           ./(2 * eps_r_z * eps_0 + dt * sigma_e_z);
 
Cezhy      = (2 * dt / dx) ...
           ./(2 * eps_r_z * eps_0 + dt * sigma_e_z);
       
 
Chyh       = (2 * mu_r_y * mu_0 - dt * sigma_m_y) ...
           ./(2 * mu_r_y * mu_0 + dt * sigma_m_y);
 
Chyez      = (2 * dt / dx) ...
           ./(2 * mu_r_y * mu_0 + dt * sigma_m_y);
 
 
% Define the Gaussian source waveform
time                  = dt * [0:number_of_time_steps-1].';
Ez_waveform           = exp(-((time-2e-10)/5e-11).^2)*1e-3/dx;
source_position_index = round(nx * source_position/domain_size)+1;

Ez_positions = [0:nx]*dx;
Hy_positions = ([0:nx-1]+0.5)*dx;

axis([0 1 -1 1]);
lez = line(Ez_positions, Ez*0, Ez, 'Color', 'b', 'linewidth', 1.5);
lhy = line(Hy_positions, 377*Hy, Hy*0, 'Color', 'r', 'LineWidth', 1.5, 'LineStyle','-.');
 
set(gca, 'fontsize', 12, 'fontweight', 'bold');
set(gcf,'Color','white');
axis square;

xlabel('x [m]');
ylabel('E [V/m]');
legend('Electric field', 'Magnetic field', 'location', 'northeast');

grid on;
 
% FDTD loop
for time_step = 1:number_of_time_steps
 
    % Update Ez for the current time step
    Ez(source_position_index) = Ez_waveform(time_step);
 
    % Update magnetic field
    Hy(1:nx) = Chyh(1:nx) .* Hy(1:nx) ...
         + Chyez(1:nx) .* (Ez(2:nx+1) - Ez(1:nx));
    % Update electric field
    Ez(2:nx) = Ceze (2:nx) .* Ez(2:nx) ...
             + Cezhy(2:nx) .* (Hy(2:nx) - Hy(1:nx-1));
              
    Ez(1)    = 0;       % Apply PEC boundary condition at x = 0 m
    Ez(nx+1) = 0;       % Apply PEC boundary condition at x = 1 m
 
    % Subroutine used to plot 1D transient field
    delete(lez);
    delete(lhy);
    lez = line(Ez_positions, Ez, 'Color', 'b', 'LineWidth', 1.5);
    lhy = line(Hy_positions, 377*Hy, 'Color', 'r', 'LineWidth', 1.5, 'LineStyle', '-.');
    ts  = num2str(time_step);
    ti  = num2str(dt*time_step*1e9);
    title(['time step = ' ts ' , time = ' ti  ' ns']);
    drawnow;
 
end

仿真结果如图：
【计算电磁学】真空中电磁波传播的MATLAB一维FDTD仿真程序

蓝色为电场，红色为磁场，右轴刻度没加。

参考
[1] Atef Elsherbeni, Veysel Demir, The Finite-Difference Time-Domain Method for ELectromagnetics with MATLAB Simulations [M]

【计算电磁学】真空中电磁波传播的MATLAB一维FDTD仿真程序

相关推荐