Hello all,
I'm trying to understand how works K-wave toolbox.
I'm trying to simulate the propagation of acoustic signal in water and I would like understand well the attenuation of acoustic signal in k-Wave simulations, because next step is simulate non linear parametric armonics in water.
With this aim I simulate a easy simulation avoiding the alfa coeff and alfa power and observing how its applied the radial divergence of the acoustic field (that I expect something related to P0/distance).
I found a little surprise when I recover the data at sensor points and if I plot pressure over distance for sensor_data.p_max and theoretic P0/distance seems that are not in good agreement (probably because I miss something at simulation code...).
Any help would be very appreciated.
Here you can find a link for plots
https://drive.google.com/open?id=0B2P-tB_iDXRidDN3SEJoV3IwOEU
https://drive.google.com/open?id=0B2P-tB_iDXRiNGNDV1BwcEQ3d2c
And here is the code used:
________________________________________________________________________________________
close all
clear all
clc
%Mallado del medio.
f_max = 50000; % Frecuencia de emision del plano radiante. [Hz]
c_max = 1500; % Velocidad maxima del medio [m/s]
long_onda = c_max/f_max; %Longitud de onda de la señal emitida
step_min = long_onda/6; % Precision minima requerida para simulacion lineal en el mallado
dx = 2e-03; % Paso espacial para realizar el mallado debe ser menos a step_min
dy = 2e-03;
% if dx > step_min
% disp('El mallado espacial debe poseer un paso menor.')
% break
% end
% if dx <= step_min
% disp('El mallado espacial es menor a long onda/6')
% end
x = 0.5; % Dimensiones de la malla en la direccion x [m]
y = 2.2; % Dimensiones de la malla en la direccion y [m]
Nx = round(x/dx); % Elementos de la malla direccion X
Ny = round(y/dy); % Elementos de la malla direccion X
kgrid = makeGrid(Nx, dx, Ny, dy);% Se define la mallado donde se realizara la simulacion.
[t_array, dt1] = makeTime(kgrid, c_max);
CFL = (c_max*dt1)/(dx);
CFL = 0.6;
t_end = y/c_max;
[kgrid.t_array, dt] = makeTime(kgrid, c_max, CFL, (t_end));
%Propiedades de los medios.
medium.sound_speed = 1500*ones(Nx, Ny); % [m/s] Velocidad del sonido en agua.
medium.density = 1000*ones(Nx, Ny); % [kg/m^3] Densidad del agua
%Parametrizacion de la atenuacion del medio (El valor del aire no esta bien definido)
% att1 = 2.17e-3; %Attenuation dB/cm·MHz2
% medium.alpha_coeff=att1*ones(Nx, Ny);
% medium.alpha_power=2;
%Propiedades del sistema emisor 2D
radio_dist = 0.0375;%Radio del disco emisor.
r = radio_dist/dx;
em_pos_y = round(0.1/dx); %El emisor se coloca en la posicion 0.5m
source.p_mask = zeros(Nx, Ny);% Se genera el espacio 2D donde colocar la linea 2D
source.p_mask(round((Nx/2)-r):round((Nx/2)+r),em_pos_y) = 1;
f_signal = 50000;
Mag_Pa = 15;
signal = Mag_Pa.*sin(2*pi*f_signal.*t_array);
%signal = signal(1:1000);
source.p = signal;
%Se genera el array de sensores.
rec_pos_y = round(0.2/dx):round(0.1/dx):round(2/dx);
rec_pos_x = round(Nx/2);
sensor.mask = zeros(Nx, Ny);
sensor.mask(rec_pos_x,rec_pos_y) = 1;
% imagesc(sensor.mask); colorbar %Comprobar la
sensor.record = {'p', 'p_max'};
input_args = { 'DisplayMask', source.p_mask|sensor.mask, 'RecordMovie', true, 'MovieType', 'image', 'MovieName', 'Atenuacion_50kHz_Prueba_0.2m_2m_01m_NOCOEFF', 'PMLInside', false};
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor,input_args{:});
D_P = [dx.*rec_pos_y' sensor_data.p_max];
% Atenuacion comprobacion
save('Atenuacion_50kHz_Prueba_0.2m_2m_01m_NOCOEFF.txt','D_P','-ascii')