光子晶体90度弯曲波导的2D FDTD（Matlab代码实现）

长安程序猿

于 2024-10-31 08:46:53 发布

阅读量823

点赞数 27

文章标签： matlab 数据结构开发语言

本文链接：https://blog.csdn.net/Yan_she_He/article/details/143380254

版权

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

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

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

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

目录

⛳️赠与读者

💥1 概述

研究背景与意义

研究方法与步骤

预期结果与结论

后续研究方向

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

⛳️赠与读者

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

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

💥1 概述

光子晶体90度弯曲波导的2D FDTD研究文档

采用4.65微米x4.65微米的结构，该结构由7x7的方形周期性结构组成，其边长为0.15微米，在域的中心形成一个7x7的方形阵列，相邻结构之间在水平和垂直方向上的间距为0.45微米。结构之间和周围的区域充满了折射率为n=1的空气。从左侧3个周期性结构中去除一条中心水平线，然后去除3个周期性结构，形成从左侧到中心再到底部的90度弯曲波导，光通过该波导被引导。这里使用了Mur的吸收边界条件。

A 4.65 micron x 4.65 micron structure is taken with 7x7 periodic structures of square dimensions (made of silicon i.e n=3.4) having 0.15 micron as the length of a side as a 7x7 square array at the centre of the domain with a separation of 0.45 micron between neighbouring structure both horizontally and vertically. The region in between and surrounding that of the structures are filled with air of refractive index, n=1. A central horizontal line from left of 3 periodic structures are removed followed by removal of 3 periodic structures down to form a 90 degree bent wave guide from left to center to bottom through which the light is guided. Here Mur's absorbing boundary condition is used.

研究背景与意义

光子晶体作为一种由不同介质材料在空间呈周期性排列的新型光学结构，具有光子带隙和光子局域化等特性，可用于制作波导、逻辑门、谐振腔、光开关等高性能光学器件。二维光子晶体因其易于制备和调控光子带隙的优势，在集成光路方面具备很高的应用潜力。90度弯曲波导作为光子晶体器件之间互联的基础，其传输性能的优化对于提升整个光子晶体电路的性能至关重要。因此，采用2D-FDTD方法对二维光子晶体90度弯曲波导进行研究具有重要意义。

研究方法与步骤

模型建立
- 根据您提供的参数，建立一个4.65微米x4.65微米的二维光子晶体结构，该结构由7x7的方形周期性结构组成，边长为0.15微米。
- 在域的中心形成一个7x7的方形阵列，相邻结构之间在水平和垂直方向上的间距为0.45微米。
- 结构之间和周围的区域填充折射率为n=1的空气。
- 从左侧3个周期性结构中去除一条中心水平线，然后去除3个周期性结构，形成从左侧到中心再到底部的90度弯曲波导。
2D-FDTD算法应用
- 采用2D-FDTD方法对光子晶体90度弯曲波导进行模拟。
- 设置Mur的吸收边界条件以吸收向外传播的电磁波，减少反射。
- 使用Yee元胞结构来安排电场和磁场分量，以保证麦克斯韦方程组的因果性和稳定性。
- 时间步长和空间步长需要满足Courant稳定性条件。
参数优化与结果分析
- 通过调整光子晶体的晶格结构、介质柱截面面积、柱体方位角等参数，研究其对光子带隙宽度和透射特性的影响。
- 模拟不同频率入射光的传输特性，分析90度弯曲波导的光传输损耗和效率。
- 对比不同参数下的模拟结果，找出最优参数组合以实现低损耗、高效率的光传输。

预期结果与结论

通过2D-FDTD方法的模拟研究，预期能够得出以下结论：

光子晶体的晶格结构、介质柱截面面积、柱体方位角等参数对光子带隙宽度和透射特性有显著影响。
通过优化参数组合，可以实现低损耗、高效率的90度弯曲波导光传输。
2D-FDTD方法是一种有效的模拟手段，可用于指导二维光子晶体90度弯曲波导的设计与优化。

