• Fb_
    了解作者
  • matlab
    开发工具
  • 2KB
    文件大小
  • zip
    文件格式
  • 1
    收藏次数
  • 1 积分
    下载积分
  • 9
    下载次数
  • 2018-12-21 15:17
    上传日期
利用特征模匹配方法计算二维声子晶体的能带图,流体和固体组合
EMMT_fluid_solid.zip
  • EMMT_fluid_solid
  • bulkchen.m
    4.1KB
  • eigenvector.m
    2.8KB
内容介绍
clc clear D_inc=8400; %brass C11_inc=16.262e10; C12_inc=8.132e10; C44_inc=4.065e10; D_mat=1000; %water C11_mat=2.0535e9; D_inc=D_inc/D_mat; C11_inc=C11_inc/C11_mat; C12_inc=C12_inc/C11_mat; C44_inc=C44_inc/C11_mat; c11m=1.0; rom=1.0; % C12_mat=0; C44_mat=0; LceLL=1; N=4;%倒格子数目:(2*N+1)*(2*N+1) nn=2*N+1; nnt=2*nn; Kx=0; vration=C44_mat/C11_mat; mm=1000; m=100; vration1=C44_inc/C11_inc; for WW=1:mm W=WW/m; hhh=0; if (W-fix(W)==0) hhh=hhh+1;continue;end%这个该如何理解 [muu,mud,mtu,mtd,kyu,kyd]=isof(N,W,rom,Kx);%计算同种材料特征值(流体) fr=1/4; [CUR,CUL,CTR,CTL,ckyr,ckyl]=eigf(N,nn,nnt,D_inc,C11_inc,W,fr,Kx,rom,c11m); [aa11,aa12,aa21,aa22]=intesmat(muu,mud,mtu,mtd,CUR,CUL,CTR,CTL);%交界面 frv=0.5/4; [bb11,bb12,bb21,bb22]=intra(N,frv,ckyr,ckyl);%内部传播 [s11,s12,s21,s22]=matlu(aa11,aa12,aa21,aa22,bb11,bb12,bb21,bb22,N); fr=1; % [CUR1,CUL1,CTR1,CTL1,ckyr1,ckyl1]=eigenvector(N,D_mat,D_inc,fr,C11_mat,C11_inc,C12_mat,C12_inc,C44_mat,C44_inc,Kx,W); [CUR1,CUL1,CTR1,CTL1,ckyr1,ckyl1]=isovector(N,W,vration1,LceLL,C12_inc,C11_inc,C44_inc); [aa11,aa12,aa21,aa22]=intesmat(CUR,CUL,CTR,CTL,CUR1,CUL1,CTR1,CTL1);%交界面 frv=1/4; [bb11,bb12,bb21,bb22]=intra(N,frv,ckyr1,ckyl1);%内部传播 [ss11,ss12,ss21,ss22]=matlu(aa11,aa12,aa21,aa22,bb11,bb12,bb21,bb22,N); [s11,s12,s21,s22]=matlu(s11,s12,s21,s22,ss11,ss12,ss21,ss22,N); [US11,US12,US21,US22]=intesmat(CUR1,CUL1,CTR1,CTL1,muu,mud,mtu,mtd);%边界 % [US11,US12,US21,US22]=intesmat(CUR,CUL,CTR,CTL,muu,mud,mtu,mtd);%边界 fr=2.5/4; [b11,b12,b21,b22]=intra(N,fr,kyu,kyd); [ss11,ss12,ss21,ss22]=matlu(US11,US12,US21,US22,b11,b12,b21,b22,N); [s11,s12,s21,s22]=matlu(s11,s12,s21,s22,ss11,ss12,ss21,ss22,N); q11=eye(2*N+1); q22=eye(2*N+1); q12=zeros(2*N+1); q21=zeros(2*N+1); for ii=1:1 [q11,q12,q21,q22]=matlu(q11,q12,q21,q22,s11,s12,s21,s22,N); end for ii=1:(2*N+1) ax=abs(kyu(ii))-W; if(abs(ax)<1.0d-10) inn=ii; end end for ii=1:(2*N+1) ax=abs(kyu(ii))-W*(vration)^0.5; if(abs(ax)<1.0d-10) innu=ii; end end utra1=muu*q11(:,innu); %横波 ttra1=mtu*q11(:,innu); uref1=mud*q21(:,innu); tref1=mtd*q21(:,innu); uin1=muu(:,innu); tin1=mtu(:,innu); utra2=muu*q11(:,inn); %纵波 ttra2=mtu*q11(:,inn); uref2=mud*q21(:,inn); tref2=mtd*q21(:,inn); uin2=muu(:,inn); tin2=mtu(:,inn); ref1=0; in1=0; tra1=0; ref2=0; in2=0; tra2=0; for k=1:(2*N+1) ref1=ref1+abs(real(conj(uref1(k))*tref1(k))); in1=in1+abs(real(conj(uin1(k))*tin1(k))); tra1=tra1+abs(real(conj(utra1(k))*ttra1(k))); ref2=ref2+abs(real(conj(uref2(k))*tref2(k))); in2=in2+abs(real(conj(uin2(k))*tin2(k))); tra2=tra2+abs(real(conj(utra2(k))*ttra2(k))); end T1(WW)=tra1/in1; R1(WW)=ref1/in1; T2(WW)=tra2/in2; R2(WW)=ref2/in2; TT(WW)=T1(WW)+T2(WW); end % hold on; % subplot(2,1,1);plot(T1,'-') % subplot(2,1,2);plot(T2,'-') % plot((1:200)/100*1433/4e-3,T1) for k=1:mm-hhh % fluidtran(k)=(k/100*(500)^0.5/4E-3)/1e3; fluidtran(k)=(C11_mat/1e3)^0.5*k/m/4e-3/1000; end c=[fluidtran;T2]; plot(c(1,:),c(2,:)) % plot(fluidtran,x0,'LineWidth', 2) xlabel('KHz','FontName','Times New Roman'); ylabel('tr','FontName','Times New Roman');%标注单位 % !title(['XYZ mode ',date]);%标注标题日期 % axis([0, 4*pi, 0, 200]); axis([0, 1000, 0, 1]);
评论
    相关推荐