Hello All
I build a cylindrical transducer in a water media with alpha being 1.5 coefficient being 0.75 and sound speed being 1540m/s as the source, and I made a sin wave with magnitude being 2 and frequency being 1.4MHz
It is a 3D simulation, I set PML size to be 2 and they are outside the computational domain.dx=dy=dz=0.25mm
For the first simulation, the computational domain size is 64(x)*64(y)*64(z) in grid. The transducer size is 64 in height(At Z axis) and 10 in radius (both in grid), with x and y coordinates being 32. For the second simulation, the computational domain size is 128*128*80 in grid and the transducer size is 80 in height and 40 in radius with x and y coordinates being 64
The sensor mask is defined as a plane at Z=32 to observe the overall distribution at the middle of the transducer.
However, some really interesting thing happened:
In the first simulation, the distribution is symmetrical in x-y plane
In the second simulation, it becomes asymmetrical and the result is very strange, one half is generally higher than the other half in the plane.
Here is the source code for the second simulation, one can change the size and center of the circle to change it into the first one. Can anybody help me to explain why this happen?
lear all;
PML_X_SIZE=2;
PML_Y_SIZE=2;
PML_Z_SIZE=2;
% The size in 3 dimension in grid
Nx=132-2*PML_X_SIZE;
Ny=132-2*PML_Y_SIZE; %
Nz=84-2*PML_Z_SIZE;
%Unit in each dimension
dx=0.25e-3; %[m],dy dz follow the same
dy=0.25e-3;
dz=0.25e-3;
%Build up the computational domain
kgrid=makeGrid(Nx,dx,Ny,dy,Nz,dz);
%Make the circle profile of the cylindrical transducer
circle=makeCircle(Nx,Ny,64,64,40); %center 1.6 cm
disk=makeDisc(Nx,Ny,64,64,40);
;
%Make the mask for the source and the sensor
mask=zeros(Nx,Ny,Nz);
mask_2=zeros(Nx, Ny, Nz);
for l=1:Nz
mask(:,:,l)=circle;
end
for l2=1:Nx
for m2=1:Ny
mask_2(l2,m2,40)=1; %This make a binary mask for the sensor array
end
end
%Visualization of the source(transducer)'s profile
voxelPlot(mask);
view([50,15]);
voxelPlot(mask_2);
view([50,15]);
%Define the media properties
medium.sound_speed = 1540; % [m/s]
medium.alpha_power = 1.5; % [dB/(MHz^y cm)]
medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)]
medium.density=1000; % [Kg/m^3]
%Define the simulation time length
t_end = 40e-6; % [s]
[kgrid.t_array,dt] = makeTime(kgrid, medium.sound_speed, [], t_end);
%Define the source
source.p_mask=mask;
source_freq = 1.5e6; % [Hz]
source_mag = 2; % [Pa]
source.p = source_mag*sin(2*pi*source_freq*kgrid.t_array); %The square wave for the source
source.p = filterTimeSeries(kgrid, medium, source.p); %Filter the frequency not supported by the grid
%Define the sensor array to record the data
sensor.mask=mask_2;
sensor.record={'p_rms','p_final','I_avg'}; % Here two types of data are recorded: root mean square of the pressure and average intensity over the simulation time
display_mask = source.p_mask;
input_args = {'DisplayMask', display_mask, 'PMLInside', false, 'PlotPML', false,'RecordMovie', true, 'MovieName','Cylindrical'}; %Set the arguments of the display
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:}); %Run the simulation function of 3D simulation
%Visualization of the Simulation Result
%Reorder the average acoustic intensity recorded here. The data returns as the
%Matlab column wise order
intensity_1=sensor_data.Iy_avg;
intensity_2=sensor_data.Ix_avg;
intensity_3=sensor_data.Iz_avg;
Map_intensity_1=zeros(Nx,Ny);
Map_intensity_2=zeros(Nx,Ny);
Map_intensity_3=zeros(Nx,Ny);
Map_intensity_4=zeros(Nx,Ny);
for p=1:Nx
for q=1:Ny
if(disk(q,p)==1)
continue;
end;
Map_intensity_1(q,p)=abs(intensity_1((p-1)*64+q,1));
Map_intensity_2(q,p)=abs(intensity_2((p-1)*64+q,1));
Map_intensity_3(q,p)=abs(intensity_3((p-1)*64+q,1));
Map_intensity_4(q,p)=( Map_intensity_1(q,p)^2+ Map_intensity_2(q,p)^2+Map_intensity_3(q,p)^2)^0.5;
end;
end;
figure;
surf(kgrid.x_vec*1000,kgrid.y_vec*1000+64,Map_intensity_1);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
title('Average Acoustic Intensity At Y direction');
figure;
surf(kgrid.x_vec*1000,kgrid.y_vec*1000,Map_intensity_2);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
title('Average Acoustic Intensity At X direction');
figure;
surf(kgrid.x_vec*1000,kgrid.y_vec*1000,Map_intensity_3);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
title('Average Acoustic Intensity At Z direction');
figure;
surf(kgrid.x_vec*1000,kgrid.y_vec*1000,Map_intensity_4);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
title('Average Acoustic Intensity ');