后续研究方向

进一步研究光子晶体90度弯曲波导的非线性效应和动态特性。
探索其他类型的二维光子晶体结构及其在不同应用场景下的性能表现。
将研究成果应用于光子晶体光开关、谐振腔等高性能光学器件的设计与制作中。

📚2 运行结果

部分代码：

%Movie type colour scaled image plot of Ez
h=imagesc(delta*(1:1:xdim)*1e+6,(1e+6*delta*(1:1:ydim))',Ez',[-1,1]);colorbar;
set(h,'AlphaData',2*epsilon'/epsilon0);
title(['\fontsize{20}Color-scaled image plot of Ez in photnic crystal 90 degree bend with PML boundary and at time = ',num2str(round(n*deltat*1e+15)),' fs']);
xlabel('x (in um)','FontSize',20);
ylabel('y (in um)','FontSize',20);
set(gca,'FontSize',20);
getframe;
end
end

%Fdtd update for Mur's absorbing boundary conditions
if abc==1

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

% Initializing electric and magnetic conductivity matrices
sigma=4e-4*ones(xdim,ydim);
sigma_star=4e-4*ones(xdim,ydim);

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

%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*sigma)./(epsilon+0.5*deltat*sigma));
D=(deltat/delta)./(epsilon+0.5*deltat*sigma);

%Mur's absorbing boundary condition parameters
p0=1;
p2=-0.5;

%Co-efficients of present and previous (regarding time-step) boundary Ez values
%in boundary update equation for forward/up and backward/down boundaries in
%the domain (as given in [1])
c0=(c/(2*S))*(1-(p0/S));
c1=-(c/(2*S))*(1+(p0/S));
c2=(c/(S^2))*(p0+(p2*S*S));
c3=-(p2*c)/2;
c0efffor=-(c0/c1);
c2efffor=-(c2/c1);
c3efffor=-(c3/c1);
c0=(c/(2*S))*(1+(p0/S));
c1=-(c/(2*S))*(1-(p0/S));
c2=-(c/(S^2))*(p0+(p2*S*S));
c3=(p2*c)/2;
c1effrev=-(c1/c0);
c2effrev=-(c2/c0);
c3effrev=-(c3/c0);

%Storage vectors for Ez boundary and boundary-1 values of previous and its
%previous timesteps
prev_xfor=zeros(1,ydim);
prev_x_minus_1for=zeros(1,ydim);
prev_yfor=zeros(xdim,1);
prev_y_minus_1for=zeros(xdim,1);
prev_xrev=zeros(1,ydim);
prev_x_minus_1rev=zeros(1,ydim);
prev_yrev=zeros(xdim,1);
prev_y_minus_1rev=zeros(xdim,1);

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

%Vector update instead of for-loop for Hy and Hx fields
Hx(2:xdim-3,2:ydim-3)=A(2:xdim-3,2:ydim-3).*Hx(2:xdim-3,2:ydim-3)-B(2:xdim-3,2:ydim-3).*(Ez(2:xdim-3,3:ydim-2)-Ez(2:xdim-3,2:ydim-3));
Hy(2:xdim-3,2:ydim-3)=A(2:xdim-3,2:ydim-3).*Hy(2:xdim-3,2:ydim-3)+B(2:xdim-3,2:ydim-3).*(Ez(3:xdim-2,2:ydim-3)-Ez(2:xdim-3,2:ydim-3));

%Vector update instead of for-loop for Ez field
Ez(3:xdim-3,3:ydim-3)=C(3:xdim-3,3:ydim-3).*Ez(3:xdim-3,3:ydim-3)+(Hy(3:xdim-3,3:ydim-3)-Hy(2:xdim-4,3:ydim-3)-Hx(3:xdim-3,3:ydim-3)+Hx(3:xdim-3,2:ydim-4)).*D(3:xdim-3,3:ydim-3);

%Mur's abc conditions obtained from Mur's difference equation for
%forward boundary
if n>1
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>1
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