【Mur吸收边界区域的2D FDTD】具有Mur吸收边界条件（ABC）的区域的二维时域有限差分法（2D FDTD）（Matlab代码实现）

本文链接：https://blog.csdn.net/qq_46009436/article/details/143432569

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

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

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

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

目录

⛳️赠与读者

💥1 概述

一、引言

二、理论基础

三、模拟设置

四、模拟过程与结果分析

五、优化与改进

六、结论

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

⛳️赠与读者

👨‍💻做科研，涉及到一个深在的思想系统，需要科研者逻辑缜密，踏实认真，但是不能只是努力，很多时候借力比努力更重要，然后还要有仰望星空的创新点和启发点。当哲学课上老师问你什么是科学，什么是电的时候，不要觉得这些问题搞笑。哲学是科学之母，哲学就是追究终极问题，寻找那些不言自明只有小孩子会问的但是你却回答不出来的问题。建议读者按目录次序逐一浏览，免得骤然跌入幽暗的迷宫找不到来时的路，它不足为你揭示全部问题的答案，但若能让人胸中升起一朵朵疑云，也未尝不会酿成晚霞斑斓的别一番景致，万一它居然给你带来了一场精神世界的苦雨，那就借机洗刷一下原来存放在那儿的“躺平”上的尘埃吧。

或许，雨过云收，神驰的天地更清朗.......🔎🔎🔎

💥1 概述

具有Mur吸收边界条件（ABC）的二维时域有限差分法（2D FDTD）研究

包含xy平面极化磁场（具有Hy和Hx分量）和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循环内，跨越整个时间网格。此外，在这里，用作更新方程的乘法因子的矩阵在循环开始之前被初始化，以避免在每个循环迭代中重复计算相同的矩阵，这是优化过程中的一个小尝试。这里的边界条件是Mur的吸收边界条件（ABC），其中网格点处的场具有使用Engquist Majda单向波动方程公式化的电场值，其中边界给人一种吸收入射到它们上的总场并且不将任何场反射回域的感觉。

A 2D TM wave containing the xy-plane polarized magnetic field having components Hy and Hx and z-polarized electric field Ez. The fields are updated at every time step, in a space, where all physical parameters of free space are not normalized to 1 but are given real and known values. The update is done using standard update equations obtained from the difference form of Maxwell's curl equations with very low electric and magnetic conductivities of 4x10^(-4) units incorporated in them. The field points are defined in a grid described by Yee's algorithm. The H fields are defined at every half coordinate of spacesteps. More precisely, the Hx part is defined at every half y-coordinate and full x-coordinate and the Hy part is defined at every half x-coordinate and full y-coordinate and E fields i.e the Ez part is defined at every full x and full y-coordinate points.Also here, the space-step length is taken as 1 micron instead of 1 unit in unitless domain assumed in previous programs. Also, the time update is done using Leapfrog time-stepping. Here, H-fields i.e Hx and Hy are updated every half time-step and E fields i.e Ez are updated every full time-step. This is shown by two alternating vector updates spanning only a part of spatial grid where the wave, starting from source, has reached at that particular time instant avoiding field updates at all points in the grid which is unnecessary at that time instant. These spatial updates are inside the main for-loop for time update, spanning the entire time grid. Also, here, the matrices used as multiplication factors for update equations are initialized before the loop starts to avoid repeated calculation of the same in every loop iteration, a minor attempt at optimization. The boundary condition here is Mur's Absorbing Boundary Condition (ABC) where the fields at the grid points have electric field values formulated using Engquist Majda one way wave equations where the boundaries give a sense of absorbing the total field incident on them and reflecting none back to the domain.

一、引言

时域有限差分法（FDTD）是一种直接求解麦克斯韦方程组的数值方法，在电磁波传播、散射和辐射等领域具有广泛应用。在二维（2D）FDTD模拟中，通常考虑TM波或TE波，其中电场和磁场分量在特定的平面内分布。本文重点研究具有Mur吸收边界条件（ABC）的2D FDTD方法，该方法用于模拟电磁波在二维空间中的传播，并减少边界处的反射。

