【折射现象的二维时域有限差分法】玻璃-空气界面折射现象的二维时域有限差分法（Matlab代码实现）-CSDN博客

本文链接：https://blog.csdn.net/weixin_67304359/article/details/148055977

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

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

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

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

目录

⛳️赠与读者

💥1 概述

一、引言

二、FDTD方法基本原理

三、模拟设置

四、模拟过程与结果

五、讨论

六、结论

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

⛳️赠与读者

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

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

💥1 概述

【折射现象的二维时域有限差分法】玻璃-空气界面折射现象的二维时域有限差分法研究文档

电场源被定义在空间域的顶部边缘（在玻璃中），位于左半部分，它是一个硬正弦源，不会因外部场的干扰而改变其值，换句话说，该源是一个完美的电导体。自由空间波长和波从源传播的平面与x轴的角度被指定为适合折射现象（相对于临界角）。随着波传播并撞击两种介质的界面，在每个时间步长上显示Ez场在整个空间域上的彩色标度图，其中波折射（弯曲）到另一种介质（即空气）中，因为α角小于两个界面之间的临界角（约42度）。

A source of electric field is defined at the top edge of the spatial domain,(in glass) on the left half, which is a hard sinusoidal source, in that it does not change its value due to interference from external fields i.e in other words,the source is a perfect electric conductor. The free space wavelength and the angle of the plane of propagation of the wave from the source with the x-axis are specified to suit refraction phenomenon (wrt critical angle). The color scaled plot of Ez field over the entire spatial domain is shown at every time step as the wave propagates and hits the interface of the two media where it refracts (bends) into the other medium (viz. air) as angle alpha would be lesser than critical angle (~42 degrees) between the two interfaces.

一、引言

光波在不同介质界面上的折射现象是经典光学的重要内容，对光学器件的设计和应用具有深远影响。精确模拟光波的传播行为对于理解和预测光学系统的性能至关重要。有限差分时域法（FDTD）作为一种强大的数值计算方法，因其能够直接求解麦克斯韦方程组，在电磁波模拟领域得到广泛应用。本文将利用二维FDTD方法模拟光波在玻璃-空气界面发生的折射现象。

二、FDTD方法基本原理

FDTD方法的核心思想是将麦克斯韦方程组在时间和空间上进行离散化，通过迭代计算得到电磁场的时空分布。对于二维情况，在横磁（TM）模式下，麦克斯韦方程组可以简化为以下形式：

∂Hz/∂t=(1/μ)(∂Ey/∂x−∂Ex/∂y)
∂Ex/∂t=(1/ε)(∂Hz/∂y)
∂Ey/∂t=−(1/ε)(∂Hz/∂x)

其中，Hz为z方向的磁场强度，Ex和Ey分别为x和y方向的电场强度，ε和μ分别为介质的介电常数和磁导率。利用中心差分格式对上述方程进行离散化，可以得到时间步长Δt和空间步长Δx、Δy下的递推公式。

三、模拟设置

电场源：定义在空间域的顶部边缘（在玻璃中），位于左半部分，是一个硬正弦源，不会因外部场的干扰而改变其值，即该源是一个完美的电导体。
自由空间波长和波的传播角度：指定自由空间波长和波从源传播的平面与x轴的角度，以适应折射现象（相对于临界角）。
介质参数：玻璃的介电常数通常大于空气的介电常数，这将导致光波在界面处发生折射。
边界条件：采用完全匹配层（Perfectly Matched Layer, PML）吸收边界条件，以有效地吸收入射波，减少边界反射。

四、模拟过程与结果

初始化电场和磁场：在计算域中初始化电场和磁场。
迭代计算：根据递推公式迭代计算电场和磁场的更新。
结果记录：在每个时间步长上记录整个空间域上Ez场的分布，以观察折射现象的发生和发展。
结果分析：通过彩色标度图展示不同时间步长上Ez场的分布，可以观察到光波在玻璃-空气界面处发生折射，折射角符合斯涅耳定律。同时，通过改变入射角、波长等参数，可以研究不同条件下光波的传播特性。

