24小时热门版块排行榜    

查看: 4971  |  回复: 11

zuzero

铁虫 (小有名气)

[交流] [交流]分享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;



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

» 收录本帖的淘帖专辑推荐

薄膜制备 表征 科研学习 matlab

» 本帖已获得的红花(最新10朵)

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

欢迎来我的个人主页,方便学术交流与合作:jijinzu.buaa.edu.cn
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

alence

木虫 (著名写手)


小木虫(金币+0.5):给个红包,谢谢回帖
nice
要是三维的就更好了
知足常乐
2楼2011-08-22 19:56:16
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

天涯的兔子

木虫 (小有名气)

谢楼主分享
走狗的路,让猫们说去吧!!!
3楼2011-08-23 08:35:30
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
4楼2011-11-25 10:50:16
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

cbyhoo

新虫 (初入文坛)


小木虫: 金币+0.5, 给个红包,谢谢回帖
送鲜花一朵
非常感谢,呵呵,多向楼主请教啊
5楼2012-05-12 13:56:22
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

kevin123581

金虫 (正式写手)


小木虫: 金币+0.5, 给个红包,谢谢回帖
送鲜花一朵
感谢啊  我做的PML完全没有吸收效果。。。头疼
6楼2012-05-24 18:18:31
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

maidou520

禁虫 (初入文坛)


小木虫: 金币+0.5, 给个红包,谢谢回帖
送红花一朵
本帖内容被屏蔽

7楼2013-12-17 16:36:36
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xttzhb

金虫 (小有名气)

马克

[ 发自小木虫客户端 ]
8楼2013-12-22 23:39:03
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hfhhfh44

银虫 (初入文坛)


小木虫: 金币+0.5, 给个红包,谢谢回帖
LZ   我就试了试咋一点反应都没呢
9楼2013-12-29 16:50:57
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

helmhlotz

金虫 (著名写手)


小木虫: 金币+0.5, 给个红包,谢谢回帖
引用回帖:
7楼: Originally posted by maidou520 at 2013-12-17 16:36:36
弱弱的问一句,程序中的笑脸是什么什么意思?

: )如果中间没有空格的话,刚好和网页的表情代码一样
10楼2014-02-20 09:39:31
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 zuzero 的主题更新
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 一志愿武汉理工材料工程专硕调剂 +5 Doleres 2026-03-19 5/250 2026-03-19 20:14 by 制度的
[考研] 一志愿南京理工大学085701资源与环境302分求调剂 +3 葵梓卫队 2026-03-18 5/250 2026-03-19 19:35 by 给你你注意休息
[考研] 一志愿中国海洋大学,生物学,301分,求调剂 +5 1孙悟空 2026-03-17 5/250 2026-03-19 18:03 by zcl123
[考研] 一志愿西安交通大学材料工程专业 282分求调剂 +5 枫桥ZL 2026-03-18 7/350 2026-03-19 14:52 by 功夫疯狂
[考研] 346求调剂[0856] +3 WayneLim327 2026-03-16 6/300 2026-03-19 11:21 by WayneLim327
[考研] 一志愿985,本科211,0817化学工程与技术319求调剂 +10 Liwangman 2026-03-15 10/500 2026-03-19 10:25 by 无际的草原
[考研] 本科郑州大学物理学院,一志愿华科070200学硕,346求调剂 +4 我不是一根葱 2026-03-18 4/200 2026-03-19 09:11 by 浮云166
[考研] 0703化学调剂,求各位老师收留 +10 秋有木北 2026-03-14 10/500 2026-03-19 05:52 by anny19840123
[考研] 一志愿华中科技大学,080502,354分求调剂 +4 守候夕阳CF 2026-03-18 4/200 2026-03-18 22:16 by li123456789.
[考研] 一志愿武理材料305分求调剂 +5 想上岸的鲤鱼 2026-03-18 6/300 2026-03-18 17:53 by 无际的草原
[考研] 311求调剂 +6 26研0 2026-03-15 6/300 2026-03-18 14:43 by haxia
[考研] 考研求调剂 +3 橘颂. 2026-03-17 4/200 2026-03-17 21:43 by 有只狸奴
[考研] 296求调剂 +5 大口吃饭 身体健 2026-03-13 5/250 2026-03-17 21:05 by 不惑可乐
[考研] 26考研求调剂 +6 丶宏Sir 2026-03-13 6/300 2026-03-17 16:13 by 醉在风里
[考研] 274求调剂 +5 时间点 2026-03-13 5/250 2026-03-17 07:34 by 热情沙漠
[考研] 机械专硕325,寻找调剂院校 +3 y9999 2026-03-15 5/250 2026-03-16 19:58 by y9999
[考研] 080500,材料学硕302分求调剂学校 +4 初识可乐 2026-03-14 5/250 2026-03-14 21:08 by peike
[考研] 中科大材料与化工319求调剂 +3 孟鑫材料 2026-03-14 3/150 2026-03-14 20:10 by ms629
[考研] 297一志愿上交085600求调剂 +5 指尖八千里 2026-03-14 5/250 2026-03-14 17:26 by a不易
[考研] 330求调剂 +3 ?酱给调剂跪了 2026-03-13 3/150 2026-03-14 10:13 by JourneyLucky
信息提示
请填处理意见