二、理论基础

麦克斯韦方程组：
- 麦克斯韦方程组描述了电磁波在介质中的传播规律，包括电场和磁场的相互作用。
- 在二维FDTD模拟中，通常考虑TM波，其中电场分量Ez和磁场分量Hx、Hy在xy平面内分布。
Yee算法：
- Yee算法是FDTD方法的核心，它定义了电场和磁场分量在空间和时间上的分布和更新规则。
- 在Yee算法中，电场分量位于网格的节点上，而磁场分量位于网格的边缘上。
Mur吸收边界条件：
- Mur吸收边界条件是一种有效的边界处理方法，用于减少FDTD模拟中边界处的反射。
- 它基于Engquist-Majda单向波动方程，通过计算边界处的电场值来模拟电磁波的吸收。

三、模拟设置

网格划分：
- 采用二维FDTD网格，网格大小根据模拟需求确定。
- 空间步长取为1微米，时间步长根据蛙跳时间步进方法确定。
场分量定义：
- 电场分量Ez位于网格的全x和全y坐标点上。
- 磁场分量Hx位于每个半y坐标和全x坐标处，Hy位于每个半x坐标和全y坐标处。
物理参数：
- 电导率和磁导率取为非常低的值，分别为4x10^(-4)单位。
- 其他物理参数（如介电常数和磁导率张量）根据模拟需求设置。
更新方程：
- 使用从麦克斯韦旋度方程的差分形式中获得的标准更新方程来更新电场和磁场分量。
- 初始化用于更新方程的乘法因子矩阵，以避免在每个循环迭代中重复计算。
边界条件：
- 采用Mur吸收边界条件，通过计算边界处的电场值来模拟电磁波的吸收。
- 边界条件应给人一种吸收入射到它们上的总场并且不将任何场反射回域的感觉。

四、模拟过程与结果分析

初始化：
- 设置网格大小、物理参数、边界条件等。
- 初始化电场和磁场分量。
时间步进：
- 使用蛙跳时间步进方法更新电场和磁场分量。
- 每半个时间步长更新一次磁场分量（Hx和Hy），每完整时间步长更新一次电场分量（Ez）。
边界处理：
- 在每个时间步长结束时，计算边界处的电场值，应用Mur吸收边界条件。
结果分析：
- 观察并分析电磁波在二维空间中的传播过程。
- 验证Mur吸收边界条件在减少边界反射方面的有效性。
- 分析电磁波的传播速度、衰减等特性。

五、优化与改进

矩阵初始化：
- 在循环开始之前初始化用于更新方程的乘法因子矩阵，以减少计算量。
并行计算：
- 考虑使用并行计算技术来加速模拟过程，如多线程或GPU加速。
边界条件优化：
- 进一步研究和优化Mur吸收边界条件，以提高其吸收效果和计算效率。
物理参数调整：
- 根据模拟需求调整电导率、磁导率等物理参数，以更准确地模拟实际物理过程。

六、结论

本文研究了具有Mur吸收边界条件的二维时域有限差分法（2D FDTD）方法，并成功模拟了电磁波在二维空间中的传播过程。通过验证Mur吸收边界条件在减少边界反射方面的有效性，本文为电磁波传播模拟提供了一种有效的边界处理方法。未来的研究可以进一步优化和改进该方法，以提高其计算效率和准确性。

📚2 运行结果

部分代码：

%Vector update instead of for-loop for Hy and Hx fields
Hx(n1:n2-1,n11:n22-1)=A(n1:n2-1,n11:n22-1).*Hx(n1:n2-1,n11:n22-1)-B(n1:n2-1,n11:n22-1).*(Ez(n1:n2-1,n11+1:n22)-Ez(n11:n2-1,n11:n22-1));
Hy(n1:n2-1,n11:n22-1)=A(n1:n2-1,n11:n22-1).*Hy(n1:n2-1,n11:n22-1)+B(n1:n2-1,n11:n22-1).*(Ez(n1+1:n2,n11:n22-1)-Ez(n11:n2-1,n11:n22-1));

