一维FDTD仿真

、一、FDTD基础知识

1、麦克斯韦方程组

E:电场强度	J:电流密度
D:电位移	M：磁流密度
H:磁场强度	：电荷密度
B:磁通量密度	：磁荷密度

 

将上述（5）-（10）方程带入（1）（2）方程得：

 2、旋度公式：

百度百科 旋度（https://baike.baidu.com/item/%E6%97%8B%E5%BA%A6/8106439）

利用旋度公式拆开（11）（12）式子的左边，再在三维空间每个矢量方向分解就可以得到六个标量方程:

 3、中心差分方程：（将导数用中心差分的式子代替）

 4、Yee胞元

Yee给出了麦克斯韦旋度方程的一组差分方程。这组方程在空间和时间上以离散的形式给出，使用的是中心差分格式。电场和磁场分量可以在时间和空间中以离散点的方式取样。FDTD将三维问题的几何结构分解为单元，以构成相应的网络。

电场分量放置在YEE单元各棱中间，平行于各棱；磁场分量放置在Yee单元各面的中心，平行于各面的法线。

 FDTF算法在离散的时间瞬间取样和计算场值，但是电场和磁场取样计算并不是在相同的时刻。时间步为∆t ​​​​

         电场的取样时刻0，∆t，2∆t，3∆t，...，n∆t

        磁场的取样时刻0.5∆t，1.5∆t，2.5∆t，...，(n+0.5)∆t

即电场取样在时间的整数步长时刻，而磁场取样时刻为半整数时间步时刻。

 （2.1）方程带入中心差分来近似

 

  

二、MATLAB仿真一维的FDTD

 一维FDTD仿真，计算由z向电流屏Jz所产生的电场和磁场
        1、电流屏位于问题空间中心，区域填充空气
        2、空间长度1m，空间步dx=1mm
        3、电流屏为高斯波形Jz(t)=exp-{[t-2*10^(-10)]/[5*10^(-11)]}^2
        4、计算Ez,Hy，并画出传播图

1、代码

%% 一维FDTD仿真，计算由z向电流屏Jz所产生的电场和磁场
% 1、电流屏位于问题空间中心，区域填充空气
% 2、空间长度1m，空间步dx=1mm
% 3、电流屏为高斯波形Jz(t)=exp-{[t-2*10^(-10)]/[5*10^(-11)]}^2
% 计算Ez,Hy


%% 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; % 电流屏，源放置的位置

%% 3、参数初始化
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
a = 2e-10; % 2*10^(-10)
b = 5e-11; % 5*10^(-11)
Jz_waveform = exp(-((time-a)/b).^2); % 高斯源
source_position_indes = round(nx*source_position/domain_size)+1; 
% 源的位置索引，一维可以不看domain_size,nx为空间胞元个数，source_position源的位置（中心）1000*0.5=500，索引第501
% plot(Jz);

%% 初始化图
% initialize_plotting_parameters;%初始化图标参数，需要自己写
hold on
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];
figure(1);
axis([0 1 -0.2 0.2 -0.2 0.2]);
lez = line(Ez_positions,Ez*0,Ez,'Color','b','LineWidth',1.5); % line 会向当前坐标区添加线条，而不删除其他图形对象或重置坐标区属性。
lhy = line(Hy_positions,377*Hy,Hy*0,'Color','r','LineWidth',1.5,'LineStyle','-.');
set(gca,'FontSize',12,'FontWeight','bold'); % bold加粗
axis square; % 将当前坐标系图形设置为方形。横轴及纵轴比例是1：1
xlabel('x_{m}');
ylabel('[A/m]');
zlabel('[V/m]');
grid on; % 开启网格
view(3);
p = patch('vertices',v,'faces',f,'facecolor','g','facealpha',0.2);
% patch()创建一个或多个填充多边形,使用 X 和 Y 的元素作为每个顶点的坐标，以创建一个或多个填充多边形,facealphal表示颜色的深浅
text(0,1,1.1,'PEC','horizontalalignment','center','fontweight','bold');
text(1,1,1.1,'PEC','horizontalalignment','center','fontweight','bold');
legend('E_{z}','H_{y}_\times_{377}','Location','NorthEast');

%% FDTD循环
for time_step = 1:number_of_time_steps
    % 更新当前时间下的JZ值
    Jz(source_position_indes)=Jz_waveform(time_step);
    
    % 更新磁场
    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);
    % 注意：Ez(2:nx+1)-Ez(1:nx)一一对应，每一次都计算了整个空间的值
    
    % 更新电场
    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);
    
    %PEC完全吸收
    Ez(1)=0;
    Ez(nx)=0;
    
    % 绘制场量动态图
    % plot_fields;
    
    delete(lez);% delete函数会删除附带的曲线标志
    delete(lhy);
 
    lez = line(Ez_positions,Ez*0,Ez,'Color','b','LineWidth',1.5); % line 会向当前坐标区添加线条，而不删除其他图形对象或重置坐标区属性。
    lhy = line(Hy_positions,377*Hy,Hy*0,'Color','r','LineWidth',1.5,'LineStyle','-.');
    
    legend('PEC','E_{z}','H_{y}_\times_{377}','Location','NorthEast');
    
    ts = num2str(time_step); % 数字转字符串
    ti = num2str(dt*time_step*1e9);
    title(('time_step_=_"ts",_time_=_"ti"_ns'));
    
    drawnow limitrate% 更新数据并处理回调,实现动图
end
drawnow

 注：

2.运行结果