www.pudn.com > pinsanquxian.rar > zzS01.m, change:2004-12-19,size:3035b

```%求解Lame波对称部分的基本模态So
format long
syms c;   %定义速度c为符号函数
xx=[];err=[];      %xx存贮不同KH值的波速，err存贮误差
eps=1e-12;         %精度
u=0.3; %泊松比u，越大对应速度越小
E=2.06e11; %弹性模量E,越大对应速度越大
p=7850;%密度越大对应速度越小
cd=((1-u)*E/(1+u)/(1-2*u)/p)^0.5;%cd为纵波波速，cs为横波波速
cs=(E/2/(1+u)/p)^0.5;
for kh=0:0.05:8       %k为波数，h为总的板厚
if kh<7.65          %利用lamea0作图法先估计出c值范围
y1=4*tanh(kh/2*(1-(c/cd)^2)^(1/2))*(1-(c/cd)^2)^(1/2)*((c/cs)^2-1)^(1/2);
y2=(2-(c/cs)^2)^2*tan(kh/2*((c/cs)^2-1)^(1/2));
if kh<0.05
c=[5365;5375];
elseif kh<2.55  c=[cs;5400];
elseif kh<3  c=[cs;xx(fix(20*kh))];
else   c=[cs;xx(fix(20*kh))];
end
for i=2:2
a=c(i-1);b=c(i);
c=(a+b)/2;
while (subs(y1-y2)~= 0) & ((b-a)>eps)
c=(a+b)/2;
fa=4*tanh(kh/2*(1-(a/cd)^2)^(1/2))*(1-(a/cd)^2)^(1/2)*((a/cs)^2-1)^(1/2)-(2-(a/cs)^2)^2*tan(kh/2*((a/cs)^2-1)^(1/2));
fb=4*tanh(kh/2*(1-(b/cd)^2)^(1/2))*(1-(b/cd)^2)^(1/2)*((b/cs)^2-1)^(1/2)-(2-(b/cs)^2)^2*tan(kh/2*((b/cs)^2-1)^(1/2));
fy=4*tanh(kh/2*(1-(c/cd)^2)^(1/2))*(1-(c/cd)^2)^(1/2)*((c/cs)^2-1)^(1/2)-(2-(c/cs)^2)^2*tan(kh/2*((c/cs)^2-1)^(1/2));
if fa*fy < 0
b=c;
else
a=c;
end
end
xx=[xx;c];
m=subs(y1-y2);
err=[err;m];
end
else
y1=4*tanh(kh/2*(1-(c/cd)^2)^(1/2))*(1-(c/cd)^2)^(1/2)*(1-(c/cs)^2)^(1/2);
y2=(2-(c/cs)^2)^2*tanh(kh/2*(1-(c/cs)^2)^(1/2));
c=[0;cs];
for i=2:2     %利用二分法求根
a=c(i-1);b=c(i);
c=(a+b)/2;
while (subs(y1-y2)~= 0) & ((b-a)>eps)
c=(a+b)/2;
fa=4*tanh(kh/2*(1-(a/cd)^2)^(1/2))*(1-(a/cd)^2)^(1/2)*(1-(a/cs)^2)^(1/2)-(2-(a/cs)^2)^2*tanh(kh/2*(1-(a/cs)^2)^(1/2));
fb=4*tanh(kh/2*(1-(b/cd)^2)^(1/2))*(1-(b/cd)^2)^(1/2)*(1-(b/cs)^2)^(1/2)-(2-(b/cs)^2)^2*tanh(kh/2*(1-(b/cs)^2)^(1/2));
fy=4*tanh(kh/2*(1-(c/cd)^2)^(1/2))*(1-(c/cd)^2)^(1/2)*(1-(c/cs)^2)^(1/2)-(2-(c/cs)^2)^2*tanh(kh/2*(1-(c/cs)^2)^(1/2));
if fa*fy < 0
b=c;
else
a=c;
end
end
xx=[xx;c];
m=subs(y1-y2);
err=[err;m];
end
end
end
cp=xx;        %计算So群速度
kh=[0:0.05:8];
for i=2:1:161
if cp(i-1) ~=cp(i)
cg(i)=cp(i)+kh(i)*(cp(i)-cp(i-1))/(kh(i)-kh(i-1)); %波的相速度cp与群速度cg的转换公式
else cg(i)=cg(i-1)-50;
end
end
cg(1)=0;
%作图
fh=[];
for kh=0:0.05:8
f=kh*cp(fix(20*kh)+1)/2/pi; %将kd与cp、cg的关系转换为fh与cp、cg的关系
fh=[fh;f];
end
figure(1)     %f与cp的图
plot(fh/2,cp);%板厚为2mm
xlabel('频率f (KHz)');ylabel('相速度Cp (m/s)');
grid on;
hold on;
figure(2)     %f与cg的图
plot(fh/2,cg);%板厚为2mm
xlabel('频率f (KHz)');ylabel('群速度Cg (m/s)');
grid on;
hold on;```