具有完美匹配层 （PML） 边界条件的区域的二维 FDTD研究（Matlab代码实现）

 👨‍🎓个人主页：研学社的博客   

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

---

💥1 概述

包含具有Hy和Hx分量的xy平面极化磁场以及z极化电场Ez的2D TM波。字段在每个时间步长更新，在一个空间中，其中自由空间的所有物理参数未规范化为 1，而是给定实值和已知值。更新是使用从麦克斯韦卷曲方程的差分形式获得的标准更新方程完成的，其中包含非常低的 4x10^（-4） 单位的电导率和磁导率。场点在Yee算法描述的网格中定义。H 字段在空间步的每个半坐标处定义。更准确地说，Hx 部分在每个半 y 坐标和完整 x 坐标处定义，Hy 部分在每半 x 坐标和完整 y 坐标和 E 场处定义，即 Ez 部分在每个完整的 x 和完整的 y 坐标点上定义。同样在这里，空间步长被取为1微米，而不是以前程序中假设的无单位域中的1个单位。此外，时间更新是使用跨越时间步长完成的。在这里，H 字段（即 Hx 和 Hy）每半个时间步更新一次，E 字段（即 Ez）每半个时间步更新一次。这表现为两个交替的矢量更新，仅跨越空间网格的一部分，其中波从源开始，已到达该特定时间时刻，避免了网格中所有点的场更新，这在当时是不必要的。这些空间更新位于时间更新的主 for 循环内，跨越整个时间网格。此外，在这里，用作更新方程乘法因子的矩阵在循环开始之前初始化，以避免在每次循环迭代中重复计算相同的矩阵，这是优化的小尝试。这里的边界条件是完全匹配层 （PML） 边界条件，其中边界附近的场在到达边界宽度的预定长度上衰减到边界处的零值，使用多项式递增的电导率值在边界宽度上，在边界宽度处最大值，并在边界宽度的每个点选择一个磁导率值以避免在该点反射。此外，这里使用Berenger的PML条件，其中在场中Ez被分成两个分量Ezx和Ezy，并且分量在两个方向上使用单独的电导率和磁导率衰减（x方向的sigmax和sigma_starx以及y方向的sigmay和sigma_stary）。

📚2 运行结果

 

 部分代码：

%Total no of time steps
time_tot=350;

%Position of the source (center of the domain)
xsource=100;
ysource=100;

%Courant stability factor
S=1/(2^0.5);

% Parameters of free space (permittivity and permeability and speed of
% light) are all not 1 and are given real values
epsilon0=(1/(36*pi))*1e-9;
mu0=4*pi*1e-7;
c=3e+8;

% Spatial grid step length (spatial grid step= 1 micron and can be changed)
delta=1e-6;
% Temporal grid step obtained using Courant condition
deltat=S*delta/c;

% Initialization of field matrices
Ez=zeros(xdim,ydim);
Ezx=zeros(xdim,ydim);
Ezy=zeros(xdim,ydim);
Hy=zeros(xdim,ydim);
Hx=zeros(xdim,ydim);

% Initialization of permittivity and permeability matrices
epsilon=epsilon0*ones(xdim,ydim);
mu=mu0*ones(xdim,ydim);

% Initializing electric conductivity matrices in x and y directions
sigmax=zeros(xdim,ydim);
sigmay=zeros(xdim,ydim);

%Perfectly matched layer boundary design
%Reference:-http://dougneubauer.com/wp-content/uploads/wdata/yee2dpml1/yee2d_c.txt
%(An adaptation of 2-D FDTD TE code of Dr. Susan Hagness)

%Boundary width of PML in all directions
bound_width=25;

%Order of polynomial on which sigma is modeled
gradingorder=6;

%Required reflection co-efficient
refl_coeff=1e-6;

%Polynomial model for sigma
sigmamax=(-log10(refl_coeff)*(gradingorder+1)*epsilon0*c)/(2*bound_width*delta);
boundfact1=((epsilon(xdim/2,bound_width)/epsilon0)*sigmamax)/((bound_width^gradingorder)*(gradingorder+1));
boundfact2=((epsilon(xdim/2,ydim-bound_width)/epsilon0)*sigmamax)/((bound_width^gradingorder)*(gradingorder+1));
boundfact3=((epsilon(bound_width,ydim/2)/epsilon0)*sigmamax)/((bound_width^gradingorder)*(gradingorder+1));
boundfact4=((epsilon(xdim-bound_width,ydim/2)/epsilon0)*sigmamax)/((bound_width^gradingorder)*(gradingorder+1));
x=0:1:bound_width;
for i=1:1:xdim
    sigmax(i,bound_width+1:-1:1)=boundfact1*((x+0.5*ones(1,bound_width+1)).^(gradingorder+1)-(x-0.5*[0 ones(1,bound_width)]).^(gradingorder+1));
    sigmax(i,ydim-bound_width:1:ydim)=boundfact2*((x+0.5*ones(1,bound_width+1)).^(gradingorder+1)-(x-0.5*[0 ones(1,bound_width)]).^(gradingorder+1));
end
for i=1:1:ydim
    sigmay(bound_width+1:-1:1,i)=boundfact3*((x+0.5*ones(1,bound_width+1)).^(gradingorder+1)-(x-0.5*[0 ones(1,bound_width)]).^(gradingorder+1))';
    sigmay(xdim-bound_width:1:xdim,i)=boundfact4*((x+0.5*ones(1,bound_width+1)).^(gradingorder+1)-(x-0.5*[0 ones(1,bound_width)]).^(gradingorder+1))';
end

%Magnetic conductivity matrix obtained by Perfectly Matched Layer condition
%This is also split into x and y directions in Berenger's model
sigma_starx=(sigmax.*mu)./epsilon;
sigma_stary=(sigmay.*mu)./epsilon;

%Choice of nature of source
gaussian=0;
sine=0;
% The user can give a frequency of his choice for sinusoidal (if sine=1 above) waves in Hz 

🎉3 参考文献

部分理论来源于网络，如有侵权请联系删除。

[1]李志. 基于FDTD的高阶PML算法及其在集成电路仿真中的应用[D].天津工业大学,2021.DOI:10.27357/d.cnki.gtgyu.2021.000394.

[2]Computational Electromagnetics At IIT Madras (2023). 2D FDTD of a region with Perfectly Matched Layer boundary.

🌈4 Matlab代码实现