24小时热门版块排行榜

返回列表

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

zuzero

铁虫 (小有名气)

应助: 0 (幼儿园)
金币: 23.5
散金: 40
红花: 2
帖子: 72
在线: 16.4小时
虫号: 1205011
注册: 2011-02-17
性别: GG
专业: 电磁场与波

[交流] [交流]分享Matlab写的PML吸收边界条件。代码很精简已有11人参与

%% 2-D FDTD PML吸收边界条件，点源辐射验证
% TE
clc;
clear;
fs = 20;
ft = 1.5;
% 边界的下标,垂直于x轴的
% 网格
Nsign = [1 -1];
diagSign = diag(Nsign);
Nx = 250;
Ny = 250;
% 吸收层厚度
Na = 20;
bxREx = 1 : Nx;
bxiREy = 2 : Nx;
byREy = 1 : Ny;
byiREx = 2 : Ny;
bxRH = 1 : Nx;
byRH = 1 : Ny;
% 场量
Hzx = zeros(Nx, Ny);
Hzy = zeros(Nx, Ny);
Ex = zeros(Nx, Ny + 1);
Ey = zeros(Nx + 1, Ny);
sigmamx = Hzx;
sigmamy = Hzy;
sigmax = Ey;
sigmay = Ex;
% 设定参数
% 最大电导率
sigmaMax = 0.2;
% 线性变化
sigma = sigmaMax * (Na : -1 : 1) / Na;
sigmam = sigmaMax * (2 * Na - 1 : -2 : 1) / 2 / Na;
sigmax([1 : Na, end : -1 : end - Na + 1],

= repmat([sigma, sigma]', 1, Ny);
sigmamx([1 : Na, end : -1 : end - Na + 1],

= repmat([sigmam, sigmam]', 1, Ny);
sigmay(:, [1 : Na, end : -1 : end - Na + 1]) = repmat([sigma, sigma], Nx, 1);
sigmamy(:, [1 : Na, end : -1 : end - Na + 1]) = repmat([sigmam, sigmam], Nx, 1);
% 计算系数
f1Ex = exp(-sigmay(bxREx, byiREx));
f1Ey = exp(-sigmax(bxiREy, byREy));
f2Ex = zeros(size(Ex));
f2Ey = zeros(size(Ey));
f2Ex(sigmay == 0) = 1 / ft;
f2Ey(sigmax == 0) = 1 / ft;
f2Ex(sigmay ~= 0) = (1 - exp(-sigmay(sigmay ~= 0))) ./ sigmay(sigmay ~= 0) / ft;
f2Ey(sigmax ~= 0) = (1 - exp(-sigmax(sigmax ~= 0))) ./ sigmax(sigmax ~= 0) / ft;
f2Ex = f2Ex(bxREx, byiREx);
f2Ey = f2Ey(bxiREy, byREy);

f1Hzx = exp(-sigmamx);
f1Hzy = exp(-sigmamy);
f2Hzx = zeros(size(Hzx));
f2Hzy = zeros(size(Hzy));
f2Hzx(sigmamx == 0) = 1 / ft;
f2Hzy(sigmamy == 0) = 1 / ft;
f2Hzx(sigmamx ~= 0) = (1 - exp(-sigmamx(sigmamx ~= 0))) ./ sigmamx(sigmamx ~= 0) / ft;
f2Hzy(sigmamy ~= 0) = (1 - exp(-sigmamy(sigmamy ~= 0))) ./ sigmamy(sigmamy ~= 0) / ft;
% 中心位置
centerX = floor((1 + Nx) / 2);
centerY = floor((1 + Ny) / 2);
% 迭代
for nn = 1 : 500
%    激励
Hzx(centerX, centerY) = Hzx(centerX, centerY) + sin(2 * pi * nn / ft / fs);
Hzy(centerX, centerY) = Hzy(centerX, centerY) + sin(2 * pi * nn / ft / fs);
%    更新E
Ex(bxREx, byiREx) = f1Ex .* Ex(bxREx, byiREx) + f2Ex .* (Hzx(bxREx, byiREx) - Hzx(bxREx, byiREx - 1) + Hzy(bxREx, byiREx) - Hzy(bxREx, byiREx - 1));
Ey(bxiREy, byREy) = f1Ey .* Ey(bxiREy, byREy) - f2Ey .* (Hzx(bxiREy, byREy) - Hzx(bxiREy - 1, byREy) + Hzy(bxiREy, byREy) - Hzy(bxiREy - 1, byREy));
%    更新H
Hzx = f1Hzx .* Hzx - f2Hzx .* (Ey(bxRH + 1,

- Ey(bxRH,

);
Hzy = f1Hzy .* Hzy + f2Hzy .* (Ex(:, byRH + 1) - Ex(:, byRH));
end
% 显示
figure('color', 'white');
imagesc(Hzx + Hzy);
colorbar;
axis image;

希望能对大家有所帮助。
下面是计算结果

回复此楼