| 查看: 5512 | 回复: 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朵）
» 猜你喜欢
40+中年人失业后，还能找到好工作吗？
已经有12人回复
-
用NBS溴化时溶剂除了CCl4外还有什么
已经有6人回复
-
申请科研项目纯靠实力能成功吗？
已经有11人回复
-
极限拉扯，如何抉择
已经有15人回复
-
二审4个多月了还是Awaiting Reviewer Scores
已经有4人回复
-
OpenAI又攻克数学难题了，离AI做有机合成还有多远？
已经有13人回复
-
课题组会帮你通过考核么？
已经有7人回复
-
售一区SCI文章T0P，我:8O.551.O54,科目全,可十急
已经有9人回复
-
售SCI一区文章，我:8O5.5.1.O5.4,科目全,可伽急
已经有10人回复
-
售SCI一区T0P文章，我:8.O55.1.O.54,科目全,可十急
已经有8人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 怎样把这段代码变成matlab代码
已经有5人回复
- PML的參數設定
已经有4人回复
- 有限元第二类边界条件 matlab仿真怎么编
已经有3人回复
- 迅雷5.8.14.706精简秒开版【转载】
已经有78人回复
- matlab数值求解边界条件微分方程组
已经有7人回复
- matlaB图像边界问题，简化程序
已经有3人回复
- 【求助】关于MATLAB中元胞自动机的边界条件和均匀分布
已经有4人回复
★
小木虫: 金币+0.5, 给个红包，谢谢回帖
送红花一朵
小木虫: 金币+0.5, 给个红包，谢谢回帖
送红花一朵 |
本帖内容被屏蔽 |
7楼2013-12-17 16:36:36
alence
木虫 (著名写手)
- 应助: 14 (小学生)
- 金币: 2929.1
- 散金: 70
- 红花: 10
- 帖子: 1153
- 在线: 652.8小时
- 虫号: 347044
- 注册: 2007-04-15
- 性别: GG
- 专业: 电磁场与波
2楼2011-08-22 19:56: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
200