matlab 蜂窝板STL计算问题

王良一

根据这篇论文，用matlab写的STL计算的程序，但是出的图和论文里的不一样，参数对了很多遍应该没有什么问题，这是写的程序的问题吗。。。。。

a=0.1524; %x边长
Ef=71e9;  %铝的弹性模量
b=0.136;  %y边长
h=0.006;  %芯层厚度
d=0.001;  %面板厚度
Vf=0.33;  %泊松比
t=0.000866;  %蜂窝壁厚
l=0.008;   %蜂窝边长
Gf=2.6692e10;  %铝剪切弹性模量
Ge=1*t*Gf/(3^0.5*l);  %芯层等效剪切弹性模量
C=Ge*h/100;   %芯层剪切刚度
C0=340;    %声速
r0=1.293;  %空气密度
rf=2810;   %面板密度

Q=0; %入射角
Z=90;    %θ
re=8*t*rf/(3*3^0.5*l);  %芯层等效密度
r=h*re+2*d*rf;  %夹层板单位面积质量
D=Ef*(h+d)^2*d/(2*(1-Vf^2));  %夹层板抗弯刚度

K=pi/a*(D/C)^0.5;   %无量纲参数

Emn11 = zeros(1,1000);
for W = 1:1000
    w = (W+1)*10;
    Emn = 0;
    Bmn = 0;
    eSum = 0;
    iSum = 0;
   
    K0=w/C0;    %空气中声波数量
    Kx=K0*sin(Q)*cos(Z);  %波数在x分量
    Ky=K0*sin(Q)*sin(Z);  %波数在y分量
    Kz=K0*cos(Q);   %波数在z分量
   
    for m = 1:70
        for n = 1:70
            
            Itop1 = 4*m*n*pi^2;
            Itop2 = (exp(-(1j*(Kx*a))))*cos(m*pi)-1;
            Itop3 = (exp(-(1j*(Ky*b))))*cos(n*pi)-1;
            Ibot1 = (Kx*a)^2-(m*pi)^2;
            Ibot2 = (Ky*b)^2-(n*pi)^2;
            Imn1 = (Itop1*Itop2*Itop3)/(Ibot1*Ibot2);   %Imn公式
            
            Bmn1 = (m*pi/a)^2+(n*pi/b)^2;
            Etop1 = 2j*r0*(w^2);
            Etop2 = (D/C)*Bmn1 + 1;
            Etop3Temp = Kz*(h+2*d);
            Etop3 = exp(1j*(Etop3Temp));
            Ebot1 = Kz;
            Ebot2 = D*(Bmn1^2);
            Ebot3 = w^2*(r-((2j*r0)/Kz))*(((D/C)*Bmn1)+1);
            
            Emn1Temp = (Etop1*Etop2*Etop3)/(Ebot1*(Ebot2-Ebot3));
            Emn1 = Emn1Temp*Imn1;   %Emn公式
            
            Itemp = (real(Imn1))^2 + (imag(Imn1))^2;
            Etemp = (real(Emn1))^2 + (imag(Emn1))^2;
            eSum=eSum+Etemp;
            iSum=iSum+Itemp;
        end
    end

    temp = iSum/eSum;
    stemp =10*log10(temp);   %STL公式

    Emn11(1,W) = stemp;
end
xp = linspace(10,10000,1000);
plot(xp,Emn11);