【物理应用】超声场可视化仿真模拟【含GUI Matlab源码 1494期】

⛄一、简介（附论文）

通过对超声场理论的数学物理方法计算，分别对圆型和矩型换能器的声轴线上声压分布、轴方向横截面的声压的分布及声场的指向性的表达式作出推导和演算，并得出结论；以及研究脉冲波声场分布特性，数值计算其声压分布，再利用数学软件Matlab进行可视化模拟进行的仿真研究，以活塞探头的各个项的参数对超声场的分布影响看作为研究的内容。最后通过Graphical User Interface设计出一个图形操作界面，有助于我们通过换能器的调参过程对声场的分布影响进行研究，同时也有助于我们更清楚地了解声场的分布理论，从而提升我们的研究效率。

⛄二、部分源代码

function varargout = gui_sound_field(varargin)
% GUI_SOUND_FIELD MATLAB code for gui_sound_field.fig
% GUI_SOUND_FIELD, by itself, creates a new GUI_SOUND_FIELD or raises the existing
% singleton*.
%
% H = GUI_SOUND_FIELD returns the handle to a new GUI_SOUND_FIELD or the handle to
% the existing singleton*.
%
% GUI_SOUND_FIELD(‘CALLBACK’,hObject,eventData,handles,…) calls the local
% function named CALLBACK in GUI_SOUND_FIELD.M with the given input arguments.
%
% GUI_SOUND_FIELD(‘Property’,‘Value’,…) creates a new GUI_SOUND_FIELD or raises the
% existing singleton*. Starting from the left, property value pairs are
% applied to the GUI before gui_sound_field_OpeningFcn gets called. An
% unrecognized property name or invalid value makes property application
% stop. All inputs are passed to gui_sound_field_OpeningFcn via varargin.
%
% *See GUI Options on GUIDE’s Tools menu. Choose “GUI allows only one
% instance to run (singleton)”.
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help gui_sound_field

