1522阅读
1回复

请教一个用MATLAB编的一维FDTD程序（有程序）

UID ：45883

注册：2009-11-05
登录：2022-04-17
发帖：117
等级：仿真二级

0楼发表于: 2009-11-14 15:42:00

要仿真一个一维的FDTD传播问题，想根据这个改动一下，在距离高斯点源15cm的地方加入一个相对介电常数为9，厚度为5CM的介质，边界条件由PEC改为吸收边界，请教大虾给条明路....
% 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 = 0.6;                  % 1D problem space length in meters
dx = 1e-3;                        % cell size in meters
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.2;            % 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);
Jz        = 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);

Cezj = (-2 * dt) ...
     ./(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);

Chym = (-2 * dt) ...
     ./(2 * mu_r_y * mu_0 + dt * sigma_m_y);

% Define the Gaussian source waveform
time       = dt*[0:number_of_time_steps-1].';
Jz_waveform = exp(-((time-8e-11)/2e-11).^2);
source_position_index = round(nx*source_position/domain_size)+1;

Ez_positions = [0:nx]*dx;
Hy_positions = ([0:nx-1]+0.5)*dx;
v = [0 -0.1 -0.1; 0 -0.1 0.1; 0 0.1 0.1; 0 0.1 -0.1; ...
     1 -0.1 -0.1; 1 -0.1 0.1; 1 0.1 0.1; 1 0.1 -0.1];
f = [1 2 3 4; 5 6 7 8];
axis([0 1 -0.2 0.2 -0.2 0.2]);
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');
axis square;
legend('E_{z}', 'H_{y} \times 377','Location','NorthEast');
xlabel('x [m]');
ylabel('[A/m]');
zlabel('[V/m]');
grid on;
p = patch('vertices', v, 'faces', f, 'facecolor', 'g', 'facealpha',0.2);

% FDTD loop
for time_step = 1:number_of_time_steps

    if (0.35=<Ez_positions<=0.4)%%%%%这样改肯定不对...额
eps_r_z =9.*ones (nx+1,1);
else
eps_r_z=ones (nx+1,1);
end
    % Update Jz for the current time step
    Jz(source_position_index) = Jz_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)) ...
         + Chym(1:nx) .* My(1:nx);

    % Update electric field
    Ez(2:nx) = Ceze (2:nx) .* Ez(2:nx) ...
             + Cezhy(2:nx) .* (Hy(2:nx) - Hy(1:nx-1)) ...
             + Cezj(2:nx) .* Jz(2:nx);

    Ez(1)    ..

未注册仅能浏览部分内容，查看全部内容及附件请先登录或注册

离线martinloong

UID ：45883

注册：2009-11-05
登录：2022-04-17
发帖：117
等级：仿真二级

1楼发表于: 2009-11-14 15:46:00

第78-82是源文件没有的，要删除

发帖回复