五、讨论

模拟精度：模拟结果的精度受到空间步长、时间步长和吸收边界条件等因素的影响。选择合适的参数可以提高模拟精度。
方法扩展：本文利用二维FDTD方法成功地模拟了光波在玻璃-空气界面发生的折射现象。该方法可以扩展到三维情况，并用于模拟更复杂的介质结构和光学器件。

六、结论

本文利用二维FDTD方法成功地模拟了光波在玻璃-空气界面发生的折射现象，并通过彩色标度图直观地展示了折射过程。模拟结果验证了FDTD方法在光学模拟中的有效性，并展示了其在研究光波传播特性方面的应用潜力。该方法为理解和预测光学系统的性能提供了有力的工具。

📚2 运行结果

部分代码：

% Defining of the permittivity profile of the region:-

% Specifying the top half of the spatial domain which is made of glass
epsilon(:,1:ydim/2)=index*index*epsilon0;

% 2D FDTD update for PML boundary as used in previous program
if pml==1

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

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


%Perfectly matched layer boundary design
%[2]-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=20;

%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;

%Multiplication factor matrices for H matrix update to avoid being calculated many times
%in the time update loop so as to increase computation speed
G=((mu-0.5*deltat*sigma_starx)./(mu+0.5*deltat*sigma_starx));
H=(deltat/delta)./(mu+0.5*deltat*sigma_starx);
A=((mu-0.5*deltat*sigma_stary)./(mu+0.5*deltat*sigma_stary));
B=(deltat/delta)./(mu+0.5*deltat*sigma_stary);

%Multiplication factor matrices for E matrix update to avoid being calculated many times
%in the time update loop so as to increase computation speed
C=((epsilon-0.5*deltat*sigmax)./(epsilon+0.5*deltat*sigmax));
D=(deltat/delta)./(epsilon+0.5*deltat*sigmax);
E=((epsilon-0.5*deltat*sigmay)./(epsilon+0.5*deltat*sigmay));
F=(deltat/delta)./(epsilon+0.5*deltat*sigmay);

% Update loop begins
for n=1:1:time_tot

%matrix update instead of for-loop for Hy and Hx fields
Hy(1:xdim-1,1:ydim-1)=A(1:xdim-1,1:ydim-1).*Hy(1:xdim-1,1:ydim-1)+B(1:xdim-1,1:ydim-1).*(Ezx(2:xdim,1:ydim-1)-Ezx(1:xdim-1,1:ydim-1)+Ezy(2:xdim,1:ydim-1)-Ezy(1:xdim-1,1:ydim-1));
Hx(1:xdim-1,1:ydim-1)=G(1:xdim-1,1:ydim-1).*Hx(1:xdim-1,1:ydim-1)-H(1:xdim-1,1:ydim-1).*(Ezx(1:xdim-1,2:ydim)-Ezx(1:xdim-1,1:ydim-1)+Ezy(1:xdim-1,2:ydim)-Ezy(1:xdim-1,1:ydim-1));

%matrix update instead of for-loop for Ez field
Ezx(2:xdim,2:ydim)=C(2:xdim,2:ydim).*Ezx(2:xdim,2:ydim)+D(2:xdim,2:ydim).*(-Hx(2:xdim,2:ydim)+Hx(2:xdim,1:ydim-1));
Ezy(2:xdim,2:ydim)=E(2:xdim,2:ydim).*Ezy(2:xdim,2:ydim)+F(2:xdim,2:ydim).*(Hy(2:xdim,2:ydim)-Hy(1:xdim-1,2:ydim));

% Source condition incorporating given free space wavelength 'free_space_wavelength'
% and having a location at the left half of the top edge (in glass) of the domain just
% after the PML boundary with the plane of propagation of the emanating wave inclined
% at an angle 'alpha' with x-axis
N_lambda=free_space_wavelength/delta;
t_start=1;
alpha=35; %Angle 'alpha' should be < 42 degrees for refraction
i=75:1:125;