该楼层疑似违规已被系统折叠 隐藏此楼查看此楼

close all; clear all;clc

a=1000;           % a为水的密度

f=1.0e6;          %   频率

c1=1.5e3;          %    水中传播速度

c2=6.4e3;          %    在铝中的传播速度

a1=0;              %   内环半径

a2=25e-3;            %   外环半径

u=5e-3;             %    平面振子的振速

r=146e-3;              %    凹面的曲率半径

k1=2*pi*f/c1;        %    在介质水中的波矢

k2=2*pi*f/c2;        %    在介质铝中的波矢

r11=sqrt(2*r*(r-sqrt(r^2-a1^2)));    %   内弦长(根据图求出

r12=sqrt(2*r*(r-sqrt(r^2-a2^2)));    %   外弦长

d1=1e-2;                      %    坐标中的x步长

d2=1e-2;                      %     坐标中的y的步长

d3=1e-3;                      %     积分中的r1步长

d4=(2*pi)/100;                %     积分中的角度z的步长

x1=0.12;        %     坐标中的x下限

x2=0.25;        %     坐标中的x上限

y1=-0.05;       %     坐标中的y下限

y2=0.05;        %     坐标中的y上限

m=0;                    %    位置x的计数

for x=x1:d1:x2;     %   第一层循环,坐标x. x代表列，注意单位为米

x;

m=m+1;

n=0;                %     坐标y的计数，注意n每次运行前要清零

for y=y1:d2:y2;     %   第二层循环,坐标y. y代表行，注意单位为米

y;

n=n+1;           %   结果计数

p=0;       %    声压的初值

is=0;       %    r1的下标

for r1=r11:d3:r12      %   第一个积分变量积分,r11,r12代表原公式中的r

r1;

is=is+1;

j=0;     %   角度z的计数

for z=0:d4:2*pi       %   第二个积分变量积分

z;

j=j+1;    %   结果计数

p=p+(k1*c1*a*u/2*pi)*((x^2+y^2-2*y*r1*(1-(r1^2/(4*r^2)))^0.5*cos(z)+(1-x/r)*r1^2)^(-0.5)*exp(-i*(r1^2/(2*r)*k2+k1*(x^2+y^2-2*y*r1*(1-(r1^2/(4*r^2)))^0.5*cos(z)+(1-x/r)*r1^2)^0.5))*r1*d3*d4);     %二重积分

end

end

v(m,n)=p   %    取声压p的幅值

end

end

last=abs(v);

col=x1:d1:x2;

ln=y1:d2:y2;

vmax=max(max(last));

v=sqrt(2)/2*vmax;

surf(ln,col,last);     %绘图

figure;

contourf(ln,col,last,[v v]);

grid;

这是圆形换能器的声场模拟程序 求矩形的！！！！