📄 f41.m
字号:
clear
%定义常量和变量
tmax=0.5e-6;
dx=0.1e-3;
dy=dx;
dt=0.5e-8;
maxgridsx=8; %x方向的网格数
maxgridsy=8; %y方向的网格数
k=0;
%金属材料常数
% media[1][1]=2799; %金属密度
% media[2][1]=0.33; %mu
% media[3][1]=0.5 %lamda
c=6342; %声速
%基体材料常数
% media[1][2]=1142; %基体的密度
% media[2][2]=0.28 %mu
% media[3][2]=0.0.48 %lamda
media=[2799 1142;
0.33 0.28;
0.48 0.46];
%初始化数组
ux=zeros(8,8,10);
uy=zeros(8,8,10);
txx1=zeros(8,8,10);
txx2=ux;
txy1=ux;
txy2=ux;
txy3=ux;
tyy1=ux;
tyy2=ux;
vx=ux;
vy=ux;
kmax=maxgridsy/4;
%定义散射体及基体
for m=1:maxgridsx
for n=1:maxgridsy
if ((mod(m-2,4)==0)|(mod(m-3,4)==0)|(mod(m-4,4)==0))&&((mod(n-2,4)==0)|(mod(n-3,4)==0)|(mod(n-4,4)==0))
media1(m,n)=2799;
media2(m,n)=0.33;
media3(m,n)=0.44;
else
media1(m,n)=1142;
media2(m,n)=0.28;
media3(m,n)=0.46;
end
end
end
disp('media1');
disp(media1);
pause(3);
disp('media2');
disp(media2);
pause(3);
disp('media3');
disp(media3);
pause(3);
disp('keyboard----1');
disp('k');
disp(k);
pause(3);
disp('media1(2,2)');
disp(media1(2,2));
pause(3);
%keyboard;
coef1=dt*dt/dx;
coef2=dt*dt/dy;
mur_coef1=(c*dt-dy)/(c*dt+dy);
mur_coef2=mur_coef1;
disp('coef1,coef2');
disp(coef1);
disp(coef2);
pause(3);
% keyboard;
for k=2:10
for n=1:maxgridsy
uy(1,n,k)=1000*sin(2*3.14*2000000*(k)*dt); %-cos(2*3.14*500000*k*dt)*exp(4*3.14*(k*dt-5.0e-6)*(k*dt-5.0e-6)/16e-12); 修改为点激励！ 再注意k的取值！
end
disp('*********************============激励uy=====0');
disp(uy);
disp('*********************============激励uy=====00');
pause(2);
% keyboard;
for m=1:maxgridsx-1
for n=2:maxgridsy-1
txx1(m,n,k)=(media3(m,n)+2*media2(m,n))*(ux(m+1,n,k)-ux(m,n,k))/dx+media3(m,n)*(uy(m,n,k)-uy(m,n-1,k))/dy; %txx1
end
end
% disp('k==');
% disp(k);
disp('*********************=================txx1==1');
disp(txx1);
disp('*********************=================txx1==1');
pause(2);
% keyboard;
for m=2:maxgridsx-1
for n=1:maxgridsy-1
txy1(m,n,k)=media2(m,n)*(ux(m,n+1,k)-ux(m,n,k))/dy+(uy(m,n,k)-uy(m-1,n,k))/dx; %txy1
end
end
disp('*********************================txy1==2');
disp(txy1);
disp('*********************================txy1==2');
pause(2);
% keyboard;
for m=1:maxgridsx-1
for n=2:maxgridsy-1
tyy1(m,n,k)=(media3(m,n)+2*media2(m,n))*(uy(m,n,k)-uy(m,n-1,k))/dy+media3(m,n)*(ux(m+1,n,k)-ux(m,n,k))/dx; %tyy1
end
end
disp('*********************==================tyy1===3');
disp(tyy1);
disp('*********************==================tyy1===3');
pause(2);
% keyboard;
for m=2:maxgridsx-1
for n=2:maxgridsy-1
txx2(m,n,k)=(media3(m,n)+2*media2(m,n))*(ux(m,n,k)-ux(m-1,n,k))/dx+media3(m,n)*(uy(m-1,n,k)-uy(m-1,n-1,k))/dy; %txx2 注意避免和txx1相等，考虑边界设定问题！
end
end
disp('*********************==================txx2====4');
disp(txx2);
disp('*********************==================txx2====4');
pause(2);
% keyboard;
for m=1:maxgridsx-1 %**************
for n=1:maxgridsy-1
txy3(m+1,n,k)=media2(m,n)*((ux(m+1,n+1,k)-ux(m+1,n,k))/dy+(uy(m+1,n,k)-uy(m,n,k))/dx); %txy3
end
end
disp('*********************===============txy3=======5');
disp(txy3);
disp('*********************===============txy3=======5');
pause(2);
% keyboard;
for m=1:maxgridsx-1
for n=1:maxgridsy-1
tyy2(m,n+1,k)=(media3(m,n+1)+2*media2(m,n+1))*(uy(m,n+1,k)-uy(m,n,k))/dy+media3(m,n+1)*(ux(m+1,n+1,k)-ux(m,n+1,k))/dx; %tyy2
end
end
disp('*********************=============tyy2=========6');
disp(tyy2);
disp('*********************=============tyy2=========6');
pause(2);
% keyboard;
for m=2:maxgridsx-1
for n=2:maxgridsy-1
txy2(m,n,k)=-txy1(m,n,k)+media2(m+1,n)*((ux(m,n,k)+ux(m,n+1,k))/dy+(uy(m,n,k)+uy(m,n-1,k)-(uy(m-1,n,k)+uy(m-1,n-1,k)))/dx); %txy2
end
end
disp('*********************=============txy2=========7');
disp(txy2);
disp('*********************=============txy2=========7');
pause(2);
% keyboard;
for m=1:maxgridsx-1
for n=1:maxgridsy-1
ux(m,n,k+1)=2*ux(m,n,k)-ux(m,n,k-1)+coef1*(txx1(m,n,k)-txx2(m,n,k))/media1(m,n)+coef2*(txy1(m,n,k)-txy2(m,n,k))/media1(m,n); %ux
uy(m,n,k+1)=2*uy(m,n,k)-uy(m,n,k-1)+coef1*(txy3(m+1,n,k)-txy1(m,n,k))/media1(m,n)+coef2*(tyy2(m,n+1,k)-tyy1(m,n,k))/media1(m,n); %uy
end
end
disp('*********************===============ux======8');
disp(ux);
disp('*********************===============ux======80');
pause(2);
% keyboard;
disp('*********************===============uy======9');
disp(uy);
disp('*********************===============uy======90');
pause(2);
keyboard;
% for m=1:maxgridsx % 边界问题****************这个问题没有解决！
% for n=1:maxgridsy
% uy(m,maxgridsy,k+1)=uy(m,maxgridsy,k)+mur_coef1*(uy(m,maxgridsy,k+1)-uy(m,n,k));
% end
% end
% disp('*********************=10');
% disp('ux=========================');
% disp(ux(m,maxgridsy,k));
% pause(3);
% disp('uy===================================');
%disp(uy(m,maxgridsy,k));
%pause(3);
% disp('ux=========================');
% disp(uy);
% pause(3);
% for k=1:5
% for m=1:maxgridsx
% for n=1:maxgridsy
% vx=(ux(m,n,k+1)-ux(m,n,k))/dt;
% vy=(uy(m,n,k+1)-uy(m,n,k))/dt;
% end
% end
% end
% disp('vx=====================****');
% disp(vx);
% disp('vy=====================****');
%disp(vy);
disp('k==');
disp(k);
disp('最后');
% disp('ux=========================');
% disp(ux);
% pause(5);
keyboard;
end
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -