Finite element time domain method for metamaterials

clear all a=150e-3; % 晶格常数 rhoA=2730; % A铝的密度 rhoB=1180; % B环氧树脂的密度 EA=7.76e10; % A铝的杨氏模量 EB=4.35e9; % B环氧树脂的杨氏模量 cA=(EA/rhoA)^(1/2); % A铝中的声速 cB=(EB/rhoB)^(1/2); % B环氧树脂中的声速 NP=5.5; % fmax=3e4; % 最大频率 dx=5e-3; % 空间步长 dt=5e-7; % 时间步长 nx=a/dx; % 一个晶格中的网格数 Nx=NP*nx+1; % 划分的网格总数 E=EB*ones(1,Nx); % 1行Nx列的EB 定义所有网格的杨氏模量为树脂的 rho=rhoB*ones(1,Nx); % 1行Nx列的rhoB 定义所有网格的密度为树脂的 for ii=0:(NP-1) % 定义声子晶体 for jj=(nx/2+2):nx E(jj+ii*nx)=EA; rho(jj+ii*nx)=rhoA; end end for ii=1:(NP*2-1) % 定义两种材料交界处密度与杨氏模量为平均值 E(ii*nx/2+1)=(EB+EA)/2; rho(ii*nx/2+1)=(rhoB+rhoA)/2; end x0=10; nt=1e5; % 时间有100000步 t0=40; spread=1/fmax/dt/sqrt(pi); v0=zeros(1,Nx); % 定义1行Nx列0数组 tau0=zeros(1,Nx); v=zeros(1,Nx); tau=zeros(1,Nx); outputpoint=Nx; output=zeros(1,nt); t=1:nt; source=exp(-((t-t0)/spread).^2); % t从1到nt，信号的数组 for ii=1:nt for jj=2:(Nx-1) v(jj)=v0(jj)+1/rho(jj)*dt/dx*(tau0(jj)-tau0(jj-1)); end v(1)=v0(2)+(cB*dt-dx)/(cB*dt+dx)*(v(2)-v0(1)); v(Nx)=v0(Nx-1)+(cB*dt-dx)/(cB*dt+dx)*(v(Nx-1)-v0(Nx)); v(x0)=v(x0)+source(ii); for jj=2:(Nx-1) tau(jj)=tau0(jj)+E(jj)*dt/dx*(v(jj+1)-v(jj)); end tau(1)=tau0(2)+(cB*dt-dx)/(cB*dt+dx)*(tau(2)-tau0(1)); tau(Nx)=tau0(Nx-1)+(cB*dt-dx)/(cB*dt+dx)*(tau(Nx-1)-tau0(Nx)); v0=v; tau0=tau; output(ii)=v0(outputpoint); end T=nt*dt; N=round(T*fmax); FDoutput=dt*fft(output); FDsource=dt*fft(source); FDoutput=abs(FDoutput(1:N)); FDsource=abs(FDsource(1:N)); f=(1:N)/T/1e3; FR=20*log10(FDoutput./FDsource); plot(f,FR) xlabel('Frequency/kHz'); ylabel('Amplitude/dB');