Hi , I have to model a "room (2 by 2 m) " in 2-D with a single acoustic source and detect pressure fields at all grid points defined . For now i am using the heterogeneous medium example , inputting density and speed of sound 330m/s,1.225kg/m3 everywhere in the grid except bounding the four edges of grid with density and speed of sound in CONCRETE . I am not able to get any reverberations recorded at sensor points infact recorded signals are decayed . Is this the correct way of constructing the "room" ?
k-Wave
A MATLAB toolbox for the time-domain
simulation of acoustic wave fields
Modelling a 2D simple room(bounded by 4 walls ) with an audible sound source .
(12 posts) (3 voices)-
Posted 10 years ago #
-
Hi Gurashish,
Have you taken into account the perfectly matched layer (PML)? This is a thin layer surrounding the computational grid that absorbs the waves as they reach the edges. In 2D, the default size of the PML is 20 grid points positioned inside the computational grid defined by the user. You can set the PML to be outside the computational grid by setting
'PMLInside'tofalse, i.e.,sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, 'PMLInside', false);Alternatively, you could make your walls fat enough that they extend through the PML and into the interior of the domain, i.e., make them thicker than 20 grid points. Keep in mind your simulation will be fastest when the total grid size including the PML is a number with small prime factors. You can use the function
checkFactorsto identify suitable values.If you're happy modifying the simulation functions, another way to model the room is suggested in this post.
Hope that helps,
Brad.
Posted 10 years ago # -
Yes ! i have disabled the PML as well as made the walls realistically thicker , but the recorded signals are not matching the source signal at all . Though what i expect is that the sound signals would hit the concrete walls and reverberate in the room. Recorded signals would atleast have similarity with source .
This is exactly what i use .
% create the computational grid [ 3m by 3m domain ]
Nx = 20; % number of grid points in the x (row) direction
Ny = 20; % number of grid points in the y (column) direction
dx = 0.15; % grid point spacing in the x direction [m]
dy = 0.15; % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);% define the properties of the propagation medium [ bounded by concrete walls 0.9m thick]
medium.sound_speed = 330*ones(Nx, Ny); % [m/s]
medium.sound_speed(1:7,:) = 3200; % [m/s]
medium.sound_speed(:,1:7) = 3200;
medium.sound_speed(:,14:Nx) = 3200;
medium.sound_speed(14:Nx,:) = 3200;
medium.density = 1.225*ones(Nx, Ny); % [kg/m^3]
medium.density(1:7,:) = 2300; % [kg/m^3]
medium.density(:,1:7) = 2300;
medium.density(14:Nx,:) = 2300;
medium.density(:,14:Nx) = 2300;% create the time array
% [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);% % define a time varying sinusoidal source
% source_freq = 7000; % [Hz]
% source_mag = 50; % [Pa]
% source.p = source_mag*sin(2*pi*source_freq*kgrid.t_array);% define a single source point
source.p_mask = zeros(Nx, Ny);
source.p_mask(10,10) = 1;% [ define sensors everywhere ]
sensor.mask = ones(Nx, Ny);% define the acoustic parameters to record
sensor.record = {'p'};% layout in addition to removing the PML from the display
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor,'PMLInside', false);Posted 10 years ago # -
It seems your source frequency is too high to be supported by the spatial grid. The grid needs to have a least two points per wavelength at the highest frequency of interest, where f_max = c0 / (2 * dx). There's more information about this in the k-Wave manual if you're stuck.
Brad.
---
clear all; % set the PML size PML_size = 20; % create the computational grid [ 3m by 3m domain ] Nx = 128 - 2*PML_size; % number of grid points in the x (row) direction Ny = 128 - 2*PML_size; % number of grid points in the y (column) direction dx = 3/Nx; % grid point spacing in the x direction [m] dy = 3/Ny; % grid point spacing in the y direction [m] kgrid = makeGrid(Nx, dx, Ny, dy); % define the properties of the propagation medium c_air = 330; % [m/s] rho_air = 1.225; % [kg/m^3] % define the properties of the wall c_wall = 1000; % [m/s] rho_wall = 1000; % [kg/m^3] thickness = 5; % [grid points] % define the position of the wall wall = zeros(Nx, Ny); wall(1:thickness, :) = 1; wall(end - thickness + 1:end, :) = 1; wall(:, 1:thickness) = 1; wall(:, end - thickness + 1:end) = 1; % assign the medium properties medium.sound_speed = c_air*ones(Nx, Ny); medium.sound_speed(wall == 1) = c_wall; medium.density = rho_air*ones(Nx, Ny); medium.density(wall == 1) = rho_wall; % set the reference sound speed used in the k-space operator medium.sound_speed_ref = c_air; % create the time array cfl = 0.1; t_end = 20e-3; % [s] [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed, cfl, t_end); % define a time varying sinusoidal source source_freq = 1.5e3; % [Hz] source_mag = 50; % [Pa] source.p = source_mag*toneBurst(1/kgrid.dt, source_freq, 3); % define a single source point source.p_mask = zeros(Nx, Ny); source.p_mask(60, 40) = 1; % define sensors everywhere sensor.mask = ones(Nx, Ny); % define the acoustic parameters to record sensor.record = {'p'}; % run simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, ... 'PMLInside', false, 'PlotPML', false, 'PMLSize', PML_size, ... 'DisplayMask', 'off', 'PlotScale', [-0.25, 0.25]*source_mag,... 'DataCast', 'single'); % plot the recorded data figure; subplot(2, 1, 1) plot(kgrid.t_array(1:length(source.p)), source.p); set(gca, 'XLim', [0, t_end]); title('input'); subplot(2, 1, 2); plot(kgrid.t_array, sensor_data.p(end/2, :)); set(gca, 'XLim', [0, t_end]); title('output in centre');Posted 10 years ago # -
to validate this model . I m trying to see if the time delay between the signals recieved at two nodes corresponds to the speed of sound in air . ie distance between nodes / time delay . I am not getting the required speed . Is there a way to validate that this model is completely equivalent to a source present in a room , Since i dont want my experimental results to differ extremely from the numerical modelling ? I tried using the Xcorr function to find lag between two signals , but its confusing .
P.S : I am Researching on Application of Time Reversal in Localising a Source in Room . So if you have any suggestions regarding some experimentation/modelling it would be interesting to know !
Posted 10 years ago # -
Hi Gurashish,
There are many ways you could do this. For example, you could replace the time varying source with an impulsive source by defining
source.p0instead ofsource.p, and then compare the time of arrival between two sensor points spaced some distance apart. If you're having trouble with defining the sensor mask or interpreting the sensor data, I'd suggest taking a looking at the k-Wave Manual.Good luck with your simulations,
Brad.
Posted 10 years ago # -
Hi brad,
Thanks so much for the help firstly . Using your above described code i m afraid of the colour of the pressure field that propagates through the walls while the simulation shows the pressure field being generated(it becomes RED as sound waves start to reflect from them, while the tone burst source has YELLOW waves propagating that makes me feel that most of the waves are entering entering the walls . I m afraid if much of waves are passing through the walls into the PML described outside the grid and hence leading to attenuation of energy. Incase you feel the same . Is there a solution to this .?clear all;
% set the PML size
PML_size = 20;% create the computational grid [ 3m by 3m domain ]
Nx = 128 - 2*PML_size; % number of grid points in the x (row) direction
Ny = 128 - 2*PML_size; % number of grid points in the y (column) direction
dx = 3/Nx; % grid point spacing in the x direction [m]
dy = 3/Ny; % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);% define the properties of the propagation medium
c_air = 330; % [m/s]
rho_air = 1.225; % [kg/m^3]% define the properties of the wall
c_wall = 3200; % [m/s] concrete walls
rho_wall = 2300; % [kg/m^3]
thickness = 5; % [grid points]% define the position of the wall
wall = zeros(Nx, Ny);
wall(1:thickness, :) = 1;
wall(end - thickness + 1:end, :) = 1;
wall(:, 1:thickness) = 1;
wall(:, end - thickness + 1:end) = 1;% assign the medium properties
medium.sound_speed = c_air*ones(Nx, Ny);
medium.sound_speed(wall == 1) = c_wall;
medium.density = rho_air*ones(Nx, Ny);
medium.density(wall == 1) = rho_wall;% set the reference sound speed used in the k-space operator
medium.sound_speed_ref = c_air;% create the time array
cfl = 0.1;
t_end = 20e-3; % [s]
[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed, cfl, t_end);% define a time varying sinusoidal source
source_freq = 1.5e3; % [Hz]
source_mag = 50; % [Pa]
source.p = source_mag*toneBurst(1/kgrid.dt, source_freq, 3);% define a single source point
source.p_mask = zeros(Nx, Ny);
source.p_mask(60, 40) = 1;% define sensors everywhere
sensor.mask = ones(Nx, Ny);% define the acoustic parameters to record
sensor.record = {'p'};% run simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, ...
'PMLInside', false, 'PlotPML', false, 'PMLSize', PML_size, ...
'DisplayMask', 'off', 'PlotScale', [-0.25, 0.25]*source_mag,...
'DataCast', 'single');% plot the recorded data
figure;
subplot(2, 1, 1)
plot(kgrid.t_array(1:length(source.p)), source.p);
set(gca, 'XLim', [0, t_end]);
title('input');subplot(2, 1, 2);
plot(kgrid.t_array, sensor_data.p(end/2, :));
set(gca, 'XLim', [0, t_end]);
title('output in centre');Posted 10 years ago # -
Hi Gurashish,
Consider a plane wave that is normally incident on a boundary. Let the characteristic acoustic impedance of the first medium be Z1 = rho1*c1 and the impedance of the second medium be Z2 = rho2*c2. When Z2 >> Z1 (as in the case of air and concrete), the pressure reflection coefficient will approach 1, and the pressure transmission coefficient will approach 2. As k-Wave plots pressure, this is the reason the colour scale looks darker in the concrete.
However, in this case, the transmission of energy into the second medium will approach zero (where T_e = 1 - R_p^2). This is clear if you consider the particle velocity of the plane wave in the second medium, which is given by u = p / Z2. As the particle velocity approaches 0, the transmitted intensity (given by I = p*u) will also approach 0.
You can find more details in any introductory book on acoustics, for example, Kinsler & Frey.
Hope that helps,
Brad.
Posted 10 years ago # -
thank you ! that certainly clears my concepts . Incase if i just want to keep the boundaries completely reflective , without having to assign a number of grid points to the wall ! could i do that ? what change do i make to the code . ?
Posted 10 years ago # -
Hi Gurashish,
There are a few possibilities. You could try imposing a boundary condition explicitly as discussed here. However, this would still require the boundary to be several grid points thick. Alternatively, if you're familiar with spectral methods, you could replace the FFTs with discrete cosine and sine transforms to directly enforce Neumann and Dirichlet boundary conditions at the edge of the grid.
Brad.
Posted 10 years ago # -
Hello Bradely Treeby,
I want to have the reflection co-effiecnt= 0.5, or 0.6 for the walls arround. The medium.sound_speed =330 ; What should be the values of medium.alpha_coeff, medium.alpha_power and PML_size. What is the numerical relationship between reflection,PLM(Thickness in grid point) and PML absorption as mention at pag 17,18 and figure 2.3.Posted 9 years ago # -
Hi Zia,
There is no way to directly define an impedance boundary condition in k-Wave. However, you could define the material properties of your walls based on the equation for the plane wave pressure reflection coefficient at normal incidence, i.e., R = (Z2 - Z1) / (Z2 + Z1), where Z1 and Z2 are the characteristic acoustic impedance of the two media.
Assuming the inner medium is air, then Z1 = 343 * 1.225 = 420. For R = 0.5, then Z2 = Z1 * (1 + R) / (1 - R) = 1260. You can then split this into the desired sound speed and density in the second medium.
We have not derived an expression between reflection coefficient and the PML parameters - the data in Fig. 2.3 was derived by running numerical experiments. Generally speaking, the default values of PML thickness (20 grid points) and PML absorption (2 Nepers per grid point) work well. If you want to test a different range it's reasonably straightforward to do - send a wave into the PML and measure the amplitude of the transmitted and reflected waves.
Brad.
Posted 9 years ago #
Reply
You must log in to post.