| 查看: 4815 | 回复: 11 | ||||||
[交流]
[交流]分享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朵）
» 猜你喜欢
Cas 72-43-5需要30g，定制合成，能接单的留言
已经有8人回复
-
求助:我三月中下旬出站，青基依托单位怎么办？
已经有6人回复
-
北京211副教授，35岁，想重新出发，去国外做博后，怎么样？
已经有8人回复
-
磺酰氟产物，毕不了业了！
已经有5人回复
-
论文终于录用啦！满足毕业条件了
已经有25人回复
-
2026年机械制造与材料应用国际会议 （ICMMMA 2026）
已经有3人回复
-
自荐读博
已经有3人回复
-
不自信的我
已经有5人回复
-
投稿Elsevier的杂志（返修），总是在选择OA和subscription界面被踢皮球
已经有8人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 怎样把这段代码变成matlab代码
已经有5人回复
- PML的參數設定
已经有4人回复
- 有限元第二类边界条件 matlab仿真怎么编
已经有3人回复
- 迅雷5.8.14.706精简秒开版【转载】
已经有78人回复
- matlab数值求解边界条件微分方程组
已经有7人回复
- matlaB图像边界问题，简化程序
已经有3人回复
- 【求助】关于MATLAB中元胞自动机的边界条件和均匀分布
已经有4人回复
alence
木虫 (著名写手)
- 应助: 14 (小学生)
- 金币: 2929.1
- 散金: 70
- 红花: 10
- 帖子: 1153
- 在线: 651.9小时
- 虫号: 347044
- 注册: 2007-04-15
- 性别: GG
- 专业: 电磁场与波
2楼2011-08-22 19:56:16
天涯的兔子
木虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 4900.2
- 散金: 201
- 红花: 1
- 帖子: 266
- 在线: 71.8小时
- 虫号: 981001
- 注册: 2010-03-24
- 性别: GG
- 专业: 电磁场与波
3楼2011-08-23 08:35:30
 |
4楼2011-11-25 10:50:16
5楼2012-05-12 13:56:22
kevin123581
金虫 (正式写手)
- 应助: 12 (小学生)
- 金币: 761.3
- 散金: 81
- 红花: 6
- 帖子: 539
- 在线: 238.7小时
- 虫号: 1658947
- 注册: 2012-03-02
- 性别: MM
- 专业: 光学
6楼2012-05-24 18:18:31
★
小木虫: 金币+0.5, 给个红包，谢谢回帖
送红花一朵
小木虫: 金币+0.5, 给个红包，谢谢回帖
送红花一朵 |
本帖内容被屏蔽 |
7楼2013-12-17 16:36:36
xttzhb
金虫 (小有名气)
- 应助: 11 (小学生)
- 金币: 242.6
- 散金: 160
- 帖子: 230
- 在线: 256.2小时
- 虫号: 1274962
- 注册: 2011-04-23
- 性别: GG
- 专业: 电磁场与波
8楼2013-12-22 23:39:03
9楼2013-12-29 16:50:57
helmhlotz
金虫 (著名写手)
- 应助: 0 (幼儿园)
- 金币: 903.2
- 帖子: 1006
- 在线: 32.6小时
- 虫号: 1021129
- 注册: 2010-05-18
- 性别: GG
- 专业: 橡胶及弹性体
10楼2014-02-20 09:39:31