%Vector update instead of for-loop for Ez field
Ez(n1+1:n2-1,n11+1:n22-1)=C(n1+1:n2-1,n11+1:n22-1).*Ez(n1+1:n2-1,n11+1:n22-1)+(Hy(n1+1:n2-1,n11+1:n22-1)-Hy(n1:n2-2,n11+1:n22-1)-Hx(n1+1:n2-1,n11+1:n22-1)+Hx(n1+1:n2-1,n11:n22-2)).*D(n1+1:n2-1,n11+1:n22-1);

%Mur's abc conditions obtained from Mur's difference equation for
%forward boundary
if n>=xdim-2-xsource
Ez(xdim-2,3:1:ydim-3)=c0efffor*(Ez(xdim-3,3:1:ydim-3)+prev_prev_xfor(1,3:1:ydim-3))-prev_prev_x_minus_1for(1,3:1:ydim-3)+c2efffor*(prev_xfor(1,3:1:ydim-3)+prev_x_minus_1for(1,3:1:ydim-3))+c3efffor*(prev_x_minus_1for(1,2:1:ydim-4)+prev_x_minus_1for(1,4:1:ydim-2)+prev_xfor(1,2:1:ydim-4)+prev_xfor(1,4:1:ydim-2));
end

%Storage vectors for boundary and boundary-1 values of previous and its
%previous time steps updated at forward boundary
prev_prev_xfor=prev_xfor;
prev_prev_x_minus_1for=prev_x_minus_1for;
prev_xfor(1,1:1:ydim)=Ez(xdim-2,1:1:ydim);
prev_x_minus_1for(1,1:1:ydim)=Ez(xdim-3,1:1:ydim);

%Mur's abc conditions obtained from Mur's difference equation for
%backward boundary
if n>=xsource-3
Ez(2,3:1:ydim-3)=-prev_prev_xrev(1,3:1:ydim-3)+c1effrev*(Ez(3,3:1:ydim-3)+prev_prev_x_minus_1rev(1,3:1:ydim-3))+c2effrev*(prev_xrev(1,3:1:ydim-3)+prev_x_minus_1rev(1,3:1:ydim-3))+c3effrev*(prev_x_minus_1rev(1,2:1:ydim-4)+prev_x_minus_1rev(1,4:1:ydim-2)+prev_xrev(1,2:1:ydim-4)+prev_xrev(1,4:1:ydim-2));
end

%Storage vectors for boundary and boundary-1 values of previous and its
%previous time steps updated at backward boundary
prev_prev_xrev=prev_xrev;
prev_prev_x_minus_1rev=prev_x_minus_1rev;
prev_xrev(1,1:1:ydim)=Ez(3,1:1:ydim);
prev_x_minus_1rev(1,1:1:ydim)=Ez(2,1:1:ydim);

%Mur's abc conditions obtained from Mur's difference equation for
%upward boundary
if n>=ydim-2-ysource
Ez(3:1:xdim-3,ydim-2)=c0efffor*(Ez(3:1:xdim-3,ydim-3)+prev_prev_yfor(3:1:xdim-3,1))-prev_prev_y_minus_1for(3:1:xdim-3,1)+c2efffor*(prev_yfor(3:1:xdim-3,1)+prev_y_minus_1for(3:1:xdim-3,1))+c3efffor*(prev_y_minus_1for(2:1:xdim-4,1)+prev_y_minus_1for(4:1:xdim-2,1)+prev_yfor(2:1:xdim-4,1)+prev_yfor(4:1:xdim-2,1));
end

%Storage vectors for boundary and boundary-1 values of previous and its
%previous time steps updated at upward boundary
prev_prev_yfor=prev_yfor;
prev_prev_y_minus_1for=prev_y_minus_1for;
prev_yfor(1:1:xdim,1)=Ez(1:1:xdim,ydim-2);
prev_y_minus_1for(1:1:xdim,1)=Ez(1:1:xdim,ydim-3);

%Mur's abc conditions obtained from Mur's difference equation for
%downward boundary
if n>=ysource-3