oz.m

计算弹性波带隙结构

 %========================z mode of BaTiO3 and Polymer ============================     
    clear;
    clc;
    t0=clock;
    
    a=1;                        
    a1=a*[1,0,0];               
    a2=a*[0,1,0];   
    a3=a*[0,0,1];
       
    roub=1150;                  
    c44b=1.6*10^9;
    e15b=0;
    kec11b=0.039825*10^(-9);             

    roua=5800;                 
    c44a=43*10^9;
    e15a=11.6;
    kec11a=11.2*10^(-9);
    
    Rc=0.333779059;
    Au=norm(cross(a1,a2));
    Pf=(pi*Rc*Rc)/Au;
    
    ra1=(2*pi/a)*cross(a2,a3)/dot(a1,cross(a2,a3));
    ra2=(2*pi/a)*cross(a3,a1)/dot(a1,cross(a2,a3));
    ra1=ra1(1:2);
    ra2=ra2(1:2);
    
    NrSquare=8;                        
    G=zeros((2*NrSquare+1)^2,2);       
    ni=1;
    
    for L=-NrSquare:NrSquare
        for m=-NrSquare:NrSquare
            Glm=L*ra1+m*ra2;
            G(ni,:)=Glm;
            ni=ni+1;
        end
    end
    
    NG=ni-1;                               
    G=G(1:NG,:);    
    
    rou=zeros(NG,NG);
    c44=zeros(NG,NG);
    e15=zeros(NG,NG);
    kec11=zeros(NG,NG);
    
    for ni=1:NG
        for nj=1:NG
            Gij=norm(G(ni,:)-G(nj,:));
            if (Gij == 0)
                rou(ni,nj)=roua*Pf+roub*(1-Pf);
                c44(ni,nj)=c44a*Pf+c44b*(1-Pf);
                e15(ni,nj)=e15a*Pf+e15b*(1-Pf);   
                kec11(ni,nj)=kec11a*Pf+kec11b*(1-Pf);
            else
                rou(ni,nj)=(roua-roub)*Pf*2*besselj(1,Gij*Rc)/(Gij*Rc);
                c44(ni,nj)=(c44a-c44b)*Pf*2*besselj(1,Gij*Rc)/(Gij*Rc);
                e15(ni,nj)=(e15a-e15b)*Pf*2*besselj(1,Gij*Rc)/(Gij*Rc);        
                kec11(ni,nj)=(kec11a-kec11b)*Pf*2*besselj(1,Gij*Rc)/(Gij*Rc);
            end
        end
    end  
    
    T=(2*pi/a)*[1e-9,0];
    M=(2*pi/a)*[1/2,1/2];
    X=(2*pi/a)*[1/2,0];
    
    THETAaa=zeros(NG,NG);          
    THETAbb=zeros(NG,NG);       
    
    
    Nkpoints=15;                   
    stepsize=0:1/(Nkpoints-1):1;   
    
    TX_eig=zeros(Nkpoints,NG);   
    XM_eig=zeros(Nkpoints,NG);    
    MT_eig=zeros(Nkpoints,NG); 
    
    for n=1:Nkpoints
        TX_step=stepsize(n)*(X-T)+T;
        for ni=1:NG
            for nj=1:NG
                kGi=TX_step+G(ni,:);
                kGj=TX_step+G(nj,:);
                
                CC44(ni,nj)=c44(ni,nj)*(kGi(1,1)*kGj(1,1)+kGi(1,2)*kGj(1,2));
                EE15(ni,nj)=e15(ni,nj)*(kGi(1,1)*kGj(1,1)+kGi(1,2)*kGj(1,2));
                KKEC11(ni,nj)=kec11(ni,nj)*(kGi(1,1)*kGj(1,1)+kGi(1,2)*kGj(1,2));
                
                THETAbb(ni,nj)=rou(ni,nj);
            end
        end
        
        THETAaa=CC44+EE15*inv(KKEC11)*EE15;        
        [VTX,DTX]=eig(THETAaa,THETAbb);
        
        for ni=1:NG
            for nj=1:NG
                if (nj == ni)
                    DDDTX(1,ni)=DTX(ni,nj);
                end
            end
        end
        
        TX_eig(n,:)=sort(sqrt(DDDTX));
 
        
        XM_step=stepsize(n)*(M-X)+X;
        for ni=1:NG
            for nj=1:NG
                kGi=XM_step+G(ni,:);
                kGj=XM_step+G(nj,:);
                
                CC44(ni,nj)=c44(ni,nj)*(kGi(1,1)*kGj(1,1)+kGi(1,2)*kGj(1,2));
                EE15(ni,nj)=e15(ni,nj)*(kGi(1,1)*kGj(1,1)+kGi(1,2)*kGj(1,2));
                KKEC11(ni,nj)=kec11(ni,nj)*(kGi(1,1)*kGj(1,1)+kGi(1,2)*kGj(1,2));
                
                THETAbb(ni,nj)=rou(ni,nj);
            end
        end
        
        THETAaa=CC44+EE15*inv(KKEC11)*EE15;        
        [VXM,DXM]=eig(THETAaa,THETAbb);
        
        for ni=1:NG
           for nj=1:NG
               if (nj == ni)
                   DDDXM(1,ni)=DXM(ni,nj);
               end
           end
       end
       
       XM_eig(n,:)=sort(sqrt(DDDXM));

       
       MT_step=stepsize(n)*(T-M)+M;
       
       for ni=1:NG
           for nj=1:NG
               kGi=MT_step+G(ni,:);
               kGj=MT_step+G(nj,:); 
               
                CC44(ni,nj)=c44(ni,nj)*(kGi(1,1)*kGj(1,1)+kGi(1,2)*kGj(1,2));
                EE15(ni,nj)=e15(ni,nj)*(kGi(1,1)*kGj(1,1)+kGi(1,2)*kGj(1,2));
                KKEC11(ni,nj)=kec11(ni,nj)*(kGi(1,1)*kGj(1,1)+kGi(1,2)*kGj(1,2));
                
                THETAbb(ni,nj)=rou(ni,nj);
            end
        end
        
        THETAaa=CC44+EE15*inv(KKEC11)*EE15;        
        [VMT,DMT]=eig(THETAaa,THETAbb);
       
       for ni=1:NG
           for nj=1:NG
               if (nj == ni)
                   DDDMT(1,ni)=DMT(ni,nj);
               end
           end
       end
       
       MT_eig(n,:)=sort(sqrt(DDDMT));
   end   

   
   kaxis=0;          
   TXaxis=kaxis:norm(T-X)/(Nkpoints-1):(kaxis+norm(T-X));
   kaxis=kaxis+norm(T-X);
   XMaxis=kaxis:norm(X-M)/(Nkpoints-1):(kaxis+norm(X-M));
   kaxis=kaxis+norm(X-M);
   MTaxis=kaxis:norm(X-T)/(Nkpoints-1):(kaxis+norm(X-T));
   kaxis=kaxis+norm(X-T);

   
   Ntraject=3;           
   EigFreq=zeros(Ntraject*Nkpoints,1);
  
   figure(1);
   hold on
   
   Nk=Nkpoints;
   
   for k=1:NG
       for ni=1:Nkpoints
           EigFreq(ni+0*Nk)=TX_eig(ni,k)/(2*pi*1180/a);
           EigFreq(ni+1*Nk)=XM_eig(ni,k)/(2*pi*1180/a);
           EigFreq(ni+2*Nk)=MT_eig(ni,k)/(2*pi*1180/a);
       end   
       plot(TXaxis(1:Nk),EigFreq(1+0*Nk:1*Nk),XMaxis(1:Nk),EigFreq(1+1*Nk:2*Nk),MTaxis(1:Nk),EigFreq(1+2*Nk:3*Nk));
   end     
   box on   
   grid on