k-Wave User Forum » Topic: Help on Attenuation !! Very Confused

k-Wave User Forum » Topic: Help on Attenuation !! Very Confused http://www.k-wave.org/forum/topic/help-on-attenuation-very-confused Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 02:31:31 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "Help on Attenuation !! Very Confused" http://www.k-wave.org/forum/topic/help-on-attenuation-very-confused#post-7855 Mon, 28 Sep 2020 19:09:03 +0000 Bradley Treeby 7855@http://www.k-wave.org/forum/ Hi RuiLiu, Regarding the first post, this is not related to the absorption properties or to k-Wave. Rather, when you model an initial value problem in 1D, the initial pressure distribution will split in two, with part of the wave propagating to the left and part travelling to the right. Each will have half the amplitude of the initial pressure distribution. You should find a more detailed explanation in most good textbooks on acoustics. Regarding the second post, the numerical model used in <code>pstdElastic2D</code> is not unconditionally stable. If the model is unstable with the chosen material parameters, you will need to reduce the size of the time step / CFL. Hope that helps, Brad RuiLiu on "Help on Attenuation !! Very Confused" http://www.k-wave.org/forum/topic/help-on-attenuation-very-confused#post-7852 Thu, 24 Sep 2020 21:49:06 +0000 RuiLiu 7852@http://www.k-wave.org/forum/ Hi, Brad. Thanks for the last reply. It help a lot. But I also confused about the attenuation is Elastic2d Model. I am trying to simulate the wave propagation in sandstone, and I set the compression alpha coefficient as 2.7[dB/MHZ^2 cm], 5.8[dB/MHZ^2 cm] for shear. But it seems the system is not stable. I don't understand how that happened. Thanks!! clc clear close all % ========================================================================= % DEFINE SIMULATION PROPERTIES % ========================================================================= % define the properties used in the simulation c0 = 3760; % compressional sound speed [m/s] shear_c0 = 2300; % shear sound speed [m/s] rho0 = 2650; % density [kg/m^3] frac_c0 = 1481; % open fracture compression sound speed [m/s], water fracshear_c0 = 1; % open fracture shear sound speed [m/s], water frac_rho0 = 1000; % open fracture density [kg/m^3], water points_per_wavelength = 30; % number of grid points per wavelength at f0 wavelength_separation = 10; % separation between the source and detector(farest) pml_size = 20; % PML size pml_alpha = 2; % PML alpha % compute corresponding grid spacing dx = 0.0001; % grid point spacing in the x direction [m] dy = dx; % grid point spacing in the y direction [m] % compute corresponding grid size Nx = wavelength_separation*points_per_wavelength + 20; % number of grid points in t he x (row) direction Ny = wavelength_separation*points_per_wavelength + 20; % number of grid points in the y (column) direction material_length = (wavelength_separation*points_per_wavelength*dx)*1e3; % [mm] grid_length = 3*dx*1e3;% [mm]fracture is three grids thick % ========================================================================= % SIMULATION % ========================================================================= % create the computational grid kgrid = makeGrid(Nx, dx, Ny, dy); % define a single sensor position an integer number of wavelengths away sensor.mask = zeros(Nx, Ny); y_loc = (10:60:wavelength_separation*points_per_wavelength + 10); x_loc = (10+60:60:wavelength_separation*points_per_wavelength + 10 - 60); sensor.mask(10, y_loc) = 1; sensor.mask(Nx-10, y_loc) = 1; sensor.mask(x_loc, 10) = 1; sensor.mask(x_loc, Ny-10) = 1; % define the crack locations start_x = randi([10 Nx-10],1,1); start_y = randi([10 Nx-10],1,1); dist = randi([50,200],1,1); angle_degree = randi([1 180],1,1); end_x = round(start_x + cos(angle_degree*3.1415926/180).* dist ); end_y = round(start_y + sin(angle_degree*3.1415926/180).* dist ); while end_x <=10 || end_x >= Nx-10 || end_y <= 10 || end_y >= Nx-10 start_x = randi([10 Nx-10],1,1); start_y = randi([10 Nx-10],1,1); dist = randi([30,200],1,1); angle_degree = randi([1 180],1,1); end_x = round(start_x + cos(angle_degree*3.1415926/180).* dist ); end_y = round(start_y + sin(angle_degree*3.1415926/180).* dist ); end [xx_new,yy_new] = XiaolinWu(start_x-1, start_y, end_x-1, end_y); [x_new,y_new] = XiaolinWu(start_x, start_y, end_x, end_y); [xxx_new,yyy_new] = XiaolinWu(start_x+1, start_y, end_x+1, end_y); % define the absorption properties medium.alpha_coeff_compression = 2.7; % [dB/(MHz^2 cm)]?????? medium.alpha_coeff_shear = 5.4; % [dB/(MHz^2 cm)]?????? % define the properties of the propagation medium medium.sound_speed_compression = c0*ones(Nx, Ny); % [m/s] medium.sound_speed_shear = shear_c0*ones(Nx, Ny); % [m/s] medium.density = rho0*ones(Nx, Ny); % [kg/m^3] for j = 1: length(x_new) medium.sound_speed_compression(x_new(j), y_new(j)) = frac_c0; % [m/s] medium.sound_speed_shear(x_new(j), y_new(j)) = fracshear_c0; % [m/s] medium.density(x_new(j), y_new(j)) = frac_rho0; medium.sound_speed_compression(xx_new(j), yy_new(j)) = frac_c0; % [m/s] medium.sound_speed_shear(xx_new(j), yy_new(j)) = fracshear_c0; % [m/s] medium.density(xx_new(j), yy_new(j)) = frac_rho0; medium.sound_speed_compression(xxx_new(j), yyy_new(j)) = frac_c0; % [m/s] medium.sound_speed_shear(xxx_new(j), yyy_new(j)) = fracshear_c0; % [m/s] medium.density(xxx_new(j), yyy_new(j)) = frac_rho0; end figure(1) imagesc(medium.density); hold on scatter (10, Nx/2, 'r.') hold on plot (y_loc, 10, 'k.') hold on plot (10, x_loc, 'k.') hold on plot (y_loc, 310, 'k.') hold on plot (310, x_loc, 'k.') % set end time t_end_cal = (material_length/3760)*10^6; t_end = 15e-6; % end at 50us % create the time array kgrid.t_array = makeTime(kgrid, max(medium.sound_speed_compression(:)) ,[], t_end); % define the initial pressure distribution source term disc_magnitude = 1e5; % [Pa]... disc_x_pos = Nx/2; % [grid points] disc_y_pos = 10; % [grid points] disc_radius = 1; % [grid points] source.p0 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius); source.p0 = smooth(kgrid, source.p0, true); % set the simulation options input_args = {'PlotFreq', 20, 'PlotScale', 'auto', 'PlotLayout', false,... 'PMLInside', false, 'PMLSize', pml_size, 'PMLAlpha', pml_alpha, 'PlotPML', true}; % define the sensor to record the maximum particle velocity everywhere sensor.record = {'p'}; % run the simulation sensor_data = pstdElastic2D(kgrid, medium, source, sensor, input_args{:}); % ========================================================================= % VISUALISATION % ========================================================================= %% [t_sc, scale, prefix] = scaleSI(max(kgrid.t_array(:))); figure(4) for n = 1 : 6 subplot(6, 1, n), plot(kgrid.t_array*scale, sensor_data.p(n,:), 'r-'); xlabel(['Time [' prefix 's]']); ylabel('Signal Amplitude-Pressure'); axis tight; % ylim([-max(sensor_data(n,:)) max(sensor_data(n,:))]) end %% figure(5) for n = 7 : 14 subplot(8, 1, n-6), plot(kgrid.t_array*scale, sensor_data.p(n,:), 'r-'); xlabel(['Time [' prefix 's]']); ylabel('Signal Amplitude-Pressure'); axis tight; % ylim([-max(sensor_data(n,:)) max(sensor_data(n,:))]) end %% figure(6) for n = 15 : 20 subplot(6, 1, n-14), plot(kgrid.t_array*scale, sensor_data.p(n,:), 'r-'); xlabel(['Time [' prefix 's]']); ylabel('Signal Amplitude-Pressure'); axis tight; % ylim([-max(sensor_data(n,:)) max(sensor_data(n,:))]) end RuiLiu on "Help on Attenuation !! Very Confused" http://www.k-wave.org/forum/topic/help-on-attenuation-very-confused#post-7849 Wed, 23 Sep 2020 21:09:02 +0000 RuiLiu 7849@http://www.k-wave.org/forum/ Thank you so much, It makes sense. But I still confuse about the attenuation. In order to make it clear, I make a 1D unit impulse. I set medium.coefficient at 0.2, power at 2. Then the pressure result seems suddenly drop to 0.5, and then not decrease anymore. I don't know how that happened. Also, even if I do not put any value for coefficient and power, It gives the same result. By the way, Is that possible to validate the impulse output wave using any general equation? I just trying to do some validation for impulse on both 1D and 2D medium. Thanks again. clc clear close all % ========================================================================= % SIMULATION % ========================================================================= % create the computational grid Nx = 500; % number of grid points in the x (row) direction dx = 0.05e-3; % grid point spacing in the x direction [m] kgrid = makeGrid(Nx, dx); % define the properties of the propagation medium medium.sound_speed = 1000; % [m/s] medium.density = 1000; % [kg/m^3] medium.alpha_coeff = 0.2; % [dB/(MHz^2 cm)] medium.alpha_power = 2; % create initial pressure distribution using a smoothly shaped sinusoid x_pos = 200; % [grid points] width = 100; % [grid points] height = 1; % [au] in = (0:pi/(width/2):2*pi).'; source.p0 = [((height/2)*sin(in-pi/2)+(height/2)); zeros(x_pos, 1); zeros(Nx - x_pos - width - 1, 1)]; plot(source.p0); % create a Cartesian sensor mask sensor.mask = [-5e-3, 0, 5e-3, 10e-3]; % [mm] % create time array t_end = 30e-6; [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed, [], t_end); % [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed); % run the simulation sensor_data = kspaceFirstOrder1D(kgrid, medium, source, sensor, 'PlotLayout', true); % ========================================================================= % VISUALISATION % ========================================================================= %% t1 = abs(-5e-3) /1000; % plot the recorded time signals figure; plot(kgrid.t_array*1e6,sensor_data(1, :), 'b-'); hold on; scatter(t1*1e6,0.5,'o') hold on; plot(kgrid.t_array*1e6,sensor_data(2, :), 'r-'); hold on; scatter(2*t1*1e6,0.5,'o') hold on; plot(kgrid.t_array*1e6,sensor_data(3, :), 'k-'); hold on; scatter(3*t1*1e6,0.5,'o') hold on; plot(kgrid.t_array*1e6,sensor_data(4, :), 'm-'); hold on; scatter(4*t1*1e6,0.5,'o') axis tight; set(gca, 'YLim', [-0.1 0.7]); ylabel('Pressure'); xlabel('Time, us'); legend('Sensor Position 1', 'Sensor Position 2', 'Sensor Position 3', 'Sensor Position 4'); Bradley Treeby on "Help on Attenuation !! Very Confused" http://www.k-wave.org/forum/topic/help-on-attenuation-very-confused#post-7842 Tue, 22 Sep 2020 13:47:06 +0000 Bradley Treeby 7842@http://www.k-wave.org/forum/ A pressure source in k-Wave corresponds to the injection of mass in free field. This is different to enforcing a boundary condition. This means that the pressure at the source position during the simulation will not match the pressure of the input signal. One way to think about it is that at each time step, you are adding the input signal to the pressure that is already at that position. At the next time point, this signal will not have had time to completely travel away, so the next sample from the input signal will be added to this. Hope that helps, Brad RuiLiu on "Help on Attenuation !! Very Confused" http://www.k-wave.org/forum/topic/help-on-attenuation-very-confused#post-7832 Wed, 16 Sep 2020 18:36:49 +0000 RuiLiu 7832@http://www.k-wave.org/forum/ Hi, I am trying to vaild my source shape. So, I put my source and sensor very close and medium-mode to no-absorption, expected to get the same signal in the sensor. But from the result even I put source and sensor in the same location, the amplitude decrease a lot. I am not sure which part I am setting wrongly. Any help will be appreciate. Thanks. clc clear close all % ========================================================================= % SIMULATION % ========================================================================= % create the computational grid Nx = 128; % number of grid points in the x (row) direction Ny = 128; % number of grid points in the y (column) direction dx = 50e-3/Nx; % grid point spacing in the x direction [m] dy = dx; % grid point spacing in the y direction [m] kgrid = makeGrid(Nx, dx, Ny, dy); % define the properties of the propagation medium medium.sound_speed = 1500; % [m/s] medium.alpha_coeff = 0; % [dB/(MHz^y cm)] medium.alpha_power = 2; medium.alpha_mode = 'no_absorption';% no-dispersion % create the time array [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed); % define a single source point source.p_mask = zeros(Nx, Ny); source.p_mask(end - Nx/4, Ny/2) = 1; % define a time varying sinusoidal source source_freq = 0.25e6; % [Hz] source_mag = 2; % [Pa] source.p = source_mag*sin(2*pi*source_freq*kgrid.t_array); % filter the source to remove high frequencies not supported by the grid % source.p = filterTimeSeries(kgrid, medium, source.p); % define a single sensor point sensor.mask = zeros(Nx, Ny); sensor.mask(end - Nx/4 - 1, Ny/2) = 1; sensor.mask(end - 3*Nx/4, Ny/2) = 1; % define the acoustic parameters to record sensor.record = {'p'}; % run the simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor); % ========================================================================= % VISUALISATION % ========================================================================= % plot the simulated sensor data figure; [t_sc, scale, prefix] = scaleSI(max(kgrid.t_array(:))); subplot(3, 1, 1), plot(kgrid.t_array*scale, source.p, 'k-'); xlabel(['Time [' prefix 's]']); ylabel('Signal Amplitude'); axis tight; title('Input Pressure Signal'); subplot(3, 1, 2), plot(kgrid.t_array*scale, sensor_data.p(2,:), 'r-'); xlabel(['Time [' prefix 's]']); ylabel('Signal Amplitude'); axis tight; title('Sensor Pressure Signal'); subplot(3, 1, 3), plot(kgrid.t_array*scale, sensor_data.p(1,:), 'r-'); xlabel(['Time [' prefix 's]']); ylabel('Signal Amplitude'); axis tight; title('Sensor Pressure Signal');