%% 电偶极辐射仿真（有限差分法）
clear; clc; close all;
%% 参数设置
lambda = 0.5; % 波长(m)
c = 3e8; % 光速(m/s)
f = c / lambda; % 频率(Hz)
omega = 2*pi*f; % 角频率
p0 = 1e-9; % 电偶极矩振幅(C·m)
sigma = 0.1; % 脉冲宽度(s)
t0 = 0.5/sigma; % 脉冲中心时间(s)
dx = lambda/50; % 空间步长(m)
dy = dx; % 二维仿真时dy=dx
dz = dx; % 三维仿真时dz=dx
dt = dx/c/1.5; % 时间步长(s)（满足CFL条件）
x_range = [-1 1]; % 计算区域x范围(m)
y_range = [-1 1]; % 计算区域y范围(m)
z_range = ; % 二维仿真时z=0
epsilon0 = 8.854187817e-12; % 真空介电常数(F/m)
mu0 = 4*pi*1e-7; % 真空磁导率(H/m)
%% 初始化网格
[x,y,z] = meshgrid(...
x_range(1):dx:x_range(2), ...
y_range(1):dy:y_range(2), ...
z_range(1):dz:z_range(2));
N_x = size(x,1); % 网格点数x方向
N_y = size(y,2); % 网格点数y方向
N_z = size(z,3); % 网格点数z方向
%% 初始化场量
Ex = zeros(N_x,N_y,N_z); % x方向电场
Ey = zeros(N_x,N_y,N_z); % y方向电场
Hz = zeros(N_x,N_y,N_z); % z方向磁场
%% 边界条件设置（PML吸收边界）
pml_order = 4; % PML阶数
pml厚 = 10; % PML厚度(网格单元数)
% 计算PML参数（这里仅简化实现）
sigma_x = zeros(N_x,N_y,N_z);
sigma_y = zeros(N_x,N_y,N_z);
for i = 1:pml厚
sigma_x(i,:,:) = i/pml厚 * 1e-3;
sigma_x(N_x-i+1,:,:) = i/pml厚 * 1e-3;
end
%% 时间迭代参数
t_max = 2*t0; % 总仿真时间(s)
n_t = round(t_max/dt); % 总时间步数
E记录 = cell(n_t,1); % 存储电场数据
%% 仿真主循环
for n = 1:n_t
t = n*dt; % 当前时间
% 电偶极激励（高斯脉冲）
p = p0 * exp(-(t - t0)^2/(2*sigma^2));
% 电场更新（FDTD公式）
Hz(2:end-1,2:end-1,2:end-1) = ...
Hz(2:end-1,2:end-1,2:end-1) ...
+ (dt/(mu0*dy)) * (Ey(2:end-1,3:end,2:end-1) - Ey(2:end-1,2:end-1,2:end-1)) ...
- (dt/(mu0*dx)) * (Ex(3:end,2:end-1,2:end-1) - Ex(2:end-1,2:end-1,2:end-1));
% 电场更新（考虑电偶极源）
Ex(floor(N_x/2),floor(N_y/2),floor(N_z/2)) = ...
Ex(floor(N_x/2),floor(N_y/2),floor(N_z/2)) ...
+ p/(epsilon0*dx*dy*dz) * dt;
% 边界处理（PML）
Hz = Hz .* exp(-sigma_x*dt/(2*mu0));
% 记录电场数据
E记录{n} = abs(Ex(:,:,1)); % 记录二维电场
% 可视化更新
if mod(n,50) == 0
figure(1);
surf(x(:,:,1), y(:,:,1), E记录{n}, 'EdgeColor','none');
shading interp;
title(['t = ',num2str(t,'%.2f'),' s']);
axis([-1 1 -1 1 0 max(E记录{n}(:))]
);
drawnow;
end
end
%% 三维场分布可视化
figure(2);
p = isosurface(x,y,z,abs(Ex),max(abs(Ex(:)))*0.5);
isonormals(x,y,z,abs(Ex),p);
patch(p,'FaceColor','interpolated','EdgeColor','none');
isonormals(x,y,z,abs(Ex),p); % 优化表面法线
daspect();
view(45,30); % 设置视角
camlight('headlight'); % 添加光源
lighting phong;
xlabel('x (m)'); ylabel('y (m)'); zlabel('z (m)');
title('三维电场分布');
colormap jet;
%% 近场特性分析
figure(3);
near_field = Ex(floor(N_x/2)-10:floor(N_x/2)+10, ...
floor(N_y/2)-10:floor(N_y/2)+10, ...
floor(N_z/2));
surf(abs(near_field));
title('近场分布');
xlabel('x方向'); ylabel('y方向'); zlabel('电场强度');中显示E记录有误，没有正确使用，请告诉我应用什么代替。显示pml厚也显示错误，用什么代替
扫一扫
492
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
