MATLAB 3D plot on dispersion curve

10 次查看（过去 30 天）
yusuf
yusuf 2015-12-9
回答： MANYINGA MARC 2023-4-12
I am working on wave propagation in multilayer structures. I wrote a code for 2-D plotting dispersion curves at specific frequency (j=1:1 condition) which is presented at the end of the question. My question is that how I can plot those curves at different frequencies (for example, j=1:3) on the same plotting. I think that it should be 3-D plotting and third axes should be frequency. Can anyone give me an idea or source for this condition ?
Thank you,
clc;
clear all;
% System Definition
N=3; %Number of layers including half spaces N:1 top half space N:bottom half space
G=zeros((N-1)*4,(N-1)*4);
%Thickness
th(1)=1000;%mm thickness of the top half space arbitrary meters
th(2)=1E-3;%mm thickness of the layer meters
th(N)=1000;%mm thickness of the layer meters
%th(4)=1E-3;%mm thickness of the layer meters
%th(N)=0;%mm thickness of the bottom half space arbitrary meters
poi=0.3;
coeff=[0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 2 3];
mass_coeff=[0.2 0.3 .4 .5 .6 .7 .8 .9 1 2];
for ii=1:10
for jj=1:11
%Material Props
%Layer 1
rho(1)=1000; %(if vacuum 0)
alfa(1)=1480; %(if vacuum 1000)%m/s longitidunal velocity
beta(1)=1; %(if vacuum 1000)%m/s shear velocity %beta=sqrt(mu/rho) %lame constant def %beta=sqrt((E)/(2*rho*(1+v)))/1000; %E and v definition
%Layer 2
rho(2)=rho(1)*mass_coeff(ii); %(if vacuum 0)
alfa(2)=alfa(1)*coeff(jj); %(if vacuum 1000)%m/s longitidunal velocity
alfa22(jj)=alfa(2);
beta(2)=alfa(2)/(sqrt((1-2*poi)/(2*(1-poi))))^-1; %(if vacuum 1000)%m/s shear velocity %beta=sqrt(mu/rho) %lame constant def %beta=sqrt((E)/(2*rho*(1+v)))/1000; %E and v definition
%Layer 2
rho(3)=rho(1); %(if vacuum 0)
alfa(3)=alfa(1); %(if vacuum 1000)%m/s longitidunal velocity
beta(3)=beta(1); %(if vacuum 1000)%m/s shear velocity %beta=sqrt(mu/rho) %lame constant def %beta=sqrt((E)/(2*rho*(1+v)))/1000; %E and v definition
% %Layer N
% rho(4)=rho(2); %(if vacuum 0)
% alfa(4)=alfa(2); %(if vacuum 1000)%m/s longitidunal velocity
% beta(4)=beta(2); %(if vacuum 1000)%m/s shear velocity %beta=sqrt(mu/rho) %lame constant def %beta=sqrt((E)/(2*rho*(1+v)))/1000; %E and v definition
%
% %Layer N
% rho(5)=rho(1); %(if vacuum 0)
% alfa(5)=alfa(1); %(if vacuum 1000)%m/s longitidunal velocity
% beta(5)=beta(1); %(if vacuum 1000)%m/s shear velocity %beta=sqrt(mu/rho) %lame constant def %beta=sqrt((E)/(2*rho*(1+v)))/1000; %E and v definition
%Velocity Range
Vstop=2000; %Top Vel
Vstart=10; %Bottom Vel
f_step=5E5; %frequency step
for j=1:1
%Function Vstart
f(j)=j*f_step;
w=2*pi*f(j);
for i=1:(Vstop-Vstart);
cp=Vstart+i-1;%km/s
k=w/cp;
for r=1:N
Ca=sqrt(((w^2)/(alfa(r)^2))-k^2);
Cb=sqrt(((w^2)/(beta(r)^2))-k^2);
B=(w^2)-(2*(beta(r)^2)*(k^2));
ga=exp(1i*th(r)*sqrt(((w^2)/(alfa(r)^2))-k^2));
gb=exp(1i*th(r)*sqrt(((w^2)/(beta(r)^2))-k^2));
%Wave vectors
L_plus=[k; Ca; B*1i*rho(r); 2*1i*rho(r)*k*(beta(r)^2)*Ca];
S_plus=[Cb; -k; -2*1i*rho(r)*k*(beta(r)^2)*Cb; B*1i*rho(r)];
L_minus=[k; -Ca; B*1i*rho(r); -2*1i*rho(r)*k*(beta(r)^2)*Ca];
S_minus=[-Cb; -k; 2*1i*rho(r)*k*(beta(r)^2)*Cb; B*1i*rho(r)];
a=-mod(r,2)*2+1;
if (r==1)
G(1:4,1:2)=[-L_minus -S_minus];
elseif (r==N)
G((N-1)*4-3:(N-1)*4,(N-1)*4-1:(N-1)*4)=[a*L_plus a*S_plus];
else
G((r-1)*4-3:(r-1)*4,(r-1)*4-1:(r-1)*4+2)=[a*L_plus a*S_plus a*L_minus a*S_minus];
G((r-1)*4+1:(r-1)*4+4,(r-1)*4-1:(r-1)*4+2)=[(a*ga)*L_plus (a*gb)*S_plus (a/ga)*L_minus (a/gb)*S_minus];
end
end
Deter_G_mag(i)=abs(det(G));
end
v_range=Vstart:(Vstop-1);
M(j,:)=Deter_G_mag;
n=1;
for m=2:(Vstop-Vstart-1)
if (m==2)
der_det(m,j)=(Deter_G_mag(m+1)-Deter_G_mag(1));
else if (m==(Vstop-Vstart-1))
der_det(m,j)=(Deter_G_mag(Vstop-Vstart-1)-Deter_G_mag(m-1));
else
der_det(m,j)=(Deter_G_mag(m+1)-Deter_G_mag(m-1));
end
end
if (sign(der_det(m-1,j))==-1 && sign(der_det(m,j))==1 && abs(v_range(m)-alfa(2))>1 && abs(v_range(m)-alfa(1))>1 && abs(v_range(m)-alfa(3))>1 && abs(v_range(m)-beta(1))>1 && abs(v_range(m)-beta(2))>1 && abs(v_range(m)-beta(3))>1);
V_p(n,j)=v_range(m);
n=n+1;
end
end
end
Vpp(jj)=V_p(1,1);
end
ratio=Vpp./alfa22;
hold on
plot(coeff,ratio,'*');
end
  1 个评论
Ghanem Alatteili
Ghanem Alatteili 2022-7-6
Hi,
I am working something similar. I am trying to plot a dispersion relation by computiong dipole intaeraction matrix. Then I will insert the translation vectors to get the dipersion relation. I would like to know how do you define your k as a number or a function. I might help.
Thank you,

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

MANYINGA MARC
MANYINGA MARC 2023-4-12

Hi everyone i can bring you some solution about your problem. My code is base of Stiffness matrix method of Rokhlin and Wang. Here i try to calculate reflexion and transmission coefficient for a multilayer a a function of frequency/incudent angle. My code didn't give me a result what i expected but it can help you for your Problem. And you can also help me to improve it

if true
  % code
end

类别

在 Help Center 和 File Exchange 中查找有关 Contour Plots 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by