k-Wave

1. cge
   

   Member
   Joined: May '14
   Posts: 16

   Hi everyone

   Recently I am helping another guy in our group to simulate the reflection at the hetrogeneous media boundary.

   I started from a simple 1D simulation, the code is provided below.

   clear all;

   % =========================================================================
   % SIMULATION
   % =========================================================================

   % create the computational grid
   Nx = 64; % 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 = 1540*ones(Nx, 1); % [m/s]
   medium.sound_speed(11:50) = 2690;
   medium.sound_speed(51:64) = 340; % [m/s]
   medium.density = 1000*ones(Nx, 1); % [kg/m^3]
   medium.density(11:50) = 1190; % [kg/m^3]
   medium.density(51:64) = 1.204;
   medium.alpha_power=2;
   medium.alpha_coeff=0.002*ones(Nx,1);
   medium.alpha_coeff(11:50)=1.19;
   medium.alpha_coeff(51:64)=1.64;
   % create initial pressure distribution using a smoothly shaped sinusoid
   source_mask=zeros(Nx,1);
   source_mask(1)=1;
   source.p0 =10^6*source_mask;

   % create a Cartesian sensor mask
   sensor_mask=zeros(Nx,1)
   sensor_mask(9)=1;
   sensor.mask = sensor_mask; % [mm]

   % set the simulation time to capture the reflections
   t_end = dx*10/1540+dx*40/2690+dx*14/340;

   % define the time array
   [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed, [], t_end);

   % run the simulation
   sensor.record={'p'};
   sensor_data = kspaceFirstOrder1D(kgrid, medium, source, sensor, 'PlotLayout', true);

   % =========================================================================
   % VISUALISATION
   % =========================================================================

   % plot the recorded time signals
   figure;
   plot(kgrid.t_array, abs(sensor_data.p));
   xlabel('time step');
   ylabel('pressure at the sensor position');

   As for the result, a interesting happens, the excitation has a amplitude only of 10^6 Pa, however, the result I got in this end is quite huge, almost 10^60~70 pa, almost infinity and impossible.

   I checked the code but can not find the problem. Can anyone helo?

   Best Regards

   Posted 9 years ago #

2. Anthony
   

   Member
   Joined: Apr '13
   Posts: 96

   Hi cge,

   I am not part of the development team but I would say that k-wave is accurate only for weakly heterogeneous media.

   What if you try other values for the second medium sound speed, density and attenuation? Is there a threshold on the difference of impedance under which this strange peak disappears?

   Intuitively, it seems natural that a spectral method is not comfortable with discontinuities. I have not investigated in details the references mentioned at the bottom of the 21st page of the user's manual but maybe you can try to play with the CFL number. You could also smooth your sound speed and density distributions (see p63 of the user's manual).

   Maybe you could also set a source mask at the place where the interface should be and use Dirichlet boundary conditions (see p34), in order to have a reflecting interface.

   Best regards,
   Anthony

   Posted 9 years ago #

3. cge
   

   Member
   Joined: May '14
   Posts: 16

   Hi Hypermoi

   Thank you very much for your suggestion. I will take a try

   Best Regards

   Posted 9 years ago #

4. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,640

   Hi cge,

   There are a few problems with your simulation:

   (1) Your source is defined within the perfectly matched layer (PML). The PML is a thin layer that surrounds the domain and absorbs the outgoing waves to simulate free-field conditions. The default PML thickness for 2D simulations is 20 grid points. You can either move your source to be outside the PML, or set the PML to be outside the domain you specify by using the optional input 'PMLInside', false.

   (2) You are using a initial pressure source at the very edge of the grid. By default, the source.p0 is smoothed inside the simulation functions. This means that parts of the source will be wrapped to the other side of the domain. To see what I mean, try running the following code:

   % define grid
   kgrid = makeGrid(20, 1);
   
   % define source and smooth
   p0 = [1; zeros(19, 1)];
   p0_smooth = smooth(kgrid, p0, true);
   
   % plot
   figure;
   stem(p0, 'bx');
   hold on;
   stem(p0_smooth, 'r');

   (3) The medium interfaces are defined within the PML (see discussion above).

   (4) As Anthony has already mentioned, the third layer has a huge impedance mismatch changing from 3.2 MRayls to 0.4 kRayl (a factor of nearly 10,000). The problem we have found with such large changes is that the transmission coefficient when waves travel from the low to high impedance material is vastly over-estimated (except using high numbers of points per wavelength), resulting in very large amplitude waves. However, in your case, this shouldn't be a problem as the source originates from within the high-impedance material. The CFL number also affects the accuracy of the simulation as mentioned above.

   A revised script is below.

   Brad.

   clear all;
   
   % =========================================================================
   % SIMULATION
   % =========================================================================
   
   % create the computational grid, with the PML outside the grid
   PML_size = 20;          % [grid points]
   Nx = 128 - 2*PML_size;  % [grid points]
   dx = 0.05e-3;           % [m]
   kgrid = makeGrid(Nx, dx);
   
   % define position of source and interfaces
   source_pos      = 15;   % [grid points]
   sensor_pos      = 30;   % [grid points]
   interface_1_pos = 35;   % [grid points]
   interface_2_pos = 60;   % [grid points]
   
   % define material properties
   c1   = 1540;  % [m/s]
   rho1 = 1000;  % [kg/m^3]
   c2   = 2690;  % [m/s]
   rho2 = 1190;  % [kg/m^3]
   c3   = 340;   % [m/s]
   rho3 = 1.2;   % [kg/m^3]
   
   % define sound speed map
   medium.sound_speed = c1*ones(Nx, 1);
   medium.sound_speed(interface_1_pos:interface_2_pos - 1) = c2;
   medium.sound_speed(interface_2_pos:end) = c3;
   
   % define density
   medium.density = rho1*ones(Nx, 1);
   medium.density(interface_1_pos:interface_2_pos - 1) = rho2;
   medium.density(interface_2_pos:end) = rho3;
   
   % set the simulation time to capture the reflections
   t_end = Nx*dx/min(medium.sound_speed);
   
   % define the time array
   CFL = 0.05;
   [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed, CFL, t_end);
   
   % define source, offset from the edge to prevent smoothing from wrapping
   % the source
   source_mag = 1e6;
   source.p0 = zeros(Nx,1);
   source.p0(source_pos) = source_mag;
   
   % create a binary sensor mask
   sensor_mask = zeros(Nx,1);
   sensor_mask(sensor_pos) = 1;
   sensor.mask = sensor_mask;
   
   % run the simulation
   sensor.record = {'p'};
   sensor_data = kspaceFirstOrder1D(kgrid, medium, source, sensor,...
       'PMLInside', false, 'PlotLayout', true, ...
       'PlotScale', [-1, 1]*source_mag, 'PlotFreq', 10);
   
   % =========================================================================
   % VISUALISATION
   % =========================================================================
   
   % plot the recorded time signals
   figure;
   plot(kgrid.t_array, sensor_data.p);
   xlabel('time step');
   ylabel('pressure at the sensor position');

   Posted 9 years ago #

5. cge
   

   Member
   Joined: May '14
   Posts: 16

   Hi Brad

   Thank you very much. Now it works and I have got the correct result

   Posted 9 years ago #

6. peet2521
   

   Member
   Joined: Sep '16
   Posts: 1

   Brad,

   Are there any case studies or something on the grid resolution we should use to get a correct transmission/reflection for a given impedance mismatch ratio?

   Posted 7 years ago #

7. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,640

   Hi peet2521,

   We don't have the data presented in exactly that form, but there are some plots in this paper that might be useful.

   Brad.

   Posted 7 years ago #

RSS feed for this topic

Reply

You must log in to post.