、一、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.运行结果