FDTD中的数值稳定性

FDTD中的数值稳定性

在FDTD的算法中，电场与磁场在空间和时间内的取样都是在离散点进行的。

% 稳定性条件，dt<= dx/(n^0.5 * c),n维度，c光速(c0=299792458m/s)
% 空间步dx=1mm,n=1,dt<=3.3356*10^(-12)s，测试稳定性dt1=3.3356ps,dt2=3.3357ps
dt1=3.3356e-12;
dt2=3.3357e-12;
% 1、设置固定值信息
eps_0 = 8.854187817e-12;
mu_0 = 4*pi*1e-7;
c = 1/sqrt(mu_0*eps_0); % 光速
% 2、定义问题空间的参数
domain_size = 1; %空间维度1
dx = 1e-3; % 空间步大小dx=1mm，胞元的大小，米为单位
dt = 3e-12; % 时间步长，以秒为单位，dt = 3*10^(-12)s
number_of_time_steps = 2000; % 时间长度2000
nx = round(domain_size/dx); % 一维问题空间胞元个数1000，% round函数：四舍五入取整函数
source_position = 0.5; % 电流屏，源放置的位置


[Ceze1,Cezhy1,Cezj1,Ez1,Jz1,Chyh1,Chyez1,Chym1,Hy1,My1,time1]=field_initial(dt1,number_of_time_steps);
a = 2e-10; % 2*10^(-10)
b = 5e-11; % 5*10^(-11)
Jz_waveform = exp(-((time1-a)/b).^2); % 高斯源
source_position_indes = round(nx*source_position/domain_size)+1; 
% FDTD循环
for time_step = 1:number_of_time_steps
    % 更新当前时间下的JZ值
    Jz1(source_position_indes)=Jz_waveform(time_step);
    
    % 更新磁场
    Hy1(1:nx) = Chyh1(1:nx).*Hy1(1:nx) + Chyez1(1:nx).*(Ez1(2:nx+1)-Ez1(1:nx)) + Chym1(1:nx).*My1(1:nx);
    % 注意：Ez(2:nx+1)-Ez(1:nx)一一对应，每一次都计算了整个空间的值
    
    % 更新电场
    Ez1(2:nx) = Ceze1(2:nx).*Ez1(2:nx) + Cezhy1(2:nx).*(Hy1(2:nx)-Hy1(1:nx-1)) + Cezj1(2:nx).*Jz1(2:nx);
    Ma_Ez1(time_step) = max(abs(Ez1(2:nx)));
    t1(time_step)=time_step*dt1*10^9;
    % 绘制场量动态图
    % plot_fields;  
end

[Ceze2,Cezhy2,Cezj2,Ez2,Jz2,Chyh2,Chyez2,Chym2,Hy2,My2,time2]=field_initial(dt2,number_of_time_steps);
a = 2e-10; % 2*10^(-10)
b = 5e-11; % 5*10^(-11)
Jz_waveform = exp(-((time2-a)/b).^2); % 高斯源
source_position_indes = round(nx*source_position/domain_size)+1; 
% FDTD循环
for time_step = 1:number_of_time_steps
    % 更新当前时间下的JZ值
    Jz2(source_position_indes)=Jz_waveform(time_step);
    
    % 更新磁场
    Hy2(1:nx) = Chyh2(1:nx).*Hy2(1:nx) + Chyez2(1:nx).*(Ez2(2:nx+1)-Ez2(1:nx)) + Chym2(1:nx).*My2(1:nx);
    % 注意：Ez(2:nx+1)-Ez(1:nx)一一对应，每一次都计算了整个空间的值
    
    % 更新电场
    Ez2(2:nx) = Ceze2(2:nx).*Ez2(2:nx) + Cezhy2(2:nx).*(Hy2(2:nx)-Hy2(1:nx-1)) + Cezj2(2:nx).*Jz2(2:nx);
    Ma_Ez2(time_step) = max(abs(Ez2(2:nx)));
    t2(time_step)=time_step*dt2*10^9;
    
    
    % 绘制场量动态图
    % plot_fields;  
end
hold on;
xlim([0,6]);%对x轴设定显示范围 
ylim([0,3]);%对Y轴设定显示范围 
grid on; % 开启网格
xlabel('时间/ns');
ylabel('电场E_{zmax}/（V/m）');
plot(t1,Ma_Ez1,'-','Linewidth',2);
plot(t2,Ma_Ez2,'-.','Linewidth',2);
legend('dt=3.3356ps','dt=3.3357ps','Location','northwest')
set(gca,'FontSize',14);
hold off;

其中用到的函数

function  [Ceze,Cezhy,Cezj,Ez,Jz,Chyh,Chyez,Chym,Hy,My,time] = field_initial(dt,number_of_time_steps)
% 1、设置固定值信息
eps_0 = 8.854187817e-12;
mu_0 = 4*pi*1e-7;

% 2、定义问题空间的参数
domain_size = 1; %空间维度1
dx = 1e-3; % 空间步大小dx=1mm，胞元的大小，米为单位
nx = round(domain_size/dx); % 一维问题空间胞元个数1000，% round函数：四舍五入取整函数

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); %eps_r_z除数不能为0
sigma_e_z = zeros(nx+1,1);

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); %mu_r_y除数不能为0
sigma_m_y = zeros(nx,1);

% 求Ez n+1的更新公式 Ez = Ceze*Ez + Cezhy*(Hy(i)-Hy(i-1)) + Cezj*Jz;
Ceze = (2*eps_r_z*eps_0 - dt*sigma_e_z)./(2*eps_r_z*eps_0 + dt*sigma_e_z);
Cezhy = (2*dt)./((2*eps_r_z*eps_0 + dt*sigma_e_z)*dx);
Cezj = (-2*dt)./(2*eps_r_z*eps_0 + dt*sigma_e_z);

% 求Hy n+1的更新公式 Hy = Chyh*Hy + Chyez*(Ez(i+1)-Ez(i)) + Chym*My;
Chyh = (2*mu_r_y*mu_0 - dt*sigma_m_y)./(2*mu_r_y*mu_0 + dt*sigma_m_y);
Chyez = 2*dt./((2*mu_r_y*mu_0 + dt*sigma_m_y)*dx);
Chym = 2*dt./(2*mu_r_y*mu_0 + dt*sigma_m_y);

% 电流屏Jz 高斯源
time = dt*(0:number_of_time_steps-1).'; % .'转置,计算每一个时间节点方便带入JZ
end

运行结果：