住房城乡建设干部学院网站,广西茶叶网站建设,软件开发五个阶段,电影免费在线观看Matlab 平面波展开法计算二维声子晶体二维声子晶体带结构计算#xff0c;材料是铅柱在橡胶基体中周期排列#xff0c;格子为正方形。采用PWE方法计算 完整程序:
%%%%%%%%%%%%%%%%%%%%%%%%% clear;clc;tic;epssys1.0e-6; %设定一个最小量#xff0c;避免系统截断误差或除零错…Matlab 平面波展开法计算二维声子晶体二维声子晶体带结构计算材料是铅柱在橡胶基体中周期排列格子为正方形。采用PWE方法计算 完整程序:
%%%%%%%%%%%%%%%%%%%%%%%%% clear;clc;tic;epssys1.0e-6; %设定一个最小量避免系统截断误差或除零错误 %%%%%%%%%%%%%%%%%%%%%%%%%%
%定义实际的正空间格子基矢 %%%%%%%%%%%%%%%%%%%%%%%%%% a0.02; a1a*[1 0]; a2a*[0 1]; %%%%%%%%%%%%%%%%%%%%%%%%%%
%定义晶格的参数 %%%%%%%%%%%%%%%%%%%%%%%%%% rho111600;E14.08e10;mju11.49e10;lambda1mju1*(E1-2*mju1)/(3*mju1-E1); %散射体的材料参数 rho21300;E21.175e5;mju24e4;lambda2mju2*(E2-2*mju2)/(3*mju2-E2); %基体的材料参数 Rc0.006; %散射体截面半径 Acpi*(Rc)^2; %散射体截面面积 Aua^2; %二维格子原胞面积 PfAc/Au; %填充率 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%生成倒格基矢 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% b12*pi/a*[1 0]; b22*pi/a*[0 1]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %选定参与运算的倒空间格矢量即参与运算的平面波数量 %设定一个l,m的取值范围变化l,m即可得出参与运算的平面波集合 NrSquare10; %选定倒空间的尺度即l,m倒格矢Gl*b1m*b2的取值范围。 %NrSquare确定后使用Bloch波数目可能为(2*NrSquare1)^2 Gzeros((2*NrSquare1)^2,2); %初始化可能使用的倒格矢矩阵 i1; for l-NrSquare:NrSquare for m-NrSquare:NrSquare G(i,:)l*b1m*b2; ii1; end; end; NGi-1; %实际使用的Bloch波数目 GG(1:NG,:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %生成k空间的rho(Gi-Gj),mju(Gi-Gj),lambda(Gi-Gj)值i,j从1到NG。 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% rhozeros(NG,NG);mjuzeros(NG,NG);lambdazeros(NG,NG); for i1:NG for j1:NG Gijnorm(G(j,:)-G(i,:)); if (Gijepssys) rho(i,j)rho1*Pfrho2*(1-Pf); mju(i,j)mju1*Pfmju2*(1-Pf); lambda(i,j)lambda1*Pflambda2*(1-Pf); else rho(i,j)(rho1-rho2)*2*Pf*besselj(1,Gij*Rc)/(Gij*Rc); mju(i,j)(mju1-mju2)*2*Pf*besselj(1,Gij*Rc)/(Gij*Rc); lambda(i,j)(lambda1-lambda2)*2*Pf*besselj(1,Gij*Rc)/(Gij*Rc); end; end; end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %定义简约布里渊区的各高对称点 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% T(2*pi/a)*[epssys 0]; M(2*pi/a)*[1/2 1/2]; X(2*pi/a)*[1/2 0]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %对于简约布里渊区边界上的每个k求解其特征频率 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% THETA_Azeros(NG,NG); %待解的本征方程A矩阵 THETA_Bzeros(NG,NG); %待解的本征方程B矩阵 Nkpoints10; %每个方向上取的点数 stepsize0:1/(Nkpoints-1):1; %每个方向上步长 TX_eigzeros(Nkpoints,NG); %沿TX方向的波的待解的特征频率矩阵 XM_eigzeros(Nkpoints,NG); %沿XM方向的波的待解的特征频率矩阵 MT_eigzeros(Nkpoints,NG); %沿MT方向的波的待解的特征频率矩阵 for n1:Nkpoints fprintf([\n k-point:,int2str(n),of,int2str(Nkpoints),.\n]); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %对于TX正方格子方向上的每个k值求解其特征频率 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TX_stepstepsize(n)*(X-T)T; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %n 求本征矩阵的元素 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i1:NG for j1:NG kGiTX_stepG(i,:); kGjTX_stepG(j,:); THETA_A(i,j)mju(i,j)*dot(kGi,kGj); THETA_B(i,j)rho(i,j); end; end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %求解TX正方格子方向上的k矩阵的特征频率 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TX_eig(n,:)sort(sqrt(eig(THETA_A,THETA_B))).; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %对于XM正方格子方向上的每个k值求解其特征频率 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% XM_stepstepsize(n)*(M-X)X; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %n 求本征矩阵的元素 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i1:NG for j1:NG kGiXM_stepG(i,:); kGjXM_stepG(j,:); THETA_A(i,j)mju(i,j)*dot(kGi,kGj); THETA_B(i,j)rho(i,j); end; end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %求解XM正方格子方向上的k矩阵的特征频率 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% XM_eig(n,:)sort(sqrt(eig(THETA_A,THETA_B))).; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %对于MT正方格子方向上的每个k值求解其特征频率 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MT_stepstepsize(n)*(T-M)M; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %n 求本征矩阵的元素 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i1:NG for j1:NG kGiMT_stepG(i,:); kGjMT_stepG(j,:); THETA_A(i,j)mju(i,j)*dot(kGi,kGj); THETA_B(i,j)rho(i,j); end; end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %求解MT正方格子方向上的k矩阵的特征频率 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MT_eig(n,:)sort(sqrt(eig(THETA_A,THETA_B))).; end; fprintf(\n Calculation Time:%d sec,toc); save pbs2D %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %绘制声子晶体能带结构图 %首先将特定方向正方格子TX,XM,MT离散化 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% kaxis0; TXaxiskaxis:norm(T-X)/(Nkpoints-1):(kaxisnorm(T-X)); kaxiskaxisnorm(T-X); XMaxiskaxis:norm(M-X)/(Nkpoints-1):(kaxisnorm(X-M)); kaxiskaxisnorm(X-M); MTaxiskaxis:norm(T-M)/(Nkpoints-1):(kaxisnorm(T-M)); kaxiskaxisnorm(T-M); Ntraject3; %所需绘制的特定方向的数目 EigFreqzeros(Ntraject*Nkpoints,1); figure(1) hold on; NkNkpoints; for k1:NG for i1:Nkpoints EigFreq(i0*Nk)TX_eig(i,k)/(2*pi); EigFreq(i1*Nk)XM_eig(i,k)/(2*pi); EigFreq(i2*Nk)MT_eig(i,k)/(2*pi); end; plot(TXaxis(1:Nk),EigFreq(10*Nk:1*Nk),b,... XMaxis(1:Nk),EigFreq(11*Nk:2*Nk),b,... MTaxis(1:Nk),EigFreq(12*Nk:3*Nk),b); end; grid on; hold off; titlestr传统平面波展开法计算得到的二维声子晶体能带结构图; title(titlestr); xlabel(波矢k); ylabel(频率f/Hz); axis([0 MTaxis(Nkpoints) 0 800]); set(gca,XTick,[TXaxis(1) TXaxis(Nkpoints) XMaxis(Nkpoints) MTaxis(Nkpoints)]); xtixlabelchar(T,X,M,T); set(gca,XTickLabel,xtixlabel);
