Matlab绘制频散曲线程序代码(20210119130722)

Matlab绘制频散曲线

程序代码

-CAL-FENGHAI-(2020YEAR-YI

function disper

%绘制平板频散曲线

%tic

clc;clear;

cl=5790;%材料纵波波速(钢板)

cs=3200;%材料横波波速(钢板) dfd=*le3;

fdO=:dfd/le3:2O)*le3;%频厚积(MHz*mm) d_Q235二6;

cps_mi n二2700;

cpa_min=100;

cp_max=10000;

mode=3;%绘制的模式数

precision=le-8;

cpa=zeros(length(fdO),mode);

cps=zeros(le ng th(fdO),mode);

for i=l:length(fdO)

fd=fdO(i);

[cpl2 n]=ss(cps_min/cp_max/fd/cl,cs,mode);

for j=l:n

cpl=cpl2(j,l);

cp2=cpl2(j,2);

cps(i,j)=serfe n(cpl,cp2,fctcl￡S'precisi on); end

[cpl2 n]=aa(cpa_min,cp_max/fd/cl/cs,mode); for j=l:n

cpl=cpl2(j,l);

cp2=cpl2(j,2);

cpa(ij)=aerfe n(cpbcp2,fd￡l‘cs,precisi on); end

end

h=zeros(mode,2);

%相速度

figure(l)

for j=l:2

ifj==l

cp=cps;

color=,b,;

else

cp=cpa;

color二T;

end

for i=l:mode

cpp=cp(:,i);

in d=fi nd(cpp==0);

if ^isempty(ind)

h(i/j)=plot((fdO(ind(end)+l:end))/d_Q235/cpp(ind(end)+l:end),color);

else

h(i/j)=plot(fdO/d_Q235,cpp/color);

end

hold on

end

ifj==2

xlabel('f/(KHz)')

ylabel('C_{p}/(km-sA{-l})')

title('6mm钢板相速度频散曲线J

set(gca,,xtick,,(0::20)*le3/d_Q235,,xticklaber/(0::20)*le3/d_Q235)

xlim([0, 1000]);%

set(gca,'ylim,,[0 cp_max],'ytick,,(0:cp_max/le3)*le3,...

'yticklaber,0:cp_max/le3)

grid on

hSGroup = hggroup;%要在子对象构建之后构建，构建后立即使用，否则将失效

hAGroup = hggroup;

set(h(:/l)/parent,,hSGroup)

set(h(:,2)/parent,/hAGroup)

set(get(get(hSGroup,'A nn otationfLege ndlnformation'),..?

'Icon DisplayStyle'/on1);

set(get(get(hAGroup,'A nn otationfLege ndlnformatiorV),..?

'Icon DisplayStyle'/on1);

legend)'对称模式T反对称模式')

end

end

%群速度

figure(2)

for j=l:2

ifj==l

cp=cps;

color=,b,;

else

cp=cpa;

color=,r,;

end

for i=l:mode

cpp=cp(:,i);

in d=fi nd(cpp==0);

if ^isempty(ind)

fd=fdO(ind(end)+l:end)';

cpp=cpp(i nd(end)+l:end);

else

fd=fdO';

end

dcdf=diff(cpp)/dfd;

cg=cpp(l:end-l).A2./(cpp(l:end-l)-fd(l:end-l).*dcdf);

h(i,j)=plot(fd(l:end-l)/d_Q235,cg/color);

hold on

end

ifj==2

xlabel('f/(KHz)')

ylabel('C_{g}/(km-sA{-l})')

title('6mm钢板群速度频散曲线，)

set(gca,,xtick,,(0::20)*le3/d_Q235,,xticklaber/(0::20)*le3/d_Q235) xlim([0,1000]);%

set(gca；ylim,/[0]*le3；ytick,/(0::*le3,,yticklabel I,0::

grid on

hSGroup = hggroup;%要在子对象构建之后构建，构建后立即使用，否则将失效hAGroup

= hggroup;

set(h(:/l)/parent'/hSGroup)

set(h(:,2)/parent,,hAGroup)

set(get(get(hSGroup,'A nn otationfLege ndlnformatiorf),…

'Icon DisplayStyle'/on1);

set(get(get(hAGroup,'A nn otationJ/Lege ndlnformatiorV),..?

Ucon DisplayStyle'/on1);

legend)'对称模式,，反对称模式')

end

end

%toc

end

%对称模式

function [cpO n]=ss(cp_min/cp_max,fd,cl,cs/mode)

cp2=cp_mi n;

deter=33;

cp0=zeros(mode,2);

n=0;

while cp2

cpl=cp2;

cp2=cpl+deter;

yl=smode(cpl,fd/cl/cs);

y2=smode(cp2,fd/cl/cs);

while yl*y2>0&&cp2

cpl=cp2;

cp2=cpl+deter;

yl=smode(cpl/fd/cl,cs);

y2=smode(cp2/fd/cl,cs);

end

讦yl*y2<0

n=n+l;

cp0(n/:)=[cpl cp2];

else

if yl==0&&y2^=0

n=n+1;

A 7d <2*d )s 8*

尸

d

d :7p 4±d )￡f F d %H S 4

① s o 二 b

晏 d

y (尸 c b 、2

*

d

l

s 花

R V d r (尸d <

E *d )s o :J *7

d

<

p

r

_

d )

c

u

_

$7

d

%d 42 尿w 眠 b

?粢

俅I?d x p y d o oa oa s o c d

。奇S a

A 7c b 、2*d )-C s o ^(

d 晏 d

d

r u 30v (l ;R V ?-+

(

d

?d <2*

d )JZ s o :>*7c b 、p 4*d )

c

u

0F d %6s 4 尿MlHItb

足 d f s f v d

。七

三uv(p

、d￡￥bs)sqe=:b

、((T R V (s

)s q f d (svqprcbopolusHSJ

u o l o u n j

pu

① P

U CD

P U CD

p u o

H+zdoHecb

,【

z d o z d o

H c c o d o "4-UHU

77d 0 T c b H r c

o d o WUHU

O

H H Z A oa od O H

上

A

4O S O

XH +e d f e d 。

CN

、(

rH d ￥e

d o N c b pue

P U CD

>l e

CN d o H H d

o

O H n ^>

4-

S

-(D

>l e

(D 」q

“d o H e d o OHHTA

七

CD S -CD

>l e

(D

」

q

0n H 0>

4O S 6

o d o H T d

end

%反对称模式

function [cpO n]=aa(cp_min,cp_max/fd z cl,cs,mode)

cp2二cp_mi n;

deter=33;

cp0=zeros(mode,2);

n=0;

while cp2

cpl=cp2;

cp2=cpl+deter;

yl=amode(cpl/fd,cl/cs);

y2=amode(cp2/fd,cl/cs);

while yl*y2>0&&cp2

cpl=cp2;

cp2=cpl+deter;

yl=amode(cpl/fd,cl/cs);

y2=amode(cp2/fd/cl/cs);

end

讦yl*y2<0

n 二n+1;

cpO(n/:)=[cpl cp2];

else

if yl==O&&y2^=O

n=n+1;

cpO(n,:)=[cpl cpl];

elseif y2==O&&yl~=O

n=n+1;

cp0(n/)=[cp2 cp2];

cp2=cp2+l;

elseif yl==O&&y2==O

n=n+l;

cpO(n,:)=[cpl cpl];

n=n+1;

cp0(n,:)=[cp2 cp2];

cp2=cp2+l;

end

end

end

end

function fs二amode(cp,fd￡l,cs)

p=abs(sqrt((cp/cs)A2-l));

q=abs(sqrt((cp/cl)A2-l));

if cp<=cs%p和q都是复数

fs=-4*p*q*sinh(pi*fd/cp*p)*cosh(pi*fd/cp*q)+(-p A2-l)A2*sinh(pi*fd/cp*q)*cosh(pi*fd/cp*p); elseif cp>cs&&cp<=cl%p 是实数，q 是复数

fs=4*p*q*sin(pi*fd/cp*p)*cosh(pi*fd/cp*q)+(p A2-l)A2*sinh(pi*fd/cp*q)*cos(pi*fd/cp*p);