% Last Modified by GUIDE v2.5 27-Oct-2021 14:38:02

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct(‘gui_Name’, mfilename, …
‘gui_Singleton’, gui_Singleton, …
‘gui_OpeningFcn’, @gui_sound_field_OpeningFcn, …
‘gui_OutputFcn’, @gui_sound_field_OutputFcn, …
‘gui_LayoutFcn’, [] , …
‘gui_Callback’, []);
if nargin && ischar(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% — Executes just before gui_sound_field is made visible.
function gui_sound_field_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% varargin command line arguments to gui_sound_field (see VARARGIN)
% Choose default command line output for gui_sound_field
set(handles.quadrate_checkbox,‘value’,0);
set(handles.circular_checkbox,‘value’,1);

%background_image1 = importdata(‘background.jpg’);
%axes(handles.background_axes);
%image(background_image1);
%alpha(0.5)
%axis off
sign1 = imread(‘tubiao.png’);
axes(handles.sign_axes);
image(sign1);
axis off
plot_image=importdata(‘plot.png’);
set(handles.plot_pushbutton,‘CDATA’,plot_image)
play_image=importdata(‘stop2.png’);
set(handles.play_pushbutton,‘CDATA’,play_image)
set(handles.text23,‘String’,[‘dB/m’]);
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes gui_sound_field wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% — Outputs from this function are returned to the command line.
function varargout = gui_sound_field_OutputFcn(hObject, eventdata, handles)
% varargout cell array for returning output args (see VARARGOUT);
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% — Executes on button press in plot_pushbutton.
function plot_pushbutton_Callback(hObject, eventdata, handles)
% hObject handle to plot_pushbutton (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

quadrate_flag=get(handles.quadrate_checkbox,‘value’);
circular_flag=get(handles.circular_checkbox,‘value’);
z0_double=str2num(get(handles.z0_edit,‘string’));
f0=str2num(get(handles.frequency_edit,‘string’)); % Transducer center frequency [MHz]
f0=f0*1e6;
c=str2num(get(handles.c_edit,‘string’)); % Speed of sound [m/s]( 水 1500 或者钢 5900)
decay=str2num(get(handles.decay_edit,‘string’)); %衰减系数
%------判断输入参数的正误--------%
if f0<20000
errorcall= errordlg( ‘The Ultrasonic Frequency is Higher than 20000Hz at least. Please input the parameter again.’ , ‘Error’ ) ;
dbcont;
end
if c<300
errorcall= errordlg( ‘The Speed of Sound is too small. Please input the parameter again.’ , ‘Error’ ) ;
dbcont;
end
if decay<0
errorcall= errordlg( ‘The Attenuation of Sound may not be less than zero. Please input the parameter again.’ , ‘Error’ ) ;
dbcont;
end
%------判断输入参数的正误--------%
%-----clear axes-----%
axes(handles.acoustic_axis_axes);
acoustic_axis_delete=get(gca,‘children’);
delete(acoustic_axis_delete);
axes(handles.sound_field_axes);
sound_field_delete=get(gca,‘children’);
delete(sound_field_delete);
axes(handles.directivity_axes);
directivity_delete=get(gca,‘children’);
delete(directivity_delete);
%-------------------%

pausetime=0.05;
flag=1;
p0=101.32510^3; %大气压强101.32510^3pa
Rp=1;
lambda=c/f0; % Wavelength
k=2pi/lambda; %波数
w=2pif0; %角频率
if quadrate_flag0 && circular_flag1
%---------------------------------------圆形活塞----------------------------------------------%
%圆形活塞参数
Rs=str2num(get(handles.radius_edit,‘string’)); %活塞半径25mm
Rs=Rs10^-3;
Fs=piRs^2; %活塞的面积
%各种材料的声阻抗
Z_air=3401.29;
Z_water=150010^3;
Z_steel=str2num(get(handles.Z_R_edit,‘string’));
%工件参数
z_d=str2num(get(handles.z_d_edit,‘string’));
z_d=z_d10^-3; %厚度
%定义网格，数值求解
n=150;
x=linspace(-1.5Rs,1.5Rs,n);
y=linspace(-1.5Rs,1.5Rs,n);
%轴向初始距离
near_field=(2Rs)^2/(4lambda); %近场距离
z0=z0_doublenear_field;
%------判断输入参数的正误--------%
if Rs<=0
errorcall= errordlg( ‘The Size of Instrument may not be less than zero. Please input the parameter again.’ , ‘Error’ ) ;
dbcont;
end
if z_d<=0
errorcall= errordlg( ‘Workpiece Size may not be less than zero. Please input the parameter again.’ , ‘Error’ ) ;
dbcont;
end
if z0_double<0 || z0>z_d
errorcall= errordlg( ‘The Initial Position may not be less than zero or more than Workpiece Size. Please input the parameter again.’ , ‘Error’ ) ;
dbcont;
end
%------判断输入参数的正误--------%
%反馈近场距离到界面
set(handles.near_field_text,‘String’,num2str(near_field10^3));
%--------------------------%
z=z0;
detaz=z;
[x,y]=meshgrid(x,y);
r=sqrt(x.^2+y.2+z.^2);
Sita=acos(z./r);
Phi=asin(y./(r.sin(Sita)));
R_acoustic_axis_incident=linspace(0,z_d,10000);
R_acoustic_axis_reflect=linspace(z_d,10(2Rs)^2/(4lambda),10000);
%-----轴线上声压-----%
P_field_incident=P_circular_acoustic_axis(lambda,Rs,R_acoustic_axis_incident).exp(-decayR_acoustic_axis_incident/8.68);
P_field_reflect=rp(Z_steel,Z_air)P_circular_acoustic_axis(lambda,Rs,R_acoustic_axis_reflect).exp(-decayR_acoustic_axis_reflect/8.68);
%由结论，直接代入参数
P=RpSound_pressure_circular( k,Rs,Sita, w,p0,c,r).exp(-decayr/8.68);
P=abs§;
%指向性
[theta,phi]=meshgrid(linspace(0,2pi,3n),linspace(0,pi/2,3n));
X=kRssin(phi);
J1=besselj(1,X);
D=abs(2J1./X);
%坐标变换
%[Dx,Dy,Dz]=[D.*((sin(Phi)).cos(Sita)),D.((sin(Sita)).sin(Phi)),D.((cos(Phi)).cos(Sita-Sita))];
Dx=D.((sin(phi)).cos(theta));
Dy=D.((sin(theta)).sin(phi));
Dz=D.((cos(phi)).*cos(theta-theta));

axes(handles.acoustic_axis_axes);
plot(R_acoustic_axis_incident,abs(P_field_incident));
hold on
plot(R_acoustic_axis_reflect,abs(P_field_reflect),'r');
plot(linspace(z_d,z_d,10),linspace(0,2,10),'--ks','LineWidth',1,'MarkerSize',2);
text(z_d,abs(P_field_incident(1000)),['reflect point']);
grid on;
title(['Circular pressure transducer acoustic axis of distribution'],'FontSize',7);
xlabel(['z/m']);
ylabel(['P/P0']);

axes(handles.sound_field_axes);
p=surf(x,y,P+0.7*((1-flag)*z_d+flag*z)*10^8);
set(p,'facealpha',z0/z-0.1);
material shiny
shading interp
colormap(cool)
light ('position',[-1 -0.5 2],'style','infinite')
hold on
cube=surf_cube( 4*Rs,4*Rs,0.7*z_d*10^8,-2*Rs,-2*Rs,0);
axis([-3*Rs 3*Rs -3*Rs 3*Rs 0 1.5*(0.7*z_d*10^8)]);
title(['When z=',num2str(((1-flag)*z_d+flag*z)*10^3),'mm,Circular pressure transducer acoustic axis of the three-dimensional distribution'],'FontSize',7);
text(0,0,0.04*z_d*10^8,['Workpiece:Plane of incidence'],'color','r');
text(0,0,0.7*z_d*10^8,['Reflective surface'],'color','r');

axes(handles.directivity_axes);
D_surf=surf(Dx,Dy,Dz);
set(D_surf,'facealpha',0.5);
axis([-0.2 0.2 -0.2 0.2 0 1])
title(['Beam directiviy'],'FontSize',7);
colorbar
material shiny
shading interp
colormap(cool)
%---------------------------------------圆形活塞----------------------------------------------%

elseif quadrate_flag1 && circular_flag0
%---------------------------------------矩形活塞----------------------------------------------%
%各种材料的声阻抗
Z_air=3401.29;
Z_water=150010^3;
Z_steel=str2num(get(handles.Z_R_edit,‘string’));

%矩形活塞换能器参数:边长a b，面积Fs
a=str2num(get(handles.length_edit,'string'))*10^-3;
b=str2num(get(handles.width_edit,'string'))*10^-3;
Fs=a*b;
%工件参数
z_d=str2num(get(handles.z_d_edit,'string'));
z_d=z_d*10^-3;                     %厚度
 %---------近场-------------%
near_field=Fs/(2*pi*lambda);
z0=z0_double*Fs/near_field;
%------判断输入参数的正误--------%
if a<=0 || b<=0
errorcall= errordlg( 'The Size of Instrument may not be less than zero. Please input the parameter again.' , 'Error'  ) ;
dbcont;
end
if z_d<=0
errorcall= errordlg( 'Workpiece Size may not be less than zero. Please input the parameter again.' , 'Error'  ) ;
dbcont;
end
if z0_double<0 || z0>z_d
errorcall= errordlg( 'The Initial Position may not be less than zero or more than Workpiece Size. Please input the parameter again.' , 'Error'  ) ;
dbcont;
end
%------判断输入参数的正误--------%
%空间网格化求数值解
n=120;                            %
x=linspace(-a,a,n);
y=linspace(-b,b,n);
Sita=linspace(0,pi/2,n);
apha=linspace(0,2*pi,n);
[Sita,apha]=meshgrid(Sita,apha);

X_acoustic_axis=linspace(-a,a,n);
Z_acoustic_axis_incident=linspace(0,z_d,n);
Z_acoustic_axis_reflect=linspace(z_d,1.5*z_d,n);
[X_acoustic_axis_incident,Z_acoustic_axis_incident]=meshgrid(X_acoustic_axis,Z_acoustic_axis_incident);
[X_acoustic_axis_reflect,Z_acoustic_axis_reflect]=meshgrid(X_acoustic_axis,Z_acoustic_axis_reflect);
R_acoustic_axis_incident=sqrt(X_acoustic_axis_incident.^2+Z_acoustic_axis_incident.^2);
R_acoustic_axis_reflect=sqrt(X_acoustic_axis_reflect.^2+Z_acoustic_axis_reflect.^2);
phi_acoustic_axis=0;
sita_acoustic_axis_incident=asin(X_acoustic_axis_incident./(R_acoustic_axis_incident.*cos(phi_acoustic_axis)));
sita_acoustic_axis_reflect=asin(X_acoustic_axis_reflect./(R_acoustic_axis_reflect.*cos(phi_acoustic_axis)));
%-------------------%
%-----------反馈近场距离到界面--------------%
set(handles.near_field_text,'String',num2str(near_field*10^3));
%------------------------------------------%
z=z0;
[x,y]=meshgrid(x,y);
r=sqrt(x.^2+y.^2+z.^2);
phi=asin(y./r);
sita=asin(x./(r.*cos(phi)));
A=k*a*cos(apha).*sin(Sita);
B=k*b*sin(apha).*sin(Sita);
%直接用结论求解
%轴线声压
P_acoustic_axis_incident=P_rectangle_acoustic_axis( Fs,a,k,lambda,R_acoustic_axis_incident ,sita_acoustic_axis_incident,phi_acoustic_axis);
P_acoustic_axis_incident=P_acoustic_axis_incident.*exp(-decay*R_acoustic_axis_incident/8.68);
P_acoustic_axis_reflect=P_rectangle_acoustic_axis( Fs,a,k,lambda,R_acoustic_axis_reflect ,sita_acoustic_axis_reflect,phi_acoustic_axis);
P_acoustic_axis_reflect=rp(Z_steel,Z_air)*P_acoustic_axis_reflect.*exp(-decay*R_acoustic_axis_reflect/8.68);
axes(handles.acoustic_axis_axes);
p_acoustic_axis_inciden=surf(X_acoustic_axis_incident,Z_acoustic_axis_incident,abs(P_acoustic_axis_incident));
hold on
p_acoustic_axis_reflect=surf(X_acoustic_axis_reflect,Z_acoustic_axis_reflect,abs(P_acoustic_axis_reflect));
plot3(linspace(0,0,10),linspace(z_d,z_d,10),linspace(0,20,10),'--ks','LineWidth',1,'MarkerSize',2);
text(0,z_d,15,['reflect point']);
title(['Recrangular pressure transducer acoustic axis of distribution'],'FontSize',7);
material shiny
shading interp
colormap(cool)
set(p_acoustic_axis_inciden,'facealpha',0.7);
set(p_acoustic_axis_reflect,'facealpha',0.7);
light ('position',[1 1 20],'style','infinite');
view([80,16]);
%声压
P=Sound_pressure_rectangle( p0,Fs,lambda,r,k,a,b,sita,phi).*exp(-decay*r/8.68);
P=abs(P);
%指向性
D_aph_Sita_omiga=(sin(A)./A).*(sin(B)./B);
D_aph_Sita_omiga=abs(D_aph_Sita_omiga);
axes(handles.sound_field_axes);
p=surf(x,y,P+0.7*((1-flag)*z_d+flag*z)*10^7);
%set(p,'facealpha',z0/z-0.2);
title(['When z=',num2str(((1-flag)*z_d+flag*z)*10^3),'mm,Recrangular pressure transducer acoustic axis of the three-dimensional distribution'],'FontSize',7);
material shiny
shading interp
colormap(cool)
light ('position',[-1 -0.5 2],'style','infinite')
hold on
cube=surf_cube( 4*a,4*b,0.7*z_d*10^7,-2*a,-2*b,0);
axis([-3*a 3*a -3*b 3*b 0 1.1*0.7*z_d*10^7]);
%坐标变换
Dx=D_aph_Sita_omiga.*(sin(Sita)).*cos(apha);
Dy=D_aph_Sita_omiga.*(sin(Sita)).*sin(apha);
Dz=D_aph_Sita_omiga.*(cos(Sita)).*cos(apha-apha);
text(0,0,0,['Workpiece:Plane of incidence'],'color','r');
text(0,0,0.7*z_d*10^7,['Reflective surface'],'color','r');

axes(handles.directivity_axes);
D_surf=surf(Dx,Dy,Dz);
set(D_surf,'facealpha',0.5)
axis([-0.2 0.2 -0.2 0.2 0 1])
title(['Beam directiviy'],'FontSize',7);
material shiny
shading interp
colormap(cool)
colorbar

end