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


%求解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<6.5        %利用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=[5333.9;5334]; 
     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;